1 Introduction↩︎

The restoration of high-quality images from corrupted or incomplete measurements is a fundamental task in modern computational imaging, encompassing applications from medical reconstruction to microscopic and astronomical imaging [1]. Mathematically, the observation process is typically modeled as a linear or nonlinear forward process subject to noise: \[y = \mathcal{N}(A(x)),\] where \(x \in \mathbb{R}^n\) is the unknown underlying image, \(y \in \mathbb{R}^m\) represents the noisy observation as a realization of the noise model, \(A : \mathbb{R}^n \to \mathbb{R}^m\) is the forward operator (e.g., a blurring kernel, subsampled Fourier transform, or Radon transform), and \(\mathcal{N}\) denotes the noise model.

Since the operator \(A\) is often ill-conditioned or possesses a nontrivial null space, recovering \(x\) from \(y\) is an ill-posed inverse problem. To overcome this, the classical framework relies on Bayesian inference. By Bayes’ theorem, the posterior probability of the image given the observation is: \[P(x|y) = \frac{P(y|x)P(x)}{P(y)} \propto P(y|x)P(x).\] The most common approach to estimate \(x\) is to find the Maximum A Posteriori (MAP) estimator, which maximizes \(P(x|y)\), or equivalently, minimizes its negative logarithm. This yields the composite optimization problem: \[\label{prob:initial} \min_{x \in \mathbb{R}^n} F(x) \equiv f_{data}(x) + \lambda \mathcal{R}(x),\tag{1}\] where \(f_{data}(x) = -\log P(y|x)\) ensures data fidelity (e.g., the least-squares \(\frac{1}{2}\left\| A x - y \right\|^2\) under Gaussian noise), and \(\mathcal{R}(x) = -\log P(x)\) is a regularization term encoding prior knowledge about the statistical distribution of images we wish to recover. The parameter \(\lambda > 0\) controls the trade-off between fidelity to the measurements and the prior. Historically, analytical priors were widely used to address this problem, with milestones such as Tikhonov regularization [2] and the Total Variation (TV) model [3].

When the regularizer \(\mathcal{R}\) is convex but nonsmooth (as is the case with \(\ell_1\)-norm sparsity priors or TV), problem 1 cannot be solved via standard gradient descent. Instead, it is traditionally addressed via proximal gradient schemes [4], [5], which alternate a forward gradient step on the smooth fidelity term with a backward proximal step on the regularizer. The iterative scheme is given by: \[\label{eq:fb95scheme} \begin{align} &z^k = x^k - \alpha_k \nabla f_{data}(x^k), \\ &x^{k+1} = \text{prox}_{\alpha_k \lambda \mathcal{R}}(z^k), \end{align}\tag{2}\] where \(\alpha_k > 0\) is the step size and the proximity operator is defined as \[\nonumber \text{prox}_{\tau \mathcal{R}}(z) = \mathop{\mathrm{arg\,min}}_{x} \left\{ \mathcal{R}(x) + \frac{1}{2\tau}\left\| x - z \right\|^2 \right\}.\] Notably, evaluating the proximity operator is mathematically equivalent to solving a pure denoising problem on \(z^k\) with assumed Gaussian noise level proportional to \(\sqrt{\tau}\).

While analytical priors are mathematically tractable, they struggle to capture the complex, high-dimensional manifolds of natural images. Then, in the last decade, the surge in the availability of large amounts of data, the development of statistical learning theory [6], and the impressive results achieved by neural network models have also influenced the inverse problems community, giving rise to data-driven approaches [7], [8]. Examples of these methods are unrolled optimization algorithms [9], [10] and Plug-and-Play (PnP) [11] schemes. We focus on the latter, which has become a dominant paradigm in computational imaging due to its modularity and empirical performance.

Framework of the contribution. The Plug-and-Play (PnP) framework was introduced as a way to leverage the powerful denoising capabilities of modern Deep Neural Networks (DNNs) within the structure of traditional optimization algorithms. In the seminal paper [12] it was observed that proximal step in 2 acts simply as an image denoiser. They proposed to replace the proximity operator with an off-the-shelf, highly expressive pre-trained deep denoiser, which we denote as \(D_{\sigma}\). The PnP forward–backward (FB) algorithm simply becomes \[\label{eq:pnp95fb95scheme} x^{k+1} =D_{\sigma}( x^k - \alpha_k \nabla f_{data}(x^k)).\tag{3}\] PnP methods have since demonstrated state-of-the-art empirical performance across numerous imaging domains. However, divorcing the denoising step from a strictly defined optimization objective means that the algorithm is no longer guaranteed to minimize a specific global function \(F(x)\), making theoretical convergence analysis extremely challenging.

To bridge the gap between empirical success and theoretical guarantees, a major breakthrough was achieved in establishing the existence of regularizers parameterized by neural networks for which the proximal operator is exactly a denoiser. A remarkable result in this framework was taken in [13], where the authors introduced the Gradient-Step (GS) denoiser which takes the form of a gradient descent step: \[D_{\sigma}= \mathrm{Id}- \nabla g_{\sigma}, \qquad g_{\sigma}= \frac{1}{2}\left\| x - N_\sigma(x) \right\|^2,\] being \(\mathrm{Id}\) the identity operator and \(N_\sigma\) a parametrized differentiable neural network. Moreover, in the same paper the following FB scheme applied to the functional \(f_{data} + \lambda g_{\sigma}\) was proposed, \[x^{k+1} = \text{prox}_{\tau f_{data}} \left(\tau \lambda D_{\sigma}(x^{k}) + (1 - \tau \lambda) x^{k}\right) = \text{prox}_{\tau f_{data}} \left(x^{k} - \tau \lambda \nabla g_\sigma(x^{k})\right),\] where the proximal step is computed with respect to the (convex) data-fidelity term. This scheme turns out to be nonpractical in many imaging inverse problems, since the proximal operator of the data-fidelity term is not computable in closed form. In the follow-up work [14], given the constraint that \(\nabla g_{\sigma}\) is \(L_{g_{\sigma}}\)–Lipschitz with \(L_{g_{\sigma}} < 1\), the same authors proved that the GS denoiser corresponds to the proximity operator of an explicit prox-regular function \(\phi_{\sigma}\), i.e., \(D_\sigma = \text{prox}_{\phi_\sigma}\), with \(\text{prox}_{\phi_\sigma}\) single valued. Hence, they consider iteration 3 with \(\alpha_k = \lambda\), showing that it corresponds to a FB step applied to the functional \(\lambda f_{data} + \phi_{\sigma}\) with a fixed step size equal to one. The convergence of this PnP-FB method is then proved for possibly nonconvex \(f_{data}\) under a constraint on the regularization parameter \(\lambda\), which must satisfy \(\lambda L_{f_{data}} <1\), where \(L_{f_{data}}\) is the Lipschitz constant of \(\nabla f_{data}\). Another step further was taken in [15], where the authors proved that the regularizer \(\phi_{\sigma}\) can be expressed for some constant \(K \in \mathbb{R}\) as \[\label{eq:hurault95phi} \phi_{\sigma}(x) = h_{\sigma}^*(x) - \frac{1}{2}\left\| x \right\|^2 + K, \qquad \forall x \in \mathop{\mathrm{Im}}(D_{\sigma}),\tag{4}\] where \(h_{\sigma}\) satisfies \(D_{\sigma}=\nabla h_{\sigma}\) and \(h_{\sigma}^*\) denotes its convex conjugate. Formula 4 shows that the function \(\phi_{\sigma}\) is weakly convex, allowing the authors to perform their convergence analysis under the assumption \(\lambda L_{f_{data}} < \frac{L_{g_{\sigma}} + 2}{L_{g_{\sigma}} + 1}\). Moreover, by formulating a relaxed version of the PnP-FB scheme, they proved convergence in the convex setting further lightening the constraints on \(\lambda\).

