# Source code for pylops.optimization.sparsity

import logging
import time

import numpy as np
from scipy.sparse.linalg import lsqr

from pylops import LinearOperator
from pylops.basicoperators import Diagonal, Identity
from pylops.optimization.eigs import power_iteration
from pylops.optimization.leastsquares import (
NormalEquationsInversion,
RegularizedInversion,
)
from pylops.optimization.solver import cgls
from pylops.utils.backend import get_array_module, get_module_name, to_numpy

try:
from spgl1 import spgl1
except ModuleNotFoundError:
spgl1 = None
spgl1_message = "Spgl1 not installed. " 'Run "pip install spgl1".'
except Exception as e:
spgl1 = None
spgl1_message = "Failed to import spgl1 (error:%s)." % e

def _hardthreshold(x, thresh):
r"""Hard thresholding.

Applies hard thresholding to vector x (equal to the proximity
operator for :math:\|\mathbf{x}\|_0) as shown in _.

..  Chen, F., Shen, L., Suter, B.W., “Computing the proximity
operator of the ℓp norm with 0 < p < 1”,
IET Signal Processing, vol. 10. 2016.

Parameters
----------
x : :obj:numpy.ndarray
Vector
thresh : :obj:float
Threshold

Returns
-------
x1 : :obj:numpy.ndarray
Tresholded vector

"""
x1 = x.copy()
x1[np.abs(x) <= np.sqrt(2 * thresh)] = 0
return x1

def _softthreshold(x, thresh):
r"""Soft thresholding.

Applies soft thresholding to vector x (equal to the proximity
operator for :math:\|\mathbf{x}\|_1) as shown in _.

..  Chen, F., Shen, L., Suter, B.W., “Computing the proximity
operator of the ℓp norm with 0 < p < 1”,
IET Signal Processing, vol. 10. 2016.

Parameters
----------
x : :obj:numpy.ndarray
Vector
thresh : :obj:float
Threshold

Returns
-------
x1 : :obj:numpy.ndarray
Tresholded vector

"""
if np.iscomplexobj(x):
# https://stats.stackexchange.com/questions/357339/soft-thresholding-
# for-the-lasso-with-complex-valued-data
x1 = np.maximum(np.abs(x) - thresh, 0.0) * np.exp(1j * np.angle(x))
else:
x1 = np.maximum(np.abs(x) - thresh, 0.0) * np.sign(x)
return x1

def _halfthreshold(x, thresh):
r"""Half thresholding.

Applies half thresholding to vector x (equal to the proximity
operator for :math:\|\mathbf{x}\|_{1/2}^{1/2}) as shown in _.

..  Chen, F., Shen, L., Suter, B.W., “Computing the proximity
operator of the ℓp norm with 0 < p < 1”,
IET Signal Processing, vol. 10. 2016.

Parameters
----------
x : :obj:numpy.ndarray
Vector
thresh : :obj:float
Threshold

Returns
-------
x1 : :obj:numpy.ndarray
Tresholded vector

.. warning::
Since version 1.17.0 does not produce np.nan on bad input.

"""
arg = np.ones_like(x)
arg[x != 0] = (thresh / 8.0) * (np.abs(x[x != 0]) / 3.0) ** (-1.5)
arg = np.clip(arg, -1, 1)
phi = 2.0 / 3.0 * np.arccos(arg)
x1 = 2.0 / 3.0 * x * (1 + np.cos(2.0 * np.pi / 3.0 - phi))
# x1[np.abs(x) <= 1.5 * thresh ** (2. / 3.)] = 0
x1[np.abs(x) <= (54 ** (1.0 / 3.0) / 4.0) * thresh ** (2.0 / 3.0)] = 0
return x1

def _hardthreshold_percentile(x, perc):
r"""Percentile Hard thresholding.

Applies hard thresholding to vector x using a percentile to define
the amount of values in the input vector to be preserved as shown in _.

..  Chen, Y., Chen, K., Shi, P., Wang, Y., “Irregular seismic
data reconstruction using a percentile-half-thresholding algorithm”,
Journal of Geophysics and Engineering, vol. 11. 2014.

Parameters
----------
x : :obj:numpy.ndarray
Vector
thresh : :obj:float
Threshold

Returns
-------
x1 : :obj:numpy.ndarray
Tresholded vector

"""
thresh = np.percentile(np.abs(x), perc)
return _hardthreshold(x, 0.5 * thresh ** 2)

def _softthreshold_percentile(x, perc):
r"""Percentile Soft thresholding.

Applies soft thresholding to vector x using a percentile to define
the amount of values in the input vector to be preserved as shown in _.

..  Chen, Y., Chen, K., Shi, P., Wang, Y., “Irregular seismic
data reconstruction using a percentile-half-thresholding algorithm”,
Journal of Geophysics and Engineering, vol. 11. 2014.

Parameters
----------
x : :obj:numpy.ndarray
Vector
perc : :obj:float
Percentile

Returns
-------
x : :obj:numpy.ndarray
Tresholded vector

"""
thresh = np.percentile(np.abs(x), perc)
return _softthreshold(x, thresh)

def _halfthreshold_percentile(x, perc):
r"""Percentile Half thresholding.

Applies half thresholding to vector x using a percentile to define
the amount of values in the input vector to be preserved as shown in _.

..  Xu, Z., Xiangyu, C., Xu, F. and Zhang, H., “L1/2 Regularization: A
Thresholding Representation Theory and a Fast Solver”, IEEE Transactions
on Neural Networks and Learning Systems, vol. 23. 2012.

Parameters
----------
x : :obj:numpy.ndarray
Vector
perc : :obj:float
Percentile

Returns
-------
x : :obj:numpy.ndarray
Tresholded vector

"""
thresh = np.percentile(np.abs(x), perc)
# return _halfthreshold(x, (2. / 3. * thresh) ** (1.5))
return _halfthreshold(x, (4.0 / 54 ** (1.0 / 3.0) * thresh) ** 1.5)

def _IRLS_data(
Op,
data,
nouter,
threshR=False,
epsR=1e-10,
epsI=1e-10,
x0=None,
tolIRLS=1e-10,
returnhistory=False,
**kwargs_solver
):
r"""Iteratively reweighted least squares with L1 data term"""
ncp = get_array_module(data)

