[Numpy-discussion] working on multiple matrices of the same shape

[Charles R Harris]

On Fri, Nov 21, 2008 at 2:57 AM, Sébastien Barthélemy
<barthelemy at crans.org>wrote:

> Hello,
>
> I would like to port a matlab library which provides functions for
> rigid body mechanics such as operations on homogeneous matrices (in
> SE(3)), twists (in se(3)) and so.
>
> In matlab, the library worked on 3d matrices: n homoneous matrices
> were stacked along the 3d dimension. This speeded up computations on
> multiple matrices,  at the price that native matrix opertators such H1
> * H2 did not work (it is not defined for 3d matrices).
>
> In this spirit, in numpy a set of rotation matrices could be built in
> the following way:
>
> def rotx(theta):
>    """
>    SE(3) matrices corresponding to a rotation around x-axis. Theta is
> a 1-d array
>    """
>    costheta = np.cos(theta)
>    sintheta = np.sin(theta)
>    H = np.zeros((theta.size,4,4))
>    H[:,0,0] = 1
>    H[:,3,3] = 1
>    H[:,1,1] = costheta
>    H[:,2,2] = costheta
>    H[:,2,1] = sintheta
>    H[:,1,2] = sintheta
>    return H
>
> I'm now seeking advices regarding an implementation with numpy (it's
> my first time with numpy).
>
> - Is there any difference between a 3d array and a 1-d array of 2-d
> arrays  (seems not)
>
> - would you use a 3d-array  or a list of 2d arrays or anything else to
> store these matrices ?
>
> - Is there a way to work easily on multiple matrices of the same shape
> ? (in single-instruction/multiple-data spirit)
>  for instance if A.shape==(3,2,4) and B.shape=(3,4,1), what is the
> best way to compute C
>  for i = (0,1,2)
>      C[i,:,:] = np.dot(A[i,:,:],B[i,:,:])
>

There is work going on along these lines called generalized ufuncs, but it
isn't there yet. The easiest current workaround is to use lists of matrices
or arrays. Or you can get clever and do the matrix multiplications yourself
using the current numpy operators, i.e.

C = (A[:,:,:,newaxis]*B[:,newaxis,:,:]).sum(axis=-2)

But this uses more memory than the list of matrices approach.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20081121/3dd98df0/attachment.html>