Sharp oracle inequalities for aggregation of affine estimators

In the statistical literature, to the best of our knowledge, theoretical foundations of aggregation procedures were first studied by Nemirovski (Nemirovski Nemirovski00 , Juditsky and Nemirovski JuditskyNemirovski00 ) and independently by a series of papers by Catoni (see Catoni04 for an account) and Yang Yang00 , Yang03 , Yang04a . For the regression model, a significant progress was achieved by Tsybakov Tsybakov03 with introducing the notion of optimal rates of aggregation and proposing aggregation-rate-optimal procedures for the tasks of linear, convex and model selection aggregation. This point was further developed in Lounici07 , RigolletTsybakov07 , BuneaTsybakovWegkamp07 , especially in the context of high dimension with sparsity constraints and in Rigollet09 for Kullback–Leibler aggregation. However, it should be noted that the procedures proposed in Tsybakov03 that provably achieve the lower bounds in convex and linear aggregation require full knowledge of design distribution. This limitation was overcome in the recent work Wang2011 .

From a practical point of view, an important limitation of the previously cited results on aggregation is that they are valid under the assumption that the aggregated procedures are deterministic (or random, but independent of the data used for aggregation). The generality of those results—almost no restriction on the constituent estimators—compensates to this practical limitation.

In the Gaussian sequence model, a breakthrough was reached by Leung and Barron LeungBarron06 . Building on very elegant but not very well-known results by George George86a ²²2Corollary 2 in George86a coincides with Theorem 1 in LeungBarron06 in the case of exponential weights with temperature $\beta=2\sigma^{2}$; cf. equation (4) below for a precise definition of exponential weights. Furthermore, to the best of our knowledge, George86a is the first reference using the Stein lemma for evaluating the expected risk of the exponentially weighted aggregate., they established sharp oracle inequalities for the exponentially weighted aggregate (EWA) for constituent estimators obtained from the data vector by orthogonally projecting it on some linear subspaces. Dalalyan and Tsybakov DalalyanTsybakov07 , DalalyanTsybakov08 showed the result of LeungBarron06 remains valid under more general (non-Gaussian) noise distributions and when the constituent estimators are independent of the data used for the aggregation. A natural question arises whether a similar result can be proved for a larger family of constituent estimators containing projection estimators and deterministic ones as specific examples. The main aim of the present paper is to answer this question by considering families of affine estimators.

Our interest in affine estimators is motivated by several reasons. First, affine estimators encompass many popular estimators such as smoothing splines, the Pinsker estimator Pinsker80 , EfromovichPinsker96 , local polynomial estimators, nonlocal means BuadesCollMorel05 , SalmonLepennec09b , etc. For instance, it is known that if the underlying (unobserved) signal belongs to a Sobolev ball, then the (linear) Pinsker estimator is asymptotically minimax up to the optimal constant, while the best projection estimator is only rate-minimax. A second motivation is that—as proved by Juditsky and Nemirovski JuditskyNemirovski09 —the set of signals that are well estimated by linear estimators is very rich. It contains, for instance, sampled smooth functions, sampled modulated smooth functions and sampled harmonic functions. One can add to this set the family of piecewise constant functions as well, as demonstrated in PolzehlSpokoiny00 , with natural application in magnetic resonance imaging. It is worth noting that oracle inequalities for penalized empirical risk minimizer were also proved by Golubev Golubev10 , and for model selection by Arlot and Bach ArlotBach09 , Baraud, Giraud and Huet BaraudGiraudHuet10 .

In the present work, we establish sharp oracle inequalities in the model of heteroscedastic regression, under various conditions on the constituent estimators assumed to be affine functions of the data. Our results provide theoretical guarantees of optimality, in terms of expected loss, for the exponentially weighted aggregate. They have the advantage of covering in a unified fashion the particular cases of frozen estimators considered in DalalyanTsybakov08 and of projection estimators treated in LeungBarron06 .

