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!

• What are the dimensions of each? It looks like $X$ is a vector, but what are $A,B,C$? – Ross Millikan Jul 8 '17 at 0:43
• Good point; I added an edit. So $X$ is also a matrix, everything is a matrix. vec(X) vectorizes X. The dimensions of everything is medium large; basically, it's easy to do A \ B, but not (kron(B.T,A) \ C). – Y. S. Jul 8 '17 at 0:46
• In some cases coarse graining can work well. If you average over, say, 100 by 100 blocks, you get a coarse grained 10 by 10 problem. You can then upscale the solution by merely copying elements to a, say 20 by 20 initial solution. If you then coarse grain the given matrices to 20 by 20 size, then you can use the differences with the upscaled 10 by 10 to 20 by 20 and the real 20 by 20 to get to an efficient way of obtaining the 20 by 20 solution. You can then repeat this doubling (or you can increase by a smaller factor) until you arrive at the original system. – Count Iblis Jul 8 '17 at 1:35
• @whyyes Sorry, I forgot to include additional conditions. We need the existence of solutions to $AY=C$ and $ZB=C$ for $X$ to work. – Sungjin Kim Jul 8 '17 at 3:19
• @whyyes This might be helpful: math.stackexchange.com/questions/585419/… – Sungjin Kim Jul 8 '17 at 7:33

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$

• Cool, this is very nice! I'm trying to draw parallels with some previous comments... the "bad" cases (the ones that don't work for the above answers) would correspond to some 0s in the diagonals of $B$ or $A$, which would complicate things? Anyway, I really like the parallel with Sylvester's equation! – Y. S. Jul 8 '17 at 20:11
• Well at some point when following through with this procedure you might run into the issue of division by zero. E.g. if $A'_{nn}=0$ but $C'_{n1}\neq 0$. That being said in this case, since the diagonal entries of the Schur decomposition are the eigenvalues, $A$ is singular and the system is possibly unsolveable to start with. As a further remark, you might be able to find implementations of these algorithms online. A quick google search gave me this dm.unibo.it/~simoncin/matrixeq.pdf which might be worth a read. – Hyperplane Jul 8 '17 at 23:35

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:

1. if $A$ has full row rank (and therefore a right-inverse), or $C = AU$ for some $U$, and
2. 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?