# Hessian Eigenspectra of More Realistic Nonlinear Models

Zhenyu Liao
ICSI and Department of Statistics
University of California, Berkeley, USA
<EMAIL_ADDRESS>Michael W. Mahoney
ICSI and Department of Statistics
University of California, Berkeley, USA
<EMAIL_ADDRESS>
###### Abstract

Given an optimization problem, the Hessian matrix and its eigenspectrum can be
used in many ways, ranging from designing more efficient second-order
algorithms to performing model analysis and regression diagnostics. When
nonlinear models and non-convex problems are considered, strong simplifying
assumptions are often made to make Hessian spectral analysis more tractable.
This leads to the question of how relevant the conclusions of such analyses
are for more realistic nonlinear models. In this paper, we exploit
deterministic equivalent techniques from random matrix theory to make a
_precise_ characterization of the Hessian eigenspectra for a broad family of
nonlinear models, including models that generalize the classical generalized
linear models, without relying on strong simplifying assumptions used
previously. We show that, depending on the data properties, the nonlinear
response model, and the loss function, the Hessian can have _qualitatively_
different spectral behaviors: of bounded or unbounded support, with single- or
multi-bulk, and with isolated eigenvalues on the left- or right-hand side of
the bulk. By focusing on such a simple but nontrivial nonlinear model, our
analysis takes a step forward to unveil the theoretical origin of many
visually striking features observed in more complex machine learning models.

##  1 Introduction

The Hessian matrix is ubiquitous in applied mathematics, statistics, and
machine learning (ML). Given a (loss) function $L(\mathbf{w})$ with some
parameters $\mathbf{w}\in{\mathbb{R}}^{p}$, the Hessian
$\mathbf{H}(\mathbf{w})\in{\mathbb{R}}^{p\times p}$ is defined as the second
derivative of the objective function with respect to the model parameter,
i.e., $\mathbf{H}(\mathbf{w})=\partial
L(\mathbf{w})/(\partial\mathbf{w}\cdot\partial\mathbf{w}^{\sf T})$. When a ML
model is being trained, it is common to parameterize that model by
$\mathbf{w}$, and then train that model by minimizing some (smooth) loss
function $L(\mathbf{w})$, with associated Hessian matrix
$\mathbf{H}(\mathbf{w})$, by backpropagating the error to improve $\mathbf{w}$
[40, 100]. Alternatively, once a ML model is trained, the Hessian matrix
$\mathbf{H}(\mathbf{w})$ (and the related Fisher information matrix [89, 91])
can be examined to identify outliers, perform diagnostics, and/or engage in
other sorts of uncertainty quantification and model validation [44, 102, 85].

For convex problems, the Hessian $\mathbf{H}(\mathbf{w})$ provides detailed
information on how to adjust the gradient to achieve improved convergence,
e.g., as in Newton-like methods. For non-convex problems, the properties of
the local loss “landscape” around a given point $\mathbf{w}$ in the parameter
space is of central significance [27, 56, 21, 61, 101, 103, 102]. In this
case, most obviously, the signs of the smallest and largest Hessian eigenvalue
can be used to test whether a given $\mathbf{w}$ is a local maximum, local
minimum, or a saddle point (i.e., the second-derivative test). More subtly,
the Hessian eigenvalue distribution characterizes the local curvature of the
loss function and provides direct access to, for instance, the fraction of
negative Hessian eigenvalues. This quantity determines the number of (local)
descent directions, a quantity that is directly connected to the rates of
convergence of various optimization algorithms [49, 23], as well as the large-
scale structure of the configuration space [13, 4].

For theoretical analysis of neural network (NN) models, Hessian eigenspectra
are often assumed to follow well-known random matrix distributions such the
Marc̆enko–Pastur law [69] or the Wigner’s semicircle law [93]. This enables
one to use Random Matrix Theory (RMT), but it involves (for NNs, at least)
making relatively strong simplifying assumptions [80, 81, 24]. A somewhat more
realistic setup involves using a so-called _spiked model_ (or a spiked
covariance model) [8, 11, 63]. In this case, the matrix follows a _signal-
plus-noise_ model and consists of _full rank_ random noise matrix and _low
rank_ statistical information structure. The “signal” eigenvalues are
generally larger than the noisy “bulk” eigenvalues; and the maximum
eigenvalues, when isolated from the bulk, are referred to as the “spikes.” A
substantial theory-practice gap exists, however. In both toy examples [41], as
well as state-of-the-art NN models [101, 103, 102, 100, 86, 30], the strong
simplifying assumptions are far from realistic or satisfactory. (A similar
theory-practice gap has been observed for other NN matrices to which RMT has
been applied, perhaps most notably weight matrices [70, 71].) A more precise
understanding of the Hessian eigenspectra (and its dependence on input data
structure, the underlying response model and model parameters, as well as the
loss function) for more realistic ML models is needed.

###  1.1 Our approach

In this article, we address these issues, in a setting that is simple enough
to be analytically tractable but complex enough to shed light on much more
realistic large-scale models. We consider a family of generalized generalized
linear models (G-GLMs) that extends the popular generalized linear model (GLM)
[29, 44]; and we show that, even for such simple models, the key simplifying
assumptions used in previous theoretical analyses of Hessian eigenspectra, can
be very inexact. In particular, apart from a few special cases (including
linear least squares and the logistic regression model in the case of
homogeneous features), most Hessians of G-GLM are _not_ close to the
Marc̆enko–Pastur and/or the semicircle law. Instead, the corresponding Hessian
depends in a complicated way on the interplay between the input feature
structure, the underlying response model, and the loss function. These
tradeoffs can be precisely characterized by the proposed analysis in the G-GLM
model we consider.

The G-GLM describes a generalized linear relation between the input feature
$\mathbf{x}_{i}\in{\mathbb{R}}^{p}$ and the corresponding response $y_{i}$, in
the sense that there exists some parameters
$\mathbf{w}_{*}\in{\mathbb{R}}^{p}$ such that for given $\mathbf{w}_{*}^{\sf
T}\mathbf{x}_{i}$, the response $y_{i}$ is independently drawn from

$y_{i}\sim f(y\mid\mathbf{w}_{*}^{\sf T}\mathbf{x}_{i}),$ (1)

for some conditional density function $f(\cdot\mid\cdot)$. This setup extends
GLMs such as the logistic model

$\mathbb{P}(y=1\mid\mathbf{w}_{*}^{\sf
T}\mathbf{x})=(1+e^{-\mathbf{w}_{*}^{\sf T}\mathbf{x}})^{-1},\quad
y\in\\{-1,1\\},$ (2)

and it covers a large family of models in many applications in statistics and
ML. Other examples of G-GLM models include:

  * •

the (noisy) nonlinear factor model [22], where
$y\sim{\mathcal{N}}(g(\mathbf{w}_{*}^{\sf T}\mathbf{x}),\sigma^{2})$ for some
nonlinear linking function $g:{\mathbb{R}}\to{\mathbb{R}}$ and $\sigma>0$;

  * •

the (noiseless) phase retrieval model [32], with $y=(\mathbf{w}_{*}^{\sf
T}\mathbf{x})^{2}$, in which case we wish to reconstruct $\mathbf{w}_{*}$ from
its magnitude measurements; and

  * •

the single-layer NN model $y=\sigma(\mathbf{w}_{*}^{\sf T}\mathbf{x})$, for
some nonlinear activation function $\sigma(t)$ such as the $\tanh$-sigmoid
$\sigma(t)=\tanh(t)$.

For a given training set $\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{n}$ of size $n$,
the standard approach to obtain/recover the parameter
$\mathbf{w}_{*}\in{\mathbb{R}}^{p}$ is to solve the following optimization
problem

$\textstyle\min_{\mathbf{w}}L(\mathbf{w})=\min_{\mathbf{w}}\frac{1}{n}\sum_{i=1}^{n}\ell(y_{i},\mathbf{w}^{\sf
T}\mathbf{x}_{i}),$ (3)

for some loss function
$\ell(y,h):{\mathbb{R}}\times{\mathbb{R}}\to{\mathbb{R}}$. This includes,
e.g., the negative log-likelihood of the observation model within the maximum
likelihood estimation framework [44], such as the logistic loss
$\ell(y,h)=\ln(1+e^{-yh})$ in the case of the logistic model (2). In many
applications, however, the optimization problem in (3) may not be convex [72,
94, 18],111As we shall see, in some non-convex cases of G-GLMs, the top
Hessian eigenvector can be shown to correlate positively with the sought-for
model parameter $\mathbf{w}_{*}$ and therefore be used as the initialization
of gradient descent methods [17, 57, 50]. This motivates our study of the
possible isolated Hessian eigenvalue-eigenvector pairs. and it can be NP-hard
in general (the noiseless phase retrieval model $y=(\mathbf{w}_{*}^{\sf
T}\mathbf{x})^{2}$ with the square loss $\ell(y,h)=(y-h^{2})^{2}$ is such an
example [17]).

###  1.2 Our main contributions

The main contribution of this work is the _exact asymptotic characterization_
of the Hessian eigenspectrum for the family of G-GLMs in (1), in the high
dimensional regime where the feature dimension $p$ and the sample size $n$ are
both large and comparable. Our main result is that we establish:

  1. 1.

the limiting eigenvalue distribution of the Hessian matrix (Theorem 1); and

  2. 2.

the behavior of (possible) isolated eigenvalues and the associated
eigenvectors (Theorem 2 and 3);

as a function of the dimension ratio $c=\lim p/n$, the feature (first- and
second-order) statistics, the loss function $\ell$ in (3), and the underlying
response model in (1). Our results are based on a technical result of
independent interest:

  1. 1.

a _deterministic equivalent_ (Theorem 4) of the random _resolvent_
$\mathbf{Q}(z)=(\mathbf{H}-z\mathbf{I}_{p})^{-1}$, for $z\in{\mathbb{C}}$ not
an eigenvalue of $\mathbf{H}$, of the generalized sample covariance:

$\mathbf{H}=\textstyle\frac{1}{n}\sum_{i=1}^{n}\ell^{\prime\prime}(y_{i},\mathbf{w}^{\sf
T}\mathbf{x}_{i})\mathbf{x}_{i}\mathbf{x}_{i}^{\sf
T}\equiv\frac{1}{n}\mathbf{X}\mathbf{D}\mathbf{X}^{\sf T},$ (4)

for
$\mathbf{D}\equiv\operatorname{diag}\\{\ell^{\prime\prime}(y_{i},\mathbf{w}^{\sf
T}\mathbf{x}_{i})\\}_{i=1}^{n}$,
$\mathbf{X}=[\mathbf{x}_{1},\ldots,\mathbf{x}_{n}]\in{\mathbb{R}}^{p\times n}$
and $\ell^{\prime\prime}(y,h)\equiv\partial^{2}\ell(y,h)/\partial h^{2}$. The
analysis is performed in the high dimensional limit as $n,p\to\infty$ with
$p/n\to c\in(0,\infty)$, under the setting of input Gaussian feature
$\mathbf{x}_{i}\sim{\mathcal{N}}(\boldsymbol{\mu},\mathbf{C})$, with some
statistical mean $\boldsymbol{\mu}\in{\mathbb{R}}^{p}$ and positive definite
covariance $\mathbf{C}\in{\mathbb{R}}^{p\times p}$. We also demonstrate our
results empirically by showing that:

  1. 1.

for a given response model (1), the Hessian eigenvalue distribution depends on
the choice of loss function and the data/feature statistics in an intrinsic
manner, e.g., bounded versus unbounded support (Figure 2) and single- versus
multi-bulk (Figure 3); and

  2. 2.

there may exist two _qualitatively_ different spikes—the well-known one due to
_data signal_ $\boldsymbol{\mu}$, and the other less-obvious one due to the
_underlying model itself_ —which may appear on different sides of the main
bulk, and their associated phase transition behaviors are different and are
characterized (Figure 6 versus 7).

To illustrate pictorially our contributions, we examine in Figure 1 the
Hessian eigenspectra as a function of the response model (1), the loss
function $\ell$, and the feature statistics $\boldsymbol{\mu},\mathbf{C}$
(that characterize the statistical property of the data/feature). In Figure
1(a) and 1(b), we compare the Hessian eigenvalues for the logistic model in
(2) with the (maximum likelihood) logistic loss $\ell(y,h)=\ln(1+e^{-yh})$,
for different choices of parameter $\mathbf{w}$. A nontrivial interplay
between the response model parameter $\mathbf{w}_{*}$, the feature statistics
$\boldsymbol{\mu},\mathbf{C}$, as well as $\mathbf{w}$ is reflected by the
range of the Hessian eigenvalue support and an additional right-hand spike in
Figure 1(b) (as opposed to Figure 1(a)), as confirmed by our theory.

For the phase retrieval model $y=(\mathbf{w}_{*}^{\sf T}\mathbf{x})^{2}$ with
square loss $\ell(y,h)=(y-h^{2})^{2}/4$, the non-convex nature of the problem
is reflected by a (relatively large) fraction of negative Hessian eigenvalues
in Figure 1(c). We also note that the top eigenvector (that corresponds to the
largest eigenvalue) contains structural information of the underlying model,
in the sense that it is positively correlated with $\mathbf{w}_{*}$, as
predicted by the theory in Figure 1(d), and that there exists an additional
negative eigenvalue. This is connected to the Hessian-based initialization
scheme discussed at the end of Section 1.1 and in Section 2.1 below widely
used in non-convex problems.

$6\cdot 10^{-2}$$8\cdot
10^{-2}$$0.1$$0.12$$0.14$$0.16$$0.18$$0.2$$0.22$$0.24$$0.26$$0.28$$0.3$$0.32$$0.34$$0.36$$0.38$$0.4$$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$5.5$$6$$6.5$EigenvaluesEmpiricalTheory
(a) Logistic, $\mathbf{w}=\mathbf{w}^{*}=\boldsymbol{\mu}$

$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$5.5$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$5.5$
(b) Logistic, $\mathbf{w}\sim{\mathcal{N}}(\mathbf{0},\mathbf{I}_{p}/p)$

$-18$$-16$$-14$$-12$$-10$$-8$$-6$$-4$$-2$$0$$2$$4$$6$$0$$0.1$$0.2$$-18$$-16$$-14$$-12$$-10$$-8$$-6$$-4$$-2$$0$$2$$4$$6$$0$$0.1$$0.2$$-18$$-16$$-14$$-12$$-10$$-8$$-6$$-4$$-2$$0$$2$$4$$6$$0$$0.1$$0.2$
(c) Phase retrieval model

Top eigenvector Theory (d) Top eigenvector in Figure 1(c)

Figure 1:  Illustration of our main results: eigenspectral properties of the
Hessian $\mathbf{H}(\mathbf{w})$ of G-GLM with $p=800$, $n=6\,000$ and
$\mathbf{C}=\mathbf{I}_{p}$. (Top): presence versus absence of a right-hand
side spike for different $\mathbf{w}$, logistic model (2) with logistic loss
$\ell(y,h)=\ln(1+e^{-yh})$ and
$\mathbf{w}_{*}=\boldsymbol{\mu}\sim{\mathcal{N}}(\mathbf{0},\mathbf{I}_{p}/p)$.
There is a nontrivial interplay between the response model (1), the feature
statistics $\boldsymbol{\mu},\mathbf{C}$, and $\mathbf{w}$, as reflected by
the “range” of the Hessian eigenvalues and an additional right-hand spike, in
Figure 1(b) versus 1(a), both captured by our theory. (Bottom): the Hessian
eigenspectrum may have a rather different shape (as opposed to the Marc̆enko-
Pastur-like in Figure 1(a) and 1(b)) for the phase retrieval model (left) and
the right-hand side isolated eigenvector is known here be an estimate of
$\mathbf{w}_{*}$ (right), as confirmed by our theory. With square loss
$\ell(y,h)=(y-h^{2})^{2}/4$,
$\mathbf{w}_{*}=[-2\cdot\mathbf{1}_{p/2};~{}2\cdot\mathbf{1}_{p/2}]/\sqrt{p}$,
$\mathbf{w}\sim{\mathcal{N}}(\mathbf{0},\mathbf{I}_{p}/p)$ and
$\boldsymbol{\mu}=\mathbf{0}$.

###  1.3 Notation and organization of the paper

##### Notations.

Throughout the paper, we follow the convention of denoting scalars by
lowercase, vectors by lowercase boldface, and matrices by uppercase boldface
letters. The notation $(\cdot)^{\sf T}$ denotes the transpose operator. The
norm $\|\cdot\|$ is the Euclidean norm for vectors and the spectral norm for
matrices. We use ${\xrightarrow{a.s.}}$ for almost sure convergence of random
variables. We use $\Im[\cdot]$ to denote the imaginary part of a complex
number and $\imath$ to present the imaginary unit.

##### Organization of the paper.

We discuss related previous efforts and preliminaries in Section 2. Our main
technical results on the limiting Hessian eigenvalue distribution, as well as
on the behavior of the Hessian isolated eigenvalues and associated
eigenvectors are presented in Section 3. In Section 4, further discussions on
the consequences of these technical results are made, together with numerical
evaluations. We conclude the article in Section 5 with some possible future
extensions.

##  2 Background and Preliminaries

In this section, we discuss related previous efforts in Section 2.1 and some
(technical) preliminaries on the spiked random matrix model and the
deterministic equivalent technique in Section 2.2 and Section 2.3,
respectively.

###  2.1 Related work

##### Random matrix theory.

Random matrices of the type (4) are related to the so-called _separable
covariance model_ [105, 26]. This model is of the form
$\frac{1}{n}\mathbf{C}^{\frac{1}{2}}\mathbf{Z}\mathbf{D}\mathbf{Z}^{\sf
T}\mathbf{C}^{\frac{1}{2}}$, for $\mathbf{Z}\in{\mathbb{R}}^{p\times n}$
having, say, i.i.d. standard Gaussian entries, and _deterministic_ or
_independent_ matrices $\mathbf{C}\in{\mathbb{R}}^{p\times p}$ and
$\mathbf{D}\in{\mathbb{R}}^{n\times n}$. Our results generalize this, in the
sense that we allow the matrix $\mathbf{D}$ to _depend_ on the random matrix
$\mathbf{Z}$, in a possibly nonlinear fashion, per (4). This is of direct
interest for the Hessian matrix of the G-GLM. We particularly show that, in
the case where $\mathbf{Z}$ is a standard Gaussian random matrix, while the
limiting spectral measure remains the same _as if $\mathbf{D}$ was
independent_ (Theorem 1), spurious isolated eigenvalue may appear due to this
dependence between $\mathbf{D}$ and $\mathbf{Z}$ (Theorem 2 and 3). This is
directly related to the popular spectral initialization method used in non-
convex optimization [67].

##### Hessian eigenspectra in ML models.

