How to do LU factorization without permutation？ The lu() is with permutation.

How to do LU factorization without permutation？ The lu() is with permutation.

Answer 1

[Bruno Luong]

It is not possible for the simple reason that such decomposition without allowed permutation might not possible. For example I claim there is no such "non-permuted lu" decomposition for

A=[0 1;
   1 0]

meaning there is no lower L and upper U such that A=L*U

[Alvin]

Sorry for the confusion, dude. What puzzled me is the case that even if it can do LU facorization without permutation, MatLab will still permute rows for larger pivots and smaller errors. For example, [1,1;2,3] will first be permuted to [2,3;1,1] and then be factorized as [1,0;0.5,1]*[2,3;0,-1]. Is there any way to factorize it as [1,0;2,1]*[1,1;0,1], just like what the elementry linear algebra textbooks tell?

[Bruno Luong]

No because MATLAB LU (or any serious LU decomposition library) has internal strategy to permute so that the pivot to be devided is the largest number in term of aboslute value. This ensures a good accuracy of the decomposition and robust to round-off error.

[Alvin]

Thanks!

Answer 2

[David Goodmanson]

Hi Alvin

Certainly you can do the decomposition without permutation, just not with Matlab lu. (Why you might want to is a different question). For example, the code below does no permutations. Not surprisingly it is less accurate than Matlab lu which is shown for comparison. In this case the no-permutation error is still quite small, though. That's true in lots of other situations, but not always.

As you can see, the algorithm divides by A(k,k), so if any of the diagonal elements (except for the last one) happen to be small, there are going to be issues. In some cases the error goes out of control. Hence the use of row and/or column permutation, to put a larger element onto the diagonal.

A = randi(10,5,5)

  A = [9    10     4    10     3
       4     2     7     7     3
       5     7     9     3     4
      10     6     1     5     2
       1     7     1     2     4]

Astore = A;     % this algorithm overwrites A, so keep a copy
n = size(A,1);
for k = 1:n-1;
    q = k+1:n;
    A(q,k) = A(q,k)/A(k,k);
    A(q,q) = A(q,q) - A(q,k)*A(k,q);
end
U = triu(A)
L = tril(A,-1) + eye(n,n)
%check
ck = L*U -Astore

% compared to
[Lmat Umat] = lu(Astore)
% check
ck2 = Lmat*Umat - Astore

% results

  U =
  9.0000   10.0000    4.0000   10.0000    3.0000
  0   -2.4444    5.2222    2.5556    1.6667
  0         0    9.8636   -1.0455    3.3182
  0         0         0  -12.9770    0.0138
  0         0         0         0    3.2717

  L =
  1.0000         0         0         0         0
  0.4444    1.0000         0         0         0        % unpermuted lower
  0.5556   -0.5909    1.0000         0         0
  1.1111    2.0909   -1.4562    1.0000         0
  0.1111   -2.4091    1.3318   -0.6502    1.0000

  ck =
  1.0e-14 *
  0         0         0         0         0
  0         0         0         0         0
  0         0         0         0   -0.0444
  0    0.0888         0   -0.1776         0
  0         0         0    0.0888   -0.0444

  % Matlab lu
  
  Lmat =
  0.9000    0.7187    0.3091    0.8345    1.0000
  0.4000   -0.0625    0.8386    1.0000         0
  0.5000    0.6250    1.0000         0         0
  1.0000         0         0         0         0
  0.1000    1.0000         0         0         0
  
  Umat =
  10.0000    6.0000    1.0000    5.0000    2.0000
  0    6.4000    0.9000    1.5000    3.8000
  0         0    7.9375   -0.4375    0.6250
  0         0         0    5.4606    1.9134
  0         0         0         0   -3.3212
  
  ck2 =
  1.0e-15 *
  0         0         0         0   -0.4441
  0         0         0         0    0.4441
  0         0         0         0         0
  0         0         0         0         0
  0         0         0         0         0
  

[Bruno Luong]

@David Goodmanson Thanks, I adapt your code to this example that show without permutation fails miserably to find correct decomposition. (You can run as many time as you wnt for slightly different inputs).

A = rand(2)*1e-20 + [0 1; 
                     1 1];

Astore = A;     % this algorithm overwrites A, so keep a copy
n = size(A,1);
for k = 1:n-1;
    q = k+1:n;
    A(q,k) = A(q,k)/A(k,k);
    A(q,q) = A(q,q) - A(q,k)*A(k,q);
end
U = triu(A);
L = tril(A,-1) + eye(n,n);
%check
ck = L*U -Astore % A(2,2) = 1 cannot be recover without permutation

ck = 2×2

0 0 0 -1

% compared to
[Lmat, Umat] = lu(Astore);
% check
ck2 = Lmat*Umat - Astore

ck2 = 2×2

0 0 0 0