Raising a large sparse matrix to a large power

I am currently implementing a recursion that involves the n-fold application of a rather large sparse matrix (ballpark of 20,000 by 20,000) to a particular vector. So, basically, my goal is to raise a large matrix to the n-th power where over the course of the recursion n reaches values >100. Unfortunately, computing the matrix power via

M^n

or alternatively storing the matrix raised to the (n-1)th power as Mm and recursively computing

M*Mm

becomes succesively slower and once n>10 takes minutes for every single step. The issue would seem to be that the number of non-zero elements grows at a rather speedy rate.

The paper proposing this very recursion, however, emphasizes that - due to the sparse nature of the matrix - the entire recursion (n=1..100) is computationally inexpensive and should run in roughly 0.17 seconds (implemented for similarly dimensioned objects).

Am I overlooking something that should be obvious? Are there more efficient ways to raise a sparse matrix to a large power? Is Matlab notoriously slow when it comes to handeling sparse matrices, so that my best course of action would be to resort to Python, Julia, etc.?

Thanks so much.

[Steven Lord]

How sparse is M? Is it a non-sparsely populated sparse matrix? What about the same question for M^2, M^3, etc.?

F1 = eye(10);
M1 = speye(10);
F2 = ones(10);
M2 = sparse(ones(10));
sparsity1 = nnz(M1)./numel(M1)

sparsity1 = 0.1000

sparsity2 = nnz(M2)./numel(M2)

sparsity2 = 1

whos F1 M1 F2 M2

  Name       Size            Bytes  Class     Attributes

  F1        10x10              800  double              
  F2        10x10              800  double              
  M1        10x10              248  double    sparse    
  M2        10x10             1688  double    sparse    

Note that even though M2 is stored as a sparse matrix, it's actually fully populated and so consumes more memory than if you had stored it as a full matrix!

Does your M matrix have special structure in the non-zero elements or are the non-zero elements more or less scattered?

spy(M1)

Are you trying to store each of those sparse matrices or are you overwriting the matrix each time to minimize the number of copies of it that are in memory?

[Jon B.]

Many thanks for your comment. The sparsity of my matrix is at 0.002, while the sparsity of M^2 is 0.025 and the sparsity of M^3 is 0.1303. So I guess that I am approaching non-sparsity fairly rapidly (which is probably at the core of my problem). In terms of its non-zero element structure, the matrix looks as follows:

I am currently overwriting the matrix in each iteration. I have, however, also tried to precompute a set of >100 matrices to be used as an ingredient in my recursion. Unfortunately, but predictably, to no avail.

Answer 1

[Matt J]

The actual recursion follows the structure x(n+1) = x(n) + AM^nBy(n) where each instance of x is a 5-by-3 matrix, A is 5-by-20,000, M is my large sparse matrix, B is 20,000-by-20,000, and each instance of y is 20,000-by-3.

C = A;
for n = 2:100 
  C = C*M;
  x(:,:,n) = x(:,:,n-1) + C*(B*y(:,:,n-1));
end

[Matt J]

Also, if y is known in advance of the loop, you should capitalize on pagemtimes,

By=pagemtimes(B,y);
C = A;
for n = 2:100 
  C = C*M;
  x(:,:,n) = x(:,:,n-1) + C*By(:,:,n-1);
end

[Jon B.]

Works perfectly! Thanks so much.

Answer 2

[David Goodmanson]

Hi if you are operating on 'a particular vector' v then you do not have to iterate as in

m = m*m             % iterate to form m^n

but rather

v = m*v               % iterate to form m^n*v

which is much much faster. In the following example the 20000x20000 matrix has 1e5 nonzero elements. One can see that m^n has very fast growth of nonzero elements whereas m^n*v does not. And even if there is growth in the latter case, there is only room for 20 thousand nonzero values compared to 400 million nonzero values when computing m^n.

N = 20000;
nnzm = 1e5;
nnzp = 1e3;
i = randi(N,nnzm,1);
j = randi(N,nnzm,1);
a = rand(nnzm,1);
m = sparse(i,j,a,N,N,1e6);
p = randi(N,nnzp,1);
q = ones(size(p));
b = rand(nnzp,1);
v = sparse(p,q,b,N,1,1e4);
for k = 1:3
  tic
  m = m*m;
  toc
end
   Elapsed time is 0.011158 seconds.
   Elapsed time is 0.145417 seconds.
   Elapsed time is 8.828610 seconds.         % m^4 is it

m = sparse(i,j,a,N,N,1e6);                   % reset m
for k = 1:6
  tic
  v = m*v;
  toc 
end
   Elapsed time is 0.000831 seconds.
   Elapsed time is 0.000562 seconds.
   Elapsed time is 0.000865 seconds.
   Elapsed time is 0.001143 seconds.
   Elapsed time is 0.001612 seconds.
   Elapsed time is 0.001126 seconds.

[Jon B.]

Thank you so much for you answer! The actual recursion follows the structure x(n+1) = x(n) + AM^nBy(n) where each instance of x is a 5-by-3 matrix, A is 5-by-20,000, M is my large sparse matrix, B is 20,000-by-20,000, and each instance of y is 20,000-by-3. So, currently, my (stylized) code looks something like this:

Mn = speye(20000);
for n = 2:100 
  Mn = M*Mn;
  x(:,:,n) = x(:,:,n-1) + A*Mn*B*y(:,:,n-1);
end

Since y changes with n, I'm afraid I don't see an immediate way to resolve the growth problem. M^nB remains 20,000-by-20,000, after all. Is there any way to do something clever about that?