domain: phys math

text:
Stephan De Bièvre
email: stephan.de-bievre@univ-lille1.fr

Carlos Mejía-Monasterio
email: carlos.mejia@upm.es

Paul E Parris
email: ‡parris@mst.edu

Dynamical mechanisms leading to equilibration in two-component gases

Keywords: numbers: 05.70.-a, 02.50.Ey, 05.20.-y

Demonstrating how microscopic dynamics cause large systems to approach thermal equilibrium remains an elusive, longstanding, and actively-pursued goal of statistical mechanics. We identify here a dynamical mechanism for thermalization in a general class of two-component dynamical Lorentz gases, and prove that each component, even when maintained in a non-equilibrium state itself, can drive the other to a thermal state with a well-defined effective temperature.

That isolated systems with many degrees of freedom evolve asymptotically in time towards thermal equilibrium lies at the heart of classical thermodynamics. Statistical mechanics teaches that for systems described by a Hamiltonian H, the thermal states are those described by the canonical Boltzmann relation ρ = Z -1 exp(-βH). This follows from original arguments of Maxwell, marginal distributions that arise from microcanonical ensembles, and the properties of the maximum entropy states to which systems thermodynamically tend [START_REF] Khinchin | Mathematical Foundations of Statistical Mechanics[END_REF].

While such statistical arguments identify the thermal state, they provide no insight into the problem of how the microscopic dynamics of diverse large systems each lead towards equilibrium from an arbitrary initial state.

Considerable progress has been made in understanding return to equilibrium [START_REF] Bach | [END_REF] in which a small system coupled to a large thermal reservoir thermalizes to the temperature of the latter. The more general problem of approach to equilibrium, in which many mutually interacting elements of an isolated system are initially out of equilibrium, is less understood [START_REF] Fermi | Studies of nonlinear problems I[END_REF]. In spite of recent progress, including identification of the thermalization mechanism in the Fermi-Pasta-Ulam problem [4], there is continuing discussion and even some controversy regarding the role played by various properties of the dynamics, such as chaos, mixing, and resonances [5,6], as well as on the emergence of non-equilibrium effective temperatures [7].

As an important step in this longstanding problem, we identify here a dynamical mechanism for thermalization in a general class of two component systems [8][9][10][11][12][13][14][15][16], in which non-interacting point particles move freely through a spatially fixed array of isolated dynamical "scatterer" particles, each of which has a few rotationally-invariant degrees of freedom with which the itinerant point particles interact and exchange energy when they are within a fixed range.

In particular, we prove here that, starting from an arbitrary initial state, an ensemble of probe subsystems (either the itinerant point particles or the "scatterer" particles) when subjected to repeated weak interactions with members of the complementary component gas, acting as a homogeneous and stationary energy reservoir, approaches a Boltzmann state at a well-defined temperature, provided merely that the reservoir is stationary-it need not be in thermal equilibrium.

While providing an understanding of the approach to equilibrium of the entire system as a whole, the result also provides insight into and allows a justification of the concept of non-equilibrium effective temperatures.

In systems of this type studied previously, the scatterers were taken to be disks [8][9][10][11]17] or needles in 2D [12] that rotate about fixed centers, and harmonic oscillators [13][14][15][16]. In each of these systems, any particle or scatterer in the gas follows a non-interacting Newtonian evolution during time intervals separating the collisions it experiences. In each collision an energy conserving interaction occurs between the particle and scatterer involved. Such systems, in which the scatterers themselves possess an internal dynamics, have been referred to as dynamical Lorentz gases [13], since they provide a natural generalization of the two-component gases originally introduced by Lorentz [START_REF] Lorentz | KNAW, Proceedings[END_REF], and subsequently studied by many others [START_REF] Dettmann | [END_REF], that feature inert (i.e., non-dynamical) scatterers, and which for over a century have been an essential tool for understanding diffusion, related equilibrium and non-equilibrium statistical properties, and the so-called Boltzmann-Grad limit.

We focus here on thermalization in a specific twodimensional system of this type, introduced in [8] as the rotating Lorentz gas model (RLG), in which the scatterers are rotating disks. The full many particle RLG exhibits realistic equilibrium and non-equilibrium behaviour [9], and has been used to study thermal rectification [10] and thermoelectricity [11]. Its time-reversible dynamics preserves phase space volume.

To understand equilibration of such RLGs as a whole, we investigate here the effect of repeated interactions on the phase space distribution of a single member of each of the two components making up the RLG. In the context of non-equilibrium temperatures, one might think of an ensemble of such single members as a "thermometer" locally probing the phase space distribution of the other component.

Thus, we first consider a single hollow disk that rotates freely about its center immersed in a gas of (and subjected to repeated impacts by) non-interacting point particles (Fig. 1a), with the particles of the gas reservoir being "probed" drawn independently from a stationary, but not necessarily thermal distribution. This obviously ignores recollisions, which would modify the momentum distribution of gas particles in the full RLG as they repeatedly encounter different scatterers. In the second part of our analysis, therefore, we follow a single particle as it passes through (i.e., probes) a spatially-fixed gas of rotating disks (Fig. 1b), where now it is the disks of the scatterer reservoir that are drawn independently from an arbitrary stationary distribution.

As we prove, in either of these situations, the probe dynamics reduce to a Markov chain that leads in the limit of small average energy exchange to an approach of the corresponding probe energy distribution to a Boltzmann state with a well-defined effective temperature.

In either case, each step of the Markov chain involves a collision (Fig. 2) between a particle of mass m = 1 and initial momentum p impinging with impact parameter b upon a hollow disk of unit radius and mass µ = M/m, rotating about its fixed center with initial angular velocity ω. During the collision, particle and disk obtain new momentum and angular velocity [8]

p ′ = p -2(p • u)u - 2µ 1 + µ (p t -ω) u ⊥ , (1)
ω ′ = ω + 2 1 + µ (p t -ω). (2)
Here, µ represents both the mass of the disk and its mo- ment of inertia, the unit vector u links the disk center to the collision point, u ⊥ is obtained by rotating u through π/2, and p t = p • u ⊥ = b p . Under these rules, particle and disk exchange energy and angular momentum, conserving both.

One rotating disk in a gas of particles -Consider a single disk (the "thermometer") with initial angular velocity ω 0 , subject to impacts by point particles (Fig. 1a) with i.i.d. momenta p n drawn from a rotationally invariant distribution ρ eff (|p|), with i.i.d. impact parameters b n uniform in [-1, 1], at a Poissonian sequence t n of impact times. Note, the marginal distribution ρ eff (|p|) ∝ ρ 0 (|p|)|p| of particles striking the disk, related to the "effusive" flux j(p) = ρ 0 (|p|)p of particles on the disk, is different from the bulk gas distribution ρ 0 (|p|).

Under these conditions p t = 0. After n collisions (2) gives

ω n = γ n ω 0 + η n-1 k=0 γ k p t,n-1-k , (3)
where η = 2/(1 + µ) and γ = 1 -η ∈ (-1, 1). As n → ∞, the first term in (3) vanishes, leaving a convergent sum

ω = lim n→+∞ η n-1 k=0 γ k p t,n-1-k (4)
of geometrically weighted i.i.d. random variables p t , ensuring that the angular velocity distribution approaches a stationary limit ρ(ω).

When the particle gas is in thermal equilibrium, the p t in (4) are Gaussian. Hence, independent of the value µ, the distribution ρ(ω) is thermal and has the same temperature as the bath (see below). Thus, in contact with thermal particles, a probe disk approaches a common thermal equilibrium with the particle bath.

Less obvious and more interesting is when the particle reservoir is not thermal [12]. Here, we characterize the limiting distribution by its low order moments. Assuming ω 0 = 0, then ω n = 0, and

ω 2 n = γ 2n ω 2 0 + η 2 p 2 t 1 -γ 2n 1 -γ 2 . ( 5
)
As n → ∞, the average disk energy then approaches On the right, the average b 2 = 1/3 has been performed, and the remaining average is over ρ eff (|p|). To relate this to the particle energy ε p = p 2 0 /2 [where . . . 0 indicates averages over ρ 0 (|p|)], we introduce an energy partitioning ratio ε = ε d /ε p , which equals 1/2 in equilibrium, satisfying the equipartition theorem. For a disk immersed in a non-equilibrium particle gas, however, ε = p 2 eff / 3p 2 0 approaches a non-universal value that depends on the particle distribution, but is independent of the mass of the disk (see Fig. 3). For a thermal particle bath, ε = 1/2 and (6) reduces to the equilibrium result µ ω 2 /2 = k B T /2.

ε d = µ 2 ω 2 = 1 2 p 2 t = 1 6 p 2 eff . (6)
To study the shape of ρ(ω) we compute the asymptotic value of the excess disk kurtosis

κ ω = β 2 -3 = η 1 + γ 1 + γ 2 p 4 t p 2 t 2 -3 , (7)
in which β 2 = ω 4 / ω 2 2 , and which vanishes when the limiting disk distribution is thermal. This is the case for a thermal particle bath, when the excess particle kurtosis κ p = p 4 t / p 2 t 2 -3 in (7) vanishes. But ρ(ω) also reduces to a thermal distribution when the particle reservoir is not thermal, whenever the coupling is sufficiently weak. This occurs when the disk is very heavy (µ ≫ 1, η → 0, γ → 1), and very light (µ ≪ 1, η → 2, γ → -1). In these two limits (7) vanishes and ρ(ω) becomes Gaussian [20]. In weak coupling, this thermalization is universal, i.e., independent of the form of ρ eff (|p|).

Thus, when placed weakly in contact with a stationary, non-equilibrium particle gas, the probe disk "thermalizes to" (or "measures") an apparent temperature such that k B T /2 = ε p 2 0 /2, where ε = p 2 eff / 3p 2 0 depends on the particle distribution and is generally not equal to 1/2. This thermalization for large and small µ is clearly seen in Fig. 3, which displays a numerical computation of the values of ε and κ ω , for 10 5 disks, each subjected to 10 4 repeated collisions for various particle distributions.

From Fig. 3 and Eq. ( 7) it is clear that for intermediate and strong coupling, with particle reservoirs for which κ p = 0, the disk does not approach a thermal state. Instead, repeated strong interactions drive it to a nonthermal state with κ ω = 0.

One particle in a gas of disks -We now consider a single particle that collides with a sequence of rotating disks (Fig. 1b) whose angular velocities are i.i.d. variables drawn from a stationary distribution ρ(|ω|), again ignoring recollisions. Denote by p s the particle's momentum before collision s, when it impinges with impact parameter b s on a disk with angular velocity ω s . According to [START_REF] Khinchin | Mathematical Foundations of Statistical Mechanics[END_REF], after this collision

p s+1 = p s -2(p s • u s )u s - 2µ 1 + µ (p s,t -ω s ) u s,⊥ . ( 8
)
Collisions occur at times t s+1 = t s + ℓ/p s , where ℓ is the collision mean free path. We assume p 0 = 0. The Markov chain (8) for p is less tractable than for a single disk in a gas of particles. In weak coupling, however, analysis of the Markov chain to determine the limiting distribution is straightforward.

Figures 4 and5 display simulation results in which 10 5 particles each undergo a sequence of collisions with rotating disks initialized as described. In Fig. 4, values of ε and κ p for the limiting particle distribution, recorded (•) to its limiting value 3 (dashed line), and kinetic energy ratio ε = ε d /εp ( ) to its limiting value 1/2 (solid line).

at long fixed simulation time, are plotted for different disk distributions ρ(|ω|).

