Pac-bayesian bounds for sparse regression estimation with exponential weights

We consider the sparse regression model where the number of parameters $p$ is larger than the sample size $n$. The difficulty when considering high-dimensional problems is to propose estimators achieving a good compromise between statistical and computational performances. The BIC estimator for instance performs well from the statistical point of view \cite{BTW07} but can only be computed for values of $p$ of at most a few tens. The Lasso estimator is solution of a convex minimization problem, hence computable for large value of $p$. However stringent conditions on the design are required to establish fast rates of convergence for this estimator. Dalalyan and Tsybakov \cite{arnak} propose a method achieving a good compromise between the statistical and computational aspects of the problem. Their estimator can be computed for reasonably large $p$ and satisfies nice statistical properties under weak assumptions on the design. However, \cite{arnak} proposes sparsity oracle inequalities in expectation for the empirical excess risk only. In this paper, we propose an aggregation procedure similar to that of \cite{arnak} but with improved statistical performances. Our main theoretical result is a sparsity oracle inequality in probability for the true excess risk for a version of exponential weight estimator. We also propose a MCMC method to compute our estimator for reasonably large values of $p$.


Introduction
We observe n independent pairs (X 1 , Y 1 ), . . . , (X n , Y n ) ∈ X × R (where X is any measurable set) such that where f : X → R is the unknown regression function and the noise variables W 1 , . . . , W n are independent of the design (X 1 , . . . , X n ), satisfy also EW i = 0 and EW 2 i σ 2 for any 1 i n and for some known σ 2 > 0. The distribution of the sample is denoted by P, the corresponding expectation is denoted by E. For any function g : X → R define g n = n i=1 g(X i ) 2 /n 1/2 and g = E g 2 n 1/2 . Let F = {φ 1 , . . . , φ p } be a set-called dictionary-of functions φ j : X → R such that φ j = 1 for any j (this assumption can be relaxed). For any θ ∈ R p define f θ = p j=1 θ j φ j , the empirical risk , and the integrated risk where {(X ′ 1 , Y ′ 1 ), . . . , (X ′ n , Y ′ n )} is an independent replication of {(X 1 , Y 1 ), . . . , (X n , Y n )}. Let us choose θ ∈ arg min θ∈R p R(θ). Note that the minimum may not be unique, however we do not need to treat the identifiability question since we consider in this paper the prediction problem, i.e., find an estimatorθ n such that R(θ n ) is close to min θ∈R p R(θ) up to a positive remainder term as small as possible.
It is a known fact that the least-square estimatorθ LSE n ∈ arg min θ∈R p r(θ) performs poorly in high-dimension p > n. Indeed, consider for instance the deterministic design case with i.i.d. noise variables N (0, σ 2 ) and a full-rank design matrix, thenθ LSE satisfies In the same context, assume now there exists a vector θ ∈ arg min θ∈R p R(θ) with a number of nonzero coordinates p 0 ≤ n. If the indices of these coordinates are known, then we can construct an estimatorθ 0 n such that The estimatorθ 0 n is called oracle estimator since the set of indices of the nonzero coordinates of θ is unknown in practice. The issue is now to build an estimator, when the set of nonzero coordinates of θ is unknown, with statistical performances close to that of the oracle estimatorθ 0 n . A possible approach is to consider solutions of penalized empirical risk minimization problems: where the penalization pen(θ) is proportional to the number of nonzero components of θ such as for instance AIC, C p and BIC criteria [1,31,37]. Bunea, Tsybakov and Wegkamp [9] established for the BIC estimatorθ BIC n the following non-asymptotic sparsity oracle inequality. For any ǫ > 0 there exists a constant C(ǫ) > 0 such that for any p 2, n 1 we have Despite good statistical properties, these estimators can only be computed in practice for p of the order at most a few tens since they are solutions of nonconvex combinatorial optimization problems. Considering convex penalty function leads to computationally feasible optimization problems. A popular example of convex optimization problem is the Lasso estimator (cf. Frank and Friedman [22], Tibshirani [39], and the parallel work of Chen, Donoho and Saunders [15] on basis pursuit) with the penalty term pen(θ) = λ|θ| 1 , where λ > 0 is some regularization parameter and, for any integer d ≥ 2, real q > 0 and vector z ∈ R d we define |z| q = ( d j=1 |z q j |) 1/q and |z| ∞ = max 1≤j≤d |z j |. Several algorithms allow to compute the Lasso for very large p, one of the most popular is known as LARS, introduced by Efron, Hastie, Johnstone and Tibshirani [21]. However, the Lasso estimator requires strong assumptions on the matrix A = (φ j (X i )) 1 i n,1 j p to establish fast rates of convergence results. Bunea, Tsybakov and Wegkamp [8] and Lounici [30] assume a mutual coherence condition on the dictionary. Bickel, Ritov and Tsybakov [7] and Koltchinskii [25] established sparsity oracle inequalities for the Lasso under a restricted eigenvalue condition. Candès and Tao [11] and Koltchinskii [26] studied the Dantzig Selector which is related to the Lasso estimator and suffers from the same restrictions. See, e.g., Bickel, Ritov and Tsybakov [7] for more details. Several alternative penalties were recently considered. Zou [45] proposed the adaptive Lasso which is the solution of a penalized empirical risk minimization problem with the penalty pen(θ) = λ p j=1 1 |ŵj | |θ j | whereŵ is an initial estimator of θ. Zou and Hastie [46] proposed the elastic net with the penalty pen(θ) = λ 1 |θ| 1 + λ 2 |θ| 2 2 , λ 1 , λ 2 > 0. Meinshausen and Bühlmann [35] and Bach [6] considered bootstrapped Lasso. See also Ghosh [23] or Cai, Xu and Zhang [10] for more alternatives to the Lasso. All these methods were motivated by their superior performances over the Lasso either from the theoretical or the practical point of view. However, strong assumptions on the design are still required to establish the statistical properties of these methods (when such results exist). A recent paper by van de Geer and Bühlmann [42] provides a complete survey and comparison of all these assumptions.
Simultaneously, the PAC-Bayesian approach for regression estimation was developed by Audibert [4,5] and Alquier [2,3], based on previous works in the classification context by Catoni [12][13][14], Mc Allester [34], Shawe-Taylor and Williamson [38], see also Zhang [44] in the context of density estimation. This framework is very well adapted for studying the excess risk R(·) − R(θ) in the regression context since it requires very weak conditions on the dictionary. However, the methods of these papers are not designed to cover the high-dimensional setting under the sparsity assumption. Dalalyan and Tsybakov [16][17][18][19] propose an exponential weights procedure related to the PAC-Bayesian approach with good statistical and computational performances. However they consider deterministic design, establishing their statistical result only for the empirical excess risk instead of the true excess risk R(·) − R(θ).
In this paper, we propose to study two exponential weights estimation procedures. The first one is an exponential weights combination of the least squares estimators in all the possible sub-models. This estimator was initially proposed by Leung and Barron [28] in the deterministic design setting. Note that in the literature on aggregation and exponential weights, the elements of the dictionary are often arbitrary preliminary estimators computed from a frozen fraction of the initial sample so that these estimators are considered as deterministic functions, the aggregate is then computed using this dictionary and the remaining data. This scheme is referred to as 'data splitting'. See for instance Dalalyan and Tsybakov [18] and Yang [43]. Leung and Barron [28] proved that data splitting is not necessary in order to aggregate least squares estimators and raised the question of computation of this estimator in high dimension. In this paper we explicit the oracle inequality satisfied by this estimator in the high-dimensional case. For the second procedure, the design may be either random or deterministic. We adapt to the regression framework PAC-Bayesian techniques developed by Catoni [14] in the classification framework to build an estimator satisfying a sparsity oracle inequality for the true excess risk. Even though we do not study the computational aspect in this paper, it should be noted that efficient Monte Carlo algorithms are available to compute these exponential weights estimators for reasonably large dimension p (p ≃ 5000), see in particular the monograph of Marin and Robert [32] for an introduction to MCMC methods. Note also that in a work parallel to ours, Rigollet and Tsybakov [36] consider also an exponential weights procedure with discrete priors and suggest a version of the Metropolis-Hastings algorithm to compute it.
The paper is organized as follows. In Section 2 we define an exponential weights procedure and derive a sparsity oracle inequality in expectation when the design is deterministic. In Section 3, the design can be either deterministic or random. We propose a modification of the first exponential weights procedure for which we can establish a sparsity oracle inequality in probability for the true excess risk. Finally Section 4 contains all the proofs of our results.

