pylops.optimization.cls_sparsity.SplitBregman

class pylops.optimization.cls_sparsity.SplitBregman(Op, callbacks=None)

Split Bregman for mixed L2-L1 norms.

Solve an unconstrained system of equations with mixed \(L_2\) and \(L_1\) regularization terms given the operator Op, a list of \(L_1\) regularization terms RegsL1, and an optional list of \(L_2\) regularization terms RegsL2.

Parameters:

Oppylops.LinearOperator

Operator to invert

Attributes:

ncpmodule

Array module used by the solver (obtained via pylops.utils.backend.get_array_module) ). Available only after setup is called.

isjaxbool

Whether the input data is a JAX array or not.

nregsL1int

Number of L1 regularization terms.

blist

Bregman variables (one per \(L_1\) regularization term).

dlist

Split variables obtained by shrinkage (one per \(L_1\) regularization term).

nregsL2int

Number of L2 regularization terms.

Regslist

List of L1 and L2 regularization terms.

epsRslist

List of L1 and L2 regularization dampings.

threshRL1slist

List of thresholds of the shrinkage steps (one per \(L_1\) regularization term).

costnumpy.ndarray, optional

History of total cost function through iterations.

iiterint

Current iteration number. Available only after setup is called and updated at each call to step.

Notes

Solve the following system of unconstrained, regularized equations given the operator \(\mathbf{Op}\) and a set of mixed norm (\(L_2\) and \(L_1\)) regularization terms \(\mathbf{R}_{2,i}\) and \(\mathbf{R}_{1,i}\), respectively:

\[J = \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 + \sum_i \epsilon_{\mathbf{R}_{1,i}}^2 \| \mathbf{R}_{1,i} \textbf{x} \|_1\]

where \(\mu\) is the reconstruction damping, \(\epsilon_{\mathbf{R}_{2,i}}\) are the damping factors used to weight the different \(L_2\) regularization terms of the cost function and \(\epsilon_{\mathbf{R}_{1,i}}\) are the damping factors of the different \(L_1\) regularization terms of the cost function. Note that the effective weight of each \(L_1\) regularization term is \(\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 inverse problems and Bregman updates.

The original system of equations is initially converted into a constrained problem:

\[J = \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 + \sum_i \epsilon_{\mathbf{R}_{1,i}}^2 \| \textbf{d}_i \|_1 \quad \text{subject to} \quad \textbf{d}_i = \mathbf{R}_{1,i} \textbf{x} \quad \forall i\]

and solved as follows:

\[\begin{split}\begin{aligned} (\textbf{x}^{k+1}, \textbf{d}_i^{k+1}) = \operatorname*{arg\,min}_{\mathbf{x}, \mathbf{d}_i} \; & \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 \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}\end{split}\]

\[\textbf{b}_i^{k+1}=\textbf{b}_i^k + \tau (\mathbf{R}_{1,i} \textbf{x}^{k+1} - \textbf{d}_i^{k+1})\]

where \(\textbf{d}_i\) are the split variables, \(\textbf{b}_i\) are the Bregman variables, \(\lambda_{\mathbf{R}_{1,i}}\) are the weights of the splitting terms, and \(\tau\) is a scaling factor of the Bregman update (\(\tau=1\) in the original algorithm [1]). The weights \(\lambda_{\mathbf{R}_{1,i}}\) do not change the solution of the problem, but affect the conditioning of the \(\textbf{x}\)-subproblem and the convergence speed of the algorithm. When not provided, \(\lambda_{\mathbf{R}_{1,i}} = \epsilon_{\mathbf{R}_{1,i}}\).

The first step is solved by alternating niter_inner times the minimization over \(\textbf{x}\) and over \(\textbf{d}_i\). The \(\textbf{x}\)-subproblem is an \(L_2\)-regularized least-squares problem, solved with scipy.sparse.linalg.lsqr (engine="scipy" and NumPy arrays) or with pylops.optimization.solver.cgls (engine="pylops", or any other array type, e.g., CuPy or JAX arrays). The \(\textbf{d}_i\)-subproblems are solved in closed form by soft thresholding \(\mathbf{R}_{1,i} \textbf{x}^{k+1} + \textbf{b}_i^k\) with threshold \(\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.

[1] (1,2)

Goldstein T. and Osher S., “The Split Bregman Method for L1-Regularized Problems”, SIAM J. on Scientific Computing, vol. 2(2), pp. 323-343. 2008.

Methods

__init__(Op[, callbacks])
callback(x, *args, **kwargs)	Callback routine
finalize([show])	Finalize solver
memory_usage([nopRegsL1, nopRegsL2, show, unit])	Compute memory usage of the solver
run(x[, engine, show, itershow, show_inner])	Run solver
setup(y, RegsL1[, x0, niter_outer, ...])	Setup solver
solve(y, RegsL1[, x0, niter_outer, ...])	Run entire solver
step(x[, engine, show, show_inner])	Run one step of solver