3
$\begingroup$

I have to calculate multiple dot products between a vector and consecutive part of another (longer) vector, shifted at certain number of elements. Naive implementation takes about 0.1 seconds.

l = 100001;
l2 = 5001;
rvec = RandomReal[{-1, 1}, l];
m = RandomReal[{-1, 1}, 2 l - 1];
ans1 = ParallelTable[rvec.m[[l - 20 (-1 + i) ;; 2 l - 20 (-1 + i) - 1]], {i, 1, l2}]; // RepeatedTiming

If I do the same calculation, but shifted vectors are gathered in matrix, and matrix-vector product is calculated, it takes 0.035 seconds, but 5000 times more memory.

M = Table[RandomReal[{-1, 1}, l], l2];(*careful,7gb matrix*)

ans2 = M.rvec; // RepeatedTiming

Is it possible to speed up the first method?

$\endgroup$
1
  • $\begingroup$ Sounds a bit like a job for ListCorrelate, but I have not figured out yet, how to use it in this case. $\endgroup$ Apr 3, 2023 at 4:23

2 Answers 2

6
$\begingroup$

Finally, I have been able to do it with ListCorrelate. The idea is to use Partition[...,c] in conjuctions with a 2-dimensional kernel ker to convince ListCorrelate to do the dot-products in steps of c = 20. Alas, the problems is that the lengths of rvec and m are not divisable by c. Thus Partition simply drops some of their last elements. But we can insert some padding to prevent that. The result is an implementation that is 15 times as fast.

l = 100001;
c = 20;
l2 = Quotient[l, c] + 1;
rvec = RandomReal[{-1, 1}, l];
m = RandomReal[{-1, 1}, 2 l - 1];

ans1 = ParallelTable[rvec . m[[l - c (i - 1) ;; 2 l - c (i - 1) - 1]], {i, 1, l2}]; // RepeatedTiming // First

ans2 = With[{
  ker = Partition[PadRight[rvec, c l2, 0.], c],
  data = Partition[PadRight[m, (2 (l2) - 1) c, 0.], c]
  },
 Reverse[Flatten[ListCorrelate[ker, data]]]
 ]; // RepeatedTiming // First

Max[Abs[ans1 - ans2]]/Max[Abs[ans1]]

0.113728

0.0074245

3.30665*10^-15

Remark:

At least on my machine, ParallelTable[..., Method -> "CoarsestGrained"] made it twice as fast.

$\endgroup$
5
  • 3
    $\begingroup$ This is good (I upvoted) but maybe you worked too hard in trying to, well, not work too hard. If you do a full ListCorrelate you avoid all partitioning and the like. Then just walk the result in steps of 20 to get the values desired, like so. In[225]:= ans1 = Table[ rvec . m[[l - 20 (-1 + i) ;; 2 l - 20 (-1 + i) - 1]], {i, 1, l2}]; // Timing Out[225]= {0.467425, Null} In[226]:= Timing[ans2 = Reverse@ListCorrelate[rvec, m];] Max[Abs[ans2[[1 ;; -1 ;; 20]] - ans1]] Out[226]= {0.014101, Null} Out[227]= 7.67386*10^-13 Speed factor is around 25-30. $\endgroup$ Apr 3, 2023 at 14:35
  • $\begingroup$ I see. Hm. I was somewhat reluctant to discard 95% of the flops... And on my M1 ListCorrelate[rvec, m] takes as long as "my" method. $\endgroup$ Apr 3, 2023 at 15:35
  • 2
    $\begingroup$ Unless there's a nice way to condense the vectors that go through FFTs under the hood, you would not be discarding 95% of the flops, just 95% of the result. $\endgroup$ Apr 3, 2023 at 23:11
  • 1
    $\begingroup$ Ah, I see. ListCorrelate does use FFT here? I tried to code that explicitly after your first post, but already Fourier[m] took 3 times longer than ListCorrelate[rvec, m]. So I assumed that ListCorrelate might recurse to something more straightforward for this data size. Anyways, I might be wrong, but wouldn't the partitioning into a matrix of size $n \times 20$ allow us to use 20 FFTs in parallel? Looking on the CPU stats, ListCorrelate does not appear to be parallelized--but that could be due to the small data size. $\endgroup$ Apr 4, 2023 at 5:02
  • $\begingroup$ Well, I thought it used a FT for this size. But I didn’t check and maybe I’m wrong (say it ain’t so!). As for some form of parallelization, I doubt this would happen explicitly in that code. So it will be fully serial unless multithreaded MKL (or equivalent) code is in use. That at least is my guess. $\endgroup$ Apr 4, 2023 at 13:48
