# Interpretation of Plug-and-Play (PnP) algorithms from a different angle

Abinash Nayak  Visiting Assistant Professor, Department of Mathematics,
University of Alabama at Birmingham, University Hall, Room 4005, 1402 10th
Avenue South, Birmingham AL 35294-1241, (p) 205.934.2154, (f) 205.934.9025
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

It’s well-known that inverse problems are ill-posed and to solve them
meaningfully, one has to employ regularization methods. Traditionally, the
most popular regularization approaches are Variational-type approaches, i.e.,
penalized/constrained functional minimization. In recent years, the classical
regularization approaches have been replaced by the so called plug-and-play
(PnP) algorithms, which copies the proximal gradient minimization processes,
such as ADMM or FISTA, but with any general denoiser. However, unlike the
traditional proximal gradient methods, the theoretical underpinnings and
convergence results have been insufficient for these PnP-algorithms. Hence,
the results from these algorithms, though empirically outstanding, are not
well-defined, in the sense of, being a minimizer of a Variational problem, or
in some other forms. In this paper, we address this question of “well-
definedness”, but from a different angle. We explain these algorithms from the
viewpoint of a semi-iterative regularization method. In addition, we expand
the family of regularized solutions, corresponding to the classical semi-
iterative methods, to a much larger class, which encompasses these algorithms,
as well as, enhance the recovery process. We conclude with several numerical
results which validate the developed theories and reflect the improvements
over the traditional PnP-algorithms, such as ADMM-PnP and FISTA-PnP.

###### Key words and phrases:

Inverse problems, Ill-posed problems, Regularization, Variational
minimization, Numerical methods, Plug-and-Play (PnP), BM3D denoiser, Computed
tomography

###### 1991 Mathematics Subject Classification:

Primary 65K05, 65K10; Secondary 65R30, 65R32

##  1\. Introduction

###  1.1. Inverse Problems and Regularization:

Mathematically, an inverse problem is often expressed as the problem of
estimating a (source) $\hat{x}$ which satisfies, for a given (effect) $b$, the
following matrix (or operator) equation

$A\hat{x}=b,$ (1.1)

where the matrix $A\in\mathbb{R}^{n\times m}$ and the vectors
$\hat{x}\in\mathbb{R}^{n}$, $b\in\mathbb{R}^{m}$ are the discrete
approximations of an infinite dimensional model describing the underlying
physical phenomenon. The inverse problem (1.1) is usually ill-posed, in the
sense of violating any one of the Hadamard’s conditions for well-posedness:
(i) Existence of a solution (for $b$ not in the range of $A$), (ii) uniqueness
of the solution (for non-trivial null-space of $A$) and, (iii) continuous
dependence on the data (for ill-conditioned A). Conditions (i) and (ii) can be
circumvented by relaxing the definition of a solution for (1.1), for example,
finding the least square solution or the minimal norm solution (i.e., the
pseudo-inverse solution $x^{\dagger}$). The most (practically) significant
condition is condition-(iii), since failing of this leads to an absurd
(unstable) solution. That is, for an (injective) A and an exact $b$
(noiseless), the solution of (1.1) can be approximated by the (LS) least-
square solution ($x^{\dagger}$), i.e., $x^{\dagger}$ is the minimizer of the
following least-square functional

$F(x)=||Ax-b||_{2}^{2},$ (1.2)

where as, for a noisy data $b_{\delta}$ (which is practically true all the
time) such that $||b-b_{\delta}||\leq\delta$, the simple least-square solution
($x^{\dagger}_{\delta}$) with respect to $b_{\delta}$ in (1.2) fails to
approximate the true solution, i.e.,
$||x^{\dagger}_{\delta}-x^{\dagger}||>>\delta$, which in turn implies,
$||x^{\dagger}_{\delta}-\hat{x}||>>\delta$, due to the ill-posedness of the
inverse problem (1.1). To counter such instabilities or ill-posedness of
inverse problems, regularization methods have to be employed, which are
broadly divided into two types.

###  1.2. Variational (or penalized) regularization

Such approaches, also known as Tikhonov-type regularization, are probably the
most well known regularization techniques for solving linear as well as
nonlinear inverse problems (see [1, 2, 3, 4, 5]), where, instead of minimizing
the simple least-square functional (1.2), one minimizes a penalized (or
constrained) functional:

$F(x;\mathcal{D},\lambda,\mathcal{R})=\mathcal{D}(Ax,b_{\delta})+\lambda\mathcal{R}(x),$
(1.3)

where $\mathcal{D}$ is called the data-fidelity term (imposing data-
consistency), $\mathcal{R}$ is the regularization term (imposing certain
structures, based on some prior knowledge of the solution $\hat{x}$) and,
$\lambda\geq 0$ is the regularization parameter that balances the trade-off
between them, depending on the noise level $\delta$, i.e.,
$\lambda=\lambda(\delta)$. The formulation (1.3) also has a Bayesian
interpretation, where the minimization of
$F(x;\mathcal{D},\lambda,\mathcal{R})$ corresponds to the maximum-a-posteriori
(MAP) estimate of $\hat{x}$ given $b_{\delta}$, where the likelihood of
$b_{\delta}$ is proportional to $exp(-\mathcal{D}(x))$ and the prior
distribution on $\hat{x}$ is proportional to $exp(-\mathcal{R}(x))$.
Classically, $\mathcal{D}(Ax,b_{\delta})=||{Ax-b_{\delta}}||_{p}^{p}$ and
$\mathcal{R}(x)=||Lx-\bar{x}||_{q}^{q}$, where $L$ is a regularization matrix
with the null spaces of $A$ and $L$ intersecting trivially, and $p$, $q$
determine the involved norms. For large scale problems, the minimization of
(1.3) is done iteratively, i.e., starting from an initial guess
$x_{0}^{\delta}$ and step-sizes $\tau_{k}>0$, the minimizer of (1.3) is
approximated via the sequence

$\displaystyle
x_{k+1}^{\delta}=x_{k}^{\delta}-\tau_{k}\nabla_{x}F(x_{k}^{\delta};\mathcal{D},\lambda,\mathcal{R}),$
(1.4)

Hence, one has to address the following two technicalities,

  1. (1)

Compute, if possible, the gradient
$\nabla_{x}F(x_{k}^{\delta};\mathcal{D},\lambda,\mathcal{R})$ at each k. It’s
easier when both $\mathcal{D}(x)$ and $\mathcal{R}(x)$ are differentiable, as
then,
$\nabla_{x}F(x_{k}^{\delta};\mathcal{D}(x_{k}^{\delta}),\lambda,\mathcal{R})=\nabla_{x}\mathcal{D}(x_{k}^{\delta})+\lambda\nabla_{x}\mathcal{R}(x_{k}^{\delta})$.
However, when they are non-differentiable, then it’s not that straight
forward. In such scenarios, for convex $\mathcal{R}$, one has to make use of
convex optimization techniques, such as proximal gradient methods, to
circumvent the differentiability issue, which is discussed in later sections.

  2. (2)

Convergence of the sequence of iterates $\\{x_{k}^{\delta}\\}$ as
$k\rightarrow\infty$. For strictly convex $\mathcal{D}$ and $\mathcal{R}$, we
have $x_{k}^{\delta}\xrightarrow{k\rightarrow\infty}x^{\delta}_{*}$, where
$x^{\delta}_{*}$ is the global minimizer of (1.3). Again, if they are non-
convex (which is the case for most modern $\mathcal{R}$), then one has to
analyze the behaviour of the sequence when $k\rightarrow\infty$, for example,
convergence to a local minimizer or a saddle-point etc.

###  1.3. Related Works

The above two queries have been studied extensively for classical
regularizers, such as TV or sparsifying wavelet transformation, and the
associated results can be found in the literature, for example in [6, 7, 8, 9,
10]. Typically, for a non-smooth $\mathcal{R}$, which is proper, closed, and
convex, the differentiability issue is circumvented by using a proximal
operator (see [11, 12, 13], and references therein), given by

$\mbox{Prox}_{\lambda\mathcal{R}}(v)=\operatorname*{arg\,min}_{x}\;\lambda\mathcal{R}(x)+\frac{1}{2}||x-v||_{2}^{2}.$
(1.5)

Basically, for smooth $\mathcal{D}(x)$ and non-smooth $\mathcal{R}(x)$, the
minimization problem corresponding to (1.3) can be solved via two first-order
iterative methods:

  1. (1)

Forward-backward splitting (FBS), also known as Iterative shrinkage/soft
thresholding algorithm (ISTA) and has a faster variant Fast ISTA (FISTA),
where each minimization step is divided into two sub-steps, given by

$\displaystyle z_{k+1}^{\delta}$
$\displaystyle=x_{k}^{\delta}-\tau_{k}\nabla_{x}\mathcal{D}(x_{k}^{\delta})\;\;\longleftarrow\text{
data-consistency step}$ (1.6) $\displaystyle x_{k+1}^{\delta}$
$\displaystyle=\text{Prox}_{\lambda\tau_{k}\mathcal{R}}(z_{k+1}^{\delta})\;\;\longleftarrow\text{
data-denoising step}$ (1.7)

where $\tau_{k}\geq 0$ is the step-size at the kth iteration.

  2. (2)

Alternating direction method of multipliers (ADMM), where three sequences are
alternatively updated as follows,

$\displaystyle x_{k+1}^{\delta}$
$\displaystyle=\text{Prox}_{\frac{1}{\rho}\mathcal{D}}(z_{k}^{\delta}-u_{k}^{\delta})\;\;\longleftarrow\text{
data-consistency step}$ (1.8) $\displaystyle z_{k+1}^{\delta}$
$\displaystyle=\text{Prox}_{\frac{\lambda}{\rho}\mathcal{R}}(x_{k+1}^{\delta}+u_{k}^{\delta})\;\;\longleftarrow\text{
data-denoising step}$ (1.9) $\displaystyle u_{k+1}^{\delta}$
$\displaystyle=u_{k}^{\delta}+x_{k+1}^{\delta}-z_{k+1}^{\delta},\;\;\longleftarrow\text{
noise update step}$ (1.10)

where $\rho>0$ is the Lagrangian parameter, which only effects the speed of
convergence and not the solution (minimizer) of (1.3).

From the above two expressions, one can observe that, each method comprises of
two fundamental steps: (1) data-consistency and (2) data-denoising. This
motivated, authors in [14], to replace the ProxσR operator in the denoising
step of ADMM by an off-the-shelf denoiser $H_{\sigma}$, with denoising
strength corresponding to noise level $\sigma$, and termed it as the PnP-
algorithm (plug-and-play method). However, note that, once the proximal
operator is replaced by any general denoiser, then the Variational problem
(1.3) breaks down, as not all denoisers can be expressed as a proximal
operator of some function $\mathcal{R}$. Hence, certain natural questions
arise, such as, where does this new sequence of iterates converge (if they
do)? And/or, how to assign a meaning to the recovered solution?

Though empirical results show the convergence of these PnP-algorithms, there
is no proof of it, for any general denoisers. However, under certain
assumptions and restrictions (such as boundedness, nonexpansiveness, etc.) on
the denoiser, there have been some convergence proof, see [15, 16, 17, 18, 19]
and references therein. There are also some other variants of such PnP-
methods, such as Regularization by Denoising (RED)[20], Regularization by
Artifact-Removal (RARE) [21], etc.

### Contribution of this paper

  * •

In this paper, we try to answer the above questions, the convergence of the
iterates and the meaning of the solution of a PnP-method, from a different
angle. We show a connection between these PnP-algorithms and semi-iterative
regularization methods. We then present a bound for the difference of a “PnP-
solution” and a “regularized LS-solution”, where the definitions of these
terms are explained below.

  * •

PnP-algorithms, like Variational regularization, are very sensitive to the
regularization/denoising parameter. Here, we present a method to significantly
reduce the sensitivity of these algorithms to the denoiser’s strength of
denoising. This is crucial, since this can transform a bad denoiser (denoisers
with too strong or too weak denoising) to a good one (an appropriate level of
denoising), without actually altering the involved denoiser.

  * •

We also address the question of when to stop/terminate the iterations to
recover “an appropriate solution”. Again, the meaning of an appropriate
solution is described below, which depends on some selection criteria.

  * •

In addition, we compare the FBS-PnP algorithm with the ADMM-PnP algorithm and
point out some of the pros and cons of these algorithms. Note that, although
in the traditional setting (1.3) both these algorithms produce a similar
result, the minimizer of (1.3), but for their PnP versions they may differ
depending on the scenarios, as the convergence is not known.

  * •

Finally, we conclude with presenting several numerical results in support of
the theories developed in this paper, which also validates the efficiency and
effectiveness of the modified algorithms suggested here. We perform
computational experiments on the X-ray computed tomography problems, and
compare the results obtained using our approach with the traditional ADMM-PnP
and FBS-PnP algorithms, with BM3D as the denoiser.

##  2\. Structure imposing descent directions

In the previous section we explained Variational or Tikhonov-type
regularization methods and their variants. The second most popular
regularization technique is known as Semi-iterative regularization method.

###  2.1. Semi-iterative regularization

Landweber, in [22], showed that minimizing the simple LS functional (1.2) via
the gradient direction, with a constant step-size, leads to the recovery of
the solution of (1.1) in a semi-convergent manner, i.e., starting from an
initial guess $x_{0}^{\delta}$, the sequence of iterates

$x_{k}^{\delta}=x_{k-1}^{\delta}+d_{k-1}^{\delta},$ (2.1)

where $d_{k-1}^{\delta}:=-\tau\nabla_{x}\mathcal{D}(x_{k-1}^{\delta})=-\tau
A^{*}(Ax_{k-1}^{\delta}-b_{\delta})$, when $\mathcal{D}(x)=||Ax-
b_{\delta}||_{2}^{2}$, for $0\leq\tau<\frac{1}{2||A||_{2}^{2}}$, initially
converges towards the true LS-solution $x^{\dagger}$, then diverges away from
$x^{\dagger}$, i.e., the recovery-errors follow a semi-convergent nature. That
is, the iteration index $k$ plays the role of a regularization parameter,
where for smaller values of $k$ the recovered solution $x_{k}^{\delta}$ is
over-regularized (over-smoothed) and for larger values of k, $x_{k}^{\delta}$
is under-regularized (over-fitted). Therefore, for “an appropriate” index
$k(\delta)$, one would obtain a regularized solution $x_{k(\delta)}^{\delta}$
which can approximate the true solution of (1.1), where an appropriate index
depends on the choices of the selection criteria, this is discussed in details
in successive sections. Also, instead of a constant step-size in (2.1), one
can have varying step-sizes $\tau_{k}\geq 0$ for each iteration and can even
generalize the simple gradient descent method to much faster methods, such as
Krylov-Subspace or Conjugate-gradient semi-iterative methods, for details and
generalizations see [23, 22, 24, 1]. However, as one can notice, through this
approach (or its extensions, which one improves the speed of convergence) one
cannot make use of any prior knowledge of the solution $\hat{x}$, i.e., no
structures can be imposed on the recovery process.

###  2.2. Regularized solution vs. Regularized solution family

Before we jump into our interpretation of the PnP-algorithms, we make a quick
inspection on the pros and cons of the aforementioned regularization methods,
as this will help in better understanding of our formulation. The similarity
between both these methods is that, first, one generates a “family of
regularized-solutions” (depending on a parameter) and then, an “appropriate
regularized solution ($x^{\delta}$)” is selected form that class, based on
some a-priori or a-posteriori selection rules or criteria, depending on the
noise-level $\delta$. In the dissimilarity, the constraint-regularization
offers a much larger set of regularized solutions which depends not only on
the choices of $\lambda$ but also on $\mathcal{R}$ (the
constraint/regularizer), i.e., the set (family) of regularized solutions in
Tikhonov-type methods is

