§ TIL: Cholesky decomposition

Not long ago I picked up a Matrix Analysis textbook and came across a matrix factorization with this entry’s namesake. The Cholesky decomposition factorizes a (Hermitian) positive-definite matrix \(\mathbf{A}\) (\( x^{T} \mathbf{A} x > 0 \) for all non-zero \(x \in \mathbb{R}^{d}\)) into a product of a lower triangular matrix \(\mathbf{L}\) and its adjoint:

\[ \mathbf{A} = \mathbf{L}\mathbf{L}^{\dagger}. \]

For a positive-definite matrix such a factorization always exists and is unique!

The numerical algorithm for constructing the lower triangular matrix \(\mathbf{L}\) is relatively straightforward. Initially, the matrix \(\mathbf{L}\) begins as a zero matrix with appropriate dimensions (say \(N \times N\)). Additionally, we take the diagonal elements of \(\mathbf{A}\) (\(A_{jj}\)) and arrange them in a non-increasing order then algorithm proceeds as follows:

  1. At step \(j\), update the diagonal elements of \(\mathbf{L}\) with those of \(\mathbf{A}\), \(L_{jj} \leftarrow A_{jj}\).
  2. Update the \(L_{ij}\), \(L_{ij} \leftarrow A_{ij} / L_{jj}\) for \(i = j+1, \ldots, N \).
  3. Update the stored diagonal elements of \(A_{ii}\), \(A_{ii} \leftarrow A_{ii} - L_{ij}^2\) for \(i = j+1, \ldots, N\)
  4. Increment \(j\) and go back to step 1 and repeat until \(j\) is equal to index of the last row of the matrix \(\mathbf{A}\).

§ Extras

  • The above algorithm can be slightly modified to obtain a low rank approximation of \(\mathbf{A}\): we can stop when all the remaining updated diagonal elements \(A_{jj}\) are below some specified tolerance \(\delta\). For some instances of \(\mathbf{A}\), e.g., electron repulsion integrals, the rank \(\nu\) of the resulting matrix after this modification can be much smaller than its maximum possible value!

  • A Python implementation of the algorithm can fit in a tweet!

import numpy as np

def cholesky(A):
  m, _ = A.shape
  for i in range(m):
      A[i,i] = np.sqrt(A[i,i])
      A[i+1:, i] /= A[i, i]
      A[i, i+1:] = 0
      A[i+1:, i+1:] = np.outer(A[i+1:, i], A[i+1:,i])
  • Animated GIF of the Cholesky decomposition of a \(20 \times 20\) matrix:
Animation of cholesky decomposition of \\(20\times20\\) matrix