arXiv:1608.00088v1 [] 30 Jul 2016

A U -classifier for high-dimensional data under
non-normality
M. Rauf Ahmad1 and Tatjana Pavlenko2
1
2

Department of Statistics, Uppsala University, Uppsala, Sweden

Department of Mathematics, KTH Royal Institute of Technology,
Stockholm, Sweden

Abstract
A classifier for two or more samples is proposed when the data are highdimensional and the underlying distributions may be non-normal. The
classifier is constructed as a linear combination of two easily computable
and interpretable components, the U -component and the P -component.
The U -component is a linear combination of U -statistics which are averages of bilinear forms of pairwise distinct vectors from two independent
samples. The P -component is the discriminant score and is a function
of the projection of the U -component on the observation to be classified.
Combined, the two components constitute an inherently bias-adjusted
classifier valid for high-dimensional data. The simplicity of the classifier
helps conveniently study its properties, including its asymptotic normal
limit, and extend it to multi-sample case. The classifier is linear but its linearity does not rest on the assumption of homoscedasticity. Probabilities
of misclassification and asymptotic properties of their empirical versions
are discussed in detail. Simulation results are used to show the accuracy
of the proposed classifier for sample sizes as small as 5 or 7 and any large
dimensions. Applications on real data sets are also demonstrated.

Keyword: bias-adjusted classifier; U -statistics; discriminant analysis;

1

Introduction

A linear classifier for two or more populations is presented when the data are
high-dimensional and possibly non-normal. Let
xik = (xi1k , . . . , xikp )0 , k = 1, . . . , ni ,
be ni independent and identically distributed random vectors from ith population with distribution function Fi , where E(xik ) = µi , Cov(xik ) = Σi > 0,
i = 1, . . . , g, g ≥ 2, are the mean vector and covariance matrix. Given this set
up, we are interested to construct a linear classifier for high-dimensional, low
1

sample size settings, i.e. p  ni , and when Fi can be non-normal.
Classification and regression are two of the most powerful tools of statistical
analysis, both as main objective of analysis on their own and also as a source of
further investigation. Due to ever growing complexity of data, classification has
particularly attracted a central place in modern statistical analysis. The wave
of large-dimensional data sets in the last few decades and their associated questions of analysis have lead the researchers to substantially think and improve
the classical framework of classification and discrimination.
This paper mainly addresses the classification problem for such a complex
data set up, particularly when the dimension of the multivariate vector, p, may
exceed the number of such vectors, ni , i.e., p  ni (see Sec. 6 for examples).
As the classical theory of classification does not work in this case, mainly due
to the singularity of empirical covariance matrix (see Sec. 2 for more details),
efforts have been made in the literature to offer some potential alternatives.
Bickel and Levina (2004) propose Independence Rule (IR), or naive Bayes rule,
by using only the diagonal of the empirical covariance matrix and compare
it to Fisher’s linear discriminant function (LDF) for the case of two normal
populations. Under certain general conditions on the eigenvalues of the scaled
covariance matrix, they show that IR under independence assumption is comparable to Fisher’s LDF under dependence when the empirical covariance matrix
is replaced with a g-inverse computed from the empirical non-zero eigenvalues
and the corresponding eigenvectors. A regularized discriminant analysis using
Fisher’s LDF is given in Witten and Tibshirani (2011). For a short but useful
review of this and other classification methods for high-dimensional data, see
Mai (2013). Further relevant approaches will be referred to in their appropriate
context in the sequel.
We begin in Sec. 2 with the two-sample U -classifier, giving detailed explanations on the its construction and justification. An extension to the multi-sample
case is given in Sec. 4. Accuracy of the classifier under different parameter settings is shown in Sec. 5, whereas the practical applications on real data sets are
demonstrated in Sec. 6. All technical proofs are deferred to the appendix.

2

The two-sample case

This section is devoted to the construction and justification of the two-sample
classifier along with its properties, asymptotic distribution and misclassification
rate. These are, respectively, the subjects of the next three subsections.

2.1

Construction and motivation of the U -classifier

Let xik = (xik1 , . . . , xikp )0 ∼ Fi be as defined above and πi denote the ith
(unknown) population, i = 1, 2 (g = 2). Let x0 be the new observation to be
classified to either of the two populations where the misclassification errors are
represented by the conditional probabilities
π(i|j) = P(x0 ∈ πi |πj ), i, j = 1, 2, i 6= j.
2

Using the information on p characteristics in each sample, i.e. (xik1 , . . . , xikp )0 ,
we aim to construct a classifier which assigns x0 to πi , i = 1, 2, optimally, i.e.
by keeping π(i|j) as small as possible, as ni , p → ∞, when (i) p may arbitrarily
exceed ni , p  ni , (ii) Fi may not necessarily be normal, and (iii) Σi ’s may
be unequal, Σ1 6= Σ2 . Note that high-dimensional or, as is frequently known,
(n, p)-asymptotic framework is kept general in that it implies both ni → ∞ and
p → ∞ but without requiring the two indices to satisfy any specific relationship
of mutual growth order. It will, however, be shown in the sequel that some of
the results hold even by assuming ni fixed and any arbitrary p.
First, to set the notations, let
xi =

ni
1 X
xik
ni

n

i
1 X
(xik − xi )(xik − xi )
ni − 1

bi =
and Σ

k=1

(1)

k=1

be the usual unbiased estimators of µi and Σi , respectively. The classical twosample linear classifier, assuming equal and known Σ with equal costs and priors
can be expressed, ignoring the constants, as (Seber, 2004, Ch. 6)
1
A(x) = (x1 − x2 )0 Σ−1 x0 − (x01 Σ−1 x1 − x02 Σ−1 x2 ),
2

(2)

where x0 is the point to be classified. Although this classifier is usually constructed under normality assumption, using a ratio of multivariate normal density functions of the two populations and substituting Σ1 = Σ2 (which makes
the classifier linear), it is well-known that Fisher constructed the same classifier without assuming normality, and hence it is also known as Fisher’s linear
discriminant function. It is the most frequently used classifier in practice and
assuming ni > p and normality, the misclassification probability can be computed using the normal distribution function. Obviously, with Σ unknown in
practice, we need to estimate A(x), replacing Σ with its usual pooled estimator,
b pooled = P2 (ni − 1)Σ
b i / P2 (ni − 1), under the homoscedasticity assumpΣ
i=1
i=1
b i are defined above, so that
tion, where Σ
0 b −1
b −1 x0 − 1 (x0 Σ
b −1
b
A(x)
= (x1 − x2 )0 Σ
pooled
1 pooled x1 − x2 Σpooled x2 ).
2

(3)

b i , hence, Σ
b pooled , are
When the data are high-dimensional, i.e., when p > ni , Σ
b
singular and can not be inverted, implying that A(x)
can not be used in this
b
case. To see how the situation develops in this framework, let us first take Σ
out of the classifier in (2) and consider
1
A(x) = (x1 − x2 )0 x0 − (x01 x1 − x02 x2 ).
2

(4)

Assuming x0 ∈ π1 , we immediately note that
E[A(x)|x0 ∈ π1 ]

=
3

1
kµ − µ2 k2 − B,
2 1

(5)

where k · k is the Euclidean norm and B = tr(Σ1 )/2n1 − tr(Σ2 )/2n2 . We note
that, removing the covariance matrix makes the resulting classifier biased with
the bias term, B, composed of the traces of the unknown covariance matrices.
If Σ1 = Σ2 , then B = (n2 − n1 )tr(Σ)/2n1 n2 , so that the classifier is positively
or negatively biased given n2 > n1 or n2 < n1 .
To inherently adjust the classifier for its bias and improve its accuracy, consider the second component of A(x) in Eqn. (4), and note that
E[(x01 x1 − x02 x2 )/2] = B + (µ01 µ1 − µ02 µ2 )/2,
where µ0i µi were used to complete the squared norm in the expectation of complete classifier in (5) and B is the same bias term. Now
x0i xi

