Return to Revisions

Here is a simple minded numerical approximation (ignoring roundoff etc) that works for the case that a factorization $A=UU^T$ is known (or $A^{1/2}$ is available):

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$ 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}$ of $A=UU^T$, then some more thought is needed. Yet another level of approximation is added if we build a subroutine to compute $A^{1/2}u$ --- such subroutines are the subject of intense study in numerical linear algebra (see e.g., Nick Higham's webpage).

EDIT. Here is a rough upper-bound, not sure if it suffices. Recall that \begin{equation*} \text{tr}(AB^{-1}) \le \langle \lambda(A), \lambda(B^{-1}) \rangle. \end{equation*} The right hand side above can be bounded using Hölder's inequality to make it easier to compute, for instance by $\lambda_{\max}(A)\text{tr}(B^{-1})$, which in turn can be easily computed using the above ideas. Alternatively, we can upper bound it using $\text{tr}(A)\lambda_{\max}(B^{-1})$, which requires computing the smallest eigenvalue of $B$, again Lanczos will help. This particular bound is much easier to obtain using existing tools in Matlab (see 'eigs').