The effect of intrinsic dimension on the Bayes-error of projected quadratic discriminant classification

High-dimensionality is a common hurdle in machine learning and pattern classification; mitigating its effects has attracted extensive research efforts. It has been found in a recent NeurIPS paper that, when the data possesses a low effective dimension, the predictive performance of a discriminative quadratic classifier with nuclear norm regularisation enjoys a reduced (logarithmic) dependence on the ambient dimension and depends on the effective dimension instead, while other regularisers are insensitive to the effective dimension. In this paper, we show that dependence on the effective dimension is also exhibited by the Bayes error of the generative Quadratic Discriminant Analysis (QDA) classifier, without any explicit regularisation, under three linear dimensionality reduction schemes. Specifically, we derive upper bounds on the Bayes error of QDA, which adapt to the effective dimension, and entirely bypass any dependence on the ambient dimension. Our findings complement previous results on compressive QDA that were obtained under compressive sensing type assumptions on the covariance structure. In contrast, our bounds make no a-priori assumptions on the covariance structure, in turn they tighten in the presence of benign traits of the covariance. We corroborate our findings with numerical experiments.


Introduction
Modern data sets often consist of large numbers of features. When the data is high dimensional, the "curse of dimensionality" typically degrades the performance of machine learning algorithms, as learning in high dimensions requires much larger training sample sizes and increased computational resources [1]. It has been observed however that many real-life data sets do not fill their ambient spaces evenly [2,3], in which case, learning form them is expected to be easier both statistically and computationally. To distinguish this structural notion of inherent dimension of data from the ambient dimension, we shall refer to it, in a general sense, as the intrinsic dimension (ID). We thus refer to such data sources as possessing a low-ID.
B Efstratios Palias exp093@bham.ac.uk Ata Kabán a.kaban@bham.ac.uk 1 A large volume of research has been dedicated to mitigating the ill effects of high dimensionality, including regularisation methods, and dimensionality reduction approaches. There is also an increasing number of studies aimed at elucidating whether some definition of ID allows learning algorithms to succeed with less resources [3][4][5][6][7], mostly focusing on complex models.
At the confluence of regularisation-based approaches and intrinsic dimension-based analyses, a recent study at NeurIPS [8] examined homogeneous discriminative quadratic classifiers, and proved that, whenever the data has a low effective dimension (a notion of ID), a nuclear norm constraint of its matrix parameter enables it to adapt to the effective dimension of the data, and have only logaritmic dependence on the ambient dimension. In other words, the classifier obtained needs less training data to achieve good generalisation whenever the data has a low effective dimension. The authors also confirmed this effect experimentally, and demonstrated that a number of other norm-constraint based regularisation schemes lack such desirable property.
Motivated by these recent results, our goal in this paper is to investigate the effect of intrinsic or effective dimension in QDA, without explicit regularisation, under dimensionality reduction schemes.
Dimensionality reduction essentially means to extract or to construct a smaller number of features than the number of observed or measured ones. Its aim is to improve computational efficiency, while potentially sacrificing the accuracy a little, or even enhancing accuracy by reducing noise. Some methods examine the data set to find its key features, others perform a general purpose compression. Linear projections are popular approaches, whereby the original high dimensional data is linearly projected into a lower dimensional subspace before pursuing a learning task. Examples include Random Projections (RP), and Principal Components Analysis (PCA).
Our interest in dimension-reducing projections stems from the fact that, for some algorithms, they can reduce or eliminate the ambient dimension dependence of their accuracy. For instance, for the popular classification algorithm of Linear Discriminant Analysis (LDA), it has been found that, under mild assumptions, the generalisation error of its randomly projected version does not depend on the ambient dimension [9,10]. In addition, a significant advantage of RP is that it can take advantage of low-ID without the need to know or compute the ID of the data.
LDA is special in that its Bayes-error has a closed form. It is natural to ask if a similar behaviour can be found in the more general classification algorithm of Quadratic Discriminant Analysis (QDA). This is what we set out to study in this paper. The importance of our work stems from the popularity of QDA as a classification algorithm, with many applications throughout several decades, including genetics [11,12], medicine [13,14] and image processing [15,16].
The remainder of this paper is structured as follows: In Sect. 2, we formally introduce the model of QDA. We then review the related literature and state our contributions. In Sect. 3, we present and discuss our results and findings. Section 4 is devoted to series of controlled experiments on synthetic data, designed to verify whether the bounds we derived reveal something real about the effect of a low-ID on the generalisation error of QDA. Section 5 concludes the paper and outlines possible future research directions.

Background
Linear and Quadratic Discriminant Analysis (LDA and QDA) were introduced by Fisher in 1936 [17] as classification algorithms for two or more classes of instances. These still represent popular methods in statistical pattern recognition due to their effectiveness, simplicity and interpretability [18,19]. Due to their success, they have been extensively studied over the past decades. However, the bulk of analytic work is focused on LDA, as its Bayes-error has a closed-form expression. No closed form is available for the Bayes error of QDA, its analysis is harder, and existing results are much more scarce. In this work we focus on QDA.

