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

Question

Alvin 2023-8-31

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2015246-how-to-do-lu-factorization-without-permutation-the-lu-is-with-permutation

编辑： David Goodmanson 2023-9-6

采纳的回答： Bruno Luong

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Bruno Luong 2023-8-31

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2015246-how-to-do-lu-factorization-without-permutation-the-lu-is-with-permutation#answer_1298051

编辑：Bruno Luong 2023-8-31

在 MATLAB Online 中打开

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

Bruno Luong 2023-9-5

编辑：Bruno Luong 2023-9-5

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 2023-9-6

Thanks!

请先登录，再进行评论。

Answer 2

David Goodmanson 2023-9-6

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2015246-how-to-do-lu-factorization-without-permutation-the-lu-is-with-permutation#answer_1302771

编辑：David Goodmanson 2023-9-6

在 MATLAB Online 中打开

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
  

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Bruno Luong 2023-9-6

编辑：Bruno Luong 2023-9-6

在 MATLAB Online 中打开

@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