8
$\begingroup$

How can I efficiently compute $\mathrm{trace}(A(B^{-1}))$ where $A$ and $B$ are both sparse symmetric PSD $n \times n$ matrices, both with $O(n)$ non-zero entries? If it helps, the pattern of non-zero entries in $A$ and $B$ can be the same.

Alternatively, is there a tight upper bound on this quantity that I can compute efficiently, e.g. in $O(n \log(n) )$ time?

Improve this question
$\endgroup$
6
  • $\begingroup$ if you just want to measure how well $B$ approximates $A$, why not use the Frobenius norm ${\rm tr}\,[(A-B)(A-B)^t]$ ? $\endgroup$
    – Carlo Beenakker
    Commented Dec 15, 2015 at 7:09
  • $\begingroup$ @CarloBeenakker Good point. I have some issues with a measure under which $A=I$ and $A=\begin{bmatrix}1 & 10^6\\ 10^6 & 1\end{bmatrix}$ both approximate equally well $B=I$. (One can play around and find examples with the same nonzero pattern, too.) $\endgroup$
    – Federico Poloni
    Commented Dec 15, 2015 at 7:36
  • $\begingroup$ @CarloBeenakker -- I removed the parenthetical about B approximating A, and quantity of interest being a measure of the quality of the approximation. It's better to think of the quantity of interest as a term in an objective function that I'm trying to minimize. $\endgroup$
    – Jeff
    Commented Dec 15, 2015 at 17:12
  • $\begingroup$ one "slightly" better than brute force approach to compute $AB^{-1}$ would be to compute only those cofactors of $B$ that would get multiplied with rows of $A$. There is work on determinants of sparse matrices that can get us started . but the worst case complexity of this method is still at least $n^2$ $\endgroup$
    – Pushpendre
    Commented Dec 15, 2015 at 17:47
  • $\begingroup$ Do you have any additional information on $A$ and $B$? Are they positive definite? Are $A$ and $B$ close in some sense? I think the general case is rather difficult. $\endgroup$
    – user35593
    Commented Dec 15, 2015 at 18:15

1 Answer 1

Reset to default
3
$\begingroup$

Simple bounds

A simple upper bound is \begin{equation*} \text{tr}(AB^{-1}) \le \min\left(\lambda_{\max}(A)\text{tr}(B^{-1}), \text{tr}(A)\lambda_{\max}(B^{-1})\right). \end{equation*} Both these bounds are numerically "easy" to compute using Lanczos. For computing $\text{tr}(B^{-1})$ a randomized trace estimator can be used (following the more general idea outlined below).

Numerical approximation

Here is a simple approach, motivated by this nice book:

  1. First compute $\alpha=\|B\|$ approximately using Lanczos
  2. Now consider $B=\alpha I - C$, so that $B^{-1}=(I-\alpha^{-1}C)^{-1}/\alpha$
  3. After that, consider \begin{equation*} \text{tr}(AB^{-1}) = \frac1\alpha\text{tr}(A^{1/2}(I-\alpha^{-1}C)^{-1}A^{1/2}) \end{equation*}
  4. Now use the von Neumann series \begin{equation*} \text{tr}(A(I-\alpha^{-1}C)^{-1})=\sum_{k\ge0} (-1)^k\alpha^{-k}\text{tr}(A^{1/2}C^kA^{1/2}) \end{equation*}
  5. Let $u \sim \mathcal{N}(0,I)$ be a mean-zero spherical Gaussian rv. Then, we approximate the above quantity by taking $m$ samples, $u_1,\ldots,u_m$ and iteratively computing $u_i^TA^{1/2}C^kA^{1/2}u_i$ for $1\le i \le m$. Observe that the key subroutine that we have is to iteratively compute $z^TC^kz = z^TC(C^{k-1}z)$.
  6. In expectation this will be an estimator for the trace in question since $E[\text{tr}(u^TA^{1/2}C^kA^{1/2}u]=E[\text{tr}(A^{1/2}C^kA^{1/2}uu^T)]$ and $E[uu^T]=I$ by assumption.

If you do not have access to $A^{1/2}$ (or a Cholesky factorization of it) then an additional level of approximation arises by building a subroutine to compute $A^{1/2}u$. Such $f(A)b)$ family of subroutines are the subject of research interest in numerical linear algebra (see e.g., Nick Higham's webpage and his book on Functions of Matrices for further information).

Improve this answer
$\endgroup$
5
  • $\begingroup$ Interesting approach. I don't think I can access to $A^{1/2}$ unfortunately. As I've convinced of it, $A$ isn't low rank and so I think $A^{1/2}$ would be dense, with $O(n^2)$ non-zero entries. $\endgroup$
    – Jeff
    Commented Dec 15, 2015 at 18:54
  • $\begingroup$ @Jeff $A^{1/2}$ does not have to be explicitly available. For instance, all you need is access to a subroutine that approximates the vector $A^{1/2}u$, and this approximation can be also done using Lanczos (using Krylov methods). But we need to do this $m$ times for the $m$ random vectors $u_i$, which can be done reasonably fast --- at no point, there is any dense matrix involved.... $\endgroup$
    – Suvrit
    Commented Dec 15, 2015 at 19:21
  • $\begingroup$ @ Suvrit , very pretty post (+1, but that will not hardly change your total!). What about the complexity ? I think that it is in $O(mn^2)$ with a "not small" constant. How do you choose $m$ ? $\endgroup$
    – loup blanc
    Commented Feb 14, 2016 at 18:26
  • $\begingroup$ @loupblanc thanks! Why the $n^2$? The subroutines use only matrix-vector calls, and under the OP's assumption, each matrix*vector call costs $O(n)$. How large $m$ has to be is harder to determine and requires more analysis (it's the same thing as figuring out how many random samples one should draw to be able to compute something to a given accuracy). But maybe you are right it won't be truly $O(mn)$ because of the nested call to $A^{1/2}u$ --- there are other fancier ideas also possible, if I get a chance I'll try to mention those. $\endgroup$
    – Suvrit
    Commented Feb 14, 2016 at 22:22
  • $\begingroup$ @ Suvrit , yes I forgot "sparsity". $\endgroup$
    – loup blanc
    Commented Feb 17, 2016 at 10:33

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .