LU decomposing a square matrix matlab gauss elimination

Question

I'm trying to create a program that takes a square (n-by-n) matrix as input, and if it is invertible, will LU decompose the matrix using Gaussian Elimination.

Here is my problem: in class we learned that it is better to change rows so that your pivot is always the largest number (in absolute value) in its column. For example, if the matrix was A = [1,2;3,4] then switching rows it is [3,4;1,2] and then we can proceed with the Gaussian elimination.

My code works properly for matrices that don't require row changes, but for ones that do, it does not. This is my code:

function newgauss(A)
    [rows,columns]=size(A);
    P=eye(rows,columns); %P is permutation matrix
    if(det(A)==0) %% determinante is 0 means no single solution
        disp('No solutions or infinite number of solutions')
        return;
    end
    U=A;
    L=eye(rows,columns);
    pivot=1;
    while(pivot<rows)
        max=abs(U(pivot,pivot));
        maxi=0;%%find maximum abs value in column pivot
        for i=pivot+1:rows
            if(abs(U(i,pivot))>max)
                max=abs(U(i,pivot));
                maxi=i;
            end
        end %%if needed then switch
        if(maxi~=0)
            temp=U(pivot,:);
            U(pivot,:)=U(maxi,:);
            U(maxi,:)=temp;
            temp=P(pivot,:);
            P(pivot,:)=P(maxi,:);
            P(maxi,:)=temp;
        end %%Grade the column pivot using gauss elimination
        for i=pivot+1:rows
            num=U(i,pivot)/U(pivot,pivot);
            U(i,:)=U(i,:)-num*U(pivot,:);
            L(i,pivot)=num;
        end
        pivot=pivot+1;
    end
    disp('PA is:');
    disp(P*A);
    disp('LU is:');
    disp(L*U);
end

Clarification: since we are switching rows, we are looking to decompose P (permutation matrix) times A, and not the original A that we had as input.

Explanation of the code:

First I check if the matrix is invertible, if it isn't, stop. If it is, pivot is (1,1)
I find the largest number in column 1, and switch rows
Grade column 1 using Gaussian elimination, turning all but the spot (1,1) to zero
Pivot is now (2,2), find largest number in column 2... Rinse, repeat

similarly instead of comparing `maxi~=0`, use a threshold. Also avoid using variable names that are also names of builtin functions (`max`), it is just confusing — Amro, Jul 22 '13 at 01:01

horchler · Accepted Answer · 2013-07-21T21:03:33.107

3

Your code seems to work fine from what I can tell, at least for the basic examples A=[1,2;3,4] or A=[3,4;1,2]. Change your function definition to:

function [L,U,P] = newgauss(A)

so you can output your calculated values (much better than using disp, but this shows the correct results too). Then you'll see that P*A = L*U. Maybe you were expecting L*U to equal A directly? You can also confirm that you are correct via Matlab's lu function:

[L,U,P] = lu(A);
L*U
P*A

Permutation matrices are orthogonal matrices, so P⁻¹ = P^T. If you want to get back A in your code, you can do:

P'*L*U

Similarly, using Matlab's lu with the permutation matrix output, you can do:

[L,U,P] = lu(A);
P'*L*U

(You should also use error or warning rather than how you're using disp in checking the determinant, but they probably don't teach that.)

edited Jul 21 '13 at 21:03

answered Jul 21 '13 at 20:17

horchler

18,384
4
37
73

I will try changing the definition but I really think there is something wrong with the code. for the matrix A=[1,2,3;4,5,6;7,8,9] PA is [7,8,9;1,2,3;4,5,6] and LU is [7,8,9;4,5.4286,6.8571;1,1.5714,2.1429]...not close – Oria Gruber Jul 21 '13 at 20:28
Maybe that's the example matrix you should have used in your question, rather than one that works in all cases. You're not performing Gaussian elimination properly -`U` isn't upper triangular in that case. – horchler Jul 21 '13 at 20:45
I'm sorry. I guess I should have used that example. can you pinpoint exactly what's wrong with the way I do gauss elimination? – Oria Gruber Jul 21 '13 at 20:49
Never mind, that example you gave is effectively a singular matrix. Do `det(A)` and look at the value -it's close to `eps`. You need to improve how you check for singular matrices (pretty easy) or not test your code with them. Try another 3-by-3 example to see if it works, but I think you still have issues. – horchler Jul 21 '13 at 20:51
1

You sir are absolutly right. Thank you very much. It works for other matrices 3x3, even ones that need to switch rows several times. just not ones where the determinante is close to zero – Oria Gruber Jul 21 '13 at 20:57
I still think that you have an issue with `L`. – horchler Jul 21 '13 at 21:02

score 1 · Answer 2 · answered Jul 22 '13 at 00:54

Note that the det function is implemented using an LU decomposition itself to compute the determinant... recursive anyone :)

Aside from that, there is a reminder towards the end of the page which suggest using cond instead of det to test for matrix singularity:

Testing singularity using abs(det(X)) <= tolerance is not recommended as it is difficult to choose the correct tolerance. The function cond(X) can check for singular and nearly singular matrices.

COND uses the singular value decomposition (see its implementation: edit cond.m)

score 0 · Answer 3 · answered Feb 11 '16 at 14:12

For anyone finding this in the future and needing a working solution:

The OP's code doesn't contain the logic for switching elements in L when creating the permutation matrix P. The adjusted code that gives the same output as Matlab's lu(A) function is:

function [L,U,P] = newgauss(A)
    [rows,columns]=size(A);
    P=eye(rows,columns); %P is permutation matrix
    tol = 1E-16; % I believe this is what matlab uses as a warning level
    if( rcond(A) <= tol) %% bad condition number
        error('Matrix is nearly singular')
    end
    U=A;
    L=eye(rows,columns);
    pivot=1;
    while(pivot<rows)
        max=abs(U(pivot,pivot));
        maxi=0;%%find maximum abs value in column pivot
        for i=pivot+1:rows
            if(abs(U(i,pivot))>max)
                max=abs(U(i,pivot));
                maxi=i;
            end
        end %%if needed then switch
        if(maxi~=0)
            temp=U(pivot,:);
            U(pivot,:)=U(maxi,:);
            U(maxi,:)=temp;
            temp=P(pivot,:);
            P(pivot,:)=P(maxi,:);
            P(maxi,:)=temp;

            % change elements in L-----
            if pivot >= 2
                temp=L(pivot,1:pivot-1);
                L(pivot,1:pivot-1)=L(maxi,1:pivot-1);
                L(maxi,1:pivot-1)=temp;
            end
        end %%Grade the column pivot using gauss elimination
        for i=pivot+1:rows
            num=U(i,pivot)/U(pivot,pivot);
            U(i,:)=U(i,:)-num*U(pivot,:);
            L(i,pivot)=num;
        end
        pivot=pivot+1;
    end
end

Hope this helps someone stumbling upon this in the future.

LU decomposing a square matrix matlab gauss elimination

3 Answers3