The docstring shows of numpy.linalg.solve:

a : (…, M, M) array_like
b : {(…, M,), (…, M, K)}, array_like

which leads to ambiguous conclusions. Let's stick to first choice. This would mean for a.shape = (M, M), b.shape = (N, M) I should expect a result of x.shape = (N, M).

Of course, it cannot really work, as there is always an uncertainty, if b should be interpreted as (..., M) or (..., M, K) (after all, what should we do if M==K?). Nevertheless, the documentation should be more explicit, as new users will be easily confused.

I personally would remove the ambiguity completely. One idea would simply be adding an argument e.g. b_is_matrix: bool to specify which signature (..., M) or (..., M, K) should apply.
The other idea would be, to remove (..., M, K) completely, as it can be obtained trivially by passing b.T (or np.swapaxes(b, -1, -2)) and using broadcasting.

Reproducing code example:

import numpy as np

a = np.ones((1, 1))
b = b = np.ones((2, 1))
np.linalg.solve(a, b)

Error message:

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

Numpy/Python version information:

1.18.1 3.6.1 (default, Nov 11 2017, 14:50:47)
[GCC 4.9.2]