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