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

回答1:

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],:);

回答2:

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

回答3:

My answer is here:

function [L, U, P] = LU_pivot(A)
[n, n1] = size(A); L=eye(n); P=eye(n); U=A;
for j = 1:n
  [pivot m] = max(abs(U(j:n, j)));     
  m = m+j-1;
  if m ~= j
    U([m,j],:) =  U([j,m], :);   % interchange rows m and j in U
    P([m,j],:) =  P([j,m], :);   % interchange rows m and j in P
    if j >= 2;    % very_important_point
      L([m,j],1:j-1) =  L([j,m], 1:j-1);   % interchange rows m and j in columns 1:j-1 of L
    end;
  end
  for i = j+1:n      
    L(i, j) = U(i, j) / U(j, j);
    U(i, :) =  U(i, :) - L(i, j)*U(j, :);
  end
end