A Cherni

E Chouzenoux

M-A Delsuc
email: delsuc@igbmc.fr

PALMA, an improved algorithm for DOSY signal processing

NMR is a tool of choice for the measure of diffusion coefficients of species in solution. The DOSY experiment, a 2D implementation of this measure, has proven to be particularly useful for the study of complex mixtures, molecular interactions, polymers, etc. However, DOSY data analysis requires to resort to inverse Laplace transform, in particular for polydisperse samples. This is a known difficult numerical task, for which we present here a novel approach. A new algorithm based on a splitting scheme and on the use of proximity operators is introduced. Used in conjunction with a Maximum Entropy and 1 hybrid regularisation, this algorithm converges rapidly and produces results robust against experimental noise. This method has been called PALMA. It is able to reproduce faithfully monodisperse as well as polydisperse systems, and numerous simulated and experimental examples are presented. It has been implemented on the server http://palma.labo.igbmc.fr where users can have their datasets processed automatically.

Introduction

Diffusion coefficients can be efficiently measured in NMR by the use of magnetic field gradients. The most classical approach consists in applying a symmetric pair of pulsed field gradients (PFG) of varying intensity, separated by a diffusion delay ∆. Random displacements of the molecule during ∆ because of Brownian motion result in modulation of the signal intensity I following the Stejskal-Tanner equation [1,2]:

I(q) = I o exp(-D∆q 2 ) ( 1
)
where D is the diffusion coefficient of the molecular species and q = γδg, the measure of the phase dispersion created by the PFG. Here, γ is the gyromagnetic ratio of the studied spin, and δ and g are the duration and intensity of the PFG respectively. A least squares fit of the experimental values to an exponential decay provides an estimate of the value of D.

The DOSY experiment, introduced by Johnson [3] is a representation of this measure as a 2D spectrum, with chemical shifts presented horizontally and diffusion coefficients vertically.

DOSY has been used intensively to analyse molecular interactions, to sort the component in complex mixtures, or to evaluate molecular size distributions [4,[START_REF] Price | NMR Studies of Translational Motion: Principles and Applications[END_REF][START_REF] Callaghan | Translational dynamics and magnetic resonance: principles of pulsed gradient spin echo NMR[END_REF].

A monodisperse sample presents a well defined diffusion coefficient, and a simple exponential adjustment of equation ( 1) allows the determination of D. When several compounds share the same chemical shift, resulting in overlapping lines in the NMR spectrum, the result of a mono-exponential fit becomes incorrect. A simple column-wise least squares fit to two or more exponentials presents instabilities in noisy data-sets which make this approach difficult to use on complex cases. Several methods have emerged in the literature for the analysis of complex mixtures of monodisperse species, where the difficult mostly arises from the presence in the spectrum of many overlapping species Approaches based on a global analysis of the whole experimental matrix have been proposed, based on a clever decomposition of the 2D spectrum matrix into multivariate model, which allows to extract the spectra of each species along with their respective diffusion profiles [START_REF] Stilbs | [END_REF]8]. Some developments on this approach have been based on a harmonic analysis of the decay [9,10,11,12,13,14]. The exponential hypothesis can even be relaxed using methods related to blind source deconvolution [15,16,17,18]. All these approaches model the sample as a mixture of species with a characteristic decay pattern. Some samples such as polydisperse polymers, dendrimers, nanoparticles, gels and aggregated species, because of the variation in size, length or aggregate state of the different molecules in the sample, present a distribution of diffusion coefficients rather than a single coefficient. Moreover, the presence of common decay patterns at different chemical shifts cannot really be assumed anymore, as in these complex systems a subtle coupling between the chemical shift and the size usually broadens the spectral line, with each spectral channel sampling a slightly different subset of the species, so that each spectral channel has to be processed independently. For these strongly polydisperse samples a precise determination of the diffusion distribution is of great analytical importance, it is however a difficult task. Because of this difficulty, polydispersity is commonly measured by a polydispersity index (PDI) defined as the ratio of the mass average molar mass M w to its number average molar mass M n : PDI = M w /M n . This quantity characterises the breadth of the distribution independently of the details of its shape, a PDI of 1.0 indicates monodispersity. PDI is commonly measured by size exclusion chromatography, (electron) microscopy, light or X-Ray scattering, or even NMR-DOSY [19].

A polydisperse sample has to be analysed with a distribution X(D) of diffusion coefficients and equation (1) becomes:

I(q) = Dmax D min X(D)exp(-D∆q 2 )dD (2)
Determining the distribution X(D) from I(q) requires to solve the Laplace inversion of the q 2 dependency of I(q). The shape of the distribution X can be modelled by a Gaussian function or by any other symmetric or asymmetric analytical shape, and the parameters for this shape fitted to the experimental data [20]. This straight-forward approach is very sensitive to the choice of the shape, and will fail if it is not well adapted to the data, or if the distribution contains several isolated massifs, and it should be used with care.

In this work, we present a general approach that solves the Laplace inversion problem presented in equation (2). A new algorithmic approach based on a splitting scheme and on the use of proximity operators is introduced. Used in conjunction with Maximum Entropy and 1 regularisations, the algorithm is stable against experimental noise, reproduces faithfully monodisperse as well as polydisperse situations, and converges rapidly.

Theory

Problem description

We assume that the diffusion experiment was performed over a series of M values of q (by varying δ, g or both) and measured as a series of intensities y m for a given chemical shift value. The problem stated by equation ( 2) can be discretised to be solved numerically:

y m = N n=1 x n exp(-D n ∆q 2 m )
with D n ranging from D min to D max . As this expression is linear in x n , it can be rewritten as follows:

Y = H H HX (3)
where ). In this work, we call X the Laplace spectrum of Y . Determining X from the knowledge of H H H and Y is an ill-posed problem as the experimental points are inevitably tainted with noise, and H H H is usually a non-square matrix, with N > M . A simple inversion does not provide a valid solution, and one has to resort to alternative approaches.

Y = {y m , 1 ≤ m ≤ M } is the experimental series, X = {x n , 1 ≤ n ≤ N }

Lagrangian formulation

A general approach for solving equation ( 3) is to generate a solution X that solves the following regularised minimisation problem: minimise

X∈R N H H HX -Y 2 + µΨ(X) (4)
The first term evaluates the distance between data and the reconstruction, while the second term is the regulariser, the Lagrangian coefficient µ > 0 acting as a weight between the two. The regularisation function Ψ is aimed at selecting among all possible distributions, the most natural one, given the experimental evidences, using some a-priori information on the problem. It is usually built as a measure of the cost of the reconstruction (in terms of energy, information, number of signals, etc. . . see below ) and tends to favour an empty spectrum. Depending on the expression chosen for Ψ, the problem can be solved by different approaches. The CONTIN method [21] solves this problem for Ψ(X) = ΓX 2 , where Γ is a matrix that contains prior assumptions about the data. Classical choices are Γ = Id which selects the solution with the least energy, or the first or second derivatives thus removing fluctuations not required for a faithful reconstruction. It has known a great success since its introduction more than 30 years ago, however it suffers from slow convergence and over-smoothed solutions. Choosing the opposite of the entropy as the regularisation function (Ψ(X) = (x n /a) log(x n /a)) allows to produce the distribution with the least information in the sense of Shannon [22]. This Maximum Entropy (MaxEnt) approach has shown to be of great efficiency and robustness for solving the DOSY problem [23] and has been widely used. However, because of the strong curvature of the entropy function, the classical implementations [24,23] of this approach are known to present slow convergence rates. Kazimierzuck et al [25] proposed recently to use Ψ(X) =

