Linear classifier design under heteroscedasticity in Linear Discriminant Analysis

TL;DR

This paper introduces the Gaussian Linear Discriminant (GLD), a new linear classifier designed for heteroscedastic data that minimizes Bayes error and is computationally efficient, outperforming traditional LDA and comparable classifiers in accuracy and speed.

Contribution

The paper derives the GLD for heteroscedastic data and proposes a local neighbourhood search algorithm to enhance robustness, significantly reducing training time compared to existing methods.

Findings

GLD outperforms traditional LDA in heteroscedastic scenarios.

GLD and LNS require up to 150 times less training time than SVM.

Proposed classifiers achieve comparable or better accuracy across diverse datasets.

Abstract

Under normality and homoscedasticity assumptions, Linear Discriminant Analysis (LDA) is known to be optimal in terms of minimising the Bayes error for binary classification. In the heteroscedastic case, LDA is not guaranteed to minimise this error. Assuming heteroscedasticity, we derive a linear classifier, the Gaussian Linear Discriminant (GLD), that directly minimises the Bayes error for binary classification. In addition, we also propose a local neighbourhood search (LNS) algorithm to obtain a more robust classifier if the data is known to have a non-normal distribution. We evaluate the proposed classifiers on two artificial and ten real-world datasets that cut across a wide range of application areas including handwriting recognition, medical diagnosis and remote sensing, and then compare our algorithm against existing LDA approaches and other linear classifiers. The GLD is shown to outperform the original LDA procedure in terms of the classification accuracy under heteroscedasticity. While it compares favourably with other existing heteroscedastic LDA approaches, the GLD requires as much as 60 times lower training time on some datasets. Our comparison with the support vector machine (SVM) also shows that, the GLD, together with the LNS, requires as much as 150 times lower training time to achieve an equivalent classification accuracy on some of the datasets. Thus, our algorithms can provide a cheap and reliable option for classification in a lot of expert systems.

Tables

Table 1. Table 1 : List and Characteristics of Datasets

Dataset	Label	$n$	$d$	$K$
D1	(a)	$2000$	$8$	$2$
D2	(b)	$2000$	$4$	$2$
Liver	(c)	$345$	$6$	$2$
Shuttle	(d)	$58000$	$9$	$7$
Vowels	(e)	$990$	$10$	$11$
Zernike Moments	(f)	$2000$	$47$	$10$
Image Segmentation (Statlog)	(g)	$2310$	$19$	$7$
Spambase	(h)	$4601$	$37$	$2$
Wine Quality (White)	(i)	$4898$	$11$	$7$
Pen Digits	(j)	$5620$	$64$	$10$
Satellite (Statlog)	(k)	$6435$	$36$	$6$
Letters	(l)	$20000$	$16$	$26$

Table 2. Table 2 : Average Bayes Error ( % percent \% )

Dataset	LDA	C-HLD	R-HLD-1	R-HLD-2	GLD
(a)	$0.0397$	$0.0382$	$0.0383$	$0.0361$	0.0360
(b)	$0.0774$	$0.0749$	$0.0749$	$0.0740$	0.0739
(c)	$0.9981$	0.9838	0.9838	0.9838	0.9838
(d)	0.0001	0.0001	0.0001	0.0001	0.0001
(e)	$0.0339$	0.0326	0.0326	0.0326	0.0326
(f)	$0.0054$	$0.0051$	0.0048	0.0048	$0.0050$
(g)	$0.0037$	0.0029	0.0029	0.0029	0.0029
(h)	$0.0253$	$\text{0.022}8$	0.0228	0.0228	0.0228
(i)	$0.0162$	$0.0201$	$0.0156$	$0.0155$	0.0154
(j)	0.0002	0.0002	0.0002	0.0002	0.0002
(k)	$0.0046$	0.0039	0.0039	0.0039	0.0039
(l)	0.0007	0.0007	0.0007	0.0007	0.0007

Table 3. Table 3 : Average Classification Accuracy ( % percent \% )

Dataset	LDA	C-HLD	R-HLD-1	R-HLD-2	GLD	LNS	SVM
(a)	$76.00$	$77.18$	$77.00$	$78.48$	78.65	$78.57$	$77.47$
(b)	$76.87$	$77.93$	$77.93$	$78.17$	78.37	$78.00$	$77.70$
(c)	$67.83$	$63.19$	$62.32$	$62.03$	$63.77$	$68.12$	68.70
(d)	$94.10$	$96.60$	$96.74$	$96.73$	$96.59$	97.91	$84.39$
(e)	$73.64$	$74.14$	$74.44$	$74.44$	$74.14$	$75.66$	76.77
(f)	$84.00$	$83.90$	$84.10$	$84.15$	84.80	$84.00$	$81.90$
(g)	$94.33$	$94.59$	$94.59$	$94.63$	$94.59$	$94.89$	96.15
(h)	$88.76$	$88.29$	$88.26$	$88.15$	$88.26$	90.28	$85.68$
(i)	$53.41$	$46.59$	$53.37$	$53.33$	$53.55$	$\text{54.1}4$	$51.88$
(j)	$96.74$	$96.99$	$96.97$	$96.98$	$97.01$	$97.41$	97.84
(k)	$85.69$	$86.06$	$86.06$	$86.03$	$86.08$	$86.65$	86.85
(l)	$81.67$	$81.87$	$81.83$	$81.78$	$81.88$	$82.25$	85.39

Table 4. Table 4 : Average Training Time (s)

Dataset	LDA	C-HLD	R-HLD-1	R-HLD-2	GLD	LNS	SVM
(a)	0.001	$0.161$	$0.140$	$0.139$	$0.002$	$0.181$	$23.192$
(b)	0.001	$0.142$	$0.121$	$0.121$	$0.002$	$0.060$	$0.721$
(c)	0.001	$0.155$	$0.1415$	$0.1337$	$0.003$	$0.028$	$2.673$
(d)	0.037	$3.531$	$3.023$	$3.012$	$0.089$	$43.32$	$4623.138$
(e)	0.036	$11.099$	$9.409$	$9.751$	$0.167$	$2.075$	$1.173$
(f)	0.387	$123.662$	$123.649$	$121.906$	$1.955$	$110.694$	$23.126$
(g)	0.128	$37.320$	$30.876$	$37.875$	$0.488$	$2.143$	$21.775$
(h)	0.101	$10.437$	$7.729$	$7.474$	$0.753$	$36.83$	$804.574$
(i)	0.017	$4.257$	$3.691$	$3.750$	$0.080$	$5.928$	$914.257$
(j)	0.638	$10.099$	$9.358$	$9.171$	$0.915$	$168.19$	$409.38$
(k)	0.304	$18.067$	$17.842$	$17.912$	$0.858$	$13.919$	$311.202$
(l)	0.835	$73.050$	$64.022$	$65.414$	$3.245$	$37.202$	$109.232$

