Yet another lu(A) question and pivoting

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.

Answer 1

[Christine Tobler]

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.

[David Wilson]

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)

Answer 2

[Christine Tobler]

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.