33
$\begingroup$

$A$ and $B$ are $n \times n$ matrices and $v$ is a vector with $n$ elements. $Av$ has $\approx 2n^2$ flops and $A+B$ has $n^2$ flops. Following this logic, $(A+B)v$ should be faster than $Av+Bv$.

Yet, when I run the following code in matlab

A = rand(2000,2000);
B = rand(2000,2000);
v = rand(2000,1);
tic
D=zeros(size(A));
D = A;
for i =1:100
    D = A + B;
    (D)*v;
end
toc
tic
for i =1:100
    (A*v+B*v);
end
toc

The opposite is true. Av+Bv is over twice as fast. Any explanations?

$\endgroup$
  • $\begingroup$ Do you see the same ratio when adding maxNumCompThreads(1); on the first line? $\endgroup$ – user7440 Jan 15 at 22:27
  • $\begingroup$ @user7440 Yes, I got the exact same times. $\endgroup$ – Sam Christensen Jan 15 at 22:41
  • 3
    $\begingroup$ Could it be that assigning $D=A+B$ requires memory transfer which is usually very slow? $\endgroup$ – vibe Jan 15 at 22:55
  • 5
    $\begingroup$ When timing things, you should make sure that you're timing the same operation. In your first block, you also time setting D to zero and D=A. Also setting D= A+B, which you don't do in the second block. $\endgroup$ – Wolfgang Bangerth Jan 16 at 1:28
  • 2
    $\begingroup$ I timed just (A+B)*v - avoids setting to D and creating it. Still, the (Av+Bv) is 3 times faster for me. Difference is smaller for 10k points, but it is still over 2 times faster. $\endgroup$ – Zizy Archer Jan 16 at 7:48
63
$\begingroup$

Except for code which does a significant number of floating-point operations on data that are held in cache, most floating-point intensive code is performance limited by memory bandwidth and cache capacity rather than by flops.

$v$ and the products $Av$ and $Bv$ are all vectors of length 2000 (16K bytes in double precision), which will easily fit into a level 1 cache. The matrices $A$ and $B$ are 2000 by 2000 or about 32 megabytes in size. Your level 3 cache might be large enough to store one of these matrices if you've got a really good processor.

Computing $Av$ requires reading 32 megabytes (for $A$) in from memory, reading in 16K bytes (for $v$) storing intermediate results in the L1 cache and eventually writing 16K bytes out to memory. Multiplying $Bv$ takes the same amount of work. Adding the two intermediate results to get the final result requires a trivial amount of work. That's a total of roughly 64 megabytes of reads and an insignificant number of writes.

Computing $(A+B)$ requires reading 32 megabytes (for A) plus 32 megabytes (for B) from memory and writing 32 megabytes (for A+B) out. Then you have to do a single matrix-vector multiplication as above which involves reading 32 megabytes from memory (if you've got a big L3 cache, then perhaps this 32 megabytes is in that L3 cache.) That's a total of 96 megabytes of reads and 32 megabytes of writes.

Thus there's twice as much memory traffic involved in computing this as $(A+B)v$ instead of $Av+Bv$.

Note that if you have to do many of these multiplications with different vectors $v$ but the same $A$ and $B$, then it will become more efficient to compute $A+B$ once and reuse that matrix for the matrix-vector multiplications.

| cite | improve this answer | |
$\endgroup$
20
$\begingroup$

Your code is limited by memory bandwidth. For trivial math, it's often better to count memory accesses rather than flops. You'll get the following table:

operation             memory reads/writes

matrix + matrix       3n²
matrix * vector       2n²+n  (if vector is not cached)
matrix * vector       n²+2n  (if vector is only read once)
vector + vector       3n

(A+B)*v (non-cached)  5n²+n       
A*v+B*v (non-cached)  4n²+5n

(A+B)*v (cached)      4n²+2n
A*v+B*v (cached)      2n²+7n

If the vector is small enough to fit into cache, Av+Bv is indeed twice as fast. Even without a cache, Av+Bv is faster, albeit not by a factor of two.

Also note that you already need 2n²+2n memory accesses just to read the arguments from memory and write back the result, so the cached variant of Av+Bv is close to optimal.

| cite | improve this answer | |
$\endgroup$
  • 2
    $\begingroup$ Note that unless special CPU instructions for streaming-stores are used, often each "write" should actually be counted as a read AND a write. This is because the CPU has to load the output into cache before writing to it. $\endgroup$ – Simon Jan 16 at 19:55
  • 1
    $\begingroup$ @Simon - Note that many CPUs have streaming write detection that avoids this behavior on sequential writes. $\endgroup$ – TLW Jan 17 at 3:31
  • 1
    $\begingroup$ @TLW Sounds interesting, I have not come across that, and fail to find relevant information with google. Can you provide a reference? $\endgroup$ – Simon Jan 17 at 15:46
  • 1
    $\begingroup$ @Simon - See e.g. infocenter.arm.com/help/index.jsp?topic=/com.arm.doc.ddi0500e/… "read allocate mode". It's a best-effort optimization - essentially if the processor detects that too many full cache lines have been written before the associated reads complete, it stops bothering to try to fetch the reads. It's an important optimization, as otherwise sequential writes are effectively at half performance (or worse). $\endgroup$ – TLW Jan 18 at 20:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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