Efficiently compute the trace of a sparse matrix times the inverse of a sparse matrix? 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? 
 A: 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:


*

*First compute $\alpha=\|B\|$ approximately using Lanczos

*Now consider $B=\alpha I - C$, so that $B^{-1}=(I-\alpha^{-1}C)^{-1}/\alpha$

*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*}

*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*}

*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)$.

*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).
