Raising a large sparse matrix to a large power
8 visualizzazioni (ultimi 30 giorni)
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
or alternatively storing the matrix raised to the (n-1)th power as Mm and recursively computing
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.
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));
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
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
m = m*m;
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
v = m*v;
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.