# On the Inherent Regularization Effects of
Noise Injection During Training

Oussama Dhifallah and Yue M. Lu O. Dhifallah and Y. M. Lu are with the John A.
Paulson School of Engineering and Applied Sciences, Harvard University,
Cambridge, MA 02138, USA (e-mails:
oussama$\\_$dhifallah@g.harvard.edu,yuelu@seas.harvard.edu).This research was
funded by the Harvard FAS Dean’s Fund for Promising Scholarship, and by the US
National Science Foundations under grants CCF-1718698 and CCF-1910410.

###### Abstract

Randomly perturbing networks during the training process is a commonly used
approach to improving generalization performance. In this paper, we present a
theoretical study of one particular way of random perturbation, which
corresponds to injecting artificial noise to the training data. We provide a
precise asymptotic characterization of the training and generalization errors
of such randomly perturbed learning problems on a random feature model. Our
analysis shows that Gaussian noise injection in the training process is
equivalent to introducing a weighted ridge regularization, when the number of
noise injections tends to infinity. The explicit form of the regularization is
also given. Numerical results corroborate our asymptotic predictions, showing
that they are accurate even in moderate problem dimensions. Our theoretical
predictions are based on a new correlated Gaussian equivalence conjecture that
generalizes recent results in the study of random feature models.

##  I Introduction

A popular approach to improving the generalization performance is to randomly
perturb the network during the training process [1, 2, 3, 4, 5]. Such random
perturbations are widely used as an implicit regularization to the learning
problem. One way that random perturbation has been used as a regularization is
by injecting it to the input data before starting the learning process [6, 7,
8]. In this paper, we provide a theoretical analysis of such learning
procedure on a random feature model [9] under Gaussian input and perturbation
vectors. Our analysis particularly shows that Gaussian noise injection
introduces a weighted ridge regularization, asymptotically.

First, we describe the models for our theoretical analysis. We are given a
collection of training data $\\{(y_{i},\boldsymbol{a}_{i})\\}_{i=1}^{n}$,
where $\boldsymbol{a}_{i}\in\mathbb{R}^{p}$ is referred to as the input vector
and $y_{i}\in\mathbb{R}$ is referred to as the label corresponding to
$\boldsymbol{a}_{i}$. In this paper, we shall assume that the labels are
generated according to the standard teacher–student model, i.e.

$\displaystyle
y_{i}=\varphi(\boldsymbol{a}_{i}^{\top}\boldsymbol{\xi}),~{}\forall
i\in\\{1,\dots,n\\},$ (1)

where $\boldsymbol{\xi}\in\mathbb{R}^{p}$ is an unknown teacher weight vector,
and $\varphi(\cdot)$ is a scalar deterministic or probabilistic function.
Here, we use the random feature model [9] to learn the model described in (1).
The random feature model considers the following class of functions

$\displaystyle\mathcal{F}_{\text{RF}}(\boldsymbol{a})=\Big{\\{}g_{\boldsymbol{w}}(\boldsymbol{a})=\boldsymbol{w}^{\top}\sigma(\boldsymbol{F}^{\top}\boldsymbol{a}),~{}\boldsymbol{w}\in\mathbb{R}^{k}\Big{\\}},$
(2)

where $\boldsymbol{a}\in\mathbb{R}^{p}$ is an input vector,
$\boldsymbol{F}\in\mathbb{R}^{p\times k}$ is a random matrix referred to as
the feature matrix, and $\sigma(\cdot)$ is a scalar function referred to as
the activation function. This model assumes that $\boldsymbol{F}$ is fixed
during the training. Note that the family in (2) can be viewed as a two–layer
neural network where the first layer weights are fixed, i.e. $\boldsymbol{F}$
is fixed.

###  I-A Learning Formulation

Before starting the learning process, $\ell$ independent perturbation vectors
are injected to each $\boldsymbol{a}_{i}$. This procedure forms the augmented
family $\\{\boldsymbol{a}_{i}+\Delta\boldsymbol{z}_{ij}\\}_{j=1}^{\ell}$ for
each $\boldsymbol{a}_{i}$, where $\\{\boldsymbol{z}_{ij}\\}_{j=1}^{\ell}$ are
independent random perturbations and $\Delta\geq 0$ denotes the noise
variance. In this paper, we study the effects of such perturbation method on
an average loss and a random feature model. Specifically, we analyze
formulations of the following form

$\displaystyle\widehat{\boldsymbol{w}}=\operatornamewithlimits{argmin}_{\boldsymbol{w}\in\mathbb{R}^{k}}$
$\displaystyle~{}\frac{1}{2n\ell}\sum_{i=1}^{n}\sum_{j=1}^{\ell}\big{(}y_{i}-\boldsymbol{w}^{\top}\sigma\big{(}\boldsymbol{F}^{\top}[\boldsymbol{a}_{i}+\Delta\boldsymbol{z}_{ij}]\big{)}\big{)}^{2}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2},$
(3)

where $\lambda>0$ denotes the regularization parameter. Note that the problem
in (3) is a standard feature formulation when $\Delta=0$. Then, we refer to
(3) as the noisy formulation, when $\Delta>0$ and the standard formulation,
otherwise.

###  I-B Performance Measure

The main objective in this paper is to study the performance of the learning
formulation in (3) on unobserved test data. For every test vector
$\boldsymbol{a}_{\text{new}}\in\mathbb{R}^{p}$, the corresponding label
$\widehat{y}$ can be predicted using the following (probabilistic) role

$\displaystyle\widehat{y}=\widehat{\varphi}[\widehat{\boldsymbol{w}}^{\top}\sigma(\boldsymbol{F}^{\top}\boldsymbol{a}_{\text{new}})],$
(4)

for some predefined function $\widehat{\varphi}(\cdot)$, where
$\widehat{\boldsymbol{w}}\in\mathbb{R}^{k}$ denotes the optimal solution of
the formulation given in (3). To measure the performance of the learning
problem in (3) on any unobserved test data
$\\{({y}_{\text{new}},\boldsymbol{a}_{\text{new}})\\}$, we use the
generalization error defined as follows

$\displaystyle\mathcal{E}_{\text{test}}=\frac{1}{4^{\upsilon}}\mathbb{E}\Big{[}\big{(}{y}_{\text{new}}-\widehat{\varphi}(\widehat{\boldsymbol{w}}^{\top}\sigma(\boldsymbol{F}^{\top}\boldsymbol{a}_{\text{new}}))\big{)}^{2}\Big{]}.$
(5)

Here, the expectation is taken over the distribution of the unobserved test
vector $\boldsymbol{a}_{\text{new}}$ and the (random) functions
$\varphi(\cdot)$ and $\widehat{\varphi}(\cdot)$. We take $\upsilon=0$ for
regression problems (e.g. $\varphi(\cdot)$ is the identity function) and
$\upsilon=1$ for binary classification problems (e.g. $\varphi(\cdot)$ is the
sign function). In this paper, we assume that the test data is generated
according to the same training model introduced in (1). Furthermore, we
measure the performance of the formulation in (3) on the training data via the
training error defined as follows

$\displaystyle\mathcal{E}_{\text{train}}=\frac{1}{2n\ell}\sum_{i=1}^{n}\sum_{j=1}^{\ell}\big{(}y_{i}-\widehat{\boldsymbol{w}}^{\top}\sigma\big{(}\boldsymbol{F}^{\top}[\boldsymbol{a}_{i}+\Delta\boldsymbol{z}_{ij}]\big{)}\big{)}^{2}.$

Note that the training error is the optimal cost value of our learning
formulation in (3) without regularization.

###  I-C Contributions

The contribution of this paper can be summarized as follows:

  * (C.1)

Our first contribution is a correlated Gaussian equivalence conjecture (cGEC).
Our conjecture considers Gaussian input and perturbation vectors. It states
that the learning formulation in (3) is asymptotically equivalent to a simpler
optimization problem that can be formulated by replacing the non–linear
vectors

$\displaystyle\boldsymbol{v}_{ij}=\sigma\big{(}\boldsymbol{F}^{\top}[\boldsymbol{a}_{i}+\Delta\boldsymbol{z}_{ij}]\big{)},$

with linear vectors with the following form

$\displaystyle\boldsymbol{q}_{ij}=\mu_{0}\boldsymbol{1}_{k}+\widetilde{\mu}_{1}\boldsymbol{F}^{\top}\boldsymbol{a}_{i}+\widehat{\mu}_{1}\boldsymbol{F}^{\top}\boldsymbol{z}_{ij}+\mu_{2}\boldsymbol{b}_{i}+\mu_{3}\boldsymbol{p}_{ij}.$

Here, $\\{\boldsymbol{b}_{i}\\}_{i=1}^{n}$ and
$\\{\boldsymbol{p}_{ij}\\}_{i,j=1}^{n,\ell}$ are independent standard Gaussian
random vectors and independent of $\\{\boldsymbol{a}_{i}\\}_{i=1}^{n}$ and
$\\{\boldsymbol{z}_{ij}\\}_{i,j=1}^{n,\ell}$. Moreover, the weights $\mu_{0}$,
$\widetilde{\mu}_{1}$, $\widehat{\mu}_{1}$, $\mu_{2}$ and $\mu_{3}$ depend on
$\sigma(\cdot)$ and $\Delta$ as follows

$\displaystyle\begin{cases}\mu_{0}=\mathbb{E}[\sigma(x_{1})],\widetilde{\mu}_{1}=\mathbb{E}[z\sigma(x_{1})],\widehat{\mu}_{1}=\mathbb{E}[v_{1}\sigma(x_{1})]\\\
\mu_{2}^{2}=\mathbb{E}[\sigma(x_{1})\sigma(x_{2})]-\mu_{0}^{2}-\widetilde{\mu}_{1}^{2}\\\
\mu_{3}^{2}=\mathbb{E}[\sigma(x_{1})^{2}]-\mathbb{E}[\sigma(x_{1})\sigma(x_{2})]-\widehat{\mu}_{1}^{2},\end{cases}$

where $x_{1}=z+\Delta v_{1}$, $x_{2}=z+\Delta v_{2}$, and $z$, $v_{1}$ and
$v_{2}$ are independent standard Gaussian random variables. Specifically, the
cGEC states that the performance of the formulation:

$\displaystyle\min_{\boldsymbol{w}\in\mathbb{R}^{k}}\frac{1}{2n\ell}\sum_{i=1}^{n}\sum_{j=1}^{\ell}\big{(}y_{i}-\widetilde{\mu}_{1}\boldsymbol{w}^{\top}\boldsymbol{F}^{\top}\boldsymbol{a}_{i}-\widehat{\mu}_{1}\boldsymbol{w}^{\top}\boldsymbol{F}^{\top}\boldsymbol{z}_{ij}$
$\displaystyle-\mu_{0}\boldsymbol{w}^{\top}\boldsymbol{1}_{k}-\mu_{2}\boldsymbol{w}^{\top}\boldsymbol{b}_{i}-\mu_{3}\boldsymbol{w}^{\top}\boldsymbol{p}_{ij}\big{)}^{2}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2},$
(6)

is asymptotically equivalent to the performance of the noisy formulation. This
conjecture is valid in the asymptotic limit (i.e. $n$, $p$ and $k$ grow to
infinity at finite ratios). More details about this equivalence is provided in
Section II. We refer to ((C.1)) as the Gaussian formulation. The cGEC is
verified by presenting multiple simulations in different scenarios.

  * (C.2)

The second contribution is a precise characterization of the training and
generalization errors of the noise injection procedure formulated in (3) for
Gaussian input and perturbation vectors. Our analysis is based on the cGEC and
valid in the high–dimensional setting (i.e. $n$, $p$ and $k$ grow to infinity
at finite ratios). Our predictions show that the asymptotic limit of the
training and generalization errors can be precisely predicted after solving a
scalar deterministic formulation. The theoretical predictions are obtained
using an extended version of the convex Gaussian min-max theorem (CGMT) [10,
11] which we refer to as the multivariate CGMT. The new version of the CGMT
accounts for the correlation introduced by injecting Gaussian noise during the
learning process. Our asymptotic results hold for a general family of feature
matrices, activation functions and generative models satisfying (1).

Figure 1: Solid line: Theoretical predictions. Circle: Numerical simulations
for (3). Black cross: Numerical simulations for ((C.1)). $\varphi(\cdot)$ is
the sign function with probability $\theta$ of flipping the sign.
$\widehat{\varphi}(\cdot)$ and $\sigma(\cdot)$ are the sign function. We set
$p=500$, $\Delta=0.5$, $\alpha=2$, $\theta=0.1$, $\lambda=10^{-5}$.
$\boldsymbol{F}$ has independent Gaussian components with zero mean and
variance $1/p$. The results are averaged over $200$ independent Monte Carlo
trials.

In Figure 1, we compare our theoretical predictions with the actual
performance of the learning problem given in (3). First, note that our
asymptotic predictions are in excellent agreement with the actual performance
of (3) and its Gaussian formulation given in ((C.1)), even for moderate values
of $p$, $n$ and $k$. This provides a first empirical validation of our
results. Figure 1 also study the effects of $\ell$ on the training and
generalization performance. Note that the generalization error follows a
double descent curve [12, 13]. Specifically, the generalization error
decreases monotonically as a function of the complexity $\eta=k/n$ after
reaching a peak known as the interpolation threshold [12, 13]. Figure 1
particularly demonstrates that the location of the interpolation threshold
depends on the number of noise samples. Specifically, the interpolation
threshold peak occurs at $\ell$ for fixed noise variance $\Delta=0.5$.
Additionally, Figure 1 shows that the interpolation threshold occurs when the
training error converges to zero. Then, we can see that perturbing the input
data with $\ell$ random noise vectors moves the interpolation threshold from
$1$ to $\ell$ and improves the generalization error for complexity $\eta$
lower than $\ell$.

  * (C.3)

The third contribution is a precise analysis of the regularization effects of
the considered noise injection procedure. Specifically, we use the asymptotic
predictions of the noisy formulation to show that the noise injection model in
(3) is equivalent to solving a standard feature formulation with an additional
weighted ridge regularization. This theoretical result is valid when the
number of noise samples $\ell$ tends to infinity. In particular, we show that
the formulation in (3) is equivalent to solving the problem

$\displaystyle\min_{\boldsymbol{w}\in\mathbb{R}^{k}}$
$\displaystyle~{}\frac{1}{2n}\sum_{i=1}^{n}\big{(}y_{i}-\boldsymbol{w}^{\top}\widehat{\sigma}(\boldsymbol{F}^{\top}\boldsymbol{a}_{i})\big{)}^{2}+\frac{1}{2}\mathinner{\\!\left\lVert\boldsymbol{R}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}^{2}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2},$
(7)

when $\ell$ grows to infinity slower than the dimensions $n$, $p$ and $k$.
Here, $\widehat{\sigma}(\cdot)$ is a new activation function and
$\boldsymbol{R}$ is defined as follows

$\displaystyle\boldsymbol{R}=\widehat{\mu}_{1}^{2}\boldsymbol{F}^{\top}\boldsymbol{F}+\mu_{3}^{2}\boldsymbol{I}_{k}.$
(8)

Finally, we provide a precise asymptotic characterization of the training and
generalization errors corresponding to (7). We refer to this formulation as
the limiting formulation.

Figure 2: Solid line: Theoretical predictions. Circle: Numerical simulations
for (3) and (7). $\varphi(\cdot)$, $\widehat{\varphi}(\cdot)$ and
$\sigma(\cdot)$ are the sign function. (a) $p=700$, $\ell=180$, $\alpha=1$,
$\Delta=1$ and $\lambda=10^{-3}$. (b) $p=600$, $\alpha=1.5$, $\Delta=1$ and
$\lambda=10^{-3}$. $\boldsymbol{F}$ has independent Gaussian components with
zero mean and variance $1/p$. The results are averaged over $100$ independent
Monte Carlo trials.

Figure 2 provides another empirical verification of our theoretical
predictions since it shows that they are in excellent agreement with the
actual performance of (3) and (7). Figure 2 shows that the noisy formulation
in (3) has approximately the same performance as the formulation in (7) for
$\ell=180$. This is aligned with our theoretical prediction which states that
the formulations in (3) and (7) are equivalent when $\ell$ grows to infinity
slower than the dimensions $n$, $p$ and $k$. Figure 2 illustrates the
convergence behavior of the generalization error corresponding to (3) for a
fixed value of $\eta$. It particularly shows that the noisy formulation has a
good convergence rate, i.e. the limit is already attained with a moderate
value of $\ell$. Moreover, we can see from Figure 2 that the convergence rate
depends on the complexity parameter $\eta$.

###  I-D Related Work

There has been significant interest in precisely characterizing the
performance of the random feature model in recent literature [14, 15, 16, 17].
The ridge regression formulation, (i.e. $\varphi(\cdot)$ is the identity
function and $\Delta=0$ in (3)) is precisely analyzed in [14] where the
feature matrix is Gaussian. In a subsequent work, [18] uses the CGMT to
accurately analyze the maximum-margin linear classifier in the
overparametrized regime. The work in [15] precisely characterizes the
performance of the standard formulation, i.e. $\Delta=0$, for general families
of feature matrices and convex loss functions. The results presented in [15]
are derived using the non–rigorous replica method [19]. The predictions in
[15] are rigiourously verified in [16] using the CGMT. All the previous work
consider an unperturbed formulation of the random feature model. In this
paper, we study the effects of adding random noise during training. Our
analysis is based on an extended version of the CGMT referred to as the
multivariate CGMT. The CGMT is first used in [20] and further developed in
[10]. It extends a Gaussian theorem first introduced in [21]. It relies on
(strong) convexity properties to prove an equivalence between two Gaussian
processes. It has been successfully applied in the analysis of convex
regression [10, 22, 16] and convex classification [23, 24, 25, 26]
formulations.

There has been significant interest in studying the effects of random noise
injection during training (see e.g. [2, 27, 3]). In particular, prior
literature [28, 29] shows that Gaussian noise injection during training
improves the robustness of the network. Moreover, several recent papers [2, 6]
show that such perturbation technique introduces some sort of regularization
to the loss function. In particular, the work in [6] shows that minimizing the
worst–case loss introduces a gradient norm regularization.

Another popular perturbation approach used in regularizing learning models is
the dropout method [1, 30]. It consists of perturbing the learning problem by
randomly dropping units from the network during the training procedure. In
this paper, we precisely analyze the Gaussian noise injection method and we
leave the analysis of the dropout technique for future work. Our empirical
studies suggest that the dropout method has a better convergence rate as
compared to the noisy formulation. Moreover, they suggest that both methods
have comparable generalization performance.

###  I-E Organization

The rest of this paper is organized as follows. Section II provides more
details about the cGEC. Section III lay out the technical assumptions under
which our results are derived. Section IV provides an asymptotic
characterization of the noisy formulation. Our theoretical predictions hold
for a general family of feature matrices, activation functions and generative
models as in (1). We provide additional simulation examples for special cases
of our results in Section V. The detailed proof of our theoretical predictions
is provided in Section VI. Section VII concludes the paper. The appendix in
Section VIII provides additional technical details.

##  II Gaussian Equivalence Conjecture with an Intuitive Explanation

Consider three independent standard Gaussian random vectors
$\boldsymbol{a}\in\mathbb{R}^{p}$, $\boldsymbol{z}_{1}\in\mathbb{R}^{p}$ and
$\boldsymbol{z}_{2}\in\mathbb{R}^{p}$. Moreover, consider the random variables
$\nu_{1}=\boldsymbol{\xi}^{\top}\boldsymbol{a}$, $\nu_{2}$ and $\nu_{3}$
defined as follows

$\displaystyle\nu_{2}=\boldsymbol{w}^{\top}\sigma(\boldsymbol{F}^{\top}[\boldsymbol{a}+\Delta\boldsymbol{z}_{1}]),\nu_{3}=\boldsymbol{w}^{\top}\sigma(\boldsymbol{F}^{\top}[\boldsymbol{a}+\Delta\boldsymbol{z}_{2}]),$

