Cholesky decomposition.

The Cholesky decomposition can be applied on any square matrix $A$ which is symmetric and positive definite. The Cholesky decomposition is one of the fastest decomposition available. It constructs a lower triangular matrix L which has the following property:

$$L \cdot L^T = A$$ (13.33)

This factorization is sometimes referred to, as ``taking the square root'' of the matrix $A$.
The Cholesky decomposition is a particular case of the $LU$ decomposition. The $LU$ decomposition is the following:

$$L \cdot U = A$$ (13.34)

where $L$ is a lower triangular matrix and $U$ is a upper triangular matrix. For example, in the case of a $4 \times 4$ matrix $A$, we have:

$$\begin{pmatrix} \alpha_{11} & 0 & 0 & 0 \\ \alpha_{21} & \alp... ...a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{pmatrix}$$ (13.35)

We can use the $LU$ decomposition to solve the linear set: $Ax=B \Leftrightarrow (LU)x=B$ by first solving for the vector $y$ such that $Ly=B$ and then solving $Ux=y$. These two systems are trivial to solve because they are triangular.

Performing $LU$ decomposition.

First let us rewrite the component $a_{ij}$ of $A$ from the equation 13.34 or 13.35. That component is always a sum beginning with

$$a_{i,j} = \alpha_{i1} \beta_{1j} + \cdots$$

The number of terms in the sum depends, however on wether $i$ or $j$ is the smallest number. We have, in fact three cases:

$$i<j: a_{ij} = \alpha_{i1} \beta_{1j} + \alpha_{i2} \beta_{2j} + \cdots + \alpha_{ii} \beta_{ij}$$ (13.36)

$$i=j: a_{ii} = \alpha_{i1} \beta_{1i} + \alpha_{i2} \beta_{2i} + \cdots + \alpha_{ii} \beta_{ii}$$ (13.37)

$$i>j: a_{ij} = \alpha_{i1} \beta_{1j} + \alpha_{i2} \beta_{2j} + \cdots + \alpha_{ij} \beta_{jj}$$ (13.38)

Equations 13.36 - 13.38, totalize $n^2$ equations for the $n^2+n$ unknown $\alpha$'s and $\beta$'s (the diagonal being represented twice). Since the number of unknowns is greater than the number of equations, we have to specify $n$ of the unknowns arbitrarily and then solve for the others. In fact, as we shall see, it is always possible to take:

$$\alpha_{ii}\equiv 1 \quad i=1, \ldots ,n$$ (13.39)

A surprising procedure, now, is Crout's algorithm, which, quite trivially, solves the set of $n^2+n$ Equations 13.36 - 13.38 for all the $\alpha$'s and $\beta$'s by just arranging the equation in a certain order! That order is as follows:

• Set $\alpha_{ii}=1, \quad i=1, \ldots,n$
• For each $j=1, \ldots, n$ do these two procedures:
  First, for $i=1, \ldots, j$ use 13.36, 13.37 and 13.39 to solve for $\beta_{ij}$, namely

  $$\beta_{ij}=a_{ij}- \sum_{k=1}^{i-1} \alpha_{ik} \beta_{kj}$$ (13.40)

  
  
  Second, for $i=j+1, \ldots, n$, use 13.38 to solve for $\alpha_{ij}$, namely

  $$\alpha_{ij}= \frac{1}{\beta_{jj}} \Bigg( a_{ij} - \sum_{k=1}^{j-1} \alpha_{ik} \beta_{kj} \Bigg)$$ (13.41)

  
  
  Be sure to do both procedures before going on to the next j.

If you work through a few iterations of the above procedure, you will see that the $\alpha$'s and $\beta$'s that occur on the right-hand side of the Equation 13.40 and 13.41 are already determined by the time they are needed.

Performing Cholesky decomposition.

We can obtain the analogs of Equations 13.40 and 13.41 for the Cholesky decomposition:

$$L_{ii}=\sqrt{a_{ii}- \sum_{k=1}^{i-1} L_{ik}^2}$$ (13.42)

and

$$L_{ji}= \frac{1}{L_{ii}} \Bigg( a_{ij} - \sum_{k=1}^{i-1} L_{ik} L_{jk} \Bigg) \quad j=i+1, \ldots, n$$ (13.43)

If you apply Equation 13.42 and 13.43 in the order $i=1, \ldots,n$, you will see the the $L$'s that occur on the right-hand side are exactly determined by the time they are needed. Also, only components $a_{ij}$ with $j>i$ are referenced.
If the matrix $A$ is not positive definite, the algorithm will stop, trying to take the square root of a negative number in equation 13.42.
What about pivoting? Pivoting (i.e., the selection of a salubrious pivot element for the division in Equation 13.43) is not really required for the stability of the algorithm. In fact, the only cause of failure is if the matrix $A$ (or, with roundoff error, another very nearby matrix) is not positive definite.