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

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.mcmorland@auckland.ac.nz 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.mcmorland@auckland.ac.nz 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.

This thread discusses one of the things highest on my wishlist for numpy. I have attached my first attempt at a code that will create a broadcastingfunction that will broadcast a function over arrays where the last N indices are assumed to be for use by the function, N=2 would be used for matrices. It is implemented for Numeric. /Jörgen 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] == 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

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

------------------------------------------------------- Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with pre-integrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.as-us.falkag.net/sel?cmd=k&kid0709&bid&3057&dat1642 _______________________________________________ Numpy-discussion mailing list Numpy-discussion@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/numpy-discussion

# -*- coding: ISO-8859-1 -*- from exceptions import ArithmeticError import Numeric,LinearAlgebra,pdb from Numeric import array,NewAxis,repeat,zeros from LinearAlgebra import inverse import operator def transpose(array): """transpose list of list """ return Numeric.transpose(array).tolist() def outer(data): result=[[]] for l in data: newresult=[] for x in result: for e in l: newresult.append(x+[e]) result=newresult return result def is_squaremultiarray(arraymatrix): try: arraymatrix.shape except AttributeError: return False if len(arraymatrix.shape)>=2: return arraymatrix.shape[-1]==arraymatrix.shape[-2] else: return False def is_multiarray(arraymatrix): try: arraymatrix.shape except AttributeError: return False if len(arraymatrix.shape)>=2: return True else: return False def broad(x): """broad(x) will pick out the highest value of the list x. There can be one or two distinct values in x one of which must be one >>> broad([1,3,3]) 3 >>> broad([1]) 1 >>> broad([1]) 1 >>> broad([]) Traceback (most recent call last): ... IndexError: Incompatible indices for broadcasting >>> broad([1,2,3]) Traceback (most recent call last): ... IndexError: Incompatible indices for broadcasting """ if len(x)==1: return x[0] elif x==[1]: return 1 q=set(x) if q==set([1]): return 1 q.discard(1) if len(q)==1: return q.pop() else: raise IndexError("Incompatible indices for broadcasting") def make_index(x,broad): if x==1: return [0 for x in range(broad)] else: return range(x) def gen_iter_idx(shapes): result=[] maxidx=map(max,zip(*shapes)) for x in zip(*shapes): result.append([make_index(q,max(x)) for q in x]) return zip(*map(outer,zip(*result))) def gen_result_shape(shapes): maxl=max(map(len,shapes)) shapes=[x for x in shapes if len(x)==maxl] return tuple(map(max,zip(*shapes))) def broadcast_shapes(matrices,dims=2): shapes=[m.shape[:-dims] for m in matrices] maxlen=max([len(x) for x in shapes]) shapes=[tuple(1 for i in range(maxlen-len(x)))+tuple(x) for x in shapes] return shapes def make_flat(mat,dims=2): """make_flat(mat) takes list of arrays and changes shape of arrays to only have three indices the last two are assumed to represent a matrix convert x.shape=(n1,n2,...,nn,i,j) to x.shape=(n1*n2...*nn,i,j) return old list of old shapes shapes """ old=[] olddict={} newdict={} for m in mat: if id(m) not in olddict: #make sure shape is changed only once oldshape=m.shape olddict[id(m)]=oldshape totlen=reduce(operator.mul,m.shape[:-dims],1) shape=(totlen,)+m.shape[-dims:] newdict[id(m)]=shape else: shape=newdict[id(m)] oldshape=olddict[id(m)] old.append(oldshape) m.shape=shape return old def make_broad(mat,dims=2): """make_broad(mat,dims) change shape of arrays in the list mat so they all have the same number of dimensions. To simplify further operations that want to do broadcasting operations assuming the last <dims> are part of the object being used for operations. dims=2 would be used for operations on matrices. """ old=[] nshape=broadcast_shapes(mat,dims=dims) olddict={} newdict={} for m,s in zip(mat,nshape): if id(m) not in olddict: #make sure shape is changed only once oldshape=m.shape if s==tuple(): newshape=(1,)+m.shape[-dims:] else: newshape=s+m.shape[-dims:] newdict[id(m)]=newshape olddict[id(m)]=oldshape else: newshape=newdict[id(m)] oldshape=olddict[id(m)] m.shape=newshape old.append(oldshape) return old def unmake_flat(mat,old): for m,s in zip(mat,old): m.shape=s def broad_fun(infunc,dims=2,retvalue=None): """broad_fun(infunc,dims,retvalue) returns a new function based on infunc that does broadcasting over the n-dims dimensions of its input variables. retvalue is used to specify which returnvalue is used for infuncs that return more than one value. example: >>> from Numeric import matrixmultiply,array,NewAxis >>> mul=broad_fun(matrixmultiply,2) >>> x=array([[1.,2],[3,4]]) >>> y=repeat(repeat(x[NewAxis,NewAxis],3,0),2,1) >>> y=y*array([1,2,3])[:,NewAxis,NewAxis,NewAxis] >>> mul(x,y) array([[[[ 7., 10.], [ 15., 22.]], [[ 7., 10.], [ 15., 22.]]], [[[ 14., 20.], [ 30., 44.]], [[ 14., 20.], [ 30., 44.]]], [[[ 21., 30.], [ 45., 66.]], [[ 21., 30.], [ 45., 66.]]]]) """ if retvalue is not None: def func(*mat): return infunc(*mat)[retvalue] else: func=infunc def newfunc(*mat): oldshape=make_flat(mat) testresult=func(*[m[0] for m in mat]) unmake_flat(mat,oldshape) oldshape=make_broad(mat) shapes=[m.shape[:-dims] for m in mat] try: resultshape=gen_result_shape(shapes)+testresult.shape utshape=gen_result_shape([m[:-dims] for m in oldshape])+testresult.shape result=zeros(resultshape,testresult.typecode()) except AttributeError: resultshape=gen_result_shape(shapes) utshape=gen_result_shape([m[:-dims] for m in oldshape]) if utshape==tuple(): return testresult result=zeros(resultshape).astype(type(testresult)) indices=gen_iter_idx(shapes+[resultshape]) for i in indices: residx=i[-1] result[residx]=func(*[m[idx] for m,idx in zip(mat,i[:-1])]) unmake_flat(mat,oldshape) result.shape=utshape return result newfunc.__name__="broadcasting_%s"%infunc.__name__ return newfunc matmul=broad_fun(Numeric.matrixmultiply) det=broad_fun(LinearAlgebra.determinant) inv=broad_fun(LinearAlgebra.inverse) geninv=broad_fun(LinearAlgebra.generalized_inverse) eigenvalues=broad_fun(LinearAlgebra.eigenvalues) linear_least_squares=broad_fun(LinearAlgebra.linear_least_squares,retvalue=0) solve_linear_equations=broad_fun(LinearAlgebra.solve_linear_equations) def _test(): import doctest doctest.testmod() if __name__=="__main__": x=array([[1.,2],[3,4]]) z=array([ [[1.,2],[3,4]], [[4.,3],[2,1]]]) y=repeat(repeat(x[NewAxis,NewAxis],3,0),2,1) a=array([[1,2],[1.1,2.05],[0.98,1.97]]) b=array([[.7],[.75],[.69]]) A=repeat(repeat(a[NewAxis,NewAxis],3,0),2,1)

