# Convergent plug-and-play with proximal denoiser and unconstrained regularization parameter

Samuel Hurault, Antonin Chambolle, Arthur Leclaire, Nicolas Papadakis

<figure><img src="image_1.png" /></figure>

## To cite this version:

Samuel Hurault, Antonin Chambolle, Arthur Leclaire, Nicolas Papadakis. Convergent plug-and-play with proximal denoiser and unconstrained regularization parameter. Journal of Mathematical Imaging and Vision, 2024, 66, pp.616-638. 10.1007/s10851-024-01195-w. hal-04269033

HAL Id: hal-04269033

https://hal.science/hal-04269033v1

Submitted on 4 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Convergent plug-and-play with proximal denoiser and unconstrained regularization parameter.

Samuel Hurault

samuel.hurault@math.u-bordeaux.fr

University of Bordeaux

Antonin Chambolle

Paris Dauphine University

Arthur Leclaire

University of Bordeaux

Nicolas Papadakis

University of Bordeaux

## Research Article

**Keywords:** Nonconvex optimization, inverse problems, plug-and-play

**Posted Date:** October 23rd, 2023

**DOI:** https://doi.org/10.21203/rs.3.rs-3462535/v1

**License:** ©① This work is licensed under a Creative Commons Attribution 4.0 International License.

Read Full License

**Additional Declarations:** No competing interests reported.

**Version of Record:** A version of this preprint was published at Journal of Mathematical Imaging and Vision on June 12th, 2024. See the published version at https://doi.org/10.1007/s10851-024-01195-w.
Convergent plug-and-play with proximal denoiser
and unconstrained regularization parameter

Samuel Hurault¹, Antonin Chambolle², Arthur Leclaire¹,
Nicolas Papadakis¹

¹*Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, F-33400,
Talence, France.

²CEREMADE, CNRS, Université Paris-Dauphine, PSL, Palaiseau,
91128, France.

Contributing authors: samuel.hurault@math.u-bordeaux.fr;
antonin.chambolle@ceremade.dauphine.fr; arthur.leclaire@math.cnrs.fr;
nicolas.papadakis@math.u-bordeaux.fr;

Abstract

In this work, we present new proofs of convergence for Plug-and-Play (PnP)
algorithms. PnP methods are efficient iterative algorithms for solving image
inverse problems where regularization is performed by plugging a pre-trained
denoiser in a proximal algorithm, such as Proximal Gradient Descent (PGD) or
Douglas-Rachford Splitting (DRS). Recent research has explored convergence by
incorporating a denoiser that writes exactly as a proximal operator. However, the
corresponding PnP algorithm has then to be run with stepsize equal to 1. The
stepsize condition for nonconvex convergence of the proximal algorithm in use
then translates to restrictive conditions on the regularization parameter of the
inverse problem. This can severely degrade the restoration capacity of the algo-
rithm. In this paper, we present two remedies for this limitation. First, we provide
a novel convergence proof for PnP-DRS that does not impose any restrictions on
the regularization parameter. Second, we examine a relaxed version of the PGD
algorithm that converges across a broader range of regularization parameters.
Our experimental study, conducted on deblurring and super-resolution experi-
ments, demonstrate that both of these solutions enhance the accuracy of image
restoration.

**Keywords:** Nonconvex optimization, inverse problems, plug-and-play
# 1 Introduction

In this work, we study the convergence of Plug-and-Play algorithms for solving image inverse problems. This is achieved by minimizing an explicit function:

$$ \hat{x} \in \arg \min_{x} \lambda f(x) + \phi(x). \qquad (1) $$

The function $f$ represents a data-fidelity term that quantifies the alignment of the estimation $x$ with the degraded observation $y$, the function $\phi$ is a nonconvex regularization function, and $\lambda > 0$ acts as a parameter, determining the balance between these two terms.

In our experimental setup, we focus on degradation models where $y$ is expressed as $y = Ax^* + \nu$, with $x^*$ representing the ground-truth signal in $\mathbb{R}^n$, $A$ a linear operator in $\mathbb{R}^{n \times m}$, and $\nu$ a white Gaussian noise in $\mathbb{R}^n$. In this context, the data-fidelity term $f(x)$ takes the form of a convex and smooth function $f(x) = \frac{1}{2}\|Ax - y\|^2$. However, our theoretical analysis extends to a wider range of noise models, as well as data-fidelity terms $f$ that may not be smooth or convex.

To find an adequate solution of the ill-posed problem of recovering $x^*$ from $y$, the choice of the regularization $\phi$ is crucial. Research has first been dedicated to find, by hand, regularizing potentials $\phi$ that are minimal when a desired property of a clean image is satisfied. Among the most classical priors, one can single out total variation [1], $L^2$ gradient norm (Tikhonov) [2], or sparsity in a given dictionary, for example in the wavelet representation [3]. These handcrafted priors are now largely outperformed by learning approaches [4–7], which may not even be associated to a closed-form regularization function $\phi$.

Estimating a local or global optima of problem (1) is classically done using proximal splitting optimization algorithms such as Proximal Gradient Descent (PGD), Douglas-Rashford Splitting (DRS) or Primal-Dual. Given an adequate stepsize $\tau > 0$, these methods alternate between explicit gradient descent steps, $\text{Id} - \tau \nabla h$ for smooth functions $h$, and/or implicit gradient steps using the proximal operator $\text{Prox}_{\tau h}(x) \in \arg \min_z \frac{1}{2\tau} \|z - x\|^2 + h(z)$, for a proper lower semi-continuous function $h$. Proximal algorithms are originally designed for convex functions, but under appropriate assumptions, PGD [8] and DRS [9] algorithms converge to a stationary point of problem (1) associated to nonconvex functions $f$ and $\phi$.

## Plug-and-Play algorithms

Plug-and-Play (PnP) [10] and Regularization by Denoising (RED) [6] methods consist in splitting algorithms in which the descent step on the regularization function is performed by an off-the-shelf image denoiser. They are respectively built from proximal splitting schemes by replacing the proximal operator (PnP) or the gradient operator (RED) of the regularization $\phi$ by an image denoiser. When used with a deep denoiser (i.e parameterized by a neural network) these approaches produce impressive results for various image restoration tasks [7, 11–16].

Theoretical convergence of PnP and RED algorithms with deep denoisers has recently been addressed by a variety of studies [16–19]. Most of these works require
specific constraints on the deep denoiser, such as nonexpansivity. However, imposing nonexpansivity of a denoiser can severely degrade its performance ([19, Table 1], [20, Figure 1], [21, Table 3]).

Another line of works [22–26] tries to address convergence by making PnP and RED algorithms exact proximal algorithms. The idea is to replace the denoiser of RED algorithms by a gradient descent operator and the one of PnP algorithm by a proximal operator. Theoretical convergence then follows from known convergence results of proximal splitting algorithms. In [23, 24], it is thus proposed to plug an explicit gradient-step denoiser of the form $D = \text{Id} - \nabla g$, for a tractable and potentially nonconvex potential $g$ parameterized by a neural network. As shown in [24], such a constrained parametrization does not harm denoising performance. The gradient-step denoiser guarantees convergence of RED methods without sacrificing performance, but it does not cover convergence of PnP algorithms. An extension to PnP has been addressed in [25, 26]: following [27], when $g$ is trained with nonexpansive gradient, the gradient-step denoiser can be written as a proximal operator $D = \text{Id} - \nabla g = \text{Prox}_{\phi}$ of a nonconvex potential $\phi$. A PnP scheme with this proximal denoiser becomes again a genuine proximal splitting algorithm associated to an explicit nonconvex functional. Following existing convergence results of the PGD and DRS algorithms in the non-convex setting, [25] proves convergence of PnP-PGD and PnP-DRS with proximal denoiser (called ProxPnP-PGD and ProxPnP-DRS).

The main limitation of this approach is that the proximal denoiser $D = \text{Prox}_{\phi}$ does not give tractability of $\text{Prox}_{\tau_{\phi}}$ for $\tau \neq 1$. Therefore, to be a provable converging proximal splitting algorithm, the stepsize of the overall PnP algorithm has to be fixed to $\tau = 1$. For instance, when used with stepsize $\tau = 1$, both the PGD algorithm [25, Theorem 4.1]

$$x_{k+1} \in \text{Prox}_{\phi} \circ (\text{Id} - \lambda \nabla f)(x_k) \quad (2)$$

and the DRS algorithm [25, Theorem 4.3]

$$x_{k+1} \in \left( \frac{1}{2} \text{Rprox}_{\lambda f} \circ \text{Rprox}_{\phi_{\sigma}} + \frac{1}{2} \text{Id} \right) (x_k) \quad (3)$$

are proved to converge to a stationary point of (1) for a regularization parameter $\lambda$ satisfying $\lambda L_f < 1$. This is an issue when the input image has mild degradations, for instance for low noise levels $\nu$, because in this case relevant solutions are obtained with a dominant data-fidelity term in (1) through high values $\lambda > 1/L_f$. Note that convergence results for the DRS algorithm without such restriction on $\lambda$ has also been proposed in [25, Theorem 4.4] but at the cost of additional technical assumptions on the geometry of the denoiser which are difficult to verify in practice.

In this work, our objective is to design convergent PnP algorithms with a proximal denoiser, and with minimal restriction on the regularization parameter $\lambda$. Contrary to previous works on PnP convergence [16, 17, 19, 23, 24, 28], we not only wish to adapt the denoiser but also the original optimization scheme of interest. We show and make use of the fact that the regularization $\phi$ (from $D = \text{Id} - \nabla g = \text{Prox}_{\phi}$) is weakly convex, with a weak convexity constant that depends only on the Lipschitz constant of the gradient of the learned potential $g$. We first recall the main result from [26] that
studies a new proximal algorithm, reminiscent to PGD, which converges for a fixed stepsize $\tau = 1$ with relaxed restriction on $\lambda$. Second, we adapt the convergence study from [9] of the DRS algorithm in the nonconvex setting and prove convergence without restriction on $\lambda$. With both results, the constraint on $\lambda$ is replaced by a constraint on the Lipschitz constant of the gradient of the learned potential $g$. This new condition can be naturally satisfied using a convex relaxation of the gradient-step denoiser.

## Contributions and outline

This work is an extension from [26]. After detailing the $\alpha$PGD relaxation of the Proximal Gradient Descent algorithm introduced in [26], we propose a new proof of convergence of the DRS algorithm, such that when used with a proximal denoiser, the corresponding PnP schemes ProxPnP-DRS can converge for any regularization parameter $\lambda$.

In section 2, we start by introducing a variety of concepts and properties that will be used in our convergence analysis.

In section 3, we explain how the gradient-step denoiser $D = \text{Id} - \nabla g$ introduced in [24] can write as the proximal operator of a $M$-weakly convex function $\phi$. This denoiser is then called proximal denoiser. The weak convexity constant $M$ depends only on the Lipschitz constant of $\nabla g$ and can be controlled with a convex relaxation of the denoiser.

In section 4, we focus on the nonconvex convergence of ProxPnP-PGD [25] i.e. the PnP version of the standard PGD algorithm, with plugged proximal denoiser. By exploiting the weak convexity of $\phi$, the condition for convergence is $\lambda L_f + M < 2$.

In section 5, we detail the $\alpha$PGD, a relaxed version of the PGD algorithm<fnref n="1" /> introduced in [26]. Its convergence is shown in Theorem 2 for a smooth convex function $f$ and a weakly convex one $\phi$. Corollary 2 then applies this result for proving convergence of ProxPnP-$\alpha$PGD, its PnP version with proximal denoiser, under the condition $\lambda L_f M < 1$. Having a multiplication, instead of an addition, between the constants $\lambda L_f$ and $M$ allows to use any regularization parameter $\lambda$, provided we decrease $M$ sufficiently.

In section 6, we propose a new convergence proof for the PnP-DRS algorithm with proximal denoiser, without any restriction on the regularization parameter $\lambda$. We extend the original nonconvex convergence result from [9] and prove that, under a new constraint on the weak convexity constant $M$, the PnP version of the DRS algorithm with proximal denoiser also converges to a stationary point of an explicit functional.

In section 7, we realize experiments for both image deblurring and image super-resolution applications. We demonstrate that, by eliminating the constraint on the regularization parameter $\lambda$, our PnP-$\alpha$PGD and PnP-DRS algorithms close the performance gap with state-of-the-art plug-and-play algorithms.

<footnote n="1">¹There are two different notions of relaxation in this paper. One is for the relaxation of the proximal denoiser and the other for the relaxation of the optimization scheme.</footnote>
# 2 Definitions and useful properties

## 2.1 Definitions

In this section, we consider $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ proper and lower semicontinuous. We recall that $f$ is proper if $\text{dom}(f) \neq \emptyset$ and lower semicontinuous (lsc) if $\forall x \in \mathbb{R}^n \liminf_{y \to x} f(y) \geq f(x)$. For $M \in \mathbb{R}$, $f$ is called $M$-weakly convex if $x \to f(x) + \frac{M}{2} \|x\|^2$ is convex. Moreover, it is coercive if $\lim_{\|x\| \to +\infty} f(x) = +\infty$.

### Nonconvex subdifferential

When $f$ is nonconvex, the standard notion of subdifferential

$$ \partial f(x) = \{\omega \in \mathbb{R}^n, \forall y \in \mathbb{R}^n, f(y) - f(x) - \langle \omega, y - x \rangle \geq 0\} \quad (4) $$

may not be informative enough. Indeed, with Fermat's rule, in both the convex and nonconvex cases, Argmin $f = \text{zeros}(\partial f)$, but for minimizing a nonconvex function, there is usually no hope to target a global minimum. Following [8], we will use as notion of subdifferential the limiting subdifferential

$$ \partial^{\text{lim}} f(x) = \{\omega \in \mathbb{R}^n, \exists x_k \to x, f(x_k) \to f(x), \omega_k \to \omega, \omega_k \in \hat{\partial} f(x_k)\} \quad (5) $$

with $\hat{\partial} f$ the Fréchet subdifferential of $f$ defined as

$$ \hat{\partial} f(x) = \{\omega \in \mathbb{R}^n, \liminf_{y \to x} \frac{f(y) - f(x) - \langle \omega, y - x \rangle}{\|x - y\|} \geq 0\}. \quad (6) $$

These three notions of subdifferential verify for $x \in \text{dom} f$

$$ \partial f(x) \subset \hat{\partial} f(x) \subset \partial^{\text{lim}} f(x). \quad (7) $$

and coincide when $f$ is convex.

A necessary (but not sufficient) condition for $x \in \mathbb{R}^n$ to be a local minimizer of a nonconvex function $f$ is $0 \in \hat{\partial} f(x)$ and thus $0 \in \partial^{\text{lim}} f(x)$. Last but not least, we give two important properties of the limiting subdifferential.

**Proposition 1** (Sum rule [29, 8.8(c)]). If $F = f + g$ with $f$ of class $C^1$ and $g$ proper. Then for $x \in \text{dom}(g)$,

$$ \partial^{\text{lim}} F(x) = \nabla F(x) + \partial^{\text{lim}} g(x). \quad (8) $$

