Compatibility with numpy arrays and other matrices

[RendersJens]

Add the ability to combine linear operators with a numpy 2D arrays (or other matrices like the ones from scipy.sparse) using VStack, HStack, Block, etc. like this:

import pylops
import numpy as np

# 2x2 identity operator
A = pylops.Identity(2)

# 2x2 identity matrix
B = np.eye(2)

C = pylops.VStack([A, B])
print(C @ np.array([1, 2]))

Currently this results in
AttributeError: 'numpy.ndarray' object has no attribute 'matvec'
while a numpy 2D array does support the @ operator (as substitute for .matvec()). The same holds for other types of matrices. The implementation of VStack, HStack, Block, etc. should use the @ operator if .matvec is not implemented.

[mrava87]

Hi,
You can easily do that by wrapping your numpy/scipy array into Pylops.MatrixMult. In your case you would just new to add another line, B=Pylops.MatrixMult(B).

Let me know if that works :)

[RendersJens]

Thank you for the swift reply. That certainly works, and this is how I use it for now, but I opened this issue as I feel like this line should be avoided, as the matrix functionality (like indexing) is then lost.

It seems like it is perfectly doable to have numpy or sparse matrices interact with linear operators, by simply using their @ operator instead of .matvec. Is it a design choice to not support this? (In that case I will probably patch it for myself)

[mrava87]

Ok, I see, that makes sense. Let me take a look where the error actually arises as it may be easy to replace matvec with @ but not so easy to replace rmatvec. Basically this would require checking if the operator is really an Pylops operator or a numpy array and behave differently.

Alternatively we could do the check once at initialization and recast numpy/scipy arrays to MatrixMult so the forward and adjoins of stacking operators would remain efficient with no further checks. This would be my preferred solution, would that make sense to use that we still use MatrixMult but we avoid the user having to do it themselves?

[RendersJens]

Yes, that makes sense. Letting the functions VStack, HStack, Block, etc. do the conversions internally is a good solution.

[mrava87]

Solved here ed193ea