$\mathcal{T}_{1}:=\left\\{x^{\delta}(\mathcal{D},\lambda,\mathcal{R}):x^{\delta}(\mathcal{D},\lambda,\mathcal{R})\mbox{
is a minimizer of \eqref{Gen. Tik. fun.}, }\lambda=\lambda(\delta)\geq
0\right\\},$ (2.2)

where as, the classical (simple) semi-iterative regularization can only
provide a set of regularized solutions depending on the iteration index $k$,
i.e.,

$\mathcal{I}_{2}:=\left\\{x^{\delta}(\mathcal{D},k):x_{k}^{\delta}=x_{k-1}^{\delta}-\tau_{k-1}\nabla_{x}\mathcal{D}(x_{k-1}^{\delta}),\;\;1\leq
k\leq k(\delta)\leq\infty\right\\}.$ (2.3)

However, on the plus side, a regularized solution in the family (2.3) can be
obtained much faster (though without any structures in it) than a regularized
solution in (2.2), as in $\mathcal{T}_{1}$ one has to minimize completely the
associated functional in (1.3).

######  Remark 2.1.

Note that, the formulations (1.3) and (2.1) only yield corresponding families
of regularized solutions $\mathcal{T}_{1}$ and $\mathcal{I}_{2}$,
respectively, and not a solution of (1.1). One estimate the true solution
$\hat{x}$ through these families of solutions depending on some “parameter
choice or selection criteria”. That is, based on a specific criterion, say
$\mathcal{S}_{0}$, one chooses a regularized solution from either the family
(2.2),
$x_{0}^{\delta}(\mathcal{D},\lambda(\delta),\mathcal{R})\in\mathcal{T}_{1}$,
such that $\delta\rightarrow 0$ implies
$\lambda(\delta;\mathcal{S}_{0})\rightarrow 0$ and
$x_{0}^{\delta}(\mathcal{D},\lambda(\delta),\mathcal{R})\rightarrow\hat{x}$,
or from the family (2.3),
$x_{0}^{\delta}(\mathcal{D},k(\delta))\in\mathcal{I}_{1}$, such that
$\delta\rightarrow 0$ implies $k(\delta;\mathcal{S}_{0})\rightarrow\infty$ and
$x_{0}^{\delta}(\mathcal{D},k(\delta))\rightarrow\hat{x}$. Of course, when the
selection criterion changes, the regularized solution changes too and hence,
the families (2.2) and (2.3) of regularized solutions get extended by the
parameter choice criterion, i.e.,

$\displaystyle\mathcal{T}_{2}:=\left\\{x^{\delta}(\mathcal{D},\lambda,\mathcal{R};\mathcal{S})\right.:$
$\displaystyle\left.x^{\delta}(\mathcal{D},\lambda,\mathcal{R})\mbox{ is a
minimizer of \eqref{Gen. Tik. fun.}, }\lambda=\lambda(\delta,\mathcal{S})\geq
0,\right.$ $\displaystyle\left.\mbox{and $\mathcal{S}$ is a selection
criterion}\right\\},$ (2.4)
$\displaystyle\mathcal{I}_{2}:=\left\\{x^{\delta}(\mathcal{D},k;\mathcal{S})\right.\;:$
$\displaystyle\;\left.x_{k}^{\delta}=x_{k-1}^{\delta}-\tau_{k-1}\nabla_{x}\mathcal{D}(x_{k-1}^{\delta}),\;\;1\leq
k\leq k(\delta,\mathcal{S})\leq\infty,\right.$ $\displaystyle\left.\mbox{ and
$\mathcal{S}$ is a selection criterion}\right\\}.$ (2.5)

Note that, for practical problems, it’s near impossible to obtain the optimal
parameters for the respective families, since the true solution $\hat{x}$ is
unknown. Although the problem of choosing an appropriate selection criterion
is non-trivial, especially for large scale problems with unknown error, there
has been a lot of studies in the literature in this regard, for example Cross-
Validation (CV), Generalized Cross-Validation (GCV), Discrepancy principle
(DP), L-curve, monotone error rule etc., see [25, 26, 27, 28, 29, 30, 31, 32].

######  Remark 2.2.

Now, one may be tempted to fuse the above two families naively, for example,
one can incorporate the iteration index $k$ as an additional regularization
parameter in (2.1) and can define a regularized solution,
$x^{\delta}(\mathcal{D},\lambda,\mathcal{R},k;\mathcal{S})$, as the k-th
iterate $x_{k}^{\delta}(\mathcal{D},\lambda,\mathcal{R})$ when minimizing
(1.3), as defined in (1.4), which satisfies the selection criterion
$\mathcal{S}$, for fixed $(\mathcal{D},\lambda,\mathcal{R};\mathcal{S})$.
However, under this definition, the new regularized solution
$x_{k}^{\delta}(\mathcal{D},\lambda,\mathcal{R})$ is neither a solution
(minimizer) corresponding to a Variational problem (1.3), nor a solution
corresponding to any classical semi-iterative regularization method, as the
gradient $\nabla_{x}F(x_{k}^{\delta};\mathcal{D},\lambda,\mathcal{R})$,
defined in (1.4) for all $k$, may not be a descent direction for minimizing
the LS-functional, defined in (1.2). That is,
$x^{\delta}(\mathcal{D},\lambda,\mathcal{R},k;\mathcal{S})$ does not belong to
either of the two classes, $\mathcal{T}_{2}$ or $\mathcal{I}_{2}$, and hence,
we need to further expand these classes. Moreover, this definition does not
work for PnP-algorithms, for any general denoiser, as not all PnP-algorithms
can be formulated as a Variational problem. Nevertheless, the advantage of
such a definition is that, one doesn’t have to worry about the convergence of
the iterates $x_{k}^{\delta}(\mathcal{D},\lambda,\mathcal{R})$, since, if it
does not belong to the family $\mathcal{I}_{3}$ then we are not looking for
the solution of (1.3) anymore. In addition, this also speeds up the recovery
process, as one does not have to iterate indefinitely for convergence. This
make sense, to a certain extent, since we are interested in the solution of
(1.1) and the formulation (1.3) only aids us to approximate the true solution
$\hat{x}$, through the regularization parameter $\lambda(\delta,\mathcal{S})$.
Hence, even when the iterates (1.4) converge (for fixed
$(\mathcal{D},\lambda,\mathcal{R})$) to $x^{\delta}_{\lambda}$, a minimizer of
(1.3), $x^{\delta}_{\lambda}$ is not the solution of (1.1), rather, one
typically generates a list of $\\{x^{\delta}_{\lambda}\\}_{\lambda}$ and an
appropriate estimate of $\hat{x}$ is then selected based on the selection
criterion $\mathcal{S}$. This naive idea lays the ground for our formulation
to interpret the PnP-algorithms, but has to be polished properly for
meaningful definitions, which is done in the following sections.

###  2.3. Structured iterations

Here, instead of expanding the family $\mathcal{T}_{2}$ to incorporate the
PnP-algorithms, we expand the family $\mathcal{I}_{2}$ to address the above
problems. First, note that, when minimizing the LS-functional (1.2),
$d_{k}^{\delta}$ in (2.1) needs to be a descent direction only and doesn’t
have to be the negative gradient (which is the steepest descent direction).
Hence, any direction with the following property can be considered as a
descent direction, for simplicity we choose $\mathcal{D}(x)=||Ax-
b_{\delta}||_{2}^{2}$,

$\Big{(}d_{k}^{\delta},-\nabla_{x}\mathcal{D}(x_{k}^{\delta})\Big{)}_{2}=\Big{(}d_{k}^{\delta},-A^{*}(Ax_{k}^{\delta}-b_{\delta})\Big{)}_{2}>0,$
(2.6)

where $\Big{(}.,.\Big{)}_{2}$ is the associated $\ell_{2}$-product. It’s also
easy to verify that if

$||e_{k}^{\delta}+\tau\nabla_{x}\mathcal{D}(x_{k}^{\delta})||_{2}^{2}\;\leq\;\epsilon_{k}(\delta),$
(2.7)

for some $\epsilon_{k}(\delta)\geq 0$, then, starting from initial points
$z_{0}^{\delta}$ and $x_{0}^{\delta}$,

$\displaystyle||z_{k}^{\delta}-x_{k}^{\delta}||_{2}^{2}\;$
$\displaystyle\leq\;||z_{0}^{\delta}-x_{0}^{\delta}||_{2}^{2}\;+\;\sum_{i=0}^{k-1}||e_{i}^{\delta}+\tau\nabla_{x}\mathcal{D}(x_{i}^{\delta})||_{2}^{2}\;\leq\;||z_{0}^{\delta}-x_{0}^{\delta}||_{2}^{2}\;+\;\sum_{i=0}^{k-1}\epsilon_{i}(\delta)$
(2.8)

where $z_{i}^{\delta}=z_{i-1}^{\delta}+e_{i}^{\delta}$ and
$x_{i}^{\delta}=x_{i-1}^{\delta}-\tau\nabla_{x}\mathcal{D}(x_{i-1}^{\delta})$.
Hence, if for some $\epsilon(\delta)\geq 0$,
$||z_{0}^{\delta}-x_{0}^{\delta}||_{2}^{2}\leq\epsilon(\delta)$ and
$\epsilon_{i}(\delta)\leq\epsilon(\delta)$, for $0\leq i\leq k-1$, then we
have

$||z_{k}^{\delta}-x_{k}^{\delta}||_{2}^{2}\;\leq\;C(k)\;\epsilon(\delta),$
(2.9)

where the (increasing) constant $C(k)$ depends only on the iteration index
$k$. In other words, if any direction $d_{k}^{\delta}$ satisfies conditions
(2.6) and (2.7), then the iterates $z_{k}^{\delta}$, formed based on it, also
remains close to the LS iterates $x_{k}^{\delta}$, upto some constant
depending on the iteration index $k$. And, since in a semi-iterative
regularization $k(\delta)$ is bounded away from infinity, the new iterates
$z_{k}^{\delta}$ are bounded at most $C(k(\delta))\epsilon(\delta)$ away from
the LS iterates $x_{k}^{\delta}$. It is also easy to check that, for the
residues or discrepancy terms $\mathcal{D}(x)=||Ax-b_{\delta}||_{2}^{2}$, we
have

$\displaystyle||Az_{k}^{\delta}-b_{\delta}||_{2}^{2}\leq\;||Ax_{k}^{\delta}-b_{\delta}||_{2}^{2}\;+\;C(k)\epsilon(\delta)||A||_{2}^{2},$
(2.10)

that is, if the residues of the LS iterates are small
($||Ax_{k}^{\delta}-b_{\delta}||_{2}^{2}\approx\delta$), then so are the
residues of $z_{k}^{\delta}$ with an error term
$C(k)\epsilon(\delta)||A||_{2}^{2}$, where $A$ is usually bounded. This is
helpful when the selection criterion, such as discrepancy principle, depends
on the residue norms, i.e., terminate the iterations when the residue norms of
the iterates $x_{k}^{\delta}$ is around certain limit
($||Ax_{k}^{\delta}-b_{\delta}||_{2}^{2}\approx\delta$), and (2.10) implies
that for those iterations the residue norms corresponding to the new iterates
$z_{k}^{\delta}$ is also around that limit, with an error term, implying that
the iterates $z_{k}^{\delta}$ have also fitted the data to that extent. Hence,
when minimizing the LS-functional (1.2), replacing the vanilla descent
directions ($-\nabla_{x}D(x_{k}^{\delta})$) with an “appropriate” descent
direction $d_{k}^{\delta}$, one can impose certain structures in the recovery
process through the iterates $z_{k}^{\delta}$, and thus, can recover a
regularized structured solution $z^{\delta}_{k(\delta)}$. With this
formulation, the family of regularized solution $\mathcal{I}_{2}$, as defined
in (2.1), is expanded to

$\displaystyle\mathcal{I}_{3}:=$
$\displaystyle\left\\{x^{\delta}(\mathcal{D},k,d_{k}^{\delta};\mathcal{S})\right.\;:\;\left.x_{k}^{\delta}=x_{k-1}^{\delta}\;+\;d_{k-1}^{\delta},\;\;1\leq
k\leq k(\delta,\mathcal{S})\leq\infty,\right.$ (2.11)
$\displaystyle\left.\;\mbox{s.t.}\;\Big{(}d_{k-1}^{\delta},-\tau_{k-1}\nabla_{x}\mathcal{D}(x_{k-1}^{\delta})\Big{)}_{2}>0\mbox{
and $\mathcal{S}$ is a selection criterion.}\right\\}.$

Figure 1 elaborates the families of solutions, $\mathcal{I}_{3}$ and
$\mathcal{T}_{2}$, visually. Note that, although
$\mathcal{I}_{3}\supseteq\mathcal{I}_{2}$, $\mathcal{I}_{3}$ is neither a
super-set, nor a subset, nor equals to $\mathcal{T}_{2}$. However, for a fixed
$\mathcal{S}_{0}$, any
$x^{\delta}(\mathcal{D},\lambda(\delta),\mathcal{R})\in\mathcal{T}_{2}$ can be
“approximated” by a
$x^{\delta}(\mathcal{D},k(\delta),d_{k}^{\delta})\in\mathcal{I}_{3}$, where
the extent of approximation depends on the choice of $\lambda(\delta)$, since,
while minimizing (1.3) using the FBS approach (1.6), the resulting direction
at the k-1 step is given by

$\displaystyle d_{k-1}^{\delta}$
$\displaystyle=\mbox{Prox}_{\lambda\mathcal{R}}(x_{k}^{\delta})\;-\;x_{k-1}^{\delta}$
(2.12)
$\displaystyle=\left(\mbox{Prox}_{\lambda\mathcal{R}}(x_{k}^{\delta})-x_{k}^{\delta}\right)+\left(-\tau\nabla_{x}\mathcal{D}(x_{k-1}^{\delta})\right),$

and hence, if

$||\mbox{Prox}_{\lambda\mathcal{R}}(x_{k}^{\delta})-x_{k}^{\delta}||_{2}\;<\;||-\tau\nabla_{x}\mathcal{D}(x_{k-1}^{\delta})||_{2},$
(2.13)

we have

$\displaystyle\Big{(}d_{k-1}^{\delta},-\tau\nabla_{x}\mathcal{D}(x_{k-1}^{\delta})\Big{)}_{2}\geq\tau||\nabla_{x}\mathcal{D}(x_{k-1}^{\delta})||\left(\tau||\nabla_{x}\mathcal{D}(x_{k-1}^{\delta})||-||\mbox{Prox}_{\lambda\mathcal{R}}(x_{k}^{\delta})-x_{k}^{\delta}||\right)>0,$

i.e., $d_{k-1}^{\delta}$ is a descent direction if (2.13) holds. But, as the
residues $\mathcal{D}(x_{k}^{\delta})$ get smaller, so does the gradient
$-\tau\ \nabla_{x}\mathcal{D}(x_{k}^{\delta})$ for a fixed $\tau$, and hence,
(2.13) won’t hold for such iterates, as well as, for large values of
$\lambda$. Therefore, for a given $\lambda(\delta)$, the largest iterate
$k(\delta)$ for which (2.13) holds, yields a regularized solution
$x^{\delta}(\mathcal{D},k(\delta),d_{k}^{\delta})\in\mathcal{I}_{3}$ which can
be close to
$x^{\delta}(\mathcal{D},\lambda(\delta),\mathcal{R})\in\mathcal{T}_{2}$,
depending on $\lambda(\delta)$.

