# Quantum baryon number fluctuations in subsystems of a hot and dense
relativistic gas of fermions

Arpan Das1 Wojciech Florkowski2 Radoslaw Ryblewski1 Rajeev Singh1 1Institute
of Nuclear Physics Polish Academy of Sciences, PL-31-342 Kraków, Poland
2Institute of Theoretical Physics, Jagiellonian University, PL-30-348 Kraków,
Poland

( )

###### Abstract

Quantum features of the baryon number fluctuations in subsystems of a hot and
dense relativistic gas of fermions are analyzed. We find that the fluctuations
in small systems are significantly increased compared to their values known
from the statistical physics, and diverge in the limit where the system size
goes to zero. The numerical results obtained for a broad range of the
thermodynamic parameters expected in heavy-ion collisions are presented. They
can be helpful to interpret and shed new light on the experimental data.

##  I Introduction

Statistical fluctuations in many-body systems (microscopic as well as
macroscopic ones) play a very important role in physics as they encode the
crucial information about possible phase transitions, dissipation, and
clustering phenomena Smoluchowski ; PhysRevLett.85.2076 ; Huang:1987asp ;
Kubo1 ; Lifshitz:1963ps ; PhysRevLett.49.1110 . An unexplored novel feature of
the fluctuations is their increase for small systems in the cases where the
quantum effects become important. Such effects have been analyzed
quantitatively in our two recent papers Das:2020ddr ; Das:2021aar  where we
addressed the fluctuations of the energy density in a hot gas of bosons and
fermions. Our results indicate limitations of the concept of the fluid
elements used in relativistic hydrodynamics applied in the description of
heavy-ion collisions. As the size of a subsystem drops below about 0.5 fm, the
energy density fluctuations (for typical values of temperatures and particle
masses) become so large that they are comparable with their mean values. In
this case, the physical picture of fluid cells with well defined energy
density becomes unjustified. An interesting feature of our calculations
Das:2020ddr ; Das:2021aar  is that not only fluctuations diverge for the
subsystem with size approaching zero, but also quantum statistical
fluctuations agree with thermodynamic fluctuations if the subsystems become
sufficiently large Das:2020ddr ; Das:2021aar . In this way, our approach
incorporates quantum features into statistical mechanics in a very natural
way.

In the present work, following similar ideas to that developed in Das:2020ddr
; Das:2021aar , we discuss fluctuations of the baryon number density in a hot
and dense relativistic gas of fermions. Our analysis is relevant for
relativistic heavy-ion physics, in particular, in the context of the beam
energy scan (BES). Hunt for the conjectured critical endpoint in the QCD phase
diagram has triggered vast theoretical and experimental studies of many
fluctuation observables. For example, one studies fluctuations of the
conserved quantum numbers in QCD, such as the baryon number, electric charge,
or strangeness. They all provide an excellent opportunity to study the
critical phenomena Berges:1998rc ; Halasz:1998qr ; Stephanov:1998dy ;
Stephanov:1999zu ; Hatta:2002sj ; Stephanov:2008qz ; Berdnikov:1999ph ;
Kitazawa:2013bta .

Following our previous framework Das:2020ddr ; Das:2021aar , we consider the
fluctuation of the baryon number in the subsystem $S_{a}$ of the thermodynamic
system $S_{V}$ described by the grand canonical ensemble characterized by the
temperature $(T)$ (or its inverse $\beta=1/T$) and the baryon chemical
potential ($\mu$). The volume $V$ of the larger system $S_{V}$ is larger than
the characteristic volume of the subsystem $S_{a}$. We derive a compact
formula that defines quantum fluctuations of the baryon number operator in
subsystems of a hot and dense relativistic gas. Then we apply this formula to
get physical insights into situations expected in relativistic heavy-ion
collisions. In the calculation we use the metric tensor with the signature
$(+1,-1,-1,-1)$. Three-vectors are shown in bold font and a dot is used to
denote the scalar product of both four and three-vectors, i.e.,
$a^{\mu}b_{\mu}=a\cdot b=a^{0}b^{0}-\bm{a}\cdot\bm{b}$.

The structure of the paper is organized as follows: We start with introducing
the basic concepts and definitions in Sec. II, where we also derive the
expression for baryon number density fluctuation. Subsequently, in Sec. III we
discuss the thermodynamic limit of baryon number density fluctuation. The
numerical results of the baryon number density fluctuations are presented in
Sec. IV. Finally, in Sec. V we interpret and summarize our results.

##  II Basic concepts and definitions

We consider a system of spin-$\nicefrac{{1}}{{2}}$ particles described by the
quantum Dirac field in thermal equilibrium. The field operator is expressed in
the standard way as Tinti:2020gyh

$\displaystyle\psi(t,\bm{x})=$
$\displaystyle\sum_{r}\int\frac{dK}{\sqrt{2\omega_{\bm{k}}}}\Big{(}U_{r}^{{\phantom{\dagger}}}(\bm{k})a_{r}^{{\phantom{\dagger}}}(\bm{k})e^{-ik\cdot
x}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+V_{r}^{{\phantom{\dagger}}}(\bm{k})b_{r}^{\dagger}(\bm{k})e^{ik\cdot
x}\Big{)},$ (1)

where we use the notation $dK\equiv{d^{3}k}/{(2\pi)^{3}}$, while
$a_{r}^{{\phantom{\dagger}}}(\bm{k})$ and $b_{r}^{\dagger}(\bm{k})$ are the
annihilation and creation operators for particles and antiparticles,
respectively. The polarization degree of freedom is represented by the index
$r$. The fermionic operators $a_{r}^{{\phantom{\dagger}}}(\bm{k})$ and
$b_{r}^{\dagger}(\bm{k})$ satisfy the canonical anticommutation relations,
$\\{a_{r}^{{\phantom{\dagger}}}(\bm{k}),a_{s}^{\dagger}(\bm{k}^{\prime})\\}=(2\pi)^{3}\delta_{rs}\delta^{(3)}(\bm{k}-\bm{k}^{\prime})$
and
$\\{b_{r}^{{\phantom{\dagger}}}(\bm{k}),b_{s}^{\dagger}(\bm{k}^{\prime})\\}=(2\pi)^{3}\delta_{rs}\delta^{(3)}(\bm{k}-\bm{k}^{\prime})$.
All the other operators anticommute with each other. The Dirac spinors
$U_{r}^{{\phantom{\dagger}}}(\bm{k})$ and
$V_{r}^{{\phantom{\dagger}}}(\bm{k})$ have normalization
${\bar{U}}_{r}^{{\phantom{\dagger}}}(\bm{k})U_{s}^{{\phantom{\dagger}}}(\bm{k})=2m\delta_{rs}$
and
${\bar{V}}_{r}^{{\phantom{\dagger}}}(\bm{k})V_{s}^{{\phantom{\dagger}}}(\bm{k})=-2m\delta_{rs}$,
and the quantity $\omega_{\bm{k}}=\sqrt{\bm{k}^{2}+m^{2}}$ is the energy of a
particle.

To perform thermal averaging of quantum operators, it is sufficient to know
the thermal expectation values of the products of two and four creation and/or
annihilation operators (for both particles and antiparticles)
CohenTannoudji:422962 ; Itzykson:1980rh ; Evans:1996bha

$\displaystyle\langle
a_{r}^{\dagger}({\bm{k}})a_{s}^{{\phantom{\dagger}}}({\bm{k}}^{\prime})\rangle=(2\pi)^{3}\delta_{rs}\delta^{(3)}({\bm{k}}-{\bm{k}}^{\prime})f(\omega_{\bm{k}}),$
(2) $\displaystyle\langle
a^{\dagger}_{r}(\bm{k})a^{\dagger}_{s}(\bm{k}^{\prime})a_{r^{\prime}}^{{\phantom{\dagger}}}(\bm{p})a_{s^{\prime}}^{{\phantom{\dagger}}}(\bm{p}^{\prime})\rangle$
$\displaystyle=(2\pi)^{6}\Big{(}\delta_{rs^{\prime}}\delta_{r^{\prime}s}\delta^{(3)}(\bm{k}-\bm{p}^{\prime})~{}\delta^{(3)}(\bm{k}^{\prime}-\bm{p})$
$\displaystyle-\delta_{rr^{\prime}}\delta_{ss^{\prime}}\delta^{(3)}({\bm{k}}-\bm{p})~{}\delta^{(3)}({\bm{k}}^{\prime}-\bm{p}^{\prime})\Big{)}f(\omega_{{\bm{k}}})f(\omega_{{\bm{k}}^{\prime}}).$
(3)

Here $f(\omega_{{\bm{k}}})=1/(\exp(\beta(\omega_{\bm{k}}-\mu))+1)$ is the
Fermi–Dirac distribution function for particles. For antiparticles, the
Fermi–Dirac distribution function differs by the sign of the baryon chemical
potential $\mu$, i.e.
$\bar{f}(\omega_{{\bm{k}}})=1/(\exp(\beta(\omega_{\bm{k}}+\mu))+1)$. Using the
anticommutation relations for $a_{r}^{{\phantom{\dagger}}}({\bm{k}})$,
$a_{r}^{\dagger}({\bm{k}})$, $b_{r}^{{\phantom{\dagger}}}({\bm{k}})$, and
$b_{r}^{\dagger}({\bm{k}})$ any other combination of two and four creation
and/or annihilation operators can be obtained from Eqs. (2) and (3).

Following Ref. Chen:2018cts , we define the baryon number density operator
$\hat{J}^{0}_{a}$, associated with the conserved baryon current in a subsystem
$S_{a}$ using a smooth Gaussian profile placed at the origin of the coordinate
system, namely

$\displaystyle\hat{J}^{0}_{a}=\frac{1}{(a\sqrt{\pi})^{3}}\int
d^{3}\bm{x}~{}\hat{J}^{0}(x)~{}\exp\left(-\frac{{\bm{x}}^{2}}{a^{2}}\right),$
(4)

where $\hat{J}^{0}=\psi^{\dagger}\psi$. To remove any possible sharp-boundary
effects we consider the smooth profile with a length scale $a$ instead of a
cubic box. The thermal expectation value of the normal ordered operator
$:\hat{J}^{0}_{a}:$ is denoted as $\langle:\hat{J}^{0}_{a}:\rangle$. Here
standard normal ordering procedure has been introduced to remove an infinite
vacuum part coming from zero-point contributions. To determine the baryon
number fluctuation of the subsystem $S_{a}$, we consider the variance

