NotImplementedError occurred when multiply a transposed operator with a vector

[anhtuan299]

Dear PyLops developers,

I was multiplying the adjoint of <256000x262144 VStack with dtype=float32> operator W with a numpy.ndarray whose length = 256000 named res. It did work very well until the following NotImplementedError came from nowhere, even when I reinstalled Scipy, PyLops and even the Anaconda enviroment:

W.T @ res

Traceback (most recent call last): File "", line 1, in File "/home/anhtuan299/anaconda3/envs/anhtuan299/lib/python3.7/site-packages/scipy/sparse/linalg/interface.py", line 429, in matmul return self.mul(other) File "/home/anhtuan299/anaconda3/envs/anhtuan299/lib/python3.7/site-packages/pylops/LinearOperator.py", line 77, in mul y = super().mul(x) File "/home/anhtuan299/anaconda3/envs/anhtuan299/lib/python3.7/site-packages/scipy/sparse/linalg/interface.py", line 393, in mul return self.dot(x) File "/home/anhtuan299/anaconda3/envs/anhtuan299/lib/python3.7/site-packages/pylops/LinearOperator.py", line 248, in dot return self.matvec(x) File "/home/anhtuan299/anaconda3/envs/anhtuan299/lib/python3.7/site-packages/pylops/LinearOperator.py", line 129, in matvec y = self._matvec(x) File "/home/anhtuan299/anaconda3/envs/anhtuan299/lib/python3.7/site-packages/pylops/LinearOperator.py", line 56, in _matvec return self.Op._matvec(x) File "/home/anhtuan299/anaconda3/envs/anhtuan299/lib/python3.7/site-packages/scipy/sparse/linalg/interface.py", line 583, in _matvec return np.conj(self.A._rmatvec(np.conj(x))) File "/home/anhtuan299/anaconda3/envs/anhtuan299/lib/python3.7/site-packages/pylops/basicoperators/VStack.py", line 109, in _rmatvec y += oper.rmatvec(x[self.nnops[iop]:self.nnops[iop + 1]]).squeeze() File "/home/anhtuan299/anaconda3/envs/anhtuan299/lib/python3.7/site-packages/pylops/LinearOperator.py", line 162, in rmatvec y = self._rmatvec(x) File "/home/anhtuan299/anaconda3/envs/anhtuan299/lib/python3.7/site-packages/pylops/LinearOperator.py", line 60, in _rmatvec return self.Op._rmatvec(x) File "/home/anhtuan299/anaconda3/envs/anhtuan299/lib/python3.7/site-packages/scipy/sparse/linalg/interface.py", line 299, in _rmatvec raise NotImplementedError NotImplementedError

I'd love to know if you could help me to solve it. Thanks a lot.

Kind regards, Anh-Tuan

[mrava87]

Hello, Thanks for the issue.

Can I ask which version of pylops you are using and if you can paste a small code snippet reproducing your error?

[anhtuan299]

Dear,

Thanks for your response. The version I am using is 1.12.0. Below is a code that returns the error on my PC:

import numpy as np
import astra
import pylops

im_size = 512
im = np.zeros((im_size, im_size))

# CT system
vol_geom = astra.create_vol_geom(im_size, im_size)
all_angles = np.linspace(0, 2*np.pi, 500, endpoint=False)
angles = np.split(all_angles, 2)
W = []
proj_geom = astra.create_proj_geom('fanflat', 1.5, im_size, angles[0], 800, 400)
proj_id = astra.create_projector('cuda', proj_geom, vol_geom)
W0 = pylops.LinearOperator(astra.OpTomo(proj_id))
W.append(W0)
W.append(W0)
p = np.append(W[0] @ im.ravel(), W[0] @ im.ravel())
W = pylops.VStack(W)
res = W @ im.ravel() - p
gradx = W.T @ res

It requires to install a toolbox named ASTRA (conda install -c astra-toolbox astra-toolbox) to create the linear operator I am working on. Sorry I couldn't make it simpler. In fact, it works on my colleague's PC and previously worked on my PC so I suspect there must be something wrong with my PC but couldn't acknowledge it though...

[mrava87]

Hello, I think I may know why. This is a fairly old version of pylops which used to work fine with scipy at that time. However at some point scipy introduced some breaking changes to their linear operator interface which we subclass.

Could you try installing scipy<1.8.0 and see if that helps.

Alternatively, could you try installing pylops=1.18.3 with a more recent version of scipy (I think for this you need at least python 3.8 - I see you are using 3.7).

In the meantime I’ll try running your code myself and see if I can reproduce your error

[anhtuan299]

Dear,

I created a new Anaconda enviroment with all the required packages in the newest versions and the problem was solved succesfully. Thanks for your help!

Kind regards, Anh-Tuan

[mrava87]

Fantastic @anhtuan299, thanks for letting us know!

Btw, we have been wanting to have an example integrating the ASTRA toolbox in PyLops. If you are interested, it would be nice if you could make a simple example and push it to https://github.com/PyLops/pylops_notebooks/tree/master/developement :)