{'fraction_non_alphanumeric': 0.05657132778741292, 'fraction_numerical': 0.015403157059290396, 'mean_word_length': 3.9321493294671983, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 61, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Theoretical Insights Into Multiclass Classification: A High-dimensional Asymptotic View

Contemporary machine learning applications often involve classification tasks with many classes. Despite their extensive use, a precise understanding of the statistical properties and behavior of classification algorithms is still missing, especially in modern regimes where the number of classes is rather large. In this paper, we take a step in this direction by providing the first asymptotically precise analysis of linear multiclass classification. Our theoretical analysis allows us to precisely characterize how the test error varies over different training algorithms, data distributions, problem dimensions as well as number of classes, inter/intra class correlations and class priors. Specifically, our analysis reveals that the classification accuracy is highly distribution-dependent with different algorithms achieving optimal performance for different data distributions and/or training/features sizes. Unlike linear regression/binary classification, the test error in multiclass classification relies on intricate functions of the trained model (e.g., correlation between some of the trained weights) whose asymptotic behavior is difficult to characterize. This challenge is already present in simple classifiers, such as those minimizing a square loss. Our novel theoretical techniques allow us to overcome some of these challenges. The insights gained may pave the way for a precise understanding of other classification algorithms beyond those studied in this paper.


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: • 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.
• 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.
• 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.
Notation. We use [k] to denote {1, . . . , k}. We use boldface lowercase letters x, y, µ, . . . to denote vectors and boldface uppercase letters X, Y , M , . . . for matrices. We write e for the -th standard basis vector in R k . We also write I k , 0 k×k and 1 k for the k × k identity and all-zeros matrices and the k × 1 all-ones vectors. For a vector c ∈ R k we write arg max c to denote the index of its largest entry, i.e., arg max c = arg max j∈[k] c i . The superscript † denotes pseudoinverse. We use Q(x) for the tail of a standard Gaussian (Q-function). Finally, we reserve variables G 0 , G 1 , . . . , G k iid ∼ N (0, 1) to denote i.i.d. standard Gaussians.

Problem formulation
We focus on multiclass classification problems with k classes. Specifically, we assume the training data consists of n feature/label pairs {(x i , Y i )} n i=1 with x i ∈ R d representing the features and Y i ∈ {1, 2, . . . , k} the associated labels representing one of k classes. It will be convenient to also model the labels as one-hot encoded vectors y i ∈ R k representing one of k classes with one-hot encoding, i.e., y i = e Yi . Therefore, when convenient we shall use {(x i , y i )} n i=1 to represent the training data. Throughout, we shall use X = x 1 x 2 . . . x n ∈ R d×n and Y = y 1 y 2 . . . y n ∈ R k×n , to denote the matrix of features and their labels aggregated into a matrix, respectively. We shall also use Y ∈ R n to denote the -th row of Y . In our analysis we focus on training linear classifiers. Specifically, we use W = w 1 w 2 ⋯ w k T ∈ R k×d and b ∈ R k to denote the weights and biases of this linear model, respectively. The overall input-output relationship of the classifier in this case is a function that maps an input vector x ∈ R d into an output of size k via x ↦ W x + b ∈ R k , where a training algorithm is used to train the corresponding weights W ∈ R k×d and biases b ∈ R k . Next we detail the data models and training algorithms that are formally studied in this paper. We end this section by discussing how the test error can be calculated for the different data models.

Data Models
In our theoretical analysis we assume the training data {(x i , Y i )} n i=1 (alternatively {(x i , y i )} n i=1 ) are generated i.i.d. according to (x, Y )/(x, y). We consider two models for the distribution of (x, y) which we detail next. In both models we shall use mean/regressor vectors {µ } k =1 ∈ R d and aggregate them into columns of a matrix of the form M ∶= µ 1 µ 2 . . . µ k ∈ R d×k .
In the first model, these vectors represent the mean of the features conditioned on the class, i.e., µ = E x Y = , whereas in the second model these vectors can be viewed as regressor coefficients. We shall refer to {µ } k =1 /M as "mean" vectors/matrix in both models. We denote the Grammian matrix of means as Σ µ,µ = M T M . Furthermore, we shall use µ ∶= µ 2 to denote the norm of the mean vector µ . Gaussian Mixture Model (GMM). In this model each example (x, Y ) belongs to class ∈ [k] with probability π , i.e., P{Y = } = π . We let π = π 1 π 2 . . . π k T ∈ R k denote the vector of priors which of course obeys π ≥ 0 and 1 T π = 1. Also, we model the class conditional density of an example in class with an isotropic Gaussian centered at a mean vector µ . In particular, we say that a data point (x, Y ) (or its one-hot encoded representation (x, y)) follows the GMM model when (2.1) We note that for a training set summarized by the feature and label matrices X and Y with columns generated i.i.d. according to the above distribution we have: Multinomial Logit Model (MLM). In this model we assume that feature vectors x are distributed i.i.d. N (0, 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 ) (or its one-hot encoded representation (x, y)) follows the multinomial logit model when

Classification algorithms
As mentioned earlier, in this paper we focus on training linear classifiers of the form x ↦ W x + b with W ∈ R k×d denoting the weights and b ∈ R k the offset values.
Least-squares (LS). In this approach we train a linear classifier x ↦ W x + b via a least-squares fit to the training data: Class averaging (Avg). This approach uses the following weight and offset values Let n be the number of training data from class then, equivalently,ŵ = n n 1 n ∑ n i∶ Yi= x i andb = n n . Therefore, this classifier picks weights according to the empirical mean of features of each class multiplied by the relative frequency of that class and the offset value as the fraction of data points from that class. We note that this algorithm has the same classification performance as the outcome of the ridge-regularized least-squares with infinite regularization.
Weighted Least-squares (WLS). This is a variation of the Least-squares approach where we fit a weighted least squares loss of the form Here, D ∈ R n×n is a diagonal matrix with the ith diagonal entry equal to D ii = ω when the i-th data point is from class (i.e. Y i = ) and ω ≥ 0, ∈ [k] denote the weights. Aggregating the weights into a vector of the form ω = ω 1 ω 2 . . . ω k T ∈ R k we can rewrite D in the form In this approach the loss associated to data points to class is weighted by a factor ω 2 . For instance, if the class priors are known, a natural choice might be ω = 1 √ π . Such a weighted approach allows the classification algorithm to focus on rare/minority classes which are not well represented in the training data.
Cross-entropy (CE). In this approach the best weight/offset values are determined by fitting a cross . Theoretical analysis for CE is substantially more involved and we defer it to future work. Nevertheless, we compare with this classifier in our numerical simulations.

Class-wise and total test classification error
Let W ,b denote the parameters of a trained classifier. Now consider a fresh data sample (x, Y ) generated according to the same distribution as the training data. Once, we have learned the parameters W ,b of the classifier, the classŶ predicted by the classifier is made by a winner takes it all strategy, as follows, Y = arg max j∈[k] ⟨ŵ j , x⟩ +b j . Therefore, the classification error condition on the the true label being c, which we shall refer to as the class-wise test error, is equal to Correspondingly, the total classification error is given by For both the GMM and MLM, the classification error depends on the vector of interceptsb ∈ R k and the following key "correlation" matrices: GMM. In model (2.1), the test error probability is explicitly given by and Y is independent of g with probability mass function P{Y = } = π , ∈ [k].

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 and Σ w,µ , the offset valuesb and the the class conditional means. For instance, as we show in the supplementary for GMM the class-wise errors are given by where z is a Gaussian random vector distributed as N (0, Similarly, based on (2.7) the total classification error in GMM is equal to As also detailed in the supplementary, the class-wise/total test errors for MLM similarly depends on quantities of the form P{Az ≥ t} with z a standard Gaussian random vector, A and t depending only on correlation matrices, conditional means and classifier offset-values; see Section D.3. There are a variety of algorithmic approaches to calculate P{Az ≥ t} once A and t are known based on Monte Carlo methods. Analytic bounds on this quantity have also been studied in the literature, e.g., [HH03, SL80]; see more details in Section D.

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 k.
Assumption 1 We focus on a double asymptotic regime where n, d → ∞ at a fixed ratio γ = d n > 0.
For the (weighted) least-squares classifier, we focus here in the overdetermined regime γ < 1. However, our approach is also directly applicable to regularized (or min-norm) LS/WLS in the overparameterized regime γ > 1.
For a sequence of random variables X n,d that converges in probability to some constant c in the limit above, we simply write X n,d P → c. For a random vector/matrix v n,d /V n,d and a deterministic vector/matrix c/C, the expressions v n,d P → c and V n,d P → C are to be understood entry-wise.

Results for Gaussian Mixture Model
In this section we discuss the asymptotics of the intercepts/correlation matrices for the averaging and the LS classifiers for the GMM. The derived formulas can be directly plugged in (2.5) and (2.7) to obtain asymptotics for the total and class-wise test error, respectively. We end this section by also characterizing the Bayes optimal estimator in this model when priors are balanced π = 1 k, ∈ [k]. Additional results on the performance of Weighted LS are deferred to the appendix.

Class averaging classifier
Proposition 3.1 Consider data generated according to GMM in an asymptotic regime with any γ > 0. For the averaging estimator discussed in Section 2.2, the following high-dimensional limits hold (3.1b) 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.
Orthogonal means, equal priors and equal energy. Next we focus on a special case with orthogonal means ⟨µ i , µ j ⟩ = 0, i ≠ j ∈ [k] of equal energy µ 2 ∶= µ ı 2 2 and of equal priors π i = π = 1 k for i ∈ [k]. In this case, the class-wise miss-classification error converges to P e c P → 1 − P{S 1 2 c z > t}, where S c = π(πµ 2 + γσ 2 )(I k−1 + 1 k−1 1 T k−1 ) and t = −πµ 2 1 k−1 . Defining after some algebraic manipulations the total classification error of the averaging estimator in this case is given by

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 Σ a diagonal positive-definite matrix and V an orthonormal matrix obeying V T V = I r .
Theorem 3.2 Consider data generated according to GMM in an asymptotic regime with γ < 1. In addition to (3.2), define the following two positive (semi)-definite matrices: P ∶= diag(π) − ππ T ⪰ 0 k×k and ∆ ∶= σ 2 I r + ΣV T P V Σ ≻ 0 r×r . Then, for the least-squares linear classifier W ,b the following limits are true asymptoticallyb (3.3b) 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ŵ i , i ∈ [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.

Bayes estimator for the balanced Gaussian Mixture Model
To check how far the above algorithms are from the lowest misclassification error achievable by any algorithm in this section, we consider a Bayesian setting with Gaussian mean vectors and we derive the Bayes-optimal risk for the case of equal priors. Recall that the Bayes estimatorŶ = arg max ∈[k] P{Y = X, Y , x} minimizes the risk Under Gaussian prior, the means are asymptotically orthogonal and equal-energy. As shown earlier, in this setting, P e,Avg P → P G 0 + max ∈[k−1] G ≥ u Avg . But, u Avg = u Bayes . Thus, the averaging method is (asymptotically) Bayes optimal for equal-norm, orthogonal means and balanced classes. An analogous result was derived in [LM19, MKLZ20], but only for binary classification.

Results for Multinomial Logit Model
In this section we discuss the asymptotics of the intercepts/correlation matrices for MLM. We present results for arbitrary mean-vectors as well as special cases where the means are mutually orthogonal. Recall the eigenvalue decomposition of the Grammian Σ µ,µ = V Σ 2 V T in (3.2). In order to state our results, it is convenient to introduce the following probability vectors in R k and R k 2 : Note that π and Π are the first and second moments of the soft-max mapping of V Σg ∼ N (0, Σ µ,µ ). In fact, for the MLM in (2.2) it holds that since M T x is distributed as V Σg. Thus, π 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).

Class averaging classifier
Proposition 4.1 Consider data generated according to MLM in an asymptotic regime with any γ > 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,µ Σ † µ,µ Σ T w,µ P → γ ⋅ diag(π) the test error obtains the following explicit form:

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.
Theorem 4.2 Consider data generated according to MLM in an asymptotic regime with 0 < γ < 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 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(π) substituted by the matrix above.

Numerical Results
This section validates our theory via numerical experiments and provides further insights on multiclass classification. See also Section A for more extensive experiments. We study the class-wise/total test misclassification error in both GMM and MLM for different sample sizes, number of classes and class priors. In line with Section 2.2 we consider four algorithms: (i) Averaging (Avg), (ii) LS, (iii) Weighted LS (WLS) with the ith class weighted by ω 2 = 1 π , (iv) Cross-Entropy (CE).
Figures 1 and 2 focus on GMM with k = 9 classes, d = 300 and µ i 2 2 = 15. To model different class prior probabilities, we use the distribution π 1 = π 2 = π 3 = 0.5, π 4 = 0.5, π 5 = 0.5, π 6 = 0.25, π 7 = 0.25, π 8 = 0.25, π 9 = 1 21. We consider three scenarios: (a) orthogonal means, equal prior (π i = 1 9); (b) orthogonal means, different prior; (c) correlated means with pairwise correlation coefficient equal to 0.5 (i.e., ⟨µ i , µ j ⟩ ( µ i 2 µ j 2 ) = 0.5 for i ≠ j) and different priors as discussed above. Figure 1 shows the test miss-classification errors as a function of γ ∶= d n. In all scenarios our theoretical predictions are a near perfect match to the empirical performance. In scenario (a), class-wise averaging achieves the lowest error as predicted by Proposition 3.4. However, in scenario (b) where the means have different norms the averaging method has higher misclassification error compared with CE, LS and WLS for large sample sizes (small γ). We note that both LS and WLS achieve lower errors compared with CE as the sample size grows. Scenario (c) is similar to (b). However, due to class correlations, the errors are uniformly higher. Figure 2 shows the corresponding class-wise miss-classification errors for the smallest γ in Figure 1 (γ = 0.117). In scenario (a), errors are equal which is expected given the equal class priors. In scenarios (b) and (c) however, due to different priors, large classes 7,8,9 achieve best accuracy. The performance difference is most visible for the averaging approach. LS mitigates this issue to some extent, while WLS creates the flattest class-wise errors suggesting that it can reduce the miss-classification error on small/minority classes. Figure 3 focuses on orthogonal classes with varying number of classes k where µ i 2 2 = 15 and d ∈ {50, 100, 200} with kd n = kγ fixed at kγ = 20 11. It plots the ratio of the empirical error probability and our theoretical prediction as k grows until k = d. Two observations are worth mentioning here. (1) The accuracy of our predictions noticeably improves as the problem dimension d, n grow as expected given the asymptotic nature of our analysis. Interestingly, the convergence appears to be noticeably faster (as a function of d) for the LS rather than the Averaging classifier. (2) Our theoretical results formally require that k is fixed while d (and n) grow large. Yet, the presented experimental results suggest that they might also hold for large k under the shown scaling. This is a fascinating research question that we believe is worth investigating further. Figure 4 provides experiments on MLM with k = 9 orthogonal classes. Unlike GMM, CE achieves the best performance in MLM. In Figure 4 (a), classes have same norms µ i 2 = 10, while in Figure 4 (b) we have quadrupled the norms of classes 7,8,9 and doubled the norms of classes 4,5,6. This disparity between the norms seems to help improve the CE accuracy, but hurt LS/averaging accuracy for small γ. Finally, Figure  4 (c) shows the class-wise probability of error associated with (b) for γ = 0.117 and demonstrates that LS outperforms averaging.     6 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,µ , and, Σ w,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 ↦ W x + b via a least-squares fit to the training data: Notice that the objective function above is separable. That is, Hence, for each ∈ [k], This decomposition is convenient for analysis as it is easier to compute the statistical properties of the simple single-output LS in (6.2) compared to the multi-output objective in (6.1). Indeed, as we show, in Step III, this simplification will eventually allow us to compute the high-dimensional behavior of the following key quantities for all ∈ [k]: (i) the interceptb , (ii) the mean-correlations ⟨ŵ , µ c ⟩, c ∈ [k], (iii) the norm ŵ 2 .
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 (x) + ∶= max(0, x) and g ∈ R n and h ∈ R d are two independent Gaussian random vectors distributed as N (0, I n ) and N (0, I d ).
Step III: Simplification of the AO and computing Σ w,µ and b. In this step we carry out a series of intricate calculations to further simplify (6.3) and characterize its various asymptotic properties. At a high-level, we follow the principled machinery introduced in [TOH15, TAH18], organizing our analysis in three intermediate steps: (a) Scalarization; (b) Convergence analysis; and (c) Deterministic analysis. We note that each one of these intermediate steps for the multiclass setting is more involved than in previously considered regression and binary classification settings. The detailed derivations are deferred to the Appendix G.1. At the end of this analysis step, we have computed the high-dimensional behavior of the interceptŝ b , ∈ [k], the mean-correlations ⟨ŵ , µ c ⟩, , c ∈ [k], the norms ŵ 2 , ∈ [k] and the LS training loss X Tŵ +b 1 n − Y 2 . In particular, for GMM these calculations allow us to conclude the following limits for all ∈ [k]: (6.5) where ∆ ∶= σ 2 I r + ΣV T P V Σ ≻ 0 r×r and P ∶= diag(π) − ππ T . Expressing (6.4) in matrix form leads to (3.3a) in Theorem 3.2. Thus, it remains to prove (3.3b), i.e., to determine the high-dimensional limit of Σ w,w . Note that (6.5) already determines the diagonal entries of Σ w,w . However, thus far, our analysis treats the optimization of each classifierŵ , ∈ [k] independently and provides no information for the cross-correlation ⟨ŵ ,ŵ c ⟩, ≠ c ∈ [k] Step IV: Computing Σ w,w and capturing cross-correlations. The final and most involved part of our analysis is characterizing the asymptotic behavior of Σ w,w . To see why this is particularly challenging note that the reduction from (6.1) to (6.2) "breaks" the dependence of allŵ 1 ,ŵ 2 , . . . ,ŵ k on the same feature matrix X. Capturing this dependence is crucial in determining the "cross-correlations" ⟨ŵ ,ŵ c ⟩, ≠ c. As noted in Section 2.3 the matrix Σ w,w is needed to calculate the class-wise and total miss-classification errors. Unfortunately, the CGMT is not directly applicable to the multi-output LS optimization in (6.1). Our idea to circumvent this challenge builds on the following simple observation: the vectorŵ ,c =ŵ +ŵ c is itself the solution to another simple single-output LS problem.
, letŵ ,ŵ c be the and c-th row of W which is the solution to the multi-output least-squares minimization (6.1). Denoteŵ ,c ∶=ŵ +ŵ c . Then,ŵ ,c is a minimizer in the following single-output least-squares problem: Thanks to Lemma 6.1, we can use the CGMT to characterize the limiting behavior of ŵ +ŵ c 2 . These calculations are similar to (but, in certain cases, such as for weighted least-squares, more involved than) those in Steps II and III above. Now note that an asymptotic characterization of ŵ +ŵ c 2 immediately yields the asymptotic characterization of ⟨ŵ ,ŵ c ⟩ as and ŵ 2 , ŵ c 2 are already computed in Step IV (cf. (6.5)). For the GMM, the analysis in this step allow us to calculate the asymptotic behavior of Σ w,w as promised in (3.3b) in Theorem 3.2:

Background on the CGMT
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: where: G ∈ R n×d , g ∈ R n , h ∈ R d , they all have entries iid Gaussian; the sets S w ⊂ R d and S u ⊂ R n are compact; and, ψ ∶ R d × R n → R. For these two processes, define the following (random) min-max optimization programs, which are refered to as the primary optimization (PO) problem and the auxiliary optimization AO: If the sets S w and S u are convex and bounded, and ψ is continuous convex-concave on S w × S u , then, for any ν ∈ R and t > 0, it holds [TOH15, Thm. 3]: (6.9) In words, concentration of the optimal cost of the AO problem around q * implies concentration of the optimal cost of the corresponding PO problem around the same value q * . Asymptotically, if we can show that Moreover, starting from (6.9) and under appropriate strict convexity conditions, the CGMT shows that concentration of the optimal solution of the AO problem implies concentration of the optimal solution of the PO around the same value. For example, if minimizers of (6.8b) satisfy w φ (g, h) 2 P → α * for some α * > 0, then, the same holds true for the minimizers of (6.8a): Thus, one can analyze the AO to infer corresponding properties of the PO, the premise being of course that the former is simpler to handle than the latter.
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 .

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,µ requires new ideas. The previously mentioned study of "extreme multiclass classification" in which the number of classes k is very large is another fascinating direction.
[   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 γ ⋆ 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 k 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 k. Following this intuition, as k increases, γ ⋆ shifts smaller due to larger sample requirement. Similarly energy µ directly controls the noise level of the problem, i.e., larger µ results in a larger signal-to-noise ratio. Thus, as we increase µ, γ ⋆ increases as well because same test accuracy can be achieved with smaller sample size.
Second, Figure 6 provides further experiments on the class-wise probabilities of the MLM model with k = 9 classes for γ = 0.1. Classes 1,2,3 have norms µ i 2 = 15, while we have quadrupled the norms of classes 7,8,9 and doubled the norms of classes 4,5,6. In scenario (a) the means are orthogonal and in scenario (b) the means are highly correlated. The WLS shown corresponds to the following choice of weights: , where π is the th entry of the vector π in (4.1). The theoretical predictions for the class-wise error probabilities are computed using formula (D.7). As was the case for the GMM in Figure 2, we see that WLS creates the flattest class-wise errors.
Finally, in Figure 7 we investigate the following question: To what extent pairwise class correlations are necessary for performance prediction? Specifically, we consider a GMM setup with k = 5 classes and orthogonal means under three scenarios: (a) π 1 = π 2 = π 3 = π 4 = π 5 , 4 µ 1 2 = 4 µ 2 2 = 2 µ 3 2 = 2 µ 4 2 = The solid lines are exact performance predictions based on our theory for averaging and least-squares estimators. The dashed lines are the theoretical upper bounds, which do not require the knowledge of cross-correlations between the classes (i.e., off-diagonal entries of Σ w,w are unknown). These bounds are calculated by applying a union bound to the class-wise probabilities P e c , c ∈ [k] in (D.1) and further appropriately bounding the off-diagonal entries of Σ w,w in terms of the self-correlations of the classes, i.e., only the diagonal entries of Σ w,w . Please see Section D.4.3 for details. Overall, the bounds shown only depend onb, Σ w,µ and diag(Σ w,w ), which can all be obtained by studying the properties of isolated least-squares on individual classes without understanding their pairwise relations. While this suggests a simpler method to calculate theoretical bounds, there is a visible gap between such upper bounds and exact bounds and this gap is particularly more visible in the third scenario (c), where the bound becomes vacuous for LS. The gap remains visible in scenarios (a) and (b). At least, in these two cases comparing the bounds for averaging and LS to each other reveals the transition in performance gain between the two estimators. However, the cross-point of the curves does not coincide with the true one. This empirical study emphasizes the fact that pairwise correlations are indeed critical for exact asymptotic analysis and naive approaches cannot reproduce, in general, the results of our sharp analysis.

B.1 WLS for GMM
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 π toπ defined in the theorem (modulo the extra additive term in the cross correlation matrix Σ w,w ). This shows that weighted LS has similar performance to an un-weighted LS applied to a model with different class priors π. 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.
Theorem B.2 Consider data generated according to MLM and γ < 1. Consider a weighted LS classifier with weights D = diag(ω 1 , . . . , ω k ) and let η be the unique solution to ∑ k =1 π ω 2 ω 2 +η = γ. Also define vector ν ∈ R k with entries given by ν = 1 γ ω 2 ω 2 +η and matrix where v ∈ R k is a random vector with entries V = e e T V Σg ∑ ′ ∈[k] e e ′ V Σg for g ∼ N (0, I r ). Then, for the WLS linear classifier W ,b the following asymptotic limits hold The corresponding formula for the asymptotic limit of the cross-correlation matrix Σ w,w is given in (J.38) in Section J.
The term E ν T v gg T in (B.2) can be computed using Monte Carlo sampling. It is also possible to slightly simplify the calcuations involved using Gaussian integration by parts as shown in Lemma C.3. As mentioned, the formula that predicts Σ w,w is given in (J.38). While somewhat more complicated than formulae (B.3), the expression that we provide is also explicit. Numerical simulations shown in Figure 6 in Section A validate the accuracy of the theoretical predictions of the theorem.

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

C.2 Gaussian integration by parts
We say that a function F ∶ R m → R is of moderate growth if for each c > 0, The following result is a direct application of Gaussian integration by parts; for instance, see [FR13, Prop. 8.29].
Lemma C.2 (Gaussian integration by parts (GIP) ) Let g ∼ N (0, I r ) and function f ∶ R r → R such that f and all its first and second order partial derivatives are of moderate growth. Then, the following statements are true: The following is a corollary of Lemma C.2 applied to the soft-max function. (C.1) Further recall the notation of π and Π in (4.1). The following statements are true: e j Π ij , and in matrix form: Then,

C.3 Block matrix inversion
Lemma C.4 (Block matrix inversion) Let T = A b b T δ be an invertible block matrix. Then D Calculating and bounding the missclassification error D.1 Proof of (2.5) and (2.6) GMM. Starting from (2.4) and using the fact that or, in matrix-form: Recall that Σ w,µ ∶= W M and note that W z is a zero-mean Gaussian vector with covariance matrix σ 2 W W T = σ 2 Σ w,w in order to conclude with the desired formula in (2.5). MLM. Recall from (2.2) that x ∼ N (0, I d ) and Y is distributed such that P{Y = x} = e ⟨µ ,x⟩ ∑ j∈[k] e ⟨µj ,x⟩ . Let g = W x and h = M T x. In this notation, (2.4) becomes P e = P arg max (g +b ) ≠ Y , To complete the proof of (2.6), it is easy to check that g h defined above is jointly Gaussian with zero-mean and covariance matrix Σ w,w Σ w,µ Σ T w,µ Σ µ,µ .

D.2 Class-wise and total miss-classification error for GMM
The class-wise miss-classification error for GMM is given by where we used that x = µ c + z. Let S c ∈ R (k−1)×(k−1) be a symmetric matrix and t c ∈ R k−1 a vector with entries: Then, we can rewrite (D.2) as where the inequality in the rightmost expression applies entry-wise. Further, by using the law of total probability we have π c P e c = k c=1 π c 1 − P S 1 2 c z ≥ t c .

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): and P{Y (h) = } = e h ∑ j∈[k] e hj , ∈ [k]. Using Gaussian decomposition we can write g =g + Σ w,µ Σ † µ,µ h whereg ∼ N (0 k , Σ w,w − Σ w,µ Σ † µ,µ Σ w,µ ). Using this, we have Further recalling the decomposition P e = ∑ c P Ŷ ≠ Y Y = c P {Y = c} and noting that we can see from (D.5) that the class-wise error probabilities can be calculated as follows: where π c is the c th entry of the vector π in (4.1) and S c , t c (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, I k−1 ), Here, A ⪰ 0 ∈ R (k−1)×(k−1) , t ∈ R k−1 are appropriate coefficient matrices (see (D.4) and (D.5)) andã i denotes the ith row of the matrix A 1 2 . For example, (D.8) maps to (D.4) for A ← S c and t ← (−t c ). 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).

D.4.1 A special case: Rank-one update of Identity
First, we discuss the case where the coefficient matrix A and vector t in (D.8) take the special form A ∝ I +11 T and t ∝ 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.
Lemma D.1 Let A = I k + 1 k 1 T k and g ∼ N (0, A). Then, for any t ∈ R, Proof For each i ∈ [n], we can decompose g i = G 0 + G i , where G 0 , G 1 , . . . , G k are iid standard normals. Indeed, it can be readily checked from this that E[g 2 i ] = 2 and E[g i g j ] = E[G 2 0 ] = 1, i ≠ j, which is consistent with g ∼ N (0, A). Thus, we can write which completes the proof.

D.4.2 Slepian's bound
When the matrix 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] A ij ≥ 0. To begin, note that A ≥ (diag(A) − aI) +a11 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: Union bound without knowledge of cross-correlations ⟨w i , w j ⟩, i ≠ j. It is worth noting that the upper bound in (D.11) requires knowledge of the cross-correlations ⟨ŵ j ,ŵ c ⟩, j ≠ c, i.e., of the off-diagonal entries of Σ w,w . Thankfully, our analysis allows predicting these values. For comparison, we ask wether it is possible to further upper bound the class-wise error probability if only the diagonal entries of Σ w,w (i.e., the norms ŵ j 2 , j ∈ [k]) were known. A simple answer to this questions is as follows. Observe that Unfortunately, this bound becomes non-trivial for the class-wise probability of error only if Intuitively, this assumes a regime wherethe weight vectorŵ c corresponding to class c aligns better with the corresponding mean vector µ c than the rest of the weight vectorsŵ j , j ≠ c. This emphasizes the important role of the cross-correlation matrix Σ w,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.
Oracle lower bound. For completeness, we briefly discuss an oracle lower bound for the class-wise probability of error in GMM. Specifically, assume that the means µ i , i ∈ [n] are known. Then the problem of classifying a new sample x is a k-ary hypothesis testing problem with Gaussian conditionals. Denote P genie,Bayes the Bayes error of this hypothesis testing problem. Clearly P genie,Bayes is a lower bound on the error of any classifier that is trained on data. The paper [WBW + 16] further lower bounds P genie,Bayes in terms of the Bayesian probability of errors between every two classes as follows: P e ≥ P genie,Bayes ≥ 2 k i≠j π i P genie,ij (D.12) where P genie,ij is the Bayesian error between classes i and j with priors πi πi+πj and πj πi+πj . For the last equality we have used the well-known formula for the Bayesian probability of binary Gaussian hypothesis testing. In the case of equal-priors the genie lower bound above simplifies to Note that, in contrast to (D.13), our analysis allows for precise evaluations of the missclassification error P e . The first statement (3.1a) follows directly from the fact that 1 n 1 T Y i = ni n P → π i . For the next two statements note thatŵ

