Joseph Salmon

Arnak Dalalyan

Optimal aggregation of affine estimators

We consider the problem of combining a (possibly uncountably infinite) set of affine estimators in non-parametric regression model with heteroscedastic Gaussian noise. Focusing on the exponentially weighted aggregate, we prove a PAC-Bayesian type inequality that leads to sharp oracle inequalities in discrete but also in continuous settings. The framework is general enough to cover the combinations of various procedures-such as the least square regression, the kernel ridge regression, the shrinkage estimators, etc.-used in the literature on statistical inverse problems. As a consequence, we show that the proposed aggregate provides an adaptive estimator in the exact minimax sense without neither discretizing the range of tuning parameters nor splitting the set of observations. We also illustrate numerically the good performance achieved by the exponentially weighted aggregate.

Introduction

There is a growing empirical evidence of superiority of aggregated statistical procedures, also referred to as blending, stacked generalization, or ensemble methods, with respect to "pure" ones. Since their introduction in the 1990's, the most famous aggregation procedures such as Boosting [START_REF] Freund | Boosting a weak learning algorithm by majority[END_REF], Bagging [START_REF] Breiman | Bagging predictors[END_REF] or Random Forest [START_REF] Amit | Shape quantization and recognition with randomized trees[END_REF] were successfully used in practice for a large variety of applications. Moreover, the most recent Machine Learning competitions such as the Pascal VOC or the Netflix challenge were won by procedures combining different types of classifiers / predictors / estimators. It is therefore of central interest to understand from a theoretical point of view what kind of aggregation strategies should be used for getting the best possible combination of the available statistical procedures.

Historical remarks and motivation

In the statistical literature, to the best of our knowledge, the lecture notes of [START_REF] Nemirovski | Topics in non-parametric statistics[END_REF] was the first work concerned by the theoretical analysis of aggregation procedures. It was followed by a paper by [START_REF] Juditsky | Functional aggregation for nonparametric regression[END_REF], as well as by a series of papers by Catoni (see [START_REF] Catoni | Statistical learning theory and stochastic optimization[END_REF] for a comprehensive account) and [START_REF] Yang | Combining different procedures for adaptive regression[END_REF][START_REF] Yang | Regression with multiple candidate models: selecting or mixing?[END_REF][START_REF] Yang | Aggregating regression procedures to improve performance[END_REF]. For the regression model, a significant progress was achieved by [START_REF] Tsybakov | Optimal rates of aggregation[END_REF] 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 by [START_REF] Lounici | Generalized mirror averaging and D-convex aggregation[END_REF], [START_REF] Ph | Linear and convex aggregation of density estimators[END_REF], [START_REF] Lecué | Optimal rates of aggregation in classification under low noise assumption[END_REF], [START_REF] Bunea | Aggregation for Gaussian regression[END_REF], especially in the context of high dimension with sparsity constraints.

From a practical point of view, an important limitation of the previously cited results on the aggregation is that they are valid under the assumption that the aggregated procedures are deterministic (or random, but independent of the data used for the aggregation). In the Gaussian sequence model, a breakthrough was reached by [START_REF] Leung | Information theory and mixing least-squares regressions[END_REF]. Building on a very elegant but not very well known result of [START_REF] George | Minimax multiple shrinkage estimation[END_REF], they established sharp oracle inequalities for the exponentially weighted aggregate (EWA) under the condition that the aggregated estimators are obtained from the data vector by orthogonally projecting it on some linear subspaces. Dalalyan andTsybakov (2007, 2008), established the validity of Leung and Barron's result under more general (non Gaussian) noise distributions provided that 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 of all, affine estimators encompass many popular estimators such as the smoothing splines, the Pinsker estimator [START_REF] Pinsker | Optimal filtration of square-integrable signals in Gaussian noise[END_REF], [START_REF] Efromovich | Sharp-optimal and adaptive estimation for heteroscedastic nonparametric regression[END_REF], the local polynomial estimators, the non-local means [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], [START_REF] Salmon | NL-Means and aggregation procedures[END_REF], 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 [START_REF] Juditsky | Nonparametric denoising of signals with unknown local structure. I. Oracle inequalities[END_REF]-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 (cf. [START_REF] Juditsky | Nonparametric denoising of signals with unknown local structure. I. Oracle inequalities[END_REF] for precise definitions). It is worth noting that oracle inequalities for the penalized empirical risk minimizer were also established by [START_REF] Golubev | On universal oracle inequalities related to high-dimensional linear models[END_REF], and for the model selection by [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF], [START_REF] Baraud | Estimator selection in the gaussian setting[END_REF].

In the present work, we establish sharp oracle inequalities in the statistical model of heteroscedastic regression, under various conditions on the constituent estimators assumed to be affine functions of the data. We assume that the design is deterministic and that the noise is Gaussian with a given covariance matrix. Our results provide theoretical guarantees of optimality, in terms of the expected loss, for the exponentially weighted aggregate. They have the advantage of covering in a unified fashion the particular cases of deterministic estimators considered by [START_REF] Dalalyan | Aggregation by exponential weighting, sharp pac-bayesian bounds and sparsity[END_REF] and of projection estimators treated by [START_REF] Leung | Information theory and mixing least-squares regressions[END_REF].

Notation

Throughout this work, we focus on the heteroscedastic regression model with Gaussian additive noise. More precisely, we assume that we are given a vector Y = (y 1 , • • • , y n ) ∈ R n obeying the model:

y i = f i + ξ i , for i = 1, . . . , n, (1)
where ξ = (ξ 1 , . . . , ξ n ) is a centered Gaussian random vector, f i = f(x i ) where f is an unknown function X → R and x 1 , . . . , x n ∈ X are deterministic points. Here, no assumption is made on the set X .

Our objective is to recover the vector f = ( f 1 , . . . , f n ), often referred to as signal, based on the data y 1 , . . . , y n . In our work, the noise covariance matrix Σ = E[ξξ ] is assumed to be diagonal (so it can be written

Σ = diag(σ 2 1 , • • • , σ 2 n ))
, with a known upper bound on the spectral norm |Σ |. In our case,

|Σ | = max i =1,••• ,n σ 2 i .
We measure the performance of an estimator f by its expected empirical quadratic loss:

r = E( f -f 2 n ) where f -f 2 n = 1 n n i =1 ( f i -fi ) 2 .
We also denote by 〈•|•〉 n the corresponding empirical inner product.

In this paper, we only focus on affine estimators f λ , i.e., estimators that can be written as affine transforms of the data Y = (y 1 , . . . , y n ) ∈ R n . Using the convention that all vectors are one-column matrices, affine estimators can be defined by

f λ = A λ Y + b λ , (2)
where the n × n real matrix A λ and the vector b λ ∈ R n are deterministic. This means that the entries of A λ and b λ may depend on the points x 1 , . . . , x n but not on the data vector Y . It is well-known that the quadratic risk of the estimator (2) is given by

r λ = E f -f λ 2 n = (A λ -I n×n ) f + b λ 2 n + Tr(A λ ΣA λ ) n (3)
and that rλ , defined by

rλ = Y -f λ 2 n + 2 n Tr(ΣA λ ) - 1 n n i =1 σ 2 i (4)
is an unbiased estimator of r λ (direct application of Stein's Lemma, cf. Appendix).

Let us describe now different families of linear and affine estimators successfully used in the statistical literature (cf., for instance, [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF]). Our results apply to all these families and lead to a procedure that behaves nearly as well as the best one of the family.

Ordinary least squares

Let {S λ : λ ∈ Λ} be a set of linear subspaces of R n . A well known family of affine estimators, successfully used in the context of model selection by [START_REF] Barron | Risk bounds for model selection via penalization[END_REF], is the set of orthogonal projections onto S λ . In the case of a family of linear regression models with design matrices X λ , one has

A λ = X λ (X λ X λ ) -1 X λ .
Diagonal filters Another set of common estimators are the so called diagonal filters f = AY , where A is a diagonal matrix A = diag(a 1 , . . . , a n ). Popular examples include:

• Ordered projections : a k = 1l (k≤λ) for some integer λ (where 1l (•) is the indicator function). Those weights are also called truncated SVD or spectral cut-off. In this case the natural parametrization is Λ = {1, . . . , n}, indexing the number of elements conserved.

• Block projections: a k = 1l (k≤w 1 ) + m-1 j =1 λ j 1l (w j ≤k≤w j +1 ) , k = 1, . . . , n, where λ j ∈ {0, 1}. Here the natural parametrization is Λ = {0, 1} m-1 , indexing subsets of {1, m -1}.

• Tikhonov-Philipps filter: a k = 1 1+(k/w) α , where w, α > 0. The set Λ = (R * + ) 2 indexes continuously the smoothing parameters.

• Pinsker filter: a k = 1 -k α w + , where x + = max(x, 0) and w, α > 0. In this case also Λ = (R * + ) 2 .

Kernel ridge regression

Assume that we have a positive definite kernel k : X × X → R and we aim at estimating the true function f in the associated reproducing kernel Hilbert space (H k , • k ). The kernel ridge estimator is obtained by minimizing the criterion Y - (Shawe-Taylor and Cristianini, 2000, page 118)).Denoting by K the n ×n kernel-matrix with element

f 2 n + λ f 2 k w.r.t. f ∈ H k (see
K i , j = k(x i , x j ), the unique solution f is a linear estimate of the data, f = A λ Y , with A λ = K (K + nλI n×n ) -1
, where I n×n is the identity matrix of size n × n.

Multiple Kernel learning

As proposed in [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF], it is also possible to handle the case of several kernels k 1 , . . . , k M , with associated positive definite matrices K 1 , . . . , K M . For a parameter

λ = (λ 1 , . . . , λ M ) ∈ Λ = R M + one can define the estimators f λ = A λ Y with A λ = M m=1 λ m K m M m=1 λ m K m + nI n×n -1 . ( 5
)
It is worth mentioning that the formulation in Eq.( 5) can be linked to the group Lasso [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] and to the multiple kernel [START_REF] Lanckriet | Learning the kernel matrix with semidefinite programming[END_REF]) -see [START_REF] Bach | Consistency of the group lasso and multiple kernel learning[END_REF], [START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF] for more details.

Organization of the paper

In Section 2, we introduce EWA and state a PAC-Bayes type bound assessing the optimality of EWA in combining affine estimators. As a consequence, we provide in Section 3 sharp oracle inequalities in various set-ups: ranging from finite to continuous families of constituent estimators and including the sparsity scenario. In Section 4, we apply our main results to prove that combining Pinsker's type filters with EWA leads to an asymptotically sharp adaptive procedure over the Sobolev ellipsoids. Section 5 is devoted to a numerical comparison of EWA with other classical filters (soft thresholding, blockwise shrinking, etc.), and illustrates the potential benefits of the aggregation. Some concluding remarks are presented in Section 6, while technical proofs are postponed to the Appendix.

Aggregation of estimators: main result

