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

Christos Thrampoulidis Samet Oymak Mahdi Soltanolkotabi Department of
Electrical and Computer Engineering, University of California, Santa Barbara,
CADepartment of Electrical and Computer Engineering, University of California,
Riverside, CAMing Hsieh Department of Electrical Engineering, University of
Southern California, Los Angeles, CA

###### Abstract

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.

##  1 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.

Notation. We use $[k]$ to denote $\\{1,\ldots,k\\}$. We use boldface lowercase
letters $\bm{x},\bm{y},\bm{\mu},\ldots$ to denote vectors and boldface
uppercase letters ${\bm{X}},{\bm{Y}},{\bm{M}},\ldots$ for matrices. We write
$\bm{e}_{\ell}$ for the $\ell$-th standard basis vector in $\mathbb{R}^{k}$.
We also write ${\bm{I}}_{k},\bm{0}_{k\times k}$ and $\bm{1}_{k}$ for the
$k\times k$ identity and all-zeros matrices and the $k\times 1$ all-ones
vectors. For a vector $\bm{c}\in\mathbb{R}^{k}$ we write $\arg\max\bm{c}$ to
denote the index of its largest entry, i.e.,
$\arg\max\bm{c}={\arg\max}_{j\in[k]}\text{ }\bm{c}_{i}$. The superscript
$\dagger$ denotes pseudoinverse. We use $Q(x)$ for the tail of a standard
Gaussian (Q-function). Finally, we reserve variables
$G_{0},G_{1},\ldots,G_{k}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}(0,1)$ to denote i.i.d. standard Gaussians.

##  2 Problem formulation

We focus on multiclass classification problems with $k$ classes. Specifically,
we assume the training data consists of $n$ feature/label pairs
$\\{(\bm{x}_{i},Y_{i})\\}_{i=1}^{n}$ with $\bm{x}_{i}\in\mathbb{R}^{d}$
representing the features and $Y_{i}\in\\{1,2,\ldots,k\\}$ the associated
labels representing one of $k$ classes. It will be convenient to also model
the labels as one-hot encoded vectors $\bm{y}_{i}\in\mathbb{R}^{k}$
representing one of $k$ classes with one-hot encoding, i.e.,
$\bm{y}_{i}=\bm{e}_{Y_{i}}$. Therefore, when convenient we shall use
$\\{(\bm{x}_{i},\bm{y}_{i})\\}_{i=1}^{n}$ to represent the training data.
Throughout, we shall use

${\bm{X}}=\begin{bmatrix}\bm{x}_{1}&\bm{x}_{2}&\ldots&\bm{x}_{n}\end{bmatrix}\in\mathbb{R}^{d\times
n}\,\quad\text{and}\quad{\bm{Y}}=\begin{bmatrix}\bm{y}_{1}&\bm{y}_{2}&\ldots&\bm{y}_{n}\end{bmatrix}\in\mathbb{R}^{k\times
n}\,,$

to denote the matrix of features and their labels aggregated into a matrix,
respectively. We shall also use $\bm{Y}_{\ell}\in\mathbb{R}^{n}$ to denote the
$\ell$-th row of ${\bm{Y}}$. In our analysis we focus on training linear
classifiers. Specifically, we use