The eigenspectra of Hessian matrices arising in ML models (in particular, for
NNs) have attracted considerable interest recently, both from a theoretical
perspective [83, 84, 19, 33, 41, 95, 80, 39, 48] and from an empirical
perspective [30, 86, 102, 100]. For example, the Hessian eigenspectrum
provides a characterization of the (local) loss curvature and therefore the
so-called “sharpness” (or “flatness”) of a trained model.222Some claim that
this corresponds to the generalization performance of the model [45, 58, 19,
28]. However, these investigations are either limited to empirical evaluation
[83, 84] or built upon somewhat unrealistic simplifying assumptions that
reduce to the “mixed” behavior of Marc̆enko–Pastur and semicircle law [80,
24], or to that of the product of independent Gaussian matrices [41]. In
contrast, we focus on the more tractable example of G-GLM in (1), and we
provide _precise_ results on the Hessian eigenspectra for structural Gaussian
feature on arbitrary loss functions. Also, instead of focusing solely on the
limiting eigenvalue distribution, our results also shed novel light on the
isolated eigenvalues (above and/or below the bulk) that are empirically
observed in modern NNs [83, 33, 41, 65, 75].

##### Hessian eigenspectra: a statistical physics perspective.

Statistical physicists are interested in the energy “landscape” of random
functions, by means of studying the Hessian of such functions [36, 38, 37, 98,
35]. This is somewhat different from the Hessian of the loss function in
optimization, and it is mainly in the context of spin glass theory. Depending
on the model, the associated Hessian eigenspectra may establish a phase
transition behavior and take either a Marc̆enko-Pastur or a semicircle shape,
but a wide range of other phenomena have also been studied. These are known to
have connections with ML systems such as NN models in high dimensions [1,
104].

##### Spectral initialization in non-convex problems.

For non-convex problems, the performance of gradient-based methods can be very
sensitive to initialization. A popular initialization method is the _spectral
initialization_ method, which can be traced back to the Principal Hessian
Directions method [64]. Similar techniques were then proposed to provide a
“good” initialization for gradient descent in non-convex problems such as
phase retrieval [32, 74, 16], matrix completion [57], low-rank matrix recovery
[50], blind deconvolution [62], and sparse coding [3]. In [67], which was
generalized by [73], the authors evaluated the eigenspectrum asymptotics of
the generalized sample covariance matrix
$\frac{1}{n}\sum_{i=1}^{n}f(y_{i})\mathbf{x}_{i}\mathbf{x}_{i}^{\sf T}$, for
some processing function $f:{\mathbb{R}}\to{\mathbb{R}}$ and
$\mathbf{x}_{i}\sim{\mathcal{N}}(\mathbf{0},\mathbf{I}_{p})$. Their technical
approach is, however, limited to the case of very homogeneous features,
$\mathbf{x}_{i}\sim{\mathcal{N}}(\mathbf{0},\mathbf{I}_{p})$, and it was
particularly focused on the phase retrieval model, instead of the Hessian
(that depends on the loss) more generally. In a sense, we generalize the
analysis in [67] to the Hessian matrix of G-GLM, by developing a systematic
approach to account for the data structures ($\boldsymbol{\mu}\neq\mathbf{0}$
and $\mathbf{C}\neq\mathbf{I}_{p}$).

##### Scalable second-order methods.

Second order methods are among the most powerful optimization methods that
have been designed, and there have been several attempts to use their many
advantages for machine learning applications, including both theoretical
developments [97, 92, 82], as well as for training and/or quantizing state-of-
the-art NNs [102, 100, 86, 30, 96]. We expect that our precise
characterization of the Hessian sheds new light on the improved design of more
efficient second-order methods. For example, insights from the proposed
analysis could shed light on improved optimization methods for particular
classes of models, as well as improved strategies for model post-processing
and diagnostics.

###  2.2 Spiked random matrix model

Define the _empirical spectral measure_ of a random matrix
$\mathbf{M}\in{\mathbb{R}}^{p\times p}$

$\mu_{\mathbf{M}}\equiv\frac{1}{p}\sum_{i=1}^{p}\delta_{\lambda_{i}(\mathbf{M})}$
(5)

as the normalized counting measure of the eigenvalues
$\lambda_{i}(\mathbf{M})$ of $\mathbf{M}$. If the random $\mu_{\mathbf{M}}$
converges to a (non-random) probability measure $\mu$ as $p\to\infty$, we say
$\mu$ is the _limiting spectral measure_ of the random matrix $\mathbf{M}$.
Examples of such limiting laws include the popular Marc̆enko-Pastur law [69]
and Wigner’s semicircle law [93].

In many applications, instead of the (limiting) eigenvalue distribution, one
may be more interested in the behavior of some particular eigenvalue-
eigenvector pairs. For a standard (noise-like) random matrix $\mathbf{M}$, the
term “ _spiked random matrix model_ ” [51, 11, 6] is used to refer to the _low
rank_ update of $\mathbf{M}$, for instance the additive perturbation of the
type $\mathbf{M}+\mathbf{A}$ for some (deterministic) $\mathbf{A}$ such that
$\operatorname{rank}(\mathbf{A})=K$ is fixed as $p\to\infty$. This type of
random matrix model plays a fundamental role in the analysis of principle
component analysis (PCA) [52, 53], sparse PCA [2, 68, 15] and Canonical
Correlation Analysis (CCA) [47, 10, 54] in high dimensions. As a concrete
example, for
$\mathbf{x}_{i}\sim{\mathcal{N}}(\boldsymbol{\mu},\mathbf{I}_{p})$ we can
write
$\mathbf{X}=[\mathbf{x}_{1},\ldots,\mathbf{x}_{n}]=\boldsymbol{\mu}\mathbf{1}_{n}^{\sf
T}+\mathbf{Z}$ for $\mathbf{Z}\in{\mathbb{R}}^{p\times n}$ having i.i.d.
standard Gaussian entries, so that the sample covariance matrix
$\frac{1}{n}\mathbf{X}\mathbf{X}^{\sf T}=\frac{1}{n}\mathbf{Z}\mathbf{Z}^{\sf
T}+\boldsymbol{\mu}\boldsymbol{\mu}^{\sf
T}+\frac{1}{n}\boldsymbol{\mu}\mathbf{1}_{n}^{\sf T}\mathbf{Z}^{\sf
T}+\frac{1}{n}\mathbf{Z}\mathbf{1}_{n}\boldsymbol{\mu}^{\sf T}$ follows a
spiked model that is an (at most) rank-two additive perturbation of the
standard Wishart random matrix $\frac{1}{n}\mathbf{Z}\mathbf{Z}^{\sf T}$. The
spiked model is particularly popular and useful in charactering signal-plus-
noise type behavior in the features, by focusing on the isolated eigenvalues
(i.e., the “spikes”) outside the main eigenspectrum bulk. Precisely, if the
signal strength $\|\boldsymbol{\mu}\|$ is larger than a certain threshold
(that depends on the dimension ratio $p/n$), an informative eigenvalue
$\hat{\lambda}$ of $\frac{1}{n}\mathbf{X}\mathbf{X}^{\sf T}$ can be observed
to “jump out” from the noisy main bulk, with the associated eigenvector
$\hat{\mathbf{u}}$ “aligned” to $\boldsymbol{\mu}$, in the sense that the
“cosine-similarly” $\hat{\mathbf{u}}^{\sf
T}\boldsymbol{\mu}/\|\boldsymbol{\mu}\|\gg 0$ as $n,p\to\infty$.333Recall that
for high dimensional random vector $\mathbf{v}\in{\mathbb{R}}^{p}$, say with
i.i.d. Gaussian entries of zero-mean and $1/p$ variance (so that
$\|\mathbf{v}\|=O(1)$ with high probability), we have, for
$\|\boldsymbol{\mu}\|=O(1)$ that $\mathbf{v}^{\sf T}\boldsymbol{\mu}\to 0$ as
$p\to\infty$; that is, a random vector is _almost orthogonal_ to any
deterministic or _independent_ random vector in high dimensions. This type of
_phase transition_ behavior has been extensively studied in the RMT literature
[8, 11, 63]. In this article, we show two _different_ types of spikes and
their corresponding phase transitions: the first well-known spike due to data
signal $\boldsymbol{\mu}$, which follows the classical BPP phase transition
[8]; and a second less-well-known spike, which arises from the underlying
(nonlinear) response model (1), even in the absence of any signal, and which
follows a very different phase transition behavior. This is discussed at
length with numerical evaluations in Section 4.

###  2.3 Deterministic equivalent

In this work, a unified analysis framework is proposed to assess
simultaneously the limiting spectral measure and the behavior of the
(possible) isolated eigenvalue-eigenvector pairs of the Hessian-type matrix in
(4). This analysis is based on the so-called deterministic equivalent
technique [43, 25] that provides asymptotic characterization of the
(expectation of the) resolvent
$\mathbf{Q}(z)\equiv(\mathbf{H}-z\mathbf{I}_{p})^{-1}$ of the random Hessian
matrix $\mathbf{H}\in{\mathbb{R}}^{p\times p}$. With additional concentration
of measure arguments [60, 76], the limiting spectral measure, the isolated
eigenvalue location, as well as the eigenvector projection onto a given
deterministic vector, while all being random, become asymptotically
deterministic and fully accessible via the analysis of $\mathbf{Q}(z)$ in the
$n,p\to\infty$ limit. More precisely, we have the following.

  1. 1.

The normalized trace $\frac{1}{p}\operatorname{tr}\mathbf{Q}(z)$ characterizes
the Stieltjes transform $m_{\mu}(z)=\int(t-z)^{-1}\mu(dt)$ of the limiting
spectral measure $\mu$ of $\mathbf{H}$, from which $\mu$ can be recovered, for
$a,b$ continuity points of $\mu$, via
$\mu([a,b])=\frac{1}{\pi}\lim_{\varepsilon\downarrow
0}\int_{a}^{b}\Im[m_{\mu}(x+\imath\varepsilon)]\,dx$.

  2. 2.

Solving the determinant equation $\det(\mathbf{H}-\lambda\mathbf{I}_{p})=0$
for $\lambda\in{\mathbb{R}}$ outside the limiting spectral support
$\operatorname{supp}(\mu)$ (that is, not in the bulk eigenspectrum) reveals
the asymptotic location of the isolated eigenvalue $\hat{\lambda}$ of
$\mathbf{H}$.

  3. 3.

For $(\hat{\lambda},\hat{\mathbf{u}})$ an eigenvalue-eigenvector pair of
interest of $\mathbf{H}$, with Cauchy’s integral formula we have
$|\mathbf{w}^{\sf
T}\hat{\mathbf{u}}|^{2}=-\frac{1}{2\pi\imath}\oint_{\Gamma_{\lambda}}\mathbf{w}^{\sf
T}\mathbf{Q}(z)\mathbf{w}\,dz$, for a deterministic vector
$\mathbf{w}\in{\mathbb{R}}^{p}$ and $\Gamma_{\lambda}$ a positively oriented
contour surrounding _only_ $\hat{\lambda}$.

As such, for $\bar{\mathbf{Q}}(z)$ a deterministic equivalent of
$\mathbf{Q}(z)$, that is,

$\mathbf{Q}(z)\leftrightarrow\bar{\mathbf{Q}}(z),\quad\text{with~{}}\frac{1}{p}\operatorname{tr}\mathbf{A}(\mathbf{Q}(z)-\bar{\mathbf{Q}}(z))\to
0,\text{~{}and~{}}\mathbf{a}^{\sf
T}(\mathbf{Q}(z)-\bar{\mathbf{Q}}(z))\mathbf{b}\to 0,$ (6)

almost surely as $n,p\to\infty$, for $\mathbf{A}\in{\mathbb{R}}^{n\times n}$
and $\mathbf{a},\mathbf{b}\in{\mathbb{R}}^{n}$ of bounded (Euclidean and
spectral) norms, the limiting spectral measure (via the associated Stieltjes
transform) and the behavior of isolated eigenpairs are directly accessible via
the study of $\bar{\mathbf{Q}}(z)$.

##  3 Main Results

In this section, we present our main results: on the limiting Hessian
eigenspectrum (in Section 3.1); and on the behavior of the (possible) isolated
eigenvalue-eigenvector(s) (in Section 3.2). These two main results depend on a
key technical deterministic equivalent result for the Hessian resolvent (in
Section 3.3), which is of independent interest. We position ourselves in the
following high dimensional regime.

######  Assumption 1 (High dimensional asymptotics).

As $n,p\to\infty$ with $p/n\to c\in(0,\infty)$, we have
$\max\\{\|\mathbf{w}\|,\|\mathbf{w}_{*}\|\\}=O(1)$ and
$\mathbf{x}_{i}\overset{i.i.d.}{\sim}{\mathcal{N}}(\boldsymbol{\mu},\mathbf{C})$
with $\max\\{\|\boldsymbol{\mu}\|,\|\mathbf{C}\|\\}=O(1)$.

###  3.1 Limiting spectral measure

Our first result concerns the limiting eigenvalue distribution of the Hessian
matrix, given under the form of its Stieltjes transform. This is a direct
consequence of our main technical Theorem 4 in Section 3.3, and it is proven
in Appendix A.2.

######  Theorem 1 (Limiting spectral measure).

Let Assumption 1 hold. We have, as $n,p\to\infty$ with $p/n\to
c\in(0,\infty)$, the empirical spectral measure $\mu_{\mathbf{H}}$ of the
Hessian matrix $\mathbf{H}$ defined in (4) converges weakly and almost surely
to a probability measure $\mu$, defined through its Stieltjes transform
$m(z)=\int(t-z)^{-1}\mu(dt)$ as the unique solution to444Uniqueness is ensured
in such a way that $\Im[m(z)]\cdot\Im[z]>0$ for $\Im[z]\neq 0$ and $zm(z)<0$
for $\Im[z]=0$, so that $(z,m(z))$ is a valid Stieltjes transform couple; see
[43].

$\textstyle
m(z)=\frac{1}{p}\operatorname{tr}\bar{\mathbf{Q}}_{b}(z),\quad\delta(z)=\frac{1}{n}\operatorname{tr}\left(\mathbf{C}\bar{\mathbf{Q}}_{b}(z)\right),$
(7)

for
$\bar{\mathbf{Q}}_{b}(z)\equiv({\mathbb{E}}[\frac{g\cdot\mathbf{C}}{1+g\cdot\delta(z)}]-z\mathbf{I}_{p})^{-1}$
and random variable

$g\equiv\partial^{2}\ell(y,h)/\partial h^{2},~{}h=\mathbf{w}^{\sf
T}\mathbf{x}\sim\mathcal{N}(\mathbf{w}^{\sf T}\boldsymbol{\mu},\mathbf{w}^{\sf
T}\mathbf{C}\mathbf{w}),$ (8)

with $y$ and $\ell$ defined respectively in (1) and (3). Moreover, if we
denote $\nu$ the law of $g$ and assume the empirical spectral measure of
$\mathbf{C}$ converges to $\tilde{\nu}$ as $p\to\infty$, then (7) can be
expressed as

$m(z)=\int\Big{(}-z+\tilde{t}\int\frac{t}{1+t\delta(z)}\,\nu(dt)\Big{)}^{-1}\,\tilde{\nu}(d\tilde{t}),\quad\delta(z)=\int\frac{c\tilde{t}}{-z+\tilde{t}\int\frac{t}{1+t\delta(z)}\,\nu(dt)}\,\tilde{\nu}(d\tilde{t}).$
(9)

A few remarks on the consequences of Theorem 1 are in order.

######  Remark 1 (Connection to other RMT models).

In the form (9), the (Stieltjes transform of the) limiting spectral measure
$\mu$ is determined by the dimension ratio $c=\lim p/n$ and the two
(probability) measures $\nu$ and $\tilde{\nu}$. As already briefly discussed
in Section 2.1, this formulation is closely connected to the separable
covariance model [59, 14, 78, 26, 99]. Moreover, if $\nu(dt)=\delta_{1}(t)$ is
a Dirac mass at one, this reduces to the classical sample covariance model
[87]; taking further $\tilde{\nu}(dt)=\delta_{1}(t)$ gives the Marc̆enko-
Pastur law. See Section 4.1 for numerical evaluations of these special cases.

In particular, the support of the (limiting) Hessian eigenvalue distribution
$\mu$ is directly linked to that of $\nu$ and $\tilde{\nu}$.

######  Remark 2 (Hessian eigen-support).

Under Assumption 1, the (limiting) spectral measure $\tilde{\nu}$ of the
covariance $\mathbf{C}$ has bounded support. However, this may not be the case
for $\nu$, the law of the random variable $g$ defined in (8). Since the
Hessian eigenvalue distribution $\mu$ is of compact support _if and only if_
both $\nu$ and $\tilde{\nu}$ have compact support [26, Porposition 3.4], $\mu$
may be of unbounded support, depending on that of $\nu$, and thus on the
underlying model and loss function.

An example of unbounded $\mu$ is the phase retrieval model with
$y=(\mathbf{w}_{*}^{\sf T}\mathbf{x})^{2}$ and square loss
$\ell(y,h)=(y-h^{2})^{2}/4$, for which we have $g=3(\mathbf{w}^{\sf
T}\mathbf{x})^{2}-(\mathbf{w}_{*}^{\sf T}\mathbf{x})^{2}$ for
$\mathbf{x}\sim{\mathcal{N}}(\boldsymbol{\mu},\mathbf{C})$. As a consequence,
with say $\mathbf{w}_{*}=\mathbf{w}$, $g$ follows a chi-square distribution
with one degree of freedom, and has thus _unbounded support_. This corresponds
to Figure 1(c), where the Hessian spectrum has a “heavier” tail compared to
Figure 1(a). In this case, the empirically observed “isolated” eigenvalue is
essentially due to a finite-dimensional effect and will be “buried” in the
noisy main bulk as $n,p$ grow large. Therefore, aiming for an (almost surely)
isolated eigenvalue-eigenvector (e.g., to recover the model parameter
$\mathbf{w}_{*}$ using the top Hessian eigenvector) for $n,p$ large, some
preprocessing function $f$ must be applied. This has been discussed in
previous work [67, 73] and corresponds to the so-called trimming strategy
widely used in phase retrieval [20] (particularly in the “linear measurements”
$n\sim p$ regime), with for instance the truncation function
$f(t)=1_{|t|\leq\epsilon}$ for some $\epsilon>0$.

Another example of unbounded $\mu$ appears when the exponential loss function
[34] is used. Precisely, consider the logistic model (2) with the choice of
loss function $\ell(y,h)=\exp(-yh)$. In this case, we have that $g=\exp(-yh)$
for $h\sim\mathcal{N}(\mathbf{w}^{\sf T}\boldsymbol{\mu},\mathbf{w}^{\sf
T}\mathbf{C}\mathbf{w})$, which follows a log-normal distribution and has
_unbounded support_. As such, the limiting Hessian eigenvalue distribution
$\mu$ has also unbounded support. On the other hand, with the logistic loss
$\ell(y,h)=\ln(1+e^{-yh})$, one has $g\leq 1/4$ and $\mu$ is guaranteed to
have bounded support. See Figure 2, where the empirical Hessian eigenvalues
and the limiting distributions are compared for logistic and exponential
losses, with a more “heavy-tailed” behavior observed for the exponential loss.

$0.1$$0.2$$0.3$$0.4$$0$$2$$4$Eigenvalues | $0$$2$$4$$6$$0$$0.2$$0.4$$0.6$
---|---
Figure 2:  Impact of loss function: bounded (left) versus unbounded (right)
Hessian eigenvalue distribution, with $p=800$, $n=6\,000$, logistic model in
(2) with $\boldsymbol{\mu}=\mathbf{0}$, $\mathbf{C}=\mathbf{I}_{p}$,
$\mathbf{w}_{*}=\mathbf{0}$ and
$\mathbf{w}=[-\mathbf{1}_{p/2},~{}\mathbf{1}_{p/2}]/\sqrt{p}$. Emphasis on the
_bounded_ Hessian eigen-support with the logistic loss (left) versus the
_unbounded_ support with the exponential loss (right).

