Residual sum of squares of block matrix with sparse elements

Question

Based on the question Residual Sum of squares in Weighted regression, a fast way to solve for $$(\mathbf{y-X\hat{\boldsymbol\beta}})^{'}\mathbf{C}^{-1}(\mathbf{y-X\hat{\boldsymbol\beta}})$$ is to transform the probelm to a regular linear regression through the following procedure

1) The cholesky decomposition of $\mathbf{C}$ is $\mathbf{C=RR^{'}}$ and $\mathbf{C^{-1}=R^{'-1}R^{-1}}$

2) The estiamtor of ${\boldsymbol\beta}$ can be written as : $$\hat{\boldsymbol\beta}=(\mathbf{X^{'}C^{-1}X})^{-1}\mathbf{X^{'}C^{-1}y}=(\mathbf{X^{'}R^{'-1}R^{-1}X})^{-1}\mathbf{X^{'}R^{'-1}R^{-1}y}$$

3) Solving $\mathbf{R^{-1}X=A}$ and $\mathbf{R^{-1}y=B}$ with backsolve we get $$\hat{\boldsymbol\beta}=(\mathbf{A^{'}A})^{-1}\mathbf{A^{'}B}$$ which is equivalent to the residual sum of squares (RSS) of unweighted regression.

In $3)$ the backsolve operator is really faster than inverting a matrix and solving for the unweighted regression estimator $\hat{\boldsymbol\beta}$ can be done really fast using the $lm$ function

My question is related to the case when $\mathbf{C}$ is a positive definite symmetric matrix of the form

$$\pmatrix{A & 0 & 0 & E \\ 0 & B & 0 & F \\ 0 & 0 & C & G \\ E^\prime & F^\prime & G^\prime & D}$$

where matrices $A,B,C,D,$ are positive definite symmetric matrices. Invering this matrix can be done using the schurr complement as in the question Inverse of block covariance matrix. However using the schur complement to find $\mathbf{C^{-1}}$ and then solving for $(\mathbf{y-X\hat{\boldsymbol\beta}})^{'}\mathbf{C}^{-1}(\mathbf{y-X\hat{\boldsymbol\beta}})$ where $\hat{\boldsymbol\beta}=(\mathbf{X^{'}C^{-1}X})^{-1}\mathbf{X^{'}C^{-1}y}$ seems far less efficient then the cholesky decomposition procedure.

Is there a way to further simplify the cholesky decomposition procedure when finding the RSS given that the matrix $\mathbf{C}$ is of the form above ? For instance is there an efficient way to calculate the cholesky decomposition of the sparse matrix $\mathbf{C}$ ?

Answer 1

The arrowhead and block arrowhead structures are preserved in Cholesky decomposition. So we can find the lower triangular Cholesky root $\mathbf{L}$ of the large sparse matrix $\mathbf{C}_0 = \mathbf{L}\mathbf{L}'$ under the form $$ \mathbf{L} = \begin{bmatrix} \mathbf{L}_A & & & \\ & \mathbf{L}_B & & \\ & & \mathbf{L}_C & \\ \mathbf{U}' & \mathbf{V}' & \mathbf{W}' & \mathbf{H} \end{bmatrix} $$ where the blocks not shown are zeros. Obviously $\mathbf{L}_A$, $\mathbf{L}_B$ and $\mathbf{L}_C$ are the Cholesky roots of the diagonal blocks $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{C}$. By identifying the first three elements of the last block line we get $$ \mathbf{U}'\mathbf{L}_A' = \mathbf{E}', \quad \mathbf{V}'\mathbf{L}_B' = \mathbf{F}', \quad \mathbf{W}'\mathbf{L}_C' = \mathbf{G}' $$ which gives the matrices $\mathbf{U}$, $\mathbf{V}$ and $\mathbf{W}$, for instance $\mathbf{U} = \mathbf{L}_A^{-1}\mathbf{E}$ would be obtained by forwardsolve(LA, E) in R. Then identifying the fourth block element of the last block row gives $$ \mathbf{U}'\mathbf{U} + \mathbf{V}'\mathbf{V} + \mathbf{W}'\mathbf{W} + \mathbf{H}\mathbf{H}' = \mathbf{D}. $$ So $\mathbf{H}$ can be the transpose of the Cholesky root $\mathbf{L}_S$ of the matrix $$ \mathbf{S} := \mathbf{D} - \mathbf{U}'\mathbf{U} - \mathbf{V}'\mathbf{V} - \mathbf{W}'\mathbf{W}. $$ To summarise, in R the computation of $\mathbf{L}$ would take the form of three forwardsolve to find $\mathbf{U}$, $\mathbf{V}$ and $\mathbf{W}$, three crossprod and a subtraction to form $\mathbf{S}$, then a chol. Of course if you only want to solve a linear system it is not necessary to form the large Cholesky root, but rather to solve systems using forwardsolve on block sub-vectors.

A problem may arise when the conditioning of the large matrix $\mathbf{C}_0$ is not good enough. It could be the case then that the matrix $\mathbf{S}$ would not be numerically positive definite.

## Change this for more tests 'n' is the size of diagonal blocks
set.seed(314159)
n <- c("A" = 2, "B" = 2, "C" = 2, "D" = 3)

## origin and end of the blocks
nB <- length(n); n0 <- sum(n)
to <- cumsum(n)
from <- 1 + c(0, to[1:(nB-1)])
names(from) <- names(n)

## build a block-triangular (but not triangular) 'R' root of 'C0' with
## the wanted structure, then get 'C0' which will be symmetric and
## positive definite by construction, and also will have the wanted
## structure.

R <- array(0, dim = c(n0, n0))

for (i in 1:nB) {
    Ri <- array(runif(n[i] * n[i]), dim = c(n[i], n[i]))
    ind <- from[i]:to[i]
    R[ind, ind] <- Ri
}

## 'indLast' contains the indices of the last block

indLast <- ind
R[indLast, 1:to[nB - 1]] <- runif(n[nB] * to[nB - 1])

## get 'C0'

C0 <- R %*% t(R)

## now retrieve the Cholesky root 'L0' of 'C0' as suggested (could be
## done with a sparse matrix as well)
## o 'C0[ind, indLast]' will succesively be 'E', 'F' and 'G'
## o The matrix 'Lind' will successively contain 'L_A', 'L_B', 'L_C'
## o The matrix 'Mind' will successively contain 'E', 'F', 'G' which
##   could better have been named say 'M_1', 'M_2', 'M_3' ...

D <- C0[indLast, indLast]
L0 <- array(0, dim = c(n0, n0))

for (i in 1:(nB-1)) {
    ind <- from[i]:to[i]
    Lind <- t(chol(C0[ind, ind]))
    Mind <- forwardsolve(Lind, C0[ind, indLast])
    D <- D - crossprod(Mind, Mind)
    L0[ind, ind] <- Lind
    L0[indLast, ind] <- t(Mind)
}
L0[indLast, indLast] <- t(chol(D))

## finally compare with the "true" value"

L0.true <- t(chol(C0))
max(abs(L0 - L0.true))