pylops.optimization.cls_sparsity.ISTA

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

Iterative Shrinkage-Thresholding Algorithm (ISTA).

Solve an optimization problem with \(L_p, \; p=0, 0.5, 1\) regularization, given the operator Op and data y. The operator can be real or complex, and should ideally be either square \(N=M\) or underdetermined \(N<M\).

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.

Opmatveccallable

Function handle to Op.matvec or Op.matmat depending on the number of dimensions of y.

Oprmatveccallable

Function handle to Op.rmatvec or Op.rmatmat depending on the number of dimensions of y.

SOpmatveccallable

Function handle to SOp.matvec or SOp.matmat depending on the number of dimensions of y.

SOprmatveccallable

Function handle to SOp.rmatvec or SOp.rmatmat depending on the number of dimensions of y.

threshfcallable

Function handle to the chosen thresholding method.

threshfloat

Threshold.

resnumpy.ndarray

Residual vector of size \([N \times 1]\) used in the solver when preallocate=True. Available only after setup is called and updated at each call to step.

gradnumpy.ndarray

Gradient vector of size \([M \times 1]\) used in the solver when preallocate=True. Available only after setup is called and updated at each call to step.

x_untheshnumpy.ndarray

Unthresholded model vector of size \([M \times 1]\) used in the solver when preallocate=True. Available only after setup is called and updated at each call to step.

xoldnumpy.ndarray

Old model vector of size \([M \times 1]\) used in the solver when preallocate=True. Available only after setup is called and updated at each call to step.

SOpx_untheshnumpy.ndarray

Old model vector pre-multiplied by the regularization operator of size \([M_S \times 1]\) used in the solver when preallocate=True. Available only after setup is called and updated at each call to step.

normresoldfloat

Old norm of the residual.

tfloat

FISTA auxiliary coefficient (not used in ISTA).

costlist

History of the L2 norm of the total objectiv function. Available only after setup is called and updated at each call to step.

iiterint

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

Raises:

NotImplementedError

If threshkind is different from hard, soft, half, soft-percentile, or half-percentile

ValueError

If perc=None when threshkind is soft-percentile or half-percentile

ValueError

If monitorres=True and residual increases

See also

OMP

Orthogonal Matching Pursuit (OMP).

FISTA

Fast Iterative Shrinkage-Thresholding Algorithm (FISTA).

SPGL1

Spectral Projected-Gradient for L1 norm (SPGL1).

SplitBregman

Split Bregman for mixed L2-L1 norms.

Notes

Solves the following synthesis problem for the operator \(\mathbf{Op}\) and the data \(\mathbf{y}\):

\[J = \|\mathbf{y} - \mathbf{Op}\,\mathbf{x}\|_2^2 + \epsilon \|\mathbf{x}\|_p\]

or the analysis problem:

\[J = \|\mathbf{y} - \mathbf{Op}\,\mathbf{x}\|_2^2 + \epsilon \|\mathbf{SOp}^H\,\mathbf{x}\|_p\]

if SOp is provided. Note that in the first case, SOp should be assimilated in the modelling operator (i.e., Op=GOp * SOp).

The Iterative Shrinkage-Thresholding Algorithms (ISTA) [1] is used, where \(p=0, 0.5, 1\). This is a very simple iterative algorithm which applies the following step:

\[\mathbf{x}^{(i+1)} = T_{(\epsilon \alpha /2, p)} \left(\mathbf{x}^{(i)} + \alpha\,\mathbf{Op}^H \left(\mathbf{y} - \mathbf{Op}\,\mathbf{x}^{(i)}\right)\right)\]

or

\[\mathbf{x}^{(i+1)} = \mathbf{SOp}\,\left\{T_{(\epsilon \alpha /2, p)} \mathbf{SOp}^H\,\left(\mathbf{x}^{(i)} + \alpha\,\mathbf{Op}^H \left(\mathbf{y} - \mathbf{Op} \,\mathbf{x}^{(i)}\right)\right)\right\}\]

where \(\epsilon \alpha /2\) is the threshold and \(T_{(\tau, p)}\) is the thresholding rule. The most common variant of ISTA uses the so-called soft-thresholding rule \(T(\tau, p=1)\). Alternatively an hard-thresholding rule is used in the case of \(p=0\) or a half-thresholding rule is used in the case of \(p=1/2\). Finally, percentile bases thresholds are also implemented: the damping factor is not used anymore an the threshold changes at every iteration based on the computed percentile.

[1]

Daubechies, I., Defrise, M., and De Mol, C., “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint”, Communications on pure and applied mathematics, vol. 57, pp. 1413-1457. 2004.

Methods

__init__(Op[, callbacks])
callback(x, *args, **kwargs)	Callback routine
finalize([show])	Finalize solver
memory_usage([show, unit])	Compute memory usage of the solver
run(x[, niter, show, itershow])	Run solver
setup(y[, x0, niter, SOp, eps, alpha, ...])	Setup solver
solve(y[, x0, niter, SOp, eps, alpha, ...])	Run entire solver
step(x[, show])	Run one step of solver

Examples using pylops.optimization.cls_sparsity.ISTA

03. Solvers

03. Solvers

03. Solvers (Advanced)

03. Solvers (Advanced)