$\bm{W}=\begin{bmatrix}\bm{w}_{1}&\bm{w}_{2}&\cdots&\bm{w}_{k}\end{bmatrix}^{T}\in\mathbb{R}^{k\times
d}\quad\text{and}\quad\bm{b}\in\mathbb{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 $\bm{x}\in\mathbb{R}^{d}$ into an output of size $k$
via $\bm{x}\mapsto\bm{W}\bm{x}+\bm{b}\in\mathbb{R}^{k},$ where a training
algorithm is used to train the corresponding weights
$\bm{W}\in\mathbb{R}^{k\times d}$ and biases $\bm{b}\in\mathbb{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.

###  2.1 Data Models

In our theoretical analysis we assume the training data
$\\{(\bm{x}_{i},Y_{i})\\}_{i=1}^{n}$ (alternatively
$\\{(\bm{x}_{i},\bm{y}_{i})\\}_{i=1}^{n}$) are generated i.i.d. according to
$(\bm{x},Y)$/$(\bm{x},\bm{y})$. We consider two models for the distribution of
$(\bm{x},\bm{y})$ which we detail next. In both models we shall use
mean/regressor vectors $\\{\bm{\mu}_{\ell}\\}_{\ell=1}^{k}\in\mathbb{R}^{d}$
and aggregate them into columns of a matrix of the form

$\bm{M}:=\begin{bmatrix}\bm{\mu}_{1}&\bm{\mu}_{2}&\ldots&\bm{\mu}_{k}\end{bmatrix}\in\mathbb{R}^{d\times
k}.$

In the first model, these vectors represent the mean of the features
conditioned on the class, i.e.,
$\bm{\mu}_{\ell}=\operatorname{\mathbb{E}}\big{[}\bm{x}|Y=\ell\big{]}$,
whereas in the second model these vectors can be viewed as regressor
coefficients. We shall refer to $\\{\bm{\mu}_{\ell}\\}_{\ell=1}^{k}$/$\bm{M}$
as “mean” vectors/matrix in both models. We denote the Grammian matrix of
means as $\bm{\Sigma}_{\bm{\mu},\bm{\mu}}={\bm{M}}^{T}{\bm{M}}.$ Furthermore,
we shall use $\mu_{\ell}:=\left\|\bm{\mu}_{\ell}\right\|_{\ell_{2}}$ to denote
the norm of the mean vector $\bm{\mu}_{\ell}$.

Gaussian Mixture Model (GMM). In this model each example $(\bm{x},Y)$ belongs
to class $\ell\in[k]$ with probability $\pi_{\ell}$, i.e.,
$\operatorname{\mathbb{P}}\\{Y=\ell\\}=\pi_{\ell}$. We let
$\bm{\pi}=\begin{bmatrix}\pi_{1}&\pi_{2}&\ldots&\pi_{k}\end{bmatrix}^{T}\in\mathbb{R}^{k}$
denote the vector of priors which of course obeys $\bm{\pi}\geq\bm{0}$ and
$\bm{1}^{T}\bm{\pi}=1$. Also, we model the class conditional density of an
example in class $\ell$ with an isotropic Gaussian centered at a mean vector
$\bm{\mu}_{\ell}$. In particular, we say that a data point $(\bm{x},Y)$ (or
its one-hot encoded representation $(\bm{x},\bm{y})$) follows the GMM model
when

$\displaystyle\operatorname{\mathbb{P}}\\{Y=\ell\\}=\pi_{\ell}\qquad\text{and}\qquad\bm{x}=\bm{\mu}_{Y}+{\bm{z}},~{}{\bm{z}}\sim\mathcal{N}(\bm{0},\sigma^{2}{\bm{I}}_{d}).$
(2.1)

We note that for a training set summarized by the feature and label matrices
${\bm{X}}$ and ${\bm{Y}}$ with columns generated i.i.d. according to the above
distribution we have: ${\bm{X}}={\bm{M}}{\bm{Y}}+\bm{Z}$ where
$\bm{Z}\in\mathbb{R}^{d\times n}$ is a Gaussian noise matrix with i.i.d.
$\mathcal{N}(0,\sigma^{2})$ entries.

Multinomial Logit Model (MLM). In this model we assume that feature vectors
$\bm{x}$ are distributed i.i.d. $\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 $(\bm{x},Y)$ (or its one-hot encoded
representation $(\bm{x},\bm{y})$) follows the multinomial logit model when

$\displaystyle\bm{x}\sim\mathcal{N}(\bm{0},{\bm{I}}_{d})\qquad\text{and}\qquad\operatorname{\mathbb{P}}\\{Y=\ell~{}|~{}\bm{x}\\}={e^{\langle\bm{\mu}_{\ell},\bm{x}\rangle}}\big{/}{\sum_{j\in[k]}e^{\langle\bm{\mu}_{j},\bm{x}\rangle}}.$
(2.2)

###  2.2 Classification algorithms

As mentioned earlier, in this paper we focus on training linear classifiers of
the form $\bm{x}\mapsto\bm{W}\bm{x}+\bm{b}$ with $\bm{W}\in\mathbb{R}^{k\times
d}$ denoting the weights and $\bm{b}\in\mathbb{R}^{k}$ the offset values.

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

$(\widehat{\bm{W}},\widehat{\bm{b}}):=\arg\min_{\bm{W},\bm{b}}\frac{1}{2n}\sum_{i=1}^{n}\left\|\bm{W}\bm{x}_{i}+\bm{b}-\bm{y}_{i}\right\|_{\ell_{2}}^{2}=\frac{1}{2n}\left\|\bm{W}\bm{X}+\bm{b}\bm{1}_{n}^{T}-\bm{Y}\right\|_{F}^{2}.$

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

$\widehat{\bm{W}}:=\frac{1}{n}{\bm{Y}}{\bm{X}}^{T}\quad\text{and}\quad\widehat{\bm{b}}:=\frac{1}{n}{\bm{Y}}\bm{1}.$

Let $n_{\ell}$ be the number of training data from class $\ell$ then,
equivalently,
$\widehat{\bm{w}}_{\ell}=\frac{n_{\ell}}{n}\left(\frac{1}{n_{\ell}}\sum_{i:\text{
}Y_{i}=\ell}^{n}\bm{x}_{i}\right)\text{ and
}\widehat{\bm{b}}_{\ell}=\frac{n_{\ell}}{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

$(\widehat{\bm{W}},\widehat{\bm{b}}):=\arg\min_{\bm{W},\bm{b}}~{}\frac{1}{2n}\left\|\left(\bm{W}\bm{X}+\bm{b}\bm{1}_{n}^{T}-\bm{Y}\right)\bm{D}\right\|_{F}^{2}.$

Here, $\bm{D}\in\mathbb{R}^{n\times n}$ is a diagonal matrix with the ith
diagonal entry equal to $D_{ii}=\omega_{\ell}$ when the i-th data point is
from class $\ell$ (i.e. $Y_{i}=\ell$) and $\omega_{\ell}\geq 0,~{}\ell\in[k]$
denote the weights. Aggregating the weights into a vector of the form
$\bm{\omega}=\begin{bmatrix}\omega_{1}&\omega_{2}&\ldots&\omega_{k}\end{bmatrix}^{T}\in\mathbb{R}^{k}$
we can rewrite $\bm{D}$ in the form

$\bm{D}=\text{diag}\left(\bm{Y}^{T}\bm{\omega}\right).$

In this approach the loss associated to data points to class $\ell$ is
weighted by a factor $\omega_{\ell}^{2}$. For instance, if the class priors
are known, a natural choice might be $\omega_{\ell}=1/\sqrt{\pi_{\ell}}$. 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 entropy loss
$(\widehat{\bm{W}},\widehat{\bm{b}}):=\arg\min_{\bm{W},\bm{b}}\frac{1}{n}\sum_{i=1}^{n}\log\Big{(}\frac{\sum_{\ell=1}^{k}e^{\langle\widehat{\bm{w}}_{\ell},\bm{x}_{i}\rangle+\bm{b}_{\ell}}}{e^{\langle\widehat{\bm{w}}_{Y_{i}},\bm{x}_{i}\rangle+\bm{b}_{Y_{i}}}}\Big{)}.$
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.

###  2.3 Class-wise and total test classification error

Let $\widehat{\bm{W}},\widehat{\bm{b}}$ denote the parameters of a trained
classifier. Now consider a fresh data sample $(\bm{x},Y)$ generated according
to the same distribution as the training data. Once, we have learned the
parameters $\widehat{\bm{W}},\widehat{\bm{b}}$ of the classifier, the class
$\widehat{Y}$ predicted by the classifier is made by a winner takes it all
strategy, as follows, $\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 $c$, which we
shall refer to as the _class-wise test error_ , is equal to

$\displaystyle\operatorname{\mathbb{P}}_{e|c}:=\operatorname{\mathbb{P}}\big{\\{}\widehat{Y}\neq
Y\big{|}Y=c\big{\\}}=\operatorname{\mathbb{P}}\big{\\{}\langle\widehat{\bm{w}}_{c},\bm{x}\rangle+\widehat{\bm{b}}_{c}\leq\max_{j\neq
c}~{}\langle\widehat{\bm{w}}_{j},\bm{x}\rangle+\widehat{\bm{b}}_{j}\big{\\}}.$
(2.3)

Correspondingly, the _total classification error_ is given by

$\displaystyle\operatorname{\mathbb{P}}_{e}:=\operatorname{\mathbb{P}}\big{\\{}\widehat{Y}\neq
Y\big{\\}}=\operatorname{\mathbb{P}}\Big{\\{}\arg\max_{j\in[k]}\\{\langle\widehat{\bm{w}}_{j},\bm{x}\rangle+\widehat{\bm{b}}_{j}\\}\neq
Y\\}\Big{\\}}=\operatorname{\mathbb{P}}\big{\\{}\langle\widehat{\bm{w}}_{Y},\bm{x}\rangle+\widehat{\bm{b}}_{Y}\leq\max_{j\neq
Y}~{}\langle\widehat{\bm{w}}_{j},\bm{x}\rangle+\widehat{\bm{b}}_{j}\big{\\}}.$
(2.4)

For both the GMM and MLM, the classification error depends on the vector of
intercepts $\widehat{\bm{b}}\in\mathbb{R}^{k}$ and the following key
“correlation" matrices:

$\bm{\Sigma}_{\bm{w},\bm{w}}:=\widehat{\bm{W}}\widehat{\bm{W}}^{T}\quad\text{and}\quad\bm{\Sigma}_{\bm{w},\bm{\mu}}:=\widehat{\bm{W}}{\bm{M}}.$

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

$\displaystyle\operatorname{\mathbb{P}}_{e}=\operatorname{\mathbb{P}}\Big{\\{}\arg\max\text{
}\left(\sigma\,{\bm{g}}+\widehat{\bm{b}}+\bm{\Sigma}_{\bm{w},\bm{\mu}}\bm{e}_{Y}\right)\neq
Y\Big{\\}},~{}~{}\text{where}~{}{\bm{g}}\sim\mathcal{N}\left(\bm{0},\bm{\Sigma}_{\bm{w},\bm{w}}\right),$
(2.5)

and $Y$ is independent of ${\bm{g}}$ with probability mass function
$\operatorname{\mathbb{P}}\\{Y=\ell\\}=\pi_{\ell},~{}~{}\ell\in[k].$

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

$\displaystyle\operatorname{\mathbb{P}}_{e}=\operatorname{\mathbb{P}}\big{\\{}\arg\max\text{
}(\,{\bm{g}}+\widehat{\bm{b}}\,)\neq
Y(\bm{h})\big{\\}}\,,\quad\text{where}~{}~{}\begin{bmatrix}{\bm{g}}\\\
\bm{h}\end{bmatrix}\sim\mathcal{N}\Big{(}\bm{0},\begin{bmatrix}\bm{\Sigma}_{\bm{w},\bm{w}}&\bm{\Sigma}_{\bm{w},\bm{\mu}}\\\
\bm{\Sigma}_{\bm{w},\bm{\mu}}^{T}&\bm{\Sigma}_{\bm{\mu},\bm{\mu}}\end{bmatrix}\Big{)}\,,$
(2.6)

and
$\operatorname{\mathbb{P}}\\{Y(\bm{h})=\ell\\}={e^{\bm{h}_{\ell}}}/{\sum_{j\in[k]}e^{\bm{h}_{j}}},~{}~{}\ell\in[k].$

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 $\bm{\Sigma}_{\bm{w},\bm{w}}$ and
$\bm{\Sigma}_{\bm{w},\bm{\mu}}$, the offset values $\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

$\displaystyle\operatorname{\mathbb{P}}_{e|c}=1-\operatorname{\mathbb{P}}\big{\\{}\bm{S}_{c}^{{1}/{2}}\,{\bm{z}}\geq\bm{t}_{c}\big{\\}},$
(2.7)

where $\bm{z}$ is a Gaussian random vector distributed as
$\mathcal{N}(\bm{0},\sigma^{2}\bm{I}_{k-1})$,
$\bm{S}_{c}\in\mathbb{R}^{{(k-1)}\times{(k-1)}}$ is a symmetric matrix such
that its $i,j$ element is given by
$[\bm{S}_{c}]_{ij}:=\langle\widehat{\bm{w}}_{c}-\widehat{\bm{w}}_{j},\widehat{\bm{w}}_{c}-\widehat{\bm{w}}_{i}\rangle$
and $\bm{t}_{c}\in\mathbb{R}^{k-1}$ a vector with entries
$[\bm{t}_{c}]_{i}:=\langle\widehat{\bm{w}}_{i}-\widehat{\bm{w}}_{c},\bm{\mu}_{c}\rangle+(\widehat{\bm{b}}_{i}-\widehat{\bm{b}}_{c})$.
Similarly, based on (2.7) the total classification error in GMM is equal to
$\operatorname{\mathbb{P}}_{e}=\sum_{\ell=1}^{k}\pi_{\ell}\operatorname{\mathbb{P}}_{e|c}=1-\sum_{\ell=1}^{k}\pi_{\ell}\operatorname{\mathbb{P}}\big{\\{}\bm{S}_{c}^{{1}/{2}}\,{\bm{z}}\geq\bm{t}_{c}\big{\\}}$.
As also detailed in the supplementary, the class-wise/total test errors for
MLM similarly depends on quantities of the form
$\operatorname{\mathbb{P}}\\{\bm{A}\bm{z}\geq\bm{t}\\}$ with $\bm{z}$ a
standard Gaussian random vector, $\bm{A}$ and $\bm{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
$\operatorname{\mathbb{P}}\\{\bm{A}\bm{z}\geq\bm{t}\\}$ once $\bm{A}$ and
$\bm{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.

###  2.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 $k$.

###### Assumption 1

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

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

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

##  3 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 $\pi_{\ell}=1/k,\ell\in[k]$. Additional results on the performance of
Weighted LS are deferred to the appendix.

###  3.1 Class averaging classifier

###### Proposition 3.1

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

$\displaystyle\widehat{\bm{b}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\bm{\pi}\,,\quad\bm{\Sigma}_{\bm{w}\bm{\mu}}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\text{diag}(\bm{\pi})\,\cdot\bm{\Sigma}_{\bm{\mu},\bm{\mu}}\,,$
(3.1a) $\displaystyle\bm{\Sigma}_{\bm{w},\bm{w}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\gamma\sigma^{2}\cdot\text{diag}(\bm{\pi})+\text{diag}(\bm{\pi})\cdot\bm{\Sigma}_{\bm{\mu},\bm{\mu}}\cdot\text{diag}(\bm{\pi})\,.$
(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.

Two classes. Consider the special case with two classes with class priors
$\pi_{1}=1-\pi_{2}=:\pi$. In this case we can compute the class-wise
misclassification probabilities $\operatorname{\mathbb{P}}_{e|1}$ and
$\operatorname{\mathbb{P}}_{e|2}$ explicitly. Specifically using (3.1), we
have
$\bm{S}_{1}=\left\|\pi\bm{\mu}_{1}-(1-\pi)\bm{\mu}_{2}\right\|_{\ell_{2}}^{2}+\gamma\sigma^{2}$
and
$t_{1}=(1-2\pi)+(1-\pi)\langle\bm{\mu}_{1},\bm{\mu}_{2}\rangle-\pi\left\|\bm{\mu}_{1}\right\|_{\ell_{2}}^{2}.$
Substituting the latter two in (2.7) we arrive at
$\operatorname{\mathbb{P}}_{e|1}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}Q\big{(}\frac{\pi\left\|\bm{\mu}_{1}\right\|_{\ell_{2}}^{2}-(1-\pi)\langle\bm{\mu}_{1},\bm{\mu}_{2}\rangle+2\pi-1}{\sqrt{\left\|\pi\bm{\mu}_{1}-(1-\pi)\bm{\mu}_{2}\right\|_{\ell_{2}}^{2}+\gamma\sigma^{2}}}\big{)}\,.$
In the case of equal priors $\pi=\pi_{1}=\pi_{2}=1/2$, antipodal and equal
energy of the means, i.e., $\bm{\mu}_{1}=-\bm{\mu}_{2}$ and
$\mu:=\left\|\bm{\mu}_{1}\right\|_{\ell_{2}}=\left\|\bm{\mu}_{2}\right\|_{\ell_{2}}$,
we can use the above to conclude that
$\operatorname{\mathbb{P}}_{e|1}=\operatorname{\mathbb{P}}_{e|2}=\frac{1}{2}\operatorname{\mathbb{P}}_{e}=\frac{1}{2}Q\big{(}\sqrt{\frac{\mu^{2}}{\mu^{2}+\gamma\sigma^{2}}}\big{)}.$
This formula recovers the result of [MKLZ20] for this special case. Also, as
mentioned in [MKLZ20], the formula matches the Bayes optimal error computed in
[LM19] for Gaussian mean vectors. This shows that the class averaging method
is Bayes optimal in this very simple setting. In Section 3.3, we generalize
this result to multiple classes: we show that the average estimator is
(asymptotically) Bayes optimal for balanced classes and equal-energy Gaussian
means for any $k\geq 2$.

Orthogonal means, equal priors and equal energy. Next we focus on a special
case with orthogonal means $\langle\bm{\mu}_{i},\bm{\mu}_{j}\rangle=0,~{}i\neq
j\in[k]$ of equal energy
$\mu^{2}:=\left\|\bm{\mu}_{\imath}\right\|_{\ell_{2}}^{2}$ and of equal priors
$\pi_{i}=\pi={1}\big{/}{k}$ for $i\in[k]$. In this case, the class-wise miss-
classification error converges to
$\operatorname{\mathbb{P}}_{e|c}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}1-\operatorname{\mathbb{P}}\\{\bm{S}_{c}^{{1}/{2}}\,{\bm{z}}>\bm{t}\\},$
where
$\bm{S}_{c}=\pi(\pi\mu^{2}+\gamma\sigma^{2})({\bm{I}}_{k-1}+\bm{1}_{k-1}\bm{1}_{k-1}^{T})$
and $\bm{t}=-\pi\mu^{2}\bm{1}_{k-1}$. Defining

$u_{\rm
Avg}:=\frac{\mu^{2}}{\sigma}\sqrt{\frac{1}{\mu^{2}+k{\gamma\sigma^{2}}}}\,,$

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

$\operatorname{\mathbb{P}}_{e|c}=\operatorname{\mathbb{P}}_{e,\rm
Avg}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\operatorname{\mathbb{P}}\big{\\{}G_{0}+\max_{j\in[k-1]}G_{j}\geq\,u_{\rm
Avg}\big{\\}},$

where $G_{0},\ldots,G_{k-1}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}(0,1).$

###  3.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

$\displaystyle\bm{\Sigma}_{\bm{\mu},\bm{\mu}}={\bm{M}}^{T}{\bm{M}}={\bm{V}}\bm{\Sigma}^{2}{\bm{V}}^{T},\qquad\bm{\Sigma}\succ\mathbf{0}_{r\times
r},~{}{\bm{V}}\in\mathbb{R}^{k\times r},~{}r\leq k.$ (3.2)

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

###### Theorem 3.2

Consider data generated according to GMM in an asymptotic regime with
$\gamma<1$. In addition to (3.2), define the following two positive
(semi)-definite matrices:
${\bm{P}}:=\text{diag}(\bm{\pi})-\bm{\pi}\bm{\pi}^{T}\succeq\bm{0}_{k\times
k}$ and
$\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
$\left(\widehat{\bm{W}},\widehat{\bm{b}}\right)$ the following limits are true
asymptotically

$\displaystyle\widehat{\bm{b}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\bm{\pi}-{\bm{P}}{\bm{V}}\bm{\Sigma}\bm{\Delta}^{-1}\bm{\Sigma}{\bm{V}}^{T}\bm{\pi}\,,\quad\bm{\Sigma}_{\bm{w},\bm{\mu}}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}{\bm{P}}{\bm{V}}\bm{\Sigma}\bm{\Delta}^{-1}\bm{\Sigma}{\bm{V}}^{T}\,,$
(3.3a) $\displaystyle\bm{\Sigma}_{\bm{w},\bm{w}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\frac{\gamma}{(1-\gamma)\sigma^{2}}{\bm{P}}+{\bm{P}}{\bm{V}}\bm{\Sigma}\bm{\Delta}^{-1}\Big{(}\bm{\Delta}^{-1}-\frac{\gamma}{(1-\gamma)\sigma^{2}}{\bm{I}}_{r}\Big{)}\bm{\Sigma}{\bm{V}}^{T}{\bm{P}}\,.$
(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 $\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.

###### Proposition 3.3

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

$\operatorname{\mathbb{P}}_{e,\rm
LS}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\operatorname{\mathbb{P}}\big{\\{}G_{0}+\max_{j\in[k-1]}G_{j}\geq
u_{\rm LS}\big{\\}}.$

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

###  3.3 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
$\hat{Y}=\arg\max_{\ell\in[k]}\operatorname{\mathbb{P}}\\{Y=\ell~{}|~{}{\bm{X}},{\bm{Y}},\bm{x}\\}$
minimizes the risk
$\operatorname{\mathbb{P}}_{e}=\operatorname{\mathbb{P}}\\{\hat{Y}\neq
Y\\}=\operatorname{\mathbb{E}}_{{\bm{X}},{\bm{Y}},\bm{x},Y}\left[\mathbb{1}[\hat{Y}\neq
Y]\right]$.

###### Proposition 3.4

Consider $\bm{\mu}_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}(\bm{0},\frac{\mu^{2}}{d}\text{\em I}_{d})$ and
$\pi_{i}={1}/{k}$ for all $i\in[k]$. Set $u_{\rm
Bayes}:=\frac{\mu^{2}}{\sigma}\frac{1}{\sqrt{\mu^{2}+k{\gamma\sigma^{2}}}}.$
Then, the Bayes risk converges to
$\operatorname{\mathbb{P}}\Big{\\{}G_{0}+\max_{\ell\in[k-1]}G_{\ell}\geq
u_{\rm Bayes}\Big{\\}}.$

Under Gaussian prior, the means are asymptotically orthogonal and equal-
energy. As shown earlier, in this setting, $\operatorname{\mathbb{P}}_{e,\rm
Avg}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\operatorname{\mathbb{P}}\big{\\{}G_{0}+\max_{\ell\in[k-1]}G_{\ell}\geq
u_{\rm Avg}\big{\\}}$. But, $u_{\rm Avg}=u_{\rm 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.

##  4 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
$\bm{\Sigma}_{\bm{\mu},\bm{\mu}}={\bm{V}}\bm{\Sigma}^{2}{\bm{V}}^{T}$ in
(3.2). In order to state our results, it is convenient to introduce the
following probability vectors in $\mathbb{R}^{k}$ and $\mathbb{R}^{k^{2}}$:

$\displaystyle\bm{\pi}:=\operatorname{\mathbb{E}}\Big{[}\frac{e^{{\bm{V}}\bm{\Sigma}{\bm{g}}}}{\bm{1}_{k}^{T}e^{{\bm{V}}\bm{\Sigma}{\bm{g}}}}\Big{]}\in\mathbb{R}^{k}~{}~{}\text{and}~{}~{}\bm{\Pi}:=\operatorname{\mathbb{E}}\Big{[}\frac{\left(e^{{\bm{V}}\bm{\Sigma}{\bm{g}}}\right)\left(e^{{\bm{V}}\bm{\Sigma}{\bm{g}}}\right)^{T}}{\left(\bm{1}_{k}^{T}e^{{\bm{V}}\bm{\Sigma}{\bm{g}}}\right)^{2}}\Big{]}\in\mathbb{R}^{k\times
k},~{}~{}\text{where}~{}{\bm{g}}\sim\mathcal{N}(\bm{0},{\bm{I}}_{r}).$ (4.1)

Note that $\bm{\pi}$ and $\bm{\Pi}$ are the first and second moments of the
soft-max mapping of
${\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

$\operatorname{\mathbb{P}}\\{Y=\ell\\}=\operatorname{\mathbb{E}}[\operatorname{\mathbb{P}}\\{Y=\ell\,|\,\bm{x}\\}]=\operatorname{\mathbb{E}}\big{[}\frac{e^{\bm{e}_{\ell}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}}{\bm{1}_{k}^{T}e^{{\bm{V}}\bm{\Sigma}{\bm{g}}}}\big{]}=\bm{\pi}_{\ell},~{}\ell\in[k]$

since ${\bm{M}}^{T}\bm{x}$ is distributed as ${{\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).

###  4.1 Class averaging classifier

###### Proposition 4.1

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

$\displaystyle\widehat{\bm{b}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\bm{\pi}\,,\quad\bm{\Sigma}_{\bm{w},\bm{\mu}}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right)\cdot{\bm{\Sigma}_{\bm{\mu},\bm{\mu}}}\,,$
(4.2a) $\displaystyle\bm{\Sigma}_{\bm{w},\bm{w}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\gamma\cdot\text{diag}(\bm{\pi})+\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right)\cdot{\bm{\Sigma}_{\bm{\mu},\bm{\mu}}}\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right)\,.$
(4.2b)

Using Gaussian decomposition in (2.6) and checking from (4.2) that
$\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:

$\displaystyle\operatorname{\mathbb{P}}_{e,\rm
Avg}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\operatorname{\mathbb{P}}\big{\\{}\arg\max\left\\{\sqrt{\gamma}\cdot\text{diag}(\sqrt{\bm{\pi}})\cdot\widetilde{{\bm{g}}}+\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right)\cdot{\bm{V}}\bm{\Sigma}\cdot{\bm{g}}+\bm{\pi}\right\\}\neq
Y({\bm{g}})\big{\\}},$ (4.3)

where $\widetilde{{\bm{g}}}\sim\mathcal{N}(\bm{0},{\bm{I}}_{k})$,
${\bm{g}}\sim\mathcal{N}(\bm{0},{\bm{I}}_{r})$ and
$\operatorname{\mathbb{P}}\\{Y({\bm{g}})=c\\}={e^{\bm{e}_{c}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}}\big{/}{\sum_{j\in[k]}e^{\bm{e}_{j}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}},c\in[k]$.

###  4.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.

###### Theorem 4.2

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

$\displaystyle\widehat{\bm{b}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\bm{\pi}\,,\quad\bm{\Sigma}_{\bm{w},\bm{\mu}}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right)\cdot{\bm{\Sigma}_{\bm{\mu},\bm{\mu}}}\,,$
(4.4a) $\displaystyle\bm{\Sigma}_{\bm{w},\bm{w}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\frac{\gamma}{1-\gamma}\cdot\left(\text{diag}(\bm{\pi})-\bm{\pi}\bm{\pi}^{T}\right)+\frac{1-2\gamma}{1-\gamma}\cdot\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right)\cdot\bm{\Sigma}_{\bm{\mu},\bm{\mu}}\cdot\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right)\,.$
(4.4b)

It is interesting to observe that (4.4a) is identical to (4.2a). However, the
cross-correlations in $\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

$\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}}\frac{\gamma}{1-\gamma}\left(\text{diag}(\bm{\pi})-\bm{\pi}\bm{\pi}^{T}-\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right)\bm{\Sigma}_{\bm{\mu},\bm{\mu}}\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right)\right).$

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

###### Proposition 4.3

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

$\gamma_{\star}=\frac{\mu^{2}k}{(k-1)^{2}}\,\Big{(}\,1-k\operatorname{\mathbb{E}}\Big{[}\frac{e^{2\mu
G_{1}}}{\left(\sum_{\ell\in[k]}e^{\mu
G_{\ell}}\,\right)^{2}}\Big{]}\,\Big{)}^{2}\in(0,1).$

Then, with probability 1 as $n\rightarrow\infty$,
$\operatorname{\mathbb{P}}_{e,\rm LS}<\operatorname{\mathbb{P}}_{e,\rm
Avg}\Longleftrightarrow\gamma<\gamma_{\star}.$

##  5 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 $i$th class
weighted by $\omega_{\ell}^{2}=1/\pi_{\ell}$, (iv) Cross-Entropy (CE).

Figures 1 and 2 focus on GMM with $k=9$ classes, $d=300$ and
$\|{\bm{\mu}_{i}}\|_{\ell_{2}}^{2}=15$. To model different class prior
probabilities, we use the distribution
$\pi_{1}=\pi_{2}=\pi_{3}=0.5,\pi_{4}=0.5,\pi_{5}=0.5,\pi_{6}=0.25,\pi_{7}=0.25,\pi_{8}=0.25,\pi_{9}=1/21.$
We consider three scenarios: (a) orthogonal means, equal prior
($\pi_{i}=1/9$); (b) orthogonal means, different prior; (c) correlated means
with pairwise correlation coefficient equal to $0.5$ (i.e.,
$\langle\bm{\mu}_{i},\bm{\mu}_{j}\rangle/(\left\|\bm{\mu}_{i}\right\|_{\ell_{2}}\left\|\bm{\mu}_{j}\right\|_{\ell_{2}})=0.5$
for $i\neq j$) and different priors as discussed above. Figure 1 shows the
test miss-classification errors as a function of $\gamma:=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 $\gamma$). 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 $\gamma$
in Figure 1 ($\gamma=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 $\|{\bm{\mu}_{i}}\|_{\ell_{2}}^{2}=15$ and $d\in\\{50,100,200\\}$ with
$kd/n=k\gamma$ fixed at $k\gamma=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 $\|{\mu_{i}}\|_{\ell_{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 $\gamma$. Finally, Figure 4 (c) shows the
class-wise probability of error associated with (b) for $\gamma=0.117$ and
demonstrates that LS outperforms averaging.

Prob. of Error$\gamma$

$\gamma$

$\gamma$

Figure 1: GMM with $k=9,d=300$. (a) orthogonal, equal prior, (b) orthogonal,
different prior, (c) correlated, different prior.

Class-wise Prob. of ErrorClass ID

Class ID

Class ID

Figure 2: Class-wise probability of errors corresponding to Figure 1 with
$\gamma=0.117$.

$\frac{\text{Empirical Prob. Error}}{\text{Theory Prob. Error}}$# of Classes /
d

# of Classes / d

# of Classes / d

Figure 3: GMM, $K/d$ is varied from $0$ to $1$ while keeping $K\gamma$
constant for (a) $d=50$, (b) $d=100$, (c) $d=200$.

Prob. of Error$\gamma$

$\gamma$

Class-wise Prob. of ErrorClass ID

Figure 4: MLM with orthogonal means for (a) equal norms and (b) different
norms. (c) Class-wise probability errors for (b).

##  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:
$\bm{b},\bm{\Sigma}_{\bm{w},\bm{\mu}},$ and, $\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
$\bm{x}\mapsto\bm{W}\bm{x}+\bm{b}$ via a least-squares fit to the training
data:

$\displaystyle(\widehat{\bm{W}},\widehat{\bm{b}}):=\frac{1}{2n}\left\|\bm{W}\bm{X}+\bm{b}\bm{1}_{n}^{T}-\bm{Y}\right\|_{F}^{2}.$
(6.1)

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

$\displaystyle\frac{1}{2n}\left\|\bm{W}\bm{X}+\bm{b}\bm{1}_{n}^{T}-\bm{Y}\right\|_{F}^{2}=\frac{1}{2n}\sum_{\ell=1}^{k}\left\|\bm{X}^{T}\bm{w}_{\ell}+b_{\ell}\bm{1}_{n}-\bm{Y}_{\ell}\right\|_{\ell_{2}}^{2}.$

Hence, for each $\ell\in[k]$,

$\displaystyle(\widehat{\bm{w}}_{\ell},\widehat{b}_{\ell})=\arg\min_{\bm{w}_{\ell},{b}_{\ell}}\text{
}\frac{1}{2n}\left\|\bm{X}^{T}\bm{w}_{\ell}+b_{\ell}\bm{1}_{n}-\bm{Y}_{\ell}\right\|_{\ell_{2}}^{2}.$
(6.2)

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 $\ell\in[k]$: (i) the
intercept $\widehat{\bm{b}}_{\ell}$, (ii) the mean-correlations
$\langle\widehat{\bm{w}}_{\ell},\bm{\mu}_{c}\rangle,~{}c\in[k]$, (iii) the
norm $\left\|\widehat{\bm{w}}_{\ell}\right\|_{\ell_{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

$\displaystyle\frac{1}{2}\left(\min_{\bm{w}_{\ell},{b}_{\ell}}\text{
}\frac{1}{\sqrt{n}}\left\|\sigma\left\|\bm{w}_{\ell}\right\|_{\ell_{2}}\bm{g}+\bm{Y}^{T}\bm{M}^{T}\bm{w}_{\ell}+b_{\ell}\bm{1}_{n}-\bm{Y}_{\ell}\right\|_{\ell_{2}}+\frac{1}{\sqrt{n}}\sigma\bm{h}^{T}\bm{w}_{\ell}\right)_{+}^{2},$
(6.3)

where $(x)_{+}:=\max(0,x)$ and $\bm{g}\in\mathbb{R}^{n}$ and
$\bm{h}\in\mathbb{R}^{d}$ are two independent Gaussian random vectors
distributed as $\mathcal{N}(\bm{0},\bm{I}_{n})$ and
$\mathcal{N}(\bm{0},\bm{I}_{d})$.

Step III: Simplification of the AO and computing
$\bm{\Sigma}_{\bm{w},\bm{\mu}}$ and $\bm{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 intercepts
$\widehat{\bm{b}}_{\ell},~{}\ell\in[k]$, the mean-correlations
$\langle\widehat{\bm{w}}_{\ell},\bm{\mu}_{c}\rangle,~{}\ell,c\in[k]$, the
norms $\left\|\widehat{\bm{w}}_{\ell}\right\|_{\ell_{2}},~{}\ell\in[k]$ and
the LS training loss
$\left\|\bm{X}^{T}\widehat{\bm{w}}_{\ell}+\widehat{\bm{b}}_{\ell}\bm{1}_{n}-\bm{Y}_{\ell}\right\|_{\ell_{2}}$.
In particular, for GMM these calculations allow us to conclude the following
limits for all $\ell\in[k]$:

$\displaystyle\widehat{b}_{\ell}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\pi_{\ell}\left(1-(\bm{e}_{\ell}-\bm{\pi})^{T}{\bm{V}}\bm{\Sigma}\bm{\Delta}^{-1}\bm{\Sigma}{\bm{V}}^{T}\bm{\pi}\right),\qquad{\bm{M}}^{T}\widehat{\bm{w}}_{\ell}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}{\pi_{\ell}}{\bm{V}}\bm{\Sigma}\bm{\Delta}^{-1}\bm{\Sigma}{\bm{V}}^{T}\left(\bm{e}_{\ell}-\bm{\pi}\right),\quad$
(6.4)

and

$\displaystyle\left\|\bm{w}_{\ell}\right\|_{\ell_{2}}^{2}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\frac{\gamma}{(1-\gamma)\sigma^{2}}\pi_{\ell}(1-\pi_{\ell})+\pi_{\ell}^{2}\left(\bm{e}_{\ell}-\bm{\pi}\right)^{T}{\bm{V}}\bm{\Sigma}\bm{\Delta}^{-1}\left(\bm{\Delta}^{-1}-\frac{\gamma}{(1-\gamma)\sigma^{2}}{\bm{I}}_{r}\right)\bm{\Sigma}{\bm{V}}^{T}\left(\bm{e}_{\ell}-\bm{\pi}\right)\,.$
(6.5)

where
$\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 ${\bm{P}}:=\text{diag}(\bm{\pi})-\bm{\pi}\bm{\pi}^{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
$\Sigma_{\bm{w},\bm{w}}$. Note that (6.5) already determines the diagonal
entries of $\Sigma_{\bm{w},\bm{w}}$. However, thus far, our analysis treats
the optimization of each classifier $\widehat{\bm{w}}_{\ell},~{}\ell\in[k]$
independently and provides no information for the cross-correlation
$\langle\widehat{\bm{w}}_{\ell},\widehat{\bm{w}}_{c}\rangle,\ell\neq c\in[k]$

Step IV: Computing $\bm{\Sigma}_{\bm{w},\bm{w}}$ and capturing cross-
correlations. The final and most involved part of our analysis is
characterizing the asymptotic behavior of $\bm{\Sigma}_{w,w}$. To see why this
is particularly challenging note that the reduction from (6.1) to (6.2)
“breaks" the dependence of all
$\widehat{\bm{w}}_{1},\widehat{\bm{w}}_{2},\ldots,\widehat{\bm{w}}_{k}$ on the
_same_ feature matrix ${\bm{X}}$. Capturing this dependence is crucial in
determining the “cross-correlations"
$\langle\widehat{\bm{w}}_{\ell},\widehat{\bm{w}}_{c}\rangle,~{}\ell\neq c$. As
noted in Section 2.3 the matrix $\bm{\Sigma}_{\bm{w},\bm{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
$\widehat{\bm{w}}_{\ell,c}=\widehat{\bm{w}}_{\ell}+\widehat{\bm{w}}_{c}$ is
itself the solution to another simple single-output LS problem.

###### Lemma 6.1

For $\ell\neq c\in[k]$, let $\widehat{\bm{w}}_{\ell}$, $\widehat{\bm{w}}_{c}$
be the $\ell$ and $c$-th row of $\widehat{\bm{W}}$ which is the solution to
the multi-output least-squares minimization (6.1). Denote
$\widehat{\bm{w}}_{\ell,c}:=\widehat{\bm{w}}_{\ell}+\widehat{\bm{w}}_{c}$.
Then, $\widehat{\bm{w}}_{\ell,c}$ is a minimizer in the following single-
output least-squares problem:

$\displaystyle\widehat{\bm{w}}_{\ell,c}=\arg\min_{\bm{w},b}\frac{1}{2n}\left\|{\bm{Y}}_{\ell}+{\bm{Y}}_{c}-{\bm{X}}^{T}\bm{w}-b\bm{1}_{n}\right\|_{\ell_{2}}^{2}\,.$

Thanks to Lemma 6.1, we can use the CGMT to characterize the limiting behavior
of $\left\|\widehat{\bm{w}}_{\ell}+\widehat{\bm{w}}_{c}\right\|_{\ell_{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
$\left\|\widehat{\bm{w}}_{\ell}+\widehat{\bm{w}}_{c}\right\|_{\ell_{2}}$
immediately yields the asymptotic characterization of
$\langle\widehat{\bm{w}}_{\ell},\widehat{\bm{w}}_{c}\rangle$ as

$\displaystyle\langle\widehat{\bm{w}}_{\ell},\widehat{\bm{w}}_{c}\rangle=\frac{\left\|\widehat{\bm{w}}_{\ell}+\widehat{\bm{w}}_{c}\right\|_{\ell_{2}}^{2}-\left\|\widehat{\bm{w}}_{\ell}\right\|_{\ell_{2}}^{2}-\left\|\widehat{\bm{w}}_{c}\right\|_{\ell_{2}}^{2}}{2},$
(6.6)

and
$\left\|\widehat{\bm{w}}_{\ell}\right\|_{\ell_{2}},\left\|\widehat{\bm{w}}_{c}\right\|_{\ell_{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 $\bm{\Sigma}_{w,w}$ as
promised in (3.3b) in Theorem 3.2:

$\displaystyle\bm{\Sigma}_{\bm{w},\bm{w}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\frac{\gamma}{(1-\gamma)\sigma^{2}}{\bm{P}}+{\bm{P}}{\bm{V}}\bm{\Sigma}\bm{\Delta}^{-1}\Big{(}\bm{\Delta}^{-1}-\frac{\gamma}{(1-\gamma)\sigma^{2}}{\bm{I}}_{r}\Big{)}\bm{\Sigma}{\bm{V}}^{T}{\bm{P}}.$

###  6.1 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:

$\displaystyle X_{\bm{w},{\bm{u}}}$
$\displaystyle:={\bm{u}}^{T}\bm{G}\bm{w}+\psi(\bm{w},{\bm{u}}),$ (6.7a)
$\displaystyle Y_{\bm{w},{\bm{u}}}$
$\displaystyle:=\left\|\bm{w}\right\|_{\ell_{2}}{\bm{g}}^{T}{\bm{u}}+\left\|{\bm{u}}\right\|_{\ell_{2}}\bm{h}^{T}\bm{w}+\psi(\bm{w},{\bm{u}}),$
(6.7b)

where: $\bm{G}\in\mathbb{R}^{n\times d}$, ${\bm{g}}\in\mathbb{R}^{n}$,
$\bm{h}\in\mathbb{R}^{d}$, they all have entries iid Gaussian; the sets
$\mathcal{S}_{\bm{w}}\subset\mathbb{R}^{d}$ and
$\mathcal{S}_{{\bm{u}}}\subset\mathbb{R}^{n}$ are compact; and,
$\psi:\mathbb{R}^{d}\times\mathbb{R}^{n}\to\mathbb{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:

$\displaystyle\Phi(\bm{G})$
$\displaystyle=\min\limits_{\bm{w}\in\mathcal{S}_{\bm{w}}}\max\limits_{{\bm{u}}\in\mathcal{S}_{{\bm{u}}}}X_{\bm{w},{\bm{u}}},$
(6.8a) $\displaystyle\phi({\bm{g}},\bm{h})$
$\displaystyle=\min\limits_{\bm{w}\in\mathcal{S}_{\bm{w}}}\max\limits_{{\bm{u}}\in\mathcal{S}_{{\bm{u}}}}Y_{\bm{w},{\bm{u}}}.$
(6.8b)

If the sets $\mathcal{S}_{\bm{w}}$ and $\mathcal{S}_{{\bm{u}}}$ are convex and
_bounded_ , and $\psi$ is continuous _convex-concave_ on
$\mathcal{S}_{\bm{w}}\times\mathcal{S}_{{\bm{u}}}$, then, for any
$\nu\in\mathbb{R}$ and $t>0$, it holds [TOH15, Thm. 3]:

$\mathbb{P}\left(\left|\Phi(\bm{G})-\nu\right|>t\right)\leq
2\,\mathbb{P}\left(\left|\phi({\bm{g}},\bm{h})-\nu\right|>t\right).$ (6.9)

In words, concentration of the optimal cost of the AO problem around
$q^{\ast}$ implies concentration of the optimal cost of the corresponding PO
problem around the same value $q^{\ast}$. Asymptotically, if we can show that
$\phi({\bm{g}},\bm{h})\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}q^{\ast}$,
then we can conclude that
$\Phi(\bm{G})\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}q^{\ast}$.
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
$\left\|\bm{w}_{\phi}({\bm{g}},\bm{h})\right\|_{\ell_{2}}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\alpha^{\ast}$
for some $\alpha^{\ast}>0$, then, the same holds true for the minimizers of
(6.8a):
$\left\|\bm{w}_{\Phi}(\bm{G})\right\|_{\ell_{2}}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\alpha^{\ast}$.
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
$\bm{\Sigma}_{\bm{w}\bm{w}}$.

##  7 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
$\bm{\Sigma}_{\bm{w},\bm{\mu}}$ 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.

## 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$, the Air Force Office of Scientific Research
Young Investigator Program (AFOSR-YIP) under award $\\#$FA$9550-18-1-0078$,
DARPA Learning with Less Labels (LwLL) and FastNICS programs, and NSF-CIF
awards $\\#1813877$ and $\\#2008443$.

## References

  * [AG82] Per Kragh Andersen and Richard D Gill.  Cox’s regression model for counting processes: a large sample study.  The annals of statistics, pages 1100–1120, 1982.
  * [ALMT13] Dennis Amelunxen, Martin Lotz, Michael B McCoy, and Joel A Tropp.  Living on the edge: A geometric theory of phase transitions in convex optimization.  arXiv preprint arXiv:1303.6672, 2013.
  * [ASH19] Ehsan Abbasi, Fariborz Salehi, and Babak Hassibi.  Universality in learning from linear measurements.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 12372–12382. Curran Associates, Inc., 2019.
  * [ASS00] Erin L Allwein, Robert E Schapire, and Yoram Singer.  Reducing multiclass to binary: A unifying approach for margin classifiers.  Journal of machine learning research, 1(Dec):113–141, 2000.
  * [BB99] Erin J Bredensteiner and Kristin P Bennett.  Multicategory classification by support vector machines.  In Computational Optimization, pages 53–79. Springer, 1999.
  * [BEH20] Anna Sergeevna Bosman, Andries Engelbrecht, and Mardé Helbig.  Visualising basins of attraction for the cross-entropy and the squared error neural network loss functions.  Neurocomputing, 2020.
  * [BKRS19] Zhiqi Bu, Jason Klusowski, Cynthia Rush, and Weijie Su.  Algorithmic analysis and statistical estimation of slope via approximate message passing.  In Advances in Neural Information Processing Systems, pages 9361–9371, 2019.
  * [BM12] Mohsen Bayati and Andrea Montanari.  The lasso risk for gaussian matrices.  Information Theory, IEEE Transactions on, 58(4):1997–2017, 2012\.
  * [CAL13] Anna Choromanska, Alekh Agarwal, and John Langford.  Extreme multi class classification.  In NIPS Workshop: eXtreme Classification, submitted, 2013.
  * [CAS16] Paul Covington, Jay Adams, and Emre Sargin.  Deep neural networks for youtube recommendations.  In Proceedings of the 10th ACM conference on recommender systems, pages 191–198, 2016.
  * [CKMY16] Corinna Cortes, Vitaly Kuznetsov, Mehryar Mohri, and Scott Yang.  Structured prediction theory based on factor graph complexity.  In Advances in Neural Information Processing Systems, pages 2514–2522, 2016.
  * [CM19] Michael Celentano and Andrea Montanari.  Fundamental barriers to high-dimensional regression with convex penalties.  arXiv preprint arXiv:1903.10603, 2019.
  * [CRPW12] Venkat Chandrasekaran, Benjamin Recht, Pablo A Parrilo, and Alan S Willsky.  The convex geometry of linear inverse problems.  Foundations of Computational Mathematics, 12(6):805–849, 2012.
  * [CS01] Koby Crammer and Yoram Singer.  On the algorithmic implementation of multiclass kernel-based vector machines.  Journal of machine learning research, 2(Dec):265–292, 2001.
  * [CS18] Emmanuel J Candès and Pragya Sur.  The phase transition for the existence of the maximum likelihood estimate in high-dimensional logistic regression.  arXiv preprint arXiv:1804.09753, 2018.
  * [DB94] Thomas G Dietterich and Ghulum Bakiri.  Solving multiclass learning problems via error-correcting output codes.  Journal of artificial intelligence research, 2:263–286, 1994.
  * [DCO20] Ahmet Demirkaya, Jiasi Chen, and Samet Oymak.  Exploring the role of loss functions in multiclass classification.  In 2020 54th Annual Conference on Information Sciences and Systems (CISS), pages 1–5. IEEE, 2020.
  * [DKT19] Zeyu Deng, Abla Kammoun, and Christos Thrampoulidis.  A model of double descent for high-dimensional binary linear classification.  arXiv preprint arXiv:1911.05822, 2019.
  * [DM16] David Donoho and Andrea Montanari.  High dimensional robust m-estimation: Asymptotic variance via approximate message passing.  Probability Theory and Related Fields, 166(3-4):935–969, 2016.
  * [DMM11] David L Donoho, Arian Maleki, and Andrea Montanari.  The noise-sensitivity phase transition in compressed sensing.  Information Theory, IEEE Transactions on, 57(10):6920–6941, 2011\.
  * [Don06] David L Donoho.  Compressed sensing.  Information Theory, IEEE Transactions on, 52(4):1289–1306, 2006\.
  * [EK18] Noureddine El Karoui.  On the impact of predictor geometry on the performance on high-dimensional ridge-regularized generalized robust regression estimators.  Probability Theory and Related Fields, 170(1-2):95–175, 2018.
  * [FR13] Simon Foucart and Holger Rauhut.  An invitation to compressive sensing.  In A mathematical introduction to compressive sensing, pages 1–39. Springer, 2013.
  * [Für02] Johannes Fürnkranz.  Round robin classification.  Journal of Machine Learning Research, 2(Mar):721–747, 2002.
  * [GCOZ17] Krzysztof Gajowniczek, Leszek J Chmielewski, Arkadiusz Orłowski, and Tomasz Ząbkowski.  Generalized entropy cost function in neural networks.  In International Conference on Artificial Neural Networks, pages 128–136. Springer, 2017.
  * [Gor88] Yehoram Gordon.  On Milman’s inequality and random subspaces which escape through a mesh in $\mathbb{R}^{n}$.  Springer, 1988.
  * [Gue02] Yann Guermeur.  Combining discriminant models with new multi-class svms.  Pattern Analysis & Applications, 5(2):168–179, 2002.
  * [HH03] Enkelejd Hashorva and Jürg Hüsler.  On multivariate gaussian tails.  Annals of the Institute of Statistical Mathematics, 55(3):507–522, 2003.
  * [HL19] Hong Hu and Yue M Lu.  Asymptotics and optimal designs of slope for sparse linear regression.  arXiv preprint arXiv:1903.11582, 2019.
  * [Hua17] Hanwen Huang.  Asymptotic behavior of support vector machine for spiked population model.  The Journal of Machine Learning Research, 18(1):1472–1492, 2017\.
  * [HYS16] Le Hou, Chen-Ping Yu, and Dimitris Samaras.  Squared earth mover’s distance-based loss for training deep neural networks.  arXiv preprint arXiv:1611.05916, 2016.
  * [JGP16] Eric Jang, Shixiang Gu, and Ben Poole.  Categorical reparameterization with gumbel-softmax.  arXiv preprint arXiv:1611.01144, 2016.
  * [JSH20] Adel Javanmard, Mahdi Soltanolkotabi, and Hamed Hassani.  Precise tradeoffs in adversarial training for linear regression.  arXiv preprint arXiv:2002.10477, 2020.
  * [KA20] Abla Kammoun and Mohamed-Slim Alouini.  On the precise error analysis of support vector machines.  arXiv preprint arXiv:2003.12972, 2020.
  * [Kar13] Noureddine El Karoui.  Asymptotic behavior of unregularized and ridge-regularized high-dimensional robust regression estimators: rigorous results.  arXiv preprint arXiv:1311.2445, 2013.
  * [KMS15] Vitaly Kuznetsov, Mehryar Mohri, and U Syed.  Rademacher complexity margin bounds for learning with a large number of classes.  In ICML Workshop on Extreme Classification: Learning with a Very Large Number of Labels, 2015.
  * [KP+02] Vladimir Koltchinskii, Dmitry Panchenko, et al.  Empirical margin distributions and bounding the generalization error of combined classifiers.  The Annals of Statistics, 30(1):1–50, 2002.
  * [KS18] Himanshu Kumar and PS Sastry.  Robust loss functions for learning multi-class classifiers.  In 2018 IEEE International Conference on Systems, Man, and Cybernetics (SMC), pages 687–692. IEEE, 2018.
  * [KT20] Ganesh Kini and Christos Thrampoulidis.  Analytic study of double descent in binary classification: The impact of loss.  arXiv preprint arXiv:2001.11572, 2020.
  * [LDBK15] Yunwen Lei, Urun Dogan, Alexander Binder, and Marius Kloft.  Multi-class svms: From tighter data-dependent generalization bounds to novel algorithms.  In Advances in Neural Information Processing Systems, pages 2035–2043, 2015.
  * [LDZK19] Yunwen Lei, Ürün Dogan, Ding-Xuan Zhou, and Marius Kloft.  Data-dependent generalization bounds for multi-class classification.  IEEE Transactions on Information Theory, 65(5):2995–3021, 2019\.
  * [LLW04] Yoonkyung Lee, Yi Lin, and Grace Wahba.  Multicategory support vector machines: Theory and application to the classification of microarray data and satellite radiance data.  Journal of the American Statistical Association, 99(465):67–81, 2004.
  * [LLY+18] Jian Li, Yong Liu, Rong Yin, Hua Zhang, Lizhong Ding, and Weiping Wang.  Multi-class learning: from theory to algorithm.  In Advances in Neural Information Processing Systems, pages 1586–1595, 2018.
  * [LM19] Marc Lelarge and Leo Miolane.  Asymptotic bayes risk for gaussian mixture in a semi-supervised setting.  arXiv preprint arXiv:1907.03792, 2019.
  * [Lol20] Panagiotis Lolas.  Regularization in high-dimensional regression and classification via random matrix theory.  arXiv preprint arXiv:2003.13723, 2020.
  * [LS20] Tengyuan Liang and Pragya Sur.  A precise high-dimensional asymptotic theory for boosting and min-l1-norm interpolated classifiers.  arXiv preprint arXiv:2002.01586, 2020.
  * [LT91] Michel Ledoux and Michel Talagrand.  Probability in Banach Spaces: isoperimetry and processes, volume 23.  Springer, 1991.
  * [Mau16] Andreas Maurer.  A vector-contraction inequality for rademacher complexities.  In International Conference on Algorithmic Learning Theory, pages 3–17. Springer, 2016.
  * [MKLZ20] Francesca Mignacco, Florent Krzakala, Yue M Lu, and Lenka Zdeborová.  The role of regularization in classification of high-dimensional noisy gaussian mixture.  arXiv preprint arXiv:2002.11544, 2020.
  * [MLC19a] X. Mai, Z. Liao, and R. Couillet.  A large scale analysis of logistic regression: asymptotic performance and new insights.  In ICASSP, 2019.
  * [MLC19b] Xiaoyi Mai, Zhenyu Liao, and Romain Couillet.  A large scale analysis of logistic regression: Asymptotic performance and new insights.  In ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3357–3361. IEEE, 2019.
  * [MM18] Léo Miolane and Andrea Montanari.  The distribution of the lasso: Uniform control over sparse balls and adaptive parameter tuning.  arXiv preprint arXiv:1811.01212, 2018.
  * [MRSY19] Andrea Montanari, Feng Ruan, Youngtak Sohn, and Jun Yan.  The generalization error of max-margin linear classifiers: High-dimensional asymptotics in the overparametrized regime.  arXiv preprint arXiv:1911.01544, 2019.
  * [MXSS20] Jincheng Mei, Chenjun Xiao, Csaba Szepesvari, and Dale Schuurmans.  On the global convergence rates of softmax policy gradient methods.  arXiv preprint arXiv:2005.06392, 2020.
  * [OH10] Samet Oymak and Babak Hassibi.  New null space results and recovery thresholds for matrix rank minimization.  arXiv preprint arXiv:1011.6326, 2010.
  * [ORS17] Samet Oymak, Benjamin Recht, and Mahdi Soltanolkotabi.  Sharp time–data tradeoffs for linear inverse problems.  IEEE Transactions on Information Theory, 64(6):4129–4158, 2017\.
  * [OTH13] Samet Oymak, Christos Thrampoulidis, and Babak Hassibi.  The squared-error of generalized lasso: A precise analysis.  arXiv preprint arXiv:1311.0830, 2013.
  * [PM05] Mahesh Pal and PM Mather.  Support vector machines for classification in remote sensing.  International journal of remote sensing, 26(5):1007–1011, 2005\.
  * [PS16] Bernardo Ávila Pires and Csaba Szepesvári.  Multiclass classification calibration functions.  arXiv preprint arXiv:1609.06385, 2016.
  * [PSG13] Bernardo Avila Pires, Csaba Szepesvari, and Mohammad Ghavamzadeh.  Cost-sensitive multiclass classification risk bounds.  In International Conference on Machine Learning, pages 1391–1399, 2013.
  * [RCY+19] Ankit Singh Rawat, Jiecao Chen, Felix Xinnan X Yu, Ananda Theertha Suresh, and Sanjiv Kumar.  Sampled softmax with random fourier features.  In Advances in Neural Information Processing Systems 32. 2019.
  * [RK04] Ryan Rifkin and Aldebaro Klautau.  In defense of one-vs-all classification.  Journal of machine learning research, 5(Jan):101–141, 2004.
  * [SAH19] Fariborz Salehi, Ehsan Abbasi, and Babak Hassibi.  The impact of regularization on high-dimensional logistic regression.  In Advances in Neural Information Processing Systems, pages 12005–12015, 2019.
  * [SC19] Pragya Sur and Emmanuel J Candès.  A modern maximum-likelihood theory for high-dimensional logistic regression.  Proceedings of the National Academy of Sciences, 116(29):14516–14525, 2019.
  * [SL80] YS Sathe and SR Lingras.  A note on the inequalities for tail probability of the multivariate normal distribution.  Communications in Statistics-Theory and Methods, 9(7):711–715, 1980\.
  * [Sto09] Mihailo Stojnic.  Various thresholds for $\ell_{1}$-optimization in compressed sensing.  arXiv preprint arXiv:0907.3666, 2009.
  * [Sto13] Mihailo Stojnic.  A framework to characterize performance of lasso algorithms.  arXiv preprint arXiv:1303.7291, 2013.
  * [SVL14] Ilya Sutskever, Oriol Vinyals, and Quoc V Le.  Sequence to sequence learning with neural networks.  In Advances in neural information processing systems, pages 3104–3112, 2014.
  * [TAH15] Christos Thrampoulidis, Ehsan Abbasi, and Babak Hassibi.  Lasso with non-linear measurements is equivalent to one with linear measurements.  In Advances in Neural Information Processing Systems, pages 3420–3428, 2015.
  * [TAH18] Christos Thrampoulidis, Ehsan Abbasi, and Babak Hassibi.  Precise error analysis of regularized $m$-estimators in high dimensions.  IEEE Transactions on Information Theory, 64(8):5592–5628, 2018\.
  * [TB07] Ambuj Tewari and Peter L Bartlett.  On the consistency of multiclass classification methods.  Journal of Machine Learning Research, 8(May):1007–1025, 2007.
  * [TOH15] Christos Thrampoulidis, Samet Oymak, and Babak Hassibi.  Regularized linear regression: A precise analysis of the estimation error.  In Proceedings of The 28th Conference on Learning Theory, pages 1683–1709, 2015.
  * [TPT20a] Hossein Taheri, Ramtin Pedarsani, and Christos Thrampoulidis.  Fundamental limits of ridge-regularized empirical risk minimization in high dimensions.  arXiv preprint arXiv:2006.08917, 2020.
  * [TPT20b] Hossein Taheri, Ramtin Pedarsani, and Christos Thrampoulidis.  Sharp asymptotics and optimal performance for inference in binary models.  In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pages 3739–3749. PMLR, 2020.
  * [TXH18] Christos Thrampoulidis, Weiyu Xu, and Babak Hassibi.  Symbol error rate performance of box-relaxation decoders in massive mimo.  IEEE Transactions on Signal Processing, 66(13):3377–3392, 2018\.
  * [WBW+16] Alan Wisler, Visar Berisha, Dennis Wei, Karthikeyan Ramamurthy, and Andreas Spanias.  Empirically-estimable multi-class classification bounds.  In 2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 2594–2598. IEEE, 2016.
  * [WW98] Jason Weston and Chris Watkins.  Multi-class support vector machines.  Technical report, Citeseer, 1998.
  * [WWM19] Shuaiwen Wang, Haolei Weng, and Arian Maleki.  Does slope outperform bridge regression?  arXiv preprint arXiv:1909.09345, 2019.
  * [YHR+16] Ian En-Hsu Yen, Xiangru Huang, Pradeep Ravikumar, Kai Zhong, and Inderjit Dhillon.  Pd-sparse: A primal and dual sparse approach to extreme multiclass and multilabel classification.  In International Conference on Machine Learning, pages 3069–3077, 2016.
  * [Zha04] Tong Zhang.  Statistical analysis of some multi-category large margin classification methods.  Journal of Machine Learning Research, 5(Oct):1225–1251, 2004.

###### Contents

  1. 1 Introduction
  2. 2 Problem formulation
    1. 2.1 Data Models
    2. 2.2 Classification algorithms
    3. 2.3 Class-wise and total test classification error
    4. 2.4 High-dimensional regime
  3. 3 Results for Gaussian Mixture Model
    1. 3.1 Class averaging classifier
    2. 3.2 Least-squares classifier
    3. 3.3 Bayes estimator for the balanced Gaussian Mixture Model
  4. 4 Results for Multinomial Logit Model
    1. 4.1 Class averaging classifier
    2. 4.2 Least-squares classifier
  5. 5 Numerical Results
  6. 6 Proof outline for least-squares: key ideas and challenges
    1. 6.1 Background on the CGMT
  7. 7 Future Directions
  8. A Additional Numerical Results
  9. B Additional Results on Weighted Least-squares classifiers
    1. B.1 WLS for GMM
    2. B.2 WLS for MLM
  10. C Preliminaries
    1. C.1 Slepian’s inequality
    2. C.2 Gaussian integration by parts
    3. C.3 Block matrix inversion
  11. D Calculating and bounding the missclassification error
    1. D.1 Proof of (2.5) and (2.6)
    2. D.2 Class-wise and total miss-classification error for GMM
    3. D.3 Class-wise and total miss-classification error for MLM
    4. D.4 Evaluating and bounding tail probabilities of multivariate Gaussians
      1. D.4.1 A special case: Rank-one update of Identity
      2. D.4.2 Slepian’s bound
      3. D.4.3 Simple bounds for GMM
  12. E The Class-averaging estimator
    1. E.1 Proofs for GMM
      1. E.1.1 GMM: Proof of Proposition 3.1
    2. E.2 Proofs for MLM
      1. E.2.1 Proof of Proposition 4.1
      2. E.2.2 Orthogonal means
  13. F On the Bayes risk of GMM: Proof of Proposition 3.4
  14. G Least-squares for GMM
    1. G.1 Proof of Theorem 3.2
      1. G.1.1 Computing $\Sigma_{w,\mu}$
      2. G.1.2 Computing $\Sigma_{w,w}$
    2. G.2 Orthogonal means
  15. H Least-squares for MLM
    1. H.1 Proof of Theorem 4.2
      1. H.1.1 Computing $\Sigma_{w,\mu}$
      2. H.1.2 Computing $\Sigma_{w,w}$
    2. H.2 Orthogonal means and equal-energy
    3. H.3 Proof of Proposition 4.3
  16. I Weighted LS for GMM (Proof of Theorem B.1)
    1. I.1 Computing $\Sigma_{w,\mu}$
    2. I.2 Computing $\Sigma_{w,w}$
  17. J Weighted LS for MLM (Proof of Theorem B.2)
    1. J.1 Computing $\Sigma_{w,\mu}$
    2. J.3 Computing $\Sigma_{w,w}$

##  Appendix A Additional Numerical Results

Transition point $\gamma_{\star}$Number of classes $k$ Figure 5: The threshold
$\gamma_{\star}$ of Proposition 4.3 as a function of the number of classes $k$
and the means’ energy $\mu$. LS provably outperforms class-averaging for
$\gamma<\gamma_{\star}$.

Class-wise Prob. of ErrorClass ID

Class ID

Figure 6:  Class-wise probabilities of error for MLM with (a) orthogonal means
and (b) correlated means.

Prob. of error$\gamma$

$\gamma$

$\gamma$

Figure 7: Precise theoretical predictions compared to a theoretical upper
bound (UB) obtained by Union bound that does _not_ require knowledge of the
off-diagonal entries of $\bm{\Sigma}_{\bm{w},\bm{w}}$. GMM with orthogonal
means and (a) equal priors, different strengths; (b) different priors,
different strengths; (c) different priors, equal strengths. See text for
details.

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 $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,
$\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.

Second, Figure 6 provides further experiments on the class-wise probabilities
of the MLM model with $k=9$ classes for $\gamma=0.1$. Classes 1,2,3 have norms
$\left\|\bm{\mu}_{i}\right\|_{\ell_{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:
$\omega_{i}^{2}=1/\bm{\pi}_{\ell},\ell\in[k]$, where $\bm{\pi}_{\ell}$ is the
$\ell$th entry of the vector $\bm{\pi}$ 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) $\pi_{1}=\pi_{2}=\pi_{3}=\pi_{4}=\pi_{5}$,
$4\left\|\bm{\mu}_{1}\right\|_{\ell_{2}}=4\left\|\bm{\mu}_{2}\right\|_{\ell_{2}}=2\left\|\bm{\mu}_{3}\right\|_{\ell_{2}}=2\left\|\bm{\mu}_{4}\right\|_{\ell_{2}}=\left\|\bm{\mu}_{5}\right\|_{\ell_{2}}=4\sqrt{3}$;
(b) $4\pi_{1}=4\pi_{2}=2\pi_{3}=2\pi_{4}=\pi_{5}$,
$4\left\|\bm{\mu}_{1}\right\|_{\ell_{2}}=4\left\|\bm{\mu}_{2}\right\|_{\ell_{2}}=2\left\|\bm{\mu}_{3}\right\|_{\ell_{2}}=2\left\|\bm{\mu}_{4}\right\|_{\ell_{2}}=\left\|\bm{\mu}_{5}\right\|_{\ell_{2}}=4\sqrt{3}$;
(c) $4\pi_{1}=4\pi_{2}=2\pi_{3}=2\pi_{4}=\pi_{5}$,
$\left\|\bm{\mu}_{1}\right\|_{\ell_{2}}=\left\|\bm{\mu}_{2}\right\|_{\ell_{2}}=\left\|\bm{\mu}_{3}\right\|_{\ell_{2}}=\left\|\bm{\mu}_{4}\right\|_{\ell_{2}}=\left\|\bm{\mu}_{5}\right\|_{\ell_{2}}=\sqrt{3}$.
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
${\bm{{\Sigma}}}_{\bm{w},\bm{w}}$ are unknown). These bounds are calculated by
applying a union bound to the class-wise probabilities
$\operatorname{\mathbb{P}}_{e|c},~{}c\in[k]$ in (D.1) and further
appropriately bounding the off-diagonal entries of
${\bm{{\Sigma}}}_{\bm{w},\bm{w}}$ in terms of the self-correlations of the
classes, i.e., only the diagonal entries of ${\bm{{\Sigma}}}_{\bm{w},\bm{w}}$.
Please see Section D.4.3 for details. Overall, the bounds shown only depend on
$\widehat{\bm{b}},{\bm{{\Sigma}}}_{\bm{w},\bm{\mu}}$ and
$\text{diag}({\bm{{\Sigma}}}_{\bm{w},\bm{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.

##  Appendix B Additional Results on Weighted Least-squares classifiers

###  B.1 WLS for GMM

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

###### Theorem B.1

Consider data generated according to GMM and $\gamma<1$. Consider a weighted
LS classifier with weights $\bm{D}=\text{diag}(\omega_{1},\ldots,\omega_{k})$
and let $\eta$ be the unique solution to
$\sum_{\ell=1}^{k}\frac{\pi_{\ell}\omega_{\ell}^{2}}{\omega_{\ell}^{2}+\eta}=\gamma$.
Also define
${\bm{P}}:=\text{diag}(\widetilde{\bm{\pi}})-\widetilde{\bm{\pi}}\widetilde{\bm{\pi}}^{T}\succeq\bm{0}_{k\times
k}$ and
$\bm{\Delta}:=\sigma^{2}{\bm{I}}_{r}+\bm{\Sigma}{\bm{V}}^{T}{\bm{P}}{\bm{V}}\bm{\Sigma}\succ\bm{0}_{r\times
r}$ with the entries of $\widetilde{\bm{\pi}}$ given by
$\widetilde{\pi}_{\ell}=\frac{1}{\gamma}\frac{\pi_{\ell}\omega_{\ell}^{2}}{\omega_{\ell}^{2}+\eta}$.
Then, for the WLS linear classifier
$\left(\widehat{\bm{W}},\widehat{\bm{b}}\right)$ the following asymptotic
limits hold

$\displaystyle\widehat{\bm{b}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\widetilde{\bm{\pi}}-{\bm{P}}{\bm{V}}\bm{\Sigma}\bm{\Delta}^{-1}\bm{\Sigma}{\bm{V}}^{T}\widetilde{\bm{\pi}}\,,\quad\bm{\Sigma}_{\bm{w},\bm{\mu}}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}{\bm{P}}{\bm{V}}\bm{\Sigma}\bm{\Delta}^{-1}\bm{\Sigma}{\bm{V}}^{T}\,,$
(B.1a) $\displaystyle\bm{\Sigma}_{\bm{w},\bm{w}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\frac{\zeta}{\sigma^{2}}{\bm{P}}+{\bm{P}}{\bm{V}}\bm{\Sigma}\bm{\Delta}^{-1}\Big{(}\bm{\Delta}^{-1}-\frac{\zeta}{\sigma^{2}}{\bm{I}}_{r}\Big{)}\bm{\Sigma}{\bm{V}}^{T}{\bm{P}}+\frac{\eta\zeta}{\sigma^{2}}\bm{Q}\,.$
(B.1b)

Here,
$\zeta:={\gamma}\big{/}{\left(\eta\,\sum_{\ell=1}^{k}\frac{\pi_{\ell}\omega_{\ell}^{2}}{(\omega_{\ell}^{2}+\eta)^{2}}\right)}$
and $\bm{Q}\in\mathbb{R}^{k\times k}$ is a known matrix depending on various
problem parameters. Its precise value is given in (I.2).

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
$\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.

###### Theorem B.2

Consider data generated according to MLM and $\gamma<1$. Consider a weighted
LS classifier with weights $\bm{D}=\text{diag}(\omega_{1},\ldots,\omega_{k})$
and let $\eta$ be the unique solution to
$\sum_{\ell=1}^{k}\frac{\pi_{\ell}\omega_{\ell}^{2}}{\omega_{\ell}^{2}+\eta}=\gamma$.
Also define vector $\bm{\nu}\in\mathbb{R}^{k}$ with entries given by
${\bm{\nu}}_{\ell}=\frac{1}{\gamma}\frac{\omega_{\ell}^{2}}{\omega_{\ell}^{2}+\eta}$
and matrix

$\displaystyle\bm{\Delta}=\operatorname{\mathbb{E}}\left[\left(\bm{\nu}^{T}\bm{v}\right){\bm{g}}{\bm{g}}^{T}\right]-\bm{\Sigma}{\bm{V}}^{T}\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right)\bm{\nu}\bm{\nu}^{T}\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right){\bm{V}}\bm{\Sigma}\succ\bm{0}_{r\times
r},$ (B.2)

where $\bm{v}\in\mathbb{R}^{k}$ is a random vector with entries
$V_{\ell}={e^{\bm{e}_{\ell}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}}\big{/}{\sum_{\ell^{\prime}\in[k]}e^{\bm{e}_{\ell^{\prime}}{\bm{V}}\bm{\Sigma}{\bm{g}}}}$
for ${\bm{g}}\sim\mathcal{N}(\bm{0},{\bm{I}}_{r})$. Then, for the WLS linear
classifier $\left(\widehat{\bm{W}},\widehat{\bm{b}}\right)$ the following
asymptotic limits hold

$\displaystyle\widehat{\bm{b}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\text{diag}(\bm{\nu})\bm{\pi}-\text{diag}(\bm{\nu})\left({\bm{I}}_{k}-\bm{\pi}\bm{\nu}^{T}\right)\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right){\bm{V}}\bm{\Sigma}\bm{\Delta}^{-1}\bm{\Sigma}{\bm{V}}^{T}\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right)\bm{\nu}.$
(B.3a) $\displaystyle\bm{\Sigma}_{\bm{w},\bm{\mu}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\text{diag}(\bm{\nu})\left({\bm{I}}_{k}-\bm{\pi}\bm{\nu}^{T}\right)\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right){\bm{V}}\bm{\Sigma}\bm{\Delta}^{-1}\bm{\Sigma}{\bm{V}}^{T}.$
(B.3b)

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

Of course, the theorem above includes Theorem 4.2 as a special case. Indeed,
we show how setting $\omega_{\ell}=1,~{}\ell\in[k]$ recovers the solution for
(un-weighted) LS. First, solving for $\eta$ simply gives
$\eta=\frac{1}{\gamma}-1.$ Thus, $\bm{\nu}=\bm{1}_{k}$. Also, observe in (4.1)
that $\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right)\bm{1}_{k}=\bm{0}$ and
$\bm{1}^{T}\bm{v}=1$. Thus, (B.2) reduces to
$\bm{\Delta}=\operatorname{\mathbb{E}}[{\bm{g}}{\bm{g}}^{T}]={\bm{I}}_{r}$.
With these, it can be readily checked that (B.3a) and (B.3b) simplify to the
expressions in (4.4a).

The term
$\operatorname{\mathbb{E}}\left[\left(\bm{\nu}^{T}\bm{v}\right){\bm{g}}{\bm{g}}^{T}\right]$
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
$\Sigma_{\bm{w},\bm{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.

##  Appendix C Preliminaries

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

###  C.1 Slepian’s inequality

######  Lemma C.1 (Slepian’s inequality [LT91])

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

$\bm{S}_{ii}={\bm{R}}_{ii},\qquad\text{and}\qquad\bm{S}_{ij}\geq{\bm{R}}_{ij}.$

Then, for any $\bm{t}\in\mathbb{R}^{k}$ it holds that

$\operatorname{\mathbb{P}}\left\\{\bigcup_{j\in[k]}\left\\{{\bm{g}}_{j}\geq\bm{t}_{j}\right\\}\right\\}\leq\operatorname{\mathbb{P}}\left\\{\bigcup_{j\in[k]}\left\\{\widetilde{{\bm{g}}}_{j}\geq\bm{t}_{j}\right\\}\right\\}\,.$

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

$1-\operatorname{\mathbb{P}}\left\\{\bm{S}^{1/2}{\bm{z}}\leq\bm{t}\right\\}\leq
1-\operatorname{\mathbb{P}}\left\\{{\bm{R}}^{1/2}{\bm{z}}\leq\bm{t}\right\\}.$

###  C.2 Gaussian integration by parts

We say that a function $F:\mathbb{R}^{m}\rightarrow\mathbb{R}$ is of moderate
growth if for each $c>0$,

$\lim_{\left\|\bm{x}\right\|_{\ell_{2}}\rightarrow\infty}~{}F(\bm{x})\exp\left(-c\left\|\bm{x}\right\|_{\ell_{2}}^{2}\right)=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 ${\bm{g}}\sim\mathcal{N}(\bm{0},{\bm{I}}_{r})$ and function
$f:\mathbb{R}^{r}\rightarrow\mathbb{R}$ such that $f$ and all its first and
second order partial derivatives are of moderate growth. Then, the following
statements are true:

(i)
$\operatorname{\mathbb{E}}\left[f({\bm{g}}){\bm{g}}\right]=\operatorname{\mathbb{E}}\left[\nabla
f({\bm{g}})\right].$

(ii)
$\operatorname{\mathbb{E}}\left[f({\bm{g}}){\bm{g}}{\bm{g}}^{T}\right]=\operatorname{\mathbb{E}}\left[f({\bm{g}})\right]{\bm{I}}_{r}+\operatorname{\mathbb{E}}\left[\nabla^{2}f({\bm{g}})\right].$

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

######  Lemma C.3 (GIP for the Softmax)

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

$\displaystyle\bm{v}=\frac{e^{{\bm{V}}\bm{\Sigma}{\bm{g}}}}{\bm{1}_{k}^{T}e^{{\bm{V}}\bm{\Sigma}{\bm{g}}}},~{}~{}~{}~{}~{}~{}V_{\ell}=\frac{e^{\bm{e}_{\ell}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}}{\sum_{j\in[k]}e^{\bm{e}_{j}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}},~{}\ell\in[k].$
(C.1)

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

(i) $\operatorname{\mathbb{E}}[\bm{v}]=\bm{\pi}.$

(ii) For all $i\in[r],\ell\in[k]$,

$\operatorname{\mathbb{E}}\left[{\bm{g}}_{i}V_{\ell}\right]=(\bm{e}_{i}^{T}\bm{\Sigma}{\bm{V}}^{T}\bm{e}_{\ell})\,\bm{\pi}_{\ell}-\bm{e}_{i}^{T}\bm{\Sigma}{\bm{V}}^{T}\,\sum_{j\in[k]}{\bm{e}_{j}\bm{\Pi}_{ij}},$

and in matrix form:

$\operatorname{\mathbb{E}}\left[\bm{v}{\bm{g}}^{T}\right]=\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right){\bm{V}}\bm{\Sigma}\,.$

(iii) For all $\ell\in[k]$ let $s_{\ell}:\mathbb{R}^{k}\rightarrow\mathbb{R}$
denote the soft-max function:
$s_{\ell}(\bm{x})=\frac{e^{x_{\ell}}}{\sum_{i\in[k]}e^{x_{i}}}.$ Then,

$\operatorname{\mathbb{E}}\left[V_{\ell}{\bm{g}}{\bm{g}}^{T}\right]=\pi_{\ell}{\bm{I}}_{r}+\bm{\Sigma}^{T}{\bm{V}}^{T}\operatorname{\mathbb{E}}\left[\nabla^{2}s_{\ell}({\bm{V}}\bm{\Sigma}{\bm{g}})\right]{\bm{V}}\bm{\Sigma}.$

###  C.3 Block matrix inversion

######  Lemma C.4 (Block matrix inversion)

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

$\displaystyle\bm{T}^{-1}\begin{bmatrix}\mathbf{f}\\\
\epsilon\end{bmatrix}=\begin{bmatrix}\bm{\Delta}^{-1}\left(\mathbf{f}-\frac{\epsilon}{\delta}\bm{b}\right)\\\
\frac{\epsilon}{\delta}-\frac{1}{\delta}\bm{b}^{T}\bm{\Delta}^{-1}\left(\mathbf{f}-\frac{\epsilon}{\delta}\bm{b}\right)\end{bmatrix}$
(C.2)

where $\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

###  D.1 Proof of (2.5) and (2.6)

GMM. Starting from (2.4) and using the fact that
$\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

$\displaystyle\operatorname{\mathbb{P}}_{e}=\operatorname{\mathbb{P}}\Big{\\{}\arg\max_{j\in[k]}\\{\langle\widehat{\bm{w}}_{j},{\bm{M}}\bm{e}_{Y}\rangle+\langle\widehat{\bm{w}}_{j},{\bm{z}}\rangle+\widehat{\bm{b}}_{j}\\}\neq
Y\\}\Big{\\}},$

or, in matrix-form:

$\displaystyle\operatorname{\mathbb{P}}_{e}=\operatorname{\mathbb{P}}\Big{\\{}\arg\max\\{\widehat{\bm{W}}{\bm{M}}\bm{e}_{Y}+\widehat{\bm{W}}{{\bm{z}}}+\widehat{\bm{b}}\\}\neq
Y\\}\Big{\\}}.$

Recall that $\bm{\Sigma}_{\bm{w},\bm{\mu}}:=\widehat{\bm{W}}{\bm{M}}$ and note
that $\widehat{\bm{W}}{{\bm{z}}}$ is a zero-mean Gaussian vector with
covariance matrix
$\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).

MLM. Recall from (2.2) that $\bm{x}\sim\mathcal{N}(\bm{0},{\bm{I}}_{d})$ and
$Y$ is distributed such that
$\operatorname{\mathbb{P}}\\{Y=\ell~{}|~{}\bm{x}\\}={e^{\langle\bm{\mu}_{\ell},\bm{x}\rangle}}\big{/}{\sum_{j\in[k]}e^{\langle\bm{\mu}_{j},\bm{x}\rangle}}.$
Let ${\bm{g}}=\widehat{\bm{W}}\bm{x}$ and $\bm{h}={\bm{M}}^{T}\bm{x}$. In this
notation, (2.4) becomes

$\operatorname{\mathbb{P}}_{e}=\operatorname{\mathbb{P}}\left\\{\arg\max\,({\bm{g}}+\widehat{\bm{b}}\,)\neq
Y\right\\},$

with
$\operatorname{\mathbb{P}}\\{Y=\ell~{}|~{}\bm{x}\\}={e^{\bm{h}_{\ell}}}\big{/}{\sum_{j\in[k]}e^{\bm{h}_{j}}}$.
To complete the proof of (2.6), it is easy to check that
$\begin{bmatrix}{\bm{g}}\\\ \bm{h}\end{bmatrix}$ defined above is jointly
Gaussian with zero-mean and covariance matrix
$\begin{bmatrix}\bm{\Sigma}_{\bm{w},\bm{w}}&\bm{\Sigma}_{\bm{w},\bm{\mu}}\\\
\bm{\Sigma}_{\bm{w},\bm{\mu}}^{T}&\bm{\Sigma}_{\bm{\mu},\bm{\mu}}\end{bmatrix}\,.$

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

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

$\displaystyle\operatorname{\mathbb{P}}_{e|c}$
$\displaystyle=\operatorname{\mathbb{P}}\left(\exists j\neq
c~{}:~{}\langle\widehat{\bm{w}}_{c}-\widehat{\bm{w}}_{j},{\bm{z}}\rangle\leq\langle\widehat{\bm{w}}_{j}-\widehat{\bm{w}}_{c},\bm{\mu}_{c}\rangle+(\widehat{\bm{b}}_{j}-\widehat{\bm{b}}_{c})\right)$
(D.1) $\displaystyle=1-\operatorname{\mathbb{P}}\left(\forall j\neq
c~{}:~{}\langle\widehat{\bm{w}}_{c}-\widehat{\bm{w}}_{j},{\bm{z}}\rangle\geq\langle\widehat{\bm{w}}_{j}-\widehat{\bm{w}}_{c},\bm{\mu}_{c}\rangle+(\widehat{\bm{b}}_{j}-\widehat{\bm{b}}_{c})\right),$
(D.2)

where we used that $\bm{x}=\bm{\mu}_{c}+{\bm{z}}$. Let
$\bm{S}_{c}\in\mathbb{R}^{{(k-1)}\times{(k-1)}}$ be a symmetric matrix and
$\bm{t}_{c}\in\mathbb{R}^{k-1}$ a vector with entries:

$\displaystyle[\bm{t}_{c}]_{j}$
$\displaystyle:=\langle\widehat{\bm{w}}_{j}-\widehat{\bm{w}}_{c},\bm{\mu}_{c}\rangle+(\widehat{\bm{b}}_{j}-\widehat{\bm{b}}_{c})~{}~{}j\neq
c\in[k]$ (D.3a) $\displaystyle[\bm{S}_{c}]_{ij}$
$\displaystyle:=\langle\widehat{\bm{w}}_{c}-\widehat{\bm{w}}_{j},\widehat{\bm{w}}_{c}-\widehat{\bm{w}}_{i}\rangle,~{}~{}i,j\neq
c\in[k].$ (D.3b)

Then, we can rewrite (D.2) as

$\displaystyle\operatorname{\mathbb{P}}_{e|c}:=1-\operatorname{\mathbb{P}}\Big{\\{}\bm{S}_{c}^{{1}/{2}}\,{\bm{z}}\geq\bm{t}_{c}\Big{\\}},$
(D.4)

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

Further, by using the law of total probability we have

$\displaystyle\operatorname{\mathbb{P}}_{e}=\sum_{c=1}^{k}\pi_{c}\operatorname{\mathbb{P}}_{e|c}=\sum_{c=1}^{k}\pi_{c}\left(1-\operatorname{\mathbb{P}}\Big{\\{}\bm{S}_{c}^{{1}/{2}}\,{\bm{z}}\geq\bm{t}_{c}\Big{\\}}\right).$

###  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):

$\displaystyle\operatorname{\mathbb{P}}_{e}=\operatorname{\mathbb{P}}\big{\\{}\arg\max\text{
}(\,{\bm{g}}+\widehat{\bm{b}}\,)\neq
Y(\bm{h})\big{\\}}\,,\quad\text{where}~{}~{}\begin{bmatrix}{\bm{g}}\\\
\bm{h}\end{bmatrix}\sim\mathcal{N}\Big{(}\bm{0},\begin{bmatrix}\bm{\Sigma}_{\bm{w},\bm{w}}&\bm{\Sigma}_{\bm{w},\bm{\mu}}\\\
\bm{\Sigma}_{\bm{w},\bm{\mu}}^{T}&\bm{\Sigma}_{\bm{\mu},\bm{\mu}}\end{bmatrix}\Big{)}\,,$

and
$\operatorname{\mathbb{P}}\\{Y(\bm{h})=\ell\\}={e^{\bm{h}_{\ell}}}/{\sum_{j\in[k]}e^{\bm{h}_{j}}},~{}~{}\ell\in[k].$
Using Gaussian decomposition we can write
${\bm{g}}=\widetilde{{\bm{g}}}+\bm{\Sigma}_{\bm{w},\bm{\mu}}\bm{\Sigma}_{\bm{\mu},\bm{\mu}}^{\dagger}\bm{h}$
where
$\widetilde{{\bm{g}}}\sim\mathcal{N}(\bm{0}_{k},\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}})$.
Using this, we have

$\displaystyle\operatorname{\mathbb{P}}_{e}$
$\displaystyle=\operatorname{\mathbb{P}}\big{\\{}\arg\max\text{
}(\widetilde{{\bm{g}}}+\bm{\Sigma}_{\bm{w},\bm{\mu}}\bm{\Sigma}_{\bm{\mu},\bm{\mu}}^{\dagger}\bm{h}+\widehat{\bm{b}})\neq
Y(\bm{h})\big{\\}}$
$\displaystyle=\operatorname{\mathbb{P}}\big{\\{}\arg\max\text{
}(\widetilde{{\bm{g}}}+\bm{\Sigma}_{\bm{w},\bm{\mu}}\bm{\Sigma}_{\bm{\mu},\bm{\mu}}^{\dagger}\bm{h}+\widehat{\bm{b}})\neq
Y(\bm{h})\big{\\}}$
$\displaystyle=\sum_{c\in[k]}\operatorname{\mathbb{E}}_{\bm{h},\widetilde{{\bm{g}}}}\left[\operatorname{\mathbb{P}}\left\\{Y(\bm{h})=c\right\\}\cdot\mathds{1}\left\\{\arg\max\text{
}(\widetilde{{\bm{g}}}+\bm{\Sigma}_{\bm{w},\bm{\mu}}\bm{\Sigma}_{\bm{\mu},\bm{\mu}}^{\dagger}\bm{h}+\widehat{\bm{b}})\neq
c\right\\}\right]$
$\displaystyle=\sum_{c\in[k]}\operatorname{\mathbb{E}}_{\bm{h},\widetilde{{\bm{g}}}}\left[\frac{e^{\bm{h}_{c}}}{{\sum_{\ell\in[k]}e^{\bm{h}_{\ell}}}}\left(1-\prod_{j\neq
c}\mathds{1}\left\\{\widetilde{{\bm{g}}}_{c}+[\bm{\Sigma}_{\bm{w},\bm{\mu}}\bm{\Sigma}_{\bm{\mu},\bm{\mu}}^{\dagger}\bm{h}]_{c}+\widehat{\bm{b}}_{c}\geq\widetilde{{\bm{g}}}_{j}+[\bm{\Sigma}_{\bm{w},\bm{\mu}}\bm{\Sigma}_{\bm{\mu},\bm{\mu}}^{\dagger}\bm{h}]_{j}+\widehat{\bm{b}}_{j}\right\\}\right)\right]$
$\displaystyle=\sum_{c\in[k]}\operatorname{\mathbb{E}}_{\bm{h}}\left[\frac{e^{\bm{h}_{c}}}{{\sum_{\ell\in[k]}e^{\bm{h}_{\ell}}}}\left(1-\operatorname{\mathbb{P}}_{{\bm{z}}\sim\mathcal{N}(\bm{0},{\bm{I}}_{k-1})}\big{\\{}\bm{S}_{c}^{1/2}{\bm{z}}\geq\bm{t}_{c}(\bm{h})\big{\\}}\right)\right],$
(D.5)

where in the last line

$\displaystyle[\bm{t}_{c}(\bm{h})]_{j}$
$\displaystyle=\widehat{\bm{b}}_{j}-\widehat{\bm{b}}_{c}+[\bm{\Sigma}_{\bm{w},\bm{\mu}}\bm{\Sigma}_{\bm{\mu},\bm{\mu}}^{\dagger}\bm{h}]_{j}-[\bm{\Sigma}_{\bm{w},\bm{\mu}}\bm{\Sigma}_{\bm{\mu},\bm{\mu}}^{\dagger}\bm{h}]_{c},~{}~{}j\neq
c\in[k]$ (D.6a) $\displaystyle[\bm{S}_{c}]_{i,j}$
$\displaystyle=\left(\bm{e}_{c}-\bm{e}_{j}\right)^{T}\left(\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}}\right)\left(\bm{e}_{c}-\bm{e}_{i}\right),~{}~{}i,j\neq
c\in[k].$ (D.6b)

Further recalling the decomposition
$\operatorname{\mathbb{P}}_{e}=\sum_{c}\mathbb{P}\Big{\\{}\widehat{Y}\neq
Y|Y=c\Big{\\}}\operatorname{\mathbb{P}}\left\\{Y=c\right\\}$ and noting that

$\displaystyle\mathbb{P}\\{Y=c\\}=\operatorname{\mathbb{E}}_{\bm{x}}\Big{[}\mathbb{P}\Big{\\{}Y=c|\bm{x}\Big{\\}}\Big{]}=\operatorname{\mathbb{E}}_{\bm{x}}\Bigg{[}\frac{1}{\left(1+\sum_{j\neq
c}e^{(\bm{\mu}_{j}-\bm{\mu}_{c})^{T}\bm{x}}\right)}\Bigg{]}=\bm{\pi}_{c}$

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

$\displaystyle\operatorname{\mathbb{P}}_{e|c}=\mathbb{P}\Big{\\{}\widehat{Y}\neq
Y|Y=c\Big{\\}}=\frac{1}{\bm{\pi}_{c}}\,\operatorname{\mathbb{E}}_{\bm{h}\sim\mathcal{N}(\bm{0}_{k},\bm{\Sigma}_{\bm{\mu},\bm{\mu}})}\left[\frac{e^{\bm{h}_{c}}}{{\sum_{\ell\in[k]}e^{\bm{h}_{\ell}}}}\left(1-\operatorname{\mathbb{P}}_{{\bm{z}}\sim\mathcal{N}(\bm{0},{\bm{I}}_{k-1})}\left\\{\bm{S}_{c}^{1/2}{\bm{z}}\geq\bm{t}_{c}(\bm{h})\right\\}\right)\right],$
(D.7)

where $\bm{\pi}_{c}$ is the $c^{\text{th}}$ entry of the vector $\bm{\pi}$ in
(4.1) and $\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 ${\bm{z}}\sim\mathcal{N}(\bm{0},{\bm{I}}_{k-1})$,

$\displaystyle
1-\operatorname{\mathbb{P}}\\{{\bm{A}}^{1/2}{\bm{z}}\leq\bm{t}\\}=\operatorname{\mathbb{P}}\left\\{\bigcup_{i\in[k-1]}\\{\widetilde{\bm{a}}_{i}^{T}{\bm{z}}\geq\bm{t}_{i}\\}\right\\}\,.$
(D.8)

Here,
${\bm{A}}\succeq\bm{0}\in\mathbb{R}^{(k-1)\times(k-1)},\bm{t}\in\mathbb{R}^{k-1}$
are appropriate coefficient matrices (see (D.4) and (D.5)) and
$\widetilde{\bm{a}}_{i}$ denotes the $i$th row of the matrix ${\bm{A}}^{1/2}$.
For example, (D.8) maps to (D.4) for ${\bm{A}}\leftarrow\bm{S}_{c}$ and
$\bm{t}\leftarrow(-\bm{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 ${\bm{A}}$ and vector
$\bm{t}$ in (D.8) take the special form
${\bm{A}}\propto{\bm{I}}+\bm{1}\bm{1}^{T}$ and $\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.

###### Lemma D.1

Let ${\bm{A}}={\bm{I}}_{k}+\bm{1}_{k}\bm{1}_{k}^{T}$ and
${\bm{g}}\sim\mathcal{N}(\bm{0},{\bm{A}})$. Then, for any $t\in\mathbb{R}$,

$\displaystyle 1-\operatorname{\mathbb{P}}\\{{\bm{g}}\leq
t\bm{1}_{k}\\}=\operatorname{\mathbb{P}}\\{G_{0}+\max_{i\in[k]}G_{i}\geq
t\\},~{}~{}G_{0},G_{1},\ldots,G_{k}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}(0,1)\,.$ (D.9)

Proof For each $i\in[n]$, we can decompose ${\bm{g}}_{i}=G_{0}+G_{i},$ where
$G_{0},G_{1},\ldots,G_{k}$ are iid standard normals. Indeed, it can be readily
checked from this that $\operatorname{\mathbb{E}}[{\bm{g}}_{i}^{2}]=2$ and
$\operatorname{\mathbb{E}}[{\bm{g}}_{i}{\bm{g}}_{j}]=\operatorname{\mathbb{E}}[G_{0}^{2}]=1,~{}i\neq
j$, which is consistent with ${\bm{g}}\sim\mathcal{N}(\bm{0},{\bm{A}})$. Thus,
we can write

$1-\operatorname{\mathbb{P}}\\{{\bm{g}}\leq
t\bm{1}_{k}\\}=\operatorname{\mathbb{P}}\\{\max_{i\in[k]}{\bm{g}}_{i}\geq
t\\}=\operatorname{\mathbb{P}}\\{G_{0}+\max_{i\in[k]}G_{i}\geq t\\},$

which completes the proof.

####  D.4.2 Slepian’s bound

When the matrix ${\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\in[k]}{\bm{A}}_{ij}\geq 0.$ To begin,
note that ${\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:

$\displaystyle
1-\operatorname{\mathbb{P}}\\{{\bm{A}}^{1/2}{\bm{g}}\leq\bm{t}\\}$
$\displaystyle\leq 1-\operatorname{\mathbb{P}}\left\\{\left(\left({\rm
diag}({\bm{A}})-a{\bm{I}}\right)\,+a\bm{1}\bm{1}^{T}\right)^{1/2}{\bm{g}}\leq\bm{t}\right\\}$
$\displaystyle\leq\operatorname{\mathbb{P}}\Big{\\{}\bigcup_{j\in[k]}\big{\\{}\,G_{0}+G_{j}\sqrt{{[{\bm{A}}]_{jj}}/{a}-1}\geq{[\bm{t}]_{j}}/{a}\,\big{\\}}\Big{\\}},\qquad
G_{0},G_{1},\ldots,G_{k}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}(0,1)\,.$

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:

$\displaystyle\operatorname{\mathbb{P}}_{e|c}$
$\displaystyle=\operatorname{\mathbb{P}}\left(\exists j\neq
c~{}:~{}\langle\widehat{\bm{w}}_{c}-\widehat{\bm{w}}_{j},{\bm{z}}\rangle\leq\langle\widehat{\bm{w}}_{j}-\widehat{\bm{w}}_{c},\bm{\mu}_{c}\rangle+(\widehat{\bm{b}}_{j}-\widehat{\bm{b}}_{c})\right)$
(D.10) $\displaystyle\leq\sum_{j\neq
c}\operatorname{\mathbb{P}}_{G\sim\mathcal{N}(0,1)}\left\\{\left\|\widehat{\bm{w}}_{j}-\widehat{\bm{w}}_{c}\right\|_{\ell_{2}}G\leq\langle\widehat{\bm{w}}_{j}-\widehat{\bm{w}}_{c},\bm{\mu}_{c}\rangle+(\widehat{\bm{b}}_{j}-\widehat{\bm{b}}_{c})\right\\}$
$\displaystyle=\sum_{j\neq
c}Q\left(\langle\frac{\widehat{\bm{w}}_{c}-\widehat{\bm{w}}_{j}}{\left\|\widehat{\bm{w}}_{c}-\widehat{\bm{w}}_{j}\right\|_{\ell_{2}}},\bm{\mu}_{c}\rangle+\frac{\widehat{\bm{b}}_{c}-\widehat{\bm{b}}_{j}}{\left\|\widehat{\bm{w}}_{c}-\widehat{\bm{w}}_{j}\right\|_{\ell_{2}}}\right)$
$\displaystyle=\sum_{j\neq
c}Q\left(\frac{-[\bm{t}_{c}]_{j}}{\sqrt{[\bm{S}_{c}]_{jj}}}\right)$ (D.11)
$\displaystyle\leq(k-1)\cdot Q\left(d_{\min}\right),$

where in (D.11) $\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{\\}}.$

Union bound without knowledge of cross-correlations
$\langle\bm{w}_{i},\bm{w}_{j}\rangle,~{}i\neq j$. It is worth noting that the
upper bound in (D.11) requires knowledge of the cross-correlations
$\langle\widehat{\bm{w}}_{j},\widehat{\bm{w}}_{c}\rangle,~{}j\neq c$, i.e., of
the off-diagonal entries of $\bm{\Sigma}_{\bm{w},\bm{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 $\bm{\Sigma}_{\bm{w},\bm{w}}$ (i.e., the norms
$\left\|\widehat{\bm{w}}_{j}\right\|_{\ell_{2}},j\in[k]$) were known. A simple
answer to this questions is as follows. Observe that
$\sqrt{[\bm{S}_{c}]_{jj}}=\left\|\widehat{\bm{w}}_{c}-\widehat{\bm{w}}_{j}\right\|_{\ell_{2}}\leq\left\|\widehat{\bm{w}}_{c}\right\|_{\ell_{2}}+\left\|\widehat{\bm{w}}_{j}\right\|_{\ell_{2}}.$
Thus, if $[\bm{t}_{c}]_{j}<0$ then the $j$th term in (D.11) is further upper
bounded by
$Q\left(\frac{-[\bm{t}_{c}]_{j}}{\left\|\widehat{\bm{w}}_{c}\right\|_{\ell_{2}}+\left\|\widehat{\bm{w}}_{j}\right\|_{\ell_{2}}}\right)$:

$\operatorname{\mathbb{P}}_{e|c}\leq\begin{cases}\sum_{j\neq
c}Q\left(\frac{-[\bm{t}_{c}]_{j}}{\left\|\widehat{\bm{w}}_{c}\right\|_{\ell_{2}}+\left\|\widehat{\bm{w}}_{j}\right\|_{\ell_{2}}}\right)&\text{if
}\bm{t}_{c}\leq\bm{0},\\\ 1&\text{otherwise.}\end{cases}$

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

$[\bm{t}_{c}]_{j}<0\Leftrightarrow\langle\widehat{\bm{w}}_{c},\bm{\mu}_{c}\rangle+\widehat{\bm{b}}_{c}\geq\langle\widehat{\bm{w}}_{j},\bm{\mu}_{c}\rangle+\widehat{\bm{b}}_{j}\text{
for all }j\neq c.$

Intuitively, this assumes a regime wherethe weight vector
$\widehat{\bm{w}}_{c}$ corresponding to class $c$ aligns better with the
corresponding mean vector $\bm{\mu}_{c}$ than the rest of the weight vectors
$\widehat{\bm{w}}_{j},~{}j\neq c.$ This emphasizes the important role of the
cross-correlation matrix $\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.

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 $\bm{\mu}_{i},~{}i\in[n]$ are known. Then the problem of classifying a
new sample $\bm{x}$ is a $k$-ary hypothesis testing problem with Gaussian
conditionals. Denote $\operatorname{\mathbb{P}}_{\rm genie,Bayes}$ the Bayes
error of this hypothesis testing problem. Clearly
$\operatorname{\mathbb{P}}_{\rm genie,Bayes}$ is a lower bound on the error of
any classifier that is trained on data. The paper [WBW+16] further lower
bounds $\operatorname{\mathbb{P}}_{\rm genie,Bayes}$ in terms of the Bayesian
probability of errors between every two classes as follows:

$\displaystyle\operatorname{\mathbb{P}}_{e}$
$\displaystyle\geq\operatorname{\mathbb{P}}_{\rm
genie,Bayes}\geq\frac{2}{k}\sum_{i\neq j}\pi_{i}\operatorname{\mathbb{P}}_{\rm
genie,ij}$ (D.12) $\displaystyle=\frac{2}{k}\sum_{i\neq
j}\frac{\pi_{i}}{\pi_{i}+\pi_{j}}\Big{\\{}\pi_{i}\cdot
Q\left(\frac{\left\|\bm{\mu}_{i}-\bm{\mu}_{j}\right\|_{\ell_{2}}}{2}+\frac{\log(\pi_{i}/\pi_{j})}{\left\|\bm{\mu}_{i}-\bm{\mu}_{j}\right\|_{\ell_{2}}}\right)+\pi_{j}\cdot
Q\left(\frac{\left\|\bm{\mu}_{i}-\bm{\mu}_{j}\right\|_{\ell_{2}}}{2}-\frac{\log(\pi_{i}/\pi_{j})}{\left\|\bm{\mu}_{i}-\bm{\mu}_{j}\right\|_{\ell_{2}}}\right)\Big{\\}}.$

where $\operatorname{\mathbb{P}}_{\rm genie,ij}$ is the Bayesian error between
classes $i$ and $j$ with priors $\frac{\pi_{i}}{\pi_{i}+\pi_{j}}$ and
$\frac{\pi_{j}}{\pi_{i}+\pi_{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

$\displaystyle\operatorname{\mathbb{P}}_{e}\geq\frac{2}{k}\sum_{i\neq
j}{\pi_{i}}\cdot
Q\left(\frac{\left\|\bm{\mu}_{i}-\bm{\mu}_{j}\right\|_{\ell_{2}}}{2}\right).$
(D.13)

Note that, in contrast to (D.13), our analysis allows for precise evaluations
of the missclassification error $\operatorname{\mathbb{P}}_{e}$.

##  Appendix E The Class-averaging estimator

###  E.1 Proofs for GMM

####  E.1.1 GMM: Proof of Proposition 3.1

The first statement (3.1a) follows directly from the fact that
$\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

$\displaystyle\widehat{\bm{w}}_{i}=\frac{1}{n}{\bm{M}}{\bm{Y}}{\bm{Y}}_{i}+\frac{1}{n}\bm{Z}{\bm{Y}}_{i}=\frac{1}{n}\sum_{j=1}^{k}\bm{\mu}_{j}({\bm{Y}}_{j}^{T}{\bm{Y}}_{i})+\frac{1}{n}\bm{Z}{\bm{Y}}_{i}=\frac{\left\|{\bm{Y}}_{i}\right\|_{\ell_{2}}^{2}}{n}\bm{\mu}_{i}+\frac{1}{n}\bm{Z}{\bm{Y}}_{i},$
(E.1)

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

$\displaystyle\langle{\bm{Y}}_{i},{\bm{Y}}_{j}\rangle=0,\forall i\neq
j\in[k].$ (E.2)

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

  1. (i)

$\frac{\left\|{\bm{Y}}_{i}\right\|_{\ell_{2}}^{2}}{n}=\frac{n_{i}}{n}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\pi_{i}$.

  2. (ii)

$\bm{Z}{\bm{Y}}_{i}\sim\sigma\left\|{\bm{Y}}_{i}\right\|_{\ell_{2}}{\bm{g}}_{i}$
with ${\bm{g}}_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}(\bm{0},{\bm{I}}_{d})$ because of (E.2).

  3. (iii)

$\frac{\left\|{\bm{g}}_{i}\right\|_{\ell_{2}}}{\sqrt{n}}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\sqrt{\gamma}$
and
$\frac{1}{\sqrt{n}}\langle{\bm{g}}_{i},\bm{\mu}_{j}\rangle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}0$.

###  E.2 Proofs for MLM

####  E.2.1 Proof of Proposition 4.1

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

$\displaystyle\bm{v}=\frac{e^{{\bm{V}}\bm{\Sigma}{\bm{g}}}}{\bm{1}_{k}^{T}e^{{\bm{V}}\bm{\Sigma}{\bm{g}}}},~{}~{}~{}~{}~{}~{}V_{i}=\frac{e^{\bm{e}_{i}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}}{\sum_{j\in[k]}e^{\bm{e}_{j}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}},~{}i\in[k].$
(E.3)

We will prove the following three statements:

$\displaystyle\widehat{\bm{b}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\operatorname{\mathbb{E}}\left[\bm{v}\right]$
(E.4a) $\displaystyle\bm{\Sigma}_{\bm{w},\bm{\mu}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\operatorname{\mathbb{E}}\left[\bm{v}{\bm{g}}^{T}\right]\bm{\Sigma}{\bm{V}}^{T}$
(E.4b) $\displaystyle\bm{\Sigma}_{\bm{w},\bm{w}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\gamma\cdot\text{diag}(\operatorname{\mathbb{E}}\left[\bm{v}\right])+\operatorname{\mathbb{E}}\left[\bm{v}{\bm{g}}^{T}\right]\cdot\operatorname{\mathbb{E}}\left[{\bm{g}}\bm{v}^{T}\right].$
(E.4c)

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

For the intercepts $\widehat{\bm{b}}_{\ell},~{}\ell\in[k]$ it holds that

$\widehat{\bm{b}}_{\ell}=\frac{1}{n}\bm{1}^{T}{\bm{Y}}_{\ell}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\operatorname{\mathbb{P}}\\{Y=\ell\\}=\operatorname{\mathbb{E}}\left[\frac{e^{\bm{h}_{\ell}}}{\sum_{j\in[k]}e^{\bm{h}_{j}}}\right],$

where
$\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
$\bm{h}\stackrel{{\scriptstyle\rm(D)}}{{=\joinrel=}}{\bm{V}}\bm{\Sigma}{\bm{g}}$.

Continuing with the vectors $\widehat{\bm{w}}_{\ell},~{}\ell\in[k]$, recall
that
$\bm{w}_{\ell}=\frac{1}{n}{\bm{X}}{\bm{Y}}_{\ell}=\frac{1}{n}\sum_{i_{\ell}\in[n]}\bm{x}_{i_{\ell}}[{\bm{Y}}_{\ell}]_{i_{\ell}}$.
Consider the singular decomposition

$\bm{M}=\bm{U}\bm{\Sigma}\bm{V}^{T}=\begin{bmatrix}{\bm{u}}_{1}&{\bm{u}}_{2}&\ldots&{\bm{u}}_{r}\end{bmatrix}\text{diag}(\sigma_{1},\sigma_{2},\ldots,\sigma_{r})\begin{bmatrix}\bm{v}^{T}_{1}\\\
\bm{v}_{2}^{T}\\\ \ldots\\\ \bm{v}_{r}^{T}\end{bmatrix},$

with $\bm{U}\in\mathbb{R}^{d\times r}$, $\bm{\Sigma}\in\mathbb{R}^{r\times
r}$, and $\bm{V}\in\mathbb{R}^{k\times r}$ where $r=$rank$(\bm{M})\leq k$.
Decompose ${\bm{X}}\in\mathbb{R}^{d\times n}$ as
${\bm{X}}={\bm{U}}{\bm{U}}^{T}{\bm{X}}+{\bm{P}}^{\perp}{\bm{X}}$ with
${\bm{P}}^{\perp}={\bm{I}}_{d}-{\bm{U}}{\bm{U}}^{T}.$ With this notation we
compute

$\displaystyle\langle\bm{w}_{\ell},\bm{\mu}_{c}\rangle=\frac{1}{n}\sum_{i\in[n]}\bm{x}^{T}_{i}\bm{\mu}_{c}[{\bm{Y}}_{\ell}]_{i}$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{r}\sum_{i\in[n]}(\bm{x}^{T}_{i}{\bm{u}}_{j})\cdot(\bm{\mu}_{c}^{T}{\bm{u}}_{j})[{\bm{Y}}_{\ell}]_{i}+\frac{1}{n}\sum_{i\in[n]}\cdot(\bm{\mu}_{c}^{T}{\bm{P}}^{\perp}\bm{x}_{i})[{\bm{Y}}_{\ell}]_{i}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\sum_{j=1}^{r}\left(\bm{e}_{c}^{T}{\bm{V}}\bm{\Sigma}\bm{e}_{j}\right)\operatorname{\mathbb{E}}\left[{\bm{g}}_{j}\frac{\left(e^{\bm{e}_{\ell}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}\right)}{\sum_{\ell^{\prime}\in[k]}e^{\bm{e}_{\ell^{\prime}}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}}\right]$
$\displaystyle=\sum_{j=1}^{r}\operatorname{\mathbb{E}}[V_{\ell}{\bm{g}}_{j}]\,\left(\bm{e}_{j}^{T}\bm{\Sigma}{\bm{V}}^{T}\bm{e}_{c}\right).$
(E.5)

Here, we have recognized that for every
$i\in[n]:{\bm{U}}^{T}\bm{x}_{i}\sim{\bm{g}}$, and also, conditioned on
$\bm{x}_{i}$: $[{\bm{Y}}_{\ell}]_{i}\sim{\rm
Bern}\left(e^{\bm{\mu}_{\ell}^{T}\bm{x}_{i}}\big{/}\sum_{\ell^{\prime}}e^{\bm{\mu}_{\ell^{\prime}}^{T}\bm{x}_{i}}\right)$
and
$\bm{\mu}_{\ell}^{T}\bm{x}_{i}=\bm{e}_{\ell}^{T}{\bm{V}}\bm{\Sigma}{\bm{U}}^{T}\bm{x}_{i}\sim\bm{e}_{\ell}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}},\ell\in[k]$.
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:

$\displaystyle\langle\bm{w}_{\ell},\bm{w}_{c}\rangle=\frac{1}{n^{2}}\sum_{i_{\ell}\in[n],i_{c}\in[n]}\bm{x}^{T}_{i_{\ell}}\bm{x}_{i_{c}}[{\bm{Y}}_{\ell}]_{i_{\ell}}[{\bm{Y}}_{c}]_{i_{c}}$
$\displaystyle=\frac{1}{n^{2}}\sum_{j=1}^{r}\sum_{i_{\ell}\in[n],i_{c}\in[n]}(\bm{x}^{T}_{i_{\ell}}{\bm{u}}_{j})(\bm{x}^{T}_{i_{c}}{\bm{u}}_{j})[{\bm{Y}}_{\ell}]_{i_{\ell}}[{\bm{Y}}_{c}]_{i_{c}}+\frac{1}{n^{2}}\sum_{i_{\ell}\in[n],i_{c}\in[n]}\left({\bm{P}}^{\perp}\bm{x}_{i_{\ell}}\right)^{T}\left({\bm{P}}^{\perp}\bm{x}_{i_{c}}\right)[{\bm{Y}}_{\ell}]_{i_{\ell}}[{\bm{Y}}_{c}]_{i_{c}}$

For $i_{\ell}=i_{c}=i\in[n]$ note that

$\displaystyle\frac{1}{n^{2}}\sum_{i\in[n]}\sum_{j=1}^{r}(\bm{x}^{T}_{i}{\bm{u}}_{j})(\bm{x}^{T}_{i}{\bm{u}}_{j})[{\bm{Y}}_{\ell}]_{i}[{\bm{Y}}_{c}]_{i}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}0,$

while, for $i_{\ell}\neq i_{c}$,

$\displaystyle\frac{1}{n^{2}}\sum_{i_{\ell}\neq
i_{c}\in[n]}\sum_{j=1}^{r}(\bm{x}^{T}_{i_{\ell}}{\bm{u}}_{j})[{\bm{Y}}_{\ell}]_{i_{\ell}}(\bm{x}^{T}_{i_{c}}{\bm{u}}_{j})[{\bm{Y}}_{c}]_{i_{c}}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\sum_{j=1}^{r}\operatorname{\mathbb{E}}\left[{\bm{g}}_{j}\cdot\frac{\left(e^{\bm{e}_{\ell}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}\right)}{\sum_{\ell^{\prime}\in[k]}e^{\bm{e}_{\ell^{\prime}}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}}\right]\operatorname{\mathbb{E}}\left[{\bm{g}}_{j}\cdot\frac{\left(e^{\bm{e}_{c}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}\right)}{\sum_{\ell^{\prime}\in[k]}e^{\bm{e}_{\ell^{\prime}}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}}\right]$
$\displaystyle=\sum_{j=1}^{r}\operatorname{\mathbb{E}}[V_{\ell}{\bm{g}}_{j}]\operatorname{\mathbb{E}}[{\bm{g}}_{j}V_{c}]=\bm{e}_{\ell}^{T}{\bm{E}}[\bm{v}{\bm{g}}^{T}]\cdot{\bm{E}}[{\bm{g}}\bm{v}^{T}]\bm{e}_{c}^{T}$

Furthermore,

$\displaystyle\frac{1}{n}\sum_{i\in[n]}\left\|{\bm{P}}^{\perp}\bm{x}_{i}\right\|_{\ell_{2}}^{2}[{\bm{Y}}_{\ell}]_{i}[{\bm{Y}}_{c}]_{i}^{2}\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\gamma\cdot\mathds{1}_{\ell,c}\cdot\operatorname{\mathbb{E}}\left[\left(\frac{e^{\bm{e}_{\ell}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}}{\sum_{\ell^{\prime}\in[k]}e^{\bm{e}_{\ell^{\prime}}^{T}{\bm{V}}\bm{\Sigma}{\bm{g}}}}\right)\right]=\gamma\cdot\bm{e}_{\ell}^{T}\text{diag}(\operatorname{\mathbb{E}}[\bm{v}])\bm{e}_{c}.$

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

####  E.2.2 Orthogonal means

Here, we specialize the general result of Proposition 4.1 to the special case
of orthogonal means: $\langle\bm{\mu}_{i},\bm{\mu}_{j}\rangle=0,\forall i\neq
j$. Recall the notation
$\mu_{i}=\left\|\bm{\mu}_{i}\right\|_{\ell_{2}},i\in[k]$. Then, in this case
the parameters in (4.1) are simply given by the following

$\displaystyle\bm{\pi}_{i}:=\operatorname{\mathbb{E}}\big{[}\frac{e^{\mu_{i}G_{i}}}{\sum_{\ell\in[k]}e^{\mu_{\ell}G_{\ell}}}\big{]},i\in[k]\quad\text{and}\quad\bm{\Pi}_{ij}:=\operatorname{\mathbb{E}}\big{[}\frac{e^{\mu_{i}G_{i}}e^{\mu_{j}G_{j}}}{\left(\sum_{\ell\in[k]}e^{\mu_{\ell}G_{\ell}}\right)^{2}}\big{]},~{}i,j\in[k].$
(E.6)

Specifically, (4.3) can be equivalently expressed as

$\displaystyle\operatorname{\mathbb{P}}_{e,\rm Avg}$
$\displaystyle\stackrel{{\scriptstyle{P}}}{{\longrightarrow}}\operatorname{\mathbb{P}}\big{(}\arg\max_{\ell\in[k]}\left\\{\gamma\cdot\text{diag}(\bm{\pi})\cdot\widetilde{{\bm{g}}}+\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right)\cdot\bm{\Sigma}{\bm{g}}\right\\}\neq
Y({\bm{g}})\big{)}$ (E.7)
$\displaystyle=\operatorname{\mathbb{P}}\big{(}\bigcup_{j\neq
Y}\left\\{\gamma\cdot\bm{\pi}_{\ell}\cdot\widetilde{{\bm{g}}}_{\ell}\geq\gamma\cdot\bm{\pi}_{Y}\cdot\widetilde{{\bm{g}}}_{Y}+(\bm{e}_{Y}-\bm{e}_{\ell})^{T}\left(\text{diag}(\bm{\pi})-\bm{\Pi}\right)\bm{\Sigma}{\bm{g}}+(\bm{\pi}_{Y}-\bm{\pi}_{\ell})\right\\}\big{)},$
(E.8)

where ${\bm{g}},\widetilde{{\bm{g}}}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}(\bm{0},{\bm{I}}_{k})$,
$\operatorname{\mathbb{P}}\left(Y({\bm{g}})=c\right)=\frac{e^{\mu_{c}{\bm{g}}_{c}}}{\sum_{\ell\in[k]}e^{\mu_{\ell}{\bm{g}}_{\ell}}}$
and $\bm{\Sigma}=\text{diag}(\mu_{1},\ldots,\mu_{k})$.

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

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

$\bm{x}_{i}={\bm{M}}\bm{y}_{i}+{\bm{z}}_{i},\quad i\in[n],$

where the matrix of means ${\bm{M}}\in\mathbb{R}^{d\times k}$ has iid Gaussian
entries with variance $\mu^{2}/d$, ${\bm{z}}_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}(\bm{0},{\bm{I}}_{d})$ and
$\bm{y}_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}\rm{Unif}\left(\mathbf{e}_{1},\ldots,\mathbf{e}_{k}\right)$ with
$\mathbf{e}_{j}$ denoting the $j^{\text{t}h}$ canonical vector in
$\mathbb{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:

$\displaystyle\hat{\bm{y}}_{n+1}$
$\displaystyle=\arg\max_{\bm{e}_{j},~{}j\in[k]}P\left\\{\bm{y}=\bm{e}_{j}~{}|~{}{\bm{X}},{\bm{Y}},\bm{x}_{n+1}\right\\}$
$\displaystyle=\arg\max_{\bm{e}_{j},~{}j\in[k]}\int
P\left\\{\bm{y}=\bm{e}_{j}~{}|~{}{\bm{M}},{\bm{X}},{\bm{Y}},\bm{x}_{n+1}\right\\}P\left\\{{\bm{M}}~{}|~{}{\bm{X}},{\bm{Y}},\bm{x}_{n+1}\right\\}\,\mathrm{d}{\bm{M}}$
$\displaystyle=\arg\max_{\bm{e}_{j},~{}j\in[k]}\int\frac{P\left\\{\bm{x}_{n+1}~{}|~{}\bm{y}=\bm{e}_{j},{\bm{M}},{\bm{X}},{\bm{Y}}\right\\}\cdot
P\left\\{\bm{y}=\bm{e}_{j}~{}|~{}{\bm{M}},{\bm{X}},{\bm{Y}}\right\\}}{P\left\\{\bm{x}_{n+1}|{\bm{M}},{\bm{X}},{\bm{Y}}\right\\}}P\left\\{{\bm{M}}~{}|~{}{\bm{X}},{\bm{Y}},\bm{x}_{n+1}\right\\}\,\mathrm{d}{\bm{M}}$
$\displaystyle=\arg\max_{\bm{e}_{j},~{}j\in[k]}\int{P\left\\{\bm{x}_{n+1}~{}|~{}\bm{y}=\bm{e}_{j},{\bm{M}}\right\\}}\,\frac{P\left\\{{\bm{M}}~{}|~{}{\bm{X}},{\bm{Y}},\bm{x}_{n+1}\right\\}}{P\left\\{\bm{x}_{n+1}|{\bm{M}},{\bm{X}},{\bm{Y}}\right\\}}\,\mathrm{d}{\bm{M}}$
(F.1)
$\displaystyle=\arg\max_{\bm{e}_{j},~{}j\in[k]}\int{P\left\\{\bm{x}_{n+1}~{}|~{}\bm{y}=\bm{e}_{j},{\bm{M}}\right\\}}\,\frac{P\left\\{{\bm{M}}~{}|~{}{\bm{X}},{\bm{Y}}\right\\}}{P\left\\{\bm{x}_{n+1}|{\bm{X}},{\bm{Y}}\right\\}}\,\mathrm{d}{\bm{M}}$
$\displaystyle=\arg\max_{\bm{e}_{j},~{}j\in[k]}\int{P\left\\{\bm{x}_{n+1}~{}|~{}\bm{y}=\bm{e}_{j},{\bm{M}}\right\\}}\,{P\left\\{{\bm{M}}~{}|~{}{\bm{X}},{\bm{Y}}\right\\}}\,\mathrm{d}{\bm{M}}$
(F.2)
$\displaystyle=\arg\max_{\bm{e}_{j},~{}j\in[k]}\int{P\left\\{\bm{x}_{n+1}~{}|~{}\bm{y}=\bm{e}_{j},{\bm{M}}\right\\}}\,{P\left\\{{\bm{X}}~{}|~{}{\bm{M}},{\bm{Y}}\right\\}}\,P\\{{\bm{M}}\\}\,\mathrm{d}{\bm{M}}$
(F.3)

To arrive in (F.1) we used that
$P\left(\bm{y}=\bm{e}_{j}~{}|~{}{\bm{M}},{\bm{X}},{\bm{Y}}\right)=\pi,~{}\forall
j\in[k]$ and
$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(\bm{x}_{n+1}~{}|~{}{\bm{X}},{\bm{Y}})>0$ is independent of the variable of
integration ${\bm{M}}$ and of the optimization variable $j$. For the same
reasons, in (F.3) we have ignored the normalizing term $P({\bm{X}}|{\bm{Y}})$.

Recalling that ${\bm{z}}_{n+1}\sim\mathcal{N}(\bm{0},{\bm{I}}_{d})$, we have
that
$P\left(\bm{x}_{n+1}~{}|~{}\bm{y}=\bm{e}_{j},{\bm{M}}\right)\propto\exp\left(-\left\|\bm{x}_{n+1}-\bm{\mu}_{j}\right\|_{\ell_{2}}^{2}\big{/}2\right)$
where $\propto$ hides constant positive terms. Moreover, the posterior
probability of the mean matrix given the training data is given by

$\displaystyle P\left({\bm{X}}~{}|~{}{\bm{M}},{\bm{Y}}\right)\cdot
P\left({\bm{M}}\right)$
$\displaystyle\propto\exp\left(-\frac{\left\|{\bm{X}}-{\bm{M}}{\bm{Y}}\right\|_{\ell_{2}}^{2}}{2}\right)\cdot\exp\left(-\frac{\left\|{\bm{M}}\right\|_{\ell_{2}}^{2}}{2(\mu^{2}/d)}\right)$
$\displaystyle\propto\prod_{c=1}^{k}\left\\{\exp\left(-\frac{\left\|\bm{\mu}_{c}\right\|_{\ell_{2}}^{2}}{2(\mu^{2}/d)}\right)\cdot\prod_{i\in\mathcal{C}_{c}}\exp\left(-\frac{\left\|\bm{x}_{i}-\bm{\mu}_{c}\right\|_{\ell_{2}}^{2}}{2}\right)\right\\},$
(F.4)

where we denote by $\mathcal{C}_{c}$ the collection of training samples that
belong to class $c\in[k]$, i.e.
$\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:

$\displaystyle\hat{\bm{y}}_{n+1}=\arg\max_{j\in[k]}~{}{\mathcal{I}}(j,\mathcal{C}_{j}\cap\\{n+1\\})\cdot\prod_{\begin{subarray}{c}c=1\\\
c\neq j\end{subarray}}^{k}{\mathcal{I}}(c,\mathcal{C}_{c})\,,$ (F.5)

where for $\ell\in[k]$ and a subset $\mathcal{A}\subset[n+1]$ we denote

${\mathcal{I}}(\ell,\mathcal{A}):=\int\mathrm{d}{\bm{\mu}_{\ell}}\,\exp\left(-\frac{\left\|\bm{\mu}_{\ell}\right\|_{\ell_{2}}^{2}}{2(\mu^{2}/d)}\right)\cdot\exp\left(-\sum_{i\in\mathcal{A}}\frac{\left\|\bm{x}_{i}-\bm{\mu}_{\ell}\right\|_{\ell_{2}}^{2}}{2}\right)\,.$

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

$\displaystyle{\mathcal{I}}(\ell,\mathcal{A})$
$\displaystyle:=\sqrt{\frac{(d/\mu^{2}+|\mathcal{A}|)}{(2{\varpi})^{d}}}\exp\left(-\frac{\left(1-\frac{1}{d/\mu^{2}+|\mathcal{A}|}\right)}{2}\sum_{i\in\mathcal{A}}\left\|\bm{x}_{i}\right\|_{\ell_{2}}^{2}+\frac{1}{2\left(d/\mu^{2}+|\mathcal{A}|\right)}\sum_{i\in\mathcal{A}}\langle\bm{x}_{i},\sum_{\begin{subarray}{c}j\in\mathcal{A}\\\
j\neq i\end{subarray}}\bm{x}_{j}\rangle\right)$
$\displaystyle:=\sqrt{\frac{(d/\mu^{2}+|\mathcal{A}|)}{(2{\varpi})^{d}}}\exp\left(-\frac{1}{2\left(\frac{d}{\mu^{2}}+{|\mathcal{A}|}\right)}\left(\left(\frac{d}{\mu^{2}}+|\mathcal{A}|-1\right)\sum_{i\in\mathcal{A}}\left\|\bm{x}_{i}\right\|_{\ell_{2}}^{2}-\sum_{i\in\mathcal{A}}\langle\bm{x}_{i},\sum_{\begin{subarray}{c}j\in\mathcal{A}\\\
j\neq i\end{subarray}}\bm{x}_{j}\rangle\right)\right)$
$\displaystyle:=\sqrt{\frac{(d/\mu^{2}+|\mathcal{A}|)}{(2{\varpi})^{d}}}\exp\left(-\frac{1}{2\left(\frac{d/n}{\mu^{2}}+{|\mathcal{A}|/n}\right)}\left(\left(\frac{d/n}{\mu^{2}}+\frac{|\mathcal{A}|}{n}-\frac{1}{n}\right)\sum_{i\in\mathcal{A}}\left\|\bm{x}_{i}\right\|_{\ell_{2}}^{2}-\frac{1}{n}\sum_{i\in\mathcal{A}}\langle\bm{x}_{i},\sum_{\begin{subarray}{c}j\in\mathcal{A}\\\
j\neq i\end{subarray}}\bm{x}_{j}\rangle\right)\right).$

Using this in (F.5) we have that

$\displaystyle\hat{\bm{y}}_{n+1}=\arg\max_{j\in[k]}~{}{\mathcal{I}}(j)\cdot\exp\left(-\frac{1}{2\left(\frac{d/n}{\mu^{2}}+\frac{n_{j}+1}{n}\right)}\left(\left(\frac{d/n}{\mu^{2}}+\frac{n_{j}}{n}\right)\left\|\bm{x}_{n+1}\right\|_{\ell_{2}}^{2}-\frac{2}{n}\langle\bm{x}_{n+1},\sum_{\begin{subarray}{c}\ell\in\mathcal{C}_{j}\end{subarray}}\bm{x}_{\ell}\rangle\right)\right),$
(F.6)

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

$\displaystyle{\mathcal{I}}(j):=\left\\{\prod_{\begin{subarray}{c}c=1\\\ c\neq
j\end{subarray}}^{k}e^{-\frac{1}{2\xi(n_{c})}\left(\left(\xi(n_{c})-\frac{1}{n}\right)\sum_{i\in\mathcal{C}_{c}}\left\|\bm{x}_{i}\right\|_{\ell_{2}}^{2}-\frac{1}{n}\sum_{i\in\mathcal{C}_{c}}\langle\bm{x}_{i},\sum_{\begin{subarray}{c}\ell\in\mathcal{C}_{c}\\\
\ell\neq i\end{subarray}}\bm{x}_{\ell}\rangle\right)}\right\\}\cdot
e^{-\frac{1}{2\left(\xi(n_{j})+\frac{1}{n}\right)}\left(\xi(n_{j})\sum_{i\in\mathcal{C}_{j}}\left\|\bm{x}_{i}\right\|_{\ell_{2}}^{2}-\frac{1}{n}\sum_{i\in\mathcal{C}_{j}}\langle\bm{x}_{i},\sum_{\begin{subarray}{c}\ell\in\mathcal{C}_{j}\\\