Suma residual de cuadrados de matriz de bloque con elementos dispersos.

Question

Basado en la pregunta de la Suma Residual de cuadrados de regresión Ponderada, una forma rápida de resolver para $$(\mathbf{y-X\hat{\boldsymbol\beta}})^{'}\mathbf{C}^{-1}(\mathbf{y-X\hat{\boldsymbol\beta}})$$ es la transformación de los probelm a una regresión lineal a través del siguiente procedimiento

1) La descomposición de cholesky de a $\mathbf{C}$ $\mathbf{C=RR^{'}}$ $\mathbf{C^{-1}=R^{'-1}R^{-1}}$

2) La estiamtor de ${\boldsymbol\beta}$ puede ser escrita como : $$\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) Resolución de $\mathbf{R^{-1}X=A}$ $\mathbf{R^{-1}y=B}$ con backsolve llegamos $$\hat{\boldsymbol\beta}=(\mathbf{A^{'}A})^{-1}\mathbf{A^{'}B}$$, que es equivalente a la suma residual de los cuadrados (RSS) de no ponderado de regresión.

En $3)$ el backsolve operador es realmente más rápido que la inversión de una matriz y de problemas para el ponderado estimador de regresión $\hat{\boldsymbol\beta}$ se puede hacer muy rápido utilizando el $lm$ función

Mi pregunta está relacionada con el caso al $\mathbf{C}$ es un simétrica positiva definida la matriz de la forma

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

donde las matrices $A,B,C,D,$ son simétrica positiva definida matrices. Invering esta matriz puede ser hecho usando el schurr complemento como en la pregunta Inversa de bloque de la matriz de covarianza. Sin embargo, utilizando el complemento de schur para encontrar $\mathbf{C^{-1}}$ y luego resolver para $(\mathbf{y-X\hat{\boldsymbol\beta}})^{'}\mathbf{C}^{-1}(\mathbf{y-X\hat{\boldsymbol\beta}})$ donde $\hat{\boldsymbol\beta}=(\mathbf{X^{'}C^{-1}X})^{-1}\mathbf{X^{'}C^{-1}y}$ parece mucho menos eficiente que la descomposición de cholesky procedimiento.

Hay una manera de simplificar aún más la descomposición de cholesky, el procedimiento de búsqueda de la RSS dado que la matriz de $\mathbf{C}$ es de la forma anterior ? Por ejemplo, ¿hay una manera eficaz para calcular la descomposición de cholesky de la matriz dispersa $\mathbf{C}$ ?

Answer 1

La punta de flecha y el bloque de punta de flecha estructuras se conservan en La descomposición de Cholesky. Así podemos encontrar el menor triangular Cholesky de la raíz $\mathbf{L}$ de la ampliación de la matriz dispersa $\mathbf{C}_0 = \mathbf{L}\mathbf{L}'$ bajo la forma $$ \mathbf{L} = \begin{bmatrix} \mathbf{L}_A & & & \\ & \mathbf{L}_B & & \\ & & \mathbf{L}_C & \\ \mathbf{U}' & \mathbf{V}' & \mathbf{W}' & \mathbf{H} \end{bmatrix} $$ en los bloques que se muestran no son ceros. Obviamente $\mathbf{L}_A$, $\mathbf{L}_B$ y $\mathbf{L}_C$ son los Cholesky raíces de la diagonal de bloques $\mathbf{A}$, $\mathbf{B}$ y $\mathbf{C}$. Mediante la identificación de los tres primeros elementos de la última línea de bloques obtenemos $$ \mathbf{U}'\mathbf{L}_A' = \mathbf{E}', \quad \mathbf{V}'\mathbf{L}_B' = \mathbf{F}', \quad \mathbf{W}'\mathbf{L}_C' = \mathbf{G}' $$ el que da las matrices $\mathbf{U}$, $\mathbf{V}$ y $\mathbf{W}$, por ejemplo $\mathbf{U} = \mathbf{L}_A^{-1}\mathbf{E}$ podrían ser obtenidos por forwardsolve(LA, E) en R. A continuación, identificar el cuarto elemento de bloque del último bloque de la fila da $$ \mathbf{U}'\mathbf{U} + \mathbf{V}'\mathbf{V} + \mathbf{W}'\mathbf{W} + \mathbf{H}\mathbf{H}' = \mathbf{D}. $$ Por lo $\mathbf{H}$ puede ser la transposición de la Cholesky raíz de $\mathbf{L}_S$ de la matriz $$ \mathbf{S} := \mathbf{D} - \mathbf{U}'\mathbf{U} - \mathbf{V}'\mathbf{V} - \mathbf{W}'\mathbf{W}. $$ Para resumir, en R el cálculo de $\mathbf{L}$ tomaría la forma de tres forwardsolve encontrar $\mathbf{U}$, $\mathbf{V}$ y $\mathbf{W}$, tres crossprod y un resta formulario de $\mathbf{S}$, chol. Por supuesto, si usted sólo desea resolver un sistema lineal no es necesario para formar la gran Cholesky de la raíz, sino para resolver sistemas el uso de forwardsolve en el bloque sub-vectores.

Un problema puede surgir cuando el condicionamiento de la matriz grande $\mathbf{C}_0$ no es suficiente. Es podría ser el caso, la matriz de $\mathbf{S}$ no sería numéricamente positiva en definitiva.

## 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))