[SciPy-Dev] extended scipy.sparse.linalg LinearOperator interface

[Lars Buitinck]

Hi all,


[TL;DR I want to extend scipy.sparse.linalg.LinearOperator and I'm
looking for current usage patterns.]


To solve a few problems that I and others [1,2] encountered with the
scipy.sparse.linalg LinearOperator interface, I decided to expand and
refactor it. In a nutshell:

* Linear operators cannot be transposed. They do have an rmatvec that
implements A^H * x, but no matrix-matrix multiplication version of
same, so this has be implemented as a loop.
* While a lot of subclasses exist in various parts of
scipy.sparse.linalg, there was no documentation on how to roll your
own operator by subclassing. Instead you have to call the constructor
with a few functions, which become the methods on the custom operator.
This doesn't scale if we want to add more methods.
* The current implementation uses monkey-patching, causing memory
leaks due to reference cycles and making the code hard to read.

I've already submitted an early PR [3] with the main parts of my proposal:

* LinearOperator is now an abstract base class. An overloaded __new__
makes sure you can still call the constructor with the old calling
conventions; it returns a subclass.
* Subclassing is possible (and documented): you must provide a method
_matvec and optionally a few more. These get used to implement the
public matvec method, which uses _matvec but adds input validation
boilerplate (the "template method pattern").
* An "adjoint" method is added, and can be overridden by supplying an
_adjoint method, to return the Hermitian adjoint (conjugate transpose)
of a linear operator. This obviates the need for an rmatmat method.

There's been discussion lately about merging in parts of PyOperators.
I decided not to wait for this to happen, but I did have a look at
their documentation and used some API conventions to make a later
merge easier (or at least be compatible).

I'd like feedback on this, and especially on the issue of backward
compatibility with the current code. In particular, I'm not sure if
the possibility to subclass LinearOperator has ever been advertised
and/or exploited in third-party code. If it was, then we have to deal
with that.

To be sure, the currently common idioms to construct a linear operator,

    A = aslinearoperator(X)  # ndarray or sparse matrix X

and

    A = LinearOperator(matvec=some_callable, shape=some_tuple)

still work and all current tests still pass.

So I'd like to ask you:

* Is my proposal useful?
* How was this API used in third-party code?


TIA,
Lars


[1] https://github.com/scipy/scipy/pull/4073
[2] https://github.com/scipy/scipy/pull/4074
[3] https://github.com/scipy/scipy/pull/4087