For thermal disks, the particle equilibrates to a thermal state with the same temperature. Moreover, when µ is very large or very small, independent of ρ(ω), the particle is also driven to a thermal state (see, e.g., Fig. 5), with a vanishing κ p . Thus, we again observe "thermalization" of a small system in weak contact with a non-thermal reservoir. However, unlike what happens to a disk in a gas of particles, the limiting particle temperature always obeys equipartition, since ε → 1/2 in this regime. For intermediate µ, provided κ ω = 0, the particle does not generally thermalize.

To analytically demonstrate the thermalization observed at large and small µ, we express using (8) the change in particle energy

E s+1 -E s = - 4µ (1 + µ) 2 E s b 2 s + 2µ (1 + µ) 2 µω 2 s + 2 √ µ 1 -µ (1 + µ) 2 2E s √ µω s b s (9)
during collision s, in terms of its energy immediately before. The coefficients multiplying the dynamical variables in (9) set the scale for the energy change in any collision, and are small for µ ≪ 1 and µ ≫ 1. In the latter case the equation for large µ follows from that for small µ by replacing µ with µ -1 .

It suffices to study this weak-coupling, small-step limit for µ ≪ 1. In this limit, (9) becomes

∆E s = -4µE s b 2 s + 4µε d k 2 s + 2 4µ ε d E s k s b s , (10)
where

ε d = 1 2 µ ω 2 s , k s = 2µ ε d ω s , k 2 s = 1. (11)
Averages here are over ρ(|ω|). The first two terms on the right hand side of (10) constitute a source of "dynamical friction" [21,22] that counterbalances the stochastic acceleration caused by the fluctuating last term [13,16,21,22]. This competition thus leads to a kind of fluctuation-dissipation like mechanism that naturally emerges from the deterministic dynamics. Inspection of ( 10) suggests ξ = 4µ as a small parameter. Introducing (scaled) continuous collision number σ = ξs, Eq. ( 10) can be described by the stochastic differential equation

dE(σ) = -α(E)dσ + λ(E)dw(σ) in which α(E) = 1 3 E(σ) -ε d , λ(E) = 2 √ 3 ε d E(σ) .
Here dw is white noise, dw = 0, and dw(σ)dw(σ ′ ) = δ(σ -σ ′ ). Clearly, the results will depend on only the first two moments of ρ(|ω|).

From the corresponding Fokker-Planck equation, one finds the stationary distribution

f ∞ (E) = Z-1 λ -2 (E) exp[ -2α(E)λ -2 (E)dE ] = 1 Z ′ E 1/2 exp(- E 2ε d )
for E(σ) at large σ. Physically, f (E, σ)dE is the fraction of particles after collision σ with energy between E and E + dE. Such particles stay in that state for time ∆t = ℓ/ √ 2E. Thus, at long times, the distribution of particles with energy E becomes thermal

f (E) = 1 Z exp(- E 2ε d ) (12)
as in a 2D non-interacting particle gas with apparent temperature given by k B T = 2ε d . Thus, rotating disks not in equilibrium drive the particle to a thermal state compatible with equipartition.

On the approach to equilibrium -In weak coupling, as we have shown, both components of the full RLG are driven to a Gaussian thermal state. From (8), and (10), it is clear that equilibration is not a simple consequence of the central limit theorem. Indeed, the momentum or energy increments that particles experience are not i.i.d. random variables; they are steps in an associated Markov chain, with the size and variance of each being a function of the dynamical variables at each step. Thus, the limiting exponential distribution for the particle energy in this case can be understood as arising from a competition between impulsive fluctuating forces of zero mean, which tend to heat the particle, and a dynamical friction that leads higher energy particles to lose energy to the disks. Thus, a fluctuation-dissipation mechanism emerges naturally from the deterministic dynamics associated with the interaction between the particles and the rotating disks.

The fact that in the present case both components of the full RLG equilibrate in the presence of a non-thermal bath allows us to identify the mechanism of approach to thermal equilibrium of the RLG as a whole.

At low particle density, e.g., each particle will scatter off many disks before any disk is likely to have interacted with more than a few particles. In weak coupling, therefore, the particle gas will equilibrate well before the gas of disks, to a thermal distribution in which equipartition of the particle and disk energies is obtained. Each disk will then be in contact with a thermal distribution of particles, and will only need therefore to undergo a return to the appropriate limiting thermal distribution. Beyond weak coupling, the situation becomes more complicated [12]. Moreover, extensive numerical simulations on large but finite systems with many particles and disks, to be reported elsewhere [12], confirm this basic picture of a two-stage equilibration process of the whole system in the weak coupling regime.

In summary, we have identified a thermalization mechanism in a rotating disk, weak-coupling version of a twocomponent dynamical Lorentz gas. We expect the basic underlying dynamical friction mechanism to be effective more generally in systems in which the individual components undergo repeated scattering events (See, for example Ref. 13)

FIG. 1 .

 1 FIG. 1. (a) A single rotating disk in a gas of point particles, and (b) a single point particle passing through a spatially fixed array of rotating disks.

FIG. 2 .

 2 FIG.2. Collision geometry for a point particle of unit mass and initial momentum p impinging on a freely rotating disk of mass µ and initial angular velocity ω.

FIG. 3 .

 3 FIG. 3. Limiting value of κω of a disk of mass µ subjected to repeated impacts by particles drawn from thermal (×), uniform ( ), and microcanonical (+) distributions, and a superposition of two thermal distributions (•) with temperatures in ratio T1/T2 = 5. Curves are theoretical predictions, symbols are simulation results. Inset: Limiting value of kinetic energy ratio ε = ε d /εp vs µ. Theoretical predictions are ε = 1/2, ε = 2/5, and ε = 1/3.

FIG. 4 .

 4 FIG.4. Limiting values of the kinetic energy ratio ε = ε d /εp and excess kurtosis κp for a particle colliding with rotating disks of mass µ and angular velocities drawn from thermal (×), uniform ( ), and microcanonical (+) distributions, and a superposition of two thermal distributions (•) with T1/T2 = 5.

1 FIG. 5 .

 15 FIG. 5. Distribution function for velocity component px for a particle undergoing repeated collisions with rotating disks of mass µ = 0.005 in a microcanonical ensemble with |ω| = √ 20, recorded (from broad to narrow) at times t = 60, 300, 600, 2000. Dashed curve corresponds to the limiting thermal distribution. Inset: convergence of the particle kurtosis β (p) 2

ACKNOWLEDGMENTS

Part of this work was performed while C.M.-M. and P.E.P. visited the Université Lille 1 and the Labex CEMPI (ANR-11-LABX-0007-01). They thank those institutions for their hospitality. C.M.-M. acknowledges partial financial support from the Spanish MICINN grants MTM2012-39101-C02-01 and MTM2015-63914-P, and from ONRG Grant N62909-15-1-C076.


Mouna Gharbi

Emilie Chouzenoux

Jean-Christophe Pesquet

An Unrolled Half-Quadratic Approach for Sparse Signal Recovery in Spectroscopy

Keywords: Unrolled architecture, Signal reconstruction, Half-Quadratic, Spectroscopy

This paper addresses the problem of sparse signal reconstruction frequently arising in analytical chemistry. Spectroscopic restoration is an important task that was tackled using two distinct approaches. First, model-based iterative methods were deployed. These methods rely on a strong theoretical background. However they face practical obstacles in terms of hyperparameter tuning. Conversely, approaches based on deep neural networks are fast and easy to implement. Nonetheless, they are often described as black-box tools, and their stability/robustness are ongoing research topics. In this paper, we propose a deep neural network based on unrolling a Majorization-Minimization (MM) Half-Quadratic (HQ) algorithm. This allows us to build interpretable layers mirroring iterations, making it possible to learn automatically data-driven hyperparameters such as regularization and stepsizes. Furthermore, we propose a dictionary of custom activation functions derived from potentials used in the original variational model. This interpretation of activations can be useful for analyzing the stability of neural networks. The efficiency of our method in comparison to iterative and learning-based methods is showcased through various experiments conducted on realistic mass spectrometry (MS) databases with various blur kernels and noise levels. All experiments show an improved reconstruction quality of chemical signals, while maintaining a low execution time at the test phase. All the results are reproducible, through the code shared on a repository available at https://github.com/GHARBIMouna/Unrolled-Half-Quadratic.

Introduction

Spectroscopy gathers data acquisition/processing pipelines aiming to decipher the chemical properties of a material or a compound. Let us cite, for instance, mass spectrometry (MS) [START_REF] Heaney | Mass spectrometry in medicine: a technology for the future?[END_REF] and nuclear magnetic resonance (NMR) [START_REF] Ernst | Principles of Nuclear Magnetic Resonance in One and Two Dimensions[END_REF]. Despite constant improvements in physical acquisition devices, the acquisition process often remains noisy and prone to degradation. More formally, let x ∈ R n be an original chemical signal we aim to estimate, and let y ∈ R m be an acquired measurement. In many spectroscopy applications (e.g., RMN DOSY [START_REF] Cherni | PALMA, an improved algorithm for DOSY signal processing[END_REF], RMN relaxometry [START_REF] Chouzenoux | Optimization of a maximum entropy criterion for 2D Nuclear Magnetic Resonance reconstruction[END_REF], MS [START_REF] Cherni | SPOQ ℓ p -over-ℓ q regularization for sparse signal recovery applied to mass spectrometry[END_REF]), both are linked by the following model

y = Hx + e, (1)
where H ∈ R m×n represents a measurement degradation process modeling for example a blur kernel [START_REF] Liu | A map-based algorithm for spectroscopic semi-blind deconvolution[END_REF] or a discrete Laplace transform [START_REF] Lin | High-resolution reconstruction for multidimensional Laplace NMR[END_REF], and e ∈ R m is an additive noise corrupting the data.

To draw relevant chemical interpretations, measurements y need to be processed to recover an estimate of the original data x, which amounts to inverting Model [START_REF] Heaney | Mass spectrometry in medicine: a technology for the future?[END_REF].

In a real-life application context, this inversion process may encounter several difficulties. On the one hand, the ill-posedness of (1), linked to the poor conditioning of H, prevents the use of basic inversion or filtering approaches. On the other hand, the necessity to solve (1) multiple times, for each new acquisition, induces hard constraints on the processing time. To cope with the aforementioned challenges, two main strategies in literature have been considered. A first set of methods relies on optimization, where a criterion combining a data fidelity term and a penalization term is defined and minimized to estimate the signal of interest. This approach was employed in many spectroscopy restoration works such as [START_REF] Cherni | SPOQ ℓ p -over-ℓ q regularization for sparse signal recovery applied to mass spectrometry[END_REF][START_REF] Mohammad-Djafari | Regularization, maximum entropy and probabilistic methods in mass spectrometry data processing problems[END_REF] for MS data processing, [START_REF] Chouzenoux | Optimization of a maximum entropy criterion for 2D Nuclear Magnetic Resonance reconstruction[END_REF][START_REF] Chouzenoux | Efficient maximum entropy reconstruction of Nuclear Magnetic Resonance T 1 -T 2 spectra[END_REF][START_REF] Song | T 1 -T 2 correlation spectra obtained using a fast twodimensional Laplace inversion[END_REF][START_REF] Rondeau-Mouro | Two dimensional IR-FID-CPMG acquisition and adaptation of a maximum entropy reconstruction[END_REF][START_REF] Venkataramanan | Solving Fredholm integrals of the first kind with tensor product structure in 2 and 2.5 dimensions[END_REF] in the context of NMR relaxometry, and in DOSY NMR [START_REF] Cherni | PALMA, an improved algorithm for DOSY signal processing[END_REF].