Despite the monumental theoretical step provided by the GS denoiser, the derived PnP algorithms based on the FB splitting in 3 still suffer from major limitations

1. when the data-fidelity is nonconvex, convergence is only guaranteed with the constraint \(\lambda L_{f_{data}} < \frac{L_{g_{\sigma}} + 2}{L_{g_{\sigma}} + 1}\), which is a strong limitation in case of low-noise levels and large Lipschitz constants of the data-fidelity gradient;

2. the constraints can be relaxed only in the convex setting, which limits the applicability of the method, and even in this case computing the exact bound is challenging;

3. due to the fact that \(\text{prox}_{\tau \phi_{\sigma}}\) is not computable in closed form for an arbitrary scaling factor \(\tau \neq 1\), the step–size of the algorithm is fixed and equal to one, which is not optimal for convergence speed.

As an additional difficulty, the function \(\phi_\sigma\) can be evaluated in closed form only on the points of the set \(\mathrm{Im}(D_\sigma)\).

Proposed approach. In this work, we propose a PnP-FB algorithm that overcomes all the limitations mentioned above while preserving the state-of-the-art convergence properties. The key point is a new way to split the objective functional \(f_{data}+\lambda\phi_\sigma\), based on the dual formulation 4 . Our splitting requires the computation of \(\text{prox}_{\tau h^*_\sigma}\), which is not available in closed form but can be approximated within a given tolerance by an inner procedure based on the evaluation of the denoiser \(D_\sigma\). The inner procedure is tailored such that it outputs a point where the function \(h^*_\sigma\) can be exactly evaluated. This allows us to design a line–search loop to adaptively compute the step length parameter.

The overall algorithm therefore takes the form of a line–search based inexact FB method, whose convergence analysis builds upon the abstract convergence scheme developed in [16]. In particular, we prove that any limit point of the iterate sequence is a stationary point of the function \(f_{data}+\lambda\phi_\sigma\) and, adopting the classical Kurdyka–Łojasiewicz (KL) [17] framework, we are able to guarantee also global convergence, without any restriction on the regularization parameter \(\lambda\) and without requiring convexity of the data-fidelity term.

Related work. The pursuit of provably convergent PnP methods has been a primary focus of the computational imaging community. Early attempts focused on strictly constraining the Lipschitz constant of the denoiser, as seen in [18]. A parallel approach emerged through Regularization by Denoising (RED) [19], which defined an explicit objective using the denoiser directly, though later works demonstrated that RED acts as a formal gradient only under highly restrictive conditions [20]. More recently, the framework has been extended to denoising diffusion models [21], [22] and flow-matching denoisers [23].

When proximity operators lack closed-form solutions, exact FB algorithms cannot be applied and the proximity operator must be computed inexactly. The foundational convergence rates of inexact proximal-gradient methods were established in [24] and later generalized in [25]. However, allowing inexactness while simultaneously applying a line–search is a nontrivial challenge because the line–search acceptance criteria generally require exact evaluation of the objective function. Crucial progress was made by Bonettini et al. [26], who designed variable-metric and line–search strategies for inexact operators. In the last years, this culminated in a generalized abstract KL convergence frameworks [16] and in the PHILA algorithm [27]. In our work, we bridge this optimization theory with the PnP imaging paradigm, tailoring the inexact KL framework to absorb the challenges posed by the neural network-based weakly-convex regularizer.

Outline of the paper. The paper is organized as follows. In Section 2, we present the general formulation of the inverse problem, highlighting the role and properties of the functions involved. In Section 3, we describe in detail our novel line–search-based Plug-and-Play FB algorithm, with particular emphasis on the inexact computation of the proximal point. The convergence analysis of the proposed method is presented in Section 4. Finally, Section 5 reports numerical experiments on image deblurring problems in the presence of different noise statistics, while our final remarks are offered in Section 6.

2 Formulation of the problem and novel splitting strategy↩︎

As explained in the introduction, we are interested in solving problem 1 with a PnP-FB scheme where the proximal step is replaced by a Gradient-Step denoiser \(D_{\sigma}\). Hence, we minimize the functional \[\label{eq:F95first} F(x) = f_{data}(x) + \lambda \phi_{\sigma}(x),\tag{5}\] where \(\phi_{\sigma}\) is defined as in 4 . In the following, we will always make the following assumption:

In the following proposition we recall the properties of \(\phi_{\sigma}\).

The choice of \(g_{\sigma}\) is crucial to ensure the existence of the implicit regularizer \(\phi_{\sigma}\) and to obtain convergence guarantees. In particular, it is required that:

1. \(g_{\sigma}(x) = \frac{1}{2}\left\| x - N_\sigma(x) \right\|^2\), where \(N_\sigma\) is a DRUNet with softplus activation functions;

2. the network is trained to impose that \(\nabla g_{\sigma}\) is \(L_{g_{\sigma}}\)–Lipschitz continuous, with \(L_{g_{\sigma}} < 1\).

We refer to [14] for more details on the training procedure and the choice of the architecture. For our purposes, we focus on the properties of \(\phi_{\sigma}\) and \(h_{\sigma}\) that are relevant for the convergence analysis. In particular, we point out that \(D_{\sigma}\) is a Lipschitz continuous operator with constant \(L_{D_{\sigma}} = L_{g_{\sigma}} + 1< 2\).