Clearly, depending on the measures $\nu$ (of the random variable $g$, which
depends on the feature statistics $\boldsymbol{\mu},\mathbf{C}$, the loss and
the underlying model) and $\tilde{\nu}$ (of the eigenvalue distribution of
feature covariance $\mathbf{C}$), the Hessian spectrum can have very different
forms. See Figure 3, where we compare the empirical Hessian eigenvalues with
their limiting behaviors, given in Theorem 1 for different feature covariance
structures.555Covariance describes the joint variability or the “correlation”
between entries of the feature vector, and it is of great significance in the
analysis of image (with local structure) and time series data. In particular,
we may observe a single main bulk with more “compact” Hessian eigenvalues as
in Figure 3-(left) or multiple bulks (two in the case of Figure 3-right) with
Hessian eigenvalues more “spread-out”, depending on the feature covariance
structure and the measure $\tilde{\nu}$. In the form of (9), the condition for
the existence of multi-bulk eigenspectrum has been throughly discussed in [26,
Section 3.2–3.4] and can be numerically evaluated with ease.

As a side remark, the “multi-bulk” behavior similar to Figure 3-(right) has
been empirically observed in Hessians of NNs in [65, 75], suggesting the
(quantitative) effectiveness of the proposed theory beyond simple G-GLM.

$0$$0.2$$0.4$$0.6$$0$$2$$4$Eigenvalues | $0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$0$$2$$4$
---|---
Figure 3:  Impact of feature covariance: Hessian spectrum of single- (left)
versus multi-bulk (right) , with $p=800$, $n=6\,000$, logistic model in (2)
with $\mathbf{w}^{*}=\mathbf{0}_{p}$,
$\mathbf{w}=\boldsymbol{\mu}\sim{\mathcal{N}}(\mathbf{0},\mathbf{I}_{p}/p)$;
for
$\mathbf{C}=\operatorname{diag}[\mathbf{1}_{p/2};~{}2\cdot\mathbf{1}_{p/2}]$
(left) and
$\mathbf{C}=\operatorname{diag}[\mathbf{1}_{p/2};~{}4\cdot\mathbf{1}_{p/2}]$
(right).

###  3.2 Isolated eigenvalues and eigenvectors

As discussed in Remark 2, the Hessian has bounded (limiting) eigen-support if
and only if $\nu$, the law of $g$, has bounded support. Under this assumption
(or, after the application of some function $f$ so that $f(g)$ is bounded), we
can then talk about the (possible) isolated Hessian eigenvalues, as in the
following result, the proof of which is in Appendix A.3.

######  Theorem 2 (Isolated eigenvalues).

In the setting of Theorem 1, assume that the probability measure $\nu$ of the
random variable $g$ defined in (8) is of bounded supported, and define

$\mathbf{G}(z)=\mathbf{I}_{3}+\boldsymbol{\Lambda}(z)\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}_{b}(z)\mathbf{V}\in{\mathbb{R}}^{3\times 3},$ (10)

with $\bar{\mathbf{Q}}_{b}(z),\delta(z)$ defined in (7),
$\mathbf{V}\equiv[\boldsymbol{\mu},~{}\mathbf{C}\mathbf{w}_{*},~{}\mathbf{C}\mathbf{w}]\in{\mathbb{R}}^{p\times
3}$,
$\mathbf{U}\equiv\mathbf{C}^{\frac{1}{2}}[\mathbf{w}_{*},~{}\mathbf{w}]\in{\mathbb{R}}^{p\times
2}$ and

$\boldsymbol{\Lambda}(z)\equiv{\mathbb{E}}\left[\frac{g}{1+g\cdot\delta(z)}\begin{bmatrix}1&(\mathbf{U}^{+}\mathbf{z})^{\sf
T}\\\
\mathbf{U}^{+}\mathbf{z}&\mathbf{U}^{+}\mathbf{z}(\mathbf{U}^{+}\mathbf{z})^{\sf
T}-(\mathbf{U}^{\sf T}\mathbf{U})^{+}\end{bmatrix}\right],$ (11)

for
$\mathbf{z}=\mathbf{C}^{-\frac{1}{2}}(\mathbf{x}-\boldsymbol{\mu})\sim\mathcal{N}(\mathbf{0},\mathbf{I}_{p})$,
where we denote $\mathbf{U}^{+}$ the Moore-Penrose pseudoinverse of
$\mathbf{U}$ with $\mathbf{U}^{+}=(\mathbf{U}^{\sf
T}\mathbf{U})^{-1}\mathbf{U}^{\sf T}$ for invertible $\mathbf{U}^{\sf
T}\mathbf{U}$. Then, for $\lambda$ such that $\mathbf{G}(\lambda)$ has a zero
eigenvalue (of multiplicity one), there exists an eigenvalue $\hat{\lambda}$
of $\mathbf{H}$ such that $\hat{\lambda}-\lambda~{}{\xrightarrow{a.s.}}~{}0$.

######  Remark 3 (Connection to isolated eigenvalue phase transition).

Theorem 2 provides an asymptotic characterization of the possible isolated
Hessian eigenvalues by computing the determinant of the much smaller (three-
by-three) deterministic matrix $\mathbf{G}$ closely related to the key
quantity $\delta(z)$. Note that Theorem 2 does not provide, at least
explicitly, the _phase transition_ condition under which these
eigenvalues/spikes become “isolated” from the main bulk. As we shall see in
more detail in Section 4, two types of quantitatively different phase
transitions can be characterized: one due to the data “signal”
$\boldsymbol{\mu}$; and one due to the structure of the underlying model (even
in the absence of any signal in the data).

We can also analyze the associated eigenvectors. First note that, for
$n\to\infty$ with $p$ fixed, we have, by the strong law of large numbers, that
$\mathbf{H}~{}{\xrightarrow{a.s.}}~{}{\mathbb{E}}[\mathbf{H}]$, with

$\displaystyle{\mathbb{E}}[\mathbf{H}]$
$\displaystyle={\mathbb{E}}[\ell^{\prime\prime}(y,\mathbf{w}^{\sf
T}\mathbf{x})\mathbf{x}\mathbf{x}^{\sf
T}]={\mathbb{E}}[g]\cdot\mathbf{C}+\mathbf{V}\begin{bmatrix}1&{\mathbb{E}}[g\cdot\mathbf{U}^{+}\mathbf{z}]^{\sf
T}\\\
{\mathbb{E}}[g\cdot\mathbf{U}^{+}\mathbf{z}]&\mathbf{U}^{+}{\mathbb{E}}[g\cdot(\mathbf{z}\mathbf{z}^{\sf
T}-\mathbf{I}_{p})](\mathbf{U}^{+})^{\sf T}\end{bmatrix}\mathbf{V}^{\sf T}$

with
$\mathbf{z}=\mathbf{C}^{-\frac{1}{2}}(\mathbf{x}-\boldsymbol{\mu})\sim\mathcal{N}(\mathbf{0},\mathbf{I}_{p})$
and $g$ defined in (8). As a consequence, it is expected that in the large
$n,p\to\infty$ limit, the top eigenvectors of $\mathbf{H}$ should be related
to the columns of $\mathbf{V}$. This is the case in Figure 1(d), where the top
eigenvector is observed to be a “noisy” version of the model parameter
$\mathbf{w}_{*}$. More precisely, for $(\hat{\lambda},\hat{\mathbf{u}})$ an
isolated eigenpair of $\mathbf{H}$, the _random_ projection $\mathbf{V}^{\sf
T}\hat{\mathbf{u}}\hat{\mathbf{u}}^{\sf T}\mathbf{V}\in{\mathbb{R}}^{3\times
3}$ can be shown to be asymptotically close to a _deterministic_ matrix as
$n,p\to\infty$. This measures the “cosine-similarly” between the Hessian
isolated eigenvector $\hat{\mathbf{u}}$ with any column of $\mathbf{V}$ and
consequently the performance of using $\hat{\mathbf{u}}$ as an estimate of,
for instance the model parameter $\mathbf{w}_{*}$ for
$\mathbf{C}=\mathbf{I}_{p}$. This result is given in the following theorem,
which is proven in Appendix A.4.

######  Theorem 3 (Isolated eigenvectors).

In the setting of Theorem 2, for an isolated eigenvalue-eigenvector pair
$(\hat{\lambda},\hat{\mathbf{u}})$ of $\mathbf{H}$ and $\lambda$ the
asymptotic position given in Theorem 2, it holds that

$\mathbf{V}^{\sf T}\hat{\mathbf{u}}\hat{\mathbf{u}}^{\sf
T}\mathbf{V}=-\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}_{b}(\lambda)\mathbf{V}\cdot\Xi(\lambda)+o(1),$

with
$\Xi(\lambda)=\frac{\mathbf{v}_{r,\mathbf{G}}\mathbf{v}_{l,\mathbf{G}}^{\sf
T}}{\mathbf{v}_{l,\mathbf{G}}^{\sf
T}\mathbf{G}^{\prime}(\lambda)\mathbf{v}_{r,\mathbf{G}}}$, for
$\bar{\mathbf{Q}}_{b}(z)$ and $\mathbf{G}(z)$ defined in (7) and (10),
$\mathbf{v}_{l,\mathbf{G}},\mathbf{v}_{r,\mathbf{G}}\in{\mathbb{R}}^{3}$ the
left and right eigenvectors of $\mathbf{G}(\lambda)$ associated with
eigenvalue zero, and $\mathbf{G}^{\prime}(z)$ the derivative of
$\mathbf{G}(z)$ with respect to $z$ evaluated at $z=\lambda$.

###  3.3 Technical tool: deterministic equivalent

As already discussed in Section 2.3, our main technical tool to derive
Theorems 1, 2 and 3 is the so-called deterministic equivalent approach [43,
25] for the Hessian resolvent
$\mathbf{Q}(z)=(\mathbf{H}-z\mathbf{I}_{p})^{-1}$, from which follows our main
results on the Hessian limiting eigenvalue distribution and the behavior of
the possible isolated eigenpairs. This result is given in the following
theorem, with the proof deferred to Appendix A.1.

######  Theorem 4 (Deterministic equivalent).

Let $\mathbf{Q}(z)\equiv(\mathbf{H}-z\mathbf{I}_{p})^{-1}$ be the resolvent of
the Hessian $\mathbf{H}$ defined in (4). Then, under Assumption 1, we have, as
$n,p\to\infty$ with $p/n\to c\in(0,\infty)$, that
$\mathbf{Q}(z)\leftrightarrow\bar{\mathbf{Q}}(z)$, with

$\bar{\mathbf{Q}}^{-1}(z)=\textstyle{\mathbb{E}}\frac{g}{1+g\cdot\delta(z)}(\mathbf{C}^{\frac{1}{2}}(\mathbf{I}_{p}-\mathbf{P}_{\mathbf{U}})\mathbf{C}^{\frac{1}{2}}+\boldsymbol{\alpha}\boldsymbol{\alpha}^{\sf
T})-z\mathbf{I}_{p},$

for random vector
$\boldsymbol{\alpha}\equiv\boldsymbol{\mu}+\mathbf{C}^{\frac{1}{2}}\mathbf{P}_{\mathbf{U}}\mathbf{z}\in{\mathbb{R}}^{p}$
and $g=\ell^{\prime\prime}(y,\mathbf{w}^{\sf
T}\boldsymbol{\mu}+\mathbf{w}^{\sf T}\mathbf{C}^{\frac{1}{2}}\mathbf{z})$ for
$\mathbf{z}\sim{\mathcal{N}}(\mathbf{0},\mathbf{I}_{p})$ as in (8), $y$
defined in (1), $\mathbf{P}_{\mathbf{U}}$ the projection onto the subspace
spanned by the columns of $\mathbf{U}$ so that
$\mathbf{P}_{\mathbf{U}}=\mathbf{U}(\mathbf{U}^{\sf
T}\mathbf{U})^{-1}\mathbf{U}^{\sf T}$ for $\mathbf{U}^{\sf T}\mathbf{U}$
invertible, and $\delta(z)$ the unique solution to (7).

##  4 Empirical Evaluation and Discussion

In this section, we provide further discussion on the consequences of Theorems
1, 2 and 3, together with numerical evaluations. The implications of Theorem 1
on the (limiting) Hessian eigenvalue distribution is discussed in Section 4.1.
Then, we discuss the consequences of Theorems 2 and 3 on possible isolated
eigenpairs, and we show that two fundamentally different phase transition
behaviors occur. Recall that classical spiked models, extensively studied in
RMT literature [8, 11, 63], exhibit the following: (i) the isolated spike
appears due to the presence of some statistical “signal” in the data; and (ii)
a “monotonic” phase transition behavior, as a function of the signal strength,
can be characterized. Here, we show that, in addition to this classical
“signal” spike, another type of spike may arise due to the underlying G-GLM
model itself (i.e., $\mathbf{w}_{*}$ and $\mathbf{w}$), and it exhibits a very
different transition behavior. In particular, in Section 4.2, we discuss
eigenvalue spikes due to data signal; and in Section 4.3, we discuss
eigenvalue spikes due to model structure, even when there is no signal in the
data.

###  4.1 Hessian eigenvalue distribution

Here, for a better interpretation of Theorem 1 on the Hessian spectrum, we
consider the special case where $\mathbf{C}=\mathbf{I}_{p}$, in which case (7)
can be further simplified as
$z\delta(z)=1-c-{\mathbb{E}}[1/(1+g\cdot\delta(z))]$, with $m(z)=\delta(z)/c$
and $c=\lim p/n$.

Let us start with the simplest setting where the random variable $g$ in (8) is
constant, say $g=1$. This happens, for instance, in the case where the square
loss $\ell(y,h)=(y-h)^{2}/2$ is employed for the logistic model (2). In this
case, the Hessian is _independent_ of the parameters
$\mathbf{w},\mathbf{w}_{*}$, and we have the Stieltjes transform $m(z)$ as the
solution to

$zcm^{2}(z)-(1-c-z)m(z)+1=0,$ (12)

that corresponds to the popular Marc̆enko-Pastur law.

As long as $g$ is _not_ constant, the limiting Hessian spectrum is, a priori,
different from the Marc̆enko-Pastur law, even in the simple
$\mathbf{C}=\mathbf{I}_{p}$ setting, since the associated Stieltjes transform
$m(z)$ is different from the solution of (12). However, we see in Figure
4-(top left) that, for the logistic model (2) with logistic loss, the Hessian
spectrum is close, at least visually, to a rescaled Marc̆enko-Pastur law. This
is due to the fact that, the distribution of $g$ is more “concentrated”
(around some constant, see Figure 4-(top right) versus (bottom right) for a
comparison between different cases) in this specific model. This is in sharp
contrast to, for example, Figure 4-(bottom) where with the exponential loss,
the law of $g$ has a much larger “spread”, and the Hessian eigenspectrum is
therefore far away from a Marc̆enko-Pastur-shape.

$0$$0.1$$0.2$$0.3$$0.4$EigenvaluesRescaled MP | $0$$0.1$$0.2$Distribution of $g$Empirical $g$Scale factor
---|---
$0$$2$$4$$6$$8$Eigenvalues | $0$$2$$4$$6$$8$$10$Distribution of $g$
Figure 4:  Comparison of Hessian eigenspectra with the Marc̆enko-Pastur law in
the setting of Figure 2. (Top): Marc̆enko-Pastur-like Hessian with logistic
loss; (left) Hessian eigenvalues and (right) empirical distribution of the
$g_{i}$s versus the scaling factor (empirically obtained from the distance
between the maximum and minimum Hessian eigenvalues to fit the Marc̆enko-
Pastur law). Comparing to Figure 2-(right) and Figure 3, we see that changes
in the choice of loss and/or feature statistics (e.g., heterogeneous features)
result in very different and non-Marc̆enko-Pastur-like Hessian spectral
behavior. (Bottom): an example of a very non-Marc̆enko-Pastur-like Hessian
with exponential loss (left) and the associated $g_{i}$s (right). Note that
the scales of the axes are different in different subfigures.

It should be mentioned that this “empirical fit” has already been observed in
[79, Figure 5], where acceleration methods proposed for a Marc̆enko-Pastur
distributed Hessian (in linear least squares) work reasonably well on logistic
regression problems. Our theory proposes a convincing explanation of this
empirical observation on logistic regression, and possibly for others more
involved models such as NNs.

However, it must be pointed out that this “visual approximation” by Marc̆enko-
Pastur law is _not robust_ , in the sense that it “visually” holds only for
(yet still is formally different from) the case of (i) logistic model with
(ii) logistic loss and (iii) identity covariance $\mathbf{C}=\mathbf{I}_{p}$.
Any change in the response model (e.g., with the phase retrieval model in
Figure 1(c)), the choice of loss function (e.g., with the exponential loss in
Figure 2-right), or going beyond the simple identity covariance (that models
very homogeneous features) setting, as in Figure 3, would induce a Hessian
spectrum that is _very_ different from the Marc̆enko-Pastur law. In this vein,
our Theorem 1 goes beyond such a “loose” Marc̆enko-Pastur approximation. It
provides a first example toward the understanding of Hessians of deep NNs and
other much more realistic models that exploit nonlinear transformations (such
as activation function) and feature statistics.

In the proofs of Theorem 1–4, we explicitly used the fact that, for jointly
Gaussian random variables, uncorrelatedness implies independence, and
therefore our main results (on the limiting Hessian eigenvalue distribution
and the possible isolated eigenpairs) are proven to hold only for Gaussian
features. On the other hand, universality phenomena commonly arise in random
matrix theory, and it has been proven for various models that the limiting
spectral measure [88], the spectrum edges [12, 99], and even involved
functionals of the eigenvalues [7] _do not_ depend on the distribution of the
random matrix, so long as its entries are independent and have zero mean, unit
variance, and finite higher, say fourth, order moment.666Counterexample
exists, though, for instance it has been shown in [31, 66] that, even in the
presence of signal, nonlinear random matrices may produce “noisy spikes” (that
are purely random and carry no information about the data signal and) that, in
particular, depends on the fourth moment of the distribution.

Figure 5 compares the Hessian eigenvalues for feature following Gaussian,
Bernoulli, and Student-t distribution, where we empirically confirm that
Theorem 1 provides consistent predictions for non-Gaussian features with
reasonably large $n,p$. We conjecture that Theorem 1, proven here only for
Gaussian distribution, holds more generally.

$0$$0.2$$0.4$$0.6$$0$$2$$4$$6$$8$Eigenvalues | $0$$0.2$$0.4$$0.6$$0$$2$$4$$6$$8$ | $0$$0.2$$0.4$$0.6$$0$$2$$4$$6$$8$
---|---|---
Figure 5:  While proved here only in the Gaussian case, the empirical Hessian
eigenvalue distribution appears rather “stable” beyond Gaussian distribution:
with $p=800$, $n=1\,200$, logistic model in (2) with
$\mathbf{w}_{*}=\mathbf{0}_{p}$,
$\mathbf{w}=\boldsymbol{\mu}\sim{\mathcal{N}}(\mathbf{0},\mathbf{I}_{p}/p)$
and $\mathbf{C}=\mathbf{I}_{p}$, for Gaussian (left), symmetric Bernoulli
(middle), and Student-t distribution with $7$ degrees of freedom (right).

