A "geometric'' infinite sum of matrices The sum $$ I + A + A^2 + A^3 + \cdots $$ equals $$(I-A)^{-1}$$ under the assumption $\rho(A)<1$, which is necessary to make the sum converge. My question: what does the sum 
$$ I + A^T A + (A^2)^T A^2 + (A^3)^T A^3 + \cdots$$ equal under the same assumption? Is there a similarly neat expression?
It is true that this last sum converges under the assumption $\rho(A)<1$. The ``obvious'' guess $(I-A^T A)^{-1}$ for the sum is not true because $ (A^2)^T A^2 \neq (A^T A)^2$; in fact, this guess does not even make sense because $\rho(A)<1$ does not rule out $\rho(A^T A)=1$. 
 A: If $F:M_n\to M_n$ is defined by $F(X)=A^\mathrm{T}XA$ and $\mathrm{id}:M_n\to M_n$ is the identity map, then your sum is $(\mathrm{id} - F)^{-1}(I)$.
A: Let $S$ be the solution to the discrete-time Lyapunov equation $A^TSA=S-I$. Note that $S$ exists and is unique because
$$
A^TSA=S-I\ \Leftrightarrow\ (I - A^T\otimes A^T)\textrm{vec}(S)=\textrm{vec}(I)
$$
and $I - A^T\otimes A^T$ is invertible owing to the fact that $\rho(A^T\otimes A^T)=\rho(A)^2<1$.
It follows that if the infinite series $I + A^T A + (A^2)^T A^2 + (A^3)^T A^3 + \cdots$ converges at all, its limit must be $S$.
But does it really converge? Note that for some matrix $A$ such as $\frac1{\sqrt{2}}\pmatrix{1&1\\ 0&1}$ we have $\rho(A)<1$ but also $\max\left\{\|A\|,\|A^T\|\right\}>1$ for every submultiplicative matrix norm. So, we cannot argue that $\|A\|\|A^T\|<1$ and pretend that $\|\sum_{k=0}^\infty(A^j)^TA^j\|\le\sum_{j=0}^\infty(\|A\|\|A^T\|)^j<\infty$ (the first inequality is always valid but the second one is not). The usual trick for proving the convergence of a Neumann series using submultiplicative norm is not applicable here.
Yet, we can still easily prove that the partial sums are convergent to $S$. Since $A^TSA=S-I$, one can prove by mathematical induction that
$$
S-\sum_{j=0}^{n-1}(A^j)^TA^j=(A^n)^TSA^n
$$
for every $n\ge1$. As $\rho(A)<1$, we have $A^n\rightarrow0$ and in turn $(A^n)^TSA^n\rightarrow0$ when $n$ approaches infinity. Hence the infinite series $\sum_{j=0}^\infty(A^j)^TA^j$ does converge to $S$.
Alternatively, we can vectorise the $n$-th partial sum of the infinite series as follows:
$$
\operatorname{vec}\left(\sum_{j=0}^n(A^j)^T A^j\right)
=\sum_{j=0}^n\left((A^j)^T \otimes (A^j)^T\right)\operatorname{vec}(I)
=\sum_{j=0}^n\left(A^T\otimes A^T\right)^j\operatorname{vec}(I).\tag{$\ast$}
$$
Now we get a geometric sum on the RHS of $(\ast)$. This effectively resurrects the Neumann series argument. As $\rho\left(A^T\otimes A^T\right)=\rho(A)^2<1$, we have $\|A^T\otimes A^T\|<1$ for some submultiplicative norm. The RHS of $(\ast)$ thus converges when $n\to\infty$ and the LHS converges too.
A: FWIW: Call $S$ the sum, provided it makes sense. Then $A^TSA=S-I$, so the series converges to a root of this equation.
If $n=2$, then the determinant of the endomorphism $S\mapsto A^TSA-S$ of $M_n(\mathbb C)$ is $(\det A-1)^2\det(A^2-I)$, so for most $A$s, the solution is unique. For larger $n$, I have no idea.