A natural FB approach to minimize the function \(F\) in 5 associated to the splitting \(f_{data}, \lambda\phi_\sigma\) requires the evaluation of the proximity operator associated to \(\tau\phi_\sigma\), where \(\tau\) is the product of the regularization parameter \(\lambda\) with the stepsize used in the forward step. However, although we can evaluate \(D_{\sigma}= \text{prox}_{\tau \phi_{\sigma}}\) for \(\tau = 1\), we cannot compute \(\text{prox}_{\tau \phi_{\sigma}}\) in closed form for \(\tau \neq 1\), because this would require inverting the denoiser \(D_{\sigma}\), which is intractable for Deep Neural Networks. Hence, we propose a novel splitting strategy that circumvents this issue and computes an approximation of \(\text{prox}_{\lambda h^*_\sigma}\) via an inner optimization procedure on the primal function \(h_{\sigma}\), whose gradient is the denoiser \(D_{\sigma}\). More precisely, in view of ?? , we express the objective function 5 as \[\label{eq:F} F(x) = f_{data}(x) + \lambda \left( h^*_{\sigma}(x) - \frac{1}{2}\left\| x \right\|^2 \right) \qquad \forall x \in \mathbb{R}^n\tag{6}\] and we split it as follows: \[\begin{align} f_0(x) &= f_{data}(x) - \frac{\lambda}{2}\left\| x \right\|^2 \quad &\text{(smooth, possibly nonconvex)} \tag{7} \\ f_1(x) &= \lambda h_{\sigma}^*(x) \quad &\text{(convex)} \tag{8} \end{align}\] Observe that by Assumption 1, \(f_0\) is smooth with \(L_{f_0}\)–Lipschitz continuous gradient, where \(L_{f_0} = L_{f_{data}} + \lambda\). Moreover, since \(h_\sigma\) is strongly convex, \(h^*_\sigma\) is differentiable with Lipschitz continuous gradient. Based on the splitting 7 –8 , we propose an inexact FB method to minimize the function \(F\) in 6 , as described in the following section.

3 Algorithm overview and inexact proximal computation↩︎

To minimize the objective function \(F(x)\) defined in our splitting strategy, we propose an iterative scheme that alternates between gradient steps, inexact proximal updates, and an adaptive line–search. We sketch the algorithmic framework: \[\begin{align} &z^k = x^k - \alpha_k \nabla f_0(x^k), \tag{9} \\ &\tilde{p}^k \approx \text{prox}_{h_{\sigma}/(\alpha_k \lambda )} \left (\frac{z^k}{\alpha_k \lambda} \right ), \tag{10} \\ &\tilde{y}^k = D_{\sigma}(\tilde{p}^k), \tag{11} \\ &x^{k+1} = x^k + \eta_k (\tilde{y}^k - x^k), \quad \eta_k \in [0, 1]. \tag{12} \end{align}\] Before diving into the technical details, we briefly describe the conceptual steps performed at each iteration:

1. Forward Step 9 . We perform a gradient step on the smooth part \(f_0\) to obtain an intermediate point \(z^k\). Our analysis is almost completely independent of the step–size sequence \(\{ \alpha_k \}_{k \in \mathbb{N}}\). For theoretical reasons, we must only impose \[\label{alphalim} 0 < \alpha_{\text{min}} \leq \alpha_k \leq \alpha_{\text{max}}, \quad \forall k \in \mathbb{N},\tag{13}\] for some positive constants \(\alpha_{\text{min}}, \alpha_{\text{max}}\). In practice, \(\alpha_k\) can be tuned to improve convergence speed.

2. Inexact Backward Step 10 –11 . We compute an inexact proximal update \(\tilde{p}^k\) for the convex part \(f_1\) by shifting the computation to the primal function \(h_{\sigma}\) via the Moreau decomposition. This step is the core focus of Sections 3.1–3.2.

3. Line–search Step 12 . We verify if the step from \(x^k\) to \(\tilde{y}^k\) provides a sufficient decrease of the objective. If not, we backtrack the step–size and repeat. Since we cannot evaluate \(f_1\) exactly, we introduce a novel surrogate merit function. This procedure is detailed in Section 3.3.

3.1 Proximal calculation via inner minimization↩︎

As already mentioned, we employ the Moreau decomposition [28] to shift the computation from the conjugate \(h_{\sigma}^*\) to the primal \(h_{\sigma}\): \[\label{eq:hatyk} \hat{y}^k := \text{prox}_{\alpha_k \lambda h_{\sigma}^*}(z^k) = z^k - \alpha_k \lambda \text{prox}_{h_{\sigma}/(\alpha_k \lambda )} \left (\frac{z^k}{\alpha_k \lambda} \right ) = z^k - \alpha_k \lambda \hat{p}^k,\tag{14}\] where \[\label{eq:inner95problem} \hat{p}^k := \text{prox}_{h_{\sigma}/(\alpha_k\lambda)}\left (\frac{z^k}{\alpha_k \lambda} \right ) = \mathop{\mathrm{arg\,min}}_{p \in \mathbb{R}^n} \underbrace{\frac{1}{2}\left\| p - \frac{z^k}{\alpha_k \lambda} \right\|^2 + \frac{1}{\alpha_k\lambda} h_{\sigma}(p)}_{:=G_k(p)}.\tag{15}\] Since \(h_{\sigma}\) is \((1 - L_{g_{\sigma}})\)–strongly convex, the function \(G_k\) in 15 is \(\mu_{G_k}\)–strongly convex with \[\label{eq:muG} \mu_{G_k} = 1 + \frac{1-L_{g_{\sigma}}}{\alpha_k\lambda}\geq 1.\tag{16}\] Moreover, the relation \(D_{\sigma}= \nabla h_{\sigma}\) leads to an explicit formula for its gradient: \[\label{eq:gradG} \nabla G_k(p) = p - \frac{z^k}{\alpha_k \lambda} + \frac{1}{\alpha_k \lambda} D_{\sigma}(p) \qquad \forall p \in \mathbb{R}^n.\tag{17}\] Hence, an approximate solution \(\tilde{p}^k\) of 15 can be computed e.g. by performing a gradient descent on \(G_k\) arrested when relation \[\label{eq:stop95crit95G} \| \nabla G_k(\tilde{p}^k) \| \leq \varepsilon_k\tag{18}\] is satisfied, where \(\varepsilon_k\) is a suitable positive tolerance.

