Hui Shi
email: hui.shi@math.u-bordeaux.fr

Yann Traonmilin
email: yann.traonmilin@math.u-bordeaux.fr

Jean-François Aujol
email: jean-francois.aujol@math.u-bordeaux.fr

Compressive learning for patch-based image denoising *

Keywords: Image denoising, Compressive learning, Sketching, Optimization, AMS subject classifications. 68U10, 94A08, 49N30

The Expected Patch Log-Likelihood algorithm (EPLL) and its extensions have shown good performances for image denoising. The prior model used by EPLL is usually a Gaussian Mixture Model (GMM) estimated from a database of image patches. Classical mixture model estimation methods face computational issues as the high dimensionality of the problem requires training on large datasets. In this work, we adapt a compressive statistical learning framework to carry out the GMM estimation. With this method, called sketching, we estimate models from a compressive representation (the sketch) of the training patches. The cost of estimating the prior from the sketch no longer depends on the number of items in the original large database. To accelerate further the estimation, we add another dimension reduction technique (low-rank modeling of the covariance matrices) to the compressing learning framework. To demonstrate the advantages of our method, we test it on real large-scale data. We show that we can produce denoising performances similar to performances obtained with models estimated from the original training database using GMM priors learned from the sketch with improved execution times.

∼ N (0, σ 2 I N ). In the last two decades, non local patch-based methods have been proven successfull for denoising. Methods such as Piecewise Linear Estimators [START_REF] Yu | Solving inverse problems with piecewise linear estimators: From gaussian mixture models to structured sparsity[END_REF][START_REF] Aguerrebere | Single shot high dynamic range imaging using piecewise linear estimators[END_REF], BM3D [START_REF] Dabov | Image denoising by sparse 3-d transformdomain collaborative filtering[END_REF][START_REF] Lebrun | An analysis and implementation of the bm3d image denoising method[END_REF] or NL-Bayes [START_REF] Lebrun | A nonlocal bayesian image denoising algorithm[END_REF][START_REF] Lebrun | Implementation of the "non-local bayes" (nl-bayes) image denoising algorithm[END_REF][START_REF] Xu | Patch group based nonlocal self-similarity prior learning for image denoising[END_REF][START_REF] Laus | A nonlocal denoising algorithm for manifold-valued images using second order statistics[END_REF] are examples of non-local methods [START_REF] Buades | A non-local algorithm for image denoising[END_REF]. In patch-based image denoising, the noisy image v is divided into small patches {v i } M i=1 . Each patch v i ∈ R P (P is the patch size) can be seen as a vector in a high dimensional space. The denoising problem is considered on each patch:

(1.2)

v i = u i + w i ,
and a corresponding denoised version u * i of the true values u i are estimated. To overcome the ill-posedness of this inverse problem, various denoising methods [START_REF] Levin | Natural image denoising: Optimality and inherent bounds[END_REF][START_REF] Lebrun | A nonlocal bayesian image denoising algorithm[END_REF][START_REF] Lebrun | Implementation of the "non-local bayes" (nl-bayes) image denoising algorithm[END_REF][START_REF] Houdard | High-dimensional mixture models for unsupervised image denoising (HDMI)[END_REF] consider patch models within a Bayesian framework. According to the Bayes' theorem, the objective is to find u * i which maximizes the posterior probability distribution f (u i |v i ) under the prior p(u i ). The Maximum A Posteriori (MAP) problem is formulated as

(1.3) u * i = arg max u i ∈R P f (u i |v i ) = arg max u i ∈R P f (v i |u i )p(u i ) ∝ arg max u i ∈R P e -u i -v i 2 2σ 2 p(u i )
where • denotes the 2 -norm. This yields (1.4) u * i = arg min

u i ∈R P u i -v i 2 2σ 2 -log(p(u i )).
Ideally, the choice of the prior distribution should be determined by the nature of the image to be estimated. In practice, Gaussian Mixture Models (GMM) [START_REF] Yu | Solving inverse problems with piecewise linear estimators: From gaussian mixture models to structured sparsity[END_REF][START_REF] Wang | Sure guided gaussian mixture image denoising[END_REF][START_REF] Houdard | High-dimensional mixture models for unsupervised image denoising (HDMI)[END_REF] have shown their effectiveness. With the GMM prior, the solution of problem (1.4) can be approximated by a Wiener filter solution (see subsection 2.2).

Among these various non-local denoising methods, the Expected Patch Log-Likelihood algorithm (EPLL) [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] occupies a central position due to its efficient denoising performance.

A large number of works build on the original EPLL formulation to deal with more general prior or go beyond the denoising problem [START_REF] Feng | Image restoration via efficient gaussian mixture model learning[END_REF][START_REF] Papyan | Multi-scale patch-based image restoration[END_REF][START_REF] Cai | Image denoising via patch-based adaptive gaussian mixture prior method[END_REF][START_REF] Luo | Adaptive image denoising by mixture adaptation[END_REF][START_REF] Ren | Example-based image synthesis via randomized patch-matching[END_REF][START_REF] Sulam | Expected patch log likelihood with a sparse prior[END_REF][START_REF] Deledalle | Image denoising with generalized gaussian mixture model patch priors[END_REF][START_REF] Saint-Dizier | A unified view on patch aggregation[END_REF]. EPLL uses a GMM prior learned from a very large set of patches extracted from clean images. The key to the success of EPLL is to find a good prior distribution. Since in practice patch sizes are typically greater than 5 × 5, estimating prior distributions in such a high-dimensional space is a difficult task. Moreover, to estimate the best possible model, we need to maximize the redundancy of structural information and use training databases as large as possible. As the traditional empirical minimization approaches require access to the whole training dataset, when the collection size is large, the learning process can be extremely costly. For instance, in the case of the classical learning method Expectation Maximization (EM), the memory consumption and computation time depend on the size of the database (see section 3).

Leveraging ideas from compressive sensing [START_REF] Foucart | A mathematical introduction to compressive sensing[END_REF] and streaming algorithms [START_REF] Cormode | An improved data stream summary: the count-min sketch and its applications[END_REF], R. Gribonval et al. propose a sketching method [START_REF] Keriven | Sketching for large-scale learning of mixture models, Information and Inference[END_REF][START_REF] Gribonval | Sketching datasets for large-scale learning (long version)[END_REF][START_REF] Gribonval | Sketching data sets for large-scale learning: Keeping only what you need[END_REF][START_REF] Gribonval | Statistical learning guarantees for compressive clustering and compressive mixture modeling[END_REF][START_REF] Gribonval | Compressive statistical learning with random feature moments[END_REF] to compress the training database.

