Graduate Student Member, IEEE, Pol Pakshal Bohra
email: pakshal.bohra@epfl.ch

Member, IEEE Del Agui Pla

Jean-François Giovannelli
email: jean-francois.giovannelli@u-bordeaux.fr

Fellow, IEEE Michael Unser
email: michael.unser@epfl.ch

Del Aguila Pol

Pla

Keywords: Inverse problems, minimum mean-square error, convolutional neural networks, sparse stochastic processes. I

We present a statistical framework to benchmark the performance of reconstruction algorithms for linear inverse problems, in particular, neural-network-based methods that require large quantities of training data. We generate synthetic signals as realizations of sparse stochastic processes, which makes them ideally matched to variational sparsity-promoting techniques. We derive Gibbs sampling schemes to compute the minimum meansquare error estimators for processes with Laplace, Student's t, and Bernoulli-Laplace innovations. These allow our framework to provide quantitative measures of the degree of optimality (in the mean-square-error sense) for any given reconstruction method. We showcase our framework by benchmarking the performance of some well-known variational methods and convolutional neural network architectures that perform direct nonlinear reconstructions in the context of deconvolution and Fourier sampling. Our experimental results support the understanding that, while these neural networks outperform the variational methods and achieve near-optimal results in many settings, their performance deteriorates severely for signals associated with heavy-tailed distributions.

A Statistical Framework to Investigate the Optimality of Signal-Reconstruction Methods prior knowledge about the signal of interest is required for the resolution of such problems.

A. Model-Based Methods

Model-based methods rely on the mathematical modeling of the signal of interest to counteract the ill-posedness of the inverse problem. We organize them in two categories.

The first category is composed of linear reconstruction methods (e.g., filtered back-projection), which are fast, well understood, and come with performance and stability guarantees [START_REF] Tikhonov | Solution of incorrectly formulated problems and the regularization method[END_REF], [START_REF] Bertero | Introduction to Inverse Problems in Imaging[END_REF]. From a variational standpoint, they can be interpreted as minimizers of a cost functional that consists of a quadratic data term to ensure consistency with the measurements, along with an additive quadratic (Tikhonov) regularization term that imposes some smoothness on the solution. Interestingly, these methods can also be derived from a statistical perspective as optimal linear reconstructors under the Gaussian hypothesis [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF].

The second category is composed of methods that exploit sparsity-the property that a signal admits a concise representation in some transform domain (e.g., wavelets) [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF], [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF], [START_REF] Baraniuk | Applications of sparse representation and compressive sensing [scanning the issue[END_REF], [START_REF] Elad | Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing[END_REF]. This powerful concept supports the theory of compressed sensing, which gives conditions under which the reconstruction of an image from a limited set of measurements is feasible [START_REF] Donoho | Compressed sensing[END_REF], [START_REF] Candès | An introduction to compressive sampling[END_REF], [START_REF] Foucart | An invitation to compressive sensing[END_REF] and stable [START_REF] Candes | Stable signal recovery from incomplete and inaccurate measurements[END_REF], [START_REF] Del Aguila Pla | Stability of imagereconstruction algorithms[END_REF]. To obtain a sparse reconstruction, one typically uses 1 -norm regularization and solves the corresponding convex optimization problem using iterative algorithms such as the fast iterative shrinkage-thresholding algorithm (FISTA) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] or the alternating direction method of multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. In practice, sparsity-promoting regularizers such as total variation (TV) [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] generally improve the quality of the image. From a statistical point of view, many of these sparsity-based methods can be interpreted as maximum a posteriori (MAP) estimators for some specific choices of stochastic models for the signal of interest [START_REF] Babacan | Bayesian compressive sensing using Laplace priors[END_REF].

B. Learning-Based Methods

Neural-network-based methods that make use of prior information learned from a large collection of training data are now the focus of much of the current research in image reconstruction [START_REF] Mccann | Convolutional neural networks for inverse problems in imaging: A review[END_REF], [START_REF] Ongie | Deep learning techniques for inverse problems in imaging[END_REF]. They shine in extreme imaging scenarios where one wishes to achieve more with fewer data, for instance when operating with short integration times, which leads to an abundance of noise, or when collecting fewer measurements to reduce either the acquisition duration and/or the radiation exposure [START_REF] Jin | Deep convolutional neural network for inverse problems in imaging[END_REF]. Here, we focus on two classes of neural-network-based methods and classify them as the counterparts of the model-based ones.

The first successful applications of deep convolutional neural networks (CNNs) in imaging build upon the classical linearreconstruction algorithms, training a CNN to correct for reconstruction artifacts in extreme imaging conditions [START_REF] Jin | Deep convolutional neural network for inverse problems in imaging[END_REF], [START_REF] Chen | Low-dose CT with a residual encoder-decoder convolutional neural network[END_REF], [START_REF] Hyun | Deep learning for undersampled MRI reconstruction[END_REF], [START_REF] Monakhova | Learned reconstructions for practical mask-based lensless imaging[END_REF], [START_REF] Perdios | CNN-based image reconstruction method for ultrafast ultrasound imaging[END_REF]. Unrolling methods [START_REF] Gregor | Learning fast approximations of sparse coding[END_REF], [START_REF] Chen | Trainable nonlinear reaction diffusion: A flexible framework for fast and effective image restoration[END_REF], [START_REF] Yang | Deep ADMM-Net for compressive sensing MRI[END_REF], [START_REF] Aggarwal | MoDL: Model-based deep learning architecture for inverse problems[END_REF], [START_REF] Adler | Learned primal-dual reconstruction[END_REF], [START_REF] Monga | Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing[END_REF] also fall into this class of direct, nonlinear reconstructions. Examples of successful applications include MRI, CT, optical imaging, and ultrasound. Their gain over the state-of-the-art is impressive and comparable in magnitude to the one afforded by a decade of refinement of the sparsity-promoting techniques.

The second class includes methods that attempt to reconstruct an image that is consistent with the measurements by replacing the proximal operator that is typically involved in the iterative sparsity-promoting methods by an appropriate denoising CNN, which then plays the role of the regularizer. They come in a variety of flavors, including plug-and-play (PnP) [START_REF] Venkatakrishnan | Plug-and-play priors for model based reconstruction[END_REF], [START_REF] Ryu | Plug-and-play methods provably converge with properly trained denoisers[END_REF], [START_REF] Zhang | Plugand-play image restoration with deep denoiser prior[END_REF], [START_REF] Sun | Scalable plugand-play ADMM with convergence guarantees[END_REF], regularization-by-denoising (RED) [START_REF] Romano | The little engine that could: Regularization by denoising (RED)[END_REF], [START_REF] Sun | Block coordinate regularization by denoising[END_REF], [START_REF] Wu | SIMBA: Scalable inversion in optical tomography using deep denoising priors[END_REF], and projected-gradient-descent [START_REF] Chang | One network to solve them all-Solving linear inverse problems using deep projection models[END_REF], [START_REF] Gupta | CNNbased projected gradient descent for consistent CT image reconstruction[END_REF] methods.

Despite their remarkable performance, CNN-based imaging methods have limitations that currently hinder their further development. Unlike the model-based methods, which are backed by sound mathematics, the development of CNN-based approaches is empirical. Expressivity is obtained through the composition of simple units, but the working of the whole is hard to comprehend and the architectural options are overwhelming (e.g., depth, number of channels, size of the filters). In practice, one usually proceeds by trial and error using the training, validation, and testing errors as quantitative criteria. Further, the training of CNNs is poorly understood and often difficult because of the underlying over-parameterization: getting a stochastic optimization algorithm to perform properly for a specific application typically requires a lot of adjustments and experimentation.