ni X
ni
ni
ni X
ni
1 X
1X
1 X
Aikr = 2
Aik + 2
Aikr = Q0i + Q1i ,
= 2
ni
ni
ni
r=1
r=1
k=1

k=1

k=1
k6=r

where Aik = x0ik xik , Aikr = x0ik xir , k 6= r, so that E(Q0i ) = (tr(Σi ) + µ0i µi )/ni
and E(Q1i ) = (1 − 1/ni )µ0i µi . Let Q0 = Q01 − Q02 , Q1 = Q11 − Q12 . Then
E(Q0 ) = 2B + R and E(Q1 ) = (µ01 µ1 − µ02 µ2 ) − R,
where R = µ01 µ1 /n1 − µ02 µ2 /n2 . As R appears with opposite signs in the
two components, adjusting each component by this amount, keeping the total
expectation same, we have E(Q0 ) + R = µ01 µ1 − µ02 µ2 and E(Q1 ) - R = 2B.
Adjusting the corresponding terms in A(x) in (4) similarly, we re-write it as
A0 (x) =

1
1
(x1 − x2 )0 x0 − (Un1 − Un2 ),
p
2

(6)

Pni
where Uni = k6=
r Aikr /pQ(ni ), Q(ni ) = ni (ni − 1), is a one-sample U -statistic
with symmetric kernel, Aikr /p = x0ik xir /p, k 6= r, which is a bilinear form of
two independent components. Now, assuming x0 ∈ π1 and independent of the
elements of both samples, we have E[(x1 − x2 )0 x0 ] = µ01 µ1 − µ01 µ2 , so that
E[A0 (x)|π1 ] = kµ1 − µ2 k2 /2,
without any bias term. Further, with x0 ∈ π1 , A0 (x) is composed of four bilinear
forms, two from sample 1, one from sample 2, and one mixed. By symmetry,
E[A0 (x)|π2 ] = −kµ1 − µ2 k2 /2, and the classifier again consists of four bilinear
forms, two from sample 2, one from sample 1 and one mixed. We, therefore,
define the classification rule for the proposed U -classifier as
Assign x0 to π1 if A0 (x) > 0, otherwise to π2 .
Before we study the properties of A0 (x) in the next section, a few important
remarks are in order. First, A0 (x) is composed of bilinear forms - and we call it
bilinear classifier - where the bi-linearity of the U -component is expressed in the
4

kernels of the two U -statistics and that of the P -component by the projection
of the new observation with respect to the difference between the empirical
centroids of the two independent samples. Further, A0 (x) is entirely composed
of empirical quantities, free of any unknown parameter, so that it can be directly
used in practice using the decision rule stated above. Note also that, A0 (x) is
linear but the linearity does not require homoscedasticity assumption, i.e., it is
linear even if Σ1 6= Σ2 .
Moreover, the first part of A0 (x) is normalized by p, and so are the kernels of
U -statistics in the second part. This will help us derive the limiting distribution
of the classifier for (n, p)-asymptotics under a general multivariate model and
mild assumptions. As a final remark, recall that the formulation of A0 (x) arises
from depriving the original classifier of empirical covariance matrix. Whereas,
this removal of an essential ingredient has its price to be paid, the resulting
classifier can still be justified from a different perspective which has its merit,
particularly for high-dimensional data.
Pni
dki /ni , where dki = kx0 −xki k2
To see this, consider d12 = d1 −d2 , di = k=1
is the Euclidean distance of x0 from sample i = 1, 2, k = 1, . . . , ni . It follows
that d12 has same bias B as for A(x). For all expressions location-invariant, write
b i = Pni D0 Dikr /ni (ni −1), Dikr = xik −xir (Ahmad, 2014b). Now A(x)
e
Σ
=
ikr
k6=r
Pni
b
b
b
d12 −[tr(Σ1 )/n1 −tr(Σ2 )/n2 ], and since k=1 di = (ni −1)tr(Σi )+ni kx0 −xi k2 ,
e
it simplifies to A(x)
= A(x) + B; compare with Eqn. (5).
This implies that A0 (x) can also be constructed using Euclidean distances
and the same bias-adjustment that lead Eqn. (5) to Eqn. (6). This distancebased approach has been discussed in Chan and Hall (2009) and the same is
further evaluated in Aoshima and Yata (2014). Our approach, however, makes
the classifier not only unbiased but also more general as well as convenient to
study and apply in practice.

2.2

Asymptotic distribution of the U -classifier

Given xik ∼ Fi , let zik = xik − µi with E(zik ) = 0, Cov(zik ) = Σi , i = 1, 2.
When we relax normality, we assume the following multivariate model
zik = Λi yik ,

(7)

where yik = (yik1 , . . . , yikp )0 has iid elements with E(yik ) = 0, Cov(yik ) = I,
and Λi is a known p × p matrix of constants such that Λ0i Λi = Ai , Λi Λ0i =
Σi > 0, i = 1, 2. To study the properties of the classifier and its asymptotic
normality, we shall supplement Model (7) with the following assumptions.
4
Assumption 1 E(yiks
) = γ < ∞, γ ∈ R+ , i = 1, 2.

Assumption 2 limp→∞ tr(Σi )/p = O(1), i = 1, 2.
Assumption 3 limp→∞ µ0i Σk µj /p = O(1), i, j, k = 1, 2.
Assumption 4 limp→∞
where

tr(Σa
Σbj )
i
b
tr(Σa
i ⊗Σj )

= 0, a, b = 1, 2, 3, a + b ≤ 4, i, j = 1, 2,

and ⊗ are Hadamard and Kronecker products.
5

Assumption 1 essentially replaces normality. Assumptions 2 is simple and mild,
and as its consequence, limp→∞ tr(Σi Σj )/p2 = O(1), so that a reference to the
assumption may also imply it consequence in the sequel. Assumptions 3 and 4
are needed only to show control of misclassification rate and consistency of the
moments of classifier. Assumption 4 ensures that the moments asymptotically
coincide with those under normality. This assumption is neither directly needed
in practical use of the classifier, nor is it required under normality whence all
terms involving the ratio vanish. The same assumptions will be extended for
multi-sample case in Sec. 4.
Now, continuing to assume x0 independent of all elements of sample 1 (where
it is already independent of all elements of sample 2), we get the following
lemma, proved in Appendix B.1, on the moments of classifier.
Lemma 1 Let the two-sample modified classifier be as given in Eqn. (6). Then,
assuming x0 ∈ πi , we have
E[A0 (x)|πi ]

=

Var[A0 (x)|πi ]

=

1 2
1
kµ − µj k2I =
∆
2p i
2p I
1
1
δi2
δ2
+ 2 kµ1 − µ2 k2Σ−1 = i2 + 2 ∆2Σ−1 ,
2
i
i
p
p
p
p

(8)
(9)

where ∆2M = kµi − µj k2M = (µi − µj )0 M−1 (µi − µj ), i, j = 1, 2, i 6= j, for any
Mp > 0, and
2
tr(Σ2i ) tr(Σi Σj ) X tr(Σ2i )
.
δi2 =
+
+
ni
nj
2ni (ni − 1)
i=1
The moments in Lemma 1 are reported using general notation for x0 ∈ πi so
that they can be easily extended to the multi-sample case later. For the present
case of g = 2, the moments easily reduce to their specific form for i = 1 or i = 2,
where the mean is obviously the same for both, i.e., kµ1 − µ1 k2I /2p.
With the picture of the proposed classifier and its moments relatively clear,
we can express our high-dimensional classification problem precisely as
Cp = {(µi , Σi ) : ∆2I , ∆2Σ−1 }p , i = 1, 2,

