Linear combinations in python/numpy

Question

Lets say I have 3 Numpy arrays, A1,A2,A3, and 3 floats, c1,c2,c3

and I'd like to evaluate B = A1*c1 + A2*c2+ A3*c3.

Will Numpy compute it like the below example?

 E1 = A1*c1
 E2 = A2*c2
 E3 = A3*c3
 D1 = E1+E2
 B = D1+E3

Or is it more clever than that? In c++ I had a neat way to abstract this kind of operation.

I defined a series of general 'LC' template functions, LC for linear combination like:

template<class T,class D>
void LC( T & R,
    T & L0,D C0,
    T & L1,D C1,
    T & L2,D C2)
{
    R = L0*C0
        +L1*C1
        +L2*C2;        
}

and then specialized this for various types,

so for instance, for an array the code looked like

for (int i=0; i<L0.length; i++)
    R.array[i] =
    L0.array[i]*C0 +
    L1.array[i]*C1 +
    L2.array[i]*C2;

thus avoiding having to create new intermediate arrays.

This may look messy but it worked really well.

I could do something similar in python, but I'm not sure if it's necessary.

Untidy

Answer 1

While numpy, in theory, could at any time always upgrade its internals to perform wondrous optimizations, at the present time it does not: B = A1*c1 + A2*c2 + A3*c3 will indeed produce and then discard intermediate temporary arrays ("spending" some auxiliary memory, of course -- nothing else).

B = A1 * c1 followed by B += A2 * c2; B += A3 * c3, again at this time, will therefore avoid spending some of that temporary memory.

Of course, you'll be able to tell the difference only if you're operating in an environment with scarce real memory (where some of that auxiliary memory is just virtual and leads to page faults) and for sufficiently large arrays to "spend" all real memory and then some. Under such extreme conditions, however, a little refactoring can buy you some performance.

Hardunn

Answer 2

This is the idea behind numexpr (A Fast numerical array expression evaluator for Python and NumPy). You might give this package a try before compiling your own routines.

Kandace

Answer 3

Some years after the question was posed, I was searching the same functionality.

The following answer was given to me including the use of np.vecdot:

A = np.array([A1, A2, A3])
c = np.array([c1, c2, c3])
B = np.vecdot(A, c, axis=0)

I post it to share.

Pudens