$\sigma^{2}(a,m,T,\mu)=\langle:\hat{J}^{0}_{a}::\hat{J}^{0}_{a}:\rangle-\langle:\hat{J}^{0}_{a}:\rangle^{2}\,$
(5)

and the normalized standard deviation as

$\sigma_{n}(a,m,T,\mu)=\frac{(\langle:\hat{J}^{0}_{a}::\hat{J}^{0}_{a}:\rangle-\langle:\hat{J}^{0}_{a}:\rangle^{2})^{1/2}}{\langle:\hat{J}^{0}_{a}:\rangle}.$
(6)

Using the Fourier mode expansion of the Dirac field operator as given by Eq.
(1), the normal ordered operator $:\hat{J}^{0}_{a}:$ can be obtained as,

$\displaystyle:\hat{J}_{a}^{0}:=\sum_{r,s}\int\frac{dK}{\sqrt{2\omega_{\bm{k}}}}\frac{dK^{\prime}}{\sqrt{2\omega_{{\bm{k}}^{\prime}}}}\times$
$\displaystyle\bigg{[}a^{\dagger}_{s}(\bm{k}^{\prime})a^{{\phantom{\dagger}}}_{r}(\bm{k})\bar{U}^{{\phantom{\dagger}}}_{s}(\bm{k}^{\prime})\gamma^{0}U_{r}(\bm{k})~{}e^{-i(\omega_{\bm{k}}-\omega_{\bm{k}^{\prime}})t}~{}e^{-\frac{a^{2}}{4}(\bm{k}-\bm{k}^{\prime})^{2}}$
$\displaystyle-b^{\dagger}_{r}(\bm{k})b^{{\phantom{\dagger}}}_{s}(\bm{k}^{\prime})\bar{V}^{{\phantom{\dagger}}}_{s}(\bm{k}^{\prime})\gamma^{0}V_{r}(\bm{k})~{}e^{i(\omega_{\bm{k}}-\omega_{\bm{k}^{\prime}})t}~{}e^{-\frac{a^{2}}{4}(\bm{k}-\bm{k}^{\prime})^{2}}$
$\displaystyle+a^{\dagger}_{s}(\bm{k}^{\prime})b^{\dagger}_{r}(\bm{k})\bar{U}^{{\phantom{\dagger}}}_{s}(\bm{k}^{\prime})\gamma^{0}V_{r}(\bm{k})~{}e^{i(\omega_{\bm{k}}+\omega_{\bm{k}^{\prime}})t}~{}e^{-\frac{a^{2}}{4}(\bm{k}+\bm{k}^{\prime})^{2}}$
$\displaystyle+b^{{\phantom{\dagger}}}_{s}(\bm{k}^{\prime})a^{{\phantom{\dagger}}}_{r}(\bm{k})\bar{V}^{{\phantom{\dagger}}}_{s}(\bm{k}^{\prime})\gamma^{0}U_{r}(\bm{k})~{}e^{-i(\omega_{\bm{k}}+\omega_{\bm{k}^{\prime}})t}~{}e^{-\frac{a^{2}}{4}(\bm{k}+\bm{k}^{\prime})^{2}}\bigg{]}.$
(7)

Using the thermal averaging of two creation and/or annihilation operators as
given by Eq. (2), the thermal expectation value of $:\hat{J}_{a}^{0}:$ has the
form

$\displaystyle\langle:\hat{J}_{a}^{0}:\rangle=2\int
dK\Big{[}f(\omega_{\bm{k}})-\bar{f}(\omega_{\bm{k}})\Big{]}.$ (8)

This expression agrees with the standard kinetic-theory definition, with the
factor of 2 accounting for the spin degeneracy.

In the next step, using the thermal averaging of four creation and/or
annihilation operators as given by Eq. (3) we obtain

$\displaystyle\sigma^{2}(a,m,T,\mu)=\langle:\hat{J}_{a}^{0}::\hat{J}_{a}^{0}:\rangle-\langle:\hat{J}_{a}^{0}:\rangle^{2}$
$\displaystyle=\int\frac{dK}{\omega_{\bm{k}}}\frac{dK^{\prime}}{\omega_{\bm{k}^{\prime}}}(\omega_{\bm{k}}\omega_{\bm{k}^{\prime}}+\bm{k}\cdot\bm{k}^{\prime}+m^{2})e^{-\frac{a^{2}}{2}(\bm{k}-\bm{k}^{\prime})^{2}}\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left[f(\omega_{\bm{k}})\left(1-f(\omega_{\bm{k}^{\prime}})\right)+\bar{f}(\omega_{\bm{k}})(1-\bar{f}(\omega_{\bm{k}^{\prime}}))\right]$
$\displaystyle-\int\frac{dK}{\omega_{\bm{k}}}\frac{dK^{\prime}}{\omega_{\bm{k}^{\prime}}}(\omega_{\bm{k}}\omega_{\bm{k}^{\prime}}+\bm{k}\cdot\bm{k}^{\prime}-m^{2})e^{-\frac{a^{2}}{2}(\bm{k}+\bm{k}^{\prime})^{2}}\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left[f(\omega_{\bm{k}})(1-\bar{f}(\omega_{\bm{k}^{\prime}}))+\bar{f}(\omega_{\bm{k}})(1-f(\omega_{\bm{k}^{\prime}}))\right].$
(9)

Although normal ordering removes unwanted vacuum contribution to the baryon
number density operator, it is not sufficient to remove all the vacuum
divergences in all composite operators. Therefore following Refs. Das:2020ddr
; Das:2021aar  we discard a divergent temperature and chemical potential
independent vacuum term that originally appears in Eq. (9). Equation (9) is
the main result which represents a fluctuation of the baryon number density in
the subsystem $S_{a}$. Note that in Eqs. (8) and (9), the spin and particle-
antiparticle degrees of freedom are included. If we take into account the
degeneracy factors $(g)$ due to other internal degrees of freedom then we
should do the following replacements:
$\langle:\hat{J}_{a}^{0}:\rangle\rightarrow g\langle:\hat{J}_{a}^{0}:\rangle$
and $\sigma^{2}\rightarrow g\sigma^{2}$ Das:2020ddr ; Das:2021aar .

Figure 1: Variation of normalized fluctuation $\sigma_{n}$ in the subsystem
$S_{a}$ with the scale $a$ for different values of the temperature $T$ and
fixed particle mass $m=1.0$ GeV and baryon chemical potential $\mu=0.5$ GeV.
Figure 2: Variation of normalized fluctuation $\sigma_{n}$ in the subsystem
$S_{a}$ with the scale $a$ for different values of the baryon chemical
potential $\mu$ and fixed particle mass $m=1.0$ GeV and temperature $T=0.15$
GeV. Figure 3: Variation of normalized fluctuation $\sigma_{n}$ in the
subsystem $S_{a}$ with the scale $a$ for different values of the particle mass
and fixed temperature $T=0.15$ GeV and baryon chemical potential $\mu=0.5$
GeV.

##  III Thermodynamic limit

Before we turn to a discussion of our numerical results based on Eq. (9), it
is important to analyze the thermodynamic limit of $\sigma^{2}(a,m,T,\mu)$.
Since $S_{a}$ is a subsystem of the system $S_{V}$, the thermodynamic limit
can be achieved by considering the $a\rightarrow\infty$ limit (still with
$a^{3}\ll V$). In the thermodynamic limit quantum fluctuation as obtained here
should reduce to the expression for the classical statistical fluctuation. To
obtain the classical limit of the fluctuations we look into the
susceptibilities describing fluctuations in the baryon number obtained for
thermal and chemical equilibrium Nahrgang:2014fza

$\displaystyle\chi_{l}^{(B)}=\frac{\partial^{l}(P/T^{4})}{\partial(\mu/T)^{l}}\bigg{|}_{T}.$
(10)

Here $P$ denotes the thermodynamic pressure. Susceptibilities can also be
related to the cumulants of the distribution of baryons, e.g.,

$\displaystyle\chi_{1}^{(B)}=\frac{1}{VT^{3}}\langle
N_{B}\rangle=\frac{1}{T^{3}}\frac{\langle
N_{B}\rangle}{V}=\frac{n_{B}}{T^{3}},$
$\displaystyle\chi_{2}^{(B)}=\frac{1}{VT^{3}}\langle(\triangle
N_{B})^{2}\rangle=\frac{1}{VT^{3}}\langle(N_{B}-\langle
N_{B}\rangle)^{2}\rangle.$ (11)

Here $n_{B}\equiv N_{B}/V$ is the net baryon number density. Equation (11)
also implies that

$\displaystyle V\langle(n_{B}-\langle
n_{B}\rangle)^{2}\rangle=T^{3}\chi_{2}^{(B)}.$ (12)

Using Eq. (10), the susceptibility $\chi_{2}^{(B)}$ can also be obtained by
taking the appropriate derivative of the thermodynamic pressure. The
thermodynamic pressure at finite temperature and baryon chemical potential can
be expressed as Nahrgang:2014fza

$\displaystyle\frac{P}{T^{4}}=$ $\displaystyle\frac{2g}{T^{3}}\int
dK\Bigg{[}\ln\bigg(1+\exp\left(-\frac{\omega_{\bm{k}}-\mu}{T}\right)\bigg{missing})$
$\displaystyle~{}~{}~{}~{}~{}+\ln\bigg(1+\exp\left(-\frac{\omega_{\bm{k}}+\mu}{T}\right)\bigg{missing})\Bigg{]}.$
(13)

Using Eqs. (10) and (13) we can easily find

$\displaystyle\chi_{1}^{(B)}$
$\displaystyle=\frac{\partial(P/T^{4})}{\partial(\mu/T)}\bigg{|}_{T}$
$\displaystyle=\frac{2g}{T^{3}}\int
dK\Big{[}f(\omega_{\bm{k}})-\bar{f}(\omega_{\bm{k}})\Big{]}$
$\displaystyle=\frac{n_{B}}{T^{3}}$ (14)

and

$\displaystyle\chi_{2}^{(B)}$
$\displaystyle=\frac{\partial^{2}(P/T^{4})}{\partial(\mu/T)^{2}}\bigg{|}_{T}$
$\displaystyle=\frac{2g}{T^{3}}\int
dK\Big{[}f(\omega_{\bm{k}})(1-f(\omega_{\bm{k}}))$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\bar{f}(\omega_{\bm{k}})(1-\bar{f}(\omega_{\bm{k}}))\Big{]}.$
(15)

Therefore, we can write