(10)

i

where the index p implies the dependence of components Cp on the dimension.
Note that, the second component in Eqn. (9) vanishes under Assumption 3.
The rest of Var[A0 (x)|πi ], and E[A0 (x)|πi ], are uniformly bounded in p, for any
fixed ni , under Assumptions 2. We can thus write
lim E[A0 (x)|πi ]

=

lim Var[A0 (x)|πi ]

=

p→∞

p→∞

1 2
∆
2 0,I
 


1
1
2
δ0,i O
+
+ o(1) ,
n1
n2

(11)
(12)

2
where ∆20,I = limp→∞ ∆2I /p ∈ (0, ∞) and δ0,
= limp→∞ δi2 /p2 ∈ (0, ∞). Now,
i
the variance obviously vanishes when we also allow ni → ∞ along with p → ∞,

6

which immediately gives consistency of the classifier, formally stated in Theorem
2 below and proved in Appendix B.3. Obviously, in practice, the consistency is
expected to hold with unknown parameters replaced by their estimators. We
need to estimate ∆2I and non-vanishing traces in δi2 to estimate the limiting
moments of the classifier. In the following, we define unbiased and consistent
plug-in estimators of these components.
For xik , i = 1, 2, let Dikr = xik − xir , k 6= r, be as defined in Sec. 2.1
with E(Dikr ) = 0, Cov(Dikr ) = 2Σi . Also let Dijkl = xik − xjl using differences of vectors from two independent samples with E(Dijkl ) = µi − µj ,
Cov(Dijkl ) = Σi + Σj , i, j = 1, 2, i 6= j. Extending the strategy of Sec. 2.1
for the estimation of tr(Σi ), we define estimators of the traces involved in δi2 .
Let A2ikrk0 r0 = (D0ikr Dik0 r0 )2 and A2ijkrls = (D0ikr Djls )2 , Aijkl = D0ik Djl , i 6= j,
using within- and between-sample independent vectors, respectively. Then, by
independence, the plug-in estimators of ∆2I /p2 , tr(Σ2i )/p2 and tr(Σi Σj )/p2 are
defined, respectively, as
E0

=

Ei

=

Uni + Unj − 2Uni nj
ni
ni X
ni X
ni X
1 2
1 X
A
2 ikrk0 r 0
12η(ni )
p
0
r=1 0

(13)
(14)

k =1 r =1
k=1
π(k,r,k0 ,r 0 )
n

n

Eij

=

nj

nj

i XX
i X
X
1 2
1
A
2 ijklrs
4η(ni nj )
p
r=1
s=1

(15)

l=1
k=1
π(k,r) π(l,s)

where η(ni ) = ni (ni − 1)(ni − 2)(ni − 3) and η(ni nj ) = ni nj (ni − 1)(nj − 1),
i, j = 1, 2, i 6= j, and π(·) implies all indices pairwise unequal. Further,
Uni nj

nj
ni X
1
1 X
Aijkl ,
=
ni nj
p
k=1 l=1

a two-sample U -statistic, so that, with Uni , i = 1, 2, as one-sample U -statistic
defined after Eqn. (6), E0 estimates ∆2I /p2 = kµi − µj k2 /p2 . E0 , Ei , Eij are
unbiased and location-invariant estimators. The following theorem, proved in
Appendix B.2, shows further that the variances of the ratios of these estimators
to the parameters they estimate are uniformly bounded in p, so that they are
consistent as p → (ni fixed) and also when ni , p → ∞.
Theorem 1 E0 , Ei , Eij , defined in Eqns. (13)-(15), are unbiased estimators
of ∆2I /p2 , tr(Σ2i )/p2 and tr(Σ1 Σ2 )/p2 . Further, under Assumptions 1-4,




1
1
E0
=
O
+
(16)
Var
∆2I
ni
nj


 
Ei
1
Var
=
O
(17)
ni
tr(Σ2i )

7




Eij
tr(Σi Σj )


Ei
Eij
Cov
,
tr(Σ2i ) tr(Σi Σj )
Var



1
1
+
ni
nj
 
1
= O
.
ni



= O

(18)
(19)