**Proposition 2** (Closedness of the limiting subdifferential). For a sequence $(x_k, \omega_k) \in \text{Graph}(\partial^{\text{lim}} f)$, if $(x_k, \omega_k) \to (x, \omega)$ and $f(x_k) \to f(x)$, then $(x, \omega) \in \text{Graph}(\partial^{\text{lim}} f)$, where $\text{Graph}(F)$ stands for the graph of a point-to-set mapping

$$ \text{Graph } F = \{(x, y) \in \mathbb{R}^n \times \mathbb{R}^m, y \in F(x)\}. \quad (9) $$
For additional details on the notion of limiting subdifferential, we refer to [29, Chapter 8]. In the rest of the paper, the subdifferential $\partial f$ corresponds to the limiting subdifferential $\partial^{\text{lim}}$ defined in (5).

### Kurdyka-Łojasiewicz property

The Kurdyka-Łojasiewicz (KL) property is a local property that characterizes the shape of the function $f$ around its critical points $\{x \in \text{int} \, \text{dom}(f), 0 \in \partial f(x)\}$. It is a tool widely used in nonconvex optimization [8, 30–33]. We use the definition from [30]: **Definition 1** (Kurdyka-Łojasiewicz (KL) property). A function $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ is said to have the Kurdyka-Łojasiewicz property at $x^* \in \text{dom}(f)$ if there exists $\eta \in (0, +\infty)$, a neighborhood $U$ of $x^*$ and a continuous concave function $\psi: [0, \eta) \to \mathbb{R}_+$ such that $\psi(0) = 0$, $\psi$ is $C^1$ on $(0, \eta)$, $\psi' > 0$ on $(0, \eta)$ and $\forall x \in U \cap [f(x^*) < f < f(x^*) + \eta]$, the Kurdyka-Łojasiewicz inequality holds:

$$ \psi'(f(x) - f(x^*)) \text{dist}(0, \partial f(x)) \ge 1. \qquad (10) $$

Proper lsc functions that satisfy the Kurdyka-Łojasiewicz inequality at each point of $\text{dom}(\partial f)$ are called KL functions.

**Remark 1.** (i) For $f$ proper lsc, the KL inequality always holds at non-critical points $x^* \in \text{dom} \, \partial f$ [30, Remark 4.(b)].

(ii) For $f$ smooth and $f(x^*) = 0$, the KL inequality (10) writes $||\nabla(\psi \circ f)(x)|| \ge 1$. The KL property can thus be interpreted as the fact that, up to a reparameterization, the function is locally sharp [8].

Following [8, Theorem 2.9], to prove the single-point convergence of iterative algorithms for minimizing the function $F = f + g$, with $f$ or $g$ nonconvex, the first requirement is to show that $F$ is KL. However, the KL condition is not stable by sum. For the purpose of the present paper, we need to choose a subclass of KL functions that, on the one hand, is large enough to encompass all our functions of interest and, on the other hand, has minimal stability properties so that inclusion to that set is easy to verify. In particular, as some of our functions will be parameterized by neural networks, stability by sum and composition are minimal requirements. In this analysis, the set of real analytic functions is large enough to include all our functions of interest.

**Definition 2** (Real analytic). A function $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ with open domain is said to be (real) analytic at $x \in \text{dom}(f)$ if $f$ may be represented by a convergent power series on a neighborhood of $x$. The function is said to be analytic if it is analytic at each point of its domain. We can extend the definition for $f: \mathbb{R}^n \to \mathbb{R}^m$ which is analytic at $x \in \mathbb{R}^n$ if for $f = (f_1, \dots, f_m)$, all $f_i: \mathbb{R}^n \to \mathbb{R}$ are analytic at $x$.

Typical real analytic functions include polynomials, exponential functions, the logarithm, trigonometric and power functions. Moreover, real analytic functions have the following stability properties

**Lemma 1** ([34]). (i) The sum, product, and composition of real analytic functions are real analytic.

(ii) A partial derivative of $f$ real analytic is real analytic.

(iii) Real Analytic Inverse Function Theorem: Let $f: \mathbb{R}^n \to \mathbb{R}^m$ be real analytic in a neighborhood of $x \in \mathbb{R}^n$ and suppose $J_f(x)$ non-singular, then $f^{-1}$ is defined and
is real analytic in a neighborhood of $f(x)$. In particular, if $f$ is real analytic and $\forall x \in \mathbb{R}^n$, $J_f(x)$ is non-singular, then $f^{-1}$ is real analytic on $\mathrm{Im}(f)$.

Finally, the important point here is that real analytic functions are KŁ:

**Lemma 2** ([35]). Let $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ be a proper real analytic function, then for any critical point $x^* \in \mathrm{dom}(f)$, there exists a neighborhood $U$ of $x^*$, an exponent $\theta \in [\frac{1}{2}, 1)$ and a constant $C$ such that the following Lojasiewicz inequality holds:

$$ \forall x \in U, \quad |f(x) - f(x^*)|^\theta \le C ||\nabla f(x)||. \qquad (11) $$

**Remark 2.** The Lojasiewicz inequality is a particular case of the Kurdyka-Lojasiewicz inequality from Definition 1 with desingularizing function $\psi(s) = cs^{1-\theta}$.

## 2.2 Useful inequalities

In this section, we give different inequalities verified by weakly convex and smooth functions. These results constitute basic blocks that will be used in our different proofs of convergence.

**Proposition 3** (Inequalities for weakly convex functions). For $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ proper lsc and M-weakly convex with $M \in \mathbb{R}$,

(i) $\forall x, y$ and $t \in [0, 1]$,

$$ f(tx + (1-t)y) \le tf(x) + (1-t)f(y) + \frac{M}{2}t(1-t)\|x-y\|^2; \qquad (12) $$

(ii) $\forall x, y$, we have $\forall z \in \partial f(y)$,

$$ f(x) \ge f(y) + \langle z, x - y \rangle - \frac{M}{2} \|x - y\|^2; \qquad (13) $$

(iii) **Three-points inequality.** For $z^+ \in \mathrm{Prox}_f(z)$, we have, $\forall x$

$$ f(x) + \frac{1}{2} \|x - z\|^2 \ge f(z^+) + \frac{1}{2} \|z^+ - z\|^2 + \frac{1-M}{2} \|x - z^+\|^2. \qquad (14) $$

*Proof.* (i) and (ii) follow from the fact that $\phi + \frac{M}{2} \|x\|^2$ is convex. We now prove (iii). Optimality conditions of the proximal operator $z^+ \in \mathrm{Prox}_\phi(z)$ gives

$$ z - z^+ \in \partial\phi(z^+). \qquad (15) $$

Hence, by (ii), we have $\forall x$,

$$ \phi(x) \ge \phi(z^+) + \langle z - z^+, x - z^+ \rangle - \frac{M}{2} \|x - z^+\|^2, \qquad (16) $$
and therefore,

$$
\begin{align}
\phi(x) + \frac{1}{2} \|x-z\|^2 &\ge \phi(z^+) + \frac{1}{2} \|x-z\|^2 + \langle z-z^+, x-z^+ \rangle - \frac{M}{2} \|x-z^+\|^2 \nonumber \\
&= \phi(z^+) + \frac{1}{2} \|x-z^+\|^2 + \frac{1}{2} \|z-z^+\|^2 - \frac{M}{2} \|x-z^+\|^2 \tag{17} \\
&= \phi(z^+) + \frac{1}{2} \|z-z^+\|^2 + \frac{1-M}{2} \|x-z^+\|^2. \nonumber
\end{align}
$$

<figure><img src="image_3.png" /></figure>

**Proposition 4** (Descent Lemma [36, Lemma 2.64]). For $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ proper lsc, if $f$ is $L_f$-smooth on an open and convex set $C \subseteq \text{int} \, \text{dom}(f)$, $\forall (x,y) \in C \times C$, we have

$$
|f(x) - f(y) - \langle \nabla f(y), x - y \rangle| \le \frac{L_f}{2} \|y - x\|^2. \quad (18)
$$

3 Relaxed Proximal Denoiser

This section introduces the denoiser used in our PnP algorithms. We first recall the definition of the Gradient Step Denoiser and show in Proposition 5 how it can be constrained to be a proximal denoiser. We finally introduce the relaxed proximal denoiser.

3.1 Gradient Step Denoiser

In this paper, we make use of the Gradient Step Denoiser introduced in [23, 24], that
writes as a gradient step over a differentiable potential g<sub>&sigma;</sub> parameterized by a neural
network:

$$
D_{\sigma} = \text{Id} - \nabla g_{\sigma}. \tag{19}
$$

In [24], g_σ is chosen to be of the form g_σ(x) = 1/2 ||x - N_σ(x)||^2 with N_σ parameterized
with a DRUNet architecture [7]. This denoiser can then be trained to denoise white
Gaussian noise v_σ of various standard deviations σ by minimizing the l_2 denoising
loss E[||D_σ(x + v_σ) - x)||^2]. It is shown that the Gradient Step Denoiser (19), despite
being constrained to be a conservative vector field (as in [6]), achieves state-of-the-art
denoising performance.

3.2 Proximal Denoiser

We first present a characterization of the Gradient Step denoiser as a proximal operator of some weakly convex potential $\phi$. The proof of this result relies on [27].

**Proposition 5 (Proximal Gradient Step denoiser).** Let $g_\sigma : \mathbb{R}^n \to \mathbb{R}$ a $C^{k+1}$ function with $k \ge 1$ and $\nabla g_\sigma L_{g_\sigma}$-Lipschitz with $L_{g_\sigma} < 1$. Let $D_\sigma := \text{Id} - \nabla g_\sigma = \nabla h_\sigma$. Then,

(i) there is $\phi_\sigma : \mathbb{R}^n \to \mathbb{R}^n \cup \{+\infty\}$, $\frac{L_{g_\sigma}}{L_{g_\sigma+1}}$-weakly convex, such that $\text{Prox}_{\phi_\sigma}$ is one-to-one and

$$
D_{\sigma} = \text{Prox}_{\phi_{\sigma}}. \tag{20}
$$
Moreover $D_\sigma$ is injective, $\mathrm{Im}(D_\sigma)$ is open and there is a constant $K \in \mathbb{R}$ such that $\phi_\sigma$ is defined on $\mathrm{Im}(D_\sigma)$ by $\forall x \in \mathrm{Im}(D_\sigma)$,

$$ \begin{aligned} \phi_\sigma(x) &= g_\sigma(D_\sigma^{-1}(x)) - \frac{1}{2} \|D_\sigma^{-1}(x) - x\|^2 + K \\ &= h_\sigma^*(x) - \frac{1}{2} \|x\|^2 + K. \end{aligned} \quad (21) $$

where $h_\sigma^*$ stands for the convex conjugate of $h_\sigma$.

(ii) $\forall x \in \mathbb{R}^n$, $\phi_\sigma(x) \ge g_\sigma(x) + K$ and for $x \in \mathrm{Fix}(D_\sigma)$, $\phi_\sigma(x) = g_\sigma(x) + K$.

(iii) $\phi_\sigma$ is $C^k$ on $\mathrm{Im}(D_\sigma)$ and $\forall x \in \mathrm{Im}(D_\sigma)$, $\nabla\phi_\sigma(x) = D_\sigma^{-1}(x) - x = \nabla g_\sigma(D_\sigma^{-1}(x))$. Moreover, $\nabla\phi_\sigma$ is $\frac{Lg_\sigma}{1-Lg_\sigma}$-Lipschitz on $\mathrm{Im}(D_\sigma)$

**Proof. (i)** This result is an extension of the Proposition 3.1 from [25]. Under the same assumptions as ours, it is shown, using the characterization of the proximity operator from [27], that for $\hat{\phi}_\sigma : \mathbb{R}^n \to \mathbb{R}$ defined as

$$ \hat{\phi}_{\sigma}(x) := \begin{cases} g_{\sigma}(D_{\sigma}^{-1}(x)) - \frac{1}{2} \|D_{\sigma}^{-1}(x) - x\|^2 & \text{if } x \in \mathrm{Im}(D_{\sigma}), \\ +\infty & \text{otherwise} \end{cases} \quad (22) $$

we have $D_\sigma = \mathrm{Prox}_{\hat{\phi}_\sigma}$.

As $\mathrm{dom}(\hat{\phi}_\sigma) = \mathrm{Im}(D_\sigma)$ is non necessarily convex, we cannot show that $\hat{\phi}_\sigma$ is weakly convex. We thus propose another choice of function $\phi_\sigma$ such that $D_\sigma = \mathrm{Prox}_{\phi_\sigma}$. As we suppose that $\nabla g_\sigma = \mathrm{Id} - D_\sigma$ is $L_{g_\sigma}$ Lipschitz, $D_\sigma$ is $L_{g_\sigma} + 1$ Lipschitz. By [27, Proposition 2], we have $\forall x \in \mathcal{X}$,

$$ D_{\sigma}(x) \in \mathrm{Prox}_{\phi_{\sigma}}(x) \quad (23) $$

with $\phi_\sigma : \mathbb{R}^n \to \mathbb{R}^n \cup \{+\infty\}$ which is $(1 - \frac{1}{L+1}) = \frac{L}{L+1}$ weakly convex. The weak convexity constant being smaller than 1, $x \to \phi_\sigma(x) + \frac{1}{2} \|x - y\|^2$ is convex and $\mathrm{Prox}_{\phi_\sigma}$ is one-to-one. According to [27, Theorem 4 (b)], there exists a constant $K \in \mathbb{R}$ such that for any $C \subseteq \mathrm{Im}(D_\sigma)$ polygonally connected $\phi_\sigma = \hat{\phi}_\sigma + K$ on $C$ i.e.

$$ \forall y \in C, \quad \phi_{\sigma}(D_{\sigma}(y)) = \langle y, D_{\sigma}(y) \rangle - \frac{1}{2} \|D_{\sigma}(y)\|^2 - h_{\sigma}(y) + K \quad (24) $$

$$ = g_{\sigma}(y) - \frac{1}{2} \|y - D_{\sigma}(y)\|^2 + K \quad (25) $$

Moreover, from the equality case of the Fenchel-Young inequality,

$$ \langle y, D_{\sigma}(y) \rangle = \langle y, \nabla h_{\sigma}(y) \rangle = h_{\sigma}(y) + h_{\sigma}^{*}(D_{\sigma}(y)). \quad (26) $$

Combining (24) and (26), we get

$$ \phi_{\sigma}(D_{\sigma}(y)) = h_{\sigma}^{*}(D_{\sigma}(y)) - \frac{1}{2} \|D_{\sigma}(y)\|^2 + K \quad (27) $$
As $\nabla g_{\sigma}$ is $\mathcal{C}^2$ and $L < 1$-Lipschitz, $\forall x \in \mathbb{R}^n$, $J_{D_{\sigma}}(x) = \text{Id} - \nabla^2 g_{\sigma}(x)$ is positive definite and $D_{\sigma}$ is injective. We can consider its inverse $D_{\sigma}^{-1}$ on $\text{Im}(D_{\sigma})$, and we get (21) from (25) Also, the inverse function theorem ensures that $\text{Im}(D_{\sigma})$ is open in $\mathbb{R}^n$. As $D_{\sigma}$ is continuous, $\text{Im}(D_{\sigma})$ is connected and open, thus polygonally connected. Therefore, (24) is true on the whole $\text{Im}(D_{\sigma})$.