(a) Semi-iterative regularization for different descent directions
$d_{k}^{\delta}$.

(b) Tikhonov-type regularization for different regularizers $\mathcal{R}$.

Figure 1. Semi-iterative vs. Tikhonov-type regularization.

######  Remark 2.3.

Note that, if a direction $d_{k}^{\delta}$ satisfies (2.6) for all $k$, then
the corresponding recovery-errors for the iterates $z_{k}^{\delta}$ will also
follow a semi-convergent nature. However, the parameter choice value
$k(\delta)$, for a fixed selection criterion $\mathcal{S}_{0}$, will be
different for the respective iterates. Also, equation (2.9) does not suggest
that
$||z_{k(\delta)}^{\delta}-\hat{x}||\leq||x_{k(\delta)}^{\delta}-\hat{x}||$,
where $\hat{x}$ is the true solution, i.e., $z_{k(\delta)}^{\delta}$ is a
better estimate of $\hat{x}$ than $x_{k}^{\delta}$, but, if we believe that
$d_{k}^{\delta}$ is a “better” direction than
$-\tau\nabla_{x}\mathcal{D}(x_{k}^{\delta})$, then we would expect that the
iterates $z_{k}^{\delta}$ to approximate the true solution better than
$x_{k}^{\delta}$. In any case, in the absence of the true solution, a “better
solution” for a fixed $\mathcal{S}_{0}$ is the one that satisfies
$\mathcal{S}_{0}$ better.

Moreover, for $z_{k(\delta)}^{\delta}$ to be a regularized solution,
$z_{k(\delta)}^{\delta}\rightarrow\hat{x}$ when $\delta\rightarrow 0$, and
this is true from (2.9) provided $C(k)\epsilon(\delta)\rightarrow 0$. That is,
for a fixed selection criterion $\mathcal{S}_{0}$, and hence $k(\delta),$
$\epsilon(\delta)=o(1/k^{2}(\delta))$. In other words, as the noise level
decreases the descent direction should get closer to the negative gradient
direction, i.e., when $\delta\rightarrow 0$,
$d_{k}^{\delta}\rightarrow-\tau\nabla_{x}\mathcal{D}(x_{k}^{\delta})\Leftrightarrow
z_{k}^{\delta}\rightarrow x_{k}^{\delta}$, and
$k_{i}(\delta)\rightarrow\infty$, for $i=1,2$, which results in
$x_{k_{1}(\delta)}^{\delta}\rightarrow\hat{x}\leftarrow
z_{k_{2}(\delta)}^{\delta}$.

The natural question that one can ask now is, when can a regularized
structured solution be obtained without encountering the semi-convergence in
the recovery errors? The answer is ‘never’ for those structure imposing
directions $d_{k}^{\delta}$ that satisfies (2.6) for all $k$, as then one
minimizes the LS functional (1.2) for noisy $b_{\delta}$, leading to the semi-
convergence of the recovery errors. However, one can relax the condition (2.6)
appropriately to circumvent the semi-convergence as well as recover a
structured regularized solution, but then, it won’t fall in the family
$\mathcal{I}_{3}$. Hence, one has to define an extension of the family
$\mathcal{I}_{3}$ to encompass such structured regularized solutions, details
in the later sections.

##  3\. PnP-algorithms as structured iterations

The core of any PnP-algorithm is the data-consistency step followed by the
data-denoising step, where the denoising operator doesn’t need to be a
proximal operator anymore. We first analyze the FBS-PnP algorithm and then the
ADMM-PnP algorithm.

###  3.1. Forward-backward splitting (FBS) - PnP

In this algorithm, for a fixed denoiser $H_{\sigma}$ (of denoising strength
corresponding to noise level $\sigma$) and starting from an initial choice
$z_{0}^{\delta}$, at any iteration step $k\geq 1$, we have

$\displaystyle z_{k-1}^{\delta}$
$\displaystyle\longmapsto\;x_{k}^{\delta}\;=\;z_{k-1}^{\delta}-\tau_{k}\nabla_{x}\mathcal{D}(z_{k-1}^{\delta})\;$
$\displaystyle\longleftarrow\mbox{ data-consistency step}$ (3.1)
$\displaystyle x_{k}^{\delta}$
$\displaystyle\longmapsto\;z_{k}^{\delta}\;=\;H_{\sigma_{k}}(x_{k}^{\delta}),\;$
$\displaystyle\longleftarrow\mbox{ data-denoising step},$ (3.2)

where $H_{\sigma_{k}}$ is the updated kth denoiser, with the denoising
strength corresponding to $\sigma_{k}=\tau_{k}\sigma$. Hence, an approach
similar to (2.12), provides a structure imposing direction $d_{k-1}^{\delta}$,
which can be defied as

$\displaystyle d_{k-1}^{\delta}$
$\displaystyle:=z_{k}^{\delta}-z_{k-1}^{\delta}$ (3.3)
$\displaystyle=\underbrace{-\tau\nabla_{x}\mathcal{D}(z_{k-1}^{\delta})}_{data-
consistency}\;+\;\underbrace{\left(H_{\sigma_{k}}(x_{k}^{\delta})-x_{k}^{\delta}\right)}_{data-
denoising}\;,$

and the direction $d_{k-1}^{\delta}$, as defined in (3.3), will be a descent
direction provided it satisfies (2.6). Therefore, for $H_{\sigma_{k}}$
satisfying

$||H_{\sigma_{k}}(x_{k}^{\delta})-x_{k}^{\delta}||\;<\;||-\tau\nabla_{x}\mathcal{D}(x_{k-1}^{\delta})||,$
(3.4)

we have $d_{k-1}^{\delta}$ satisfying (2.6), since

$\displaystyle\Big{(}d_{k-1}^{\delta},-\tau\nabla_{x}\mathcal{D}(z_{k-1}^{\delta})\Big{)}_{2}\geq\tau||\nabla_{x}\mathcal{D}(z_{k-1}^{\delta})||\left(\tau||\nabla_{x}\mathcal{D}(z_{k-1}^{\delta})||-||H_{\sigma_{k}}(x_{k}^{\delta})-x_{k}^{\delta}||\right)>0.$

In other words, (3.4) suggests that if the change (denoising) in the denoising
step is smaller than the change (improvement) in the data-consistency step,
than the resulting direction $d_{k-1}^{\delta}$ serves as a structure imposing
descent direction, as defined in §2.3, and hence, the iterates
$z_{k}^{\delta}$ recovers a structured regularized-solution
$z_{k(\delta)}^{\delta}$, depending on the selection criteria $\mathcal{S}$.
This makes sense, since, if the strength of the denoiser exceeds the changes
arising in the data-consistency step then it will over-denoise or over-smooth
the iterates, leading to an over-smoothed recovery. Of course, if
$||H_{\sigma_{k}}(x_{k}^{\delta})-x_{k}^{\delta}||$ is very small then, from
(3.3), $d_{k-1}^{\delta}$ will be very close to
$-\tau\nabla_{x}\mathcal{D}(z_{k-1}^{\delta})$ and thus, the structured
iterates $z_{k}^{\delta}$ will also be close to the unstructured (noisy)
iterates $x_{k}^{\delta}$, leading to an under-denoised or under-smoothed
recovery. Note that, the condition (3.4) is sufficient but not necessary for
the directions $d_{k}^{\delta}$ to be a descent direction, i.e.,
$d_{k}^{\delta}$ can satisfy (2.6) even when it does not satisfy (3.4).

Now, similar to (2.9), we would like to estimate a bound on the difference of
the new denoised iterates $z_{k}^{\delta}$ with the original LS-iterates
$y_{k}^{\delta}$ (generated from the negative gradients), i.e., starting from
the same initial point $z_{0}^{\delta}=y_{0}^{\delta}\equiv 0$ and with the
same step-size $\tau>0$, we have $x_{1}^{\delta}=y_{1}^{\delta}$ and

$\displaystyle||z_{1}^{\delta}-y_{1}^{\delta}||$
$\displaystyle=||H_{\sigma_{1}}(x_{1}^{\delta})-x_{1}^{\delta}||,$
$\displaystyle||z_{2}^{\delta}-y_{2}^{\delta}||$
$\displaystyle\leq||H_{\sigma_{2}}(x_{2}^{\delta})-x_{2}^{\delta}||+(1+\tau||A^{*}A||)||H_{\sigma_{1}}(x_{1}^{\delta})-x_{1}^{\delta}||$
$\displaystyle||z_{3}^{\delta}-y_{3}^{\delta}||$
$\displaystyle\leq||H_{\sigma_{3}}(x_{3}^{\delta})-x_{3}^{\delta}||+(1+\tau||A^{*}A||)||z_{2}^{\delta}-x_{2}^{\delta}||$
$\displaystyle\leq||H_{\sigma_{3}}(x_{3}^{\delta})-x_{3}^{\delta}||+(1+\tau||A^{*}A||)||H_{\sigma_{2}}(x_{2}^{\delta})-x_{2}^{\delta}||$
$\displaystyle\hskip
85.35826pt+(1+\tau||A^{*}A||)^{2}||H_{\sigma_{1}}(x_{1}^{\delta})-x_{1}^{\delta}||$
$\displaystyle\vdots$ $\displaystyle||z_{k}^{\delta}-y_{k}^{\delta}||$
$\displaystyle\leq\sum_{i=0}^{k-1}(1+\tau||A^{*}A||)^{i}||H_{\sigma_{k-i}}(x_{k-i}^{\delta})-x_{k-i}^{\delta}||.$
(3.5)

That is, for the $k(\delta)$ iteration the denoised iterate
$z_{k(\delta)}^{\delta}$ will be bounded away from the LS iterate
$y_{k(\delta)}^{\delta}$ at most by the factor given in (3.1), which also
depends on the denoising strength of the denoiser $H_{\sigma_{k}}$ at each
iteration $k$. Again, this does not guarantee that
$||z_{k(\delta)}^{\delta}-\hat{x}||\leq||y_{k(\delta)}^{\delta}-\hat{x}||$, as
it can happen that the iterates may be over-smoothed (for stronger denoiser
$H_{\sigma_{k}}$) and might have crossed the optimal solution. However, if we
trust (tune) our denoiser to be “appropriate” and directing the iterates
$x_{k}^{\delta}$ towards the true solution $\hat{x}$ through $z_{k}^{\delta}$,
i.e., towards a “better” direction, then we would hope the shift in (3.1) to
be in the proper direction, and hence, leading to
$||z_{k(\delta)}^{\delta}-\hat{x}||\leq||y_{k(\delta)}^{\delta}-\hat{x}||$,
i.e., a better recovery. In §3.1.1 and §3.1.2 we provide some suggestions to
tune a bad denoiser $H_{\sigma}$ appropriately, i.e., when the original
denoiser $H_{\sigma}$ is too strong or too weak, then how can it be improved
to perform the appropriate level of denoising to recover better estimates of
the true solution, without actually altering the original denoiser.

In this section we showed that, if the denoiser $H_{\sigma_{k}}$ satisfies the
descent condition (3.4) at each iteration k, then the FBS-PnP algorithm will
fall in the family $\mathcal{I}_{3}$. The weakness, which also the strength,
of the condition (3.4) is that, the recovery errors corresponding to the
iterates $z_{k}^{\delta}$ will also follow a semi-convergent path, and hence,
the iterations have to be terminated at an appropriate early instance
($k(\delta)$), as explained in remark 2.3. To overcome this inconvenience, as
stated in remark 2.3, the family $\mathcal{I}_{3}$ needs to be further
extended, which is done later.

####  3.1.1. Reinforcing weaker denoisers

In practice, one usually starts with a denoiser of certain strength, say
$\sigma$, and at every iteration $k$, the denoising of $x_{k}^{\delta}$ is
done via $H_{\sigma_{k}}$, where $\sigma_{k}=\sigma\tau_{k}$, where $\tau_{k}$
is the associated step-size for that iteration. As discussed in §3.1 the FBS-
PnP algorithm is guaranteed to fall in the regularized family
$\mathcal{I}_{3}$ provided the denoiser $H_{\sigma_{k}}$ satisfy (3.4). Now,
for a weaker denoiser (i.e.,
$||H_{\sigma_{k}}(x_{k}^{\delta})-x_{k}^{\delta}||$ is small), though the
structure imposing direction $d_{k}^{\delta}$, as defined in (3.3), will
satisfy (3.4), the iterates $x_{k}^{\delta}$ won’t be sufficiently denoised,
i.e., the denoised iterates $z_{k}^{\delta}$ will still contain certain amount
of noise. Of course, during the next minimization process, one can opt for an
“appropriate” stronger denoiser $H_{\sigma_{k}}$, if possible, to provide
stronger denoising of the iterates $x_{k}^{\delta}$. There is also a second
approach to augment the denoising strength of $d_{k}^{\delta}$, without
altering the denoiser, for example, when there is only one denoiser to work
with. Observe that, the satisfaction of the condition (3.4), by a denoiser
$H_{\sigma_{k}}$, is also dependent on the step-size $\tau$. That is, even for
a weaker denoiser $H_{\sigma_{k}}$, if the value of $\tau$ is small enough,
then $H_{\sigma_{k}}$ will fail to satisfy the condition (3.4), i.e., the
denoising strength of the weaker denoiser will increase implicitly. Hence, by
decreasing the step-size $\tau$ smaller than the previous choice, such that
$||H_{\sigma_{k}}(x_{k}^{\delta})-x_{k}^{\delta}||$ is only slightly smaller
than $||-\tau\nabla_{x}\mathcal{D}(z_{k-1}^{\delta})||$ will, not only lead to
a descent direction but also, the amount of denoising will be increased for
the original denoiser $H_{\sigma}$, see Example 5.1. However, this does
increase the computational time, as with smaller values of $\tau$ the descent
rate also decreases.

####  3.1.2. Attenuating stronger denoisers

In contrast to a weak denoiser, a stronger denoisers $H_{\sigma}$ may not lead
to the satisfaction of the condition (3.4), since, the extent of denoising
$||H_{\sigma_{k}}(x_{k}^{\delta})-x_{k}^{\delta}||$ may exceed the
improvements coming from the data-consistency step
$||x_{k}^{\delta}-z_{k-1}^{\delta}||$. Again, the obvious way to overcome this
issue is by opting for a weaker denoiser $H_{\sigma}$, the next time, such
that the condition (3.4) is satisfied. However, this can lead to several
attempts of numerous minimization process to find an appropriate denoiser,
which can substantially increase the computational time. Note that, unlike the
previous tweak, here one cannot increase the step-size to compensate the
strong denoising in (3.4), as then the descent process will break down, for
larger step-sizes. A cleverer approach of handling this issue is through
relaxing the denoising strength of the strong denoiser $H_{\sigma}$ by
incorporating an additional relaxation parameter $\alpha\in(0,1]$, via

$\displaystyle z_{k}^{\delta}=z_{k}^{\delta}(\alpha)$
$\displaystyle:=\;x_{k}^{\delta}\;+\;\alpha\left(H_{\sigma_{k}}(x_{k}^{\delta})\;-\;x_{k}^{\delta}\right)$
(3.6)
$\displaystyle=(1-\alpha)\;x_{k}^{\delta}\;+\;\alpha\;H_{\sigma_{k}}(x_{k}^{\delta}).$

Hence, if
$||H_{\sigma_{k}}(x_{k}^{\delta})-x_{k}^{\delta}||>||x_{k}^{\delta}-z_{k-1}^{\delta}||=||-\tau\nabla_{x}\mathcal{D}(z_{k-1}^{\delta})||$,
then for any

