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]:

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][2] and his book on *Functions of Matrices* for further information).


  [1]: http://press.princeton.edu/titles/9104.html
  [2]: http://www.maths.manchester.ac.uk/~higham/