In this section we describe the statistical framework for aggregating estimators and we also introduce the exponentially weighted aggregate. The task of aggregation consists in estimating f by a suitable combination of the elements of a family of constituent estimators

F Λ = ( f λ ) λ∈Λ ∈ R n .
The target objective of the aggregation is to build an aggregate f aggr , not necessarily in the family F Λ , that mimics the performance of the best constituent estimator. It is called oracle because of its dependence on the unknown function f . We assume that Λ is a measurable subset of R M , for some M ∈ N.

The theoretical tool commonly used for evaluating the quality of an aggregation procedure is the oracle inequality (OI), generally written in the following form:

E f aggr -f 2 n ≤ C n inf λ∈Λ E f λ -f 2 n + R n , (6)
with residual term R n tending to zero, and leading constant C n being bounded. The OIs with leading constant one are of central theoretical interest since they allow to bound the excess risk and to assess the aggregation-rate-optimality. The residual term R n depends on the complexity (size) of the family F Λ , as on the amount of noise, measured in term of variance in our context.

Exponentially Weighted Aggregate (EWA)

Let r λ = E( f λf 2 n ) denote the risk of the estimator f λ , for any λ ∈ Λ, and let rλ be an estimator of r λ . The precise form of rλ strongly depends on the nature of the constituent estimators. For any probability distribution π over the set Λ and for any β > 0, we define the probability measure of exponential weights, π, by the following formula:

π(d λ) = θ(λ)π(d λ) with θ(λ) = exp(-n rλ /β) Λ exp(-n rω /β)π(d ω) . ( 7
)
The corresponding exponentially weighted aggregate, henceforth denoted by f EWA , is the expectation of the f λ w.r.t. the probability measure π:

f EWA = Λ f λ π(d λ) . ( 8
)
It is convenient and customary to use the terminology of Bayesian statistics: the measure π is called prior, the measure π is called posterior and the aggregate f EWA is then the posterior mean. The parameter β 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 [START_REF] George | Minimax multiple shrinkage estimation[END_REF].

The interpretation of the weights θ(λ) is simple: they up-weight estimators all the more that their performance, measured in terms of the risk estimate rλ , is good. The temperature parameter reflects the confidence we have in this criterion: if the temperature is small (β ≈ 0) the distribution concentrates on the estimators achieving the smallest value for rλ , assigning almost zero weights to the other estimators. On the other hand, if β → +∞ then the probability distribution over Λ is simply the prior π, and the data do not modify our confidence in the estimators. It should also be noted that averaging w.r.t. the posterior π is not the only way of constructing an estimator of f , some alternative estimators based on π have been studied, see for instance [START_REF] Zhang | Information-theoretic upper and lower bounds for statistical estimation[END_REF][START_REF] Audibert | Fast learning rates in statistical inference through aggregation[END_REF].

Main result

To state our main result, we denote by P Λ the set of all probability measures on Λ and by K (p, p ) the Kullback-Leibler divergence between two probability measures p, p ∈ P Λ :

K (p, p ) = Λ log d p d p (λ) p(d λ) if p p , +∞ otherwise.
.

Theorem 1 (PAC Bayesian Bound) If either one of the following conditions is satisfied:

C 1 :
The matrices A λ are orthogonal projections (i.e., symmetric and idempotent) and the vectors b λ satisfy A λ b λ = 0, for all λ ∈ Λ.

C 2 : The matrices A λ are all symmetric, positive semidefinite and satisfy

A λ A λ = A λ A λ , A λ Σ = ΣA λ for all λ, λ ∈ Λ. All the vectors b λ are zero.
Then, the risk of the aggregate f EWA defined by Equations ( 7), ( 8) and ( 4) satisfies the inequality

r EWA = E( f EWA -f 2 n ) ≤ inf p∈P Λ Λ E f λ -f 2 n p(d λ) + β n K (p, π) (9) provided that β ≥ α |Σ |, where α = 4 if C 1 holds true and α = 8 if C 2 holds true.
All the proofs of our results are given in the appendix, at the end of the paper. Note also that the result of Theorem 1 applies to the estimator f EWA that uses the full knowledge of the covariance matrix Σ. Indeed, even if for the choice of β only an upper bound on the spectral norm of Σ is required, the entire matrix Σ enters in the definition of the unbiased risks rλ that is used for defining f EWA . The exponentially weighted aggregate f EWA is easily extended to handle the more realistic situation where an unbiased estimate Σ, independent of Y , of the covariance matrix Σ is available. Simply replace Σ by Σ in the definition of the unbiased risk estimate (4). When the estimators f λ satisfy π-a.e. condition C 1 or C 2 , choosing β = α | Σ |, it can be checked that a claim similar to Theorem 1 remains valid. Another observation is that using the extension of Stein's lemma presented in [START_REF] Dalalyan | Aggregation by exponential weighting, sharp pac-bayesian bounds and sparsity[END_REF], Lemma 1), a result similar to Theorem 1 can be established for some specific non Gaussian noise distributions, provided that the components of the noise vector are independent.

Sharp oracle inequalities

In this section, we discuss consequences of the main result for specific choices of prior measures. Some of them are closely related to the oracle inequalities presented in Dalalyan andTsybakov (2007, 2008), [START_REF] Alquier | Pac-bayesian bounds for sparse regression estimation with exponential weights[END_REF]Lounici (2010), Rigollet and[START_REF] Ph | Exponential screening and optimal rates of sparse estimation[END_REF] especially when dealing with the sparsity scenario in the high dimensional framework.

Discrete oracle inequality

In order to demonstrate that Inequality (9) can be reformulated in terms of an OI as defined by ( 6), let us consider the simple case when the prior π is discrete. That is, we assume that π(Λ 0 ) = 1 for a countable set Λ 0 ⊂ Λ. Without loss of generality, we assume that Λ 0 = N. Then, the following result holds true.

Proposition 1 If either one of the conditions C 1 and C 2 (cf. Theorem 1) is fulfilled and π is supported by N, then the aggregate f EWA defined by Equations ( 7), ( 8) and ( 4) satisfies the inequality

E( f EWA -f 2 n ) ≤ inf j ∈N:π j >0 E f j -f 2 n + β log(1/π j ) n (10) provided that β ≥ α |Σ |, where α = 4 if C 1 holds true and α = 8 if C 2 holds true.
Proof: It suffices to apply Theorem 1 and to bound the RHS from above by the minimum over all Dirac measures p = δ j with j such that π j > 0.

Continuous oracle inequality

It may be useful in practice to combine a family of affine estimators indexed by an open subset of R M , for some integer M > 0, for instance when the aim is to build an estimator that is nearly as accurate as the best kernel estimator with fixed kernel and varying bandwidth. In order to state an oracle inequality in such a "continuous" setup, let us denote by d 2 (λ, Λ) the largest real τ > 0 such that the ball centered at λ with radius τ is included in Λ. In what follows, Leb(•) stands for the Lebesgue measure.

Proposition 2 Let Λ ⊂ R M be an open and bounded set and let π be the uniform probability on Λ. Assume that the mapping

λ → r λ is Lipschitz continuous, i.e., |r λ -r λ | ≤ L r λ -λ 2 , ∀λ, λ ∈ Λ. Under the conditions C 1 or C 2 aggregate f EWA satisfies the inequality E f EWA -f 2 n ≤inf λ∈Λ E f λ -f 2 n + βM n log M 2 min(n -1 , d 2 (λ, Λ)) + L r + β log Leb(Λ) n . ( 11
)
for every β ≥ α |Σ | where α = 4 if C 1 holds true and α = 8 if C 2 holds true.

Sparsity oracle inequality

The continuous oracle inequality stated in previous subsection is well adapted to the case where the dimension M of Λ is small compared to the sample size n (or, more precisely, the signal to noise ratio n/ max i σ 2 i ). If this is not the case, the choice of the prior should be done more carefully. For instance, consider the case of a set Λ ⊂ R M with large M under the sparsity scenario: there is a sparse vector λ * ∈ Λ such that the risk of f λ * is small. Then, it is natural to choose a prior π that promotes the sparsity of λ. This can be done in the same vein as in Dalalyan andTsybakov (2007, 2008), by means of the heavy tailed prior:

π(d λ) ∝ M j =1 1 (1 + |λ j /τ| 2 ) 2 1l Λ (λ)d (λ), ( 12
)
where τ > 0 is a tuning parameter.

Proposition 3 Let Λ = R M and let π be defined by (12). Assume that the mapping λ → r λ is continuously differentiable and, for some M × M matrix M , satisfies:

r λ -r λ -∇r λ (λ -λ ) ≤ (λ -λ ) M (λ -λ ), ∀λ, λ ∈ Λ. ( 13
)
If either one of the conditions C 1 and C 2 (cf. Theorem 1) is fulfilled, then the aggregate f EWA defined by Equations ( 7), ( 8) and ( 4) satisfies the inequality

E f EWA -f 2 n ≤ inf λ∈R M E f λ -f 2 n + 4β n M j =1 log 1 + |λ j | τ + Tr(M )τ 2 (14) provided that β ≥ α |Σ |, where α = 4 if C 1 holds true and α = 8 if C 2 holds true.
Let us discuss here some consequences of this sparsity oracle inequality. First of all, let us remark that in most cases Tr(M ) is on the order of M and the choice τ = β/(nM ) ensures that the last term in the RHS of Eq. ( 14) decreases at the parametric rate 1/n. This is the choice we recommend for practical applications.

Assume now that we are given a large number of linear estimators ĝ 1 = G 1 Y , . . . , ĝ M = G M Y satisfying, for instance, condition C 2 . We will focus on matrices G j having a spectral norm bounded by one (it is well known that the failure of this condition makes the linear estimator inadmissible, cf. [START_REF] Cohen | All admissible linear estimates of the mean vector[END_REF]). Assume furthermore that our aim is to propose an estimator that mimics the behavior of the best possible convex combination of a pair of estimators chosen among ĝ 1 , . . . , ĝ M . This task can be accomplished in the framework of the present paper by setting Λ = R M and f λ = λ 1 ĝ 1 + . . . λ M ĝ M , where λ = (λ 1 , . . . , λ M ). If the collection { ĝ i } satisfies condition C 2 , then it is also the case for the collection of their linear combinations { f λ }. Moreover, the mapping λ → r λ is quadratic with the Hessian matrix ∇ 2 r λ given by the entries 2〈G j f |G j f 〉 n + 2 n Tr(G j ΣG j ), j , j = 1, . . . , M . This implies that Inequality (13) holds with M being the Hessian divided by 2. Therefore, setting σ = (σ 1 , . . . , σ n ), we get

Tr(M ) ≤ | M j =1 G 2 j |( f 2 n + σ 2 n ) ≤ M ( f 2 n + σ 2 n )
, where the norm of a matrix is understood as its largest singular value. Applying Proposition 3 with τ = β/(nM ), we get for

