Cholesky decomposition
In linear algebra, the Cholesky decomposition or Cholesky factorization (pronounced /ʃəˈlɛski/ shə-LES-kee) is a decomposition of a Hermitian, positive-definite matrix into the product of a lower triangular matrix and its conjugate transpose, which is useful for efficient numerical solutions, e.g., Monte Carlo simulations.
[3] The converse holds trivially: if A can be written as LL* for some invertible L, lower triangular or otherwise, then A is Hermitian and positive definite.
However, if the rank of A is r, then there is a unique lower triangular L with exactly r positive diagonal elements and n − r columns containing all zeroes.
The main advantage is that the LDL decomposition can be computed and used with essentially the same algorithms, but avoids extracting square roots.
decomposition exists where the number of non-zero elements on the diagonal D is exactly the rank of A.
For linear systems that can be put into symmetric form, the Cholesky decomposition (or its LDL variant) is the method of choice, for superior efficiency and numerical stability.
[2] Systems of the form Ax = b with A symmetric and positive definite arise quite often in applications.
[14] The Cholesky decomposition is commonly used in the Monte Carlo method for simulating systems with multiple correlated variables.
[15] The following simplified example shows the economy one gets from the Cholesky decomposition: suppose the goal is to generate two correlated normal variables
Unscented Kalman filters commonly use the Cholesky decomposition to choose a set of so-called sigma points.
The columns of L can be added and subtracted from the mean x to form a set of 2N vectors called sigma points.
These sigma points completely capture the mean and covariance of the system state.
The explicit inverse of a Hermitian matrix can be computed by Cholesky decomposition, in a manner similar to solving linear systems, using
Generally, the first algorithm will be slightly slower because it accesses the data in a less regular manner.
The Cholesky algorithm, used to calculate the decomposition matrix L, is a modified version of Gaussian elimination.
Note that bi b*i is an outer product, therefore this algorithm is called the outer-product version in (Golub & Van Loan).
This is repeated for i from 1 to n. After n steps, A(n+1) = I is obtained, and hence, the lower triangular matrix L sought for is calculated as
For complex and real matrices, inconsequential arbitrary sign changes of diagonal and associated off-diagonal elements are allowed.
The computation is usually arranged in either of the following orders: The above algorithm can be succinctly expressed as combining a dot product and matrix multiplication in vectorized programming languages such as Fortran as the following, where conjg refers to complex conjugate of the elements.
The above algorithm can be succinctly expressed as combining a dot product and matrix multiplication in vectorized programming languages such as Fortran as the following, where conjg refers to complex conjugate of the elements.
Either pattern of access allows the entire computation to be performed in-place if desired.
Specifically, if Ax = b, and y denotes the computed solution, then y solves the perturbed system (A + E)y = b, where
Here ||·||2 is the matrix 2-norm, cn is a small constant depending on n, and ε denotes the unit round-off.
[17] While this might lessen the accuracy of the decomposition, it can be very favorable for other reasons; for example, when performing Newton's method in optimization, adding a diagonal matrix can improve stability when far from the optimum.
Again, the pattern of access allows the entire computation to be performed in-place if desired.
This involves matrix products and explicit inversion, thus limiting the practical block size.
Notice that the equations above that involve finding the Cholesky decomposition of a new matrix are all of the form
, which allows them to be efficiently calculated using the update and downdate procedures detailed in the previous section.
The argument is not fully constructive, i.e., it gives no explicit numerical algorithms for computing Cholesky factors.
The Cholesky factorization can be generalized [citation needed] to (not necessarily finite) matrices with operator entries.