This scalable technique compress the whole training collection into a fixed-size representation (a vector): a sketch of the training dataset before learning. The sketch captures the necessary information for the considered learning task. For certain mixture model estimation, it is then possible to learn their parameters directly from the sketch, without access to the original dataset. Hence the space and time complexity of the learning algorithm no longer depends on the original database size, but only on the size of the sketch which is linked to the dimensionality of the model. Sketching has been already used successfully in machine learning [START_REF] Schellekens | Compressive classification (machine learning without learning)[END_REF][START_REF] Gribonval | Compressive statistical learning with random feature moments[END_REF][START_REF] Keriven | Compressive k-means[END_REF][START_REF] Chatalic | Large-scale high-dimensional clustering with fast sketching[END_REF][START_REF] Byrne | Sketched clustering via hybrid approximate message passing[END_REF][START_REF] Permiakova | Sketched stochastic dictionary learning for large-scale data and application to high-throughput mass spectrometry[END_REF], generative networks [START_REF] Schellekens | Compressive learning of generative networks[END_REF], source localization [START_REF] Fontaine | Scalable source localization with multichannel α-stable distributions[END_REF][START_REF] Fontaine | Sketching for nearfield acoustic imaging of heavy-tailed sources[END_REF], independent component analysis [START_REF] Sheehan | Compressive independent component analysis[END_REF] and depth imaging [START_REF] Sheehan | A sketching framework for reduced data transfer in photon counting lidar[END_REF]. In [START_REF] Keriven | Sketching for large-scale learning of mixture models, Information and Inference[END_REF], the sketching is implemented and evaluated on synthetic data to estimate a GMM with diagonal covariances. It is shown that on large synthetic data, for the estimation of GMM, the sketching produces precise results while requiring fewer memory space and computations. In this work, we explore the sketching method in the image patches context where GMM with full covariance must be estimated from the compressed database.

Due to the curse of dimensionality, it is computationally expensive to manipulate the GMMs' covariance matrices. In [START_REF] Renna | Reconstruction of signals drawn from a gaussian mixture via noisy compressive measurements[END_REF], the authors show that most natural images and videos can be represented by a GMM with low-rank covariance matrices. The experiments have also shown the efficiency of low-rank covariance matrices applied to image denoising [START_REF] Parameswaran | Accelerating gmm-based patch priors for image restoration: Three ingredients for a 100× speed-up[END_REF], image inpainting, high-speed video and hyperspectral imaging [START_REF] Yang | Video compressive sensing using gaussian mixture models[END_REF]. This motivates us to use such low-rank covariances in the GMM modeling of patches and extend the sketching framework accordingly to gain computational speedup and to manage the modeling of the image patches in the most possible flexible way.

1.1. Contributions. A preliminary and short version of this work has appeared in [START_REF] Shi | Sketched learning for image denoising[END_REF]. In this paper, we provide a more detailed version of this work with a final consolidated version of the proposed learning algorithm, validated by extended numerical experiments.

Figure 1 summarizes the principle of our approach. We first construct a sketch by averaging random Fourier features computed over the whole image patch database. Then the model parameters are learned directly from the sketch by our Low-rank Continuous Orthogonal

Matching Pursuit (LR-COMP) algorithm without access to the original database. Finally, the learned model is used with a Bayesian method (EPLL) for the denoising task. Our the low-rank approximation of the covariance matrices, we lighten the computation burden in the denoising process while keeping good denoising performances.

• We demonstrate the performance of our approach on real large-scale data (over 4

millions training samples of patch size of 7 × 7) for the task of patch-based image denoising. We show that using models trained with the compressed database, we can obtain similar denoising performances compared to the models obtained with the classical EM algorithm. To the best of our knowledge, this is also the first time that the sketching framework has been applied with such high dimensional GMMs.

• Computationally, we estimate the model from a compressed database which is about 1000 times smaller than the original patch database. It leads to running time approximately two times faster compared to the EM method.

The paper is organized as follows. Section 2 is a reminder of the EPLL framework. Then we review the EM algorithm in section 3. In section 4, we explain the compressive learning method. In section 5, we focus on explaining how to adapt the sketching framework to learn a GMM in the image patch context. We also interpret the extension to low-rank covariances and the implementation details of the adapted learning algorithm LR-COMP. In section 6, we provide numerical experiments that demonstrate the performance of our approach. Some conclusions and tracks for further works follow in section 7.

2. Image denoising with EPLL. We review in this section the Expected Patch Log-Likelihood (EPLL) framework for image denoising. EPLL is a patch-based image restoration algorithm introduced by Zoran and Weiss [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]. The EPLL framework restores an image u by performing the following maximum a posteriori (MAP) estimation over all N patches:

(2.1) u * = arg min

u∈R N P 2σ 2 u -v 2 - N i=1 log(p(P i u))
where P i : R N -→ R P is a linear operator that extracts a patch of P pixels centered at the position i, typically P = 7 × 7. The function p(•) is the density of the prior probability distribution of the patches. Note that in practice, we assume that patches are distributed independently.

2.1. Optimization. Due to the non-convexity of p(•), direct optimization of the problem may be difficult. The authors of EPLL propose to perform the optimization with "halfquadratic splitting" [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF]. By introducing N auxiliary unknown vectors z i ∈ R P and a denoising parameter β > 0, the problem is then considered as:

(2.2)

u * = arg min u∈R N z 1 ,...,z N ∈R P P 2σ 2 u -v 2 + β 2 N i=1 P i u -z i 2 - N i=1 log(p(z i )).
The optimization (2.2) is accomplished by alternating the minimization of u and z i .

• Solving u for fixed z i -Problem (2.2) turns into a linear inverse problem with the Tikhonov regularization. It has a closed form solution: û = arg min

u∈R N P 2σ 2 u -v 2 + β 2 N i=1 P i u -z i 2 = (I + βσ 2 P N i=1 P T i P i ) -1 (v + βσ 2 P N i=1 P T i z i ) (2.3)
with N i=1 P i P i = P I, where P is the number of patches overlapping each pixel.

Hence we have

(2.4) û = (I + σ 2 βI) -1 (v + σ 2 β zi )
where zi := ( N i=1 P i P i ) -1 N i=1 P i z i = 1 P N i=1 P i z i is the image after averaging all overlapping patches z i .

• Solving z i for fixed u -The minimization problem (2.2) is separable with respect to the latent variable z i . It means that for each z i we solve a patch MAP estimation under the patch prior p(z i ), i.e. for all i,

(2.5) ẑi = arg min

z i ∈R P β 2 P i û -z i 2 -log(p(z i )).
The solution of this problem depends on the choice of patch prior p(•).

2.2. Denoising with a GMM prior. EPLL assumes that the prior is a finite Gaussian mixture model (GMM) with zero-mean on centered patches: the empirical mean estimated from noisy patches are removed before the denoising process (2.5) and added back in the end.

We consider that a zero-mean patch x ∈ R P is a random vector generated from a distribution with density p(x) defined as

(2.6) p(x) = K k=1 α k N P (x; 0, Σ k )
where K is the number of Gaussian components and α k ≥ 0 are weights of each component such that K k=1 α k = 1. The function N P (x; 0, Σ k ) denotes the density of a Gaussian distribution with zero-mean with covariance Σ k ∈ R P ×P . Recall that the zero-mean Gaussian distribution density is:

(2.7)

N P (x; 0, Σ k ) = 1 (2π) P/2 |Σ k | 1/2 e -1 2 x T Σ -1 k x .
Hence, under the GMM prior, Problem (2.5) turns to:

(2.8) ẑi = arg min

