High-dimensional classification by sparse logistic regression

TL;DR

This paper introduces a new sparse logistic regression method for high-dimensional binary classification, providing theoretical bounds and a computationally feasible estimator that is rate-optimal under certain conditions.

Contribution

It proposes a novel model selection procedure with a complexity penalty related to VC-dimension and extends the Slope estimator for logistic regression to high-dimensional settings.

Findings

Derived non-asymptotic bounds for misclassification risk.

Established the rate-optimality of the extended Slope estimator.

Linked the complexity penalty to VC-dimension of sparse classifiers.

Abstract

We consider high-dimensional binary classification by sparse logistic regression. We propose a model/feature selection procedure based on penalized maximum likelihood with a complexity penalty on the model size and derive the non-asymptotic bounds for the resulting misclassification excess risk. The bounds can be reduced under the additional low-noise condition. The proposed complexity penalty is remarkably related to the VC-dimension of a set of sparse linear classifiers. Implementation of any complexity penalty-based criterion, however, requires a combinatorial search over all possible models. To find a model selection procedure computationally feasible for high-dimensional data, we extend the Slope estimator for logistic regression and show that under an additional weighted restricted eigenvalue condition it is rate-optimal in the minimax sense.

Equations

$$Y_{i}\sim Bin(1,p_{i}),\;\;\;\;\;\ln\frac{p_{i}}{1-p_{i}}=\mbox{\boldmath$\beta$}^{t}{\bf x}_{i}$$

$$\ell(\mbox{\boldmath$\beta$})=\sum_{i=1}^{n}\left\{\mbox{\boldmath$\beta$}^{t}{\bf x}_{i}~{}Y_{i}-\ln(1+\exp(\mbox{\boldmath$\beta$}^{t}{\bf x}_{i})\right\}$$

$$\widehat{\mbox{\boldmath$\beta$}}_{M}=\arg\max_{\widetilde{\mbox{\boldmath$\beta$}}\in{\cal B}_{M}}\sum_{i=1}^{n}\left\{\widetilde{\mbox{\boldmath$\beta$}}^{t}{\bf x}_{i}Y_{i}-\ln\left(1+\exp(\widetilde{\mbox{\boldmath$\beta$}}^{t}{\bf x}_{i})\right)\right\}$$

$$\widehat{M}=\arg\min_{M\in\mathfrak{M}}\left\{\sum_{i=1}^{n}\left(\ln\left(1+\exp(\widehat{\mbox{\boldmath$\beta$}}_{M}^{t}{\bf x}_{i})\right)-\widehat{\mbox{\boldmath$\beta$}}_{M}^{t}{\bf x}_{i}Y_{i}\right)+Pen(|M|)\right\},$$

$$KL(\mathbf{p},{\hat{\mathbf{p}}}_{\hat{M}})=\frac{1}{n}\sum _{i=1}^n\left\{p_i\ln \left(\frac{p_i}{{\hat{p}}_{\hat{M}i}}\right)+(1-p_i)\ln \left(\frac{1-p_i}{1-{\hat{p}}_{\hat{M}i}}\right)\right\},$$

$$Pen(|M|)=c|M|\ln \frac{de}{|M|},|M|=1,\dots ,r-1\mathrm{and}Pen(r)=c\cdot r$$

$$\sup_{\mbox{\boldmath$\beta$}:||\mbox{\boldmath$\beta$}||_{0}\leq d_{0}}EKL({{\bf p}},\widehat{{\bf p}}_{\widehat{M}})\leq C\frac{1}{\delta(1-\delta)}~{}\frac{\min\left(d_{0}\ln\frac{de}{d_{0}},r\right)}{n}$$

$$d_0{\log }_2\left(\frac{2d}{d_0}\right)\le V(\mathcal{C}(d_0))\le 2d_0{\log }_2\left(\frac{de}{d_0}\right)$$

$$\widehat{\eta}_{\widehat{M}}({\bf x})=I\{\widehat{\mbox{\boldmath$\beta$}}_{\widehat{M}}^{t}{\bf x}\geq 0\}$$

$${\mathcal{E}}_X({\hat{\eta }}_{\hat{M}},{\eta }^{\ast })\le \sqrt{2EKL(\mathbf{p},{\hat{\mathbf{p}}}_{\hat{M}})}$$

$$\underset{{\eta }^{\ast }\in \mathcal{C}(d_0)}{\sup }{\mathcal{E}}_X({\hat{\eta }}_{\hat{M}},{\eta }^{\ast })\le C_1\sqrt{\frac{1}{\delta (1-\delta )}\frac{\min \left(d_0\ln \frac{de}{d_0},r\right)}{n}}$$

$$\underset{\tilde{\eta }}{\inf }\underset{{\eta }^{\ast }\in \mathcal{C}(d_0)}{\sup }{\mathcal{E}}_{X_0}(\tilde{\eta },{\eta }^{\ast })\ge C_2\sqrt{\frac{d_0\ln \frac{de}{d_0}}{n}}$$

$$\sqrt{\frac{d_0\ln \left(\frac{de}{d_0}\right)}{n}}\sim \sqrt{\frac{V(\mathcal{C}(d_0))}{n}}$$

$$R_X({\hat{\eta }}_{\hat{M}})-R_X({\eta }^{\ast })=\left(R_X({\hat{\eta }}_{\hat{M}})-R_X({\eta }_L^{\ast })\right)+\left(R_X({\eta }_L^{\ast })-R_X({\eta }^{\ast })\right),$$

$$|p_{i}-1/2|\geq h\;\;{\rm or,\;equivalently,}\;\;|\mbox{\boldmath$\beta$}^{t}{\bf x}_{i}|\geq\ln\left(\frac{1+2h}{1-2h}\right)$$

$$h\le |p_i-1/2|\le \Delta$$

$$\underset{{\eta }^{\ast }\in {\mathcal{C}}_X(d_0,h)}{\sup }{\mathcal{E}}_X({\hat{\eta }}_{\hat{M}},{\eta }^{\ast })\le C_1\min \left(\sqrt{\frac{1-4h^2}{1-4{\Delta }^2}\frac{\min \left(d_0\ln \frac{de}{d_0},r\right)}{n}},\frac{1-4h^2}{1-4{\Delta }^2}\frac{\min \left(d_0\ln \frac{de}{d_0},r\right)}{nh}\right)$$

$$\underset{\tilde{\eta }}{\inf }\underset{{\eta }^{\ast }\in {\mathcal{C}}_{X_0}(d_0,h)}{\sup }{\mathcal{E}}_{X_0}(\tilde{\eta },{\eta }^{\ast })\ge C_2\min \left(\sqrt{\frac{d_0\ln \frac{de}{d_0}}{n}},\frac{d_0\ln \frac{de}{d_0}}{nh}\right)$$

$$\widehat{\mbox{\boldmath$\beta$}}_{Lasso}=\arg\min_{\widetilde{\mbox{\boldmath$\beta$}}}\left\{\sum_{i=1}^{n}\left(\ln\left(1+\exp(\widetilde{\mbox{\boldmath$\beta$}}^{t}{\bf x}_{i})\right)-\widetilde{\mbox{\boldmath$\beta$}}^{t}{\bf x}_{i}Y_{i}\right)+\lambda||\widetilde{\mbox{\boldmath$\beta$}}||_{1}\right\}$$

$$\widehat{\mbox{\boldmath$\beta$}}_{Slope}=\arg\min_{\widetilde{\mbox{\boldmath$\beta$}}}\left\{||{\bf Y}-X\widetilde{\mbox{\boldmath$\beta$}}||_{2}^{2}+\sum_{j=1}^{d}\lambda_{j}|\widetilde{\mbox{\boldmath$\beta$}}|_{(j)}\right\},$$

$$\widehat{\mbox{\boldmath$\beta$}}_{Slope}=\arg\min_{\widetilde{\mbox{\boldmath$\beta$}}}\left\{\sum_{i=1}^{n}\left(\ln\left(1+\exp(\widetilde{\mbox{\boldmath$\beta$}}^{t}{\bf x}_{i})\right)-\widetilde{\mbox{\boldmath$\beta$}}^{t}{\bf x}_{i}Y_{i}\right)+\sum_{j=1}^{d}\lambda_{j}|\widetilde{\mbox{\boldmath$\beta$}}|_{(j)}\right\},$$

$$\kappa (d_0,c_0)=\underset{\mathbf{u}\in \mathcal{S}(d_0,c_0);\mathbf{u}\ne 0}{\min }\frac{||X\mathbf{u}|{|}_2}{||\mathbf{u}|{|}_2}>0$$

$${\lambda }_j=A\frac{c_0+1}{c_0-1}\sqrt{\ln (2d/j)},j=1,\dots ,d$$

$$\sup_{\mbox{\boldmath$\beta$}:||\mbox{\boldmath$\beta$}||_{0}\leq d_{0}}\mathbb{E}KL({{\bf p}},\widehat{{\bf p}}_{Slope})\leq 8A^{2}~{}\frac{c_{0}^{2}}{(c_{0}-1)^{2}}~{}\frac{1}{\delta(1-\delta)}\left(\frac{2\pi+8}{\ln(2d)}+\frac{1}{\kappa^{2}(d_{0},c_{0})}\right)\frac{d_{0}}{n}\ln\left(\frac{2de}{d_{0}}\right)$$

$${\mathcal{E}}_X({\hat{\eta }}_{Slope},{\eta }^{\ast })=O\left(\sqrt{\frac{d_0\ln \frac{de}{d_0}}{n}}\right)$$

$$Y|({\bf X}={\bf x})\sim B(1,p({\bf x})),\;\;\;p({\bf x})=E(Y|{\bf X}={\bf x})=\frac{\exp(\mbox{\boldmath$\beta$}^{t}{\bf x})}{1+\exp(\mbox{\boldmath$\beta$}^{t}{\bf x})}\;\;{\rm and}\;\;{\bf X}\sim q({\bf x}),$$

$$\underset{\tilde{\eta }}{\inf }\underset{{\eta }^{\ast }\in \mathcal{C}(d_0),q}{\sup }\mathcal{E}(\tilde{\eta },{\eta }^{\ast })\ge C\sqrt{\frac{V(\mathcal{C}(d_0))}{n}}\ge \tilde{C}\sqrt{\frac{d_0\ln \frac{de}{d_0}}{n}}$$

$$Pen(|M|)=C|M|\ln \frac{de}{|M|},|M|=1,\dots ,\min (d,n),$$

$$\underset{{\eta }^{\ast }\in \mathcal{C}(d_0)}{\sup }\mathcal{E}({\hat{\eta }}_{\hat{M}},{\eta }^{\ast })\le C\sqrt{\ln \left(\frac{1}{\delta {\lambda }_{\min }(G))}\right)\frac{d_0\ln \frac{de}{d_0}}{n}}$$

