Sufficient Dimension Reduction for Classification

TL;DR

This paper introduces the maximal mean variance (MMV) method for high-dimensional classification, which is model-free, consistent, and effective even when the number of variables exceeds the sample size.

Contribution

The paper proposes the MMV approach, a novel sufficient dimension reduction technique tailored for high-dimensional classification, with proven consistency and asymptotic properties.

Findings

MMV outperforms existing methods in simulations.

The method maintains effectiveness when p > n.

Real data analysis confirms practical utility.

Abstract

We propose a new sufficient dimension reduction approach designed deliberately for high-dimensional classification. This novel method is named maximal mean variance (MMV), inspired by the mean variance index first proposed by Cui, Li and Zhong (2015), which measures the dependence between a categorical random variable with multiple classes and a continuous random variable. Our method requires reasonably mild restrictions on the predicting variables and keeps the model-free advantage without the need to estimate the link function. The consistency of the MMV estimator is established under regularity conditions for both fixed and diverging dimension (p) cases and the number of the response classes can also be allowed to diverge with the sample size n. We also construct the asymptotic normality for the estimator when the dimension of the predicting vector is fixed. Furthermore, our method works pretty well when n < p. The surprising classification efficiency gain of the proposed method is demonstrated by simulation studies and real data analysis.

Tables

Table 1. Table 1: Average classification error (percentage)

Method	MMV+LDA		LDA
	$p=50$
Model I	9.95 (3.43)		24.15 (6.21)
	$p=200$
Model I	13.87 (4.42)		19.83 (6.33)

Table 2. Table 2: Average classification error (percentage)

Method	MMV+Logistic regression		Logistic regression
	$p=20$
Model II	9.33 (3.36)		13.31 (4.16)
	$p=50$
Model II	14.85 (4.36)		31.71 (6.27)

Table 3. Table 3: Average classification error (percentage)

Method	MMV+SVM		SVM		MMV+KKN		KKN
	$p=50$
Model III	16.65 (3.61)		20.49 (4.01)		17.80 (3.75)		34.52 (4.51)
Model IV	16.66 (4.09)		25.22 (4.46)		18.33 (4.65)		22.75 (5.43)
	$p=200$
Model III	25.18 (4.05)		27.08 (4.45)		25.63 (3.99)		42.21 (4.40)
Model IV	16.20 (5.02)		19.35 (5.89)		14.09 (4.06)		22.70 (6.45)

Table 4. Table 4: Average classification error (percentage)

	MMV+LDA		LDA
	11.24 (1.16)		16.74 (3.00)
	MMV+SVM		SVM
	12.40 (1.52)		16.53 (2.41)
	MMV+KKN		KKN
	14.37 (2.14)		18.73 (1.43)

Equations

$$MV(Z|Y)={\mathrm{E}}_Z[{\mathrm{Var}}_Y(F(Z|Y))]=\sum _{r=1}^R{p}_r\int [F_r(z)-F(z){]}^{2}dF(z)$$

$$\widehat{MV}(Z|Y):=M{V}_n(Z|Y)=\frac{1}{n}\sum _{r=1}^R\sum _{i=1}^n{\hat{p}}_r[\hat{F}(Z_i)-{\hat{F}}_r(Z_i){]}^2$$

$${\mathbf{\beta }}_{01}=\arg \underset{{\mathbf{\beta }}_1}{\max }MV({\mathbf{\beta }}_1^{\mathrm{T}}\mathbf{X}|Y)$$

$${\mathbf{\beta }}_{0k}=\arg \underset{{\mathbf{\beta }}_k}{\max }MV({\mathbf{\beta }}_k^{\mathrm{T}}\mathbf{X}|Y)$$

$$S_{Y|\mathbf{X}}\subseteq \mathrm{span}\{{\mathbf{\beta }}_{01},\dots ,{\mathbf{\beta }}_{0d}\},$$

$$\mathrm{span}\{{\mathbf{\beta }}_{01},\dots ,{\mathbf{\beta }}_{0k}\}⊉S_{Y|\mathbf{X}}.$$

$$\hat{F}(z):={\hat{F}}_h(z)={\int }_{-\infty }^z{\hat{f}}_h(u)\mathrm{d}u=n^{-1}\sum _{i=1}^n{\int }_{-\infty }^z{K}_h(Z_i-u)\mathrm{d}u$$

$${\hat{F}}_r(z):={\hat{F}}_{hr}(z)=n_r^{-1}\sum _{j=1}^{n_r}{\int }_{-\infty }^z{K}_{h_r}(Z_j-u)\mathrm{d}u$$

$${\mathbf{X}}_i\sim N(Y_i\cdot \mathbf{\mu },\mathbf{\Sigma }),$$

$$\mathbf{w}\propto \mathbf{\Omega }\mathbf{\mu }.$$

$$Y=f({\mathbf{\beta }}_1^{\mathrm{T}}\mathbf{X},\dots ,{\mathbf{\beta }}_k^{\mathrm{T}}\mathbf{X},\epsilon )(k<p),$$

$$Y\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}{\mathbf{X}}|{\mathbf{B}}{\mathbf{X}}$$

$${\mathbf{B}}\boldsymbol{\gamma}=\mathbf{0},\quad\boldsymbol{\gamma}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}{\mathbf{B}}{\mathbf{X}},$$

$$L^{(1)}({\mathbf{\theta }}_1)=:L^{(1)}({\mathbf{\beta }}_1,{\lambda }_1)=MV({\mathbf{\beta }}_1^{\mathrm{T}}\mathbf{X}|Y)+{\lambda }_1({\mathbf{\beta }}_1^{\mathrm{T}}{\mathbf{\beta }}_1-1),$$

$$L_{nh}^{(1)}({\mathbf{\theta }}_1)=:L_{nh}^{(1)}({\mathbf{\beta }}_1,{\lambda }_1)=M{V}_n({\mathbf{\beta }}_1^{\mathrm{T}}\mathbf{X}|Y)+{\lambda }_1({\mathbf{\beta }}_1^{\mathrm{T}}{\mathbf{\beta }}_1-1),$$

$$L^{(2)}({\mathbf{\theta }}_2)=:L^{(2)}({\mathbf{\beta }}_2,{\mathbf{\lambda }}_2)=MV({\mathbf{\beta }}_2^{\mathrm{T}}\mathbf{X}|Y)+{\lambda }_{21}({\mathbf{\beta }}_2^{\mathrm{T}}{\mathbf{\beta }}_2-1)+{\lambda }_{22}({\mathbf{\beta }}_{01}^{\mathrm{T}}{\mathbf{\beta }}_2),$$

$$L_{nh}^{(2)}({\mathbf{\theta }}_2)=:L_{nh}^{(2)}({\mathbf{\beta }}_2,{\mathbf{\lambda }}_2)=M{V}_n({\mathbf{\beta }}_2^{\mathrm{T}}\mathbf{X}|Y)+{\lambda }_{21}({\mathbf{\beta }}_2^{\mathrm{T}}{\mathbf{\beta }}_2-1)+{\lambda }_{22}({\hat{\mathbf{\beta }}}_1^{\mathrm{T}}{\mathbf{\beta }}_2),$$

$$\underset{{\mathbf{\beta }}_1\in \partial B({\kappa }_1)\cap {\Gamma }_1}{\sup }MV({\mathbf{\beta }}_1^{\mathrm{T}}\mathbf{X}|Y)<MV({\mathbf{\beta }}_{01}^{\mathrm{T}}\mathbf{X}|Y).$$

$$\underset{{\mathbf{\beta }}_2\in \partial B({\kappa }_2)\cap {\Gamma }_2\cap {\rm Y}}{\sup }MV({\mathbf{\beta }}_2^{\mathrm{T}}\mathbf{X}|Y)<MV({\mathbf{\beta }}_{02}^{\mathrm{T}}\mathbf{X}|Y).$$

$$|\frac{\partial L(\mathbf{\theta })}{\partial {\theta }_i\partial {\theta }_j\partial {\theta }_k}|\le M_1<\infty |\frac{\partial L(\mathbf{\theta })}{\partial {\theta }_i\partial {\theta }_j\partial {\theta }_k}-\frac{\partial L_{nh}(\mathbf{\theta })}{\partial {\theta }_i\partial {\theta }_j\partial {\theta }_k}|\le M_2<\infty$$

$$\mathbf{X}={\mathbf{\beta }}_{1}Y+\mathbf{\Delta }\mathbf{ϵ},$$

$$Y=\mathrm{I}\left(1/\{1+\exp ({\mathbf{\beta }}_1^{\mathrm{T}}\mathbf{X})\}\ge 0.5\right),$$

$$Y=\mathrm{I}\left({\mathbf{\beta }}_1^{\mathrm{T}}X/\{0.5+({\mathbf{\beta }}_2^{\mathrm{T}}X+1.5{)}^2\}+0.2ϵ\ge 0\right),$$

$$Y=\mathrm{I}\left(({\mathbf{\beta }}_1^{\mathrm{T}}\mathbf{X}{)}^2+({\mathbf{\beta }}_2^{\mathrm{T}}\mathbf{X}{)}^2+0.2ϵ\ge 1\right),$$

$$\begin{array}{rl}F(Z)& =p_1\cdot F_1(Z)+p_{-1}\cdot F_{-1}(Z)\\ & =p_1\cdot \Phi ((Z-{\mathbf{\beta }}^{\mathrm{T}}\mathbf{\mu })/\sqrt{{\mathbf{\beta }}^{\mathrm{T}}\mathbf{\Sigma }\mathbf{\beta }})+p_{-1}\cdot \Phi ((Z+{\mathbf{\beta }}^{\mathrm{T}}\mathbf{\mu })/\sqrt{{\mathbf{\beta }}^{\mathrm{T}}\mathbf{\Sigma }\mathbf{\beta }}),\end{array}$$

