Pac-bayesian bounds for sparse regression estimation with exponential weights

Introduction

The estimator $\hat{\theta}^{0}_{n}$ is called oracle estimator since the set of indices of the nonzero coordinates of $\overline{\theta}$ is unknown in practice. The issue is now to build an estimator, when the set of nonzero coordinates of $\overline{\theta}$ is unknown, with statistical performances close to that of the oracle estimator $\hat{\theta}^{0}_{n}$.

A possible approach is to consider solutions of penalized empirical risk minimization problems:

where the penalization $pen(\theta)$ is proportional to the number of nonzero components of $\theta$ such as for instance AIC, $C_{p}$ and BIC criteria . Bunea, Tsybakov and Wegkamp established for the BIC estimator $\hat{\theta}_{n}^{BIC}$ the following non-asymptotic sparsity oracle inequality. For any $\epsilon>0$ there exists a constant $C(\epsilon)>0$ such that for any $p\geqslant 2,n\geqslant 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 non-convex combinatorial optimization problems.

Simultaneously, the PAC-Bayesian approach for regression estimation was developed by Audibert and Alquier , based on previous works in the classification context by Catoni , Mc Allester , Shawe-Taylor and Williamson , see also Zhang in the context of density estimation. This framework is very well adapted for studying the excess risk $R(\cdot)-R(\overline{\theta})$ 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 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(\cdot)-R(\overline{\theta})$.

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 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 and Yang . Leung and Barron 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 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\simeq 5000$), see in particular the monograph of Marin and Robert for an introduction to MCMC methods. Note also that in a work parallel to ours, Rigollet and Tsybakov 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 $W_{1},\ldots,W_{n}$ are i.i.d. gaussian $N(0,\sigma^{2})$.

For any $J\subset\{1,\ldots,p\}$ and $K>0$ define

For the sake of simplicity we will write $\Theta_{K}=\Theta_{K}(\{1,\ldots,p\})$.

For any subset $J\subset\{1,\ldots,p\}$ define

where $r(\theta)=\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-f_{\theta}(X_{i}))^{2}=\|Y-f_{\theta}\|_{n}^{2}$ with $Y=(Y_{1},\ldots,Y_{n})^{T}$. Denote by $\mathcal{P}_{n}(\{1,\ldots,p\})$ the set of all subsets of $\{1,\ldots,p\}$ containing at most $n$ elements. The aggregate $\hat{f}_{n}$ is defined as follows

where $\lambda>0$ is the temperature parameter, $\pi$ is the prior probability distribution on $\mathcal{P}(\{1,\ldots,p\})$, the set of all subsets of $\{1,\ldots,p\}$, that is, for any $J\in\{1,\ldots,p\}$, $\pi_{J}\geq 0$ and $\sum_{J\in\mathcal{P}(\{1,\ldots,p\})}\pi_{J}=1$.

The next result is a reformulation in our context of Theorem 8 of .

Assume that the noise variables $W_{1},\ldots,W_{n}$ are i.i.d. $N(0,\sigma^{2})$. Then the aggregate $\hat{\theta}_{n}$ defined by (2.3) with $0<\lambda\leqslant\frac{n}{4\sigma^{2}}$ satisfies

Proposition 2.1 holds true for any prior $\pi$. We suggest using the following prior. Fix $\alpha\in(0,1)$ and define

As a consequence, we obtain the following immediate corollary of Proposition 2.1.

Assume that the noise variables $W_{1},\ldots,W_{n}$ are i.i.d. $N(0,\sigma^{2})$. Then the aggregate $\hat{f}_{n}=f_{\hat{\theta}_{n}}$, with $\lambda=\frac{n}{4\sigma^{2}}$ and $\pi$ taken as in (2.5), satisfies

This result improves upon which established in Theorem 6 for gaussian noise and deterministic design

where $\log_{+}x=\max\left\{\log x,0\right\}$ and we recall that $A=(\phi_{j}(X_{i}))_{1\leq i\leq n,1\leq j\leq p}$. Note that our bound is faster by a factor $\log_{+}|\theta|_{1}$. Note also that the bound in the above display grows worse for large values of $|\theta|_{1}$.

In order to evaluate the performance of these exponential weights procedures, developed a notion of optimal rate of sparse prediction. In particular, established that there exists a numerical constant $c^{*}>0$ such that for all estimator $T_{n}$

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:

We assume in this section that the noise variables are subgaussian. More precisely we have the following condition.

The noise variables $W_{1},\ldots,W_{n}$ are independent and independent of $X_{1},\ldots,X_{n}$. We assume also that there exist two known constants $\sigma>0$ and $\xi>0$ such that

The estimation method is a version of the Gibbs estimator introduced in . Fix $K\geq 1$. First we define the prior probability distribution as follows. For any $J\subset\{1,\ldots,p\}$ let ${\rm u}_{J}$ denote the uniform measure on $\Theta_{K+1}(J)$. We define

We can now state the main result of this section.

The choice $\lambda=\lambda^{*}$ 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 4.9). However this choice is not necessarily the best choice in practice even though it gives the good order of magnitude for $\lambda$. 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(\bar{\theta})$ is equal to $1$ and we require only $\max_{j}\|\phi_{j}\|<\infty$ 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. ). For instance, imposes the design matrix $A=(\phi_{j}(X_{i}))_{1\leq i\leq n,1\leq j\leq p}$ to be deterministic and to satisfy the following restricted eigenvalue condition for the Lasso

where $\eta,C(\eta)>0$ and $C(\eta)$ increases to $+\infty$ as $\eta$ tends to .

