Denoising: from classical methods to deep CNNs
Jean-Eric Campagne

Comparison of adaptive mesh refinement techniques for numerical weather prediction.
    Isidora Jankovcor4
    April 29, 2024
===================================================================================






This paper aims to explore the evolution of image denoising in a pedagological way. We briefly review classical methods such as Fourier analysis and wavelet bases, highlighting the challenges they faced until the emergence of neural networks, notably the U-Net, in the 2010s. The remarkable performance of these networks has been demonstrated in studies such as <cit.>. They exhibit adaptability to various image types, including those with fixed regularity, facial images, and bedroom scenes, achieving optimal results and biased towards geometry-adaptive harmonic basis. The introduction of score diffusion has played a crucial role in image generation. In this context, denoising becomes essential as it facilitates the estimation of probability density scores. We discuss the prerequisites for genuine learning of probability densities, offering insights that extend from mathematical research to the implications of universal structures.


Keywords. denoising, Fourier analysis, wavelet analysis, diffusion models, U-Net




§ INTRODUCTION


Noise reduction is a crucial step in signal processing, aiming to restore a signal of interest disturbed by unwanted interferences. Noise perception is often subjective, as what constitutes a relevant signal for some may be considered noise for others, and vice versa. Thus, the primary objective of noise reduction is to recover the useful signal while minimizing the impact of stochastic disturbances.
Recent advancements in generative models, particularly those based on score matching <cit.>, have sparked renewed interest in noise reduction, opening new perspectives in this field.

In this pedagogical article, while it's impossible to cover all aspects of noise reduction, we focus on a noise reduction method that exploits deep neural networks to reveal what the network learns. For comprehensive reviews, one can refer, for example, to <cit.>. However, to contextualize these advancements, we also examine traditional noise reduction approaches based on classic orthonormal bases. We highlight the differences between Fourier transform-based and wavelet-based analysis, as well as the limitations of each approach. Furthermore, we explore in detail how deep neural networks, such as  <cit.> and  <cit.>, push the boundaries of traditional methods. We will see how they achieve this and what they learn. To support our claims, we rely on previous works, notably that of <cit.>.

For the sake of brevity, we limit our study to the case of additive white Gaussian noise, leaving aside issues related to observation acquisition.



§ NOISE REDUCTION IN "CLASSICAL" ORTHONORMAL BASES



To illustrate, let x(u) denote an original signal with real values, where u is the low-dimensional variable indexing the sample: for example, u=i could be a "time" index in the case of an audio frame, or u=(i,j) a double index to locate a pixel in an image. We assume that x(u) is not affected by quantization effects. To this signal, we add noise z(u), which we consider to be uncorrelated Gaussian white noise with variance σ^2, so the observations y(u) can be written as

    y(u) = x(u) + z(u)    z∼𝒩(0,σ^2)

The task of a denoiser is to restore the original signal x from the observations y. Even if σ is a priori unknown, we will see how to obtain an estimate with the data. To simplify, in many cases, we assume that this estimation is sufficiently relevant.

To begin, let's consider the 1D case depicted in Figure <ref>, where the original signal and observations {y(i)}_i≤ N (N=2028 with Δ t=1 sec.) are superimposed. This type of signal is chosen to represent transient phenomena. Indeed, the vast majority of signals (music pieces, speeches, photos of landscapes, indoor scenes, or faces) are structured and exhibit transitions between regular intervals.



 §.§ Fourier Analysis

A first idea that naturally arises is to consider the Fourier spectra of such signals, as presented in Figure <ref> for those of the noisy signal |ŷ(ω)| (ω=2π f) and of the noise alone |ẑ(ω)| obtained by synthesizing the noise from an estimate of the variance obtained through preprocessing[We will see later for the wavelet coefficient thresholding. In the case of the 1D signal, by taking the constant frame from the beginning, we use σ = 1.4826×MAD(y(u)-med(y(u))) where MAD is the Median Absolute Deviation. One can also use an iterative method that allows determining an initial estimation of x(u), which is subtracted from y(u) to apply the same formula.]. Note that with Δ t=1 sec, the maximum accessible frequency is 0.5 Hz, or ω_max=π, according to the Nyquist theorem.

We observe that the spectral power of the noise approximately dominates for frequencies beyond f_cut=40 mHz or ω_cut=0.245. If we depict the decomposition on the orthonormal Fourier basis as

    y(k) = ∑_ω=0^ω_maxŷ(ω) e^iω k

then we can establish a "Hard" thresholding strategy on the Fourier coefficients as follows:

    ρ_ω_cut(ŷ,ω) = ŷ(ω)    ω≤ω_cut
    
    0    ω>ω_cut
 
which allows obtaining the low-frequency estimation of the original signal by performing the following summation:

    x̃(k) = ∑_ωρ_ω_cut(ŷ,ω) e^iω k = ∑_ω<ω_cutŷ(ω) e^iω k

The result of the thresholding in Eq. <ref> is shown in Figure <ref>. Note in passing that this frequency thresholding mechanically translates, so to speak, into a thresholding on the Fourier coefficients. This is inherent to the nature of the Fourier transform; we will see the difference with wavelet analysis.

What do we observe? Firstly, we have indeed reduced the noise when compared to Figure <ref>, but upon closer inspection, there are artifacts in the form of small oscillations present everywhere, with more pronounced trends at the transition levels. This is where we perceive the flaws and limitations of the Fourier analysis approach.

To begin with, let's recall a result concerning the relationship between the decay of the Fourier spectrum and the regularity of a function:

If g (1D) is bounded and there exist K and ε>0 such that 

    |ĝ(ω)| ≤K/1+|ω|^n+1+ε
 
then g∈ C^n, meaning that g is continuously differentiable n times.

This has three consequences: (1) we can "read" from the power spectrum the degree of regularity of the function, and the faster the decay, the more regular the function, (2) a function whose power spectrum has compact support ("bounded") is infinitely differentiable, (3) the decay of |ĝ(ω)| is governed by the worst regularity of the function over its entire domain. This latter property helps us understand the shape of the power spectrum of the observed signal (Fig. <ref>): since there are discontinuities, the spectrum is constrained to decay only as 1/|ω|, which is the shape of the Fourier spectrum of a "step-function." This constraint masks the fact that the function is regular outside the transition phases[Note a technical point: notice that the chosen signal returns to the same level (0) at the beginning and end of the frame. If this were not the case, due to a simple FFT considering the periodic spectrum with a period equal to its length, we would then have a discontinuity that would generate a spectrum in 1/|ω|. To overcome this technical defect, one would need to use a Discrete Cosine Transform (DCT) that symmetrizes the signal by making it even before considering it periodic (this time with a period twice its length).].

Therefore, as the observed spectrum decays slowly, when placing a cutoff to reduce white noise whose spectrum is constant, ω_cut (f_cut) must be relatively low frequency, depriving us of components resulting from transients. Moreover, by truncating the Fourier series, the residue is an oscillating function composed of high-frequency contributions, hence the oscillations. This phenomenon, which occurs due to thresholding, regardless of whether there is noise or not, is what we call the Gibbs phenomenon.

