6
$\begingroup$

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.

| cite | improve this question | |
$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – Kirill Feb 2 '15 at 22:45
  • 1
    $\begingroup$ 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}$ $\endgroup$ – RogueDodecahedron Feb 3 '15 at 12:43
  • 2
    $\begingroup$ 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? $\endgroup$ – Wolfgang Bangerth Feb 4 '15 at 13:04
2
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.