# Classifying Overlapping Gaussian Mixtures in High Dimensions: From Optimal
Classifiers to Neural Nets

Khen Cohen∗, Yaron Oz
Raymond and Beverly Sackler School of Physics and Astronomy
Tel-Aviv University
Tel-Aviv 69978, Israel
<EMAIL_ADDRESS>
&Noam Levi
École Polytechnique Fédérale de Lausanne (EPFL)
Switzerland
<EMAIL_ADDRESS>
Equal Contribution

###### Abstract

We derive closed-form expressions for the Bayes optimal decision boundaries in
binary classification of high dimensional overlapping Gaussian mixture model
(GMM) data, and show how they depend on the eigenstructure of the class
covariances, for particularly interesting structured data. We empirically
demonstrate, through experiments on synthetic GMMs inspired by real-world
data, that deep neural networks trained for classification, learn predictors
which approximate the derived optimal classifiers. We further extend our study
to networks trained on authentic data, observing that decision thresholds
correlate with the covariance eigenvectors rather than the eigenvalues,
mirroring our GMM analysis. This provides theoretical insights regarding
neural networks’ ability to perform probabilistic inference and distill
statistical patterns from intricate distributions.

##  1 Introduction

There is a well-accepted understanding in machine learning that inherent
correlations or statistical structure within the data play a crucial role in
enabling effective modeling. The presence of structure in the data provides
important context that machine learning algorithms can leverage for accurate
and meaningful knowledge extraction and generalization. Harnessing any
inherent structure is widely seen as an important factor for achieving
successful learning outcomes.

Determining the characteristics of natural datasets that contribute to
effective training and good generalization has tremendous theoretical and
practical significance. Regrettably, real-world datasets are typically viewed
as samples drawn from some unknown, intricate, high-dimensional underlying
population, rendering analyses that depend on accurately modeling the true
distribution enormously challenging. Connecting properties of the data to
learning outcomes is difficult because real data does not readily yield its
generative process or population distribution.

While fully characterizing natural data distributions remains difficult,
progress can still be made by considering simplified models that capture key
aspects of the true distribution. One approach is to model the data via the
moment expansion of the underlying population. For a random vector
$X\in\mathbb{R}^{d}$ representing a single sample, its distribution can be
described by the moment generating function
$M_{X}(t)=\sum_{n=0}^{\infty}t^{n}\mathbb{E}(X^{n})/n!$. Approximating this by
retaining a finite number of moments provides a tractable model of the data
distribution.

A simple starting point is the Gaussian model, which posits that a dataset’s
statistical properties can be fully captured by its first two moments - the
mean vector and sample covariance matrix. This assumes higher-order
dependencies are negligible such that the distribution is characterized solely
by its first and second central statistical moments.

In the last decade, analyzing the behavior of machine learning algorithms on
idealized i.i.d. Gaussian datasets has emerged as an important area of
research in high-dimensional statistics (Donoho and Tanner, 2009; Korada and
Montanari, 2011; Monajemi et al., 2013; Candès et al., 2020; Bartlett et al.,
2020). Studying these simplified synthetic distributions has provided valuable
insights, such as neural network scaling laws (Maloney et al., 2022; Kaplan et
al., 2020), universal convergence properties (Seddik et al., 2020) and an
improved understanding of the "Double Descent" phenomenon (Montanari and
Saeed, 2022). Examining learning on Gaussian approximations of real data
distributions has helped establish foundational understandings of algorithmic
behavior in high-dimensional settings.

A particular case of interest is that of data which is divided into a set
number of classes. In certain instances, the data can be described by a
mixture model, where each sample is generated separately for each class. The
simplest example of such distributions is that of a Gaussian Mixture Model
(GMM), on which we focus in this work.

Gaussian mixtures are a popular model in high-dimensional statistics since,
besides being an universal approximator, they often lead to mathematically
tractable problems. Indeed, a recent line of work has analyzed the asymptotic
performance of a large class of machine learning problems in the proportional
high-dimensional limit under the Gaussian mixture data assumption, see e.g.
Mai and Liao (2019); Mignacco et al. (2020a); Taheri et al. (2020); Kini and
Thrampoulidis (2021); Wang and Thrampoulidis (2021); Refinetti et al. (2021);
Loureiro et al. (2021a).

In light of the insights gained by studying the unsupervised, perspective on
GMM classification, we focus here on a complementary direction. Namely, we
consider the supervised setting, where a single neural network is trained on a
dataset generated from a GMM, to perform binary classification where the true
labels are given. To connect the GMM with real-world data, we restrict
ourselves to the case of strictly non-linearly separable GMMs, such that the
class means difference is negligible compared to the class covariance
difference. We refer to this setup as overlapping GMM classification.

Our main contributions are as follows:

  * •

In Section 3, we derive the Bayes Optimal Classifiers (BOCs) and decision
boundaries for overlapping GMM classification problems, both in the population
and the empirical limits.

  * •

In Section 4, under the assumptions of correlated features and Haar
distributed eigenvectors, we are able to provide approximate closed form
equations for the decision boundaries and discriminator distributions, as a
function of the eigenvalues and eigenvectors of the different class
covariances.

  * •

In Section 5, we present empirical evidence, demonstrating that some networks
trained to perform binary classification on a GMM, approximate the BOC. We
further relate these observations to existing results, namely convergence of
homogeneous networks to a KKT point.

  * •

In Section 6, we demonstrate empirically that for high dimensional data, the
covariance eigenvectors and not the eigenvalues determine the classification
threshold for deep neural networks trained on real-world datasets, and provide
an explanation inspired by our results on GMMs. In Section 7 we summarize our
conclusions and outlook.

##  2 Background and Related Work

There has been significant work aimed at understanding the classification
capabilities of Gaussian mixture models (GMMs), that is, recovering the
cluster label for each data point rather than using labels provided by a
teacher. For binary classification problems, examples include methods proposed
by Mai et al. (2019), Mignacco et al. (2020b), and Deng et al. (2022), with
the latter also demonstrating an equivalence between classification and
single-index models. In multi-class settings, (Thrampoulidis et al., 2020)
analyzed the performance of ridge regression classifiers. The most general
results in this area are from Gaussian Mix Group (Loureiro et al., 2021b),
which considers the use of GMMs with any convex loss function for
classification. Additional results regarding kernels and GMMs include
(Couillet et al., 2018; Liao and Couillet, 2019; Kammoun and Couillet, 2023;
Refinetti et al., 2021), precise asymptotics (Loureiro et al., 2021c) and the
connection between GMMs and real-world vision data (Ingrosso and Goldt, 2022).
Additionally, work has been done to understand Quadratic Discriminant Analysis
(QDA) in high dimensions (Elkhalil et al., 2017; Ghojogh and Crowley, 2019;
Das and Geisler, 2021), which we rely heavily upon in this work.

##  3 Overlapping Gaussian Mixtures in High Dimensions

Consider the task of performing binary classification on samples drawn from a
two class GMM, where the underlying Gaussian distributions have similar means
but different covariance matrices. We generate two datasets
$\mathcal{D}_{A},\mathcal{D}_{B}$ of equal size $N$, and the same data
dimensions $x\in\mathbb{R}^{d}$. Here, each sample vector is drawn from a
jointly normal distribution, either $x\sim\mathcal{N}(\mu_{A},\Sigma_{A})$ or
$x\sim\mathcal{N}({\mu}_{B},\Sigma_{B})$, where the two means and covariance
matrices distinguish between the two classes. In this setup, we assign a label
value $y=1$ to samples drawn from class $A$ and $y=-1$ for the $B$ class
samples. Throughout this work, we consider the large number of features and
large number of samples limit $d,N\to\infty$, while the ratio
$d/N=\gamma\in\mathbb{R^{+}}$ is constant. We denote $\gamma\to 0$ as the
population limit, while $\gamma\to 1$ as the empirical limit.

In the following sections, we discuss the optimal classifiers obtained first
from the population and then the empirical data distributions. We then apply
our results to a simplified model dataset which captures some of the
properties of real-world datasets, and distinguish the roles of the covariance
eigenvectors and eigenvalues for this case.

Figure 1: Probability density function of $\beta(x)$, evaluated on GMMs with
different class covariance matrices. Blue and red bins indicate samples drawn
from the classes $A$ and $B$, respectively. Solid curves represent the
numerical evaluation of Equation 3 as a generalized $\chi^{2}$ distribution.
Dots indicate the values given by Equation 5, as well as by Equation 9. Left:
$\beta$ distribution for covariances with the same basis but different
spectra. Center: $\beta$ distribution for covariances with the same spectrum
but different random bases. Right: $\beta$ distribution comparison for
covariances with both different spectra and different bases, shown in gray.
Here, we take $\alpha_{A}=0.5,\Delta\alpha=-0.3,d=100$.

###  3.1 Optimal Classification on Population Data

The Bayes-optimal classifier (BOC) assigns each sample to the class that
maximizes the expected value gain or its logarithm defined as

$\displaystyle
G(X=x)=\arg\max_{Y=y}p(X=x|Y=y)p(Y=y)=\arg\max_{Y=y}\log(p_{X|Y})=\arg\max_{Y=y}\beta_{Y},$
(1)

where $p(Y=y)=1/C$ is the class density with $C$ the number of classes in the
case of a balanced dataset. Here, $p(X=x|Y=y)=\exp{(\beta_{Y})}$ is the
conditional in-class density. The decision boundary between any two classes
$y,y^{\prime}\in C$ is given when
${\beta_{Y=y}}(X=x)={\beta_{Y=y^{\prime}}}(X=x)$, where the probabilities of a
sample being classified as $y$ or $y^{\prime}$ are equal. For a two-class GMM,
this is obtained by simply equating the class conditional probability
densities, i.e., $p_{A}=p_{B}$, leading to the quadratic decision rule,
derived in Elkhalil et al. (2017) as well as Das and Geisler (2021), to pick
class $A$ if

