Theoretical Insights Into Multiclass Classification: A High-dimensional Asymptotic View

Christos Thrampoulidis, Samet Oymak, Mahdi Soltanolkotabi

Introduction

Multiclass classification is fundamental to a large number of real-world machine learning applications that demand the ability to automatically distinguish between thousands of different classes. Applications include essentially any problem with categorical outputs spanning natural language processing [SVL14], where a seq2seq decoder has to choose the correct word token, reinforcement learning [JGP16, MXSS20], where the agent has to choose the correct action, to recommendation systems, where the model should recommend the correct movie out of many other options. For instance, YouTube’s recommendation system is modeled as an extreme multiclass problem with more than a million classes where each video corresponds to a viable class [CAS16].

The growing list of applications motivate an in-depth exploration of multiclass classification algorithms. Despite their extensive use however, a precise understanding of the statistical properties and behavior of classification algorithms is still missing with many open questions: What is the total and per class test accuracy? How does this quantity depend on various problem parameters such as data distributions, problem dimensions, etc.? What is the highest test accuracy achievable by any algorithm? What is the best algorithm for each scenario? Which algorithm achieves the highest accuracy on rare or minority classes? How does the answer to the above question change in modern regimes where the number of classes is large?

Asymptotic analysis in modern high-dimensional regimes where the number of training data and feature sizes grow in tandem with each other provides a promising setting for precisely quantifying the accuracy of classification algorithms as a function of problem variables and resolving the questions above. However, despite the rich literature on precise high-dimensional estimation and more recently binary classification, multiclass classification is an under-explored venue possibly due to the difficulty of capturing the intricate dependencies between the classes even for relatively simple linear classifiers.

Contributions. We initiate a precise asymptotic study of linear multiclass classification in the modern high-dimensional regime, where the sizes of the training data and of the feature vectors grow large at a proportional rate. A key promise of such a precise analysis is that it allows us to accurately compare between different classification algorithms and data models. Compared to linear regression/binary classification, we identify the following crucial challenge: the test accuracy in multiclass classification relies on intricate cross-correlations between the trained weights of the classifier. This has two consequences that drive our analysis. First, in order to obtain sharp asymptotics on the test error of any classifier, it is a prerequisite to precisely quantify the asymptotics of these cross-correlations. Second, the test error does not depend on the correlations in closed-form expressions. Thus, to compare between different classifiers, we need efficient numerical and analytic means to evaluate the test error in terms of the correlation matrices. Interestingly, we show that these challenges are already present in simple classifiers, such as minimizing the square loss, and in stylized distributional settings, such as Gaussian features. Our contributions are as follows:

∙\bullet We study two different data models: a Gaussian Mixtures Model (GMM) and a Multinomial Logit Model (MLM) with Gaussian features. For each one of them, we provide a precise characterization of total and class-wise test accuracy for three different training algorithms: (i) a least-squares (LS) based classifier, (ii) a weighted least-squares (WLS) based classifier, and (iii) a simple per class averaging (Avg) estimator. For the least-squares based classifiers, we develop a new technique to overcome the technical challenge of characterizing the limiting behavior of the weights’ cross-correlations. For the per class averaging classifier, we show that it is Bayes optimal for a GMM with equal priors.

∙\bullet We discuss efficient means of evaluating the test accuracy as a function of the weights’ cross-correlations. This, together with the derived asymptotic formulae for the latter, lead to the first precise high-dimensional characterization of how the total/class-wise accuracy varies for different algorithms, data distributions, problem dimensions as well as number of classes, the inter/intra class correlations and class priors. For special problem geometries, we derive precise conditions on the data distribution and on the relative size of the training set over which each of the two studied algorithms dominates.

∙\bullet We present and discuss numerical simulations that corroborate our theoretical findings. For instance, with an eye towards making classification algorithms more fair/equitable, we use our precise characterization of the class-wise accuracy to demonstrate how different algorithms behave in the presence of rare/minority classes. We also empirically compare the algorithms studied in this paper to other popular losses such as cross-entropy minimization. This allows us to better understand the performance of various algorithms in modern regimes of large number of classes.

Related Work. There is a classical body of algorithmic work on multiclass classification, e.g., [CS01, LLW04, WW98, BB99, DB94] and several empirical studies of their comparative performance [RK04, Für02, ASS00, PM05]. A more recent extension of this line of work investigates the effect of the loss function in deep neural networks [HYS16, GCOZ17, KS18, BEH20, DCO20]. Algorithms for extreme multiclass problems with huge number of classes has also been studied in several [CAL13, YHR+16, RCY+19, KMS15] works. On the theory front, numerous works have investigated consistency [Zha04, LLW04, TB07, PSG13, PS16] and finite-sample behavior [KP+02, Gue02, ASS00, LLY+18, CKMY16, LDBK15, Mau16, LDZK19] of multiclass classification algorithms. Our work differs from this literature in that we are interested in precise characterizations of the test accuracy rather than order-wise bounds. Here we focus on linear classifiers, but we consider the modern high-dimensional regime in which both the sample size and the features’ dimension are large.

Specifically, our theoretical approach to linear multiclass classification fits in the rapidly growing literature on sharp high-dimensional asymptotics of convex optimization-based estimators [Don06, Sto09, OH10, CRPW12, ALMT13, DMM11, BM12, ALMT13, Sto13, OTH13, TOH15, Kar13, EK18, DM16, ORS17, TXH18, TAH18, MM18, WWM19, CM19, HL19, BKRS19, ASH19, JSH20]. Most of this line of work studies linear models and regression problems. More recently there has been a surge of interest in sharp analysis of a variety of methods tailored to binary classification models [TAH15, Hua17, CS18, SC19, MLC19b, MLC19a, KA20, SAH19, TPT20b, DKT19, MRSY19, LS20, KT20, MKLZ20, Lol20, TPT20a]. Nevertheless, none of these prior works have yet considered multiclass classification settings. Our paper unveils the salient features of the multiclass setting and shows that corresponding results from the binary setting do not directly apply here. We emphasize that this is the case even for seemingly simple one-vs-all (OVA) classifiers, such as minimizing the square-loss, that involve training a single binary classifier per class [RK04]. The key technical tool behind our sharp analysis is the convex Gaussian min-max Theorem (CGMT) [TOH15, Sto13]. However, a “naive" application of the CGMT on the original optimization of the classifier does not allow us to compute all the necessary correleations between the classfier’s weights to precisely capture the total/class-wise errors. Instead, our key idea is to formulate an artificial optimization problem, which captures the missing correlations and at the same time conveniently allows us to leverage the CGMT.