$\alpha\;\leq\;\frac{||x_{k}^{\delta}-z_{k-1}^{\delta}||}{||z_{k}^{\delta}-x_{k}^{\delta}||}\;=\;\frac{||-\tau\nabla_{x}\mathcal{D}(z_{k-1}^{\delta})||}{||H_{\sigma_{k}}(x_{k}^{\delta})-x_{k}^{\delta}||}\;<1,$
(3.7)

we will have the resulting direction

$\displaystyle d_{k}^{\delta}$ $\displaystyle=z_{k}^{\delta}-z_{k-1}^{\delta}$
$\displaystyle=(x_{k}^{\delta}-z_{k-1}^{\delta})\;+\;\alpha\left(H_{\sigma_{k}}(x_{k}^{\delta})-x_{k}^{\delta}\right)$
(3.8)

satisfying the condition (2.6). Therefore, if $\alpha<<1$ (small) then
$d_{k}^{\delta}$ is close to $-\tau\nabla_{x}\mathcal{D}(z_{k-1}^{\delta})$,
which leads to weak denoising, and if $\alpha\approx 1$ (close to 1) then it
leads to strong denoising. This can be seen in Table 2, for Example 5.2.
Usually, one can choose
$\alpha=\gamma\;\frac{||-\tau\nabla_{x}\mathcal{D}(z_{k-1}^{\delta})||}{||H_{\sigma_{k}}(x_{k}^{\delta})-x_{k}^{\delta}||}$,
where the proximity parameter $\gamma\leq 1$ controls the closeness of
$d_{k}^{\delta}(\alpha(\gamma))$ to
$-\tau\nabla_{x}\mathcal{D}(z_{k-1}^{\delta})$ through the parameter
$\alpha(\gamma)$.

An alternative way of estimating a proper $\alpha$ value is through minimizing
$z_{k}^{\delta}(\alpha)$ based on the selection criterion, i.e.,

$\displaystyle\alpha_{0}:=$
$\displaystyle\operatorname*{arg\,min}_{\alpha\in(0,1]}\;\;\;\mathcal{S}(z_{k}^{\delta}(\alpha))$
(3.9) $\displaystyle\mbox{such
that,}\;\;z_{k}^{\delta}(\alpha)=x_{k}^{\delta}+\alpha\left(H_{\sigma_{k}}(x_{k}^{\delta})\;-\;x_{k}^{\delta}\right).$

Note that, the minimization problem (3.9) may not be a strictly convex, i.e.,
there might not be a global minimizer $\alpha_{0}$. Nevertheless, this is
simply a sub-problem intended to find an appropriate $\alpha$ value between 0
and 1, depending on the selection criterion $\mathcal{S}$, and hence, even if
the best $\alpha_{0}$ is not obtained (which is not a requirement), any
$\alpha\in(0,1)$ will reduce the denoising strength and empirical experiments
(see Example 5.2) show that it works fine. Also, the alpha value from (3.9)
doesn’t guarantee that
$d_{k}^{\delta}(\alpha_{0}):=z_{k}^{\delta}(\alpha_{0})-z_{k-1}^{\delta}$ will
satisfy (2.6), but, it does attenuates the strong denoising of
$H_{\sigma_{k}}$ and provides an iterate $z_{k}^{\delta}$ which satisfies the
selection criterion the best.

######  Remark 3.1.

Note that, inspired from (3.9), one can also choose the denoising strength of
the denoiser $H_{\sigma_{k}}$, at every iteration, depending on the selection
criterion, i.e., by defining the denoised iterates as a function of the
denoising strength $\sigma$

$z_{k}^{\delta}(\sigma_{k}):=H_{\sigma_{k}}(x_{k}^{\delta}),$ (3.10)

one can choose the denoising strength $\sigma_{k}$, within a certain range,
which best satisfies the selection criterion, i.e.,

$\displaystyle\sigma_{k}:=$
$\displaystyle\operatorname*{arg\,min}_{\sigma\in[\sigma_{{min}},\;\sigma_{max}]}\;\;\;\mathcal{S}(z_{k}^{\delta}(\sigma))$
(3.11) $\displaystyle\mbox{such that,
}\;\;z_{k}^{\delta}(\sigma)=H_{\sigma}(x_{k}^{\delta}).$

However, the denoising step can be computationally expensive step during the
iterative process, and hence, performing a minimization over the denoising
strength can notably increase the computational time. In addition, similar to
(3.9), the minimization in (3.11) need not be even convex.

###  3.2. Fast Forward-backward splitting (Fast FBS) - PnP

The plain FBS-PnP method tends to show very slow convergence rate, even with
the semi-iterative approach, i.e., $k(\delta)$ is very large. To speed up the
descent rate, faster algorithms are usually employed such as, the fast
proximal gradient methods (FPGM) or FISTA, see [11]. The iterates in this
algorithm have an intermediate momentum step to accelerate the descent
process, which is given by, starting from initial
$z_{0}^{\delta}=x_{0}^{\delta}$ and $t_{0}=1$, for $k\geq 1$,

$\displaystyle z_{k-1}^{\delta}$
$\displaystyle\longmapsto\;x_{k}^{\delta}\;=\;z_{k-1}^{\delta}-\tau_{k}\nabla_{x}\mathcal{D}(z_{k-1}^{\delta})\;$
$\displaystyle\longleftarrow\mbox{ data-consistency step}$ (3.12)
$\displaystyle x_{k}^{\delta}$
$\displaystyle\longmapsto\;z_{k}^{\delta}\;=\;H_{\sigma_{k}}(x_{k}^{\delta}),\;$
$\displaystyle\longleftarrow\mbox{ data-denoising step}$ (3.13) $\displaystyle
z_{k}^{\delta}$
$\displaystyle\longmapsto\;z_{k}^{\delta}\;=\;z_{k}^{\delta}\;+\;\alpha_{k}(z_{k}^{\delta}-z_{k-1}^{\delta}),\;$
$\displaystyle\longleftarrow\mbox{ momentum-step}$ (3.14)

where $\alpha_{k}=\frac{t_{k-1}-1}{t_{k}}$ with
$t_{k}=\frac{(1+\sqrt{1+4t_{k-1}^{2}})}{2}$.

Now, one can also perform similar analysis as done for the vanilla FBS-PnP
method, we skip the details.

##  4\. Expansion of the regularization family

Before we proceed to attach a meaning to the ADMM-PnP algorithms we would like
to address the question raised in the last part of remark 2.3, i.e., can the
family of regularized solutions $\mathcal{I}_{3}$ be generalized to soften the
condition (2.6), which leads to the semi-convergence of the recovery errors?
To answer this question let’s first interpret the family $\mathcal{I}_{3}$
from a different perspective. Note that, the formulation of $\mathcal{I}_{3}$
in (2.11) is also equivalent to

$\displaystyle\mathcal{I}_{3}^{\prime}:=\left\\{x^{\delta}(\mathcal{D},k,d_{k}^{\delta};\mathcal{S}):\;x_{k}^{\delta}\;=\;x_{k-1}^{\delta}\;+\;d_{k-1}^{\delta},\mbox{
subject to }\mathcal{D}(x_{k}^{\delta})\geq\epsilon(\delta,\mathcal{S})\geq
0,\right.$ $\displaystyle\left.\mbox{such that $d_{k}^{\delta}$ satisfies
\eqref{LS descend prop.} and $\mathcal{S}$ is a selection
criterion.}\right\\},$ (4.1)

where the lower bound $\epsilon(\delta,\mathcal{S})$ for the discrepancy term
$\mathcal{D}$ now serves the purpose of $k(\delta)$ and determines the extent
of regularization, for example, if $\mathcal{S}$ is the Morozov’s discrepancy
principle (DP), then $\epsilon(\delta)\approx\eta\delta$, for $\eta>1$, see
[5, 33, 25]. Equivalently, the formulation in $\mathcal{I}_{3}^{\prime}$
suggests to minimize the data-discrepancy functional $\mathcal{D}$, via a
(descent) direction $d_{k}^{\delta}$ satisfying (2.6), but constraint to
$\mathcal{D}(x_{k}^{\delta})\geq\epsilon(\delta,\mathcal{S})$, i.e., not to
minimize $\mathcal{D}$ completely, which avoids overfitting to the noisy data
$b_{\delta}$. Therefore, the key is to not overfit the noisy data, which
motivates us to generalize the constraint
$\mathcal{D}(x_{k}^{\delta})\geq\epsilon(\delta)$, for larger iterations
($k>k(\delta)$), to

$0\leq\epsilon_{1}(\delta)\leq\mathcal{D}(x_{k}^{\delta})\leq\epsilon_{2}(\delta),$
(4.2)

such that $\epsilon_{2}(\delta)\xrightarrow{\delta\rightarrow 0}0$ (so that
$x_{k}^{\delta}\xrightarrow{\delta\rightarrow 0}x^{\dagger}$), where the lower
bound $\epsilon_{1}(\delta)$ avoids overfitting of the noisy data and the
upper bound $\epsilon_{2}(\delta)$ avoids underfitting of the data. Now, if we
enforce the directions $d_{k}^{\delta}$ to satisfy (2.6) for $x_{k}^{\delta}$
satisfying $\mathcal{D}(x_{k}^{\delta})>\epsilon_{2}(\delta)$, we minimize the
discrepancy-term $\mathcal{D}$, where as, removing the constraint (2.6) on
those directions $d_{k}^{\delta}$ for which
$\mathcal{D}(x_{k}^{\delta})\leq\epsilon_{2}(\delta)$, but constraining them
implicitly through $x_{k}^{\delta}$ such that
$\epsilon_{1}(\delta)\leq\mathcal{D}(x_{k}^{\delta}+d_{k}^{\delta})\leq\epsilon_{2}(\delta)$,
one is able to relax the descent criterion while (still) avoiding the
overfitting to the noisy data. Hence, the new conditions for the directions
$d_{k}^{\delta}$ are

$\displaystyle\Big{(}d_{k}^{\delta},-\nabla_{x}\mathcal{D}(x_{k}^{\delta})\Big{)}_{2}>0,\;\;\;\mbox{
for }\;\mathcal{D}(x_{k}^{\delta})>\epsilon_{2}(\delta)\;\mbox{ or }\;k\leq
k(\delta),\mbox{ and}$ (4.3)
$\displaystyle\mathcal{D}(x_{k}^{\delta}+d_{k}^{\delta})\geq\epsilon_{1}(\delta),\;\;\;\
\mbox{ for }\;\mathcal{D}(x_{k}^{\delta})\leq\epsilon_{2}(\delta)\;\mbox{ or
}\;k>k(\delta).$

Note that, with this modification, the recovery errors for the iterates
$x_{k}^{\delta}$ does not necessarily follow a semi-convergence trail, nor a
convergence, anymore; it may end up oscillating/fluctuating around the optimal
solution, see Example 5.2. Again, the answer to “does it improve the recovery
process?”, depends on the choice of the recovery directions $d_{k}^{\delta}$
and the selection criteria $\mathcal{S}$. The formal extension of the family
$\mathcal{I}_{3}^{\prime}$ can be expressed as

$\displaystyle\mathcal{I}_{4}:=$
$\displaystyle\left\\{x^{\delta}(\mathcal{D},k,d_{k}^{\delta};\mathcal{S})\right.\;:\;\left.x_{k}^{\delta}=x_{k-1}^{\delta}\;+\;d_{k-1}^{\delta},\;\mbox{
s.t. $d_{k-1}^{\delta}$ satisfies \eqref{descent dir. gen.},}\right.$ (4.4)
$\displaystyle\left.\;\mbox{ and $\mathcal{S}$ is a selection
criterion.}\right\\}.$

Observe that, by increasing the number of constraints from one, in equation
(4), to two, in equation (4.4), we are able to expand the family of
regularized solutions, i.e., for $\epsilon_{1}(\delta)<\epsilon_{2}(\delta)$,
we have $\mathcal{I}_{4}\supset\mathcal{I}_{3}^{\prime}$ and when
$\epsilon_{1}(\delta)=\epsilon_{2}(\delta)$, then
$\mathcal{I}_{4}=\mathcal{I}_{3}^{\prime}$. Moreover, the family
$\mathcal{I}_{4}\supseteq\mathcal{T}_{2}$ (the family of Tikhonov-
regularization), since, the discrepancy of a regularized solution in
$\mathcal{T}_{2}$ is always positive, for any $\lambda(\delta)>0$, i.e.,
$\mathcal{D}(x^{\delta}(\mathcal{D},\lambda(\delta),\mathcal{R}))>0$, and the
directions $d_{k}^{\delta}$ in $\mathcal{I}_{4}$ need not satisfy (2.6) for
all $k$. Moreover, the advantage of defining a regularized solution as in
$\mathcal{I}_{4}$ over $\mathcal{T}_{2}$ is that, one does not have to worry
about the convergence of $\\{x_{k}^{\delta}\\}$ anymore, since we are not
interested in finding the minimizer of a Variational problem (1.3), rather, we
would like to estimate the solution of the inverse problem (1.1) iteratively
via a sequence of well-defined iterates $x_{k}^{\delta}$, depending on the
directions $d_{k}^{\delta}$, where the best estimate
$x_{k(\delta,\mathcal{S})}^{\delta}$ to the true solution $\hat{x}$ is
determined via a selection criterion $\mathcal{S}$.

Observe that, although the initial iterate $x_{0}^{\delta}$ is not explicitly
mentioned in the families of the regularized solutions
$\mathcal{I}_{1},\mathcal{I}_{2},\mathcal{I}_{3}$ and $\mathcal{I}_{4}$, it is
inherently embedded in the recovery process, unlike in the families
$\mathcal{T}_{1}$ and $\mathcal{T}_{2}$, for convex regularizers
$\mathcal{R}$, where $x_{0}^{\delta}$ does not matter. To explicitly state the
importance of the initial iterates in these (non-Variational or (semi-)
iterative) families, $\mathcal{I}_{4}$ can be expanded to

$\displaystyle\mathcal{I}_{5}:=$
$\displaystyle\left\\{x^{\delta}(\mathcal{D},k,d_{k}^{\delta},x_{0}^{\delta};\mathcal{S})\right.:\left.x_{k}^{\delta}=x_{k-1}^{\delta}\;+\;d_{k-1}^{\delta},\;\mbox{
starting from $x_{0}^{\delta}$},\right.$ $\displaystyle\left.\;\mbox{s.t.
$d_{k-1}^{\delta}$ satisfies \eqref{descent dir. gen.} and $\mathcal{S}$ is a
selection criterion.}\right\\}.$ (4.5)

The significance of the initial iterate $x_{0}^{\delta}$ can be prominently
seen in the ADMM-PnP algorithm, which is discussed later. However, if the
descent directions depends continuously, as a function, on its iterates, then
for a small change in the initial iterates the difference in the recovered
regularized solutions will be bounded by a factor similar to (2.9), i.e., if
$||d_{k}^{\delta}(x_{k}^{\delta})-d_{k}^{\delta}(z_{k}^{\delta}))||$ is small
for small $||x_{k}^{\delta}-z_{k}^{\delta}||$ and
$||x_{0}^{\delta}-z_{0}^{\delta}||\leq\epsilon$, then
$||z_{k}^{\delta}-x_{k}^{\delta}||\leq C(k)\epsilon$. This implies stability
in the recovery process with respect to any small perturbations in the initial
choices.

###  4.1. Differential equation’s solutions as a regularization family

As one might have guessed by now, all the above families of solutions can be
categorized as the solution states corresponding to a discrete approximation
of some differential equations, in particular, the evolution equation given by

