How to write a class that makes numpy's dot() use my __rmatmul__()?

Question

What is the simplest way to make numpy's dot() method use my object's __rmatmul__ method?

The issue is easily demonstrated with scipy.sparse:

import numpy as np
from scipy import sparse

A = np.array([[1, 2], [3, 4]]) 
B = sparse.diags([1, 2])

BA = B.dot(A)          # Works fine
AB = B.__rmatmul__(A)  # Works: This is what I want to happen
AB = B.T.dot(A.T).T    # Works but not simple
AB = A @ B             # Preferred way for Python 3.5+
AB = A.dot(B)          # This is what people do, but it gives very non-intuitive results

The desired result here is B.__rmatmul__(A) but I do not want users used to using numpy.dot to get tripped up when the reasonable expression A.dot(B) returns a 2x2 array where each element is a sparse matrix (a copy of B):

>>> AB
array([[<2x2 sparse matrix of type '<class 'numpy.float64'>'
    with 2 stored elements (1 diagonals) in DIAgonal format>,
        <2x2 sparse matrix of type '<class 'numpy.float64'>'
    with 2 stored elements (1 diagonals) in DIAgonal format>],
       [<2x2 sparse matrix of type '<class 'numpy.float64'>'
    with 2 stored elements (1 diagonals) in DIAgonal format>,
        <2x2 sparse matrix of type '<class 'numpy.float64'>'
    with 2 stored elements (1 diagonals) in DIAgonal format>]], dtype=object)

If I am writing my own "sparse" class, how can I get numpy to use my objects __rmatmul__ function? From the answer here, it seems like the only solution might be to make my object a subclass of np.ndarray, but this seems way too complicated (and may still fail since I doubt I can provide sensible view-like access to a sparse array).

(The particular application I would like to use is a custom sparse representation of an array as a sum over dyads, but the details are not important.)

References

This issue has been pointed out several times, and the result explained in terms of the behaviour of numpy.dot(), but the solutions are all workarounds. None answer this question.

People running into this issue:

• scipy.sparse dot extremely slow in Python
• numpy.dot() gives TypeError: can't multiply sequence by non-int of type 'float'
• Dot product between scipy sparse matrix and numpy array give ValueError

Ugly workarounds:

• numpy: formal definition of "array_like" objects?: Formal definition of "array-like" objects as far as numpy is concerned, but this seems to require that my object behave as a dense array of numbers (defeating the purpose of a sparse implementation).
• Array and __rmul__ operator in Python Numpy: Subclass ndarray (with warnings).
• Commute numpy sparse matrix dot product: Suggests essentially B.T.dot(A.T).T.

(Update)

I think in Python 3.5 and above it is clear that one should used the notation A @ B instead of dot(). The question still remains for earlier versions of python, but is less relevant.

• Why is a.dot(b) faster than a@b although Numpy recommends a@b: It is not significantly slower for large enough matrices.
• Difference between numpy dot() and Python 3.5+ matrix multiplication @:
• PEP 465 -- A dedicated infix operator for matrix multiplication
• NumPy dot documentation: Recommends A @ B or matmul() for 2-D arrays.