Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

LU decomposition with partial pivoting Matlab

I am trying to implement my own LU decomposition with partial pivoting. My code is below and apparently is working fine, but for some matrices it gives different results when comparing with the built-in [L, U, P] = lu(A) function in matlab

Can anyone spot where is it wrong?

function [L, U, P] = lu_decomposition_pivot(A)
    n = size(A,1);
    Ak = A;
    L = zeros(n);
    U = zeros(n);
    P = eye(n);
    for k = 1:n-1
        for i = k+1:n
            [~,r] = max(abs(Ak(:,k)));

            Ak([k r],:) = Ak([r k],:);
            P([k r],:) = P([r k],:);

            L(i,k) = Ak(i,k) / Ak(k,k);
            for j = k+1:n
                U(k,j-1) = Ak(k,j-1);
                Ak(i,j) = Ak(i,j) - L(i,k)*Ak(k,j);
            end
        end
    end
    L(1:n+1:end) = 1;
    U(:,end) = Ak(:,end);
return

Here are the two matrices I've tested with. The first one is correct, whereas the second has some elements inverted.

A = [1 2 0; 2 4 8; 3 -1 2];

A = [0.8443 0.1707 0.3111;
     0.1948 0.2277 0.9234;
     0.2259 0.4357 0.4302];

UPDATE

I have checked my code and corrected some bugs, but still there's something missing with the partial pivoting. In the first column the last two rows are always inverted (compared with the result of lu() in matlab)

function [L, U, P] = lu_decomposition_pivot(A)
    n = size(A,1);
    Ak = A;
    L = eye(n);
    U = zeros(n);
    P = eye(n);
    for k = 1:n-1
        [~,r] = max(abs(Ak(k:end,k)));
        r = n-(n-k+1)+r;
        Ak([k r],:) = Ak([r k],:);
        P([k r],:) = P([r k],:);
        for i = k+1:n
            L(i,k) = Ak(i,k) / Ak(k,k);
            for j = 1:n
                U(k,j) = Ak(k,j);
                Ak(i,j) = Ak(i,j) - L(i,k)*Ak(k,j);
            end
        end
    end
    U(:,end) = Ak(:,end);
return
like image 725
BRabbit27 Avatar asked Mar 08 '13 21:03

BRabbit27


People also ask

What is LU decomposition with partial pivoting?

LU Decomposition with Partial PivotingPA=LU. L is an n×n lower-triangular matrix with all diagonal entries equal to 1. U is an n×n upper-triangular matrix. P is an n×n permutation matrix.

How do you write LU decomposition in Matlab?

[ L , U ] = lu( A ) factorizes the full or sparse matrix A into an upper triangular matrix U and a permuted lower triangular matrix L such that A = L*U . [ L , U , P ] = lu( A ) also returns a permutation matrix P such that A = P'*L*U . With this syntax, L is unit lower triangular and U is upper triangular.

What is pivoting in LU factorization?

Pivoting for LU factorization is the process of systematically selecting pivots for Gaussian elimina- tion during the LU factorization of a matrix. The LU factorization is closely related to Gaussian elimination, which is unstable in its pure form.


2 Answers

I forgot that If there was a swap in matrix P I had to swap also the matrix L. So just add the next line after after swapping P and everything will work excellent.

L([k r],:) = L([r k],:);
like image 64
BRabbit27 Avatar answered Oct 13 '22 20:10

BRabbit27


Both functions are not correct. Here is the correct one.

function [L, U, P] = LU_pivot(A)
    [m, n] = size(A); L=eye(n); P=eye(n); U=A;
    for k=1:m-1
        pivot=max(abs(U(k:m,k)))
        for j=k:m
            if(abs(U(j,k))==pivot)
                ind=j
                break;
            end
        end
        U([k,ind],k:m)=U([ind,k],k:m)
        L([k,ind],1:k-1)=L([ind,k],1:k-1)
        P([k,ind],:)=P([ind,k],:)
        for j=k+1:m
            L(j,k)=U(j,k)/U(k,k)
            U(j,k:m)=U(j,k:m)-L(j,k)*U(k,k:m)
        end
        pause;
    end
end
like image 35
panch Avatar answered Oct 13 '22 20:10

panch