$0$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$0.55$$0.6$$0.65$$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$5.5$EigenvaluesEigenvaluesLimiting
lawRight edge$0$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$0.55$$0.6$$0.65$$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$5.5$EigenvaluesEigenvaluesLimiting
lawRight edge

$\sqrt{c}$11.5$0$$0.02$$0.04$$0.06$$0.08$$0.1$$0.12$$0.14$$0.16$$0.18$$0.2$$\|\boldsymbol{\mu}\|^{2}$EigengapEmpirical
Theoretical

$\sqrt{c}$11.5$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$\|\boldsymbol{\mu}\|^{2}$Eigen-
alignment

Figure 6:  Spike (on right-hand side of bulk) due to data signal in Corollary
1: theory versus practice of (left) Hessian eigenspectrum with
$\|\boldsymbol{\mu}\|^{2}=0.8$, (middle) eigengap
$\operatorname{dist}(\lambda_{\boldsymbol{\mu}},\operatorname{supp}(\mu))$,
and (right) top eigenvector alignment $\alpha$ in (14), as a function of the
signal strength $\|\boldsymbol{\mu}\|^{2}$, on logistic model with logistic
loss, for $\boldsymbol{\mu}\propto[-\mathbf{1}_{p/2},~{}\mathbf{1}_{p/2}]$,
$\mathbf{w}=\mathbf{w}_{*}=\mathbf{0}$, $\mathbf{C}=\mathbf{I}_{p}$, $p=512$
and $n=2\,048$. Results averaged over $50$ runs.

$0$$2\cdot 10^{-2}$$4\cdot 10^{-2}$$6\cdot 10^{-2}$$8\cdot
10^{-2}$$0.1$$0.12$$0.14$$0.16$$0.18$$0.2$$0.22$$0.24$$0.26$$0.28$$0.3$$0.32$$0.34$$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$EigenvaluesEigenvaluesLimiting
lawLeft edge$0$$2\cdot 10^{-2}$$4\cdot 10^{-2}$$6\cdot 10^{-2}$$8\cdot
10^{-2}$$0.1$$0.12$$0.14$$0.16$$0.18$$0.2$$0.22$$0.24$$0.26$$0.28$$0.3$$0.32$$0.34$$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$EigenvaluesEigenvaluesLimiting
lawLeft edge

$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$5.5$$6$$6.5$$7$$7.5$$8$$0$$0.01$$0.01$$0.02$$0.02$$0.03$$0.03$$0.04$$\|\mathbf{w}\|$EigengapEmpirical
Theoretical

$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$5.5$$6$$6.5$$7$$7.5$$8$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$\|\mathbf{w}\|$Eigen-
alignment

Figure 7:  Spike (on left-hand side of bulk) due to response model in
Corollary 2 _in the absence of data signal_ : (left) Hessian eigenspectrum for
$\|\mathbf{w}\|=2$, with an emphasis on the _left_ isolated eigenvalue
$\hat{\lambda}_{\mathbf{w}}$, (middle) eigengap
$\operatorname{dist}(\lambda_{\mathbf{w}},\operatorname{supp}(\mu))$, and
(right) dominant eigenvector alignment, as a function $\|\mathbf{w}\|$ with
$\mathbf{w}\propto[-\mathbf{1}_{p/2},~{}\mathbf{1}_{p/2}]$,
$\mathbf{w}_{*}=\boldsymbol{\mu}=\mathbf{0}$, $\mathbf{C}=\mathbf{I}_{p}$,
$p=800$ and $n=8\,000$. Results averaged over $50$ runs.

###  4.2 Isolated eigenvalues: Spike due to data signal

To study the spike due to data “signal” $\boldsymbol{\mu}$, as well as its
phase transition behavior, we focus here on the case
$\mathbf{w}_{*}=\mathbf{w}=\mathbf{0}$. This, in the case of logistic model
(2), gives rise to a much simpler form of limiting eigenspectrum (per Theorem
1) and possible isolated eigenpairs (per Theorem 2 and 3), as summarized in
the following corollary, the proof of which is given in Appendix B.3.

######  Corollary 1 (Spike due to data signal: logistic model).

Consider the logistic model in (2) with logistic loss, for
$\mathbf{w}=\mathbf{w}_{*}=\mathbf{0}$ and $\mathbf{C}=\mathbf{I}_{p}$, we
have $g=1/4$ and the limiting Hessian eigenvalue distribution is the
Marc̆enko-Pastur law, but rescaled by a factor of $1/4$. Moreover, it follows
from Theorem 2 and 3 that there is _at most one_ isolated eigenpair
$(\hat{\lambda}_{\boldsymbol{\mu}},\hat{\mathbf{u}}_{\boldsymbol{\mu}})$ of
$\mathbf{H}$ and it satisfies

