Use chol to factorize a symmetric coefficient matrix, and then solve a linear system using the Cholesky factor.

Create a symmetric matrix with positive values on the diagonal.

A = [1 0 1; 0 2 0; 1 0 3]

A = 3×3

     1     0     1
     0     2     0
     1     0     3

Calculate the Cholesky factor of the matrix.

R = chol(A)

R = 3×3

    1.0000         0    1.0000
         0    1.4142         0
         0         0    1.4142

Create a vector for the right-hand side of the equation $\mathrm{Ax}=\mathit{b}$.

b = sum(A,2);

Since $\mathit{A}={\mathit{R}}^{\mathit{T}}\mathit{R}$ with the Cholesky decomposition, the linear equation becomes ${\mathit{R}}^{\mathit{T}}\mathit{R}\mathit{x}=\mathit{b}$. Solve for x using the backslash operator.

x = R\(R'\b)

x = 3×1

    1.0000
    1.0000
    1.0000

Calculate the upper and lower Cholesky factorizations of a matrix and verify the results.

Create a 6-by-6 symmetric positive definite test matrix using the gallery function.

A = gallery('lehmer',6);

Calculate the Cholesky factor using the upper triangle of A.

R = chol(A)

R = 6×6

    1.0000    0.5000    0.3333    0.2500    0.2000    0.1667
         0    0.8660    0.5774    0.4330    0.3464    0.2887
         0         0    0.7454    0.5590    0.4472    0.3727
         0         0         0    0.6614    0.5292    0.4410
         0         0         0         0    0.6000    0.5000
         0         0         0         0         0    0.5528

Verify that the upper triangular factor satisfies R'*R - A = 0, within roundoff error.

norm(R'*R - A)

ans = 2.5801e-16

Now, specify the 'lower' option to calculate the Cholesky factor using the lower triangle of A.

L = chol(A,'lower')

L = 6×6

    1.0000         0         0         0         0         0
    0.5000    0.8660         0         0         0         0
    0.3333    0.5774    0.7454         0         0         0
    0.2500    0.4330    0.5590    0.6614         0         0
    0.2000    0.3464    0.4472    0.5292    0.6000         0
    0.1667    0.2887    0.3727    0.4410    0.5000    0.5528

Verify that the lower triangular factor satisfies L*L' - A = 0, within roundoff error.

norm(L*L' - A)

ans = 2.5801e-16

Use chol with two outputs to suppress errors when the input matrix is not symmetric positive definite.

Create a 5-by-5 matrix of binomial coefficients. This matrix is symmetric positive definite, so subtract 1 from the last element to ensure it is no longer positive definite.

A = pascal(5);
A(end) = A(end) - 1

A = 5×5

     1     1     1     1     1
     1     2     3     4     5
     1     3     6    10    15
     1     4    10    20    35
     1     5    15    35    69

Calculate the Cholesky factor for A. Specify two outputs to avoid generating an error if A is not symmetric positive definite.

[R,flag] = chol(A)

R = 4×4

     1     1     1     1
     0     1     2     3
     0     0     1     3
     0     0     0     1

flag = 5

Since flag is nonzero, it gives the pivot index where the factorization fails. chol is able to calculate q = flag-1 = 4 rows and columns correctly before failing when it encounters the part of the matrix that changed.

Verify that R'*R returns four rows and columns that agree with A(1:q,1:q).

q = flag-1;
R'*R

ans = 4×4

     1     1     1     1
     1     2     3     4
     1     3     6    10
     1     4    10    20

A(1:q,1:q)

ans = 4×4

     1     1     1     1
     1     2     3     4
     1     3     6    10
     1     4    10    20

Calculate the Cholesky factor of a sparse matrix, and use the permutation output to create a Cholesky factor with fewer nonzeros.

Create a sparse positive definite matrix based on the west0479 matrix.

load west0479
A = west0479;
S = A'*A;

Calculate the Cholesky factor of the matrix two different ways. First specify two outputs, and then specify three outputs to enable row and column reordering.

[R,flag] = chol(S);
[RP,flagP,P] = chol(S);

For each calculation, check that flag = 0 to confirm the calculation is successful.

if ~flag && ~flagP
    disp('Factorizations successful.')
else
    disp('Factorizations failed.')
end

Factorizations successful.

Compare the number of nonzeros in chol(S) vs. the reordered matrix chol(P'*S*P). Best practice is to use the three output syntax of chol with sparse matrices, since reordering the rows and columns can greatly reduce the number of nonzeros in the Cholesky factor.

subplot(1,2,1)
spy(R)
title('Nonzeros in chol(S)')
subplot(1,2,2)
spy(RP)
title('Nonzeros in chol(P''*S*P)')

Use the 'vector' option of chol to return the permutation information as a vector rather than a matrix.

Create a sparse finite element matrix.

S = gallery('wathen',10,10);
spy(S)

Calculate the Cholesky factor for the matrix, and specify the 'vector' option to return a permutation vector p.

[R,flag,p] = chol(S,'vector');

Verify that flag = 0, indicating the calculation is successful.

if ~flag
    disp('Factorization successful.')
else
    disp('Factorization failed.')
end

Factorization successful.

Verify that S(p,p) = R'*R, within roundoff error.

norm(S(p,p) - R'*R,'fro')

ans = 2.1039e-13