$$\begin{array}{rl}MV(\mathbf{\beta })& =p_1\int [F_1(z)-F(z){]}^2\mathrm{d}F(z)+p_{-1}\int [F_{-1}(z)-F(z){]}^2\mathrm{d}F(z)\\ & =\left(p_1{p}_{-1}^2+p_{-1}{p}_1^2\right)\left\{p_1\int [F_1(z)-F_{-1}(z){]}^2\mathrm{d}{F}_1(z)+p_{-1}\int [F_1(z)-F_{-1}(z){]}^2\mathrm{d}{F}_{-1}(z)\right\}\\ & =p_1{p}_{-1}\{p_1\int {\left[\Phi (t)-\Phi (t+2{\mathbf{\beta }}^{\mathrm{T}}\mathbf{\mu }/\sqrt{{\mathbf{\beta }}^{\mathrm{T}}\mathbf{\Sigma }\mathbf{\beta }})\right]}^2\mathrm{d}\Phi (t)\\ & +p_{-1}\int {\left[\Phi (t)-\Phi (t-2{\mathbf{\beta }}^{\mathrm{T}}\mathbf{\mu }/\sqrt{{\mathbf{\beta }}^{\mathrm{T}}\mathbf{\Sigma }\mathbf{\beta }})\right]}^2\mathrm{d}\Phi (t)\}.\end{array}$$

$$L=({\mathbf{\beta }}^{\mathrm{T}}\mathbf{\Sigma }\mathbf{\beta }{)}^{-1/2}{\mathbf{\beta }}^{\mathrm{T}}\mathbf{\mu }+\lambda ({\mathbf{\beta }}^{\mathrm{T}}\mathbf{\beta }-1)+\pi {\mathbf{\beta }}^{\mathrm{T}}{\mathbf{\Sigma }}^{-1}\mathbf{\mu }.$$

$$\frac{\partial L}{\partial \mathbf{\beta }}=({\mathbf{\beta }}^{\mathrm{T}}\mathbf{\Sigma }\mathbf{\beta }{)}^{-1/2}{\mathbf{\mu }}^{\mathrm{T}}-{\mathbf{\beta }}^{\mathrm{T}}\mathbf{\mu }({\mathbf{\beta }}^{\mathrm{T}}\mathbf{\Sigma }\mathbf{\beta }{)}^{-3/2}{\mathbf{\beta }}^{\mathrm{T}}\mathbf{\Sigma }+2\lambda {\mathbf{\beta }}^{\mathrm{T}}+\pi {\mathbf{\mu }}^{\mathrm{T}}{\mathbf{\Sigma }}^{-1}=0.$$

$$\begin{array}{rl}0& =\frac{\partial L}{\partial \mathbf{\beta }}\mathbf{\beta }\\ & =({\mathbf{\beta }}^{\mathrm{T}}\mathbf{\Sigma }\mathbf{\beta }{)}^{-1/2}{\mathbf{\mu }}^{\mathrm{T}}\mathbf{\beta }-({\mathbf{\beta }}^{\mathrm{T}}\mathbf{\Sigma }\mathbf{\beta }{)}^{-1/2}{\mathbf{\beta }}^{\mathrm{T}}\mathbf{\mu }+2\lambda {\mathbf{\beta }}^{\mathrm{T}}\mathbf{\beta }+\pi {\mathbf{\mu }}^{\mathrm{T}}{\mathbf{\Sigma }}^{-1}\beta \\ & =2\lambda ,\end{array}$$

$$0=({\mathbf{\beta }}^{\mathrm{T}}\mathbf{\Sigma }\mathbf{\beta }){\mathbf{\mu }}^{\mathrm{T}}-{\mathbf{\beta }}^{\mathrm{T}}\mathbf{\mu }{\mathbf{\beta }}^{\mathrm{T}}\mathbf{\Sigma }+({\mathbf{\beta }}^{\mathrm{T}}\mathbf{\Sigma }\mathbf{\beta }{)}^{3/2}\pi {\mathbf{\mu }}^{\mathrm{T}}{\mathbf{\Sigma }}^{-1}.$$

Full text

Keywords: Classification; Mean variance index; Sufficient dimension reduction; Consistency.

MSC2010 subject classifications: Primary 62H30; secondary 62G05

1 Introduction

Sufficient dimension reduction fits into what is currently quite a hot area in research of high dimensional data. Large quantities of related articles and studies have appeared in recent decades. However, most of the literature focuses on the regression problem where the responsible variable $Y$ is a continuous scalar, while little is designed specially for the problem of classification with a categorical response. The slice-based methods, including but not limited to the seminal sliced inverse regression (SIR; Li (1991)), sliced average variance estimation (SAVE; Cook and Weisberg (1991)), directional regression (DR; Li and Wang (2007)) and sliced regression (SR; Wang and Xia (2008)), can be naturally applied to the classification problem with the slices determined directly by the classes of the response. It seems to work nicely but the number of the slices is restricted by the number of the classes which can be problematic when there are only a few categories. More specifically, when the response is a binary variable, the number of the slices is imposed as $2$ and the number of effective dimension-reduction directions is correspondingly forced to be $1$, which would directly reduce the accuracy of classification. Moreover, almost all of the above methods require the linearity condition or constant covariance condition, or both, which are difficult to verify in practice, and the results may be misleading if these conditions are violated. Other popular sufficient dimension reduction methods, like minimum average variance estimation (MAVE; Xia et al., (2002)), inverse regression (IR; Cook and Ni (2005)) and distance covariance based sufficient dimension reduction (DC; Sheng and Yin (2013), Sheng and Yin (2016)) either require the responsible variable to be continuous or treat the response as a numeric variable. Therefore, from the perspective of sufficient dimension reduction, it’s necessary to do some work deliberately into the case of the categorical responsible variable or in other words, the classification.

More importantly, from the perspective of classification itself, dimension reduction is of importance for constantly emerging high-dimensional classification problems considering the “curse of dimensionality” which appears in most classification approaches, such as the frequently used linear and quadratic classifiers, support vectors machines, k-nearest neighbours, decision trees, neural networks and new methods like distance weighted discrimination (Marron, Todd and Ahn (2007), Marron (2015) and Wang and Zou (2018)) and so on. The usual practice is to conduct projection or variable selection to reduce dimensionality as a first step. The projection methods have been widely applied to classification for the gene expression data. Related research includes but is not limited to principle component analysis in Zou, Hastie and Tibshirani (2006), Bair et al.  (2006) and Shao et al.  (2014); sliced inverse regression in Bura and Pfeiffer (2003) and Antoniadis et al.  (2003); and partial least squares in Boülesteix (2004). Of note is the fact that sufficient dimension reduction is indeed a kind of projection methods where all the information related to classification is preserved. Variable selection is another line of dimension reduction approach. It can deal with classification problems with extremely large dimensionality. See, for example, the nearest shrunken centroids method in Tibshirani et al.  (2002), the features annealed independence rules in Fan and Fan (2008) and the latest research, including Andrews and McNicholas (2014), Stefanski, Wu and White (2014), and Partovi and Davison (2015), among others. Furthermore, Cui, Li and Zhong (2015) used a novel mean variance (MV) index to implement model-free feature screening for high dimensional discriminant analysis, which is related closely to the proposed method.

In addition to the classification methods aforementioned, there is also a great deal of other research work on classification. To name but a few, Dettling and Bühlmann (2003) and Buhlmann and Yu (2003) studied boosting under logit loss, and $L_{2}$ loss in the high-dimensional setting, respectively; Greenshtein and Ritov (2004) introduced the concept of persistence, which is weaker than consistency and pays more attention to misclassification error; Donoho and Jin (2008) employed higher criticism thresholding for feature screening when the useful features are both rare and weak; and following their work, Fan, Jin and Yao (2013) proposed a two-stage classification procedure based on innovated thresholding and high criticism thresholding in the sparse Gaussian graphical model.

In this article, we propose a novel sufficient dimension reduction approach – maximum mean variance(MMV), designed deliberately for a high-dimensional classification problem based on the mean variance index first proposed by Cui, Li and Zhong (2015). This method is not slice-based, thus it circumvents the restriction on slice number. Moreover, the approach does not require the linearity condition or constant variance condition, nor does it require any special distributions on the predicting vector ${\mathbf{X}}$, ${\mathbf{X}}|Y$ or $Y|{\mathbf{X}}$, which is essential in the methods of Zhu and Zeng (2006), Cook and Forzani (2009), Cook and Li (2009), Bura and Forzani (2015), Bura, Duarte and Forzani (2016) and Zhang, Chen and Zhou (2018). In addition, our method keeps the model-free advantage without the need to estimate the link function. These advantages broaden the scope of the application of our method. The consistency of the MMV estimator is established for both fixed and diverging dimension cases and asymptotic normality is constructed for the case of fixed dimension. The relationship between MMV and classification is more than the usual stepwise heuristics of dimension reduction first and classification next, which is elaborated upon by taking the examples of linear discriminant analysis (LDA) and index models. Moreover, the asymptotic theory of MMV estimator is quite challenging to set up, because the empirical MV index includes the kernel estimation of conditional and unconditional distribution functions, thus it can not be directly expressed by the sum of independent and identically distributed random variables.

The rest of the paper is organised as follows. In Section 2, we propose a new sufficient dimension reduction approach (MMV) to high-dimensional classification. In Section 3, we elaborate on the delicate relationship between the MMV method and classification. Consistency and the asymptotic normality of the MMV estimator are studied in Section 4. Several simulation studies together with numerical comparisons and a real data example are conducted to illustrate the efficiency and priority of the proposed method in Section 5. Section 6 concludes the article and the technical proofs are deferred to the Appendix.

2 Method

2.1 Mean variance index

Let $Y$ be a categorical response with $R$ classes $\{y_{1},y_{2},\dots,y_{R}\}$ and $Z$ be a continuous covariate. The Mean Variance index (Cui, Li and Zhong (2015)) is defined by:

[TABLE]

