How to solve a certain recurrent matrix equation I am trying to find the steady-state covariance of a vector that repeatedly undergoes a linear transformation with noise. Specifically, the following recurrent relation should hold:
$$cov(AX + \epsilon) = cov(X)$$
where $A \in \mathbb{R}^{m\times m}$ is a constant matrix and $X$ and $\epsilon$ are independent random vectors; $\epsilon \sim \cal{N}(\bf{0}, \Sigma)$. Due to the independence, we can write:
$$ Acov(X)A^T + \Sigma = cov(X)$$
From this, is it possible to solve for $cov(X)$? 
 A: Nobody's answered yet, but I found a solution:
Let $C := cov(X)$ and let $\odot$ mean elementwise multiplication.
This identity holds trivially:
$$\mathbf{v}^TC\mathbf{w} = \sum_{i, j} (\mathbf{vw}^T \odot C)_{ij} = \sum_{i,j} \mathbf{v}_i\mathbf{w}_jC_{ij}$$
Since each entry of $ACA^T$ takes this form, where v is the ith row and w is the jth row of $A^T$, we can make a system of $m^2$ equations to solve for each entry of $C$.
From the last equation in the question:
$$C_{ij} = \Sigma_{ij} + \sum_{k=1}^m \sum_{l=1}^m A_{ik}A_{jl}C_{kl}$$
$$\Sigma_{ij} = C_{ij} - \sum_{k=1}^m \sum_{l=1}^m A_{ik}A_{jl}C_{kl}$$
Let $r(i, j) := im + j$, $\mathbf{\sigma} \in \mathbb{R}^{m^2}$ where $\mathbf{\sigma}_{r(i, j)} = \Sigma_{ij}$, and $\mathbf{c} \in \mathbb{R}^{m^2}$ where $\mathbf{c}_{r(i, j)} = C_{ij}$. Then we can write:
$$\mathbf{\sigma} = D\mathbf{c}$$
The entries of $D$ are:
$$D_{r(i,j)r(k,l)} = \delta_{r(i,j)r(k,l)} - A_{ik}A_{jl}$$
Then we can solve for c and reformat it as a matrix.
(I actually implemented this in Python using numpy.linalg.tensorsolve, which allowed me to avoid reshaping the matrices into vectors and back again.)