As concerns the computation of an approximate solution \(\tilde{y}^k\) of 14 , we remark that the straightforward approach of applying the explicit formula \(\tilde{y}^k = z^k - \alpha_k \lambda \tilde{p}^k\) in 14 is not practical, since \(\tilde{y}^k\) would not be assured to belong to \(\mathrm{Im}(D_{\sigma})\); therefore, we would be missing an explicit way to compute \(f_1(\tilde{y}^k)\), which, as we will see later, is needed to define a sufficient decrease condition. To overcome this issue, first we observe that from equations 14 –17 and since \(\nabla G_k(\hat{p}^k) = 0\) one obtains that \[\nabla h_{\sigma}(\hat{p}^k) = \hat{y}^k.\] Therefore, the approximation \(\tilde{y}^k\) can be computed as \[\label{eq:ytildeDenoised} \tilde{y}^k = \nabla h_{\sigma}(\tilde{p}^k) = D_{\sigma}(\tilde{p}^k).\tag{19}\] Using this new formulation, we show that we can control the distance between the exact prox and its approximation: \[\nonumber \left\| \tilde{y}^k - \hat{y}^k \right\| = \left\| D_{\sigma}( \tilde{p}^k ) - D_{\sigma}( \hat{p}^k ) \right\| \leq L_{D_{\sigma}} \left\| \tilde{p}^k - \hat{p}^k \right\| \leq \frac{L_{D_{\sigma}}}{\mu_{G_k}} \left\| \nabla G_k(\tilde{p}^k) \right\| \leq \frac{L_{D_{\sigma}}}{\mu_{G_k}}\varepsilon_k,\] where the second inequality follows from the strong convexity of \(G_k\) and the last one from the stopping criterion 18 . Using the fact that \(\mu_{G_k} > 1\) and \(L_{D_\sigma} < 2\), this leads to \[\label{eq:distyk95ine} \left\| \tilde{y}^k - \hat{y}^k \right\| \leq 2\varepsilon_k.\tag{20}\] In addition, if 19 holds we are able to explicitly compute \(f_1(\tilde{y}^k)\) and \(h_\sigma(\tilde{p}^k)\). Indeed, from 19 and ?? , we derive the following explicit formula \[\begin{align} f_1(\tilde{y}^k) = \lambda h_{\sigma}^*(\tilde{y}^k) &= \lambda\Big( g_\sigma(D_\sigma^{-1}(\tilde{y}^k))-\frac{1}{2} \|D_\sigma^{-1}(\tilde{y}^k) -\tilde{y}^k\|^2 +\frac{1}{2} \|\tilde{y}^k\|^2 \Big)\\ &= \lambda \Big( g_\sigma(\tilde{p}^k) +\langle \tilde{y}^k, \tilde{p}^k \rangle - \frac{1}{2} \|\tilde{p}^k\|^2 \Big). \end{align}\] Moreover, by the Fenchel–Young inequality and the relation 19 we obtain [28] \[\nonumber h_{\sigma}^*(\tilde{y}^k) = \langle \tilde{y}^k, \tilde{p}^k \rangle - h_{\sigma}(\tilde{p}^k)\] which, combined with the previous relation, yields \(h_\sigma (\tilde{p}^k)= \frac{1}{2} \|\tilde{p}^k\|^2 - g_\sigma(\tilde{p}^k)\).

3.2 Properties of the inexact proximal gradient point↩︎

In this section, we state some crucial properties of the point \(\tilde{y}^k\) obtained as described in the previous section. To this end, we follow the notation in [27] and observe that \(\hat{y}^k = \text{prox}_{\alpha_k f_1}(z^k)\) is the unique point such that \[\label{eq:Hkfirstdef} \hat{y}^k = \mathop{\mathrm{arg\,min}}_{y \in \mathbb{R}^n} H_k(y), \qquad H_k(y) = \langle \nabla f_0(x^k), y - x^k \rangle + \frac{1}{2 \alpha_k } \left\| y - x^k \right\|^2 + f_1(y) - f_1(x^k).\tag{21}\] Notice that \(H_k(x^k) = 0\), therefore \(H_k(\hat{y}^k) \leq 0\), and \(H_k(\hat{y}^k) < 0\) whenever \(x^k\neq \hat{y}^k\). On the other hand, \(x^k=\hat{y}^k\) implies that \(x^k\) is a stationary point for the function \(F\), i.e., \(\nabla F(x^k) = 0\). In general, when \(f_1\) is convex, \(H_k\) is \(\theta\)–strongly convex, with \(\theta = 1/\alpha_{\max}\). Moreover, since in our case \(h_{\sigma}\) is \((1 - L_{g_{\sigma}})\)–strongly convex, the gradient of \(f_1=\lambda h_{\sigma}^*\) is Lipschitz continuous with constant \(\frac{\lambda}{1 - L_{g_{\sigma}}}\). Then, the gradient of \(H_k\) is \(\mu\)–Lipschitz continuous with \(\mu \geq 1/{\alpha_{\min}} + \lambda/(1-L_{g_\sigma })\).

The function \(H_k\) plays a crucial role in our algorithm since if \(H_k(\tilde{y} ^k) < 0\), the vector \(d^k = \tilde{y}^k - x^k\) computed in 11 is a descent direction for the objective function \(F\) [29]. However, the negative sign of \(H_k(\tilde{y}^k)\) alone does not guarantee suitable properties of the descent direction \(d^k\), needed to perform the convergence analysis of iteration 9 –12 . Such properties are actually obtained if the value \(H_k(\tilde{y}^k)\) is sufficiently close to the optimal one, \(H_k(\hat{y}^k)\). For this reason, it is important to estimate the distance from the optimal value. In our settings, if \(\tilde{y}^k\) is the point defined in 19 , with \(\tilde{p}^k\) satisfying 17 , the \(\mu\)–smoothness of \(H_k\) and the Descent Lemma imply \[\label{eq:ineqH} H_k (\tilde{y}^k) - H_k (\hat{y}^k) \leq \frac{\mu}{2} \left\| \tilde{y}^k - \hat{y}^k \right\|^2 \leq 2{\mu}{\varepsilon}_k^2.\tag{22}\] Hence, up to some constants, we can control the distance \(H_k (\tilde{y}^k) - H_k (\hat{y}^k)\) by means of the tolerance parameter \(\varepsilon_k\), without knowing the optimal value \(H_k (\hat{y}^k)\). Before giving more details on how the tolerance parameter \(\varepsilon_k\) is chosen, we explain our line–search approach, because the two issues are closely related.

3.3 Line–search procedure↩︎

To ensure convergence in the nonconvex setting, we employ an Armijo-like line search. Considering the descent direction \(d^k = \tilde{y}^k - x^k\), the new point is set as \[x^{k+1} = x^k + \eta_k d^k, \quad \eta_k \in (0,1],\] where the line–search parameter \(\eta_k\) must satisfy a sufficient decrease condition. A very common example of such conditions is based on the objective function values and writes as \[\label{eq:armijo0} F (x^{k}+\eta_k d^k) \leq F (x^k) + \omega \eta_k \Delta_k,\tag{23}\] where \(\Delta_k\) is a negative quantity expressing the expected decrease and \(\omega\in(0,1)\) is a fixed parameter [27]. The value of \(\eta_k\) is typically found with a backtracking procedure, i.e., by defining \(\eta_k = \delta^{m_k}\), with \(\delta \in (0,1)\) and \(m_k\) equal to the smallest non–negative integer for which 23 is satisfied. Moreover, borrowing the ideas in [26], [29], [30] the sufficient decrease quantity is usually set as \(\Delta_k = H_k(\tilde{y}^k)\). However, this approach cannot be directly adopted in our settings, due to the lack of an explicit way to compute \(f_1(x^k +\eta d^k)\) except for the case \(\eta = 1\), which is needed to obtain both the quantities \(F(x^k+\eta d^k)\) and \(H_k(\tilde{y}^k)\). We overcome this difficulty by introducing an upper estimate for \(f_1(x^k)\) in the following way. From the convexity of \(f_1\), for \(\eta \in (0,1]\) it holds that \[\nonumber f_1(x^k + \eta d^k) \leq \eta f_1(\tilde{y}^k) + (1 - \eta) f_1(x^k).\] Initializing \(U^0 = f_1(x^0)\) for some initial guess \(x^0\), at each step we can obtain an upper estimate of \(f_1(x^{k+1})\) as follows \[U^{k+1} := \eta_k f_1(\tilde{y}^k) + (1 - \eta_k) U^k.\] We next show that the sequence \(\{U^k\}_{k\in\mathbb{N}}\) provides an upper estimate of \(f_1(x^k)\) for all \(k\).