Recently, a second set of approaches has emerged, which relies on supervised machine learning techniques applied on real or synthetic chemical signal databases to learn mappings (typically, deep neural networks) capable of producing an estimate of the sought signal from the degraded data. For example, dense convolutional neural networks (CNNs) were used to remove artifacts from high-dimensional NMR signals in [START_REF] Qu | Accelerated Nuclear Magnetic Resonance spectroscopy with deep learning[END_REF] and, in [START_REF] Kim | Compressive sensing spectroscopy using a residual convolutional neural network[END_REF], a residual convolutional neural network, was designed for compressive sensing spectroscopy. The work [START_REF] Debus | Deep learning in analytical chemistry[END_REF] reviews various deep networks such as LSTMs, CNNs and DNNs, used in the context of analytical chemistry. Let us also mention the survey paper [START_REF] Chen | Review and prospect: deep learning in Nuclear Magnetic Resonance spectroscopy[END_REF], presenting an overview of various deep learning methods for NMR spectra processing, including the recent works from [START_REF] Karunanithy | FID-Net: A versatile deep neural network architecture for NMR spectral reconstruction and virtual decoupling[END_REF][START_REF] Hansen | Using deep neural networks to reconstruct non-uniformly sampled NMR spectra[END_REF][START_REF] Wu | Improvement in signal-to-noise ratio of liquid-state NMR spectroscopy via a deep neural network DN-Unet[END_REF]. Other spectroscopy modalities, such as Raman and FTIR, have also been considered, from the side of deep learning, for instance in [START_REF] Horgan | High-Throughput molecular imaging via deep-learning-enabled Raman spectroscopy[END_REF] with residual Unets and [START_REF] Brandt | Deep learning for reconstructing low-quality FTIR and Raman spectra a case study in microplastic analyses[END_REF] via autoencoders.

Supervised learning-based methods can lead to high-quality results, as soon as a sufficiently large labeled database is available. Furthermore, they have proven to be fast since they benefit from the advances of GPU-based deep learning programming frameworks. However, the interpretability and explainability of these methods remain open issues (see [START_REF] Haar | An analysis of explainability methods for convolutional neural networks[END_REF][START_REF] Neacsu | Design of robust complex-valued feed-forward neural networks[END_REF][START_REF] Gupta | Multivariate lipschitz analysis of the stability of neural networks[END_REF] for recent works on these aspects). They might suffer from instability issues [START_REF] Antun | On instabilities of deep learning in image reconstruction and the potential costs of AI[END_REF][START_REF] Combettes | Lipschitz certificates for layered network structures driven by averaged activation operators[END_REF], and guaranteeing their robustness remains a challenge [START_REF] Genzel | Solving inverse problems with deep neural networks -robustness included?[END_REF]. In contrast, optimization-based methods offer the required theoretical guarantees, but tuning their hyperparameters might be a tedious task, and their execution time for CPU based implementations may be too high in large scale applicative contexts. For instance, the Matlab implementation1 of the MS signal restoration method from [START_REF] Cherni | SPOQ ℓ p -over-ℓ q regularization for sparse signal recovery applied to mass spectrometry[END_REF] requires more than 1 minute to restore an MS spectra on a grid of n = 10 3 mass values.

Lately, the two previous approaches have been combined, yielding an appealing class of techniques based on 'unrolling' (or 'unfolding'). The core idea is to re-interpret existing iterative optimization algorithms as neural network architectures [START_REF] Hershey | Deep unfolding: Modelbased inspiration of novel deep architectures[END_REF], where layers mirror iterations. A recent survey of unrolling techniques can be found in [START_REF] Monga | Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing[END_REF]. This new paradigm presents several advantages. First, new architectures can be proposed, which enjoy an interpretable network structure. Second, prior and/or physical knowledge about the problem at hand can be naturally embedded into the architecture. Third, a reduced execution time is obtained during the test phase, thanks to deep learning GPU-based frameworks. Lastly, unrolled architectures are more flexible than optimization algorithms, and in particular, they offer the possibility to learn some hyperparameters and/or operators of the original algorithm. For instance, stepsizes and regularization weights are learnt in [START_REF] Bertocchi | Deep Unfolding of a Proximal Interior Point Method for Image Restoration[END_REF], operators and penalty function parameters in [START_REF] Ito | Trainable ISTA for sparse signal recovery[END_REF][START_REF] Savanier | Deep Unfolding of the DBFB Algorithm With Application to ROI CT Imaging With Limited Angular Density[END_REF], and shrinkage functions and biases are learnt in [START_REF] Sun | Color image denoising via discriminatively learned iterative shrinkage[END_REF].

In various applications [START_REF] Monga | Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing[END_REF], unrolling techniques have been shown to outperform (sometimes, by far) existing methods, with a noticeable gain in explainability and robustness with respect to plain deep learning techniques, while leading to similar (fast) computing times [START_REF] Le | The faster proximal algorithm, the better unfolded deep learning architecture ? the study case of image denoising[END_REF]. This paper investigates unrolled algorithms dedicated to the resolution of spectroscopy inverse problems. We explore specifically the unrolling of a popular class of optimization algorithms called Majorization-Minimization (MM) approaches [START_REF] Sun | Majorization-Minimization algorithms in signal processing, communications, and machine learning[END_REF]. Namely, we focus on the MM half-quadratic (HQ) approach [START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF] which minimizes sequentially quadratic majorizing surrogates of the objective function.

HQ was initially developed for image restoration applications [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF], then deployed for chemical signal recovery in [START_REF] Mazet | Baseline spectrum estimation using halfquadratic minimization[END_REF][START_REF] Liu | Goldindec: A novel algorithm for Raman spectrum baseline correction[END_REF]. Processing large-scale data with MM HQ has been explored in [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF] by relying on subspace acceleration strategy [START_REF] Chouzenoux | A majorize-minimize strategy for subspace optimization applied to image restoration[END_REF][START_REF] Chouzenoux | A Majorize-Minimize subspace approach for ℓ 2 -ℓ 0 image regularization[END_REF], or parallel CPU-based implementations [START_REF] Cadoni | A block parallel Majorize-Minimize memory gradient algorithm[END_REF][START_REF] Chalvidal | Block distributed 3MG algorithm and its application to 3D image restoration[END_REF][START_REF] Hong | A unified algorithmic framework for block-structured optimization involving big data: With applications in machine learning and signal processing[END_REF][START_REF] Breloy | Block Majorization-Minimization algorithms for low-rank clutter subspace estimation[END_REF].

Our contributions in this work are as follows:

• The proposition of an unrolled MM HQ algorithm, for which an efficient GPU-based implementation is available.

• The construction of new activation functions, with established connection with choices of penalty functions.

• The design of a supervised architecture allowing us to learn the involved regularization and stepsize parameters.

• The validation of the proposed approach in the restoration of signals arising from realistic databases of MS signals.

Link with existing works: First, as we mentioned, deep learning strategies or optimization-based ones, have already been considered for tackling spectroscopy inverse problems. Up to our knowledge, deep unrolling has been only scarcely explored in the field of spectroscopy, and its deployment for MS signal processing remains unprecedented. We can mention [START_REF] Bobin | Learning to unmix from Poisson measurements with application to γ-spectroscopy[END_REF][START_REF] Ahmadi | Photothermal-SR-Net: A customized deep unfolding neural network for photothermal super resolution imaging[END_REF][START_REF] Wang | A sparse model-inspired deep thresholding network for exponential signal reconstruction-application in fast biological spectroscopy[END_REF] proposing deep unfolded neural networks to solve inverse problems in γ-ray spectroscopy, photothermal radiometry, and non-uniformly sampled NMR, respectively. In the broader context of sparse 1D signal recovery, we can cite [START_REF] Yang | A robust deep unfolded network for sparse signal recovery from noisy binary measurements[END_REF][START_REF] Li | Jointly sparse signal recovery via deep autoencoder and parallel coordinate descent unrolling[END_REF][START_REF] Huijben | Learning sampling and model-based signal recovery for compressed sensing MRI[END_REF]. Although the deep unrolling paradigm in these papers is common to our study, the methods highly differ in the physical model they are considering, in the signal processing problem that is tackled, and in the iterative algorithm that is unrolled. Second, unrolling of MM methods is also not much explored in the literature. We can only mention [START_REF] Li | An algorithm unrolling approach to deep image deblurring[END_REF] which proposes to unroll an HQ method for an application to image deblurring. Third, our analysis connecting penalty functions and activation layers generalizes the one presented in [START_REF] Alt | Translating diffusion, wavelets, and regularisation into residual networks[END_REF].

The paper is organized as follows, Section 2 introduces the optimization problem, and the MM HQ algorithm to solve it. Section 3 details our contributions, mainly the unrolling of HQ algorithm, the proposition of new activation functions, and the design of architectures for hyperparameter learning. Section 4 presents our experimental results and discussions. Lastly, Section 5 concludes the work.

Background on MM-based restoration methods

Problem formulation