$$P\left(|p(\mathbf{X})-1/2|\le h\right)\le C{h}^{\alpha }$$

Full text

Keywords: Complexity penalty; feature selection; high-dimensionality; misclassification excess risk; sparsity; VC-dimension.

1 Introduction

Classification is one of the most important setups in statistical learning and has been studied in various contexts. Theoretical foundations of classification are presented in the books of Devroye, Györfi and Lugosi (1996) and Vapnik (2000), while the surveys of the state-of-the-art can be found in Boucheron, Bousquet and Lugosi (2005) and Giraud (2015, Section 9).

Consider a general (binary) classification with a (high-dimensional) vector of features ${\bf x}\in\mathbb{R}^{d}$ and the outcome class label $Y|{\bf x}\sim Bin(1,p({\bf x}))$. The accuracy of a classifier $\eta$ is defined by a misclassification error $R(\eta)=P(Y\neq\eta({\bf x}))$. It is well-known that $R(\eta)$ is minimized by the Bayes classifier $\eta^{*}({\bf x})=I\{p({\bf x})\geq 1/2\}$. However, the probability function $p({\bf x})$ is unknown and the resulting classifier $\widehat{\eta}({\bf x})$ should be designed from the data $D$: a random sample of $n$ independent observations $({\bf x}_{1},Y_{1}),\ldots,({\bf x}_{n},Y_{n})$. The design points ${\bf x}_{i}$ may be considered as fixed or random. The corresponding (conditional) misclassification error of $\widehat{\eta}$ is $R(\widehat{\eta})=P(Y\neq\widehat{\eta}({\bf x})|D)$ and the goodness of $\widehat{\eta}$ w.r.t. $\eta^{*}$ is measured by the misclassification excess risk ${\cal E}(\widehat{\eta},\eta^{*})=ER(\widehat{\eta})-R(\eta^{*})$.

A common general (nonparametric) approach for finding a classifier $\widehat{\eta}$ from the data is empirical risk minimization (ERM), where minimization of a true misclassification error $R(\eta)$ is replaced by minimization of the corresponding empirical risk $\widehat{R}_{n}(\eta)=\frac{1}{n}\sum_{i=1}^{n}I\{Y_{i}\neq\eta({\bf x}_{i})\}$ over a given class of classifiers. Misclassification excess risk of ERM classifiers has been intensively studied in the literature (see, e.g., Boucheron, Bousquet and Lugosi, 2005 and Giraud, 2015, Section 9 for surveys and references therein). However, ERM can be hardly used directly in practice due to its computational cost and is typically relaxed by some related convex minimization surrogate (e.g., SVM).

Another possibility to obtain $\widehat{\eta}$ is to estimate $p({\bf x})$ from the data by some $\widehat{p}({\bf x})$ and use a plug-in classifier of the form $\widehat{\eta}({\bf x})=I\{\widehat{p}({\bf x})\geq 1/2\}$. A standard approach is to assume some (parametric or nonparametric) model for $p({\bf x})$. In this paper we consider one of the most commonly used models – logistic regression, where it is assumed that p({\bf x})=\frac{\exp(\mbox{\boldmath\beta}^{t}{\bf x})}{1+\exp(\mbox{\boldmath\beta}^{t}{\bf x})} and \mbox{\boldmath\beta}\in\mathbb{R}^{d} is a vector of unknown regression coefficients. The corresponding Bayes classifier is a linear classifier \eta^{*}({\bf x})=I\{p({\bf x})\geq 1/2\}=I\{\mbox{\boldmath\beta}^{t}{\bf x}\geq 0\}. One then estimates $\beta$ from the data by the maximum likelihood estimator (MLE) \widehat{\mbox{\boldmath\beta}}, plugs-in \widehat{\mbox{\boldmath\beta}} (or, equivalently, $\widehat{p}({\bf x})$) and the resulting (linear) classifier is \widehat{\eta}({\bf x})=I\{\widehat{p}({\bf x})\geq 1/2\}=I\{\widehat{\mbox{\boldmath\beta}}^{t}{\bf x}\geq 0\}. Unlike ERM, the MLE \widehat{\mbox{\boldmath\beta}} though not available in the closed form, can be nevertheless obtained numerically by the fast iteratively reweighted least squares algorithm (McCullagh and Nelder, 1989, Section 2.5). Nonparamertric plug-in classifiers were considered in Yang (1999), Koltchinskii and Beznosova (2005), Audibert and Tsybakov (2007).

In the era of “Big Data”, however, the number of features $d$ describing the objects for classification might be very large and even larger than the sample size $n$ (large $d$ small $n$ setups) that raises a severe “curse of dimensionality” problem. Reducing the dimensionality of a feature space by selecting a sparse subset of “significant” features becomes essential. Thus, Bickel and Levina (2004), Fan and Fan (2008) showed that even in simple cases, high-dimensional classification without feature selection might be as bad as just pure guessing.

Nevertheless, unlike model selection in high-dimensional Gaussian regression that has been intensively studied in 2000s (see Birgé and Massart, 2001, 2007; Abramovich and Grinshtein, 2010; Rigollet and Tsybakov, 2011; Verzelen, 2012 among many others), there are much less theoretical results on model/feature selection in classification. Devroy, Györfi and Lugosi (1996, Chapter 18) and Vapnik (2000, Chapter 4) considered selection from a sequence of classifiers within a sequence of classes by penalized ERM with the structural penalty depending on the Vapnik-Chervonenkis (VC) dimension of a class. They established the oracle inequalities and the upper bounds for the misclassification excess risk of the selected classifier but did not provide the lower bound to assess its optimality. See also Boucheron, Bousquet and Lugosi (2005, Section 8) for related penalized ERM approaches and references therein. Recall, however, that a computational cost (even for a given model) is a serious drawback of any ERM-based procedure.

The main goal of the paper is to fill the gap. In particular, we investigate feature selection in sparse logistic regression classification. Although logistic regression is widely used in various classification problems, its rigorous theoretical ground has not been yet properly established. Model selection in a general framework of generalized linear models (GLM) and in logistic regression in particular was studied in Abramovich and Grinshtein (2016). They proposed model selection procedure based on penalized maximum likelihood with a complexity penalty on the model size and investigated the goodness-of-fit of the resulting estimator in terms of the Kullback-Leibler risk. They derived the nonasymptotic bounds for this risk and showed that the resulting estimator is asymptotically minimax and adaptive to the unknown sparsity. In this paper we utilize their approach for classification and consider the corresponding plug-in classifier. In particular, we show that the considered complexity penalty is remarkably related to the VC-dimension of a set of sparse linear classifiers. We establish the non-asymptotic upper bound for misclassification excess risk of the resulting classifier and construct explicitly the design for which it is sharp in the minimax sense. We also show that the excess risk bounds can be improved under the additional low-noise assumption.

Any model selection criterion based on a complexity penalty requires, however, a combinatorial search over all possible models that makes its usefulness problematic for high-dimensional data. A common remedy is to replace the original complexity penalty by a related convex surrogate. The probably most well-known techniques is Lasso. However, it can achieve only the sub-optimal rates under some extra conditions on the design matrix $X$ (van de Geer, 2008). Recently, for Gaussian linear regression Bogdan et al. (2015) proposed a Slope estimator. Bellec, Lecué and Tsybakov (2018) showed that under certain additional conditions on $X$, Slope is rate-optimal for linear regression. We adapt it to the logistic regression (and, in fact, to a general GLM) setup and extend the results of Bellec, Lecué and Tsybakov (2018).