$\displaystyle T^{3}\chi_{2}^{(B)}$ $\displaystyle=V\langle(n_{B}-\langle
n_{B}\rangle)^{2}\rangle,$ $\displaystyle=2g\int
dK\Big{[}f(\omega_{\bm{k}})(1-f(\omega_{\bm{k}}))$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\bar{f}(\omega_{\bm{k}})(1-\bar{f}(\omega_{\bm{k}}))\Big{]}.$
(16)

In the large volume limit Eq. (9) should be consistent with Eq. (16). This can
be verified using the Gaussian representation of the three-dimensional Dirac
delta function

$\displaystyle\delta^{(3)}({\bm{k}}-{\bm{k}^{\prime}})=\lim_{a\to\infty}\frac{a^{3}}{(2\pi)^{3/2}}e^{-\frac{a^{2}}{2}({\bm{k}}-{\bm{k}^{\prime}})^{2}}.$
(17)

This leads us to the following expression

$\displaystyle\lim_{a\to\infty}a^{3}(2\pi)^{3/2}\bigg{[}\langle:\hat{J}_{a}^{0}::\hat{J}_{a}^{0}:\rangle-\langle:\hat{J}_{a}^{0}:\rangle^{2}\bigg{]}$
$\displaystyle=\int\frac{d^{3}k}{(2\pi)^{3}}d^{3}k^{\prime}\frac{1}{\omega_{\bm{k}}}\frac{1}{\omega_{\bm{k}^{\prime}}}(\omega_{\bm{k}}\omega_{\bm{k}^{\prime}}+\bm{k}\cdot\bm{k}^{\prime}+m^{2})$
$\displaystyle\left[f(\omega_{\bm{k}})(1-f(\omega_{\bm{k}^{\prime}}))+\bar{f}(\omega_{\bm{k}})(1-\bar{f}(\omega_{\bm{k}^{\prime}}))\right]\delta^{(3)}(\bm{k}-\bm{k}^{\prime})$
$\displaystyle-\int\frac{d^{3}k}{(2\pi)^{3}}d^{3}k^{\prime}\frac{1}{\omega_{\bm{k}}}\frac{1}{\omega_{\bm{k}^{\prime}}}(\omega_{\bm{k}}\omega_{\bm{k}^{\prime}}+\bm{k}\cdot\bm{k}^{\prime}-m^{2})$
$\displaystyle\left[f(\omega_{\bm{k}})(1-\bar{f}(\omega_{\bm{k}^{\prime}}))+\bar{f}(\omega_{\bm{k}})(1-f(\omega_{\bm{k}^{\prime}}))\right]\delta^{(3)}(\bm{k}+\bm{k}^{\prime})$
$\displaystyle=2\int
dK\Big{[}f(\omega_{\bm{k}})(1-f(\omega_{\bm{k}}))+\bar{f}(\omega_{\bm{k}})(1-\bar{f}(\omega_{\bm{k}}))\Big{]},$
(18)

hence,

$\displaystyle\lim_{a\to\infty}V_{a}\bigg{[}\langle:\hat{J}_{a}^{0}::\hat{J}_{a}^{0}:\rangle-\langle:\hat{J}_{a}^{0}:\rangle^{2}\bigg{]}=T^{3}\chi_{2}^{(B)},$
(19)

where $V_{a}=a^{3}(2\pi)^{3/2}$ may be identified as the volume of the
“Gaussian” subsystem $S_{a}$.

Figure 4: Variation of normalized fluctuation
$V_{a}\sigma^{2}/(T^{3}\chi_{2}^{(B)})$ for different values of temperature
($T$) but with fixed baryon chemical potential ($\mu$) and particle mass
($m$). Figure 5: Variation of normalized fluctuation
$V_{a}\sigma^{2}/(T^{3}\chi_{2}^{(B)})$ for different values of baryon
chemical potential ($\mu$) but with fixed temperature ($T$) and particle mass
($m$). Figure 6: Variation of normalized fluctuation
$V_{a}\sigma^{2}/(T^{3}\chi_{2}^{(B)})$ for different values of particle mass
($m$) but with fixed temperature ($T$) and baryon chemical potential ($\mu$).

##  IV Numerical results

Our main result describing fluctuations of the baryon number density in a hot
and dense gas of baryons is represented by Eq. (9). By straightforward
numerical integration, we can obtain the results for any subsystem of size
$a$, temperature $T$, baryon chemical potential $\mu$, and particle mass $m$.
In Figs. (1), (2) and (3) we show the variation of the normalized fluctuation
$\sigma_{n}$ as defined in Eq. (6) with the subsystem size $a$ for different
values of temperature, baryon chemical potential and mass, respectively.
Keeping in mind the relativistic heavy ion collisions in the regime where
large baryon number densities are expected, we consider the following range of
the parameters: $0.1$ GeV $\leq T\leq 0.25$ GeV, $0.5$ GeV $\leq\mu\leq 1.0$
GeV and $0.0$ $\leq m\leq 1.0$ GeV. The internal degeneracy factor is taken to
be $g=10$.

Figs. (1), (2) and (3) show that the normalized fluctuation $\sigma_{n}$
decreases with the increase in the size of the subsystem $S_{a}$, which is the
expected behavior of fluctuations. Furthermore, for small system size, the
fluctuation is substantially large, which is the manifestation of the quantum
mechanical behavior. From these figures, it is clear that the normalized
fluctuation $\sigma_{n}$ decreases with an increase in temperature and baryon
chemical potential, however, $\sigma_{n}$ is large if the particle mass is
large. We should mention that variation of $\sigma$ with temperature, baryon
chemical potential and mass is opposite to the variation of $\sigma_{n}$ (not
shown here).

In Figs. (4), (5) and (6) we demonstrate the thermodynamic limit of the volume
scaled fluctuation $V_{a}\sigma^{2}$. From Eq. (19) we conclude that in the
thermodynamic limit, i.e., in the $a\rightarrow\infty$ limit,
$V_{a}\sigma^{2}/(T^{3}\chi_{2}^{(B)})$ should approach unity. This property
can be seen clearly in Figs. (4), (5) and (6), where for $a>1$ fm quantum
fluctuations approaches the classical limit which also depend upon the
temperature ($T$), baryon chemical potential ($\mu$) and particle mass ($m$).
Furthermore, for small system size the quantum fluctuations can be
significant.

##  V conclusion

In this work, we have analyzed quantum baryon-number fluctuations in
subsystems of a hot and dense relativistic gas of fermions and found that they
diverge for small system sizes. On the other hand, our results agree with the
results known from statistical physics for sufficiently large system size $a$.
In this way, we have delivered a useful formula that accounts for both
statistical and quantum features of the fluctuations. The numerical results
have been obtained for a broad range of thermodynamic parameters. They may be
useful to interpret and shed new light on the heavy-ion experimental data.

###### Acknowledgements.

This research was supported in part by the Polish National Science Centre
Grants No. 2016/23/B/ST2/00717 and No. 2018/30/E/ST2/00432.

## References

  * (1) M. Smoluchowski, “Beitrag zur Theorie der Opaleszenz von Gasen im kritischen Zustande,” Bulletin international de l’Académie des sciences de Cracovie (1911) 493–502.
  * (2) S. Jeon and V. Koch, “Charged particle ratio fluctuation as a signal for quark-gluon plasma,” Phys. Rev. Lett. 85 (Sep, 2000) 2076–2079.
  * (3) K. Huang,  Statistical Mechanics .  Wiley New York, 2nd ed., 1987.
  * (4) R. Kubo, “The Fluctuation-Dissipation Theorem,” Rep. Prog. Phys. 29 (1966) 255.
  * (5) E. Lifshitz and I. Khalatnikov, “Investigations in relativistic cosmology,” Adv. Phys. 12 (1963) 185–249.
  * (6) A. H. Guth and S.-Y. Pi, “Fluctuations in the new inflationary universe,” Phys. Rev. Lett. 49 (Oct, 1982) 1110–1113.
  * (7) A. Das, W. Florkowski, R. Ryblewski, and R. Singh, “Quantum fluctuations of energy in subsystems of a hot relativistic gas,” arXiv:2012.05662 [hep-ph].
  * (8) A. Das, W. Florkowski, R. Ryblewski, and R. Singh, “Pseudo-gauge dependence of quantum fluctuations of energy in a hot relativistic gas of fermions,” arXiv:2103.01013 [nucl-th].
  * (9) J. Berges and K. Rajagopal, “Color superconductivity and chiral symmetry restoration at nonzero baryon density and temperature,” Nucl. Phys. B 538 (1999) 215–232, arXiv:hep-ph/9804233.
  * (10) A. M. Halasz, A. D. Jackson, R. E. Shrock, M. A. Stephanov, and J. J. M. Verbaarschot, “On the phase diagram of QCD,” Phys. Rev. D 58 (1998) 096007, arXiv:hep-ph/9804290.
  * (11) M. A. Stephanov, K. Rajagopal, and E. V. Shuryak, “Signatures of the tricritical point in QCD,” Phys. Rev. Lett. 81 (1998) 4816–4819, arXiv:hep-ph/9806219.
  * (12) M. A. Stephanov, K. Rajagopal, and E. V. Shuryak, “Event-by-event fluctuations in heavy ion collisions and the QCD critical point,” Phys. Rev. D 60 (1999) 114028, arXiv:hep-ph/9903292.
  * (13) Y. Hatta and T. Ikeda, “Universality, the QCD critical / tricritical point and the quark number susceptibility,” Phys. Rev. D 67 (2003) 014028, arXiv:hep-ph/0210284.
  * (14) M. Stephanov, “Non-Gaussian fluctuations near the QCD critical point,” Phys. Rev. Lett. 102 (2009) 032301, arXiv:0809.3450 [hep-ph].
  * (15) B. Berdnikov and K. Rajagopal, “Slowing out-of-equilibrium near the QCD critical point,” Phys. Rev. D 61 (2000) 105017, arXiv:hep-ph/9912274.
  * (16) M. Kitazawa, M. Asakawa, and H. Ono, “Non-equilibrium time evolution of higher order cumulants of conserved charges and event-by-event analysis,” Phys. Lett. B 728 (2014) 386–392, arXiv:1307.2978 [nucl-th].
  * (17) L. Tinti and W. Florkowski, “Particle polarization, spin tensor and the Wigner distribution in relativistic systems,” arXiv:2007.04029 [nucl-th].
  * (18) C. Cohen-Tannoudji, B. Diu, F. Laloë, and S. R. Hemley, Quantum mechanics: Vol. 3: fermions, bosons, photons, correlations and entanglement.  A Wiley-Interscience publication. Wiley, New York, 1977.
  * (19) C. Itzykson and J. Zuber, Quantum Field Theory.  International Series In Pure and Applied Physics. McGraw-Hill, New York, 1980.
  * (20) T. Evans and D. A. Steer, “Wick’s theorem at finite temperature,” Nucl. Phys. B 474 (1996) 481–496, arXiv:hep-ph/9601268.
  * (21) S. Coleman, Lectures of Sidney Coleman on Quantum Field Theory.  WSP, Hackensack, 12, 2018.
  * (22) M. Nahrgang, M. Bluhm, P. Alba, R. Bellwied, and C. Ratti, “Impact of resonance regeneration and decay on the net-proton fluctuations in a hadron resonance gas,” Eur. Phys. J. C 75 no. 12, (2015) 573, arXiv:1402.1238 [hep-ph].