We focus on the theoretical guarantees expressed in terms of oracle inequalities for the expected squared loss. Interestingly, although several recent papers ArlotBach09 , BaraudGiraudHuet10 , GoldenshlugerLepski08 discuss the paradigm of competing against the best linear procedure from a given family, none of them provide oracle inequalities with leading constant equal to one. Furthermore, most existing results involve some constants depending on different parameters of the setup. In contrast, the oracle inequality that we prove herein is with leading constant one and admits a simple formulation. It is established for (suitably symmetrized, if necessary) exponentially weighted aggregates George86a , Catoni04 , DalalyanTsybakov07 with an arbitrary prior and a temperature parameter which is not too small. The result is nonasymptotic but leads to an asymptotically optimal residual term when the sample size, as well as the cardinality of the family of constituent estimators, tends to infinity. In its general form, the residual term is similar to those obtained in the PAC-Bayes setting Mcallester98 , Langford02 , Seeger03 in that it is proportional to the Kullback–Leibler divergence between two probability distributions.

The problem of competing against the best procedure in a given family was extensively studied in the context of online learning and prediction with expert advice KivinenWarmuth99 , Cesa-BianchiLugosi06 . A connection between the results on online learning and statistical oracle inequalities was established by Gerchinovitz Gerchinovitz11 .

Throughout this work, we focus on the heteroscedastic regression model with Gaussian additive noise. We assume we are given a vector $\mathbf{Y}=(y_{1},\ldots,y_{n})^{\top}\in\mathbb{R}^{n}$ obeying the model

where $\bolds{\xi}=(\xi_{1},\ldots,\xi_{n})^{\top}$ is a centered Gaussian random vector, $f_{i}=\mathbf{f}(x_{i})$ where $\mathbf{f}\dvtx\mathcal{X}\rightarrow\mathbb{R}$ is an unknown function and $x_{1},\ldots,x_{n}\in\mathcal{X}$ are deterministic points. Here, no assumption is made on the set $\mathcal{X}$. Our objective is to recover the vector $\mathbf{f}=(f_{1},\ldots,f_{n})^{\top}$, often referred to as signal, based on the data $y_{1},\ldots,y_{n}$. In our work, the noise covariance matrix $\Sigma=\mathbb{E}[\bolds{\xi}\bolds{\xi}^{\top}]$ is assumed to be finite with a known upper bound on its spectral norm $|\!|\!|{\Sigma}|\!|\!|$. We denote by $\langle\cdot|\cdot\rangle_{n}$ the empirical inner product in $\mathbb{R}^{n}$: $\langle\mathbf{u}|\mathbf{v}\rangle_{n}=(1/n)\sum_{i=1}^{n}u_{i}v_{i}$. We measure the performance of an estimator $\hat{\mathbf{f}}$ by its expected empirical quadratic loss: $r=\mathbb{E}[\|\mathbf{f}-\hat{\mathbf{f}}\|_{n}^{2}]$ where $\|\mathbf{f}-\hat{\mathbf{f}}\|_{n}^{2}=\frac{1}{n}\sum_{i=1}^{n}(f_{i}-\hat{f}_{i})^{2}$.

We only focus on the task of aggregating affine estimators $\hat{\mathbf{f}}_{\lambda}$ indexed by some parameter $\lambda\in\Lambda$. These estimators can be written as affine transforms of the data $\mathbf{Y}=(y_{1},\ldots,y_{n})^{\top}\in\mathbb{R}^{n}$. Using the convention that all vectors are one-column matrices, we have $\hat{\mathbf{f}}_{\lambda}=A_{\lambda}\mathbf{Y}+\mathbf{b}_{\lambda}$, where the $n\times n$ real matrix $A_{\lambda}$ and the vector $\mathbf{b}_{\lambda}\in\mathbb{R}^{n}$ are deterministic. It means the entries of $A_{\lambda}$ and $\mathbf{b}_{\lambda}$ may depend on the points $x_{1},\ldots,x_{n}$ but not on the data $\mathbf{Y}$. Let us describe now different families of linear and affine estimators successfully used in the statistical literature. Our results apply to all these families, leading to a procedure that behaves nearly as well as the best (unknown) one of the family.