where $\sigma(\cdot)$, $\boldsymbol{\xi}\in\mathbb{R}^{p}$ and
$\boldsymbol{F}\in\mathbb{R}^{p\times k}$ satisfy some regularity assumptions,
and where $\boldsymbol{w}\in\mathbb{R}^{k}$. Moreover, define the joint
probability distribution of $\nu_{1}$, $\nu_{2}$ and $\nu_{3}$ as
$\mathbb{P}(\nu_{1},\nu_{2},\nu_{3})$. The cGEC states that the joint
distribution $\mathbb{P}(\nu_{1},\nu_{2},\nu_{3})$ is asymptotically Gaussian,
i.e.
$d(\mathbb{P}(\nu_{1},\nu_{2},\nu_{3}),\mathbb{P}(\nu_{g,1},\nu_{g,2},\nu_{g,3}))$
converges in probability to zero where $\nu_{g,1}$, $\nu_{g,2}$ and
$\nu_{g,3}$ are jointly Gaussian with the same first and second moments of
$\nu_{1}$ , $\nu_{2}$ and $\nu_{3}$ and $d(\cdot,\cdot)$ is some probability
distance that metrizes the convergence in distribution (e.g maximum-sliced
(MS) distance [31, 32]). To have the same first two moments, the random
variables ${\nu}_{g,1}$, ${\nu}_{g,2}$ and ${\nu}_{g,3}$ are selected as
follows $\nu_{g,1}=\nu_{1}$ and

$\displaystyle{\nu}_{g,2}=\boldsymbol{w}^{\top}\big{(}\mu_{0}\boldsymbol{1}_{k}+\boldsymbol{F}^{\top}[\widetilde{\mu}_{1}\boldsymbol{a}+\widehat{\mu}_{1}\boldsymbol{z}_{1}]+\mu_{2}\boldsymbol{b}+\mu_{3}\boldsymbol{p}_{1}\big{)},$
$\displaystyle{\nu}_{g,3}=\boldsymbol{w}^{\top}\big{(}\mu_{0}\boldsymbol{1}_{k}+\boldsymbol{F}^{\top}[\widetilde{\mu}_{1}\boldsymbol{a}+\widehat{\mu}_{1}\boldsymbol{z}_{2}]+\mu_{2}\boldsymbol{b}+\mu_{3}\boldsymbol{p}_{2}\big{)},$

where $\boldsymbol{1}_{k}$ represents the all $1$ vector with size $k$. Here,
$\boldsymbol{b}\in\mathbb{R}^{k}$, $\boldsymbol{p}_{1}\in\mathbb{R}^{k}$ and
$\boldsymbol{p}_{2}\in\mathbb{R}^{k}$ are three independent standard Gaussian
random vectors and they are independent of $\boldsymbol{a}$,
$\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$. The weights $\mu_{0}$,
$\widetilde{\mu}_{1}$, $\widehat{\mu}_{1}$, $\mu_{2}$ and $\mu_{3}$ are as
defined in Section I-C.

In the standard setting, i.e. $\Delta=0$, the cGEC is equivalent to the
uniform Gaussian equivalence theorem (uGET), observed and used in many earlier
papers [18, 15, 33, 16]. Recently, the work in [17] provided a rigorous proof
of the uGET. Specifically, the work in [17] proves a special case of cGEC when
$\Delta=0$, the feature matrix is Gaussian and the activation functions have
bounded first three derivatives. However, similar to previous literature [33],
we conjecture that the cGEC is valid under more general settings. We believe
that the analysis in [17] can be extended to prove the cGEC and we leave the
full technical details for future work.

Our theoretical results are based on this conjecture. It is thus useful to
also provide an intuitive explanation for the plausibility of the cGEC. Assume
that $\boldsymbol{f}_{i}$ is the $i$th column of $\boldsymbol{F}$. The
nonlinear term
$I_{1}=\sigma(\boldsymbol{f}_{i}^{\top}(\boldsymbol{a}+\Delta\boldsymbol{z}_{1}))$
can be decomposed by projecting on the basis
$(1,\boldsymbol{f}_{i}^{\top}\boldsymbol{a},\boldsymbol{f}_{i}^{\top}\boldsymbol{z}_{1})$,
i.e.
$I_{1}=\mu_{0}+\widetilde{\mu}_{1}\boldsymbol{f}_{i}^{\top}\boldsymbol{a}+\widehat{\mu}_{1}\boldsymbol{f}_{i}^{\top}\boldsymbol{z}_{1}+\sigma^{\perp}_{i}$.
The term $\sigma^{\perp}_{i}$ is selected so that we match the variance of
$I_{1}$ and the correlation with
$I_{2}=\sigma(\boldsymbol{f}_{i}^{\top}(\boldsymbol{a}+\Delta\boldsymbol{z}_{2}))$.
We note that the cGEC makes sense when the columns of $\boldsymbol{F}$ are
independent and have the same norm. These are the regularity assumptions for
the feature matrix in [33]. The same intuition also appears in the analysis of
the unperturbed random kernel models, in particular, the random feature model
[18]. In this paper, we suppose that the feature matrix and the activation
function satisfy the regularity assumptions in [33] and conjecture that the
Gaussian equivalence is valid for $(\nu_{1},\nu_{2},\dots,\nu_{\ell})$ for
$\ell\geq 1$ and uniformly in $\boldsymbol{w}\in\mathbb{R}^{k}$. Using the
cGEC, the performance of the formulation in (3) can be fully characterized by
asymptotically analyzing the Gaussian formulation given in ((C.1)). We verify
this conjecture by performing multiple simulation examples in different
settings (regression and classification).

##  III Technical Assumptions

In this paper, we precisely characterize the noisy formulation under the
following technical assumptions.

######  Assumption 1 (Gaussian Vectors).

The input vectors $\\{\boldsymbol{a}_{i}\\}_{i=1}^{n}$ and the perturbation
vectors $\\{\boldsymbol{z}_{ij}\\}_{i=1,j=1}^{n,\ell}$ are known and drawn
independently from a standard Gaussian distribution. Without loss of
generality, we assume that the hidden vector
$\boldsymbol{\xi}\in\mathbb{R}^{p}$ has unit norm. Also, it is independent of
the input vectors, the noise vectors and $\boldsymbol{F}$.

Our theoretical predictions are valid in the high-dimensional setting where
$n$, $p$ and $k$ grow to infinity at finite ratios.

######  Assumption 2 (Asymptotic Limit).

The number of samples and the number of hidden neurons satisfy $n=n(p)$ and
$k=k(p)$, respectively. We assume that $\alpha_{p}=n(p)/p\to\alpha>0$ and
$\eta_{p}=k(p)/n(p)\to\eta>0$ as $p\to\infty$. Also, the number of noise
injections $\ell$ is independent of $p$.

Moreover, we consider the following assumption to ensure that the
generalization error defined in (5) concentrates in the high–dimensional
limit.

######  Assumption 3 (Generative Model).

The data generating function $\varphi(\cdot)$ introduced in (1) is independent
of the input vectors, the noise vectors and the feature matrix. Moreover, the
following conditions are satisfied.

  * (a)

$\varphi(\cdot)$ and $\widehat{\varphi}(\cdot)$ are continuous almost
everywhere in $\mathbb{R}$. For every $h>0$ and $z\sim\mathcal{N}(0,h)$, we
have $\mathbb{E}[\varphi^{2}(z)]<+\infty$, $\mathbb{E}[z\varphi(z)]\neq 0$ and
$0<\mathbb{E}[\widehat{\varphi}^{2}(z)]<+\infty$.

  * (b)

For any compact interval $[c,C]$, there exists a function $g(\cdot)$ such that

$\sup_{h,\chi\in[c,C]}\mathinner{\\!\left\lvert\widehat{\varphi}(\chi+hx)\right\rvert}^{2}\leq
g(x)\quad\text{for all }x\in\mathbb{R}.$

Additionally, the function $g(\cdot)$ satisfies
$\mathbb{E}[{{g}^{2}(z)}]<+\infty$, where $z\sim\mathcal{N}(0,1)$.

In addition to the assumptions in Section II, we consider the following
regularity conditions for the activation function.

######  Assumption 4 (Activation Function).

The activation function $\sigma(\cdot)$ is independent of the input vectors,
the noise vectors and the feature matrix. It also satisfies
$\mathbb{E}[\sigma(z)^{2}]<+\infty$ and $\mathbb{E}[z\sigma(z)]\neq 0$, where
$z\sim\mathcal{N}(0,1)$.

In addition to the assumptions discussed in Section II, we consider a family
of feature matrices that satisfy the following assumption to guarantee that
the Gaussian formulation converges to a deterministic problem.

######  Assumption 5 (Feature Matrix).

The SVD decomposition of the feature matrix can be expressed as
$\boldsymbol{F}=\boldsymbol{U}\boldsymbol{S}\boldsymbol{V}$, where
$\boldsymbol{U}\in\mathbb{R}^{p\times p}$ and
$\boldsymbol{V}\in\mathbb{R}^{k\times k}$ are random orthogonal matrices and
$\boldsymbol{S}\in\mathbb{R}^{p\times k}$ is a diagonal matrix formed by the
singular values of $\boldsymbol{F}$. Define the matrix $\boldsymbol{M}$ as
$\boldsymbol{M}=\boldsymbol{F}^{\top}\boldsymbol{F}$.

  * (a)

We assume that $\boldsymbol{U}$ is a Haar-distributed random unitary matrix.

  * (b)

We also assume that the empirical distribution of the eigenvalues of the
matrix $\boldsymbol{M}$ converges weakly to a probability distribution
$\mathbb{P}_{\kappa}(\cdot)$ supported in $[0~{}\zeta_{\text{max}}]$, where
$\zeta_{\text{max}}>0$ is a constant independent of $p$ and $\ell$.

  * (c)

We finally assume that $\mathbb{E}_{\kappa}[\kappa]>0$, where the expectation
is taken over the distribution $\mathbb{P}_{\kappa}(\cdot)$.

Based on Assumption 2, we also have the following property
$\delta_{p}=k(p)/p\to\delta>0$ as $p$ grows to infinity. Our theoretical
predictions use the weak convergence in Assumption 5 to fully characterize the
noisy formulation.

##  IV Precise Analysis of the Noisy Formulation

In this section, we asymptotically analyze the noise injection procedure
introduced in (3). Specifically, we provide a precise asymptotic
characterization of the training and generalization errors corresponding to
(3).

###  IV-A Precise Asymptotic Analysis

Before stating our technical results, we start with few definitions. Define
the following two deterministic functions