The QDA model
Quadratic Discriminant Analysis (QDA) is a generative classifier. Let X × Y be the instance space, where X is the feature space and Y = {0, 1} is the set of labels. We will write each instance as a pair (x, y), where x ∈ X is the feature vector and y ∈ Y is its label. Let C 0 and C 1 denote the two data classes. QDA assumes class-conditional Gaussian models on both classes. That is, the data instances are assumed to be generated i.i.d. as follows.
Here, N (μ, ) denotes the multi-variate Gaussian distribution with mean vector μ and covariance matrix . LDA assumes further that 0 = 1 , while QDA does not make this assumption. Thus, the parameters of QDA are the following: In practice, the true parameters are not available, so we usually compute estimates from a finite sample or training set. However, this paper is concerned with the Bayes-error of QDA, which is a theoretical quantity representing the expected misclassification error as a function of the true parameters. We will therefore assume that the true parameters in (1) are known.
For a given feature vector x ∈ X , and the true parameter values, the prediction rule of QDA is a function f : X → Y defined as where p(x|μ, ) is the density of N (μ, ) evaluated at the point x. Expanding out the definition of the Gaussian density it is easy to see that this prediction rule takes the form of a quadratic decision boundary in general. The model reduces to LDA when 0 = 1 , in which case the decision boundary is linear [20].
The generalisation error of any classifier is the probability of misclassification. For QDA, this is where D := π 0 · N (μ 0 , 0 ) + π 1 · N (μ 1 , 1 ) is the data generating distribution under the QDA model, I{·} is the indicator function, which equals 1 if its argument is true, and 0 otherwise, and f is defined as in (2). This is the probability that a future instance from either class is incorrectly classified. The ultimate goal of any machine learning algorithm to minimize this quantity. Under the QDA model, there is a well-known upper bound on the Bayes error of QDA, the Bhattacharyya bound, which will be our starting point in this paper. To simplify the setting, we will consider the common case where the class priors π 0 and π 1 are equal.
Then, an upper bound for the generalisation error of (3) is given by The expression (4) is called the Bhattacharyya distance (hence the name of the bound), as it is a measure of statistical distance between two Gaussian distributions, being zero if and only if the two distributions have the same parameters [21].

Related work
Our work is focused on the generalisation error of QDA under projection based dimensionality reduction. We want to find out whether it adapts to some notion of low-ID in the data. This question was previously studied in the simpler setting of LDA [9,10], where the authors give high-probability upper bounds on the error of LDA under Gaussian random projection. In the finite dimensional case, if the true covariance matrix has a low rank, they find that the ambient dimension is replaced with this rank. Moreover, in the infinite-dimensional case, if the covariance has a finite trace, the dimension is replaced with another notion of ID known as the effective dimension. This is the ratio of the trace of the covariance matrix and its largest eigenvalue. Interestingly, the same notion of ID appears in the analysis of nuclear-norm regularised discriminative quadratic classifier in the recent study of [8]. This particular notion of ID will also appear in our results, so we will define it formally in Sect. 4.
As we already mentioned, similar results are sparse for QDA. The recent study at Neurips [8] considered a specific form of quadratic classifiers, namely homogeneous discriminative quadratic classifiers of the form x T Ax ≷ 0, where A can be an indefinite matrix parameter, and is subjected to a nuclear norm regularisation constraint. Their results improve on the then-best known bounds from [22] by reducing the dependence on the ambient dimension to logarithmic, and showing the ability of the nuclear norm regulariser to adapt to the intrinsic dimension. They also demonstrate experimentally that several other regularisers, such as the Frobenius norm regulariser, lack this ability.
The main analytic tool used in [8] is the Rademacher complexity, which has the advantage of distribution-freeness (no class-conditional Gaussian assumption). However, they do not consider the full generative QDA model, and do not consider any dimensionality reduction method, but only the discriminative classifier that operates in the full ambient space. Moreover, their bound still has a logarithmic dependence on the ambient dimension. It is therefore of interest to us to consider the full generative QDA model, and consider dimensionality reduction rather than regularisation. In QDA, the matrix that appears in the decision boundary can be positive-or negative-definite, and the quadratic expression in the decision boundary is not necessarily homogeneous.
Another thread of work similar to ours appeared in [23], where a multi-class Gaussian classifier very similar to QDA is analysed under randomly projected measurement data. However, the authors assume that the original high dimensional covariances are low-rank, and have some Gaussian noise added after being randomly projected. Furthermore, the covariance matrices are assumed to be known and span random subspaces of equal dimensions. These subspaces are assumed to have been sampled randomly uniformly, in order to guarantee that the average of two covariance matrices has double the rank of the individual matrices almost surely. These assumptions originate from the field of compressed sensing, and help obtaining results though techniques inspired from that field.
Among other results, they find that, under their assumptions, the minimum number of measurements to guarantee a reliable classification performance depends on the rank of the covariance matrices instead of the ambient dimension. This is another step towards quantifying the dependence of classification on the ID. However, the restrictions on the covariance matrices being low-rank and also of equal ranks are rather unrealistic, and most of their results degrade rapidly when the matrices are of "almost" low-rank. The equal-and zeromeans case, in the ambient space is also quite restrictive for real-world data.
The imposition of such strong assumptions means that the results only hold when the assumptions hold. Instead, we aim to forego such a-priori structural assumptions where possible, and let our bounds reveal benign structural traits where the bounds tighten. Interestingly, we shall see good agreement between the benign structures identified this way from our bounds and those assumed a-priori in [23]; however our bounds are still valid without those a-priori assumptions.
We consider three dimensionality-reducing projections owing to their popularity, simplicity, and their diversity on how they operate on the data. The first two are random Gaussian and random orthogonal projections; these both are oblivious to the parameters and data, while the third one is principal component analysis (PCA), a deterministic function of the data parameters. To the best of our knowledge, there are no results on the error of QDA under projections that would completely eliminate dependence on the ambient dimension without a strong low-rank assumption. We are also unaware of work quantifying the effect of PCA on the error of QDA.
Furthermore, it is well known that QDA in a large ambient space suffers from the curse of dimensionality. It is therefore imperative to study its classification performance under several different dimensionality reduction schemes, and unveil its degree of dependence from the ambient dimension, at least theoretically. This would shed light and improve our understanding of how we can expect the algorithm to perform when faced with high-dimensional data sets, especially when they possess a low-ID. Motivated by this goal, we revisit QDA under projections, without imposing any assumptions other than the QDA model itself.