where $F(z|Y)=\mathbb{P}(Z\leq z|Y)$, $F(z)=\mathbb{P}(Z\leq z)$, $F_{r}(z)=\mathbb{P}(Z\leq z|Y=y_{r})$ and $p_{r}=\mathbb{P}(Y=y_{r})$ for $r=1,\dots,R$. It has been verified that $MV(Z|Y)=0$ if and only if $Y$ and $Z$ are independent. It is worth noting that the MV index characterizes both linear and nonlinear correlations between categorical variable $Y$ and continuous variable $Z$.

Let $\{(Y_{i},Z_{i}):1\leq i\leq n\}$ be an i.i.d random sample of size $n$. Let $\hat{F}(Z)$ and $\hat{F}_{r}(Z)$ be some sample estimations of $F(Z)$ and $F_{r}(Z)$. The MV index can be estimated by:

[TABLE]

where $\hat{p}_{r}=1/n\sum_{i=1}^{n}\mathrm{I}\{Y_{i}=y_{r}\}$ with $\mathrm{I}(\cdot)$ being the indicator function. Cui, Li and Zhong (2015) used the empirical distributions of $Z$ and $Z|Y$ as their sample estimators in a screening procedure.

2.2 Maximum Mean Variance method for sufficient dimension reduction

Classification is a crucial statistical problem which has attracted interest for decades. Let $Y$ be the categorial response defined above and ${\mathbf{X}}\in\mathbb{R}^{p}$ be continuous predictors. In this paper, we propose a novel sufficient reduction approach designed deliberately for high-dimensional classification problems based on the mean variance index. The idea is to find a few linear combinations (or indexes) of original predictors that contribute to classification without a loss of information, so that these derived indexes can then be utilized for classification. This is achieved through a stepwise maximization procedure of the MV index. Recall that $MV(Z|Y)=0$ if and only if $Z$ and $Y$ are statistically independent. Thanks to this property, the MV index is used for marginal feature screening in discriminant analysis (Cui, Li and Zhong (2015)). Our novel idea is to abandon this and, on the contrary, we seek a $\boldsymbol{\beta}\in\mathbb{R}^{p}$ such that $MV(\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}}|Y)$ achieves its maximum under some constraints. This is why we named this method Maximum Mean Variance. Hereafter, we refer to the MMV together with the following classification approach as ”MMV+.”. A sequential algorithm is elaborated as follows: we find the first linear combination of the predictors from

[TABLE]

subject to $\boldsymbol{\beta}_{1}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\beta}_{1}=1$. Then the $k$th linear combination can be calculated from

[TABLE]

subject to $\boldsymbol{\beta}_{k}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\beta}_{k}=1$ and $[\boldsymbol{\beta}_{01},\cdots,\boldsymbol{\beta}_{0(k-1)}]^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\beta}_{k}=\mathbf{0}$. We continue this process till the MV index reaches [math]. This procedure is indeed conducting sufficient dimension reduction, which can be seen clearly from the following theorem.

Theorem 1**.**

Suppose there exists a positive integer $d<p$ such that $MV(\boldsymbol{\beta}_{01}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)\geq MV(\boldsymbol{\beta}_{02}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)\geq\dots\geq MV(\boldsymbol{\beta}_{0d}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)>0=MV(\boldsymbol{\beta}_{0(d+1)}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)=\dots=MV(\boldsymbol{\beta}_{0p}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)$ where $\boldsymbol{\beta}_{0i}^{{\mathrm{\scriptscriptstyle T}}}\boldsymbol{\beta}_{0i}=1$ and $\boldsymbol{\beta}_{0i}^{{\mathrm{\scriptscriptstyle T}}}\boldsymbol{\beta}_{0j}=0$ for $i,j=1,\dots,p$ and $i\neq j$. Then

[TABLE]

where $S_{Y|{\mathbf{X}}}$ denotes the central dimension reduction subspace (Cook (1994), Cook (1996)) and for any $k<d$,

[TABLE]

Note that the subscript [math] in $\boldsymbol{\beta}_{0i}$’s in Theorem 1 is used to indicate that these variables are in the population level. The existence of $d<p$ is validated in the classical linear discriminant analysis (LDA) and in the index model setting. See more details in Section 3.

In practice, we don’t know the population MV index for any given $\boldsymbol{\beta}$ and have to use $\widehat{MV}$ to replace it in the sequential optimization procedure. It is natural to estimate $F(z)$ in (2.1)by its empirical counterpart: $\hat{F}(z)=n^{-1}\sum_{i=1}^{n}\mathrm{I}(Z_{i}\leq z)$, as is done in Cui, Li and Zhong (2015). However, the empirical distribution is a step function which makes the optimization algorithm problematic. Instead, we use

[TABLE]

where $\hat{f}_{h}$ is a kernel density estimator of $f$ and $K_{h}(\cdot)=1/hK(\cdot/h)$ with $K$ being a kernel function and $h=h_{n}$ the bandwidth converging to 0 as $n\to\infty$. Similarly, $\hat{F}_{r}(z)$ can be estimated by

[TABLE]

for $r=1,2,\dots,R$, where $n_{r}$ is the sample size of the $r$th category and $Z_{j}$, $j=1,\dots,n_{r}$, are the sample points in this category. Then, we can use the estimator $\widehat{MV}$ given in $(\ref{mv})$ together with $\hat{F}_{h}(z)$ and $\hat{F}_{hr}(z)$ to implement the optimization algorithm. Note that $\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}}$ is of one dimension, $\widehat{MV}(\boldsymbol{\beta}):=\widehat{MV}(\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}}|Y)$ can be estimated by the approach given above.

3 MMV in classification

3.1 Fisher’s LDA

Consider a two-class classification problem. Suppose we have $n$ labeled i.i.d. training samples $(Y_{i},{\mathbf{X}}_{i}),1\leq i\leq n$, where ${\mathbf{X}}_{i}$ is a $p$-dimensional feature vector and $Y_{i}\in\{-1,1\}$ is the corresponding class label. Let $p_{1}=\mathbb{P}(Y_{i}=1)$, $p_{-1}=\mathbb{P}(Y_{i}=-1)$ and assume

[TABLE]

where $\boldsymbol{\mu}$ is the contrast mean vector between the two classes, and $\boldsymbol{\Sigma}$ is the $p\times p$ covariance matrix. Given a new independent feature vector from the same population, i.e. ${\mathbf{X}}\sim N(Y\cdot\boldsymbol{\mu},\,\boldsymbol{\Sigma})$, our goal is to train $({\mathbf{X}}_{i},Y_{i})$ to decide whether $Y=-1$ or $Y=1$. Here we use the contrast mean in model (3.1), but the method and result below also applies to a more general model with mean vectors $\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2}\in\mathbb{R}^{p}$ with no extra difficulty.

Linear Discriminant Analysis (LDA), namely Fisher’s LDA, is a well-known method for classification, which is essentially based on a weighted average of the test features $L({\mathbf{X}})=\sum_{j=1}^{p}w(j)X(j)$ and predicts $Y=\pm 1$ if $L({\mathbf{X}})><0$. Here, ${\mathbf{w}}=(w_{1},\ldots,w_{p})^{\mathrm{\scriptscriptstyle T}}$ is a preselected weight vector. Fisher showed that the optimal weight vector satisfies

[TABLE]

where $\boldsymbol{\Omega}=\boldsymbol{\Sigma}^{-1}$. In the classical setting where $n\gg p$, $\mu$ and $\Omega$ can be conveniently estimated and Fisher’s LDA is approachable. Unfortunately, in the modern regime where $p\gg n$, Fisher’s LDA faces immediate challenges. To bypass the difficulty of estimating $\Omega$ in LDA, we propose a classifier based on the transformed variables $\boldsymbol{\beta}_{0k}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}}$, $1\leq k\leq d$, where $\boldsymbol{\beta}_{0k}$s are the optimizers of the MMV procedure. The rationale of proposing this classifier is the intimate relationship given in the next theorem between the MMV method and the LDA.

Theorem 2**.**

Under model (3.1), $d=1$ and $\boldsymbol{\beta}_{01}\propto\boldsymbol{\Omega}\boldsymbol{\mu}$, where $\boldsymbol{\beta}_{01}$ is defined in (2.3).

In practice, if we have an estimator of $MV(Z|Y)$, say $\widehat{MV}(Z|Y)$ given in (2.2), then we can find the maximizer, denoted as $\hat{\boldsymbol{\beta}}_{1}$, of $\widehat{MV}(Z|Y)$. For a new given feature vector ${\mathbf{X}}$, we classify it as $Y=1$ if ${\hat{\boldsymbol{\beta}}_{1}}^{\top}{\mathbf{X}}>0$ and $Y=-1$ if ${\hat{\boldsymbol{\beta}}_{1}}^{\top}{\mathbf{X}}<0$. Since the $\hat{\boldsymbol{\beta}}_{1}$ is an estimator of $\boldsymbol{\beta}_{01}$, by Theorem 2 it is also an estimator of $\boldsymbol{\Omega}\boldsymbol{\mu}$. This classifier has the following benefits.

It does not have to estimate the precision matrix $\boldsymbol{\Omega}$ and can be implemented efficiently in the case where $n<p$. 2. 2.

It requires minimal conditions. Given $Y$, the transformed variables ${\boldsymbol{\beta}}_{01}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}},\cdots,{\boldsymbol{\beta}}_{0d}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}}$ have distributions closer to normal than the original variable ${\mathbf{X}}$ when $d\ll p$ (Hall and Li (1993)). It is well known that the conditional normality is an important assumption for linear discriminant analysis. 3. 3.

It has better performance in classification than the traditional LDA. We will show in our simulation study that it reduces the classification error significantly compared with the LDA.

This theorem states that for model (3.1), the true $d=1$ and $\boldsymbol{\beta}_{01}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}}$ contains all the information for classification. This means at the population level (assuming population functions are all known), the MMV procedure gives exactly the LDA classifier. The theorem also justifies the efficiency of the LDA for the normal model (3.1) in terms of maximum mean variance. When $p>n$, the LDA needs to estimate the inverse of the covariance matrix; thus it is unsolvable or requires extra sparsity assumption. The MMV is an efficient alternative to circumvent this problem.