Ordinary least squares. Let $\{\mathcal{S}_{\lambda}\dvtx\lambda\in\Lambda\}$ be a set of linear subspaces of $\mathbb{R}^{n}$. A well-known family of affine estimators, successfully used in the context of model selection BarronBirgeMassart99 , is the set of orthogonal projections onto $\mathcal{S}_{\lambda}$. In the case of a family of linear regression models with design matrices $X_{\lambda}$, one has $A_{\lambda}=X_{\lambda}(X_{\lambda}^{\top}X_{\lambda})^{+}X_{\lambda}^{\top}$, where $(X_{\lambda}^{\top}X_{\lambda})^{+}$ stands for the Moore–Penrose pseudo-inverse of $X_{\lambda}^{\top}X_{\lambda}$.

Diagonal filters. Other common estimators are the so-called diagonal filters corresponding to diagonal matrices $A=\operatorname{diag}(a_{1},\ldots,a_{n})$. Examples include the following:

• •

  Ordered projections: $a_{k}=\mathbh{1}_{(k\leq\lambda)}$ for some integer $\lambda$ [$\mathbh{1}_{(\cdot)}$ is the indicator function]. Those weights are also called truncated SVD (Singular Value Decomposition) or spectral cutoff. In this case a natural parametrization is $\Lambda=\{1,\ldots,n\}$, indexing the number of elements conserved.

• •

  Block projections: $a_{k}=\mathbh{1}_{(k\leq w_{1})}+\sum_{j=1}^{m-1}\lambda_{j}\mathbh{1}_{(w_{j}\leq k\leq w_{j+1})}$, $k=1,\ldots,n$, where $\lambda_{j}\in\{0,1\}$. Here the natural parametrization is $\Lambda=\{0,1\}^{m-1}$, indexing subsets of $\{1,\ldots,m-1\}$.

• •

  Tikhonov–Philipps filter: $a_{k}=\frac{1}{1+(k/w)^{\alpha}}$, where $w,\alpha>0$. In this case, $\Lambda=(\mathbb{R}_{+}^{*})^{2}$, indexing continuously the smoothing parameters.

• •

  Pinsker filter: $a_{k}=(1-\frac{k^{\alpha}}{w})_{+}$, where $x_{+}=\max(x,0)$ and $(w,\alpha)=\lambda\in\Lambda=(\mathbb{R}_{+}^{*})^{2}$.

Kernel ridge regression. Assume that we have a positive definite kernel $k\dvtx\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ and we aim at estimating the true function $f$ in the associated reproducing kernel Hilbert space ($\mathcal{H}_{k},\|\cdot\|_{k}$). The kernel ridge estimator is obtained by minimizing the criterion $\|\mathbf{Y}-\mathbf{f}\|^{2}_{n}+\lambda\|\mathbf{f}\|_{k}^{2}$ w.r.t. $f\in\mathcal{H}_{k}$ (see Shawe-TaylorChristianini00 , page 118). Denoting by $K$ the $n\times n$ kernel-matrix with element $K_{i,j}=k(x_{i},x_{j})$, the unique solution $\hat{\mathbf{f}}$ is a linear estimate of the data, $\hat{\mathbf{f}}=A_{\lambda}\mathbf{Y}$, with $A_{\lambda}=K(K+n\lambda{I_{n\times n}})^{-1}$, where ${I_{n\times n}}$ is the $n\times n$ identity matrix.

Multiple Kernel learning. As described in ArlotBach09 , it is possible to handle the case of several kernels $k_{1},\ldots,k_{M}$, with associated positive definite matrices $K_{1},\ldots,K_{M}$. For a parameter $\lambda=(\lambda_{1},\ldots,\lambda_{M})\in\Lambda=\mathbb{R}_{+}^{M}$, one can define the estimators $\hat{\mathbf{f}}_{\lambda}=A_{\lambda}\mathbf{Y}$ with

It is worth mentioning that the formulation in equation (2) can be linked to the group Lasso YuanLin06 and to the multiple kernel learning introduced in LanckrietCristianiniBartlettElGhaouiJordan03 —see ArlotBach09 for more details.

