Suppose I have full rank $n\times n$ matrix $A$ with $\rho(A) < 1$ and I want to find an expression for
$$S = X + A^\top X A + A^{2\top} X A^2 + A^{3\top} X A^3 + \dots$$
where $X$ is an $n\times n$ positive definite matrix. Thus,
$$ S = X + A^\top S A$$
Following this answer to a similar problem, we can make an eigenvalue decomposition of $A$ such that $A=UDU^{-1}$ with $D = \text{diag}(\lambda_1,\dots,\lambda_n)$. Then
$$ S = X + U^{-\top}DU^\top S UDU^{-1}$$
and
$$ Z = U^\top S U = U^\top X U + DU^\top S UD$$
If we then define $T = U^\top XU$ then
$$ Z = T + DZD = T + DTD + D^2TD^2 + D^3TD^3 + \dots$$
which implies that $(i,j)$'th entry of $Z$ is
$$ Z_{ij} = \frac{t_{ij}}{1 - \lambda_i\lambda_j}$$
Having obtained $Z$ then we can find S through
$$S = U^{-\top}ZU^{-1}$$
My questions are
- Is the above correct?
- Is the it a good (fast and numerically stable way) way to compute it?
R code confirming the above in one case
# assign matrices
A <- matrix(c(.8, .4, .1, .5), 2, 2)
X <- matrix(c( 1, .5, .5, 2), 2)
# compute with above formulas
eg <- eigen(A)
U <- eg$vectors
U_t <- t(U)
las <- eg$values
T. <- crossprod(U, X %*% U)
Z <- T. / (1 - tcrossprod(las))
S <- solve(U_t, t(solve(U_t, t(Z))))
# naive solution
nai <- X
for(i in 1:1000)
nai <- crossprod(A, nai %*% A) + X
nai - S
#R [,1] [,2]
#R [1,] -3.55e-15 -4.44e-16
#R [2,] -8.88e-16 0.00e+00
