Source code for pyproximal.optimization.sr3

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

def _lsqr(Op, data, iter_lim, z_old, x0, kappa, eps, Reg):
r"""LSQR

This function uses LSQR to solve the inner iteration of SR3, given by

.. math::
\min_x \dfrac{1}{2}\Vert\begin{bmatrix}A\\ \sqrt{kappa}L\end{bmatrix}x
- \begin{bmatrix}b\\ \sqrt{kappa}w\Vert_2^2

LSQR is stopped when the maximum amount of iterations are reached or when
the stopping criterion is satisfied. There are no other stopping criteria,
hence the solver should not be used outside of the context of solving the
inner iteration of SR3.

Parameters
----------
Op: :obj:pylops.LinearOperator
Forward operator. This is the stacked operator
:math:\begin{bmatrix}A\\ \sqrt{kappa}L\end{bmatrix}
data: :obj:numpy.ndarray
Data
iter_lim: :obj:int
Maximum number of iterations
z_old: :obj:np.ndarray
The previous outer iteration
x0: :obj:numpy.ndarray
initial guess
kappa: :obj:float
The regularization parameter for the inner iteration
eps: :obj:float
The regularization parameter for the outer iteration
Reg: :obj:np.ndarray
The regularization operator L

Returns
-------
x: :obj:numpy.ndarray
Approximate solution

"""
data -= Op.matvec(x0)
x = x0
beta = np.linalg.norm(data)
u = data / beta
v = Op.rmatvec(u)
alpha = np.linalg.norm(v)
v = v / alpha
w = v
phi_bar = beta
rho_bar = alpha
for _ in range(iter_lim):
u = Op.matvec(v) - alpha*u
beta = np.linalg.norm(u)
u = u / beta
v = Op.rmatvec(u) - beta*v
alpha = np.linalg.norm(v)
v = v / alpha
rho = np.sqrt(rho_bar**2 + beta**2)
c = rho_bar/rho
s = beta/rho
theta = s*alpha
rho_bar = -c*alpha
phi = c*phi_bar
phi_bar = s*phi_bar
x += (phi/rho)*w
temp = Reg.matvec(x)
z = np.sign(temp) * np.maximum(abs(temp) - eps/kappa, 0)
if np.linalg.norm(z - z_old) < 1e-6 * np.linalg.norm(z_old):
return x
w = v - (theta/rho)*w
z_old = z
return x

[docs]def SR3(Op, Reg, data, kappa, eps, x0=None, adaptive=True,
iter_lim_outer=int(1e2), iter_lim_inner=int(1e2)):
r"""Sparse Relaxed Regularized Regression

Applies the Sparse Relaxed Regularized Regression (SR3) algorithm to
an inverse problem with a sparsity constraint of the form

.. math::

\min_x \dfrac{1}{2}\Vert \mathbf{Ax} - \mathbf{b}\Vert_2^2 +
\lambda\Vert \mathbf{Lx}\Vert_1

SR3 introduces an auxiliary variable :math:\mathbf{z} = \mathbf{Lx},

.. math::

\min_{\mathbf{x},\mathbf{z}} \dfrac{1}{2}\Vert \mathbf{Ax} -
\mathbf{b}\Vert_2^2 + \lambda\Vert \mathbf{z}\Vert_1 +
\dfrac{\kappa}{2}\Vert \mathbf{Lx} - \mathbf{z}\Vert_2^2

Parameters
----------
Op : :obj:pylops.LinearOperator
Forward operator
Reg : :obj:numpy.ndarray
Regularization operator
data : :obj:numpy.ndarray
Data
kappa : :obj:float
Parameter controlling the difference between :math:\mathbf{z}
and :math:\mathbf{Lx}
eps : :obj:float
Regularization parameter
x0 : :obj:numpy.ndarray, optional
Initial guess
adaptive : :obj:bool, optional
Use adaptive SR3 with a stopping criterion for the inner iterations
or not
iter_lim_outer : :obj:int, optional
Maximum number of iterations for the outer iteration
iter_lim_inner : :obj:int, optional
Maximum number of iterations for the inner iteration

Returns
-------
x: :obj:numpy.ndarray
Approximate solution.

Notes
-----
SR3 uses the following algorithm:

.. math::
\mathbf{x}_{k+1} = (\mathbf{A}^T\mathbf{A} + \kappa
\mathbf{L}^T\mathbf{L})^{-1}(\mathbf{A}^T\mathbf{b} +
\kappa \mathbf{L}^T\mathbf{y}_k) \\
\mathbf{y}_{k+1} = \prox_{\lambda/\kappa\mathcal{R}}
(\mathbf{Lx}_{k+1})

"""
(m, n) = Op.shape
if x0 is None:
x0 = np.zeros(n)
x = x0.copy()
p = Reg.shape[0]
v = np.zeros(p)
w = v
eta = 1/kappa
theta = 1
AL = pylops.VStack([Op, np.sqrt(kappa) * Reg])
for _ in range(iter_lim_outer):
# Compute the inner iteration
x = _lsqr(AL, np.concatenate((data, np.sqrt(kappa)*v)),
iter_lim_inner, v, x, kappa, eps, Reg)
else:
x = sp_lsqr(AL, np.concatenate((data, np.sqrt(kappa)*v)),
iter_lim=iter_lim_inner, x0=x)[0]
# Compute the outer iteration
w_old = w
temp = Reg.matvec(x)
w = np.sign(temp) * np.maximum(abs(temp) - eta*eps, 0)
err1 = np.linalg.norm(v - w) / max(1, np.linalg.norm(w))
if err1 < 1e-6:
return x
theta = 2/(1 + np.sqrt(1+4/(theta**2)))
v = w + (1-theta)*(w - w_old)
return x