Moving averages. If we think of coordinates of $\mathbf{f}$ as some values assigned to the vertices of an undirected graph, satisfying the property that two nodes are connected if the corresponding values of $\mathbf{f}$ are close, then it is natural to estimate $f_{i}$ by averaging out the values $Y_{j}$ for indices $j$ that are connected to $i$. The resulting estimator is a linear one with a matrix $A=(a_{ij})_{i,j=1}^{n}$ such that $a_{ij}=\mathbh{1}_{V_{i}}(j)/n_{i}$, where $V_{i}$ is the set of neighbors of the node $i$ in the graph and $n_{i}$ is the cardinality of $V_{i}$.

In Section 2 we introduce EWA and state a PAC-Bayes type bound in expectation assessing optimality properties of EWA in combining affine estimators. The strengths and limitations of the results are discussed in Section 3. The extension of these results to the case of grouped aggregation—in relation with ill-posed inverse problems—is developed in Section 4. As a consequence, we provide in Section 5 sharp oracle inequalities in various setups: ranging from finite to continuous families of constituent estimators and including sparse scenarii. In Section 6 we apply our main results to prove that combining Pinsker’s type filters with EWA leads to asymptotically sharp adaptive procedures over Sobolev ellipsoids. Section 7 is devoted to numerical comparison of EWA with other classical filters (soft thresholding, blockwise shrinking, etc.) and illustrates the potential benefits of aggregating. The conclusion is given in Section 8, while the proofs of some technical results (Propositions 2–6) are provided in the supplementary material DalalyanSalmonsupp .

Let $r_{\lambda}=\mathbb{E}[\|\hat{\mathbf{f}}_{\lambda}-\mathbf{f}\|_{n}^{2}]$ denote the risk of the estimator $\hat{\mathbf{f}}_{\lambda}$, for any $\lambda\in\Lambda$, and let $\hat{r}_{\lambda}$ be an estimator of $r_{\lambda}$. The precise form of $\hat{r}_{\lambda}$ strongly depends on the nature of the constituent estimators. For any probability distribution $\pi$ over $\Lambda$ and for any $\beta>0$, we define the probability measure of exponential weights, $\hat{\pi}$, by

The corresponding exponentially weighted aggregate, henceforth denoted by ${\hat{\mathbf{f}}}_{\mathrm{EWA}}$, is the expectation of $\hat{\mathbf{f}}_{\lambda}$ w.r.t. the probability measure $\hat{\pi}$:

We will frequently use the terminology of Bayesian statistics: the measure $\pi$ is called prior, the measure $\hat{\pi}$ is called posterior and the aggregate ${\hat{\mathbf{f}}}_{\mathrm{EWA}}$ is then the posterior mean. The parameter $\beta$ will be referred to as the temperature parameter. In the framework of aggregating statistical procedures, the use of such an aggregate can be traced back to George George86a .

The interpretation of the weights $\theta(\lambda)$ is simple: they up-weight estimators all the more that their performance, measured in terms of the risk estimate $\hat{r}_{\lambda}$, is good. The temperature parameter reflects the confidence we have in this criterion: if the temperature is small ($\beta\approx 0$), the distribution concentrates on the estimators achieving the smallest value for $\hat{r}_{\lambda}$, assigning almost zero weights to the other estimators. On the other hand, if $\beta\rightarrow+\infty$, then the probability distribution over $\Lambda$ is simply the prior $\pi$, and the data do not influence our confidence in the estimators.

In this paper we only focus on affine estimators

where the $n\times n$ real matrix $A_{\lambda}$ and the vector $\mathbf{b}_{\lambda}\in\mathbb{R}^{n}$ are deterministic. Furthermore, we will assume that an unbiased estimator $\widehat{\Sigma}$ of the noise covariance matrix $\Sigma$ is available. It is well known (cf. Appendix for details) that the risk of the estimator (6) is given by

and that $\hat{r}_{\lambda}^{\mathrm{unb}}$, defined by

is an unbiased estimator of $r_{\lambda}$. Along with $\hat{r}_{\lambda}^{\mathrm{unb}}$, we will use another estimator of the risk that we call the adjusted risk estimate and define by

One can notice that the adjusted risk estimate $\hat{r}_{\lambda}^{\mathrm{adj}}$ coincides with the unbiased risk estimate $\hat{r}_{\lambda}^{\mathrm{unb}}$ if and only if the matrix $A_{\lambda}$ is an orthogonal projector.