Our contributions
We summarise our contributions as follows.
1. We upper bound the Bayes error of QDA under three popular projection schemes: Gaussian random projection; random orthogonal projection; and principal components analysis.
• Our bounds eliminate the dependence on the ambient dimension and replace it with a specific notion of ID, i.e. the effective dimension. • Our bounds highlight some benign structural traits of the problem, including favourable covariance structures, without imposing any such structures a-priori.
2. We corroborate our theoretical findings with controlled numerical experiments on synthetic data sets, to confirm that our bounds reflect the behaviour of observed test error of QDA under these projections.

Notation 1 We denote scalars and vectors with lowercase letters and matrices with uppercase letters. The notations
R n and R m×n for the set of n-dimensional real vectors and the set of m ×n real matrices respectively are standard. Also, S n + denotes the set of symmetric positive-semidefinite n × n matrices. For vectors, · denotes the Euclidean norm. For matrices, | · | denotes the determinant. We let λ i (·) denote the i-th largest eigenvalue and tr(·) denote the trace of a matrix. We also use λ max (·) and λ min (·) to denote the largest and smallest eigenvalues of a matrix respectively.
Suppose the data distribution resides in R d where the ambient dimension d is too large. To reduce dimension, consider a linear projection mapping We are interested in the effect of such dimensionality reduction on the Bayes error of QDA.
Recall the Bhattacharyya bound given in Theorem 1. First, the effect of dimensionality reduction on the parameters (means and covariances) is as follows. The class means will be mapped as μ 0 → Mμ 0 and μ 1 → Mμ 1 (5) and the class covariance matrices will be mapped as Let us see what the Bhattacharyya bound becomes as a result of these mappings. To simplify notation, from this point on, let us write where δ is the Bhattacharyya distance, defined in (4). Then, the bound of Theorem 1 becomes exp(−δ M ), where This projected form of the Bhattacharyya bound is the object of our analysis. To upper bound exp(−δ M ), we need to lower bound δ M that appears in (8). We do so, under the three projection methods we discussed, each of which corresponds to a different choice for the matrix M. These are the three main results of our paper, and we devote one subsection to each.
In what follows, we denote the ambient dimension of the original data space by d, and the projection dimension with k, where k ≤ d. We only consider Euclidean vector spaces, so all norms and distances will be taken to be Euclidean. For completeness, we include in the appendix any supporting lemmas used for obtaining our results.

QDA under Gaussian random projection
The first dimensionality reduction method we consider is Gaussian random projection. Take a k × d matrix R filled with i.i.d. randomly sampled elements from N (0, 1), where k ≤ d, and pre-multiply the d-dimensional data source to reduce dimension. While not a projection in a strict geometric sense, its name of random projection (RP) is established, due to the near-isometry property of such linear mappings. This projection method has thoroughly been studied, including its effects on LDA classification, as discussed in Sect. 2.
The Bhattacharyya bound under the mapping x → Rx becomes exp(−δ R ), where δ R is given in (8). In this section, we are therefore interested in upper bounding exp(−δ R ); we do this in the following theorem.

Theorem 2 (Bhattacharyya bound under random projection)
Let R ∈ R k×d be a matrix whose elements are sampled i.i.d. from N (0, 1). Then, the following is an upper bound for From Theorem 2 we see that the error bound does not directly depend on the ambient dimension d but only through the traces of the covariance matrices relative to their largest eigenvalues in conjunction with the projection dimension k. This means that even though in the worst case the bound grows with d, whenever the data distribution does not fill the entire ambient space evenly, the error bound adapts to structure, essentially through a notion of intrinsic dimension. Indeed, the Bhattacharyya bound stays constant no matter how large d is, as long as these quantities, along with the distance between class centers stay the same. We will verify this experimentally in Sect. 4. Furthermore, from Theorem 2 we see that our bounds tightens when traces of the individual class covariances 0 , 1 get small and the trace of the pooled covariance gets large. The trace of a covariance relative to its largest eigenvalue quantifies to what extent the distribution fills the ambient space (and termed as the effective rank [24])-hence we see that one instance where our bound is tightest is when the individual covariances have lowest rank while their sum has highest rank-that is exactly the structural trait assumed in [23]. However, in contrast to that approach, our bound holds without this assumption while it highlights and lets us read off such favourable traits.

