Improving performance in matrix operations

[Oscar Benjamin]

On 9 March 2016 at 20:09, Drimades <e.zhupa at gmail.com> wrote:
> I'm doing some tests with operations on numpy matrices in Python. As an example, it takes about 3000 seconds to compute eigenvalues and eigenvectors using scipy.linalg.eig(a) for a matrix 6000x6000. Is it an acceptable time?

I don't know really but you need to understand that numpy delegates
this kind of operation to the underlying BLAS library. It's possible
to have different BLAS libraries depending on how you installed numpy.
For example if you install numpy from
    http://www.lfd.uci.edu/~gohlke/pythonlibs/
then you will have a numpy that is linked with the Intel MKL library
for BLAS which I think is that same as used in e.g. Matlab and many
other things. Alternatively if you installed from the numpy
sourceforge page then you'll have the ATLAS BLAS library. If you're
using e.g. Ubuntu and installed numpy from the Ubuntu repos it's
possible that you're using numpy's vendored unoptimised BLAS library.

Each of these different BLAS libraries has different characteristics
in terms of accuracy and speed so it's worth knowing which one you're
actually using.

> Any suggestions to improve? Does C++ perform better with matrices?

If you were working in C++ you would still want to link to a BLAS
library to do this so I don't see why it would make any difference
except that it would require you to work out how to compile and use
BLAS directly and then link to it from your C++ code.

> Another thing to consider is that matrices I'm processing are heavily sparse.

Then you should definitely use something that is targeted at sparse
matrices (as suggested by Fabien). This can give a massive boost in
performance.

> Do they implement any parallelism? While my code is running, one of my cores is 100% busy, the other one 30% busy.

It sounds like the particular BLAS library you are using is not using
several cores for this workload. Different BLAS libraries have
different capabilities. Again you need to figure out which one you've
got and how it's compiled. It's possible that e.g. MKL has a parallel
eig function but that it is compiled with that behaviour disabled in
your setup.

--
Oscar