To state our main results, we denote by $\mathcal{P}_{\Lambda}$ the set of all probability measures on $\Lambda$ and by $\mathcal{K}(p,p^{\prime})$ the Kullback–Leibler divergence between two probability measures $p,p^{\prime}\in\mathcal{P}_{\Lambda}$:

We write $S_{1}\preceq S_{2}$ (resp., $S_{1}\succeq S_{2}$) for two symmetric matrices $S_{1}$ and $S_{2}$, when $S_{2}-S_{1}$ (resp., $S_{1}-S_{2}$) is semi-definite positive.

The simplest setting in which all the conditions of part (i) of Theorem 1 are fulfilled is when the matrices $A_{\lambda}$ and $\Sigma$ are all diagonal, or diagonalizable in a common base. This result, as we will see in Section 6, leads to a new estimator which is adaptive, in the exact minimax sense, over the collection of all Sobolev ellipsoids. It also suggests a new method for efficiently combining varying-block-shrinkage estimators, as described in Section 5.4.

However, part (i) of Theorem 1 leaves open the issue of aggregating affine estimators defined via noncommuting matrices. In particular, it does not allow us to evaluate the MSE of EWA when each $A_{\lambda}$ is a convex or linear combination of a fixed family of projection matrices on nonorthogonal linear subspaces. These kinds of situations may be handled via the result of part (ii) of Theorem 1. One can observe that in the particular case of a finite collection of projection estimators (i.e., $A_{\lambda}=A_{\lambda}^{2}$ and $\mathbf{b}_{\lambda}=0$ for every $\lambda$), the result of part (ii) offers an extension of LeungBarron06 , Corollary 6, to the case of general noise covariances (LeungBarron06 deals only with i.i.d. noise).

An important situation covered by part (ii) of Theorem 1, but not by part (i), concerns the case when signals of interest $\mathbf{f}$ are smooth or sparse in a basis $\mathcal{B}_{\mathrm{sig}}$ which is different from the basis $\mathcal{B}_{\mathrm{noise}}$ orthogonalizing the covariance matrix $\Sigma$. In such a context, one may be interested in considering matrices $A_{\lambda}$ that are diagonalizable in the basis $\mathcal{B}_{\mathrm{sig}}$ which, in general, do not commute with $\Sigma$.

The results presented so far concern the situation when the matrices $A_{\lambda}$ are symmetric. However, using the last part of Theorem 1, it is possible to propose an estimator of $\mathbf{f}$ that is almost as accurate as the best affine estimator $A_{\lambda}\mathbf{Y}+\mathbf{b}_{\lambda}$ even if the matrices $A_{\lambda}$ are not symmetric. Interestingly, the estimator enjoying this property is not obtained by aggregating the original estimators $\hat{\mathbf{f}}_{\lambda}=A_{\lambda}\mathbf{Y}+\mathbf{b}_{\lambda}$ but the “symmetrized” estimators $\tilde{\mathbf{f}}_{\lambda}=\tilde{A}_{\lambda}\mathbf{Y}+\mathbf{b}_{\lambda}$, where $\tilde{A}_{\lambda}=A_{\lambda}+A_{\lambda}^{\top}-A_{\lambda}^{\top}A_{\lambda}$. Besides symmetry, an advantage of the matrices $\tilde{A}_{\lambda}$, as compared to the $A_{\lambda}$’s, is that they automatically satisfy the contraction condition $\tilde{A}_{\lambda}\preceq I_{n\times n}$ required by part (ii) of Theorem 1. We will refer to this method as Symmetrized Exponentially Weighted Aggregates (or SEWA) DalalyanSalmon11b .

To understand the scope of condition (C), let us present several cases of widely used linear estimators for which this condition is satisfied:

• •

  The simplest class of matrices $A_{\lambda}$ for which condition (C) holds true are orthogonal projections. Indeed, if $A_{\lambda}$ is a projection matrix, it satisfies $A_{\lambda}^{\top}A_{\lambda}=A_{\lambda}$ and, therefore, $\operatorname{Tr}(\widehat{\Sigma}A_{\lambda})=\operatorname{Tr}(\widehat{\Sigma}A_{\lambda}^{\top}A_{\lambda})$.