# On a perturbation analysis of Higham squared maximum Gaussian elimination
growth matrices

Alan Edelman  Department of Mathematics, Massachusetts Institute of
Technology, Cambridge, MA, 02139 USA<EMAIL_ADDRESS>,  John Urschel
<EMAIL_ADDRESS>and  Bowen Zhu<EMAIL_ADDRESS>School of Engineering
and Applied Sciences, Harvard University, Cambridge, MA, 02138 USA.

###### Abstract.

Gaussian elimination is the most popular technique for solving a dense linear
system. Large errors in this procedure can occur in floating point arithmetic
when the matrix’s growth factor is large. We study this potential issue and
how perturbations can improve the robustness of the Gaussian elimination
algorithm. In their 1989 paper, Higham and Higham characterized the complete
set of real n by n matrices that achieves the maximum growth factor under
partial pivoting. This set of matrices serves as the critical focus of this
work. Through theoretical insights and empirical results, we illustrate the
high sensitivity of the growth factor of these matrices to perturbations and
show how subtle changes can be strategically applied to matrix entries to
significantly reduce the growth, thus enhancing computational stability and
accuracy.

###### 2020 Mathematics Subject Classification:

Primary 65F05, 15A23.

##  1\. Introduction to Partial Pivoting Growth

It is well-known, even to beginning students in engineering and science, that
the most popular method for solving the dense linear system $Ax=b$ for $x$ is
Gaussian elimination. Indeed, this method, employed using partial pivoting, is
widely available through interfaces from high level languages such as Julia
[1], Mathematica [12], Matlab [11], Python NumPy [7], R [13], etc. Error
estimates for the stability of the Gaussian elimination algorithm in floating
point arithmetic are governed by the bits of precision used, the condition
number of the matrix, and the growth factor (i.e., the largest magnitude entry
encountered during the Gaussian elimination algorithm) [5, Theorem 3.3.2].
Many researchers have studied and continue to study the question of why
Gaussian elimination with partial pivoting has been so very effective [4, 8,
9, 10, 14, 15, 16, 17]. In contrast to complete pivoting, where the existence
of matrices with even super-linear growth remains an open problem [2, 3, 6],
it has been known since Wilkinson’s classic text The Algebraic Eigenvalue
Problem [18, p.212] that, for partial pivoting, the growth factor is bounded
above by $2^{n-1}$ and that this quantity can be achieved by the matrix

$A=\begin{pmatrix}[r]1&0\phantom{.}&\cdots&0&1\\\
-1&\ddots&\ddots&\vdots&\vdots\\\ \vdots&\ddots&1&0&1\\\ -1&\cdots&-1&1&1\\\
-1&\cdots&-1&-1&1\end{pmatrix}.$ (1.1)

Much later, Higham and Higham (from now on denoted Higham222We denote “Higham
and Higham” as Higham2, to be read as “Higham squared,” yet we realize this
has the appearance of a footnote, so for readers who saw it this way, we have
included this footnote. ) identified the complete set of $n$ by $n$ real
matrices that achieve the maximal growth of $2^{n-1}$ [9]. We call such
matrices Higham2 matrices (see Proposition 2.2 for a description). A scalar
quantity of interest is the last pivot of a Higham2 matrix, which is a
differentiable function of the matrix entries. We can therefore ask for the
gradient of this last pivot or, even better, to have a full (non-
infinitesimal) perturbation analysis of the last pivot (for Gaussian
elimination without pivoting). We provide such a perturbation analysis in
Theorem 2.3. The last pivot is an ideal quantity to measure in order to
understand the growth factor, as every entry of $U$ is the last pivot of the
LU factorization of some submatrix of $A$. We observe that generically, large
growth does not last very long in the sense that often a small perturbation
can dramatically reduce a large pivot. We have a mental image that the Higham2
matrices live on a kind of “ridge” that one can easily fall off of. This
picture is consistent with the smoothed analysis of Sankar, Spielman, and Teng
[14]. The structure of the Higham2 matrices provide an ideal setting to better
understand the ridge and its profile. Perhaps unsurprisingly, not all
directions of descent are created equal. We provide numerical experiments (in
Section 3) to visualize the effects of perturbing Higham2 matrices and confirm
the conclusions gleaned from the theoretical results of Section 2.

##  2\. Entrywise Perturbations & Higham2 Matrices

Here we provide mathematical estimates for the effects of entrywise
perturbations on the last pivot of the LU factorization (without pivoting) of
a matrix (Lemma 2.1), recall a characterization of Higham2 matrices
(Proposition 2.2), and consider the effects of entrywise perturbations on this
class (Theorem 2.3, Corollaries 2.4 & 2.5). These theoretical results give
insight into the experimental observations in Section 3.

######  Lemma 2.1.

Let

$A=\begin{pmatrix}\hat{L}&0\\\
\bm{\ell}^{T}&1\end{pmatrix}\begin{pmatrix}\hat{U}&\bm{u}\\\
0&p\end{pmatrix}\in\mathrm{GL}_{n}(\mathbb{R}),$

where $\hat{L}\in\mathrm{SL}_{n-1}(\mathbb{R})$ is lower unitriangular,
$\hat{U}\in\mathrm{GL}_{n-1}(\mathbb{R})$ is upper triangular, and
$\bm{\ell},\bm{u}\in\mathbb{R}^{n-1}$. Then the LU factorization, if it
exists, of $A+\epsilon\,\bm{e}_{i}\bm{e}_{j}^{T}$, where $\bm{e}_{i}$ is the
$i^{th}$ standard basis vector, has last pivot

$p^{(i,j)}_{\epsilon}=p+\epsilon\frac{(\hat{U}^{-1}\bm{u})_{j}(\bm{\ell}^{T}\hat{L}^{-1})_{i}}{1+\epsilon(\hat{L}\hat{U})^{-1}_{ji}}\qquad\text{for}\quad
i,j<n,$

$p^{(i,n)}_{\epsilon}=p-\epsilon(\bm{\ell}^{T}\hat{L}^{-1})_{i}$ for $i<n$,
$p^{(n,j)}_{\epsilon}=p-\epsilon(\hat{U}^{-1}\bm{u})_{j}$ for $j<n$, and
$p^{(n,n)}_{\epsilon}=p+\epsilon$.

###### Proof.

We first consider the case $i,j<n$. The last pivot is a ratio of determinants

$p^{(i,j)}_{\epsilon}=\frac{\det(A+\epsilon\,\bm{e}_{i}\bm{e}_{j}^{T})}{\det(\hat{L}\hat{U}+\epsilon\,\bm{\hat{e}}_{i}\bm{\hat{e}}_{j}^{T})}=\frac{\det(A)}{\det(\hat{L}\hat{U})}\,\bigg{(}\frac{1+\epsilon\bm{e}_{j}^{T}A^{-1}\bm{e}_{i}}{1+\epsilon\bm{\hat{e}}_{j}^{T}(\hat{L}\hat{U})^{-1}\bm{\hat{e}}_{i}}\bigg{)}=p\bigg{(}\frac{1+\epsilon
A^{-1}_{ji}}{1+\epsilon(\hat{L}\hat{U})^{-1}_{ji}}\bigg{)},$

where $\bm{\hat{e}}_{i}$ is the $i^{th}$ standard basis vector in
$\mathbb{R}^{n-1}$. The matrix $A^{-1}$ has block form

$A^{-1}=\begin{pmatrix}\hat{U}^{-1}&-p^{-1}\hat{U}^{-1}\bm{u}\\\
0&p^{-1}\end{pmatrix}\begin{pmatrix}\hat{L}^{-1}&0\\\
-\bm{\ell}^{T}\hat{L}^{-1}&1\end{pmatrix}=\begin{pmatrix}(\hat{L}\hat{U})^{-1}+p^{-1}\hat{U}^{-1}\bm{u}\bm{\ell}^{T}\hat{L}^{-1}&-p^{-1}\hat{U}^{-1}\bm{u}\\\
-p^{-1}\bm{\ell}^{T}\hat{L}^{-1}&p^{-1}\end{pmatrix},$

and so
$A^{-1}_{ji}=(\hat{L}\hat{U})_{ji}^{-1}+p^{-1}(\hat{U}^{-1}\bm{u})_{j}(\bm{\ell}^{T}\hat{L}^{-1})_{i}$.
Therefore,

$p^{(i,j)}_{\epsilon}-p=p\bigg{[}\frac{1+\epsilon
A^{-1}_{ji}}{1+\epsilon(\hat{L}\hat{U})^{-1}_{ji}}-1\bigg{]}=\frac{\epsilon
p(A^{-1}_{ji}-(\hat{L}\hat{U})^{-1}_{ji})}{1+\epsilon(\hat{L}\hat{U})^{-1}_{ji}}=\frac{\epsilon(\hat{U}^{-1}\bm{u})_{j}(\bm{\ell}^{T}\hat{L}^{-1})_{i}}{1+\epsilon(\hat{L}\hat{U})^{-1}_{ji}}.$

When $i=n$ or $j=n$,
$\displaystyle{p^{(i,j)}_{\epsilon}=\frac{\det(A+\epsilon\,\bm{e}_{i}\bm{e}_{j}^{T})}{\det(\hat{L}\hat{U})}=p(1+\epsilon
A^{-1}_{ji}).}$ Noting that $A_{jn}^{-1}=-p^{-1}(\hat{U}^{-1}\bm{u})_{j}$ for
$j<n$, $A_{ni}^{-1}=-p^{-1}(\bm{\ell}^{T}\hat{L}^{-1})_{i}$ for $i<n$, and
$A_{nn}^{-1}=p^{-1}$ completes the proof. ∎