if x0 is not None:
data = data - Op * x0
if returnhistory:
xinv_hist = ncp.zeros((nouter + 1, int(Op.shape)))
rw_hist = ncp.zeros((nouter + 1, int(Op.shape)))

# first iteration (unweighted least-squares)
xinv = NormalEquationsInversion(
Op, None, data, epsI=epsI, returninfo=False, **kwargs_solver
)
r = data - Op * xinv
if returnhistory:
xinv_hist = xinv
for iiter in range(nouter):
# other iterations (weighted least-squares)
xinvold = xinv.copy()
if threshR:
rw = 1.0 / ncp.maximum(ncp.abs(r), epsR)
else:
rw = 1.0 / (ncp.abs(r) + epsR)
rw = rw / rw.max()
R = Diagonal(rw)
xinv = NormalEquationsInversion(
Op, [], data, Weight=R, epsI=epsI, returninfo=False, **kwargs_solver
)
r = data - Op * xinv
# save history
if returnhistory:
rw_hist[iiter] = rw
xinv_hist[iiter + 1] = xinv
# check tolerance
if ncp.linalg.norm(xinv - xinvold) < tolIRLS:
nouter = iiter
break

if x0 is not None:
xinv = x0 + xinv
if returnhistory:
xinv_hist = x0 + xinv_hist

if returnhistory:
return xinv, nouter, xinv_hist[: nouter + 1], rw_hist[: nouter + 1]
else:
return xinv, nouter

def _IRLS_model(
Op,
data,
nouter,
threshR=False,
epsR=1e-10,
epsI=1e-10,
x0=None,
tolIRLS=1e-10,
returnhistory=False,
**kwargs_solver
):
r"""Iteratively reweighted least squares with L1 model term"""
ncp = get_array_module(data)

if x0 is not None:
data = data - Op * x0
if returnhistory:
xinv_hist = ncp.zeros((nouter + 1, int(Op.shape)))
rw_hist = ncp.zeros((nouter + 1, int(Op.shape)))

Iop = Identity(data.size, dtype=data.dtype)
# first iteration (unweighted least-squares)
if ncp == np:
xinv = Op.H @ lsqr(Op @ Op.H + (epsI ** 2) * Iop, data, **kwargs_solver)
else:
xinv = (
Op.H
[docs]            @ cgls(
Op @ Op.H + (epsI ** 2) * Iop,
data,
ncp.zeros(int(Op.shape), dtype=Op.dtype),
**kwargs_solver
)
)
if returnhistory:
xinv_hist = xinv
for iiter in range(nouter):
# other iterations (weighted least-squares)
xinvold = xinv.copy()
rw = np.abs(xinv)
rw = rw / rw.max()
R = Diagonal(rw, dtype=rw.dtype)
if ncp == np:
xinv = (
R
@ Op.H
@ lsqr(Op @ R @ Op.H + epsI ** 2 * Iop, data, **kwargs_solver)
)
else:
xinv = (
R
@ Op.H
@ cgls(
Op @ R @ Op.H + epsI ** 2 * Iop,
data,
ncp.zeros(int(Op.shape), dtype=Op.dtype),
**kwargs_solver
)
)
# save history
if returnhistory:
rw_hist[iiter] = rw
xinv_hist[iiter + 1] = xinv
# check tolerance
if np.linalg.norm(xinv - xinvold) < tolIRLS:
nouter = iiter
break

if x0 is not None:
xinv = x0 + xinv
if returnhistory:
xinv_hist = x0 + xinv_hist

if returnhistory:
return xinv, nouter, xinv_hist[: nouter + 1], rw_hist[: nouter + 1]
else:
return xinv, nouter

