ldl

Block LDL' factorization for Hermitian indefinite matrices

Description

Full Data

example

[L,D] = ldl(A) factorizes full matrix A into a permuted lower triangular matrix L and a block diagonal matrix D satisfying A = L*D*L'.

example

[L,D] = ldl(A,triangle), where triangle is "upper", uses the upper triangle of A to compute the factorization. By default, triangle is "lower", which uses the lower triangle of A to compute the factorization.

example

[L,D,P] = ldl(___) also returns the permutation matrix satisfying P'*A*P = L*D*L'. You can use any of the input argument combinations in previous syntaxes.

example

[L,D,P] = ldl(___,outputForm) returns the permutation information in the form specified by outputForm. Specify outputForm as "vector" to return the permutation information as a vector.

Sparse Data

[L,D,P] = ldl(S) returns the lower triangular factor, block diagonal factor, and permutation information satisfying P'*S*P = L*D*L' for real sparse matrix S.

example

[L,D,P,C] = ldl(S) also returns the scaling matrix satisfying P'*C*S*C*P = L*D*L'. This syntax is available only for sparse matrix inputs.

[L,D,P,C] = ldl(S,tol) specifies the pivot tolerance in the algorithm for factoring S.

[___] = ldl(___,triangle), where triangle is "upper", uses the upper triangle of real sparse S to compute the factorization. By default, triangle is "lower", which uses the lower triangle of S to compute the factorization. You can use any of the input and output argument combinations in previous syntaxes for sparse matrices.

[___] = ldl(___,outputForm) returns the permutation information in the form specified by outputForm. Specify outputForm as "vector" to return the permutation information as a vector.

Examples

collapse all

Calculate the full LDL factorization of a 3-by-3 matrix by using the ldl function with three output arguments.

A = [2 -1 0; -1 2 -1; 0 -1 1];
[L1,D1,P1] = ldl(A)
L1 = 3×3

1.0000         0         0
-0.5000    1.0000         0
0   -0.6667    1.0000

D1 = 3×3

2.0000         0         0
0    1.5000         0
0         0    0.3333

P1 = 3×3

1     0     0
0     1     0
0     0     1

Incorporate the permutation matrix P into the L factor by specifying two output arguments. Here, L2 is equal to P1*L1.

[L2,D2] = ldl(A)
L2 = 3×3

1.0000         0         0
-0.5000    1.0000         0
0   -0.6667    1.0000

D2 = 3×3

2.0000         0         0
0    1.5000         0
0         0    0.3333

P1*L1
ans = 3×3

1.0000         0         0
-0.5000    1.0000         0
0   -0.6667    1.0000

Solve a linear system by performing an LDL factorization and using the factors to simplify the problem. Compare solutions obtained by using the backslash operator, the ldl function, and a decomposition object.

Create a 3-by-3 matrix and 3-by-1 vector and solve the linear system A*x = b using the backslash operator.

A = [2 -3 4; -3 2 -3; 4 -3 1];
b = [1; 0; 0];
x1 = A\b
x1 = 3×1

-0.4118
-0.5294
0.0588

Compute the LDL decomposition of A. Use the factors to solve the system. Because P'*A*P = L*D*L', the linear system A*x = b can be rewritten as inv(P')*L*D*L'*inv(P)*x = b. Solving for x, the linear system can be rewritten as x = P*inv(L')*inv(D)*inv(L)*P'*b or x = P*(L'\(D\(L\(P'*b)))).

Precomputing the matrix factors prior to solving the linear system can improve performance when solving a linear system with several different values of b because the factorization occurs only once and does not need to be repeated.

[L,D,P] = ldl(A)
L = 3×3

1.0000         0         0
0    1.0000         0
-0.6429   -0.4286    1.0000

D = 3×3

2.0000    4.0000         0
4.0000    1.0000         0
0         0   -1.2143

P = 3×3

1     0     0
0     0     1
0     1     0

x2 = P*(L'\(D\(L\(P'*b))))
x2 = 3×1

-0.4118
-0.5294
0.0588

The decomposition object is useful for solving linear systems using specialized factorizations because you get the benefits of precomputing the matrix factors but do not need to know how to use the factors. Use the decomposition object with the "ldl" type to recreate the same results.

dA = decomposition(A,"ldl");
x3 = dA\b
x3 = 3×1

-0.4118
-0.5294
0.0588

The larger the matrix to factor, the more efficient it is to use a permutation vector. Using a permutation vector can also save on execution time in subsequent operations.

Create a 1000-by-1000 symmetric random matrix, and calculate the LDL factorization of the matrix. Specify three output arguments for the ldl function to return permutation information. Compare the size of the permutation information stored as a matrix and as a vector.

A = rand(1000);
A = A + A';
[L,D,P] = ldl(A);
[Lv,Dv,v] = ldl(A,"vector");
whos P v
Name         Size                Bytes  Class     Attributes

P         1000x1000            8000000  double
v            1x1000               8000  double

When using a permutation matrix, the factors satisfy the identity P'*A*P = L*D*L'. When using a permutation vector, the factors instead satisfy the identity A(v,v) = Lv*Dv*Lv'.

Calculate the LDL factorization of a 3-by-3 matrix using its lower triangle.

A = [2 -1 0; -3 2 -1; 0 -3 1];
[L,D] = ldl(A)
L = 3×3

1.0000         0         0
-1.5000    1.0000         0
0    1.2000    1.0000

D = 3×3

2.0000         0         0
0   -2.5000         0
0         0    4.6000

Calculate the LDL factorization of the same matrix using its upper triangle by specifying triangle as "upper".