$\displaystyle\frac{\partial}{\partial t}x^{\delta}(.,t)$
$\displaystyle=f(x^{\delta}(.,t))$ (4.6) $\displaystyle x^{\delta}(.,t_{0})$
$\displaystyle=x_{0}^{\delta},$

where the function $f(x^{\delta}(.,t))$ determines the flow or the evolution
of a family of estimated solutions $x^{\delta}(.,t)$ to the inverse problem
(1.1), and an appropriate estimate $x^{\delta}(.,T(\delta))$ for $\hat{x}$ is
determined via a selection criterion. When
$f(x^{\delta}(.,t))=-\tau(t)\nabla_{t}\mathcal{D}(x^{\delta}(.,t))$, for
$\tau(t)\geq 0$, we have the classical steepest-descent method in the
continuous settings. For $f(x^{\delta}(.,t))$ satisfying
$\Big{(}f(x^{\delta}(.,t)),-\nabla_{t}\mathcal{D}(x^{\delta}(.,t))\Big{)}_{2}>0$,
it falls in the family $\mathcal{I}_{3}$, where as, for $f(x^{\delta}(.,t))$
satisfying (4.3), we end in the family $\mathcal{I}_{5}$, in their respective
continuous versions. Here, one can observe that a regularized solution is not
only a function of $f(x^{\delta}(.,t))$ (the flow direction), but also, on
$x_{0}^{\delta}$ (the initial/starting state) and $T(\delta)$ (the
stopping/terminal time), i.e.,
$x^{\delta}=x^{\delta}(.;f,x_{0}^{\delta},T(\delta))$. In contrast, the
minimization of the (convex) Variational problem (1.3) is independent of the
starting point and the descent direction, when minimized in an iterative
manner, since a good starting point ($x_{0}^{\delta}$) and a faster descent
direction $(f(x^{\delta}(.,t)))$ only improves the descent rate to reach the
minimizer of the (1.3), but not the final solution. It only depends on the
definition of the penalized/constrained functional ($\mathcal{R}$) and the
penalty/regularization parameter ($\lambda$) in (1.3), i.e.,
$x^{\delta}(.;\mathcal{R},\lambda)$. Therefore, based on a fixed selection
criterion $\mathcal{S}$, comparing a Variational regularized solution
$x^{\delta}(.;\mathcal{R},\lambda)$ to a (semi-) iterative regularized
solution $x^{\delta}(.;f,T,x_{0}^{\delta})$ we see that the analogous of
$\mathcal{R}$ is $f$ and $\lambda$ is $T$. However, the (true) dependence of
the initial point $x_{0}^{\delta}$ is still missing in the Variational
solution (assuming one minimizes the convex functional completely and the
iterates converge to the global minimizer) and has an importance in the
(semi-) iterative method. Nevertheless, the recovery will be stable to a small
perturbation in the initial condition, if $f(x^{\delta}(.,t))$ depends
continuously on $x^{\delta}(.,t)$, i.e., one can estimate a bound for the
difference of the solution flow at any time $t$, which is given by

$||x^{\delta}(.,t)-z^{\delta}(.,t)||\leq
c(t)\int_{0}^{t}||f(x^{\delta}(.,s))-f(z^{\delta}(.,s))||ds\;+\;||x_{0}^{\delta}-z_{0}^{\delta}||,$
(4.7)

for some constant $c(t)$, which implies

$||x^{\delta}(.,t)-z^{\delta}(.,t)||\leq
c(t)\int_{0}^{t}\epsilon^{\delta}(s)ds+||x_{0}^{\delta}-z_{0}^{\delta}||\;\leq\;C(t)\epsilon,$
(4.8)

if $||f(x^{\delta}(.,s))-f(z^{\delta}(.,s))||\leq\epsilon^{\delta}(t)$, for
$t_{0}\leq s\leq t$, and $||x_{0}^{\delta}-z_{0}^{\delta}||\leq\epsilon$.

Note that, if $f(x^{\delta}(.,t))$ satisfies (2.6) or (4.3) then it is
implicitly dependent on the discrepancy term $\mathcal{D}(x)$ through the
inner product with its gradient $-\nabla_{x}\mathcal{D}(x_{k}^{\delta})$,
i.e., $f(x^{\delta}(.,t))=f(x^{\delta}(.,t);\mathcal{D})$. Hence, to further
generalize the class $\mathcal{I}_{5}$, one can even remove any dependence
(implicit or explicit) on the discrepancy term $\mathcal{D}$, i.e., no
knowledge of the forward operator $A$ is used. For example, one can generate
the flow of solutions $x^{\delta}(.,t)$ driven solely by a collection of
informative dataset, i.e., using machine learning or deep learning algorithms,
for example, the AUTOMAP network in [34]. The advantage of this formulation is
that, any short-comings or mismatch in the model $A$ can be compensated by
replacing it with a data based model (say neural networks). However, this
strength also becomes its weakness for certain (adversarial) noisy data and
pushes the recovery process further into the realm of instabilities (now being
generated by the data model), see [35], if no “proper regularization” is
employed, where the proper regularization implies regularizing the recovery
process rather than regularizing the neural network, this is a topic of
discussion for another paper. Nevertheless, we end this section by defining
the largest family of regularized solution to the inverse problem (1.1) as

$\displaystyle\mathcal{I}_{6}:=$
$\displaystyle\left\\{x^{\delta}(f,x_{0}^{\delta},t;\mathcal{S})\right.\;:\;x^{\delta}(f,x_{0}^{\delta},t)\;\mbox{
is a solution state of \eqref{evolution equation} at time t,}$ (4.9)
$\displaystyle\left.\mbox{s.t.
$f(x^{\delta}(.,t))\xrightarrow{\delta\rightarrow
0}-\tau(t)\nabla_{x}\mathcal{D}(x(.,t))$ and $\mathcal{S}$ is a selection
criterion.}\right\\}.$

###  4.2. Alternating direction method of multipliers (ADMM) - PnP

Now we have developed all the theories needed to attach a meaning to the ADMM-
PnP algorithms. Note that, a typical step in the ADMM-PnP algorithm is given
by

$\displaystyle x_{k+1}^{\delta}$
$\displaystyle=\text{Prox}_{\frac{1}{\rho}\mathcal{D}}(z_{k}^{\delta}-u_{k}^{\delta})\;\;\longleftarrow\text{
data-consistency step}$ (4.10)
$\displaystyle=\operatorname*{arg\,min}_{x}\;\;\mathcal{D}(x)+\rho||x-(z_{k}^{\delta}-u_{k}^{\delta})||_{2}^{2}$
$\displaystyle z_{k+1}^{\delta}$
$\displaystyle=H_{\sigma_{k}}(x_{k+1}^{\delta}+u_{k}^{\delta})\;\;\longleftarrow\text{
data-denoising step}$ (4.11) $\displaystyle u_{k+1}^{\delta}$
$\displaystyle=u_{k}^{\delta}+x_{k+1}^{\delta}-z_{k+1}^{\delta}\;\;\longleftarrow\text{
noise update step}$ (4.12)

where if $H_{\sigma_{k}}$, for $\sigma_{k}=\sigma/\rho$, corresponds to the
proximal map of a proper, closed and convex regularizer $\mathcal{R}(x)$, then
the Lagrangian parameter ($\rho>0$) only effects the speed of the convergence
and not the solution of (1.3), i.e., the iterates $x_{k}^{\delta}$ and
$z_{k}^{\delta}$ will converge to the minimzer of the functional (1.3),
irrespective of the $\rho$-value. However, this can not be guaranteed for any
general denoiser $H_{\sigma_{k}}$, in fact, there might not even exists an
associated proximal map corresponding to any general denoiser, and hence, a
regularization function $\mathcal{R}(x)$ in (1.3). In this section, we connect
ADMM-PnP algorithms to the theories developed in the previous sections and
also point out the significance of the Lagrangian parameter $\rho$ in this
setting. Note that, at any step $k$, the minimizer ($x_{k+1}^{\delta}$) in the
data-consistency step (4.10) has a closed form solution, for
$\mathcal{D}(x)=||Ax-b_{\delta}||_{2}^{2}$, given by

$\displaystyle x_{k}^{\delta}=(A^{*}A+\rho
I)^{-1}(A^{*}b_{\delta}+\rho(z_{k-1}^{\delta}-u_{k-1}^{\delta}))\;\;\longleftarrow\;\mbox{Tikhonov-$\ell_{2}$
solution with $\lambda=\rho$}.$ (4.13)

Therefore, the resulting direction from $x_{k}^{\delta}$ to
$x_{k+1}^{\delta}$, for $k\geq 1$, is given by

$d_{k}^{\delta}:=x_{k+1}^{\delta}-x_{k}^{\delta}=(A^{*}A+\rho
I)^{-1}\rho\left[(z_{k}^{\delta}-z_{k-1}^{\delta})-(u_{k}^{\delta}-u_{k-1}^{\delta})\right],$
(4.14)

and it’s a descent direction if it satisfy (2.6), i.e.,

$\Big{(}d_{k}^{\delta},-\nabla_{x}\mathcal{D}(x_{k})\Big{)}_{2}=\Big{(}(A^{*}A+\rho
I)^{-1}\rho\left[(z_{k}^{\delta}-z_{k-1}^{\delta})-(u_{k}^{\delta}-u_{k-1}^{\delta})\right],-A^{*}(Ax_{k}^{\delta}-b_{\delta})\Big{)}_{2}>0$
(4.15)

which, though is trivial for $k=0$ with
$x_{0}^{\delta}=z_{0}^{\delta}=u_{0}^{\delta}\equiv 0$, as

$\Big{(}d_{0}^{\delta},-\nabla_{x}\mathcal{D}(x_{0})\Big{)}_{2}=\Big{(}(A^{*}A+\rho
I)^{-1}A^{*}b_{\delta},A^{*}b_{\delta}\Big{)}_{2}>0,$ (4.16)

is non-trivial for $k\geq 1$. However, this is where we can either use the
regularization family $\mathcal{I}_{5}$, as $d_{0}^{\delta}$ satisfies (2.6),
or $\mathcal{I}_{6}$, where $d_{k}^{\delta}$ doesn’t need to satisfy any
conditions, to attach a meaning to the iterative process of generating a
regularized family of structured-solutions $\\{x_{k}^{\delta}\\}$, where an
appropriate estimate ($x_{k(\delta)}^{\delta}$) to the true solution
($\hat{x}$) solution is determined by some selection criteria. Note that, for
$\epsilon_{1}(\delta)\leq\epsilon_{2}(\delta)$, there is a $\rho(\delta)$ such
that
$\epsilon_{1}(\delta)\leq\mathcal{D}(x_{1}^{\delta}(\rho(\delta)))\leq\epsilon_{2}(\delta)$,
where
$x_{1}^{\delta}(\rho(\delta))=(A^{*}A+\rho(\delta)I)^{-1}A^{*}b_{\delta}$, and
for $\delta\rightarrow 0$, we need $\rho(\delta)\rightarrow 0$, since
$\epsilon_{2}(\delta)\xrightarrow{\delta\rightarrow 0}0$.

As mentioned earlier, the ADDM-minimization corresponding to a traditional
regularizer $\mathcal{R}(x)$ is independent of the Lagrangian parameter
$\rho$, as the iterates converge. But, for any general denoiser $H_{\sigma}$
there is no proof of convergence and hence, the choice of $\rho$ will affect
the recovered regularized solution. In addition, the initial choices
$(x_{0}^{\delta},z_{0}^{\delta},u_{0}^{\delta})$ also influence the recovery
process in this case, if we consider the iterative process to be in the family
$\mathcal{I}_{5}$ or $\mathcal{I}_{6}$ of regularized solutions. To show the
nature of dependence on $\rho$, we first transform the minimization in (4.10)
to an equivalent constraint minimization problem, given by

$\displaystyle
x_{k+1}^{\delta}:=\begin{cases}\operatorname*{arg\,min}_{x}\;\;&\mathcal{D}(x)=||Ax-
b_{\delta}||_{2}^{2},\\\ \mbox{subject to
}\;&||x-(z_{k}^{\delta}-u_{k}^{\delta})||_{2}^{2}\;\leq\;\epsilon(\rho),\end{cases}$
(4.17)

for some associated $\epsilon(\rho)\sim 1/\rho$. That is, for larger values of
$\rho$, $x_{k+1}^{\delta}$ is closer to $z_{k}^{\delta}-u_{k}^{\delta}$, and
for smaller values of $\rho$, $x_{k+1}^{\delta}$ will be farther away from
$z_{k}^{\delta}-u_{k}^{\delta}$, and closer to the noise corrupted solution
$(A^{*}A+\rho I)^{-1}(A^{*}b_{\delta}+\rho(z_{k}^{\delta}-u_{k}^{\delta}))$.
Hence, we divide the analysis into two parts:

####  4.2.1. When Lagrangian parameter is large

Assuming we start from $x_{0}^{\delta}=z_{0}^{\delta}=u_{0}^{\delta}\equiv 0$,
after the $1^{st}$ step, we have a slightly improved $x_{1}^{\delta}$ towards
the LS-solution (without being much corrupted by the noise), and hence, even
for a weaker denoiser $H_{\sigma_{1}}$, one can denoise it to get
$z_{1}^{\delta}$. On the $2^{nd}$ step, $x_{2}^{\delta}$ will again be
improved only slightly towards the noisy LS-solution, because of the stronger
constraint in (4.17), which can then be effectively denoised by a weaker
denoiser and so on. Therefore, starting from smooth initial choices, one can
recover well denoised iterates $\\{z_{k}^{\delta}\\}$, even for weaker
denoisers, see Example 5.3. This process is very similar to FBS-PnP algorithm,
though the iterates here are not explicitly generated by any descent
directions, rather, through a constrained minimization at each step.

####  4.2.2. When Lagrangian parameter is small

In contrast to the above recovery process, here, even when the initial choices
$(x_{k}^{\delta},z_{k}^{\delta},u_{k}^{\delta})$ are smooth (say 0s), after
the $1^{st}$ iteration we will have $x_{1}^{\delta}$ close to the noisy LS-
solution ($x^{\dagger}_{\delta}$), i.e., $x_{1}^{\delta}$ too noisy, and
hence, we need a much stronger denoiser to clean the corruptions; which in
turn can move $z_{1}^{\delta}$ far away from $x_{1}^{\delta}$, because of the
weaker constraint in (4.17). Thus, the denoisers $H_{\sigma_{k}}$ need to be
really effective in cleaning the heavily noised iterates $x_{k}^{\delta}$ at
every step, to produce well-denoised iterates $z_{k}^{\delta}$, only then we
may hope to recover a clean regularized solution. However, for weaker to
moderately strong denoisers, unlike the FBS-PnP algorithm, here one might
recover a noisy estimate to the solution of (1.1), see Example 5.3.

######  Remark 4.1.

From the above two analysis, one can notice that when $\rho$ is large the
recovery process is very similar to the FBS-PnP method and one can clean the
noise in the iterates $x_{k}^{\delta}$ gradually towards the “improved”
iterates $z_{k}^{\delta}$, of course, the improvement depends on the denoiser
$H_{\sigma_{k}}$. But, the rate of descent can be very slow, similar to ISTA,
and hence, it will be preferable to opt for the Fast FBS-PnP methods. Where
as, if $\rho$ is small then the rate of descent corresponding to the
discrepancy term $\mathcal{D}(x)$ is fast, however, this doesn’t ensure that
the rate of recovery for a denoised solution to be fast too, since, the faster
$\mathcal{D}$ decreases, the faster the noise increases in the iterates.
Therefore, here, it greatly depends on the effectiveness of the denoisers to
obtain a cleaner recovery.

