Why these two similar operations (a sparse matrix w/o transpose times a vector) take different time to finish?

3 views (last 30 days)
n = 1e5;
a = sprand(n,n,5/n);
v = rand(n,1);
tic; for i = 1:1e3, b = a*v; end; toc
tic; for i = 1:1e3, c = a'*v; end; toc
I wanted to test if sparse matrix transpose takes a lot of time, so I used a sparse matrix times a full vector for an example (this is common in engineering problems). The only difference between the two loops is the matrix transpose.
I suprises me that the first loop (a*v) takes 1.28s whereas the second loop (a'*v) takes ONLY 0.18s. Switching the order of the two loops does not change the time cost.
What makes the higu difference? Why is matrix transpose faster than non-transpose?

Accepted Answer

Riccardo Scorretti
Riccardo Scorretti on 6 May 2022
Interesting. I tried to see what happens with matrix of different sizes:
n = round(logspace(1, 6, 30));
t = [];
figure
for k = 1 : numel(n)
a = sprand(n(k),n(k),5/n(k));
v = rand(n(k),1);
tic; for i = 1:1e3, b = a*v; end; t(k,1) = toc;
tic; for i = 1:1e3, c = a'*v; end; t(k,2) = toc;
loglog(n(1:k), t(:,1), 'o-', n(1:k), t(:,2), 's-') ; grid on ;
xlabel('n') ; ylabel('T_{CPU}') ;
legend('a*v', 'a''*v');
drawnow;
end
Clearly something happens close to n = 10000 (more precisely, on my PC between n = 10002 and n = 10003, guess why this value!?):
My very personal hypothesis is that until n = 10002 the two multiplications are executed in the same way, by using a trivial algorithm. After, the transposition and multiplication are combined in a single operation.
The interest of this approach is that, as in MATLAB both full and sparse matrix are stored column-wisely (= like in Fortran), multiplicating the transpose of A by a vector v boils up into executing a kind of column-column multiplication, which is more efficient because if reduces cache misses.
  1 Comment
Zhuoran Han
Zhuoran Han on 10 May 2022
Thank you very much for the explaination!
It might be possible to save a lot of computation time with pre-saved matrix transposes!

Sign in to comment.

More Answers (1)

Sandeep Shrivastava
Sandeep Shrivastava on 6 May 2022
This is an interesting problem to look into! When I calculated the transpose separately and performed only the multiplication in loop, it gave me the same elapsed time as the earlier case. This probably means that it's not the operand that is affecting the performance, but rather:
1. the size of the matrix/vector, and;
2. the operator
They are possibly triggering BLAS routines which handle the large data sets and certain types of operations more efficiently. I think this behavior is not release-dependent.
Script:
n = 1e5;
a = sprand(n,n,5/n);
aT = a';
v = rand(n,1);
tic; for i = 1:1e3, b = a*v; end; toc
tic; for i = 1:1e3, c = aT*v; end; toc
tic; for i = 1:1e3, c = a'*v; end; toc
Output:
Elapsed time is 1.790306 seconds.
Elapsed time is 1.824447 seconds.
Elapsed time is 0.354630 seconds.
Here are some other posts on this forum where BLAS routines may have had a role to play:
  1 Comment
Riccardo Scorretti
Riccardo Scorretti on 6 May 2022
@Sandeep Shrivastava Indeed, in the first link James Tursa writes:
D' * bsxfun(@times,W',D)
"The latter will likely result in D' not being explicitly calculated. Rather, a transpose flag will simply be passed to the dgemm routine in the BLAS call. And the W' in the latter will not result in a data copy, it will be just be a simple reshaped shared data copy of the W vector."
This seems to support the idea that the difference in performances can be explained by the better usage of the cache memory. By explicitly computing and storing aT = a' MATLAB is forced to transpose the matrix physically, and the interpreter (or the JIT compiler) cannot be smart enough to realize that aT is the transposed of a. This would explains also your results.

Sign in to comment.

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!