Now, one might consider that the employed thresholding is very abrupt, and we can gain insight into what happens by taking a filtering perspective. Indeed, let h be a real filter with spectral response ĥ(ω); its usage proceeds with the following operations:

    x̂̃̂(ω)    = ŷ(ω) ĥ(ω)    ⇔  x̃(k) = ∑_ωŷ(ω) ĥ(ω) e^iω k = ∑_ωρ_h(ŷ,ω) e^iω k

where ρ defines a function generalizing the operation <ref> by defining the associated filter as that of a square function used by Cl. Shannon. However, the operation of multiplying spectral responses in Eq. <ref> translates to a convolution in the real domain:

    x̃(k) = (y∗ h)(k)= (h∗ y)(k) = ∫ h(k-u) y(u)du

However, in the case of the square function, its counterpart in the real domain is the sinc function (i.e., the function sin(x)/x), which is notoriously oscillatory and has long tails. We understand the effects of signal averaging that impairs signal reconstruction, especially in transition zones.

However, one could consider changing the filter h (or a function ρ_h) with a smoother spectral response beyond the cutoff frequency while still aiming to reduce the noise portion. There are many filters in the literature as they are widely used in electronics, notably Stephen Butterworth's eponymous filter (1930). Figure <ref> presents the spectral response of this type of filter and its denoising action.



We observe that the result is slightly improved compared to Figure <ref>, but it is clear that the improvement is marginal.
Certainly, other types of thresholding can be applied <cit.>, but the result is generally the same: we would need a tool that can analyze the function locally and concentrate the description into a sparse representation where the largest coefficients reveal the transients. One might consider performing what is called a windowed Fourier transform, but the major drawback lies in defining these windows. In fact, the tool should adapt automatically to the studied function. This tool has existed since the 1990s: it is the wavelet multi-resolution analysis discussed in the following section. The concepts developed for Fourier analysis prove useful in designing such analyses.




 §.§ Wavelet Analysis





  §.§.§ Some Basic Notions


Wavelet theory is elaborated in detail in the reference book by <cit.>. Here are some elements within the scope of this article.

A wavelet ψ(x) is an oscillating function such that

    ∫_-∞^∞ψ(x) dx = 0

thus, in the Fourier domain, it is a band-pass filter. This function undergoes translations and dyadic dilations/contractions[In the case of audio signal processing, the base scale factor is rather 2^1/Q where Q acts like the notion of an octave.] such that, up to a normalization factor:

    ψ_n,j(x) ∝ψ(2^j x - n),    j,n∈ℤ

As j increases, the "size" of the support of the wavelet ψ_n,j(x) contracts, and vice versa for j<0. This allows for a local analysis of a function f(x) by considering the inner products:

    W_j,n(f) = ⟨ψ_n,j, f ⟩  = ∫_-∞^∞ f(v)ψ_n,j^∗(v)  dv

If the variations of f are slow compared to the support of ψ_n,j(x), then W_j,n(f) is zero. The difference between a sine wave and a wavelet for the analysis of signals with transients is clearly seen (Fig. <ref>): the wavelet provides a local introspection of the function unlike the sine wave.

In the Fourier domain, when considering dilations/contractions, we have the opposite effect compared to that in the real space:

    ψ̂_n,j(ω) = 1/2^jψ̂(ω/2^j)

so for j>0, the spectrum ψ_n,j shifts towards high frequencies, while for j<0, it is the opposite.
Unlike a sine wave at a fixed frequency, which is delocalized in real space, a wavelet is localized in both spaces: the wavelet's support in the "time-frequency" plane is low-frequency and well-localized in frequency and quite wide in time, and conversely, at high frequency, the localization in time is good and less so in frequency. The total "area" of the support in the time-frequency plane satisfies the Heisenberg inequality (Fig. <ref>).



To determine whether W_j,n(f) is significant or not, one can reason in the Fourier domain: slow variations of f imply that its spectrum is concentrated at low frequencies, and small scales of ψ_n,j imply that the support of its Fourier spectrum is rather concentrated at high frequencies, so there is no overlap, hence the result W_j,n(f)≈ 0. On the contrary, where there are rapid variations of f during a transient, non-zero coefficients W_j,n(f) can be found. Through the translations at a fixed scale j, the number of significant coefficients is typically given by the width of the wavelet. Hence the sought-after sparse representation at small scales that detect transients.



In fact, instead of considering the wavelets ψ_n,j(x) with j≤ 0, to cover the entire Fourier domain allowing for a complete representation of a function f (Littlewood-Paley condition)[Note: In the Fourier section, we saw that ω_max=π when considering a sampling frequency of 1, so the support of ψ̂(ω) must cover this high part of the spectrum.], one can use a low-pass filter denoted ϕ̂(ω). We only subject ϕ(x) to the translation operation, to obtain inner products:

    V_n(f) = ⟨ϕ_n, f ⟩

An example of the complete collection of filters is given in Figure <ref>.

By inverse transformation, the coefficients V_n(f) provide a large-scale approximation (low frequency) of f, while the coefficients W_j,n(f) provide detailed information about the function.
Thus, the complete characterization of a function f is done through the set of scalar products:

    C(f) = {V_n(f), W_j,n(f)}_j∈ℕ,n∈ℤ


In the specific case of multi-resolution analyses developed by S. Mallat, there is a relationship between ϕ̂ and ψ̂, allowing the design of families of orthonormal bases of L^2(ℝ) (the space of square-integrable functions)[The first orthonormal basis of L^2(ℝ) is that of A. Haar (1909), which is part of the Daubechies wavelet family "db1". D. Gabor and J. Morlet designed in 1948 and 1984 (following the "sniffer planes" affair) what would later be recognized as complex wavelets for their respective research, the first orthonormal basis of L^2(ℝ) and C^∞ being that of Y. Meyer (1986), thus initiating multi-resolution analysis and the development of fast transformation algorithms, including the work of S. Mallat (1989).]. Thus, I. Daubechies (1992) developed wavelets with compact support noted "". Figure <ref> shows an example of ϕ and ψ ("db2"). The design of these wavelets is presented in <cit.>. The small spikes of the "db2" functions are not artifacts. The higher the n order of the wavelet, the more regular it is, as shown by the following theorem:


(Daubechies)

If ψ has p zero moments[The moment of order k is the result of ∫ x^k ψ(x) dx.], then its support is at least equal to 2p-1, and the wavelets of type "" have a minimum support of [-p+1,p], and the corresponding support of ϕ is [0,2p-1].

If ϕ is also constrained to have zero moments, then we obtain Daubechies wavelets called "Coiflets"[The name is a tribute to R. Coiffman who asked I. Daubechies to find wavelets for quadratures.]. Other versions exist in libraries as  <cit.> and  <cit.>.



