Procedure: Pivoted Cholesky decomposition¶
Description and Background¶
The Cholesky decomposition is named after Andre-Louis Cholesky, who found that a symmetric positive-definite matrix can be decomposed into a lower triangular matrix and the transpose of the lower triangular matrix. Formally, the Cholesky method decomposes a symmetric positive definite matrix \(A\) uniquely into a product of a lower triangular matrix \(L\) and its transpose, i.e. \(A = LL^T\), or equivalently \(A = R^TR\) where \(R\) is a upper triangular matrix.
The Pivoted Cholesky decomposition, or the Cholesky decomposition with complete pivoting, of a matrix \(A\) returns a permutation matrix \(P\) and the unique upper triangular matrix \(R\) such that \(P^TAP = R^TR\). The permutation matrix is an orthogonal matrix, so the matrix \(A\) can be rewritten as \(A = (RP^T)^T RP^T\)
An arbitrary permutation of rows and columns of matrix \(A\) can be decomposed by the Cholesky algorithm, \(P^TAP=R^TR\), where \(P\) is a permutation matrix and \(R\) is an upper triangular matrix. The permutation matrix is an orthogonal matrix, so the matrix \(A\) can be rewritten as \(A = P R^T RP^T\), so that \(C =P R^T\) is the pivoted Cholesky decomposition of \(A\).
Pivoted Cholesky algorithm: This algorithm computes the Pivoted Cholesky decomposition \(P^TAP=R^TR\) of a symmetric positive semidefinite matrix \(A \in \Re^{n \times n}\). The nonzero elements of the permutation matrix \(P\) are given by \(P(piv(k),k)=1\), \(k=1,\ldots,n\)
More details about the numerical analysis of the pivoted Cholesky decomposition can be found in Higham (2002).
Inputs¶
- Symmetric positive-definite matrix \(A\)
Outputs¶
- Permutation matrix \(P\)
- Unique upper triangular matrix \(R\)
Procedure¶
Initialise
\(R = \mathrm{zeros}(\mathrm{size}(A))\)
\(piv(k) = k,\, \forall k\in[1,n]\)
Repeat for \(k=1:n\)
{Finding the pivot}\(\quad B = A(k:n,k:n)\)\(\quad l = \left\{ i : A(i,i) == \max diag \left( B \right) \right\}\){Swap rows and columns}\(\quad A(:,k) <=> A(:,l)\)\(\quad R(:,k) <=> R(:,l)\)\(\quad A(k,:) <=> A(l,:)\)\(\quad piv(k) <=> piv(l)\){Cholesky decomposition}\(\quad R(k,k) = \sqrt{A(k,k)}\)\(\quad R(k,k+1:n) = R(k,k)^{-1} A(k,k+1:n)\){Updating \(A\)}\(\quad A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - R(k,k+1:n)^TR(k,k+1:n)\)
End repeat
Additional Comments¶
If \(A\) is a variance matrix, the pivoting order given by the permutation matrix \(P\) has the following interpretation: the first pivot is the index of the element with the largest variance, the second pivot is the index of the element with the largest variance conditioned on the first element, the third pivot is the index of the element with the largest variance conditioned on the first two elements element, and so on.
References¶
- Higham, N. J. Accuracy and Stability of Numerical Algorithms. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2002. ISBN 0-89871-521-0. Second Edition.