def IRLS(
Op,
data,
nouter,
threshR=False,
epsR=1e-10,
epsI=1e-10,
x0=None,
tolIRLS=1e-10,
returnhistory=False,
kind="data",
**kwargs_solver
):
r"""Iteratively reweighted least squares.

Solve an optimization problem with :math:L^1 cost function (data IRLS)
or :math:L^1 regularization term (model IRLS) given the operator Op
and data y.

In the *data IRLS*, the cost function is minimized by iteratively solving a
weighted least squares problem with the weight at iteration :math:i
being based on the data residual at iteration :math:i-1. This IRLS solver
is robust to *outliers* since the L1 norm given less weight to large
residuals than L2 norm does.

Similarly in the *model IRLS*, the weight at at iteration :math:i
is based on the model at iteration :math:i-1. This IRLS solver inverts
for a sparse model vector.

Parameters
----------
Op : :obj:pylops.LinearOperator
Operator to invert
data : :obj:numpy.ndarray
Data
nouter : :obj:int
Number of outer iterations
threshR : :obj:bool, optional
Apply thresholding in creation of weight (True)
or damping (False)
epsR : :obj:float, optional
Damping to be applied to residuals for weighting term
espI : :obj:float, optional
Tikhonov damping
x0 : :obj:numpy.ndarray, optional
Initial guess
tolIRLS : :obj:float, optional
Tolerance. Stop outer iterations if difference between inverted model
at subsequent iterations is smaller than tolIRLS
returnhistory : :obj:bool, optional
Return history of inverted model for each outer iteration of IRLS
kind : :obj:str, optional
Kind of solver (data or model)
**kwargs_solver
Arbitrary keyword arguments for
:py:func:scipy.sparse.linalg.cg solver for data IRLS and
:py:func:scipy.sparse.linalg.lsqr solver for model IRLS when using
numpy data(or :py:func:pylops.optimization.solver.cg and
:py:func:pylops.optimization.solver.cgls when using cupy data)

Returns
-------
xinv : :obj:numpy.ndarray
Inverted model
nouter : :obj:int
Number of effective outer iterations
xinv_hist : :obj:numpy.ndarray, optional
History of inverted model
rw_hist : :obj:numpy.ndarray, optional
History of weights

Notes
-----
*Data IRLS* solves the following optimization problem for the operator
:math:\mathbf{Op} and the data :math:\mathbf{d}:

.. math::
J = \|\mathbf{d} - \mathbf{Op}\,\mathbf{x}\|_1

by a set of outer iterations which require to repeatedly solve a
weighted least squares problem of the form:

.. math::
\DeclareMathOperator*{\argmin}{arg\,min}
\mathbf{x}^{(i+1)} = \argmin_\mathbf{x} \|\mathbf{d} -
\mathbf{Op}\,\mathbf{x}\|_{2, \mathbf{R}^{(i)}}^2 +
\epsilon_\mathbf{I}^2 \|\mathbf{x}\|_2^2

where :math:\mathbf{R}^{(i)} is a diagonal weight matrix
whose diagonal elements at iteration :math:i are equal to the absolute
inverses of the residual vector :math:\mathbf{r}^{(i)} =
\mathbf{y} - \mathbf{Op}\,\mathbf{x}^{(i)} at iteration :math:i.
More specifically the :math:j-th element of the diagonal of
:math:\mathbf{R}^{(i)} is

.. math::
R^{(i)}_{j,j} = \frac{1}{\left| r^{(i)}_j \right| + \epsilon_\mathbf{R}}

or

.. math::
R^{(i)}_{j,j} = \frac{1}{\max\{\left|r^{(i)}_j\right|, \epsilon_\mathbf{R}\}}

depending on the choice threshR. In either case,
:math:\epsilon_\mathbf{R} is the user-defined stabilization/thresholding
factor _.

Similarly *model IRLS* solves the following optimization problem for the
operator :math:\mathbf{Op} and the data :math:\mathbf{d}:

.. math::
\mathbf{d} = \mathbf{Op}\,\mathbf{x}

by a set of outer iterations which require to repeatedly solve a
weighted least squares problem of the form _:

.. math::
\mathbf{x}^{(i+1)} = \operatorname*{arg\,min}_\mathbf{x}
\mathbf{d} = \mathbf{Op}\,\mathbf{x}

where :math:\mathbf{R}^{(i)} is a diagonal weight matrix
whose diagonal elements at iteration :math:i are equal to the absolutes
of the model vector :math:\mathbf{x}^{(i)} at iteration
:math:i. More specifically the :math:j-th element of the diagonal of
:math:\mathbf{R}^{(i)} is

.. math::
R^{(i)}_{j,j} = \left|x^{(i)}_j\right|.

..  https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares
..  Chartrand, R., and Yin, W. "Iteratively reweighted algorithms for
compressive sensing", IEEE. 2008.

"""
if kind == "data":
solver = _IRLS_data
elif kind == "model":
solver = _IRLS_model
else:
raise NotImplementedError("kind must be data or model")
return solver(
Op,
data,
nouter,
threshR=threshR,
epsR=epsR,
epsI=epsI,
x0=x0,
tolIRLS=tolIRLS,
returnhistory=returnhistory,
**kwargs_solver
)

[docs]def OMP(
Op,
data,
niter_outer=10,
niter_inner=40,
sigma=1e-4,
normalizecols=False,
show=False,
):
r"""Orthogonal Matching Pursuit (OMP).

Solve an optimization problem with :math:L^0 regularization function given
the operator Op and data y. The operator can be real or complex,
and should ideally be either square :math:N=M or underdetermined
:math:N<M.

Parameters
----------
Op : :obj:pylops.LinearOperator
Operator to invert
data : :obj:numpy.ndarray
Data
niter_outer : :obj:int, optional
Number of iterations of outer loop
niter_inner : :obj:int, optional
Number of iterations of inner loop. By choosing niter_inner=0, the
Matching Pursuit (MP) algorithm is implemented.
sigma : :obj:list
Maximum :math:L^2 norm of residual. When smaller stop iterations.
normalizecols : :obj:list, optional
Normalize columns (True) or not (False). Note that this can be
expensive as it requires applying the forward operator
:math:n_{cols} times to unit vectors (i.e., containing 1 at
position j and zero otherwise); use only when the columns of the
operator are expected to have highly varying norms.
show : :obj:bool, optional
Display iterations log

Returns
-------
xinv : :obj:numpy.ndarray
Inverted model
iiter : :obj:int
Number of effective outer iterations
cost : :obj:numpy.ndarray, optional
History of cost function

--------
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
:math:\mathbf{Op} and the data :math:\mathbf{d}:

.. math::
\|\mathbf{Op}\,\mathbf{x}-\mathbf{b}\|_2^2 \leq \sigma^2,

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

.. math::
\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{b}\right\|_2^2

Note that by choosing niter_inner=0 the basic Matching Pursuit (MP)
optimization at each iteration to find the best :math:\mathbf{x} for the
currently selected basis functions, the vector :math:\mathbf{x} is just
updated at the new basis function by taking directly the value from
the inner product :math:\mathbf{Op}_j^H\,\mathbf{r}_k.

In this case it is highly reccomended 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.

"""
ncp = get_array_module(data)

Op = LinearOperator(Op)
if show:
tstart = time.time()
algname = "OMP optimization\n" if niter_inner > 0 else "MP optimization\n"
print(
algname
+ "-----------------------------------------------------------------\n"
"The Operator Op has %d rows and %d cols\n"
"sigma = %.2e\tniter_outer = %d\tniter_inner = %d\n"
"normalization=%s"
% (Op.shape, Op.shape, sigma, niter_outer, niter_inner, normalizecols)
)
# find normalization factor for each column
if normalizecols:
ncols = Op.shape
norms = ncp.zeros(ncols)
for icol in range(ncols):
unit = ncp.zeros(ncols, dtype=Op.dtype)
unit[icol] = 1
norms[icol] = np.linalg.norm(Op.matvec(unit))
if show:
print("-----------------------------------------------------------------")

if niter_inner == 0:
x = []
cols = []
res = data.copy()
cost = ncp.zeros(niter_outer + 1)
cost = np.linalg.norm(data)
iiter = 0
while iiter < niter_outer and cost[iiter] > sigma:
# compute inner products
cres = Op.rmatvec(res)
cres_abs = np.abs(cres)
if normalizecols:
cres_abs = cres_abs / norms
# choose column with max cres
cres_max = np.max(cres_abs)
imax = np.argwhere(cres_abs == cres_max).ravel()
nimax = len(imax)
if nimax > 0:
imax = imax[np.random.permutation(nimax)]
else:
imax = imax
# update active set
if imax not in cols:
cols.append(int(imax))
else:
imax_in_cols = cols.index(imax)

