# Raising a large sparse matrix to a large power

8 visualizzazioni (ultimi 30 giorni)
Jon B. il 29 Giu 2023
Commentato: Jon B. il 30 Giu 2023
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.
##### 2 CommentiMostra 1 commento meno recenteNascondi 1 commento meno recente
Jon B. il 30 Giu 2023
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.

Accedi per commentare.

### Risposta accettata

Matt J il 30 Giu 2023
Modificato: Matt J il 30 Giu 2023
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
##### 2 CommentiMostra 1 commento meno recenteNascondi 1 commento meno recente
Jon B. il 30 Giu 2023
Works perfectly! Thanks so much.

Accedi per commentare.

### Più risposte (1)

David Goodmanson il 30 Giu 2023
Modificato: David Goodmanson il 30 Giu 2023
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.
##### 1 CommentoMostra NessunoNascondi Nessuno
Jon B. il 30 Giu 2023
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?

Accedi per commentare.

### Categorie

Scopri di più su Sparse Matrices in Help Center e File Exchange

R2022a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!