How to do LU factorization without permutation？ The lu() is with permutation.
显示 更早的评论
How to do LU factorization without permutation？ The lu() is with permutation.
采纳的回答
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 个评论
Alvin
2023-9-5
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
更多回答（1 个）
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 个评论
@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
% compared to
[Lmat, Umat] = lu(Astore);
% check
ck2 = Lmat*Umat - Astore
类别
在 帮助中心 和 File Exchange 中查找有关 Linear Algebra 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