$\hat{\lambda}_{\boldsymbol{\mu}}~{}{\xrightarrow{a.s.}}~{}\left\\{\begin{array}[]{ll}\lambda_{\boldsymbol{\mu}}=\frac{1}{4}(1+\rho+c\cdot\frac{\rho+1}{\rho})&\rho>\sqrt{c},\\\
\frac{1}{4}(1+\sqrt{c})^{2}&\rho\leq\sqrt{c};\end{array}\right.$ (13)

and

$\frac{|\boldsymbol{\mu}^{\sf
T}\hat{\mathbf{u}}_{\boldsymbol{\mu}}|^{2}}{\|\boldsymbol{\mu}\|^{2}}~{}{\xrightarrow{a.s.}}~{}\left\\{\begin{array}[]{ll}\alpha=\frac{\rho^{2}-c}{\rho^{2}+c\rho}&\rho>\sqrt{c},\\\
0&\rho\leq\sqrt{c};\end{array}\right.$ (14)

with the signal strength $\rho=\lim\|\boldsymbol{\mu}\|^{2}$ and $c=\lim p/n$.

The isolated eigenvalue-eigenvector behavior described in Corollary 1 follows
the classical phase transition [9, 8, 77]. In particular, it exhibits the
following properties: (i) the isolated eigenvalue always appears on the right-
hand side of the main (Marc̆enko-Pastur) bulk; and (ii) the eigenvalue
location and associated eigenvector alignment is “monotonic” with respect to
the signal strength $\|\boldsymbol{\mu}\|^{2}$. By “monotonic,” we mean that,
for a fixed dimension ratio $c$, the largest Hessian eigenvalue is bound to
become ‘isolated” once $\|\boldsymbol{\mu}\|^{2}$ exceeds $\sqrt{c}$, and also
that its value, as well as the eigenvector alignment, increase monotonically
as $\|\boldsymbol{\mu}\|^{2}$ grows. This behavior is confirmed in Figure 6.

As we shall see below in Section 4.3, this behavior is _not_ exhibited for a
spike due to the model parameter $\mathbf{w}$.

###  4.3 Isolated eigenvalues: Spike due to model structure

To study the spike due to the underlying model structure (i.e.,
$\mathbf{w}_{*}$ and $\mathbf{w}$), as well as its phase transition behavior,
we position ourselves in the situation where $\boldsymbol{\mu}=\mathbf{0}$,
that is, in the absence of any data signal, and we consider again the logistic
model. This leads to the following corollary, the proof of which is given in
Appendix B.4.

######  Corollary 2 (Spike due to model: logistic model).

Consider the logistic model in (2) with logistic loss,
$\boldsymbol{\mu}=\mathbf{0}$ and $\mathbf{C}=\mathbf{I}_{p}$, then the
Stieltjes transform $m(z)$ satisfies $m(z)=1/({\mathbb{E}}[f(r,z)]-z)$ for
$f(t,z)=1/(cm(z)+2+e^{-t}+e^{t})$ and $r\sim\mathcal{N}(0,\|\mathbf{w}\|^{2})$
that _depends_ on $\mathbf{w}$ but _not_ on $\mathbf{w}_{*}$. Moreover, there
is _at most one_ isolated eigenvalue $\hat{\lambda}_{\mathbf{w}}$ of
$\mathbf{H}$ that is due to $\mathbf{w}$ and satisfies
$\hat{\lambda}_{\mathbf{w}}-\lambda_{\mathbf{w}}~{}{\xrightarrow{a.s.}}~{}0$
with $\lambda_{\mathbf{w}}$ solution to
$0=\det\mathbf{G}(\lambda_{\mathbf{w}})=1+m(\lambda_{\mathbf{w}})\frac{{\mathbb{E}}[f(r,\lambda_{\mathbf{w}})(r^{2}-\|\mathbf{w}\|^{2})]}{\|\mathbf{w}\|^{2}}$.

The situation here (for a spike due to model) is more subtle than for the
spike due to data signal (discussed in Section 4.2): as $\mathbf{w}$ changes
(e.g., as the “energy” $\|\mathbf{w}\|$ grows), both the Hessian (limiting)
eigenvalue distribution and the possible spike location are impacted. Figure 7
illustrates the behavior of the spike due to $\mathbf{w}$ in the setting of
Corollary 2. First, note that, different from the case of spike due to data
signal $\boldsymbol{\mu}$, the spike in Figure 7-(left) appears on the _left-
hand side_ of the main bulk. In particular, this means that the Hessian may
admit an eigenvalue that is _significantly smaller_ than all the other
eigenvalues. Second, note from Figure 7-(middle) that, different from the
spike due to $\boldsymbol{\mu}$, the spike due to $\mathbf{w}$ exhibits here a
“non-monotonic” behavior in the sense that, it is absent for small values of
$\|\mathbf{w}\|$ (as for small $\|\boldsymbol{\mu}\|$ in Figure 6-middle) and
it becomes “isolated” as $\|\mathbf{w}\|$ increases, but then again gets
closer to the main bulk, as the norm $\|\mathbf{w}\|$ continues to increase,
resulting a unimodal eigengap.

Finally, it is perhaps even more surprising to observe in Figure 7-(right)
that, the alignment between the associated isolated eigenvector and the
parameter $\mathbf{w}$ is, unlike the eigengap in Figure 7-(middle),
_monotonically increasing_ as $\|\mathbf{w}\|$ grows large. This is as in the
case of spike due to data signal $\boldsymbol{\mu}$, illustrated in Figure
6-(right). This particularly suggests that, in the case of spike due to model,
a _smaller_ eigengap may not always imply _less_ statistical “information”
contained in the associated eigenvector, which is somehow contrary to the
conventional eigengap heuristic [90, 55].

We already see that there are spikes due to model on the _left-hand side_ of
the main bulk (see, e.g., Figure 7-left). One may wonder if it is always the
case that the “model spikes” are smaller than the main bulk of eigenvalues. It
turns out that this is _not_ true, and these spikes may appear on the right-
hand side of the bulk.

In Figure 8, we evaluate the noiseless phase retrieval model
$y=(\mathbf{w}_{*}^{\sf T}\mathbf{x})^{2}$ with the square loss
$\ell(y,h)=(y-h^{2})^{2}/4$. As already mentioned in Remark 2, the Hessian
eigenvalue distribution of phase retrieval model has unbounded support, due to
the unboundedness of the support of the chi-square distribution. As a result,
some preprocessing function must be applied for a spike to be observed. Here
in Figure 8 we particularly focus on the preprocessing function

$f(t)=\frac{\max(t,0)-1}{\max(t,0)+\sqrt{2/c}-1}$ (15)

proposed in [73] in the case of $\mathbf{w}\propto\mathbf{w}_{*}$,
$\boldsymbol{\mu}=\mathbf{0}$ and $\mathbf{C}=\mathbf{I}_{p}$. Here, we
observe a spike on the _right-hand side_ of the main bulk, the “eigengap” of
which follows again a non-monotonic behavior as a function of the norm
$\|\mathbf{w}_{*}\|$. Perhaps more surprisingly, in contrast to Figure
7-(right) where the eigenvector alignment is monotonically increasing, here we
see in Figure 8-(right) that the eigenvector alignment follows a similar _non-
monotonic_ trend as the eigengap.

$-0.2$$-0.1$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$1.4$$1.6$$1.8$$2$$2.2$$2.4$Eigenvalue
distribution$-0.2$$-0.1$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$1.4$$1.6$$1.8$$2$$2.2$$2.4$Eigenvalue
distribution

$0.2$$0.4$$0.6$$0.8$$1$$1.2$$1.4$$1.6$$1.8$$2$$0$$0.02$$0.04$$0.06$$0.08$$0.1$$0.12$$0.14$$0.16$$0.18$$0.2$$\|\mathbf{w}_{*}\|$Eigengap

$0.2$$0.4$$0.6$$0.8$$1$$1.2$$1.4$$1.6$$1.8$$2$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$\|\mathbf{w}_{*}\|$Eigen-
alignment

Figure 8:  Spike (on right-hand side of bulk) due to response model _in the
absence of data signal_ : (left) Hessian eigenspectrum for
$\|\mathbf{w}_{*}\|=1$, with an emphasis on the _right_ isolated eigenvalue
$\hat{\lambda}_{\mathbf{w}_{*}}$, (middle) eigengap, and (right) dominant
eigenvector alignment, as a function $\|\mathbf{w}_{*}\|$ with
$\mathbf{w}_{*}\propto[-\mathbf{1}_{p/2},~{}\mathbf{1}_{p/2}]$,
$\mathbf{w}=\sqrt{2/3}\cdot\mathbf{w}_{*}$, $\boldsymbol{\mu}=\mathbf{0}$,
$\mathbf{C}=\mathbf{I}_{p}$, $p=800$ and $n=4\,000$. Results averaged over
$50$ runs.

As a side remark, note in Figure 8 that negative eigenvalues are observed for
$\mathbf{w}\propto\mathbf{w}_{*}$: this is not surprising since the modified
matrix $\mathbf{H}$ with the preprocessing function $f(t)$ in (15) _no longer_
corresponds to the Hessian of the original optimization problem.

##  5 Conclusions

In this paper, we precisely characterized the Hessian eigenspectra of the
nonlinear G-GLMs, in the high dimensional regime where both feature dimension
and sample size are large and comparable. By focusing on this more realistic,
yet still tractable, family of models, our results provide theoretical
justifications for many empirical observations in large-scale machine learning
models such as deep neural networks. The deterministic equivalent technique
that we employ is very powerful, and we expect that it will enable analysis of
a range of more realistic machine learning models. Of particular future
interest would be to use this technique to understand heavy-tailed structure
in the correlations of state-of-the-art-neural network models [70, 71], and
the connections with recent results on heavy-tailed behavior in stochastic
optimization algorithms that are used to train these models [46, 42].

#### Acknowledgments

We would like to acknowledge DARPA, IARPA (contract W911NF20C0035), NSF, and
ONR via its BRC on RandNLA for providing partial support of this work. Our
conclusions do not necessarily reflect the position or the policy of our
sponsors, and no official endorsement should be inferred.

## References

  * [1] Madhu Advani, Subhaneil Lahiri and Surya Ganguli  “Statistical mechanics of complex neural systems and high dimensional data”  In _Journal of Statistical Mechanics: Theory and Experiment_ 2013.03, 2013, pp. P03014  DOI: 10.1088/1742-5468/2013/03/p03014
  * [2] Arash A. Amini and Martin J. Wainwright  “High-dimensional analysis of semidefinite relaxations for sparse principal components”  In _The Annals of Statistics_ 37.5B, 2009, pp. 2877–2921
  * [3] Sanjeev Arora, Rong Ge, Tengyu Ma and Ankur Moitra  “Simple, efficient, and neural algorithms for sparse coding”  In _Journal of Machine Learning Research_ 40.2015 Microtome Publishing, 2015
  * [4] Antonio Auffinger, Gérard Ben Arous and Jiří Černý  “Random Matrices and Complexity of Spin Glasses”  In _Communications on Pure and Applied Mathematics_ 66.2, 2013, pp. 165–201  DOI: 10.1002/cpa.21422
  * [5] Zhidong Bai and Jack W Silverstein  “Spectral analysis of large dimensional random matrices”  Springer, 2010
  * [6] Zhidong Bai and Jianfeng Yao  “On sample eigenvalues in a generalized spiked population model”  In _Journal of Multivariate Analysis_ 106 Elsevier, 2012, pp. 167–177
  * [7] Zhidong D Bai and Jack W Silverstein  “CLT for linear spectral statistics of large-dimensional sample covariance matrices”  In _Advances In Statistics_ World Scientific, 2008, pp. 281–333
  * [8] Jinho Baik, Gérard Ben Arous and Sandrine Péché  “Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices”  In _The Annals of Probability_ 33.5 Institute of Mathematical Statistics, 2005, pp. 1643–1697
  * [9] Jinho Baik and Jack W Silverstein  “Eigenvalues of large sample covariance matrices of spiked population models”  In _Journal of multivariate analysis_ 97.6 Elsevier, 2006, pp. 1382–1408
  * [10] Zhigang Bao, Jiang Hu, Guangming Pan and Wang Zhou  “Canonical correlation coefficients of high-dimensional Gaussian vectors: Finite rank case”  In _The Annals of Statistics_ 47.1, 2019, pp. 612–640  DOI: 10.1214/18-aos1704
  * [11] Florent Benaych-Georges and Raj Rao Nadakuditi  “The eigenvalues and eigenvectors of finite, low rank perturbations of large random matrices”  In _Advances in Mathematics_ 227.1 Elsevier, 2011, pp. 494–521
  * [12] Paul Bourgade, László Erdös and Horng-Tzer Yau  “Edge universality of beta ensembles”  In _Communications in Mathematical Physics_ 332.1 Springer, 2014, pp. 261–353
  * [13] Alan J. Bray and David S. Dean  “Statistics of Critical Points of Gaussian Fields on Large-Dimensional Spaces”  In _Physical Review Letters_ 98.15, 2007, pp. 150201  DOI: 10.1103/physrevlett.98.150201
  * [14] Zdzisław Burda, Jerzy Jurkiewicz and Bartłomiej Wacław  “Spectral moments of correlated Wishart matrices”  In _Physical Review E_ 71.2 APS, 2005, pp. 026111
  * [15] Tony Cai, Zongming Ma and Yihong Wu  “Optimal estimation and rank detection for sparse spiked covariance matrices”  In _Probability Theory and Related Fields_ 161.3-4, 2015, pp. 781–815
  * [16] Emmanuel J Candes, Thomas Strohmer and Vladislav Voroninski  “Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming”  In _Communications on Pure and Applied Mathematics_ 66.8 Wiley Online Library, 2013, pp. 1241–1274
  * [17] Emmanuel J. Candès, Xiaodong Li and Mahdi Soltanolkotabi  “Phase Retrieval via Wirtinger Flow: Theory and Algorithms”  In _IEEE Transactions on Information Theory_ 61.4, 2015, pp. 1985–2007  DOI: 10.1109/tit.2015.2399924
  * [18] Olivier Chapelle, Choon Teo, Quoc Le and Alex Smola  “Tighter bounds for structured estimation”  In _Advances in neural information processing systems_ 21, 2008, pp. 281–288
  * [19] Pratik Chaudhari et al.  “Entropy-SGD: Biasing gradient descent into wide valleys”  In _Journal of Statistical Mechanics: Theory and Experiment_ 2019.12 IOP Publishing, 2019, pp. 124018
  * [20] Yuxin Chen and Emmanuel Candes  “Solving random quadratic systems of equations is nearly as easy as solving linear systems”  In _Advances in Neural Information Processing Systems_ , 2015, pp. 739–747
  * [21] Yuejie Chi, Yue M. Lu and Yuxin Chen  “Nonconvex Optimization Meets Low-Rank Matrix Factorization: An Overview”  In _IEEE Transactions on Signal Processing_ 67.20, 2018, pp. 5239–5269  eprint: 1809.09573
  * [22] Dennis Child  “The essentials of factor analysis”  Cassell Educational, 1990
  * [23] Yacine Chitour, Zhenyu Liao and Romain Couillet  “A Geometric Approach of Gradient Descent Algorithms in Neural Networks”  In _arXiv_ , 2018  eprint: 1811.03568
  * [24] Anna Choromanska et al.  “The loss surfaces of multilayer networks”  In _Artificial intelligence and statistics_ , 2015, pp. 192–204
  * [25] Romain Couillet and Merouane Debbah  “Random matrix methods for wireless communications”  Cambridge University Press, 2011
  * [26] Romain Couillet and Walid Hachem  “Analysis of the limiting spectral measure of large random matrices of the separable covariance type”  In _Random Matrices: Theory and Applications_ 3.04 World Scientific, 2014, pp. 1450016
  * [27] Yann N Dauphin et al.  “Identifying and attacking the saddle point problem in high-dimensional non-convex optimization”  In _Advances in Neural Information Processing Systems_ 27, NIPS’14  Curran Associates, Inc., 2014, pp. 2933–2941
  * [28] Laurent Dinh, Razvan Pascanu, Samy Bengio and Yoshua Bengio  “Sharp Minima Can Generalize For Deep Nets”  In _International Conference on Machine Learning_ , 2017, pp. 1019–1028
  * [29] Annette J Dobson and Adrian G Barnett  “An introduction to generalized linear models”  CRC press, 2018
  * [30] Zhen Dong et al.  “HAWQ: Hessian aware quantization of neural networks with mixed-precision”  In _Proceedings of the IEEE International Conference on Computer Vision_ , 2019, pp. 293–302
  * [31] Zhou Fan and Andrea Montanari  “The spectral norm of random inner-product kernel matrices”  In _Probability Theory and Related Fields_ 173.1-2 Springer, 2019, pp. 27–85
  * [32] James R. Fienup  “Phase retrieval algorithms: a comparison”  In _Applied Optics_ 21.15, 1982, pp. 2758  DOI: 10.1364/ao.21.002758
  * [33] Stanislav Fort and Surya Ganguli  “Emergent properties of the local geometry of neural loss landscapes”  In _arXiv preprint arXiv:1910.05929_ , 2019
  * [34] Yoav Freund, Robert Schapire and Naoki Abe  “A short introduction to boosting”  In _Journal-Japanese Society For Artificial Intelligence_ 14.771-780 JAPANESE SOC ARTIFICIAL INTELL, 1999, pp. 1612
  * [35] Yan V Fyodorov and Pierre Le Doussal  “Hessian spectrum at the global minimum of high-dimensional random landscapes”  In _Journal of Physics A: Mathematical and Theoretical_ 51.47, 2018, pp. 474002  DOI: 10.1088/1751-8121/aae74f
  * [36] Yan V. Fyodorov  “Complexity of Random Energy Landscapes, Glass Transition, and Absolute Value of the Spectral Determinant of Random Matrices”  In _Physical Review Letters_ 92.24, 2004, pp. 240601  DOI: 10.1103/physrevlett.92.240601
  * [37] Yan V. Fyodorov and Pierre Le Doussal  “Manifolds Pinned by a High-Dimensional Random Landscape: Hessian at the Global Energy Minimum”  In _Journal of Statistical Physics_ 179.1, 2020, pp. 176–215  DOI: 10.1007/s10955-020-02522-2
  * [38] Yan V. Fyodorov and Ian Williams  “Replica Symmetry Breaking Condition Exposed by Random Matrix Calculation of Landscape Complexity”  In _Journal of Statistical Physics_ 129.5-6, 2007, pp. 1081–1116  DOI: 10.1007/s10955-007-9386-x
  * [39] Mario Geiger et al.  “Jamming transition as a paradigm to understand the loss landscape of deep neural networks”  In _Physical Review E_ 100.1 APS, 2019, pp. 012115
  * [40] Ian Goodfellow, Yoshua Bengio and Aaron Courville  “Deep Learning”  MIT Press, 2016
  * [41] Diego Granziol  “Beyond Random Matrix Theory for Deep Networks”  In _arXiv preprint arXiv:2006.07721_ , 2020
  * [42] Mert Gurbuzbalaban, Umut Simsekli and Lingjiong Zhu  “The Heavy-Tail Phenomenon in SGD”  In _arXiv_ , 2020  eprint: 2006.04740
  * [43] Walid Hachem, Philippe Loubaton and Jamal Najim  “Deterministic equivalents for certain functionals of large random matrices”  In _The Annals of Applied Probability_ 17.3 Institute of Mathematical Statistics, 2007, pp. 875–930
  * [44] Trevor Hastie, Robert Tibshirani and Jerome Friedman  “The elements of statistical learning: data mining, inference, and prediction”  Springer Science & Business Media, 2009
  * [45] Sepp Hochreiter and Jürgen Schmidhuber  “Flat minima”  In _Neural Computation_ 9.1 MIT Press, 1997, pp. 1–42
  * [46] Liam Hodgkinson and Michael W. Mahoney  “Multiplicative noise and heavy tails in stochastic optimization”  In _arXiv_ , 2020  eprint: 2006.06293
  * [47] Harold Hotelling  “Relations Between Two Sets of Variates”  In _Biometrika_ 28.3/4, 1936, pp. 321
  * [48] Arthur Jacot, Franck Gabriel and Clement Hongler  “The asymptotic spectrum of the Hessian of DNN throughout training”  In _International Conference on Learning Representations_ , 2019
  * [49] Prateek Jain and Purushottam Kar  “Non-convex Optimization for Machine Learning”  In _Foundations and Trends® in Machine Learning_ 10.3-4 Now Publishers Inc. Hanover, MA, USA, 2017, pp. 142–336
  * [50] Prateek Jain, Praneeth Netrapalli and Sujay Sanghavi  “Low-rank matrix completion using alternating minimization”  In _Proceedings of the forty-fifth annual ACM symposium on Theory of computing_ , 2013, pp. 665–674
  * [51] Iain M Johnstone  “On the distribution of the largest eigenvalue in principal components analysis”  In _Annals of statistics_ JSTOR, 2001, pp. 295–327
  * [52] Iain M Johnstone and Debashis Paul  “PCA in high dimensions: An orientation”  In _Proceedings of the IEEE_ 106.8 IEEE, 2018, pp. 1277–1292
  * [53] Iain M. Johnstone and Arthur Yu Lu  “Sparse Principal Components Analysis”  In _arXiv_ , 2009  URL: https://arxiv.org/abs/0901.4392
  * [54] Iain M. Johnstone and Alexei Onatski  “Testing in high-dimensional spiked models”  In _Annals of Statistics_ 48.3, 2020, pp. 1231–1254
  * [55] Antony Joseph and Bin Yu  “Impact of regularization on spectral clustering”  In _The Annals of Statistics_ 44.4 Institute of Mathematical Statistics, 2016, pp. 1765–1791
  * [56] Kenji Kawaguchi  “Deep Learning without Poor Local Minima”  In _Advances in Neural Information Processing Systems_ 29, NIPS’16  Curran Associates, Inc., 2016, pp. 586–594
  * [57] Raghunandan H Keshavan, Andrea Montanari and Sewoong Oh  “Matrix completion from a few entries”  In _IEEE transactions on information theory_ 56.6 IEEE, 2010, pp. 2980–2998
  * [58] Nitish Shirish Keskar et al.  “On large-batch training for deep learning: Generalization gap and sharp minima”  In _arXiv preprint arXiv:1609.04836_ , 2016
  * [59] A Khorunzhy  “Eigenvalue distribution of large random matrices with correlated entries”  In _Mat. Fiz. Anal. Geom_ 3.1-2, 1996, pp. 80–101
  * [60] Michel Ledoux  “The concentration of measure phenomenon”  American Mathematical Soc., 2001
  * [61] Jason D. Lee et al.  “First-order methods almost always avoid strict saddle points”  In _Mathematical Programming_ 176.1-2, 2019, pp. 311–337
  * [62] Kiryung Lee, Yanjun Li, Marius Junge and Yoram Bresler  “Blind recovery of sparse signals from subsampled convolution”  In _IEEE Transactions on Information Theory_ 63.2 IEEE, 2016, pp. 802–821
  * [63] Thibault Lesieur, Florent Krzakala and Lenka Zdeborová  “Phase transitions in sparse PCA”  In _2015 IEEE International Symposium on Information Theory (ISIT)_ , 2015, pp. 1635–1639  IEEE
  * [64] Ker-Chau Li  “On principal Hessian directions for data visualization and dimension reduction: Another application of Stein’s lemma”  In _Journal of the American Statistical Association_ 87.420 Taylor & Francis, 1992, pp. 1025–1039
  * [65] Xinyan Li et al.  “Hessian based analysis of SGD for deep nets: Dynamics and generalization”  In _Proceedings of the 2020 SIAM International Conference on Data Mining_ , 2020, pp. 190–198  SIAM
  * [66] Zhenyu Liao, Romain Couillet and Michael W Mahoney  “Sparse quantized spectral clustering”  In _arXiv preprint arXiv:2010.01376_ , 2020
  * [67] Yue M Lu and Gen Li  “Phase transitions of spectral initialization for high-dimensional non-convex estimation”  In _Information and Inference: A Journal of the IMA_ Oxford University Press, 2019
  * [68] Zongming Ma  “Sparse principal component analysis and iterative thresholding”  In _The Annals of Statistics_ 41.2, 2013, pp. 772–801  eprint: 1112.2432
  * [69] Vladimir Alexandrovich Marchenko and Leonid Andreevich Pastur  “Distribution of eigenvalues for some sets of random matrices”  In _Matematicheskii Sbornik_ 114.4 Russian Academy of Sciences, Steklov Mathematical Institute of Russian …, 1967, pp. 507–536
  * [70] C.. Martin and M.. Mahoney  “Implicit Self-Regularization in Deep Neural Networks: Evidence from Random Matrix Theory and Implications for Learning”, 2018
  * [71] C.. Martin, T.. Peng and M.. Mahoney  “Predicting trends in the quality of state-of-the-art neural networks without access to training or testing data” Accepted for publication, _Nature Communications_., 2020
  * [72] Llew Mason, Jonathan Baxter, Peter L Bartlett and Marcus R Frean  “Boosting algorithms as gradient descent”  In _Advances in neural information processing systems_ , 2000, pp. 512–518
  * [73] Marco Mondelli and Andrea Montanari  “Fundamental Limits of Weak Recovery with Applications to Phase Retrieval”  In _Conference On Learning Theory_ , 2018, pp. 1445–1450
  * [74] Praneeth Netrapalli, Prateek Jain and Sujay Sanghavi  “Phase retrieval using alternating minimization”  In _Advances in Neural Information Processing Systems_ , 2013, pp. 2796–2804
  * [75] Vardan Papyan  “Traces of Class/Cross-Class Structure Pervade Deep Learning Spectra”  In _Journal of Machine Learning Research_ 21.252, 2020, pp. 1–64
  * [76] Leonid A Pastur  “A simple approach to the global regime of Gaussian ensembles of random matrices”  In _Ukrainian Mathematical Journal_ 57.6 Springer, 2005, pp. 936–966
  * [77] Debashis Paul  “Asymptotics of sample eigenstructure for a large dimensional spiked covariance model”  In _Statistica Sinica_ JSTOR, 2007, pp. 1617–1642
  * [78] Debashis Paul and Jack W Silverstein  “No eigenvalues outside the support of the limiting empirical spectral distribution of a separable covariance matrix”  In _Journal of Multivariate Analysis_ 100.1 Elsevier, 2009, pp. 37–57
  * [79] Fabian Pedregosa and Damien Scieur  “Average-case Acceleration Through Spectral Density Estimation”  In _arXiv preprint arXiv:2002.04756_ , 2020
  * [80] Jeffrey Pennington and Yasaman Bahri  “Geometry of neural network loss surfaces via random matrix theory”  In _Proceedings of the 34th International Conference on Machine Learning-Volume 70_ , 2017, pp. 2798–2806  JMLR. org
  * [81] Jeffrey Pennington and Pratik Worah  “The spectrum of the fisher information matrix of a single-hidden-layer neural network”  In _Advances in Neural Information Processing Systems_ 31, 2018, pp. 5410–5419
  * [82] Farbod Roosta-Khorasani and Michael W. Mahoney  “Sub-sampled Newton methods”  In _Mathematical Programming_ 174.1-2, 2019, pp. 293–326  DOI: 10.1007/s10107-018-1346-5
  * [83] Levent Sagun, Leon Bottou and Yann LeCun  “Eigenvalues of the hessian in deep learning: Singularity and beyond”  In _arXiv preprint arXiv:1611.07476_ , 2016
  * [84] Levent Sagun et al.  “Empirical analysis of the hessian of over-parametrized neural networks”  In _arXiv preprint arXiv:1706.04454_ , 2017
  * [85] Adepu Ravi Sankar, Yash Khasbage, Rahul Vigneswaran and Vineeth N Balasubramanian  “A Deeper Look at the Hessian Eigenspectrum of Deep Neural Networks and its Applications to Regularization”  In _arXiv preprint arXiv:2012.03801_ , 2020
  * [86] Sheng Shen et al.  “Q-BERT: Hessian Based Ultra Low Precision Quantization of BERT”  In _AAAI_ , 2020, pp. 8815–8821
  * [87] Jack W Silverstein and ZD Bai  “On the empirical distribution of eigenvalues of a class of large dimensional random matrices”  In _Journal of Multivariate analysis_ 54.2 Elsevier, 1995, pp. 175–192
  * [88] Terence Tao, Van Vu and Manjunath Krishnapur  “Random matrices: Universality of ESDs and the circular law”  In _The Annals of Probability_ 38.5 Institute of Mathematical Statistics, 2010, pp. 2023–2065
  * [89] Aad W. der Vaart  “Asymptotic Statistics” 3, Cambridge Series in Statistical and Probabilistic Mathematics  Cambridge University Press, 2000  DOI: 10.1017/cbo9780511802256
  * [90] Ulrike Von Luxburg  “A tutorial on spectral clustering”  In _Statistics and computing_ 17.4 Springer, 2007, pp. 395–416
  * [91] Martin J. Wainwright  “High-Dimensional Statistics: : A Non-Asymptotic Viewpoint”, Cambridge Series in Statistical and Probabilistic Mathematics  Cambridge University Press, 2019  DOI: 10.1017/9781108627771
  * [92] Shusen Wang, Farbod Roosta-Khorasani, Peng Xu and Michael W. Mahoney  “GIANT: Globally Improved Approximate Newton Method for Distributed Optimization”  In _Advances in Neural Information Processing Systems_ 31 Curran Associates, Inc., 2018, pp. 2332–2342
  * [93] Eugene P. Wigner  “Characteristic Vectors of Bordered Matrices With Infinite Dimensions”  In _Annals of Mathematics_ 62.3 Annals of Mathematics, 1955, pp. 548–564
  * [94] Yichao Wu and Yufeng Liu  “Robust truncated hinge loss support vector machines”  In _Journal of the American Statistical Association_ 102.479 Taylor & Francis, 2007, pp. 974–983
  * [95] Yikai Wu et al.  “Dissecting Hessian: Understanding Common Structure of Hessian in Neural Networks”  In _arXiv preprint arXiv:2010.04261_ , 2020
  * [96] Peng Xu and Michael W. Mahoney  “Second-order Optimization for Non-convex Machine Learning: an Empirical Study”  In _Proceedings of the 2020 SIAM International Conference on Data Mining_ , 2020, pp. 199–207  DOI: 10.1137/1.9781611976236.23
  * [97] Peng Xu et al.  “Sub-sampled Newton Methods with Non-uniform Sampling”  In _Advances in Neural Information Processing Systems_ 29 Curran Associates, Inc., 2016, pp. 3000–3008
  * [98] Masaki Yamada and Alexander Vilenkin  “Hessian eigenvalue distribution in a random Gaussian landscape”  In _Journal of High Energy Physics_ 2018.3, 2018, pp. 29  DOI: 10.1007/jhep03(2018)029
  * [99] Fan Yang  “Edge universality of separable covariance matrices”  In _Electronic Journal of Probability_ 24 The Institute of Mathematical Statisticsthe Bernoulli Society, 2019
  * [100] Zhewei Yao et al.  “ADAHESSIAN: An Adaptive Second Order Optimizer for Machine Learning”  In _arXiv preprint arXiv:2006.00719_ , 2020
  * [101] Zhewei Yao et al.  “Hessian-based Analysis of Large Batch Training and Robustness to Adversaries”  In _arXiv_ , 2018  eprint: 1802.08241
  * [102] Zhewei Yao, Amir Gholami, Kurt Keutzer and Michael Mahoney  “PyHessian: Neural networks through the lens of the Hessian”  In _arXiv preprint arXiv:1912.07145_ , 2019
  * [103] Zhewei Yao et al.  “Trust Region Based Adversarial Attack on Neural Networks”  In _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ 00, CVPR’19, 2019, pp. 11342–11351  DOI: 10.1109/cvpr.2019.01161
  * [104] Lenka Zdeborová and Florent Krzakala  “Statistical physics of inference: thresholds and algorithms”  In _Advances in Physics_ 65.5, 2016, pp. 453–552  DOI: 10.1080/00018732.2016.1211393
  * [105] Lixin Zhang  “Spectral analysis of large dimensional random matrices”  In _National University of Singapore PHD Thesis_ , 2006

##  Appendix A Proofs

In this section, we present the proofs of our main technical results, Theorem
1-4. We will start with the key technical Theorem 4, which, as discussed in
Section 2.3, will serve as the cornerstone in proving Theorem 1, 2 and 3.

In the sequel, the big $O$ notation $O(u_{n})$ is understood in an almost sure
sense when the objects are random. When multidimensional objects are
concerned, $O(u_{n})$ is understood entry-wise and $O_{\|\cdot\|}(\cdot)$ is
understood, for a vector $\mathbf{v}=O_{\|\cdot\|}(u_{n})$, as its Euclidean
norm satisfying $\|\mathbf{v}\|=O(u_{n})$ and for a square matrix
$\mathbf{M}$, as its spectral norm being of the order $O(u_{n})$. The small
$o$ notation is understood likewise. Note that under Assumption 1 it is
equivalent to using either $O(u_{n})$ or $O(u_{p})$ as $n,p$ scales linearly,
and we will use constantly $O(u_{n})$ for simplicity of exposition.

###  A.1 Proof of Theorem 4

The proof of the deterministic equivalent (recall from Section 2.3) result of
the type $\mathbf{Q}(z)\leftrightarrow{\mathbb{E}}[\mathbf{Q}(z)]$ generally
comes in two steps:

  1. 1.

establish the convergence in mean, that is, to show that the proposed
$\bar{\mathbf{Q}}(z)$ is an asymptotic approximation for the expected
resolvent ${\mathbb{E}}[\mathbf{Q}(z)]$ in the sense that

$\|\bar{\mathbf{Q}}(z)-{\mathbb{E}}[\mathbf{Q}(z)]\|\to 0,$ (16)

as $n,p\to\infty$; and

  2. 2.

show that the random quantities
$\frac{1}{p}\operatorname{tr}\mathbf{A}\mathbf{Q}(z)$ and $\mathbf{a}^{\sf
T}\mathbf{Q}(z)\mathbf{b}$ of interest concentrate around their expectations
in the sense that

$\frac{1}{p}\operatorname{tr}\mathbf{A}(\mathbf{Q}(z)-{\mathbb{E}}\mathbf{Q}(z))\to
0,\quad\text{~{}and~{}}\mathbf{a}^{\sf
T}(\mathbf{Q}(z)-{\mathbb{E}}\mathbf{Q}(z))\mathbf{b}\to 0,$ (17)

almost surely as $n,p\to\infty$. This can be done with concentration of
measure arguments [60] or, in the Gaussian case, the Nash–Poincaré inequality
to bound the variance (see for example [76, Proposition 2.4]) and then
applying Markov’s inequality and the Borel Cantelli lemma.

Here we will focus on the first part of the proof and show for
$\mathbf{Q}(z)=(\mathbf{H}-z\mathbf{I}_{p})^{-1}$,
$\mathbf{H}\in{\mathbb{R}}^{p\times p}$ defined in (4) that

$\|{\mathbb{E}}[\mathbf{Q}(z)]-\bar{\mathbf{Q}}(z)\|\to
0,\quad\bar{\mathbf{Q}}(z)\equiv\left({\mathbb{E}}\left[\frac{g}{1+g\cdot\delta(z)}(\mathbf{C}^{\frac{1}{2}}(\mathbf{I}_{p}-\mathbf{P}_{\mathbf{U}})\mathbf{C}^{\frac{1}{2}}+\boldsymbol{\alpha}\boldsymbol{\alpha}^{\sf
T})\right]-z\mathbf{I}_{p}\right)^{-1},$ (18)

as $n,p\to\infty$, for random vector
$\boldsymbol{\alpha}\equiv\boldsymbol{\mu}+\mathbf{C}^{\frac{1}{2}}\mathbf{P}_{\mathbf{U}}\mathbf{z}\in{\mathbb{R}}^{p}$
and $g=\ell^{\prime\prime}(y,\mathbf{w}^{\sf
T}\boldsymbol{\mu}+\mathbf{w}^{\sf T}\mathbf{C}^{\frac{1}{2}}\mathbf{z})$ for
$\mathbf{z}\sim{\mathcal{N}}(\mathbf{0},\mathbf{I}_{p})$ as in (8), $y$
defined in (1), $\mathbf{P}_{\mathbf{U}}$ the projection onto the subspace
spanned by the columns of
$\mathbf{U}=\mathbf{C}^{\frac{1}{2}}[\mathbf{w}_{*},~{}\mathbf{w}]\in{\mathbb{R}}^{p\times
2}$ so that $\mathbf{P}_{\mathbf{U}}=\mathbf{U}(\mathbf{U}^{\sf
T}\mathbf{U})^{-1}\mathbf{U}^{\sf T}$ for $\mathbf{U}^{\sf T}\mathbf{U}$
invertible, and $\delta(z)$ the unique solution to (7).

Before diving into the detailed proof of Theorem 4, we will illustrate the
main idea leading to Theorem 4, notably the form of $\bar{\mathbf{Q}}$ in
(18), in the simple setting of $\boldsymbol{\mu}=\mathbf{0}$ and
$\mathbf{C}=\mathbf{I}_{p}$. As discussed above in Section 2.1, the main
technical difficulty lies in the fact that the diagonal matrix
$\mathbf{D}=\operatorname{diag}\\{\ell^{\prime\prime}(y_{i},\mathbf{w}^{\sf
T}\mathbf{x}_{i})\\}_{i=1}^{n}$ depends on the $\mathbf{x}_{i}$’s in a rather
involved nonlinear manner per (1). Nonetheless, note that standard Gaussian
vector $\mathbf{x}\sim{\mathcal{N}}(\mathbf{0},\mathbf{I}_{p})$ can decomposed
as

$\mathbf{x}=\mathbf{U}(\mathbf{U}^{\sf T}\mathbf{U})^{-1}\mathbf{U}^{\sf
T}\mathbf{x}+\mathbf{x}^{\perp}\equiv\mathbf{P}_{\mathbf{U}}\mathbf{x}+\mathbf{x}^{\perp},\quad\text{for}\quad\mathbf{x}^{\perp}\equiv(\mathbf{I}_{p}-\mathbf{P}_{\mathbf{U}})\mathbf{x},$
(19)

the sum of its projection onto the subspace spanned by the columns of
$\mathbf{U}=[\mathbf{w}_{*},~{}\mathbf{w}]$ (recall that
$\mathbf{C}=\mathbf{I}_{p}$ and assume for the moment that
$\mathbf{w},\mathbf{w}_{*}$ are linearly independent) and the complementary
subspace. Due to the orthogonal invariance of standard multivariate Gaussian,
we have that both $\mathbf{P}_{\mathbf{U}}\mathbf{x}$ and
$\mathbf{x}^{\perp}=(\mathbf{I}_{p}-\mathbf{P}_{\mathbf{U}})\mathbf{x}$ are
Gaussian, and they are _independent_ , for example by checking that
${\mathbb{E}}[\mathbf{P}_{\mathbf{U}}\mathbf{x}\mathbf{x}^{\sf
T}(\mathbf{I}_{p}-\mathbf{P}_{\mathbf{U}})]=\mathbf{0}$.

Since $\mathbf{D}$ depends on $\mathbf{X}$ only through $\mathbf{U}^{\sf
T}\mathbf{X}$, we have that $\mathbf{D}$ is _independent_ of
$\mathbf{X}^{\perp}\equiv[\mathbf{x}_{1}^{\perp},\ldots,\mathbf{x}_{n}^{\perp}]=(\mathbf{I}_{p}-\mathbf{P}_{\mathbf{U}})\mathbf{X}\in{\mathbb{R}}^{p\times
n}$ so that

$\mathbf{H}=\frac{1}{n}\mathbf{X}\mathbf{D}\mathbf{X}^{\sf
T}=\frac{1}{n}(\mathbf{P}_{\mathbf{U}}\mathbf{X}+\mathbf{X}^{\perp})\mathbf{D}(\mathbf{P}_{\mathbf{U}}\mathbf{X}+\mathbf{X}^{\perp})^{\sf
T}=\frac{1}{n}\mathbf{X}^{\perp}\mathbf{D}(\mathbf{X}^{\perp})^{\sf
T}+\mathbf{A}$ (20)

for some $\mathbf{A}\in{\mathbb{R}}^{p\times p}$ of low rank (since
$\operatorname{rank}(\mathbf{P}_{\mathbf{U}})\leq 2$). As a consequence, our
object of interest $\mathbf{H}$ is shown to follow a spiked random matrix
model as discussed in Section 2.2, with the main bulk
$\frac{1}{n}\mathbf{X}^{\perp}\mathbf{D}(\mathbf{X}^{\perp})^{\sf T}$ now more
accessible due to the independence between $\mathbf{X}^{\perp}$ and
$\mathbf{D}$. Moreover, since $\mathbf{X}^{\perp}$ is also a low rank update
of the original $\mathbf{X}$, the above derivation shows that, asymptotically,
the limiting spectral measure of
$\mathbf{H}=\frac{1}{n}\mathbf{X}\mathbf{D}\mathbf{X}^{\sf T}$ in (4) is the
same as that of
$\frac{1}{n}\mathbf{X}^{\perp}\mathbf{D}(\mathbf{X}^{\perp})^{\sf T}$ and thus
the same as if $\mathbf{D}$ was _independent_ of $\mathbf{X}$. While the
analysis of limiting eigenvalue distribution of $\mathbf{H}$ is largely
simplified based on the argument above, the characterization of the low rank
matrix $\mathbf{A}$ requires more effort, so as to determine the behavior of
the possible isolated eigenvalue-eigenvector pairs of $\mathbf{H}$.

The idea above can be generalized, under Assumption 1, to arbitrary profile of
the statistical mean $\boldsymbol{\mu}$ and covariance $\mathbf{C}$, leading
to Theorem 4, the proof of which is given as follows.

###### Proof of Theorem 4.

For mathematical convenience, for the resolvent
$\mathbf{Q}(z)=(\mathbf{H}-z\mathbf{I}_{p})^{-1}$, we will take $z<0$ outside
the spectral support of $\mathbf{H}$ in what follows. Since $\mathbf{Q}(z)$
and $\bar{\mathbf{Q}}(z)$ from the theorem statement are complex analytic
functions for $z\notin{\mathbb{R}}^{+}$ (matrix-valued Stieltjes transforms
are analytic [43]), obtaining the convergence results on ${\mathbb{R}}^{-}$ is
equivalent to obtaining the result on all of
${\mathbb{C}}\setminus{\mathbb{R}}^{+}$.

Following the orthogonal decomposition idea in (19), for
$\mathbf{U}=\mathbf{C}^{\frac{1}{2}}[\mathbf{w}_{*},~{}\mathbf{w}]\in{\mathbb{R}}^{p\times
2}$ and $\mathbf{x}_{i}\sim{\mathcal{N}}(\boldsymbol{\mu},\mathbf{C})$ we have
$\mathbf{z}_{i}=\mathbf{C}^{-\frac{1}{2}}(\mathbf{x}_{i}-\boldsymbol{\mu})\sim{\mathcal{N}}(\mathbf{0},\mathbf{I}_{p})$
and

$\mathbf{z}_{i}=\mathbf{P}_{\mathbf{U}}\mathbf{z}_{i}+\mathbf{z}_{i}^{\perp},\quad\text{~{}for~{}}\mathbf{z}_{i}^{\perp}\equiv(\mathbf{I}_{p}-\mathbf{P}_{\mathbf{U}})\mathbf{z}_{i}.$
(21)

Note that (i) $\mathbf{U}^{\sf T}\mathbf{z}_{i}^{\perp}=\mathbf{0}$, (ii)
$\mathbf{U}^{\sf T}\mathbf{z}_{i}$ is independent of $\mathbf{z}_{i}^{\perp}$
due to the Gaussianity of $\mathbf{z}_{i}$ and

${\mathbb{E}}[\mathbf{z}_{i}^{\perp}]=\mathbf{0},\quad{\mathbb{E}}[\mathbf{z}_{i}^{\perp}(\mathbf{z}_{i}^{\perp})^{\sf
T}]=\mathbf{I}_{p}-\mathbf{P}_{\mathbf{U}}.$ (22)

With these results, we now move on to the proof of (18). Recall the definition
$\mathbf{Q}=(\frac{1}{n}\mathbf{X}\mathbf{D}\mathbf{X}^{\sf
T}-z\mathbf{I}_{p})^{-1}$ with
$\mathbf{D}=\operatorname{diag}\\{g_{i}\\}_{i=1}^{n}$ and

$g_{i}\equiv\ell^{\prime\prime}(y_{i},\mathbf{w}^{\sf T}\mathbf{x}_{i}),$ (23)

we first introduce the following intermediate deterministic quantities

$\theta\equiv\frac{1}{n}\operatorname{tr}\mathbf{C}{\mathbb{E}}[\mathbf{Q}_{-1}]>0,\quad\mathbf{M}_{\theta}\equiv{\mathbb{E}}\left[\frac{g_{1}\cdot\mathbf{x}_{1}\mathbf{x}_{1}}{1+g_{1}\cdot\theta}\right],\quad\bar{\mathbf{Q}}_{\theta}\equiv\left(\mathbf{M}_{\theta}-z\mathbf{I}_{p}\right)^{-1},$
(24)

where we denote, for i.i.d. $\mathbf{x}_{i},g_{i}$ and without loss of
generality, $\mathbf{Q}_{-1}\equiv(\frac{1}{n}\sum_{i\neq
1}g_{i}\mathbf{x}_{i}\mathbf{x}_{i}^{\sf
T}-z\mathbf{I}_{p})^{-1}\in{\mathbb{R}}^{p\times p}$ the “leave-one-out”
version of $\mathbf{Q}$ that is _independent_ of $\mathbf{x}_{1}$ and $g_{1}$.
Our proof strategy is to show that

$\|{\mathbb{E}}[\mathbf{Q}]-\bar{\mathbf{Q}}_{\theta}\|\to
0,\quad\|\bar{\mathbf{Q}}_{\theta}-\bar{\mathbf{Q}}\|\to 0,$

and thus the conclusion of (18).

With the resolvent identity
$\mathbf{A}^{-1}-\mathbf{B}^{-1}=\mathbf{A}^{-1}(\mathbf{B}-\mathbf{A})\mathbf{B}^{-1}$
for invertible $\mathbf{A},\mathbf{B}$, we write,

$\displaystyle{\mathbb{E}}[\mathbf{Q}-\bar{\mathbf{Q}}_{\theta}]$
$\displaystyle={\mathbb{E}}\left[\mathbf{Q}\left(\mathbf{M}_{\theta}-\frac{1}{n}\mathbf{X}\mathbf{D}\mathbf{X}^{\sf
T}\right)\right]\bar{\mathbf{Q}}_{\theta}={\mathbb{E}}[\mathbf{Q}]\mathbf{M}_{\theta}\bar{\mathbf{Q}}_{\theta}-\frac{1}{n}\sum_{i=1}^{n}{\mathbb{E}}[g_{i}\mathbf{Q}\mathbf{x}_{i}\mathbf{x}_{i}^{\sf
T}]\bar{\mathbf{Q}}_{\theta}$
$\displaystyle={\mathbb{E}}[\mathbf{Q}]\mathbf{M}_{\theta}\bar{\mathbf{Q}}-\frac{1}{n}\sum_{i=1}^{n}{\mathbb{E}}\left[\frac{g_{i}\cdot\mathbf{Q}_{-i}\mathbf{x}_{i}\mathbf{x}_{i}^{\sf
T}}{1+g_{i}\cdot\frac{1}{n}\mathbf{x}_{i}^{\sf
T}\mathbf{Q}_{-i}\mathbf{x}_{i}}\right]\bar{\mathbf{Q}}_{\theta},$

where we used the Sherman–Morrison formula.

Since the random variable $\frac{1}{n}\mathbf{x}_{i}^{\sf
T}\mathbf{Q}_{-i}\mathbf{x}_{i}$ is expected to concentrate around its
expectation
$\theta=\frac{1}{n}\operatorname{tr}\mathbf{C}{\mathbb{E}}[\mathbf{Q}_{-i}]$,
we write

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}{\mathbb{E}}\left[\frac{g_{i}\cdot\mathbf{Q}_{-i}\mathbf{x}_{i}\mathbf{x}_{i}^{\sf
T}}{1+g_{i}\cdot\frac{1}{n}\mathbf{x}_{i}^{\sf
T}\mathbf{Q}_{-i}\mathbf{x}_{i}}\right]$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}{\mathbb{E}}\left[\frac{g_{i}\cdot\mathbf{Q}_{-i}\mathbf{x}_{i}\mathbf{x}_{i}^{\sf
T}}{1+g_{i}\cdot\theta}\right]+\frac{1}{n}\sum_{i=1}^{n}{\mathbb{E}}\left[\frac{g_{i}^{2}\cdot\mathbf{Q}_{-i}\mathbf{x}_{i}\mathbf{x}_{i}^{\sf
T}(\theta-\frac{1}{n}\mathbf{x}_{i}^{\sf
T}\mathbf{Q}_{-i}\mathbf{x}_{i})}{(1+g_{i}\cdot\frac{1}{n}\mathbf{x}_{i}^{\sf
T}\mathbf{Q}_{-i}\mathbf{x}_{i})(1+g_{i}\cdot\theta)}\right]$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}{\mathbb{E}}[\mathbf{Q}_{-i}]\mathbf{M}_{\theta}+\frac{1}{n}\sum_{i=1}^{n}{\mathbb{E}}\left[\frac{g_{i}^{2}\cdot\mathbf{Q}\mathbf{x}_{i}\mathbf{x}_{i}^{\sf
T}(\theta-\frac{1}{n}\mathbf{x}_{i}^{\sf
T}\mathbf{Q}_{-i}\mathbf{x}_{i})}{1+g_{i}\cdot\theta}\right]$
$\displaystyle\equiv{\mathbb{E}}[\mathbf{Q}_{-1}]\mathbf{M}_{\theta}+\frac{1}{n}{\mathbb{E}}[\mathbf{Q}\mathbf{X}\mathbf{D}\mathbf{D}_{1}\mathbf{X}^{\sf
T}]$