β ≥ 8 |Σ |, E f EWA -f 2 n ≤ inf α, j , j E α ĝ j + (1 -α) ĝ j -f 2 n + 8β n log 1 + Mn β + β n ( f 2 n + σ 2 n ), ( 15
)
where the inf is taken over all α ∈ [0, 1] and j , j ∈ {1, . . . , M }. This shows that, using EWA with a sufficiently large temperature, one can achieve the best possible risk over the convex combinations of a pair of linear estimators-selected from a large (but finite) family-at the price of a residual term that decreases at the parametric rate up to a log factor.

Oracle inequalities for varying-block-shrinkage estimators

Let us consider now the problem of aggregation of two-block shrinkage estimators. It means that the constituent estimators have the following form:

for λ = (a, b, k) ∈ [0, 1] 2 × {1, . . . , n} := Λ, f λ = A λ Y where A λ = diag a1l(i ≤ k) + b1l(i > k), i = 1, • • • , n . Let
E( f EWA -f 2 n ) ≤ E( fλ -f 2 n ) + β n 1 + log n 2 f 2 n + n Tr(Σ) 12β . ( 16
)
The proof of this result can be found in [START_REF] Dalalyan | Sharp oracle inequalities for aggregation of affine estimators[END_REF].

In the case Σ = I n×n , this result is comparable to [START_REF] Leung | Information Theory and Mixing Least Squares Regression[END_REF], page 20, Theorem 2.49), which states that in the model of homoscedastic regression (Σ = I n×n ), the EWA acting on two-block positive-part James-Stein shrinkage estimators satisfies, for any k = 3, • • • , n -3, and for β = 8, the oracle inequality

E( f Leung -f 2 n ) ≤ E( fλ -f 2 n ) + 9 n + 8 n min K >0 K ∨ log n -6 K -1 . ( 17
)

Application to minimax adaptive estimation

In the celebrated paper [START_REF] Pinsker | Optimal filtration of square-integrable signals in Gaussian noise[END_REF]proved that in the model (1) the minimax risk over ellipsoids can be asymptotically attained by a linear estimator. Let us denote by

θ k ( f ) = 〈 f |ϕ k 〉 n the coeffi- cients of the (orthogonal) discrete sine transform of f , hereafter denoted by D f . Pinsker's result- restricted to Sobolev ellipsoids F (α, R) = f ∈ R n : n k=1 k 2α θ k ( f ) 2 ≤ R and to the homoscedastic noise (Σ = σ 2 I n×n )-states that, as n → ∞, the equivalences inf f sup f ∈F (α,R) E f -f 2 n ∼ inf A sup f ∈F (α,R) E AY -f 2 n (18) ∼ inf w>0 sup f ∈F (α,R) E A α,w Y -f 2 n ( 19
)
hold (Tsybakov, 2009, Theorem 3.2), where the first inf is taken over all possible estimators f and A α,w = D diag (1k α /w) + ; k = 1, . . . , n D is the Pinsker filter in the discrete sine basis. In simple words, this implies that the (asymptotically) minimax estimator can be chosen from the quite narrow class of linear estimators with Pinsker's filter. However, it should be emphasized that the minimax linear estimator depends on the parameters α and R, that are generally unknown. An (adaptive) estimator, that does not depend on (α, R) and is asymptotically minimax over a large scale of Sobolev ellipsoids has been proposed by [START_REF] Efromovich | A self-training algorithm for nonparametric filtering[END_REF]. The next result, that is a direct consequence of Theorem 1, shows that EWA with linear constituent estimators is also asymptotically sharp adaptive over Sobolev ellipsoids.

Proposition 5 Let λ = (α, w) ∈ Λ = R 2 + and consider the prior

π(d λ) = 2n -α/(2α+1) σ 1 + n -α/(2α+1) σ w 3 e -α d αd w, ( 20
)
where n σ = n/σ 2 . Then, in model ( 1) with homoscedastic errors, the aggregate f EWA based on the temperature β = 8σ 2 and the constituent estimators f α,w = A α,w Y (with A α,w being the Pinsker filter) is adaptive in the exact minimax sense1 on the family of classes {F (α, R)

: α > 0, R > 0}.
It is worth noting that the exact minimax adaptivity property of our estimator f EWA is achieved without any tuning parameter. All previously proposed methods that are provably adaptive in exact minimax sense depend on some parameters such as the lengths of blocks for blockwise Stein and Efromovich-Pinsker estimators or the step of discretization and the maximal value of bandwidth [START_REF] Cavalier | Oracle inequalities for inverse problems[END_REF]. Another nice property of the estimator f EWA is that it does not require any pilot estimator based on the data splitting device [START_REF] Efromovich | On nonparametric regression for IID observations in a general setting[END_REF], [START_REF] Yang | Aggregating regression procedures to improve performance[END_REF].

Experiments

In this section we present some numerical experiments on synthetic data, by focusing only on the case of homoscedastic Gaussian noise (Σ = σ 2 I n×n ) with known variance. Following the philosophy of reproducible research, a toolbox is made available freely for download at: www.math.jussieu.fr/~salmon/code/index_codes.php

We evaluate different estimation routines on several 1D signals, introduced by Donoho and Johnstone (1994,1995) and considered as benchmark in literature on signal processing. The six signals we retained for our experiments because of their diversity are depicted in Figure 1. Since all these signals are non-smooth, we have also carried out experiments on their smoothed versions obtained by taking an antiderivative, see Figure 1. In what follows, the experiment on non-smooth signals will be referred to as Experiment I, whereas the experiment on their smoothed counterparts will be referred to as Experiment II. In both cases, prior to applying estimation routines, we normalize the (true) sampled signal to have an empirical norm equal to one and use the Discrete Sine Transform (DST) denoted by θ(Y ) = θ 1 (Y ), . . . , θ n (Y ) .

The four estimation routines-including EWA-used in our experiments are detailed below:

Soft Thresholding (ST), Donoho and Johnstone (1994): For a given threshold parameter t , the soft thresholding estimator of the vector of DST coefficients θ k ( f ) is defined by

θ k = sgn θ k (Y ) |θ k (Y )| -σt + . ( 21
)
In our experiments, we use the threshold minimizing the estimated unbiased risk defined via Stein's lemma. This procedure is referred to as SURE-shrink in [START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF].

Blockwise James-Stein (BJS) shrinkage, Cai (1999):

The set of indices {1, • • • , n} is partitioned into N = [n/ log(n)] non-overlapping blocks B 1 , B 2 , • • • B N of equal size L. (If n
is not a multiple of N , the last block may be of smaller size than all the others.) The corresponding blocks of true coefficients

θ B k ( f ) = θ j ( f ) j ∈B k
are estimated by shrinking the blocks of noisy coefficients θ B k (Y ):

θ B k = 1 - λLσ 2 S 2 k (Y ) σ + θ B k (Y ), k = 1, • • • , N (22)
where

S 2 k (Y ) = θ B k (Y ) 2
2 and λ = 4.50524 as in [START_REF] Cai | Adaptive wavelet estimation: a block thresholding and oracle inequality approach[END_REF].

Unbiased risk estimate (URE) minimization, Golubev (1992), Cavalier et al. (2002):

it consists in using a Pinsker filter, as defined in Section 4, with a data-driven choice of parameters α and w. This choice is done by minimizing an unbiased estimate of the risk over a suitably chosen grid for the values of α and w. Here, we use geometric grids ranging from 0.1 to 100 for α and from 1 to n for w. The bi-dimensional grid used in all the experiments has 100×100 elements. We refer to [START_REF] Cavalier | Oracle inequalities for inverse problems[END_REF] for the closed-form formula of the unbiased risk estimator.

EWA on Pinsker's filters:

We consider the same finite family of linear smoothers-defined by Pinsker's filters-as in the URE routine described above. According to Proposition 1, this leads to an estimator which is nearly as accurate as the best Pinsker's estimator in the given finite family.

To report the result of our experiments, we have also computed the best linear smoother based on a Pinsker filter chosen among the candidates that we used for defining URE and EWA routines. By best smoother we mean the one minimizing the squared error, which can be computed since we know the ground truth. This pseudo-estimator will be referred to as oracle. The results summarized in Table 1 for Experiment I and Table 2 for Experiment II correspond to the average over 1000 trials of the mean squared error (MSE) from which we subtract the MSE of the oracle and multiply the resulting difference by the sample size. We report the results for σ = 0.33 and for n ∈ {2 8 , 2 9 , 2 10 , 2 11 }.

Simulations show that EWA and URE have very comparable performances and are significantly more accurate than Soft Thresholding and Block James-Stein (see Table 1) for every size n of signals considered. The improvement is particularly important when the signal has large peaks (cf. Figure 2) or discontinuities (cf. Figure 3). In most cases, the EWA method also outperforms URE, but this difference is much less pronounced. One can also observe that in the case of smooth signals, the difference of the MSEs between EWA and the oracle, multiplied by n, remains nearly constant when n varies. This is in perfect agreement with our theoretical results in which the residual term decreases to zero inversely proportionally to the sample size.

Of course, soft thresholding and blockwise James-Stein procedures have been designed for being applied to the wavelet transform of a Besov smooth function, rather than to the Fourier transform of a Sobolev-smooth function. However, the point here is not to demonstrate the superiority of EWA as compared to ST and BJS procedures. The point is to stress the importance of having sharp adaptivity up to optimal constant and not simply adaptivity in the sense of rate of convergence. Indeed, the procedures ST and BJS are provably rate-adaptive when applied to Fourier transform of a Sobolev-smooth function, but they are not sharp adaptive-they do not attain the optimal constant-whereas EWA and URE do attain.

Summary and future work

In this paper, we have addressed the problem of aggregating a set of affine estimators in the context of regression with fixed design and heteroscedastic noise. Under some assumptions on the constituent estimators, we have proven that the EWA with a suitably chosen temperature parameter satisfies PAC-Bayesian type inequality, from which different types of oracle inequalities have been deduced. All these inequalities are with leading constant one and with rate-optimal residual term. As a by-product of our results, we have shown that EWA applied to the family of Pinsker's estimators produces an estimator, which is adaptive in the exact minimax sense. Next in our agenda is carrying out an experimental evaluation of the proposed aggregate using the approximation schemes described by [START_REF] Dalalyan | Sparse regression learning by aggregation and Langevin Monte-Carlo[END_REF], [START_REF] Ph | Exponential screening and optimal rates of sparse estimation[END_REF] and [START_REF] Alquier | Pac-bayesian bounds for sparse regression estimation with exponential weights[END_REF]. It will also be interesting to extend the results of this work to the case of the unknown noise variance in the same vein as in [START_REF] Ch | Mixing least-squares estimators when the variance is unknown[END_REF].

r = Y -f 2 n + 2 n n i =1 σ 2 i ∂ y i fi - 1 n n i =1 σ 2 i , (23)
is an unbiased estimate of r , ie. E r = r .

Proof:

