Azzera filtri
Azzera filtri

Raising a large sparse matrix to a large power

16 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
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 Commenti
Steven Lord
Steven Lord il 30 Giu 2023
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.
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
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 Commenti
Matt J
Matt J il 30 Giu 2023
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.
Jon B. il 30 Giu 2023
Works perfectly! Thanks so much.

Accedi per commentare.

Più risposte (1)

David Goodmanson
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 Commento
Jon B.
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

Prodotti


Release

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by