Enriched-nonexpansive Plug-and-Play image restoration

1 Introduction

Image restoration seeks a clean image from degraded and partial measurements, and it underlies denoising, deblurring, super-resolution, and tomographic reconstruction. The standard formulation is variational: one minimizes the sum of a data-fidelity term and a regularizer, the former tying the estimate to the measurements and the latter encoding prior knowledge about natural images. Proximal splitting algorithms solve such problems by alternating a gradient step on the smooth fidelity with a proximal step on the regularizer. The proximity operator that realizes this second step was introduced by Moreau [1], and its use inside forward–backward splitting for signal recovery was developed by Combettes and Wajs [2]. Acceleration of these schemes has a long history: Nesterov [3] introduced the optimal first-order momentum method for smooth convex minimization, and Beck and Teboulle [4] turned this idea into the fast iterative shrinkage-thresholding algorithm that is now standard in imaging. Inertial variants that also handle nonconvex regularizers were given by Ochs et al. [5], whose iPiano method combines a forward–backward step with a heavy-ball extrapolation.

A decisive shift came when Venkatakrishnan et al. [6] observed that the proximal step of the alternating direction method of multipliers acts exactly as a Gaussian denoiser, and may therefore be replaced by any off-the-shelf denoiser, giving the Plug-and-Play paradigm. Because the substituted denoiser need not be the proximity operator of any explicit regularizer, this move gains considerable empirical power at the cost of the variational interpretation, and hence of the classical convergence guarantees. The idea was quickly broadened. Sreehari et al. [7] placed it in a consensus-equilibrium framework and identified nonexpansiveness of the denoiser as a sufficient condition for the fixed-point iteration to be well behaved. Romano et al. [8] proposed the related regularization-by-denoising framework, which builds an explicit regularizer from a denoiser, and the exact conditions under which that construction is consistent were clarified by Reehorst and Schniter [9]. Chan et al. [10] analyzed Plug-and-Play ADMM and proved fixed-point convergence under a bounded-denoiser condition together with a continuation scheme. An extensive account of the many physical and learned models that the paradigm now unifies was given in the survey of Kamilov et al. [11]. On the network side, the denoisers that make the paradigm effective are deep: Zhang et al. [12] introduced the residual denoising network DnCNN, and Zhang et al. [13] developed the DRUNet denoiser prior tailored to Plug-and-Play restoration.

Outside deblurring, the size of the admissible denoiser class becomes a practical question rather than a formal one. Ahmad et al. [14] developed Plug-and-Play reconstruction for magnetic resonance imaging, where the forward operator is a subsampled Fourier transform and the prior must supply everything the missing frequencies do not; Sun et al. [15] treated the large-scale tomographic and compressive regime through an online variant that touches only a random subset of the measurements at each step; and Wei et al. [16] learned the algorithm parameters themselves so that a single denoiser serves several tasks. A parallel line replaces the iteration by a finite unrolled network, surveyed by Monga et al. [17] and instantiated for compressive sensing by Zhang and Ghanem [18], while deep equilibrium reconstruction, introduced for imaging by Gilton et al. [19] on the equilibrium networks of Bai et al. [20], drives a learned operator to a fixed point directly. In all of these settings the operator that is iterated is trained, the measurement operator is ill-conditioned or rank deficient, and the quality of the reconstruction is governed by how expressive the prior is allowed to be. A theory that enlarges the certified class therefore acts where these applications are most constrained.

The central difficulty is that inserting an arbitrary network destroys the structure that guarantees convergence. Progress has therefore taken the form of constraining the denoiser so that the governing operator recovers a nonexpansive character. Ryu et al. [21] proved that Plug-and-Play forward–backward and ADMM converge when the denoiser residual has Lipschitz constant strictly below one, and enforced this with a spectral penalty during training. Pesquet et al. [22] constructed learned operators that are firmly nonexpansive, equivalently maximally monotone, so that the associated iteration is provably convergent, and a network architecture realizing this constraint was given by Terris et al. [23]. Hurault et al. [24] obtained convergence with a gradient-step denoiser whose residual is a conservative field, and later [25] extended this to a proximal denoiser with an unconstrained regularization parameter. Cohen et al. [26] secured convergence of regularization by denoising by projecting onto the fixed-point set of a demicontractive denoiser, one for which the distance to a fixed point decreases by an amount controlled by the residual. Each of these guarantees rests on the same lever, a Lipschitz or firm-nonexpansiveness constraint on the denoiser, and it is well documented that tightening this constraint removes exactly the sharp, high-curvature responses that make deep denoisers effective. Inertial acceleration has been imported into this setting as well: Tang and Wu [27] proposed an inertial forward–backward Plug-and-Play scheme for Poisson deconvolution, while the rate behavior of inexact nonexpansive iterations was analyzed in generality by Liang et al. [28]. What unites all of these works, and what we identify as the central research gap, is that the operator whose fixed points are sought is required to be nonexpansive or firmly nonexpansive; there is no convergence theory for Plug-and-Play that admits an expansive composite operator.

The fixed-point machinery on which the present work builds is long established. The averaging idea originates with Krasnosel'skiĭ [29], and the mean-value iteration with Mann [30]; the two-step process is due to Ishikawa [31]. Existence of fixed points of nonexpansive self-maps of bounded closed convex sets in uniformly convex spaces was established independently by Browder [32] and Göhde [33]. Weak convergence of the iterates hinges on the property isolated by Opial [34] and on the demiclosedness principle; Reich [35] proved weak convergence theorems for nonexpansive mappings in Banach spaces, and Tan and Xu [36] gave the sharp weak and strong convergence analysis of the Ishikawa process that we rely on. The uniformly convex estimate that drives asymptotic regularity is Schu's lemma [37]. On the smoothness side, the cocoercivity of the gradient of a smooth convex function, which makes the gradient step averaged, is the theorem of Baillon and Haddad [38], and the encyclopedic reference for the operator theory we use is the monograph of Bauschke and Combettes [39]. Convergence under summable perturbations is handled by the quasi-Fejérian analysis of Combettes [40], and linear rates follow from the regularity notions cataloged by Rockafellar and Wets [41] and by Dontchev and Rockafellar [42].

The operator class that removes the gap was introduced by Berinde. In [43] he defined enriched nonexpansive mappings in Hilbert spaces and showed that they are approximated by the Krasnosel'skiĭiteration, and in [44] he extended the study to Banach spaces through a retraction-displacement condition. An enriched nonexpansive map need not be nonexpansive; it becomes nonexpansive only after averaging with the identity, so the class strictly contains the nonexpansive maps. Bejenaru and Postolache [45] placed this class inside a single displacement hierarchy, averaged inside nonexpansive inside enriched, and exhibited a transfer–return dynamics that moves any operator to a nonexpansive generator and back while preserving the fixed-point set; they also introduced the inertial Mann and Ishikawa iterations adapted to displaced operators. Their study, however, is abstract: it contains no application and, by its own data-availability statement, no numerical experiments. The gap we address is therefore twofold. On the imaging side, no Plug-and-Play convergence theory admits an expansive composite operator, so the denoiser is over-constrained. On the fixed-point side, the enriched hierarchy and its transfer–return dynamics have never been connected to a concrete computational problem.

This paper closes both gaps. We ask only that the composite Plug-and-Play operator be enriched nonexpansive, a strictly weaker requirement than nonexpansiveness, and we recover convergence and acceleration through the transfer–return dynamics. The development rests on four theorems, each stated in the main text with its hypotheses written out and supported by lemmas placed where they are used or, when technical, in the appendices.

Theorem 17 performs the reduction. If the composite operator T = D_θ°G_λ is b-enriched nonexpansive, the inertial Mann and Ishikawa iterations for T coincide, term by term, with the classical Mann and Ishikawa iterations of the nonexpansive generator N=(1-ᾱ)I+ᾱT, where ᾱ=1/(b+1), and both converge strongly to a Plug-and-Play solution. The inertial correction turns out to be exactly the quantity that effects this coincidence, so the classical convergence theory applies unchanged to an expansive operator.

Whether the hypothesis of that theorem holds for a given network is decided by Theorem 22. A continuously differentiable T is b-enriched precisely when the displacement I−T is infinitesimally cocoercive with modulus 1/(2(b+1)), so the smallest admissible constant is suph−1 for a Rayleigh-type quotient h of that displacement taken over the constraint set. Two computable relaxations follow, among them the field-of-values formula b=(LT2-1)/(2(1-σ)) in terms of the Lipschitz constant L_*T* and the largest eigenvalue σ of the symmetric part of the Jacobian. The exact quotient is attained on non-scalar operators, it stays bounded as T approaches the identity, where the coarser formula degenerates, and it is evaluated by automatic differentiation.

Within the window that the certificate fixes, the step size may be chosen adaptively. Theorem 32 shows that any policy with outputs in a compact subinterval of that window preserves the conclusion of Theorem 17, whatever the policy has learned, and identifies in closed form the relaxation that minimizes the first-order residual, whose contraction is governed by the angle between the residual and its image under I − JN.

Theorem 36 quantifies the resulting behavior: the fixed-point residual decays as O(n^−1/2) with a constant proportional to b+1, the decay is linear under metric subregularity of I−N, and the conclusions survive denoiser evaluation errors of summable norm.

Around these four results we place a training objective that promotes the certified inequality and whose eigenvalue terms are differentiated through the Rayleigh–Danskin identity [46], with the spectral-normalization technique of Miyato et al. [47] and the singular-value analysis of convolutional layers by Sedghi et al. [48] supplying the practical Lipschitz control; a step-size policy trained by unrolled differentiation in the spirit of Andrychowicz et al. [49], whose greedy limit is a Barzilai–Borwein-type quantity [50]; a perturbation estimate that turns spectral measurements on a finite sample into a bound over a segment or a covered set; and an extension to arbitrary splittings that covers Plug-and-Play ADMM. Section 8 reports two experiments, one on a linear surrogate whose constants are known in closed form and one on image restoration with trained denoisers. Figure 1 summarizes the pipeline and Table 1 positions the framework against existing approaches.

Table 1. Comparison of representative Plug-and-Play methods by the denoiser condition required for convergence and by the acceleration provided.

Since the argument draws on an established body of fixed-point theory, it is worth separating the two strands at the outset. The enriched class and its defining inequality are due to Berinde [43, 44], and the displacement calculus, the mutual inversion of displacement by α and by 1/α, the invariance of the fixed-point set, the hierarchy averaged ⊂ nonexpansive ⊂ enriched, and the inertial Mann and Ishikawa iterations for displaced operators are due to Bejenaru and Postolache [45]; Section 2 assembles this material with its sources. The Plug-and-Play specialization and the apparatus that certification and control require are developed here: the identification of the composite operator as an enriched map together with the strictness statement of Proposition 15, the differential characterization of the enrichment constant and its two relaxations in Theorem 22, the near-identity and normality analysis of Propositions 23 and 25, the sampling estimate of Proposition 26, the training objective (Equation 23), the safeguarded controller and its greedy rule in Theorem 32, the rate and stability statements of Theorem 36, the splitting-agnostic extension of Theorem 40, and the experiments of Section 8.

Convergence guarantees for learned iterative solvers are usually obtained by forcing the learned map into a rigid class, nonexpansive or firmly nonexpansive, and then paying for that rigidity in reconstruction quality. The transfer–return viewpoint shows that the rigidity can be moved off the network and onto the iteration: one lets the network be as expressive as an enriched map allows, and one repairs convergence by averaging, which the inertial scheme performs automatically. Whether the same separation of concerns transfers to deep equilibrium models [19, 20] and to other settings in which a learned operator is driven to a fixed point is a question the present experiments do not settle, and Section 9 records it as such.

