numpy-discussion

[Numpy-discussion] Dot product threading?

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

Is there a way to specify which dimensions I want dot to work over?

For example, if I have two arrays:

    In [78]:ma = array([4,5,6]) # shape = (1,3)
    In [79]:mb = ma.transpose()   # shape = (3,1)
    In [80]:dot(mb,ma)
    Out[80]:
    array([[16, 20, 24],
           [20, 25, 30],
           [24, 30, 36]])

No problems there.
Now I want to do that multiple times, threading over the first dimension
of larger arrays:

    In [85]:mc = array([[[4,5,6]],[[7,8,9]]])        # shape = (2, 1, 3)
    In [86]:md = array([[[4],[5],[6]],[[7],[8],[9]]]) #shape = (2, 3, 1)

and I want to calculate the two (3, 1) x (1, 3) dot products to get a
shape = (2, 3, 3) result, so that

    res[i,...] == dot(md[i,...], mc[i,...])

From my example above res[0,...] would be the same as dot(mb,ma) and
res[1,...] would be

    In [108]:dot(md[1,...],mc[1,...])
    Out[108]:
    array([[49, 56, 63],
           [56, 64, 72],
           [63, 72, 81]])

I could do it by explicitly looping over the first dimension (as
suggested by my generic example), but it seems like there should be a
better way by specifying to the dot function the dimensions over which
it should be 'dotting'.

Cheers,

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.

[Numpy-discussion] Re: Dot product threading?

[Robert K. <rob...@gm...>]

Angus McMorland wrote:
> Is there a way to specify which dimensions I want dot to work over?

Use swapaxes() on the arrays to put the desired axes in the right places.

In [2]: numpy.swapaxes?
Type:           function
Base Class:     <type 'function'>
String Form:    <function swapaxes at 0x5d6d30>
Namespace:      Interactive
File:
/Library/Frameworks/Python.framework/Versions/2.4/lib/python2.4/site-packages/numpy-0.9.7.2476-py2.4-macosx-10.4-ppc.egg/numpy/core/oldnumeric.py
Definition:     numpy.swapaxes(a, axis1, axis2)
Docstring:
    swapaxes(a, axis1, axis2) returns array a with axis1 and axis2
    interchanged.


In [3]: numpy.dot?
Type:           builtin_function_or_method
Base Class:     <type 'builtin_function_or_method'>
String Form:    <built-in function dot>
Namespace:      Interactive
Docstring:
    matrixproduct(a,b)
    Returns the dot product of a and b for arrays of floating point types.
    Like the generic numpy equivalent the product sum is over
    the last dimension of a and the second-to-last dimension of b.
    NB: The first argument is not conjugated.

-- 
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

Re: [Numpy-discussion] Dot product threading?

[Angus M. <a.m...@au...>]

Robert Kern wrote:
> Angus McMorland wrote:
>
>>Is there a way to specify which dimensions I want dot to work over?
>
> Use swapaxes() on the arrays to put the desired axes in the right places.

Thanks for your reply, Robert. I've explored a bit further, and have
made sense of what's going on, to some extent, but have further questions.

My interpretation of the dot docstring, is that the shapes I need are:
    a.shape == (2,3,1) and b.shape == (2,1,3)
so that the sum is over the 1s, giving result.shape == (2,3,3)
but:
    In [85]:ma = array([[[4],[5],[6]],[[7],[8],[9]]]) #shape = (2, 3, 1)
    In [86]:mb = array([[[4,5,6]],[[7,8,9]]])         #shape = (2, 1, 3)
so
    In [87]:res = dot(ma,mb).shape
    In [88]:res.shape
    Out[88]:(2, 3, 2, 3)
such that
    res[i,:,j,:] == dot(ma[i,:,:], mb[j,:,:])
which means that I can take the results I want out of res by slicing
(somehow) res[0,:,0,:] and res[1,:,1,:] out.

