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:

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

Answer 1

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.)

Answer 2

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)

Answer 3

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.