8
$\begingroup$

I'm looking at speeding up matrix-vector products but everything I read is about how to do it for very large matrices. My case, the matrices are small but the number of times it must be done is very large.

What methods, if any, are there to optimize this? Would it be faster to construct a really big diagonal block matrix out of the small matrices and one big vector made of the smaller vectors and use the techniques for the large matrix-vector speedups? Or would setting up the global matrix and vector kill any benefit there?

Improve this question
$\endgroup$
4
  • $\begingroup$ Do you have to multiply the same matrix times many vectors or is there no reuse of the matrices? $\endgroup$
    – Brian Borchers
    Commented Jan 11, 2014 at 5:20
  • $\begingroup$ I doubt that using the techniques for large matrices for a large block-diagonal will give you a significant speedup. How small is 'small' in your case and how many matrices are we typically talking about? Do you know anything else about these matrices, e.g., do they describe rotations etc.? $\endgroup$
    – Christian Waluga
    Commented Jan 11, 2014 at 6:36
  • $\begingroup$ @BrianBorchers There's no re-using the matrix, it's different for every point every time step $\endgroup$
    – tpg2114
    Commented Jan 11, 2014 at 12:16
  • $\begingroup$ @ChristianWaluga They are 5x5 up to maybe 10x10 sometimes, dense, not symmetric and not diagonally dominant generally speaking. How many times it needs to be done varies on the case, but typically 10000 to 60000 times per time step $\endgroup$
    – tpg2114
    Commented Jan 11, 2014 at 12:17

3 Answers 3

Reset to default
7
$\begingroup$

Before trying to optimize your code, it is worth asking if there is anything to optimize to begin with. Libraries that optimize matrix-vector products do so by working around two issues: limitations on the size of the cache and the latency to load data from memory. The first is done by using the data that's currently in cache to its fullest extent for everything it needs to be used for before replacing it by other data, the latter is done by prefetching data into the cache before actually making use of it.

In your case, you have relatively little arithmetic intensity of your data -- you load data from memory, you use it exactly once, and then you move on to the next matrix. This leaves only the second avenue to optimize for: prefetch data before you use it.

But, as I said, before trying to optimize things, it may be worth figuring out what you already have: time how many matrix-vector products you are doing per second, calculate how many bytes this requires to load from memory onto your processor, and then compare this with the bandwidth of the processor you happen to have in your machine. You may find that there is nothing you can do to make things faster.

Improve this answer
$\endgroup$
5
$\begingroup$

Generating C++ code and use Eigen/Armadillo is a possibility, but this depends on the application.

The solution for us was to just write out the result explicitly for $N<8$. Without the loops, the code is very fast with modern compilers and vector support (sse2, avx2, and avx512 in 64 bit).

Do take care of memory alignment of your data (up to 64 bytes aligned) and do restrict the pointers to make life easier for the compiler. No need to use multicore support with these matrix sizes, the overhead is larger than the gain.

We use scripting to automatically generate separate functions for each possible combination, and cache the function pointers for consecutive calls.

There is a nice cheap library for small matrix operations that is highly optimized by hand, called OptiVec, and it works very well in our case. We use it for $N \geq 8$.

Improve this answer
$\endgroup$
4
$\begingroup$

It may not actually matter since your matrices are cache contained already, but you should be calling dgemv() or sgemv() or the equivalent from the best BLAS library you can get your hands on. You should try the Intel MKL if you can get access to it, and also BLIS or ATLAS or one of the many other optimized BLAS libraries out there.

Improve this answer
$\endgroup$
7
  • 1
    $\begingroup$ Interestingly the BLAS routines run slower than the MATMUL intrinsic function in Fortran, even with MKL and hardware-specific implementations. $\endgroup$
    – tpg2114
    Commented Jan 12, 2014 at 6:04
  • 1
    $\begingroup$ I'm not entirely surprised by that, but I'd have to see the code to know for sure. There are a number of issues to worry about. I'd would have to suggest checking the alignment of your arrays before writing off the MKL, but at these small sizes, the MKL MATVEC may not be higly optimized. $\endgroup$
    – Bill Barth
    Commented Jan 12, 2014 at 14:52
  • $\begingroup$ Bill, a co-worker of mine ran into this same problem. The conclusion was that either there was some non-negligible overhead in the MKL call or it was otherwise not well optimized for small matrices. Either way, a hand-written matmul was considerably faster when doing a very large number of 5x5 matrix multiplications. $\endgroup$
    – Aurelius
    Commented Jan 14, 2014 at 2:38
  • 2
    $\begingroup$ The matmul for a NxN multiplying an Nx1 should have $O(N^2)$ operations. How did you get to $O(N^3)$? $\endgroup$
    – Bill Barth
    Commented Jan 14, 2014 at 22:38
  • 3
    $\begingroup$ If the matrices are very small (e.g. 4x4), try giving one of the templated libraries a go - it may remove a lot of function call overhead. Eigen is a good candidate. $\endgroup$
    – Damien
    Commented Jan 19, 2014 at 4:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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