domain: stat math

text:
Jérémie Bigot

Charles Deledalle

Delphine Féral

Generalized SURE for optimal shrinkage of singular values in low-rank matrix denoising

Keywords: matrix denoising, singular value decomposition, low-rank model, Gaussian spiked population model, spectral estimator, Stein's unbiased risk estimate, random matrix theory, exponential family, optimal shrinkage rule, degrees of freedom AMS classifications: 62H12, 62H25

We consider the problem of estimating a low-rank signal matrix from noisy measurements under the assumption that the distribution of the data matrix belongs to an exponential family. In this setting, we derive generalized Stein's unbiased risk estimation (SURE) formulas that hold for any spectral estimators which shrink or threshold the singular values of the data matrix. This leads to new data-driven spectral estimators, whose optimality is discussed using tools from random matrix theory and through numerical experiments. Under the spiked population model and in the asymptotic setting where the dimensions of the data matrix are let going to infinity, some theoretical properties of our approach are compared to recent results on asymptotically optimal shrinking rules for Gaussian noise. It also leads to new procedures for singular values shrinkage in finite-dimensional matrix denoising for Gamma-distributed and Poisson-distributed measurements.

Introduction

Low rank matrix denoising in an exponential family

In various applications, it is of interest to estimate a signal matrix from noisy data. Typical examples include the case of data that are produced in a matrix form, while others are concerned with observations from multiple samples that can be organized in a matrix form. In such setting, a typical inference problem involves the estimation of an unknown (non-random) signal matrix X ∈ R n×m from a noisy data matrix Y satisfying the model:

Y = X + W , (1.1)
where W is an n × m noise matrix with real entries W ij assumed to be independent random variables with E[W ij ] = 0 and Var(W ij ) = τ 2 ij for 1 ≤ i ≤ n and 1 ≤ j ≤ m. In this paper, we focus on the situation where the signal matrix X is assumed to have a low rank structure, and we consider the general setting where the distribution of Y belongs to a continuous exponential family parametrized by the entries of the matrix X = E[Y ]. For discrete observations (count data), we also consider the specific case of Poisson noise.

The low rank assumption on X is often met in practice when there exists a significant correlation between the columns of X. This can be the case when the columns of X represent 2D images at different wavelength of hyperspectral data, since images at nearby wavelengths are strongly correlated [START_REF] Candès | Unbiased risk estimates for singular value thresholding and spectral estimators[END_REF]. Further applications, where low-rank modeling of X is relevant, can be found in genomics [START_REF] Wall | Svdman-singular value decomposition analysis of microarray data[END_REF][START_REF] Alter | Singular value decomposition for genomewide expression data processing and modeling[END_REF], NMR spectroscopy [START_REF] Nguyen | Spatiotemporal denoising of mr spectroscopic imaging data by low-rank approximations[END_REF], collaborative filtering [START_REF] Candès | Exact matrix completion via convex optimization[END_REF] or medical imaging [BD06, LBH + 12], among many others.

Low-rank matrix estimation is classically done in the setting where the additive noise is Gaussian with homoscedastic variance. The more general case of observations sampled from an exponential family is less developed, but there exists an increasing research interest in the study of low rank matrix recovery beyond the Gaussian case. Examples of low-rank matrix recovering from Poisson distributed observations can be found in applications with count data such as network traffic analysis [START_REF] Andres Bazerque | Inference of Poisson count processes using low-rank tensor data[END_REF] or call center data [START_REF] Shen | Analysis of call centre arrival data using singular value decomposition: Research articles[END_REF]. A theory for low-rank matrix recovery and completion in the case of Poisson observations has also been recently proposed in [START_REF] Cao | Poisson matrix recovery and completion[END_REF]. Matrix completion under a low rank assumption with additive errors having a subexponential distribution and belonging to an exponential family has also been considered in [START_REF] Lafond | Low rank matrix completion with exponential family noise[END_REF]. The recent work [START_REF] Udell | Generalized low rank models[END_REF] proposes a novel framework to approximate, by a low rank matrix, a tabular data set made of numerical, Boolean, categorical or ordinal observations.

The class of spectral estimators

A standard approach to estimate a low rank matrix relies on the singular value decomposition (SVD) of the data matrix

Y = min(n,m) k=1 σk ũk ṽt k , (1.2)
where σ1 ≥ σ2 ≥ . . . ≥ σmin(n,m) ≥ 0 denote its singular values, and ũk , ṽk denote the associated singular vectors. In this paper, we propose to consider the class of spectral estimators Xf = f (Y ), where f : R n×m → R n×m is a (possibly data-dependent) mapping that acts on the singular values of the data matrix Y while leaving its singular vectors unchanged. More precisely, these estimators take the form

Xf = f (Y ) = min(n,m) k=1 f k (Y )ũ k ṽt k , (1.3)
where, for each 1 ≤ k ≤ min(n, m), f k (Y ) are real positive values that may depend only on σk (hence we write f k (σ k )) or on the whole matrix Y .

Investigated spectral estimators

Typical examples of spectral estimators include the classical principal component analysis (PCA) applied to matrix denoising defined, for some 1 ≤ r ≤ min(n, m), as

Xr = r k=1 σk ũk ṽt k with σk = f k (σ k ) = σk (1.4)
for all 1 ≤ k ≤ r and where it is implicitely understood that f k (σ k ) = 0 for k ≥ r + 1.

Another typical spectral estimator in matrix denoising with Gaussian measurements is the softthresholding [START_REF] Candès | Unbiased risk estimates for singular value thresholding and spectral estimators[END_REF] which corresponds to the choice

Xsoft = min(m,n) k=1 σk ũk ṽt k with σk = f k (Y ) = 1 - λ(Y ) σk + σk , (1.5)
for all 1 ≤ k ≤ min(n, m) and where λ(Y ) > 0 is a possibly data-dependent threshold parameter, and (x) + = max(x, 0) for any x ∈ R. Finaly, we will consider a more general class of shrinkage estimators, encompassing the PCA and the soft-thresholding, that perform

Xw = min(m,n) k=1 σk ũk ṽt k with σk = f k (Y ) = w k (Y )σ k , (1.6)
where w k (Y ) ∈ [0, 1] is a possibly data-dependent shrinking weight.

Main contributions

Under the assumption that the distribution of Y belongs to an exponential family, the goal of this paper is to derive data-driven choices for the weights w k (Y ) in (1.3).We construct estimators via a two-step procedure. First, an active set of non-zero singular values is defined. Then, in a second step, weights w k (Y ) associated with non-zero singular values are optimized, and shown to reach desired asymptotical properties in the Gaussian spiked population model. The main contributions of the paper are then the following ones.

An AIC inspired criterion for rank and singular values locations estimation

When no a priori is available on the rank of the signal matrix X, optimizing for the weights w k , for all 1 ≤ k ≤ min(m, n), can lead to estimators with large variance (i.e., overfitting the noise).

We propose an automatic rule to prelocalize the subset of non-zero singular values. An active set s ⊆ I = {1, 2, . . . , min(n, m)} of singular values is defined as the minimizer of a penalized log-likelihood criterion that is inspired by the Akaike information criterion (AIC)

s * ∈ arg min s⊆I -2 log q(Y ; Xs ) + 2|s|p n,m with p n,m = 1 2 √ m + √ n 2 , (1.7)
where Xs = k∈s σk ũk ṽt k , |s| is the cardinal of s, and q(Y ; Xs ) is the likelihood of the data in a given exponential family with estimated parameter Xs . For the case of Gaussian measurements with homoscedastic variance τ 2 , one has that q(Y ;

Xs ) = Y -Xs 2 F /2τ 2 , where • F denotes the Frobenius norm of a matrix, and we show that the active set of singular values boils down to

s = {k ; σk > c n,m + }, (1.8)
where c n,m

+ = τ ( √ m + √ n).
For Gamma and Poisson measurements, we resort to a greedy optimization procedure described in Section 4.

Once the active set has been determined, the subsequent shrinkage estimator is obtained by optimizing only for the weights within this subset while setting the other ones to zero.

Novel data-driven shrinkage rules minimizing SURE-like formulas

We use the principle of Stein's unbiased risk estimation (SURE) [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF] to derive unbiased estimation formulas for the mean squared error (MSE) risk and mean Kullback-Leibler (MKL) risks of spectral estimators. Minimizing such SURE-like formulas over an appropriate class of spectral estimators is shown to lead to novel data-driven shrinkage rules of the singular values of the matrix Y . In particular, our approach leads to novel spectral estimators in situations where the variances τ 2 ij of the entries W ij of the noise matrix are not necessarily equal, and may depend on the signal matrix X.

As an illustrative example, let us consider spectral estimators of the form

X1 w = f (Y ) = w 1 (Y )σ 1 ũ1 ṽt 1 , (1.9)
which only act on the first singular value σ1 of the data while setting all the other ones to zero. In this paper, examples of data-driven choices for the weight w 1 (Y ) are the following ones:

• for Gaussian measurements with n ≤ m and known homoscedastic variance

τ 2 w 1 (Y ) = 1 - τ 2 σ2 1 1 + |m -n| + 2 n =2 σ2 1 σ2 1 -σ2 + 1 1 {σ1>c n,m + } , (1.10)
• for Gamma measurements with τ 2 ij = X 2 ij /L and L > 2 (see Section 2.1 for a precise definition),

w 1 (Y ) = min   1,   L -1 Lmn n i=1 m j=1 X1 ij Y ij + 1 Lmn   1 + |m -n| + 2 min(n,m) =2 σ2 1 σ2 1 -σ2     -1   1 1 {1∈s * } , (1.11)
• for Poisson measurements with τ 2 ij = X ij (see Section 2.1 for a precise definition)

w 1 (Y ) = min   1, n i=1 m j=1 Y ij n i=1 m j=1 X1 ij   1 1 {1∈s * } .
(1.12)

Beyond the case of rank one, closed-form solutions for the weights cannot be obtained, except for the case of Gaussian measurements with homoscedastic variance τ 2 . In this latter case, the rule for w 1 (Y ) in (1.10) generalizes to other eigenvalues w k (Y ) as

w k (Y ) =   1 - τ 2 σ2 k   1 + |m -n| + 2 min(n,m) =1; =k σ2 k σ2 k -σ2     + 1 1 {σk>c n,m + } . (1.13)
For Gamma or Poisson distributed measurements, we propose fast algorithms to get numerical approximations of the weights w k (Y ) (see Section 5.2 for more details).

Asymptotic properties in the Gaussian spiked population model

Another contribution of the paper is to discuss the optimality of the shrinking weights (1.13) for Gaussian noise in the asymptotic setting where the dimensions of the matrix Y are let going to infinity. These theoretical results are obtained for the so-called spiked population model that has been introduced in the literature on random matrix theory and high-dimensional covariance matrix estimation (see e.g. [BS06, BN12, DS07, SN13]). All the theoretical and asymptotic results of the paper (other than derivation of proposed estimators) assume this model.

Definition 1.1. The Gaussian spiked population model corresponds to the following setting:

• the W ij in (1.1) are iid Gaussian random variables with zero mean and variance τ 2 = 1/m,

• the X ij 's in (1.1) are the entries of an unknown n×m matrix X that has a low rank structure, meaning that it admits the SVD X =

r * k=1 σ k u k v t k
, where u k and v k are the left and right singular vectors associated to the singular value σ k > 0, for each 1 ≤ k ≤ r * , with σ 1 > σ 2 > . . . > σ r * ,

• the rank r * of the matrix X is assumed to be fixed,

• the dimensions of the data matrix Y = X + W are let going to infinity in the asymptotic framework where the sequence m = m n ≥ n is such that lim n→+∞ n m = c with 0 < c ≤ 1. In the Gaussian spiked population model, the asymptotic locations of the empirical singular values σ1 ≥ . . . ≥ σmin(n,m) are well understood in the random matrix theory (further details are given in Section 3.1). Note that the setting where the rank r * is not held fixed but allowed to grow with min(n, m) is very different, see e.g. [START_REF] Ledoit | Nonlinear shrinkage estimation of large-dimensional covariance matrices[END_REF] and references therein.

Under the Gaussian spiked population model, our contributions are then as follows:

• we prove the convergence of the SURE formula when the dimensions of Y tend to infinity,

