Azzera filtri
Azzera filtri

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

26 visualizzazioni (ultimi 30 giorni)
How to do LU factorization without permutation? The lu() is with permutation.

Risposta accettata

Bruno Luong
Bruno Luong il 31 Ago 2023
Modificato: Bruno Luong il 31 Ago 2023
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 Commenti
Bruno Luong
Bruno Luong il 5 Set 2023
Modificato: Bruno Luong il 5 Set 2023
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.

Accedi per commentare.

Più risposte (1)

David Goodmanson
David Goodmanson il 6 Set 2023
Modificato: David Goodmanson il 6 Set 2023
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 Commento
Bruno Luong
Bruno Luong il 6 Set 2023
Modificato: Bruno Luong il 6 Set 2023
@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

Accedi per commentare.

Categorie

Scopri di più su Linear Algebra in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by