Problem formulation

Multinomial Logit Model (MLM). In this model we assume that feature vectors x\bm{x} are distributed i.i.d. N(0,Id)\mathcal{N}(\bm{0},\bm{I}_{d}) and that the conditional density of the class labels is given by the soft-max function. Concretely, we say that a data point (x,Y)(\bm{x},Y) (or its one-hot encoded representation (x,y)(\bm{x},\bm{y})) follows the multinomial logit model when

2 Classification algorithms

Least-squares (LS). In this approach we train a linear classifier x↦Wx+b\bm{x}\mapsto\bm{W}\bm{x}+\bm{b} via a least-squares fit to the training data:

Class averaging (Avg). This approach uses the following weight and offset values

Weighted Least-squares (WLS). This is a variation of the Least-squares approach where we fit a weighted least squares loss of the form

3 Class-wise and total test classification error

Let W^,b^\widehat{\bm{W}},\widehat{\bm{b}} denote the parameters of a trained classifier. Now consider a fresh data sample (x,Y)(\bm{x},Y) generated according to the same distribution as the training data. Once, we have learned the parameters W^,b^\widehat{\bm{W}},\widehat{\bm{b}} of the classifier, the class Y^\widehat{Y} predicted by the classifier is made by a winner takes it all strategy, as follows, Y^=arg⁡max⁡j∈[k] ⟨w^j,x⟩+b^j.\widehat{Y}={\arg\max}_{j\in[k]}\text{ }\langle\widehat{\bm{w}}_{j},\bm{x}\rangle+\widehat{\bm{b}}_{j}. Therefore, the classification error condition on the the true label being cc, which we shall refer to as the class-wise test error, is equal to

Correspondingly, the total classification error is given by

GMM. In model (2.1), the test error probability is explicitly given by

MLM. In model (2.2), the test error probability is explicitly given by

Calculating the class-wise/total misclassifcation errors. The identities (2.5) and (2.6) (see Section D.1 for a proof) as well as similar ones for the class-wise test error demonstrate that the total/class-wise errors only depend on the correlation matrices Σw,w\bm{\Sigma}_{\bm{w},\bm{w}} and Σw,μ\bm{\Sigma}_{\bm{w},\bm{\mu}}, the offset values b^\widehat{\bm{b}} and the the class conditional means. For instance, as we show in the supplementary for GMM the class-wise errors are given by

4 High-dimensional regime

This paper derives sharp asymptotic formulae for the class-wise and total classification error of averaging and (weighted) LS algorithms for GMM and MLM. We defer all our proofs to the appendix. All our results hold in the following high-dimensional regime with finite kk.

We focus on a double asymptotic regime where n,d→∞n,d\rightarrow\infty at a fixed ratio γ=d/n>0\gamma=d/n>0.

For the (weighted) least-squares classifier, we focus here in the overdetermined regime γ<1\gamma<1. However, our approach is also directly applicable to regularized (or min-norm) LS/WLS in the overparameterized regime γ>1\gamma>1.

For a sequence of random variables Xn,d\mathcal{X}_{n,d} that converges in probability to some constant cc in the limit above, we simply write Xn,d⟶Pc\mathcal{X}_{n,d}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}c. For a random vector/matrix vn,d\bm{v}_{n,d}/Vn,d{\bm{V}}_{n,d} and a deterministic vector/matrix c\bm{c}/C{\bm{C}}, the expressions vn,d⟶Pc\bm{v}_{n,d}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\bm{c} and Vn,d⟶PC{\bm{V}}_{n,d}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}{\bm{C}} are to be understood entry-wise.

Results for Gaussian Mixture Model

Consider data generated according to GMM in an asymptotic regime with any γ>0\gamma>0. For the averaging estimator discussed in Section 2.2, the following high-dimensional limits hold

The above result allows us to precisely characterize the behavior of the averaging estimator in the high-dimensional regime. Let us consider a few special cases.

after some algebraic manipulations the total classification error of the averaging estimator in this case is given by

where G0,…,Gk−1∼iidN(0,1).G_{0},\ldots,G_{k-1}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,1).

2 Least-squares classifier

This section focuses on characterizing the intercepts and correlation matrices for the least-squares classifier. To present our results, we assume that the Grammian matrix has eigenvalue decomposition

with Σ\bm{\Sigma} a diagonal positive-definite matrix and V{\bm{V}} an orthonormal matrix obeying VTV=Ir{\bm{V}}^{T}{\bm{V}}={\bm{I}}_{r}.

Consider data generated according to GMM in an asymptotic regime with γ<1\gamma<1. In addition to (3.2), define the following two positive (semi)-definite matrices: P:=diag(π)−ππT⪰0k×k{\bm{P}}:=\text{diag}(\bm{\pi})-\bm{\pi}\bm{\pi}^{T}\succeq\bm{0}_{k\times k} and Δ:=σ2Ir+ΣVTPVΣ≻0r×r.\bm{\Delta}:=\sigma^{2}{\bm{I}}_{r}+\bm{\Sigma}{\bm{V}}^{T}{\bm{P}}{\bm{V}}\bm{\Sigma}\succ\bm{0}_{r\times r}. Then, for the least-squares linear classifier (W^,b^)\left(\widehat{\bm{W}},\widehat{\bm{b}}\right) the following limits are true asymptotically