Proof of Theorem 2
Since we need to upper bound exp(−δ R ), we lower bound δ R and take the exponential of the negative of the bound we derive. To this end, we have to lower or upper bound each of the four terms in δ R separately, depending on their sign. For consistency of notations, we use the labels introduced in (8) for these four terms.
We lower bound E 1 with high probability as follows: with probability at least 1 − exp(−k 2 /4) − exp(− 2 /2λ max ( )). We used the lower bound of Theorem 1 to obtain (11), Theorem 2 to obtain (13) and the upper bound of Theorem 4 to obtain (14). Next, we lower bound E 2 with high probability: with probability at least 1 − exp(− 2 /2λ max ( )). Here we used the lower bound given in Theorem 4 to obtain (19). We upper bound E 3 with high probability: with probability at least 1 − exp(− 2 /2λ max ( 0 )). To obtain (25), we used the upper bound in Theorem 4. Finally, we can upper bound E 4 in the same way as E 3 , switching 0 to 1 , which gives with probability at least 1 − exp(− 2 /2λ max ( 1 )).
Combining the bounds on E 1 , E 2 , E 3 and E 4 by using the union bound, and plugging back into the exponential, we obtain the result stated in Theorem 2.

QDA under random orthogonal projection
Next, we consider random orthogonal projection to reduce dimensionality. This is similar to random projection but uses a random matrix that is semi-orthogonal 1 with a uniformly random orientation. It is expected that this will perform very similarly to the Gaussian random projection when d is large due to a well-known measure concentration effect by which the rows of R have nearly the same norms, and pairwise rows are near-orthogonal when d is large. However, R o allows us to derive a slightly tighter bound.
To generate R o , first we generate a k × d Gaussian RP matrix R as in Sect. 3.1, and then take R o := (R R T ) −1/2 R ∈ R k×d ; this is indeed semi-orthogonal. Now, we are interested in upper bounding exp(−δ R o ), which we do in the following theorem.
Theorem 3 (Bhattacharyya bound under random orthogonal projection) Let R ∈ R k×d be a matrix whose elements are sampled i.i.d. from N (0, 1) and define R k×d R o := (R R T ) −1/2 R. Then, the following is an upper bound for The bound in Theorem 3 is tighter than the one we obtained in the case of Gaussian random projection, due to the random projection being orthogonal and thus offering more control over the eigenvalues of the projected covariance matrices, and indeed three of the four terms are bounded with probability 1.
Most importantly, as before in the case of Theorem 2, the ambient dimension d does not appear in the error bound of Theorem 3. The error depends on it only through the traces of the covariance matrices relative to their largest eigenvalues, displaying adaptivity to intrinsic dimension without the user needing to know its value.

Proof of Theorem 3
As in the proof of Theorem 2, we first lower bound δ R o , considering each term separately, and using the labels in (8).
For the term E 1 , note that expanding out the definition of Therefore, this term admits the same lower bound as we had in Sect. 3.1, namely We lower bound E 2 using the lower side of Poincaré's separation theorem given in Theorem 3 to obtain To bound E 3 , we use the upper side of Poincaré's separation theorem (Theorem 3), Finally, we can upper bound E 4 in the same way as E 3 , by switching 0 to 1 , Plugging everything back into the exponential, and noting that the bounds on E 2 , E 3 , E 4 are deterministic, we obtain the result of Theorem 3.

