I want to get an efficient way of computing $P_k$ from $P_0$ that satisfies the following recursion:
$P_k = FP_{k-1}F^T+Q$
Where $P_k$, $F$ and $Q$ are matrices ($Q$ is diagonal, if it changes anything).
We can get the following explicit formula:
$P_k = F^kP_0(F^T)^k+\sum _{i=0}^{k-1}F^iQ(F^T)^i$
Where the first term can be easily computed by diagonalization, but the series in the second term seems tough !
As you noted, the first term can be easily computed by diagonalization. If $F = T\Lambda T^{-1}$, then $$ P_k = T\Lambda^k T^{-1} P_0 T^{-T}\Lambda^k T^T + \sum_{i=0}^k T\Lambda^i T^{-1} Q T^{-T}\Lambda^i T^T $$ You can simplify this slightly to $$ P_k = T\Lambda^k T^{-1} P_0 T^{-T}\Lambda^k T^T + T \left[ \sum_{i=0}^k \Lambda^i \tilde{Q} \Lambda^i \right] T^T $$ where $\tilde Q = T^{-1} Q T^{-T}$. Now note that each term of the series is a matrix $\tilde Q$ with row and column scalings governed by $\Lambda^{i}$. In particular, each element of that matrix is of the form $\lambda_p^i \lambda_q^{i} \tilde{Q}_{pq}$. You can apply the geometric series summation formula to each element separately.