(Asymmetric) matrix power series in closed form: $\sum_{i=0}^{\infty} A^i \left(A^i\right)^{\top}={?}$

Question

Let $A\in \mathbb{S}^{N\times N}$ be a symmetric, real and stable matrix, i.e., $\rho(A)<1$, where $\rho(A)$ stands for the spectral radius of $A$. Then, $$\sum\limits_{i=0}^{\infty} A^{2i}=\left( I-A^2 \right)^{-1}.\tag{$\star$}\label{star}$$

Now, let $A$ be asymmetric. It is not clear whether the power series $$F(A):=\sum_{i=0}^{\infty} A^i (A^i)^{\top}$$ admits a compact closed-form expression as in equation \eqref{star} (when it converges). If the matrix $A$ is normal — i.e., it commutes with its transpose — and stable, then we have that $F(A)=\left(I-AA^{\top}\right)^{-1}$. In general, that is not the case.

Framework. Let $A$ be a random matrix generated as follows. Let $\widetilde{A}\in \mathcal{G}(N,p)$ be the adjacency matrix associated with the realization of a random graph over $N$ nodes with probability $p$ of arrow drawing, i.e., the entries $\widetilde{A}_{ij}$ of $\widetilde{A}$ are i.i.d. $\mathsf{Bernoulli}(p)$ for all $i\neq j$ — this is often referred to as a binomial random graph. Set $\widetilde{A}_{ii}=0$ for all $i$. Let $A$ be obtained from $\widetilde{A}$ by: i) [off-diagonal] normalizing the entries of $\widetilde{A}$ as $A_{ij}:=\alpha_1 \widetilde{A}_{ij}/d_{\max}$, for all $i\neq j$, where $d_{\max}$ is the maximum in-flow degree of the random graph and $0<\alpha_1<1$ is some constant; ii) [diagonal] setting the diagonal terms as $A_{ii}:=\alpha-\sum_{k\neq i} \widetilde{A}_{ik}$, where $0<\alpha_1\leq\alpha<1$. In other words, the rows of $A$ sum to $\alpha$ and its support is given by the realization of the binomial random graph on $N$ nodes.

Question. Is it true that $$\mathbb{P}\left(\lVert\sum_{i=0}^{\infty}A^i \left(A^i\right)^{\top} - \left(I-AA^{\top}\right)^{-1}\rVert_\text{max}>\varepsilon\right)\overset{N\rightarrow \infty}\longrightarrow 0?$$

Is there any known result in this direction? Any reference would be appreciated. The specifics of the random nature of $A$ are not important other than that of $A$ being stable and randomly generated (so that it is more natural to consider the limit $N\rightarrow\infty$).

Some numerical experiments are pointing in the direction of a concentration of $\sum_{i=0}^{\infty} A^i \left(A^i\right)^{\top}$ about $\left(I-AA^{\top}\right)^{-1}$ as the dimension $N$ grows large (I may attach them later, if necessary) under the referred framework.

Context. $F(A)$ is the covariance matrix associated with the time-series reflecting the state-evolution of certain linear stochastic dynamical systems. Having a closed form expression for the covariance in terms of $A$ is helpful for model identification purposes.

Answer 1

Let me answer the first question:

Q: Does $F=\sum_{i=0}^{\infty} A^i (A^i)^{\top}$ for non-symmetric $A$ have a closed-form solution analogous to the solution $\sum_{i=0}^{\infty} A^{2i}=\left( I-A^2 \right)^{-1}$ for a symmetric $A$?

A: Yes, in terms of the vectorization operation: $$\operatorname{vec}(F)=(I-{A} \otimes A)^{-1}\cdot\operatorname{vec}(I).\tag{$\ast$}\label{ast}$$ This does require the inversion of an $N^2\times N^2$ matrix, rather than the inversion of an $N\times N$ matrix as in the case of a symmetric $A$.

The vectorization formula follows because $F$ solves the Lyapunov equation $$AFA^\top=F-I.$$

I do not have an answer to the second question, whether \eqref{ast} converges for large $N$ to $(I-AA^\top)^{-1}$.