In the field of high dimensional classification, we are usually faced with rare and weak signals, i.e. important features are sparse and each contributes weakly to the classification decision. In this setting, reasonable classification becomes quite difficult. If rare and weak signals imply sparse optimal weights, the proposed MMV estimator can be adjusted by adding a penalty term to the objective function or by some kinds of thresholding. We conjecture that higher criticism thresholding (HCT) might also be applied in this context. It has been verified that HCT performs quite well theoretically and practically when the signals are rare and weak; see Donoho and Jin (2008) and Fan, Jin and Yao (2013) for details. Hence, the extension of HCT to the MMV estimator seems rather direct and reasonable. Nevertheless, a theoretical investigation of this would be very complicated, and we leave it for further research.

3.2 Index model

The index model enjoys a lot of popularity in regression and classification. The logistic model and probit model are special cases where the link function is known with a single index. A general index model can be expressed as the following semi-parametric model. Let $Y\in\mathbb{R}$ denote the response variable and ${\mathbf{X}}\in\mathbb{R}^{p}$ denote the covariates. Assume there exist orthogonal $p$-dimensional vectors $\boldsymbol{\beta}_{1},\dots,\boldsymbol{\beta}_{k}$ with unit norm such that

[TABLE]

where $f$ is an arbitrary unknown link function and $\varepsilon$ is independent of ${\mathbf{X}}$. With a bit of an abuse of the notation, the notation $k$ in (3.2) can be seen as a fixed integer indicating the number of the indexes. The column space spanned by $\{\boldsymbol{\beta}_{1},\dots,\boldsymbol{\beta}_{k}\}$ is defined as the central dimension reduction subspace by Cook (1994) and Cook (1998). Under the setting of index model (3.2), we can detail Theorem 1 to some extent. Let ${\mathbf{B}}=(\boldsymbol{\beta}_{1},\dots,\boldsymbol{\beta}_{k})^{\mathrm{\scriptscriptstyle T}}$. Assume

[TABLE]

and there exists a $p$-dimensional vector $\boldsymbol{\gamma}$ such that

[TABLE]

then by Lemma 4.3 in Dawid (1979) and Proposition 4.6 in Cook (1998), we can get $MV(\boldsymbol{\gamma}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}}|Y)=0$. This implies that under mild conditions, the MMV method can recover all the information in ${\mathbf{X}}$ related to classification with $d<p$ indexes in Theorem 1. Specifically, when ${\mathbf{X}}\sim N(\boldsymbol{\mu},\boldsymbol{\Sigma})$, we can obtain the following theorem.

Theorem 3**.**

Assume ${\mathbf{X}}\sim N(\boldsymbol{\mu},\boldsymbol{\Sigma})$ and Y\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}{\mathbf{X}}|{\mathbf{B}}{\mathbf{X}} where ${\mathbf{B}}=(\boldsymbol{\beta}_{1},\dots,\boldsymbol{\beta}_{k})^{\mathrm{\scriptscriptstyle T}}$ with $k<p$. For $d$ defined in Theorem 1, if $2k\leq p$, then $d\leq 2k$. Specifically, if $\boldsymbol{\Sigma}={\mathbf{I}}$, then $d=k$.

Theorems 1 and 3 indicate that when $\boldsymbol{\Sigma}={\mathbf{I}}$, the MMV procedure can exactly recover the central subspace with $d=k$ steps under the setting of index model and normal covariates. To be specific, in the logistic (probit) model with normal covariates, $\boldsymbol{\beta}_{01}\propto\boldsymbol{\beta}_{1}$, where $\boldsymbol{\beta}_{1}$ denotes the coefficient vector of the logistic (probit) model and $\boldsymbol{\beta}_{01}$ is defined in Theorem 1. This implies that $d=1$ is enough for the logistic (probit) model. The advantage of our method is that it is a nonparametric method, thus we do not need to assume any specific link functions.

3.3 Other classification algorithms

Other popular classification methods such as K-Nearest Neighbours (KNN), neural networks and support vector machine (SVM), can be connected to MMV by a two step procedure, i.e. dimension reduction first and classification next. We conjecture that such a two step procedure will improve the accuracy of the classification because high dimensionality causes problems in the classification algorithms mentioned above. Simulations in Section 5 validate this conjecture.

4 Consistency and asymptotic normality

We introduce the following notations to simplify the description of the theory. $\{{\mathbf{X}}_{i},Y_{i}\}_{i=1}^{n}$ are i.i.d. samples. $n_{r}$ indicates the number of samples in the class $Y=y_{r}$ for $r=1,\dots,R$. For the case $d=2$ where $d$ is defined in Theorem 1, denote

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

where $\boldsymbol{\lambda}_{2}=(\lambda_{21},\lambda_{22})^{{\mathrm{\scriptscriptstyle T}}}$. Let $\boldsymbol{\theta}_{0i}=(\boldsymbol{\beta}_{0i}^{\mathrm{\scriptscriptstyle T}},\boldsymbol{\lambda}_{0i})^{\mathrm{\scriptscriptstyle T}}$ and $\hat{\boldsymbol{\theta}}_{i}=(\hat{\boldsymbol{\beta}}_{i}^{\mathrm{\scriptscriptstyle T}},\hat{\boldsymbol{\lambda}}_{i})^{\mathrm{\scriptscriptstyle T}}$ be the maximizers of $L^{(i)}(\boldsymbol{\theta}_{i})$ and $L_{nh}^{(i)}(\boldsymbol{\theta}_{i})$ respectively for $i=1,2$. Let $\Omega_{i}$ denote the parameter space of $\boldsymbol{\beta}_{i}$ and $B(\kappa_{i})=\{\boldsymbol{\beta}_{i}:\|\boldsymbol{\beta}_{i}-\boldsymbol{\beta}_{0i}\|\leq\kappa_{i}\}$ be a ball with center $\boldsymbol{\beta}_{0i}$ and radius $\kappa_{i}$ for $i=1,2$. The boundary of the ball is denoted by $\partial B(\kappa_{i})$. Denote $\Gamma_{i}=\{\boldsymbol{\beta}_{i}:\boldsymbol{\beta}_{i}^{{\mathrm{\scriptscriptstyle T}}}\boldsymbol{\beta}_{i}=1\}$ for $i=1,2$ and $\Upsilon=\{\boldsymbol{\beta}_{2}:\boldsymbol{\beta}_{01}^{{\mathrm{\scriptscriptstyle T}}}\boldsymbol{\beta}_{2}=0\}$. Let $C(\kappa_{i})=\{\boldsymbol{\beta}_{i}:\|\boldsymbol{\beta}_{i}-\boldsymbol{\beta}_{0i}\|\leq\kappa_{i}\}$ be a complex ball in $\mathbb{C}^{p}$ with center $\boldsymbol{\beta}_{0i}$ and radius $\kappa_{i}$ for $i=1,2$. Note that each element of $\boldsymbol{\beta}_{i}$ is complex. Denote $\Gamma_{Ci}=\{\boldsymbol{\beta}_{i}:\boldsymbol{\beta}_{i}^{{\mathrm{\scriptscriptstyle T}}}\boldsymbol{\beta}_{i}=1\}$ for $i=1,2$ and $\Upsilon_{C}=\{\boldsymbol{\beta}_{2}:\boldsymbol{\beta}_{01}^{{\mathrm{\scriptscriptstyle T}}}\boldsymbol{\beta}_{2}=0\}$. For $i=1,2$, let $\Omega_{Ci}$ be the parameter space of the complex $\boldsymbol{\beta}_{i}$.

The conditions below are required to establish the consistency and asymptotic normality.

$(1)$ $c_{1}/R\leq\min_{1\leq r\leq R}p_{r}\leq\max_{1\leq r\leq R}p_{r}\leq c_{2}/R$ and $R=O(n^{\delta})$ with $0<\delta\leq 1/2$.

$(2)$ a. There exists an open subset $\omega_{1}$ of $\Omega_{1}\cap\Gamma_{1}$ that contains the true parameter $\boldsymbol{\beta}_{01}$ for almost all $({\mathbf{X}},Y)$ and $\sup_{\boldsymbol{\beta}_{1}\in B(\kappa_{01})}MV(\boldsymbol{\beta}_{1}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)<\infty$, and for any $\kappa_{1}\in(0,\kappa_{01}]$

[TABLE]

b. There exists an open subset $\omega_{2}$ of $\Omega_{2}\cap\Gamma_{2}\cap\Upsilon$ that contains the true parameter $\boldsymbol{\beta}_{02}$ for almost all $({\mathbf{X}},Y)$ and $\sup_{\boldsymbol{\beta}_{2}\in B(\kappa_{02})}MV(\boldsymbol{\beta}_{2}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)<\infty$, and for any $\kappa_{2}\in(0,\kappa_{02}]$

[TABLE]

$(3)$ $\int uK(u)du=0$, $\int u^{2}K(u)du<\infty$ and $nh^{4}\to 0$ where $h=h_{1}=\dots=h_{R}$; $|K|_{\infty}=\sup_{u\in\mathbb{R}}|1/hK(u/h)|<\infty$.

