[Numpy-discussion] Vectorizing "dot" on the two last axis

[Anne Archibald]

2008/10/10 Gael Varoquaux <gael.varoquaux at normalesup.org>:

> I have been unable to vectorize the following operation::
>
>    window_size = 10
>    nb_windows = 819
>    nb_clusters = 501
>    restricted_series = np.random.random(size=(window_size, nb_clusters,
>                                                    nb_windows))
>    this_cov = np.zeros((nb_windows, nb_clusters, nb_clusters))
>    print >>sys.stderr, "alive"
>    for index, window in enumerate(restricted_series):
>        this_cov[index, ...] = np.dot(window, window.T)
>
>
> The last for loop is the one I am unhappy with.
>
> Note that it is fairly easy to get huge arrays when trying the vectorize
> this through np.dot or np.tensordot.
>
> I am unhappy with myself: I feel that it should be possible to vectorize
> this operation, but I cannot figure it out.

I am pretty sure that, unfortunately, there is no way to vectorize
this without an intermediate array substantially larger than either
the inputs or the outputs. (Add one to the tally of people wishing for
ufunc-like linear algebra.) From a computational point of view, this
isn't particularly a problem: the intermediate array cannot be large
enough to be a problem without the slices along at least one axis
being large enough that for loop overhead is pretty minor. So you can
just loop over this axis. Coding-wise, of course, this is a royal
pain: unless you know ahead of time which axis is the smallest, you
need to code several versions of your routine, and each version must
include a for loop (just what numpy is supposed to help avoid).

So: ufunc-like linear algebra would be great, but in the meantime, I
guess we all learn to live with for loops.

Anne