Note: This page has been translated by MathWorks. Click here to see

To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

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. `ldl`

uses MA57 for sparse real symmetric `A`

.

`[L,D,P,S] = LDL(A,THRESH)`

uses `THRESH`

as
the pivot tolerance in MA57. `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;

`ldl`

uses the MA57 routines in the Harwell
Subroutine Library (HSL) for real sparse matrices.

[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.

Was this topic helpful?