$\displaystyle
T_{2,\lambda}(\boldsymbol{t},\boldsymbol{\tau})=\frac{\delta}{T_{1}^{2}}\mathbb{E}_{\kappa}\Big{[}\frac{\kappa}{g_{\kappa,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}\Big{]}/\Big{(}1-\frac{\widetilde{\mu}_{1}^{2}t_{1}\delta}{\tau_{1}}\mathbb{E}_{\kappa}\Big{[}\frac{\kappa}{g_{\kappa,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}\Big{]}\Big{)}$
$\displaystyle
T_{3,\lambda}(\boldsymbol{t},\boldsymbol{\tau})=\frac{t_{1}^{2}}{\ell}\mathbb{E}_{\kappa}\Big{[}\frac{\widetilde{\mu}_{1}^{2}\kappa+\mu_{2}^{2}}{g_{\kappa,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}\Big{]}+\frac{t_{1}^{2}+t_{2}^{2}}{\ell^{2}}\mathbb{E}_{\kappa}\Big{[}\frac{\widehat{\mu}_{1}^{2}\kappa+\mu_{3}^{2}}{g_{\kappa,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}\Big{]},$

where the expectations are taken over the probability distribution
$\mathbb{P}_{\kappa}(\cdot)$ defined in Assumption 5 and where
$\boldsymbol{t}=[t_{1},t_{2}]^{\top}$ and
$\boldsymbol{\tau}=[\tau_{1},\tau_{2}]^{\top}$. Here, the function
$g_{\kappa,\lambda}(\cdot,\cdot)$ is defined as follows

$\displaystyle g_{\kappa,\lambda}(\boldsymbol{t},\boldsymbol{\tau})$
$\displaystyle=\frac{t_{1}}{\tau_{1}}\big{(}\widetilde{\mu}_{1}^{2}\kappa+\mu_{2}^{2}\big{)}+\Big{(}\frac{t_{1}}{\tau_{1}\ell}+\frac{t_{2}(\ell-1)}{\tau_{2}\ell}\Big{)}\big{(}\widehat{\mu}_{1}^{2}\kappa+\mu_{3}^{2}\big{)}+\lambda.$
(9)

Furthermore, define the following four-dimensional deterministic optimization
problem

$\displaystyle\inf_{\begin{subarray}{c}\tau_{1}>0\\\
\tau_{2}>0\end{subarray}}\max_{\begin{subarray}{c}t_{1}\geq 0\\\ t_{2}\geq
0\end{subarray}}$
$\displaystyle~{}\frac{t_{1}}{2\tau_{1}}\big{(}\gamma_{1}-2\widetilde{\mu}_{1}T_{1}q_{t,\tau}^{\star}\gamma_{2}+\widetilde{\mu}_{1}^{2}T_{1}^{2}(q_{t,\tau}^{\star})^{2}+\mu_{0}^{2}(\vartheta^{\star})^{2}-2\mu_{0}\vartheta^{\star}\gamma_{3}\big{)}$
$\displaystyle~{}+\frac{\tau_{1}t_{1}+\tau_{2}t_{2}}{2\ell}-\frac{t_{1}^{2}+t_{2}^{2}}{2\ell}+\frac{(q_{t,\tau}^{\star})^{2}}{2T_{2,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}-\frac{\eta
T_{3,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}{2},$ (10)

where the constant $\vartheta^{\star}$ satisfies $\vartheta^{\star}=0$ if
$\mu_{0}=0$ and $\vartheta^{\star}=\gamma_{3}/\mu_{0}$ otherwise, and
$T_{1}=\sqrt{\delta\mathbb{E}_{\kappa}[\kappa]}$. Here, $\gamma_{1}$,
$\gamma_{2}$ and $\gamma_{3}$ depend on the data distribution and are defined
as $\gamma_{1}=\mathbb{E}[y^{2}]$, $\gamma_{2}=\mathbb{E}[ys]$,
$\gamma_{3}=\mathbb{E}[y]$, where $y=\varphi(s)$, and $s$ is a standard
Gaussian random variable. Note that the problem defined in (IV-A) depends on
$q_{t,\tau}^{\star}$ which is given by

$\displaystyle
q_{t,\tau}^{\star}=\frac{\gamma_{2}t_{1}\widetilde{\mu}_{1}T_{1}T_{2,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}{\tau_{1}+t_{1}\widetilde{\mu}_{1}^{2}T_{1}^{2}T_{2,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}.$
(11)

Now, we summarize our main theoretical results in the following theorem.

######  Theorem 1 (Noisy Formulation Characterization).

Suppose that the assumptions in Section III are all satisfied and the cGEC
introduced in Section II is valid. Then, the training error converges in
probability as follows

$\displaystyle\mathcal{E}_{\text{train}}\xrightarrow{~{}p\to+\infty~{}}C^{\star}(\Delta,\lambda)-\frac{\lambda}{2}\left((q^{\star})^{2}+h^{\prime}(\lambda)\right),$

where $C^{\star}(\Delta,\lambda)$ is the optimal cost of the deterministic
problem in (IV-A). Here, the function $h(\cdot)$ is defined as follows

$\displaystyle
h(\lambda)=-(q^{\star})^{2}\Big{(}\lambda-\frac{1}{T_{2,\lambda}(\boldsymbol{t}^{\star},\boldsymbol{\tau}^{\star})}\Big{)}-\eta
T_{3,\lambda}(\boldsymbol{t}^{\star},\boldsymbol{\tau}^{\star}).$

Moreover, the generalization error defined in (5) converges in probability to
a deterministic function as follows

$\displaystyle{\mathcal{E}}_{\text{test}}\xrightarrow{~{}p\to+\infty~{}}\frac{1}{4^{\upsilon}}\mathbb{E}\left[\left(\varphi(g_{1})-\widehat{\varphi}(g_{2})\right)^{2}\right],$
(12)

where $g_{1}$ and $g_{2}$ have a bivariate Gaussian distribution with mean
vector $[0,\mu_{0s}\vartheta^{\star}]$ and covariance matrix $\boldsymbol{C}$,
defined as follows

$\displaystyle\boldsymbol{C}=\begin{bmatrix}1&\mu_{1s}T_{1}q^{\star}\\\
\mu_{1s}T_{1}q^{\star}&\mu_{1s}^{2}\beta^{\star}+\mu_{2s}^{2}\left((q^{\star})^{2}+h^{\prime}(\lambda)\right)\end{bmatrix}.$

The constant $\vartheta^{\star}$ satisfies $\vartheta^{\star}=0$ if
$\mu_{0}=0$ and $\vartheta^{\star}=\gamma_{3}/\mu_{0}$ otherwise. Here, the
constants $\mu_{0s}$, $\mu_{1s}$ and $\mu_{2s}$ are defined as
$\mu_{0s}=\mathbb{E}[\sigma(z)]$, $\mu_{1s}=\mathbb{E}[z\sigma(z)]$ and
$\mu_{2s}^{2}=\mathbb{E}[\sigma(z)^{2}]-\mu_{0s}^{2}-\mu_{1s}^{2}$, where $z$
is a standard Gaussian random variable. Additionally, the constant
$\beta^{\star}$ can be computed via the following expression

$\displaystyle\beta^{\star}$
$\displaystyle=\frac{1}{V_{1}+V_{3}}\Big{(}V_{1}T_{1}^{2}-V_{2}-V_{4}-\lambda+\frac{1}{T_{2,\lambda}(\boldsymbol{t}^{\star},\boldsymbol{\tau}^{\star})}\Big{)}(q^{\star})^{2}$
$\displaystyle+\frac{\eta
T_{3,\lambda}(\boldsymbol{t}^{\star},\boldsymbol{\tau}^{\star})}{V_{1}+V_{3}}-\frac{V_{2}+V_{4}+\lambda}{V_{1}+V_{3}}h^{\prime}(\lambda),$
(13)

where the constants $V_{1}$, $V_{2}$, $V_{3}$ and $V_{4}$ are defined as
follows

$\displaystyle
V_{1}=\frac{t_{1}^{\star}\widetilde{\mu}_{1}^{2}}{\tau_{1}^{\star}},~{}V_{3}=\widehat{\mu}_{1}^{2}\Big{(}\frac{t_{1}^{\star}}{\tau_{1}^{\star}\ell}+\frac{t_{2}^{\star}(\ell-1)}{\tau_{2}^{\star}\ell}\Big{)}$
$\displaystyle
V_{2}=\frac{t_{1}^{\star}\mu_{2}^{2}}{\tau_{1}^{\star}},V_{4}\textit{\textbullet}=\mu_{3}^{2}\Big{(}\frac{t_{1}^{\star}}{\tau_{1}^{\star}\ell}+\frac{t_{2}^{\star}(\ell-1)}{\tau_{2}^{\star}\ell}\Big{)}.$

Here, $q^{\star}=q^{\star}_{t^{\star},\tau^{\star}}$ is given in (11),
$\boldsymbol{t}^{\star}=[t_{1}^{\star},t_{2}^{\star}]^{\top}$ and
$\boldsymbol{\tau}^{\star}=[\tau_{1}^{\star},\tau_{2}^{\star}]^{\top}$.
Moreover,
$\\{t_{1}^{\star},t_{2}^{\star},\tau_{1}^{\star},\tau_{2}^{\star}\\}$ denotes
the optimal solution of the problem defined in (IV-A). Also, we treat
$q^{\star}$, $\boldsymbol{t}^{\star}$ and $\boldsymbol{\tau}^{\star}$ as
constants independent of $\lambda$ when we compute the derivative of the
function $h(\cdot)$.

To streamline our presentation, we postpone the proof of Theorem 1 to Section
VI. Note that Theorem 1 provides a full asymptotic characterization of the
training and generalization errors corresponding to the formulation given in
(3). Specifically, it shows that the performance of (3) can be fully
characterized after solving a deterministic scalar formulation where the cost
function depends on the parameters $\ell$ and $\Delta$. The theoretical
predictions stated in Theorem 1 are valid for any fixed noise variance
$\Delta\geq 0$ and number of noise samples $\ell\geq 1$. Additionally, it is
valid for a general family of feature matrices, activation functions and
generative models satisfying (1). The analysis presented in Section VI shows
that the deterministic problem in (IV-A) is strictly convex-concave. This
implies the uniqueness of the optimal solutions of the optimization in (IV-A).
Next, we study the properties of the noise injection method in (3) when $\ell$
grows to infinity slower than $n$, $p$ and $k$.

###  IV-B Noise Regularization Effects

Now, we consider the setting where $\ell$ grows to infinity slower than the
dimensions $n$, $p$ and $k$. We use the theoretical predictions stated in
Theorem 1 to study the regularization effects of the noise injection method in
(3). Our first theoretical result is introduced in the following theorem.

######  Theorem 2 (Regularization Effects).

Suppose that the assumptions in Theorem 1 are all satisfied. Moreover, define
the following formulation

$\displaystyle\min_{\boldsymbol{w}\in\mathbb{R}^{k}}$
$\displaystyle~{}\frac{1}{2n}\sum_{i=1}^{n}\big{(}y_{i}-\boldsymbol{w}^{\top}\widehat{\sigma}(\boldsymbol{F}^{\top}\boldsymbol{a}_{i})\big{)}^{2}+\frac{1}{2}\mathinner{\\!\left\lVert\boldsymbol{R}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}^{2}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2}.$
(14)

Here, the regularization matrix $\boldsymbol{R}$ is defined as follows

$\displaystyle\boldsymbol{R}=\widehat{\mu}_{1}^{2}\boldsymbol{F}^{\top}\boldsymbol{F}+\mu_{3}^{2}\boldsymbol{I}_{k},$
(15)

and the new activation function $\widehat{\sigma}(\cdot)$ satisfies the
properties

$\displaystyle\mathbb{E}[\widehat{\sigma}(z)]=\mathbb{E}[\sigma(x_{1})],~{}\mathbb{E}[z\widehat{\sigma}(z)]=\mathbb{E}[z\sigma(x_{1})]$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathbb{E}[\widehat{\sigma}(z)^{2}]=\mathbb{E}[\sigma(x_{1})\sigma(x_{2})],$
(16)

where $x_{1}=z+\Delta v_{1}$, $x_{2}=z+\Delta v_{2}$ and $z$, $v_{1}$ and
$v_{2}$ are independent standard Gaussian random variables. Also, define
$\widehat{\mathcal{E}}_{\text{train}}$ and
$\widehat{\mathcal{E}}_{\text{test}}$ as the training and generalization
errors corresponding to the problem in (14). Then, for any $\zeta>0$, we have
the following convergence results

$\displaystyle\begin{cases}\lim\limits_{\ell\to+\infty}\lim\limits_{p\to+\infty}\mathbb{P}\Big{(}\mathinner{\\!\left\lvert{\mathcal{E}}_{\text{train}}-\widehat{\mathcal{E}}_{\text{train}}\right\rvert}<\zeta\Big{)}=1\\\
\lim\limits_{\ell\to+\infty}\lim\limits_{p\to+\infty}\mathbb{P}\Big{(}\mathinner{\\!\left\lvert{\mathcal{E}}_{\text{test}}-\widehat{\mathcal{E}}_{\text{test}}\right\rvert}<\zeta\Big{)}=1,\end{cases}$
(17)

where ${\mathcal{E}}_{\text{test}}$ and ${\mathcal{E}}_{\text{train}}$ are the
training and generalization errors corresponding to the noisy formulation.

To streamline our presentation, we postpone the proof of Theorem 2 to Section
VI. The above theorem shows that the noisy formulation given in (3) is
equivalent to a standard formulation with a new activation function and an
additional weighted ridge regularization, when $\ell\to+\infty$. It also
provides the explicit form of the regularization. This shows that inserting
Gaussian noise during the training procedure introduces a regularization that
depend on the activation function, the feature matrix and the noise variance.
Now, we provide a precise asymptotic characterization of the formulation in
(14). Before stating our asymptotic result, we define the following
deterministic problem

$\displaystyle\inf_{\begin{subarray}{c}\tau_{1}>0\end{subarray}}\sup_{\begin{subarray}{c}t_{1}\geq
0\end{subarray}}~{}\frac{t_{1}}{2\tau_{1}}\big{(}\gamma_{1}-2\widetilde{\mu}_{1}T_{1}\widehat{q}_{t,\tau}^{\star}\gamma_{2}+\widetilde{\mu}_{1}^{2}T_{1}^{2}(\widehat{q}_{t,\tau}^{\star})^{2}+\mu_{0}^{2}(\vartheta^{\star})^{2}$
(18)
$\displaystyle-2\mu_{0}\vartheta^{\star}\gamma_{3}\big{)}+\frac{\tau_{1}t_{1}}{2}-\frac{t_{1}^{2}}{2}+\frac{(\widehat{q}_{t,\tau}^{\star})^{2}}{2\widehat{T}_{2,\lambda}(t_{1},\tau_{1})}-\frac{\eta\widehat{T}_{3,\lambda}(t_{1},\tau_{1})}{2}.$

Here, the constant $\vartheta^{\star}$ satisfies $\vartheta^{\star}=0$ if
$\mu_{0}=0$ and $\vartheta^{\star}=\gamma_{3}/\mu_{0}$ otherwise, and $T_{1}$
is defined in Section IV-A. Moreover, the functions
$\widehat{q}_{t,\tau}^{\star}$ and $\widehat{T}_{2,\lambda}(\cdot,\cdot)$ are
defined as follows

$\displaystyle\widehat{q}_{t,\tau}^{\star}=\frac{\gamma_{2}t_{1}\widetilde{\mu}_{1}T_{1}\widehat{T}_{2,\lambda}(t_{1},\tau_{1})}{\tau_{1}+t_{1}\widetilde{\mu}_{1}^{2}T_{1}^{2}\widehat{T}_{2,\lambda}(t_{1},\tau_{1})},~{}\text{and}~{}\widehat{T}_{2,\lambda}(t_{1},\tau_{1})=\frac{\delta}{T_{1}^{2}}\mathbb{E}_{\kappa}\Big{[}\frac{\kappa}{\widehat{g}_{\kappa,\lambda}(t_{1},\tau_{1})}\Big{]}/\Big{(}1-\frac{\widetilde{\mu}_{1}^{2}t_{1}\delta}{\tau_{1}}\mathbb{E}_{\kappa}\Big{[}\frac{\kappa}{\widehat{g}_{\kappa,\lambda}(t_{1},\tau_{1})}\Big{]}\Big{)}.$
(19)

Here, the functions $\widehat{T}_{3,\lambda}(\cdot,\cdot)$ and
$\widehat{g}_{\kappa,\lambda}(\cdot,\cdot)$ are defined as follows

$\displaystyle\widehat{T}_{3,\lambda}(t_{1},\tau_{1})=t_{1}^{2}\mathbb{E}_{\kappa}\Big{[}{(\widetilde{\mu}_{1}^{2}\kappa+\mu_{2}^{2}})/{\widehat{g}_{\kappa,\lambda}(t_{1},\tau_{1})}\Big{]},$
$\displaystyle\widehat{g}_{\kappa,\lambda}(t_{1},\tau_{1})=\frac{t_{1}}{\tau_{1}}\Big{(}\widetilde{\mu}_{1}^{2}\kappa+\mu_{2}^{2}\Big{)}+\Big{(}\widehat{\mu}_{1}^{2}\kappa+\mu_{3}^{2}\Big{)}+\lambda,$

where the expectations are taken over the probability distribution
$\mathbb{P}_{\kappa}(\cdot)$ defined in Assumption 5. Now, we summarize the
asymptotic properties of the limiting formulation in (14) in the following
theorem.

######  Lemma 1 (Limiting Formulation Characterization).

Suppose that the assumptions in Theorem 1 are all satisfied. Then, the
training error corresponding to the limiting formulation in (14) converges in
probability as follows

$\displaystyle\widehat{\mathcal{E}}_{\text{train}}\xrightarrow{~{}p\to+\infty~{}}\widehat{C}^{\star}(\Delta,\lambda)-\frac{\lambda}{2}\left((\widehat{q}^{\star})^{2}+\widehat{h}^{\prime}(\lambda)\right),$

where $\widehat{C}^{\star}(\Delta,\lambda)$ is the optimal cost of the
deterministic problem in (18). Here, the function $\widehat{h}(\cdot)$ is
defined as follows

$\displaystyle\widehat{h}(\lambda)=-(\widehat{q}^{\star})^{2}\Big{(}\lambda-\frac{1}{\widehat{T}_{2,\lambda}(t_{1}^{\star},\tau_{1}^{\star})}\Big{)}-\eta\widehat{T}_{3,\lambda}(t_{1}^{\star},\tau_{1}^{\star}).$

Moreover, the generalization error corresponding to the limiting formulation
in (14) converges in probability to a deterministic function as follows

$\displaystyle\widehat{\mathcal{E}}_{\text{test}}\xrightarrow{~{}p\to+\infty~{}}\frac{1}{4^{\upsilon}}\mathbb{E}\left[\left(\varphi(g_{1})-\widehat{\varphi}(g_{2})\right)^{2}\right],$
(20)

where $g_{1}$ and $g_{2}$ have a bivariate Gaussian distribution with mean
vector $[0,\mu_{0s}\vartheta^{\star}]$ and covariance matrix $\boldsymbol{C}$,
defined as follows

$\displaystyle\boldsymbol{C}=\begin{bmatrix}1&\mu_{1s}T_{1}\widehat{q}^{\star}\\\
\mu_{1s}T_{1}\widehat{q}^{\star}&\mu_{1s}^{2}\widehat{\beta}^{\star}+\mu_{2s}^{2}\left((\widehat{q}^{\star})^{2}+\widehat{h}^{\prime}(\lambda)\right)\end{bmatrix}.$

The constant $\vartheta^{\star}$ satisfies $\vartheta^{\star}=0$ if
$\mu_{0}=0$ and $\vartheta^{\star}=\gamma_{3}/\mu_{0}$ otherwise. Here, the
constants $\mu_{0s}$, $\mu_{1s}$ and $\mu_{2s}$ are defined as
$\mu_{0s}=\mathbb{E}[\sigma(z)]$, $\mu_{1s}=\mathbb{E}[z\sigma(z)]$ and
$\mu_{2s}^{2}=\mathbb{E}[\sigma(z)^{2}]-\mu_{0s}^{2}-\mu_{1s}^{2}$, where $z$
is a standard Gaussian random variable. Additionally, the constant
$\widehat{\beta}^{\star}$ can be computed via the following expression

$\displaystyle\widehat{\beta}^{\star}$
$\displaystyle=\frac{1}{V_{1}+V_{3}}\Big{(}V_{1}T_{1}^{2}-V_{2}-V_{4}-\lambda+\frac{1}{\widehat{T}_{2,\lambda}(t_{1}^{\star},\tau_{1}^{\star})}\Big{)}(\widehat{q}^{\star})^{2}$
$\displaystyle+\frac{\eta\widehat{T}_{3,\lambda}(t_{1}^{\star},\tau_{1}^{\star})}{V_{1}+V_{3}}-\frac{V_{2}+V_{4}+\lambda}{V_{1}+V_{3}}\widehat{h}^{\prime}(\lambda),$
(21)

where the constants $V_{1}$, $V_{2}$, $V_{3}$ and $V_{4}$ are defined as
follows

$\displaystyle
V_{1}=\frac{t_{1}^{\star}\widetilde{\mu}_{1}^{2}}{\tau_{1}^{\star}},~{}V_{3}=\widehat{\mu}_{1}^{2},~{}V_{2}=\frac{t_{1}^{\star}\mu_{2}^{2}}{\tau_{1}^{\star}},V_{4}=\mu_{3}^{2}.$

Here, $\widehat{q}^{\star}=\widehat{q}^{\star}_{t^{\star},\tau^{\star}}$ is
given in (19). Moreover, $\\{t_{1}^{\star},\tau_{1}^{\star}\\}$ denotes the
optimal solution of the problem defined in (18). Also, we treat
$\widehat{q}^{\star}$, $t_{1}^{\star}$ and $\tau_{1}^{\star}$ as constants
independent of $\lambda$ when we compute the derivative of the function
$\widehat{h}(\cdot)$.

The proof of Lemma 1 is provided in Section VI. The results in Theorem 2 and
Lemma 1 are based on the asymptotic predictions stated in Theorem 1.
Specifically, we show in Section VI that the asymptotic problem corresponding
to the noisy formulation in (IV-A) converges to the deterministic problem in
(18), when $\ell$ grows to infinity. Then, we show that the deterministic
problem in (18) is the asymptotic limit of the formulation in (14) using the
CGMT framework. The analysis presented in Section VI shows that the
deterministic problem in (18) is strictly convex-concave. This implies the
uniqueness of its optimal solutions.

##  V Simulation Results

In this part, we provide additional simulation examples to verify our
asymptotic results stated in Theorem 1, Theorem 2 and Lemma 1. Our predictions
stated in Section IV are valid for a general family of feature matrices,
activation functions and generative models satisfying (1). We specialize our
general results to popular learning models.

In particular, we consider two families of feature matrices. We consider
feature matrices that can be expressed as $\boldsymbol{F}=d\boldsymbol{V}$,
where: (a) The scalar $d$ satisfies $d=1/\sqrt{p}$ and the matrix
$\boldsymbol{V}$ has independent standard Gaussian components. We refer to
this matrix as the Gaussian feature matrix. (b) The scalar $d$ satisfies
$d=\sqrt{3/p}$ and the matrix $\boldsymbol{V}$ has independent uniformly
distributed components in $[-1~{}1]$. We refer to this matrix as the uniform
feature matrix.

Also, we consider two popular regression and classification models. For the
regression model, we assume that $\varphi(\cdot)$ is the ReLu function and
$\widehat{\varphi}(\cdot)$ is the identity function. For the classification
model, we assume that $\varphi(\cdot)$ is the sign function with possible sign
flip with probability $\theta$ and $\widehat{\varphi}(\cdot)$ is the sign
function.

###  V-A Correlated Gaussian Equivalence

In this simulation example, we consider a particular case of the cGEC to
illustrate the Gaussian equivalence. Here, we consider a projected version of
the probability density function $\mathbb{P}(\nu_{1},\nu_{2},\nu_{3})$,
defined in Section II. Specifically, we consider the probability density
function of the random variable $\nu=\nu_{1}+\nu_{2}+\nu_{3}$.

Figure 3: Original: measures the distribution of
$\nu=\nu_{1}+\nu_{2}+\nu_{3}$. Gaussian: measures the distribution of
$\nu=\nu_{g,1}+\nu_{g,2}+\nu_{g,3}$. Gaussian density: denotes the Gaussian
pdf with the corresponding first two moments of $\nu$. (a) $\boldsymbol{F}$ is
the Gaussian feature matrix and $\sigma(\cdot)$ is the tanh function. (b)
$\boldsymbol{F}$ is the uniform feature matrix in $[-\sqrt{3}~{}\sqrt{3}]$ and
$\sigma(\cdot)$ is the ReLu function. We take $p=10^{4}$, $\alpha=1.1$,
$\Delta=1$ and $\eta=1.3$.

Figure 3 considers two different feature matrices and activation functions. It
also compares the probability density function of the random variable
$\nu=\nu_{1}+\nu_{2}+\nu_{3}$ with the probability density function of the
random variable $\nu_{g}=\nu_{g,1}+\nu_{g,2}+\nu_{g,3}$. Figure 3 shows that
the random variable $\nu$ is Gaussian with the same first two moments of
$\nu_{g}$ for particular choices of $\boldsymbol{F}$, $\sigma(\cdot)$,
$\boldsymbol{\xi}$ and $\boldsymbol{w}$. This provides a particular
illustration of the cGEC.

###  V-B Limiting Performance

Our fourth simulation considers the non–linear regression model. Figure 4
compares the numerical predictions and our predictions stated in Theorem 2 and
Lemma 1.

Figure 4: Solid line: Theoretical predictions. Circle: Numerical simulations
for (3) in 4 and for both (3) and (7) in 4. Black cross: Numerical simulations
for (7). (a) $p=500$, $\Delta=0.4$, $\alpha=1.5$ and $\lambda=10^{-2}$. (b)
$p=500$, $\Delta=0.6$, $\alpha=2$, $\lambda=10^{-4}$, $\eta=1$ and
$\sigma(\cdot)$ is the SoftPlus. Binary denotes the binary step activation.
$\boldsymbol{F}$ is the Gaussian feature matrix. The number of Monte Carlo
trials is $100$.

This simulation example first provides an empirical verification of the
theoretical predictions in Theorem 2 and Lemma 1. It particularly shows that
our predictions are in excellent agreement with the empirical results for (3)
and (7). Furthermore, note that the performance of the deterministic
formulation given in (18) is achieved with a moderate number of noise samples,
i.e. $\ell=70$ and $\ell=100$. This verifies the results stated in Theorem 2
and Lemma 1 and provides an empirical verification of the cGEC introduced in
Section II. Figure 4 further shows that the considered noisy formulation can
asymptotically mitigate the double descent in the generalization error for an
appropriately selected activation function and fixed noise variance.
Specifically, note that the ReLu and binary activation functions lead to a
decreasing generalization performance which is not the case for the SoftPlus
activation. Figure 4 illustrates the convergence behavior of the
generalization error corresponding to (3) for the SoftPlus activation and
fixed $\eta$. It particularly shows that the generalization error of (3)
converges to the generalization error of (7) when $\ell$ grows to infinity.
Moreover, note that the limit is already achieved with a small value of
$\ell$. This verifies the predictions in Theorem 2.

###  V-C Impact of the Noise Variance

In this simulation example, we study the effects of the noise variance
$\Delta$ on the generalization error corresponding to the noisy formulation
and the limiting formulation. Here, we consider the binary classification
model. Figure 5 compares the numerical predictions and our theoretical
predictions stated in Theorem 1, Theorem 2 and Lemma 1.

Figure 5: Solid line: Theoretical predictions. Circle: Numerical simulations
for (3) in 5 and 5 and for (7) in 5. Black cross: Numerical simulations for
((C.1)). (a) $\boldsymbol{F}$ is the Gaussian feature matrix and
$\sigma(\cdot)$ is the tanh activation function. We set $p=400$, $\ell=50$,
$\alpha=2$, $\theta=0.1$ and $\lambda=10^{-5}$. (b) $\boldsymbol{F}$ is the
uniform feature matrix and $\sigma(\cdot)$ is the SoftPlus activation. We set
$p=1500$, $\alpha=1.5$, $\theta=0$ and $\lambda=10^{-4}$. (c) $\boldsymbol{F}$
is the uniform feature matrix. We set $p=600$, $\alpha=2$, $\theta=0$,
$\ell=20$, $\eta=1.5$ and $\lambda=10^{-4}$. The number of Monte Carlo trials
is $100$.

It provides another empirical verification of our theoretical predictions
since our results are in excellent agreement with the actual performance of
the considered formulations. It also provides an empirical verification of the
cGEC discussed in Section II. Figure 5 studies the effects of the noise
variance $\Delta$ on the generalization error corresponding to the noisy
formulation for fixed $\ell$. Note that increasing the noise variance improves
the generalization error especially at low $\eta$. Figure 5 also suggests that
an optimized noise variance can reduce the effects of the double descent
phenomenon. Now, Figure 5 studies the effects of the noise variance $\Delta$
on the generalization error corresponding to the limiting formulation. We can
see that the generalization error increases after reaching a minimum for
$\Delta=0.5$. For $\Delta=4$, observe that the generalization error is
deceasing. This suggests that the double descent phenomenon can be mitigated
for an appropriately selected noise variance. Figure 5 also shows that the
generalization error corresponding to the noisy formulation has a unique
minimum as a function of $\Delta$ for the considered activation functions.
This suggests that the optimization problem over the noise variance $\Delta$
has a unique solution. This can simplify the design of an efficient
optimization scheme of the generalization error in terms of $\Delta$.

###  V-D Alternative Formulations

Now, we consider the binary classification model, where $\theta=0$. We compare
the performance of the noisy formulation given in (3) and the dropout
technique. In this paper, we consider the following version of the dropout
method

$\displaystyle\min_{\boldsymbol{w}\in\mathbb{R}^{k}}$
$\displaystyle\frac{1}{2n\ell}\sum_{i=1}^{n}\sum_{j=1}^{\ell}\big{(}y_{i}-\boldsymbol{w}^{\top}\sigma\big{(}\boldsymbol{D}_{ij}\boldsymbol{F}^{\top}\boldsymbol{a}_{i}\big{)}\big{)}^{2}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2},$

where $\\{\boldsymbol{D}_{ij}\\}_{i,j}^{n,\ell}$ are diagonal matrices with
independent and identically distributed diagonal entries drawn from the
distribution, $\mathbb{P}(d=1)=1-\epsilon$ and $\mathbb{P}(d=0)=\epsilon$,
where $\epsilon$ denotes the probability of dropping a unit. The above
formulation is similar to the one considered in [1, 30].

####  V-D1 Performance Comparison

In Figure 6, we compare the performance of the noisy formulation and the
dropout formulation for four different activation functions.

Figure 6: Solid line: Theoretical predictions. Circle: Numerical simulations
for (3) and (7). Hexagram: Numerical simulations for the dropout formulation.
(a) Sigmoid activation and $p=600$, $\alpha=1.4$, $\Delta=1.5$,
$\lambda=10^{-4}$, $\ell=60$ and $\epsilon=0.3$. (b) Erf activation and
$p=600$, $\alpha=1.4$, $\Delta=2$, $\lambda=10^{-3}$, $\ell=40$ and
$\epsilon=0.3$. (c) ReLu activation and $p=600$, $\alpha=1.4$, $\Delta=3.5$,
$\lambda=10^{-3}$, $\ell=60$ and $\epsilon=0.3$. (d) Sign activation and
$p=600$, $\alpha=1.6$, $\Delta=1.5$, $\lambda=10^{-4}$, $\ell=60$ and
$\epsilon=0.4$. $\boldsymbol{F}$ is the uniform feature matrix. The results
are averaged over $35$ independent Monte Carlo trials.

First, we can notice that our asymptotic results provided in Theorem 1,
Theorem 2 and Lemma 1 match with the actual performance of the noisy
formulation and its limiting formulation. This gives another empirical
verification of our theoretical predictions. Figure 6 suggests that the noisy
and dropout formulations significantly improve the generalization performance
of the standard formulation for an appropriately selected activation function,
$\Delta$, $\ell$ and $\epsilon$. Figure 6 shows that the noisy and dropout
formulations have a similar generalization performance for fixed $\ell=60$ and
for the sigmoid activation function. Moreover, it shows that the dropout
formulation approaches the generalization performance of the limiting
formulation at high $\eta$. Moreover, Figures 6 and 6 shows that the dropout
method provides a largely better performance as compared to the noisy
formulation for fixed $\ell=40$ and $\ell=60$ and for the Erf and the ReLu
activation functions, respectively. It also suggests that the limiting and
dropout formulations have a similar generalization performance where the
dropout method is better at high $\eta$ for the ReLu activation function. Now,
Figure 6 considers the sign activation function and shows that the noisy
formulation provides a better generalization performance as compared to the
dropout formulation at low $\eta$ and for a fixed number of noise injections,
$\ell=60$. We can also see that the limiting formulation generalizes better
than the dropout method for the considered parameters. This simulation example
particularly suggests that the performance of the noisy and dropout
formulations depends on the activation function. Moreover, the dropout
formulation have a similar generalization performance as compared to the
limiting performance for the considered parameters.

####  V-D2 Convergence Behavior

In the last simulation example, we study the convergence behavior of the noisy
and dropout formulations for different activation functions. Figure 7 first
shows that our theoretical predictions stated in Theorem 1, Theorem 2 and
Lemma 1 match with the actual performance of the noisy formulation and its
limiting formulation. This gives another empirical verification of our
predictions.

Figure 7: Solid line: Theoretical predictions. Circle: Numerical simulations
for (3) and (7). Hexagram: Numerical simulations for the dropout formulation.
(a) Sigmoid activation function and we set the parameters as $p=500$,
$\alpha=1.4$, $\Delta=1.5$, $\lambda=10^{-4}$, $\eta=2.5$ and $\epsilon=0.3$.
(b) Erf activation function and we set $p=600$, $\alpha=1.4$, $\Delta=2$,
$\lambda=10^{-3}$, $\eta=2$ and $\epsilon=0.3$. (c) ReLu activation function
and we set the parameters as $p=500$, $\alpha=1.4$, $\Delta=3.5$,
$\lambda=10^{-3}$, $\eta=1.5$ and $\epsilon=0.3$. (d) Sign activation function
and we set the parameters as $p=500$, $\alpha=1.6$, $\Delta=1.5$,
$\lambda=10^{-4}$, $\eta=2$ and $\epsilon=0.4$. The feature matrix
$\boldsymbol{F}$ is the uniform feature matrix. The results are averaged over
$35$ independent Monte Carlo trials.

Figure 7 studies the convergence properties of both approaches as a function
of $\ell$ and for different activation functions. It particularly suggests
that the dropout method has a better convergence rate as compared to the noisy
formulation. Moreover, Figure 7 suggests that the noisy and dropout
formulations have comparable generalization performance when the number of
noise injections grows to infinity.

##  VI Technical Details: Analysis of the Noisy Formulation

In this part, we provide a rigorous proof of the predictions stated in Theorem
1, Theorem 2 and Lemma 1. To this end, we suppose that the assumptions
considered in Sections II and III are all satisfied. We derive our theoretical
results using an extended version of the CGMT framework which we refer to as
the multivariate CGMT.

###  VI-A Multivariate Convex Gaussian Min–Max Theorem

To derive the asymptotic results stated in Theorem 1, Theorem 2 and Lemma 1,
we use an extended version of the CGMT framework introduced in [10]. The CGMT
is used to accurately analyze a generally hard primary formulation by
introducing an asymptotically equivalent auxiliary optimization problem. In
this paper, we consider primary optimization problems of the following form

$\Phi_{k}=\min\limits_{\boldsymbol{w}\in\mathcal{S}_{\boldsymbol{w}}}\max\limits_{\boldsymbol{u}\in\mathcal{S}_{\boldsymbol{u}}}\sum_{i=1}^{\ell}\boldsymbol{u}_{i}\boldsymbol{G}_{i}\boldsymbol{w}_{i}+\psi(\boldsymbol{w},\boldsymbol{u}),$
(22)

where $\boldsymbol{u}_{i}\in\mathbb{R}^{n_{i}}$ and
$\boldsymbol{w}_{i}\in\mathbb{R}^{k_{i}}$ are optimization variables and
$\boldsymbol{G}_{i}\in\mathbb{R}^{n_{i}\times k_{i}}$ has independent standard
Gaussian random components, for any $i\in\\{1\dots,\ell\\}$. Additionally, the
vectors $\boldsymbol{w}$ and $\boldsymbol{u}$ are formed by the concatenation
of the vectors $\\{\boldsymbol{w}_{i}\\}_{i=1}^{\ell}$ and
$\\{\boldsymbol{u}_{i}\\}_{i=1}^{\ell}$, respectively. We refer to the
formulation in (22) as the multivariate primary optimization (multivariate
PO). We show that the corresponding multivariate auxiliary optimization
(multivariate AO) is given by

$\displaystyle\phi_{k}=\min\limits_{\boldsymbol{w}\in\mathcal{S}_{\boldsymbol{w}}}\max\limits_{\boldsymbol{u}\in\mathcal{S}_{\boldsymbol{u}}}$
$\displaystyle\sum_{i=1}^{\ell}\mathinner{\\!\left\lVert\boldsymbol{u}_{i}\right\rVert}\boldsymbol{g}_{i}^{\top}\boldsymbol{w}_{i}+\sum_{i=1}^{\ell}\mathinner{\\!\left\lVert\boldsymbol{w}_{i}\right\rVert}\boldsymbol{h}_{i}^{\top}\boldsymbol{u}_{i}+\psi(\boldsymbol{w},\boldsymbol{u}),$
(23)

where $\boldsymbol{g}_{i}\in\mathbb{R}^{k_{i}}$ and
$\boldsymbol{h}_{i}\in\mathbb{R}^{n_{i}}$ are independent standard Gaussian
random vectors, for any $i\in\\{1\dots,\ell\\}$. Here, we assume that
$\boldsymbol{G}_{i}\in\mathbb{R}^{{n_{i}}\times{k_{i}}}$,
$\boldsymbol{g}_{i}\in\mathbb{R}^{{k_{i}}}$ and
$\boldsymbol{h}_{i}\in\mathbb{R}^{{n_{i}}}$, are all independent, the
feasibility sets $\mathcal{S}_{\boldsymbol{w}}\subset\mathbb{R}^{{k}}$ and
$\mathcal{S}_{\boldsymbol{u}}\subset\mathbb{R}^{{n}}$ are convex and compact,
and the function $\psi\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{{k}}\times\mathbb{R}^{{n}}\to\mathbb{R}$ is continuous
_convex-concave_ on
$\mathcal{S}_{\boldsymbol{w}}\times\mathcal{S}_{\boldsymbol{u}}$, where
$k=\sum_{i=1}^{\ell}k_{i}$ and $n=\sum_{i=1}^{\ell}n_{i}$. Now, we summarize
our theoretical result in the following theorem.

######  Theorem 3 (Multivariate CGMT).

Assume that the above assumptions are all satisfied. For any fixed $\ell\geq
1$ and $\epsilon>0$, consider an open set $\mathcal{S}_{p,\epsilon}$.
Moreover, define the set
$\mathcal{S}^{c}_{p,\epsilon}={\mathcal{S}}_{\boldsymbol{w}}\setminus\mathcal{S}_{p,\epsilon}$.
Let $\phi_{k}$ and $\phi^{c}_{k}$ be the optimal cost values of the
multivariate AO formulation in (23) with feasibility sets
${\mathcal{S}}_{\boldsymbol{w}}$ and $\mathcal{S}^{c}_{p,\epsilon}$,
respectively. Assume that the following properties are all satisfied

  * (1)

There exists a constant $\phi$ such that the optimal cost $\phi_{k}$ converges
in probability to $\phi$ as $k$ goes to $+\infty$.

  * (2)

There exists a constant $\phi^{c}$ such that the optimal cost $\phi^{c}_{k}$
converges in probability to $\phi^{c}$ as $k$ goes to $+\infty$, for any fixed
$\epsilon>0$.

  * (3)

There exists a positive constant $\zeta>0$ such that $\phi^{c}\geq\phi+\zeta$,
for any fixed $\epsilon>0$.

Then, the following convergence in probability holds

$\mathinner{\\!\left\lvert\Phi_{k}-\phi_{k}\right\rvert}\overset{k\to+\infty}{\longrightarrow}0,~{}\text{and}~{}\mathbb{P}(\widehat{\boldsymbol{w}}_{k}\in\mathcal{S}_{p,\epsilon})\overset{k\to\infty}{\longrightarrow}1,$

for any fixed $\epsilon>0$, where $\Phi_{k}$ and
$\widehat{\boldsymbol{w}}_{k}$ are the optimal cost and the optimal solution
of the multivariate PO formulation in (22).

The above theorem allows us to analyze the generally easy multivariate AO
formulation given in (23) to infer asymptotic properties of the generally hard
multivariate PO formulation in (22). The proof of Theorem 3 follows by showing
that the formulation in (23) and the following formulation

$\displaystyle\widehat{\Phi}_{k}$
$\displaystyle=\min\limits_{\boldsymbol{w}\in\mathcal{S}_{\boldsymbol{w}}}\max\limits_{\boldsymbol{u}\in\mathcal{S}_{\boldsymbol{u}}}\sum_{i=1}^{\ell}\boldsymbol{u}_{i}\boldsymbol{G}_{i}\boldsymbol{w}_{i}+\psi(\boldsymbol{w},\boldsymbol{u})+\sum_{i=1}^{\ell}\mathinner{\\!\left\lVert\boldsymbol{u}_{i}\right\rVert}\mathinner{\\!\left\lVert\boldsymbol{w}_{i}\right\rVert}g_{i},$

satisfy all the assumptions in [21], where $\\{g_{i}\\}_{i=1}^{\ell}$ are
independent standard Gaussian random variables. Then, following the same
analysis in [11] and [10], we can show that for any $\chi\in\mathbb{R}$ and
$c>0$, it holds

$\mathbb{P}\left(\mathinner{\\!\left\lvert\Phi_{k}-\chi\right\rvert}>c\right)\leq
2^{\ell}~{}\mathbb{P}\left(\mathinner{\\!\left\lvert\phi_{k}-\chi\right\rvert}>c\right).$
(24)

Combining this result with the assumptions of Theorem 3 completes the proof.
We omit the detailed proof since it is similar to the analysis in [11] and
[10]. We refer to Theorem 3 as the multivariate convex Gaussian min-max
theorem (multivariate CGMT).

Next, we use the multivariate CGMT to rigorously prove the technical results
provided in Theorem 1, Theorem 2 and Lemma 1. Our approach is to reformulate
the Gaussian formulation in ((C.1)) in the form of the multivariate PO problem
given in (22). Then, use the multivariate CGMT framework to show that the
formulation in (3) is asymptotically equivalent to an easier formulation that
can be written in the form of the multivariate AO problem given in (23). The
next step is to show that the multivariate AO formulation converges in
probability to a deterministic problem that can be expressed in the form of
the formulation given in (IV-A).

###  VI-B Asymptotic Analysis of the Noisy Formulation

In this part, we provide the technical steps to obtain the theoretical results
stated in Theorem 1. Specifically, we use the multivariate CGMT framework to
precisely analyze the noisy formulation introduced in (3). Next, we suppose
that the assumptions introduced in Section III are all satisfied.

####  VI-B1 Formulating the Multivariate Primary Formulation

Based on the cGEC introduced in Section II, it suffices to precisely analyze
the Gaussian formulation in the large system limit. Then, it suffices to
analyze the following formulation

$\displaystyle\min_{\boldsymbol{w}\in\mathbb{R}^{k}}$
$\displaystyle\frac{1}{2n\ell}\sum_{i=1}^{n}\sum_{j=1}^{\ell}\big{(}y_{i}-\widetilde{\mu}_{1}\boldsymbol{w}^{\top}\boldsymbol{F}^{\top}\boldsymbol{a}_{i}-\widehat{\mu}_{1}\boldsymbol{w}^{\top}\boldsymbol{F}^{\top}\boldsymbol{z}_{ij}$
$\displaystyle-\mu_{0}\boldsymbol{w}^{\top}\boldsymbol{1}_{k}-\mu_{2}\boldsymbol{w}^{\top}\boldsymbol{b}_{i}-\mu_{3}\boldsymbol{w}^{\top}\boldsymbol{p}_{ij}\big{)}^{2}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2}.$
(25)

Note that the formulation in (VI-B1) is strongly convex with a strong
convexity parameter equals to $\lambda$. This means that it has a unique
optimal solution. Note that the multivariate CGMT framework assumes that the
feasibility sets of the multivariate PO formulation in (22) are compact. The
following lemma shows that this assumption is satisfied by our formulation.

######  Lemma 2 (Primal Compactness).

Assume that $\widehat{\boldsymbol{w}}_{p}$ is the unique optimal solution of
the formulation given in (VI-B1). Then, there exist two positive constants
$C_{w}>0$ and $C_{\vartheta}>0$ such that

$\displaystyle\mathbb{P}\big{(}\mathinner{\\!\left\lVert\widehat{\boldsymbol{w}}_{p}\right\rVert}\leq
C_{w}\big{)}\overset{p\to\infty}{\longrightarrow}1,~{}\mathbb{P}\big{(}\mathinner{\\!\left\lvert\boldsymbol{1}_{k}^{\top}\widehat{\boldsymbol{w}}_{p}\right\rvert}\leq
C_{\vartheta}\big{)}\overset{p\to\infty}{\longrightarrow}1,$ (26)

where the second asymptotic result is valid only when $\mu_{0}\neq 0$.

Given that the loss function in (VI-B1) is proper and strongly convex, one can
use the results in [16, Lemma 1] to prove Lemma 2. This asymptotic result
follows using Assumptions 3, 4, and 5 and [34, Theorem 2.1]. Combining this
result with the theoretical result stated in [16, Proposition 1], the Gaussian
formulation is asymptotically equivalent to the following formulation

$\displaystyle\min_{\begin{subarray}{c}c_{w}\leq\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}\leq
C_{w}\\\ \mathinner{\\!\left\lvert\vartheta\right\rvert}\leq
C_{\vartheta}\end{subarray}}$
$\displaystyle\frac{1}{2n\ell}\sum_{j=1}^{\ell}\mathinner{\\!\left\lVert\boldsymbol{y}-[\widetilde{\mu}_{1}\boldsymbol{A}+\widehat{\mu}_{1}\boldsymbol{Z}_{j}]\boldsymbol{F}\boldsymbol{w}-\mu_{0}\vartheta\boldsymbol{1}_{n}-\mu_{2}\boldsymbol{B}\boldsymbol{w}-\mu_{3}\boldsymbol{P}_{j}\boldsymbol{w}\right\rVert}^{2}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2},$
(27)

where the matrices $\boldsymbol{A}\in\mathbb{R}^{n\times p}$,
$\boldsymbol{Z}_{j}\in\mathbb{R}^{n\times p}$,
$\boldsymbol{B}\in\mathbb{R}^{n\times k}$ and
$\boldsymbol{P}_{j}\in\mathbb{R}^{n\times k}$ are formed by the concatenation
of the vectors $\boldsymbol{a}_{i}$, $\boldsymbol{z}_{ij}$,
$\boldsymbol{b}_{i}$ and $\boldsymbol{p}_{ij}$, respectively. Here, the label
vector $\boldsymbol{y}\in\mathbb{R}^{n}$ is formed by the labels
$\\{y_{i}\\}_{i=1}^{n}$ and can be expressed as follows
$\boldsymbol{y}=\varphi(\boldsymbol{A}\boldsymbol{\xi})$. Based on the
analysis in [16, Proposition 2], it suffices to precisely analyze the problem
in (27) for fixed feasible $\vartheta$. Then, minimize its asymptotic limit
over $\vartheta$ to infer the asymptotic properties of (27). Next, we start by
analyzing the formulation in (27) for fixed feasible $\vartheta$. To express
(27) in the form of the multivariate PO introduced in (22), we introduce
additional dual optimization variables. Specifically, the formulation in (27)
can be equivalently formulated as follows

$\displaystyle\min_{\begin{subarray}{c}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}\leq
C_{w}\end{subarray}}\max_{\boldsymbol{u}\in\mathbb{R}^{\ell n}}$
$\displaystyle-\frac{\mathinner{\\!\left\lVert\boldsymbol{u}\right\rVert}^{2}}{2n\ell}+\frac{1}{n\ell}\sum_{j=1}^{\ell}\boldsymbol{u}_{j}^{\top}\Big{(}\mu_{0}\vartheta\boldsymbol{1}_{n}+\mu_{2}\boldsymbol{B}\boldsymbol{w}-\boldsymbol{y}$
(28)
$\displaystyle+[\widetilde{\mu}_{1}\boldsymbol{A}+\widehat{\mu}_{1}\boldsymbol{Z}_{j}]\boldsymbol{F}\boldsymbol{w}+\mu_{3}\boldsymbol{P}_{j}\boldsymbol{w}\Big{)}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2},$

where the dual optimization vector $\boldsymbol{u}\in\mathbb{R}^{\ell n}$ can
be decomposed as follows
$\boldsymbol{u}^{\top}=[\boldsymbol{u}_{1}^{\top},\dots,\boldsymbol{u}_{\ell}^{\top}]^{\top}$,
where $\boldsymbol{u}_{i}\in\mathbb{R}^{n}$, for any $1\leq i\leq\ell$. Note
that the optimization problem given in (28) has a unique optimal solution. The
multivariate CGMT also assumes that the feasibility set of the maximization
problem in (28) is compact. The following lemma shows that this assumption is
also satisfied by our formulation.

######  Lemma 3 (Dual Compactness).

Assume that $\widehat{\boldsymbol{u}}_{p}$ is the unique optimal solution of
the formulation given in (28). Then, there exists a positive constants
$C_{u}>0$ such that

$\displaystyle\mathbb{P}\big{(}\mathinner{\\!\left\lVert\widehat{\boldsymbol{u}}_{p}\right\rVert}/\sqrt{n}\leq
C_{u}\big{)}\overset{p\to\infty}{\longrightarrow}1.$ (29)

This result can also be proved using similar steps as in [16, Lemma 2].
Specifically, we can use the result in [35, Proposition 11.3] to show the
compactness of the optimal dual vector $\widehat{\boldsymbol{u}}_{p}$. The
results in Lemma 2 and Lemma 3 show that the Gaussian formulation is
asymptotically equivalent to the following formulation

$\displaystyle\min_{\begin{subarray}{c}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}\leq
C_{w}\end{subarray}}\max_{\frac{\mathinner{\\!\left\lVert\boldsymbol{u}\right\rVert}}{\sqrt{n}}\leq
C_{u}}$
$\displaystyle-\frac{\mathinner{\\!\left\lVert\boldsymbol{u}\right\rVert}^{2}}{2n\ell}+\frac{1}{n\ell}\sum_{j=1}^{\ell}\boldsymbol{u}_{j}^{\top}\big{(}\mu_{0}\vartheta\boldsymbol{1}_{n}+\mu_{2}\boldsymbol{B}\boldsymbol{w}-\boldsymbol{y}$
(30)
$\displaystyle+[\widetilde{\mu}_{1}\boldsymbol{A}+\widehat{\mu}_{1}\boldsymbol{Z}_{j}]\boldsymbol{F}\boldsymbol{w}+\mu_{3}\boldsymbol{P}_{j}\boldsymbol{w}\big{)}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2}.$