|x n | = X 1 where Ψ is the 1 norm of X. Their approach relies on previous works that have shown that this is equivalent to select the spectrum with the less non-null values. The principal advantage of their approach is to rely on recent major advances in the field of convex minimisation and compressed sensing. The algorithm ITAMeD they developed is based on the soft thresholding approach and allows a rapid convergence toward the solution. Enforcing a minimum number of non-null values in X is a good approach when the sample is a mixture of monodisperse compounds, and the Laplace spectrum a set of sharp lines, it is not appropriate however for the analysis of polydisperse samples that may present very large distributions. Urbańczyk et al. [26] recently extended this work to minimising the p norm (

X p = ( |x n | p ) 1 p ) with 1 ≤ p ≤ 2.
As expected, the authors show that adapting the p parameter to the kind of data allows to reconstruct spectra of various widths. However, the choice of p is somewhat ad hoc and has to be adapted to each situation, additionally the authors rely here on the IRLS [27] algorithm that is slower than ITAMeD. The TRAIn method [28] that has been proposed recently is not explicitly based on a regularisation approach, but rather on an early stopping strategy, in conjunction with a trust region algorithm. This method is claimed to be efficient on polydisperse distributions. The early stopping approach creates a bias in the final solution which can be assimilated to an implicit regularisation albeit with no analytical definition, and this could be seen as a deficiency [START_REF]Bayesian Approach to Inverse Problems[END_REF].

Constrained formulation

The Lagrangian coefficient µ involved in equation ( 4) may be difficult to adjust in practice. However, one has often some precise knowledge about the level of noise corrupting the data. We propose to adopt a more practical formulation by solving the following constrained optimisation problem: minimise

X∈R N Ψ(X) subject to H H HX -Y ≤ η (5)
where η > 0 is related to the expected quality of the fit, based on an estimate of the experimental noise. This has the advantage to shift the burden of determining the adequate value of a Lagrangian coefficient to the much simpler task of estimating a noise level.

Hybrid regularisation

In order to favour both smooth (polydisperse) and sparse (monodisperse) shapes in the estimated signal, we propose a novel regularisation defined as follows:

Ψ(X) = λ ent(X, a) + (1 -λ) 1 (X) (6)
where ent (X, a) is given by ent(X, a) =

     N n=1 xn a log xn a if x n > 0 0 if x n = 0 +∞ elsewhere,
and is the opposite of the Shannon entropy with a flat prior a > 0, typically chosen here from an estimate value of x n . 1 (X) is the 1 norm of the vector X defined as [START_REF] Elad | Sparse and Redundant Representations[END_REF]:

1 (X) = N n=1 |x n |
and λ ∈ [0, 1] allows to control the balance between the sparsity prior and the entropy prior.

Proximity operator

The general problem expressed in equation ( 5) can be elegantly solved by the convex optimisation algorithm PPXA+, based on the use of proximity operators [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF]. This allows a generic solver for any choice of convex function Ψ while insuring a very rapid convergence toward the optimal solution. The proximal approach is also the basis of the ITAMeD method. The FISTA algorithm on which it is based [START_REF] Beck | [END_REF] is derived from the soft-thresholding operator, the proximity operator of the 1 norm operator. To our knowledge, it has never been applied to the MaxEnt penalty nor to a hybrid approach as proposed here.

It can be shown [START_REF] Cherni | Proc. 8th Int. Symp. Signal, Image, Video and Commun[END_REF] that the proximity operator of the functional Ψ used in equation ( 6) can be expressed as follows (see ESI sections 1 and 2 † for details):

prox Ψ (X) = (p(x n )) 1≤n≤N (7)
where,

p(x n ) = λ a W a λ exp( axn-a(1-λ) λ + log(a) -1) if λ ∈]0, 1] sign(x n ) max (|x n | -(1 -λ), 0) if λ = 0 (8)
In the above expression, W states for the Lambert function, defined as the inverse function of f : z → ze z for all z ∈ C, [START_REF] Corless | [END_REF] i.e:

z = we w ⇔ w = W(z).
In the current context, only a restriction of W to R + is required. It should be noted that for pure Maximum Entropy (λ = 1) we recover the result from Combettes and Pesquet [START_REF] Combettes | Fixed-Point Algorithms for Inverse Problems in Science and Engineering[END_REF]:

p ent (x) = 1 a W (a exp(x) + log(a) -1)
Similarly, pure 1 regularisation (λ = 0) brings the soft thresholding operator:

p 1 (x) = sign(x) max (|x| -1, 0)

Algorithm

With the expression of the proximity operator given in equations ( 7)-( 8), the convex optimisation problem ( 5) can be easily solved using a proximal splitting algorithm. At each iteration, such a method alternates between the proximity operator of Ψ, and the proximity operator associated to the constraint H H HX -Y ≤ η, (i.e. the projection operator onto this constrained set).

In order to ensure good convergence properties of our algorithm, we adopt the PPXA+ approach from Pustelnik et al. [START_REF] Pustelnik | [END_REF], generalizing the PPXA method from Combettes and Pesquet [37]. These algorithms both rely on the Douglas-Rachford scheme [38], which consists in replacing the involved proximity operators by their reflections (see ESI sections 3, 4 and 5 † for details). This leads us to the so-called PALMA algorithm, standing for "Proximal Algorithm for L 1 combined with MAxent prior". This algorithm is fully detailed in the Electronic Supplementary Information †.

Material & Methods

Simulations.

Several simulated data-sets, chosen to represent various analytical situations, were used for the evaluation of the algorithm. Set A consists in three monodisperse components with diffusion coefficients 16 µm 2 /s , 63 µm 2 /s and 230 µm 2 /s, with respective intensities 1.0, 0.33, 0.66. This data-set is equivalent to the simulation used in Urbańczyk et al. [25]. Set B is a wide distribution, simulated as a log-normal distribution centred at 35 µm 2 /s, it presents a PDI estimated to 6.26. Sets C1 and C2 are asymmetric distributions built from 15 lognormal components, ranging from 18 to 85 µm 2 /s, with intensities ranging from 0.1 to 10, they have PDI estimated respectively to 1.79 and 1.32. In all simulations, 64 gradient values were simulated, and a Gaussian noise equal to 0.001%, 0.01%, 0.1%, or 1% of the initial point was added. The gradient values were chosen with a harmonic progression for set A, and with linear increments for sets B, C1, and C2. All Laplace spectra were reconstructed on 256 logarithmically sampled points. Other simulations with varying distributions and noise levels are also presented in ESI (see section 6 †).

NMR measure.

A set of PEO standards were purchased from American Polymer Standards Corporation (Mentor, OH, USA), and 3 samples were prepared. Sample a) is a standard PEO with M w = 2343.3 g mol -1 and PDI = 1.07; sample b) is a standard PEO with M w = 4051.2 g mol -1 and PDI = 1.28; sample c) is a mixture prepared from standard PEOs ranging from 350 to 5250 g mol -1 for a theoretical M w of 3238.5 g mol -1 and a theoretical PDI of 2.01. Each sample was prepared and measured as described in Viéville et al. [19]. The crude plant extract was obtained from the brown algae Sargassum muticum as described in Vonthron-Sénécheau et al. [39]. Two equivalent samples were prepared by dissolving each time 18.3 mg in 0.75 mL MeOD plus 0.4 mL D 2 O, chloroquine was added to one sample at a concentration of 0.16 mg mL -1 (1%w/w of plant dry extract). The NMR experiments were run on a Avance III Bruker spectrometer, operating at 700 MHz, and equipped with a TXI cryo-probe. DOSY were acquired with the convection compensated experiment using bipolar pulses [40] (dstebpgp3spr pulse program). 50 gradient increments from 0.5 G/cm to 52.5 G/cm were used, with a cosine roll-off PFG shape. Each elementary PFG had a duration of δ = 1.1 ms and the diffusion delay ∆ was set to 150 ms. For each gradient intensity, 64 scans were acquired with a relaxation of 1.5 s, for a total experimental time of 1 hour 50 minutes. Each 1D spectrum was apodised with an unshifted sine-bell, zerofilled once and Fourier transformed. A spline baseline correction was applied, as well as a correction of small shifts caused by possible instabilities of the temperature control [41].

Processing.

The PALMA algorithm was implemented with the programming language python version 2.7 using the numpy/scipy libraries. The algorithm was then packaged in a plugin to the SPIKE program developed in our group [42]. All programs are available from the authors. All computations were performed on a Macintosh Mac Pro dual Xeon with a total of 8 cores, equipped with 32 GB of memory and running MacOsX 10.7.

