BUG: Missing ``matmat`` coverage in ``scipy.sparse.linalg.aslinearoperator``

[mfeuerle]

Describe your issue.

When calling scipy.sparse.linalg.aslinearoperator with an matrix A that implements all of the methods matvec, matmat, rmatvec and rmatmat the resulting scipy.sparse.linalg.LinearOperator uses the provided methods for matvec, rmatvec and rmatmat but not for matmat where the default implementation is used instead of the user-provided implementation. This is annoying since the default implementation is inefficient and falls back to several matvec calls.

This should be changed such that the user provided implementation is used.

Reproducing Code Example

from scipy.sparse.linalg import aslinearoperator
class Matrix:
    shape = (1,1)
    dtype = float
    
    def matvec(self, x):
        print("matvec")

    def matmat(self, x):
        print("matmat")
    
    def rmatvec(self, x):
        print("rmatvec")
        
    def rmatmat(self, x):
        print("rmatmat")
        
A = aslinearoperator(Matrix())
A._matvec(1)
A._rmatvec(1)
A._rmatmat(1)
A._matmat(1)

Error message

matvec
rmatvec
rmatmat
Traceback (most recent call last):
  File *** in <module>
    A._matmat(1)
  File "/venv/lib/python3.9/site-packages/scipy/sparse/linalg/_interface.py", line 527, in _matmat
    return super()._matmat(X)
  File "/venv/lib/python3.9/site-packages/scipy/sparse/linalg/_interface.py", line 187, in _matmat
    return np.hstack([self.matvec(col.reshape(-1,1)) for col in X.T])
AttributeError: 'int' object has no attribute 'T'

SciPy/NumPy/Python version information

1.8.0

[lobpcg]

Is not it needed to return something meaningful in the definition of matmat so that it doesn't fall back to matvec?

[mfeuerle]

If you call aslinearoperator with a object that is not a (sparse) matrix or a LinearOperator the following code is executed:

if hasattr(A, 'shape') and hasattr(A, 'matvec'):
    rmatvec = None
    rmatmat = None
    dtype = None
    
    if hasattr(A, 'rmatvec'):
        rmatvec = A.rmatvec
    if hasattr(A, 'rmatmat'):
        rmatmat = A.rmatmat
    if hasattr(A, 'dtype'):
        dtype = A.dtype
    return LinearOperator(A.shape, A.matvec, rmatvec=rmatvec,
                          rmatmat=rmatmat, dtype=dtype)

Obviously the matmul function is completely ignored. There are two major problems with that:

1. Your object must implement the matvec function. A object that only implements matmul is not accepted even though matvec could be implemented by simply falling back to matmul as the default implementation of LinearOperator already does anyways.
2. Even if your object supplies a efficient matmul implementation it will always fall back to a probably inefficient default implementation based on matvec since the matmu function is never considered.

[lobpcg]

I do not see this example as a proof of a bug. It's documented at https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.aslinearoperator.html : "An object with .shape and .matvec attributes"

It is reasonable to request that aslinearoperator takes objects with other attributes in which case the request is for an enhancement, rather than a bug report., - just like your already opened https://github.com/scipy/scipy/issues/16114