$(4)$ There exists a constant $\kappa_{01}$ such that $MV(\boldsymbol{\beta}_{1}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)$ is an analytic function of each coordinate of $\boldsymbol{\beta}_{1}$ in $C(\kappa_{01})\subseteq\Omega_{C1}\cap\Gamma_{C1}$ and $\sup_{\boldsymbol{\beta}_{1}\in C(\kappa_{01})}MV(\boldsymbol{\beta}_{1}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)<\infty$. For any $\boldsymbol{\beta}_{1}\in C(\kappa_{01})$, $MV^{{}^{\prime}}(\boldsymbol{\beta}_{1})=MV^{{}^{\prime}}(\boldsymbol{\beta}_{1}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)$ and $MV^{{}^{\prime\prime}}(\boldsymbol{\beta}_{1})=MV^{{}^{\prime\prime}}(\boldsymbol{\beta}_{1}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)$ exist. $\sup_{\boldsymbol{\beta}_{1}\in C(\kappa_{01})}\|MV^{{}^{\prime}}(\boldsymbol{\beta}_{1})\|_{\infty}<\infty$ and $\sup_{\boldsymbol{\beta}_{1}\in C(\kappa_{01})}\|MV^{{}^{\prime\prime}}(\boldsymbol{\beta}_{1})\|_{\infty}<\infty$, where $\|MV^{{}^{\prime}}(\boldsymbol{\beta}_{1})\|_{\infty}=\inf\{C>0:\mathbb{P}(\|MV^{{}^{\prime}}(\boldsymbol{\beta}_{1})\|\leq C)=1\}$ and $\|MV^{{}^{\prime\prime}}(\boldsymbol{\beta}_{1})\|_{\infty}=\inf\{C>0:\mathbb{P}(\|MV^{{}^{\prime\prime}}(\boldsymbol{\beta}_{1})\|\leq C)=1\}$. (Note that we only give the condition for $i=1$. The conditions for the $i=2$ case are a bit more complex but quite similar to the $i=1$ case.)

$(5)$ $L^{{}^{\prime\prime}}(\boldsymbol{\theta}_{0})$ is nonsingular and $\rm E(\|{\mathbf{X}}/\sqrt{p}\|)<\infty$.

Condition $(1)$ requires the proportion of each response class to be moderate, not too small nor too large. $R=O(n^{\delta})$ allows the number of the classes to grow with the sample size. This condition is also imposed by Cui, Li and Zhong (2015). Condition $(2)$ is assumed to ensure the existence of the MMV optimizers. A similar condition is assumed in Chen, Ma and Zhou (2017) for a likelihood function. Simulation results below also verify this assumption. Condition $(3)$ is widely used in the literature of kernel density estimators and is assumed for the consistency of the estimator; see Cheng (2017) for reference. Condition $(4)$ is not as strict as it seems. Recall that the MV index is defined on the cumulative distribution functions which are, of course, bounded. Therefore, we only need the bounded constraints for the corresponding density function and the first derivative of the density function, which are both quite mild. Condition $(5)$ is required to guarantee the root-$n$ consistency of the proposed estimator; it is in the spirit of the Von Mises proposition (Serfling (1980), Section $6.1$).

Theorem 4** (Consistency).**

Let $\{{\mathbf{X}}_{i},Y_{i}\}_{i=1}^{n}$ be i.i.d. samples. For both fixed and diverging $R$, under Conditions $(1)$ – $(3)$, we have $\widehat{\boldsymbol{\beta}}_{i}\to\boldsymbol{\beta}_{0i}$ in probability as $n\to\infty$ for $i=1,2$. Moreover, when $p$ satisfies $p^{p/2}n^{-\alpha(1-\delta)}=o(1)$ for any $\alpha<1/2$, under Conditions $(1)$ – $(3)$, $\widehat{\boldsymbol{\beta}}_{i}\to\boldsymbol{\beta}_{0i}$ in probability as $n\to\infty$ for $i=1,2$.

Theorem 4 indicates that the MMV estimators are consistent for both fixed and diverging $p$ and $R$ cases under regular conditions. Here, we only present the theorem for the case $d=2$ for the convenience of statement. When $p$ is fixed, the $\sqrt{n}$ consistency and asymptotic normality can be further proved.

Theorem 5** (Asymptotic normality).**

Let $\{{\mathbf{X}}_{i},Y_{i}\}_{i=1}^{n}$ be i.i.d. samples. Under Conditions $(1)$ – $(5)$, $\sqrt{n}(\widehat{\boldsymbol{\beta}}_{i}-\boldsymbol{\beta}_{0i})$ is asymptotically normal with mean zero and covariance matrix $V_{i}$, where $V_{i}$ is defined in the proof for $i=1,2$.

When $p$ is diverging with $n$, the asymptotic normality given above may or may not hold, which depends on the property of the MV index, thus the distribution of ${\mathbf{X}}$ and ${\mathbf{X}}|Y$. When the uniform convergence stays for the first and second derivatives of the estimated MV index with a properly fast convergence rate, the convergence rate of the MMV estimator can be proved via techniques similar to those employed in Fan and Peng (2004) for the non-concave penalized likelihood when the number of parameters is diverging with the sample size. Consider the simple case $d=1$. More conditions are needed (the subscript “1” of $\boldsymbol{\beta}_{01}$, $\boldsymbol{\theta}_{1}$ and $\boldsymbol{\theta}_{01}$ and the superscript “(1)” of $L^{(1)}$ and $L^{(1)}_{nh}$ are omitted for simplicity):

$(6)$ $\|\nabla MV_{n}(\boldsymbol{\beta}_{0})-\nabla MV(\boldsymbol{\beta}_{0})\|=O_{p}(\alpha_{n})$ where $\alpha_{n}=\alpha_{n}(n,p,h)$ and $p^{3/2}\alpha_{n}=o_{p}(1)$.

$(7)$ $\lambda_{\max}\{\nabla^{2}MV_{n}(\boldsymbol{\beta}_{0})-\nabla^{2}MV(\boldsymbol{\beta}_{0})\}=o_{p}(1)$.

$(8)$ There is a large enough ball centering at $\boldsymbol{\theta}_{0}$ such that for any $\boldsymbol{\theta}$ in the ball,

[TABLE]

for $1\leq j,k,l\leq p$ where $M_{1}$ and $M_{2}$ are positive numbers.

Theorem 6**.**

If Conditions $(1)$-$(8)$ hold, there exists a local maximizer $\hat{\boldsymbol{\beta}}$ such that $\|\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}_{0}\|=O_{p}(\alpha_{n})$ where $\alpha_{n}$ is defined in Condition $(6)$.

5 Simulations

Although the MMV method can be readily used under settings of $n<p$, the computation cost for a sequential algorithm like ours is quite high. It is usually better to use cross validation to choose the proper number of indexes $d$ and empirical bandwidth $h$ in practice. However, the computation cost of doing so is also very high. For simplicity, instead, in the following simulation studies, we use $d$ as the true dimension of the central subspace and $h=3\cdot sd(\tilde{\boldsymbol{\beta}}_{1})n^{-1/3}$ where $sd$ stands for standard deviation and $\tilde{\boldsymbol{\beta}}_{1}$ is a good initial estimate of $\boldsymbol{\beta}_{1}$. If the predictors’ dimension is ultra-high, our suggestion is to conduct feature screening first to reduce dimensionality $p$ (say, $\exp(O(n^{\xi}))$ for some $\xi>0$) to a relatively large scale $d^{\prime}$ (e.g. , o(n)) by fast methods such as those of Fan and Fan (2008) and Cui, Li and Zhong (2015). When the size of $p$ is comparable to $n$, our method is quite fast and effective. We use ten-fold cross validation to calculate the classification error in both simulation and real data analysis. We repeat the experiment 400 times and the average classification error and the corresponding standard deviation (in parentheses) are then calculated. Let $\boldsymbol{\beta}_{1}=(1,1,1,1,0,\ldots,0)^{\mathrm{\scriptscriptstyle T}}$ and $\boldsymbol{\beta}_{2}=(1,-1,1,-1,0,\ldots,0)^{{\mathrm{\scriptscriptstyle T}}}$. The calculation for LDA, logistics regression, SVM and KKN is based on the corresponding Matlab (R2015a) packages using default settings.

5.1 Fisher’s linear discriminant analysis

In this study, we set $p$ be 50 and 200 respectively, with the sample size $n=80$. We generate ${\mathbf{Y}}=(1,\ldots,1,-1,\ldots,-1)^{\sf T}$ first, then generate ${\mathbf{X}}$ as follows. It is an ordinary LDA model which is in fact an inverse model with a one-dimensional central subspace.

Model I

[TABLE]

where $\boldsymbol{\epsilon}\sim N({\mathbf{0}},{\bf I}_{n})$ and $\boldsymbol{\Delta}_{ij}=0.5^{|i-j|}$ for $1\leq i,j\leq 10$.

Table 1 shows that MMV+LDA outperforms LDA significantly in both settings ($p=50$ and $p=200$). By applying MMV, the classification error is decreased by 50 percent or so. Although MMV+LDA equals LDA in the population level, the former does work better in the finite sample settings. The reason is that dimension reduction through MMV increases estimation efficiency.

5.2 Logistic regression

In this study, we set $p$ be 20 and 50 respectively, with the sample size $n=80$. Since logistic regression utilizes likelihood estimation, the sample size is required to be larger than the dimension of the predictors. We generate data using the logistic model as follows.

Model II

[TABLE]

where $\mathrm{I}$ is the indictor function and ${\mathbf{X}}=(X_{1},\ldots,X_{p})^{\ T}\sim N({\bf 0},\boldsymbol{\Psi})$ with $\boldsymbol{\Psi}_{ij}=0.5^{|i-j|}$ for $1\leq i,j\leq p$. In this study, the central subspace is spanned by the direction $\boldsymbol{\beta}_{1}$.

It can been seen clearly from Table 2 that MMV, as a dimension reduction technique, improves estimation efficiency, and thus reduces classification error remarkably when it is combined with Logistic regression.

5.3 More complex models

We compare our method with more advanced algorithms like SVM and KKN in this study. Models III and IV are multiple index models. We set $p$ be 50 and 200 respectively, while the sample size $n=160$. In these two models, the central subspace is spanned by the directions $\boldsymbol{\beta}_{1}$ and $\boldsymbol{\beta}_{2}$.

Model III

[TABLE]

where $\epsilon\sim N(0,1)$, ${\mathbf{X}}=(X_{1},\ldots,X_{p})^{{\mathrm{\scriptscriptstyle T}}}\sim N({\bf 0},\boldsymbol{\Psi})$ with $\boldsymbol{\Psi}_{ij}=0.5^{|i-j|}$ for $1\leq i,j\leq p$, and {\mathbf{X}}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}\epsilon.