######  Remark 4.2.

In addition to the issues mentioned in the above remark, note that, in the
ADMM minimization process, step (4.10) involves an additional minimization
problem at each iteration $k$. This further increases the computational time
and complexity, for example, for large scale problem the matrix inversion
$(A^{*}A+\rho I)^{-1}$ is impossible and one has to approximate the exact
minimizer $x_{k}^{\delta}$ via an iterative scheme (such as conjugate-
gradient) for a certain number of fixed iterations, and hence, there is a
further degradation to the proper definition of ADMM-minimization. Therefore,
unless one has some prior knowledge of an “appropriate” Lagrangian parameter
$\rho$ and “appropriate” denoisers $H_{\sigma_{k}}$, which can balance the
speed of $\mathcal{D}$-descent and the cleaning of the noisy iterates, we
prefer the fast FBS-PnP algorithms, where even inappropriate denoisers can be
tweaked appropriately to produce better results, as it involves a gradual
recovery process. However, for appropriately strong $H_{\sigma_{k}}$ and
$\rho$-value, one can get better performance (empirically) in ADMM-PnP than in
FBS-PnP. Nevertheless, considering both these algorithms as iterative
regularization processes in either family $\mathcal{I}_{5}$ or
$\mathcal{I}_{6}$, the better one is the one that satisfies the selection
criterion the best.

######  Remark 4.3.

Note that, the first iterate of the ADMM process (starting from initial
$x_{0}^{\delta}=z_{0}^{\delta}=u_{0}^{\delta}\equiv 0$) is an
$\ell_{2}$-Tikhonov solution, i.e.,
$x_{1}^{\delta}=\operatorname*{arg\,min}_{x}\;\mathcal{D}(x)+\rho||x||_{2}^{2}$.
Hence, one can generalize it by replacing the standard $\ell_{2}$-Tikhonov
solution with any general $\ell_{2}$-Tikhonov solution, i.e.,

$x_{1}^{\delta}\;=\;\operatorname*{arg\,min}_{x}\;\mathcal{D}(x)\;+\;\rho||Lx||_{2}^{2},$
(4.18)

where $L$ can be either any wavelet transformation or the gradient operator.
This may serve as a better initial $x_{1}^{\delta}$ for the ADMM-PnP process,
than the standard $\ell_{2}$-solution and can lead to an improved regularized
solution.

##  5\. Numerical Results

In this section we provide computational results validating the theories
developed in this paper. Note that, the goal here is not to show the
effectiveness of PnP-algorithms, which has been empirically shown in numerous
papers, rather, we try to present some of the caveats of these algorithms,
when do they fail, and how to improve them. Hence, in most of the examples we
won’t focus on repeating the same experiment over and over with different
parameter values (i.e., tuning the denoising parameter) to attain the best
solution, with the prior knowledge of the true solution. On contrary, for a
fixed selection criterion $\mathcal{S}_{0}$ and a fixed denoiser $H_{\sigma}$,
we present: (1) how can one attain the best solution based on
$\mathcal{S}_{0}$, during the iterative process, and (2) how to further
improve the recovery for the same $H_{\sigma}$ and $\mathcal{S}_{0}$, for
example, when the denoiser is too strong or too weak. We also compare the Fast
FBS-PnP algorithm with the ADMM-PnP algorithm, i.e., the pros and cons of both
these algorithms and, when do they fall in the family $\mathcal{I}_{3}$ or in
$\mathcal{I}_{5}$.

All the experiments are computed in MATLAB, where we consider the selection
criterion ($\mathcal{S}_{0}$) as the cross-validation criterion, for some
leave out set, and the denoiser $H_{\sigma}$ as the BM3D denoiser, where the
denoising parameter $\sigma>0$ corresponds to the standard deviation of the
noise and, the MATLAB code for the BM3D denoiser is obtained from
http://www.cs.tut.fi/ foi/GCF-BM3D/, which is based on [36, 37]. Here, we kept
all the attributes of the code in their original (default) settings, except
the standard deviation, which depends on $\sigma$. Also, the BM3D denoiser
code provided at the above link works for (grayscale or color) images with
intensities in the range [0,1]. Hence, for our problems when the image
intensities fall outside this range, we transform the original denoiser with a
simple rescaling trick, i.e.,
$\hat{H}_{\sigma}(x)=S^{-1}(H_{\sigma}(S(x)),x_{-},x_{+})$, where for any $x$
we have $x_{-}=min(x)$, $x_{+}=max(x)$, $S(x)=\frac{x-x_{-}}{x_{+}-x_{-}}$ and
$S^{-1}(y,x_{-},x_{+})=(x_{+}-x_{-})y+x_{-}$. From now on, we denote
$H_{\sigma}(x)=\hat{H}_{\sigma}(x)$, unless otherwise stated. We consider the
$\ell_{2}$-discrepancy term, $\mathcal{D}(x)=||Ax-b_{\delta}||_{2}^{2}$, where
the noise level ($\delta$) in the data is such that the relative error
$\left(\frac{||b-b_{\delta}||_{2}}{||b||_{2}}\%\right)$ is around 1%, which
corresponds to a 40 dB SNR (signal to noise ratio). Moreover, in a typical
FBS-PnP or ADMM-PnP algorithm, the denoiser is denoted as $H_{\hat{\sigma}}$,
where $\hat{\sigma}=\tau\sigma$ ($\tau$ being the step-size) or
$\hat{\sigma}=\frac{\sigma}{\rho}$ ($\rho$ being the Lagrangian parameter),
respectively. However, to analyze the behavior of the iterative process with
respect to the original denoiser, we instead fix the denoiser $H_{\sigma}$ in
the iterative process, i.e., $H_{\hat{\sigma}}=H_{{\sigma}}$, unless otherwise
stated. As for the evaluation metrics, we compare the PSNR (peak signal-to-
noise ratio) values of the recovered solutions using the MATLAB inbuilt
function $psnr(x_{k}^{\delta},\hat{x},\mbox{max}(\hat{x}))$, the SSIM
(structure similarity index measure) values using MATLAB function
$ssim(S(x_{k}^{\delta}),S(\hat{x}))$, where $S$ is the scaling operator
(defined above), since the MATLAB $ssim()$ function is well-defined for image
intensities in [0,1], and the relative MSE (mean square error) computed as
$\frac{||x_{k}^{\delta}-\hat{x}||_{2}}{||\hat{x}||_{2}}$. The numerical values
are shown in Tables 1, 2, 3 and 4, where $`\mathcal{D}$-err.’ stands for the
discrepancy error
$\left(\frac{||Ax_{k}^{\delta}-b_{D,\delta}||_{2}}{||b_{D,\delta}||_{2}}\right)$
and $`\mathcal{S}$-err.’ stands for the cross-validation error
$\left(\frac{||Ax_{k}^{\delta}-b_{S,\delta}||_{2}}{||b_{S,\delta}||_{2}}\right)$,
where $b_{S,\delta}\subset b_{\delta}$ is the left-out set and
$b_{D,\delta}=b_{\delta}\backslash b_{S,\delta}$, and the recoveries are shown
in Figures 3, 4, 5 and 7.

In the following examples, the matrix equation (1.1) corresponds to the
discretization of a radon transformation, which is associated with the X-ray
computed tomography (CT) reconstruction problem, where we generate the matrix
$A\in\mathbb{R}^{m\times n}$, $\hat{x}\in\mathbb{R}^{n}$ and
$b\in\mathbb{R}^{m}$ from the MATLAB codes presented in [38]. The dimension
$n$ corresponds to the size of a $N\times N$ image, i.e., $n=N^{2}$, and the
dimension $m$ is related to the number of rays per projection angle and the
number of projection angles, i.e., $m=M_{1}\times M_{2}$, where $M_{1}$
implies the number of rays/angle and $M_{2}$ implies the number of angles.

######  Example 5.1.

[Semi-convergence + Weaker denoiser + Boosting]
Here, we show the semi-convergence nature in the recovery errors shown by PnP-
algorithms when an inappropriate denoiser (in this case a weak denoiser) is
used during the minimization process and, how can one retrieve the “best”
possible solution in this case. We also present how to boost the denoising
strength of a weak denoiser $H_{\sigma}$, without altering the $\sigma$ value,
to produce better results. We start with a standard ($128\times 128$) Shepp-
Logan phantom (true image $\hat{x}\in\mathbb{R}^{16384}$ and
$\hat{x}_{i}\in[0,1]$), which is then re-scaled between [-1,1]. The purpose of
rescaling is that, real phantoms are not necessarily restricted to the non-
negativity constraint. Moreover, the constraint ($x_{i}\geq 0$) minimization
corresponding to the standard Shepp-Logan phantom is relatively robust to
noisy data, and since the point of this example is to show the instabilities
occurring in PnP-algorithms, we re-scaled the original phantom. The matrix
$A\in\mathbb{R}^{8064\times 16384}$ is generated using the $PRtomo()$ code
from [38], corresponding to a ‘fancurved’ CT problem with only 45 view angles
(which are evenly spread over $360^{o}$). The noiseless data is generated by
$b:=A\hat{x}\in\mathbb{R}^{181\times 45=8064}$, which is then contaminated by
Gaussian noise ($\epsilon_{\delta}$) to produce noisy data $b_{\delta}$ such
that the relative error is around 1%. Now, we leave out 1% of the noisy data
$b_{\delta}$ for the cross-validation error (the selection criterion
$\mathcal{S}_{0}$), and start the Fast FBS-PnP algorithm with the modified
BM3D denoiser $H_{\sigma}$, with $\sigma=0.0005$, for all the iterations. We
perform the minimization process for two values of the step-size, i.e., (1)
$\tau_{1}=2\times 10^{-4}$ and (2) $\tau_{2}=10^{-5}$, over 1000 iterations
each, to check the convergence or semi-convergence of the recovery errors.

The results are shown in Table 1 and the figures in Figure 3. From Figure
3(e), one can see that the relative errors, when $\tau=2\times 10^{-4}$,
follow a semi-convergence trail, and hence, if the final iterate
($x_{N}^{\delta}$) is considered as the recovered solution, then it’s much
worse than the CV-solution ($x_{k(\delta,\mathcal{S})}^{\delta}$), where
$k(\delta,\mathcal{S})=55$ in this case. Also, note that, by decreasing the
step-size ($\tau=10^{-5}$), we are able to recover a much better (smoother)
estimate, for the reasons explained in §3.1.1. However, it leads to a much
slower descent rate and the recovery errors are still following a convergence
path at the final iteration, i.e., the recovery will improve upon further
iterations. The reason for such a slower recovery rate is the usage of a weak
denoiser, for which, one has to use a smaller step-size (to augment the
denoising).

######  Example 5.2.

[Stronger denoiser + Attenuation + Family $\mathcal{I}_{3}$ or
$\mathcal{I}_{5}$]
Now, we repeat the experiment performed in Example 5.1 but, using a stronger
denoiser $H_{\sigma}$ instead, i.e., a larger value of $\sigma$. Here we set
$\sigma=0.02$ and $\lambda=2\times 10^{-4}$, and perform the Fast FBS-PnP
algorithm for 250 iterations, since $H_{\sigma}$ is a stronger denoiser. Then
we repeat the process by attenuating the denoising strength of $H_{\sigma}$
through the process described in §3.1.2, i.e., using (3.1.2), (3.9) and a
combination of them. Note that, while using (3.1.2), the descent direction is
a function of $\gamma$, i.e., greater values of $\gamma$ implies greater
denoising and vice-verse, as can be seen in Figures 4, in addition to,
$d_{k}^{\delta}(\alpha(\gamma))$ satisfying (2.6). Figure 2(a) shows the graph
of
$\Big{(}d_{k}^{\delta},-\tau\nabla_{x}\mathcal{D}(x_{k}^{\delta})\Big{)}_{2}$
vs. $k$, and Figure 2(b) shows the graph of
$\Big{(}d_{k}^{\delta}(\alpha(\gamma)),-\tau\nabla_{x}\mathcal{D}(x_{k}^{\delta})\Big{)}_{2}$
vs. $k$, for different values of $\gamma$. One can see that $d_{k}^{\delta}$
does not satisfy (2.6) for all values of $k$, where as,
$d_{k}^{\delta}(\alpha(\gamma))$ does satisfy (2.6) for all values of $k$ and
$\gamma$, i.e., without the attenuation the structured directions
$d_{k}^{\delta}$ are not descent directions (due to the strong denoising), but
after the attenuation, the directions $d_{k}^{\delta}(\alpha(\gamma))$ are
descent directions. Where as, the same is not true for
$d_{k}^{\delta}(\alpha_{0})$ generated from (3.9), which can be seen in Figure
2(b), as it does not enforce the condition (2.6) on $d_{k}^{\delta}$. However,
one can see from Figure 2(a), which shows the graph of
$\Big{(}d_{k}^{\delta},-\tau\nabla_{x}\mathcal{D}(x_{k}^{\delta})\Big{)}_{2}$
vs. k, that the structure imposing directions are completely opposite to the
data-consistency directions for larger values of $k$, as the inner-product is
close to -1, suggesting that iterative process is close to an equilibrium
state (i.e., data-consistency and data-denoising steps are negating each
other), which is also visible from Figure 2(b), which shows the graph of
$\Big{(}d_{k}^{\delta}(\alpha_{0}),-\tau\nabla_{x}\mathcal{D}(x_{k}^{\delta})\Big{)}_{2}$
vs. k (highly oscillating between 1 and -1), and Figure 4(g), the MSE curve.
In addition, note that, from Figure 2(c), the $\alpha(\gamma)$ (attenuating)
values for the $d_{k}^{\delta}(\alpha(\gamma))$ steadily decreases to zero,
since the discrepancy values of the iterates $\mathcal{D}(x_{k}^{\delta})$
(and hence, the gradient values) decreases to zero; where as, the $\alpha_{0}$
values does decreases over the iterations but doesn’t tend to zero, as it is
not associated with the condition (2.6), in fact, the oscillations of the
(small) $\alpha_{0}$ values help the iterative process to attain an
equilibrium state; in contrast, the decreasing $\alpha(\gamma)$ values lead to
the semi-convergent nature, where the sharpness of the semi-convergence nature
depends to the $\gamma$ values (smaller $\gamma\Rightarrow$ sharper semi-
convergence). This phenomena is also reflected in the $\mathcal{S}$-error
curves, Figure 2(d), where for the directions $d_{k}^{\delta}(\alpha(\gamma))$
the CV-errors are semi-convergent, but for $d_{k}^{\delta}(\alpha_{0})$, it is
not.

Table 2 shows the error metrics of the recoveries for different descent
directions and different stopping iterations. Figure 4 shows the different
recoveries, corresponding to the $k(\delta,\mathcal{S})$ iterations, the graph
of the PSNR values and the MSE curve. As explained in §3.1, when using
$d_{k}^{\delta}(\alpha(\gamma))$, the iterative process falls in the
regularization family $\mathcal{I}_{3}$, i.e., semi-convergence of the
recovery errors, as can be seen from Figures 4(g), for different values of
$\alpha(\gamma)$. Where as, when using $d_{k}^{\delta}(\alpha_{0})$, the
recovery process falls in the regularization family $\mathcal{I}_{5}$, i.e.,
may not exhibit semi-convergence of the recovery errors, as can be seen from
Figures 4(g). Therefore, it suggests that, when the recovery process belongs
to the family $\mathcal{I}_{3}$, i.e., all the $d_{k}^{\delta}$s are descent
directions, then the minimization process can be terminated at a much earlier
instance and a recovered solution $x_{k(\delta)}^{\delta}$ can correspond to
the one satisfying the selection criterion the best, where as, if recovery
process falls in family $\mathcal{I}_{5}$, then one can wait for longer period
of iterations and can recover better solutions, as the fluctuating/oscillating
iterates $z_{k}^{\delta}$s can produce better approximations over time.

