# Solving for the trace of an unknown matrix in a linear matrix equation

I would like an efficient way to solve for $\mathrm{tr}(X)$, where $$XA=B$$ and I know the matrices $A$ and $B$, but not $X$. $A$ and $B$ are Hermitian and sparse and I would like to avoid computing $X$ in its entirety, only to take its trace and nothing else.

There is a detailed answer to a similar question involving the Lyapunov equation, but I hoped that in this simpler case that there may be less involved or more efficient options.

I have been impressed recently by the performance of iterative techniques, such as those of Yousef Saad in solving matrix-vector equations and I wonder whether there is a similar technique that can be applied here? Sadly, the best I can think of at the moment, however, is to solve row-by-row and summing the element corresponding to the diagonal. Of course, that would be counterproductive when an LU decomposition approach would be much more efficient (assuming I have to compute the whole matrix).

Thanks,

Joly

• I am assuming $A$ is square and that you do not want to compute $A^{-1}$ (otherwise, $\operatorname{tr}(X)=\operatorname{tr}(BA^{-1})$). – parsiad Nov 7 '16 at 16:51
• Thanks. Yes, $A$ is square and I would very much like to avoid computing $A^{-1}$, since $A$ is very sparse and its inverse is dense. – holyjoly Nov 7 '16 at 16:56
• A simple iteration would be $X^{n+1}=B+X^n(I-A)$ (you would have to also check to see the iteration is contractive), and you could compute the trace of an iterate that is desirably close to your solution. You could also do a more sophisticated iteraiton. Saad's book is actually a pretty good source to look at if you want to speed things up beyond the above simple iteration. – parsiad Nov 7 '16 at 17:00
• Note also that you can right-multiply the equation by a diagonal matrix $D$ (of appropriate dimension) without changing the set of solutions. – parsiad Nov 7 '16 at 17:11
• My concern with this is that I know that $X$ is approximately dense and I assume that the diagonal elements in this iteration will not converge more rapidly than the others in general? If not, then I will be faced with a dense matrix in memory before the trace converges. – holyjoly Nov 7 '16 at 17:30

For some matrix function, $\mathbf{f}(\mathbf{X})$, the trace may be approximated as: $$\mathrm{tr}(\mathbf{f}(\mathbf{X}))\approx\sum_{i=1}^N \vec{v}^T_i \mathbf{f}(\mathbf{X})\vec{v}_i,$$ where $\vec{v}_i$ are some set of vectors, which may be random or otherwise. For my problem, $\mathbf{f}(\mathbf{X})=\mathbf{X}=\mathbf{B}\mathbf{A}^{-1}$ and so the trace may be approximated as $$\mathrm{tr}(\mathbf{X})\approx\sum_{i=1}^N \vec{u}^T_i\vec{w}_i,$$ where $\vec{u}_i^T=\vec{v}_i^T\mathbf{B}$ and $\vec{w}_i$ may be obtained by solving $\vec{v}_i=\mathbf{A}\vec{w}_i$. In my case, if $\vec{v}_i$ are normally distributed random numbers, N must be greater than is practical. Instead, if I take vectors from a Hadamard matrix as suggested in those papers, then for the problem I want to solve I seem to require at most 512 vectors to converge to $10^{-8}$ in the trace. This is remarkably cheap considering the dimension of the matrices which I am dealing with.
I mocked this up in Matlab:  acc=0; N=2048; H=hadamard(2^ceil(log2(size(A,1)))); exactr=trace(B*inv(A)); for i=1:N; v=H(i,1:size(A,1))'; acc=acc+((v'*B)*(A\v)); if(i-(2^floor(log2(i)))==0) fprintf('%d %20.16f %20.16f %20.16f\n',i,acc/i,exactr,abs((acc/i)-exactr)); end; end; trest=acc/N