For any i = 1, • • • , n, one has E(Y i -fi ) 2 = E(Y i -f i ) 2 + E( f i -fi ) 2 + 2E (Y i -f i )( f i -fi ) ,
The following identity is the classical Stein Lemma (cf. [START_REF] Tsybakov | Introduction to nonparametric estimation[END_REF] p.157), based on integration by parts:

E (Y i -f i ) fi = σ 2 i E ∂ y i f i . ( 24
)
where the differentiation is according to Y i . Using the last two displays, one has:

E Y -f 2 n = E Y -f 2 n + E f -f 2 n - 2 n E n i =1 σ 2 i ∂ y i f i , (25)
leading to the announced unbiased risk estimate.

B Main Result

Now, we can apply Stein's Lemma for any estimator f λ , so that we can build rλ for any λ ∈ Λ. In this paper, we only focus on affine estimators f λ , i.e., estimators that can be written as affine transforms of the data

Y = (y 1 , • • • , y n ) ∈ R n .
Affine estimators can be defined by

f λ = A λ Y + b λ , (26)
where the n × n real matrix A λ and the vector b λ ∈ R n are deterministic. This means that the entries of A λ and b λ may depend on the design points x 1 , • • • , x n but not on the data vector Y . It is easy to check that the the divergence term n i =1 σ 2 i ∂ y i f i in Stein's Lemma is simply Tr(ΣA λ ) for affine estimators. Then rλ , defined by

rλ = Y -f λ 2 n + 2 n Tr(ΣA λ ) - 1 n n i =1 σ 2 i ( 27
)
is an unbiased estimator of r λ .

In order to state our main result, we denote by P Λ the set of all probability measures on Λ and by K (p, p ) the Kullback-Leibler divergence between two probability measures p, p ∈ P Λ .

K (p, p ) = Λ log d p d p (λ) p(d λ) if p p , +∞ otherwise.
Theorem 1 If either one of the following conditions is satisfied:

C 1 :
The matrices A λ are orthogonal projections (i.e., symmetric and idempotent) and the vectors b λ satisfy A λ b λ = 0, for all λ ∈ Λ.

C 2 : The matrices A λ are all symmetric, positive semidefinite and satisfy

A λ A λ = A λ A λ , A λ Σ = ΣA λ for all λ, λ ∈ Λ. All the vectors b λ are zero.
Then, the aggregate f EWA defined by Equations ( 7), ( 8) and ( 4) satisfies the inequality

E( f EWA -f 2 n ) ≤ inf p∈P Λ Λ E f λ -f 2 n p(d λ) + β n K (p, π) provided that β ≥ α |Σ |, where α = 4 if C 1 holds true and α = 8 if C 2 holds true.
Proof:[when C 2 is satisfied] According to Stein's lemma, the quantity

rEWA = Y -f EWA 2 n + 2 n n i =1 σ 2 i ∂ y i fEWA,i - 1 n n i =1 σ 2 i (28)
is an unbiased estimate of the risk r EWA = E( f EWAf 2 n ). Using simple algebra, one checks that

Y -f EWA 2 n = Λ Y -f λ 2 n -f λ -f EWA 2 n θ(λ)π(d λ). ( 29
)
By interchanging the integral and differential operators, we get the following expression for the derivatives of fEWA,i :

∂ y i fEWA,i = Λ ∂ y i fλ,i θ(λ) π(d λ) + Λ fλ,i ∂ y i θ(λ) π(d λ). ( 30
)
Let us defined A EWA Λ A λ θ(λ)π(d λ). With this notation, the last equality, combined with Equations ( 4), ( 28), ( 29) and the fact that

n i =1 σ 2 i ∂ y i fλ,i = Tr(ΣA λ ), implies that rEWA = Λ rλ -f λ -f EWA 2 n θ(λ) π(d λ) + 2 n n i =1 σ 2 i Λ fλ,i ∂ y i θ(λ) π(d λ).
Taking into account that Λ fEWA,i

∂ y i θ(λ) π(d λ) = fEWA,i ∂ y i Λ θ(λ)π(d λ)
= 0, we come up with the following expression for the unbiased risk estimate:

rEWA = Λ rλ -f λ -f EWA 2 n + 2 ∇ Y log θ(λ)|Σ( f λ -f EWA ) n θ(λ) π(d λ) (31) = Λ rλ -f λ -f EWA 2 n -2nβ -1 ∇ Y rλ |Σ( f λ -f EWA ) n θ(λ) π(d λ). ( 32
)
Note that, so far, the precise form of the constituent estimators has not been exploited. This form is important for computing ∇ Y rλ . In view of Equations ( 26) and ( 4), as well as the assumptions A λ = A λ and b λ ≡ 0 (holding thanks to C 2 ), we get

∇ Y rλ = 2 n (I n×n -A λ ) (I n×n -A λ )Y - 2 n (I n×n -A λ ) b λ = 2 n (I n×n -A λ ) 2 Y . ( 33
)
In what follows, we use the shorthand I = I n×n . Using this notation and Eq. ( 33), we get

rEWA = Λ rλ -f λ -f EWA 2 n - 4 β (I -A λ ) 2 Y |Σ(A λ -A EWA )Y n θ(λ) π(d λ). ( 34
)
Recall now that for any pair of commuting matrices P and Q the identity (I -P ) 2 = (I -Q) 2 +2 I -P +Q 2 (Q-P ) holds true. Applying this formula to P = A λ and Q = A EWA we get the following expression: (I -

A λ ) 2 Y |Σ(A λ -A EWA )Y n = (I -A EWA ) 2 Y |Σ(A λ -A EWA )Y n -2 I - Â+A λ 2 (A EWA -A λ )Y |Σ(A EWA -A λ )Y n .
When one integrates over Λ with respect to the measure θ • π, the term of the first scalar product in the RHS of the last equation vanishes. On the other hand, positive semidefiniteness of matrices A λ implies that of the matrix A EWA and, therefore, I -

Â+A λ 2 (A EWA -A λ )Y |Σ(A EWA -A λ )Y n ≤ 〈(A EWA -A λ )Y |Σ(A EWA - A λ )Y n . This inequality, in conjunction with (34) implies that rEWA ≤ Λ rλ -f λ -f EWA 2 n + 8 β 〈(A EWA -A λ )Y |Σ(A EWA -A λ )Y n θ(λ) π(d λ) = Λ rλ -f λ -f EWA 2 n + 8 β 〈 f EWA -f λ |Σ( f EWA -f λ ) n θ(λ) π(d λ) ≤ Λ rλ -1 - 8 max i σ 2 i β f λ -f EWA 2 n θ(λ) π(d λ).
Taking into account the fact that β ≥ 8 max i σ 2 i , we get rEWA

≤ Λ rλ θ(λ) π(d λ) ≤ Λ rλ π(d λ)+ β n K ( π, π).
To conclude, it suffices to remark that π is the probability measure minimizing the criterion Λ rλ p(d λ) + β n K (p, π) among all p ∈ P Λ (see for instance [START_REF] Catoni | Statistical learning theory and stochastic optimization[END_REF] p.160). Thus, for every p ∈ P Λ , it holds that rEWA

≤ Λ rλ p(d λ) + β n K (p, π).
Taking the expectation of both sides, the desired result follows with Fatou's lemma.

Proof:[ when C 1 is satisfied] We can do the same calculation as when C 2 is satisfied until (32). In view of Equations ( 2) and ( 4), as well as the assumptions A 2 λ = A λ = A λ and A λ b λ ≡ 0, we get

∇ Y rλ = 2 n (I n×n -A λ ) (I n×n -A λ )Y - 2 n (I n×n -A λ ) b λ = 2 n (I n×n -A λ )Y - 2 n b λ . ( 35
)
Using the same shorthand I = I n×n with Eq. ( 35) we come up with

rEWA = Λ rλ -f λ -f EWA 2 n - 4 β Y -f λ |Σ( f λ -f EWA ) n θ(λ) π(d λ). ( 36
)
Now, since f is the expectation of f λ with respect to the measure θ • π, we have

rEWA = Λ rλ -f λ -f EWA 2 n + 4 β 〈Y -f EWA + f EWA -f λ |Σ( f λ -f EWA ) n θ(λ) π(d λ) = Λ rλ -f λ -f EWA 2 n + 4 β 〈 f EWA -f λ |Σ( f EWA -f λ ) n θ(λ) π(d λ) ≤ Λ rλ -1 - 4 max i σ 2 i β f λ -f EWA 2 n θ(λ) π(d λ).
Taking into account the fact that β ≥ 4 max i σ 2 i , we get the same results as with condition

C 2 : rEWA ≤ Λ rλ θ(λ) π(d λ) ≤ Λ rλ π(d λ) + β n K ( π, π).
The end of the proof is unchanged and leads to the same general result as with condition C 2 , except for the choice of α.

C Continuous oracle inequality

Proposition 2 Let Λ ⊂ R M be an open and bounded set and let π be the uniform probability on Λ. Assume that the mapping λ → r λ is Lipschitz continuous, i.e., |r λ -

r λ | ≤ L r λ -λ 2 , ∀λ, λ ∈ Λ. Under the conditions C 1 or C 1 aggregate f EWA satisfies the inequality E f EWA -f 2 n ≤inf λ∈Λ E f λ -f 2 n + βM n log M 2 min(n -1 , d 2 (λ, Λ)) + L r + β log Leb(Λ) n .
for every β ≥ α |Σ | where α = 4 if C 1 holds true and α = 8 if C 2 holds true.

Proof: It suffices to apply Theorem 1 and to bound from above the RHS of inequality ( 9)

E( f EWA -f 2 n ) ≤ inf p∈P Λ Λ r λ p(d λ) + β n K (p, π) E( f EWA -f 2 n ) ≤ inf p∈P Λ Λ |r λ -r λ 0 | + r λ 0 p(d λ) + β n K (p, π)
Then, the RHS of the last inequality can be bounded from above by the minimum over all measures having as density p λ 0 ,τ 0 (λ) = 1l B λ 0 (τ 0 ) (λ)/Leb(B λ 0 (τ 0 )), with λ 0 ∈ Λ and τ 0 = min(1/n, d 2 (λ 0 , Λ)) (hence B λ 0 (τ 0 ) ⊂ Λ). Using the Lipschitz condition on r λ , the bound on the risk becomes

E( f EWA -f 2 n ) ≤ Λ |r λ -r λ 0 | + r λ 0 p λ 0 ,τ 0 (d λ) + β n K (p λ 0 ,τ 0 , π) E( f EWA -f 2 n ) ≤ r λ 0 + L r Λ λ -λ 0 2 p λ 0 ,τ 0 (d λ) + β n K (p λ 0 ,τ 0 , π) E( f EWA -f 2 n ) ≤ r λ 0 + L r τ 0 + β n K (p λ 0 ,τ 0 , π) (37)
Now, since λ 0 is such that B λ 0 (τ 0 ) ⊂ Λ, the measure p λ 0 ,τ 0 (λ) d λ is absolutely continuous w.r.t. π and the Kullback-Leibler divergence between these measures equals log Leb(Λ)/Leb B λ 0 (τ 0 ) . By the simple inequality x 2 2 ≤ M x 2 ∞ for any x ∈ R M , one can see that the Euclidean ball of radius τ 0 contains the hypercube of width 2τ 0 M . So we have the following lower bound for the volume B λ 0 : Leb B λ 0 (τ 0 ) ≥ (2τ 0 / M ) M . By combining this with inequality (37) the results of Proposition 2 is straightforward.