######  Example 5.3.

[ADMM - small vs. large $\rho$ \+ with vs. without $\sigma$ update]
In this example we show the importance of the Lagrangian parameter ($\rho$)
and the influence of the update or no-update of the denoiser’s denoising
strength ($\sigma$-value), during the iterative process, i.e., (1st) we
perform the traditional ADMM-PnP algorithm, where $H_{\sigma}$ is updated to
$H_{\hat{\sigma}}$, where $\hat{\sigma}=\frac{\sigma}{\rho}$, during the
iterative process, and (2nd) we perform the iterative process without updating
$H_{\sigma}$ during the descent process. In addition, for each of the two
experiments, we repeat the processes for different values of the Lagrangian
parameter ($\rho$), to observe its effect on the recovered solution. Again,
for the classical Tikhonov-type regularization, where the regularization
function $\mathcal{R}(x)$ in (1.3) is proper, closed and convex, the (1st)
method is guaranteed to converge to the minimizer of (1.3), irrespective of
the $\rho$-value, although, the minimizer might not be the best estimate of
the true solution $\hat{x}$ during the minimization process, for a fixed
$\lambda$ value. Note that, as described in §4.2, for smaller $\rho$-value the
first iterate $x_{1}^{\delta}$ is close to the noisy LS-solution
($x^{\dagger}_{\delta}$), associated with the noisy data $b_{\delta}$. Hence,
the denoiser $H_{\sigma}$ needs to be particularly good so that, after
updating to $H_{\hat{\sigma}}$, it should be able to clean the large
corruptions generated by the noisy data. Where as, for large $\rho$,
$x_{1}^{\delta}$ is close to $x_{0}^{\delta}$, and thus, will be smooth for
smoother $x_{0}^{\delta}$. Hence, one may clean the artifacts (corruptions)
steadily via a moderately strong denoiser. However, if $H_{\sigma}$ is updated
to $H_{\hat{\sigma}=\frac{\sigma}{\rho}}$, then the denoising ability of
$H_{\sigma}$ is also reduced, by a factor corresponding to $\frac{1}{\rho}$,
and thus, a moderately strong denoiser becomes a weaker denoiser, and may not
be able to clean the corruptions. Therefore, for the fixed denoiser
$H_{\sigma}$, (2nd) method suggests to incorporate the denoiser $H_{\sigma}$
without updating it, to preserve the denoising strength of a moderately strong
denoiser and provide effective denoising. Moreover, since the convergence of
ADMM-PnP algorithm for any general denoiser $H_{\sigma}$ is not known, it does
not hurt to implement it in either way; after-all, both these processes can be
thrown in either $\mathcal{I}_{5}$ or $\mathcal{I}_{6}$ family and the better
process is the one for which the recovered solution satisfies the selection
criterion ($\mathcal{S}$) the best.

Here we keep the experimental settings similar to Example 5.1, but the
iterative process follows the ADMM-PnP algorithm, using the stronger denoiser
$H_{\sigma}$, with $\sigma=0.02$. For each of the above two strategies, we
repeat the experiments for four values of the Lagrangian parameter,
$\rho=0.01,0.1,10,100$, for over 250 iterations. Also, the minimizer of each
$x$-update, i.e., $x_{k}^{\delta}$ as defined in (4.13), is only estimated by
a certain number of conjugate-gradient (CG) steps and not through the matrix
inversion. Here, we considered 100 CG-steps, initiated from
$z_{k-1}^{\delta}-u_{k-1}^{\delta}$, to approximate the k-th step x-minimizer,
i.e., to generate the $x_{k}^{\delta}$ iterate. Hence, the 100 CG-steps, in
combination with the 250 ADMM-steps, leads to a total of 250,000 iterations,
although, there are only 250 denoising steps (which can be computationally
expensive). Also, note that, the denoiser $H_{\sigma}$ for $\sigma=0.02$ acted
as a strong denoiser in Example 5.2, where as, here it behaves like a weak
denoiser, except when $\rho$ is large and $\sigma$ is not updated to
$\hat{\sigma}$. Interestingly, even when $\sigma$ is updated to
$\hat{\sigma}=\frac{\sigma}{\rho}>>\sigma$ (for smaller $\rho$), the denoiser
$H_{\hat{\sigma}}$ is still not adequately strong to clean the corruptions in
the iterates as effectively as $H_{\sigma}$ with larger $\rho$ value, which
validates the aforementioned explanations. Table 3 shows the error metrics of
the recoveries corresponding to the minimum CV-error, Figure 5 shows the
recoveries and, Figure 6 shows the MSE curves of the recoveries.

######  Example 5.4.

In this example we compare the recoveries obtained via the precondition matrix
$L=|\nabla|$, for the first iterate ($x_{1}^{\delta}$) in (4.18), with the
recoveries using $L=I$. We repeat the experiment in Example 5.3, for
$\rho=100$, but using $L=|\nabla|$ to calculate $x_{1}^{\delta}$. The CV-
solution ($x_{k(\delta,\mathcal{S})}^{\delta}$), with and without $\sigma$
update, for $L=I$ is shown in Figures 5(g) and 5(h), and the CV-solutions for
$L=|\nabla|$ is shown in Figure 7. Note that, the recoveries corresponding to
$H_{\hat{\sigma}}$ is bad in both the cases and the recovery using
$L=|\nabla|$ is slightly better than $L=I$, see Table 4.

##  6\. Conclusion and Future Research

In this paper we tried to explain the PnP-algorithms from a different angle.
We defined certain families of regularized solutions and showed that the
solution of a PnP-algorithm will fall in one of them. We started with an
extension of the Landweber iterations to generate structure imposing descent
directions, i.e., any direction $d_{k}^{\delta}$ satisfying (2.6). This lead
to the formulation of the regularization family $\mathcal{I}_{3}$, which
contains all the semi-iterative processes. This is further generalized to the
family $\mathcal{I}_{5}$, where the directions $d_{k}^{\delta}$ satisfies
(2.6) for certain number of iterations and then the restriction is relaxed,
i.e., once the iterates $x_{k}^{\delta}$ have approximated the data
$b_{\delta}$ to a certain extent
($\mathcal{D}(x_{k}^{\delta})\leq\epsilon_{1}(\delta)$), we don’t want
$d_{k}^{\delta}$ to be a descent direction anymore, as then, it will lead to
semi-convergence in the recovery errors. The advantage of such a formulation
is that, in addition to avoiding the semi-convergence of the recovery errors,
one doesn’t even have to worry about the convergence of the iterates
$x_{k}^{\delta}$ anymore, i.e., the iterative process doesn’t necessarily need
to be associated with a penalized/constrained Variational-minimization
problem, and thus, the iterates $x_{k}^{\delta}$ do not need to converge to a
minimizer. Hence, it compensates the short-comings of both the classical
regularization methods, i.e., the semi-iterative methods and the Tikhonov-type
methods. In other words, any iterate $x_{k}^{\delta}$ during the minimization
process is an “approximate estimate” to the solution of the inverse problem
($x^{\dagger}$), and the appropriateness of the estimate is determined via a
selection criterion $\mathcal{S}_{0}$ for the solution of the inverse problem
(1.1). Therefore, one can terminate the iterative process at an early
instance, if the error in the selection criterion $\mathcal{S}_{0}$ starts
getting worse, and hence, can reduce the computational time significantly; or
can wait for longer period of iterations, if $\mathcal{S}(x_{k}^{\delta})$
values is decreasing, assuming the selection criterion $\mathcal{S}_{0}$ is
appropriate for the problem. Note that, one does not necessarily need a
Variational formulation to regularize the solution of an inverse problem. A
“regularized solution $x^{\delta}$” is a solution which “approximates
$\hat{x}$”, depending on the noise level $\delta$, and “gets better” as “noise
vanishes”, i.e., $\delta\rightarrow 0$ should imply
$||x^{\delta}-\hat{x}||\rightarrow 0$, and “$x^{\delta}$ is well-defined” if
there is a “unique process” of obtaining it, for example, either via a
Variational formulation or via an iterative formulation. Hence, this leads to
the formulation of our largest family of regularized solution
$\mathcal{I}_{5}$, as defined in (4.9), i.e., a solution state
$x^{\delta}(.,t)$ corresponding to an initial value differential equation
problem.

Furthermore, we showed how to improve the recoveries when dealing with an
inappropriate denoiser $H_{\sigma}$, without altering the denoiser, i.e., when
the denoiser is too strong to oversmooth the recovered solution, then one can
attenuate the denoising strength via a relaxing parameter, as described in
§3.1.2, or, when the denoiser is very weak to remove the corruptions properly,
then one can reduce the step-size to boost the denoising strength, as
described in §3.1.1. We also showed the importance of the Lagrangian parameter
($\rho$) in the ADMM-PnP algorithm and the influence of updating or not
updating the denoiser $H_{\sigma}$ to $H_{\hat{\sigma}}$, where
$\hat{\sigma}=\frac{\sigma}{\rho}$, during the iterative process, as well, as
the significance of the preconditioned matrix $L$ (for the first iterate
$x_{1}^{\delta}$) in (4.18). The validations of these improving techniques can
be seen in Examples 5.1, 5.2, 5.3 and 5.4.

In an upcoming paper, we extend this idea to the setting of training examples
or data, i.e., for a neural network denoiser or a denoising iterative scheme
(such as unrolled neural network). We believe that by incorporating this idea
of regularization, i.e., monitoring the recovery process via an appropriate
selection criterion, one can understand and control the instabilities arising
in deep learning based reconstruction algorithms, as shown in [35].

Fast FBS-PnP + a weak denoiser $H_{\sigma}$ ($\sigma=0.0005$), for
$\tau=2\times 10^{-4}$ vs. $10^{-5}$
---
step-size |  iter.(k) |  MSE |  $\mathcal{D}$-err. |  $\mathcal{S}$-err. |  PSNR |  SSIM |  Min.MSE
$\tau=2\times 10^{-4}$ |  1000(N) |  0.3432 |  0.3264 |  0.0290 |  10.51 |  0.2926 |  0.0208 (50)
|  55 ($\mathcal{S}$) |  0.2088 |  0.4891 |  0.0248 |  14.83 |  0.2540 |
$\tau=10^{-5}$ |  1000(N) |  0.0899 |  0.0075 |  0.0134 |  22.15 |  0.3932 |  0.0899(1000)
|  1000($\mathcal{S}$) |  0.0899 |  0.0075 |  0.0134 |  22.15 |  0.3932 |
Table 1. Comparing solution $z_{N}^{\delta}$ vs.
$z_{k(\delta,\mathcal{S})}^{\delta}$ for Example 5.1. Fast FBS-PnP + a strong
denoiser $H_{\sigma}$ ($\sigma=0.01$) + $d_{k}^{\delta}$ vs.
$d_{k}^{\delta}(\alpha(\gamma))$ vs. $d_{k}^{\delta}(\alpha_{0})$
---
directions |  iter.(k) |  MSE |  $\mathcal{D}$-err. |  $\mathcal{S}$-err. |  PSNR |  SSIM |  Min.MSE
$d_{k}^{\delta}$ |  250(N) |  0.1491 |  0.0174 |  0.0182 |  17.75 |  0.3991 |  0.1477(36)
|  36($\mathcal{S}$) |  0.1477 |  0.0171 |  0.0176 |  17.83 |  0.3750 |
$d_{k}^{\delta}(\alpha(\gamma=0.1))$ |  250(N) |  0.2312 |  0.0033 |  0.0271 |  13.94 |  0.2284 |  0.1660(52)
|  36($\mathcal{S}$) |  0.1682 |  0.0064 |  0.0171 |  16.71 |  0.2903 |
$d_{k}^{\delta}(\alpha(\gamma=0.5))$ |  250(N) |  0.1846 |  0.0034 |  0.0245 |  15.90 |  0.2730 |  0.1181(75)
|  46($\mathcal{S}$) |  0.1242 |  0.0055 |  0.0156 |  19.34 |  0.3395 |
$d_{k}^{\delta}(\alpha(\gamma=0.9))$ |  250(N) |  0.1433 |  0.0036 |  0.0215 |  18.10 |  0.3046 |  0.0823(69)
|  63($\mathcal{S}$) |  0.0832 |  0.0057 |  0.0122 |  22.82 |  0.3697 |
$d_{k}^{\delta}(\alpha(\gamma=1))$ |  250(N) |  0.0729 |  0.0048 |  0.0142 |  23.97 |  0.3673 |  0.0542(102)
|  71($\mathcal{S}$) |  0.0648 |  0.0073 |  0.0121 |  25.00 |  0.3928 |
$d_{k}^{\delta}(\alpha_{0})$ |  250(N) |  0.0455 |  0.0066 |  0.0122 |  28.07 |  0.4106 |  0.0440(131)
|  192($\mathcal{S}$) |  0.0467 |  0.0066 |  0.0118 |  27.83 |  0.4039 |
Table 2. Comparing solution $z_{N}^{\delta}$ vs.
$z_{k(\delta,\mathcal{S})}^{\delta}$ for Example 5.2. ADMM-PnP using
$H_{\sigma}$, for $\sigma=0.02$, vs. $H_{\hat{\sigma}}$, where
$\hat{\sigma}=\sigma/\rho$
---
denoiser and $\rho$ |  iter.(k) |  MSE |  $\mathcal{D}$-err. |  $\mathcal{S}$-err. |  PSNR |  SSIM |  Min.MSE
$H_{\hat{\sigma}}$ & $\rho=0.01$ |  250(N) |  0.2100 |  0.0044 |  0.0239 |  14.78 |  0.2467 |  0.2088(3)
|  219($\mathcal{S}$) |  0.2108 |  0.0043 |  0.0232 |  14.75 |  0.2472 |
$H_{{\sigma}}$ & $\rho=0.01$ |  250(N) |  01.07 |  0.0026 |  0.0776 |  0.6266 |  0.1307 |  0.2088(3)
|  3($\mathcal{S}$) |  0.2088 |  0.0048 |  0.0249 |  14.83 |  0.2533 |
$H_{\hat{\sigma}}$ & $\rho=0.1$ |  250(N) |  0.1347 |  0.0064 |  0.0149 |  18.64 |  0.3289 |  0.1346(172)
|  248($\mathcal{S}$) |  0.1349 |  0.0064 |  0.0148 |  18.63 |  0.3284 |
$H_{{\sigma}}$ & $\rho=0.1$ |  250(N) |  0.9147 |  0.0026 |  0.0671 |  2.00 |  0.2115 |  0.2087(3)
|  3($\mathcal{S}$) |  0.2087 |  0.0048 |  0.0249 |  14.83 |  0.2537 |
$H_{\hat{\sigma}}$ & $\rho=10$ |  250(N) |  0.4225 |  0.0029 |  0.0373 |  08.71 |  0.1589 |  0.2087(4)
|  4($\mathcal{S}$) |  0.2087 |  0.0047 |  0.0247 |  14.83 |  0.2533 |
$H_{{\sigma}}$ & $\rho=10$ |  250(N) |  0.1428 |  0.0044 |  0.0159 |  18.13 |  0.3434 |  0.1379(74)
|  214($\mathcal{S}$) |  0.1403 |  0.0044 |  0.0153 |  18.28 |  0.3388 |
$H_{\hat{\sigma}}$ & $\rho=100$ |  250(N) |  0.2754 |  0.0032 |  0.0294 |  12.42 |  0.1954 |  0.2082(8)
|  10($\mathcal{S}$) |  0.2083 |  0.0049 |  0.0246 |  14.85 |  0.2548 |
$H_{{\sigma}}$ & $\rho=100$ |  250(N) |  0.0406 |  0.0087 |  0.0121 |  29.06 |  0.4130 |  0.0405(96)
|  38($\mathcal{S}$) |  0.0540 |  0.0087 |  0.0120 |  26.58 |  0.4074 |
Table 3. Comparing solution $z_{N}^{\delta}$ vs.
$z_{k(\delta,\mathcal{S})}^{\delta}$ for Example 5.3. ADMM-PnP, for
$\sigma=0.02$ and $\rho=100$, $H_{\sigma}$ vs. $H_{\hat{\sigma}}$ and $L=I$
vs. $L=|\nabla|$
---
denoiser and $L$ |  iter.(k) |  MSE |  $\mathcal{D}$-err. |  $\mathcal{S}$-err. |  PSNR |  SSIM |  Min.MSE
$H_{\hat{\sigma}}$ & $L=I$ |  250(N) |  0.2754 |  0.0032 |  0.0294 |  12.42 |  0.1954 |  0.2082(8)
|  10($\mathcal{S}$) |  0.2083 |  0.0049 |  0.0246 |  14.85 |  0.2548 |
$H_{{\sigma}}$ & $L=I$ |  250(N) |  0.0406 |  0.0087 |  0.0121 |  29.06 |  0.4130 |  0.0405(96)
|  38($\mathcal{S}$) |  0.0540 |  0.0087 |  0.0120 |  26.58 |  0.4074 |
$H_{\hat{\sigma}}$ & $L=|\nabla|$ |  250(N) |  0.2750 |  0.0032 |  0.0291 |  12.44 |  0.1957 |  0.2073(8)
|  8($\mathcal{S}$) |  0.2073 |  0.0051 |  0.0235 |  14.89 |  0.2572 |
$H_{{\sigma}}$ & $L=|\nabla|$ |  250(N) |  0.0404 |  0.0087 |  0.0121 |  29.10 |  0.4131 |  0.0401(79)
|  67($\mathcal{S}$) |  0.0404 |  0.0087 |  0.0120 |  29.09 |  0.4074 |
Table 4. Comparing solution $z_{N}^{\delta}$ vs.
$z_{k(\delta,\mathcal{S})}^{\delta}$ for Example 5.4.