The DOSY experiments presented in Fig. 2 and Fig. 3 were processed column-wise with λ = 0.01. The noise level along the gradient axis was estimated by the difference of the actual data to the result of a polynomial smoothing, this value being used to estimate η in equation ( 5) (See ESI, section 5.2 † for more details). The diffusion axis was logarithmically sampled from 50 to 10.000 µm 2 s -1 over 256 points. Experiments on PEO were processed with a maximum of 20.000 iterations, for a total time of about 30 seconds for 110 DOSY profiles computed. Experiments on plant extract were processed taking all signals with an estimated SNR above 20 (26 dB) in the first 1D spectrum of the experimental matrix, with a maximum of 200.000 iterations, for a total time of about 1 hour for 1200 to 1400 DOSY profiles. The algorithms ITAMeD, p tailored-ITAMeD, and TRAIn utilised in the ESI were used as downloaded from the respective web-sites, using MATLAB program version R 2013b on MacOs.

Server.

A Web server is available at http://palma.labo.igbmc.fr, where users may submit datasets for automatic processing. The python code of the PALMA algorithm is available on the same server as well as on https://github.com/delsuc/PALMA.

Results

Tests on Simulated Data

The PALMA algorithm described above was first tested on a series of simulated data-sets. Fig. 1a presents the results obtained on the simulated experiment A consisting in the superposition of three monodisperse species, separated by less than a factor of 4 in diffusion coefficients, equivalent to the test used in Urbańczyk et al. [25]. When analysed with PALMA using a null λ value, indicating a pure 1 regularisation, a Laplace spectrum consisting of 3 sharp peaks is produced as expected. Using the pure MaxEnt mode (λ = 1) on the same dataset, a broader spectrum is reconstructed. PALMA allows the weight between the two approaches to be freely varied. When doing so it can be observed a narrowing of the MaxEnt distribution for decreasing value of λ, characteristic of a bias toward monodisperse distributions, with a sudden transition to sharp lines for λ = 0. In Fig. 1b the same procedure was applied on a broad Gaussian line simulating a polydisperse polymer with a PDI of about 6, corresponding to experiment B. Again it can be observed that a null λ gives rise to sharp, sparse lines, this time only sampling the broad line in an inadequate manner. However, in contrast with the previous case, all the non-null values of λ lead to a nearly perfect reconstruction of the line-shape, with a correct determination of its width. In this example, the minimal reconstruction error was obtained for λ = 0.05. Tests performed on asymmetric distributions (see ESI figures S5 and S6 and table S1 †) show the same tendency, with stable results for all non-null values of λ. It should be noted that, because of the algorithm that maintains the analysis within the noise distance of the data, i.e. H H HX -Y η (see ( 5)) all the reconstructed Laplace spectra fit equally well the data. They differ only in how they match the regularisation term, a term which holds and expresses the a priori information we have on the data-set. The PALMA algorithm was tested against ITAMeD [25], p tailored-ITAMeD [26], and TRAIn [28] algorithms, using the same simulated data as above. In our hands, PALMA and TRAIn present the best results in terms of faithfulness and robustness, in particular for polydisperse data-sets, with TRAIN showing better results on the C2 data-set, while PALMA behaving better for data sampled with a small number of data-points.

Application to polydisperse polymers.

It is acknowledged that experimental data are quite different from simulated one, with a mixture of sharp and large diffusion distribution, tainted with instrumental artifacts and nonstationary noise. To the behaviour of the method on polydisperse systems, the program was first applied on DOSY experiments measured from poly-ethyleneoxide (PEO) polymers in water with calibrated chain lengths and polydispersity. Fig. 2 presents the results obtained for three PEO samples with polydispersity ranging from 1.07 to 2.0 measured in a standard manner, and processed with PALMA. The polydispersity of the different samples can be clearly seen in the profile widthes. Sample a) is a standard PEO polymer, with a rather low polydispersity. Sample c) was prepared from a set of rather monodisperse polymers, in order to cover regularly a large range of chain lengths. Sample b) on the other hand is a PEO polymer given to have a standard polydispersity, however the details of the composition is not known.

Application to plant extract.

To test the robustness of the approach, it was applied to crude ethanolic plant extract obtained from brown algae. Fig. 3 shows the aliphatic region of the DOSY experiment performed on this algae extract. Only the more abundant species are visible at this plot level. This kind of analysis on complex mixture has been extensively used to analyse natural products, plant extract [43], or even adulteration of herbal and dietary supplements [44]. Here the presence of a particular molecule or family of molecule is characterised by the alignment at the same diffusion coefficient of the characteristic lines located at their corresponding chemical shift positions. For instance, in Fig. 3, the polyol and fatty acid signals are outlined. The fatty acids are 3, but plotted at a level four times lower; b) the same brown algae extract with 0.16 mg mL -1 of chloroquine added, showing the aromatic signals and the methyl signal outlined with an arrow in the spectrum, the dashed horizontal label is at 285µm 2 s -1 ; c) the 1D spectrum of the brown algae extract with added chloroquine. In all three spectra, the aromatic panel is plotted four times lower than the aliphatic panel.

certainly partly aggregated in this sample, as indicated by the width of the line, both along the spectral and Laplace axes. The DOSY is useful in this context because it provides a high dynamic analysis, where the most intense lines do not "hide" less intense ones. In Fig. 4, the chloroquine molecule was added at a low concentration (1% w/w) to the same sample. While the aromatic signals of chloroquine, located in a rather empty region of the spectrum, are easily detected, the aliphatic chain signals fall in crowed region where they are difficult to observe. The methyl groups that fall at position 1.3 ppm and 1.45 ppm are completely buried under the fatty acid signals and the diffusion coefficient are not different enough. In contrast, the signals from the methylene moieties are observed around 1.9 ppm, in a relatively free spectral region. Despite being about ×100 smaller than the larger signals (methanol, polyol or fatty acid chains), the signal is well separated in the DOSY spectrum, and is aligned with the aromatic signals. One can also observed the reproducibility of the PALMA-processed DOSY spectra, as the main features of the spectra are nearly identical for both samples.

Discussion

