# Using Cayley-Hamilton to calculate powers of arbitrary matrices of fixed size

For a symmetric matrix $A$ with eigenvalues $\lambda_i<1$ (with $i=0,\dots,n$), one can write the resolvent ($R(z,A):=(A-zI)^{-1}$) at $z=1$ as a series of powers of $A$ $[1]$:

$R(1,A)=(A-1)^{-1}=-(I+A+A^2+A^3+\dots)$

From Cayley-Hamilton, I know that one can express the power $k\geq n$ of an $n\times n$ matrix $A$ as a linear combination of $A^0,A^1,A^2,\dots, A^{n-1}$. I would like to use this to convert the infinite series of matrix powers above into a finite series of matrix powers (with possibly very ugly coefficients). For this, I need a formula for the $k$-th power of an $n\times n$ matrix of the form

$A^k=\sum_{j=0}^{n-1}\xi_jA^j$

for $k\geq n$. The $\xi_j=\xi_j(c_0,c_1,c_2,\dots,c_n)$ should be function s of the coefficients $c_i$ (with $i=0,\dots,n$) of the $A$'s characteristic polynomial. I have tried to derive this formula by hand and so far it seems to be possible, but very tedious. I would assume that other people have worked on finding such a formula before, but I couldn't find it on the internet.

$[1]$ Barnett, Lionel, Christopher L. Buckley, and Seth Bullock. "Neural complexity and structural connectivity." Physical Review E 79.5 (2009): 051914.

• Well I guess you can express the coefficients of the infinite series $(A-1)^{-1}$ as $\sum_{k=0}^{n-1} \alpha_k A^{k}$ where $\alpha_k$ is the solution to the system of equations $(\lambda_i - 1)^{-1} = \sum_{k=0}^{n-1} \alpha_k \lambda_i^k, i=1,\ldots,n$. Then rewrite each $A^k$ in terms of the characteristic polynomial starting from $A^{n-1}$ and rearranging as you go and finally collecting terms - sounds messy though – Nadiels Aug 24 '17 at 12:27

Obviously the expansion must be unique, so we can find it however we like. If $$p(t) = \sum_{k=0}^n c_k t^k$$ is the characteristic polynomial of $A$, then if we write $B=A-I$, then $A=I+B$ satisfies $$\sum_{k=0}^n c_k (I+B)^k = 0$$ by Cayley–Hamilton. We can then expand to get $$0 = \sum_{k=0}^n c_k \sum_{m=0}^k \binom{k}{m} B^m = \sum_{m=0}^n B^m \sum_{k=m}^n c_k \binom{k}{m} = \sum_{m=0}^n \alpha_m B^m$$ by changing the order of summation, and the $\alpha_m = \sum_{k=m}^n c_k \binom{k}{m}$ are constants. So $$\alpha_0 = -\sum_{m=1}^n B^m \alpha_m.$$ $\alpha_0$ is nonzero since $B$ is supposed to be invertible, so multiplying both sides by $B^{-1}\alpha^{-1}$, we have $$B^{-1} = \sum_{m=1}^n -\frac{\alpha_m}{\alpha_0} B^{m-1} .$$ Now we can go back to $A$ using $B=A-I$, which gives $$(A-I)^{-1} = \sum_{m=0}^{n-1} -\frac{\alpha_{m+1}}{\alpha_0} (A-I)^{m} = \sum_{m=0}^{n-1} -\frac{\alpha_{m+1}}{\alpha_0} \sum_{j=0}^m (-1)^{m-j} \binom{m}{j} A^j \\ = \sum_{j=0}^{n-1} A^j \frac{-1}{\alpha_0}\sum_{m=j}^{n-1} (-1)^{m-j} \binom{m}{j} \alpha_{m+1} \\ = \sum_{j=0}^{n-1} A^j \frac{-1}{\alpha_0}\sum_{m=j}^{n-1} (-1)^{m-j} \binom{m}{j} \sum_{k=m+1}^n c_k \binom{k}{m+1}.$$ Therefore $$(A-I)^{-1} = \sum_{j=0}^{n-1} \beta_j A^j,$$ with coefficients $$\beta_j = -\left( \sum_{k=0}^n c_k \right)^{-1} \sum_{m=j}^{n-1} (-1)^{m-j} \binom{m}{j} \sum_{k=m+1}^n c_k \binom{k}{m+1}$$ It may be possible to simplify this expression further by changing the order of summation again and using some binomial identities.
$$\sum_{m=j}^{n-1} (-1)^{m-j+1} \binom{m}{j} \sum_{k=m+1}^n c_k \binom{k}{m+1} = (-1)^j \sum_{k=j+1}^n c_k \sum_{M=j+1}^{k} (-1)^{M} \binom{M-1}{j} \binom{k}{M} = -\sum_{k=j+1}^n c_k$$ by Euler's formula on sums of $k$th differences of powers, since $j<k$.
Actually, since everything commutes, we can do some formal matrix division to get this result more simply: $$0 = (A-I)^{-1}(c_n A^n+ c_{n-1} A^{n-1} + \dotsb + c_1 A + c_0 I ).$$ Then $$(A-I)^{-1}A^k = A^{k-1}+A^k+\dotsb+I + (A-I)^{-1},$$ and so we have $$\left( \sum_{k=0}^n c_k \right) (A-I)^{-1} = -\sum_{m=0}^{n} A^{m} \sum_{k=m+1}^n c_k,$$ and again we find $$(A-I)^{-1} = \sum_{m=0}^{n-1} A^m \frac{\sum_{k=m+1}^n c_k}{\sum_{k=0}^n c_k}$$