$\displaystyle\beta(x)=\frac{1}{2}(x^{T}Qx-2q^{T}x+c)>0,\quad
Q=(\Sigma_{B}^{-1}-\Sigma_{A}^{-1}),\quad
q=(\Sigma_{B}^{-1}\mu_{B}-\Sigma_{A}^{-1}\mu_{A}),$ (2) $\displaystyle
c=\mu_{B}^{T}\Sigma_{B}^{-1}\mu_{B}-\mu_{A}^{T}\Sigma_{A}^{-1}\mu_{A}-\log\left(\frac{|{\Sigma_{A}}|}{|\Sigma_{B}|}\right)\
,$

where $|O|$ is the determinant of the matrix $O$. This quadratic $\beta(x)$ is
the Bayes classifier, or the Bayes decision variable that, when compared to
zero, maximizes expected gain.

Note that the $q^{T}x$ term in Equation 2 determines the linear decision
boundary, while $x^{T}Qx$ determines the nonlinear part. For the rest of this
work, we consider the nonlinear, overlapping mixture configuration, where
$\mu_{A}=\mu_{B}=\mu$, while $\Sigma_{A}\neq\Sigma_{B}$, as a proxy for the
behavior of many real world datasets (Schilling et al., 2021). Without loss of
generality, we may set $\mu=0$, since we can always shift away the mean of the
distributions, resulting in the simpler form of $\beta(x)$

$\displaystyle\beta(x)=\frac{1}{2}(x^{T}Qx+c)=\frac{1}{2}(\textrm{Tr}(Qxx^{T})+c)=\frac{1}{2}\sum_{i}q_{i}\tilde{\chi}^{2}_{1}-\frac{a}{2},$
(3)

where the last transition in Equation 3 is valid since $\beta$ is given by a
sum of $\chi^{2}$ distributed random variables, known as the generalized chi-
squared distribution (Das and Geisler, 2021), where
$a\sim\mathcal{N}(\log(|\Sigma_{B}^{-1}\Sigma_{A}|),0)$, $q_{i}$ are the
eigenvalues of $\Sigma_{Y}^{1/2}Q\Sigma_{Y}^{1/2}$, and $Y=A,B$. The
probability density function (PDF) of $\beta$ does not have a closed analytic
form, but can be obtained numerically by various methods (Ferrari, 2019; Das
and Geisler, 2021). In Figure 1, we show the $\beta$-distribution for high-
dimensional Gaussian data with different covariances, obtained both
empirically from Equation 3 and by drawing from a sum of $\chi^{2}$ variables
matching Equation 3.

While the full PDF of $\beta$ cannot be written in a simple form, the
empirical expectation value on a given set of measurements can be derived.
Namely, given a set of $N$ samples $x_{i},i=1,\ldots,N$, taken from a normal
distribution $\mathcal{N}(0,M),M\in\mathbb{R}^{d\times d}$, the expectation
value of the Bayes classifier $\beta$ on these samples is simply

$\displaystyle\langle\beta_{M}\rangle=\frac{1}{2}\left[\textrm{Tr}(QM)-\log\left({|\Sigma_{B}^{-1}|{|{\Sigma_{A}}|}}\right)\right],$
(4)

where $\langle\cdot\rangle$ indicates averaging over the data distribution. It
can be used to define what it means for a sample to belong to class $A$ or
$B$, by computing the average distance between the expectation value of
$\beta$ on any dataset and the class values, given by

$\displaystyle\langle\beta_{A}\rangle=-d+\textrm{Tr}(\Sigma_{B}^{-1}\Sigma_{A})+c>0,\quad\langle\beta_{B}\rangle=d-\textrm{Tr}(\Sigma_{A}^{-1}\Sigma_{B})+c<0,$
(5)

where we averaged over $N\to\infty$ samples from the $A,B$ distributions, and
$c=-\log(|\Sigma_{B}^{-1}\Sigma_{A}|)$. We show these values as black dots in
Figure 1.

Figure 2: Comparison between BOC and quadratic network trained on data with
different covariance spectrum. Left: accuracy, measured as the sigmoid
function acting on the network/BOC output. Center: network/BOC output
distribution. Blue indicates network results while red is the BOC prediction.
Here, the network is trained with $d_{h}=100$ and $N=200k$ to $91.27\%$
training accuracy and reaches $91.1\%$ with $\alpha_{A}=0.2$ and
$\alpha_{B}=0.3$. Right: Results for KKT convergence of a quadratic NN trained
on an $\alpha=0.2$ correlated GMM, where the class covariances share the same
spectrum but are given a different random basis. Red indicates the BOC
prediction, blue is the quadratic network output before the softmax function,
when trained with $d_{h}=100$, and green are the KKT predictions with
$\lambda_{a}\propto 1/N$ and $d_{h}=100$. The data dimensions are similar to
the main text, i.e., $d=100$, using $N=100k$ samples, the network reached
$100\%$ training and $99.9\%$ test accuracy.

###  3.2 Empirical Optimal Classification

So far, we have assumed that we have access to the population covariances
$\Sigma_{A},\Sigma_{B}$. Here, we comment on the classification problem
changes when relaxing this assumption and consider the empirical covariances
which are obtained at finite dimension to sample number ratio $\gamma=const.$

In particular, we consider the case where we have access only to finitely
sampled versions of the population moments, but have infinitely many samples
of this noisy distribution. Here, the distributions in question for classes
$A,B$ are still gaussian, but their respective moments are not the population
ones, implying that the BOC is deformed. The covariance matrices we have
access to are given by measurements, that define the empirical covariances for
each class
$\Sigma_{N,Y}=\frac{1}{N}X_{Y}X_{Y}^{T}\simeq\Sigma_{Y}+\sqrt{\gamma}R_{Y}$,
where $X_{Y}\in\mathbb{R}^{d\times N}$ is the design matrix, and $R_{Y}$ is a
random matrix with $\mathcal{O}(1)$ eigenvalues (Biroli and Mézard, 2023)111
This result holds in the limit of $d^{2}\ll N$. In the case of $d\ll N$,
$Q_{N}$ is a Wishart matrix, which requires more careful treatment such as
given in Liao and Couillet (2019); Tiomoko et al. (2019); Zhang et al. (2018),
which we leave for future work.. We can therefore rewrite the empirical BOC as

$\displaystyle\\!\\!\\!\beta_{N}(x)$
$\displaystyle\simeq\frac{x^{T}\left(Q_{N}\right)x}{2}+\frac{c_{N}}{2},~{}Q_{N}\simeq
Q+\sqrt{\gamma}({R_{B}}-{R_{A}}),~{}c_{N}\\!=c-\log\left(\frac{|\mathbb{I}+\sqrt{\gamma}{\Sigma_{B}R_{B}}|}{|\mathbb{I}+\sqrt{\gamma}{\Sigma_{A}R_{A}}|}\right),$
(6)

where we expand to leading order in $\gamma=\frac{d}{N}$. This implies that
the empirical $\beta_{N}$ distribution remains a generalized $\chi^{2}$ with
the substitution of coefficients to be the eigenvalues of
$\Sigma_{Y}^{1/2}Q_{N}\Sigma_{Y}^{1/2}$ and the constant $c$ is shifted by the
$\log$ term. Therefore, the empirical deviation from the BOC
$\Delta\beta_{N}=\|\beta_{N}-\beta\|$ is estimated by

$\displaystyle\\!\\!\\!\\!\Delta\beta_{N}$
$\displaystyle\simeq\sqrt{\gamma}|\textrm{Tr}({R_{B}}(xx^{T}-\Sigma_{B})-{R_{A}}(xx^{T}-\Sigma_{A}))|,$
(7)

which decreases with $\sqrt{\gamma}$. This setting is exhibited in Figure 1,
which is meant to mimic the real-world setting, where we do not have access to
the population covariance, and we can only estimate how close the empirical
moments are to the population ones.

##  4 Analysis for a Toy Model of Complex Data

In order to connect the BOC with real world classification tasks, we focus on
a specific modeling of natural data. Here, we consider a simple model of a
correlated GMM inspired by the neural scaling law literature (Maloney et al.,
2022; Montanari and Saeed, 2022; Levi and Oz, 2023a) as well as the signal
recovery literature (Loureiro et al., 2021d, 2022). Concretely, we assume a
power law scaling spectrum, and a different basis matrix for each class

$\displaystyle\Lambda_{Y,ij}=i^{-1-\alpha_{Y}}\delta_{ij},\qquad\Sigma_{Y}=O_{Y}^{T}\Lambda
O_{Y},$ (8)

where $\delta_{ij}$ is the $d$ dimensional identity matrix and $O_{Y}$ is a
random orthogonal matrix. In the following sections, we distinguish the role
of eigenvalues and eigenvectors in determining the BOC for overlapping GMMs.

###  4.1 Diagonal Correlated Covariances

As a first example, we consider the case where $\Sigma_{A},\Sigma_{B}$ are
both diagonal matrices with different eigenvalues, such that
$O_{A}=O_{B}=\mathbb{I}$, and $\Sigma_{A}=\Lambda_{A},\Sigma_{B}=\Lambda_{B}$.
Given a model of power law scaling dataset $A$ with $\alpha_{A}$, we can
define the power law exponent of the dataset $B$ to be
$\alpha_{B}=\alpha_{A}+\Delta\alpha$. Here, we can obtain explicit expressions
for $\beta$ as

$\displaystyle\\!\\!\\!\langle\beta_{A}\rangle\\!=\\!\frac{1}{2}\left(H_{d}^{(-\Delta\alpha)}\\!\\!-\\!d-\Delta\alpha\log(\Gamma(d+1))\right),~{}\langle\beta_{B}\rangle\\!=\\!\frac{1}{2}\left(d-H_{d}^{(\Delta\alpha)}\\!\\!-\Delta\alpha\log\left(\Gamma(d+1)\right)\right),$
(9)

where $H_{d}^{i}$ is the Harmonic number, and $\Gamma(x)$ is the Gamma
function. In the limit of $d\to\infty$, one obtains that
$\langle\beta_{A}\rangle\propto\frac{1}{4}\left(2d\left(\frac{d^{\Delta\alpha}}{\Delta\alpha+1}-1\right)-(2d\Delta\alpha+\Delta\alpha)\log(d)\right)$
at leading order. This is a rapidly growing function in $d$, showing that even
a small difference in spectra can be magnified simply by dimensionality,
leading to correct classification.