Next, we focus on precisely analyzing the formulation in (30). Now, note that
the label vector $\boldsymbol{y}$ depend on the Gaussian matrix
$\boldsymbol{A}$. Then, we decompose $\boldsymbol{A}$ as follows

$\displaystyle\boldsymbol{A}=\boldsymbol{A}\boldsymbol{P}_{\boldsymbol{\xi}}+\boldsymbol{A}\boldsymbol{P}^{\perp}_{\boldsymbol{\xi}}=\boldsymbol{A}{\boldsymbol{\xi}}{\boldsymbol{\xi}}^{\top}+\boldsymbol{A}\boldsymbol{P}^{\perp}_{\boldsymbol{\xi}},$
(31)

where $\boldsymbol{P}_{\boldsymbol{\xi}}\in\mathbb{R}^{p\times p}$ denotes the
projection matrix onto the space spanned by the vector
$\boldsymbol{\xi}\in\mathbb{R}^{p}$ and
$\boldsymbol{P}^{\perp}_{\boldsymbol{\xi}}=\boldsymbol{I}_{p}-{\boldsymbol{\xi}}{\boldsymbol{\xi}}^{\top}$
denotes the projection matrix onto the orthogonal complement of the space
spanned by the vector $\boldsymbol{\xi}$. Note that the random matrix
$\boldsymbol{A}{\boldsymbol{\xi}}{\boldsymbol{\xi}}^{\top}$ is independent of
the random matrix $\boldsymbol{A}\boldsymbol{P}^{\perp}_{\boldsymbol{\xi}}$.
Then, we can express $\boldsymbol{A}$ as follows without changing its
statistics