(ii) We have

$$ 
\begin{aligned}
\phi_{\sigma}(x) &= \frac{1}{2} \|x - x\|^2 + \phi_{\sigma}(x) \\
&\geq \frac{1}{2} \|x - D_{\sigma}(x)\|^2 + \phi_{\sigma}(D_{\sigma}(x)) \\
&= g_{\sigma}(x) + K
\end{aligned} 
\quad (28) $$

where the first inequality comes from the definition of the proximal operator $D_{\sigma} = \text{Prox}_{\phi_{\sigma}}$ and the last equality is given by (21).

(iii) This is the application of [27, Corollary 6]. For the Lipschitz property, let $x, y \in \text{Im}(D_{\sigma})$, there exists $u, v \in \mathbb{R}^n$ such that $x = D_{\sigma}(u)$ and $y = D_{\sigma}(v)$. Hence, we have

$$ 
\begin{aligned}
\|\nabla\phi_{\sigma}(x) - \nabla\phi_{\sigma}(y)\| &= \|D_{\sigma}^{-1}(x) - D_{\sigma}^{-1}(y) - (x - y)\| \\
&= \|u - D_{\sigma}(u) - (v - D_{\sigma}(v))\| \\
&= \|\nabla g_{\sigma}(u) - \nabla g_{\sigma}(v)\| \\
&\leq L \|u - v\|
\end{aligned} 
\quad (29) $$

because $g_{\sigma}$ is $L$-Lipschitz. Moreover, as $J_{D_{\sigma}}(x) = \nabla^2 h_{\sigma}(x) = \text{Id} - \nabla^2 g_{\sigma}(x)$, $\forall u \in \mathbb{R}^n$, $\langle \nabla^2 h_{\sigma}(x)u, u \rangle = \|u\|^2 - \langle \nabla^2 g_{\sigma}(x)u, u \rangle \geq (1 - L) \|u\|^2$ and $h_{\sigma}$ is $(1 - L)$-strongly convex. Thus continuing from last inequalities

$$ 
\begin{aligned}
\|\nabla\phi_{\sigma}(x) - \nabla\phi_{\sigma}(y)\| &\leq \frac{L}{1-L} \|\nabla h_{\sigma}(u) - \nabla h_{\sigma}(v)\| \\
&= \frac{L}{1-L} \|D_{\sigma}(u) - D_{\sigma}(v)\| \\
&= \frac{L}{1-L} \|x - y\|.
\end{aligned} 
\quad (30) $$

□

This result states that a proximal denoiser can be defined from the denoiser (19), if the gradient of the learned potential $g_{\sigma}$ is contractive. In [25] the Lipschitz constant of $\nabla g_{\sigma}$ is softly constrained to satisfy $L_{g_{\sigma}} < 1$, by penalizing the spectral norm $\|\nabla^2 g_{\sigma}(x + \nu_{\sigma})\|_{S}$ in the denoiser training loss. The resulting proximal denoiser has a fixed weak-convexity constant that depends only on $L_{g_{\sigma}}$. In the next section, we propose a way to control this value without retraining $g_{\sigma}$.
### 3.3 More details on the regularization $\phi_{\sigma}$

We keep the parameterization of the potential $g_{\sigma}$ proposed by [25]:

$$g_{\sigma}(x) = \frac{1}{2} \|x - N_{\sigma}(x)\|^2, \qquad (31)$$

where $N_{\sigma}: \mathbb{R}^n \to \mathbb{R}^n$ is a $C^{\infty}$ neural network with softplus activations

$$s_{\epsilon}(x) = \frac{1}{\epsilon} \log(1 + e^{\epsilon x}) \qquad (32)$$

As explained in [24], $g_{\sigma}$ can be made coercive by choosing a convex compact set $C \subset \mathbb{R}^n$ where the iterates should stay and by adding an extra term

$$g_{\sigma}(x) = \frac{1}{2} \|x - N_{\sigma}(x)\|^2 + \frac{1}{2} \|x - \Pi_{C}(x)\|^2 \qquad (33)$$

with for example $C = [-1, 2]^n$ for images taking their values in $[0; 1]^n$. Assuming that $L_{g_{\sigma}} < 1$, the potential $\phi_{\sigma}$ obtained with Proposition 5 then verifies:

**Proposition 6 (Properties of $\phi_{\sigma}$).** With the above parameterization of $g_{\sigma}$, let $\phi_{\sigma}$ obtained from $g_{\sigma}$ via Proposition 5, then $\phi_{\sigma}$ verifies

(i) $\phi_{\sigma}$ is lower-bounded.

(ii) $\phi_{\sigma}$ is $\frac{L}{L+1}$-weakly convex.

(iii) $\phi_{\sigma}$ is $C^{\infty}$ and $\nabla\phi_{\sigma}$ is $\frac{L}{1-1}$-Lipschitz continuous on $\text{Im}(D_{\sigma})$.

(iv) $\phi_{\sigma}$ is coercive when $g_{\sigma}$ is coercive.

(v) $\phi_{\sigma}$ is real analytic on $\text{Im}(D_{\sigma})$.

*Proof.* (i) As $g_{\sigma}$ is non-negative, this follows directly from Proposition 5(ii).

(ii) Proposition 5(i)

(iii) Proposition 5(iii)

(iv) This follows directly from $\phi_{\sigma}(x) \ge g_{\sigma}(x)$ (Proposition 5(ii)).

(v) $g_{\sigma}$ is parameterized with the softplus activation, which is a real analytic function. By composition and sum of real analytic functions (Lemma 1(i)) $N_{\sigma}$ and then subsequently $g_{\sigma}$ are real analytic functions. Then, by Lemma 1(ii), $\nabla g_{\sigma}$ and thus $D_{\sigma}$ are also real analytic. By the inverse function theorem (Lemma 1(iii)), as $\forall x \in \mathbb{R}^n$, $J_{D_{\sigma}}(x) > 0$, $D_{\sigma}^{-1}$ is then real analytic on $\text{Im}(D_{\sigma})$. We finally obtain, using the expression of $\phi_{\sigma}$ on $\text{Im}(D_{\sigma})$ (21), again by sum and composition, that $\phi_{\sigma}$ is real analytic on $\text{Im}(D_{\sigma})$. $\square$

### 3.4 Relaxed Denoiser

Once trained, the Gradient Step Denoiser $D_{\sigma} = \text{Id} - \nabla g_{\sigma}$ can be relaxed as in [24] with a parameter $\gamma \in [0, 1]$

$$D_{\sigma}^{\gamma} = \gamma D_{\sigma} + (1 - \gamma) \text{Id} = \text{Id} - \gamma \nabla g_{\sigma}. \qquad (34)$$
Applying Proposition 5 with $g_{\sigma}^{\gamma} = \gamma g_{\sigma}$ which has a $\gamma L_{g_{\sigma}}$-Lipschitz gradient, we get that if $\gamma L_g < 1$, there exists a $\frac{\gamma L_{g_{\sigma}}}{\gamma L_{g_{\sigma}}+1}$-weakly convex $\phi_{\sigma}^{\gamma}$ such that

$$D_{\sigma}^{\gamma} = \text{Prox}_{\phi_{\sigma}^{\gamma}}, \qquad (35)$$

satisfying $\phi_{\sigma}^{0} = 0$ and $\phi_{\sigma}^{1} = \phi_{\sigma}$. Hence, one can control the weak convexity of the regularization function by relaxing the proximal denoising operator $D_{\sigma}^{\gamma}$.

# 4 PnP Proximal Gradient Descent (PnP-PGD)

In this section, we give convergence results for the ProxPnP-PGD algorithm

$$ \begin{aligned} x_{k+1} &= D_{\sigma} \circ (\text{Id} - \lambda f)(x_k) \\ &= \text{Prox}_{\phi_{\sigma}} \circ (\text{Id} - \lambda f)(x_k) \end{aligned} \qquad (36) $$

which is the PnP version of PGD, with plugged Proximal Denoiser (20). The convergence proof proposed in [25] for this algorithm is suboptimal as the weak convexity of $\phi_{\sigma}$ is not exploited. Doing so, we improve here the condition on the regularization parameter $\lambda$ for convergence.

## 4.1 Proximal Gradient Descent with a weakly convex function

We consider the general minimization problem

$$ \min_{x} F(x) := f(x) + \phi(x), \qquad (37) $$

for a smooth nonconvex function $f$ and a weakly convex function $\phi$ that are both bounded from below. We study under which conditions the classical Proximal Gradient Descent

$$ x_{k+1} \in \text{Prox}_{\tau\phi} \circ (\text{Id} - \tau \nabla f)(x_k) \qquad (38) $$

converges to a stationary point of (37). We first show convergence of the function values, and then convergence of the iterates, if $F$ verifies the Kurdyka-Łojasiewicz (KL) property (Definition 1).

**Theorem 1** (Convergence of the PGD algorithm (38)). Assume $f$ and $\phi$ proper lsc, bounded from below with $f$ differentiable with $L_f$-Lipschitz gradient, and $\phi$ M-weakly convex. Then for $\tau < \max(\frac{2}{L_f+M}, \frac{1}{L_f})$, the iterates (38) verify

(i) $(F(x_k))$ is non-increasing and converges.

(ii) The sequence has finite length i.e. $\sum_{k=0}^{+\infty} ||x_{k+1} - x_k||^2 < +\infty$ and $||x_{k+1} - x_k||$ converges to 0 at rate $\min_{k<K} ||x_{k+1} - x_k|| = O(1/\sqrt{K})$.

(iii) All cluster points of the sequence $x_k$ are stationary points of $F$.

(iv) If the sequence $(x_k)$ is bounded and if $F$ verifies the KL property at the cluster points of $(x_k)$, then $(x_k)$ converges, with finite length, to a stationary point of $F$.

**Remark 3.** (i) The boundedness of the iterates $(x_k)$ is verified as soon as the objective $F$ is coercive. Indeed, it ensures that $\{F(x) \le F(x_0)\}$ is bounded and, since $F(x_k)$ is non-increasing (Theorem 1(i)), that the iterates remain bounded.
(ii) We explained in Section 2 that, in practice, the KL property is verified as soon as F is real analytic.

Proof. The proof follows standard arguments of the convergence analysis of the PGD in the nonconvex setting [8, 31, 37].

(i) Relation (38) leads to $\frac{x_k - x_{k+1}}{\tau} - \nabla f(x_k) \in \partial\phi(x_{k+1})$, by definition of the proximal operator. As $\phi$ is M-weakly convex, we have from Proposition 3 (ii)

$$ \phi(x_k) \geq \phi(x_{k+1}) + \frac{\|x_k - x_{k+1}\|^2}{\tau} - \frac{M}{2} \|x_k - x_{k+1}\|^2 + \langle \nabla f(x_k), x_{k+1} - x_k \rangle. \quad (39) $$

The descent Lemma (Proposition 4) gives for $f$:

$$ f(x_{k+1}) \leq f(x_k) + \langle \nabla f(x_k), x_{k+1} - x_k \rangle + \frac{L_f}{2} \|x_k - x_{k+1}\|^2. \quad (40) $$

Combining both inequalities, for $F = f + \phi$, we obtain

$$ F(x_k) \geq F(x_{k+1}) + \left( \frac{1}{\tau} - \frac{M + \lambda L_f}{2} \right) \|x_k - x_{k+1}\|^2. \quad (41) $$

Therefore, if $\tau < 2/(M + \lambda L_f)$, $(F(x_k))$ is non-increasing. As $F$ is assumed lower-bounded, $(F(x_k))$ converges. We call $F^*$ its limit.

(ii) Summing (41) over $k = 0, 1, \dots, m$ gives

$$ \begin{aligned} \sum_{k=0}^{m} \|x_{k+1} - x_k\|^2 &\leq \frac{1}{\frac{1}{\tau} - \frac{M+L_f}{2}} (F(x_0) - F(x_{m+1})) \\ &\leq \frac{1}{\frac{1}{\tau} - \frac{M+L_f}{2}} (F(x_0) - F^*). \end{aligned} \quad (42) $$

Therefore, $\lim_{k \to \infty} \|x_{k+1} - x_k\| = 0$ with the convergence rate

$$ \gamma_k = \min_{0 \leq i \leq k} \|x_{i+1} - x_i\|^2 \leq \frac{1}{k} \frac{F(x_0) - F^*}{\frac{1}{\tau} - \frac{M+L_f}{2}}. \quad (43) $$

(iii) Suppose that a subsequence $(x_{k_i})_i$ is converging towards $x$. Let us show that $x$ is a critical point of $F$. By optimality of the proximal operator of $\phi$, for all $k \geq 0$

$$ \frac{x_{k_i+1} - x_{k_i}}{\tau} - \nabla f(x_{k_i+1}) \in \partial\phi(x_{k_i+1}), \quad (44) $$

where $\partial\phi$ stands for the limiting subdifferential (5). In other words, for $\omega_{k_i} = \frac{x_{k_i} - x_{k_i-1}}{\tau} - \nabla f(x_{k_i})$, $(\omega_{k_i}, x_{k_i}) \in \text{Graph } \partial\phi$. From the continuity of $\nabla f$, we have $\nabla f(x_{k_i}) \to \nabla f(x)$. As $\|x_{k_i+1} - x_{k_i}\| \to 0$, we get

$$ \frac{x_{k_i} - x_{k_i-1}}{\tau} - \nabla f(x_{k_i}) \to -\nabla f(x). \quad (45) $$
By closeness of $\partial\phi$ (Proposition 2), it remains to show that $\phi(x_{k_i}) \to \phi(x)$, to get $-\nabla f(x) \in \partial\phi(x)$, i.e. $x$ is a critical point of $F$ (by Proposition 1). To do so, we first use the fact that $g$ is lsc to obtain

$$ \liminf_{i \to \infty} \phi(x_{k_i}) \geq \phi(x). \quad (46) $$

On the other hand, by optimality of the Prox,

$$ \begin{aligned} & \phi(x_{k_i+1}) + \langle x_{k_i+1} - x_{k_i}, \nabla f(x_{k_i}) \rangle \frac{1}{2\tau} \|x_{k_i+1} - x_{k_i}\|^2 \\ & \leq \phi(x) + \langle x - x_{k_i}, \nabla f(x_{k_i}) \rangle \frac{1}{2\tau} \|x - x_{k_i}\|^2. \end{aligned} \quad (47) $$

Using that $x_{k_i} \to x$ and $x_{k_i+1} - x_{k_i} \to 0$ when $i \to +\infty$, we get

$$ \limsup_{i \to \infty} \phi(x_{k_i}) \leq \phi(x), \quad (48) $$

and

$$ \lim_{i \to \infty} \phi(x_{k_i}) = \phi(x). \quad (49) $$

(iv) We wish to apply Theorem 2.9 from [8]. We need to verify that the sequence $(x_k)$ satisfies the three conditions H1, H2, and H3 specified in [8, Section 2.3]:

* H1: Sufficient decrease condition $\forall k \in \mathbb{N}$,

$$ f(x_{k+1}) + a \|x_{k+1} - x_k\|^2 \leq F(x_k). \quad (50) $$

* H2: Relative error condition $\forall k \in \mathbb{N}$, there exists $\omega_{k+1} \in \partial F(x_{k+1})$ such that

