Issue with current documentation:

The numpy.linalg docs claim, "Returned shape is identical to b [the second argument]." It isn't, though. Returned ndim is identical to b.

e.g.

import numpy as np

A = np.random.normal(size=(2,3,3))
b = np.random.normal(size=(1,3))

x = np.linalg.solve(A, b)
x.shape

gives (2, 3). More generally, the shape of b can be "{(…, M,), (…, M, K)}", which I think is what requires linalg to break the normal rules for broadcasting, which would begin from the right and pad missing axes. E.g.

A = np.random.normal(size=(2,2,3,3))
b = np.random.normal(size=(1,3))

x = np.linalg.solve(A, b)

ValueError: solve: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (m,m),(m,n)->(m,n) (size 1 is different from 3)

I would expect the (1,) to be broadcast to (2, 2) as it would in normal array broadcasting, but linalg doesn't expand missing axes. In another case, if I want b to be of shape (…, M, K), I would expect the following to have shape (3, 3, 3)

A = np.random.normal(size=(3,3,3))
b = np.random.normal(size=(3,3))

x = np.linalg.solve(A, b)

but it has shape (3, 3).

So it feels clear what's happening here.

Idea or request for content:

It would be great if the wording

the last 1 or 2 dimensions of a multidimensional array are interpreted as vectors or matrices, as appropriate for each operation.

in the linalg section on broadcasting could be made more explicit that the ... are broadcast from the left without padding missing axes.