Beside the strain that this empirical approach exerts on developers, the performance greatly depends on the quality, cardinality, and representability of the training dataset, while the outcome is not necessarily transposable to other applications. The bottleneck with biomedical imaging is often a limited access to large, representative datasets. This is mostly because of legal issues in medical imaging and because of the lack of standardized protocols in biomicroscopy. Another issue is the chickenand-egg nature of the training process because the desired image (the physical object that corresponds to the measurements) is not known precisely-in practice, the goldstandard is an image produced by a state-of-the-art model-based method with highdensity/low-noise measurements. This is adequate for developing methods for compressed sensing, but not otherwise. This explains why the works that demonstrate the superiority of the CNN-based approaches over the more traditional model-based methods for image reconstruction have used limited benchmarks so far.

C. Contribution

In this work, we present an objective environment to benchmark the performance of reconstruction algorithms for linear inverse problems. Our proposed framework offers quantitative measures of the degree of optimality (in the mean-square-error sense) for any given reconstruction method. Further, it provides access to large amounts of training data, which enables the benchmarking of CNN-based approaches.

We synthesize ground-truth signals and then simulate the measurement process (e.g., convolution for deconvolution microscopy, Fourier sampling for MRI) in the presence of noise. Specifically, we consider a statistical framework where the underlying signals are realizations of 1D sparse stochastic processes (SSPs) [START_REF] Unser | An Introduction to Sparse Stochastic Processes[END_REF]. The motivation there is that these processes are ideally matched to model-based methods, the most prominent of which can be interpreted as their MAP estimators [START_REF] Kamilov | MMSE estimation of sparse Lévy processes[END_REF]. Since the true statistical distribution of the signal is known exactly in our framework, the minimum-mean-squareerror (MMSE) estimator is indeed optimal in the mean-squareerror (MSE) sense. Therefore, we are able to provide statistical guarantees of optimality by specifying an upper limit on the reconstruction performance.

Our framework also provides training data for CNN-based approaches. Indeed, we can produce any desired number of training pairs for a given reconstruction task and some chosen stochastic signal model, which allows for an informed comparison of network architectures. Thus, the availability of the goldstandard (MMSE estimator) and training data make our benchmark a good ground for the tuning of CNN architectures and for the identification of the best designs in a tightly controlled environment.

The MAP estimates of SSPs are solutions of optimization problems that resemble the ones used in model-based methods, and can be computed efficiently. However, it has been observed that these MAP estimators are suboptimal in the MSE sense [START_REF] Kamilov | MMSE estimation of sparse Lévy processes[END_REF], [START_REF] Amini | Bayesian estimation for continuous-time sparse stochastic processes[END_REF], except in the Gaussian scenario where the MAP and MMSE estimators (generalized Wiener filter) coincide [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]. In this work, we focus on non-Gaussian signal models. In principle, the MMSE estimator involves the calculation of high-dimensional integrals, which are not numerically tractable in general. Thus, we develop efficient Gibbs-sampling-based algorithms to compute the MMSE estimator for specific classes of SSPs, with innovations following the Laplace, Student's t, and Bernoulli-Laplace distributions. To the best of our knowledge, no such working solution for generic linear inverse problems with SSPs has been presented in the literature.

Finally, we present experimental results that illustrate the usefulness of our framework. Specifically, we benchmark the performance of some well-known model-based methods and CNNs that perform direct nonlinear reconstructions, in the context of deconvolution and Fourier sampling for first-order SSPs. The CNNs that we consider are optimized by minimizing the MSE loss for training datasets. On one hand, when the innovations follow a Bernoulli-Laplace distribution, we observe that CNNs (with sufficient capacity and training data) outperform the sparsity-promoting methods, which are well-suited to these piecewise-constant signals. In fact, some of these CNNs achieve near-optimal MSE performance. On the other hand, our experiments with Student's t innovations indicate regimes where CNNs fail to reconstruct the signals well. More specifically, we observe that, when the tails of the Student's t distribution are made heavier (i.e., when we move towards a Cauchy distribution), CNNs perform rather poorly.

D. Roadmap

In Section II, we describe a continuous-domain model for the measurement process along with a way to discretize it. In Section III, we introduce Lévy processes as stochastic models for our signals and we derive the probability distribution for samples of such processes. We then discuss MAP and MMSE estimation in Section IV before we develop Gibbs samplers for Lévy processes associated with Laplace, Student's t, and Bernoulli-Laplace distributions in Section V. Finally, we present experimental results in Section VI.

II. MEASUREMENT MODEL

In the proposed framework, we consider the recovery of a continuous-domain signal s : R → R from a finite number M of measurements y = (y m ) M m=1 .

A. Continuous-Domain Measurement Model

We model the measurements y = (y m ) M m=1 as

y m = R s(t)ν m (t) dt + n[m], (1)
where (ν m ) M m=1 are linear functionals that describe the physics of the acquisition process and n[•] is an additive white Gaussian noise (AWGN) with variance σ 2 n . By choosing appropriate functionals (ν m ) M m=1 , we can study a variety of linear inverse problems such as denoising, deconvolution, inpainting, and Fourier sampling.

B. Discrete Measurement Model

We need to discretize (1) to obtain a computationally feasible model for the measurements. To that end, we consider a finite region of interest Ω = (0, T ) of the signal and approximate it with

s h (t) = K k=1 s(kh)sinc t h -k , ( 2
)
where h is the sampling step and K = ( T h -1). When h is small enough, s h is a good approximation of s within the interval Ω [START_REF] Unser | Sampling-50 Years after Shannon[END_REF]. On introducing (2) into (1), we get that

y = Hs + n, (3)
where s = (s(kh)) K k=1 ∈ R K contains equidistant samples of the signal, H : R K → R M is the discrete system matrix with

[H] m,k = R sinc t h -k ν m (t) dt, (4)
and n ∈ R M is the noise. Thus, for any signal samples s ∈ R K , we can simulate noisy measurements using [START_REF] Bertero | Introduction to Inverse Problems in Imaging[END_REF]. Next, we derive the discrete system matrices for deconvolution and Fourier sampling. Hereafter, we assume for simplicity that h = 1.

C. Deconvolution

In deconvolution, the measurements are acquired by sampling the result of the convolution between the signal and the point-spread function (PSF) ψ of the acquisition system, which we model by letting the measurement functionals be ν m = ψ(m -•). We assume that the cutoff frequency of ψ is ω 0 ≤ π, as this allows us to sample (s * ψ) on an integer grid without aliasing effects. In this case, The entries of the resulting system matrix H are given by

[H] m,k = R sinc(t -k)ψ(m -t) dt = ψ(m -k).
(

) 5
Here, H is a discrete convolution matrix whose entries are samples of the bandlimited PSF ψ.

D. Fourier Sampling

In Fourier sampling, the measurements are acquired by sampling the Fourier transform of the signal at arbitrary frequencies {ω m } M m=1 . Accordingly, the measurement functionals are the complex exponentials

ν m = e -jω m • . Assuming that |ω m | ≤ π, we get that [H] m,k = R sinc(t -k)e -jω m t dt = e -jω m k . ( 6
)
Here, H is a discrete Fourier-like matrix, except that the frequencies ω m do not necessarily lie on an uniform grid.