###  4.2 Rotated Correlated Covariances

Next, we consider the case where $\Sigma_{A},\Sigma_{B}$ share the same
spectrum, but are rotated with respect to one another, such that
$\Sigma_{A}=O_{A}^{T}\Lambda O_{A},\Sigma_{B}=O_{B}^{T}\Lambda O_{B}$. The
sample expectation for the BOC reads

$\displaystyle\langle\beta_{A}\rangle=-\langle\beta_{B}\rangle=\frac{1}{2}\left(\textrm{Tr}(\mathcal{O}^{T}\Lambda^{-1}\mathcal{O}\Lambda)-d\right),$
(10)

where $\mathcal{O}=O_{A}O_{B}^{T}$ is itself a random orthogonal matrix. In
the limit of $d\to\infty$, we may replace the expectation value over the
dataset with its integral with respect to the Haar measure on the rotation
group, and obtain a closed form expression as

$\displaystyle\langle\beta_{A}\rangle_{\mathcal{O}}$
$\displaystyle=\frac{1}{2}\left[\frac{1}{d}\textrm{Tr}(\Lambda^{-1})\textrm{Tr}(\Lambda)-d\right]=\frac{1}{2}\left(\frac{H_{d}^{(-1-\alpha)}H_{d}^{(1+\alpha)}}{d}-d\right),$
(11)

where the final expression is given for a power law scaling dataset. In the
limit of $d\to\infty$, Equation 11 grows as
$\langle\beta_{A}\rangle_{\mathcal{O}}\propto\frac{1}{2}\left(\frac{1}{\alpha(\alpha+2)}+1\right)d-\frac{d^{\alpha+1}\zeta(\alpha+1)}{2(\alpha+2)}$,
which is a faster growing function than Equation 9 for $\Delta\alpha/\alpha<1$
indicating that it requires a large difference in the relative eigenvalue
scaling exponents $\Delta\alpha/\alpha$ to shift the BOC distribution, when
compared to the effect of a basis difference. We show this effect explicitly
in the different panels of Figure 1. This implies that the eigenvectors are
likely to play a more significant role than eigenvalues in datasets which are
well modeled by the above setup.

Figure 3: Two Gaussian classification task, with the same eigenvectors but
different eigenvalues bulk: the first one with $\alpha=0.5$ and the second
with $\alpha+\Delta\alpha$. The model is FC, averaged over 5 training runs -
the solid line represents the mean accuracy and shaded area represents 1-sigma
error bars.

##  5 Neural Networks as Nearly Optimal Classifiers

Having established the BOC behavior on overlapping GMMs, we can now ask how
these conclusions translate to the learning process in some neural network
architectures. Namely, NNs have been shown to achieve optimal classification
in some cases (Radhakrishnan et al., 2023). We demonstrate that NNs quite
generically approximate the quadratic BOC, and provide some intuition as to
how it may occur by appealing to the directional convergence of gradient
descent to a KKT point.

###  5.1 Neural Network Classifier

Let
$\mathcal{D}=\\{(\mathbf{x}_{a},y_{a})\\}_{a=1}^{N}\subseteq{\mathbb{R}}^{d}\times\\{-1,1\\}$
be a binary classification training dataset. Let
$\Phi({\boldsymbol{\theta}};\cdot):{\mathbb{R}}^{d}\to{\mathbb{R}}$ be a
neural network parameterized by ${\boldsymbol{\theta}}\in{\mathbb{R}}^{p}$.
For a loss function $\ell:{\mathbb{R}}\to{\mathbb{R}}$ the empirical loss of
$\Phi({\boldsymbol{\theta}};\cdot)$ on the dataset $\mathcal{D}$ is ${\cal
L}({\boldsymbol{\theta}}):=\sum_{a=1}^{N}\ell(y_{a}\Phi({\boldsymbol{\theta}};\mathbf{x}_{a}))$.
We focus on the _logistic loss_ (a.k.a. _binary cross entropy_), namely,
$\ell(q)=\log(1+e^{-q})$.

The network function $\Phi({\boldsymbol{\theta}};\cdot)$, can therefore be
interpreted as the equivalent of the BOC variable $\beta$ in the previous
sections. With this intuition in mind, we expect that neural networks trained
to perform binary classification should converge to the BOC, i.e.,
$\Phi({\boldsymbol{\theta}};\cdot)\to\beta$ with sufficient with sufficient
network expressivity, number of samples and successful optimization dynamics.

The simplest network architecture that one may employ is a linear classifier,
where $\Phi(z)=z$. Such linear networks of any depth cannot approximate the
BOC given by Equation 3, as it is fully nonlinear.

Next, we consider a non-trivial case, by taking a two layer network, with the
activation function $\sigma(z)=z^{2}$, i.e. the quadratic activation.
Quadratic activations are natural in the context of signal processing, where
often detectors can only measure the amplitude (and not phase) of a signal,
e.g. the phase retrieval problem (Jaganathan et al., 2015; Dong et al., 2023).
It has also gained in popularity recently as a prototypical non-convex
optimization problem with strict saddles, e.g. Candes et al. (2014); Chen et
al. (2019); Arnaboldi et al. (2024); Martin et al. (2023). It has also been
shown such two-layer blocks can be used to simulate higher-order polynomial
neural networks and sigmoidal activated neural networks (Livni et al., 2014;
Soltani and Hegde, 2018). Concretely, a sufficiently expressive network
function to fully realize the BOC is

$\displaystyle\Phi=v^{T}(Wx)^{2}+b=\sum_{i}\sum_{\mu,\nu}v_{i}W_{i\mu}W_{i\nu}x_{\mu}x_{\nu}+b,$
(12)

where $W\in\mathbb{R}^{d_{h}\times d}$, $v^{T}\in\mathbb{R}^{d_{h}}$ are
weight matrices at the input and hidden layers, $d_{h}$ is the number of
hidden units, and $b$ is the bias at the last layer. In this setup, a
predictor function can be easily matched to the BOC by matching

$\displaystyle\sum_{i=1}^{d}v_{i}W_{i\mu}W_{i\nu}=\frac{1}{2}Q_{\mu\nu},\qquad
b=\frac{c}{2},$ (13)

where the number of hidden neurons correspond to the minimal number required
to match the rank of the BOC matrix $Q$, hence $d_{h}\geq d$. Since this
network is expressive enough to fully reproduce the BOC, our analysis in
Section 3.1 holds, and the importance of eigenvalues and eigenvectors is
manifest. We show the result of a quadratic network converging to the BOC in
Figure 2 where $d=d_{h}=100$.

###  5.2 Karush–Kuhn–Tucker (KKT) Convergence

Another source of intuition for the value of $d_{h}$ may be found in the KKT
convergence equation (Ji and Telgarsky, 2020; Lyu and Li, 2020), which states,
in a nutshell, that the weights of a homogeneous neural network
222Homogeneity, i.e. $f(L\theta)=L^{a}f(\theta)$, where $f$ is the output of
the network, $\theta$ are the network parameters, and $L,a\in\mathbb{R}$, a
quadratic network satisfies this requirement as it is composed of monomials.,
trained with the logistic loss for binary classification, using an
infinitesimal learning rate $\eta\to 0$ and infinite iterations $t\to\infty$,
will converge in direction to a KKT fixed point given by

$\displaystyle\tilde{{\boldsymbol{\theta}}}=\sum_{a=1}^{N}\lambda_{a}y_{a}\nabla_{\boldsymbol{\theta}}\Phi(\tilde{{\boldsymbol{\theta}}};\mathbf{x}_{a}),$
(14)

where $\tilde{{\boldsymbol{\theta}}}$ are the network weights at convergence,
$\mathbf{x}_{a}$ are the samples and $y_{a}$ are the sample labels. Here,
$\lambda_{a}\in\mathbb{R}^{+}$, and are nonzero only for samples on the
decision boundary. The full theorem and requirements are provided in Appendix
B. We can apply Equation 14 to the quadratic network as

$\displaystyle\\!\\!\\!\\!\tilde{v}=\sum_{a=1}^{N}\lambda_{a}y_{a}(\tilde{W}x_{a})^{2},~{}\tilde{W}=\sum_{a=1}^{N}\lambda_{a}y_{a}\tilde{v}^{T}x_{a}(Wx_{a}),$
(15)

where we assume no biases (consistent with the case of similar spectrum and
different basis). We posit, that in high dimensions $\lambda_{a}\propto
N^{-1}$, since at large $d$ every sample in the dataset should lie on the
margin and contribute an equal amount, and the decision hyper-surface area
should be comparable to the entire volume (Hsu et al., 2022). Under this
ansatz, the above equations can be solved numerically, and indeed approximate
the BOC as $d_{h}$ increases. We present evidence that the Bayes optimal
decision boundary can be realized by a two-layer neural network with quadratic
activation, and increasing the number of hidden units increases the
probability of converging to this minima under gradient flow in Figure 2.

##  6 Results Extending to Realistic Data and Networks

Figure 4: The classification flipping test on the FMNIST and CIFAR10 datasets
between class 0 and other classes. Top row: training and tests on GMMs
generated from the class covariance matrices. Middle row: training and tests
on real images after whitening, rescaling and coloring. Bottom row:
eigenvector threshold at the flipping point, for FMNIST and CIFAR10, using
GMMs and real images, when either training with a NN or predicting with the
BOC. Blue columns indicate the flipping point on GMM classification, red
indicates the same test on real images, and green shows the predictions of the
BOC given only by information on the covariance matrices, as given by Equation
3. The BOC is only shown for CIFAR10, due to numerical instabilities.

In the previous section, we showed that small differences in eigenvalues and
eigenvectors are magnified by the dimension, making classification easier. We
further demonstrated that differences in the eigenvectors are faster growing
in the dimension of the data, suggesting eigenvectors play a more significant
role in the Bayesian decision boundary. Next, we will provide clear evidence
that these conclusions hold, and extend to real-world datasets, by performing
a set of tests related to the various properties of the covariance matrices.
We train two common network architectures to perform binary classification,
both on GMMs and on real images, exploring their performance as we change the
covariance structure of each class.

