Skip to content

Feat: added lambdaRL1s to SplitBregman #810

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 3 commits into from
Oct 7, 2026
Merged

Feat: added lambdaRL1s to SplitBregman #810

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2499a2f
feat: added lambdaRL1s to SplitBregman
mrava87 Oct 6, 2026
0b33b54
minor: move lambdaRL1s next to epsRL1s/epsRL2s and fix docstring inde…
mrava87 Oct 6, 2026
2f7e00d
test: remove lambdaRL1s convergence comparison from SplitBregman test
mrava87 Oct 6, 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
96 changes: 68 additions & 28 deletions pylops/optimization/cls_sparsity.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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`
Expand All @@ -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
Expand All @@ -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}
Expand All @@ -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
Expand All @@ -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``.
Expand All @@ -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)
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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``
Expand Down Expand Up @@ -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
Expand All @@ -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
]
Expand All @@ -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
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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``
Expand Down Expand Up @@ -3363,6 +3402,7 @@ def solve(
mu=mu,
epsRL1s=epsRL1s,
epsRL2s=epsRL2s,
lambdaRL1s=lambdaRL1s,
tol=tol,
tau=tau,
restart=restart,
Expand Down
21 changes: 15 additions & 6 deletions pylops/optimization/sparsity.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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``
Expand Down Expand Up @@ -850,6 +858,7 @@ def splitbregman(
mu=mu,
epsRL1s=epsRL1s,
epsRL2s=epsRL2s,
lambdaRL1s=lambdaRL1s,
tol=tol,
tau=tau,
restart=restart,
Expand Down
41 changes: 41 additions & 0 deletions pytests/test_sparsity.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"""
Expand Down Expand Up @@ -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)
Loading