# Source code for pyproximal.optimization.primal

import time
import warnings
import numpy as np

from math import sqrt
from pylops.optimization.leastsquares import regularized_inversion
from pylops.utils.backend import to_numpy
from pyproximal.proximal import L2
from pyproximal.utils.bilinear import BilinearOperator

def _backtracking(x, tau, proxf, proxg, epsg, beta=0.5, niterback=10):
r"""Backtracking

Line-search algorithm for finding step sizes in proximal algorithms when
the Lipschitz constant of the operator is unknown (or expensive to
estimate).

"""
def ftilde(x, y, f, tau):
xy = x - y
return f(y) + np.dot(f.grad(y), xy) + \
(1. / (2. * tau)) * np.linalg.norm(xy) ** 2

iiterback = 0
while iiterback < niterback:
z = proxg.prox(x - tau * proxf.grad(x), epsg * tau)
ft = ftilde(z, x, proxf, tau)
if proxf(z) <= ft:
break
tau *= beta
iiterback += 1
return z, tau

[docs]def ProximalPoint(prox, x0, tau, niter=10, callback=None, show=False):
r"""Proximal point algorithm

Solves the following minimization problem using Proximal point algorithm:

.. math::

\mathbf{x} = \argmin_\mathbf{x} f(\mathbf{x})

where :math:f(\mathbf{x}) is any convex function that has a known
proximal operator.

Parameters
----------
prox : :obj:pyproximal.ProxOperator
Proximal operator
x0 : :obj:numpy.ndarray
Initial vector
tau : :obj:float
Positive scalar weight
niter : :obj:int, optional
Number of iterations of iterative scheme
callback : :obj:callable, optional
Function with signature (callback(x)) to call after each iteration
where x is the current model vector
show : :obj:bool, optional
Display iterations log

Returns
-------
x : :obj:numpy.ndarray
Inverted model

Notes
-----
The Proximal point algorithm can be expressed by the following recursion:

.. math::

\mathbf{x}^{k+1} = \prox_{\tau f}(\mathbf{x}^k)

"""
if show:
tstart = time.time()
print('Proximal point algorithm\n'
'---------------------------------------------------------\n'
'Proximal operator: %s\n'
'tau = %10e\tniter = %d\n' % (type(prox), tau, niter))
head = '   Itn       x[0]        f'
x = x0.copy()
for iiter in range(niter):
x = prox.prox(x, tau)

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

if show:
if iiter < 10 or niter - iiter < 10 or iiter % (niter // 10) == 0:
msg = '%6g  %12.5e  %10.3e' % \
(iiter + 1, x[0], prox(x))
print(msg)
if show:
print('\nTotal time (s) = %.2f' % (time.time() - tstart))
print('---------------------------------------------------------\n')
return x

tau=None, beta=0.5, eta=1.,
niter=10, niterback=100,
acceleration=None,
callback=None, show=False):

Solves the following minimization problem using (Accelerated) Proximal

.. math::

\mathbf{x} = \argmin_\mathbf{x} f(\mathbf{x}) + \epsilon g(\mathbf{x})

where :math:f(\mathbf{x}) is a smooth convex function with a uniquely
defined gradient and :math:g(\mathbf{x}) is any convex function that
has a known proximal operator.

Parameters
----------
proxf : :obj:pyproximal.ProxOperator
Proximal operator of f function (must have grad implemented)
proxg : :obj:pyproximal.ProxOperator
Proximal operator of g function
x0 : :obj:numpy.ndarray
Initial vector
epsg : :obj:float or :obj:np.ndarray, optional
Scaling factor of g function
tau : :obj:float or :obj:numpy.ndarray, optional
Positive scalar weight, which should satisfy the following condition
to guarantees convergence: :math:\tau  \in (0, 1/L] where L is
the Lipschitz constant of :math:\nabla f. When tau=None,
backtracking is used to adaptively estimate the best tau at each
iteration. Finally, note that :math:\tau can be chosen to be a vector
when dealing with problems with multiple right-hand-sides
beta : :obj:float, optional
Backtracking parameter (must be between 0 and 1)
eta : :obj:float, optional
Relaxation parameter (must be between 0 and 1, 0 excluded).
niter : :obj:int, optional
Number of iterations of iterative scheme
niterback : :obj:int, optional
Max number of iterations of backtracking
acceleration : :obj:str, optional
Acceleration (None, vandenberghe or fista)
callback : :obj:callable, optional
Function with signature (callback(x)) to call after each iteration
where x is the current model vector
show : :obj:bool, optional
Display iterations log