$\displaystyle\boldsymbol{A}=\boldsymbol{s}{\boldsymbol{\xi}}^{\top}+\boldsymbol{A}\boldsymbol{P}^{\perp}_{\boldsymbol{\xi}},$
(32)

where $\boldsymbol{s}\in\mathbb{R}^{n}$ has independent standard Gaussian
components and the two random quantities $\boldsymbol{s}$ and $\boldsymbol{A}$
are independent. This shows that the optimization problem formulated in (30)
is statistically equivalent to the following formulation

$\displaystyle\min_{\begin{subarray}{c}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}\leq
C_{w}\end{subarray}}\max_{\frac{\mathinner{\\!\left\lVert\boldsymbol{u}\right\rVert}}{\sqrt{n}}\leq
C_{u}}$
$\displaystyle-\frac{\mathinner{\\!\left\lVert\boldsymbol{u}\right\rVert}^{2}}{2n\ell}+\frac{1}{n\ell}\sum_{j=1}^{\ell}\boldsymbol{u}_{j}^{\top}\big{(}-\boldsymbol{y}+\widetilde{\mu}_{1}\boldsymbol{s}{\boldsymbol{\xi}}^{\top}\boldsymbol{F}\boldsymbol{w}$
(33)
$\displaystyle+\mu_{0}\vartheta\boldsymbol{1}_{n}+\boldsymbol{G}\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{w}+\boldsymbol{T}_{j}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}\big{)}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2},$

where $\boldsymbol{G}$ and $\\{\boldsymbol{T}_{j}\\}_{j=1}^{\ell}$ are
independent matrices with independent and identically distributed standard
Gaussian components. Here, $\boldsymbol{G}\in\mathbb{R}^{n\times k}$,
$\boldsymbol{T}_{j}\in\mathbb{R}^{n\times k}$. Moreover,
$\boldsymbol{\Sigma}\in\mathbb{R}^{k\times k}$ and
$\boldsymbol{\Gamma}\in\mathbb{R}^{k\times k}$ are positive definite matrices
and defined as follows

$\displaystyle\begin{cases}\boldsymbol{\Sigma}=\widetilde{\mu}_{1}^{2}\boldsymbol{F}^{\top}\boldsymbol{P}^{\perp}_{\boldsymbol{\xi}}\boldsymbol{F}+\mu_{2}^{2}\boldsymbol{I}_{k}\\\
\boldsymbol{\Gamma}=\widehat{\mu}_{1}^{2}\boldsymbol{F}^{\top}\boldsymbol{F}+\mu_{3}^{2}\boldsymbol{I}_{k}.\end{cases}$
(34)

The above results show that it suffices to precisely analyze the formulation
given in (33). Moreover, note that (33) can be equivalently formulated as
follows

$\displaystyle\min_{\begin{subarray}{c}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}\leq
C_{w}\end{subarray}}\max_{\frac{\mathinner{\\!\left\lVert\boldsymbol{u}\right\rVert}}{\sqrt{n}}\leq
C_{u}}$
$\displaystyle~{}\frac{\mu_{0}\vartheta}{n\ell}\boldsymbol{u}^{\top}\boldsymbol{1}_{\ell
n}-\frac{\mathinner{\\!\left\lVert\boldsymbol{u}\right\rVert}^{2}}{2n\ell}-\frac{\boldsymbol{u}^{\top}\widehat{\boldsymbol{y}}}{n\ell}+\frac{\widetilde{\mu}_{1}\bar{\boldsymbol{\xi}}^{\top}\boldsymbol{F}\boldsymbol{w}}{\ell
n}\boldsymbol{u}^{\top}\widehat{\boldsymbol{s}}$ (35)
$\displaystyle+\frac{1}{n\ell}\boldsymbol{u}^{\top}\boldsymbol{M}^{\top}\boldsymbol{G}\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{w}+\frac{1}{n\ell}\boldsymbol{u}^{\top}\boldsymbol{T}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2},$

where $\boldsymbol{1}_{\ell n}$ denotes the all one vector of size $\ell n$,
$\widehat{\boldsymbol{y}}\in\mathbb{R}^{\ell n}$ and
$\widehat{\boldsymbol{s}}\in\mathbb{R}^{\ell n}$ are formed after performing
an $\ell$ times concatenation of the vectors $\boldsymbol{y}$ and
$\boldsymbol{s}$, respectively. Here, the matrices
$\boldsymbol{T}\in\mathbb{R}^{n\ell\times k}$ and
$\boldsymbol{M}\in\mathbb{R}^{n\times\ell n}$ are given as follows

$\displaystyle\boldsymbol{T}^{\top}=[\boldsymbol{T}_{1}^{\top},\boldsymbol{T}_{2}^{\top},\dots,\boldsymbol{T}_{\ell}^{\top}]^{\top},~{}\boldsymbol{M}=[\boldsymbol{I}_{n},\boldsymbol{I}_{n},\dots,\boldsymbol{I}_{n}].$

We can notice that the optimization problem formulated in (35) is in the form
of the multivariate PO problem given in (22). Therefore, applying the
multivariate CGMT, the corresponding multivariate AO problem can be expressed
as follows

$\displaystyle\min_{\begin{subarray}{c}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}\leq
C_{w}\end{subarray}}\max_{\frac{\mathinner{\\!\left\lVert\boldsymbol{u}\right\rVert}}{\sqrt{n}}\leq
C_{u}}$
$\displaystyle-\frac{\mathinner{\\!\left\lVert\boldsymbol{u}\right\rVert}^{2}}{2n\ell}+\frac{\widetilde{\mu}_{1}\bar{\boldsymbol{\xi}}^{\top}\boldsymbol{F}\boldsymbol{w}}{\ell
n}\boldsymbol{u}^{\top}\widehat{\boldsymbol{s}}+\frac{\mu_{0}\vartheta}{n\ell}\boldsymbol{u}^{\top}\boldsymbol{1}_{\ell
n}$ (36)
$\displaystyle-\frac{\boldsymbol{u}^{\top}\widehat{\boldsymbol{y}}}{n\ell}+\frac{1}{n\ell}\mathinner{\\!\left\lVert\boldsymbol{M}\boldsymbol{u}\right\rVert}\boldsymbol{g}_{1}^{\top}\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{w}+\frac{1}{\ell
n}\mathinner{\\!\left\lVert\boldsymbol{u}\right\rVert}\boldsymbol{g}_{2}^{\top}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}$
$\displaystyle+\frac{\mathinner{\\!\left\lVert\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}}{\ell
n}\boldsymbol{h}_{1}^{\top}\boldsymbol{M}\boldsymbol{u}+\frac{\mathinner{\\!\left\lVert\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}}{\ell
n}\boldsymbol{h}_{2}^{\top}\boldsymbol{u}+\frac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2},$

where the vectors $\boldsymbol{g}_{1}\in\mathbb{R}^{k}$,
$\boldsymbol{g}_{2}\in\mathbb{R}^{k}$, $\boldsymbol{h}_{1}\in\mathbb{R}^{n}$
and $\boldsymbol{h}_{2}\in\mathbb{R}^{\ell n}$ are independent standard
Gaussian random vectors. First, observe that the convexity assumption in
Theorem 3 is satisfied by our multivariate PO formulation in (35).
Furthermore, note that the compactness assumptions in the multivariate CGMT
framework are also satisfied by our primary problem in (35). Then, following
the multivariate CGMT framework, we focus on analyzing the multivariate AO
formulation introduced in (36). Specifically, the objective is to simplify the
multivariate AO problem and study its asymptotic properties.

####  VI-B2 Simplifying the Multivariate Auxiliary Formulation

In this part, our objective is to simplify the multivariate AO problem given
in (36). Specifically, the main objective is to express the formulation in
(36) in terms of scalar optimization variables. First, observe that the
singular value decomposition (SVD) of the matrix $\boldsymbol{M}$ can be
expressed as
$\boldsymbol{M}=\boldsymbol{U}\boldsymbol{S}\boldsymbol{V}^{\top}$, where
$\boldsymbol{U}\in\mathbb{R}^{n\times n}$ and
$\boldsymbol{V}\in\mathbb{R}^{\ell n\times\ell n}$ are two orthogonal matrices
and $\boldsymbol{S}\in\mathbb{R}^{n\times\ell n}$ is given by
$\boldsymbol{S}=[\sqrt{\ell}\boldsymbol{I}_{n}~{}~{}\boldsymbol{0}_{n\times(\ell-1)n}]$.
Therefore, the optimization problem expressed in (36) can be formulated as
follows

$\displaystyle\min_{\begin{subarray}{c}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}\leq
C_{w}\end{subarray}}\max_{\frac{\mathinner{\\!\left\lVert\boldsymbol{u}\right\rVert}}{\sqrt{n}}\leq
C_{u}}$
$\displaystyle~{}\frac{\widetilde{\mu}_{1}\bar{\boldsymbol{\xi}}^{\top}\boldsymbol{F}\boldsymbol{w}}{\ell
n}\boldsymbol{u}^{\top}\boldsymbol{V}^{\top}\widehat{\boldsymbol{s}}-\frac{\mathinner{\\!\left\lVert\boldsymbol{u}\right\rVert}^{2}}{2n\ell}+\frac{\mu_{0}\vartheta}{n\ell}\boldsymbol{u}^{\top}\boldsymbol{V}^{\top}\boldsymbol{1}_{\ell
n}$ (37)
$\displaystyle-\frac{\boldsymbol{u}^{\top}\boldsymbol{V}^{\top}\widehat{\boldsymbol{y}}}{n\ell}+\frac{\sqrt{\ell}}{n\ell}\mathinner{\\!\left\lVert\boldsymbol{u}_{1}\right\rVert}\boldsymbol{g}_{1}^{\top}\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{w}+\frac{1}{\ell
n}\mathinner{\\!\left\lVert\boldsymbol{u}\right\rVert}\boldsymbol{g}_{2}^{\top}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}$
$\displaystyle+\frac{\sqrt{\ell}}{\ell
n}\mathinner{\\!\left\lVert\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}\boldsymbol{h}_{1}^{\top}\boldsymbol{u}_{1}+\frac{\mathinner{\\!\left\lVert\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}}{\ell
n}\boldsymbol{h}_{2}^{\top}\boldsymbol{u}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2},$

where we perform the change of variable
$\boldsymbol{u}_{\text{new}}=\boldsymbol{V}^{\top}\boldsymbol{u}$, we
decompose the new vector as
$\boldsymbol{u}_{\text{new}}^{\top}=[\boldsymbol{u}_{1}^{\top},\dots,\boldsymbol{u}_{\ell}^{\top}]$
and we replace $\boldsymbol{u}_{\text{new}}$ by $\boldsymbol{u}$. Now, we
denote by $t_{1}$ and $t_{2}$ the norms of the independent vectors
$\boldsymbol{u}_{1}$ and
$\boldsymbol{u}_{-1}^{\top}=[\boldsymbol{u}_{2}^{\top},\dots,\boldsymbol{u}_{\ell}^{\top}]$,
i.e. $t_{1}=\mathinner{\\!\left\lVert\boldsymbol{u}_{1}\right\rVert}$ and
$t_{2}=\mathinner{\\!\left\lVert\boldsymbol{u}_{-1}\right\rVert}$.
Additionally, we decompose the orthogonal matrix $\boldsymbol{V}$ as follows
$\boldsymbol{V}=[\boldsymbol{V}_{1}~{}~{}\boldsymbol{V}_{2}]$, where
$\boldsymbol{V}_{1}\in\mathbb{R}^{\ell n\times n}$ and
$\boldsymbol{V}_{2}\in\mathbb{R}^{\ell n\times(\ell-1)n}$. Define the vector
$\boldsymbol{v}\in\mathbb{R}^{k}$ as
$\boldsymbol{v}=\boldsymbol{F}^{\top}{\boldsymbol{\xi}}$ and the scalar $q$ as
$q=\bar{\boldsymbol{v}}^{\top}\boldsymbol{w}$, where $\bar{\boldsymbol{v}}$ is
defined as follows
$\bar{\boldsymbol{v}}=\boldsymbol{v}/\mathinner{\\!\left\lVert\boldsymbol{v}\right\rVert}$.
Also, define the scalar $T_{p,1}$ as
$T_{p,1}=\mathinner{\\!\left\lVert\boldsymbol{v}\right\rVert}$.

Now, we are ready to further simplify the multivariate AO formulation. The
first step is to fix $t_{1}$ and $t_{2}$ and solve the formulation in (37)
over the direction of the independent vectors $\boldsymbol{u}_{1}$ and
$\boldsymbol{u}_{-1}$. Specifically, based on the result in Lemma 3, the
formulation given in (37) can be simplified as follows

$\displaystyle\min_{\begin{subarray}{c}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}\leq
C_{w}\end{subarray}}\max_{\begin{subarray}{c}0\leq{t_{1}}/{\sqrt{n}}\leq
C_{t_{1}}\\\ 0\leq{t_{2}}/{\sqrt{n}}\leq
C_{t_{2}}\end{subarray}}~{}\frac{\sqrt{\ell}t_{1}}{n\ell}\boldsymbol{g}_{1}^{\top}\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{w}+\frac{\sqrt{t_{1}^{2}+t_{2}^{2}}}{\ell
n}\boldsymbol{g}_{2}^{\top}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}-\frac{t_{1}^{2}+t_{2}^{2}}{2n\ell}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2}$
(38) $\displaystyle+\frac{t_{1}}{\ell
n}\mathinner{\\!\left\lVert\sqrt{\ell}\mathinner{\\!\left\lVert\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}\boldsymbol{h}_{1}+\mathinner{\\!\left\lVert\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}\widehat{\boldsymbol{h}}_{2}-\boldsymbol{V}_{1}^{\top}\widehat{\boldsymbol{y}}+\mu_{0}\vartheta\boldsymbol{V}_{1}^{\top}\boldsymbol{1}_{\ell
n}+\widetilde{\mu}_{1}T_{p,1}q\boldsymbol{V}_{1}^{\top}\widehat{\boldsymbol{s}}\right\rVert}$
$\displaystyle+\frac{t_{2}}{\ell
n}\mathinner{\\!\left\lVert\mathinner{\\!\left\lVert\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}\widetilde{\boldsymbol{h}}_{2}-\boldsymbol{V}_{2}^{\top}\widehat{\boldsymbol{y}}+\mu_{0}\vartheta\boldsymbol{V}_{2}^{\top}\boldsymbol{1}_{\ell
n}+\widetilde{\mu}_{1}T_{p,1}q\boldsymbol{V}_{2}^{\top}\widehat{\boldsymbol{s}}\right\rVert},$

where we decompose the Gaussian vector $\boldsymbol{h}_{2}$ as
$\boldsymbol{h}_{2}=[\widehat{\boldsymbol{h}}_{2}^{\top}~{}~{}\widetilde{\boldsymbol{h}}_{2}^{\top}]$,
where $\widehat{\boldsymbol{h}}_{2}\in\mathbb{R}^{n}$ and
$\widetilde{\boldsymbol{h}}_{2}\in\mathbb{R}^{(\ell-1)n}$. Here, $C_{t_{1}}$
and $C_{t_{2}}$ are sufficiently large positive constants that ensure the
asymptotic result in Lemma 3. Note that it remains to solve over the primal
vector $\boldsymbol{w}$ to obtain a scalar formulation of the multivariate AO
problem. We continue our analysis by defining the following optimization
problem

$\displaystyle\min_{\begin{subarray}{c}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}\leq
C_{w}\end{subarray}}\max_{\begin{subarray}{c}0\leq{t_{1}}/{\sqrt{n}}\leq
C_{t_{1}}\\\ 0\leq{t_{2}}/{\sqrt{n}}\leq
C_{t_{2}}\end{subarray}}~{}\frac{\sqrt{\ell}t_{1}}{n\ell}\boldsymbol{g}_{1}^{\top}\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{w}+\frac{\sqrt{t_{1}^{2}+t_{2}^{2}}}{\ell
n}\boldsymbol{g}_{2}^{\top}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}-\frac{t_{1}^{2}+t_{2}^{2}}{2n\ell}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2}$
(39) $\displaystyle+\frac{t_{1}}{\ell
n}\sqrt{{\ell}\mathinner{\\!\left\lVert\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}^{2}\mathinner{\\!\left\lVert\boldsymbol{h}_{1}\right\rVert}^{2}+\mathinner{\\!\left\lVert\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}^{2}\mathinner{\\!\left\lVert\widehat{\boldsymbol{h}}_{2}\right\rVert}^{2}+\mathinner{\\!\left\lVert\mu_{0}\vartheta\boldsymbol{V}_{1}^{\top}\boldsymbol{1}_{\ell
n}+\widetilde{\mu}_{1}T_{p,1}q\boldsymbol{V}_{1}^{\top}\widehat{\boldsymbol{s}}-\boldsymbol{V}_{1}^{\top}\widehat{\boldsymbol{y}}\right\rVert}^{2}}$
$\displaystyle+\frac{t_{2}}{\ell
n}\sqrt{\mathinner{\\!\left\lVert\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}^{2}\mathinner{\\!\left\lVert\widetilde{\boldsymbol{h}}_{2}\right\rVert}^{2}+\mathinner{\\!\left\lVert\mu_{0}\vartheta\boldsymbol{V}_{2}^{\top}\boldsymbol{1}_{\ell
n}-\boldsymbol{V}_{2}^{\top}\widehat{\boldsymbol{y}}+\widetilde{\mu}_{1}T_{p,1}q\boldsymbol{V}_{2}^{\top}\widehat{\boldsymbol{s}}\right\rVert}^{2}}.$

Note that the difference between the cost functions of the formulations in
(38) and (39) are terms that converge in probability to zero. Before showing
the asymptotic equivalence between the formulations in (38) and (39), we
provide important convexity properties of the optimization problem in (39) as
given in the following lemma.

