feat: added ResidualNormToInitialCallback and rtol1 to solvers

Author: mrava87

docs/source/addingsolver.rst

@@ -189,8 +189,7 @@ input and returns some of the most valuable properties of the class-based solver
 
     def cg(Op, y, x0, niter=10, tol=1e-4, rtol=0.0,
            show=False, itershow=(10, 10, 10), callback=None):
-        rcallback = ResidualNormCallback(rtol)
-        cgsolve = CG(Op, callbacks=[rcallback, ])
+        cgsolve = CG(Op, callbacks=[ResidualNormToInitialCallback(rtol), ])
         if callback is not None:
             cgsolve.callback = callback
         x, iiter, cost = cgsolve.solve(

pylops/optimization/basic.py

@@ -6,7 +6,10 @@
 
 from typing import TYPE_CHECKING, Callable, Optional, Tuple
 
-from pylops.optimization.callback import ResidualNormCallback
+from pylops.optimization.callback import (
+    ResidualNormToDataCallback,
+    ResidualNormToInitialCallback,
+)
 from pylops.optimization.cls_basic import CG, CGLS, LSQR
 from pylops.utils.decorators import add_ndarray_support_to_solver
 from pylops.utils.typing import NDArray
@@ -22,7 +25,8 @@ def cg(
     x0: Optional[NDArray] = None,
     niter: int = 10,
     tol: float = 1e-4,
-    rtol: bool = 0.0,
+    rtol: float = 0.0,
+    rtol1: float = 0.0,
     show: bool = False,
     itershow: Tuple[int, int, int] = (10, 10, 10),
     callback: Optional[Callable] = None,
@@ -47,8 +51,12 @@ def cg(
         Absolute tolerance on residual norm. Stops the solver when the
         residual norm is below this value.
     rtol : :obj:`float`, optional
-        Relative tolerance on residual norm. Stops the solver when the
-        ratio of the current residual norm to the initial residual norm is
+        Relative tolerance on residual norm wrt initial residual norm. Stops
+        the solver when the ratio of the current residual norm to the initial
+        residual norm is below this value.
+    rtol1 : :obj:`float`, optional
+        Relative tolerance on residual norm wrt to data. Stops the solver
+        when the ratio of the current residual norm to the data norm is
         below this value.
     show : :obj:`bool`, optional
         Display iterations log
@@ -78,12 +86,15 @@ def cg(
     See :class:`pylops.optimization.cls_basic.CG`
 
     """
-    rcallback = ResidualNormCallback(rtol)
+    callbacks = []
+    if rtol > 0.0:
+        callbacks.append(ResidualNormToInitialCallback(rtol))
+    if rtol1 > 0.0:
+        callbacks.append(ResidualNormToDataCallback(rtol1))
+
     cgsolve = CG(
         Op,
-        callbacks=[
-            rcallback,
-        ],
+        callbacks=callbacks if len(callbacks) > 0 else None,
     )
     if callback is not None:
         cgsolve.callback = callback
@@ -108,6 +119,7 @@ def cgls(
     damp: float = 0.0,
     tol: float = 1e-4,
     rtol: float = 0.0,
+    rtol1: float = 0.0,
     show: bool = False,
     itershow: Tuple[int, int, int] = (10, 10, 10),
     callback: Optional[Callable] = None,
@@ -134,8 +146,12 @@ def cgls(
         Absolute tolerance on residual norm. Stops the solver when the
         residual norm is below this value.
     rtol : :obj:`float`, optional
-        Relative tolerance on residual norm. Stops the solver when the
-        ratio of the current residual norm to the initial residual norm is
+        Relative tolerance on residual norm wrt initial residual norm. Stops
+        the solver when the ratio of the current residual norm to the initial
+        residual norm is below this value.
+    rtol1 : :obj:`float`, optional
+        Relative tolerance on residual norm wrt to data. Stops the solver
+        when the ratio of the current residual norm to the data norm is
         below this value.
     show : :obj:`bool`, optional
         Display iterations log
@@ -180,12 +196,14 @@ def cgls(
     See :class:`pylops.optimization.cls_basic.CGLS`
 
     """
-    rcallback = ResidualNormCallback(rtol)
+    callbacks = []
+    if rtol > 0.0:
+        callbacks.append(ResidualNormToInitialCallback(rtol))
+    if rtol1 > 0.0:
+        callbacks.append(ResidualNormToDataCallback(rtol1))
     cgsolve = CGLS(
         Op,
-        callbacks=[
-            rcallback,
-        ],
+        callbacks=callbacks if len(callbacks) > 0 else None,
     )
     if callback is not None:
         cgsolve.callback = callback

pylops/optimization/callback.py

@@ -1,11 +1,14 @@
 __all__ = [
     "Callbacks",
     "MetricsCallback",
-    "ResidualNormCallback",
+    "ResidualNormToDataCallback",
+    "ResidualNormToInitialCallback",
 ]
 
 from typing import TYPE_CHECKING, Dict, List, Optional, Sequence
 
+import numpy as np
+
 from pylops.utils.metrics import mae, mse, psnr, snr
 from pylops.utils.typing import NDArray
 
@@ -191,8 +194,37 @@ def on_step_end(self, solver: "Solver", x: NDArray) -> None:
             self.metrics["psnr"].append(psnr(self.xtrue, x))
 
 
-class ResidualNormCallback(Callbacks):
-    """Residual norm callback
+class ResidualNormToDataCallback(Callbacks):
+    """Residual norm to data callback
+
+    This callback can be used to stop the solver when the residual norm
+    is below a certain threshold defined as a percentage of the
+    initial residual norm.
+
+    Parameters
+    ----------
+    rtol : :obj:`float`
+        Percentage of the initial residual norm below which the solver
+        will stop iterating. For example, if `rtol` is 0.1, the solver
+        will stop when the residual norm is below 10% of the initial
+        residual norm.
+
+    """
+
+    def __init__(self, rtol: float) -> None:
+        self.rtol = rtol
+        self.stop = False
+
+    def on_setup_end(self, solver: "Solver", x: NDArray) -> None:
+        self.ynorm = self.ncp.linalg.norm(self.y)
+
+    def on_step_end(self, solver: "Solver", x: NDArray) -> None:
+        if solver.cost[-1] < self.rtol * self.ynorm:
+            self.stop = True
+
+
+class ResidualNormToInitialCallback(Callbacks):
+    """Residual norm to initial callback
 
     This callback can be used to stop the solver when the residual norm
     is below a certain threshold defined as a percentage of the

pylops/optimization/sparsity.py

@@ -9,7 +9,10 @@
 
 from typing import TYPE_CHECKING, Any, Callable, Dict, List, Optional, Tuple
 
-from pylops.optimization.callback import ResidualNormCallback
+from pylops.optimization.callback import (
+    ResidualNormToDataCallback,
+    ResidualNormToInitialCallback,
+)
 from pylops.optimization.cls_sparsity import FISTA, IRLS, ISTA, OMP, SPGL1, SplitBregman
 from pylops.utils.decorators import add_ndarray_support_to_solver
 from pylops.utils.typing import NDArray, SamplingLike
@@ -142,6 +145,7 @@ def omp(
     niter_inner: int = 40,
     sigma: float = 1e-4,
     rtol: float = 0.0,
+    rtol1: float = 0.0,
     normalizecols: bool = False,
     Opbasis: Optional["LinearOperator"] = None,
     optimal_coeff: bool = False,
@@ -172,9 +176,13 @@ def omp(
     sigma : :obj:`float`, optional
         Maximum :math:`L_2` norm of residual. When smaller stop iterations.
     rtol : :obj:`float`, optional
-        Relative tolerance on residual. Stops the solver when the
-        ratio of the current residual norm to the initial residual norm
-        is below this value.
+        Relative tolerance on residual norm wrt initial residual norm. Stops
+        the solver when the ratio of the current residual norm to the initial
+        residual norm is below this value.
+    rtol1 : :obj:`float`, optional
+        Relative tolerance on residual norm wrt to data. Stops the solver
+        when the ratio of the current residual norm to the data norm is
+        below this value.
     normalizecols : :obj:`list`, optional
         Normalize columns (``True``) or not (``False``). Note that this can be
         expensive as it requires applying the forward operator
@@ -229,12 +237,15 @@ def omp(
     See :class:`pylops.optimization.cls_sparsity.OMP`
 
     """
-    rcallback = ResidualNormCallback(rtol)
+    callbacks = []
+    if rtol > 0.0:
+        callbacks.append(ResidualNormToInitialCallback(rtol))
+    if rtol1 > 0.0:
+        callbacks.append(ResidualNormToDataCallback(rtol1))
+
     ompsolve = OMP(
         Op,
-        callbacks=[
-            rcallback,
-        ],
+        callbacks=callbacks if len(callbacks) > 0 else None,
     )
     if callback is not None:
         ompsolve.callback = callback
@@ -264,7 +275,8 @@ def ista(
     alpha: Optional[float] = None,
     eigsdict: Optional[Dict[str, Any]] = None,
     tol: float = 1e-10,
-    rtol: bool = 0.0,
+    rtol: float = 0.0,
+    rtol1: float = 0.0,
     threshkind: str = "soft",
     perc: Optional[float] = None,
     decay: Optional[NDArray] = None,
@@ -309,9 +321,13 @@ def ista(
         Absolute tolerance on model update. Stop iterations if difference between inverted model
         at subsequent iterations is smaller than ``tol``
     rtol : :obj:`float`, optional
-        Relative tolerance on total cost function. Stops the solver when the
-        ratio of the current cost function to the initial cost function
-        is below this value.
+        Relative tolerance on residual norm wrt initial residual norm. Stops
+        the solver when the ratio of the current residual norm to the initial
+        residual norm is below this value.
+    rtol1 : :obj:`float`, optional
+        Relative tolerance on residual norm wrt to data. Stops the solver
+        when the ratio of the current residual norm to the data norm is
+        below this value.
     threshkind : :obj:`str`, optional
         Kind of thresholding ('hard', 'soft', 'half', 'hard-percentile',
         'soft-percentile', or 'half-percentile' - 'soft' used as default)
@@ -370,12 +386,15 @@ def ista(
     See :class:`pylops.optimization.cls_sparsity.ISTA`
 
     """
-    rcallback = ResidualNormCallback(rtol)
+    callbacks = []
+    if rtol > 0.0:
+        callbacks.append(ResidualNormToInitialCallback(rtol))
+    if rtol1 > 0.0:
+        callbacks.append(ResidualNormToDataCallback(rtol1))
+
     istasolve = ISTA(
         Op,
-        callbacks=[
-            rcallback,
-        ],
+        callbacks=callbacks if len(callbacks) > 0 else None,
     )
     if callback is not None:
         istasolve.callback = callback
@@ -410,6 +429,7 @@ def fista(
     eigsdict: Optional[Dict[str, Any]] = None,
     tol: float = 1e-10,
     rtol: float = 0.0,
+    rtol1: float = 0.0,
     threshkind: str = "soft",
     perc: Optional[float] = None,
     decay: Optional[NDArray] = None,
@@ -513,12 +533,15 @@ def fista(
     See :class:`pylops.optimization.cls_sparsity.FISTA`
 
     """
-    rcallback = ResidualNormCallback(rtol)
+    callbacks = []
+    if rtol > 0.0:
+        callbacks.append(ResidualNormToInitialCallback(rtol))
+    if rtol1 > 0.0:
+        callbacks.append(ResidualNormToDataCallback(rtol1))
+
     fistasolve = FISTA(
         Op,
-        callbacks=[
-            rcallback,
-        ],
+        callbacks=callbacks if len(callbacks) > 0 else None,
     )
     if callback is not None:
         fistasolve.callback = callback
@@ -673,6 +696,7 @@ def splitbregman(
     epsRL2s: Optional[SamplingLike] = None,
     tol: float = 1e-10,
     rtol: float = 0.0,
+    rtol1: float = 0.0,
     tau: float = 1.0,
     restart: bool = False,
     engine: str = "scipy",
@@ -726,9 +750,13 @@ def splitbregman(
         Tolerance. Stop the solver if difference between inverted model
         at subsequent iterations is smaller than ``tol``
     rtol : :obj:`float`, optional
-        Relative tolerance on total cost function. Stops the solver when the
-        ratio of the current cost function to the initial cost function
-        is below this value.
+        Relative tolerance on residual norm wrt initial residual norm. Stops
+        the solver when the ratio of the current residual norm to the initial
+        residual norm is below this value.
+    rtol1 : :obj:`float`, optional
+        Relative tolerance on residual norm wrt to data. Stops the solver
+        when the ratio of the current residual norm to the data norm is
+        below this value.
     tau : :obj:`float`, optional
         Scaling factor in the Bregman update (must be close to 1)
     restart : :obj:`bool`, optional
@@ -772,12 +800,14 @@ def splitbregman(
     See :class:`pylops.optimization.cls_sparsity.SplitBregman`
 
     """
-    rcallback = ResidualNormCallback(rtol)
+    callbacks = []
+    if rtol > 0.0:
+        callbacks.append(ResidualNormToInitialCallback(rtol))
+    if rtol1 > 0.0:
+        callbacks.append(ResidualNormToDataCallback(rtol1))
     sbsolve = SplitBregman(
         Op,
-        callbacks=[
-            rcallback,
-        ],
+        callbacks=callbacks if len(callbacks) > 0 else None,
     )
     if callback is not None:
         sbsolve.callback = callback

pytests/test_sparsity.py

@@ -14,7 +14,7 @@
 import pytest
 
 from pylops.basicoperators import FirstDerivative, Identity, MatrixMult
-from pylops.optimization.callback import ResidualNormCallback
+from pylops.optimization.callback import ResidualNormToInitialCallback
 from pylops.optimization.cls_sparsity import IRLS
 from pylops.optimization.sparsity import fista, irls, ista, omp, spgl1, splitbregman
 
@@ -226,7 +226,7 @@ def test_IRLS_model_stopping(par):
     rtol = 6e-1
     kwars_solver = dict(iter_lim=5) if backend == "numpy" else dict(niter=5)
 
-    rcallback = ResidualNormCallback(rtol)
+    rcallback = ResidualNormToInitialCallback(rtol)
     irlssolve = IRLS(
         Aop,
         callbacks=[