diff --git a/pylops/optimization/cls_sparsity.py b/pylops/optimization/cls_sparsity.py index ce30d3be..3f14cc13 100644 --- a/pylops/optimization/cls_sparsity.py +++ b/pylops/optimization/cls_sparsity.py @@ -2758,6 +2758,9 @@ class SplitBregman(Solver): List of L1 and L2 regularization terms. epsRs : :obj:`list` List of L1 and L2 regularization dampings. + threshRL1s : :obj:`list` + List of thresholds of the shrinkage steps (one per :math:`L_1` + regularization term). cost : :obj:`numpy.ndarray`, optional History of total cost function through iterations. iiter : :obj:`int` @@ -2782,11 +2785,8 @@ class SplitBregman(Solver): are the damping factors used to weight the different :math:`L_2` regularization terms of the cost function and :math:`\epsilon_{\mathbf{R}_{1,i}}` are the damping factors of the different :math:`L_1` regularization - terms of the cost function. Note that :math:`\epsilon_{\mathbf{R}_{1,i}}` is used - both as weight of the augmented :math:`L_2` term and as threshold of the - shrinkage step of the Split-Bregman algorithm (see below): as a consequence, - the effective weight of each :math:`L_1` regularization term is - :math:`\epsilon_{\mathbf{R}_{1,i}}^2`. + terms of the cost function. Note that the effective weight of each + :math:`L_1` regularization term is :math:`\epsilon_{\mathbf{R}_{1,i}}^2`. The generalized Split-Bregman algorithm [1]_ is used to solve such cost function: the algorithm is composed of a sequence of unconstrained @@ -2810,7 +2810,7 @@ class SplitBregman(Solver): \; & \frac{\mu}{2} \|\textbf{y} - \textbf{Op}\,\textbf{x}\|_2^2 \\ & + \frac{1}{2}\sum_i \epsilon_{\mathbf{R}_{2,i}} \|\mathbf{y}_{\mathbf{R}_{2,i}} - \mathbf{R}_{2,i} \textbf{x}\|_2^2 \\ - & + \frac{1}{2}\sum_i \epsilon_{\mathbf{R}_{1,i}} \|\textbf{d}_i - + & + \frac{1}{2}\sum_i \lambda_{\mathbf{R}_{1,i}} \|\textbf{d}_i - \mathbf{R}_{1,i} \textbf{x} - \textbf{b}_i^k\|_2^2 \\ & + \sum_i \epsilon_{\mathbf{R}_{1,i}}^2 \| \textbf{d}_i \|_1 \end{aligned} @@ -2820,8 +2820,14 @@ class SplitBregman(Solver): \tau (\mathbf{R}_{1,i} \textbf{x}^{k+1} - \textbf{d}_i^{k+1}) where :math:`\textbf{d}_i` are the split variables, :math:`\textbf{b}_i` - are the Bregman variables, and :math:`\tau` is a scaling factor of the - Bregman update (:math:`\tau=1` in the original algorithm [1]_). + are the Bregman variables, :math:`\lambda_{\mathbf{R}_{1,i}}` are the + weights of the splitting terms, and :math:`\tau` is a scaling factor of the + Bregman update (:math:`\tau=1` in the original algorithm [1]_). The + weights :math:`\lambda_{\mathbf{R}_{1,i}}` do not change the solution of + the problem, but affect the conditioning of the + :math:`\textbf{x}`-subproblem and the convergence speed of the algorithm. + When not provided, :math:`\lambda_{\mathbf{R}_{1,i}} = + \epsilon_{\mathbf{R}_{1,i}}`. The first step is solved by alternating ``niter_inner`` times the minimization over :math:`\textbf{x}` and over :math:`\textbf{d}_i`. The @@ -2832,7 +2838,8 @@ class SplitBregman(Solver): any other array type, e.g., CuPy or JAX arrays). The :math:`\textbf{d}_i`-subproblems are solved in closed form by soft thresholding :math:`\mathbf{R}_{1,i} \textbf{x}^{k+1} + \textbf{b}_i^k` - with threshold :math:`\epsilon_{\mathbf{R}_{1,i}}`. The entire + with threshold :math:`\epsilon_{\mathbf{R}_{1,i}}^2 / + \lambda_{\mathbf{R}_{1,i}}`. The entire procedure is repeated ``niter_outer`` times, or until the norm of the difference between the models of two subsequent outer iterations is smaller than ``tol``. @@ -2848,7 +2855,8 @@ def _print_setup(self, xcomplex: bool = False) -> None: strpar = ( f"niter_outer = {self.niter_outer:3d} niter_inner = {self.niter_inner:3d} tol = {self.tol:2.2e}\n" - f"mu = {self.mu:2.2e} epsL1 = {self.epsRL1s}\t epsL2 = {self.epsRL2s}" + f"mu = {self.mu:2.2e} epsL1 = {self.epsRL1s}\t epsL2 = {self.epsRL2s}\n" + f"lambdaL1 = {self.lambdaRL1s}" ) print(strpar) print("-" * 65) @@ -2930,6 +2938,7 @@ def setup( mu: float = 1.0, epsRL1s: SamplingLike | None = None, epsRL2s: SamplingLike | None = None, + lambdaRL1s: SamplingLike | None = None, tol: float = 1e-10, tau: float = 1.0, restart: bool = False, @@ -2962,14 +2971,21 @@ def setup( of ``RegsL2`` or equal to ``None`` to use a zero data for every regularization operator in ``RegsL2``) mu : :obj:`float`, optional - Data term damping + Data term damping epsRL1s : :obj:`list` - :math:`L_1` Regularization dampings (must have the same number of elements - as ``RegsL1``). Note that the effective weight of each :math:`L_1` - regularization term in the cost function is ``epsRL1s[i]**2`` + :math:`L_1` Regularization dampings (must have the same number of elements + as ``RegsL1``). Note that the effective weight of each :math:`L_1` + regularization term in the cost function is ``epsRL1s[i]**2`` epsRL2s : :obj:`list` - :math:`L_2` Regularization dampings (must have the same number of elements - as ``RegsL2``) + :math:`L_2` Regularization dampings (must have the same number of elements + as ``RegsL2``) + lambdaRL1s : :obj:`list`, optional + .. versionadded:: 2.9.0 + + Weights of the splitting terms (must have the same number of + elements as ``RegsL1``). They do not change the solution of the + problem, but affect the convergence speed of the algorithm. If + ``None``, they are set equal to ``epsRL1s`` tol : :obj:`float`, optional Tolerance. Stop outer iterations if difference between inverted model at subsequent iterations is smaller than ``tol`` @@ -3004,6 +3020,7 @@ def setup( self.mu = mu self.epsRL1s = list(epsRL1s) if epsRL1s is not None else [] self.epsRL2s = list(epsRL2s) if epsRL2s is not None else [] + self.lambdaRL1s = list(lambdaRL1s) if lambdaRL1s is not None else self.epsRL1s self.tol = tol self.tau = tau self.restart = restart @@ -3014,6 +3031,20 @@ def setup( # L1 regularizations self.nregsL1 = len(RegsL1) + if lambdaRL1s is not None: + if len(self.lambdaRL1s) != self.nregsL1: + msg = ( + f"lambdaRL1s must have the same number of elements as RegsL1 " + f"({len(self.lambdaRL1s)} != {self.nregsL1})" + ) + raise ValueError(msg) + if any(lambdaRL1 <= 0 for lambdaRL1 in self.lambdaRL1s): + msg = "lambdaRL1s must be strictly positive" + raise ValueError(msg) + self.threshRL1s = [ + epsRL1 * (epsRL1 / lambdaRL1) + for epsRL1, lambdaRL1 in zip(self.epsRL1s, self.lambdaRL1s, strict=True) + ] self.b = [ self.ncp.zeros(RegL1.shape[0], dtype=self.Op.dtype) for RegL1 in RegsL1 ] @@ -3039,10 +3070,9 @@ def setup( self.epsRs += [ sqrt(epsRL2s[ireg] / 2) / sqrt(mu / 2) for ireg in range(self.nregsL2) ] - if epsRL1s is not None: - self.epsRs += [ - sqrt(epsRL1s[ireg] / 2) / sqrt(mu / 2) for ireg in range(self.nregsL1) - ] + self.epsRs += [ + sqrt(lambdaRL1 / 2) / sqrt(mu / 2) for lambdaRL1 in self.lambdaRL1s + ] self.x0 = x0 x = self.ncp.zeros(self.Op.shape[1], dtype=self.Op.dtype) if x0 is None else x0 @@ -3122,14 +3152,15 @@ def step( if not self.preallocate: for ireg in range(self.nregsL1): self.d[ireg] = _softthreshold( - self.RegsL1[ireg].matvec(x) + self.b[ireg], self.epsRL1s[ireg] + self.RegsL1[ireg].matvec(x) + self.b[ireg], + self.threshRL1s[ireg], ) else: for ireg in range(self.nregsL1): self.ncp.add( self.RegsL1[ireg].matvec(x), self.b[ireg], out=self.d[ireg] ) - self.d[ireg] = _softthreshold(self.d[ireg], self.epsRL1s[ireg]) + self.d[ireg] = _softthreshold(self.d[ireg], self.threshRL1s[ireg]) # Bregman update for ireg in range(self.nregsL1): @@ -3262,6 +3293,7 @@ def solve( mu: float = 1.0, epsRL1s: SamplingLike | None = None, epsRL2s: SamplingLike | None = None, + lambdaRL1s: SamplingLike | None = None, tol: float = 1e-10, tau: float = 1.0, restart: bool = False, @@ -3298,14 +3330,21 @@ def solve( of ``RegsL2`` or equal to ``None`` to use a zero data for every regularization operator in ``RegsL2``) mu : :obj:`float`, optional - Data term damping + Data term damping epsRL1s : :obj:`list` - :math:`L_1` Regularization dampings (must have the same number of elements - as ``RegsL1``). Note that the effective weight of each :math:`L_1` - regularization term in the cost function is ``epsRL1s[i]**2`` + :math:`L_1` Regularization dampings (must have the same number of elements + as ``RegsL1``). Note that the effective weight of each :math:`L_1` + regularization term in the cost function is ``epsRL1s[i]**2`` epsRL2s : :obj:`list` - :math:`L_2` Regularization dampings (must have the same number of elements - as ``RegsL2``) + :math:`L_2` Regularization dampings (must have the same number of elements + as ``RegsL2``) + lambdaRL1s : :obj:`list`, optional + .. versionadded:: 2.9.0 + + Weights of the splitting terms (must have the same number of + elements as ``RegsL1``). They do not change the solution of the + problem, but affect the convergence speed of the algorithm. If + ``None``, they are set equal to ``epsRL1s`` tol : :obj:`float`, optional Tolerance. Stop outer iterations if difference between inverted model at subsequent iterations is smaller than ``tol`` @@ -3363,6 +3402,7 @@ def solve( mu=mu, epsRL1s=epsRL1s, epsRL2s=epsRL2s, + lambdaRL1s=lambdaRL1s, tol=tol, tau=tau, restart=restart, diff --git a/pylops/optimization/sparsity.py b/pylops/optimization/sparsity.py index 8dc9ea2c..13278b70 100644 --- a/pylops/optimization/sparsity.py +++ b/pylops/optimization/sparsity.py @@ -715,6 +715,7 @@ def splitbregman( mu: float = 1.0, epsRL1s: SamplingLike | None = None, epsRL2s: SamplingLike | None = None, + lambdaRL1s: SamplingLike | None = None, tol: float = 1e-10, rtol: float = 0.0, rtol1: float = 0.0, @@ -760,14 +761,21 @@ def splitbregman( of ``RegsL2`` or equal to ``None`` to use a zero data for every regularization operator in ``RegsL2``) mu : :obj:`float`, optional - Data term damping + Data term damping epsRL1s : :obj:`list` - :math:`L_1` Regularization dampings (must have the same number of elements - as ``RegsL1``). Note that the effective weight of each :math:`L_1` - regularization term in the cost function is ``epsRL1s[i]**2`` + :math:`L_1` Regularization dampings (must have the same number of elements + as ``RegsL1``). Note that the effective weight of each :math:`L_1` + regularization term in the cost function is ``epsRL1s[i]**2`` epsRL2s : :obj:`list` - :math:`L_2` Regularization dampings (must have the same number of elements - as ``RegsL2``) + :math:`L_2` Regularization dampings (must have the same number of elements + as ``RegsL2``) + lambdaRL1s : :obj:`list`, optional + .. versionadded:: 2.9.0 + + Weights of the splitting terms (must have the same number of + elements as ``RegsL1``). They do not change the solution of the + problem, but affect the convergence speed of the algorithm. If + ``None``, they are set equal to ``epsRL1s`` tol : :obj:`float`, optional Tolerance. Stop the solver if difference between inverted model at subsequent iterations is smaller than ``tol`` @@ -850,6 +858,7 @@ def splitbregman( mu=mu, epsRL1s=epsRL1s, epsRL2s=epsRL2s, + lambdaRL1s=lambdaRL1s, tol=tol, tau=tau, restart=restart, diff --git a/pytests/test_sparsity.py b/pytests/test_sparsity.py index 69b6e998..d35b1668 100644 --- a/pytests/test_sparsity.py +++ b/pytests/test_sparsity.py @@ -604,6 +604,20 @@ def test_SPGL1(par): assert_array_almost_equal(x, xinv, decimal=1) +@pytest.mark.parametrize("par", [(par1)]) +def test_SplitBregman_lambdaRL1s_wrong(par): + """Check errors for lambdaRL1s with wrong number of elements or values""" + nx = 3 * par["nx"] + Iop = Identity(nx) + Dop = FirstDerivative(nx, edge=True) + y = np.ones(nx) + + with pytest.raises(ValueError, match="same number of elements"): + splitbregman(Iop, y, [Dop], epsRL1s=[0.3], lambdaRL1s=[0.3, 0.3]) + with pytest.raises(ValueError, match="strictly positive"): + splitbregman(Iop, y, [Dop], epsRL1s=[0.3], lambdaRL1s=[0.0]) + + @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j)]) def test_SplitBregman(par): """Invert denoise problem with SplitBregman""" @@ -718,3 +732,30 @@ def test_SplitBregman_cost(par): + epsRL1**2 * np.linalg.norm(D1op @ xinv, ord=1) ) assert_array_almost_equal(cost[0], J, decimal=6) + + +@pytest.mark.parametrize("par", [(par1), (par1j)]) +def test_SplitBregman_lambdaRL1s(par): + """Check that lambdaRL1s=None is equivalent to lambdaRL1s=epsRL1s""" + np.random.seed(42) + nx = 3 * par["nx"] + Iop = Identity(nx) + Dop = FirstDerivative(nx, edge=True) + + x = np.zeros(nx) + x[: nx // 2] = 10 + x[nx // 2 : 3 * nx // 4] = -5 + y = x + np.random.normal(0, 1, nx) + epsRL1 = 1.0 + + kwars_solver = ( + dict(iter_lim=20, damp=0) if backend == "numpy" else dict(niter=20, damp=0) + ) + kwargs_sb = dict( + niter_outer=2000, niter_inner=5, mu=1.0, epsRL1s=[epsRL1], tol=1e-10 + ) + xdef = splitbregman(Iop, y, [Dop], **kwargs_sb, **kwars_solver)[0] + xeps = splitbregman( + Iop, y, [Dop], lambdaRL1s=[epsRL1], **kwargs_sb, **kwars_solver + )[0] + assert_array_almost_equal(xdef, xeps, decimal=12)