QDA under principal component analysis
The third dimensionality reduction method we consider is PCA. We write the best rank-k approximation of as where U k ∈ R k×d is semi-orthogonal and k ∈ S k + is diagonal. The matrix U k is our projection matrix in this case, so we have to upper bound exp(−δ U k ). This is provided in the following theorem.
Theorem 4 (Bhattacharyya bound under PCA) Let U k ∈ R k×d be a semi-orthogonal matrix whose rows consist of the k principal eigenvectors of , and˜ k be the best rank-k approximation of . Then, the following is an upper bound for exp(−δ U k ), where˜ + k is the Moore-Penrose pseudoinverse of˜ k .
Proof Unlike previously, here the dimensionality reducing transformation U k is a deterministic function of ; however, as before, we aim to express the bound in terms of the parameters, so we work to eliminate isolated occurrences of U k . As in the proofs of Theorems 2 and 3, we lower bound δ U k , considering each term separately and using the labels in (8).
The term E 1 is can be written as follows: where the semi-orthogonality of U k was used to obtain (44).
The term E 2 can also be computed exactly as We upper bound E 3 as follows: In the last line we used Theorem 3. Swapping 0 and 1 in E 3 gives an upper bound on E 4 : Combining the bounds of E 1 , E 2 , E 3 and E 4 and plugging back, we obtain the result of Theorem 4.
Note that, due to the deterministic nature of PCA (given the parameters), Theorem 4 holds with probability 1-unlike Theorems 2 and 3, which only hold with high probability.
As another observation, we note that the second term of Theorem 4 is almost the same as the second term of Theorem 3, the only difference being that the largest k eigenvalues of appear in place of the smallest ones. This makes the bound for PCA strictly tighter than the bound for random orthogonal projection when the true means are equal. Remark 1 PCA is classic, but there may be other viable alternatives to finding the rank-k approximation of in conjunction with QDA. For instance, one could use the best rank-k approximation of either 0 or 1 , or some convex combination thereof. Another natural choice would be to find the projection that maximises the Bhattacharyya distance under the projected parameters, that is, find where U ∈ R k×d is semi-orthogonal. To maximize δ U , we need to maximize E 1 and E 2 and minimize E 3 and E 4 . Ignoring the term E 1 (assuming that it does not change significantly, or that the means are equal) we are left with maximizing the following three expressions Apart from the constants and the signs, the only major difference is the matrix that appears in these three expressions. This simultaneous optimisation problem appears to be not straightforward though, apart from the equal-means case, for which there exists a method using convex optimization [25]. In this context, in the general case, maximising only the first term from (52)-i.e. PCA-may then be interpreted or viewed as a simplification of the simultaneous optimisation problem (38). Crucially, as in Theorems 2 and 3, again the ambient dimension does not appear in the bound for the PCA-QDA combination, which means that letting d grow indefinitely, while keeping other parameters unchanged, will prevent the bound from increasing. Interestingly, here the trace of covariances does not appear either, so it is less clear whether the error actually depends on a notion of ID.

Numerical experiments
We now proceed to testing experimentally the theoretical results we obtained, in conjunction with the out-of-sample test error. As in Sect. 3, we devote one subsection to each projection method. For each method, we are interested in the following questions: 1. How is the test error affected by the intrinsic dimension of the data distribution compared to the ambient dimension? 2. Under what conditions on the true parameters does the predictive performance withstand the arbitrary increase of the ambient dimension?
One should reckon that it would be rather difficult to design the appropriate experiments for studying these questions on a purely experimental basis. However, our theoretical analysis gives us insights about these questions already, which serves as a guide for gaining further insights from experimentation. We would like to highlight this as a nice example demonstrating the practical value of theoretical analysis, and indeed a test of it. A second sense in which we put our theoretical results to a test is by assessing to what extent the variations of our upper bounds, when varying certain parameters of the problem, agree with variations in empirical estimates of the predictive error. An upper bound only tells us that the error does not exceed a certain threshold with high probability, but this still permits a lot of variation below that threshold, which we assess through the forthcoming experiments.
Regarding our two research questions above, we have insights pertaining to the second question directly, while for the first one we can only conjecture a monotonic relation between the error and the ID prior to conducting the numerical work.

Technical preliminaries
Part of our experimental setup involves increasing the ambient dimension of the data, while keeping all other quantities that appear in the bounds fixed, most of which involve the eigenvalues of the covariance matrices 0 , 1 and . In this section we give the procedures necessary to achieve this.
The trace of a covariance relative to its largest eigenvalue is known as the effective rank. This is actually a notion of ID that reflects to what extent the distribution fills the ambient space.

Definition 1 (Effective rank [24]) Let
The effective rank of A is defined as It is straightforward to show that 1 ≤ r (A) ≤ rank(A) ≤ d, for all A ∈ S d + . One advantage of the effective rank is that it is more detailed than the rank, as it considers the magnitude of all eigenvalues, rather than just counting the non-zero ones, and can also take on non-integer values. As we will soon see, controlling the parameters that appear in some of our bounds amounts to controlling the effective rank of a covariance matrix. Also, since it is a notion of ID, it makes sense to vary the effective rank (apart from the rank) of covariances to see its effects.
The second ingredient we need is the following result that allows us to control the largest eigenvalues of a matrix sum. Thus, there exists an orthogonal matrix V ∈ R d×d that "separates" these two complements, that is we have

