### 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

_{11}= √ a

_{11}

l-components in the same column and

*below*l

_{11}are found by scaling the a-components by l

_{11}:

l

_{21}= a

_{21}/ l

_{11}

l

_{31}= a

_{31}/ l

_{11}

Moving down to the diagonal in the second row,

l

_{22}= √ (a

_{22}- l

_{21}

^{2})

The component below l

_{22}is

l

_{32}= (a

_{32}- l

_{31}l

_{21}) / l

_{22}

Finally, moving down to the third row:

l

_{33}= √ (a

_{33}- (l

_{32}

^{2}+ l

_{31}

^{2}))

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

_{11}= √ a

_{11}

u-components in the same row and

*to the right of*u

_{11}are found by scaling the a-components by u

_{11}:

u

_{12}= a

_{12}/ u

_{11}

u

_{13}= a

_{13}/ u

_{11}

Moving down to the diagonal in the second row,

u

_{22}= √ (a

_{22}- u

_{12}

^{2})

The component to the right of u

_{22}is

u

_{23}= (a

_{23}- u

_{12}u

_{13}) / u

_{22}

Finally, moving down to the third row

u

_{33}= √ (a

_{33}- (u

_{13}

^{2}+ u

_{23}

^{2}))

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: calculator, cholesky, decomposition, factorization, math, mathematics

## 0 Comments:

Post a Comment

<< Home