[Numpy-discussion] array of matrices

[Robert Kern]

On Sat, Mar 28, 2009 at 23:15, Anne Archibald <peridot.faceted at gmail.com> wrote:
> 2009/3/28 Geoffrey Irving <irving at naml.us>:
>> On Sat, Mar 28, 2009 at 12:47 AM, Robert Kern <robert.kern at gmail.com> wrote:
>>> 2009/3/27 Charles R Harris <charlesr.harris at gmail.com>:
>>>>
>>>> On Fri, Mar 27, 2009 at 4:43 PM, Robert Kern <robert.kern at gmail.com> wrote:
>>>>>
>>>>> On Fri, Mar 27, 2009 at 17:38, Bryan Cole <bryan at cole.uklinux.net> wrote:
>>>>> > I have a number of arrays of shape (N,4,4). I need to perform a
>>>>> > vectorised matrix-multiplication between pairs of them I.e.
>>>>> > matrix-multiplication rules for the last two dimensions, usual
>>>>> > element-wise rule for the 1st dimension (of length N).
>>>>> >
>>>>> > (How) is this possible with numpy?
>>>>>
>>>>> dot(a,b) was specifically designed for this use case.
>>>>
>>>> I think maybe he wants to treat them as stacked matrices.
>>>
>>> Oh, right. Sorry. dot(a, b) works when a is (N, 4, 4) and b is just
>>> (4, 4). Never mind.
>>
>> It'd be great if this operation existed as a primitive.  What do you
>> think would be the best way in which to add it?  One option would be
>> to add a keyword argument to "dot" giving a set of axes to map over.
>> E.g.,
>>
>>    dot(a, b, map=0) = array([dot(u,v) for u,v in zip(a,b)]) # but in C
>>
>> "map" isn't a very good name for the argument, though.
>
> I think the right long-term solution is to make dot (and some other
> linear algebra functions) into "generalized ufuncs", so that when you
> dot two multidimensional objects together, they are treated as arrays
> of two-dimensional arrays, broadcasting is done on all but the last
> two dimensions, and then the linear algebra is applied "elementwise".
> This covers basically all "stacked matrices" uses in a very general
> way, but would require some redesigning of the linear algebra system -
> for example, dot() currently works on both two- and one-dimensional
> arrays, which can't work in such a setting.
>
> The infrastructure to support such generalized ufuncs has been added
> to numpy, but as far as I know no functions yet make use of it.

I don't think there is a way to do it in general with dot(). Some
cases are ambiguous. I think you will need separate matrix-matrix,
matrix-vector, and vector-vector gufuncs, to coin a term.

-- 
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