• •

  When the matrix $\widehat{\Sigma}$ is diagonal, then a sufficient condition for (C) is $a_{ii}\leq\sum_{j=1}^{n}a_{ji}^{2}$. Consequently, (C) holds true for matrices having only zeros on the main diagonal. For instance, the $k$NN filter in which the weight of the observation $Y_{i}$ is replaced by zero, that is, $a_{ij}=\mathbf{1}_{j\in\{j_{i,1},\ldots,j_{i,k}\}}/k$ satisfies this condition.

• •

  Under a little bit more stringent assumption of homoscedasticity, that is, when $\widehat{\Sigma}=\widehat{\sigma}^{2}I_{n\times n}$, if the matrices $A_{\lambda}$ are such that all the nonzero elements of each row are equal and sum up to one (or a quantity larger than one), then $\operatorname{Tr}(A_{\lambda})=\operatorname{Tr}(A_{\lambda}^{\top}A_{\lambda})$ and (C) is fulfilled. A notable example of linear estimators that satisfy this condition are Nadaraya–Watson estimators with rectangular kernel and nearest neighbor filters.

In some rare situations, the matrix $\Sigma$ is known and it is natural to use $\widehat{\Sigma}=\Sigma$ as an unbiased estimator. Besides this not very realistic situation, there are at least two contexts in which it is reasonable to assume that an unbiased estimator of $\Sigma$, independent of $\mathbf{Y}$, is available.

The first case corresponds to problems in which a signal can be recorded several times by the same device, or once but by several identical devices. For instance, this is the case when an object is photographed many times by the same digital camera during a short time period. Let $\mathbf{Z}_{1},\ldots,\mathbf{Z}_{N}$ be the available signals, which can be considered as i.i.d. copies of an $n$-dimensional Gaussian vector with mean $\mathbf{f}$ and covariance matrix $\Sigma_{Z}$. Then, defining $\mathbf{Y}=(\mathbf{Z}_{1}+\cdots+\mathbf{Z}_{N})/N$ and $\widehat{\Sigma}_{Z}=(N-1)^{-1}(\mathbf{Z}_{1}\mathbf{Z}_{1}^{\top}+\cdots+\mathbf{Z}_{N}\mathbf{Z}_{N}^{\top}-N\mathbf{Y}\mathbf{Y}^{\top})$, we find ourselves within the framework covered by previous theorems. Indeed, $\mathbf{Y}\sim\mathcal{N}_{n}(\mathbf{f},\Sigma_{Y})$ with $\Sigma_{Y}=\Sigma_{Z}/N$ and $\widehat{\Sigma}_{Y}=\widehat{\Sigma}_{Z}/N$ is an unbiased estimate of $\Sigma_{Y}$, independent of $\mathbf{Y}$. Note that our theory applies in this setting for every integer $N\geq 2$.

The second case is when the dominating part of the noise comes from the device which is used for recording the signal. In this case, the practitioner can use the device in order to record a known signal, $\mathbf{g}$. In digital image processing, $\mathbf{g}$ can be a black picture. This will provide a noisy signal $\mathbf{Z}$ drawn from Gaussian distribution $\mathcal{N}_{n}(\bm{g},\Sigma)$, independent of $\mathbf{Y}$ which is the signal of interest. Setting $\widehat{\Sigma}=(\mathbf{Z}-\mathbf{g})(\mathbf{Z}-\mathbf{g})^{\top}$, one ends up with an unbiased estimator of $\Sigma$, which is independent of $\mathbf{Y}$.

All the results stated in this work provide sharp nonasymptotic bounds on the expected risk of EWA. It would be insightful to complement this study by risk bounds that hold true with high probability. However, it was recently proved in DaiRigZha12 that EWA is deviation suboptimal: there exist a family of constituent estimators and a constant $C>0$ such that the difference between the risk of EWA and that of the best constituent estimator is larger than $C/\sqrt{n}$ with probability at least $0.06$. Nevertheless, several empirical studies (see, e.g., DaiZhang11 ) demonstrated that EWA has often a smaller risk than some of its competitors, such as the empirical star procedure Audibert07 , which are provably optimal in the sense of OIs with high probability. Furthermore, numerical experiments carried out in Section 7 show that the standard-deviation of the risk of EWA is of the order of $1/n$. This suggests that under some conditions on the constituent estimators it might be possible to establish OIs for EWA that are similar to (10) but hold true with high probability. A step in proving this kind of result was done in LecueMendelson10 , Theorem C, for the model of regression with random design.