The above result allows us to precisely characterize the behavior of the least-squares classifier in the high-dimensional regime. In Section G.2, we specialize (3.3) to the case of orthogonal means. Compared to the weight vectors w^i,i∈[k]\widehat{\bm{w}}_{i},i\in[k] of the class averaging classifier that are also (asymptotically) orthogonal when means are orthogonal, this is not the case for LS. We show next that these spurious correlations only hurt the classification error when classes are balanced.

Consider the case of orthogonal, equal energy-means Σμ,μ=μIk\bm{\Sigma}_{\bm{\mu},\bm{\mu}}=\mu{\bm{I}}_{k}, balanced priors πi=1/k, i∈[k]\pi_{i}=1/k,~{}i\in[k] and γ<1\gamma<1. Setting uLS:=μ2σ1−γμ2+kγσ2 ,u_{\rm LS}:=\frac{\mu^{2}}{\sigma}\sqrt{\frac{1-\gamma}{\mu^{2}+k{\gamma\sigma^{2}}}}\,, it holds that

Specifically, since uLS=uAvg1−γ < uAvgu_{\rm LS}=u_{\rm Avg}\sqrt{1-\gamma}\,<\,u_{\rm Avg}, the averaging estimator strictly outperforms LS for all 0<γ<10<\gamma<1 and k≥2k\geq 2 in this setting.

3 Bayes estimator for the balanced Gaussian Mixture Model

Results for Multinomial Logit Model

Note that π\bm{\pi} and Π\bm{\Pi} are the first and second moments of the soft-max mapping of VΣg∼N(0,Σμ,μ){\bm{V}}\bm{\Sigma}{\bm{g}}\sim\mathcal{N}\left(\bm{0},\bm{\Sigma}_{\bm{\mu},\bm{\mu}}\right). In fact, for the MLM in (2.2) it holds that

since MTx{\bm{M}}^{T}\bm{x} is distributed as VΣg{{\bm{V}}\bm{\Sigma}{\bm{g}}}. Thus, π\bm{\pi} is the vector of class priors (which explains the slight abuse of notation here in relation to our notation for the class priors of the GMM).

Consider data generated according to MLM in an asymptotic regime with any γ>0\gamma>0. For the averaging classifier, the following high-dimensional limits hold

Using Gaussian decomposition in (2.6) and checking from (4.2) that Σw,w−Σw,μΣμ,μ†Σw,μT⟶Pγ⋅diag(π)\bm{\Sigma}_{\bm{w},\bm{w}}-\bm{\Sigma}_{\bm{w},\bm{\mu}}\bm{\Sigma}_{\bm{\mu},\bm{\mu}}^{\dagger}\bm{\Sigma}_{\bm{w},\bm{\mu}}^{T}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\gamma\cdot\text{diag}(\bm{\pi}) the test error obtains the following explicit form:

2 Least-squares classifier

This section focuses on characterizing the intercepts and correlation matrices for the least-squares classifier. We also use the result to characterize conditions under which LS outperforms averaging.

Consider data generated according to MLM in an asymptotic regime with 0<γ<10<\gamma<1. Recall the notation in (4.1). For the LS classifier, the following high-dimensional limits hold.

It is interesting to observe that (4.4a) is identical to (4.2a). However, the cross-correlations in Σw,w\bm{\Sigma}_{\bm{w},\bm{w}} differ. We prove below that this leads to an improved performance of the LS classifier for large sample sizes. First, Theorem 4.2 can be used to check that

Thus, the only change in the test-error formula compared to (4.3) is the term γ⋅diag(π)\gamma\cdot\text{diag}(\bm{\pi}) substituted by the matrix above.

Assume orthogonal, equal-energy means Σμ,μ=μ2Ik\bm{\Sigma}_{\bm{\mu},\bm{\mu}}=\mu^{2}{\bm{I}}_{k}, k≥2k\geq 2. Let

Numerical Results

Proof outline for least-squares: key ideas and challenges

In this section, we provide a proof sketch for the analysis of the multiclass least-squares (LS) classifier.

Specifically, we discuss our approach towards specifying the high-dimensional limits of the key quantities needed to evaluate the classification error: b,Σw,μ,\bm{b},\bm{\Sigma}_{\bm{w},\bm{\mu}}, and, Σw,w\bm{\Sigma}_{\bm{w},\bm{w}}. For simplicity, we focus here on the performance of the LS classifier GMM. We note that our proofs for the MLM and the Weighted Least-Squares (WLS) classifiers follow the same general strategy, but in some parts require more involved and intricate analysis and derivations. Our proof follows the following general steps; see the appendix for complete details and derivations.

Step I: Decomposing the loss across classes. Recall from Section 2.2 that the multiclass LS classifier produces a linear classifier x↦Wx+b\bm{x}\mapsto\bm{W}\bm{x}+\bm{b} via a least-squares fit to the training data:

Notice that the objective function above is separable. That is,

Step II: Reduction to an Auxiliary Optimization (AO) problem via CGMT. To calculate the high-dimensional statistical behavior of (6.2) we use the Convex Gaussian min-max Theorem (CGMT) [Sto13, TOH15] framework. We provide a brief introduction of the CGMT machinery in Section 6.1. Roughly stated, this framework allows us to replace a Primary Optimization (PO) problem of the form (6.2) with an Auxiliary Optimization (AO) problem that is simpler to analyze, but is predictive of the behavior of the latter. For instance, for the PO in (6.2) in the GMM, after some algebraic manipulations, the AO problem takes the form

where Δ:=σ2Ir+ΣVTPVΣ≻0r×r\bm{\Delta}:=\sigma^{2}{\bm{I}}_{r}+\bm{\Sigma}{\bm{V}}^{T}{\bm{P}}{\bm{V}}\bm{\Sigma}\succ\bm{0}_{r\times r} and P:=diag(π)−ππT{\bm{P}}:=\text{diag}(\bm{\pi})-\bm{\pi}\bm{\pi}^{T}.

