[Numpy-discussion] Repeated dot products

Anne Archibald peridot.faceted@gmail....
Sat Dec 12 16:30:46 CST 2009

2009/12/12 T J <tjhnson@gmail.com>:
> Hi,
> Suppose I have an array of shape:  (n, k, k).  In this case, I have n
> k-by-k matrices.  My goal is to compute the product of a (potentially
> large) user-specified selection (with replacement) of these matrices.
> For example,
>   x = [0,1,2,1,3,3,2,1,3,2,1,5,3,2,3,5,2,5,3,2,1,3,5,6]
> says that I want to take the 0th matrix and dot it with the 1st matrix
> and dot that product with the 2nd matrix and dot that product with the
> 1st matrix again and so on...
> Essentially, I am looking for efficient ways of doing this.  It seems
> like what I *need* is for dot to be a ufunc with a reduce() operator.
> Then I would construct an array of the matrices, as specified by the
> input x.  For now, I am using a python loop and this is unbearable.
>>>> prod = np.eye(k)
>>>> for i in x:
> ...  prod = dot(prod, matrices[i])
> ...
> Is there a better way?

You are right that numpy/scipy should have a matrix product ufunc. In
fact it should have ufunc versions of all the linear algebra tools. It
even has the infrastructure (generalized ufuncs) to support such a
thing in a very general way. Sadly this infrastructure is basically

Your best workaround depends on the sizes of the various objects
involved. If k is large enough, then a k by k by k matrix multiply
will be slow enough that it doesn't really matter what else you do. If
x is much longer than n, you will want to avoid constructing a len(x)
by k by k array. If n is really small and x really large, you might
even win by some massively clever Knuthian operation that found
repeated substrings in x and evaluated each corresponding matrix
product only once. If on the other hand len(x) is much shorter than n,
you could perhaps benefit from forming an intermediate len(x) by k by
k array. This array could then be used in a vectorized reduce
operation, if we had one. Normally, one can work around the lack of a
vectorized matrix multiplication by forming an (n*k) by k matrix and
multiplying it, but I don't think this will help any with reduce. If k
is small, you can multiply two k by k matrices by producing the k by k
by k elementwise product and then adding along the middle axis, but I
don't think this will work well with reduction either.

So I have just two practical solutions, really:
(a) use python reduce() and the matrix product, possibly taking
advantage of the lack of need to produce an intermediate len(x) by k
by k matrix, and live with python in that part of the loop, or
(b) use the fact that in recent versions of numpy there's a
quasi-undocumented ufuncized matrix multiply built as a test of the
generalized ufunc mechanism. I have no idea whether it supports
reduce() (and if it doesn't adding it may be an unpleasant

Good luck,

P.S. if you come up with a good implementation of the repeated
substring approach I'd love to hear it! -A

More information about the NumPy-Discussion mailing list