numpy : efficient sum computations

[TG]

Hi there.

I want to do some intensive computations with numpy, and I'm
struggling a bit to find myyyyy wayyyyyy. Here is the problem :

m and d are two matrices :

m.shape = (x,y,a,b)
d.shape = (a,b)

I want to return

i.shape = (x,y) with
i[x,y] = sum(m[x,y] * d)

I already found that

m[:,:] * d

will give me a matrix of shape (x,y,a,b) containing the products.

Now I want to sum up on axis 2 and 3. If I do :

(m[:,:] * d).sum(axis=3).sum(axis=2)

it seems like I get my result.

I'm wondering : is this syntax leading to efficient computation, or is
there something better ?

thanks

Thomas

 

[TG]

Okay, another one which I don't have answer for.

it is the reverse case, sort of :

phi.shape (x,y)
d.shape (a,b)

I want to return m :

m.shape = (x,y,a,b)

with
m[x,y] = d * phi[x,y]

currently, my code is :

m = empty(phi.shape + d.shape)
m[:,:] = d

this repeats the matrix d x*y times, which is what I want. But now, if
I do :

m[:,:] = d * phi[:,:]

in order to multiply each "d" by a specific value of "phi", it doesn't
work.

<type 'exceptions.ValueError'>: shape mismatch: objects cannot be
broadcast to a single shape

Right now I'm stuck here. If anyone can help, I would be glad.

 

[Robert Kern]

Hi there.

I want to do some intensive computations with numpy, and I'm
struggling a bit to find myyyyy wayyyyyy.

The best place to ask numpy questions is on the numpy mailing list.

http://www.scipy.org/Mailing_Lists

Here is the problem :

m and d are two matrices :

m.shape = (x,y,a,b)
d.shape = (a,b)

I want to return

i.shape = (x,y) with
i[x,y] = sum(m[x,y] * d)

I already found that

m[:,:] * d

will give me a matrix of shape (x,y,a,b) containing the products.

Now I want to sum up on axis 2 and 3. If I do :

(m[:,:] * d).sum(axis=3).sum(axis=2)

it seems like I get my result.

I'm wondering : is this syntax leading to efficient computation, or is
there something better ?

That's about it. There's no need for the [:,:], though.

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
that is made terrible by our own mad attempt to interpret it as though it had
an underlying truth."
-- Umberto Eco

 

[Robert Kern]

Okay, another one which I don't have answer for.

it is the reverse case, sort of :

phi.shape (x,y)
d.shape (a,b)

I want to return m :

m.shape = (x,y,a,b)

with
m[x,y] = d * phi[x,y]

currently, my code is :

m = empty(phi.shape + d.shape)
m[:,:] = d

this repeats the matrix d x*y times, which is what I want. But now, if
I do :

m[:,:] = d * phi[:,:]

in order to multiply each "d" by a specific value of "phi", it doesn't
work.

<type 'exceptions.ValueError'>: shape mismatch: objects cannot be
broadcast to a single shape

Right now I'm stuck here. If anyone can help, I would be glad.

Use broadcasting to get the desired behavior. It's a widely used feature of
numpy that takes some getting used to, but it is quite powerful. Here is some
documentation:

http://www.scipy.org/EricsBroadcastingDoc

It still references Numeric, numpy's predecessor, but the concepts are still
valid. Use the numpy.newaxis object to add new axes to phi and d in the right
places, then multiply them together:


In [8]: from numpy import arange, newaxis

In [9]: phi = arange(5*4).reshape((5,4))

In [10]: phi.shape
Out[10]: (5, 4)

In [11]: d = arange(3*2).reshape((3,2))

In [12]: d.shape
Out[12]: (3, 2)

In [13]: m = phi[:,:,newaxis,newaxis] * d[newaxis,newaxis,:,:]

In [14]: m.shape
Out[14]: (5, 4, 3, 2)

In [17]: m[1, 2]
Out[17]:
array([[ 0, 6],
[12, 18],
[24, 30]])

In [18]: d * phi[1, 2]
Out[18]:
array([[ 0, 6],
[12, 18],
[24, 30]])

In [19]: for i in range(phi.shape[0]):
....: for j in range(phi.shape[1]):
....: assert (m[i,j] == d*phi[i,j]).all()
....:
....:

In [20]:

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
that is made terrible by our own mad attempt to interpret it as though it had
an underlying truth."
-- Umberto Eco