D Sparsity oracle inequality

Let us choose a prior π that promotes the sparsity of λ. This can be done in the same vein as in Dalalyan and Tsybakov Dalalyan andTsybakov (2007, 2008), by means of the heavy tailed prior (Student t (3) distribution):

π(d λ) ∝ M j =1 1 (1 + |λ j /τ| 2 ) 2 1l Λ (λ), (38)
where τ > 0 is a tuning parameter, that takes small values.

Proposition 3 Let Λ = R M and let π be defined by ( 12). Assume that the mapping λ → r λ is continuously differentiable and, for some M × M matrix M , satisfies:

r λ -r λ -∇r λ (λ -λ ) ≤ (λ -λ ) M (λ -λ ), ∀λ, λ ∈ Λ.
If either one of the conditions C 1 and C 2 (cf. Theorem 1) is fulfilled, then the aggregate f EWA defined by Equations ( 7) and ( 4) satisfies the inequality

E f EWA -f 2 n ≤ inf λ∈R M E f λ -f 2 n + 4β n M j =1 log 1 + |λ j | τ + Tr(M )τ 2
provided that β ≥ α max i =1,...,n σ 2 i , where α = 4 if C 1 holds true and α = 8 if C 2 holds true. Proof: The proof is a simplified version of proofs given in Dalalyan andTsybakov (2007, 2008), since Λ is the whole space, Λ = R M instead of a bounded subset of R M .

We begin the proof as for the previous proposition, but pushing the development of the function λ → r λ up to second order. So, for any λ * ∈ R M , we have

E( f EWA -f 2 n ) ≤ inf λ * ∈R M r λ * + Λ ∇r λ * (λ -λ * ) + (λ -λ * ) M (λ -λ * ) p λ * (d λ) + β n K (p λ * , π)
By choosing p λ * (λ) = π(λλ * ) for any λ ∈ R M , the second term in the last display vanishes since the distribution π is symmetric. The third term is computed thanks to the moment of order 2 of a scaled Student t (3) distribution. Recall that if T is drawn from the scaled Student t (3) distribution, its distribution function is u → 2/[π(1 + u 2 ) 2 ], and that ET 2 = 1. Thus, we have that Λ λ 2 1 π(λ)d λ = τ 2 . We can then bound the risk of the EWA estimator by

E( f EWA -f 2 n ) ≤ inf λ * ∈R M r λ * + Tr(M )τ 2 + β n K (p λ * , π) (39)
So far, the particular choice of heavy tailed prior has not been used. This choice is important to control the Kullback-Leibler divergence between two translated versions of the same distribution

K (p λ * , π) = Λ log M j =1 (τ 2 + λ 2 j ) 2 (τ 2 + (λ j -λ * j ) 2 ) 2 p λ * (d λ) K (p λ * , π) = 2 M j =1 Λ log τ 2 + λ 2 j τ 2 + (λ j -λ * j ) 2 p λ * (d λ)
. We bound the quotient in the above equality by

τ 2 + λ 2 j τ 2 + (λ j -λ * j ) 2 = 1 + 2τ(λ j -λ * j ) τ 2 + (λ j -λ * j ) 2 λ * j τ + (λ * j ) 2 τ 2 + (λ j -λ * j ) 2 τ 2 + λ 2 j τ 2 + (λ j -λ * j ) 2 ≤ 1 + λ * j τ + λ * j τ 2 ≤ 1 + λ * j τ 2 .
Since the last inequality is independent of λ, the integral disappears (p λ * is a probability measure) in the previous bound on the Kullback-Leibler divergence, so we eventually get

K (p λ * , π) ≤ 4 M j =1 log 1 + λ * j τ ,
and combine with Inequality (39), this ends the proof of the proposition.

E Application to minimax estimation

Let us denote by θ k ( f ) = 〈 f |ϕ k 〉 n the coefficients of the (orthogonal) discrete sine transform of f and the Sobolev ellipsoids

F (α, R) = { f ∈ R n : n k=1 k 2α θ k ( f ) 2 ≤ R}.
Assume in this section that the noise is homoscedastic (Σ = σ 2 I n×n ) and A α,w = D diag (1k α /w) + ; k = 1, . . . , n D is the Pinsker filter in the discrete sine basis.

Proposition 5 Let λ = (α, w) ∈ Λ = R 2
+ and consider the prior

π(d λ) = 2n -α/(2α+1) σ 1 + n -α/(2α+1) σ w 3 e -α d αd w,
where n σ = n/σ 2 . Then, in model ( 1) with homoscedastic errors, the aggregate f EWA based on the temperature β = 8σ 2 and the constituent estimators f α,w = A α,w Y (with A α,w being the Pinsker filter) is adaptive in the exact minimax sense on the family of classes {F (α, R) : α > 0, R > 0}.

Proof: We assume, without loss of generality, that the matrix n 1/2 D coincides with the identity matrix.

First, let us fix α 0 > 0 and R 0 > 0, such that n -1/2 f ∈ F (α 0 , R 0 ) and define λ 0 = (α 0 , w 0 ) ∈ Λ with w 0 chosen such that the Pinsker estimator f α 0 ,w 0 is minimax over the ellipsoid F (α 0 , R 0 ).

In what follows, we set n σ = n/σ 2 and denote by p π the density of π w.r.t. the Lebesgue measure on R 2 + : p π (α, w) = e -α n -α/(2α+1) σ p w (wn -α/(2α+1) σ ), where p w is a probability density function supported by (0, ∞) such that up w (u) d u = 1. One easily checks that

R 2 α p π (α, w) d αd w = 1, R 2 w p π (α, w) d αd w ≤ n 1/2 σ . ( 40
)
Let τ be a positive number such that τ ≤ min(1, α 0 /(2 log w 0 )) and choose p λ 0 ,τ as a translation/dilatation of π, concentrating on λ 0 when τ → 0:

p λ 0 ,τ (d λ) = p π λ -λ 0 τ d λ τ 2 . In view of Theorem 1, E( f EWA -f 2 n ) ≤ r λ 0 + R 2
r α,wr α 0 ,w 0 p λ 0 ,τ (d λ) + β n K (p λ 0 ,τ , π) . ( 41)

Let us decompose the term r α,wr α 0 ,w 0 into two pieces: r α,wr α 0 ,w 0 = {r α,wr α,w 0 } + {r α,w 0r α 0 ,w 0 } and find upper bounds for the resulting terms. With our choice of estimator, the difference between the risk functions at (α, w) and (α, w 0 ) is:

n(r α,w -r α,w 0 ) = n k=1 (1 -k α /w) + -1 2 -(1 -k α /w 0 ) + -1 2 f 2 k + n k=1 (1 -k α /w) + 2 -(1 -k α /w 0 ) + 2 σ 2
Since the weights of the Pinsker estimators are in [0, 1], we have

n|r α,w -r α,w 0 | ≤2 n k=1 ( f 2 k + σ 2 ) (1 -k α /w) + -(1 -k α /w 0 ) + . ( 42
)
For any x, y ∈ R, the inequality |x +y + | ≤ |x -y| is obvious. Combined with α 0 ≤ α and w 0 ≤ w, we have that

1 - k α w + -1 - k α w 0 + ≤ k α w - k α w 0 1l {k α ≤w} ≤ w -w 0 w 0 . ( 43
)
By virtue of Inequalities ( 42) and ( 43) we get

|r α,w -r α,w 0 | ≤2n -1 n k=1 ( f 2 k + σ 2 ) (w -w 0 ) w 0 ≤ 2(R 0 + σ 2 ) w -w 0 w 0 . ( 44
)
Similar calculations lead to a bound for the other absolute difference between risk functions:

|r α,w 0 -r α 0 ,w 0 | ≤ 2n -1 n k=1 ( f 2 k + σ 2 ) k α -k α 0 w 0 1l {k α 0 ≤w 0 } ≤ 2(R 0 + σ 2 ) w α-α 0 α 0 0 -1 . ( 45
)
Recall that we aim to bound the second term in the RHS of (41). To this end, we need an accurate upper bound on the integrals of the RHSs of ( 44) and (45) w.r.t. the probability measure p λ 0 ,τ . For the first one, we get r α,wr α,w 0 p λ 0 ,τ (d λ) ≤ 2(R 0 + σ 2 )w -1 0 R 2 (ww 0 )p λ 0 ,τ (d λ)

≤ 4n 1/2 σ w -1 0 τ(R 0 + σ 2 ). ( 46)

Similar arguments apply to bound the integral of the second difference between risk functions: where we used the inequality τ ≤ α 0 /(2 log w 0 ). The last term to bound in inequality (41) requires the evaluation of the Kullback-Leibler divergence between p λ 0 ,τ and π. It can be done as follows: where the third equality is derived thanks to Eq. ( 40 Eventually, we can reformulate our bound on the risk of the EWA given in (41), leading to

K (p λ 0 ,τ , π) =
E( f EWA -f 2 n ) ≤ r λ 0 +4τ(R 0 + σ 2 ) n 1/2 σ w 0 + log w 0 α 0 + 8σ 2 (α 0 + 3 log 2+w 0 τ ) n . ( 48
)
To conclude the proof of the proposition, we set

r λ 0 = (1 + o n (1)) min f max f ∈F (α 0 ,R 0 ) E( f -f 2 n ).
In view of this result, taking the max over f ∈ F (α 0 , R 0 ) in (48), we get max f ∈F (α 0 ,R)

E( f EWA -f 2 n ) ≤ (1 + o n (1)) min f max f ∈F (α 0 ,R) E( f -f 2 n ) + O log n n .
This leads to the desired result in view of the relation lim inf

n→∞ min f max f ∈F (α 0 ,R) n 2α 0 2α 0 +1 E( f -f 2 n ) > 0,
which follows from (Tsybakov, 2009, Theorem 3.1).

Figure 1 :

 1 Figure 1: Test signals used in our experiment: Piece-Regular, Ramp, Piece-Polynomial, HeaviSine, Doppler and Blocks. (a) non-smooth (Experiment I) and (b) smooth (Experiment II).

Figure 2 :Figure 3 :

 23 Figure 2: Heavisine. The first row is the true signal (left) and a noisy version corrupted by Gaussian noise with standard deviation σ = 0.33 (right). The second row gives denoised version obtained by EWA (left), BJS, ST and URE (right). The PSNR is computed by the formula PSNR = 10 log 10 max( f ) 2 /MSE .

Figure 4 :

 4 Figure 4: Piece-Regular. The first row is the true signal (left) and a noisy version corrupted by Gaussian noise with standard deviation σ = 1 (right). The second row gives denoised version obtained by EWA (left) and by BJS, ST and URE (right). The PSNR is computed by the formula PSNR = 10 log 10 max( f ) 2 /MSE .

  w 0r α 0 ,w 0 |p λ 0 ,τ (d λ) ≤ 2(R 0 + σ 2 )

  0 ,τ (d λ) -2 log(τ) ≤ α 0 + (τ -1) + τ (d λ) -2 log(τ).

  ) and the obvious relation p w ∞ = 2. Now, making the change of variable w = w 0 + τn α/(2α+1) σ u and using the fact that w 0 + τn α/w 0 + u p w (u)d u≤ 3 log 1 + w 0 + R + up w (u)d u = 3 log(2 + w 0 ).