# estimate model for current set of columns
if niter_inner == 0:
# MP update
Opcol = Op.apply_columns(
[
int(imax),
]
)
res -= Opcol.matvec(cres[imax] * ncp.ones(1))
x.append(cres[imax])
else:
x[imax_in_cols] += cres[imax]
else:
# OMP update
Opcol = Op.apply_columns(cols)
if ncp == np:
x = lsqr(Opcol, data, iter_lim=niter_inner)
else:
x = cgls(
Opcol,
data,
ncp.zeros(int(Opcol.shape), dtype=Opcol.dtype),
niter=niter_inner,
)
res = data - Opcol.matvec(x)
iiter += 1
cost[iiter] = np.linalg.norm(res)
if show:
if iiter < 10 or niter_outer - iiter < 10 or iiter % 10 == 0:
msg = "%6g        %12.5e" % (iiter + 1, cost[iiter])
print(msg)
xinv = ncp.zeros(int(Op.shape), dtype=Op.dtype)
xinv[cols] = ncp.array(x)
if show:
print(
"\nIterations = %d        Total time (s) = %.2f"
% (iiter, time.time() - tstart)
)
print("-----------------------------------------------------------------\n")
return xinv, iiter, cost

[docs]def ISTA(
Op,
data,
niter,
eps=0.1,
alpha=None,
eigsiter=None,
eigstol=0,
tol=1e-10,
monitorres=False,
returninfo=False,
show=False,
threshkind="soft",
perc=None,
callback=None,
decay=None,
SOp=None,
x0=None,
):
r"""Iterative Shrinkage-Thresholding Algorithm (ISTA).

Solve an optimization problem with :math: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 :math:N=M
or underdetermined :math:N<M.

Parameters
----------
Op : :obj:pylops.LinearOperator
Operator to invert
data : :obj:numpy.ndarray
Data
niter : :obj:int
Number of iterations
eps : :obj:float, optional
Sparsity damping
alpha : :obj:float, optional
Step size. To guarantee convergence, ensure
:math:\alpha \le 1/\lambda_\text{max}, where :math:\lambda_\text{max}
is the largest eigenvalue of :math:\mathbf{Op}^H\mathbf{Op}.
If None, the maximum eigenvalue is estimated and the optimal step size
is chosen as :math:1/\lambda_\text{max}. If provided, the
convergence criterion will not be checked internally.
eigsiter : :obj:float, optional
Number of iterations for eigenvalue estimation if alpha=None
eigstol : :obj:float, optional
Tolerance for eigenvalue estimation if alpha=None
tol : :obj:float, optional
Tolerance. Stop iterations if difference between inverted model
at subsequent iterations is smaller than tol
monitorres : :obj:bool, optional
Monitor that residual is decreasing
returninfo : :obj:bool, optional
Return info of CG solver
show : :obj:bool, optional
Display iterations log
threshkind : :obj:str, optional
Kind of thresholding ('hard', 'soft', 'half', 'hard-percentile',
'soft-percentile', or 'half-percentile' - 'soft' used as default)
perc : :obj:float, optional
Percentile, as percentage of values to be kept by thresholding (to be
provided when thresholding is soft-percentile or half-percentile)
callback : :obj:callable, optional
Function with signature (callback(x)) to call after each iteration
where x is the current model vector
decay : :obj:numpy.ndarray, optional
Decay factor to be applied to thresholding during iterations
SOp : :obj:pylops.LinearOperator, optional
Regularization operator (use when solving the analysis problem)
x0: :obj:numpy.ndarray, optional
Initial guess

Returns
-------
xinv : :obj:numpy.ndarray
Inverted model
niter : :obj:int
Number of effective iterations
cost : :obj:numpy.ndarray, optional
History of cost function

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

--------
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
:math:\mathbf{Op} and the data :math:\mathbf{d}:

.. math::
J = \|\mathbf{d} - \mathbf{Op}\,\mathbf{x}\|_2^2 +
\epsilon \|\mathbf{x}\|_p

or the analysis problem:

.. math::
J = \|\mathbf{d} - \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) _ is used, where
:math:p=0, 0.5, 1. This is a very simple iterative algorithm which
applies the following step:

.. math::
\mathbf{x}^{(i+1)} = T_{(\epsilon \alpha /2, p)} \left(\mathbf{x}^{(i)} +
\alpha\,\mathbf{Op}^H \left(\mathbf{d} - \mathbf{Op}\,\mathbf{x}^{(i)}\right)\right)

or

.. math::
\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{d} -
\mathbf{Op} \,\mathbf{x}^{(i)}\right)\right)\right\}

where :math:\epsilon \alpha /2 is the threshold and :math:T_{(\tau, p)}
is the thresholding rule. The most common variant of ISTA uses the
so-called soft-thresholding rule :math:T(\tau, p=1). Alternatively an
hard-thresholding rule is used in the case of :math:p=0 or a half-thresholding
rule is used in the case of :math: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.

..  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.

"""
if not threshkind in [
"hard",
"soft",
"half",
"hard-percentile",
"soft-percentile",
"half-percentile",
]:
raise NotImplementedError(
"threshkind should be hard, soft, half,"
"hard-percentile, soft-percentile, "
"or half-percentile"
)
if (
threshkind in ["hard-percentile", "soft-percentile", "half-percentile"]
and perc is None
):
raise ValueError(
"Provide a percentile when choosing hard-percentile,"
"soft-percentile, or half-percentile thresholding"
)

# choose thresholding function
if threshkind == "soft":
threshf = _softthreshold
elif threshkind == "hard":
threshf = _hardthreshold
elif threshkind == "half":
threshf = _halfthreshold
elif threshkind == "hard-percentile":
threshf = _hardthreshold_percentile
elif threshkind == "soft-percentile":
threshf = _softthreshold_percentile
else:
threshf = _halfthreshold_percentile

