Large-sample properties of unsupervised estimation of the linear discriminant using projection pursuit

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Large-sample properties of unsupervised estimation of the linear discriminant using projection pursuit © Authors, 2021 Published version Radojičić, Una; Nordhausen, Klaus; Virta, Joni Radojičić, U., Nordhausen, K., & Virta, J. (2021). Large-sample properties of unsupervised estimation of the linear discriminant using projection pursuit. Electronic Journal of Statistics, 15(2), 6677-6739. https://doi.org/10.1214/21-EJS1956 2021 Electronic Journal of Statistics Vol. 15 (2021) 6677–6739 ISSN: 1935-7524 https://doi.org/10.1214/21-EJS1956 Large-sample properties of unsupervised estimation of the linear discriminant using projection pursuit∗ Una Radojiˇci´c1, Klaus Nordhausen1,2and Joni Virta3 1Institute of Statistics & Mathematical Methods in Economics, Vienna University of Technology, Austria. e-mail: [email protected] 2Department of Mathematics and Statistics, University of Jyv¨askyl¨a, Finland. e-mail: [email protected] 3Department of Mathematics and Statistics, University of Turku, Finland. e-mail: [email protected] Abstract: We study the estimation of the linear discriminant with projection pursuit, a method that is unsupervised in the sense that it does not use the class labels in the estimation. Our viewpoint is asymptotic and, as our main contribution, we derive central limit theorems for estimators based on three different projection indices, skewness, kurtosis, and their convex combination. The results show that in each case the limiting covariance matrix is proportional to that of linear discriminant analysis (LDA), a supervised estimator of the discriminant. An extensive comparative study between the asymptotic variances reveals that projection pursuit gets arbitrarily close in efficiency to LDA when the distance between the groups is large enough and their proportions are reasonably balanced. Additionally, we show that consistent unsupervised estimation of the linear discriminant can be achieved also in high-dimensional regimes where the dimension grows at a suitable rate to the sample size, for example, pn=o(n1/3) is sufficient under skewness-based projection pursuit. We conclude with a real data example and a simulation study investigating the validity of the obtained asymptotic formulas for finite samples. Keywords and phrases: Clustering, Kurtosis, Skewness, Linear discriminant analysis, Projection pursuit. Received March 2021. Contents 1 Introduction................................ 6678 2 Estimationofthelineardiscriminant.................. 6681 2.1 Supervised estimation of the linear discriminant . . . . . . . . . 6682 2.2 Unsupervised estimation of the linear discriminant . . . . . . . 6683 2.3 Asymptotic comparison of the three estimators . . . . . . . . . 6685 3 Convex combination of skewness and kurtosis . . . . . . . . . . . . . 6687 3.1 Theoretical properties . . . . . . . . . . . . . . . . . . . . . . . 6687 ∗The work of Joni Virta was supported by the Academy of Finland (Grant 335077). 6677 6678 U. Radojiˇci´cetal. 3.2 Asymptotic comparisons . . . . . . . . . . . . . . . . . . . . . . 6688 4 High-dimensional projection pursuit . . . . . . . . . . . . . . . . . . 6692 5 Simulations ................................ 6694 6 Realdataexample ............................ 6698 7 Discussion................................. 6699 A EquivalentresultsforPCA ....................... 6701 B Proofsoftechnicalresults ........................ 6702 C Additionalsimulationresults ...................... 6730 Acknowledgments .............................. 6734 References................................... 6734 1. Introduction Classification and clustering are two central themes in modern data analysis and can be seen, respectively, as the supervised and unsupervised versions of the same problem: In classification, the group memberships, or labels, of the training data points are known and the objective is to use the training data to form a classification rule for future observations, some of the standard methods including, e.g., linear discriminant analysis, support vector machines, and random forests, see [20]. Whereas in clustering, no labels for the data points are known but we postulate that a reasonable grouping exists and aim to find it, with, e.g., k-means clustering or spectral clustering, see [20,64]. In this paper, we work under a clustering context and the assumption that the data admit a natural grouping but that their labels are indeed unknown to us. In their seminal work, [51] studied in this setting the use of projection pursuit (PP), a general family of methods searching for a projection direction that maximizes the value of the so-called projection index, see, e.g., [26,11,9,17] and the references therein. Namely, denoting the within-class covariance matrix by Σand the two group means by μ1,μ2,[51] established that using kurtosis as the projection index in projection pursuit allows the unsupervised estimation of the projection direction θ:= Σ−1(μ2−μ1) that is used in linear discriminant analysis to construct the optimal Bayes classifier, in the full absence of any label information. In other words, projection pursuit essentially allows conducting LDA in a unsupervised fashion to recover the subspace that optimally separates the two groups. Afterward, various clustering methods can then be applied to the projected data to conduct efficient clustering. While very interesting, the result of [51] raises a natural question regarding the efficiency of the procedure. Namely, how much does one lose by not knowing the labels and relying on projection pursuit compared to using LDA to recover the same direction θwhen the group memberships are known? This is the main question we study in the current paper, working, for simplicity, under the assumption of two-group normal mixtures. Our approach is asymptotic in nature and we perform the comparison through the limiting covariance matrices of the estimators in question. In particular, we show that the limiting covariance matrices of projection pursuit and LDA are proportional, allowing us to conduct Properties of PP-based estimators of the linear discriminant 6679 the comparisons simply through the corresponding constants of proportionality. Interestingly, the ratios of these constants depend on the model parameters only through the mixing proportion and the squared Mahalanobis distance (MD) between the group means, τ:= (μ2−μ1)Σ−1(μ2−μ1). In particular, the ratios do not depend directly on the dimension pof the data. As our second contribution, we show that projection pursuit can be used to consistently estimate the optimal direction even in the high-dimensional regime where p≡pn→∞as long as the dimension grows at a suitable rate to the sample size n. The present work can thus be seen as a continuation to the classical literature on confirmatory (or inferential) projection pursuit. However, unlike the present study that focuses on asymptotic efficiencies, the main question in confirmatory PP is typically to assess whether the found projections reflect actual structures in the population and are not simply artifacts caused by noise. For example, to assess the significance of the obtained results, [21] uses a bootstraplike procedure and compares the observed value of the projection index to the one obtained by applying the procedure to the corresponding Gaussian data. Similar problems are discussed also in, e.g., [58,46,36]. To illustrate the importance of assessing the significance of the obtained results, [14] give cautionary examples in which they show how exploratory projection pursuit can always find structures when applied to sparse, high-dimensional data. Besides the confirmatory projection pursuit, asymptotic results for general projection indices have been derived earlier also in the context of independent component analysis (ICA), see, e.g., [50,16,45,63]. In ICA, one assumes that the observed random p-vector xis an independent component (IC) model, i.e., there exists a full rank p×pmatrix Γsuch that Γ{x−E(x)}has independent components. The IC model is a rather wide family of distributions and, in particular, contains our model of choice, the multivariate normal mixture. This, apparently novel, result is given as Lemma B.1 in the Appendix Aand reveals that the multivariate normal mixture decomposes into p−1 independent standard normals and a non-Gaussian component which corresponds, up to the scale and sign, to the optimal linear discriminant projection of the data. This connection between the two models implies that the results of the current paper are intimately related to [63] who considered (in the context of ICA) the same projection indices as we do here. However, we remark that our contributions surpass those of [63] in two critical regards: 1) [63] derived only the asymptotic variances of the ICA parameters (ignoring their covariances), whereas we give the full limiting distribution of the estimated projection direction. Besides completing the asymptotic story, knowledge of the full distribution is crucial concerning the comparison of PP and LDA as it reveals that the limiting covariance matrix of PP is exactly proportional to the limiting covariance of LDA, see Theorems 1–4later on. 2) From a technical viewpoint, the derivation of the convergence rates of the estimators was in [63] left implicit and our proofs provide a rigorous treatment of this. In particular, to guarantee well-defined Taylor expansions of the objective functions, we need to establish the almost sure convergence of the projection pursuit estimates and, as far as we are aware, such results have not been given previously either in ICA or PP-literature. Finally, we note that the 6680 U. Radojiˇci´cetal. link between the two models furthermore implies that any projection index capable of recovering independent components in an IC model can, in this setting, be used to recover the optimal linear discriminant (assuming the index yields distinct values for the mixture and Gaussian noise). For example, [54] show that any increasing function of a subadditive squared dispersion measure can be used for such a purpose. While kurtosis is the most popular choice for the projection index in projection pursuit, also several alternatives are commonly used. In particular, skewness is a somewhat standard choice, see, for example, [40], and was shown in [38]to have the same property of being able to find the optimal projection direction without the label information as possessed by kurtosis. As such, we study also skewness-based projection pursuit in the current work. We also note that, prior to their use in projection pursuit, both skewness and kurtosis have a rich history as test statistics when testing for (multivariate) normality. For example, [43] first introduced the maximal kurtosis and skewness obtained by projections as test statistics when testing for multivariate normality and, using Monte Carlo methods, compared the power of the obtained tests to other common tests for normality against various alternatives, including several two-dimensional Gaussian mixtures. Distributional properties of the statistics introduced by [43]were later studied in [42]and[7], under the null hypotheses of multivariate normality and elliptical symmetry, respectively. See also the conjecture by [40] that the limiting distribution of the maximal skewness attainable by a linear combination of normal variables is skew-normal. Despite their ubiquitousness, as shown by [51,38], for both kurtosis and skewness there exist particular values of the mixing proportion under which the two indices are unable to recover the optimal projection direction (for example, skewness fails to produce a consistent estimate of θwhen the two groups have equal proportions). These drawbacks can be mitigated by combining both cumulants into a single projection index, in a form of a weighted linear combination. This combined projection index was first proposed in [29] and some of its distributional properties were discussed in [36]. Our results show that with a proper choice of weighting, a rather efficient unsupervised competitor for the LDA-based supervised estimator can be obtained with the combined index. Indeed, in the extreme case where the distance between the group means is large enough and the group sizes are reasonably balanced, projection pursuit is able to achieve efficiency arbitrarily close to LDA. As remarked in the previous paragraph, the asymptotic properties of the hybrid index have been studied also earlier, in the context of independent component analysis, in [63]. We note that despite the theoretical guarantees of projection pursuit, the most common unsupervised method for revealing clusters is still arguably PCA, see, e.g., [28]. However, it is also well known that PCA does not, in general, yield a consistent estimator of the linear discriminant direction. A standard example demonstrating this is the extreme case where the within-group covariance matrix Σis heavily concentrated on a direction orthogonal to the difference of the group means μ2−μ1. In such a case, the projections of the two group means onto the first principal component direction overlap, making clustering based Properties of PP-based estimators of the linear discriminant 6681 on the direction impossible. Hence, due to its unreliability in estimating the linear discriminant, PCA cannot truly be seen as a unsupervised estimator of the separating direction and, as such, we do not include it in the comparisons in the current paper. However, we have still included, for completeness, equivalent asymptotic results for PCA as we state for the other methods, and these are given in Appendix A. In recent years there has been a large amount of work on parameter estimation in Gaussian mixture models [68,35,33,1,24,30,23], particularly in high-dimensional settings and by the EM-algorithm, see for example [71,65,69, 70,12] and references therein. It is worth mentioning how, in general, methods for parameter estimation in Gaussian mixture models can also be used for unsupervised estimation of the linear discriminant in this setting, with the potential estimator being the plug-in estimator in which the group means and the common covariance are estimated by the method of choice. Some statistical guarantees and properties of EM-based estimators of parameters in Gaussian mixtures can be found in, e.g., [6,48], and are beyond the scope of this paper. However, under the real data example of Section 6, we compare the results obtained by the projection pursuit approach to those obtained using the EM-algorithm when estimating the optimal linear discriminant. The rest of the manuscript is organized as follows. In Section 2 we derive the asymptotic behavior of three estimators of the linear discriminant direction: LDA and kurtosisand skewness-based projection pursuit. A short comparison of the results is also presented. In Section 3, we give the corresponding results for projection pursuit based on a weighted combination of skewness and kurtosis and conduct a more extensive set of asymptotic comparisons between all considered methods. In Section 4we show that both kurtosisand skewness-based PP, as well as the PP based on the convex combination of those produce consistent estimators of the linear discriminant in the high-dimensional setting, where both the sample size nand the dimension pdiverge to infinity in a suitable ratio. Simulation studies exploring both the finite-sample performance of the methods and the applicability of our asymptotic results to practice are given in Section 5, while the performance and the applicability of presented methods to a real data example, as well as the comparison to the PCA, are given in Section 6. Finally, we conclude with some discussion in Section 7. All proofs of the technical results are postponed to Appendix B. 2. Estimation of the linear discriminant Let (Ω,F,P) be a probability space. Throughout the following, we assume that the (p+ 1)-dimensional pair (x,y) obeys the following model: y∼Ber(α1)andx|y∼N p{yμ1+(1−y)μ2,Σ},(1) for 0 <α 1<1, μ1,μ2∈Rp,μ1=μ2, and a full rank Σ∈Rp×p.Themarginal distribution of xis then the multivariate normal mixture, x∼α1Np(μ1,Σ)+α2Np(μ2,Σ), 6682 U. Radojiˇci´cetal. where α2:= 1 −α1. Under model (1), the classification of xis usually based on its projection onto the linear discriminant direction θ=Σ−1(μ2−μ1). This projection direction is optimal in the sense that the optimal Bayes classifier (having the minimal miss-classification rate out of all classifiers) depends on the data only through the projection θx, see, e.g., [44]. Our objective throughout the paper is the estimation of the standardized projection direction θ/θ(the scale of the projection direction is irrelevant, meaning that the unit length constraint is without loss of generality). As described in Section 1, we will consider two types of estimators, unsupervised ones which use only the random vector x(a sample from its distribution) in the estimation, and an supervised one which bases the estimation on the full pair (x,y). The supervised method is allowed more information in the estimation and is, naturally, expected to provide a more efficient estimator, a fact that is verified by our comparisons later on. 2.1. Supervised estimation of the linear discriminant If we have a sample (x1,y 1),...,(xn,y n) from the distribution of the full pair (x,y) available, the standard estimator of Σ−1(μ2−μ1) is the plug-in estimator (which is also its MLE, up to the scaling of the pooled covariance matrix) used in standard LDA. That is, using the notation, ¯ xn1:= 1 n i=1 yi n  i=1 yixi,¯ xn2:= 1 n i=1(1 −yi) n  i=1 (1 −yi)xi, Sn:= 1 n−2n  i=1 yi(xi−¯ x1)(xi−¯ x1)+ n  i=1 (1 −yi)(xi−¯ x2)(xi−¯ x2), we consider the estimator, wn:= S−1 n(¯ xn2−¯ xn1). Asymptotic results for LDA are very standard in the literature, see for example [4]. However, these results are usually given in the case of fixed group sizes, whereas in our model the group sizes are determined by the indicator variables y1,...,y nand are, as such, random. Hence, as far as we know, the following theorem is, if not particularly groundbreaking in its conclusions, a novel one. Theorem 1. Under model (1), we have, as n→∞, √n(wn/wn−θ/θ)Np(0,ΨU), where ΨU:= 1+βτ θ2βIp−θθ θ2Σ−1Ip−θθ θ2, β:= α1α2and τ:= (μ2−μ1)Σ−1(μ2−μ1). Properties of PP-based estimators of the linear discriminant 6683 The form of the limiting covariance matrix in Theorem 1is rather simple and inspection of the proof of the result reveals that the involved projection matrices onto the orthogonal complement of the direction θ/θare simply consequences of the standardization of the estimator to unit length. Note also that the scalar factor in front can be written as 1/(θ2β)+(θ/θ)Σ(θ/θ), the two summands of which have the following rough interpretations: If the groups are imbalanced, βis small, making the first summand large and inflating the asymptotic variance. Similarly, if the data exhibit a large amount of variation in the direction of the optimal discriminant direction, i.e., (θ/θ)Σ(θ/θ)is large, the second term increases the magnitude of the asymptotic variance. 2.2. Unsupervised estimation of the linear discriminant Kurtosis-based projection pursuit Let δ1:= 1/2−1/√12, δ2:= 1/2+1/√12 and ˜ x:= x−E(x). The kurtosis κ:Sp−1→Rof the projection of xon a given direction u∈Sp−1is then defined as, κ(u)= E{(u˜ x)4} [E{(u˜ x)2}]2. The fact that projection pursuit based on kurtosis is Fisher consistent for the linear discriminant under normal mixtures was first shown in [51, Corollary 2]. However, the successful use of their result in practice requires knowing something about the mixing proportion α1. Namely, if α1∈(δ1,δ 2) then the linear discriminant θ/θis found as the minimizer of κ, whereas if α1∈(0,δ 1)∪(δ2,1) then θ/θis found as the maximizer of κ. Naturally, as a workaround, one could in practice always search for both the minimizer and the maximizer of κbut, even in this case, it might be non-trivial to recognize the linear discriminant amongst the two. Thus, to obtain a truly unsupervised estimator, we propose instead using the squared excess kurtosis {κ(u)−3}2as an objective function. Indeed, the next lemma reveals that the squared excess kurtosis yields a Fisher consistent estimate of the linear discriminant, apart from the degenerate cases α1∈{δ1,δ 2}where excess kurtosis vanishes, without the need to choose between minimization and maximization. Lemma 1. Given model (1), 1) if α1/∈{δ1,δ 2}, then the function u→{κ(u)−3}2is uniquely maximized by ±θ/θ, 2) if α1∈{δ1,δ 2},then{κ(u)−3}2=0for all u∈Sp−1. Moving next to study the asymptotic properties of κ,letx1,...,xnbe a random sample from the marginal distribution of xin the model (1). The sample counterpart of κis κn:Sp−1→R,κ n(u)= (1/n)n i=1(u˜ xi)4 {(1/n)n i=1(u˜ xi)2}2, 6684 U. Radojiˇci´cetal. where ˜ xi:= xi−¯ x.Ifn≥pthe denominator of the random function κnis a.s. positive, making κnwell-defined and an estimator for θ/θis then obtained as any maximizer of u→{κn(u)−3}2(note that if n≥p, a maximizer exists almost surely due to the compacity of Sp−1). The following theorem shows that any sequence of such maximizers has a limiting normal distribution. Note that the need to include the “corrective” signs snin Theorem 2stems from the sign-invariance of the objective function (which also causes the existence of two maximizers in Lemma 1). Theorem 2. Given model (1), assume that α1/∈{δ1,δ 2}and let unbe any sequence of maximizers of u→{κn(u)−3}2. Then, there exists a sequence of signs sn∈{−1,1}such that, as n→∞, 1) snun→θ/θ, almost surely. 2) √n(snun−θ/θ)Np(0,Ψκ),where Ψκ:= CκΨU, and Cκ:= 6+24βτ +9β(1 −2β)τ2+β(1 −3β)τ3 βτ3(1 −6β)2, with ΨU,βand τas in Theorem 1. The limiting covariance matrices in Theorems 1and 2are proportional, the only difference being the factor Cκ. This makes their comparisons in Subsection 2.3 particularly straightforward. However, even without the formal comparisons, it is evident that the kurtosis-based estimator has a clear flaw in that it fails to be consistent for the mixing proportions α1∈{δ1,δ 2}(for these values of α1, we have 1 −6β= 0 in the denominator of Cκin Theorem 2). And even though these are only two points in the continuum (0,1), the continuity of Cκin α1 outside of these points implies that the estimator is highly inefficient for values of α1near δ1or δ2. Hence, we will next discuss an alternative estimator that is consistent when α1∈{δ1,δ 2}(at the price of lacking consistency in another point). Skewness-based projection pursuit To complement the kurtosis-based projection pursuit, we next consider skewnessbased projection pursuit. Note that, despite its dependency on lower moments, this form of PP is less studied in the literature (see the references in Section 1). The skewness of the projection of xon a given direction u∈Sp−1is measured by the objective function γ:Sp−1→Rdefined as, γ(u)= E{(u˜ x)3} [E{(u˜ x)2}]3/2, Properties of PP-based estimators of the linear discriminant 6691 Fig 6. Relative asymptotic efficiencies of the hybrid estimator as a function of w1when τ=5 and α1=δ1+ε,whereε∈{0.001,0.002,0.005,0.010}. The value of the weight w1achieving the maximal efficiency can be seen to approach zero when ε→0. still in more detail in Figure 6which plots the relative asymptotic efficiency C−1 η as a function of w1for τ=5andα1−δ1=: ε∈{0.001,0.002,0.005,0.010}. The weight achieving the maximal efficiency indeed approaches zero as ε→0. Algebraically, it is easy to see what is happening: For τ→∞and β≈1/6, the approximation of Cη, obtained by ignoring the terms of order (1 −6β)2and higher is Cη≈1 1−4β+8 3 1−w1 w1 1−6β β(1−4β).Forα1→δ1from the inside of the interval we have 1 −6β<0, which yields that Cηis minimized for w1→0. Similarly, for α1→δ1from the outside of the interval we have 1 −6β>0, which yields that Cη≥0 and shows that it is minimized for w1= 1. Also, as β→1/6, no matter from which side, Cηconverges to constant in w1, implying that there is no discontinuity in the efficiency value itself. No such behavior is observed for α1→0.5 and the reason for this is that 1 −4β≥0, implying that no sign change occurs when passing the critical value α1=1/2. The efficiencies achieved by the optimal weighting are shown by the solid black line in Figure 5and indicate that the hybrid estimator is able to reach satisfying levels of efficiency, particularly when the mixing proportion lies in the interval (δ1,δ 2). Indeed, in the limit τ→∞, within the interval there always exists a weighting that reaches efficiency equal to LDA, as evidenced by the right-most panel of Figure 5. On the other hand, outside of the interval (δ1,δ 2), LDA is still, even in the limit τ→∞, a superior choice. We conjecture that the reason for this critical difference in behavior inside and outside of the interval is that when the value of α1is extreme, one of the groups is small, making the pin-pointing of the optimal direction difficult in general, but even more so for the unsupervised methods which have no class information available. However, it is not clear why the particular points δ1,δ 2serve as the cut-off values for this 6692 U. Radojiˇci´cetal. behavior. Finally, note that the discontinuities make the use of the optimal choice of weighting somewhat difficult in practice, as, if one’s prior information/guess on the value of the mixing proportion α1is even slightly off, relying on the seemingly optimal choice can in the worst case lead to relative efficiency close to zero. Moreover, recall that the previous experiments were asymptotical in nature and do not necessarily reflect the behavior of the method under sample sizes encountered in practical situations. Hence, we suggest using a “safe” universal value of w1, most preferably falling in the interval (0.725,0.825) identified in conjunction with Figure 4. For example, Figure 3shows that the value w1=0.80 delivers, for finite τ, performance not far behind the optimal choice for any α1. However, if one is reasonably certain about the value of α1(which, optimally, is far away from δ1,δ 2) and has nsufficiently large, resorting to the optimal choice is, of course, also possible. 4. High-dimensional projection pursuit In this section, we study projection pursuit in a high-dimensional regime where both the sample size nand the dimension p≡pndiverge to infinity in a suitable rate. That is, we will work with the n-indexed sequence of high-dimensional centered normal mixtures, xn∼α1Npn(−α2hn,Σn)+α2Npn(α1hn,Σn),(2) where the parameters pn∈N,hn∈Rpnand Σn∈Rpn×pnare all functions of the sample size n,andwhereα1,α 2∈(0,1) are taken to be fixed. The notation Sp−1refers to the unit sphere in Rpand Σn2denotes the spectral norm of the matrix Σn. The specific forms for the two locations guarantee that xnhas zero mean and is without loss of generality. For simplicity, we work with the assumption that the location of the data is known (and equals zero), allowing us to consider non-centered quantities in our objective functions. For a fixed n, the population and sample skewness-based objective functions are thus γ2 n0:Spn−1→Rand γ2 n:Spn−1→R,with, γn0(u)= E{(uxn)3} [E{(uxn)2}]3/2and γn(u)= 1 nn i=1(uxni)3 {1 nn i=1(uxni)2}3/2, while the population and sample kurtosis-based objective functions are (κn0− 3)2:Spn−1→Rand (κn−3)2:Spn−1→R,with, κn0(u)= E{(uxn)4} [E{(uxn)2}]2and κn(u)= 1 nn i=1(uxni)4 {1 nn i=1(uxni)2}2. Furthermore, denote for a fixed n, the population and sample hybrid objective functions ηn0:Spn−1→Rand ηn:Spn−1→R,with, ηn0(u;w1)=w1γ2 n0(u)+(1−w1){κn0(u)−3}2and Properties of PP-based estimators of the linear discriminant 6693 ηn(u;w1)=w1γ2 n(u)+(1−w1){κn(u)−3}2, respectively. Note that, indeed, also the population objective functions are now indexed by nsince our model evolves with the growing sample size. By Theorems 3,4and Lemma 3, the unique maximizers of u→ γ2 n0(u), u→{κn0(u)− 3}2,andu→ ηn0(u;w1), w1∈(0,1), for α1=α2,α1=1/2±1/√12 and α1∈(0,1), respectively, are now ±θnwhere θn:= Σ−1 nhn/Σ−1 nhn. As our main results in this section, we show that projection pursuits using all three considered projection indices produce consistent estimates of θnin the high-dimensional regime (2) where the dimension pn→∞, as long as its growth rate is sufficiently slow compared to nand the model parameters are bounded in size from both above and below. Note that, since the target parameter θn is indexed by n, by its “consistent estimate” we mean that the angle between our estimator unand θngets arbitrarily small in the sense of convergence of probability, |u nθn|−1→p0asn→∞. Our technique of proof is essentially a high-dimensional version of the standard M-estimator argument where additional care has been taken to accommodate the n-indexed model (2). Theorem 5. Let xn1,...,xnn be a random sample (a triangular array) from the model (2)where α1=α2and assume that there exists C1,C 2>0such that, for all n, 1/C1≤hn≤C1,Σn2≤C2,Σ−1 n2≤C2. Assume further that, pn→∞ and pn=o(n1/3). Then any sequence unof maximizers of u→ γ2 n(u)satisfies, |u nθn|−1→p0, as n→∞. The equivalent statement to Theorem 5holds for kurtosis-based projection pursuit as well and, since the proof of this is exactly analogous to that of Theorem 5, we have refrained from including it in Appendix B. Theorem 6. Let xn1,...,xnn be a random sample (a triangular array) from the model (2)where α1=α2and assume that there exists C1,C 2>0such that, for all n, 1/C1≤hn≤C1,Σn2≤C2,Σ−1 n2≤C2. Assume further that, pn→∞ and pn=o(n1/4). Then any sequence unof maximizers of u→{κn(u)−3}2satisfies, |u nθn|−1→p0, as n→∞. 6694 U. Radojiˇci´cetal. Theorems 5and 6imply the equivalent statement for the hybrid estimator as well, whose proof we again omit. Theorem 7. Let xn1,...,xnn be a random sample (a triangular array) from the model (2)where α1=α2and assume that there exists C1,C 2>0such that, for all n, 1/C1≤hn≤C1,Σn2≤C2,Σ−1 n2≤C2. Assume further that, pn→∞ and pn=o(n1/4). Then any sequence unof maximizers of u→ ηn(u;w1),0<w 1<1, satisfies, |u nθn|−1→p0, as n→∞. We note that while some of the assumptions of Theorems 5-7may seem counter-intuitive given the nature of the problem (e.g., one might expect that the problem would get easier when the distance hnbetween the groups increases), they are in fact needed for controlling the third and fourth moments of the distribution. Similarly, the growth rates pn=o(n1/3)andpn=o(n1/4) are indeed a consequence of using third and fourth moment based objective functions, respectively. Finally, a natural continuation to the above consistency results would be to derive limiting distributions for the corresponding quantities |u nθn|−1, in order to allow efficiency comparisons also in the high-dimensional case. However, this task is beyond our current scope and thus left for future work. 5. Simulations The three projection pursuit estimators considered here have been discussed in the context of ICA in detail in [63] where also fixed point algorithms for their computations are described. For our purpose here we can use their deflationbased algorithms when only one direction is to be extracted. Projection pursuit is considered notoriously prone to local optima and therefore it is known that good initial values for such algorithms are crucial. [63] suggest to use initial values based on a simple ICA method called FOBI [13]. This is also suitable in our context as the normal mixture is a sub-model of the IC model, see Lemma B.1 in Appendix B. The algorithms of [63] are implemented in the package ICtest [49], which we will use in the following together with R 3.6.1 [59]. Further details about the software used are contained in the appendix. Let unbe any of the estimators of θ/θdiscussed in Section 2. The accuracy of the estimator can in simulations be measured through the inner product u nθ/θ, which, by the Cauchy-Schwarz inequality, achieves the absolute value one if and only if the two vectors are parallel. In the continuation, we call the presented inner product the “Maximal similarity index” (MSI). The following lemma presents the limiting distribution of this performance measure. Properties of PP-based estimators of the linear discriminant 6695 Lemma 4. Let the unit length vector unsatisfy √n(snun−θ/θ)Np(0,Ψ) for some sequence of signs sn∈{−1,1}and limiting covariance matrix Ψ.Then, as n→∞, 2n(1 −snu nθ/θ)zΨz,(3) where the random vector zobeys the p-variate standard normal distribution. Moreover, the expected value of the right-hand side of (3)is tr(Ψ). Note that the sign correction in Lemma 4can be incorporated in practice by choosing the sign snsuch that the quantity snu nθ/θis positive. In the simulations, we will evaluate the performances of methods through the left-hand side of (3). By Lemma 4the average of this criterion over several replicates should be close to the trace of the limiting covariance matrix of the corresponding estimator, for sample size nlarge enough. Hence, the simulations also serve to “verify” our asymptotic results. In the following simulations four projection pursuit (PP) directions have been calculated: kurtosis based (obtained by maximization of (κn−3)2), skewness based (obtained by maximization of γ2 n), “safe” hybrid estimator (obtained by maximization of ηnfor w1=0.8) and “optimal” hybrid estimator (obtained by maximization of ηnfor w1=w1(α1,τ) which maximizes the relative asymptotic efficiency of the hybrid estimator w.r.t. LDA). Sign snis chosen such that snu nθ/||θ|| ≥ 0. As discussed, the performances of the four presented PP directions unare evaluated using the maximal similarity index u nθ/||θ|| from above. For the first simulation setting, the means of maximal similarity indices snu nθ/||θ|| in the simulations are obtained using 1000 random samples in each setting. In each (τ,α1,n)-setting, where the Mahalanobis distance between the group means τ=1,2,...,20, mixing proportion α1=0.05,0.1,...,0.45,0.5and sample size n= 500,1000,2000,4000,8000, 16000,32000, m= 1000 random samples are generated from a 10-dimensional normal mixture α1N10(0,Σ)+ (1 −α1)N10(μ2,Σ), where Σis a covariance matrix with autoregression AR(1) structure with ρ=0.6, and μis in each setting chosen randomly such that μΣ−1μ=τ. The heatmaps of the MSI-values in Figure 7show that for moderate sample sizes (n≥2000), both hybrid estimators estimate the optimal LDA direction very well. It is also visible that kurtosis and skewness based PP directions perform very badly when α1is near their corresponding discontinuity points. The hybrid estimators suffer from the same problem when α1is near 1/2−1/√12. The hybrid estimator with optimal w1is performing worse when α1is approaching 1/2−1/√12 from the inside the interval (1/2−1/√12,1/2+1/√12). It is important to recall, that the criterion for choosing the optimal weight w1is an asymptotic one and thus might not perform well in small sample settings. Furthermore, to calculate the optimal weight w1for the hybrid estimator, one needs to know both the Mahalanobis distance τbetween the group means and the mixing proportion α1, which is a rather unrealistic requirement in practice. Luckily, the hybrid estimator with the “safe” weight w1=0.8showsavery 6696 U. Radojiˇci´cetal. Fig 7. Average values of the MSI snu nθ/||θ|| as a function of Mahalanobis distance τbetween the group means and mixing proportion α1,whereunis one of the four estimators discussed above. good performance in this simulation study and is therefore recommended in cases where knowledge of α1and τis lacking. Another observation based on this simulation is that for small sample sizes skewness based PP seems to be preferable. This might be due to the fact that moments of order three are easier to estimate than moments of order four. The corresponding heatmap of the standard deviation of MSI can be found in the Appendix, Figure C.1, and shows that for sample size and distance between the group means moderately large, deviation of the MSI to the corresponding mean, which is very close to the optimal value of 1, is negligible, for most of the values of the mixing proportion α1. Heatmaps of mean and standard deviation of the MSI in Figures C.3 and C.2 show that for large sample sizes (n= 8000,16000,32000) MSI is virtually 1. In the next simulation, the theoretic results are to be confirmed by exploiting the results of Lemma 4. For that purpose we select three values for τ∈{1,5,10}, to represent hardly, moderately, and clearly separated clusters, respectively. Then we simulate, for sample sizes n= 500,1000,2000,4000,8000, 16000,32000, from a three-variate Gaussian mixture model as specified above with Σ=I3+13, where 13is matrix of ones, μ1=0and μ2=μ2(τ) is for each τchosen such that Mahalanobis distance between the means is equal to τ, i.e., μ2(1) = (0.68,−0.55,0.6),μ2(5) = (0.81,−2.24,−0.36),μ2(10) = (3.06,1.6,−1.11). Properties of PP-based estimators of the linear discriminant 6697 Fig 8. Average values of the MSI snu nθ/||θ|| as a function of n, for mixing proportion α1∈{0.1,0.2,0.3,0.4,0.5}and Mahalanobis distance between the group means τ∈{1,5,10}, where unis one of the four estimators discussed above. For α1=0.1,0.2,...,0.5 we compute then for sample sizes n= 500,1000,2000, 4000,8000,16000,32000 the means of the 2n(1−snu n)θ/θbased on 2000 repetitions. The question is then whether those averages for the presented methods stabilize for the cases they are expected to work. In order not to clutter the figure we computed only for α1=0.1 trace tr(Ψ), for the corresponding matrix Ψ.Figure8shows the results of this simulation and confirms the corresponding theoretic findings from above. The less separated the clusters are the more difficult the estimation and even for n= 32000 observations there is no stabilization visible. But the more separated the two clusters are, the faster the stabilization. Also, the closer α1is to the critical value 1/2−1/√12, the worse is kurtosis based PP. Skewness based PP similarly is better the more skewed the distribution and clearly does not work in the symmetric case. Both hybrid estimators show excellent performance in this setting. It is also clearly visible that the empirical lines for the mixing proportion α1=0.1 correspond to the theoretically computed dashed lines given that the groups are separated enough and we assume for τ= 1 the line would be reached for much larger sample sizes. Though we show the theoretic lines only for one mixing proportion the behavior is similar for all others naturally with the exception of skewness not working in the symmetric case. 6698 U. Radojiˇci´cetal. Fig 9. PP direction unbased on LDA (here denoted as “LDA, sample”), PCA, and one of the four estimators discussed above. Principal component analysis (PCA) can also be seen as a projection method where the variance is maximized. PCA is arguably the most popular dimension reduction method and is often used before clustering. While skewness and kurtosis can be related to mixtures the variance does not have the same connection with the discriminant direction as the other cumulants. A theoretic consideration of when PCA can be used to estimate the discriminant is in Appendix A. Here we show an example where PCA fails. Figure 9visualizes a sample of size n= 100 from our Gaussian mixture model with μ1=(0,0),μ2= (0,5),Σ=10 0.3 0.31  and mixing proportion α1=0.3. The figure contains then the direction θ/θof the population LDA as well as its estimate together with our four PP methods considered in this section and the direction of the first principal component. As it is clearly visible here, there is no big difference between the methods except for PCA which points in a direction that contains no information for the separation of the two groups. 6. Real data example To evaluate the performance of the hybrid estimator in a real data set we consider the finance data set available in the R package Rmixmod [37], which consists of 889 records of companies where based on four numeric summary statistics, it should be decided if the company is financially healthy or not, where the information is provided in the data set. Properties of PP-based estimators of the linear discriminant 6699 The scatter plot matrix is given in the Appendix as Figure C.4 and shows no clear clusters. As a reference we compute for the data set LDA and then compare this supervised estimate via the estimate snu nθn/θnof the MSI to our hybrid estimator for different weights, to PCA and to an estimate obtained by fitting a two-component Gaussian mixture model via the expectation - maximization (EM) algorithm. Namely, since the nature of the presented problem is essentially clustering, it is only natural to consider also the classification by EM-algorithm [18]. We consider the EM-algorithm implemented in R package mclust [57], where it is being initialized using the initial partitions from model-based hierarchical agglomerative clustering. By assuming the model (1), the EM-algorithm estimates the parameters of the Gaussian mixture; covariance matrix ΣEM and group means μ1,EM and μ2,EM . Then, the estimated mean and the covariance matrix can be used in order to estimate the linear discriminant direction uEM =Σ−1 EM(μ2,EM −μ1,EM )/ΣEM(μ2,EM −μ1,EM ). We further refer to the estimator uEM of the LDA direction as the mclust-based estimate. Figure 10 shows obtained MSI values for the discussed estimators. The figure clearly shows that as long as enough weight is given to kurtosis, the hybrid estimator based PP clearly outperforms both PCA and mclust based estimators, while its performance is poor if skewness gets too much weight. This is not surprising as the amount of healthy (457) and bankrupt (432) companies is almost equal. The weight of 0.8 gives again a good performance. Nevertheless, even though for most values of w1, and especially for suggested w1∈[0.7,0.8], hybrid estimators clearly outperform both PCA and mclust based estimators, the achieved MSI values of around 0.5 are not ideal. Such performance can be explained by the low sample size and that the cluster centers are not that far apart, as is shown in the boxplots of Figure C.5 in the Appendix C,which also indicates that the results obtained by mclust based estimation are not satisfactory. 7. Discussion In this paper, we conducted an asymptotic comparison of two popular estimators of the linear discriminant direction, LDA, and projection pursuit based on skewness and kurtosis. For the latter, we proposed using the weighted combination of kurtosis and skewness as the projection index (giving the individual cumulants as special cases). Both the theoretical results and simulations indicate that, with a suitable choice of weighting, such projection pursuit achieves reasonably good performance compared to LDA (e.g., around 15% relative asymptotic efficiency if the Mahalanobis distance between the groups is 5, see Figure 3), considering it operates in the complete absence of any label information. Moreover, in the extreme case of balanced and infinitely well-separated groups, projection pursuit is able to reach asymptotic efficiency equal to LDA with an optimal choice of weighting. 6700 U. Radojiˇci´cetal. Fig 10.MSIsnunθn/θnfor the finance data, where unis the direction based on PCA and PP estimators for w1=0,0.1,0.2,...,0.9,1.θnis the direction based on LDA. The use of our optimal weighting results is difficult in practice by the discontinuities around the mixing proportions δ1,δ 2observed in Section 3, see Figure 5. As such, unless one is sure that the mixing proportion is not in these regions, our recommendation is to use a universal choice of weighting, anything between 0.7and0.8 (as the weight for skewness) likely being a good choice. At first, we thought that the discontinuities, and the surprising recommendation to favor kurtosis just outside the interval (δ1,δ 2), might be caused by the uneven robustness properties of skewness and kurtosis in the objective function. Namely, being based on fourth moments, kurtosis is more affected by outliers than skewness (despite the standardization with second moments). Hence, we also considered using the “balanced” objective function, η∗(u):=w1γ(u)8/3+w2{κ(u)−3}2, in an attempt to put skewness and kurtosis on an equal footing. However, the asymptotic properties of η∗(not shown here) turn out to be essentially the same as for η, including the discontinuities which are also observed for it. Note also that the discontinuous behavior was observed also in [63], where the normal mixture model was studied using independent component analysis. Similarly one could extend our considerations here to many other PP indices as well, which often are modifications of skewness or kurtosis (see e.g. [25]) or otherwise motivated to be useful in clustering or structure detection, see, for example, [14,17] and references therein for alternative indices. These indices are however often computationally expensive and therefore much less popular than skewness and kurtosis. Finally, besides projection pursuit, there exist also other unsupervised estimators of the linear discriminant. For example, it is known that the linear discriminant can be reconstructed using invariant coordinate selection (ICS) Properties of PP-based estimators of the linear discriminant 6707 ProofofLemmaB.3.The conditional distribution of zgiven vzis z|vz=s∼N(v−2sv,Ip−v−2vv). Thus, E{(vz)kz}= E[E{(vz)kz|vz}] = E[(vz)kE{z|vz}]=v−2E{(vz)k+1}v. The second claim is shown analogously and by using the fact that E(zz|vz)= Cov(z|vz)+E(z|vz)E(z|vz). ProofofLemma1.The distribution of the projection u˜ xis u˜ x∼α1N(−α2t, g)+α2N(α1t, g), where t:= uh,g:= uΣu and h:= μ2−μ1. By the moment formulas of univariate normal distribution, E{(u˜ x)2}=g+βt2,whereβ:= α1α2. Similarly, E{(u˜ x)4}=β(α3 1+α3 2)t4+6βt2g+3g2which can be further simplified by noting that α3 1+α3 2=1−3β. Hence, {κ(u)−3}2=β2(1 −6β)2f4 (1 + βf)4,(B.2) where f:= t2/g ≥0. If α1∈{δ1,δ 2},then1−6β=0making{κ(u)−3}2= 0. Assume then that α1/∈{δ1,δ 2}, implying that (1 −6β)2>0. The derivative of the map x→ x4/(1 + βx)4is 4x3/(1 + βx)5, showing that the map is strictly increasing in (0,∞). Hence, {κ(u)−3}2is maximal when fis at its largest. Now, f=t2 g= Σ1/2u Σ1/2u Σ−1/2h2 , showing that, by the Cauchy-Schwarz inequality, fis maximal if and only if Σ1/2u∝Σ−1/2h, i.e., when u=±θ/θ(where θ=Σ−1h). ProofofTheorem2.The objective functions are translation invariant, meaning that we may, without loss of generality, assume that E(x)=0. This makes the marginal distribution of xbe x∼α1Np(−α2h,Σ)+α2Np(α1h,Σ), where h:= μ2−μ1. The strong consistency of the estimator can be shown in the usual way by establishing that the objective function is strongly uniformly convergent in the compact parameter set Sp−1(or, more precisely, in its subset where the sign of the estimator is fixed), that is, sup u∈Sp−1|{κn(u)−3}2−{κ(u)−3}2|→0,a.s.(B.3) For simplicity, we give the proof of the uniform convergence only in Theorem 3, in the context of skewness (having lower moments than kurtosis), and similar (but lengthier) arguments can be used to show (B.3). 6708 U. Radojiˇci´cetal. To show the limiting normality, note that the Largrangian corresponding to the optimization problem is n(u)={κn(u)−3}2+λn(uu−1) where λnis the Lagrangian multiplier. Using some matrix calculus, the corresponding gradient is seen to be ∇n(u)= 8 ˜sn2(u)3{κn(u)−3}{˜sn2(u)˜ mn3(u)−˜sn4(u)˜ mn1(u)}−2λnu, where ˜snk(u):=(1/n)i(u˜ xi)kand ˜ mnk(u):=(1/n)i(u˜ xi)k˜ xi. The gradient vanishes at the (sign-adjusted) sample maximum snunand multiplication of the gradient from the left with snu nthus yields that 0 = snu n∇n(snun)= −2λn, showing that λn=0. We next work on the level of individual probability elements ω∈Ω. By Lemma 1, LLN and the strong consistency of snun, there exists a probability one set Hsuch that snun→u0and κn(un)−3→t=0forallω∈H. Thus, for each ω∈H, the maximizer unsatisfies, for nlarge enough, the estimating equation ˜sn2(snun)˜ mn3(snun)−˜sn4(snun)˜ mn1(snun)=0. Using Lagrangian multipliers we can similarly show that the population maximizer u0:= θ/θsatisfies s2(u0)m3(u0)−s4(u0)m1(u0) = 0, where sk(u)=E{(ux)k}and mk(u)= E{(ux)kx}. Let gnκ :Rp\{0}→Rbe such that gnκ(u)=˜sn2(u)˜ mn3(u)−˜sn4(u)˜ mn1(u). For each ω∈H,wehave,fornlarge enough, the Taylor expansion gnκ(snun)=gnκ(u0)+∇gnκ(u0)(snun−u0) +{(snun−u0)×∇∇gnκ(˜ un)}(snun−u0), where ∇∇gnκ(˜ un) is the third order tensor of second derivatives of g, the symbol ×denotes the vector-by-tensor multiplication (producing a matrix) and ˜ un satisfies ˜ un−u0≤un−u0, implying that ˜ un→u0. Multiplying the expansion by √nand using the fact that gnκ(snun)=0gives that {(snun−u0)×∇∇gnκ(˜ un)+∇gnκ(u0)}√n(snun−u0)=√ngnκ(u0). (B.4) Now, the elements of ∇∇gnκ(u) are polynomials of the sample moments of ˜ xiand the elements of uimplying that, by LLN, ∇∇gnκ(˜ un) converges to a constant and (snun−u0)×∇∇gnκ(˜ un) converges to zero, for any ω∈H.Now, by the unit lengths of snunand u0,wehavec0h(snun+u0)√n(snun−u0)= 0, where c0:= (1/2){3s2(u0)2−s4(u0)}θ(hΣ−1h)−1(the inclusion of the constant c0simplifies things later on). Summing this with equation (B.4)gives {(snun−u0)×∇∇gnκ(˜ un)+∇gnκ(u0)+c0h(snun+u0)}√n(snun−u0) =√ngnκ(u0). Assume now for a moment that ∇gnκ(u0)+c0h(snun+u0)converges to a full-rank matrix G∈Rp×p. Then, for nlarge enough, we have, √n(snun−u0)={(snun−u0)×∇∇gnκ(˜ un)+∇gnκ(u0)(B.5) Properties of PP-based estimators of the linear discriminant 6709 +c0h(snun+u0)}−1√ngnκ(u0).(B.6) Hence, assuming further that we have √ngnκ(u0)Np(0,Π), then the limiting distribution of snunis, by Slutsky’s theorem, √n(snun−u0)Np{0,G−1Π(G−1)}.(B.7) Thus, to complete the proof, we next derive expressions for Gand Π(and show that the former has indeed full rank). The Jacobian of gnκ is, ∇gnκ(u0)=2˜ mn1(u0)˜ mn3(u0)+3˜sn2(u0)˜ Gn2(u0)−4˜ mn3(u0)˜ mn1(u0) −˜sn4(u0)˜ Gn0(u0), where ˜ Gnk(u):=(1/n)i(u˜ xi)k˜ xi˜ x i. Thus, by LLN and using the population level estimating equation, s2(u0)m3(u0)=s4(u0)m1(u0), we get ∇gnκ(u0)→p=−2s4(u0) s2(u0)m1(u0)m1(u0)+3s2(u0)G2(u0)−s4(u0)G0(u0), (B.8) where Gk(u):=E{(ux)kxx}. Denote next τ:= hΣ−1h. To compute the moments mk(u0)andGk(u0), we use Lemma B.3.Theformer satisfies Σ−1/2mk(u0)=θ−kE{(vz)kz},wherev:= Σ−1/2hand z∼ α1Np(−α2h,Ip)+α2Np(α1h,Ip). Denoting the components of the mixture by z1and z2, we have, by the first part of Lemma B.3,forz1that E{(vz1)kz1}= v−2E{(vz1)k+1}v, and similarly for z2. Hence, Σ−1/2mk(u0)=θ−kv−2E{(vz)k+1}v. Finally, since sk(u0)=θ−kE{(vz)k}, we get mk(u0)=θv−2sk+1(u0)Σ1/2v=θτ−1sk+1(u0)h.(B.9) For Gk(u0), we have, using the same notation, that Σ−1/2Gk(u0)Σ−1/2=θ−kE{(vz)kzz}. The second part of Lemma B.3 then shows that Gk(u0)=θ−kΣ1/2[E{(vz)k}(Ip−v−2vv)+v−4E{(vz)k+2}vv]Σ1/2 =θ−k[θksk(u0)(Σ−τ−1hh)+τ−2θk+2sk+2(u0)hh] =sk(u0)Σ+τ−1{τ−1θ2sk+2(u0)−sk(u0)}hh. (B.10) Plugging in the expressions to (B.8), we get ∇gnκ(u0)→p(3s2 2−s4)Σ1/2(Ip− ww)Σ1/2,wherew:= Σ−1/2h/Σ−1/2hand sk≡sk(u0). Moreover, we also 6710 U. Radojiˇci´cetal. have c0h(snun+u0)→p(3s2 2−s4)Σ1/2wwΣ−1Σ1/2.NowGis the sum of these two, giving, G= (3s2 2−s4)Σ1/2(Ip−ww+wwΣ−1)Σ1/2. The invertibility of Gnow follows from Lemma B.2, which also gives (Ip−ww+wwΣ−1)−1=Ip+(wΣ−1w)−1ww(Ip−Σ−1) =Ip+θ−2Σ−1/2hhΣ−1/2(Ip−Σ−1). Finally, this makes the inverse of Gbe, G−1=1 3s2 2−s4Σ−1+1 θ2θθ(Ip−Σ−1). The fact that 3s2 2−s4= 0 follows from the formulas for skgiven later in the proof. We next obtain the limiting distribution of √ngnκ(u0)=√n{˜sn2(u0)˜ mn3(u0)−˜sn4(u0)˜ mn1(u0)}. Define non-centered counterparts for the sample moments as snk(u):=(1/n) i(uxi)kand mnk(u):=(1/n)i(uxi)kxi. Then, LLN together with the calculus of op(1) and Op(1) sequences shows that ˜sn2(u)=sn2(u)+op(1/√n) and ˜ mn1(u)=mn1(u)+op(1/√n). However, the same equivalence does not hold for the terms ˜ mn3(u)and˜sn4(u) but we instead have ˜ mn3(u)=mn3(u)−3sn1(u)m2(u)−s3(u)mn0(u)+op(1/√n), and ˜sn4(u)=sn4(u)−4s3(u)sn1(u)+op(1/√n). Using these, we expand √ngnκ(u0) to be (dropping u0from the notation), √ngnκ =√n(sn2−s2)m3+s2√n(mn3−m3)−√n(sn4−s4)m1 −s4√n(mn1−m1)+(4s3m1−3s2m2)√nsn1−s2s3√nmn0. (B.11) Hence, by CLT, √ngnκ has a limiting normal distribution with the covariance matrix, Π=Cov{(u 0x)2m3+s2(u 0x)3x−(u 0x)4m1−s4(u 0x)x + (4s3m1−3s2m2)(u 0x)−s2s3x}. This matrix consists of 36 terms, which we next present and simplify using (B.9) and (B.10). We use the notation ψ=θτ−1. Note that s1=0,m0=0and f:= 4s3m1−3s2m2=ψs2s3h. Properties of PP-based estimators of the linear discriminant 6711 (1, 1): (s4−s2 2)m3m 3=ψ2s2 4(s4−s2 2)hh. (1, 2): s2(m3m 5−s2m3m 3)=ψ2s2s4(s6−s2s4)hh. (1, 3): −(s6−s2s4)m3m 1=−ψ2s2s4(s6−s2s4)hh. (1, 4): −s4(m3m 3−s2m3m 1)=−ψ2s2 4(s4−s2 2)hh. (1, 5): s3m3f=ψ2s2s2 3s4hh. (1, 6): −s2s3m3m 2=−ψ2s2s2 3s4hh. (2, 1): ψ2s2s4(s6−s2s4)hh. (2, 2): s2 2(G6−m3m 3)=s2 2s6Σ+s2 2{ψ2(s8−s2 4)−τ−1s6}hh. (2, 3): −s2(m7m 1−s4m3m 1)=−ψ2s2 2(s8−s2 4)hh. (2, 4): −s2s4(G4−m3m 1)=−s2s2 4Σ−s2s4{ψ2(s6−s2s4)−τ−1s4}hh. (2, 5): s2m4f=ψ2s2 2s3s5hh. (2, 6): −s2 2s3G3=−s2 2s2 3Σ−s2 2s3{ψ2s5−τ−1s3}hh. (3, 1): −ψ2s2s4(s6−s2s4)hh. (3, 2): −ψ2s2 2(s8−s2 4)hh. (3, 3): (s8−s2 4)m1m 1=ψ2s2 2(s8−s2 4)hh. (3, 4): s4(m1m 5−s4m1m 1)=ψ2s2s4(s6−s2s4)hh. (3, 5): −s5m1f=−ψ2s2 2s3s5hh. (3, 6): s2s3m1m 4=ψ2s2 2s3s5hh. (4, 1): −ψ2s2 4(s4−s2 2)hh. (4, 2): −s2s2 4Σ−s2s4{ψ2(s6−s4s2)−τ−1s4}hh. (4, 3): ψ2s2s4(s6−s2s4)hh. (4, 4): s2 4(G2−m1m 1)=s2s2 4Σ+s2 4{ψ2(s4−s2 2)−τ−1s2}hh. (4, 5): −s4m2f=−ψ2s2s2 3s4hh. (4, 6): s2s3s4G1=ψ2s2s2 3s4hh. (5, 1): ψ2s2s2 3s4hh. (5, 2): ψ2s2 2s3s5hh. (5, 3): −ψ2s2 2s3s5hh. (5, 4): −ψ2s2s2 3s4hh. (5, 5): s2ff=ψ2s3 2s2 3hh. (5, 6): −s2s3fm 1=−ψ2s3 2s2 3hh. (6, 1): −ψ2s2s2 3s4hh. (6, 2): −s2 2s2 3Σ−s2 2s3{ψ2s5−τ−1s3}hh. (6, 3): ψ2s2 2s3s5hh. (6, 4): ψ2s2s2 3s4hh. (6, 5): −ψ2s3 2s2 3hh. (6, 6): s2 2s2 3G0=s2 2s2 3Σ+s2 2s2 3(ψ2s2−τ−1)hh. Summation of the previous 36 terms results in Π=s2(s2s6−s2s2 3−s2 4)(Σ− τ−1hh). Thus, from the reasoning preceding (B.7), we have that √n(snun− u0) has a limiting normal distribution and with the covariance matrix Ψκ= G−1Π(G−1). Plugging now in the values of Gand Πand simplifying, we obtain, Ψκ=s2(s2s6−s2s2 3−s2 4) (3s2 2−s4)2Ip−θθ θ2Σ−1Ip−θθ θ2.(B.12) Now, recall that sk≡sk(u0)=E{(u 0x)k}=θ−kE{(θx)k}where θx∼ 6712 U. Radojiˇci´cetal. α1N1(−α2τ,τ)+α2N1(α1τ,τ)andτ=θh=θΣθ. Using the moment formulas for univariate normal distribution we now obtain that s2=θ−2τ(1 + βτ),s 3=θ−3(α1−α2)βτ3, s4=θ−4τ2{βτ2(1 −6β) + 3(1 + βτ)2}, and s6=θ−6τ3{β(1 −5β+5β2)τ3+15β(1 −3β)τ2+ 45βτ + 15} where β:= α1α2and we have used the identities α3 1+α3 2=1−3βand α5 1+α5 2= 1−5β+5β2. Plugging these in to (B.12) and simplifying (using (α1−α2)2= 1−4β), shows that the constant in front is (1 + βτ)(6 + 24βτ +9βτ2+βτ3−18β2τ2−3β2τ3) τ3β2(6β−1)2θ2 Proof of Lemma 2.The distribution of the projection u˜ xis u˜ x∼α1Np(−α2t, g)+α2Np(α1t, g), where t:= uh,g:= uΣu and h:= μ2−μ1. By the moment formulas of univariate normal distribution, E{(u˜ x)2}=g+βt2,whereβ:= α1α2. Similarly, E{(u˜ x)3}=(α1−α2)βt3. Hence, γ(u)=β2(1 −4β)f3 (1 + βf)3, where f:= t2/g ≥0. Now, if α1=α2=1/2, then clearly γ(u)2=0.Ifα1=1/2, the derivative of the map x→ x3/(1 + βx)3is 3x2/(1 + βx)4, showing that the map is strictly increasing outside of the origin. The conclusion now follows as in the proof of Lemma 1. Proof of Theorem 3.The strong consistency follows as soon as we show the strong uniform consistency, sup u∈Sp−1|γn(u)2−γ(u)2|→0,a.s.. (B.13) By Theorem 2 and Lemma 1 in [5], (B.13) holds if, 1) the parameter space is compact, 2) we have γn(u)2→γ(u)2, a.s., for all u∈Sp−1(this holds by LLN and the continuous mapping theorem), 3) γ2is uniformly continuous in uand, 4) γ2 nis Lipschitz continuous in the sense that |γn(u1)2−γn(u2)2|≤Knu1−u2 for all u1,u2∈Sp−1and some random variable Knconverging almost surely to a constant. We now verify condition 4) above. Using the notation of the proof of Theorem 2,wehave |γn(u1)2−γn(u2)2|= ˜s2 n3(u1) ˜s3 n2(u1)−˜s2 n3(u2) ˜s3 n2(u2) Properties of PP-based estimators of the linear discriminant 6713 ≤|˜s2 n3(u1)−˜s2 n3(u2)|˜s3 n2(u2)−˜s2 n3(u2)|˜s3 n2(u1)−˜s3 n2(u2)| ˜s3 n2(u1)˜s3 n2(u2). Now, ˜sn2(u)is,forallu∈Sp−1, lower bounded by the smallest eigenvalue of the sample covariance matrix, which by the continuity of the eigenvalues and the positive-definiteness of the covariance matrix converges almost surely to a positive constant. Moreover, we have |˜sn3(u)|≤1 n n  i=1 |u˜ xi|3≤1 n n  i=1 xi−¯ x3≤1 n n  i=1 (xi+¯ x)3, which converges, by LLN, almost surely to a constant, and similar result can be shown for |˜sn2(u)|. Finally, |˜s2 n3(u1)−˜s2 n3(u2)|≤|˜sn3(u1)−˜sn3(u2)|2 n n  i=1 xi−¯ x3, and |˜sn3(u1)−˜sn3(u2)|≤1 n n  i=1 |(u1−u2)˜ xi||(u 1˜ xi)2+u 1˜ xiu 2˜ xi+(u 2˜ xi)2| ≤u1−u23 n n  i=1 xi−¯ x3, and putting everything above together, we conclude that the Lipschitz continuity 4) holds. What remains to be verified is then condition 3), which can be shown similarly to 4) after recalling that Lipschitz continuity implies uniform continuity. Hence, the strong consistency of the estimator follows. Also, the proof of the limiting normality has exactly the same steps as in the proof of Theorem 2and we only provide the key steps and expressions, using the same notation as in the proof of Theorem 2. The gradient of γnis ∇γn(u)= 6 ˜sn2(u)5/2γn(u){˜sn2(u)˜ mn2(u)−˜sn3(u)˜ mn1(u)}, leading to the estimating equation gnγ(un)=0,forgnγ(u):=˜sn2(u)˜ mn2(u)− ˜sn3(u)˜ mn1(u). The Jacobian of gnγ at u0satisfies (after simplification via the estimating equation) ∇gnγ(u0)→p−s3(u0) s2(u0)m1(u0)m1(u0)+2s2(u0)G1(u0)−s3(u0)G0(u0). By formulas for mk(u0)andGk(u0), the limit equals −s3Σ1/2(Ip−ww)Σ1/2. Using the same trick as in the proof of Theorem 2to make the Jacobian full rank (addition of c0h(snun+u0)√n(snun−u0) = 0 for suitably chosen c0to the Taylor expansion), we obtain the corresponding matrix Gto be G=−s3Σ1/2(Ip−ww+wwΣ−1)Σ1/2, 6714 U. Radojiˇci´cetal. with the inverse, G−1=−1 s3Σ−1+1 θ2θθ(Ip−Σ−1). Moving to study the limiting distribution of √ngnγ, we note that ˜ mn2(u)=mn2(u)−2sn1(u)m1(u)−s2(u)mn0(u)+op(1/√n), and ˜sn3(u)=sn3(u)−3s2(u)sn1(u)+op(1/√n). With these, we expand √ngnγ (u0)tobe, √ngnγ =√n(sn2−s2)m2+s2√n(mn2−m2)−√n(sn3−s3)m1 −s3√n(mn1−m1)+s2m1√nsn1−s2 2√nmn0.(B.14) Hence, by CLT, √ngnγ has a limiting normal distribution with the covariance matrix, Π=Cov{(u 0x)2m2+s2(u 0x)2x−(u 0x)3m1−s3(u 0x)x+s2m1(u 0x)−s2 2x}. The covariance matrix has the following 36 terms. (1, 1): (s4−s2 2)m2m 2=ψ2s2 3(s4−s2 2)hh. (1, 2): s2(m2m 4−s2m2m 2)=ψ2s2s3(s5−s2s3)hh. (1, 3): −(s5−s2s3)m2m 1=−ψ2s2s3(s5−s2s3)hh. (1, 4): −s3(m2m 3−s2m2m 1)=−ψ2s2 3(s4−s2 2)hh. (1, 5): s2s3m2m 1=ψ2s2 2s2 3hh. (1, 6): −s2 2m2m 2=−ψ2s2 2s2 3hh. (2, 1): ψ2s2s3(s5−s2s3)hh. (2, 2): s2 2(G4−m2m 2)=s2 2s4(Σ−τ−1hh)+ψ2s2 2(s6−s2 3)hh. (2, 3): −s2(m5m 1−s3m2m 1)=−ψ2s2 2(s6−s2 3)hh. (2, 4): −s2s3(G3−m2m 1)=−s2s2 3(Σ−τ−1hh)−ψ2s2s3(s5−s2s3)hh. (2, 5): s2 2m3m 1=ψ2s3 2s4hh. (2, 6): −s3 2G2=−s4 2(Σ−τ−1hh)−ψ2s3 2s4hh. (3, 1): −ψ2s2s3(s5−s2s3)hh. (3, 2): −ψ2s2 2(s6−s2 3)hh. (3, 3): (s6−s2 3)m1m 1=ψ2s2 2(s6−s2 3)hh. (3, 4): s3(m1m 4−s3m1m 1)=ψ2s2s3(s5−s2s3)hh. (3, 5): −s2s4m1m 1=−ψ2s3 2s4hh. (3, 6): s2 2m1m 3=ψ2s3 2s4hh. (4, 1): −ψ2s2 3(s4−s2 2)hh. (4, 2): −s2s2 3(Σ−τ−1hh)−ψ2s2s3(s5−s2s3)hh. (4, 3): ψ2s2s3(s5−s2s3)hh. (4, 4): s2 3(G2−m1m 1)=s2s2 3(Σ−τ−1hh)+ψ2s2 3(s4−s2 2)hh. (4, 5): −s2s3m2m 1=−ψ2s2 2s2 3hh. Properties of PP-based estimators of the linear discriminant 6715 (4, 6): s2 2s3G1=ψ2s2 2s2 3hh. (5, 1): ψ2s2 2s2 3hh. (5, 2): ψ2s3 2s4hh. (5, 3): −ψ2s3 2s4hh. (5, 4): −ψ2s2 2s2 3hh. (5, 5): s3 2m1m 1=ψ2s5 2hh. (5, 6): −s3 2m1m 1=−ψ2s5 2hh. (6, 1): −ψ2s2 2s2 3hh. (6, 2): −s4 2(Σ−τ−1hh)−ψ2s3 2s4hh. (6, 3): ψ2s3 2s4hh. (6, 4): ψ2s2 2s2 3hh. (6, 5): −ψ2s5 2hh. (6, 6): s4 2G0=s4 2(Σ−τ−1hh)+ψ2s5 2hh. Summing the terms gives Π=s2(s2s4−s3 2−s2 3)(Σ−τ−1hh). This yields the limiting covariance, Ψκ=G−1Π(G−1)=s2(s2s4−s3 2−s2 3) s2 3Ip−θθ θ2Σ−1Ip−θθ θ2. (B.15) Finally, simplifying the constant in front shows that it equals (1 + βτ)(2 + 6βτ +βτ2) τ2β2(1 −4β)θ2. ProofofTheorem4.Again, the strong consistency follows as in Theorem 2and we omit its proof. For the limiting distribution, we give in the following the key steps of the proof (and use the same notation as in the proofs of Theorems 2 and 3). The gradient of ηnis ∇ηn(u)= 1 ˜s3 n2 [6w1γn(u)˜s1/2 n2(u)gnγ(u)+8w2{κn(u)−3}gnκ(u)], where gnγ(u)=˜sn2(u)˜ mn2(u)−˜sn3(u)˜ mn1(u)andgnκ(u)=˜sn2(u)˜ mn3(u)− ˜sn4(u)˜ mn1(u) were used in the proofs of Theorems 3and 2, respectively. Thus, unsolves the estimating equation gnη(un)=0,where gnη(u):=3w1γn(u)˜s1/2 n2(u)gnγ(u)+4w2{κn(u)−3}gnκ(u). The Jacobian of gnη at u0satisfies ∇gnη(u0)=3w1[∇{γn(u0)˜s1/2 n2(u0)}gnγ(u0)+γn(u0)˜s1/2 n2(u0)∇gnγ(u0)] +4w2[∇{κn(u0)−3}gnκ(u0)+{κn(u0)−3}∇gnκ(u0)]. 6716 U. Radojiˇci´cetal. Recalling that mk(u0)=ψsk+1hand Gk(u0)=sk(Σ−τ−1hh)+ψ2sk+2hh, where ψ=θτ−1, LLN now gives that gnγ(u0)→p0and gnκ(u0)→p0, implying that ∇gnη(u0)→p−s−2 2{3w1s2s2 3+4w2(s4−3s2 2)2}Σ1/2(Ip−ww)Σ1/2. Completing now this matrix to full rank through the unit length constraint on un(as in the proofs of Theorems 3and 2), we now obtain that, G−1=−s2 2 3w1s2s2 3+4w2(s4−3s2 2)2Σ−1+1 θ2θθ(Ip−Σ−1). We then derive the limiting distribution of √ngnη(u0)=3w1√n˜s1/2 n2(u0)gnγ(u0)+4w2√ngnκ(u0). Recalling that the population version satisfies 3w1γ(u0)s1/2 2(u0)gγ(u0)+4w2{κ(u0)−3}gκ(u0)=0, we get the expansion, √ngnη =3w1{√n(γn˜s1/2 n2−γs1/2 2)(s2m2−s3m1)+s−1 2s3√ngnγ} +4w2{√n(κn−κ)(s2m3−s4m1)+s−2 2(s4−3s2 2)√ngnκ}+op(1), where √ngnκ has the expansion given in (B.11). Now, s2m2−s3m1=0and s2m3−s4m1=0, implying that √ngnη =3w1s−1 2s3√ngnγ +4w2s−2 2(s4−3s2 2)√ngnκ +op(1), where √ngnγ has the expansion given in (B.14). Consequently, by CLT, √ngnη has a limiting normal distribution. By the proof of Theorem 2, the limiting covariance matrix of 4w2s−2 2(s4−3s2 2)√ngnκ is 16w2 2s−3 2(s4−3s2 2)2(s2s6−s2s2 3− s2 4)(Σ−τ−1hh) and, by the proof of Theorem 3, the limiting covariance matrix of 3w1s−1 2s3√ngnγ is 9w2 1s−1 2s2 3(s2s4−s3 2−s2 3)(Σ−τ−1hh). Thus, the limiting covariance matrix of √ngnη is {9w2 1s−1 2s2 3(s2s4−s3 2−s2 3)+16w2 2s−3 2(s4− 3s2 2)2(s2s6−s2s2 3−s2 4)}(Σ−τ−1hh)+24w1w2s−3 2s3(s4−3s2 2)Cov(y1,y2), where y1:= (u 0x)2m2+s2(u 0x)2x−(u 0x)3m1−s3(u 0x)x+s2(u 0x)m1−s2 2x, y2:= (u 0x)2m3+s2(u 0x)3x−(u 0x)4m1−s4(u 0x)x+s3(u 0x)m1−s2s3x. The matrix Cov(y1,y2) consists of the following 36 terms: (1, 1): ψ2s3s4(s4−s2 2)hh. (1, 2): ψ2s2s3(s6−s2s4)hh. (1, 3): −ψ2s2s3(s6−s2s4)hh. (1, 4): −ψ2s3s4(s4−s2 2)hh. (1, 5): ψ2s2s3 3hh. Properties of PP-based estimators of the linear discriminant 6723 ProofofTheoremB.1.Let Rn:= (hjk)=(1/n) n  i=1{xnijxnikxni −E(xnjxnkxn)} be the pn×pn×pnthird-order symmetric tensor containing all centered sample third moments. Given u1,u2,u3∈Rpnwe denote by Rn×(u1⊗u2⊗u3) the scalar pn j=1 pn k=1 pn =1 u1ju2ku3hjk. Using this notation, our quantity of interest is, Rn2:= sup u∈Spn−1|Rn×(u⊗u⊗u)|, i.e., the spectral norm of the tensor Rn, see, e.g., [19, Lemma 1]. Fix next ε>0andletNnε ⊆Spn−1be an ε-net of Spn−1. That is, for each u∈Spn−1there exists v∈N nε such that u−v≤ε.From([62], Lemma 5.2) we know that Nnε can be chosen such that its cardinality |Nnε|is at most (1+2/ε)pn.Fixnowu∈Spn−1and let v≡vube its ε-neighbour in Nnε. Then, we have, |Rn×(u⊗u⊗u)−Rn×(v⊗v⊗v)| ≤|Rn×(u⊗u⊗{u−v})| +|Rn×(u⊗{u−v}⊗v)| +|Rn×({u−v}⊗v⊗v)| ≤3εRn2, (B.22) where the final inequality follows as |Rn×(a1⊗a2⊗a3)|≤R n2a1a2a3 for any a1,a2,a3∈Rpn,see[19]. By the reverse triangle inequality, (B.22)gives, |Rn×(u⊗u⊗u)|≤|R n×(v⊗v⊗v)|+3εRn2 ≤max w∈Nnε |Rn×(w⊗w⊗w)|+3εRn2. Since the above holds for all u∈Spn−1, we further get, Rn2=sup u∈Spn−1|Rn×(u⊗u⊗u)|≤ max w∈Nnε |Rn×(w⊗w⊗w)|+3εRn2, that is, Rn2≤1 1−3εmax u∈Nnε |Rn×(u⊗u⊗u)|, whenever ε<1/3. We next apply this bound with the choice ε=2/9toget P(Rn2≥t)≤P( max u∈Nn(2/9) |Rn×(u⊗u⊗u)|≥t/3) ≤ u∈Nn(2/9) P(|Rn×(u⊗u⊗u)|≥t/3).(B.23) 6724 U. Radojiˇci´cetal. Now, for a fixed u∈N n(2/9), the quantity Rn×(u⊗u⊗u) is distributed as 1 nn i=1 X3 i−E(X3)whereX1,...,X nis a random sample from the distribution of X∼α1N(−α2uhn,uΣnu)+α2N(α1uhn,uΣnu). Consequently, Lemma B.5 implies that P(|Rn×(u⊗u⊗u)|≥t/3) ≤2exp−√nt 3K{1+(Σn2+hn2)3/2}2/3 whenever t≥3K[1 + {uΣnu+(uhn)2}3/2]n−1/2. I.e., in particular when t≥ 3K{1+(Σn2+hn2)3/2}n−1/2.As|Nnε|≤10pn, plugging in to (B.23), we finally get, P(Rn2≥t)≤2exp−√nt 3K{1+(Σn2+hn2)3/2}2/3 +pnlog 10, for all t≥3K{1+(Σn2+hn2)3/2}n−1/2,whichis,forafixedt>0, satisfied for all large enough n, thanks to our assumptions that pn→∞and pn(Σn2+ hn2)=o(n1/3). Denoting Hn:= Σn2+hn2, the same assumptions, in conjunction with the requirement that Hn≥Cfor all nlarge enough, guarantee that we have, for each fixed t>0, that −n1/3 Hnt2/3 {3K}2/3{H−3/2 n+1}2/3+pnHn n1/3log 10→−∞, as n→∞. Thus the claim follows. We next show equivalent results for the second moment instead of the third, beginning with the finiteness of the ψ1-Orlicz “norm” of X2−E(X2) for univariate normal mixtures. The proof of this is exactly analogous to that of Lemma B.4 andweomitit. Lemma B.6. Let X∼α1N(−α2h, σ2)+α2N(α1h, σ2)where h∈R,σ2>0. Then, X2−E(X2)ψ1≤10(σ2+h2). Lemma B.6 allows us to derive a concentration bound for the sample-centered second moment. Again, the proof exactly follows Lemma B.5, causing us to leave it out. Lemma B.7. Let X1,...,X nbe a random sample from the distribution of X∼ α1N(−α2h, σ2)+α2N(α1h, σ2)where h∈R,σ2>0.Then, P 1 n n  i=1 X2 i−E(X2)≥ε≤2exp−√nε K(1 + σ2+h2), for all ε≥K(1 + σ2+h2)n−1/2where K>0is a constant not depending on any of the parameters. Properties of PP-based estimators of the linear discriminant 6725 Finally, Lemmas B.6 and B.7 give us the uniform law of large numbers for the second sample moment over all projections. Theorem B.2. Let xn1,...,xnn be a random sample from the model (2)and assume that pn→∞ and pn(Σn2+hn2)=o(n1/2). Assume further that for some C>0, we have Σn2+hn2>Cfor all n large enough. Then, sup u∈Spn−1 1 n n  i=1 (uxni)2−E{(uxn)2}→p0, as n→∞. ProofofTheoremB.2.We provide only the main steps of the proof (which is very similar to that of Theorem B.1). Letting Sn:= (1/n)n i=1 xnix ni − E{xnx n}, our claim is that Sn2→p0. Arguing as in Theorem B.1 we obtain, Sn2≤1 1−2εmax u∈Nnε |uSnu|, whenever ε<1/2, where Nnε is an ε-net of Spn−1.Thus, P(Sn2≥t)≤ u∈Nn(1/4) P(|uSnu|≥t/2),(B.24) where the right-hand side probabilities each satisfy, by Lemma B.7, P(|uSnu|≥t/2) ≤2exp−√nt 2K(1 + Σn2+hn2), when t≥2K(1 + Σn2+hn2)n−1/2. Plugging in to (B.24), we get, P(Sn2≥t)≤2exp−√nt 2K(1 + Σn2+hn2)+pnlog 9, for all t≥2K(1 + Σn2+hn2)n−1/2, a condition which is, for a fixed t>0, satisfied for all large enough nby our assumptions. The desired result now follows. Auxiliary results for the proof of Theorem 5,part2 In this subsection, we prove the two main parts required in the M-estimator argument, i.e., uniform identifiability and uniform law of large numbers. However, before them, we first give a lemma that quantifies the extent to which invertible linear transformations preserve the non-parallelity of vectors. 6726 U. Radojiˇci´cetal. Lemma B.8. Assume that a,b∈Rp,a,b=0, are such that  ab ab≤1−ε, for some ε>0. Moreover, let M∈Rp×psatisfy max{M2,M−12}≤Cfor some C>0.Then,  aM2b MaMb≤1−εC−4. Proof of Lemma B.8.We have aM2b MaMb=1−1 2    Ma Ma−Mb Mb    2 ≤1−1 2M−12 2    a Ma−b Mb    2 ≤1−1 2C2a Ma−b Mb2 +2εab MaMb ≤1−εC−4, where the last line uses Ma≤M2a≤Caand similarly for Mb. The desired bound for the other direction can be obtained by starting from the equation, −aM2b MaMb=1−1 2    Ma Ma+Mb Mb    2 , and proceeding analogously. Lemma B.9 then next gives sufficient conditions for the model parameters under which the sequence of population level maximizers is identifiable uniformly in n. Lemma B.9. Let xn1,...,xnn be a random sample from the model (2)with α1=α2and assume that there exists C1,C 2>0such that, for all n, 1/C1≤hn≤C1,Σn2≤C2,Σ−1 n2≤C2. Then, for all fixed ε∈(0,C2 2),thereexistsδ≡δ(ε, C1,C 2,α 1)>0such that (a) below implies (b) for all n. (a) u∈Spn−1is such that |uθn|≤1−ε. (b) γ2 n0(θn)−γ2 n0(u)≥δ. Proof of Lemma B.9.By the proof of Lemma 2, γ2 n0(u)=(α1α2)4(1−4α1α2)2fn(u) 1+α1α2fn(u)6 =: (α1α2)4(1−4α1α2)2g6 n(u), Properties of PP-based estimators of the linear discriminant 6727 where fn(u):=(uhn)2/uΣnuand the constant ι:= (α1α2)4(1 −4α1α2)2is strictly positive. Moreover, gn(θn)≥gn(u)≥0, implying that we have γ2 n0(θn)−γ2 n0(u)≥ι{gn(θn)3−gn(u)3}{gn(θn)3+gn(u)3} ≥ι{gn(θn)−gn(u)}g5 n(θn). Now, fn(θn)=h nΣ−1 nhn≥hn2Σn−1 2≥C−2 1C−1 2for all n, implying that there exists C3>0 such that gn(θn)≥C3for all n. Thus, the claim of the lemma holds once we show that gn(θn)−gn(u)≥C4for some C4>0not depending on n. Now, gn(θn)−gn(u)= fn(θn)−fn(u) {1+α1α2fn(θn)}{1+α1α2fn(u)} ≥fn(θn)−fn(u) {1+α1α2fn(θn)}2 ≥fn(θn)−fn(u) {1+α1α2hn2Σ−1 n2}2 ≥fn(θn)−fn(u) {1+α1α2C2 1C2}2, showing that it remains to lower bound fn(θn)−fn(u) is a manner not depending on n. Some algebra reveals that, fn(u)=uΣnθn Σ1/2 nuΣ1/2 nθn2 h nΣ−1 nhn≤(1 −εC−2 2)2h nΣ−1 nhn, where the inequality follows from applying Lemma B.8 with a=u,b=θnand M=Σ1/2 n. Consequently, fn(θn)−fn(u)≥C−2 1C−3 2ε(2 −εC−2 2)≥C−2 1C−3 2ε. The lower bound does not depend on n, finally implying the desired claim. Before showing the uniform law of large numbers for our objective function, we first establish an auxiliary result on some basic properties of op(1) and Op(1) sequences of random variables. Lemma B.10. Assume that a sequence of random variables Yn>0has Yn≥ ε+Rnfor some ε>0and some sequence of random variables Rn=op(1). Then, 1 Yn≤1 ε+Tn, for some sequence of random variables Tn=op(1). 6728 U. Radojiˇci´cetal. Proof of Lemma B.10.We first show that 1/Yn=Op(1). To see this, we write, P(1/Yn>2/ε)=P(Yn<ε/2) ≤P(ε+Rn<ε/2) = P(Rn<−ε/2), where the final probability goes to zero as n→∞. Hence 1/Yn=Op(1), and we can write, 1 Yn−1 ε=ε−Yn Ynε≤−Rn Ynε=op(1)Op(1) = op(1), concluding the proof. Lemma B.11. Let xn1,...,xnn be a random sample from the model (2)and assume that there exists C1,C 2>0such that, for all n, hn≤C1,Σn2≤C2,Σ−1 n2≤C2. Assume further that, pn→∞ and pn=o(n1/3). Then, sup u∈Spn−1γ2 n(u)−γ2 n0(u)→p0, as n→∞. Proof of Lemma B.11.We use the notation rn0(u):=E{(uxn)3},rn(u):= (1/n)n i=1(uxni)3,sn0(u):=E{(uxn)2}and sn(u):=(1/n)n i=1(uxni)2. Now, rn0(u) is uniformly upper bounded in uand n, as can be seen by writing |rn0(u)|=|α1α2(α1−α2)(uhn)3|≤C3 1. Similarly, we have |sn0(u)|=uΣnu+ α1α2(uhn)2≥Σ−1 n−1 2≥C−1 2and |sn0(u)|≤Σn2+hn2≤C2+C2 1. And, by writing, |sn0(u)|≤|sn0(u)−sn(u)|+|sn(u)|≤ sup u∈Spn−1|sn0(u)−sn(u)|+|sn(u)|, Theorem B.2 further gives that |sn(u)|≥|sn0(u)|+qn≥C−1 2+qn,where qn=op(1). We also note that this implies that |sn(u)|3≥(C−1 2+qn)3= C−3 2+op(1), as taking the third power is an increasing mapping. This further gives 1/|sn(u)|3≤C3 2+op(1), almost surely, by Lemma B.10, as the second moment sn(u) is almost surely positive for all large enough n(in the sequel, we implicitly restrict to this almost sure set). We now write, sup u∈Spn−1γ2 n(u)−γ2 n0(u) ≤sup u∈Spn−1 1 s3 n0(u)|rn(u)−rn0(u)||rn(u)+rn0(u)| +sup u∈Spn−1 r2 n(u) 1 s3 n(u)−1 s3 n0(u) , (B.25) Properties of PP-based estimators of the linear discriminant 6729 and treat the two supremums on the right-hand side separately. Denoting Bn:= supu∈Spn−1|rn(u)−rn0(u)|, the first one has the upper bound, C3 2Bn(Bn+2C3 1), which is of the order op(1) by Theorem B.1. Whereas, the second supremum has the upper bound, {Bn(Bn+2C3 1)+C6 1}sup u∈Spn−1 s3 n0(u)−s3 n(u) s3 n(u)s3 n0(u) ≤{Bn(Bn+2C3 1)+C6 1}C3 2sup u∈Spn−1 s3 n0(u)−s3 n(u) s3 n(u) ≤{Bn(Bn+2C3 1)+C6 1}C3 2{C3 2+op(1)}sup u∈Spn−1|s3 n0(u)−s3 n(u)|. (B.26) Denoting, En:= sup u∈Spn−1|sn(u)−sn0(u)|, the left-over supremum in the final expression of (B.26) has the upper bound, Ensup u∈Spn−1|s2 n0(u)+sn0(u)sn(u)+s2 n(u)| ≤En{3(C2+C2 1)2+ 3(C2+C2 1)En+E2 n}, which converges in probability to zero by Theorem B.2. Hence, both terms on the right-hand side of (B.25) are upper bounded by sequences of random variables converging in probability to zero, and we finally obtain the desired claim. Proof of Theorem 5 ProofofTheorem5.We first note that a sequence of maximizers exists almost surely for all nlarge enough as γ2 nis a rational function whose denominator is a positive power of a quantity of the form u{(1/n)n i=1 xnix ni}uwhere the matrix n i=1 xnix ni is almost surely positive-definite when pn≤n(again, we restrict implicitly to the corresponding almost sure set). Fix now ε∈(0,C2 2). Then, by Lemma B.9,wehave, P(1 −|u nθn|≥ε)≤P{γ2 n0(θn)−γ2 n0(un)≥δ}, for some δ≡δ(ε, C1,C 2,α 1)>0. Consequently, by the triangle inequality, P(1 −|u nθn|≥ε)≤P{|γ2 n0(θn)−γ2 n(θn)|≥δ/3} +P{γ2 n(θn)−γ2 n(un)≥δ/3} +P{|γ2 n(un)−γ2 n0(un)|≥δ/3}.(B.27) The first and the third term on the right-hand side in (B.27) are by Lemma B.11 o(1), while the second term equals zero as unis a maximizer of u→ γ2 n(u). Hence, the claim follows after noting that we always have |u nθn|≤1 due to the unit lengths of unand θn. 6730 U. Radojiˇci´cetal. Appendix C: Additional simulation results In this section, we give supporting plots as supplementary material to claims made and plots presented in the article. Simulations and the corresponding plots are done using R 3.6.1 [59] together with R packages ICtest [49], mvtnorm [22], MASS [61], GGally [56], ggpubr [31], dplyr [67], tidyr [66] and RColorBrewer [47]. Figures C.1 and C.2 show the standard deviation of maximal similarity index snu nθ/θwhere unis one of PP estimators discussed in the article, as a function of the Mahalanobis distance between the group means τand mixing proportion α1, for sample sizes n∈{500,1000,2000,4000}and n∈ {8000,16000,32000}respectively. Fig C.1. The heatmaps show standard deviation of the MSI snu nθ/||θ|| as a function of Mahalanobis distance between the group means τand mixing proportion α1,whereunis one of estimators of θ/||θ|| obtained by maximizing (κn−3)2,γ2 n,ηn(·;0.8) and ηn(·;w1), where in the latter case, w1=w1(α1,τ)maximizes asymptotic relative efficiency of hybrid estimator w.r.t. LDA, for sample sizes n∈{500,1000,2000,4000}.Signsnis chosen such that snu nθ/||θ|| ≥ 0. In each setting, mean is calculated using m= 1000 replicates and the data is randomly generated from 10-dimensional normal mixtures with covariance matrix Σ having AR(1) structure with ρ=0.6., while μ1=0and μ2chosen in each setting such that μ 2Σμ2=τ. Properties of PP-based estimators of the linear discriminant 6731 Fig C.2. The heatmaps show standard deviation of the MSI snu nθ/||θ|| as a function of Mahalanobis distance between the group means τand mixing proportion α1,whereunis one of estimators of θ/||θ|| obtained by maximizing (κn−3)2,γ2 n,ηn(·;0.8) and ηn(·;w1), where in the latter case, w1=w1(α1,τ)maximizes asymptotic relative efficiency of hybrid estimator w.r.t. LDA, for sample sizes n∈{8000,16000,32000}.Signsnis chosen such that snu nθ/||θ|| ≥ 0. In each setting, mean is calculated using m= 1000 replicates and the data is randomly generated from 10-dimensional normal mixtures with covariance matrix Σ having AR(1) structure with ρ=0.6, while μ1=0and μ2chosen in each setting such that μ 2Σμ2=τ. Figure C.3 shows mean of the maximal similarity index snu nθ/θwhere un is one of PP estimators discussed in the article, as a function of the Mahalanobis distance between the group means τand mixing proportion α1, for large sample sizes, n∈{8000,16000,32000}. 6732 U. Radojiˇci´cetal. Fig C.3. The heatmaps show means of the MSI snu nθ/||θ|| as a function of Mahalanobis distance between the group means τand mixing proportion α1,whereunis one of estimators of θ/||θ|| obtained by maximizing (κn−3)2,γ2 n,ηn(·;0.8) and ηn(·;w1), where in the latter case, w1=w1(α1,τ)maximizes asymptotic relative efficiency of hybrid estimator w.r.t. LDA, for sample sizes n∈{500,1000,2000,4000}.Signsnis chosen such that snu nθ/||θ|| ≥ 0.In each setting, mean is calculated using m= 1000 replicates and the data is randomly generated from 10-dimensional normal mixtures with covariance matrix Σhaving AR(1) structure with ρ=0.6,whileμ1=0and μ2chosen in each setting such that μ 2Σμ2=τ. Figure C.4 shows a scatter matrix plot of the finance data set from the Rpackage Rmodmix, where the point in the plot is being colored red if the company is being bankrupt, and blue otherwise, as well as the marginal densities for both groups which are given at the diagonal. Properties of PP-based estimators of the linear discriminant 6739 Processing Systems (C. Cortes,N. Lawrence,D. Lee,M. Sugiyama and R. Garnett,eds.)28 1567–1575. Curran Associates, Inc. [70] Zhao, Y.,Shrivastava, A. K. and Tsui, K.-L. (2019). Regularized Gaussian mixture model for high-dimensional clustering. IEEE Transactions on Cybernetics 49 3677-3688. [71] Zhu, R.,Wang, L.,Zhai, C. and Gu, Q. (2017). High-dimensional variance-reduced stochastic gradient Expectation-Maximization algorithm. In Proceedings of the 34th International Conference on Machine Learning (D. Precup and Y. W. Teh,eds.).Proceedings of Machine Learning Research 70 4180–4188. PMLR.