Broadcasting rules for linalg.solve are different to those for matmul

[subhylahiri]

The issue is with the line

numpy/numpy/linalg/linalg.py

Line 402 in d35cd07

if b.ndim == a.ndim - 1:

matmul uses the matrix-matrix interpretatin most of the time. It expects the last axis of a to match the second-last axis of b:

>>> import numpy as np
>>> rng = np.random.default_rng()
>>> a, b, c = rng.random((3, 2, 2)), rng.random((2, 3)), rng.random((3, 2))
>>> a @ b
array([[[...]]])
>>> a @ c
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0...

matmul only uses the vector version if an array has exactly one dimension.

>>> d = rng.random(2)
>>> a @ d
array([[...]])

With solve, on the other hand, it will use the vector version when b has one fewer dimension than a:

>>> from numpy.linalg import solve
>>> solve(a, b)
ValueError: solve1: Input operand 1 has a mismatch in its core dimension 0...
>>> solve(a, c)
array([[...]])
>>> solve(a, d)
ValueError: solve: Input operand 1 does not have enough dimensions...

I get the argument for the current behavior (as far as I can tell it means that solve(a, a @ b) == b always works), but it catches me by surprise whenever the errors above happen. I can also concoct examples where matmul and solve fail to be inverses:

>>> e, f = rng.random((2, 2)), rng.random((2, 2, 2))
>>> np.allclose(f @ solve(f, e), e)
False
>>> np.allclose(solve(f, f @ e), e)
True

though I'm not sure that there's any way to ensure they are always inverses.
This means we cannot always use solve instead of inv as is often recommended. In any case, there are ways around this using newaxis and squeeze`.
At the very least, it would be nice if the documentation were more specific about when the vector and matrix versions are used. Currently, it says " the last 1 or 2 dimensions of a multidimensional array are interpreted as vectors or matrices, as appropriate for each operation". It doesn't say which version has priority when both interpretations are possible.

numpy version: 1.26.4
python version 3.12.3

[jakevdp]

Thanks! This issue is discussed in #15349, and will be fixed in NumPy 2.0, where the batch semantics of solve will match those of matmul.

[seberg]

If it's changed in NumPy 2, no reason to keep this open. Thanks both.