[Lu,Du] = ldl(A,"upper")
Lu = 3×3

1.0000   -0.5000         0
0    1.0000   -0.6667
0         0    1.0000

Du = 3×3

2.0000         0         0
0    1.5000         0
0         0    0.3333

Calculate the LDL factorization of a random Hermitian matrix. Specify four output arguments for the ldl function to return a scaling matrix. Values in scaling matrix C1 have different scales, while the block diagonal matrix D1 is more balanced.

A = randn(7);
A = A + A';
d = logspace(0,-20,7);
A = A.*d.*d';
[L1,D1,P1,C1] = ldl(sparse(A));
full(C1)
ans = 7×7
1019 ×

0.0000         0         0         0         0         0         0
0    0.0000         0         0         0         0         0
0         0    0.0000         0         0         0         0
0         0         0    0.0000         0         0         0
0         0         0         0    0.0000         0         0
0         0         0         0         0    0.0046         0
0         0         0         0         0         0    6.7627

full(D1)
ans = 7×7

1.0000         0         0         0         0         0         0
0   -0.2722         0         0         0         0         0
0         0    1.2133         0         0         0         0
0         0         0   -3.6488         0         0         0
0         0         0         0    0.2950         0         0
0         0         0         0         0   -0.7251         0
0         0         0         0         0         0   13.1577

Calculate the LDL factorization using ldl with three output arguments. If you do not specify a scaling matrix output, matrices L2 and D2 have varied scales.

[L2,D2,P2] = ldl(sparse(A));
full(D2)
ans = 7×7

0.0000    0.0010         0         0         0         0         0
0.0010    1.0753         0         0         0         0         0
0         0    0.0000    0.0000         0         0         0
0         0    0.0000   -0.0000         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0

Input Arguments

collapse all

Input matrix, specified as a real or complex square matrix. ldl assumes that A is symmetric (Hermitian if complex) and uses only the lower triangle of A by default. You can check if A is Hermitian by using ishermitian(A).

Data Types: single | double
Complex Number Support: Yes

Sparse input matrix, specified as a real sparse matrix. ldl assumes that S is symmetric.

Data Types: double

Triangular factor of the input matrix, specified as "lower" or "upper". Use this option to specify that ldl should use the lower or upper triangle of the input matrix to compute the factorization. ldl assumes that the input matrix is symmetric (Hermitian if complex) and uses only the lower or upper triangle to perform computations.

If you specify three output arguments and triangle is "lower", then the decomposition satisfies P'*A*P = L*D*L' and L is lower triangular. If triangle is "upper", then the decomposition satisfies P'*A*P = L'*D*L and L is upper triangular.

Shape of permutation output, specified as "matrix" or "vector". If outputForm is "matrix" and the lower triangle is used, then the permutation matrix satisfies P'*A*P = L*D*L'. If outputForm is "vector" and the lower triangle is used, then the permutation vector satisfies A(P,P) = L*D*L'.

Pivot tolerance for real sparse matrices, specified as a scalar in the range [0,0.5]. By default, tol is 0.01. Using smaller values may result in faster factorization times and lead to sparser factors but can result in a less stable factorization.

Data Types: double

Output Arguments

collapse all

Lower triangular factor, returned as a square matrix. By default, L is a lower unit triangular matrix, which is a lower triangular matrix with ones on the diagonal. If triangle is "upper", then L is an upper unit triangular matrix.

Block diagonal factor, returned as a square matrix. D has 1-by-1 and 2-by-2 blocks on its diagonal.

Permutation information, returned as a matrix or, if outputForm is "vector", as a vector.

Scaling factor, returned as a square matrix. C is a positive sparse diagonal matrix that satisfies P'*C*S*C*P = L*D*L' for real sparse matrix S.

collapse all

Hermitian Matrix

• A square matrix, A, is Hermitian if it is equal to its complex conjugate transpose, A = A'.

In terms of the matrix elements,

${a}_{i,\text{\hspace{0.17em}}j}={\overline{a}}_{j,\text{\hspace{0.17em}}i}\text{\hspace{0.17em}}\text{\hspace{0.17em}}.$

• The entries on the diagonal of a Hermitian matrix are always real. Because real matrices are unaffected by complex conjugation, a real matrix that is symmetric is also Hermitian. For example, this matrix is both symmetric and Hermitian.

$A=\left[\begin{array}{cc}\begin{array}{c}1\\ 0\end{array}& \begin{array}{cc}\begin{array}{c}0\\ 2\end{array}& \begin{array}{c}1\\ 0\end{array}\end{array}\\ 1& \begin{array}{cc}0& 1\end{array}\end{array}\right]$

• The eigenvalues of a Hermitian matrix are real.

References

[1] Ashcraft, Cleve, Roger G. Grimes, and John G. Lewis. “Accurate Symmetric Indefinite Linear Equation Solvers.” SIAM Journal on Matrix Analysis and Applications 20, no. 2 (January 1998): 513–61. https://doi.org/10.1137/S0895479896296921.

[2] Anderson, E., Z. Bai, C. Bischof, L. S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenny, and D. Sorensen. LAPACK Users’ Guide. Philadelphia: Society for Industrial and Applied Mathematics, 1999. https://doi.org/10.1137/1.9780898719604.

[3] Duff, Iain S. “MA57---a Code for the Solution of Sparse Symmetric Definite and Indefinite Systems.” ACM Transactions on Mathematical Software 30, no. 2 (June 2004): 118–44. https://doi.org/10.1145/992200.992202.

Version History

Introduced before R2006a