0
$\begingroup$

I tried the normal Table construct, takes half the time as ParallelTable. The expressiion is a a running ListConvolve, but its tedious to find the exact form. So I took the obvious forms of Array contructs

In[1]:=

ClearAll["Global`*"]

In[2]:= l = 100001;
l2 = 5001;
rvec = RandomReal[{-1, 1}, l];
m = RandomReal[{-1, 1}, 2 l - 1];

In[6]:= Quiet[ ans[0] = Table[rvec . m[[l - 20 (-1 + i) ;; 2 l - 20 (-1 + i) - 1]], {i, 1, l2}]; // RepeatedTiming]

Out[6]= {0.375008, Null}

In[7]:= Quiet[ans[1] =  ParallelTable[ rvec . m[[l - 20 (-1 + i) ;; 2 l - 20 (-1 + i) - 1]], {i, 1, l2}]; // RepeatedTiming]

Out[7]= {0.595901, Null}

Now the functional kernel of the tables: just replace

{ i -> # , {i,1,l2}->  {l2}} 

In[8]:= ff = ( rvec . Take[m, {l - 20 (# - 1), 2 l - 1 - 20 (# - 1) }] &);

ParallelArray:

In[9]:=   RepeatedTiming[ans[2] = ParallelArray[ff[#] &, {l2}];]

Out[9]= {0.466023, Null}

Parallelizing Array:

In[10]:= RepeatedTiming[Parallelize[ ans[3] = Array[ ff, {l2}]];]

Out[10]= {0.492969, Null}

Parallelizing mapping to Range

In[11]:= RepeatedTiming[Parallelize[ ans[4] =  ff /@ Range[l2]];]

Out[11]= {0.817287, Null}

Just mapping is the fastest :-)

In[12]:= RepeatedTiming[ ans[5] =  ff /@ Range[l2];]

Out[12]= {0.158288, Null}

These is done on a I5 10xxx notebook.

As is widely known now, notebook systems never utilize the full power of muliticore or GPU processing capabilities above some 25% of the theoretical timings because thermal management has absolute priority.

$\endgroup$
4
  • $\begingroup$ I don't get the statement about "notebook systems" and "thermal management ", and I think it is just wrong. The speed difference is more about interpreted vs. compiled language -- or rather whether the interpreted language has already a good compiled and parallelized implementation of an algorithm suitable for the problem. And about runtime security checks (e.g. against segmentation faults) in code that is linked to the notebook interface. Moreover, the desire for immutable data types is also a (minor) cause of slowdown as some superfluous copies may have to be made. $\endgroup$ Apr 3, 2023 at 18:01
  • $\begingroup$ As they say in the specifacations of any CPU, the maximal speed of around 3.5 GHz is switched off to very low frequencies if the core temperature reaches 70C. Often notebooks stall completely because even the system bus system is pu on Hold. For the same procedures on a gaming desktop in these cases the cooling system switches to maximum speed and sometimes it is working $\endgroup$
    – Roland F
    Apr 4, 2023 at 8:35
  • $\begingroup$ Ah, with "notebook system" you meant a portable computer, not interpreted languages with a notebook interface like Mathematica, Python, or Matlab. My fault. Anyways, that must be a very, very old laptop for which thermal throttling would be so relevant for this tiny computational problem at hand. Also be assured that some non-Intel notebooks have significantly less problems in holding up to their marketed performance promises... ;) $\endgroup$ Apr 4, 2023 at 9:13
  • $\begingroup$ I do Mathematica with physics problems always on the newest and most expensive nootebooks since 1995. No one solved any of the problems, that took really time and memory. Put your hand on the cooling fan output and compare with a desktop system. $\endgroup$
    – Roland F
    Apr 5, 2023 at 11:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.