######  Lemma 4 (Strong–convexity of (39)).

Define $\widehat{f}_{p,2}$ as the cost function of the problem in (39). Then,
$\widehat{f}_{p,2}$ is strongly convex in the vector $\boldsymbol{w}$ where
$\lambda$ is a strong convexity parameter. Moreover, it is strongly concave in
the variables $t_{1}$ and $t_{2}$ in the feasibility sets where $-1$ is a
strong concavity parameter.

###### Proof.

The strong convexity can be proved by observing that the cost function of (39)
is a positive sum of convex and strongly convex functions in terms of
$\boldsymbol{w}$ for fixed feasible $t_{1}$ and $t_{2}$. Moreover, note that
the term
${\sqrt{t_{1}^{2}+t_{2}^{2}}}\boldsymbol{g}_{2}^{\top}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}$
can be replaced with
$t_{1}\boldsymbol{g}_{21}^{\top}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}+t_{2}\boldsymbol{g}_{22}^{\top}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}$
without changing the statistics of our formulation, where
$\boldsymbol{g}_{21}$ and $\boldsymbol{g}_{22}$ are two independent Gaussian
vectors. Then, one can see that the cost function of (39) is strongly concave
in the variables $t_{1}$ and $t_{2}$ where $-1$ is a strong concavity
parameter. ∎

Lemma 4 provides important convexity properties of the optimization problem
formulated in (39). These properties are essential to prove the equivalence
between (38) and (39) as stated in the following lemma.

######  Lemma 5 (High-dimensional Equivalence I).

Define $\widehat{\mathcal{S}}^{\star}_{p,1}$ and
$\widehat{\mathcal{S}}^{\star}_{p,2}$ as the sets of optimal solutions of the
minimization problems in (38) and (39), respectively. Moreover, let
$\widehat{O}^{\star}_{p,1}$ and $\widehat{O}^{\star}_{p,2}$ be the optimal
objective values of the optimization problems in (38) and (39), respectively.
Then, the following convergence in probability holds

$\mathinner{\\!\left\lvert\widehat{O}^{\star}_{p,1}-\widehat{O}^{\star}_{p,2}\right\rvert}\overset{p\to+\infty}{\longrightarrow}0,~{}\text{and}~{}\mathbb{D}(\widehat{\mathcal{S}}^{\star}_{p,1},\widehat{\mathcal{S}}^{\star}_{p,2})\overset{p\to+\infty}{\longrightarrow}0,$
(40)

where $\mathbb{D}(\mathcal{A},\mathcal{B})$ denotes the deviation between the
sets $\mathcal{A}$ and $\mathcal{B}$ and is defined as
$\mathbb{D}(\mathcal{A},\mathcal{B})=\sup_{\boldsymbol{x}_{1}\in\mathcal{A}}\inf_{\boldsymbol{x}_{2}\in\mathcal{B}}\mathinner{\\!\left\lVert\boldsymbol{x}_{1}-\boldsymbol{x}_{2}\right\rVert}_{2}$.

The detailed proof of Lemma 5 is deferred to Appendix VIII-A. Lemma 5
particularly shows that the optimization problems in (38) and (39) are
asymptotically equivalent. Then, it suffices to precisely analyze the
formulation in (39). To solve over the primal vector $\boldsymbol{w}$, we
introduce two independent scalar optimization variables $\tau_{1}$ and
$\tau_{2}$ where they both solve optimization problems of the following form

$\displaystyle\sqrt{x}=\inf_{\tau>0}~{}\frac{\tau}{2}+\frac{x}{2\tau},~{}\text{for
any}~{}x\geq 0.$ (41)

Here, note that the optimal solution of the problem in (41) can be expressed
as $\tau^{\star}=\sqrt{x}$. Next, we use the identity in (41) to transform the
non-smooth square roots in the cost function of the formulation given in (39)
to smooth terms. This is an essential step to solve over the primal vector
$\boldsymbol{w}$. Specifically, based on the result in Lemma 2, our
multivariate AO formulation given in (39) can be expressed as follows

$\displaystyle\min_{\begin{subarray}{c}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}\leq
C_{w}\end{subarray}}\max_{\begin{subarray}{c}0\leq{t_{1}}\leq C_{t_{1}}\\\
0\leq{t_{2}}\leq
C_{t_{2}}\end{subarray}}\inf_{\begin{subarray}{c}\tau_{1}>0\\\
\tau_{2}>0\end{subarray}}~{}\frac{\sqrt{\ell}t_{1}}{\sqrt{n}\ell}\boldsymbol{g}_{1}^{\top}\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{w}+\frac{\sqrt{t_{1}^{2}+t_{2}^{2}}}{\ell\sqrt{n}}\boldsymbol{g}_{2}^{\top}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}+\frac{\tau_{1}t_{1}}{2\ell}+\frac{\tau_{2}t_{2}}{2\ell}-\frac{t_{1}^{2}+t_{2}^{2}}{2\ell}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2}$
(42)
$\displaystyle+\frac{t_{1}\mathinner{\\!\left\lVert\boldsymbol{h}_{1}\right\rVert}^{2}}{2\tau_{1}n}\mathinner{\\!\left\lVert\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}^{2}+\frac{\zeta_{p,t,\tau}}{2n}\mathinner{\\!\left\lVert\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}^{2}+\frac{t_{1}}{2\tau_{1}\ell
n}\mathinner{\\!\left\lVert\mu_{0}\vartheta\boldsymbol{V}_{1}^{\top}\boldsymbol{1}_{\ell
n}-\boldsymbol{V}_{1}^{\top}\widehat{\boldsymbol{y}}+\widetilde{\mu}_{1}T_{p,1}q\boldsymbol{V}_{1}^{\top}\widehat{\boldsymbol{s}}\right\rVert}^{2}$
$\displaystyle+\frac{t_{2}}{2\tau_{2}\ell
n}\mathinner{\\!\left\lVert\mu_{0}\vartheta\boldsymbol{V}_{2}^{\top}\boldsymbol{1}_{\ell
n}-\boldsymbol{V}_{2}^{\top}\widehat{\boldsymbol{y}}+\widetilde{\mu}_{1}T_{p,1}q\boldsymbol{V}_{2}^{\top}\widehat{\boldsymbol{s}}\right\rVert}^{2},$

where
$\zeta_{p,t,\tau}={t_{1}\mathinner{\\!\left\lVert\widehat{\boldsymbol{h}}_{2}\right\rVert}^{2}}/{(\ell\tau_{1})}+{t_{2}\mathinner{\\!\left\lVert\widetilde{\boldsymbol{h}}_{2}\right\rVert}^{2}}/{(\ell\tau_{2})}$.
Here, we also perform the change of variable $t_{1,\text{new}}=t_{1}/\sqrt{n}$
and $t_{2,\text{new}}=t_{2}/\sqrt{n}$, then, replace $t_{1,\text{new}}$ and
$t_{2,\text{new}}$ by $t_{1}$ and $t_{2}$. Note that the feasibility sets of
the optimization variables $\tau_{1}$ and $\tau_{2}$ are open unbounded sets.
To simplify the analysis, we show that the feasibility sets of the variables
$\tau_{1}$ and $\tau_{2}$ can be restricted to compact sets with probability
going to $1$ as $p$ grows to $+\infty$ as stated in the following lemma.

######  Lemma 6 (Additional Compactness).

There exist positive constants independent of $p$, $c_{\tau_{1}}>0$,
$C_{\tau_{1}}>0$, $c_{\tau_{2}}>0$ and $C_{\tau_{2}}>0$, such that the
following convergence in probability holds

$\displaystyle\mathbb{P}(c_{\tau_{1}}\leq\widehat{\tau}_{1}\leq
C_{\tau_{1}})\xrightarrow{p\to\infty}1,~{}\mathbb{P}(c_{\tau_{2}}\leq\widehat{\tau}_{2}\leq
C_{\tau_{2}})\xrightarrow{p\to\infty}1,$ (43)

where $\widehat{\tau}_{1}$ and $\widehat{\tau}_{2}$ are the optimal solutions
of the formulation in (LABEL:ana_fm7).

The detailed proof of Lemma 6 is provided in Appendix VIII-B. Based on Lemmas
4 and 6, the optimization problem given in (LABEL:ana_fm7) is asymptotically
equivalent to the following problem

$\displaystyle\min_{\begin{subarray}{c}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}\leq
C_{w}\end{subarray}}\max_{\begin{subarray}{c}0\leq{t_{1}}\leq C_{t_{1}}\\\
0\leq{t_{2}}\leq
C_{t_{2}}\end{subarray}}\min_{\begin{subarray}{c}c_{\tau_{1}}\leq{\tau_{1}}\leq
C_{\tau_{1}}\\\ c_{\tau_{2}}\leq{\tau_{2}}\leq
C_{\tau_{2}}\end{subarray}}~{}\frac{\sqrt{\ell}t_{1}}{\sqrt{n}\ell}\boldsymbol{g}_{1}^{\top}\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{w}+\frac{\sqrt{t_{1}^{2}+t_{2}^{2}}}{\ell\sqrt{n}}\boldsymbol{g}_{2}^{\top}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}+\frac{\tau_{1}t_{1}}{2\ell}+\frac{\tau_{2}t_{2}}{2\ell}-\frac{t_{1}^{2}+t_{2}^{2}}{2\ell}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{w}\right\rVert}^{2}$
(44)
$\displaystyle+\frac{t_{1}\mathinner{\\!\left\lVert\boldsymbol{h}_{1}\right\rVert}^{2}}{2\tau_{1}n}\mathinner{\\!\left\lVert\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}^{2}+\frac{\zeta_{p,t,\tau}}{2n}\mathinner{\\!\left\lVert\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{w}\right\rVert}^{2}+\frac{t_{1}}{2\tau_{1}\ell
n}\mathinner{\\!\left\lVert\mu_{0}\vartheta\boldsymbol{V}_{1}^{\top}\boldsymbol{1}_{\ell
n}-\boldsymbol{V}_{1}^{\top}\widehat{\boldsymbol{y}}+\widetilde{\mu}_{1}T_{p,1}q\boldsymbol{V}_{1}^{\top}\widehat{\boldsymbol{s}}\right\rVert}^{2}$
$\displaystyle+\frac{t_{2}}{2\tau_{2}\ell
n}\mathinner{\\!\left\lVert\mu_{0}\vartheta\boldsymbol{V}_{2}^{\top}\boldsymbol{1}_{\ell
n}-\boldsymbol{V}_{2}^{\top}\widehat{\boldsymbol{y}}+\widetilde{\mu}_{1}T_{p,1}q\boldsymbol{V}_{2}^{\top}\widehat{\boldsymbol{s}}\right\rVert}^{2}.$

Now, we are ready to simplify the formulation in (LABEL:ana_fm7p) over the
optimization vector $\boldsymbol{w}$. We start our analysis by decomposing the
optimization vector $\boldsymbol{w}\in\mathbb{R}^{k}$ as follows

$\displaystyle\boldsymbol{w}=q\bar{\boldsymbol{v}}+\boldsymbol{B}^{\perp}_{\boldsymbol{v}}\boldsymbol{r},$
(45)

where $\boldsymbol{r}\in\mathbb{R}^{k-1}$ and
$\boldsymbol{B}^{\perp}_{\boldsymbol{v}}\in\mathbb{R}^{k\times(k-1)}$ is
formed by an orthonormal basis orthogonal to the vector
$\boldsymbol{v}\in\mathbb{R}^{k}$. Based on the result in Lemma 2, one can
equivalently formulate the problem in (LABEL:ana_fm7p) as follows

$\displaystyle\min_{\begin{subarray}{c}\mathinner{\\!\left\lvert
q\right\rvert}\leq C_{q}\\\
\mathinner{\\!\left\lVert\boldsymbol{r}\right\rVert}\leq
C_{r}\end{subarray}}\max_{\begin{subarray}{c}0\leq{t_{1}}\leq C_{t_{1}}\\\
0\leq{t_{2}}\leq
C_{t_{2}}\end{subarray}}\min_{\begin{subarray}{c}c_{\tau_{1}}\leq{\tau_{1}}\leq
C_{\tau_{1}}\\\ c_{\tau_{2}}\leq{\tau_{2}}\leq
C_{\tau_{2}}\end{subarray}}~{}\frac{\sqrt{\ell}t_{1}}{\sqrt{n}\ell}\boldsymbol{g}_{1}^{\top}\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{B}^{\perp}_{\boldsymbol{v}}\boldsymbol{r}+\frac{\sqrt{t_{1}^{2}+t_{2}^{2}}}{\ell\sqrt{n}}\boldsymbol{g}_{2}^{\top}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{B}^{\perp}_{\boldsymbol{v}}\boldsymbol{r}+\frac{\tau_{1}t_{1}}{2\ell}+\frac{\tau_{2}t_{2}}{2\ell}-\frac{t_{1}^{2}+t_{2}^{2}}{2\ell}$
(46)
$\displaystyle+\frac{t_{1}\mathinner{\\!\left\lVert\boldsymbol{h}_{1}\right\rVert}^{2}}{2\tau_{1}n}\mathinner{\\!\left\lVert\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{B}^{\perp}_{\boldsymbol{v}}\boldsymbol{r}+q\boldsymbol{\Sigma}^{\frac{1}{2}}\bar{\boldsymbol{v}}\right\rVert}^{2}+\frac{\zeta_{p,t,\tau}}{2n}\mathinner{\\!\left\lVert\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{B}^{\perp}_{\boldsymbol{v}}\boldsymbol{r}+q\boldsymbol{\Gamma}^{\frac{1}{2}}\bar{\boldsymbol{v}}\right\rVert}^{2}+\tfrac{\lambda}{2}(q^{2}+\mathinner{\\!\left\lVert\boldsymbol{r}\right\rVert}^{2})$
$\displaystyle+\frac{t_{1}}{2\tau_{1}\ell
n}\mathinner{\\!\left\lVert\mu_{0}\vartheta\boldsymbol{V}_{1}^{\top}\boldsymbol{1}_{\ell
n}-\boldsymbol{V}_{1}^{\top}\widehat{\boldsymbol{y}}+\widetilde{\mu}_{1}T_{p,1}q\boldsymbol{V}_{1}^{\top}\widehat{\boldsymbol{s}}\right\rVert}^{2}+\frac{t_{2}}{2\tau_{2}\ell
n}\mathinner{\\!\left\lVert\mu_{0}\vartheta\boldsymbol{V}_{2}^{\top}\boldsymbol{1}_{\ell
n}-\boldsymbol{V}_{2}^{\top}\widehat{\boldsymbol{y}}+\widetilde{\mu}_{1}T_{p,1}q\boldsymbol{V}_{2}^{\top}\widehat{\boldsymbol{s}}\right\rVert}^{2},$

where $C_{q}>0$ and $C_{r}>0$ are two positive constants selected to satisfy
the asymptotic result in Lemma 2. Here, we also drop terms that converge in
probability to zero. One way to justify this step is using similar analysis as
in Lemma 5. Note that the convexity results in Lemma 4 are still satisfied by
the formulation in (LABEL:ana_fm8). Specifically, the cost function in
(LABEL:ana_fm8) is jointly strongly convex in the minimization variables and
jointly strongly concave in the maximization variables.

Now, it remains to solve over the optimization vector
$\boldsymbol{r}\in\mathbb{R}^{k-1}$. To solve over $\boldsymbol{r}$, we
interchange the minimization over $\boldsymbol{r}$ and the maximization over
$t_{1}$ and $t_{2}$. This step is justified using the result in [36]. The cost
function of the optimization vector $\boldsymbol{r}$ can then be expressed as
follows

$\displaystyle g(\boldsymbol{r})$
$\displaystyle=\frac{\sqrt{\ell}t_{1}}{\sqrt{n}\ell}\boldsymbol{g}_{1}^{\top}\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{B}^{\perp}_{\boldsymbol{v}}\boldsymbol{r}+\frac{\sqrt{t_{1}^{2}+t_{2}^{2}}}{\ell\sqrt{n}}\boldsymbol{g}_{2}^{\top}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{B}^{\perp}_{\boldsymbol{v}}\boldsymbol{r}+\frac{t_{1}\mathinner{\\!\left\lVert\boldsymbol{h}_{1}\right\rVert}^{2}}{2\tau_{1}n}\mathinner{\\!\left\lVert\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{B}^{\perp}_{\boldsymbol{v}}\boldsymbol{r}+q\boldsymbol{\Sigma}^{\frac{1}{2}}\bar{\boldsymbol{v}}\right\rVert}^{2}+\tfrac{\lambda}{2}\mathinner{\\!\left\lVert\boldsymbol{r}\right\rVert}^{2}$
(47)
$\displaystyle+\frac{\zeta_{p,t,\tau}}{2n}\mathinner{\\!\left\lVert\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{B}^{\perp}_{\boldsymbol{v}}\boldsymbol{r}+q\boldsymbol{\Gamma}^{\frac{1}{2}}\bar{\boldsymbol{v}}\right\rVert}^{2},$

where we ignore the terms independent of $\boldsymbol{r}$. Note that the
function $g(\cdot)$ is convex and smooth. Before solving the minimization
problem of the function $g(\cdot)$, we define the matrix
$\boldsymbol{G}\in\mathbb{R}^{k\times k}$ and the vectors
$\boldsymbol{f}\in\mathbb{R}^{k-1}$ and $\boldsymbol{z}\in\mathbb{R}^{k-1}$ as
follows

$\displaystyle\begin{cases}\boldsymbol{G}=\frac{t_{1}\mathinner{\\!\left\lVert\boldsymbol{h}_{1}\right\rVert}^{2}}{\tau_{1}n}\boldsymbol{\Sigma}+\frac{\zeta_{p,t,\tau}}{n}\boldsymbol{\Gamma}\\\
\boldsymbol{f}=\frac{t_{1}\mathinner{\\!\left\lVert\boldsymbol{h}_{1}\right\rVert}^{2}}{\tau_{1}n}\bar{\boldsymbol{B}}^{\perp}_{\boldsymbol{v}}\boldsymbol{\Sigma}\bar{\boldsymbol{v}}+\frac{\zeta_{p,t,\tau}}{n}\bar{\boldsymbol{B}}^{\perp}_{\boldsymbol{v}}\boldsymbol{\Gamma}\bar{\boldsymbol{v}}\\\
\boldsymbol{z}=\frac{\sqrt{\ell}t_{1}}{\sqrt{n}\ell}\bar{\boldsymbol{B}}^{\perp}_{\boldsymbol{v}}\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{g}_{1}+\frac{\sqrt{t_{1}^{2}+t_{2}^{2}}}{\ell\sqrt{n}}\bar{\boldsymbol{B}}^{\perp}_{\boldsymbol{v}}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{g}_{2},\end{cases}$
(48)

where
$\bar{\boldsymbol{B}}^{\perp}_{\boldsymbol{v}}=({\boldsymbol{B}}^{\perp}_{\boldsymbol{v}})^{\top}$.
After computing the derivative of the function $g(\cdot)$ and setting it to
zero, the optimal solution of the unconstrained version of minimizing the
function $g(\cdot)$ can be expressed as follows

$\displaystyle\widetilde{\boldsymbol{r}}^{\star}=-\left[\bar{\boldsymbol{B}}^{\perp}_{\boldsymbol{v}}\boldsymbol{G}\boldsymbol{B}^{\perp}_{\boldsymbol{v}}+\lambda\boldsymbol{I}_{k-1}\right]^{-1}\left[q\boldsymbol{f}+\boldsymbol{z}\right].$
(49)

Similar to the analysis in Lemmas 2, 3 and 6, one can show that the norm of
the optimal vector $\widetilde{\boldsymbol{r}}^{\star}$ is bounded. This means
that $\widetilde{\boldsymbol{r}}^{\star}$ is an optimal solution of the
formulation in (LABEL:ana_fm8). Then, the optimal loss function can be
expressed as follows

