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

[Sébastien Barthélemy]

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,:,:])

- is it a good idea/not to hard to subclass ndarray with a
HomogenousMatrix class ? (one could then redefine inv() for instance)

- Is my implementation of rotx efficient ?

Thank you for any help, and please feel free to send me back reading
anything I missed.

regards
-- 
Sébastien Barthélemy