The rest of the paper is organized as follows. In Section 2 we present the model (feature) selection procedure for sparse logistic regression with fixed design based on a general procedure of Abramovich and Grinshtein (2016) and provide the upper bounds for the resulting estimator in terms of Kullback-Leibler risk. In Section 3 we apply it for classification to establish the non-asymptotic upper bound for its misclassification excess risk and derive the corresponding minimax lower bounds. The improvement of the obtained risk bounds under the additional low-noise assumption is given in Section 4. In Section 5 we consider the logistic Slope classifier as a convex surrogate for the proposed feature selection procedure and show that its misclassification excess risk is still rate-optimal under an extra weighted restricted eigenvalue condition on the design matrix $X$. The random design case is considered in Section 6 . Section 7 provides a short real-data example. All the proofs are given in the Appendix.

2 Notation and preliminaries

Consider a sparse logistic regression model

[TABLE]

with deterministic design points ${\bf x}_{i}\in\mathbb{R}^{d}\;\;\;i=1,\ldots,n$, where we assume that the unknown vector of regression coefficients \mbox{\boldmath\beta}\in\mathbb{R}^{d} is sparse.

Let d_{0}=||\mbox{\boldmath\beta}||_{0} be the size of true (unknown) model, where the $l_{0}$ (quasi)-norm of regression coefficients ||\mbox{\boldmath\beta}||_{0} is the number of nonzero entries. Let $X\in\mathbb{R}^{n\times d}$ be the design matrix of rows ${\bf x}_{i}$, $r=rank(X)$ and assume that any $r$ columns of $X$ are linearly independent.

For the model (1) the log-likelihood function is

[TABLE]

Let $\mathfrak{M}$ be the set of all $2^{d}$ possible models $M\subseteq\{1,\ldots,d\}$. For a given model $M$, define {\cal B}_{M}=\{\mbox{\boldmath\beta}\in\mathbb{R}^{d}:\beta_{j}=0\;{\rm if}\;j\not\in M\}. The MLE \widehat{\mbox{\boldmath\beta}}_{M} of $\beta$ is then

[TABLE]

The \widehat{\mbox{\boldmath\beta}}_{M} in (2) is not available in the closed form but can be obtained numerically by the iteratively reweighted least squares algorithm (see McCullagh and Nelder, 1989, Section 2.5). The corresponding MLE for probabilities $p_{i}$ are \widehat{p}_{Mi}=\frac{\exp(\widehat{\mbox{\boldmath\beta}}_{M}^{t}{\bf x}_{i})}{1+\exp(\widehat{\mbox{\boldmath\beta}}_{M}^{t}{\bf x}_{i})}.

Select the model $\widehat{M}$ by the penalized maximum likelihood model selection criterion of the form

[TABLE]

where $Pen(|M|)$ is a complexity penalty on the model size $|M|$. In fact, we can restrict $\mathfrak{M}$ in (3) to models with sizes at most $r$ since for any $\beta$ with ||\mbox{\boldmath\beta}||_{0}>r, there necessarily exists another \mbox{\boldmath\beta}^{\prime} with ||\mbox{\boldmath\beta}^{\prime}||_{0}\leq r such that X\mbox{\boldmath\beta}=X\mbox{\boldmath\beta}^{\prime}.

Within general GLM framework, Abramovich and Grinshtein (2016) investigated the goodness-of-fit of the resulting estimator $\widehat{{\bf p}}_{\widehat{M}}$. They considered the Kullback-Leibler divergence $KL({\bf p},\widehat{{\bf p}}_{\widehat{M}})$ between the data distribution with the true probabilities ${\bf p}=(p_{1},\ldots,p_{n})$ and its empirical distribution generated by $\widehat{{\bf p}}_{\widehat{M}}$ given by

[TABLE]

and measured the accuracy of $\widehat{{\bf p}}_{\widehat{M}}$ by the corresponding Kullback-Leibler risk $EKL({\bf p},\widehat{{\bf p}}_{\widehat{M}})$ (in fact, the Kullback-Leibler divergence $KL(\cdot,\cdot)$ in Abramovich and Grinshtein, 2016 was defined as $n$ times $KL(\cdot,\cdot)$ in this paper).

Assumption** **(A).

Assume that there exists $0<\delta<1/2$ such that $\delta<p_{i}<1-\delta$ or, equivalently, there exists $C_{0}>0$ such that |\mbox{\boldmath\beta}^{t}{\bf x}_{i}|<C_{0} in (1) for all $i=1,\ldots,n$.

Assumption (A) prevents the variances $Var(Y_{i})=p_{i}(1-p_{i})$ to be infinitely close to zero.

Consider a set of models of size at most $d_{0}$, where $1\leq d_{0}\leq r$. Obviously, $|M|\leq d_{0}$ iff ||\mbox{\boldmath\beta}||_{0}\leq d_{0}. Abramovich and Grinshtein (2016) showed that for the complexity penalty

[TABLE]

in (3), where $c>\frac{4}{\delta(1-\delta)}$, under Assumption (A), the upper bound of the Kullnack-Leibler risk is given by

[TABLE]

for some $C>0$. They also derived the corresponding minimax lower bounds for the Kullback-Leibler risk and showed that for weakly-collinear design, the upper bound in (6) is of the optimal order (in the minimax sense).

The above results on the Kullback-Leibler risk can be extended to model selection under additional structural constraints on the set of admissible models $\mathfrak{M}$ (see Secition 4.1 of Abramovich and Grinshtein, 2016).

In what follows we utilize (6) to derive the upper bounds for the misclassification excess risk of the corresponding plug-in classifier \widehat{\eta}_{\widehat{M}}({\bf x})=I\{\widehat{\mbox{\boldmath\beta}}_{\widehat{M}}^{t}{\bf x}\geq 0\}.

To gain more insight into the complexity penalty (5) within classification framework we present the following lemma on the Vapnik-Chervonenskis (VC) dimension of the set of all $d_{0}$-sparse linear classifiers :

Lemma 1**.**

Let {\cal C}(d_{0})=\{\eta({\bf x})=I\{\mbox{\boldmath\beta}^{t}{\bf x}\geq 0\}:\mbox{\boldmath\beta}\in\mathbb{R}^{d},\;||\mbox{\boldmath\beta}||_{0}\leq d_{0}\} be the set of all $d_{0}$-sparse linear classifiers and $V({\cal C}(d_{0}))$ its VC-dimension. Then,

[TABLE]

Thus, the complexity penalty $Pen(|M|)$ in (5) is essentially proportional to the VC-dimension of the corresponding class of $|M|$-sparse linear classifiers ${\cal C}(|M|)$.

3 Misclassification excess risk bounds

In this section we apply the selected model $\widehat{M}$ in (3) for classification and derive the bounds for the corresponding misclassification exceess risk.

We consider first the fixed design. For a given design matrix $X$ the misclassification error of a classifier $\eta$ is $R_{X}(\eta)=\frac{1}{n}\sum_{i=1}^{n}P(Y_{i}\neq\eta({\bf x}_{i}))$. Following our previous arguments define a (linear) plug-in classifier

[TABLE]

and consider its misclassification excess risk ${\cal E}_{X}(\widehat{\eta}_{\widehat{M}},\eta^{*})=ER_{X}(\widehat{\eta}_{\widehat{M}})-R_{X}(\eta^{*})$, where recall that the (ideal) Bayes classifier \eta^{*}({\bf x})=I\{\mbox{\boldmath\beta}^{t}{\bf x}\geq 0\} with the true (unknown) $\beta$ in (1). The remarkable results of Zhang (2004) and Bartlett, Jordan and McAuliffe (2006) establish the relations between the Kullback-Leibler risk $EKL({{\bf p}},\widehat{{\bf p}}_{\widehat{M}})$ and the misclassification excess risk ${\cal E}_{X}(\widehat{\eta}_{\widehat{M}},\eta^{*})$:

[TABLE]

Thus, (6) and (9) imply immediately the following upper bound for ${\cal E}_{X}(\widehat{\eta}_{\widehat{M}},\eta^{*})$:

Theorem 1**.**

Consider a sparse logistic regression model (1) with ||\mbox{\boldmath\beta}||_{0}\leq d_{0}. Let ${\widehat{M}}$ be a model selected in (3) with the complexity penalty (5) and consider the corresponding plug-in classifier $\widehat{\eta}_{\widehat{M}}({\bf x})$ in (8). Under Assumption (A),

[TABLE]

for some $C_{1}>0$, simultaneously for all $1\leq d_{0}\leq r$.

The constant ${\sqrt{\frac{1}{\delta(1-\delta)}}}$ in (10) is a result of the direct application of (6) and (9). It can be improved by establishing the similar relations between misclassification excess risk and other losses rather than Kullback-Leibler in (9) and deriving the corresponding upper bounds for their risks. See, for example, the results and the proof of Theorem 6 below for the random design. In a way, the Kullback-Leibler loss implies the most conservative upper bound (Painsky and Wornell, 2018).

