Is there an "enhanced" numpy/scipy dot method?

Question

Problem

I would like to compute the following using numpy or scipy:

Y = A**T * Q * A

where A is a m x n matrix, A**T is the transpose of A and Q is an m x m diagonal matrix.

Since Q is a diagonal matrix I store only its diagonal elements as a vector.

Ways of solving for Y

Currently I can think of two ways of how to calculate Y:

1. Y = np.dot(np.dot(A.T, np.diag(Q)), A) and
2. Y = np.dot(A.T * Q, A).

Clearly option 2 is better than option 1 since no real matrix has to be created with diag(Q) (if this is what numpy really does...)
However, both methods suffer from the defect of having to allocate more memory than there really is necessary since A.T * Q and np.dot(A.T, np.diag(Q)) have to be stored along with A in order to calculate Y.

Question

Is there a method in numpy/scipy that would eliminate the unnecessary allocation of extra memory where you would only pass two matrices A and B (in my case B is A.T) and a weighting vector Q along with it?

Answer 1

(w/r/t the last sentence of the OP: i am not aware of such a numpy/scipy method but w/r/t the Question in the OP Title (i.e., improving NumPy dot performance) what's below should be of some help. In other words, my answer is directed to improving performance of most of the steps comprising your function for Y).

First, this should give you a noticeable boost over the vanilla NumPy dot method:

>>> from scipy.linalg import blas as FB >>> vx = FB.dgemm(alpha=1., a=v1, b=v2, trans_b=True) 

Note that the two arrays, v1, v2 are both in C_FORTRAN order

You can access the byte order of a NumPy array through an array's flags attribute like so:

>>> c = NP.ones((4, 3)) >>> c.flags       C_CONTIGUOUS : True          # refers to C-contiguous order       F_CONTIGUOUS : False         # fortran-contiguous       OWNDATA : True       MASKNA : False       OWNMASKNA : False       WRITEABLE : True       ALIGNED : True       UPDATEIFCOPY : False 

to change the order of one of the arrays so both are aligned, just call the NumPy array constructor, pass in the array and set the appropriate order flag to True

>>> c = NP.array(c, order="F")  >>> c.flags       C_CONTIGUOUS : False       F_CONTIGUOUS : True       OWNDATA : True       MASKNA : False       OWNMASKNA : False       WRITEABLE : True       ALIGNED : True       UPDATEIFCOPY : False 

You can further optimize by exploiting array-order alignment to reduce excess memory consumption caused by copying the original arrays.

But why are the arrays copied before being passed to dot?

The dot product relies on BLAS operations. These operations require arrays stored in C-contiguous order--it's this constraint that causes the arrays to be copied.

On the other hand, the transpose does not effect a copy, though unfortunately returns the result in Fortran order:

Therefore, to remove the performance bottleneck, you need to eliminate the predicate array-copying step; to do that just requires passing both arrays to dot in C-contiguous order*.

So to calculate dot(A.T., A) without making an extra copy:

>>> import scipy.linalg.blas as FB >>> vx = FB.dgemm(alpha=1.0, a=A.T, b=A.T, trans_b=True) 

In sum, the expression just above (along with the predicate import statement) can substitute for dot, to supply the same functionality but better performance

you can bind that expression to a function like so:

>>> super_dot = lambda v, w: FB.dgemm(alpha=1., a=v.T, b=w.T, trans_b=True) 

Answer 2

I just wanted to put that up on SO, but this pull request should be helpful and remove the need for a separate function for numpy.dot https://github.com/numpy/numpy/pull/2730 This should be available in numpy 1.7

In the meantime, I used the example above to write a function that can replace numpy dot, whatever the order of your arrays are, and make the right call to fblas.dgemm. http://pastebin.com/M8TfbURi

Hope this helps,

Answer 3

numpy.einsum is what you're looking for:

numpy.einsum('ij, i, ik -> jk', A, Q, A)

This shall not need any additional memory (though usually einsum works slowlier than BLAS operations)