diff --git a/examples/plot_tvreg.py b/examples/plot_tvreg.py index 30692eb0..352f9cac 100755 --- a/examples/plot_tvreg.py +++ b/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)) diff --git a/pylops/optimization/cls_sparsity.py b/pylops/optimization/cls_sparsity.py index 33eff99b..ac91c85a 100644 --- a/pylops/optimization/cls_sparsity.py +++ b/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 diff --git a/pylops/optimization/sparsity.py b/pylops/optimization/sparsity.py index 3cbf6fce..611e102b 100644 --- a/pylops/optimization/sparsity.py +++ b/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 diff --git a/pytests/test_sparsity.py b/pytests/test_sparsity.py index 7d4e61dd..0e6db86b 100644 --- a/pytests/test_sparsity.py +++ b/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) diff --git a/pytests/test_swts.py b/pytests/test_swts.py index c2c8c00e..d1a1833c 100644 --- a/pytests/test_swts.py +++ b/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"])