Sparsity oracle inequality in expectation
Throughout this section, we assume that the design is deterministic and the noise variables For any J ⊂ {1, . . . , p} and K > 0 define For the sake of simplicity we will write Θ . . , p}) the set of all subsets of {1, . . . , p} containing at most n elements. The aggregatef n is defined as followŝ where λ > 0 is the temperature parameter, π is the prior probability distribution on P({1, . . . , p}), the set of all subsets of {1, . . . , p}, that is, for any J ∈ {1, . . . , p}, π J ≥ 0 and J∈P({1,...,p}) π J = 1. The next result is a reformulation in our context of Theorem 8 of [28].
Proposition 2.1 holds true for any prior π. We suggest using the following prior. Fix α ∈ (0, 1) and define

P. Alquier and K. Lounici
As a consequence, we obtain the following immediate corollary of Proposition 2.1.
Then the aggregatef n = fθ n , with λ = n 4σ 2 and π taken as in (2.5), satisfies This result improves upon [17] which established in Theorem 6 for gaussian noise and deterministic design Note that our bound is faster by a factor log + |θ| 1 . Note also that the bound in the above display grows worse for large values of |θ| 1 . In order to evaluate the performance of these exponential weights procedures, [40] developed a notion of optimal rate of sparse prediction. In particular, [36] established that there exists a numerical constant c * > 0 such that for all estimator T n sup θ∈R p \{0}: The above display combined with Theorem 2.1 shows that the exponential weights procedure (2.3) with the prior (2.5) achieves the optimal rate of sparse prediction for any vector θ satisfying |J(θ)| ≤ rank(A) log(1+ep) .