$\displaystyle g^{\star}$
$\displaystyle=-\frac{q^{2}}{2}\boldsymbol{f}^{\top}\left[\bar{\boldsymbol{B}}^{\perp}_{\boldsymbol{v}}\boldsymbol{G}\boldsymbol{B}^{\perp}_{\boldsymbol{v}}+\lambda\boldsymbol{I}_{k-1}\right]^{-1}\boldsymbol{f}+\frac{q^{2}}{2}\bar{\boldsymbol{v}}^{\top}\boldsymbol{G}\bar{\boldsymbol{v}}$
(50)
$\displaystyle-\frac{1}{2}\boldsymbol{z}^{\top}\left[\bar{\boldsymbol{B}}^{\perp}_{\boldsymbol{v}}\boldsymbol{G}\boldsymbol{B}^{\perp}_{\boldsymbol{v}}+\lambda\boldsymbol{I}_{k-1}\right]^{-1}\boldsymbol{z}.$

Based on the SVD decomposition of the matrix $\boldsymbol{M}$, it can be
checked that the last term in the multivariate AO formulation given in
(LABEL:ana_fm8) is zero. Then, the formulation given in (LABEL:ana_fm8) can be
expressed as follows

$\displaystyle\min_{\begin{subarray}{c}\mathinner{\\!\left\lvert
q\right\rvert}\leq
C_{q}\end{subarray}}\max_{\begin{subarray}{c}0\leq{t_{1}}\leq C_{t_{1}}\\\
0\leq{t_{2}}\leq
C_{t_{2}}\end{subarray}}\min_{\begin{subarray}{c}c_{\tau_{1}}\leq{\tau_{1}}\leq
C_{\tau_{1}}\\\ c_{\tau_{2}}\leq{\tau_{2}}\leq C_{\tau_{2}}\end{subarray}}$
$\displaystyle\frac{\tau_{1}t_{1}}{2\ell}-\frac{t_{1}^{2}+t_{2}^{2}}{2\ell}+\frac{\tau_{2}t_{2}}{2\ell}+\frac{\lambda+V_{p,2}(\boldsymbol{t},\boldsymbol{\tau})-V_{p,3}(\boldsymbol{t},\boldsymbol{\tau})}{2}q^{2}$
(51)
$\displaystyle+\frac{t_{1}}{2\tau_{1}n}\mathinner{\\!\left\lVert\mu_{0}\vartheta\boldsymbol{1}_{n}-{\boldsymbol{y}}+\widetilde{\mu}_{1}T_{p,1}q{\boldsymbol{s}}\right\rVert}^{2}-\frac{1}{2}V_{p,4}(\boldsymbol{t},\boldsymbol{\tau}),$

where $\boldsymbol{t}=[t_{1},t_{2}]^{\top}$ and
$\boldsymbol{\tau}=[\tau_{1},\tau_{2}]^{\top}$. Here, $\boldsymbol{1}_{n}$
denotes the vector of all one with size $n$ and the functions
$V_{p,2}(\cdot,\cdot)$, $V_{p,3}(\cdot,\cdot)$ and $V_{p,4}(\cdot,\cdot)$
depend on the optimization variables and are given by

$\displaystyle\begin{cases}V_{p,2}(\boldsymbol{t},\boldsymbol{\tau})=\bar{\boldsymbol{v}}^{\top}\boldsymbol{G}\bar{\boldsymbol{v}}\\\
V_{p,3}(\boldsymbol{t},\boldsymbol{\tau})=\boldsymbol{f}^{\top}\left[\bar{\boldsymbol{B}}^{\perp}_{\boldsymbol{v}}\boldsymbol{G}\boldsymbol{B}^{\perp}_{\boldsymbol{v}}+\lambda\boldsymbol{I}_{k-1}\right]^{-1}\boldsymbol{f}\\\
V_{p,4}(\boldsymbol{t},\boldsymbol{\tau})=\boldsymbol{z}^{\top}\left[\bar{\boldsymbol{B}}^{\perp}_{\boldsymbol{v}}\boldsymbol{G}\boldsymbol{B}^{\perp}_{\boldsymbol{v}}+\lambda\boldsymbol{I}_{k-1}\right]^{-1}\boldsymbol{z}.\end{cases}$
(52)

Note that we simplified the multivariate AO formulation given in (36) to a
scalar optimization problem as given in (51). Then, it remains to study the
asymptotic properties of the scalar formulation in (51). We refer to this
problem as the scalar formulation.

####  VI-B3 Asymptotic Analysis of the Scalar Formulation

In this part, we study the asymptotic properties of the scalar formulation in
(51) corresponding to the multivariate AO problem. Based on Assumption 5 and
the result in [37, Proposition 3], the random variable ${T}_{p,1}$ converges
pointwisely in probability to the scalar $T_{1}$ defined as follows

$\displaystyle
T_{p,1}=\sqrt{{\boldsymbol{\xi}}^{\top}\boldsymbol{F}\boldsymbol{F}^{\top}{\boldsymbol{\xi}}}\xrightarrow{~{}p\to\infty~{}}T_{1}=\sqrt{\delta\mathbb{E}_{\kappa}[\kappa]},$
(53)

where the expectations are over the probability distribution
$\mathbb{P}_{\kappa}(\cdot)$ defined in Assumption 5. Furthermore, the random
function $V_{p,2}(\cdot,\cdot)$ can be expressed as follows

$\displaystyle V_{p,2}(\boldsymbol{t},\boldsymbol{\tau})$
$\displaystyle=\frac{1}{T_{p,1}^{2}}{\boldsymbol{\xi}}^{\top}\boldsymbol{F}\Big{(}\frac{t_{1}\mathinner{\\!\left\lVert\boldsymbol{h}_{1}\right\rVert}^{2}}{\tau_{1}n}\boldsymbol{\Sigma}+\frac{\zeta_{p,t,\tau}}{n}\boldsymbol{\Gamma}\Big{)}\boldsymbol{F}^{\top}{\boldsymbol{\xi}}$
(54)
$\displaystyle=\frac{1}{T_{p,1}^{2}}{\boldsymbol{\xi}}^{\top}\boldsymbol{F}\Big{[}\Big{(}\frac{t_{1}\mathinner{\\!\left\lVert\boldsymbol{h}_{1}\right\rVert}^{2}\widetilde{\mu}_{1}^{2}}{\tau_{1}n}+\frac{\zeta_{p,t,\tau}\widehat{\mu}_{1}^{2}}{n}\Big{)}\boldsymbol{F}^{\top}\boldsymbol{F}-\frac{t_{1}\mathinner{\\!\left\lVert\boldsymbol{h}_{1}\right\rVert}^{2}\widetilde{\mu}_{1}^{2}}{\tau_{1}n}$
$\displaystyle\times\boldsymbol{F}^{\top}{\boldsymbol{\xi}}{\boldsymbol{\xi}}^{\top}\boldsymbol{F}+\Big{(}\frac{t_{1}\mathinner{\\!\left\lVert\boldsymbol{h}_{1}\right\rVert}^{2}{\mu}_{2}^{2}}{\tau_{1}n}+\frac{\zeta_{p,t,\tau}{\mu}_{3}^{2}}{n}\Big{)}\boldsymbol{I}_{k}\Big{]}\boldsymbol{F}^{\top}{\boldsymbol{\xi}},$

where $\boldsymbol{t}=[t_{1},t_{2}]^{\top}$ and
$\boldsymbol{\tau}=[\tau_{1},\tau_{2}]^{\top}$. Then, using the theoretical
results in [37] and based on Assumption 5, the random function
$V_{p,2}(\cdot,\cdot)$ converges pointwisely in probability as follows

$\displaystyle
V_{p,2}(\boldsymbol{t},\boldsymbol{\tau})\xrightarrow{~{}p\to+\infty~{}}V_{2}(\boldsymbol{t},\boldsymbol{\tau})$
$\displaystyle=\frac{\delta}{T_{1}^{2}}\Big{(}\frac{t_{1}\widetilde{\mu}_{1}^{2}}{\tau_{1}}+\zeta_{t,\tau}\widehat{\mu}_{1}^{2}\Big{)}\mathbb{E}_{\kappa}[\kappa^{2}]-\delta^{2}\frac{t_{1}\widetilde{\mu}_{1}^{2}}{\tau_{1}}\mathbb{E}_{\kappa}[\kappa]^{2}$
(55)
$\displaystyle+\Big{(}\frac{t_{1}{\mu}_{2}^{2}}{\tau_{1}}+\zeta_{t,\tau}{\mu}_{3}^{2}\Big{)}\delta\mathbb{E}_{\kappa}[\kappa],$

where
$\zeta_{t,\tau}={t_{1}}/{(\ell\tau_{1})}+{t_{2}(\ell-1)}/{(\ell\tau_{2})}$ and
the expectations are over the probability distribution
$\mathbb{P}_{\kappa}(\cdot)$ defined in Assumption 5. The theoretical results
in [37, Proposition 3] also show that the random function
$V_{p,4}(\cdot,\cdot)$ satisfies the following asymptotic result

$\displaystyle
V_{p,4}(\boldsymbol{t},\boldsymbol{\tau})-\widehat{V}_{p,4}(\boldsymbol{t},\boldsymbol{\tau})\xrightarrow{~{}p\to+\infty~{}}0,$
(56)

where the random function $\widehat{V}_{p,4}(\cdot,\cdot)$ is defined as
follows

$\displaystyle\widehat{V}_{p,4}(\boldsymbol{t},\boldsymbol{\tau})$
$\displaystyle=\frac{t_{1}^{2}}{{n}\ell}\text{Tr}\Big{[}\widehat{\boldsymbol{\Sigma}}^{\frac{1}{2}}\left[\widehat{\boldsymbol{G}}+\lambda\boldsymbol{I}_{k-1}\right]^{-1}\widehat{\boldsymbol{\Sigma}}^{\frac{1}{2}}\Big{]}$
(57)
$\displaystyle+\frac{{t_{1}^{2}+t_{2}^{2}}}{\ell^{2}{n}}\text{Tr}\Big{[}\boldsymbol{\Gamma}^{\frac{1}{2}}\left[\widehat{\boldsymbol{G}}+\lambda\boldsymbol{I}_{k-1}\right]^{-1}\boldsymbol{\Gamma}^{\frac{1}{2}}\Big{]}.$

Here, $\text{Tr}[.]$ represents the trace. Additionally, the matrix
$\widehat{\boldsymbol{\Sigma}}$ is given as
$\widehat{\boldsymbol{\Sigma}}=\widetilde{\mu}_{1}^{2}\boldsymbol{F}^{\top}\boldsymbol{F}+\mu_{2}^{2}\boldsymbol{I}_{k}$
and the matrix $\widehat{\boldsymbol{G}}$ has the following expression
$\widehat{\boldsymbol{G}}={t_{1}\mathinner{\\!\left\lVert\boldsymbol{h}_{1}\right\rVert}^{2}}/{(\tau_{1}n)}\widehat{\boldsymbol{\Sigma}}+{\zeta_{p,t,\tau}}/{n}\boldsymbol{\Gamma}$.
Using again [37, Proposition 3] and Assumption 5, we can also see that the
random function $V_{p,4}(\cdot,\cdot)$ converges in probability to the
function $V_{4}(\cdot,\cdot)$ defined as follows

$\displaystyle V_{4}(\boldsymbol{t},\boldsymbol{\tau})$
$\displaystyle=\frac{\eta
t_{1}^{2}}{\ell}\mathbb{E}_{\kappa}\Big{[}\frac{\widetilde{\mu}_{1}^{2}\kappa+\mu_{2}^{2}}{\frac{t_{1}}{\tau_{1}}(\widetilde{\mu}_{1}^{2}\kappa+\mu_{2}^{2})+\zeta_{t,\tau}(\widehat{\mu}_{1}^{2}\kappa+\mu_{3}^{2})+\lambda}\Big{]}$
$\displaystyle+\frac{\eta{(t_{1}^{2}+t_{2}^{2})}}{\ell^{2}}\mathbb{E}_{\kappa}\Big{[}\frac{\widehat{\mu}_{1}^{2}\kappa+\mu_{3}^{2}}{\frac{t_{1}}{\tau_{1}}(\widetilde{\mu}_{1}^{2}\kappa+\mu_{2}^{2})+\zeta_{t,\tau}(\widehat{\mu}_{1}^{2}\kappa+\mu_{3}^{2})+\lambda}\Big{]}.$

Now, it remains to study the asymptotic properties of the random function
$V_{p,3}(\cdot,\cdot)$. Based on the block matrix inversion lemma, it can be
checked that the random function $V_{p,3}(\cdot,\cdot)$ satisfies the
following

$\displaystyle
V_{p,3}(\boldsymbol{t},\boldsymbol{\tau})=V_{p,2}(\boldsymbol{t},\boldsymbol{\tau})+\lambda-\frac{1}{T_{p,2}(\boldsymbol{t},\boldsymbol{\tau})}.$
(58)

Here, the random function $T_{p,2}(\cdot,\cdot)$ is defined as follows

$\displaystyle
T_{p,2}(\boldsymbol{t},\boldsymbol{\tau})=\bar{\boldsymbol{v}}^{\top}\left[\boldsymbol{G}+\lambda\boldsymbol{I}_{k}\right]^{-1}\bar{\boldsymbol{v}}.$
(59)

Using the matrix inversion lemma, it can be checked that the random function
$T_{p,2}(\cdot,\cdot)$ converges in probability to the function
$T_{2,\lambda}(\cdot,\cdot)$ defined as follows

$\displaystyle
T_{2,\lambda}(\boldsymbol{t},\boldsymbol{\tau})=\frac{\delta}{T_{1}^{2}}\mathbb{E}\Big{[}\frac{\kappa}{g_{\kappa,\lambda}(t,\tau)}\Big{]}/\Big{(}1-\frac{\widetilde{\mu}_{1}^{2}t_{1}\delta}{\tau_{1}}\mathbb{E}\Big{[}\frac{\kappa}{g_{\kappa,\lambda}(t,\tau)}\Big{]}\Big{)},$

where the function $g_{\kappa,\lambda}(\cdot,\cdot)$ is defined in Section IV.
Additionally, using the weak law of large numbers (WLLN), we have the
following convergence property

$\displaystyle\frac{1}{n}\mathinner{\\!\left\lVert\mu_{0}\vartheta\boldsymbol{1}_{n}-{\boldsymbol{y}}+\widetilde{\mu}_{1}T_{p,1}q{\boldsymbol{s}}\right\rVert}^{2}\xrightarrow{~{}p\to+\infty~{}}\mu_{0}^{2}\vartheta^{2}+\gamma_{1}+\widetilde{\mu}_{1}^{2}T_{1}^{2}q^{2}-2\widetilde{\mu}_{1}T_{1}q\gamma_{2}-2\mu_{0}\vartheta\gamma_{3}.$
(60)

Here, $\gamma_{1}$, $\gamma_{2}$ and $\gamma_{3}$ depend on the data
distribution and are defined as $\gamma_{1}=\mathbb{E}[y^{2}]$,
$\gamma_{2}=\mathbb{E}[ys]$, $\gamma_{3}=\mathbb{E}[y]$, where $y=\varphi(s)$,
and $s$ is a standard Gaussian random variable. This proves that the cost
function of the following deterministic problem