Equations

$${\mathcal{C}}^{\ast }(\text{x})=\arg \underset{{\mathcal{C}}_k}{\max }p({\mathcal{C}}_k|\text{x}),k\in \{1,2,...,K\}$$

$$p({\mathcal{C}}_k|\text{x})=\frac{p(\text{x}|{\mathcal{C}}_k)\times p({\mathcal{C}}_k)}{p(\text{x})},k\in \{1,2\}$$

$$\lambda (\text{x})=\frac{p(\text{x}|{\mathcal{C}}_1)}{p(\text{x}|{\mathcal{C}}_2)}$$

$$p(\textbf{x}|\mathcal{C}_{k})=\frac{1}{\sqrt{(2\pi)^{d}\det(\bm{\Sigma}_{k})}}\exp\bigg{[}-\frac{1}{2}(\textbf{x}-\bar{\textbf{x}}_{k})^{T}\bm{\Sigma}_{k}^{-1}(\textbf{x}-\bar{\textbf{x}}_{k})\bigg{]}$$

$$\ln \lambda (\text{x})$$

$$\displaystyle=\frac{1}{2}\ln\frac{\det\bm{\Sigma}_{2}}{\det\bm{\Sigma}_{1}}+\frac{1}{2}\bigg{[}(\textbf{x}-\bar{\textbf{x}}_{2})^{T}\bm{\Sigma}_{2}^{-1}(\textbf{x}-\bar{\textbf{x}}_{2})-(\textbf{x}-\bar{\textbf{x}}_{1})^{T}\bm{\Sigma}_{1}^{-1}(\textbf{x}-\bar{\textbf{x}}_{1})\bigg{]}$$

$$(\text{x}-{\bar{\text{x}}}_2{)}^T{\mathbf{\Sigma }}_2^{-1}(\text{x}-{\bar{\text{x}}}_2)-(\text{x}-{\bar{\text{x}}}_1{)}^T{\mathbf{\Sigma }}_1^{-1}(\text{x}-{\bar{\text{x}}}_1)\ge \ln \frac{{\tau }^2\det {\mathbf{\Sigma }}_1}{\det {\mathbf{\Sigma }}_2}$$

$${\text{x}}^T{\mathbf{\Sigma }}_x^{-1}({\bar{\text{x}}}_1-{\bar{\text{x}}}_2){⋛}_{{\mathcal{C}}_2}^{{\mathcal{C}}_1}\ln \tau +\frac{1}{2}({\bar{\text{x}}}_1^T{\mathbf{\Sigma }}_x^{-1}{\bar{\text{x}}}_1-{\bar{\text{x}}}_2^T{\mathbf{\Sigma }}_x^{-1}{\bar{\text{x}}}_2)$$

$$S=\frac{{\text{w}}^T({\bar{\text{x}}}_1-{\bar{\text{x}}}_2)({\bar{\text{x}}}_1-{\bar{\text{x}}}_2{)}^T\text{w}}{{\text{w}}^T{\mathbf{\Sigma }}_x\text{w}}$$

$$p_e={\pi }_{1}p(y<w_0|{\mathcal{C}}_1)+{\pi }_{2}p(y\ge w_0|{\mathcal{C}}_2)$$

$$\displaystyle\textbf{w}=\big{[}s_{1}\bm{\Sigma}_{1}+s_{2}\bm{\Sigma}_{2}\big{]}^{-1}(\bar{\textbf{x}}_{1}-\bar{\textbf{x}}_{2})$$

$$w_0={\mu }_1-s_1{\sigma }_1^2={\mu }_2+s_2{\sigma }_2^2$$

$$\displaystyle\textbf{w}=\big{[}s\bm{\Sigma}_{2}+(1-s)\bm{\Sigma}_{1}\big{]}^{-1}(\bar{\textbf{x}}_{1}-\bar{\textbf{x}}_{2})$$

$$w_0={\mu }_1-(1-s){\sigma }_1^2={\mu }_2+s{\sigma }_2^2$$

$$\textbf{w}=\big{[}s\bm{\Sigma}_{1}+(1-s)\bm{\Sigma}_{2}\big{]}^{-1}(\bar{\textbf{x}}_{1}-\bar{\textbf{x}}_{2})$$

$$w_0=\frac{s{\mu }_2{\sigma }_1^2+(1-s){\mu }_1{\sigma }_2^2}{s{\sigma }_1^2+(1-s){\sigma }_2^2}$$

$$\mathrm{d}\textbf{w}=\big{(}s\bm{\Sigma}_{2}+(1-s)\bm{\Sigma}_{1}\big{)}^{-1}(\bm{\Sigma}_{1}-\bm{\Sigma}_{2})\big{(}s\bm{\Sigma}_{2}+(1-s)\bm{\Sigma_{1}}\big{)}^{-1}(\bar{\textbf{x}}_{1}-\bar{\textbf{x}}_{2})\mathrm{d}s$$