The hybrid regularisation proposed in eq [START_REF] Callaghan | Translational dynamics and magnetic resonance: principles of pulsed gradient spin echo NMR[END_REF] implements two well known approaches, namely the 1 regularisation, which tends to minimise the number of signals required to explain the data, and the MaxEnt regularisation usually presented as a way to maximise the posterior probability of the analysis while preserving the positivity of the retrieved spectra [45]. While both regularisation are well established and based on clear principles, the MaxEnt is known to be somewhat more difficult to implement. The hybrid regularisation proposed here allows to produce very robust results even in the case of complex signals such as the one presented in the simulated examples. The constrained problem is solved using a new convex optimisation algorithm, based on the use of proximity operators and a split version of the Douglas-Rachford procedure. The use of the proximity operators allows to implement a simple incremental step, requiring no inner line-search minimisation step, where the main burden is three applications of the linear operator H H H or of its generalised inverse B B B = (I I I + H H H T H H H) -1 . This algorithm allows a rapid convergence even with the hybrid regularisation used here. In the simulation presented here, an approximate solution is obtained very rapidly (in less than a second). The results presented in this work were obtained with longer convergences, using typically 10.000 to 100.000 iterations, however thanks to the rapidity of the iterative step, this corresponds typically to a few seconds on a laptop. Because of the constrained approach used here, there is no need to determine a Lagrangian parameter as in most other techniques (sometimes called a smoothing parameter ). Nevertheless, the approach requires some parameters, of which the prior a and the noise η can readily be estimated from the experimental data-set, using respectively the first point of the decay and an estimate of the noise from a polynomial smoothing of the experimental data-set (see ESI section for details †). The weight λ between the MaxEnt and the 1 regularisations embodies a prior assumption on the presence of sparse component in the Laplace spectrum. In the simulations of a sparse theoretical spectrum (Fig. 1a) λ = 0 corresponding to a pure 1 regularisation provides the best reconstruction as expected. Simulations performed on several wide distributions (Fig. 1b and ESI figures S4, S6, S15, and S16 †) show that the method recovers faithfully the position of the signal and the theoretical profile, for most non-null values of the λ parameter, even in the presence of noise. Results are more contrasted for the sparse spectrum in Fig. 1a, where the pure MaxEnt analysis (λ = 1.0) presents large features located at the position of the three components, and it can be observed a narrowing of the MaxEnt distribution for decreasing value of λ. The results of pure MaxEnt analysis produced by the PALMA algorithm were checked to be equivalent to the results obtained with the original algorithm [23] based on a fixed point approach, however the convergence is much faster, and processing times are about 10 times shorter for the same results. It should be recalled that the MaxEnt analysis produces a statistical analysis of the data, where the final spectrum is the density distribution which maximises the posterior probability of finding a signal [START_REF] Jaynes | Probability Theory: The Logic of Science[END_REF]45]. The width of the actual signals can thus be considered as an uncertainty on the position of the monodisperse components, uncertainty that is present in the pure 1 case as errors in the position of the lines, but not directly manifest. These results suffer however from a lack of resolutive power, and the possibility to bias toward a sparser result is certainly a plus. In a general approach, the optimal value for λ should be chosen from assumptions on the data based on explicit previous knowledge, however on a practical point, this is not feasible. Even on a sample known to be composed solely of monodisperse species, the choice of a null λ is problematic. On the one hand, some polydisperse impurities might be present with the risk of overlooking them as we observe in Fig. 2 and Fig. 3; on the other hand, most instrument imperfections such as temperature drift, gradient non-linearity, phase distortions, convection, etc. will distort the pure exponential decay and create some apparent polydispersity. Confronted with the same difficulty, Urbańczyk et al. [26] chose to vary the p parameter of their tailored-ITAMeD algorithm, somewhat similar to λ and chose the larger value which allows a minimal residual. The same approach could easily be used here, however considering the fact that polydisperse samples are correctly analysed for most non-null values of λ, with the better results obtained for small values, we suggest using values between 0.01 and 0.05 as monodisperse data are well described with these values. The quality of these results are in sharp contrast with equivalent analyses presented in the literature. The λ = 1 mode reproduces the classical MaxEnt regularisation [23] although in less processing time. The λ = 0 mode, compared to the 1 based ITAMeD approach [25], the final resolution resulting from PALMA reconstruction is much higher, as observed in our simulations, in agreement to what has been published. Intermediate values of λ, creating a bias of the MaxEnt solution toward more sparse data, produce more resolved spectra, which usually better match the patterns expected in solution NMR. We recommend using a value of λ in the range 0.01 to 0.05 for safe results, with the possibility to adapt this value in particular cases (for instance extreme polydispersity or spectral superposition). We do not recommend to use the pure 1 mode (λ = 0), even in the case of monodisperse samples, because the instrumental fluctuations mentioned above certainly disturb this pure behaviour, and because of the difficulty to estimate the uncertainty of the result from a sparse spectrum alone.

Conclusion

The DOSY experiment holds a special position among the different experiments available to the NMR spectroscopist. Whilst it provides invaluable information on the size and the interaction of the molecules in solution, it presents an important challenge for its acquisition and its analysis. he acquired data are usually corrupted by many artefacts produced by the very high sensitivity of this experiment to instrument imperfections such as temperature drift, non-linearity of the gradients and of the detection electronics, phase distortions, convection, etc. In parallel, while the parameter dependency expressed in the basic evolution equations is extremely simple, it is well established that a simple fit of the data to this equation usually fails in providing a faithful analysis of the data in the general case, and one has to resort to inverse Laplace transform for the analysis step, a problem well known to be of extreme noise sensitivity. For these reasons, the stability and robustness of the acquisition and processing schemes are of great importance for the quality of DOSY spectroscopy, and many acquisition schemes and many processing procedures have already been proposed in the literature for this purpose. In this work, we have introduced a general method to solve the inverse problem as found in the analysis of DOSY experiments that we believe provides an unequalled level of quality and robustness in the processing step. The method is based on a constrained regularisation of the least square problem, and we showed that an hybrid regularisation, combining maximum sparsity and maximum entropy provides optimal results in the general case. This approach is controlled by a single parameter λ weighting between these two criteria, and the results are not very sensitive to the exact value of this parameter as long as extreme values are not chosen. This method, that we called PALMA, is faster and more robust than previous MaxEnt implementation and provides better results. It requires only weak assumptions from the user, and can be run in a fully automatic manner. It has been implemented on a server freely available at http://palma.labo.igbmc.fr, where users may submit their data-sets for automatic processing. The code of the algorithm is available at the same address.

Let us first consider the case when

λ ∈]0, 1]. Then, φ(u) = 0 ⇔ u -x n + λ a log(u) + λ a - λ a log(a) + 1 -λ = 0 ⇔ log(u) = a λ -u + x - λ a + λ a log(a) -1 + λ ⇔ u = exp - a λ u + ax -a(1 -λ) λ + log(a) -1 .
Finally:

u = λ a W a λ exp ax -a(1 -λ) λ + log(a) -1 . (6)
In the above expression, W states for the Lambert function, also called Omega, defined as the inverse function of f : x → x exp(x) for all x ∈ C:

z = x exp(x) ⇔ x = W(z).
If λ = 0, we have:

x n -u ∈ ∂|u| ⇐⇒ x n -u ∈ sign(u) if u = 0 [-1, 1] elsewhere ⇐⇒ u =    x n -1 if u > 1 0 if u ∈ [-1, 1] x n + 1 elsewhere.
Finally, we can conclude that (∀X = (x n ) 1≤n≤N ∈ R N ), prox Ψ (X) is given by (4), where, for all n ∈ {1, . . . , N },

p(x n ) = λ a W a λ exp axn-a(1-λ) λ + log(a) -1 if λ ∈]0, 1] sign(x n ) (max(|x n | -1, 0)) if λ = 0. ( 7
)

Asymptotic development

Let us emphasize that the Lambert function has many interesting properties. In particular:

W(e x ) -→ x→+∞ x 1 - log(x) 1 + x (8)
This asymptotic development is of main interest in our context. Indeed, when λ ∈]0, 1], for all n ∈ {1, . . . , N },

p(x n ) = λ a W (exp(c n )) , (9)
with c n = axn-a(1-λ) λ -1 + 2 log(a)log(λ). When c n is large (typically c n > 10 2 ), standard numerical implementation of the Lambert function yields infinity as an output.

We thus propose to use, for large c n , the following approximation of the proximity operator, which is a consequence of the asymptotic development (8):

p(x n ) = λ a c n - c n 1 + c n log(c n ) . (10)
Figure 8 illustrates relative error of the approximation, when λ = a = 1.

PALMA algorithm

Using the PPXA+ algorithm, we propose a new proximal algorithm to solve (1) that makes use of the proximal formula [START_REF] Stilbs | [END_REF]. The so-called PALMA algorithm, standing for "Proximal Algorithm for L 1 combined with MAxent prior", is given below:

Initialization V (0,1) ∈ R N , V (0,2) ∈ R M X (0) = (I N + H H H H H H) -1 (V (0,1) + H H H V (0,2) ) B B B = (I N + H H H H H H) -1 γ ∈ (0, 2), Minimization For k = 0, 1, . . .            Z (k,1) = prox Ψ (V (k,1) ) Z (k,2) = proj •-Y ≤η V (k,2) U (k) = B B B(Z (k,1) + H H H Z (k,2) ) X (k+1) = X (k) + γ(U (k) -X (k) ) V (k+1,1) = V (k,1) + γ 2U (k) -X (k) -Z (k,1) V (k+1,2) = V (k,2) + γ H H H(2U (k) -X (k) ) -Z (k,2)
where I N is the identity matrix of R N and the projection proj •-Y ≤η is defined as follows:

proj •-Y ≤η (Z) = Z + (Z -Y ) min η Z -Y , 1 -Y (∀(Y, Z) ∈ (R N ) 2 ).

Choice of processing parameters

As it is detailed in the previous theoretical section, the PALMA algorithm is controlled by several scaling parameters which have to be adapted to the current problem.

Choice of a

The entropy approach requires an expression of the prior knowledge on the system, which in the present case corresponds to a priori spectrum in the absence of experimental evidences. A general expression is thus a flat spectrum of intensity a, the value of which has to be adapted to the scaling of the experimental measurement. A natural approach consists in estimating the area under the signal expected signal: n x n and scaling the prior with this value. Due to the properties of the Laplace transform, we have chosen here a = y 0 ≈ n x n .

Choice of η

If we assume that the experimental values are tainted by an additive random noise: y m = ŷm + ε m , and if this noise is supposed to be centered and Gaussian i.i.d. with variance σ 2 , then we can expect the residual of the fit to be η ≈ σ √ M . The value of σ can be measured here either from additional measures or from a simple polynomial fit of the y m curve as it is done in the current implementation. The case of a correlated Gaussian noise ε, with covariance matrix Σ, is encompassed by our method. Indeed, the least squares constraint becomes:

(∀X ∈ R N ) (H H HX -Y ) Σ -1 (H H HX -Y ) 1/2 ≤ η,
which is equivalent to H H HX -Y 2 ≤ η, up to any change of variable with the form H H H ←-θΣ -1/2 H H H and Y ←-θΣ -1/2 Y , η ←θη, for some θ > 0. In order to facilitate the choice of η, we propose to take θ = σ max , i.e. the maximal singular value of Σ, so that a suitable choice for η is σ max √ M .

Choice of λ

In order to evaluate the influence of the λ factor on the balance between sparsity and the amount of information in the processed signal, we launch our PALMA algorithm for different signal of monodisperse and polysdisperse species (A, B, C1 and C2) and we trace on Figures S2 to S6, the recovered signal by variying λ from 0 to 1.

The associated reconstruction quality, measured in terms of signal to noise ratio (SNR) = 10 log 10 X X-X of each tested case, is reported in Table S1. Larger value of quality of reconstruction correspondents to the best reconstruction.

• Signal A Signal A consists in three monodisperse components with diffusion coefficients 16 µm 2 /s, 63 µm 2 /s, and 230 µm 2 /s, with respective intensities 1.0, 0.33 and 0.66. This data-set is equivalent to the simulation used in Kazimierczuk et al.

• Signal B Signal B is a wide monodisperse distribution, simulated as a symmetric log-normal distribution centered at 35 µm 2 /s and with variance 25.

• Signals C1 & C2 C1 and C2 signals are asymmetric distributions built from 15 log-normal components, ranging from 18 to 85 µm 2 /s, with intensities ranging from 0.1 to 10 s. They have PDI estimated respectively to 1.79 and 1.32. For C1, the intensities are chosen as 10.0/(1.4 i ), i = 1, . . . , 15, and for C2, as 10.0/(1.4 15-i ), i = 1, . . . , 15. In all simulations, an additive zero-mean white Gaussian noise and standard deviation σ equals 0.1% of the initial point of H H HX was added to the observed data. The number of measurement is M = 64 and the original signal has the dimension N = 256.

Figure 1 :

 1 Figure 1: PALMA Reconstruction for various values of λ of the simulated experiments, with an added 0.1% Gaussian noise. a) experiment A with 3 monodisperse species, indicated by the black vertical lines, λ = 0 spectrum was divided by 3 for clarity; b) experiment B with a large polydisperse Gaussian profile λ = 0 spectrum was divided by 8 for clarity.

Figure 2 :

 2 Figure 2: DOSY profile of the main NMR signal of different standard PEO samples, a) a reference PEO with PDI of 1.07, b) a reference PEO with PDI of 1.28, c) a mixtures of reference PEOs with a global PDI of 2.01,

Figure 3 :

 3 Figure 3: DOSY experiment on a brown algae methanol/water extract showing only the major constituents.Dashed horizontal label lines are indicated for fatty acid chains (220 µm 2 s -1 ), glycerol and short polyol (380 µm 2 s -1 ), methanol (1100 µm 2 s -1 ).

Figure 4 :

 4 Figure 4: Comparison of two DOSY experiments on a brown algae extracts. a) the aromatic (empty) and aliphatic regions of the same experiment as shown in Fig.3, but plotted at a level four times lower; b) the same brown algae extract with 0.16 mg mL -1 of chloroquine added, showing the aromatic signals and the methyl signal outlined with an arrow in the spectrum, the dashed horizontal label is at 285µm 2 s -1 ; c) the 1D spectrum of the brown algae extract with added chloroquine. In all three spectra, the aromatic panel is plotted four times lower than the aliphatic panel.

Figure S1 :

 S1 Figure S1: Relative error value between W(e cn ) and the approximation of p(x n ) in 9. The relative error becomes lower than 10 -3 for c n > 20.

Figure S3 :

 S3 Figure S2: Measurement of A, B, C1 and C2 signals

Figure S4 :

 S4 Figure S4: Reconstruction of signal B using PALMA with different λ. The minimal error is obtained for λ = 0.05 (Spectra obtained with λ = 0 have their intensities divided by 8 for clarity).

Figure S5 :

 S5 Figure S5: Reconstruction of signal C1 using PALMA with different λ. The minimal error is obtained for λ = 0.01 (Spectra obtained with λ = 0 have their intensities divided by 8 for clarity).

Figure S6 :

 S6 Figure S6: Reconstruction of signal C2 using PALMA with different λ. The minimal error is obtained for λ = 0.01 (Spectra obtained with λ = 0 have their intensities divided by 8 for clarity).

Figure S7 :

 S7 Figure S7: Reconstruction of signal B, using different algorithms. Here, M = 64 and noise level = 1%.

Figure S8 :

 S8 Figure S8: Reconstruction of signal B, using different algorithms. Here, M = 64 and noise level = 0.1%.

Figure S9 :

 S9 Figure S9: Reconstruction of signal B, using different algorithms. Here, M = 64 and noise level = 0.01%.

Figure S10 :

 S10 Figure S10: Reconstruction of signal B, using different algorithms. Here, M = 64 and noise level = 0.001%.

Figure S11 :

 S11 Figure S11: Reconstruction of signal C2, using different algorithms. Here, M = 64 and noise level = 1%.

Figure S12 :

 S12 Figure S12: Reconstruction of signal C2, using different algorithms. Here, M = 64 and noise level = 0.1%.

Figure S13 :

 S13 Figure S13: Reconstruction of signal C2, using different algorithms. Here, M = 64 and noise level = 0.01%.

Figure S14 :

 S14 Figure S14: Reconstruction of signal C2, using different algorithms. Here, M = 64 and noise level = 0.001%.

Figure S16 :

 S16 Figure S16: Reconstruction of signal C2 with λ optimal for different noise levels.

Table 1 :

 1 Table 1 presents synthetic results, extensive results are presented in ESI (see figures S7 to S14 and tables S2 and S3 †). Quality of reconstruction of signal B with different algorithms for various noise levels. Quality is computed as

	Algorithm		noise level
		1%	0.1% 0.01% 0.001%
	ITAMeD	3.37 18.65 29.04 29.40
	ITAMeD with p 6.06 25.26 36.69 37.08
	TRAIn	24.75 28.63 26.53 19.47
	PALMA λ = 0.01 20.54 28.57 41.69 53.25
	PALMA λ = 0.05 24.01 32.51 48.28 51.37

Xsim

Xsim-X calc expressed in dB. For each noise level, the highest quality results are outlined.

