Calculating the trace of the product of two matrices

I have to calculated $\mbox{trace}(A^{-1}B)$ where $A$ is a symmetric positive definite matrix and $B$ is a symmetric matrix, very sparse with only two elements non zero. I want to find a way that I could calculate the above expression efficiently specially when A and B are high dimensional like $10000\times 10000.$ What is the best way to do it.

I have a bunch of Bs each very sparse with only two non zero values. I cannot store $A^{-1}$ since it is dense and I won't have enough memory. Any efficient ways/tricks to do it efficiently, like trace properties or something?

• Start by decomposing $A = L L^\top$ using a Cholesky decomposition with $L$ a lower diagonal matrix, and $B = \sum e_i e_j^\top$ with $e_i$ a unit vector with a 1 in the i-th element and zero elsewhere. – ja72 Feb 13 '14 at 20:06
• When you say "a bunch of B's", how many are we talking about? It's going to be difficult to come up with as solution here that doesn't require $O(n^2)$ storage; certainly, neither the solution offered below or in the comments avoids this. – Michael Grant Feb 18 '14 at 13:46
• How is $A$ currently stored? In double-precision it requires 80GiB of memory (40GiB if you use a packed format that leverages symmetry). Do you actually have the fully formed version of $A$ stored somewhere, or is it determined algorithmically? – Michael Grant Feb 18 '14 at 13:48
• I have computed and stored inverses of dense 10K x 10K matrices using MATLAB or C++/LAPACK, and the storage is 8 * 10^4 * 10^4 = 800MB (not 80GB), so if you have a decent computer the storage should not be a problem. If you take advantage of the symmetry and use float instead of double the storage is only 200MB to store the inverses. Using @MichaelC.Grant comment below it seems you should have no problems. – Peder Feb 20 '14 at 3:42
• Ha ha! My bad, I think I did 100K x 100K in my calculations. – Michael Grant Feb 20 '14 at 4:17

2 Answers

First, let me offer some purely mathematical derivation, then we will attempt to address the storage problem once I get answers to the questions I posed in the comments above. I will edit this answer as needed.

Since $A$ is symmetric and positive definite, it admits a Cholesky factorization $A=LL^T$, where $L$ is lower triangular; and $A=L^{-T}L^{-1}$. Let us define $M=L^{-1}$, which is also a lower triangular matrix, so $A^{-1}=M^TM$; and let $m_k$ denote the $k$th column of $M$.

Furthermore, you say that $B$ is symmetric with two non-zero elements. This means that $B$ can take one of two forms: $$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$$ where $e_k$ is a vector with a one in the $k$th position and zeros elsewhere. Let's consider the first form for a moment: \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} So as you can see, the trace requires exactly one element of $A^{-1}$, or the inner product of two columns of $M$. A similar derivation for the second case yields $$\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$$

So hopefully now it is clear why I asked: how many $B$ matrices are there? How is $A$ stored? What kinds of operations can we perform with $A$? Those questions are essential for determining what to do in this case. For instance, if there are only a handful of unique indices $i,j$ above, then one approach is to compute each $f_i\triangleq A^{-1}e_i$ using some sort of iterative method, then use $e_j^TA^{-1}e_I=e_j^Tf_i$.

But if most of the indices $i=1,2,\dots,10000$ are represented, it may be more expedient to do some sort of Cholesky factorization on the matrix. Yes, you may not have enough memory---to do an in-core factorization. But Cholesky factorizations can be done out-of-core. This involves performing the calculations in blocks, reading in only enough data into memory to solve that particular block, and writing each block to disk before proceeding with the next.

• Nice approach. Wouldn't storing $f_i$ for all $i$ be as costly as storing $A^{-1}$ itself? – Peder Feb 20 '14 at 3:48
• Indeed. At one extreme, you compute all of the $f_i$ vectors at once; this minimizes computation but at the expense of memory. At the other extreme, you keep only 2 at a time in memory, but you have to do a lot of recomputation. Somewhere in the middle might be the most practical choice! – Michael Grant Feb 20 '14 at 4:01
• That makes a lot of sense. – Peder Feb 20 '14 at 15:55
• @Daniel Mahler Your method of calculating two cofactors and the determinant could be expensive too. I mean calculation of cofactors includes finding determinant too. So you have to calculate three determinants. That is expensive. – user34790 Feb 21 '14 at 3:47
• You probably meant for this comment to go with the other answer; hopefully he will see it anyway... – Michael Grant Feb 21 '14 at 4:29

If there really are only 2 non zero values in $$B$$ then you can compute $$tr(A^{-1}B)$$ from $$A$$'s determinant and 2 of its minors. A 2-element matrix is a sum of 2 1-element matrices and a 1-element matrix is the outer product of 1-element vectors, using bra-ket notation: $$B = X + Y = \left| r_1 \right> x \left< c_1 \right| + \left| r_2 \right> y \left< c_2 \right|$$ Since trace is a linear operator $$tr(A^{-1}B) = tr(A^{-1}(X+Y)) = tr(A^{-1}X) + tr(A^{-1}Y)$$ Let $$C$$ be the matrix of cofactors of $$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}$$ This saves having to invert $$A$$. Only the determinant and 2 specific cofactors of $$A$$ need to be computed, so $$tr(A^{-1}B)$$ can be computed within a small constant factor of the cost of $$det\,A$$.

There has been progress in recent years on practical algorithms for determinants of large sparse matrices. This is not my area of expertise but here are some references:

• Erlend Aune, Daniel P. Simpson: Parameter estimation in high dimensional Gaussian distributions, particularly section 2.1 (arxiv:1105.5256) (longer version published version)
• Ilse C.F. Ipsen, Dean J. Lee: Determinant Approximations (arxiv:1105.0437)
• Arnold Reusken: Approximation of the determinant of large sparse symmetric positive definite matrices (arxiv:hep-lat/0008007)
• notes for an implementation in the shogun library

These methods seem to be primarily approximation methods that can compute the determinant to arbitrary accuracy at the cost of increased running time, so you can choose the balance between speed & accuracy. The also seem to avoid materializing large dense matrices in intermediate calculations

• How would you propose to compute the determinant without $O(n^2)$ storage? – Michael Grant Feb 18 '14 at 13:47
• I have attached some references to the answer. – Daniel Mahler Feb 19 '14 at 5:52