Fix: separate Split Bregman inner solver options

examples/plot_tvreg.py

@@ -75,7 +75,9 @@
 plt.tight_layout()
 
 ###############################################################################
-# Now we impose blockiness in the solution using the Split Bregman solver
+# Now we impose blockiness in the solution using the Split Bregman solver.
+# Note that ``kwargs_solver`` can be used to supply parameters to the inner
+# solver.
 Dop = pylops.FirstDerivative(nx, edge=True, kind="backward")
 mu = 0.01
 lamda = 0.3
@@ -92,7 +94,7 @@
     epsRL1s=[lamda],
     tol=1e-4,
     tau=1.0,
-    **dict(iter_lim=30, damp=1e-10)
+    kwargs_solver=dict(iter_lim=30, damp=1e-10),
 )[0]
 
 plt.figure(figsize=(10, 5))
@@ -175,7 +177,7 @@
     tol=1e-4,
     tau=1.0,
     show=False,
-    **dict(iter_lim=5, damp=1e-4)
+    **dict(iter_lim=5, damp=1e-4),
 )[0]
 xinv = np.real(xinv.reshape(ny, nx))
 

pylops/optimization/cls_sparsity.py

@@ -3355,7 +3355,8 @@ def solve(
         show: bool = False,
         itershow: tuple[int, int, int] = (10, 10, 10),
         show_inner: bool = False,
-        **kwargs_lsqr,
+        kwargs_solver: dict | None = None,
+        **kwargs,
     ) -> tuple[NDArray, int, NDArray]:
         r"""Run entire solver
 
@@ -3428,11 +3429,19 @@ def solve(
             three element of the list.
         show_inner : :obj:`bool`, optional
             Display iteration logs of the solver of the :math:`\mathbf{x}`-subproblem
-        **kwargs_lsqr
-            Arbitrary keyword arguments for the solver of the
-            :math:`\mathbf{x}`-subproblem of the Split Bregman algorithm
-            (:py:func:`scipy.sparse.linalg.lsqr` or
+        kwargs_solver : :obj:`dict`, optional
+            .. versionadded:: 2.9.0
+
+            Keyword arguments for the solver of the
+            :math:`\mathbf{x}`-subproblem (:py:func:`scipy.sparse.linalg.lsqr` or
             :py:func:`pylops.optimization.solver.cgls`, see ``engine``).
+            Entries override duplicate values in ``**kwargs``.
+        **kwargs
+            *Deprecated*, will be removed in v3.0.0. Simply kept for
+            back-compatibility with previous implementation.
+
+            Another way to pass keyword arguments to the solver of the
+            :math:`\mathbf{x}`-subproblem.
 
         Returns
         -------
@@ -3462,13 +3471,14 @@ def solve(
             preallocate=preallocate,
             show=show,
         )
+        kwargs_solver = {**kwargs, **(kwargs_solver or {})}
         x = self.run(
             x,
             engine=engine,
             show=show,
             itershow=itershow,
             show_inner=show_inner,
-            **kwargs_lsqr,
+            **kwargs_solver,
         )
         self.finalize(show)
         return x, self.iiter, self.cost

pylops/optimization/sparsity.py

@@ -727,7 +727,8 @@ def splitbregman(
     show_inner: bool = False,
     callback: Callable | None = None,
     preallocate: bool = False,
-    **kwargs_lsqr,
+    kwargs_solver: dict | None = None,
+    **kwargs,
 ) -> tuple[NDArray, int, NDArray]:
     r"""Split Bregman for mixed L2-L1 norms.
 
