A truncated "geometric" matrix series Let $A$ be an $n\times n$ matrix, $B$ be an $n\times m$ matrix, $C$ an $m \times m$ matrix, and consider the sum
$$\sum_{k = 0}^{N-1} A^k B C^k.$$
Is there any smart way to rewrite this sum in a way similar to the partial sums of the geometric series; namely for $a,b \in \mathbb{R},$
$$\sum_{k = 0}^{N-1} b a^k = b \frac{1-a^N}{1-a}?$$
 A: This is a long comment. Write $S$ for this sum. If we just directly imitate the usual geometric series argument we are led to consider
$$ASC = \sum_{k=1}^N A^k B C^k = S - B + A^N B C^N$$
so $ASC - S = A^N B C^N - B$. And here we're a bit stuck: $ASC - S$ just isn't a multiplication of $S$ by any matrix on either the left or the right so we can't hope to perform a simple matrix division here. We just need to talk about the inverse, if it exists, of the linear map $S \mapsto ASC - S$. But at least this is a system of linear equations whose solution, which we can hope is unique, gives the coefficients of $S$ so that's not too bad.
A: With meager loss of generality, let's first assume that $A,C$ are diagonalizable so that $A=RDR^{-1},C=SES^{-1}$ where $R,S$ are invertible and $D,E$ are diagonal.
Then $$\sum_{k=0}^{N-1}A^kBC^k=\sum_{k=0}^{N-1}RD^kR^{-1}BSE^kS^{-1}$$
$$=R(\sum_{k=0}^{N-1}D^kR^{-1}BSE^k)S^{-1}.$$
Suppose now that $R^{-1}BS=(\alpha_{i,j})_{i,j}$ and that the diagonal entries in $D$ are $(\delta_1,\dots,\delta_n)$ and the diagonal entries in $E$ are $(\epsilon_1,\dots,\epsilon_m)$.
Then $$\sum_{k=0}^{N-1}A^kBC^k=R(\sum_{k=0}^{N-1}\delta_i^k\alpha_{i,j}\epsilon_j^k)_{i,j}S^{-1}$$
$$=R\cdot\Big(\frac{(\delta_i\epsilon_j)^N-1}{\delta_i\epsilon_j-1}\cdot\alpha_{i,j}\Big)_{i,j}\cdot S^{-1}.$$
It looks like there is a not quite as nice expression for our sum when $A,C$ are not diagonalizable, and $D,E$ are in Jordan normal form.
Generalization:
Suppose that $f(x,y)$ is a polynomial (or an analytic function when we assume convergence), and $f(x,y)=\sum_{k,l}b_{k,l}x^ky^l$. Then by using the same reasoning and matrices as before, we have
$$\sum_{k,l}b_{k,l}A^kBC^l=R\cdot(f(\delta_i,\epsilon_j)\cdot\alpha_{i,j})_{i,j}\cdot S^{-1}.$$
