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.