Table 1 :

 1 Comparison of several adaptive methods on the six (non-smooth) signals of interest. For each signal length n and each method, we give the average value of n×(MSE -MSE Oracle ) and the corresponding standard deviation below, for 1000 replications of the experiment. Negative values indicate that in some cases the EWA procedure has a smaller risk than that of the best linear estimator used for the aggregation, which is possible since the EWA itself is not a linear estimator.

	n	EWA	URE	BJS	ST	EWA	URE	BJS	ST
			Blocks			Doppler
	256	0.051	0.245	9.617	4.846	0.062	0.212 13.233	6.036
		(0.42) (0.39)	(1.78)	(1.29)	(0.35) (0.31)	(2.11)	(1.23)
	512 -0.052	0.302 13.807	9.256	-0.100	0.205 17.080 12.620
		(0.35) (0.50)	(2.16)	(1.70)	(0.30) (0.39)	(2.29)	(1.75)
	1024 -0.050	0.299 19.984 17.569	-0.107	0.270 21.862 23.006
		(0.36) (0.46)	(2.68)	(2.17)	(0.35) (0.41)	(2.92)	(2.35)
	2048 -0.007	0.362 28.948 30.447	-0.150	0.234 28.733 38.671
		(0.42) (0.57)	(3.31)	(2.96)	(0.34) (0.42)	(3.19)	(3.02)
			HeaviSine			Piece-Regular
	256 -0.060	0.247	1.155	3.966	-0.069	0.248	8.883	4.879
		(0.19) (0.42)	(0.57)	(1.12)	(0.32) (0.40)	(1.76)	(1.20)
	512 -0.079	0.215	2.064	5.889	-0.105	0.237 12.147	9.793
		(0.19) (0.39)	(0.86)	(1.36)	(0.30) (0.37)	(2.28)	(1.64)
	1024 -0.059	0.240	3.120	8.685	-0.092	0.291 15.207 16.798
		(0.23) (0.36)	(1.20)	(1.64)	(0.34) (0.46)	(2.18)	(2.13)
	2048 -0.051	0.278	4.858 12.667	-0.059	0.283 21.543 27.387
		(0.25) (0.48)	(1.42)	(2.03)	(0.34) (0.54)	(2.47)	(2.77)
			Ramp			Piece-Polynomial
	256	0.038	0.294	6.933	5.644	0.017	0.203 12.201	3.988
		(0.37) (0.47)	(1.54)	(1.20)	(0.37) (0.37)	(1.81)	(1.19)
	512	0.010	0.293	9.712	9.977	-0.078	0.312 17.765	9.031
		(0.36) (0.51)	(1.76)	(1.67)	(0.35) (0.49)	(2.72)	(1.62)
	1024 -0.002	0.300 13.656 16.790	-0.026	0.321 23.321 17.565
		(0.30) (0.45)	(2.25)	(2.06)	(0.38) (0.48)	(2.96)	(2.28)
	2048	0.007	0.312 19.113 27.315	-0.007	0.314 31.550 29.461
		(0.34) (0.50)	(2.68)	(2.61)	(0.41) (0.49)	(3.05)	(2.95)

Table 2 :

 2 Comparison of several adaptive methods on the six smoothed signals of interest. For each signal length n and each method, we give the average value of n(MSE -MSE Oracle ) and the corresponding standard deviation below, for 1000 replications of the experiment. With the model (1), if the estimator f is almost everywhere differentiable in Y and if each ∂ y i f i has finite first moment, then

	Stein's Lemma 1

see(Tsybakov,

2009, Definition 3.8)

Acknowledgments

The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR) under the grant PARCIMONIE.

Appendix

A Stein's Lemma with heteroscedastic noise

To define the EWA estimator, we first need to determine an unbiased risk estimate for any of the constituent estimators. We adapt a systematic method based on Stein's Lemma to the heteroscedastic framework. We recall this lemma given in [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF], for our setting:


Mikel Rodriguez

Ivan Laptev

Josef Sivic

Jean-Yves Audibert

Density-aware person detection and tracking in crowds

We address the problem of person detection and tracking in crowded video scenes. While the detection of individual objects has been improved significantly over the recent years, crowd scenes remain particularly challenging for the detection and tracking tasks due to heavy occlusions, high person densities and significant variation in people's appearance. To address these challenges, we propose to leverage information on the global structure of the scene and to resolve all detections jointly. In particular, we explore constraints imposed by the crowd density and formulate person detection as the optimization of a joint energy function combining crowd density estimation and the localization of individual people. We demonstrate how the optimization of such an energy function significantly improves person detection and tracking in crowds. We validate our approach on a challenging video dataset of crowded scenes.

Introduction

Detecting and tracking people in crowded scenes is a crucial component for a wide range of applications including surveillance, group behavior modeling and crowd disaster prevention. The reliable person detection and tracking in crowds, however, is a highly challenging task due to heavy occlusions, view variations and varying density of people as well as the ambiguous appearance of body parts, e.g. the head of one person could be similar to a shoulder of a near-by person. High-density crowds, such as illustrated in Figure 1, present particular challenges due to the difficulty of isolating individual people with standard low-level methods of background subtraction and motion segmentation typically applied in low-density surveillance scenes.

In recent years significant progress has been made in the field of object detection and recognition [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF][START_REF] Everingham | The PASCAL Visual Object Classes Challenge[END_REF][START_REF] Felzenszwalb | Object detection with discriminatively trained part based models[END_REF]. While standard "scanning-window" methods attempt to localize objects independently, several recent approaches extend this work and exploit scene context as well as relations among objects for improved object detection [START_REF] Desai | Discriminative models for multi-class object layout[END_REF][START_REF] Rabinovich | Objects in context[END_REF][START_REF] Torralba | Contextual priming for object detection[END_REF][START_REF] Yao | Modeling mutual context of object and human pose in human-object interaction activities[END_REF]. Related ideas have been investigated for human motion analysis where incorporating scene-level and behavioral factors effecting the spatial arrangement and movement of people have been shown beneficial for achieving improved detection and tracking accuracy. Examples of explored cues include: the destination of a pedestrian within the scene [START_REF] Pellegrini | You'll never walk alone: Modeling social behavior for multi-target tracking[END_REF], repulsion from near-by agents due to the preservation of personal space and social grouping behavior [START_REF] Breitenstein | Robust tracking-by-detection using a detector confidence particle filter[END_REF], as well as the speed of an agent in the group [START_REF] Johansson | Specification of the social force pedestrian model by evolutionary adjustment to video tracking data[END_REF].

We follow this line of work and extend it to the detection and tracking of people in high-density crowds. Rather than modeling individual interactions of people, this work exploits information at a more global level provided by the crowd density and scene geometry. Crowd density estimation has been addressed in a number of recent works which often pose it as a regression problem [START_REF] Chan | Privacy preserving crowd monitoring: Counting people without people models or tracking[END_REF][START_REF] Kong | A viewpoint invariant approach for crowd counting[END_REF][START_REF] Lempitsky | Learning to count objects in images[END_REF]. Such methods avoid the hard detection task and attempt to infer person counts directly from low-level image measurements, e.g. histograms of feature responses. Such methods, hence, provide person counts in image regions but are uncertain 1 about the location of people in these regions. This information is complementary to the output of standard person detectors which optimize precise localization of individual people but lack the global knowledge on the crowd structure. Our goal and contribution is to combine these two sources of complementary information for improved person detection and tracking. The intuition behind our method is illustrated in Figure 1 where the constraints of person counts in local image regions help improving the standard head detector.

We formulate our method in an energy minimization framework which combines crowd density estimates with the strength of individual person detections. We minimize this energy by jointly optimizing the density and the location of individual people in the crowd. We demonstrate how such optimization leads to significant improvements of state-of-the-art person detection in crowded scenes with varying densities. In addition to crowd density cues, we explore constraints provided by scene geometry and temporal continuity of person tracks in the video and demonstrate further improvements for person tracking in crowds. We validate our approach on challenging crowded scenes from multiple video datasets.

The rest of the paper is organized as follows. We put our method in the context of related work in Section 2. Section 3 describes our model, its individual components as well as the optimization procedure for our objective function. Finally, the proposed model is experimentally evaluated in Section 4.

Related work

Object recognition and localization has been studied extensively in computer vision and substantial progress has been achieved over the last decade [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF][START_REF] Everingham | The PASCAL Visual Object Classes Challenge[END_REF][START_REF] Felzenszwalb | Object detection with discriminatively trained part based models[END_REF]. Recent studies demonstrate that the explicit modeling of scene-level information can aid performance of individual object-level detectors. For example, [START_REF] Torralba | Contextual priming for object detection[END_REF] uses scene priming to constrain locations of objects in the image, [START_REF] Barinova | On detection of multiple object instances using hough transforms[END_REF][START_REF] Wu | Detection of multiple, partially occluded humans in a single image by bayesian combination of edgelet part detectors[END_REF] seek detections that explain a set of detected parts, whereas [START_REF] Desai | Discriminative models for multi-class object layout[END_REF][START_REF] Rabinovich | Objects in context[END_REF] model spatial and co-occurence interactions between objects in the scene while [START_REF] Yao | Modeling mutual context of object and human pose in human-object interaction activities[END_REF] attempts to explore relations between human poses and manipulated objects. Our detection method relates to this previous work and extends it by exploring the new scene-level constraint provided by the crowd density.

