# A "geometric''-ish infinite sum of matrices

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

1. Is the above correct?
2. 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

• For an alternate approach: if we use vectorization and the Kronecker product as in user1551's answer to the question, then we can compute the desired matrix explicitly as $$\operatorname{vec}^{-1}[(I - A^T \otimes A^T)^{-1}\operatorname{vec}(X)]$$ Jul 17, 2018 at 10:00
• Another useful observation (should you want to search existing literature) is that $S = X + A^TSA$ is an instance of the discrete Lyapunov equation. Jul 17, 2018 at 10:04
• @Omnomnomnom Are there advantages of this method? I end up with a $n^2 \times n^2$ matrix but I guess that there may be an easy way to do the computation? Jul 18, 2018 at 9:37
• I really just meant that, at worst, we can do the problem by inverting an $n^2 \times n^2$ matrix. I don't know enough about the complexity of computing eigendecompositions to know which of our two methods will be faster Jul 18, 2018 at 10:43