Returns
-------
x : :obj:numpy.ndarray
Inverted model

Notes
-----
The Proximal point algorithm can be expressed by the following recursion:

.. math::

\mathbf{x}^{k+1} = \mathbf{y}^k + \eta (\prox_{\tau^k \epsilon g}(\mathbf{y}^k -
\tau^k \nabla f(\mathbf{y}^k)) - \mathbf{y}^k) \\
\mathbf{y}^{k+1} = \mathbf{x}^k + \omega^k
(\mathbf{x}^k - \mathbf{x}^{k-1})

where at each iteration :math:\tau^k can be estimated by back-tracking
as follows:

.. math::

\begin{aligned}
&\tau = \tau^{k-1} &\\
&repeat \; \mathbf{z} = \prox_{\tau \epsilon g}(\mathbf{x}^k -
\tau \nabla f(\mathbf{x}^k)), \tau = \beta \tau \quad if \;
f(\mathbf{z}) \leq \tilde{f}_\tau(\mathbf{z}, \mathbf{x}^k) \\
&\tau^k = \tau, \quad \mathbf{x}^{k+1} = \mathbf{z} &\\
\end{aligned}

where :math:\tilde{f}_\tau(\mathbf{x}, \mathbf{y}) = f(\mathbf{y}) +
\nabla f(\mathbf{y})^T (\mathbf{x} - \mathbf{y}) +
1/(2\tau)||\mathbf{x} - \mathbf{y}||_2^2.

Different accelerations are provided:

- acceleration=None: :math:\omega^k = 0;
- acceleration=vandenberghe [1]_: :math:\omega^k = k / (k + 3) for
- acceleration=fista: :math:\omega^k = (t_{k-1}-1)/t_k for  where
:math:t_k = (1 + \sqrt{1+4t_{k-1}^{2}}) / 2 [2]_

.. [1] Vandenberghe, L., "Fast proximal gradient methods", 2010.
.. [2] 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.

"""
# check if epgs is a vector
if np.asarray(epsg).size == 1.:
epsg = epsg * np.ones(niter)
epsg_print = str(epsg[0])
else:
epsg_print = 'Multi'

if acceleration not in [None, 'None', 'vandenberghe', 'fista']:
raise NotImplementedError('Acceleration should be None, vandenberghe '
'or fista')
if show:
tstart = time.time()
'---------------------------------------------------------\n'
'Proximal operator (f): %s\n'
'Proximal operator (g): %s\n'
'tau = %s\tbeta=%10e\n'
'epsg = %s\tniter = %d\n'
''
'niterback = %d\tacceleration = %s\n' % (type(proxf), type(proxg),
'Adaptive' if tau is None else str(tau), beta,
epsg_print, niter, niterback, acceleration))
head = '   Itn       x[0]          f           g       J=f+eps*g       tau'

backtracking = False
if tau is None:
backtracking = True
tau = 1.

# initialize model
t = 1.
x = x0.copy()
y = x.copy()

# iterate
for iiter in range(niter):
xold = x.copy()

# proximal step
if not backtracking:
if eta == 1.:
x = proxg.prox(y - tau * proxf.grad(y), epsg[iiter] * tau)
else:
x = x + eta * (proxg.prox(x - tau * proxf.grad(x), epsg[iiter] * tau) - x)
else:
x, tau = _backtracking(y, tau, proxf, proxg, epsg[iiter],
beta=beta, niterback=niterback)
if eta != 1.:
x = x + eta * (proxg.prox(x - tau * proxf.grad(x), epsg[iiter] * tau) - x)

# update internal parameters for bilinear operator
if isinstance(proxf, BilinearOperator):
proxf.updatexy(x)

# update y
if acceleration == 'vandenberghe':
omega = iiter / (iiter + 3)
elif acceleration == 'fista':
told = t
t = (1. + np.sqrt(1. + 4. * t ** 2)) / 2.
omega = ((told - 1.) / t)
else:
omega = 0
y = x + omega * (x - xold)

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

if show:
if iiter < 10 or niter - iiter < 10 or iiter % (niter // 10) == 0:
pf, pg = proxf(x), proxg(x)
msg = '%6g  %12.5e  %10.3e  %10.3e  %10.3e  %10.3e' % \
(iiter + 1, np.real(to_numpy(x[0])) if x.ndim == 1 else np.real(to_numpy(x[0, 0])),
pf, pg,
pf + np.sum(epsg[iiter] * pg),
tau)
print(msg)
if show:
print('\nTotal time (s) = %.2f' % (time.time() - tstart))
print('---------------------------------------------------------\n')
return x

[docs]def AcceleratedProximalGradient(proxf, proxg, x0, tau=None, beta=0.5,
epsg=1., niter=10, niterback=100,
acceleration='vandenberghe',
callback=None, show=False):

This is a thin wrapper around :func:pyproximal.optimization.primal.ProximalGradient with
vandenberghe or fistaacceleration. See :func:pyproximal.optimization.primal.ProximalGradient
for details.

"""
'from v0.5.0. It is recommended to start using ProximalGradient by selecting the '
'appropriate acceleration parameter as this behaviour will become default in '
'version v1.0.0 and AcceleratedProximalGradient will be removed.', FutureWarning)
return ProximalGradient(proxf, proxg, x0, tau=tau, beta=beta,
epsg=epsg, niter=niter, niterback=niterback,
acceleration=acceleration,
callback=callback, show=show)