The bounds in Theorem 1 suffice for consistency of estimators so that exact
variances and covariances are not needed. These exact moments follow from
Theorem 5 and Lemma 3; see Ahmad (2014b).
P
From Theorem 1, it immediately follow that E0 /E(E0 ) −
→ 1 which gives
1
b 0 (x)|πi ] = E0 , as a consistent estimator of the true mean of
empirical mean, E[A
2
the classifier. Using similar probability convergence of the other two estimators
Ei and Eij , we obtain, by Slutsky’s lemma (van der Vaart, 1998, p 11), the
d 0 (x)|πi ], i.e. δb2 , as a consistent
first component of the empirical variance, Var[A
i
P
estimator of δi2 , i = 1, 2, using plug-in estimators such that δbi2 /δi2 −
→ 1. The
limiting empirical moments, parallel to Eqns. (11)-(12), thus follow as
b 0 (x)|πi ]
lim E[A

p→∞

d 0 (x)|πi ]
lim Var[A

p→∞

1 2
∆
2 0,I
 


1
1
2
= δ0,i
O
+
+ oP (1)
n1
n2
=

(20)
(21)

2
with ∆20,I = limp→∞ E0 /p ∈ (0, ∞), δ0,i
= limp→∞ δbi2 /p2 ∈ (0, ∞). Hence
i
h
b 0 (x)|πi ] − E[A0 (x)|πi ] = oP (1)
(22)
lim E[A
p→∞
i
h
d 0 (x)|πi ] − Var[A0 (x)|πi ] = oP (1).
lim Var[A
(23)
p→∞

The following theorem, proved in Appendix B.3, summarizes both true and
empirical consistency of the classifier.
Theorem 2 Given A0 (x) in Eqn. (6) with its moments as in Lemma 1. Let
x0 ∈ πi . Under Assumptions 1-3, as ni , p → ∞, i = 1, 2,
A0 (x) P (−1)i
−
→
+ oP (1),
∆2I /p
2
with δi defined above. Further, the consistency holds when the moments of the
classifier are replaced with their empirical estimators, given in Eqns. (20)-(21).
The same arguments help us establish asymptotic normality of the classifier as
stated in the following theorem, proved in Appendix B.4.
Theorem 3 Given A0 (x) in Eqn. (6) with its moments as in Lemma 1. Let
x0 ∈ πi . Under Assumptions 1-3, as ni , p → ∞, i = 1, 2,
A0 (x) − E[A0 (x)] D
p
−→ N (0, 1),
Var[A0 (x)]
Further, the normal limit holds when the moments are replaced with their empirical estimators, given in Eqns. (20)-(21).
8

The construction of A0 (x) is of great benefit in proving Theorem 3. Its composition of two parts, each of which in turn a linear combination of two independent
components, reduces the bulk of computational burden. Moreover, the optimality property of U -statistics ensures the minimum variance (efficiency) of the
classifier. A further verification of these properties of the classifier through
simulations is demonstrated in Sec. 5.

3

Estimation of misclassification probabilities

Consider the classification problem Cp in (10) again. Using notation introduced
in the beginning of Sec. 2.1, the optimality of the proposed classifier can be
evaluated by the misclassification rates π(1|2) and π(2|1) with
Z
π(i|j) = P (x ∈ πi |πj ) =
dFj (x),
(24)
Ri

where Fj denotes the distribution function and Ri = {x : x ∈ πi } is the region
of observed data from ith population with R1 ∪ R2 = X , R1 ∩ R2 = ∅, where X
denotes the space of observed x and ∅ is the empty set. Under the assumption
of equal probabilities and equal costs, we are interested to minimize the total
probability of misclassification, say o , given the observed data, i.e.
arg min 0 = arg min[π(1|2) + π(2|1)]/2,
x

x

where the subscript o stands for optimal. Obviously, the ideal minimum can only
be achieved when the parameters are known in which case the (ideal) classifier
takes the form
Aideal
(x ∈ π1 ) = x00 (µ1 − µ2 )/p − (µ01 µ1 − µ02 µ2 )/2p,
0
using the fact that xi and Uni are unbiased estimators of µi and µ0i µi , respectively, i = 1, 2. Now, if Fi are known, say multivariate normal, i.e. xik ∼
Np (µi , Σi ), then, under the homoscedasticity assumption Σ1 = Σ2 = Σ, the
error rate of Aideal
can be expressed as
0


2
kµ
−
µ
k
2 I
,
ideal
= Φ − q 1
o
2 kµ1 − µ2 k2Σ−1
where Φ is the standard normal distribution function. Assuming equal priors,
the best possible performance in this ideal setting, i.e. with µ1 , µ2 , Σ known,
is achieved by Fisher’s linear classifier (equivalently, Bayes rule)
AFisher
(x) = x00 Σ−1 (µ1 − µ2 ) −
0

1
0
(µ + µ2 ) Σ−1 (µ1 − µ2 )
2 1

(Anderson, 2003) with the corresponding misclassification rate given by


q
1
Fisher
2
o
=Φ −
kµ1 − µ2 kΣ ,
2
9

where kµ1 −µ2 k2Σ is the Mahalanobis distance. Denoting Fisher
as a benchmark,
o
the relative performance of Aideal
(x)
can
be
theoretically
evaluated
by using the
0
ratio of the arguments of Φ, say q, where
kµ1 − µ2 k2I
q=
1/2 .
kµ1 − µ2 k2Σ · kµ1 − µ2 k2Σ−1
Bickel and Levina (2004) put forth a nice strategy to compute a bound for an
expression like q, based on Kantorovich inequality (Bernstein, 2009). Following
the same idea, let M be any positive definite symmetric p × p matrix. Then for
any vector υ
4λmin (M) · λmax (M)
kυk2I
≥
2,
2
kυkM · kυk2M−1
[λmin (M) + λmax (M)]
where λmin (M) and λmax (M) denote the smallest and the largest eigenvalues of
(Σ)
=
M, respectively. Applying this inequality to q and denoting the ratio λλmax
min (Σ)
κ (assuming the two extreme eigenvalues bounded away from 0 and ∞), we get
√
2 κ
q≥
,
(25)
1+κ
so that the upper bound on the misclassification probability of Aideal
(x) is
0


√

2 κ −1
ideal
Fisher
o
≤Φ −
Φ
−o
.
1+κ
which essentially depends on κ, the range of non-zero eigenvalues of Σ. We
note that, for moderate κ, the increase in the misclassification rate, induced
by taking the covariance matrix away while constructing Aideal
(x), is not large
0
(see
Fig.
1).
Further, the
relative to the best possible performance, i.e. Fisher
o
upper bound in (25) represents the worst-case scenario so that the empirical
results are expected to be better.
Now, for an alternative flavor of the evaluation of the classifier, while still
continuing to assume normality, let us condition the classifier on the data, i.e.,
A0 (x) = A0 (x)|(xi , Uni ), say, and it immediately follows that


1 0
1
A0 (x) ∼ N
µi (x1 − x2 ) − (Un1 − Un2 ), kx1 − x2 k2Σ−1
i
p
2
i = 1, 2. This, using the standardized version of classifier (Theorem 3), leads to
the actual error rate
"
!
1
µ01 (x1 − x2 ) − (Un1 − Un2 )/2
n =
Φ − p
2
(x1 − x2 )0 Σ1 (x1 − x2 )
!#
µ02 (x2 − x1 ) − (Un1 − Un2 )/2
+ Φ − p
,
(x1 − x2 )0 Σ2 (x1 − x2 )
10

0.5

0.45

0.4

Upper bound of eIdeal
0

0.35
κ=80
0.3

κ=3

0.25
κ=10
0.2

0.15

0.1

0.05

0
0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

eFisher
0

Figure 1: Upper bound on the misclassification probability of Aideal
(x) as a function
0
of Fisher
for normal (thick line) and t5 (dashed line) distributions with κ = 3, 10, 80.
o

where the subscript n denotes the dependence on the observed sample. Using
Theorem 2,


1 0
1
1
P
µ (x1 − x2 ) − (Un1 − Un2 ) − ∆20,I −
→ 0
p i
2
2
kx1 − x2 k2Σ−1 − ∆2Σ−1
i

P

−
→ 0,

i

so that by Slutsky’s lemma (van der Vaart, 1998, p 11) it follows that, as
ni , p → ∞,
"
!
!#
∆20,I
∆20,I
P 1
n −
→
Φ −
+Φ
,
2
∆Σ−1
∆Σ−1
1

2

where the convergence remains (asymptotically) true even for the sample based
classifier A0 (x) since ideal
is the limiting value of n .
o
Finally, we consider a similar evaluation of the classifier under non-normality
which we discuss for the general class of elliptically contoured distributions (including multivariate normal). Let Fi ≡ Ep (h, µi , Σi ), i = 1, 2, where Ep (h, µ, Σ)
denotes a p-dimensional elliptical distribution with density function

|Σ|−1/2 h kx − µk2Σ ,
(26)
where h is a monotone (decreasing) function on [0, ∞) and parameters µ and Σ
are specified as in model (7). Assume now that xik ∼ Ep (h, µi , Σi ) and denote
11

the first two conditional moments of the classifier A0 (x) by
Ei = E [A0 (x)|xi , Uni ] ,

Vi = Var [A0 (x)|xi , Uni ] ,

respectively, i = 1, 2. Wakaki (1994) has discussed Fisher’s linear discriminant function for elliptical distributions and also its robustness particularly in
the context of M -estimation. Our assertion in the following closely follows his
structure.
By conditioning the classifier on the data, i.e. by considering A0 (x) =
A0 (x)|(xi , Uni ) and using its standardized version, we get (see Wakaki, 1994,
Theorem 1.1, p 260)


A0 (x) − Ei
√
P
≤ z =: Q(z)
(27)
Vi
for any z ∈ R, where Γ denotes the gamma function, h is defined in (26), and
Q is the distribution function whose density function is given by
Z
π (p−1)/2 ∞ (p−3)/2

q(z) =
s
h(z 2 + s)ds.
Γ p−1
0
2
The normal distribution is a special case of Ep , for if h(s) = (2π)−p/2 exp(−s/2),
then q(z) reduces to the standard normal. Now, using (27), we can express the
conditional (or actual) misclassification probability as



 
1
E2
E1
n =
+Q √
.
Q −√
2
V1
V2

4

Multi-sample U -classifier

It is obvious from the construction of the two-sample classifier that it can be
easily extended to the multi-sample case. Let
xik ∼ Fi , k = 1, . . . , ni
be iid vectors with E(xik ) = µi , Cov(xik ) = Σi , i = 1, . . . , g ≥ 2. The multisample version of classifier in Eqn. (6) can be expressed as
A0 (x) =

1 0
1
x0 (xi − xk ) − (Uni − Unk ),
p
2

(28)

i, k = 1, . . . , g, k 6= i. Alternatively, to write it in a more explicit form, let
A0i (x) =

1
1 0
x0 xi − Uni
p
2

be the discriminant function for population i, so that the classifier is
A0g (x) = A0i (x) − A0k (x)
12

(29)

for any distinct pair (i, k). The classification rule modifies to
Assign x0 to πi if A0g (x) > 0, i.e., if A0i (x) > A0k (x); otherwise to πk .
To study the properties of the multi-sample classifier and its asymptotic behavior, we first extend the two-sample assumptions for the general case.
Assumption 5 E(x4iks ) = γ < ∞, γ ∈ R+ , i = 1, . . . , g.
Assumption 6 limp→∞ tr(Σi )/p = O(1), i = 1, . . . , g.
Assumption 7 limp→∞ µ0i Σk µj /p = O(1), i, j, k = 1, . . . , g.
Assumption 8 limp→∞
where

tr(Σa
Σbj )
i
b)
tr(Σa
⊗Σ
i
j

= 0, a, b = 1, 2, 3, a + b ≤ 4, i, j = 1, . . . , g,

and ⊗ are Hadamard and Kronecker products.

We begin by the following lemma which generalizes Lemma 1 on the moments
of the two-sample classifier.
Lemma 2 Given A0 (x) in Eqn. (28) or (29). Let x0 ∈ πi . Then
E[A0 (x)|πi ]
Var[A0 (x)|πi ]

= kµi − µk k2 /2p = ∆2I /2p
=

δi2 /p2

+ kµi −

µk k2Σ−1 /p2
i

(30)
=

δi2 /p2

+

∆2Σ−1 /p2 ,
i

(31)

where ∆2M−1 = (µi − µk )0 M(µi − µk ), M = Σ or M = I, and
2

δi2 =

tr(Σ2i ) tr(Σi Σk ) X tr(Σ2i )
+
+
, i, k = 1, . . . , g, i 6= k.
ni
nk
2ni (ni − 1)
i=1

The moment estimators follow obviously from those of the two-sample case given
in Eqns. (13)-(15). Likewise, the consistency of these estimators follows from
Lemma 1. This helps us extend Theorem 2 on the consistency and asymptotic
normality of A0 (x) for the general case as following.
Theorem 4 Given A0 (x) in Eqn. (28) or (29) with its moments in Lemma 1.
Let x0 ∈ πi . Under Assumptions 5-7, as ni , p → ∞, i = 1, . . . , g,
A0 (x)
∆2I /p
A0 (x) − E[A0 (x)]
p
Var[A0 (x)]

P

−
→
D

(−1)i
+ oP (1),
2

−→ N (0, 1),

(32)
(33)

Further, the limits hold when the moments are replaced with their empirical
estimators given in Eqns. (20)-(21).
As the multi-sample case is a straightforward extension of its two-sample counterpart in Sec. 2, we skip many detailed proofs to avoid unnecessary repetitions.
13

5

Simulations

We use simulation results to evaluate the performance of A0 (x) under practical
scenarios, mainly focusing on consistency, asymptotic normality and control of
misclassification under high-dimensional framework. We consider g = 2 case
and generate data from multivariate normal and t distributions, i.e. Fi is either
Np (µi , Σi ), i = 1, 2, or tν (µi , Σi ), ν = 10, i = 1, 2. For each distribution, we
set µ1 = 0 with bp/3c elements of µ2 also 0 and the rest as 1 where b·c denotes
the smallest integer. For Σi , we consider two cases: (1) Both populations have
AR(1) structure, Cov(Xk , Xl ) = κρ|k−l| , ∀ k, l, with σ 2 = 1 for i = 1 and 2,
where ρ = 0.3 for i = 1 and ρ = 0.7 for i = 2, to represent both low and high
correlation structures; (2) The same AR(1) structure for i = 1 with σ 2 = 1,
ρ = 0.5, whereas an unstructured (UN) Σi for i = 2, defined as Σ = (σij )pi,j=1
with σii = 1(1)p and ρij = (i − 1)/p, i 6= j.
For finite-sample performance of the classifier under arbitrarily growing dimension, emphasizing the p  ni , we generate samples of sizes n1 = 5, n2 = 7,
and combine with p = {10, 20, 50, 100, 300, 500, 700, 1000, 3000, 5000, 10000}.
Finally, all results are averages of 1000 simulation runs for each combination
of parameters mentioned above. To additionally observe the effect of large ni ,
the misclassification rates are also presented for n1 = 10, n2 = 12. Similarly, to
assess the classifier for very different sample sizes, we also used n1 = 5, n2 = 25
and n1 = 10, n2 = 50 and observed very similar results, hence not reported here.
Fig. 2 shows the results of asymptotic normality of A0 (x), where the first
two rows are for normal distribution, respectively for AR-AR and AR-UN structures. Each row gives a histogram of A0 (x), with empirical density added to it,
for p = 100, 500 and 1000 (left to right). The last two rows are for multivariate t
distribution with ν = 10. As stated above, the results are carried out for several
other dimensions as well, up to p = 10000, and also for other sample sizes, but
due to similarity of the graphs, only a selection is reported here.
We observe close normal approximation for ni as small as 5 or 7, and the
results for t distribution depict small sample robustness of the classifier to nonnormality. To make the results of the two distributions comparable, the density
axes are scaled to the same height, so that the heavy-tailed behavior of t distribution can be witnessed from a slightly extended range on the x-axis. In general,
it is observed that the increasing dimension does not damage the asymptotic
normality of the classifier even if the data are non-normal.
A similar performance is observed for the control of misclassification rate,
shown in Figs. 3 for n1 = 5, n2 = 7, and 4 for n1 = 10, n2 = 12. The thick
line represents the
p actual error rate under asymptotic normality of the classifier,
i.e. Φ(−E(A0 )/ Var(A0 )), where Φ(·) is the (univariate) normal distribution
function and E(·), Var(·) are the moments of the classifier. This actual error
rate is used as a reference to assess the estimated error rate shown in dashed
(dotted) line for normal (t10 ) distribution. Further, the upper (lower) panel in
each figure is for AR-AR (AR-UN) pair of covariances.
The estimated error closely follows the actual error for n1 = 5, n2 = 7, and
the error rate also converges to zero, showing consistency of the classifier. For t
14

Density

Density

Density

0
A0(x)

0
A0(x)

0
A0(x)

0
A0(x)

1

1

2

2

2

2

4

3

3

4

0.2

0.3

0.4

0.0

0.1

0.2

0.3

0.4

0.0

0.1

0.2

0.3

0.4

0.0

0.0

−2

í1

−1

0.0

−2

í2

−2

0.1

−4

−4

í3

−3

0.1

0.2

0.3

0.4

0.1

0.2

0.3

0.4

0.0

0.1

0.2

0.3

0.4

0.0

0.1

0.2

0.3

0.4

0.0

0.1

0.2

0.3

0.4

−4

−3

−3

−4

−2

−2

−2

−2

−1

−1

0
A0(x)

0
A0(x)

0
A0(x)

0
A0(x)

1

1

2

2

2

2

4

3

3

4

0.0

0.1

0.2

0.3

0.4

0.0

0.1

0.2

0.3

0.4

0.0

0.1

0.2

0.3

0.4

0.0

0.1

0.2

0.3

0.4

−4

−4

−3

−3

−2

−2

−2

−2

−1

−1

0

0
A0(x)

0
A0(x)

A0(x)

A0(x)

0

1

1

2

2

2

2

n1 = 5, n2 = 7, p = 100, 500, 1000 (L to R each row) and covariance structures AR-AR (Rows 1, 3) & AR-UN (Rows 2, 4).

Figure 2: Asymptotic normality of A0 for two-class case with data from multivariate normal (Rows 1-2) & t (Rows 3-4) distributions,

Density

Density
Density
Density
Density

Density
Density
Density
Density

15

3

4

4

3

0.35

0.35

0.30

0.30

0.25

0.25
Theoretical
Estimated: MVN
Estimated: MVT

Theoretical
Estimated: MVN
Estimated: MVT

e(1|0)

0.20

e(1|0)

0.20

0.15

0.15

0.10

0.10

0.05

0.05

0.00

0.00
200

400

600

800

1000

200

400

p

0.30

800

1000

600

800

1000

0.30

0.25

0.25
Theoretical
Estimated: MVN
Estimated: MVT

Theoretical
Estimated: MVN
Estimated: MVT
0.20

e(0|1)

0.20

e(1|0)

600
p

0.15

0.15

0.10

0.10

0.05

0.05

0.00

0.00
200

400

600

800

1000

200

p

400
p

Figure 3: Theoretical (thick line) and Estimated error rates of A0 (x) for twoclass case with data from multivariate normal (dashed line) and t (dotted line)
distributions, n1 = 5, n2 = 7, p = 10, 20, 50, 100, 200, 300, 500 and covariance structures
AR-AR (upper panel) & AR-UN (lower penal).

distribution with n1 = 5, n2 = 7, the estimated error rates are relatively
higher than under normality, but with ni increased only by 5, a discernable
difference in the performance of the classifier is observed in Fig. 4. Note that,
the x-axis in Figs. 3-4 is truncated at p = 500 since the misclassification rates
already converge to 0 by this value and remain so for larger p.

6

Applications

We apply A0 (x) on two large data sets for g = 2 and 3. With moderate sample
sizes (77 and 102), we use K = 3-fold CV for evaluation (see Dudoit et al, 2002).
16

0.35

0.35

0.30

0.30

0.25

0.25
Theoretical
Estimated: MVN
Estimated: MVT

Theoretical
Estimated: MVN
Estimated: MVT

e(1|0)

0.20

e(1|0)

0.20

0.15

0.15

0.10

0.10

0.05

0.05

0.00

0.00
0

100

200

300

400

500

0

100

200

p

300

0.30

500

0.30

0.25

0.25
Theoretical
Estimated: MVN
Estimated: MVT

Theoretical
Estimated: MVN
Estimated: MVT
0.20

e(0|1)

0.20

e(1|0)

400

p

0.15

0.15

0.10

0.10

0.05

0.05

0.00

0.00
0

100

200

300

400

500

0

p

100

200

300

400

500

p

Figure 4: Theoretical (thick line) and Estimated error rates of A0 (x) for twoclass case with data from multivariate normal (dashed line) and t (dotted line)
distributions, n1 = 10, n2 = 12, p = 10, 20, 50, 100, 200, 300, 500 and covariance structures AR-AR (upper panel) & AR-UN (lower penal).

Let L and T denote the learning and test sets. We randomly divide data set
into K classes of roughly equal size where T consists of K − 1 classes with Kth
class held out as test data. The procedure is repeated K times, each time with a
different test class, and a misclassification rate is computed for each repetition.
The evaluation criterion is the average misclassification rate over all repetitions.
For kth fold of CV, let nki (L), nki (T ) and mkij (T ) be, respectively, the sample
sizes for learning and test data in sample i and the number of misclassified
observations from class i into class j, i, j = 1, . . . , g = 2 or 3, k = 1, . . . , K = 3.
Let ek (i|j) be the estimated misclassification rate, an estimator of π(i|j) in (24),
for kth rotation, i.e. ek (i|j) = mkij (T )/nk (T ), where nk (T ) = nki (T ) + nkj (T ).

17

For g = 3, we do the same procedure for each of three pairs and compute overall
misclassification rate. For details on the used and other data sets, see Statnikov
et al (2005) and also Shipp et al (2002) and Armstrong et al (2002).
Example 1: DLBCL Data The Diffuse Large B-cell Lymphoma (DLBCL)
data belongs to a study of lymphoid malignancy in adults. The analysis reported
here consists of p = 5469 gene expressions studied on pre-treatment biopsies
from two independent groups of 77 patients, one with DLBCL (n1 = 58), the
other with follicular lymphoma (FL) (n2 = 19).
For a 3-fold CV, we randomly divide the data into three groups of sizes 26,
26, 25 with n1 (L) = 52, n1 (T ) = 25 and nk (L) = 51, nk (T ) = 26 for k = 2, 3.
By coding the populations as 1 (DLBCL) and 2 (FL), the misclassifications
observed from the three rotations of CV, i.e. mk12 and mk21 , for k = 1, 2, 3, are
m112 = 3, m121 = 1;

m212 = 6, m221 = 0;

m312 = 2, m321 = 3,

so that the overall misclassification rate is computed as 15/77. Although of
relatively less importance, due to randomly sampled folds for cross-validation,
we also report the sample sizes for each fold as following:
K=1:

n11 (L) = 38, n12 (L) = 14, n11 (T ) = 20, n12 (T ) = 5

K=2:

n21 (L) = 40, n22 (L) = 11, n21 (T ) = 18, n22 (T ) = 8

K=3:

n31 (L) = 38, n32 (L) = 13, n31 (T ) = 20, n32 (T ) = 6

Example 2: Leukemia Data The data set pertains to a study of patients
with acute lymphoblastic leukemia (ALL) carrying a chromosomal translocation involving mixed-lineage leukemia (MLL) gene. The analysis reported here
consists of p = 11225 gene expression profiles of leukemia cells from n2 = 24
patients diagnosed with B-precursor ALL carrying an MLL translocation and
compared to a group of n3 = 20 individual diagnosed with conventional Bprecursor without MLL translocation. In addition, there is a third group of a
random sample of n1 = 28 with acute myelogenous leukemia (AML).
For a 3-fold cross-validation, we randomly divide the data into three equal
groups of size 24 each and use K −1 = 2 classes of total nk (L) = 48 observations
in learning set and nk (T ) = 24 in the test set every time, k = 1, 2, 3. The rest
of the procedure is the same as explained in the Example 1 above. We obtain
the following misclassifications for the three folds of cross-validation:
K=1:
K=2:

m112 = 1, m121 = 0;
m212 = 0, m221 = 1;

m113 = 2, m131 = 1;
m213 = 2, m231 = 0;

m123 = 0, m132 = 1
m223 = 1, m232 = 0

K=3:

m312 = 0, m321 = 0;

m313 = 0, m331 = 0;

m323 = 0, m332 = 0.

This gives an overall misclassification rate 9/72. The sample sizes used in each
rotation are also reported below.
K=1:

n11 (L) = 21, n12 (L) = 15, n13 (L) = 12; n11 (T ) = 7, n12 (T ) = 9, n13 (T ) = 8

K=2:

n21 (L) = 18, n22 (L) = 14, n23 (L) = 16; n21 (T ) = 10, n22 (T ) = 10, n23 (T ) = 4

K=3:

n31 (L) = 17, n32 (L) = 19, n33 (L) = 12; n31 (T ) = 11, n32 (T ) = 5, n33 (T ) = 8

18

7

Discussion and conclusions

A U -classifier for high-dimensional and possibly non-normal data is proposed.
The threshold part of the classifier, called U -component, is a linear combination
of two bivariate U -statistics of computed from the two independent samples.
The discriminant function part, called P -component, forms an inner product
between the observation to be classified and the difference of the mean vectors
of the corresponding independent samples. It results into a computationally
simple classifier which is linear without requiring the underlying covariance matrices to be equal. A multi-class extension with same properties is also given.
The classifier is unbiased, consistent and asymptotically normal, under a general multivariate model, including (but not necessarily) the multivariate normal
distribution. Rapid convergence of the misclassification rate of the classifier is
shown for very small sample sizes and non-normal distributions, under mild and
practically justifiable assumptions. The performance of the classifier, in terms
of its consistency, asymptotic normality and control of misclassification rate, is
shown through simulations for normal and non-normal distributions with sample sizes as small as 5 or 7 and for arbitrary large dimension.
We apply the classifier to genetics and microarray data sets, some of the
most popular areas for classification analysis. To emphasize the role of highdimensionality, we demonstrate that the use, accuracy, and validity of the classifier does not rest on any form of data pre-processing as is usually shown in the
literature. In other words, a data set measured in large dimension can be directly used for the classifier without any pre-requisites of reducing the dimension
through sorting or clustering or other means.

A

Moments of Quadratic and Bilinear Forms

0
0
For zik in (7), let Aik = z0ik Σi zik = yik
A2 yik and Aijkl = z0ik zjl = yik
Λi Λj yjl ,
k 6= l, be a quadratic and a bilinear form of independent components with
0
Aik = yik
yik = Qik for Σi = I. As all terms involving Ai eventually vanish
under Assumption 4, we write Ai = A for simplicity. Theorem 5 gives basic
moments of quadratic and bilinear forms which we extend in Lemma 3. Proofs
of these results are tedious but simple, therefore skipped (see Ahmad, 2014b).

Theorem 5 For Aik and Aijkl , as defined above, we have
E Q2ik

2

=

2tr(Σ2i ) + [tr(Σi )]2 + M1

(34)

E A2ik

2

=

2tr(Σ4i ) + [tr(Σ2i )]2 + M2

(35)

E (Aik Ajk )

E A4ijkl

E Qik Qjk A2ijkl

=

2tr(Σ3i Σj ) + tr(Σ2i )tr(Σi Σj ) + M2

=
=

2

(36)

2

6tr(Σi Σj ) + 3 [tr(Σi Σj )] + M3
2

4tr(Σi Σj ) +

4tr(Σ3i )tr(Σj )

19

+ [tr(Σi )]

(37)
2

tr(Σ2j )

+ M4

(38)

P
P
with M1 = γtr(A A), M2 = γtr(A2 A2 ), M3 = 6γtr(A2 A2 )+γ 2 ps=1 pt=1 A4st ,
2
3
M4 = 2γtr(Σi )tr(A
A) + 4γtr(A
A) + γtr(A ADA) and D = diag(A).
Moreover, E(Aik ) = tr(Σi ), E(A2ikr ) = tr(Σ2i ) and Cov(Aik , Aikr ) = 0.

Lemma 3 Let zit be as given above with zit , ziu independent if t 6= u. Then
E[z0it ziu z0it ziv z0iu Σi ziv ]

= tr(Σ4i )

E[z0it ziu z0iw ziu z0it ziv z0iw ziv ]
E[(z0it ziu )2 z0it Σi zit ]
Var(z0it ziu z0iv ziu )
Cov[(z0it ziu )2 , (z0it ziv )2 ]
E[(z0it zju )2 z0it Σj zit ]

=

Var(z0it zju z0iv zju )

=

Cov[(z0jt ziu )2 , (z0jt ziv )2 ]
E(z0jt ziu z0jt ziv z0iu Σj ziv )
Cov[(z0it ziu )2 , z0it Σj zit ]
Cov(z0it Σi ziu , z0it Σj ziu )
E[(z0iu ziv )2 z0it Σj zit ]

=
=
=
=

(39)

tr(Σ4i )
2tr(Σ4i )
2tr(Σ4i )
2tr(Σ4i )

(40)
+
+

[tr(Σ2i )]2
[tr(Σ2i )]2

+ M2

(41)

+ M2

(42)

+ M2

(43)

2



2

2tr{(Σi Σj ) } + tr(Σi Σj ) + M2

2
2tr{(Σi Σj )2 } + tr(Σi Σj ) + M2

(44)
(45)

2

(46)

=

tr{(Σi Σj ) }

(47)

=

2tr(Σ3i Σj ) +
tr{(Σi Σj )2 }

=

=
=

2tr{(Σi Σj ) } + M2
2

M2

(48)
(49)

2

tr(Σi Σj )tr(Σi ) ,

(50)

where E[(z0it ziu )2 z0it ziu z0it ziv ], E[z0it ziu z0it ziv z0it Σi zit ], E[z0it ziu z0it ziv z0it Σi ziu ],
E[(z0it ziu )2 z0it Σi ziu ], E[(z0it ziu )2 z0it Σj ziu ], E[(z0it ziu )2 z0it ziv z0iu ziv ] all vanish.

B

Main Proofs

B.1

Proof of Lemma 1

Let x0 ∈ π1 . With x0 independent of both samples, E[A0 (x)] is trivial. For
variance, ignoring p for simplicity, we begin with
Var[x00 (x1 − x2 )] = E[x00 (x1 − x2 )]2 − [µ01 (µ1 − µ2 )]2 .
Since E[x00 (x1 − x2 )]2 = E[x00 (x1 − x2 )(x1 − x2 )0 x0 ], we immediately get



Σ1
Σ2
E[x00 (x1 − x2 )]2 = tr (Σ1 + µ1 µ01 )
+
+ (µ1 − µ2 )(µ1 − µ2 )0 ,
n1
n2
so that
tr(Σ21 ) tr(Σ1 Σ2 ) µ01 Σ1 µ1 µ01 Σ2 µ1
+
+
+
+(µ1 −µ2 )0 Σ1 (µ1 −µ2 ).
n1
n2
n1
n2
P2
0
Now Var(Un1 − Un2 ) =
i=1 Var(Uni ). For Uni with h(xik , xir ) = xik xir
0
0
(Serfling, 1980, Ch. 5)), h1 (xik ) = xik µi with ξ1 = Var[h1 (xik )] = µi Σi µi , and
h2 (·) = h(·) with ξ2 = Var(Aik ) = tr(Σ2i ) + 2µ0i Σi µi , so that