III. STOCHASTIC SIGNAL MODEL

In this section, we describe a continuous-domain stochastic model for the signal. We also derive the probability distribution for the discrete signal vector s = (s(k)) K k=1 .

A. Lévy Processes

In our framework, the underlying signals are realizations of a well-known class of first-order sparse stochastic processes: the Lévy processes [START_REF] Unser | An Introduction to Sparse Stochastic Processes[END_REF], [START_REF] Ken-Iti | Lévy Processes and Infinitely Divisible Distributions[END_REF].

Definition 1 (Lévy process):

A stochastic process s = {s(t) : t ∈ R + } is a Lévy process if 1) s(0) = 0 almost surely; 2) (independent increments) for any N ∈ N \ {0, 1} and 0 ≤ t 1 < t 2 • • • < t N < ∞, the increments (s(t 2 ) - s(t 1 )), (s(t 3 ) -s(t 2 )), . . . , (s(t N ) -s(t N -1
)) are mutually independent; 3) (stationary increments) for any given step h, the increment process

u h = {s(t) -s(t -h) : t ∈ R + } is stationary; 4) (stochastic continuity) for any > 0 and t ≥ 0 lim h→0 Pr{|s(t + h) -s(t)| > } = 0.
Lévy processes are closely linked to infinitely divisible (id) distributions.

Definition 2 (Infinite divisibility): A random variable X is infinitely divisible if, for any N ∈ N \ {0}, there exist independent and identically distributed (i.i.d.) random variables X

1 , . . . , X N such that X = X 1 + • • • + X N .
For any Lévy process s, the random variable s(t) for some t > 0 is infinitely divisible. Moreover, its probability density function (pdf) is given by

p s(t) (x) = R R p s(1) (y)e jωy dy t e -jωx dω 2π . ( 7
)
Conversely, for any id distribution with pdf p id , it is possible to construct a Lévy process s such that p s(1) = p id . Thus, there is a one-to-one correspondence between Lévy processes and id distributions [START_REF] Ken-Iti | Lévy Processes and Infinitely Divisible Distributions[END_REF]. Among all id distributions, the pdf of the Gaussian distribution exhibits the fastest rate of decay at infinity. In this sense, we refer to the non-Gaussian, heavier-tailed members (e.g., Laplace, Bernoulli-Laplace, Student's t, symmetric-alpha-stable) of the class of id distributions as sparse [START_REF] Amini | Sparsity and infinite divisibility[END_REF]. Indeed, some of these sparse distributions have a mass at the origin in their probability distribution (e.g., Bernoulli-Laplace) and some of them are strongly compressible (e.g., Student's t, symmetric-alphastable) [START_REF] Amini | Compressibility of deterministic and random infinite sequences[END_REF].

The stochastic model of Lévy processes allows us to consider a variety of signals with different types of sparsity. In our framework, we focus on the subclass of Lévy processes associated with the Gaussian, Laplace, Bernoulli-Laplace and Student's t distributions. Some realizations of these processes are shown in Fig. 1.

B. Discrete Stochastic Model

Now, we derive the pdf of the random vector s = (s(k)) K k=1 , which contains uniform samples of a Lévy process. Consider the stationary increment process u(t) = {s(t)s(t -1) : t ∈ R + } whose first-order pdf p u is the same as p s(1) and so is infinitely divisible. Its samples u = (u(k)) K k=1 can be expressed as

u = Ds, (8)
where D is a finite-difference matrix of the form

D = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 0 • • • 0 -1 1 0 • • • 0 0 -1 1 • • • 0 . . . . . . 0 0 • • • -1 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 9
)
Using [START_REF] Elad | Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing[END_REF] and the fact that the increments are independent, we obtain the pdf of the discrete signal as

p s (s) = K k=1 p u ([Ds] k ) . ( 10
)
Note that (8) can also be written as which gives us a direct way to generate samples of Lévy processes.

[s] k = k n=1 [u] n , k = 1, . . . , K, (11)

C. Extensions

In this work, we have considered inverse problems involving 1D signals that are modelled as realizations of Lévy processes with increments that follow the Gaussian, Laplace, Bernoulli-Laplace and Student's t distributions. Our framework can further be extended in a straightforward manner to include the more general signal model of continuous-domain first-order autoregressive processes [START_REF] Unser | An Introduction to Sparse Stochastic Processes[END_REF]Chapter 7] driven by white noises associated with the aforementioned distributions. These AR [START_REF] Unser | Biomedical image reconstruction: From the foundations to deep neural networks[END_REF] processes yield a discrete stochastic model that is similar to the Authorized licensed use limited to: Jean-Francois Giovannelli. Downloaded on November 02,2023 at 16:03:19 UTC from IEEE Xplore. Restrictions apply.

one described in [START_REF] Candès | An introduction to compressive sampling[END_REF]. There, the application of a suitable transformation matrix to the discrete signal vector, which contains equidistant samples of the process, decouples it and generates a random vector (called the innovation or generalized increments) with i.i.d. entries. Thus, the MMSE estimation methods presented in Section V can be readily adapted for such AR(1) processes.

We can also directly extend the proposed framework to handle multidimensional signals for the particular stochastic model of continuous-domain AR Lévy sheets [47, Chapter 3], [START_REF] Unser | An Introduction to Sparse Stochastic Processes[END_REF] associated with the Gaussian, Laplace, Bernoulli-Laplace and Student's t distributions. These are higher-dimensional generalizations (based on separable whitening operators) of the corresponding AR(1) processes and they result in desirable discrete models of the form [START_REF] Candès | An introduction to compressive sampling[END_REF]. Unfortunately, the vectorized discrete signal for other ("non-separable") higher-dimensional stochastic processes described in [START_REF] Unser | An Introduction to Sparse Stochastic Processes[END_REF] cannot be fully decoupled by applying a linear transformation. This makes the task of designing schemes to compute their MMSE estimators very challenging. An alternate way of extending our framework could be to define a new class of continuous-domain multidimensional stochastic models using the spline-operator-based framework of [START_REF] Unser | Ridges, neural networks, and the Radon transform[END_REF]. However, this approach would require substantial development of novel mathematical ideas and is thus not discussed further in this article.

IV. BAYESIAN INFERENCE

So far, we have introduced the signal and measurement models that allow us to generate our ground-truth signals and simulate their noisy measurements for a certain acquisition setup. Next, we focus on statistical estimators for the reconstruction problem at hand, which is to recover the signal s from the measurements y.

In Bayesian inference, the goal is to characterize the posterior distribution p s|y and derive estimators based on it. Using Bayes' rule and (10), we get p s|y (s|y) = p y|s (y|s)p s (s)

R K p y|s (y|s)p s (s) ds ∝ exp - y -Hs 2 2 2σ 2 n K k=1 p u ([Ds] k ) . ( 12
)

A. Maximum a Posteriori Estimator

The MAP estimator calculates the mode of the posterior distribution p s|y and is given by

s MAP (y) = arg max s∈R K p s|y (s|y) = arg min s∈R K 1 2σ 2 n y -Hs 2 2 + K k=1 Φ u ([Ds] k ) , ( 13
)
where Φ u (x) =log(p u (x)). The cost functional in (13) consists of a quadratic data-fidelity term and a penalty term that encodes the prior signal model. The optimization task in [START_REF] Del Aguila Pla | Stability of imagereconstruction algorithms[END_REF] resembles the one formulated in the variational model-based methods. For instance, if p u is a Gaussian pdf, then the penalty term is proportional to Ds 2 2 , which is a classical Tikhonov regularizer [START_REF] Tikhonov | Solution of incorrectly formulated problems and the regularization method[END_REF]. However, if p u is a Laplace pdf, then we have a sparsity-promoting 1 -norm penalty term Ds 1 , which corresponds to the popular TV regularizer [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF].