The OI of the previous section requires various conditions on the constituent estimators $\hat{\mathbf{f}}_{\lambda}=A_{\lambda}\mathbf{Y}+\mathbf{b}_{\lambda}$. One may wonder how general these conditions are and is it possible to extend these OIs to more general $\hat{\mathbf{f}}_{\lambda}$’s. Although this work does not answer this question, we can sketch some elements of response.

First of all, we stress that the conditions of the present paper relax significantly those of previous results existing in statistical literature. For instance, Kneip Kneip94 considered only linear estimators, that is, $\mathbf{b}_{\lambda}\equiv 0$ and, more importantly, only ordered sets of commuting matrices $A_{\lambda}$. The ordering assumption is dropped in Leung and Barron LeungBarron06 , in the case of projection matrices. Note that neither of these assumptions is satisfied for the families of Pinsker and Tikhonov–Philipps estimators. The present work strengthens existing results in considering more general, affine estimators extending both projection matrices and ordered commuting matrices.

Despite the advances achieved in this work, there are still interesting cases that are not covered by our theory. We now introduce a family of estimators commonly used in image processing that do not satisfy our assumptions. In recent years, nonlocal means (NLM) became quite popular in image processing BuadesCollMorel05 . This method of signal denoising, shown to be tied in with EWA SalmonLepennec09b , removes noise by exploiting signals self-similarities. We briefly define the NLM procedure in the case of one-dimensional signals.

Assume that a vector $\mathbf{Y}=(y_{1},\ldots,y_{n})^{\top}$ given by (1) is observed with $f_{i}=\mathbf{f}(i/n)$, $i=1,\ldots,n$, for some function $\mathbf{f}\dvtx[0,1]\to\mathbb{R}$. For a fixed “patch-size” $k\in\{1,\ldots,n\}$, let us define $\mathbf{f}_{[i]}=(f_{i},f_{i+1},\ldots,f_{i+k-1})^{\top}$ and $\mathbf{Y}_{[i]}=(y_{i},y_{i+1},\ldots,y_{i+k-1})^{\top}$ for every $i=1,\ldots,n-k+1$. The vectors $\mathbf{f}_{[i]}$ and $\mathbf{Y}_{[i]}$ are, respectively, called true patch and noisy patch. The NLM consists in regarding the noisy patches $\mathbf{Y}_{[i]}$ as constituent estimators for estimating the true patch $\mathbf{f}_{[i_{0}]}$ by applying EWA. One easily checks that the constituent estimators $\mathbf{Y}_{[i]}$ are affine in $\mathbf{Y}_{[i_{0}]}$, that is, $\mathbf{Y}_{[i]}=A_{i}\mathbf{Y}_{[i_{0}]}+\mathbf{b}_{i}$ with $A_{i}$ and $\mathbf{b}_{i}$ independent of $\mathbf{Y}_{[i_{0}]}$. Indeed, if the distance between $i$ and $i_{0}$ is larger than $k$, then $\mathbf{Y}_{[i]}$ is independent of $\mathbf{Y}_{[i_{0}]}$ and, therefore, $A_{i}=0$ and $\mathbf{b}_{i}=\mathbf{Y}_{[i]}$. If $|i-i_{0}|<k$, then the matrix $A_{i}$ is a suitably chosen shift matrix and $\mathbf{b}_{i}$ is the projection of $\mathbf{Y}_{[i]}$ onto the orthogonal complement of the image of $A_{i}$. Unfortunately, these matrices $\{A_{i}\}$ and vectors $\{b_{i}\}$ do not fit our framework, that is, the assumption $A_{i}b_{i}=A_{i}^{\top}b_{i}=0$ is not satisfied.