The remainder of the paper is organized as follows. Section 2 assembles the operator-theoretic background, including the material on monotone operators and on the Mann and Ishikawa processes that the later arguments use. Section 3 introduces the enriched hypothesis, the iteration, and Theorem 17. Section 4 develops the certificate, Theorem 22, and the training objective, and Section 5 the controller and Theorem 32. Section 6 collects the rate and stability statements in Theorem 36, and Section 7 the splitting-agnostic extension. Section 8 reports the numerical study and Section 9 concludes. The appendices contain the auxiliary lemmas, the proofs deferred from the main text, and proofs of the classical results on which the arguments rest, given for Hilbert spaces, which is the setting of every application.

2 Background and notation

Throughout, H is a real Hilbert space with inner product 〈·, ·〉 and norm ∥·∥; for the imaging application H = ℝ^*n* with the Euclidean norm. For a self-map S of a set C⊆H we write Fix(S) = {x ∈ C:Sx = x}, and I denotes the identity. A subset C is called admissible if it is nonempty, closed, convex, and bounded. We write x_*n*⇀x for weak convergence and x_*n*→x for strong convergence; in ℝ^*n* the two coincide. For α > 0 and S:C→C, the α-displacement is S_α = (1−α)I+αS, and ᾱ:=1/α is the return coefficient. Table 2 collects the recurring symbols.

Table 2. Recurring notation used throughout the paper.

2.1 Nonexpansive and monotone operators

The two languages used below, the language of nonexpansive maps and the language of monotone operators, describe the same objects through the Cayley correspondence, and both are needed: the iteration is stated in the first, the certificate of Section 4 is naturally stated in the second.

Definition 1 ([43]). Let C ⊆ H be nonempty and S:C → C. The map S is nonexpansive if ∥Sx−Sy∥ ≤ ∥x−y∥ for all x, y ∈ C; it is κ-averaged, with κ ∈ (0, 1), if S = (1−κ)I + κN for some nonexpansive N; and it is b-enriched nonexpansive, with b ∈ [0, ∞), if

∥b(x-y)+Sx-Sy∥≤(b+1)∥x-y∥, for all x,y∈C.(1)

A set-valued operator B:H⇉H is monotone if 〈u−v, x−y〉 ≥ 0 whenever u ∈ Bx and v ∈ By, and maximally monotone if its graph is not properly contained in the graph of another monotone operator. A single-valued B is ν-cocoercive, ν > 0, if

〈Bx-By,x-y〉 ≥ ν∥Bx-By∥2, for all x,y,(2)

which is monotonicity strengthened by a quantitative bound and is equivalent to νB being firmly nonexpansive. The dictionary between the two languages is the elementary equivalence

S nonexpansive⇔I-S is 12-cocoercive,(3)

proved in Lemma 49. Theorem 22 extends Equation 3 to the enriched class: S is b-enriched exactly when I−S is cocoercive with modulus 1/(2(b+1)), so the enrichment constant measures by how much the displacement of S falls short of the cocoercivity modulus that nonexpansiveness would require. That reading is what makes the certificate computable.

Lemma 2 ([39]). *Let κ ∈ (0, 1) and S:C → C. Then S is κ-averaged if and only if*

∥Sx-Sy∥2≤∥x-y∥2-1-κκ∥(I-S)x-(I-S)y∥2, for all x,y∈C.(4)

Lemma 3 ([38]). Let f(x)=12∥Ax-y∥2, so that ∇f(x) = A^⊤(Ax − y) is Lipschitz with constant L = ∥A∥². For λ ∈ (0, 2/L) the operator G_λ = I−λ∇f is λL2-averaged, and for λ ∈ (0, 2/L] it is nonexpansive. For λ ∈ (0, 1/L] the Jacobian I−λA^⊤A is symmetric with spectrum in [0, 1].

Appendix E proves Lemmas 2 and 3. The two analytic theorems that drive every convergence statement follow.

Theorem 4 ([32]). Let C *be a nonempty, bounded, closed, convex subset of a uniformly convex Banach space, in particular a Hilbert space, and let N:C → C be nonexpansive. Then Fix(N) is nonempty, closed, and convex*.

Appendix C proves Theorem 4 for Hilbert spaces; the uniformly convex Banach case is the one quoted from [32, 33].

2.2 The Mann and Ishikawa processes

Successive approximation x_*n*+1 = Nx_*n* need not converge for a nonexpansive N: the rotation of the plane by a fixed angle is nonexpansive with the origin as its only fixed point, yet every orbit stays on a circle. Convergence is restored by moving only part of the way, which is the content of the mean-value process of Mann [30],

xn+1=(1-μn)xn+μnNxn, μn∈[0,1],(5)

and of the two-step process of Ishikawa [31], in which the argument of N is itself relaxed,

yn=(1-γn)xn+γnNxn, xn+1=(1-μn)xn+μnNyn.(6)

The mechanism behind both is Fejér monotonicity: the distance to any fixed point never increases, and the amount by which it decreases at step n is proportional to μn(1-μn)∥xn-Nxn∥2. Summing this over n makes the residuals square-summable as soon as ∑nμn(1-μn)=∞, which is why every hypothesis on the step sizes below has that divergence at its center; the requirement limsup_*n*γ_*n* < 1 in Equation 6 plays the same role for the inner step. Once ∥x_*n*−Nx_*n*∥ → 0, the demiclosedness principle of Appendix C identifies every weak cluster point of {x_*n*} as a fixed point, and the property isolated by Opial [34] makes that cluster point unique. Lemma 5 and Theorems 6 and 7 record the outcome, and Appendices B to D prove them from the definitions.

Lemma 5 ([39]). Let N:C → C be nonexpansive with Fix(N) ≠ ∅, and let {x_*n*} follow Equation 5 with μ_*n* ∈ [0, 1]. Then for every p ∈ Fix(N),

∥xn+1-p∥2≤∥xn-p∥2-μn(1-μn)∥xn-Nxn∥2.(7)

In particular {∥x_*n*−p∥} is nonincreasing, {x_*n*} is bounded, and ∑nμn(1-μn)∥xn-Nxn∥2≤∥x0-p∥2.

Theorem 6 ([39]). Let C be nonempty, closed, and convex in a Hilbert space, let N:C → C be nonexpansive with Fix(N) ≠ ∅, and let {x_*n*} follow Equation 5 with μ_*n* ∈ [0, 1] and ∑nμn(1-μn)=∞. Then ∥x_*n*−Nx_*n*∥ → 0 and {x_*n*} converges weakly to a point of Fix(N); if dimH < ∞ the convergence is strong.

Theorem 7 ([36]). Let C be a bounded, closed, convex subset of a uniformly convex Banach space and N:C → C nonexpansive. Consider the Ishikawa iteration (Equation 6) with γ_*n*, μ_*n* ∈ [0, 1] satisfying ∑nγn(1-μn)<∞, ∑nμn(1-μn)=∞, and limsup_*n*γ_*n* < 1. If the norm is Fréchet differentiable, or if the space has the Opial property, that is if x_*n*⇀x implies liminf_*n*∥x_*n*−x∥ < liminf_*n*∥x_*n*−y∥ for every y ≠ x, then {x_*n*} converges weakly to a fixed point of N; if the range of C under N *is relatively compact, or if N satisfies condition (I), that is if* ∥x−Nx∥ ≥ ω(d(x, FixN)) for some nondecreasing ω with ω(0) = 0 and ω(r) > 0 for r > 0, the convergence is strong. All hypotheses hold in H = ℝ^*n*.

2.3 The enriched class and the displacement calculus

The next three statements form the calculus on which Section 3 rests: Lemma 8 is the return identity, Proposition 9 the hierarchy and Lemma 10 the invariance of the fixed-point set.

Lemma 8 ([43]). Let b ≥ 0, α: = b + 1, and ᾱ:=1/α∈(0,1]. A map E:C → C is b-enriched nonexpansive if and only if its return map N:=Eᾱ=(1-ᾱ)I+ᾱE is nonexpansive, equivalently E = N_α = (1 − α)I + αN.

Proof: Write d = x−y and w = Ex−Ey. Since ᾱ=1/(b+1) we have 1-ᾱ=b/(b+1), so

Nx-Ny=(1-ᾱ)d+ᾱw=1b+1(bd+w)=1b+1(b(x-y)+Ex-Ey).

Taking norms, ∥Nx-Ny∥=1b+1∥b(x-y)+Ex-Ey∥, so ∥Nx−Ny∥ ≤ ∥x−y∥ holds for all x, y exactly when Equation 1 holds. Solving N=(1-ᾱ)I+ᾱE for E gives E=1ᾱN-1-ᾱᾱI=(1-α)I+αN, using 1/ᾱ=α.

Proposition 9 ([45]). Let α > 0, ᾱ=1/α, and S:C → C. Then (Sᾱ)α=S and (Sα)ᾱ=S, so displacement by α and by ᾱ are mutually inverse. Consequently every averaged map, every nonexpansive map, and every enriched map is a displacement of a nonexpansive generator, and with A, N, E the averaged, nonexpansive, and enriched classes, A ⊂ N ⊂ E, both inclusions being strict.

Proof: Since αᾱ=1,

(Sᾱ)α=(1-α)I+α((1-ᾱ)I+ᾱS)=(1-αᾱ)I+αᾱS=S,

and the second identity follows by exchanging the roles of α and ᾱ. An averaged map is nonexpansive by Lemma 50, and a nonexpansive map satisfies Equation 1 with b = 0, so A⊆N⊆E. For strictness, S = −I on ℝ is nonexpansive, while S = (1−κ)I+κN with κ ∈ (0, 1) forces N =((κ−2)/κ)I, whose Lipschitz constant (2−κ)/κ exceeds one, so −I is not averaged; and Sx = −cx with c > 1 is b-enriched for b ≥ (c−1)/2 by Proposition 15 (ii) without being nonexpansive.

Lemma 10 ([45]). For every α > 0 and S:C → C, one has Fix(S_α) = Fix(S).

Proof: The equivalences x = S_αx = (1−α)x+αSx ⇔ αx = αSx ⇔ x = Sx hold because α ≠ 0.

Remark 11. The displacement coefficient α is a coordinate along the hierarchy of Proposition 9, with α < 1 averaged, α = 1 nonexpansive and α > 1 enriched, and Lemma 10 says that traveling along it fixes the fixed-point set. Certification, iteration and analysis may therefore all be carried out on the nonexpansive generator, which is the strategy adopted below, and the set of Plug-and-Play solutions is unaffected by the passage between T and N.

The background closes with the inertial iterations of the displacement calculus. Let α ≥ 1 and T:C → C satisfy Tᾱ(C)⊆C.

Definition 12 ([45]). With step-size sequences α_*n*, β_*n* ∈ [0, 1/α], the inertial Ishikawa iteration for T is

