Series of Specific Matrix Products I'm currently facing a problem where I'm not sure whether a closed form solution exists. Suppose you have two real matrices $A \in \mathbb{R}^{nxn}$ and $B \in \mathbb{R}^{nxm}$ (which arise from a time-discrete linear dynamical system: $x_{t+1}=Ax_t+Bu_t$). Now I would like to have a closed-form expression for
$$
\sum_{i=0}^{N-1}A^{i}BB^{T}\left(A^{i}\right)^{T}
$$ 
where
$$
A^2 = AA\\
A^3 = AAA\\
\vdots
$$
I already thought about looking at the decomposition of $BB^T$ and $A$
$$
BB^{T}=D_{B}\varLambda_{B}D_{B}^{T},\ \ \Lambda_{B}=\left(\begin{array}{cccc}
\lambda_{B_{1}} & 0 & 0 & 0\\
0 & \lambda_{B_{2}} & 0 & 0\\
0 & 0 & \ddots & 0\\
0 & 0 & 0 & \lambda_{B_{n}}
\end{array}\right)
$$
$$
A = D_A \Lambda_A D_A^{-1},\ \ \Lambda_{A}=\left(\begin{array}{cccc}
\lambda_{A_{1}} & 0 & 0 & 0\\
0 & \lambda_{A_{2}} & 0 & 0\\
0 & 0 & \ddots & 0\\
0 & 0 & 0 & \lambda_{A_{n}}
\end{array}\right)
$$
where $\Lambda_A$ and $\Lambda_B$ contain the eigenvalues of $A$ and $BB^T$. Note that while we know that $BB^T$ is symmetric we only know that $A$ is diagonizable. With that we could rewrite
$$
\sum_{i=0}^{N-1}A^{i}BB^{T}\left(A^{i}\right)^{T} \\
 =\sum_{i=0}^{N-1}D_{A}\Lambda_{A}^{i}D_{A}^{-1}D_{B}\varLambda_{B}D_{B}^{T}D_{A}^{-T}\Lambda_{A}^{i}D_{A}^{T}\\
 =D_{A}\left(\sum_{i=0}^{N-1}\Lambda_{A}^{i}D_{A}^{-1}D_{B}\varLambda_{B}D_{B}^{T}D_{A}^{-T}\Lambda_{A}^{i}\right)D_{A}^{T}\\
 =D_{A}\left(\sum_{i=0}^{N-1}\Lambda_{A}^{i}F\Lambda_{A}^{i}\right)D_{A}^{T},\ \ F=D_{A}^{-1}D_{B}\varLambda_{B}D_{B}^{T}D_{A}^{-T}
$$
However, from there on I do not know any tricks that could bring me further. Does anyone have an idea or maybe knows a theorem that could help here? Maybe it is also beneficial to know that the term
$$
\sum_{i=0}^{N-1}A^{i}BB^{T}\left(A^{i}\right)^{T}
$$
arises from the product of the reachability matrix $R$ and its transpose $R^T$ where
$$
R=\left[A^{N-1}B\ A^{N-2}B\ \ldots\ B\right]
$$
Best
Pascal
 A: Here's my go at it. My answer isn't complete, but it may be on the right track or at least food for thought. You have already shown that $A$ and $B$ are simultaneously diagonalizable. This means that there exists a common eigenbasis $\{|\lambda \rangle \}$, which can be made orthonormal. Once you find this basis you can insert unity in the form of $\sum_\lambda |\lambda \rangle \langle \lambda| = 1$. Doing this yields that 
$$ \sum_{i=0}^{N-1}\sum_{\lambda}A^{i}B|\lambda \rangle \langle \lambda|B^{T}\left(A^{i}\right)^{T} = \sum_{i=0}^{N-1}\sum_{\lambda}A^{i}|\lambda \rangle \lambda_B^2\langle \lambda|\left(A^{i}\right)^{T} = \sum_{i=0}^{N-1}\sum_{\lambda}|\lambda \rangle \lambda_{A}^i\lambda_B^2\langle \lambda|$$
Now you can swap the finite sums and get 
$$\sum_{\lambda}|\lambda \rangle \lambda_B^2\langle \lambda|  \sum_{i=0}^{N-1}\lambda_{A}^i = \sum_{\lambda}Tr(A^{N-1})|\lambda \rangle \lambda_B^2\langle \lambda|$$
Now you just have to sum over the outer product of eigenvectors and eigenvalues of $B$, which in principle should be able to be done. 
A: I do not think that you can have a closed form solution in matrix form, but component-wise representation can be obtained pretty easily. 
Let $\Phi^i=\Lambda_A^i F\Lambda_A^i$, where $F$ and $\Lambda_A$ are defined above. Then the $(j,k)$th element of $\Phi^i$ is $\phi_{j,k}^i=\lambda_{A_j}\lambda_{A_k}f_{j,k}$ and we have
$$\sum_{i=0}^N\phi^i_{j,k}=\frac{1-(\lambda_{A_j}\lambda_{A_k})^N}{1-\lambda_{A_j}\lambda_{A_k}}f_{j,k}$$ 
If the system is stable, i.e., $|\lambda_{A_i}|<1$ and $N$ is sufficiently large, we can approximate the finite sum by an infinite one
$$\sum_{i=0}^N \phi^i_{j,k}\approx \sum_{i=0}^\infty \phi^i_{j,k} = \frac{1}{1-\lambda_{A_j}\lambda_{A_k}}f_{j,k}.$$ 