epsg=1., weights=None,
eta=1., niter=10,
acceleration=None,
callback=None, show=False):

Solves the following minimization problem using Generalized Proximal

.. math::

\mathbf{x} = \argmin_\mathbf{x} \sum_{i=1}^n f_i(\mathbf{x})
+ \sum_{j=1}^m \epsilon_j g_j(\mathbf{x}),~~n,m \in \mathbb{N}^+

where the :math:f_i(\mathbf{x}) are smooth convex functions with a uniquely
defined gradient and the :math:g_j(\mathbf{x}) are any convex function that
have a known proximal operator.

Parameters
----------
proxfs : :obj:list of pyproximal.ProxOperator
Proximal operators of the :math:f_i functions (must have grad implemented)
proxgs : :obj:list of pyproximal.ProxOperator
Proximal operators of the :math:g_j functions
x0 : :obj:numpy.ndarray
Initial vector
tau : :obj:float
Positive scalar weight, which should satisfy the following condition
to guarantees convergence: :math:\tau  \in (0, 1/L] where L is
the Lipschitz constant of :math:\sum_{i=1}^n \nabla f_i.
epsg : :obj:float or :obj:np.ndarray, optional
Scaling factor(s) of g function(s)
weights : :obj:float, optional
Weighting factors of g functions. Must sum to 1.
eta : :obj:float, optional
Relaxation parameter (must be between 0 and 1, 0 excluded). Note that
this will be only used when acceleration=None.
niter : :obj:int, optional
Number of iterations of iterative scheme
acceleration:  :obj:str, optional
Acceleration (None, vandenberghe or fista)
callback : :obj:callable, optional
Function with signature (callback(x)) to call after each iteration
where x is the current model vector
show : :obj:bool, optional
Display iterations log

Returns
-------
x : :obj:numpy.ndarray
Inverted model

Notes
-----
The Generalized Proximal point algorithm can be expressed by the
following recursion:

.. math::
\text{for } j=1,\cdots,n, \\
~~~~\mathbf z_j^{k+1} = \mathbf z_j^{k} + \eta
\left[prox_{\frac{\tau^k \epsilon_j}{w_j} g_j}\left(2 \mathbf{x}^{k} - \mathbf{z}_j^{k}
- \tau^k \sum_{i=1}^n \nabla f_i(\mathbf{x}^{k})\right) - \mathbf{x}^{k} \right] \\
\mathbf{x}^{k+1} = \sum_{j=1}^n w_j \mathbf z_j^{k+1} \\

where :math:\sum_{j=1}^n w_j=1. In the current implementation, :math:w_j=1/n when
not provided.