Concretely, we consider two popular network architectures: fully connected
(FC) and convolutional neural networks (CNN), tested on the CIFAR10
(Krizhevsky, 2012) and Fashion-MNIST (Xiao et al., 2017) datasets. For each
training procedure, the model is optimized to classify samples from two
classes in the given dataset (class 0 versus classes 1-9 aggregated).

The optimization procedure proceeds as follows: for each class, the samples
are split into training and evaluation subsets. We then compute the covariance
matrix of the training and evaluation subsets separately for the first and
second classes. New synthetic data is generated by sampling from a
multivariate Gaussian distribution with zero mean and the corresponding
covariance matrix. This synthetic data is used to train the model to
distinguish between the two classes (i.e. classify the Gaussians). When
trained on real images, the model is tested on real images and not gaussians.

The optimization objective is the binary cross-entropy loss. The FC
architecture consists of 3 dense layers with 2048 units each utilizing ReLU
activations, followed by a softmax output layer. For CNNs, we employ ResNet-18
(He et al., 2015). All models are optimized with the Adam optimizer, a
learning rate of 0.001, and Cosine Annealing LR scheduling over 50 epochs with
a batch size of 256. The synthetic datasets contain 1,280 images per class,
split 80-20 for training and validation. Our computational resource was a
single NVIDIA GeForce GTX 1650 GPU, with 16GB RAM.

###  6.1 $\Delta\alpha$ Tests on GMMs

We first perform an experiment to quantitatively evaluate a network’s ability
to utilize covariance matrix eigenvalue scaling for classification tasks.
Specifically, two random covariance matrices are constructed using Equation 8,
with one covariance has the spectral scaling of $\alpha$ and the other of
$\alpha+\Delta\alpha$. Critically, both matrices share the same basis of
eigenvectors. Approximately 500 samples are then drawn from each covariance
matrix.

A model is trained to classify between the two classes defined by the distinct
covariance matrices. Figure 3 shows the results for a fully connected network,
as well as for a ResNet18 architecture (He et al., 2015). To assess the
performance variability, the results are averaged over an ensemble of five
independent training runs.

By varying only the scaling exponent between the covariances while keeping the
eigenvectors fixed, this experiment protocol directly tests a network’s
capability to leverage the eigenvalue scaling imparted by the covariance
matrix for discriminative learning. The stability of the classifications
indicate the degree to which these deep models can extract and utilize the
intrinsic geometrical information encoded in the covariance descriptors.

###  6.2 Flip Tests on GMMs Constructed from Real Data

Next, we consider the role of eigenvectors and eigenvalues in classifying
overlapping GMMs, as well as real images. It is known that the eigenvectors
and eigenvalues of the data covariance matrix reveal different aspects of the
data. In order to understand which of the two contains pertinent information
for NN classification performance, we consider the following setup: we study
two datasets, with covariance matrices $C_{1}$ and $C_{2}$, and decompose each
one of them to a set of eigenvectors and eigenvalues:
$\left\\{\lambda^{1}_{i},v^{1}_{i}\right\\}$,
$\left\\{\lambda^{2}_{i},v^{2}_{i}\right\\}$. We then construct a new
covariance matrix by combining some of the eigenvalues and eigenvectors of
$C_{1}$ and $C_{2}$. To quantify the ratio of eigenvectors/eigenvalues coming
from each class, we define two threshold indices $\tau_{\lambda}$ for the
eigenvalues and $\tau_{v}$ for the eigenvectors , where each of them is an
integer between 0 and the dimension of $C_{1},C_{2}$, given by $d$. The new
covariance matrix reads:

$\displaystyle C_{\tau_{\lambda},\tau_{v}}=V^{T}\Lambda V\ ,$ (16)

where $V=V_{\tau_{v}}$ and $\Lambda=\Lambda_{\tau_{\lambda}}$ are the new
basis, and eigenvalues, respectively, which are composed by a mixture between
the two classes:

$\displaystyle v_{i}=\begin{cases}{v^{1}}_{i}&i\leq\tau\\\
{v^{2}}_{i}&\tau<i\\\ \end{cases},\hskip
28.45274pt\lambda_{i}=\begin{cases}{\lambda^{1}}_{i}&i\leq\tau\\\
{\lambda^{2}}_{i}&\tau<i\ .\end{cases}$ (17)

This process is denoted as a flip test, illustrated in Figure 8.

It is worth noting that when composing the rotation matrix $V$, we must ensure
that it is nearly orthogonal to maintain the properties of a true basis, which
is not guaranteed in our scheme. However, we find that this matrix is very
close to orthogonal in the following sense: we define the following error, and
require $e\ll 1$

$\displaystyle
e\equiv\max_{\tau}{\tfrac{1}{d}\left||V^{\dagger}V-I|\right|_{F}},$ (18)

where $\|\cdot\|_{F}$ is the Frobenius norm. Equation 18 sets an upper bound
for the error over all of the thresholds $\tau\in[0,d]$. We practically found
that $e\sim 0.018$, while computing it for a randomly generated matrix
composed of unit norm vectors gives $e\sim 1$.

In Figure 4 (top row), we show the results of training a FC and ResNet18
networks to classify GMMs constructed from different $C_{1},C_{2}$ and tested
on samples drawn from $C_{\tau_{\lambda},\tau_{v}}$. The baseline results are
that the network classifies samples as coming from $C_{1}$ if
$C_{\tau_{\lambda},\tau_{v}}$ has both the eigenvalues and eigenvectors of
$C_{1}$. When constructing $C_{\tau_{\lambda},\tau_{v}}$ with the eigenvectors
of $C_{1}$ and the eigenvalues of $C_{2}$, the network still classifies
samples as coming from $C_{1}$. This is not the case when taking the
eigenvalues from $C_{1}$ but some of the eigenvectors from $C_{2}$. Evidently,
it requires a $\mathcal{O}(100)$ eigenvectors from $C_{2}$ to convince the
network that the samples came from $C_{2}$.

In Figure 4 (bottom row), we show that this result is architecture dependent,
but the existence of an eigenvector flipping point persists throughout our
experiments. We further compare this result with the BOC prediction, showing
that indeed the BOC requires more eigenvectors to be flipped before
determining that samples belong in $C_{2}$ rather than $C_{1}$. This is shown
in detail in Figure 11.

We further analyze the above flip tests for different data sizes used to
construct $C_{1}$ and $C_{2}$. We find that the results do not change above a
certain number of samples, as the covariance matrices are well approximated
beyond that point. Figure 12 shows the results.

###  6.3 Flip Tests on Real Data

To ensure that the overlapping mixture model holds for real-data, we examined
the effect of adding the vector class means computed from the CIFAR10 and
FMNIST datasets. We observed no qualitative change in the results, indicating
that the assumption of zero mean is justified in these cases.

In order to perform the same flipping tests on real images, we simply perform
a whitening transformation, followed by a rescaling, and then rotating to the
basis of choice (Belhumeur et al., 1997). For an illustration of the resulting
images see Figure 9 and Figure 10. This process ensures that the class
covariance matrices of the data have been matched to our design, but
inherently affects all the higher moments of the underlying distribution in an
unpredictable way. Interestingly, as can be seen in Figure 4 (top row), the
occurrence of a clear flipping point which depends on eigenvectors and not
eigenvalues is manifest in real images. In Figure 4 (bottom row) we show a
comparison between the thresholds for GMMs and real images, where the results
indicate a close similarity, particularly for ResNet18 trained on CIFAR10,
which can be attributed to its high expressivity.

##  7 Conclusions and Limitations

We studied binary classification of high dimensional overlapping GMM data as a
framework to model real-world image datasets and to quantify the importance
for the classification tasks of the eigenvalues and the eigenvectors of the
data covariance matrix. We showed that deep neural networks trained for
classification, learn predictors that approximate the Bayes optimal
classifiers, and demonstrated that the decision thresholds for networks
trained on authentic data correlate with the covariance eigenvectors rather
than the eigenvalues, compatible with the GMM analysis. Our results reveal new
theoretical insights to the neural networks’ ability to perform probabilistic
inference and distill statistical patterns from complex distributions.

Limitations: Firstly, we focused on the $\gamma\ll 1$ regime of the estimated
class covariances, while in many real-world cases, the empirical limit may be
more applicable Levi and Oz (2023b). In Section 5.2, our ansatz neglected the
dependence of $\lambda$ on $d/N$ and the spectral density, which nonlinearly
determines the number of points which lie on the decision surface. It might be
possible to derive a similar result to Hsu et al. (2022) by taking a feature
map approach to the quadratic net, but we postpone this to future work.
Finally, while we quantified the significance of the detailed structure of the
data covariance matrix, i.e. its eigenvalues and eigenvectors and their
relative importance for classification, we are still left with the open
question regarding the relative importance of the higher moments of the
distribution, which we leave for a future study.

##  8 Acknowledgements

We would like to thank Amir Globerson, Yohai Bar-Sinai, Yue Lu and Bruno
Loureiro for useful discussions. NL would like to thank G-Research for the
award of a research grant, as well as the CERN-TH department for their
hospitality during various stages of this work. The work of Y.O. is supported
in part by Israel Science Foundation Center of Excellence.