These MAP estimators can be computed efficiently with the help of iterative algorithms such as gradient descent, FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], and ADMM [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF].

B. Minimum Mean-Square Error Estimator

The MMSE estimator is given by

s MMSE (y) = arg min s∈R K R K s -s 2 2 p s|y (s|y) ds = R K s p s|y (s|y) ds, (14)
which is the mean of the posterior distribution p s|y . For a fixed stochastic model, the MMSE estimator is the optimal reconstructor in the MSE sense and thus serves as the goldstandard in our benchmarking framework.

In the Gaussian case, the MMSE estimator is known to coincide with the MAP estimator and is straightforward to calculate [START_REF] Unser | Biomedical image reconstruction: From the foundations to deep neural networks[END_REF], [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]. However, in the non-Gaussian case, we need to numerically evaluate the high-dimensional integral in [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], which is computationally challenging.

V. MMSE ESTIMATORS FOR SPARSE LÉVY PROCESSES

In this section, we present efficient methods to compute the MMSE estimator for sparse Lévy processes with increments that follow the Laplace, Student's t, and Bernoulli-Laplace distributions, which constitutes a key contribution of this article.

A. Markov Chain Monte Carlo Methods

The MMSE estimator s MMSE involves the calculation of the integral [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. The high dimensionality of this integral makes its approximation by simple techniques such as uniform-grid-based Riemann sums infeasible. Instead, one can use Markov Chain Monte Carlo (MCMC) methods [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF], [START_REF] Geyer | Practical Markov chain Monte Carlo[END_REF], [START_REF] Gilks | Markov chain Monte Carlo in Practice[END_REF], [START_REF] Gamerman | Markov chain Monte Carlo: Stochastic Simulation for Bayesian Inference[END_REF] for the numerical approximation of ( 14) in a tractable manner.

MCMC methods are designed for generating random samples from nontrivial high-dimensional probability distributions. Broadly speaking, the idea in MCMC is to design a Markov chain such that the distribution that one wishes to draw samples from is its stationary distribution. The desired samples can be obtained by simulating the Markov chain and recording its states after convergence.

In order to compute the integral in [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], we first generate samples {s (q) } Q q=1 from p s|y using an MCMC method. We then approximate s MMSE by the empirical mean (q) . Although the collected samples are correlated, the Markov chain central limit theorem [START_REF] Gordin | The central limit theorem for stationary Markov processes[END_REF] guarantees that s Q is a good approximation of s MMSE for a large-enough Q.

s Q = 1 Q Q q=1 s
Algorithm 1: Gibbs Sampling.

1: Input: Q (number of samples), B (burn-in period) 2: Initialization: ( x (0) , y (0) ) 3: for q = 1, . . . , B + Q do 4: Generate x (q) ∼ p x|y (x| y (q-1) ) 5:

Generate y (q) ∼ p y|x (y| x (q) ) 6: end for 7: Output: {(x (q) , y (q) )} Q q=1 = {( x (q+B) , y (q+B) )} Q q=1

B. Gibbs Sampling

In this work, we propose to use the MCMC method called Gibbs sampling [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF], [START_REF] Casella | Explaining the Gibbs sampler[END_REF] to generate samples {u (q) } Q q=1 from the posterior distribution p u|y . These can then be transformed in accordance with [START_REF] Foucart | An invitation to compressive sensing[END_REF] to obtain samples {D -1 u (q) } Q q=1 from p s|y . We now give the gist of this algorithm.

Let x and y be two random variables. Consider the task of generating samples from their joint distribution p x,y . Gibbs sampling is advantageous whenever it is computationally difficult to sample from the joint distribution directly but the conditional distributions p x|y and p y|x are easy to sample from. The steps involved in this method are presented in Algorithm 1. They yield a Markov chain whose stationary distribution is indeed p x,y [START_REF] Casella | Explaining the Gibbs sampler[END_REF]. In practice, one discards some of the initial samples (burn-in period) to allow the chain to converge. Moreover, quantities (expectation integrals) based on the marginal distributions p x and p y can be computed from the individual samples {x (q) } Q q=1 and {y (q) } Q q=1 , respectively. Next, we present Gibbs sampling schemes for Lévy processes with Laplace, Student's t, and Bernoulli-Laplace increments. Our strategy is to introduce an auxiliary vector w and perform Gibbs sampling for the joint distribution p u,w|y [START_REF] Tanner | The calculation of posterior distributions by data augmentation[END_REF], [START_REF] Mira | On the use of auxiliary variables in Markov chain Monte Carlo sampling[END_REF]. The key is to choose w such that the conditional distributions p u|w,y and p w|u,y can be sampled from in an efficient manner.

Hereafter, we assume that the noise variance σ 2 n and the parameters of the signal model are known.

C. Laplace Increments

For Lévy processes with Laplace increments, we adapt the approach that was developed in [START_REF] Park | The Bayesian LASSO[END_REF].

The pdf for the Laplace distribution is

p u (u) = b 2 exp (-b|u|) , ( 15
)
where b is the scale parameter. The density in (15) can be expressed as a scale mixture of normal distributions [START_REF] Andrews | Scale mixtures of normal distributions[END_REF], as

p u (u) = R p u|w (u|w)p w (w) dw, (16)
where

p u|w (u|w) = 1 √ 2πw exp - u 2 2w ( 17
)
is the Gaussian pdf and

p w (w) = b 2 2 exp - b 2 w 2 1 + (w) (18)
is a mixing exponential pdf 1 with λ = 2/b 2 . This property allows us to define an auxiliary random vector w ∈ R K with i.i.d. entries following the distribution p w in [START_REF] Mccann | Convolutional neural networks for inverse problems in imaging: A review[END_REF], such that

p u|w (u|w) = K k=1 p u|w ([u] k |[w] k ) , ( 19
)
where p u|w is shown in [START_REF] Babacan | Bayesian compressive sensing using Laplace priors[END_REF]. Due to the chain rule of probability (or the general product rule), the full joint distribution p y,u,w can be written as p y,u,w (y, u, w) = p y|u,w (y|u, w)p u,w (u, w)

= p y|u (y|u)p u|w (u|w)p w (w). (20)
Consequently, the distribution p u,w|y takes the form p u,w|y (u, w|y) ∝ exp -

1 2σ 2 n y -Au 2 2 × K k=1 [w] -1 2 k exp - [u] 2 k 2[w] k × K k=1 b 2 2 exp - b 2 [w] k 2 1 + ([w] k ) , (21)
where A := HD -1 .

Based on [START_REF] Chen | Low-dose CT with a residual encoder-decoder convolutional neural network[END_REF], the conditional distribution p u|w,y is then obtained as

p u|w,y (u|w, y) ∝ exp - 1 2 1 σ 2 n y -Au 2 2 + u T C L (w)u , ( 22
)
where C L (w) is a diagonal matrix with elements ([w] -1 k ) K k=1 . Specifically, p u|w,y is a multivariate Gaussian pdf with mean

u = σ -2 n (σ -2 n A T A + C L (w)) -1 A T y and covariance matrix R = (σ -2 n A T A + C L (w)) -1 .
There exist several methods for the efficient generation of samples from a multivariate Gaussian density [START_REF] Rue | Fast sampling of Gaussian Markov random fields[END_REF], [START_REF] Orieux | Sampling high-dimensional Gaussian distributions for general linear inverse problems[END_REF], [START_REF] Gilavert | Efficient Gaussian sampling for solving large-scale inverse problems using MCMC[END_REF], [START_REF] Vono | High-dimensional Gaussian sampling: A review and a unifying approach based on a stochastic proximal point algorithm[END_REF].

The conditional distribution p w|u,y is p w|u,y (w|u, y)

∝ K k=1 p w|u,y ([w] k |[u] k , y) , ( 23
)
where p w|u,y (w|u, y) ∝ exp -1 2

u 2 w + b 2 w × w -1 2 1 + (w) (24)
1 The pdf of the exponential distribution is

p exp (x) = (1/λ)e -x/λ 1 + (x),
where λ > 0 is the scale parameter.

Authorized licensed use limited to: Jean-Francois Giovannelli. Downloaded on November 02,2023 at 16:03:19 UTC from IEEE Xplore. Restrictions apply.

belongs to the family of generalized inverse Gaussian distributions 2 with λ 1 = b 2 , λ 2 = u 2 and a = 0.5. We use the method proposed in [START_REF] Devroye | Random variate generation for the generalized inverse Gaussian distribution[END_REF] to draw samples from the pdf in [START_REF] Perdios | CNN-based image reconstruction method for ultrafast ultrasound imaging[END_REF].

To summarize, at each iteration q of the constructed blocked Gibbs sampler, we generate u (q) ∼ p u|w,y (u|w (q-1) , y) and [w (q) ] k ∼ p w|u,y (w|[u (q) ] k , y) for all k ∈ {1, . . . , K}. The collected samples {u (q) } q follow the desired distribution p u|y .

D. Student's t Increments

The case of Student's t increments can be handled by adapting the method shown in [START_REF] Fevotte | A Bayesian approach for blind separation of sparse sources[END_REF], which is in fact similar to the one we described for Laplace increments.

The Student's t pdf is given by

p u (u) = Γ( α+1 2 ) Γ α 2 1 √ π(1 + u 2 ) α+1 2 , ( 25
)
where α is the number of degrees of freedom and controls the tail of the distribution, and where Γ denotes the gamma function.

It can also be expressed as

p u (u) = R p u|w (u|w)p w (w) dw, (26)
where

p u|w (u|w) = w 2π exp - wu 2 2 ( 27
)
is a Gaussian pdf and p w (w) = (0.5)

α 2 Γ( α 2 ) w α 2 -1 exp - w 2 1 + (w) (28)
is the pdf of a gamma 3 distribution. Again, we introduce an auxiliary vector w ∈ R K whose i.i.d. entries follow p w defined in [START_REF] Aggarwal | MoDL: Model-based deep learning architecture for inverse problems[END_REF]. It is such that

p u|w (u|w) = K k=1 p u|w ([u] k |[w] k ) , ( 29
)
where p u|w is defined in [START_REF] Yang | Deep ADMM-Net for compressive sensing MRI[END_REF].

Here, the distribution p u,w|y is given by

p u,w|y (u, w|y) ∝ exp - 1 2σ 2 n y -Au 2 2 × K k=1 [w] 1 2 k exp - [w] k [u] 2 k 2 2
The pdf of the generalized inverse Gaussian distribution is

p gig (x) = (λ 1 /λ 2 ) a/2 2K a ( √ λ 1 λ 2 ) x a-1 e -(λ 1 x+λ 2 /x)/2 1 + (x),
where K a is the modified Bessel function of the second kind, λ 1 > 0, λ 2 > 0, and a ∈ R. 3 The pdf of the gamma distribution is

p gam (x) = 1 λ λ 1 2 Γ(λ 1 ) x λ 1 -1 e -x/λ 2 1 + (x),
where λ 1 > 0 and λ 2 > 0 are the shape and scale parameters, respectively.

× K k=1 [w] α 2 -1 k exp - [w] k 2 1 + ([w] k ) , ( 30
)
where A := HD -1 . Now, the conditional distribution p u|w,y (u|w, y) turns out to be

p u|w,y (u|w, y) ∝ exp - 1 2 1 σ 2 n y -Au 2 2 + u T C T (w)u , ( 31
)
where C T (w) is a diagonal matrix with entries ([w] k ) K k=1 . Similar to the Laplace case, p u|w,y is a multivariate Gaussian density with mean

u = σ -2 n (σ -2 n A T A + C T (w)) -1 A T y and covariance matrix R = (σ -2 n A T A + C T (w)) -1 .
The distribution p w|u,y is again separable and takes the form

p w|u,y (w|u, y) ∝ K k=1 p w|u,y ([w] k |[u] k , y) , ( 32
)
where

p w|u,y (w|u, y) ∝ exp - (1 + u) 2 w 2 × w α-1 2 1 + (w). ( 33
)
is a gamma distribution with λ 1 = α+1 2 and λ 2 = 2 (1+u) 2 , which can easily be sampled from.

E. Bernoulli-Laplace Increments

In [START_REF] Ge | Enhanced sampling schemes for MCMC based blind Bernoulli-Gaussian deconvolution[END_REF], Gibbs sampling schemes have been designed for a deconvolution problem where the underlying signal is an i.i.d. spike train that follows the Bernoulli-Gaussian distribution. Unfortunately, the Bernoulli-Gaussian distribution is not infinitely divisible and so is not compatible with our framework of Lévy processes. While there exists some work [START_REF] Chaari | Sparse Bayesian regularization using Bernoulli-Laplacian priors[END_REF] on Bernoulli-Laplace priors, according to the analysis presented in [START_REF] Ge | Enhanced sampling schemes for MCMC based blind Bernoulli-Gaussian deconvolution[END_REF], their proposed sampler would have a tendency to get stuck in certain configurations. Thus, we build upon the method in [START_REF] Ge | Enhanced sampling schemes for MCMC based blind Bernoulli-Gaussian deconvolution[END_REF] and develop a novel Gibbs sampler for Lévy processes with Bernoulli-Laplace increments.

The Bernoulli-Laplace pdf is

p u (u) = λδ(u) + (1 -λ) b 2 exp (-b|u|), (34)
where λ ∈ (0, 1) denotes the mass probability at the origin and b is a scale parameter. We can represent this same density as

p u (u) = R 1 v=0 p u|v,w (u|v, w)p(v) p(w) dw, (35)
where

p v (v) = (λ) 1-v (1 -λ) v for v ∈ {0, 1} (36)
is a Bernoulli distribution,

p w (w) = b 2 2 exp - b 2 w 2 1 + (w) (37)
is an exponential pdf, and p u|v,w is defined such that

p u|v,w (u|v = 0, w) = δ(u) ( 38
)
p u|v,w (u|v = 1, w) = 1 √ 2πw exp - u 2 2w . ( 39
)
Based on this representation, we introduce two independent auxiliary vectors v ∈ R K and w ∈ R K . Their elements are i.i.d. and follow the distributions p v and p w , as defined in ( 36) and [START_REF] Wu | SIMBA: Scalable inversion in optical tomography using deep denoising priors[END_REF], respectively. Further, these vectors satisfy

p u|v,w (u|v, w) = K k=1 p u|v,w ([u] k |[v] k , [w] k ) , ( 40
)
where p u|v,w is defined in [START_REF] Chang | One network to solve them all-Solving linear inverse problems using deep projection models[END_REF] and [START_REF] Gupta | CNNbased projected gradient descent for consistent CT image reconstruction[END_REF].

Here, the full joint distribution p y,u,v,w is given by p y,u,v,w (y, u, v, w) = p y|u,v,w (y|u, v, w)p u,v,w (u, v, w)

= p y|u (y|u)p u|v,w (u|v, w)

× p v (v)p w (w). (41)
As a result, the distribution p u,v,w|y takes the form

p u,v,w|y (u, v, w|y) ∝ exp - 1 2σ 2 n y -Au 2 2 × K k=1 p u|v,w ([u] k |[v] k , [w] k ) × K k=1 λ 1-[v] k (1 -λ) [v] k × K k=1 b 2 2 exp - b 2 [w] k 2 1 + ([w] k ), (42)
where A = HD -1 .

Let us now introduce some notations. For any binary vector q ∈ R K , let I q,0 and I q,1 denote sets of indices such that [q] k = 0 for k ∈ I q,0 and [q] k = 1 for k ∈ I q,1 . Further, let A(q) be the matrix constructed by taking the columns of A corresponding to the indices in I q,1 . We then define the matrix B(q, r) = σ 2 n I + A(q)C BL (q, r)A(q) T , where r ∈ R K is a vector with positive entries and C BL (q, r) is a diagonal matrix with entries ([r] k ) k∈I q,1 . Here, we also introduce the vector q (-k) ∈ R K-1 that contains all the entries of q except the kth one, so that q (-k) = ([q] 1 , . . . , [q] k-1 , [q] k+1 , . . . , [q] K ) T . Lastly, for q ∈ {0, 1}, we define the vector q q (-k) ∈ R K such that q q (-k) = ([q] 1 , . . . , [q] k-1 , q, [q] k+1 , . . . , [q] K ) T . First, we look at the conditional distribution p u|v,w,y . From [START_REF] Chang | One network to solve them all-Solving linear inverse problems using deep projection models[END_REF] and [START_REF] Amini | Bayesian estimation for continuous-time sparse stochastic processes[END_REF], we deduce that any sample from p u|v,w,y takes the value of zero at the indices in

I v,0 . If we define u 1 = ([u] k ) k∈I v,1 , then we get p u 1 |v,w,y (u 1 |v, w, y) ∝ exp - 1 2 1 σ 2 n y -A(v)u 1 2 2 + u T 1 C BL (v, w)u 1 . (43)
Thus, p u 1 |v,w,y is a multivariate Gaussian density with mean

u 1 = σ -2 n (σ -2 n A(v) T A(v) + C BL (v, w)) -1 A(v) T y and co- variance matrix R = (σ -2 n A(v) T A(v) + C BL (v, w)) -1
. The conditional distribution p w|u,v,y takes the form [START_REF] Ken-Iti | Lévy Processes and Infinitely Divisible Distributions[END_REF] where p w|u,v,y is given by

p w|u,v,y (w|u, v, y) ∝ K k=1 p w|u,v,y ([w] k |[u] k , [v] k , y) ,
p w|u,v,y (w|u, v = 0, y) ∝ b 2 2 exp - b 2 w 2 1 + (w) (45) p w|u,v,y (w|u, v = 1, y) ∝ exp - 1 2 u 2 w + b 2 w × w -1 2 1 + (w). (46)
The densities in ( 45) and ( 46) correspond to the exponential distribution with λ = 2/b 2 and the generalized inverse Gaussian distribution with λ 1 = b 2 , λ 2 = u 2 , and a = 0.5.

Next, inspired by the work in [START_REF] Ge | Enhanced sampling schemes for MCMC based blind Bernoulli-Gaussian deconvolution[END_REF], we consider sampling from the marginalized conditional distribution of [v] k in a sequential manner as this can allow for a more efficient exploration of configurations of v. More specifically, at each iteration q, we draw [v (q) ] k from the distribution p [v] k |v (-k) ,w,y (v|v (q) (-k) , w (q) , u (q-1) ), where v

(q) (-k) = ([v (q) ] 1 , . . . , [v (q) ] k-1 , [v (q-1) ] k+1 , . . . , [v (q-1) ] K ) and k ∈ {1, . . . , K}.
The marginalized posterior distribution p v,w|y is given by p v,w|y (v, w|y) ∝ p y|v,w (y|v, w)p v (v)p w (w),

where

p y|v,w (y|v, w) = R K
p y|u,v,w (y|u, v, w)p u|v,w (u|v, w) du.

(48) It can be shown that ( 47) and ( 48) lead to

p v,w|y (v, w|y) ∝ |B(v, w)| -1 2 exp - 1 2 y T B(v, w) -1 y × K k=1 λ 1-[v] k (1 -λ) [v] k × K k=1 b 2 2 exp - b 2 [w] k 2 1 + ([w] k ). (49)
From ( 49), we see that p = 1 + exp -

1 2 h 1 -v; v (-k) , w, y -h v; v (-k) , w, y -1 , ( 50
)
where

h v; v (-k) , w, y = y T B v v (-k) , w -1 y + log B(v v (-k) , w) + 2v log λ 1 -λ . ( 51
)
To summarize, in each iteration q of the above-described sampler, we generate w (q) ∼ p w|u,v,y (w|u (q-1) , v (q-1) , y), (q) , u (q-1) ) for all k and u (q) ∼ p u|v,w,y (u|v (q) , w (q) , y). This particular order of updates is important as it yields a partially collapsed Gibbs sampler [START_REF] Van Dyk | Partially collapsed Gibbs samplers: Theory and methods[END_REF] where the stationary distribution is still p u,v,w|y .