"""
# check if weights sum to 1
if weights is None:
weights = np.ones(len(proxgs)) / len(proxgs)
if len(weights) != len(proxgs) or np.sum(weights) != 1.:
raise ValueError(f'omega={weights} must be an array of size {len(proxgs)} '
f'summing to 1')
print(weights)

# check if epgs is a vector
if np.asarray(epsg).size == 1.:
epsg_print = str(epsg)
epsg = epsg * np.ones(len(proxgs))
else:
epsg_print = 'Multi'

if acceleration not in [None, 'None', 'vandenberghe', 'fista']:
raise NotImplementedError('Acceleration should be None, vandenberghe '
'or fista')
if show:
tstart = time.time()
'---------------------------------------------------------\n'
'Proximal operators (f): %s\n'
'Proximal operators (g): %s\n'
'tau = %10e\nepsg = %s\tniter = %d\n' % ([type(proxf) for proxf in proxfs],
[type(proxg) for proxg in proxgs],
0 if tau is None else tau,
epsg_print, niter))
head = '   Itn       x[0]          f           g       J=f+eps*g'

if tau is None:
tau = 1.

# initialize model
t = 1.
x = x0.copy()
y = x.copy()
zs = [x.copy() for _ in range(len(proxgs))]

# iterate
for iiter in range(niter):
xold = x.copy()

for i, proxf in enumerate(proxfs):

# proximal step
x = np.zeros_like(x)
for i, proxg in enumerate(proxgs):
ztmp = 2 * y - zs[i] - tau * grad
ztmp = proxg.prox(ztmp, tau * epsg[i] / weights[i])
zs[i] += eta * (ztmp - y)
x += weights[i] * zs[i]

# update y
if acceleration == 'vandenberghe':
omega = iiter / (iiter + 3)
elif acceleration == 'fista':
told = t
t = (1. + np.sqrt(1. + 4. * t ** 2)) / 2.
omega = ((told - 1.) / t)
else:
omega = 0
y = x + omega * (x - xold)

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

if show:
if iiter < 10 or niter - iiter < 10 or iiter % (niter // 10) == 0:
pf, pg = np.sum([proxf(x) for proxf in proxfs]), np.sum([proxg(x) for proxg in proxgs])
msg = '%6g  %12.5e  %10.3e  %10.3e  %10.3e' % \
(iiter + 1, x[0] if x.ndim == 1 else x[0, 0],
pf, pg[0] if epsg_print == 'Multi' else pg,
pf + np.sum(epsg * pg))
print(msg)
if show:
print('\nTotal time (s) = %.2f' % (time.time() - tstart))
print('---------------------------------------------------------\n')
return x

[docs]def HQS(proxf, proxg, x0, tau, niter=10, z0=None, gfirst=True,
callback=None, callbackz=False, show=False):

Solves the following minimization problem using Half Quadratic splitting
algorithm:

.. math::

\mathbf{x},\mathbf{z}  = \argmin_{\mathbf{x},\mathbf{z}}
f(\mathbf{x}) + g(\mathbf{z}) \\
s.t. \; \mathbf{x}=\mathbf{z}

where :math:f(\mathbf{x}) and :math:g(\mathbf{z}) are any convex
function that has a known proximal operator.

Parameters
----------
proxf : :obj:pyproximal.ProxOperator
Proximal operator of f function
proxg : :obj:pyproximal.ProxOperator
Proximal operator of g function
x0 : :obj:numpy.ndarray
Initial vector
tau : :obj:float or :obj:numpy.ndarray, optional
Positive scalar weight, which should satisfy the following condition
to guarantees convergence: :math:\tau  \in (0, 1/L] where L is
the Lipschitz constant of :math:\nabla f. Finally note that
:math:\tau can be chosen to be a vector of size niter such that
different :math:\tau is used at different iterations (i.e., continuation
strategy)
niter : :obj:int, optional
Number of iterations of iterative scheme
z0 : :obj:numpy.ndarray, optional
Initial z vector (not required when gfirst=True
gfirst : :obj:bool, optional
Apply Proximal of operator g first (True) or Proximal of
operator f first (False)
callback : :obj:callable, optional
Function with signature (callback(x)) to call after each iteration
where x is the current model vector
callbackz : :obj:bool, optional
Modify callback signature to (callback(x, z)) when callbackz=True
show : :obj:bool, optional
Display iterations log

Returns
-------
x : :obj:numpy.ndarray
Inverted model
z : :obj:numpy.ndarray
Inverted second model

Notes
-----
The HQS algorithm can be expressed by the following recursion [1]_:

.. math::

\mathbf{z}^{k+1} = \prox_{\tau g}(\mathbf{x}^{k}) \\
\mathbf{x}^{k+1} = \prox_{\tau f}(\mathbf{z}^{k+1})

for gfirst=False, or

.. math::

\mathbf{x}^{k+1} = \prox_{\tau f}(\mathbf{z}^{k}) \\
\mathbf{z}^{k+1} = \prox_{\tau g}(\mathbf{x}^{k+1})

for gfirst=False. Note that x and z converge to each other,
however if iterations are stopped too early x is guaranteed to belong to
the domain of f while z is guaranteed to belong to the domain of g.
Depending on the problem either of the two may be the best solution.

.. [1] D., Geman, and C., Yang, "Nonlinear image recovery with halfquadratic
regularization", IEEE Transactions on Image Processing,
4, 7, pp. 932-946, 1995.

"""
# check if epgs is a ve
if np.asarray(tau).size == 1.:
tau_print = str(tau)
tau = tau * np.ones(niter)
else:
tau_print = 'Variable'