Sparsity oracle inequality in probability
In Section 2 we assumed the design is deterministic and we established an oracle inequality in expectation with the optimal rate of sparse prediction. We want now to establish an oracle inequality in probability that holds true for deterministic and random design. From now on, the design can be either deterministic or random. We make the following mild assumption:

133
We assume in this section that the noise variables are subgaussian. More precisely we have the following condition.
Assumption 3.1. The noise variables W 1 , . . . , W n are independent and independent of X 1 , . . . , X n . We assume also that there exist two known constants σ > 0 and ξ > 0 such that The estimation method is a version of the Gibbs estimator introduced in [13,14]. Fix K ≥ 1. First we define the prior probability distribution as follows. with π taken as in (2.5).
We are now ready to define our estimator. For any λ > 0 we consider the probability measureρ λ admitting the following density w.r.t. the probability measure m dρ λ dm The aggregatef n is defined as follows f n = fθ n ,θ n =θ n (λ, m) = ΘK θρ λ (dθ). We can now state the main result of this section.
The choice λ = λ * comes from the optimization of a (rather pessimistic) upper bound on the risk R (see Inequality (4.9) in the proof of this theorem, page 142). However this choice is not necessarily the best choice in practice even though it gives the good order of magnitude for λ. The practitioner may use cross-validation to properly tune the temperature parameter. Theorem 3.1 improves upon previous results on the following points: 1) Our oracle inequality is sharp in the sense that the leading factor in front of R(θ) is equal to 1 and we require only max j φ j < ∞ whereas l 1 -penalized empirical risk minimization procedures such as Lasso or Dantzig selector have a leading factor strictly larger than 1 and impose in addition stringent conditions on the dictionary (cf. [7,11]). For instance, [7] imposes the design matrix A = (φ j (X i )) 1≤i≤n,1≤j≤p to be deterministic and to satisfy the following restricted eigenvalue condition for the Lasso where for any ∆ ∈ R p and J ⊂ {1, . . . , p}, we denote by ∆ J the vector in R p which has the same components as ∆ on J and zero coordinates on the complement J c . Assuming in addition that the noise is gaussian N (0, σ 2 ) and taking λ = Aσ log p n , A > 8, [7] proved that the Lassoθ L satisfies with probability at where η, C(η) > 0 and C(η) increases to +∞ as η tends to 0. On the downside, our estimator requires the additional condition |θ| 1 ≤ K. This condition is common in the PAC-bayesian literature. Removing this condition is a difficult problem and does not seem possible with the actual techniques of proof where this condition is needed in order to apply Bernstein's inequality.
2) We establish a sparsity oracle inequality in probability for the integrated risk R(·) whereas previous results on the exponential weights are given in expectation [16,17,19,24,29,36].
3) Unlike mirror averaging or progressive mixture rules, satisfying similar inequalities in expectation, our estimator does not involve an averaging step [18,24,29]. As a consequence, its computational complexity is significantly reduced as compared to those procedures with averaging step. For instance [29] considered the model (1.1) with random design and i.i.d. observations (X 1 , Y 1 ), . . . , (X n , Y n ), n ≥ 2. The studied estimator is the following mirror averaging scheme whereθ 0 = Θ(K) θdΠ and for any 1 ≤ k ≤ n − 1,θ k is defined similarly as These estimators can be implemented for example by MCMC. In this case, computing the integral Θ(K) θρ λ (dθ) is the most time-consuming part of the procedure.
The procedure (3.1)-(3.2) requires computing this integral only once whereas the mirror averaging procedure of [29] requires computing integrals of this form n times. 4) Under the assumption |θ| 1 ≤ K for some absolute constant K and taking ǫ = n −1 we have with probability at least 1 − n −1 for some absolute constant C > 0. In [36] a minimax lower bound in expectation is established for deterministic design and Gaussian noise. A similar result holds in probability with the same proof combined with Theorem 2.7 of [41]. There exists absolute constants c 1 , c 2 > 0 such that inf Tn sup θ∈R p : If s ≤ rank(A) log(1+ep) then we observe that the upper bound in (3.3) is optimal up to the additional logarithmic factor log n. Note however that if p ≥ n 1+δ for some δ > 0, which is relevant with the high-dimensional setting we consider in this paper, then our bound is rate optimal.