A classical strategy for solving inverse problem (1) and producing an estimate x ∈ R n of x from observed data y, consists in solving a penalized least-squares minimization problem [START_REF] Willoughby | Solutions of ill-posed problems (A. N. tikhonov and V. Y. arsenin[END_REF] of the form minimize

x∈R n F (x) = 1 2 ∥Hx -y∥ 2 + λ 1 Ψ 1 (x) + λ 2 Ψ 2 (x) . (2)
Hereabove, the regularization functions Ψ i : R n → R, i ∈ {1, 2}, aim at ensuring the stability of the solution x, by enforcing a priori assumptions related to sparsity. A generic expression for such sparsity-promoting regularization functions is given, for i ∈ {1, 2}, by

(∀x ∈ R n ) Ψ i (x) = n s=1 ψ i ([x] s ), (3)
where ψ i : R → [0, +∞) is some positive-valued potentials, and [•] s with s ∈ {1, . . . , n} denotes the s-th component of a vector in R n . λ 1 > 0 and λ 2 > 0 are regularization hyperparameters balancing the weight between both regularization terms (Ψ i ) i∈{1,2} , and the least-squares fidelity in the objective function F . Let us provide the rationale for adopting such a hybrid two-part penalty term. Sparsity-based priors are widely used in the signal processing literature so as to promote few nonzero entries in the sought signal (or in its representation in a transformed domain [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF]). Two main families of such priors can be distinguished, namely (i) convex ones, (ii) non-convex ones. Convex sparse penalties include, as a notorious example, the ℓ 1 norm (related to the LASSO approach in machine learning [START_REF] Park | The Bayesian Lasso[END_REF]). An extension of this penalty is the elastic net regularizer [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]. Smooth approximations for it have also been very popular, such as the Huber penalty or the Fair potential [START_REF] Guitton | Robust inversion of seismic data using the Huber norm[END_REF]. Convex penalties are convenient since they lead to more tractable optimization problems associated with efficient solvers with sounded convergence guarantees. However, they may lead to a bias in the estimated non-zero coefficients, which might be detrimental to the final quality of the results or to their quantitative interpretation. There is a plethora of non-convex priors, starting with the popular ℓ 0 quasi-norm up to its various approximations through norm ratios [START_REF] Repetti | Euclid in a Taxicab: Sparse blind deconvolution with smoothed ℓ 1 /ℓ 2 regularization[END_REF], piecewise defined penalties [START_REF] Soubies | A continuous exact ℓ 0 penalty (CEL0) for least-squares regularized problem[END_REF], and smoothed ℓ 0 potentials [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF]. Such choices are closer, in essence, to the true goal of sparse restoration, as they directly aim at minimizing the number of significant components in the sought signal, without altering (too much) their amplitudes. Thus, they can lead to much better restoration quality, especially for very sparse signals, and avoid the biasing effect of the convex penalties. However, non-convexity has a significant cost in terms of difficulty of the optimization procedure, making the algorithms sometimes unstable and prone to getting stuck in spurious local minima. In this work, we propose to get the best of both worlds by mixing the two strategies, setting ψ 1 as a convex sparsity-enhancing potential and ψ 2 as a non-convex one. Furthermore, we opt for smooth (i.e., differentiable) regularizers, so as to preserve an efficient optimization procedure. It can be understood that, in such a context, the choice of proper regularization parameters (λ 1 , λ 2 ) will play a critical role. The unrolling approach will serve to decipher, in an automatic, supervised manner, the setting of them, as we will explain in Section 3.

Quadratic majorization-minimization algorithm

The estimate x of x is obtained by minimizing the penalized least squares formulation (2) by assuming that this optimization problem admits a solution. This is, in particular, guaranteed if the cost function in (2) is coercive, which is satisfied, for standard choices of function Ψ 1 , provided that λ 1 > 0. Given our choices for ψ 1 and ψ 2 , a popular class of minimization schemes rely on the key concept of quadratic majorization-minimization we introduce hereafter.

Before presenting our algorithm, let us be more specific about the assumption on our potential functions ψ 1 and ψ 2 . Assumption 1. For every i ∈ {1, 2}, ψ i : R → R is (i) a differentiable, even function, (ii) increasing on [0, +∞),

(iii) such that ψ i ( √ •) is concave on [0, +∞).
In addition, ψ 1 is convex.

Assumption 1 is satisfied by a wide class of convex, and non-convex potentials, listed in the second column of Table 1 (see also [START_REF] Chouzenoux | A stochastic Majorize-Minimize subspace algorithm for online penalized least squares estimation[END_REF]).

In the context of Problem (2) involving penalties satisfying Assumption 1, an efficient strategy consists of adopting a half-quadratic (HQ) optimizer [START_REF] Nikolova | Analysis of half-quadratic methods for signal and image recovery[END_REF][START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF]. HQ belongs to the class of Majorization-Minimization (MM) algorithms, which is a prominent family of methods in the field of optimization [START_REF] Zhang | Surrogate maximization/minimization algorithms and extensions[END_REF][START_REF] Hunter | A tutorial on MM algorithms[END_REF], with successful applications in the context of signal/image processing, telecommunications, and machine learning [START_REF]Incremental majorization-minimization optimization with application to large-scale machine learning[END_REF][START_REF] Sun | Majorization-Minimization algorithms in signal processing, communications, and machine learning[END_REF]. MM techniques substitute the often complex problem of minimizing F by a sequence of more tractable problems, based upon the minimization of sequences of majorant functions (also called surrogate functions or auxiliary functions), satisfying the following conditions: at every iteration k ∈ N of the algorithm,

(∀x ∈ R n ) F (x) ≤ Q(x, x k ), F (x k ) = Q(x k , x k ), (4)
where x k ∈ R n is the current iterate. Under Assumption 1, a quadratic majorant function for F at x k ∈ R n is given by

(∀x ∈ R n ) Q(x, x k ) = F (x k ) + ∇F (x k ) ⊤ (x -x k ) + 1 2 (x -x k ) ⊤ A(x k )(x -x k ), (5)
where A(x k ) is symmetric definite positive matrix in R n×n , defined as

(∀x ∈ R n ) A(x) = H ⊤ H + λ 1 Diag{(ω 1 ([x] s )) 1≤s≤n } + λ 2 Diag{(ω 2 ([x] s )) 1≤s≤n },
where, for every i ∈ {1, 2},

(∀u ∈ R \ {0}) ω i (u) = ϱ i (u) u , ω i (0) = lim u→0 ϱ i (u) u , (6)
with ϱ i : R → R the derivative of the scalar potential function ψ i . Note that the above limit is well-defined, thanks to Assumption 1 [START_REF] Nikolova | Analysis of half-quadratic methods for signal and image recovery[END_REF]. The expression of the derivatives for a bunch of potential functions satisfying Assumption 1 is provided in the third column of Table 1. Additional properties for the potential functions, also given in the table, will be discussed in Subsection 3.2. Solving Problem (2) using the MM algorithm associated with quadratic majorant function [START_REF] Cherni | SPOQ ℓ p -over-ℓ q regularization for sparse signal recovery applied to mass spectrometry[END_REF] yields

(∀k ∈ N) x k+1 = x k -γ k A(x k ) -1 ∇F (x k ), (7)
where x 0 ∈ R n and, for every

k ∈ N, A(x k ) -1 is the inverse of A(x k ), ∇F (x k ) = H ⊤ (Hx k -y) + λ 1 Ω 1,k x k + λ 2 Ω 2,k x k ,
with, for i ∈ {1, 2},

Ω i,k = Diag{(ω i ([x k ] s )) 1≤s≤n }. ( 8
)
{γ k } k∈N ⊂ (0, 2
) is a sequence of stepsizes. Algorithm [START_REF] Lin | High-resolution reconstruction for multidimensional Laplace NMR[END_REF] can be viewed as a preconditioned gradient method, or a quasi-Newton method, where the scaling matrix multiplying the gradient descent term is built following the MM principle, in such a way that (F (x k )) k∈N decreases monotonically. It can also be interpreted as an alternating minimization method acting in an extended space, where the function to minimize is quadratic with respect to some auxiliary vector variable. This interpretation explains the name given to HQ algorithm, and is actually a key element for proving its convergence [START_REF] Allain | On global and local convergence of half-quadratic algorithms[END_REF]. The performance of HQ method has been illustrated in the context of signal restoration [START_REF] Nikolova | Analysis of half-quadratic methods for signal and image recovery[END_REF].

The matrix inversion required in [START_REF] Lin | High-resolution reconstruction for multidimensional Laplace NMR[END_REF] can be detrimental to a stable and fast practical implementation of this method. Moreover, setting the stepsize sequence, and the penalization hyperparameters (λ 1 , λ 2 ) in the hybrid regularizer might be tedious. To this end, we propose in the next section an unrolling procedure that helps reducing computational time by creating an architecture with a given number of layers, benefiting from the advantages of implementations on GPUs, and with automatized hyperparameter tuning.

Proposed unrolled network

Unrolled HQ architecture

The first important step for algorithm unrolling consists in translating the algorithm of interest into a multi-layer network architecture. Let us consider a given number K > 0 of iterations for HQ algorithm. We will show in this section that each iteration k ∈ {0, . . . , K -1} of it can be mirrored by a layer L k of a feedforward network with a residual connection [START_REF] Combettes | Deep neural network structures solving variational inequalities[END_REF]. Let us rewrite the update rule of HQ algorithm [START_REF] Lin | High-resolution reconstruction for multidimensional Laplace NMR[END_REF] as

(∀k ∈ N) x k+1 = x k -W k ∇F (x k ), (9)
where

x 0 ∈ R n , W k = γ k A(x k ) -1 , (10)
and the gradient expression is given in [START_REF] Mohammad-Djafari | Regularization, maximum entropy and probabilistic methods in mass spectrometry data processing problems[END_REF]. We now introduce linear and non linear operators that will play the role of linear layers and non-linear functions in the unrolled architecture. Let

V 0 = H ⊤ H, V 1 = I n , V 2 = I n , V 3 = [I n λ 1 I n λ 2 I n ], (11)
b 0 = -H ⊤ y, b 1 = 0, b 2 = 0, b 3 = 0, (12)
where I n is the identity matrix of size n × n. Moreover, let R 0 to R 3 be defined as, for every

x ∈ R n , R 0 (x) = R 3 (x) = x, R 1 (x) = ϱ 1 ([x] s ) 1≤s≤n , R 2 (x) = ϱ 2 ([x] s ) 1≤s≤n . (13
) Note these activations operate componentwise. Activations R 0 and R 3 , both reducing to identity operators, are introduced for the sake of completeness of the presentation, so they play no role in the final output expression. Given the expressions ( 8) and ( 9), and the above definitions, we can rewrite the K first iterations of the HQ algorithm [START_REF] Lin | High-resolution reconstruction for multidimensional Laplace NMR[END_REF] as

(∀k ∈ {0, . . . , K-1}) x k+1 = x k -W k   R 3   V 3   R 0 (V 0 x k + b 0 ) R 1 (V 1 x k + b 1 ) R 2 (V 2 x k + b 2 )   + b 3     , (14)
with x 0 ∈ R n . Equation ( 14) has a remarkable structure, where we can identify fundamental bricks of a standard deep neural network architecture. For every k ∈ {0, . . . , K -1}, the input x k of layer L k first goes simultaneously through the 3-channel linear layer with weights V 0 , V 1 , V 2 , and biases b 0 , b 1 , b 2 , then through the activation operators R 0 , R 1 , R 2 . The three outputs of these activation blocks are concatenated and go through a linear layer with weight matrix V 3 , aggregating them. This layer has bias b 3 , and activation R 3 (here reducing to identity). The output x k+1 is finally obtained by computing the difference between the input x k and the output of R 3 composed with the linear operator W k . Such a last step is similar to a residual connection. Each layer of the network depends on the noisy observation. For every k ∈ {0, . . . , K -1} the k-th layer is denoted as L y k , so that, for every k ∈ {0, . . . , K -1},

x k+1 = L y k (x k ). ( 15
)
Figure 1 summarizes the structure of layer L y k for some k ∈ {0, . . . , K -1}. This equivalent formulation of the HQ algorithm presents several advantages. As we will demonstrate in Subsection 3.3, it is a first step towards the supervised learning of its hyperparameters such as the regularization pa-rameters and the stepsizes. The optimization of these parameters is based on a dedicated loss function directly related to the application at hand (e.g., mean squared error on the restored signal). This deep unrolling strategy leads to the design of a deep architecture with interpretable structure, which contrasts with traditional deep networks. Namely, the involved weight and bias operators have specific forms derived from the original inverse problem formulation. The nonlinear activation functions are not chosen in an ad hoc manner, but they are fully interpretable in an optimization sense. This will be further discussed in Subsection 3.2. Another specificity of our architecture with respect to conventional DL ones is that the weight operator W k depends on the k-th layer input. This weight operator thus exhibits some similarity to a multivariate attention mechanism. Moreover, relying on the unrolling paradigm offers certain mathematical guarantees of stability and robustness as it has been proven in [START_REF] Bertocchi | Deep Unfolding of a Proximal Interior Point Method for Image Restoration[END_REF][START_REF] Valle | Stability of unfolded forward-backward to perturbations in observed data[END_REF]. Finally, from a practical viewpoint, the neural network-based interpretation of the HQ method allows the deployment of GPU-friendly tools to implement the different bricks of the optimizer. This is beneficial in terms of computational efficiency, during the training process, when supervised unrolling is performed. This is also useful for computational time reduction at the inference stage, with or without trained parameters, as we will demonstrate in our numerical experiments.

Activation functions

The architecture displayed in Fig. 1 makes use of several activation functions whose properties are further analysed in this section. Since R 0 and R 3 reduce to identity functions, we focus our discussion on R 1 and R 2 , which are directly related to our choice for the penalization functions Ψ 1 and Ψ 2 in (2). As can be seen in ( 13), R 1 and R 2 are applied componentwise. Their expressions involve the scalar potentials ϱ 1 and ϱ 2 that are the derivatives of ψ 1 and ψ 2 , respectively. We list in the third column of Table 1 the derivative ϱ for typical choices of function ψ satisfying Assumption 1. Interestingly, one can recognize that some correspond to commonly used activation functions in neural networks, while others are less familiar. For example, the hyperbolic tangent is obtained for the Green potential, Eliott activation function corresponds to the Fair potential, and the Huber penalty leads to a shifted version of capped ReLU. As shown in [START_REF] Combettes | Lipschitz certificates for layered network structures driven by averaged activation operators[END_REF], proximal properties of activation functions play an important role in the mathematical understanding of neural networks, in particular in their stability analysis. We analyse further the properties of ϱ in that respect. Let Γ 0 (R) denote the class of lower-semicontinuous convex functions from R to R ∪ {+∞} which are proper (i.e., at least finite in one point). We recall that, if φ belongs to Γ 0 (R), its proximity operator prox φ : R → R is defined as

(∀x ∈ R) prox φ (x) = argmin t∈R φ(t) + 1 2 (t -x) 2 . ( 16
)
We are now ready to state the following property.

Proposition 1. Let ψ : R → R be a differentiable even function and let ϱ be its derivative. Assume that ϱ is 1-Lipschitz on R, i.e.

(∀(x, t) ∈ R 2 ) |ϱ(x) -ϱ(t)| ≤ |x -t|.
1. There exists α ∈ [1/2, 1] and an even function φ ∈ Γ 0 (R) such that

(∀x ∈ R) ϱ(x) = x + 2α(prox φ (x) -x). (17) 2.
Let

(∀x ∈ R) φ(x) = φ(x) + x 2 2 ( 18
)
and let φ * ∈ Γ 0 (R) be the Fenchel-Young conjugate of φ. 2 Then

(∀x ∈ R) ψ(x) c = (1 -2α) x 2 2 + 2α φ * (x), (19)
where c = designates equality up to an additive constant.

3. If ψ is convex, then α = 1/2.
Proof. See Section S1 of our supplementary material.

The interplay between penalty function ψ, the activation ρ, and the underlying convex function φ is summarized in Fig. 2. The previous proposition ensures that, when ψ is a convex penalty function satisfying Assumption 1, there exists φ ∈ Γ 0 (R) such that

ϱ = prox φ . ( 20
) φ ϱ ψ (17)
( 1 9 ) D e r i v a t i o n In this case, ϱ enjoys the property of being firmly nonexpansive, i.e. increasing and 1-Lipschitz [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Proposition 4.4] . As emphasized in [START_REF] Combettes | Lipschitz certificates for layered network structures driven by averaged activation operators[END_REF], it turns out that most of the existing activation functions used in neural networks are proximity operators of convex functions (e.g. ReLU activation function). On the other hand, if ψ is non convex, then [START_REF] Karunanithy | FID-Net: A versatile deep neural network architecture for NMR spectral reconstruction and virtual decoupling[END_REF] shows that ϱ is the overrelaxation of the proximity operator of a function in Γ 0 (R). According to (E1) in Section S1 of the supplementary material, ϱ is then α-averaged, which means that it satisfies the inequality [71, Proposition 4.35] [72],

(∀(x, t) ∈ R) 2 |ϱ(x) -ϱ(t)| 2 ≤ |x -t| 2 - 1 -α α |x -ϱ(x) -t + ϱ(t)| 2 . ( 21
)
The smallest α, the tightest the upper bound above. In other words, a small value of α induces more stability. Equation [START_REF] Wu | Improvement in signal-to-noise ratio of liquid-state NMR spectroscopy via a deep neural network DN-Unet[END_REF] gives another interpretation of α. It shows that ψ is a weakly-convex function [73, Proposition 4.3] with modulus 2α-1. Thus, α can be viewed as a measure of the degree of nonconvexity of potential ψ. To conclude, the activation functions in our unrolled architecture benefit from specific stability properties inherited from their optimization-based interpretation. This is in line with recent works [START_REF] Combettes | Lipschitz certificates for layered network structures driven by averaged activation operators[END_REF][START_REF] Alt | Translating diffusion, wavelets, and regularisation into residual networks[END_REF].

In [START_REF] Alt | Translating diffusion, wavelets, and regularisation into residual networks[END_REF], diffusion methods, wavelet approaches, and variational regularization are employed to perform a denoising task. The link with CNN architectures is made through residual neural blocks. Diffusitivies, shrinkage functions, and regularizers are reinterpreted as activation functions. The authors emphasize the prospects of nonmonotone functions, rarely used in the deep learning community. Through our method, we obtain a dictionary of novel activation functions, associated to usual convex and non-convex sparsity promoting penalizations. We provide, the values for α in the last column of Table 1. The convex function φ associated to a penalty ψ is implicitly defined by ( 18) and [START_REF] Wu | Improvement in signal-to-noise ratio of liquid-state NMR spectroscopy via a deep neural network DN-Unet[END_REF]. As shown in Section S3 of our supplementary material, it is possible for most of the convex penalties in Table 1 to derive a closed form expression for φ.

Penalty Name Penalization ψ(t), t ∈ R Activation ϱ(t), t ∈ R α Convex penalties Fair potential δ(|t| -δ log( |t| δ + 1)) δt |t|+δ 1 2 Huber t 2 2 if |t| < δ δ|t| -δ 2 /2 otherwise t if |t| < δ δsign(t) otherwise 1 2
Green log(cosh(t)) tanh(t)

1 2
Hyperbolic

δ 2 κ (1 + t 2 δ 2 ) κ 2 -1 t(1 + t 2 δ 2 ) κ 2 -1 1 2

Nonconvex penalties

Welsch

δ 2 1 -exp(-t 2 2δ 2 ) t exp(-t 2 2δ 2 ) exp(-3 2 ) + 1 2 Geman-McClure δ 2 t 2 2δ 2 +t 2 4δ 4 t (2δ 2 +t 2 ) 2 5 8 Tukey biweight δ 2 (1 -(1 -t 2 6δ 2 ) 3 ) if |t| ≤ √ 6δ 1 otherwise t(1 -t 2 6δ 2 ) 2 if |t| ≤ √ 6δ 0 otherwise 9 10
Hyperbolic tangent

δ 2 tanh( t 2 2δ 2 ) t (cosh( t 2 2δ 2 ) 2 0.9581 Cauchy δ 2 2 log(1 + t 2 δ 2 ) δ 2 t t 2 +δ 2 9 16
Table 1: Examples of penalties ψ satisfying Assumption 1. In the third column, we provide also the expression of their derivatives (i.e., nonlinear activation in our architecture) ϱ and, in the last column, the averaging constants α (see Section 3.2). All expressions are valid for every t ∈ R, (λ, δ) ∈]0, +∞[ 2 and κ ∈ [START_REF] Heaney | Mass spectrometry in medicine: a technology for the future?[END_REF][START_REF] Ernst | Principles of Nuclear Magnetic Resonance in One and Two Dimensions[END_REF]. The calculations for the averaging constants in non-convex cases are provided in Section S2 of our supplementary material.

Learning strategy

As aforementioned, the rewriting of HQ algorithm as a multi-layer neural network structure (see Fig. 1) paves the way for supervised learning of its hyperparameters. In this subsection, we showcase our proposed strategy to tune automatically such parameters.

First, regularization parameters (λ 1 , λ 2 ) are untied into two sequences (λ i,k ) 0≤k≤K-1 , i ∈ {1, 2}, varying with each layer of the proposed architecture. This is in contrast with classical optimization-based restoration procedures, which rely on setting those hyperparameters over all iterations, relying on empirical or statistical tuning techniques often difficult to implement [START_REF] Bauer | Comparing parameter choice methods for regularization of ill-posed problems[END_REF].

Here, we propose instead to learn automatically data-driven parameters. For a given number of layers K > 0, each layer L y k indexed by k ∈ {0, . . . , K -1} is attached to learnable units defined as

(∀i ∈ {1, 2}) λ i,k = ReLU FC i,k ([Hx k -y] ℓ ) 2 ) 1≤ℓ≤m ) . (22)
Hereabove, ReLU is the rectified linear unit [START_REF] Zeiler | On rectified linear units for speech processing[END_REF] applied componentwise, aiming at enforcing positivity on regularization parameters, and defined as follows:

