MATLAB Answers

Yet another lu(A) question and pivoting

20 views (last 30 days)
David Wilson
David Wilson on 4 Apr 2019
Answered: Christine Tobler on 8 Apr 2019
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 Comments

Sign in to comment.

Answers (2)

Christine Tobler
Christine Tobler on 4 Apr 2019
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.

  1 Comment

David Wilson
David Wilson on 8 Apr 2019
Yes, I understand that mldivide does check the matrix for L or PL, but exactly how? And is this fast algorithm and efficient? After all it needs to do this everytime! (according to the flow chart)
My crude attempt is below. I first load a large, sparse matrix & generate a permuted L. But suppose I want to check this.
Below I scan row by row, and pick off the position of the final non-zero element, then sort & then finally check if actually lower triangular. Is this the approach, or is there something more cunning that I have overlooked?
load west0479
A = west0479;
[pL,U] = lu(A)
% Now re-generate p
n = size(pL,1); % # of rows
p = NaN(n,1);
for i=1:n
r = pL(i,:);
%p(i) = max(find(r~=0,1,'last'));
p(i) = find(r~=0,1,'last');
end
[sp,p] = sort(p);
Lr = L(sp,:); spy(Lr)
istril(Lr)

Sign in to comment.


Christine Tobler
Christine Tobler on 8 Apr 2019
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 Comments

Sign in to comment.

Tags