• it is shown that minimizing the asymptotic value of SURE leads to the same estimator as the limiting value of the estimator obtained by minimizing the SURE,

• this model allows to show that the novel data-driven spectral estimators derived in this paper are asymptotically connected to existing optimal shrinkage rules [SN13, GD14a, Nad14] for low-rank matrix denoising,

• in this setting, we are also able to connect the choice of the penalty function 2|s|p n,m in (1.7) with Stein's notion of degrees of freedom (see e.g. [START_REF] Efron | The estimation of prediction error: Covariance penalties and crossvalidation[END_REF]) for spectral estimators.

Numerical experiments and publicly available source code

As the theoretical properties of our estimators are studied in an asymptotic setting, we report the results of various numerical experiments to analyze the performances of the proposed estimators for finite-dimensional matrices. These experiments allow the comparison with existing shrinkage rules for Gaussian-distributed measurements and they are also used to shed some lights on the finite sample properties of the method for Gamma-distributed or Poisson-distributed measurements. We also exhibit the settings where the signal matrix X is either easy or more difficult to recover. From these experiments, the main findings are the following ones:

• the use of an appropriate active set s of singular values is an essential step for the quality of shrinkage estimators whose weights are data-driven by SURE-like estimators; taking s = {1, . . . , min(n, m)} leads to poor results while the choice of s = s * minimizing the AIC criterion (1.7) appears to yield the best performances,

• for Gaussian noise, the performances of our approach are similar to those obtained by the asymptotically optimal spectral estimator proposed in [START_REF] Gavish | Optimal shrinkage of singular values[END_REF] when the true rank r * of the signal matrix X is sufficiently small,

• for Gamma or Poisson distributed measurements, the spectral estimators proposed in this paper give better results than estimators based on PCA (restricted to the active set s * ) or soft-thresholding of singular values.

Beyond the case of Gaussian noise, the implementation of the estimators is not straightforward, and we thus provide publicly available source code at https://www.math.u-bordeaux.fr/~cdeledal/gsure_low_rank to reproduce the figures and the numerical experiments of this paper.

Related results in the literature

Early work on singular value thresholding began with the work in [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF] on the best approximation of fixed rank to the data matrix Y . Spectral estimators with different amounts of shrinkage for each singular value of the data matrix have then been proposed in [START_REF] Efron | Empirical Bayes on Vector Observations: An Extension of Stein's Method[END_REF][START_REF] Efron | Multivariate empirical bayes and estimation of covariance matrices[END_REF]. In the case of Gaussian measurements with homoscedastic variance, the problem of estimating X under a low-rank assumption has recently received a lot of attention in the literature on high-dimensional statistics, see e.g. [CSLT13, DG14, JS15, SN13]. Recent works [START_REF] Gavish | Optimal shrinkage of singular values[END_REF][START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF] also consider the more general setting where the distribution of the additive noise matrix W is orthogonally invariant, and such that its entries are iid random variables with zero mean and finite fourth moment. In all these papers, the authors have focused on spectral estimators which shrink or threshold the singular values of Y , while its singular vectors are left unchanged. In this setting, the main issue is to derive optimal shrinkage rules that depends on the class of spectral estimators that is considered, on the loss function used to measure the risk of an estimator of X, and on appropriate assumptions for the distribution of the additive noise matrix W .

Organization of the paper

Section 2 is devoted to the analysis of a data matrix whose entries are distributed according to a continuous exponential family. SURE-like formula are first given for the mean squared error risk, and then for the Kullback-Leibler risk. As an example of discrete exponential family, we also derive such risk estimators for Poisson distributed measurements. The computation of datadriven shrinkage rules is then discussed for Gaussian, Gamma and Poisson noises. In Section 3, we restrict our attention to the Gaussian spiked population model in order to derive asymptotic properties of our approach. We study the asymptotic behavior of the SURE formula proposed in [START_REF] Candès | Unbiased risk estimates for singular value thresholding and spectral estimators[END_REF][START_REF] Donoho | Minimax risk of matrix denoising by singular value thresholding[END_REF] for spectral estimators using tools from RMT. This result allows to make a connection between data-driven spectral estimators minimizing the SURE for Gaussian noise, and the asymptotically optimal shrinkage rules proposed in [START_REF] Shabalin | Reconstruction of a low-rank matrix in the presence of Gaussian noise[END_REF][START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF] and [START_REF] Gavish | Optimal shrinkage of singular values[END_REF]. In Section 4, we study the penalized log-likelihood criterion (1.7) used to select an active set of singular values. Its connection to the degrees of freedom of spectral estimators and rank estimation in matrix denoising is discussed. Various numerical experiments are finally proposed in Section 5 to illustrate the usefulness of the approach developed in this paper for low-rank denoising and to compare its performances with existing methods. The proofs of the main results of the paper are gathered in a technical Appendix A, and numerical implementation details are described in Appendix B.

SURE-like formulas in exponential families

For an introduction to exponential families, we refer to [START_REF] Brown | Fundamentals of statistical exponential families: with applications in statistical decision theory[END_REF]. The idea of unbiased risk estimation in exponential families dates back to [START_REF] Hudson | A natural identity for exponential families with applications in multiparameter estimation[END_REF]. More recently, generalized SURE formulas have been proposed for the estimation of the MSE risk, for denoising under various continuous and discrete distributions in [START_REF] Raphan | Learning to be Bayesian without supervision[END_REF], and for inverse problems whithin the continuous exponential families in [START_REF] Eldar | Generalized sure for exponential families: Applications to regularization[END_REF]. In [START_REF] Deledalle | Estimation of Kullback-Leibler losses for noisy recovery problems within the exponential family[END_REF], SURE-like formula are derived for the estimation of the Kullback-Leibler risk that applies to both continuous and discrete exponential families. In what follows, we borrow some ideas and results from these works. We first treat the case of continuous exponential families, and then we focus on Poisson data in the discrete case.

Data sampled from a continuous exponential family

We recall that Y is an n×m matrix with independent and real entries Y ij . For each 1 ≤ i ≤ n and 1 ≤ j ≤ m, we assume that the random variable Y ij is sampled from a continuous exponential family, in the sense that each Y ij admits a probability density function (pdf) q(y; X ij ) with respect to the Lebesgue measure dy on the real line Y = R. The pdf q(y; X ij ) of Y ij can thus be written in the general form:

q(y; X ij ) = h(y) exp (η(X ij )y -A(η(X ij ))) , y ∈ Y, (2.1)
where η (the link function) is a one-to-one and smooth function, A (the log-partition function) is a twice differentiable mapping, h is a known function, and X ij is an unknown parameter of interest belonging to some open subset X of R. Throughout the paper, we will suppose that the following assumption holds:

Assumption 2.1. The link function η and the log-partition function A are such that

A (η(x)) = x for all x ∈ X ,
where A denotes the first derivative of A.

Since E[Y ij ] = A (η(X ij )) for exponential families in the general form (2.1), Assumption 2.1 implies that E[Y ij ] = X ij , and thus the data matrix satisfies the relation Y = X + W where W is a centered noise matrix, which is in agreement with model (1.1). Now, if we let Θ = η(X ), it will be also convenient to consider the expression of the pdf of Y ij in the canonical form:

p(y; θ ij ) = h(y) exp (θ ij y -A(θ ij )) , y ∈ Y, (2.2)
where θ ij = η(X ij ) ∈ Θ is usually called the canonical parameter of the exponential family. Finally, we recall the relation

Var(Y ij ) = A (θ ij ) = A (η(X ij ))
where A denotes the second derivative of A. Then, we denote by θ the n × m matrix whose entries are the θ ij 's.

Examples of data satisfying model (2.1) are the following ones:

Gaussian noise with known variance τ 2 :

q(y; X ij ) = 1 √ 2π exp - (y -X ij ) 2 2τ 2 , E[Y ij ] = X ij , Var(Y ij ) = τ 2 , Y = R, X = R, Θ = R, h(y) = 1 √ 2πτ exp - y 2 2τ 2 , η(x) = x τ 2 , A(θ) = τ 2 θ 2 2 .
Gamma-distributed measurements with known shape parameter L > 0:

