arXiv CS open access

A geometry-based deep equilibrium model for image restoration under multiplicative Gamma noise

2026 · arxiv_cs

arXiv CS · Papers · 2026

Conceptio Open Archive

Analyze & explore related research in Conceptio

Connect this record directly into your research workflow or AI agent. Explore citations across open repositories, query full text, and export structured proof bundles.

Explore in Search →Connect via MCP / API ↗

A geometry-based deep equilibrium model for image restoration under multiplicative Gamma noise Shengkun Yang ∗

Luca Ratti †

Zhichang Guo ‡

August 6, 2026

arXiv:2608.04944v1 [math.NA] 5 Aug 2026

Abstract We propose a deep learning framework for image restoration from images degraded by both multiplicative Gamma noise and blur. Unlike conventional deep equilibrium (DEQ) models that rely on implicit neural regularization, the proposed method learns an explicit and interpretable regularizer parameterized by geometric priors associated with surface area and mean curvature. To minimize the resulting variational model, we develop a mirror descent algorithm tailored to the commonly used Gamma-noise fidelity terms. Leveraging the Kurdyka–Lojasiewicz property for functions defined in o-minimal structures, we establish the global convergence of the generated iterates to a critical point. Experimental results on both grayscale and color image restoration demonstrate that the proposed method consistently outperforms representative model-based approaches while achieving performance comparable to state-of-the-art DEQ models based on implicit regularization, despite requiring substantially fewer trainable parameters.

Keywords: Deep equilibrium models; Learned geometric regularization; Gamma noise; Nonconvex optimization

Introduction

Multiplicative Gamma noise commonly arises in image restoration problems associated with coherent imaging systems. Representative applications include Synthetic Aperture Radar (SAR), laser imaging, ultrasound imaging, and tomographic modalities, all of which involve coherent interference or exponential attenuation phenomena [50, 27, 9]. Mathematically, the corresponding inverse problem aims to recover the unknown image x from degraded observations y = (Ax)η, where A denotes a known ill-conditioned measurement operator and η is a multidimensional Gamma random variable whose components have mean 1 and variance L1 , L being the equivalent number of looks [22]. The signal-dependent nature of the noise intensity, coupled with the complexity of the blurring operator, renders this inverse problem intrinsically ill-posed and challenging for both model formulation and algorithm design [47]. ∗

School of Mathematics, Harbin Institute of Technology, Harbin, 150001, Heilongjiang, China. Email: 23B912006@stu.hit.edu.cn. † Department of Mathematics, University of Bologna, Bologna, 40126, Emilia-Romagna, Italy. Email: luca.ratti5@unibo.it. Corresponding author. ‡ School of Mathematics, Harbin Institute of Technology, Harbin, 150001, Heilongjiang, China. Email: mathgzc@hit.edu.cn.

A geometry-based deep equilibrium model

Early studies on multiplicative noise removal were predominantly variational, incorporating prior image information through carefully designed regularization and fidelity terms. One of the earliest variational models for multiplicative denoising is the RLO model [50], which combines total variation (TV) regularization with fidelity constraints on the image mean and variance. To better align the model with the statistical description of Gamma noise, Aubert and Aujol [6] proposed a new fidelity term derived from a maximum a posteriori (MAP) perspective. The resulting model, however, is nonconvex, making convergence analysis considerably more challenging. Subsequent research introduced a variety of improvements, including variable transformation methods [51], variable splitting techniques [34, 65], and formulations involving auxiliary fidelity terms [25]. While TV-based regularization achieves satisfactory denoising performance, it is prone to artifacts such as contrast attenuation, corner smearing, and staircase effects [67]. Higher-order geometric regularization has since been proposed to mitigate these drawbacks: Zhang [64] introduced an adaptive Euler’s elastica model within a level-set framework that suppresses staircase artifacts while preserving geometric structure, and Yang et al. [57] developed a curvature-based elastica functional that further improves the preservation of edges and fine details. Beyond variational models, nonlocal approaches have also been extensively investigated. Exploiting sparse representations among similar patches, Deledalle et al. [24] developed the MuLoG+BM3D framework for speckle reduction in synthetic aperture radar images. Moreover, Liu et al. [44] constructed a low-rank distance regularization based on nonlocal patch similarity and proposed a convergent optimization algorithm for the resulting model. For a more comprehensive overview of recent advances in multiplicative Gamma noise removal, the interested reader is referred to [28, 53]. Despite the aforementioned model-based methods have shown promising performance, they rely heavily on handcrafted prior assumptions to characterize image structures through regularization. Consequently, when the prescribed priors fail to accurately capture the underlying image properties, the restoration quality may deteriorate significantly. Moreover, these methods typically require iterative stage-wise optimization and careful manual parameter tuning, which substantially increases the computational and practical burden. Model-based deep learning (MBDL) methods have recently emerged as an effective alternative to purely model-driven approaches by integrating physical measurement models with data-driven learned priors. Among them, the plug-and-play (PnP) framework has attracted considerable attention due to its flexibility in incorporating learned Gaussian denoisers into a wide range of imaging tasks [62, 54, 41]. Another representative MBDL paradigm is deep unfolding (DU), which interprets each iteration of an optimization algorithm as a network layer and trains the resulting architecture in an end-to-end supervised manner [21, 61, 45]. Despite their success, DU methods are often limited by the memory consumption and numerical instability associated with backpropagation through deep iterative architectures. Consequently, practical implementations usually involve only a small number of unfolding iterations, which restricts flexibility during inference. To address this limitation, deep equilibrium (DEQ) models [8], which can be interpreted as infinitely deep networks, have recently been introduced for solving imaging problems [31]. Rather than explicitly stacking network layers, DEQ learns a fixed point whose equilibrium state approximates the ground-truth image while requiring substantially less memory during training. Under suitable smoothness assumptions on the learned operator together with standard optimization conditions, the convergence of the corresponding fixedpoint iterations can be rigorously established. These theoretical guarantees have enabled successful applications to Poisson denoising [23] and Gaussian image restoration [33]. Nevertheless, extending the DEQ framework to multiplicative Gamma noise is substantially more involved, since the corresponding data-fidelity terms exhibit markedly stronger non-Lipschitz behavior, rendering convergence analyses based on the standard Lipschitz smoothness assumption no longer applicable. In this work, we develop a deep equilibrium (DEQ) framework for image restoration under joint blur and multiplicative Gamma noise degradations. From a modeling perspective, we propose an

A geometry-based deep equilibrium model

explicit variational regularizer parameterized by learnable convolution filters and influence functions. The resulting regularizer naturally falls within the Field-of-Experts (FoE) framework [49] while explicitly incorporating geometric priors based on surface area and mean curvature, thereby combining the interpretability of variational models with the flexibility of data-driven learning. From an algorithmic perspective, we introduce a novel Bregman potential together with a mirror descent inference scheme tailored to Gamma noise. From a theoretical perspective, we establish global convergence of the inference iterations to a critical point of the underlying objective functional. Finally, extensive experiments on both grayscale and color image restoration demonstrate that the proposed method consistently outperforms representative variational and deep learning approaches. The remainder of this paper is organized as follows. Section 2 introduces the parameterized variational model and the proposed mirror descent inference algorithm. Section 3 presents the convergence analysis of the lower-level problem. Section 4 describes the forward fixed-point computation and backward training strategy within the DEQ framework. Experimental results are reported in Section 5. Finally, Section 6 draws some conclusions of the proposed studies. To keep the manuscript selfcontained, the basic definitions and auxiliary theoretical results are provided in Appendices A and B.

Deep equilibrium modeling with mirror descent

Deep equilibrium models learn fixed points of parametric functions represented by neural network, and are equivalent to infinitely deep recurrent neural networks whose equilibrium states provide accurate approximations of the ground truth. By drawing connections to classical optimization algorithms, such as gradient descent and proximal gradient methods, DEQ models can be naturally adapted to (linear) inverse problems with additive Gaussian noise [31, 32, 33]. However, in the presence of multiplicative Gamma noise, the commonly used data-fidelity terms are generally not L-smooth, which introduces additional difficulties in both algorithm design and convergence analysis. Motivated by [68, 23], we incorporate a mirror descent strategy into the DEQ framework to address image restoration problems involving blur operators and Gamma noise. Under suitable constraints on the optimization parameters, convergence guarantees of the proposed method can be established. To systematically formulate the proposed DEQ framework, we consider a dataset {(x∗i , yi )}N i=1 , where x∗i ∈ Rn denotes the ground-truth image and yi ∈ Rm represents the corresponding degraded observation. Let fθ (·; yi ) : Rn → Rn be an image-to-image neural network parameterized by θ ∈ Θ ⊂ Rp . The proposed DEQ model is formulated as the bilevel problem: N ∗ 1 X 0 θ̂ ∈ arg min L(θ) = ℓ x∞ i θ; (xi , yi ) , xi , N θ∈Θ i=1 ∞ 0 s.t. xi θ; (xi , yi ) ∈ Fix(fθ (·; yi )).

(1)

Here, x0i denotes the initial estimate for the i-th sample, ℓ is the upper-level loss function (e.g., the mean squared error), and Fix(fθ (·; yi )) denotes the set of fixed points of the mapping fθ (·; yi ). Consequently, the lower-level problem amounts to computing a fixed point of fθ (·; yi ). To ensure consistency with the underlying variational model, we construct fθ such that its fixed points coincide with the critical points of a regularized objective functional Fθ . Specifically, under a box-constrained variational formulation [20, 59], we consider the composite energy Fθ (x; yi ) = D(x; yi ) + λRθ (x) + ι[0,a]n (x), where D(x; yi ) denotes the data-fidelity term enforcing consistency between the reconstruction x and the observed data yi , Rθ is a learnable regularizer equipped with explicit geometric priors, λ > 0 is

A geometry-based deep equilibrium model