E The Class-averaging estimator
where in the last line we used orthogonality of the rows Y j of the matrix Y : To conclude simply use the facts that for all i ∈ [k]: (iii) .

(E.3)
We will prove the following three statements: These lead to (4.2) using Lemma C.3. Therefore, in what follows, we prove (E.4) For the interceptsb , ∈ [k] it holds that where h ∼ N (0 k , Σ µµ ). To deduce the first statement in (E.4a), note that h with U ∈ R d×r , Σ ∈ R r×r , and V ∈ R k×r where r =rank(M ) ≤ k. Decompose X ∈ R d×n as X = U U T X + P ⊥ X with P ⊥ = I d − U U T . With this notation we compute Here, we have recognized that for every i ∈ [n] ∶ U T x i ∼ g, and also, conditioned on . This shows the second statement in (E.4b) when expressed in matrix form.
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.
F On the Bayes risk of GMM: Proof of Proposition 3.4 Without loss of generality in this proof we assume σ = 1. The general result follows by simply replacing (µ, σ) with ( µ σ , 1) and using the proof for σ = 1. Recall that the feature vectors x 1 , . . . , x n of the training data set are given by: where the matrix of means M ∈ R d×k has iid Gaussian entries with variance µ 2 d, z i iid ∼ N (0, I d ) and y i iid ∼ Unif (e 1 , . . . , e k ) with e j denoting the j th canonical vector in R k . By definition here the Bayes estimator is the maximum-likelihood (ML) estimator. By applying the law of total probability and by successive application of the Bayes rule we have the following chain of reformulations of the ML: To arrive in (F.1) we used that P (y = e j M , X, Y ) = π, ∀j ∈ [k] and P (x n+1 y = e j , M , X, Y ) = P (x n+1 y = e j , M ). Also, (F.2) follows by recognizing that P (x n+1 X, Y ) > 0 is independent of the variable of integration M and of the optimization variable j. For the same reasons, in (F.3) we have ignored the normalizing term P (X Y ).
where ∝ hides constant positive terms. Moreover, the posterior probability of the mean matrix given the training data is given by By completing the squares and invoking a gaussian integral it can be shown that Using this in (F.5) we have that where ξ(n c ) ∶= d n µ 2 + nc n , c ∈ [k] and We conclude that log (I(j)) − 1 log (I(j)) + 1 Next, we evaluate the objective in (F.7) in the asymptotic limit n, d → ∞, n d = γ. First, since n c n P → π, note that I(j) − I( ) For each one of the four terms in (F.8), we have the following by the CLT: where in the last line we used the fact that 1 → N (0, γ). Therefore, in the asymptotic limit, the Bayes estimator is the solution to: πµ 2 ⟨y n+1 , e j ⟩ + π (πµ 2 + γ) G j , G 1 , . . . , G k iid ∼ N (0, 1). (F.9) As such, the probability of error is The LS classifier solves: Identifying the AO. To continue further note that by duality we have Note that the above is jointly convex in (w , b ) and concave in s and the Gaussian matrix Z is independent of everything else. Thus, the objective is in the form of (6.7a) and so we consider the corresponding Auxiliary Optimization (AO) problem: where g ∈ R n and h ∈ R d are independent Gaussian random vectors with i.i.d. N (0, 1) entries. Maximizing over the direction of s and setting its norm β = s 2 we arrive at To continue, consider the singular value decomposition with U ∈ R d×r , Σ ∈ R r×r , and V ∈ R k×r where r =rank(M ) ≤ k. We further decompose w in its projections on the orthogonal columns u 1 , . . . , u r of U : We also define α = α 1 α 2 . . . α k T . In this notation, we havē where in the second line we denote h ⊥ the projection of h onto the complement subspace of the span of u 1 , . . . , u r and we recalled that w ⊥ 2 = 1 and ⟨w ⊥ , At this point, observe that we have reduced the AO to an optimization problem over only r + 2 scalar variables.
Using the law of large numbers, the fact that h ⊥ 2 concentrates around √ d − r and (d − r) n P → γ, as well as the limit calculation above, it is not hard to see that for fixed α 0 , b and α = [α 1 , . . . , α r ] T ∈ R r , the objective function in (G.5) converges to the following: We will show in the next paragraph that the argument inside the square-root in (G.6) is a convex quadratic over (α 0 , α, b ) (see (G.19)). Thus, the function D (α 0 , α, b ) is jointly convex. Using uniform convergence of convex functions over compact sets [AG82, Cor.. II.1], we arrive at Deterministic Analysis. Here, we analyze the deterministic scalar minimization on the RHS of (G.7). Define A ∶= σ 2 I r + ΣV T diag(π)V Σ ΣV T π π T V Σ 1 and c = ΣV T e 1 , (G.8) and observe that we can write First, note that the matrix A is positive definite. This can be checked by computing the Schur complement of A: Positive definiteness above holds because P ∶= diag(π) − ππ T ⪰ 0 k×k . Thus the term under the square-root in (G.9) is a strictly convex quadratic. Thus, D is jointly convex in its arguments.
To simplify the RHS of (G.7) we proceed by minimizing .
(G.11) Thus, the minimum value attained is Using the above, (G.7) reduces tō Setting the derivative with respect to α 0 to zero we arrive at Thus,α Plugging the latter into (G.12) we arrive at Asymptotic predictions. First, from (G.11) the bias term converges as follows: Recall from (G.4) that α = U T w . Thus, the correlations ⟨µ i , w ⟩, i ∈ [k] converge as follows: (G.14) Here, convergence applies element-wise to the entries of the involved random vectors. Moreover, from the analysis above we can predict the limit of the norm w 2 . For this, note that w 2 2 =α 2 0 +α Tα . Thus, In the previous section we used the CGMT to predict the biasb , the correlations ⟨µ i ,ŵ ⟩, i[k] and the norm ŵ 2 for all ∈ [k] members of the multi-output classifier. Here, we show how to compute the limits of the cross-correlations ⟨ŵ ,ŵ j ⟩, ≠ j ∈ [k].
Lemma G.1 For ≠ j ∈ [k], letŵ ŵ j be solutions to the least-squares minimization Denoteŵ ,j ∶=ŵ +ŵ j andb ,j ∶=b +b j . Then, (ŵ ,j ,b ,j ) is a minimizer in the following least-squares problem: (ŵ ,j ,b ,j ) = arg min Proof Clearly the minimization in (G.16) is convex. Thus, it suffices to prove thatŵ +ŵ j satisfies the KKT conditions. First, by optimality ofŵ , we have that Similarly, forŵ j : X Y j − X Tŵ j −b j 1 n = 0. 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.
Thanks to Lemma G.1, we can use the CGMT to characterize the limiting behavior of ŵ +ŵ j 2 . Observe that this immediately gives the limit of ⟨ŵ ,ŵ j ⟩ since 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 we have decomposed with w ⊥ ,j 2 = 1 and U T w ⊥ ,j = 0 r . Using a calculation similar to the one leading to (G.9) we can show that (G.18) converges point-wise in β 0 , β = β 1 , . . . , β r , b ,j to the following: where A is as in (G.8) and we have further defined d ,j ∶= π ΣV T e + π j ΣV T e j π + π j .

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.