## References

  * Donoho and Tanner (2009) David Donoho and Jared Tanner.  Observed universality of phase transitions in high-dimensional geometry, with implications for modern data analysis and signal processing.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 367(1906):4273–4293, 2009.
  * Korada and Montanari (2011) Satish Babu Korada and Andrea Montanari.  Applications of the lindeberg principle in communications and statistical learning.  _IEEE Transactions on Information Theory_ , 57(4):2440–2450, 2011.  doi: 10.1109/TIT.2011.2112231.
  * Monajemi et al. (2013) Hatef Monajemi, Sina Jafarpour, Matan Gavish, David L. Donoho, Sivaram Ambikasaran, Sergio Bacallado, Dinesh Bharadia, Yuxin Chen, Young Choi, Mainak Chowdhury, Soham Chowdhury, Anil Damle, Will Fithian, Georges Goetz, Logan Grosenick, Sam Gross, Gage Hills, Michael Hornstein, Milinda Lakkam, Jason Lee, Jian Li, Linxi Liu, Carlos Sing-Long, Mike Marx, Akshay Mittal, Albert No, Reza Omrani, Leonid Pekelis, Junjie Qin, Kevin Raines, Ernest Ryu, Andrew Saxe, Dai Shi, Keith Siilats, David Strauss, Gary Tang, Chaojun Wang, Zoey Zhou, and Zhen Zhu.  Deterministic matrices matching the compressed sensing phase transitions of gaussian random matrices.  _Proceedings of the National Academy of Sciences_ , 110(4):1181–1186, 2013.  doi: 10.1073/pnas.1219540110.
  * Candès et al. (2020) Emmanuel J Candès, Pragya Sur, et al.  The phase transition for the existence of the maximum likelihood estimate in high-dimensional logistic regression.  _The Annals of Statistics_ , 48(1):27–42, 2020\.
  * Bartlett et al. (2020) Peter L. Bartlett, Philip M. Long, Gábor Lugosi, and Alexander Tsigler.  Benign overfitting in linear regression.  _Proceedings of the National Academy of Sciences_ , 117(48):30063–30070, 2020.  ISSN 0027-8424.  doi: 10.1073/pnas.1907378117.
  * Maloney et al. (2022) Alexander Maloney, Daniel A. Roberts, and James Sully.  A solvable model of neural scaling laws, 2022.
  * Kaplan et al. (2020) Jared Kaplan, Sam McCandlish, Tom Henighan, Tom B. Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei.  Scaling laws for neural language models, 2020.
  * Seddik et al. (2020) Mohamed El Amine Seddik, Cosme Louart, Mohamed Tamaazousti, and Romain Couillet.  Random matrix theory proves that deep learning representations of gan-data behave as gaussian mixtures, 2020.
  * Montanari and Saeed (2022) Andrea Montanari and Basil N. Saeed.  Universality of empirical risk minimization.  In _Proceedings of Thirty Fifth Conference on Learning Theory_ , volume 178 of _Proceedings of Machine Learning Research_ , pages 4310–4312. PMLR, 02–05 Jul 2022.
  * Mai and Liao (2019) Xiaoyi Mai and Zhenyu Liao.  High-dimensional classification via empirical risk minimization: Improvements and optimality.  _arXiv: 1905.13742_ , 2019.
  * Mignacco et al. (2020a) Francesca Mignacco, Florent Krzakala, Yue Lu, Pierfrancesco Urbani, and Lenka Zdeborova.  The role of regularization in classification of high-dimensional noisy gaussian mixture.  In _International Conference on Machine Learning_ , pages 6874–6883. PMLR, 2020a.
  * Taheri et al. (2020) Hossein Taheri, Ramtin Pedarsani, and Christos Thrampoulidis.  Optimality of least-squares for classification in gaussian-mixture models.  In _2020 IEEE International Symposium on Information Theory (ISIT)_ , pages 2515–2520. IEEE, 2020.
  * Kini and Thrampoulidis (2021) Ganesh Ramachandra Kini and Christos Thrampoulidis.  Phase transitions for one-vs-one and one-vs-all linear separability in multiclass gaussian mixtures.  In _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 4020–4024. IEEE, 2021.
  * Wang and Thrampoulidis (2021) Ke Wang and Christos Thrampoulidis.  Benign overfitting in binary classification of gaussian mixtures.  In _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 4030–4034. IEEE, 2021.
  * Refinetti et al. (2021) Maria Refinetti, Sebastian Goldt, Florent Krzakala, and Lenka Zdeborová.  Classifying high-dimensional gaussian mixtures: Where kernel methods fail and neural networks succeed, 2021.
  * Loureiro et al. (2021a) Bruno Loureiro, Gabriele Sicuro, Cedric Gerbelot, Alessandro Pacco, Florent Krzakala, and Lenka Zdeborova.  Learning gaussian mixtures with generalized linear models: Precise asymptotics in high-dimensions.  In _Advances in Neural Information Processing Systems_ , volume 34, pages 10144–10157, 2021a.
  * Mai et al. (2019) 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, May 2019.  doi: 10.1109/ICASSP.2019.8683376.  ISSN: 2379-190X.
  * Mignacco et al. (2020b) Francesca Mignacco, Florent Krzakala, Yue M. Lu, and Lenka Zdeborová.  The role of regularization in classification of high-dimensional noisy gaussian mixture.  2020b.  doi: 10.48550/ARXIV.2002.11544.  URL https://arxiv.org/abs/2002.11544.
  * Deng et al. (2022) Zeyu Deng, Abla Kammoun, and Christos Thrampoulidis.  A model of double descent for high-dimensional binary linear classification.  _Information and Inference: A Journal of the IMA_ , 11(2):435–495, June 2022.  ISSN 2049-8772.  doi: 10.1093/imaiai/iaab002.  URL https://doi.org/10.1093/imaiai/iaab002.
  * Thrampoulidis et al. (2020) Christos Thrampoulidis, Samet Oymak, and Mahdi Soltanolkotabi.  Theoretical insights into multiclass classification: a high-dimensional asymptotic view.  In _Proceedings of the 34th International Conference on Neural Information Processing Systems_ , NIPS’20, pages 8907–8920, Red Hook, NY, USA, December 2020. Curran Associates Inc.  ISBN 9781713829546.
  * Loureiro et al. (2021b) Bruno Loureiro, Gabriele Sicuro, Cédric Gerbelot, Alessandro Pacco, Florent Krzakala, and Lenka Zdeborová.  Learning gaussian mixtures with generalised linear models: Precise asymptotics in high-dimensions.  2021b.  doi: 10.48550/ARXIV.2106.03791.  URL https://arxiv.org/abs/2106.03791.
  * Couillet et al. (2018) Romain Couillet, Zhenyu Liao, and Xiaoyi Mai.  Classification asymptotics in the random matrix regime.  In _2018 26th European Signal Processing Conference (EUSIPCO)_ , pages 1875–1879, 2018.  doi: 10.23919/EUSIPCO.2018.8553034.
  * Liao and Couillet (2019) Zhenyu Liao and Romain Couillet.  A large dimensional analysis of least squares support vector machines.  _IEEE Transactions on Signal Processing_ , 67(4):1065–1074, February 2019.  ISSN 1941-0476.  doi: 10.1109/tsp.2018.2889954.  URL http://dx.doi.org/10.1109/TSP.2018.2889954.
  * Kammoun and Couillet (2023) Abla Kammoun and Romain Couillet.  Covariance discriminative power of kernel clustering methods.  _Electronic Journal of Statistics_ , 17(1):291 – 390, 2023.  doi: 10.1214/23-EJS2107.  URL https://doi.org/10.1214/23-EJS2107.
  * Loureiro et al. (2021c) Bruno Loureiro, Gabriele Sicuro, Cedric Gerbelot, Alessandro Pacco, Florent Krzakala, and Lenka Zdeborová.  Learning gaussian mixtures with generalized linear models: Precise asymptotics in high-dimensions.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, _Advances in Neural Information Processing Systems_ , volume 34, pages 10144–10157. Curran Associates, Inc., 2021c.  URL https://proceedings.neurips.cc/paper_files/paper/2021/file/543e83748234f7cbab21aa0ade66565f-Paper.pdf.
  * Ingrosso and Goldt (2022) Alessandro Ingrosso and Sebastian Goldt.  Data-driven emergence of convolutional structure in neural networks.  _Proceedings of the National Academy of Sciences_ , 119(40), September 2022.  ISSN 1091-6490.  doi: 10.1073/pnas.2201854119.  URL http://dx.doi.org/10.1073/pnas.2201854119.
  * Elkhalil et al. (2017) Khalil Elkhalil, Abla Kammoun, Romain Couillet, Tareq Y Al-Naffouri, and Mohamed-Slim Alouini.  Asymptotic performance of regularized quadratic discriminant analysis based classifiers.  In _2017 IEEE 27th International Workshop on Machine Learning for Signal Processing (MLSP)_ , pages 1–6. IEEE, 2017.
  * Ghojogh and Crowley (2019) Benyamin Ghojogh and Mark Crowley.  Linear and quadratic discriminant analysis: Tutorial, 2019.
  * Das and Geisler (2021) Abhranil Das and Wilson S. Geisler.  A method to integrate and classify normal distributions.  _Journal of Vision_ , 21(10):1, September 2021\.  ISSN 1534-7362.  doi: 10.1167/jov.21.10.1.  URL http://dx.doi.org/10.1167/jov.21.10.1.
  * Schilling et al. (2021) Achim Schilling, Andreas Maier, Richard Gerum, Claus Metzner, and Patrick Krauss.  Quantifying the separability of data classes in neural networks.  _Neural Networks_ , 139:278–293, July 2021.  ISSN 0893-6080.  doi: 10.1016/j.neunet.2021.03.035.  URL http://dx.doi.org/10.1016/j.neunet.2021.03.035.
  * Ferrari (2019) Alberto Ferrari.  A note on sum and difference of correlated chi-squared variables, 2019\.
  * Biroli and Mézard (2023) Giulio Biroli and Marc Mézard.  Generative diffusion in very large dimensions.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2023(9):093402, September 2023.  ISSN 1742-5468.  doi: 10.1088/1742-5468/acf8ba.  URL http://dx.doi.org/10.1088/1742-5468/acf8ba.
  * Tiomoko et al. (2019) Malik Tiomoko, Romain Couillet, Eric Moisan, and Steeve Zozor.  Improved estimation of the distance between covariance matrices.  In _ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 7445–7449, 2019.  doi: 10.1109/ICASSP.2019.8682621.
  * Zhang et al. (2018) Guodong Zhang, Chaoqi Wang, Bowen Xu, and Roger B. Grosse.  Three mechanisms of weight decay regularization.  _CoRR_ , abs/1810.12281, 2018.  URL http://arxiv.org/abs/1810.12281.
  * Levi and Oz (2023a) Noam Levi and Yaron Oz.  The underlying scaling laws and universal statistical structure of complex datasets.  _arXiv preprint arXiv:2306.14975_ , 2023a.
  * Loureiro et al. (2021d) Bruno Loureiro, Cedric Gerbelot, Hugo Cui, Sebastian Goldt, Florent Krzakala, Marc Mezard, and Lenka Zdeborova.  Learning curves of generic features maps for realistic datasets with a teacher-student model.  In _Advances in Neural Information Processing Systems_ , volume 34, 2021d.
  * Loureiro et al. (2022) Bruno Loureiro, Cedric Gerbelot, Hugo Cui, Sebastian Goldt, Florent Krzakala, Marc Mezard, and Lenka Zdeborova.  Learning curves of generic features maps for realistic datasets with a teacher-student model.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2022(11):114001, nov 2022.  doi: 10.1088/1742-5468/ac9825.  URL https://doi.org/10.1088%2F1742-5468%2Fac9825.
  * Radhakrishnan et al. (2023) Adityanarayanan Radhakrishnan, Mikhail Belkin, and Caroline Uhler.  Wide and deep neural networks achieve consistency for classification.  _Proceedings of the National Academy of Sciences_ , 120(14), March 2023.  ISSN 1091-6490.  doi: 10.1073/pnas.2208779120.  URL http://dx.doi.org/10.1073/pnas.2208779120.
  * Jaganathan et al. (2015) Kishore Jaganathan, Yonina C. Eldar, and Babak Hassibi.  Phase retrieval: An overview of recent developments, 2015.
  * Dong et al. (2023) Jonathan Dong, Lorenzo Valzania, Antoine Maillard, Thanh-an Pham, Sylvain Gigan, and Michael Unser.  Phase retrieval: From computational imaging to machine learning: A tutorial.  _IEEE Signal Processing Magazine_ , 40(1):45–57, 2023.  doi: 10.1109/MSP.2022.3219240.
  * Candes et al. (2014) Emmanuel J. Candes, Xiaodong Li, and Mahdi Soltanolkotabi.  Phase retrieval via wirtinger flow: Theory and algorithms.  _CoRR_ , abs/1407.1065, 2014.  URL http://arxiv.org/abs/1407.1065.
  * Chen et al. (2019) Yuxin Chen, Yuejie Chi, Jianqing Fan, and Cong Ma.  Gradient descent with random initialization: fast global convergence for nonconvex phase retrieval.  _Mathematical Programming_ , 176(1–2):5–37, February 2019.  ISSN 1436-4646.  doi: 10.1007/s10107-019-01363-6.  URL http://dx.doi.org/10.1007/s10107-019-01363-6.
  * Arnaboldi et al. (2024) Luca Arnaboldi, Florent Krzakala, Bruno Loureiro, and Ludovic Stephan.  Escaping mediocrity: how two-layer networks learn hard generalized linear models with sgd, 2024.
  * Martin et al. (2023) Simon Martin, Francis Bach, and Giulio Biroli.  On the impact of overparameterization on the training of a shallow neural network in high dimensions, 2023.
  * Livni et al. (2014) Roi Livni, Shai Shalev-Shwartz, and Ohad Shamir.  On the computational efficiency of training neural networks, 2014.
  * Soltani and Hegde (2018) Mohammadreza Soltani and Chinmay Hegde.  Towards provable learning of polynomial neural networks using low-rank matrix estimation.  In Amos Storkey and Fernando Perez-Cruz, editors, _Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics_ , volume 84 of _Proceedings of Machine Learning Research_ , pages 1417–1426. PMLR, 09–11 Apr 2018.  URL https://proceedings.mlr.press/v84/soltani18a.html.
  * Ji and Telgarsky (2020) Ziwei Ji and Matus Telgarsky.  Directional convergence and alignment in deep learning, 2020.
  * Lyu and Li (2020) Kaifeng Lyu and Jian Li.  Gradient descent maximizes the margin of homogeneous neural networks, 2020\.
  * Hsu et al. (2022) Daniel Hsu, Vidya Muthukumar, and Ji Xu.  On the proliferation of support vectors in high dimensions, 2022.
  * Krizhevsky (2012) Alex Krizhevsky.  Learning multiple layers of features from tiny images.  _University of Toronto_ , 05 2012.  URL https://www.cs.toronto.edu/~kriz/learning-features-2009-TR.pdf.
  * Xiao et al. (2017) Han Xiao, Kashif Rasul, and Roland Vollgraf.  Fashion-mnist: A novel image dataset for benchmarking machine learning algorithms.  _arXiv: 1708.07747_ , 2017.  doi: 10.48550/ARXIV.1708.07747.
  * He et al. (2015) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition, 2015.
  * Belhumeur et al. (1997) Peter N. Belhumeur, João P. Hespanha, and David J. Kriegman.  Eigenfaces vs. fisherfaces: Recognition using class specific linear projection.  _IEEE Trans. Pattern Anal. Mach. Intell._ , 19(7):711–720, jul 1997.  ISSN 0162-8828.  doi: 10.1109/34.598228.  URL https://doi.org/10.1109/34.598228.
  * Levi and Oz (2023b) Noam Levi and Yaron Oz.  The universal statistical structure and scaling laws of chaos and turbulence.  _arXiv preprint arXiv:2311.01358_ , 2023b.
  * Papyan et al. (2020) Vardan Papyan, X. Y. Han, and David L. Donoho.  Prevalence of neural collapse during the terminal phase of deep learning training.  _Proceedings of the National Academy of Sciences_ , 117(40):24652–24663, September 2020.  ISSN 1091-6490.  doi: 10.1073/pnas.2015509117.  URL http://dx.doi.org/10.1073/pnas.2015509117.