Theorem 5 Let A, B ∈ S d + with orthogonal eigenvalue decompositions
where A 0 , B 0 ∈ S m + contain the m largest eigenvalues of A and B respectively, contain the d − m smallest eigenvalues of A and B respectively and 0 ∈ R (d−m)×m is the zero matrix. Note that neither of A 0 , A 1 , B 0 , B 1 is necessarily diagonal, although V can always be chosen to make one pair diagonal (either A 0 and A 1 , or B 0 and B 1 ). The matrix C will therefore be For the eigenvalues, we can ignore V and V T and only consider the block matrix in this case. Since C ∈ S d + , we know that m of its eigenvalues lie in the upper-left block and d − m lie in the lower-right block. To show that the largest m eigenvalues of C only depend on the largest m eigenvalues of A and B, we just need to show that they lie in the upper-left block. One way is to bound the quadratic forms of A 0 + B 0 and A 1 + B 1 . To this end, it suffices to only consider vectors in W or W ⊥ . Let x ∈ W be a unit vector. We have Now let y ∈ W ⊥ be a unit vector. We have Since λ m (A) + λ m (B) ≥ λ m+1 (A) + λ m+1 (B), the m largest eigenvalues of C lie in its upper-left block, whereas its d − m smallest eigenvalues lie in its lower-right block. Thus, its m largest eigenvalues only depend on the m largest eigenvalues of A and B, as they lie in the block that contains the block matrices A 0 and B 0 . This completes the proof.
Theorem 5 basically tells us that if we increase d by adding eigenvalues to A and B that are smaller than their respective m-th largest eigenvalues, and the condition in (55) continues to hold, then the m largest eigenvalues of their sum do not change. We will later give more details on how this result can be used in each experimental setup.

Remark 2
To obtain a matrix with a fixed rank or effective rank, we need to control its eigenvalues. While for the rank of a matrix this is straightforward, for the effective rank this is not always the case. To obtain a matrix with a fixed effective rank, we need control over both the sum of its eigenvalues and also its largest eigenvalue. While there might be more sophisticated methods of doing so, ours relies on a simple algebraic trick.
To sample the eigenvalues so that the effective rank equals r , we freely choose the largest eigenvalue λ 1 , and then we sample λ 2 , . . . , λ d such that the following two conditions are satisfied: and where (60) ensures that the effective rank equals r , which can be confirmed by comparing it with Definition 1 after doing a little algebra, and (61) ensures that all eigenvalues are nonnegative and that none of them is greater than λ 1 . Clearly we must have 1 ≤ r ≤ d for these two conditions to hold simultaneously. Sampling numbers such that they all fall in the same interval (in this case [0, λ 1 ]) and equal a fixed sum (in this case (r − 1)λ 1 ) is a well-studied problem and routines exist in several programming languages to implement it (e.g. [26]). If one later needs to change the largest eigenvalue, while keeping the effective rank the same, they can just multiply all eigenvalues by the same positive constant. This is why the initial choice of λ 1 can be arbitrary.
After sampling the eigenvalues, we let the matrix in question take the following form where V ∈ R d×d is an orthogonal matrix which is sampled randomly. The problem of sampling random orthogonal matrices has also been well-studied (e.g. [27]), and several routines exist in most programming languages. This setup makes the covariance matrices as random as possible, apart from a fixed effective rank. The form of (62) is also used to obtain a matrix with a fixed rank, after determining its eigenvalues.

Experimental setup
Each of the three projection methods we considered in the earlier sections has its own characteristics. Part of our setup is the same for all three, as we seek to answer our research questions and observe similarities and differences between the methods. We describe our setup here. The first step is to fix the true parameters (means, covariances) of the two Gaussian classes, as well as the projection dimension. Then, we sample a test set in the ambient space from one of these. We then choose the projection matrix (randomly or deterministically, depending on the method) and project the test set using this matrix. Finally, we run QDA on the projected set using the projected means and covariances and record the empirical test error fraction measured on the test sample. Exploiting symmetry, we will use a test set sampled from class C 0 , and count the proportion of this test set that was classified to C 1 .
We repeat this strategy as we vary the values of the ambient and projection dimensions, and also repeat for different choices of the covariance matrices. In one type of experiments, we vary the rank or the effective rank of 0 and 1 together, that is, they will always have equal rank and effective rank. This reduces the number of experiments for each experimental setup. In another type of experiments, we vary the ambient dimension, while controlling the ranks and effective ranks of the covariance matrices. In both cases, we need to fix the largest eigenvalues of 0 and 1 . We will set both of them to 1.
Equipped with these preliminaries and tools, we are now ready to present results obtained for each projection scheme separately. For completeness, we include the formal steps for all experiment designs we used in Appendix 2.

