Matrix geometric series with multiplying from both sides

Question

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?

Answer 1

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

Answer 2

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)$$