Use _matmul_ operator (@) for matrix multiplication, map composition, tensor contraction

[mkoeppe]

#22760 added support for __matmul__ in the coercion model.

We should start using it.

First step: Review the semantics of this operator in major Python software for matrix and tensor computation (NumPy, Numba, TensorFlow, PyTorch, ...) so that we do not paint ourselves into a corner.

Follow-up tickets:

• sage.geometry.hyperbolic_space: __matmul__ for composition of isometries #32212 sage.geometry.hyperbolic_space

CC: @tscrim @egourgoulhon

Component: linear algebra

Author: Matthias Koeppe

Branch/Commit: u/mkoeppe/use__matmul__operator____ @ da5104c

Issue created by migration from https://trac.sagemath.org/ticket/30244

[nbruin]

comment:1

I think in most cases we don't need it. In a matrix algebra, it's pretty clear that multiplication means matrix multiplication; not Hadamard product.

The first place I'd see a possible application is for self-maps of a ring, where one may want to consider both composition and point-wise products.

[mkoeppe]

comment:2

Replying to @nbruin:

The first place I'd see a possible application is for self-maps of a ring, where one may want to consider both composition and point-wise products.

+1. Wondering now if we need a category of @-itive semigroups.

[jhpalmieri]

comment:3

I think using @ for tensor products would be nicer than tensor([M1,M2,M3]). Aside from the fact that it's shorter, Sage's current implementation is a bit of a mess. For example, I can never remember that it should be tensor([M1, M2, M3]), not tensor(M1, M2, M3), or you could use M1.tensor(M2, M3) but not M1.tensor([M2, M3]).

[tscrim]

comment:4

Replying to @jhpalmieri:

I think using @ for tensor products would be nicer than tensor([M1,M2,M3]). Aside from the fact that it's shorter, Sage's current implementation is a bit of a mess. For example, I can never remember that it should be tensor([M1, M2, M3]), not tensor(M1, M2, M3), or you could use M1.tensor(M2, M3) but not M1.tensor([M2, M3]).

+1 for this. This would be a good place to start providing an implementation.

[nbruin]

comment:6

I think '@' for tensor products would be very confusing for other uses in python. While it was introduced non-prescriptively as another binary operator, the implementation name __matmul__ suggests otherwise, and the numpy reason for requesting it was very solidly for matrix multiplication, i.e., map composition.

It would seem to me that the place where having operator notation for tensor products might be interesting is for Kronecker products of matrices. That's exactly the place where it directly clashes with numpy notation!! I also think it's very inefficient for multivariate Kronecker products, because the infix notation necessitates the construction of intermediate results.

For modules, I would expect tensor products and homs to be equally prevalent, so keeping a symmetry in their notation seems like a desirable thing to have.

So, I think "@" is a poor fit for tensor product.

[mkoeppe]

comment:7

+1 on giving the "map composition" priority for this operator.

In the realm of tensors, matrix multiplication would generalize to tensor contraction, not tensor product.

[mkoeppe]

comment:8

Replying to @mkoeppe:

Wondering now if we need a category of @-itive semigroups.

Any thoughts on this? Do we need a category SemigroupsWithRespectToMatmul?

[nbruin]

comment:9

I'd say no, for the same reason why we don't have this for "|", "&", "<<", ">>", and other binary operators that python provides. I think we first need a convincing use-case before we start building infrastructure. Code is a burden, not an asset :-).

[mkoeppe]

comment:10

Note that these other binary operators do not actually participate in the coercion framework...

But great point, of course, that we should start with something concrete first.

[mkoeppe]

Branch: u/mkoeppe/use__matmul__operator____

[nbruin]

comment:12

I think you may want to reconsider #22760 in the light of possible usage scenarios. The coercion framework is particularly designed to figure out common parents into which the operands can be coerced so that the operation can be applied. This does not apply to all cases; for instance, for actions there is no appropriate common parent. For actions, other procedures are followed (and generally, less powerful ones; which is appropriate).

If @ is going to be composition, it's going to be mostly a partial operation if regarded as, say, an operation on homomorphisms between modules over a field. Alternatively, it's an operation that combines objects in DIFFERENT parents (e.g., a pairing $\mathrm{Hom}(A,B) \times \mathrm{Hom}(B,C) \to \mathrm{Hom(}A,C)$) in which case the coercion framework probably doesn't have an appropriate setting yet. I think you really want to know what the coercion framework is supposed to accomplish for you before you try and hook @ into it. If the only scenarios where things work are the cases where _matmul_ already knows what to do, then there's no benefit from an extra indirection layer: you could just put the login into __matmul__ directly.

---

New commits:

ea74b6f	Add support for `__matmul__` in the coercion model
d939b2c	Update doctests for py3
e8d7924	Merge branch 't/22760/add_support_for___matmul___in_the_coercion_model' into t/30244/use__matmul__operator____
a982dd8	sage.categories.map.Map: Add __matmul__

[nbruin]

Commit: a982dd8

[mkoeppe]

comment:13

This is a great point.

On this branch, as you perhaps saw, I am already overriding the (double-underscore) __matmul__ operator, so the coercion framework is not actually involved at all when it comes to the @ operation between two Maps. So in this case, I don't think there is actually any overhead/indirection.

[mkoeppe]

comment:14

But a concern is that by overriding it, it is disabling coerce actions...

[tscrim]

comment:15

I don't think there is any harm in #22760 by adding the hook for it. However, I thought the coercion framework was designed to also handle actions. For example, you need coercion for z * A when z is an integer but A only has a QQ-action. I also thought the coercion framework was what called the different actions to see what was appropriate too.

That being said, I see your point about _matmul_ potentially not being so useful by itself. However, we might as well include it in case someone does have a use case for it.

[egourgoulhon]

comment:16

Replying to @mkoeppe:

In the realm of tensors, matrix multiplication would generalize to tensor contraction, not tensor product.

+1

[nbruin]

comment:17

Replying to @tscrim:

I don't think there is any harm in #22760 by adding the hook for it. However, I thought the coercion framework was designed to also handle actions. For example, you need coercion for z * A when z is an integer but A only has a QQ-action. I also thought the coercion framework was what called the different actions to see what was appropriate too.

Indeed, some coercion steps are possible for actions, but not nearly as much as for operations internal to structures. For instance, for addition between ZZ[x,y] and QQ[z], the coercion system will construct a common covering structure QQ[x,y,z] by combining a sequence of "construction functors" according to certain (heuristic!) rules. I don't think the rules for actions are nearly as advanced -- probably no coercions on the acted-upon set are tried at all; and probably shouldn't.

I'd hope there are better tools available for exploring what coercion to take than to query and see "what works". If that's actually what happens, then applying it to partial operators such as composition is definitely inappropriate.

In general, I'm not so sure coercion will help for @ if it goes the composition route, and then having the hook in the system is going to be counterproductive, because people will stumble on it and do unhelpful things with it. So I disagree with the idea that putting hooks just in case someone finds a use for it is harmless. Not making a design decision on it now can also mean not making a design mistake now. Avoiding mistakes has benefits.

[mkoeppe]

comment:18

I think I will try out an alternative implementation of composition as a coerce-action (keeping the inherited __mul__ operator). Then we can experiment with this a little to gather more insights.

[mkoeppe]

comment:20

Replying to @jhpalmieri:

I think using @ for tensor products would be nicer than tensor([M1,M2,M3]). Aside from the fact that it's shorter, Sage's current implementation is a bit of a mess. For example, I can never remember that it should be tensor([M1, M2, M3]), not tensor(M1, M2, M3), or you could use M1.tensor(M2, M3) but not M1.tensor([M2, M3]).

As notation for tensor products, perhaps we can use a different operator ... how about the bitwise-and operator &?

[mkoeppe]

comment:22

Setting new milestone based on a cursory review of ticket status, priority, and last modification date.

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

c08e5b5

sage.categories.map.Map: Add __matmul__

[sagetrac-git]

Changed commit from c08e5b5 to 63b8892

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

9588f52

FiniteSetEndoMap*: Define `__matmul__`, delegate to it from __mul__

[sagetrac-git]

Changed commit from 63b8892 to 9588f52

[mkoeppe]

Author: Matthias Koeppe

[sagetrac-git]

Changed commit from 9588f52 to f626394

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

f626394

TensorWithIndices: Make `__matmul__` an alias of __mul__

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

3006833	FiniteSetEndoMap*: Define `__matmul__`, delegate to it from __mul__
de14c34	TensorWithIndices: Make `__matmul__` an alias of __mul__
451eb95	Map.__mul__: Add doctest output

[sagetrac-git]

Changed commit from f626394 to 451eb95

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

348f680

TensorWithIndices: Update doctests

[sagetrac-git]

Changed commit from 451eb95 to 348f680

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

878e923	CoercionModel.verify_coercion_maps: Use @ instead of * for composition
da5104c	PrecomposedAction.__init__: Use @ instead of * for composition

[sagetrac-git]

Changed commit from 348f680 to da5104c

[mkoeppe]

comment:32

Not sure if @ should also be used for matrix-vector multiplication. The current code on the branch does not do this; but scipy seems to think so (see for example https://docs.scipy.org/doc/scipy/reference/optimize.linprog-highs.html)