[Numpy-discussion] arrays of matrices

[Geoffrey Irving]

On Thu, Feb 28, 2008 at 06:55:11PM -0600, Robert Kern wrote:
> On Thu, Feb 28, 2008 at 6:43 PM, Geoffrey Irving <irving at pixar.com> wrote:
> >  > The magic is in In[27]. We reshape the array of vectors to be
> >  > compatible with the shape of the array of matrices. When we multiply
> >  > the two together, it is as if we multiplied two (n,3,3) matrices, the
> >  > latter being the vectors repeated 3 times. Then we sum along the rows
> >  > of each of the product matrices to get the desired dot product.
> >
> >  Thanks!  That'll do nicely.
> >
> >  For large matrices, that could be problematic due to the blowup in
> >  intermediate memory, but on the other hand for large matrices a loop
> >  through the toplevel index wouldn't add much cost.
> 
> If you really want to save memory and you can destroy A, then you
> could do the multiplication in-place. If you really want to get fancy
> and can destroy b, you can use it as storage for the summation output,
> too.
> 
> In [11]: A *= b.reshape([n,1,3])
> 
> In [12]: c = A.sum(axis=-1, out=b)
> 
> In [13]: b
> Out[13]:
> array([[   50,   140,   230],
>        [ 1220,  1580,  1940],
>        [ 4010,  4640,  5270],
>        [ 8420,  9320, 10220],
>        [14450, 15620, 16790]])
> 
> In [14]: c is b
> Out[14]: True

By large I meant if the array shapes are (n,i,j), (n,j,k), where all indices
are large.  The permanent memory cost is O(nij+njk), but O(nijk) flops are
required, and the broadcast/sum solution would require intermediate storage
for each one.

It doesn't matter in my case though, so it's just a curiosity.

> >  > PS: Are you perchance the Geoffrey Irving I knew at CalTech, class of '03?
> >
> >  Yep.  That would answer the question I had when I started reading this email.
> >  However, it's spelled Caltech, not CalTech!
> 
> Yeah, yeah, yeah. The Wikis, they have taken over my finger ReFlexes.
> 
> NumPy Rudds += 1. Take that, Tim Hochberg!  :-)

There are a bunch of techers here (from further back than us), but I may be the
only one actively using numpy at the moment.  It would be very painful to have
to lug large mesh data in python around without it.  And I wouldn't have the
pleasure of writing stuff like this:

    # triangulate a polygon mesh, given as (counts,vertices,points)

    # counts is the number of vertices in each polygon, and vertices
    # are the concatenated vertex indices of each polygon

    triangleFace = numpy.arange(len(counts)).repeat(counts-2)
    vertex1 = (numpy.add.accumulate(counts)-counts)[triangleFace]
    vertex2 = numpy.arange(len(triangleFace)) + 2*triangleFace + 1
    triangles = vertices[numpy.vstack([vertex1,vertex2,vertex2+1]).transpose()]

Geoffrey