##  Appendix A Neural Collapse on Overlapping GMM Data

Here, we explore the concept of Neural Collapse (NC) on overlapping GMMs
generated from real data. Practically, we employ the metrics given in Papyan
et al. [2020] to measure the NC when training a ResNet18 on gaussian versions
of CIFAR10 images. In Figure 5, we show that NC occurs on overlapping GMMs
even when the class mean vector of the data has been set to zero, indicating
that expressive networks are able to perform nonlinear transformations that
generate class means, where none were present, and at the final layer still
perform linear classification.

Various metrics of Neural Collapse, as given by Papyan et al. [2020], are
shown for GMMs generated from two classes of CIFAR10. Here, we use 5000
samples per class, training a ResNet18 to perform binary classification using
cross-entropy loss, on GMMs generated from greyscale versions of CIFAR10
images. Each sample has dimension $1024$, and the training is done using the
default specifications provided in the supplementary code of Papyan et al.
[2020]

Figure 5: Results for neural collapse convergence for training ResNet18 on two
classes (0,7) sampled from Gaussian versions of CIFAR10 data. Left: setting
$\mu_{A}=\mu_{B}=0$. Right: $\mu_{A}\neq\mu_{B}$, as given by the class means.
The network attains 99.9% test accuracy in both cases. We see no substantial
difference in convergence to neural collapse metrics, affirming our claims in
the main text.

##  Appendix B KKT formulation

As discussed in the main text, we posit that the approach towards a BOC for a
NN can be explained by convergence to a KKT point.

Our reasoning follows Theorem Theorem B.1 below, which holds for _gradient
flow_ (i.e., gradient descent with an infinitesimally small step size). Before
stating the theorem, we need the following definitions: (1) We say that
gradient flow converges in direction to $\tilde{{\boldsymbol{\theta}}}$ if
$\lim_{t\to\infty}\frac{{\boldsymbol{\theta}}(t)}{\|{\boldsymbol{\theta}}(t)\|}=\frac{\tilde{{\boldsymbol{\theta}}}}{\|\tilde{{\boldsymbol{\theta}}}\|}$,
where ${\boldsymbol{\theta}}(t)$ is the parameter vector at time $t$; (2) We
say that a network $\Phi$ is _homogeneous_ w.r.t. the parameters
${\boldsymbol{\theta}}$ if there exists $L>0$ such that for every $\lambda>0$
and ${\boldsymbol{\theta}},\mathbf{x}$ we have
$\Phi(\alpha{\boldsymbol{\theta}};\mathbf{x})=\lambda^{L}\Phi({\boldsymbol{\theta}};\mathbf{x})$.
Thus, scaling the parameters by any factor $\lambda>0$ scales the outputs by
$\lambda^{L}$. We note that essentially any fully-connected or convolutional
neural network with ReLU activations is homogeneous w.r.t. the parameters
${\boldsymbol{\theta}}$ if it does not have any skip-connections (i.e.,
residual connections) or bias terms, except possibly for the first layer.

######  Theorem B.1 (Paraphrased from Lyu and Li [2020], Ji and Telgarsky
[2020]).

Let $\Phi({\boldsymbol{\theta}};\cdot)$ be a homogeneous ReLU neural network.
Consider minimizing the logistic loss over a binary classification dataset
$\\{(\mathbf{x}_{a},y_{a})\\}_{a=1}^{N}$ using gradient flow. Assume that
there exists time $t_{0}$ such that ${\cal
L}({\boldsymbol{\theta}}(t_{0}))<1$333This ensures that
$\ell(y_{a}\Phi({\boldsymbol{\theta}}(t_{0});\mathbf{x}_{a}))<1$ for all $i$,
i.e. at some time $\Phi$ classifies every sample correctly.. Then, gradient
flow converges in direction to a first order stationary point (KKT point) of
the following maximum-margin problem:

$\min_{{\boldsymbol{\theta}}^{\prime}}\frac{1}{2}\|{\boldsymbol{\theta}}^{\prime}\|^{2}\;\;\;\;\text{s.t.
}\;\;\;\forall
i\in[n]\;\;y_{a}\Phi({\boldsymbol{\theta}}^{\prime};\mathbf{x}_{a})\geq 1~{}.$
(19)

Moreover, ${\cal L}({\boldsymbol{\theta}}(t))\to 0$ as $t\to\infty$.

The above theorem guarantees directional convergence to a first order
stationary point (of the optimization problem (Equation 19)), which is also
called _Karush–Kuhn–Tucker point_ , or _KKT point_ for short. The KKT approach
allows inequality constraints, and is a generalization of the method of
_Lagrange multipliers_ , which allows only equality constraints.