[v (q) ] k ∼ p [v] k |v (-k) ,w,y (v|v (q) (-k) , w

VI. EXPERIMENTAL RESULTS

In our experiments, we benchmark the performance of some popular signal reconstruction schemes, including a CNN-based method, on deconvolution and Fourier sampling problems with Lévy processes associated with the Bernoulli-Laplace and Student's t distributions.

A. Signal Models

We consider a signal vector s ∈ R 100 that contains samples of a Lévy process whose increments follow the Bernoulli-Laplace or Student's t distribution.

1) Bernoulli-Laplace Increments: The Bernoulli-Laplace pdf [START_REF] Sun | Scalable plugand-play ADMM with convergence guarantees[END_REF] is characterized by the parameters λ and b, where λ determines the mass probability at the origin and b represents the scale of the Laplace component. We perform experiments for models corresponding to λ ∈ {0.6, 0.7, 0.8, 0.9}. The scale parameter is set to b = 1 for each case.

2) Student't t Increments: The Student's t pdf ( 25) is parameterized by α, which controls the tails of the distribution. We conduct experiments for α ∈ {1, 3, 5, 39}.

B. Measurement Models

We consider both deconvolution and Fourier sampling problems for each of the above-described signal models.

1) Deconvolution: As shown in Section II-C, the system matrix H for deconvolution is a discrete convolution matrix. Accordingly, we construct H : R 100 → R 88 such that

H = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ [h] 13 • • • [h] 1 0 • • • 0 0 . . . . . . . . . . . . . . . . . . . . . . . . 0 0 • • • 0 [h] 13 • • • [h] 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 52
)
where h ∈ R 13 consists of the central samples of a truncated Gaussian PSF with variance σ 2 0 = 4.

2) Fourier Sampling: For Fourier sampling in 1D, which is reminiscent of MRI, the forward model H resembles a discrete Fourier matrix (see Section II-D). Thus, in order to construct H, we first sample M = 16 rows of the DFT matrix. The first row of the DFT matrix (DC component) is always kept, while the remaining ones are selected in a quasi-random fashion with a denser sampling at low frequencies. We then create the real system matrix H : R 100 → R M , where M = 2M -1, by separating the real and imaginary parts.

In both measurement models, the AWGN variance σ 2 n is chosen such that the (average) measurement SNR is around 30 dB.

C. Reconstruction Algorithms

For each combination of the signal and measurement models, we compare the performance of some (variational) model-based techniques, a CNN-based scheme and the MMSE estimator. We generate validation and test datasets, each consisting of 1,000 pairs of ground-truth signals and their noisy measurements. Further, in order to train the CNNs, we also synthesize a repository R containing a large number of training examples.

1) Model-Based methods: We consider the model-based methods