H Least-squares for MLM
H. Assume that X, Y are generated from the MLM. Fix any ∈ [k]. The classifier parametersŵ ,b minimize the following objective function Identifying the AO. To continue further note that by duality we have and the optimization is jointly convex in (w , b ) and concave in s. Here, note that Y depends on the Gaussian matrix X. Thus, before applying the CGMT, we need to break this dependence as follows. Consider the singular value decomposition with U ∈ R d×r , Σ ∈ R r×r , and V ∈ R k×r where r =rank(M ) ≤ k. For every i ∈ [n], we decompose x i in its projection on the subspace spanned orthogonal columns u 1 , . . . , u r as follows: where P ⊥ = I r − U U T , and we denotẽ Recalling that x i ∼ N (0, I d ) note thatg where we used (H.3) and the SVD decomposition of M . In this notation, we can rewrite the PO as follows: From (H.4) and (H.5) notice that Y depends only onG andG is independent of X T P ⊥ . Therefore, the corresponding Auxiliary Optimization (AO) problem becomes where g ∈ R n and h ∈ R d are iid Gaussian vectors independent of everything else.
Scalarization of the AO. Maximizing over the direction of s and denoting its norm β = s 2 ≥ 0 we arrive at In the remaining, we focus in the inner minimization above. Let us denote a ∶= U T w and α 0 = P ⊥ w 2 . Notice that a ⊥ P ⊥ w and thus the orthogonal decomposition w = U a + P ⊥ w . With this observation, we can optimize over the direction of P T w in (H.7) by aligning it with −P T h. With this, the minimization in (H.7) reduces to the following min a,α0≥0,b Convergence of the AO. First, we argue on point-wise convergence of the objective function in (H.8).
Fix a, α 0 and b . From the WLLN, 1 √ n P ⊥ h 2 P → √ γ and where the expectation is over g ∼ N (0 r , I r ) (with some abuse of notation) and . (H.10) Therefore, point-wise on a, α 0 and b , the objective of the AO converges to 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, Deterministic analysis of the AO. Here, we solve the deterministic minimization problem in (H.12). Optimization over b is straightforward. By setting we now have to optimize min α0≥0,α By direct differentiation and first-order optimality, we compute the optimal values as follows: , (H.14) Asymptotic Predictions. From the analysis above, we conclude with the following limits about the solutionb ,ŵ of the PO: (H.16d) Recall the notation in (4.1). Note that E[V ] = π . Moreover, using Gaussian integration by parts Lemma C.3, it can be shown that E[V g] = ΣV T (diag(π) − Π) e . Using these and writing and in matrix form, we arrive at (4.4a).