Var[x00 (x1 −x2 )] =

Var(Uni ) =

2
ni (ni − 1)

[2(ni − 2)ξ1 + ξ2 ] =
20

2tr(Σ2i )
4µ0i Σi µi
+
, i = 1, 2.
ni (ni − 1)
ni

For Cov[x00 (x1 − x2 ), Un1 − Un2 ], Cov(x00 x2 , Un1 ) = 0 = Cov(x00 x1 , Un2 ), by
independence, where it immediately follows that Cov(x00 xi , Uni ) = 2µ0i Σi µi /ni ,
i = 1, 2. Combining all results and simplifying gives Eqn. (9).

B.2

Proof of Theorem 1

The unbiasedness is trivial. For Var(E0 ), Var(Uni ), i = 1, 2, are given in Sec.
B.1. For Var(Unij ), h(xik , xjl ) = x0ik xjl with h10 = µ0j xik , h01 = µ0i xjl so that
ξ10 = Var[h10 (·)] = µ0j Σi µj and ξ10 = Var[h10 (·)] = µ0i Σj µi . Also h11 (·) = h(·)
with ξ11 = Var[h11 (·)] = µ0i Σj µi + µ0j Σi µj + tr(Σi Σj ). Hence (Lee, 1990)
Var(Uni nj ) =


