Can anyone suggest a way to numerically compute the following matrix vector product?


$$u=A^{-1}b=(AA\otimes BB + AB \otimes BA)^{-1}\operatorname{vec}(C)$$


Here $AA,BB,AB,BA$ and $C$ are $d\times d$ matrices with $d\approx 1000$. A naive method is to expand Kronecker products and use linear solver for an answer in $O(d^6)$ time, but I need to do it in $O(d^3)$ time.


Not sure if helps, but $AA,AB,BA,BB$ are blocks of a covariance matrix, ie for partitioned random variable $(x,y)$
$$
\text{cov}\left(
\begin{array}{ll}
x\\y
\end{array}\right)
=\left(\begin{array}{cc}AA& AB\\BA& BB\end{array}\right)$$


I tried Jacobi-like iterative method, writing $A=M-N$ with $M=AA\otimes BB$ and $N=-AB\otimes BA$ and iterating $u=M^{-1}Nu+M^{-1}b$ which is cheap because of Kronecker products, but this sometimes diverges. Empirically, divergence seems to happen when $x$ and $y$ are highly correlated with each other. Any suggestions/literature pointers are appreciated!



Computing $(AA\otimes BB + AB \otimes BA)^{-1}$
Your equation is equivalent to
$$
(B^{-1}A)X + X(B^{-1}A)^T = B^{-1}A^{-1}CB^{-T}B^{-T},
$$
where $X$ is the matricization of $u$. This is a [Sylvester equation](https://en.wikipedia.org/wiki/Sylvester_equation), which can be solved for $X$ in $O(d^3)$ time.


You'll get some numerical ugliness if $A$ or $B$ are close to non-invertible, but that might be unavoidable in some sense: if either of $A$ or $B$ are ill-conditioned then so is the original coefficient matrix $AA \otimes BB + AB \otimes BA$.