We recall the following characterization of Higham2 matrices from the original
paper of Higham and Higham [9], where we have slightly adjusted the
normalization and notation to suit our needs.

######  Proposition 2.2.

[9, Theorem 2.2] Every matrix $A\in\mathrm{GL}_{n}(\mathbb{R})$,
$\|A\|_{\max}=1$, with growth factor under partial pivoting equal to $2^{n-1}$
is of the form

$DPA=\begin{pmatrix}\hat{L}&0\\\
-\bm{1}^{T}&1\end{pmatrix}\begin{pmatrix}\hat{U}&\bm{u}\\\
0&2^{n-1}\end{pmatrix}=\begin{pmatrix}\hat{L}\hat{U}&\bm{1}\\\
-\bm{1}^{T}\hat{U}&1\end{pmatrix},$ (2.1)

where $D\in\mathrm{GL}_{n}(\mathbb{R})$ is a $\pm 1$ diagonal matrix,
$P\in\mathrm{GL}_{n}(\mathbb{R})$ is a permutation matrix associated with a
partial pivoting of $A$, $\hat{L}\in\mathrm{GL}_{n-1}(\mathbb{R})$ is lower
unitriangular with $\hat{L}_{ij}=-1$ for all $i>j$, $\bm{1}$,
$\bm{u}=(1,2,...,2^{n-2})^{T}\in\mathbb{R}^{n-1}$, and
$\hat{U}\in\mathrm{GL}_{n-1}(\mathbb{R})$ is upper triangular, with entries
satisfying $\|\hat{L}\hat{U}\|_{\max}\leq 1$ and
$\|\bm{1}^{T}\hat{U}\|_{\infty}\leq 1$.

Applying Lemma 2.1 to Higham2 matrices (described in Proposition 2.2) gives
the following theorem.

######  Theorem 2.3.

Let $A\in\mathrm{GL}_{n}(\mathbb{R})$, $\|A\|_{\max}=1$, be a Higham2 matrix
of the form in Equation 2.1 with $P=D=I$. Then the LU factorization, if it
exists, of $A+\epsilon\,\bm{e}_{i}\bm{e}_{j}^{T}$ has last pivot

$p_{\epsilon}^{(i,j)}=\begin{cases}\displaystyle{\frac{1}{2^{1-n}+\epsilon\,2^{-i}\sum_{\ell=1}^{n-j}2^{-\ell}\hat{U}_{j,n-\ell}^{-1}}}&\qquad\text{for}\quad
i<j<n\\\\[12.91663pt]
\displaystyle{\frac{1+\tfrac{1}{2}\epsilon\,(\hat{U}_{ji}^{-1}-\sum_{\ell=1}^{i-j}2^{-\ell}\hat{U}_{j,i-\ell}^{-1})}{2^{1-n}+\epsilon(2^{1-n}\hat{U}_{ji}^{-1}+2^{-i}\sum_{\ell=1}^{n-i-1}2^{-\ell}\hat{U}_{j,n-\ell}^{-1})}}&\qquad\text{for}\quad
j\leq i<n\\\\[4.30554pt] \end{cases}\quad,$ (2.2)

$p_{\epsilon}^{(i,n)}=2^{n-1}(1+\epsilon\,2^{-i})$ for $i<n$,
$p_{\epsilon}^{(n,j)}=2^{n-1}(1-\epsilon\sum_{\ell=1}^{n-j}2^{-\ell}\hat{U}_{j,n-\ell}^{-1})$
for $j<n$, and $p_{\epsilon}^{(n,n)}=2^{n-1}+\epsilon$.

###### Proof.

By Proposition 2.2, $A=\begin{pmatrix}\hat{L}&0\\\
-\bm{1}^{T}&1\end{pmatrix}\begin{pmatrix}\hat{U}&\bm{u}\\\
0&2^{n-1}\end{pmatrix}$, where $\hat{L}\in\mathrm{SL}_{n-1}(\mathbb{R})$ is
lower unitriangular with $\hat{L}_{ij}=-1$ for all $i>j$,
$\bm{u}=(1,2,...,2^{n-2})^{T}$, and $\hat{U}\in\mathrm{GL}_{n-1}(\mathbb{R})$
is upper triangular. The matrix $\hat{L}$ has a simple structure, and its
inverse has entries given by $\hat{L}^{-1}_{ij}=\phi(i-j)$, where $\phi(k)$
equals zero for $k<0$, one for $k=0$, and $2^{k-1}$ for $k>0$. Therefore,

$(\bm{1}^{T}\hat{L}^{-1})_{i}=\sum_{k=i}^{n-1}\phi(k-i)=2^{n-i-1}\qquad\text{and}\qquad(\hat{L}\hat{U})_{ji}^{-1}=\sum_{k=i}^{n-1}\phi(k-i)\hat{U}_{jk}^{-1}=\hat{U}_{ji}^{-1}+\sum_{k=i+1}^{n-1}2^{k-i-1}\hat{U}_{jk}^{-1}.$

Applying Lemma 2.1 to $A$ and noting that
$(\hat{U}^{-1}\bm{u})_{j}=\sum_{k=j}^{n-1}2^{k-1}\hat{U}^{-1}_{jk}$, we have
$p_{\epsilon}^{(n,n)}=2^{n-1}+\epsilon$,
$p_{\epsilon}^{(i,n)}=2^{n-1}(1+\epsilon\,2^{-i})$ for $i<n$,
$p_{\epsilon}^{(n,j)}=2^{n-1}(1-\epsilon\sum_{\ell=1}^{n-j}2^{-\ell}\hat{U}_{j,n-\ell}^{-1})$
for $j<n$, and

$p_{\epsilon}^{(i,j)}=2^{n-1}-\frac{\epsilon
2^{n-i-1}\sum_{k=j}^{n-1}2^{k-1}\hat{U}^{-1}_{jk}}{1+\epsilon(\hat{U}_{ji}^{-1}+\sum_{k=i+1}^{n-1}2^{k-i-1}\hat{U}_{jk}^{-1})}\qquad\text{for}\quad
i,j<n.$

When $i<j<n$, $\hat{U}_{ji}^{-1}=0$ and
$\sum_{k=i+1}^{n-1}2^{k-i-1}\hat{U}_{jk}^{-1}=\sum_{k=j}^{n-1}2^{k-i-1}\hat{U}^{-1}_{jk}$,
and so $p_{\epsilon}^{(i,j)}$ equals
$(2^{1-n}+\epsilon\,2^{-i}\sum_{\ell=1}^{n-j}2^{-\ell}\hat{U}_{j,n-\ell}^{-1})^{-1}$.
Finally, in the case $j\leq i<n$, we have

$p_{\epsilon}^{(i,j)}=\frac{1+\tfrac{1}{2}\epsilon\,(\hat{U}_{ji}^{-1}-\sum_{\ell=1}^{i-j}2^{-\ell}\hat{U}_{j,i-\ell}^{-1})}{2^{1-n}+\epsilon(2^{1-n}\hat{U}_{ji}^{-1}+2^{-i}\sum_{\ell=1}^{n-i-1}2^{-\ell}\hat{U}_{j,n-\ell}^{-1})}.$

∎

Equation 2.2 of Theorem 2.3 deserves a number of observations. First, we note
a connection between the rate at which $p_{\epsilon}^{(i,j)}$ decreases and
the condition number of the matrix: if the last pivot under $\epsilon$
perturbation is much smaller than $1/|\epsilon|$, then $A$ is ill-conditioned.

######  Corollary 2.4.

Let $A\in\mathrm{GL}_{n}(\mathbb{R})$, $\|A\|_{\max}=1$, be a Higham2 matrix
of the form in Equation 2.1 with $P=D=I$. If $|\epsilon|<1$ and
$\sqrt{n}<|p_{\epsilon}^{(i,j)}|<2^{n-5}$, then
$\kappa_{2}(A)=\|A\|_{2}\|A^{-1}\|_{2}\geq\sqrt{n}\,\big{|}3\epsilon
p_{\epsilon}^{(i,j)}\big{|}^{-1}$.

###### Proof.

Let $|\epsilon|<1$ and $\sqrt{n}<|p_{\epsilon}^{(i,j)}|<2^{n-5}$ for some
$i,j$. We first relate $\kappa_{2}(A)$ to the entries of $U^{-1}$:

$\|A\|_{2}\|A^{-1}\|_{2}\geq\|A\mathbf{e}_{n}\|_{2}\|A^{-1}\|_{2}=\sqrt{n}\max_{y\neq
0}\frac{\|U^{-1}y\|_{2}}{\|Ly\|_{2}}\geq\sqrt{n}\max_{\ell}\frac{\|U^{-1}e_{n-\ell}\|_{2}}{\|Le_{n-\ell}\|_{2}}\geq\max_{k,\ell}\frac{\sqrt{n}\,|U^{-1}_{k,n-\ell}|}{\sqrt{\ell+1}}.$

What remains is to show that
$|U^{-1}_{k,n-\ell}|\geq\sqrt{\ell+1}\,\big{|}3\epsilon
p_{\epsilon}^{(i,j)}\big{|}^{-1}$ for some $k,\ell$. We proceed by
contradiction. If this is not the case, then

$\left|\epsilon\sum_{\ell=1}^{k}\frac{\hat{U}_{j,n-\ell}^{-1}}{2^{\ell}}\right|<\frac{1}{\big{|}p_{\epsilon}^{(i,j)}\big{|}}\sum_{\ell=1}^{\infty}\frac{\sqrt{\ell+1}}{3\times
2^{\ell}}<\frac{0.57}{\big{|}p_{\epsilon}^{(i,j)}\big{|}}\quad\text{and}\quad\left|\epsilon\,(\hat{U}_{ji}^{-1}-\sum_{\ell=1}^{i-j}2^{-\ell}\hat{U}_{j,i-\ell}^{-1})\right|<\sum_{\ell=0}^{\infty}\frac{1}{3\times
2^{\ell}}=\frac{2}{3}.$

These two bounds, combined with the formulae of Theorem 2.3, lead to a
contradiction for all choices of $(i,j)$. For example, if $j\leq i<n$, then