Table S1 :

 S1 Quality of reconstruction of A, B, C1 and C2 signals for various λ values.6 Comparison with state-of-the-art algorithmsSeveral algorithms have been developed to solve the ill-posed problem in DOSY experience. In TableS2, we present the comparison results between our approach and several recent algorithms, namely ITAMeD, ITAMeD with p , and TRAIN.

		λ	Qlty reconstruction of Signal A B C1 C2
	0.0 0.01 1.05 5.65 -11.04 -9.87 -12.58 28.57 15.21 16.15 0.02 0.80 31.87 13.76 As an illustration, we present in Figures S7 to S10, the reconstruction of B and C2 14.89 signals with different algorithms for 4 different noise levels.
	0.05 0.58 The following settings have been used, which lead to the best performance in terms of 32.51 12.92 14.06
	0.1 both reconstruction quality and computational cost: 0.47 31.11 12.71 0.2 0.39 31.17 12.37 PALMA: λ = 0.01. 0.5 0.30 31.44 11.65 ITAMED : Regularization parameter = 10 -6 . ITAMeD with p : Smoothing parameter τ = 10 -7 , ration between 1st and 2nd term 13.66 13.35 12.91
	= 10. TRAIn: τ = 1.02	1	0.26	31.55	11.21	12.44
	Signal	Algorithm	Noise level	1%	Qlty reconstruction in dB 0.1% 0.01% 0.001%
		ITAMeD			3.37		18.65	29.04	29.40
	B	ITAMeD with p		6.06		25.26	36.69	37.08
		TRAIn			24.75 28.63	26.53	19.47
		PALMA with λ = 0.01	20.54	28.57	41.69	53.25
		PALMA with λ = 0.05	24.01	32.51	48.28 51.37
		ITAMeD			1.57		15.13	14.36	14.06
	C2	ITAMeD with p TRAIn		3.37 8.62		6.30 15.39	6.29 23.36	6.29 25.5
		PALMA with λ = 0.01	10.6	12.72	17.72	23.24
		PALMA with λ = 0.05	7.62		10.97	16.59	20.75

Table S2 :

 S2 Quality of reconstruction of signals B and C2 with different algorithms for various noise levels. Here, M = 64

	M	Algorithm	Noise level	1%	Qlty reconstruction in dB 0.1% 0.01% 0.001%
		ITAMeD		1.57		15.13	14.36	14.06
	64	ITAMeD with p	3.37		6.30	6.29	6.29
		TRAIn		8.62		15.39	23.36	25.5
		PALMA with λ = 0.01	10.6	12.72	17.72	23.24
		PALMA with λ = 0.05	7.62		10.97	16.59	20.75
	32	ITAMeD ITAMeD with p	-3.72 8.89 0.57 6.18	12.89 6.27	13.91 6.27
		TRAIn		3.68		6.60	9.396	18.04
		PALMA with λ = 0.01	11.09	14.60	19.48	23.04
		PALMA with λ = 0.05	8.69		13.09	18.92	20.56

Table S3 :

 S3 Quality of reconstruction of signal C2 for various noise levels and problem sizes.

Acknowledgements

This work was supported by the CNRS MASTODONS project under grant 2016TABASCO, by the Agence Nationale pour la Recherche, grant ANR2014 ONE SHOT 2D FT ICR. The authors thank Jean-Christophe Pesquet for the initial idea and for discussion all over this project. We thank Catherine Vonthron-Sénécheau, Mélanie Bourgeot and Laure Marguerite for the plant extract sample. We thank Samuel Nicaise for the first version of the WEB server interface and Julien Seiler for helping in the deployment of the WEB server.

* madelsuc@unistra.fr

Variational formulation

A well-known and very efficient strategy for solving ill-posed inverse problems is to adopt a penalization approach that provides an estimate X ∈ R N of the original signal X ∈ R N , that is the solution of the constrained minimization problem: minimize

where Ψ is the regularization function and η is an estimate of the experimental noise.

In this work, we choose to define Ψ for every X ∈ R N and λ ∈ [0, 1], as:

where the Shannon entropy with a flat prior a > 0 and the 1 norm are defined, respectively, as:

Proximity operator

To solve Problem (1), we propose to rely on a proximal optimization method, that makes use of the so-called proximity operator. Since Ψ is convex on R N , its proximity operator at a given point X ∈ R N is defined as the unique minimizer of 1 2

• -X 2 + Ψ. Since Ψ in (2) takes a separable form, i.e. it can be written as:

its proximity operator is given by:

Hereabove, for every n ∈ {1, . . . , N }, p(x n ) ∈ R is the unique minimizer of:

Robustness against noise

We propose to evaluate the quality of reconstruction of signal C2 for different λ values and for 4 different noise levels. Figure S15 shows that the best quality reconstruction is obtained when λ = 0.01 whatever the different noise level. For this optimal λ PALMA algorithm ensures a great quality of reconstruction as it is illustrated in Figure S16.


Cherni Afef

Chouzenoux Émilie

Delsuc Marc-André

PROXIMITY OPERATORS FOR A CLASS OF HYBRID SPARSITY + ENTROPY PRIORS. APPLICATION TO DOSY NMR SIGNAL RECONSTRUCTION

Keywords: Inverse problems, Laplace inversion, entropy regularization, sparsity prior, proximity operator, spectroscopy, nuclear magnetic resonance

Inverse problems arising from Laplace transform inversion are ill-posed, and require suitable regularization strategies. Although the maximum entropy regularization approach usually appears as an adequate strategy due to its ability to recover regular positive valued signals, it was observed to lead to poor reconstruction results when the sought signal contains narrow peaks. In that case, a sparsity promoting penalty such as the 1 norm, combined with a positivity constraint, is more suitable. In order to derive a flexible resolution method, hybrid approaches combining both entropy and sparsity regularization strategies should be envisaged. However, the choice of an efficient optimization algorithm remains a challenging task. Among available optimization techniques, proximal methods have shown their efficiency in solving large scale possibly nonsmooth problems. This paper provides an extensive list of new proximity operators for the sum of entropy and sparsity penalties. The applicability of these results is illustrated by means of experiments, in the context of DOSY NMR signal reconstruction.

INTRODUCTION

