Re: [Numpy-discussion] Dot product threading?

[Angus M. <am...@gm...>]

Pau Gargallo wrote:
> I think that the operation Angus wants is the following (at least I
> would like that one ;-)
> 
[snip]
> 
> But sometimes it would be useful (at least for me) to have:
>     a.shape = I+(n,k)
>     b.shape = I+(k,m)
> and to get only:
>     dot2( a,b ).shape == I+(n,m)
>     dot2( a,b )[i] == dot2( a[i], b[i] )
> 
> This would be a natural extension of the scalar product (a*b)[i] ==
> a[i]*b[i]
> If dot2 was a kind of ufunc, this will be the expected behaviour,
> while the current dot's behaviour will be obtained by dot2.outer(a,b).
> 
> Does this make any sense?

Thank-you Pau, this looks exactly like what I want. For further
explanation, here, I believe, is an implementation of the desired
routine using a loop. It would, however, be great to do this using
quicker (ufunc?) machinery. Pau, can you confirm that this is the same
as the routine you're interested in?

def dot2(a,b):
    '''Returns dot product of last two dimensions of two 3-D arrays,
    threaded over first dimension.'''
    try:
        assert a.shape[1] == b.shape[2]
        assert a.shape[0] == b.shape[0]
    except AssertionError:
        print "incorrect input shapes"
    res = zeros( (a.shape[0], a.shape[1], a.shape[1]), dtype=float )
    for i in range(a.shape[0]):
        res[i,...] = dot( a[i,...], b[i,...] )
    return res

I think the 'arrays of 2-D matrices' comment (which I've snipped out, oh
well) captures the idea well.

Angus.
-- 
Angus McMorland
email a.m...@au...
mobile +64-21-155-4906

PhD Student, Neurophysiology / Multiphoton & Confocal Imaging
Physiology, University of Auckland
phone +64-9-3737-599 x89707

Armourer, Auckland University Fencing
Secretary, Fencing North Inc.

-- 
Angus McMorland
email a.m...@au...
mobile +64-21-155-4906

PhD Student, Neurophysiology / Multiphoton & Confocal Imaging
Physiology, University of Auckland
phone +64-9-3737-599 x89707

Armourer, Auckland University Fencing
Secretary, Fencing North Inc.