Numerical results with QDA under Gaussian RP
We start our experiments with Gaussian random projection. First we test how the empirical test error is affected by a low ID. To this end, we consider a setting where d remains fixed and the ID of 0 and 1 (measured by their rank and effective rank) is varied.
For each choice of the value of the ID, we randomly generate and fix the class mean parameters μ 0 , μ 1 such that they have fixed distance of 1. We then generate and fix the class covariances 0 , 1 using the procedure described in Remark 2. We also generate an out-of-sample test set of 1000 points.
We vary k, ensuring that it is always no higher than the rank (it can be higher than the effective rank). For each choice of k, we compute the test error. We report the average and standard deviation of the test errors from 100 independent random draws of the random projection matrix.
The complete steps of this set of experiments are given in Appendix 2 (Algorithm 1). The results are presented in Fig. 1. We see that, the lower the ID, in comparison with d, the lower the empirical test error. In addition, we also observe that, when we vary the rank, the empirical test error drops significantly faster than when we vary the effective rank.
In the next round of experiments we will vary the ambient dimension d. Our theoretical bound indicates that the error should not increase with d as long as other parameters in the bound, notably the effective ranks of covariance matrices, remain fixed. Our second set of experiments aims to test this behaviour on the actual empirical test errors. Therefore, we need to ensure that, as d increases, all other parameters that appear in Theorem 2 stay the same. Apart from k and , these are the following: The mean distance μ 0 − μ 1 can be controlled by resampling the means at a fixed distance for each d. To fix λ max ( 0 ) and λ max ( 1 ), we just need to re-determine the covariance matrices such that their largest eigenvalues stay the same. After that, tr( 0 ) and tr( 1 ) can be fixed by keeping the effective ranks r ( 0 ) and r ( 1 ) the same for all d, according to its definition. We therefore re-sample the eigenvalues of 0 and 1 for each d, ensuring that they satisfy this property. Recall that Remark 2 described our method of how to do this.
The most difficult parameter to control is λ max ( ). This is because, unlike the trace, λ max is not a linear function over S d + . To control it, we use the result we established in Theorem 5 applied to 0 /2 and 1 /2 in place of A and B, and choosing m = 1. In this specific case, Theorem 5 says that, as d increases (and thus more eigenvalues are added), as long as the first principal eigenvectors of 0 and 1 remain collinear, the largest eigenvalue of stays the same.
Having ensured the conditions of Theorem 5, we will create 0 and 1 to have the following form: where D 0 ∈ S d−1 + is diagonal with all eigenvalues of 0 expect its largest, D 1 ∈ S d−1 + is diagonal with all eigenvalues is orthogonal and 0 ∈ R d−1 is the zero vector. The setup ensures that the eigenvectors corresponding to λ max ( 0 ) and λ max ( 1 ) both lie along the direction of the first axis, thus satisfying the condition of Theorem 5 for m = 1. The setup of (64) comes without loss of generality, and applying the same rotation to both 0 and 1 will result in the condition still holding, albeit for another direction.
As d increases, we re-sample the eigenvalues of 0 and 1 , keeping both of their largest eigenvalues, as well as their effective ranks the same. The covariance matrices are then created anew using (64), which ensures that the conditions of Theorem 5 are still satisfied for m = 1. That is, the matrices D 0 , D 1 and V grow in size but the first principal eigenvectors of 0 and 1 remain collinear.
The steps of this set of experiments are summarised in Appendix 2 (Algorithm 1). The results are given in Fig. 2. For all values of k, we see that, as d increases, the mean empirical test error fluctuates a little for small d, but then it stabilises for larger d of around 750. Most importantly, k = 5 which corresponds to an upper bound on the error in comparison to all larger values of k, the empirical test error stays almost constant for all values of d tested. This corroborates our theoretical finding that the error does not increase with d directly but only though the effective rank; in other words, compressive QDA is suitable in arbitrary highdimensional problems as long as the effective rank is low.
The sudden drop of the empirical test error when d increases from 50 to 100 can be attributed to the fact that, since r ( 0 ) and r ( 1 ) remain fixed at 50, then at d = 50 the matrices are isotropic. Therefore, all of their eigenvalues are equal and the random projection results in a significant loss of information, regardless of the resulting low-dimensional subspace. When the dimension increases, however, some eigenvalues will have to decrease drastically to keep the effective rank fixed. Therefore, there are less directions of greatest variance, thus resulting in a much lower information loss, if the projection lands on a suitable subspace.
Finally, we remark that, intuitively, we would expect that the empirical test error will be higher when k is smaller. Our theoretical upper bound does not give us information on this, and it is therefore interesting to see from the empirical results in Fig. 2 that this is indeed the case in the examples tested.
We conclude from the results of this section that random projection is a suitable dimensionality reduction method for QDA when the underlying distribution has a low-ID covariance. The bound we derived for the generalisation error might be not particularly tight, but it captures most of the behaviour observed throughout our empirical investigations. In particular, the error does not grow with d but only through the effective rank.

Numerical results with QDA under random orthogonal projection
With random orthogonal projection-based dimensionality reduction we expect very similar results as those observed with Gaussian random projection. As argued in [9], and proved formally in [28], when d is large, the columns of R with entries i.i.d. from N (0, 1), are close to being orthonormal.
The setup will be identical to the one used in random projection; for full reproducibility all steps are summarised in Appendix 2 (Algorithm 2). The first set of experiments aims to find out whether the empirical test error adapts to a low-ID, and the results are presented in Fig. 3.
As expected, the empirical test error is almost the same as in Fig. 1 in both of their respective subfigures, for all choices of ID and all values of k. This confirms experimentally that random projection essentially plays the role of random orthogonal projection, or equivalently, the matrix R used for random projection, is close to being semi-orthogonal. Therefore, not much difference is seen when orthonormalizing R.
For the second set of experiments, we test the extent to which the empirical test error does not depend on d. As before, we would like to create the class parameters in a way to ensure that all quantities that appear in the bound of Theorem 3 are fixed. While the first factor is the same as in Theorem 2, and can thus be controlled in the same way, for the second factor we would need to control the k largest eigenvalues of 0 and 1 and the k smallest eigenvalues of . However, Theorem 5 only guarantees fixing the largest eigenvalues of and we have not found a way to control its smallest.
Nevertheless, we expect random projection to behave very similarly to orthogonal projection, therefore we instead opt to use the same setup as we used there in Sect. 4.3, using the steps of Algorithm 2. The results are presented in Fig. 4.
As expected, the results appear nearly identical to those in Fig. 2. The conclusions are therefore the same as in Sect. 4.3, namely, increasing d does not blow up the errors as long as the quantities that appear in the bound stay unchanged.
As a final remark, we should mention that, despite the seemingly different bounds in Theorems 2 and 3, the empirical performance is very similar. Of course for smaller d, this might have not been the case, but the significance of our work is in the high-dimensional settings.