$\big{|}p_{\epsilon}^{(i,j)}\big{|}\geq\frac{1-\tfrac{1}{2}\left|\epsilon\,(\hat{U}_{ji}^{-1}-\sum_{\ell=1}^{i-j}2^{-\ell}\hat{U}_{j,i-\ell}^{-1})\right|}{2^{1-n}+\left|\epsilon(2^{1-n}\hat{U}_{ji}^{-1}+2^{-i}\sum_{\ell=1}^{n-i-1}2^{-\ell}\hat{U}_{j,n-\ell}^{-1})\right|}>\frac{2/3}{2^{1-n}+0.57/\big{|}p_{\epsilon}^{(i,j)}\big{|}}>\frac{2/3}{2^{-4}+0.57}\,\big{|}p_{\epsilon}^{(i,j)}\big{|},$

a contradiction. The remaining cases are similar, and are left to the reader.
∎

Now, let us restrict our attention to the case $i<j<n$, with $i$ relatively
small. As long as
$\epsilon\,2^{-i}\sum_{\ell=1}^{n-j}2^{-\ell}\hat{U}_{j,n-\ell}^{-1}$ is not
exponentially small in $n$, $p_{\epsilon}^{(i,j)}$ is approximately
$2^{i}\,(\epsilon\,\sum_{\ell=1}^{n-j}2^{-\ell}\hat{U}_{j,n-\ell}^{-1})^{-1}$.
One would expect this to be the case for “most” Higham2 matrices when
$\epsilon$ is only polynomially small, though exceptions certainly exist
(e.g., the Wilkinson matrix $\hat{U}=I$). This intuition is supported by
experimental results in Section 3. The case $i=1$ and $j=n-1$ is particularly
striking, and gives guaranteed improvement for all Higham2 matrices.

######  Corollary 2.5.

Let $A\in\mathrm{GL}_{n}(\mathbb{R})$, $\|A\|_{\max}=1$, be a Higham2 matrix
of the form in Equation 2.1 with $P=D=I$. Then

$|p^{(1,n-1)}_{\epsilon}|\leq\frac{4+o_{n}(1)}{|\epsilon|},\qquad\text{where}\qquad|o_{n}(1)|<2^{-(n-6)}|\epsilon|^{-1}\;\text{for}\;|\epsilon|>2^{-(n-4)}.$
(2.3)

###### Proof.

By Equation 2.2,

$p^{(1,n-1)}_{\epsilon}=\frac{4\hat{U}_{n-1,n-1}}{\epsilon+2^{-(n-3)}\hat{U}_{n-1,n-1}}=\frac{\hat{U}_{n-1,n-1}}{\epsilon}\bigg{(}4-\frac{2^{-(n-5)}\hat{U}_{n-1,n-1}}{\epsilon+2^{-(n-3)}\hat{U}_{n-1,n-1}}\bigg{)}.$

Noting that
$|\hat{U}_{n-1,n-1}|=\big{|}\big{[}(\hat{L}\hat{U})_{n-1,n-1}+(\bm{1}^{T}\hat{U})_{n-1}\big{]}/2\big{|}\leq
1$ gives the desired result. ∎

The entry $(1,n-1)$ is an example of a perturbation direction that always
produces a small last pivot when $\epsilon$ is only polynomially small.
Finally, we note that, while Lemma 2.1 and Theorem 2.3 apply only to the last
pivot, this general framework holds for arbitrary entries of $U$, as every
entry of $U$ is the last pivot of the LU factorization of some submatrix of
$A$. For instance, Corollary 2.5 implies that a only polynomially small
$\epsilon$ perturbation to the $(1,n-k)^{th}$ entry of a Higham2 matrix leads
to $|U_{n-k+1,n}|\leq(4+o_{n}(1))/|\epsilon|$ for any fixed $k$ (with $n$
growing). In Section 3, we make use of the insights gained from Theorem 2.3 to
suggest perturbations tailored to the most influential components of $A$, and
compare their effect to perturbations applied uniformly to $A$.

##  3\. Experimental Results

Here we perform two experiments, illustrated in Figures 1 and LABEL:fig:data.
First, in Figure 1, we consider the effects of an $\epsilon=10^{-8}$
perturbation to the last pivot of Higham2 matrices of dimension $n=100$.
Though $n=100$ is quite a small test-case, it is already more than sufficient
for our purposes, as, in double precision, the last pivot is nearly equal to
the inverse of machine epsilon squared. The heatmap on the left is of a
Higham2 matrix generated by taking the
$\hat{U}\in\mathrm{GL}_{n-1}(\mathbb{R})$ of Proposition 2.2 to have
independent standard normal entries, scaled so that $\|A\|_{\max}=1$, the map
on the right is of the Wilkinson matrix (see Equation 1.1), and the map in the
middle is of a matrix in between the two (w.r.t. choice of $\hat{U}$).
Perturbations in the top left portion of a random Higham2 matrix appear to be
most impactful. This is consistent with Theorem 2.3 and the fact that the
inverse of an upper triangular matrix with normal entries tends to have
exponentially large entries (in $n$) near the upper-right corner. However, the
Wilkinson matrix provides a clear reminder that this is not always the case.
The inverse of $\hat{U}=I$ has zero entries above the diagonal, rendering
perturbations to the top-left entries of $A$ relatively useless. Our results
are also consistent with Inequality 2.3: perturbing the top-right $(1,n-1)$
entry with a sufficiently large $\epsilon$ is a reliable way to always
decrease the last pivot size.

Figure 1. The effects of an $\epsilon=10^{-8}$ perturbation at different
entries of three $100\times 100$ Higham2 matrices for different choices of
$\hat{U}$ ($\hat{U}=I$ is the Wilkinson matrix). We plot a heatmap for each,
where the $(i,j)^{th}$ grid point is the value of
$\log_{10}|p^{(i,j)}_{\epsilon}|$ (in exact arithmetic).

Of course, we are not merely interested in the last pivot, but in the quality
of solution to $Ax=b$ we obtain using Gaussian elimination. Our theoretical
and experimental results in Theorem 2.3, Corollary 2.5, and Figure 1 give key
insights into the stability of the last pivot, with implications for the
growth factor itself, as every entry of $U$ is the last pivot of a sub-matrix
of $A$. In Figure LABEL:fig:data, we examine the effects of matrix
perturbation on the numerical solution to $Ax=b$ for Higham2 matrices using
Gaussian elimination with no pivoting. In the left plot, we observe that the
ill-conditioning of a random Higham2 matrix is a major barrier to a reasonable
solution. This is consistent with the theoretical observation that extremely
fast decay in the last pivot implies ill-conditioning (Corollary 2.4). In both
the middle and right plots, we observe that the perturbation to the first row
is the superior strategy. It is quite possible that this observation holds
more broadly than the class of maximum growth factor Higham2 matrices
considered here, and that there may be benefits to considering perturbations
tailored specifically to Gaussian elimination with partial pivoting.

## Acknowledgements

This material is based upon work supported by the National Science Foundation
under grant no. OAC-1835443, grant no. SII-2029670, grant no. ECCS-2029670,
grant no. OAC-2103804, and grant no. PHY-2021825. We also gratefully
acknowledge the U.S. Agency for International Development through Penn State
for grant no. S002283-USAID. The information, data, or work presented herein
was funded in part by the Advanced Research Projects Agency-Energy (ARPA-E),
U.S. Department of Energy, under Award Number DE-AR0001211 and DE-AR0001222.
The views and opinions of authors expressed herein do not necessarily state or
reflect those of the United States Government or any agency thereof. This
material was supported by The Research Council of Norway and Equinor ASA
through Research Council project “308817 - Digital wells for optimal
production and drainage”. Research was sponsored by the United States Air
Force Research Laboratory and the United States Air Force Artificial
Intelligence Accelerator and was accomplished under Cooperative Agreement
Number FA8750-19-2-1000. The views and conclusions contained in this document
are those of the authors and should not be interpreted as representing the
official policies, either expressed or implied, of the United States Air Force
or the U.S. Government. The U.S. Government is authorized to reproduce and
distribute reprints for Government purposes notwithstanding any copyright
notation herein. The third author’s masters thesis [19] contains some of the
results presented here, as well as some additional analysis and figures.

## References

  * [1] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah.  Julia: A fresh approach to numerical computing.  SIAM review, 59(1):65–98, 2017.
  * [2] Ankit Bisain, Alan Edelman, and John Urschel.  A new upper bound for the growth factor in Gaussian elimination with complete pivoting.  arXiv preprint arXiv:2312.00994, 2023.
  * [3] Alan Edelman and John Urschel.  Some new results on the maximum growth factor in Gaussian elimination.  SIAM Journal on Matrix Analysis and Applications, 45(2):967–991, 2024.
  * [4] Leslie V Foster.  Gaussian elimination with partial pivoting can fail in practice.  SIAM Journal on Matrix Analysis and Applications, 15(4):1354–1362, 1994.
  * [5] Gene H Golub and Charles F Van Loan.  Matrix computations.  JHU press, 2013.
  * [6] Nick Gould.  On growth in Gaussian elimination with complete pivoting.  SIAM Journal on Matrix Analysis and Applications, 12(2):354–361, 1991.
  * [7] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  Nature, 585(7825):357–362, September 2020.
  * [8] Desmond J Higham, Nicholas J Higham, and Srikara Pranesh.  Random matrices generating large growth in lu factorization with pivoting.  SIAM Journal on Matrix Analysis and Applications, 42(1):185–201, 2021.
  * [9] Nicholas J Higham and Desmond J Higham.  Large growth factors in Gaussian elimination with pivoting.  SIAM Journal on Matrix Analysis and Applications, 10(2):155–164, 1989.
  * [10] Han Huang and Konstantin Tikhomirov.  Average-case analysis of the Gaussian elimination with partial pivoting.  arXiv preprint arXiv:2206.01726, 2022.
  * [11] The MathWorks Inc.  Matlab version: 9.13.0 (r2022b), 2022.
  * [12] Wolfram Research, Inc.  Mathematica, Version 14.0.  Champaign, IL, 2024.
  * [13] R Core Team.  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria, 2021.
  * [14] Arvind Sankar, Daniel A Spielman, and Shang-Hua Teng.  Smoothed analysis of the condition numbers and growth factors of matrices.  SIAM Journal on Matrix Analysis and Applications, 28(2):446–476, 2006.
  * [15] Daniel A Spielman and Shang-Hua Teng.  Smoothed analysis of algorithms: Why the simplex algorithm usually takes polynomial time.  Journal of the ACM (JACM), 51(3):385–463, 2004.
  * [16] Lloyd N Trefethen and David Bau.  Numerical linear algebra, volume 181.  Siam, 2022.
  * [17] Lloyd N Trefethen and Robert S Schreiber.  Average-case stability of Gaussian elimination.  SIAM Journal on Matrix Analysis and Applications, 11(3):335–360, 1990.
  * [18] J.H. Wilkinson.  The Algebraic Eigenvalue Problem.  Clarendon Press, Oxford, 1965.
  * [19] Bowen Zhu.  Perturbation Analysis of Higham Squared Maximum Growth Matrices.  Master’s thesis, Harvard University, Cambridge, MA, May 2024.


