--- title: "Extended Cholesky decomposition in R" author: "Stéphane Laurent" date: "2017-12-03" output: md_document: variant: markdown preserve_yaml: true html_document: keep_md: no tags: R, maths highlighter: "pandoc-solarized" --- {r setup, include=FALSE} knitr::opts_chunk$set(collapse = TRUE)  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. {r, attr.source='.numberLines'} 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: {r, attr.source='.numberLines'} 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)) #~ 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) ##~~ 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)) #~ 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) ##~~ 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)) #~ 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)