# identify backend to use
ncp = get_array_module(data)

# prepare decay (if not passed)
if perc is None and decay is None:
decay = ncp.ones(niter)

if show:
tstart = time.time()
print(
"ISTA optimization (%s thresholding)\n"
"-----------------------------------------------------------\n"
"The Operator Op has %d rows and %d cols\n"
"eps = %10e\ttol = %10e\tniter = %d"
% (threshkind, Op.shape, Op.shape, eps, tol, niter)
)
# step size
if alpha is None:
if not isinstance(Op, LinearOperator):
Op = LinearOperator(Op, explicit=False)
# compute largest eigenvalues of Op^H * Op
Op1 = LinearOperator(Op.H * Op, explicit=False)
if get_module_name(ncp) == "numpy":
maxeig = np.abs(
Op1.eigs(
neigs=1,
symmetric=True,
niter=eigsiter,
**dict(tol=eigstol, which="LM")
)
)
else:
maxeig = np.abs(
power_iteration(
Op1, niter=eigsiter, tol=eigstol, dtype=Op1.dtype, backend="cupy"
)
)
alpha = 1.0 / float(maxeig)

# define threshold
thresh = eps * alpha * 0.5

if show:
if perc is None:
print("alpha = %10e\tthresh = %10e" % (alpha, thresh))
else:
print("alpha = %10e\tperc = %.1f" % (alpha, perc))
print("-----------------------------------------------------------\n")
head1 = "   Itn       x        r2norm     r12norm     xupdate"

# initialize model and cost function
if x0 is None:
if data.ndim == 1:
xinv = ncp.zeros(int(Op.shape), dtype=Op.dtype)
else:
xinv = ncp.zeros((int(Op.shape), data.shape), dtype=Op.dtype)
else:
if data.ndim != x0.ndim:
# error for wrong dimensions
raise ValueError("Number of columns of x0 and data are not the same")
elif x0.shape != Op.shape:
# error for wrong dimensions
raise ValueError("Operator and input vector have different dimensions")
else:
xinv = x0.copy()

if monitorres:
normresold = np.inf
if returninfo:
cost = np.zeros(niter + 1)

# iterate
for iiter in range(niter):
xinvold = xinv.copy()

# compute residual
res = data - Op @ xinv
if monitorres:
normres = np.linalg.norm(res)
if normres > normresold:
raise ValueError(
"ISTA stopped at iteration %d due to "
"residual increasing, consider modifying "
"eps and/or alpha..." % iiter
)
else:
normresold = normres

grad = alpha * (Op.H @ res)

# update inverted model
if SOp is not None:
xinv_unthesh = SOp.H @ xinv_unthesh
if perc is None:
xinv = threshf(xinv_unthesh, decay[iiter] * thresh)
else:
xinv = threshf(xinv_unthesh, 100 - perc)
if SOp is not None:
xinv = SOp @ xinv

# model update
xupdate = np.linalg.norm(xinv - xinvold)

if returninfo or show:
costdata = 0.5 * np.linalg.norm(res) ** 2
costreg = eps * np.linalg.norm(xinv, ord=1)
if returninfo:
cost[iiter] = costdata + costreg

# run callback
if callback is not None:
callback(xinv)

if show:
if iiter < 10 or niter - iiter < 10 or iiter % 10 == 0:
msg = "%6g  %12.5e  %10.3e   %9.3e  %10.3e" % (
iiter + 1,
to_numpy(xinv[:2]),
costdata,
costdata + costreg,
xupdate,
)
print(msg)

# check tolerance
if xupdate < tol:
logging.warning("update smaller that tolerance for " "iteration %d" % iiter)
niter = iiter
break

# get values pre-threshold at locations where xinv is different from zero
# xinv = np.where(xinv != 0, xinv_unthesh, xinv)

if show:
print(
"\nIterations = %d        Total time (s) = %.2f"
% (niter, time.time() - tstart)
)
print("---------------------------------------------------------\n")
if returninfo:
return xinv, niter, cost[:niter]
else:
return xinv, niter

[docs]def FISTA(
Op,
data,
niter,
eps=0.1,
alpha=None,
eigsiter=None,
eigstol=0,
tol=1e-10,
returninfo=False,
show=False,
threshkind="soft",
perc=None,
callback=None,
decay=None,
SOp=None,
x0=None,
):
r"""Fast Iterative Shrinkage-Thresholding Algorithm (FISTA).

Solve an optimization problem with :math: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
:math:N=M or underdetermined :math:N<M.

Parameters
----------
Op : :obj:pylops.LinearOperator
Operator to invert
data : :obj:numpy.ndarray
Data
niter : :obj:int
Number of iterations
eps : :obj:float, optional
Sparsity damping
Step size. To guarantee convergence, ensure
:math:\alpha \le 1/\lambda_\text{max}, where :math:\lambda_\text{max}
is the largest eigenvalue of :math:\mathbf{Op}^H\mathbf{Op}.
If None, the maximum eigenvalue is estimated and the optimal step size
is chosen as :math:1/\lambda_\text{max}. If provided, the
convergence criterion will not be checked internally.
eigsiter : :obj:int, optional
Number of iterations for eigenvalue estimation if alpha=None
eigstol : :obj:float, optional
Tolerance for eigenvalue estimation if alpha=None
tol : :obj:float, optional
Tolerance. Stop iterations if difference between inverted model
at subsequent iterations is smaller than tol
returninfo : :obj:bool, optional
Return info of FISTA solver
show : :obj:bool, optional
Display iterations log
threshkind : :obj:str, optional
Kind of thresholding ('hard', 'soft', 'half', 'soft-percentile', or
'half-percentile' - 'soft' used as default)
perc : :obj:float, optional
Percentile, as percentage of values to be kept by thresholding (to be
provided when thresholding is soft-percentile or half-percentile)
callback : :obj:callable, optional
Function with signature (callback(x)) to call after each iteration
where x is the current model vector
decay : :obj:numpy.ndarray, optional
Decay factor to be applied to thresholding during iterations
SOp : :obj:pylops.LinearOperator, optional
Regularization operator (use when solving the analysis problem)
x0: :obj:numpy.ndarray, optional
Initial guess

Returns
-------
xinv : :obj:numpy.ndarray
Inverted model
niter : :obj:int
Number of effective iterations
cost : :obj:numpy.ndarray, optional
History of cost function

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

--------
OMP: Orthogonal Matching Pursuit (OMP).
ISTA: Iterative Shrinkage-Thresholding Algorithm (ISTA).
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
:math:\mathbf{Op} and the data :math:\mathbf{d}:

.. math::
J = \|\mathbf{d} - \mathbf{Op}\,\mathbf{x}\|_2^2 +
\epsilon \|\mathbf{x}\|_p

or the analysis problem:

.. math::
J = \|\mathbf{d} - \mathbf{Op}\,\mathbf{x}\|_2^2 +
\epsilon \|\mathbf{SOp}^H\,\mathbf{x}\|_p

if SOp is provided.

The Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) _ is used,
where :math:p=0, 0.5, 1. This is a modified version of ISTA solver with
improved convergence properties and limited additional computational cost.
Similarly to the ISTA solver, the choice of the thresholding algorithm to
apply at every iteration is based on the choice of :math:p.

..  Beck, A., and Teboulle, M., “A Fast Iterative Shrinkage-Thresholding
Algorithm for Linear Inverse Problems”, SIAM Journal on
Imaging Sciences, vol. 2, pp. 183-202. 2009.

"""
if not threshkind in [
"hard",
"soft",
"half",
"hard-percentile",
"soft-percentile",
"half-percentile",
]:
raise NotImplementedError(
"threshkind should be hard, soft, half,"
"hard-percentile, soft-percentile, "
"or half-percentile"
)
if (
threshkind in ["hard-percentile", "soft-percentile", "half-percentile"]
and perc is None
):
raise ValueError(
"Provide a percentile when choosing hard-percentile,"
"soft-percentile, or half-percentile thresholding"
)