On 5/18/06, Bill Baxter <wbaxter@gmail.com> wrote:

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

And you want: c_[r]ik = 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 to be # interpolated into the C++ source mat_ten = ['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 = """ for (int i=0;i<order;i++) { for (int j=0;j<order;j++) { for (int k=0;k<order;k++) { double sum=0; for (int m=0;m<order;m++) { sum += %s; } out(i,j,k) = sum; } } } """ self.code = [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. The returned array has the following structure: out(i,j,k) = sum_m(mat(i,m)*ten(m,j,k)) if idx==0 out(i,j,k) = sum_m(mat(j,m)*ten(i,m,k)) if idx==1 out(i,j,k) = sum_m(mat(k,m)*ten(i,j,m)) if idx==2 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)==2,\ "mat must be a rank 2 array, shape=%s" % mat.shape assert mat.shape[0]==mat.shape[1],\ "Only square matrices are supported: mat shape=%s" % mat.shape assert len(ten.shape)==3,\ "mat must be a rank 3 array, shape=%s" % ten.shape assert ten.shape[0]==ten.shape[1]==ten.shape[2],\ "Only equal-dim tensors are supported: ten shape=%s" % ten.shape order = mat.shape[0] out = zeros_like(ten) inline(self.code[idx],('mat','ten','out','order'), type_converters = converters.blitz) return out # Make actual instance mt3_dot = mt3_dot_factory()