$\displaystyle\mathrm{Var}_{\mathrm{post}}\left[f(x)\right]-\lim_{N\to\infty}\mathrm{Var}_{\mathrm{post}}\left[f(x)\right]$
$\displaystyle=\frac{c}{N}\nu\Sigma_{\perp}\propto\frac{1}{N}\propto\frac{1}{P}\propto\frac{1}{L}.$

##  11 Proof of Theorem 4.2

To prove Theorem 4.2 we begin with the bias-variance decomposition:

$\left\langle
f(x)-V_{0}x\right\rangle^{2}=\left(\theta_{*}x_{||}-V_{0}x\right)^{2}+\frac{\left|\left|x_{\perp}\right|\right|^{2}\left|\left|\theta_{*}\right|\right|^{2}}{P},$

where $x_{\perp}=x-x_{||}$ and

$x_{||}=\mathrm{im}(X)\mathrm{im}(X)^{T}x.$

Consider first the case when $\alpha_{0}<1$. Then

$\theta_{*}=V_{0}+\epsilon(X^{T}X)^{-1}X^{T}.$

Thus, the error introduced by the bias is

$\displaystyle\mathbb{E}\left[\left(\theta_{*}x_{||}-V_{0}x\right)^{2}\right]$
$\displaystyle=\mathbb{E}\left[(V_{0}x_{\perp})^{2}\right]+\mathbb{E}\left[(\epsilon(X^{T}X)^{-1}X^{T}x_{||})^{2}\right]$
$\displaystyle=\mathbb{E}\left[(V_{0}x)^{2}-2(V_{0}x)(V_{0}x_{||})+(V_{0}x_{||})^{2}\right]+\sigma_{\epsilon}^{2}\mathbb{E}\left[\|(X^{T}X)^{-1}X^{T}x_{||}\|^{2}\right]$
$\displaystyle=1-\alpha_{0}+\sigma_{\epsilon}^{2}\mathbb{E}\left[\|(X^{T}X)^{-1}X^{T}x_{||}\|^{2}\right].$

We observe that, for $X^{\dagger}=(X^{T}X)^{-1}X^{T}$,

$\displaystyle\mathbb{E}\left[\|(X^{T}X)^{-1}X^{T}x_{||}\|^{2}\right]$
$\displaystyle=\mathbb{E}\left[\mathrm{Tr}\left(x_{||}^{T}\left(X^{\dagger}\right)^{T}X^{\dagger}x_{||}\right)\right]$
$\displaystyle=\mathbb{E}\left[\mathrm{Tr}\left(\left(X^{\dagger}\right)^{T}X^{\dagger}x_{||}x_{||}^{T}\right)\right]$
$\displaystyle=\frac{\alpha_{0}}{1-\alpha_{0}}.$

Hence, the bias is

$\displaystyle\mathbb{E}\left[\left(\theta_{*}x_{||}-V_{0}x\right)^{2}\right]$
$\displaystyle=1-\alpha_{0}+\frac{\alpha_{0}}{1-\alpha_{0}}\sigma_{\epsilon}^{2}.$
(11.1)

The error introduced by the variance is

$\displaystyle\mathbb{E}\left[\frac{\left|\left|x_{\perp}\right|\right|^{2}\left|\left|\theta_{*}\right|\right|^{2}}{P}\right]$
$\displaystyle=\frac{1}{P}\mathbb{E}\left[\left(\|V_{0}\|^{2}+\|\epsilon(X^{T}X)^{-1}X^{T}\|^{2}\right)\left(\|x\|^{2}+\|x_{||}\|^{2}-2\|\mathrm{im}(X)^{T}x\|^{2}\right)\right]$
$\displaystyle=\frac{1-\alpha_{0}}{P}\mathbb{E}\left[\|x\|^{2}\left(1+\sigma_{\epsilon}^{2}\|(X^{T}X)^{-1}X^{T}\|^{2}\right)\right]$
$\displaystyle=\left(\frac{1}{\alpha_{0}}-1\right)\left(1+\sigma_{\epsilon}^{2}\mathbb{E}\left[\|(X^{T}X)^{-1}X^{T}\|^{2}\right]\right).$

Observing that

$\displaystyle\mathbb{E}\left[\left\|(X^{T}X)^{-1}X^{T}\right\|^{2}\right]$
$\displaystyle=\mathbb{E}\left[\mathrm{Tr}\left((X^{T}X)^{-1}X^{T}X(X^{T}X)^{-1}\right)\right]=\mathbb{E}\left[\mathrm{Tr}\left((X^{T}X)^{-1}\right)\right]=\frac{\alpha_{0}}{1-\alpha_{0}}$

from the $-1$st moment of the Marchenko-Pastur distribution for
$\alpha_{0}<1$, we obtain

$\displaystyle\left\langle f(x)-V_{0}x\right\rangle^{2}$
$\displaystyle=1-\alpha_{0}+\frac{\alpha_{0}}{1-\alpha_{0}}\sigma_{\epsilon}^{2}+\left(\frac{1}{\alpha_{0}}-1\right)\left(1+\frac{\alpha_{0}}{1-\alpha_{0}}\sigma_{\epsilon}^{2}\right)$
$\displaystyle=\frac{1}{\alpha_{0}}-\alpha_{0}+\frac{\sigma_{\epsilon}^{2}}{1-\alpha_{0}}.$

In the case of $\alpha_{0}>1$, the variance is zero since
$\|x_{\perp}\|^{2}=0$. Given

$\displaystyle\theta_{*}$ $\displaystyle=V_{0}+\epsilon X^{T}(XX^{T})^{-1},$

the total error originates from the bias, i.e.,

$\displaystyle\left\langle f(x)-V_{0}x\right\rangle^{2}$
$\displaystyle=\mathbb{E}\left[(\epsilon(X^{T}X)^{-1}X^{T}x)^{2}\right]$
$\displaystyle=\frac{\sigma_{\epsilon}^{2}}{\alpha_{0}-1}$ (11.2)

similarly to the bias computation for $\alpha_{0}<1$.

