pylops.optimization.cls_sparsity.OMP

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

Orthogonal Matching Pursuit (OMP).

Solve an optimization problem with \(L_0\) regularization function 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.

normsnumpy.ndarray

Vector of size \([Nbasis \times 1]\) containing the norms of each column of the Opbasis operator. Available only after setup is called.

resnumpy.ndarray

Residual vector of size \([N \times 1]\). Available only after setup is called and updated at each call to step.

costlist

History of the L2 norm of the residual. 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.

See also

ISTA

Iterative Shrinkage-Thresholding Algorithm (ISTA).

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 optimization problem for the operator \(\mathbf{Op}\) and the data \(\mathbf{y}\):

\[\|\mathbf{x}\|_0 \quad \text{subject to} \quad \|\mathbf{Op}\,\mathbf{x}-\mathbf{y}\|_2^2 \leq \sigma^2,\]

using Orthogonal Matching Pursuit (OMP). This is a very simple iterative algorithm which applies the following step:

\[\begin{split}\DeclareMathOperator*{\argmin}{arg\,min} \DeclareMathOperator*{\argmax}{arg\,max} \Lambda_k = \Lambda_{k-1} \cup \left\{\argmax_j \left|\mathbf{Op}^{j H}\,\mathbf{r}_k\right| \right\} \\ \mathbf{x}_k = \argmin_{\mathbf{x}} \left\|\mathbf{Op}_{\Lambda_k}\,\mathbf{x} - \mathbf{y}\right\|_2^2\end{split}\]

where \(\mathbf{Op}^j\) is the \(j\)-th column of the operator, \(\mathbf{r}_k\) is the residual at iteration \(k\), and \(\mathbf{Op}_{\Lambda_k}\) is the operator restricted to the columns in the set \(\Lambda_k\).

Note that by choosing niter_inner=0 the basic Matching Pursuit (MP) algorithm is implemented instead. In other words, instead of solving an optimization at each iteration to find the best \(\mathbf{x}\) for the currently selected basis functions, either the vector \(\mathbf{x}\) is just updated at the new basis function by adding the value from the inner product \(\mathbf{Op}_j^H\,\mathbf{r}_k\) to the current value (optimal_coeff=False) or the optimal coefficient that minimizes the norm of the residual \(\mathbf{r} - c * \mathbf{Op}^j\) is estimated (optimal_coeff=True) and added to the current value.

In the case the MP solver is used, it is highly recommended to provide a normalized basis function. If different basis have different norms, the solver is likely to diverge. Similar observations apply to OMP, even though mild unbalancing between the basis is generally properly handled. Two possible ways to handle the scenario fo non-normalized basis functions are:

• Find the normalization factor of the the basis functions before running the solver (this is done by choosing normalizecols=True);

• Find the optimal coefficient that minimizes the norm of the residual \(\mathbf{r} - c * \mathbf{Op}^j\) at every iteration (this is done by choosing optimal_coeff=True).

Finally, when the operator is a chain of operators, with the rigth-most representing the basis function, if the operator of the basis function is provided in the Opbasis parameter, the solver will use this operator to find the normalization factor for each column of the operator.

Methods

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