We now show that there exists a design matrix $X_{0}$ for which the upper bound for the misclassification excess risk (10) is essential sharp (up to a probably different constant).

Consider the set of all possible $d_{0}$-sparse linear classifiers ${\cal C}(d_{0})$ defined in Lemma 1 and the case, where a Bayes classifier $\eta^{*}({\bf x})$ is not perfect, that is, $R(\eta^{*})>0$ (aka an agnostic model). Then, the following result holds:

Theorem 2**.**

Consider a $d_{0}$-sparse agnostic logistic regression model (1), where $2\leq d_{0}\log_{2}\left(\frac{2d}{d_{0}}\right)\leq n$.

Then, there exists a design matrix $X_{0}\in\mathbb{R}^{n\times d}$ such that

[TABLE]

for some constant $C_{2}>0$, where the infimum is taken over all classifiers $\tilde{\eta}$ based on the data $(X_{0},{\bf Y})$.

Theorem 2 is a particular case of Theorem 4 from Section 4 below.

The upper and lower bounds established in Theorem 1 and Theorem 2 allow one to derive the asymptotic minimax rate for misclassification excess risk in sparse logistic regression model as $n$ increases. We allow the number of features $d$ to increase with $n$ as well and even faster than $n$ ($d\gg n$ setup). The following immediate Corollary 1 shows that the proposed classifier $\widehat{\eta}_{\widehat{M}}$ is asymptotically minimax in terms of “the worst case” design and adaptive to the unknown sparsity:

Corollary 1**.**

Consider a $d_{0}$-sparse logistic regression agnostic model (1), where $d_{0}$ satisfies $2\leq d_{0}\log_{2}\left(\frac{2d}{d_{0}}\right)\leq n$. Then, as $n$ and $d$ increase, for a fixed $\delta$ in Assumption (A),

The asymptotic minimax misclassification excess risk $\sup_{X}\inf_{\tilde{\eta}}\sup_{\eta^{*}\in{\cal C}(d_{0})}{\cal E}_{X}(\tilde{\eta},\eta^{*})$ is of the order

[TABLE] 2. 2.

The classifier $\widehat{\eta}_{\widehat{M}}$ defined in (8), where the model $\widehat{M}$ was selected by (3) with the complexity penalty (5), attains the minimax rates simultaneously for all $2\leq d_{0}\log_{2}\left(\frac{2d}{d_{0}}\right)\leq n$.

Finally, we note that if the considered logistic regression model is misspecified and the Bayes classifier $\eta^{*}$ is not linear, we still have the following risk decomposition

[TABLE]

where $\eta^{*}_{L}=\arg\min_{\eta\in{\cal C}(d)}R_{X}(\eta)$ is the best (ideal) linear classifier. Our previous arguments can then be applied to the first term in the RHS of (12) representing the estimation error, while the second term is an approximation error and measures the ability of linear classifiers to perform as good as $\eta^{*}$. Enriching the class of classifiers may improve the approximation error but will necessarily increase the estimation error in (12). In a way, it is similar to the variance/bias tradeoff in regression.

4 Tighter risk bounds under low-noise condition

The main challenges for any classifier occur near the the boundary $\{{\bf x}:p({\bf x})=1/2\}$ (equivalently, a hyperplane \mbox{\boldmath\beta}^{t}{\bf x}=0 for the logistic regression model), where it is hard to predict the class label accurately. However, for regions, where $p({\bf x})$ is bounded away from $1/2$ (margin or aka low-noise condition), the bounds for misclassification excess risk established in the previous Section 3 can be improved. Following Massart and Nédélec (2006) introduce the following low-noise assumption:

Assumption** **(B).

Consider the logistic regression model (1) and assume that there exists $0\leq h<1/2$ such that

[TABLE]

for all $i=1,\ldots,n$.

Assumption (B) essentially assumes the existence of the “corridor” of width $2\ln\left(\frac{1+2h}{1-2h}\right)$ that separates the two sets \{{\bf x}_{i}:\mbox{\boldmath\beta}^{t}{\bf x}_{i}>0,\;i=1,\ldots,n\} and \{{\bf x}_{i}:\mbox{\boldmath\beta}^{t}{\bf x}_{i}<0,\;i=1,\ldots,n\}.

For a given design matrix $X$, define {\cal C}_{X}(d_{0},h)=\{\eta:\eta\in{\cal C}(d_{0}),\;|\mbox{\boldmath\beta}^{t}{\bf x}_{i}|\geq\ln\left(\frac{1+2h}{1-2h}\right),\;i=1,\ldots,n\}. Evidently, ${\cal C}_{X}(d_{0},0)={\cal C}(d_{0})$ for any $X$.

Theorem 3 below establishes the upper bound for the misclassification excess risk of the proposed classifier $\widehat{\eta}_{\widehat{M}}$ under the additional low noise Assumption (B):

Theorem 3**.**

Consider a sparse logistic regression model (1), where ||\mbox{\boldmath\beta}||_{0}\leq d_{0}. Assume that there exist $0<h<\Delta<1/2$ such that

[TABLE]

for all $i=1,\ldots,n$.

Let ${\widehat{M}}$ be a model selected in (3) with the complexity penalty (5) and consider the corresponding classifier $\widehat{\eta}_{\widehat{M}}({\bf x})$ in (8). Then, for all $1\leq d_{0}\leq r$,

[TABLE]

for some $C_{1}>0$.

Thus, if the margin parameter $h$ is large enough, namely, $h>\sqrt{\frac{d_{0}\ln\frac{de}{d_{0}}}{n}}$, the misclassification excess risk bound (10) is reduced. The classifier $\widehat{\eta}_{\widehat{M}}({\bf x})$ does not depend on $h$ and the procedure is inherently adaptive to its value.

Similar to the previous Section 3, one can construct a design matrix for which the upper bound (15) is sharp:

Theorem 4**.**

Consider a $d_{0}$-sparse agnostic logistic regression model (1) with $2\leq d_{0}\log_{2}\frac{2d}{d_{0}}\leq n$.

There exists a design matrix $X_{0}\in\mathbb{R}^{n\times d}$ such that under Assumption (B)

[TABLE]

for some $C_{2}>0$.

The design matrix $X_{0}$ is constructed explicitly in the proof of Theorem 4 in the Appendix. Note that Theorem 2 may be viewed as a particular case of Theorem 4 for $h=0$.

5 Logistic Slope classifier

Solving for $\widehat{M}$ in (3) requires generally a combinatorial search over all possible models in $\mathfrak{M}$ that makes the use of complexity penalties to be computationally problematic when the number of features is large. Greedy algorithms (e.g., forward selection) approximate the global solution of (3) by a stepwise sequence of local ones. However, they require strong constraints on the design matrix $X$ that can hardly hold for high-dimensional data. A more reasonable approach is convex relaxation, where the original combinatorial problem is replaced by a related convex surrogate. Thus, for linear-type complexity penalties of the form Pen(|M|)=\lambda|M|=\lambda||\mbox{\boldmath\beta}||_{0}, the celebrated Lasso replaces the $l_{0}$-(quasi) norm by $l_{1}$-norm:

[TABLE]

Assume that all the columns of the design matrix $X$ are normalized to have unit norms. From the results of van de Geer (2008) it follows that under an assumption similar to Assumption (A) and certain extra conditions on $X$, the logistic Lasso with a tuning parameter $\lambda$ of the order $\sqrt{\ln d}$ results in sub-optimal Kullback-Leibler risk $O\left(\frac{d_{0}}{n}\ln d\right)$ and, therefore, sub-optimal misclassification excess risk $O\left(\sqrt{\frac{d_{0}}{n}\ln d}\right)$. For Gaussian regression, Bellec, Lecué and Tsybakov (2018) showed that under certain conditions on $X$, Lasso can achieve the optimal rate with adaptively chosen $\lambda$ by Lepski procedure.

Recently, for Gaussian regression, Bogdan et al. (2015) suggested the Slope estimator – a penalized maximum likelihood estimator with a sorted $l_{1}$-norm penalty defined as follows:

[TABLE]

where $||\cdot||_{2}$ denotes the Euclidean norm in $\mathbb{R}^{n}$, |\widetilde{\mbox{\boldmath\beta}}|_{(1)}\geq\ldots\geq|\widetilde{\mbox{\boldmath\beta}}|_{(d)} are the descendingly ordered absolute values of \widetilde{\mbox{\boldmath\beta}}_{j}’s and $\lambda_{1}\geq\ldots\geq\lambda_{d}>0$ are the tuning parameters. It is a convex minimization problem. Unlike the constant $\lambda$ in Lasso, there is a sequence of decreasing $\lambda_{j}$ in Slope. Bellec, Lecué and Tsybakov (2018) proved that under a weighted restricted eigenvalue condition on $X$ with normalized columns, the quadratic risk of the Slope estimator (17) with $\lambda_{j}=A\sqrt{\ln(2d/j)}$ for a certain constant $A$ is of the (rate-optimal) order $O\left(\frac{d_{0}}{n}\ln(\frac{de}{d_{0}})\right)$.

