Matrix vector multiplication performance

Question

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?

Answer 1

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.

Answer 2

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).

Answer 3

An addition to the other two (excellent) answers:

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 :)