if show:
tstart = time.time()
print('HQS\n'
'---------------------------------------------------------\n'
'Proximal operator (f): %s\n'
'Proximal operator (g): %s\n'
'tau = %s\tniter = %d\n' % (type(proxf), type(proxg),
tau_print, niter))
head = '   Itn       x[0]          f           g       J = f + g'

x = x0.copy()
if z0 is not None:
z = z0.copy()
else:
z = np.zeros_like(x)
for iiter in range(niter):
if gfirst:
z = proxg.prox(x, tau[iiter])
x = proxf.prox(z, tau[iiter])
else:
x = proxf.prox(z, tau[iiter])
z = proxg.prox(x, tau[iiter])

# run callback
if callback is not None:
if callbackz:
callback(x, z)
else:
callback(x)

if show:
if iiter < 10 or niter - iiter < 10 or iiter % (niter // 10) == 0:
pf, pg = proxf(x), proxg(x)
msg = '%6g  %12.5e  %10.3e  %10.3e  %10.3e' % \
(iiter + 1, np.real(to_numpy(x[0])),
pf, pg, pf + pg)
print(msg)
if show:
print('\nTotal time (s) = %.2f' % (time.time() - tstart))
print('---------------------------------------------------------\n')
return x, z

[docs]def ADMM(proxf, proxg, x0, tau, niter=10, gfirst=False,
callback=None, callbackz=False, show=False):
r"""Alternating Direction Method of Multipliers

Solves the following minimization problem using Alternating Direction
Method of Multipliers (also known as Douglas-Rachford splitting):

.. math::

\mathbf{x},\mathbf{z}  = \argmin_{\mathbf{x},\mathbf{z}}
f(\mathbf{x}) + g(\mathbf{z}) \\
s.t. \; \mathbf{x}=\mathbf{z}

where :math:f(\mathbf{x}) and :math:g(\mathbf{z}) are any convex
function that has a known proximal operator.

ADMM can also solve the problem of the form above with a more general
constraint: :math:\mathbf{Ax}+\mathbf{Bz}=c. This routine implements
the special case where :math:\mathbf{A}=\mathbf{I}, :math:\mathbf{B}=-\mathbf{I},
and :math:c=0, as a general algorithm can be obtained for any choice of
:math:f and :math:g provided they have a known proximal operator.

On the other hand, for more general choice of :math:\mathbf{A}, :math:\mathbf{B},
and :math:c, the iterations are not generalizable, i.e. thye depends on the choice of
:math:f and :math:g functions. For this reason, we currently only provide an additional
solver for the special case where :math:f is a :class:pyproximal.proximal.L2
operator with a linear operator  :math:\mathbf{G} and data :math:\mathbf{y},
:math:\mathbf{B}=-\mathbf{I} and :math:c=0,
called :func:pyproximal.optimization.primal.ADMML2. Note that for the very same choice
of :math:\mathbf{B} and :math:c, the :func:pyproximal.optimization.primal.LinearizedADMM
can also be used (and this does not require a specific choice of :math:f).

Parameters
----------
proxf : :obj:pyproximal.ProxOperator
Proximal operator of f function
proxg : :obj:pyproximal.ProxOperator
Proximal operator of g function
x0 : :obj:numpy.ndarray
Initial vector
tau : :obj:float, optional
Positive scalar weight, which should satisfy the following condition
to guarantees convergence: :math:\tau  \in (0, 1/L] where L is
the Lipschitz constant of :math:\nabla f.
niter : :obj:int, optional
Number of iterations of iterative scheme
gfirst : :obj:bool, optional
Apply Proximal of operator g first (True) or Proximal of
operator f first (False)
callback : :obj:callable, optional
Function with signature (callback(x)) to call after each iteration
where x is the current model vector
callbackz : :obj:bool, optional
Modify callback signature to (callback(x, z)) when callbackz=True
show : :obj:bool, optional
Display iterations log

Returns
-------
x : :obj:numpy.ndarray
Inverted model
z : :obj:numpy.ndarray
Inverted second model

--------

Notes
-----
The ADMM algorithm can be expressed by the following recursion:

.. math::

\mathbf{x}^{k+1} = \prox_{\tau f}(\mathbf{z}^{k} - \mathbf{u}^{k})\\
\mathbf{z}^{k+1} = \prox_{\tau g}(\mathbf{x}^{k+1} + \mathbf{u}^{k})\\
\mathbf{u}^{k+1} = \mathbf{u}^{k} + \mathbf{x}^{k+1} - \mathbf{z}^{k+1}

Note that x and z converge to each other, however if iterations are
stopped too early x is guaranteed to belong to the domain of f
while z is guaranteed to belong to the domain of g. Depending on
the problem either of the two may be the best solution.

"""
if show:
tstart = time.time()
'---------------------------------------------------------\n'
'Proximal operator (f): %s\n'
'Proximal operator (g): %s\n'
'tau = %10e\tniter = %d\n' % (type(proxf), type(proxg),
tau, niter))
head = '   Itn       x[0]          f           g       J = f + g'

