9

What is the best way to calculate sum of matrices such as A^i + A^(i+1) + A^i+2........A^n for very large n?

I have thought of two possible ways:

1) Use logarithmic matrix exponentiation(LME) for A^i, then calculate the subsequent matrices by multiplying by A.

Problem: Doesn't really take advantage of the LME algorithm as i am using it only for the lowest power!!

2)Use LME for finding A^n and memoize the intermediate calculations.

Problem: Too much space required for large n.

Is there a third way?

ishan3243
  • 1,870
  • 4
  • 30
  • 49

3 Answers3

13

Notice that:

A + A^2 = A(I + A)
A + A^2 + A^3 = A(I + A) + A^3
A + A^2 + A^3 + A^4 = (A + A^2)(I + A^2)
                    = A(I + A)(I + A^2)

Let

B(n) = A + ... + A^n

We have:

B(1) = A
B(n) = B(n / 2) * (I + A^(n / 2)) if n is even
B(n) = B(n / 2) * (I + A^(n / 2)) + A^n if n is odd

So you will be doing a logarithmic number of steps and there is no need to compute inverses.

While a direct implementation will lead to a (log n)^2 factor, you can keep it at log n by computing the powers of A as you compute B.

IVlad
  • 43,099
  • 13
  • 111
  • 179
  • Calculating the inverse is about O(m^2.4) where m is the size of the matrix. Multiplication cost is O(m^2.4) too, and you are doing way more multiplication then I do. The only advantage of this approach is it works even if you can't invert I-A, but I'm quite sure that it's worth checking first if you can invert it, and if so use my solution. – Save Sep 19 '13 at 09:50
  • @Save - finding the inverse is also significantly more difficult to implement correctly. Even checking if it exists requires computing the determinant, which is not a trivial problem. – IVlad Sep 19 '13 at 09:52
  • There are tons of libraries to do taht for you. And a trivial implementation of both multiplication and inversion costs O(n^3) and it's quite easy to do (Gaussian Elimination), so, again, I'd suggest about checking first and then decide the method. – Save Sep 19 '13 at 09:54
  • It might work, but inverses will require floating point numbers. Consider that you might not get exact results. – IVlad Sep 19 '13 at 09:57
  • Correct me if I'm wrong, but for large n, calculating A^n is still going to be a real efficiency killer at O(m^(2.4*n)) where m is the matrix size (assuming the 2.4 above is correct) – fresidue Sep 19 '13 at 10:50
  • @IVlad basically in this solution you are computing A^n then A^(n/2) then A^(n/4) and so on.....Wouldn't this lead to O(log(n)+log(n/2)+log(n/4).....) complexity? Ok i get that it is O((log n)^2) but i don't understand how to convert to log n – ishan3243 Sep 19 '13 at 10:56
  • @user2545792 - even if you leave it at `O((log n)^2)` it should be very fast. To reduce it, notice that if you have `A^(n/2)` computed and stored in, say `t`, then `A^n` is just `t*t`, so an `O(1)` operation. – IVlad Sep 19 '13 at 11:01
  • 1
    @fresidue - The complexity is `O(m^3 * log n)` if using the naive matrix multiplication and `O(m^2.4 * log n)` using the optimal multiplication algorithm. As long as the matrices are small enough, this is going to be very efficient, regardless of how large `n` is. – IVlad Sep 19 '13 at 11:02
  • What does $B(n/2)$ mean when $n$ is odd? – user32849 May 31 '18 at 12:47
  • @user32849 integer division, `5 / 2 = 2`. – IVlad May 31 '18 at 21:20
4

You can use the fact that the sum to n of a geometric series of matrices equals:

S_n = (I-A)^(-1) (I-A^n)

and, since you don't start from 0, you can simply calculate:

result = S_n - S_i

where i is your starting index.

Save
  • 11,450
  • 1
  • 18
  • 23
1

Why not just diagonalize the matrix to make the multiplication cheap.

edit:

As long as the matrix is nonsingular you should be able to find a diagonal representation D of matrix A such that A = PDP^-1 where P is made up of the eigenvectors of A, and D has the eigenvalues of A along the diagonal. Getting D^m = D*D^(m-1) is cheap since it's you're only multiplying along the diagonal (i.e. the same number of multiplications as the dimension of the matrix)

Getting S(m)=S(m-1)+D^m is also cheap since you're only adding diagonal elements.

Then you have

A^i + A^(i+1) + A^i+2........A^n = P(D^i + D^(i+1) + D^i+2........D^n)P^-1 = P( S(n) - S(i) )P^-1

The only difficult thing is finding P and P^-1

fresidue
  • 888
  • 8
  • 10