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

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

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.

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