Extended Cholesky decomposition in R

Posted on December 3, 2017 by Stéphane Laurent

Tags: R, maths

Let \(S\) be a symmetric positive semidefinite matrix of order \(d\) having rank \(r\). An extended Cholesky decomposition of \(S\) is a triplet \((L,M,P)\) consisting of a lower triangular \(r\times r\)-matrix \(L\), a \((d-r) \times r\)-matrix \(M\), and a permutation matrix \(P\) of order \(d\) such that, setting \[ C = \begin{pmatrix} L & 0 \\ M & 0 \end{pmatrix}, \] one has \(PSP' = CC'\). Besides, setting \[ \widetilde{C} = \begin{pmatrix} L & 0 \\ M & I_{d-r} \end{pmatrix}, \] one has \(S={(\widetilde{C}'P)}'I_d^r \widetilde{C}'P\) where \(I_d^r\) is the \(d\times d\)-matrix \(\begin{pmatrix} I_r & 0 \\ 0 & 0 \end{pmatrix}\).

The R function below calculates an extended Cholesky decomposition.

extendedCholesky <- function(S){
  C <- suppressWarnings(chol(S, pivot=TRUE))
  d <- nrow(C)
  P <- matrix(0, d, d)
  P[cbind(1:d, attr(C,"pivot"))] <- 1
  r <- attr(C, "rank")
  return(list(L = t(C[seq_len(r), seq_len(r), drop=FALSE]), 
              M = t(C[seq_len(r), seq_len(d-r)+r, drop=FALSE]), 
              P = P))
}

Let’s check:

d <- 3
##~~ check for a rank 1 matrix ~~##
S <- tcrossprod(c(1:d))
#~ extended Cholesky of S ~#
EC <- extendedCholesky(S); P <- EC$P; L <- EC$L; M <- EC$M
#~ C matrix ~#
C <- cbind(rbind(L,M), matrix(0, d, d-ncol(L)))
all.equal(P %*% S %*% t(P), C%*%t(C))
## [1] TRUE

#~ C tilde matrix ~#
Ctilde <- cbind(rbind(L,M), 
                rbind(matrix(0, nrow(L), d-nrow(L)), diag(d-nrow(L))))
all.equal(
  t(t(Ctilde)%*%P) %*% 
    diag(c(rep(1, nrow(L)), rep(0, d-nrow(L)))) %*% 
    (t(Ctilde)%*%P), 
  S)
## [1] TRUE

##~~ check for a rank 2 matrix ~~##
S <- tcrossprod(c(1:d)) + tcrossprod(d:1)
#~ extended Cholesky of S ~#
EC <- extendedCholesky(S); P <- EC$P; L <- EC$L; M <- EC$M
#~ C matrix ~#
C <- cbind(rbind(L,M), matrix(0, d, d-ncol(L)))
all.equal(P %*% S %*% t(P), C%*%t(C))
## [1] TRUE

#~ C tilde matrix ~#
Ctilde <- cbind(rbind(L,M), 
                rbind(matrix(0, nrow(L), d-nrow(L)), diag(d-nrow(L))))
all.equal(
  t(t(Ctilde)%*%P) %*% 
    diag(c(rep(1, nrow(L)), rep(0, d-nrow(L)))) %*% 
    (t(Ctilde)%*%P), 
  S)
## [1] TRUE

##~~ check for a rank 3 matrix ~~##
S <- toeplitz(d:1)
#~ extended Cholesky of S ~#
EC <- extendedCholesky(S); P <- EC$P; L <- EC$L; M <- EC$M
#~ C matrix ~#
C <- cbind(rbind(L,M), matrix(0, d, d-ncol(L)))
all.equal(P %*% S %*% t(P), C%*%t(C))
## [1] TRUE

#~ C tilde matrix ~#
Ctilde <- cbind(rbind(L,M), 
                rbind(matrix(0, nrow(L), d-nrow(L)), diag(d-nrow(L))))
all.equal(
  t(t(Ctilde)%*%P) %*% 
    diag(c(rep(1, nrow(L)), rep(0, d-nrow(L)))) %*% 
    (t(Ctilde)%*%P), 
  S)
## [1] TRUE