Yet another lu(A) question and pivoting
5 ビュー (過去 30 日間)
古いコメントを表示
Does Matlab automatically detect "psychologically" lower triangular matrices when solving Ax=b?
I.e. if you do an lu decomposition to get a (permuted) L and an upper triangular U, and solve the system explicitly with something like:
% solve A*x=b
[L,U] = lu(A)
x = U\(L\b) % solve two triangular systems, where one is not truely triangular.
Does Matlab "know" that when solving L\b, L is a (permuted) lower triangular, and take advantage of that?
If so, how does it detect that? Note that in this case the "P" orthognal permutation is not exported.
BTW: I'm glad the help file has removed the "psychologically" term, I never thought that helped understanding in the slightest.
0 件のコメント
回答 (2 件)
Christine Tobler
2019 年 4 月 4 日
Yes, MATLAB checks if L is a permuted triangular matrix. See the doc for mldivide - Algorithm for full inputs. However, it's still cheaper if you get the third output P from LU and use it directly - this way, backslash does not have to reconstruct the permutation vector and triangular matrix from L.
Christine Tobler
2019 年 4 月 8 日
I'm afraid I can't post the algorithm used in mldivide. Note that it's a bit more general: it also works if both the rows and the columns of a matrix have been permuted.
Here's a slight variation on your algorithm:
% Set up lower triangular A with dense diagonal:
A = tril(sprandn(100, 100, 0.1)) + speye(100);
A = A(randperm(100), :);
% Find permutation vector:
[i, j] = find(A); % find row and column indices of all non-zeros in A
p = accumarray(i, j, [], @max); % find maximal column index for each row
perm(p) = 1:length(p); % compute inverse permutation of p
% Check the result
istril(A(perm, :))
Note this will only work if all diagonal elements of the triangular matrix A are non-zero. The same restriction also applies to the algorithm used in mldivide: It does not find a triangular form for matrices that are structurally singular.
0 件のコメント
参考
カテゴリ
Help Center および File Exchange で Operating on Diagonal Matrices についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!