$$ \|\omega_{k+1}\| \leq b \|x_{k+1} - x_k\|. \quad (51) $$

* H3: Continuity condition There exists a subsequence $(x_{k_i})_{i \in \mathbb{N}}$ and $\hat{x} \in \mathbb{R}^n$ such that

$$ x_{k_i} \to \hat{x} \text{ and } F(x_{k_i}) \to F(\hat{x}) \text{ as } i \to +\infty. \quad (52) $$

Condition H1 corresponds to the sufficient decrease condition shown in (i). For condition H3, we use that, as $(x_k)$ is bounded, there exists a subsequence $(x_{k_i})$ converging towards $y$. Then $F(x_{k_i}) \to F(y)$ has been shown in (iii). Finally, for condition H2, from (39), we had that

$$ z_{k+1} + \lambda \nabla f(x_{k+1}) = \frac{x_{k+1} - x_k}{\tau} \quad (53) $$

where $z_{k+1} \in \partial\phi(x_{k+1})$. Therefore,

$$ \|z_{k+1} + \lambda \nabla f(x_{k+1})\| = \frac{\|x_{k+1} - x_k\|}{\tau} \quad (54) $$

which gives condition H2. $\square$
## 4.2 ProxPnP Proximal Gradient Descent (ProxPnP-PGD)

Equipped with the convergence of PGD, we can now study the convergence of ProxPnP-PGD, the PnP-PGD algorithm with plugged Proximal Denoiser (20):

$$x_{k+1} = D_{\sigma} \circ (\mathrm{Id} - \lambda f)(x_k) = \mathrm{Prox}_{\phi_{\sigma}} \circ (\mathrm{Id} - \lambda f)(x_k). \quad (55)$$

This algorithm corresponds to the PGD algorithm (38) with fixed stepsize $\tau = 1$ and targets stationary points of the functional $F$ defined by:

$$F := \lambda f + \phi_{\sigma}. \quad (56)$$

At each iteration, $x_k \in \mathrm{Im}(D_\sigma)$, and thus $\phi_\sigma$ defined in (21) is tractable along the algorithm. The value of the objective function (56) at $x_k$ is

$$F(x_k) = \lambda f(x_k) + g_\sigma(z_k) - \frac{1}{2} \|z_k - x_k\|^2 + K \quad (57)$$

where $z_k = (\mathrm{Id} - \lambda \nabla f)(x_{k-1}) = D_\sigma^{-1}(x_k)$.

The following result, obtained from Theorem 1, improves [25] using the fact that the potential $\phi_\sigma$ is not any nonconvex function but a weakly convex one.

**Corollary 1 (Convergence of ProxPnP-PGD).** Let $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ be differentiable with $L_f$-Lipschitz gradient, bounded from below. Assume that $L_{g_\sigma}$ the Lipschitz constant of $\nabla g_\sigma$ verifies $L_{g_\sigma} < 1$. Then, for $\lambda L_f < \frac{L_{g_\sigma}+2}{L_{g_\sigma}+1}$, the iterates $(x_k)$ given by the iterative scheme (55) verify

(i) $(F(x_k))$ is non-increasing and converges.

(ii) The sequence has finite length i.e. $\sum_{k=0}^{+\infty} \|x_{k+1} - x_k\|^2 < +\infty$ and $\|x_{k+1} - x_k\| \le \mathcal{O}(1/\sqrt{K})$ converges to 0 at rate $\min_{k<K} \|x_{k+1} - x_k\| = \mathcal{O}(1/\sqrt{K})$.

(iii) If $f$ is real analytic, and $g_\sigma$ is coercive, then the iterates $(x_k)$ converge towards a critical point of $F$.

*Proof.* This is a direct application of Theorem 1 with $\lambda f$ the smooth function and with $\phi_\sigma$ the weakly convex function. The stepsize condition

$$\tau < \max\left(\frac{2}{\lambda L_f + \frac{L}{L+1}}, \frac{1}{\lambda L_f}\right) \quad (58)$$

becomes, with $\tau = 1$,

$$\Leftrightarrow \lambda L_f + \frac{L_{g_\sigma}}{L_{g_\sigma} + 1} < 2 \text{ or } \lambda L_f < 1 \quad (59)$$

$$\Leftrightarrow \lambda L_f < \frac{L_{g_\sigma} + 2}{L_{g_\sigma} + 1}. \quad (60)$$

For (iii), as shown in Section 3.3, $\phi_\sigma$ is real analytic on $\mathrm{Im}(D_\sigma)$, lower-bounded and coercive if $g_\sigma$ is coercive. As $f$ is assumed real analytic, by Lemma 1, $F = \lambda f + \phi_\sigma$ is
also real analytic and thus KL on Im($D_σ$). This is enough for the proof of Theorem 1 to hold because any limit point of ($x_k$) is a fixed point of the PnP-PGD operator $D_σ \circ (\text{Id} - \lambda f)$ and thus belongs to Im($D_σ$). □

By exploiting the weak convexity of $\phi_σ$, the convergence condition $\lambda L_f < 1$ of [25] is here replaced by $\lambda L_f < \frac{L_{g_σ}+2}{L_{g_σ}+1}$. Even if the bound is improved, we are still limited to regularization parameters satisfying $\lambda L_f < 2$.

In the same way, it was proven in [25, Theorem 4.3] that the PnP version of Douglas-Rachford Splitting (DRS) scheme with Proximal denoiser (ProxPnP-DRS), which writes

$$ \left\{ \begin{array}{l} y_{k+1} \in \text{Prox}_{\lambda f}(x_k) \end{array} \right. \qquad (61a) $$

$$ z_{k+1} = \text{Prox}_{\phi_σ}(2y_{k+1} - x_k) \qquad (61b) $$

$$ x_{k+1} = x_k + 2\beta(z_{k+1} - y_{k+1}), \qquad (61c) $$

converges towards a critical point of $\lambda f + \phi_σ$ provided $\lambda L_f < 1$ (see Section 6 for more details).

In both cases, the value of the regularization trade-off parameter $\lambda$ is then limited. This is an issue when restoring an image with mild degradations for which relevant solutions are obtained with a low amount of regularization and a dominant data-fidelity term in $F = \lambda f + g$. We could argue that this is a not a problem as the regularization strength is also regulated by the $\sigma$ parameter, which we are free to tune manually. However, it is observed in practice, for instance in [24], that the performance of PnP method greatly benefits from the ability to tune the two regularization parameters.

Given this limitation, our objective is to design new convergent PnP algorithms with proximal denoiser, and with minimal restriction on the regularization parameter $\lambda$. We not only wish to adapt the denoiser but also the original optimization schemes of interest.

# 5 PnP Relaxed Proximal Gradient Descent (PnP-αPGD)

In this section, we study a relaxation of the general Proximal Gradient Descent algorithm, called αPGD, such that when used with the proximal denoiser, the corresponding PnP scheme ProxPnP-αPGD can converge for a wider range of regularization parameters $\lambda$.

## 5.1 αPGD algorithm

We first introduce the αPGD algorithm for solving the general problem

$$ \min_{x \in \mathbb{R}^n} f(x) + \phi(x) \qquad (62) $$
for $f, \phi: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ with $f$ convex and smooth and $\phi$ weakly convex. The algorithm writes, for $\alpha \in (0, 1)$,

$$ q_{k+1} = (1 - \alpha)y_k + \alpha x_k \qquad (63a) $$

$$ x_{k+1} = \operatorname{Prox}_{\tau\phi}(x_k - \tau\nabla f(q_{k+1})) \qquad (63b) $$

$$ y_{k+1} = (1 - \alpha)y_k + \alpha x_{k+1}. \qquad (63c) $$

Algorithm (63) with $\alpha = 1$ exactly corresponds to the PGD algorithm (38). This scheme is reminiscent of [38] (taking $\alpha = \theta_k$ and $\tau = \frac{1}{\theta_k L_f}$ in Algorithm 1 of [38]), which generalizes Nesterov-like accelerated proximal gradient methods [37, 39]. Inspired by [40], this scheme was derived from the Primal-Dual algorithm [41] with Bregman proximal operator [42]. We describe this derivation in Appendix A.

### $\alpha$ PGD convergence

In the convex setting, using the convergence result from [42] of the Bregman Primal-Dual algorithm (from which $\alpha$PGD was derived), one can prove convergence of $\alpha$PGD for $\tau\lambda L_f < 2$ and small enough $\alpha$. However, in our case, the objective $F = \lambda f + \phi$ is nonconvex and a new convergence result needs to be derived. This is done with the following theorem.

**Theorem 2 (Convergence of $\alpha$PGD (63)).** Assume $f$ and $\phi$ proper, lsc, lower-bounded, with $f$ convex and $L_f$-smooth and $\phi$ M-weakly convex. Then for $\alpha \in (0, 1)$ and $\tau < \min(\frac{1}{\alpha L_f}, \frac{\alpha}{M})$, the updates (63) verify

(i) $F(y_k) + \frac{\alpha}{2\tau}(1 - \frac{1}{\alpha})^2 ||y_k - y_{k-1}||^2$ is non-increasing and converges.

(ii) the sequence $y_k$ has finite length i.e. $\sum_{k=0}^{+\infty} ||y_{k+1} - y_k||^2 < +\infty$ and $||y_{k+1} - y_k||$ converges to 0 at rate $\min_{k<K} ||y_{k+1} - y_k|| = O(1/\sqrt{K})$.

(iii) All cluster points of the sequence $y_k$ are stationary points of $F$. In particular, if $g$ is coercive, then $y_k$ has a subsequence converging towards a stationary point of $F$.

**Remark 4.** With this theorem, $\alpha$PGD is shown to verify convergence of the iterates and of the norm of the residual to 0. Note that we do not have here the analog of Theorem 1(iv) on the iterates convergence using the KL hypothesis. Indeed, as we detail at the end of the proof, the nonconvex convergence analysis with KL functions from [8] or [31] does not extend to our case.

**Remark 5.** As shown in Appendix B, a better bound on $\tau$ could be found, but with little numerical gain. Moreover, when $\alpha = 1$, algorithms $\alpha$PGD (63) and PGD (38) are equivalent, but we get a slightly worst bound in Theorem 2 than in Theorem 1 ($\tau < \min(\frac{1}{L_f}, \frac{1}{M}) \le \frac{2}{L_f+M}$). Nevertheless, when used with $\alpha < 1$, we next show that the relaxed algorithm is more relevant in the perspective of PnP with proximal denoiser.

**Remark 6.** The data-fidelity term $f$ is here assumed convex with Lipschitz gradient. This is verified by the classical $L^2$ data-fidelity term that we use in our experiments of Section 7. However, this excludes several data-fidelity terms such as the $L^1$ or the Kullback-Leiber distances.

*Proof.* The proof relies on Lemma 4 and Proposition 3. It follows the general strategy of the proofs in [38], and also requires the convexity of $f$.
(i) and (ii) : We can write (63b) as

$$
\begin{align}
x_{k+1} &\in \arg\min_{y \in \mathbb{R}^n} \phi(y) + \langle \nabla f(q_{k+1}), y - x_k \rangle + \frac{1}{2\tau} \|y - x_k\|^2 \nonumber \\
&\in \arg\min_{y \in \mathbb{R}^n} \phi(y) + f(q_{k+1}) + \langle \nabla f(q_{k+1}), y - q_{k+1} \rangle + \frac{1}{2\tau} \|y - x_k\|^2 \tag{64} \\
&\in \arg\min_{y \in \mathbb{R}^n} \Phi(y) + \frac{1}{2\tau} \|y - x_k\|^2 \nonumber
\end{align}
$$

with $\Phi(y) := \phi(y) + f(q_{k+1}) + \langle \nabla f(q_{k+1}), y - q_{k+1} \rangle$. As $\phi$ is $M$-weakly convex, so is $\Phi$. The three-points inequality of Proposition 3 (iii) applied to $\Phi$ thus gives $\forall y \in \mathbb{R}^n$,

$$
\Phi(y) + \frac{1}{2\tau} \|y - x_k\|^2 \geq \Phi(x_{k+1}) + \frac{1}{2\tau} \|x_{k+1} - x_k\|^2 + \left(\frac{1}{2\tau} - \frac{M}{2}\right) \|x_{k+1} - y\|^2 \quad (65)
$$

that is to say,

$$
\begin{align}
& \phi(y) + f(q_{k+1}) + \langle \nabla f(q_{k+1}), y - q_{k+1} \rangle + \frac{1}{2\tau} \|y - x_k\|^2 \ge \nonumber \\
& \qquad \phi(x_{k+1}) + f(q_{k+1}) + \langle \nabla f(q_{k+1}), x_{k+1} - q_{k+1} \rangle + \frac{1}{2\tau} \|x_{k+1} - x_k\|^2 + \left(\frac{1}{2\tau} - \frac{M}{2}\right) \|x_{k+1} - y\|^2. \tag{66}
\end{align}
$$

Using relation (63c), the descent Lemma (Proposition 4) as well as the convexity on $f$,
we obtain

$$
\begin{align*}
& f(q_{k+1}) + \langle \nabla f(q_{k+1}), x_{k+1} - q_{k+1} \rangle \\
= & f(q_{k+1}) + \left\langle \nabla f(q_{k+1}), \frac{1}{\alpha}y_{k+1} + \left(1 - \frac{1}{\alpha}\right)y_k - q_{k+1} \right\rangle \\
= & \frac{1}{\alpha}\left(f(q_{k+1}) + \langle \nabla f(q_{k+1}), y_{k+1} - q_{k+1} \rangle\right) + \left(1 - \frac{1}{\alpha}\right)\left(f(q_{k+1}) + \langle \nabla f(q_{k+1}), y_k - q_{k+1} \rangle\right) \\
\ge & \frac{1}{\alpha}\left(f(y_{k+1}) - \frac{L_f}{2}\|y_{k+1} - q_{k+1}\|^2\right) + \left(1 - \frac{1}{\alpha}\right)\left(f(q_{k+1}) + \langle \nabla f(q_{k+1}), y_k - q_{k+1} \rangle\right) \\
\ge & \frac{1}{\alpha}\left(f(y_{k+1}) - \frac{L_f}{2}\|y_{k+1} - q_{k+1}\|^2\right) + \left(1 - \frac{1}{\alpha}\right)f(y_k). \tag{67}
\end{align*}
$$

From relations (63a) and (63c), we have $y_{k+1} - q_{k+1} = \alpha(x_{k+1} - x_k)$. Combining this expression with (66) and (67), we have for all $y \in \mathbb{R}^n$

$$
\begin{align}
& \phi(y) + \left(\frac{1}{\alpha} - 1\right) f(y_k) + \frac{1}{2\tau} \|y - x_k\|^2 f(q_{k+1}) + \langle \nabla f(q_{k+1}), y - q_{k+1} \rangle \tag{68} \\
& \geq \phi(x_{k+1}) + \frac{1}{\alpha} f(y_{k+1}) + \left(\frac{1}{2\tau} - \frac{M}{2}\right) \|x_{k+1} - y\|^2 + \left(\frac{1}{2\tau} - \frac{\alpha L_f}{2}\right) \|x_{k+1} - x_k\|^2. \nonumber
\end{align}
$$
Using the convexity of $f$ we get for all $y \in \mathbb{R}^n$,

$$ \begin{aligned} & \phi(y) + f(y) + \left(\frac{1}{\alpha} - 1\right) f(y_k) + \frac{1}{2\tau} \|y - x_k\|^2 \\ & \geq \phi(x_{k+1}) + \frac{1}{\alpha} f(y_{k+1}) + \left(\frac{1}{2\tau} - \frac{\alpha L_f}{2}\right) \|x_{k+1} - x_k\|^2 + \left(\frac{1}{2\tau} - \frac{M}{2}\right) \|x_{k+1} - y\|^2, \end{aligned} \quad (69) $$

while weak convexity of $\phi$ with relation (63c) gives

$$ \phi(x_{k+1}) \geq \frac{1}{\alpha}\phi(y_{k+1}) + \left(1 - \frac{1}{\alpha}\right)\phi(y_k) - \frac{M}{2}(1 - \alpha)\|y_k - x_{k+1}\|^2. \quad (70) $$

Using (69), (70) and $F = f + \phi$, we obtain $\forall y \in \mathbb{R}^n$

$$ \begin{aligned} & \left(\frac{1}{\alpha} - 1\right) (F(y_k) - F(y)) + \frac{1}{2\tau} \|y - x_k\|^2 \geq \frac{1}{\alpha} (F(y_{k+1}) - F(y)) \\ & + \left(\frac{1}{2\tau} - \frac{\alpha L_f}{2}\right) \|x_{k+1} - x_k\|^2 + \left(\frac{1}{2\tau} - \frac{M}{2}\right) \|x_{k+1} - y\|^2 - \frac{M}{2}(1 - \alpha) \|y_k - x_{k+1}\|^2. \end{aligned} \quad (71) $$

For $y = y_k$, we get

$$ \begin{aligned} \frac{1}{\alpha} (F(y_k) - F(y_{k+1})) &\geq \left(\frac{1}{2\tau} - \frac{\alpha L_f}{2}\right) \|x_{k+1} - x_k\|^2 \\ &\quad - \frac{1}{2\tau} \|y_k - x_k\|^2 \left(\frac{1}{2\tau} - \frac{M(2-\alpha)}{2}\right) \|x_{k+1} - y_k\|^2. \end{aligned} \quad (72) $$

For $\alpha \in (0, 1)$, using that

$$ y_k - x_{k+1} = \frac{1}{\alpha}(y_k - y_{k+1}) \quad (73) $$

$$ y_k - x_k = \left(1 - \frac{1}{\alpha}\right) (y_k - y_{k-1}) \quad (74) $$

we get

$$ \begin{aligned} F(y_k) - F(y_{k+1}) &\geq -\frac{\alpha}{2\tau} \left(1 - \frac{1}{\alpha}\right)^2 \|y_k - y_{k-1}\|^2 \\ &\quad + \alpha \left(\frac{1}{2\tau} - \frac{\alpha L_f}{2}\right) \|x_{k+1} - x_k\|^2 \\ &\quad + \frac{1}{\alpha} \left(\frac{1}{2\tau} - \frac{M(2-\alpha)}{2}\right) \|y_{k+1} - y_k\|^2. \end{aligned} \quad (75) $$
With the assumption $\tau\alpha L_f < 1$, the second term of the right-hand side is non-negative and therefore,

$$ 
\begin{aligned}
F(y_k) - F(y_{k+1}) &\geq -\frac{\alpha}{2\tau} \left(1 - \frac{1}{\alpha}\right)^2 \|y_k - y_{k-1}\|^2 \\
&\quad + \frac{1}{\alpha} \left(\frac{1}{2\tau} - \frac{M(2-\alpha)}{2}\right) \|y_{k+1} - y_k\|^2. \qquad (76) \\
&= -\delta \|y_k - y_{k-1}\|^2 + \delta \|y_{k+1} - y_k\|^2 + (\gamma - \delta) \|y_{k+1} - y_k\|^2
\end{aligned}
$$

with

$$ \delta = \frac{\alpha}{2\tau} \left(1 - \frac{1}{\alpha}\right)^2 \quad \text{and} \quad \gamma = \frac{1}{\alpha} \left(\frac{1}{2\tau} - \frac{M(2-\alpha)}{2}\right). \qquad (77) $$

We now make use of the following lemma.

**Lemma 3** ([43]). Let $(a_n)_{n \in \mathbb{N}}$ and $(b_n)_{n \in \mathbb{N}}$ be two real sequences such that $b_n \geq 0$ $\forall n \in \mathbb{N}$, $(a_n)$ is bounded from below and $a_{n+1} + b_n \leq a_n$ $\forall n \in \mathbb{N}$. Then $(a_n)_{n \in \mathbb{N}}$ is a monotonically non-increasing and convergent sequence and $\sum_{n \in \mathbb{N}} b_n < +\infty$.

To apply Lemma 3 we look for a stepsize satisfying

$$ \gamma - \delta > 0 \quad \text{i.e.} \quad \tau < \frac{\alpha}{M}. \qquad (78) $$

Therefore, hypothesis $\tau < \min\left(\frac{1}{\alpha L_f}, \frac{\alpha}{M}\right)$ gives that $(F(y_k) + \delta \|y_k - y_{k-1}\|^2)$ is a non-increasing and convergent sequence and that $\sum_k \|y_k - y_{k+1}\|^2 < +\infty$.

(iii) The proof of this result is an extension of the proof of Theorem 1(ii) in the context of the classical PGD. Suppose that a subsequence $(y_{k_i})$ is converging towards $y$. Let us show that $y$ is a critical point of $F$. From (63b), we have

$$ \frac{x_{k+1} - x_k}{\tau} - \nabla f(q_{k+1}) \in \partial\phi(x_{k+1}). \qquad (79) $$

First we show that $x_{k_i+1} - x_{k_i} \to 0$. We have $\forall k > 1$,

$$ 
\begin{aligned}
\|x_{k+1} - x_k\| &= \left\| \frac{1}{\alpha}y_{k+1} + \left(1 - \frac{1}{\alpha}\right)y_k - \frac{1}{\alpha}y_k - \left(1 - \frac{1}{\alpha}\right)y_{k-1} \right\| \\
&\leq \frac{1}{\alpha} \|y_{k+1} - y_k\| + \left(\frac{1}{\alpha} - 1\right) \|y_k - y_{k-1}\| \\
&\to 0.
\end{aligned}
\qquad (80)
$$
From (63a), we also get $||q_{k+1} - q_k|| \to 0$. Now, let us show that $x_{k_i} \to y$ and $q_{k_i} \to y$. First using (63c), we have

$$ \begin{aligned} ||x_{k_i} - y|| &\le ||x_{k_i+1} - y|| + ||x_{k_i+1} - x_{k_i}|| \\ &\le \frac{1}{\alpha} ||y_{k_i+1} - y|| + (\frac{1}{\alpha} - 1) ||y_{k_i} - y|| + ||x_{k_i+1} - x_{k_i}|| \quad (81) \\ &\to 0. \end{aligned} $$

Second, from (63a), we get in the same way $q_{k_i} \to y$. From the continuity of $\nabla f$, we obtain $\nabla f(q_{k_i}) \to \nabla f(y)$ and therefore

$$ \frac{x_{k_i} - x_{k_i-1}}{\tau} - \nabla f(q_{k_i}) \to -\nabla f(y). \quad (82) $$

If we can also show that $\phi(x_{k_i}) \to \phi(y)$, we get from the closeness of the limiting subdifferential (Proposition 2) that $-\nabla f(y) \in \partial\phi(y)$ i.e. $y$ is a critical point of $F$.

Using the fact that $\phi$ is lsc and $x_{k_i} \to y$, we have

$$ \liminf_{i \to \infty} \phi(x_{k_i}) \ge \phi(y). \quad (83) $$

On the other hand, with Equation (69) for $k+1 = k_i$, taking $i \to +\infty$, $||y - x_{k_i+1}|| \to 0$, $||y - x_{k_i}|| \to 0$, $f(y_{k_i}) \to f(y)$, $f(y_{k_i+1}) \to f(y)$ and we get

$$ \limsup_{i \to \infty} \phi(x_{k_i}) \le \phi(y), \quad (84) $$

and therefore

$$ \lim_{i \to \infty} \phi(x_{k_i}) = \phi(y). \quad (85) $$

As $f$ is lower-bounded, if $\phi$ is coercive, so is $F$ and by (i) the iterates $(y_k)$ remain bounded and $(y_k)$ admits a converging subsequence.

### Remark on the convergence of the iterates with the KL hypothesis

In order to prove a result similar to Theorem 1 (iv) on the convergence of the iterates with the KL hypothesis, we cannot directly apply Theorem 2.9 from [8] on $F$ as the objective function $F(x_k)$ by itself does not decrease along the sequence but $F(x_k) + \delta ||x_{k+1} - x_k||^2$ does (where $\delta = \frac{\alpha}{2\tau} (1 - \frac{1}{\alpha})^2$).

Our situation is more similar to the variant of this result presented in [31, Theorem 3.7]. Indeed, denoting $\mathcal{F}: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}$ defined as $\mathcal{F}(x, y) = F(x) + \delta ||x - y||^2$ and considering $\forall k \ge 1$, the sequence $z_k = (y_k, y_{k-1})$ with $y_k$ following our algorithm, we can easily show that $z_k$ verifies the conditions H1 and H3 specified in [31, Section 3.2]. However, condition H2 does not extend to our algorithm. $\square$