Once we have defined the sequence \(\{ U^k \}_{k \in \mathbb{N}}\), we introduce a new figure of merit as follows \[\label{eq:Phi95def} \Phi(x^k, U^k) = f_0(x^k) + U^k\tag{24}\] and we define a computable approximate value for \(\Delta_k=H_k(\tilde{y}^k)\) as \[\label{eq:def95tildeDeltak} \tilde{\Delta}_k = \langle \nabla f_0(x^k), \tilde{y}^k - x^k \rangle + \frac{1}{2 \alpha_k } \left\| \tilde{y}^k - x^k \right\|^2 + f_1(\tilde{y}^k) - U^k,\tag{25}\] which is, in practice, \(H_k(\tilde{y}^k)\) with \(f_1(x^k)\) replaced by \(U^k\). Clearly, the approximation 25 satisfies the inequality \[\label{eq:deltatildeineq} \tilde{\Delta}_k \leq H_k(\tilde{y}^k).\tag{26}\] Based on the previous definitions, we set the backtracking loop to compute the smallest non–negative integer \(m_k\) such that \[\label{eq:armijo} \min\{\Phi (x^{k}+\delta^{m_k} d^k, \delta^{m_k}f_1(\tilde{y}^k)+ (1-\delta^{m_k})U^k), F(\tilde{y}^k) \}\leq \Phi (x^k,U^k) + \omega \delta^{m_k} \tilde{\Delta}_k.\tag{27}\] The well posedness of the above condition is stated in Lemma 5 to follow. Its proof relies on a specific choice of the tolerance parameter \(\varepsilon_k\) in 18 , which is detailed in the next section. Finally, we adopt the following criterion to select the actual steplength parameter to be used in the updating rule 12 \[\label{eq:updating95complete} \eta_k = \begin{cases} \delta^{m_k} &if\Phi (x^{k}+\delta^{m_k} d^k, \delta^{m_k}f_1(\tilde{y}^k)+ (1-\delta^{m_k})U^k)< F(\tilde{y}^k)\\ 1 &otherwise \end{cases}.\tag{28}\] The above choice is motivated by the fact that it implies the following two properties on the point \(x^{k+1}\): \[\label{eq:line95search95cons} \Phi(x^{k+1},U^{k+1} )\leq \Phi(x^k,U^k) + \omega \delta^{m_k}\tilde{\Delta}_k\;\;and\;\;\Phi(x^{k+1},U^{k+1} )\leq F(\tilde{y}^k),\tag{29}\] which are both needed in the convergence analysis performed in Section 4.

3.4 Choice of the tolerance \({\varepsilon}_k\)↩︎

In this section, we explain how the tolerance parameter can be chosen in order to: (a) make the vector \(d^k\) a descent direction for the merit function \(\Phi\); (b) guarantee that the sufficient decrease condition 27 is satisfied for all sufficiently large values of \(m_k\). Instead of selecting a prefixed sequence \(\{{\varepsilon}_k \}_{k \in \mathbb{N}}\), we borrow the ideas from [25], [27], [29] and adopt an adaptive strategy that relates \(\epsilon_k\) with \(\tilde{\Delta}_k\). In fact, we formally set \[\nonumber {\varepsilon}_k = \frac{\sqrt{- \tilde{\Delta}_k}}{2}.\] The above setting must be interpreted in the sense of Algorithm 2, where the tolerance value is dynamically adjusted through iterations performed to construct \(\tilde{y}^k\). More precisely, Algorithm 2 consists of gradient descent steps applied to minimize the function \(G_k\) in 15 , which generate a sequence of dual vectors \(\{p^{k,j}\}_{j\in\mathbb{N}}\) and the corresponding approximate proximal gradient points \(\{\tilde{y}^{k,j} := D_\sigma(p^{k,j})\}_{j\in\mathbb{N}}\). The inner steps are then stopped at iteration \(j\) so that \[\left\| \nabla G_k({p^{k,j}}) \right\|^2 \leq - \frac{\tilde{\Delta}_{k,j}}{4} = - \frac{1}{4} \Big(\langle \nabla f_0(x^k), \tilde{y}^{k,j} - x^k \rangle + \frac{1}{2 \alpha_k } \left\| \tilde{y}^{k,j} - x^k \right\|^2 + \lambda h^*(\tilde{y}^{k,j}) - U^k \Big).\] To show that the above inequality is well posed, we first notice that \(\left\| \nabla G_k(p^{k,j}) \right\| \underset{j\to\infty}{\longrightarrow} \left\| \nabla G_k(\hat{p}^k) \right\| = 0\), while, if \(x^k\) is not a stationary point so that \(H_k(\hat{y}^k)<0\), by continuity we also have \[H_k(\tilde{y}^{k,j}) \underset{j\to\infty}{\longrightarrow} H_k (\hat{y}^k) < 0.\] Since \(\tilde{\Delta}_{k,j}\leq H_k(\tilde{y}^{k,j})\), for sufficiently large \(j\), we have \(\tilde{\Delta}_{k,j} < 0\). Therefore, the inner stopping rule is well defined and Algorithm 2 outputs \(\tilde{y}^k\) and the corresponding value \(\tilde{\Delta}_k < 0\) such that \(\|\nabla G(\tilde{p}^k)\|^2 \leq -\tilde{\Delta}_k/4\). Combining this inequality with 22 , we obtain \[\label{eq:basicHYtilde} H_k(\tilde{y}^k) - H_k(\hat{y}^k) \leq - \frac{\mu}{2} \tilde{\Delta}_k.\tag{30}\] We observe that our arguments still hold if replace the inner gradient step in Algorithm 2 with any other converging optimization method, e.g. accelerated or inertial gradient methods, since 30 depends only on the termination rule.

3.5 PnP-IPA pseudocode↩︎

We summarize the whole proposed algorithm in the following pseudocode. The main routine is PnP-IPA (which stands for Plug-and-Play Inexact Proximal Algorithm) and the inner routine is inexact_den_prox, which computes the inexact proximal update via an inner optimization method on the primal function \(h_{\sigma}\).

4 Convergence analysis↩︎

In this section, we perform the convergence analysis for Algorithm 1 by adopting a more general point of view. In particular, we refer to any composite optimization problem of the form \[\label{eq:problem} \min_{x\in\mathbb{R}^n} F(x), \;\; withF(x)=f_0(x) + f_1(x),\tag{31}\] under the following assumptions on the functions \(f_0, f_1\):