## References

  * [ACGH19] Sanjeev Arora, Nadav Cohen, Noah Golowich, and Wei Hu.  A convergence analysis of gradient descent for deep linear neural networks.  ICLR, 2019.
  * [ACH18] Sanjeev Arora, Nadav Cohen, and Elad Hazan.  On the optimization of deep networks: Implicit acceleration by overparameterization.  ICML, 2018.
  * [ALP22] Ben Adlam, Jake Levinson, and Jeffrey Pennington.  A random matrix perspective on mixtures of nonlinearities for deep learning.  AISTATS, 2022.
  * [AP20] Ben Adlam and Jeffrey Pennington.  The neural tangent kernel in high dimensions: Triple descent and a multi-scale theory of generalization.  In International Conference on Machine Learning, pages 74–84. PMLR, 2020.
  * [APH+21] Moloud Abdar, Farhad Pourpanah, Sadiq Hussain, Dana Rezazadegan, Li Liu, Mohammad Ghavamzadeh, Paul Fieguth, Xiaochun Cao, Abbas Khosravi, U Rajendra Acharya, et al.  A review of uncertainty quantification in deep learning: Techniques, applications and challenges.  Information Fusion, 76:243–297, 2021.
  * [AS64] Milton Abramowitz and Irene A Stegun.  Handbook of mathematical functions with formulas, graphs, and mathematical tables, volume 55.  US Government printing office, 1964.
  * [ASS20] Madhu S Advani, Andrew M Saxe, and Haim Sompolinsky.  High-dimensional dynamics of generalization error in neural networks.  Neural Networks, 132:428–446, 2020.
  * [BDK+21] Yasaman Bahri, Ethan Dyer, Jared Kaplan, Jaehoon Lee, and Utkarsh Sharma.  Explaining neural scaling laws.  arXiv preprint arXiv:2102.06701, 2021.
  * [BH89] Pierre Baldi and Kurt Hornik.  Neural networks and principal component analysis: Learning from examples without local minima.  Neural networks, 2(1):53–58, 1989.
  * [BHMM19] Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal.  Reconciling modern machine-learning practice and the classical bias–variance trade-off.  Proceedings of the National Academy of Sciences, 116(32):15849–15854, 2019.
  * [BHX20] Mikhail Belkin, Daniel Hsu, and Ji Xu.  Two models of double descent for weak features.  SIAM Journal on Mathematics of Data Science, 2(4):1167–1180, 2020\.
  * [BLLT20] Peter L Bartlett, Philip M Long, Gábor Lugosi, and Alexander Tsigler.  Benign overfitting in linear regression.  Proceedings of the National Academy of Sciences, 2020.
  * [Bor14] Emile Borel.  Introduction géométrique à quelques théories physiques.  Gauthier-Villars, 1914.
  * [CB18] Lenaic Chizat and Francis Bach.  On the global convergence of gradient descent for over-parameterized models using optimal transport.  In Advances in neural information processing systems, pages 3036–3046, 2018.
  * [DF87] Persi Diaconis and David Freedman.  A dozen de finetti-style results in search of a theory.  In Annales de l’IHP Probabilités et statistiques, volume 23, pages 397–423, 1987.
  * [DGST18] Matthew M Dunlop, Mark A Girolami, Andrew M Stuart, and Aretha L Teckentrup.  How deep are deep gaussian processes?  Journal of Machine Learning Research, 19(54):1–46, 2018.
  * [DH19] Simon Du and Wei Hu.  Width provably matters in optimization for deep linear neural networks.  In International Conference on Machine Learning, pages 1655–1664. PMLR, 2019.
  * [DL13] Andreas Damianou and Neil D Lawrence.  Deep gaussian processes.  In Artificial intelligence and statistics, pages 207–215. PMLR, 2013.
  * [DLT+18] Simon S Du, Jason D Lee, Yuandong Tian, Barnabas Poczos, and Aarti Singh.  Gradient descent learns one-hidden-layer cnn: Don’t be afraid of spurious local minima.  ICML, 2018.
  * [FSH21] Gianluca Finocchio and Johannes Schmidt-Hieber.  Posterior contraction for deep gaussian process priors.  arXiv preprint arXiv:2105.07410, 2021.
  * [GG16] Yarin Gal and Zoubin Ghahramani.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  In international conference on machine learning, pages 1050–1059. PMLR, 2016.
  * [GRM+18] Marta Garnelo, Dan Rosenbaum, Christopher Maddison, Tiago Ramalho, David Saxton, Murray Shanahan, Yee Whye Teh, Danilo Rezende, and SM Ali Eslami.  Conditional neural processes.  In International Conference on Machine Learning, pages 1704–1713. PMLR, 2018.
  * [Han18] Boris Hanin.  Which neural net architectures give rise to exploding and vanishing gradients?  In Advances in Neural Information Processing Systems, 2018.
  * [Han21] Boris Hanin.  Correlation functions in random fully connected neural networks at finite width.  arXiv preprint arXiv:2204.01058, 2021.
  * [HBM+22] Jordan Hoffmann, Sebastian Borgeaud, Arthur Mensch, Elena Buchatskaya, Trevor Cai, Eliza Rutherford, Diego de Las Casas, Lisa Anne Hendricks, Johannes Welbl, Aidan Clark, et al.  Training compute-optimal large language models.  arXiv preprint arXiv:2203.15556, 2022.
  * [HHLMF18] Marton Havasi, José Miguel Hernández-Lobato, and Juan José Murillo-Fuentes.  Inference in deep gaussian processes using stochastic gradient hamiltonian monte carlo.  Advances in neural information processing systems, 31, 2018.
  * [HLT20] Bobby He, Balaji Lakshminarayanan, and Yee Whye Teh.  Bayesian deep ensembles via the neural tangent kernel.  Advances in neural information processing systems, 33:1010–1022, 2020.
  * [HMRT22] Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan J Tibshirani.  Surprises in high-dimensional ridgeless least squares interpolation.  The Annals of Statistics, 50(2):949–986, 2022.
  * [HN19] Boris Hanin and Mihai Nica.  Finite depth and width corrections to the neural tangent kernel.  ICLR 2020 and arXiv:1909.05989, 2019.
  * [HN20] Boris Hanin and Mihai Nica.  Products of many large random matrices and gradients in deep neural networks.  Communications in Mathematical Physics, 376(1):287–322, 2020.
  * [HNPSD22] Jiri Hron, Roman Novak, Jeffrey Pennington, and Jascha Sohl-Dickstein.  Wide bayesian neural networks have a simple weight posterior: theory and accelerated sampling.  In International Conference on Machine Learning, pages 8926–8945. PMLR, 2022.
  * [HP21] Boris Hanin and Grigoris Paouris.  Non-asymptotic results for singular values of gaussian matrix products.  Geometric and Functional Analysis, 31(2):268–324, 2021.
  * [HR18] Boris Hanin and David Rolnick.  How to start training: The effect of initialization and architecture.  In Advances in Neural Information Processing Systems, pages 571–581, 2018.
  * [IVHW21] Pavel Izmailov, Sharad Vikram, Matthew D Hoffman, and Andrew Gordon Gordon Wilson.  What are bayesian neural network posteriors really like?  In International conference on machine learning, pages 4629–4640. PMLR, 2021.
  * [JGH18] Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  In Advances in neural information processing systems, pages 8571–8580, 2018.
  * [JT19] Ziwei Ji and Matus Telgarsky.  Gradient descent aligns the layers of deep linear networks.  ICLR, 2019.
  * [Kaw16] Kenji Kawaguchi.  Deep learning without poor local minima.  In Advances in Neural Information Processing Systems, pages 586–594, 2016.
  * [KMH+20] 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.  arXiv preprint arXiv:2001.08361, 2020.
  * [LB18] Thomas Laurent and James Brecht.  Deep linear networks with arbitrary loss: All local minima are global.  In International conference on machine learning, pages 2902–2907. PMLR, 2018.
  * [LBN+18] Jaehoon Lee, Yasaman Bahri, Roman Novak, Samuel S Schoenholz, Jeffrey Pennington, and Jascha Sohl-Dickstein.  Deep neural networks as gaussian processes.  ICML 2018 andarXiv:1711.00165, 2018.
  * [LG18] Andrew K Lampinen and Surya Ganguli.  An analytic theory of generalization dynamics and transfer learning in deep linear networks.  arXiv preprint arXiv:1809.10374, 2018.
  * [LIB+22] Sanae Lotfi, Pavel Izmailov, Gregory Benton, Micah Goldblum, and Andrew Gordon Wilson.  Bayesian model selection, the marginal likelihood, and generalization.  arXiv preprint arXiv:2202.11678, 2022.
  * [LK17] Haihao Lu and Kenji Kawaguchi.  Depth creates no bad local minima.  arXiv preprint arXiv:1702.08580, 2017.
  * [LNR22] Mufan Bill Li, Mihai Nica, and Daniel M Roy.  The neural covariance sde: Shaped infinite depth-and-width networks at initialization.  arXiv preprint arXiv:2206.02768, 2022.
  * [LS21] Qianyi Li and Haim Sompolinsky.  Statistical mechanics of deep linear neural networks: The backpropagating kernel renormalization.  Physical Review X, 11(3):031059, 2021.
  * [Mac92] David JC MacKay.  Bayesian interpolation.  Neural computation, 4(3):415–447, 1992.
  * [Mei36] CS Meijer.  Über whittakersche bzw. besselsche funktionen und deren produkte.  Nieuw Archief voor Wiskunde, 18(2):10–29, 1936.
  * [MM19] Song Mei and Andrea Montanari.  The generalization error of random features regression: Precise asymptotics and the double descent curve.  Communications on Pure and Applied Mathematics, 2019.
  * [MMM21] Song Mei, Theodor Misiakiewicz, and Andrea Montanari.  Generalization error of random feature and kernel methods: hypercontractivity and kernel matrix concentration.  Applied and Computational Harmonic Analysis, 2021.
  * [MMN18] Song Mei, Andrea Montanari, and Phan-Minh Nguyen.  A mean field view of the landscape of two-layer neural networks.  Proceedings of the National Academy of Sciences, 115(33):E7665–E7671, 2018.
  * [MRS22] Alexander Maloney, Daniel A Roberts, and James Sully.  A solvable model of neural scaling laws.  arXiv preprint arXiv:2210.16859, 2022.
  * [MZ22] Andrea Montanari and Yiqiao Zhong.  The interpolation phase transition in neural networks: Memorization and generalization under lazy training.  The Annals of Statistics, 50(5):2816–2847, 2022.
  * [NBR+21] Lorenzo Noci, Gregor Bachmann, Kevin Roth, Sebastian Nowozin, and Thomas Hofmann.  Precise characterization of the prior predictive distribution of deep relu networks.  Advances in Neural Information Processing Systems, 34:20851–20862, 2021.
  * [NDSR21] Gadi Naveh, Oded Ben David, Haim Sompolinsky, and Zohar Ringel.  Predicting the outputs of finite deep neural networks trained with noisy gradients.  Physical Review E, 104(6):064301, 2021.
  * [Nea96] Radford M Neal.  Priors for infinite networks.  In Bayesian Learning for Neural Networks, pages 29–53. Springer, 1996.
  * [NR21] Gadi Naveh and Zohar Ringel.  A self consistent theory of gaussian processes captures feature learning effects in finite cnns.  Advances in Neural Information Processing Systems, 34, 2021.
  * [PC21] Geoff Pleiss and John P Cunningham.  The limitations of large width in neural networks: A deep gaussian process perspective.  Advances in Neural Information Processing Systems, 34:3349–3363, 2021.
  * [PN21] Huy Tuan Pham and Phan-Minh Nguyen.  Global convergence of three-layer neural networks in the mean field regime.  ICLR, 2021.
  * [RVE18] Grant M Rotskoff and Eric Vanden-Eijnden.  Neural networks as interacting particle systems: Asymptotic convexity of the loss landscape and universal scaling of the approximation error.  arXiv preprint arXiv:1805.00915, 2018.
  * [RYH22] Daniel A Roberts, Sho Yaida, and Boris Hanin.  The Principles of Deep Learning Theory: An Effective Theory Approach to Understanding Neural Networks.  Cambridge University Press, 2022.
  * [SMG14] Andrew M Saxe, James L McClelland, and Surya Ganguli.  Exact solutions to the nonlinear dynamics of learning in deep linear neural networks.  ICLR, 2014.
  * [SR21] Inbar Seroussi and Zohar Ringel.  Separation of scales and a thermodynamic description of feature learning in some cnns.  arXiv preprint arXiv:2112.15383, 2021.
  * [SS20] Justin Sirignano and Konstantinos Spiliopoulos.  Mean field analysis of neural networks: A central limit theorem.  Stochastic Processes and their Applications, 130(3):1820–1852, 2020\.
  * [SS21] Justin Sirignano and Konstantinos Spiliopoulos.  Mean field analysis of deep neural networks.  Mathematics of Operations Research, 2021.
  * [WI20] Andrew G Wilson and Pavel Izmailov.  Bayesian deep learning and a probabilistic perspective of generalization.  Advances in neural information processing systems, 33:4697–4708, 2020.
  * [WRV+20] Florian Wenzel, Kevin Roth, Bastiaan S Veeling, Jakub Swikatkowski, Linh Tran, Stephan Mandt, Jasper Snoek, Tim Salimans, Rodolphe Jenatton, and Sebastian Nowozin.  How good is the bayes posterior in deep neural networks really?  arXiv preprint arXiv:2002.02405, 2020.
  * [Yai20] Sho Yaida.  Non-gaussian processes and neural networks at finite widths.  MSML, 2020.
  * [YH21] Greg Yang and Edward J Hu.  Tensor programs iv: Feature learning in infinite-width neural networks.  In International Conference on Machine Learning, pages 11727–11737. PMLR, 2021.
  * [YS20] Jiahui Yu and Konstantinos Spiliopoulos.  Normalization effects on shallow neural networks and related asymptotic expansions.  arXiv preprint arXiv:2011.10487, 2020.
  * [ZVP21] Jacob Zavatone-Veth and Cengiz Pehlevan.  Exact marginal prior distributions of finite bayesian neural networks.  Advances in Neural Information Processing Systems, 34, 2021.
  * [ZVTP22] Jacob A. Zavatone-Veth, William L. Tong, and Cengiz Pehlevan.  Contrasting random and learned features in deep bayesian linear regression.  Phys. Rev. E, 105:064118, 2022.