# Is there a faster method to compute the geometric series of a matrix?

I want to calculate the geometric series of a matrix $A$:

$$S=I+A+A^2+\dots+A^n$$

and then apply to a vector $v$, $Sv$.

I've done it in Matlab with a loop and I think it's quite efficient applying the matrix $A$ to $v$ at each step instead of calculating $S$ first.

Anyway, I tried to compute the series using $S=(I-A^n)(I-A)^{-1}$, but when calculating the inverse I get an error because $A$ is nearly singular.

Is the an efficient way to compute the inverse in this case? I thought it would be better than a loop since Matlab is quite good manipulating matrices.

I know that $A$ has dominant diagonal and tridiagonal.

• Is $A$ singular or $I-A$? Also, Horner's rule takes $n$ matrix-vector multiplications here, which is the same it takes to multiply by $I-A^n$ (unless you compute $A^n$ explicitly in $O(\log n)$ multiplications, but that matrix would have much wider bands), so I'm not yet sure why $(I-A^n)(I-A)^{-1}$ could be expected to be faster. – Kirill Feb 2 '15 at 22:45
• Do you calculate the inverse in Matlab with inv(A)? If so, this is a bad idea, as explained here: de.mathworks.com/help/matlab/ref/inv.html . It is better to calculate x=inv(1-A)*v by x=(1-A)\v. Edit: there is a mistake in the formula for $S$. It is $S=(1-A^{n+1}) \cdot (1-A)^{-1}$ – RogueDodecahedron Feb 3 '15 at 12:43
• You say that you compute the application of $S$ by a loop, but it's not clear to me how you do it. Do you factor your polynomial $I+A+A^2+\ldots$ so that you only need $n$ matrix-vector products, or do you multiply your vector with each term individually, requiring $\frac 12 n^2$ matrix-vector products? – Wolfgang Bangerth Feb 4 '15 at 13:04

Perhaps it would help to decompose $v$ in terms of the (normalized) eigenvectors $w_i$ of the matrix $A$. Let the corresponding eigenvalues be $\lambda_i$, so that $A \cdot w_i = \lambda_i w_i$.
If $c_i = <w_i, v>$ then $(I + A + A^2 + ... + A^n) \cdot v = \sum_i c_i (1+\lambda_i + \lambda_i^2 + ... + \lambda_i^n) w_i$
For the eigenvalues that have $|\lambda_i - 1| > \epsilon$ you can use $(1-\lambda^{n+1}_i)/(1-\lambda_i)$ to quickly compute your sum, and for the others you can for example do an explicit summation.