x = x0.copy()
u = z = np.zeros_like(x)
for iiter in range(niter):
if gfirst:
z = proxg.prox(x + u, tau)
x = proxf.prox(z - u, tau)
else:
x = proxf.prox(z - u, tau)
z = proxg.prox(x + u, tau)
u = u + x - z

# run callback
if callback is not None:
if callbackz:
callback(x, z)
else:
callback(x)
if show:
if iiter < 10 or niter - iiter < 10 or iiter % (niter // 10) == 0:
pf, pg = proxf(x), proxg(x)
msg = '%6g  %12.5e  %10.3e  %10.3e  %10.3e' % \
(iiter + 1, np.real(to_numpy(x[0])),
pf, pg, pf + pg)
print(msg)
if show:
print('\nTotal time (s) = %.2f' % (time.time() - tstart))
print('---------------------------------------------------------\n')
return x, z

[docs]def ADMML2(proxg, Op, b, A, x0, tau, niter=10, callback=None, show=False, **kwargs_solver):
r"""Alternating Direction Method of Multipliers for L2 misfit term

Solves the following minimization problem using Alternating Direction
Method of Multipliers:

.. math::

\mathbf{x},\mathbf{z}  = \argmin_{\mathbf{x},\mathbf{z}}
\frac{1}{2}||\mathbf{Op}\mathbf{x} - \mathbf{b}||_2^2 + g(\mathbf{z}) \\
s.t. \; \mathbf{Ax}=\mathbf{z}

where :math:g(\mathbf{z}) is any convex function that has a known proximal operator.

Parameters
----------
proxg : :obj:pyproximal.ProxOperator
Proximal operator of g function
Op : :obj:pylops.LinearOperator
Linear operator of data misfit term
b : :obj:numpy.ndarray
Data
A : :obj:pylops.LinearOperator
Linear operator of regularization term
x0 : :obj:numpy.ndarray
Initial vector
tau : :obj:float, optional
Positive scalar weight, which should satisfy the following condition
to guarantees convergence: :math:\tau \in (0, 1/\lambda_{max}(\mathbf{A}^H\mathbf{A})].
niter : :obj:int, optional
Number of iterations of iterative scheme
callback : :obj:callable, optional
Function with signature (callback(x)) to call after each iteration
where x is the current model vector
show : :obj:bool, optional
Display iterations log
**kwargs_solver
Arbitrary keyword arguments for :py:func:scipy.sparse.linalg.lsqr used
to solve the x-update

Returns
-------
x : :obj:numpy.ndarray
Inverted model
z : :obj:numpy.ndarray
Inverted second model

--------

Notes
-----
The ADMM algorithm can be expressed by the following recursion:

.. math::

\mathbf{x}^{k+1} = \argmin_{\mathbf{x}} \frac{1}{2}||\mathbf{Op}\mathbf{x}
- \mathbf{b}||_2^2 + \frac{1}{2\tau} ||\mathbf{Ax} - \mathbf{z}^k + \mathbf{u}^k||_2^2\\
\mathbf{z}^{k+1} = \prox_{\tau g}(\mathbf{Ax}^{k+1} + \mathbf{u}^{k})\\
\mathbf{u}^{k+1} = \mathbf{u}^{k} + \mathbf{Ax}^{k+1} - \mathbf{z}^{k+1}

"""
if show:
tstart = time.time()
'---------------------------------------------------------\n'
'Proximal operator (g): %s\n'
'tau = %10e\tniter = %d\n' % (type(proxg), tau, niter))
head = '   Itn       x[0]          f           g       J = f + g'

sqrttau = 1. / sqrt(tau)
x = x0.copy()
u = z = np.zeros(A.shape[0], dtype=A.dtype)
for iiter in range(niter):
# create augumented system
x = regularized_inversion(Op, b, [A, ], x0=x,
dataregs=[z - u, ], epsRs=[sqrttau, ],
**kwargs_solver)[0]
Ax = A @ x
z = proxg.prox(Ax + u, tau)
u = u + Ax - z

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

