# Why can Mathematica compute numerical sums more efficiently when they are written as matrix operations?

Let $f(n)$ and $K(n,m)$ be functions such that the double sum, which we wish to evaluate numerically, $$\sum_{n=1}^a \sum_{m=1}^a f(n) f(m) K(n,m)$$ exists when $a$ is some large positive number. I have noticed that Mathematica evaluates these sums much more quickly when they are expressed using matrix multiplication than when they are done the 'naive' way.

For example,

a=1000;
f = Table[RandomReal[], {n, 1, a}];
k = Table[RandomReal[], {n, 1, a}, {m, 1, a}];

ParallelSum[f[[n]] f[[m]] k[[n, m]], {n, 1, a}, {m, 1, a}] // AbsoluteTiming
{5.73706, 121605.}

(({f}.k).Transpose[{f}])[[1, 1]] // AbsoluteTiming
{0.0392919, 121605.}


This is not news to many of you; I found a few previous questions where askers were advised to convert sums to an equivalent expression using matrix multiplication (here and here) precisely because it is so much more efficient.

I'm curious as to what's going on 'under the hood.' Why is Mathematica so much faster at evaluating equivalent expressions using matrix algebra rather than Sum[]? Is Mathematica using a well-known algorithm that somehow does better than just summing the terms? Or is it that the Sum[] function itself, or calling the elements from the arrays are the culprit?

• You could, BTW, simply type f.k.f and get the same. MMA figures out orientation of the vectors fine on its own. – LLlAMnYP Apr 29 '16 at 11:53
• Matrix multiplication using Dot will call BLAS GEMM routines, which use optimized machine instructions and are tuned for memory bandwidth utilization efficiency. A manual loop over the elements at the top level, especially using ParallelSum (which necessitates copying the operands to the subkernels, performing the sum, and then copying the result back to the main kernel), cannot hope to compete with that. You can probably do better with Compile, but you still won't match BLAS. – Oleksandr R. Apr 29 '16 at 12:18