(∀z ∈ R m ) ReLU(z) = (max(0, [z] ℓ )) 1≤ℓ≤m . (23)
Moreover, FC i,k , i ∈ {1, 2}, k ∈ {0, . . . , K -1}, are fully connected layers acting on component-wise squares of the data fit residual at the entry of each layer L y k . In such way, the regularization parameters automatically adapt to the noise level at hand, without any required prior knowledge. A similar strategy was employed in [START_REF] Bertocchi | Deep Unfolding of a Proximal Interior Point Method for Image Restoration[END_REF] in the context of image restoration. The rationale behind this choice is derived from the Bayesian interpretation of the objective function, which leads to having regularization hyperparameters proportional to the Gaussian noise variance.

Second, we suggest to also learn, in a supervised fashion, the stepsize sequence (γ k ) 0≤k≤K-1 involved in the linear weight layer sequence (W k ) 0≤k≤K-1 . The rationale is to tune optimally the contribution of each layer for delivering the best estimate x at the output of the unrolled architecture. In order to reduce spurious local minima in the training phase, we adopt an overparametrization strategy [START_REF] Bell-Kligler | Blind super-resolution kernel estimation using an internal-GAN[END_REF], that consists in defining

(∀k ∈ {0, . . . , K -1}) γ k = ReLU(c k,1 × c k,2 × . . . × c k,N ), (24)
with N ≥ 1 (typically, N = 10), and (c k,j ) 0≤k≤K-1,1≤j≤N are scalars quantities that are learned. Vector θ = (θ k ) 0≤k≤K-1 ∈ R K(2m+N ) gathers the parameters to be learnt, namely 2mK weights for the fully-connected layers in [START_REF] Haar | An analysis of explainability methods for convolutional neural networks[END_REF], and N K values to learn the stepsize quantities as in [START_REF] Gupta | Multivariate lipschitz analysis of the stability of neural networks[END_REF]. Simple constraints, such as range or sum-to-one, on the estimates can be easily imposed by appending a dedicated activation function at the end of the unrolled architecture. Such strategy will be employed in our experimental section, using a ReLU as last layer, to enforce the positivity of the restored MS spectra. The global multi-layer architecture, called U-HQ (Unrolled Half-Quadratic), is summarized in Fig. 3, with a detailed overview of a layer L θ k k structure displayed in Fig. 4, responsible for learning regularization and stepsize hyperparameters to be fed to the layer L y k simulating one iteration of the HQ scheme, as shown in Fig. 1. Given a training set {(x (i) , y (i) ), i ∈ {1, . . . , S}} comprised of S pairs of groundtruth and degraded data samples, the prediction feedforward model for one instance i ∈ {1, . . . , S} is given by

f θ (x (i) 0 ; y (i) ) = L y (i) K-1 • • • • • L y (i) k • • • • • L y (i) 0 (x (i) 0 ), (25)
x (i) 0 denoting the input for the first layer of the network for the i-th sample. Then, the learnt hyperparameter vector θ is defined as the solution to the optimization problem minimize

θ∈R K(2m+1) 1 S S i=1 L(f θ (x (i) 0 ; y (i) ), x (i) ), ( 26
)
with L a predefined loss, for instance, the mean squared error. Problem ( 26) can be solved end-to-end through backpropagation with standard neural network training algorithms such as stochastic gradient or ADAM [START_REF] Kingma | A method for stochastic optimization[END_REF]. At the inference stage, for a given degraded signal y, and an initial prediction x 0 , the inferred output is given by x = f θ (x 0 ; y). The proposed approach benefits from automatic parameter tuning with layer-dependent and data driven structure, while having an interpretable structure.

Restoration of mass spectrometry signals

We now evaluate the performance of U-HQ on the restoration of MS sparse signals. We focus on the resolution of inverse problem (1) where x ∈ R n is a ground truth signal with length n = 2000, H ∈ R m×n is a Toeplitz matrix mimicking a circular padded convolution with a blur kernel, and e ∈ R m corresponds to an additive i.i.d. Gaussian noise with zero-mean and standard deviation σ > 0.

