[Numpy-discussion] ndarray of matrices

[Pau Gargallo]

On 6/8/06, Alexandre Guimond <alexandre.guimond at mirada-solutions.com> wrote:
>
>
>
>
> Hi all.
>
>
>
> i work mainly with "volume" (3d) images, and numpy.ndarray answers most of
> my needs (addition of images, etc.). The problem I'm faced now with is that
> I have images of matrices and vectors and would like that when I do
> image_of_matrices * image_of_vector is does the dot product of each of my
> matrices with all of my vectors, and when I do image_of_matrices.mean() it
> gives me the mean matrix. Basically, I want the same functionalities that
> are currently provided with scalars, but applied to matrices.
>
>
>
> It seems that a nice way of doing this is to have and ndarray of
> numpy.matrix, but this isn't supported it seems. Can anyone recommend a good
> way of implementing this? I'm new with the numpy thing and I'm not sure if
> subclassing ndarray is a good idea since I'll have to overload all the
> operators and i don't believe this will result in a very fast
> implementation, but I might be mistaken. Another possibility may be to
> create a new dtype for numpy.matrix, but I don't know if this is possible.
> Anyone have recommandations?
>
>
>
> Thx for any help.
>

We are several of us wondering which is the best way to do this kind
of things. We were discussing this before
(http://aspn.activestate.com/ASPN/Mail/Message/numpy-discussion/3130104),
and some solutions were proposed, but we still don't have the definitive answer.

Building arrays of matrices objects will be too inefficient.
For me the best thing would be to have n-dimensional universal
functions, but this don't exist yet.

Meanwhile, I am using the following code (which is not *the* solution):

from numpy import *

nz,ny,nx = 1,1,1

im_of_mat = rand( nz, ny, nx, 3,3 )
im_of_vec = rand( nz, ny, nx, 3 )

im_of_products = ( im_of_mat * im_of_vec[...,newaxis,:] ).sum(axis=-1)


# test that everything it's ok
for m,v,p in zip(im_of_mat.reshape(-1,3,3),
		 im_of_vec.reshape(-1,3),
		 im_of_products.reshape(-1,3)):
	assert allclose( dot(m,v), p )



pau