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.