The fast resolution of large size ill-posed linear inverse problems presents a big challenge in the context of biophysical data processing [START_REF] Pence | Big data and chemical education[END_REF]. A difficult inverse problem, arising for instance in the context of diffusion ordered nuclear magnetic resonance spectroscopy (DOSY NMR) [START_REF] Johnson | Diffusion ordered nuclear magnetic resonance spectroscopy: principles and applications[END_REF], consists in the estimation of an original unknown positive-valued signal x = (x n ) 1≤n≤N ∈ [0, +∞[ N , from degraded measurements y = (y m ) 1≤m≤M ∈ R M related to x through a linear model

y = Kx + w, (1)
where K = (K m,n ) 1≤m≤M,1≤n≤N ∈ R M ×N is the observation matrix corresponding to a discretized version of a Laplace transform, and w = (w m ) 1≤m≤M ∈ R M represents some additive acquisition noise. An efficient resolution strategy is to employ a penalized approach that defines an estimate x ∈ R N of x as a solution to the following constrained minimization problem [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF]:

minimize x∈R N Ψ(x) subject to Kx -y ≤ τ, (2)
where τ > 0 is a parameter depending on the noise characteristics, and Ψ is a regularization function that allows to add prior information on the sought signal and to impose the positivity of its entries. The so-called maximum entropy reconstruction approach [START_REF] Eggermont | Maximum entropy regularization for Fredholm integral equations of the first kind[END_REF], corresponding to the choice of the Shannon entropy function [START_REF] Gray | Entropy and Information Theory[END_REF][START_REF] Basseville | On entropies, divergences, and mean values[END_REF] for the penalization term Ψ, has been at the core of several papers dealing with regularized inverse Laplace transform [START_REF] Laue | Maximum entropy method in nuclear magnetic resonance spectroscopy[END_REF][START_REF] Livesey | Maximum entropy analysis of quasielastic light scattering from colloidal dispersions[END_REF][START_REF] Livesey | Analysing the distribution of decay constants in pulse fluorimetry using the maximum entropy method[END_REF][START_REF] Delsuc | Maximum entropy processing of DOSY NMR spectra[END_REF][START_REF] Mariette | Continuous relaxation time distribution decomposition by MEM[END_REF][START_REF] Daniell | Maximum entropy and bayesian methods: Cambridge, england[END_REF]. A more recent approach consists in adopting for Ψ a criterion enforcing both sparsity and positivity, with the aim to improve the resolution of narrow peaks possibly present in the sought signal [START_REF] Urbanczyk | Iterative thresholding algorithm for multiexponential decay applied to pgse nmr data[END_REF][START_REF] Kaipin | Trust-region algorithm for the inversion of molecular diffusion nmr data[END_REF][START_REF] Worley | Convex accelerated maximum entropy reconstruction[END_REF][START_REF] Berman | Laplace inversion of low-resolution NMR relaxometry data using sparse representation methods[END_REF], but this strategy may be at the price of loosing the smoothness of the solution. Hybrid regularization approaches combining both entropy and sparsity terms in Ψ should thus be envisaged so as to derive a flexible resolution method. However, the choice of an efficient optimization algorithm to solve the resulting constrained problem (2) (or its Lagrangian formulation) remains a challenging task. On the one hand, in the case of entropy regularization, the optimization techniques proposed in the literature are usually ad hoc and difficult to extend to other types of priors [START_REF] Gorodnitsky | An extension of an interior-point method for entropy minimization[END_REF][START_REF] Chouzenoux | Efficient maximum entropy reconstruction of nuclear magnetic resonance t1-t2 spectra[END_REF][START_REF] Willis | Maximum entropy image restoration revisited[END_REF]. On the other hand, in the case of sparsity regularization, the minimization step is usually handled by proximal optimization methods [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. These methods, grounded on the use of the proximity operator [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF] are highly flexible, and particularly efficient for handling functionals for which the proximity operators have a closed form expression. In this paper, we provide a comprehensive list of formula for the proximity operators of new hybrid regularization functions combining a non necessarily convex sparsity enhancing term and an entropy penalty. This list constitutes a very helpful tool for deriving fast resolution methods for the regularized inversion of [START_REF] Pence | Big data and chemical education[END_REF]. The applicability of our method is illustrated by means of an example in the context of DOSY NMR.

The paper is organized as follows: In section 2, we introduce our notation and recall the definition of the proximity operator. Section 3 presents our main contribution, that is the expression of the proximity operators of several combination of entropy and sparsity promoting penalties. Then, Section 4 presents an application of our results to the resolution of an inverse problem arising in DOSY NMR signal reconstruction. Finally, some conclusions are drawn in Section 5.

PROBLEM STATEMENT

Hybrid regularization

Let us define the following family of hybrid regularization functions of the form:

(∀x ∈ R N ) Ψ(x) = αΨ 1 (x) + βΨ 2 (x), (3)
where Ψ 1 : R N →] -∞, +∞] and Ψ 2 : R N → R are lower semi continuous (lsc) and proper functions and (α, β) are some positive weights. In this paper, we focus on the case when Ψ 1 is either the Shannon or the Burg entropy [START_REF] Nityananda | Maximum entropy image reconstruction-a practical non-information-theoretic approach[END_REF] and Ψ 2 is a sparsity promoting prior. An efficient strategy to promote the sparsity of the sought signal is to choose Ψ 2 as the 0 penalty that counts the number of nonzeros in x. However, the later function is neither differentiable nor continuous at 0, so that its convex approximation, the 1 norm, that sums the absolute value of the signal entries is often used instead [START_REF] Candès | Enhancing sparsity by reweighted l1 minimization[END_REF]. Continuous, but non convex, approximations have also been proposed for the 0 penalty, namely the log-sum penalty [START_REF] Wipf | Iterative reweighted 1 and 2 methods for finding sparse solutions[END_REF] and the Cauchy penalty [START_REF] Antoniadis | Wavelet thresholding for some classes of non-Gaussian noise[END_REF]. Note that all the aforementioned regularization terms are separable, so that (3) can be rewritten as

(∀x ∈ R N ) Ψ(x) = N n=1 ψ(x n ), (4)
with ψ = αψ 1 +βψ 2 , ψ 1 : R →]-∞, +∞] and ψ 2 : R → R.

Figures 1 and2 illustrate the considered functions, for the entropy and sparse penalty terms ψ 1 and ψ 2 , respectively.

Proximity Operator

When Ψ is convex on R N , its proximity operator at x ∈ R N is defined as the unique minimizer of Ψ+ 1 2

•-x 2 [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF]. This operator has been generalized for lsc proper functions that are not necessarily convex in [27, Sec.XV-4], as the multi-valued operator: Fig. 2. Examples of sparsity promoting penalties, in the scalar case: 0 (continuous black line), 1 (dashed-dotted green line), log-sum (dashed red line) and Cauchy (dotted blue line) priors.

prox Ψ : x → Argmin y∈R N 1 2 y -x 2 + Ψ(y) . (5)
Since Ψ in (4) takes a separable form, its proximity operator ( 5) is given by ( [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF]):

prox Ψ (x) = (p n (x n )) 1≤n≤N , (6)
where, for every n ∈ {1, . . . , N },