The great virtue of Theorem B.1 is that it characterizes the _implicit bias_
of gradient flow with the logistic loss for homogeneous networks. Namely, even
though there are many possible directions of
$\frac{{\boldsymbol{\theta}}}{\|{\boldsymbol{\theta}}\|}$ that classify the
dataset correctly, gradient flow converges only to directions that are KKT
points of Problem (Equation 19). In particular, if the trajectory
${\boldsymbol{\theta}}(t)$ of gradient flow under the regime of Theorem B.1
converges in direction to a KKT point $\tilde{{\boldsymbol{\theta}}}$, then we
have the following: There exist
$\lambda_{1},\ldots,\lambda_{n}\in{\mathbb{R}}$ such that

$\displaystyle\tilde{{\boldsymbol{\theta}}}=\sum_{a=1}^{N}\lambda_{a}y_{a}\nabla_{\boldsymbol{\theta}}\Phi(\tilde{{\boldsymbol{\theta}}};\mathbf{x}_{a})~{}$
(stationarity) (20) $\displaystyle\forall
a\in[N],\;\;y_{a}\Phi(\tilde{{\boldsymbol{\theta}}};\mathbf{x}_{a})\geq 1$
(primal feasibility) (21) $\displaystyle\lambda_{1},\ldots,\lambda_{N}\geq 0$
(dual feasibility) (22) $\displaystyle\forall
a\in[N],~{}~{}\lambda_{a}=0~{}\text{if}~{}y_{a}\Phi(\tilde{{\boldsymbol{\theta}}};\mathbf{x}_{a})\neq
1$ (complementary slackness) (23)

Our main insight is based on Equation 20, which implies that the parameters
$\tilde{{\boldsymbol{\theta}}}$ are a linear combinations of the derivatives
of the network at the training data points. We say that a data point
$\mathbf{x}_{a}$ is _on the margin_ if
$y_{a}\Phi(\tilde{{\boldsymbol{\theta}}};\mathbf{x}_{a})=1$ (i.e.
$|\Phi(\tilde{{\boldsymbol{\theta}}};\mathbf{x}_{a})|=1$) . Note that Equation
23 implies that only samples which are on the margin affect Equation 20, since
samples not on the margin have a coefficient $\lambda_{a}=0$. In sufficiently
high dimension we expect no samples to be off the margin in GMM classification
therefore all $\lambda_{a}$ should be nonzero and contribute a similar amount.

##  Appendix C Dataset Properties

We show the eigenvalues scaling structure and the Inverse Participation Ratio
(IPR) of the eigenvectors. We separate the eigenvalues to three parts: I. The
large eigenvalues, II. The bulk eigenvalues, and III. The small eigenvalues,
as can be seen in Figure 6. These three regimes are seen also in real datasets
as shown in previous work [Levi and Oz, 2023a, b]. Figure 6 shows the
eigenvalues for each of the classes in the FMNIST and CIFAR10 datasets,
respectively.

Figure 6: Left: Partition the eigenvector indices to three regimes: the large
eigenvalues, the bulk and the small eigenvalues. The top figure shows the
eigenvalues versus their index, and the bottom show the corresponding IPR
versus the eigenvalues index. Center/Right: The eigenvalues of the covariance
matrix for the different classes in CIFAR10/FMNIST dataset, with a comparison
to uncorrelated data (blue). The orange line represents the eigenvalues of the
covariance matrix of the whole dataset and the green dashed line shows the
average of the different classes. We see the three regimes of eigenvalues: the
large eigenvalues - up to $\pm 10$ first eigenvalues, the bulk scaling -
between $\pm 10$ to $\pm 300$, and the small eigenvalues - from $\pm 350$.

The IPR defined as:

$\displaystyle IPR\equiv\left(\sum_{n}\left|v_{n}\right|^{4}\right)^{-1},$
(24)

is used to measure the localization rate (entropy) of a vector. Figure 7 shows
the IPR of the eigenvectors of the covariance matrices for the different
classes in the CIFAR10 dataset and FMNIST. As can be seen from the Figure, the
changes in the eigenvalues structure and the IPR is small.

Figure 7: The histogram of the IPR of the covariance matrix eigenvectors for
the different classes in the CIFAR10 dataset (left), and FMNIST (right), with
a comparison to noise (blue).

##  Appendix D Appendix - Flipping tests examples

We include the different figures concerning the classification flipping tests
referred to in the main text. The diagram in Figure 8 outlines the structure
of the flipping tests. In Figure 9 and Figure 10 we see rotated images of
FMNIST and CIFAR10, respectively. In Figure 11, Figure 12 and Figures 13 and
14 we plot the eigenvalues and eigenvectors thresholds.

The code for our experiments is provded in
https://github.com/khencohen/FlippingsTests.

Figure 8: The classification flipping tests. We evaluate the model
classification performance on different covariance matrices. Left - building a
new covariance matrix by flipping the eigenvectors up to a certain threshold.
Right - building a new covariance matrix by flipping the eigenvalues up to a
certain threshold. This method provides information about the relative
importance of the different eigenvectors and eigenvalues for the model to make
a classification decision. Figure 9: Examples of rotated images in the FMNIST
dataset. Each of the samples of class 0 (shirt, top row), was rotated using
the covariance matrix of another class (bottom row), and the middle row shows
the result of the rotation. In all of our tests, for any rotation of class 0
sample, with class $c$ basis, the model classified the rotated image as a
sample from class $c$. Figure 10: Examples of rotated images in the CIFAR
dataset. Each of the samples of class 0 (airplane, top row), was rotated using
the covariance matrix of another class (bottom row), and the middle row shows
the result of the rotation. In all of our tests, for any rotation of class 0
sample, with class $c$ basis, the model classified the rotated image as a
sample from class $c$. Figure 11: The classification eigenvector (left) and
eigenvalues (right) flipping test on the optimal classifier. CIFAR10 dataset,
comparing class 0 with other classes.

Figure 12: The classification eigenvector flipping test on the Gaussian
dataset between class 0 and class 5, for different number of samples. FC model
(left) and ResNet18 (right), CIFAR10 dataset.

Figure 13: The classification eigenvalues flipping test on the rotated FMNIST
dataset between class 0 and other classes for the FC model.

Figure 14: The classification eigenvalues flipping test on the rotated FMNIST
dataset between class 0 and other classes for the ResNet 18 model.


and read/RMW labels agree with current buffer and memory contents. The first
condition holds because propagation occurs only after write transition that
puts an item into a buffer, and propagation is the only way to clear a buffer.
The second condition holds because enumeration respects the program order
between writes and RMW, so by the moment of such transition’s execution all
previous thread’s writes will be already propagated.

To prove the latter condition we’ll consider the write event that a given read
event reads from in $G$ and show that this write coincides with the one
determined by the $\mathsf{TSO}$ memory subsystem999Here we reason about
events in the context of operational semantics. Formally, it only deals with
event labels. But every value stored in memory or buffer is placed there due
to some transition execution which in turn is determined by some graph event,
since the trace under consideration is built from $G$. So we implicitly assume
such mapping between memory values and graph events that produce them. . Let
$r$ be the read event and $loc$ the location it reads from, $w_{tr}$ the write
that has been observed during the transition and
$w_{{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}}$
the
${\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}$-predecessor
of $r$. Suppose by contradiction that $w_{tr}\neq
w_{{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}}$.

  * •

If the thread buffer has non-propagated writes to the location $loc$, the
$w_{tr}$ will be the latest non-propagated write to $loc$ of the same thread.

    * –

If
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}(w_{{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}},w_{tr})$,
then
${\color[rgb]{.75,0,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,0,.25}\mathtt{fr}}(r,w_{tr})$.
Therefore, there is an
${\color[rgb]{0.3,0.3,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0.3,0.3,0.3}\pgfsys@color@gray@stroke{0.3}\pgfsys@color@gray@fill{0.3}\mathtt{po}};{\color[rgb]{.75,0,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,0,.25}\mathtt{fr}}$
cycle prohibited by $\mathsf{TSO}$.

    * –

Suppose that
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}(w_{tr},w_{{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}})$.

      * *

If
$w_{{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}}$
belongs to the same thread, it can only appear after $r$ in the program order,
since $w_{tr}$ is the latest write. So there is a
${\color[rgb]{0.3,0.3,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0.3,0.3,0.3}\pgfsys@color@gray@stroke{0.3}\pgfsys@color@gray@fill{0.3}\mathtt{po}};{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}$
cycle prohibited by $\mathsf{TSO}$.

      * *

If
$w_{{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}}$
belongs to the other thread, it has been already been propagated since
${\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}\mathtt{e}$
is respected by the enumeration. But it means that $w_{tr}$ is also
propagated.

  * •

If all previous writes to $loc$ from the same thread have been already
propagated, the $w_{tr}$ will be the
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$-latest
(in the enumeration) write to $loc$.

    * –

If
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}(w_{{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}},w_{tr})$,
then
${\color[rgb]{.75,0,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,0,.25}\mathtt{fr}}(r,w_{tr})$.
Then, since the enumeration respects
${\color[rgb]{.75,0,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,0,.25}\mathtt{fr}}$,
$w_{tr}$ will be preceded by $r$ in the trace order and thus cannot be
observed during the transition.

    * –

Suppose that
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}(w_{tr},w_{{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}})$.
Then
$w_{{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}}$
is not propagated yet. The only case when $r$ can still read from it is when
$w_{{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}}$
belongs to the same thread and is buffered. But it contradicts the condition
on the empty buffer. ∎

This concludes the proof of Lemma A.22. ∎

##  Appendix B Termination Checking

Just like termination, non-termination can be proven by considering only the
iteration of a spinloop which reads only from
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$-maximal
stores. If in a finite graph some threads are blocked by spinloop iterations
which read from only
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$-maximal
stores, and all other threads have terminated, then the graph can be extended
to a fair infinite graph with infinite spinloops.

###### Definition B.1.

A thread ${\tau}$ is terminated in execution graph $G$ if in the (unique)
behavior $\beta$ with $G.{\mathtt{E}}=\mathsf{Event}(\beta)$,
$\mu_{\tau}(\beta)$ is terminating.