The CGMT is an extension of Gordon’s Gaussian min-max inequality (GMT) [Gor88]. In the context of high-dimensional inference problems, Gordon’s inequality was first successfully used in the study oh sharp phase-transitions in noiseless Compressed Sensing [Sto09, CRPW12, ALMT13, Sto09]. More recently, [Sto13] (see also [ALMT13, Sec. 10.3]) discovered that Gordon’s inequality is essentially tight for certain convex problems. A concrete and general formulation of this idea was given by [TOH15] and was called the CGMT.

In order to summarize the essential ideas, consider the following two Gaussian processes:

In [TAH18], the authors introduce a principled machinery that allows to (a) express a quite general family of convex inference optimization problems in the form of the PO and (b) properly analyze the corresponding AO. In particular, the analysis of the AO is performed in three intermediate steps. First, the (random) optimization over vector variables is simplified to an easier optimization over only few scalar variables, termed the “scalarized AO". After the scalarization step, it is possible to establish (uniform) convergence of the scalarized AO to a deterministic min-max optimization problem over only a few scalar variables. The convergence step is followed by the analysis of the latter deterministic problem, which leads to the desired asymptotic characterizations. Our proofs outlined in Section 6 follow this general strategy, but the new idea introduced in Step IV therein is key to capture the asymptotic behavior of the off-diagonal entries of Σww\bm{\Sigma}_{\bm{w}\bm{w}}.

Future Directions

This work aims at initiating a precise asymptotic study of multiclass classifiers that provides a promising setting for resolving a rich set of open questions regarding the (comparative) performance of classification algorithms as a function of the involved problem variables. As mentioned, even understanding the statistical performance of one-vs-all multiclass classifiers does not follow directly from the existing literature on binary classifiers. Extending the results of this paper to the one-vs-all logistic and SVM classifiers would allow for a principled comparison among these different choices. A possibly more challenging, albeit mathematically intriguing and practically relevant task, is characterizing the asymptotics of more complicated (non-separable) losses, such as the cross-entropy loss. For this, even characterizing the asymptotic behavior of the correlations Σw,μ\bm{\Sigma}_{\bm{w},\bm{\mu}} requires new ideas. The previously mentioned study of “extreme multiclass classification" in which the number of classes kk is very large is another fascinating direction.

Acknowledgments

C. Thrampoulidis is partially supported by the NSF under Grant Numbers CCF-2009030 and HDR-1934641. S. Oymak is partially supported by the NSF award CNS-1932254. M. Soltanolkotabi is supported by the Packard Fellowship in Science and Engineering, a Sloan Research Fellowship in Mathematics, an NSF-CAREER under award #1846369\#1846369, the Air Force Office of Scientific Research Young Investigator Program (AFOSR-YIP) under award #\#FA9550−18−1−00789550-18-1-0078, DARPA Learning with Less Labels (LwLL) and FastNICS programs, and NSF-CIF awards #1813877\#1813877 and #2008443\#2008443.

References

Appendix A Additional Numerical Results

In this section, we provide further numerical experiments.

First, in Figure 5 we investigate the question: When does least-squares provably outperform averaging? Our Proposition 4.3 provides a fundamental transition point in sample complexity above which least-squares is provably better than averaging under MLM. In Figure 5, we visualize γ⋆\gamma_{\star} as a function of different number of classes as well as different levels of mean energy. Least-squares outperform averaging in the region below the lines displayed in Figure 5. Our key message is that least-squares work better when the sample complexity is higher and the problem is less noisy. As the number of classes kk increase, the problem becomes more difficult/noisy and we require a larger sample complexity to ensure classifier achieves a similar amount of accuracy as small kk. Following this intuition, as kk increases, γ⋆\gamma_{\star} shifts smaller due to larger sample requirement. Similarly energy μ\mu directly controls the noise level of the problem, i.e., larger μ\mu results in a larger signal-to-noise ratio. Thus, as we increase μ\mu, γ⋆\gamma_{\star} increases as well because same test accuracy can be achieved with smaller sample size.

Appendix B Additional Results on Weighted Least-squares classifiers

We now focus on characterizing the intercepts/correlation matrices for the WLS classifier.

Surprisingly, the effect of the weights is essentially equivalent to adjusting the class priors from π\bm{\pi} to π~\widetilde{\bm{\pi}} defined in the theorem (modulo the extra additive term in the cross correlation matrix Σw,w\bm{\Sigma}_{\bm{w},\bm{w}}). This shows that weighted LS has similar performance to an un-weighted LS applied to a model with different class priors π~\widetilde{\bm{\pi}}. This characterization allows us to precisely understand how different weighting schemes can alter test accuracy for rare/minority classes.

B.2 WLS for MLM

Theorem B.2 predicts the asymptotic performance of weighted least-squares for data generated according to MLM.

The corresponding formula for the asymptotic limit of the cross-correlation matrix Σw,w\bm{\Sigma}_{\bm{w},\bm{w}} is given in (J.38) in Section J.

Appendix C Preliminaries

In this section we gather a few preliminary results that will be used later on in our proofs.

Let g∼N(0,S){\bm{g}}\sim\mathcal{N}(\bm{0},\bm{S}) and g~∼N(0,R)\widetilde{{\bm{g}}}\sim\mathcal{N}(\bm{0},{\bm{R}}) such that for all i,j∈[k]i,j\in[k]:

Equivalently, letting z∼N(0,Ik){\bm{z}}\sim\mathcal{N}(\bm{0},{\bm{I}}_{k}),

C.2 Gaussian integration by parts

The following result is a direct application of Gaussian integration by parts; for instance, see [FR13, Prop. 8.29].

The following is a corollary of Lemma C.2 applied to the soft-max function.

Let g∼N(0r,Ir){\bm{g}}\sim\mathcal{N}(\bm{0}_{r},{\bm{I}}_{r}) and random vector v=[V1,V2,…,Vk]T\bm{v}=[V_{1},V_{2},\ldots,V_{k}]^{T} with entries:

Further recall the notation of π\bm{\pi} and Π\bm{\Pi} in (4.1). The following statements are true:

C.3 Block matrix inversion

Let T=[AbbTδ]\bm{T}=\begin{bmatrix}{\bm{A}}&\bm{b}\\ \bm{b}^{T}&\delta\end{bmatrix} be an invertible block matrix. Then

where Δ=A−1δbbT≻0\bm{\Delta}={\bm{A}}-\frac{1}{\delta}\bm{b}\bm{b}^{T}\succ\bm{0} is the Schur complement.

Appendix D Calculating and bounding the missclassification error

GMM. Starting from (2.4) and using the fact that xi=μY+z=MeY+z, z∼N(0,σ2Ik)\bm{x}_{i}=\bm{\mu}_{Y}+{\bm{z}}={\bm{M}}\bm{e}_{Y}+{\bm{z}},~{}{\bm{z}}\sim\mathcal{N}(\mathbf{0},\sigma^{2}{\bm{I}}_{k}), we have that

Recall that Σw,μ:=W^M\bm{\Sigma}_{\bm{w},\bm{\mu}}:=\widehat{\bm{W}}{\bm{M}} and note that W^z\widehat{\bm{W}}{{\bm{z}}} is a zero-mean Gaussian vector with covariance matrix σ2W^W^T=σ2Σw,w\sigma^{2}\widehat{\bm{W}}\widehat{\bm{W}}^{T}=\sigma^{2}\bm{\Sigma}_{\bm{w},\bm{w}} in order to conclude with the desired formula in (2.5).

D.2 Class-wise and total miss-classification error for GMM

The class-wise miss-classification error for GMM is given by

where the inequality in the rightmost expression applies entry-wise.

Further, by using the law of total probability we have

D.3 Class-wise and total miss-classification error for MLM

In this section, we derive an explicit formula for the class-wise error for MLM. Recall (2.6):

we can see from (D.5) that the class-wise error probabilities can be calculated as follows:

where πc\bm{\pi}_{c} is the cthc^{\text{th}} entry of the vector π\bm{\pi} in (4.1) and Sc,tc(h)\bm{S}_{c},\bm{t}_{c}(\bm{h}) are defined in (D.6).

D.4 Evaluating and bounding tail probabilities of multivariate Gaussians

In Sections D.3 and D.2, we expressed the class-wise probability of missclassification error for both GMM and MLM in the following convenient form for z∼N(0,Ik−1){\bm{z}}\sim\mathcal{N}(\bm{0},{\bm{I}}_{k-1}),

The formulation above is convenient both in our theoretical analysis, as well as, in simulations. In the rest of this section, we briefly discuss some relevant tools that allow to further simplify or bound expressions in the form of (D.8).

First, we discuss the case where the coefficient matrix A{\bm{A}} and vector t\bm{t} in (D.8) take the special form A∝I+11T{\bm{A}}\propto{\bm{I}}+\bm{1}\bm{1}^{T} and t∝1\bm{t}\propto\bm{1}. This special case appears in some of the stylized symmetric problem settings studied in this paper, such as classification problems with orthogonal and equally-balanced means.

D.4.2 Slepian’s bound

When the matrix A{\bm{A}} does not have the special structure assumed by Lemma D.1, it is not possible in general to provide simple expressions as the one in (D.9). Yet, it might be possible to obtain upper bounds of the same simple form. Such simple bounds can be useful for theoretical interpretations of otherwise complicated formulae, or can provide efficient means for quick (but, non-tight) implementations.

In this section, we discuss Slepian’s inequality (see C.1) as a useful tool in this direction. Assume that a=min⁡i,j∈[k]Aij≥0.a=\min_{i,j\in[k]}{\bm{A}}_{ij}\geq 0. To begin, note that A≥(diag(A)−aI) +a11T,{\bm{A}}\geq\left({\rm diag}({\bm{A}})-a{\bm{I}}\right)\,+a\bm{1}\bm{1}^{T}, where the inequality holds element-wise and equality is true for the diagonal elements. Then, one can apply Slepian’s Lemma C.1 to upper bound the conditional probability of error in (D.8) with the following simple bound:

In the second line above, we used the Gaussian decomposition of Lemma D.1.

D.4.3 Simple bounds for GMM

Union bound. Of course, it is also possible to apply (a simpler) union bound to upper bound the tail probability in (D.8). Here, we show explicitly the result of applying union bound to the class-wise error probabilities of the GMM. Specifically, consider (D.1). An application of the union bound leads to the following:

where in (D.11) Sc,tc\bm{S}_{c},\bm{t}_{c} are defined in (D.3) and in the last line we denote d_{\min}:=\min_{j\neq c}\big{\{}{-[\bm{t}_{c}]_{j}}\Big{/}{\sqrt{[\bm{S}_{c}]_{j,j}}}\big{\}}.

Unfortunately, this bound becomes non-trivial for the class-wise probability of error only if

Intuitively, this assumes a regime wherethe weight vector w^c\widehat{\bm{w}}_{c} corresponding to class cc aligns better with the corresponding mean vector μc\bm{\mu}_{c} than the rest of the weight vectors w^j, j≠c.\widehat{\bm{w}}_{j},~{}j\neq c. This emphasizes the important role of the cross-correlation matrix Σw,w\bm{\Sigma}_{\bm{w},\bm{w}} (including the off-diagonals) for accurate performance prediction. For an illustration, we have implemented this bound and have compared it to our sharp predictions in Figure 7.

Appendix E The Class-averaging estimator

The first statement (3.1a) follows directly from the fact that 1n1TYi=nin⟶Pπi.\frac{1}{n}\bm{1}^{T}{\bm{Y}}_{i}=\frac{n_{i}}{n}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\pi_{i}. For the next two statements note that

where in the last line we used orthogonality of the rows Yj{\bm{Y}}_{j} of the matrix Y{\bm{Y}}:

To conclude simply use the facts that for all i∈[k]i\in[k]:

E.2 Proofs for MLM