$${\mathcal{C}}^{\ast }(\text{x})=\{\begin{array}{ll}{\mathcal{C}}_1& \text{if}y={\text{w}}^T\text{x}\ge w_0\\ {\mathcal{C}}_2& \text{if}y={\text{w}}^T\text{x}<w_0\end{array}$$

$${\mu }_1={\text{w}}^T{\bar{\text{x}}}_1{\mu }_2={\text{w}}^T{\bar{\text{x}}}_2{\sigma }_1^2={\text{w}}^T{\mathbf{\Sigma }}_1\text{w}{\sigma }_2^2={\text{w}}^T{\mathbf{\Sigma }}_2\text{w}$$

$$p(y<w_0|{\mathcal{C}}_1)$$

$$\displaystyle=\int_{-\infty}^{w_{0}}\frac{1}{\sqrt{2\pi}\sigma_{1}}\exp\bigg{[}-\frac{(\zeta-\mu_{1})^{2}}{2\sigma_{1}^{2}}\bigg{]}d\zeta=1-Q\bigg{(}\frac{w_{0}-\mu_{1}}{\sigma_{1}}\bigg{)}$$

$$p(y\geq w_{0}|\mathcal{C}_{2})=\int_{w_{0}}^{\infty}\frac{1}{\sqrt{2\pi}\sigma_{2}}\exp\bigg{[}-\frac{(\zeta-\mu_{2})^{2}}{2\sigma_{2}^{2}}\bigg{]}d\zeta=Q\bigg{(}\frac{w_{0}-\mu_{2}}{\sigma_{2}}\bigg{)}$$

$$p_{e}=\pi_{1}\big{[}1-Q(z_{1})\big{]}+\pi_{2}\big{[}Q(z_{2})\big{]}$$

$$z_1=\frac{w_0-{\mu }_1}{{\sigma }_1}\text{and}{z}_2=\frac{w_0-{\mu }_2}{{\sigma }_2}$$

$$\nabla p_{e}(\textbf{w},w_{0})=\bigg{[}\frac{\partial p_{e}}{\partial\textbf{w}^{T}},\frac{\partial p_{e}}{\partial w_{0}}\bigg{]}^{T}=\textbf{0}$$

$$\frac{\partial p_{e}}{\partial\textbf{w}}=\pi_{1}\bigg{(}\frac{1}{\sqrt{2\pi}}e^{-z_{1}^{2}/2}\frac{\partial z_{1}}{\partial\textbf{w}}\bigg{)}-\pi_{2}\bigg{(}\frac{1}{\sqrt{2\pi}}e^{-z_{2}^{2}/2}\frac{\partial z_{2}}{\partial\textbf{w}}\bigg{)}$$

$$\frac{\partial z_1}{\partial \text{w}}=\frac{-{\sigma }_1{\bar{\text{x}}}_1-z_1{\mathbf{\Sigma }}_1\text{w}}{{\sigma }_1^2}\text{and}\frac{\partial z_2}{\partial \text{w}}=\frac{-{\sigma }_2{\bar{\text{x}}}_2-z_2{\mathbf{\Sigma }}_2\text{w}}{{\sigma }_2^2}$$

$$\frac{\partial p_{e}}{\partial\textbf{w}}=\frac{1}{\sqrt{2\pi}}\bigg{[}-\pi_{1}e^{-z_{1}^{2}/2}\bigg{(}\frac{\sigma_{1}\bar{\textbf{x}}_{1}+z_{1}\bm{\Sigma}_{1}\textbf{w}}{\sigma_{1}^{2}}\bigg{)}+\pi_{2}e^{-z_{2}^{2}/2}\bigg{(}\frac{\sigma_{2}\bar{\textbf{x}}_{2}+z_{2}\bm{\Sigma}_{2}\textbf{w}}{\sigma_{2}^{2}}\bigg{)}\bigg{]}$$

$$\frac{\partial p_{e}}{\partial w_{0}}=\pi_{1}\bigg{(}\frac{1}{\sqrt{2\pi}}e^{-z_{1}^{2}/2}\frac{\partial z_{1}}{\partial w_{0}}\bigg{)}-\pi_{2}\bigg{(}\frac{1}{\sqrt{2\pi}}e^{-z_{2}^{2}/2}\frac{\partial z_{2}}{\partial w_{0}}\bigg{)}$$

$$\frac{\partial z_1}{\partial w_0}=\frac{1}{{\sigma }_1}\text{and}\frac{\partial z_2}{\partial w_0}=\frac{1}{{\sigma }_2}$$

$$\frac{\partial p_{e}}{\partial w_{0}}=\frac{\pi_{1}}{\sqrt{2\pi}}\bigg{(}\frac{1}{\sigma_{1}}e^{-z_{1}^{2}/2}\bigg{)}-\frac{\pi_{2}}{\sqrt{2\pi}}\bigg{(}\frac{1}{\sigma_{2}}e^{-z_{2}^{2}/2}\bigg{)}$$

$$\bigg{(}\frac{\pi_{2}z_{2}}{\sigma_{2}^{2}}e^{-z_{2}^{2}/2}\bm{\Sigma}_{2}-\frac{\pi_{1}z_{1}}{\sigma_{1}^{2}}e^{-z_{1}^{2}/2}\bm{\Sigma}_{1}\bigg{)}\textbf{w}=\bigg{(}\frac{\pi_{1}}{\sigma_{1}}e^{-z_{1}^{2}/2}\bigg{)}\bar{\textbf{x}}_{1}-\bigg{(}\frac{\pi_{2}}{\sigma_{2}}e^{-z_{2}^{2}/2}\bigg{)}\bar{\textbf{x}}_{2}$$

Full text

††journal: Expert Systems with Applications

1 Introduction

In many applications one encounters the need to classify a given object under one of a number of distinct groups or classes based on a set of features known as the feature vector. A typical example is the task of classifying a machine part under one of a number of health states. Other applications that involve classification include face detection, object recognition, medical diagnosis, credit card fraud prediction and machine fault diagnosis.

A common treatment of such classification problems is to model the conditional density functions of the feature vector [Ng & Jordan, 2002]. Then, the most likely class to which a feature vector belongs can be chosen as the class that maximises the a posteriori probability of the feature vector. This is known as the maximum a posteriori (MAP) decision rule.

Let $K$ be the number of classes, $\mathcal{C}_{k}$ be the $k$th class, x be a feature vector and $\mathcal{D}_{k}$ be training samples belonging to the $k$th class $(k\in\{1,2,...,K\})$. The MAP decision rule for the classification task is then to choose the most likely class of x, $\mathcal{C}^{*}(\textbf{x})$ given as:

[TABLE]

We assume for the moment that there are only $K=2$ classes, i.e. binary classification (we consider multi-class classification in a later section). Then, using Bayes’ rule, the two posterior probabilities can be expressed as:

[TABLE]

It is often the case that the prior probabilities $p(\mathcal{C}_{1})$ and $p(\mathcal{C}_{2})$ are known, or else they may be estimable from the relative frequencies of $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ in $\mathcal{D}$ where $\mathcal{D}=\mathcal{D}_{1}\cup\mathcal{D}_{2}$. Let these priors be given by $\pi_{1}$ and $\pi_{2}$ respectively for class $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$. Then, the likelihood ratio defined as:

[TABLE]

is compared against a threshold defined as $\tau=\pi_{2}/\pi_{1}$ so that one decides on class $\mathcal{C}_{1}$ if $\lambda(\textbf{x})\geq\tau$ and class $\mathcal{C}_{2}$ otherwise.

Linear Discriminant Analysis (LDA) proceeds from here with two basic assumptions [Izenman, 2009, Chapter 8]:

The conditional probabilities $p(\textbf{x}|\mathcal{C}_{1})$ and $p(\textbf{x}|\mathcal{C}_{2})$ have multivariate normal distributions. 2. 2.

The two classes have equal covariance matrices, an assumption known as homoscedasticity.

Let $\bar{\textbf{x}}_{1}$, $\bm{\Sigma}_{1}$ be the mean and covariance matrix of $\mathcal{D}_{1}$ and $\bar{\textbf{x}}_{2}$, $\bm{\Sigma}_{2}$ be the mean and covariance of $\mathcal{D}_{2}$ respectively. Then, for $k\in\{1,2\}$,

[TABLE]

where $d$ is the dimensionality of $\mathcal{X}$, which is the feature space of x. Given the above definitions of the conditional probabilities, one may obtain a log-likelihood ratio given as:

[TABLE]

which is then compared against $\ln\tau$ so that $\mathcal{C}_{1}$ is chosen if $\ln\lambda(\textbf{x})\geq\ln\tau$, and $\mathcal{C}_{2}$ otherwise. Thus, the decision rule for classifying a vector x under class $\mathcal{C}_{1}$ can be rewritten as:

[TABLE]

In general, this result is a quadratic discriminant. However, a linear classifier is often desired for the following reasons:

A linear classifier is robust against noise since it tends not to overfit [Mika et al., 1999]. 2. 2.

A linear classifier has relatively shorter training and testing times [Yuan et al., 2012]. 3. 3.

Many linear classifiers allow for a transformation of the original feature space into a higher dimensional feature space using the kernel trick for better classification in the case of a non-linear decision boundary [Bishop, 2006, Chapter 6].

By calling on the assumption of homoscedasticity, i.e. $\bm{\Sigma}_{1}=\bm{\Sigma}_{2}=\bm{\Sigma}_{x}$, the original quadratic discriminant given by (6) for classifying a given vector x decomposes into the following linear decision rule:

[TABLE]

Here, $\bm{\Sigma}_{x}^{-1}(\bar{\textbf{x}}_{1}-\bar{\textbf{x}}_{2})$ is a vector of weights denoted by w and $\ln\tau+\frac{1}{2}(\bar{\textbf{x}}_{1}^{T}\bm{\Sigma}_{x}^{-1}\bar{\textbf{x}}_{1}-\bar{\textbf{x}}_{2}^{T}\bm{\Sigma}_{x}^{-1}\bar{\textbf{x}}_{2})$ is a threshold denoted by $w_{0}$. This linear classifier is also known as Fisher’s Linear Discriminant. If only the weight vector w is required for dimensionality reduction, w may be obtained by maximising Fisher’s criterion [Fisher, 1936], given by:

[TABLE]

where $\bm{\Sigma}_{x}=n_{1}\bm{\Sigma}_{1}+n_{2}\bm{\Sigma}_{2}$ and $n_{1}$, $n_{2}$ are the cardinalities of $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ respectively.

LDA is the optimal Bayes’ classifier for binary classification if the normality and homoscedasticity assumptions hold [Hamsici & Martinez, 2008] [Izenman, 2009, Chapter 8]. It demands only the computation of the dot product between w and x, which is a relatively computationally inexpensive operation.

As a supervised learning algorithm, LDA is performed either for dimensionality reduction (usually followed by classification) [Barber, 2012, Chapter 16], [Buturovic, 1994, Duin & Loog, 2004, Sengur, 2008], or directly for the purpose of statistical classification [Fukunaga, 2013, Chapter 4][Izenman, 2009, Mika et al., 1999]. LDA has been applied to several problems such as medical diagnosis e.g. [Sharma & Paliwal, 2008, Coomans et al., 1978, Sengur, 2008, Polat et al., 2008], face and object recognition e.g. [Song et al., 2007, Chen et al., 2000, Liu et al., 2007, Yu & Yang, 2001] and credit card fraud prediction e.g. [Mahmoudi & Duman, 2015]. The widespread use of LDA in these areas is not because the datasets necessarily satisfy the normality and homoscedasticity assumptions, but mainly due to the robustness of LDA against noise, being a linear model [Mika et al., 1999]. Since the linear Support Vector Machine (SVM) can be quite expensive to train, especially for large values of $K$ or $n$ ($n=n_{1}+n_{2}$), LDA is often relied upon [Hariharan et al., 2012].

Yet, practical implementation of LDA is not without problems. Of note is the small sample size (SSS) problem that LDA faces with high-dimensional data and much smaller training data [Sharma & Paliwal, 2015, Lu et al., 2003]. When $d\gg n$, the scatter matrix $\Sigma_{x}$ is not invertible, as it is not full-rank. Since the decision rule as given by (7) requires the computation of the inverse of $\Sigma_{x}$, the singularity of $\Sigma_{x}$ makes the solution infeasible. In works by, for example, [Liu et al., 2007, Paliwal & Sharma, 2012], this problem is overcome by taking the Moore-Penrose pseudo-inverse of the scatter matrix, rather than the ordinary matrix inverse. Sharma & Paliwal [2008] use a gradient descent approach where one starts from an initial solution of w and moves in the negative direction of the gradient of Fisher’s criterion (8). This method avoids the computation of an inverse altogether. Another approach to solving the SSS problem involves adding a scalar multiple of the identity matrix to the scatter matrix to make the resulting matrix non-singular, a method known as regularised discriminant analysis [Friedman, 1989, Lu et al., 2003].

However, for a given dataset that does not satisfy the homoscedasticity or normality assumption, one would expect that modifications to the original LDA procedure accounting for these violations would yield an improved performance. One such modification, in the case of a non-normal distribution, is the mixture discriminant analysis [Hastie & Tibshirani, 1996, McLachlan, 2004, Ju et al., 2003] in which a non-normal distribution is modelled as a mixture of Gaussians. However, the parameters of the mixture components or even the number of mixture components, are usually not known a priori. Other non-parametric approaches to LDA that remove the normality assumption involve using local neighbourhood structures [Cai et al., 2007, Fukunaga & Mantock, 1983, Li et al., 2009] to construct a similarity matrix instead of the scatter matrix $\Sigma_{x}$ used in LDA. However, these approaches aim at linear dimensionality reduction, rather than linear classification. Another modification, in the case of a non-linear decision boundary between $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$, is the Kernel Fisher Discriminant (KFD) [Mika et al., 1999, Zhao et al., 2009, Polat et al., 2008]. KFD maps the original feature space $\mathcal{X}$ into some other space $\mathcal{Y}$ (usually higher dimensional) via the kernel trick [Mika et al., 1999]. While the main utility of the kernel is to guarantee linear separability in the transformed space, the kernel may also be employed to transform non-normal data into one that is near-normal.

Our proposed method differs from the above approaches in that we primarily consider violation of the homoscedasticity assumption, and do not address the SSS problem. We seek to provide a linear approximation to the quadratic boundary given by (6) under heteroscedasticity without any kernel transformation; we note that several heteroscedastic LDA approaches have been proposed to this effect. Nevertheless, for reasons which we highlight in the next section, our contributions in this paper are stated explicitly as follows:

We propose a novel linear classifier, which we term the Gaussian Linear Discriminant (GLD), that directly minimises the Bayes error under heteroscedasticity via an efficient optimisation procedure. This is presented in Section 3. 2. 2.

We propose a local neighbourhood search method to provide a more robust classifier if the data has a non-normal distribution (Section 4).

2 Related Work

Under the heteroscedasticity assumption, many LDA approaches have been proposed among which we mention [Fukunaga, 2013, Chapter 4],[Zhang & Liu, 2008, McLachlan, 2004, Duin & Loog, 2004, Decell Jr & Marani, 1976, Decell & Mayekar, 1977, Malina, 1981, Loog & Duin, 2002]. As it is known that Fisher’s criterion (whose maximisation is equivalent to the LDA derivation described in the Introduction section) only takes into account the difference in the projected class means, existing heteroscedastic LDA approaches tend to obtain a generalisation on Fisher’s criterion. In the work of [Loog & Duin, 2002], for instance, a directed distance matrix (DDM) known as the Chernoff distance, which takes into account the difference in covariance matrices between the two classes as well as the projected class means, is maximised instead of Fisher’s criterion (8). The same idea employing the Chernoff criterion is used by [Duin & Loog, 2004]. A wider class of Bregman divergences including the Bhattacharya distance [Decell Jr & Marani, 1976] and the Kullback-Leibler divergence [Decell & Mayekar, 1977] have also been used for heteroscedastic LDA, as Fisher’s criterion can be considered a special case of these measures when the covariance matrices of the classes are equal.

However, most of these approaches aim at linear dimensionality reduction, which involves finding a linear transformation that transforms the original data into one of reduced dimensionality, while at the same time maximising the discriminatory information between the classes. Our focus with this paper, however, is not on dimensionality reduction, but on obtaining a Bayes optimal linear classifier for binary classification assuming that the covariance matrices are not equal. As far as we know, the closest work to ours in this regard are the works by [Marks & Dunn, 1974, Anderson & Bahadur, 1962, Peterson & Mattson, 1966, Fukunaga, 2013]

Obtaining the Bayes optimal linear classifier involves minimising the probability of misclassification $p_{e}$ as given by:

[TABLE]

where $y=\textbf{w}^{T}\textbf{x}$. Unfortunately, there is no closed-form solution to the minimisation of (9) [Anderson & Bahadur, 1962]. Thus, an iterative procedure is inevitable in order to obtain the Bayes optimal linear classifier.

In the work of [Marks & Dunn, 1974], for example, the iterative procedure described is to solve for w and $w_{0}$ as given by

[TABLE]

by obtaining the optimal values of $s_{1}$ and $s_{2}$ via systematic trial and error. We denote this heteroscedastic LDA procedure by R-HLD-2, for the reason that the two parameters $s_{1}$ and $s_{2}$ are chosen at random.

[Anderson & Bahadur, 1962] makes the observation that if the weight vector w and the threshold $w_{0}$ are both multiplied by the same positive scalar, the decision boundary remains unchanged. Therefore, by multiplying (10) through by the scalar $s_{1}+s_{2}$, w and $w_{0}$ can be put in the form of:

[TABLE]

Still, the optimal value of $s$ has to be chosen by systematic trial and error. We denote this heteroscedastic LDA approach by R-HLD-1, for the reason that only one parameter $s$ is chosen at random. As we show in the next section, $s$ is unbounded. Therefore, the difficulty faced by this approach is that $s$ has to be chosen from the interval $(-\infty,\infty)$, so that the probability of finding the optimal $s$ for a given dataset is low, without extensive trial and error to limit the choice of $s$ to some finite interval $[a,b]$.

To avoid the unguided trial and error procedure in [Marks & Dunn, 1974, Anderson & Bahadur, 1962], [Peterson & Mattson, 1966] and [Fukunaga, 2013, Chapter 4] propose a theoretical approach described below:

Change $s$ from [math] to $1$ with small step increments $\Delta s$. 2. 2.

Evaluate w as given by:

[TABLE] 3. 3.

Evaluate $w_{0}$ as given by:

[TABLE] 4. 4.

Compute the probability of misclassification $p_{e}$. 5. 5.

Choose w and $w_{0}$ that minimise $p_{e}$.

We refer to this procedure as C-HLD, for the reason that the optimal $s$ is constrained in the interval $[0,1]$.

However, we highlight two main problems with the above C-HLD procedure:

There is no obvious choice of the step rate $\Delta s$. Too small a value of $\Delta s$ will demand too many matrix inversions in Step 2, as there will be too many $s$ values. On the other hand, if $\Delta s$ is too large, the optimal $s$ may not be refined enough, and the w obtained may not be optimal. Specifically, the change in w that results from a small change in $s$ is given as:

[TABLE]

which can affect the classification accuracy. 2. 2.

The solution obtained this way is only locally optimal as $s$ is bounded in the interval $[0,1]$. As we show in the next section, $s$ is actually unbounded. When there is a class imbalance [Xue & Titterington, 2008], the optimal $s$ may be found outside the interval $[0,1]$ which can lead to poor classification accuracy.

Our proposed algorithm, which is described in the next section, unlike the trial and error approach by [Marks & Dunn, 1974, Anderson & Bahadur, 1962], has a principled optimisation procedure, and unlike [Fukunaga, 2013, Peterson & Mattson, 1966] does not encounter the problem of choosing an inappropriate $\Delta s$, nor restricts $s$ to the interval $[0,1]$. Consequently, our proposed algorithm achieves a far lower training time than the C-HLD, R-HLD-1 and R-HLD-2, for roughly the same classification accuracy.

3 Gaussian Linear Discriminant

Let $\textbf{w}\in\mathbb{R}^{d}$ be a vector of weights, and $w_{0}\in\mathbb{R}$, a threshold such that:

[TABLE]

Since x is assumed to have a multivariate normal distribution in classes $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$, $y$ has a mean of $\mu_{1}$ and a variance of $\sigma_{1}^{2}$ for class $\mathcal{C}_{1}$ and a mean of $\mu_{2}$ and a variance of $\sigma_{2}^{2}$ for class $\mathcal{C}_{2}$ given as:

[TABLE]

With reference to the Bayes error of (9), the individual misclassification probabilities can be expressed as:

[TABLE]

and

[TABLE]

where $Q(^{.})$ is the Q-function. Therefore, the Bayes error to be minimised may be rewritten as:

[TABLE]

where

[TABLE]

Our aim is to find a local minimum of $p_{e}$. A necessary condition is for the gradient of $p_{e}$ to be zero, i.e.,

[TABLE]

From (9), it can be shown that:

[TABLE]

From (20), however, we obtain the following:

[TABLE]

Therefore,

[TABLE]

It can similarly be shown from (9) that,

[TABLE]

Again, from (20),

[TABLE]

Therefore,

[TABLE]

Now, equating the gradient $\nabla p_{e}(\textbf{w},w_{0})$ to zero, the following set of equations are obtained:

[TABLE]

[TABLE]

Substituting (29) into (28) yields:

[TABLE]

Then the vector w can be given by:

[TABLE]

It will be noted however that (31) is still in terms of $w_{0}$, so that an explicit representation of $w_{0}$ in terms of w is needed from (29) to substitute in $z_{1}$ and $z_{2}$ in (31). This is where our approach most significantly differs from [Fukunaga, 2013]. Solving for $w_{0}$ from (29) results in the following quadratic:

[TABLE]

which can be simplified to:

[TABLE]

where $\tau$ is given as before as $\tau=\pi_{2}/\pi_{1}$. If $\tau$ is defined and not equal to zero, and $\sigma_{1}^{2}\not=\sigma_{2}^{2}$ (since $\bm{\Sigma}_{1}\not=\bm{\Sigma}_{2}$ for heteroscedastic LDA), (33) can be shown to have the following solutions:

[TABLE]

Nevertheless, since there are two solutions to $w_{0}$ in (34), a choice has to be made as to which of them is substituted into (31). To eliminate one of the solutions, we consider the second-order partial derivative of $p_{e}$ with respect to $w_{0}$ evaluated at $w_{0}$ as given by (34), and determine under what condition it is greater than or equal to zero. This is a second-order necessary condition for $p_{e}$ to be a local minimum. From (27), it can be shown that:

[TABLE]

We denote this second-order derivative by $h$. We then consider all possibilities of $z_{1}$ and $z_{2}$ (which are the variables in (35) that depend on $w_{0}$) under three cases, and analyse the sign of $h$ in each.

Case 1

$z_{2}\leq 0$ and $z_{1}\geq 0$: then $h$ is trivially non-positive.

Case 2

$z_{2}\geq 0$ and $z_{1}\leq 0$: then $h$ is trivially non-negative.

Case 3

$z_{2}>0$ and $z_{1}>0$ or $z_{2}<0$ and $z_{1}<0$: then $h$ is non-negative if and only if

[TABLE]

i.e.,

[TABLE]

It will be noted that the right-hand side of the inequality (37) is identically zero, as can be seen from (32). Therefore, the condition under which $h$ is greater than or equal to zero is when:

[TABLE]

Note also that Case 2 necessarily satisfies (38) so that we consider (38) as the general inequality for the non-negativity of $h$ for all cases, and thus for $w_{0}$ to be a local minimum.

Now, when one considers the two solutions of $w_{0}$ in (35), only the solution given by:

[TABLE]

satisfies the inequality of (38), i.e., only this choice of $w_{0}$ corresponds to a local minimum. The proof of this is given in the appendix.

We may then substitute this expression of $w_{0}$ into (31) so that (31) is in terms of w only. Even so, w has to be solved for iteratively. This is because (31) has no closed-form solution since $\mu_{1},\mu_{2},\sigma_{1},\sigma_{2}$ are themselves functions of w. As the iterative procedure requires an initial choice of w, we use Fisher’s choice of the weight vector as given by:

[TABLE]

as our initial solution. Again, we mention that $n_{1}$ and $n_{2}$ are the cardinalities of $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$. After a number of such iterative updates, the optimal $w_{0}$ is then solved for from (39). This algorithm, known as the Gaussian Linear Discriminant (GLD), is described in detail in Algorithm 1.

Note that by multiplying both w of (31) and $w_{0}$ proportionally by $c=(\sigma_{1}z_{2}-\sigma_{2}z_{1})/\sigma_{1}\sigma_{2}$ (due to (38), $c$ is non-negative and hence the discrimination criterion given by (15) is not changed), the GLD may be viewed in terms of the optimal solution of (11), where

[TABLE]

which is unbounded given the inequality of (38). However, unlike [Anderson & Bahadur, 1962, Marks & Dunn, 1974], $s$ is not chosen by systematic trial and error, and unlike [Fukunaga, 2013], $s$ is not varied between [math] and $1$ at small step increments. Instead, since $s$ is a function of w and $w_{0}$, our algorithm may be interpreted as obtaining increasingly refined values of $s$ by improving upon w and $w_{0}$ starting from Fisher’s solution, as is described in Algorithm 1.

3.1 Stopping Criteria

The GLD algorithm may be terminated under any of the following conditions:

When the change in the objective function $p_{e}$ remains within a certain tolerance $\epsilon_{1}$ for a number of consecutive iterations. 2. 2.

When the change in the norm of w remains within a certain tolerance $\epsilon_{2}$ for a number of consecutive iterations. 3. 3.

When the gradient of $p_{e}$ as given by (21) remains within a certain tolerance $\epsilon_{3}$ for a number of consecutive iterations. 4. 4.

After a fixed number of iterations $I$, if convergence is slow.

At the end of the algorithm, the final solution may be chosen either as the solution to which the iterations converge, or the solution corresponding to the minimum $p_{e}$ found in the iterative updates.

3.2 Multiclass Classification

Suppose now that there are $K>2$ classes in the dataset $\mathcal{D}$, then the classification problem may be reduced to a number of binary classification problems. The two main approaches usually taken for this reduction are the One-vs-All (OvA) and One-vs-One (OvO) strategies [Bishop, 2006, Hsu & Lin, 2002].

3.2.1 One-vs-All (OvA)

In OvA, one trains a classifier to discriminate between one class and all other classes. Thus, there are $K$ different classifiers. An unknown vector x is then tested on all $K$ classifiers so that the class corresponding to the classifier with the highest discriminant score is chosen. However, with respect to the proposed GLD algorithm, this is an ill-suited approach. This is because the collection of all other classes on one side of the discriminant will not necessarily have a normal distribution, and could in fact be multimodal, if the means are well-separated. Since our algorithm is built on strong normality assumptions of the data on each side of the discriminant, the GLD, as has been formulated, is expected to perform poorly.

3.2.2 One-vs-One

In OvO, a classifier is trained to discriminate between every pair of classes in the dataset, ignoring the other $K-2$ classes. Thus, there are $K(K-1)/2$ unique classifiers that may be constructed. Again, an unknown vector x is tested on all $K(K-1)/2$ classifiers. The predicted classes for all the classifiers are then tallied so that the class that occurs most frequently is chosen. This is equivalent to a majority vote decision. In a lot of cases, however, there is no clear-cut winner, as more than one class may have the highest number of votes. In such a case, the most likely class is often chosen randomly between those most frequently occurring classes. The GLD provides a more appropriate means for breaking such ties, by making use of the minimised Bayes error $p_{e}$ for each classifier. Specifically, one may instead use a weighted voting system, where the count of every predicted class is weighted by $1-p_{e}$, since $p_{e}$ provides an appropriate measure of uncertainty associated with each classifier output. Thus, the decision rule is reduced to choosing the maximum weighted vote among the $K$ classes.

Note that even though the GLD minimises the Bayes error for each classifier, the overall Bayes error for a multiclass problem may not be minimised by using multiple binary classifiers.

4 Non-normal Distributions

So far, the fundamental assumption that has been used to derive the GLD is that the data in each class has a normal distribution. Thus, for an unknown non-normal distribution, the linear classifier we have obtained does not minimise the Bayes error for that unknown distribution. We argue, however, that if this unknown distribution is nearly-normal [Mudholkar & Hutson, 2000], then a more robust linear classifier may be found in some neighbourhood of the GLD. For this reason, we use a local neighbourhood search algorithm to explore the region in $\mathbb{R}^{d+1}$ around the GLD to obtain the classifier that minimises the number of misclassifications on the training dataset. We do this by perturbing each of the $d+1$ vector elements in the optimal $\tilde{\textbf{w}}=[w_{0},\textbf{w}^{T}]^{T}$ obtained from the GLD procedure by a small amount $\delta\tilde{w}_{i}$. After every perturbation, the resulting classifier is evaluated on the test dataset. This procedure is repeated as described in Algorithm 2 until the stopping criterion is satisfied.

The algorithm is terminated after a certain maximum number of iterations $R$ is reached. Additionally, one may perform an early termination if after a predefined number of iterations $r_{max}$, there is no improvement in the minimum number of misclassifications on the training dataset that has been found in the search.

5 Experimental Validation

We validate our proposed algorithm on two artificial datasets denoted by D1 and D2, as well as on ten real-world datasets taken from the University of California, Irvine (UCI) Machine Learning Repository. These datasets are shown in Table 1, and cut across a wide range of applications including handwriting recognition, medical diagnosis, remote sensing and spam filtering. D1 and D2 are normally distributed with different covariance matrices. For D1, we generate $1000$ samples for class $\mathcal{C}_{1}$ and $2000$ samples for class $\mathcal{C}_{2}$ using the following Gaussian parameters:

[TABLE]

For D2, we generate $2000$ samples for class $\mathcal{C}_{1}$ and $4000$ samples for class $\mathcal{C}_{2}$ using the following Gaussian parameters:

[TABLE]

The above Gaussian parameters are slightly modified from the two class data used by [Fukunaga, 2013] and [Xue & Titterington, 2008] in order to make the sample means less separated.

For each dataset in Table 1, we perform $10$-fold cross validation. We run $20$ different trials. On each training dataset, we evaluate the minimum Bayes error achievable by our proposed algorithm averaged over all $10$ folds and $20$ trials. If there are more than two classes, we use OvO, and calculate the mean Bayes error over all $K(K-1)/2$ discriminants. As we are interested only in linear classification, we compare the performance of the GLD with the original LDA as well as the heteroscedastic LDA procedures by [Fukunaga, 2013],[Anderson & Bahadur, 1962] and [Marks & Dunn, 1974] as described in Section 2 in terms of the Bayes error (9). For the sake of brevity, we denote these three heteroscedastic LDA algorithms by the annotations earlier introduced: C-HLD, R-HLD-1 and R-HLD-2 respectively. These results are shown in Table 2.

Moreover, for each of the test datasets, we evaluate the average classification accuracy for each of LDA, C-HLD, R-HLD-1, R-HLD-2, GLD and GLD with local neighbourhood search (LNS). We also compare the performance of these LDA approaches to the SVM. These results are shown in Table 3, while the average training times of the algorithms are shown in Table 4.

We estimate the prior probabilities based on the relative frequencies of the data in each class in the dataset, and the stopping criterion for the GLD is thus: we stop if the gradient of w change is less than or equal to $\epsilon_{3}=10^{-6}$, or else we terminate our algorithm after $I=20$ iterations and choose the solution corresponding to the minimum $p_{e}$. Also, for the LNS procedure, we perturb each vector element by $10\%$ of its absolute value, i.e. $\delta\tilde{w}_{i}=0.1|\tilde{w}_{i}|$, and we run for $R$=1000 iterations, terminating prematurely if $r_{max}=0.1R$. We use a step size of $\Delta s=0.001$ for the C-HLD algorithm, and run $1000$ trials for R-HLD-1 and R-HLD-2. All the parameters used in the experiments are optimised via cross-validation. Note that if the sample covariance matrix is singular, we use the Moore-Penrose pseudo-inverse.

6 Results and Discussion

For real-world datasets, the covariance matrices of the classes are rarely equal, therefore the homoscedasticity assumption in LDA does not hold. Our results in Table 2 confirm that LDA does not minimise the Bayes error under heteroscedasticity, as none of the datasets used has equal covariance matrices. With the exception of datasets (d), (j) and (l), where LDA achieves an equal Bayes error as the other heteroscedastic LDA approaches, LDA is outperformed by the GLD on all remaining datasets in terms of minimising the Bayes error. It will be noted that the other three heteroscedastic LDA approaches algorithms achieve a performance comparable to the GLD on all the datasets in terms of the Bayes error. However, R-HLD-1 and R-HLD-2 require a lot of trials ($1000$ in our experiments) in order to obtain the optimal parameters $s$ and $s_{1}$, $s_{2}$ respectively, while C-HLD requires a step size of $\Delta s=0.001$ which translates to $1001$ trials. Consequently, the training time for these algorithms far exceed that of the GLD, as can be seen in Table 4. For example, the gain in training time of the GLD over C-HLD, R-HLD-1 and R-HLD-2 is over $62$ folds for dataset (g), and about $20$ folds for dataset (l). Moreover, since C-HLD, R-HLD-1 and R-HLD-2 all require matrix inversions, performing a matrix inversion for each of the $1000$ trials can be a computationally demanding task especially for high-dimensional data, which have large covariance matrices. Instead, since the GLD follows a principled optimisation procedure, the number of matrix inversions required is far lower. For example, on dataset (f), which has a dimensionality of $47$, the GLD requires over $60$ times less time to train than the other heteroscedastic LDA approaches.

It is conceivable that the minimisation of the Bayes error would translate into a good performance in terms of the classification accuracy, if the normality assumption of LDA holds. For this reason, it can be seen in Table 3 that the GLD achieves the best classification accuracy on datasets (a) and (b), which are generated from known normal distributions. Thus, the proposed GLD algorithm is particularly suited for applications with datasets that tend to be normally distributed in each class e.g. in machine fault diagnosis, or accelerometer-based human activity recognition [Ojetola et al., 2015], as it also requires far less training time than the existing heteroscedastic LDA approaches.

However, for datasets (c) through to (l), the classes do not have any known normal distribution. Therefore, minimising the Bayes error under the normality assumption would not necessarily result in a classifier that has the best classification accuracy, even if the difference in covariance matrices has been accounted for. For this reason, it is not surprising that LDA achieves a superior classification accuracy than C-HLD, R-HLD-1, R-HLD-2 and the GLD on datasets (c) and (h) as can be seen in Table 3. However, by searching around the neighbourhood of the GLD, the Local Neighbourhood Search (LNS) algorithm is able to account for the non-normality and obtain a more robust classifier. Thus, the GLD, together with the LNS procedure, achieves a higher classification accuracy than all the LDA approaches on all the real-world datasets (i.e. (c) to (l)) with the exception of dataset (f) which has the GLD showing superior classification accuracy.

While the SVM outperforms the LDA approaches on half of the datasets, its training time can be rather long for large datasets. For instance, for dataset (d) which has $58000$ elements, the SVM takes about $1.3$ hours to train whereas the GLD with LNS, which achieves the best classification accuracy on this dataset, takes $43$ seconds to train, representing over $100$ fold savings in computational time over the SVM. Similar patterns can be seen in other datasets like (i), where the GLD with LNS achieves a superior classification accuracy with over $150$ times shorter training time than the SVM. This suggests that for such large datasets, the GLD with Local Neighbourhood Search is a low-complexity alternative to the SVM, as it requires far less computational time than the SVM.

We, however, make note of two weaknesses our proposed algorithms have. For the GLD, the procedure as described in Algorithm 1, may converge to a saddle point, instead of a local minimum. Even if it were to converge to a local minimum, there is no guarantee that is the global optimum solution due to the fact that the objective function $p_{e}$ is known to be non-convex Anderson & Bahadur [1962]. Also, since the Local Neighbourhood Search involves evaluating the misclassification rate on the training set for every perturbation, the procedure does not scale well with large amounts of training data. Because of this, it is important to have a good initial solution like the GLD, so that an early termination may be performed if there is no improvement after some number of iterations.

7 Conclusion

In this paper, we have presented the Gaussian Linear Discriminant (GLD), a novel and computationally efficient method for obtaining a linear discriminant for heteroscedastic Linear Discriminant Analysis (LDA) for the purpose of binary classification. Our algorithm minimises the Bayes error via an iterative optimisation procedure that uses Fisher’s Linear Discriminant as the initial solution. Moreover, the GLD does not require any parameter adjustments. We have also proposed a local neighbourhood search method by which a more robust linear classifier may be obtained for non-normal distributions. Our experimental results on two artificial and ten real world applications show that when the covariance matrices of the classes are unequal, LDA is unable to minimise the Bayes error. Thus, under heteroscedasticity, our proposed algorithm achieves superior classification accuracy to the LDA for normally distributed classes. While the proposed GLD algorithm compares favourably with other heteroscedastic LDA approaches, the GLD requires a far less training time. Moreover, the GLD, together with the LNS, has been shown to be particularly robust, comparing favourably with the SVM, but requiring far less training time on our datasets. Thus, for expert systems like machine fault diagnosis or human activity monitoring that require linear classification, the proposed algorithms provide a low-complexity, high-accuracy solution.

While this work has focused on linear classification, on-going work is focused on modifying the GLD procedure for the purpose of linear dimensionality reduction. Moreover, it is of particular interest to us to be able to derive the Bayes error for some known non-normal distributions. An alternative to this is to be able to obtain a kernel that implicitly transforms some data of a known non-normal distribution into a feature space where the classes are normally distributed. Finally, like all local search algorithms, the performance and complexity of the LNS procedure depends on the choice of the initial solution. Therefore, further work that explores the use initial solutions (including the heteroscedastic LDA approaches discussed) other than the GLD for the LNS procedure is being done.

Appendix A

Theorem 1**.**

Let $w_{0}^{+}$ and $w_{0}^{-}$ be the two distinct solutions of (34), then $w_{0}^{+}$ and $w_{0}^{-}$ cannot both satisfy (38) given that $\sigma_{1}\neq\sigma_{2}$.

Proof.

Let

[TABLE]

and let

[TABLE]

Then

[TABLE]

Suppose that $w_{0}^{+}$ satisfies (38), then

[TABLE]

i.e.,

[TABLE]

This implies that $\sigma_{1}^{2}/(\sigma_{1}^{2}-\sigma_{2}^{2})>\sigma_{2}^{2}/(\sigma_{1}^{2}-\sigma_{2}^{2})$ since $\beta$ is a positive scalar.

Consider now $w_{0}^{-}$ given as:

[TABLE]

Then

[TABLE]

In order for (38) to be satisfied, it can be shown, similar to (A.5), that

[TABLE]

which can be simplified to give $1\leq 0$. Since this conclusion is false, only $w_{0}^{+}$ satisfies (26). ∎