We will now extend the above results for Slope for logistic regression and, in fact, for a general GLM (see the Appendix C). Naturally modifying the definition of the Slope estimator for the considered logistic regression model (1), define

[TABLE]

where $\lambda_{1}\geq\ldots\geq\lambda_{d}>0$. Note that (18) is also a convex program that makes the logistic Slope estimator computationally feasible for high-dimensional data. The corresponding estimated probabilities \widehat{p}_{Slope,i}=\frac{\exp(\widehat{\mbox{\boldmath\beta}}_{Slope}^{t}{\bf x}_{i})}{1+\exp(\widehat{\mbox{\boldmath\beta}}_{Slope}^{t}{\bf x}_{i})},\;i=1,\ldots,n.

As usual, any convex relaxation requires certain extra conditions on the restricted minimal eigenvalue of the design matrix $X$ over some set of vectors. In particular, similar to Gaussian regression considered in Bellec, Lecué and Tsybakov (2018), we assume the following Weighted Restricted Eigenvalue (WRE) condition for Slope estimator (18) :

Assumption**.**

($WRE(d_{0},c_{0})$** condition)*** Consider the sparse logistic regression model (1) with ||\mbox{\boldmath\beta}||_{0}\leq d_{0}, where the columns of the design matrix $X$ are normalized to have unit norms. Consider the set ${\cal S}(d_{0},c_{0})=\{{\bf u}\in\mathbb{R}^{d}:\sum_{j=1}^{d}\sqrt{\ln(2d/j)}|u|_{(j)}\leq(1+c_{0})||{\bf u}||_{2}~{}\sqrt{\sum_{j=1}^{d_{0}}\ln(2d/j)}\}$ and assume that $X{\bf u}\neq{\bf 0}$ for any ${\bf u}\neq{\bf 0}\in{\cal S}(d_{0},c_{0})$.*

An interesting discussion on the relations between the WRE condition and the restricted eigenvalue condition (RE) required for Lasso is given in Section 8 of Bellec, Lecué and Tsybakov (2018).

Define a restricted minimal eigenvalue $\kappa(d_{0},c_{0})$ as follows :

[TABLE]

Theorem 5**.**

Consider a sparse logistic regression model (1), where ||\mbox{\boldmath\beta}||_{0}\leq d_{0}, the columns of the design matrix $X$ are normalized to have unit norms and, in addition, $X$ satisfies the $WRE(d_{0},c_{0})$ condition for some $c_{0}>1$. Assume that Assumption (A) holds.

Let the tuning parameters

[TABLE]

with the constant $A\geq 20\sqrt{6}$.

Then,

[TABLE]

for all $1\leq d_{0}\leq r$.

Note that $\lambda_{j}$’s in (19) are of the same form as those in Bellec, Lecué and Tsybakov (2018) for Gaussian regression but differ in a constant $A$.

Theorem 5 is a particular case of Theorem 8 for a general GLM (see Appendix C).

Using (9) one immediately gets the corresponding result for the misclassification exceess risk of the logistic Slope classifier:

Corollary 2**.**

Assume all the conditions of Theorem 5 and choose $\lambda_{j}$ according to (19). Consider the logistic Slope classifier \widehat{\eta}_{Slope}({\bf x})=I\{\widehat{\mbox{\boldmath\beta}}^{t}_{Slope}{\bf x}\geq 0\}. Then,

[TABLE]

Thus, the logistic Slope estimator is computationally feasible and yet achieves the optimal rates under the additional $WRE(d_{0},c_{0})$ condition on the design for all but very dense models for which $d_{0}\ln(\frac{de}{d_{0}})>r$ (see Theorem 1).

Furthermore, following the arguments in the proof of Theorem 3, one can show that the bound (21) for ${\cal E}_{X}(\widehat{\eta}_{Slope},\eta^{*})$ may be reduced under the additional low noise Assumption (B).

6 Random design

The results above have been obtained for the fixed design. In machine learning, it is more common to consider classification with random design. In this section we show that our main previous results for the fixed design can be extended for the random design.

Consider the following model:

[TABLE]

where $q(\cdot)$ is a marginal density of ${\bf X}$ with a bounded support ${\cal X}\subset\mathbb{R}^{d}$. By re-scaling we can assume without loss of generality that $||{\bf x}||_{2}\leq 1$ for all ${\bf x}\in{\cal X}$, where recall that $||\cdot||_{2}$ is the Euclidean norm.

We assume that all $X_{j}$ are linearly independent. Hence, the minimal eigenvalue $\lambda_{\min}(G)$ of the matrix $G=E({\bf X}{\bf X}^{t})$ is strictly positive.

Recall that the misclassification excess risk of a classifier $\widehat{\eta}$ designed from a random sample $({\bf X}_{1},Y_{1}),\ldots,({\bf X}_{n},Y_{n})$ from the joint distribution $({\bf X},Y)$ is ${\cal E}(\widehat{\eta}_{\widehat{M}},\eta^{*})=ER(\widehat{\eta})-ER(\eta^{*})$, where the Bayes classifier \eta^{*}({\bf x})=I(\mbox{\boldmath\beta}^{t}{\bf x}\geq 0). Adapting a general result on the minimax lower bound for the misclassification excess risk for random design (see, e.g., Devroye, Györfi and Lugosi, 1996, Chapter 14 and Boucheron, Bousquet and Lugosi, 2005, Section 5.5) for a $d_{0}$-sparse agnostic logistic regression model (22), by Lemma 1 we have:

[TABLE]

Similar to the fixed design setup, consider the penalized maximum likelihood model selection procedure (3) with the complexity penalty

[TABLE]

where the exact choice for the constant $C$ will follow from the proof of Theorem 7 below.

The following Assumption (A1) is a direct analog of Assumption (A) for random design:

Assumption** **(A1).

Assume that there exists $0<\delta<1/2$ such that $\delta<p({\bf x})<1-\delta$ or, equivalently, there exists $C_{0}>0$ such that |\mbox{\boldmath\beta}^{t}{\bf x}|<C_{0} in (22) for all ${\bf x}\in{\cal X}$.

The following Theorem 6 extends the results of Theorem 1 for random design:

Theorem 6**.**

Consider a sparse logistic regression model (22), where ||\mbox{\boldmath\beta}||_{0}\leq d_{0}.

Let ${\widehat{M}}$ be a model selected in (3) with the complexity penalty (23) and consider the corresponding plug-in classifier $\widehat{\eta}_{\widehat{M}}({\bf x})$ in (8). Under Assumption (A1),

[TABLE]

for some positive $C>0$ , simultaneously for all $1\leq d_{0}\leq\min(d,n)$.

Thus, the classifier $\widehat{\eta}_{\widehat{M}}({\bf x})$ is adaptively rate-optimal (in the minimax sense) for random design as well.

Theorem 6 is a particular case of Theorem 7 stated below.

We should note that similar upper bounds can be obtained for model selection by penalized ERM utilizing general results of Devroye, Györfi and Lugosi (1996, Chapter 18) and Vapnik (2000, Chapter 4) on structural penalties depending on a VC-dimension and applying Lemma 1 for their adaptation to sparse logistic regression. See also Boucheron, Bousquet and Lugosi, (2005, Section 8) for related ERM approaches and references therein. Recall, however, that a computational cost is a crucial drawback of any ERM-based procedure.

The misclassification excess risk of $\widehat{\eta}_{\widehat{M}}({\bf x})$ can be again improved under the low-noise condition which can even be formulated in a more general form for random design (Mammen and Tsybakov, 1999 and Tsybakov, 2004):

Assumption** **(B1).

Assume that there exist $C>0$ and $\alpha\geq 0$ such that

[TABLE]

for all $0<h<h^{*}$, where $h^{*}<1/2$.

Assumption (B) in Section 4 for the fixed design can be viewed as a limiting case $\alpha=\infty$.

Theorem 7**.**

Consider a sparse logistic regression model (22), where ||\mbox{\boldmath\beta}||_{0}\leq d_{0}. Let ${\widehat{M}}$ be a model selected in (3) with the complexity penalty (23) and consider the corresponding plug-in classifier $\widehat{\eta}_{\widehat{M}}({\bf x})$ in (8).

Under Assumptions (A1) and (B1), there exists $C>0$ such that

[TABLE]

for all $1\leq d_{0}\leq\min(d,n)$.

Theorem 6 (no low-noise assumption) corresponds to the extreme case $\alpha=0$. For another extreme case $\alpha=\infty$ (complete separation from 1/2), the upper bound (25) is $O\left(\frac{d_{0}\ln\frac{de}{d_{0}}}{n}\right)$ similar to the results of Section 4 for the fixed design.

The rates (25) can be reduced further under additional conditions on the support ${\cal X}$ and the density $q({\bf x})$ using the arguments of Koltchinskii and Beznosova (2005) and Audibert and Tsybakov (2007) but this is beyond the scope of the paper.

To extend the results of Theorem 5 for Slope estimator for random design one needs the $WRE(d_{0},c_{0})$ condition to be held with high probability. It evidently depends on the marginal distribution $q({\bf x})$. Thus, Bellec, Lecué and Tsybakov (2018, Theorem 8.3) showed that it is satisfied for multivariate Gaussian and even sub-Gaussian ${\bf X}$ when $d_{0}\ln^{2}(de/d_{0})\leq cn$ for some constant $c>0$, and under mild conditions on the covariance matrix.