for
$\mathbf{D}_{1}\equiv\operatorname{diag}\\{\frac{g_{i}}{1+g_{i}\cdot\theta}(\theta-\frac{1}{n}\mathbf{x}_{i}^{\sf
T}\mathbf{Q}_{-i}\mathbf{x}_{i})\\}_{i=1}^{n}$ and
$\mathbf{D}=\operatorname{diag}\\{g_{i}\\}_{i=1}^{n}$, where we used the fact
that $\mathbf{Q}_{-i}$ is independent of both $g_{i}$ and $\mathbf{x}_{i}$. As
such,

${\mathbb{E}}[\mathbf{Q}-\bar{\mathbf{Q}}_{\theta}]={\mathbb{E}}[\mathbf{Q}-\mathbf{Q}_{-1}]\mathbf{M}_{\theta}\bar{\mathbf{Q}}_{\theta}-\frac{1}{n}{\mathbb{E}}[\mathbf{Q}\mathbf{X}\mathbf{D}\mathbf{D}_{1}\mathbf{X}^{\sf
T}]\bar{\mathbf{Q}}_{\theta}.$ (25)

For the first right-hand side term in (25), again by the resolvent identity

$\displaystyle{\mathbb{E}}[\mathbf{Q}_{-1}-\mathbf{Q}]$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}{\mathbb{E}}\left[\mathbf{Q}_{-i}-\mathbf{Q}\right]=\frac{1}{n}\sum_{i=1}^{n}{\mathbb{E}}\left[\mathbf{Q}_{-i}\frac{1}{n}g_{i}\mathbf{x}_{i}\mathbf{x}_{i}^{\sf
T}\mathbf{Q}\right]=\frac{1}{n}\sum_{i=1}^{n}{\mathbb{E}}\left[\frac{1}{n}\frac{g_{i}\cdot\mathbf{Q}_{-i}\mathbf{x}_{i}\mathbf{x}_{i}^{\sf
T}\mathbf{Q}_{-i}}{1+g_{i}\cdot\frac{1}{n}\mathbf{x}_{i}^{\sf
T}\mathbf{Q}_{-i}\mathbf{x}_{i}}\right]$
$\displaystyle=\frac{1}{n}{\mathbb{E}}\left[\mathbf{Q}_{-1}\frac{g_{1}\mathbf{x}_{1}\mathbf{x}_{1}^{\sf
T}}{1+g_{i}\cdot\theta}\mathbf{Q}_{-1}\right]+\frac{1}{n}{\mathbb{E}}\left[\mathbf{Q}_{-1}\frac{g_{1}^{2}\mathbf{x}_{1}\mathbf{x}_{1}^{\sf
T}(\theta-\frac{1}{n}\mathbf{x}_{1}^{\sf
T}\mathbf{Q}_{-1}\mathbf{x}_{1})}{(1+g_{1}\cdot\theta)(1+g_{1}\cdot\frac{1}{n}\mathbf{x}_{1}^{\sf
T}\mathbf{Q}_{-1}\mathbf{x}_{1})}\mathbf{Q}_{-1}\right]$
$\displaystyle=\frac{1}{n}{\mathbb{E}}[\mathbf{Q}_{-1}\mathbf{M}_{\theta}\mathbf{Q}_{-1}]+\frac{1}{n}{\mathbb{E}}\left[\mathbf{Q}\frac{1}{n}\mathbf{X}\mathbf{D}\mathbf{D}_{1}\mathbf{X}^{\sf
T}\mathbf{Q}_{-1}\right]$