Proofs of Section 2
This proof uses an argument from Leung and Barron [28].
Proof of Proposition 2.1. The mapping Y →f n (Y ) △ = (f n (X 1 , Y ), . . . ,f n (X n , Y )) T is clearly continuously differentiable by composition of elementary differentiable functions. For any subset J ⊂ {1, . . . , p} define A J = (φ j (X i )) 1≤i≤n,j∈J , Σ J = 1 n A T J A J , Φ J (·) = (φ j (·)) j∈J and and Σ + J denotes the pseudo-inverse of Σ J . Denote by ∂ i the derivative w.r.t. Y i . Simple computations give Thus we have We have Consider the following estimator of the risk Using an argument based on Stein's identity as in [27] we now prove that

137
We have In view of (4.1) we can apply Fubini's Theorem to the right-hand-side of (4.4).
We obtain under the assumption W ∼ N(0, σ 2 ) that A Similar equality holds for the integral over R − . Thus we obtain Combining (4.2), (4.3) and the above display gives Sincef n (·, Y ) is the expectation of f J (·, Y ) w.r.t. the probability distribution ∝ g · π, we have For the sake of simplicity set f J = f J (·, Y ) andf n =f n (·, Y ). Combining (4.2), the above display and λ n 4σ 2 yieldŝ By definition of g J we have Integrating the above inequality w.r.t. the probability distribution 1 C g ·π (where C = J∈Pn({1,...,p}) g J π J is the normalization factor) and using the fact that , π 0 as well as a convex duality argument (cf., e.g., [20], p. 264) we get for all probability measure π ′ on P({1, . . . , p}). Taking the expectation in the last inequality we get for any π ′ Proof of Theorem 2.1. First note that any minimizer θ ∈ R p of the right-handside in (2.6) is such that |J(θ)| rank(A) n where we recall that A = (φ j (X i )) 1 i n,1 j p . Indeed, for any θ ∈ R p such that |J(θ)| > rank(A) we can construct a vector θ ′ ∈ R p such that f θ = f θ ′ and |J(θ ′ )| rank(A) and the mapping x → x log epα x is nondecreasing on (0, p]. Next for any J ∈ P n ({1, . . . , p}) we have Combining the above display with Proposition 2.1 and our definition of the prior π gives the result.