Datasets

Let us first introduce the three datasets we employ to evaluate our method as well as its competitors. In all cases, the dataset is made of original signals x paired with degraded signals y obtained through Model [START_REF] Heaney | Mass spectrometry in medicine: a technology for the future?[END_REF]. Moreover, the blur support size is 50, so that m = 2049. Dataset 1. Signals x are real MS spectra extracted from the mass bank inventory of MS spectra [78], considering a monoisotopic mass range of [0,200] ppm and a charge state +1. The spectra intensities are rescaled to reach a maximum value of 100. The blur kernel varies for each signal. It has a "Ricker" shape with spreading factor chosen randomly uniformly in (0.25, 1) [START_REF] Ricker | The form and nature of seismic waves and the structure of seismograms[END_REF]. The noise standard variation is also varying throughout the dataset, and is chosen randomly uniformly between 0 and 0.5. The blur kernel is known, while the noise level is assumed to be unknown. Training, validation, and test sets contain 900, 100 and 100 pairs (x, y), respectively. Dataset 2. Signals x are also real MS spectra extracted from the same mass spectrometry bank inventory [78], we consider again the monoisotopic mass range [0,200] ppm and a charge state +1. Hereagain, the spectra intensities are rescaled to reach a maximum value of 100. The blur kernel is chosen as for Dataset 1. The noise standard deviation now varies uniformly between 0.5 and 1. Training, validation, and test sets contain 900, 100 and 100 pairs (x, y), respectively. Dataset 3. Groundtruth signals x have a monoisotopic mass between 0 and 200 ppm, similarly to Dataset 1 and Dataset 2. We now use blur kernel with shape "Fraser Suzuki" [START_REF] Felinger | Data analysis and signal processing in chromatography[END_REF], whose spreading factor is chosen randomly uniformly between 0.25 and 1 and the asymmetry factor between 0.2 and 0.6. The noise standard deviation varies between 0 and 0.5. Training, validation, and test sets contain 900, 100 and, 100 pairs (x, y), respectively.

Experimental settings

Evaluation of the results. All compared supervised models are trained so as to solve [START_REF] Combettes | Lipschitz certificates for layered network structures driven by averaged activation operators[END_REF], where L is the averaged mean squared error (MSE) between the estimated signal x and the ground truth x, MSE(x, x) = (1/n)∥x-x∥ 2 . MSE metric is also used to finetune regularization hyperparameters, as explained in a forthcoming section. In our comparative tables, we provide the signalto-noise ratio (SNR, in dB) defined as SNR(x, x) = 20 log 10 (∥x∥/∥x -x∥) , as well as the thresholded SNR metric (TSNR) in dB, defined as the SNR computed only on the support (i.e., indices of non-zero entries) of the sought signal x. We display both the mean and standard deviation (STD) of these quantitative metrics, computed over the test sets.

Penalties. The regularization functions Ψ 1 and Ψ 2 in (2), are defined as in (3), with ψ 1 the (convex) Fair potential and ψ 2 the (non-convex) Cauchy penalty, whose mathematical expressions and properties are listed in Table 1. They depend on smoothing parameters δ 1 and δ 2 , respectively, which are set through a procedure precised hereafter. Within an ablation study, we will also compare to cases when only one of these penalty functions is employed (i.e., when we remove the network branch related to either Ψ 1 or Ψ 2 ). Results when using a quadratic (i.e., Tikhonov) penalty for ψ 1 (while ψ 2 = 0) will also be presented. This experimental plan allows, in particular, to study the role and effect of the layer W k , and the activations R 1 and R 2 , in our U-HQ architecture.

Initialization. All the compared approaches are initialized by the null vector, in particular, all the compared neural networks take x 0 = 0 as an input.

Programming framework. The methods are implemented using Pytorch (version 1.8.0) under Python (version 3.8.5). The GPU-based code is run on an Nvidia DGX server using a V100 SXM2 GPU card (32 GB of RAM). The CPU is an Intel@ Xeon(R) W-2135 CPU @ 3.70 GHz with 12 active cores.

Comparison settings

The proposed U-HQ approach is compared to various methods, that can be classified into three categories, namely (i) optimization-based methods which rely on the HQ scheme to minimize (2), (ii) deep unrolling techniques relying on similar or distinct architecture than ours, with various training strategies, and (iii) deep neural networks inspired from state-of-the-art architectures deployed for the restoration of sparse and/or spectroscopic signals.

Let us specify hereafter the tuning for all the competing methods, including ours.

Optimization-based methods. We implement the HQ algorithm (7) described in Section 2.2 to minimize (2). We run HQ until convergence, that is either until it reaches a maximal number of iterations (here, 100), or until it shows a stabilization of its iterates through ∥x k+1 -x k ∥ 2 /∥x k ∥ 2 < 10 -5 . We denote HQ-SC this approach. We also consider a modified variant of HQ algorithm, where the number of iterations is fixed to a value K identical to the number of layers retained in the unrolled architecture. This method is inspired from early stopping (ES), therefore we name it HQ-ES.

For both methods, the stepsize is set to γ k ≡ 1 for simplicity. Moreover, the hyperparameters (λ 1 , λ 2 ) are finetuned manually. Specifically, a gridsearch is performed on 200 signal randomly extracted from the training set to minimize the averaged MSE. The resulting parameters are retained later to assess the performance of these methods on the test set. Specifically, we perform the gridsearch using the convex and non-convex potentials separately to finetune both the regularization and smoothness hyperparameters (λ 1 , λ 2 , δ 1 , δ 2 ). Then, we run HQ algorithm using the hybrid penalty, with (δ 1 , δ 2 ) fixed according to the previous disjoint gridsearch and we finetune the regularization hyperparameters (λ 1 , λ 2 ). The grid used for λ 1 and λ 2 is [0, 0.01, 0.1, 1, 2, 4, 6, 8, 10, 20], while for δ 1 and δ 2 it is [0.01, 0.1, 0.5, 1, 2, 4, 6, 8].

Learning-based methods. In order to assess the benefits of the proposed architecture in U-HQ, we provide comparisons with simplified versions of it, with an increasing number of free parameters, so as to perform an ablation study. We start with a low parametrized model inspired from deep equilibrium (DE) models [START_REF] Gilton | Deep equilibrium architectures for inverse problems in imaging[END_REF]. Namely, the stepsize is set as γ k ≡ 1, and the trained architecture shares the same two parameters (λ 1 , λ 2 ) for all layers. So, for a fixed number of layers K, we have

L (θ 0 ) 0 = L (θ 1 ) 1 = . . . = L (θ K-1 ) K-1 = L (θ) . ( 27
)
This yields the so-called U-HQ-DE method. We then implement two architectures entitled U-HQ-FixS (U-HQ Fix Stepsize) and U-HQ-FixN (U-HQ Fixed Noise). In both cases, the noise estimation layer from Figure 4 is discarded, and the regularization hyperparameters (λ i,k ) 0≤k≤K-1 , are simply defined as outputs of ReLU unit with 2K learnable inputs (a i,k ) 0≤k≤K-1 , i ∈ {1, 2}. Moreover, in U-HQ-FixS, we furthermore impose γ k ≡ 1, so that we learn only (λ i,k ) 0≤k≤K-1 , while U-HQ-FixN learns stepsize hyperparameters (γ k ) 0≤k≤K-1 with enforced positivity using ReLU. Finally, we implement U-HQ-FixN-OverP. This last architecture uses U-HQ-FixN as a baseline, and includes the overparametrization strategy similar as the one described in [START_REF] Gupta | Multivariate lipschitz analysis of the stability of neural networks[END_REF], for the learning of both the regularization and stepsize hyperparameters per layer. The unrolled schemes U-HQ-FixS, U-HQ-FixN, U-HQ-FixN-OverP, and the proposed U-HQ, are experimented using either ψ 1 (convex) or ψ 2 (non-convex), or (ψ 1 , ψ 2 ) (hybrid), in their activation layers. Note that using only one of these two penalizations practically results in training half the number of weights. Regarding U-HQ-DE we present the results only for the hybrid penalty pair (ψ 1 , ψ 2 ), for the sake of conciseness. Furthermore, we make comparisons with two unrolled methods recently introduced for sparse signal recovery in [START_REF] Gharbi | Unrolled primal-dual deep network for sparse signal restoration[END_REF], namely U-ISTA and U-PD, based on an unrolling of the ISTA algorithm [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm with application to wavelet-based image deblurring[END_REF] and the primal dual algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], respectively. Both methods rely on an ℓ 1 (thus convex) penalization on the signal. Regularization and stepsize parameters of both methods are learned, with enforced positivity using ReLU activations.

Deep learning methods. We compare our proposed approach U-HQ to three deep neural network architectures. First, we implement a fully connected network FCNet, based on the concatenation of 4 blocks of layers, each composed of a linear map and a ReLU activation layer. Second, we create an autoencoder structure, inspired from [START_REF] Brandt | Deep learning for reconstructing low-quality FTIR and Raman spectra a case study in microplastic analyses[END_REF], named AE. Similarly to FCNet, the AE composes linear maps and ReLU layers. The difference is that the linear maps are such that the encoder reduces, with a factor 2 per layer, the dimension of the input then the decoder maps, through up-samplings with factor 2, the low-dimensional signals to the output. We set 4 layers in both encoder and decoder parts. We finally deploy a handcrafted Residual UNet, named ResUNet, reminiscent from the one proposed in [START_REF] Horgan | High-Throughput molecular imaging via deep-learning-enabled Raman spectroscopy[END_REF], in the context of Raman spectra recovery. The latter architecture relies on an autoencoder structure like UNet [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF]. We adapt the code available in the github of [START_REF] Horgan | High-Throughput molecular imaging via deep-learning-enabled Raman spectroscopy[END_REF], by reducing the number of parameters of the network, to limit over-fitting effects. We provide the details of the resulting ResUNet architecture in Section S6 of the supplementary material as well as the two comparative deep learning methods.

The acronyms of all the compared approaches are summarized in For reproducibility purposes, we share a repository available at https:// github.com/GHARBIMouna/Unrolled-Half-Quadratic with our implementation in PyTorch of all the compared algorithms.

Results

Table 4 summarizes the results in terms of mean and standard deviation of the SNR and TSNR scores computed on test sets, for the three datasets. We recall that the higher the SNR (resp. the TSNR), the better the estimation. Best quantitative results are highlighted in bold style. At a first glance, this confirms the competitiveness of our proposed architecture U-HQ with respect to all other compared methods. Let us inspect and discuss more deeply the results.

U-HQ vs optimization-based methods. First, let us inspect the results of optimization-based methods. Using the proposed hybrid penalty term with iterative methods HQ-SC and HQ-ES reaches better metrics compared to using convex or non-convex potentials separately in most cases on all datasets. This demonstrates the advantage of a hybrid penalty in this applicative problem, as it brings flexibility in the regularization model. However, it is at the price of a tedious grid search. Second, from a general glance, we can notice that optimization-based methods are outperformed by most learning-based ones. This is because learning-based competitors perform a supervised data-driven tuning of its hyperparameters to optimize the quality metric (here, MSE, which is directly related to SNR). In particular, our proposed approach U-HQ yields better restoration quality compared to the use of a standard iterative HQ algorithm HQ-SC, and to its early stopped variant HQ-ES (up to 3dB improvement in both SNR/TSNR).

