[Numpy-discussion] arrays of matrices

[Robert Kern]

On Thu, Feb 28, 2008 at 4:34 PM, Geoffrey Irving <irving at pixar.com> wrote:
> Hello,
>
>  I have a large number of points (shape (n,3)), and a matching
>  number of 3x3 matrices (shape (n,3,3)), and I want to compute
>  the product of each matrix times the corresponding point.
>
>  I can't see a way to do this operation with dot or tensordot,
>  since these routines either sum across an index or treat it
>  as independent between the two arguments.
>
>  Is this case, I can use the fact that 3 is small to do it, but
>  is there a clean way for handle this kind of "array of matrix"
>  situation in general?

For dot-products, yes. We can use broadcasting followed by axis summation.


In [20]: from numpy import *

In [21]: n = 5

In [22]: A = arange(n*3*3).reshape([n,3,3])

In [23]: b = 10*arange(n*3).reshape([n,3])

In [24]: A
Out[24]:
array([[[ 0,  1,  2],
        [ 3,  4,  5],
        [ 6,  7,  8]],

       [[ 9, 10, 11],
        [12, 13, 14],
        [15, 16, 17]],

       [[18, 19, 20],
        [21, 22, 23],
        [24, 25, 26]],

       [[27, 28, 29],
        [30, 31, 32],
        [33, 34, 35]],

       [[36, 37, 38],
        [39, 40, 41],
        [42, 43, 44]]])

In [25]: b
Out[25]:
array([[  0,  10,  20],
       [ 30,  40,  50],
       [ 60,  70,  80],
       [ 90, 100, 110],
       [120, 130, 140]])

In [26]: for i in range(n):
   ....:     print dot(A[i], b[i])
   ....:
   ....:
[ 50 140 230]
[1220 1580 1940]
[4010 4640 5270]
[ 8420  9320 10220]
[14450 15620 16790]

In [27]: (A * b.reshape([n, 1, 3])).sum(axis=-1)
Out[27]:
array([[   50,   140,   230],
       [ 1220,  1580,  1940],
       [ 4010,  4640,  5270],
       [ 8420,  9320, 10220],
       [14450, 15620, 16790]])


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.


PS: Are you perchance the Geoffrey Irving I knew at CalTech, class of '03?

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
  -- Umberto Eco