Suppose that we have two operators $A$ and $B$ which satisfy $[A,B]\neq0, A^{T}=A, B^{T}=B$. I'm going to keep these operators vague to be concise, but I have precise definitions of $A$ and $B$ in mind.
In addition, the operator $A$ is diagonal. Diagonalizing $B$ is not trivial but it can be done with some effort, so let us assume that we have the eigenvalues and eigenvectors of both $A$ and $B$.
What I would like to do is compute $\text{Tr}[(AB)^n]$, where $n$ is a positive integer. In general this is computationally expensive unless we can diagonalize the product $AB$. Unfortunately I don't know how to do that, so my question is if it is possible to use our knowledge of the spectrum of $A$ and $B$ to get a computational speedup? Again I think the answer is usually "no", but since $A$ is diagonal maybe this case is manageable. Any suggestions for how to proceed are welcome, these matrices are big and any speedup helps.
EDIT: It's been suggested that I (1) explain what $A$ and $B$ are and (2) explain what I mean by computationally expensive, so here we go.
Let $t$ be a positive integer and consider a spin-1/2 chain of length $t$. The states of the chain are in a Hilbert space $\mathcal{H}_t$ of dimension $2^t$. We will want to "double" the Hilbert space by defining a second spin chain, also of length $t$, and we denote the doubled Hilbert space as $\mathcal{H}_{2t}=\mathcal{H}_{t}\otimes\mathcal{H}_t$.
Next we define the following two Ising-like operators which act on $\mathcal{H}_t$: $$U_X = \exp\left[-i(-\frac{\pi}{4}+ig)\sum_{\tau=1}^{t}X_{\tau}\right]$$ $$U_Z = \exp\left[-i(-\frac{\pi}{4}+iJ)\sum_{\tau=1}^{t}Z_{\tau}Z_{\tau+1}\right]$$
where $g,J$ are real numbers that we can tune and $X, Z$ are Pauli matrices, and periodic boundaries are assumed. Note that, in general these matrices are not unitary, but they are symmetric.
We are almost ready to define $A$ and $B$, but first we need to define a projector $P$ that acts on $\mathcal{H}_{2t}$. Let $|\vec{n}_1,\vec{n}_2\rangle\in\mathcal{H}_{2t}$ and denote the total $S_z$ of $\vec{n}_i$ by $N_i$. Then
$$P = \sum_{\vec{n}_1,\vec{n}_2:N_1=N_2}|\vec{n}_1,\vec{n}_2\rangle\langle\vec{n}_1,\vec{n}_2|$$
is a Hermitian projector which projects onto the subspace $\mathcal{H}_P\subset\mathcal{H}_{2t}$ where both spin chains have the same total $S_z$. With this projector, we finally have
$$A = P(U_{Z}\otimes U_{Z}^{\ast})P, \quad B = P(U_{X}\otimes U_{X}^{\ast})P$$
It's easy to see that the operator $A$ is diagonal in the computational basis of $\mathcal{H}_P$. The operator $B$ is much more complicated since we cannot form eigenstates of the Pauli X operator in $\mathcal{H}_P$. However, I have a quasi-analytic way of diagonalizing $B$ that reduces it to a numerically tractable problem.
Now, let us return to the issue of computational complexity. In my original post I was looking to compute the trace of $(AB)^n$ in an efficient manner. To understand my concern, note that the dimension of the projected Hilbert space is big: one can check that $\dim(\mathcal{H}_P) = \binom{2t}{t}\sim 4^t/\sqrt{t}$. In contrast, my method of diagonalizing $B$ has complexity $O(\text{poly}(t))$ (I think it's $O(t^3)$ but there may be a hidden log term or something). I was hoping to find a comparable reduction for computing the desired trace, but that might be asking too much.
Another option would be to simply approximate the trace, but this seems tricky. I can prove that $\text{Tr}[(AB)^n]\in\mathbb{R}_{+}$ for all $n$, but the eigenvalues of $AB$ are not generically real or positive. So it is not clear how we could truncate the spectrum in any way to get a reasonable approximation.