Is there an easier way, which would make dot only calculate the dot
products for the cases where i==j (which is what I call threading over
the first dimension)?

Since the docstring makes no mention of what happens over other
dimensions, should that be added, or is this the conventional numpy
behaviour that I need to get used to?

Cheers,

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.

[Numpy-discussion] Re: Dot product threading?

[Robert K. <rob...@gm...>]

Angus McMorland wrote:
> Robert Kern wrote:
> 
>>Angus McMorland wrote:
>>
>>>Is there a way to specify which dimensions I want dot to work over?
>>
>>Use swapaxes() on the arrays to put the desired axes in the right places.
> 
> Thanks for your reply, Robert. I've explored a bit further, and have
> made sense of what's going on, to some extent, but have further questions.
> 
> My interpretation of the dot docstring, is that the shapes I need are:
>     a.shape == (2,3,1) and b.shape == (2,1,3)
> so that the sum is over the 1s, giving result.shape == (2,3,3)

I'm not sure why you think you should get that resulting shape. Yes, it will
"sum" over the 1s (in this case there is only one element in those axes so there
is nothing really to sum). What exactly are the semantics of the operation that
you want? I can't tell from just the input and output shapes.

> but:
>     In [85]:ma = array([[[4],[5],[6]],[[7],[8],[9]]]) #shape = (2, 3, 1)
>     In [86]:mb = array([[[4,5,6]],[[7,8,9]]])         #shape = (2, 1, 3)
> so
>     In [87]:res = dot(ma,mb).shape
>     In [88]:res.shape
>     Out[88]:(2, 3, 2, 3)
> such that
>     res[i,:,j,:] == dot(ma[i,:,:], mb[j,:,:])
> which means that I can take the results I want out of res by slicing
> (somehow) res[0,:,0,:] and res[1,:,1,:] out.
> 
> Is there an easier way, which would make dot only calculate the dot
> products for the cases where i==j (which is what I call threading over
> the first dimension)?

I'm afraid I really don't understand the operation that you want.

> Since the docstring makes no mention of what happens over other
> dimensions, should that be added, or is this the conventional numpy
> behaviour that I need to get used to?

It's fairly conventional for operations that reduce values along an axis to a
single value. The remaining axes are left untouched. E.g.

In [1]: from numpy import *

In [2]: a = random.randint(0, 10, size=(3,4,5))

In [3]: s1 = sum(a, axis=1)

In [4]: a.shape
Out[4]: (3, 4, 5)

In [5]: s1.shape
Out[5]: (3, 5)

In [6]: for i in range(3):
   ...:     for j in range(5):
   ...:         print i, j, (sum(a[i,:,j]) == s1[i,j]).all()
   ...:
   ...:
0 0 True
0 1 True
0 2 True
0 3 True
0 4 True
1 0 True
1 1 True
1 2 True
1 3 True
1 4 True
2 0 True
2 1 True
2 2 True
2 3 True
2 4 True

-- 
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

Re: [Numpy-discussion] Re: Dot product threading?

[Pau G. <pau...@gm...>]

> I'm afraid I really don't understand the operation that you want.

I think that the operation Angus wants is the following (at least I
would like that one ;-)

if you have two 2darrays of shapes:
    a.shape =3D (n,k)
    b.shape =3D (k,m)
you get:
    dot( a,b ).shape =3D=3D (n,m)

Now, if you have higher dimensional arrays  (kind of "arrays of matrices")
    a.shape =3D I+(n,k)
    b.shape =3D J+(k,m)
where I and J are tuples, you get
    dot( a,b ).shape =3D=3D I+J+(n,m)
    dot( a,b )[ i,j ] =3D=3D dot( a[i],b[j] )          #i,j represent tuple=
s
That is the current behaviour, it computes the matrix product between
every possible pair.
For me that is similar to 'outer' but with matrix product.

But sometimes it would be useful (at least for me) to have:
     a.shape =3D I+(n,k)
     b.shape =3D I+(k,m)