Model IV

[TABLE]

where $\epsilon\sim N(0,1)$, ${\mathbf{X}}=(X_{1},\ldots,X_{p})^{\sf T}\sim N({\bf 0},\boldsymbol{\Psi})$ with $\boldsymbol{\Psi}_{ij}=0.5^{|i-j|}$ for $1\leq i,j\leq p$, and {\mathbf{X}}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}\epsilon.

Table 3 indicates that even for complex classification techniques such as SVM and KNN, employing MMV before classification still enjoys a significant decrease of the classification error. This further conforms the efficiency and priority of our proposed method.

5.4 Real data analysis

We apply our method to human colon cancer data with $n=62$ and $p=2000$, which is available in R. There are 40 samples from tumors “t”, and 22 samples are from normal “n” biopsies. The data was originally collected on microarrays with 6500 probes. 2000 of them were selected apparently, randomly, to be used for demonstrating statistical methods. We first screen the number of the predictors to 100 by the method of Cui, Li and Zhong (2015), and compare our methods with LDA, SVM and KKN. We apply the same bandwidth selection and cross validation methods as those in the simulation study. A few choices of the dimension $d=1,2,3$ are tried and they have rather similar results. Here we only present the result with fixed $d=1$. To get a fair comparison, we repeat the permutation 100 times for cross validation results. The table below summarizes the average classification errors and the corresponding standard deviations.

Table 4 demonstrates that “MMV+.” performs much better than the original classification method. It seems that the performances of the three methods are comparable to each other, with SVM performing a little better, while MMV+LDA performs the best among the three “MMV+.” methods. We conjecture that the relationship between MMV and the original classification techniques may be more than simple addition.

6 Conclusion

In this paper, we propose a new sufficient dimension reduction approach - Maximal Mean Variance - which is designed deliberately for high-dimensional classification. Our method requires fairly mild restrictions on the predicting variables and keeps the model-free advantage without the need to estimate the link function; thus it can be widely used. The relationships between MMV and popular classification methods are discussed in detail. The consistency and convergence rate of the MMV estimator are established for both fixed and diverging dimension ($p$) cases and the number of the response classes can also be allowed to grow with the sample size $n$. When the dimension of the predicting vector is fixed, we further construct the asymptotic normality for the MMV estimator. Simulations and real data analysis validate the efficiency and priority of the proposed method. While cross validation can be employed to choose a proper $d$ in practice, it would be quite challenging to derive an optimal $d$ theoretically under some kind of criterion because the algorithm of MMV is a stepwise procedure. We leave it for further research. Besides, the MMV method can be readily applied to the ultra high dimensional setting by conducting a screening procedure first. This two-scale learning framework is in the sprit of Fan and Lv (2008) for sure independence screening.

7 Appendix

7.1 Proof of Theorem 1

Proof.

By the property of MV index, we have \boldsymbol{\beta}_{0i}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}\not\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y for $i=1,\dots,d$ and \boldsymbol{\beta}_{0i}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y for $i=d+1,\dots,p$. Recalling the definition of central subspace, this argument means $\boldsymbol{\beta}_{0i}^{{\mathrm{\scriptscriptstyle T}}}\notin S_{Y|{\mathbf{X}}}^{\perp}$ for $i=1,\dots,d$ and $\boldsymbol{\beta}_{0i}^{{\mathrm{\scriptscriptstyle T}}}\in S_{Y|{\mathbf{X}}}^{\perp}$ for $i=d+1,\dots,p$ where $S_{Y|{\mathbf{X}}}^{\perp}$ is the orthogonal complement space of $S_{Y|{\mathbf{X}}}$. Noting that $\boldsymbol{\beta}_{0i}^{{\mathrm{\scriptscriptstyle T}}}\boldsymbol{\beta}_{0i}=1$, $\boldsymbol{\beta}_{0i}^{{\mathrm{\scriptscriptstyle T}}}\boldsymbol{\beta}_{j0}=0$ for $i,j=1,\dots,p$ and $i\neq j$ and ${\rm span}\{\boldsymbol{\beta}_{(d+1)0},\dots,\boldsymbol{\beta}_{p0}\}\subseteq S_{Y|{\mathbf{X}}}^{\perp}$, we easily obtain $S_{Y|{\mathbf{X}}}\subseteq{\rm span}\{\boldsymbol{\beta}_{01},\dots,\boldsymbol{\beta}_{0d}\}$.

For any $k<d$, if ${\rm span}\{\boldsymbol{\beta}_{01},\dots,\boldsymbol{\beta}_{0k}\}\supseteq S_{Y|{\mathbf{X}}}$, we can obtain for any $i>k$, \boldsymbol{\beta}_{i}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y which leads to a contradiction at least for $i=d>k$.

∎

7.2 Proof of Theorem 2

Proof.

Let $Z$ denote the variable $\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}}$ and $Z_{i}$ denote the variable $\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}}|Y=i$ ($i=1,-1$) for $\boldsymbol{\beta}\in\mathbb{R}^{p}$ satisfying $\|\boldsymbol{\beta}\|=1$. Then, $Z_{1}\sim N(\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\mu},\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}\boldsymbol{\beta})$, $Z_{-1}\sim N(-\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\mu},\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}\boldsymbol{\beta})$ and

[TABLE]

where $F_{i}(Z)$ denotes the conditional distribution function of $Z|Y=i$ for $i=1,-1$. Then,

[TABLE]

Hence, we only need to maximize $\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\mu}/\sqrt{\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}\boldsymbol{\beta}}$. The solution is exactly the optimal weight of LDA, i.e. $\boldsymbol{\beta}_{01}\propto\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}$, where $\boldsymbol{\beta}_{01}$ is defined in Theorem 1. Thus, if we choose $d=1$ in the MMV procedure, then MMV+LDA equals LDA at the population level.

Then we try to find $\boldsymbol{\beta}_{02}$. We always assume $\boldsymbol{\Sigma}$ is positive definite, and thus for any non-zero $\boldsymbol{\beta}$, $\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}\boldsymbol{\beta}>0$. Since now we maximize $(\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}\boldsymbol{\beta})^{-1/2}\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\mu}$ subject to $\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\beta}=1$ and $\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}=0$, let

[TABLE]

Take first partial derivative w.r.t. $\boldsymbol{\beta}$ and set it to zero we have

[TABLE]

Then

[TABLE]

i.e. $\lambda=0$. Now $(\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}\boldsymbol{\beta})^{3/2}\frac{\partial L}{\partial\boldsymbol{\beta}}=0$ is reduced to

[TABLE]

Multiply both sides of above equation by $\boldsymbol{\Sigma}^{-1}\boldsymbol{\beta}$ from right, together with $\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\beta}=1$ and $\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}=0$, we have

[TABLE]

and thus

[TABLE]

Note that none of $\pi$, $\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}^{-2}\boldsymbol{\beta}$ and $\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\beta}$ can be 0. If $\pi=0$, then (7.2) gives $\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\beta}=0$. Plugging this and $\pi=0$ back into (7.1) gives $\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}\boldsymbol{\beta}\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}=0$, i.e. $\boldsymbol{\mu}=0$, which contradicts the assumption $\boldsymbol{\mu}\neq 0$ and thus $\pi\neq 0$. If $\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}^{-2}\boldsymbol{\beta}=0$, then (7.2) gives $\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\beta}=0$. Multiplying $\boldsymbol{\Sigma}^{-2}\boldsymbol{\beta}$ from the right to (7.1) and plugging in $\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}^{-2}\boldsymbol{\beta}=0$ and $\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\beta}=0$, we have $\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}^{-3}\boldsymbol{\beta}=0$ as $\pi\neq 0$. If we repeat this deduction then we have $\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}^{-i}\boldsymbol{\beta}=0$, $i=0,1,\cdots$. This can not hold for general $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ unless $\boldsymbol{\beta}=0$, which contradicts $\boldsymbol{\beta}=1$, and thus $\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}^{-2}\boldsymbol{\beta}\neq 0$. If $\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\beta}=0$, then $\pi=0$ which contradicts $\pi\neq 0$. Therefore (7.3) holds, is well-defined and $\pi\neq 0$. Now plug the $\pi$ back into (7.1) and we have

[TABLE]

and thus $\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}U=\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}U^{\mathrm{\scriptscriptstyle T}}-W$. By multiplying this equation from the right by $\boldsymbol{\beta}$ and applying it again, we have

[TABLE]

i.e. $0=W\boldsymbol{\beta}+\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}W^{\mathrm{\scriptscriptstyle T}}=\left(\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}^{-2}\boldsymbol{\beta}\right)^{-1}\left(\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\beta}\right)\left[\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}^{-1}\boldsymbol{\beta}+\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}\right]$. Since $\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\beta}\neq 0$, we have $\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}^{-1}\boldsymbol{\beta}+\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}=0$. Its solution is of the form $\boldsymbol{\beta}=\boldsymbol{\Sigma}(V-V^{\mathrm{\scriptscriptstyle T}})\boldsymbol{\mu}$ with $V$ a $p\times p$ matrix. Plugging this $\boldsymbol{\beta}$ back into (7.4), we have

[TABLE]

By the arbitrariness of $\boldsymbol{\mu}$, the above equation gives

[TABLE]

Multiplying it by $\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}^{-2}$ from the left, we have

[TABLE]

from which, by the arbitrariness of $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ again, we have

[TABLE]

Thus

[TABLE]

If we take transpose it follows that

[TABLE]