# choose thresholding function
if threshkind == "soft":
threshf = _softthreshold
elif threshkind == "hard":
threshf = _hardthreshold
elif threshkind == "half":
threshf = _halfthreshold
elif threshkind == "hard-percentile":
threshf = _hardthreshold_percentile
elif threshkind == "soft-percentile":
threshf = _softthreshold_percentile
else:
threshf = _halfthreshold_percentile

# identify backend to use
ncp = get_array_module(data)

# prepare decay (if not passed)
if perc is None and decay is None:
decay = ncp.ones(niter)

if show:
tstart = time.time()
print(
"FISTA optimization (%s thresholding)\n"
"-----------------------------------------------------------\n"
"The Operator Op has %d rows and %d cols\n"
"eps = %10e\ttol = %10e\tniter = %d"
% (threshkind, Op.shape, Op.shape, eps, tol, niter)
)
# step size
if alpha is None:
if not isinstance(Op, LinearOperator):
Op = LinearOperator(Op, explicit=False)
# compute largest eigenvalues of Op^H * Op
Op1 = LinearOperator(Op.H * Op, explicit=False)
if get_module_name(ncp) == "numpy":
maxeig = np.abs(
Op1.eigs(
neigs=1,
symmetric=True,
niter=eigsiter,
**dict(tol=eigstol, which="LM")
)
)
else:
maxeig = np.abs(
power_iteration(
Op1, niter=eigsiter, tol=eigstol, dtype=Op1.dtype, backend="cupy"
)
)
alpha = 1.0 / float(maxeig)

# define threshold
thresh = eps * alpha * 0.5

if show:
if perc is None:
print("alpha = %10e\tthresh = %10e" % (alpha, thresh))
else:
print("alpha = %10e\tperc = %.1f" % (alpha, perc))
print("-----------------------------------------------------------\n")
head1 = "   Itn       x        r2norm     r12norm     xupdate"

# initialize model and cost function
if x0 is None:
if data.ndim == 1:
xinv = ncp.zeros(int(Op.shape), dtype=Op.dtype)
else:
xinv = ncp.zeros((int(Op.shape), data.shape), dtype=Op.dtype)
else:
if data.ndim != x0.ndim:
# error for wrong dimensions
raise ValueError("Number of columns of x0 and data are not the same")
elif x0.shape != Op.shape:
# error for wrong dimensions
raise ValueError("Operator and input vector have different dimensions")
else:
xinv = x0.copy()

zinv = xinv.copy()
t = 1
if returninfo:
cost = np.zeros(niter + 1)

# iterate
for iiter in range(niter):
xinvold = xinv.copy()

# compute residual
resz = data - Op @ zinv

grad = alpha * (Op.H @ resz)

# update inverted model
if SOp is not None:
xinv_unthesh = SOp.H @ xinv_unthesh
if perc is None:
xinv = threshf(xinv_unthesh, decay[iiter] * thresh)
else:
xinv = threshf(xinv_unthesh, 100 - perc)
if SOp is not None:
xinv = SOp @ xinv

# update auxiliary coefficients
told = t
t = (1.0 + np.sqrt(1.0 + 4.0 * t ** 2)) / 2.0
zinv = xinv + ((told - 1.0) / t) * (xinv - xinvold)

# model update
xupdate = np.linalg.norm(xinv - xinvold)

if returninfo or show:
costdata = 0.5 * np.linalg.norm(data - Op @ xinv) ** 2
costreg = eps * np.linalg.norm(xinv, ord=1)
if returninfo:
cost[iiter] = costdata + costreg

# run callback
if callback is not None:
callback(xinv)

if show:
if iiter < 10 or niter - iiter < 10 or iiter % 10 == 0:
msg = "%6g  %12.5e  %10.3e   %9.3e  %10.3e" % (
iiter + 1,
to_numpy(xinv[:2]),
costdata,
costdata + costreg,
xupdate,
)
print(msg)

# check tolerance
if xupdate < tol:
niter = iiter
break

# get values pre-threshold  at locations where xinv is different from zero
# xinv = np.where(xinv != 0, xinv_unthesh, xinv)

if show:
print(
"\nIterations = %d        Total time (s) = %.2f"
% (niter, time.time() - tstart)
)
print("---------------------------------------------------------\n")
if returninfo:
return xinv, niter, cost[:niter]
else:
return xinv, niter