Let us define g∼N(0,Ir){\bm{g}}\sim\mathcal{N}(\bm{0},{\bm{I}}_{r}) and random vector v=[V1,V2,…,Vk]T\bm{v}=[V_{1},V_{2},\ldots,V_{k}]^{T} with entries:

We will prove the following three statements:

These lead to (4.2) using Lemma C.3. Therefore, in what follows, we prove (E.4)

where h∼N(0k,Σμμ)\bm{h}\sim\mathcal{N}\left(\bm{0}_{k},\bm{\Sigma}_{\bm{\mu}\bm{\mu}}\right). To deduce the first statement in (E.4a), note that h=\joinrel=(D)VΣg\bm{h}\stackrel{{\scriptstyle\rm(D)}}{{=\joinrel=}}{\bm{V}}\bm{\Sigma}{\bm{g}}.

We proceed similarly with the proof of the last statement in (E.4c) as follows:

Combining the last two displays results in (E.4b), as desired.

E.2.2 Orthogonal means

Specifically, (4.3) can be equivalently expressed as

Appendix F On the Bayes risk of GMM: Proof of Proposition 3.4

Without loss of generality in this proof we assume σ=1\sigma=1. The general result follows by simply replacing (μ,σ)(\mu,\sigma) with (μσ,1)(\frac{\mu}{\sigma},1) and using the proof for σ=1\sigma=1. Recall that the feature vectors x1,…,xn\bm{x}_{1},\ldots,\bm{x}_{n} of the training data set are given by:

To arrive in (F.1) we used that P(y=ej ∣ M,X,Y)=π, ∀j∈[k]P\left(\bm{y}=\bm{e}_{j}~{}|~{}{\bm{M}},{\bm{X}},{\bm{Y}}\right)=\pi,~{}\forall j\in[k] and P(xn+1 ∣ y=ej,M,X,Y)=P(xn+1 ∣ y=ej,M)P\left(\bm{x}_{n+1}~{}|~{}\bm{y}=\bm{e}_{j},{\bm{M}},{\bm{X}},{\bm{Y}}\right)=P\left(\bm{x}_{n+1}~{}|~{}\bm{y}=\bm{e}_{j},{\bm{M}}\right). Also, (F.2) follows by recognizing that P(xn+1 ∣ X,Y)>0P(\bm{x}_{n+1}~{}|~{}{\bm{X}},{\bm{Y}})>0 is independent of the variable of integration M{\bm{M}} and of the optimization variable jj. For the same reasons, in (F.3) we have ignored the normalizing term P(X∣Y)P({\bm{X}}|{\bm{Y}}).

where we denote by Cc\mathcal{C}_{c} the collection of training samples that belong to class c∈[k]c\in[k], i.e. Cc={i∈[n] ∣ yi=ec}.\mathcal{C}_{c}=\{i\in[n]~{}|~{}\bm{y}_{i}=\bm{e}_{c}\}.

With these the objective function of the ML rule in (F.3) becomes:

By completing the squares and invoking a gaussian integral it can be shown that

where ξ(nc):=d/nμ2+ncn, c∈[k]\xi(n_{c}):=\frac{d/n}{\mu^{2}}+\frac{n_{c}}{n},~{}c\in[k] and

For each one of the four terms in (F.8), we have the following by the CLT:

Therefore, in the asymptotic limit, the Bayes estimator is the solution to:

Appendix G Least-squares for GMM

Identifying the AO. To continue further note that by duality we have

Scalarization of the AO. For convenience, define

To continue, consider the singular value decomposition

We also define α=[α1α2…αk]T\bm{\alpha}=\begin{bmatrix}\alpha_{1}&\alpha_{2}&\ldots&\alpha_{k}\end{bmatrix}^{T}. In this notation, we have

Deterministic Analysis. Here, we analyze the deterministic scalar minimization on the RHS of (G.7). Define

First, note that the matrix A{\bm{A}} is positive definite. This can be checked by computing the Schur complement of A{\bm{A}}:

Setting the derivative with respect to α0\alpha_{0} to zero we arrive at

Plugging the latter into (G.12) we arrive at

Asymptotic predictions. First, from (G.11) the bias term converges as follows:

Adding the equations on the above displays we find that

Recognize that this coincides with the optimality condition for (G.16). Thus, the proof is complete.

The analysis of (G.16) is very similar to that of (G.1); thus, most details are omitted. Similar to (G.5) we can relate (G.16) with the following AO problem:

where A{\bm{A}} is as in (G.8) and we have further defined

Thus, similar to (G.11) we can compute the minimizer of the deterministic objective in (G.19):

where recall that Δ\bm{\Delta} is as in (G.10).

Finally, using (G.1.2) and (G.15) in (G.17) it follows that

G.2 Orthogonal means

Here, we specialize the asymptotic predictions of Theorem 3.2 to the case of orthogonal means ⟨μi,μj⟩=0, i≠j\langle\bm{\mu}_{i},\bm{\mu}_{j}\rangle=0,~{}i\neq j.

Then, the following asymptotic limits hold for the least-squares classifier, for all i,j∈[k]i,j\in[k]:

Furthermore, if the means have equal norms μ:=μi\mu:=\mu_{i} and the classes are balanced: πi=1/k, i∈[k]\pi_{i}=1/k,~{}i\in[k], then, setting uLS:=μ2σ1−γμ2+kγσ2u_{\rm LS}:=\frac{\mu^{2}}{\sigma}\sqrt{\frac{1-\gamma}{\mu^{2}+k{\gamma\sigma^{2}}}}, it holds that

Proof This is a direct corollary of Theorem 3.2. Indeed, (G.24) can be derived from (3.3) after substituting V=Ik,Σ=diag(μ1,μ2,…,μk){\bm{V}}={\bm{I}}_{k},\bm{\Sigma}=\text{diag}(\mu_{1},\mu_{2},\ldots,\mu_{k}) and some algebra steps that we omit for brevity.

and applying Lemma (D.1), the probability of error is given by the advertised expression.

Appendix H Least-squares for MLM