s 2 = arg min s∈R K y -Hs 2 2 + τ Ds 2 2 , (53)
s 1 = arg min s∈R K y -Hs 2 2 + τ Ds 1 , (54)
and

s log = arg min s∈R K y -Hs 2 2 + τ K k=1 log 1 + ([Ds] k ) 2 , ( 55
) where τ ∈ R + . ( 53), ( 54) and (55) resemble the MAP estimators of Lévy processes associated with Gaussian, Laplace, and Student's t distributions, respectively. However, unlike the MAP estimators, these include an adjustable hyperparameter τ . For each of these methods, the same regularization parameter τ is used for the entire test dataset. This particular value of τ is the one that yields the lowest MSE for the validation dataset.

These estimators are implemented in MATLAB using Global-BioIm [START_REF] Soubies | Pocket guide to solve inverse problems with Global-BioIm[END_REF]-a library for solving inverse problems. Specifically, the 2 estimator is expressed in closed-form as

s 2 = H T H + τ D T D -1 H T y. (56)
The 1 and log estimators are computed iteratively using ADMM. Since the cost functional in ( 55) is non-convex, we initialize ADMM for s log with s 1 so that it can reach a better local minimum.

2) CNN-Based Method:

The concept here is to train a CNN as a regressor that maps an initial low-quality reconstruction s 0 to a high-quality one s CNN [START_REF] Jin | Deep convolutional neural network for inverse problems in imaging[END_REF], [START_REF] Chen | Low-dose CT with a residual encoder-decoder convolutional neural network[END_REF], [START_REF] Hyun | Deep learning for undersampled MRI reconstruction[END_REF], [START_REF] Monakhova | Learned reconstructions for practical mask-based lensless imaging[END_REF], [START_REF] Perdios | CNN-based image reconstruction method for ultrafast ultrasound imaging[END_REF]. The architecture of the CNN used in our experiments is based on the well-known denoising network DnCNN [START_REF] Zhang | Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising[END_REF] and is described in Fig. 2 andTable I.

First, we build a training dataset of cardinality M T by taking the first M T examples {s m , y m } M T m=1 from the repository R. We

TABLE I CONVOLUTION LAYERS

then train the model by minimizing the MSE loss function

L(θ) = 1 M T M T m=1 s m -s CNN (θ; s 0 (y m )) 2 2 , (57)
where θ represents the learnable parameters of the network, with the help of the ADAM optimizer [START_REF] Kingma | ADAM: A. method for stochastic optimization[END_REF]. The CNN is trained for 1,000 epochs with a batch size of 256 and a weight decay of γ.

The initial learning rate is set as 10 -2 . For some duration of the training (first 600 epochs for deconvolution and first 750 epochs for Fourier sampling), it is decreased by a factor of 0.5 every 50 epochs. We choose the initial low-quality reconstruction to be s 0 (y) = H T y for the deconvolution problems. In the case of Fourier sampling, s 0 (y) is the zero-filled reconstruction. All the CNN-based reconstruction schemes are implemented in PyTorch.

3) Goldstandard (MMSE Estimator): Our MMSE estimators are implemented in MATLAB, according to the methods detailed in Section V. There, we set the number of samples as Q = 8,000 and the burn-in period as B = 3,000 for signals with Bernoulli-Laplace increments. We use Q = 15,000 and B = 5,000 for signals associated with the Student's t distribution.

D. Results

We present our results for all the test datasets in Figs. 3,4, 5 and 6. For the sake of clarity, instead of the MSE, we display the "MSE optimality gap" which is the difference between the MSE obtained by a specific method and the MSE attained by the MMSE estimator. In these figures, the CNNs are labelled using the tuple (F, C, L, M T , γ), where F is the filter size, C is the number of channels, L is the number of layers, M T is the cardinality of the training dataset and γ is the weight decay. For the interested reader, we also provide information about the computation times required by all the methods in the supplementary material.

1) Lévy Processes with Bernoulli-Laplace Increments: Here, we summarize our observations for both the deconvolution and Fourier sampling experiments (Figs. 3 and4). The sparsitypromoting 1 estimator, which corresponds to the popular TV regularization, is known to be well-suited to piecewise-constant Lévy processes with Bernoulli-Laplace increments. As the value of λ increases, these signals become sparser and exhibit fewer jumps. Consequently, we observe that the 1 estimator performs better than the 2 estimator. The log estimator also promotes sparse solutions [START_REF] Wipf | Iterative reweighted 1 and 2 methods for finding sparse solutions[END_REF] and we see that it performs well for these piecewise-constant signals. However, despite the good fit, there is still some gap between the MSE attained by the 1 and log, and MMSE estimators.

The performance of the CNN-based method improves as we increase the capacity of the CNN and the amount of training data. In fact, with sufficient capacity and training data, they outperform the 1 and log estimators. Remarkably, some of the CNNs achieve a near-optimal MSE.

2) Lévy Processes With Student's t Increments: The parameter α allows us to consider a wide range of signals. As α → ∞, we approach the Gaussian regime. The other extreme is α = 1, which corresponds to the super heavy-tailed (sparse) Cauchy distribution. This scenario can be particularly challenging for the correct setting of algorithm parameters. Due to the heavy tails of the Cauchy distribution, the validation and test datasets may contain signals with a vastly different range of values. Consequently, for a given model-based method, the regularization parameter τ that is chosen to yield the lowest MSE for the validation dataset may differ significantly from the value τ * that achieves the lowest MSE on the test dataset. Thus, in Figs. 5 and6, we also include the performance of model-based methods when their regularization parameter is tuned for optimal MSE performance on the test dataset directly. These "boosted" model-based methods are labelled as * 2 , * 1 and log * . In Figs. 5 and6, we can see that the 2 estimator is optimal for a large value of α. As the value of α decreases, the performance of the 2 estimator deteriorates and becomes worse than that of the 1 estimator. For all the cases, the log estimator, which corresponds to a tunable MAP estimator for the Student's t distribution, attains reasonable MSE values. Note that for the deconvolution experiment involving Cauchy signals, there is a significant gap between the MSE values obtained by the 2 and 1 and * 2 and * 1 estimators, respectively. Interestingly, the log estimator is less affected by this issue.

Finally, for both deconvolution and Fourier sampling problems, CNNs with sufficient capacity and training data perform well up to α = 3, after which there seems to be a steep transition and their performance drops sharply. In fact, for Cauchy signals, we observe that the training process for these CNNs is quite unstable-the training loss marginally decreases and seems to converge, and the networks do not generalize to the validation (or test) datasets. We believe that this last example poses an open challenge for designing robust neural-network-based schemes that can handle signals following (super) heavy-tailed distributions.

