[1]Giulia Ferrandi

1]Department of Mathematics and Computer Science, TU Eindhoven, PO Box 513,
Eindhoven, 5600MB, The Netherlands

2]CEMAT and Department of Mathematics, Instituto Superior Técnico,
Universidade de Lisboa, Av. Rovisco Pais, Lisbon, 1049-001, Portugal

# On the trace ratio method and Fisher’s discriminant analysis for robust
multigroup classification

<EMAIL_ADDRESS>Igor V. Kravchenko<EMAIL_ADDRESS>Michiel E. Hochstenbach<EMAIL_ADDRESS>M. Rosário Oliveira
<EMAIL_ADDRESS>[  [

###### Abstract

We compare two different linear dimensionality reduction strategies for the
multigroup classification problem: the trace ratio method and Fisher’s
discriminant analysis. Recently, trace ratio optimization has gained in
popularity due to its computational efficiency, as well as the occasionally
better classification results. However, a statistical understanding is still
incomplete. We study and compare the properties of the two methods. Then, we
propose a robust version of the trace ratio method, to handle the presence of
outliers in the data. We reinterpret an asymptotic perturbation bound for the
solution to the trace ratio, in a contamination setting. Finally, we compare
the performance of the trace ratio method and Fisher’s discriminant analysis
on both synthetic and real datasets, using classical and robust estimators.

###### keywords:

Trace ratio method, Fisher’s discriminant analysis, multigroup classification,
linear dimensionality reduction, regularized MCD, robust statistics.

###### pacs:

[

MSC Classification]62H30, 15A18

##  1 Introduction

Classification is an important statistical task where an observation is
assigned to one of the non-overlapping known groups, based on the statistical
properties of the data characterizing these groups. Statistics, machine
learning, data science, and pattern recognition are some of the areas that use
this family of methods to solve practical problems. In this paper, we take a
closer look at two linear methods for multigroup classification: the trace
ratio method and Fisher’s discriminant analysis. By both theoretical analysis
and numerical simulations, we attempt to increase our understanding of the
properties of these methods. We also introduce and study a robust version of
the trace ratio method.

Suppose we have a population described by the random vector
$\mathbf{X}=(X_{1},\dots,X_{p})^{T}$, with mean
$\operatorname*{\mathbb{E}}(\mathbf{X})={\boldsymbol{\mu}}$ and covariance
$\operatorname*{{\rm Var}}(\mathbf{X})=\mathbf{\Sigma}$. We make the common
assumption that $\mathbf{\Sigma}$ is nonsingular. The population is divided
into $g$ groups (or classes), modeled by the class variable
$Y\in\\{1,\dots,g\\}$. Denote the expectation of
$\mathbf{X}_{j}=(\mathbf{X}\mid Y=j)$ by
$\operatorname*{\mathbb{E}}(\mathbf{X}_{j})={\boldsymbol{\mu}}_{j}$, and its
covariance matrix by $\operatorname*{{\rm
Var}}(\mathbf{X}_{j})=\mathbf{\Sigma}_{j}$. The probability of being in group
$j$ is $p_{j}=P(Y=j)$. Then

${\boldsymbol{\mu}}=\operatorname*{\mathbb{E}}(\mathbf{X})=\operatorname*{\mathbb{E}}(\operatorname*{\mathbb{E}}(\mathbf{X}\mid
Y))=\sum_{j=1}^{g}p_{j}\,{\boldsymbol{\mu}}_{j}.$ (1)

Given a new observation $\mathbf{x}\in\mathbb{R}^{p}$, we would like to assign
it to one of the $g$ classes, via a so-called _classification rule_ , that is
built on top of the available data. We are interested in methods where the
classification task includes a _linear_ dimensionality reduction step, where
the data are projected onto a lower-dimensional subspace.

In particular, we study the well-known Fisher’s discriminant analysis (FDA;
see, e.g., [1, Ch. 11]) and a more recent dimensionality reduction method, the
trace ratio method (TR; see, e.g., [2]). Both strategies aim at finding a
lower-dimensional representation of the data where the separation between the
groups is maximized in a certain sense. They both involve a measure of
distance between and within the groups. This is inspired by the decomposition
of the covariance matrix $\operatorname*{{\rm Var}}(\mathbf{X})$ as

$\operatorname*{{\rm
Var}}(\mathbf{X})=\operatorname*{\mathbb{E}}(\operatorname*{{\rm
Var}}(\mathbf{X}\mid Y))+\operatorname*{{\rm
Var}}(\operatorname*{\mathbb{E}}(\mathbf{X}\mid
Y))=\sum_{j=1}^{g}\,p_{j}\,\mathbf{\Sigma}_{j}+\sum_{j=1}^{g}\,p_{j}\,({\boldsymbol{\mu}}_{j}-{\boldsymbol{\mu}})\,({\boldsymbol{\mu}}_{j}-{\boldsymbol{\mu}})^{T},$
(2)

where $\mathbf{W}=\sum_{j=1}^{g}\,p_{j}\,\mathbf{\Sigma}_{j}$ is the _within
covariance_ matrix and
$\mathbf{B}=\sum_{j=1}^{g}\,p_{j}\,({\boldsymbol{\mu}}_{j}-{\boldsymbol{\mu}})\,({\boldsymbol{\mu}}_{j}-{\boldsymbol{\mu}})^{T}$
the _between covariance_ matrix. If projecting the data leads to a good
separation among the groups, a favorable effect on the performance of the
classifiers is expected, since the present noise may be damped or filtered. In
addition, the dimension reduction step tends to increase the efficiency of
subsequent steps and may help the analyst to interpret the classification
output.

The within and between covariance matrices can be estimated in different ways.
We start our discussion from the classical estimates in Section 2, reviewing
FDA and the trace ratio method. Next, we are interested in assessing the
classification performances of the two methods when robust estimates are used.
These are very useful when the data is contaminated, i.e., there is a certain
proportion $\varepsilon>0$ of outliers that do not follow the distribution of
$\mathbf{X}$. As we will observe in Section 5.2, while some nonlinear
classifiers seem to be less affected by outliers, linear classifiers may
benefit from the use of robust statistics.

FDA with robust estimates has already been studied mainly by exploring the
connection between FDA and linear discriminant analysis (LDA); see, e.g., [3].
Several methods to robustify discriminant analysis have been proposed. The
most straightforward proposals suggest employing robust estimators of group
means and robust covariance matrices [3, 4, 5, 6, 7]. These methods assume
that the number of observations per group, compared to the number of
variables, is sufficiently large to guarantee good estimates of the unknown
parameters. When we are in the opposite situation (which happens in many
practical problems), we say that our data is high-dimensional. The first
attempts to deal with this issue were done by Vanden Branden and Hubert [8],
who provide a robust classifier for high-dimensional data based on SIMCA [9],
and by Pires and Branco [10] who propose robust classifiers based on
projection-pursuit. Later on, Hoffmann and coauthors [11] suggest a sparse and
robust classification method based on partial least squares for binary
classification problems, while [12] presents a robust version of optimal
scoring, for both the non-sparse and sparse case.

In this paper, we propose a new robust TR method, obtained by exploiting MCD
robust estimates of the within and between covariance matrices. The method can
deal with high-dimensional data since it uses regularized MCD estimates (MRCD,
[13]); however, as we will see in Section 5.3, when a high number of
irrelevant variables has a negative impact on the performance of the
estimation methods, even on the robust ones.

We derive theoretical expressions showing that the between and within
theoretical covariance matrices, under a certain contamination scenario (cf.
Section 3.1), can be written as the sum of their non-contaminated counterparts
and a perturbation term. In the TR case, this formulation allows us to derive
an upper bound for the distance between the non-contaminated and the
contaminated solutions of the method.

We compare TR and FDA on synthetic and real datasets. Synthetic scenarios
consider cases where one method performs better than the others. In this case,
FDA and TR are used as classifiers and are compared by means of two different
criteria. While the first one is based on the performance of the associated
classification rules, the second criterion is related to the proximity between
the true solution of one method and the estimated one. Real datasets have been
chosen from the UCI [14] and KEEL [15] platforms, and they illustrate the
performance of FDA and TR as dimensionality reduction methods, used before the
construction of a classifier. In many of these datasets, FDA is as good as or
better than TR. Moreover, robust TR shows clear improvements compared to
classical TR in several datasets.

This paper is organized as follows. In Section 2 we review FDA, LDA, and the
trace ratio method. In Section 3 we derive the theoretical within and between
covariance matrices under a specific contamination scenario. For TR we derive
an upper bound quantifying the effect of an infinitesimal percentage of
contamination into the TR matrix of coefficients. In Section 4 we propose a
robust TR estimator. In Section 5 we carry out several simulation experiments
and in Section 6 we explore FDA and TR as dimensionality reduction strategies
for real datasets. The main conclusions are presented in Section 7.

We conclude this section by introducing some notations. We denote the $k\times
k$ identity matrix by $\mathbf{I}_{k}$. While $\|\cdot\|$ denotes the standard
Euclidean norm for vectors and the induced spectral norm for matrices (equal
to its largest singular value),
$\|\mathbf{x}\|_{\mathbf{A}}^{2}=\mathbf{x}^{T}\\!\mathbf{A}\mathbf{x}$ is
used for the $\mathbf{A}$-weighted norm for a symmetric positive definite
(SPD) matrix $\mathbf{A}$. Given the orthogonal bases of two subspaces (of
possibly unequal dimension) $\mathbf{V}$ and $\widehat{\mathbf{V}}$, the
(largest) angle between the two subspaces is defined as
$\angle(\mathbf{V},\widehat{\mathbf{V}})=\arcsin\|\mathbf{V}\mathbf{V}^{T}-\widehat{\mathbf{V}}\widehat{\mathbf{V}}^{T}\|$.
This quantity is well defined since
$\|\mathbf{V}\mathbf{V}^{T}-\widehat{\mathbf{V}}\widehat{\mathbf{V}}^{T}\|\in[0,\,1]$
(see, e.g., [16]).

##  2 FDA, LDA, and TR for multigroup classification

We present and compare FDA [17] with the more recent trace ratio optimization
method (see, e.g., [2] for a review). We discuss the two methods for the case
where classical estimates for the covariance matrices $\mathbf{B}$ and
$\mathbf{W}$ are used.

Let $\mathbf{x}_{j1},\dots,\mathbf{x}_{jn_{j}}$ represent $n_{j}$ realizations
of a random sample from $\mathbf{X}_{j}$; the total sample size is
$n=\sum_{j=1}^{g}n_{j}$. The sample mean and the sample covariance matrix of
group $j$ are
$\overline{\mathbf{x}}_{j}=n_{j}^{-1}\,\sum_{h=1}^{n_{j}}{\mathbf{x}_{jh}}$
and
$\mathbf{S}_{j}=(n_{j}-1)^{-1}\,\sum_{h=1}^{n_{j}}\,(\mathbf{x}_{jh}-\overline{\mathbf{x}}_{j})\,(\mathbf{x}_{jh}-\overline{\mathbf{x}}_{j})^{T}$,
respectively. The overall sample mean and covariance are then estimated by
$\overline{\mathbf{x}}=\sum_{j=1}^{g}\tfrac{n_{j}}{n}\,\overline{\mathbf{x}}_{j}$
and ${\mathbf{S}}_{{\rm
pooled}}=\sum_{j=1}^{g}\tfrac{n_{j}-1}{n-g}\,\mathbf{S}_{j}$. We define the
classical between-class scatter matrix, $\mathbf{B}$, and the classical
within-class scatter matrix, $\mathbf{W}$, as follows:

$\mathbf{B}=n^{-1}\,\sum_{j=1}^{g}n_{j}\,(\overline{\mathbf{x}}_{j}-\overline{\mathbf{x}})\,(\overline{\mathbf{x}}_{j}-\overline{\mathbf{x}})^{T},\qquad\mathbf{W}=n^{-1}\,\sum_{j=1}^{g}\sum_{h=1}^{n_{j}}\
(\mathbf{x}_{jh}-\overline{\mathbf{x}}_{j})\,(\mathbf{x}_{jh}-\overline{\mathbf{x}}_{j})^{T}.$
(3)

To keep the notation simple, we use $\mathbf{B}$ ($\mathbf{W})$ for both the
theoretical and classical (i.e., non-robust) estimated between (within)
covariance matrices. It is easy to check that
$\mathbf{W}=n^{-1}\,(n-g)\,{\mathbf{S}}_{{\rm pooled}}$. Let ${\rm
rank}(\cdot)$ indicate the rank of a matrix. Then $r={\rm
rank}(\mathbf{B})\leq\min\\{g-1,\,p\\}$ and ${\rm
rank}(\mathbf{W})\leq\min\\{n-g,\,p\\}$. Throughout the paper, we assume the
common situation that there are sufficiently many data points such that
$n-g\geq p$. Since $\mathbf{\Sigma}$ is assumed to be nonsingular, we also
assume that $\mathbf{W}$ is also nonsingular, and thus SPD. However, if some
variables are highly correlated, numerical singularity of $\mathbf{W}$ may
occur. One possible way to handle this situation is described in Section 6.

###  2.1 Fisher’s discriminant analysis

The original idea of Fisher [17] is to find the linear combination of the
data, $\mathbf{v}^{T}\mathbf{x}$, that leads to the highest separability
between the groups. By assuming that the population group covariance matrices
are equal and of full rank, we look for a projection that maximizes the
separation between the groups and minimizes the variability within the groups.
This corresponds to finding a projection vector, $\mathbf{v}_{\rm FDA}$, that
maximizes a separability index: $\mathbf{v}_{\rm
FDA}=\operatorname*{\arg\\!\max}_{\mathbf{v}\neq\mathbf{0}}\rho_{\rm
FDA}(\mathbf{v})$, where

$\rho_{\rm
FDA}(\mathbf{v})=\frac{\sum_{j=1}^{g}n_{j}\,(\mathbf{v}^{T}\overline{\mathbf{x}}_{j}-\mathbf{v}^{T}\overline{\mathbf{x}})^{2}}{\sum_{j=1}^{g}\sum_{h=1}^{n_{j}}(\mathbf{v}^{T}\mathbf{x}_{jh}-\mathbf{v}^{T}\overline{\mathbf{x}}_{j})^{2}}=\frac{\mathbf{v}^{T}\mathbf{B}\mathbf{v}}{\mathbf{v}^{T}\mathbf{W}\mathbf{v}}.$
(4)

This quantity is a Rayleigh quotient of the pair $(\mathbf{B},\mathbf{W})$.
Suppose we want to project our data onto a $k$-dimensional subspace spanned by
the columns of $\mathbf{V}\in\mathbb{R}^{p\times k}$, with $k\leq r$. FDA (4)
can be formalized as a generalized eigenvalue problem; see, e.g., [1, p. 623].

######  Proposition 1.

Let $\lambda_{1}\geq\dots\geq\lambda_{k}>0$ be the $k\leq r$ largest nonzero
eigenvalues of the generalized eigenvalue problem
$\mathbf{B}\mathbf{v}=\lambda\,\mathbf{W}\mathbf{v}$ and let
$\mathbf{v}_{1},\dots,\mathbf{v}_{k}$ be the associated eigenvectors (usually
scaled such that $\mathbf{v}_{i}^{T}{\mathbf{S}}_{{\rm
pooled}}\mathbf{v}_{i}=1$). Then the vector of coefficients $\mathbf{v}_{\rm
FDA}$ that maximizes (4) is given by $\mathbf{v}_{\rm FDA}=\mathbf{v}_{1}$.
The $i$th eigenvector is the solution to

$\mathbf{v}_{i}=\operatorname*{\arg\\!\max}_{\begin{subarray}{c}\mathbf{v}_{j}^{T}{\mathbf{S}}_{{\rm
pooled}}\mathbf{v}=0,\ j<i\\\ \mathbf{v}\neq\mathbf{0}\end{subarray}}\
\frac{\mathbf{v}^{T}\mathbf{B}\mathbf{v}}{\mathbf{v}^{T}\mathbf{W}\mathbf{v}}.$
(5)

The linear combination $\mathbf{v}_{1}^{T}\mathbf{x}$ is called _sample first
discriminant_ , and likewise the linear combination
$\mathbf{v}_{i}^{T}\mathbf{x}$, for $i\leq k$, is the sample $i$th
discriminant; Fisher’s discriminants generate a lower-dimensional space whose
span optimizes a sequence of separation criteria $\rho_{\rm FDA}(\cdot)$
defined in (4).

Problem (5) is equivalent to finding the first $k$ generalized eigenvectors
$\mathbf{V}_{\rm FDA}=[\mathbf{v}_{1},\dots,\mathbf{v}_{k}]$ of the pencil
$(\mathbf{B},\mathbf{W})$, with scaling $\mathbf{V}_{\rm
FDA}^{T}{\mathbf{S}}_{{\rm pooled}}\mathbf{V}_{\rm FDA}=\mathbf{I}_{k}$. Since
$\mathbf{W}$ and ${\mathbf{S}}_{{\rm pooled}}$ are the same matrix up to a
multiplicative factor, the generalized eigenvalue problem for FDA is
equivalent to finding the solution $\mathbf{V}_{\rm FDA}$ to the constrained
optimization problem (cf., e.g., [2])

$\max_{\mathbf{V}^{T}\mathbf{W}\,\mathbf{V}=\mathbf{I}_{k}}{\rm
tr}(\mathbf{V}^{T}\mathbf{B}\mathbf{V}),$ (6)

where ${\rm tr}(\cdot)$ indicates the trace of a matrix. Another equivalent
formulation is (cf. [18, Prop. 2.1.1]) $\max_{\text{\scriptsize
rank}(\mathbf{V})=k}{\rm
tr}((\mathbf{V}^{T}\mathbf{W}\,\mathbf{V})^{-1}\,\mathbf{V}^{T}\mathbf{B}\mathbf{V})$,
where the objective function may be viewed as a _generalized Rayleigh
quotient_.

###  2.2 Relation between linear and Fisher’s discriminant analysis

When we assume that all the group covariance matrices are equal, i.e.,
$\mathbf{\Sigma}_{j}\equiv\mathbf{\Sigma}$, FDA has a strong connection with
LDA, which assumes that, within groups, the random variables have a
multivariate normal distribution with equal covariance matrices. The
associated classification rule is derived such that it minimizes the _total
probability of misclassification_ , which is the probability of assigning a
given data point to an incorrect class (cf. [1, Result 11.5]). According to
LDA, an observation $\mathbf{x}$ is assigned to class $j$ if (see, e.g., [1,
p. 630])

$j=\operatorname*{\arg\\!\min}_{i}\|\mathbf{x}-{\boldsymbol{\mu}}_{i}\|_{\mathbf{\Sigma}^{-1}}^{2}-2\,\log
p_{i},$ (7)

and $\|\mathbf{x}-{\boldsymbol{\mu}}_{i}\|_{\mathbf{\Sigma}^{-1}}$ is the
Mahalanobis distance between $\mathbf{x}$ and the group mean
${\boldsymbol{\mu}}_{i}$. LDA is also the classifier that maximizes the total
probability of correctly assigning an observation (called the accuracy) among
all possible classifiers for multivariate normal mixtures with equal
covariance matrices.

When $k=r$, LDA with the classical estimates $\overline{\mathbf{x}}_{i}$ and
${\mathbf{S}}_{{\rm pooled}}$ coincides with projecting a data point
$\mathbf{x}$ onto $\mathbf{V}_{\rm FDA}\in\mathbb{R}^{p\times r}$ and
classifying it to the group with the closest projected mean, according to the
Euclidean distance. This result is reported in Proposition 2, which is a
slightly different formulation of [1, Result 11.6]. Recall that $\mathbf{B}$
has $r\leq\min{\\{g-1,p\\}}$ positive eigenvalues.

######  Proposition 2.

Let $[\mathbf{V}_{\rm FDA}\;\mathbf{V}_{\perp}]$ contain the
${\mathbf{S}}_{{\rm pooled}}$-orthogonal generalized eigenvectors of the
pencil $(\mathbf{B},\mathbf{W})$, with $\mathbf{V}_{\rm
FDA}=[\mathbf{v}_{1},\dots,\mathbf{v}_{r}]$ and
$\mathbf{V}_{\perp}^{T}{\mathbf{S}}_{{\rm pooled}}\mathbf{V}_{\rm
FDA}=\mathbf{0}$. Then only the terms $\|\mathbf{V}_{\rm
FDA}^{T}(\mathbf{x}-\overline{\mathbf{x}}_{i})\|^{2}$ contribute to the
classification. In mathematical terms,

$\operatorname*{\arg\\!\min}_{i}\|\mathbf{x}-\overline{\mathbf{x}}_{i}\|_{{\mathbf{S}}_{{\rm
pooled}}^{-1}}^{2}-2\,\log\tfrac{n_{i}}{n}=\operatorname*{\arg\\!\min}_{i}\|\mathbf{V}_{\rm
FDA}^{T}(\mathbf{x}-\overline{\mathbf{x}}_{i})\|^{2}-2\,\log\tfrac{n_{i}}{n}.$

######  Proof 2.3.

Since ${\rm rank}(\mathbf{B})=r$, we have
$\mathbf{B}\mathbf{V}_{\perp}=\mathbf{0}$. In particular,
$\mathbf{V}_{\perp}^{T}(\overline{\mathbf{x}}_{i}-\overline{\mathbf{x}})=\bf
0$ for all $i$. This can be deduced from the fact that
$0=\text{tr}(\mathbf{V}^{T}_{\perp}\mathbf{B}\mathbf{V}_{\perp})=\sum_{i=1}^{g}\frac{n_{i}}{n}\,\|\mathbf{V}^{T}_{\perp}(\overline{\mathbf{x}}_{i}-\overline{\mathbf{x}})\|^{2}$.
Since ${\mathbf{S}}_{{\rm pooled}}^{-1}=\mathbf{V}_{\rm FDA}\mathbf{V}_{\rm
FDA}^{T}+\mathbf{V}_{\perp}\mathbf{V}_{\perp}^{T}$, we have

$\displaystyle\|\mathbf{x}-\overline{\mathbf{x}}_{i}\|_{{\mathbf{S}}_{{\rm
pooled}}^{-1}}^{2}$ $\displaystyle=\|\mathbf{V}_{\rm
FDA}^{T}(\mathbf{x}-\overline{\mathbf{x}}_{i})\|^{2}+\|\mathbf{V}_{\perp}^{T}\,[\,(\mathbf{x}-\overline{\mathbf{x}})+(\overline{\mathbf{x}}-\overline{\mathbf{x}}_{i})\,]\,\|^{2}$
$\displaystyle=\|\mathbf{V}_{\rm
FDA}^{T}(\mathbf{x}-\overline{\mathbf{x}}_{i})\|^{2}+\|\mathbf{V}_{\perp}^{T}(\mathbf{x}-\overline{\mathbf{x}})\|^{2},$

where the second term in the right-hand side is common to all groups.

From this last proposition, it is easy to see that Rule (7) is equivalent to
classifying a data point $\mathbf{x}$ to the class with the closest projected
mean $\mathbf{V}_{\rm FDA}^{T}\overline{\mathbf{x}}_{j}$, up to an additive
weight related to the class proportions.

For a reduced dimension $k<r$, there is still a connection between Rule (7)
and the classification via the Euclidean distance. This requires an estimation
of the group means and the common covariance matrix by assuming that the group
means lie in a $k$-dimensional space. The result is the following.

######  Proposition 2.4.

[19, Prop. 1] Consider a family of $g$ multivariate normal distributions, with
mean vectors ${\boldsymbol{\mu}}_{i}$ and common covariance matrix
$\mathbf{\Sigma}$. Let $\mathbf{V}_{\rm FDA}$ be the $p\times k$ solution to
the generalized eigenvalue problem $(\mathbf{B},\mathbf{W})$ (6), scaled such
that $\mathbf{V}_{\rm FDA}^{T}{\mathbf{S}}_{{\rm pooled}}\mathbf{V}_{\rm
FDA}=\mathbf{I}_{k}$. Let $\widehat{\boldsymbol{\mu}}_{i}$ and
$\widehat{\mathbf{\Sigma}}$ be the maximum likelihood estimators of
${\boldsymbol{\mu}}_{i}$ and $\mathbf{\Sigma}$, with the additional constraint
that the centroids ${\boldsymbol{\mu}}_{i}$ lie in a $k$-dimensional space.
Then $\widehat{\boldsymbol{\mu}}_{i}={\mathbf{S}}_{{\rm
pooled}}\mathbf{V}_{\rm FDA}\mathbf{V}_{\rm
FDA}^{T}\,\overline{\mathbf{x}}_{i}+(\mathbf{I}_{p}-{\mathbf{S}}_{{\rm
pooled}}\mathbf{V}_{\rm FDA}\mathbf{V}_{\rm FDA}^{T})\,\overline{\mathbf{x}}$
and $\widehat{\mathbf{\Sigma}}={\mathbf{S}}_{{\rm
pooled}}+(\mathbf{I}_{p}-{\mathbf{S}}_{{\rm pooled}}\mathbf{V}_{\rm
FDA}\mathbf{V}_{\rm FDA}^{T})\,\mathbf{B}\,(\mathbf{I}_{p}-{\mathbf{S}}_{{\rm
pooled}}\mathbf{V}_{\rm FDA}\mathbf{V}_{\rm FDA}^{T})^{T}.$ Moreover,

$\|\mathbf{x}-\widehat{\boldsymbol{\mu}}_{i}\|_{\widehat{\mathbf{\Sigma}}^{-1}}^{2}=\|\mathbf{V}_{\rm
FDA}^{T}\mathbf{x}-\mathbf{V}_{\rm
FDA}^{T}\overline{\mathbf{x}}_{i}\|^{2}+(\mathbf{x}-\overline{\mathbf{x}})^{T}\,(\widehat{\mathbf{\Sigma}}^{-1}-\mathbf{V}_{\rm
FDA}\mathbf{V}_{\rm FDA}^{T})\,(\mathbf{x}-\overline{\mathbf{x}}),$

where the rightmost term is common to all groups.

We can conclude that classifying a projected data point according to the
Euclidean distance from the projected means is equivalent to classifying it in
the original space according to the Mahalanobis distance, where
${\boldsymbol{\mu}}_{i}$ and $\mathbf{\Sigma}$ are estimated appropriately
based on the rank of $\mathbf{V}_{\rm FDA}$, i.e., by
$\overline{\mathbf{x}}_{i}$ and ${\mathbf{S}}_{{\rm pooled}}$ when ${\rm
rank}(\mathbf{V}_{\rm FDA})=r$ (Proposition 2), and by
$\widehat{\boldsymbol{\mu}}_{i}$ and $\widehat{\mathbf{\Sigma}}$ when ${\rm
rank}(\mathbf{V}_{\rm FDA})=k<r$ (Proposition 2.4). In other words, regardless
of the rank of $\mathbf{V}_{\rm FDA}$, $\mathbf{x}$ is assigned to class $j$
if

$j=\operatorname*{\arg\\!\min}_{i}\|\mathbf{V}_{\rm
FDA}^{T}(\mathbf{x}-\overline{\mathbf{x}}_{i})\|^{2}-2\,\log\tfrac{n_{i}}{n},$
(8)

where $\frac{n_{i}}{n}$ estimates the prior probability $p_{i}$ of group $i$.

We remark that the classification rule (8) does not require any distributional
assumption. In view of this, in general, it is not a Bayesian classification
rule and thus does not minimize the total probability of misclassification.
Nevertheless, FDA with $k=r$ is optimal under the homoscedastic assumption and
with multivariate normally distributed data; in this case FDA and LDA
coincide.

###  2.3 Trace ratio optimization

Recently, a different alternative to FDA has been studied (see, e.g., [20, Ch.
10] and [21, 2, 22]). The essence is that instead of the trace of a ratio (see
Section 2.1), a ratio of traces is considered. If we consider the projected
data $\mathbf{V}^{T}\mathbf{x}$, another possible generalization of $\rho_{\rm
FDA}$ (4) comes from considering the Euclidean distance between the projected
vectors. The between covariance matrix of the projected data can be written as
a function of the distances between $\mathbf{V}^{T}\overline{\mathbf{x}}_{j}$
and $\mathbf{V}^{T}\overline{\mathbf{x}}$ (see, e.g., [2]) in the following
way:

$\displaystyle\sum_{j=1}^{g}n_{j}\,\|\mathbf{V}^{T}\overline{\mathbf{x}}_{j}-\mathbf{V}^{T}\overline{\mathbf{x}}\|^{2}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{g}n_{j}\,(\mathbf{V}^{T}\overline{\mathbf{x}}_{j}-\mathbf{V}^{T}\overline{\mathbf{x}})^{T}\,(\mathbf{V}^{T}\overline{\mathbf{x}}_{j}-\mathbf{V}^{T}\overline{\mathbf{x}})$
$\displaystyle=$ $\displaystyle\sum_{j=1}^{g}{\rm
tr}(n_{j}\,\mathbf{V}^{T}(\overline{\mathbf{x}}_{j}-\overline{\mathbf{x}})\,(\overline{\mathbf{x}}_{j}-\overline{\mathbf{x}})^{T}\,\mathbf{V})\
=\ n\ {\rm tr}(\mathbf{V}^{T}\mathbf{B}\mathbf{V}).$

Similarly, for the common variability within the projected groups we have

$\sum_{j=1}^{g}\,\sum_{h=1}^{n_{j}}\,\|\mathbf{V}^{T}\mathbf{x}_{jh}-\mathbf{V}^{T}\overline{\mathbf{x}}_{j}\|^{2}=n\
{\rm tr}(\mathbf{V}^{T}\mathbf{W}\,\mathbf{V}).$

The new optimization problem aims at finding the orthogonal projection
$\mathbf{V}_{\rm TR}\in\mathbb{R}^{p\times k}$ that maximizes the trace ratio
(TR):

$\mathbf{V}_{\rm
TR}=\operatorname*{\arg\\!\max}_{\mathbf{V}^{T}\mathbf{V}=\mathbf{I}_{k}}\
\rho(\mathbf{V})=\frac{{\rm tr}(\mathbf{V}^{T}\mathbf{B}\mathbf{V})}{{\rm
tr}(\mathbf{V}^{T}\mathbf{W}\,\mathbf{V})}.$ (9)

Other variants of the problem (9) have been considered in [20, Ch. 10] and
[21, 2, 22]; these include a different set of constraints on $\mathbf{V}$ or
consider matrices different from $\mathbf{W}$ in the denominator of $\rho$.

In the one-dimensional case, when $k=1$, problems (4), (6), and (9) coincide,
up to the scaling of the solution vector. For $k>1$, the trace ratio problem
(9) does not have an explicit solution and thus relies on iterative
approaches; see, e.g., [23, 2, 24, 25]. Still, the authors of [22] report that
the TR problem reduces to a least squares problem for the high-dimensional
case $p\geq n$, and provide a closed form solution.

The existence of $\mathbf{V}_{\rm TR}$ is guaranteed under a milder hypothesis
on the within covariance matrix. We hereby report [2, Prop. 3.2] for
completeness.

######  Proposition 2.5.

Let $\mathbf{B}$ be symmetric and $\mathbf{W}$ symmetric positive
semidefinite, with ${\rm rank}(\mathbf{W})\geq p-k+1$. Then (9) admits a
finite maximum.

In particular, Proposition 2.5 implies that for a finite maximum, $k$ may be
larger than $r={\rm rank}(\mathbf{B})\leq\min\\{g-1,p\\}$, and $\mathbf{W}$
does not need to be positive definite. Several authors (see, e.g., [2])
mention the constraint $\mathbf{V}^{T}\mathbf{C}\mathbf{V}=\mathbf{I}_{k}$ for
a certain SPD matrix $\mathbf{C}$. Since the eigenvectors of the pair
$(\mathbf{B},\mathbf{W})$ are both $\mathbf{B}$-orthogonal and
$\mathbf{W}$-orthogonal, we remark that for the convex combination
$\mathbf{C}=(1-\gamma_{1})\,\mathbf{B}+\gamma_{1}\mathbf{W}$, with
$0<\gamma_{1}\leq 1$ such that $\mathbf{C}$ is SPD, the TR method is
equivalent to a generalized eigenvalue problem (6), which renders FDA.

Algorithm 1 shows the pseudocode for a natural iterative method (see, e.g.,
[2]) to find the solution to (9).

Algorithm 1: Iterative method to solve the trace ratio problem (9)

Input: $\mathbf{B}$, $\mathbf{W}\in\mathbb{R}^{p\times p}$, dimension $k$,
initial $(p\times k)$ matrix $\mathbf{V}$ with orthogonal columns.
Output: Solution $\mathbf{V}\in\mathbb{R}^{p\times k}$ to (9).
1: If $\mathbf{V}$ is not provided, choose a random $p\times k$ matrix with
orthogonal columns while not converged (as determined by $\rho(\mathbf{V})$ or
$\mathbf{V}$) 2: M $\rho=\rho(\mathbf{V})=\frac{{\rm
tr}(\mathbf{V}^{T}\mathbf{B}\mathbf{V})}{{\rm
tr}(\mathbf{V}^{T}\mathbf{W}\,\mathbf{V})}$ 3: M Determine (orthogonal)
$\mathbf{V}$ as largest $k$ eigenvectors of $\mathbf{B}-\rho\mathbf{W}$ end

We would like to point out that we may also start Algorithm 1 with the
orthogonalized FDA solution $\mathbf{V}=\mathbf{V}_{\rm
FDA}\in\mathbb{R}^{p\times k}$.

The global convergence of Algorithm 1 is stated in, e.g., [26, Thm. 5.1].
There the uniqueness of the solution $\mathbf{V}_{\rm TR}$ is also discussed,
and it is related to line 3 of Algorithm 1: if there is a gap between the
largest $k$ eigenvalues of $\mathbf{B}-\rho(\mathbf{V}_{\rm TR})\,\mathbf{W}$
and the $p-k$ smallest ones, then the solution is uniquely described by the
set $\\{\mathbf{V}_{\rm TR}\mathbf{Q}\,:\,\mathbf{Q}\in\mathbb{R}^{k\times
k},\,\mathbf{Q}^{T}\mathbf{Q}=\mathbf{I}_{k}\\}$.

The stopping criterion may be based on the (relative) error between two
consecutive $\rho$-values, or the angle between two consecutive subspaces
spanned by the columns of $\mathbf{V}$; in our implementation, we use the
first approach.

Unlike FDA, TR does not come with a natural classifier. In [2], the method of
the $k$-nearest neighbors is used to classify new data points in the projected
space. We propose to classify a data point to the group with the closest
projected sample mean vectors, according to the Euclidean distance and the
class proportions:

$j=\operatorname*{\arg\\!\min}_{i}\|\mathbf{V}_{\rm
TR}^{T}(\mathbf{x}-\overline{\mathbf{x}}_{i})\|^{2}-2\,\log\tfrac{n_{i}}{n},$

as it is done in FDA, although we do not have a connection with the
Mahalanobis distance in the original space and the Euclidean distance in the
projected space.

###  2.4 FDA, TR, and the reduced dimension

We now study the solution to FDA and TR, and some relevant quantities, as a
function of $k$, the selected reduced dimension. We recall that for FDA,
$1\leq k\leq r={\rm rank}(\mathbf{B})$; since $\mathbf{W}$ is assumed
nonsingular, the reduced dimension for TR is not restricted, thus $1\leq k\leq
p={\rm rank}(\mathbf{W})$.

The FDA solution $\mathbf{V}_{\rm FDA}^{(k)}$ corresponds to the first $k$
generalized eigenvectors of $(\mathbf{B},\mathbf{W})$; therefore it holds
$\text{span}(\mathbf{V}_{\rm FDA}^{(k)})\subset\text{span}(\mathbf{V}_{\rm
FDA}^{(k+1)})$ for $k=1,\dots,r-1$. If we are interested in choosing the best
value for $k$, it is enough to compute the eigenvalue decomposition of
$(\mathbf{B},\mathbf{W})$ once and retain the desired number of eigenvectors.
Another consequence is the following result. For simplicity, we choose
$\mathbf{V}_{\rm FDA}^{(k)}$ to have $\mathbf{W}$-orthogonal eigenvectors
(instead of ${\mathbf{S}}_{{\rm pooled}}$-orthogonal).

######  Proposition 2.6.

Let the eigenvectors of $(\mathbf{B},\mathbf{W})$ be scaled
$\mathbf{W}$-orthogonal. Then ${\rm tr}((\mathbf{V}_{\rm
FDA}^{(k)})^{T}\,\mathbf{B}\mathbf{V}_{\rm
FDA}^{(k)})=\lambda_{1}+\cdots+\lambda_{k}$, ${\rm tr}((\mathbf{V}_{\rm
FDA}^{(k)})^{T}\,\mathbf{W}\,\mathbf{V}_{\rm FDA}^{(k)})={\rm
tr}(\mathbf{I}_{k})=k$, and therefore both quantities are increasing as a
function of $k$.

The trace ratio $\rho(\mathbf{V}_{\rm FDA}^{(k)})$ for the FDA subspace is
equal to $(\lambda_{1}+\cdots+\lambda_{k})/k$, and therefore is a non-
increasing function of $k$.

The quantity ${\rm tr}((\mathbf{V}_{\rm
FDA}^{(k)})^{T}\,\mathbf{B}\mathbf{V}_{\rm FDA}^{(k)})$ is mentioned in [1,
(11-68)] as a theoretical criterion to choose $k$.

Now we turn our attention to the TR problem. Algorithm 1 needs the value of
$k$ in advance; the processes and outputs for different values of $k$ are
generally unrelated, and therefore we expect $\text{span}(\mathbf{V}_{\rm
TR}^{(k)})\not\subset\text{span}(\mathbf{V}_{\rm TR}^{(k+1)})$. The following
example illustrates this statement.

######  Example 2.7.

By applying an orthogonal transformation, one may choose the within scatter
matrix $\mathbf{W}$ to be diagonal without loss of generality. Therefore, we
generate $1000$ random problems, where $\mathbf{W}={\rm diag}(1,2,\dots,20)$
is fixed, while the between scatter matrix is given by
$\mathbf{B}=\mathbf{A}\mathbf{A}^{T}$, where $\mathbf{A}$ is a $20\times 10$
matrix with normally distributed entries, after which $\mathbf{A}$ is column
centered. The columns represent $10$ groups. For each problem, we run TR with
$k=1,2,\dots,8$ and compute the angle between two consecutive solutions,
$\mathbf{V}_{\rm TR}^{(k)}$ and $\mathbf{V}_{\rm TR}^{(k+1)}$. Figure 1(a)
reports the sorted angles $\angle(\mathbf{V}_{\rm TR}^{(2)},\,\mathbf{V}_{\rm
TR}^{(3)})$ for each problem. Their values range from rather small to
approximately $\pi/5$, meaning that the two subspaces can be rather different.
In Figure 1(b), we carry out similar experiments, but now we display the
median $\angle(\mathbf{V}_{\rm TR}^{(k)},\,\mathbf{V}_{\rm TR}^{(k+1)})$ over
1000 experiments for various values of $k$.

(a) (b)

Figure 1: (a): Subspace angle between two-dimensional and three-dimensional TR
spaces $\mathbf{V}_{\rm TR}^{(2)}$ and $\mathbf{V}_{\rm TR}^{(3)}$. (b):
Median angle $\angle(\mathbf{V}_{\rm TR}^{(k)},\,\mathbf{V}_{\rm TR}^{(k+1)})$
over 1000 random problems, for $k=1,\dots,8$.

Since $\mathbf{V}_{\rm TR}^{(k)}\not\subset\mathbf{V}_{\rm TR}^{(k+1)}$, the
quantities ${\rm tr}((\mathbf{V}_{\rm
TR}^{(k)})^{T}\,\mathbf{B}\mathbf{V}_{\rm TR}^{(k)})$ and ${\rm
tr}((\mathbf{V}_{\rm TR}^{(k)})^{T}\,\mathbf{W}\mathbf{V}_{\rm TR}^{(k)})$ do
not have an easy expression and interpretation. Nevertheless, we can show that
the optimal trace ratio $\rho^{(k)}=\rho(\mathbf{V}_{\rm TR}^{(k)})$ is a non-
increasing function of $k$. To the best of our knowledge, this result is new.

######  Proposition 2.8.

Let $\rho^{(k)}$ be the optimal function value of the TR, for a
$k$-dimensional subspace of $\mathbb{R}^{p}$ and the pencil
$(\mathbf{B},\,\mathbf{W})$, where $\mathbf{W}$ is SPD. Then $\rho^{(k)}$ is a
non-increasing function of $k$, for $k=1,\dots,p$.

######  Proof 2.9.

Let $\mathbf{V}_{\rm TR}^{(k)}\in\mathbb{R}^{p\times k}$ be the solution to TR
with optimal value $\rho^{(k)}$, for $k=1,\dots p-1$. Because of the linearity
of the trace, we can write

$\rho^{(k+1)}-\rho^{(k)}=\frac{{\rm tr}((\mathbf{V}_{\rm
TR}^{(k+1)})^{T}\,(\mathbf{B}-\rho^{(k)}\mathbf{W})\,\mathbf{V}_{\rm
TR}^{(k+1)})}{{\rm tr}((\mathbf{V}_{\rm
TR}^{(k+1)})^{T}\,\mathbf{W}\mathbf{V}_{\rm TR}^{(k+1)})},$

where the denominator is positive. Denote by $\lambda_{i}(\mathbf{A})$ the
$i$th largest eigenvalue of a matrix $\mathbf{A}$. We recall that

$0={\rm tr}((\mathbf{V}_{\rm
TR}^{(k)})^{T}\,(\mathbf{B}-\rho^{(k)}\mathbf{W})\,\mathbf{V}_{\rm
TR}^{(k)})=\sum_{i=1}^{k}\lambda_{i}(\mathbf{B}-\rho^{(k)}\mathbf{W}).$

The first equality comes from the first-order necessary condition for the
maximum (see, e.g., [2]), while the second follows from the fact that
$\mathbf{V}_{\rm TR}^{(k)}$ contains the $k$ largest eigenvectors of
$\mathbf{B}-\rho^{(k)}\mathbf{W}$. Since the sum of the $k$ largest
eigenvalues of $\mathbf{B}-\rho^{(k)}\mathbf{W}$ is zero, the $(k+1)$st
largest eigenvalue must be nonpositive:
$\lambda_{k+1}(\mathbf{B}-\rho^{(k)}\mathbf{W})\leq 0$. Then

${\rm tr}((\mathbf{V}_{\rm
TR}^{(k+1)})^{T}(\mathbf{B}-\rho^{(k)}\mathbf{W})\,\mathbf{V}_{\rm
TR}^{(k+1)})\leq\max_{\mathbf{V}^{T}\mathbf{V}=\mathbf{I}_{k+1}}{\rm
tr}(\mathbf{V}^{T}(\mathbf{B}-\rho^{(k)}\mathbf{W})\,\mathbf{V})=\lambda_{k+1}(\mathbf{B}-\rho^{(k)}\mathbf{W})\leq
0.$

This proves that $\rho^{(k+1)}-\rho^{(k)}\leq 0$.

The monotonicity from this proposition does not need to be strict, in the case
of a multiple generalized eigenvalue. An easy example of this is
$\mathbf{B}={\rm diag}(2,4)$ and $\mathbf{W}={\rm diag}(1,2)$, when
$\rho^{(1)}=\rho^{(2)}=2$. We also note that, for this result, the hypothesis
that $\mathbf{W}$ is SPD can be weakened; one only needs that ${\rm
tr}((\mathbf{V}_{\rm TR}^{(k)})^{T}\,\mathbf{W}\mathbf{V}_{\rm TR}^{(k)})>0$
for all $k$.

By rather extensive experimentation, we have reason to believe that ${\rm
tr}((\mathbf{V}_{\rm TR}^{(k)})^{T}\mathbf{B}\mathbf{V}_{\rm TR}^{(k)})$ is a
non-decreasing function of $k$, just as for FDA. A proof or counterexample is
currently lacking.

######  Conjecture 2.10.

Let $\mathbf{B}$ be symmetric positive semidefinite and $\mathbf{W}$ be SPD,
and let $\mathbf{V}_{\rm TR}^{(k)}$ be the quantity determined by the trace
ratio method. The quantity ${\rm tr}((\mathbf{V}_{\rm
TR}^{(k)})^{T}\mathbf{B}\mathbf{V}_{\rm TR}^{(k)})$ is a non-decreasing
function of $k$.

If this conjecture is true, the following related result follows as well.

######  Corollary 2.11.

If Conjecture 2.10 holds, then ${\rm tr}((\mathbf{V}_{\rm
TR}^{(k)})^{T}\,\mathbf{W}\mathbf{V}_{\rm TR}^{(k)})$ is also non-decreasing
in $k$.

######  Proof 2.12.

This follows from Proposition 2.8. In fact, since
$\rho^{(k+1)}\leq\rho^{(k)}$,

$\frac{{\rm tr}((\mathbf{V}_{\rm TR}^{(k+1)})^{T}\mathbf{B}\mathbf{V}_{\rm
TR}^{(k+1)})}{{\rm tr}((\mathbf{V}_{\rm
TR}^{(k)})^{T}\mathbf{B}\mathbf{V}_{\rm TR}^{(k)})}\leq\frac{{\rm
tr}((\mathbf{V}_{\rm TR}^{(k+1)})^{T}\,\mathbf{W}\mathbf{V}_{\rm
TR}^{(k+1)})}{{\rm tr}((\mathbf{V}_{\rm
TR}^{(k)})^{T}\,\mathbf{W}\mathbf{V}_{\rm TR}^{(k)})}.$

Therefore, if Conjecture 2.10 holds, since the trace ratio is decreasing, the
quantity ${\rm tr}((\mathbf{V}_{\rm TR}^{(k)})^{T}\,\mathbf{W}\mathbf{V}_{\rm
TR}^{(k)})$ increases faster than ${\rm tr}((\mathbf{V}_{\rm
TR}^{(k)})^{T}\mathbf{B}\mathbf{V}_{\rm TR}^{(k)})$ as $k$ increases.

From the descriptions of FDA and TR methods, it is clear that any pencil
$(\mathbf{B},\mathbf{W})$, with $\mathbf{B}$ symmetric and $\mathbf{W}$ SPD,
estimating the between and the within covariance matrices, may be exploited in
the contexts of FDA and TR. This motivates the use of robust estimators in
Section 4 and of the true covariance matrices in the simulation study (Section
5).

Although FDA seems to have a stronger statistical motivation than TR, in
Section 5 we illustrate some situations in which TR performs better than FDA
in the classification task.

##  3 Effect of contamination in the within and between covariance matrices

In the following sections, we express the contaminated theoretical within and
between covariance matrices as the original matrices plus a perturbation term.
These perturbation terms allow us to derive an upper limit of a proximity
measure between the TR solution when no contamination is present, and their
equivalent under contamination.

###  3.1 Contaminated theoretical between and within covariance matrices

We model the presence of atypical data as a mixture of two random vectors: the
majority of the data comes from a certain (central, clean, or non-
contaminated) mixture model $\mathbf{X}$, and, in each group, there is a
probability $\varepsilon$ for getting data generated by a different
distribution, characterizing the atypical scenario. Let
${}^{c}{\mathbf{X}}_{j}$ be the contaminating random variable, where
$\operatorname*{\mathbb{E}}({}^{c}{\mathbf{X}}_{j})={}^{c}{{\boldsymbol{\mu}}}_{j}$
and $\operatorname*{{\rm
Var}}({}^{c}{\mathbf{X}}_{j})={}^{c}{\mathbf{\Sigma}}_{j}$. Let
$\mathbf{Z}_{j}$ be the contaminated variable that describes the distribution
of group $j$:

$\mathbf{Z}_{j}=(1-\Omega_{j})\,\mathbf{X}_{j}+\Omega_{j}\,{}^{c}{\mathbf{X}}_{j},$

with $\Omega_{j}\sim{\rm Bernoulli}(\varepsilon)$. In addition, we assume that
$\Omega_{j}$, $\mathbf{X}_{j}$, and ${}^{c}{\mathbf{X}}_{j}$ are independent
random variables for any value of $j$. Then, we can derive the mean vector and
covariance matrix of $\mathbf{Z}_{j}$:

$\displaystyle\operatorname*{\mathbb{E}}(\mathbf{Z}_{j})$
$\displaystyle=(1-\varepsilon)\,{\boldsymbol{\mu}}_{j}+\varepsilon\
{}^{c}{{\boldsymbol{\mu}}}_{j}={\boldsymbol{\mu}}_{j}+\varepsilon\,({}^{c}{{\boldsymbol{\mu}}}_{j}-{\boldsymbol{\mu}}_{j})={\boldsymbol{\mu}}_{j}+\varepsilon\,\Delta_{j},\quad{\rm
where}\quad\Delta_{j}={}^{c}{{\boldsymbol{\mu}}}_{j}-{\boldsymbol{\mu}}_{j}.$
(10)

The covariance matrix of $\mathbf{Z}_{j}$ is given by

$\displaystyle\operatorname*{{\rm Var}}(\mathbf{Z}_{j})$
$\displaystyle=\operatorname*{\mathbb{E}}(\mathbf{Z}_{j}\mathbf{Z}_{j}^{T})-\operatorname*{\mathbb{E}}(\mathbf{Z}_{j})\operatorname*{\mathbb{E}}(\mathbf{Z}_{j})^{T}$
$\displaystyle=\varepsilon\,\operatorname*{\mathbb{E}}({}^{c}{\mathbf{X}}_{j}({}^{c}{\mathbf{X}}_{j})^{T})+(1-\varepsilon)\,\operatorname*{\mathbb{E}}(\mathbf{X}_{j}(\mathbf{X}_{j})^{T})$
$\displaystyle\phantom{MMM}-\varepsilon^{2}\
{}^{c}{{\boldsymbol{\mu}}}_{j}\,({}^{c}{{\boldsymbol{\mu}}}_{j})^{T}-(1-\varepsilon)^{2}\,{\boldsymbol{\mu}}_{j}\,{\boldsymbol{\mu}}_{j}^{T}-\varepsilon\,(1-\varepsilon)\,\left({}^{c}{{\boldsymbol{\mu}}}_{j}\,{\boldsymbol{\mu}}_{j}^{T}+{\boldsymbol{\mu}}_{j}\,({}^{c}{{\boldsymbol{\mu}}}_{j})^{T}\right)$
$\displaystyle=\varepsilon\
{}^{c}{\mathbf{\Sigma}}_{j}+(1-\varepsilon)\,\mathbf{\Sigma}_{j}+\varepsilon\,(1-\varepsilon)\,\left({}^{c}{{\boldsymbol{\mu}}}_{j}\,({}^{c}{{\boldsymbol{\mu}}}_{j})^{T}+{\boldsymbol{\mu}}_{j}\,{\boldsymbol{\mu}}_{j}^{T}-{}^{c}{{\boldsymbol{\mu}}}_{j}{\boldsymbol{\mu}}_{j}^{T}-{\boldsymbol{\mu}}_{j}({}^{c}{{\boldsymbol{\mu}}}_{j})^{T}\right)$
$\displaystyle=\varepsilon\
{}^{c}{\mathbf{\Sigma}}_{j}+(1-\varepsilon)\,\mathbf{\Sigma}_{j}+\varepsilon\,(1-\varepsilon)\,\Delta_{j}\,\Delta_{j}^{T}=\mathbf{\Sigma}_{j}+\varepsilon\left({}^{c}{\mathbf{\Sigma}}_{j}-\mathbf{\Sigma}_{j}\right)+\varepsilon\,(1-\varepsilon)\,\Delta_{j}\,\Delta_{j}^{T}.$
(11)

Let $\mathbf{Z}$ be the random vector describing the mixture of the
contaminated random vectors $\mathbf{Z}_{j}$. From (1), the expectation of
$\mathbf{Z}$ is

$\operatorname*{\mathbb{E}}(\mathbf{Z})=\sum_{j=1}^{g}p_{j}\operatorname*{\mathbb{E}}(\mathbf{Z}_{j})={\boldsymbol{\mu}}+\varepsilon\,\overline{\Delta},\quad{\rm
where}\quad\overline{\Delta}=\sum_{j=1}^{g}p_{j}\,\Delta_{j}.$ (12)

From (2) and (11), we can compute the contaminated theoretical within
covariance matrix:

$\mathbf{W}_{Z}=\sum_{j=1}^{g}p_{j}\operatorname*{{\rm
Var}}(\mathbf{Z}_{j})=\mathbf{W}+\varepsilon\,\sum_{j=1}^{g}p_{j}\,({}^{c}{\mathbf{\Sigma}}_{j}-\mathbf{\Sigma}_{j})+\varepsilon\,(1-\varepsilon)\,\sum_{j=1}^{g}p_{j}\,\Delta_{j}\,\Delta_{j}^{T},$
(13)

where $\mathbf{W}=\sum_{j=1}^{g}p_{j}\,\mathbf{\Sigma}_{j}$ is the theoretical
within covariance matrix. We notice that when
${}^{c}{\mathbf{\Sigma}}_{j}=\mathbf{\Sigma}_{j}$ for all $j$, the
contaminated theoretical within covariance matrix is the original
$\mathbf{W}$, plus a low-rank matrix.

By using (2), we compute the theoretical between covariance matrix of
$\mathbf{Z}$ by substituting (10) and (12):

$\displaystyle\mathbf{B}_{Z}$
$\displaystyle=\sum_{j=1}^{g}p_{j}\,(\operatorname*{\mathbb{E}}(\mathbf{Z}_{j})-\operatorname*{\mathbb{E}}(\mathbf{Z}))\,(\operatorname*{\mathbb{E}}(\mathbf{Z}_{j})-\operatorname*{\mathbb{E}}(\mathbf{Z}))^{T}$
$\displaystyle=\sum_{j=1}^{g}p_{j}\,({\boldsymbol{\mu}}_{j}+\varepsilon\Delta_{j}-{\boldsymbol{\mu}}-\varepsilon\,\overline{\Delta})\,({\boldsymbol{\mu}}_{j}+\varepsilon\Delta_{j}-{\boldsymbol{\mu}}-\varepsilon\,\overline{\Delta})^{T}.$

Let us define
$\boldsymbol{\delta}_{j}={\boldsymbol{\mu}}_{j}-{\boldsymbol{\mu}}$. It
follows that

$\mathbf{B}_{Z}=\mathbf{B}+\varepsilon\sum_{j=1}^{g}p_{j}\,[\boldsymbol{\delta}_{j}(\Delta_{j}-\overline{\Delta})^{T}+(\Delta_{j}-\overline{\Delta})\,\boldsymbol{\delta}_{j}^{T}]+\varepsilon^{2}\,\sum_{j=1}^{g}p_{j}\,(\Delta_{j}-\overline{\Delta})\,(\Delta_{j}-\overline{\Delta})^{T}.$
(14)

Now suppose that only the first group has been contaminated. This means that
$\Delta_{j}=0$ for $j\neq 1$, and $\overline{\Delta}=p_{1}\Delta_{1}$.
Exploiting the fact that
$\sum_{j=2}^{g}p_{j}\boldsymbol{\delta}_{j}=-p_{1}\boldsymbol{\delta}_{1}$, we
obtain

$\displaystyle\sum_{j=1}^{g}p_{j}\,[\boldsymbol{\delta}_{j}\big{(}\Delta_{j}-\overline{\Delta}\big{)}^{T}+\big{(}\Delta_{j}-\overline{\Delta}\big{)}\boldsymbol{\delta}_{j}^{T}]$
$\displaystyle=p_{1}\,(1-p_{1})\,(\boldsymbol{\delta}_{1}\Delta_{1}^{T}+\Delta_{1}\boldsymbol{\delta}_{1}^{T})-p_{1}^{2}(\boldsymbol{\delta}_{1}\Delta_{1}^{T}+\Delta_{1}\boldsymbol{\delta}_{1}^{T})$
$\displaystyle=p_{1}\,(\boldsymbol{\delta}_{1}\Delta_{1}^{T}+\Delta_{1}\boldsymbol{\delta}_{1}^{T}),$

and
$\sum_{j=1}^{g}p_{j}\left(\Delta_{j}-\overline{\Delta}\right)\left(\Delta_{j}-\overline{\Delta}\right)^{T}=p_{1}\,(1-p_{1})\,\Delta_{1}\,\Delta_{1}^{T}$.
Then, up to second-order terms in $\varepsilon$, we can write the contaminated
covariance matrices as

$\displaystyle\mathbf{W}_{Z}$
$\displaystyle=\mathbf{W}+\Delta\mathbf{W}=\mathbf{W}+\varepsilon\,p_{1}\,[({}^{c}{\mathbf{\Sigma}}_{1}-\mathbf{\Sigma})+\Delta_{1}\,\Delta_{1}^{T}]+o(\varepsilon),$
$\displaystyle\mathbf{B}_{Z}$
$\displaystyle=\mathbf{B}\;+\Delta\mathbf{B}\;=\mathbf{B}+\varepsilon\,p_{1}\,(\boldsymbol{\delta}_{1}\Delta_{1}^{T}+\Delta_{1}\boldsymbol{\delta}_{1}^{T})+o(\varepsilon).$

In general, all the contaminated parameters
$\operatorname*{\mathbb{E}}(\mathbf{Z}_{j})$, $\operatorname*{{\rm
Var}}(\mathbf{Z}_{j})$, $\operatorname*{\mathbb{E}}(\mathbf{Z})$,
$\mathbf{W}_{Z}$, and $\mathbf{B}_{Z}$ can be written as their original
counterparts ($\operatorname*{\mathbb{E}}(\mathbf{X}_{j})$,
$\operatorname*{{\rm Var}}(\mathbf{X}_{j})$,
$\operatorname*{\mathbb{E}}(\mathbf{X})$, $\mathbf{W}$, and $\mathbf{B}$,
respectively) plus an additional term each, that quantifies the effect of the
contamination. The additional terms are functions of the percentage of
contamination, $\varepsilon$, the probability of the contaminated classes,
$p_{j}$, the parameters of the uncontaminated (${\boldsymbol{\mu}}_{j}$ and
$\mathbf{\Sigma}_{j}$) and contaminating distribution
(${}^{c}{{\boldsymbol{\mu}}}_{j}$ and ${}^{c}{\mathbf{\Sigma}}_{j}$).

###  3.2 Perturbation analysis of TR

In Section 3.1, we have given explicit expressions for the contaminated
theoretical between and within covariance matrices, when the contamination
happens in one class (LABEL:eq:cont-oneclass). We have stressed that they are
of the form $\mathbf{B}_{Z}=\mathbf{B}+\Delta\mathbf{B}$ and
$\mathbf{W}_{Z}=\mathbf{W}+\Delta\mathbf{W}$, where $\|\Delta\mathbf{B}\|$,
$\|\Delta\mathbf{W}\|\to 0$ when $\varepsilon\to 0$. We use this formulation
to understand how the TR solution, $\mathbf{V}_{\rm TR}$, changes when we
perturb the pair $(\mathbf{B},\mathbf{W})$ to
$(\mathbf{B}_{Z},\,\mathbf{W}_{Z})$. A generalization of this problem is
addressed by Zhang and Yang [27], for arbitrary perturbations such that
$\|\Delta\mathbf{B}\|$, $\|\Delta\mathbf{W}\|\to 0$. We report their upper
bound for the resulting change in $\mathbf{V}_{\rm TR}$. By
$[\Delta\mathbf{B}\ \,\Delta\mathbf{W}]$ we denote the $p\times 2p$ matrix
constructed by appending $\Delta\mathbf{W}$ to $\Delta\mathbf{B}$.

######  Theorem 3.13.

[27, Thm. 4.1] Let $\mathbf{B}$ be symmetric and $\mathbf{W}$ SPD. Let the
$p\times k$ matrix $\mathbf{V}_{\rm TR}$ be the solution to (9) for the pair
$(\mathbf{B},\mathbf{W})$ and suppose that there is a gap $\gamma>0$ between
the $k$th and the $(k+1)$st largest eigenvalues of
$\mathbf{B}-\rho(\mathbf{V}_{\rm TR})\,\mathbf{W}$. Then, when $\eta\to 0$,
for all $\|[\Delta\mathbf{B}\ \,\Delta\mathbf{W}]\|\leq\eta$, there exists
${}^{c}{\mathbf{V}_{\rm TR}}$ solution to (9) for the pencil
$(\mathbf{B}+\Delta\mathbf{B},\,\mathbf{W}+\Delta\mathbf{W})$ such that

$\|\mathbf{V}_{\rm TR}\mathbf{V}_{\rm TR}^{T}-{}^{c}{\mathbf{V}_{\rm
TR}}{}^{c}{\mathbf{V}_{\rm
TR}^{T}}\|\leq\frac{4\sigma}{\gamma}\,\Big{(}1+\frac{k}{\tau}\,\|\mathbf{W}+\Delta\mathbf{W}\|\Big{)}+o(\eta),$

where $\sigma=\|\Delta\mathbf{B}\|+\rho(\mathbf{V}_{\rm
TR})\,\|\Delta\mathbf{W}\|\leq(1+\rho(\mathbf{V}_{\rm TR}))\,\eta$, and $\tau$
is the sum of the $k$ smallest eigenvalues of $\mathbf{W}$.

Note that, as stated in the introduction, we always have $\|\mathbf{V}_{\rm
TR}\mathbf{V}_{\rm TR}^{T}-{}^{c}{\mathbf{V}_{\rm TR}}{}^{c}{\mathbf{V}_{\rm
TR}^{T}}\|\leq 1$. This theorem gives an asymptotic upper bound that is linear
in $\eta$, for $\eta\to 0$.

When the contamination level tends to zero, i.e., $\varepsilon\to 0$, we can
express $\Delta\mathbf{B}$ and $\Delta\mathbf{W}$ as in (LABEL:eq:cont-
oneclass). Let us assume that
${}^{c}{\mathbf{\Sigma}}_{1}=\mathbf{\Sigma}_{1}$. Then

$\|\Delta\mathbf{B}\|\leq
2\,\varepsilon\,p_{1}\,\|\boldsymbol{\delta}_{1}\|\,\|\Delta_{1}\|+o(\varepsilon),\qquad\|\Delta\mathbf{W}\|\leq\,\varepsilon\,p_{1}^{2}\,\|\Delta_{1}\|^{2}+o(\varepsilon).$

Thus,

$\|[\Delta\mathbf{B}\
\,\Delta\mathbf{W}]\|\leq\,\varepsilon\,p_{1}\,\|\Delta_{1}\|\,(2\,\|\boldsymbol{\delta}_{1}\|\,+\,p_{1}\,\|\Delta_{1}\|)+o(\varepsilon).$

If $p_{1}$, $\|\boldsymbol{\delta}_{1}\|$, and $\|\Delta_{1}\|$ are fixed,
then we have just shown that $\|[\Delta\mathbf{B}\
\,\Delta\mathbf{W}]\|\leq\varepsilon\,C$, where
$C=p_{1}\,\|\Delta_{1}\|\,(2\,\|\boldsymbol{\delta}_{1}\|\,+\,p_{1}\,\|\Delta_{1}\|)>0$.
If we set $\eta=\varepsilon\,C$ in Theorem 3.13, then
$o(\eta)=o(\varepsilon)$. Thus we can write the upper bound for the
dissimilarity measure between $\mathbf{V}_{\rm TR}$ and
${}^{c}{\mathbf{V}_{\rm TR}}$:

$\displaystyle\|\mathbf{V}_{\rm TR}\mathbf{V}_{\rm
TR}^{T}-{}^{c}{\mathbf{V}_{\rm TR}}{}^{c}{\mathbf{V}_{\rm TR}^{T}}\|$
$\displaystyle\leq\frac{4}{\gamma}\,\varepsilon\,p_{1}\,\|\Delta_{1}\|\,\big{(}2\,\|\boldsymbol{\delta}_{1}\|+p_{1}\,\rho(\mathbf{V}_{\rm
TR})\,\|\Delta_{1}\|\big{)}$
$\displaystyle\qquad\qquad\cdot\Big{(}1+\frac{k}{\tau}\,\left(\|\mathbf{W}\|+2\,\varepsilon\,p_{1}^{2}\,\|\Delta_{1}\|^{2}\right)\Big{)}+o(\varepsilon)$
$\displaystyle=\frac{4}{\gamma}\,\varepsilon\,p_{1}\,\|\Delta_{1}\|\,\big{(}2\,\|\delta_{1}\|+p_{1}\,\rho(\mathbf{V}_{\rm
TR})\,\|\Delta_{1}\|\big{)}\,\Big{(}1+\frac{k}{\tau}\,\|\mathbf{W}\|\Big{)}+o(\varepsilon),$

where in the last step we discard a new $o(\varepsilon)$-term. Let
$\lambda_{p-k+1}(\mathbf{W})\geq\cdots\geq\lambda_{p}(\mathbf{W})$ be the $k$
smallest eigenvalues of $\mathbf{W}$; then
$\lambda_{1}(\mathbf{W})=\|\mathbf{W}\|$ and
$\tau=\sum_{i=1}^{k}\lambda_{p-k+i}(\mathbf{W})$. Since $\tau\geq
k\,\lambda_{p}(\mathbf{W})$, we can write:

$\|\mathbf{V}_{\rm TR}\mathbf{V}_{\rm TR}^{T}-{}^{c}{\mathbf{V}_{\rm
TR}}{}^{c}{\mathbf{V}_{\rm
TR}^{T}}\|\leq\frac{4}{\gamma}\,\varepsilon\,p_{1}\,\|\Delta_{1}\|\left(2\,\|\delta_{1}\|+p_{1}\,\rho(\mathbf{V}_{\rm
TR})\,\|\Delta_{1}\|\right)\,\Big{(}1+\frac{\lambda_{1}(\mathbf{W})}{\lambda_{p}(\mathbf{W})}\Big{)}+o(\varepsilon).$
(16)

The quotient $\lambda_{1}(\mathbf{W})/\lambda_{p}(\mathbf{W})$ is the
condition number of $\mathbf{W}$, which is an indication of the ill-
conditioning of $\mathbf{W}$. In the situation of a high condition number,
even small perturbations in the $\mathbf{W}$ matrix can lead to big changes in
the subspaces generated by TR.

The upper limit (16) is linear in the distance between
${\boldsymbol{\mu}}_{1}$ and ${\boldsymbol{\mu}}$, i.e.,
$\|\boldsymbol{\delta}_{1}\|$, but quadratic in $\|\Delta_{1}\|$, the distance
between the non-contaminated mean vector and the contaminated mean vector.
Thus $\|\Delta_{1}\|$, which is related to the contamination process, is the
term that contributes the most to increasing the upper limit of this
dissimilarity index. Moreover, the upper limit depends explicitly on the
theoretical within covariance matrix, $\mathbf{W}$, but only implicitly on
$\mathbf{B}$ (via $\rho(\mathbf{V}_{\rm TR})$ and $\gamma$), assigning
different levels of importance to these two matrices.

A way to control the impact of contamination in the estimated TR and FDA
coefficients is to consider robust estimators, which are discussed in the next
section.

##  4 Robust trace ratio and robust Fisher discriminant analysis

In this section, we propose a robust version of the TR algorithm, and also
study a robust version of FDA. Instead of considering the classical estimators
of the between and within covariance matrices defined in (3), we use their
robust versions based on the minimum covariance determinant (MCD) estimators
of location and covariance matrix, described in [28]. Let us denote the robust
estimators by $\widetilde{{\boldsymbol{\mu}}}_{j}$ and
$\widetilde{\mathbf{S}}_{j}$ for the location and the covariance matrix,
respectively, of the $j$th group, $j=1,\dots,g$. The global location estimator
is defined as
$\widetilde{{\boldsymbol{\mu}}}=\sum_{j=1}^{g}\tfrac{n_{j}}{n}\,\widetilde{{\boldsymbol{\mu}}}_{j}$.
The robust versions of the between and within scatter matrices are then
defined as

$\widetilde{\mathbf{B}}=\sum_{j=1}^{g}\tfrac{n_{j}}{n}\,(\widetilde{{\boldsymbol{\mu}}}-\widetilde{{\boldsymbol{\mu}}}_{j})\,(\widetilde{{\boldsymbol{\mu}}}-\widetilde{{\boldsymbol{\mu}}}_{j})^{T},\qquad\widetilde{\mathbf{W}}=\sum_{j=1}^{g}\tfrac{n_{j}-1}{n-g}\,\widetilde{\mathbf{S}}_{j}.$
(17)

We now propose a robust version of the trace ratio method, which we denote by
rTR. For this method, we determine $\widetilde{\mathbf{V}}_{\rm TR}$ from
Algorithm 1 in Subsection 2.3, but now using as input the robust estimates of
the between and within scatter matrices, $\widetilde{\mathbf{B}}$ and
$\widetilde{\mathbf{W}}$, respectively.

We compare rTR with a robust version of FDA (rFDA), where a robust coefficient
matrix, $\widetilde{\mathbf{V}}_{\rm FDA}$, is obtained by finding the first
$k$ generalized eigenvectors of the pencil
$(\widetilde{\mathbf{B}},\widetilde{\mathbf{W}})$. For the classical
estimates, the statistical literature usually recommends the derivation of the
eigenvectors of the matrix ${\mathbf{W}}^{-1}{\mathbf{B}}$, as can be
confirmed in, e.g., [1] or [29]. Nevertheless, we note that standard methods
working on this pair (such as the ones implemented in the R package geigen,
see [30] for details) are numerically preferable over approaches based on the
transformed matrix $\mathbf{W}^{-1}\mathbf{B}$ (see, e.g., [31, Ch. 8]).
Therefore, we use the function geigen from [30] to derive the classical and
robust FDA matrix of coefficients.

The procedure to classify a test set observation, $\mathbf{x}_{0}$, into one
of the groups, for both rTR and rFDA, is as follows.

  1. (i)

Project the training set observations on the subspace spanned by
$\widetilde{\mathbf{V}}$ (estimated by rTR or rFDA).

  2. (ii)

Compute the MCD estimates of the group projected locations
$\widetilde{{\boldsymbol{\mu}}}_{Uj}$ and the covariance matrices, and from
these the associated $k\times k$ pooled covariance matrix
$\widetilde{\mathbf{W}}_{U}$.

  3. (iii)

Compute $\mathbf{u}_{0}=\widetilde{\mathbf{V}}^{T}\mathbf{x}_{0}$ and the
Mahalanobis distance between $\mathbf{u}_{0}$ and
$\widetilde{{\boldsymbol{\mu}}}_{Uj}$,
$\|\mathbf{u}_{0}-\widetilde{{\boldsymbol{\mu}}}_{Uj}\|_{\widetilde{\mathbf{W}}_{U}^{-1}}$.

  4. (iv)

Assign $\mathbf{x}_{0}$ to the group $j$ that minimizes the quantity
$\|\mathbf{u}_{0}-\widetilde{{\boldsymbol{\mu}}}_{Uj}\|_{\widetilde{\mathbf{W}}_{U}^{-1}}^{2}-2\,\log(n_{j}/n)$.
The second term balances the class prior information in the decision rule.

Step (iv) is equivalent to performing robust LDA over the projected data on
the reduced subspace, spanned by a robust estimate of $\mathbf{V}$. The robust
LDA is obtained by plugging the MCD estimators of location and covariance in
the Mahalanobis distances, computed in the reduced subspace. This is a robust
way to mitigate the potentially harmful effect of outliers in the
classification rules based on the Euclidean distance in the projected space.
We remark that, in case of high-dimensional data, we may also use regularized
MCD [13] or any other robust estimators of the groups location and covariance
matrix.

In the next sections, we compare FDA and TR when either classical or robust
estimates for the covariance matrices are considered.

##  5 Simulation study

We first explore TR and FDA as classifiers. Here the dimensionality reduction
is embedded in the classifier: first, we build the projection matrix, either
$\mathbf{V}_{\rm FDA}$ or $\mathbf{V}_{\rm TR}$, and then we classify the
projected data points based on the distances from the projected group means.
The simulation scheme presented in this section has two main objectives: (i)
compare TR with FDA classical and robust classifiers, under clean and
contaminated data, and (ii) analyze the performance of TR and FDA classifiers
under the presence of irrelevant variables for the associated classification
problem.

###  5.1 Scenarios

In all scenarios, we consider four groups ($g=4$), each one characterized by a
random vector of dimension $p=3$, except for the last scenario, where we
continue to have three relevant variables for the classification problem, but
we add $q\in\\{0,10,20,50,100\\}$ other irrelevant variables. The
contamination scheme is the one described in Section 3.1. More precisely, in
the $\ell$th scenario, the clean $j$th group is
$\mathbf{X}_{j}=(\mathbf{X}\mid
Y=j)\sim\mathcal{N}_{p}\big{(}{\boldsymbol{\mu}}_{j}^{(\ell)}\,,\,\mathbf{\Sigma}_{j}^{(\ell)}\big{)}$
and the contaminating variable is distributed as
${}^{c}{\mathbf{X}}_{j}\sim\mathcal{N}_{p}\big{(}{}^{c}{\boldsymbol{\mu}}_{j}^{(\ell)}\,,\,\mathbf{\Sigma}_{j}^{(\ell)}\big{)}$,
for $j=1,\dots,g$,
$\ell=\mbox{\scriptsize{\rm{I}}},\dots,\mbox{\scriptsize{\rm{IV}}}$.
Contamination is absent when $\varepsilon=0$.

For all scenarios, we choose the group mean vectors to suggest that two
discriminant variables are enough to separate the groups, thus we choose $k=2$
as the reduced dimension. To make the scenarios more realistic, the covariance
matrices challenge this separation.

As a benchmark for each scenario, we use the associated parameters of the
regular population (${\boldsymbol{\mu}}_{j}^{(\ell)}$ and
$\mathbf{\Sigma}_{j}^{(\ell)}$) to write the theoretical between and within
scatter matrices. Using these as input, we find the solutions to (6) and (9),
and refer to them as theoretical matrices of coefficients, which are indicated
by ${\bf V}_{\rm FDA}$ and ${\bf V}_{\rm TR}$, respectively.

Table 1: Parameters characterizing regular and contaminated distributions for
each simulation scenario
($\ell=\mbox{\scriptsize{\rm{I}}},\dots,\mbox{\scriptsize{\rm{IV}}}$), where
$\mathbf{M}^{(\ell)}=[{\boldsymbol{\mu}}_{1}^{(\ell)},\dots,{\boldsymbol{\mu}}_{4}^{(\ell)}]$,
${}^{c}\mathbf{M}^{(\ell)}=[{}^{c}{\boldsymbol{\mu}}_{1}^{(\ell)},\dots,{}^{c}{\boldsymbol{\mu}}_{4}^{(\ell)}]$,
and $\mathbf{\Sigma}_{j}^{(\ell)}$ is the covariance matrix of the $j$th group
in the $\ell$th scenario.

$\ell$ | $\mathbf{M}^{(\ell)}$ | ${}^{c}\mathbf{M}^{(\ell)}$ | $\mathbf{\Sigma}_{1}^{(\ell)}$ | $\mathbf{\Sigma}_{2}^{(\ell)}$ | $\mathbf{\Sigma}_{3}^{(\ell)}$ | $\mathbf{\Sigma}_{4}^{(\ell)}$
---|---|---|---|---|---|---
I | $\left[\\!\\!\begin{array}[]{rrrr}15&15&&\\\ 3&-3&&\\\ &&2&-2\end{array}\\!\\!\right]$ | $\left[\\!\\!\begin{array}[]{rrrr}15&15&&\\\ -27&-3&&\\\ &&2&-2\end{array}\\!\\!\right]$ | $\left[\\!\\!\begin{array}[]{lll}1&&\\\ &\phantom{-}1&-0.25\\\ &-0.25&\phantom{-}1\end{array}\\!\\!\right]$ | $\mathbf{\Sigma}_{1}^{(\mbox{\scriptsize{\rm{I}}})}$ | $\mathbf{\Sigma}_{1}^{(\mbox{\scriptsize{\rm{I}}})}$ | $\mathbf{\Sigma}_{1}^{(\mbox{\scriptsize{\rm{I}}})}$
II | $\left[\\!\\!\begin{array}[]{rrrr}&&3&-3\\\ -3&3&&\\\ &&1&1\end{array}\\!\\!\right]$ | $\left[\\!\\!\begin{array}[]{rrrr}&&3&-3\\\ -3&-27&&\\\ &&1&1\end{array}\\!\\!\right]$ | $\left[\\!\\!\begin{array}[]{rrrr}1&&\\\ &3&\\\ &&1\end{array}\\!\\!\right]$ | $\mathbf{\Sigma}_{1}^{(\mbox{\scriptsize{\rm{II}}})}$ | $\left[\\!\\!\begin{array}[]{rrrr}3&&\\\ &1&\\\ &&3\end{array}\\!\\!\right]$ | $\mathbf{\Sigma}_{3}^{(\mbox{\scriptsize{\rm{II}}})}$
III | $\left[\\!\\!\begin{array}[]{rrrr}&10&&-10\\\ -3&&3&\\\ 1&&1&\end{array}\\!\\!\right]$ | $\left[\\!\\!\begin{array}[]{rrrr}&10&&-10\\\ -3&&-27&\\\ 1&&1&\end{array}\\!\\!\right]$ | $\mathbf{\Sigma}_{1}^{(\mbox{\scriptsize{\rm{II}}})}$ | $\mathbf{I}_{3}$ | $\left[\\!\\!\begin{array}[]{lll}1&&\\\ &\phantom{-}3&-0.5\\\ &-0.5&\phantom{-}1\end{array}\\!\\!\right]$ | $\mathbf{I}_{3}$
IV | $\left[\\!\\!\begin{array}[]{c}\mathbf{M}^{(\mbox{\scriptsize{\rm{II}}})}\\\ \mathbf{0}_{q\times 4}\end{array}\\!\\!\right]$ | $\left[\\!\\!\begin{array}[]{c}{}^{c}\mathbf{M}^{(\mbox{\scriptsize{\rm{II}}})}\\\ \mathbf{0}_{q\times 4}\end{array}\\!\\!\right]$ | $\left[\\!\\!\begin{array}[]{rrrr}\mathbf{\Sigma}_{1}^{(\mbox{\scriptsize{\rm{II}}})}&\mathbf{0}_{3\times q}\\\ \mathbf{0}_{q\times 3}&\mathbf{I}_{q}\end{array}\\!\\!\right]$ | $\mathbf{\Sigma}_{1}^{(\mbox{\scriptsize{\rm{IV}}})}$ | $\left[\\!\\!\begin{array}[]{rrrr}\mathbf{\Sigma}_{3}^{(\mbox{\scriptsize{\rm{II}}})}&\mathbf{0}_{3\times q}\\\ \mathbf{0}_{q\times 3}&\mathbf{I}_{q}\end{array}\\!\\!\right]$ | $\mathbf{\Sigma}_{3}^{(\mbox{\scriptsize{\rm{IV}}})}$

Note: $\mathbf{0}_{q\times q^{\prime}}$ represents a $q\times q^{\prime}$
matrix of zeros. In scenario IV, $q\in\\{0,\,10,\,20,\,50,\,100\\}$ represents
the number of irrelevant variables. Omitted matrix entries represent zeros.

(a) Scenario I.

(b) Scenario III.

(c) Scenario II.

Figure 2: Mean vectors and 95% confidence ellipses for the projected data
assuming that the projected data follows a bivariate normal distribution. The
projection matrices are the theoretical $\mathbf{V}_{\rm TR}$ and
$\mathbf{V}_{\rm FDA}$, corresponding to the contamination levels
$\varepsilon=0$ (solid line) and $\varepsilon=0.1$ (dashed line).

When all group covariance matrices are equal, i.e., in the first scenario
($\ell=\mbox{\scriptsize{\rm{I}}}$), and the underlying group population
follows a multivariate normal distribution, the classification rule that
minimizes (maximizes) the total probability of misclassification (probability
of correctly assign an observation to a group, called theoretical accuracy) is
linear (LDA, cf. Section 2.2) and would coincide with FDA, if we considered
$k=3$. In the present work, we only consider $k=2$, which does not enjoy of
the same optimality property. In addition, we see that TR leads to better
results than FDA in this scenario, showing that the missing dimension has a
relevant impact on the total probability of misclassification. The solutions
to (6) and (9) associated with the parameters characterizing the regular
distribution are:

${\mathbf{V}}^{(\mbox{\scriptsize{\rm{I}}})}_{\rm
FDA}={\mbox{\scriptsize$\left[\\!\\!\begin{array}[]{cl}1&0\\\ 0&1\\\
0&0.415\end{array}\\!\\!\right]$}},\qquad{\mathbf{V}}^{(\mbox{\scriptsize{\rm{I}}})}_{\rm
TR}={\mbox{\scriptsize$\left[\\!\\!\begin{array}[]{cl}1&0\\\ 0&0.757\\\
0&0.654\end{array}\\!\\!\right]$}}.$ (18)

In this case, the first two FDA discriminant variables are $X_{1}$ and a
linear combination of $X_{2}$ and $X_{3}$, where $X_{2}$ has a higher weight
than $X_{3}$. The TR matrix of coefficients has a similar structure, but gives
a more balanced weight to $X_{2}$ and $X_{3}$. It should be noted that these
weights cannot be compared one-to-one because the classification is based on a
different metric: while the Mahalanobis distance is used for FDA, the
Euclidean norm is exploited for classification based on the trace ratio
method.

To get a better understanding, we visualize the effect of the contamination on
the projected mean vector and covariance matrix in each scenario. In Figure 2
we draw 95% confidence ellipses for the projected data, assuming that the
projected data follows a bivariate normal distribution with mean vector and
covariance matrix equal to the theoretical values under each scenario. We
consider $\varepsilon=0$ and $\varepsilon=0.10$; in the latter case, the
projected contaminated group does not follow a bivariate normal distribution,
so the 95% confidence ellipse would not be an appropriate confidence region.
Nevertheless, this approximation serves its purpose.

For scenario $\ell=\mbox{\scriptsize{\rm{I}}}$, the contamination strongly
affects the first and second group mean vectors, leading to an overlap between
these two groups in the projected space spanned by $\mathbf{V}_{\rm
FDA}^{\mbox{\scriptsize{\rm{(I)}}}}$ and $\mathbf{V}_{\rm
TR}^{\mbox{\scriptsize{\rm{(I)}}}}$. As a side effect, the separation between
the third and fourth group increases for both FDA and TR, as shown in Figure
2(a).

In the following scenarios, we consider different class covariance matrices.
The second scenario ($\ell=\mbox{\scriptsize{\rm{II}}}$) describes a case
where TR and FDA show a similar performance. The covariance matrices are all
diagonal (see Table 1 for more details). The theoretical matrices of
coefficients are

${\mathbf{V}}^{(\mbox{\scriptsize{\rm{II}}})}_{\rm
FDA}={\mathbf{V}}^{(\mbox{\scriptsize{\rm{II}}})}_{\rm
TR}={\mbox{\scriptsize$\left[\\!\\!\begin{array}[]{cc}1&0\\\ 0&1\\\ 0&0\\\
\end{array}\\!\\!\right]$}}.$ (19)

As illustrated in Figure 2(c), in the projected space we cannot distinguish
between groups 1 and 2, but the contamination also has a considerable impact
on the other two groups, mainly on the shapes of the covariance matrices, as
illustrated by the changes in the 95% confidence ellipses.

In the third scenario ($\ell=\mbox{\scriptsize{\rm{III}}}$), we illustrate an
example where FDA leads to better results than TR. Like in the previous
scenario, by looking at the group means (cf. Table 1), $X_{1}$ and $X_{2}$
seem to be enough to separate the groups. However, the correlation between
$X_{2}$ and $X_{3}$ in the third group affects $\mathbf{v}_{2}$, where both
variables have a nonzero weight. We can anticipate that FDA performs better
than TR since it gives higher weight to $X_{2}$ than to $X_{3}$, as follows:

${\mathbf{V}}^{(\mbox{\scriptsize{\rm{III}}})}_{\rm
FDA}={\mbox{\scriptsize$\left[\\!\\!\begin{array}[]{cl}1&0\\\ 0&1\\\
0&0.125\end{array}\\!\\!\right]$}},\qquad{\mathbf{V}}^{(\mbox{\scriptsize{\rm{III}}})}_{\rm
TR}={\mbox{\scriptsize$\left[\\!\\!\begin{array}[]{cl}1&0\\\ 0&0.147\\\
0&0.989\end{array}\\!\\!\right]$}}.$ (20)

This is reinforced by Figure 2(b), where the 95% confidence region of the
first and third group show a high overlap in the subspace spanned by TR, when
no contamination is present. Contamination increases this overlap but does not
seem to affect groups 2 and 4. For FDA, we observe a considerable separation
between groups 1 and 3 in the projected space, when no contamination occurs.
The contamination scheme substantially increases the overlap between groups 1
and 3 in the projected space, as can be confirmed by Figure 2(b).
Additionally, contamination affects the variability of groups 2 and 4.

In the fourth scenario ($\ell=\mbox{\scriptsize{\rm{IV}}}$) we test the
performances of FDA and TR in presence of irrelevant variables in the data. In
[12] a simulation scheme is designed to study the effect of increasing the
number of irrelevant variables for the classification problem, when keeping
the number of relevant variables and the group sample sizes fixed. The authors
have shown that, for clean and contaminated data, the accuracy of both
classical and robust linear discriminant classifiers decreases considerably,
as the number of irrelevant variables increases. To verify whether this
behavior also occurs in TR and FDA, we design the fourth scenario
($\ell=\mbox{\scriptsize{\rm{IV}}}$) by adding $q\in\\{0,10,20,50,100\\}$
irrelevant variables to the second scenario. We choose scenario II as starting
point, since it is the one with no significant differences in the median
accuracy between theoretical TR and FDA (cf. Section 5.2). The irrelevant
variables are independent and identically distributed with zero mean and
unitary variance; they are also uncorrelated from the potentially relevant
variables, i.e., $X_{1}$, $X_{2}$, and $X_{3}$. The parameters associated with
the fourth scenario are listed in Table 1. From the theoretical viewpoint, the
relevant variables are well captured by the matrices of coefficients:

${\mathbf{V}}^{(\mbox{\scriptsize{\rm{IV}}})}_{\rm
FDA}={\mathbf{V}}^{(\mbox{\scriptsize{\rm{IV}}})}_{\rm
TR}={\mbox{\scriptsize$\left[\\!\\!\begin{array}[]{cc}1&0\\\ 0&1\\\
0&0\\\\[-2.84526pt] \vdots&\vdots\\\ 0&0\\\ \end{array}\\!\\!\right]$}}.$ (21)

In the following section, we are going to generate data from the scenarios in
Table 1 to study the sample behavior of $\mathbf{V}_{\rm TR}$ and
$\mathbf{V}_{\rm FDA}$ under different contamination levels and different
estimators for the between and within covariance matrices.

###  5.2 Experimental setting

Given a scenario, for each experiment we generate a training set of
$n_{j}=400$ data points for each class, and a test set composed by $n_{j}=40$
observations per group, according to the parameters defined in Table 1. For
scenario IV we also consider $n_{j}=40$ for the training set. Training sets
and test sets are independently generated. We study the different
contamination levels
$\varepsilon\in\\{0,\,0.05,\,0.10,\,0.15,\,0.20,\,0.25,\,0.30\\}$. We test the
performances of the following classifiers.

  * •

cTR and cFDA: TR and FDA with classical estimators for the covariance
matrices. The associated classification rules are described in Section 2.

  * •

rTR and rFDA: TR and FDA with robust estimators for the covariance matrices,
as described in Section 4. These are obtained by the MRCD robust location and
covariance matrix estimators via the CovMrcd routine from the R package rrcov,
where the parameter ${\sf alpha}=0.75$ controls the size of the subsets over
which the covariance determinant is minimized (see [6, 32] for further
details). The classification rules are described in Section 4. Given that the
MRCD estimates select $75\%$ of the data, we expect that they are able to deal
with datasets with at most $25\%$ of outliers. Finally, we remark that the
regularization parameter of the MRCD estimates is chosen adaptively, so that
MRCD estimates coincide with MCD estimates when regularization is not needed
[13].

  * •

tTR and tFDA: theoretical TR and FDA, where the covariance matrices are
computed analytically (cf. (2)). In this case, no training set is needed. The
associated classification rules are described in Section 2 and make use of the
theoretical group means.

  * •

tQDA: quadratic discriminant classification rule (see, e.g., [29]) for
multivariate normal data. The parameters of the distributions are known, so we
do not need a training set.

  * •

AdaBoost: nonlinear classifier (with decision trees as the weak learners)
[33]. This is abbreviated by ab in the figures. We have selected AdaBoost due
to its popularity (and proven performance) in the Machine Learning community.

We remark that tQDA and AdaBoost do not involve any linear dimensionality
reduction strategy. When the groups have a multivariate normal distribution
with different covariance matrices, tQDA is the optimal classifier, because it
maximizes the accuracy (or, equivalently, it minimizes the total probability
of misclassification). In scenario I, the covariance matrices are all equal,
so QDA coincides with LDA, which is the optimal linear classifier.

In cTR, cFDA, rTR, and rFDA, the training set is used to find a projection
matrix for each method, while we assess the performances of the classification
rules on the test set. Each experiment is repeated $200$ times, meaning that
for each contamination level $\varepsilon$, we generate $200$ training sets
and $200$ test sets. To get trustful performance measures, the test set is
never contaminated.

For each algorithm and scenario, we report the sample accuracy. We also
evaluate the goodness of a sample subspace $\widehat{\mathbf{V}}$ compared to
its theoretical counterpart $\mathbf{V}$. To do so, we compute the angles
between these two subspaces $\angle(\mathbf{V},\widehat{\mathbf{V}})$. The
matrices $\mathbf{V}$ and $\widehat{\mathbf{V}}$ need to have orthogonal
columns: these can be obtained from the solution to FDA via a QR
decomposition. The solution to TR already has orthogonal columns. Angles
compare only subspaces: they indicate the quality of the estimates of the
theoretical coefficients obtained with FDA and TR.

###  5.3 Results and discussion

We now report the results of the simulation study. Table 2 summarizes the
ranking of the algorithms based on the median accuracy over the generated test
datasets. Relative positions in Table 2 are tested via the Wilcoxon signed-
rank test for paired samples. Except for AdaBoost, the comparison is made only
under the hypothesis that the mean vectors and covariance matrices are known
and not contaminated. As we will see, the simulation study agrees with the
considerations that we have made in Section 5.1, when no contamination occurs:
in the first scenario, tTR shows a better accuracy than tFDA; in the second
and fourth scenario, the two methods are equivalent; in the third scenario,
tFDA gives better results than tTR.

Table 2: Ranking of theoretical methods and AdaBoost, with $\varepsilon=0$, according to the conclusions of the Wilcoxon signed-rank tests for paired samples, with the False Discovery Rate correction [34]. The median accuracy is reported in brackets. I |  tQDA (0.994) | $>$ |  AdaBoost (0.988) | $>$ |  tTR (0.969) | $>$ |  tFDA (0.900)
---|---|---|---|---|---|---|---
II, IV |  tQDA (0.894) | $>$ |  tFDA (0.881) | $=$ |  tTR (0.881) | $>$ |  AdaBoost (0.875)
III |  tQDA (0.981) | $=$ |  tFDA (0.981) | $>$ |  AdaBoost (0.975) | $>$ |  tTR (0.838)

Figure 3 shows the angles between the benchmark matrices $\mathbf{V}_{\rm TR}$
or $\mathbf{V}_{\rm FDA}$ and the solution to TR or FDA when $\mathbf{B}$ and
$\mathbf{W}$ are estimated from a sample. We consider both classical and
robust estimates. Since we are able to compute the theoretical contaminated
matrices ${}^{c}{\mathbf{B}}$ and ${}^{c}{\mathbf{W}}$ from Section 3, we can
show the theoretical angles between the original subspaces and the
contaminated ones (represented in gray in Figure 3).

Figure 3: Angle between the solution to TR (FDA) when the $\mathbf{B}$ and
$\mathbf{W}$ are the theoretical ones, and the solution to TR (FDA) in
presence of contamination when $\mathbf{B}$ and $\mathbf{W}$ are estimated
from a sample. We consider both classical (blue) and robust (orange)
estimates. The gray lines correspond to the theoretical angle between the
original subspace with and without contamination. The horizontal axis shows
different levels of contamination in the sample. For scenario IV, we only
consider $q=10$.

Figures 4 to 6 report the boxplots of the accuracy for each method and level
of contamination. We indicate as benchmark results those related to tTR, tLDA,
and tQDA, since these are not affected by the contamination level. Given the
theoretical findings of the contaminated scatter matrices, derived in Section
3.1, it is possible to obtain the accuracy for the theoretical tFDA and tTR
for each level of contamination. The results based on the theoretical
derivations and solely on simulations are quite similar, so only cFDA and cTR
(simulation-based) are shown.

In the first scenario ($\ell=\mbox{\scriptsize{\rm{I}}}$), and when no
contamination is considered, the Friedman test [35] concludes that there are
statistically significant differences between tQDA, AdaBoost, tTR, and tFDA.
Moreover, Table 2 reports the following ranking of the methods, based on
Wilcoxon signed-rank test for the median accuracy: tQDA, AdaBoost, tTR, and
tFDA. This ranking is also supported by the plots in Figure 4. In this case,
it is known that tQDA reduces to LDA, which would coincide with FDA if we
considered $k=3$ instead of $k=2$. FDA with $k=2$ loses almost 10% in median
accuracy by not considering a third dimension on the projected space, while
tTR leads to higher median accuracy than tFDA, when $k=2$. In addition,
AdaBoost turns out to be the second best method. Under contamination, we
observe that the robust methods can accommodate until 25% of contamination
while the classical counterparts break down, presenting an increase in
variability. From Figure 3, we see that the subspaces obtained with cFDA and
cTR change as soon as $\varepsilon\geq 0.05$; cFDA changes slightly more than
cTR. The subspaces generated by rFDA and rTR stay close to their theoretical
counterparts, as expected, although the robust methods also start breaking
down at $\varepsilon=0.25$.

Figure 4:  Distribution of the accuracy over 200 simulation runs of Scenario
I. The accuracy of tTR, tFDA, and tQDA is reported as benchmark in the
leftmost plot. Classical/robust TR and FDA classifiers, and AdaBoost, are
considered under several levels of contamination.

In the second scenario ($\ell=\mbox{\scriptsize{\rm{II}}}$), when
$\varepsilon=0$, the Friedman tests followed by the Wilcoxon signed-rank test
(in Table 2) lead to the following order of the methods: tQDA, tFDA and tTR
with (statistically) the same median accuracy, AdaBoost; this is agrees with
Figure 5. We stress that this is an example where linear approximations (tTR
and tFDA) of the optimal quadratic rule lead to better results than a
nonlinear classification algorithm (AdaBoost). When the data is not
contaminated, classical and robust linear estimators perform similarly to
their theoretical counterparts, and better than AdaBoost. When we increase the
proportion of contamination, Figure 5 shows that cTR is immediately affected
by the contamination, while the accuracy of cFDA gradually decreases, but ends
in the lowest median accuracy for $\varepsilon=0.30$, compared to all other
methods. Both robust alternatives are unaffected by small to moderate
percentages of contamination, but break down for $\varepsilon=0.30$, as
expected. AdaBoost does not seem to be affected by the contamination. In
Figure 3, the subspaces obtained with cFDA and cTR change as soon as
$\varepsilon\geq 0.05$; cTR changes more than cFDA, which also shows more
variability than cTR for $\varepsilon=0.05,\,0.1$. As expected, the subspaces
generated by rFDA and rTR are close to their theoretical counterparts; we
observe their breakdown starting from $\varepsilon=0.25$.

Figure 5: Distribution of the accuracy over 200 simulation runs of Scenario
II. The accuracy of tTR, tFDA, and tQDA is reported as benchmark in the
leftmost plot. Classical/robust TR and FDA classifiers, and AdaBoost, are
considered under several levels of contamination.

The third scenario ($\ell=\mbox{\scriptsize{\rm{III}}}$) describes a situation
where there is no statistically significant difference between tQDA and tFDA,
both having a higher median accuracy than AdaBoost; tTR is the method with the
lowest median accuracy. This is confirmed by the Friedman test followed by the
Wilcoxon signed-rank tests, and illustrated in Figure 6. The cTR method is
highly affected by the presence of contamination, showing a lower accuracy,
while the performance degradation of cFDA is slower with the increase of
$\varepsilon$. The increasing contamination level also causes an immediate
increase in the angle between the subspace generated by cFDA and its
theoretical counterpart. This does not hold for cTR, where the contamination
causes only a little change in the angle. We find the same behavior of cTR in
Figure 6, where even with clean data the accuracy is already quite low, and
the contamination slightly worsens the performances of the classifier. Again
the subspaces of rFDA and rTR are close to the theoretical ones, except for
the contamination level $\varepsilon=0.3$.

Figure 6: Distribution of the accuracy over 200 simulation runs of Scenario
III. The accuracy of tTR, tFDA, and tQDA is reported as benchmark in the
leftmost plot. Classical/robust TR and FDA classifiers, and AdaBoost, are
considered under several levels of contamination.

The fourth scenario ($\ell=\mbox{\scriptsize{\rm{IV}}}$) combines the
distribution of scenario II with $q\in\\{0,10,20,50,100\\}$ additional
variables, which are irrelevant for the classification problem. Under this
scenario, in the training set we consider two group sample sizes
$n_{j}\in\\{40,\,400\\}$, clean data ($\varepsilon=0$) and 10% of
contamination ($\varepsilon=0.10$). We compare only cFDA, cTR, rFDA, and rTR
for different numbers of irrelevant variables. The median accuracy is shown in
Figure 7 for $n_{j}=400$ and Figure 8 for $n_{j}=40$. In the latter case, when
$p>n_{j}$, the MRCD estimators are needed.

When $n_{j}=400$, Figure 7 shows that the number of irrelevant variables has a
small impact on the accuracy, if the data is clean ($\varepsilon=0$). As
expected, the contamination affects mainly cFDA and cTR; cFDA median accuracy
seems to decrease faster than the one of cTR, but for every $q$ the median
accuracy of cTR is lower than the one of cFDA.

Figure 7: Median accuracy of 200 simulations in Scenario IV, for different
numbers of irrelevant variables ($q$), with $n_{j}=400$. Figure 8: Median
accuracy of 200 simulations in Scenario IV, for different numbers of
irrelevant variables ($q$), with $n_{j}=40$.

The most relevant finding is the impact of the sample size. Even if both tTR
and tFDA find the right theoretical matrices of coefficients (see (21)),
smaller sample sizes ($n_{j}=40$) lead to poorer classification results when
irrelevant variables exist, regardless the presence of contamination. In the
uncontaminated case, the median accuracy for cTR (cFDA) is 0.881 (0.888) when
$q=0$, and decreases up to 0.603 (0.569) when $q=100$; in contrast, when
$n_{j}=400$, we go from 0.881 to 0.869 for both cFDA and cTR. When
contamination is present and $n_{j}=40$, robust estimators achieve a similar
accuracy of their classical counterparts, when no contamination is present.
For classical estimators we also observe a faster decrease in median accuracy
as the number of irrelevant variables increases.

The pattern for the classical estimators under contamination is similar to the
one of Figure 7, but characterized by lower median accuracies. cFDA median
accuracy becomes lower than the one observed for cTR when $q\geq 25$, and
decreases much faster, meaning that cFDA is clearly more sensitive to the
increase of irrelevant variables than cTR. The degradation in the median
accuracies of the robust estimators is quite similar to the one registered for
the classical counterparts when $\varepsilon=0$, as expected.

So far we have discussed three scenarios in which TR and FDA either perform
equally well, or show some differences. Interestingly, scenario III highlights
the idea that a nonlinear and more complicated model is not always better than
a linear one in terms of accuracy. As the contamination level increases, in
scenarios II and III the accuracy of cFDA decreases gradually; in all three
scenarios, cTR is immediately affected by the presence of contamination but
its performance does not deteriorate as the contamination level increases.
Even if FDA seems to be more resistant to outliers, it may end up being worst
than TR for high contamination levels. Robustification helps both methods up
to approximately $25\,\%$ of contamination, by construction, leading to the
same accuracies of cFDA and cTR when $\varepsilon=0$. In scenario IV, the
increasing number of irrelevant variables affects the accuracy of both TR and
FDA, regardless the chosen estimators for the covariance matrices. This
behavior is more pronounced when the group sizes are small compared to the
number of irrelevant variables.

In the next section, we present FDA and TR as linear dimensionality reduction
strategies. We study their performances on real datasets.

##  6 Real datasets

When TR and FDA are regarded as dimensionality reduction techniques, the
projection step can be viewed as a pre-processing step, and a classifier is
built on top of the projected data. We choose rLDA as classifier, but other
alternatives (e.g., logistic regression, classification trees, etc.) are also
possible.

Table 3: Real datasets, with size ($n\times p$), ranges of group sizes ($n_{j}$), numbers of groups ($g$), and the considered maximum reduced dimension $r=\min\\{p,\,g-1\\}$; $p$ represents the number of variables with $Q_{n}>0$. Dataset |  Size $n\times p$ |  Range $n_{j}$ | $g$ | $r$ | Complete name | Source
---|---|---|---|---|---|---
ecoli | $\phantom{1}327\times\phantom{1}5$ | $[\phantom{1}20,\,143]$ | 5 |  14 | Ecoli | [14]
floral buds | $\phantom{1}550\times\phantom{1}6$ | $[\phantom{1}44,\,363]$ | 4 |  14 | Floral Buds | [36]
libras | $\phantom{1}360\times 90$ | $[\phantom{1}24,\,\phantom{1}24]$ | 15 | 14 | Movement Libras | [14]
mfeat | $2000\times 64$ | $[200,\,200]$ | 10 |  19 | Multiple Features (Karhunen-Love coefs.) | [14]
penguins | $\phantom{1}333\times\phantom{1}4$ | $[\phantom{1}34,\,\phantom{1}73]$ | 6 |  13 | Palmer Penguins | [37]
plant margin | $1600\times 53$ | $[\phantom{1}16,\,\phantom{1}16]$ | 100 | 53 | One-hundred plant species leaves (margin) | [38]
plant shape | $1600\times 64$ | $[\phantom{1}16,\,\phantom{1}16]$ | 100 | 64 | One-hundred plant species leaves (shape) | [38]
plant texture | $1599\times 49$ | $[\phantom{1}15,\,\phantom{1}16]$ | 100 | 49 | One-hundred plant species leaves (texture) | [38]
segment | $2310\times 16$ | $[330,\,330]$ | 7 |  16 | Image segmentation | [14]
texture | $5500\times 40$ | $[500,\,500]$ | 11 | 10 | Texture | [15]

Experiments are conducted on the real datasets listed in Table 3. For each
dataset, we discard observations containing missing values, and remove those
variables for which the $Q_{n}$ is zero (we recall that $Q_{n}$ [39] is a
high-breakdown point robust estimate of the scale of a random variable). This
step aims at removing those variables that are almost constant on the entire
dataset. For penguins, we combine sex and species to get 6 classes, while we
discard the information about the penguins’ home islands. In ecoli, we do not
consider classes with less than 20 observations.

The experiments consist of two steps: first, we reduce the dimension of the
problem from $n\times p$ to $n\times k$, for $k=1,\dots,r$, via either
classical or robust TR and FDA. Secondly, we use the projected data to build
the robust LDA (rLDA) classifier: this is equivalent to computing the robust
Mahalanobis distances between points and means in the $k$-dimensional subspace
spanned by the chosen dimensionality reduction method, and assigning each data
point to its closest group (we refer to Section 4 for more details).

We compute the accuracy to evaluate the performances of rLDA combined with the
various dimensionality reduction strategies. Each dataset is split into a
training set, $\mathbf{X}_{\rm train}$, and a test set, $\mathbf{X}_{\rm
test}$. The training set is used to compute the projection matrix,
$\mathbf{V}$, and to estimate the projected group means and pooled covariance
matrix for rLDA. The test set data points are then classified according to
rLDA. The median accuracy is estimated via 10-fold cross-validation. We remark
that each fold contains the classes in the same proportions as the original
dataset.

When estimating the covariance matrices, we may have the following problems:
the (numerical) singularity of $\mathbf{W}$ estimates, and the singularity of
one or more individual group covariance matrices. The first issue may be
caused by possibly high correlations among the variables. The nullspace of
$\mathbf{W}$ may intersect the nullspace of $\mathbf{B}$. In FDA (which is a
generalized eigenvalue problem) this fact leads to the so-called singular
pencils. In TR we may lose the uniqueness of the solution, because for some
reduced dimension $k$ there may be no gap between the $k$th and the $(k+1)$st
eigenvalue of $\mathbf{B}-\rho^{(k)}\mathbf{W}$. To avoid this situation, we
compute the eigenvectors of $\mathbf{W}$ and project both $\mathbf{W}$ and
$\mathbf{B}$ onto the subspace spanned by the range of $\mathbf{W}$, i.e., by
the eigenvectors corresponding to its positive eigenvalues. Numerically, we
select only the eigenvalues of $\mathbf{W}$ satisfying
$\lambda_{i}(\mathbf{W})>10^{-10}\,\lambda_{1}(\mathbf{W})$ and the
corresponding eigenspace ${\bf\Gamma}$. We solve TR and FDA for the pencil
projected onto the span of the eigenvectors
$({\bf\Gamma}^{T}\mathbf{B}{\bf\Gamma},\,{\bf\Gamma}^{T}\mathbf{W}{\bf\Gamma})$
and thus find a projection matrix $\mathbf{V}$ with fewer rows. In the end,
this is equivalent to work with $\mathbf{X}_{\rm train}{\bf\Gamma}$ as input
of the training task; test data is projected accordingly. We adopt the same
strategy for the robust covariance matrices $\widetilde{\mathbf{B}}$ and
$\widetilde{\mathbf{W}}$.

The singularity of the individual group covariance matrices can be related to
the fact that, for some $j$’s, the group size is smaller than the number of
variables $n_{j}<p$. This is not problematic for the classical estimates
$\mathbf{S}_{j}$, and, in the worst-case scenario, it results in a singular
within covariance matrix, which can be tackled as we have already suggested.
As mentioned in Section 4, in the robust case we use the MRCD estimates [13]
for mean and covariance, in place of the MCD estimates. These prevent the
determinant of $\widetilde{\mathbf{S}}_{j}$ from being zero.

When $n_{j}\ll p$, the MRCD estimates could give very different results
depending on the selected subset of points. Therefore, we have repeated the
same experiments with a different regularization strategy for dimensionality
reduction. This consists in centering the groups with respect to their
$L_{1}$-median (see, e.g., [40]) and then computing the MCD estimate of the
common covariance matrix. We note that this approach implicitly assumes that
all groups share the same covariance matrix. In general, we did not see any
improvement compared to the use of MRCD estimates. In some cases, the
performance of this other strategy is even less satisfactory.

We mention in passing that we have also compared the rLDA classifier with the
classical LDA method. As expected, rLDA results in equal or a better
classification performance, except for the plant shape set. The uncommon
behavior of plant shape is also noted when using rLDA and will be discussed
later.

In Figure 9, we report the median accuracy of each method as a function of the
reduced dimension $k$. The results are also compared with rLDA where no prior
dimensionality reduction is performed. The rLDA median accuracy is also
computed using a 10-fold cross-validation strategy.

Figure 9: Median accuracy for each classification problem, as a function of
the reduced dimension $k$. Methods that use robust estimates (including rLDA)
are indicated with solid lines; dashed lines are for classical estimates.

In all problems, FDA with classical or robust estimates reaches the median
accuracy of rLDA (on the original data); this means that we can reduce the
dimension of the problem without losing important information for the
classification problem. In addition, the lower-dimensional representation of
the data enhances interpretability and eases further possible analyses. TR
does not always reach the same accuracy, and the increase in the performance
looks slower than FDA. In some cases, e.g., libras, mfeat and segment, it
seems that a higher dimension than $k=r$ is needed to get to the same accuracy
as rLDA and FDA. Finally, there are some situations where FDA and TR give a
higher median accuracy than the complete model: cTR and cFDA in ecoli; cFDA,
rFDA, and cTR in floral buds; rFDA in libras and cFDA in segment.

In general, FDA seems to perform as good as or better than TR, with two small
exceptions: ecoli for $k\geq 3$, penguins for $k=2$. In the latter case, FDA
is still better when $k>2$. While robustification seems to be of little help
for FDA, with the exception of libras, we remark that rTR can perform much
better than cTR. This is the case of libras, mfeat, plant shape, and Texture.

There is an interesting and uncommon behavior in plant shape. For large values
of $k$, the median accuracy of rFDA and cFDA decreases, while it keeps
increasing for rTR and cTR. One would expect that when the reduced dimension
is actually the original size of the problem, i.e., for $k=p$, rFDA coincides
with rLDA. This might not be the case for two reasons: $\mathbf{V}_{\rm
FDA}\in\mathbb{R}^{p\times p}$ is not an orthogonal basis for $\mathbb{R}^{p}$
and thus it does not preserve distances. Secondly, the MRCD covariance matrix
is not affine equivariant [13], meaning that the estimator of the variance of
the projected data is not necessarily equal to $\mathbf{V}_{\rm
FDA}^{T}\widetilde{\mathbf{W}}\mathbf{V}_{\rm FDA}$. We also point out that
this is the dataset where the largest difference between FDA and TR is
observed: while FDA reaches its maximum median accuracy around $k=10$, TR
needs $k=p$ to compete with rLDA.

##  7 Conclusions

We have studied and tested FDA and TR, two linear dimensionality reduction
methods that are based on the within and between covariance matrices. We have
shown that, under a particular contamination scenario, the two contaminated
covariance matrices can be written as the sum of the clean covariance
matrices, and a perturbation term that depends on the contamination level of
the data. This allowed us to reinterpret the asymptotic perturbation bounds
for the solution to TR [27] in terms of the parameters of the original
variables, the contaminating variables, and the probability of contamination.

We have introduced a new robust version of TR, by running its algorithm with
the (regularized) MCD estimates [13] of the between and within covariance
matrices. In the simulation study of Section 5, the performance of the method
has been compared with its classical counterpart, FDA, QDA, and AdaBoost. In
all scenarios, robust estimates have a good impact on both TR and FDA
performances. In some cases, FDA shows a higher accuracy than AdaBoost: this
reinforces the idea that a nonlinear model does not necessarily have a higher
accuracy than a linear one. We have also studied FDA and TR for an increasing
number of irrelevant variables, and noticed that the accuracy of the
classification decays for both FDA and TR, even when the covariance matrices
are estimated robustly. This decay is faster when the group sizes are much
smaller than the number of irrelevant variables.

Finally, we have assessed the performance of FDA and TR as dimensionality
reduction methods on real datasets, in combination with robust LDA as
classifier. In general, FDA performs as well as or better than TR; in most
problems, both methods reach the performance of the full model, where the
classifier is trained without any dimensionality reduction; in a few cases,
they perform even slightly better. For the tested real datasets, the robust
estimators have little impact on the accuracy of FDA, while, in general, a big
improvement is observed in the TR method.

Acknowledgments This work has received funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Sklodowska-
Curie grant agreement No 812912. It has also received support from Fundação
para a Ciência e Tecnologia, Portugal, through the project UIDB/04621/2020.

## References

  *   * Johnson and Wichern [2007] Johnson, R.A., Wichern, D.W.: Applied Multivariate Statistical Analysis, 6th edn. Prentice Hall, Upper Saddle River, NJ, USA (2007)
  * Ngo et al. [2012] Ngo, T.T., Bellalij, M., Saad, Y.: The trace ratio optimization problem. SIAM Rev. 54(3), 545–569 (2012)
  * Todorov and Pires [2007] Todorov, V., Pires, A.M.: Comparative performance of several robust linear discriminant analysis methods. Revstat Stat. J. 5(1), 63–83 (2007)
  * Hawkins and McLachlan [1997] Hawkins, D.M., McLachlan, G.J.: High-breakdown linear discriminant analysis. J. Am. Stat. Assoc. 92(437), 136–143 (1997)
  * Croux and Dehon [2001] Croux, C., Dehon, C.: Robust linear discriminant analysis using S-estimators. Can. J. Stat. 29(3), 473–493 (2001)
  * Hubert and Van Driessen [2004] Hubert, M., Van Driessen, K.: Fast and robust discriminant analysis. Comput. Stat. Data Anal. 45(2), 301–320 (2004)
  * Bashir and Carter [2005] Bashir, S., Carter, E.M.: High breakdown mixture discriminant analysis. J. Multivar. Anal. 93(1), 102–111 (2005)
  * Vanden Branden and Hubert [2005] Vanden Branden, K., Hubert, M.: Robust classification in high dimensions based on the SIMCA method. Chemometr. Intell. Lab. Syst. 79(1-2), 10–21 (2005)
  * Wold [1976] Wold, S.: Pattern recognition by means of disjoint principal components models. Pattern Recognit. 8(3), 127–139 (1976)
  * Pires and Branco [2010] Pires, A.M., Branco, J.A.: Projection-pursuit approach to robust linear discriminant analysis. J. Multivar. Anal. 101(10), 2464–2485 (2010)
  * Hoffmann et al. [2016] Hoffmann, I., Filzmoser, P., Serneels, S., Varmuza, K.: Sparse and robust PLS for binary classification. J. Chemom. 30(4), 153–162 (2016)
  * Ortner et al. [2020] Ortner, I., Filzmoser, P., Croux, C.: Robust and sparse multigroup classification by the optimal scoring approach. Data Min. Knowl. Discov. 34(3), 723–741 (2020)
  * Boudt et al. [2020] Boudt, K., Rousseeuw, P.J., Vanduffel, S., Verdonck, T.: The minimum regularized covariance determinant estimator. Stat. Comput. 30(1), 113–128 (2020)
  * Dua and Graff [2017] Dua, D., Graff, C.: UCI Machine Learning Repository (2017). http://archive.ics.uci.edu/ml
  * Derrac et al. [2015] Derrac, J., Garcia, S., Sanchez, L., Herrera, F.: KEEL data-mining software tool: data set repository, integration of algorithms and experimental analysis framework. J. Mult.-valued Log. Soft. Comput. 17 (2015)
  * Knyazev and Argentati [2002] Knyazev, A.V., Argentati, M.E.: Principal angles between subspaces in an A-based scalar product: algorithms and perturbation estimates. SIAM J. Sci. Comp. 23(6), 2008–2040 (2002)
  * Fisher [1936] Fisher, R.A.: The use of multiple measurements in taxonomic problems. Annals of eugenics 7(2), 179–188 (1936)
  * Absil et al. [2009] Absil, P.A., Mahony, R., Sepulchre, R.: Optimization Algorithms on Matrix Manifolds. Princeton University Press, Princeton, NJ, USA (2009)
  * Hastie and Tibshirani [1996] Hastie, T., Tibshirani, R.: Discriminant analysis by Gaussian mixtures. J. R. Stat. Soc.: Series B (Methodological) 58(1), 155–176 (1996)
  * Fukunaga [1990] Fukunaga, K.: Introduction to Statistical Pattern Recognition, 2nd edn. Academic Press Professional, Inc., Cambridge, MA, USA (1990)
  * Park and Park [2008] Park, C.H., Park, H.: A comparison of generalized linear discriminant analysis algorithms. Pattern Recognit. 41, 1083–1097 (2008)
  * Shi and Wu [2021] Shi, W., Wu, G.: New algorithms for trace-ratio problem with application to high-dimension and large-sample data dimensionality reduction. Mach. Learn., 1–28 (2021)
  * Jia et al. [2009] Jia, Y., Nie, F., Zhang, C.: Trace ratio problem revisited. IEEE Trans. Neural Netw. 20, 729–735 (2009)
  * Wang et al. [2007] Wang, H., Yan, S., Xu, D., Huang, X.: Trace-ratio vs. ratio-trace for dimensionality reduction. Proc. IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognit., 1–8 (2007)
  * Zhao et al. [2013] Zhao, M., Chan, R.H., Tang, P., Chow, T.W., Wong, S.W.: Trace ratio linear discriminant analysis for medical diagnosis: A case study of dementia. IEEE Signal Process. Lett. 20(5), 431–434 (2013)
  * Zhang et al. [2010] Zhang, L.H., Liao, L.Z., Ng, M.K.: Fast algorithms for the generalized Foley–Sammon discriminant analysis. SIAM J. Matrix Anal. Appl. 31(4), 1584–1605 (2010)
  * Zhang and Yang [2013] Zhang, L.H., Yang, W.H.: Perturbation analysis for the trace quotient problem. Linear Multilinear Algebra 61(12), 1629–1640 (2013)
  * Rousseeuw and Driessen [1999] Rousseeuw, P.J., Driessen, K.V.: A fast algorithm for the minimum covariance determinant estimator. Technometrics 41(3), 212–223 (1999)
  * Hastie et al. [2009] Hastie, T., Tibshirani, R., Friedman, J.: The Elements of Statistical Learning. Springer Series in Statistics. Springer, New York, NY, USA (2009)
  * Hasselman [2019] Hasselman, B.: Geigen: Calculate Generalized Eigenvalues of a Matrix Pair. (2019). R package 1.5. http://CRAN.R-project.org/package=geigen
  * Golub and Van Loan [2013] Golub, G.H., Van Loan, C.F.: Matrix Computations. Johns Hopkins University Press, Baltimore, MD, USA (2013)
  * Todorov and Filzmoser [2009] Todorov, V., Filzmoser, P.: An object-oriented framework for robust multivariate analysis. J. Stat. Softw. 32(3), 1–47 (2009)
  * Freund and Schapire [1997] Freund, Y., Schapire, R.E.: A decision-theoretic generalization of on-line learning and an application to boosting. J. Comput. Syst. Sci. 55(1), 119–139 (1997)
  * Benavoli et al. [2016] Benavoli, A., Corani, G., Mangili, F.: Should we really use post-hoc tests based on mean-ranks? J. Mach. Learn. Res. 17(5), 1–10 (2016)
  * Demšar [2006] Demšar, J.: Statistical comparisons of classifiers over multiple data sets. J. Mach. Learn. Res. 7(1), 1–30 (2006)
  * Raymaekers and Rousseeuw [2022] Raymaekers, J., Rousseeuw, P.: Classmap: Visualizing Classification Results. (2022). R package version 1.2.0. https://CRAN.R-project.org/package=classmap
  * Horst et al. [2020] Horst, A.M., Hill, A.P., Gorman, K.B.: Palmerpenguins: Palmer Archipelago (Antarctica) Penguin data. (2020). R package version 0.1.0. https://allisonhorst.github.io/palmerpenguins/
  * Mallah et al. [2013] Mallah, C., Cope, J., Orwell, J.: Plant leaf classification using probabilistic integration of shape, texture and margin features. Signal Process. Pattern Recognit. Appl. 5(1), 45–54 (2013)
  * Rousseeuw and Croux [1993] Rousseeuw, P.J., Croux, C.: Alternatives to the median absolute deviation. J. Am. Stat. Assoc. 88(424), 1273–1283 (1993)
  * Fritz et al. [2012] Fritz, H., Filzmoser, P., Croux, C.: A comparison of algorithms for the multivariate ${L}_{1}$-median. Comput. Stat. 27, 393–410 (2012)