and to get only:
     dot2( a,b ).shape =3D=3D I+(n,m)
     dot2( a,b )[i] =3D=3D dot2( a[i], b[i] )

This would be a natural extension of the scalar product (a*b)[i] =3D=3D 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?

Cheers,
pau

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.

Re: [Numpy-discussion] Dot product threading?

[Pau G. <pau...@gm...>]

> 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] =3D=3D b.shape[2]
>         assert a.shape[0] =3D=3D b.shape[0]
>     except AssertionError:
>         print "incorrect input shapes"
>     res =3D zeros( (a.shape[0], a.shape[1], a.shape[1]), dtype=3Dfloat )
>     for i in range(a.shape[0]):
>         res[i,...] =3D dot( a[i,...], b[i,...] )
>     return res
>

yes, that is what I would like. I would like it even with more
dimensions and with all the broadcasting rules ;-)
These can probably be achieved by building actual 'arrays of matrices'
(an array of matrix objects) and then using the ufunc machinery.
But I think that a simple dot2 function (with an other name of course)
will still very useful.

pau

Re: [Numpy-discussion] Dot product threading?

[<jor...@bo...>]

This thread discusses one of the things highest on my wishlist for=20
numpy. I have attached my first attempt at a code that will create a=20
broadcastingfunction that will broadcast a function over arrays where=20
the last N indices are assumed to be for use by the function, N=3D2 would=
=20
be used for matrices. It is implemented for Numeric.


/J=F6rgen


Pau Gargallo skrev:
>> 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] =3D=3D b.shape[2]
>>         assert a.shape[0] =3D=3D b.shape[0]
>>     except AssertionError:
>>         print "incorrect input shapes"
>>     res =3D zeros( (a.shape[0], a.shape[1], a.shape[1]), dtype=3Dfloat=
 )
>>     for i in range(a.shape[0]):
>>         res[i,...] =3D dot( a[i,...], b[i,...] )
>>     return res
>>
>=20
> yes, that is what I would like. I would like it even with more
> dimensions and with all the broadcasting rules ;-)
> These can probably be achieved by building actual 'arrays of matrices'
> (an array of matrix objects) and then using the ufunc machinery.
> But I think that a simple dot2 function (with an other name of course)
> will still very useful.
>=20
> pau
>=20
>=20
> -------------------------------------------------------
> Using Tomcat but need to do more? Need to support web services, securit=
y?
> Get stuff done quickly with pre-integrated technology to make your job=20
> easier
> Download IBM WebSphere Application Server v.1.0.1 based on Apache Geron=
imo
> http://sel.as-us.falkag.net/sel?cmd=3Dk&kid=120709&bid&3057&dat=121642
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>=20

Re: [Numpy-discussion] Re: Dot product threading?

[Bill B. <wb...@gm...>]

In Einstein summation notation, what numpy.dot() does now is:
c_riqk =3D a_rij * b_qjk

And you want:
c_[r]ik =3D a_[r]ij * b_[r]jk

where the brackets indicate a 'no summation' index.
Writing the ESN makes it clearer to me anyway.  :-)

--bb



On 5/18/06, Pau Gargallo <pau...@gm...> wrote:
>
> > I'm afraid I really don't understand the operation that you want.
>
> I think that the operation Angus wants is the following (at least I
> would like that one ;-)
>
> if you have two 2darrays of shapes:
>     a.shape =3D (n,k)
>     b.shape =3D (k,m)
> you get:
>     dot( a,b ).shape =3D=3D (n,m)
>
> Now, if you have higher dimensional arrays  (kind of "arrays of matrices"=
)
>     a.shape =3D I+(n,k)
>     b.shape =3D J+(k,m)
> where I and J are tuples, you get
>     dot( a,b ).shape =3D=3D I+J+(n,m)
>     dot( a,b )[ i,j ] =3D=3D dot( a[i],b[j] )          #i,j represent tup=
les
> That is the current behaviour, it computes the matrix product between
> every possible pair.
> For me that is similar to 'outer' but with matrix product.
>
> But sometimes it would be useful (at least for me) to have:
>      a.shape =3D I+(n,k)
>      b.shape =3D I+(k,m)
> and to get only:
>      dot2( a,b ).shape =3D=3D I+(n,m)
>      dot2( a,b )[i] =3D=3D dot2( a[i], b[i] )
>
> This would be a natural extension of the scalar product (a*b)[i] =3D=3D
> 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?
>
> Cheers,
> pau
>
>

Re: [Numpy-discussion] Re: Dot product threading?

[Fernando P. <fpe...@gm...>]

On 5/18/06, Bill Baxter <wb...@gm...> wrote:
> In Einstein summation notation, what numpy.dot() does now is:
> c_riqk =3D a_rij * b_qjk
>
> And you want:
> c_[r]ik =3D a_[r]ij * b_[r]jk
>
> where the brackets indicate a 'no summation' index.
> Writing the ESN makes it clearer to me anyway.  :-)

I recently needed something similar to this, and being too lazy to
think up the proper numpy ax-swapping kung-fu, I just opened up weave
and was done with it in a hurry.   Here it is, in case anyone finds
the basic idea of any use.

Cheers,

f

###

class mt3_dot_factory(object):
    """Generate the mt3t contract function, holding necessary state.

    This class only needs to be instantiated once, though it doesn't try to
    enforce this via singleton/borg approaches at all."""

    def __init__(self):
        # The actual expression to contract indices, as a list of strings t=
o be
        # interpolated into the C++ source
        mat_ten =3D ['mat(i,m)*ten(m,j,k)', # first index
                   'mat(j,m)*ten(i,m,k)', # second
                   'mat(k,m)*ten(i,j,m)', # third
                   ]
        # Source template
        code_tpl =3D """
        for (int i=3D0;i<order;i++) {
            for (int j=3D0;j<order;j++) {
                for (int k=3D0;k<order;k++) {
                    double sum=3D0;
                    for (int m=3D0;m<order;m++) {
                        sum +=3D %s;
                    }
                    out(i,j,k) =3D sum;
                }
            }
        }
        """
        self.code =3D [code_tpl % mat_ten[idx] for idx in (0,1,2)]

    def __call__(self,mat,ten,idx):
        """mt3s_contract(mat,ten,idx) -> tensor.

        A special type of matrix-tensor contraction over a single index.  T=
he
        returned array has the following structure:

              out(i,j,k) =3D sum_m(mat(i,m)*ten(m,j,k)) if idx=3D=3D0
              out(i,j,k) =3D sum_m(mat(j,m)*ten(i,m,k)) if idx=3D=3D1
              out(i,j,k) =3D sum_m(mat(k,m)*ten(i,j,m)) if idx=3D=3D2

        Inputs:

         - mat: an NxN array.
         - ten: an NxNxN array.
         - idx: the position of the index to contract over, 0 1 or 2."""

        # Minimal input validation - we use asserts so they don't kick in
        # under a -O run of python.
        assert len(mat.shape)=3D=3D2,\
               "mat must be a rank 2 array, shape=3D%s" % mat.shape
        assert mat.shape[0]=3D=3Dmat.shape[1],\
               "Only square matrices are supported: mat shape=3D%s" % mat.s=
hape
        assert len(ten.shape)=3D=3D3,\
               "mat must be a rank 3 array, shape=3D%s" % ten.shape
        assert ten.shape[0]=3D=3Dten.shape[1]=3D=3Dten.shape[2],\
               "Only equal-dim tensors are supported: ten shape=3D%s" % ten=
.shape

        order =3D mat.shape[0]
        out =3D zeros_like(ten)
        inline(self.code[idx],('mat','ten','out','order'),
               type_converters =3D converters.blitz)
        return out

# Make actual instance
mt3_dot =3D mt3_dot_factory()