¿Hay una manera más rápida para calcular algunos elementos diagonales de la inversa de una matriz definida positiva simétrica enorme?

Question

Pregunté esto en LO primero, pero decidió mover la parte de matemáticas de mi pregunta aquí.

Considere la posibilidad de una $p \times p$ simétrica y positiva definida la matriz de $\bf A$ (p=70000, es decir, $\bf A$ es de aproximadamente 40 GB con 8 bytes dobles). Queremos calcular los tres primeros elementos de la diagonal de la matriz inversa $({\bf A^{-1}})_{11}$, $({\bf A^{-1}})_{22}$ y $({\bf A^{-1}})_{33}$.

He encontrado este artículo por James R. Montón que parece resolver este problema exacto sin calcular el total inverso $\bf A^{-1}$. Si he entendido correctamente, se le calcula primero la descomposición de Cholesky, es decir, la parte superior triangular de la matriz de $\bf R$ que satisface $\bf A=R^T R$, el cual necesita de $\frac16p^2+\frac12p^2-\frac23p$ operaciones de punto flotante (multiplicaciones/divisiones) utilizando el LINPACK función SPOFA. Luego procede a calcular individual de los elementos de la diagonal de la matriz inversa $({\bf A^{-1}})_{ii}$ utilizando una expresión que explota la dispersión de ${\bf R}^T{\bf y}={\bf e}_j$ y que requiere de $\frac12(p-i)^2+\frac52(p-i)+2$ operaciones de punto flotante. (No entiendo los detalles completos de esta, por lo que no puede en la actualidad suma correctamente).

El documento se basa en LINPACK; no es citado por nadie, por lo que parece que nadie se preocupaba por los últimos 23 años? Después de leer esto, me pregunto si esta es la mejor manera de hacer las cosas, o si un moderno LAPACK enfoque basado en evitar la descomposición de Cholesky?

En resumen, hay una forma más rápida para calcular los elementos de la diagonal de la matriz inversa?

Answer 1

Yo sugiero que busque en el libro de las Matrices, los Momentos, y en Cuadratura con Aplicaciones. La idea básica es que usted tiene un $n \times n$ positiva semi-definida la matriz de $A$, y se desea calcular la forma cuadrática $u^T f(A)u$ donde $u \in \mathbb{R}^n$. En el caso de que usted desee calcular $$ A^{-1}_{ii} = e_i^T f(A)e_i, $$ donde $e_i \in \mathbb{R}^n$ es un vector de ceros, con un 1 en el $i$th la posición y la función de $f$ se define como $f(\lambda) = \lambda^{-1}$.

Puede activar la forma cuadrática en una integral de Stieltjes, ya que $$ u^T f(A)u = u^T Q^T f(\Lambda) Qu = \alpha^T f(\Lambda) \alpha = \sum_{i=1}^n f(\lambda_i) \alpha_i^2 = \int_a^b f(\lambda) d\alpha(\lambda), $$ donde el autovalor de la descomposición de $A$ está dado por $A = Q\Lambda Q^T$ $\Lambda$ es una matriz diagonal que contiene los autovalores de a $A$, y el vector $\alpha = Qu$. La integral puede ser estimada a partir de la cuadratura de Gauss a través de $$ I[f] = \int_a^b f(\lambda) d \alpha(\lambda) = \sum_{i=1}^N \omega_i f(t_i) + R[f], $$ donde $\{\omega_i\}_{i=1}^N$ son un conjunto de pesos y $\{t_i\}_{i=1}^N$ son un conjunto de nodos en el que evaluar la función de $f$, e $R[f]$ es un resto o término de error. Los valores de la $\omega$'s y $t$'s son desconocidos y deben ser resueltos. Los valores para los pesos pueden ser calculadas a través de un algoritmo iterativo similar al algoritmo de Lanczos para el cómputo de los autovalores. Los valores de los nodos puede ser obtenida a partir de los componentes de un vector propio de una matriz de derivadas de $A$. Este cálculo se puede hacer de manera eficiente. Para más detalles, ver el libro, así como de esta conferencia sobre el tema por James Lambers.

Los subyacentes de las matemáticas y álgebra lineal, puede parecer un poco de miedo al principio, pero te aseguro que esto conduce a un modo bastante sencillo y eficiente algoritmo. Escribí el código de Matlab para calcular los pesos y los nodos para una clase en la escuela de posgrado. No fue muy difícil. Echa un vistazo a el libro. La buena suerte.

Answer 2

Me sugirió que un parcial de descomposición de Cholesky, debería ser posible aquí, pero en base a un comentario que he recibido, empecé a pensar en mi respuesta de nuevo y se dio cuenta de que mi suposición era incorrecta. Disculpas si puedo llevar a nadie por el camino equivocado.

---

