Timothy Hochberg has proposed a generalization of the matrix mechanism
to support manipulating arrays of linear algebra objects. For example,
one might have an array of matrices one wants to apply to an array of
vectors, to yield an array of vectors:
In [88]: A = np.repeat(np.eye(3)[np.newaxis,...],2,axis=0)
In [89]: A
Out[89]:
array([[[ 1., 0., 0.],
[ 0., 1., 0.],
[ 0., 0., 1.]],
[[ 1., 0., 0.],
[ 0., 1., 0.],
[ 0., 0., 1.]]])
In [90]: V = np.array([[1,0,0],[0,1,0]])
Currently, it is very clumsy to handle this kind of situation even
with arrays, keeping track of dimensions by hand. For example if one
wants to multiply A by V "elementwise", one cannot simply use dot:
In [92]: np.dot(A,V.T)
Out[92]:
array([[[ 1., 0.],
[ 0., 1.],
[ 0., 0.]],
[[ 1., 0.],
[ 0., 1.],
[ 0., 0.]]])
tensordot() suffers from the same problem. It arises because when you
combine two multiindexed objects there are a number of ways you can do
it:
A: Treat all indices as different and form all pairwise products
(outer product):
In [93]: np.multiply.outer(A,V).shape
Out[93]: (2, 3, 3, 2, 3)
B: Contract the outer product on a pair of indices:
In [98]: np.tensordot(A,V,axes=(-1,-1)).shape
Out[98]: (2, 3, 2)
C: Take the diagonal along a pair of indices:
In [111]: np.diagonal(np.multiply.outer(A,V),axis1=0,axis2=3).shape
Out[111]: (3, 3, 3, 2)
What we want in this case is a combination of B and C:
In [106]: np.diagonal(np.tensordot(A,V,axes=(-1,-1)),axis1=0,axis2=2).T.shape
Out[106]: (2, 3)
but it cannot be done without constructing a larger array and pulling
out its diagonal.
If this seems like an exotic thing to want to do, suppose instead we
have two arrays of vectors, and we want to evaluate the array of dot
products. None of no.dot, np.tensordot, or np.inner produce the
results you want. You have to multiply elementwise and sum. (This also
precludes automatic use of BLAS, if indeed tensordot uses BLAS.)
Does it make sense to implement a generalized tensordot that can do
this, or is it the sort of thing one should code by hand every time?
Is there any way we can make it easier to do simple common operations
like take the elementwise matrix-vector product of A and V?
The more general issue, of making linear algebra natural by keeping
track of which indices are for elementwise operations, and which
should be used for dots (and how), is a whole other kettle of fish. I
think for that someone should think hard about writing a full-featured
multilinear algebra package (might as well keep track of coordinate
transformations too while one was at it) if this is intended.
Anne
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