Finally, our proof technique is specific to affine estimators. Its extension to estimators defined as a more complex function of the data will certainly require additional tools and is a challenging problem for future research. Yet, it seems unlikely to get sharp OIs with optimal remainder term for a fairly general family of constituent estimators (without data-splitting), since this generality inherently increases the risk of overfitting.

In order to demonstrate that inequality (14) can be reformulated in terms of an OI as defined by (3), let us consider the case when the prior $\pi$ is discrete, that is, $\pi(\Lambda_{0})=1$ for a countable set $\Lambda_{0}\subset\Lambda$, and w.l.o.g $\Lambda_{0}=\mathbb{N}$. Then, the following result holds true.

It suffices to apply Theorem 3 and to upper-bound the right-hand side by the minimum over all Dirac measures $p=\delta_{\ell}$ such that $\pi_{\ell}>0$. This inequality can be compared to Corollary 2 in BaraudGiraudHuet10 , Section 4.3. Our result has the advantage of having factor one in front of the expectation of the left-hand side, while in BaraudGiraudHuet10 a constant much larger than 1 appears. However, it should be noted that the assumptions on the (estimated) noise covariance matrix are much weaker in BaraudGiraudHuet10 .

It may be useful in practice to combine a family of affine estimators indexed by an open subset of $\mathbb{R}^{M}$ for some $M\in\mathbb{N}$ (e.g., to build an estimator nearly as accurate as the best kernel estimator with fixed kernel and varying bandwidth). To state an oracle inequality in such a “continuous” setup, let us denote by $d_{2}(\bolds{\lambda},\partial\Lambda)$ the largest real $\tau>0$ such that the ball centered at $\bolds{\lambda}$ of radius $\tau$—hereafter denoted by $B_{\bolds{\lambda}}(\tau)$—is included in $\Lambda$. Let $\operatorname{Leb}(\cdot)$ be the Lebesgue measure in $\mathbb{R}^{M}$.

It suffices to apply assertion (i) of Theorem 1 and to upper-bound the right-hand side in inequality (10) by the minimum over all measures having as density $p_{\bolds{\lambda}^{*},\tau^{*}}(\bolds{\lambda})=\mathbh{1}_{B_{\bolds{\lambda}^{*}}(\tau^{*})}(\bolds{\lambda})/\operatorname{Leb}(B_{\bolds{\lambda}^{*}}(\tau^{*}))$. Choosing $\tau^{*}=\min(n^{-1},d_{2}(\bolds{\lambda}^{*},\partial\Lambda))$ such that $B_{\bolds{\lambda}^{*}}(\tau^{*})\subset\Lambda$, the measure $p_{\bolds{\lambda}^{*},\tau^{*}}(\bolds{\lambda})\,d\bolds{\lambda}$ is absolutely continuous w.r.t. the uniform prior $\pi$ and the Kullback–Leibler divergence between these two measures equals $\log\{\operatorname{Leb}(\Lambda)/\operatorname{Leb}(B_{\bolds{\lambda}^{*}}(\tau^{*}))\}$. Using $\operatorname{Leb}(B_{\bolds{\lambda}^{*}}(\tau^{*}))\geq({2\tau^{*}}/{\sqrt{M}})^{M}$ and the Lipschitz condition, we get the desired inequality.

Note that it is not very stringent to require the risk function $r_{\bolds{\lambda}}$ to be Lipschitz continuous, especially since this condition needs not be satisfied uniformly in $\mathbf{f}$. Let us consider the ridge regression: for a given design matrix $X\in\mathbb{R}^{n\times p}$, $A_{\lambda}=X(X^{\top}X+\gamma_{n}\lambda I_{n\times n})^{-1}X^{\top}$ and $\mathbf{b}_{\lambda}=0$ with $\lambda\in[{\lambda_{*}},{\lambda^{*}}]$, $\gamma_{n}$ being a given normalization factor typically set to $n$ or $\sqrt{n}$, $\lambda_{*}>0$ and $\lambda^{*}\in[\lambda_{*},\infty]$. One can easily check the Lipschitz property of the risk function with $L_{r}=L_{r}(f)=4\lambda_{*}^{-1}\|\mathbf{f}\|_{n}^{2}+(2/n)\operatorname{Tr}(\Sigma)$.