# How to solve for the matrix $X$ in the following equation $AXB + X = CD$

How to solve for the matrix $X$ in the following equation $AXB + X = CD$?

$A$ and $B$ are full rank symmetric matrices, and there is no structure to $CD$. $CD$ just could be $C$.

• Any info on what properties A, B, C, and D have? Why is it CD and not just C?
– sheß
Aug 31, 2015 at 22:02
• I have updated the question. Aug 31, 2015 at 22:08
• It's a system of linear equations--solve it as you would any system. Aug 31, 2015 at 22:38
• Is $X$ necessarily square (such as assumed in Mark L. Stone's answer)? Note that if $X$ is $r\times c$ the conditions in your equation all still seem to work fine; you'd usually just get either an over-determined or under-determined system. In the first case there's generally no exact solution (but you may be interested in something like, say, a least-squares approximation in that case) and in the second case you should get a set of solutions on a linear subspace. Sep 1, 2015 at 0:57
• @Glen_b - Even if $r \ne c$, the system may still be well-posed. Whether the system is well-posed or not depends entirely on the spectra of $A$ and $B$. Sep 1, 2015 at 1:31

Without loss of generality $$CD = C$$. By properties of the Kronecker product, the problem is equivalent to $$(B^T \otimes A) vec(X) + vec(X) = vec(C)$$ with solution $$vec(X) = \left( B^T \otimes A + I \right)^{-1} vec(C)$$ assuming the inverse exists. Here $$vec(A)$$ is the vector obtained from the matrix $$A$$ by stacking its columns.

Let $$A$$ and $$B$$ have eigenvalues $$\mu_i, \lambda_j$$. Then the inverse exists iff $$\mu_i \lambda_j \ne -1$$ for all $$i,j$$. In particular $$A$$ or $$B$$ do not need to have full rank.

• While this is a "correct" solution, explicitly forming the matrix inverse is not a numerically good way of solving the problem. Sep 1, 2015 at 2:56
• It should be $(B^{\rm T}\otimes A)vec(X)+vec(X)=vec(C)$, right?minor error! Aug 17, 2016 at 6:15
• @user252783 - Thank you, I fixed it. Aug 2, 2020 at 22:32

As whuber wrote in a comment "It's a system of linear equations--solve it as you would any system."

Here's how you could do it in Octave (or MATLAB) with YALMIP (free) plus the free solver GLPK.

>> n=5;A=randn(n);B=rand(n);C=randn(n);D=rand(n); % generate random data
>> X=sdpvar(n,n,'full') % declare X to be an n by n matrix variable
>> optimize(A*X*B+X==C*D,[],sdpsettings('solver','glpk')) % find a solution to the constraint
>> value(X) % here is the solution
-0.6236   -2.2800   -1.7939   -0.5188   -1.0156
0.1853    1.4420    0.9496    0.4852   -0.3270
-0.7491   -2.6045   -1.8847   -1.2238   -0.8824
-0.8494    1.9423   -0.1014   -1.7763    0.1937
0.4692    0.9236    0.4308    0.9344   -0.0396

• Following up on Glen_b's comment above: If the equation has a unique solution, the above technique will find it. If there's an infinity of solutions, this method will find one of them. By inserting an objective function in lieu of [] in optimize, a solution which optimizes the objective function given satisfying the equation will be found. If the problem is overdetermined and no solution exists, that will be reported (a least squares solution could be found if desired.) I gave an example using all n by n matrices, but in fact this technique is general to any compatiblly dimensioned matrices. Sep 1, 2015 at 2:52