why is A*v+B*v faster than (A+B)*v?

Question

$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?

Answer 1

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.

Answer 2

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.