Numerical results with QDA under principal components analysis
We finally experiment with PCA-based dimensionality reduction for QDA. In this set of experiments, we again inspect the dependence of the empirical test error on the ID. The setup will be the same as in Sects. 4.3 and 4.4, except for the fact that, since the projection is deterministic given the parameters, there is no need to sample random matrices. We will therefore record the empirical test error directly for each pair (k, d). This reduces the required computational resources, albeit performing the eigen-decomposition is significantly more costly than a random projection. For reproducibility, all steps of this experiment are given in Algorithm 3, and the results are presented in Fig. 5. From Fig. 5, we observe a completely different picture from the previously obtained Figs. 1 and 3. Here, when the rank is lower, the empirical test error seems to be unaffected by it. This is unlike random projection and random orthogonal projection, where the test error shows a clear dependence on the rank of the covariance matrices. However, for the effective rank, we see that a lower effective rank gives a lower empirical test error. For the second round, we need to fix all parameters in Theorem 4 as d increases. Apart from k, we need to control the following: To control the k largest eigenvalues of 0 , 1 and , Theorem 5 for m ≥ k provides a sufficient condition. Namely, we determine 0 and 1 to be of the following form: where D 0 ∈ S m + and F 0 ∈ S d−m + are diagonal with the m largest and d −m smallest eigenvalues of 0 respectively, and similarly D 1 ∈ S m + and F 1 ∈ S d−m + for 1 , V ∈ R m×m and W ∈ R (d−m)×(d−m) are orthogonal and 0 ∈ R (d−m)×m is the zero matrix. This is similar to the setting in (64), except that for each k we need the condition of Theorem 5 to hold with m ≥ k. This is to ensure that, as d increases, only specified Along with this condition, we need to fix the first k eigenvalues of the two covariance matrices. Therefore, we first create 50 eigenvalues, and then progressively add more as d increases, ensuring that all of them are smaller than the first 50. The covariance matrices are then created with (66). As d increases, the matrix V must remain fixed, in order to keep the upper-left block of unchanged, whereas W grows in size and has to be re-determined.
To control the norm in (65), we opt to fix the vector appearing in the norm. To this end, we again require that the condition of Theorem 5 holds for k ≤ m. This way, as d increases k remains the same and U k obtains more columns that consist of zeros. Therefore, we just require that as d increases and more coordinates of μ 0 − μ 1 appear, its existing coordinates do not change. In other words, we determine all 1000 coordinates of μ 0 and μ 1 from the start and unveil d of them progressively as we increase d. Note that the equal-means case is just a specific case of this.
The complete sequence of steps for running these experiments is given in Algorithm 3. We present the results in Fig. 6. From this figure we can clearly see that the empirical test error remains very similar for all d, but only when k ≤ 30. When k > 30, the empirical test error increases with d slightly, but it remains under an upper bound (e.g. that corresponds to k = 5). This is expected because the condition of Theorem 5 is satisfied for m = 30.
Comparing Figs. 2, 4 and 6, we observe that PCA tends to outperform both versions of random projection by about 15%, even when the means are closer to each other (since truncating coordinates reduces their distance for smaller d).
The bound of PCA-QDA is also numerically tighter. These results also show the some precise conditions on the data, under which PCA outperforms the other two projection methods, i.e. the conditions of Theorem 5.

Conclusions
We derived upper bounds for the generalisation error of QDA when the data is subjected to three different alternative projection schemes. Specifically, we considered random Gaussian projection, random orthogonal projection and PCA. The first two make use of random matrices to reduce dimensionality, whereas the third one is a deterministic function of the model parameters. The bounds we derived are analysing the effect of these projections on the wellknown Bhattacharyya bound. All three bounds turned out to be independent of the ambient dimension and instead depend on the effective rank -resembling an interesting dimension-adaptive behaviour recently found in the case of nuclear-norm regularised discriminative QDA [8]. We also confirmed our findings by extensive empirical simulations on synthetic data.
Future work can focus on deriving similar results for other dimensionality reduction methods (stochastic or deterministic) for QDA. The goal is to find the conditions under which the empirical test error adapts to a low-ID and the extend on which this is true.
Another possible research direction would be to examine whether other classifiers, apart from LDA and QDA, also enjoy a lower classification error in the presence of some low-ID assumption on the data generator, and to find out which notions of ID appear to capture this behaviour.

Appendix B: Experimentation details
We include here all steps involved in our experiments reported in Sect. 4.