David's Blog

Living a quiet life in Coquitlam, B.C.

Name:
Location: Coquitlam, British Columbia, Canada

Sunday, September 17, 2017

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 * LT, or
(ii) an upper triangular matrix, U, such that A = UT * 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.

l11 = √ a11

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

l21 = a21 / l11
l31 = a31 / l11

Moving down to the diagonal in the second row,

l22 = √ (a22 - l212)

The component below l22 is

l32 = (a32 - l31 l21) / l22

Finally, moving down to the third row:

l33 = √ (a33 - (l322 + l312))

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,

u11 = √ a11

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

u12 = a12 / u11
u13 = a13 / u11

Moving down to the diagonal in the second row,

u22 = √ (a22 - u122)

The component to the right of u22 is

u23 = (a23 - u12 u13) / u22

Finally, moving down to the third row

u33 = √ (a33 - (u132 + u232))

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.

Labels: , , , , ,

0 Comments:

Post a Comment

<< Home