[docs]def SPGL1(Op, data, SOp=None, tau=0, sigma=0, x0=None, **kwargs_spgl1):

Solve a constrained system of equations given the operator Op
and a sparsyfing transform SOp aiming to retrive a model that
is sparse in the sparsyfing domain.

This is a simple wrapper to :py:func:spgl1.spgl1
which is a porting of the well-known
SPGL1 <https://www.cs.ubc.ca/~mpf/spgl1/>_ MATLAB solver into Python.
In order to be able to use this solver you need to have installed the
spgl1 library.

Parameters
----------
Op : :obj:pylops.LinearOperator
Operator to invert
data : :obj:numpy.ndarray
Data
SOp : :obj:pylops.LinearOperator
Sparsyfing transform
tau : :obj:float
Non-negative LASSO scalar. If different from 0,
SPGL1 will solve LASSO problem
sigma : :obj:list
BPDN scalar. If different from 0,
SPGL1 will solve BPDN problem
x0 : :obj:numpy.ndarray
Initial guess
**kwargs_spgl1
Arbitrary keyword arguments for
:py:func:spgl1.spgl1 solver

Returns
-------
xinv : :obj:numpy.ndarray
Inverted model in original domain.
pinv : :obj:numpy.ndarray
Inverted model in sparse domain.
info : :obj:dict
Dictionary with the following information:

- tau, final value of tau (see sigma above)

- rnorm, two-norm of the optimal residual

- rgap, relative duality gap (an optimality measure)

- gnorm, Lagrange multiplier of (LASSO)

- stat, final status of solver
* 1: found a BPDN solution,
* 2: found a BP solution; exit based on small gradient,
* 3: found a BP solution; exit based on small residual,
* 4: found a LASSO solution,
* 5: error, too many iterations,
* 6: error, linesearch failed,
* 7: error, found suboptimal BP solution,
* 8: error, too many matrix-vector products.

- niters, number of iterations

- nProdA, number of multiplications with A

- nProdAt, number of multiplications with A'

- n_newton, number of Newton steps

- time_project, projection time (seconds)

- time_matprod, matrix-vector multiplications time (seconds)

- time_total, total solution time (seconds)

- niters_lsqr, number of lsqr iterations (if subspace_min=True)

- xnorm1, L1-norm model solution history through iterations

- rnorm2, L2-norm residual history through iterations

- lambdaa, Lagrange multiplier history through iterations

Raises
------
ModuleNotFoundError
If the spgl1 library is not installed

Notes
-----
Solve different variations of sparsity-promoting inverse problem by
imposing sparsity in the retrieved model _.

The first problem is called *basis pursuit denoise (BPDN)* and
its cost function is

.. math::
\left\|\mathbf{Op}\,\mathbf{S}^H\mathbf{x}-\mathbf{b}\right\|_2^2
\leq \sigma,

while the second problem is the *ℓ₁-regularized least-squares or LASSO*
problem and its cost function is

.. math::
\left\|\mathbf{Op}\,\mathbf{S}^H\mathbf{x}-\mathbf{b}\right\|_2^2

..  van den Berg E., Friedlander M.P., "Probing the Pareto frontier
for basis pursuit solutions", SIAM J. on Scientific Computing,
vol. 31(2), pp. 890-912. 2008.

"""
if spgl1 is None:
raise ModuleNotFoundError(spgl1_message)
pinv, _, _, info = spgl1(
Op if SOp is None else Op * SOp.H,
data,
tau=tau,
sigma=sigma,
x0=x0,
**kwargs_spgl1
)

xinv = pinv.copy() if SOp is None else SOp.H * pinv
return xinv, pinv, info

[docs]def SplitBregman(
Op,
RegsL1,
data,
niter_outer=3,
niter_inner=5,
RegsL2=None,
dataregsL2=None,
mu=1.0,
epsRL1s=None,
epsRL2s=None,
tol=1e-10,
tau=1.0,
x0=None,
restart=False,
show=False,
**kwargs_lsqr
):
r"""Split Bregman for mixed L2-L1 norms.

Solve an unconstrained system of equations with mixed :math:L^2 and :math:L^1
regularization terms given the operator Op, a list of :math:L^1
regularization terms RegsL1, and an optional list of :math:L^2
regularization terms RegsL2.

Parameters
----------
Op : :obj:pylops.LinearOperator
Operator to invert
RegsL1 : :obj:list
:math:L^1 regularization operators
data : :obj:numpy.ndarray
Data
niter_outer : :obj:int
Number of iterations of outer loop
niter_inner : :obj:int
Number of iterations of inner loop of first step of the Split Bregman
algorithm. A small number of iterations is generally sufficient and
for many applications optimal efficiency is obtained when only one
iteration is performed.
RegsL2 : :obj:list
Additional :math:L^2 regularization operators
(if None, :math:L^2 regularization is not added to the problem)
dataregsL2 : :obj:list, optional
:math:L^2 Regularization data (must have the same number of elements
of RegsL2 or equal to None to use a zero data for every
regularization operator in RegsL2)
mu : :obj:float, optional
Data term damping
epsRL1s : :obj:list
:math:L^1 Regularization dampings (must have the same number of elements
as RegsL1)
epsRL2s : :obj:list
:math:L^2 Regularization dampings (must have the same number of elements
as RegsL2)
tol : :obj:float, optional
Tolerance. Stop outer iterations if difference between inverted model
at subsequent iterations is smaller than tol
tau : :obj:float, optional
Scaling factor in the Bregman update (must be close to 1)
x0 : :obj:numpy.ndarray, optional
Initial guess
restart : :obj:bool, optional
The unconstrained inverse problem in inner loop is initialized with
the initial guess (True) or with the last estimate (False)
show : :obj:bool, optional
Display iterations log
**kwargs_lsqr
Arbitrary keyword arguments for
:py:func:scipy.sparse.linalg.lsqr solver used to solve the first
subproblem in the first step of the Split Bregman algorithm.

Returns
-------
xinv : :obj:numpy.ndarray
Inverted model
itn_out : :obj:int
Iteration number of outer loop upon termination

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