H.1.2 Computing Σ w,w
Here, we prove (4.4b). Specifically, we compute the correlations ⟨ŵ ,ŵ c ⟩, ≠ c ∈ [k] by following the strategy of Section G.1.2. Specifically, in view of Lemma G.1 we need to study the following PO: which is minimized byŵ +ŵ c . Thus the analysis will lead us to an asymptotic formula for ŵ +ŵ c 2 . This when combined with the formulae for ŵ 2 and ŵ c 2 in (H.16d) will give the desired. 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 G 0 ∼ N (0, 1), g ∼ N (0 r , I r ), only now (H.10) is modified to: Y ,c ∼ Bern(V c + V ) and as before: With these, it can be shown that Combining this with (H.16d), we conclude that for ≠ c ∈ [k]: 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 equalenergy means: LS outperforms the averaging classifier for large enough sample sizes. Assuming orthogonal means of equal energy µ: Then, where we defined p ∶= π 1 − µ 2 (π 1 − Π 11 + Π 12 ) 2 and q ∶= (π 2 1 + Π 2 12 µ 2 k − 2µ 2 Π 12 (π 1 − Π 11 + Π 12 )). (H.23) Thus, similar to (4.3) and with the same notation, In (H.24) (as well as in (4.3)), note that the matrices multiplyingg and g have all the form of a rank one update of a (scaled) identity matrix. It turns out that we can exploit this structure to simplify the formulae for the test error even further. Importantly, this lets us directly compare P e,LS and P e,Avg of the two classifiers. These are detailed in Section H.3.
To see this note that p 2 > 0 from positive semi-definiteness of the Schur matrix in (H.23). It takes simple algebra to conclude that p 2 > 0 ⇒ γ ⋆ < 1.
such that pI k − q1 k 1 T k ≻ 0 and fixed αI k − β1 k 1 T k ≻ 0. Then, the following statements are true.
Proof Fix any p > 0, q ≤ p k . Denote T ∶= pI k − q1 k 1 T k 1 2 and S ∶= αI k − β1 k 1 T k 1 2 for convenience. It can be checked that This shows the first statement. Next, we show the second statement. Using the distribution of c(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 g 1 , . . . , g k−1 are independent. Next, we will show that the function G(⋅) defined above is strictly decreasing in (0, ∞). Towards this goal, using Q ′ (x) = − 1 √ 2π e −x 2 2 = −φ(x) and using the shorthand we may compute the derivative of G at any s > 0 as follows: where in the last line we have defined Next, we use Gaussian integration by parts (GIBP) to further simplify the expression in (H.26). Fix any i ∈ [k − 1]. Then, by (GIBP): where in the second line, we used the fact that h ∼ N (0, I k ) to compute We now compute the derivatives in (H.28). First, for any ∈ [k − 1], ≠ i, where we defined  For the second term we observe that: Moreover, using again GIBP, E[g 2 k ] = 1, E[g k H kj ] = 0, j ∈ [k] and the fact that φ ′ (x) = −xφ(x), We are now ready to put things together: The WLS estimator solves: Identifying the AO. By duality we have Note that the above is jointly convex in (u, w , b ) and concave in s. Thus, we consider the Auxiliary Optimization (AO) problem where g ∈ R n and h ∈ R d are independent Gaussian random vectors with i.i.d. N (0, 1) entries. Moreover, we carry out a change of variable s → Ds to arrive at Simplification of the AO. Maximizing over the direction of s and setting its norm β = s 2 above we arrive at To continue, consider the singular value decomposition with r ∶= rank(M ) ≤ k and define the variable α = U T w and α ⊥ = U T ⊥ w where U ⊥ is the orthogonal complement of the columns of U . With these definitions the above optimization problem reduces to Decomposing the optimization over α ⊥ in terms of its direction and norm α 0 = α ⊥ 2 we arrive at Since U T h is r ≤ k dimensional in our asymptotic regime the term h T U α n can be ignored. Also replacing β with β √ n we thus arrive at Setting the derivative with respect to u to zero we arrive at Plugging the latter into the above the AO simplifies to n To continue note that in our asymptotic regime we have and the cross terms can be ignored so that in an asymptotic sense Therefore we arrive at which can be rewritten in the form To continue further we shall assume D = diag Y T ω . Note that in this case Also, Next define We thus arrive at Deterministic Analysis of the AO. Setting the derivative of the above with respect to α 0 to zero we arrive at Note that the above objective has the form Thus setting the derivatives with respect to β and τ to zero, we have Combining the latter two we conclude that π T ν β τ = 1. Thus, η = β τ is the solution to π T ν (η) = 1. To calculate τ and hence α 0 we calculate f ′ which is equal to Now note that at the optimal point we have Thus from the above we can conclude that Thus the AO optimization problem reduces to where η is the solution to k =1 π ω 2 ω 2 + η = γ , and A(η) ∶= First, note that the matrix A is positive definite. This can be checked by computing the Schur complement of A: Positive definiteness above holds because P ∶= diag(π ⊙ ν) − (π ⊙ ν) (π ⊙ ν) T ⪰ 0 k×k . Thus the objective is a strictly convex quadratic and is jointly convex in its arguments. We proceed by minimizing the objective over (α, b ) which is equal to Thus, the minimum value attained is Thus the objective reduces to Therefore, To continue note that Thus, Asymptotic predictions. First, from (I.7) the bias term converges as follows: Thus,b Recall that α = U T w . Thus, the correlations ⟨µ i , w ⟩, i ∈ [k] converge as follows: (I.11) Here, convergence applies element-wise to the entries of the involved random vectors. Moreover, from the analysis above we can predict the limit of the norm w 2 . For this, note that w 2 2 =α 2 0 +α Tα . Thus, w 2 2 P → ζ σ 2 π ν (1 − π ν ) + π 2 ν 2 (π ⊙ ν − e ) T V Σ∆ −1 ∆ −1 − ζ σ 2 I r ΣV T (π ⊙ ν − e ) + ηζ σ 2 π ν ′ (η) + π 2 ν 2 c T A −1 A ′ A −1 c − 2π 2 ν ′ ν c T A −1 c . (I.12)

I.2 Computing Σ w,w
In the previous section we used the CGMT to predict the biasb , the correlations ⟨µ i ,ŵ ⟩, i[k] and the norm ŵ 2 for all ∈ [k] members of the multi-output classifier. Here, we show how to compute the limits of the cross-correlations ⟨ŵ ,ŵ j ⟩, ≠ j ∈ [k].
To continue note that Using the above and recallingπ ′ = π ⊙ ν ′ we arrive at Q =diag(π ′ ) Using the above the cross-correlation matrix Σ w,w is given by where D is as in (J.1). In fact, it is convenient to rewrite the above as follows: where we defined Finally, we show how to compute α 0 using (J.15). The RHS in (J.15) can be computed as where a, b are as in (J.23) and (J.24). Also, note that Put together, we have the following expression for α 0 : where a, b are as in (J.23), (J.24) and we have also defined Also, recalling that e T Σ w,µ =ŵ T U ΣV T P → a T ΣV T and using (J.23): Finally, for the magnitudes of the weight vectors, recall that ŵ ŵ 2 2 P → ∆ −1 ΣV T (diag(π) − Π) ⋅ ν ⋅ (e − π ν) Remark J.2 Consider the special case ω i = 1, i ∈ [k]. We show how the above recovers the solution for (un-weighted) LS. First, note that in this case (J.16) simply gives η = 1 γ − 1. Thus, ν = 1 k andπ = π. Also, recall that (diag(π) − Π) 1 k = 0 and 1 T v = 1. Thus, (J.25) simply gives ∆ = E[gg T ] = I r . With these, it can be readily checked that (J.28) and (J.23) simplify to the expressions in (4.4a). Similarly, A = A ′ = I r+1 and (J.26) reduces in this case to (H.1.1). For general weight coefficients, such simplifications do not seem possible and one needs to compute the matrix E ν T v gg T that appears in the definitions of ∆, A and A ′ . We note that this calculation can be somewhat simplified by applying Gaussian integration by parts similar to lemma C.3.

J.3 Computing Σ w,w
In this section, we use Lemma G.1 to compute the cross-correlations ⟨ŵ ,ŵ j ⟩, j ≠ ∈ [k]. Specifically, the analysis of (I.13) is almost identical to the analysis of (J.3) in the previous section. Specifically, without repeating all the details for brevity, it can be shown that the AO of (I.13) converges to min a,α0≥0,b ,τ >0 max β≥0 D(a, α 0 , b , τ, β) where D(a, α 0 , b , τ, β) is as in (J.8) only with Y substituted by Y ,c : 1 Note the slight abuse of notation compared to the definitions in (J.30) and (J.30). This "renaming" should not be confusing as the constant γ ⋅ η (that is different between the two definitions) cancels when computing a b = A −1 c (see (J.12)).