7 Numerical example

We now demonstrate the performance of the proposed feature selection and classification procedures on a short numerical real-data study.

The online marketing data contains various information about 500 customers registered at the site: their personal data (e.g., gender, country of living, etc.) and purchase activities during the past period (e.g., frequencies and types of purchases, purchase amounts, currencies, etc.). Overall, there were 61 explanatory variables. Given such data, one of the main goals is to predict customers who are about to become inactive in order to incentivize them to remain by various discounts. The output then is a binary variable indicating whether a customer was still active during the next time period or not.

The data was randomly split into training ($n_{1}=400$) and test ($n_{2}=100$) sets. The number of possible features $d=61$ was too large to perform a complete combinatorial search for penalized maximum likelihood model selection procedure (3) with the complexity penalty (5). Instead we used its forward selection version, logistic Lasso and logistic Slope classifiers. The corresponding tuning constants were chosen by 5-fold cross-validation and the resulting three classifiers were then applied to the test set.

The best misclassification rate was achieved by Slope (17%), followed by Lasso (19%) and forward selection (22%). In addition, we compared the sizes of the models selected by the three classifiers. The conservative forward selection procedure yielded a very sparse model with only 2 predictors – the time since last purchase and purchase amount. On the other hand, Lasso with the CV-chosen tuning parameter is known to tend to select too many variables (see, e.g., Bühlmann and van de Geer, 2011, Section 2) and resulted in the model of size 13 by adding 11 other variables. The Slope classifier with a decreasing sequence of tuning parameters commonly implies even larger models (32 in the considered example). Note, however, that prediction and model identification are two different problems and, in particular, the choices of tuning parameters for them should be different.

Acknowledgments

The work was supported by the Israel Science Foundation (ISF), grants ISF-820/13 and ISF-589/18. The authors would like to thank Noga Alon for his help in the proof of Lemma 1, Alexander Tsybakov for valuable remarks and Roi Granot for the real-data example.

Appendix

Throughout the proofs we use various generic positive constants, not necessarily the same each time they are used even within a single equation.

Appendix A: Proof of Lemma 1

Denote for brevity $V=V({\cal C}(d_{0}))$. For any fixed subset of $d_{0}$ $\beta_{j}$’s the VC of the corresponding set of $d_{0}$-dimensional linear classifiers is known to be $d_{0}$ (e.g., Giraud, 2015, Exercise 9.5.2). Then, by Sauer’s lemma the maximal number of different labelling of $V$ points in $\mathbb{R}^{d_{0}}$ that such set of classifiers can produce is $\sum_{k=0}^{d_{0}}\binom{V}{k}\leq\left(\frac{Ve}{d_{0}}\right)^{d_{0}}$ (see, e.g., Giraud, 2015, Section 9.2.2). The overall number of different labelling is, therefore, $\binom{d}{d_{0}}\sum_{k=0}^{d_{0}}\binom{V}{k}$, and by the definition of $V({\cal C}(d_{0}))$ we have

[TABLE]

that implies an upper bound $V\leq 2~{}d_{0}\log_{2}\left(\frac{de}{d_{0}}\right)$.

On the other hand, take $k=\log_{2}(2d/d_{0})$ and let $K$ be the $k\times 2^{k-1}$ matrix whose columns are all possible vectors with entries $\pm 1$ and the first entry $1$. Note that $d_{0}2^{k-1}=d$. Let $W$ be the $d_{0}k\times d$ block-wise matrix consisting of $d_{0}\times d_{0}$ blocks, each being a $k\times 2^{k-1}$ matrix, where the diagonal matrices are copies of $K$, while all others are zero matrices. Thus, $W$ has $d_{0}k=d_{0}\log_{2}(2d/d_{0})$ rows. It is easy to verify that these rows are shattered by half-spaces whose supporting vectors $w$ have a single non-zero $\pm 1$ entry in each of the $d_{0}$ blocks and, therefore, $V\geq d_{0}\log_{2}(2d/d_{0})$.

Appendix B: Tighter bounds for low-noise condition

B1: Proof of Theorem 3

Assumption (14) obviously implies Assumption (A) with $\delta=1/2-\Delta$. In addition, under (14), $Var(Y_{i})=p_{i}(1-p_{i})\leq(1/2-h)(1/2+h)=(1-4h^{2})/4$. Hence, adapting the results of Abramovich and Grinshtein (2016) on Kullback-Leibler risk in general GLM framework for logistic regression, the upper bound (6) for $EKL({{\bf p}},\widehat{{\bf p}}_{\widehat{M}})$ can be improved:

[TABLE]

and, therefore, from (9) we have

[TABLE]

On the other hand, we can adapt the general Theorem 3 of Bartlett, Jordan and McAuliffe (2006) for $\psi(f)=(1/2)\left((1-f)\ln(1-f)+(1+f)\ln(1+f)\right)\geq f^{2}/2$ corresponding to the Kullback-Leibler risk (Zhang, 2004, Section 3.5), $\alpha=1$ corresponding to (14) and $c=1/(2h)$ to get

[TABLE]

Applying (26) implies then

[TABLE]

B2: Proof of Theorem 4

For any $\tilde{\eta}$ and $\eta^{*}\in{\cal C}(d_{0},h)$ we have

[TABLE]

for any $X$.

As we have mentioned, the worst case scenario for classification is when $p_{i}=1/2\pm h$ or, equivalently, |\mbox{\boldmath\beta}^{t}{\bf x}_{i}|=\ln\left(\frac{1+2h}{1-2h}\right). Let $V=d_{0}\log_{2}(2d/d_{0})$. In the proof of Lemma 1 we constructed explicitly the matrix $W_{V\times d}$ whose rows ${\bf w}_{1},\ldots,{\bf w}_{V}$ are shattered by ${\cal C}(d_{0})$. Then, for any ${\bf p}=\{\frac{1}{2}\pm h\}^{V}$ there exists \mbox{\boldmath\beta}\in\mathbb{R}^{d} such that ||\mbox{\boldmath\beta}||_{0}\leq d_{0} and \mbox{\boldmath\beta}^{t}{\bf w}_{i}=\ln\frac{p_{i}}{1-p_{i}}=\pm\ln\frac{1+2h}{1-2h} for all $i=1,\ldots,V$. Define also the corresponding binary vector ${\bf b}$ with b_{i}=I\{\mbox{\boldmath\beta}^{t}{\bf w}_{i}\geq 0\}, that is, $b_{i}=1$ if $p_{i}=\frac{1}{2}+h$ and $b_{i}=0$ if $p_{i}=\frac{1}{2}-h$. Obviously, the set of all ${\bf b}$’s is a hypercube $H^{V}=\{0,1\}^{V}$.

Define now a $n\times d$ design matrix $X_{0}$ with $\varkappa$ rows of ${\bf w}_{1}$, $\varkappa$ rows of ${\bf w}_{2}$, …, $\varkappa$ rows of ${\bf w}_{V-1}$ and the remaining $n-(V-1)\varkappa$ rows of ${\bf w}_{V}$, where an integer $1\leq\varkappa\leq\lfloor\frac{n}{V-1}\rfloor$ will be defined later.

The proof will now follow the general scheme of the proof of Theorem 4 of Massart and Nédélec (2006) but with necessary modifications for the fixed design.

For any ${\bf p}\in\{\frac{1}{2}\pm h\}^{V}$ and the corresponding ${\bf b}\in H^{V}$ define an $n$-dimensional indicator vector $\eta_{\bf b}=(\underbrace{b_{1},\ldots,b_{1}}_{\varkappa},\ldots,\underbrace{b_{V-1},\ldots,b_{V-1}}_{\varkappa},\underbrace{b_{V},\ldots,b_{V}}_{n-(V-1)\varkappa})$ and let $\tilde{C}_{X_{0}}(d_{0},h)=\{\eta_{\bf b},\;{\bf b}\in H^{V}\}$. By its design, \tilde{{\cal C}}_{X_{0}}(d_{0},h)\subseteq\{\eta:\eta\in{\cal C}(d_{0}),\;|\mbox{\boldmath\beta}^{t}{\bf x}_{0i}|=\ln\frac{1+2h}{1-2h},\;i=1,\ldots,n\}\subseteq C_{X_{0}}(d_{0},h).

Hence, we can reduce the minimax risk over the entire ${\cal C}_{X_{0}}(d_{0},h)$ to $\tilde{{\cal C}}_{X_{0}}(d_{0},h)$:

[TABLE]

Furthermore, for a given $\tilde{\eta}$, define $\tilde{\eta}^{*}=\arg\min_{\eta\in\widetilde{{\cal C}}(d_{0},h)}||\tilde{\eta}-\eta||_{1}$. Then, for any $\eta^{*}\in\tilde{{\cal C}}_{X_{0}}(d_{0},h)$ we have

[TABLE]

and, therefore, from (27)-(29)

[TABLE]