[A1]

\(f_1 : \mathbb{R}^n \to {\mathbb{R}}\cup\{+\infty\}\) is a proper, lower semicontinuous, convex function.

[A2]

\(f_0 : \mathbb{R}^n \to \mathbb{R}\) is continuously differentiable on an open set \(\Omega_0 \supset \overline{\text{dom}(f_1)}\).

[A3]

\(f_0\) has \(L_{f_0}\)–Lipschitz continuous gradient on \(\text{dom}(f_1)\).

[A4]

\(F\) is bounded from below.

Clearly, the specific settings in 7 –8 satisfy [A1]–[A4]. For nonconvex nonsmooth optimization problems, the stationary points are characterized as follows.

Any (local) minimum point of \(F\) is a stationary point, while the converse, in general, is not true, except the case where \(f_1\) is also convex. The subdifferential of the objective function of problem 31 under assumption [A2] can be expressed as stated in the following Lemma.

In particular, for convex functions, the definition of subdifferential reduces as follows.

In the sequel, instead of focusing on Algorithm 1, we will analyze the convergence from an abstract point of view, i.e., through the properties of the sequences of the iterates. Therefore, we consider any sequences \(\{x^k\}_{k\in\mathbb{N}}\), \(\{\tilde{y}^k\}_{k\in\mathbb{N}}\), \(\{\hat{y}^k\}_{k\in\mathbb{N}}\), \(\{U_k\}_{k\in\mathbb{N}}\) where \(x^{k+1} = x^k+\eta_k(\tilde{y}^k-x^k)\), with \(\eta_k\in (0,1]\) and the following holds:

[S1]

\(U^0=f_1(x^0)\) and \(U^{k+1} = \eta_kf_1(\tilde{y}^k) + (1-\eta_k) U^k\) for \(k\geq 0\);

[S2]

\(\hat{y}^k = \text{prox}_{\alpha_kf_1}(x^k-\alpha_k\nabla f_0(x^k))\) and \(H_k(\tilde{y}^k) - H_k(\hat{y}^k) \leq -\frac{\mu}{2} \tilde{\Delta}_k\) for some \(\mu > 0\), with \(\tilde{\Delta}_k\) defined as in 25 and \(\alpha_k\in [\alpha_{\min},\alpha_{\max}]\);

[S3]

\(\eta_k\) is defined as in 27 –28 , with \(\Phi\) as in 24 .

The sequences generated by Algorithm 1 clearly satisfy [S1]–[S3].

4.1 Preliminary results and well posedness of the line–search↩︎

Now we start our analysis by three preliminary lemmata, whose proofs can be found in Appendix 7.

Next, we show that the line search procedure is well-defined and that the parameter \(\eta_k\) computed at each step is bounded from below by a positive constant depending only on \(L_{f_0}, \theta, \mu\).

Besides the well posedness of the steplength selection rule, Lemma 5 has further relevant theoretical consequences, stated in the following propositions, whose proof is given in the Appendix 7.

4.2 Stationarity of the limit points↩︎

The first convergence result is stated in the following theorem.

4.3 Convergence analysis in the Kurdyka–Łojasiewicz framework↩︎

In this section we give more insight on the convergence of the sequences associated to conditions [S1]–[S3] adopting in particular the framework of the KL inequality, which is stated in the following definition.

The analysis performed in this section runs under the following assumption:

[A5]

The function \(\mathcal{F} : \mathbb{R}^n \times \mathbb{R}\to \mathbb{R}\cup\{+\infty\}\) defined as \[\label{eq:merit95func} \mathcal{F}(x, \rho) = F(x) + \frac{1}{2} \rho^2\tag{33}\] is a KL function.

Assumption [A5] is satisfied by a wide class of functions, including real-analytic and semi-algebraic mappings, and functions definable in an o-minimal structure. In particular, if \(f_0\) and \(f_1\) fall both into one of the aforementioned classes, then \(\mathcal{F}\) belongs to the same class and thus satisfies the KL property. Notice that the function \(f_1\) defined in 8 is real-analytic [15]. This assumption is crucial to ensure the convergence of the algorithm to a stationary point, as stated in the following theorem.

The boundedness of the sequence \(\{x^{k}\}_{k \in \mathbb{N}}\) is a common assumption in convergence analysis of nonconvex optimization algorithms. The assumption is naturally satisfied if the objective function \(F\) is coercive, i.e., \(F(x) \to +\infty\) as \(\left\| x \right\| \to +\infty\).

The proof of Theorem 2 is given by showing that [S1]–[S3] are a special case of the abstract scheme defined in [16], which is restated below.

We will show that the sequences associated to [S1]–[S2], with a proper setting of the surrogate function \(\Phi\) and all the auxiliary sequences, satisfies the assumptions [H1]–[H5], following the same convergence analysis adopted in [27]. Then, we will conclude the proof of Theorem 2 by applying the aforementioned abstract convergence result.

5 Numerical experiments↩︎

In this section, we evaluate the empirical performance of the proposed PnP-IPA algorithm. We focus on the deblurring inverse problem, using the 8 real-world camera shake kernels from [34] as done in previous works [35], [36], and benchmark our approach against state-of-the-art PnP methods. Results are obtained on a subset of the first 10 \(256 \times 256\) center-cropped color images of the CBSD68 dataset [37], hereafter referred to as CBSD10. In particular, our experiments are designed to demonstrate three aspects of our method:

1. the effectiveness of PnP-IPA in a high Gaussian noise regime;

2. the superiority of our algorithm over previous PnP methods in a low Gaussian noise regime, where we can leverage the flexibility of the line–search to optimally tune the regularization parameter;

3. the convergence of PnP-IPA in a deblurring scenario with Cauchy noise, which is a heavy-tailed distribution which leads to a nonconvex data-fidelity term.

With regards to the Gaussian noise, for which the data-fidelity term is convex, we compare the method to three aforementioned state-of-the-art PnP algorithms with convergence guarantees: Prox-PnP [14], relaxed-Prox-PnP and alpha-Prox-PnP [15]. Moreover, we also compare with DPIR [36], which is a widely used PnP method based on the HQS algorithm and a DRUNet denoiser that achieves excellent visual results but lacks convergence guarantees. For the Cauchy noise, we compare with the prox-PnP, for which convergence is guaranteed even in the nonconvex setting (unlike alpha-Prox-PnP), and with a Gradient Descent method plus backtracking using the Gradient-Step prior in a RED framework, as described in [38].

5.1 Implementation details↩︎