q(y; X ij ) = L L y L-1 Γ(L)X L ij exp -L y X ij 1 1 ]0,+∞[ (y), E[Y ij ] = X ij , Var(Y ij ) = X 2 ij L , Y = R, X =]0, +∞[, Θ =]-∞, 0[, h(y) = L L y L-1 Γ(L) 1 1 ]0,+∞[ (y), η(x) = - L x , A(θ) = -L log - θ L .
The matrix θ = η(X) can then be estimated via the n × m matrix θf = θf (Y ) whose entries are given by

θf ij (Y ) = η Xf ij , for all 1 ≤ i ≤ n, 1 ≤ j ≤ m, (2.3)
where Xf ij is a spectral estimator as defined in eq. (1.3). In the rest of this section, we follow the arguments in [START_REF] Eldar | Generalized sure for exponential families: Applications to regularization[END_REF] and [START_REF] Deledalle | Estimation of Kullback-Leibler losses for noisy recovery problems within the exponential family[END_REF] to derive SURE-like formulas under the exponential family for the estimators θf and Xf , using either the meansquared error (MSE) risk or the Kullback-Leibler (KL) risk.

Unbiased estimation of the MSE risk

We consider the following MSE risk which provides a measure of discrepancy in the space Θ of natural parameters, and then indirectly in the space of interest X .

Definition 2.1. The squared error (SE) risk of θf is SE( θf , θ) = θfθ 2 F , and the meansquared error (MSE) risk of θf is defined as

MSE( θf , θ) = E SE( θf , θ) = E θf -θ 2 F .
Using the above MSE risk to compare θf and θ implies that the discrepancy between the estimator Xf and the matrix of interest X is measured by the quantity

MSE η ( Xf , X) = MSE(η( Xf ), η(X)) which is different from MSE( Xf , X). For Gaussian noise, MSE η ( Xf , X) = 1 τ 2 E Xf -X 2
F , while for Gamma distributed measurements with known shape parameter L > 0, it follows that

MSE η ( Xf , X) = L 2 n i=1 m j=1   X ij - Xf ij X ij Xf ij   2 .
The following proposition gives a SURE formula for the MSE risk introduced in Definition 2.1.

Proposition 2.1. Suppose that the data are sampled from a continuous exponential family. Assume that the function h, in the definition (2.2) of the exponential family, is twice continuously differentiable on Y = R. If the following condition holds

E θf ij (Y ) < +∞, for all 1 ≤ i ≤ n, 1 ≤ j ≤ m, (2.4)
then, the quantity

GSURE( θf ) = θf (Y ) 2 + n i=1 m j=1 2 h (Y ij ) h(Y ij ) θf ij (Y ) + h (Y ij ) h(Y ij ) + 2 div θf (Y ), (2.5
)

where div θf (Y ) = n i=1 m j=1 ∂ θf ij (Y ) ∂Y ij , is an unbiased estimator of MSE( θf , θ)
Note that GSURE( θf ) is an unbiased estimator of MSE( θf , θ) and not of MSE( Xf , X). It is shown in Section 3.3 that the results of Proposition 2.4 coincide with the approach in [START_REF] Candès | Unbiased risk estimates for singular value thresholding and spectral estimators[END_REF] on the derivation of a SURE formula in the case of Gaussian noise for smooth spectral estimators.

In the case of Gamma noise, assuming L > 2 implies that the conditions on the function h in Proposition 2.1 is satisfied, hence assuming that conditions (2.4) holds as well, and using that

θf ij (Y ) = - L f ij (Y ) and ∂ θf ij (Y ) ∂Y ij = L |f ij (Y )| 2 ∂f ij (Y ) ∂Y ij , it follows that GSURE( θf ) = n i=1 m j=1 L 2 |f ij (Y )| 2 - 2L(L -1) Y ij f ij (Y ) + 2L |f ij (Y )| 2 ∂f ij (Y ) ∂Y ij - (L -1)(L -2) |Y ij | 2 .
(2.6)

Unbiased estimation of KL risks

Following the terminology in [START_REF] Deledalle | Estimation of Kullback-Leibler losses for noisy recovery problems within the exponential family[END_REF], let us now introduce two different notions of Kullback-Leibler risk, which arise from the non-symmetry of this discrepancy measure.

Definition 2.2. Let f : R n×m → R n×m be a smooth spectral function. Consider the estimator θf defined by (2.3), where Y is a matrix whose entries Y ij are independent random variables sampled from the exponential family (2.2) in canonical form:

• the Kullback-Leibler synthesis (KLS) risk of θf is defined as

KLS( θf , θ) = n i=1 m j=1 R log   p(y; θf ij ) p(y; θ ij )   p(y; θf ij ) dy,
and the mean KLS risk of θf is defined as

MKLS( θf , θ) = E KLS( θf , θ) ,
• the Kullback-Leibler analysis (KLA) risk of θf is defined as

KLA( θf , θ) = n i=1 m j=1 R log   p(y; θ ij ) p(y; θf ij )   p(y; θ ij ) dy,
and the mean KLA risk of θf is defined as MKLA( θf , θ) = E KLA( θf , θ) .

A key advantage of the Kullback-Leibler risk is that it measures the discrepancy between the unknown distribution p(y; θ ij ) and its estimate p(y; θf ij ). It is thus invariant with respect to the reparametrization θf = η( Xf ) (unlike the MSE risk), and we may also write MKLS( θf , θ) = MKLS( Xf , X) and MKLA( θf , θ) = MKLA( Xf , X). As suggested in [START_REF] Deledalle | Estimation of Kullback-Leibler losses for noisy recovery problems within the exponential family[END_REF], the MKLA risk represents how well the distribution p(y; θf ij ) explain a random variable Y ij sampled from the pdf p(y; θ ij ). The MKLA risk is a natural loss function in many statistical problems since it takes as a reference measure the true distribution of the data, see e.g. [START_REF] Hall | On Kullback-Leibler loss and density estimation[END_REF]. The MKLS risk represents how well one may generate an independent copy of Y ij by sampling a random variable from the pdf p(y; θf ij ). The MKLS risk has also been considered in various inference problems in statistics [START_REF] Hannig | On Poisson signal estimation under Kullback-Leibler discrepancy and squared risk[END_REF][START_REF] Yanagimoto | The Kullback-Leibler risk of the Stein estimator and the conditional MLE[END_REF].

By simple calculation, it follows that

MKLS( θf , θ) = n i=1 m j=1 E θf ij -θ ij A ( θf ij ) + A(θ ij ) -E A( θf ij ) , (2.7) and MKLA( θf , θ) = n i=1 m j=1 E θ ij - θf ij A (θ ij ) + E A( θf ij ) -A(θ ij ). (2.8)
Hence, in the case of Gaussian measurements with known variance τ 2 , we easily retrieve that

MKLS( θf , θ) = MKLA( θf , θ) = τ 2 2 MSE( θf , θ) = 1 2τ 2 E Xf -X 2 F .
In the case of Gamma distributed measurements with known shape parameter L > 0, it follows that

MKLS( θf , θ) = L n i=1 m j=1 E   Xf ij X ij -log   Xf ij X ij   -1   , MKLA( θf , θ) = L n i=1 m j=1 E   X ij Xf ij -log   X ij Xf ij   -1   .
Below, we use some of the results in [START_REF] Deledalle | Estimation of Kullback-Leibler losses for noisy recovery problems within the exponential family[END_REF] whose main contributions are the derivation of new unbiased estimators of the MKLS and MKLA risks. For continuous exponential family, the risk estimate derived in [START_REF] Deledalle | Estimation of Kullback-Leibler losses for noisy recovery problems within the exponential family[END_REF] is unbiased for the MKLS risk, while it is only asymptotically unbiased for the MKLA risk with respect to the signal-to-noise ratio. For data sampled from a continuous exponential family, this makes simpler the use of the MKLS risk to derive data-driven shinkage in low rank matrix denoising, and we have therefore chosen to concentrate our study on this risk in this setting. The following proposition establishes a SURE formula to estimate the MKLS risk in the continuous case.

Proposition 2.2. Suppose that the data are sampled from a continuous exponential family.

Assume that the function h, in the definition (2.2) of the exponential family, is continuously differentiable on Y = R. Suppose that the function A, in the definition (2.2) of the exponential family, is twice continuously differentiable on Θ. If the following condition holds

E A ( θf ij (Y )) < +∞, for all 1 ≤ i ≤ n, 1 ≤ j ≤ m, (2.9)
then, the quantity

SUKLS( θf ) = n i=1 m j=1 θf ij (Y ) + h (Y ij ) h(Y ij ) A ( θf ij (Y )) -A( θf ij (Y )) + div f (Y ),
(2.10)

where div f (Y ) = m i=1 n j=1 ∂f ij (Y ) ∂Y ij , is an unbiased estimator of MKLS( θf , θ) - n i=1 m j=1 A(θ ij ).
A key difference in the formula of unbiased estimates for the MSE and the KL risks is the computation of the divergence term in (2.5) and (2.10), when

Xf = min(n,m) k=1 f k (σ k )ũ k ṽt
k is a smooth spectral estimator in the sense where each function f k : R + → R + is assumed to be (almost everywhere) differentiable for 1 ≤ k ≤ min(n, m). In this setting, the divergence term in the expression of GSURE( θf ) depends upon the matrix θf (Y ) = η( Xf ). Therefore, when η is a nonlinear mapping, it is generally not possible to obtain a simpler expression for div θf (Y ). To the contrary, for SUKLS( θf ), the divergence term is div f (Y ) which has the following closed-form expression for any smooth spectral estimators

div f (Y ) = |m -n| min(n,m) k=1 f k (σ k ) σk + min(n,m) k=1 f k (σ k ) + 2 min(n,m) k=1 f k (σ k ) min(n,m) =1; =k σk σ2 k -σ2 , (2.11)
thanks to the results from Theorem IV.3 in [START_REF] Candès | Unbiased risk estimates for singular value thresholding and spectral estimators[END_REF].

Note that SUKLS( Xr w ) = τ 2 2 SURE( Xr w ) for Gaussian measurements, hence, the GSURE and SUKLS strategies match in this case. In the case of Gamma measurements, assuming that L > 2 implies that the conditions on the function h in Proposition 2.2 is satisfied, and by assuming that condition (2.9) holds as well, it follows that

SUKLS( θf ) = n i=1 m j=1 (L -1) f ij (Y ) Y ij -L log (f ij (Y )) -Lmn + div f (Y ),
where the expression of div f (Y ) is given by (2.11). Note that it is implicitly understood in the definition of div f (Y ) that each mapping f ij : R n×m → R is differentiable. The differentiability of the spectral function f (and thus of its components f ij ) is a consequence of the assumption that the functions f 1 , . . . , f min(n,m) (acting on the singular values) are supposed to be differentiable. For further details, on the differentiability of f and the f ij 's, we refer to Section IV in [START_REF] Candès | Unbiased risk estimates for singular value thresholding and spectral estimators[END_REF]. From the arguments in [START_REF] Candès | Unbiased risk estimates for singular value thresholding and spectral estimators[END_REF], it follows that formula (2.11) for the divergence of f is also valid under the assumption that each function f k is differentiable on R + except on a set of Lebesgue measure zero.

The case of Poisson data

For Poisson data, the key result to obtain unbiased estimate of a given risk is the following lemma which dates back to the work in [START_REF] Hudson | A natural identity for exponential families with applications in multiparameter estimation[END_REF].

Lemma 2.1. Let f : Z n×m → R n×m be a measurable mapping. Let 1 ≤ i ≤ n and 1 ≤ j ≤ m, and denote by f ij : Z n×m → R a measurable function. Let Y ∈ Z n×m be a matrix whose entries are independently sampled from a Poisson distribution on Z. Then,

E   n i=1 m j=1 X ij f ij (Y )   = E   n i=1 m j=1 Y ij f ij (Y -e i e t j )   ,
where, for each 1 ≤ i ≤ n and 1 ≤ j ≤ m, f ij (Y ) denotes the (i, j)-th entry of the matrix f (Y ), and e i (resp. e j ) denotes the vector of Z n (resp. Z m ) with the i-th entry (resp. j-th entry) equals to one and all others equal to zero.

Hudson's lemma provides a way to estimate (in an unbiased way) the expectation of the Frobenius inner product between the matrix X and the matrix f (Y ). To see the usefulness of this result, one may consider the following mean-squared error

MSE( Xf , X) = E Xf -X 2 F = E   Xf 2 F -2 n i=1 m j=1 X ij Xf ij (Y ) + X 2 F   .
Therefore, by Lemma 2.1, one immediately obtains that

PURE( θf ) = Xf 2 F -2 n i=1 m j=1 Y ij f ij (Y -e i e t j ), (2.12)
is an unbiased estimate for the quantity MSE( Xf , X) -X 2 F . For Poisson data, one may also define the following KL risks

MKLS( θf , θ) = n i=1 m j=1 E   X ij - Xf ij - Xf ij log   X ij Xf ij     , (2.13) MKLA( θf , θ) = n i=1 m j=1 E   Xf ij -X ij -X ij log   Xf ij X ij     , ( 2
PUKLA( θf ) = n i=1 m j=1 Xf ij -Y ij log f ij (Y -e i e t j ) , (2.15) is an unbiased estimator of MKLA( θf , θ) + n i=1 m j=1 X ij -X ij log (X ij ).

Data-driven shrinkage in low-rank matrix denoising

For a matrix X with entries X ij ∈ X = R, we consider shrinkage estimators of the form

Xs w = f (Y ) = k∈s w k σk ũk ṽt k , (2.16) with s ⊆ I = {1, 2, . . . , min(n, m)} and w k ∈ [0, 1], for all k ∈ s.
When the underlying matrix X is constrained to have positive entries, e.g. X =]0, +∞[ in the Gamma and Poisson cases, we consider instead estimators of the form

Xs w = f (Y ) = max k∈s w k σk ũk ṽt k , ε , (2.17)
where ε > 0 is an a priori lower bound on the smallest value of X ij , where for any matrix X,

max[X, ε] ij = max[X ij , ε], for all 1 ≤ i ≤ m and 1 ≤ j ≤ n.
The construction of the subset s is postponed to Section 4, and we focus here in selecting the weights in a data-driven way for a fixed given s. In the following, we denote by s c the complementary set of s in I, i.e., s c = I \ s, and we let θs w = η Xs w . When X =]0, +∞[, we have found that considering estimators of the form (2.17) is more appropriate than trying to find shrinking weights (w k ) k∈s such that all the entries of the matrix k∈s w k σk ũk ṽt k are positive, for a given subset s.

Gaussian noise with known homoscedastic variance τ 2

By applying the GSURE formula (2.5) for Gaussian distributed measurements and thanks to the expression (2.11) for the divergence of smooth spectral estimators, we obtain for Xs w , as defined in (2.16), the SURE expression given by

SURE( Xs w ) = -mnτ 2 + k∈s (w k -1) 2 σ2 k + k∈s c σ2 k + 2τ 2 s k=1   1 + |m -n| + 2 min(n,m) =1; =k σ2 k σ2 k -σ2   w k
which unbiasedly estimate MSE( Xs w , X). Hence, for each k ∈ s, by differentiating the above expression with respect to w k , it follows that a data-driven weight for the k-th empirical singular value is given by

w k (Y ) =   1 - τ 2 σ2 k   1 + |m -n| + 2 min(n,m) =1; =k σ2 k σ2 k -σ2     + , (2.18)
which fullfils the requirement that w k (Y ) ∈ [0, 1]. Note that as SUKLS( Xs w ) = τ 2 2 SURE( Xs w ) for Gaussian measurements, the exact same data-driven weight would be obtained by minimizing an estimate of the MKLS( Xs w , X).

The case of estimators with rank one. Consider the case of estimators with rank 1, i.e., let s = {1}. It follows that

X1 w = X{1} w = w 1 X1 where w 1 ∈ [0, 1] is given by w 1 (Y ) =   1 - τ 2 σ2 1   1 + |m -n| + 2 min(n,m) =1; =1 σ2 1 σ2 1 -σ2     + .

Gamma and Poisson distributed measurements

In Gamma and Poisson cases, it is not possible to follow the same strategy as in the Gaussian case to derive optimal weights for (2.17) in a closed-form using the established SURE-like formulas.

We shall investigate how data-driven shinkage can be approximated in Section 5 on numerical experiments using fast algorithms. Nevertheless, when the estimator is restricted to rank 1, optimizing KL risk estimators lead to closed-form expressions under the assumption that all the entries of the data matrix Y are strictly positive.

The case of estimators with rank one under Gamma noise. Consider again the case of estimators with rank 1, i.e., let s = {1}, and let X1 = σ1 ũ1 ṽt 1 denote the PCA approximation of rank 1 of X. If all the entries of the matrix Y are strictly positive, by the Perron-Frobenius theorem, all the entries of the first singular vectors ũ1 and ṽ1 are strictly positive. Therefore, all the entries of X1 belong to the set X =]0, +∞[, and we can consider

X1 w = X{1} w = w 1 σ1 ũ1 ṽt 1
as defined in (2.16) instead of (2.17). Assuming L > 2 for the SUKLS formula to hold, it follows by simple calculations that

SUKLS( θ1 w ) = n i=1 m j=1 (L -1)w 1 X1 ij Y ij -mnL log (w 1 ) -L log X1 ij Y ij -Lmn + (1 + |m -n|)w 1 + 2w 1 min(n,m) =2 σ2 1 σ2 1 -σ2 .
Hence, by differentiating the above expression with respect to w 1 and as it is monotonic on both sides of its unique minimum, the optimal value of w 1 ∈ [0, 1] minimizing SUKLS( θ1 w ) is given by

w 1 (Y ) = min   1,   L -1 Lmn n i=1 m j=1 X1 ij Y ij + 1 Lmn   1 + |m -n| + 2 min(n,m) =2 σ2 1 σ2 1 -σ2     -1   ,
which yields the shrinking rule (1.11) stated in the introduction of this paper. Note that it is not possible to obtain, in a closed-form, the optimal value of the weight w 1 that minimizes the criterion GSURE( θ1 w ).

The case of estimators with rank one under Poisson noise. Using again that all the assumption that the entries of Y are positive, we can consider (by the Perron-Frobenius theorem) X1 w = X{1} w = w 1 σ1 ũ1 ṽt 1 as defined in (2.16) instead of (2.17). Then, the PURE formula (2.12) and Proposition 2.3 apply to the estimator θ1 w = log X1 w which yield to

PURE( θ1 w ) = w 2 1 σ2 1 -2 n i=1 m j=1 Y ij w 1 σ(ij) 1 ũ(ij) 1,i ṽ(ij) 1,j ,
and

PUKLA( θ1 w ) = n i=1 m j=1 w 1 X1 ij -Y ij log (w 1 ) + log σ(ij) 1 ũ(ij) 1,i ṽ(ij) 1,j
, where

X1 ij = σ1 ũ1,i ṽ1,j , σ(ij) 1
is the largest singular value of the matrix Ye i e t j , and ũ(ij) 1 (resp. ṽ(ij) 1 ) denotes its left (resp. right) singular vectors. Therefore, by differentiating the above expression with respect to w 1 and as it is monotonic on both sides of its unique minimum, an optimal value for w 1 ∈ [0, 1] which minimizes PURE( θ1 w ) is given by

w 1 (Y ) = min   1, 1 σ2 1 n i=1 m j=1 Y ij σ(ij) 1 ũ(ij) 1,i ṽ(ij) 1,j   .
However, this optimal shrinking rule cannot be used in practice since evaluating the values of

σ(ij) 1 , ũ(ij) 1 , ṽ(ij) 1
for all 1 ≤ i ≤ n and 1 ≤ j ≤ m is not feasible from a computational point of view for large values of n and m. Nevertheless, a fast algorithm to find a numerical approximation of the optimal value w 1 (Y ) is proposed in Section 5.

To the contrary, using again that all the X1 ij are positive by the Perron-Frobenius theorem, the value of

w 1 ∈ [0, 1] minimizing PUKLA( θ1 w ) is w 1 (Y ) = min   1, n i=1 m j=1 Y ij n i=1 m j=1 X1 ij   ,
which is straightforward to compute. This corresponds to the shrinkage rule (1.12) given in the introduction.

3 Gaussian spiked population model

In this section, we restrict our analysis to the Gaussian spiked population model and the asymptotic setting introduced in Definition 1.1.

Asymptotic location of empirical singular values

We

c ± = 1 ± √ c
where c + is the so-called bulk (right) edge. When X = 0 has a low rank structure, the asymptotic behavior of the singular values of Y = X + W is also well understood [BN12, DS07, SN13], and generalizations to noise matrix W whose distribution is orthogonally invariant have also been recently considered in [START_REF] Benaych-Georges | The singular values and vectors of low rank perturbations of large rectangular random matrices[END_REF]. Below, we recall some of these results that will be needed in this paper. To this end, let us introduce the real-valued function ρ defined by

ρ (σ) = (1 + σ 2 )(c + σ 2 ) σ 2
for any σ > 0.

Then, the following result holds (see e.g. Theorem 2.8 in [START_REF] Benaych-Georges | The singular values and vectors of low rank perturbations of large rectangular random matrices[END_REF] and Proposition 9 in [START_REF] Shabalin | Reconstruction of a low-rank matrix in the presence of Gaussian noise[END_REF]).

Proposition 3.1. Assume that Y = X + W is a random matrix sampled from the Gaussian spiked population model with τ = 1 √ m and X = r * k=1 σ k u k v t k . Then, for any fixed k ≥ 1, one has that, almost surely,

lim n→+∞ σk = ρ (σ k ) if k ≤ r * and σ k > c 1/4 , c + otherwise.
Moreover,

lim n→+∞ σmin(n,m) = c -.
In what follows, we shall also use the relation

1 σ 2 = ρ 2 (σ) -(c + 1) -(ρ 2 (σ) -(c + 1)) 2 -4c 2c that holds for any σ > c 1/4 , (3.1)
which is a consequence of e.g. the results in Section 3.1 in [START_REF] Benaych-Georges | The singular values and vectors of low rank perturbations of large rectangular random matrices[END_REF].

Existing asymptotically optimal shrinkage rules

Below, we briefly summarize some results in [START_REF] Gavish | Optimal shrinkage of singular values[END_REF] and [START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF] on the construction of asymptotically optimal spectral estimators. Let

Xf = f (Y ) = min(n,m) k=1 f k (σ k )ũ k ṽt k (3.2)
be a given smooth spectral estimator, and consider the standard squared error SE( Xf , X) = Xf -X 2 F as a measure of risk. The set of spectral functions minimizing SE( Xf , X) is given by f k (σ k ) = ũt k X ṽk , for 1 ≤ k ≤ min(n, m). However, it cannot be used in practice since X is obviously unknown. A first alternative suggested in [START_REF] Gavish | Optimal shrinkage of singular values[END_REF] and [START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF] is to rather study the asymptotic risk

SE ∞ ( Xf ) = lim n→∞ SE( Xf , X) (in the almost sure sense) (3.3)
in the Gaussian spiked population model. Then, it is proposed in [START_REF] Gavish | Optimal shrinkage of singular values[END_REF] and [START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF] to find an asymptotically optimal choice of f by minimizing SE ∞ ( Xf ) among a given class of smooth spectral functions. The results in [START_REF] Gavish | Optimal shrinkage of singular values[END_REF] show that, among spectral estimators of the form Xη = min(n,m) k=1

η(σ k )ũ k ṽt k , where η : R + → R + is a continuous shrinker such that η(σ) = 0 whenever σ ≤ c + , an asymptotically optimal shrinkage rule is given by the choice

η * (σ) = 1 σ (σ 2 -(c + 1)) 2 -4c if σ > c + , 0 otherwise.
(3.4)

In [START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF], it is proposed to consider spectral estimators of the form Xδ = r k=1 δ k ũk ṽt k where δ 1 , . . . , δ r are positive weights. By Theorem 2.1 in [START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF], it follows that, if σ k > c 1/4 for all 1 ≤ k ≤ r with r ≤ r * , then the weights which minimize SE ∞ ( Xδ ) over R r + are given by

δ * k = δ k (σ k ) = σ 4 k -c σ k (1 + σ 2 k )(c + σ 2 k ) , for all 1 ≤ k ≤ r. (3.5)
In what follows, the shrinkage rules (3.4) and (3.5) are shown to be equivalent, and they will serve as a reference of asymptotic optimality. It should be stressed that the estimators in [START_REF] Gavish | Optimal shrinkage of singular values[END_REF] and [START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF] are not equivalent. Indeed, the method in [START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF] requires an estimate of the rank, while the approach in [START_REF] Gavish | Optimal shrinkage of singular values[END_REF] applies the same shrinker to all empirical singular values. Nevertheless, the shrinkage function that is applied to significant singular values (either above the bulk edge in [START_REF] Gavish | Optimal shrinkage of singular values[END_REF] or up to a given rank in [START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF]) is the same.

Asymptotic behavior of data-driven estimators based on SURE

Following the principle of SURE, a second alternative to choose a smooth spectral estimator of the form (3.2) is to study the problem of selecting a set of functions

(f k ) 1≤k≤min(n,m) that minimize an unbiased estimate of MSE( Xf , X) = E Xf -X 2 F .
For any 1 ≤ i ≤ m and 1 ≤ j ≤ n, we recall that f ij (Y ) denotes the (i, j)-th entry of the matrix Xf = f (Y ). Under the condition that

E |Y ij f ij (Y )| + ∂f ij (Y ) ∂Y ij < +∞, for all 1 ≤ i ≤ n, 1 ≤ j ≤ m. (3.6)
it follows from the results in [START_REF] Candès | Unbiased risk estimates for singular value thresholding and spectral estimators[END_REF] (or equivalently from Proposition 2.1 for Gaussian noise with τ 2 = 1/m) that

SURE Xf = -n + f (Y ) -Y 2 F + 2 m div f (Y ), (3.7)
is an unbiased estimate of MSE( Xf , X), where the divergence div f (Y ) admits the closed-form expression (2.11). The SURE formula (3.7) has been used in [START_REF] Candès | Unbiased risk estimates for singular value thresholding and spectral estimators[END_REF] to find a data-driven value for λ = λ(Y ) in the the case of singular values shrinkage by soft-thresholding which corresponds to the choice

f k (σ k ) = (σ k -λ) + , for all 1 ≤ k ≤ min(n, m).
We study now the asymptotic behavior of the SURE formula (3.7). To this end, we shall use Proposition 3.1, but we will also need the following result (whose proof can be found in the Appendix) to study some of the terms in expression (2.11) of the divergence of f (Y ).

Proposition 3.2. Assume that Y = X + W is a random matrix sampled from the Gaussian spiked population model with τ = 1 √ m and X = r * k=1 σ k u k v t k .
Then, for any fixed 1 ≤ k ≤ r * such that σ k > c 1/4 , one has that, almost surely,

lim n→+∞ 1 n n =1; =k σk σ2 k -σ2 = 1 ρ (σ k ) 1 + 1 σ 2 k .
In what follows, we restrict our analysis to the following class of spectral estimators (the terminology in the definition below is borrowed from [START_REF] Gavish | Optimal shrinkage of singular values[END_REF]).

Definition 3.1. Let Xf = f (Y ) = min(n,m) k=1 f k (σ k )ũ k
ṽt k be a smooth spectral estimator. For a given 1 ≤ r ≤ min(n, m), the estimator f is said to be a spectral shrinker of order r that collapses the bulk to 0 if f k (σ) = 0 whenever σ ≤ c + and 1 ≤ k ≤ r, f k (σ) = 0 for all σ ≥ 0 and k > r.

The reason for restricting the study to spectral estimators such that f k (σ k ) = 0 whenever σk < c + is linked to the choice of the active set s * (1.8) of singular values in the Gaussian case, as detailed in Section 4. Now, for a spectral shrinker Xf of order r that collapses the bulk to 0, we study the asymptotic behavior of the terms in expression (3.7) that only depend on f , namely

SURE Xf = r k=1 (f k (σ k ) -σk ) 2 + 2 1 - n m r k=1 f k (σ k ) σk + 2 m r k=1 f k (σ k ) +4 n m r k=1 f k (σ k )   1 n n =1; =k σk σ2 k -σ2   (3.8)
The reason for studying SURE Xf is that finding an optimal shrinkage rule that minimizes

SURE

Xf is equivalent to minimizing expression (3.8) over spectral shrinkers of order r that collapses the bulk to 0, since SURE Xf -SURE Xf = -n + n k=r+1 σ2 k for such Xf .

Then, using Proposition 3.1, Proposition 3.2, and the assumption that the f k 's are continuously differentiable functions on R + , we immediately obtain the following result.

Lemma 3.1. Assume that Y = X + W is a random matrix sampled from the Gaussian spiked population model with τ = 1 √ m and X = r * k=1 σ k u k v t k . Let Xf be a spectral shrinker of order r ≤ r * that collapses the bulk to 0, such that each function f k , for 1 ≤ k ≤ r, is continuously differentiable on ]c + , +∞[. Moreover, assume that σ k > c 1/4 for all 1 ≤ k ≤ r. Then, one has that, almost surely,

lim n→+∞ SURE Xf = r k=1 (f k (ρ(σ k )) -ρ(σ k )) 2 + 2f k (ρ(σ k )) σ 2 k (1 + c) + 2c σ 2 k ρ(σ k ) (3.9)
Asymptotically optimal shrinkage of singular values. Thanks to Lemma 3.1, one may determine an asymptotic optimal spectral shrinker as the one minimizing lim n→+∞ SURE Xf .

For this purpose, let us define the class of estimators

Xr w = r k=1 w k σk 1 1 {σ k >c + } ũk ṽt k , (3.10)
where 1 ≤ r ≤ r * is a given integer, and the w k 's are positive weights. In practice, the estimator Xr w is computed by replacing the bulk edge c + by its approximation c n,m + = 1+ n m in eq. (3.10).

For moderate to large values of n and m, the quantities c + and c n,m + are very close, and this replacement does not change the numerical performances of Xr w . Then, provided that σ k > c 1/4 for all 1 ≤ k ≤ r, it follows from Lemma 3.1 that

lim n→+∞ SURE Xr w = r k=1 ρ 2 (σ k )(w k -1) 2 + 2w k σ 2 k (1 + c) + 2c σ 2 k .
Differentiating the above expression with respect to each weight w k leads to the following choice of asymptotically optimal weights

w * k = 1 - σ 2 k (1 + c) + 2c σ 2 k ρ 2 (σ k ) for all 1 ≤ k ≤ r. (3.11)
Therefore, if the singular values of the matrix X to be estimated are sufficiently large (namely σ k > c 1/4 for all 1 ≤ k ≤ r), by using Proposition 3.1 and eq. (3.11), one has that an asymptotically optimal spectral shrinker of order r ≤ r * is given by the choice of functions

f * k (ρ(σ k )) = 1 - σ 2 k (1+c)+2c σ 2 k ρ 2 (σ k ) ρ(σ k ) if ρ(σ k ) > c + , 0 otherwise, for all 1 ≤ k ≤ r. (3.12)
Using, the relation (3.1) one may also express the asymptotically optimal shrinking rule (3.12) either as a function of ρ(σ k ) only,

f * k (ρ(σ k )) = 1 ρ(σ k ) (ρ 2 (σ k ) -(c + 1)) 2 -4c if ρ(σ k ) > c + , 0 otherwise. (3.13)
or as function of σ k only (using that ρ(σ k ) > c + is equivalent to σ k > c 1/4 ),

f * k (ρ(σ k )) = σ 4 k -c σ k √ (1+σ 2 k )(c+σ 2 k ) if σ k > c 1/4 , 0 otherwise.
(3.14) Therefore, for spectral shrinker of order r, we remark that the shrinkage rule (3.13) coincides with the rule (3.4) which has been obtained in [START_REF] Gavish | Optimal shrinkage of singular values[END_REF]. Similarly, when the quantity f * k (ρ(σ k )) is expressed as a function of σ k only in (3.14), then we retrieve the shrinking rule (3.5) derived in [START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF]. Therefore, it appears that minimizing either the asymptotic behavior of the SURE, that is lim n→+∞ SURE Xf , or the limit of SE risk (3.3) leads to the same choice of an asymptotically optimal spectral estimator.

Data-driven shrinkage of empirical singular values. From the results in Section 2.3, the principle of SURE minimisation leads to the following data-driven choice of spectral shrinker of order r that collapses the bulk to 0

Xr w = r k=1 f k (σ k )ũ k ṽt k , (3.15) where f k (σ k ) = w k (Y )σ k 1 1 {σ k >c + } , for all 1 ≤ k ≤ r, with w k (Y )
given by (2.18). From Proposition 3.1 and Proposition 3.2 it follows that, if σ k > c 1/4 , then, almost surely,

lim n→+∞ f k (σ k ) = 1 - σ 2 k (1 + c) + 2c σ 2 k ρ 2 (σ k ) ρ(σ k ), for all 1 ≤ k ≤ r ≤ r * .
Therefore, the data-driven spectral estimator Xr w (3.15) asymptotically leads to the optimal shrinking rule of singular values given by (3.12) which has been obtained by minimizing the asymptotic behavior of the SURE.

Note that when τ = 1/ √ m, it suffices to replace the condition σk > c + by σk > τ (

√ m + √ n)
in the definition of Xr w , which yields the shrinking rule (1.13) stated in the introduction of this paper.

Estimating active sets of singular values in exponential families

In this section, we propose to formulate a new Akaike information criterion (AIC) to select an appropriate set of singular values over which a shrinkage procedure might be applied. To this end, we shall consider the estimator Xs = k∈s σk ũk ṽt k defined for a subset s ⊆ I = {1, 2, . . . , min(n, m)}, and we address the problem of selecting an optimal subset s from the data Y .

In the case of Gaussian measurements, the shrinkage estimators that we use in our numerical experiments are of the form

Xf = k∈s f k (σ k )ũ k ṽt k where s = {k ; σk > c n,m + } with c n,m + = 1 + n m ,
for some (possibly data-dependent) shrinkage functions f k . The set s is based on the knowledge of an approximation c n,m + of the bulk edge c + . Thanks to Proposition 3.1, the bulk edge c + is interpreted as the threshold which allows to distinguish the locations of significant singular values in the data from those due to the presence of additive noise. Interestingly, the following result shows that the active set s may be interpreted through the prism of model selection using the minimisation of a penalized log-likelihood criterion.

Proposition 4.1. Assume that Y = X + W where the entries of W are iid Gaussian variables with zero mean and standard deviation τ = 1/ √ m. Then, we have

s * = arg min s⊆I m Y - Xs 2 F + 2|s|p n,m with p n,m = 1 2 √ m + √ n 2 , (4.1)
where Xs = k∈s σk ũk ṽt k for s ∈ I = {1, 2, . . . , min(n, m)}, and |s| is the cardinal of s.

Proof. We remark that Y -

Xs = k / ∈s σk ũk ṽt k . It results that m Y - Xs 2 F + 2|s|p n,m = m k / ∈s σ2 k + 2|s|p n,m = n k=1 mσ 2 k if k / ∈ s 2p n,m otherwise . (4.2)
Using that 2p n,m /m = c n,m + , it follows that the set s = {k ; σk > c n,m + } is by definition such that k ∈ s if and only if 2p n,m < mσ 2 k . Therefore, by (4.2), the criterion s → m Y -Xs 2 F + 2|s|p n,m is minimum at s = s which concludes the proof.

In the model Y = X + W , where the entries of W are iid Gaussian variables with zero mean and variance τ 2 , it is well known that the degrees of freedom (DOF) of a given estimator X is defined as

DOF( X) = 1 τ 2 n i=1 m j=1 Cov( Xij , Y ij ) = 1 τ 2 n i=1 m j=1 E[ Xij W ij ].
The DOF is widely used in statistics to define various criteria for model selection among a collection of estimators, see e.g. [START_REF] Efron | The estimation of prediction error: Covariance penalties and crossvalidation[END_REF]. In low rank matrix denoising, the following proposition shows that it is possible to derive the asymptotic behavior of the DOF of spectral estimators.

Proposition 4.2. Assume that Y = X + W is a random matrix sampled from the Gaussian spiked population model with

τ = 1 √ m and X = r * k=1 σ k u k v t k . Let
Xf be a spectral shrinker of order r ≤ r * that collapses the bulk to 0, such that each function

f k , for 1 ≤ k ≤ r, is continuously differentiable on ]c + , +∞[. Moreover, assume that σ k > c 1/4 for all 1 ≤ k ≤ r.
Then, one has that, almost surely,

lim n→+∞ 1 m DOF( Xf ) = r k=1 f k (ρ(σ k )) ρ(σ k ) 1 + c + 2c σ 2 k .
Proof. Thanks to the derivation of the SURE in [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF] and formula (2.11) on the divergence of spectral estimators, one has that

DOF( Xf ) = E div Xf = E   |m -n| r k=1 f k (σ k ) σk + r k=1 f k (σ k ) + 2 r k=1 f k (σ k ) n =1; =k σk σ2 k -σ2   .
By Proposition 3.1, Proposition 3.2, and our assumptions on the f k 's, one has that, almost surely,

lim n→+∞ 1 m div Xf = r k=1 f k (ρ(σ k )) ρ(σ k ) 1 + c + 2c σ 2 k .
which completes the proof.

Hence, in the Gaussian spiked population model, by Proposition 4.2 and using that

σ 2 k > √ c
for all 1 ≤ k ≤ r, it follows that if s ⊆ {1, . . . , r} then

lim n→+∞ 1 m DOF( Xs ) = |s| 1 + c + 2c σ 2 k ≤ |s| 1 + √ c 2 = |s|c 2 + . (4.3)
Hence, the quantity 2|s| 1 2 ( √ m + √ n) 2 is asymptotically an upper bound of DOF( Xs ) (when normalized by 1/m) for any given set s ⊆ {1, . . . , r}.

Let us now consider the more general case where the entries of Y are sampled from an exponential family. To the best of our knowledge, extending the notion of the bulk edge to non-Gaussian data sampled from an exponential family has not been considered so far in the literature on random matrices and low rank perturbation model. Therefore, except in the Gaussian case, it is far from being trivial to find an appropriate threshold value c to define an active set in the form s = {k ; σk > c}.

Nevertheless, to select an appropriate active set of singular values, we introduce the following criterion that is inspired by the previous results on the DOF of the estimator Xs in the Gaussian case and the statistical literature on the well known AIC for model selection [START_REF] Akaike | A new look at the statistical model identification[END_REF].

Definition 4.1. The AIC associated to Xs = k∈s σk ũk ṽt k is

AIC( Xs ) = -2 log q(Y ; Xs ) + 2|s|p n,m with p n,m = 1 2 √ m + √ n 2 . (4.4)
where |s| is the cardinal of s, and q(Y ;

Xs ) = n i=1 m j=1 q(Y ij ; Xs ij )
is the likelihood of the data in the general form (2.1) at the estimated parameters

X ij = Xs ij .
In the above definition of AIC( Xs ), the quantity 2|s|p n,m is an approximation of the degree of freedom of Xs , i.e., of the numbers of its free parameters as it is justified by Proposition 4.2 in the case of Gaussian measurements. The AIC allows us to define an optimal subset of active variables as

s * = arg min s⊆I AIC( Xs ).
For Gaussian measurements, Proposition 4.1 gives the value of the optimal set s * in a closed-form. Following the arguments in Section 2.3, for Gamma or Poisson measurements and for a given subset s, we consider the estimator

Xs = max k∈s σk ũk ṽt k , ε , (4.5)
when > 0 is an a priori value to satisfy the positivity constraint on the entries of an estimator in this setting. However, contrary to the case of Gaussian noise, the search of an optimal subset s ⊂ arg min s⊆I AIC( Xs ) becomes a combinatorial problem in this context. In our numerical experiments, we thus choose to construct an approximation s of s with a greedy search strategy that reads as follows

s = I \ k ∈ I ; AIC( XI\{k} ) ≤ AIC( XI ) . (4.6)
For Gaussian measurements, s = s since the optimisation problem (4.6) becomes separable. In our numerical experiments, we have found that s selects a relevant set of active singular values which separates well the structural content of X while removing most of the noise component. Further details are given in Section 5 below.

For Gaussian noise, the computation of the active set s * of singular values may also be interpreted as a way to estimate the unknown rank r * of the signal matrix X. In this setting, one has that s = {k ; σk > c n,m + } which suggests the choice r = max{k ; σk > c n,m + }, (4.7)

as an estimator of r * .

There exists an abundant literature of the problem of estimating the rank of an empirical covariance matrix for the purpose of selecting the appropriate number of significant components to be kept in PCA or factor analysis. It is much beyond the scope of this paper to give an overview of this topic. We point to the review in [START_REF] Jolliffe | Principal component analysis[END_REF] for a summary of existing methods to determine the number of components in PCA that are grouped into three categories: subjective methods, distribution-based test tools, and computational procedures. For recent contributions in the matrix denoising model (1.1) with Gaussian noise, we refer to the works [START_REF] Choi | Selecting the number of principal components: Estimation of the true rank of a noisy matrix[END_REF][START_REF] Gavish | The optimal hard threshold for singular values is \(4/\sqrt {3}\)[END_REF] and references therein. For example for Gaussian data with know variance τ 2 = 1/m, Eq. ( 11) in [START_REF] Gavish | The optimal hard threshold for singular values is \(4/\sqrt {3}\)[END_REF] on optimal hard thresholding of singular values suggest to take

r = max{k ; σk > λ(c)}, with λ(c) = 2(c + 1) + 8c (c + 1) + √ c 2 + 14c + 1 , (4.8)
as a simple method to estimate the rank. It should be remarked that the problem of estimating the true rank r * of X in model (1.1) is somewhat ill-posed as, in the Gaussian spiked population model, Proposition 3.1 implies that one may only expect to estimate the so-called effective rank r eff = max{k ; σ k > c 1/4 } (see e.g. Section II.D in [START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF]).

In our numerical experiments, we shall compare different choices for the active set of singular values of the form ŝ = {1, . . . , r} where r is either given by (4.7), (4.8), or by the "oracle choices" r = r * and r = r eff .

Other methods based on hypothesis testing [START_REF] Choi | Selecting the number of principal components: Estimation of the true rank of a noisy matrix[END_REF] could be used for rank estimation in the Gaussian model (1.1), but it is beyond the purpose of this paper to give a detailed comparison.

For Poisson or Gamma noise, it is more difficult to interpret the computation of s * as a way to estimate the rank of X since, in our numerical experiments, we have found that the cardinality of s * is generally not equal to max{k ; k ∈ s * }. Moreover, to the best of our knowledge, there is not so much work on the estimation of the true rank of a noisy matrix beyond the Gaussian case. Therefore, we have not included a numerical comparison with other methods for the choice of the active set of singular values in these two cases.

The case of a signal matrix of rank one

We consider the simple setting where the rank r * of the matrix X is known and equal to one meaning that

X = σ 1 u 1 v t 1 ,
where u 1 ∈ R n and v 1 ∈ R m are vectors with unit norm that are fixed in this numerical experiment, and σ 1 is a positive real that we will let varying. We also choose to fix n = m = 100, and so to take c = n m = 1 and c + = 2. For the purpose of sampling data from Gamma and Poisson distribution, we took singular vectors u 1 and v 1 with positive entries. The i-th entry (resp. j-th entry) of u 1 (resp. v 1 ) is chosen to be proportional to 1

-(i/n -1/2) 2 (resp. 1 -(j/m -1/2) 2 ). Let Y = min(n,m) k=1
σk ũk ṽt k be an n × m matrix whose entries are sampled from model (2.1) and then satisfying E[Y ] = X.

Gaussian measurements

We first consider the case of Gaussian measurements, where

Y = X + W with E[W ij ] = 0, Var(W ij ) = τ 2 with τ = 1 √ m .
In this context, we compare the following spectral shrinkage estimators:

• Rank-1 PCA shrinkage X1 = σ1 ũ1 ṽt 1 ,
• Rank-1 SURE-driven soft-thresholding

X1 soft = σ1 ũ1 ṽt 1 with σ1 = (σ 1 -λ(Y )) + ,
• Rank-1 asymptotically optimal shrinkage proposed in [START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF] and [GD14a]

X1 * = σ1 ũ1 ṽt 1 with σ1 = σ2 1 -4 1 1 {σ 1 >2} ,
• Rank-1 SURE-driven weighted estimator that we have derived in Section 2.3

X1 w = σ1 ũ1 ṽt 1 with σ1 = 1 - 1 σ2 1 1 m + 2 m n =2 σ2 1 σ2 1 -σ2 + σ1 1 1 {σ 1 >2} ,
where the above formula follows from the results in Section 3.3 using that c = 1 and c + = 2 in these numerical experiments, and where, for the soft-thresholding, the value λ(Y ) > 0 is obtained by a numerical solver in order to minimize the SURE. As a benchmark, we will also consider the oracle estimator X 1 * that performs shrinkage by using the knowledge of the true singular-value σ 1 defined as

X 1 * = σ1 ũ1 ṽt 1 with σ1 = ρ(σ 1 ) 2 -4 1 1 {ρ(σ 1 )>2}
which corresponds to the asymptotically optimal shrinking rule (3.13) as a function of ρ(σ k ) in the setting c = 1 and c + = 2. Note that form the formula above ŵ1 = σ1 / σ1 is necessary in the range [0, 1] for all considered estimators.

In Figure 1, we compare the estimated singular-values σ1 and the estimated weights ŵ1 = σ1 /σ 1 as functions of σ 1 for the four aforementionned estimators. Because all estimators are subject to noise variance, we display, for all estimators, the median values and the 80% confidence intervals obtained from M = 100 noise realizations. It can be seen that the median curves for the eigenvelues and the weights of X1 w and X1 * coincide (up to variations that are slightly larger for the former) which is in agreement with the asymptotic analysis of shrinkage rules that has been carried out in Section 3.3. Spectral estimator obtained by SURE-driven soft-thresholding also leads to an optimal shrinkage rule.

In Figure 2, for each of the four spectral estimators above, we display for M = 100 noise realizations, as functions of σ 1 , the following normalized MSE

NMSE( X) = X -X 2 F X 2 F .
The normalized MSE of the estimators X1 soft , X1 * and X1 w are the same for values of σ 1 larger than c 1/4 = 1, and they only differ for values of σ 1 close or below the threshold c 1/4 = 1 (corresponding to values of ρ(σ 1 ) below the bulk edge c + = 2). More remarkably, above c 1/4 = 1, they offer similar NMSE values to the oracle shrinkage estimator X 1 * , not only in terms of median but also in terms of variability, as assessed by the confidence intervals. The performances of the estimator X1

(standard PCA) are clearly poorer. These numerical experiments also illustrate that, for finitedimensional low rank matrix denoising with r * = 1, data-driven spectral estimators obtained by minimizing a SURE criterion achieve performances that are similar to asymptotically optimal shrinkage rules.

Gamma and Poisson distributed measurements

Let us now consider the case where the entries of Y ij ≥ 0 of the data matrix Y are independently sampled from a Gamma or Poisson distribution with mean X ij > 0. To satisfy the constraint that the estimators must be matrices with positive entries, we consider estimators of the form (2.17). In this context, we compare the following spectral shrinkage estimators, set for ε = 10 -6 , as:

• Rank-1 PCA shrinkage X1 = max σ1 ũ1 ṽt 1 , ε ,
• Rank-1 GSURE/SUKLS/PURE/PUKLA-driven soft-thresholding • Rank-1 GSURE/SUKLS/PURE/PUKLA-driven weighted estimator

X1 soft = max σ1 ũ1 ṽt 1 , ε with σ1 = (σ 1 -λ(Y )) + ,
X1 w = max σ1 ũ1 ṽt 1 , ε with σ1 = w 1 (Y )σ 1 1 1 {1∈s} ,
where s is the approximated active subset as defined in Section 4. For the soft-thresholding, the value λ(Y ) > 0 is obtained by a numerical solver in order to minimize either the GSURE or the SUKLS criterion (in the Gamma case) and either the PURE or the PUKLA criterion (in the Poisson case). The weight w 1 (Y ) ∈ [0, 1] is obtained by a numerical solver in order to minimize the GSURE and the PURE, as described in Section B. According to Section 2.3, the weight w 1 (Y ) ∈ [0, 1], minimizing the SUKLS criterion, has the following closed-form formula

w 1 (Y ) = min   1,   L -1 Lmn n i=1 m j=1 X1 ij Y ij + 1 Lmn 1 + 2 n =2 σ2 1 σ2 1 -σ2   -1   ,
and for the PUKLA criterion, we have

w 1 (Y ) = min   1, n i=1 m j=1 Y ij n i=1 m j=1 X1 ij   .
To evaluate the performances of these estimators, we perform again a study involving M = 100 noise realizations.

In the Gamma case with shape parameter L = 3, results are reported in Figure 3 where σ 1 ranges from 0.1 to 5. In the Poisson case, results are reported in Figure 4. To generate data from a Poisson distribution with mean value X = σ 1 u 1 v t 1 , we took σ 1 ranging from 25 to 400. In this context, the entries X i,j are in average ranging from 0.25 to 4. When σ 1 = 25, about 78% of the entries of Y are 0 and 20% are equals to 1 which correspond to an extreme level of noise, while when σ 1 = 400, the entries of Y concentrate around 4 with a standard deviation of 2 which correspond to a simpler noisy setting.

In these experiments, it can be seen that all the data-dependent spectral estimators achieve comparable results with really small errors in terms of MSE and MKL risks. Their performances are similar to X1 = σ1 ũ1 ṽt 1 meaning that optimizing either SURE-like criteria leads to a spectral estimator closed to correspond to matrix denoising by ordinary PCA. However, unlike the Gamma case, it might be observed in the Poisson case that when reaching a stronger noise level, i.e, for small value of σ 1 , the NMSE of all estimator increases as the denoising problem becomes more challenging. Nevertheless, only the weight of X1 w driven by PUKLA does not present a drop wich allows reaching a slightly smaller MKLA. In the Gamma case, the noise level being proportional to the signal level, the NMSE/MKLS remain constant for all σ 1 . Finally, it should be remarked that a comparison with the asymptotically optimal shrinkage rule X1 * proposed in [START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF] and [START_REF] Gavish | Optimal shrinkage of singular values[END_REF] for Gaussian noise with variance τ 2 = 1 m is not realistic in the case of Gamma or Poisson measurements. Indeed, as remark in [START_REF] Gavish | Optimal shrinkage of singular values[END_REF], to use the estimator X1 * in a Gaussian model with homoscedastic variance τ 2 = 1 m , one may take the estimator

X1 * = √ mτ f * 1 (σ 1 /( √ mτ ))ũ 1 ṽt 1 .
However, this clearly means that the approach in [START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF] and [START_REF] Gavish | Optimal shrinkage of singular values[END_REF] requires the knowledge of the variance of the entries Y ij of the data matrix Y which is not possible for Gamma or Poisson measurements as either Var(

Y ij ) = X 2 ij /L or Var(Y ij ) = X ij in these settings.
While all estimators behave similarly in the rank 1 setting, we will see in the next section that they can significantly differ when the rank is let to be larger than 2.

The case of a signal matrix of rank larger than two

We now consider the more complex an realistic setting where the rank r * of the matrix X is unknown and potentially larger than two, i.e.,

X = r * k=1 σ k u k v t k ,
where u k ∈ R n and v k ∈ R m are vectors with unit norm that are fixed in this numerical experiment, and σ k are positive real values also fixed in this experiment. We also choose to fix n = 100 and m = 200, while the true rank is r * = 9 as shown by the red curve in Figure 5 . (i) Their corresponding singular values. (j) NMSE of the various approximations as a function of the rank r. (k) Same but without knowledge the bulk edge, namely c + = 0. (l,m,n) Same when the active set of singular values is of the form ŝ = {1, . . . , r} where r is either given by r = r * (oracle/true rank), r = r eff (effective rank) or by (4.8). In all the figures, the solid curves correspond to oracle estimators and the dashed curves correspond to data-driven estimators, obtained over M = 1, 000 noise realizatrions. The grey areas represent a 80% confidence interval. • GSURE/SUKLS/PURE/SUKLA driven weighted estimator

• (a) • (b) • (c) • (d) • (e) • (f) • (g) • (h) • (i) • (j)
Xr w = r k=1 max σk ũk ṽt k , ε with σk = w k (Y )σ k 1 1 {k∈s} ,
where r ∈ [1, min(n, m)], and s is the approximated active subset as defined in Section 4. For the soft-thresholding, the value λ(Y ) > 0 is obtained by a numerical solver in order to minimize either the GSURE or the SUKLS criterion (in the Gamma case) and either the PURE or the PUKLA criterion (in the Poisson case). As shown in Section 2.3, in the case of Gamma (resp. Poisson) measurements, the value of w k (Y ) for k ∈ s which minimizes the GSURE (resp. PURE) or the SUKLS (resp. PUKLA), cannot be obtained in closed form. As an alternative, we adopt a greedy one-dimensional optimization strategy starting from the matrix σ1 ũ1 ṽt 1 and next updating the weights w sequentially by starting = 1 to = min(n, m), with the constraint that, for all / ∈ s, the weight w is set to zero. To this end, we resort to one-dimensional optimization techniques in the interval [0, 1] using Matlab's command fminbnd. This strategy is used for GSURE, SUKLS, PURE and PUKLA by evaluating them as described in Section B. As in the Gaussian setting, we compare this spectral estimators with their oracle counterparts given by

X soft = min(m,n) k=1 max σk ũk ṽt k , ε with σk = (σ k -λ oracle (Y )) + , and
X r w = r k=1 max σk ũk ṽt k , ε with σk = w oracle k (Y )σ k 1 1 {k∈s} .
where w oracle k (X) and λ oracle (Y ) minimizes one of the objective SE η , KLS, SE or KLA (nonexpected risks) over the set of matrices sharing with Y the same r first left and right singular vectors, and soft-thresholding approximations respectively. Note again that X r w and X soft are ideal approximations of X that cannot be used in practice but serve as benchmarks to evaluate the performances of the data-driven estimators Xr w and Xsoft . The results for the Gamma noise are reported on Figure 6. As in the Gaussian setting, it can be observed that Xr w and X r w achieve comparable performances, as well as Xsoft and X soft showing that the GSURE (resp. SUKLS) accurately estimates the MSE η (resp. KLS). Visual inspection of the restored matrices tends to show that the estimators driven by MSE η or GSURE produce less relevant results compared to KLS or SUKLS, as confirmed by the curves of NMSE and MKLS. Performance in terms of NMSE also illustrates that minimizers of SE η do not coincides with those of SE. As in the Gaussian setting, Xr w and X r outperform Xsoft , X soft and standard PCA Xr provided that r is large enough. Moreover, the performance of Xr w obtained with KL objectives plateaus to its optimum when the rank r becomes large. Again, this allows us to choose r = min(n, m) when we do not have a priori on the true rank r .

The results for the Poisson noise are reported on Figure 7. The conclusions are similar to the Gaussian and Gamma cases. Obviously, the NMSE is smaller for approximations that minimizes SE (or PURE) than for those minimizing KLA (or PUKLA). However, visual inspection of the obtained matrices tends to demonstrate that minimizing such objectives might be less relevant than minimizing KL objectives. In this setting, the performance of Xr w is on a par with the one

A.2 A technical result to prove SURE-like formulas

We recall the key lemma needed to prove the SURE-like formulas in an exponential family in the continuous case. Similar results have already been formulated in different papers in the literature, see e.g. the review proposed in [START_REF] Deledalle | Estimation of Kullback-Leibler losses for noisy recovery problems within the exponential family[END_REF].

Lemma A.1. Let Y ∈ R n×m be a random matrix whose entries Y ij are independently sampled from the continuous exponential family (2.2) in canonical form (that is the distribution of Y ij is absolutely continuous with respect to the Lebesgue measure dy on R). Suppose that the function h is continuously differentiable on Y = R. Let 1 ≤ i ≤ n and 1 ≤ j ≤ m, and denote by F ij : R n×m → R a continuously differentiable function such that

E [|F ij (Y )|] < +∞. (A.4)
Then, the following relation holds

E [θ ij F ij (Y )] = -E h (Y ij ) h(Y ij ) F ij (Y ) + ∂F ij (Y ) ∂Y ij .
Proof. Using the expression (2.2) of the pdf of the random varibles Y ij , one has that , an integration by part and eq. (A.5) imply that

E [θ ij F ij (Y )] = R n×m F ij (Y )h(y ij )θ ij exp (θ ij y ij -A(θ ij )) dy ij
E [θ ij F ij (Y )] = - R n×m ∂F ij (Y )h(y ij ) ∂y ij exp (θ ij y ij -A(θ ij )) dy ij n 1≤k≤n 1≤ ≤m (k, ) =(i,j) p(y k ; θ k ) dy k . Now, since ∂F ij (Y )h(y ij ) ∂y ij = h (y ij )F ij (Y ) + ∂F ij (Y )
∂y ij h(y ij ), we finally obtain that

E [θ ij F ij (Y )] = -E h (Y ij ) h(Y ij ) F ij (Y ) + ∂F ij (Y ) ∂Y ij ,
which completes the proof.

A.3 Proof of Proposition 2.1

We remark that

E θ ij θf ij (Y ) = E h (Y ij ) h(Y ij ) θf ij (Y ) + E   ∂ θf ij (Y ) ∂Y ij   . (A.7)
Then, by definition (2.2) of the exponential family, we remark that

E h (Y ij ) h(Y ij ) = R h (y ij ) exp (θ ij y ij -A(θ ij )) dy ij .
Hence, using an integration by parts twice, we arrive at

E h (Y ij ) h(Y ij ) = θ 2 ij R h(y ij ) exp (θ ij y ij -A(θ ij )) dy ij = θ 2 ij . (A.8)
To complete the proof, it suffices to insert equalities (A.7) and (A.8) into (A.6).

A.4 Proof of Proposition 2.2

Thanks to eq. (2.8), one has that

MKLS( θf , θ) = n i=1 m j=1 E θf ij (Y )A ( θf ij (Y )) -θ ij A ( θf ij (Y )) -A( θf ij (Y )) + A(θ ij ).
(A.9) Using Lemma A.1 with F ij (Y ) = A ( θf ij (Y )) and condition (2.9), it follows that

E θ ij A ( θf ij (Y )) = -E h (Y ij ) h(Y ij ) A ( θf ij (Y )) -E   ∂ θf ij (Y ) ∂Y ij A ( θf ij (

B Implementation details

We discuss below an algorithmic approach to find data-driven spectral estimators. First, we discuss on how to compute data-driven spectral estimators from the expression of risk estimators. For SUKLS in continuous exponential families, and for SURE in the Gaussian case only, eq. (3.7) and (2.10) provide respectively a closed-form solution that can be evaluated in linear time O(nm). On the contrary, the computations of GSURE (beyond the Gaussian case), PURE and PUKLA, given respectively in eq. (2.6), (2.12) and (2.15), cannot be evaluated in reasonable time. They rely respectively on the computation of the divergence div θf (Y ),

Y ij f ij (Y -e i e t j
) and Y ij log f ij (Ye i e t j ) . Without further assumptions, such quantities requires O(n 2 m 2 ) operations in general. A standard approach for the computation of the divergence, suggested in [START_REF] Girard | A fast monte-carlo cross-validation procedure for large least squares problems with noisy data[END_REF][START_REF] Ramani | Monte-Carlo SURE: a black-box optimization of regularization parameters for general denoising algorithms[END_REF], is to unbiasedly estimate it with Monte-Carlo simulations by sampling the following relation

Y ij f ij (Y -e i e t j ) ≈ Y ij f ij (Y ) -δ i,j ∂f (Y ) ∂Y δ i,j
, and

Y ij log f ij (Y -e i e t j ) ≈ Y ij log f ij (Y ) -δ i,j ∂f (Y ) ∂Y δ i,j
where the entries of δ should be chosen Bernoulli distributed with parameter p = 0.5. The advantage of these three approximations is that they can be computed where Ũ and Ṽ are the matrices whose columns are ũk and ṽk , and D, S and A are n × m matrices defined, for all 1 ≤ i ≤ n and 1 ≤ j ≤ m, as

D i,j = δi,j × f i (σ i ) if i = j 0 otherwise, S i,j = δi,j + δj,i 2 × 0 if i = j f i (σ i )-f j (σ j ) σi -σ j otherwise,
A i,j = δi,jδj,i 2 × 0 if i = j f i (σ i )+f j (σ j ) σi +σ j otherwise,

where σk and f k (σ k ) are extended to 0 for k > min(n, m) and δ = Ũ t δ Ṽ ∈ R n×m .

  summarize below the asymptotic behavior of the singular values of the data matrix Y = min(n,m) k=1 σk ũk ṽt k in the Gaussian spiked population model. In the case where X = 0, it is well known [AGZ10, BS10] that the empirical distribution of the singular values of Y = W (with τ = 1 √ m ) converges, as n → +∞, to the quarter circle distribution if c = 1 and to its generalized version if c < 1. This distribution is supported on the compact interval [c -, c + ] with

Figure 1 :

 1 Figure 1: The case of Gaussian measurements with m = n = 100. Estimated first singular value σ1 as a function of the true underlying one σ 1 , for (a) our proposed estimator X1 w , (b) the soft-thresholding X1 soft and (c) the asymptotical one X1 * . All of them are compared to the first singular value σ1 of Y 1 and the one of the oracle asymptotical estimator X 1 * . (c,d,e) Same but for the corresponding weight ŵ1 = σ1 /σ 1 . Curves have been computed on M = 100 noise realizations, only the median and an 80% confidence interval are represented respectively by a stroke and a shadded area of the same color.

Figure 2 :

 2 Figure 2: Same as Fig. 1 but for the normalized MSE of the corresponding estimators.

Figure 3 :

 3 Figure 3: The case of Gamma measurements with m = n = 100. (a) Estimated first eigenvalue σ1 as a function of the true underlying one σ 1 for our proposed estimator X1 w and the soft-thresholding X1 soft when both are guided by the GSURE. Both of them are compared to the first singular value σ1 of Y 1 . Same but for (b) the corresponding weights ŵ1 = σ1 /σ 1 , (c) the NMSE risk and (d) the MKLS risk. (e-h) Exact same esperiments but when our proposed estimator and the soft-thresholding are both guided by SUKLS. Curves have been computed on M = 100 noise realizations, only the median and an 80% confidence interval are represented respectively by a stroke and a shadded area of the same color.

Figure 4 :

 4 Figure 4: The case of Poisson measurements with m = n = (a) Estimated first eigenvalue σ1 as a function of the true underlying one σ 1 for our proposed estimator X1 w and the soft-thresholding X1 soft when both are guided by the PURE. Both of them are compared to the first singular value σ1 of Y 1 . Same but for (b) the corresponding weights ŵ1 = σ1 /σ 1 , (c) the NMSE risk and (d) the MKLA risk. (e-h) Exact same esperiments but when our proposed estimator and the soft-thresholding are both guided by PUKLA. Curves have been computed on M = 100 noise realizations, only the median and an 80% confidence interval are represented respectively by a stroke and a shadded area of the same color.

  (i). Again, let Y = min(n,m) k=1 σk ũk ṽt k be an n × m matrix whose entries are sampled from model (2.1) and then satisfying E[Y ] = X.

•

Figure 5 :

 5 Figure 5: (a) Zoom on a 100 × 200 matrix and (e) a single realization of corrupted version by Gaussian noise (τ = 80). (b,c) Oracle soft-thresholding X soft and data-driven soft-thresholding Xsoft . (d) PCA full rank Xrmax , i.e., r max = min(n, m). (f,g,h) Oracle full rank approximation X rmax w , and data-driven full rank estimation Xrmax w

Figure 6 :

 6 Figure 6: (a) A single realization of version by Gamma noise (L = 80) with zoom on a 100 × 200 matrix. (b,c,d,e) Oracle soft-thresholding X soft and datadriven soft-thresholding Xsoft respectively for SE η , GSURE, KLS and SUKLS. (f) PCA Xrmax with full rank approximation i.e. r max = min(n, m). (g,h,i,j) Oracle full rank approximation X rmax w , and data-driven full rank estimation Xrmax w respectively for SE η , GSURE, KLS and SUKLS. (k) Their corresponding singular values averaged over M = 100 noise realizations. (l,m) NMSE averaged over M = 100 noise realizations as a function of the rank r with and without using the active set. (n,o) Same but with respect to MKLS.

•Figure 7 :

 7 Figure 7: (a) A single realization of corrupted version by Poisson noise with zoom on a 100 × 200 noise-free matrix (b,c,d,e) Oracle soft-thresholding X soft and datadriven soft-thresholding Xsoft respectively for SE, PURE, KLA and PUKLA. (f) PCA Xrmax with full rank approximation i.e. r max = min(n, m). (g,h,i,j) Oracle full rank approximation X rmax w , and data-driven full rank estimation Xrmax w respectively for SE, PURE, KLA and PUKLA. (k) Their corresponding singular values averaged over 200 noise realizations. (l,m) NMSE averaged over 200 noise realizations as a function of the rank r with and without using the active set. (n,o) Same but with respect to MKLA. (Matrix entries are displayed in log-scale for better visual assessment.)

  ) =(i,j) p(y k ; θ k ) dy k .where Y = (y k ) 1≤k≤n,1≤ ≤m . Thanks to condition (A.4), it follows thatR n×m F ij (Y )h(y ij ) exp (θ ij y ij -A(θ ij )) dy ij n 1≤k≤n 1≤ ≤m (k, ) =(i,j) p(y k ; θ k ) dy k < +∞. (A.5) Therefore, given that θ ij exp (θ ij y ij -A(θ ij )) = ∂ exp(θ ij y ij -A(θ ij )) ∂y ij

  ij (Y )| 2 -2θ ij θf ij (Y ) + θ 2 ij . (A.6) Using Lemma A.1 with F ij (Y ) = θf ij(Y ) and condition (2.4), it follows that

XY

  θ ij ). Now recall that f ij (Y ) = η -1 θf ij (Y ) and that A ( θf ij (Y )) = η -1 θf ij (Y )by Assumption 2.1. Therefore, Thanks to the expression (2.14) of the MKLA risk for data sampled from a Poisson distribution, it follows thatMKLA( θf , θ) + ij -X ij log (X ij ) = Xf ij -X ij log Xf ijIn the case of Poisson data, one has that exp (θ ij ) = X ij andh(Y ij -1) h(Y ij ) = Y ij .Therefore, by applying Hudson's Lemma 2.1 with F ij (Y ) = log Xf ij , it follows that ij log f ij (Ye i e t j )

  div θf (Y ) = E δ tr δ t ∂ θf (Y ) ∂Y δ at random directions δ ∈ R n×m satisfying E[δ] = 0, E[δ i δ i ] = 1 and E[δ i δ j ] = 0. Following [Del15], a similar first order approximation can be used for the other two quantities as

  in linear time O(nm) by making use of the results of [LS01, SS03, Ede05, CSLT13, DVP + 12] that provide an expression for the directional derivative given by ∂f (Y ) ∂Y δ = Ũ (D + S + A) Ṽ t (B.1)

Numerical experimentsIn this section, we assess of the performance of data-driven srhinkage rules under various numerical experiments involving Gaussian, Gamma and Poisson measurements.

Gaussian distributed measurements

We first consider the case of Gaussian measurements, where Y = X + W with E[W ij ] = 0, Var(W ij ) = τ 2 with τ = 1 √ m . In the following numerical experiments, we study the behavior of the spectral estimator:

• SURE-driven soft-thresholding

• Asymptotically optimal shrinkage proposed in [START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF] and [START_REF] Gavish | Optimal shrinkage of singular values[END_REF]

• SURE-driven weighted estimator that we have derived in Section 2.3

σk ũk ṽt k with σk = 1 -

where r ∈ [1, min(n, m)], and for the soft-thresholding, the value λ(Y ) > 0 is obtained by a numerical solver in order to minimize the SURE. Otherwise specified, we consider r = max{k ; σk > c n,m + }, i.e., an estimator of the rank using knowledge of the bulk edge c + ≈ c n,m + , hence, 1 1 {k ≤ r} 1 1 {σk>c n,m + } . As discussed in Section 4, we compare, in these experiments, the influence of rank estimation by analyzing the performances of the same estimators when either r = r max = min(n, m) (i.e. without knowledge the bulk edge, namely c + = 0), r = r * (oracle/true rank), r = r eff (effective rank [START_REF] Nadakuditi | OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[END_REF]) or by (4.8) (from hard-thresholding of singular values in [START_REF] Gavish | The optimal hard threshold for singular values is \(4/\sqrt {3}\)[END_REF]).

In order to assess the quality of SURE as an estimator of the MSE, we also compare the aforementioned approach with their oracle counterparts given by

where λ oracle (Y ) minimizes the squared error SE (non-expected risk) over the sets and softthresholding approximations respectively. Note that X r w and X soft are ideal approximations of X that cannot be used in practice but serve as benchmarks to evaluate the performances of the data-driven estimators Xr , Xsoft , Xr * and Xr w . In order to shed some light on the variance of these estimators, and indirectly on the variance of the SURE, we perform this experiments over M = 1000 independent realizations of Y .

The results are reported on Figure 5. For an estimator of the rank given either by r = max{k ; σk > c n,m + } (knowledge of the bulk edege), r = r * (oracle/true rank), r = r eff (effective rank) or by (4.8), it can be observed that Xr w , Xr * and X r w achieve comparable performances for all r ∈ [1, min(m, n)] even though the two first do not rely on the unknown matrix X. Similarly Xsoft and X soft achieve also comparable performances showing again that the SURE accurately estimates the MSE. In terms of error bands for the NMSE, Xr w , Xr * and X r w outperform Xsoft and X soft provided that r is large enough. Moreover, the performance of Xr w plateaus to its optimum when the rank r becomes large. This allows us to choose r = min(n, m) when we do not have a priori on the true or effective rank.

Interestingly, Fig. 5.(k) shows that when the above estimators are used without the knowledge of the bulk edge (i.e. by taking c n,m + = 0 in their computation instead of c n,m + = 1 + n m , which corresponds to the choice r = r max = min(n, m)), the performance of Xr w actually decreases when the rank r becomes too large. Indeed, it is clear from Fig. 5.(k), that the the error band of the NMSE of Xr w becomes much larger as the rank r increases. This illustrates that the SURE suffers from estimation variance in the case of over parametrization when r becomes too large, and thus it cannot be used to estimate jointly a too large number of weights. Therefore, the knowledge of an appropriate estimator r of the rank (e.g. using the bulk edge) seems to provide a relevant upper bound on the number of weights that can be jointly and robustly estimated with the SURE.

Gamma and Poisson measurements

Let us now consider the case where the entries of Y ij > 0 of the data matrix Y are independently sampled from a Gamma or Poisson distribution with mean X ij > 0. We again consider estimators of the form (2.17). In this context, we compare the following spectral shrinkage estimators, set for ε = 10 -6 , as:

• GSURE/SUKLS/PURE/SUKLA driven soft-thresholding

of Xsoft based on PUKLA. In fact, for other choices of matrices X, Xr w based on PUKLA might improve, in terms of MKLS, much more on Xsoft , and might improve not as much on Xr w based on PURE. Nevertheless, whatever X, we observed that Xr w driven by PUKLA always reaches at least as good performance in terms of MKLS as the best of Xr w driven by SE and Xsoft . Xr w actually decreases when the rank r becomes too large. As in the Gaussian setting, this can be explained by the fact that the GSURE, SUKLS, PURE and PUKLA suffer from estimation variance in the case of over parametrization, hence, they cannot be used to estimate jointly a too large number of weights. The active set s (in the same manner as the bulk edge) seems to provide a relevant selection of the weights that can be jointly and robustly estimated in a data driven way.

Signal matrix with equal singular values and increasing rank

We finally propose to highlight potential limitations of our approach in the situation where the rank r * of the matrix X = r * k=1 σ k u k v t k is let growing and all positive singular values σ k of X are equal, namely

where u k ∈ R n and v k ∈ R m are vectors with unit norm that are fixed, c n,m = n m and W is centered random matrix whose entries are iid Gaussian variables with variance τ 2 = 1/m. We again choose to fix n = 100 and m = 200, while the true rank is r * let growing from 1 to min(n, m) in the following numerical experiments. The constant γ is chosen to be larger than 1. Hence, eq. (5.1) corresponds to the Gaussian spiked population model in the setting where all positive singular values are equal and larger than the threshold c . All these estimators have been introduced in Section 5.2.

In Figure 8, we report the results of numerical experiments by displaying errors bars of the NMSE of these estimators as functions of the true rank r * . For low values of the true rank (r * ≤ 20), the data-driven estimators Xrmax w (our approach) and Xrmax * (shrinkage rule from [START_REF] Gavish | Optimal shrinkage of singular values[END_REF]) achieve the best performances that are similar in term of median value of the NMSE. However, our approach has some limitations with respect to the performances of the estimator from [START_REF] Gavish | Optimal shrinkage of singular values[END_REF] or data-driven soft-thresholding [START_REF] Candès | Unbiased risk estimates for singular value thresholding and spectral estimators[END_REF] in the setting where the signal matrix has equal positive singular values and when its rank is increasing. Moreover, the error bands of the NMSE for our approach becomes significantly larger than those of the other data-driven estimators when the true rank r * increases. This illustrates that SURE minimization leads to estimators with a high variance in the case of over parametrization, that is, when there exists a large number of significant singular values in the signal matrix.

A Proof of the main results

A.1 Proof of Proposition 3.2

Let us first introduce some notation and definitions to be used in the proof. For all 1 ≤ ≤ n, let λ be the eigenvalues of Y Y t namely λ = σ2 . For a fixed 1 ≤ k ≤ r * such that σ k > c 1/4 , let us introduce the complex-valued function g k defined by

where supp(µ k ) = λ ; 1 ≤ ≤ n, = k is the support of the random measure µ k = 1 n n =1; =k δ λ on R + , where δ λ denotes the Dirac measure at λ. It is clear that

The main difficulty in the proof is to show that, almost surely,

which is the purpose of what follows. For a matrix A ∈ R n×m (with n ≤ m), we denote its singular values by σ 1 (A) ≥ σ 2 (A) ≥ . . . ≥ σ n (A) ≥ 0. Hence, one has that σ = σ (Y ) for all 1 ≤ ≤ n. Now, we recall that Y = X + W where X is a fixed matrix of rank r * and W is a random matrix with iid entries sampled from a Gaussian distribution with zero mean and variance 1 m . The first step in the proof is to show that the random measure µ k behaves asymptotically as the almost sure limit of the empirical spectral measure µ W W t of the Wishart matrix W W t . By definition, the eigenvalues of W W t are λ (W ) = σ 2 (W ) for all 1 ≤ ≤ n and µ W W t is thus defined as

It is well know (see e.g. Theorem 3.6 in [START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF]) that, once m = m n ≥ n and lim n→+∞ n m = c with 0 < c ≤ 1, then, almost surely, the empirical spectral measure µ W W t converges weakly to the so-called Marchenko-Pastur distribution µ M P which is deterministic and has the following

We recall that such a convergence can also be characterized through the so-called Cauchy or Stieltjes transform which is defined for any probability measure µ on R as ∀z ∈ C outside the support of µ, g µ (z) = 1 zλ dµ(λ).

By eq. (3.3.2) in [START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF], one obtains that, almost surely,

where g M P is the Cauchy transform of µ M P and

Moreover, by Proposition 6 in [START_REF] Pastur | Matrices aléatoires: statistique asymptotique des valeurs propres[END_REF], the convergence (A.1) is uniform over any compact subset of C \ R.

Then, it follows from the so-called Weyl's interlacing inequalities (see e.g. Theorem 3.1.2 in [START_REF] Horn | Topics in matrix analysis[END_REF]) that for all 1 ≤ ≤ n

Thanks to the results that have been recalled above on the asymptotic properties of µ W W t , one may use inequalities (A.2) to prove that, almost surely, the random measure µ k converges weakly to the Marchenko-Pastur distribution µ M P . Under the assumptions of Proposition 3.2 and using Proposition 3.1, it can be shown that there exists η k > 0 such that, almost surely and for all sufficiently large n λ Therefore, using the uniform convergence (A.3) of g k to g M P and the continuity of g M P at z = ρ 2 (σ k ), one obtains that, almost surely,