(a)
$\Big{(}d_{k}^{\delta}(\alpha(\gamma)),-\tau\nabla\mathcal{D}(x_{k}^{\delta})\Big{)}_{2}$
before attenuating the denoiser’s strength

(b)
$\Big{(}d_{k}^{\delta}(\alpha(\gamma)),-\tau\nabla\mathcal{D}(x_{k}^{\delta})\Big{)}_{2}$
after attenuating the denoiser’s strength

(c) $\alpha(\gamma)$ (attenuating) values over iterations

(d) CV-errors for different $\gamma$ values and $\alpha_{0}$

Figure 2.  Attenuation of denoising strength, see Example 5.2.

(a) True Phantom.

(b) $z_{k(\delta,\mathcal{S})}^{\delta}$, $k(\delta,\mathcal{S})=1000$ and
$\tau=10^{-5}$.

(c) $z_{N}^{\delta}$, $N=1000$ and $\tau=2\times 10^{-4}$.

(d) $z_{k(\delta,\mathcal{S})}^{\delta}$, $k(\delta,\mathcal{S})=55$ and
$\tau=2\times 10^{-4}$.

(e) MSE vs. $k$, for $\tau=2\times 10^{-4}$ and $10^{-5}$.

(f) PSNR vs. $k$, for $\tau=2\times 10^{-4}$ and $10^{-5}$.

Figure 3. $z_{N}^{\delta}$ vs. $z_{k,\mathcal{S}}^{\delta}$ and performance
curves for Example 5.1.

(a) $z_{k(\delta,\mathcal{S})}$, for $k(\delta,\mathcal{S})=36$ and
$d_{k}^{\delta}$

(b) $z_{k(\delta,\mathcal{S})}$, for $k(\delta,\mathcal{S})=36$ &
$d_{k}^{\delta}(\alpha(\gamma=0.1))$

(c) $z_{k(\delta,\mathcal{S})}$, for $k(\delta,\mathcal{S})=36$ &
$d_{k}^{\delta}(\alpha(\gamma=0.5))$

(d) $z_{k(\delta,\mathcal{S})}$, for $k(\delta,\mathcal{S})=36$ &
$d_{k}^{\delta}(\alpha(\gamma=0.9))$

(e) $z_{k(\delta,\mathcal{S})}$, for $k(\delta,\mathcal{S})=36$ &
$d_{k}^{\delta}(\alpha(\gamma=1))$

(f) $z_{k(\delta,\mathcal{S})}$, for $k(\delta,\mathcal{S})=36$ &
$d_{k}^{\delta}(\alpha_{0})$

(g) MSE curves for $d_{k}^{\delta}$, $d_{k}^{\delta}(\alpha_{0})$ and
$d_{k}^{\delta}(\gamma)$

(h) PSNR curves for $d_{k}^{\delta}$, $d_{k}^{\delta}(\alpha_{0})$ and
$d_{k}^{\delta}(\alpha(\gamma))$

Figure 4. CV-recoveries and performance curves for Example 5.2.

(a) $x_{k(\delta)}^{\delta}$ for $H_{\hat{\sigma}}$ and $\rho=0.01$

(b) $x_{k(\delta)}^{\delta}$ for $H_{{\sigma}}$ and $\rho=0.01$

(c) $x_{k(\delta)}^{\delta}$ for $H_{\hat{\sigma}}$ and $\rho=0.1$

(d) $x_{k(\delta)}^{\delta}$ for $H_{\sigma}$ and $\rho=0.1$

(e) $x_{k(\delta)}^{\delta}$ for $H_{\hat{\sigma}}$ and $\rho=10$

(f) $x_{k(\delta)}^{\delta}$ for $H_{\sigma}$ and $\rho=10$

(g) $x_{k(\delta)}^{\delta}$ for $H_{\hat{\sigma}}$ and $\rho=100$

(h) $x_{k(\delta)}^{\delta}$ for $H_{\sigma}$ and $\rho=100$

Figure 5. $H_{\hat{\sigma}=\sigma/\rho}$ vs. $H_{\sigma}$ for ADMM-PnP, see
Example 5.3.

(a) MSE curves for $H_{\hat{\sigma}}$, where $\hat{\sigma}=\sigma/\rho$

(b) MSE curves for $H_{{\sigma}}$

Figure 6. Recovery errors for $H_{\hat{\sigma}}$ vs. $H_{\sigma}$, see Example
5.3.

(a) $x_{k(\delta,\mathcal{S})}^{\delta}$, for $H_{\hat{\sigma}}$ and
$L=|\nabla|$ for $\rho=100$

(b) $x_{k(\delta,\mathcal{S})}^{\delta}$, for $H_{{\sigma}}$ and $L=|\nabla|$
for $\rho=100$

Figure 7. $x_{k(\delta,\mathcal{S})}^{\delta}$ with precondition matrix
$L=|\nabla|$, see Example 5.4.

## References

  * [1] H. W. Engl, M. Hanke, and A. Neubauer, Regularization of inverse problems, vol. 375 of Mathematics and its Applications.  Kluwer Academic Publishers Group, Dordrecht, 1996.
  * [2] A. Bakushinsky and A. Goncharsky, Ill-posed problems: theory and applications, vol. 301 of Mathematics and its Applications.  Kluwer Academic Publishers Group, Dordrecht, 1994.  Translated from the Russian by I. V. Kochikov.
  * [3] C. W. Groetsch, The theory of Tikhonov regularization for Fredholm equations of the first kind, vol. 105 of Research Notes in Mathematics.  Pitman (Advanced Publishing Program), Boston, MA, 1984.
  * [4] J. Baumeister, Stable solution of inverse problems.  Advanced Lectures in Mathematics, Friedr. Vieweg & Sohn, Braunschweig, 1987.
  * [5] V. A. Morozov, Methods for solving incorrectly posed problems.  Springer-Verlag, New York, 1984.  Translated from the Russian by A. B. Aries, Translation edited by Z. Nashed.
  * [6] S. Osher, M. Burger, D. Goldfarb, J. Xu, and W. Yin, “An iterative regularization method for total variation-based image restoration,” Multiscale Modeling & Simulation, vol. 4, no. 2, pp. 460–489, 2005.
  * [7] S. Ma, W. Yin, Y. Zhang, and A. Chakraborty, “An efficient algorithm for compressed mr imaging using total variation and wavelets,” in 2008 IEEE Conference on Computer Vision and Pattern Recognition, pp. 1–8, 2008.
  * [8] J. Yang, Y. Zhang, and W. Yin, “A fast alternating direction method for tvl1-l2 signal reconstruction from partial fourier data,” IEEE Journal of Selected Topics in Signal Processing, vol. 4, no. 2, pp. 288–297, 2010.
  * [9] S. Ravishankar and Y. Bresler, “Mr image reconstruction from highly undersampled k-space data by dictionary learning,” IEEE Transactions on Medical Imaging, vol. 30, no. 5, pp. 1028–1041, 2011.
  * [10] H. Y. Liao and G. Sapiro, “Sparse representations for limited data tomography,” in 2008 5th IEEE International Symposium on Biomedical Imaging: From Nano to Macro, pp. 1375–1378, 2008.
  * [11] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sciences, vol. 2, pp. 183–202, 01 2009.
  * [12] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Found. Trends Mach. Learn., vol. 3, p. 1–122, Jan. 2011\.
  * [13] A. Chambolle and T. Pock, “A first-order primal-dual algorithm for convex problems with applications to imaging.,” Journal of Mathematical Imaging and Vision, vol. 40, no. 1, pp. 120–145, 2011.
  * [14] S. V. Venkatakrishnan, C. A. Bouman, and B. Wohlberg, “Plug-and-play priors for model based reconstruction,” in 2013 IEEE Global Conference on Signal and Information Processing, pp. 945–948, 2013.
  * [15] S. H. Chan, X. Wang, and O. A. Elgendy, “Plug-and-play admm for image restoration: Fixed-point convergence and applications,” IEEE Transactions on Computational Imaging, vol. 3, no. 1, pp. 84–98, 2017.
  * [16] G. T. Buzzard, S. H. Chan, S. Sreehari, and C. A. Bouman, “Plug-and-play unplugged: Optimization-free reconstruction using consensus equilibrium,” SIAM Journal on Imaging Sciences, vol. 11, no. 3, pp. 2001–2020, 2018.
  * [17] E. Ryu, J. Liu, S. Wang, X. Chen, Z. Wang, and W. Yin, “Plug-and-play methods provably converge with properly trained denoisers,” in Proceedings of the 36th International Conference on Machine Learning (K. Chaudhuri and R. Salakhutdinov, eds.), vol. 97 of Proceedings of Machine Learning Research, pp. 5546–5557, PMLR, 09–15 Jun 2019.
  * [18] A. M. Teodoro, J. M. Bioucas-Dias, and M. A. T. Figueiredo, “A convergent image fusion algorithm using scene-adapted gaussian-mixture-based denoising,” IEEE Transactions on Image Processing, vol. 28, no. 1, pp. 451–463, 2019.
  * [19] Y. Sun, B. E. Wohlberg, and U. Kamilov, “An online plug-and-play algorithm for regularized image reconstruction,” IEEE Transactions on Computational Imaging, vol. 5, 1 2019.
  * [20] Y. Romano, M. Elad, and P. Milanfar, “The little engine that could: Regularization by denoising (red),” SIAM Journal on Imaging Sciences, vol. 10, no. 4, pp. 1804–1844, 2017.
  * [21] J. Liu, Y. Sun, C. Eldeniz, W. Gan, H. An, and U. S. Kamilov, “Rare: Image reconstruction using deep priors learned without groundtruth,” IEEE Journal of Selected Topics in Signal Processing, vol. 14, no. 6, pp. 1088–1099, 2020.
  * [22] L. Landweber, “An iteration formula for fredholm integral equations of the first kind,” American Journal of Mathematics, vol. 73, no. 3, pp. 615–624, 1951.
  * [23] M. Hanke, “Accelerated landweber iterations for the solution of ill-posed equations,” Numerische Mathematik, vol. 60, pp. 341–373, Dec 1991.
  * [24] M. Hanke, A. Neubauer, and O. Scherzer, “A convergence analysis of the landweber iteration for nonlinear ill-posed problems,” Numerische Mathematik, vol. 72, pp. 21–37, Nov 1995.
  * [25] V. A. Morozov, “On the solution of functional equations by the method of regularization,” Soviet Math. Dokl., vol. 7, pp. 414–417, 1966.
  * [26] G. M. Vaĭnikko, “The principle of the residual for a class of regularization methods,” Zh. Vychisl. Mat. i Mat. Fiz., vol. 22, no. 3, pp. 499–515, 763, 1982.
  * [27] H. Gfrerer, “An a posteriori parameter choice for ordinary and iterated Tikhonov regularization of ill-posed problems leading to optimal convergence rates,” Math. Comp., vol. 49, no. 180, pp. 507–522, S5–S12, 1987.
  * [28] P. C. Hansen, “Analysis of discrete ill-posed problems by means of the l-curve,” SIAM Rev., vol. 34, no. 4, pp. 561–580, 1992.
  * [29] C. L. Lawson and R. J. Hanson, Solving least squares problems.  Prentice-Hall, Inc., Englewood Cliffs, N.J., 1974.  Prentice-Hall Series in Automatic Computation.
  * [30] U. Tautenhahn and U. Hämarik, “The use of monotonicity for choosing the regularization parameter in ill-posed problems,” Inverse Problems, vol. 15, no. 6, pp. 1487–1505, 1999.
  * [31] F. Bauer and T. Hohage, “A Lepskij-type stopping rule for regularized Newton methods,” Inverse Problems, vol. 21, no. 6, pp. 1975–1991, 2005\.
  * [32] P. Mathé, “The Lepskiĭ principle revisited,” Inverse Problems, vol. 22, no. 3, pp. L11–L15, 2006.
  * [33] V. Morozov, “The error principle in the solution of operational equations by the regularization method,” USSR Computational Mathematics and Mathematical Physics, vol. 8, no. 2, pp. 63 – 87, 1968.
  * [34] B. Zhu, J. Z. Liu, B. Rosen, and M. Rosen, “Image reconstruction by domain transform manifold learning,” Nature, vol. 555, 03 2018.
  * [35] V. Antun, F. Renna, C. Poon, B. Adcock, and A. C. Hansen, “On instabilities of deep learning in image reconstruction and the potential costs of ai,” Proceedings of the National Academy of Sciences, vol. 117, no. 48, pp. 30088–30095, 2020.
  * [36] Y. Mäkinen, L. Azzari, and A. Foi, “Exact transform-domain noise variance for collaborative filtering of stationary correlated noise,” in 2019 IEEE International Conference on Image Processing (ICIP), pp. 185–189, 2019.
  * [37] Y. Mäkinen, L. Azzari, and A. Foi, “Collaborative filtering of correlated noise: Exact transform-domain variance for improved shrinkage and patch matching,” IEEE Transactions on Image Processing, vol. 29, pp. 8339–8354, 2020.
  * [38] S. Gazzola, P. Hansen, and J. Nagy, “Ir tools - a matlab package of iterative regularization methods and large-scale test problems,” Numerical Algorithms, 2018.