In the area of tracking, [START_REF] Breitenstein | Robust tracking-by-detection using a detector confidence particle filter[END_REF][START_REF] Khan | MCMC-based particle filtering for tracking a variable number of interacting targets[END_REF][START_REF] Nillius | Multi-target tracking: Linking identities using bayesian network inference[END_REF][START_REF] Pellegrini | You'll never walk alone: Modeling social behavior for multi-target tracking[END_REF] demonstrate how the tracking accuracy can be improved by the joint modeling of motion of multiple people in the scene. For example, [START_REF] Nillius | Multi-target tracking: Linking identities using bayesian network inference[END_REF] tracks football players and explores the roles of players to resolve tracking ambiguity. Others study tracking in crowded scenes containing, for example, groups of ants [START_REF] Khan | MCMC-based particle filtering for tracking a variable number of interacting targets[END_REF] or biological cells [START_REF] Li | Cell population tracking and lineage construction using multiple-model dynamics filters and spatiotemporal optimization[END_REF]. Given the difficulty of tracking individuals in crowds, several approaches leverage common motion patterns as constraints for the track-ing problem. In [START_REF] Ali | Floor Fields for Tracking in High Density Crowd Scenes[END_REF] global motion patterns are learned and participants of the crowd are assumed to behave in a similar manner. Overlapping motion patterns have been studied in [START_REF] Rodriguez | Tracking in unstructured crowded scenes[END_REF] as a means of coping with multi-modal crowd behaviors and [START_REF] Rodriguez | Data-driven crowd analysis in videos[END_REF] uses a large database of crowds as a datadriven prior for tracking. Our method is complementary to this prior work and can benefit from its sophisticated motion models. Differently to the prior work on high-density crowd tracking, our method builds on the tracking-by-detection framework in [START_REF] Everingham | Hello! My name is... buffy-automatic naming of characters in tv video[END_REF] and does not require manual initialization.

The problem of determining the density of objects in a scene has been also studied extensively. Several density estimation methods are based on aggregating counts obtained from local object detectors. These methods typically assume that image regions of target objects have distinct motion and/or appearance [START_REF] Descombes | Object extraction using a stochastic birth-and-death dynamics in continuum[END_REF][START_REF] Rabaud | Counting crowded moving objects[END_REF]. While this assumption may be applicable to low-and medium-density crowds, it usually does not hold for high-density scenes. An alternative commonly used density estimation approach is based on regression. These methods [START_REF] Chan | Privacy preserving crowd monitoring: Counting people without people models or tracking[END_REF][START_REF] Cho | A neural-based crowd estimation by hybrid global learning algorithm[END_REF][START_REF] Kong | A viewpoint invariant approach for crowd counting[END_REF][START_REF] Lempitsky | Learning to count objects in images[END_REF][START_REF] Marana | Estimation of crowd density using image processing[END_REF] typically forego the challenges of detecting individual agents and instead focus on learning a mapping between density and a set of feature responses. In our work we apply the regression-style crowd density estimation approach and use it to improve localization of individual people.

Crowd Model

This section presents our method for detecting and tracking people in crowds. Section 3.1 formulates the new model for density-informed person detection. The components of this model, i.e. person detection and crowd density estimation are described in Sections 3.2 and 3.3 respectively. Section 3.4 describes our tracking framework.

Energy formulation

We formulate the density-informed person detection as follows. We assume to have a confidence score s(p) of a person detector for each location p i , i = 1...N in an image. In addition, we assume we are given a person density, i.e. the number of people per pixel, D(p i ) estimated in a window of size σ at each location p i . Our goal is to identify locations of people in the image such that the sum of detector confidence scores at those locations is maximized while respecting the density of people given by D and preventing significantly overlapping detections, i.e. detections with the area overlap greater than a certain threshold. Using similar notation as in [START_REF] Desai | Discriminative models for multi-class object layout[END_REF], we encode detections in the entire image by a single N -vector x ∈ {0, 1} N , where x i = 1 if the detection at p i is "switched on" and 0 otherwise. The detection problem can be then formulated as a minimization of the following cost function

min x∈{0,1} N -s x E S + x W x E P + α D -Ax 2 2 E D . ( 1
)
Minimizing the first term, E S , in (1) ensures the high confidence values of the person detector at locations of detected people (indicated by x i = 1). The second, pairwise, term E P ensures that only valid configurations of non-overlapping detections are selected. This is achieved by setting W ij = ∞, if detections at locations p i and p j have significant area overlap ratio, and 0 otherwise. The first two terms of the cost function are similar in spirit to the formulation used in [START_REF] Desai | Discriminative models for multi-class object layout[END_REF] and implement a variation of the standard non-maximum suppression. In addition, we introduce a new term, E D , that models the crowd density constraints by penalizing the differences between the density values (i) measured with a regression-based density estimator D and (ii) obtained by counting the "switched on" (or active) detections in x. The evaluation of the density of active detections in x is performed by matrix multiplication Ax, where A is a N × N matrix with rows

A i A i (q j ) = 1 2πσ 2 exp - p i -q j 2 2σ 2 (2)
corresponding to Gaussian windows of size σ centered at positions p i . The idea of minimizing the term E D is illustrated in Figure 2. Intuitively, optimizing the cost (1) including the third, density, term E D enables improving person detection by penalizing confident detections in low person density image regions while promoting low-confidence detections in high person density regions.

In the simplest form, when σ is equal to the width of the entire image, the density term E D reduces (up to a multiplicative constant) to penalizing the difference between the total number of switched on detections in x and the total number of people in the image obtained by the density estimator D.

Note that results of the optimization will be naturally dependent on the quality of the used detector s and density estimator D. Moreover, the quality of D is expected to decrease for small window sizes σ, since if the density estimates D were reliable for arbitrary small σ, having D would solve the detection problem. To balance the contributions of person detection and density estimation, we introduce in (1) a weighting parameter α, which we set manually during training. We also select appropriate window size σ resulting in reasonable density estimates as will be explained in Section 3.3.

Optimization. Optimizing (1) with respect to a discrete vector x is a NP-hard problem in general. As the cost function (1) has a quadratic form, one could attempt optimizing it using quadratic programming by relaxing the binary values of x to x ∈ [0, 1] N . As the quadratic form is not convex (because of the E P term), there is no real gain in the computational complexity to solve the relaxed problem. We thus follow a different approach and optimize (1) using a greedy search procedure similar to [START_REF] Desai | Discriminative models for multi-class object layout[END_REF]. We initialize x = 0, then at each new iteration we update x with x i = 1 by activating a detection at p i for i which decreases the cost function by the largest amount. We continue iterations until the cost in (1) cannot be decreased further by activating more detections. Similar to [START_REF] Desai | Discriminative models for multi-class object layout[END_REF], we found this minimization strategy to work well in practice.

Person detection

In this section we describe the multi-scale person detector, which provides the dense detection score map s(p) used in the cost function (1). Person detection is a widely studied problem with many methods and evaluation benchmarks available [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF][START_REF] Everingham | The PASCAL Visual Object Classes Challenge[END_REF][START_REF] Felzenszwalb | Object detection with discriminatively trained part based models[END_REF]. Most of the methods consider fullbody (pedestrians) or upper-body detection as these views are typical for streetview scenes, movies and consumer photographs. In dense crowd scenes, however, the body of the person does not provide reliable cues due to frequent occlusions by other people. In this work we use the stateof-the-art object detector [START_REF] Felzenszwalb | Object detection with discriminatively trained part based models[END_REF] and train it on extended head regions obtained by manual annotation of head rectangles in our training images of crowds. To improve discrimination, we sample negatives from the same images by allowing a small overlap with positive samples. Once trained, the detector is applied to all positions and scales in the image to provide a dense score map s(p). To localize individual people a non-maximum suppression procedure is usually applied to s(p) [START_REF] Felzenszwalb | Object detection with discriminatively trained part based models[END_REF]. In our work we use s(p) in combination with density estimation D(p) in (1). Results of person detections in crowds are illustrated in Figures 1 and4, and are quantitatively evaluated in Section 4.

Geometric consistency. Footage of crowd scenes typically involves people moving on a ground plane filmed with an overhead camera (see Figure 3). This imposes strong prior on the plausible positions and sizes of people in the frame.

We use this prior to improve results of the basic person detector described above. For this purpose, we seek to roughly estimate the camera viewpoint, which in turn will give us constraints on person sizes observed across the image. We follow the approximation of scene-camera geometry developed in [START_REF] Hoiem | Putting Objects in Perspective[END_REF] and [START_REF] Gallagher | Understanding images of groups of people[END_REF] for generic object detection and adapt it to constrain head detection in crowd videos. We assume people stand on the ground plane, and further assume that (bottom tip of) their heads lie in a single "head plane", parallel with the ground plane. The camera position with respect to the scene is parameterized by two variables: the height of the camera, y c , and the camera tilt, parameterized by the position of the horizon in the image v 0 . We assume there is no horizontal or in-plane rotation of the camera. Given this simplified model, the 3D head size y i of a person i is related to its measured image height h i at vertical image co-ordinate v i as

y i = h i y c v i -v 0 . ( 3
)
Assuming known 3D head size, equal to the average person head (22.6cm [START_REF] Marieb | Human anatomy & physiology[END_REF]), we can estimate the two unknown camera parameters (y c , v 0 ) from positions and heights of two detections in the image. In our experiments, we estimate the camera parameters using least squares fit from the top m head detection hypothesis with the highest confidences (see figure 4). While the scene-camera geometry model is only approximate, we show that the resulting geometric constraints can significantly improve the performance of the detector in Section 4.

Having obtained an estimate for the camera parameters, we can in turn use them to prune head detection hypotheses which are not consistent with the geometry of the scene. In particular, given the combined detection score map obtained over a set of geometrically plausible scales, a detection hypothesis for a given location in the frame is compared the reported height of the head detection in real-world coordinates to that of the average person head. When the absolute difference between the head heights of the detection hypothesis and the average head size is too large, the detection hypothesis is suppressed.

Crowd density estimation

Here we outline the density estimator D(p) used to drive the density-informed detection given by cost function (1). Crowd density estimation has been addressed in a number of works [START_REF] Chan | Privacy preserving crowd monitoring: Counting people without people models or tracking[END_REF][START_REF] Cho | A neural-based crowd estimation by hybrid global learning algorithm[END_REF][START_REF] Lempitsky | Learning to count objects in images[END_REF][START_REF] Marana | Estimation of crowd density using image processing[END_REF][START_REF] Rabaud | Counting crowded moving objects[END_REF]. While some of this work is based on counting detected people, here we seek an alternative approach to complement our person detector. We therefore resort to regression-based density estimation and follow the recent method of [START_REF] Lempitsky | Learning to count objects in images[END_REF]. We assume a set of training images with dense feature map φ(p) ∈ R m at each pixel p and ground-truth annotations of head positions ξ i . We define ground truth density F 0 (p) as a kernel density estimate based on positions of annotated points as

F 0 (p) = 1 2πτ 2 i exp - ξ i -p 2 2τ 2 ,
where we use a 2D Gaussian kernel with a small bandwidth τ = 4 corresponding to the size of heads in our feature map φ. Following [START_REF] Lempitsky | Learning to count objects in images[END_REF], we learn a linear transformation of feature responses that approximates the density function at each pixel F (p|w) = w φ(p) where w ∈ R m is a weight vector. We learn w from k training samples by minimizing the regularized MESA distance D introduced in [START_REF] Lempitsky | Learning to count objects in images[END_REF] as

ŵ = argmin w∈R m w w + λ k D(F 0 k (•), F k (•|w)) . ( 4
)
The density estimator above can be used with different types of feature maps φ. In this paper we define φ using the dense scores s(p) of the person detector described in Section 3.2. In particular, we define the feature map φ i (p), i = 1...m by smoothing detector scores s(p) that are obtained from a range of geometrically plausible scales with m Gaussian filters G(•; t i ) of different size t i .