Usted tendrá que utilizar una completa descomposición de Cholesky, pero el problema para deducir la inversa puede ser reducida en escala para guardar redundante de la computación.

Si el SPD matriz $\mathbf{A}$ está escrito en forma de bloque como

$\mathbf{A}=\begin{bmatrix} \mathbf{A}_{00} & \mathbf{A}_{10}^T \\ \mathbf{A}_{10} & \mathbf{A}_{11} \end{bmatrix}$

con $\mathbf{A}_{00}$ $3\times3$ SPD bloque, su inversa será dada por

$\mathbf{A}^{-1}=\begin{bmatrix} \mathbf{B}_{00} & \mathbf{B}_{10}^T \\ \mathbf{B}_{10} & \mathbf{B}_{11} \end{bmatrix}$

El equivalente a la descomposición de Cholesky de a $\mathbf{A}$ es el dado por

$\mathbf{R}=\begin{bmatrix} \mathbf{R}_{00} & \\ \mathbf{R}_{10} & \mathbf{R}_{11} \end{bmatrix}$

La matriz resultante ecuación a resolver por el inverso del bloque de matriz $\mathbf{B}_{00}$ puede ser reducida a

$\begin{bmatrix} \mathbf{R}_{00} & \\ \mathbf{R}_{10} & \mathbf{R}_{11} \end{bmatrix} \begin{bmatrix} \mathbf{R}_{00}^T & \mathbf{R}_{10}^T\\ & \mathbf{R}_{11}^T \end{bmatrix} \begin{bmatrix} \mathbf{B}_{00} \\ \mathbf{B}_{10} \end{bmatrix}=\begin{bmatrix}\mathbf{I}_{00} \\ \mathbf{0} \end{bmatrix}$

y resuelto a través de

$\begin{bmatrix} \mathbf{R}_{00} & \\ \mathbf{R}_{10} & \mathbf{R}_{11} \end{bmatrix} \begin{bmatrix} \mathbf{X}_{0} \\ \mathbf{X}_{1} \end{bmatrix} =\begin{bmatrix}\mathbf{I}_{00} \\ \mathbf{0} \end{bmatrix}$

y

$\begin{bmatrix} \mathbf{R}_{00}^T & \mathbf{R}_{10}^T\\ & \mathbf{R}_{11}^T \end{bmatrix} \begin{bmatrix} \mathbf{B}_{00} \\ \mathbf{B}_{10} \end{bmatrix}=\begin{bmatrix}\mathbf{X}_{0} \\ \mathbf{X}_{1} \end{bmatrix}$

La solución que usted está buscando va a ser la diagonal entradas de $\mathbf{B}_{00}$. El lado derecho de esta expresión sólo es $n\times 3$, y esto, combinado con la baja de la estructura triangular de $\mathbf{R}$ debe ser de alrededor de la manera más eficiente para obtener la solución que usted requiere.

Answer 3

Una alternativa podría ser el método de gradiente conjugado (previo) para la solución de Ax = b con b = (1.0.0). Es la primera entrada de x $(A^{-1})_{11}$

Answer 4

La primera parte de aquí abajo no es una "respuesta", pero un comentario extendido/seguimiento-pregunta a talonmies anteriores de respuesta/comentario

Me he extendido su notación un poco y consiguió una más clara exposición: en Primer lugar vamos a extender la Xde la matriz a ser la inversa de la completa R-matriz, de tal manera que $\small \mathbf A^{-1}= \mathbf B =\mathbf X^\tau \mathbf X $ y también vamos a ampliar los índices de X , entonces
$\qquad \small \mathbf X_{00} =\mathbf R_{00}^{-1} $, $\small \qquad \mathbf X_{11} =\mathbf R_{11}^{-1} \quad $ pero lamentablemente $\small \mathbf X_{10} \ne \mathbf R_{10}^{-1} $ .

Por la inversa de la relación de la $\small \mathbf X =\mathbf R^{-1} $ tenemos, que $\small \mathbf X_{10} \mathbf R_{00}+ \mathbf X_{11} \mathbf R_{10} =\mathbf 0 $ e lo $\small \mathbf X_{10} = - \mathbf X_{11} \mathbf R_{10} \mathbf R_{00}^{-1} $ (a partir de la cual queremos recibir la parte superior izquierda de la diagonal de elementos de B por a $\small \mathbf X^\tau_{00} \mathbf X_{00} + \mathbf X^\tau_{10} \mathbf X_{10} = \mathbf B_{00} $ )

Aquí $\small \mathbf X_{11} $ representa la parte del problema que queremos para reducir el esfuerzo computacional, debido a que esta es la inversión de la (gran) parte restante de la cholesky-factor.Pensé en introducir el concepto de la pseudoinverse pinv(X) por el que podríamos continuar:

$\pequeño \begin{array} {rcl} \mathbf Y &=& - ( \mathbf X_{11} \mathbf R_{10} \mathbf R_{00}^{-1} )^{-1} \\ &=& - \mathbf R_{00} \operatorname{pinv} (\mathbf R_{10}) \mathbf R_{11} \\ \mathbf X_{10} &=& \operatorname{pinv} ( \mathbf Y) \end{array}$

Esto permitiría invertir pequeñas matrices sólo, pero a pesar de la pseudoinverse-operaciones funciona bien en muchos casos, no podría hacer el trabajo correctamente para el pinv(R10)-parte; para todas las variantes de mis cálculos siempre necesitaba la versión completa R11 en esta o aquella versión. ¿Ves algún general de la razón, por qué el pinv(R10) no está trabajando lo suficiente aquí?

---

Práctica de cálculo

---

Hmm, por lo que aceptar que necesitamos todo el cholesky-descomposición de todos modos, el proceso podría incluso ser descrito más corto:

a) considerar la simetría del SPD-matriz de Una , buscamos la diagonal-elementos de la 3x3-submatriz $\small \mathbf B_{00} $ $ \small \mathbf B = \mathbf A^{-1} $

b) realizar el cholesky-descomposición de abajo a arriba en la matriz R; denotar $\small \mathbf R_{00} $ la parte superior izquierda de 3x3 triangular superior submatrix (pensando en términos de una correlación de la matriz de Un esto representa el parcial o "inexplicable" de varianza/covarianza)

c) invertir $\small \mathbf R_{00} $ conseguir $\small \mathbf X_{00} $ (esto también es barato porque ya triangular)

d) en $ \small \mathbf B_{00} = \mathbf X_{00}^\tau \mathbf X_{00} $ nos encontramos con los necesarios elementos de la diagonal. (Esta última operación puede incluso ser sustituido por una simple suma de los cuadrados a lo largo de las columnas de X )

---

Aquí está un ejemplo completo utilizable para Pari/GP. A tenerlo tan claro como sea posible sin optimizaciones, errorchecks etc:

Cultivo(M,dim) = matriz(dim dim,r,c,M[r,c]) \\ reduce el tamaño de una matriz

 \\ procedimiento para mostrar la primera dim diagonal-entradas de A^-1
{invdiag(A,dim=3)=local(rs=filas(Una),lV,R,X);

 \\ reducir Una por una cholesky-proceso de dim x dim-residual de la matriz A_00
 \\ la cholesky proceso que va de abajo hacia arriba
 forstep( d = rs, dim+1, -1, 
 lV = [, d]/sqrt(a[d,d]);
 A = de Cultivo( Una - lV*lV~ , d-1 );
);
 \\después de esta es el dim x dim residual de la matriz de 

 \\ calcular el cholesky-factor de que Un residual en la matriz R
 \\ acabamos de continuar trabajando de abajo hacia arriba, pero no es más recortada ahora
 R = matriz(dim dim); 
 forstep( d = dim, 1, -1, 
 R[,d] = [d]/sqrt(a[d,d]);
 A = a - R[,d] * R[,d]~;
);

 X = R^-1; \\ inversa de R de tamaño dim dim x
 \\ y para extraer los valores de diag (^-1) =diag( R^-1~ * R^-1 )
 \\ basta la suma de los cuadrados a lo largo de las columnas
 lv = vector(dim,c,sum(r=1,c, X[r,c]^2) );
retorno(lv);
}

\\ crear una matriz de la muestra Un
A = matriz(8,8,r,c,binomial(r-1,c-1)/2^(c-1))
A = a* a~ \\ hacer un simétrica positiva definida una

de impresión(invdiag(A,3)) \\ muestran los 3 primeros resultados parciales de cholesky-método
de impresión(diag (^-1)) \\ mostrar la verdadera resultados por inversión de la matriz completa

\\ de salida:
 [21845.0000000, 980612.000000, 8259152.00000]
 [21845, 980612, 8259152, 21815360, 21017856, 7373824, 806912, 16384]~

Answer 5

He aquí una página y el paquete de la solución de este problema exacto. Se hace uso de la descomposición de Cholesky, pero tenga en cuenta que la descomposición de Cholesky de una gran escasa matriz se puede calcular rápidamente por un paquete como CHOLMOD.

Los documentos que cita son:

• Takahashi K, Fagan J, y Chen M-S (1973). La formación de una escasa autobús de la impedancia de la matriz y su aplicación al estudio de corto circuito. En El Poder De La Industria De Equipo De Aplicación De Los Procedimientos De La Conferencia. IEEE Power Engineering Society.

• Vanhatalo J y Vehtari Un (2008). Modelado local y global de los fenómenos con escasa Gaussiano procesos. En los Procedimientos de la 24ª Conferencia de la Incertidumbre en Inteligencia Artificial. AU IA Prensa.