where $\tilde{\bf b},{\bf b}^{*}\in H^{V}$ are the binary vectors corresponding to $\tilde{\eta}^{*}$ and $\eta^{*}$ respectively (see above).

By a simple calculus one can verify that the square Hellinger distance $H^{2}\left(Bin(1,\frac{1}{2}+h),Bin(1,\frac{1}{2}-h)\right)$ between two Bernoulli distributions $Bin(1,\frac{1}{2}+h)$ and $Bin(1,\frac{1}{2}-h)$ is $1-\sqrt{1-4h^{2}}$. For any ${\bf b}\in H^{V}$ and the corresponding $\eta_{\bf b}$ define ${\bf p}_{\bf b}\in\mathbb{R}^{n}$ as follows: $p_{{\bf b}i}=\frac{1}{2}+h$ if $\eta_{{\bf b}i}=1$ and $p_{{\bf b}i}=\frac{1}{2}-h$ if $\eta_{{\bf b}i}=0,\;i=1,\ldots,\varkappa(V-1)$, and $p_{{\bf b}i}=0,\;i=\varkappa(V-1)+1,\ldots,n$. Then, for any ${\bf b}_{1},{\bf b}_{2}\in H^{V}$ and the corresponding ${\bf p}_{{\bf b}_{1}}$ and ${\bf p}_{{\bf b}_{2}}$ we have

[TABLE]

Hence, applying the version of Assouad’s lemma given in Lemma 7 of Barron, Birgé and Massart (1999) yields

[TABLE]

that together with (30) implies

[TABLE]

Consider two cases.

Case 1. $h\leq\frac{1}{6}$.

For $h\geq\sqrt{\frac{V-1}{18n}}$, apply (31) for $\varkappa=\lfloor\frac{1}{18h^{2}}\rfloor$ (note that $2\leq\varkappa\leq\lfloor\frac{n}{V-1}\rfloor$), to get

[TABLE]

For $h<\sqrt{\frac{V-1}{18n}}$, one can follow all the above arguments for $\tilde{h}=\sqrt{\frac{V-1}{18n}}$ and the corresponding $\varkappa=\lfloor\frac{n}{V-1}\rfloor$ to have

[TABLE]

Case 2. $h>\frac{1}{6}$.

Set $\varkappa=1$ and note that $C_{X_{0}}(d_{0},\frac{1}{2})\subseteq C_{X_{0}}(d_{0},h)$ for any $0\leq h\leq\frac{1}{2}$. Hence, (30) implies

[TABLE]

By obvious combinatoric calculus, for any (binary) vector $\tilde{\bf b}$, $\frac{1}{2^{V}}\sum_{j=1}^{2^{V}}P(\tilde{b}_{i}\neq b_{ji})=\frac{1}{2}$ for any $i$ and, therefore,

[TABLE]

for large $h>\frac{1}{6}$ (in fact, larger than any fixed $h_{0}$).

Appendix C: Slope estimator for a general GLM

Consider a GLM setup with a response variable $Y$ and a set of $d$ predictors $x_{1},...,x_{d}$. We observe a series of independent observations $({\bf x}_{i},Y_{i}),\;i=1,\ldots,n$, where the design points ${\bf x}_{i}\in\mathbb{R}^{p}$ are deterministic. The distribution $f_{\theta_{i}}(y)$ of $Y_{i}$ belongs to a (one-parameter) natural exponential family with a natural parameter $\theta_{i}$ and a scaling parameter $a$:

[TABLE]

The function $b(\cdot)$ is assumed to be twice-differentiable. In this case $\mathbb{E}(Y_{i})=b^{\prime}(\theta_{i})$ and $Var(Y_{i})=ab^{\prime\prime}(\theta_{i})$. To complete GLM we assume the canonical link \theta_{i}=\mbox{\boldmath\beta}^{t}{\bf x}_{i} or, equivalently, in the matrix form, \mbox{\boldmath\theta}=X\mbox{\boldmath\beta}, where $X_{n\times p}$ is the design matrix and \mbox{\boldmath\beta}\in\mathbb{R}^{p} is a vector of the unknown regression coefficients. The logistic regression (1) is a particular case of a general GLM (32) for the Bernoulli distribution $Bin(1,p_{i})$, where the natural parameter is $\theta=\ln\frac{p}{1-p},b(\theta_{0})=\ln(1+e^{\theta})$ and $a=1$.

Following Abramovich and Grinshtein (2016) assume the extended version of Assumption (A) for GLM :

Assumption** **(A’).

Assume that $\theta_{i}\in\Theta$, where the parameter space $\Theta\subseteq\mathbb{R}$ is a closed (finite or infinite) interval. 2. 2.

Assume that there exist constants $0<{\cal L}\leq{\cal U}<\infty$ such that the function $b^{\prime\prime}(\cdot)$ satisfies the following conditions:

1. (a)

$\sup_{t\in\mathbb{R}}b^{\prime\prime}(t)\leq{\cal U}$** 2. (b)

$\inf_{t\in\Theta}b^{\prime\prime}(t)\geq{\cal L}$**

Conditions on $b^{\prime\prime}(\cdot)$ in Assumption (A’) are intended to exclude two degenerate cases, where the variance $Var(Y)$ is infinitely large or small. They also ensure strong convexity of $b(\cdot)$ over $\Theta$. For the binomial distribution, ${\cal U}=1/4$ and Assumption (A’) reduces to Assumption (A) with ${\cal L}=\delta(1-\delta)$.

Recall that the Slope estimator is a penalized maximum likelihood with an ordered $l_{1}$-norm penalty and, therefore, defined for a GLM as follows:

[TABLE]

for $\lambda_{1}\geq\cdots\geq\lambda_{d}>0$. The corresponding Kullback-Leibler risk

[TABLE]

where \mbox{\boldmath\theta}=X\mbox{\boldmath\beta} and \widehat{\mbox{\boldmath\theta}}_{Slope}=X\widehat{\mbox{\boldmath\beta}}_{Slope} (see Abramovich and Grinshtein, 2016).

Theorem 8**.**

Consider a GLM (32), where ||\mbox{\boldmath\beta}||_{0}\leq d_{0}, the columns of the design matrix $X$ are normalized to have unit norms and $X$ satisfies the $WRE(d_{0},c_{0})$ condition for some $c_{0}>1$. Assume that Assumption (A’) holds.

Let

[TABLE]

in (33) with the constant $A\geq 40\sqrt{6}$.

Then, simultaneously for all \mbox{\boldmath\beta}\in\mathbb{R}^{d} such that ||\mbox{\boldmath\beta}||_{0}\leq d_{0},

[TABLE]

for any $0<\Delta<1$. 2. 2.

[TABLE]

Proof.

Since \widehat{\mbox{\boldmath\beta}}_{Slope} is the minimizer of (33),

[TABLE]

From (29) of Abramovich and Grinshtein (2016) one has

[TABLE]

(recall that the Kullback-Leibler divergence $KL(\cdot,\cdot)$ in Abramovich and Grinshtein, 2016 was defined as $n$ times $KL(\cdot,\cdot)$ in this paper). Thus,

[TABLE]

Let {\bf u}=\widehat{\mbox{\boldmath\beta}}_{Slope}-\mbox{\boldmath\beta}. Applying Lemma A.1 of Bellec, Lecué and Tsybakov (2018) with $\tau=0$ implies

[TABLE]

Consider now the first term of the RHS in (38). Since the distribution of $Y$ belongs to the exponential family with the bounded variance $ab^{\prime\prime}(\theta)\leq a{\cal U}$ (Assumption (A’)), a centered zero mean random variable $Y-b^{\prime}(\theta)$ is sub-Gaussian with the scale factor $\sqrt{a{\cal U}}$, that is, $\mathbb{E}e^{t(Y-b^{\prime}(\theta))}\leq e^{a{\cal U}t^{2}/2}$ and, therefore, $\mathbb{E}e^{(Y-b^{\prime}(\theta))^{2}/(6{\cal U}a)}\leq e$. Applying Theorem 9.1 of Bellec, Lecué and Tsybakov (2018) (adapted to our normalization conditions on the columns of $X$) yields

[TABLE]

with probability at least $1-\Delta$.

Set

[TABLE]

and

[TABLE]

The proof will now go along the lines of the proof of Theorem 6.1 of Bellec, Lecué and Tsybakov (2018) for Gaussian regression with necessary adaptations to GLM and different normalization conditions on the columns of $X$.

To prove (36) consider two cases.

$Case~{}1.~{}\tilde{H}({\bf u})\leq G({\bf u})$. In this case

[TABLE]

and, therefore, combining (35) and (38)-(42) with probability at least $1-\Delta$ yields

[TABLE]

for any $\epsilon>0$.

Lemma 1 of Abramovich and Grinshtein (2016) established the equivalence of the Kullback-Leibler divergence KL(\mbox{\boldmath\theta},\widehat{\mbox{\boldmath\theta}}_{Slope}) and the squared quadratic norm ||\widehat{\mbox{\boldmath\theta}}_{Slope}-\mbox{\boldmath\theta}||^{2} under Assumption (A’):

[TABLE]