An execution graph $G$ is a non-termination witness if there is a thread which
is not terminated in $G$ and all threads ${\tau}$ which are not terminated in
$G$ have finite runs $\mu_{\tau}(\beta)$ where $\beta$ is the (unique)
behavior with $G.{\mathtt{E}}=\mathsf{Event}(\beta)$, and these runs end in a
spinloop iteration reading from
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$-maximal
stores, i.e., for all ${\tau}$ there is a spinloop iteration
$[n_{\tau},{n_{\tau}}^{\prime}]$ of ${\tau}$ in $\beta$ where
${n_{\tau}}^{\prime}$ is the length of $\mu_{\tau}(\beta)$ and

$[\\{{e\in
G.{\mathtt{E}}^{\tau}\;|\;n_{\tau}\leq{\mathtt{sn}}(e)\leq{n_{\tau}}^{\prime}}\\}];G.{\color[rgb]{.75,0,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,0,.25}\mathtt{fr}}=\emptyset$

Every non-termination witness proves the existence of a fair graph with
infinite spinloops.

###### Theorem B.2.

Let $G$ be a non-termination witness with
$G.{\mathtt{E}}=\mathsf{Event}(\beta)$ where $\beta$ is a behavior of a
deterministic program. Then there is a fair execution graph $G^{\prime}$ with
$G^{\prime}.{\mathtt{E}}=\mathsf{Event}(\beta^{\prime})$ where
$\beta^{\prime}$ is a behavior of the program and every thread which is not
terminated in $G$ has an infinite spinloop in $\beta^{\prime}$.

###### Proof.

We first define $\beta^{\prime}$. For each ${\tau}$ we defined
$\beta^{\prime}({\tau})$ by case analysis on whether ${\tau}$ is terminated in
$G$. If ${\tau}$ is terminated, then $\beta^{\prime}({\tau})=\beta({\tau})$.
Otherwise, since $G$ is a non-termination witness, the run $\mu_{\tau}(\beta)$
of ${\tau}$ is finite. Let $K_{\tau}$ be the length of this run. For steps
$k<K_{\tau}$ inside this run, the behavior is unchanged

$k<K_{\tau}\ \implies\ \beta^{\prime}({\tau},k)=\beta({\tau},k)$

For later steps, we repeat the last spinloop iteration infinitely often. Let
$[n_{\tau},{n_{\tau}}^{\prime}]$ be the final spinloop iteration of ${\tau}$
in $\beta$ which reads from
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$-maximal
stores, and $L_{\tau}={n_{\tau}}^{\prime}-n_{\tau}+1$ the length of that
iteration (i.e., the number of transitions in the iteration). Then the
infinite repetition of this spinloop iteration is defined by

$k\geq K_{\tau}\ \implies\
\beta^{\prime}({\tau},k)=\beta({\tau},K_{\tau}-L_{\tau}+(k-K_{\tau})\mod
L_{\tau})$

Due to determinism, each of these spinloop iterations in the program reads
exactly the same values from the same locations and hence traverses exactly
the same states. This shows that $\beta^{\prime}$ is a behavior of the
program.

We proceed to define $G^{\prime}$. The set of events is defined by
$G^{\prime}.{\mathtt{E}}=\mathsf{Event}(\beta^{\prime})$. Modification order
is unchanged
($G^{\prime}.{\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}=G.{\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$).
For the reads-from relation we let each of the new events read from the
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$-maximal
store that the final spinloop iteration read from. We first define a mapping
$l({\tau},k)$ which maps each serial number $k$ of an event in $G^{\prime}$ to
the serial number of the corresponding event in $G$

$l({\tau},k)=\begin{cases}k&{\tau}\ \text{is terminated}\ \lor\ k<K_{\tau}\\\
K_{\tau}-L_{\tau}+(k-K_{\tau})\mod L_{\tau}&\text{o.w.}\end{cases}$

With this definition, we have
$\beta^{\prime}({\tau},k)=\beta({\tau},l({\tau},k))$. We lift this mapping to
events in the obvious way:

$l(e)={\langle{l({\mathtt{tid}}(e),{\mathtt{sn}}(e)),{\mathtt{tid}}(e):{\mathtt{elab}}(e)}\rangle}$

Observe that we have $l(e)\in G.{\mathtt{E}}$ for all $e\in
G^{\prime}.{\mathtt{E}}$. Furthermore $G.{\mathtt{E}}\subseteq
G^{\prime}.{\mathtt{E}}$; hence $l$ is a relation on
$G^{\prime}.{\mathtt{E}}$. It relates new events to the corresponding old
event. Its inverse $l^{-1}$ relates old events to the corresponding new
events. This allows defining the reads-from relation of $G^{\prime}$ as

$G^{\prime}.{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}\
=\
G.{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}};l^{-1}$

Since
$G^{\prime}.{\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$
is
$G.{\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$,
it remains prefix-finite. Furthermore, the
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$-maximal
stores in $G^{\prime}$ are exactly the
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$-maximal
stores in $G$. Since all the new reads added read from
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$-maximal
stores, they have no
${\color[rgb]{.75,0,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,0,.25}\mathtt{fr}}$-successors.
Thus
${\color[rgb]{.75,0,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,0,.25}\mathtt{fr}}$
remains also prefix finite. Therefore $G^{\prime}$ is fair, and the claim
follows. ∎

This infinite extension of the graph does not in general preserve
$\mathcal{G}$-consistency. For example, consider an (absurd) declarative
memory system $\mathcal{G}$ which only admits finite execution graphs. We
capture declarative memory systems which allow the infinite extension by the
following predicate.

###### Definition B.3.

A declarative memory system $\mathcal{G}$ is infinitely-read-extensible if
adding infinitely many reads to
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$-maximal
elements preserves $\mathcal{G}$-consistency, i.e., if $G\in\mathcal{G}$ then
so is every $G^{\prime}$ with

$G^{\prime}.{\mathtt{E}}=G.{\mathtt{E}}\cup E$

with $E\subseteq{\mathtt{R}}$ satisfying

$[E];G.{\color[rgb]{.75,0,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,0,.25}\mathtt{fr}}=\emptyset$

Practically relevant declarative memory systems are all infinitely-read-
extensible. This holds in particular for
$\mathcal{G}_{\\{{\mathsf{SC},\mathsf{TSO},\mathsf{RA}}\\}}$.

###### Lemma B.4.

$\mathcal{G}_{\\{{\mathsf{SC},\mathsf{TSO},\mathsf{RA}}\\}}\quad\text{are
infinitely-read-extensible}$

###### Proof.

Omitted. ∎

###### Theorem B.5.

If $\mathcal{G}$ is infinitely-read-extensible and $G$ is
$\mathcal{G}$-consistent, so is $G^{\prime}$ (where $G$, $G^{\prime}$ are the
graphs from Theorem B.2).

###### Proof.

By assumption, the events in the final spinloop iterations read from
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$-maximal
events. The newly added events

$E=G^{\prime}.{\mathtt{E}}-G.{\mathtt{E}}$

are by definition all reading from those same events. Observe first that the
newly added events are all mapped by $l$ to events in spinloop iterations

$[E];l=[\\{{e\in
G.{\mathtt{E}}\;|\;\exists{\tau}.\;n_{\tau}\leq{\mathtt{sn}}(e)\leq{n_{\tau}}^{\prime}}\\}]$

Since $G$ is a non-termination witness by assumption, these events read from
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$
maximal stores

$[\\{{e\in
G.{\mathtt{E}}\;|\;\exists{\tau}.\;n_{\tau}\leq{\mathtt{sn}}(e)\leq{n_{\tau}}^{\prime}}\\}];G.{\color[rgb]{.75,0,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,0,.25}\mathtt{fr}}=\emptyset$

We conclude by definition of $G^{\prime}$ that the newly added events $E$ read
from
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$-maximal
stores:

$\displaystyle[E];G^{\prime}.{\color[rgb]{.75,0,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,0,.25}\mathtt{fr}}\
$ $\displaystyle=\
[E];G^{\prime}.{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}^{-1};G^{\prime}.{\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$
$\displaystyle=\
[E];(G.{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}};l^{-1})^{-1};G.{\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$
$\displaystyle=\
[E];l;G.{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\mathtt{rf}}^{-1};G.{\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$
$\displaystyle=\
[E];l;G.{\color[rgb]{.75,0,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,0,.25}\mathtt{fr}}$
$\displaystyle=\ [\\{{e\in
G.{\mathtt{E}}\;|\;\exists{\tau}.\;n_{\tau}\leq{\mathtt{sn}}(e)\leq{n_{\tau}}^{\prime}}\\}];G.{\color[rgb]{.75,0,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,0,.25}\mathtt{fr}}$
$\displaystyle=\ \emptyset$

Since $\mathcal{G}$ is infinitely-read-extensible, the claim follows. ∎

###### Lemma B.6.

If $G$ from Theorem B.2 is
$\mathcal{G}_{\\{{\mathsf{SC},\mathsf{TSO},\mathsf{RA}}\\}}$-consistent, so is
$G^{\prime}$.

###### Proof.

Follows directly from Theorems B.5 and B.4. ∎

Using Theorems B.5 and 6.3, we can use model checking to automatically decide
the termination of spinloops for any client which does not generate infinite
modification orders. The model checker must simply generate all finite
executions of a program in which spinloops are executed for at most one
iteration. This enumerates all non-termination witnesses (modulo staying in
spinloops longer than necessary). If any non-termination witness exists, then
by Theorem B.5 there is an infinite, fair extension of the graph in which the
spinloop never terminates. If there is no non-termination witness, then all
spinloops must terminate after reading from
${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\mathtt{mo}}$-maximal
stores. Thus by Theorem 6.3, all spinloops terminate.

In (Oberhauser et al., 2021a) we have implemented this scheme as an extension
of GenMC (Kokologiannakis et al., 2019), a model checker which can be
parametrized with declarative memory systems satisfying certain conditions
(MM1–MM4 in (Kokologiannakis et al., 2019)). We applied the extended model
checker to several locks. Among others, the MCS lock, waiter-count based
locks, and the Linux qspinlock can get stuck in a spinloop when using
incorrect barrier combinations. Many of these barrier combinations ensure
safety of the locks. Since non-termination problems under weak consistency can
not be detected by existing model checkers, and these barrier combinations are
otherwise correct, existing model checkers may give verification engineers an
unwarranted level of confidence for programs that have spinloops.