These wavelets with compact support are very practical and allow the implementation of a fast wavelet transformation algorithm for a finite-size signal (and the inverse reconstruction) in O(N) time, therefore faster than the FFT. Starting from the N samples of the signal, during a transformation, we obtain N coefficients[nb. except for adaptations at the edges] where half of them constitute the low-frequency approximation coefficients, and the other half contains the detail coefficients (high-frequency). The process can then be iterated where at each step, the signal is replaced by the low-frequency coefficients from the previous step. For example, in three steps, an introspection of a signal is performed with filters {ϕ̂(2^jω),ψ̂(2^jω)} with j=0,1,2 where only the wavelet coefficients obtained with the lowest frequency filter ϕ̂(4ω) and all detail filters ψ̂(2^jω) are retained, as shown in Figure <ref>. The schematization of these wavelet transformation steps (DWT for Discrete Wavelet Transform) and its inverse is given in Figure <ref>.



  §.§.§ Application to 1D Denoising




On the right part of Figure <ref> at the top, next to the graph of the noiseless signal x (Fig. <ref>), the coefficients W_j,n(x) (the "details") concerning the largest scale at the bottom and the smallest scale at the top obtained by a fast wavelet transformation[nb. In fact, a "stationary" or "translation-invariant" version of the transformation was used ("undecimated wavelet transform" or "à trous algorithm").] using the "db2" wavelets are shown. It can be observed that the significant coefficients are few (sparse) and concentrate as expected around the discontinuities of the signal.

Since the wavelet transformation is linear, the coefficients of the observed (noisy) signal y decompose as the sum of the coefficients of the noiseless (original) signal and the coefficients of the noise:

    W_j,n(y) = W_j,n(x) + W_j,n(z)

However, since the noise z comes from a random process (Gaussian white noise), the same applies to the coefficients W_j,n(z). The "noisy" nature of W_j,n(y) is observed in the middle part of Figure <ref> on the right. However, it is also noticeable that the significant coefficients of the noiseless signal are hardly affected. One may then consider thresholding of the coefficients defined by[nb. This is a "hard" thresholding in the literature and libraries.] 

    ρ_T(W) = 
    1     |W|≥ T 
    
    0     |W|< T


If we denote by e_ℓ the elements of the orthonormal wavelet basis and C_ℓ(f) the wavelet coefficients of a function f (nb. ℓ is an index combining the pairs (j,n) from the representation Eq. <ref>):

    f = ∑_ℓ C_ℓ(f)  e_ℓ

then the estimation of the signal x from the observed signal y is performed as follows:

    x̃ = ∑_ℓρ_T(C_ℓ(y)) e_ℓ

Notice the similarity to expression <ref> used in the context of Fourier analysis. However, there is a major difference, as the sum above can be written in the form:

    x̃ = ∑_ℓ∈ℰ_T C_ℓ(y) e_ℓ

with

    ℰ_T = {ℓ such that   |C_ℓ(y)|≥ T }

The difference lies in the fact that the membership condition to ℰ_T is conditioned by C_ℓ(y), which depends on the (observed) signal y, it is not the form ℓ<ℓ_cut of equation <ref> used in the Fourier framework. The condition adapts automatically to the shape of y, notably to its transient phases. This is what was sought with sparsity to overcome the shortcomings identified by Fourier analysis.

The choice to clarify is that of the threshold T. It is owed to D. L. Donoho and I. M. Johnstone (1990s) (Theorem <ref> Sec. <ref>) who found that for a large class of problems

    T = σ×√(2log N)        (Donoho & Johnstone)

achieves nearly optimal risk[Indeed, one can find a slightly smaller value that performs better.] and moreover, the wavelet coefficients of the noise are such that

    lim_N→∞ℙ(
    T - σloglog N/log N≤max_0≤ℓ<N |C_ℓ(z)| ≤ T
    ) = 1

The fact that T increases with N reflects the statistics of Gaussian white noise, which makes more and more excursions into the tails of the Gaussian as N increases. Moreover, the value of σ can be obtained by using the Median Absolute Deviation calculated on the detail coefficients of the highest frequency (smallest scale) and adjusted by a coefficient that links it to the variance of a Gaussian[The coefficient is equal to 1.4826, it is the inverse of the quantile of the ℕ(0,1) distribution equal to 3/4.].

The result of such thresholding for the case of the noisy signal is shown in Figure <ref> (bottom left plot). A noticeable improvement can be observed compared to Figure <ref> obtained by Fourier analysis. To provide a quantitative expression for this visual impression, I have reported the signal-to-noise ratio (SNR[nb. its definition is 
snr(x,x_ref) = 10 log_10(x_ref^2/x_ref-x^2)
with x_ref being the noise-free signal]). In this case, if the noisy signal has an SNR of approximately 19 dB, the signal denoised by Fourier analysis achieves an SNR of 22 dB, while by wavelet analysis, an SNR of 39 dB is achieved. Therefore, wavelet analysis presents significant advantages over Fourier analysis in the case of signals containing transients. Let's see how this applies to images.



  §.§.§ Application to 2D denoising


The wavelet theory described briefly in the previous section for the 1D case adapts quite easily to the multi-dimensional case <cit.>. Essentially, one needs to have a low-pass filter and a series of band-pass filters that together cover the necessary Fourier domain to decompose a function on an orthonormal basis.

In the 2D case, one can construct a representation from the 1D described functions ϕ and ψ, defining 3 wavelets that cover 3 Fourier zones as follows:

    ψ^1(u_1,u_2)    = ψ(u_1) ϕ(u_2)
    ψ^2(u_1,u_2)    = ϕ(u_1) ψ(u_2) 
    ψ^3(u_1,u_2)    = ψ(u_1) ψ(u_2)

As this needs to be done at all scales, these three wavelets are dilated/contracted and translated according to

    ψ^k_j,n(u) = 1/2^jψ^k(u-2^j n/2^j)      u=(u_1,u_2), n =(n_1,n_2)

The three types of wavelets detect horizontal contours, i.e., rapidly varying along the vertical axis hence vertical high frequencies, vertical contours, i.e., rapidly varying along the horizontal axis hence horizontal high frequencies, and finally contours that involve both types of variations. The low-pass filter

    ϕ(u_1,u_2) = ϕ(u_1) ϕ(u_2)

allows obtaining an averaged approximation. It just needs to be translated following the same principle as for the 1D case to cover the support of the function to be processed. An example is given with the Haar wavelets ("db1") in Figure <ref>.


The same type of progressive decomposition of an image into low-frequency approximation and detail coefficients is identical to the case described in the section dealing with 1D, using Daubechies wavelets, for example. An example illustrating this is shown in Figure <ref>. The detection of contours marking transients, as in 1D, can be observed. The gray areas indicate nearly zero coefficients, while significant amplitude coefficients are represented by white or black pixels depending on their sign. As observed, we indeed have a sparse representation. Now, as in 1D, progressively, the low-frequency approximation at step j of size N_j× N_j is decomposed into 4 sub-images of sizes N_j/2× N_j/2 with a new, even more averaged approximation and 3 sub-images of "details". The reconstruction is done by adding the "details" obtained at different scales to the obtained low-frequency image.


