Describe the issue:

The Python docs on dunder methods with reflected operands says:
These functions are only called if the left operand does not support the corresponding operation and the operands are of different types.
“Does not support” here means that the class has no such method, or the method returns NotImplemented. Do not set the method to None if you want to force fallback to the right operand’s reflected method—that will instead have the opposite effect of explicitly blocking such fallback.
For operands of the same type, it is assumed that if the non-reflected method – such as add() – fails then the overall operation is not supported, which is why the reflected method is not called.

This behaviour seems to be violated when performing an operation like np.ndarray * np.matrix. This code evaluates like np.matrix.__rmul__(np.ndarray).
Confusingly, np.ndarray.__mul__(np.matrix) exists and can be called explicitly, leading to a different result.
AFAIK dispatching to np.matrix rmul is expected, as both left- and right-handed multiplicaitons on np.matrix are evaluated as matrix multiplication. However, according to the standard, calling through mul should also follow this behaviour, or follow the "Does not support" path described above.

I believe this bug happens because a*b bypasses __getattribute("__mul__")__ due to special method lookup rules - x.__mul__ goes through class __getattribute__, type(x).__mul__ goes through metaclass __getattribute__, while x*foo goes directly to x's __mul__ implementation. Unfortunately I do not have enough experience with the numpy codebase to check if this is really the case here.

Side note:
np.ndarray *= np.matrix leads to behaviour similar to np.ndarray.__mul__(np.matrix), just with the result being an array - this is confusing, but described in the docs and does not violate the standard. Debugging this is how I found the main issue I'm reporting.

Yes, I know np.matrix use is discouraged.
This might be related to #5844

Reproduce the code example:

import numpy as np
nda = np.ones((3,3))
mat = np.matrix((np.ones((3,1))))
infix_op = A*mat
print(infix_op)
# [[3.]
#  [3.]
#  [3.]]
print(type(infix_op))
# <class 'numpy.matrix'>

dunder_instance = A.__mul__(mat)
print(dunder_instance)
# [[1. 1. 1.]
#  [1. 1. 1.]
#  [1. 1. 1.]]
print(type(dunder_instance))
# <class 'numpy.matrix'>

dunder_class = np.ndarray.__mul__(A, mat)
print(dunder_class)
# [[1. 1. 1.]
#  [1. 1. 1.]
#  [1. 1. 1.]]
print(type(dunder_class))
# <class 'numpy.matrix'>

Error message:

No response

Runtime information:

1.23.5
3.10.12 (main, Jun 11 2023, 05:26:28) [GCC 11.4.0]

Context for the issue:

No response