Hence, taking $\epsilon={\cal L}/(2a)$ in (43) after a straightforward calculus yields

[TABLE]

with probability at least $1-\Delta$.

$Case~{}2.~{}\tilde{H}({\bf u})>G({\bf u})$. Using the definition of $\lambda_{j}$’s in (35) and (38)-(42), with probability at least $1-\Delta$ we have

[TABLE]

The KL(\mbox{\boldmath\theta},\widehat{\mbox{\boldmath\theta}}_{Slope})\geq 0 and, therefore, the RHS of (46) is necessarily positive. Thus,

[TABLE]

and, therefore, by $WRE(d_{0},c_{0})$ condition, (46) implies

[TABLE]

for any $\epsilon>0$. Taking $\epsilon={\cal L}/(4a)$ and exploiting the equivalence between KL(\mbox{\boldmath\theta},\widehat{\mbox{\boldmath\theta}}_{Slope}) and ||\widehat{\mbox{\boldmath\theta}}_{Slope}-\mbox{\boldmath\theta}||_{2}^{2} in (44) imply that with probability at least $1-\Delta$,

[TABLE]

where we used the definition (35) of $\lambda_{j}$’s and the upper bound $\sum_{j=1}^{d_{0}}\ln(2d/j)\leq d_{0}\ln(2ed/d_{0})$ (see, e.g., (2.7) of Bellec, Lecué and Tsybakov, 2018).

To prove the second statement (37) of the theorem denote $C^{*}=8A^{2}\frac{c_{0}^{2}}{(c_{0}-1)^{2}}$ and note that

[TABLE]

Then, by integrating, (36) and (47) after a straightforward calculus yield

[TABLE]

∎

Appendix D: Proof of Theorem 7

We first introduce several notations. Let $||g||_{L_{2}}=(\int_{\cal X}g^{2}({\bf x})d{\bf x})^{1/2}$ be a standard $L_{2}$-norm of a function $g$ and $||g||_{L_{2}(q)}=(\int_{\cal X}g^{2}({\bf x})q({\bf x})d{\bf x})^{1/2}$ be the $L_{2}$-norm of $g$ weighted by the marginal distribution $q$ of ${\bf X}$. In addition, the $L_{\infty}$-norm $||g||_{\infty}=\sup_{{\bf x}\in{\cal X}}|g({\bf x})|$.

Applying Theorem 3 of Bartlett, Jordan and McAlliffe (2006) for $\psi(g)=g^{2}/2$, Assumption (B1) implies that there exists $C>0$ such that

[TABLE]

Furthermore, let f_{\mbox{\boldmath\beta}}({\bf x},y) be the joint distribution of $({\bf X},Y)$ for a given $\beta$, i.e. f_{\mbox{\boldmath\beta}}({\bf x},y)=p({\bf x})^{y}(1-p({\bf x}))^{1-y}q({\bf x}), where p({\bf x})=\frac{\exp\{\mbox{\boldmath\beta}^{t}{\bf x}\}}{1+\exp\{\mbox{\boldmath\beta}^{t}{\bf x}\}}. Consider the square Hellinger distance $H^{2}(Bin(1,p_{1}),Bin(1,p_{2}))$ between two Bernoulli distributions with success probabilities $p_{1}$ and $p_{2}$. It is easy to verify that $H^{2}(Bin(1,p_{1}),Bin(1,p_{2}))\geq\frac{(p_{1}-p_{2})^{2}}{4(1-\delta)}\geq\frac{(p_{1}-p_{2})^{2}}{2}$. Then for the square Hellinger distance d^{2}_{H}(f_{\mbox{\boldmath\beta}_{1}},f_{\mbox{\boldmath\beta}_{2}}) between f_{\mbox{\boldmath\beta}_{1}} and f_{\mbox{\boldmath\beta}_{2}} we have

[TABLE]

and from (48) it is, therefore, sufficient to bound the Hellinger risk Ed^{2}_{H}(f_{\widehat{}\mbox{\boldmath\beta}_{\widehat{M}}},f_{\mbox{\boldmath\beta}}).

We will show that the penalty (23) falls within a general class of penalties considered in Yang and Barron (1998) and then apply their Theorem 1 to find an upper bound for Ed^{2}_{H}(f_{\widehat{}\mbox{\boldmath\beta}_{\widehat{M}}},f_{\mbox{\boldmath\beta}}).

Using the standard inequality $\ln(1+t)\leq t$, under Assumption (A1) we have

[TABLE]

Define \rho(f_{\mbox{\boldmath\beta}_{1}},f_{\mbox{\boldmath\beta}_{2}})=||\ln f_{\mbox{\boldmath\beta}_{2}}-\ln f_{\mbox{\boldmath\beta}_{1}}||_{\infty}. Thus,

[TABLE]

For a given model $M$ consider the set of coefficients ${\cal B}_{M}$ defined in Section 2. One can easily verify that under Assumption (A1), for any \mbox{\boldmath\beta}_{1},\mbox{\boldmath\beta}_{2}\in{\cal B}_{M} and the corresponding $p_{1}({\bf x}),p_{2}({\bf x})$

[TABLE]

for any ${\bf x}\in{\cal X}$.

In particular, (51) implies

[TABLE]

where recall that $G=E({\bf X}{\bf X}^{t})$ and $\lambda_{min}(G)>0$ is its minimal eigenvalue.

For each \mbox{\boldmath\beta}_{0}\in{\cal B}_{M} consider the corresponding Hellinger ball {\cal H}_{f_{\mbox{\boldmath\beta}_{0}},r}=\{f_{\mbox{\boldmath\beta}}:d_{H}(f_{\mbox{\boldmath\beta}},f_{\mbox{\boldmath\beta}_{0}})\leq r,\;\mbox{\boldmath\beta}\in{\cal B}_{M}\}. From (49) and (52) it then follows that if f_{\mbox{\boldmath\beta}}\in{\cal H}_{f_{\mbox{\boldmath\beta}_{0}},r}, the corresponding \mbox{\boldmath\beta}\in{\cal B}_{M} lies in the Euclidean ball {\cal B}_{\mbox{\boldmath\beta}_{0},r^{\prime}}=\{\mbox{\boldmath\beta}\in\mathbb{R}^{|M|}:||\mbox{\boldmath\beta}-\mbox{\boldmath\beta}_{0}||_{2}\leq r^{\prime}\} with $r^{\prime}=\frac{\sqrt{2}r}{\delta(1-\delta)\sqrt{\lambda_{min}(G)}}$.

Furthermore, for any $||{\bf x}||_{2}\leq 1$, (51) and Cauchy–-Schwarz inequality imply that \left|p_{2}({\bf x})-p_{1}({\bf x})\right|\leq\frac{1}{4}||\mbox{\boldmath\beta}_{2}-\mbox{\boldmath\beta}_{1}||_{2} and, therefore, by (50)

[TABLE]

Let N({\cal B}_{\mbox{\boldmath\beta}_{0},r^{\prime}},l_{2},\epsilon) be the $\epsilon$-covering number of {\cal B}_{\mbox{\boldmath\beta}_{0},r^{\prime}} w.r.t. $l_{2}$-distance. It is well-known that N({\cal B}_{\mbox{\boldmath\beta}_{0},r^{\prime}},l_{2},\epsilon)\leq\left(1+\frac{2r^{\prime}}{\epsilon}\right)^{|M|}\leq\left(\frac{3r^{\prime}}{\epsilon}\right)^{|M|} for any $\epsilon<r^{\prime}$.

Thus, for the $\epsilon$-covering number N({\cal H}_{f_{\mbox{\boldmath\beta}_{0}},r},\rho,\epsilon) of {\cal H}_{f_{\mbox{\boldmath\beta}_{0}},r} w.r.t. the distance \rho(f_{\mbox{\boldmath\beta}_{1}},f_{\mbox{\boldmath\beta}_{2}}), from (53) we have

[TABLE]

The considered family of sparse logistic regression models satisfies then Assumption 1 of Yang and Barron (1998) with $A_{M}=\frac{c}{\delta^{2}(1-\delta)\sqrt{\lambda_{\min}(G)}}$ for some $c>0$ and $m_{M}=|M|$. Apply now their Theorem 1 for a penalized maximum likelihood model selection procedure (3) with a complexity penalty $Pen(|M|)=C_{1}~{}m_{M}\ln A_{M}+C_{2}\cdot C_{M}\leq\tilde{C}_{1}\ln\left(\frac{1}{\delta\lambda_{\min}(G)}\right)~{}|M|+C_{2}|M|\ln\frac{de}{|M|}$, where $C_{M}=|M|\ln\frac{de}{|M|}$, and the exact positive constants $C_{1}$ and $C_{2}$ are given in the paper. Thus,

[TABLE]

To complete the proof note that one can always find a constant $C$ in the penalty (23) such that the resulting $Pen(|M|)=C\ln\left(\frac{1}{\delta\lambda_{\min}(G)}\right)|M|\ln\frac{de}{|M|}\geq\tilde{C}_{1}\left(\frac{1}{\delta\lambda_{\min}(G)}\right)~{}|M|+C_{2}|M|\ln\frac{de}{|M|}$.