Fig. 1 .

 1 Fig. 1. Realizations of different Lévy processes as characterized by the corresponding infinitely divisible pdfs.

  [v] k |v (-k) ,w,y is a Bernoulli distribution withp [v] k |v (-k) ,w,y (v|v (-k) , w, y)Authorized licensed use limited to: Jean-Francois Giovannelli. Downloaded on November 02,2023 at 16:03:19 UTC from IEEE Xplore. Restrictions apply.

Fig. 2 .

 2 Fig. 2. Architecture of the CNN, where BN denotes the operation of batch normalization.

Fig. 3 .

 3 Fig. 3. Deconvolution of Lévy processes with Bernoulli-Laplace increments.

Fig. 4 .

 4 Fig. 4. Fourier sampling of Lévy processes with Bernoulli-Laplace increments.

Fig. 5 .

 5 Fig. 5. Deconvolution of Lévy processes with Student's t increments. The figure at the bottom is a zoomed-in version of the dotted rectangular box shown in the figure at the top.

Fig. 6 .

 6 Fig. 6. Fourier sampling of Lévy processes with Student's t increments. The figure at the bottom is a zoomed-in version of the dotted rectangular box shown in the figure at the top.

Authorized licensed use limited to: Jean-Francois Giovannelli. Downloaded on November 02,2023 at 16:03:19 UTC from IEEE Xplore. Restrictions apply.

VII. CONCLUSION

We have introduced a controlled environment, based on sparse stochastic processes (SSPs), for the objective benchmarking of reconstruction algorithms, including CNN-based methods that require lots of training data, in the context of linear inverse problems. We have developed efficient posterior sampling schemes to compute the minimum-mean-square-error estimators for specific classes of SSPs. These yield the upper limit on reconstruction performance and allow us to provide a measure of statistical optimality. We have highlighted the abilities of our framework by benchmarking some popular variational methods and convolutional neural-network (CNN) architectures for deconvolution and Fourier-sampling problems. In particular, we have observed that, while CNNs outperform the variational methods and achieve a near-optimal performance in terms of mean-square error for a wide range of conditions, they can sometimes fail too, especially for signals with heavy-tailed innovations.

ACKNOWLEDGEMENT

We acknowledge access to the facilities and expertise of the CIBM Center for Biomedical Imaging, a Swiss research center of excellence founded and supported by Lausanne University Hospital (CHUV), University of Lausanne (UNIL), École polytechnique fédérale de Lausanne (EPFL), University of Geneva (UNIGE), and Geneva University Hospitals (HUG).

manuscript and approving it for publication was Dr. Xiao Fu. This work was supported in part by the Swiss National Science Foundation under Grant 200020_184646 / 1 and in part by the European Research Council (ERC Project FunLearn) under Grant 101020573. (


Nelson Fonseca

Cyril Mangenot

Sami Hebib

Hervé Aubert

Olivier Pascal

Interleaved multi-band pyramidal antennas combining radio navigation and telemetry satellite applications

published or not. The documents may come


Ashwin Raman
email: ashwin.raman@my.utsa.edu

Hrishitva Patel
email: hrishitva.patel@my.utsa.edu

Web Application for Visualizing Bitcoin Prices Over the Years

Keywords: Python, Flask, Chart.js, Google App Engine

This project involves the development of a web application to visualize historical Bitcoin prices using Python, Flask, and Chart.js. The application retrieves data from the CoinGecko API and presents it in a graphical format on a web page. The web page was then deployed on Google App Engine

I. INTRODUCTION

Cryptocurrency markets, especially Bitcoin, have garnered significant attention in recent years. Analyzing historical price trends is crucial for understanding market behavior and making informed investment decisions. This project aims to create a user-friendly web application for visualizing Bitcoin prices over a specified time period.

II. SYSTEM ARCHITECTURE

A. Backend

The Flask backend handles data retrieval from the CoinGecko API and serves as the interface between the API and the frontend.

1. Data Retrieval: Utilizes the requests library in python to interact with the CoinGecko API in order to get the data points for the graph. We also implemented a function (get_bitcoin_data()) to fetch Bitcoin price data over a specified time period in our case it was over a span of 10 years.

2. Web Server: Used Flask to create a simple web server and implemented a route to render the HTML template with Bitcoin Price data over the 10 year span.

B. Frontend (HTML and Chart.js)

The frontend is responsible for rendering the user interface and displaying the Bitcoin price data in a line graph.

HTML Template:

Created an HTML template (index.html) for the webpage structure. Within this HTML template we integrated the Chart.js library for creating an interactive line chart.

JavaScript (Chart.js):

Configured the Chart.js library to visualize the Bitcoin price data. We also went about and defined the structure of the line chart, including labels, data points, and styling of the graph to our required specifications.

III. IMPLEMENTATION

In order to deploy the website, the following implementation was made. We first set up a Flask project using Visual Studio Code as that was the IDE we had the most experience deploying projects in. Next was to import and install the necessary dependencies used for this project which include Flask and Requests. The next step was to define and function to fetch the bitcoin price data from the CoinGecko API and configure Flask to create a simple web server and established a route for rendering the HTML template. On the front-end side of the Identify applicable funding agency here. If none, delete this text box.

website, we created an HTML template labeled (index.html) and structured it to accommodate the Chart.js library. We then proceded to integrate the Chart.js into the template in order to visualize the Bitcoin price data given by the API allowing for the users to interact with the graph and see the exact price on that given timeframe. During the development of this project, we subjected the code to local deployment in order to see what else should be modified for the frontend and backend of the application. This also allowed us to make sure that we have proper functionality of the webpage and address any issues that arise from the testing phase.

A. Some Mistakes

 One issue that we ran into while going through the testing phase was properly setting up the index.html file as it wasn't being read into the python code  Another issue that was a big problem was when deploying the application onto Google App Engine there was a 502 bad gateway error that we could not figure out how to circumvent.

 A final mistake that was made during this project was not deploying the code properly and not having the proper implementations on my computer in order to deploy it over to the Google App Engine.

IV. RESULTS

The final results were a successful deployment of the web application showcasing the prices of Bitcoin over the span of 10 years. https://bitcoin-website-405714.uc.r.appspot.com/ Here is a link to the deployed webpage. In this section as well we have provided the source code for this project and each section is the file name. This project showcases the feasibility of developing a web application to visualize Bitcoin prices, catering to the increasing demand for accessible and insightful cryptocurrency tools. The integration of external APIs, such as CoinGecko, enhances the application's capabilities by providing real-time data and ensuring its relevance in the everevolving cryptocurrency landscape.

A. App

VI. FUTURE WORK

While the current implementation meets the project's objectives, several avenues for future enhancement exist:

User Authentication: Implement user authentication mechanisms to personalize user experiences and provide secure access to additional features.

Advanced Analysis Tools: Introduce additional charting options and technical indicators, allowing users to perform advanced analyses on Bitcoin price data. Improved User Interface: Enhance the overall user interface for improved user experience, incorporating user feedback and design principles. VII. PANOPTO LINK FOR DEMO https://utsa.hosted.panopto.com/Panopto/Pages/Viewer.aspx?i d=832b7d56-f67e-49d2-a3bb-b0ca001b1ce3