where we used again the Sherman–Morrison formula. The first term is easily
seen to be of order $O(n^{-1})$ in spectral norm (with
$\|\mathbf{Q}(z)\|\leq\operatorname{dist}(z,\operatorname{supp}\mu_{\mathbf{H}})^{-1}$
and $\|\mathbf{M}_{\theta}\|=O(1)$ under Assumption 1). For the second term,
note that with the additional $1/n$ in front, it suffices to bound the
spectral norm of the second right-hand side term in (25). Since
$\|\mathbf{Q}\frac{1}{n}\mathbf{X}\mathbf{D}\mathbf{X}^{\sf T}\|\leq 1$, it
thus remains to control

${\mathbb{E}}\left[\max_{i}\left|\theta-\frac{1}{n}\mathbf{x}_{i}^{\sf
T}\mathbf{Q}_{-i}\mathbf{x}_{i}\right|\right],$ (26)

that appears in $\mathbf{D}_{1}$. To this end, note that

$\displaystyle{\mathbb{E}}\left[\left|\theta-\frac{1}{n}\mathbf{x}_{1}^{\sf
T}\mathbf{Q}_{-1}\mathbf{x}_{1}\right|^{4}\right]\leq{\mathbb{E}}\left[\left|\frac{1}{n}\mathbf{x}_{1}^{\sf
T}\mathbf{Q}_{-1}\mathbf{x}_{1}-\frac{1}{n}\operatorname{tr}\mathbf{C}\mathbf{Q}_{-1}+\frac{1}{n}\operatorname{tr}\mathbf{C}(\mathbf{Q}_{-1}-{\mathbb{E}}[\mathbf{Q}_{-1}])\right|^{4}\right]$
$\displaystyle={\mathbb{E}}\left[\left|\frac{1}{n}\mathbf{x}_{1}^{\sf
T}\mathbf{Q}_{-1}\mathbf{x}_{1}-\frac{1}{n}\operatorname{tr}\mathbf{C}\mathbf{Q}_{-1}\right|^{4}\right]+{\mathbb{E}}\left[\left|\frac{1}{n}\operatorname{tr}\mathbf{C}(\mathbf{Q}_{-1}-{\mathbb{E}}[\mathbf{Q}_{-1}])\right|^{4}\right]$
$\displaystyle+6\cdot{\mathbb{E}}\left[\left|\frac{1}{n}\mathbf{x}_{1}^{\sf
T}\mathbf{Q}_{-1}\mathbf{x}_{1}-\frac{1}{n}\operatorname{tr}\mathbf{C}\mathbf{Q}_{-1}\right|^{2}\cdot\left|\frac{1}{n}\operatorname{tr}\mathbf{C}(\mathbf{Q}_{-1}-{\mathbb{E}}[\mathbf{Q}_{-1}])\right|^{2}\right]+O(n^{-4})=O(n^{-2}),$

with
$\alpha\equiv\frac{1}{n}\operatorname{tr}\mathbf{C}{\mathbb{E}}[\mathbf{Q}_{-1}]$
by definition, where we used, for the last line, [5, Lemma B.26], the
Burkholder inequality (see e.g., [5, Lemma 2.13]), as well as the fact that
$\boldsymbol{\mu}^{\sf
T}\mathbf{Q}_{-i}\boldsymbol{\mu}\leq\|\boldsymbol{\mu}\|^{2}\|\mathbf{Q}_{-i}\|=O(1)$
under Assumption 1.

As a result, by Markov’s inequality
$\mathbb{P}(X>t)\leq{\mathbb{E}}[X^{k}]/t^{k}$ for every $k$ (for $X,t>0$) and
the union bound, we have that

$\displaystyle{\mathbb{E}}\left[\max_{i}\left|\theta-\frac{1}{n}\mathbf{x}_{i}^{\sf
T}\mathbf{Q}_{-i}\mathbf{x}_{i}\right|\right]$
$\displaystyle=\int_{0}^{C^{\prime}n^{-1/4}}\mathbb{P}\left(\max_{i}\left|\theta-\frac{1}{n}\mathbf{x}_{i}^{\sf
T}\mathbf{Q}_{-i}\mathbf{x}_{i}\right|>t\right)\,dt$
$\displaystyle+n\int_{C^{\prime}n^{-1/4}}^{\infty}\mathbb{P}\left(\left|\theta-\frac{1}{n}\mathbf{x}_{1}^{\sf
T}\mathbf{Q}_{-1}\mathbf{x}_{1}\right|>t\right)\,dt\leq Cn^{-\frac{1}{4}}$

holds for some constant $C,C^{\prime}>0$ that depend on $p,n$ only via the
ratio $p/n$.

Putting all together, we obtain the first approximation

$\|{\mathbb{E}}[\mathbf{Q}]-\bar{\mathbf{Q}}_{\theta}\|=O(n^{-\frac{1}{4}}),\quad\text{~{}for~{}}\bar{\mathbf{Q}}_{\theta}\equiv\left(\mathbf{M}_{\theta}-z\mathbf{I}_{p}\right)^{-1},\quad\mathbf{M}_{\theta}\equiv{\mathbb{E}}\left[\frac{g_{1}\cdot\mathbf{x}_{1}\mathbf{x}_{1}}{1+g_{1}\cdot\theta}\right],$
(27)

with
$\theta\equiv\frac{1}{n}\operatorname{tr}\mathbf{C}{\mathbb{E}}[\mathbf{Q}_{-1}]>0$.

Note in passing we also proved that
$\|{\mathbb{E}}[\mathbf{Q}-\mathbf{Q}_{-1}]\|=O(n^{-1})$.

It thus remains to establish the second approximation by showing

$\|\bar{\mathbf{Q}}_{\theta}-\bar{\mathbf{Q}}\|\to 0,$ (28)

as $n,p\to\infty$, with $\bar{\mathbf{Q}}=(\mathbf{M}-z\mathbf{I}_{p})^{-1}$
and
$\mathbf{M}={\mathbb{E}}\left[\frac{g}{1+g\cdot\delta(z)}(\mathbf{C}^{\frac{1}{2}}(\mathbf{I}_{p}-\mathbf{P}_{\mathbf{U}})\mathbf{C}^{\frac{1}{2}}+\boldsymbol{\alpha}\boldsymbol{\alpha}^{\sf
T})\right]$ as defined in (18). To this end, we first compute the expectation
in $\mathbf{M}_{\theta}$. Note that
$g_{i}=\ell^{\prime\prime}(y_{i},\mathbf{w}^{\sf T}\mathbf{x}_{i})$ depends on
$\mathbf{x}_{i}$ _only_ via the projection $\mathbf{U}^{\sf T}\mathbf{z}_{i}$
under (1). We thus decompose $\mathbf{x}_{i}\mathbf{x}_{i}^{\sf T}$ as

$\displaystyle\mathbf{x}_{i}\mathbf{x}_{i}^{\sf T}$
$\displaystyle=(\boldsymbol{\mu}+\mathbf{C}^{\frac{1}{2}}\mathbf{z}_{i})(\boldsymbol{\mu}+\mathbf{C}^{\frac{1}{2}}\mathbf{z}_{i})^{\sf
T}=(\boldsymbol{\mu}+\mathbf{C}^{\frac{1}{2}}\mathbf{P}_{\mathbf{U}}\mathbf{z}_{i}+\mathbf{C}^{\frac{1}{2}}\mathbf{z}_{i}^{\perp})(\boldsymbol{\mu}+\mathbf{C}^{\frac{1}{2}}\mathbf{P}_{\mathbf{U}}\mathbf{z}_{i}+\mathbf{C}^{\frac{1}{2}}\mathbf{z}_{i}^{\perp})^{\sf
T}$ $\displaystyle=\boldsymbol{\alpha}_{i}\boldsymbol{\alpha}_{i}^{\sf
T}+\mathbf{C}^{\frac{1}{2}}\mathbf{z}_{i}^{\perp}(\mathbf{z}_{i}^{\perp})^{\sf
T}\mathbf{C}^{\frac{1}{2}}+\boldsymbol{\alpha}_{i}(\mathbf{z}_{i}^{\perp})^{\sf
T}\mathbf{C}^{\frac{1}{2}}+\mathbf{C}^{\frac{1}{2}}\mathbf{z}_{i}^{\perp}\boldsymbol{\alpha}_{i}^{\sf
T},$

where we denote
$\boldsymbol{\alpha}_{i}\equiv\boldsymbol{\mu}+\mathbf{C}^{\frac{1}{2}}\mathbf{P}_{\mathbf{U}}\mathbf{z}_{i}$.
As such

$\displaystyle\mathbf{M}_{\theta}$
$\displaystyle\equiv{\mathbb{E}}\left[\frac{g_{1}\cdot\mathbf{x}_{1}\mathbf{x}_{1}}{1+g_{1}\cdot\theta}\right]={\mathbb{E}}\left[\frac{g_{1}\cdot(\boldsymbol{\alpha}_{1}\boldsymbol{\alpha}_{1}^{\sf
T}+\mathbf{C}^{\frac{1}{2}}\mathbf{z}_{i}^{\perp}(\mathbf{z}_{i}^{\perp})^{\sf
T}\mathbf{C}^{\frac{1}{2}})}{1+g_{1}\cdot\theta}\right]$
$\displaystyle={\mathbb{E}}\left[\frac{g_{1}}{1+g_{1}\cdot\theta}\left(\boldsymbol{\alpha}_{1}\boldsymbol{\alpha}_{1}^{\sf
T}+\mathbf{C}^{\frac{1}{2}}(\mathbf{I}_{p}-\mathbf{P}_{\mathbf{U}})\mathbf{C}^{\frac{1}{2}}\right)\right]$

where we used the crucial fact that $\mathbf{z}_{i}^{\perp}$ is independent of
$g_{i},\mathbf{P}_{\mathbf{U}}\mathbf{z}_{i}$ and $\boldsymbol{\alpha}_{i}$
due to the Gaussianity of $\mathbf{z}_{i}$, with
${\mathbb{E}}[\mathbf{z}_{i}^{\perp}]=\mathbf{0}$ and
${\mathbb{E}}[\mathbf{z}_{i}^{\perp}(\mathbf{z}_{i}^{\perp})^{\sf
T}]=\mathbf{I}_{p}-\mathbf{P}_{\mathbf{U}}$ from (22). As a result, we have

$\mathbf{M}-\mathbf{M}_{\theta}={\mathbb{E}}\left[\frac{g^{2}(\delta-\theta)}{(1+g\cdot\delta)(1+g\cdot\theta)}\cdot(\mathbf{C}^{\frac{1}{2}}(\mathbf{I}_{p}-\mathbf{P}_{\mathbf{U}})\mathbf{C}^{\frac{1}{2}}+\boldsymbol{\alpha}\boldsymbol{\alpha}^{\sf
T})\right]$ (29)

for $\delta=\delta(z)=\frac{1}{p}\operatorname{tr}\bar{\mathbf{Q}}_{b}(z)$ as
defined in (7).

Further note that

$\bar{\mathbf{Q}}_{\theta}-\bar{\mathbf{Q}}=\bar{\mathbf{Q}}_{\theta}\left(\mathbf{M}-\mathbf{M}_{\theta}\right)\bar{\mathbf{Q}}$

so that with
$\bar{\mathbf{Q}}_{b}(z)=\left({\mathbb{E}}\left[\frac{g\cdot\mathbf{C}}{1+g\delta(z)}\right]-z\mathbf{I}-p\right)^{-1}$
as defined in (7),

$\displaystyle|\delta-\theta|=\left|\frac{1}{n}\operatorname{tr}\mathbf{C}\left(\bar{\mathbf{Q}}_{b}-{\mathbb{E}}[\mathbf{Q}_{-1}]\right)\right|$
$\displaystyle\leq\left|\frac{1}{n}\operatorname{tr}\mathbf{C}(\bar{\mathbf{Q}}_{b}-\bar{\mathbf{Q}})\right|+\left|\frac{1}{n}\operatorname{tr}\mathbf{C}(\bar{\mathbf{Q}}-\bar{\mathbf{Q}}_{\theta})\right|$
$\displaystyle+\left|\frac{1}{n}\operatorname{tr}\mathbf{C}(\bar{\mathbf{Q}}_{\theta}-{\mathbb{E}}[\mathbf{Q}])\right|+\left|\frac{1}{n}\operatorname{tr}\mathbf{C}{\mathbb{E}}[\mathbf{Q}-\mathbf{Q}_{-1}]\right|$

where we used the fact that
$\|{\mathbb{E}}[\mathbf{Q}]-\bar{\mathbf{Q}}_{\theta}\|=O(n^{-\frac{1}{4}})$
and $\|{\mathbb{E}}[\mathbf{Q}-\mathbf{Q}_{-1}]\|=O(n^{-1})$ to bound the
third and fourth right-hand side term, respectively, and the following lemma
to bound the first right-hand side term.

######  Lemma 1 ([87, Lemma 2.6]).

For $\mathbf{A},\mathbf{M}\in{\mathbb{R}}^{n\times n}$ symmetric,
$\mathbf{u}\in{\mathbb{R}}^{n}$, $\tau\in{\mathbb{R}}$ and
$z\in{\mathbb{C}}\setminus{\mathbb{R}}$,

$\operatorname{tr}\mathbf{A}(\mathbf{M}+\tau\mathbf{u}\mathbf{u}^{\sf
T}-z\mathbf{I}_{n})^{-1}-\operatorname{tr}\mathbf{A}(\mathbf{M}-z\mathbf{I}_{n})^{-1}\leq\frac{\|\mathbf{A}\|}{|\Im[z]|}.$

The inequality further extends to the case of symmetric and positive
semidefinite $\mathbf{A}$, with $\tau>0$ and $z<0$ so long that the inverses
are well-defined, by replacing the right-hand-side $|\Im[z]|$ by $|z|$.

We thus obtain from (29) that

$\displaystyle|\delta-\theta|$
$\displaystyle\leq\left|\frac{1}{n}\operatorname{tr}\mathbf{C}(\bar{\mathbf{Q}}-\bar{\mathbf{Q}}_{\theta})\right|+O(n^{-\frac{1}{4}})$
$\displaystyle=|\delta-\theta|\cdot\frac{1}{n}\operatorname{tr}\left(\mathbf{C}\bar{\mathbf{Q}}_{\theta}{\mathbb{E}}\left[\frac{g^{2}}{(1+g\cdot\delta)(1+g\cdot\delta)}(\mathbf{C}^{\frac{1}{2}}(\mathbf{I}_{p}-\mathbf{P}_{\mathbf{U}})\mathbf{C}^{\frac{1}{2}}+\boldsymbol{\alpha}\boldsymbol{\alpha}^{\sf
T})\bar{\mathbf{Q}}\right]\right)+O(n^{-\frac{1}{4}})$

so that

$|\delta-\theta|=O(n^{-\frac{1}{4}}),\quad\text{and}\quad\|\bar{\mathbf{Q}}_{\theta}-\bar{\mathbf{Q}}\|=O(n^{-\frac{1}{4}}).$
(30)

This concludes the proof of Theorem 4. ∎

Note in passing that we can alternatively write

$\bar{\mathbf{Q}}=\left({\mathbb{E}}\left[\frac{g\cdot\mathbf{C}}{1+g\cdot\delta(z)}\right]-z\mathbf{I}_{p}+\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{\sf
T}\right)^{-1}$ (31)

with
$\mathbf{V}\equiv[\boldsymbol{\mu},~{}\mathbf{C}^{\frac{1}{2}}\mathbf{U}]\in{\mathbb{R}}^{p\times
3}$, $\mathbf{U}=\mathbf{C}^{\frac{1}{2}}[\mathbf{w}_{*},~{}\mathbf{w}]$ and

$\boldsymbol{\Lambda}(z)\equiv{\mathbb{E}}\left[\frac{g}{1+g\cdot\delta(z)}\begin{bmatrix}1&(\mathbf{U}^{+}\mathbf{z})^{\sf
T}\\\
\mathbf{U}^{+}\mathbf{z}&\mathbf{U}^{+}\mathbf{z}(\mathbf{U}^{+}\mathbf{z})^{\sf
T}-(\mathbf{U}^{\sf T}\mathbf{U})^{+}\end{bmatrix}\right]$ (32)

for $\mathbf{U}^{+}$ the Moore–Penrose pseudoinverse of $\mathbf{U}$ as
defined in Theorem 2.

###  A.2 Proof of Theorem 1

The proof of Theorem 1 follows straightforwardly from Theorem 4 by writing the
corresponding Stieltjes transform

$\displaystyle\frac{1}{p}\operatorname{tr}\mathbf{Q}(z)$
$\displaystyle=\frac{1}{p}\operatorname{tr}\left(\mathbf{Q}(z)-{\mathbb{E}}[\mathbf{Q}(z)]+{\mathbb{E}}[\mathbf{Q}(z)]-\bar{\mathbf{Q}}(z)+\bar{\mathbf{Q}}(z)-\bar{\mathbf{Q}}_{b}(z)+\bar{\mathbf{Q}}_{b}(z)\right)$
$\displaystyle=O(n^{-\frac{1}{4}})+O(n^{-1})+\frac{1}{p}\operatorname{tr}\bar{\mathbf{Q}}_{b}(z)$

where we used the fact that
$\|{\mathbb{E}}[\mathbf{Q}]-\bar{\mathbf{Q}}\|=O(n^{-\frac{1}{4}})$ from the
proof of Theorem 4 above and Lemma 1 for the approximations, respectively. ∎

###  A.3 Proof of Theorem 2

For $\bar{\mathbf{Q}}(z)$ defined in Theorem 4, we have, with the Woodbury
identity that

$\bar{\mathbf{Q}}(z)=\left({\mathbb{E}}\left[\frac{g\cdot\mathbf{C}}{1+g\cdot\delta(z)}\right]-z\mathbf{I}_{p}+\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{\sf
T}\right)^{-1}=\left(\bar{\mathbf{Q}}^{-1}_{b}(z)+\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{\sf
T}\right)^{-1}$