Computational complexity and stability of U-HQ. The interpretation of HQ algorithm as layers of a feed-forward neural network has allowed us to implement it efficiently on the GPU-based Pytorch environment, which considerably reduces its execution time, when compared to a more standard CPU implementation. Since it requires, in average, a lower number of layers/iterations than its optimization-based counterpart, U-HQ benefits from a reduced execution time compared to HQ-SC, as illustrated in Table 3. For the sake of comparison, we also display the test time of the second best performing method, U-HQ-FixN (i.e., U-HQ without noise estimation layer). One can notice that this architecture has similar test time than U-HQ, which shows that introducing data-driven fully connected layers (see Fig. 7) does not increase the network complexity. On top of its good qualitative and computational performance, the training stability of U-HQ, and the absence of overfitting phenomenon can be assessed on the training loss curves on Dataset 1 given by Fig. 5. An example of restored signal in Dataset 1 is displayed in Fig. 6. We refer the reader to Section S5 of our supplementary file for similar plots on Datasets 2 and 3.

Ablation study on U-HQ. Third, the proposed architecture U-HQ appears superior to all its simplified variants, namely U-HQ-DE, U-HQ-FixS, U-HQ-FixN, and U-HQ-FixN-OverP. This illustrates the importance of the learnable stepsize and hyper-parameter sequences, as well as the benefits of overparametrization. A visual representation of learnt parameters is given in Fig. 7 for Dataset 1 (see Section S5 of our supplementary file for the other datasets). Interestingly, several regimes can be observed, namely boosting phases with higher regularization parameters and higher stepsizes, and stabilization phases where all parameters stabilize to constant values. Table 3: Averaged iteration/layer numbers and execution time on GPU per signal, for inference on test set on Dataset 1 to Dataset 3, using HQ-SC,U-HQ-FixN, U-HQ-FixN-OverP, U-HQ and DL methods (FCNet, AE, and ResUNet, respectively). For the former, we indicate also CPU times, for the sake of comparison.

Finally, it can be seen that, like for optimization-based methods, employing both penalty branches in the architecture of U-HQ has clear benefit on the performance when compared with the use of only the convex or non-convex term. The choice of good penalty terms is crucial to get high performance. A large class of alternative penalties, such as the Tikhonov-based penalty used in NMR spectroscopy in [START_REF] Venkataramanan | Solving Fredholm integrals of the first kind with tensor product structure in 2 and 2.5 dimensions[END_REF] for instance, can be easily encompassed in our framework. However, as our datasets are made of sparse signals, such penalty is not well adapted, and considerably degrades the quality of the results, as it can be seen in the table. The aforementioned observations confirm our choices for the final U-HQ architecture with a hybrid convex/non-convex penalty model.

U-HQ vs state-of-the-art benchmarks. The proposed U-HQ appears superior, in terms of both qualitative metrics, with respect to its competitors U-PD and U-ISTA, which again shows the benefit of considering a flexible and learnable regularization strategy in this applicative context. Furthermore, as shown in the bottom lines of Table 4, the deep learning models FCNet and AE perform very poorly on this task. Our approach displays consistently better restoration scores compared to with ResUNet. In terms of inference time, DL methods are however faster than the other competitors.

The main reason is that U-HQ and its variants, require at each layer (i.e., each iteration in its iterative form), the inversion of an n × n matrix to build the operator W k .

Discussion

Let us now take a step back, to discuss the strengths of the proposed unrolled architecture U-HQ, and the remaining ways for improvement. A key feature of unrolled architectures is that they are model-based. Specifically, the proposed approach U-HQ explicitly integrates the knowledge of the degradation operator H and of the statistical noise model, as well as the prior knowledge on the sought signals. This is clearly beneficial, as it brings interpretability, flexibility, and good performance to our model, compared to the DL methods. U-HQ is able to perform consistently well, even though our datasets are very challenging, with varying shapes for H, varying signal sparsity ratios, and varying noise levels, and a relatively low number of training samples. First, an avenue for future research is to go one step further, by considering an even more flexible architecture, for instance able to learn the degradation operator H (i.e., blind deconvolution context), or to cope with non-Gaussian noise. This would certainly require complexifying the problem formulation and the associated iterative algorithm method, for instance as in [START_REF] Huang | Unrolled variational Bayesian algorithm for image blind deconvolution[END_REF], and larger datasets for training. Second, as illustrated in the experiments, the performance of U-HQ depend on the choice for the potential functions retained, to build the activations R 1 , R 2 , as well as the linear block W k . The hybrid strategy allows to combine optimally two classes of penalties, and thus a wide class of priors on the sought signal. In a future work, we plan to investigate a higher number of penalty branches, with the aim to encompass a broader class of datasets encountereed in other applications of signal processing (e.g., seismic data [START_REF] Repetti | Euclid in a Taxicab: Sparse blind deconvolution with smoothed ℓ 1 /ℓ 2 regularization[END_REF]). Third, as discussed in the introduction, the unrolling paradigm paves the way for a thorough stability analysis, using advanced tools from fixed point analysis [START_REF] Valle | Stability of unfolded forward-backward to perturbations in observed data[END_REF]. Such analysis remains an open question for U-HQ, and would be an interesting future work perspective.

Conclusion

In this work, we propose a novel approach to solve sparse recovery problems, motivated by a class of inverse problems arising in chemistry. The proposed method relies on merging a model-based half-quadratic iterative method and a data driven supervised approach, by adopting the unrolling paradigm. Unrolling the half-quadratic algorithm leads to an original interpretable deep network architecture. Leveraging on automatic differentiation allows us to efficiently tune the algorithm hyperparameters. Furthermore, the resulting architecture turns out to be efficient in terms of reconstruction quality and computational budget, both of which are important features in any experimental setting. Our simulations illustrate the capability of our unrolling technique on various realistic databases of mass spectrometry spectra degraded by different blur kernels and noise levels.

Abstract

This supplementary file is organized as follows. In Section S1, we present our proof for Proposition 1. In Section S2, we provide the calculation details for α-averaging constants of the activation functions ϱ of Table 1. Then, in Section S3, we derive the expression of the convex function φ which was introduced in Proposition 1. We finally provide our training settings in Section S4 and display in Section S5 additional figures, complementing the Section 4 of the paper.

S1. Proof of Proposition 1

• If ϱ is 1-Lipschitz then there exists α ∈ [ 1 2 , 1] and a 1-Lipschitz function Since ψ is even by assumption, its derivative ϱ is odd, and we deduce from ( 17) that prox φ is odd too. For every p ∈ R, let u ∈ ∂φ(p) and let x = u + p. Then

ϑ such that (∀x ∈ R) ϱ(x) = (1 -α)x + αϑ(x). ( E1
x -p ∈ ∂φ(p) ⇔ p = prox φ (x) ⇔ -p = prox φ (-x) ⇔ -u = -x + p ∈ ∂φ(-p), (E2)
where ∂φ denotes the subdifferential of φ (see [R3]). This shows that ∂φ is odd.

According to [R3, Proposition 16.17],

∅ ̸ = int(dom φ) ⊂ dom ∂φ, (E3)
where

dom φ = {x ∈ R | φ(x) < +∞} (E4) dom ∂φ = {x ∈ R | ∂φ(x) ̸ = ∅}, (E5)
and int(S) denotes the interior of a set S. There thus exists p ∈ R and ϵ > 0 such that

p ∈ int(dom φ) ⊂ dom ∂φ = dom ∂ φ (E6) ⊂ dom φ, (E7)
where φ : x → φ(-x). Equality in (E6) follows from the fact that, for every p ∈ R, ∂ φ(p) = -∂φ(-p) = ∂φ(p) since ∂φ is odd. Then, we deduce from (E7) that p ∈ int(dom φ) ∩ dom φ.

Let us consider the even part φ e of φ defined as

φ e = 1 2 (φ + φ). ( E8
)
Since we have shown that int(dom φ) ∩ dom φ is nonempty, it follows from [R3, Theorem 16.47] that

∂φ e = 1 2 ∂φ + ∂ φ = ∂φ. (E9)
As φ e and φ have the same subdifferential, prox φe = prox φ . This shows that, in [START_REF] Karunanithy | FID-Net: A versatile deep neural network architecture for NMR spectral reconstruction and virtual decoupling[END_REF], we can assume that φ is even.

• It follows from [R3, Corollary 24.5] that φ * is a differentiable function such that prox φ = ( φ * ) ′ (E10)
where (•) ′ denotes the derivative of the function in argument. We deduce from ( 17) that

(∀x ∈ R) ψ ′ (x) = ϱ(x) = (1 -2α)x + 2α( φ * ) ′ (x). (E11)
By integrating, this yields

(∀x ∈ R) ψ(x) = (1 -2α) x 2 2 + 2α φ * (x) + C, (E12)
where C ∈ R.

• If ψ is convex, its derivative ϱ is increasing. Since ϱ is defined on R, it follows from [R2, Proposition 3.1(iii)] that α = 1/2.

S2. α-averaging constants of activation functions

In this section, we present calculation details for α-averaging constants of the activation functions ϱ of Table 1, as defined in Proposition 1. Firmly nonexpansive activations ϱ have an α constant equal to 1 2 . This applies to the first five rows of Table 1. Let us now calculate α for the non-monotonic activation functions ϱ listed in lines 5 -10 of Table 1

. It follows from (E1) that α ∈]1/2, 1] is the smallest constant such that ϑ = Id +α -1 (ϱ -Id) is 1-Lipschitz, i.e., it satifies the condition sup t∈R |ϑ ′ (t)| = 1 (E13)
• Welsch:

It follows from the expression of Welsch activation function that (E13) is equivalent to

sup t∈R 1 + α -1 exp -t 2 2δ 2 1 - t 2 δ 2 -1 = 1. (E14)
The zeros of the derivative of ϑ ′ are t = 0, t = √ 3δ, and t = -√ 3δ.

(E15)

Condition (E14) is then equivalent to checking the latter two values, namely

2 exp(-3/2) + 1) α -1 ≤ 1 (E16)
The smallest averaging constant is thus α = exp(-3/2) + 1/2.

• Geman-McClure:

From the expression of Geman-McClure activation function (E13) reads

sup t∈R 1 + α -1 4δ 4 (2δ 2 -3t 2 ) (2δ 2 + t 2 ) 3 -1 = 1. (E17)
The zeros of the derivative of ϑ ′ are t = 0, t = √ 2δ, and t = -√ 2δ

Thus, the maximum in (E17) is reached for the latter two values, which yield the minimum value of the averaging constant α = 5/8.

• Tukey biweight:

By using the piecewise form of the activation function, plugging in the expression of ϱ on [-√ 6δ, √ 6δ], and deriving ϑ, we obtain the equivalent formulation of (E13):

sup t∈[- √ 6δ, √ 6δ] 1 + α -1 1 - t 2 6δ 2 1 - 5t 2 6δ 2 -1 = 1 (E18)
The maximum is reached for |t| = 3 2 5 δ and the averaging constant of Tukey biweight activation function is equal to α = 9/10.

• Cauchy:

Condition (E13) reads sup t∈R 1 + α -1 δ 2 (δ 2 -t 2 ) (δ 2 + t 2 ) 2 -1 = 1,
the maximum being reached when |t| = √ 3δ. This results in α = 9 16 .

• Hyperbolic Tangent:

Condition (E13) is equivalent to sup t∈R 1 + α -1 δ 2 cosh( t 2 2δ 2 ) -2t 2 sinh( t 2 2δ 2 ) δ 2 cosh 3 ( t 2 2δ 2 ) -1 = 1. (E19)
The zeros of the derivative of ϑ ′ satisfy 2t 3t 2 sinh 2 ( t 2 2δ 2 ) -3δ 2 cosh(

t 2 2δ 2 ) sinh( t 2 2δ 2 ) -t 2 cosh 2 ( t 2 2δ 2 ) = 0.
The maximum in (E13) is obtained for the nonzero solutions to this equation. To determine the latter values, a numerical resolution is required. By using Newton-Raphson's method, we obtain numerically that the solutions are t/δ ≃ ±1.5355. We deduce that α ≃ 0.9581.