## 5.2 ProxPnP-αPGD algorithm

The αPGD algorithm (63) gives birth to the PnP-αPGD algorithm by replacing the proximity operator $\text{Prox}_{\tau g}$ by a Gaussian denoiser $D_\sigma$. Now, similar to what was done
in Section 4.2 with the PGD algorithm, in this section, we study the convergence of this PnP-αPGD with our particular Proximal Gradient-Step Denoiser (20) $D_σ = \text{Prox}_{\phi_σ}$ and stepsize $τ = 1$. This corresponds to the following algorithm, which we refer to ProxPnP-αPGD.

$$ 
\left\{ 
\begin{aligned}
q_{k+1} &= (1 - \alpha)y_k + \alpha x_k && (86a) \\
x_{k+1} &= D_σ(x_k - \lambda \nabla f(q_{k+1})) && (86b) \\
y_{k+1} &= (1 - \alpha)y_k + \alpha x_{k+1} && (86c)
\end{aligned}
\right. 
$$

Like ProxPnP-PGD (55), the ProxPnP-αPGD scheme targets the critical points of the explicit functional $F = \lambda f + \phi_σ$ where $\phi_σ$ is obtained from the deep potential $g_σ$ via Proposition 5. Applying the previous Theorem 2 on the convergence of the αPGD algorithm with $τ = 1$ and $g = \phi_σ$ which is $M = \frac{L_{g_σ}}{L_{g_σ}+1}$ weakly convex, we get the following convergence result for ProxPnP-αPGD.

**Corollary 2** (Convergence of ProxPnP-αPGD (86)). Let $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ convex and differentiable with $L_f$-Lipschitz gradient, bounded from below. Assume that $L_{g_σ}$ the Lipschitz constant of $\nabla g_σ$ verifies $L_{g_σ} < 1$. Let $M = \frac{L_{g_σ}}{L_{g_σ}+1}$ be the weak convexity constant of $\phi_σ$ obtained from $g_σ$ via Proposition 5. Then, for any $\alpha \in (0, 1)$ such that

$$ \lambda L_f M < 1 \text{ and } M < \alpha < 1/(\lambda L_f) \quad (87) $$

the iterates given by the iterative scheme (86) verify

(i) $F(y_k) + \frac{\alpha}{2}(1 - \frac{1}{\alpha})^2 ||y_k - y_{k-1}||^2$ is non-increasing and converges.

(ii) The sequence $y_k$ as finite length i.e. $\sum_{k=0}^{+\infty} ||y_{k+1} - y_k||^2 < +\infty$ and $||y_{k+1} - y_k||$ converges to 0 at rate $\min_{k<K} ||y_{k+1} - y_k|| = O(1/\sqrt{K})$.

(iii) If $g_σ$ is coercive, then $(y_k)$ has a subsequence converging towards a critical point of $F$.

The existence of $\alpha \in (0, 1)$ satisfying relation (87) is ensured as soon as $\lambda L_f M < 1$. As a consequence, when $M$ gets small (i.e $\phi_σ$ gets "more convex") $\lambda L_f$ can get arbitrarily large. This is a major advance compared to ProxPnP-PGD, which was limited to $\lambda L_f < 2$ in Corollary 1, even for convex $\phi$ ($M = 0$). To further exploit this property, we now consider the relaxed denoiser $D_σ^γ$ (34) that is associated to a function $\phi_σ^γ$ with a tunable weak convexity constant $M^γ$.

**Corollary 3** (Convergence of ProxPnP-αPGD with relaxed denoiser). Let $F := \lambda f + \phi_σ^γ$, with the $M^γ = \frac{\gamma L_{g_σ}}{\gamma L_{g_σ}+1}$-weakly convex potential $\phi_σ^γ$ introduced in (35) and $L_{g_σ} < 1$. Then, for $M^γ < \alpha < 1/(\lambda L_f)$, the iterates $x_k$ given by the ProxPnP-αPGD (86) with $\gamma$-relaxed denoiser $D_σ^γ$ defined in (34) verify the convergence properties (i)-(iii) of Theorem 2.

Therefore, using the $\gamma$-relaxed denoiser $D_σ^γ = \gamma D_σ + (1 - \gamma)$ Id, the overall convergence condition on $\lambda$ is now

$$ \lambda < \frac{1}{L_f} \left( 1 + \frac{1}{\gamma L} \right). \quad (88) $$

Provided $\gamma$ gets small, $\lambda$ can be arbitrarily large. Small $\gamma$ means small amount of regularization brought by denoising at each step of the PnP algorithm. Moreover, for
small $\gamma$, the targeted regularization function $\phi_{\sigma}^{\gamma}$ gets close to a convex function. It has already been observed that deep convex regularization can be suboptimal compared to more flexible nonconvex ones [23]. Depending on the inverse problem, and on the necessary amount of regularization, the choice of the couple $(\gamma, \lambda)$ will be of paramount importance for efficient restoration.

# 6 PnP Douglas-Rachford Splitting (PnP-DRS)

We now focus on the ProxPnP version of the Douglas-Rachford Splitting (DRS). It corresponds to a classical DRS algorithm, with fixed stepsize $\tau = 1$, for optimizing the function $\lambda f + \phi_{\sigma}$. The algorithm writes, for $\beta \in (0, 1]$,

$$y_{k+1} = \text{Prox}_{\phi_{\sigma}}(x_k) \qquad (89a)$$

$$z_{k+1} = \text{Prox}_{\lambda f}(2y_{k+1} - x_k) \qquad (89b)$$

$$x_{k+1} = x_k + 2\beta(z_{k+1} - y_{k+1}) \qquad (89c)$$

or equivalently

$$x_{k+1} \in (\beta \text{Rprox}_{\lambda f} \circ \text{Rprox}_{\phi_{\sigma}} \circ + (1 - \beta) \text{Id})(x_k). \qquad (90)$$

The case $\beta = \frac{1}{2}$ corresponds to the standard version of Douglas-Rachford Splitting, and is equivalent to ADMM, while $\beta = 1$ is generally referred to as Peaceman-Rachford splitting.

The convergence of this algorithm in the nonconvex setting, studied in [44] and [9], requires one of the two functions to be differentiable with globally Lipschitz gradient. Assuming that the data-fidelity term $f$ satisfies this property, we get convergence of the DRS algorithm provided that $\lambda L_f < 1$ [25, Theorem 4.3]. In the experiments, this version of DRS will be referred to *DRSdiff*.

Convergence of this algorithm without the restriction on the parameter $\lambda$ is possible if $\phi_{\sigma}$ has a globally Lipschitz gradient. However, $\phi_{\sigma}$ obtained from Proposition 5 is differentiable only on $\text{Im}(D_{\sigma}) = \text{Im}(\text{Prox}_{\phi_{\sigma}})$. In [25], convergence of the iterates was proved assuming $\text{Im}(D_{\sigma})$ to be convex in order to keep the descent lemma for $\phi_{\sigma}$ valid on $\text{Im}(D_{\sigma})$ (Proposition 4). However, the convexity of $\text{Im}(D_{\sigma})$ is difficult to verify and unlikely to be true in practice.

Instead, we propose a new adaptation of the original convergence results from [9, 44] on the convergence of the DRS algorithm with nonconvex functions. We directly write the convergence result in the ProxPnP case of interest. In particular, we consider a very general data-fidelity term $f$ nonconvex and non-differentiable and $\phi_{\sigma}$ obtained from Proposition 5 via the proximal Gradient step denoiser $D_{\sigma} = \text{Id} - \nabla g_{\sigma} = \text{Prox}_{\phi_{\sigma}}$. As detailed in Section 3.3, $\phi_{\sigma}$ is $\frac{L_{g_{\sigma}}}{1+L_{g_{\sigma}}}$-weakly convex on $\mathbb{R}^n$ and has $\frac{L_{g_{\sigma}}}{1-L_{g_{\sigma}}}$ Lipschitz gradient on $\text{Im}(D_{\sigma})$. Our convergence result makes use of these two properties. The main adaptation concerns the proof of the sufficient decrease
of the Douglas-Rachford envelope

$$F^{DR}(x, y, z) = \lambda f(z) + \phi_{\sigma}(y) + \frac{1}{\tau} \langle y - x, y - z \rangle + \frac{1}{2} \|y - z\|^2. \quad (91)$$

**Theorem 3** (Convergence of ProxPnP-DRS). Let $f$ proper, $lsc$, bounded from below. Assume that $L_{g\sigma}$ the Lipschitz constant of $\nabla g\sigma$ and $\beta \in (0, 1)$ verify

$$\beta(2L_{g\sigma}^3 - 3L_{g\sigma}^2 + 1) + (2L_{g\sigma}^2 + L_{g\sigma} - 1) < 0. \quad (92)$$

Then, for all $\lambda > 0$, the iterates $(y_k, z_k, x_k)$ given by the iterative scheme (89) verify (i) $(F^{DR}(x_{k-1}, y_k, z_k))$ is non-increasing and converges.

(ii) $x_k - x_{k-1} = \beta(y_k - z_k)$ vanishes with rate $\min_{k \le K} \|y_k - z_k\| = O(\frac{1}{\sqrt{K}})$.

(iii) If $f$ is real analytic, and $g\sigma$ is coercive, then the iterates $(x_k)$ converge towards a critical point of $F$.

**Remark 7.** The standard Douglas-Rachford algorithm is for $\beta = 1/2$. The restriction on $L_{g\sigma}$ is then

$$2L_{g\sigma}^3 + L_{g\sigma}^2 + 2L_{g\sigma} - 1 < 0 \quad (93)$$

which is verified for approximately $L_{g\sigma} < 0.376$. Using a smaller $\beta$ relaxes the constraint on $L_{g\sigma}$ until $L_{g\sigma} < 0.5$ when $\beta \to 0$. For example for $\beta = 0.25$, the constraint is verified up to approximately $L_{g\sigma} < 0.45$.

**Remark 8.** In addition to deal with non-differentiable data-fidelity terms, the main advantage of this result, compared to [25, Theorem 4.3], is that the condition on $\lambda$ disappears, and the convergence remains true for all regularization parameter $\lambda$. However, the restriction on $L_{g\sigma}$ the Lipschitz constant of $\nabla g\sigma$ is stronger than just $L_{g\sigma} < 1$, which should harm the denoising capacity of the denoiser.

*Proof. (i)* As explained above, this point represents the main adaptation. Following the initial proof from [9, Theorem 1], we derive a new sufficient decrease property for $F^{DR}$. We first rewrite the envelope $F^{DR}$ from (91) as

$$\begin{aligned} F^{DR}(x, y, z) &= \lambda f(z) + \phi_{\sigma}(y) + \langle y - x, y - z \rangle + \frac{1}{2} \|y - z\|^2 \\ &= \lambda f(z) + \phi_{\sigma}(y) + \frac{1}{2} \|(2y - x) - z\|^2 - \frac{1}{2} \|x - y\|^2. \end{aligned} \quad (94)$$

As $z_{k+1} = \text{Prox}_f(2y_{k+1} - x_k)$ (89b), denoting $\forall k > 0$ $F_k^{DR} = F^{DR}(x_{k-1}, y_k, z_k)$,

$$\begin{aligned} F_{k+1}^{DR} &= \min_{z \in \mathbb{R}^n} \lambda f(z) + \phi_{\sigma}(y_{k+1}) + \frac{1}{2} \|(2y_{k+1} - x_k) - z\|^2 - \frac{1}{2} \|x_k - y_{k+1}\|^2 \\ &= \min_{z \in \mathbb{R}^n} \lambda f(z) + \phi_{\sigma}(y_{k+1}) + \langle y_{k+1} - x_k, y_{k+1} - z \rangle + \frac{1}{2} \|y_{k+1} - z\|^2. \end{aligned} \quad (95)$$

The optimality condition for (89a) is

$$x_k - y_{k+1} \in \partial\phi_{\sigma}(y_{k+1}). \quad (96)$$
As $\phi_{\sigma}$ is differentiable on $\mathrm{Im}(D_{\sigma})$ and $y_{k+1} \in \mathrm{Im}(D_{\sigma})$,

$$x_k - y_{k+1} = \nabla \phi_{\sigma}(y_{k+1}). \qquad (97)$$

The Douglas-Rachford envelope $F^{DR}$ then writes

$$F_{k+1}^{DR} = \min_{z \in \mathbb{R}^n} \lambda f(z) + \phi_{\sigma}(y_{k+1}) + \langle \nabla \phi_{\sigma}(y_{k+1}), z - y_{k+1} \rangle + \frac{1}{2} \|y_{k+1} - z\|^2. \quad (98)$$

We get, by comparing with the right term when $z = z_k$,

$$\begin{aligned} F_{k+1}^{DR} & \leq \lambda f(z_k) + \phi_{\sigma}(y_{k+1}) + \langle \nabla \phi_{\sigma}(y_{k+1}), z_k - y_{k+1} \rangle + \frac{1}{2} \|y_{k+1} - z_k\|^2 \\ & = \phi_{\sigma}(y_{k+1}) + \langle \nabla \phi_{\sigma}(y_{k+1}), y_k - y_{k+1} \rangle + \langle \nabla \phi_{\sigma}(y_{k+1}), z_k - y_k \rangle \\ & \quad + \lambda f(z_k) + \frac{1}{2} \|y_{k+1} - z_k\|^2. \end{aligned} \qquad (99)$$

From Proposition 5, $\phi_{\sigma}$ is $M = \frac{L_{g\sigma}}{L_{g\sigma}+1}$ weakly convex. Then with Proposition 3 (ii) we have

$$\phi_{\sigma}(y_{k+1}) + \langle \nabla \phi_{\sigma}(y_{k+1}), y_k - y_{k+1} \rangle \leq \phi_{\sigma}(y_k) + \frac{M}{2} \|y_k - y_{k+1}\|^2 \quad (100)$$

and

$$\begin{aligned} F_{k+1}^{DR} & \leq \phi_{\sigma}(y_k) + \lambda f(z_k) + \langle \nabla \phi_{\sigma}(y_{k+1}), z_k - y_k \rangle \\ & \quad + \frac{M}{2} \|y_k - y_{k+1}\|^2 + \frac{1}{2} \|y_{k+1} - z_k\|^2. \end{aligned} \qquad (101)$$

Moreover,

$$\langle \nabla \phi_{\sigma}(y_{k+1}), z_k - y_k \rangle = \langle \nabla \phi_{\sigma}(y_k), z_k - y_k \rangle + \langle \nabla \phi_{\sigma}(y_{k+1}) - \nabla \phi_{\sigma}(y_k), z_k - y_k \rangle \quad (102)$$

and using

$$\begin{aligned} F_k^{DR} &= \lambda f(z_k) + \phi_{\sigma}(y_k) + \langle y_k - x_{k-1}, y_k - z_k \rangle + \frac{1}{2} \|y_k - z_k\|^2 \\ &= \lambda f(z_k) + \phi_{\sigma}(y_k) + \langle \nabla \phi_{\sigma}(y_k), z_k - y_k \rangle + \frac{1}{2} \|y_k - z_k\|^2 \end{aligned} \qquad (103)$$

we get

$$\begin{aligned} F_{k+1}^{DR} & \leq F_k^{DR} + \langle \nabla \phi_{\sigma}(y_{k+1}) - \nabla \phi_{\sigma}(y_k), z_k - y_k \rangle \\ & \quad - \frac{1}{2} \|y_k - z_k\|^2 + \frac{1}{2} \|y_{k+1} - z_k\|^2 + \frac{M}{2} \|y_k - y_{k+1}\|^2. \end{aligned} \qquad (104)$$
Now using

$$ \frac{1}{2} \|y_{k+1} - z_k\|^2 = \frac{1}{2} \|y_k - y_{k+1}\|^2 + \frac{1}{2} \|y_k - z_k\|^2 + \langle y_{k+1} - y_k, y_k - z_k \rangle \quad (105) $$

we get

$$ F_{k+1}^{DR} \le F_k^{DR} + \langle (\nabla \phi_{\sigma}(y_{k+1}) - \nabla \phi_{\sigma}(y_k)) - (y_{k+1} - y_k), z_k - y_k \rangle + \frac{1+M}{2} \|y_k - y_{k+1}\|^2. \quad (106) $$

From (89c) and (97)

$$ z_k - y_k = \frac{1}{2\beta}(x_k - x_{k-1}) = \frac{1}{2\beta}(y_{k+1} - y_k) + \frac{1}{2\beta}(\nabla\phi_{\sigma}(y_{k+1}) - \nabla\phi_{\sigma}(y_k)), \quad (107) $$

so that

$$ F_{k+1}^{DR} - F_k^{DR} \le \left(1 + M - \frac{1}{\beta}\right) \frac{1}{2} \|y_k - y_{k+1}\|^2 + \frac{1}{2\beta} \|\nabla\phi_{\sigma}(y_{k+1}) - \nabla\phi_{\sigma}(y_k)\|^2. \quad (108) $$

Using the fact that $\nabla\phi_{\sigma}$ is $\frac{L_{g_{\sigma}}}{1-L_{g_{\sigma}}}$ Lipschitz on $\mathrm{Im}(D_{\sigma})$ and that $\forall k > 0, y_k \in \mathrm{Im}(D_{\sigma})$, it simplifies to

$$ F_{k+1}^{DR} - F_k^{DR} \le \frac{1}{2} \left( 1 + M + \frac{1}{\beta} \left( \frac{L_{g_{\sigma}}^2}{(1 - L_{g_{\sigma}})^2} - 1 \right) \right) \|y_k - y_{k+1}\|^2. \quad (109) $$

Recalling that $M = \frac{L_{g_{\sigma}}}{1+L_{g_{\sigma}}}$, the condition on the stepsize for $F_{k+1}^{DR} \le F_k^{DR}$ is thus

$$ \begin{aligned} & 1 + M + \frac{1}{\beta} \left( \frac{L_{g_{\sigma}}^2}{(1 - L_{g_{\sigma}})^2} - 1 \right) < 0 \\ \Leftrightarrow & \frac{2L_{g_{\sigma}} + 1}{L_{g_{\sigma}} + 1} + \frac{1}{\beta} \frac{2L_{g_{\sigma}} - 1}{(1 - L_{g_{\sigma}})^2} < 0 \\ \Leftrightarrow & \beta(2L_{g_{\sigma}}^3 - 3L_{g_{\sigma}}^2 + 1) + (2L_{g_{\sigma}}^2 + L_{g_{\sigma}} - 1) < 0. \end{aligned} \quad (110) $$

For completeness, recalling that $D_{\sigma} = \nabla h_{\sigma}$ with $h_{\sigma} = 1 - L_{g_{\sigma}}$ strongly convex, we have

$$ \|y_k - y_{k+1}\| = \|D_{\sigma}(x_k) - D_{\sigma}(x_{k-1})\| \ge (1 - L_{g_{\sigma}}) \|x_k - x_{k-1}\| \quad (111) $$

and the sufficient decrease writes

$$ F_k^{DR} - F_{k+1}^{DR} \ge \delta \|x_k - x_{k-1}\|^2 \quad (112) $$
with

$$ \delta = \frac{1}{2}(1 - L_{g_{\sigma}}) \left( \frac{1}{\beta} \left( 1 - \frac{L_{g_{\sigma}}^2}{(1 - L_{g_{\sigma}})^2} \right) - (1 + M) \right) \quad (113) $$

(ii) This point directly follows from the above sufficient decrease condition.

(iii) For the original DRS algorithm, this point was proven in [44, Theorem 2]. Given the sufficient decrease property, assuming that (a) $F^{DR}$ verifies the KL property on $\mathrm{Im}(D_{\sigma})$ and (b) that the iterates are bounded, the proof follows equally. With the same arguments as in the proof of Corollary 1, both assumptions (a) and (b) are verified when $f$ is real analytic and $g_{\sigma}$ is coercive. □

# 7 Experiments

In this section, we apply, with the proximal denoiser Prox-DRUNET, the PnP algorithms ProxPnP-PGD (55), ProxPnP-αPGD (86), ProxPnP-DRSdiff (61) (diff specifies that this PnP-DRS is dedicated to differentiable data-fidelity terms $f$) and ProxPnP-DRS (89) for deblurring and super-resolution with Gaussian noise. These four algorithms target critical points of the objective

$$ F = \lambda f + \phi_{\sigma} \quad (114) $$

with $f$ the data-fidelity term and $\phi_{\sigma}$ the regularization function detailed in Section 3.3.

For both applications, we consider a degraded observation $y = Ax^* + \nu \in \mathbb{R}^m$ of a clean image $x^* \in \mathbb{R}^n$ that is estimated by solving problem (1) with $f(x) = \frac{1}{2} \|Ax - y\|^2$. Its gradient $\nabla f(x) = A^T(Ax - y)$ is thus Lipschitz with constant $\|A^T A\|_S$. The blur kernels are normalized so that the Lipschitz constant $L_f = 1$. We use for evaluation and comparison the 68 images from the CBSD68 dataset, center-cropped to $n = 256 \times 256$ and Gaussian noise with 3 noise levels $\nu \in \{0.01, 0.03, 0.05\}$.

For deblurring, the degradation operator $A = H$ is a convolution performed with circular boundary conditions. As in [7, 12, 18, 24], we consider the 8 real-world camera shake kernels of [45], the 9 × 9 uniform kernel and the 25 × 25 Gaussian kernel with standard deviation 1.6.

For single image super-resolution (SR), the low-resolution image $y \in \mathbb{R}^m$ is obtained from the high-resolution one $x \in \mathbb{R}^n$ via $y = SHx + \nu$ where $H \in \mathbb{R}^{n \times n}$ is the convolution with anti-aliasing kernel. The matrix $S$ is the standard $s$-fold down-sampling matrix of size $m \times n$ and $n = s^2 \times m$. As in [7], we evaluate SR performance on 4 isotropic Gaussian blur kernels with standard deviations 0.7, 1.2, 1.6 and 2.0; and consider downsampled images at scale $s = 2$ and $s = 3$.

## Hyperparameter selection

As explained in Section 3.2, the proximal denoiser $D_{\sigma}$ defined in Proposition 5 is trained following [25] with a penalization encouraging $L_g < 1$. Convergence of ProxPnP-PGD, ProxPnP-αPGD, ProxPnP-DRSdiff and ProxPnP-DRS are then respectively guaranteed by Corollary 1, Corollary 2, [25, Theorem 4.3] and Theorem 3.
For each algorithm, in Table 1 we propose default values for the involved hyperpa-
rameters. Note that we use the same choice of hyperparameters for both deblurring
and super-resolution. The $\lambda$ parameter always satisfies the corresponding constraint
required for convergence. As each algorithm uses its own set of parameters depending
on the constraint on $\lambda$, each algorithm targets critical points of a different functional.

For ProxPnP-PGD, ProxPnP-αPGD and ProxPnP-DRS algorithms, we use the
$\gamma$-relaxed version of the denoiser (34). In practice, we found that the same choice
of parameters $\gamma$ and $\sigma$ is optimal for both PGD and $\alpha$PGD, with values depending
on the amount of noise $\nu$ in the input image. We thus choose $\lambda \in [0, \lambda_{lim}]$ where,
following Corollary 1 and Corollary 2, for ProxPnP-PGD $\lambda_{lim}^{\text{PGD}} = \frac{1}{L_f} \frac{\gamma+2}{\gamma+1}$ and for
ProxPnP-αPGD $\lambda_{lim}^{\text{αPGD}} = \frac{1}{L_f} \frac{\gamma+1}{\gamma} \geq \lambda_{lim}^{\text{PGD}}$. For both $\nu = 0.01$ and $\nu = 0.03$, $\lambda$ is set
to its maximal allowed value $\lambda_{lim}$. As $\lambda_{lim}^{\text{αPGD}} \geq \lambda_{lim}^{\text{PGD}}$, ProxPnP-αPGD is expected
to outperform ProxPnP-PGD at these noise levels. Finally, for ProxPnP-αPGD, $\alpha$ is
set to its maximum possible value $1/(\lambda L_f)$.

For ProxPnP-DRSdiff, $\lambda$ is also set to its maximal possible $\frac{1}{L_f}$ for theoretical
convergence [25, Theorem 4.3].

Eventually, for ProxPnP-DRS, Theorem 3 requires $L_{g_{\sigma}} < L_{max}(\beta)$ via the con-
straint (92). As $D_{\sigma}$ is trained to ensure $L_{g_{\sigma}} < 1$, we do not retrain the denoiser for a
specific $L_{max}(\beta)$ value, but we use again the $\gamma$-relaxed version of the denoiser $D_{\sigma}^{\gamma}$ (34)
with $\gamma = L_{max}(\beta)$. $\gamma \nabla g_{\sigma}$ is then $L_{max}(\beta)$-Lipschitz and $D_{\sigma}^{\gamma} = \text{Prox}_{\phi_{\sigma}^{\gamma}}$. In practice,
we find $\beta = 0.25$ to be a good compromise. For this choice of $\beta$, (92) is satisfied for
$L_{max}(\beta) \le 0.45$.

Numerical performance analysis

We numerically evaluate in Table 2, for deblurring and for super-resolution, the PSNR
performance of our four ProxPnP algorithms. We give comparisons with the deep
state-of-the-art PnP methods IRCNN [12] and DPIR [7] which both apply the PnP-
HQS algorithm with decreasing stepsize but without convergence guarantees. We also
provide comparisons with the GS-PnP (referring to PnP with Gradient Step denoiser)
method presented in [24].

Observe that, among ProxPnP methods, ProxPnP-DRS and ProxPnP-αPGD give
the best performance over the variety of kernels and noise levels. Indeed, ProxPnP-
DRS convergence is guaranteed for any of λ, which can thus be tuned to optimize
performance. ProxPnP-DRSdiff is, on the other hand, constrained to λ < 1. Simi-
larly, ProxPnP-PGD is constrained to λ < 2 when ProxPnP-αPGD restriction on λ
is relaxed. When a low amount of regularization (i.e. a large λ value) is necessary,
an upper bound on λ can severely limit the restoration capacity of the algorithm.
Indeed, we observe that when the input noise is low (ν = 2.55), ProxPnP-DRSdiff and
ProxPnP-PGD perform significantly worse. However, when the input noise is high,
a stronger regularization is necessary (i.e. a small λ value) and all methods perform
comparably.

We also provide visual comparisons for deblurring in Figure 1 and super-resolution
in Figure 2. Along with the output images, for each ProxPnP algorithm, we plot the
evolution of the corresponding provably-decreasing function. We recall in Table 3,
<table><caption>**Table 1:** Choice of the different hyperparameters involved for each ProxPnP algorithm. The same set of hyperparameters are used for both deblurring and super-resolution experiments.</caption><tr><td>ν(/255)</td><td></td><td>2.55</td><td>7.65</td><td>12.75</td></tr><tr><td rowspan="3">PGD</td><td>γ</td><td>0.6</td><td>1</td><td>1</td></tr><tr><td>λ = γ+2 / γ+1</td><td>1.625</td><td>1.5</td><td>1.5</td></tr><tr><td>σ/ν</td><td>1.25</td><td>0.75</td><td>0.5</td></tr><tr><td rowspan="4">αPGD</td><td>γ</td><td>0.6</td><td>1</td><td>1</td></tr><tr><td>λ = γ+1 / γ</td><td>2.66</td><td>2</td><td>2</td></tr><tr><td>α = 1 / λ</td><td>0.37</td><td>0.5</td><td>0.5</td></tr><tr><td>σ/ν</td><td>1.25</td><td>0.75</td><td>0.5</td></tr><tr><td rowspan="4">DRS</td><td>β</td><td>0.25</td><td>0.25</td><td>0.25</td></tr><tr><td>γ = L<sub>max</sub>(β)</td><td>0.45</td><td>0.45</td><td>0.45</td></tr><tr><td>λ</td><td>5</td><td>1.5</td><td>0.75</td></tr><tr><td>σ/ν</td><td>2</td><td>1</td><td>0.5</td></tr><tr><td rowspan="3">DRSdiff</td><td>β</td><td>0.5</td><td>0.5</td><td>0.5</td></tr><tr><td>λ</td><td>1.</td><td>1.</td><td>1.</td></tr><tr><td>σ/ν</td><td>0.75</td><td>0.5</td><td>0.5</td></tr></table>

<table><caption>**Table 2:** PSNR (dB) results on CBSD68 for deblurring (left) and super-resolution (right). PSNR are averaged over 10 blur kernels for deblurring (left) and 4 blur kernels along various scales s for super-resolution (right).</caption><thead><tr><th rowspan="2">Noise level ν</th><th colspan="3">Deblurring</th><th colspan="6">Super-resolution</th></tr><tr><th>0.01</th><th>0.03</th><th>0.05</th><th colspan="3">scale s = 2</th><th colspan="3">scale s = 3</th></tr><tr><th></th><th>0.01</th><th>0.03</th><th>0.05</th><th>0.01</th><th>0.03</th><th>0.05</th><th>0.01</th><th>0.03</th><th>0.05</th></tr></thead><tbody><tr><td>IRCNN [12]</td><td>31.42</td><td>28.01</td><td>26.40</td><td>26.97</td><td>25.86</td><td>25.45</td><td>25.60</td><td>24.72</td><td>24.38</td></tr><tr><td>DPIR [7]</td><td><strong>31.93</strong></td><td><strong>28.30</strong></td><td><u>26.82</u></td><td>27.79</td><td>26.58</td><td><u>25.83</u></td><td><u>26.05</u></td><td><u>25.27</u></td><td><u>24.66</u></td></tr><tr><td>GS-PnP [24]</td><td><u>31.70</u></td><td><u>28.28</u></td><td><strong>26.86</strong></td><td>27.88</td><td><strong>26.81</strong></td><td><strong>26.01</strong></td><td>25.97</td><td><strong>25.35</strong></td><td><strong>24.74</strong></td></tr><tr><td>ProxPnP-PGD</td><td>30.91</td><td>27.97</td><td>26.66</td><td>27.68</td><td>26.57</td><td>25.81</td><td>25.94</td><td>25.20</td><td>24.62</td></tr><tr><td>ProxPnP-αPGD</td><td>31.55</td><td>28.03</td><td>26.66</td><td><u>27.92</u></td><td><u>26.61</u></td><td>25.80</td><td><u>26.03</u></td><td>25.26</td><td>24.61</td></tr><tr><td>ProxPnP-DRSdiff</td><td>30.56</td><td>27.78</td><td>26.61</td><td>27.44</td><td>26.58</td><td>25.82</td><td>25.75</td><td>25.19</td><td>24.63</td></tr><tr><td>ProxPnP-DRS</td><td>31.51</td><td>28.01</td><td>26.63</td><td><strong>27.95</strong></td><td>26.58</td><td>25.81</td><td><strong>26.13</strong></td><td><u>25.27</u></td><td>24.65</td></tr></tbody></table>

for each algorithm, the formula of each function that is proved to decrease and converge along the iterations. We also plot the evolution of the norm of the residuals and of the PSNR along the iterates. These plots empirically confirm the theoretical convergence results. Observe that, despite being trained with additional constraints to guarantee convergence, ProxPnP-αPGD and ProxPnP-DRS globally compare with the performance of the state-of-the-art DPIR method.
<table><caption>**Table 3:** For each ProxPnP algorithm, the corresponding convergence theorem exhibits a function that is proved to decrease along the iterates.</caption><thead><tr><th>Algorithm</th><th>Decreasing function</th></tr></thead><tbody><tr><td>ProxPnP-PGD</td><td>F(x) = λf(x) + φ<sub>σ</sub>(x)</td></tr><tr><td>ProxPnP-αPGD</td><td>F<sup>α</sup>(y) = λf(y) + φ<sub>σ</sub>(y) + <sup>α</sup>&frasl;<sub>2</sub>(1 − <sup>1</sup>&frasl;<sub>α</sub>)<sup>2</sup>||y − y<sub>−</sub>||<sup>2</sup></td></tr><tr><td>ProxPnP-DRSdiff</td><td>F<sup>DR</sup>(x, y, z) = λf(y) + φ<sub>σ</sub>(z) + &lang;y − x, y − z&rang; + <sup>1</sup>&frasl;<sub>2</sub>||y − z||<sup>2</sup></td></tr><tr><td>ProxPnP-DRS</td><td>F<sup>DR</sup>(x, y, z) = λf(z) + φ<sub>σ</sub>(y) + &lang;y − x, y − z&rang; + <sup>1</sup>&frasl;<sub>2</sub>||y − z||<sup>2</sup></td></tr></tbody></table>

## 8 Conclusion

In this paper, we propose new convergent plug-and-play algorithms with minimal restrictions on the parameters of the inverse problem. We study two distinct schemes, the first one is a relaxed version of the Proximal Gradient Descent (PGD) algorithm and the second one is the Douglas Rachford Splitting (DRS) scheme. When used with a proximal denoiser, while the original PnP-PGD and PnP-DRS convergence results imposed restrictive conditions on the regularization parameter of the problem, the proposed algorithms converge, with milder conditions, towards stationary points of an explicit functional. Numerical experiments exhibit the convergence and the efficiency of the methods on deblurring and super-resolution problems.

A limitation of these approaches is that the γ-relaxation of the Gradient Step denoiser adds again an extra hyperparameter that needs to be tuned. For instance, PnP-αPGD has in total four hyperparameters to tune: σ, λ, γ and α. Even though we proposed default values in Section 7, this could be challenging and time-consuming for the user.

Different from αPGD, another possible modification of the PGD algorithm is the Primal-Dual Davis-Yin (PDDY) algorithm [46]. It is a Primal Dual version of the Davis-Yin splitting method [47]. As shown by Condat and Richtárik [48] (Algorithm 3), it can be written as a relaxation of the PGD algorithm with a parameter α ∈ (0, 1). From [48], convex convergence of the ProxPnP version of this algorithm could be established, for any regularization parameter λ, by choosing α small enough. However, the convergence of the PDDY algorithm in a nonconvex setting has not been established and remains an open research avenue for future investigations.
<figure><img src="image_4.png" /><figcaption>**Fig. 1:** Deblurring with various methods of “starfish” degraded with the indicated blur kernel and input noise level $\nu = 0.03$. We also plot for each ProxPnP algorithm the evolution of the respective decreasing Lyapunov functions, the residual $\gamma_k = \min_{0 \le i \le k} \frac{\|x_{i+1} - x_i\|^2}{\|x_0\|^2}$ and the PSNR.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 2:** Super-resolution with various methods on “leaves” downsampled by 2, with the indicated blur kernel and input noise level $\nu = 0.03$. We plot the evolution of the respective Lyapunov functions, the residual $\gamma_k = \min_{0 \le i \le k} \frac{\|x_{i+1} - x_i\|^2}{\|x_0\|^2}$ and the PSNR.</figcaption></figure>
# Appendix A Primal-Dual derivation of the αPGD algorithm

Suppose we target a solution of the following minimization problem:

$$ \min_{x \in \mathbb{R}^n} f(x) + \phi(x) \qquad (A1) $$

for $f$ $L_f$-smooth and strongly convex and $\phi$ convex. The method in [42] targets a minimizer with a Bregman Primal-Dual algorithm

$$ \left\{ \begin{array}{l} y_{k+1} = \arg \min_y \frac{1}{\sigma} D_{h^Y}(y, y_k) + f^*(y) - \langle \bar{x}_k, y \rangle \\ x_{k+1} = \arg \min_x \frac{1}{\tau} D_{h^X}(x, x_k) + \phi(x) + \langle x, y_{k+1} \rangle \\ \bar{x}_{k+1} = x_{k+1} + \beta(x_{k+1} - x_k), \end{array} \right. \qquad (A2) $$

where $D_h$ here denotes the Bregman divergence associated to a convex potential $h$. For convex $f$ and $\phi$, provided that the potentials $h^X$ and $h^Y$ are 1-convex with respect to the norm $||.||^2$, convergence of (A2) towards a solution of (A1) is ensured as long as $\tau\sigma ||K^*K|| < 1$ [42, Remark 3].

Notice that $L_f f^*$ is 1-strongly convex with respect to the norm $||.||^2$. Indeed, $f$ smooth and convex with a $L_f$-Lipschitz gradient implies $f^* 1/L_f$ strongly-convex, i.e. $L_f f^*$ is 1-strongly convex. Then we can use $h^Y = L_f f^*$ (and $h^X = \frac{1}{2} ||.||^2$) and the algorithm becomes

$$ \left\{ \begin{array}{l} y_{k+1} = \arg \min_y \left( \frac{L_f}{\sigma} + 1 \right) f^*(y) - \langle \bar{x}_k + \frac{L_f}{\sigma} \nabla f^*(y^k), y \rangle \\ x_{k+1} \in \arg \min_x \frac{1}{2\tau} \|x - x^k\|^2 + \phi(x) + \langle x, y_{k+1} \rangle \\ \bar{x}_{k+1} = x_{k+1} + \beta(x_{k+1} - x_k). \end{array} \right. \qquad (A3) $$

The optimality condition for the first update is

$$ \left( \frac{L_f}{\sigma} + 1 \right) \nabla f^*(y_{k+1}) = \bar{x}_k + \frac{L_f}{\sigma} \nabla f^*(y^k). \qquad (A4) $$

Using the fact that $(\nabla f)^{-1} = \nabla f^*$, with the change of variable $y_k \leftarrow \nabla f^*(y_k)$, the previous algorithm can then be rewritten in a fully primal formulation:

$$ \left\{ \begin{array}{l} y_{k+1} = \frac{\bar{x}_k + \frac{L_f}{\sigma} y_k}{1 + \frac{L_f}{\sigma}} \\ x_{k+1} \in \text{Prox}_{\tau\phi}(x_k - \tau \nabla f(y_{k+1})) \\ \bar{x}_{k+1} = x_{k+1} + \beta(x_{k+1} - x_k). \end{array} \right. \qquad (A5) $$

For $\alpha = \frac{1}{1 + \frac{L_f}{\sigma}}$, $\beta = 1 - \alpha$, the algorithm writes

$$ \left\{ \begin{array}{l} y_{k+1} = \alpha \bar{x}_k + (1 - \alpha) y_k \\ x_{k+1} \in \text{Prox}_{\tau\phi}(x_k - \tau \nabla f(y_{k+1})) \\ \bar{x}_{k+1} = x_{k+1} + (1 - \alpha)(x_{k+1} - x_k). \end{array} \right. \qquad (A6) $$
or

$$ \begin{cases} y_{k+1} = \alpha x_k + (1-\alpha)y_k + \alpha(1-\alpha)(x_k - x_{k-1}) \\ x_{k+1} \in \text{Prox}_{\tau\phi}(x_k - \tau \nabla f(y_{k+1})) \end{cases} \qquad (A7) $$

The latter is equivalent to (63) with $y_k$ in place of $q_k$.

## Appendix B Stepsize condition in Theorem 2.

In Theorem 2, the stepsize condition

$$ \tau < \min \left( \frac{1}{\alpha L_f}, \frac{\alpha}{M} \right) \qquad (B8) $$

can be replaced by the slightly better bound

$$ \tau < \min \left( \frac{1}{\alpha L_f}, \frac{2\alpha}{\alpha^3 L_f + (2-\alpha)M} \right). \qquad (B9) $$

which give little numerical gain. For sake of completeness, we develop the proof here.

*Proof.* We keep the second term in (75) instead of just using its non-negativity. From the relation

$$ x_{k+1} - x_k = \frac{1}{\alpha}(y_{k+1} - y_k) + \left(1 - \frac{1}{\alpha}\right)(y_k - y_{k-1}), \qquad (B10) $$

by convexity of the squared $\ell_2$ norm, for $0 < \alpha < 1$, we have

$$ \|y_{k+1} - y_k\|^2 \le \alpha \|x_{k+1} - x_k\|^2 + (1-\alpha) \|y_k - y_{k-1}\|^2 \qquad (B11) $$

and

$$ \|x_{k+1} - x_k\|^2 \ge \frac{1}{\alpha} \|y_{k+1} - y_k\|^2 + \left(1 - \frac{1}{\alpha}\right) \|y_k - y_{k-1}\|^2. \qquad (B12) $$

which gives finally

$$ \begin{aligned} F(y_k) - F(y_{k+1}) &\ge \left( \alpha \left(1 - \frac{1}{\alpha}\right) \left( \frac{1}{2\tau} - \frac{\alpha L_f}{2} \right) - \frac{\alpha}{2\tau} \left(1 - \frac{1}{\alpha}\right)^2 \right) \|y_k - y_{k-1}\|^2 \\ &\quad + \left( \frac{1}{2\tau} - \frac{\alpha L_f}{2} + \frac{1}{\alpha} \left( \frac{1}{2\tau} - \frac{M(2-\alpha)}{2} \right) \right) \|y_{k+1} - y_k\|^2 \\ &\quad = -\frac{1-\alpha}{2\alpha\tau} \left(1 - \alpha^2 \tau L_f\right) \|y_k - y_{k-1}\|^2 \\ &\quad + \frac{1}{2\alpha\tau} \left(1 + \alpha - \alpha^2 \tau L_f - \tau M(2-\alpha)\right) \|y_{k+1} - y_k\|^2 \\ &\quad = -\delta \|y_k - y_{k-1}\|^2 + \delta \|y_{k+1} - y_k\|^2 + (\gamma - \delta) \|y_{k+1} - y_k\|^2 \end{aligned} \qquad (B13) $$

with

$$ \delta = \frac{1-\alpha}{2\alpha\tau} \left(1 - \alpha^2 \tau L_f\right) \qquad (B14) $$
$$ \gamma = \frac{1}{2\alpha\tau} (1 - \alpha^2\tau L_f + \alpha - \tau M(2 - \alpha)). \qquad (B15) $$

The condition on the stepsize becomes

$$ \begin{gathered} \gamma - \delta > 0 \\ \Leftrightarrow \tau < \frac{2\alpha}{\alpha^3 L_f + (2 - \alpha)M}, \end{gathered} \qquad (B16) $$

and the overall condition is

$$ \tau < \min \left( \frac{1}{\alpha L_f}, \frac{2\alpha}{\alpha^3 L_f + (2 - \alpha)M} \right). \qquad (B17) $$

□

## References

[1] Rudin, L.I., Osher, S., Fatemi, E.: Nonlinear total variation based noise removal algorithms. Phys. D 60, 259–268 (1992)

[2] Tikhonov, A.N.: On the solution of ill-posed problems and the method of regularization. In: Doklady Akademii Nauk, vol. 151, pp. 501–504 (1963). Russian Academy of Sciences

[3] Mallat, S.: A Wavelet Tour of Signal Processing, The Sparse Way, 3rd edition edn. Academic Press, Elsevier, ??? (2009)

[4] Zoran, D., Weiss, Y.: From learning models of natural image patches to whole image restoration. In: IEEE ICCV, pp. 479–486 (2011)

[5] Bora, A., Jalal, A., Price, E., Dimakis, A.G.: Compressed sensing using generative models. arXiv preprint arXiv:1703.03208 (2017)

[6] Romano, Y., Elad, M., Milanfar, P.: The little engine that could: Regularization by denoising (red). SIAM Journal Imaging Sci. 10(4), 1804–1844 (2017)

[7] Zhang, K., Li, Y., Zuo, W., Zhang, L., Van Gool, L., Timofte, R.: Plug-and-play image restoration with deep denoiser prior. IEEE Transactions on Image Processing (2021)

[8] Attouch, H., Bolte, J., Svaiter, B.F.: Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods. Math. Programming 137(1-2), 91–129 (2013)
[9] Themelis, A., Patrinos, P.: Douglas–rachford splitting and ADMM for nonconvex optimization: Tight convergence results. SIAM Journal Optimization 30(1), 149–181 (2020)

[10] Venkatakrishnan, S.V., Bouman, C.A., Wohlberg, B.: Plug-and-play priors for model based reconstruction. In: IEEE Glob. Conf. Signal Inf. Process., pp. 945–948 (2013)

[11] Meinhardt, T., Moller, M., Hazirbas, C., Cremers, D.: Learning proximal operators: Using denoising networks for regularizing inverse imaging problems. In: Proceedings of the IEEE International Conference on Computer Vision, pp. 1781–1790 (2017)

[12] Zhang, K., Zuo, W., Gu, S., Zhang, L.: Learning deep CNN denoiser prior for image restoration. In: IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 3929–3938 (2017)

[13] Sun, Y., Liu, J., Kamilov, U.S.: Block coordinate regularization by denoising. arXiv preprint arXiv:1905.05113 (2019)

[14] Ahmad, R., Bouman, C.A., Buzzard, G.T., Chan, S., Liu, S., Reehorst, E.T., Schniter, P.: Plug-and-play methods for magnetic resonance imaging: Using denoisers for image recovery. IEEE signal processing magazine 37(1), 105–116 (2020)

[15] Yuan, X., Liu, Y., Suo, J., Dai, Q.: Plug-and-play algorithms for large-scale snapshot compressive imaging. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) (2020)

[16] Sun, Y., Wu, Z., Xu, X., Wohlberg, B., Kamilov, U.S.: Scalable plug-and-play admm with convergence guarantees. IEEE Transactions Computational Imaging 7, 849–863 (2021)

[17] Ryu, E., Liu, J., Wang, S., Chen, X., Wang, Z., Yin, W.: Plug-and-play methods provably converge with properly trained denoisers. In: International Conference on Machine Learning, pp. 5546–5557 (2019)

[18] Pesquet, J.-C., Repetti, A., Terris, M., Wiaux, Y.: Learning maximally monotone operators for image recovery. SIAM Journal Imaging Sci. 14(3), 1206–1237 (2021)

[19] Hertrich, J., Neumayer, S., Steidl, G.: Convolutional proximal neural networks and plug-and-play algorithms. Linear Algebra and its Applications 631, 203–234 (2021)

[20] Bohra, P., Goujon, A., Perdios, D., Emery, S., Unser, M.: Learning lipschitz-controlled activation functions in neural networks for plug-and-play image reconstruction methods. In: NeurIPS Workshop on Deep Learning and Inverse
[21] Nair, P., Chaudhury, K.N.: On the construction of averaged deep denoisers for image regularization. arXiv preprint arXiv:2207.07321 (2022)

[22] Sreehari, S., Venkatakrishnan, S.V., Wohlberg, B., Buzzard, G.T., Drummy, L.F., Simmons, J.P., Bouman, C.A.: Plug-and-play priors for bright field electron tomography and sparse interpolation. IEEE Transactions on Computational Imaging 2(4), 408–423 (2016)

[23] Cohen, R., Blau, Y., Freedman, D., Rivlin, E.: It has potential: Gradient-driven denoisers for convergent solutions to inverse problems. In: Neural Information Processing Systems, vol. 34 (2021)

[24] Hurault, S., Leclaire, A., Papadakis, N.: Gradient step denoiser for convergent plug-and-play. In: International Conference on Learning Representations (2022)

[25] Hurault, S., Leclaire, A., Papadakis, N.: Proximal denoiser for convergent plug-and-play optimization with nonconvex regularization. In: International Conference on Machine Learning (2022)

[26] Hurault, S., Chambolle, A., Leclaire, A., Papadakis, N.: A relaxed proximal gradient descent algorithm for convergent plug-and-play with proximal denoiser. In: International Conference on Scale Space and Variational Methods in Computer Vision, pp. 379–392 (2023). Springer

[27] Gribonval, R., Nikolova, M.: A characterization of proximity operators. Journal of Mathematical Imaging and Vision 62(6), 773–789 (2020)

[28] Terris, M., Repetti, A., Pesquet, J.-C., Wiaux, Y.: Building firmly nonexpansive convolutional neural networks. In: IEEE International Conference on Acoustics, Speech, and Signal Processing, pp. 8658–8662 (2020)

[29] Rockafellar, R.T., Wets, R.J.-B.: Variational Analysis vol. 317. Springer, ??? (2009)

[30] Attouch, H., Bolte, J., Redont, P., Soubeyran, A.: Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the kurdyka-lojasiewicz inequality. Mathematics of operations research 35(2), 438–457 (2010)

[31] Ochs, P., Chen, Y., Brox, T., Pock, T.: ipiano: Inertial proximal algorithm for nonconvex optimization. SIAM Journal Imaging Sci. 7(2), 1388–1419 (2014)

[32] Bolte, J., Sabach, S., Teboulle, M., Vaisbourd, Y.: First order methods beyond convexity and lipschitz gradient continuity with applications to quadratic inverse problems. SIAM Journal on Optimization 28(3), 2131–2151 (2018)
[33] Zeng, J., Lau, T.T.-K., Lin, S., Yao, Y.: Global convergence of block coordinate descent in deep learning. In: International Conference on Machine Learning, pp. 7313–7323 (2019). PMLR

[34] Krantz, S.G., Parks, H.R.: A Primer of Real Analytic Functions, (2002)

[35] Lojasiewicz, S.: Ensembles semi-analytiques. IHES notes, 220 (1965)

[36] Bauschke, H.H., Combettes, P.L.: Convex Analysis and Monotone Operator Theory in Hilbert Spaces, p. 468 (2011)

[37] Beck, A., Teboulle, M.: Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems. IEEE Transactions on Image Processing **18**(11), 2419–2434 (2009)

[38] Tseng, P.: On accelerated proximal gradient methods for convex-concave optimization. submitted to SIAM Journal Optimization **2**(3) (2008)

[39] Nesterov, Y.: Gradient methods for minimizing composite functions. Mathematical programming **140**(1), 125–161 (2013)

[40] Lan, G., Zhou, Y.: An optimal randomized incremental gradient method. Mathematical programming **171**(1), 167–215 (2018)

[41] Chambolle, A., Pock, T.: A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Math. Imaging and Vision **40**(1), 120–145 (2011)

[42] Chambolle, A., Pock, T.: On the ergodic convergence rates of a first-order primal-dual algorithm. Mathematical Programming **159**(1), 253–287 (2016)

[43] Bauschke, H.H., Combettes, P.L.: Convex Analysis and Monotone Operator Theory in Hilbert Spaces vol. 408. Springer, ??? (2011)

[44] Li, G., Pong, T.K.: Douglas–rachford splitting for nonconvex optimization with application to nonconvex feasibility problems. Math. Progr. **159**, 371–401 (2016)

[45] Levin, A., Weiss, Y., Durand, F., Freeman, W.T.: Understanding and evaluating blind deconvolution algorithms. In: IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 1964–1971 (2009)

[46] Salim, A., Condat, L., Mishchenko, K., Richtárik, P.: Dualize, split, randomize: Toward fast nonsmooth optimization algorithms. Journal of Optimization Theory and Applications **195**(1), 102–130 (2022)

[47] Davis, D., Yin, W.: A three-operator splitting scheme and its optimization applications. Set-valued and variational analysis **25**, 829–858 (2017)
[48] Condat, L., Richtárik, P.: Randprox: Primal-dual optimization algorithms with randomized proximal updates. arXiv preprint arXiv:2207.12891 (2022)