Our experiments¹ are implemented in Python using the DeepInverse [39] library, which provides a modular framework for implementing algorithms for inverse problems. For a fair comparison, the Gradient-Step denoiser and the DRUNet used in the DPIR algorithm are implemented using the neural network architectures and pre-trained weights given by the library. Regarding the Prox-PnP algorithm, its relaxed version and the alpha-Prox-PnP (as well as PnP-IPA), we use the official weights provided by the authors of [14], which can be downloaded from their GitHub repository. The parameters of the high-noise Gaussian setting are set accordingly to the guidelines provided in the respective references, while for the low-noise Gaussian setting and for the Cauchy noise an extensive grid search has been performed for all the methods to ensure the best possible PSNR. The hyperparameters configuration for the baseline methods is given in Appendix 8. All methods are stopped when the relative error \(\left\| x^k - x^{k-1} \right\| / \left\| x^{k-1} \right\|\) falls below a threshold of \(10^{-4}\). Moreover, we fixed the seed to 35. All experiments ran on a machine equipped with an Intel Core i7-12700H CPU, 16 GB of RAM, and an NVIDIA GeForce RTX 4050 (Laptop) GPU with 6 GB of VRAM.

Regarding our proposed PnP-IPA algorithm, the line-search hyperparameters are set to \(\delta = 0.5\) and \(\omega = 10^{-4}\). The step–size \(\alpha_k\), which we remind that it can be freely chosen without any theoretical constraint (apart from lower and upper bounds), is initially set to \(10^{6}\) and reduced by a factor of 3 every \(N_{\alpha}\) iterations until \(\alpha_k \leq 10^{2}\), and then it remains constant. The rationale behind this choice is given in Section 5.5. We set \(N_{\alpha} = 10\) for the Gaussian-noise experiments and \(N_{\alpha} = 25\) for the Cauchy-noise experiments, where a slower decrease of \(\alpha_k\) was found to be beneficial. The regularization parameter \(\lambda\) is tuned for each noise level and each method to achieve the best possible PSNR.

Our results have been obtained using the BenchOpt framework [40], which ensures a standardized evaluation protocol and fair comparison across different algorithms. In particular, the reported computation time is handled directly by the Benchopt library and corresponds to the wall-clock from solver start until the stopping criterion is met.

5.2 Image deblurring under high Gaussian noise↩︎

We first consider the standard image deblurring task under additive white Gaussian noise with standard deviation \(\nu = 0.05\), which is a high level of noise. We remind, as a classical result, that the data-fidelity term related to Gaussian noise is given by the standard least–squares function \(f_{data}(x) = \frac{1}{2} \left\| Ax - y \right\|^2\), which is convex. The Lipschitz constant of its gradient is \(L_{f_{data}} = \left\| A^T A \right\|\), equal to 1 for normalized kernels (as the ones we test the methods on). Regarding PnP-IPA, we set the regularization parameter \(\lambda\) to \(3\) and the denoising level \(\sigma\) to \(0.05 = 1 \nu\).

In this aggressive noise regime a high level of regularization is required to stabilize the optimization process; hence, the constraints on the regularization parameter imposed by Prox-PnP do not limit the performance of the algorithm, and all the methods perform well. In particular, as shown in Table ¿tbl:tab:high95noise95table?, Prox-PnP, alpha-Prox-PnP and PnP-IPA achieve comparable PSNR values, while DPIR performs better (but without theoretical guarantees). Our method, nevertheless, is faster than Prox-PnP and its variants.

Lastly, we remark that relaxed-Prox-PnP is not included in the comparison since, in this high noise regime, the optimal value of the relaxation parameter provided in the reference paper is exactly 1, which makes the algorithm equivalent to Prox-PnP.

In Figure 3, we plot the evolution of PSNR and relative error along time, as well as the reconstructed images from an example image provided in the CBSD68 dataset.

5.3 Image deblurring under weak Gaussian noise↩︎

We compare the same methods under weak Gaussian noise. Since in this setting we can trust more the data fidelity term, the optimal regularization parameter \(\lambda\) can be higher than in the aggressive noise case. Then, the constraints on \(\lambda\) imposed by Prox-PnP and its variants become relevant, as they force the algorithm to use a suboptimal \(\lambda\) that is too small for this noise level. In contrast, PnP-IPA is completely agnostic to constraints. This allows us to apply less regularization (tuning \(\lambda\) to its true optimal value for image quality) without compromising convergence.

PnP-IPA is run with \(\lambda = 10\), which largely exceeds the upper bound on \(\lambda\) imposed by Prox-PnP, and the denoising level \(\sigma\) is set to \(0.02 = 2 \nu\). Quantitative results are reported in Table ¿tbl:tab:low95noise95table?. We observe that PnP-IPA consistently outperforms Prox-PnP and its variants. In Figure 4 we show a visual comparison of the restored images.

5.4 Deblurring under Cauchy noise↩︎

We move beyond convex data fidelity and evaluate the algorithms on a nonconvex task: deblurring under Cauchy noise. Heavy-tailed Cauchy noise is used to model impulsive noise or outliers, and finds applications in SAR imaging [41], radar signal processing [42], underwater acoustic signal processing [43], as well as many other domains where extreme deviations from the mean are common. It is parametrized by a scale parameter \(\gamma\) that controls the tail heaviness: higher \(\gamma\) corresponds to heavier tails and more extreme outliers. Its probability density function is given by \[p(x) = \frac{1}{\pi} \frac{\gamma}{\gamma^2 + x^2}.\] In practice, the noisy image is clipped to [0,1] before reconstruction, a common practice for displayable images. The corresponding Maximum A Posteriori (MAP) estimation yields, up to additive and positive multiplicative constants that do not affect the minimization, the following data-fidelity term, which is nonconvex: \[f_{data}(x) = \frac{1}{2} \sum_{i=1}^m \log\left( \gamma^2 + (Ax - y)_i^2 \right). \label{eq:cauchy95data95fidelity}\tag{34}\] Furthermore, the Lipschitz constant of its gradient is exactly \(L_{f_{data}} = \left\| A^T A \right\| / \gamma^2\), which is extremely large for small \(\gamma\) (i.e., heavy-tailed noise). For instance, with \(\gamma = 0.01\), we have \(L_{f_{data}} = 10^4\) for normalized kernels.

Under these conditions, the theoretical bounds for Prox-PnP and its relaxed version become extremely restrictive, since the parameter \(\lambda\) should be set equal to \(1.49 \cdot 10^{-4}\) to satisfy \(\lambda L_{f_{data}} < 1.5\). Moreover, alpha-Prox-PnP cannot be applied, since it requires the data fidelity to be convex and the relaxed version of Prox-PnP does not significantly relax the constraints on \(\lambda\), which would be set in the same order of magnitude. Furthermore, DPIR is omitted from this scenario: despite an extensive hyperparameters tuning, we observed that the algorithm consistently diverges without producing any viable intermediate result. In contrast, PnP-IPA converges without any constraint.