S3. Associated convex functions

In this appendix, we derive the expression of the convex function φ which was introduced in Proposition 1.

• Fair potential: By using ( 20) and (E10), we have Therefore, we deduce from (18) that, for every t ∈ R,

φ(t) = -δ|t| -δ 2 log(δ -|t|) -t 2 2 if |t| < δ +∞ otherwise (E26)
• Huber:

(∀t ∈ R) ϱ(t) = t if |t| < δ δ sign(t) otherwise = proj [-δ,δ] (t) (E27)
where proj S denotes the projection onto a nonempty closed convex set S. Thus

(∀t ∈ R) φ(t) = ι [-δ,δ] (t) = 0 if t ∈ [-δ, δ] +∞ otherwise.
(E28)

• Green: (E30)

(∀t ∈
• Hyperbolic:

We will just consider the cases when κ ∈ {1, 2} for which closed form expressions can be obtained.

Figure 1 :

 1 Figure 1: Unrolled HQ architecture: Layer L y k mimicks the k-th iteration of the HQ algorithm (7), for k ∈ {0, . . . , K -1}, and K > 0, as described in (14). V 0 , V 1 , V 2 and V 3 and W k are weight operators while b 0 , b 1 , b 2 , b 3 are biases. R 0 , R 1 , R 2 and R 3 are activation functions. Black dotted arrows showcase the dependency of the pointed block to the mentioned value. For instance, block V 3 requires the knowledge of (λ 1 , λ 2 ).

Figure 2 :

 2 Figure 2: Interplay between penalty ψ, activation ϱ and convex function φ.

Figure 3 :

 3 Figure 3: Supervised U-HQ network: For every k ∈ {0, . . . , K -1}, the detailed view of layer L y k is given by Fig. 1. Details of the learning block L θ k k are given in (22) and (24).

Figure 4 :

 4 Figure 4: Overview of one layer of proposed architecture U-HQ.

Figure 5 :

 5 Figure 5: Evolution of training loss (solid black) and validation loss (dashed blue) along epochs for U-HQ for Dataset 1.

Dataset

  FCNet GPU time (s) 3.42 × 10 -4 3.43 × 10 -4 3.28 × 10 -4 AE GPU time (s) 6.53 × 10 -4 7.89 × 10 -4 6.67 × 10 -4 ResUNet GPU time (s) 4.72 × 10 -3 4.55 × 10 -3 4.36 × 10 -3

Figure 6 :

 6 Figure 6: Groundtruth signal x (top left), degraded observation y (top middle) and reconstruction x using proposed method U-HQ (top right), ResUNet (bottom left), U-ISTA (bottom middle) and U-PD (bottom right) for one example from the test set of Dataset 1. x-axis indicates the mass-to-charge (m/z) ratio of chemical compounds and y-axis the abundance of the corresponding compound. SNR/TSNR scores of restored signal (in dB) are 31.75/34.96, 29.97/31.84, 23.08/25.50 and 26.13/29.50 for U-HQ, ResUNet, U-ISTA, and U-PD, respectively.

Figure 7 :

 7 Figure 7: Learnt hyperparameter sequences, averaged on test set, (λ 1,k ) 0≤k≤K-1 (left), (λ 2,k ) 0≤k≤K-1 (middle) and (γ k ) 0≤k≤K-1 (right) with respect to layer index k (x-axis) for the proposed method U-HQ, for Dataset 1.

)(

  In the worst case, one can simply choose α = 1.) A function ϱ satisfying (E1) is said to be α-averaged. The existence of φ ∈ Γ 0 (R) such that (17) holds then follows from [R2, Proposition 3.1(ii)].

(

  ∀t ∈ R) ϱ(t) = ( φ * ) ′ (t) = δt |t| + δ . (E20)This function can be seen as a scaled version of the Elliot activation function investigated in [R1,Example 2.15]. We have thus(∀t ∈ R) φ * (t) = δ(|t| -δ log(|t| + δ)) + η (E21) with η ∈ R. Since φ ∈ Γ 0 (R), (∀t ∈ R) φ(t) = φ * * (t) = sup u∈R tu -φ * (u). (E22)For every t ∈ R, let us define(∀u ∈ R) g(u) = tu -φ * (u) = tu -δ(|u| -δ log(|u| + δ)) -η. (E23) If t ≥ δ, then g(u) -→ u→+∞ +∞, and if t ≤ δ then g(u) -→ u→-∞+∞. Fort ∈] -δ, δ[, g reaches its maximum at u such thatt = δu |u| + δ ⇔ u = δt δ -|t| . (E24)For this value of u, by setting η = 2δ 2 log δ, g(u) = -δ|t| -δ 2 log(δ -|t|).(E25)

Tikhonov 8 K(m + N ) 10 - 5 , 6 ,

 81056 8 K(m + N ) 10 -5 , 15, 10 8 K(m + N ) 10 -5 , 15, 10 8 K(m + N ) 10 -m + N ) 10 -6 , 15, 10 8 2K(m + N ) 10 -6 , 15, 10 8 2K(m + N ) 10 -6 , 5, 5 10 -2 , 15, 10 10 -2 , 15, 10 10 -2 , 5, 5Table T1: Summary of training settings for experimented deep unrolled architectures: number of layers K, number of learnt parameters p, optimizer, learning rate (lr) and training bst and validation bsv batch sizes respectively. For U-HQ, we indicate lr λ (top) and lr γ (bottom).

Figure F1 :

 F1 Figure F1: Evolution of training loss (solid black) and validation loss (dashed blue) along epochs for U-HQ for Dataset 2 (left) and Dataset 3 (right).

Figure F2 :

 F2 Figure F2: Groundtruth signal x (left up), degraded observation y (middle up) and reconstruction x (right up) using proposed method U-HQ, ResUNet (left bottom), U-ISTA (middle bottom) and U-PD (right bottom), for an example from Dataset 2. x-axis indicates the mass-to-charge (m/z) ratio of chemical compounds and y-axis the abundance of corresponding compound. SNR/TSNR scores of restored signals (in dB) are 26.26/28.09, 25.05/26.23, 22.92/25.21, and 22.11/25.83 for the respective methods.

Figure F3 :

 F3 Figure F3: Groundtruth signals x (left), degraded observations y (middle) and reconstruction x (right) using proposed method U-HQ, ResUNet (left bottom), U-ISTA (middle bottom) and U-PD (right bottom) for an example from Dataset 3. x-axis indicates the mass-to-charge (m/z) ratio of chemical compounds and y-axis the abundance of corresponding compound. SNR/TSNR scores of restored signals (in dB) are 26.85/28.03, 28.67/29.83, 17.17/19.04 and 21.35/25.19 respectively.

Figure F4 :

 F4 Figure F4: Learnt hyperparameter sequences (λ 1,k ) 0≤k≤K-1 , (λ 2,k ) 0≤k≤K-1 and (γ k ) 0≤k≤K-1 with respect to layer index k (x-axis) for the proposed method U-HQ, for Dataset 2 (top) and Dataset 3 (bottom). Learnt parameters are averaged on test set.

Figure F5 :

 F5 Figure F5: ResUNet architecture. Borders of the input signal y are truncated, to reach an input size of n = 2000.

Table 2 .

 2

	Method's name	Acronym definition
	U-HQ	Unrolled Half Quadratic approach
	HQ-SC	Half Quadratic algorithm with Stopping Cri-
		teria
	HQ-ES	Half Quadratic algorithm with Early Stop-
		ping rule
	U-HQ-DE	Unrolled Half Quadratic approach with Deep
		Equilibrium model.
	U-HQ-FixS	Unrolled Half Quadratic approach with
		Fixed Stepsize
	U-HQ-FixN	Unrolled Half Quadratic approach with
		Fixed Noise model
	U-HQ-FixN-OverP Unrolled	Half	Quadratic	approach
		with Fixed Noise model using over-
		parametrization
	U-PD	Unrolled Primal Dual approach
	U-ISTA	Unrolled Iterative Soft Threshholding Algo-
		rithm
	FCNet	Fully Connected Network
	AE	AutoEncoder
	ResUNet	Residual UNet

Table 2 :

 2 Glossary of termsTraining settings. Let us now present the training specifications used for the learning-based methods described in the previous subsection. The number of layers is set, for each dataset and each method, through empirical search, to maintain a trade-off between a low validation loss, a comparable number of parameters between all tested methods, a reasonable computational training time and a stable learning. For architectures U-HQ, U-HQ-FixS, U-HQ-FixN, U-HQ-FixN-OverP and U-HQ-DE, we set (δ 1 , δ 2 ) = (0.1, 1) for Dataset 1 and Dataset 2 and (δ 1 , δ 2 ) = (0.01, 0.1) for Dataset 3 in the unrolled HQ methods (while these parameters were finetuned in optimization-based methods). We set initial values of learnable weights/parameters to 1, and the parameter N to 10 in all our overparametrized settings. Training is performed using Adam optimizer[START_REF] Kingma | A method for stochastic optimization[END_REF], with mini-batch size 5. Regarding the training strategy, we set a unique learning rate lr to update the weights in most architectures, at the exception of U-HQ. As this latter uses completely different architectures for the stepsize (i.

e., simple ReLU layer) and the regularization weights (i.e., fully connected layers), we use two different learning rates, namely lr γ and lr λ , respectively. Regarding the deep learning methods, Adam optimizer is also used to minimize the MSE loss. Training parameters, including learning rates, batch sizes, and layer number, for all architectures, are summarized in Sections S4 (for unrolled models) and S6 (for DL models), in our supplementary material.

Table 4 :

 4 Summary of results in terms of mean (std) of SNR/TSNR metrics, in dB, computed on test sets (see acronyms definitions in Table2). Best scores are highlighted in bold. Second best scores are underlined.

  An Unrolled Half-Quadratic Approach for Sparse Signal Recovery in Spectroscopy (supplementary) Mouna Gharbi a , Emilie Chouzenoux a , Jean-Christophe Pesquet a a CVN, CentraleSupélec, Inria Saclay, University Paris Saclay, Gif-sur-Yvette, France

Preprint submitted to Signal ProcessingDecember 22, 2023

https://fr.mathworks.com/matlabcentral/fileexchange/ 88897-spoq-smooth-sparse-p-over-q-ratio-regularization-toolbox

The conjugate of φ is defined as (∀u ∈ R) φ * (u) = sup x∈R xu -φ(x).

Preprint submitted to ElsevierDecember 22, 2023

Acknowledgments

This work has been supported by the ITN-ETN project TraDE-OPT funded by the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 861137 and by the European Research Council Starting Grant MAJORIS ERC-2019-STG850925.

We can thus set

By defining g as in (E23), if t ≥ δ, then g(u) -→ u→+∞ +∞ and, if t ≤ δ, then g(u) -→ u→-∞ +∞. For t ∈] -δ, δ[, the maximum of g is reached at

We deduce that, for every t ∈ R,

(E34)

• If κ = 2, then ϱ reduces to the identity function and φ = 0.

S4. Training settings for Section 4

We summarize in the Table T1 the training settings used for the unrolled methods experimented in the paper.

S5. Additional figures for Section 4

We display in Figure F1

S6. Settings of DL methods

In this section, we display the ResUNet architecture used for the sake of comparison on Fig. F5. A detailed overview of each of this network's blocks is provided in Table . T2. It is important to note that each residual convolutional block has a skip connection parallel to the convolutional branch. Moreover, we display the details of both AE and FCNet architectures in Table T3. Finally, we provide the DL methods training settings in Table T4.