How to solve $AXB = C$ efficiently (numerically) This might be a simple question but I can't think of an obvious answer. How do I solve for $X$ in $AXB=C$ efficiently? I don't want to use $(B^T\otimes A) \operatorname{vec}(X) = C$ because $B^T\otimes A$ is very large. Is there a better 1-step solution? Or should I just use gradient descent? 
$A$, $B$, $C$, and $X$ are all matrices, with dimensions in the order of 1000 x 1000
Thanks!
 A: To simplify, you should do a Schur decomposition of $A$ and $B$, i.e. $A=UA'U^*$ and $B=VB'V^*$ with $U,V$ unitary and $A',B'$ upper triangular. Then the system becomes
$$ AXB=C \iff  U A' U^* X V B' V^* = C \iff  A' Y B' =C'$$
with $Y =U^*XV$ and $C' = U^*CV$. The resulting system is efficiently solvable by carefully carrying out the matrix multiplications. Observe that, since $A'$ and $B'$ are upper triangular,
$$ A'YB' = 
\begin{bmatrix}*&*&*\\*&*&* \\ A'_{nn}Y_{n1} & \cdots & A'_{nn}Y_{nn}&\end{bmatrix}\cdot 
\begin{bmatrix}B'_{11}&\cdots& B'_{1n}\\ &\ddots&\vdots\\&&B'_{nn}\end{bmatrix} $$
Thus $(A'YB')_{n1} = A'_{nn}Y_{n1}B'_{11} = C'_{n1}$ determines $Y_{n1}$. But then we can continue with 
$$(A'YB')_{n2} = A'_{nn}Y_{n1}B'_{12} + A'_{nn}Y_{n2}B'_{22} = C'_{n2}$$
which subsequently determines $Y_{n2}$ and so on and so forth.
This idea is basically a variation on the Bartels-Stewart Algorithm which solves Sylvester's equation $AX - XB = C$
A: To summarize the comments, there are a number of ways of approximating a solution efficiently (using some optimization method, or using a coarse graining / upscaling procedure). 
For an exact solution, we can use $X = A^\dagger C B^\dagger$ (or, in matlab, $X = (A \backslash C) / B$) if:


*

*if $A$ has full row rank (and therefore a right-inverse), or $C = AU$ for some $U$, and

*if $B$ has full column rank (and therefore a left-inverse), or $C = VB$ for some $V$.


Then $AXB = C$. 
For exact solutions, this seems the best we can do?
