Cholesky Decomposition Program

A Cholesky Factorization of a real, symmetric, positive-definite, matrix A is the decomposition of the matrix into either
(i) a lower triangular matrix, L, such that A = L * L^T, or
(ii) an upper triangular matrix, U, such that A = U^T * U.

(Note that the terms matrix factorization and matrix decomposition are interchangeable.)

Computing the Lower Triangular Matrix, L

Consider a symmetric 3 x 3 matrix, A, on which we would like to perform a Cholesky decomposition, finding the lower triangular matrix. i.e.,

Expanding the right hand side (RHS) of this equation, it becomes

Working through several of the components, a couple patterns can be detected.

l₁₁ = √ a₁₁

l-components in the same column and below l₁₁ are found by scaling the a-components by l₁₁:

l₂₁ = a₂₁ / l₁₁
l₃₁ = a₃₁ / l₁₁

Moving down to the diagonal in the second row,

l₂₂ = √ (a₂₂ - l₂₁²)

The component below l₂₂ is

l₃₂ = (a₃₂ - l₃₁ l₂₁) / l₂₂

Finally, moving down to the third row:

l₃₃ = √ (a₃₃ - (l₃₂² + l₃₁²))

From this point, a rough algorithm can be developed.
Given an N x N symmetric matrix, A, work through A from the first row (k = 1) to the last row (k = N), first computing the diagonal of that row, and then elements of the column under it.

1) Compute the diagonal element of the k-th row according to the following formula:

2) Compute elements below the diagonal for the k-th column according to the following formula (where i > k):

Computing the Upper Triangular Matrix, U

Repeating the process to find the upper triangular matrix:

Expanding the RHS of this equation, it becomes

Again,

u₁₁ = √ a₁₁

u-components in the same row and to the right of u₁₁ are found by scaling the a-components by u₁₁:

u₁₂ = a₁₂ / u₁₁
u₁₃ = a₁₃ / u₁₁

Moving down to the diagonal in the second row,

u₂₂ = √ (a₂₂ - u₁₂²)

The component to the right of u₂₂ is

u₂₃ = (a₂₃ - u₁₂ u₁₃) / u₂₂

Finally, moving down to the third row

u₃₃ = √ (a₃₃ - (u₁₃² + u₂₃²))

Again, from this point, a rough algorithm can be developed (but with different summations).
Given an N x N symmetric matrix, A, work through A from the first row (k = 1) to the last row (k = N), first computing the diagonal of that row, and then elements of the row to the right of it.

1) Compute the diagonal element of the k-th row according to the following formula:

2) Compute elements to the right of the diagonal for the k-th row according to the following formula (where i > k):

Whichever approach is taken for the decomposition, computing the L matrix or the U matrix, the numbers are the same; the transpose of L is the same as U, and vice versa.

Significance of the Cholesky Decomposition

Among its several applications, a Cholesky Factorization can be used to determine if a symmetric matrix is positive-definite. (The positive-definiteness of a matrix could be checked by computing its eigenvalues: if they are all greater than zero, the matrix is positive-definite. However, computing the eigenvalues of a matrix is computationally a much more expensive process than a Cholesky decomposition.) Given a symmetric matrix, A, we initially assume A is positive-definite and pass it into a Cholesky Decomposition routine. If the routine fails, the initial assumption is proven wrong. On the other hand, if the decomposition succeeds, the positive-definiteness of A is confirmed.

Knowing whether a matrix is positive-definite or not is an important bit of information to have when working in the field of Quadratic Programming (minimizing a specific type of function, F.) When F is stated as a matrix-vector equation, the matrix involved is a symmetric matrix, called the Hessian matrix. If the Hessian matrix is not only symmetric, but also positive-definite, F has a global minimum. If the Hessian matrix is not positive-definite, F does not have a global minimum. Because Quadratic Programming problems can be extremely difficult to solve, knowledge about the existence of a global minimum-- before starting the task of solving for it-- would indicate if efforts would be better spent directed toward a different goal. In other words, if a function is known to not have a global minimum, but does have local minima, it might be desirable to direct the search program to find a local minimum. Parameters and inputs to the Quadratic Programming routine could be modified so that the algorithm is not seeking a global minimum that does not exist.

The Cholesky Decomposition program to which the title of this post hot-links computes U. Note that entry of the full symmetric matrix, A, is expected, even though only the upper triangular half of it is used by the decomposition routine; the strictly lower triangular half of A is not referenced. The symmetry of A is confirmed by the main program before being passed into the decomposition routine. This JavaScript program is a translation of the LAPACK routine DPOTRF.F, which is posted on the NETLIB website. DPOTRF.F includes a multitude of options to customize the program: whether the code is to be blocked or unblocked, whether to compute the upper triangular U or lower triangular L matrix, whether to use loop unrolling, and whether loop increments are the same, to name a few. This JavaScript translation offers none of these options. It computes U. That's it.

A C++ translation of this program is available on GitHub.