z i ∈R P β 2 P i û -z i 2 -log( K k=1 α k N P (z i ; m i , Σ k ))
where we supposed that the mean m i are correctly estimated by the empirical mean of noisy patches.

This problem cannot be solved in closed form as the second term is the logarithm of a sum of exponential. In [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF], the authors proposed to solve this problem by keeping only one Gaussian component. For a given centered patch zi = P i û -m i , we chose the component k * i that maximizes the posterior probability p(k i |z i ). This leads to computationally efficient implementations. [START_REF] Tran | Number of useful components in gaussian mixture models for patch-based image denoising[END_REF] also justified that only one component is required for good reconstructions.

The index k * i is chosen by

k * i = arg max 1≤k i ≤K p(k i |z i ) = arg max 1≤k i ≤K p(k i )p(z i |k i ) = arg min 1≤k i ≤K -2 log α k i + log Σ k i + 1 β I P + z i (Σ k i + 1 β I P ) -1 zi (2.9)
where α k i and Σ k i are the weights and the covariance matrices of the kth Gaussian component for the given patch zi . With k * i (instead of a sum of K components), the solution of (2.8) is then a Wiener filtering solution:

(2.10) ẑi = (Σ k * i + 1 β I P ) -1 Σ k * i zi +m i .

Eigenspace implementation of EPLL.

The matrix inversions in (2.9) and (2.10) can be done efficiently by using the singular value decomposition over the covariance matrices.

We denote (2.11)