the regularization parameter that balances the influence of the data term and the regularizer, and ι is the indicator function associated with the admissible set [0, a]n , defined by ( 0, x ∈ [0, a]n , ι[0,a]n (x) = +∞, otherwise. For imaging problems corrupted by multiplicative Gamma noise, the data-fidelity term D(x; yi ) is typically chosen as one of the following two formulations. • (KL-type fidelity) A generalized Kullback–Leibler fidelity, originally introduced for Poisson noise and subsequently adapted to Gamma noise [52, 25]: D

(x; yi ) =

m X

((Ax)q − (yi )q log(Ax)q ) .

q=1

• (AA-type fidelity) From the maximum a posteriori (MAP) perspective, the Aubert–Aujol fidelity term [6] is given by: m X (yi )q AA + log(Ax)q . D (x; yi ) = (Ax)q q=1

In the following subsections, we introduce the parameterized regularization term Rθ and the fixedpoint map fθ associated with the minimization of Fθ employed in this work.

2.1

Parameterized explicit variational regularization

A variety of model-based neural networks built upon implicit architectures have recently been proposed for imaging inverse problems [37, 31]. Although these implicit frameworks demonstrate promising empirical performance, they are often treated as black-box models and lack rigorous mathematical interpretability. To improve the interpretability of the proposed framework, we explicitly incorporate geometric priors into the network parameterization. Such geometric structures have been shown to play a fundamental role in network design and image reconstruction; see, e.g., [58, 40]. Inspired by the image area and mean curvature regularizers in [57], we consider the following parameterized regularization model equipped with learnable influence functions and linear convolution kernels: Rθ (x) =

C X n X

Na h X

c=1 j=1

l=1

Nc h X

i h i A x(c) ψla kl ∗ x(c) j

i h (c) i ψlc′ K x(c) A x

l′ =1

(2)

Here, C denotes the number of image channels, Na and Nc denote the numbers of learnable areaand curvature-related influence functions, respectively, kl are learnable convolution kernels, ∗ is the convolution operator, and b is a learnable scalar. Let Dh+ (Dh− ) and Dv+ (Dv− ) denote the forward (backward) finite difference operators along the horizontal and vertical directions, respectively. Based on these operators, we define the following two geometric quantities associated with x(c) . (i) The area-related geometric quantity is defined by q i h 2 2 A x(c) = Dh+ x(c) j + Dv+ x(c) j + ϵ, j

where ϵ is a learnable positive parameter.

A geometry-based deep equilibrium model

(ii) The curvature-related geometric quantity is defined as D+ x(c) D+ x(c) K x(c) = Dh− h (c) + Dv− v (c) , A(x ) A(x ) where the division is understood element-wise. The functions ψla : Rn → Rn and ψlc′ : Rn → Rn denote learnable influence functions associated with the area and curvature terms, respectively. Since the parameterizations of the two influence functions are analogous, we only describe the construction of ψla . Following the strategy of [21], the influence function ψla acts separately on each component and is described using a radial basis function (RBF) expansion as follows M X |zj − µk | a a , [ψl (z)]j = ωlk ϕ γk k=1

a are learnable coefficients, µ are uniformly spaced centers, and γ are the corresponding where ωlk k k scaling parameters. The learnable coefficients associated with ψlc′ are denoted by ωlc′ k . The basis

function ϕ is chosen as a Gaussian one, ϕ

|zj − µk | γk

(zj − µk )2 = exp − 2γk2

! .

This choice guarantees the analyticity of the learned influence functions and the resulting regularizer Rθ , which is required for the convergence analysis in Proposition 5. Finally, each convolution kernel is parameterized using a discrete cosine transform (DCT) basis [21], Nr X wlr br , kl = ∥wl ∥ r=1

where br ∈ RNr are fixed DCT basis filters and wl is a vector of Nr learnable coefficients. Under this parameterization, the resulting kernels naturally satisfy the unit-norm property and characterize the linear component of the regularization model, while the influence functions capture the corresponding a , ω c , w , ϵ, b} endows the nonlinear behavior. Consequently, the learnable parameter set θ = {ωlk lr l′ k proposed regularization model with considerable flexibility and expressive power. Notably, this flexibility is not merely qualitative: by appropriately specializing θ, the regularizer Rθ recovers a range of classical variational models as special cases, indicating that it constitutes a genuine generalization rather than an ad hoc parameterization. To illustrate this property, we present several representative special cases in the single-channel setting with Na = Nc = 1. • Total Variation (TV). If [ψla (·)]j ≡ 1, ∀j, ϵ = 0, b = 0, then (2) reduces to the classical total variation (TV) regularizer [50]. • Minimal Surface regularization. If [ψla (·)]j ≡ 1, ∀j, ϵ = 1, b = 0, then it becomes the minimal surface regularizer [46].

A geometry-based deep equilibrium model

• MC-Lp . If [ψla (·)]j ≡ 0, ∀j, [ψlc′ (z)]j = |zj |p , ϵ = 1, where |·| denotes the absolute value, with fixed constants p and b, then (2) reduces to the MC-Lp regularizer proposed in [66]. • Euler’s elastica regularization. If

p

 k l = Gσ ,

[ψla (z)]j = 

zj  , [ψlc′ (z)]j = (zj )2 , ϵ = 0, max zj j

where Gσ denotes the Gaussian convolution kernel with standard deviation σ, and b is a fixed constant, then (2) reduces to the Euler’s elastica regularizer proposed in [64]. With these choices, the first regularization term in (2) provides an estimate of the gray level of the (degraded) image, which scales the intensity of the regularization according to the intensity of each pixel. This is a desirable effect when dealing with multiplicative noise. • Hybrid geometric regularization. If

p

 kl = Gσ ,

[ψla (z)]j = 

zj  , [ψlc′ (z)]j = (zj )2 , ϵ = 1, max zj j

then the resulting regularizer coincides with the hybrid geometric regularization proposed in [57].

2.2

Mirror descent algorithms for deep equilibrium models

This section derives the fixed-point mapping fθ associated with the minimization of the parametric energy Fθ . In convex optimization, it is common to rely on iterative schemes to progressively refine an approximation of a minimizer: as a consequence, such minimizers can be seen as fixed points of the map employed to produce the next iterate given the previous one. For the multiplicative Gamma noise restoration problems considered in this work, the AA-type fidelity DAA (x; yi ) is not globally convex, while DAA (x; yi ) and DKL (x; yi ) both fail to satisfy the standard L-smoothness condition. Consequently, classical optimization methods, such as proximal gradient descent, cannot be directly applied with rigorous convergence guarantees. Motivated by [23], we adopt a mirror descent framework to address these difficulties. By introducing a Bregman potential function h, mirror descent extends the standard proximal gradient method and relaxes the conventional smoothness and convexity requirements imposed on the objective functional. To present the mirror descent framework in a general setting, we consider the following nonconvex optimization problem: inf Ψ(x) := f (x) + g(x), x∈C

where C = int dom h and h ∈ G(C) is a Bregman potential satisfying the kernel generating distance property (see Definition 2 for the precise definition). To establish the well-posedness of the subsequent algorithm in the considered nonconvex case, the following assumptions are needed (see [18]): Assumption 1. (i) h ∈ G(C) and C = dom h.

A geometry-based deep equilibrium model

(ii) f : Rn → (−∞, +∞] is proper, lower semicontinuous and convex, with dom f ∩ C ̸= ∅. (iii) g : Rn → (−∞, +∞] is proper and lower semicontinuous such that dom h ⊆ dom g. Moreover, g is continuously differentiable on C. (iv) The objective function Ψ is bounded below on C, i.e., inf Ψ(x) > −∞. x∈C

(v) For any θ > 0, the function h(x) + θf (x) is supercoercive, i.e., h(x) + θf (x) = +∞. ∥x∥ ∥x∥→+∞ lim

Starting from an initialization x0 ∈ int (dom h), the mirror descent iteration is recursively updated by n o 1 xk+1 = arg min Ψh (x) := f (x) + ⟨x − xk , ∇g(xk )⟩ + Dh (x, xk ) . τ x∈Rn

(3)

Here, τ > 0 denotes the step size and Dh is the Bregman distance induced by h; see Definition 4 for its definition. The update can be interpreted as minimizing a local model obtained by linearizing the differentiable component of the objective fuction g at the current iterate and regularizing the resulting subproblem through the non-Euclidean metric Dh . Remark 1. When h = ∥ · ∥2 , the above scheme reduces to the classical proximal gradient method. Hence, mirror descent may be viewed as a non-Euclidean generalization of proximal gradient descent [10]. Remark 2. Although the objective function Ψ(x) = f (x) + g(x) is generally nonconvex due to the nonconvexity of g, the local surrogate function Ψh arising from the linearization of g becomes convex under Assumption 1. Consequently, the resulting update mapping from int(dom h) to itself is well defined and single valued. Returning to the restoration problem considered in this work, we set f as the indicator function ι[0,a]n (x), while g consists of the parameterized regularization network λRθ (x) and the data-fidelity term D(x; y), with domain C = (0, +∞). Consequently, the generic objective Ψ introduced above reduces to the parameterized functional Fθ employed in the proposed DEQ model. Under Assumption 1, iteration (3) admits the following closed-form expression as the lower-level subproblem of the DEQ model: xk+1 = Proxhτι[0,a]n ∇h∗ ∇h(xk ) − τ λ∇Rθ (xk ) + ∇D(xk ; y) (4) where h∗ denotes the Fenchel conjugate of h [48]. The choice of the Bregman potential h depends on the specific form of the fidelity term. For the KL-type fidelity, a natural choice is Burg’s entropy [10], hKL (x) = −

n X

log(xj ).

j=1

which inherently enforces the positivity constraint on the solution. In contrast, for the AA-type fidelity term, this choice becomes unsuitable since it fails to satisfy the smoothness conditions required in the subsequent convergence analysis; see Lemma 2. In this case, we instead choose the potential function hAA (x) =

n X 1 j=1

A geometry-based deep equilibrium model

Although this potential hAA (x) does not explicitly preserve positivity, the constraint is still effectively enforced through the indicator function ι associated with the admissible set [0, a]n , making it an appropriate choice within the proposed mirror descent framework. As a particular instance of the Bregman proximal gradient framework, the Bregman proximal step in mirror descent is generally not equivalent to an Euclidean projection. Nevertheless, under the specific choices of the potential functions hKL (x) and hAA (x), corresponding respectively to the KL-type and AA-type fidelity terms, this equivalence still holds. Pn −1 n Proposition 1. Let h(x) = j=1 xj . Then, for any x ∈ R++ , the Bregman proximal operator associated with the indicator function ι[0,a]n coincides with the Euclidean projection onto [0, a]n . Proof. The Bregman proximal operator associated with ι[0,a]n is defined by Proxhτι[0,a]n (x) = arg minn τ ι[0,a]n (u) + Dh (u, x) u∈R " !# n X uj − xj 1 1 = arg minn τ ι[0,a] (uj ) + − + . u∈R uj xj x2j

(5)

j=1

Since the objective function is separable, the minimization decouples componentwise. Define D̃h (ũ, x̃) =

1 ũ − x̃ 1 − + , ũ x̃ x̃2

x̃ > 0.

Then (5) reduces to Proxhτι[0,a]n (x) = Proxhτι[0,a] (x1 ), . . . , Proxhτι[0,a] (xn ) , where Proxhτι[0,a] (xj ) = arg min D̃h (uj , xj ). 0≤uj ≤a

Differentiating D̃h with respect to ũ yields ∂ũ D̃h (ũ, x̃) =

ũ2 − x̃2 . ũ2 x̃2

Hence, D̃h (·, x̃) is strictly decreasing on (0, x̃) and strictly increasing on (x̃, ∞), so its unique minimizer is attained at ũ = x̃. Therefore, Proxhτι[0,a] (xj ) = min{xj , a}, which is precisely the Euclidean projection onto [0, a]. Consequently, Proxhτι[0,a]n (x) = Π[0,a]n (x), where Π[0,a]n denotes the Euclidean projection onto [0, a]n .

□

We conclude this section by presenting the explicit mirror descent updates for image restoration under multiplicative Gamma noise. The resulting closed-form iterations are derived for the two fidelity terms most commonly employed in applications. For the KL-type fidelity term, the corresponding update is given by xk k+1 x = Π[0,a]n . (6) 1 + τ xk (λ∇Rθ (xk ) + ∇DKL (xk ; y))

A geometry-based deep equilibrium model

For the AA-type fidelity term, the mirror descent iteration takes the form x

k+1

= Π[0,a]n

xk p 1 + τ (xk )2 (λ∇Rθ (xk ) + ∇DAA (xk ; y))

! .

(7)

These closed-form expressions constitute the computational core of the proposed method and will serve as the forward operator in the subsequent DEQ framework. Remark 3. The update formulas (6) and (7) are derived under the assumption that the argument of the projection operator Π[0,a]n belongs to Rn++ . To ensure that this assumption is satisfied throughout the iterative process, we introduce in Section 4 a backtracking strategy that preserves the positivity of the projection argument.

Convergence analysis of mirror descent with smooth-adaptable functionals

In this section, we establish the convergence of the mirror descent iterations to a critical point of the objective functional Ψ, where Ψ = Fθ . A common assumption in the convergence analysis of firstorder optimization methods is the Lipschitz continuity of the gradient; see, e.g., [5, 17]. However, this assumption is violated by both the AA-type and KL-type data-fidelity terms arising in multiplicative Gamma noise restoration. Moreover, the objective functional Ψ is nonconvex, so classical convex optimization theory is not applicable. To overcome these difficulties, we follow the analytical framework of [18, 23] and establish the convergence of the mirror descent scheme (4) under weaker assumptions satisfied by the proposed model. Specifically, the classical Lipschitz smoothness requirement is replaced by the notion of L-smooth adaptability with respect to the Bregman potential h, while the convexity assumption is replaced by the Kurdyka–Lojasiewicz property. Since the convergence analysis for the KL-type fidelity follows directly from the arguments in [23], we restrict our attention to the AA-type fidelity, which requires additional technical developments. We begin by recalling the notion of L-smooth adaptability (L-smad). Definition 1. Let h ∈ G(C) with C = int(dom h), and let g : Rn → (−∞, +∞] be a proper lower semicontinuous function such that dom h ⊂ dom g. Assume that g is continuously differentiable on C. The pair (g, h) is said to satisfy the L-smooth adaptability (L-smad) condition on C if: there exists L > 0 such that Lh − g is convex on C. Remark 4. Strictly speaking, the complete L-smad condition additionally requires Lh+g to be convex. However, for the convergence analysis considered here, it is sufficient to assume only the convexity of Lh − g [18]. Let aq ∈ Rn+ denote the q-th row of the forward operator A ∈ Rm×n . Assume that aq ̸= 0 for all q = 1, . . . , m. Then, for every x ∈ Rn++ , ⟨aq , x⟩ > 0,

q = 1, . . . , m.

This is a standard assumption for blur operators in imaging inverse problems; see, e.g., [25, 12]. Under this condition, one can derive a mild sufficient criterion ensuring the L-smad property for the AA-type fidelity term. P Lemma 2. Let g(x) = DAA (x; y) and hAA (x) = nj=1 x1j . Then, for any L≥M

m X q=1

yq ,

A geometry-based deep equilibrium model

where

max1≤j≤n {a2qj }

M = max

, 1≤q≤m min1≤j≤n {a3 qj , aqj > 0}

the function LhAA − g is convex on Rn++ . Proof. Since g and hAA are C 2 on Rn++ , it suffices to find a constant L > 0 such that, for any x ∈ Rn++ and any direction d ∈ Rn , L⟨∇2 hAA (x)d, d⟩ ≥ ⟨∇2 g(x)d, d⟩. Firstly, it is straightforward to verify that, for any d ∈ Rn , L⟨∇2 hAA (x)d, d⟩ = 2L

n X d2j j=1

x3j

On the other hand, by the definition of g, we have ∇g(x) =

m X

yq 1 aq − aq , ⟨aq , x⟩ ⟨aq , x⟩2

q=1

and ⟨∇2 g(x)d, d⟩ =

m X q=1

−

⟨aq , d⟩2 2yq ⟨aq , d⟩2 X ⟨aq , d⟩2 + ≤ 2y . q ⟨aq , x⟩2 ⟨aq , x⟩3 ⟨aq , x⟩3 q=1

By the Cauchy-Schwarz inequality and the assumption that aq ∈ Rn+ with aq ̸= 0 for all q = 1, . . . , m, we obtain Pn n 2 3 X d2j ⟨aq , d⟩2 j=1 aqj xj ≤ . ⟨aq , x⟩3 ⟨aq , x⟩3 x3j j=1

Let Mq =

maxj {a2qj } minj {a3qj :aqj >0}

. Then it follows that Pn

P 2 3 maxj {a2qj } nj=1 x3j j=1 aqj xj P ≤ ≤ Mq . ⟨aq , x⟩3 minj {a3qj : aqj > 0} nj=1 (xj )

Consequently, we have 2

m X

⟨aq , d⟩2 ⟨aq , x⟩3 q=1   m n X X d2j  ≤ 2yq Mq  3 x q=1 j=1 j    m n 2 X X dj , ≤M  2yq   3 x q=1 j=1 j

⟨∇ g(x)d, d⟩ ≤

2yq

where M = maxq∈{1,...,m} Mq > 0 is a constant that depends only on the matrix A and is independent of x and d. Therefore it suffices to choose m X L≥M yq , q=1

which guarantees the convexity of Lh − g.

□

A geometry-based deep equilibrium model

Remark 5. Burg’s entropy is the standard choice of Bregman potential for the KL-type fidelity and plays a central role in the existing convergence theory. It is therefore natural to consider using the same potential for the AA-type fidelity. Nevertheless, the corresponding L-smad condition generally fails because the AA-type fidelity possesses a stronger singularity near the origin, induced by its 1/xtype structure. The proposed potential can thus be viewed as an extension of Burg’s entropy tailored to the stronger singular behavior of the AA-type fidelity, while retaining the analytical properties required for our convergence analysis. For the KL-type fidelity term, the L-smad condition associated with the Bregman potential hKL admits a considerably simpler characterization than in the AA case. In fact, the corresponding result is classical; see [10, Lemma 7], which provides an explicit condition on L ensuring the convexity of Lh − g. Lemma 3. Let g(x) = DKL (x; y), and let hKL (x) denote the Burg’s entropy. If the constant L satisfies m X L≥ yq , q=1

then the function Lh − g is convex on Rn++ . Based on Lemmas 2 and 3, establishing the L-smad property for the objective function D(x; y)+λRθ reduces to verifying it for the pair (λRθ , h), since the sum of convex functions remains convex. As observed in [35], verifying the L-smad condition over the entire domain int(dom h) is generally difficult. A practical simplification is to restrict the analysis to a compact subset of conv(dom Rθ ) ∩ int(dom h). Motivated by this observation, we introduce the indicator function ι[0,a]n to constrain the admissible set to [0, a]n . Such a restriction is natural in imaging applications, where image intensities are typically normalized to ranges such as [0, 1]n or [0, 255]n . Under this bounded-domain setting, the selection rule for the constant L′ ensuring the L-smad property of (λRθ , h) can be provided. LemmaP4 (L-smad property of (λRθ , h)). Let Rθ : (0, a]n → R be a real analytic function and let h(x) = nj=1 x1j . Then there exists a constant L′ > 0 such that the pair (λRθ , h) satisfies the L-smad property on (0, a]n . In particular, letting ρ(·) denote the spectral radius, M := supx∈(0,a]n ρ ∇2 Rθ (x) < 3

∞, then any constant L′ ≥ a2 λM is sufficient. Proof. Since both Rθ and h are both C 2 -functions on their domains, the L-smad property of the pair (λRθ , h) can be characterized through their Hessian matrices. In particular, it suffices to show that there exists L′ > 0 such that [10] L′ ∇2 h(x) − λ∇2 Rθ (x) ⪰ 0,

∀x ∈ (0, a]n .

where ⪰ 0 denotes positive semidefiniteness. Equivalently, one seeks L′ > 0 satisfying inf λmin L′ ∇2 h(x) − λ∇2 Rθ (x) ≥ 0, x∈(0,a]n

with λmin (·) the minimum eigenvalue. Therefore, following the argument of [23, Lemma B.4], it remains to establish a positive lower bound for the eigenvalues of ∇2 h(x) and a uniform upper bound for those of ∇2 Rθ (x) on (0, a]n . Since Rθ is analytic, it is Lipschitz smooth on the bounded set (0, a]n ⊂ [0, a]n . Consequently, there exists a constant M > 0 such that sup ρ ∇2 Rθ (x) ≤ M, x∈(0,a]n

A geometry-based deep equilibrium model

where ρ(·)P denotes the spectral radius. On the other hand, consider the potential function h(x) = AA h (x) = nj=1 x1j . For x ∈ (0, a]n , we have ∇2 h(x) = 2diag

1 x3

with

1 = x3

1 1 ,..., 3 xn x31

Hence, the Hessian matrix of h satisfies inf

x∈(0,a]n

2 2 ≥ 3. 3 a x∈(0,a] x

λmin ∇2 h(x) ≥ inf

Therefore, whenever L′ ≥ a2 λM > 0, we obtain inf

x∈(0,a]n

λmin L′ ∇2 h(x) − λ∇2 Rθ (x) ≥ 0,

This proves that the pair (λRθ , h) satisfies the L-smad condition on (0, a]n .

□

Finally, combining Lemmas 2 and 4, we obtain an estimate of the L-smad constant for the objective functional F = DAA + λRθ . In particular, an admissible choice is L̄ = L + L′ , yielding L̄hAA − F

is convex on int(dom hAA ).

(8)

Hence, the objective functional satisfies the L-smooth adaptability condition. To establish the convergence of the mirror descent algorithm in this nonconvex setting, in addition to the L-smad property, we also rely on the Kurdyka-Lojasiewicz (KL) property of the minimized functional. Directly verifying the KL property from its definition is generally challenging. Fortunately, a broad class of functions arising in optimization and imaging are known to satisfy the KL property through o-minimal geometry. A detailed review of KL functions and o-minimal structures is reported in Appendix B. We therefore establish the KL property of the objective functional Ψ by showing that it is definable in an o-minimal structure, which in turn guarantees the desired KL property and provides the remaining key ingredient for proving convergence to a critical point. Proposition 5. Let Ψ = DAA (x; y) + ι[0,a]n (x) + λRθ (x). Then Ψ satisfies the KL property at any x ∈ crit(Ψ). Moreover, Ψ is coercive, and the set S∞ = {x ∈ dom Ψ : Ψ(x) = +∞} is closed. Proof. Following [23, Proposition 4.3], it suffices to verify the properness, lower semicontinuity, continuity on the domain, and definability of Ψ to establish the KL property. 1. Ψ is proper, lower semicontinuous, and continuous on its domain. Indeed, for any x ∈ (0, a)n , one has Ψ(x) < +∞, implying that Ψ is proper. Since [0, a]n is closed, the indicator function ι[0,a]n is lower semicontinuous [48]. Consequently, Ψ = DAA (x; b) + ι[0,a]n (x) + λRθ (x) being the sum of lower semicontinuous functions, is itself lower semicontinuous. Moreover, for x ∈ dom(Ψ), one has Ψ = DAA (x; y) + λRθ (x), and continuity follows directly from the continuity of DAA and Rθ . 2. Ψ is a definable fuction in Ran,exp . We first recall the AA-type fidelity term DAA (x; y) =

m X q=1

X yq + log(Ax)q = DqAA (x). (Ax)q q=1

A geometry-based deep equilibrium model

To establish its definability, it suffices to consider the elementary function f (x) = 1/x on the domain {x > 0}. Its graph is given by Γ(f ) = {(x, t) ∈ R2 : tx − 1 = 0, x > 0}, which is a semialgebraic set. Hence, f is definable in the semialgebraic structure (R, +, ·, (r)r∈R ) and, consequently, definable in the larger o-minimal structure Ran,exp . Since the logarithm is also definable in Ran,exp , each component in DqAA (x) , being the sum of a logarithmic term and a reciprocal term, is definable. By closure of definable functions under finite sums, the fidelity term DAA is therefore a definable function in Ran,exp . We next show that ι[0,a]n is semialgebraic, and hence definable in Ran,exp . Consider the polynomials Pj0 (x) = 0, Pj1 (x) = −1, Pj2 (x) = xj (xj − a), Pj3 (x) = xj , Pj4 (x) = xj − a, ∀j = 1 . . . , n. Then the cube [0, a]n can be expressed as a finite union of finite intersections of semialgebraic sets: \ [ [0, a]n = {x ∈ Rn : Pj3 (x) = 0, Pj1 (x) < 0} S0 ,Sa ⊆{1,...,n} S0 ∩Sa =∅

j∈S0

∩

{x ∈ Rn : Pj4 (x) = 0, Pj1 (x) < 0}

j∈Sa

∩

{x ∈ Rn : Pj0 (x) = 0, Pj2 (x) < 0} .

j ∈S / 0 ∪Sa

Hence [0, a]n is a semialgebraic set, and its indicator function ι[0,a]n is a semialgebraic function. The network Rθ consists of compositions of affine mappings, convolution operators, and nonlinear activation functions constructed from Gaussian basis expansions. Therefore, Rθ is definable in Ran,exp . Since all three components of Ψ are definable in the o-minimal structure Ran,exp , the objective functional Ψ is also definable in Ran,exp . Consequently, Theorem 9 implies that Ψ satisfies the KurdykaLojasiewicz property. 3. Ψ is coercive. Owing to the presence of the indicator function ι[0,a]n , one has Ψ(x) = +∞ whenever x ∈ / [0, a]n . Hence ∥x∥ → +∞ ⇒ Ψ(x) → +∞, which establishes coercivity. 4. S∞ is closed. In the present setting, dom(Ψ) = [0, a]n . Moreover, for x ∈ dom(Ψ), the condition x ∈ S∞ is equivalent to the existence of some q ∈ {1, . . . , m} such that (Ax)q = 0. Let aq denote the q-th row of A, and let Ker(aq ) denote its null space. Then n S∞ = ∪m q=1 (Ker(aq ) ∩ [0, a ]) .

Since the mapping fq (x) = ⟨aq , x⟩ is continuous and {0} ⊂ R is closed, each set Ker(aq ) = fq−1 ({0}) is closed in Rn . Therefore, S∞ , being a finite union and intersection of closed sets, is itself closed. □

A geometry-based deep equilibrium model

Remark 6. An alternative proof of the second component can be obtained by extending the argument of [23, Proposition 4.3]. For this reason, in the previous result, we discussed the properties of the set S∞ , which are necessary to prove [23, Proposition 4.3] but not Theorem 6. It should be noted, however, that the objective functional Ψ considered here is not globally subanalytic, but only locally subanalytic. Nevertheless, by [23, Corollary 4.5], this local characterization is still sufficient to establish that Ψ is a regularized KL functional. In the present work, we instead adopt the o-minimal framework, since it covers a substantially broader class of functions and avoids the technical construction of auxiliary real-analytic functions required in [23, SM 2.4]. This provides a more direct and unified proof of the KL property. Theorem 6 (Convergence of mirror descent with AA-type fidelity). Let Ψ = DAA (x; y) + ι[0,a]n (x) + λRθ (x) be the objective functional. Assume that the step size τ satisfies 0 < τ L̄ < 1, where L̄ is defined in (8). Then the sequence {xk } produced by the iteration (7) converges to a critical point of Ψ. Proof. The result follows by verifying the assumptions of [35, Theorem 3] (see also [18, Theorem 4.1]). Let Ch = dom(hAA ). The required conditions are: (i) hAA : Ch ⊂ Rn → R ∪ {+∞} is a proper, convex, C 2 Legendre function, and is strongly convex on every bounded convex subset of Ch . (ii) F (x) := DAA (x; y) + λRθ (x) is proper, continuously differentiable, and bounded from below on conv(dom(ι[0,a]n )) ∩ int(dom(hAA )). (iii) The functional L̄hAA − F is convex on conv(dom(ι[0,a]n )) ∩ int(dom(hAA )). (iv) Ψ = F + ι[0,a]n is bounded from below, coercive, and satisfies the KL property. (v) For every α > 0, both ∇hAA and ∇F are Lipschitz continuous on the sublevel set {x ∈ conv(dom(ι[0,a]n )) ∩ int(dom(hAA )) : Ψ(x) ≤ α}. (vi) The set [0, a]n is nonempty, closed, convex, and semialgebraic. (vii) For every xk ∈ int(dom(hAA )), the solution set n o 1 Mτ (xk ) := arg minn ι[0,a]n (x) + ⟨x − xk , ∇F (xk )⟩ + DhAA (x, xk ) x∈R τ is nonempty and contained in int(dom(hAA )). Proposition 5 establishes condition (vi). In conjunction with the lower bound DAA (x; y) ≥ m +

m X

log yq ,

q=1

it also yields conditions (ii) and (iv). Condition (iii) follows from Lemmas 2 and P 4. Furthermore, by Definition 3 and the explicit form of the chosen Bregman potential, hAA (x) = nj=1 x1j is a Legendre function and satisfies condition (i). It therefore remains to verify (v) and (vii). For (v), let Ωα := {x ∈ conv(dom(ι[0,a]n )) ∩ int(dom(hAA )) : Ψ(x) ≤ α}. Since dom ι[0,a]n = [0, a]n and dom(hAA ) = (0, +∞)n , we have Ωα ⊂ (0, a]n .

A geometry-based deep equilibrium model

Hence Ωα is bounded. Because Ψ(x) ≤ α on Ωα , there exists some β > 0, such that Ωα ⊂ [β, a]n . Since both hAA and F are C 2 on [β, a]n , their Hessians are bounded on Ωα . Consequently, ∇hAA and ∇F are Lipschitz continuous on Ωα , establishing (v). To verify (vii), note that ι[0,a]n is convex and hAA is strictly convex on int(dom(hAA )). Therefore, the surrogate objective 1 x 7→ ι[0,a]n (x) + ⟨x − xk , ∇F (xk )⟩ + DhAA (x, xk ) τ is strictly convex on [0, a]n . Since ι[0,a]n restricts the feasible set to the compact set [0, a]n , the above minimization problem admits a unique global minimizer. Furthermore, for any xk ∈ int(dom(hAA )), the Bregman distance satisfies DhAA (x, xk ) = +∞ whenever x ∈ ∂ dom(hAA ). Hence the minimizer cannot lie on the boundary of dom(hAA ), implying that Mτ (xk ) ⊂ int(dom(hAA )). Therefore Mτ is a well-defined single-valued mapping from int(dom(hAA )) into itself, consistent with [18, Remark 3.1]. Similarly, one can see from the iterative update (7) that when x0 ∈ int(dom(hAA )), the pre-projection result still lies in Rn++ rather than in ∂ dom(hAA ). All assumptions of [35, Theorem 3] are thus satisfied, and the conclusion follows immediately. □ Proposition 7 (Existence and convergence of fixed points). Define fθ (x; y) = Π[0,a]n

x p 1 + τ x2 (∇Rθ (x) + ∇DAA (x; y))

(9)

Then the AA-type update rule (7) can be written as xk+1 = fθ (xk ; y). Suppose that 0 < τ L̄ < 1, where L̄ is defined in (8), and let x0 ∈ int(dom(hAA )). Then the fixed-point set Fix(fθ (·; y)) = {x ∈ Rn : x = fθ (x; y)} is nonempty. Moreover, the sequence {xk } generated by (7) converges to a point in Fix(fθ (·; y)). Proof. By Theorem 6, the sequence {xk } converges to some limit point x∞ ∈ int(dom(hAA )). Since Rθ and DAA are analytic, the mapping fθ defined in (9) is continuous. Passing to the limit in xk+1 = fθ (xk ; y) gives x∞ = lim xk+1 = lim fθ (xk ; y) = fθ ( lim xk ; y) = fθ (x∞ ; y). k→∞

k→∞

Hence x∞ ∈ Fix(fθ (·; y)) showing that the fixed-point set is nonempty. the proof is complete.

Since {xk } converges to x∞ ,

□

Forward and gradient computation in DEQs

Having established the theoretical framework of the proposed DEQ model at inference stage, two computational challenges remain to be addressed in practice. First, during the forward pass, one must efficiently compute a fixed point of the mapping fθ for a given observation and network parameters θ ∈ Θ. Second, since a DEQ model can be interpreted as an infinitely deep network, standard backpropagation is generally infeasible due to its prohibitive memory requirements. Consequently, after solving the lower-level fixed-point problem in (1) and evaluating the upper-level loss, an efficient and memory-intensive strategy is required to optimize the network parameters θ. In the remainder of this section, we describe the numerical techniques employed to address these two challenges.

A geometry-based deep equilibrium model

4.1

Fixed-point computation via backtracking mirror descent

Within the DEQ framework, both training and inference require the computation of a fixed point of the mapping fθ . A straightforward approach is to apply the fixed-point iteration induced by (6) and (7). The efficiency and convergence behavior of this iteration, however, depend critically on the choice of the step size τ . The convergence of the AA iteration (7) is guaranteed by Theorem 6, provided that τ < L̄1 , where L̄ = L + L′ with the constants L and L′ given in Lemmas 2 and 4, respectively. For the KL iteration (6), the same step-size condition is required in [35, Theorem 3], where L is characterized in Lemma 3 and L′ is derived in [23, SM 2.3.3]. Nevertheless, such a strategy is impractical for several reasons. First, as discussed in Lemma 4, estimating the L-smad constant L′ of the network requires evaluating bounds on the spectral radius of the Hessian of Rθ , which can be computationally expensive. Second, the resulting estimates are often overly conservative [35], leading to unnecessarily small step sizes and consequently slow convergence. Finally, the validity of the update itself imposes additional restrictions on τ . Indeed, inspection of the update formula (7) shows that the denominator terms must remain strictly positive; otherwise, the square-root operation becomes undefined. Therefore, an appropriate step-size selection mechanism is essential not only for convergence but also for the welldefinedness of the iteration. Remark 7. (Step size for the KL fidelity). Although the update formula (6) associated with the KL-type fidelity term does not explicitly involve a square-root operation, positivity constraints remain necessary. Indeed, for Burg’s entropy h(x) = − log x, the Fenchel conjugate is given by h∗ (y) = − log(−y) − 1, whose domain is (−∞, 0) [10]. Comparing (4) with the update formula (6) shows that the corresponding denominator terms must still remain positive in order for the iteration to be well defined. Motivated by the backtracking principle introduced in [11], we develop an adaptive step-size strategy specifically for the proposed mirror descent scheme. Starting from an initial step size τ0 > 0, the step size is decreased adaptively until the sufficient decrease condition (10) and the positivity condition (11) are both satisfied. The former is formulated in terms of the Bregman distance Dh , while the latter guarantees the well-definedness of the subsequent iterate. Consequently, the explicit estimation of the L-smad constants is avoided without sacrificing the convergence guarantees of the mirror descent scheme. The complete procedure is summarized in Algorithm 1. Having described the algorithm, we now present the following proposition regarding its convergence and the finite termination of the backtracking procedure. Theorem 8. The backtracking procedure in Algorithm 1 terminates after finitely many iterations. Moreover, the sequence {xk } generated by (4) with step sizes determined by the backtracking strategy converges to a critical point of Ψ. Proof. The backtracking strategy preserves the sufficient decrease property required in the proof of Theorem 6. Therefore, the convergence analysis extends to the backtracking scheme, and the convergence of the generated sequence follows by the same argument; see [35, Proposition 4]. It therefore remains to show that the backtracking procedure terminates in finitely many steps. By [35, Lemma 2], for any step size τ > 0, the following sufficient decrease estimate holds: 1 Ψ(xk ) − Ψ(Mτ (xk )) ≥ − L̄ Dh (Mτ (xk ), xk ), τ where L̄ > 0 denotes the Lipschitz constant of λRθ + D. Consequently, whenever τ < 1−γ , we have L γ 1 τ − L̄ > τ , which implies γ Ψ(xk ) − Ψ(Mτ (xk )) > Dh (Mτ (xk ), xk ). τ

A geometry-based deep equilibrium model

Algorithm 1 Adaptive Bregman Backtracking Step-Size Selection Input: Current iterate xk ∈ int(dom h), initial step size τ0 > 0, sufficient decrease parameter γ ∈ (0, 1), backtracking factor η ∈ (0, 1). Output: Updated iterate xk+1 . 1: Set τ ← τ0 . 2: Compute Mτ (xk ) = Proxhτι[0,a]n ∇h∗ ∇h(xk ) − τ ∇Rθ (xk ) + ∇D(xk ; y) . 3: while

γ Dh (Mτ (xk ), xk ) τ

(10)

∇h(xk ) − τ ∇Rθ (xk ) + ∇D(xk ; y) ≥ 0

(11)

Ψ(xk ) − Ψ(Mτ (xk )) ≤ or do 4: τ ← ητ . 5: Recompute Mτ (xk ). 6: end while 7: Set

xk+1 ← Mτ (xk ). 8: return xk+1 .

Hence the first acceptance condition in the while-loop of Algorithm 1 is automatically satisfied for τ chosen in this way. Next, define Th (τ ) = ∇h(xk ) − τ ∇Rθ (xk ) + ∇D(xk ; y) . At τ = 0, one has ( − x1k , Th (0) = − (x1k )2 ,

for the KL-type fidelity, for the AA-type fidelity,

where 1 ∈ Rn denotes the vector of ones. When xk ∈ (0, a]n , Th (0) is strictly negative componentwise in both cases. By continuity of Th , there exists δD > 0 such that Th (τ ) < 0 for all 0 < τ < δD . Combining the two observations, all acceptance conditions are fulfilled for τ < min{δD , 1−γ }. Since L̄ the backtracking strategy decreases the step size geometrically, it must terminate after finitely many reductions. This completes the proof. □ Finally, for the implementation of the forward pass, it is important to specify how to evaluate the gradient ∇Rθ appearing in (4). To do so, two strategies are available. The most direct approach is to employ automatic differentiation to evaluate ∇Rθ and substitute it into the fixed-point mapping fθ . Alternatively, following the idea of [21], one may exploit the explicit form of the regularizer in (2) to derive ∇Rθ analytically, which can then be incorporated into the forward iteration (4). Although an analytical implementation may offer additional computational efficiency, it is generally less flexible and requires problem-dependent derivations. In particular, explicitly differentiating the terms involved in Rθ , such as the surface area A and curvature K, often results in cumbersome expressions. Moreover, certain components of the functional may vary depending on the specific application, including the normalization strategies adopted in ψla and the treatment of convolutions across multiple image channels. By contrast, automatic differentiation requires no model-specific derivation, is straightforward to implement, and benefits from the optimized differentiation routines provided by modern deep learning frameworks such as PyTorch. For these reasons, we adopt the automatic differentiation implementation throughout this work.

A geometry-based deep equilibrium model

4.2

Gradient computation via implicit differentiation

After computing the lower-level solution of the bilevel problem (1) via mirror descent, we evaluate the upper-level objective function L(θ) =

N ∗ 1 X 0 ℓ x∞ i θ; (xi , yi ) , xi . N i=1

The network parameters θ are then updated using standard optimization algorithms, such as first-order stochastic methods (e.g., Adam [38]) or quasi-Newton schemes (e.g., L-BFGS [43]). Consequently, an efficient procedure for computing gradients with respect to θ is essential in the DEQ framework, whose infinite depth precludes the direct application of conventional backpropagation. A well-known method to overcome this issue is to rely on implicit differentiation (see [31]). For clarity of exposition, we consider a single pair consisting of an observation and its ground truth, denoted by (y, x∗ ). Denote the loss by ℓθ := ℓ x∞ θ; (x0 , y) , x∗ and write x∞ θ; (x0 , yi ) as shorthand for x∞ . Since x∞ satisfies the fixed-point relation x∞ = fθ (x∞ ; y), we can treat x∞ as an implicit function of θ. By the chain rule, the gradient of the loss with respect to the network parameters is given by ∞ T ∂x ∂ℓθ ∂ℓθ = . (12) ∂θ ∂θ ∂x∞ ∂ℓθ For commonly used losses such as the mean squared error, the term ∂x ∞ admits a simple closed-form ∞ expression. The main computational challenge therefore lies in evaluating ∂x ∂θ . This difficulty can be substantially simplified by exploiting the fixed-point relation satisfied by the equilibrium state. Since x∞ is a fixed point of fθ (·; y), it satisfies

x∞ = fθ (x∞ ; y). Differentiating both sides with respect to θ and rearranging the resulting expression yields −1 ∂x∞ ∂fθ ∂fθ (x∞ ; y) = I− , ∂θ ∂x x=x∞ ∂θ ∞

;y) where ∂fθ (x denotes the Jacobian of fθ with respect to θ, evaluated at the equilibrium point x∞ . ∂θ Substituting this identity into the chain rule (12) leads to the final gradient formula −T ∂ℓθ ∂ℓθ ∂fθ (x∞ ; y) T ∂fθ = I− . ∂θ ∂θ ∂x x=x∞ ∂x∞ ∞

;y) The Jacobian ∂fθ (x can be computed via automatic differentiation, once the fixed point x∞ is ∂θ obtained using Algorithm 1 for a fixed θ. Consequently, the only remaining task is the evaluation −T θ of I − ∂f . Rather than solving the associated linear system explicitly, as in [8, 31], we ∂x ∞ x=x

adopt a Jacobian-free backpropagation (JFB) strategy [29]. This approach avoids the computational overhead of large-scale linear system solves, requires only fixed memory, and has been observed to provide competitive gradient accuracy in practice. −T θ At an implementation level, JFB approximates the inverse-transpose Jacobian term I − ∂f ∂x x=x∞ by the identity operator I. Although this approximation may appear heuristic, it admits several rigorous interpretations. One viewpoint is to regard the equilibrium state fθ (x∞ ; y) as independent of

A geometry-based deep equilibrium model

the network input x∞ during differentiation, effectively treating the fixed point as if it were obtained after a single converged iteration. For this reason, JFB is often referred to as one-step differentiation [16]. An attractive consequence of this interpretation is that the resulting gradient can be computed directly within standard automatic differentiation frameworks. From an optimization perspective, JFB can also be interpreted as a preconditioned gradient method in which the exact gradient is replaced by an alternative descent direction obtained through a suitable preconditioning strategy; see [29] for a detailed analysis. Another interpretation views JFB as the zeroth-order truncation of the Neumannseries expansion of the inverse Jacobian term, as discussed in [23]. These theoretical justifications, together with its favorable computational efficiency, have made JFB a widely adopted approximation strategy in the training of equilibrium models [68, 30].

Experimental results

To quantitatively assess reconstruction quality, we first introduce three widely used quantitative metrics for evaluating reconstruction quality, each characterizing a different aspect of image fidelity. • PSNR (Peak Signal-to-Noise Ratio) [36]: Measures the pixel-wise mean squared error between the reconstructed image and the ground truth, expressed in decibels, reflecting fidelity at the pixel level. • SSIM (Structural Similarity Index) [56]: Computes a weighted combination of luminance, contrast, and structural components, yielding a value in [0, 1] that quantifies structural similarity. • LPIPS (Learned Perceptual Image Patch Similarity) [63]: Based on a pretrained deep neural network (AlexNet in our experiments), computes a distance in feature space between the reconstruction and the ground truth. This metric, typically nonnegative, correlates well with human perceptual judgment.

5.1

Training details

The proposed model is trained on the DIV2K dataset [2], which contains a diverse collection of highresolution natural images. The 720 training images and 80 validation images are cropped to a size of (256×256) before training. All experiments are conducted on a workstation equipped with an Intel Core i9-12900K CPU, 125 GB RAM, and an NVIDIA RTX A5000 GPU with 24 GB memory, running Ubuntu 24.04.4 LTS.

5.2

Parameter and function initialization

The area-related influence function is initialized by the smooth approximation ψ a (z) = log(1 + z 2 ), while the curvature-related influence function is initialized as the quadratic penalty ψ c (z) = z 2 . The corresponding RBF coefficients are obtained by least-squares fitting, where the inputs are normalized and rescaled to ensure that they lie within the effective support of the basis functions. The implementation details are available in the accompanying code repository https: // github. com/ Shengkunyang/ DEQ4Gamma . To ensure positivity, we parameterize ϵ ← exp(ϵ̃) and b = exp(b̃), with ϵ̃ and b̃ initialized such that ϵ = 10−6 and b = 0.01. Since the learned influence functions are expected to adapt the relative scaling between the data fidelity and regularization terms, the balancing parameter λ is fixed to 1 throughout all experiments. Given a degraded observation y, the initial estimate of the DEQ model is chosen as x0 = A⊤ y. Each convolution kernel is initialized as a random linear combination of DCT basis filters. During training, the parameters θ are learned by minimizing the mean squared error loss L(θ). For grayscale image restoration, the network employs 12 learnable 5 × 5 convolution kernels, each parameterized by 24 DCT basis filters (Nr = 24), together with 33 Gaussian

A geometry-based deep equilibrium model

radial basis functions (M = 33). The regularizer contains a total of 13 learnable influence functions (Na + Nc = 13). For color image restoration, to reduce GPU memory consumption, the numbers of area- and curvature-related influence functions are set to Na = 8 and Nc = 3, respectively, while all other network configurations remain unchanged. The resulting parameterized architecture is illustrated in Fig. 1.

Figure 1: Overall architecture of the proposed Mirror Descent-based Deep Equilibrium (DEQ) model. (a) Overall pipeline of the proposed framework. (b) Mirror Descent module for equilibrium computation. (c) Geometry-aware learning architecture for estimating the regularization operator.

5.3

Gray image restoration

We evaluate the proposed method on grayscale image restoration under two blur settings: a Gaussian blur with standard deviation 2 and a motion blur with length 5 and angle 30◦ . After blurring, the images are contaminated with multiplicative Gamma noise with equivalent numbers of looks L = 4 and L = 10. To the best of our knowledge, although a deep learning method has recently been proposed for joint deblurring and multiplicative Gamma denoising [41], no publicly available implementation is currently available, preventing a direct experimental comparison. Fortunately, the literature provides a number of well-established model-based methods for this challenging restoration task. Accordingly, we compare the proposed method with five representative model-based approaches, namely the AA [6], RLO

A geometry-based deep equilibrium model

[50], DZ [25], ZWN [65], and MG [57] models. For each approach, the regularization parameters are tuned independently for each image to achieve the best reconstruction performance. For the proposed method, the variants employing the KL-type and AA-type data fidelity terms are denoted by DEQKL and DEQ-AA, respectively. The same notation is adopted throughout the color image restoration experiments. The experiments are performed on three representative grayscale images: the synthetic image “Dart”, the SAR image “SAR”, and the natural image “Leaves”, which are shown in Figure 2.

(a) Dart

(b) SAR

Figure 2: Three test images for grayscale image restoration experiments.

(a) Blurred & Noisy

(b) AA

(d) DZ

(e) ZWN

(f) MG

(g) DEQ-KL

(h) DEQ-AA

Figure 3: Comparison of restoration results obtained by different methods on the Dart image degraded by Gaussian blur and noise with noise level L = 4. Figures 3-6 present representative restoration results on the grayscale test images. Overall, the proposed method consistently produces superior visual quality compared with the variational-based baselines, particularly for SAR and natural images. For example, in Fig. 4, our method recovers sharper road structures while effectively suppressing isolated dark artifacts. Similarly, the fine leaf textures are more faithfully reconstructed in Fig. 5. For the synthetic Dart image, however, variational methods

A geometry-based deep equilibrium model

(a) Blurred & Noisy

(b) AA

(d) DZ

(e) ZWN

(f) MG

(g) DEQ-KL

(h) DEQ-AA

Figure 4: Comparison of restoration results obtained by different methods on the SAR image degraded by Gaussian blur and noise with noise level L = 4.

(a) Blurred & Noisy

(b) AA

(d) DZ

(e) ZWN

(f) MG

(g) DEQ-KL

(h) DEQ-AA

Figure 5: Comparison of restoration results obtained by different methods on the Leaves image degraded by Gaussian blur and noise with noise level L = 10.

A geometry-based deep equilibrium model

Table 1: Comparison of PSNR, SSIM, and LPIPS achieved by different restoration methods on grayscale test images degraded by Gaussian blur and noise under various noise levels. Image

Method

PSNR↑

SSIM↑

LPIPS↓

AA RLO DZ ZWN MG DEQ-KL DEQ-AA

13.52 13.56 14.49 17.09 15.47 16.05 16.07

0.49 0.44 0.50 0.78 0.70 0.79 0.76

0.44 0.51 0.52 0.16 0.21 0.15 0.14

AA RLO DZ ZWN MG DEQ-KL DEQ-AA

14.97 14.97 15.15 17.34 16.19 18.95 18.80

0.51 0.52 0.53 0.88 0.73 0.86 0.85

0.41 0.39 0.47 0.11 0.21 0.09 0.09

AA RLO DZ ZWN MG DEQ-KL DEQ-AA

19.41 19.59 20.27 20.45 20.28 20.65 21.45

0.46 0.48 0.57 0.61 0.55 0.57 0.64

0.36 0.38 0.33 0.25 0.32 0.32 0.23

AA RLO DZ ZWN MG DEQ-KL DEQ-AA

20.92 20.98 21.18 21.79 21.54 22.07 22.69

0.58 0.58 0.60 0.66 0.63 0.70 0.72

0.34 0.33 0.30 0.23 0.28 0.21 0.18

AA RLO DZ ZWN MG DEQ-KL DEQ-AA

16.46 16.47 16.93 17.57 17.62 17.81 18.07

0.44 0.43 0.46 0.56 0.57 0.60 0.61

0.51 0.57 0.54 0.40 0.36 0.37 0.36

AA RLO DZ ZWN MG DEQ-KL DEQ-AA

17.93 17.77 18.14 17.87 18.47 19.91 20.02

0.51 0.53 0.54 0.65 0.62 0.71 0.72

0.54 0.46 0.47 0.30 0.31 0.28 0.24

Dart

SAR

Leaves

A geometry-based deep equilibrium model

Table 2: Comparison of PSNR, SSIM, and LPIPS achieved by different restoration methods on grayscale test images degraded by motion blur and noise under various noise levels. Image

Method

PSNR↑

SSIM↑

LPIPS↓

AA RLO DZ ZWN MG DEQ-KL DEQ-AA

12.13 14.09 15.22 17.20 16.35 17.72 17.64

0.49 0.51 0.56 0.69 0.68 0.83 0.83

0.51 0.38 0.40 0.22 0.22 0.17 0.17

AA RLO DZ ZWN MG DEQ-KL DEQ-AA

15.75 15.88 17.37 17.97 16.87 18.17 18.77

0.59 0.56 0.79 0.88 0.78 0.83 0.85

0.32 0.34 0.17 0.11 0.17 0.14 0.12

AA RLO DZ ZWN MG DEQ-KL DEQ-AA

20.03 19.94 20.69 20.52 21.33 20.77 21.53

0.53 0.51 0.62 0.63 0.63 0.62 0.65

0.37 0.37 0.31 0.25 0.26 0.39 0.28

AA RLO DZ ZWN MG DEQ-KL DEQ-AA

21.57 21.64 22.02 22.63 22.34 22.69 23.16

0.63 0.63 0.66 0.72 0.69 0.72 0.74

0.28 0.29 0.28 0.21 0.24 0.25 0.20

AA RLO DZ ZWN MG DEQ-KL DEQ-AA

16.42 16.49 16.94 17.51 17.84 18.30 18.40

0.44 0.45 0.48 0.53 0.59 0.63 0.64

0.56 0.52 0.59 0.47 0.34 0.41 0.38

AA RLO DZ ZWN MG DEQ-KL DEQ-AA

17.32 18.41 18.60 18.44 19.16 19.77 20.01

0.55 0.56 0.58 0.68 0.67 0.67 0.73

0.47 0.45 0.43 0.34 0.28 0.34 0.30

Dart

SAR

Leaves

A geometry-based deep equilibrium model

(a) Blurred & Noisy

(b) AA

(d) DZ

(e) ZWN

(f) MG

(g) DEQ-KL

(h) DEQ-AA

Figure 6: Comparison of restoration results obtained by different methods on the Dart image degraded by motion blur and noise with noise level L = 10. remain competitive under severe degradations, especially the ZWN model. As shown in Fig. 3, ZWN produces higher-contrast dartboard patterns, whereas our method achieves more effective deblurring. Under a lower noise level, the advantage of the proposed method becomes more pronounced. As illustrated in Fig. 6, our reconstruction preserves sharper cross patterns in the four corners as well as finer structures around the dartboard center. Quantitative comparisons are reported in Tables 1 and 2 for Gaussian and motion blur, respectively, in terms of PSNR, SSIM, and LPIPS, where the best result in each setting is highlighted in bold. For the synthetic Dart image with Gaussian blur and noise level L = 4, the proposed method achieves a slightly lower PSNR than the ZWN model, while yielding superior SSIM and LPIPS values, indicating a better balance between structural fidelity and perceptual quality. Apart from this single case, our method consistently outperforms all competing variational models on the remaining grayscale test images. In particular, for the challenging natural image Leaves, the proposed method improves the PSNR by at least 0.46 dB over the best competing baseline under all considered degradation settings. Finally, DEQ-AA consistently achieves slightly better performance than DEQ-KL across the test set, as reflected, for example, by the quantitative results on Leaves in Tables 1 and 2. A possible explanation is that the joint deblurring and Gamma denoising task considered in this work exhibits a more complex degradation model, under which the AA-type fidelity provides a more accurate statistical characterization of Gamma noise than the KL-type fidelity. A deeper theoretical and experimental investigation of this observation is left for future work.

A geometry-based deep equilibrium model

5.4

Color image restoration

We evaluate the proposed method on a diverse collection of color images, including 3 images from Set3C 1 , 24 from Kodak24 2 , 17 from BSD500 [3], and 11 from Set14 [60]. Degraded observations are generated by first applying either an 15 × 15 Gaussian blur kernel with standard deviation 2 or a realistic motion blur kernel [39], followed by multiplicative Gamma noise with equivalent numbers of looks L = 4 and L = 10. Similar to the grayscale restoration setting, a method for joint deblurring and multiplicative Gamma denoising of color images has been proposed in [55]. Unfortunately, its source code is not publicly available, preventing us from including it in our experimental comparison. As a representative learningbased baseline, we therefore adopt the DEQ framework of [23]. Although originally designed for Poisson image restoration, its KL-type data-fidelity model has also been successfully applied to multiplicative Gamma noise removal [52, 64], and coincides with the fidelity term employed in our DEQ-KL model. To ensure a fair comparison, we regenerate the training data according to the degradation model adopted in this work and retrain the network of [23] using the corresponding blurred Gamma-noise observations. The retrained network follows the DEQ framework and adopts a Regularization by Denoising (RED) parameterization; accordingly, we denote it by DEQ-RED. The principal difference between DEQ-RED and DEQ-KL lies in the regularization model: DEQ-RED employs an implicit CNN regularizer, whereas DEQ-KL incorporates an explicit geometry-aware regularizer together with the proposed positivity-preserving backtracking strategy. Table 3: Average quantitative results (PSNR, SSIM, and LPIPS) on the color test set under different blur types and noise levels. L

Method

PSNR↑

SSIM↑

LPIPS↓

Gaussian Blur

DEQ-RED [23]

22.87

0.58

0.55

DEQ-KL

23.09

0.59

0.48

DEQ-AA

23.17

0.60

0.47

DEQ-RED [23]

23.67

0.60

0.51

DEQ-KL

23.79

0.61

0.44

DEQ-AA

24.03

0.63

0.41

Motion Blur

DEQ-RED [23]

20.55

0.46

0.63

DEQ-KL

20.62

0.48

0.60

DEQ-AA

20.99

0.49

0.55

DEQ-RED [23]

21.08

0.49

0.58

DEQ-KL

21.92

0.54

0.49

DEQ-AA

21.97

0.54

0.48

Table 3 reports the quantitative results of the proposed methods and the competing approaches on the color image test set, where the best performance in each setting is highlighted in bold. Under both 1 2

https://huggingface.co/datasets/deepinv/set3c https://www.kaggle.com/datasets/drxinchengzhu/kodak24

A geometry-based deep equilibrium model

(a) Clean

(b) Blurred & Noisy

(PSNR↑, SSIM↑, LPIPS↓)

(8.25,0.10,1.25)

(d) DEQ-KL

(e) DEQ-AA

(22.17, 0.75, 0.36)

(22.25, 0.78, 0.24)

(22.26, 0.79, 0.24)

Figure 7: Comparison of restoration results obtained by different methods on the Magazine image degraded by Gaussian blur and noise with noise level L = 4. Gaussian and motion blur degradations, DEQ-KL achieves performance comparable to that of DEQRED, indicating that incorporating an explicit variational regularizer does not compromise restoration quality relative to implicit regularization. A further advantage of the proposed framework is its parameter efficiency. Owing to the incorporation of explicit geometric priors, our model contains only 8 × 24 + 8 × 33 + 3 × 33 + 2 ≈ 600 trainable parameters, whereas DEQ-RED requires approximately 105 parameters. This substantial reduction in model complexity is expected to improve training efficiency while maintaining competitive performance. Moreover, DEQ-AA consistently outperforms DEQ-KL across most test cases, which is consistent with the observations from the grayscale image experiments. Figures 7-9 present representative visual comparisons. Under Gaussian blur, the proposed DEQAA method reconstructs sharper geometric structures and finer textures, as illustrated by the butterfly patterns in Figure 8. Another notable observation is that DEQ-RED, which relies on an implicit neural prior for image reconstruction, occasionally introduces visually noticeable artifacts. This phenomenon is also reflected by its relatively inferior LPIPS scores reported in Table 3. For example, although DEQ-RED and the proposed method achieve comparable PSNR values in the Gaussian blur experiment shown in Figure 7, DEQ-RED produces visible shadow-like artifacts that degrade perceptual quality. The difference becomes more pronounced under the more challenging real motion blur degradation. As illustrated in Figure 9, DEQ-RED generates noticeable blotchy artifacts in the Peppers image, whereas the proposed method preserves cleaner structures and more natural visual appearance. We attribute this behavior to the explicit geometry-aware variational regularization, which provides stronger structural priors and

A geometry-based deep equilibrium model

(a) Clean

(b) Blurred & Noisy

(PSNR↑, SSIM↑, LPIPS↓)

(11.47,0.23,0.79)

(d) DEQ-KL

(e) DEQ-AA

(21.53, 0.72, 0.27)

(21.62, 0.74, 0.23)

(21.74, 0.75, 0.18)

Figure 8: Comparison of restoration results obtained by different methods on the Butterfly image degraded by Gaussian blur and noise with noise level L = 4. effectively suppresses the artifacts that may arise in purely data-driven restoration models.

5.5

On the convergence of DEQ models at inference

In all numerical experiments, the iterative solver is terminated once the relative error satisfies the stopping criterion ∥xk+1 − xk ∥ ≤ ϵS . ∥xk+1 ∥ To better illustrate the convergence behavior of the proposed method, we adopt a stringent stopping tolerance of ϵS = 10−6 in this subsection. In all other experiments, we use the empirical choice ϵS = 10−5 , which provides a good compromise between reconstruction accuracy and computational efficiency. To illustrate the convergence behavior of the proposed DEQ models, we report the evolution of the objective value Ψ and the PSNR for image Magazine under Gaussian blur with Gamma noise of level L = 4; the corresponding results for DEQ-KL and DEQ-AA are shown in Fig. 10. It can be observed that the objective values of both models decrease monotonically and eventually stabilize. In terms of reconstruction quality, the PSNR of DEQ-AA reaches a stable value after approximately 700 iterations, whereas DEQ-KL requires about 1000 iterations to converge. This observation indicates that the minimization of the proposed objective functional is well aligned with improvements in reconstruction quality, suggesting that convergence of the optimization process is accompanied by convergence to a

A geometry-based deep equilibrium model

(a) Clean

(b) Blurred & Noisy

(PSNR↑, SSIM↑, LPIPS↓)

(14.62,0.09,0.94)

(d) DEQ-KL

(e) DEQ-AA

(24.11, 0.53, 0.43)

(24.12, 0.55, 0.36)

(24.52, 0.60, 0.30)

Figure 9: Comparison of restoration results obtained by different methods on the Peppers image degraded by motion blur and noise with noise level L = 10. high-quality reconstruction. Furthermore, since a backtracking strategy is employed, the step size is adaptively adjusted during the iterations, which leads to discontinuities in the convergence curves of both DEQ-KL and DEQ-AA. The iterations at which the step size is updated are indicated by red vertical lines in Fig. 10. It can be observed that the objective energy continues to decrease and the PSNR keeps increasing after these adjustments, demonstrating the effectiveness of the adopted backtracking strategy.

5.6

Learned filters and characteristic influence functions in neural networks

One of the main contributions of the developed DEQ model is the possibility to explain the reconstructions as critical points of an energy functional with interpretable learned regularization term (2). To gain further insight into the learned image priors, we visualize representative convolution filters and influence functions in this subsection. Fig. 11 compares the convolution filters before and after training, where the color intensity indicates the magnitude and sign of the filter coefficients. Although the randomly initialized filters exhibit no discernible structure, the learned filters evolve into organized patterns that capture meaningful image P a local a (k ∗ y), features. To further reveal the learned representation, Fig. 12 visualizes the responses N ψ l l=1 l obtained by applying the learned filters and influence functions to the degraded observation y of the grayscale Parrot image corrupted by Gaussian blur and multiplicative Gamma noise with L = 4. As training progresses, these responses become progressively closer to the corresponding ground-truth im-

A geometry-based deep equilibrium model

(a) DEQ-KL

(b) DEQ-AA

Figure 10: Convergence verification of (a) DEQ-KL and (b) DEQ-AA methods, showing PSNR and energy variations with iterations. age, indicating that the learned representation effectively captures the underlying image structure. Initialization

After training

Figure 11: Representative convolution filters before and after training. This observation is consistent with the purpose of the area-related component in the proposed regularizer. Since the variance of multiplicative Gamma noise is proportional to the underlying signal intensity, PNa a brighter regions are generally subject to stronger noise corruption. Therefore, the quantity l=1 ψl (kl ∗ y) is expected to provide a rough estimate of the latent clean image, enabling the strength of the area-related regularization to be adapted according to the local imagecontent. A related strategy p Gσ ∗y was adopted in [64, 57], where the normalized Gaussian-smoothed image, max Gσ ∗y , was employed as a gray-level indicator. As illustrated in Fig. 12, the learned gray-level indicator provides a substantially more accurate approximation of the clean image than the fixed Gaussian indicator. These results demonstrate the benefit of learning the convolution filters and influence functions directly from data, thereby providing further support for the proposed extension of the model in [57]. The learned influence functions subsequently perform nonlinear modulation of these extracted local features. As shown in Fig. 13(a), the curvature-related penalty function ψ c consistently attains its minimum around zero curvature, indicating that the proposed model favors locally flat or straight image structures while suppressing unnecessary geometric oscillations. This observation is consistent with the

A geometry-based deep equilibrium model

Clean

Indicator [64]

Epoch 0

Epoch 5

Final epoch

Recovery

Figure 12: Evolution of the learned area-related influence function throughout training and the corresponding reconstruction From left to right: the ground-truth Parrot image; the performance. p Gσ ∗y from [64], where the recommended parameters σ = 0.8, p = 1.2 are gray-level indicator max Gσ ∗y P a a adopted; the learned area-related indicator N l=1 ψl (kl ∗ y) at epochs 0, 5, and the final epoch; and the reconstruction produced by the final learned regularizer. geometric interpretation of curvature regularization reported in [67], where zero curvature corresponds to locally smooth structures. In contrast, the influence functions ψ a associated with the learned convolutional filters exhibit substantially richer behaviors, as illustrated in Figs. 13(b)-13(d). Besides the conventional unimodal penalties centered at zero, double-well and multi-well potentials are frequently observed, which have been widely adopted in image restoration models [13, 42]. This suggests that different filters favor different response magnitudes rather than uniformly penalizing all nonzero responses. Since different filters capture distinct local image characteristics, the corresponding learned penalties naturally adapt to these heterogeneous statistics. Moreover, several learned influence functions are asymmetric, implying that positive and negative filter responses are treated differently and thereby providing greater flexibility for modeling directional image structures. We conduct an additional experiment to investigate the dependence of the learned regularizer Rθ on the training data. Specifically, we train the network on a dataset of piecewise-constant images with smooth boundaries, for which we expect the curvature component of Rθ to play a less significant role. The dataset consists of 300 images containing ellipses with random positions, aspect ratios, orientations, and intensities, generated using the Operator Discretization Library [1]. Taking the AA-type fidelity as an example, we found that the learned area-related influence functions consistently exhibit the pattern shown as Activation II in Fig. 13, whereas the curvature-related influence functions remain nearly identical to their quadratic initialization throughout training (Fig. 14(a)). Meanwhile, the learned parameter b gradually decreases during training (Fig. 14(b)). These results indicate that, for piecewise-constant images, the proposed framework naturally emphasizes the arearelated regularization and learns a regularizer resembling the classical TV model.

(a) Activation I

(b) Activation II

(d) Activation IV

Figure 13: Representative learned influence functions. (a) The learned curvature influence function ψ c . (b)–(d) Representative learned activation influence functions ψ a corresponding to different convolutional filters, illustrating truncated convex, double-well, and multiple-well structures.

A geometry-based deep equilibrium model

(a) Curvature-aware influence function

(b) Parameter b evolution

Figure 14: Ellipses dataset results. (a) Learned curvature-aware influence function. (b) Evolution of parameter b across training epochs.

Conclusions

In this work, we proposed a deep equilibrium framework for joint image deblurring and multiplicative Gamma noise removal. By explicitly incorporating geometry-aware regularization with learnable filters and influence functions, the proposed method provides an interpretable balance between model-based and data-driven approaches. Furthermore, based on a suitable Bregman mirror descent scheme and the o-minimal framework, the convergence of the proposed algorithm to a critical point is established. Extensive experiments on grayscale and color images demonstrate the effectiveness of the proposed method compared with existing variational-based, and deep learning approaches.

Bregman geometry and proximal operators

In this section, we recall the definitions of the fundamental concepts from nonconvex optimization used throughout Sections 2.2 and 3. Definition 2 (Kernel generating distance [7]). Let Ω ⊆ Rn be a nonempty, open, and convex set. A function ϕ : Rn → (−∞, +∞] is called a kernel generating distance associated with Ω if it satisfies the following conditions: (i) ϕ is proper, lower semicontinuous, and convex. Moreover, dom ϕ ⊆ Ω, int(dom ϕ) = Ω. (ii) ϕ is continuously differentiable on Ω, i.e., ϕ ∈ C 1 (Ω). We denote the class of kernel generating distances by G(Ω). Definition 3. [48] Let h : X ⊂ Rn → (−∞, +∞] be a proper lower semicontinuous function. (i) (essential smoothness). A differentiable function h on int(dom(h)) is said to be essentially smooth if, for every sequence {xk } ⊂ int(dom(h)) converging to a boundary point of dom(h), we have ∥∇h(xk )∥ → +∞. (ii) (Legendre function). The function h is called a Legendre function if it is essentially smooth and strictly convex on int(dom(h)).

A geometry-based deep equilibrium model

Definition 4. [19]. Given a Legendre function h ∈ G(Ω) : X ⊂ Rn → (−∞, +∞], the associated Bregman distance between x ∈ dom h and y ∈ int(dom(h)) is defined by Dh (x, y) = h(x) − h(y) − ⟨∇h(y), x − y⟩.

Remark 8. Owing to the convexity of h, the Bregman distance is strictly convex with respect to its first argument. Moreover, Dh (x, y) ≥ 0 for all (x, y) ∈ dom(h) × int(dom(h)), with equality if and only if x = y. Consequently, Dh serves as a generalized proximity measure. Unlike the Euclidean distance, however, it is generally non-symmetric. The special choice h = ∥ · ∥2 yields the classical squared Euclidean distance. Definition 5. Let R : X ⊂ Rn → (−∞, +∞] be a proper lower semicontinuous function and h : X ⊂ Rn → (−∞, +∞] be a Legendre function. The Bregman proximal operator associated with R (w.r.t. h) is defined by ProxhτR (x) = arg min λR(u) + Dh (u, x), x ∈ int(dom(h)), u∈X

where λ > 0 is the step-size. Remark 9. Since R may be nonconvex, the corresponding proximal mapping ProxhτR is generally setvalued. If, in addition, R is convex and coercive, then the Bregman proximal operator is single-valued [23].

Kurdyka-Lojasiewicz functions

This section specifically focuses on the Kurdyka-Lojasiewicz (KL), which plays a central role in the convergence analysis of nonconvex optimization algorithms. We first recall its definition and then discuss its verification in more general settings, such as the one considered in this project, where simpler tools based on global subanalyticity cannot be directly applied. Definition 6. Let Φ : Rn → (−∞, +∞] be a proper lower semicontinuous function. (a) The function Φ is said to satisfy the KL property at a point x̄ ∈ dom Φ if there exist η ∈ (0, +∞], a neighborhood U of x̄, and a continuous concave function ψ : [0, η) → R+ such that: (i) ψ(0) = 0; (ii) ψ is continuously differentiable on (0, η) and continuous at 0; (ii) ψ ′ (t) > 0 for all t ∈ (0, η); (iv) for every x ∈ U ∩ {x ∈ Rn : Φ(x̄) < Φ(x) < Φ(x̄) + η}, the KL inequality ψ ′ (Φ(x) − Φ(x̄)) dist(0, ∂Φ(x)) ≥ 1 holds, where dist(0, ∂Φ(x)) denotes the distance from 0 to the subgradient set ∂Φ. (b) The function Φ is called a KL function if it satisfies the KL property at every point in dom Φ. Although the KL property is generally difficult to verify directly from its definition, it is known to hold for several broad classes of functions frequently encountered in imaging and machine learning. Two particularly important examples are semialgebraic functions and functions definable in an ominimal structure, which we briefly review below.

A geometry-based deep equilibrium model

B.1

Semialgebraic functions

We next recall the class of semialgebraic sets and functions. Definition 7. (Semialgebraic sets and semialgebraic functions). • A set S ⊂ Rn is called semialgebraic if there exist finitely many polynomial functions fij , gij : Rn → R such that p \ q [ S= {x ∈ Rn : fij (x) = 0, gij (x) < 0}. i=1 j=1

• A function f : Rn → R ∪ {+∞} is called semialgebraic if its graph {(x, t) ∈ Rn+1 : f (x) = t} is a semialgebraic subset of Rn+1 . The class of semialgebraic sets is closed under finite unions, finite intersections, complements, and Cartesian products. Correspondingly, semialgebraic functions form a rich and flexible class, which is closed under finite sums, products, and compositions. Important examples include real polynomial functions and indicator functions of semialgebraic sets.

B.2

Functions definable in an o-minimal structure

We then recall the notion of an o-minimal structure, introduced by Lou van den Dries [26], which provides a powerful framework for analyzing functions arising in optimization; see also [15]. Definition 8. (o-minimal structure) Let M = {Mn }n≥1 be a sequence such that each Mn is a family of subsets of Rn . We say that M is an o-minimal structure on the field (R, +, ·) if it satisfies the following conditions: (i) For every n ≥ 1, the collection Mn is a Boolean algebra of subsets of Rn . Equivalently, it contains ∅ and Rn , and is closed under finite unions, finite intersections, and complements. (ii) If S ∈ Mn , then both cylinder sets S×R

and

R×S

belong to Mn+1 . (iii) Whenever T ∈ Mn+1 , its projection onto the first n coordinates, π(T ) := {(x1 , . . . , xn ) ∈ Rn : ∃ xn+1 ∈ R, (x1 , . . . , xn , xn+1 ) ∈ T } , belongs to Mn . (iv) For every polynomial function p : Rn → R, its zero set belongs to Mn ; that is, {x ∈ Rn : p(x) = 0} ∈ Mn . (v) The sets in M1 are exactly the finite unions of intervals and points. A set A is said to be definable in M if A ∈ M. This notion naturally extends to functions through their graphs, leading to the following definition of definable functions.

A geometry-based deep equilibrium model

Definition 9. (Definable functions) A function f : Rn → R ∪ {+∞} is called definable in M if its graph gph f := {(x, f (x)) : x ∈ dom f } is a definable subset of Rn+1 , i.e., gph f ∈ Mn+1 . As discussed in [4], functions definable in an o-minimal structure inherit many of the favorable properties enjoyed by semialgebraic functions. In particular, the class of definable functions is closed under several fundamental operations, including finite sums, composition, generalized inversion, and indicator functions of definable sets. The o-minimal framework encompasses a remarkably broad class of functions that commonly arise in optimization and imaging applications, including (see [26]): • Semialgebraic functions, which can be interpreted as functions definable in the o-minimal structure M(R, +, ·, (r)r∈R ), namely the smallest structure in (R, +, ·) containing all singletons. Indeed, by the Tarski-Seidenberg theorem, the collection of semialgebraic sets coincides with M(R, +, ·, (r)r∈R ). • Globally subanalytic functions [14], which are definable in the structure M(Ran ), where Ran = (R, +, ·, (f )) and f ranges over restricted analytic functions. • Log-exp functions, which are definable in the structure M(Ran,exp ), related to the field Ran,exp = (R, +, ·, (f ), exp), where f ranges over restricted analytic functions and exp(x) = ex . These definable structures satisfy the following inclusion relationship: M(R, +, ·, (r)r∈R ) ⊆ M(Ran ) ⊆ M(Ran,exp ). For a comprehensive introduction to o-minimal structures, we refer the reader to [26, 15]. Among their numerous applications in optimization theory, one of the most important is the extension of the Kurdyka-Lojasiewicz framework. Theorem 9. Suppose that f : Rn → R ∪ {+∞} is a proper lower semicontinuous function definable in an o-minimal structure M. Then, for every x ∈ dom ∂f , the function f satisfies the Kurdyka– Lojasiewicz property at x.

Acknowledgments The research of Luca Ratti was partially supported by the PNRR project M4C2-PE00000013 (FAIR), Spoke 8 “Pervasive AI”, funded by the NextGeneration EU programme (CUP J33C22002830006), and by the INDAM-GNCS 2026 project (CUP E53C25002010001). The research of Zhichang Guo was supported by the National Key Research and Development Program of China (Grant No. 2024YFE0216800).

References [1] Adler, J., Kohr, H., Öktem, O., 2017. Operator discretization library (odl). Zenodo . [2] Agustsson, E., Timofte, R., 2017. Ntire 2017 challenge on single image super-resolution: Dataset and study, in: Proceedings of the IEEE conference on computer vision and pattern recognition workshops, pp. 126–135. [3] Arbelaez, P., Maire, M., Fowlkes, C., Malik, J., 2010. Contour detection and hierarchical image segmentation. IEEE transactions on pattern analysis and machine intelligence 33, 898–916.

A geometry-based deep equilibrium model

[4] Attouch, H., Bolte, J., Redont, P., Soubeyran, A., 2010. Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the kurdyka-lojasiewicz inequality. Mathematics of operations research 35, 438–457. [5] Attouch, H., Bolte, J., Svaiter, B.F., 2013. Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward–backward splitting, and regularized gauss–seidel methods. Mathematical programming 137, 91–129. [6] Aubert, G., Aujol, J.F., 2008. A variational approach to removing multiplicative noise. SIAM journal on applied mathematics 68, 925–946. [7] Auslender, A., Teboulle, M., 2006. Interior gradient and proximal methods for convex and conic optimization. SIAM Journal on Optimization 16, 697–725. [8] Bai, S., Kolter, J.Z., Koltun, V., 2019. Deep equilibrium models. Advances in neural information processing systems 32. [9] Barnett, J., Li, W., Resmerita, E., Vese, L., 2025. Multiscale hierarchical decomposition methods for images corrupted by multiplicative noise. Journal of Mathematical Imaging and Vision 67, 12. [10] Bauschke, H.H., Bolte, J., Teboulle, M., 2017. A descent lemma beyond lipschitz gradient continuity: first-order methods revisited and applications. Mathematics of Operations Research 42, 330–348. [11] Beck, A., 2017. First-order methods in optimization. SIAM. [12] Bertero, M., Boccacci, P., Desiderà, G., Vicidomini, G., 2009. Image deblurring with poisson data: from cells to galaxies. Inverse Problems 25, 123006. [13] Bertozzi, A.L., Flenner, A., 2016. Diffuse interface models on graphs for classification of high dimensional data. siam REVIEW 58, 293–328. [14] Bolte, J., Daniilidis, A., Lewis, A., 2007a. The lojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems. SIAM Journal on Optimization 17, 1205–1223. [15] Bolte, J., Daniilidis, A., Lewis, A., Shiota, M., 2007b. Clarke subgradients of stratifiable functions. SIAM Journal on Optimization 18, 556–572. [16] Bolte, J., Pauwels, E., Vaiter, S., 2023. One-step differentiation of iterative algorithms. Advances in Neural Information Processing Systems 36, 77089–77103. [17] Bolte, J., Sabach, S., Teboulle, M., 2014. Proximal alternating linearized minimization for nonconvex and nonsmooth problems. Mathematical Programming 146, 459–494. [18] Bolte, J., Sabach, S., Teboulle, M., Vaisbourd, Y., 2018. First order methods beyond convexity and lipschitz gradient continuity with applications to quadratic inverse problems. SIAM Journal on Optimization 28, 2131–2151. [19] Bregman, L.M., 1967. The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming. USSR computational mathematics and mathematical physics 7, 200–217. [20] Chan, R.H., Ma, J., 2012. A multiplicative iterative algorithm for box-constrained penalized likelihood image restoration. IEEE transactions on image processing 21, 3168–3181.

A geometry-based deep equilibrium model

[21] Chen, Y., Pock, T., 2016. Trainable nonlinear reaction diffusion: A flexible framework for fast and effective image restoration. IEEE transactions on pattern analysis and machine intelligence 39, 1256–1272. [22] Cui, Y., Zhou, G., Yang, J., Yamaguchi, Y., 2011. Unsupervised estimation of the equivalent number of looks in sar images. IEEE Geoscience and remote sensing letters 8, 710–714. [23] Daniele, C., Villa, S., Vaiter, S., Calatroni, L., 2026. Deep equilibrium models for poisson imaging inverse problems via mirror descent. SIAM Journal on Imaging Sciences 19, 1077–1109. [24] Deledalle, C.A., Denis, L., Tabti, S., Tupin, F., 2017. Mulog, or how to apply gaussian denoisers to multi-channel sar speckle reduction? IEEE Transactions on Image Processing 26, 4389–4403. [25] Dong, Y., Zeng, T., 2013. A convex variational model for restoring blurred images with multiplicative noise. SIAM Journal on Imaging Sciences 6, 1598–1625. [26] Van den Dries, L., Miller, C., 1996. Geometric categories and o-minimal structures . [27] El Hadji, S.D., 2026. Nonlocal multiplicative denoising method based on variational approach and topological shape optimization. Computers & Mathematics with Applications 204, 1–23. [28] Feng, X., Zhu, X., 2021. Models for multiplicative noise removal, in: Handbook of Mathematical Models and Algorithms in Computer Vision and Imaging: Mathematical Imaging and Vision. Springer, pp. 1–34. [29] Fung, S.W., Heaton, H., Li, Q., McKenzie, D., Osher, S., Yin, W., 2022. Jfb: Jacobian-free backpropagation for implicit networks, in: Proceedings of the AAAI Conference on Artificial Intelligence, pp. 6648–6656. [30] Gan, W., Ying, C., Boroojeni, P.E., Wang, T., Eldeniz, C., Hu, Y., Liu, J., Chen, Y., An, H., Kamilov, U.S., 2023. Self-supervised deep equilibrium models with theoretical guarantees and applications to mri reconstruction. IEEE Transactions on Computational Imaging 9, 796–807. [31] Gilton, D., Ongie, G., Willett, R., 2021. Deep equilibrium architectures for inverse problems in imaging. IEEE Transactions on Computational Imaging 7, 1123–1133. [32] Gkillas, A., Ampeliotis, D., Berberidis, K., 2023a. A highly interpretable deep equilibrium network for hyperspectral image deconvolution, in: ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE. pp. 1–5. [33] Gkillas, A., Ampeliotis, D., Berberidis, K., 2023b. An optimization-based deep equilibrium model for hyperspectral image deconvolution with convergence guarantees. arXiv preprint arXiv:2306.06378 . [34] Huang, Y.M., Ng, M.K., Wen, Y.W., 2009. A new total variation method for multiplicative noise removal. SIAM Journal on imaging sciences 2, 20–40. [35] Hurault, S., Kamilov, U., Leclaire, A., Papadakis, N., 2023. Convergent bregman plug-andplay image restoration for poisson inverse problems. Advances in Neural Information Processing Systems 36, 27251–27280. [36] Jähne, B., 2005. Digital image processing. Springer. [37] Kadkhodaie, Z., Simoncelli, E., 2021. Stochastic solutions for linear inverse problems using the prior implicit in a denoiser. Advances in Neural Information Processing Systems 34, 13242– 13254.

A geometry-based deep equilibrium model

[38] Kingma, D.P., Ba, J., 2015. Adam: A method for stochastic optimization. International Conference on Learning Representations (ICLR) . [39] Levin, A., Weiss, Y., Durand, F., Freeman, W.T., 2009. Understanding and evaluating blind deconvolution algorithms, in: 2009 IEEE conference on computer vision and pattern recognition, IEEE. pp. 1964–1971. [40] Li, Y., Chang, H., Duan, Y., 2025. Curvpnp: Plug-and-play blind image restoration with deep curvature denoiser. Signal Processing 233, 109951. [41] Li, Y., Hu, J., Ni, G., Zeng, T., 2023. Deep cnn denoiser prior for blurred images restoration with multiplicative noise. Inverse Problems and Imaging 17, 726–745. [42] Li, Y., Kim, J., 2011. Multiphase image segmentation using a phase-field model. Computers & Mathematics with Applications 62, 737–745. [43] Liu, D.C., Nocedal, J., 1989. On the limited memory bfgs method for large scale optimization. Mathematical programming 45, 503–528. [44] Liu, X., Lu, J., Shen, L., Xu, C., Xu, Y., 2020. Multiplicative noise removal: nonlocal low-rank model and its proximal alternating reweighted minimization algorithm. SIAM Journal on Imaging Sciences 13, 1595–1629. [45] Muckley, M.J., Riemenschneider, B., Radmanesh, A., Kim, S., Jeong, G., Ko, J., Jun, Y., Shin, H., Hwang, D., Mostapha, M., et al., 2021. Results of the 2020 fastmri challenge for machine learning mr image reconstruction. IEEE transactions on medical imaging 40, 2306–2317. [46] Pang, Z., Guo, L., Duan, Y., Lu, J., 2018. Image restoration based on the minimized surface regularization. Computers & Mathematics with Applications 76, 1893–1905. [47] Ran, Y., Guo, Z., Shao, J., Li, Y., Wu, B., 2025. The bv solution for a nonlocal diffusion equation and its application on multiplicative noise removal. Mathematical Methods in the Applied Sciences 48, 14367–14384. [48] Rockafellar, R.T., 1997. Convex analysis. volume 28. Princeton university press. [49] Roth, S., Black, M.J., 2005. Fields of experts: A framework for learning image priors, in: 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), IEEE. pp. 860–867. [50] Rudin, L., Lions, P.L., Osher, S., 2003. Multiplicative denoising and deblurring: theory and algorithms. Springer. [51] Shi, J., Osher, S., 2008. A nonlinear inverse scale space method for a convex multiplicative noise model. SIAM Journal on imaging sciences 1, 294–321. [52] Steidl, G., Teuber, T., 2010. Removing multiplicative noise by douglas-rachford splitting methods. Journal of Mathematical Imaging and Vision 36, 168–184. [53] Ullah, A., Chen, W., Khan, M.A., Sun, H., 2017. A new variational approach for multiplicative noise and blur removal. PloS one 12, e0161787. [54] Wang, G., Pan, Z., Zhang, Z., 2019. Deep cnn denoiser prior for multiplicative noise removal. Multimedia Tools and Applications 78, 29007–29019.

A geometry-based deep equilibrium model

[55] Wang, W., Yao, M., Ng, M.K., 2021. Color image multiplicative noise and blur removal by saturation-value total variation. Applied Mathematical Modelling 90, 240–264. [56] Wang, Z., Bovik, A.C., Sheikh, H.R., Simoncelli, E.P., 2004. Image quality assessment: from error visibility to structural similarity. IEEE transactions on image processing 13, 600–612. [57] Yang, S., Guo, Z., Li, J., Song, F., Yao, W., 2025. Mixed geometry information regularization for image deblurring with multiplicative noise. Inverse Problems 41, 095005. [58] Ye, Z., Liu, K.S., Ma, T., Gao, J., Chen, C., 2019. Curvature graph network, in: International conference on learning representations. [59] Zeng, C., Wu, C., 2019. On the discontinuity of images recovered by noncovex nonsmooth regularized isotropic models with box constraints. Advances in Computational Mathematics 45. [60] Zeyde, R., Elad, M., Protter, M., 2010. On single image scale-up using sparse-representations, in: International conference on curves and surfaces, Springer. pp. 711–730. [61] Zhang, J., Ghanem, B., 2018. Ista-net: Interpretable optimization-inspired deep network for image compressive sensing, in: Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 1828–1837. [62] Zhang, K., Li, Y., Zuo, W., Zhang, L., Van Gool, L., Timofte, R., 2021. Plug-and-play image restoration with deep denoiser prior. IEEE Transactions on Pattern Analysis and Machine Intelligence 44, 6360–6376. [63] Zhang, R., Isola, P., Efros, A.A., Shechtman, E., Wang, O., 2018. The unreasonable effectiveness of deep features as a perceptual metric, in: Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 586–595. [64] Zhang, Y., Li, S., Guo, Z., Wu, B., Du, S., 2022. Image multiplicative denoising using adaptive euler’s elastica as the regularization. Journal of Scientific Computing 90, 69. [65] Zhao, X.L., Wang, F., Ng, M.K., 2014. A new convex optimization model for multiplicative noise and blur removal. SIAM Journal on Imaging Sciences 7, 456–475. [66] Zhu, W., 2020. Image denoising using l p-norm of mean curvature of image surface. Journal of Scientific Computing 83, 32. [67] Zhu, W., Chan, T., 2012. Image denoising using mean curvature of image surface. SIAM Journal on Imaging Sciences 5, 1–32. [68] Zou, Z., Liu, J., Wohlberg, B., Kamilov, U.S., 2023. Deep equilibrium learning of explicit regularization functionals for imaging inverse problems. IEEE Open Journal of Signal Processing 4, 390–398.

A geometry-based deep equilibrium model for image restoration under multiplicative Gamma noise

Analyze & explore related research in Conceptio

Related documents