{yn=(1−βn)xn+βnTxn,xn+1=(1−αn)xn+αnTyn+(1−α)αnβn(xn−Txn),(8)

subject to ∑nβn(1α-αn)<∞, ∑nαn(1α-αn)=∞, and limsupnβn<1α. Setting β_*n* ≡ 0 gives the inertial Mann iteration x_{*n* + 1} = (1−α_*n*)x_*n* + α_*n*Tx_*n* with α_*n* ∈ [0, 1/α] and ∑nαn(1α-αn)=∞.

3 Enriched Plug-and-Play and its convergence

The variational model for recovering x ∈ ℝ^*n* from a measurement y = Ax+ε corrupted by noise ε is minxf(x)+g(x), in which f(x)=12∥Ax-y∥2 is the data fidelity and g is a regularizer that the Plug-and-Play substitution never has to make explicit, and the forward–backward Plug-and-Play iteration is x_*k*+1 = D_θ(x_*k*−λ∇f(x_*k*)) = Tx_*k* with T = D_θ°G_λ; a Plug-and-Play solution is any fixed point x^⋆ = Tx^⋆. Two hypotheses are standing throughout.

Assumption 13. H = ℝ^*n*; f(x)=12∥Ax-y∥2 with L = ∥A∥² and λ ∈ (0, 2/L), so G_λ is λL2-averaged by Lemma 3; C⊆ℝ^*n* is admissible, for instance the image hypercube [0, 1]^*n*; and D_θ(G_λ(C)) ⊆ C, so T maps C into C.

Assumption 14. There exists b ≥ 0 such that T = D_θ°G_λ is b-enriched nonexpansive on C; equivalently, with α: = b + 1 ≥ 1 and ᾱ:=1/α, the return operator

N:=Tᾱ=(1-ᾱ)I+ᾱT(9)

is nonexpansive on C.

Assumption 14 is a statement about every pair of points of C, and this is the hypothesis on which all the convergence statements below rest. In practice b is not known in advance: Section 4 produces it from spectral measurements of JT, and those measurements are taken on a finite sample. Proposition 26 supplies the estimate that closes this gap, converting a measurement on a δ-net into a bound valid on all of C at the cost of a term proportional to the Lipschitz modulus ϰ of JT; without a bound of that kind, what Section 4 certifies is Assumption 14 restricted to the sampled points, and

If D_θ is nonexpansive then, under Assumption 13, T is a composition of nonexpansive maps and hence nonexpansive, so b = 0: the classical setting is the special case b = 0 of Assumption 14. Allowing b > 0 admits denoisers whose composite operator is not nonexpansive. The next proposition shows that this relaxation is strict, and that no bound on the Lipschitz constant implies it.

Proposition 15. The following hold for T:C → C.

• (i) Every nonexpansive T is 0-enriched.
• *(ii) There are b-enriched operators with* b > 0 that are not nonexpansive; on ℝ the map Tx = −cx with c > 1 is b-enriched exactly for b ≥ (c−1)/2 and is expansive.
• (iii) An ℓ-Lipschitz operator with ℓ > 1 *need not be b-enriched for any* b ≥ 0: on ℝ² the scaled rotation T = ℓR_ϑ with cosϑ ≥ 1/ℓ satisfies no inequality Equation 1.
• (iv)T is b-enriched if and only if, writing d = x−y, w = Tx−Ty, t = ∥w∥/∥d∥ and χ = 〈d, w〉/(∥d∥∥w∥) for x ≠ y,

t2+2bχt≤2b+1,(10)

*so enrichment restricts the increment w through the pair* (t, χ), that is through its length and its alignment with d, and never through the length alone. In particular, if χ = 1 at some pair with t > 1, then Equation 10 fails for every b ≥ 0.

Part (iv) is the quantitative form of parts (ii) and (iii) and explains them: an expansive increment is admissible when it points away from d, since then χ < 0 and the term 2bχt can absorb t², and is inadmissible when it points along d, since then no amount of averaging removes the excess. One mechanism lies behind every obstruction met later, and it is the reason the certificate of Section 4 is built from the symmetric part of JT, which is the infinitesimal form of the alignment coefficient χ, and not from ∥JT∥₂, which is the infinitesimal form of t alone. The proof is in Appendix F.

Under Assumption 14, T is exactly an enriched operator N_α in the sense of Proposition 9, so Definition 12 applies. We call the resulting scheme enriched-inertial Plug-and-Play.

Definition 16. Let T = D_θ°G_λ be b-enriched with α = b + 1. The EI-PnP-Mann iteration is

xn+1=(1-αn)xn+αnDθ(xn-λ∇f(xn)), αn∈[0,1α],(11)

and the EI-PnP-Ishikawa iteration is Equation 8 with this T.

The following theorem is the first of the four announced in Section 1. Its data are those of Assumptions 13 and 14 and Definition 16, and the remaining symbols are fixed in the statement.

Theorem 17 (Reduction and convergence). Let Assumptions 13 and 14 hold with enrichment constant b ≥ 0, put α = b+1 and ᾱ=1/α, and let N=(1-ᾱ)I+ᾱT be the generator of Equation 9. Set μ_*n*: = αα_*n* and γ_*n*: = αβ_*n*, both of which lie in [0, 1] whenever α_*n*, β_*n* ∈ [0, 1/α]. Then the following hold.

• (i) Displacement identity. Tx = (1−α)x+αNx and x−Tx = α(x−Nx) for every x ∈ C, and Fix(T) = Fix(N) ≠ ∅.
• (ii) Reduction. The EI-PnP-Mann iteration (Equation 11) is, term by term, the Mann iteration x_*n*+1 = (1−μ_*n*)x_*n*+μ_*n*Nx_*n* of N; the EI-PnP-Ishikawa iteration (Equation 8) is, term by term, the Ishikawa iteration (Equation 6) of N with parameters γ_*n*, μ_*n*. The inertial correction (1−α)α_*n*β_*n*(x_*n*−Tx_*n*) in Equation 8 is exactly the term that cancels the cross term produced by the displacement, and without it the reduction fails.
• (iii) Step-size dictionary. αn(1α-αn)=μn(1-μn)/α2 and βn(1α-αn)=γn(1-μn)/α2, and limsupnβn<1α is equivalent to limsup_*n*γ_*n* < 1; hence the restrictions of Definition 12 hold if and only if those of Theorems 6 and 7 hold for μ_*n*, γ_*n*.
• (iv) Convergence. If α_*n* ∈ [0, 1/α] with ∑nαn(1α-αn)=∞, then the EI-PnP-Mann iterates satisfy ∥x_*n*−Tx_*n*∥ → 0 and converge strongly to some x^⋆ ∈ Fix(T), that is to a Plug-and-Play solution. If α_*n*, β_*n* ∈ [0, 1/α] satisfy the restrictions of Definition 12, the same conclusions hold for the EI-PnP-Ishikawa iterates.

Proof: Part (i) is Lemma 8 together with Lemma 10, and Fix(N) ≠ ∅ follows from Theorem 4, since N is nonexpansive by Assumption 14 and maps the admissible set C into itself, C being convex and T(C)⊆C. Parts (ii) and (iii) are Lemmas 54 and 55 of Appendix F. For part (iv), the Mann case follows from (ii) and (iii) and Theorem 6, which applies because ∑nμn(1-μn)=∞ and dimH < ∞; then ∥x_*n*−Tx_*n*∥ = α∥x_*n*−Nx_*n*∥ → 0 by (i). The Ishikawa case follows in the same way from Theorem 7, whose hypotheses hold in ℝ^*n*: the norm is Fréchet differentiable, C is bounded closed convex, and N(C) is relatively compact.

Corollary 18. Let C be an admissible subset of a uniformly convex Banach space, let T:C → C *be b-enriched nonexpansive with α = b+1 and generator* N=Tᾱ, and let α_*n*, β_*n* ∈ [0, 1/α] obey the restrictions of Definition 12. If N satisfies condition (I) of Theorem 7, then the EI-PnP-Ishikawa iterates converge strongly to a point of Fix(T).

Proof: The reduction of Theorem 17 (ii) and the dictionary of Theorem 17 (iii) are algebraic identities and hold in any normed space, so the iterates coincide with the Ishikawa iterates Equation 6 of N with parameters γ_*n*, μ_*n*. The set Fix(N) is nonempty by Theorem 4, and Theorem 7 gives strong convergence to a point of Fix(N) under condition (I). Lemma 10 identifies that point as an element of Fix(T).

Remark 19. Part (ii) of Theorem 17 has a geometric reading: passing from N to T = N_α stretches every increment by the factor α along the direction x−Nx, and the correction removes the resulting overshoot. The admissible window [0, 1/α] on T is the image of the classical window [0, 1] on N under contraction by ᾱ, so no step outside it can be certified.

Remark 20. The window [0, 1/α] in Equation 11 narrows as b grows, and the cost is smaller than the shrinkage suggests. By part (iii) of Theorem 17 the quantity that governs the motion is μ_*n* = αα_*n*, which ranges over the whole classical interval [0, 1] however large α is: the contraction of the window is exactly compensated by the stretching of the increments of T, since x−Tx = α(x−Nx). One iteration of Equation 11 at the largest admissible α_*n* = 1/α moves the iterate by ∥x_*n*−Nx_*n*∥, the same displacement as an unrelaxed step of the generator, for every b. The enrichment costs a constant in the residual bound: by Theorem 36 the guaranteed decay of ∥x_*n*−Tx_*n*∥ *carries the factor α = b+1, because the residual of* T *is α times the residual of N while the underlying iteration is unchanged. So a larger b inflates the measured fixed-point residual by a known factor and leaves the order of convergence and the per-iteration progress on the generator untouched. The practical limitation is elsewhere, namely that b is finite only when the certificate of Section 4 returns a finite value, which requires the aligned part of JT to stay below one; Proposition 23 describes when this is easy and Section 8 reports the values that trained denoisers attain*.

Setting b = 0 collapses T = N, the admissible window [0, 1/α] becomes [0, 1], and Theorem 17 reduces to the standard guarantee for Plug-and-Play with a nonexpansive denoiser. Corollary 18 records what survives outside a Hilbert space: the reduction is algebraic, so the conclusion extends to a uniformly convex Banach space whenever the generator satisfies condition (I).

4 Certifying and training enriched denoisers

Assumption 14 is the sole hypothesis behind Theorem 17, so a usable method must produce, for a given trained network, a number b for which the assumption holds, and must produce it from quantities a program can evaluate. This section does that. Throughout, T ∈ C¹(C; ℝ^*n*) with C convex, JT(x) is the Jacobian, JsT(x)=12(JT(x)+JT(x)⊤) its symmetric part, and

F:=I-T, JF(x)=I-JT(x), JsF(x)=I-JsT(x)(12)

the displacement of T and its Jacobian. We first record the differential criterion for nonexpansiveness, then state the characterization of the enrichment constant, then examine two regimes in which it behaves in ways the coarser bounds do not predict, and finally address the passage from a finite sample to the whole of C.

Lemma 21. Let C be convex and S ∈ C¹(C; ℝ^*n*). Then S is nonexpansive on C if and only if ∥JS(x)∥₂ ≤ 1 for all x ∈ C.

Proof: If the bound holds, then for x, y ∈ C the segment [y, x] lies in C and Sx-Sy=∫01JS(y+t(x-y))(x-y)dt, so ∥Sx−Sy∥ ≤ ∥x−y∥. Conversely, if S is nonexpansive then ∥JS(x)v∥=limh→0∥S(x+hv)-Sx∥/|h|≤∥v∥ for all v, so ∥JS(x)∥₂ ≤ 1.

A bound on ∥JT∥₂ alone cannot certify enrichment, by part (iv) of Proposition 15: at a pair where the increment is positively aligned with x−y and longer than it, no b works, and spectral normalization [47] constrains only the length. The symmetric part is the quantity that decides the matter, and Theorem 22 identifies it exactly.

Theorem 22 (Exact certificate). Let C⊆ℝ^*n* be convex, let T ∈ C¹(C; ℝ^*n*), and let F = I−T be as in Equation 12. For x ∈ C and a unit vector v set

s(x,v)=v⊤JsT(x)v, q(x,v)=∥JT(x)v∥2.

Then the following hold.

• (i) Differential criterion. For ᾱ∈(0,1] the return operator Nᾱ=(1-ᾱ)I+ᾱT is nonexpansive on C if and only if

ᾱ ≤ ψ(s(x,v),q(x,v)):=2(1-s(x,v))1+q(x,v)-2s(x,v)for all x∈C, ∥v∥=1,(13)

whenever the denominator is positive, and unconditionally where it vanishes.

• (ii) Displacement form. The two entries of Equation 13 are the displacement quantities

1+q(x,v)-2s(x,v)=∥JF(x)v∥2, 1-s(x,v)=〈v,JsF(x)v〉,(14)

so that ψ=2〈v,JsF(x)v〉/∥JF(x)v∥2.

• (iii) Characterization. T is b-enriched nonexpansive on C if and only if

〈v,JsF(x)v〉 ≥ 12(b+1)∥JF(x)v∥2 for all x∈C, ∥v∥=1,(15)

*that is, if and only if the displacement F = I−T is cocoercive with modulus 1/(2(b+1)) in the infinitesimal sense. The smallest admissible constant is*

b⋆=max{0, supx∈C, ∥v∥=1, JF(x)v≠0h(x,v)-1},h(x,v):=∥JF(x)v∥22〈v,JsF(x)v〉,(16)

*with the convention h = +∞ where the denominator is not positive. The case b = 0 of Equation 15 is the infinitesimal form of the equivalence Equation 3*.

• (iv) Computable relaxations. Put σ=supx∈Cλmax(JsT(x)), the supremum over C of the field of values of JT, that is of the set {〈v, JT(x)v〉:∥v∥ = 1}; put LT=supx∈C∥JT(x)∥2, mF=infx∈Cλmin(JsF(x))=1-σ and LF=supx∈C∥JF(x)∥2. If σ < 1 then

b⋆ ≤ bfov:=max{0,LT2-12(1-σ)},b⋆ ≤ bdis:=max{0,LF22(1-σ)-1},(17)

and both ᾱ=1/(bfov+1) and ᾱ=1/(bdis+1) are admissible in Equation 13.

Proof: (i) By Lemma 21, Nᾱ is nonexpansive on the convex set C if and only if v⊤JNᾱ(x)⊤JNᾱ(x)v≤1 for every x ∈ C and every unit v. With JNᾱ=(1-ᾱ)I+ᾱJT and ∥v∥ = 1,

v⊤JNᾱ⊤JNᾱv=(1-ᾱ)2+2ᾱ(1-ᾱ)s+ᾱ2q,

using v^⊤(JT+JT^⊤)v = 2s. Subtracting 1 and factoring,

v⊤JNᾱ⊤JNᾱv-1=ᾱ[ᾱ(1+q-2s)-2(1-s)],(18)

which for ᾱ>0 is nonpositive exactly when ᾱ(1+q-2s)≤2(1-s), that is when Equation 13 holds where 1+q−2s > 0 and unconditionally where 1+q−2s = 0.

(ii) Since JFv = v−JTv, expanding gives ∥JFv∥² = 1 − 2〈v, JTv〉+q = 1 − 2s+q, and 〈v, Js_*F*v〉 = 〈v, JFv〉 = 1−〈v, JTv〉 = 1−s.

(iii) By Lemma 8, T is b-enriched if and only if Nᾱ with ᾱ=1/(b+1) is nonexpansive, and by (i) and (ii) this is ᾱ∥JFv∥2≤2〈v,JsFv〉 for all x, v, which is Equation 15. Where JF(x)v = 0 the inequality reads 〈v, Js_*F*(x)v〉 ≥ 0 and holds automatically, since then 〈v, Js_*F*(x)v〉 = 〈v, JF(x)v〉 = 0. Solving Equation 15 for b gives b+1 ≥ h(x, v) at every pair with JF(x)v ≠ 0, whence Equation 16. Taking b = 0 in Equation 15 gives 2〈v,JsFv〉≥∥JFv∥2, the infinitesimal statement of Equation 3.

(iv) For the first bound, 1+q−2s > 0 and ψ is decreasing in q; and ∂sψ=2(1-q)/(1+q-2s)2, which is nonpositive precisely when q ≥ 1; the infimum over q is therefore at q=LT2, and since L_*T* ≥ 1 the function s↦ψ(s,LT2) is nonincreasing, so on the reachable set {s≤σ, q≤LT2, q≥s2} the infimum of ψ is at (s,q)=(σ,LT2) when L_*T* ≥ 1, giving ᾱ≤2(1-σ)/(1+LT2-2σ) and hence b≤(LT2-1)/(2(1-σ)); when L_*T* < 1 Lemma 21 already gives b = 0. The constraint q ≥ s² used here is Cauchy–Schwarz, q = ∥JTv∥² ≥ 〈v, JTv〉² = s², and it makes 1+q−2s ≥ (1−s)² ≥ 0. For the second bound, h≤supv∥JFv∥2/infv2〈v,JsFv〉≤LF2/(2mF) with m_*F* = λmin(Js_*F*) = 1−λmax(Js_*T*) ≥ 1−σ. Both relaxations dominate b^⋆ by Equation 16.

The characterization Equation 15 is the working form of the certificate. It reads enrichment as a cocoercivity property of the displacement, with the classical requirement as the single value b = 0, and it explains why the certificate is a quotient rather than a bound on one spectral quantity: the numerator of h measures how far the displacement moves a direction, the denominator how much of that motion is aligned with the direction itself, and only the ratio decides. The next two propositions describe the two regimes in which reading only the numerator, or only the denominator, gives the wrong answer.

Proposition 23 (Near-identity regime). Let T = I−εB on C with ε > 0 and B ∈ C¹(C; ℝ^*n*), and suppose there are m > 0 and β < ∞ with

λmin(JsB(x))≥m, ∥JB(x)∥2≤β for all x∈C.

Then h(x, v) ≤ εβ²/(2m) for all x, v, so

b⋆ ≤ max{0, εβ22m-1},(19)

and T is nonexpansive, that is b^⋆ = 0, as soon as ε ≤ 2m/β². In particular b^⋆ → 0 as ε → 0.

Proof: Here F = εB, so JF = εJB and Js_*F* = εJs_*B*. Hence

h=ε2∥JBv∥22ε〈v,JsBv〉=ε∥JBv∥22〈v,JsBv〉≤εβ22m,

and Equation 19 follows from Equation 16.

Remark 24. *Proposition 23 settles the behavior of the certificate for denoisers close to the identity, the regime in which a learned denoiser is most accurate. Writing the denoiser in residual form D_θ = I−R_θ, one has Js_{D_θ} = I − Js_{R_θ} and therefore λmax(Js_{D_θ}) = 1−λmin(Js_{R_θ}), so a small residual pushes σ toward 1, and the first relaxation in Equation 17 then appears to diverge. It remains bounded: as σ → 1 through a shrinking residual, the numerator* LT2-1 *tends to zero at the same time as the denominator 2(1−σ), and the ratio of the two is exactly what Proposition 23 evaluates. The limit is b^⋆ = 0, so the theory holds up where the denoiser is closest to the identity. The divergence visible in Figure 2 is a property of the relaxation b_fov along a path on which L_*T* is held fixed while σ increases, a path that a shrinking residual does not follow. A large b^⋆ comes instead from a displacement whose motion is poorly aligned with the direction it moves, which is a statement about the non-normality of JT and not about its distance from the identity; Proposition 25 makes this precise*.

Proposition 25 (Normal and symmetric operators, and attainment). Let T ∈ C¹(C; ℝ^*n*) and F = I−T.

• If JF(x) is symmetric and positive semidefinite for every x ∈ C, then b⋆=max{0,LF/2-1} with LF=supx∥JF(x)∥2; in particular T is nonexpansive precisely when L_*F* ≤ 2, and no enrichment is needed.
• In general h(x, v) = ∥JF(x)v∥/(2cosϕ), where ϕ is the angle between v and JF(x)v. Enrichment is therefore required only to the extent that the displacement turns the direction it acts on: the length factor ∥JFv∥/2 is the symmetric contribution of part (i), and the factor 1/cosϕ is what departure from symmetry adds.
• The relaxation b_fov of Equation 17 is attained, that is b⋆=bfov, whenever JT(x)≡cQ with c > 0 and Q orthogonal and cλmax(Q^s) < 1. This family contains non-scalar operators, for instance the scaled rotations T = cR_ϑ of ℝ² with ϑ ∈ (0, π) and cosϑ < 1/c, so attainment is not confined to the one-parameter scalar family Tx = −cx. For operators outside it the relaxation is in general strict; Section 8 exhibits a non-normal instance on ℝ¹²⁰ where b_fov exceeds b^⋆ by a factor close to 1.74.

Proof: (i) With JF symmetric positive semidefinite, diagonalize JF(x)=∑iςieiei⊤ with ς_*i* ≥ 0 and write v=∑iuiei with ∑iui2=1. Then h=∑iςi2ui2/(2∑iςiui2), a ratio of two convex combinations weighted by ςiui2, so h≤maxiςi/2=∥JF(x)∥2/2, with equality when v is a top eigenvector. Taking the supremum over x gives b⋆=max{0,LF/2-1}.

(ii) Writing ϕ for the angle between v and JF(x)v, one has 〈v, Js_*F*v〉 = 〈v, JFv〉 = ∥JFv∥cosϕ, so h = ∥JFv∥/(2cosϕ), which reduces to ∥JFv∥/2 exactly when the displacement moves v along v.

(iii) If JT = cQ with Q orthogonal then q = ∥cQv∥² = c² for every unit v, so L_*T* = c and q attains its maximum at every v; the maximum of s = cv^⊤Q^sv is σ = cλmax(Q^s), attained at a top eigenvector v_* of Q^s. At v_* both extremes occur together, so the infimum of ψ over (s, q) is attained on the reachable set and coincides with the value at (σ,LT2) used in Equation 17, whence b⋆=bfov. For T = cR_ϑ on ℝ² one has Q^s = cosϑI, so σ = ccosϑ and b⋆=bfov=(c2-1)/(2(1-ccosϑ)), and R_ϑ is not a multiple of the identity for ϑ ∈ (0, π).

The last question is the passage from measurement to guarantee. Both σ and L_*T*, and the quotient h itself, are suprema over C, whereas any program evaluates them at finitely many points. The following estimate closes the gap under a Lipschitz condition on the Jacobian, and is stated in the segment-wise form that the experiments of Section 8 use.

Proposition 26 (From a sample to the constraint set). Let C⊆ℝ^*n* be convex, let T ∈ C¹(C; ℝ^*n*) satisfy ∥JT(u)-JT(u′)∥2≤ϰ∥u-u′∥ for all u, u′ ∈ C, fix ᾱ∈(0,1] and write Nᾱ=(1-ᾱ)I+ᾱT.

• (i) If S⊆C is such that every point of C lies within δ of S, and

m^:=maxu∈S∥JNᾱ(u)∥2 satisfies m^+ᾱϰδ≤1,(20)

then Nᾱ is nonexpansive on C and T is b-enriched on C with b=1/ᾱ-1.

• Let x, y ∈ C and let u₀, …, u_*K* be equally spaced on the segment [y, x], with spacing δ = ∥x−y∥/K. If max0≤k≤K∥JNᾱ(uk)∥2+ᾱϰδ/2≤1, then ∥Nᾱx-Nᾱy∥≤∥x-y∥.

Proof: (i) For u ∈ C choose u′ ∈ S with ∥u−u′∥ ≤ δ. Since JNᾱ=(1-ᾱ)I+ᾱJT,

∥JNᾱ(u)∥2≤∥JNᾱ(u′)∥2+ᾱ∥JT(u)-JT(u′)∥2≤m^+ᾱϰδ≤1,

and Lemma 21 applies on the convex set C. The last assertion is Lemma 8. (ii) Every point of [y, x] lies within δ/2 of some u_*k*, so the same estimate bounds ∥JNᾱ∥2 by 1 along the segment, and Nᾱx-Nᾱy=∫01JNᾱ(y+t(x-y))(x-y)dt gives the claim.

Remark 27. Part (i) is the statement one wants, but its cost is a covering of C, whose cardinality grows exponentially with n, so it is informative for low-dimensional constraint sets and not for the image hypercube. Part (ii) is within reach of a program: Fejér monotonicity in Lemma 5 is applied to the pairs (x_n*, p) with p a fixed point, so it is enough to certify* Nᾱ along the segments joining the iterates to the limit, and each such segment is one-dimensional. Section 8 reports this segment-wise measurement for the trained denoisers, alongside the pointwise values, and the guarantees are to be read as holding on the region so certified. For a network built from smooth activations the modulus ϰ is finite and bounded by a product of layer quantities; for a network built from piecewise-affine activations JT* is piecewise constant, ϰ = 0 within each activation cell, and the pointwise measurement is exact on the cell containing the sample*.

4.1 Computing the certificate

Because f is quadratic, Gλx=x-λA⊤(Ax-y) is affine with constant Jacobian I−λA^⊤A, so by the chain rule JT(x)=JDθ(Gλx)(I-λA⊤A), and no Jacobian is ever formed. Every quantity in Theorem 22 is reached through Jacobian-vector and vector-Jacobian products, both supplied by automatic differentiation at the cost of one network evaluation each. With J = JT(x), the two relaxations need the leading eigenvalues

LT(x)2=λmax(J⊤J), σ(x)=λmax(12(J+J⊤)),LF(x)2=λmax((I-J)⊤(I-J)),(21)

computed by the power iterations v↦J^⊤Jv/∥J^⊤Jv∥ and v↦12(Jv+J⊤v), the second with a spectral shift 12(J+J⊤)+ρI, ρ > L_*T*, to make the leading eigenvalue positive.

The exact constant Equation 16 costs one further step. Maximizing h(x, ·) over the unit sphere is the generalized symmetric eigenproblem

JF(x)⊤JF(x)v=2ξJsF(x)v, hmax(x)=ξmax(x),(22)

whose pencil is symmetric with positive definite right-hand side whenever σ(x) < 1; it is solved either by a preconditioned power iteration, in which the solve against Js_*F* is performed by conjugate gradients using only the products JFv and JF^⊤v, or by direct ascent of the quotient h over the sphere. Both need only the same two automatic-differentiation primitives, and the resulting ξ_max is a lower bound on the supremum for every iterate of the ascent, so the run may be stopped at any point and the value used, with the proviso that a lower bound on h is an optimistic certificate and has to be combined with the margin of Proposition 26. In the reported experiments the ascent is run to stationarity and its value is compared with both relaxations of Equation 17.

Example 28. For Tx = −cx on ℝ with c > 1 one has JT≡−c, so σ = −c and L_*T* = c, and Equation 17 gives ᾱ=2(1+c)/(1+c2+2c)=2/(1+c) and bfov=(c2-1)/(2(1+c))=(c-1)/2, the constant found by hand in Proposition 15. Here JF* = 1+c and h = (1+c)²/(2(1+c)) = (1+c)/2, so b^⋆ = (c−1)/2 as well: the two agree, as Proposition 25 (iii) predicts, since JT is c times an orthogonal map of ℝ¹. Figure 2 plots the certified b_fov as a function of L_*T* for several values of σ*.

Example 29. Let T = ℓR_ϑ on ℝ² with ℓ = 1.2 and ϑ=π2, a scaled quarter-turn. Then JT = ℓR_ϑ has L_*T* = ℓ = 1.2 > 1, so T is expansive, while its symmetric part is ℓ2(Rϑ+Rϑ⊤)=ℓcosϑI=0, so σ = 0. Theorem 22 (iv) gives ᾱ=2/(1+ℓ2)≈0.82 and bfov=(ℓ2-1)/2=0.22, and the return operator satisfies ∥JN∥2=(1-ᾱ)2+ᾱ2ℓ2=1, *so the certificate is attained, in agreement with Proposition 25 (iii). A pure rotation with ℓ = 1 needs no enrichment, whereas the scaled rotation does, because scaling introduces the expansion that averaging must absorb, even though the symmetric part vanishes and no argument based on L_*T* alone distinguishes the two*.

4.2 A training objective that promotes the certificate

Remark 30. *Spectral normalization [47] caps the Lipschitz constant of each layer and hence bounds L_*T* through the chain rule, but by Proposition 15 (iv) it does not bound σ, so it cannot by itself certify enrichment. Controlling σ means limiting the aligned expansion of the composite Jacobian. Writing D_θ = I−R_θ, one has Js_{D_θ} = I − Js_{R_θ}, so λmax(Js_{D_θ}) = 1−λmin(Js_{R_θ}), and a nonnegativity constraint on the symmetric part of the residual Jacobian yields λmax(Js_{D_θ}) ≤ 1. The penalty Equation 23 enforces the composite analog of this condition directly, which is why it acts on the composite operator rather than on the layers of the network.*

Let Lden(θ)=E(u,w)∥Dθ(u+w)-u∥2 be the Gaussian-denoising loss and augment it with one-sided spectral penalties driving (σ, L_*T*) toward a target (σ_*, L_*) with σ_* < 1,

L(θ)=Lden(θ)+η1Ex[(λmax(JsT(x))-σ)+2]+η2Ex[(∥JT(x)∥2-L)+2].(23)

The penalties are one-sided and act on samples, so Equation 23 promotes the certified inequality rather than imposing it; the value of b that enters the algorithm is the one measured after training by the procedure of Section 4.1, together with the margin of Proposition 26, and never the target (σ_*, L_*) itself. When several forward operators share one denoiser, the expectation in Equation 23 is taken over the operators as well, and the composite is recertified for each of them. The gradients of the eigenvalue terms follow from first-order perturbation theory.

Lemma 31 ([46]). Let M(θ) be a symmetric matrix depending smoothly on θ with a simple largest eigenvalue λmax and unit eigenvector v_*. Then λmax is differentiable and ∇θλmax(M(θ))=v⊤(∂θM(θ))v, obtained by differentiating the Rayleigh quotient v⊤M(θ)v with v_* held fixed.

Proof: Differentiating λmax=v⊤Mv and using Mv_* = λmaxv_* gives ∇θλmax=2λmaxv⊤∂θv+v⊤(∂θM)v, and differentiating v⊤v=1 shows v⊤∂θv=0, so the first term vanishes.

Because each penalty needs only the leading eigenpair, a few power-iteration steps per minibatch suffice, and Lemma 31 lets the resulting scalar flow through ordinary back-propagation.

5 Safeguarded step-size control

Theorem 17 holds for every admissible step-size sequence, so the steps may be chosen adaptively provided the choice never leaves the admissible set. We work with the Mann form through the dictionary μ_*n* = αα_*n* ∈ [0, 1] of Theorem 17 (iii); the residual map is r(x) = x−Nx with Jacobian M(x) = I − JN(x). The next theorem is the third of the four announced in Section 1, and it separates what the controller may change from what it may not.

Theorem 32 (Safeguarded control). Let Assumptions 13 and 14 hold with α = b+1 and ᾱ=1/α, fix a∈(0,12), and let πθ:ℝd→[0,1] be measurable. Define

μn=a+(1-2a)πθ(zn)∈[a,1-a],αn=ᾱμn∈[0,1/α],(24)

where z_*n* is any feature computed from the history up to x_*n*. Then the following hold.

• (i) Guarantee. For every realization of π_θ, and whatever θ has been trained to do, the EI-PnP-Mann iterates Equation 11 with steps Equation 24 converge strongly to a Plug-and-Play solution and ∥x_*n*−Tx_*n*∥ → 0. The value of θ affects the speed and never the conclusion.
• (ii) Greedy step. Let J = JN(x_*n*), M = I−J, r_*n* = x_*n*−Nx_*n* ≠ 0 and Mr_*n* ≠ 0. Under the first-order model Nx≈Nx_*n*+J(x−x_*n*) the update x_*n*+1 = x_*n*−μr_*n* produces the residual r_*n*+1 = (I−μM)r_*n*, and the relaxation minimizing ∥r_*n*+1∥ over μ ∈ ℝ is

μn⋆=〈rn,Mrn〉∥Mrn∥2,∥rn+1∥∥rn∥(μn⋆=1-〈rn,Mrn〉2∥rn∥2∥Mrn∥2=|sin∠(rn,Mrn)|.(25)

One has μn⋆≥0, and the clamped value μn=min{max{μn⋆,a},1-a} satisfies Equation 24 and is therefore covered by (i). To first order, ∥rn+1(μn⋆)∥≤∥rn+1(μ)∥ for every μ ∈ [a, 1−a], with strict inequality unless μ=μn⋆.

• (iii) Linear contraction on affine generators. If Nx = Cx+e with M = I−C and λmax(Js_*C*) =:σ_*N* < 1, so that λmin(Js_*M*) = 1−σ_*N* > 0, the greedy iteration contracts the residual at every step,

∥rn+1∥∥rn∥ ≤ 1-(1-σN∥M∥2)2 < 1,(26)

so it converges R-linearly with this factor, meaning that ∥r_n∥ is dominated by a geometric sequence with that ratios.

Proof: (i) By construction μ_*n* ∈ [a, 1−a], so μ_*n*(1−μ_*n*) ≥ a(1−a) > 0 and ∑nμn(1-μn)=∞; Theorem 17 (ii) identifies the iteration with the Mann iteration of N with these μ_*n*, and Theorem 17 (iv) gives the conclusion, together with ∥x_*n*−Tx_*n*∥ = α∥x_*n*−Nx_*n*∥ → 0.

(ii) From x_*n*+1 = x_*n*−μr_*n* and the first-order model, Nx_*n*+1≈Nx_*n*−μJr_*n*, so r_*n*+1 = r_*n*−μr_*n*+μJr_*n* = (I−μM)r_*n* and ∥rn+1∥2=∥rn∥2-2μ〈rn,Mrn〉+μ2∥Mrn∥2, a strictly convex quadratic in μ with the stated minimizer and minimum value ∥rn∥2-〈rn,Mrn〉2/∥Mrn∥2; dividing by ∥rn∥2 gives the sine expression. Since ∥J∥₂ ≤ 1 by Lemma 21, 〈rn,Mrn〉=∥rn∥2-〈rn,Jrn〉≥0. The comparison with a fixed μ follows because projection onto [a, 1−a] preserves the ordering of the values of a convex function around its minimizer.

(iii) By (ii) the per-step factor is 1-cos2∠(rn,Mrn), and 〈rn,Mrn〉=〈rn,JsMrn〉≥(1-σN)∥rn∥2 while ∥Mr_*n*∥ ≤ ∥M∥₂∥r_*n*∥, so cos∠(r_*n*, Mr_*n*) ≥ (1−σ_*N*)/∥M∥₂, which gives Equation 26.

Equation 25 is a Barzilai–Borwein-type step [50] for the fixed-point residual, and the attainable one-step gain is the sine of the angle between the residual and its image under M. Since Mr_*n* = r_*n*−JN(x_*n*)r_*n* is computed by one Jacobian-vector product, or by the finite difference Mr_*n*≈r_*n*−(N(x_*n*+ϵr_*n*)−Nx_*n*)/ϵ, the greedy step costs one extra evaluation of the denoiser per iteration.

Remark 33. The recursion r_*n*+1 = (I−μ_*n*M)r_*n* identifies the greedy iteration with a Richardson method for the linear system Mr = 0, and Equation 25 is its exact steepest-descent step for the quadratic energy 12∥r∥2. What the greedy rule cannot do is look ahead: it optimizes the next residual only, and on problems where the residual must first move into a slowly contracting subspace this is not the best available choice. That is the gap an amortized policy is meant to close, and it is the reason the feature vector below is built from the same inner products that appear in Equation 25.

5.1 An amortized policy

A learned controller π_θ can amortize multi-step look-ahead by minimizing the horizon-K cost J(θ)=E[∑n=1Kωnφ(xn)], with φ(x) = ∥x−Nx∥² and weights ω_*n* > 0, over problem instances. Since the unrolled solver is differentiable, its gradient is given by the adjoint recursion, proved in Appendix F.

Proposition 34. Write Φ_*n*(x) = x−μ_*n*r(x) with μ_*n* = a+(1 − 2a)π_θ(z(x)), so x_*n*+1 = Φ_*n*(x_*n*). Let p_*K* = ω_*K*∇φ(x_*K*) and run, for n = K−1, …, 1,

pn=ωn∇φ(xn)+(∂xΦn(xn))⊤pn+1,∂xΦn(xn)=I-μnM(xn)-r(xn)(∇xμn)⊤,(27)

with ∇xμn=(1-2a)Jz(xn)⊤∇πθ(zn), where Jz is the Jacobian of the feature map and ∇π_θ the gradient of the policy in its argument. Then

∇θJ=-∑n=1K-1〈pn+1,r(xn)〉(1-2a)∇θπθ(zn).(28)

Every quantity in Equations 27 and 28 is a Jacobian-vector or vector-Jacobian product with N, so the controller is trained by truncated back-propagation through the solver in the manner of [49]. The feature vector used in Section 8 is

zn=(log∥rn∥, μn⋆, ∥Mrn∥/∥rn∥, 〈rn,Mrn〉/∥rn∥2),(29)

whose last three entries are the statistics of Theorem 32 (ii), and π_θ is a multilayer perceptron with two hidden layers of thirty-two units and hyperbolic-tangent activations, followed by a sigmoid output, so that Equation 24 holds by construction. In Section 8.6 it is trained for ninety Adam updates at learning rate 3 × 10⁻³ on a horizon of K = 8 steps with weights ωn=0.85K-1-n, the denoiser Jacobian being held fixed within each unrolled window, which is the usual approximation for learned optimizers. The horizon objective also admits a Markov-decision-process reading with state z_*n*, action μ_*n* ∈ [a, 1−a], and reward R_*n* = log(∥r_*n*∥/∥r_*n*+1∥), whose undiscounted return is the total log-reduction of the residual; a policy-gradient update on a stochastic policy supported on [a, 1−a] optimizes it while preserving convergence for every trajectory.

Remark 35. Theorem 32 separates two layers of the design. The certificate of Section 4 fixes the admissible window [0, 1/α] from the enrichment constant, and the controller of the present section learns, within that fixed window, a step-size policy that accelerates convergence. The two layers are independent: retraining the controller cannot enlarge the window, and recertifying the denoiser cannot violate the policy's guarantee, so the safety of the method is decided by the certificate and its speed by the controller.

6 Rates and stability

By Theorem 17 (ii) the EI-PnP-Mann iteration is Mann on the generator N, and ∥x_*n*−Tx_*n*∥ = α∥x_*n*−Nx_*n*∥, so every rate for N transfers to T with the explicit factor α = b+1. The fourth theorem collects the resulting statements together with the behavior under inexact evaluation of the denoiser, which is the situation in practice.

Theorem 36 (Rates and stability). Let Assumptions 13 and 14 hold with α = b+1, let N be the generator of Equation 9, and let {x_*n*} follow the EI-PnP-Mann iteration with steps μ_*n* = αα_*n* on the generator. Fix p ∈ Fix(T).

• (i) Sublinear residual rate. If μ_*n*≡μ ∈ (0, 1), then for all n ≥ 0

∥xn-Txn∥ ≤ α∥x0-p∥μ(1-μ)(n+1),1n+1∑k=0n∥xk-Txk∥2 ≤ α2∥x0-p∥2μ(1-μ)(n+1).(30)

• (ii) Linear rate under metric subregularity. If μ_*n* ∈ [a, 1−a] with a∈(0,12) and I−N is metrically subregular on the Fejér-invariant ball B=B¯(p0,∥x0-p0∥) for some p₀ ∈ Fix(N), that is d(x, FixN) ≤ τ∥x−Nx∥ on B, then with ρ:=max{0,1-a(1-a)τ-2}<1 the distances d_*n*: = d(x_*n*, FixN) obey d_*n*+1 ≤ ρd_*n*, and

∥xn-x⋆∥≤2ρn1-ρd0, ∥xn-Txn∥≤2αρnd0,(31)

for a solution x^⋆ ∈ Fix(T).

• (iii) Stability under inexact evaluation. If instead x_*n*+1 = (1−μ_*n*)x_*n*+μ_*n*(Nx_*n*+e_*n*) with μ_*n* ∈ [a, 1−a] and ∑n∥en∥<∞, then {x_*n*} is bounded, ∥x_*n*−Tx_*n*∥ → 0, and x_*n* converges strongly to a fixed point of T. If moreover the hypotheses of (ii) hold and ∥en∥≤Cρ̄n with ρ̄∈[0,1), the convergence is R-linear with rate max{ρ,ρ̄}.
• (iv) Optimal constant relaxation for affine generators. If Nx = Cx+e with C symmetric and spectrum in [m₋, m₊], −1 < m₋ ≤ m₊ < 1, and μ_*n*≡μ ∈ (0, 1], then ∥xn+1-x⋆∥≤ϱ(μ)∥xn-x⋆∥ with x^⋆ = (I−C)⁻¹e and

ϱ(μ)=max{|1-μ(1-m+)|, |1-μ(1-m-)|}.(32)

Over all μ > 0 this is minimized at μ†=2/(2-m+-m-), with value (m₊−m₋)/(2−m₊−m₋) < 1. The admissible range is (0, 1], and μ^† ≤ 1 exactly when m₊+m₋ ≤ 0; when m₊+m₋ > 0 the minimum over (0, 1] is at the endpoint μ = 1, with value max{|m₊|, |m₋|}. In both cases the minimizer over the admissible range is min{μ^†, 1}.

Figure 3 plots the factor Equation 32 of part (iv) against the relaxation for three symmetric spectra. The proof is assembled in Appendix F from Lemmas 56, 45 and the argument of Lemma 5; part (iii) is a quasi-Fejérian argument, in which the distance to a fixed point grows by at most a summable amount at each step, in the style of [40].

Remark 37. *The factor α = b+1 in Equation 30 is the only trace the enrichment leaves on the rate, in agreement with Remark 20. At b = 0 the bound reduces to the classical residual rate of the Krasnoselskii–Mann iteration, so the enriched theory degrades smoothly to the nonexpansive one as the denoiser is tightened toward the boundary of the admissible class*.

Remark 38. *The subregularity hypothesis in part (ii) is not restrictive for restoration. It holds automatically when N is affine, as in the surrogate of Section 8, and more generally at any isolated Plug-and-Play solution, where the residual map has a nonsingular generalized derivative in the sense cataloged in [41, 42]. When the fidelity term is strongly convex, for instance when A is injective, the gradient step is already a contraction and part (ii) applies with an explicit modulus, giving a linear rate with no regularity assumption on the denoiser beyond enrichment*.

Remark 39. *Part (iii) tolerates any error sequence of summable norm, which covers finite-precision evaluation and early-stopped inner solvers. If the denoiser is itself defined by an inner optimization, a convenient rule is to solve the n-th inner problem to accuracy* ∥en∥≤C/n2; the resulting errors are summable, so convergence is preserved while the inner tolerance is allowed to relax as the outer iteration proceeds, which saves computation in the early iterations where high inner accuracy is wasted.

7 A general splitting perspective

The forward–backward operator is one instance of a broader principle: the framework applies to any splitting whose Plug-and-Play operator is enriched.

Theorem 40. Let C⊆ℝ^*n* be admissible and let S:C → C be any operator whose fixed points are the Plug-and-Play solutions of a given splitting. Suppose S is b-enriched nonexpansive, so that NS=(1-ᾱ)I+ᾱS is nonexpansive. Then Fix(S) ≠ ∅, *and both the inertial Mann and Ishikawa iterations for S converge to a point of* Fix(S); *the certificate of Theorem 22 applies verbatim to S through JS*.

Proof: Every step of Section 3 used only that T is enriched with T(C)⊆C and C admissible, not the specific form D_θ°G_λ. Replacing T by S and N by N_*S* throughout Theorem 17 yields the claim; existence follows from Theorem 4, and Theorem 22 depends only on the Jacobian of the iterated operator.

Example 41. Plug-and-Play ADMM iterates vk+1=proxf/ρ(zk-uk), zk+1=Dθ(vk+1+uk), u^*k*+1 = u^*k*+v^*k*+1−z^*k*+1. Eliminating v and z *presents this as a Douglas–Rachford iteration on the splitting variable w = z−u, driven by an operator S_DR built from the reflections of prox_*f*/ρ and D_θ, with a relaxation parameter multiplying the reflection; the correspondence between the fixed points of that operator and the solutions of the scheme is the standard one for Douglas–Rachford splitting [39, Ch. 26], and it is unaffected by the substitution of D_θ for a proximity operator, since only the fixed-point relation is used. When D_θ is only enriched, one certifies S_DR directly through Theorem 22 and applies Theorem 40, so the over-relaxed inertial Plug-and-Play ADMM inherits the same guarantee and admissible-window structure, with the reflection relaxation playing the role of α_*n*. The firm nonexpansiveness of the proximity operator that makes S_DR well defined is Lemma 52*.

Remark 42. Theorem 40 covers, beyond forward–backward and ADMM, the three-operator splitting used when the fidelity decomposes into a smooth and a nonsmooth part, and the primal–dual schemes used when the forward operator is not readily inverted. In each case the object to certify is the fixed-point operator of the scheme, whose Jacobian is assembled from those of the constituent resolvents and the linear maps by the chain rule, so Theorem 22 applies without modification, and the admissible window is read off the certified displacement exactly as in the forward–backward case.

Table 1 positions the framework against representative Plug-and-Play methods. The distinguishing feature is the admissible denoiser class: existing convergent methods require the relevant operator to be nonexpansive or firmly nonexpansive.

The second column of Table 1 places the present work. Two of the entries buy convergence by tightening the denoiser itself: the residual-Lipschitz condition of [21], the sharpest of the strictly nonexpansive guarantees, constrains the network globally and, as its authors observe, degrades denoising away from the training noise level, while [22] confines the operator to the half-averaged regime, the most conservative point of the averaged scale. Two more buy it by fixing the form of the residual, which [24] and [25] require to be a conservative field, recovering an explicit though nonconvex energy in exchange. Elsewhere the price is paid by the algorithm rather than the operator: [10] needs a continuation scheme that strengthens the denoiser as the iteration proceeds, and [26], which does admit the wider demicontractive class, replaces the iteration by a projection onto the fixed-point set. The remaining two entries constrain the setting instead, [27] to Poisson noise and [8] to denoisers of symmetric Jacobian, a condition rarely met by deep networks [9]. Enrichment of the composite, strictly weaker than nonexpansiveness by Proposition 15, asks less than any of these, and it is compatible with both inertial and learned acceleration inside one admissible-window theory. What it costs is the shrinkage of that window from [0, 1] to [0, 1/α] described in Remark 20, which the controller of Section 5 is designed to absorb.

8 Numerical study

Two experiments are reported. The first is a linear surrogate on ℝ¹²⁰ for which every constant of Theorem 22 is available in closed form, so that the certified values can be checked against ground truth and the admissible window against the true stability boundary. The second is image restoration with trained denoisers, which asks what the enlarged operator class contains in practice and what it costs in reconstruction quality at an unchanged guarantee. Everything reported is produced by one script with fixed seeds.

8.1 A certified linear surrogate

We build a composite operator T(x) = Wx+c on ℝ¹²⁰ that is enriched nonexpansive and not nonexpansive. The matrix W is assembled from two-by-two upper-triangular blocks with double eigenvalue d_*i* on an equally spaced grid in [0.4730, 0.6883] and off-diagonal t_*i* = 2(0.92−d_*i*), placed in a fixed orthonormal basis; each block has spectral radius d_*i* < 1 and symmetric-part maximal eigenvalue d_*i*+t_*i*/2 = 0.92, while the blocks with the smallest d_*i* have largest singular value above one, so W is non-normal with ∥W∥₂ = 1.0978 > 1. This is the situation isolated by Proposition 15 (iv). Table 3 reports the measured quantities; the closed-form block analysis is Appendix G.

Table 3. Certification of the linear surrogate via Theorem 22.

Although ||W||₂ > 1, so that the operator is not nonexpansive, the certified averaging yields a return operator with norm at most one.

The decisive lines of Table 3 are the last two. The relaxation predicts, and direct computation confirms, that averaging the expansive operator W, with ∥W∥₂ = 1.098, by ᾱ=0.438 yields a nonexpansive generator with norm 0.990. Solving the generalized eigenproblem (Equation 22) for the same operator returns b^⋆ = 0.736 against b_fov = 1.282, a factor of 1.74, and the resulting return operator has norm exactly 1.0000, so the exact certificate is attained here while both relaxations are strictly conservative. For a non-normal operator on ℝ¹²⁰ this confirms what Proposition 25 (iii) establishes in closed form for scaled orthogonal maps, and it widens the admissible window from (0, 0.438] to (0, 0.576], an increase of 31%. The constant b = 1.28 is in any case inaccessible to a method that controls only the Lipschitz constant.

With the operator certified, we compare iterations toward the fixed point x^⋆ = (I−W)⁻¹c: Picard x_*k*+1 = Tx_*k*; EI-PnP-Mann with constant relaxation μ = 0.5 on the generator; and EI-PnP-Mann with the greedy relaxation (Equation 25) clamped to the admissible window, each run with the window from b_fov and with the window from b^⋆. Figure 4 shows the trajectories and Table 4 the iteration counts to ∥xk-x⋆∥≤10-10.

Table 4. Iterations to reach ||xk-x⋆||≤10-10 on the certified surrogate.

The second column gives the relaxation applied to T, written s for the coefficient in x_*k*+1 = (1−s)x_*k*+sTx_*k* and α_*n* for the admissible steps of Equation 11. The greedy step roughly halves the guaranteed constant-step count, and the window of Theorem 22 (iii) removes a further quarter.

The safeguarded iterations converge exactly as Theorem 17 predicts, and the greedy rule reduces the constant-step count by 52% with the relaxed window and by 52% again with the exact window, confirming Theorem 32 (ii). The exact certificate is worth about a quarter of the remaining iterations, 173 against 129, because a wider window admits a larger clamped step. The position of the Picard iteration on this instance is particular to it. Picard reaches the tolerance in 68 iterations, faster than any safeguarded run, and this is a property of the instance and not of the method: the composite here happens to have spectral radius 0.688 < 1, so successive approximation converges, while ∥W∥₂ > 1 means that no nonexpansive theory certifies it. For a nonlinear enriched denoiser the unrelaxed iteration can diverge, and the safeguarded schemes cannot. The acceleration reported above is therefore a comparison among safeguarded schemes on an instance on which the unguaranteed iteration is the fastest. Figure 5 plots the per-iteration contraction factor of the safeguarded runs. The greedy factor lies below the constant one over the part of the run in which the error is still decreasing geometrically, and rises once the iterate is within rounding distance of the fixed point, where the ratio ceases to be informative. The greedy step is stable across the run, its clamped value remaining at the upper edge of the admissible window, so the acceleration comes from consistently using the largest certified step and not from delicate tuning; the per-step bound of Theorem 32 (iii) is far from tight on this instance, since (1−σ_*N*)/∥I − JN∥₂ = 0.070 there and the guaranteed factor 1-0.0702 accounts for only a small part of the observed decrease; the bound is a worst case over residual directions, and the greedy rule does better because the residual aligns with the slow subspace.

Figure 6 sweeps the relaxation s applied to T and plots the spectral radius of the resulting affine iteration matrix (1−s)I+sW; convergence holds where this radius is below one. The eigenvalues of W are real and lie in (0, 1), so the exact stability boundary is scrit=2/(1-minidi)=3.80. The certified windows, (0, 0.438] and (0, 0.576], lie well inside it. The step-size restrictions of Definition 12 therefore never certify a divergent relaxation, and, as is usual for worst-case guarantees, they are conservative, by a factor of 8.7 for the relaxation and 6.6 for the exact certificate. The controller of Section 5 operates inside that margin.

Because the surrogate is affine, the metric subregularity hypothesis of Theorem 36 (ii) holds globally, which is why every method converges R-linearly in Figure 4. The surrogate tests the operator theory, the certificate, the reduction, the window and the greedy step under conditions where every constant is known exactly, and it says nothing about restoration quality, which is the subject of the rest of this section.

8.2 Setting of the restoration study

The denoisers are residual networks of the DnCNN type [12], seven convolutional layers of forty-eight channels with smooth activations, about 1.05 × 10⁵ parameters, written as D_θ = I−R_θ. A piecewise-affine network would make JT piecewise constant and the pointwise certificate exact on each activation cell, as noted in Remark 27; we use a smooth activation instead, so that the Jacobian is Lipschitz and Proposition 26 applies as stated. Training uses 44 × 44 patches drawn from the four hundred training images of the Berkeley segmentation database [51] in the split distributed with [12], additive Gaussian noise of level 10/255, three thousand Adam [52] updates with batch sixteen at learning rate 10⁻³, reduced by a factor 0.2 after 60% and 85% of the schedule, and the objective Equation 23 with weights η₁ = η₂ = 60 and targets (σ_*, L_*) that differ between models and are stated with the results. The activation is the exponential linear unit, which is C¹ so that Proposition 26 applies. The spectral terms of Equation 23 are evaluated on a sub-batch of four patches of side forty by warm-started power iterations and differentiated through Lemma 31. Because one denoiser serves three forward models, the expectation in Equation 23 is taken over the three operators, drawn in turn across updates.

Three inverse problems are used, all on gray-scale images in [0, 1] with additive Gaussian measurement noise of standard deviation 0.01: Gaussian deblurring with a 9 × 9 kernel of width 1.6; motion deblurring with a 9 × 9 diagonal streak; and single-image super-resolution by a factor of two, with a Gaussian antialiasing filter of width 1.0 before decimation. In every case ∥A∥ ≤ 1, so L = ∥A∥² ≤ 1, and the gradient step size is selected for each operator from {1.2, 1.6, 1.9}⊂(0, 2/L) on two held-out images; every selection returned λ = 1.9, and Lemma 3 keeps G_λ nonexpansive at that value. The safeguard floor of Equation 24 is a = 10⁻³. The iteration is started from the observation, upsampled bicubically in the super-resolution case, and stopped when the relative fixed-point residual falls below 10⁻⁵ or after eighty iterations; the relaxation sweep of Section 8.7 uses a budget of one hundred iterations instead. Reconstruction is measured by peak signal-to-noise ratio and by structural similarity [53] on four test images of size 128 × 128 from the classical set distributed with [12]; the certificates are measured on two held-out crops of size 64 × 64 from the sixty-eight test images of [51], which are disjoint from the training split.

Three quantities are reported for each operator: the two relaxations b_fov and b_dis of Theorem 22 (iv), and the constant b^⋆ of Theorem 22 (iii). All spectral quantities are computed by a Lanczos iteration on the corresponding symmetric operator, whose matrix-vector products are the automatic-differentiation primitives of Section 4.1; the leading eigenvalues of a deep composite are clustered, and a plain power iteration is unreliable there. The constant b^⋆ is located by bisection on ᾱ, with the criterion of Theorem 22 (i) checked at each trial value, and the returned ᾱ is verified by measuring ∥JN∥₂ at that value. The bisection has resolution 2⁻⁵, so the reported b^⋆ is an admissible constant and an upper bound for the smallest one, and every comparison with the relaxations below is conservative in that direction. Each quantity is a maximum over the sample, and by Proposition 26 it certifies the whole of C only up to the perturbation term of that statement; Section 8.7 reports the segment-wise measurement.

8.3 Certified constants of trained denoisers

Four denoisers were trained, identical except for the targets in Equation 23: one with no spectral penalty, one with (σ_*, L_*) = (0.85, 0.90), one with (0.85, 1.15), and one in which the Lipschitz term is dropped and only the symmetric part is constrained, with σ_* = 0.90. Table 5 reports the constants of the composite operator for each of them; the denoiser with the looser Lipschitz target was evaluated on Gaussian deblurring, and the remaining three on all three forward models.

Table 5. Certified constants of the composite operator for denoisers of identical architecture and identical training budget, differing only in the targets of Equation 23.

Measurements are on 2 held-out crops. The unconstrained model has σ > 1, so no finite constant exists for it and it carries no guarantee. The last column is the spectral norm of the Jacobian of the return operator at the certified ᾱ=1/(b⋆+1), which confirms the criterion of Theorem 22 (i) directly.

The unconstrained denoiser is outside the theory. Its composite has σ between 1.29 and 1.48 on every task, and by Theorem 22 (iii) no finite b exists once σ exceeds one, so the operator is not enriched and no relaxation of Equation 11 converges by Theorem 17. The penalty Equation 23 brings σ below one, and it is the symmetric part and not the Lipschitz constant that has to be moved: the fourth denoiser, whose objective contains no Lipschitz term at all, is certified on every task.

The enriched class is occupied. On Gaussian deblurring the denoiser trained with L_* = 1.15 produces a composite with L_*T* = 2.08, so the operator is expansive and satisfies none of the nonexpansive, firmly nonexpansive or contractive-residual conditions of Table 1, while its symmetric-part maximum is 0.823 and Theorem 22 certifies it with b^⋆ = 0.600, giving the admissible window (0, 0.625]. The last column of Table 5 verifies the criterion of Theorem 22 (i) directly: the return operator at the certified ᾱ has Jacobian norm 0.952 on the sample, safely below one.

The certified constant is far below both relaxations. For that same operator b_fov = 9.39 and b_dis = 25.49, against b^⋆ = 0.600, a ratio of 15.6 between the field-of-values constant and the certified one, and of 42.5 between the displacement constant and it. The relaxation of Equation 17 would shrink the admissible window from (0, 0.625] to (0, 0.096], a factor of 6.5 in the step size and, by Theorem 36 (i), the same factor in the residual constant α = b+1. The gap is a direct measurement of the non-normality of the composite Jacobian described in Proposition 25: a deep denoiser expands strongly in directions poorly aligned with the residual, which the quotient Equation 16 sees and the pair (σ, L_*T*) does not. Figure 7 shows the spread of the three spectral quantities over the certification sample; σ separates the certified models from the unconstrained one cleanly, while L_*T* does not.

8.4 Crossing from the nonexpansive class into the enriched class

Table 5 compares operators that differ in their training. To isolate the effect of the operator class with everything else held fixed, we take the residual of a single certified denoiser and rescale it, forming D_*c* = I−cR_θ for a range of gains c. Raising c has two effects that pull in opposite directions on the two relaxations: it raises ∥JT∥₂, and it lowers λmax(Js_*T*), because the symmetric part of the residual Jacobian enters with the factor c. Table 6 records the family.

Table 6. A family of denoisers D_*c* = I−cR_θ built from one certified residual.

The family leaves the nonexpansive class and enters the enriched one with every other ingredient held fixed. Denoising and restoration are peak signal-to-noise ratios in dB.

The family crosses the boundary between the two classes at a gain between 1.3 and 1.6. Below the crossing the composite is nonexpansive and b^⋆ = 0; above it the composite is expansive, with L_*T* reaching 2.71 at c = 2.5, and Theorem 22 still certifies it, with b^⋆ growing from 0.231 to 0.882 and the window contracting from (0, 0.812] to (0, 0.531]. The last three rows of the table are operators that the nonexpansive and firmly nonexpansive conditions of Table 1 exclude and that Theorem 17 covers. The ratio between b_fov and b^⋆ grows along the family, from 5.3 at c = 1.6 to 6.7 at c = 2.5, so the coarser certificate becomes less informative exactly where the operator becomes more expansive. Restoration falls slowly along the family, from 25.85 to 24.93 dB, which is expected since c = 1 is the gain the denoiser was trained for; what the family isolates is the class transition and the behavior of the certificate, not a gain in quality.

8.5 Restoration

Table 7 reports restoration on the three tasks, and Figure 8 shows one deblurring instance.

Table 7. Restoration averaged over 4 test images.

Certified operators are run with the greedy relaxation inside their certified window; the unconstrained operator has no window and is run unrelaxed. The iteration column counts steps to a relative fixed-point residual below 10⁻⁵.

The certificate is not obtained at the price of quality on these problems. The unconstrained denoiser, which carries no guarantee, has the lowest structural similarity on all three tasks and the lowest peak signal-to-noise ratio on two of them: 25.24 dB against 25.85 to 26.12 on Gaussian deblurring and 25.88 against 26.52 to 26.92 on super-resolution. On motion deblurring it sits between the two certified models in peak signal-to-noise ratio, 25.58 against 25.42 and 25.87, while its structural similarity, 0.7812, is below both. Constraining the composite acts here as a useful regularizer of the fixed point rather than as a tax on it, which is the opposite of what the received account of the trade-off predicts. The statement concerns small networks on modest crops rather than the state of the art in restoration.

The comparison the theory is built for is the next row. On Gaussian deblurring the enriched operator, with L_*T* = 2.08 and b^⋆ = 0.600, reaches 25.89 dB and 0.8156 structural similarity, against 25.85 and 0.8155 for the nonexpansive operator of identical architecture. The two are within four hundredths of a decibel, so what the numbers support is not that enrichment improves restoration but that it costs nothing: an operator that the nonexpansive and firmly nonexpansive conditions of Table 1 reject is certified here, and it restores as well as the operator those conditions accept. The larger class is therefore available at no measured cost, and which operator to place in it is a question for the training objective.

The last row concerns the training objective itself. The denoiser whose loss constrains only the symmetric part, with no Lipschitz term, is the best certified model on all three tasks, at 26.12, 25.87 and 26.92 dB. Proposition 15 (iv) and Theorem 22 single out precisely this quantity, and constraining it alone is enough to certify while leaving the network more freedom than a Lipschitz constraint does. On these runs that freedom was not used to become expansive, so the resulting composite is nonexpansive and b^⋆ = 0; the enriched class was needed for the operator of the second reading and not for this one.

8.6 Step-size rules

Table 8 compares the constant relaxation, the greedy rule (Equation 25) and the trained policy inside the certified window, and Figure 9 shows the residual and the reconstruction quality along one run.

Table 8. Step-size rules inside the certified window.

The learned policy is trained by unrolled differentiation through the iteration with the adjoint recursion of Proposition 34 and the features (Equation 29). All three rules satisfy Equation 24 and are covered by Theorem 32 (i).

The greedy rule reproduces on images what Table 4 shows on the surrogate. For the nonexpansive operator on super-resolution, where both rules reach the tolerance on every image, it takes 36 iterations against 72 for the constant relaxation, a reduction of 50%. On Gaussian deblurring, again for the nonexpansive operator, the constant relaxation reaches the tolerance on two of the four images and the greedy rule on all four, so the entry 39 of Table 7 is a mean over four images and is not directly comparable with the constant-step entry; restricted to the two images on which both rules stop, the counts are 36.5 for the greedy rule against 72 for the constant one. For the enriched operator the constant relaxation does not reach the tolerance within the iteration budget while the greedy rule does. The clamped greedy value sits near the upper end of [a, 1−a] throughout, with mean between 0.976 and 0.999, so the acceleration comes from using the largest certified step consistently. The trained policy matches the greedy rule, both reaching the tolerance at 67 iterations to the nearest step, and does not improve on it here; the horizon it was trained on is short and the instances are close to affine, which is the regime in which Theorem 32 (ii) says the greedy step is already the one-step optimum. What Theorem 32 (i) guarantees, and what the run confirms, is that no realization of the policy leaves the convergent regime.

8.7 The window and the sampled certificate

Table 9 sweeps the return coefficient across and beyond the certified window on Gaussian deblurring with the enriched operator.

Table 9. Relaxation sweep across and beyond the certified window on Gaussian deblurring.

The sweep is run on a single test image with a budget of one hundred iterations, so its reconstruction values are not comparable with the averages of Table 7. The return coefficient is set to the stated multiple of the certified ᾱ; the multiple 1 is the certified value. At the certified value and at one and a half times it the residual falls below the stopping tolerance of 10⁻⁵; at half the certified value the iteration is convergent but slower and has not reached the tolerance after 100 steps; at twice the certified value and beyond it does not converge.

The iteration converges at the certified ᾱ and at 1.5 times it, and diverges at twice it, so on this instance the certificate is conservative by a factor between one and a half and two. That is a much tighter margin than the factor of 6.6 separating the certified window of Theorem 22 (iii) from the stability boundary on the linear surrogate of Section 8.1, which is what one expects: the surrogate was built to have a large gap between the worst-case field of values and the spectrum, while a trained denoiser is closer to the worst case that the certificate anticipates.

Finally, the certificate was measured on a finite sample, and Proposition 26 (ii) says what has to be checked for the guarantee to apply where Fejér monotonicity uses it, namely along the segments joining the iterates to the limit. We ran that check on one deblurring instance: sampling seven equally spaced points on the segment from the limit to the initial iterate and measuring ∥JN∥₂ at each, the largest value is 0.9929, below one at every point. The criterion that Proposition 26 (ii) tests pointwise is therefore satisfied at every sampled point of that segment, with a margin of 0.0071. Turning this into a statement about the whole segment needs the modulus ϰ, which we do not estimate for these networks, so what the measurement establishes is that the certificate is not violated where the iteration actually travels, and the entries of Table 7 are to be read in that sense rather than as claims over the whole hypercube.

9 Conclusion

The displacement and transfer–return calculus for nonexpansive-type operators supplies a working language for convergent Plug-and-Play image restoration. The central message is a relaxation of hypothesis: a learned denoiser need not be nonexpansive, only enriched nonexpansive, that is nonexpansive after averaging, for the Plug-and-Play iteration to converge, and the averaging is supplied by the inertial iteration, whose correction term is the exact quantity that reduces the enriched iteration on the composite operator to a classical iteration on its nonexpansive generator. Around this reduction, Theorem 22 identifies enrichment with infinitesimal cocoercivity of the displacement I−T, and from that identification the enrichment constant is read off a Rayleigh quotient. The characterization is exact rather than sufficient, it is attained on non-scalar operators, and it costs two automatic-differentiation primitives to evaluate. It also decides what happens as the denoiser approaches the identity, a regime that Proposition 23 shows to be favorable. Theorem 32 places a safeguarded controller inside the certified window, Theorem 36 gives rates and stability under inexact evaluation, and Theorem 40 extends the construction to arbitrary splittings. Section 8 confirms the certificates on a linear surrogate whose constants are known in closed form and, on deblurring, motion deblurring, and single-image super-resolution, exhibits, on the deblurring task, a trained composite operator that lies in the enriched class and in no smaller one, shows the characterization to be more than an order of magnitude tighter there than the field-of-values relaxation, and finds that this operator restores within four hundredths of a decibel of a certified nonexpansive composite of identical architecture.

Three limitations bound the scope of these results. The certificate is first order and is evaluated on a finite sample; Proposition 26 converts sampled measurements into a bound over a segment or over a covered set, but the covering version is impractical in high dimension, so the guarantees reported in Section 8 hold on the region actually certified and not, without further work, on the whole hypercube. The acceleration is established per step and measured among safeguarded schemes; on the linear surrogate the unrelaxed Picard iteration, which carries no guarantee there because the composite is expansive, is still the fastest, and a global complexity bound for the amortized policy is open. Finally, the restoration study uses a small residual network on modest crops, so the reported reconstruction values are not a claim about the state of the art in restoration, and they do not show that enrichment improves reconstruction. What they establish is that the enlarged class is occupied by operators the existing guarantees reject and that using it costs nothing in quality; whether it can be turned to advantage is a question for the training objective, on which Section 8.5 reports a single observation.

The results suggest a program whose scope the present experiments do not settle. Wherever a learned operator is iterated, in equilibrium networks [19, 20], in learned optimizers, or in unrolled reconstruction [17], the tension between expressivity and guaranteed convergence recurs, and the transfer–return mechanism offers a uniform way to address it: certify the composite as enriched, transfer to the nonexpansive generator, iterate with the inertial correction, and transfer back. Whether the mechanism reaches those settings is a conjecture, not a result: it requires the corresponding fixed-point operator to admit a finite enrichment constant, which is exactly what Theorem 22 makes checkable and what only a study of those operators can decide. The certificate is the sole problem-specific ingredient, and it reduces to eigenvalue computations that automatic differentiation already supports, so the checking is within reach. The imaging application developed here instantiates a principle in which convergence is certified and repaired rather than imposed on the network by force.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fams.2026.1945678/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Lixin Shen, Syracuse University, United States

Reviewed by: Zhenghan Fang, Johns Hopkins University, United States

Mary John, Abu Dhabi Polytechnic, United Arab Emirates