$\displaystyle\min_{\begin{subarray}{c}\mathinner{\\!\left\lvert
q\right\rvert}\leq
C_{q}\end{subarray}}\max_{\begin{subarray}{c}0\leq{t_{1}}\leq C_{t_{1}}\\\
0\leq{t_{2}}\leq
C_{t_{2}}\end{subarray}}\min_{\begin{subarray}{c}c_{\tau_{1}}\leq{\tau_{1}}\leq
C_{\tau_{1}}\\\ c_{\tau_{2}}\leq{\tau_{2}}\leq C_{\tau_{2}}\end{subarray}}$
$\displaystyle\frac{\tau_{1}t_{1}+\tau_{2}t_{2}}{2\ell}-\frac{t_{1}^{2}+t_{2}^{2}}{2\ell}+\frac{q^{2}}{2T_{2,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}-\frac{\eta
T_{3,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}{2}$ (61)
$\displaystyle+\frac{t_{1}}{2\tau_{1}}\Big{(}\gamma_{1}-2\widetilde{\mu}_{1}T_{1}q\gamma_{2}+\widetilde{\mu}_{1}^{2}T_{1}^{2}q^{2}+\mu_{0}^{2}\vartheta^{2}-2\mu_{0}\vartheta\gamma_{3}\Big{)},$

is the converging limit of the cost function the scalar formulation in (51),
where the function $T_{3,\lambda}(\cdot,\cdot)$ is given by
$V_{4}(\boldsymbol{t},\boldsymbol{\tau})/\eta$. Before continuing our
analysis, we summarize convexity properties of the cost function of (61) in
the following lemma.

######  Lemma 7 (Strong–convexity of (61)).

Define $\widehat{f}$ as the cost function of the problem in (61) defined in
the feasibility set. Then, $\widehat{f}$ is jointly strongly convex in the
variables $(q,\vartheta,\tau_{1},\tau_{2})$ for fixed feasible
$(t_{1},t_{2})$. Moreover, it is jointly strongly concave in the variables
$(t_{1},t_{2})$ for fixed feasible $(q,\vartheta,\tau_{1},\tau_{2})$.

This result can be proved by observing that the strong convexity parameters in
Lemma 4 are independent of $p$ and that the operations performed after Lemma 4
preserve the convexity properties. Another property of the scalar formulation
is that its set of optimal solutions concentrates around the set of optimal
solutions of the formulation in (61) as summarized in the following lemma.

######  Lemma 8 (Consistency of the Scalar Formulation).

Define $\tau_{p,1}^{\star}$, $\tau_{p,2}^{\star}$, $t_{p,1}^{\star}$,
$t_{p,2}^{\star}$ and $q_{p}^{\star}$ as the optimal solutions of the scalar
formulation given in (51). Additionally, define $\tau_{1}^{\star}$,
$\tau_{2}^{\star}$, $t_{1}^{\star}$, $t_{2}^{\star}$ and $q^{\star}$ as the
optimal solutions of the deterministic optimization problem given in (61).
Therefore, the following convergence in probability holds

$\displaystyle\tau_{p,1}^{\star}\xrightarrow{p\to+\infty}\tau_{1}^{\star},~{}\tau_{p,2}^{\star}\xrightarrow{p\to+\infty}\tau_{2}^{\star},~{}t_{p,1}^{\star}\xrightarrow{p\to+\infty}t_{1}^{\star},$
(62)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}t_{p,2}^{\star}\xrightarrow{p\to+\infty}t_{2}^{\star},~{}q_{p}^{\star}\xrightarrow{p\to+\infty}q^{\star}.$

Moreover, define $\vartheta_{p}^{\star}$ and $\vartheta^{\star}$ as the
optimal solutions of the minimization problems of (51) and (61) over
$\vartheta$ in the feasibility set defined in (27). Then, we also have the
following convergence in probability

$\displaystyle\vartheta_{p}^{\star}\xrightarrow{p\to+\infty}\vartheta^{\star}.$
(63)

The convergence result in (LABEL:sop_conv) follows using [38, Theorem 2.1]. We
can see that all the assumptions in [38, Theorem 2.1] are satisfied by the
formulations in (51) and (61). Moreover, the result in (63) follows using [16,
Proposition 2]. The detailed proof is omitted since it follows similar ideas
as in Proposition $4$ and Proposition $5$ in [16]. Based on [36], we can
further simplify the formulation in (61) by solving the minimization problem
over the variables $q$ and $\vartheta$. Note that the optimal
$\vartheta^{\star}$ satisfies $\vartheta^{\star}=0$ if $\mu_{0}=0$ and
$\vartheta^{\star}=\gamma_{3}/\mu_{0}$ otherwise. Furthermore, the optimal $q$
denoted by $q_{t,\tau}^{\star}$ can be expressed as follows

$\displaystyle
q_{t,\tau}^{\star}=\frac{\gamma_{2}t_{1}\widetilde{\mu}_{1}T_{1}T_{2,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}{\tau_{1}+t_{1}\widetilde{\mu}_{1}^{2}T_{1}^{2}T_{2,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}.$
(64)

Observe that the optimal solutions, $\vartheta^{\star}$ and
$q_{t,\tau}^{\star}$, satisfy the boundedness constraints. Moreover, note that
our results are valid for any bounds that satisfy the results in Lemmas 2, 3
and 6. Now that we obtained the asymptotic scalar optimization problem, it
remains to study the asymptotic behavior of the training and generalization
errors.

####  VI-B4 Asymptotic Analysis of the Training and Generalization Errors

First, the generalization error is given by

$\displaystyle\mathcal{E}_{\text{test}}$
$\displaystyle=\frac{1}{4^{\upsilon}}\mathbb{E}\left[\left(\varphi(\boldsymbol{a}_{\text{new}}^{\top}\boldsymbol{\xi})-\widehat{\varphi}(\widehat{\boldsymbol{w}}^{\top}\sigma(\boldsymbol{F}^{\top}\boldsymbol{a}_{\text{new}}))\right)^{2}\right],$

where $\boldsymbol{a}_{\text{new}}$ is an unseen data sample and
$\widehat{\boldsymbol{w}}$ is the optimal solution of the noisy formulation.
Based on the uniform Gaussian equivalence theorem (uGET), observed and proved
in many earlier papers [17, 18, 15, 33, 16], the asymptotic properties of the
generalization error are equivalent to the asymptotic properties of
$\overline{\mathcal{E}}_{\text{test}}$ defined as follows

$\displaystyle\overline{\mathcal{E}}_{\text{test}}=\frac{1}{4^{\upsilon}}\mathbb{E}\Big{[}$
$\displaystyle\Big{(}\varphi(\boldsymbol{a}_{\text{new}}^{\top}\boldsymbol{\xi})-\widehat{\varphi}(\mu_{0s}\widehat{\vartheta}_{p}+\mu_{1s}\widehat{\boldsymbol{w}}^{\top}\boldsymbol{F}^{\top}\boldsymbol{a}_{\text{new}}+\mu_{2s}\widehat{\boldsymbol{w}}^{\top}\boldsymbol{z})\Big{)}^{2}\Big{]}.$

Here, $\widehat{\boldsymbol{w}}$ and $\widehat{\vartheta}_{p}$ are the optimal
solutions of our primary formulation given in (27). The expectation is taken
over the distribution of the random vector $\boldsymbol{a}_{\text{new}}$, the
random vector $\boldsymbol{z}$ and the possibly random functions
$\varphi(\cdot)$ and $\widehat{\varphi}(\cdot)$, where $\boldsymbol{z}$ is
independent of $\boldsymbol{a}_{\text{new}}$ and drawn from a standard
Gaussian distribution. Moreover, the constants $\mu_{0s}$, $\mu_{1s}$ and
$\mu_{2s}$ are defined as $\mu_{0s}=\mathbb{E}[\sigma(z)]$,
$\mu_{1s}=\mathbb{E}[z\sigma(z)]$ and
$\mu_{2s}^{2}=\mathbb{E}[\sigma(z)^{2}]-\mu_{0s}^{2}-\mu_{1s}^{2}$, where $z$
is a standard Gaussian random vector. Now, consider the following two random
variables

$\displaystyle
g_{1}=\boldsymbol{a}_{\text{new}}^{\top}\boldsymbol{\xi},~{}\text{and}~{}g_{2}=\mu_{0s}\widehat{\vartheta}_{p}+\mu_{1s}\widehat{\boldsymbol{w}}^{\top}\boldsymbol{F}^{\top}\boldsymbol{a}_{\text{new}}+\mu_{2s}\widehat{\boldsymbol{w}}^{\top}\boldsymbol{z}.$

Given the optimal solutions $\widehat{\boldsymbol{w}}$ and
$\widehat{\vartheta}_{p}$, the random variables $g_{1}$ and $g_{2}$ have a
bivariate Gaussian distribution with mean vector
$[0,\mu_{0s}\widehat{\vartheta}_{p}]^{\top}$ and covariance matrix given by

$\displaystyle\boldsymbol{C}_{p}=\begin{bmatrix}\mathinner{\\!\left\lVert\boldsymbol{\xi}\right\rVert}^{2}&\mu_{1s}\boldsymbol{\xi}^{\top}\boldsymbol{F}\widehat{\boldsymbol{w}}\\\
\mu_{1s}\boldsymbol{\xi}^{\top}\boldsymbol{F}\widehat{\boldsymbol{w}}&\mu_{1s}^{2}\mathinner{\\!\left\lVert\boldsymbol{F}\widehat{\boldsymbol{w}}\right\rVert}^{2}+\mu_{2s}^{2}\mathinner{\\!\left\lVert\widehat{\boldsymbol{w}}\right\rVert}^{2}\end{bmatrix}.$

Define the random variables $\widehat{q}_{p}^{\star}$,
$\widehat{\beta}_{p}^{\star}$ and $\widehat{r}_{p}^{\star}$ as follows

$\displaystyle\widehat{q}_{p}^{\star}=\bar{\boldsymbol{v}}^{\top}\widehat{\boldsymbol{w}},~{}\widehat{\beta}_{p}^{\star}=\mathinner{\\!\left\lVert\boldsymbol{F}\widehat{\boldsymbol{w}}\right\rVert}^{2},~{}\text{and}~{}\widehat{r}_{p}^{\star}=\mathinner{\\!\left\lVert\widehat{\boldsymbol{w}}\right\rVert}^{2},$
(65)

where
$\bar{\boldsymbol{v}}=\boldsymbol{v}/\mathinner{\\!\left\lVert\boldsymbol{v}\right\rVert}$
and the vector $\boldsymbol{v}$ is defined as
$\boldsymbol{v}=\boldsymbol{F}^{\top}{\boldsymbol{\xi}}$. Then, the covariance
matrix $\boldsymbol{C}_{p}$ can be expressed as follows

$\displaystyle\boldsymbol{C}_{p}=\begin{bmatrix}1&\mu_{1s}T_{p,1}\widehat{q}^{\star}_{p}\\\
\mu_{1s}T_{p,1}\widehat{q}^{\star}_{p}&\mu_{1s}^{2}\widehat{\beta}_{p}^{\star}+\mu_{2s}^{2}\widehat{r}_{p}^{\star}\end{bmatrix}.$
(66)

Hence, to study the asymptotic properties of the generalization error, it
suffices to study the asymptotic properties of
$\widehat{\vartheta}^{\star}_{p}$, $\widehat{q}^{\star}_{p}$,
$\widehat{\beta}^{\star}_{p}$ and $\widehat{r}^{\star}_{p}$. The following
lemma summarizes the asymptotic properties of our primal formulation given in
(27).

######  Lemma 9 (Primal Consistency).

The random variables $\widehat{\vartheta}^{\star}_{p}$,
$\widehat{q}^{\star}_{p}$, $\widehat{\beta}^{\star}_{p}$ and
$\widehat{r}^{\star}_{p}$ converge in probability as follows

$\displaystyle\widehat{q}_{p}^{\star}\xrightarrow{p\to+\infty}q^{\star},~{}\widehat{\vartheta}_{p}^{\star}\xrightarrow{p\to+\infty}\vartheta^{\star},\widehat{\beta}^{\star}_{p}\xrightarrow{p\to+\infty}\beta^{\star}$
(67)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\widehat{r}_{p}^{\star}\xrightarrow{p\to+\infty}r^{\star}=(q^{\star})^{2}+h^{\prime}(\lambda),$

where $q^{\star}$ and $\vartheta^{\star}$ are optimal solutions of the
deterministic scalar formulation in (61). Moreover, the function $h(\cdot)$
and $\beta^{\star}$ are defined in Theorem 1.

###### Proof.

Note that the analysis in Section VI-B2 shows that the scalar formulation
given in (51) is a simplified version of the multivariate AO formulation given
in (36). Define the random variable $\widetilde{\vartheta}_{p}^{\star}$ as the
optimal solution of the minimization of the problem (36) over $\vartheta$ in
the feasibility set defined in (27). Moreover, define the random variables
$\widetilde{q}_{p}^{\star}$, $\widetilde{\beta}_{p}^{\star}$ and
$\widetilde{r}_{p}^{\star}$ as follows

$\displaystyle\widetilde{q}_{p}^{\star}=\bar{\boldsymbol{v}}^{\top}\widetilde{\boldsymbol{w}},~{}\widetilde{\beta}_{p}^{\star}=\mathinner{\\!\left\lVert\boldsymbol{F}\widetilde{\boldsymbol{w}}\right\rVert}^{2},~{}\text{and}~{}\widetilde{r}_{p}^{\star}=\mathinner{\\!\left\lVert\widetilde{\boldsymbol{w}}\right\rVert}^{2},$
(68)

where $\widetilde{\boldsymbol{w}}$ is the optimal solution of the multivariate
AO formulation given in (36). Based on the decomposition in (45), note that
$\widetilde{\beta}^{\star}_{p}$ satisfies the following expression

$\displaystyle\widetilde{\beta}^{\star}_{p}=\mathinner{\\!\left\lVert\boldsymbol{F}\widetilde{\boldsymbol{w}}\right\rVert}^{2}$
$\displaystyle=(\widetilde{q}^{\star}_{p})^{2}\bar{\boldsymbol{v}}^{\top}\boldsymbol{F}^{\top}\boldsymbol{F}\bar{\boldsymbol{v}}+2\widetilde{q}^{\star}_{p}\bar{\boldsymbol{v}}^{\top}\boldsymbol{F}^{\top}\boldsymbol{F}\boldsymbol{B}^{\perp}_{\boldsymbol{v}}\widetilde{\boldsymbol{r}}^{\star}$
(69)
$\displaystyle+(\widetilde{\boldsymbol{r}}^{\star})^{\top}\bar{\boldsymbol{B}}^{\perp}_{\boldsymbol{v}}\boldsymbol{F}^{\top}\boldsymbol{F}\boldsymbol{B}^{\perp}_{\boldsymbol{v}}\widetilde{\boldsymbol{r}}^{\star},$

where $\widetilde{\boldsymbol{r}}^{\star}$ is defined in (49) and is the
optimal solution of minimizing the function $g(\cdot)$ introduced in (47).
Substituting the expression of $\widetilde{\boldsymbol{r}}^{\star}$ given in
(49), performing the same analysis as in Section VI-B3 and using the
convergence result in (LABEL:sop_conv), it can be shown that the random
quantity $\widetilde{\beta}^{\star}_{p}$ converges in probability to
$\beta^{\star}$ defined in (1). Additionally, observe that
$\widetilde{r}_{p}^{\star}$ can be expressed as follows

$\displaystyle\widetilde{r}_{p}^{\star}=\mathinner{\\!\left\lVert\widetilde{\boldsymbol{w}}\right\rVert}^{2}=(\widetilde{q}^{\star}_{p})^{2}+\mathinner{\\!\left\lVert\widetilde{\boldsymbol{r}}^{\star}\right\rVert}^{2}.$
(70)

Define the function $h_{p}\mathrel{\mathop{\mathchar
58\relax}}\lambda\to-(\widetilde{q}^{\star}_{p})^{2}V_{p,3}(\boldsymbol{t}^{\star}_{p},\boldsymbol{\tau}^{\star}_{p})-V_{p,4}(\boldsymbol{t}^{\star}_{p},\boldsymbol{\tau}^{\star}_{p})$,
where the random functions $V_{p,3}(\cdot,\cdot)$ and $V_{p,4}(\cdot,\cdot)$
are defined in (52) and where
$\boldsymbol{t}^{\star}_{p}=[t_{p,1}^{\star},t_{p,2}^{\star}]^{\top}$ and
$\boldsymbol{\tau}^{\star}_{p}=[\tau_{p,1}^{\star},\tau_{p,2}^{\star}]^{\top}$.
Here,
$\\{t_{p,1}^{\star},t_{p,2}^{\star},\tau_{p,1}^{\star},\tau_{p,2}^{\star}\\}$
are defined in Lemma 8. Given the expression of
$\widetilde{\boldsymbol{r}}^{\star}$ in (49), we can see that
$\widetilde{r}_{p}^{\star}$ can be expressed as follows

$\displaystyle\widetilde{r}_{p}^{\star}=(\widetilde{q}^{\star}_{p})^{2}+h_{p}^{\prime}(\lambda),$
(71)

where the optimal solutions are treated as constants independent of $\lambda$.
Performing the same analysis as in Section VI-B3 and using the convergence
result in (LABEL:sop_conv), it can be shown that the random quantity
$\widetilde{r}_{p}^{\star}$ converges in probability as follows

$\displaystyle\widetilde{r}_{p}^{\star}\xrightarrow{p\to+\infty}r^{\star}=(q^{\star})^{2}+h^{\prime}(\lambda),$
(72)

where $q^{\star}$ is the optimal solution of (61) and the function $h(\cdot)$
is defined in Theorem 1. Given that the scalar formulation given in (51) is a
simplified version of the multivariate AO formulation and based on Lemma 8, we
obtain the following asymptotic properties

$\displaystyle\widetilde{q}_{p}^{\star}\xrightarrow{p\to+\infty}q^{\star},~{}\widetilde{\vartheta}_{p}^{\star}\xrightarrow{p\to+\infty}\vartheta^{\star},$
(73)

where $q^{\star}$ and $\vartheta^{\star}$ are the optimal solutions of the
deterministic scalar formulation in (61). Following a similar analysis as in
[39], we can show that the assumptions in Theorem 3 are all satisfied. The
main idea is to define the set $\mathcal{S}_{p,\epsilon}$ introduced in
Theorem 3 as

$\displaystyle\\{\boldsymbol{w}\mathrel{\mathop{\mathchar
58\relax}}\mathinner{\\!\left\lvert\mathinner{\\!\left\lVert\boldsymbol{F}{\boldsymbol{w}}\right\rVert}^{2}-\beta^{\star}\right\rvert}<\epsilon\\},~{}\text{and}~{}\\{\boldsymbol{w}\mathrel{\mathop{\mathchar
58\relax}}\mathinner{\\!\left\lvert\mathinner{\\!\left\lVert{\boldsymbol{w}}\right\rVert}^{2}-r^{\star}\right\rvert}<\epsilon\\}.$
(74)

Then, use the strong convexity properties of the formulation in
(LABEL:ana_fm8) to prove that the assumptions in Theorem 3 are satisfied. This
means that $\widehat{\vartheta}^{\star}_{p}$, $\widehat{q}^{\star}_{p}$,
$\widehat{\beta}^{\star}_{p}$ and $\widehat{r}^{\star}_{p}$ defined in (65)
concentrates around the same values as $\widetilde{\vartheta}^{\star}_{p}$,
$\widetilde{q}^{\star}_{p}$, $\widetilde{\beta}^{\star}_{p}$ and
$\widetilde{r}^{\star}_{p}$ defined in (68). ∎

Now, to show the convergence of the generalization error in Theorem 1, it
suffices to show that $\overline{\mathcal{E}}_{\text{test}}$ is a continuous
function in $\widehat{\vartheta}^{\star}_{p}$, $\widehat{q}^{\star}_{p}$,
$\widehat{\beta}^{\star}_{p}$ and $\widehat{r}^{\star}_{p}$. Based on
Assumption 3, the functions $\varphi(\cdot)$ and $\widehat{\varphi}(\cdot)$
are square integrable over Gaussian distributions. Moreover, the optimal
solutions $\widehat{\vartheta}^{\star}_{p}$, $\widehat{q}^{\star}_{p}$,
$\widehat{\beta}^{\star}_{p}$ and $\widehat{r}^{\star}_{p}$ are bounded. Based
on Assumption 3 and the continuity under integral sign property [40], the
continuity of $\overline{\mathcal{E}}_{\text{test}}$ follows. These properties
lead to the convergence result given in (20) in Theorem 1. Based on the
analysis in Lemma 9 and Theorem 3, the optimal cost value of the noisy
formulation converges in probability to the optimal cost value of the
deterministic formulation in (61). Combining this result with the asymptotic
property stated in (72) shows the convergence of the training error stated in
Theorem 1.

###  VI-C Large Number of Noise Injections

Note that the analysis in the previous parts studies the properties of the
training and generalization errors corresponding to (3) when $p$ grows to
infinity. In this part, we study the properties of the noisy formulation when
$\ell$ grows to infinity slower than the dimensions $n$, $p$ and $k$. Note
that in this regime and based on the analysis in the previous part, the noisy
formulation converges to the deterministic formulation in (61). Then, the
objective is to analyze the deterministic formulation in (61) when $\ell$
grows to infinity. Here, we note that the bounds on the feasibility sets of
the deterministic formulation in (61) depends on $\ell$ as follows

$\displaystyle
C_{t_{1}}=\sqrt{\ell}\overline{C}_{t_{1}},~{}C_{t_{2}}=\sqrt{\ell}\overline{C}_{t_{2}},~{}c_{\tau_{1}}=\sqrt{\ell}\overline{c}_{\tau_{1}},$
$\displaystyle
c_{\tau_{2}}=\sqrt{\ell}\overline{c}_{\tau_{2}},~{}C_{\tau_{1}}=\sqrt{\ell}\overline{C}_{\tau_{1}},~{}C_{\tau_{2}}=\sqrt{\ell}\overline{C}_{\tau_{2}}.$
(75)

We start our analysis by performing the change of variable
${\tau}_{1}=\tau_{1}/\sqrt{\ell}$, ${\tau}_{2}=\tau_{2}/\sqrt{\ell}$,
${t}_{1}=t_{1}/\sqrt{\ell}$ and ${t}_{2}=t_{2}/\sqrt{\ell}$. This means that
the deterministic scalar optimization problem given in (61) can be expressed
as follows

$\displaystyle\max_{\begin{subarray}{c}0\leq{t_{1}}\leq\overline{C}_{t_{1}}\\\
0\leq{t_{2}}\leq\overline{C}_{t_{2}}\end{subarray}}\min_{\begin{subarray}{c}\overline{c}_{\tau_{1}}\leq{\tau_{1}}\leq\overline{C}_{\tau_{1}}\\\
\overline{c}_{\tau_{2}}\leq{\tau_{2}}\leq\overline{C}_{\tau_{2}}\end{subarray}}~{}\frac{\tau_{1}t_{1}+\tau_{2}t_{2}}{2}-\frac{t_{1}^{2}+t_{2}^{2}}{2}+\frac{(q_{t,\tau}^{\star})^{2}}{2T_{2,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}-\frac{\eta
T_{3,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}{2}$ (76)
$\displaystyle+\frac{t_{1}}{2\tau_{1}}\Big{(}\gamma_{1}-2\widetilde{\mu}_{1}T_{1}q_{t,\tau}^{\star}\gamma_{2}+\widetilde{\mu}_{1}^{2}T_{1}^{2}(q_{t,\tau}^{\star})^{2}+\mu_{0}^{2}(\vartheta^{\star})^{2}-2\mu_{0}\vartheta^{\star}\gamma_{3}\Big{)},$

where $\boldsymbol{t}=[t_{1},t_{2}]^{\top}$ and
$\boldsymbol{\tau}=[\tau_{1},\tau_{2}]^{\top}$. Here, the functions
$T_{2,\lambda}(\cdot,\cdot)$ and $q_{t,\tau}^{\star}$ are the same as the ones
provided in Section IV. Furthermore, the functions
$T_{3,\lambda}(\cdot,\cdot)$ and $g_{\kappa,\lambda}(\cdot,\cdot)$ are given
as follows

$\displaystyle T_{3,\lambda}(\boldsymbol{t},\boldsymbol{\tau})$
$\displaystyle=t_{1}^{2}\mathbb{E}\Big{[}\frac{\widetilde{\mu}_{1}^{2}\kappa+\mu_{2}^{2}}{g_{\kappa,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}\Big{]}+\frac{t_{1}^{2}+t_{2}^{2}}{\ell}\mathbb{E}\Big{[}\frac{\widehat{\mu}_{1}^{2}\kappa+\mu_{3}^{2}}{g_{\kappa,\lambda}(\boldsymbol{t},\boldsymbol{\tau})}\Big{]},$
$\displaystyle g_{\kappa,\lambda}(\boldsymbol{t},\boldsymbol{\tau})$
$\displaystyle=\frac{t_{1}}{\tau_{1}}\Big{(}\widetilde{\mu}_{1}^{2}\kappa+\mu_{2}^{2}\Big{)}+\Big{(}\frac{t_{1}}{\tau_{1}\ell}+\frac{t_{2}(\ell-1)}{\tau_{2}\ell}\Big{)}$
$\displaystyle\times\Big{(}\widehat{\mu}_{1}^{2}\kappa+\mu_{3}^{2}\Big{)}+\lambda.$

Now, we focus on analyzing the formulation in (76) when the number of noise
injections grows to infinity. The following lemma summarizes our main
technical results.

######  Lemma 10 (Large Number of Noise Injections).

When $\ell$ grows to infinity, the asymptotic limit of the formulation in (76)