Cálculo de la traza del producto de dos matrices

Question

Tengo que calcular $\mbox{trace}(A^{-1}B)$ donde $A$ es una matriz simétrica positiva definida y $B$ es una matriz simétrica, muy dispersa con sólo dos elementos no nulos. Quiero encontrar una manera de calcular la expresión anterior de manera eficiente, especialmente cuando A y B son de alta dimensión como $10000\times 10000.$ ¿Cuál es la mejor manera de hacerlo?

Tengo un puñado de B, cada una muy dispersa, con sólo dos valores no nulos. No puedo almacenar $A^{-1}$ ya que es denso y no tendré suficiente memoria. Alguna forma o truco eficiente para hacerlo de forma eficiente, como propiedades de rastreo o algo así?

Answer 1

En primer lugar, permítanme ofrecer una derivación puramente matemática, y luego intentaremos abordar el problema del almacenamiento una vez que obtenga respuestas a las preguntas que planteé en los comentarios anteriores. Editaré esta respuesta según sea necesario.

Desde $A$ es simétrica y definida positiva, admite una factorización Cholesky $A=LL^T$ , donde $L$ es triangular inferior; y $A=L^{-T}L^{-1}$ . Definamos $M=L^{-1}$ que también es una matriz triangular inferior, por lo que $A^{-1}=M^TM$ ; y dejar que $m_k$ denotan el $k$ columna de $M$ .

Además, usted dice que $B$ es simétrica con dos elementos no nulos. Esto significa que $B$ puede adoptar una de estas dos formas: $$B=\alpha(e_ie_j^T+e_je_i^T) \quad \text{or}\quad B=\alpha e_ie_i^T + \beta e_je_j^T$$ donde $e_k$ es un vector con un uno en el $k$ y ceros en el resto. Consideremos por un momento la primera forma: $$\begin{aligned} \mathop{\textrm{Tr}}(A^{-1}B)&=\alpha\mathop{\textrm{Tr}}(A^{-1}(e_ie_j^T+e_je_i^T))\\ &=\alpha\mathop{\textrm{Tr}}(A^{-1}e_ie_j^T)+\alpha\mathop{\textrm{Tr}}(A^{-1}e_je_i^T) \\ &=\alpha e_j^TA^{-1}e_i+\alpha e_i^TA^{-1}e_j = 2\alpha\left[A^{-1}\right]_{ij} \\ &= 2\alpha e_j^TM^TMe_i = 2\alpha \langle m_i,m_j \rangle \end{aligned}$$ Así que, como puede ver, el trazado requiere exactamente un elemento de $A^{-1}$ , o el producto interior de dos columnas de $M$ . Una derivación similar para el segundo caso da como resultado $$\mathop{\textrm{Tr}}(A^{-1}B)=\alpha\left[A^{-1}\right]_{ii}+\beta\left[A^{-1}\right]_{jj}+\alpha\langle m_i,m_i\rangle+\beta\langle m_j,m_j\rangle$$

Así que espero que ahora quede claro por qué he preguntado: ¿cuántos $B$ ¿hay matrices? ¿Cómo es $A$ ¿Guardado? ¿Qué tipo de operaciones podemos realizar con $A$ ? Estas preguntas son esenciales para determinar qué hacer en este caso. Por ejemplo, si sólo hay un puñado de índices únicos $i,j$ anterior, entonces un enfoque es calcular cada $f_i\triangleq A^{-1}e_i$ utilizando algún tipo de método iterativo, y luego utilizar $e_j^TA^{-1}e_I=e_j^Tf_i$ .

Pero si la mayoría de los índices $i=1,2,\dots,10000$ se representan, puede ser más conveniente hacer algún tipo de factorización Cholesky en la matriz. Sí, puede que no tengas suficiente memoria para hacer una en el núcleo factorización. Pero las factorizaciones Cholesky se pueden hacer fuera del núcleo . Esto implica realizar los cálculos en bloques, leer en memoria sólo los datos suficientes para resolver ese bloque concreto y escribir cada bloque en el disco antes de proceder al siguiente.

Answer 2

Si realmente sólo hay 2 valores no nulos en $B$ entonces se puede calcular $tr(A^{-1}B)$ de $A$ y 2 de su menores . Una matriz de 2 elementos es la suma de 2 matrices de 1 elemento y una matriz de 1 elemento es el producto exterior de vectores de 1 elemento, utilizando bra-ket notación: $$B = X + Y = \left| r_1 \right> x \left< c_1 \right| + \left| r_2 \right> y \left< c_2 \right|$$ Dado que la traza es un operador lineal $$tr(A^{-1}B) = tr(A^{-1}(X+Y)) = tr(A^{-1}X) + tr(A^{-1}Y)$$ Dejemos que $C$ sea el matriz de cofactores de $A$ $$tr(A^{-1}X) = tr( A^{-1} \left| r_1 \right> x \left< c_1 \right| ) = x \left< c_1 \right| A^{-1} \left| r_1 \right> = x (A^{-1})_{c_1r_1} = x \left(\frac{C^{\top}}{det\,A}\right)_{c_1r_1} = x \frac {C_{r_1c_1}}{det\,A} \\ \therefore tr(A^{-1}B) = \frac {x\,C_{r_1c_1} + y\,C_{r_2c_2}}{det\,A}$$ Esto ahorra tener que invertir $A$ . Sólo el determinante y 2 cofactores específicos de $A$ deben ser calculados, por lo que $tr(A^{-1}B)$ puede calcularse con un pequeño factor constante del coste de $det\,A$ .

En los últimos años se ha avanzado en los algoritmos prácticos para los determinantes de grandes matrices dispersas. Esta no es mi área de experiencia, pero aquí hay algunas referencias:

• Erlend Aune, Daniel P. Simpson: Estimación de parámetros en distribuciones gaussianas de alta dimensión En particular, la sección 2.1 ( arxiv:1105.5256 ) (versión más larga publicada versión )
• Ilse C.F. Ipsen, Dean J. Lee: Aproximaciones de determinantes ( arxiv:1105.0437 )
• Arnold Reusken: Aproximación del determinante de grandes matrices simétricas positivas definidas dispersas ( arxiv:hep-lat/0008007 )
• notas para una implementación en la biblioteca shogun

Estos métodos parecen ser principalmente métodos de aproximación que pueden calcular el determinante con una precisión arbitraria a costa de aumentar el tiempo de ejecución, así que puedes elegir el equilibrio entre velocidad y precisión. También parecen evitar la materialización de grandes matrices densas en los cálculos intermedios