# Matrix vector multiplication performance

I have been learning about the impact of cache size on code performance. I wrote a small code to see how using a column major loop in MATLAB would be better than using a row major loop, since MATLAB stores matrices in column major like FORTRAN. I also compared against MATLAB's internal multiplication routine. Here is the code and the results:

% row major access
tic
for i=1:n
for j=1:n
b1(i)=b1(i)+A(i,j)*x(j);
end
end
t1(count)=toc;

% column major access
tic
for j=1:n
for i=1:n
b2(i)=b2(i)+A(i,j)*x(j);
end
end
t2(count)=toc;

% column major vector ops
tic
for i=1:n
b3(i)=b3(i)+A(i,:)*x(:);
end
t3(count)=toc;

% row major vector ops
for j=1:n
b4(:)=b4(:)+A(:,j)*x(j);
end
t4(count)=toc;

% MATLAB built in
tic
b5=A*x;
t5(count)=toc;

% double vop
tic
b6=A(:,:)*x(:);
t6(count)=toc;


The column major loop is faster than the row major loop, but how is MATLAB so much faster?

MATLAB loops are slow. This is known. If you use something like C++ or Julia you'll get much closer (like 5x-10x or so IIRC). So that's why it's super crazy.

Still, the algorithms that MATLAB is calling are actually BLAS libraries. Specifically, MATLAB is calling Intel MKL under the hood here, so it's a blazing fast multithreaded code. In fact, it does all sorts of things. It first loads all inputs into contiguous stack-allocated buffers Then it does some bit twiddling on pointers to make things move faster, SIMD's a few things, and does this all with parallelism. It's wild stuff to try and understand. If you're curious OpenBLAS is open-source (instead of MKL which is not), so you can take a look and try to understand what's going on, but it's crazy!

## Edit:

I originally said "and then does operations via stencils so that way it's actually an order $\mathcal{O}(n^{log_2 7})$ instead of the naive $\mathcal{O}(n^{3})$ that using the definition of matrix multiplication gives". I guess it doesn't actually use Strassen in practice. Also, that's for matrix multiplication, not specifically matrix-vector multiplication like the question asked. Sorry.

• Pretty amazing stuff!! – EternusVia Jan 18 '18 at 1:12
• MKL uses Strassen? I had never heard about it. Actually, in this recent discussion it seems that an Intel employee claims the opposite. – Federico Poloni Jan 18 '18 at 22:50
• Strassen isn't useful for matrix-vector multiplication (it's a matrix-matrix multiplication algorithm), and matrix-vector takes only $O(n^{2})$ time, so those parts of this answer aren't correct. – Brian Borchers Jan 19 '18 at 2:23
• MATLAB does use Intel MKL on Intel processors, and it is vastly faster than MATLAB loops. – Brian Borchers Jan 19 '18 at 2:32

Matlab is probably handing off to a vendor BLAS (MKL, etc). Those libraries are totally different animals than a Matlab double-for loop (thread-parallel, SIMD intrinsics, algorithmically reblocked, the works).

Also, for loops are generally pokey in Matlab anyway. You'll probably see speedups (in both cases) just by rewriting the calculation in some form of vector notation (handle whole columns/rows by writing A(:,j)/A(i,:), that kind of thing).

• Good point! I'll try the vector notation variants as well. – EternusVia Jan 17 '18 at 22:57

Understanding the source code of something like OpenBLAS can be a daunting task. As an alternative starting point, you could read through

How to Write Fast Numerical Code: A Small Introduction by Chellappa et al.

Section 5 deals with the MMM (matrix-matrix multiplication) operation. You start from a "naïve" triple-loop MMM kernel and examine its memory data access patterns. Then the authors show you how to:

1. optimize cache memory usage via blocking and buffering (section 5.1)
2. optimize MMM at the CPU and register level using blocking and loop unrolling (section 5.2)
3. "tune" the performance by selecting the optimal parameters (cache/register block size, etc.) for your CPU architecture (section 5.3). You do this with a tool like ATLAS (section 5.4).

This won't get you anywhere near the performance of a highly-optimized library like Intel MKL, but you should nevertheless see a significant improvement over the naïve kernel. And hey, it's a fun exercise :)

• I think another good resource for the adventurous implementor is the BLIS documentation, see cs.utexas.edu/users/flame/pubs/blis1_toms_rev3.pdf (it's more advanced than the Chellappa citation, but more pedagogical perhaps than jumping straight to source). I think the key points are (i) implement all BLAS3 in terms of a "fat dot product" microkernel that maps onto platform-specific SIMD instructions (ii) judiciously copy from input buffers to aligned and SIMD-reshaped workspace, blocking and reordering loop nests to maximize potential to reuse that workspace. – rchilton1980 Jan 19 '18 at 4:09