6.7. Cholesky Decomposition for Positive Definite Systems

Systems of linear equations of the form \mathbf{A}\,\bm{x} = \bm{x} are often symmetric and positive definite. Equations deriving from physical applications such as electrical, mechanical, and structural systems are often positive definite. We will see in Over-determined Systems and Vector Projections that the solutions to over-determined systems, such as found with least squares regression problems (Least Squares Regression), make use positive definite matrices of the form \mathbf{A}^T\mathbf{A}. The diagonal values of \mathbf{A}^T\mathbf{A} are positive because each is the dot product of a column of \bf{A} with itself. Likewise, the diagonal values of a positive definite matrix are always positive. The special symmetry of positive definite matrices allows for a simple and numerically accurate factoring called the Cholesky decomposition that returns an upper triangular matrix, \bf{R}. The transpose of \bf{R} gives us a lower triangular matrix. The product of the lower triangular matrix and upper triangular is the decomposition of \bf{A} that we will use as we did with LU decomposition to find the solution using a triangular solver.

(6.17)\mathbf{A} = \mathbf{R}^T\mathbf{R}.

The Cholesky decomposition can only be used for positive definite matrices. A positive definite matrix \bf{A} is real, symmetric, and satisfies the property

(6.18)\bm{x}^T\mathbf{A}\,\bm{x} > 0, \; \text{for all nonzero } \bm{x} \in \mathbb{R}^n.

The left side of equation (6.18) reduces to a dot product that yields a real, scalar number. Unfortunately, equation (6.18) is a difficult test to apply because it allows for an infinite set of \bm{x} vectors. An alternate criteria for evaluating if a matrix is positive definite consists of three tests.

  1. Is the matrix symmetric?
  2. Are all of the elements along the diagonal of the matrix positive?
  3. Are all of the eigenvalues of the matrix positive?

The procedure for finding the eigenvalues of a matrix is addressed in Application of Eigenvalues and Eigenvectors, but we will not use the third test. A faster final test to see if the Cholesky decomposition can be used is to attempt to use it. The decomposition fails if the algorithm needs to take the square root of a number that is zero or negative. So we can attempt to use the Cholesky decomposition algorithm if the matrix is symmetric and has only positive values along the diagonal. The decomposition function listed in

The cholesky function below returns a boolean flag indicating if the algorithm was able to complete the decomposition, which also tells us if the matrix is positive definite or not.

The development of the Cholesky decomposition algorithm comes from the algebra of how \bf{A} is computed as the product of \mathbf{R}^T and \bf{R}.

\mathbf{A} = \mat{r_{11} {} {} {}; r_{12} r_{22} {} {}; r_{13} r_{23} r_{33} {}; r_{14} r_{24} r_{34} r_{44}} \, \mat{r_{11} r_{12} r_{13} r_{14}; {} r_{22} r_{23} r_{24}; {} {} r_{33} r_{34}; {} {} {} r_{44}}

The first row is simple because each element of \bf{A} comes from only one product.

(6.19)\begin{array}{lcl} a_{11} = r_{11}\,r_{11} &\quad &r_{11} = \sqrt{a_{11}} \\ a_{12} = r_{11}\,r_{12} &\quad &r_{12} = a_{12} / r_{11} \\ a_{13} = r_{11}\,r_{13} &\quad &r_{13} = a_{13} / r_{11} \\ a_{14} = r_{11}\,r_{14} &\quad &r_{14} = a_{14} / r_{11} \end{array}

In the subsequent rows, we find the value of the diagonal element using previous found values in the column of the diagonal element.

\begin{array}{lcl} a_{22} = r_{12}^2 + r_{22}^2 &\quad &r_{22} = \sqrt{a_{22} - r_{12}^2} \\ a_{33} = r_{13}^2 + r_{23}^2 + r_{33}^2 &\quad &r_{22} = \sqrt{a_{33} - r_{13}^2 - r_{23}^2} \end{array}

(6.20)r_{ii} = \sqrt{a_{ii} - \sum_{k = 1}^{i-1}r_{ki}^2}, \quad \text{for } i = 2 \ldots n

The algorithm uses equation (6.20) to check if the matrix still appears to be positive definite. If the argument to square root is zero or a negative number, then the algorithm is finished with the matrix marked as not being positive definite. The remaining values of each row are found after finding the diagonal element.

\begin{array}{lcl} a_{23} = r_{12}\,r_{13} + r_{22}\,r_{23} &\quad &r_{23} = \left(a_{23} - r_{12}\,r_{13}\right) / r_{22} \\ a_{34} = r_{13}\,r_{14} + r_{23}\,r_{24} + r_{33}\,r_{34} &\quad &r_{34} = \left(a_{34} - \left(r_{13}\,r_{14} + r_{23}\,r_{24}\right)\right) / r_{33} \end{array}

(6.21)r_{ij} = \left( a_{ij} - \sum_{k = 1}^{i-1}r_{ki}\,r_{kj}\right) / r_{ii}, \quad \text{for } i = 2 \ldots n, \text{ and } j = i+1 \ldots n

The cholesky function uses equations (6.19) to (6.21) to find the Cholesky decomposition. MATLAB has a function called chol that finds the same.

function [R, isPositiveDefinite] = cholesky(A)
% CHOLESKY decomposition. If A is posive definite, A = R'*R.
% isPositiveDefinite is a boolean variable. Use R only if
% isPositiveDefinite is true.

[m,n] = size(A);
if m ~= n
    error("Matrix is not square")
elseif ~issymmetric(A)
    error("Matrix is not positive definite or symmetric")
end
if any(diag(A<=0))
    msg = "Matrix is not positive definite/n";
    msg = strcat(msg, "Diagonal not all positive");
    error(msg);
end
isPositiveDefinite = true;
for i = 1:n
    for k = 1:i-1 % skip i = 1
        A(i,i) = A(i,i) - A(k,i)^2;
    end
    if A(i,i) <= 0
        isPositiveDefinite = false;
        R = 0;
        return;
    end
    A(i,i) = sqrt(A(i,i));
    for j = i+1:n
        for k = 1:i-1
            A(i,j) = A(i,j) - A(k,i)*A(k,j);
        end
        A(i,j) = A(i,j)/A(i,i);
    end
end
R = triu(A);