Main Content

ldl

Block LDL' factorization for Hermitian indefinite matrices

Syntax

L = ldl(A)
[L,D] = ldl(A)
[L,D,P] = ldl(A)
[L,D,p] = ldl(A,'vector')
[U,D,P] = ldl(A,'upper')
[U,D,p] = ldl(A,'upper','vector')
[L,D,P,S] = ldl(A)
[L,D,P,S] = LDL(A,THRESH)
[U,D,p,S] = LDL(A,THRESH,'upper','vector')

Description

L = ldl(A) returns only the permuted lower triangular matrix L as in the two-output form. The permutation information is lost, as is the block diagonal factor D. By default, ldl references only the diagonal and lower triangle of A, and assumes that the upper triangle is the complex conjugate transpose of the lower triangle. Therefore [L,D,P] = ldl(TRIL(A)) and [L,D,P] = ldl(A)both return the exact same factors. Note, this syntax is not valid for sparse A.

[L,D] = ldl(A) stores a block diagonal matrix D and a permuted lower triangular matrix in L such that A = L*D*L'. The block diagonal matrix D has 1-by-1 and 2-by-2 blocks on its diagonal. Note, this syntax is not valid for sparse A.

[L,D,P] = ldl(A) returns unit lower triangular matrix L, block diagonal D, and permutation matrix P such that P'*A*P = L*D*L'. This is equivalent to [L,D,P] = ldl(A,'matrix').

[L,D,p] = ldl(A,'vector') returns the permutation information as a vector, p, instead of a matrix. The p output is a row vector such that A(p,p) = L*D*L'.

[U,D,P] = ldl(A,'upper') references only the diagonal and upper triangle of A and assumes that the lower triangle is the complex conjugate transpose of the upper triangle. This syntax returns a unit upper triangular matrix U such that P'*A*P = U'*D*U (assuming that A is Hermitian, and not just upper triangular). Similarly, [L,D,P] = ldl(A,'lower') gives the default behavior.

[U,D,p] = ldl(A,'upper','vector') returns the permutation information as a vector, p, as does [L,D,p] = ldl(A,'lower','vector'). A must be a full matrix.

[L,D,P,S] = ldl(A) returns unit lower triangular matrix L, block diagonal D, permutation matrix P, and scaling matrix S such that P'*S*A*S*P = L*D*L'. This syntax is only available for real sparse matrices, and only the lower triangle of A is referenced.

[L,D,P,S] = LDL(A,THRESH) uses THRESH as the pivot tolerance in the algorithm. THRESH must be a double scalar lying in the interval [0, 0.5]. The default value for THRESH is 0.01. Using smaller values of THRESH may give faster factorization times and fewer entries, but may also result in a less stable factorization. This syntax is available only for real sparse matrices.

[U,D,p,S] = LDL(A,THRESH,'upper','vector') sets the pivot tolerance and returns upper triangular U and permutation vector p as described above.

Examples

These examples illustrate the use of the various forms of the ldl function, including the one-, two-, and three-output form, and the use of the vector and upper options. The topics covered are:

Before running any of these examples, you will need to generate the following positive definite and indefinite Hermitian matrices:

A = full(delsq(numgrid('L', 10)));
B = gallery('uniformdata',10,0);
M = [eye(10) B; B' zeros(10)]; 

The structure of M here is very common in optimization and fluid-flow problems, and M is in fact indefinite. Note that the positive definite matrix A must be full, as ldl does not accept sparse arguments.

Example 1 — Two-Output Form of ldl

The two-output form of ldl returns L and D such that A-(L*D*L') is small, L is permuted unit lower triangular, and D is a block 2-by-2 diagonal. Note also that, because A is positive definite, the diagonal of D is all positive:

[LA,DA] = ldl(A); 
fprintf(1, ...
'The factorization error ||A - LA*DA*LA''|| is %g\n', ...
norm(A - LA*DA*LA'));
neginds = find(diag(DA) < 0)

Given a b, solve Ax=b using LA, DA:

bA = sum(A,2);
x = LA'\(DA\(LA\bA));
fprintf(...
'The absolute error norm ||x - ones(size(bA))|| is %g\n', ...
norm(x - ones(size(bA)))); 

Example 2 — Three Output Form of ldl

The three output form returns the permutation matrix as well, so that L is in fact unit lower triangular:

[Lm, Dm, Pm] = ldl(M);
fprintf(1, ...
'The error norm ||Pm''*M*Pm - Lm*Dm*Lm''|| is %g\n', ...
norm(Pm'*M*Pm - Lm*Dm*Lm'));
fprintf(1, ...
'The difference between Lm and tril(Lm) is %g\n', ...
norm(Lm - tril(Lm)));

Given b, solve Mx=b using Lm, Dm, and Pm:

bM = sum(M,2);
x = Pm*(Lm'\(Dm\(Lm\(Pm'*bM))));
fprintf(...
'The absolute error norm ||x - ones(size(b))|| is %g\n', ...
norm(x - ones(size(bM)))); 

Example 3 — The Structure of D

D is a block diagonal matrix with 1-by-1 blocks and 2-by-2 blocks. That makes it a special case of a tridiagonal matrix. When the input matrix is positive definite, D is almost always diagonal (depending on how definite the matrix is). When the matrix is indefinite however, D may be diagonal or it may express the block structure. For example, with A as above, DA is diagonal. But if you shift A just a bit, you end up with an indefinite matrix, and then you can compute a D that has the block structure.

figure; spy(DA); title('Structure of D from ldl(A)');
[Las, Das] = ldl(A - 4*eye(size(A)));
figure; spy(Das); 
title('Structure of D from ldl(A - 4*eye(size(A)))');

Example 4 — Using the 'vector' Option

Like the lu function, ldl accepts an argument that determines whether the function returns a permutation vector or permutation matrix. ldl returns the latter by default. When you select 'vector', the function executes faster and uses less memory. For this reason, specifying the 'vector' option is recommended. Another thing to note is that indexing is typically faster than multiplying for this kind of operation:

[Lm, Dm, pm] = ldl(M, 'vector');
fprintf(1, 'The error norm ||M(pm,pm) - Lm*Dm*Lm''|| is %g\n', ...
  norm(M(pm,pm) - Lm*Dm*Lm'));

% Solve a system with this kind of factorization.
clear x;
x(pm,:) = Lm'\(Dm\(Lm\(bM(pm,:))));
fprintf('The absolute error norm ||x - ones(size(b))|| is %g\n', ...
  norm(x - ones(size(bM))));

Example 5 — Using the 'upper' Option

Like the chol function, ldl accepts an argument that determines which triangle of the input matrix is referenced, and also whether ldl returns a lower (L) or upper (L') triangular factor. For dense matrices, there are no real savings with using the upper triangular version instead of the lower triangular version:

Ml = tril(M);
[Lml, Dml, Pml] = ldl(Ml, 'lower'); % 'lower' is default behavior.
fprintf(1, ...
'The difference between Lml and Lm is %g\n', norm(Lml - Lm));
[Umu, Dmu, pmu] = ldl(triu(M), 'upper', 'vector');
fprintf(1, ...
'The difference between Umu and Lm'' is %g\n', norm(Umu - Lm'));

% Solve a system using this factorization.
clear x;
x(pm,:) = Umu\(Dmu\(Umu'\(bM(pmu,:))));
fprintf(...
'The absolute error norm ||x - ones(size(b))|| is %g\n', ...
norm(x - ones(size(bM))));

When specifying both the 'upper' and 'vector' options, 'upper' must precede 'vector' in the argument list.

Example 6 — linsolve and the Hermitian indefinite solver

When using the linsolve function, you may experience better performance by exploiting the knowledge that a system has a symmetric matrix. The matrices used in the examples above are a bit small to see this so, for this example, generate a larger matrix. The matrix here is symmetric positive definite, and below we will see that with each bit of knowledge about the matrix, there is a corresponding speedup. That is, the symmetric solver is faster than the general solver while the symmetric positive definite solver is faster than the symmetric solver:

Abig = full(delsq(numgrid('L', 30)));
bbig = sum(Abig, 2);
LSopts.POSDEF = false;
LSopts.SYM = false;
tic; linsolve(Abig, bbig, LSopts); toc;
LSopts.SYM = true;
tic; linsolve(Abig, bbig, LSopts); toc;
LSopts.POSDEF = true;
tic; linsolve(Abig, bbig, LSopts); toc;

References

[1] Ashcraft, Cleve, Roger G. Grimes, and John G. Lewis. “Accurate Symmetric Indefinite Linear Equation Solvers.” SIAM Journal on Matrix Analysis and Applications 20, no. 2 (January 1998): 513–561. https://doi.org/10.1137/S0895479896296921.

[2] Anderson, E., ed. LAPACK Users’ Guide. 3rd ed. Software, Environments, Tools. Philadelphia: Society for Industrial and Applied Mathematics, 1999. https://doi.org/10.1137/1.9780898719604.

[3] Duff, Iain S. “MA57---a Code for the Solution of Sparse Symmetric Definite and Indefinite Systems.” ACM Transactions on Mathematical Software 30, no. 2 (June 2004): 118–144. https://doi.org/10.1145/992200.992202.

Extended Capabilities

See Also

| |