.. math::
J = \frac{\mu}{2} \|\textbf{d} - \textbf{Op}\,\textbf{x} \|_2^2 +
\frac{1}{2}\sum_i \epsilon_{\mathbf{R}_{2,i}} \|\mathbf{d}_{\mathbf{R}_{2,i}} -
\mathbf{R}_{2,i} \textbf{x} \|_2^2 +
\sum_i \epsilon_{\mathbf{R}_{1,i}} \| \mathbf{R}_{1,i} \textbf{x} \|_1

where :math:\mu is the reconstruction damping, :math:\epsilon_{\mathbf{R}_{2,i}}
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 used to weight the different :math:L^1 regularization
terms of the cost function.

The generalized Split-Bergman algorithm _ is used to solve such cost
function: the algorithm is composed of a sequence of unconstrained

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

.. math::
J = \frac{\mu}{2} \|\textbf{d} - \textbf{Op}\,\textbf{x}\|_2^2 +
\frac{1}{2}\sum_i \epsilon_{\mathbf{R}_{2,i}} \|\mathbf{d}_{\mathbf{R}_{2,i}} -
\mathbf{R}_{2,i} \textbf{x}\|_2^2 +
\textbf{y}_i = \mathbf{R}_{1,i} \textbf{x} \quad \forall i

and solved as follows:

.. math::
\DeclareMathOperator*{\argmin}{arg\,min}
\begin{align}
(\textbf{x}^{k+1}, \textbf{y}_i^{k+1}) =
\argmin_{\mathbf{x}, \mathbf{y}_i}
\|\textbf{d} - \textbf{Op}\,\textbf{x}\|_2^2
&+ \frac{1}{2}\sum_i \epsilon_{\mathbf{R}_{2,i}} \|\mathbf{d}_{\mathbf{R}_{2,i}} -
\mathbf{R}_{2,i} \textbf{x}\|_2^2 \\
&+ \frac{1}{2}\sum_i \epsilon_{\mathbf{R}_{1,i}} \|\textbf{y}_i -
\mathbf{R}_{1,i} \textbf{x} - \textbf{b}_i^k\|_2^2 \\
&+ \sum_i \| \textbf{y}_i \|_1
\end{align}

.. math::
\textbf{b}_i^{k+1}=\textbf{b}_i^k +
(\mathbf{R}_{1,i} \textbf{x}^{k+1} - \textbf{y}^{k+1})

The :py:func:scipy.sparse.linalg.lsqr solver and a fast shrinkage
algorithm are used within a inner loop to solve the first step. The entire
procedure is repeated niter_outer times until convergence.

..  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.

"""
ncp = get_array_module(data)

if show:
tstart = time.time()
print(
"Split-Bregman optimization\n"
"---------------------------------------------------------\n"
"The Operator Op has %d rows and %d cols\n"
"niter_outer = %3d     niter_inner = %3d   tol = %2.2e\n"
"mu = %2.2e         epsL1 = %s\t  epsL2 = %s     "
% (
Op.shape,
Op.shape,
niter_outer,
niter_inner,
tol,
mu,
str(epsRL1s),
str(epsRL2s),
)
)
print("---------------------------------------------------------\n")
head1 = "   Itn          x           r2norm          r12norm"

# L1 regularizations
nregsL1 = len(RegsL1)
b = [ncp.zeros(RegL1.shape, dtype=Op.dtype) for RegL1 in RegsL1]
d = b.copy()

# L2 regularizations
nregsL2 = 0 if RegsL2 is None else len(RegsL2)
if nregsL2 > 0:
Regs = RegsL2 + RegsL1
if dataregsL2 is None:
dataregsL2 = [ncp.zeros(Reg.shape, dtype=Op.dtype) for Reg in RegsL2]
else:
Regs = RegsL1
dataregsL2 = []

# Rescale dampings
epsRs = [
np.sqrt(epsRL2s[ireg] / 2) / np.sqrt(mu / 2) for ireg in range(nregsL2)
] + [np.sqrt(epsRL1s[ireg] / 2) / np.sqrt(mu / 2) for ireg in range(nregsL1)]
xinv = ncp.zeros(Op.shape, dtype=Op.dtype) if x0 is None else x0
xold = ncp.full(Op.shape, ncp.inf, dtype=Op.dtype)

itn_out = 0
while ncp.linalg.norm(xinv - xold) > tol and itn_out < niter_outer:
xold = xinv
for _ in range(niter_inner):
# Regularized problem
dataregs = dataregsL2 + [d[ireg] - b[ireg] for ireg in range(nregsL1)]
xinv = RegularizedInversion(
Op,
Regs,
data,
dataregs=dataregs,
epsRs=epsRs,
x0=x0 if restart else xinv,
**kwargs_lsqr
)
# Shrinkage
d = [
_softthreshold(RegsL1[ireg] * xinv + b[ireg], epsRL1s[ireg])
for ireg in range(nregsL1)
]
# Bregman update
b = [b[ireg] + tau * (RegsL1[ireg] * xinv - d[ireg]) for ireg in range(nregsL1)]
itn_out += 1

if show:
costdata = mu / 2.0 * ncp.linalg.norm(data - Op.matvec(xinv)) ** 2
costregL2 = (
0
if RegsL2 is None
else [
epsRL2 * ncp.linalg.norm(dataregL2 - RegL2.matvec(xinv)) ** 2
for epsRL2, RegL2, dataregL2 in zip(epsRL2s, RegsL2, dataregsL2)
]
)
costregL1 = [
ncp.linalg.norm(RegL1.matvec(xinv), ord=1)
for epsRL1, RegL1 in zip(epsRL1s, RegsL1)
]
cost = (
costdata + ncp.sum(ncp.array(costregL2)) + ncp.sum(ncp.array(costregL1))
)
msg = "%6g  %12.5e       %10.3e        %9.3e" % (
ncp.abs(itn_out),
ncp.real(xinv),
costdata,
cost,
)
print(msg)

if show:
print(
"\nIterations = %d        Total time (s) = %.2f"
% (itn_out, time.time() - tstart)
)
print("---------------------------------------------------------\n")
return xinv, itn_out