if show:
if iiter < 10 or niter - iiter < 10 or iiter % (niter // 10) == 0:
pf, pg = 0.5 * np.linalg.norm(Op @ x - b) ** 2, proxg(Ax)
msg = '%6g  %12.5e  %10.3e  %10.3e  %10.3e' % \
(iiter + 1, np.real(to_numpy(x[0])),
pf, pg, pf + pg)
print(msg)
if show:
print('\nTotal time (s) = %.2f' % (time.time() - tstart))
print('---------------------------------------------------------\n')
return x, z

[docs]def LinearizedADMM(proxf, proxg, A, x0, tau, mu, niter=10,
callback=None, show=False):
r"""Linearized Alternating Direction Method of Multipliers

Solves the following minimization problem using Linearized Alternating
Direction Method of Multipliers (also known as Douglas-Rachford splitting):

.. math::

\mathbf{x} = \argmin_\mathbf{x} f(\mathbf{x}) + g(\mathbf{A}\mathbf{x})

where :math:f(\mathbf{x}) and :math:g(\mathbf{x}) are any convex
function that has a known proximal operator and :math:\mathbf{A} is a
linear operator.

Parameters
----------
proxf : :obj:pyproximal.ProxOperator
Proximal operator of f function
proxg : :obj:pyproximal.ProxOperator
Proximal operator of g function
A : :obj:pylops.LinearOperator
Linear operator
x0 : :obj:numpy.ndarray
Initial vector
tau : :obj:float, optional
Positive scalar weight, which should satisfy the following
condition to guarantee convergence: :math:\mu \in (0,
\tau/\lambda_{max}(\mathbf{A}^H\mathbf{A})].
mu : :obj:float, optional
Second positive scalar weight, which should satisfy the following
condition to guarantees convergence: :math:\mu \in (0,
\tau/\lambda_{max}(\mathbf{A}^H\mathbf{A})].
niter : :obj:int, optional
Number of iterations of iterative scheme
callback : :obj:callable, optional
Function with signature (callback(x)) to call after each iteration
where x is the current model vector
show : :obj:bool, optional
Display iterations log

Returns
-------
x : :obj:numpy.ndarray
Inverted model
z : :obj:numpy.ndarray
Inverted second model

--------

Notes
-----
The Linearized-ADMM algorithm can be expressed by the following recursion:

.. math::

\mathbf{x}^{k+1} = \prox_{\mu f}(\mathbf{x}^{k} - \frac{\mu}{\tau}
\mathbf{A}^H(\mathbf{A} \mathbf{x}^k - \mathbf{z}^k + \mathbf{u}^k))\\
\mathbf{z}^{k+1} = \prox_{\tau g}(\mathbf{A} \mathbf{x}^{k+1} +
\mathbf{u}^k)\\
\mathbf{u}^{k+1} = \mathbf{u}^{k} + \mathbf{A}\mathbf{x}^{k+1} -
\mathbf{z}^{k+1}

"""
if show:
tstart = time.time()
'---------------------------------------------------------\n'
'Proximal operator (f): %s\n'
'Proximal operator (g): %s\n'
'Linear operator (A): %s\n'
'tau = %10e\tmu = %10e\tniter = %d\n' % (type(proxf),
type(proxg),
type(A),
tau, mu, niter))
head = '   Itn       x[0]          f           g       J = f + g'
x = x0.copy()
Ax = A.matvec(x)
u = z = np.zeros_like(Ax)
for iiter in range(niter):
x = proxf.prox(x - mu / tau * A.rmatvec(Ax - z + u), mu)
Ax = A.matvec(x)
z = proxg.prox(Ax + u, tau)
u = u + Ax - z

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

