Skip to content

Fix: scale the operator when dividing it by a scalar #811

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mrava87 merged 2 commits into from
Oct 11, 2026
Merged

Fix: scale the operator when dividing it by a scalar #811

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
64052f2
Fix: scale the operator when dividing it by a scalar
raashish1601 Oct 7, 2026
f386b51
minor: use np.isscalar for scalar division and update type hints
raashish1601 Oct 9, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 9 additions & 5 deletions pylops/linearoperator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -835,14 +835,15 @@ def dot(self, x: NDArray) -> NDArray:

def div(
self,
y: NDArray,
y: NDArray | float,
niter: int = 100,
densesolver: str = "scipy",
) -> NDArray:
) -> NDArray | LinearOperator:
r"""Solve the linear problem :math:`\mathbf{y}=\mathbf{A}\mathbf{x}`.
Comment thread
mrava87 marked this conversation as resolved.

Overloading of operator ``/`` to improve expressivity of `Pylops`
when solving inverse problems.
when solving inverse problems. Dividing by a scalar instead scales
Comment thread
mrava87 marked this conversation as resolved.
the operator, so ``Op / a`` is the same as ``Op * (1 / a)``.

Parameters
----------
Expand All @@ -864,10 +865,13 @@ def div(

def __truediv__(
self,
y: NDArray,
y: NDArray | float,
niter: int = 100,
densesolver: str = "scipy",
) -> NDArray:
) -> NDArray | LinearOperator:
if np.isscalar(y):
# A scalar divisor scales the operator, symmetric with ``Op * a``
return self * (1 / y)
if self.explicit is True:
if sp.sparse.issparse(self.A):
# use scipy solver for sparse matrices
Expand Down
20 changes: 16 additions & 4 deletions pylops/signalprocessing/fft.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,10 @@ def _rmatvec(self, x: NDArray) -> NDArray:
y = y.astype(self.rdtype)
return y

def __truediv__(self, y: NDArray) -> NDArray:
def __truediv__(self, y: NDArray | float) -> NDArray | LinearOperator:
if np.isscalar(y):
# A scalar divisor scales the operator (see LinearOperator.div)
return super().__truediv__(y)
if self.norm is not _FFTNorms.ORTHO:
return self._rmatvec(y) / self._scale
return self._rmatvec(y)
Expand Down Expand Up @@ -238,7 +241,10 @@ def _rmatvec(self, x: NDArray) -> NDArray:
y = scipy.fft.fftshift(y, axes=self.axis)
return y

def __truediv__(self, y: NDArray) -> NDArray:
def __truediv__(self, y: NDArray | float) -> NDArray | LinearOperator:
if np.isscalar(y):
# A scalar divisor scales the operator (see LinearOperator.div)
return super().__truediv__(y)
if self.norm is not _FFTNorms.ORTHO:
return self._rmatvec(y) / self._scale
return self._rmatvec(y)
Expand Down Expand Up @@ -403,7 +409,10 @@ def _rmatvec(self, x: NDArray) -> NDArray:
y = np.real(y)
return y

def __truediv__(self, y: NDArray) -> NDArray:
def __truediv__(self, y: NDArray | float) -> NDArray | LinearOperator:
if np.isscalar(y):
# A scalar divisor scales the operator (see LinearOperator.div)
return super().__truediv__(y)
if self.norm is _FFTNorms.ORTHO:
return self._rmatvec(y)
return self._rmatvec(y) / self._scale
Expand Down Expand Up @@ -495,7 +504,10 @@ def _rmatvec(self, x: NDArray) -> NDArray:
y = scipy.fft.fftshift(y, axes=self.axis)
return y

def __truediv__(self, y: NDArray) -> NDArray:
def __truediv__(self, y: NDArray | float) -> NDArray | LinearOperator:
if np.isscalar(y):
# A scalar divisor scales the operator (see LinearOperator.div)
return super().__truediv__(y)
if self.norm is not _FFTNorms.ORTHO:
return self._rmatvec(y) / self._scale
return self._rmatvec(y)
Expand Down
15 changes: 12 additions & 3 deletions pylops/signalprocessing/fft2d.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,10 @@ def _rmatvec(self, x):
y = ncp.fft.fftshift(y, axes=self.axes[self.ifftshift_before])
return y

def __truediv__(self, y: NDArray) -> NDArray:
def __truediv__(self, y: NDArray | float) -> NDArray | LinearOperator:
if np.isscalar(y):
# A scalar divisor scales the operator (see LinearOperator.div)
return super().__truediv__(y)
if self.norm is not _FFTNorms.ORTHO:
return self._rmatvec(y) / self._scale
return self._rmatvec(y)
Expand Down Expand Up @@ -243,7 +246,10 @@ def _rmatvec(self, x):
y = scipy.fft.fftshift(y, axes=self.axes[self.ifftshift_before])
return y

def __truediv__(self, y: NDArray) -> NDArray:
def __truediv__(self, y: NDArray | float) -> NDArray | LinearOperator:
if np.isscalar(y):
# A scalar divisor scales the operator (see LinearOperator.div)
return super().__truediv__(y)
if self.norm is not _FFTNorms.ORTHO:
return self._rmatvec(y) / self._scale / self._scale
return self._rmatvec(y)
Expand Down Expand Up @@ -357,7 +363,10 @@ def _rmatvec(self, x):
y = scipy.fft.fftshift(y, axes=self.axes[self.ifftshift_before])
return y

def __truediv__(self, y: NDArray) -> NDArray:
def __truediv__(self, y: NDArray | float) -> NDArray | LinearOperator:
if np.isscalar(y):
# A scalar divisor scales the operator (see LinearOperator.div)
return super().__truediv__(y)
if self.norm is not _FFTNorms.ORTHO:
return self._rmatvec(y) / self._scale / self._scale
return self._rmatvec(y)
Expand Down
15 changes: 12 additions & 3 deletions pylops/signalprocessing/fftnd.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,7 +122,10 @@ def _rmatvec(self, x: NDArray) -> NDArray:
y = ncp.fft.fftshift(y, axes=self.axes[self.ifftshift_before])
return y

def __truediv__(self, y: NDArray) -> NDArray:
def __truediv__(self, y: NDArray | float) -> NDArray | LinearOperator:
if np.isscalar(y):
# A scalar divisor scales the operator (see LinearOperator.div)
return super().__truediv__(y)
if self.norm is not _FFTNorms.ORTHO:
return self._rmatvec(y) / self._scale
return self._rmatvec(y)
Expand Down Expand Up @@ -220,7 +223,10 @@ def _rmatvec(self, x: NDArray) -> NDArray:
y = sp_fft.fftshift(y, axes=self.axes[self.ifftshift_before])
return y

def __truediv__(self, y: NDArray) -> NDArray:
def __truediv__(self, y: NDArray | float) -> NDArray | LinearOperator:
if np.isscalar(y):
# A scalar divisor scales the operator (see LinearOperator.div)
return super().__truediv__(y)
if self.norm is not _FFTNorms.ORTHO:
return self._rmatvec(y) / self._scale
return self._rmatvec(y)
Expand Down Expand Up @@ -320,7 +326,10 @@ def _rmatvec(self, x: NDArray) -> NDArray:
y = scipy.fft.fftshift(y, axes=self.axes[self.ifftshift_before])
return y

def __truediv__(self, y: NDArray) -> NDArray:
def __truediv__(self, y: NDArray | float) -> NDArray | LinearOperator:
if np.isscalar(y):
# A scalar divisor scales the operator (see LinearOperator.div)
return super().__truediv__(y)
if self.norm is not _FFTNorms.ORTHO:
return self._rmatvec(y) / self._scale
return self._rmatvec(y)
Expand Down
26 changes: 26 additions & 0 deletions pytests/test_ffts.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -318,6 +318,32 @@ def test_unknown_engine(par):
)


@pytest.mark.parametrize("engine", ["numpy", "scipy", "fftw", "mkl_fft"])
def test_scalar_division(engine):
"""Check that dividing by a scalar scales the operator, like ``Op * (1 / a)``,
instead of applying its inverse to the scalar
"""
if engine == "mkl_fft" and not mkl_fft_enabled:
pytest.skip("mkl_fft is not installed")
if engine == "fftw" and backend == "cupy":
pytest.skip("fftw does not work with CuPy arrays")
if engine == "scipy" and backend == "cupy":
pytest.skip("scipy does not work with CuPy arrays")
Ops = [FFT(dims=(8,), nfft=8, engine=engine)]
if engine != "fftw":
# FFT2D and FFTND have no fftw engine
Ops += [
FFT2D(dims=(4, 6), engine=engine),
FFTND(dims=(2, 4, 6), engine=engine),
]
for Op in Ops:
x = np.arange(Op.shape[1], dtype=np.float64)
for scalar in (2, 2.0):
Sop = Op / scalar
assert Sop.shape == Op.shape
assert_array_almost_equal(Sop @ x, (Op @ x) / 2, decimal=10)


dtype_precision = [
(np.float16, 1),
(np.float32, 4),
Expand Down
18 changes: 18 additions & 0 deletions pytests/test_linearoperator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,24 @@ def test_scaled(par):
assert S3op.dtype == dtype


@pytest.mark.parametrize("par", [(par1), (par2), (par1j)])
def test_scalar_division(par):
"""Verify that dividing by a scalar scales the operator, like ``Op * (1 / a)``,
instead of solving ``a = Op x``
"""
shape = (par["ny"], par["nx"])
A = np.random.normal(0, 10, shape) + par["imag"] * np.random.normal(0, 10, shape)
x = np.ones(par["nx"]) + par["imag"] * np.ones(par["nx"])
Mop = MatrixMult(A, dtype=par["dtype"])
# an explicit operator and one that is not
for Op in (Mop, VStack([Mop])):
for scalar in (2, 2.0):
Sop = Op / scalar
assert isinstance(Sop, LinearOperator)
assert Sop.shape == Op.shape
assert_array_almost_equal(Sop @ x, (Op @ x) / 2, decimal=6)


@pytest.mark.parametrize("par", [(par1), (par1j)])
def test_scipyop(par):
"""Verify interaction between pylops and scipy Linear operators"""
Expand Down
Loading