1 
ni µ0i Σj µi + nj µ0j Σi µj + tr(Σi Σj )
ni nj p2

where Cov(Uni , Uni nj ) = 2µ0j Σi µi /ni p2 , Cov(Unj , Uni nj ) = 2µ0i Σj µj /nj p2 and
Cov(Uni , Unj ) = 0 by independence. Var(E0 /p) can now be approximated as
Var(E0 /p) = 2tr(Σ20ij )/p2 + 4(µi − µj )0 Σ0ij (µi − µj )/p2 , i, j = 1, 2, i 6= j.
With second term vanishing under Assumption 3 and first term bounded in p
under Assumption 2, Var(E0 ) reduces to O(1/n1 +1/n2 )2 as p → ∞, so that the
consistency follows immediately as ni , p → ∞. The bound in (16) also follows
trivially. As Ei and Eij , are also one- and two-sample U -statistics with higher
order kernels, we essentially follow the same strategy as for E0 . First, from
Theorem 5 and Lemma 3, it can be shown that (see Ahmad, 2014b, Ch. 2)
Var(Ei )

=

4 h 3
(2ni − 9n2i + 9ni − 16)tr(Σ4i )
η(ni )p4

i
+ (n2i − 3ni + 8)[tr(Σ2i )]2 + M2 O(n3i ) + M3 O(n2i ) (51)
h
2
Var(Eij ) =
(n − 1)tr{(Σ1 Σ2 )2 } + [tr(Σ1 Σ2 )]2
(ni − 1)(nj − 1)p4
i
+ M2 O(n) + M3 O(1)
(52)
h
i
4
ni tr(Σ3i Σj ) + M2 O(ni )
(53)
Cov(Ei , Eij ) =
ni (ni − 1)p4
n = ni + nj , i, j = 1, 2, i 6= j, M2 , M3 are as in Theorem 5 and Cov(Ei , Ej ) = 0.
As terms involving M ’s vanish under Assumption 4, the consistency and the
bounds (by Cauchy-Schwarz inequality) follow the same way as for E0 . Note
also that, the terms involving M ’s are exactly zero under normality in which
case the same results follow even more conveniently.