While F (p) provides per-pixel estimate of the density function, here we aggregate it in Gaussian windows of size σ and define D(p) in (1) as D(p) = F (p) * G(•; σ). As mentioned earlier, σ is an important parameter related to the quality of the density estimator, small values of σ result in stronger density prior in our framework, but also have less reliable results of density estimation. To balance this trade-off, we choose σ by cross-validation and fix its value to 90 pixels in our experiments. We note that adapting the value of σ to the estimated size of people in the image may provide further improvements.

Tracking detections

In this section we describe how the person detections, obtained in individual frames by minimizing the densityaware cost function (1), are associated into person-tracks over multiple video frames. The objective here is to to associate head detections in individual frames into a set of head tracks corresponding to the same person within the crowd across time, thereby generating tracks of individuals that conform to the crowd. We follow the tracking by detection approach of [START_REF] Everingham | Hello! My name is... buffy-automatic naming of characters in tv video[END_REF], which demonstrated excellent performance in tracking faces in TV footage, but here apply it to track heads in crowded video scenes. The method uses local point tracks throughout the video to associate detections of the same person obtained in individual frames. For each crowd video sequence, we obtain point tracks using the Kanade-Lucas-Tomasi tracker [START_REF] Shi | Good features to track[END_REF]. The point tracks are used to establish correspondence between pairs of heads that have been detected within the crowd.

For a given pair of head detections in different frames, the number of point tracks which pass through both heads is counted, and if this number is large relative to the number of point tracks which are not in common to both heads, a match is declared. We use a single-link agglomerative grouping strategy as in [START_REF] Everingham | Hello! My name is... buffy-automatic naming of characters in tv video[END_REF] which gradually merges head detections into larger groups starting from the closest (most connected) detections.

This simple tracking procedure is extremely robust and can establish matches between head detections where the head has not been continuously detected due to pose variation or partial occlusions due to other members of the crowd. In Section 4 we demonstrate the improvement in detection performance using this type of tracking by association: missing detection below detection threshold can be filled-in, and short tracks corresponding to false positive detections can be discarded.

Evaluation

This section evaluates results of our method on the task of detecting and tracking people in crowded video scenes. For comparison, the plot also shows performance of the density-aware detector using the ground truth density, obtained by smoothing ground truth detections by a Gaussian with different sigmas. Note the improvement in performance for smaller sigmas. In this case, for sigma approaching zero, the density would approach the ground truth and hence the perfect performance.

Dataset. For the training and testing we have collected videos of crowd scenes from the Internet and recorded our own set of videos with scenes of political rally. These videos contain considerable variation in terms of viewing angle, scale, crowd motion and crowd density. The camera remained static in all the videos. For the training and evaluation of detection accuracy, we have annotated heads of all visible people in randomly selected frames. All the videos were divided into three subsets. The first subset with 1200 annotated head bounding boxes was used to train a person detector [START_REF] Felzenszwalb | Object detection with discriminatively trained part based models[END_REF] (Section 3.2). The second subset with all people annotated in 60 frames was used to train a density estimator [START_REF] Lempitsky | Learning to count objects in images[END_REF] (Section 3.3). The last test subset with 1009 annotated head bounding boxes was deployed to evaluate the performance of the detection.

In order to evaluate the task of tracking, we employed a subset of 13 video clips captured at a large political rally. Examples of the video frames from this tracking testing set are depicted in Figures 1 and8. On average, each video clip is roughly two minutes long with a frame size of 720 × 480. As can be seen, this collection of crowd video clips spans a wide range of densities, viewpoints and zoom levels. We have annotated 122 person tracks to evaluate the tracking performance on videos from the test subset. This data is available on the project website [1].

Detection. To test and compare the detection performance, we follow PASCAL VOC evaluation protocol [START_REF] Everingham | The PASCAL Visual Object Classes Challenge[END_REF]. A predicted bounding box is considered correct if it overlaps more than 50% with a ground-truth bounding box, only Figure 5. Tracking RMS error (in pixels) for different tracks with and without using the proposed density-aware detector. The error is measured with respect to manually labelled tracks. Note that in many cases the error is significantly reduced using the proposed method.

one detection overlapping with ground-truth bounding box is considered correct, other overlapping detections are declared as false positives. The Performance is measured in terms of precision-recall and average precision (AP) values. To demonstrate the advantage of our method on the detection task, we have compared it to three alternative detectors. Our first baseline detector is [START_REF] Felzenszwalb | Object detection with discriminatively trained part based models[END_REF] trained on our training data. The second detector integrates the baseline detector with geometric filtering imposing the prior on the size of detections as described in Section 3.2. The third detector further integrates temporal consistency constraints using agglomerative clustering described in Section 3.4. Finally, our density-aware detector optimizes the introduced cost function (1) and integrates geometric filtering and temporal consistency constraints as in the case of other detectors. The comparative evaluation is shown in Figure 6. As can be observed, the density-aware detector (red curve) outperforms all three other detectors by a large margin. Qualitative detection results are illustrated in Figure 8.

To gain understanding of the density prior introduced in this work, Figure 6 also shows detection results for the density-aware detector using ground truth density estimates obtained by smoothing ground truth person detections with Gaussians of varying σ (green curves). Interestingly, the detection performance increases significantly in this case, suggesting that our detector can further benefit from future better performing density estimation methods. As expected, the performance of the detector increases for better localized ground truth density estimators with smaller values of σ.

We finally illustrate detection results over scene regions with different person densities in Figure 7. We observe that the density-aware detector (red bars) obtains highest gain for the medium range of density values (2-7 people per 4m 2 ). This can be explained by the increasing error of the density estimator for low-density scenes, and decreasing performance of the person detector for the very high density scenes.

Tracking. The objective of this set of experiments is to assess the improvement that can be attained in tracking accuracy using the proposed density-aware crowd model in the presence of a range of crowd densities.

Quantitative analysis of the proposed tracking algorithm was performed by generating ground-truth trajectories for 122 people, which were selected randomly from the set of all people in the crowd. The ground-truth was generated by manually tracking the centroid of each selected person across the video. In our experiments we evaluate tracks independently, by measuring tracking error (measured in pixels) which is achieved by comparing the tracker position at each frame with respect to the position indicated by the ground truth. When our system does not detect a person that has been labeled in the ground truth, this corresponding track is considered lost. In total, our system was able to detect and track 89 out of the 122 labelled individuals.

A set of trajectories generated by our tracking algorithm is shown in Figure 8. The average tracking error for the annotated heads that were detected obtained using the proposed model was 52.61 pixels (Figure 5). In order to assess the contribution of density estimation in tracking accuracy, a baseline tracking procedure consisting of detection, geometric filtering and agglomerative clustering was evaluated. The mean tracking error of this baseline algorithm was 64.64 pixels. This increase in tracking accuracy can be attributed to the fact the association of individual head detections into tracks contains fewer false positives.

Detection results

Tracking results

We further evaluated the ability to track people over a span of frames by measuring the difference in the length of the generated tracks in relation to the manually annotated tracks. The mean absolute difference between the length of the ground truth tracks and the tracks generated by our system was 18.31 frames, whereas the baseline (which does not incorporate density information) resulted in a mean difference of 30.48 frames. It can be observed from these results that our tracking was very accurate, in most cases, and able to maintain correct track labels over time.

Conclusion

We have shown that automatically obtained person density estimates can be used to improve person localization and tracking performance in crowded scenes. We have formulated the person detection task as a minimization of a joint energy function incorporating scores of individual detections, pair-wise non-overlap constraints, and constraints imposed by the estimated person density over the scene. We have demonstrated significant gains in detection and tracking performance on challenging videos of crowded scenes with varying density. Currently, video frames are processed individually and obtained detections are tracked in postprocessing. Next, we plan to extend the proposed cost function to spatio-temporal video volume to jointly detect and track people in multiple frames.

Figure 1 .

 1 Figure 1. Individual head detections provided by state-of-the-art object detector [12] (bottom-left; green: true positives; red: false positives) are improved significantly by our method (bottom-right; yellow: new true positives) using the crowd density estimate (topright) obtained from the original frame (top-left).

Figure 2 .

 2 Figure 2. Illustration of the energy term ED of (1). Minimizing ED implies reducing the difference in person density estimates obtained by the estimator D(p) (blue color coding) and by local counting person detections (red color coding).

Figure 3 .

 3 Figure 3. Estimating approximate camera-scene geometry from head detections. Assuming known average head size in the scene, the camera height and tilt (parametrized by the location of the horizon line) can be computed from imaged positions and scales of two head detections in the image. See text for more details.

Figure 4 .

 4 Figure 4. (a) The top most confident head detections (shown in red) used to estimate the camera parameters (height and tilt). The camera parameters are in turn used to constrain the heights of detected heads in the scene. (b) Visualization of the metric rectification of a frame of a crowd video using the estimated scene geometry.

Figure 6 .

 6 Figure 6. Evaluation of person detection performance. Precisionrecall curves for the (a) baseline detector, (b) after geometric filtering, (c) tracking by agglomerative clustering and (d) using the proposed density-aware person detector. Note the significant improvement in detection performance obtained by the density-aware detector.For comparison, the plot also shows performance of the density-aware detector using the ground truth density, obtained by smoothing ground truth detections by a Gaussian with different sigmas. Note the improvement in performance for smaller sigmas. In this case, for sigma approaching zero, the density would approach the ground truth and hence the perfect performance.

Figure 7 .

 7 Figure 7. Person detection performance, measured by average precision (y-axis) for different person densities (x-axis). Note that proposed density-aware detector outperforms the baseline detector across almost all density-levels, and specially medium-level densities of 2-7 people in 2x2 meters.

Figure 8 .

 8 Figure 8. Examples of detection and tracking results for different crowded scenes and levels of person density. For detection (top), true positives are shown in green, false positives are shown in red. The additional true positives detected by the proposed density-aware detector (compared to the baseline with geometric filtering and tracking) are shown in yellow. False positives removed by the proposed method are shown by red dashed line. For tracking results (bottom), head detections are depicted with green bounding boxes along with their associated tracks over 100 frames. Additional results and videos are available on the project web page [1].

WILLOW project, Laboratoire d'Informatique de l' École Normale Supérieure, ENS/INRIA/CNRS UMR 8548.

Acknowledgements. This work was partly supported by the Quaero, OSEO, MSR-INRIA, ANR DETECT (ANR-09-JCJC-0027-01) and the CROWDCHECKER project. We thank Pierre Bernas, Philippe Drabczuk, and Guillaume Ne from Evitech for the helpful discussions and the testing videos.


Mathilde Munier

Fredéric Law

Géri Meduri

Damien Le Menuet

Marc Lombes

Mineralocorticoid receptor

teaching and research institutions in France or abroad, or from public or private research centers.

  Figure1