if show:
if iiter < 10 or niter - iiter < 10 or iiter % (niter // 10) == 0:
pf, pg = proxf(x), proxg(Ax)
msg = '%6g  %12.5e  %10.3e  %10.3e  %10.3e' % \
(iiter + 1, np.real(to_numpy(x[0])),
pf, pg, pf + pg)
print(msg)
if show:
print('\nTotal time (s) = %.2f' % (time.time() - tstart))
print('---------------------------------------------------------\n')
return x, z

[docs]def TwIST(proxg, A, b, x0, alpha=None, beta=None, eigs=None, niter=10,
callback=None, show=False, returncost=False):
r"""Two-step Iterative Shrinkage/Threshold

Solves the following minimization problem using Two-step Iterative
Shrinkage/Threshold:

.. math::

\mathbf{x} = \argmin_\mathbf{x} \frac{1}{2}
||\mathbf{b} - \mathbf{Ax}||_2^2 + g(\mathbf{x})

where :math:\mathbf{A} is a linear operator and :math:g(\mathbf{x})
is any convex function that has a known proximal operator.

Parameters
----------
proxg : :obj:pyproximal.ProxOperator
Proximal operator of g function
A : :obj:pylops.LinearOperator
Linear operator
b : :obj:numpy.ndarray
Data
x0 : :obj:numpy.ndarray
Initial vector
alpha : :obj:float, optional
Positive scalar weight (if None, estimated based on the
eigenvalues of :math:\mathbf{A}, see Notes for details)
beta : :obj:float, optional
Positive scalar weight (if None, estimated based on the
eigenvalues of :math:\mathbf{A}, see Notes for details)
eigs : :obj:tuple, optional
Largest and smallest eigenvalues of :math:\mathbf{A}^H \mathbf{A}.
If passed, computes alpha and beta based on them.
niter : :obj:int, optional
Number of iterations of iterative scheme
callback : :obj:callable, optional
Function with signature (callback(x)) to call after each iteration
where x is the current model vector
show : :obj:bool, optional
Display iterations log
returncost : :obj:bool, optional
Return cost function

Returns
-------
x : :obj:numpy.ndarray
Inverted model
j : :obj:numpy.ndarray
Cost function

Notes
-----
The TwIST algorithm can be expressed by the following recursion:

.. math::

\mathbf{x}^{k+1} = (1-\alpha) \mathbf{x}^{k-1} +
(\alpha-\beta) \mathbf{x}^k +
\beta \prox_{g} (\mathbf{x}^k + \mathbf{A}^H
(\mathbf{b} - \mathbf{A}\mathbf{x}^k)).

where :math:\mathbf{x}^{1} = \prox_{g} (\mathbf{x}^0 + \mathbf{A}^T
(\mathbf{b} - \mathbf{A}\mathbf{x}^0)).

The optimal weighting parameters :math:\alpha and :math:\beta are
linked to the smallest and largest eigenvalues of
:math:\mathbf{A}^H\mathbf{A} as follows:

.. math::

\alpha = 1 + \rho^2 \\
\beta =\frac{2 \alpha}{\Lambda_{max} + \lambda_{min}}

where :math:\rho=\frac{1-\sqrt{k}}{1+\sqrt{k}} with
:math:k=\frac{\lambda_{min}}{\Lambda_{max}} and
:math:\Lambda_{max}=max(1, \lambda_{max}).

Experimentally, it has been observed that TwIST is robust to the
choice of such parameters. Finally, note that in the case of
:math:\alpha=1 and :math:\beta=1, TwIST is identical to IST.

"""
# define proxf as L2 proximal
proxf = L2(Op=A, b=b)

# find alpha and beta
if alpha is None or beta is None:
if eigs is None:
emin = A.eigs(neigs=1, which='SM')
emax = max([1, A.eigs(neigs=1, which='LM')])
else:
emax, emin = eigs
k = emin / emax
rho =  (1 - sqrt(k)) / (1 + sqrt(k))
alpha = 1 + rho ** 2
beta = 2 * alpha / (emax + emin)

# compute proximal of g on initial guess (x_1)
xold = x0.copy()
x = proxg.prox(xold - proxf.grad(xold), 1.)

if show:
tstart = time.time()
print('TwIST\n'
'---------------------------------------------------------\n'
'Proximal operator (g): %s\n'
'Linear operator (A): %s\n'
'alpha = %10e\tbeta = %10e\tniter = %d\n' % (type(proxg),
type(A),
alpha, beta, niter))
head = '   Itn       x[0]          f           g       J = f + g'

# iterate
j = None
if returncost:
j = np.zeros(niter)
for iiter in range(niter):
# compute new x
xnew = (1 - alpha) * xold + \
(alpha - beta) * x + \
beta * proxg.prox(x - proxf.grad(x), 1.)
# save current x as old (x_i -> x_i-1)
xold = x.copy()
# save new x as current (x_i+1 -> x_i)
x = xnew.copy()

# compute cost function
if returncost:
j[iiter] = proxf(x) + proxg(x)

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

if show:
if iiter < 10 or niter - iiter < 10 or iiter % (niter // 10) == 0:
pf, pg = proxf(x), proxg(x)
msg = '%6g  %12.5e  %10.3e  %10.3e  %10.3e' % \
(iiter + 1, np.real(to_numpy(x[0])),
pf, pg, pf + pg)
print(msg)
if show:
print('\nTotal time (s) = %.2f' % (time.time() - tstart))
print('---------------------------------------------------------\n')
if returncost:
return x, j
else:
return x
`