[SciPy-Dev] Implementation choices in scipy.sparse.linalg.LinearOperator

[matteo ravasi]

Dear All,
just trying to revive this as I didn’t get any reply.

I can see that there is more interest in the community with respect to ad-hoc matrix-vector multiplications (aka Linear operators) - see also the discussion in Re: [SciPy-Dev] Efficient Toeplitz Matrix-Vector Multiplication and related PR.

Hopefully some of the core developers or anyone involved in this area can comment here ;)

Thank you!
MR



> On 8 Jan 2020, at 20:31, matteo ravasi <matteoravasi at gmail.com> wrote:
> 
> Dear Scipy Developers,
> some of you may be aware of the PyLops project (https://github.com/equinor/pylops <https://github.com/equinor/pylops>), which heavily builds on the scipy.sparse.linalg.LinearOperator class of scipy and provides an extensive library of general purpose and domain specific linear operators.
> 
> Recently a users asked for a functionality that would create the equivalent dense matrix representation of an operator (to be mostly used for teaching purposes). While implementing it, we came across two internal behaviours in the interface.py file that we would like to mention and discuss about:
> 
> - The method _matmat and _rmatmat as implemented per today requires operators to work on (N, 1)/(M, 1) arrays. Whilst the first sentence in the Notes of LinearOperator states "The user-defined matvec() function must properly handle the case where v has shape (N,) as well as the (N,1) case", I tend to believe that the first option should be the encouraged one, as I do not really see any advantage with carrying over trailing dimensions in this case. Moreover, some of the linear solvers in scipy that can be used together with linear operator such as lsqr or lsmr, require the RHS b vector to have a single dimension (b : (m,) ndarray). Again, although a quick test showed that they work also when providing a b vector of shape (m, 1), this shows to me that the extra trailing dimension is just a nuisance and should be discouraged. Assuming that you agree with this argument, an alternative solution that works on (N, ) input is actually possible and currently implemented in PyLops for _matmat (see here https://github.com/equinor/pylops/blob/master/pylops/LinearOperator.py <https://github.com/equinor/pylops/blob/master/pylops/LinearOperator.py>). This avoids having implementations that work for both (N, ) and (N, 1), which sometimes requires adding additional checks and squeezes within the implementation of _matvec and _rmatvec in situations where trailing dimensions in the input vector are not really bringing any benefit.
> 
> - The above mentioned new feature has been currently added as a method to the PyLops LinearOperator class (and named todense in a similar fashion to scipy.sparse.xxx_matrix.todense methods). Doing so, we realized that every time two operators are summed/subtracted/multiplied… (pretty much when one of the overloading of scipy.sparse.linalg.LinearOperator is invoked), the returned object is not of LinearOperator type anymore but of the type of the private class in scipy.sparse.linalg.interface used to perform the requested operator (e.g, _SumLinearOperator). I wonder if returning an object of a ‘private’ type is actually something recommended to do in scipy (or even if this happens anywhere else in scipy). Again, in PyLops we need to overload the __add__ method (and similar methods) and make it into a thin wrappers of the scipy equivalent, returning an object of our parent class LinearOperator. This is to ensure that the returned object would maintain the additional properties and methods of our Linearoperator class, but also to avoid ‘leaking’ any of the private classes to the end users. I wonder if also the scipy API should perhaps avoid returning to end users objects of some private classes.
> 
> This email is not intended to criticise any of current behaviours in scipy.sparse.linalg.LinearOperator, rather to gain an understanding of some of the decisions made by the core developers of the scipy.sparse.linalg.interface submodule, which may all have their own right to be the way they are today :)
> 
> Looking forward to hearing from you.
> 
> Best wishes,
> MR

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20200124/2902d29e/attachment.html>