Assume that X,Y{\bm{X}},{\bm{Y}} are generated from the MLM.

Identifying the AO. To continue further note that by duality we have

where P⊥=Ir−UUT{\bm{P}}^{\perp}={\bm{I}}_{r}-{\bm{U}}{\bm{U}}^{T}, and we denote

Recalling that xi∼N(0,Id)\bm{x}_{i}\sim\mathcal{N}(\bm{0},{\bm{I}}_{d}) note that

Further recall that for all i∈[n]i\in[n], conditioned on xi\bm{x}_{i}

where we used (H.3) and the SVD decomposition of M{\bm{M}}. In this notation, we can rewrite the PO as follows:

In the remaining, we focus in the inner minimization above. Let us denote

where the expectation is over g∼N(0r,Ir){\bm{g}}\sim\mathcal{N}(\bm{0}_{r},{\bm{I}}_{r}) (with some abuse of notation) and

Next, with an argument based on convexity and compactness similar to that in “Convergence analysis of the AO" in Section G it can be argued that the convergence above is uniform. Thus,

By direct differentiation and first-order optimality, we compute the optimal values as follows:

The analysis of (H.17) is almost identical to the analysis of (H.1) in the previous section. Specifically, without repeating all the details for brevity, it can be shown that the AO of (H.17) converges to the following (cf. (H.11):

where as before G0∼N(0,1),g∼N(0r,Ir)G_{0}\sim\mathcal{N}(0,1),{\bm{g}}\sim\mathcal{N}(\bm{0}_{r},{\bm{I}}_{r}), only now (H.10) is modified to:

This shows (4.4b) after applying Gaussian integration by parts and expressing it in matrix form; see Lemma C.3.

H.2 Orthogonal means and equal-energy

Here, we use Theorem 4.2 to prove that, in contrast to the GMM, in the MLM under orthogonal and equal-energy means: LS outperforms the averaging classifier for large enough sample sizes. Assuming orthogonal means of equal energy μ\mu:

Thus, similar to (4.3) and with the same notation,

H.3 Proof of Proposition 4.3

In (E.7) and (H.24), we showed the following limits for orthogonal means of equal-energy μ>0\mu>0:

We compare the expression on the RHS in the above display by applying Lemma H.1 below with the following substitutions

such that pIk−q1k1kT≻0p{\bm{I}}_{k}-q\bm{1}_{k}\bm{1}_{k}^{T}\succ 0 and fixed αIk−β1k1kT≻0\alpha{\bm{I}}_{k}-\beta\bm{1}_{k}\bm{1}_{k}^{T}\succ 0. Then, the following statements are true.

For 0<p2<p10<p_{2}<p_{1} and any q1<p1k,q2<p2kq_{1}<\frac{p_{1}}{k},q_{2}<\frac{p_{2}}{k}, it holds that F(p2,q2)<F(p1,q1)F(p_{2},q_{2})<F(p_{1},q_{1}).

Proof Fix any p>0,q≤pkp>0,q\leq\frac{p}{k}. Denote T:=(pIk−q1k1kT)1/2{\bm{T}}:=\left(p{\bm{I}}_{k}-q\bm{1}_{k}\bm{1}_{k}^{T}\right)^{1/2} and S:=(αIk−β1k1kT)1/2\bm{S}:=\left(\alpha{\bm{I}}_{k}-\beta\bm{1}_{k}\bm{1}_{k}^{T}\right)^{1/2} for convenience. It can be checked that T:=(pIk+p−qk−pk1k1kT){\bm{T}}:=\left(\sqrt{p}{\bm{I}}_{k}+\frac{\sqrt{p-qk}-\sqrt{p}}{k}\bm{1}_{k}\bm{1}_{k}^{T}\right) and S:=(αIk+α−βk−αk1k1kT)\bm{S}:=\left(\sqrt{\alpha}{\bm{I}}_{k}+\frac{\sqrt{\alpha-\beta k}-\sqrt{\alpha}}{k}\bm{1}_{k}\bm{1}_{k}^{T}\right). From these, it follows directly that

Next, we show the second statement. Using the distribution of c(h)c(\bm{h}) and symmetry we have the following chain of equalities:

where in the last line we used the rotational symmetry of the Gaussian distribution:

and the fact that g1,…,gk−1{\bm{g}}_{1},\ldots,{\bm{g}}_{k-1} are independent.

Next, we will show that the function G(⋅){\cal{G}}(\cdot) defined above is strictly decreasing in (0,∞)(0,\infty). Towards this goal, using Q′(x)=−12πe−x2/2=−ϕ(x)Q^{\prime}(x)=-\frac{1}{\sqrt{2\pi}}e^{-x^{2}/2}=-\phi(x) and using the shorthand

we may compute the derivative of G{\cal{G}} at any s>0s>0 as follows:

Next, we use Gaussian integration by parts (GIBP) to further simplify the expression in (H.26). Fix any i∈[k−1]i\in[k-1]. Then, by (GIBP):

where in the second line, we used the fact that h∼N(0,Ik)\bm{h}\sim\mathcal{N}(\bm{0},{\bm{I}}_{k}) to compute

where in the penultimate line we used the fact that ϕ′(x)=−xϕ(x)\phi^{\prime}(x)=-x\phi(x). Consider the two terms in (H.33). Clearly,

where, we have recalled (H.35) and (H.31). Using (H.37) in (H.36), we find that

From this, (H.26) and (H.27), we have shown that G{\cal{G}} is strictly decreasing in (0,∞)(0,\infty). Recalling the definition of G{\cal{G}} in (H.25), this implies that F(p,q)F(p,q) is strictly increasing in p>0p>0, as desired to complete the proof.

Appendix I Weighted LS for GMM (Proof of Theorem B.1)

To continue, consider the singular value decomposition

Since UTh\bm{U}^{T}\bm{h} is r≤kr\leq k dimensional in our asymptotic regime the term hTUαn\frac{\bm{h}^{T}\bm{U}\bm{\alpha}}{n} can be ignored. Also replacing β\beta with β/n\beta/\sqrt{n} we thus arrive at

Setting the derivative with respect to u\bm{u} to zero we arrive at

Plugging the latter into the above the AO simplifies to

To continue note that in our asymptotic regime we have

and the cross terms can be ignored so that in an asymptotic sense

To continue further we shall assume D=diag(YTω)\bm{D}=\text{diag}\left(\bm{Y}^{T}\bm{\omega}\right). Note that in this case

Deterministic Analysis of the AO. Setting the derivative of the above with respect to α0\alpha_{0} to zero we arrive at

Note that the above objective has the form

Thus setting the derivatives with respect to β\beta and τ\tau to zero, we have

Combining the latter two we conclude that πTν(βτ)=1\bm{\pi}^{T}\bm{\nu}\left(\frac{\beta}{\tau}\right)=1. Thus, η=βτ\eta=\frac{\beta}{\tau} is the solution to πTν(η)=1\bm{\pi}^{T}\bm{\nu}\left(\eta\right)=1. To calculate τ\tau and hence α0\alpha_{0} we calculate f′f^{\prime} which is equal to

Now note that at the optimal point we have

Thus the AO optimization problem reduces to

First, note that the matrix A{\bm{A}} is positive definite. This can be checked by computing the Schur complement of A{\bm{A}}:

Asymptotic predictions. First, from (I.7) the bias term converges as follows:

Adding the equations on the above displays we find that

Recognize that this coincides with the optimality condition for (I.13). Thus, the proof is complete.

The analysis of (I.13) is very similar to that of (G.1). In particular we use the following decomposition

Setting the derivative of the above with respect to α0\alpha_{0} to zero we arrive at

Note that the above objective has the form

Thus, the derivatives with respect to β\beta and τ\tau to zero we have

Combining the latter two we conclude that πTν(βτ)=1\bm{\pi}^{T}\bm{\nu}\left(\frac{\beta}{\tau}\right)=1. Thus, η=βτ\eta=\frac{\beta}{\tau} is the solution to πTν(η)=1\bm{\pi}^{T}\bm{\nu}\left(\eta\right)=1. To calculate τ\tau and hence β0\beta_{0} we calculate g′g^{\prime} which is equal to

Now note that at the optimal point we have

Thus, similar to (I.10) we can compute the minimizer of the deterministic

Finally, using (I.2) and (I.1) in (I.14) it follows that

Let us end by simplifying Q\bm{Q} to this aim

Thus, defining π~′=π⊙ν′\widetilde{\bm{\pi}}^{\prime}=\bm{\pi}\odot\bm{\nu}^{\prime} we have

Using the above and recalling π~′=π⊙ν′\widetilde{\bm{\pi}}^{\prime}=\bm{\pi}\odot\bm{\nu}^{\prime} we arrive at

Using the above the cross-correlation matrix Σw,w\bm{\Sigma}_{\bm{w},\bm{w}} is given by

Appendix J Weighted LS for MLM (Proof of Theorem B.2)

Let D:=D(n):=diag(D1,…,Dn){\bm{D}}:={\bm{D}}^{(n)}:=\text{diag}(D_{1},\ldots,D_{n}) be a diagonal matrix with non-zero diagonal entries. In particular, assume that the diagonal entries of D{\bm{D}} are distributed Di∼iidDD_{i}\stackrel{{\scriptstyle iid}}{{\sim}}D where the random variable DD may depend on the entries of the matrix of response variables Y{\bm{Y}}. Here, we focus on the following setting:

MMT=VΣ2VT{\bm{M}}{\bm{M}}^{T}={\bm{V}}\bm{\Sigma}^{2}{\bm{V}}^{T}, and g∼N(0,Ir).{\bm{g}}\sim\mathcal{N}(\bm{0},{\bm{I}}_{r}).

where D{\bm{D}} is as in (J.1). In fact, it is convenient to rewrite the above as follows:

Identifying the AO. The PO in (J.3) is very similar to (H.1). In particular, following step by step the same decomposition trick as in Section G.1.1, it can be shown that the AO corresponding to (J.3) becomes (cf. (H.6))

Note that the resulting minimization is convex in u{\bm{u}} and concave in β\beta. Also, by considering the bounded AO (such that β\beta is bounded; see [DKT19, Sec. A]), we can flip the order of min-max and optimize over u{\bm{u}} first. In particular, u{\bm{u}} minimizes the following strictly convex quadratic

Putting things together, the new objective function of (J.5) becomes

where the expectation is over g∼N(0r,Ir){\bm{g}}\sim\mathcal{N}(\bm{0}_{r},{\bm{I}}_{r}) (with some abuse of notation) and

Thus, at optimality either β=0\beta=0 or τ=α0γ\tau=\alpha_{0}\sqrt{\gamma}. In what follows, consider the solution τ=α0γ\tau=\alpha_{0}\sqrt{\gamma}. We will show that this leads to the true saddle point of D{\cal{D}}.

Rearranging (J.13) and using τ=α0γ\tau=\alpha_{0}\sqrt{\gamma} gives the following equation for η\eta:

where we have also used the RHS of (J.14). Next, we specialize these findings to the special structure of the weighting matrix D{\bm{D}} in (J.1).

Applying weighting (J.1). Assume (J.1) holds. In this case, Equation (J.14) that determines the value of η>0\eta>0 becomes

Also, in this case we can write (J.11) in the following more convenient form:

Because of (J.16), notice that π~\widetilde{\bm{\pi}} is a probability vector, i.e.

Using (J.30) and (J.31), we conclude from (J.12) the following expressions for a\bm{a} and b\bm{b}:

Finally, we show how to compute α0\alpha_{0} using (J.15). The RHS in (J.15) can be computed as

Put together, we have the following expression for α0\alpha_{0}:

Asymptotic Predictions. Writing (J.24) in vector form we find that

Further recall the matrix A′{\bm{A}}^{\prime} in (J.27).

Thus, what changes in the calculations above is in (J.18) and (J.35), where we now have instead