p n (x n ) = prox αψ1+βψ2 (x n ). (7)
In the sequel, we will focus our attention on the resolution of the scalar optimization problem [START_REF] Laue | Maximum entropy method in nuclear magnetic resonance spectroscopy[END_REF], when α ∈]0, +∞[ and β ∈ [0, +∞[.

PROXIMITY OPERATORS FOR HYBRID SPARSE + ENTROPY PRIORS

Shannon entropy

The Shannon entropy regularization is defined, for all x ∈ R, as

ψ 1 (x) =      x log x if x > 0 0 if x = 0 +∞ elsewhere. (8)
According to [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF], the proximity operator of αψ 1 at x ∈ R N , with α ∈]0, +∞[, reads:

prox αψ1 (x) = αW 1 α exp x α -1 , (9)
where W states for the W Lambert function [START_REF] Corless | On the lambert W function[END_REF]. The remainder of this section shows how to generalize this expression in the case when β ∈]0, +∞[.

3.1.1. Shannon entropy + 1 Set (α, β) ∈]0, +∞[ 2 .
The hybrid Shannon entropy + 1 penalty is defined, for every x ∈ R, as

ψ(x) =      αx log x + βx if x > 0 0 if x = 0 +∞ elsewhere. (10)
Its proximity operator reads, for every x ∈ R:

prox ψ (x) = αW 1 α exp x -β α -1 . (11)

Shannon entropy + 0

Set (α, β) ∈]0, +∞[ 2 . The hybrid Shannon entropy + 0 penalty is defined, for every x ∈ R, as

ψ(x) =      αx log x + β if x > 0 0 if x = 0 +∞ elsewhere. ( 12
)
Its proximity operator reads, for every x ∈ R:

prox ψ (x) =      p if β < β {0, p} if β = β 0 elsewhere, (13)
where β = 1 2 p 2 + αp ∈]0, +∞[ and

p = αW 1 α exp x α -1 . ( 14
)
Note that, due to the discontinuity of the 0 penalty in 0, prox ψ is multi-valued in the limit case when β = β.

Shannon + log-sum

Set (α, β, δ) ∈]0, +∞[ 3 . The hybrid Shannon entropy + logsum penalty is defined, for every x ∈ R, as

ψ(x) =      αx log x + β log(δ + x) if x > 0 β log(δ) if x = 0 +∞ elsewhere. ( 15
)
For every x ∈ R, its proximity operator is given by

prox ψ (x) = Argmin p∈]0,+∞[ s.t. ϕ(p)=0 1 2 (x -p) 2 + ψ(p) ,
with ϕ(p) = p 2 +(δ-x+α)p+α(δ+p) log(p)+δ(α-x)+β.

Shannon + Cauchy

Set (α, β, δ) ∈]0, +∞[ 3 . The hybrid Shannon entropy + Cauchy penalty is defined, for every x ∈ R, as

ψ(x) =      αx log x + β log(δ + x 2 ) if x > 0 β log(δ) if x = 0 +∞ elsewhere. ( 16
)
For every x ∈ R, its proximity operator is given by

prox ψ (x) = Argmin p∈]0,+∞[ s.t. ϕ(p)=0 1 2 (x -p) 2 + ψ(p) ,
with ϕ(p) = p 3 +(α-x)p 2 +(δ+2β)p+α(δ+p 2 ) log(p)+δ(α-x).

Case of Burg entropy

The Burg entropy regularization is defined, for all x ∈ R, as

ψ 1 (x) = -log x if x > 0 +∞ elsewhere. ( 17
)
According to [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF], the proximity operator of the Burg entropy function reads, for every x ∈ R N , for every α ∈]0, +∞[:

prox αψ1 (x) = x + √ x 2 + 4α 2 . ( 18
)
Let us now present the extension of this result to the case when β ∈]0, +∞[.

Burg entropy + 1

Set (α, β) ∈]0, +∞[ 2 . The hybrid Burg entropy + 1 penalty is defined, for every x ∈ R, as

ψ(x) = -α log x + βx if x > 0 +∞ elsewhere. ( 19
)
Its proximity operator reads, for every x ∈ R ( [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF]):

prox ψ (x) = x -β + (β -x) 2 + 4α 2 . ( 20
)

Burg + 0

Set (α, β) ∈]0, +∞[ 2 . The hybrid Burg entropy + 0 penalty is defined, for every x ∈ R, as

ψ(x) = -α log x + β if x > 0 +∞ elsewhere. ( 21
)
Its proximity operator reads, for every x ∈ R:

prox ψ (x) = x + √ x 2 + 4α 2 . ( 22
)

Burg + log-sum

Set (α, β, δ) ∈]0, +∞[ 3 . The hybrid Burg entropy +log-sum penalty is defined, for every x ∈ R, as

ψ(x) = -α log x + β log(δ + x) if x > 0 +∞ elsewhere. ( 23
)
For every x ∈ R, its proximity operator is given by

prox ψ (x) = Argmin p∈]0,+∞[ s.t. ϕ(p)=0 1 2 (x -p) 2 + ψ(p) , with ϕ(p) = p 3 + (δ -x)p 2 + p(β -δx -α) -δα.

Burg + Cauchy

Set (α, β, δ) ∈]0, +∞[ 3 . The hybrid Burg entropy + Cauchy penalty is defined, for every x ∈ R, as

ψ(x) = -α log x + β log(δ + x 2 ) if x > 0 +∞ elsewhere. (24)
For every x ∈ R, its proximity operator is given by

prox ψ (x) = Argmin p∈]0,+∞[ s.t. ϕ(p)=0 1 2 (x -p) 2 + ψ(p) , with ϕ(p) = p 4 -xp 3 + (δ + 2β -α)p 2 -δxp -δα.

APPLICATION TO NMR SIGNAL RECONSTRUCTION

Problem formulation

During the DOSY (Diffusion Order SpectroscopY) experiment, proposed by [START_REF] Morris | Diffusion-ordered 2D NMR spectroscopy[END_REF] to analyze the properties of complex chemical mixtures, a series of measurements are acquired for different pulsed field gradient strengths, the data are then processed with the aim to separate different species according to their diffusion coefficient. The DOSY NMR data y ∈ R M gathers the results of M ≥ 1 experiments corresponding to different acquisition settings characterized by the set of parameters t = (t m ) 1≤m≤M . The problem is then to estimate, from these measurements, the values of the true diffusion distribution χ(T ) at given positions T = (T n ) 1≤n≤N . The relation between y and the sought signal x = (χ(T n )) 1≤n≤N ∈ R N can be written under the form (1) where K ∈ R M ×N is given, for every m ∈ {1, . . . , M } and every n ∈ {1, . . . , N } by K m,n = exp(-T n t m ),

where w ∈ R M is a perturbation noise. In practice, the noise is assumed to be zero-mean Gaussian, i.i.d, with known standard deviation σ > 0 so that we propose to find an estimate x ∈ R N of x by solving (2) with τ = η √ M σ, where η > 0 is a weight closed to 1 [START_REF] Trussell | Convergence criteria for iterative restoration methods[END_REF].

Experimental results

We now present the experimental results obtained when, for every n ∈ {1, . . . , N }, for every m ∈ {1, . . . , M },

t m = t min + m-1 M t max , T n = T min exp -(n-1) N log Tmin Tmax , (26)
with M = 50, N = 200, T min = 1, T max = 10 3 , t min = 0, t max = 1.5. We consider two synthetic signals A and B. The synthetic signal A models a monodisperse distribution with symmetric log-normal shape located at T = 16 while signal B corresponds to a polydisperse distribution that is the sum of two log-normal patterns located at T 1 = 4 and T 2 = 32. Table 1 presents the resulting signal to noise ratio (SNR) defined as 10 log 10 x 2 / x -x 2 which was obtained when solving Problem (2) for both datasets, for different values of σ, and several choices for the regularization term Ψ, namely (8), ( 10), ( 17) and [START_REF] Willis | Maximum entropy image restoration revisited[END_REF]. Here, we choose to focus only on convex priors so that the optimization problem can be efficiently solved with the PPXA+ algorithm [START_REF] Pesquet | A parallel inertial proximal optimization method[END_REF]. Parameter η is set equals to 1.2, since it was observed to lead to the best reconstruction results. Moreover, when hybrid penalties are considered, we fix β = 1 -α, and we optimize α ∈]0, 1] manually so as to obtain a minimal reconstruction error. One can observe that the addition of the 1 norm in the penalty functions improves significantly the quality of restoration results, especially in the case of the polydisperse signal B. The best results are obtained with the combination of Shannon entropy and 1 penalties. This is also confirmed by visual inspection, as it can be noticed in Figure 3

CONCLUSION

This work expands the extensive list of proximity operators available in the litterature [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF][START_REF] Gheche | A proximal approach for optimization problems involving kullback divergences[END_REF] by considering the case of separable functions combining entropy and sparsity promoting terms. Through numerical experiments, we show that these functions may serve as efficient hybrid penalties for solving ill-posed inverse problems in the context of DOSY NMR spectroscopy. Since only convex penalties were tested in our experiments, the PPXA+ algorithm was retained. It should however be emphasized that the novel proximity operators we derive could be applied in a variety of proximal algorithms, in the convex [START_REF] Komodakis | Playing with duality: An overview of recent primal-dual approaches for solving largescale optimization problems[END_REF] or the non-convex case [START_REF] Chouzenoux | Variable metric forward-backward algorithm for minimizing the sum of a differentiable function and a convex function[END_REF][START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF]. The latter case should be of particular interest in the context of blind signal restoration problems such as those encountered in [START_REF] Repetti | Euclid in a taxicab: Sparse blind deconvolution with smoothed 1/ 2 regularization[END_REF][START_REF] Toumi | A review of blind source separation in nmr spectroscopy[END_REF] where the proposed hybrid penalties could be beneficial.

Fig. 1 .

 1 Fig. 1. Examples of entropy based penalties, in the scalar case: Shannon (continuous red line) and Burg (dashed blue line) entropy priors. Note that ψ 1 = +∞ for negative values of x.

  displaying an example of reconstruction result we obtained in the case of dataset B, with σ = 10 -5 .

Fig. 3 .

 3 Fig. 3. Reconstruction results for dataset B and σ = 10 -5 .

Table 1 .

 1 SNR in dB of the restored signals A and B for various choices of the penalization function Ψ.

			Shannon Shannon	Burg	Burg
		σ	prior + 1 prior prior + 1 prior
	Dataset A	10 -2 12.45 10 -3 18.16 10 -4 20.87	13.16 20.86 25.95	12.92 12.11 12.03	12.92 13.44 15.53
	Dataset B	10 -3 11.14 10 -4 18.11 10 -5 19.05	18.52 20.23 26.30	6.20 7.41 7.54	9.54 10.95 10.98

This work was supported by the CNRS MASTODONS project under grant 2016TABASCO and by the Agence Nationale pour la Recherche (ANR, France) under grant 2010FT-ICR2D and grant Défi de tous savoirs 2014, ONE-SHOT-FT-ICR-MS-2D).