# Solve Matrix equation AXB+C=D

I have 4 known matrices:

1. mA = 4x4 -> {{-4, 1, 5, 0}, {1, -1, 1, 2}, {-1, 0, 1, 0}, {1, 17, 0, 1}}
2. mB = 3x3 -> {{-4, 7, 10}, {2, -1, 1}, {1, 0, 1}}
3. mC = 4x3 -> {{1, 2, 2}, {0, -1, 0}, {1, 0, 1}, {0, 0, 0}}
4. mD = 4x3 -> {{1, 2, 3}, {0, 0, 0}, {1, 0, 2}, {0, 0, -1}}

And one unknown:

1. X = 4x3

I want to solve an equation for X in the form mA.X.B+C=D. How can I do this in mathematica? I'm brand new to this language. =)

[EDIT] What I tried and my thought process:

1. AXB+C=D -> AXB = D-C -> D-C=E
2. AXB=E
3. A^-1x(AXB)=A^-1xE
4. I(XB)=A^-1xE
5. XB=A^-1xE
matrixE = matrixD-MatrixC
invA = Inverse[matrixA]

right_equation = invA.matrixE
LinearSolve[matrixB, right_equation]

• Please show the code you have tried and display it along with your matrices in Mathematica format. May 16, 2021 at 16:30
• I'm sorry there was a typo; x is 4x3. I'll update with what I've tried.
– OLGJ
May 16, 2021 at 16:39
• mX = Array[x, {4, 3}]; solution = mX /. First[Solve[(mA . mX . mB + mC) == mD]] May 16, 2021 at 16:39
• The mX statement creates a 4x3 matrix. But what happens afterwards?
– OLGJ
May 16, 2021 at 16:49
• @OLGJ Solve finds the solution(s) and returns lists of lists of rules mapping x[i,j] to some number - in this case there's only one list of rules so we take the First. The mX /. then takes the mX array of x[i,j] variables and applies those rules leaving you with just a 4x3 matrix of numbers. May 16, 2021 at 20:09

mA = {{-4, 1, 5, 0}, {1, -1, 1, 2}, {-1, 0, 1, 0}, {1, 17, 0, 1}};
mB = {{-4, 7, 10}, {2, -1, 1}, {1, 0, 1}};
mC = {{1, 2, 2}, {0, -1, 0}, {1, 0, 1}, {0, 0, 0}};
mD = {{1, 2, 3}, {0, 0, 0}, {1, 0, 2}, {0, 0, -1}};

mX = Inverse[mA] . (mD - mC) . Inverse[mB]
(*    {{-4/7, -141/35, 202/35},
{0,    1/35,    -2/35},
{-3/7, -106/35, 152/35},
{3/7,  89/35,   -118/35}}    *)


test:

mA . mX . mB + mC == mD
(*    True    *)

• I was not aware you could multiply with the inverse of both matrix a & b. But your explanation of it makes sense. Thank you very much Roman!
– OLGJ
May 17, 2021 at 8:01
• The trick that is so trivial that we keep forgetting it is that LEFT and RIGHT inverses are the same: mA . Inverse[mA] == Inverse[mA] . mA == IdentityMatrix[4]. See here for details. May 17, 2021 at 8:29

LyapunovSolve can be used:

zeros = SparseArray[{}, #] &; (* create array of zeros *)
mA = {{-4, 1, 5, 0}, {1, -1, 1, 2}, {-1, 0, 1, 0}, {1, 17, 0, 1}};
mB = {{-4, 7, 10}, {2, -1, 1}, {1, 0, 1}};
mC = {{1, 2, 2}, {0, -1, 0}, {1, 0, 1}, {0, 0, 0}};
mD = {{1, 2, 3}, {0, 0, 0}, {1, 0, 2}, {0, 0, -1}};

mX = LyapunovSolve[{mA, zeros[{4, 4}]}, {zeros[{3, 3}], mB}, mD - mC]
(*
{{-(4/7), -(141/35),  202/35},
{  0,        1/35,   -(2/35)},
{-(3/7), -(106/35),  152/35},
{  3/7,     89/35, -(118/35)}}
*)

mA . mX . mB + mC == mD
(*  True  *)


Update: Alternative

mX = LinearSolve[mA,
Transpose@LinearSolve[Transpose@mB, Transpose[mD - mC]]]

• I have not encountered Lyapunov's method before, but it seems to do the job. Thanks for taking the time Michael!
– OLGJ
May 17, 2021 at 8:05