It is usually faster to do matrix exponentiation A ^ n
rather than repeated multiplication A * A * A * ... n times
where the exponent n
is some positive integer in this case. This is indeed the case for full matrices, but sparse matrix exponentiation is unexpectedly slower than repeated multiplication.
I have an adjacency matrix representing a large and sparse state transition graph:
>> whos adj
Variables visible from the current scope:
variables in scope: top scope
Attr Name Size Bytes Class
==== ==== ==== ===== =====
s adj 25515x25515 1047824 double
Total is 651015225 elements using 1047824 bytes
>> numel (adj)
ans = 651015225
>> nnz (adj)
ans = 52731
>> unique (adj(:))
ans =
0
1
(It is represented as double
instead of logical
because it will be multiplied soon.)
It is a wellknown property of 0/1 adjacency matrices that the (i,j)th
element of adj ^ n
gives the number of paths of length n
from i
to j
.
If we call adj ^ n
and measure the time:
>> t = zeros(1,15); for i = 1:15, tic; tmp = adj ^ i; t(i) = toc; disp([i t(i)]); end
1 0.07019686698913574
2 0.003432989120483398
3 0.0135948657989502
4 0.03022098541259766
5 0.1219630241394043
6 0.2582740783691406
7 0.5015640258789062
8 0.906404972076416
9 3.269675016403198
10 5.212830066680908
11 7.401607036590576
12 11.35441994667053
13 16.15549778938293
14 27.15614080429077
15 51.56035304069519
If we do adj * adj * adj ... n times
and measure the time:
>> u = zeros(1,15); for i = 1:15, tic; tmp = eval(["adj" repmat(" * adj", 1, i1)]); u(i) = toc; disp([i u(i)]); end
1 0.07631683349609375
2 0.003623008728027344
3 0.01030302047729492
4 0.02698016166687012
5 0.06217288970947266
6 0.1419768333435059
7 0.3206641674041748
8 0.7241461277008057
9 1.282343149185181
10 2.265716075897217
11 3.962852001190186
12 6.107408046722412
13 7.560194969177246
14 9.875128030776978
15 14.05766201019287
Compare the times t
and u
, and it’s at least 3 times faster to do repeated multiplication instead of exponentiation, even with eval
.
Questions:

Why does this happen?

Should we intercept the case where a sparse matrix is being exponentiated to a positive integer and do repeated multiplication behind the scenes for performance?