[Numpy-discussion] Preliminary thoughts on implementing __matmul__

Nathaniel Smith njs at pobox.com
Wed Aug 6 08:31:36 EDT 2014

On Wed, Aug 6, 2014 at 2:19 AM, Charles R Harris
<charlesr.harris at gmail.com> wrote:
> Hi All,
> I've been looking to implement the "@" operator from Python 3.5. Looking at
> the current implementation of the dot function, it only uses a vector inner
> product, which is either that defined in arraytypes.c.src or a version using
> cblas defined in _dotblas for the float, cfloat, double, cdouble types. I
> note that the versions defined in arraytypes.c.src include all the numeric
> types plus boolean, datetime, timedelta, and object. I'm not clear why
> datetime and timedelta should have dot products, except perhaps for scalar
> multiplication.

I guess numeric @ timedelta is at least well-defined, but dot products
on datetime make no sense -- datetimes do not support +!

One thing we should keep in mind as well is how to allow user-defined
dtypes to provide efficient matmul implementations.

> The boolean version has the advantage that it can short
> circuit. I also note that all the operations proposed for "@" can easily be
> done with einsum except for objects. So I'm wondering if one easy way to
> implement the functions is to extend einsum to work with objects and make it
> use blas when available.

Those do seem like nice features regardless of what we do for @ :-).

I think the other obvious strategy to consider, is defining a 'dot'
gufunc, with semantics identical to @. (This would be useful for
backcompat as well: adding/dropping compatibility with older python
versions would be as simple as mechanically replacing a @ b with
newdot(a, b) or vice-versa.) This would require one new feature in the
gufunc machinery: support for "optional core axes", to get the right
semantics for 1d arrays. OTOH this would also be useful in general
because there are other gufuncs that want to handle 1d arrays the same
way @ does -- e.g., 'solve' variants.

This would automatically solve both the user-defined dtype problem
(ufuncs already allow for new loops to be registered) and the
third-party array type problem (via __numpy_ufunc__).

> Another thing that may be worth looking into would be some way to multiply
> by the complex conjugate, as that is easy to implement at the low level. I'd
> welcome any thoughts as to how that might be done.

One idea that's come up before was to define a complex-conjugate
dtype, which would allow .H to be a view on the original array.

A simpler solution would be to define a specialized conjdot gufunc.


Nathaniel J. Smith
Postdoctoral researcher - Informatics - University of Edinburgh

More information about the NumPy-Discussion mailing list