Combining the above two equations gives $\boldsymbol{\Sigma}^{-1}(V^{\mathrm{\scriptscriptstyle T}}-V)\boldsymbol{\Sigma}S=S\boldsymbol{\Sigma}(V^{\mathrm{\scriptscriptstyle T}}-V)\boldsymbol{\Sigma}^{-1}$ which has the solutions $S=\boldsymbol{\Sigma}^{-2}$ and $S=0$. If $S=\boldsymbol{\Sigma}^{-2}$, then by the fact $S^{\mathrm{\scriptscriptstyle T}}=-S$ we have $\boldsymbol{\Sigma}^{-2}=-\boldsymbol{\Sigma}^{-2}$ which is not valid. Thus we have $S=0$. From it we have $\boldsymbol{\mu}\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}(V^{\mathrm{\scriptscriptstyle T}}-V)\boldsymbol{\Sigma}^{2}=\left[\boldsymbol{\mu}\boldsymbol{\mu}^{\mathrm{\scriptscriptstyle T}}(V^{\mathrm{\scriptscriptstyle T}}-V)\boldsymbol{\Sigma}^{2}\right]^{\mathrm{\scriptscriptstyle T}}$ which has solutions $V^{\mathrm{\scriptscriptstyle T}}-V=\boldsymbol{\Sigma}^{-2}$ and $V^{\mathrm{\scriptscriptstyle T}}-V=0$. If $V^{\mathrm{\scriptscriptstyle T}}-V=\boldsymbol{\Sigma}^{-2}$, then $(V^{\mathrm{\scriptscriptstyle T}}-V)^{\mathrm{\scriptscriptstyle T}}=(\boldsymbol{\Sigma}^{-2})^{\mathrm{\scriptscriptstyle T}}$, i.e. $V-V^{\mathrm{\scriptscriptstyle T}}=\boldsymbol{\Sigma}^{-2}$, i.e. $-\boldsymbol{\Sigma}^{-2}=\boldsymbol{\Sigma}^{-2}$ a contradiction. Thus $V^{\mathrm{\scriptscriptstyle T}}-V=0$ which gives $\boldsymbol{\beta}=\boldsymbol{\Sigma}(V-V^{\mathrm{\scriptscriptstyle T}})\boldsymbol{\mu}=0$.

Therefore $\boldsymbol{\beta}_{02}=0$ for general $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$. For some specific $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$, there might exist nonzero $\boldsymbol{\beta}_{02}$, but always $\boldsymbol{\beta}_{02}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\mu}=0$. This gives the conclusion that for multivariate normal, the true $d=1$ and only the first $\boldsymbol{\beta}_{01}$ contributes to the variation among different classes.

∎

7.3 Proof of Theorem 3

Proof.

From Y\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}{\mathbf{X}}|{\mathbf{B}}{\mathbf{X}} and ${\mathbf{B}}\boldsymbol{\gamma}=\mathbf{0}$, we get Y\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}\boldsymbol{\gamma}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}}|{\mathbf{B}}{\mathbf{X}} (Proposition 4.3 in Cook (1998)). If we also have \boldsymbol{\gamma}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}{\mathbf{B}}{\mathbf{X}}, then by Lemma 4.3 of Dawid (1979) and Proposition 4.6 of Cook (1998), we obtain \boldsymbol{\gamma}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y which implies $MV(\boldsymbol{\gamma}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}}|Y)=0$ according to the property of the MV index. When ${\mathbf{X}}\sim N(\boldsymbol{\mu},\boldsymbol{\Sigma})$, then \boldsymbol{\gamma}^{\mathrm{\scriptscriptstyle T}}{\mathbf{X}}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}{\mathbf{B}}{\mathbf{X}} if and only if

[TABLE]

If $\boldsymbol{\Sigma}=\mathbf{I}$, then (7.5) holds by ${\mathbf{B}}\boldsymbol{\gamma}=\mathbf{0}$. Let $\boldsymbol{\Gamma}$ denote the vector space spanned by all $\boldsymbol{\gamma}$ satisfying (7.5) and ${\mathbf{B}}\boldsymbol{\gamma}=\mathbf{0}$, then we can easily get $\mathrm{dim}(\boldsymbol{\Gamma})=p-k$ which implies $d=k$ for $d$ defined in Theorem 1. If $\boldsymbol{\Sigma}\neq\mathbf{I}$ and $2k\leq p$, $\mathrm{dim}(\boldsymbol{\Gamma})\geq p-2k$ (note that $\boldsymbol{\Sigma}$ is positive definite), thus $d\leq 2k$. ∎

7.4 Proof of Theorem 4

Lemma 1**.**

*Under Condition $(1)$, for any $\boldsymbol{\beta}\in\mathbb{R}^{p}$, we have

$(1)$ if $R$ is fixed, $MV_{n}(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)\to MV(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)$ in probability as $n\to\infty$.

$(2)$ if $R$ is diverging with $n$ and satisfies $R=O(n^{\delta})$ with $\delta\leq 1/2$, $MV_{n}(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)\to MV(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)$ in probability as $n\to\infty$.*

Proof.

Denote $\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}$ by $X$ with support $\mathbb{R}_{X}$ and the transformed samples $\{\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}_{j}\}_{j=1}^{n}$ by $\{X_{j}\}_{j=1}^{n}$. By the definitions of $MV(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)$ and $MV_{n}(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)$, we have

[TABLE]

For the first term $A_{1}$,

[TABLE]

where the second inequality is obtained by $\int d\hat{F}(x)=1$. We then consider the term $B_{1}$,

[TABLE]

where the second equality is implied by Theorem 2.2 of Cheng (2017) with any $\alpha\leq 1/2$ under mild conditions and the last equality is given by Condition $(1)$ and the fact that $|\hat{p}_{r}-p_{r}|=O_{p}(n^{-1/2})$. For the second term, also by Theorem 2.2 of Cheng (2017), we obtain $B_{2}=\sup_{x\in\mathbb{R}_{X}}|\hat{F}_{h}(x)-F(x)|=O_{p}(n^{-\alpha})$.

We turn to the term $A_{2}$,

[TABLE]

where the last equality is based on the extended Glivenko-Cantelli lemma (in Fabian (1985)) with $\alpha$ defined above.

For the last term $A_{3}$, we have $|A_{3}|=O_{p}(n^{-\alpha})$ with any $\alpha\leq 1/2$ by Lemma A.4 of Cui, Li and Zhong (2015). To sum up, $|MV_{n}(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)-MV(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)|=|A_{1}+A_{2}+A_{3}|\leq O_{p}(n^{-\alpha(1-\delta)})+O_{p}(n^{-\alpha})=O_{p}(n^{-\alpha(1-\delta)})$. Thus, we complete the proof. ∎

Here, we give a proof for the $d=2$ case of Theorem 4.

Proof.

We first prove $\widehat{\boldsymbol{\beta}}_{1}\to_{p}\boldsymbol{\beta}_{01}$.

From Lemma 1, we obtain that for any $\boldsymbol{\beta}_{1}\in B(\kappa_{1})$ with $\kappa_{1}\in(0,\kappa_{01}]$, $MV_{n}(\boldsymbol{\beta}_{1}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)\to_{p}MV(\boldsymbol{\beta}_{1}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)$. For any $\epsilon>0$, let $\{\boldsymbol{\beta}^{1},\dots,\boldsymbol{\beta}^{M}\}$ be an $\epsilon/\sqrt{p}$-net of $B(\kappa_{01})$ with $M=(2\kappa_{01}\sqrt{p}/\epsilon+1)^{p}$. Since $M$ is fixed, by Lemma 1 we obtain

[TABLE]

as $n\to\infty$. For any $\boldsymbol{\beta}\in B(\kappa_{01})$, there exists a $m\in\{1,2,\dots,M\}$ such that $\|\boldsymbol{\beta}-\boldsymbol{\beta}^{m}\|\leq\epsilon/\sqrt{p}$. Then by Conditions $(4)$ and $(5)$, $|MV_{n}(\boldsymbol{\beta})-MV_{n}(\boldsymbol{\beta}^{m})|\leq 1/n\sum_{i=1}^{n}|K_{h}|_{\infty}\|{\mathbf{X}}_{i}/\sqrt{p}\|\epsilon\to_{p}|K_{h}|_{\infty}\rm E\|{\mathbf{X}}/\sqrt{p}\|\epsilon$ and $|MV(\boldsymbol{\beta})-MV(\boldsymbol{\beta}^{m})|\to_{p}0$. Therefore, combining the above three equations, we obtain

[TABLE]

Then by Condition $(2.a)$, it is easy to see that for any sufficiently small $\kappa_{1}>0$,

[TABLE]

as $n\to\infty$. Thus, there exists a local maximum local point $\widehat{\boldsymbol{\beta}}_{1}\in\partial B(\kappa_{1})\cap\Gamma_{1}$ with probability approaching to $1$ which means $\mathbb{P}(\|\widehat{\boldsymbol{\beta}}_{1}-\boldsymbol{\beta}_{01}\|<\kappa_{1})\to 1$.

Next, Let $p$ be diverging with $n$ and $p^{p/2}n^{-\alpha(1-\delta)}=o(1)$. Note that Lemma 1 still holds when $p$ is diverging, because the dimension of $\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}$ stays to be $1$ whether the dimension of ${\mathbf{X}}$ diverges or not. Recall that $|MV_{n}(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)-MV(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)|=O_{p}(n^{-\alpha(1-\delta)})$ (See Lemma 1). Then

[TABLE]

By the same arguments for the fixed $p$ case, we can complete the first part of the proof.

Noting that $\widehat{\boldsymbol{\beta}}_{1}\to_{p}\boldsymbol{\beta}_{01}$ and Condition $(2.b)$, we can easily obtain $\widehat{\boldsymbol{\beta}}_{2}\to_{p}\boldsymbol{\beta}_{02}$ by the same arguments given above.

∎

7.5 Proof of Theorem 5

Lemma 2**.**

*Under Condition $(1)$, for any $\boldsymbol{\beta}\in\mathbb{C}^{p}$, we have

$(1)$ if $R$ is fixed, $MV_{n}(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)\to MV(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)$ in probability as $n\to\infty$.

