What is the fastest way to quadratic form numpy array multiplication?

Question

I have tried those two alternatives

objective = lambda A, x : (np.dot(x.T ,np.dot(A, x)))[0,0]
objective = lambda A, x : (np.matrix(x).T * np.matrix(A) * np.matrix(x))[0,0]

With primary one I got 5 sec of running time with my algorithm With secondary I got 14 sec

With MATLAB I got 2 sec

I want to go with the Numpy but obviously I need a way to ameliorate this crummy results. How can I get faster quadratic form matrix, vector product?

Note: I profiled the code and this lambda function drinks the juice of all. Improvemnt: I merely remove the native Ubuntu package of scipy and numpy then installed with followings

sudo pip install numpy
sudo apt-get install libatlas-base-dev gfortran
sudo pip install scipy
sudo apt-get install libpng-dev libfreetype6-dev
sudo pip install matplotlib 

I improves the performance a bit, however still lower than Matlab

Answer 1

I have both NumPy and Matlab installed and they both take around 45 ms for a 10000x10000 matrix.

Considering your timings, I suspect that x is not a single column vector. If you want to do this computation for multiple column vectors all at once, look at my answer to this question: Calculate "v^T A v" for a matrix of vectors v . The timings you are listing are terribly slow if x is just a single column vector (in NumPy or Matlab).

I suspect, however, that the difference could also be coming from how your NumPy installation was compiled. This is really a timing of the BLAS functions used by NumPy and Matlab. I believe both are really calling the same underlying library on my machine since I have NumPy linked against Intel's MKL. If NumPy is built against a well-optimized BLAS like the Intel MKL, large vector operations like this should run at roughly the same speed as Matlab since they are both likely to be calling the same lower level BLAS functions. If your version of NumPy is not compiled using an optimized BLAS, performance will be worse.

If you know your installation of NumPy is already linked against the MKL, you can try setting the MKL_NUM_THREADS environment variable to match the number of processors on your system.

One simple way to get the properly compiled version of NumPy is to use a pre-built distribution. Anaconda and Enthought are very good but they will require a subscription to get the optimized versions. Academic licenses are available for free. You can also look here: http://www.lfd.uci.edu/~gohlke/pythonlibs/

Answer 2

Finally what I have done is to change the bounded library of numpy for linear algebra functions. It was using ATLAS for default but I push harder ( like 4 hours ) to change this to OpenBlas. I found that guide Compiling numpy with OpenBLAS integration and followed bit by bit. Result is working with faster time. It is still lacking in relation to Matlab (Intel MLK) 2.5 sec but tolerable with 3 sec execution.