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

49 次查看（过去 30 天）
Alvin
Alvin 2023-8-31
编辑： David Goodmanson 2023-9-6
How to do LU factorization without permutation？ The lu() is with permutation.

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

采纳的回答

Bruno Luong
Bruno Luong 2023-8-31
编辑：Bruno Luong 2023-8-31
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更早的评论
Bruno Luong
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.

请先登录，再进行评论。

更多回答（1 个）

David Goodmanson
David Goodmanson 2023-9-6
编辑：David Goodmanson 2023-9-6
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 个评论
Bruno Luong
Bruno Luong 2023-9-6
编辑：Bruno Luong 2023-9-6
Ran in:
@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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Decomposition 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by