On the downside, our estimator requires the additional condition $|\bar{\theta}|_{1}\leq 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(\cdot)$ whereas previous results on the exponential weights are given in expectation .

3) Unlike mirror averaging or progressive mixture rules, satisfying similar inequalities in expectation, our estimator does not involve an averaging step . As a consequence, its computational complexity is significantly reduced as compared to those procedures with averaging step. For instance considered the model (1.1) with random design and i.i.d. observations $(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$, $n\geq 2$. The studied estimator is the following mirror averaging scheme

4) Under the assumption $|\bar{\theta}|_{1}\leq K$ for some absolute constant $K$ and taking $\epsilon=n^{-1}$ we have with probability at least $1-n^{-1}$

for some absolute constant $C>0$. In 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 . There exists absolute constants $c_{1},c_{2}>0$ such that

Proofs

This proof uses an argument from Leung and Barron .

The mapping $Y\rightarrow\hat{f}_{n}(Y)\stackrel{{\scriptstyle\triangle}}{{=}}(\hat{f}_{n}(X_{1},Y),\ldots,\hat{f}_{n}(X_{n},\break Y))^{T}$ is clearly continuously differentiable by composition of elementary differentiable functions. For any subset $J\subset\{1,\ldots,p\}$ define $A_{J}=(\phi_{j}(X_{i}))_{1\leq i\leq n,j\in J}$, $\Sigma_{J}=\frac{1}{n}A^{T}_{J}A_{J}$, $\Phi_{J}(\cdot)=(\phi_{j}(\cdot))_{j\in J}$ and

and $\Sigma_{J}^{+}$ denotes the pseudo-inverse of $\Sigma_{J}$. Denote by $\partial_{i}$ the derivative w.r.t. $Y_{i}$. Simple computations give

Consider the following estimator of the risk

Using an argument based on Stein’s identity as in we now prove that

Combining (4.2), (4.3) and the above display gives

Since $\hat{f}_{n}(\cdot,Y)$ is the expectation of $f_{J}(\cdot,Y)$ w.r.t. the probability distribution $\propto g\cdot\pi$, we have

For the sake of simplicity set $f_{J}=f_{J}(\cdot,Y)$ and $\hat{f}_{n}=\hat{f}_{n}(\cdot,Y)$. Combining (4.2), the above display and $\lambda\leqslant\frac{n}{4\sigma^{2}}$ yields

Integrating the above inequality w.r.t. the probability distribution $\frac{1}{C}g\cdot\pi$ (where $C=\sum_{J\in\mathcal{P}_{n}(\{1,\ldots,p\})}g_{J}\pi_{J}$ is the normalization factor) and using the fact that

as well as a convex duality argument (cf., e.g., , p. 264) we get

for all probability measure $\pi^{\prime}$ on $\mathcal{P}(\{1,\ldots,p\})$. Taking the expectation in the last inequality we get for any $\pi^{\prime}$

Next for any $J\in\mathcal{P}_{n}(\{1,\ldots,p\})$ we have

Combining the above display with Proposition 2.1 and our definition of the prior $\pi$ gives the result. ∎

2 Proof of Theorem 3.1

We state below a version of Bernstein’s inequality useful in the proof of Theorem 3.1. See Proposition 2.9 page 24 in , more precisely Inequality (2.21).

Let $T_{1}$, …, $T_{n}$ be independent real valued random variables. Let us assume that there is two constants $v$ and $w$ such that

For any $\theta\in\Theta_{K+1}$ define the random variables

Note that these variables are independent. We have

where we have used in the last line Pythagore’s Theorem to prove $\|f_{\theta}-f_{\bar{\theta}}\|^{2}=R(\theta)-R(\bar{\theta})$. Next we have, for any integer $k\geq 3$, that

with $w:=8(\xi+[\|f\|_{\infty}+L(K+1/2)])L(K+1/2)$.

Next, for any $\lambda\in(0,n/w)$ and $\theta\in\Theta_{K+1}$, applying Lemma 4.1 with $\zeta=\lambda/n$ gives

Set $C=8\left(\sigma^{2}+[\|f\|_{\infty}+L(K+1/2)]^{2}\right)$. For the sake of simplicity let us put

For any $\varepsilon>0$ the last display yields

Integrating w.r.t. the probability distribution ${\rm m}(\cdot)$ we get

Combining the last two displays we obtain

Now, using Lemma 1.1.3 in Catoni we obtain that

Combining (4.8) and (4.6) with a union bound argument gives

where $\mathcal{M}_{+}^{1}(\Theta_{K+1})$ is the set of all probability measures over $\Theta_{K+1}$.

Now for any $\delta\in(0,1]$ taking $\rho$ as the uniform probability measure on the set $\{t\in\Theta(J(\bar{\theta})):|t-\bar{\theta}|_{1}\leq\delta\}\subset\Theta_{K+1}(J(\bar{\theta}))$ gives

Taking $\delta=n^{-1}$ and the inequality $\log{p\choose|J(\bar{\theta})|}\leq|J(\bar{\theta})|\log\frac{pe}{|J(\bar{\theta})|}$ gives

where we replaced $\gamma$ and $\beta$ by their definitions, see (4.5) and (4.7). Taking now $\lambda=n/(2\mathcal{C}_{1})$ (where we recall that $\mathcal{C}_{1}=C\vee w$) in (4.9) gives

where we have used that $1-\frac{\lambda C}{2(n-w\lambda)}\geq 1/2$ and $1+\frac{\lambda C}{2(n-w\lambda)}\leq 3/2$. ∎

References