Prox-PnP runs with the largest admissible proximal–gradient step–size \(\tau = 1.49 / L_{f_{data}} = 1.49 \cdot 10^{-4}\) and \(\sigma = 0.006\). For the Gradient Descent method, we set the initial step–size to \(10^{-3}\) and let it be handled by the backtracking procedure, while the regularization parameter is set to \(\lambda = 10^{-3}\) and the denoising level to \(\sigma = 0.015\). The initial step size is divided by a factor \(\beta =1.5\) until the descent condition is satisfied. For PnP-IPA, we set \(\lambda = 1 / 300\) and \(\sigma = 0.03\). These parameters have been tuned to achieve the best possible PSNR for each method after an extensive grid search.

Intuitively, the heavy tails of the Cauchy distribution result in a data-fidelity term that, locally, is steep enough to make the Lipschitz constant of its gradient extremely large; convergence of Prox-PnP therefore requires a tiny \(\lambda\) (i.e., a strong regularization weight \(1/\lambda\) in our splitting convention of equation ?? ). PnP-IPA, having no such constraint, can afford a regularization weight roughly one order of magnitude smaller, granting more credit to the data-fidelity term and yielding sharper reconstructions.

Table ¿tbl:tab:cauchy95table? reports the numerical results. PnP-IPA outperforms the other methods, both in terms of PSNR and computational time. In Figure 5 we report a visual comparison of the restored images, as well as the evolution of PSNR and residual along iterations.

5.5 Further remarks on PnP-IPA↩︎

In this final subsection, we provide some further remarks on the observed convergence properties of our proposed algorithm.

First, both the number of inner steps and the number of backtracking reductions \(m_k\) performed by the line search are low. More precisely, in Gaussian noise problems the inner solver performs 2 iterations on average, and always no more than 4, while for Cauchy noise problems only one inner iteration is sufficient to meet the stopping criterion. As for the backtracking loop, the number of backtracking reduction ranges from 0,1 for Gaussian low noise problems to 2 for Gaussian high and Cauchy noise. In particular, the latter information is important since a huge number of inner steps would have represented a bottleneck for the method, since every step requires an evaluation of the denoiser. This behaviour of the algorithm can be explained by the strong convexity of the inner objective function.

Secondly, our objective is to motivate the choice of the step size \(\alpha_k\), as detailed in Section 5.1. To begin with, we can rewrite the forward step of our method as \[\nonumber z^k = x^k - \alpha_k \nabla f_0(x^k) = (1 + \alpha_k \lambda) x^k - \alpha_k \nabla f_{data} (x^k).\] Hence, the inexact proximal operator of \(h_{\sigma}/ (\alpha_k \lambda)\) is computed at the point \[\nonumber \frac{z^k}{\alpha_k \lambda} = \left( 1 + \frac{1}{\alpha_k \lambda} \right) x^k - \frac{1}{\lambda} \nabla f_{data} (x^k).\] When \(\alpha_k\) is large, the term \(\frac{1}{\alpha_k \lambda}\) is small. Therefore, the proximal operator \(\mathrm{prox}_{h_{\sigma}/(\alpha_k \lambda)}(\frac{z^k}{\alpha_k \lambda})\) approaches the identity mapping, as the function \(h_{\sigma}/(\alpha_k \lambda)\) tends to zero. Consequently, we obtain \(\tilde{p}^k \approx \frac{z^k}{\alpha_k \lambda}\), where we can also neglect the term \(\frac{1}{\alpha_k \lambda} x^k\).

As a result, the update \(\tilde{y}^k = D_\sigma(\tilde{p}^k)\) of our method essentially reduces to a pure PnP proximal gradient step with step size \(\tau = 1 / \lambda\): \[\label{eq:PnP95approximation} \tilde{y}^k \approx D_\sigma \left( x^k - \frac{1}{\lambda} \nabla f_{data}(x^k) \right).\tag{35}\] Naturally, the overall algorithm is not equivalent to Prox-PnP, also due to the presence of the line-search procedure. Moreover, our method theoretically converges for any choice of \(\lambda\). We observe that for a small \(\lambda\), a larger \(\alpha_k\) is needed to reduce the term \(\frac{1}{\alpha_k\lambda}\) and obtain the approximation in 35 . Thus, even if the choice of a large initial value \(\alpha_0 = 10^6\) might appear unusual, in these specific settings the parameter \(\alpha_k\) acts less like a traditional step size and more as a parameter regulating how closely the inexact proximal step approximates a pure PnP step of the form 35 . The order of magnitude of \(\alpha_0\) in the numerical experiments was set by manual adjustment. We then let \(\alpha_k\) decrease to a fixed value, which we found beneficial to stabilize the method in its final iterations.

6 Conclusions↩︎

In this paper, we introduced PnP-IPA (Plug-and-Play Inexact Proximal Algorithm), a novel and provably convergent optimization scheme for solving imaging inverse problems using a deep neural network-based regularizer. While the Plug-and-Play framework has empirically dominated the field of computational imaging, previous approaches providing rigorous theoretical guarantees have suffered from severe limitations, such as strict upper bounds on the regularization parameter, fixed step–sizes, and the requirement of convex data-fidelity terms.

To overcome these bottlenecks, we built upon the analytical properties of the Gradient-Step (GS) denoiser and designed a novel splitting strategy which allows for the inexact proximal computation of the scaled regularizer, hence enabling the use of variable step–sizes. Furthermore, to enable an adaptive step–size selection without exact evaluations of the objective function, we introduced a novel surrogate merit function that successfully guides an Armijo-like backtracking line–search.

From a theoretical standpoint, we embedded our algorithm within the abstract KL framework, proving global convergence of the generated sequence to a stationary point of the nonconvex objective. Crucially, this result is achieved without any assumption on the regularization parameter \(\lambda\) and fully encompasses nonconvex data-fidelity terms.

Numerical experiments on image deblurring tasks confirmed the practical superiority of PnP-IPA. By lifting the theoretical constraints on \(\lambda\), our method allows for optimal parameter tuning, yielding state-of-the-art restoration quality and competitive computational times under both aggressive and weak Gaussian noise. Moreover, the robustness and flexibility of PnP-IPA were strongly highlighted in the nonconvex setting of Cauchy noise, where our algorithm significantly outperformed existing provable PnP methods that struggle with the extremely large Lipschitz constants of the data-fidelity gradient.

Future research directions will focus on proving that a single inner step is sufficient to satisfy the inexact proximal point condition, which would lead to a much more efficient algorithm. Moreover, we will explore the extension of our method to other denoisers beyond the Gradient-Step, such as Flow-Matching denoisers, thereby integrating a renoising procedure into our framework.

7 Proofs for the convergence analysis↩︎

8 Hyperparameters of the baseline algorithms↩︎

For the sake of reproducibility, this appendix collects in Table 1 the hyperparameter setting used for the three Hurault-based baselines compared against PnP-IPA in Section 5. The values mostly follow the guidelines of the original references [14], [15], except for the weak Gaussian noise for which a different combination of hyperparameters has been found to produce better results. We remind that Relaxed-Prox-PnP is omitted from the high-noise regime since its optimal \(\gamma_{\text{relax}}\) in that setting is \(1\), making the algorithm coincide with Prox-PnP.

References↩︎