## Documentation Center |

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 "psychologically 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 "psychologically
lower triangular matrix" (i.e a product of unit lower triangular and
permutation matrices) 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
"psychologically unit lower triangular" (i.e., a permuted unit lower
triangular matrix), 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.

Was this topic helpful?