Σ k = U k Λ k U T k , with U k ∈ R P ×P an unitary matrix and Λ k = diag(λ (k) 1 , ..., λ (k)
k * i = arg min 1≤k≤K -2 log α k + P j=1 log(λ (k) j + 1 β ) + [ṽ (k) i ] 2 j λ k j + 1 β where (2.12) ṽ(k) i = U T k zi .
Then (2.10) leads to

(2.13) ẑi = U k * i S k * i U T k * i zi +m i = U k * i S k * i ṽ(k * i ) i +m i with (2.14) S k * i = diag   λ (k * i ) j λ (k * i ) j + 1 β   j=1,...,P .
3. Learning a GMM with EM. The Expectation-Maximization (EM) algorithm is a classical mixture estimation approach. This algorithm starts with some initial estimates of model parameters and then iteratively updates the estimate until the the estimates are not changing much. See Appendix B for the details of the EM algorithm. In each iteration, it carries out two steps: the E-Step (expectation step) and the M-Step (maximization step). In E-Step, using the current estimate of the parameters, we evaluate the posterior probabilities. In the M-Step we compute parameters that maximize the probabilities found on the E-Step. These estimated parameters are then used to determine the distribution of the latent variables in the next E-Step.

As for the time complexity of one iteration of this algorithm, it is linear in the number of model components K and the number of elements in the database n. However it is cubic with respect to the dimensions P due to the fact that we need to inverse the covariance matrix when calculating the density in E-Step. Thus, when estimating a K-components GMM on a database of n elements of dimension P , the computational complexity of one iteration of the EM algorithm is O(nKP 2 + KP 3 ). The major criticism of the EM algorithm is that when dealing with a large dataset, it often converges slowly. To address this problem, researchers have developed various variations of the traditional EM algorithm [START_REF] Mclachlan | The EM algorithm and extensions[END_REF][START_REF] Varadhan | Simple and globally convergent methods for accelerating the convergence of any em algorithm[END_REF]. Learning parameters using EM technique face computational issues linked to the size of the dataset and the number of parameters to estimate, which would make the use of (very) large image patches databases impractical. In the next section we will see an alternative manner to learn parameters using compressive learning.

Sketching.

Sketching is a dimensionality reduction method. The principle is to compress the whole dataset massively before learning. First, the dataset χ

= {x i } n i=1 is summa- rized into a vector y ∈ C m (m n) called the sketch: (4.1) y := Sketch(χ).
Note that the computation of the sketch can be performed in a distributed manner. Then we apply a learning procedure Υ that allows us to learn an estimate Ψ * of some statistical parameters Ψ of the dataset directly from the sketch y, namely

(4.2) Ψ * = Υ(y) = Υ(Sketch(χ))
More specifically, learning from the sketch corresponds to a minimization problem

(4.3) Ψ * ∈ arg min Ψ E(y, Ψ)
where the energy of the model E(•, •) quantifies the fit between the sketch y and the parameter Ψ. In the context of statistical learning, the energy E can be seen as a proxy of the empirical risk. The principle of sketching is summarized in Fig. 2. vector space while keeping all the necessary information of the dataset. Mathematically, given a linear sketching operator S:

(4.4) S :D -→ C m y = Sf
and for some finite K ∈ N * , we define a K-sparse model f Θ,α ∈ D:

(4.5) f Θ,α = K k=1 α k f θ k
where f θ k ∈ D are elementary measures parametrized by θ k , α k ≥ 0 for all components and K k=1 α k = 1. We can express the vector y as

(4.6) y = Sf Θ,α = K k=1 α k Sf θ k .
In practice we only have access to the empirical probability distribution f =

1 n n i=1 δ x i
where δ x i is a unit mass at x i . So we can define the empirical sketch as ŷ = 1 n S n i=1 δ x i .

The goal of the sketching framework is to recover f Θ,α from y, hence we do the following minimization to estimate the parameters

(4.7) (Θ * , α * ) ∈ arg min Θ=(θ k ) K k=1 α∈R K ,α k ≥0, K k=1 α k =1 Sf Θ,α -ŷ 2 2 .
The objective of sketched learning algorithms is to minimize a datafit functional between the compressed database and the sketch of the estimation. In other words, our aim is to find parameters α, Θ such that the sketch of the probability distribution parameterized by α, Θ is the closest to the empirical sketch ŷ.

Recovery guarantees.

It was shown in [START_REF] Gribonval | Statistical learning guarantees for compressive clustering and compressive mixture modeling[END_REF] that we can guarantee theoretically the success of this estimation with a condition on the sketch size. These guarantees necessitate a "Lower Restricted Isometry Property" (LRIP) of the sketching operator. This property, is verified with high probability, for GMM with sufficiently separated mean and random Fourier sketching as long as the sketch size m ≥ O(K 2 dpolylog(K, d)), i.e. when the size of the sketch essentially depends on the parameters K (the number of components) and d (the model dimension). Empirical results seem to indicate that for d tot the total number of parameters, a database size of the order of d tot is sufficient: in the case of estimating a GMM with diagonal covariance matrices [START_REF] Keriven | Sketching for large-scale learning of mixture models, Information and Inference[END_REF], the authors observe that the quality of the reconstruction exhibits a sharp phase-transition with respect to the sketch size m. This phase transition happens for m proportional to d tot . In our model, d tot = K(P r + 1). The excess risk of the GMM learning task is then controlled by the sum of an empirical error term and a modeling error term. This guarantees that the estimated GMM approximates well the distribution of the data [START_REF] Gribonval | Sketching datasets for large-scale learning (long version)[END_REF].

Note that EPLL uses a zero-mean GMM as the patch prior, therefore, during the learning process, the patches are centered before sketching and we do not estimate the mean of Gaussians. In our case, the sketched GMM learning problem reduces to the estimation of the sum

of k zero-mean Gaussians with covariances Θ = (Σ k ) K k=1 , i.e f Θ,α = K k=1 α k g Σ k where g Σ is
the zero-mean Gaussian measure with covariance Σ. In this context, the notion of separation used to prove guarantees in [START_REF] Gribonval | Statistical learning guarantees for compressive clustering and compressive mixture modeling[END_REF] does not hold. We still show empirically that the sketching process is successful without this separation assumption.

4.3. Design of sketching operator: randomly sampling the characteristic function. In [START_REF] Keriven | Sketching for large-scale learning of mixture models, Information and Inference[END_REF], the sketch is a sampling of the characteristic function (i.e the Fourier transform of the probability distribution f ). Recall that the characteristic function ψ f of a measure f is defined as:

(4.8) ψ f (ω) = R d e -iω T x df (x) ∀ω ∈ R d .
The sketching operator is therefore expressed as:

(4.9) Sf = [ψ(ω 1 ), ..., ψ(ω m )] T
where {ω 1 , ..., ω m } is a set of well chosen frequencies. In the spirit of Random Fourier Sampling, the authors of [START_REF] Keriven | Sketching for large-scale learning of mixture models, Information and Inference[END_REF] propose to draw the frequencies from a probability distribution, i.e.

(ω 1 , ..., ω m )

i.i.d.
∼ ∆. The choice of frequencies is essential to the success of sketching, and we will discuss it in details in subsection 5.1.

Sketching image patches.

In this section, we adapt the sketching framework to the context of image patches. As when using the classical EM algorithm, the GMM learning from sketch is performed under the assumption that the training patches are i. e -iω T x j with ω ∈ R P .

Thus the empirical sketch ŷ is expressed as

(5.2) ŷ = [ ψ(ω 1 ), ..., ψ(ω m )] T = 1 n   n j=1 e -iω T 1 x j , ..., n j=1 e -iω T m x j   T .
In other words, a sample of the sketched database is a P -dimensional frequency component calculated by averaging over patches (not to be mixed with usual 2D Fourier components of images). Thanks to the properties of the Fourier transform of Gaussians, the sketch of a single zero-mean Gaussian component g Σk at frequency ω l is

(5.3) (S(g Σk )) l = ψ g Σ k (ω l ) = e -1 2 ω T l Σ k ω l .
Thus, given the weights α = (α k ) K k=1 and the covariance matrices Σ = (Σ k ) K k=1 , the sketch of

a zero-mean GMM f Σ,α = K k=1 α k g Σ k is (5.4) y = [S(f Σ,α ) l ] l=1,...,m = K k=1 α k e -1 2 ω T l Σ k ω l l=1,...,m
.

As a consequence, denoting P SD P the set of P × P positive symmetric definite matrices, the problem (4.7) of estimating GMM parameters becomes

(5.5) (Σ * , α * ) ∈ arg min Σ=(Σ k ) K k=1 ,Σ k ∈P SD P α∈R K ,α k ≥0, K k=1 α k =1 ŷ -Sf Σ,α 2 2 , i.e.
(5.6)

(Σ * k ) K k=1 , (α * k ) K k=1 ∈ arg min Σ k ∈P SD P ,∀k α k ≥0, K k=1 α k =1 m l=1 1 n n j=1 e -iω T l x j - K k=1 α k e -1 2 ω T l Σ k ω l 2 .
In practice, the positive definite constraint in the optimization problem is hard to enforce directly on the space of P ×P matrices (as previous work only considered diagonal covariances, it was not an issue). Our method, based on the Burer-Monteiro, permits us to respect the PSD constraint by recasting the covariance estimation problem as an optimization over R P ×P without constraint (see subsection 5.2 for more details).

5.1. Frequency sampling. The design of the probability distribution ∆ for sampling the frequencies {ω 1 , ..., ω m } is essential to the success of sketching. In our work, we draw frequencies from the adapted radius frequency distribution proposed in [START_REF] Keriven | Sketching for large-scale learning of mixture models, Information and Inference[END_REF]. The adapted radius heuristic proposes to sample ω as

(5.7) ω = Rϕ where R ∈ R + is the norm of ω and ϕ ∈ R P is the random direction. The radius R is chosen with a radius distribution R ∼ p R (R; η) = ((ηR) 2 + 1 4 (ηR) 4 ) 1 2 e -1 2 (

Extension to low-rank covariances.

Bayesian MAP theory permits to use a GMM with degenerate covariance matrices as a denoising prior. In this case, the prior is given only in the union of subspaces spanned by the r leading eigenvectors of the K covariance matrices of the GMM. The experiments [START_REF] Parameswaran | Accelerating gmm-based patch priors for image restoration: Three ingredients for a 100× speed-up[END_REF][START_REF] Renna | Reconstruction of signals drawn from a gaussian mixture via noisy compressive measurements[END_REF] have shown that we can use low-rank covariance matrices for denoising while keeping good performance. This motivates us to approximate the covariance matrices in the GMM prior by low-rank matrices.

Following classical Burer-Monteiro method [START_REF] Burer | Local minima and convergence in low-rank semidefinite programming[END_REF][START_REF] Chi | Nonconvex optimization meets low-rank matrix factorization: An overview[END_REF] in low-rank matrix estimation, we parameterize Σ k by its factors

X k : Σ k = X k X T k . We define f X,α the density function of a zero-mean GMM with X = (X k ) K k=1
, where X k is a factor of a covariance matrix.

Supposing that ŷ -Sf X,α 2 2 has a minimizer, we approximate the minimization (5.5) by

(5.8) ( X, α) ∈ arg min X=(X k ) K k=1 ,X k ∈R P ×r α∈R K ,α k ≥0, K k=1 α k =1 ŷ -Sf X,α 2 2 , i.e.
(5.9) ( Xk ) K k=1 , (α k ) K k=1 ∈ arg min

X k ∈R P ×r ,∀k α k ≥0, K k=1 α k =1 m l=1 ŷ - K k=1 α k e -1 2 ω T l X k X T k ω l 2
where X = { X1 , ..., XK } is the collection of factorized rank reduced covariances. With the following proposition, we show that the difference between the costs minimized in (5.5) and

(5.9) (the full rank and the low-rank cases, respectively) is associated with the smallest eigenvalues of the covariance matrices. We qualitatively validate this approximation since these eigenvalues are typically small.

Sf Φ -ŷ 2 -Sf Φ * -ŷ 2 ≤ C max k (σ r+1 (Σ * k ))
where the σ r+1 (Σ * k ) is the (r + 1)-th singular value of Σ * k sorted by decreasing order.

The proof is detailed in Appendix D.

Ideally, we would like to obtain a similar bound for Σ * k -Xk XT k F . We conjecture that a RIP (Restricted Isometry Property) would be needed for such a result. As the verification of RIP remains an open theoretical question in the zero-mean GMM case, we leave this theoretical question for further work.

Related work.

There exist other methods to learn GMMs by incorporating a dimensionality reduction technique. In the mixture models of probabilistic PCAs (MPPCA) [START_REF] Tipping | Mixtures of probabilistic principal component analyzers[END_REF],

the covariance matrix Σ k is parameterized as

(5.10) Σ k = X k X k + γ 2 k I.
The authors use the EM algorithm to optimize the parameters γ k and X k . In the high dimensional data clustering (HDDC) model [START_REF] Bouveyron | High-dimensional data clustering[END_REF], the authors generalized the MPPCA by setting the covariance matrix as

(5.11) Σ k = X k diag(λ k )X T k + γ 2 k I, λ k > 0.
As for the MPPCA, the parameters are learned by the EM algorithm. In the HDMI model [START_REF] Houdard | High-dimensional mixture models for unsupervised image denoising (HDMI)[END_REF],

a model selection algorithm for the intrinsic dimension of each mixture component is proposed.

In this model, the noise variance γ k is a priori fixed and the other parameters are optimized by the EM algorithm. Finally, in the PCA-GMM model [START_REF] Hertrich | PCA reduced gaussian mixture models with applications in superresolution[END_REF], the covariance matrix Σ k is expressed as:

(5.12)

Σ k = 1 γ 2 k (I -X k X T k ) + X k S-1 X T k -1
, S ∈ SP D(P ).

This model is a more general one than HDDC. The parameter γ k can either be fixed a priori or optimized simultaneously with the other parameters by the EM algorithm. Unlike the above models, our model doesn't estimate γ k , we rather set a similar user-defined parameter (called µ in our case) at the denoising step. In our model, we also assume that the intrinsic dimension of each Gaussian component is the same and a priori fixed. Using automatically estimated ranks for covariances is a possible future work.

5.

3. An algorithm for learning patch prior from a sketch : LR-COMP (Low-Rank Continuous Orthogonal Matching Pursuit). Problem (4.7) can be solved approximately using the greedy Compressive Learning OMP called CL-OMP and a variation of CL-OMP called CL-OMP with Replacement (CL-OMPR) [START_REF] Keriven | Sketching for large-scale learning of mixture models, Information and Inference[END_REF][START_REF] Keriven | Blind source separation using mixtures of alpha-stable distributions[END_REF]. These algorithms are based on the Matching Pursuit [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF], Orthogonal Matching Pursuit [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] and Orthogonal Matching Pursuit with Replacement [START_REF] Jain | Orthogonal matching pursuit with replacement[END_REF] for classical compressive sensing, which handle sparse approximation problems.

It starts from an empty support and it expands the support by greedily adding new atoms to the current support. Each new atom θ is found by maximizing the correlation Sf θ , r where r is the current residual. Then it updates the weights and reduces the cost function with a descent algorithm initialized with the current parameters. For better feasible recovery, the algorithm using the replacement method was proposed. The approach increases the number of iterations of CL-OMP and extends the size of support more than the desired sparsity. Then it deletes the extra atoms by using the hard thresholding operator. We adapt these algorithms in the GMMs context with our low-rank approximation. Several modifications are detailed below:

• No Replacement. Although the algorithms using the replacement method show better results on synthetic data, our results tested on image patches show that the replacement step has a negligible effect. Therefore, we run our algorithm without this hard thresholding operator to decrease the computation time.

• Estimation of the factors of covariance instead of the covariance matrices.

As we approximate the covariance matrices with their factors, in each step of the algorithm, we do operations directly on the factorized rank reduced covariance X instead of the covariance matrix Σ to lighten the computations.

• Non-negativity and the normalization of the weights. In Step 1 of the algorithm, we compute the real part of the correlation between the normalized atom and the residual as done in CL-OMP(R). This avoids a negative correlation and negative weights in practice. No matter how Step 3 was computed, using the projected gradient descent or the gradient descent, or with a direct calculation, there's negligible difference in the result and the running time. The weights are not forced to be sumto-one at each iteration. However, after transforming the negative weights to zero, an l 1 -normalization of the weights is performed at the end of the algorithm.

Algorithm 5.1 LR-COMP: Compressive GMM estimation with low-rank covariances [START_REF] Shi | [END_REF].

Input Empirical sketch ŷ, sketching operator S, sparsity K, rank r r ← ŷ; X ← ∅ for t = 1 to K do

Step 1: Perform a descent initialized with a P ×r matrix of normally distributed random numbers:

X * k ← arg max X k Re Sf X k Sf X k 2 , r 2 , init = rand
Step 2: Extend the support:

X ← X ∪ {X * k } Step 3: Find weights: α ← arg min α ŷ - |X| k=1 α k Sf X k 2 2
Step 4: Perform a descent initialized with current parameters:

X, α ← arg min X,α ŷ - |X| k=1 α k Sf X k 2 2 , init = (X, α) Step 5: Update residual: r ← ŷ - |X| k=1 α k Sf X k ; end for Normalize the weights α k such that k α k = 1 return Support X, weights α
The proposed algorithm is summarized in Algorithm 5.1. In practice, we perform Step 4 with a descent algorithm (L-BFGS). We use more iteration in the ultimate Step 4 (for t = K) as a "final adjustment". With this "final adjustment" step, we could reduce the running time by using fewer iterations in Step 4 for t < K. Our algorithm was implemented by extending the MATLAB toolbox [START_REF] Keriven | SketchMLbox -A MATLAB toolbox for large-scale mixture learning[END_REF]. The Matlab implementation of our approach is available at [START_REF] Shi | [END_REF]. The main tool for the implementation of Algorithm 5.1 is to compute the necessary gradients for the optimization problems in Steps 1, 3, and 4.

For the following section, denote the vector v(X) = Sf X ∈ R m .

Expression of the gradient for

Step 1. In step 1, we have the optimization problem (5.13)

X * k ∈ arg max X k ∈R P ×r Re Sf X k Sf X k 2 , r 2 , r ∈ C m . Let F (X k ) = -Re Sf X k Sf X k 2 , r 2 = -v(X k ) T Re(r) v(X k ) 2 , then problem (5.13) turns to (5.14) X k * ∈ arg min X k ∈R P ×r F (X k ).
In practice, with W = [ω 1 , ..., ω m ] ∈ R P ×m the frequency matrix, we compute the gradient of F (X k ) with the following operation:

(

G = - 1 v(X k ) 2 W W T X k * v (X k ) * F (X k ) v (X k ) v (X k ) 2 -Re (r) . 5.15)
Here the symbol * represents the multiplication element by element in MATLAB. A matrix of size m×r multiplied using dot * with a m×1 vector leads to a matrix of size m×r (multiplying all columns of the left side by the same column vector of the right side). The result G is a matrix of size P × r. We need to reshape all the elements of the matrix G into a single column vector, whose result is the gradient ∇ X k F (X k ). The detailed computation is in Appendix E.

The solution of

Step 3. The problem is (5.16)

α * = arg min α∈R |X| y - |X| k=1 α k Sf X k 2 2 y ∈ C m . Denote V (X) = [v(X 1 ), ..., v(X |X| )] ∈ R m×|X| , α = [α, ..., α |X| ] T ∈ R |X| , then the problem can
be expressed as a least-squares minimization (5.17)

α * = arg min α∈R |X| g(α) = arg min α∈R |X| y -V α 2 2 .
We thus have

(5.18) α * = (V T V ) -1 V T ŷ.

Expression of the gradient for

Step 4. The problem is

(5.19) (X * , α) ∈ arg min X∈R | X|,X k ∈R P +P ×r α∈R | X| ŷ - |X[ k=1 α k Sf X k 2 2 . Denote V = [v(X 1 ), ..., v(X | X|)], α = [α 1 , ...α K ] T , we express (5.20) h(X, α) = ŷ -V α 2 2 ,
so we have the gradients

(5.21) ∇ α h(X, α) = 2V T (V α -ŷ) and (5.22) ∇ X k h(X, α) = 2α k ∇ X k v(X k ) T (V α -ŷ).
In practice, as in Step 1, we compute the second gradient by calculating the matrix

(5.23) G 2 = -2α k W (W T X k * v(X k ) * (V α -y)).
The gradient ∇ X k h(X, α) corresponds the vector after reshaping G 2 .

As the function minimized here is smooth; the descent will be guaranteed to converge to a local minimum. Recent works suggest that if all the OMP steps fall close enough to the Gaussian of the global optimum [START_REF] Traonmilin | The basins of attraction of the global minimizers of nonconvex inverse problems with low-dimensional models in infinite dimension[END_REF], this step will converge to the global optimum under a restricted isometry condition.

Complexity of LR-OMP.

When estimating a K-components GMM, the proposed algorithm LR-OMP has a computational cost of the order of O(mP 2 rK 2 ). In each iteration, the computational cost is dominated by the matrix-vector product W (W T X) where W is a matrix of size P ×m and W T X is a matrix of size m×r. As m n, the computational cost of our algorithm is lower than that of the EM. Moreover, it is possible to exploit the advantages of GPU computing, the matrix multiplication can be performed by using multiple GPUs in parallel [START_REF] Zhang | Matrix multiplication on high-density multi-gpu architectures: Theoretical and experimental investigations[END_REF]. This could result in a speed-up, especially for the "final adjustment" step.

5.5. Denoising with low-rank covariance matrices. In this section, we describe some modifications required in EPLL to use our estimated model. The estimated parameters are Φ = { X1 , ..., XK , α1 , ..., αK } with Xk ∈ R P ×r and α k ∈ R + . A singular value decomposition of Xk is given by Xk = Ûk Ŝk V T k . Ûk , Vk ∈ R P ×P are orthogonal matrices and

Ŝk = diag(ŝ k 1 , • • • , ŝkr , 0, • • • , 0) ∈ R P ×P is a diagonal matrix. The r-rank covariance matrix can be expressed with Σkr = Xk XT k = Ûk Ŝ2 k Û T k . We approximate the covariance matrix Σ k with Σ k Σk = Ûk Λk Û T k
where Λk is formed as:

(5.24) Λk =           ŝ2 k 1 . . . ŝ2 kr 0 0 µ . . . µ          
where µ is a user parameter. Denoting Û r k ∈ R P ×r the matrix formed by the first r columns of Ûk and Λr k the matrix formed with the first r rows and r columns of Λk , we have:

(5.25) Σ k + 1 β I P -1 = Û r k ( Λr k + 1 β I r ) -1 Û rT k + β βµ + 1 (I p -Û r k Û rT k ) and (5.26) Σ k + 1 β I P -1 Σ k = Û r k ( Λr k + 1 β I r ) -1 Λr k Û rT k + βµ βµ + 1 (I p -Û r k Û rT k ).
The detailed computation of (5.25) is in Appendix F. Then the Gaussian selection step of EPLL (2.11) becomes

k * i = arg min 1≤k≤K -2 log α k + r j=1 log(ŝ 2 k j + 1 β ) + [v (k) i ] 2 j ŝ2 k j + 1 β - β βµ + 1 [v (k) i ] 2 j (5.27)
where

(5.28) v(k) i = Û rT k zi .
With the optimal component k * i , the estimated patch (2.10) becomes (recall that zi are centered patches and that m i are the estimated mean of patches ẑi = (Σ

k * i + 1 β I P ) -1 Σ k * i zi = Ûk r i * ( Λr k * i + 1 β I r ) -1 Λr k * i Û rT k * i zi + βµ βµ + 1 (I p -Û r k * i Û rT k * i )z i = Û r k * i Λ k * i v(k * i ) i + βµ βµ + 1 (z i -Ûk * r v(k * i ) i )+m i (5.29) with (5.30) Λ k * i = ( Λr k * + 1 β I r ) -1 Λr k * = diag   ŝ2 k * ij ŝ2 k * ij + 1 β   j=1,...,r .
6. Experimental Results. In this section we present several numerical experiments to illustrate the benefits of our approach.

We randomly extract n = 4 × 10 6 patches of size P = 7 × 7 from the training images of the Berkeley Segmentation Database (BSDS) [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF]. Then the patches are compressed into a sketch.

Based on observations from numerical simulations, the scale parameter η must be adjusted for each task and dataset [START_REF] Schellekens | When compressive learning fails: blame the decoder or the sketch?[END_REF]. In [START_REF] Keriven | Sketching for large-scale learning of mixture models, Information and Inference[END_REF], the authors propose to estimate this parameter with a small sketch on a small subset from the dataset. In our work, we choose the optimal parameter η by hand. We then learn a mixture model of K = 20 Gaussian components with low-rank covariance matrices. We compare the denoised results with the results obtained with a GMM (full-rank) prior model learned by the EM algorithm. For the comparison, we train the prior from the same image patches dataset. The denoising is performed with EPLL1 . To evaluate the quality of denoised images, we use two measures: PSNR (Peak Signal to Noise Ratio) and SSIM (Structural Similarity) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]. For the test images, we use two datasets: Set12 [START_REF] Zhang | Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising[END_REF] and BSD68 [START_REF] Roth | Fields of experts: A framework for learning image priors[END_REF] for a thorough evaluation. The code is available at [START_REF] Shi | [END_REF] to reproduce the results below.

Figure 4 shows the denoising performance on 6 images of the Set12 dataset. The noisy images are obtained by adding zero-mean Gaussian noise with standard deviations σ = 20 to the test images. The covariance matrices of the model learned by the sketching have the rank r = 20. We observe that for most of images, we obtain similar or better values of PSNR and SSIM.

Another evaluation was carried out on the images from the BSD68 dataset. The test images have been corrupted by adding white Gaussian noise with standard deviations σ = 15, 60. Table 1 shows the average PSNR and SSIM values on the dataset. On average, our approach results are 0.2dB below the results with EM in terms of PSNR. However, our approach is about 2 times faster than the EM. Moreover, a loss of 0.2dB does not affect the visual quality in most natural images. We also evaluate the similarity of the models learned via EM and LR-OMP. In Figure 5 we visualize the leading eigenvectors of the learned covariance whose weight is the largest. The represented eigenvectors are ordered decreasingly with respect to the eigenvalues. The figure shows that the learned components have rich and similar structures except for the smallest eigenvalues where we observe differences. As we also observe that the eigenvalues decay much faster with our method than with EM, it is hard to interpret further the difference with the result of EM. We can still say that these different "dictionaries" lead to similar denoising results.

6.1. Influence of realization of sketching operator. Our approach performs stable performances with different initialization. Table 2 shows the variability of the PSNR/SSIM over different random sketch realizations. The evaluation is carried out on the classical images: cameraman, house, etc. The noisy images are obtained by adding white Gaussian noise with standard deviations σ = 20 to the test images.

6.2. Influence of sketch size and the compression rate. Theoretically, we can successfully estimate a GMM with sufficiently separated mean and random Fourier sketching with high probability as long as the sketch size m ≥ O(K 2 P polylog(K, P )). In our case, we learn zero-mean Gaussians. From [START_REF] Keriven | Sketching for large-scale learning of mixture models, Information and Inference[END_REF], empirical results indicate that a sketch size of the order of the number of parameters is sufficient (i.e. it is conjectured that K 2 P could be reduced to 6 shows the denoising performance and estimation time with models learned by using different sketch sizes (c = 1, 5, 10, 20). Our experiments show that a larger sketch size doesn't improve the denoising performance necessarily, and indeed it causes more learning time. (D.8)

(∇ X k v(X k )) T r v(X k ) 2 - v(X k ) T r(∇ X k v(X k )) T v(X k ) v(X k ) 2 = - (∇ X k v(X k )) T v(X k ) 2 r + v(X k ) T rv(X k ) v(X k ) 2 2 = (∇ X k v(X k )) T v(X k ) 2 F (X k )v(X k ) v(X k ) 2 -r .
(E.1)

For each component v l (X k ) = e -1 2 ω T l X k X k T ω l , we have (E.2) ∂v l (X k ) ∂X k = -v l (X k )X k T ω l ω T l .
Then for a given vector z ∈ R m

(E.3) ∇ X k v(X k ), z = - m l=1 z l v l (X k )X k T ω l ω T l .
In practice, we compute the scalar product with ∇ X k v(X k ), z = -W (W T X k * (v(X k ) * z)) (E. [START_REF] Burer | Local minima and convergence in low-rank semidefinite programming[END_REF] where W = [ω 1 , ..., ω m ] ∈ M P,m (R) the frequency matrix and * the multiplication element by element in MATLAB. We compute the matrix G as

(E.5) G = - 1 v(X k ) 2 W W T X k * v(X k ) * F (X k )v(X k ) v(X k ) 2 -r .
As a consequence, we reshape the matrix G to a vector which corresponds the gradient

∇ X k F (X k ).
Appendix F. The expression of (5.25). Denoting Û r k ∈ R P ×r (resp. Û c k ∈ R P ×(P -r) )the matrix formed by the first r (resp. the last P -r) columns of Ûk and Λr k the matrix formed with the first r rows and r columns of Λk . We use the bloc matrix multiplication:

1 .

 1 Introduction. We consider the classical noisy observation model of a clean natural image u ∈ R N (composed of N pixels): (1.1) v = u + w where v is the observed degraded version of u. The acquisition noise w is usually assumed to be an additive white Gaussian noise of variance σ, i.e. w i.i.d.

Figure 1 .

 1 Figure 1. A summary of our method

  i of Λ k are the singular values of Σ k . Then we can compute (2.9) by:

4. 1 .Figure 2 .

 12 Figure 2. Schema of sketching

  i.d. Given a training set of n centered patches χ = {x 1 , ..., x n } ⊂ R P , we define the empirical characteristic function with

ηR) 2 Figure 3 .

 23 Figure 3. Curve of the radius distribution density

Proposition 5 . 1 . 2 m l=1 ω l 4 2.

 5124 Let Φ * = {Σ * 1 , ..., Σ * K , α * 1 , ..., α * K } be a minimizer of (5.5). Suppose that there exists a minimizer Φ = { X1 , ..., XK , α1 , ..., αK } for the problem (5.8). Let C = 1 Then we have:

Figure 4 .

 4 Figure 4. From left to right: Original images, noisy images with noise σ = 20, results with EM model, results with LR-COMP model. The denoising results are evaluated with PSNR/SSIM. Similar denoising performances are obtained with LR-COMP with a 1000 times smaller compressed database. To estimate the prior model, our method is 2 times faster than the EM algorithm.

Figure 5 .

 5 Figure 5. The first 20 eigenvectors of the covariance matrices (for the heaviest weight) learned by LR-OMP with rank r = 20 (a), r = 49(b) and EM (c). The decay of the corresponding eigenvalues (d).

Figure 6 .

 6 Figure 6. Denoising performance (PSNR/SSIM) and estimation time (hours) of models learned with different sketch size. c =1 (a), 5 (b), 10 (c), 20 (d).

6. 3 .

 3 Influence of the rank r.

Figure 7

 7 shows the denoising performance of models estimated with different ranks. Our experiments show that the model with reduced rank results in a minor PSNR/SSIM drop compared to the full-rank model. However, the learning time is much faster. According to the experiments, we cannot reduce the rank further (less than 20) to keep good denoising performance.

Figure 7 .

 7 Figure 7. Denoising performance (PSNR/SSIM) of models learned with different intrinsic dimensions. r = 10 (left), 20 (middle), 49(right).

6. 4 .ω l -ŷl 2 . 2 = 2 ≤ 2 . 5 ) 2 2 2 2

 42222522 Learning time. In terms of time complexity, the running time depends on the number of components K and the complexity of the descent algorithm. In our approach, we use the Limited-memory BFGS algorithm to handle the optimization problems in Step 1 and 4. The latter is the most time-consuming part of the algorithm. To get the model (c = 10, r = 20) that achieves the denoising performance of our experiments (Figure4), it takes less than 2 hours on a computer with 2 * 32 cores AMD EPYC 7452 @ 2,35 GHz. With the same environment, our learning algorithm is about 2 times faster than the EM algorithm (3.68h)2 . and suppose that there exists a minimizer Φk = ( Xk , αk ) for the problem (5Let Σk be the best rank-r approximation of Σ * k with the rank r i.e.(D.3) Σk ∈ arg min Σ,rank(Σ)=r Σ * k -Σ 2 F .Define Φ = ( Σk , α * k ). According to the definition (D.2) and the triangle inequality, we haveSf Φ -y 2 ≤ Sf Φ -ŷ Sf Φ -Sf Φ * + Sf Φ * -ŷ Sf Φ -Sf Φ * 2 + Sf Φ * -yUsing the convexity inequality |1 -e -x | ≤ |x| and Cauchy-Schwarz inequality, we havee -1 2 ω T l Σk ω l (1 -e -1 2 ω T l (Σ * k -Σk )ω l ) ≤ 1 -e -(Σ * k -Σk )ω l . (D.6)By the Eckart-Young-Mirsky theorem, we have that the largest singular value of Σ * k -Σk isσ r+1 (Σ * k ) and (D.7) ω l , (Σ * k -Σk )ω l ≤ ω l (Σ * k -Σk ) op = ω l σ r+1 (Σ * k ).Therefore, using K k=1 α k = α 1 = 1 and Hölder's inequality, and upper bound on (D.

Table 1

 1 The average PSNR and SSIM on the BSD68 dataset with 2 different levels of noise.

	σ	Sketching	EM
	15 31.8 / .876 32.0 / .879
	50 24.4 / .637 24.6 / .646

Table 2

 2 Image denoising performance comparison of models estimated over different random sketch realization.KP ). In our experiments, we set m = cK(P × r + 1) = 10K(P × r + 1) ≈ 2 × 10 5 , i.e the compressed database is approximately 1000 times smaller than the original patch database. The gains in terms of memory is approximately n m times compared to the EM approach. Figure

		Realization 1 Realization 2 Realization 3 Realization 4
	cameraman	33.1 / .930	33.3 / .930	33.4 / .930	33.4 / .930
	house	35.4 / .926	35.4 / .925	35.4 / .926	35.4 / .925
	jetplane	32.1 / .936	32.3 / .937	32.4 / .937	32.4 / .937
	lena	32.0/ .931	32.2 / .931	32.2 / .931	32.2 / .931
	pirate	29.8 / .907	30.0 / .908	30.0 / .908	30.0 / .908

Matlab implementation based on the code of[START_REF] Parameswaran | Accelerating gmm-based patch priors for image restoration: Three ingredients for a 100× speed-up[END_REF].

Mo Chen (2021). EM Algorithm for Gaussian Mixture Model (EM GMM) (https://www.mathworks.com/matlabcentral/fileexchange/26184-em-algorithm-for-gaussian-mixture-modelem-gmm), MATLAB Central File Exchange. Retrieved October 11, 2021.

Acknowledgments. The authors acknowledge the support of the French National Research Agency (ANR) under reference ANR-20-CE40-0001 (EFFIREG project). Experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr/). The authors are indebted to anonymous reviewers for providing insightful comments which have resulted in this paper.

This work was funded by the French National Research Agency (ANR) under reference ANR-20-CE40-0001 (EFFIREG project).

Conclusions.

In this work, we adapt the sketching framework in the context of image patches. We propose an algorithm LR-COMP to estimate a GMM with low-rank approximation and provide an implementation of the algorithm. Experiments illustrate that a highdimensional GMM can be learned from a compressed database and then used for patch-based denoising. We achieve denoising performances close to state-of-the art model based methods while the learning procedure uses less memory and time than the classical EM algorithm.

In future works, we can generalize our approach to other models such as GGMM (Generalized Gaussian Mixture Model) for a better denoising performance [START_REF] Deledalle | Image denoising with generalized gaussian mixture model patch priors[END_REF]. We also aim to adapt the sketching to more inverse problems such as image super-resolution, image deblurring, etc.

Another perspective is to extend our model to the study of video denoising method as the potential of the technique for video restoration remains unexplored. As mentioned earlier, the scale parameter η must be adjusted for each task and dataset. In our work, we choose this parameter by hand. Moreover, in our model, the intrinsic dimension of each Gaussian component is assumed the same and a prior fixed. It could be useful to design a procedure to estimate the hyper-parameters automatically in future work. In our work, we estimate a GMM with zero-mean. In this context, the notion of separation used to prove the restricted isometry property which in turn proves identifiability of the GMM in [START_REF] Gribonval | Statistical learning guarantees for compressive clustering and compressive mixture modeling[END_REF] and convergence of gradient descent in [START_REF] Traonmilin | The basins of attraction of the global minimizers of nonconvex inverse problems with low-dimensional models in infinite dimension[END_REF]. Proving a RIP on zero-mean would require a new notion of separation. We conjecture that an angular separation between Gaussian might enable us to prove such RIP.

Such separation could e.g. compare the angle between eigenvectors of covariances by decreasing eigenvalue amplitude and weight the separation accordingly. The low-rank model should add even more separation as the Gaussian are supported on different sub-spaces. We still show empirically that the sketching process is successful without this separation assumption. This opens interesting new theoretical questions for the study of the success of compressive learning in patch-based image processing.

Definition A.2. Frobenius norm. For a matrix A ∈ C m×n , the Frobenius norm of A is defined as

where σ i (A) are the singular values of A.

Definition A.3. Operator norm. For a continuous linear operator A : V → W , the operator norm of A is defined as

the singular value decomposition of D with Σ = diag(σ 1 , ..., σ m ). Let U r (resp. V r ) be the matrix formed by the first r columns of U (resp. V ) and Σ r = diag(σ 1 , ..., σ r ). Then the r-rank matrix, obtained from the truncated singular value decomposition:

the minimizer of the low-rank approximation:

The minimizer D * is unique if and only if σ r+1 < σ r .

Appendix B. EM algorithm.

Given a data set of n clean training patches χ = {x 1 , ..., x n } ⊂ R P ×n , the EM algorithm for estimating a GMM can be summarized as follows:

1. Define the number of components K. For each component k, we initialize the parameters Θ k = (µ k , Σ k , α k ) randomly, and we compute the log likelihood

2

. E-Step

Compute the posterior function Γ i,k with the current parameters Θ k :

M-Step

Re-estimate the parameters Θ new k with the Γ i,k obtained in the E-Step:

4. Re-evaluate the log likelihood. Iterate E-Step and M-Step until the log likelihood or the parameters are not changing much.

Appendix C. Design of the adapted radius distribution. The choice of the probability distribution ∆ to draw the frequencies is a key ingredient in designing the sketching operator. In our work, we choose the frequency distribution called the adapted radius [START_REF] Keriven | Compressive k-means[END_REF] with a heuristic. In this appendix, we describe the design of the adapted radius distribution.

Assuming that we want to estimate a P -dimensional Gaussian g = N (0, I P ), we can compute the characteristic function ψ g (ω) associated with g:

The adapted radius heuristic proposes not to sample ω directly but rather to sample the radius of the P -dimensional Gaussian R = √ ω T ω. Thus, we draw the frequency ω ∈ R P as

where the radius R ∈ R + is chosen with a radius distribution R ∼ p R (R; η), and the direction ϕ ∈ R P is uniformly generated on the l 2 unit sphere S P -1 , i.e. ϕ ∼ U(S P -1 ). Then, the characteristic function ψ g (ω) reduces to

We obtain a one-dimensional Gaussian function for R. To design the radius distribution, we consider the estimation of a Gaussian g = N (0, 1). We aim at sampling the radius R leading to large variations of the characteristic function when the parameters are closed to the true parameters. In other words, when parameters (µ, σ 2 ) are closed to (0,1), we want have a large

This can be accomplished by promoting the radius R which makes the norm of the gradient ∇ψ (µ,σ 2 ) (R) 2 large. Recall that ψ (µ,σ 2 ) (R) = e -iµR e -1 2 σ 2 R 2 and the norm of the gradient is:

It yields the density of a radius distribution :

Here the scale parameter η should be chosen to probe the spectral content of the true GMM model. Appendix E. Calculation of the gradient. Denote F (X k ) = -v(X k ) T r v(X k ) 2 , where r ∈ R m is the real part of r. We compute the gradient of F as follows: