Editing Cholesky decomposition (section)

=== The Cholesky–Banachiewicz and Cholesky–Crout algorithms ===
[[File:Chol.gif|thumb|Access pattern (white) and writing pattern (yellow) for the in-place Cholesky—Banachiewicz algorithm on a 5×5 matrix]]
If the equation
<math display=block>\begin{align}
\mathbf{A} = \mathbf{LL}^T & =
\begin{pmatrix}   L_{11} & 0 & 0 \\
   L_{21} & L_{22} & 0 \\
   L_{31} & L_{32} & L_{33}\\
\end{pmatrix}
\begin{pmatrix}   L_{11} & L_{21} & L_{31} \\
   0 & L_{22} & L_{32} \\
   0 & 0 & L_{33}
\end{pmatrix} \\[8pt]
& =
\begin{pmatrix}   L_{11}^2 &   &(\text{symmetric})   \\
   L_{21}L_{11} & L_{21}^2 + L_{22}^2& \\
   L_{31}L_{11} & L_{31}L_{21}+L_{32}L_{22} & L_{31}^2 + L_{32}^2+L_{33}^2
\end{pmatrix},
\end{align}</math>

is written out, the following is obtained:

<math display=block>\begin{align}
\mathbf{L} = 
\begin{pmatrix} \sqrt{A_{11}} &  0 & 0  \\
A_{21}/L_{11} & \sqrt{A_{22} - L_{21}^2} & 0 \\
A_{31}/L_{11} &  \left( A_{32} - L_{31}L_{21} \right) /L_{22}  &\sqrt{A_{33}- L_{31}^2 - L_{32}^2}
\end{pmatrix}
\end{align}</math>

and therefore the following formulas for the entries of {{math|'''L'''}}:

<math display=block> L_{j,j} = (\pm)\sqrt{ A_{j,j} - \sum_{k=1}^{j-1} L_{j,k}^2 }, </math>
<math display=block> L_{i,j} = \frac{1}{L_{j,j}} \left( A_{i,j} - \sum_{k=1}^{j-1} L_{i,k} L_{j,k} \right) \quad \text{for } i>j. </math>

For complex and real matrices, inconsequential arbitrary sign changes of diagonal and associated off-diagonal elements are allowed. The expression under the [[square root]] is always positive if {{math|'''A'''}} is real and positive-definite. 

For complex Hermitian matrix, the following formula applies:

<math display=block> L_{j,j} = \sqrt{ A_{j,j} - \sum_{k=1}^{j-1} L_{j,k}^*L_{j,k} }, </math>
<math display=block> L_{i,j} = \frac{1}{L_{j,j}} \left( A_{i,j} - \sum_{k=1}^{j-1} L_{j,k}^* L_{i,k} \right) \quad \text{for } i>j. </math>

So it now is possible to compute the {{math|(''i'', ''j'')}} entry if the entries to the left and above are known. The computation is usually arranged in either of the following orders:
* The '''Cholesky–Banachiewicz algorithm''' starts from the upper left corner of the matrix {{mvar|L}} and proceeds to calculate the matrix row by row. 
<syntaxhighlight lang="C">
for (i = 0; i < dimensionSize; i++) {
    for (j = 0; j <= i; j++) {
        float sum = 0;
        for (k = 0; k < j; k++)
            sum += L[i][k] * L[j][k];

        if (i == j)
            L[i][j] = sqrt(A[i][i] - sum);
        else
            L[i][j] = (1.0 / L[j][j] * (A[i][j] - sum));
    }
}
</syntaxhighlight>
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,
<syntaxhighlight lang="Fortran">
do i = 1, size(A,1)
    L(i,i) = sqrt(A(i,i) - dot_product(L(i,1:i-1), L(i,1:i-1)))
    L(i+1:,i) = (A(i+1:,i) - matmul(conjg(L(i,1:i-1)), L(i+1:,1:i-1))) / L(i,i)
end do
</syntaxhighlight>
where <code>conjg</code> refers to complex conjugate of the elements.

* The '''Cholesky–Crout algorithm''' starts from the upper left corner of the matrix {{mvar|L}} and proceeds to calculate the matrix column by column. <syntaxhighlight lang="C">
for (j = 0; j < dimensionSize; j++) {
    float sum = 0;
    for (k = 0; k < j; k++) {
        sum += L[j][k] * L[j][k];
    }
    L[j][j] = sqrt(A[j][j] - sum);

    for (i = j + 1; i < dimensionSize; i++) {
        sum = 0;
        for (k = 0; k < j; k++) {
            sum += L[i][k] * L[j][k];
        }
        L[i][j] = (1.0 / L[j][j] * (A[i][j] - sum));
    }
}
</syntaxhighlight>
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,
<syntaxhighlight lang="Fortran">
do i = 1, size(A,1)
    L(i,i) = sqrt(A(i,i) - dot_product(L(1:i-1,i), L(1:i-1,i)))
    L(i,i+1:) = (A(i,i+1:) - matmul(conjg(L(1:i-1,i)), L(1:i-1,i+1:))) / L(i,i)
end do
</syntaxhighlight>
where <code>conjg</code> refers to complex conjugate of the elements.

Either pattern of access allows the entire computation to be performed in-place if desired.