# Matrix Solve for a particular form

Suppose that we have two matrices $\mathbf A$ and $\mathbf B$ which are known to us and both of them are square matrices of the same dimensions.

Now we want to find a square matrix $\mathbf C$ that solves following equation:

$$\mathbf C=\mathbf A\mathbf C\mathbf A^\top+\mathbf B$$

How can we find this $\mathbf C$ matrix in Mathematica?

Just for additional information, it is for finding the state-space unconditional covariance matrix

Thanks

LyapunovSolve[] is designed for this:

amat = Array[a, {2, 2}]; bmat = Array[b, {2, 2}];
cmat = LyapunovSolve[{amat, -IdentityMatrix[Length[amat]]},
{IdentityMatrix[Length[amat]], Transpose[amat]}, -bmat];


Check:

amat.cmat.Transpose[amat] + bmat - cmat // Simplify
{{0, 0}, {0, 0}}


Alternatively, you can reformulate as a Kronecker product linear system:

cmat2 = Partition[LinearSolve[KroneckerProduct[amat, amat] -
IdentityMatrix[Length[amat]^2],
-Flatten[bmat]], Length[amat]];


Check:

amat.cmat2.Transpose[amat] + bmat - cmat2 // Simplify
{{0, 0}, {0, 0}}

• Thanks for your help. I checked Mathematica documentation and you are right. I should use LyapunovSolve[], but since it is a discrete version, I should use DiscreteLyapunovSolve[]. – George Sep 11 '17 at 20:42