# 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, 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.