with
$\bar{\mathbf{Q}}_{b}(z)\equiv\left({\mathbb{E}}\left[\frac{g\cdot\mathbf{C}}{1+g\cdot\delta(z)}\right]-z\mathbf{I}_{p}\right)^{-1}$
that characterizes the main bulk,
$\mathbf{V}\equiv[\boldsymbol{\mu},~{}\mathbf{C}^{\frac{1}{2}}\mathbf{U}]\in{\mathbb{R}}^{p\times
3}$ and

$\boldsymbol{\Lambda}(z)\equiv{\mathbb{E}}\frac{g}{1+g\cdot\delta(z)}\begin{bmatrix}1&(\mathbf{U}^{+}\mathbf{z})^{\sf
T}\\\
\mathbf{U}^{+}\mathbf{z}&\mathbf{U}^{+}\mathbf{z}(\mathbf{U}^{+}\mathbf{z})^{\sf
T}-(\mathbf{U}^{\sf T}\mathbf{U})^{+}\end{bmatrix}.$ (33)

Recall from Theorem 1 that $z\in{\mathbb{R}}$ is an eigenvalue within the main
bulk if the Stieltjes transform is not defined

$m(z)=\frac{1}{p}\operatorname{tr}\left({\mathbb{E}}_{g}\frac{g\cdot\mathbf{C}}{1+g\cdot\delta(z)}-z\mathbf{I}_{p}\right)^{-1}=\frac{1}{p}\operatorname{tr}\bar{\mathbf{Q}}_{b}(z).$

To characterize the possible isolated eigenvalues (i.e., the “spikes”), we aim
to find $z\in{\mathbb{R}}$ such that $\bar{\mathbf{Q}}(z)$ is not defined but
$\bar{\mathbf{Q}}_{b}(z)$ is (and thus, outside the main bulk), or
equivalently

$\det\bar{\mathbf{Q}}^{-1}(z)=0,\quad\text{such that
$\det\bar{\mathbf{Q}}_{b}^{-1}(z)\neq 0$ }$

by solving for $z\in{\mathbb{R}}$ such that

$\displaystyle 0$
$\displaystyle=\det\bar{\mathbf{Q}}^{-1}(z)=\det\left(\bar{\mathbf{Q}}^{-1}_{b}(z)+\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{\sf
T}\right)=\det\bar{\mathbf{Q}}^{-1}_{b}(z)\cdot\det(\mathbf{I}_{p}+\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}_{b}(z))$
$\displaystyle=\det(\mathbf{I}_{2K+1}+\boldsymbol{\Lambda}\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}_{b}(z)\mathbf{V})$

with Sylvester’s identity. This concludes the proof of Theorem 2. ∎

###  A.4 Proof of Theorem 3

Similar to the behavior of the isolated eigenvalues studied in Theorem 2, the
corresponding eigenvector can be characterized via the Cauchy’s integration
technique. More precisely, from Theorem 4, recall that
$\mathbf{V}\equiv[\boldsymbol{\mu},~{}\mathbf{C}\mathbf{w}_{*},~{}\mathbf{C}\mathbf{w}]\in{\mathbb{R}}^{p\times
3}$ captures the low rank structure of the Hessian matrix $\mathbf{H}$. Thus,
we can characterize the projection of a given isolated eigenvector
$\hat{\mathbf{u}}$ onto the subspace spanned by the columns of $\mathbf{V}$
(that are assumed to be linearly independent for simplicity of presentation).
To this end, we need to evaluate, for $\lambda$ the asymptotic approximation
given in Theorem 2 of an isolated eigenvalue $\hat{\lambda}$ of $\mathbf{H}$
with corresponding eigenvector $\hat{\mathbf{u}}$, the following matrix

$\mathbf{V}^{\sf T}\hat{\mathbf{u}}\hat{\mathbf{u}}^{\sf
T}\mathbf{V}=-\frac{1}{2\pi\imath}\oint_{\Gamma_{\lambda}}\mathbf{V}^{\sf
T}\mathbf{Q}(z)\mathbf{V}\,dz=-\frac{1}{2\pi\imath}\oint_{\Gamma_{\lambda}}\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}(z)\mathbf{V}\,dz+o(1)$ (34)

where $\Gamma_{\lambda}$ is a positively oriented contour surrounding the
isolated eigenvalue only.

Since the columns of $\mathbf{V}$, denoted $\mathbf{v}_{a}$ for
$a\in\\{1,2,3\\}$, are assumed to be linearly independent, one can write

$\hat{\mathbf{u}}=\sum_{a=1}^{3}\theta_{a}\mathbf{v}_{a}+\boldsymbol{\omega}_{a}$

where $\boldsymbol{\omega}_{a}\in{\mathbb{R}}^{p}$ is a random vector of zero
mean and orthogonal to $\mathbf{v}_{a}$. As a consequence, the projection in
(34) allows to retrieve the coefficient vector
$\boldsymbol{\theta}\in{\mathbb{R}}^{3}$ via

$\mathbf{V}^{\sf T}\hat{\mathbf{u}}\hat{\mathbf{u}}^{\sf
T}\mathbf{V}=\mathbf{V}^{\sf
T}\mathbf{V}\boldsymbol{\theta}\boldsymbol{\theta}^{\sf T}\mathbf{V}^{\sf
T}\mathbf{V}+o(1)$ (35)

with a central limit theorem argument.

We now present the proof of Theorem 3 as follows.

###### Proof of Theorem 3.

Let $(\hat{\lambda},\hat{\mathbf{u}})$ be an isolated eigenvalue-eigenvector
pair of $\mathbf{H}$, and let $\lambda$ denote the asymptotic approximation
given in Theorem 2. We have

$\displaystyle\mathbf{V}^{\sf T}\hat{\mathbf{u}}\hat{\mathbf{u}}^{\sf
T}\mathbf{V}=-\frac{1}{2\pi\imath}\oint_{\Gamma_{\lambda}}\mathbf{V}^{\sf
T}\mathbf{Q}(z)\mathbf{V}\,dz=-\frac{1}{2\pi\imath}\oint_{\Gamma_{\lambda}}\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}(z)\mathbf{V}\,dz+o(1)$
$\displaystyle=-\frac{1}{2\pi\imath}\oint_{\Gamma_{\lambda}}\mathbf{V}^{\sf
T}\left(\bar{\mathbf{Q}}^{-1}_{b}(z)+\mathbf{V}\boldsymbol{\Lambda}(z)\mathbf{V}^{\sf
T}\right)^{-1}\mathbf{V}\,dz+o(1)$
$\displaystyle=-\frac{1}{2\pi\imath}\oint_{\Gamma_{\lambda}}\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}_{b}(z)\mathbf{V}\left(\mathbf{I}_{3}+\boldsymbol{\Lambda}(z)\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}_{b}(z)\mathbf{V}\right)^{-1}\,dz+o(1)$

with Woodbury identity, where we recall
$\bar{\mathbf{Q}}_{b}(z)\equiv\left({\mathbb{E}}\left[\frac{g\cdot\mathbf{C}}{1+g\cdot\delta(z)}\right]-z\mathbf{I}_{p}\right)^{-1}$
defined in (7) that characterizes the main bulk. With Cauchy’s residue
theorem, it remains to evaluate the associated residue and thus the values of
$z$ lying within $\Gamma_{\lambda}$ such that
$\mathbf{I}_{3}+\boldsymbol{\Lambda}\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}_{b}(z)\mathbf{V}$ is singular. As a consequence,

$\displaystyle\mathbf{V}^{\sf T}\hat{\mathbf{u}}\hat{\mathbf{u}}^{\sf
T}\mathbf{V}=-{\rm Res}\left(\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}_{b}(z)\mathbf{V}\cdot\left(\mathbf{I}_{3}+\boldsymbol{\Lambda}(z)\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}_{b}(z)\mathbf{V}\right)^{-1}\right)+o(1)$
$\displaystyle=-\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}_{b}(\lambda)\mathbf{V}\cdot\lim_{z\to\lambda}\left(\mathbf{I}_{3}+\boldsymbol{\Lambda}(z)\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}_{b}(z)\mathbf{V}\right)^{-1}+o(1)$
$\displaystyle=-\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}_{b}(\lambda)\mathbf{V}\cdot\frac{\mathbf{v}_{r,\mathbf{G}}\mathbf{v}_{l,\mathbf{G}}^{\sf
T}}{\mathbf{v}_{l,\mathbf{G}}^{\sf
T}\mathbf{G}^{\prime}(\lambda)\mathbf{v}_{r,\mathbf{G}}}+o(1)\equiv-\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}_{b}(\lambda)\mathbf{V}\cdot\Xi(\lambda)+o(1)$

with $\mathbf{G}(z)$ defined in (10),
$\mathbf{v}_{l,\mathbf{G}},\mathbf{v}_{r,\mathbf{G}}\in{\mathbb{R}}^{3}$ the
left and right eigenvectors of $\mathbf{G}(\lambda)$ associated with
eigenvalue zero, and $\mathbf{G}^{\prime}(z)$ the derivative of
$\mathbf{G}(z)$ with respect to $z$ evaluated at $z=\lambda$.

This concludes the proof of Theorem 3. ∎

Note in particular that

$\mathbf{G}^{\prime}(z)=\boldsymbol{\Lambda}^{\prime}(z)\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}_{b}(z)\mathbf{V}+\boldsymbol{\Lambda}(z)\mathbf{V}^{\sf
T}\bar{\mathbf{Q}}_{b}^{\prime}(z)\mathbf{V}$

with

$\displaystyle\boldsymbol{\Lambda}^{\prime}(z)$
$\displaystyle=-{\mathbb{E}}\frac{g^{2}\cdot\delta^{\prime}(z)}{(1+g\cdot\delta(z))^{2}}\begin{bmatrix}1&(\mathbf{U}^{+}\mathbf{z})^{\sf
T}\\\
\mathbf{U}^{+}\mathbf{z}&(\mathbf{U}^{+}\mathbf{z})(\mathbf{U}^{+}\mathbf{z})^{\sf
T}-(\mathbf{U}^{\sf T}\mathbf{U})^{+}\end{bmatrix},$
$\displaystyle\bar{\mathbf{Q}}_{b}^{\prime}(z)$
$\displaystyle=\bar{\mathbf{Q}}_{b}(z)\left({\mathbb{E}}\frac{g^{2}\cdot\delta^{\prime}(z)}{(1+g\cdot\delta(z))^{2}}\mathbf{C}+\mathbf{I}_{p}\right)\bar{\mathbf{Q}}_{b}(z),$

and
$\delta^{\prime}(z)=\frac{\frac{1}{n}\operatorname{tr}\bar{\mathbf{Q}}_{b}\mathbf{C}\bar{\mathbf{Q}}_{b}}{1-{\mathbb{E}}\left[\frac{g^{2}}{(1+g\cdot\delta)^{2}}\right]\frac{1}{n}\operatorname{tr}\mathbf{C}\bar{\mathbf{Q}}_{b}\mathbf{C}\bar{\mathbf{Q}}_{b}}$.
These _explicit_ derivatives will be used in the numerical evaluations of
Theorem 3 in Section 4.

##  Appendix B Detailed derivations for specific models

###  B.1 Logistic model

Consider the logistic model as in (2)

$\mathbb{P}(y=1\mid\mathbf{w}_{*}^{\sf
T}\mathbf{x})=\sigma(\mathbf{w}_{*}^{\sf
T}\mathbf{x})=1-\mathbb{P}(y=-1\mid\mathbf{w}_{*}^{\sf
T}\mathbf{x}),\quad\sigma(t)=(1+e^{-t})^{-1}$ (36)

with the logistic loss $\ell(y,h)=\ln(1+e^{-yh})$ and $h=\mathbf{w}^{\sf
T}\mathbf{x}\sim\mathcal{N}(\mathbf{w}^{\sf T}\boldsymbol{\mu},\mathbf{w}^{\sf
T}\mathbf{C}\mathbf{w})$, we have for $g\equiv\ell^{\prime\prime}(y,h)$
defined in (8) that

$g=\frac{\partial^{2}\ell(y,h)}{\partial h^{2}}=\frac{e^{yh}}{(1+e^{yh})^{2}}$
(37)

so that for any fixed $\delta(z)$,

${\mathbb{E}}\left[\frac{g}{1+g\delta(z)}\right]={\mathbb{E}}\left[\frac{1}{\delta(z)+2+e^{-yh}+e^{yh}}\right]={\mathbb{E}}[1/(\delta(z)+2+e^{-h}+e^{h})].$
(38)

###  B.2 Phase retrieval model

Consider the noiseless phase retrieval model for which $y=(\mathbf{w}_{*}^{\sf
T}\mathbf{x})^{2}$ with square loss $\ell(y,h)=(y-h^{2})^{2}/4$ and
$h=\mathbf{w}^{\sf T}\mathbf{x}\sim\mathcal{N}(\mathbf{w}^{\sf
T}\boldsymbol{\mu},\mathbf{w}^{\sf T}\mathbf{C}\mathbf{w})$, we have for
$g\equiv\ell^{\prime\prime}(y,h)$ defined in (8) that

$g=\frac{\partial^{2}\ell(y,h)}{\partial h^{2}}=3h^{2}-y$ (39)

so that for any fixed $\delta(z)$,

${\mathbb{E}}\left[\frac{g}{1+g\cdot\delta(z)}\right]={\mathbb{E}}\left[\frac{3h^{2}-y}{1+(3h^{2}-y)\delta(z)}\right]$
(40)

where we note that $3h^{2}-y=\mathbf{x}^{\sf T}\Delta\mathbf{x}$ with
$\Delta\equiv 3\mathbf{w}\mathbf{w}^{\sf T}-\mathbf{w}_{*}\mathbf{w}_{*}^{\sf
T}\in{\mathbb{R}}^{p\times p}$ so that $3h^{2}-y$ is a weighted sum of (at
most) two independent non-central chi-square random variable. (Recall that,
for $\mathbf{x}\sim\mathcal{N}(\boldsymbol{\mu},\mathbf{I}_{p})$,
$\|\mathbf{x}\|^{2}$ follows a non-central $\chi^{2}$ with $p$ degree of
freedom and denotes
$\|\mathbf{x}\|^{2}\sim\chi^{2}(p,\|\boldsymbol{\mu}\|^{2})$.) In particular,
for $\mathbf{w}=\alpha\mathbf{w}_{*}$ with $\alpha\in{\mathbb{R}}$, we obtain
$\Delta=(3\alpha^{2}-1)\mathbf{w}_{*}\mathbf{w}_{*}^{\sf T}$ and $3h^{2}-y$ is
thus a (non-centered) chi-square random variable.

###  B.3 Derivation of Corollary 1

For the logistic model defined in (2) with
$\mathbf{w}=\mathbf{w}_{*}=\mathbf{0}$ and $\mathbf{C}=\mathbf{I}_{p}$, it
follows from Appendix B.1 that $g=1/4$, so that, by Theorem 1, the Stieltjes
transform $m(z)$ satisfies

$zcm^{2}(z)-(1-c-4z)m(z)+4=0$ (41)

which corresponds to the Marc̆enko-Pastur law in (12) rescaled by $1/4$. It
then follows from Theorem 2 that

$\det\mathbf{G}(\lambda)=1+\frac{\|\boldsymbol{\mu}\|^{2}m(\lambda)}{4+cm(\lambda)}=0\Leftrightarrow
m(\lambda)=-\frac{4}{c+\|\boldsymbol{\mu}\|^{2}}$ (42)

which, by plugging back into the equation above gives the expression of
$\lambda_{\boldsymbol{\mu}}$, and the condition
$\|\boldsymbol{\mu}\|^{2}>\sqrt{c}$ follows from the property $zm(z)<0$ for
$\Im[z]=0$ of the Stieltjes transform. With Theorem 3, we conclude the proof
of Corollary 1.

###  B.4 Derivation of Corollary 2

For the logistic model defined in (2) with $\boldsymbol{\mu}=\mathbf{0}$ and
$\mathbf{C}=\mathbf{I}_{p}$, it follows from Appendix B.1 and Theorem 1 that
the Stieltjes transform $m(z)$ satisfies

$m(z)=\frac{\delta(z)}{c}=\frac{1}{{\mathbb{E}}[(2+\exp(-r)+\exp(r)+cm(z))^{-1}]-z}$
(43)

for $r\sim\mathcal{N}(0,\|\mathbf{w}\|^{2})$ that _depends_ on $\mathbf{w}$
but not on $\mathbf{w}_{*}$. Moreover, Theorem 2 gives
$\bar{\mathbf{Q}}_{b}(z)=m(z)\mathbf{I}_{p}$ and

$\displaystyle 0$
$\displaystyle=\det\mathbf{G}(z)=\det\left(\mathbf{I}_{2}+m(z){\mathbb{E}}[f(r,z)(r^{2}-\|\mathbf{w}\|^{2})](\mathbf{U}^{\sf
T}\mathbf{U})^{+}\mathbf{U}^{\sf T}\frac{\mathbf{w}\mathbf{w}^{\sf
T}}{\|\mathbf{w}\|^{4}}\mathbf{U}\right)$
$\displaystyle=1+m(z)\frac{{\mathbb{E}}[f(r,z)(r^{2}-\|\mathbf{w}\|^{2})]}{\|\mathbf{w}\|^{2}}$

with $\mathbf{U}=[\mathbf{w}_{*},~{}\mathbf{w}]$ and
$f(t,z)=\frac{1}{cm(z)+2+\exp(-t)+\exp(t)}$, where we use the fact that

$\displaystyle{\mathbb{E}}\left[\frac{g\cdot\mathbf{z}\mathbf{z}^{\sf
T}}{1+g\cdot\delta(z)}\right]$
$\displaystyle={\mathbb{E}}[f(r,z)]\cdot\mathbf{I}_{p}+{\mathbb{E}}[f(r,z)(r^{2}-\|\mathbf{w}\|^{2})]\cdot\frac{\mathbf{w}\mathbf{w}^{\sf
T}}{\|\mathbf{w}\|^{4}},$

and ${\mathbb{E}}[\frac{g}{1+g\cdot\delta(z)}]={\mathbb{E}}[f(r,z)]$ from
Appendix B.1. As such, there exists _no_ spike due to $\mathbf{w}_{*}$.

Moreover, taking $\mathbf{w}_{*}=\mathbf{0}$, we further obtain

$\mathbf{G}(z)=1+m(z){\mathbb{E}}[f(r,z)(r^{2}/\|\mathbf{w}\|^{2}-1)]$ (44)

and

$\mathbf{G}^{\prime}(z)=m^{\prime}(z)\left({\mathbb{E}}[f(r,z)(r^{2}/\|\mathbf{w}\|^{2}-1)]-cm(z){\mathbb{E}}[f^{2}(r,z)(r^{2}/\|\mathbf{w}\|^{2}-1)]\right)$
(45)

where we have
$m^{\prime}(z)=\frac{m^{2}(z)}{1-cm^{2}(z){\mathbb{E}}[f^{2}(r,z)]}$ by
differentiating (43). Therefore, we obtain from Theorem 3 that

$\frac{|\mathbf{w}^{\sf
T}\hat{\mathbf{u}}|^{2}}{\|\mathbf{w}\|^{2}}=-\frac{m(\lambda_{\mathbf{w}})}{\mathbf{G}^{\prime}(\lambda_{\mathbf{w}})}+o(1).$
(46)