Matrix Solve for a particular form

Question

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

Answer 1

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}}