B.3

Proof of Theorem 2

The proof essentially follows from that of Theorem 1 without much new computations. In particular, the first part, assuming true parameters known, is
21

trivial. For the second part with empirical estimators, the (ni , p)-consistency
P
of estimators proved in Sec. B.2 implies that E0 /E(E0 ) −
→ 1, and the same
holds for Ei , Eij . Plugging these estimators in the moments of A0 (x) and using
P
\
Slutsky’s lemma, δbi2 /δi2 −
→ 1 so that Var[A
0 (x)] = Var[A0 (x] + oP (1), and the
consistency follows similarly as with known parameters.

B.4

Proof of Theorem 3

e0 (x) = [A0 (x)|x0 ∈ π1 ] − E[A0 (x)|x0 ∈ π1 ], where
Write A
e0 (x) = [x00 (x1 − x2 ) − µ01 (µ1 − µ2 )] − [(Un − µ01 µ1 ) − (Un − µ02 µ2 )]/2,
A
1
2
bn be the projection of U
en = Un −µ0 µi , i = 1, 2.
ignoring p for simplicity. Let U
i
i
i
i
0
Then (Serfling, 1980, Ch. 5) g1 (x1k ) = h1 (·) − µ1 µ1 = (X1k − µ1 )0 µ1 for Un1 ,
and similarly g1 (x2l ) for Un2 , with E[g(·)] = 0 in both cases, so that
n1
n2
X
2 X
0
bn − U
bn = 2
U
(X
−
µ
)
µ
−
(X2l − µ2 )0 µ2 ,
1k
1
1
1
2
n1
n2
k=1

l=1

en = U
bn + oP (1), i = 1, 2. With x0 ∈ π1 and independence of samples,
where U
i
i
e0 (x) results into a sum of two independent components, each
this projection of A
an average of iid variables (van der Vaart, 1998). Taking p into account, the
asymptotic normality follows by the CLT under Assumptions 1-3 as ni , p → ∞.
Acknowledgement We are sincerely thankful to Prof. Thomas Mikosch,
Department of Mathematics, University of Copenhagen, for careful perusal of
the paper, discussion and encouragement. The research of Tatjana Pavlenko is
partially supported by the Swedish Research Council, Grant No. C0595201.

References
Ahmad, MR (2014a). A U -statistic approach for a high-dimensional two-sample
mean testing problem under non-normality and Behrens-Fisher setting. Annals of the Institute of Statistical Mathematics, 66(1), 33-61.
Ahmad, MR (2014b). Location-invariant and non-invariant tests for highdimensional covariance matrices under normality and non-normality. Technical report, 2014:4, Uppsala University, Sweden.
Anderson, TW (2003). An introduction to multivariate statistical analysis. Third
edition. Wiley, New York.
Aoshima, M, and K Yata (2014). A distance-based, misclassification rate adjusted classifier for multiclass, high-dimensional data. Annals of the Institute
of Statistical Mathematics, 66, 983-1010.

22

Armstrong, SA, JE Staunton, LB Silverman et al. (2002). MLL translocations
specify a distinct gene expression profile that distinguishes a unique leukemia.
Nature Genetics, 30, 41-47.
Bernstein, DS (2009). Matrix mathematics: Theory, facts and formulas. Second
edition. Princeton University Press, Princeton.
Bickel, P and E Levina (2004). Some theory for Fisher’s linear discriminant
function, ’naive Bayes’, and some alternatives when there are many more
variables than observations. Bernoulli, 10(6), 989-1010.
Chan, Y-B and P Hall (2009). Scale adjustments for clasifiers in highdimensional, low sample size settings. Biometrika, 96(2), 469-478.
Dehling, H, M Denker, and W Philipp (1984). Invarianee principles for yon
Mises and U -statistics. Z. Wahrsch. verw. Gebiete 67, 139-167.
Dudoit, S, J Fridlyand, and TP Speed (2002). Comparison of discriminant methods for the classification of tumors using gene expression data. Journal of
Americal Statistical Association,97, 77-87.
Johnson, RA and DW Wichern (2007). Applied multivariate statistical analysis,
6th ed. Pearson Education, NJ.
Koroljuk, VS, and YV Borovskich (1994). Theory of U -statistics. Kluwer Academic Press, Dordrecht.
Lee, AJ (1990). U-statistics: Theory and Practice. CRC Press, Boca Raton.
Mai, Q (2013). A review of discriminant analysis in high dimensions. WIREs
Computational Statistics, 5, 190-197.
Mikosch, T (1993). Weak invariance principle for weighted U -statistics with
varying kernels. Journal of Multivariate Analysis,47, 82-102.
Mikosch, T (1994). Weak invariance principle for weighted U -statistics. Journal
of Theoretical Probability,7, 147-173.
Seber, GAF (2004). Multivariate observations. Wiley, New York.
Serfling, RJ (1980). Approximation theorems of mathematical statistics. Wiley,
Weinheim.
Shipp, MA, KN Ross, P Tamayo, et al. (2002). Diffuse large B-cell lymphoma outcome prediction by gene-expression profiling and supervised machine learning. Nature Medicine, 8(1), 68-74.
Statnikov, A, I Tsamardinos, Y Dosbayev, and CF Aliferis (2005). GEMS: A
system for automated cancer diagnosis and biomarker discovery from microarray gene expression data. International Journal of Medical Informatics, 74,
491-503.
23

van der Vaart, AW (1998). Asymptotic statistics. Cambridge University Press,
Cambridge.
Wakaki, H (1994). Discriminant analysis under elliptical populations. Hiroshima
Journal of Mathematics, 24, 257-298.
Witten, D and R Tibshirani (2011). Penalized classification using Fisher’s linear
discrimininat. JRSS B, 73, 753-772.

24