$(2)$ if $R$ is diverging with $n$ and satisfies $R=O(n^{\delta})$ with $\delta\leq 1/2$, $MV_{n}(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)\to MV(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)$ in probability as $n\to\infty$.*

Proof.

The proof is similar to that of Lemma 1. The only difference comes from the complex $\boldsymbol{\beta}$ and the corresponding complex random variable $X=:\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}$ with support $\mathbb{C}_{X}$. Denote $X=:a+ib$ where ${\mathbf{Z}}:=(a,b)^{{\mathrm{\scriptscriptstyle T}}}$ is a real random vector. By the definition of the cumulative distribution function of a complex variable, that is, the joint distribution function of the real part and the imaginary part of the variable, under several mild conditions we obtain

[TABLE]

This convergence rate comes from Theorem 3 of Liu and Yang (2008). Then, by the same arguments of Lemma 1, we complete the proof. ∎

The following is the proof of Theorem 5.

Proof.

We first present the proof for $i=1$. For simplicity, we omit $i=1,2$ in the subscripts and superscripts. Denote $\boldsymbol{\theta}=(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}},\lambda)^{{\mathrm{\scriptscriptstyle T}}}\in\mathbb{R}^{p+1}$. Let $\widehat{\boldsymbol{\theta}}=(\widehat{\boldsymbol{\beta}}^{\mathrm{\scriptscriptstyle T}},\hat{\lambda})^{\mathrm{\scriptscriptstyle T}}$ be the maximizer of $L_{nh}(\boldsymbol{\theta})$, then $\widehat{\boldsymbol{\theta}}=(\widehat{\boldsymbol{\beta}}^{\mathrm{\scriptscriptstyle T}},\hat{\lambda})^{\mathrm{\scriptscriptstyle T}}$ is a stationary point of $L_{nh}(\boldsymbol{\theta})$, that is, $L_{nh}^{{}^{\prime}}(\widehat{\boldsymbol{\theta}})=0$. Similarly, let $\boldsymbol{\theta}_{0}=(\boldsymbol{\beta}_{0}^{\mathrm{\scriptscriptstyle T}},\lambda_{0})^{\mathrm{\scriptscriptstyle T}}$ be the maximizer of $L(\boldsymbol{\theta})$, then $\boldsymbol{\theta}_{0}=(\boldsymbol{\beta}_{0}^{\mathrm{\scriptscriptstyle T}},\lambda_{0})^{\mathrm{\scriptscriptstyle T}}$ is a stationary point of $L(\boldsymbol{\theta})$. Denote $MV(\boldsymbol{\beta})=MV(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)$ and $MV_{n}(\boldsymbol{\beta})=MV_{n}(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{X}}|Y)$ for the simplicity of notation.

We first prove $n^{1/2}(\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0})\to_{d}N(0,V)$ where the covariance matrix $V$ will be given in the proof below. By the Taylor expansion, we have

[TABLE]

where $\boldsymbol{\theta}^{*}$ satisfies $\|\boldsymbol{\theta}^{*}-\boldsymbol{\theta}_{0}\|\leq\|\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\|$ with $\|\cdot\|$ being the Frobenius norm and $\boldsymbol{\theta}^{*}=(\boldsymbol{\beta}^{*{\mathrm{\scriptscriptstyle T}}},\lambda^{*})^{{\mathrm{\scriptscriptstyle T}}}$. With regular calculation, we obtain

[TABLE]

and

[TABLE]

where $I_{p}$ denotes the identity matrix of dimension $p\times p$. The remainder term $R(\boldsymbol{\theta}^{*})$ contains the third derivative of $L_{nh}(\boldsymbol{\theta})$ at $\boldsymbol{\theta}=\boldsymbol{\theta}^{*}$. Let $T_{n}=L_{nh}^{{}^{\prime\prime\prime}}(\boldsymbol{\theta}^{*})$, where $T_{n}$ is an array of dimension $(p+1)\times(p+1)\times(p+1)$ and for each $j=1,\dots,(p+1)$ and $T_{n}(j,;,;)$ is a matrix of dimension $(p+1)\times(p+1)$. Hence, we can write

[TABLE]

Then, based on the explicit expressions of the derivatives given above and (7.6), we obtain

[TABLE]

Next, our proof is divided into two parts:

Part 1:

[TABLE]

Part 2:

[TABLE]

where $=_{d}$ indicates convergence in distribution.

Proof of Part 1. If $\boldsymbol{\theta}=(\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}},\lambda)^{{\mathrm{\scriptscriptstyle T}}}$ is a complex vector, the above arguments still stand. By Condition $(4)$ and Lemma 2, it follows from the convergence of analytic functions that

[TABLE]

Denote $\boldsymbol{\beta}=(b_{1},\dots,b_{p})^{{\mathrm{\scriptscriptstyle T}}}$ and $T_{jn}=\sqrt{n}\partial{MV_{n}(\boldsymbol{\beta})}/\partial{b_{j}}|_{\boldsymbol{\beta}=\boldsymbol{\beta}_{0}}$. Note that $\boldsymbol{\beta}_{0}^{T}\boldsymbol{\beta}_{0}-1=0$, then it suffices to prove the asymptotic normality of $T_{jn}$ for each $j=1,\dots,p$. It is worth noting the elements of $-A\sqrt{n}L_{nh}^{{}^{\prime}}(\boldsymbol{\theta}_{0})$ are linear combinations of $T_{jn}$s. We consider the $j=1$ case.

By the definition of $MV_{n}(\boldsymbol{\beta})$, we have

[TABLE]

By Lemma 2,

[TABLE]

where $\boldsymbol{\beta}\in C(\kappa_{0})$, $i^{2}=-1$, and $u_{n}(\boldsymbol{\beta})+iv_{n}(\boldsymbol{\beta})=o_{p}(1)$ is uniform in $\boldsymbol{\beta}\in C(\kappa_{0})$ when $n\to\infty$ with $u_{n}(\boldsymbol{\beta})$ and $v_{n}(\boldsymbol{\beta})$ being real functions of $\boldsymbol{\beta}$. By Cauchy’s residue theorem, we have

[TABLE]

where $\tilde{\boldsymbol{\beta}}=(b_{1},b_{02},\dots,b_{0p})^{{\mathrm{\scriptscriptstyle T}}}$ with $\boldsymbol{\beta}_{0}=(b_{01},\dots,b_{0p})^{\mathrm{\scriptscriptstyle T}}$, and $C_{1}$ satisfies $\{b_{1}\in\mathbb{C}:\|b_{1}-b_{01}\|=r\}$ with $r<\kappa_{0}$. Define

[TABLE]

Then,

[TABLE]

Let

[TABLE]

Noting that the left-hand side of (7.8) is real, we consider the real part of the other hand, that is

[TABLE]

where the arguments of $R_{n},I_{n},u_{n},v_{n}$ are $b_{01}+re^{i\mu}$. By the mean value theorem, we obtain

[TABLE]

where the arguments of $R_{0n},I_{0n},u_{0n},v_{0n}$ are all $b_{01}+re^{i\mu_{0}}$ with $\mu_{0}\in[0,2\pi]$. By the definition of $S_{1n}$, we have

[TABLE]

By central limit theorem, $S_{1n}$ is asymptotically normally distributed. Then, noticing $u_{n}(\boldsymbol{\beta})+iv_{n}(\boldsymbol{\beta})=o_{p}(1)$ for $\boldsymbol{\beta}\in C(\kappa_{0})$ and letting $r\to 0$, $|T_{1n}-S_{1n}|\leq 2/r(|R_{0n}|+|I_{0n}|)(|u_{0n}|+|v_{0n}|)=o_{p}(1)$ as $n\to\infty$. Using Slutsky’s theorem, we complete the proof for the asymptotic normality of $T_{1n}$ as well as Part 1.

Proof of Part 2. By Condition $(4)$–$(5)$, we can easily obtain

[TABLE]

Hence, by Slutsky’s theorem,

[TABLE]

Therefore, combining Part 1 and Part 2, we have $\sqrt{n}(\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0})\to_{d}N(0,V)$. Because $\widehat{\boldsymbol{\theta}}=(\widehat{\boldsymbol{\beta}}^{{\mathrm{\scriptscriptstyle T}}},\hat{\lambda})^{{\mathrm{\scriptscriptstyle T}}}$, by the property of multivariate normal distribution, we complete the proof for the case $i=1$ with the covariance matrix $V_{1}$ being the $p\times p$ sub-matrix at the top right-hand corner of $V$.

For the case $i=2$, the proof is a direct extension of the $i=1$ case. We avoid presenting it here.∎

7.6 Proof of Theorem 6

Proof.

It’s enough to show that for any given $\varepsilon$ there exists a positive constant $C$ such that, for $n$ large enough, the following inequality holds:

[TABLE]

where $\boldsymbol{\theta}_{0}$ is the maximizer of $L(\boldsymbol{\theta})=MV(\boldsymbol{\beta})+\lambda(\boldsymbol{\beta}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\beta}-1)$. This means that with probability tending to $1$ there exists a local maximum $\hat{\boldsymbol{\theta}}$ in the ball $\{\boldsymbol{\theta}_{0}+\alpha_{n}{\mathbf{u}}:\|{\mathbf{u}}\|\leq C\}$ so that $\|\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\|=O_{p}(\alpha_{n})$. This implies $\|\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}_{0}\|=O_{p}(\alpha_{n})$. Then,

[TABLE]

where $\boldsymbol{\theta}^{*}$ lies between $\boldsymbol{\theta}_{0}+\alpha_{n}{\mathbf{u}}$ and $\boldsymbol{\theta}_{0}$.

By Condition $(6)$,

[TABLE]

For $A_{2}$, by Condition (7)

[TABLE]

By the Cauchy-Schwarz inequality and Conditions $(6)$ and $(8)$,

[TABLE]

Allowing $\|{\mathbf{u}}\|$ to be large enough, $B_{1}$ and $B_{3}$ are dominated by $B_{2}$ which is less than [math]. This proves (7.12).

∎