Now let's consider image denoising. The process indicated for the 1D case is applied as-is in the 2D case. An example using the original image from Figure <ref> is shown in Figure <ref>. The initial image (512× 512) on the left has pixel values in the interval [0,1]. The middle image is the result of adding Gaussian white noise with a moderate σ of 0.1, resulting in a PSNR[PSNR: Peak Signal to Noise Ratio defined as 10log_10(I_max^2/x-x_ref^2) with I_max=1 in our case.] of 20 dB. The image on the right shows the result of denoising using the thresholding method[nb. other thresholding functions may give similar or worse results.] "hard" as described in Section <ref>, taking into account that in the threshold formula, N is replaced by N^2, i.e., the total number of pixels. The wavelet used is the Haar wavelet, and the transformation was pushed to level 3. The PSNR of the denoised image is 28 dB.

 


Therefore, we have effectively denoised the image; however, the same zoomed portion of the images highlights that denoising "blurs" the contours and reduces the contrasts. This is where the limitation of 2D wavelet transformation with wavelets like "db" is pointed out. To study this limitation, we will use a certain type of images in the next section that control the regularity of contours and background.



  §.§.§ Limit of classical wavelet analysis: C^α images



As observed during the denoising of a natural image, there is a phenomenon of blurring contours. To investigate this phenomenon, during the 2000s, "synthetic" images were generated with a background featuring a certain level of regularity controlled by a parameter α, akin to a typical contour within the image. The higher the value of α, the greater the regularity observed in both the background and the contour. These are called C^α images. Examples are shown in Figure <ref>.
 

The denoising of a C^2 image (i.e., α=2) is shown in Figure <ref>. It can be noted that the PSNR is much better than that obtained for a natural image of the same dimension. The evolution of the PSNR of the denoised image is depicted as a function of either the number of pixels or the regularity factor α in Figure <ref>. The log N scaling law, where N represents the number of pixels along one dimension of the image, can be elucidated as follows: when the size is increased by a factor of × 2, the image can be filtered using, for example, a Gaussian kernel of size 2× 2, resulting in the original image size but with the noise variance reduced by a factor of 4. As PSNR is given by -log_10(σ), this indeed yields a log_10 N scaling law. However, determining the exact law is more intricate, as is discerning the shape of the curves illustrating the denoised PSNR versus the noisy PSNR across the parameter α. We will revisit this topic in the section discussing denoising via deep neural networks.



Now, even though denoising is at work, it remains observed as for the natural image, a deterioration of the boundary between the two parts of the C^α image as can be seen in the zooms of Figure <ref> (second row) and the difference between the original image and the denoised image. In fact, as shown by <cit.>, this fundamentally stems from the fact that the support of the wavelets defined by equations <ref>, <ref> does not adapt to the geometry of the image, especially that of the contours. A strategy developed in the 1995-2015 period roughly consisted of finding "directional" wavelets and establishing adapted thresholding algorithms.



  §.§.§ Directional wavelets



As briefly illustrated in the previous section, the classical 2D Wavelet Transform (Eq. <ref>) does not adapt well to the geometry of images, especially to contours. An improvement would be to design directional analyses. The subject is very vast, so I will only take a few cases to illustrate.

One possible solution is to use the Morlet/Gabor wavelet transposed in 2D:

    ψ(u_1,u_2)   =(-e^-ξ^2/2+e^i ξ  u_1) e^-u_1^2+u_2^2/2 σ^2
    ψ(ω_1,ω_2)    =  πσ ^2 e^-1/2σ^2 ω_2^2(e^-1/2σ^2 (ω_1-ξ)^2-e^-1/2(ξ^2+σ^2 w_1^2))

The parameter σ plays a role in the characteristic width of the wavelet, while ξ sets the position of the power maximum in Fourier. An illustration is given in Figure <ref>.
  






Now, one can change the scale λ as we have already done previously (similarly for translation) which we take dyadic (λ=2^j with j∈ℤ), but we can also apply a rotation by an angle θ. It follows[nb. r_θ=[  cosθ -sinθ;  sinθ  cosθ ].]

    ψ_j,θ(u) = 2^-2jψ(2^-jr_θ.u) ⟹ψ_j,θ(ω) = ψ(2^jr_-θ.ω)


  
Then, by combining a collection of wavelets with different parameters j and θ, we can cover the Fourier plane. Figure <ref> illustrates this point. Of course, one can adapt the base wavelet and the number of rotations to fill in any holes left here just to make the different associated "blobs" visible. However, as in 1D, to cover the low frequencies (the central zone, ω∼(0,0)), we resort to using a low-pass filter ϕ(ω) such as a simple Gaussian.  

So finally, the image (signal) is passed through a set of bandpass filters and a low-pass filter (θ_k=π k/K)

    Wx= ((ψ_j,θ_k∗ X)_j∈{j_min,…,j_max},k∈{0,…,K-1}, ϕ∗ X)

which allows the signal to be analyzed at all scales and orientations. It is also necessary to take into account translations in real space in order to have a complete view of the transients throughout the entire image. Note that a particular orientation of a wavelet is sensitive to transients along the perpendicular axis.

From an image, using different orientations and spatial translations, it is "split" into all these filtering channels, and we obtain sub-images that can be represented as in Figure <ref>. It can be noticed that at each scale, the orientation of the wavelets makes orthogonal transitions visible in the original image. Thus, the boundaries between areas of equal intensity are detected, which are the transient phenomena sought.


This scale/orientation decomposition is implemented in the so-called scattering networks architecture <cit.> and their developments as in references <cit.>. For applications in cosmology, see for example <cit.>.


  §.§.§ Wavelets aiming to adapt to geometry


In the same vein as the previous orientable wavelets, orientable wavelets of the steerable wavelet type were developed in the 1990s (see, for example, <cit.>). This marked the second generation of wavelets. An example is given in Figure <ref> for a single decomposition scale. Given an input image, it is first divided into a high-frequency part and a low-frequency part, then band-pass filters are applied to the low-frequency part of the image to detect edges at different angles. Regarding denoising, it is necessary to estimate the noise of the different coefficients in order to apply effective thresholding and restore a denoised image (see, for example, the  method by <cit.>).


Another type of orientable wavelets uses the concept of ridgelets <cit.>, defined as follows: if ψ(x) is a "traditional" wavelet, then the family of functions

    r_s,b,θ(x_1,x_2) = s^-1/2ψ((x_1cosθ + x_2 sinθ - b)/s)

are constant along the line x_1cosθ + x_2 sinθ = constant. A version that contributed to the development of edge detection, compression, and denoising algorithms was introduced later <cit.>. The idea is as follows. In the classical case to obtain the wavelet coefficient (at a certain scale, without translation) of a function f, convolution is performed as follows:

    ∫ f(x_1,x_2) 2^j ψ(2^j x_1,2^j x_2) dx_1dx_2 = ∫ 2^-j f_2^-j(x_1,x_2) ψ(x_1,x_2) dx_1dx_2

with

    f_s(x,y) = f(s x,s  y)