@@ -814,11 +815,19 @@ def splitbregman(
         Pre-allocate all variables used by the solver. Note that if ``y``
         is a JAX array, this option is ignored and variables are not
         pre-allocated since JAX does not support in-place operations.
-    **kwargs_lsqr
-        Arbitrary keyword arguments for the solver of the
-        :math:`\mathbf{x}`-subproblem of the Split Bregman algorithm
-        (:py:func:`scipy.sparse.linalg.lsqr` or
+    kwargs_solver : :obj:`dict`, optional
+        .. versionadded:: 2.9.0
+
+        Keyword arguments for the solver of the
+        :math:`\mathbf{x}`-subproblem (:py:func:`scipy.sparse.linalg.lsqr` or
         :py:func:`pylops.optimization.solver.cgls`, see ``engine``).
+        Entries override duplicate values in ``**kwargs``.
+    **kwargs
+        *Deprecated*, will be removed in v3.0.0. Simply kept for
+        back-compatibility with previous implementation.
+
+        Another way to pass keyword arguments to the solver of the
+        :math:`\mathbf{x}`-subproblem.
 
     Returns
     -------
@@ -867,6 +876,7 @@ def splitbregman(
         show=show,
         itershow=itershow,
         show_inner=show_inner,
-        **kwargs_lsqr,
+        kwargs_solver=kwargs_solver,
+        **kwargs,
     )
     return xinv, itn_out, cost

pytests/test_sparsity.py

@@ -737,3 +737,57 @@ def test_SplitBregman_lambdaRL1s(par):
         Iop, y, [Dop], lambdaRL1s=[epsRL1], **kwargs_sb, **kwars_solver
     )[0]
     assert_array_almost_equal(xdef, xeps, decimal=12)
+
+
+def test_SplitBregman_kwargs_solver_cgls():
+    """Check that CGLS parameters can be passed through kwargs_solver"""
+    np.random.seed(42)
+    nx = 10
+    Iop = Identity(nx)
+    Dop = FirstDerivative(nx, edge=True)
+    y = np.random.normal(size=nx)
+
+    # Compute expected solution with dense solver
+    derivative = Dop.todense(backend=backend)
+    expected = np.linalg.solve(np.eye(nx) + derivative.T @ derivative, y)
+
+    # Computed solution with SplitBregman using CGLS with tight inner tolerance
+    # to ensure full convergence to the expected solution.
+    result, _, _ = splitbregman(
+        Iop,
+        y,
+        [Dop],
+        niter_outer=1,
+        niter_inner=1,
+        mu=1.0,
+        epsRL1s=[1.0],
+        tol=1e-10,
+        engine="pylops",
+        damp=0,
+        kwargs_solver=dict(niter=100, tol=1e-14, rtol=0.0, rtol1=0.0),
+    )
+    assert_array_almost_equal(result, expected, decimal=6)
+
+
+@pytest.mark.parametrize("rtol,rtol1", [(0.99, 0.0), (0.0, 0.99)])
+def test_SplitBregman_inner_relative_tolerance(rtol, rtol1):
+    """Check that rtol/rtol1 do not leak into the inner CGLS solver."""
+    np.random.seed(42)
+    nx = 10
+    Iop = Identity(nx)
+    Dop = FirstDerivative(nx, edge=True)
+    y = np.linspace(-1.0, 1.0, nx)
+
+    # CGLS rtol/rtol1 are left at their defaults: if the outer ones leaked,
+    # CGLS would stop early and the solutions below would differ
+    kwargs = dict(
+        niter_outer=1,
+        niter_inner=1,
+        epsRL1s=[0.2],
+        tol=0.0,
+        engine="pylops",
+        kwargs_solver=dict(niter=20, tol=1e-14),
+    )
+    xref, _, _ = splitbregman(Iop, y, [Dop], **kwargs)
+    x, _, _ = splitbregman(Iop, y, [Dop], rtol=rtol, rtol1=rtol1, **kwargs)
+    assert_array_almost_equal(x, xref, decimal=10)

pytests/test_swts.py

@@ -78,6 +78,7 @@ def test_wrong_level(par):
 @pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4)])
 def test_SWT_1dsignal(par):
     """Dot-test and inversion for SWT operator for 1d signal"""
+    np.random.seed(42)
     dtype = np.empty(0, dtype=par["dtype"]).real.dtype
 
     SWTop = SWT(dims=[par["nt"]], axis=0, wavelet="haar", level=3, dtype=par["dtype"])
@@ -110,6 +111,7 @@ def test_SWT_1dsignal(par):
 @pytest.mark.parametrize("axis", [0, 1])
 def test_SWT_2dsignal(par, axis):
     """Dot-test and inversion for SWT operator for 2d signal"""
+    np.random.seed(42)
     dtype = np.empty(0, dtype=par["dtype"]).real.dtype
 
     SWTop = SWT(
@@ -148,6 +150,7 @@ def test_SWT_2dsignal(par, axis):
 @pytest.mark.parametrize("axis", [0, 1, 2])
 def test_SWT_3dsignal(par, axis):
     """Dot-test and inversion for SWT operator for 3d signal"""
+    np.random.seed(42)
     dtype = np.empty(0, dtype=par["dtype"]).real.dtype
     SWTop = SWT(
         dims=(par["nt"], par["nx"], par["ny"]),
@@ -186,6 +189,7 @@ def test_SWT_3dsignal(par, axis):
 @pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4)])
 def test_SWT2D_2dsignal(par):
     """Dot-test and inversion for SWT2D operator for 2d signal"""
+    np.random.seed(42)
     dtype = np.empty(0, dtype=par["dtype"]).real.dtype
 
     SWTop = SWT2D(
@@ -224,6 +228,7 @@ def test_SWT2D_2dsignal(par):
 @pytest.mark.parametrize("axes", [(0, 1), (0, 2), (1, 2)])
 def test_SWT2D_3dsignal(par, axes):
     """Dot-test and inversion for SWT operator for 3d signal"""
+    np.random.seed(42)
     dtype = np.empty(0, dtype=par["dtype"]).real.dtype
 
     SWTop = SWT2D(
@@ -263,6 +268,7 @@ def test_SWT2D_3dsignal(par, axes):
 @pytest.mark.parametrize("par", [(par1), (par2), (par3), (par4)])
 def test_SWTND_3dsignal(par):
     """Dot-test and inversion for SWTND operator for 3d signal"""
+    np.random.seed(42)
     dtype = np.empty(0, dtype=par["dtype"]).real.dtype
 
     SWTop = SWTND(
@@ -305,6 +311,7 @@ def test_SWTND_3dsignal(par):
 )
 def test_SWTND_4dsignal(par, axes):
     """Dot-test and inversion for SWTND operator for 4d signal"""
+    np.random.seed(42)
     dtype = np.empty(0, dtype=par["dtype"]).real.dtype
 
     SWTop = SWTND(
@@ -346,6 +353,7 @@ def test_SWTND_4dsignal(par, axes):
 @pytest.mark.parametrize("kind", ["SWT", "SWT2D", "SWTND"])
 def test_SWTs_wavelets(par, wavelet, kind):
     """Dot-test for SWT, SWT2D and SWTND operators with non-haar wavelets"""
+    np.random.seed(42)
     dtype = np.empty(0, dtype=par["dtype"]).real.dtype
     dims = (par["nt"], par["nx"], par["ny"])