Matrix geometric series with multiplying from both sides I am interested in solving following infinite series $J + M^T J M + M^T M^T J M M$ + ...
I know that this series converges in my case.
I would like to derive a theorem analogous to following one.

Let A be a square matrix. If $|A|<1$, the series $$S=I+A+A^2+\cdots$$
  converges to $(I-A)^{-1}$.

This is what comes to my mind:
Let us denote $S$ the sum of above mentioned series.
Then the following equation holds: $$S = J + M^T S M$$
I know that I can solve this equation by treating every element of S as a variable. If shape of $S$ is $n*n$ then I would need to solve the system of equations with $n^2$ variables.
My question is: Is there a faster way how to compute this? And what additional assumptions do I need?
$J$ is symetric possitive definite matrix and I already have it factored to Cholesky decomposition. So the complexity of factoring this matrix doesn't count. I also don't mind if you make assumptions about the matrix $M$ (e.g. that it is diagonalizable and that I already have it factored.)
What else do I need to to find the sum quickly? And how to do it?
 A: What the series converges to, when it converges, is the solution of your equation $S = J + M^T S M$.  And the whole point of it is that if
the spectral radius of $M$ is small, this series converges quickly and provides a fast way to approximate that solution.  
If you have already diagonalized $M$, it is very fast to compute terms of the series: if $M = V^{-1} D V$ where $D$ is diagonal with diagonal elements $d_j$, then 
$$ (M^T)^n J M^n = V^T D^n (V^{-1})^T J V^{-1} D^n V $$
The $(i,j)$ matrix element of this is 
$$ d_i^n ((V^{-1})^T J V^{-1})_{ij} d_j^n $$
Summing over all $n$, assuming $|d_i d_j| < 1$, we get
$$ \dfrac{((V^{-1})^T J V^{-1})_{ij}}{1-d_i d_j} $$
Thus you can quickly compute the result by
 $$ \eqalign{Q &= (V^{-1})^T J V^{-1} \cr
             R_{ij} &= \frac{Q_{ij}}{1-d_i d_j} \ \text{for each $i,j$}\cr
             S &= V^T R V\cr} $$
A: You could vectorize the equation (utilizing the Kronecker product $\otimes$)
$${\rm vec}(S) = {\rm vec}(J) + (M^T\otimes M^T)\,{\rm vec}(S)$$
Then solve it
$${\rm vec}(S) = (I\otimes I - M^T\otimes M^T)^{-1}\,{\rm vec}(J)$$
And finally de-vectorize it to recover a matrix-valued result
$$S = {\rm devec}\Big((I\otimes I - M^T\otimes M^T)^{-1}\,{\rm vec}(J)\Big)$$