thus there is identical scaling on both coordinates. However, if we consider the function f(x,y) = 1_y≥ x^2, the scaling transformation f(s  x, s^2 y) seems most appropriate (parabolic scaling). In fact, this type of scaling transformation needs to be performed locally in the "tangent-normal" frame. This is where the curvelets are developed, an example of which is shown in Figure <ref> in both the real space and the Fourier space. Fast transformation algorithms have emerged and are described in the articles by <cit.> (e.g.,  [<http://www.curvelet.org>], and  [<https://fr.mathworks.com/matlabcentral/fileexchange/112215-act>]). In astro-cosmology, there are numerous applications such as <cit.>.




  §.§.§ Generalization of Curvelets and Bandlets


A generalization of curvelets was implemented by <cit.> to adapt to contour geometries and exploit the regularity of the image (or function in general) along contours. In terms of shape regularity, the classical transformation exploits possible isotropy regularity over square domains at different scales. Curvelets exploit regularity along quadratic contours, while bandlets generalize the concept <cit.>. An example is shown in Figure <ref>. Many applications have been developed, including image compression <cit.> and denoising, as in the articles by <cit.> (French) and <cit.>. The article by <cit.> provides an application in astronomy, and <cit.> in the field of medical imaging.



Regarding the denoising of C^2 images, what about it? Figure <ref> presents results in addition to those reported in Figure <ref>, obtained with curvelets and steerable wavelets, as well as the thresholding methods used in the codes[ (MATLAB code) and  (<https://www.ipol.im/pub/art/2014/86/>)]. I couldn't find a code implementing bandlets which would have achieved an optimum (see Sec. <ref> Eq. <ref>).



It is clear that while the PSNR of denoised images is better for directional wavelets, the "border zone" is still the least well denoised. This was the state of the art around 2015. Neural networks will change the game, as they did for image classification, then extended to image processing in general, as well as sound processing, speech, and then language-related tasks, etc.


§ DENOISING WITH DEEP NEURAL NETWORKS





 §.§ Architecture: Some Reminders

I won't delve into the entire zoo of neural networks, nor their history. The reader can consult the following survey by <cit.> even not up-to-date for a basic overview of some Deep  neural networks concepts. It's worth noting that denoising with neural networks <cit.> emerged shortly before the advent of AlexNet <cit.>. Networks used for denoising, for example, include the  architecture by <cit.> and [Note: This network was originally designed for medical image segmentation. Its multi-scale structure is also relevant for denoising.] by <cit.> (Fig. <ref>).


These two networks are composed of layers with specific actions:

    
  * Convolutional layer: If we denote (H_in, W_in, C_in) as the dimensions of the input tensor and (H_out, W_out, C_out) for the output tensor, and consider a batch of convolutional kernels 𝒦 of size K × K, then the values of the output tensor are calculated as follows[To simplify, I haven't included a possible bias b_k.]:
    
    T^out_i,j,k = ∑_c=1^C_in∑_(u,v)∈Δ_K𝒦_k,c,u+s,v+s× T^in_i+u,j+v,c  with   s=⌊ K/2 ⌋,  Δ_K={-s,…,s}^2

    One can observe at first glance that not only does this operation act in the two spatial dimensions but also mixes the input channels. It's an operation to be visualized in 3D (Fig. <ref>), where the part concerning the dimensions (H,W) resembles a wavelet decomposition in every aspect except that the filters are learned. However, it's important to note that the operation on the channels leads to the remarkable properties of convolutional networks. If this operation is removed, the performance degrades significantly <cit.>[see "Scat alone" results in Tab. 1.].
    

    
    
  * Batch Normalization layer: This operation conditions the Hessian matrix of the convolution with respect to its parameters by performing the substitution 
    
    x_i → x^'_i =  x_i - M/S

    where M and S^2 are the empirical mean and variance.
    
    
  * Pointwise Non-linear Activation: Such as a Rectified Linear Unit (ReLU) ρ(x) = max(0,x). These nonlinearities are essential; otherwise, the network operation would be purely linear, incapable of separating "XOR" type samples. One other argument concerns the non-extinction of gradients.
    
    
  * Pooling layer: Average, max, etc.
    
    
  * Deconvolution and upsampling layer in the case of contracting and expansive architectures like in the case of .

If the architecture of  presents a linear sequence of blocks (convolution, batch-norm, ReLU), that of  proceeds in two parts: the first (contracting) chains convolutions and subsamplings as for a wavelet decomposition, while the second (expansive) performs the inverse transformation. There are other architectures that proceed in this way, such as Variational Auto-Encoders <cit.>. What distinguishes  is the fusion at each scale of the results of the two phases. This mechanism, called "skip connection", merges information from higher resolution with information from lower resolution. This feature allows correlating information at different scales.



 §.§ Learning the noise

Regarding the denoising task, the model (e.g., ) denoted D_θ (with θ being the set of parameters to optimize) is used to learn the noise. That is, if we denote y as the sample at the input of the network corresponding to the noisy signal and D_θ(y) as the output, then

    x̃(y) = y - D_θ(y)

is constrained to match the noise-free signal x. This means that the cost function to minimize is the mean square error, such that if we have n samples {x_i,y_i}_i≤ n then:

    θ^∗  = argmin_θ ∑_i=1^nx_i-x̃(y_i)^2

The models are trained with all levels of noise σ∈ [0,1], knowing that the pixels of the x_i are also in [0,1]. It should be noted that training with only one level of noise gives poorer results. In the following, I will only consider the results concerning the  architecture, but they are similar to those obtained with  <cit.>. A peculiarity used by the authors of this article is to remove bias terms from the layers (convolution and batchnorm), which gives the network the property of being first-order homogeneous in the sense of differentiable functions from ℝ^n→ℝ^n. Euler's theorem then tells us that

    D_θ(y) = (∇_y D_θ(y))  y

where the Jacobian of D_θ(y) appears.



The denoising result of a C^2 image (80× 80) to which Gaussian noise with a PSNR of 20 dB (σ = 0.1) was added is shown in Figure <ref>. The network was trained with 100,000 images. These results should be compared to those obtained in Section <ref> taking into account the reduction in size[Refer to the curves in Figure <ref> to appreciate the change in PSNR with size.] required for training the . The quality of the reconstruction by the neural network can also be appreciated as a function of the regularity order α of C^α images and the level of added noise (input PSNR) (Fig. <ref>). The curves from Figure <ref> are reported here for comparison. It is noted that the network performs better especially at the boundary, and it behaves better as the regularity order increases.


Questions arise: why do networks like  perform significantly better than previously developed methods? What do they learn? Is it generalizable to other types of images? What about the size of the training image batch? To address these questions, I will heavily rely on recent work by <cit.>.


 §.§ What does the network learn? Orthogonal basis and thresholding



If we consider the definition (Eq. <ref>) of the estimated x̃(y) of the noise-free signal x from the noisy signal y, and take into account the homogeneous nature of the network's response (Eq. <ref>), then we can write

    x̃(y) = (Id - ∇_y D_θ(y)) y = 𝒪_θ(y)  y

By diagonalizing the operator 𝒪_θ(y), denoting (λ_k(y),e_k(y))_k the set of eigenvalues and eigenvectors, we can express y in this basis by denoting c_k(y) = ⟨ y, e_k(y)⟩:

    y = ∑_k c_k(y) e_k(y)

Similarly, the denoised signal decomposes as

    x̃(y) = ∑_k λ_k(y) c_k(y) e_k(y)

This relationship is reminiscent of the thresholding formulas in Eqs. <ref>, <ref> while maintaining adaptivity to the geometry of y as seen for wavelets because the thresholding function ρ(c_k(y))=λ_k(y) c_k(y) depends on y. Thus, the network learns not only the basis but also the thresholding to apply to remove the noise.

In practice, what does this give us? Obtaining the eigenvalues and eigenvectors of 𝒪_θ(y) is done via the Singular Value Decomposition (SVD) from an image y with, for example, noise such that σ=0.1 (20 dB). The results are shown in Figures <ref> and <ref> for C^α images with α=2 and 5. First, we notice that the dot products decrease rapidly[Reminder: they are organized according to the decreasing values of the eigenvalues.], indicating that we are indeed dealing with a sparse decomposition. The eigenvalues also decrease rapidly, exploiting sparsity well, especially since their clear decrease is noticed for k≳ 10^2, i.e., as soon as the noise becomes noticeable. The most surprising aspect is the images of the eigenvectors: they are oscillating functions both along the contours and in the uniformly regular regions of the image (without crossing the boundary). In short, they adapt to the geometry.


These oscillating functions resemble solutions of a Laplace equation; they are called geometry-adaptive harmonic basis (GAHB). It can be seen that these bases are quite different from those used in Section <ref>. It must be understood that GAHB bases are adaptive; the basis vector depends on the signal (e_k(y)), whereas wavelet bases (or other "curvelets/bandlets, etc.) are not adaptive; the basis vector is known a priori (e_k). This is to be differentiated from the sparse nature, where the decomposition coefficients concentrate on local peculiarities (transients) of each signal: wavelet bases are sparse, while the Fourier basis is not for signals with transients. 







 §.§ Prediction of oracles and practical case



To understand the behavior of the network when performing denoising, let's take a small detour through a classical analysis called "oracles". If we take the expression of the denoised signal (Eq. <ref>), then the mean squared error (fixed signal x, noise z∼𝒩(0,σ^2), y=x+z) is given by:

    ℰ = 𝔼_z[x-x̃(y)^2] 
    = ∑_k 𝔼_z[(⟨ x, e_k(y)⟩ - λ_k(y) ⟨ y, e_k(y)⟩)^2]

To simplify, let's consider that e_k(y) is independent of y, as if we were in a fixed orthonormal basis like the wavelet bases in Section <ref>. Let c_k(x)=⟨ x, e_k⟩ and similarly for  
c_k(y) and c_k(z). Note that by linearity c_k(y)=c_k(x)+c_k(z), so we have

    ℰ    = ∑_k 𝔼_z[(c_k(x)(1-λ_k(y))-λ_k(y) c_k(z))^2]
       = ∑_k{ c^2_k(x)𝔼_z[(1-λ_k(y))^2] - 2 c_k(x)𝔼_z[c_k(z)(1-λ_k(y))λ_k(y)] + 𝔼_z[(λ_k(y) c_k(z))^2]}

At this stage, we have a problem to continue the calculation because λ_k(y) implicitly depends on z. So, we will use an oracle which will give us λ_k(x) instead. The result will be necessarily optimistic. It should also be noted that any decomposition of white Gaussian noise 𝒩(0,σ^2) over any orthonormal basis leads to random values of the coefficients according to 𝒩(0,σ^2). Therefore, 𝔼_z[c_k(z)]=0 and 𝔼_z[c^2_k(z)]=σ^2.   Then we have

    ℰ    =∑_k{ 
    c^2_k(x) (1-λ_k(x))^2 -2 c_k(x)(1-λ_k(x))λ_k(x) 𝔼_z[c_k(z)]
    + λ^2_k(x) 𝔼_z[c^2_k(z)]
    }
        = ∑_k{ c^2_k(x) (1-λ_k(x))^2 + σ^2 λ^2_k(x) }

Minimizing such an error amounts to minimizing each term, which leads to the prediction of the oracle:

    λ^om_k(x)    = c^2_k(x)/c^2_k(x)+σ^2  ⇒  ℰ^om =  ∑_k σ^2 c^2_k(x)/c^2_k(x)+σ^2 = σ^2 ∑_k λ^om_k(x)


Looking at equation <ref>, one might think of a harsher oracle that would use an adaptive "hard" thresholding (Eq. <ref>)

    λ^oh_k(x) = 
    1     c^2_k(x)≥σ^2 
    
    0    otherwise

Let's re-index the coefficients c_k(x) in a decreasing order such that

    c_1(x) ≥ c_2(x) ≥…≥ c_M(x)≥σ^2  ≥ c_M+1(x)…

The expression of the mean squared error of this oracle is equal to 

    ℰ^oh   = ∑_1≤ k≤ Mσ^2 + ∑_k≥ M+1^k_max c_k^2(x)  
    = ∑_k min(σ^2,c_k^2(x)) 
    = Mσ^2 + x-x_M ^2

where for the last expression, we simply expressed that x_M is the approximation of x with the first M coefficients (re-indexed), and k_max=N^2 in the case of an N× N image. Note that M is the cardinality of the following set:

    M=Card{k /  c^2_k(x)≥σ^2 }


We can relate ℰ^oh and ℰ^om, because ∀ (a,b)∈ℝ^+ we have

    min(a,b)/2 ≤ab/a+b≤min(a,b)

thus

    ℰ^oh/2≤ℰ^om≤ℰ^oh


Therefore, up to a factor of 2, we can write that

    ℰ^om∼ℰ^oh = Mσ^2 + x-x_M ^2

This is a classical result of the behavior of the mean squared error of denoising in an orthonormal basis <cit.>.

The first observation is that if the coefficients c_k(x) are zero for k>M, then the estimation error x-x_M ^2 is zero, and the denoising error is minimal, equal to Mσ^2 (corresponding to the variance of M independent random variables). This indicates that the more sparse the representation of the (noise-free) signal is (M small), the better denoising can be achieved.

We can also provide an estimate of the convergence rate of ℰ^oh (ℰ^om) in the case of a power-law decay of the coefficients <cit.>:

If there exist three strictly positive constants c_1,c_2,p such that

    c_1 k^-(p+1)≤ c_k^2(x) ≤ c_2 k^-(p+1)  ⇒  ℰ^oh∼σ^2p/(1+p)




















Now, to address the real case, we can't use the coefficients c_k(x) as the oracle, since the signal x is unknown to us. Therefore, we can only use thresholding on the coefficients c_k(y):

    λ^T_k(y) = 
    1     c^2_k(y)≥ T^2 
    
    0    otherwise

while still considering a fixed orthonormal basis. In fact, σ is replaced by T, and we can follow the previous developments, leading to 

    ℰ^T = M T^2 + x-x_M ^2

with M = Card{k /  c^2_k(y)≥ T^2 }, similar to expression <ref>. The problem then boils down to finding the best threshold in the given basis. This is where the result of D. L. Donoho and I. M. Johnstone <cit.> provides a practical solution that we have already used (e.g., Sec. <ref>):
(Donoho & Johnstone)


If T=σ√(2log d) with d being the size of the signal, then[the error by the oracle uses a "hard" thresholding given by σ.]

    ℰ^T≤ C_d(σ^2 + ℰ^oh)

with C_d=(1+2log d) and ℰ^oh = ∑_k min(σ^2,c_k^2(x)). Moreover, C_d and T are asymptotically optimal.

Therefore, when d≫ 1 (which is the case), then up to a constant, with the above thresholding, we can achieve oracle performance. Thus, we can write (C being a constant)

    ℰ^T≤ C  σ^2p/1+p
 

Now, <cit.> have shown that for C^α type images, the optimal quadratic error is such that

    ℰ^optim≤ C σ^2α/1+α

which means that one can try to find a basis such that 

    |⟨ x, e_k⟩|^2 ∼ k^-(1+α)

Moreover, if we translate the scaling of ℰ^optim in terms of PSNR, the slope of the PSNR of the optimally denoised image as a function of the PSNR of the noisy image[incidentally equal to -20log_10σ.] is α/(1+α).

This result motivated the search for orthonormal bases described in section <ref>. However, it is noteworthy that the result for "bandlets" (Sec. <ref>), which form the best sparse orthonormal basis for C^α images, gives <cit.> (σ^2 ≫ 1/d)

    ℰ^bandlet≤ C (|logσ|σ^2)^α/(1+α)

There is an additional |logσ| factor compared to the optimal result of <cit.>. 

Regarding the denoising result by the neural network, it is optimal: the decay of the coefficients c^2_k(x) is as expected (Fig. <ref>) as well as the slope of the PSNR as illustrated in Figure <ref>. Since theory cannot access the offset, we cannot go further to know if we can do better in exploiting these results. However, by analyzing the result of <cit.>, it is indeed the detection of contours in the presence of noise, and the exploitation of the regularity of the function along these contours that are the keys <cit.>. It is clear that the network manages to capture these regularities. The next section will show another example.



 §.§ The case of "celebrity" portraits

The next question arises: what are the performance of a  for "natural" images? <cit.> used face images from the  database[<https://mmlab.ie.cuhk.edu.hk/projects/CelebA.html>] <cit.>.

After training with 100,000 80× 80 images, Figure <ref> gives the results in a similar form to Figure <ref>. We can see that the eigenvectors again show the adaptation of the basis to the different contours present in the image. The decay (not shown) of the coefficients c^2_k(x) is of the same nature as those obtained with C^α images. Other examples are presented in the article, including pictures of bedrooms.



 §.§ Influence of the Database Size


In the previous sections, the  networks were trained with 100,000 80× 80 images. What if the size of the training batch is smaller? Figure <ref> answers the question for an extreme case where the database consists of 10 images. The image used for denoising is an image from this same batch[Otherwise, the results are even worse: the denoised image of a test image corresponds to a mixture that does not resemble the original image at all. However, for networks trained with 100,000 images, test images were used to produce denoising results and orthonormal bases.].





It can be observed that the denoising is almost perfect. However, it can be noticed that the eigenvectors do not exhibit the same property of oscillating functions respecting the contours of the image as observed in Figure <ref>. Similarly, Figure <ref> shows that the decay of the coefficients c^2_k(x) (i.e., |⟨ x, e_k(y)⟩|^2) as well as that of the eigenvalues λ_k(y) is not satisfactory at all when N=10: the decay is much less pronounced (lack of sparsity) and the thresholding much less effective.

Clearly, a sufficient number of training images is needed to obtain the qualities stated in the previous sections. Can we understand this?



 §.§ What does the network learn? The score!


In the previous sections, the signal x is considered as a deterministic process. We seek an approximation from the observation y which is the result of the corruption of x by additive white Gaussian noise: y=x+z with z∼𝒩(0,σ^2). In this section, we will consider a probabilistic framework where x is the realization of a random variable with probability density p_X(x).
By the way, it might seem odd to consider that the images used in the previous sections are the result of random processes. However, we have been accustomed to the idea that systems with a large number of degrees of freedom are better understood through Statistical Physics since the beginning of the 20th century and the Boltzmann-Gibbs theory. And in this theory, a distribution of Gibbs is generically formulated underlying the process as follows:

    p_X(x) = Z^-1 e^-U(x)

Moreover, the processes studied in Statistical Physics are certainly in most cases non-Gaussian (U(x) does not have just a quadratic term in x) but they are ergodic. Thus, it is easy to make the implicit connection between an ergodic process and a probability distribution. However, there are non-ergodic processes, at least because they are not stationary, like face images. However, since the underlying process is not ergodic, this may indeed challenge our intuition shaped by Statistical Physics. Therefore, let's try to break free from this constraint, without completely disregarding it.

In this context, the search for x̃(y), which best approximates the signal x, is done by minimizing the quadratic error in expectation over the two probability densities, that of y (via z) and that of x:

    x̃(y) = argmin_𝒪𝔼_x,y[x-𝒪(y)^2]

where 𝒪 is an operator from ℝ^d to ℝ^d (it transforms an image d=N× N into another image in the same space).

This is a Bayesian estimation whose following theorem gives the solution:
 
The Bayesian estimator of the denoised signal is given by (i=1,…,d)

    x̃(y)_i = 𝔼_x[x_i|y=(y_1,…,y_d)] = ∫ x_i p_i(x_i|y) dx_i
























The main problem, even if we have a compact form, lies in the fact that p_i(x_i|y), whether we like it or not, depends on the a priori p_i(x_i) (for all i=1,…,d). Moreover, the integral calculation must be performed in high dimension. This is why initially <cit.> sought a linear operator of the form x̃(y) = L.y.

But as pointed out by <cit.> and more recently by <cit.>, for Gaussian white noise 𝒩(0,σ^2) independent of the unnoisy signal x, we have the following result independent of the prior on x:

The estimator x̃(y) of the signal x, such that y=x+z (z∼𝒩(0,σ^2)), is given by

    x̃(y) = y + σ^2 ∇_y log p(y)



















This theorem could be just as frustrating as the Bayesian estimator theorem invoking conditional probability, but it also tells us that if we have a good denoiser, like the well-trained  network, then knowing x̃(y) and the noise variance[It could be estimated with the signal itself as we did in the context of denoising in wavelet bases for example.], then we have an estimator of the score ∇_y log p(y) of the distribution of observables. And it is this ingredient that proves to be crucial for generating new samples from training images.



 §.§ Score-Based Generation

To generate samples from p(x), we proceed as follows (see e.g., <cit.>). Starting from a sample x assumed to be drawn from the distribution p(x), we progressively transform it (steps indexed by t∈ [0,T]) until the underlying distribution is very simple, namely a Gaussian with σ=1. This is the forward process in Figure <ref>. The idea is then to reverse the process (the backward mode) by generating a new sample from a new instance of noise.


The transport process used in the Score Diffusion algorithm is the Ornstein-Uhlenbeck equation <cit.>. The forward process follows the following stochastic differential equation:

    dx_t = -x_t dt + √(2) dB_t    t∈[0,T]

where dB_t is a Wiener process (Brownian motion). The discretization gives (constant time step δ):

    x_t+δ = (1-δ)x_t + z̃    z̃iid∼𝒩(0,2δ Id)

It is a Gaussian and Markovian process whose solution is given by

    x_t = x_0 e^-t + (1-e^-2t)^1/2_σ_t z      z iid∼𝒩(0,Id)

The associated probability density function is given by the expression 

    p_t(x_t) = ∫ p_t(x_t,x_0) dx_0 = ∫ p_t(x_t|x_0) p_0(x_0)dx_0

where according to the Fokker-Planck representation of the Ornstein-Uhlenbeck process, we have

    p_t(x_t|x_0) = 1/(2πσ_t^2)^d/2 e^ -x_t-x_0e^-t^2/2σ_t^2= 
    g_σ_t(x_t-x_0e^-t)

In fact, p_t(x_t) is the result of the convolution of the initial distribution p_0(x_0) with this Gaussian which gradually tends towards 𝒩(0,1), producing a progressive smoothing that also makes p_t(x_t) tend towards the normal distribution.

Now, let's move on to the backward process, which is also the result of a stochastic differential equation (damped Langevin): 

    dx_T-t = (x_T-t + 2∇_x log p_T-t(x_T-t))dt + √(2) dB_t    t∈[0,T]

which in discrete form (δ>0 constant time step) gives the relation:

    x_t-δ = (1-δ)x_t + 2δ s_t(x_t) + √(2δ) z      z iid∼𝒩(0,Id)

where s_t(x_t) = ∇_x log p_t(x_t) is the score.

Here, there are two possible approaches: either we have a model for p_t (e.g., <cit.> in the case of generating φ^4 fields, for example) that can exploit regularity patterns and assumptions about correlations to provide a constrained parametric framework for learning the score, or we use a denoising neural network for the same purpose. So, imagine that we have a batch of input images (x_0), at each step t of the forward pass, we have new noisy images with variance σ_t, which can be used to train a network[This is the method of score matching <cit.>.] like a  to denoise them and thereby obtain the score s_t(x_t) at step t. Then, starting from a pure noise instance (x_T), we perform the iterations of the backward process, and at t=0 we obtain a new sample from p_0.

The real question when we see such images produced by these generative models is: are we sampling from a probability density? In other words, has the model learned p(x) with x being a face image, and what is the nature of this distribution p? The alternative would probably be that the new images are sophisticated patchworks of database images? A first answer was given by <cit.>.



 §.§ Memorization vs Generalization



Generating new samples from a distribution p(x) using the Score Diffusion method briefly described in the previous section raises questions. In particular, the quality of the generated images and their real novelty may be questionable. To try to answer both questions, consider the following experiment: prepare two disjoint sets of N images ((B_1(N),B_2(N))) from the same database (i.e., CelebA), and train two  networks (S_1(N),S_2(N)), S_i(N) with B_i(N). This is repeated for batch sizes N=1,10,100,…,10^5. Then, for a fixed N, draw a single image of pure noise (the image x_T), and perform the backward process described in the previous section (Eq. <ref>), where the score s_t is estimated using S_1(N) and S_2(N). This yields 2 samples denoted x_0^(i,N) to indicate that they are obtained by the two models S_i(N).

Figure <ref> shows a result of the experimentation. The questions that arise then are: does x_0^(i,N) resemble closely or remotely an image from B_i? does x_0^(1,N) resemble x_0^(2,N)? The first question raises the issue of whether the obtained sample is possibly a mixture of images from the training database, or is it truly new. The second question concerns the uniqueness of the transformation process from the normal distribution to a hypothetical distribution p_0.

As long as N<O(10,000), S_1(N) and S_2(N) produce different images but ones that resemble an image from the database more closely as N decreases. For N=10,000, the generated images x_0^(1,N) and x_0^(2,N) are different from the database. However, they are different from each other, which raises questions because it seems that two transformation processes are at work (recall x_T is the same noise image given to both backward processes). However, as soon as N=100,000, the synthesized images x_0^(1,N) and x_0^(2,N) by S_1(N) and S_2(N) are not only different from the images closest to the two databases as before, but they are also very close to each other (red rectangle), indicating that a single process is indeed learned.

We are witnessing a transition from models that behave like large memories to models that effectively learn a probability distribution as long as the number of training images is sufficient. For 80× 80 images, for a  with 7.6 10^6 parameters, the transition seems to occur around N ∼ 10^5, and it occurs around N =O(10^4) for 40× 40 images. However, it is not certain that a simple extrapolation of the training batch size would suffice for high-resolution images. That being said, one should not stop at the criterion of similarity between the generated sample and the closest sample from the training database to assess the quality of a generator.



So, a properly trained denoising neural network indeed learns a distribution p(x), whether it's for face photos, indoor photos, etc., for which it was not a priori imaginable that there would be such an underlying probability density. In high dimensions, there are phenomena of probability concentration on typical sets <cit.>, so it is possible that the description of the structures of this set in the case of photos is not as complex as it may seem a priori. Indeed, if  networks manage to overcome the curse of dimensionality, it is because they have captured something, a form of regularity of the underlying distribution: "what is it?" that is the question.



§ SUMMARY


In this article, after briefly reviewing Fourier analysis and classical wavelet analysis, as well as studying the denoising of synthetic images with fixed regularity (C^α), we examined why, in the period from approximately 1995 to 2005, research focused on constructing orthonormal wavelet bases attempting to adapt to the signal's geometry. However, approximation errors in these bases did not reach the theoretical optimum. The advent of neural networks around the 2010s changed the landscape. For example, the use of deep networks such as  and  (2015) is an example. On the other hand, the technique of generating samples from a density p(x) by Score Diffusion has shed new light. The score ∇_x log p(x), a concept developed by R. Fisher, can be directly learned by a network performing perfect denoising. This score, and thus this network, allows to invert a Markov process which, starting from a sample drawn from a simple d-dimensional normal distribution, produces a sample from p(x). Studying the performance of such networks on C^α images shows that they achieve the desired optimum. The wavelet bases produced adapt perfectly to the geometry of the images. This property is preserved when applied to face photos, or even indoor photos. They are therefore more efficient than "traditional" methods. Moreover, by analyzing according to the size of the training samples, it has been shown that from a minimum size, the networks can learn an underlying probability density for natural images, as surprising as it may seem. These results open perspectives ranging from the field of research in mathematics to philosophical questions about the existence of universal underlying structures, through very practical considerations on the development of such models so that they are indeed able to generalize.



§ ACKNOWLEDGEMENTS


Zahra Kadkhodaie and Florentin Guth gave me access to their pre-trained  models as well as the code for generating C^α images. With them and Gadriel Peyré, I was able to exchange ideas to develop this note. Nasser Eslahi gave me access to the  library (v2.1.3) and to his own  algorithm for using curvelets.



§ CODES


I used the  <cit.> and  <cit.> environments. An implementation of the wavelet steerable library by <cit.> is available on the following  site: <https://github.com/TetsuyaOdaka/SteerablePyramid>. The  algorithm by <cit.> is available on the site <https://www.ipol.im/pub/art/2014/86/>.



tocsectionRéférences


elsarticle-harv 




















