Block LDL' factorization for Hermitian indefinite matrices
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')
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.
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.
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))));
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))));
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)))');
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))));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.
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;[1] Ashcraft, C., R.G. Grimes, and J.G. Lewis. “Accurate Symmetric Indefinite Linear Equations Solvers.” SIAM J. Matrix Anal. Appl. Vol. 20. Number 2, 1998, pp. 513–561.
[2] Duff, I. S. "MA57 — A new code for the solution of sparse symmetric definite and indefinite systems." Technical Report RAL-TR-2002-024, Rutherford Appleton Laboratory, 2002.