Invertir una matriz de covarianza escasa

Question

Tengo un postive definitiva simétrica matriz de covarianza que se parece a esto:

Tenga en cuenta que todos los a,B,C,D,E,F,G son también poitive definitiva simétrica de las matrices de covarianza quiero encontrar una manera fácil en que se puede invertir la matriz de partes en lugar de invertir todo en un momento. Realmente agradecería un código R o de cualquier método de R que puede resolver mi problema Mi objetivo principal es encontrar una eficiente-de forma rápida, para invertir esta matriz y encontrar su determinante

Answer 1

Podría usted hacer uso de un simple bloque de la matriz de inversión?

$$\left[\begin{array}{cc} A & B \\C & D \end{array}\right]^{-1}$$

$$=\begin{bmatrix} (A - BD^{-1}C)^{-1} & -A^{-1}B(D - CA^{-1}B)^{-1} \\ -D^{-1}C(A - BD^{-1}C)^{-1} & (D - CA^{-1}B)^{-1} \end{bmatrix}$$

(donde aquí "$A$" probablemente correspondan mejor a su matriz diagonal con o bloque-diagonal $(A,B,C)$)

Si la covarianza es realmente ordinaria de varianza-covarianza de la matriz, hay más adecuada (aunque relacionados) los cálculos. También podría considerar la matriz de descomposición, tales como Choleski, que se simplifica en gran medida por la estructura particular que usted tiene.

--

user11852 la convierte en un excelente punto en los comentarios; hay muchas situaciones en las que usted realmente no necesita para calcular la inversa de la misma.

Answer 2

Para el general de las matrices dispersas, muy eficiente Choleski basado en la inversión se puede hacer uso de la SparseM paquete. Los pasos básicos para crear un vector distinto de cero de los elementos y de los vectores de números de fila y columna a la que estos elementos se corresponden. SparseM se lo lleva de allí. La he usado en el R rms paquete orm función que permite el manejo de miles de intersecciones en un modelo de regresión ordinal. El método para obtener el equivalente a $(X'X)^{-1}X'Y$ es increíblemente rápido. Independiente inversa no es calculada.

Answer 3

Debido a que esta es una matriz de covarianza, $e=h$, $f=l$, $g=m$, y $A,B,C,D$ son todas positivas. El cambio de notación para manejar una generalización de este problema, considere la matriz simétrica

$$\mathbb{A} = \pmatrix{ \sigma_1^2 & 0 & 0 & \cdots & x_1 \sigma_1 \sigma_n \\ 0 & \sigma_2^2 & 0 & \cdots & x_2 \sigma_2 \sigma_n \\ 0 & 0 & \sigma_3^2 & \cdots & x_3 \sigma_3 \sigma_n \\ \vdots & \vdots & \vdots & \ddots& \vdots \\ x_1 \sigma_1 \sigma_n & x_2 \sigma_2 \sigma_n & x_3 \sigma_3\sigma_n &\cdots &\sigma_n^2}.$$

By comparing this to the question it is apparent that $n=4$, $\sigma_1 = \sqrt{A}$, $\sigma_2 = \sqrt{B}$, $\sigma_3 = \sqrt{C}$, $\sigma_4=\sigma_n=\sqrt{D}$, and $x_1 = e/(\sigma_1\sigma_n)$ etc. There are $n-1$ such $x_i$.

It is evident that

$$\mathbb{A} = \mathbb{U}\mathbb{X}\mathbb{U}$$

where $\mathbb{U}$ is the diagonal matrix with entries $(\sigma_1, \sigma_2, \ldots, \sigma_n)$ and (therefore)

$$\mathbb{X} = \pmatrix{ 1 & 0 & 0 & \cdots & x_1 \\ 0 & 1 & 0 & \cdots & x_2 \\ 0 & 0 & 1 & \cdots & x_3 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ x_1 & x_2 & x_3 & \cdots & 1}.$$

It is straightforward to demonstrate (inductively) that $\mathbb{X}$ is invertible if and only if $x_1^2 + x_2^2 + \cdots + x_{n-1}^2 \ne 1$. In this case set

$$\Delta = \det(\mathbb{X}) = 1 - (x_1^2 + x_2^2 + \cdots + x_{n-1}^2)$$

and observe (using, for instance, the block-matrix formulae in Glen_b's answer) that

$$\mathbb{X}^{-1} = \frac{1}{\Delta}\pmatrix{ \Delta + x_1^2 & x_1, x_2 & x_1 x_3 & \cdots & -x_1 \\ x_1, x_2 & \Delta + x_2^2 & x_2 x_3 & \cdots & -x_2 \\ x_3 x_1 & x_3 x_2 & \Delta + x_3^2 & \cdots & -x_3 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -x_1 & -x_2 & -x_3 & \cdots & 1}.$$

The computation of $\mathbb{U}^{-1}$ is simple--invert each of its diagonal entries--and the subsequent construction of

$$\mathbb{A}^{-1} = \mathbb{U}^{-1} \mathbb{X}^{-1} \mathbb{U}^{-1}$$

is equally easy.

---

To illustrate this algorithm, here is an R implementation. It creates a random symmetric matrix with positive diagonal entries and arbitrary values in the last row and column, inverts it with the algorithm, multiplies that by the original, and compares the result to the identity matrix. The arguments to the inversion function are the vector of $n$ diagonal elements a and the vector of $n-1$ elements in the bottom row of $\mathbb{X}$. Thus $\mathbb{X}$ nunca realmente se necesita para ser construido.

El tiempo es un poco más lento que usar el sparseMatrix de la clase en la Matrix paquete, porque sólo los nativos R operaciones (por ejemplo, outer y t) son empleados y estos no están optimizados para la velocidad. Debido a la inversa de la matriz resultante es nada pero dispersas, parece ser que no hay inherente ventaja del uso de la matriz dispersa cálculos.

invert <- function(a, y) {
  n <- length(a)    # The diagonal elements
  sigma <- sqrt(a)
  x <- (y / (sigma[-n] * sigma[n]))
  Delta <- 1 - sum(x*x)
  X.inv <- outer(c(x, -1), c(x, -1)) + diag(c(rep(Delta, n-1), 0))
  X.inv <- t(t(X.inv) / sigma) / sigma
  return (X.inv / Delta)
}
#
# Create a matrix.
#
n <- 4000
a <- rexp(n)
y <- rnorm(n-1)

X <- diag(a)
n <- length(a)
X[-n, n] <- X[n, -n] <- y
#
# Invert it with the algorithm.
#
system.time(X.inv <- invert(a, y))
#
# Check that it is the inverse.
#
if (n <= 1000) {
  One <- zapsmall(X %*% X.inv)
  if(!all.equal(One, diag(rep(1,n))))
    warning("Not an inverse!") else 
      message("Inverse is correct.")
}

El uso de la dispersa-matriz de problemas, acceso a la solve función en el Matrix biblioteca. Aquí es una ilustración que sigue el código anterior.

library(Matrix)
i <- c(1:n, 1:(n-1), rep(n, n-1))
j <- c(1:n, rep(n, n-1), 1:(n-1))
X.0 <- sparseMatrix(i, j, x=c(a, y, y))
system.time(X.0.inv <- solve(X.0))