asymptotic form of the result follows from the limit in the proportional
asymptotics:

$E_{\mu}=\lim_{n_{0},m\to\infty}E_{\Sigma^{*}}=\frac{\sigma_{\epsilon}^{2}\phi}{1-\phi}\mathbb{E}_{\mu}[r/\lambda]\,.$
(S26)

Next we prove Prop. 3.0 using the Harris inequality before proving a slightly
more general result that properly handles the case where $\Sigma$ may have
repeated eigenvalues.

Proof [Proof of Prop. 3.0] Using the Harris inequality,

$E_{\Sigma^{*}_{1}}^{\text{LR}}=\frac{\sigma_{\varepsilon}^{2}n_{0}}{m-n_{0}-1}\frac{\bar{\mathrm{tr}}(\Sigma^{*}_{2})}{n_{0}}\sum_{i=1}^{n_{0}}\frac{r_{i,2}}{\bar{\mathrm{tr}}(\Sigma^{*}_{2})}\frac{r_{i,1}}{r_{i,2}}\frac{1}{\lambda_{i}}\leq\frac{\sigma_{\varepsilon}^{2}n_{0}}{m-n_{0}-1}\frac{\bar{\mathrm{tr}}(\Sigma^{*}_{1})}{\bar{\mathrm{tr}}(\Sigma^{*}_{2})}\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\frac{r_{i,2}}{\lambda_{i}}\leq
E_{\Sigma^{*}_{2}}^{\text{LR}},$ (S27)

since $1/\lambda_{i}$ and ${r_{i,1}}/{r_{i,2}}$ are nonincreasing and
nondecreasing in $i$ respectively.

Prop. 3.0 can be strengthened, and doing so motivates the occurrence of the
conditional expectation in Def. 4.0. Note that we assume that the eigenvalues
of $\Sigma$ are in nondecreasing order, and each eigenvalue has associated to
it two overlap coefficients, $r_{i,1}$ and $r_{i,2}$. Prop. 3.0 assumes that
$r_{i,1}/r_{i,2}$ form a nondecreasing sequence. However, what happens when
$\Sigma$ has repeated eigenvalues? In this case, the ordering of the
$\lambda_{i}$ can be changed, which in turn changes the associated $r_{i,1}$
and $r_{i,2}$. Therefore, the assumption on $r_{i,1}/r_{i,2}$ is too
strong—reordering the repeated eigenvalues might be sufficient to satisfy the
condition even if it is violated for the original ordering. Instead, we can
introduce the conditional expectation to handle this more gracefully.

In the following, we use $\tilde{\lambda}_{1},\ldots,\tilde{\lambda}_{k}$ to
denote the $k$ non-repeated eigenvalues of the training covariance $\Sigma$
and the sets $S_{j}$ for $j\in\\{1,\ldots,k\\}$ to denote the indices of
eigenvalues in $\\{1,\ldots,n_{0}\\}$ associated to the $j$th repeated
eigenvalue of $\Sigma$. Analogously, we define the corresponding non-repeated
overlap coefficients as,

$\displaystyle\tilde{r}_{j}=\frac{1}{\lvert S_{j}\rvert}\sum_{i\in
S_{j}}r_{i},\quad j\in\\{1,\ldots,k\\}$ (S28)

This is equivalent to the conditional expectation discussed in Eq. S22. In
this case, the measure-theoretic definition of _hardness_ in Def. 4.0 becomes
equivalent to stating that the sequence
$\tfrac{\tilde{r}_{j,2}}{\tilde{r}_{j,1}}$ is nondecreasing as $j$ ranges from
$\\{1,\ldots,k\\}$ when $\Sigma^{*}_{2}$ is _harder_ then $\Sigma^{*}_{1}$
(here $\tilde{r}_{j,1}$ and $\tilde{r}_{j,2}$ denote the non-repeated overlap
coefficients of $\Sigma^{*}_{1}$ and $\Sigma^{*}_{2}$ respectively).

###### Proposition D.0

Let $\tilde{r}_{j,1}$ and $\tilde{r}_{j,2}$ denote the non-repeated overlap
coefficients999Recall the definition in Eq. S28. of $\Sigma^{*}_{1}$ and
$\Sigma^{*}_{2}$ relative to $\Sigma$. If
$\mathrm{tr}(\Sigma^{*}_{2})\geq\mathrm{tr}(\Sigma^{*}_{1})$ and the ratios
${\tilde{r}_{j,1}}/{\tilde{r}_{j,2}}$ form a nondecreasing sequence, then
$E_{\Sigma^{*}_{2}}^{\text{LR}}\geq E_{\Sigma^{*}_{1}}^{\text{LR}}$

Proof  From Prop. 3.0 it is enough to show
$\mathrm{tr}(\Sigma^{*}_{2}\Sigma^{-1})\geq\mathrm{tr}(\Sigma^{*}_{1}\Sigma^{-1})$.
First, note that $\sum_{i=1}^{n_{0}}r_{i,1}=\sum_{j=1}^{k}\lvert
S_{j}\rvert\tilde{r}_{j,1}=\mathrm{tr}(\Sigma^{*}_{1})$ and
$\sum_{i=1}^{n_{0}}r_{i,2}=\sum_{j=1}^{k}\lvert
S_{j}\rvert\tilde{r}_{j,2}=\mathrm{tr}(\Sigma^{*}_{2})$. Then,

$\displaystyle\mathrm{tr}(\Sigma^{*}_{1}\Sigma^{-1})$
$\displaystyle=\sum_{i=1}^{n_{0}}\frac{r_{i,1}}{\lambda_{i}}=\sum_{j=1}^{k}\lvert
S_{j}\rvert\frac{\tilde{r}_{j,1}}{\tilde{\lambda}_{j}}=\mathrm{tr}(\Sigma^{*}_{2})\sum_{j=1}^{k}\lvert
S_{j}\rvert\frac{\tilde{r}_{j,2}}{\mathrm{tr}(\Sigma^{*}_{2})}\frac{\tilde{r}_{j,1}}{\tilde{r}_{j,2}}\frac{1}{\tilde{\lambda}_{j}}$
(S29) $\displaystyle\leq\mathrm{tr}(\Sigma^{*}_{2})\left(\sum_{j=1}^{k}\lvert
S_{j}\rvert\frac{\tilde{r}_{j,2}}{\mathrm{tr}(\Sigma^{*}_{2})}\frac{\tilde{r}_{j,1}}{\tilde{r}_{j,2}}\right)\left(\sum_{j=1}^{k}\lvert
S_{j}\rvert\frac{\tilde{r}_{j,2}}{\mathrm{tr}(\Sigma^{*}_{2})}\frac{1}{\tilde{\lambda}_{j}}\right)$
(S30)
$\displaystyle=\frac{\mathrm{tr}(\Sigma^{*}_{1})}{\mathrm{tr}(\Sigma^{*}_{2})}\mathrm{tr}(\Sigma^{*}_{2}\Sigma^{-1})$
(S31) $\displaystyle\leq\mathrm{tr}(\Sigma^{*}_{2}\Sigma^{-1}),$ (S32)

where the inequality is a consequence of the Harris inequality (Grimmett,
1999, Section 2.2): since $\lvert
S_{j}\rvert{\tilde{r}_{j,1}}/{\tilde{\lambda}_{j}}$ is a normalized measured
with respect to $j$, the sequence ${1}/{\tilde{\lambda}_{j}}$ is nonincreasing
in $j$, while the sequence ${\tilde{r}_{j,1}}/{\tilde{r}_{j,2}}$ is
nondecreasing in $j$.

###  D.2 Linear regression limit of random feature regression

In this section, we show that taking appropriate limits of Cor. 5.0 and Cor.
5.0 recovers existing results for linear regression in high-dimensions. We
first examine the ridgeless limit, for which the results fall into two cases:
$\phi<1$, which we discuss in Sec. D.2.1 and connect to Prop. 3.0 from the
main text; and $\phi>1$, which we discuss in Sec. D.2.2 and connect to Hastie
et al. (2019, Cor. 2). We then consider non-zero regularization, and in Sec.
D.2.3 we connect Cor. 5.0 to Wu and Xu (2020, Theorem 1).

To recover these results, we let $\sigma$ approach the identity activation
function and $\psi\to 0$101010The order of these limits does not change the
result, but for concreteness we take take the limit as $\psi\to 0$ first.,
since as more random features are added the kernel concentrate around its
limit. Moreover, since we are taking the limit $\psi\to 0$ and we assume
$\phi>0$, we may assume $\phi>\psi$ in all calculations for simplicity. We
also note that as $\sigma$ approaches the identity, the constants associated
to $\sigma$ approach the following limits:

$\eta,\zeta,\eta_{*},\zeta_{*},\xi\to
1\quad\text{and}\quad\omega,\omega_{*}\to 0.$ (S33)

####  D.2.1 Recovering asymptotic form of Prop. 3.0 for $\phi<1$

From Cor. 5.0, the expression for the total error is

$\displaystyle E_{\mu}$
$\displaystyle=(1-\xi)^{2}s_{*}+2(1-\xi)\xi\mathcal{I}^{*}_{1,1}+\xi^{2}\phi\mathcal{I}^{*}_{1,2}+\xi^{2}\frac{\psi}{\phi-\psi}x(\sigma_{\varepsilon}^{2}+\mathcal{I}_{1,1})(\omega_{*}+\mathcal{I}^{*}_{1,1})$
(S34)
$\displaystyle\qquad+\xi^{2}x\Big{(}1-\frac{x(\omega-\sigma_{\varepsilon}^{2})}{1-x^{2}\mathcal{I}_{2,2}}\Big{)}\mathcal{I}^{*}_{2,2},$
(S35)

where $x=1/(\omega+\mathcal{I}_{1,1})$ since we are assuming $\psi<\phi$.
Taking the limit $\psi\to 0$, the total error $E_{\mu}$ converges to

$(1-\xi)^{2}s_{*}+2(1-\xi)\xi\mathcal{I}^{*}_{1,1}+\xi^{2}\phi\mathcal{I}^{*}_{1,2}+\xi^{2}x\Big{(}1-\frac{x(\omega-\sigma_{\varepsilon}^{2})}{1-x^{2}\mathcal{I}_{2,2}}\Big{)}\mathcal{I}^{*}_{2,2}.$
(S36)

Recall that Prop. 3.0 assumes that $\phi<1$. We begin by analyzing the
solution to the self-consistent equation for $x$ in the ridgeless limit when
the activation function $\sigma$ becomes linear.

###### Lemma D.0

Suppose $0<\phi<1$ and $\psi<\phi$. In the ridgeless limit, the solution $x$
to the self-consistent equation in Cor. 5.0,

$x=\frac{1}{\omega+\mathcal{I}_{1,1}}\,,$ (S37)

satisfies $\lim_{\omega\to 0}x=\infty$.

Proof  From the definition of $\mathcal{I}_{1,1}$, we see

$\displaystyle\mathcal{I}_{1,1}=\phi\mathbb{E}_{\mu}\left[\frac{\lambda}{\phi+\lambda
x}\right]=\frac{\phi}{x}\mathbb{E}_{\mu}\left[\frac{\lambda}{\phi/x+\lambda}\right]\leq\frac{\phi}{x}.$
(S38)

Using Eq. S37, we find $x\geq({1-\phi})/{\omega}$. Taking $\omega\to 0$
completes the proof.

From Lem. D.0, the solution to the self-consistent equation for $x$ diverges,
i.e. $x\to\infty$ as $\sigma$ becomes linear. In this case, by the dominated
convergence theorem, we have that

$\displaystyle\lim_{\omega\to 0}x\mathcal{I}^{*}_{1,1}$
$\displaystyle=\lim_{x\to\infty}\mathbb{E}_{\mu}\left[{\frac{xr}{\phi+x\lambda}}\right]=\phi\mathbb{E}_{\mu}\left[{r/\lambda}\right],$
(S39) $\displaystyle\lim_{\omega\to 0}x^{2}\mathcal{I}_{2,2}$
$\displaystyle=\lim_{x\to\infty}\phi\mathbb{E}_{\mu}\left[{\frac{x^{2}\lambda^{2}}{(\phi+x\lambda)^{2}}}\right]=\phi,$
(S40) $\displaystyle\lim_{\omega\to 0}x^{2}\mathcal{I}^{*}_{2,2}$
$\displaystyle=\lim_{x\to\infty}\phi\mathbb{E}_{\mu}\left[{\frac{x^{2}r\lambda}{(\phi+x\lambda)^{2}}}\right]=\phi\mathbb{E}_{\mu}\left[{r/\lambda}\right],$
(S41) $\displaystyle\lim_{\omega\to 0}\mathcal{I}^{*}_{1,1}$
$\displaystyle=\lim_{x\to\infty}\phi\mathbb{E}_{\mu}\left[{\frac{r}{(\phi+x\lambda)^{2}}}\right]=0,$
(S42) $\displaystyle\text{and}\quad\lim_{\omega\to 0}\mathcal{I}^{*}_{1,2}$
$\displaystyle=\lim_{x\to\infty}\phi\mathbb{E}_{\mu}\left[{\frac{r}{(\phi+x\lambda)^{2}}}\right]=0\,.$
(S43)

As such, the total error in Eq. S36 converges to

$\displaystyle\sigma_{\varepsilon}^{2}\frac{\phi}{1-\phi}\mathbb{E}_{\mu}[r/\lambda]$
(S44)

as $\omega\to 0$ and $\psi\to 0$, as desired.

####  D.2.2 Recovering results from Hastie et al. (2019) when
$\Sigma=\Sigma^{*}$ and $\phi>1$

The minimum-norm solution for under-determined linear regression is the same
as the limiting ridge-regularized solution as the ridge constant converges to
0. As such, studying the ridgeless limit as in Cor. 5.0 allows for a
comparison to prior results on minimum-norm interpolants (Hastie et al.,
2019). To do so, we again take the limit as $\sigma$ becomes linear; however,
in contrast to the previous section, we now assume $\phi>1$ in order to
compare with Hastie et al. (2019, Cor. 2).

As Hastie et al. (2019) examines the setting in which the training and test
covariate distributions are equal, $\Sigma=\Sigma^{*}$, we have that
$\mathcal{I}^{*}_{a,b}=\mathcal{I}_{a,b}$. Using this relation, the total
error Eq. S36 becomes

$E_{\mu}=(1-\xi)^{2}s+2(1-\xi)\xi\mathcal{I}_{1,1}+\xi^{2}\phi\mathcal{I}_{1,2}+\xi^{2}x\Big{(}1-\frac{x(\omega-\sigma_{\varepsilon}^{2})}{1-x^{2}\mathcal{I}_{2,2}}\Big{)}\mathcal{I}_{2,2}.$
(S45)

Next, letting $\sigma$ become linear as in Eq. S33, we obtain
$x={1}/{\mathcal{I}_{1,1}}$ and

$\displaystyle E_{\mu}\,\stackrel{{\scriptstyle\omega\to
0}}{{\longrightarrow}}\,\mathcal{I}_{1,1}+\sigma_{\varepsilon}^{2}\frac{x^{2}\mathcal{I}_{2,2}}{1-x^{2}\mathcal{I}_{2,2}},$
(S46)

where we used the identity Eq. S5.

We now simplify and relate the result for the limiting error from Hastie et
al. (2019, Cor. 2) in the case of isotropic $\beta$ satisfying
$\lVert\beta\rVert_{2}=1$, where the total error is written as111111Hastie et
al. (2019, Cor. 2) provides expressions for a bias and variance term
separately, but the decomposition is defined slightly differently than the one
we use, so we compare directly to the total error. Note that Eq. S47 corrects
a typo in Hastie et al. (2019).

$\displaystyle\frac{1}{\phi^{2}c_{0}(H,\phi)}+\sigma_{\varepsilon}^{2}c_{0}\phi\frac{\int{\frac{s^{2}}{(1+c_{0}\phi
s)^{2}dH(s)}}}{\int{\frac{s}{(1+c_{0}\phi s)^{2}dH(s)}}}\,.$ (S47)

The measure $H$ is the limiting empirical spectral density of the covariance,
which is equivalent to the marginal distribution of $\lambda$, and $c_{0}$
satisfies the equation

$1-\frac{1}{\phi}=\int{\frac{1}{1+c_{0}\phi s}dH(s)}.$ (S48)

Substituting $c_{0}=x/\phi^{2}$ and rearranging terms, we see

$\displaystyle 0$ $\displaystyle=1-\frac{1}{\phi}-\int{\frac{1}{1+c_{0}\phi
s}dH(s)}$ (S49)
$\displaystyle=-\frac{1}{\phi}+1-\int{\frac{1}{1+(x/\phi)s}dH(s)}$ (S50)
$\displaystyle=-\frac{1}{\phi}+x\int{\frac{s}{\phi+xs}dH(s)}$ (S51)
$\displaystyle=-\frac{x}{\phi}\big{(}\frac{1}{x}-\mathcal{I}_{1,1}\big{)},$
(S52)

which is satisfied by $x=1/\mathcal{I}_{1,1}$, implying the two self-
consistent equations are equivalent and validating the identification
$c_{0}=x/\phi^{2}$. Using the same substitution $c_{0}=x/\phi^{2}$ in Eq. S47
gives

$\displaystyle\mathcal{I}_{1,1}+\sigma_{\varepsilon}^{2}\frac{1}{x\phi}\frac{x^{2}\mathcal{I}_{2,2}}{\mathcal{I}_{1,2}}=\mathcal{I}_{1,1}+\sigma_{\varepsilon}^{2}\frac{x^{2}\mathcal{I}_{2,2}}{1-x^{2}\mathcal{I}_{2,2}},$
(S53)

which matches our expression in Eq. S46.

####  D.2.3 Recovering results from Wu and Xu (2020) when $\Sigma=\Sigma^{*}$
and $\gamma>0$

Viewing the expressions in Cor. 5.0 as functions of $\gamma$ and $\omega$, the
total error satisfies the self-consistency relation
$E_{\mu}(\gamma,\omega)=E_{\mu}(\gamma_{\text{eff}},0)$, i.e. the error with
regularizer $\gamma$ and nonlinearity coefficient $\omega$ (from Eq. 17) is
identical to the error with regularizer $\gamma_{\text{eff}}$ and nonlinearity
constant $\omega=0$ (achieved by taking a linear activation function). Given
this relationship, we seek to connect Cor. 5.0 with the results of Wu and Xu
(2020, Theorem 1) for linear ridge regression, in the case where
$\Sigma=\Sigma^{*}$ and their underlying $\beta$ has an isotropic prior and
their ridge-weighting matrix is the identity. This connection establishes that
the total test error in Cor. 5.0 with regularizer $\gamma$ and nonlinearity
constant $\omega$ matches the total test error of linear ridge regression with
regularizer $\gamma_{\text{eff}}$.

To begin, recall in the case that $\Sigma=\Sigma^{*}$,
$\mathcal{I}^{*}_{a,b}=\mathcal{I}_{a,b}$ and $\xi=1$. In this case, our total
test error can be written as:

$\displaystyle
E_{\mu}=\phi\mathcal{I}_{1,2}+\frac{\sigma_{\varepsilon}^{2}+\phi\mathcal{I}_{1,2}}{\gamma_{\text{eff}}+\phi\mathcal{I}_{1,2}}\cdot
x\mathcal{I}_{2,2}.$ (S54)

with $x=1/(\gamma_{\text{eff}}+\mathcal{I}_{1,1})$.

Simplified to the aforementioned setting, the results of Wu and Xu (2020,
Theorem 1) can be written as121212Note this expression subtracts the
contribution of the irreducible label noise on the test point to match our
definition of the test error. Additionally, the generative model in Wu and Xu
(2020) takes the form
$y(\mathbf{x}_{i})=\beta^{\top}\mathbf{x}_{i}/\sqrt{m}+\varepsilon_{i}$
instead of the convention
$y(\mathbf{x}_{i})=\beta^{\top}\mathbf{x}_{i}/\sqrt{n_{0}}+\varepsilon_{i}$ we
take. Hence the term $\int{\frac{\lambda}{(\lambda
m(-\tilde{\gamma})+1)^{2}}dH(\lambda)}$ is rescaled by a relative factor of
$1/\phi$ relative to the comparable term in Wu and Xu (2020, Theorem 1).,

$\displaystyle
E^{\text{LR}}=\frac{m^{\prime}(-\tilde{\gamma})}{m(-\tilde{\gamma})^{2}}\cdot\left(\int{\frac{\lambda}{(\lambda
m(-\tilde{\gamma})+1)^{2}}dH(\lambda)}+\sigma_{\varepsilon}^{2}\right)-\sigma_{\varepsilon}^{2}$
(S55)

where measure $H$ is the limiting empirical spectral density of the
covariance, which is equivalent to the marginal distribution of $\lambda$.
Additionally, $m(z)$ (the Stieltjes transform of the limiting spectral density
corresponding to the Gram matrix $X^{\top}X/m$) and its derivative are given
by Wu and Xu (2020, Theorem 1) as,

$\displaystyle\tilde{\gamma}$
$\displaystyle=\frac{1}{m(-\tilde{\gamma})}-\phi\int{\frac{\lambda}{1+\lambda
m(-\tilde{\gamma})}dH(\lambda)}$ (S56) $\displaystyle 1$
$\displaystyle=\left(\frac{1}{m^{2}(-\tilde{\gamma})}-\phi\int{\frac{\lambda^{2}}{(\lambda\cdot
m(-\tilde{\gamma})+1)^{2}}dH(\lambda)}\right)m^{\prime}(-\tilde{\gamma})\,,$
(S57)

which can be rearranged to read,

$\displaystyle\phi m(-\tilde{\gamma})$
$\displaystyle=\frac{1}{\tilde{\gamma}/\phi+\phi\int{\frac{\lambda}{\phi+\lambda\cdot\phi
m(-\tilde{\gamma})}dH(\lambda)}}$ (S58)
$\displaystyle\phi^{2}m^{\prime}(-\tilde{\gamma})$
$\displaystyle=\frac{1}{\frac{1}{\phi^{2}m^{2}(-\tilde{\gamma})}-\phi\int{\frac{\lambda^{2}}{(\lambda\cdot\phi
m(-\tilde{\gamma})+\phi)^{2}}dH(\lambda)}}\,.$ (S59)

Whereas we define the Gram matrix as $X^{\top}X/n_{0}$, it is defined as
$X^{\top}X/m$ in Wu and Xu (2020), leading to some normalization differences.
To resolve these discrepancies, it is necessary to identify
$m(-\tilde{\gamma})=x/\phi$ and $\tilde{\gamma}=\phi\gamma_{\text{eff}}$ in
order for the regression vector $\hat{\beta}$ to be consistent between the two
conventions. Under these identifications, Eq. S58 can be written as,

$\displaystyle x$
$\displaystyle=\frac{1}{\gamma_{\text{eff}}+\mathcal{I}_{1,1}}\,,$ (S60)

which agrees with the expression in Cor. 5.0, and Eq. S59 can be written as,

$\displaystyle\frac{\partial x}{\partial\gamma_{\text{eff}}}$
$\displaystyle=\frac{1}{1/x^{2}-\mathcal{I}_{2,2}}$ (S61)
$\displaystyle=-\frac{x}{1/x-\mathcal{I}_{1,1}+\phi\mathcal{I}_{1,2}}$ (S62)
$\displaystyle=-\frac{x}{\gamma_{\text{eff}}+\phi\mathcal{I}_{1,2}}\,,$ (S63)

which, using $\tfrac{\partial
x}{\partial\gamma_{\text{eff}}}=\rho\tfrac{\partial x}{\partial\gamma}$ agrees
with the expression in Eq. S829. Next, we use these identifications to
simplify the expression for the test error in Eq. S57,

$\displaystyle E^{\text{LR}}$
$\displaystyle=\frac{\phi^{2}m^{\prime}(-\tilde{\gamma})}{(\phi
m(-\tilde{\gamma}))^{2}}\cdot\left(\phi^{2}\int{\frac{\lambda}{(\lambda\cdot\phi
m(-\tilde{\gamma})+\phi)^{2}}dH(\lambda)}+\sigma_{\varepsilon}^{2}\right)-\sigma_{\varepsilon}^{2}$
(S64)
$\displaystyle=\frac{1}{1-x^{2}\mathcal{I}_{2,2}}\cdot\left(\phi\mathcal{I}_{1,2}+\sigma_{\varepsilon}^{2}\right)-\sigma_{\varepsilon}^{2}\,.$
(S65)

Invoking the identity
$\frac{1}{1-x^{2}\mathcal{I}_{2,2}}=\frac{1/x}{1/x-x\mathcal{I}_{2,2}}=\frac{\gamma_{\text{eff}}+\mathcal{I}_{1,1}}{\gamma_{\text{eff}}+\phi\mathcal{I}_{1,2}}$,
which follows from Eq. S5 and Eq. S60, we can rewrite the error as,

$\displaystyle
E^{\text{LR}}=\frac{\phi\mathcal{I}_{1,2}\cdot(\gamma_{\text{eff}}+\mathcal{I}_{1,1})}{\gamma_{\text{eff}}+\phi\mathcal{I}_{1,2}}+\frac{\sigma_{\varepsilon}^{2}\cdot(\mathcal{I}_{1,1}-\phi\mathcal{I}_{1,2})}{\gamma_{\text{eff}}+\phi\mathcal{I}_{1,2}}\,.$
(S66)

Finally, we invoke Eq. S5 and Eq. S60 once more to obtain,

$\displaystyle
E^{\text{LR}}=\phi\mathcal{I}_{1,2}+\phi\mathcal{I}_{1,2}\cdot\frac{x\mathcal{I}_{2,2}}{\gamma_{\text{eff}}+\phi\mathcal{I}_{1,2}}+\frac{\sigma_{\varepsilon}^{2}x\mathcal{I}_{2,2}}{\gamma_{\text{eff}}+\phi\mathcal{I}_{1,2}}\,,$
(S67)

which matches the expression in Eq. S54, confirming the equivalence of the two
expressions.

##  E Harder shifts increase the bias and the total error

To begin, we state Prop. 5.0 in the general setting where the scales in
$\mu_{1}$ and $\mu_{2}$ need not match. In this setting, the marginals of the
two test distribution may have two scales $s_{1}$ and $s_{2}$ distinct from
the scale of the training covariance $s$, and the parameters in Thm. 5.0 are
subscripted by a $1$ or $2$ to indicate it whether they are induced by
$\mu_{1}$ or $\mu_{2}$. As discussed in Sec. 4 and Sec. 7, additional
constraints on these parameters are necessary to guarantee an ordering of the
bias and error.

######  Proposition E.0 (Restatement of Prop. 5.0)

Consider two LJSDs $\mu_{1}$ and $\mu_{2}$ such that $\mu_{1}\leq\mu_{2}$.
Then, in the setting of Thm. 5.0,

$\displaystyle B_{\mu_{1}}\leq B_{\mu_{2}}$
$\displaystyle\quad\text{if}\quad\xi_{2}\leq\xi_{1}\leq 1$ (S68)
$\displaystyle E_{\mu_{1}}\leq E_{\mu_{2}}$
$\displaystyle\quad\text{if}\quad\xi_{2}\leq\xi_{1}\leq
1\leq\zeta_{2}/\zeta_{1}\leq\eta_{2}/\eta_{1}\quad\text{and}\quad\sigma_{\varepsilon}^{2}\leq\omega\,.$
(S69)

Proof [Proof of Prop. 5.0] Consider two LJSDs $\mu_{1}$ and $\mu_{2}$ such
that $\mu_{1}\leq\mu_{2}$, and assume $\xi_{2}\leq\xi_{1}\leq 1$. Then,

$\displaystyle B_{\mu_{2}}-B_{\mu_{1}}$
$\displaystyle=(1-\xi_{2})^{2}s_{2}+2(1-\xi_{2})\xi_{2}(\mathcal{I}^{*}_{1,1})_{2}+\xi_{2}^{2}\phi(\mathcal{I}^{*}_{1,2})_{2}$
(S70)
$\displaystyle\qquad-(1-\xi_{1})^{2}s_{1}-2(1-\xi_{1})\xi_{1}(\mathcal{I}^{*}_{1,1})_{1}-\xi_{1}^{2}\phi(\mathcal{I}^{*}_{1,2})_{1}$
(S71)
$\displaystyle\geq(1-\xi_{2})^{2}s_{2}+2(1-\xi_{2})\xi_{2}(\mathcal{I}^{*}_{1,1})_{2}+\xi_{2}^{2}\phi(\mathcal{I}^{*}_{1,2})_{2}$
(S72)
$\displaystyle\qquad-(1-\xi_{1})^{2}s_{2}-2(1-\xi_{1})\xi_{1}(\mathcal{I}^{*}_{1,1})_{2}-\xi_{1}^{2}\phi(\mathcal{I}^{*}_{1,2})_{2}$
(S73)
$\displaystyle=\mathbb{E}_{\mu_{2}}\left[r\bigg{(}(1-\xi_{2})+\xi_{2}\frac{\phi}{\phi+x\lambda}\bigg{)}^{2}\right]-\mathbb{E}_{\mu_{2}}\left[r\bigg{(}(1-\xi_{1})+\xi_{1}\frac{\phi}{\phi+x\lambda}\bigg{)}^{2}\right]$
(S74)
$\displaystyle=\mathbb{E}_{\mu_{2}}\left[r\left(\bigg{(}(1-\xi_{2})+\xi_{2}\frac{\phi}{\phi+x\lambda}\bigg{)}^{2}-\bigg{(}(1-\xi_{1})+\xi_{1}\frac{\phi}{\phi+x\lambda}\bigg{)}^{2}\right)\right]$
(S75)
$\displaystyle=\mathbb{E}_{\mu_{2}}\left[r\left((\xi_{1}-\xi_{2})\Big{(}2-(\xi_{1}+\xi_{2})\frac{x\lambda}{\phi+x\lambda}\Big{)}\frac{x\lambda}{\phi+x\lambda}\right)\right]$
(S76) $\displaystyle\geq 0\,,$ (S77)

where we have again used Lem. A.0 and the inequalities in Eq. S107 which
follow from Cor. A.0.

Now, additionally assume that $\sigma_{\varepsilon}^{2}\leq\omega$ and
$1\leq\zeta_{2}/\zeta_{1}\leq\eta_{2}/\eta_{1}$. Note that the latter
condition implies $s_{2}\xi_{2}^{2}\geq s_{1}\xi_{1}^{2}$ and
${\omega}_{2}/s_{2}\geq{\omega}_{1}/s_{1}$. Reorganizing the terms specifying
the limiting variance, we can rewrite the total error as,

$\displaystyle E_{\mu}$ $\displaystyle=B_{\mu}+V_{\mu}$ (S78)
$\displaystyle=(1-\xi)^{2}s_{*}+2(1-\xi)\xi(\mathcal{I}^{*}_{1,1})+\xi^{2}\phi(\mathcal{I}^{*}_{1,2})$
$\displaystyle\quad-\rho\xi^{2}\frac{\psi}{\phi}\frac{\partial
x}{\partial\gamma}\bigg{(}(\sigma_{\varepsilon}^{2}+\mathcal{I}_{1,1})(\omega+\phi\mathcal{I}_{1,2})(\omega_{*}+\mathcal{I}^{*}_{1,1})+\gamma\tau\mathcal{I}_{2,2}(\omega_{*}+\phi\mathcal{I}^{*}_{1,2})$
$\displaystyle\qquad\qquad\qquad\quad+\frac{\phi}{\psi
x}\gamma\bar{\tau}(\sigma_{\varepsilon}^{2}+\phi\mathcal{I}_{1,2})(\mathcal{I}^{*}_{1,1}-\phi\mathcal{I}^{*}_{1,2})\bigg{)}$
(S79)
$\displaystyle\equiv(1-\xi)^{2}s_{*}+2(1-\xi)\xi(\mathcal{I}^{*}_{1,1})+\xi^{2}\big{(}C_{1}\omega_{*}+C_{2}\mathcal{I}^{*}_{1,1}+C_{3}\mathcal{I}^{*}_{1,2}\big{)}\,,$
(S80)

where the $C_{i}\geq 0$ and depend on $\mu$ only through the marginal
$\lambda$ (i.e. they only depend on the training distribution):

$\displaystyle C_{1}$ $\displaystyle=-\rho\frac{\psi}{\phi}\frac{\partial
x}{\partial\gamma}\big{(}(\sigma_{\varepsilon}^{2}+\mathcal{I}_{1,1})(\omega+\phi\mathcal{I}_{1,2})+\gamma\tau\mathcal{I}_{2,2}\big{)}\geq
0$ (S81) $\displaystyle C_{2}$ $\displaystyle=-\rho\xi^{2}\frac{\partial
x}{\partial\gamma}\left(\frac{\psi}{\phi}(\sigma_{\varepsilon}^{2}+\mathcal{I}_{1,1})(\omega+\phi\mathcal{I}_{1,2})+\frac{\gamma\bar{\tau}}{x}(\sigma_{\varepsilon}^{2}+\phi\mathcal{I}_{1,2})\right)\geq
0$ (S82) $\displaystyle C_{3}$ $\displaystyle=-\rho\frac{\partial
x}{\partial\gamma}\left(\psi\gamma\tau\mathcal{I}_{2,2}-\frac{\phi\gamma\bar{\tau}}{x}(\sigma_{\varepsilon}^{2}+\phi\mathcal{I}_{1,2})-\frac{\phi}{\rho\tfrac{\partial
x}{\partial\gamma}}\right)$ (S83) $\displaystyle=-\rho\frac{\partial
x}{\partial\gamma}\left(\psi\gamma\tau\mathcal{I}_{2,2}-\frac{\phi\gamma\bar{\tau}}{x}(\sigma_{\varepsilon}^{2}+\phi\mathcal{I}_{1,2})+\frac{\phi}{\rho
x}(\gamma+\rho\gamma(\tau\psi/\phi+\bar{\tau})(\omega+\phi\mathcal{I}_{1,2}))\right)$
(S84) $\displaystyle=-\rho\gamma\frac{\partial
x}{\partial\gamma}\left(\psi\tau\mathcal{I}_{2,2}+\frac{\phi\bar{\tau}}{x}(\omega-\sigma_{\varepsilon}^{2})+\frac{\phi}{\rho
x}\big{(}1+\rho\tau\frac{\psi}{\phi}(\omega+\phi\mathcal{I}_{1,2})\big{)}\right)$
(S85) $\displaystyle\geq 0\,.$ (S86)

where the inequalities follow from Lem. A.0 and from the assumption
$\sigma_{\varepsilon}^{2}\leq\omega$. It is now straightforward to see that,

$\displaystyle E_{\mu_{2}}-E_{\mu_{1}}$
$\displaystyle=(1-\xi_{2})^{2}s_{2}+2(1-\xi_{2})\xi_{2}(\mathcal{I}^{*}_{1,1})_{2}-(1-\xi_{1})^{2}s_{1}-2(1-\xi_{1})\xi_{1}(\mathcal{I}^{*}_{1,1})_{1}$
(S87)
$\displaystyle\quad+\xi_{2}^{2}\big{(}C_{1}\omega_{2}+C_{2}(\mathcal{I}^{*}_{1,1})_{2}+C_{3}(\mathcal{I}^{*}_{1,2})_{2}\big{)}-\xi_{1}^{2}\big{(}C_{1}\omega_{1}+C_{2}(\mathcal{I}^{*}_{1,1})_{1}+C_{3}(\mathcal{I}^{*}_{1,2})_{1}\big{)}$
$\displaystyle\geq(1-\xi_{2})^{2}s_{2}+2(1-\xi_{2})\xi_{2}(\mathcal{I}^{*}_{1,1})_{2}-(1-\xi_{1})^{2}s_{2}-2(1-\xi_{1})\xi_{1}(\mathcal{I}^{*}_{1,1})_{2}$
(S88)
$\displaystyle\quad+\xi_{2}^{2}\big{(}C_{1}\omega_{2}+C_{2}(\mathcal{I}^{*}_{1,1})_{2}+C_{3}(\mathcal{I}^{*}_{1,2})_{2}\big{)}-\xi_{1}^{2}\big{(}C_{1}\omega_{1}+C_{2}(\mathcal{I}^{*}_{1,1})_{1}+C_{3}(\mathcal{I}^{*}_{1,2})_{1}\big{)}$
$\displaystyle=(\xi_{1}-\xi_{2})\big{(}(2-\xi_{1}-\xi_{2})s_{2}-2(1-\xi_{1}-\xi_{2})(\mathcal{I}^{*}_{1,1})_{2}\big{)}$
(S89)
$\displaystyle\quad+s_{2}\xi_{2}^{2}\big{(}C_{1}\frac{\omega_{2}}{s_{2}}+C_{2}\frac{(\mathcal{I}^{*}_{1,1})_{2}}{s_{2}}+C_{3}\frac{(\mathcal{I}^{*}_{1,2})_{2}}{s_{2}}\big{)}-s_{1}\xi_{1}^{2}\big{(}C_{1}\frac{\omega_{1}}{s_{1}}+C_{2}\frac{(\mathcal{I}^{*}_{1,1})_{1}}{s_{1}}+C_{3}\frac{(\mathcal{I}^{*}_{1,2})_{1}}{s_{1}}\big{)}$
$\displaystyle\geq(\xi_{1}-\xi_{2})\big{(}(2-\xi_{1}-\xi_{2})(\mathcal{I}^{*}_{1,1})_{2}-2(1-\xi_{1}-\xi_{2})(\mathcal{I}^{*}_{1,1})_{2}\big{)}$
(S90)
$\displaystyle\quad+s_{2}\xi_{2}^{2}\big{(}C_{1}\frac{\omega_{2}}{s_{2}}+C_{2}\frac{(\mathcal{I}^{*}_{1,1})_{2}}{s_{2}}+C_{3}\frac{(\mathcal{I}^{*}_{1,2})_{2}}{s_{2}}\big{)}-s_{2}\xi_{2}^{2}\big{(}C_{1}\frac{\omega_{1}}{s_{1}}+C_{2}\frac{(\mathcal{I}^{*}_{1,1})_{1}}{s_{1}}+C_{3}\frac{(\mathcal{I}^{*}_{1,2})_{1}}{s_{1}}\big{)}$
$\displaystyle\geq(\xi_{1}-\xi_{2})(\xi_{1}+\xi_{2})(\mathcal{I}^{*}_{1,1})_{2}$
(S91)
$\displaystyle\quad+s_{2}\xi_{2}^{2}\big{(}C_{1}\frac{\omega_{2}}{s_{2}}+C_{2}\frac{(\mathcal{I}^{*}_{1,1})_{2}}{s_{2}}+C_{3}\frac{(\mathcal{I}^{*}_{1,2})_{2}}{s_{2}}\big{)}-s_{2}\xi_{2}^{2}\big{(}C_{1}\frac{\omega_{2}}{s_{2}}+C_{2}\frac{(\mathcal{I}^{*}_{1,1})_{2}}{s_{2}}+C_{3}\frac{(\mathcal{I}^{*}_{1,2})_{2}}{s_{2}}\big{)}$
$\displaystyle\geq 0\,,$ (S92)

where we have used the inequalities

$\displaystyle\frac{(\mathcal{I}^{*}_{1,1})_{2}}{s_{2}}-\frac{(\mathcal{I}^{*}_{1,1})_{1}}{s_{1}}\geq
0\,,\quad\frac{(\mathcal{I}^{*}_{1,2})_{2}}{s_{2}}-\frac{(\mathcal{I}^{*}_{1,2})_{1}}{s_{1}}\geq
0\,,\quad\text{and}\quad(\mathcal{I}^{*}_{1,1})_{2}\leq s_{2}\,.$ (S93)

The former two inequalities again follow from Cor. A.0.

###  E.1 Pure-scale shifts

######  Proposition E.0 (Restatement of Prop. 5.0)

Consider two LJSDs $\mu_{1}$ and $\mu_{2}$ such that $\mu_{1}\leq\mu_{2}$ and
$\mathbb{E}_{\mu_{1}}[r|\lambda]/\mathbb{E}_{\mu_{1}}[r|\lambda]=\mathbb{E}_{\mu_{1}}[r]/\mathbb{E}_{\mu_{2}}[r]$,
i.e. a _pure-scale shift_. Then, in the setting of Thm. 5.0,

$\displaystyle B_{\mu_{1}}\leq B_{\mu_{2}}$
$\displaystyle\quad\text{if}\quad\xi_{2}\leq\xi_{1}\leq 1$ (S94)
$\displaystyle V_{\mu_{1}}\leq V_{\mu_{2}}$ $\displaystyle\quad\text{if}\quad
1\leq\zeta_{2}/\zeta_{1}\leq\eta_{2}/\eta_{1}$ (S95) $\displaystyle
E_{\mu_{1}}\leq E_{\mu_{2}}$
$\displaystyle\quad\text{if}\quad\xi_{2}\leq\xi_{1}\leq
1\leq\zeta_{2}/\zeta_{1}\leq\eta_{2}/\eta_{1}\quad\text{and}\quad\sigma_{\varepsilon}^{2}\leq\omega\,.$
(S96)

The result for the bias and total error follows immediately from Prop. 5.0.
For the variance, we can reorganize the terms specifying the limiting variance
as,

Proof

$\displaystyle V_{\mu}$
$\displaystyle=-\rho\xi^{2}\frac{\psi}{\phi}\frac{\partial
x}{\partial\gamma}\bigg{(}(\sigma_{\varepsilon}^{2}+\mathcal{I}_{1,1})(\omega+\phi\mathcal{I}_{1,2})(\omega_{*}+\mathcal{I}^{*}_{1,1})+\frac{\phi}{\psi}\gamma\bar{\tau}(\sigma_{\varepsilon}^{2}+\phi\mathcal{I}_{1,2})\mathcal{I}^{*}_{2,2}$
$\displaystyle\qquad\qquad\qquad\quad+\gamma\tau\mathcal{I}_{2,2}(\omega_{*}+\phi\mathcal{I}^{*}_{1,2})\bigg{)}$
(S97) $\displaystyle\equiv
C_{1}(\omega_{*}+\mathcal{I}^{*}_{1,1})+C_{2}\mathcal{I}^{*}_{2,2}+C_{3}(\omega_{*}+\phi\mathcal{I}^{*}_{1,2})\,,$
(S98)

where the $C_{i}\geq 0$ and depend on $\mu$ only through the marginal
$\lambda$ (i.e. they only depend on the training distribution):

$\displaystyle C_{1}$
$\displaystyle=-\rho\xi^{2}\frac{\psi}{\phi}\frac{\partial
x}{\partial\gamma}(\sigma_{\varepsilon}^{2}+\mathcal{I}_{1,1})(\omega+\phi\mathcal{I}_{1,2})\geq
0$ (S99) $\displaystyle C_{2}$ $\displaystyle=-\rho\xi^{2}\frac{\partial
x}{\partial\gamma}\gamma\bar{\tau}(\sigma_{\varepsilon}^{2}+\phi\mathcal{I}_{1,2})\geq
0$ (S100) $\displaystyle C_{3}$
$\displaystyle=-\rho\xi^{2}\frac{\psi}{\phi}\frac{\partial
x}{\partial\gamma}\gamma\tau\mathcal{I}_{2,2}\geq 0\,,$ (S101)

where the inequalities follow directly from Lem. A.0. From this expression, it
is now straightforward to see that,

$\displaystyle V_{\mu_{2}}-V_{\mu_{1}}$
$\displaystyle=\xi_{2}^{2}\left(C_{1}(\omega_{2}+(\mathcal{I}^{*}_{1,1})_{2})+C_{2}(\mathcal{I}^{*}_{2,2})_{2}+C_{2}(\omega_{2}+\phi\mathcal{I}^{*}_{1,2})_{2}\right)$
$\displaystyle\qquad-\xi_{1}^{2}\left(C_{1}(\omega_{1}+(\mathcal{I}^{*}_{1,1})_{1})+C_{2}(\mathcal{I}^{*}_{2,2})_{1}+C_{3}(\omega_{1}+\phi\mathcal{I}^{*}_{1,2})_{1}\right)$
(S102)
$\displaystyle=s_{2}\xi_{2}^{2}\left(C_{1}(\frac{\omega_{2}}{s_{2}}+\frac{(\mathcal{I}^{*}_{1,1})_{2}}{s_{2}})+C_{2}\frac{(\mathcal{I}^{*}_{2,2})_{2}}{s_{2}}+C_{3}(\frac{\omega_{2}}{s_{2}}+\phi\frac{(\mathcal{I}^{*}_{1,2})_{2}}{s_{2}})\right)$
$\displaystyle\qquad-
s_{1}\xi_{1}^{2}\left(C_{1}(\frac{\omega_{1}}{s_{1}}+\frac{(\mathcal{I}^{*}_{1,1})_{1}}{s_{1}})+C_{2}\frac{(\mathcal{I}^{*}_{2,2})_{1}}{s_{1}}+C_{3}(\frac{\omega_{1}}{s_{1}}+\phi\frac{(\mathcal{I}^{*}_{1,2})_{1}}{s_{1}})\right)$
(S103) $\displaystyle\geq
s_{2}\xi_{2}^{2}\left(C_{1}(\frac{\omega_{2}}{s_{2}}+\frac{(\mathcal{I}^{*}_{1,1})_{2}}{s_{2}})+C_{2}\frac{(\mathcal{I}^{*}_{2,2})_{2}}{s_{2}}+C_{3}(\frac{\omega_{2}}{s_{2}}+\phi\frac{(\mathcal{I}^{*}_{1,2})_{2}}{s_{2}})\right)$
$\displaystyle\qquad-
s_{2}\xi_{2}^{2}\left(C_{1}(\frac{\omega_{1}}{s_{1}}+\frac{(\mathcal{I}^{*}_{1,1})_{1}}{s_{1}})+C_{2}\frac{(\mathcal{I}^{*}_{2,2})_{1}}{s_{1}}+C_{3}(\frac{\omega_{1}}{s_{1}}+\phi\frac{(\mathcal{I}^{*}_{1,2})_{1}}{s_{1}})\right)$
(S104) $\displaystyle\geq
s_{2}\xi_{2}^{2}\left(C_{1}(\frac{\omega_{2}}{s_{2}}+\frac{(\mathcal{I}^{*}_{1,1})_{2}}{s_{2}})+C_{2}\frac{(\mathcal{I}^{*}_{2,2})_{2}}{s_{2}}+C_{3}(\frac{\omega_{2}}{s_{2}}+\phi\frac{(\mathcal{I}^{*}_{1,2})_{2}}{s_{2}})\right)$
$\displaystyle\qquad-
s_{2}\xi_{2}^{2}\left(C_{1}(\frac{\omega_{2}}{s_{2}}+\frac{(\mathcal{I}^{*}_{1,1})_{2}}{s_{2}})+C_{2}\frac{(\mathcal{I}^{*}_{2,2})_{2}}{s_{2}}+C_{3}(\frac{\omega_{2}}{s_{2}}+\phi\frac{(\mathcal{I}^{*}_{1,2})_{2}}{s_{2}})\right)$
(S105) $\displaystyle=0\,,$ (S106)

where we have used the inequalities

$\displaystyle\frac{(\mathcal{I}^{*}_{1,1})_{2}}{s_{2}}-\frac{(\mathcal{I}^{*}_{1,1})_{1}}{s_{1}}\geq
0\,,\qquad\frac{(\mathcal{I}^{*}_{2,2})_{2}}{s_{2}}-\frac{(\mathcal{I}^{*}_{2,2})_{1}}{s_{1}}\geq
0\,,\quad\text{and}\quad\frac{(\mathcal{I}^{*}_{1,2})_{2}}{s_{2}}-\frac{(\mathcal{I}^{*}_{1,2})_{1}}{s_{1}}\geq
0\,,$ (S107)

which are immediate consequences of Cor. A.0 and the assumption
$\mathbb{E}_{\mu_{1}}[r|\lambda]/\mathbb{E}_{\mu_{1}}[r|\lambda]=\mathbb{E}_{\mu_{1}}[r]/\mathbb{E}_{\mu_{2}}[r]-s_{1}/s_{2}$.

###  E.2 Incomparable shifts of equal scale

In this section, we state a bound on the bias and total error of shifts that
are incomparable to the null shift $\mu_{\emptyset}$ in terms of the in-
distribution bias and error.

###### Proposition E.0

Let $\mu$ be a LJSD such that
$\mathbb{E}_{\mu}[\lambda]=s=s_{*}=\mathbb{E}_{\mu}[r]$, and define the
functions

$L(x)\mathrel{\mathop{:}}=\inf_{0\leq\lambda\leq
x}\frac{\mathbb{E}_{\mu}[r|\lambda]}{\lambda}\quad\text{and}\quad
U(x)\mathrel{\mathop{:}}=\sup_{0\leq\lambda\leq
x}\frac{\mathbb{E}_{\mu}[r|\lambda]}{\lambda}.$ (S108)

Denote the bias and total error of $\mu$ by $B_{\mu}$ and $E_{\mu}$
respectively, then

$\frac{\mathbb{E}_{\mu_{\text{train}}}\left[{\lambda
L(\lambda)}\right]}{\mathbb{E}_{\mu_{\text{train}}}\left[{\lambda}\right]}B_{\mu_{\emptyset}}\leq
B_{\mu}\leq\frac{\mathbb{E}_{\mu_{\text{train}}}\left[{\lambda
U(\lambda)}\right]}{\mathbb{E}_{\mu_{\text{train}}}\left[{\lambda}\right]}B_{\mu_{\emptyset}},$
(S109)

and further if $\sigma_{\varepsilon}^{2}\leq\omega$, then

$\frac{\mathbb{E}_{\mu_{\text{train}}}\left[{\lambda
L(\lambda)}\right]}{\mathbb{E}_{\mu_{\text{train}}}\left[{\lambda}\right]}E_{\mu_{\emptyset}}\leq
E_{\mu}\leq\frac{\mathbb{E}_{\mu_{\text{train}}}\left[{\lambda
U(\lambda)}\right]}{\mathbb{E}_{\mu_{\text{train}}}\left[{\lambda}\right]}E_{\mu_{\emptyset}}.$
(S110)

Proof  Since $s=s_{*}$, $\xi=1$. Thus, recall
$B_{\mu}=\phi\mathcal{I}^{*}_{1,2}$. Note that the function
$\lambda\mapsto(\phi+x\lambda)^{-2}$ is nonincreasing, thus Lem. A.0 implies
(S109), since $B_{\mu_{\emptyset}}=\phi\mathcal{I}_{1,2}$.

For the total error, recall from Eqs. (S81)-(S86) that
$E_{\mu}=C_{1}\omega+C_{2}\mathcal{I}^{*}_{1,1}+C_{3}\mathcal{I}^{*}_{1,2}$,
where the coefficients $C_{1}$, $C_{2}$, and $C_{3}$ are nonnegative.
Moreover, these coefficients are the same for $\mu_{\emptyset}$ as $s=s_{*}$.
Finally, we can again apply Lem. A.0 as the functions
$\lambda\mapsto(\phi+x\lambda)^{-1}$ and $\lambda\mapsto(\phi+x\lambda)^{-2}$
are nonincreasing, to find Eq. S110.

##  F The benefit of overparameterization

###  F.1 The bias is nonincreasing

We begin by examining the behavior of the bias as a function of the
overparameterization ratio $\phi/\psi$ in a more general setting than
discussed in the main text. In particular, we relax Sec. 4 and allow for
unequal scales, $s\neq s_{*}$. As discussed in Sec. 7, monotonicity of the
bias is not guaranteed if $\xi>1$. As such, we adopt the additional condition
$\xi\leq 1$.

######  Proposition F.0 (Restatement of Prop. 5.0)

In the setting of Thm. 5.0, if $\xi\leq 1$, the bias $B_{\mu}$ is a
nonincreasing function of the overparameterization ratio $\phi/\psi$.

Proof [Proof of Prop. F.0] Recall from Thm. 5.0 that the bias is given by

$B_{\mu}=(1-\xi)^{2}s_{*}+2(1-\xi)\xi\mathcal{I}^{*}_{1,1}+\xi^{2}\phi\mathcal{I}^{*}_{1,2}\,,$
(S111)

where $x$ is the unique positive real root of the self-consistent equation,

$x=\frac{1-\gamma\tau}{\omega+\mathcal{I}_{1,1}}\,.$ (S112)

Differentiating Eq. S111 with respect to $\phi/\psi$ gives,

$\displaystyle\frac{\partial
B_{\mu}}{\partial{(\phi/\psi)}}=-\frac{\psi^{2}}{\phi}\frac{\partial
B_{\mu}}{\partial{\psi}}=2\frac{\psi^{2}}{\phi}\frac{\partial
x}{\partial\psi}\cdot\left((1-\xi)\xi\mathcal{I}^{*}_{2,2}+\xi^{2}\phi\mathcal{I}^{*}_{2,3}\right)\,.$
(S113)

Since by assumption $\xi\leq 1$, and Lem. A.0 gives $\mathcal{I}^{*}_{a,b}\geq
0$, it is sufficient to show $\frac{\partial x}{\partial{\psi}}\leq 0$, which
immediately follows by implicitly differentiating Eq. S112 and simplifying the
expression:

$\frac{\partial x}{\partial\psi}=-\frac{\rho
x\tau(\omega+\mathcal{I}_{1,1})}{\phi\big{(}1+\rho(\bar{\tau}+\frac{\psi}{\phi}\tau)(\omega+\phi\mathcal{I}_{1,2})\big{)}}\leq
0\,.$ (S114)

Therefore we conclude that $\frac{\partial B_{\mu}}{\partial{(\phi/\psi)}}\leq
0$.

###  F.2 The variance is nonincreasing

Next, we turn our attention to the variance. We note that Prop. 5.0 does not
rely on Sec. 4: the variance is a nonincreasing function of the
overparameterization ratio $\phi/\psi$ whenever $\psi<\phi$, regardless of the
magnitude of the overall covariance scales $s$ and $s_{*}$. However, the
proposition only focuses on the ridgeless limit, whereas Fig. 1(d) and Fig.
8(a) show that result may in fact hold for nonzero ridge constant. As the
proof is considerably simpler in the ridgeless limit, we defer the analysis of
the nonzero ridge setting to future work.

######  Proposition F.0 (Restatement of Prop. 5.0)

In the setting of Cor. 5.0 and in the overparameterized regime (i.e.
$\psi<\phi$), the variance $V_{\mu}$ is a nonincreasing function of the
overparameterization ratio $\phi/\psi$.

Proof [Proof of Prop. 5.0] Using the chain rule we have that, $\frac{\partial
V_{\mu}}{\partial{(\phi/\psi)}}=\frac{\partial
V_{\mu}}{\partial{\psi}}\left[\frac{\partial(\phi/\psi)}{\partial\psi}\right]^{-1}=-\frac{\partial
V_{\mu}}{\partial{\psi}}\frac{\psi^{2}}{\phi}$, so it is sufficient to show
that $\frac{\partial V_{\mu}}{\partial{\psi}}\geq 0$.

From Cor. 5.0 in the overparameterized regime, the self-consistent equation
for $x$ reads $x=\tfrac{1}{\omega+\mathcal{I}_{1,1}}$ and is independent of
$\psi$. Therefore, $\partial x/\partial\psi=0$ and the expression for
$\partial V_{\mu}/\partial\psi$ follows directly from Eq. 25,

$\displaystyle\frac{\partial
V_{\mu}}{\partial\psi}=\xi^{2}\frac{\phi}{(\phi-\psi)^{2}}x(\sigma_{\varepsilon}^{2}+\mathcal{I}_{1,1})(\omega_{*}+\mathcal{I}^{*}_{1,1})\geq
0\,,$ (S115)

where the inequality follows from Lem. A.0.

###  F.3 The generalization gap is nonincreasing

Finally, we prove that overparameterization also confers enhanced robustness,
specifically that the generalization gap between shifted and unshifted test
error is a nonincreasing function of the overparameterization ratio in the
overparameterized regime.

######  Proposition F.0 (Restatement of Prop. 5.0)

Consider two LJSDs $\mu_{1}$ and $\mu_{2}$ such that $\mu_{1}\leq\mu_{2}$ and
$1\leq\zeta_{2}/\zeta_{1}\leq\eta_{2}/\eta_{1}$. Then, in the setting of Cor.
5.0 and in the overparameterized regime (i.e. $\psi<\phi$), the generalization
gap $E_{\mu_{2}}-E_{\mu_{1}}$ is a nonincreasing function of the
overparameterization ratio $\phi/\psi$.

Proof [Proof of Prop. 5.0] The chain rule gives
$\frac{\partial(E_{\mu_{2}}-E_{\mu_{1}})}{\partial(\phi/\psi)}=-\frac{\psi^{2}}{\phi}\frac{\partial(E_{\mu_{2}}-E_{\mu_{1}})}{\partial\psi}$
so it is sufficient to show
${\frac{\partial(E_{\mu_{2}}-E_{\mu_{1}})}{\partial\psi}\geq 0}$. To that end,
recall from Cor. 5.0 in the overparameterized regime that the self-consistent
equation for $x$ reads $x=\tfrac{1}{\omega+\mathcal{I}_{1,1}}$ and is
independent of $\psi$. Therefore, $\partial x/\partial\psi=0$, which implies
$\partial B_{\mu_{1}}/\partial\psi=0$ and $\partial
B_{\mu_{2}}/\partial\psi=0$, so that
$\frac{\partial(E_{\mu_{2}}-E_{\mu_{1}})}{\partial\psi}=\frac{\partial(V_{\mu_{2}}-V_{\mu_{1}})}{\partial\psi}$,
the expression for which follows directly from Eq. 25. Finally, note that, as
in Sec. E, the conditions $\mu_{1}\leq\mu_{2}$ and
$1\leq\zeta_{2}/\zeta_{1}\leq\eta_{2}/\eta_{1}$ imply $s_{2}\xi_{2}^{2}\geq
s_{1}\xi_{1}^{2}$ and ${\omega}_{2}/s_{2}\geq{\omega}_{1}/s_{1}$, so we have

$\displaystyle\frac{\partial(E_{\mu_{2}}-E_{\mu_{1}})}{\partial\psi}$
$\displaystyle=\xi_{2}^{2}\left(C_{E}(\omega_{2}+(\mathcal{I}^{*}_{1,1})_{2})\right)-\xi_{1}^{2}\left(C_{E}(\omega_{1}+(\mathcal{I}^{*}_{1,1})_{1}\right)$
(S116)
$\displaystyle=s_{2}\xi_{2}^{2}\left(C_{E}\left({\frac{\omega_{2}}{s_{2}}+\frac{(\mathcal{I}^{*}_{1,1})_{2}}{s_{2}}}\right)\right)-s_{1}\xi_{1}^{2}\left(C_{E}\left({\frac{\omega_{1}}{s_{1}}+\frac{(\mathcal{I}^{*}_{1,1})_{1}}{s_{1}}}\right)\right)$
(S117) $\displaystyle\geq
s_{2}\xi_{2}^{2}\left(C_{E}\left({\frac{\omega_{2}}{s_{2}}+\frac{(\mathcal{I}^{*}_{1,1})_{2}}{s_{2}}}\right)\right)-s_{2}\xi_{2}^{2}\left(C_{E}\left({\frac{\omega_{1}}{s_{1}}+\frac{(\mathcal{I}^{*}_{1,1})_{1}}{s_{1}}}\right)\right)$
(S118) $\displaystyle\geq
s_{2}\xi_{2}^{2}\left(C_{E}\left({\frac{\omega_{2}}{s_{2}}+\frac{(\mathcal{I}^{*}_{1,1})_{2}}{s_{2}}}\right)\right)-s_{2}\xi_{2}^{2}\left(C_{E}\left({\frac{\omega_{2}}{s_{2}}+\frac{(\mathcal{I}^{*}_{1,1})_{2}}{s_{2}}}\right)\right)$
(S119) $\displaystyle=0\,,$ (S120)

where we have introduced the shorthand
$C_{E}=\frac{\phi}{(\phi-\psi)^{2}}x(\sigma_{\varepsilon}^{2}+\mathcal{I}_{1,1})\geq
0$ and in the last inequality we have used Cor. A.0 to argue
$\frac{(\mathcal{I}^{*}_{1,1})_{2}}{s_{2}}\geq\frac{(\mathcal{I}^{*}_{1,1})_{1}}{s_{1}}$.

##  G Linear trends between in-distribution and out-of-distribution
generalization

In Sec. 5.4, we investigated the linear relationship between the shifted and
unshifted test error in the ridgeless, overparameterized regime, arriving at
Prop. 5.0. Here we provide a short proof of that proposition, which we restate
for clarity.

######  Proposition G.0 (Restatement of Prop. 5.0)

Consider two LJSDs $\mu_{1}$ and $\mu_{2}$. In the setting of Cor. 5.0 and in
the overparameterized regime (i.e. $\psi<\phi$), $E_{\mu_{2}}$ is linear in
$E_{\mu_{1}}$, parametrically in $\phi/\psi$:

$E_{\mu_{2}}=E_{0}+\underbrace{\frac{\xi_{2}^{2}}{\xi_{1}^{2}}\left({\frac{\omega_{2}+(\mathcal{I}^{*}_{1,1})_{2}}{\omega_{1}+(\mathcal{I}^{*}_{1,1})_{1}}}\right)}_{\textsc{slope}}E_{\mu_{1}}\,,$
(S121)

where $E_{0}$ and $\textsc{slope}\geq 0$ are constants independent of the
overparameterization ratio $\phi/\psi$. Moreover, if $\mu_{1}\leq\mu_{2}$ and
$1\leq\zeta_{2}/\zeta_{1}\leq\eta_{2}/\eta_{1}$, then $\textsc{slope}\geq 1$.

Proof [Proof of Prop. 5.0] In the overparameterized regime (i.e. $\phi>\psi$)
and in the ridgeless limit, the limiting test error for LJSD $\mu_{1}$ is
given by Cor. 5.0 as

$E_{\mu_{1}}=B_{\mu_{1}}+\xi_{1}^{2}\frac{1}{\phi/\psi-1}x(\sigma_{\varepsilon}^{2}+\mathcal{I}_{1,1})(\omega_{1}+(\mathcal{I}^{*}_{1,1})_{1})+\xi_{1}^{2}x\Big{(}1-\frac{x(\omega-\sigma_{\varepsilon}^{2})}{1-x^{2}\mathcal{I}_{2,2}}\Big{)}(\mathcal{I}^{*}_{2,2})_{1}\,.$
(S122)

Similarly, the limiting test error for LJSD $\mu_{2}$ is given by Cor. 5.0 as

$E_{\mu_{2}}=B_{\mu_{2}}+\xi_{2}^{2}\frac{1}{\phi/\psi-1}x(\sigma_{\varepsilon}^{2}+\mathcal{I}_{1,1})(\omega_{2}+(\mathcal{I}^{*}_{1,1})_{2})+\xi_{2}^{2}x\Big{(}1-\frac{x(\omega-\sigma_{\varepsilon}^{2})}{1-x^{2}\mathcal{I}_{2,2}}\Big{)}(\mathcal{I}^{*}_{2,2})_{2}\,.$
(S123)

Note that in both expressions, the self-consistent equation for $x$ reads
$x=\frac{1}{\omega+\mathcal{I}_{1,1}}$, which has no dependence on $\psi$,
and, as such, $x$, $\mathcal{I}_{a,b}$, $(\mathcal{I}^{*}_{a,b})_{1}$,
$(\mathcal{I}^{*}_{a,b})_{2}$, $B_{\mu_{1}}$, and $B_{\mu_{2}}$ do not depend
on the overparameterization ratio $\phi/\psi$. Hence we can simply eliminate
the quantity $\tfrac{1}{\phi/\psi-1}$ from Eqs. S122, LABEL: and S123 to
obtain

$E_{\mu_{2}}=E_{0}+\frac{\xi_{2}^{2}}{\xi_{1}^{2}}\Big{(}\frac{\omega_{2}+(\mathcal{I}^{*}_{1,1})_{2}}{\omega_{1}+(\mathcal{I}^{*}_{1,1})_{1}}\Big{)}E_{\mu_{1}}\,,$
(S124)

where $E_{0}$ does not depend on the overparameterization ratio $\phi/\psi$
and is given by

$\begin{split}E_{0}&\mathrel{\mathop{:}}=B_{\mu_{2}}-\frac{\xi_{2}^{2}}{\xi_{1}^{2}}\frac{\omega_{2}+(\mathcal{I}^{*}_{1,1})_{2}}{\omega_{1}+(\mathcal{I}^{*}_{1,1})_{1}}B_{\mu_{1}}+\xi_{2}^{2}x\Big{(}1-\frac{x(\omega-\sigma_{\varepsilon}^{2})}{1-x^{2}\mathcal{I}_{2,2}}\Big{)}\Big{(}(\mathcal{I}^{*}_{2,2})_{2}-\frac{\omega_{2}+(\mathcal{I}^{*}_{1,1})_{2}}{\omega_{1}+(\mathcal{I}^{*}_{1,1})_{1}}(\mathcal{I}^{*}_{2,2})_{1}\Big{)}\,.\end{split}$
(S125)

The fact that $\textsc{slope}\geq 0$ follows immediately from the expression
and the inequalities in Lem. A.0. To establish the final conclusion, recall
from Sec. E that the conditions $\mu_{1}\leq\mu_{2}$ and
${1\leq\zeta_{2}/\zeta_{1}\leq\eta_{2}/\eta_{1}}$ imply $s_{2}\xi_{2}^{2}\geq
s_{1}\xi_{1}^{2}$ and ${\omega}_{2}/s_{2}\geq{\omega}_{1}/s_{1}$, so we have

slope
$\displaystyle=\frac{\xi_{2}^{2}}{\xi_{1}^{2}}\frac{\omega_{2}+(\mathcal{I}^{*}_{1,1})_{2}}{\omega_{1}+(\mathcal{I}^{*}_{1,1})_{1}}$
(S126)
$\displaystyle=\frac{s_{2}\xi_{2}^{2}}{s_{1}\xi_{1}^{2}}\frac{\omega_{2}/s_{2}+(\mathcal{I}^{*}_{1,1})_{2}/s_{2}}{\omega_{1}/s_{1}+(\mathcal{I}^{*}_{1,1})_{1}/s_{1}}$
(S127)
$\displaystyle\geq\frac{s_{2}\xi_{2}^{2}}{s_{2}\xi_{2}^{2}}\frac{\omega_{2}/s_{2}+(\mathcal{I}^{*}_{1,1})_{2}/s_{2}}{\omega_{2}/s_{2}+(\mathcal{I}^{*}_{1,1})_{1}/s_{1}}$
(S128)
$\displaystyle\geq\frac{s_{2}\xi_{2}^{2}}{s_{2}\xi_{2}^{2}}\frac{\omega_{2}/s_{2}+(\mathcal{I}^{*}_{1,1})_{2}/s_{2}}{\omega_{2}/s_{2}+(\mathcal{I}^{*}_{1,1})_{2}/s_{2}}$
(S129) $\displaystyle=1\,,$ (S130)

where in the last inequality we have used Eq. S107.

##  H Optimal regularization

As the overparameterization ratio $\phi/\psi$ goes to infinity, the expression
for the test error simplifies and is given by the expressions in Cor. 5.0. In
this case, it is straightforward to derive the value of the regularization
constant that minimizes the test error.

######  Proposition H.0 (Restatement of Prop. 6.0)

In the setting of Cor. 5.0 and under Sec. 4, the optimal regularization
$\gamma^{\text{opt}}_{\mu}$ is independent of $\mu$ and satisfies,

$\displaystyle\gamma^{\text{opt}}_{\mu}\mathrel{\mathop{:}}=\operatorname*{arg\,min}_{\gamma\geq-\rho\omega}E_{\mu}(\gamma)=\rho(\sigma_{\varepsilon}^{2}-\omega)\qquad\text{and}\qquad
E_{\mu}^{\text{opt}}\mathrel{\mathop{:}}=E_{\mu}(\gamma^{\text{opt}}_{\mu})=\mathcal{I}^{*}_{1,1}\,.$
(S131)

Moreover, the optimal error is ordered with respect to shift strength, i.e.
for two LJSDs $\mu_{1}\leq\mu_{2}$, $E_{\mu_{1}}^{\text{opt}}\leq
E_{\mu_{2}}^{\text{opt}}$.

Proof [Proof of Prop. 6.0] From the expressions in Sec. I.7 when
$\psi/\phi=0$, we have for $\gamma_{\text{eff}}=\gamma/\rho+\omega$,

$\displaystyle
x=\frac{1}{\gamma_{\text{eff}}+\mathcal{I}_{1,1}}>0\,,\qquad\frac{\partial
x}{\partial\gamma}=-\frac{\tau}{\gamma_{\text{eff}}+\phi\mathcal{I}_{1,2}}<0\,,\qquad\tau=\frac{x}{\rho}\,,\qquad\text{and}\qquad\bar{\tau}=\frac{1}{\gamma}\,,$
(S132)

where the inequalities follow from the constraint $\gamma\geq-\rho\omega$,
which implies that $\gamma_{\text{eff}}\geq 0$.

Differentiation shows that,

$\displaystyle\frac{\partial^{2}x}{\partial\gamma^{2}}$
$\displaystyle=\frac{2}{x}\left(\left(\frac{\partial
x}{\partial\gamma}\right)^{2}-\rho\frac{\partial
x}{\partial\gamma}\mathcal{I}_{2,3}\right)>0\,.$ (S133)

Under Sec. 4, these expressions and those in Cor. 5.0 and Eq. S5 give,

$\displaystyle E_{\mu}$ $\displaystyle=\mathcal{I}^{*}_{1,1}+\frac{\partial
x}{\partial\gamma}\rho(\gamma_{\text{eff}}-\sigma_{\varepsilon}^{2})\mathcal{I}^{*}_{2,2}\,,$
(S134)

so that,

$\displaystyle\frac{\partial E_{\mu}}{\partial\gamma}$
$\displaystyle=-\frac{\partial
x}{\partial\gamma}\mathcal{I}^{*}_{2,2}+\frac{\partial^{2}x}{\partial\gamma^{2}}\rho(\gamma_{\text{eff}}-\sigma_{\varepsilon}^{2})\mathcal{I}^{*}_{2,2}+\frac{\partial
x}{\partial\gamma}\mathcal{I}^{*}_{2,2}-2\left(\frac{\partial
x}{\partial\gamma}\right)^{2}\rho(\gamma_{\text{eff}}-\sigma_{\varepsilon}^{2})\mathcal{I}^{*}_{3,3}$
(S135)
$\displaystyle=\rho(\gamma_{\text{eff}}-\sigma_{\varepsilon}^{2})\left(\frac{\partial^{2}x}{\partial\gamma^{2}}\mathcal{I}^{*}_{2,2}-2\left(\frac{\partial
x}{\partial\gamma}\right)^{2}\mathcal{I}^{*}_{3,3}\right)$ (S136)
$\displaystyle=\rho(\gamma_{\text{eff}}-\sigma_{\varepsilon}^{2})\left(\frac{\partial^{2}x}{\partial\gamma^{2}}\left(\phi\mathcal{I}^{*}_{2,3}+x\mathcal{I}^{*}_{3,3}\right)-2\left(\frac{\partial
x}{\partial\gamma}\right)^{2}\mathcal{I}^{*}_{3,3}\right)$ (S137)
$\displaystyle=\rho(\gamma_{\text{eff}}-\sigma_{\varepsilon}^{2})\left(\frac{\partial^{2}x}{\partial\gamma^{2}}\phi\mathcal{I}^{*}_{2,3}-2\rho\frac{\partial
x}{\partial\gamma}\mathcal{I}_{2,3}\mathcal{I}^{*}_{3,3}\right)\,.$ (S138)

Noting that the term in parentheses on the last line is positive, we find that
the unique minimizer131313Note that $\partial E_{\mu}/\partial\gamma<0$ for
$\gamma_{\text{eff}}<\sigma^{2}$ and $\partial E_{\mu}/\partial\gamma>0$ for
$\gamma_{\text{eff}}>\sigma^{2}$. of the error is given by
$\gamma_{\text{eff}}=\sigma_{\varepsilon}^{2}$. Eq. S131 then follows
immediately from the definition of $\gamma_{\text{eff}}=\gamma/\rho+\omega$
and from Eq. S134.

Finally, to establish the last statement, namely that
$E_{\mu_{1}}^{\text{opt}}\leq E_{\mu_{2}}^{\text{opt}}$ for
$\mu_{1}\leq\mu_{2}$, simply note that
$E_{\mu}(\gamma^{\text{opt}}_{\mu})=\mathcal{I}^{*}_{1,1}$ and apply Cor. A.0.

##  I Proof of Thm. 5.0

As discussed in Sec. 2, we consider predictive functions $\hat{y}$ defined by
random feature kernel ridge regression,

$\hat{y}(\mathbf{x})\mathrel{\mathop{:}}=YK^{-1}K_{\mathbf{x}}\,$ (S139)

for $K\mathrel{\mathop{:}}=K(X,X)+\gamma I_{m}$,
$K_{\mathbf{x}}\mathrel{\mathop{:}}=K(X,\mathbf{x})$, with $\gamma$ is a ridge
regularization constant. Here

$\displaystyle K$ $\displaystyle=\frac{1}{n_{1}}F^{\top}F+\gamma I_{m}\,,$
(S140)

and, as in the main text, we have introduced the abbreviations
$F\mathrel{\mathop{:}}=\sigma(WX/\sqrt{n_{0}})$. The labels are generated by a
linear function parameterized by $\beta\in\mathbb{R}^{n_{0}}$, whose entries
are drawn independently from $\mathcal{N}(0,1)$, i.e.
$Y=\beta^{\top}X/\sqrt{n_{0}}+\varepsilon$. In this section, we develop the
techniques and detailed calculations needed to determine the high-dimensional
limit of the test error,

$E_{\Sigma^{*}}=\mathbb{E}_{\mathbf{x},\beta}\mathbb{E}[(y(\mathbf{x})-\hat{y}(\mathbf{x}))^{2}]-\sigma_{\varepsilon}^{2}=\mathbb{E}_{\mathbf{x},\beta}\mathbb{E}[(\beta^{\top}\mathbf{x}/\sqrt{n_{0}}-YK^{-1}K_{\mathbf{x}})^{2}]\,,$
(S141)

as well as its decomposition into bias and variance terms,

$E_{\Sigma^{*}}=\underbrace{\mathbb{E}_{\mathbf{x},\beta}\mathbb{E}[(\mathbb{E}[\hat{y}(\mathbf{x})]-y(\mathbf{x}))^{2}]}_{B_{\Sigma^{*}}}+\underbrace{\mathbb{E}_{\mathbf{x},\beta}[\mathbb{V}[\hat{y}(\mathbf{x})]]}_{V_{\Sigma^{*}}}\,.$
(S142)

We may also write the test loss as,

$E_{\Sigma^{*}}=\mathbb{E}_{(\mathbf{x},y)}(y-\hat{y}(\mathbf{x}))^{2}=E_{1}+E_{2}+E_{3}$
(S143)

with

$\displaystyle E_{1}$
$\displaystyle=\mathbb{E}_{(\mathbf{x},\beta,\varepsilon)}y(\mathbf{x})^{2}=\mathbb{E}_{(\mathbf{x},\beta,\varepsilon)}\mathrm{tr}(y(\mathbf{x})y(\mathbf{x})^{\top})$
(S144) $\displaystyle E_{2}$
$\displaystyle=-2\mathbb{E}_{(\mathbf{x},\beta,\varepsilon)}(K_{\mathbf{x}}^{\top}K^{-1}Y)y(\mathbf{x})=-2\mathbb{E}_{(\mathbf{x},\varepsilon)}\mathrm{tr}(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top}y(\mathbf{x}))$
(S145) $\displaystyle E_{3}$
$\displaystyle=\mathbb{E}_{(\mathbf{x},\beta,\varepsilon)}(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top})^{2}=\mathbb{E}_{(\mathbf{x},\varepsilon)}\mathrm{tr}(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top}YK^{-1}K_{\mathbf{x}})\,.$
(S146)

###  I.1 Reducing to the mean-zero case

In this section, we aim to show that we may assume without loss of generality
that the activation function is centered. While many of the ideas in the proof
are similar to those of Mei and Montanari (2021), at the technical level many
distinctions arise owing to our differing assumptions on the data matrix $X$
and weight matrix $W$. In particular, their approach utilizes a decomposition
into Gegenbauer polynomials and exact orthogonality among these orthogonal
basis functions. Since this approach is not feasible in our setting of iid
Gaussian matrices, we develop an alternative proof that we present below.

To proceed, we define centered versions of the the random feature matrices $F$
and $f$ as

$\bar{F}_{ij}\mathrel{\mathop{:}}=F_{ij}-\mathbb{E}_{Z}\sigma\left({\sqrt{\bar{\mathrm{tr}}(\Sigma)}Z}\right)\quad\text{and}\quad\bar{f}_{i}\mathrel{\mathop{:}}=f_{i}-\mathbb{E}_{Z}\sigma\left({\sqrt{\bar{\mathrm{tr}}(\Sigma^{*})}Z}\right)$
(S147)

for all $i$ and $j$ and $Z\sim\mathcal{N}(0,1)$. Note that this centering
operation is $n_{0}$-dependent, but in the limit
$\bar{\mathrm{tr}}(\Sigma),\bar{\mathrm{tr}}(\Sigma^{*})\to s,s_{*}$, which
appear in the limiting self-consistent equation. To show replacing $F$ and $f$
by $\bar{F}$ and $\bar{f}$ does not alter the test error in the limit, we have
to show that the $E_{1}$, $E_{2}$, and $E_{3}$ terms of Equations
(S144)-(S146) are not changed in the limit. Clearly this is true for $E_{1}$
as it contains neither $F$ nor $f$. For $E_{2}$ we must show

$\mathbb{E}\left[{(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top}-\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}Y^{\top})y(\mathbf{x})}\right]\to
0$ (S148)

and for $E_{3}$ we must show

$\mathbb{E}\left[{(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top})^{2}-(\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}Y^{\top})^{2}}\right]\to
0,$ (S149)

where

$\bar{K}\mathrel{\mathop{:}}=\frac{1}{n_{1}}\bar{F}^{\top}\bar{F}+\gamma
I\quad\text{and}\quad\bar{K}_{\mathbf{x}}\mathrel{\mathop{:}}=\frac{1}{n_{1}}\bar{F}^{\top}\bar{f}.$
(S150)

Define the random variable

$\Delta\mathrel{\mathop{:}}=(K_{\mathbf{x}}^{\top}K^{-1}-\bar{K}_{\mathbf{x}}\bar{K}^{-1})Y^{\top}.$
(S151)

To control the typical behavior of $\Delta$, we will define an event
$\mathcal{E}$ (see Def. I.0), and let $\mathbf{1}_{\mathcal{E}}$ be an
indicator for $\mathcal{E}$. In Lem. I.0 we show
$\mathbb{P}[\mathcal{E}^{\mathsf{c}}]\to 0$.

We can argue

$\left\lvert\mathbb{E}\left[{\Delta
y(\mathbf{x})}\right]\right\rvert\leq\left\lvert\mathbb{E}\left[{\mathbf{1}_{\mathcal{E}}\Delta
y(\mathbf{x})}\right]\right\rvert+\left\lvert\mathbb{E}\left[{(1-\mathbf{1}_{\mathcal{E}})\Delta
y(\mathbf{x})}\right]\right\rvert,$ (S152)

for Equation (S148), and

$\displaystyle\left\lvert\mathbb{E}\left[{\Delta
y(\mathbf{x})}\right]\right\rvert$
$\displaystyle\leq\left\lvert\mathbb{E}\left[{\mathbf{1}_{\mathcal{E}}\Delta(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top}+\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}Y^{\top})}\right]\right\rvert$
(S153)
$\displaystyle\quad+\left\lvert\mathbb{E}\left[{(1-\mathbf{1}_{\mathcal{E}})\left({(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top})^{2}-(\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}Y^{\top})^{2}}\right)}\right]\right\rvert,$
(S154)

for Equation (S149). We demonstrate Equations (S152) and (S153) are $o(1)$ in
two steps. The first terms on the right-hand sides of Equations (S152) and
(S153) represent the typical behavior of the random variables. To bound them,
we show that, given $\mathcal{E}$, $\Delta\to 0$ in Lem. I.0. The second terms
on the right-hand sides of Equations (S152) and (S153) represent the atypical
behavior of the random variables. To bound them, we use the Cauchy-Schwarz
inequality and the fact that $\mathbb{P}[\mathcal{E}^{\mathsf{c}}]\to 0$ in
Lem. I.0.

To briefly outline the structure of this section: Sec. I.1.1 proves
concentration of some quadratic forms of the underlying random matrices and
bounds their operator norms with high probability; Sec. I.1.2 controls the
atypical behavior mentioned above; Sec. I.1.3 applies the Schur complement
formula to derive an expression for $\Delta$ that we can bound; Sec. I.1.4
defines an event where the expression for $\Delta$ from the Schur complement
formula can be easily bounded; and Sec. I.1.5 completes the argument by
bounding $\Delta$ on this typical event.

####  I.1.1 Prerequisites for concentration and linearization

Just as in Adlam et al. (2019), the concentration of the linear data kernel
about its expectation is a key ingredient to the analysis. Define the random
variables

$\Upsilon\mathrel{\mathop{:}}=\frac{1}{n_{0}}X^{\top}X\quad\text{and}\quad\upsilon_{*}\mathrel{\mathop{:}}=\lVert\mathbf{x}\rVert_{2}^{2}/n_{0}.$
(S155)

###### Remark I.0

Throughout this section and the next we will use $C$ to denote an arbitrarily
large, $n_{0}$-independent, positive constant, which can increase from line to
line. For example, if $X\leq C$ and $Y\leq C$, we will simply write $XY\leq
C$, since $C^{2}$ is still some $n_{0}$-independent constant. This approach is
valid as such replacements only occur a finite number of times. Similarly, $c$
denotes some arbitrarily small, $n_{0}$-independent, positive constant that
can decrease from line to line.

Then we note

$\mathbb{E}\Upsilon=\bar{\mathrm{tr}}(\Sigma)I_{n_{0}}\quad\text{and}\quad\mathbb{E}\upsilon_{*}=\bar{\mathrm{tr}}(\Sigma^{*}).$
(S156)

For the variance, we see by Sec. 2.2 that

$\mathbb{V}\Upsilon_{ij}=\frac{1+\mathbb{I}(i=j)}{n_{0}}\bar{\mathrm{tr}}(\Sigma)\leq
C/n_{0}\quad\text{and}\quad\mathbb{V}\upsilon_{*}=\frac{2}{n_{0}}\bar{\mathrm{tr}}((\Sigma^{*})^{2})\leq
C/n_{0}.$ (S157)

This motivates the following lemma.

###### Lemma I.0

The event

$\mathcal{E}_{\text{data}}\mathrel{\mathop{:}}=\left\\{{\lVert\Upsilon\rVert_{\infty}\leq
C}\right\\}\cap\left\\{{\left\lvert\upsilon_{*}\\!-\\!\bar{\mathrm{tr}}(\Sigma^{*})\right\rvert\leq
Cn_{0}^{c-1/2}}\right\\}\cap\bigcap_{i,j}\left\\{{\left\lvert\Upsilon_{ij}\\!-\\!\delta_{ij}\bar{\mathrm{tr}}(\Sigma)\right\rvert\leq
Cn_{0}^{c-1/2}}\right\\}$ (S158)

occurs with high-probability, that is, for some positive, $n_{0}$-independent
constant $C$

$\mathbb{P}[\mathcal{E}_{\text{data}}^{\mathsf{c}}]\leq Cn_{0}^{-c}$ (S159)

for any positive, $n_{0}$-independent constant $c$.

Proof  Tighter concentration than is obtained directly from the variance can
be shown by observing that $\Upsilon$ and $\upsilon_{*}$ are equal in
distribution to

$\frac{1}{n_{0}}Z^{T}\Sigma
Z\quad\text{and}\quad\frac{1}{n_{0}}z^{\top}\Sigma^{*}z$ (S160)

for $Z$ an $(n_{0},m)$-dimensional matrix and $z$ an $n_{0}$-dimensional
vector both containing i.i.d. standard Gaussian random variables. Then using
the assumption
$\lVert\Sigma\rVert_{\infty},\lVert\Sigma^{*}\rVert_{\infty}\leq C$ and
applying the results of Erdős et al. (2012, Sec. B) (or similar concentration
results), we find that for some positive constant $C$ that

$\mathbb{P}[\left\lvert\Upsilon_{ij}-\delta_{ij}\bar{\mathrm{tr}}(\Sigma)\right\rvert\geq
Cn_{0}^{c-1/2}]\leq C\exp({-n_{0}^{c}}),$ (S161)

where $c$ is any positive constant. An identical concentration result holds
for $\upsilon_{*}-\bar{\mathrm{tr}}(\Sigma^{*})$. Then, by the union bound, we
see

$\displaystyle\mathbb{P}\left[{\left\\{{\left\lvert\upsilon_{*}-\bar{\mathrm{tr}}(\Sigma^{*})\right\rvert\geq
Cn_{0}^{-c}}\right\\}\cup\bigcup_{i,j}\left\\{{\left\lvert\Upsilon_{ij}-\delta_{ij}\bar{\mathrm{tr}}(\Sigma)\right\rvert\geq
Cn_{0}^{c-1/2}}\right\\}}\right]$ (S162) $\displaystyle\qquad\leq
Cn_{0}^{2}C\exp\left({-n_{0}^{c}}\right)$ (S163) $\displaystyle\qquad\leq
C\exp(-n_{0}^{c}).$ (S164)

The operator norm of $\Upsilon$ can also be bounded probabilistically. For
some $n_{0}$-independent constant $C$,

$\lVert\Upsilon\rVert_{\infty}\leq\lVert\Sigma\rVert_{\infty}\lVert
Z/\sqrt{n_{0}}\rVert_{\infty}^{2}\leq
C\lVert\Sigma\rVert_{\infty}\left({1+\sqrt{\frac{m}{n_{0}}}+\sqrt{\frac{\log(2/\delta)}{n_{0}}}}\right)\leq
C$ (S165)

with probability at least $1-\delta$, where $Z$ is i.i.d. Gaussian (see for
example Vershynin, 2018, Theorem 4.4.5). Thus, setting $\delta=n_{0}^{-c}$, we
see $\mathbb{P}[\lVert\Upsilon\rVert_{\infty}\geq C]\leq n_{0}^{-c}$. Applying
the union bound again completes the proof.

###### Lemma I.0

The (squared) expected operator norm of $\Upsilon$ is also bounded:

$\mathbb{E}\lVert\Upsilon\rVert_{\infty}^{2}\leq C$ (S166)

for some positive, $n_{0}$-independent constant $C$.

Proof  This result follows from a tail integration argument. That is using
Vershynin (2018, Theorem 4.4.5)) we have that

$\lVert\Upsilon\rVert_{\infty}\leq\lVert\Sigma\rVert_{\infty}\lVert
Z/\sqrt{n_{0}}\rVert_{\infty}^{2}\leq
C\lVert\Sigma\rVert_{\infty}\left({1+\sqrt{\frac{m}{n_{0}}}+u}\right)$ (S167)

with probability at least $1-2\exp(-n_{0}u^{2})$, where $Z$ is i.i.d.
Gaussian. Hence we have that $\mathbb{P}[\lVert\Upsilon\rVert_{\infty}\geq
C\lVert\Sigma\rVert_{\infty}\left({1+\sqrt{\frac{m}{n_{0}}}+u}\right)]\leq
2\exp(-n_{0}u^{2}).$ Now by the tail integration identity,

$\displaystyle\mathbb{E}[\lVert\Upsilon\rVert_{\infty}^{2}]$
$\displaystyle=2\int_{0}^{\infty}u\mathbb{P}[\lVert\Upsilon\rVert_{\infty}\geq
u]du$ (S168) $\displaystyle\leq
C^{2}\lVert\Sigma\rVert_{\infty}^{2}\\!\left({1\\!+\\!\sqrt{\frac{m}{n_{0}}}}\right)^{2}\\!+\\!C^{2}\lVert\Sigma\rVert_{\infty}^{2}\int_{0}^{\infty}\\!u\mathbb{P}\left[{\lVert\Upsilon\rVert_{\infty}\\!\geq\\!C\lVert\Sigma\rVert_{\infty}\\!\left({1\\!+\\!\sqrt{\frac{m}{n_{0}}}\\!+\\!u}\right)}\right]du$
(S169) $\displaystyle\leq
C^{2}\lVert\Sigma\rVert_{\infty}^{2}\left({1+\sqrt{\frac{m}{n_{0}}}}\right)^{2}+C^{2}\lVert\Sigma\rVert_{\infty}^{2}\cdot\frac{1}{n_{0}}$
(S170) $\displaystyle\leq C.$ (S171)

by integrating the Gaussian tail.

Next, we state similar results for the matrix
$\tilde{\Upsilon}\mathrel{\mathop{:}}=W\Sigma^{*}W^{\top}/n_{0}$ as were
obtained for $\Upsilon$ above.

###### Lemma I.0

The event

$\mathcal{E}_{W}\mathrel{\mathop{:}}=\left\\{{\left\lVert\tilde{\Upsilon}\right\rVert_{\infty}\leq
C}\right\\}\cap\bigcap_{i,j}\left\\{{\left\lvert\tilde{\Upsilon}_{ij}\\!-\\!\delta_{ij}\bar{\mathrm{tr}}(\Sigma^{*})\right\rvert\leq
Cn_{0}^{c-1/2}}\right\\}$ (S172)

occurs with high-probability, that is, for some positive, $n_{0}$-independent
constant $C$

$\mathbb{P}[\mathcal{E}_{W}^{\mathsf{c}}]\leq Cn_{0}^{-c}$ (S173)

for any positive, $n_{0}$-independent constant $c$.

Proof  The claims follow identically to Lem. I.0.

###### Lemma I.0

Moreover on the event $\mathcal{E}_{W}$, we have

$\mathbb{E}_{\mathbf{x}}\bar{\sigma}(W\mathbf{x}/\sqrt{n_{0}})_{i}\leq
Cn_{0}^{c-1/2}$ (S174)

and

$\mathbb{E}_{\mathbf{x}}\bar{\sigma}(W\mathbf{x}/\sqrt{n_{0}})\bar{\sigma}(W\mathbf{x}/\sqrt{n_{0}})^{\top}=\zeta_{*}\tilde{\Upsilon}+(\eta_{*}-\zeta_{*})I_{m}+\Delta,$
(S175)

where $\lVert\Delta\rVert_{\infty}\leq Cn_{0}^{c}$.

Proof  Note that conditional on $W$,
$W\mathbf{x}/\sqrt{n_{0}}\sim\mathcal{N}(0,\tilde{\Upsilon})$. Using a Taylor
expansion in
${\epsilon\mathrel{\mathop{:}}=\tilde{\Upsilon}_{ii}-\bar{\mathrm{tr}}(\Sigma*)}$
below with Sec. 2.2 shows that elementwise

$\displaystyle\left\lvert\mathbb{E}_{\mathbf{x}}\bar{\sigma}(W\mathbf{x}/\sqrt{n_{0}})_{i}\right\rvert$
$\displaystyle=\left\lvert\mathbb{E}_{Z}\bar{\sigma}\left({\sqrt{\tilde{\Upsilon}_{ii}}Z}\right)\right\rvert$
(S176)
$\displaystyle=\left\lvert\mathbb{E}_{Z}{\sigma}\left({\sqrt{\epsilon+\bar{\mathrm{tr}}(\Sigma^{*})}Z}\right)-\mathbb{E}_{Z}{\sigma}\left({\sqrt{\bar{\mathrm{tr}}(\Sigma^{*})}Z}\right)\right\rvert$
(S177) $\displaystyle\leq C\epsilon$ (S178) $\displaystyle\leq
Cn_{0}^{c-1/2},$ (S179)

where $Z$ is a standard Gaussian. See Adlam et al. (2019) for additional
details. This allows us to assume that $\sigma(W\mathbf{x})$ is centered as a
conditional expectation over $\mathbf{x}$ at the expense of an elementwise
$Cn_{0}^{c-1/2}$ error. Denoting
$\bar{f}\mathrel{\mathop{:}}=\bar{\sigma}(W\mathbf{x}/\sqrt{n_{0}})$, we note
this can be extended to the matrix $\mathbb{E}\bar{f}\bar{f}^{\top}$. We see

$\mathbb{E}_{\mathbf{x}}\bar{f}\bar{f}^{\top}=\mathbb{E}_{\mathbf{x}}[(\bar{f}-\mathbb{E}_{\mathbf{x}}\bar{f})(\bar{f}-\mathbb{E}_{\mathbf{x}}\bar{f})^{\top}]+(\mathbb{E}_{\mathbf{x}}\bar{f})(\mathbb{E}_{\mathbf{x}}\bar{f})^{\top}.$
(S180)

However, the second term on the right-hand side of Equation (S180) is small is
operator norm, since it is rank-1 and

$\left\lVert(\mathbb{E}_{\mathbf{x}}\bar{f})(\mathbb{E}_{\mathbf{x}}\bar{f})^{\top}\right\rVert_{\infty}=\left\lVert\mathbb{E}_{\mathbf{x}}\bar{f}\right\rVert_{2}^{2}\leq
Cn_{0}^{c}$ (S181)

by Equation (S179). Next, we see for $i\neq j$ that

$\mathbb{E}_{\mathbf{x}}\bar{f}_{i}\bar{f}_{j}=\mathbb{E}_{Z_{1},Z_{2}}\bar{\sigma}(c_{1}Z_{1})\bar{\sigma}(c_{12}Z_{1}+c_{2}Z_{2})$
(S182)

for constants $c_{1}$, $c_{2}$, and $c_{12}$ depending on $W$ and
$\Sigma^{*}$, where $Z_{1}$ and $Z_{2}$ are independent, standard Gaussians.
Again, see Adlam et al. (2019) for details. The $i=j$ terms can be handled
similarly. Moreover,
${c_{12}=\tilde{\Upsilon}_{ij}/\sqrt{\tilde{\Upsilon}_{ii}}\leq
Cn_{0}^{c-1/2}}$, so Taylor expanding in $c_{12}Z_{1}$, we find

$\mathbb{E}_{\mathbf{x}}\bar{f}\bar{f}^{\top}=\zeta_{*}\Upsilon+(\eta_{*}-\zeta_{*})I+\tilde{\Delta}+(\mathbb{E}_{\mathbf{x}}\bar{f})(\mathbb{E}_{\mathbf{x}}\bar{f})^{\top}.$
(S183)

As a consequence of the Taylor expansion and Sec. 2.2, $\tilde{\Delta}$ has
entries that are bounded by $Cn_{0}^{c-1}$ in absolute value (Adlam et al.,
2019). The final conclusion for
${\Delta=\tilde{\Delta}+(\mathbb{E}_{\mathbf{x}}\bar{f})(\mathbb{E}_{\mathbf{x}}\bar{f})^{\top}}$
follows by upper bounding
$\lVert\tilde{\Delta}\rVert_{\infty}\leq\lVert\tilde{\Delta}\rVert_{F}\leq
Cn_{0}^{c}$ and using the previous bound in Equation (S181).

Finally we include a bound on the operator norm of the random feature matrix.
The argument follows Mei and Montanari (2021, Lemma C.3) closely.

###### Lemma I.0

Under Sec. 2.2, the event

$\mathcal{E}_{F}\mathrel{\mathop{:}}=\left\\{{\left\lVert\bar{F}/\sqrt{n_{0}}\right\rVert_{\infty}\leq
Cn_{0}^{c}}\right\\}$ (S184)

occurs with high-probability, that is, for some positive, $n_{0}$-independent
constant $C$

$\mathbb{P}[\mathcal{E}_{F}^{\mathsf{c}}]\leq Cn_{0}^{-c}$ (S185)

for any positive, $n_{0}$-independent constant $c<10$.

Proof  Consider the matrix

$\displaystyle\bar{R}_{ij}=1_{i\neq
j}\bar{\sigma}(z_{i}^{\top}\Sigma^{1/2}z_{j}/\sqrt{n_{1}})/\sqrt{n_{1}}$
(S186)

for $z_{i}=\Sigma^{-1/2}X_{i}$ (i.e. columns of $\Sigma^{-1/2}X$) for $1\leq
i\leq m$ and $z_{m+i}=W_{i}$ (i.e. rows of $W$) for $1\leq i\leq n_{0}$. By
construction, we note that $\bar{F}/\sqrt{n_{1}}$ is a minor of $\bar{R}$.
Thus, bounding the operator norm of $\bar{R}$ suffices to bound the operator
norm of $\bar{F}$. Moreover, $z_{i}$ are independent Gaussians distributed as
$\mathcal{N}(0,I_{n_{0}})$.

Next we show that the entries to the activation function cannot be too large
with high probability. For convenience throughout we define $M=m+n_{0}$ and
let $\Omega$ be a symmetric matrix with $\lVert\Omega\rVert_{\infty}\leq C$.
Note that, for $i\neq j$, the random variable $z_{i}^{\top}\Omega
z_{j}/\sqrt{n_{0}}=\sum_{k=1}^{n_{0}}\lambda_{k}(\Omega)a_{k}b_{k}/\sqrt{n_{0}}$
where $a_{k}$ and $b_{k}$ are i.i.d. from $\sim\mathcal{N}(0,1)$. So the
moment generating function can be bounded as,

$\displaystyle\mathbb{E}\left[{\exp\left({\sum_{i=1}^{n_{0}}t\lambda_{i}(\Omega)a_{i}b_{i}/\sqrt{n_{0}}}\right)}\right]$
$\displaystyle=\prod_{i=1}^{n_{0}}\mathbb{E}[\exp(a_{i}\cdot
t\lambda_{i}(\Omega)b_{i}/\sqrt{n_{0}})]$ (S187)
$\displaystyle=\prod_{i=1}^{n_{0}}\mathbb{E}[\exp(t^{2}b_{i}^{2}\lambda_{i}(\Omega)^{2}/(2n_{0})]$
(S188)
$\displaystyle\leq\prod_{i=1}^{n_{0}}\mathbb{E}[\exp(t^{2}b_{i}^{2}C^{2}/(2n_{0})]$
(S189) $\displaystyle\leq\mathbb{E}[\exp(t^{2}b_{i}^{2}C^{2}/2)]$ (S190)
$\displaystyle=\frac{1}{\sqrt{1-C^{2}t^{2}}}$ (S191)

for $\lvert t\rvert\leq\frac{1}{C}$. Further
$\frac{1}{\sqrt{1-C^{2}t^{2}}}\leq\exp(C^{2}t^{2})$ for $\lvert
t\rvert\leq\frac{1}{2C}$. This establishes the original random variable is
subexponential. Thus this random variables satisfies

$\displaystyle\mathbb{P}[\lvert z_{i}^{\top}\Omega
z_{j}/\sqrt{n_{0}}\rvert\geq u]\leq
2\max\\{\exp(-u^{2}/2C^{2}),\exp(-u/2C)\\}$ (S192)

(see for example Wainwright, 2019, Proposition 2.9). Define the event

$\mathcal{G}\mathrel{\mathop{:}}=\bigcap_{1\leq i<j\leq
M}\left\\{{\left\lvert\frac{z_{i}^{\top}\Sigma^{1/2}z_{j}}{\sqrt{n_{0}}}\right\rvert\geq
C\sqrt{\log M}}\right\\}\cap\bigcap_{1\leq i<j\leq
M}\left\\{{\left\lvert\frac{z_{i}^{\top}\Sigma
z_{j}}{\sqrt{n_{0}}}\right\rvert\geq C^{2}\sqrt{\log M}}\right\\}.$ (S193)

Then by a union bound and Markov’s inequality, we have that

$\mathbb{P}[\mathcal{G}^{\mathsf{c}}]\leq{4}/{M^{20}}$ (S194)

for large enough $C$.

We now define a modified version of $\bar{\sigma}$, that is the same up to a
constant factor on $\mathcal{G}$ but is truncated outside of $\mathcal{G}$.
Let $\bar{u}=C\sqrt{\log M}$ and taking the constants from Sec. 2.2, we see

$\tilde{\sigma}(u)\mathrel{\mathop{:}}=\begin{cases}\bar{\sigma}(u)\exp(-c_{1}\lvert\bar{u}\rvert)/c_{0}&\text{for
}\lvert u\rvert\leq\bar{u}\\\
\bar{\sigma}(\bar{u})\exp(-c_{1}\lvert\bar{u}\rvert)/c_{0}&\text{for
}u>\bar{u}\\\
\bar{\sigma}(-\bar{u})\exp(-c_{1}\lvert\bar{u}\rvert)/c_{0}&\text{for
}u<-\bar{u}\end{cases}.$ (S195)

Just as we did with $\sigma$ in Equation (S147), we center $\tilde{\sigma}$
with its mean $\tilde{a}$. Note $\tilde{\sigma}$ is a 1-bounded and
1-Lipschitz function so $\lvert\tilde{a}\rvert\leq 1$.

Now consider the matrix

$\displaystyle\tilde{R}_{ij}\mathrel{\mathop{:}}=1_{i\neq
j}(\tilde{\sigma}(z_{i}^{\top}\Sigma^{1/2}z_{j}/\sqrt{n_{0}})-\tilde{a})/\sqrt{n_{1}}.$
(S196)

By controlling the operator norm of $\tilde{R}$, we can control the operator
norm of $\bar{R}$ at the end of the proof. By a covering argument it suffices
to control

$\displaystyle\lVert\tilde{R}\rVert_{\infty}\leq\max_{v\in
S}10\underbrace{\lvert v^{\top}\tilde{R}v\rvert}_{F_{v}(Z)},$ (S197)

where $S$ is a $1/4$-covering of the $M$-dimensional sphere with cardinality
$\exp(cM)$ and the matrix $Z$ is of all the variables $z_{i}$.

We now seek to apply Deshpande and Montanari (2016, Lemma 9, Lemma 20). To
this end we wish to show that $F_{v}(Z)$ is Lipschitz in $Z$, that is, if we
define $\mathcal{Z}\mathrel{\mathop{:}}=\sqrt{t}Z+\sqrt{1-t}Z^{\prime}$ for
($Z,Z^{\prime})\in\mathcal{G}\times\mathcal{G}$, we need to show that

$\displaystyle\max_{v\in S}\max_{t\in[0,1]}\left\lVert\nabla
F_{v}(\sqrt{t}Z+\sqrt{1-t}Z^{\prime})\right\rVert_{F}\leq L$ (S198)

for suitable $L$. Consider the gradient with respect to a column of
$\mathcal{Z}$ (denoted by $\zeta_{l}$):

$\displaystyle\nabla_{\zeta_{l}}F_{v}(\mathcal{Z})$
$\displaystyle=2\frac{v_{l}}{\sqrt{n_{0}n_{1}}}\sum_{i\neq
l}\Sigma^{1/2}\zeta_{i}\underbrace{v_{i}\tilde{\sigma}^{\prime}(\zeta_{i}^{\top}\Sigma^{1/2}\zeta_{l}^{\top}/\sqrt{n_{0}})}_{\xi_{i}}$
(S199)
$\displaystyle=2\frac{v_{l}}{\sqrt{n_{0}n_{1}}}\Sigma^{1/2}\mathcal{Z}\mathbf{\xi},$
(S200)

where $\mathbf{\xi}$ is the vector with coordinates $\xi_{i}$ except at $l$
where it is zero. Continuing, on the set $\mathcal{G}$, we find

$\displaystyle\lVert\nabla_{w_{l}}F_{v}(W)\rVert_{2}^{2}$ $\displaystyle\leq
C^{2}v_{l}^{2}/n_{0}^{2}\sum_{i\neq l,j\neq
l}\left\lvert\xi_{i}\xi_{j}w_{i}^{\top}\Sigma w_{j}\right\rvert$ (S201)
$\displaystyle\leq Cv_{l}^{2}\sqrt{\log M}/n_{0}^{2}\sum_{i}\sum_{j\neq
i}\lvert\xi_{i}\xi_{j}\rvert$ (S202) $\displaystyle\leq Cv_{l}^{2}\sqrt{\log
M}/n_{0}^{2}\sum_{i}\sum_{j\neq i}(\xi_{i}^{2}+\xi_{j}^{2})$ (S203)
$\displaystyle\leq Cv_{l}^{2}\sqrt{\log M}/n_{0}$ (S204)

where the last line follows since
$\lvert\tilde{\sigma}^{\prime}(\cdot)\rvert_{2}\leq 1$ implies that
$\lVert\xi\rVert\leq 1$. So finally, we obtain

$\displaystyle\lVert\nabla_{W}F_{v}(W)\rVert_{F}^{2}\leq C\sqrt{\log
M}/n_{0}=:L^{2}.$ (S205)

Now Deshpande and Montanari (2016, Lemma 9) shows for a constant $D$ that

$\displaystyle\mathbb{P}[{\lvert v^{\top}\tilde{R}v\rvert}>D]\leq
C\exp\left(Cn_{0}-\frac{D^{2}}{L^{2}}\right)+\frac{C}{D^{2}}\mathbb{E}[\max_{v\in
S}(F_{v}(Z)-F_{v}(Z^{\prime}))^{2}\cdot\mathbf{1}_{\mathcal{G}^{\mathsf{c}}}].$
(S206)

We use a crude bound on the complement. First note that

$\max_{v\in S}(F_{v}(Z)-F_{v}(Z^{\prime}))^{2}\leq
C\lVert\tilde{\sigma}(Z^{\top}\Sigma^{1/2}Z/\sqrt{n_{0}})\rVert_{F}^{2}\leq
C\lVert Z^{\top}Z/\sqrt{n_{0}}\rVert_{F}^{2}+Cn_{0}^{2},$ (S207)

since removing the 1-bounded-Lipschitz $\tilde{\sigma}(\cdot)$ can be done by
centering each activation with $\tilde{\sigma}(0)$. Then,

$\displaystyle\mathbb{E}\left[{\max_{v\in
S}(F_{v}(Z)-F_{v}(Z^{\prime}))^{2}\cdot\mathbf{1}_{\mathcal{G}^{\mathsf{c}}}}\right]\leq\sqrt{\frac{C}{n_{1}}(\mathbb{E}[\lVert
Z^{\top}\Sigma^{1/2}Z/\sqrt{n_{0}}\rVert_{F}^{4}]+n_{0}^{4})\mathbb{P}[\mathcal{G}^{\mathsf{c}}]}.$
(S208)

A short computation shows that $\mathbb{E}[\lVert
Z^{\top}\Sigma^{1/2}Z/\sqrt{n_{0}}\rVert_{F}^{4}]\leq Cn_{0}^{8}$. Combining
all terms then shows that

$\displaystyle\mathbb{P}[\lvert v^{\top}\tilde{R}v\rvert>D]\leq
C\exp(Cn_{0})\cdot\exp\left({-\frac{n_{0}D^{2}}{\log
M}}\right)+\frac{C}{D^{2}}\cdot\frac{C}{M^{5}}.$ (S209)

Choosing $D=c_{3}\sqrt{\log M}$ for sufficiently large $c_{3}$ shows that

$\displaystyle\mathbb{P}[\lvert v^{\top}\tilde{R}v\rvert>D]\leq
C\exp(-cn_{0})+\frac{C}{n_{0}^{10}}\leq\frac{C}{n_{0}^{10}},$ (S210)

where $c>0$.

Recalling the original $\epsilon$-net covering, this implies that

$\displaystyle\lVert\tilde{R}\rVert_{\infty}\leq C\sqrt{\log n_{0}}$ (S211)

with probability at least $C/n^{10}$. We can now finish the argument by
relating the operator norm of $\bar{R}$ and $\tilde{R}$. Recalling the
rescaling between the modified and unmodified activations, we see

$\displaystyle\mathbb{P}[\lVert\bar{R}\rVert_{\infty}\geq C\sqrt{\log
n_{0}}\cdot
c_{0}\exp(c_{1}\lvert\bar{u}\rvert)]\leq\mathbb{P}[\lVert\tilde{R}\rVert_{\infty}\geq
C\sqrt{\log(n_{0})},\mathcal{G}]+\mathbb{P}[\mathcal{G}^{\mathsf{c}}]\leq
C/n_{0}^{10}.$ (S212)

####  I.1.2 Controlling the atypical behavior

Since the argument for the atypical event is the most straightforward, we
provide that first.

###### Lemma I.0

Suppose $\mathbb{P}[\mathcal{E}^{\mathsf{c}}]=Cn_{0}^{-c}$ for some
$n_{0}$-independent constants $C>0$ and $c>0$. Then,

$\left\lvert\mathbb{E}\left[{(1-\mathbf{1}_{\mathcal{E}})\Delta
y(\mathbf{x})}\right]\right\rvert\to 0$ (S213)

and

$\left\lvert\mathbb{E}\left[{(1-\mathbf{1}_{\mathcal{E}})\left({(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top})^{2}-(\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}Y^{\top})^{2}}\right)}\right]\right\rvert\to
0$ (S214)

as $n_{0}\to\infty$.

Proof  Applying the Cauchy-Schwarz inequality to the left-hand side of
Equation (S213), we may bound it by

$\mathbb{E}\left\lvert\left[{(1-\mathbf{1}_{\mathcal{E}})\Delta
y(\mathbf{x})}\right]\right\rvert\leq\left({\mathbb{E}\left[{\mathbf{1}_{\mathcal{E}^{\mathsf{c}}}}\right]}\right)^{1/2}\left({\mathbb{E}\left\lvert\Delta
y\right\rvert^{2}}\right)^{1/2}.$ (S215)

Since the first term $\mathbb{P}[\mathcal{E}_{c}]\to 0$ it also follows that
$(\mathbb{P}[\mathcal{E}_{c}])^{1/2}\to 0$. So it suffices to show
${\mathbb{E}[\lvert\Delta y\rvert^{2}]=O(1)}$.

We now apply the Cauchy-Schwarz inequality again to bound the second term by

$\displaystyle\mathbb{E}\left\lvert(\Delta
y(\mathbf{x}))^{2}\right\rvert\leq\sqrt{\mathbb{E}[\Delta^{4}]\mathbb{E}[y(\mathbf{x})^{4}]}.$
(S216)

A simple argument shows that $\mathbb{E}[y(\mathbf{x})^{4}]\leq
C(\mathbb{E}[(\beta^{\top}\mathbf{x})^{4}/n_{0}^{2}]+\mathbb{E}[\varepsilon^{4}])\leq
C\mathbb{E}[\lVert\Sigma^{1/2}\mathbf{z}\rVert^{4}/n_{0}^{2}]+C\sigma_{\varepsilon}^{4}\leq
C$ for $\mathbf{z}\sim\mathcal{N}(0,I_{d})$ since
$\mathbb{E}[\lVert\mathbf{z}\rVert_{2}^{4}]/n_{0}^{2}]\leq C$. Hence it
suffices to show $\mathbb{E}[\Delta^{4}]\leq C$. To this end note that

$\displaystyle\mathbb{E}[\Delta^{4}]\leq
C(\mathbb{E}[(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top})^{4}]+\mathbb{E}[(\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}Y^{\top})^{4}])$
(S217)

An identical argument can be applied to bound both terms. We outline the
argument for the first term. We see

$\displaystyle\lvert\mathbb{E}[(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top})^{4}]\rvert\leq\mathbb{E}\lvert\mathbb{E}_{\beta,\varepsilon}[(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top})^{4}|W,X,x]\rvert.$
(S218)

Now note that by the definition of $Y=\beta^{\top}X/\sqrt{n_{0}}+\varepsilon$
we have that,

$\displaystyle(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top})^{4}\leq
C\left((K_{\mathbf{x}}^{\top}K^{-1}\frac{X}{\sqrt{n_{0}}}\beta)^{4}+(K_{\mathbf{x}}^{\top}K^{-1}\varepsilon)^{4}\right)$
(S219)

Recalling that $\beta\sim\mathcal{N}(0,I_{n_{0}})$ and
$\varepsilon\sim\mathcal{N}(0,I_{m})$, we can compute the Gaussian moments
(marginally in $\beta,\varepsilon$) as,

$\displaystyle\mathbb{E}_{\beta,\varepsilon}[(K_{\mathbf{x}}^{\top}K^{-1}\frac{X}{\sqrt{n_{0}}}\beta)^{4}+(K_{\mathbf{x}}^{\top}K^{-1}\varepsilon)^{4}|W,X,x]$
(S220) $\displaystyle\qquad\leq
C\left((K_{\mathbf{x}}^{\top}K^{-1}\left(\frac{X^{\top}X}{n_{0}}\right)K^{-1}K_{\mathbf{x}}^{\top})^{2}+\sigma_{\varepsilon}^{2}(K_{\mathbf{x}}^{\top}K^{-2}K_{\mathbf{x}}^{\top})^{2}\right)$
(S221)

Continuing, using the definition of $K_{\mathbf{x}}$,

$\displaystyle\mathbb{E}[(K_{\mathbf{x}}^{\top}K^{-1}\left(\frac{X^{\top}X}{n_{0}}\right)K^{-1}K_{\mathbf{x}}^{\top})^{2}]$
$\displaystyle\leq\mathbb{E}\left[\left\|\frac{f^{\top}}{\sqrt{n_{0}}}\frac{FK^{-1}}{\sqrt{n_{0}}}(\frac{X^{\top}X}{n_{0}})\frac{K^{-1}F^{\top}}{\sqrt{n_{0}}}\frac{f}{\sqrt{n_{0}}}\right\|_{\infty}^{2}\right]$
(S222)
$\displaystyle\leq\mathbb{E}\left[{\left\lVert\frac{f}{\sqrt{n_{0}}}\right\rVert_{2}^{4}\cdot\left\lVert\frac{FK^{-1}}{\sqrt{n_{0}}}\right\rVert_{\infty}^{4}\cdot\left\lVert\frac{X^{\top}X}{n_{0}}\right\rVert_{\infty}^{2}}\right]$
(S223)

Noting that $K=FF^{\top}/n_{0}+\gamma I_{m}$ an application of the SVD to
$\frac{F}{\sqrt{n_{0}}}$ (denoting the corresponding singular values as
$s_{i}$) shows that
$\lVert\frac{FK^{-1}}{\sqrt{n_{0}}}\rVert_{\infty}=\max_{i}\frac{s_{i}}{s_{i}^{2}+\gamma}\leq\max_{s\geq
0}\frac{s}{s^{2}+\gamma}=1/({2\sqrt{\gamma}})$. Hence, we can continue
bounding Equation (S223) by

$\displaystyle
C\mathbb{E}[\left\lVert{f}/{\sqrt{n_{0}}}\right\rVert_{\infty}^{4}]\cdot\mathbb{E}[\|X^{\top}X/n_{0}\|_{\infty}^{2}],$
(S224)

where we exploited the independence of $f$ and $X$. Lem. I.0 ensures that
$\mathbb{E}[\lVert
X^{\top}X/n_{0}\rVert_{\infty}^{2}]=\mathbb{E}[\lVert\Upsilon\rVert_{\infty}^{2}]\leq
C$. The former term becomes
$\mathbb{E}[\|{f}/{\sqrt{n_{0}}}\|_{2}^{4}]=\mathbb{E}[\sigma(\sqrt{\bar{\mathrm{tr}}({\Sigma^{*}})}z)^{4}]\leq
C$ for $z\sim\mathcal{N}(0,1)$ by Sec. 2.2. Thus, the previous displays are
bounded by some constant $C$. An entirely analogous argument shows that
$\mathbb{E}[(K_{\mathbf{x}}^{\top}K^{-2}K_{\mathbf{x}}^{\top})^{2}]\leq C$.
Together these two results along with the previous computations establish
that,

$\displaystyle\lvert\mathbb{E}[(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top})^{4}]\rvert\leq
C$ (S225)

Note that our argument did not exploit any explicit properties of the centered
vs uncentered activation function $\sigma$ vs. $\bar{\sigma}$. Hence an
identical argument shows that,
$\lvert\mathbb{E}[(\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}Y^{\top})^{4}]\rvert\leq
C$. These two results imply ${\mathbb{E}[\Delta^{4}]}\leq C$ as desired.

Now we turn to the second term. Using an identical application of the Cauchy-
Schwarz inequality as used for the first term, we can bound Equation (S214) by

$\left({\mathbb{E}\left[{\mathbf{1}_{\mathcal{E}^{\mathsf{c}}}}\right]}\right)^{1/2}\left({\mathbb{E}\left\lvert(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top})^{2}-(\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}Y^{\top})^{2}\right\rvert^{2}}\right)^{1/2}.$
(S226)

Hence we can show Equation (S214) is $o(1)$ if
$\mathbb{E}\left\lvert(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top})^{2}-(\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}Y^{\top})^{2}\right\rvert^{2}=O(1)$.
This result follows immediately from the computations shown for the previous
term since we have already established that
$\mathbb{E}[(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top})^{4}]\leq C$ and
$\mathbb{E}[(\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}Y^{\top})^{4}]\leq C$.
These two previously shown results imply

$\displaystyle\mathbb{E}\left\lvert(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top})^{2}-(\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}Y^{\top})^{2}\right\rvert^{2}\leq
C\left(\mathbb{E}[(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top})^{4}]+\mathbb{E}[(\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}Y^{\top})^{4}]\right)\leq
C.$ (S227)

####  I.1.3 Schur complement formula for bounding $\Delta$

The centering procedure is a rank-1 change to $F$ and $f$, so applying the
Schur complement formula to understand its effect is natural. Write
${a\mathrel{\mathop{:}}=\mathbb{E}_{Z}\sigma\left({\sqrt{\bar{\mathrm{tr}}(\Sigma)}Z}\right)}$,
${a_{*}\mathrel{\mathop{:}}=\mathbb{E}_{Z}\sigma\left({\sqrt{\bar{\mathrm{tr}}(\Sigma^{*})}Z}\right)}$,
${\mathbf{v}\mathrel{\mathop{:}}=1/n_{1}\bar{F}^{\top}1_{n_{1}}}$,
${\mathbf{u}\mathrel{\mathop{:}}=a1_{m}}$,
${U\mathrel{\mathop{:}}=[\mathbf{u},\mathbf{v}]^{\top}}$, and
${C\mathrel{\mathop{:}}=\left(\begin{smallmatrix}1&1\\\
1&0\end{smallmatrix}\right)}$. Then,

$\displaystyle K^{-1}$
$\displaystyle=(\bar{K}+\mathbf{u}\mathbf{v}^{\top}+\mathbf{v}\mathbf{u}^{\top}+\mathbf{u}\mathbf{u}^{\top})^{-1}$
(S228) $\displaystyle=(\bar{K}+U^{\top}CU)^{-1}$ (S229)
$\displaystyle=\bar{K}^{-1}-\bar{K}^{-1}U^{\top}(C^{-1}+U\bar{K}^{-1}U^{\top})^{-1}U\bar{K}^{-1}\,,$
(S230)

and, for $\delta\mathrel{\mathop{:}}=1/n_{1}\bar{f}^{\top}1_{n_{1}}$ and
$P\mathrel{\mathop{:}}=(\delta+a_{*},a_{*})^{\top}$,

$\displaystyle K_{\mathbf{x}}$
$\displaystyle=\bar{K}_{\mathbf{x}}+\frac{1}{n_{1}}(F-\bar{F})^{\top}\bar{f}+\frac{1}{n_{1}}\bar{F}^{\top}(f-\bar{f})+\frac{1}{n_{1}}(F-\bar{F})^{\top}(f-\bar{f})$
(S231)
$\displaystyle=\bar{K}_{\mathbf{x}}+(\delta+a_{*})\mathbf{u}+a_{*}\mathbf{v}$
(S232) $\displaystyle=\bar{K}_{\mathbf{x}}+U^{\top}P\,.$ (S233)

Combining these expressions,

$\displaystyle K^{-1}K_{\mathbf{x}}$
$\displaystyle=K^{-1}(\bar{K}_{\mathbf{x}}+U^{\top}P)$ (S234)
$\displaystyle=\bar{K}^{-1}(\bar{K}_{\mathbf{x}}+U^{\top}P)-\bar{K}^{-1}U^{\top}(C^{-1}+U\bar{K}^{-1}U^{\top})^{-1}U\bar{K}^{-1}(\bar{K}_{\mathbf{x}}+U^{\top}P)$
(S235) $\displaystyle=\bar{K}^{-1}\bar{K}_{\mathbf{x}}+T_{1}+T_{2}\,,$ (S236)

with

$\displaystyle T_{1}$
$\displaystyle=\bar{K}^{-1}U^{\top}\left(I_{2}-(C^{-1}+U\bar{K}^{-1}U^{\top})^{-1}U\bar{K}^{-1}U^{\top}\right)P$
(S237)
$\displaystyle=\bar{K}^{-1}U^{\top}\left(I_{2}+CU\bar{K}^{-1}U^{\top}\right)^{-1}P$
(S238) $\displaystyle T_{2}$
$\displaystyle=-\bar{K}^{-1}U^{\top}(C^{-1}+U\bar{K}^{-1}U^{\top})^{-1}U\bar{K}^{-1}\bar{K}_{\mathbf{x}}$
(S239)
$\displaystyle=-\bar{K}^{-1}U^{\top}\left(I_{2}+CU\bar{K}^{-1}U^{\top}\right)^{-1}CU\bar{K}^{-1}\bar{K}_{\mathbf{x}}\,.$
(S240)

Furthermore, writing
${c_{0}\mathrel{\mathop{:}}=\tfrac{1}{m}\mathbf{u}^{\top}\bar{K}^{-1}\mathbf{u}}$,
${c_{1}\mathrel{\mathop{:}}=1+\mathbf{u}^{\top}\bar{K}^{-1}\mathbf{v}}$,
${c_{2}\mathrel{\mathop{:}}=1-\mathbf{v}^{\top}\bar{K}^{-1}\mathbf{v}}$,

$\displaystyle\left(I_{2}+CU\bar{K}^{-1}U^{\top}\right)^{-1}$
$\displaystyle=\frac{1}{c_{1}^{2}+mc_{0}c_{2}}\begin{pmatrix}c_{1}&c_{2}-c_{1}\\\
-mc_{0}&mc_{0}+c_{1}\end{pmatrix}\,,$ (S241)

so that,

$\displaystyle\left(I_{2}+CU\bar{K}^{-1}U^{\top}\right)^{-1}P$
$\displaystyle=\frac{1}{c_{1}^{2}+mc_{0}c_{2}}\begin{pmatrix}c_{1}\delta+c_{2}a_{*}\\\
-mc_{0}\delta+c_{1}a_{*}\end{pmatrix}\,,$ (S242)

and,

$\displaystyle\left(I_{2}+CU\bar{K}^{-1}U^{\top}\right)^{-1}C$
$\displaystyle=\frac{1}{c_{1}^{2}+mc_{0}c_{2}}\begin{pmatrix}c_{2}&c_{1}\\\
c_{1}&-mc_{0}\end{pmatrix}\,.$ (S243)

####  I.1.4 Concentration with high-probability

###### Definition I.0

Recall the definitions from Sec. I.1.3. To characterize the typical behavior
of the random variables, we define the event

$\displaystyle\mathcal{E}\mathrel{\mathop{:}}=\mathcal{E}_{\text{data}}\cap\mathcal{E}_{W}\cap\mathcal{E}_{F}\cap\left\\{{|\delta|\leq
Cn_{0}^{c-1/2}}\right\\}\cap\left\\{{cn_{0}^{-c}\leq c_{0}\leq
Cn_{0}^{c}}\right\\}\cap\left\\{{c_{1}\leq Cn_{0}^{c+1/2}}\right\\}$ (S244)
$\displaystyle\quad\cap\left\\{{cn_{0}^{-c}\leq c_{2}\leq
Cn_{0}^{c}}\right\\}\cap\left\\{{Y\bar{K}^{-1}\mathbf{v}^{\top}\leq
Cn_{0}^{c}}\right\\}\cap\left\\{{Y\bar{K}^{-1}\mathbf{u}^{\top}\leq
Cn_{0}^{c+1/2}}\right\\}.$ (S245)
$\displaystyle\cap\left\\{{\bar{K}_{\mathbf{x}}\bar{K}^{-1}\mathbf{u}^{\top}\leq
Cn_{0}^{c}}\right\\}\cap\left\\{{\bar{K}_{\mathbf{x}}\bar{K}^{-1}\mathbf{v}^{\top}\leq
Cn_{0}^{c-1/2}}\right\\}\cap\left\\{{y(\mathbf{x})\leq Cn_{0}^{c}}\right\\}$
(S246)

###### Lemma I.0

The event $\mathcal{E}$ is high-probability, that is,
$\mathbb{P}[\mathcal{E}^{\mathsf{c}}]\leq Cn_{0}^{-c}$ for some constant $C>0$
and any constant $0<c<10$.

Proof  First, we already know $\mathcal{E}_{\text{data}}$, $\mathcal{E}_{W}$,
and $\mathcal{E}_{F}$ all occur with high-probability (see Lems. I.0, I.0 and
I.0). When bounding the other events in Equations (S244)-(S246), it will be
useful to introduce conditioning on one of $\mathcal{E}_{\text{data}}$,
$\mathcal{E}_{W}$, or $\mathcal{E}_{F}$. This will suffice since for any event
$\mathcal{A}$, we have

$\mathbb{P}[\mathcal{A}^{\mathsf{c}}]\leq\mathbb{P}[\mathcal{A}^{\mathsf{c}}|\mathcal{E}_{\text{data}}]+\mathbb{P}[\mathcal{E}_{\text{data}}^{\mathsf{c}}]\leq\mathbb{P}[\mathcal{A}^{\mathsf{c}}|\mathcal{E}_{\text{data}}]+Cn_{0}^{-c}$
(S247)

for example. We will explicitly denote this conditioning, but when taking
expectation over some random variables (e.g. $x$ or $W$) we will not
explicitly denote that we are conditioning on the remaining independent random
variables (e.g. $X$). Once we have bounded the probability of the complement
of each event in Equations (S244)-(S246), we can apply the union bound to
complete the proof.

For future reference recall
${\mathbf{v}\mathrel{\mathop{:}}=1/n_{1}\bar{F}^{\top}1_{n_{1}}}$,
${\mathbf{u}\mathrel{\mathop{:}}=a1_{m}}$. We deal with each event in
Equations (S244)-(S246) in turn.

##### Controlling $\delta$:

Here, we introduce conditioning on the event $\mathcal{E}_{\text{data}}$. To
control $\delta$, we want to calculate its mean and variance to apply
Chebyshev’s inequality. Note that $\mathcal{E}_{\text{data}}$ is independent
of $W$, so the expectations over $W$ below are unchanged by this conditioning.
Recall

$\delta=\frac{1}{n_{1}}\sum_{k=1}^{n_{1}}{\sigma}\left({\sum_{j=1}^{n_{0}}W_{kj}x_{j}/\sqrt{n_{0}}}\right)-a_{*}$
(S248)

and note that conditional on $x$ the sum $\sum_{j}W_{kj}x_{j}/\sqrt{n_{0}}$ is
distributed as $\mathcal{N}(0,\upsilon_{*})$ for all $k$.

Expanding in $\upsilon_{*}-\bar{\mathrm{tr}}(\Sigma^{*})$, we see

$\displaystyle\mathbb{E}_{W}[\delta|\mathcal{E}_{\text{data}}]$
$\displaystyle=\mathbb{E}[\mathbb{E}_{Z\sim\mathcal{N}(0,1)}[{\sigma}(\sqrt{\bar{\mathrm{tr}}(\Sigma^{*})+(\upsilon_{*}-\bar{\mathrm{tr}}(\Sigma^{*}))}Z)]-a_{*}|\mathcal{E}_{\text{data}}]$
(S249)
$\displaystyle=\mathbb{E}[\mathbb{E}_{Z}[{\sigma}(\sqrt{\bar{\mathrm{tr}}(\Sigma^{*})}Z)]-a_{*}|\mathcal{E}_{\text{data}}]$
(S250)
$\displaystyle\qquad+\mathbb{E}\left[{\frac{1}{2}\frac{(\upsilon_{*}-\bar{\mathrm{tr}}(\Sigma^{*}))}{\bar{\mathrm{tr}}(\Sigma^{*})}\mathbb{E}_{Z}[\sqrt{\bar{\mathrm{tr}}(\Sigma^{*})}Z{\sigma}(\sqrt{\bar{\mathrm{tr}}(\Sigma^{*})}Z)]+\Delta\Big{|}\mathcal{E}_{\text{data}}}\right]$
(S251)
$\displaystyle=\mathbb{E}\left[{\frac{1}{2}\frac{(\upsilon_{*}-\bar{\mathrm{tr}}(\Sigma^{*}))}{\bar{\mathrm{tr}}(\Sigma^{*})}\mathbb{E}_{Z}[\sqrt{\bar{\mathrm{tr}}(\Sigma^{*})}Z{\sigma}(\sqrt{\bar{\mathrm{tr}}(\Sigma^{*})}Z)]+\Delta|\mathcal{E}_{\text{data}}}\right],$
(S252)

where $Z$ is a standard Gaussian. The first term of Equation (S252) is less
than $Cn_{0}^{c-1/2}$ and the remainder term $\Delta$ is less than
$Cn_{0}^{c-1}$ on $\mathcal{E}_{\text{data}}$. More detail on this argument
can be found in Adlam et al. (2019). Taking expectation over the remaining
randomness, we see $\mathbb{E}[\delta|\mathcal{E}_{\text{data}}]\leq
Cn_{0}^{c-1/2}$.

Similarly,

$\mathbb{V}_{W}[\delta|\mathcal{E}_{\text{data}}]=\frac{1}{n_{1}^{2}}\sum_{k=1}^{n_{1}}\mathbb{V}_{W}\left[{{\sigma}\left({\sum_{j=1}^{n_{0}}W_{kj}x_{j}/\sqrt{n_{0}}}\right)\Big{|}\mathcal{E}_{\text{data}}}\right]=\frac{1}{n_{1}}\mathbb{V}_{Z\sim\mathcal{N}(0,\upsilon_{*})}{\sigma}(Z),$
(S253)

where we can again Taylor expand in
$\upsilon_{*}-\bar{\mathrm{tr}}(\Sigma^{*})$ to get the bound
$\mathbb{V}_{W}[\delta|\mathcal{E}_{\text{data}}]\leq C/n_{0}$ on
$\mathcal{E}_{\text{data}}$. Using the law of total variance, we can then
bound $\mathbb{V}[\delta|\mathcal{E}_{\text{data}}]\leq C/n_{0}$.

Finally applying Chebyshev’s inequality implies for all $c>0$ that

$\mathbb{P}[\lvert\delta-\mathbb{E}[\delta|\mathcal{E}_{\text{data}}]\rvert>Cn_{0}^{c-1/2}|\mathcal{E}_{\text{data}}]\leq
Cn_{0}^{-c},$ (S254)

and so

$\mathbb{P}[\lvert\delta\rvert>Cn_{0}^{c-1/2}|\mathcal{E}_{\text{data}}]\leq
Cn_{0}^{-c}.$ (S255)

##### Controlling $c_{0}$, $c_{1}$, and $c_{2}$:

Here, we introduce conditioning on the event $\mathcal{E}_{F}$. The results
for $c_{0}$, $c_{1}$, and $c_{2}$ can be found using a similar argument to
terms in Mei and Montanari (2021, Lemma 9.6/Step 3). The identifications
$c_{0}\leftrightarrow K_{11}$, $c_{1}\leftrightarrow K_{12}$,
$c_{2}\leftrightarrow 1-K_{22}$ hold.

For $c_{0}=\frac{1}{m}\mathbf{u}^{\top}\bar{K}^{-1}\mathbf{u}$ we have that

$c_{0}\leq\frac{a^{2}\lVert
1_{m}\rVert^{2}}{m}\lVert\bar{K}^{-1}\rVert_{\infty}\leq\frac{C}{\gamma}$
(S256)

and

$c_{0}\geq\frac{a^{2}\lVert
1_{m}\rVert^{2}}{m}\lambda_{\min}(\bar{K}^{-1})\geq\frac{C}{(\gamma+\lVert\bar{F}\rVert_{\infty}^{2}/n_{0})}\geq
cn_{0}^{-c}$ (S257)

using the operator norm bound on the event $\mathcal{E}_{F}$.

For $c_{1}=1+\mathbf{u}^{\top}\bar{K}^{-1}\mathbf{v}$ we have

$\lvert c_{1}\rvert=\left\lvert
1+a1_{m}^{\top}\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}\frac{1_{n_{1}}}{\sqrt{n_{1}}}\right\rvert\leq
1+a\lVert 1_{m}\rVert_{2}\frac{\lVert
1_{n_{1}}\rVert_{2}}{\sqrt{n_{1}}}\lVert\bar{K}^{-1}\frac{\bar{F}}{\sqrt{n_{1}}}\rVert\leq
Cn_{0}^{1/2},$ (S258)

where the bound $\lVert\bar{K}^{-1}\frac{\bar{F}}{\sqrt{n_{1}}}\rVert\leq C$
follows by considering the SVD of $\bar{F}$.

For $c_{2}=1-\mathbf{v}^{\top}\bar{K}^{-1}\mathbf{v}$ we obviously have that
$c_{2}\leq 1$ since $\mathbf{v}^{\top}\bar{K}^{-1}\mathbf{v}\geq 0$.
Additionally, we have that

$\displaystyle c_{2}$
$\displaystyle=1-\frac{1}{n_{1}}1_{n_{1}}^{\top}\frac{\bar{F}}{\sqrt{n_{1}}}\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}1_{n_{1}}$
(S259)
$\displaystyle=\frac{1_{n_{1}}^{\top}}{\sqrt{n_{1}}}\left({I_{n_{1}}-\frac{\bar{F}}{\sqrt{n_{1}}}\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}}\right)\frac{1_{n_{1}}}{\sqrt{n_{1}}}$
(S260) $\displaystyle\geq
1-\lVert\frac{\bar{F}}{\sqrt{n_{1}}}\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}\rVert_{\infty}$
(S261)
$\displaystyle=1-\frac{\lVert\bar{F}/\sqrt{n_{1}}\rVert_{\infty}^{2}}{\gamma+\lVert\bar{F}/\sqrt{n_{1}}\rVert_{\infty}^{2}}\geq
1-\frac{cn_{0}^{c}}{\gamma+cn_{0}^{c}}$ (S262)
$\displaystyle\geq\frac{C\gamma}{cn_{0}^{c}}$ (S263) $\displaystyle\geq
cn_{0}^{-c}$ (S264)

as desired once again using the conditioning on the event $\mathcal{E}_{F}$.

Finally,

$\mathbf{v}^{\top}\bar{K}^{-1}\mathbf{v}=\frac{1}{n_{1}}1_{n_{1}}^{\top}\frac{\bar{F}}{\sqrt{n_{1}}}\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}1_{n_{1}}\leq\frac{\lVert
1_{n_{1}}\rVert_{2}^{2}}{n_{1}}\left\lVert\frac{\bar{F}}{\sqrt{n_{1}}}\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}\right\rVert_{\infty}\leq
C$ (S265)

since
$\frac{\bar{F}}{\sqrt{n_{1}}}\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}\leq
C$ follows once again by considering the SVD of $\bar{F}$.

##### Controlling $Y\bar{K}^{-1}\mathbf{v}^{\top}$:

Here, we introduce conditioning on the event $\mathcal{E}_{\text{data}}$ and
note its independence from $\beta$ and $\varepsilon$. Recalling that
$Y=\beta^{\top}X/\sqrt{n_{0}}+\varepsilon$, we note that $Y$ has expectation
zero since
$\mathbb{E}_{\beta,\varepsilon}[Y\bar{K}^{-1}\mathbf{v}^{\top}|\mathcal{E}_{\text{data}}]=0$.
Similarly the conditional variance can be bounded as

$\displaystyle\mathbb{V}_{\beta,\varepsilon}[Y\bar{K}^{-1}\mathbf{v}^{\top}|\mathcal{E}_{\text{data}}]$
$\displaystyle=\mathbf{v}^{\top}\bar{K}^{-1}(\Upsilon+\sigma_{\varepsilon}^{2}I_{m})\bar{K}^{-1}\mathbf{v}$
(S266)
$\displaystyle\leq\left\lVert\Upsilon+\sigma_{\varepsilon}^{2}I_{m}\right\rVert_{\infty}\left\lVert\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}\right\rVert_{\infty}^{2}\left\lVert\frac{1_{n_{1}}}{\sqrt{n_{1}}}\right\rVert_{2}^{2}$
(S267) $\displaystyle\leq C.$ (S268)

The bound $\|\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}\|_{\infty}\leq C$
follows by considering the SVD of $F$ as in the proof of Lem. I.0, while the
bound $\|\Upsilon\|\leq C$ holds on $\mathcal{E}_{\text{data}}$. The law of
total variance then shows
${\mathbb{V}[Y\bar{K}^{-1}\mathbf{v}^{\top}|\mathcal{E}_{\text{data}}]\leq C}$
also. Hence an application of Chebyshev’s inequality shows that

$\displaystyle\mathbb{P}[\lvert Y\bar{K}^{-1}\mathbf{v}^{\top}\rvert\geq
Cn_{0}^{c}|\mathcal{E}_{\text{data}}]\leq Cn_{0}^{-c}$ (S269)

for some $C>0$ and any $c>0$.

##### Controlling $Y\bar{K}^{-1}\mathbf{u}^{\top}$:

Again, we introduce conditioning on the event $\mathcal{E}_{\text{data}}$ and
note its independence from $\beta$ and $\varepsilon$. A nearly identical
argument to the previous one suffices, so we only outline the argument. Again
$\mathbb{E}[Y\bar{K}^{-1}\mathbf{v}^{\top}|\mathcal{E}_{\text{data}}]=0$ and
the conditional variance can be bounded as

$\displaystyle\mathbb{V}_{\beta,\varepsilon}[Y\bar{K}^{-1}\mathbf{u}^{\top}|\mathcal{E}_{\text{data}}]$
$\displaystyle=\mathbf{u}^{\top}\bar{K}^{-1}(\Upsilon+\sigma_{\varepsilon}^{2}I_{m})\bar{K}^{-1}\mathbf{u}$
(S270)
$\displaystyle\leq\lVert\bar{K}^{-1}\rVert_{\infty}^{2}\|\Upsilon+\sigma_{\varepsilon}^{2}I_{m}\|_{\infty}\|\mathbf{u}\|_{2}^{2}$
(S271) $\displaystyle\leq\frac{1}{\gamma^{2}}\cdot
a^{2}m\cdot\|\Upsilon+\sigma_{\varepsilon}^{2}I_{m}\|_{\infty}$ (S272)
$\displaystyle\leq Cn_{0}.$ (S273)

Finally, Chebyshev’s inequality suffices to show

$\displaystyle\mathbb{P}[\lvert Y\bar{K}^{-1}\mathbf{u}^{\top}\rvert\geq
Cn_{0}^{c+1/2}|\mathcal{E}_{data}]\leq Cn_{0}^{-c}$ (S274)

for some $C>0$ and any $c>0$.

##### Controlling $\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}\mathbf{u}^{\top}$:

Here, we introduce conditioning on the event $\mathcal{E}_{W}$ and note its
independence from $\mathbf{x}$. The overall argument then uses Chebyshev’s
inequality exploiting the randomness in $\mathbf{x}$. First,

$\displaystyle\left\lvert\mathbb{E}_{\mathbf{x}}[\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}\mathbf{u}^{\top}|\mathcal{E}_{W}]\right\rvert$
$\displaystyle=\left\lvert\mathbb{E}_{\mathbf{x}}[\bar{f}^{\top}|\mathcal{E}_{W}]\frac{\bar{F}}{\sqrt{n_{1}}}\bar{K}^{-1}\mathbf{u}^{\top}/\sqrt{n_{1}}\right\rvert$
(S275)
$\displaystyle\leq\left\lVert\mathbb{E}_{\mathbf{x}}[\bar{f}|\mathcal{E}_{W}]\right\rVert_{2}\left\lVert\frac{\bar{F}}{\sqrt{n_{1}}}\bar{K}^{-1}\right\rVert_{\infty}\left\lVert\mathbf{u}/\sqrt{n_{1}}\right\rVert_{2}$
(S276) $\displaystyle\leq Cn_{0}^{c}.$ (S277)

The bound $\lVert\frac{\bar{F}}{\sqrt{n_{1}}}\bar{K}^{-1}\rVert_{\infty}\leq
C$ follows by an SVD as before, while the first bound
${\lVert\mathbb{E}_{\mathbf{x}}[\bar{f}^{\top}|\mathcal{E}_{W}]\rVert_{2}\leq
Cn_{0}^{c}}$ follows from the fact
$\mathbb{E}_{\mathbf{x}}\bar{\sigma}(W\mathbf{x})_{i}\leq Cn_{0}^{c-1/2}$ on
$\mathcal{E}_{W}$ (see Lem. I.0). Next we turn to the conditional variance. We
see

$\displaystyle\mathbb{V}_{\mathbf{x}}[\bar{K}_{\mathbf{x}}\bar{K}^{-1}\mathbf{u}|\mathcal{E}_{W}]$
$\displaystyle=\mathbf{u}^{\top}\bar{K}^{-1}\mathbb{E}_{\mathbf{x}}[\bar{K}_{\mathbf{x}}\bar{K}_{\mathbf{x}}^{\top}|\mathcal{E}_{W}]\bar{K}^{-1}\mathbf{u}$
(S278)
$\displaystyle=\frac{1}{n_{1}^{2}}\mathbf{u}^{\top}\bar{K}^{-1}\bar{F}^{\top}\big{(}\frac{\rho_{*}}{n_{0}}W\Sigma^{*}W^{\top}+(\eta_{*}-\zeta_{*})I_{n_{1}}+\Delta\big{)}\bar{F}\bar{K}^{-1}\mathbf{u}$
(S279)
$\displaystyle\leq\left\lVert\frac{\mathbf{u}}{\sqrt{n_{1}}}\right\rVert_{2}^{2}\left\lVert\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}\right\rVert_{2}^{2}\cdot\left({\left\lVert\rho_{*}\tilde{\Upsilon}\right\rVert_{\infty}\\!+\left\lVert(\eta_{*}-\zeta_{*})I_{n_{1}}\right\rVert_{\infty}\\!+\left\lVert\Delta\right\rVert_{\infty}}\right)$
(S280) $\displaystyle\leq C\cdot(C+C+n_{0}^{c})$ (S281) $\displaystyle\leq
Cn_{0}^{c}.$ (S282)

Once again the bound
$\lVert\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}\rVert_{2}^{2}$ follows
by considering the SVD, while the nontrivial operator norm bounds follow from
Lem. I.0. Using the law of total variance, we see

$\mathbb{V}[\bar{K}_{\mathbf{x}}\bar{K}^{-1}\mathbf{u}|\mathcal{E}_{W}]=\mathbb{E}[\mathbb{V}_{\mathbf{x}}[\bar{K}_{\mathbf{x}}\bar{K}^{-1}\mathbf{u}|\mathcal{E}_{W}]|\mathcal{E}_{W}]+\mathbb{V}[\mathbb{E}_{\mathbf{x}}[\bar{K}_{\mathbf{x}}\bar{K}^{-1}\mathbf{u}|\mathcal{E}_{W}]|\mathcal{E}_{W}]\leq
Cn_{0}^{c}$ (S283)

by Equations (S277) and (S282). Applying Chebyshev’s inequality yields

$\mathbb{P}[\lvert\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}\mathbf{u}-\mathbb{E}_{\mathbf{x}}[\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}\mathbf{u}]\rvert\geq
Cn_{0}^{c}|\mathcal{E}_{W}]\leq Cn_{0}^{-c},$ (S284)

and so

$\mathbb{P}[\lvert\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}\mathbf{u}\rvert\geq
Cn_{0}^{c}|\mathcal{E}_{W}]\leq Cn_{0}^{-c}$ (S285)

by the triangle inequality since
$\lvert\mathbb{E}_{\mathbf{x}}[\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}\mathbf{u}]\rvert\leq
Cn_{0}^{c}$.

##### Controlling $\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}\mathbf{v}^{\top}$:

Again, we introduce conditioning on the event $\mathcal{E}_{W}$ and note its
independence from $\mathbf{x}$. A nearly identical argument to the previous
one shows the result for this term. Hence we only outline the argument. First,
we see

$\displaystyle\lvert\mathbb{E}_{\mathbf{x}}[\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}\mathbf{v}^{\top}|\mathcal{E}_{W}]\rvert$
$\displaystyle=\left\lvert\mathbb{E}_{\mathbf{x}}[\bar{f}^{\top}|\mathcal{E}_{W}]\frac{\bar{F}}{\sqrt{n_{1}}}\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}\frac{1_{n_{1}}}{n_{1}}\right\rvert$
(S286)
$\displaystyle\leq\left\lVert\mathbb{E}_{\mathbf{x}}[\bar{f}^{\top}|\mathcal{E}_{W}]\right\rVert_{2}\left\lVert\frac{\bar{F}}{\sqrt{n_{1}}}\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}\right\rVert_{\infty}\left\lVert\frac{1_{n_{1}}}{n_{1}}\right\rVert_{2}$
(S287) $\displaystyle\leq Cn_{0}^{c-1/2}.$ (S288)

The bound
$\lVert\frac{\bar{F}}{\sqrt{n_{1}}}\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}\rVert_{\infty}\leq
1$ follows by an SVD as before, while the first bound again follows from Lem.
I.0. Next we compute the conditional variance and bound as in Equation (S282)
to find

$\displaystyle\mathbb{V}_{\mathbf{x}}[\bar{K}_{\mathbf{x}}\bar{K}^{-1}\mathbf{v}^{\top}|\mathcal{E}_{W}]$
$\displaystyle=\mathbf{v}^{\top}\bar{K}^{-1}\mathbb{E}_{\mathbf{x}}[\bar{K}_{\mathbf{x}}\bar{K}_{\mathbf{x}}^{\top}|\mathcal{E}_{W}]\bar{K}^{-1}\mathbf{v}$
(S289)
$\displaystyle=\frac{1}{n_{1}^{4}}1_{n_{1}}^{\top}\bar{F}\bar{K}^{-1}\bar{F}^{\top}\left({\frac{\rho_{*}}{n_{0}}W\Sigma^{*}W^{\top}+(\eta_{*}-\zeta_{*})I_{n_{1}}+\Delta}\right)\bar{F}\bar{K}^{-1}\bar{F}^{\top}1_{n_{1}}$
(S290)
$\displaystyle\leq\frac{1}{n_{1}}\left\lVert\frac{1_{n_{1}}}{\sqrt{n_{1}}}\right\rVert_{2}^{2}\\!\cdot\left\lVert\frac{\bar{F}}{\sqrt{n_{1}}}\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}\right\rVert_{2}^{2}\\!\cdot\left({\left\lVert{\rho_{*}}\tilde{\Upsilon}\right\rVert_{\infty}\\!+\\!\left\lVert(\eta_{*}\\!-\\!\zeta_{*})I_{n_{1}}\right\rVert_{\infty}\\!+\\!\left\lVert\Delta\right\rVert_{\infty}}\right)$
(S291) $\displaystyle\leq\frac{C\cdot(C+C+n_{0}^{c})}{n_{1}}$ (S292)
$\displaystyle\leq Cn_{0}^{c-1}.$ (S293)

Once again the bound
$\lVert\frac{\bar{F}}{\sqrt{n_{1}}}\bar{K}^{-1}\frac{\bar{F}^{\top}}{\sqrt{n_{1}}}\rVert_{2}\leq
1$ follows by considering the SVD while the nontrivial operator norm bounds
follow on $\mathcal{E}_{W}$. Now applying Chebyshev’s inequality as before
shows

$\displaystyle\mathbb{P}[\lvert\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}\mathbf{v}^{\top}\rvert\geq
Cn_{0}^{c-1/2}|\mathcal{E}_{W}]\leq Cn_{0}^{-c}.$ (S294)

##### Controlling $y(\mathbf{x})$:

Finally we show $y(\mathbf{x})\leq Cn_{0}^{c}$ with probability at least
$Cn_{0}^{-c}$. Note that $\mathbb{E}[y(\mathbf{x})]=0$ and a short computation
shows that
$\mathbb{V}[y(\mathbf{x})]=\bar{\mathrm{tr}}(\Sigma^{*})+\sigma_{\varepsilon}^{2}\leq
C$. So a direct application of Chebyshev’s inequality shows that

$\displaystyle\mathbb{P}[\lvert y(\mathbf{x})\rvert\geq Cn_{0}^{c}]\leq
Cn_{0}^{-c}.$ (S295)

####  I.1.5 Completing the argument for the typical behavior

###### Lemma I.0

For some $n_{0}$-independent constants $c>0$ and $C>0$, we have

$\left\lvert\mathbb{E}\left[{\mathbf{1}_{\mathcal{E}}\Delta
y(\mathbf{x})}\right]\right\rvert\leq Cn_{0}^{-c}$ (S296)

and

$\left\lvert\mathbb{E}\left[{\mathbf{1}_{\mathcal{E}}\Delta(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top}+\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}Y^{\top})}\right]\right\rvert\leq
Cn_{0}^{-c}$ (S297)

Proof  Obviously, $\left\lvert\mathbb{E}\left[{\mathbf{1}_{\mathcal{E}}\Delta
y(\mathbf{x})}\right]\right\rvert$ is bounded by
$\left\lvert\mathbb{E}\left[{\Delta
y(\mathbf{x})}\right]|\mathcal{E}\right\rvert$, so we can remove the indicator
$\mathbf{1}_{\mathcal{E}}$ from the expectations in Equations (S296) and
(S297) by conditioning on $\mathcal{E}$.

We want to show that conditional on $\mathcal{E}$ and for some
$n_{0}$-independent constants $c>0$ and $C>0$, the bounds
$\left\lvert\Delta\right\rvert\leq Cn_{0}^{-2c}$, $\left\lvert
y(\mathbf{x})\right\rvert\leq Cn_{0}^{c}$, and
$\left\lvert(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top}+\bar{K}_{\mathbf{x}}^{\top}\bar{K}^{-1}Y^{\top})\right\rvert\leq
Cn_{0}^{c}$ hold.

Recall,
${a\mathrel{\mathop{:}}=\mathbb{E}_{Z}\sigma\left({\sqrt{\bar{\mathrm{tr}}(\Sigma)}Z}\right)}$,
${a_{*}\mathrel{\mathop{:}}=\mathbb{E}_{Z}\sigma\left({\sqrt{\bar{\mathrm{tr}}(\Sigma^{*})}Z}\right)}$,
${\mathbf{v}\mathrel{\mathop{:}}=1/n_{1}\bar{F}^{\top}1_{n_{1}}}$,
${\mathbf{u}\mathrel{\mathop{:}}=a1_{m}}$,
${U\mathrel{\mathop{:}}=[\mathbf{u},\mathbf{v}]^{\top}}$, and
${C\mathrel{\mathop{:}}=\left(\begin{smallmatrix}1&1\\\
1&0\end{smallmatrix}\right)}$.

First, consider the term $\Delta=YT_{1}+YT_{2}$. To show that
$\left\lvert\Delta\right\rvert\leq Cn_{0}^{-2c}$ on $\mathcal{E}$, we write
out $YT_{1}$ and $YT_{2}$ more explicitly:

$\displaystyle YT_{1}$
$\displaystyle=Y\bar{K}^{-1}U^{\top}\left(I_{2}+CU\bar{K}^{-1}U^{\top}\right)^{-1}P$
(S298)
$\displaystyle=\frac{1}{c_{1}^{2}+mc_{0}c_{2}}\begin{pmatrix}Y\bar{K}^{-1}\mathbf{u}\\\
Y\bar{K}^{-1}\mathbf{v}\end{pmatrix}^{\top}\begin{pmatrix}c_{1}\delta+c_{2}a_{*}\\\
-mc_{0}\delta+c_{1}a_{*}\end{pmatrix}$ (S299)
$\displaystyle=\frac{1}{c_{1}^{2}+mc_{0}c_{2}}\left(Y\bar{K}^{-1}\mathbf{u}\cdot(c_{1}\delta+c_{2}a_{*})+Y\bar{K}^{-1}\mathbf{v}\cdot(-mc_{0}\delta+c_{1}a_{*})\right)$
(S300)

and

$\displaystyle YT_{2}$
$\displaystyle=-Y\bar{K}^{-1}U^{\top}\left(I_{2}+CU\bar{K}^{-1}U^{\top}\right)^{-1}CU\bar{K}^{-1}\bar{K}_{\mathbf{x}}$
(S302)
$\displaystyle=\frac{1}{c_{1}^{2}+mc_{0}c_{2}}\begin{pmatrix}Y\bar{K}^{-1}\mathbf{u}\\\
Y\bar{K}^{-1}\mathbf{v}\end{pmatrix}^{\top}\begin{pmatrix}c_{2}&c_{1}\\\
c_{1}&-mc_{0}\end{pmatrix}\begin{pmatrix}\mathbf{u}^{\top}\bar{K}^{-1}\bar{K}_{\mathbf{x}}\\\
\mathbf{v}^{\top}\bar{K}^{-1}\bar{K}_{\mathbf{x}}\end{pmatrix}$ (S303)
$\displaystyle\leq\frac{1}{c_{1}^{2}+mc_{0}c_{2}}\Big{(}Y\bar{K}^{-1}\mathbf{u}\cdot(c_{2}\mathbf{u}^{\top}\bar{K}^{-1}\bar{K}_{\mathbf{x}}+c_{1}\mathbf{v}^{\top}\bar{K}^{-1}\bar{K}_{\mathbf{x}})$
(S304)
$\displaystyle\qquad+Y\bar{K}^{-1}\mathbf{v}\cdot(c_{1}\mathbf{u}^{\top}\bar{K}^{-1}\bar{K}_{\mathbf{x}}-mc_{0}\mathbf{v}^{\top}\bar{K}^{-1}\bar{K}_{\mathbf{x}})\Big{)}.$
(S305)

Now using Equation (S298) and the definition of $\mathcal{E}$, we see that on
$\mathcal{E}$

$\displaystyle\left\lvert YT_{1}\right\rvert$
$\displaystyle\leq\left\lvert\frac{1}{c_{1}^{2}+mc_{0}c_{2}}\right\rvert\left({\left\lvert
Y\bar{K}^{-1}\mathbf{u}\right\rvert(\left\lvert
c_{1}\right\rvert\left\lvert\delta\right\rvert+\left\lvert
c_{2}\right\rvert\left\lvert a_{*}\right\rvert)+\left\lvert
Y\bar{K}^{-1}\mathbf{v}\right\rvert(m\left\lvert
c_{0}\right\rvert\left\lvert\delta\right\rvert+\left\lvert
c_{1}\right\rvert\left\lvert a_{*}\right\rvert)}\right)$ (S306)
$\displaystyle\leq
Cn_{0}^{2c-1}\left({Cn_{0}^{c+1/2}(Cn_{0}^{c}Cn_{0}^{c-1/2}+Cn_{0}^{c})+Cn_{0}^{c}(Cn_{0}Cn_{0}^{c}Cn_{0}^{c-1/2}+Cn_{0}^{c})}\right)$
(S307) $\displaystyle\leq Cn_{0}^{5c-1/2}.$ (S308)

Since the $c$s used in the bounds above can be arbitrarily small, we may
replace write $\left\lvert YT_{1}\right\rvert\leq Cn_{0}^{c-1/2}$ for
arbitrarily small $c>0$. Finally, a similar argument for Equation (S302)
yields the bound

$\displaystyle\left\lvert YT_{2}\right\rvert$
$\displaystyle\leq\left\lvert\frac{1}{c_{1}^{2}+mc_{0}c_{2}}\right\rvert\Big{(}\lvert
Y\bar{K}^{-1}\mathbf{u}\rvert\cdot(\lvert
c_{2}\rvert\lvert\mathbf{u}^{\top}\bar{K}^{-1}\bar{K}_{\mathbf{x}}\rvert+\lvert
c_{1}\rvert\lvert\mathbf{v}^{\top}\bar{K}^{-1}\bar{K}_{\mathbf{x}}\rvert)$
(S309) $\displaystyle\qquad+\lvert Y\bar{K}^{-1}\mathbf{v}\rvert\cdot(\lvert
c_{1}\rvert\lvert\mathbf{u}^{\top}\bar{K}^{-1}\bar{K}_{\mathbf{x}}\rvert-m\lvert
c_{0}\rvert\lvert\mathbf{v}^{\top}\bar{K}^{-1}\bar{K}_{\mathbf{x}}\rvert)\Big{)}$
(S310) $\displaystyle\leq
Cn_{0}^{2c-1}\cdot\Big{(}Cn_{0}^{c+1/2}\cdot(Cn_{0}^{c}Cn_{0}^{c}+Cn_{0}^{c+1/2}Cn_{0}^{c-1/2})$
(S311)
$\displaystyle\qquad+Cn_{0}^{c}\cdot(Cn_{0}^{c+1/2}Cn_{0}^{c}+Cn_{0}Cn_{0}^{c}Cn_{0}^{c-1/2})\Big{)}$
(S312) $\displaystyle\leq Cn_{0}^{2c-1}\cdot(Cn_{0}^{3c+1/2}+Cn_{0}^{3c+1/2})$
(S313) $\displaystyle\leq Cn_{0}^{5c-1/2}.$ (S314)

and as before since $c$ can be chosen arbitrarily small we can write
$\left\lvert YT_{2}\right\rvert\leq Cn_{0}^{c-1/2}$.

###  I.2 Gaussian equivalents

Before describing how (in the high-dimensional limit) Gaussian equivalents can
be used to compute the test error we make the following observation: the
limiting test error is invariant to the centering of the activation function
across the training and test distributions. That is, by Sec. I.1 the limiting
test error is unchanged by replacing,

$\bar{F}_{ij}\mathrel{\mathop{:}}=F_{ij}-\mathbb{E}_{z\sim\mathcal{N}(0,s)}\sigma\left({z}\right)\quad\text{and}\quad\bar{f}_{i}\mathrel{\mathop{:}}=f_{i}-\mathbb{E}_{z\sim\mathcal{N}(0,s_{*})}\sigma\left({z}\right)$
(S315)

for all $i$ and $j$.

With this simplification in mind, the proof of Thm. 5.0 relies on the concept
of Gaussian equivalents and the linearization analysis developed in Adlam and
Pennington (2020a, b); Adlam et al. (2019), which we briefly review here,
though we refer the reader to these works for a more detailed description. To
proceed, we first define the moment-matched Gaussian linearizations,

$\displaystyle\bar{F}\to$
$\displaystyle\sqrt{\frac{\bar{\rho}}{n_{0}}}WX+\sqrt{\bar{\eta}-\bar{\zeta}}\Theta$
(S316) $\displaystyle\bar{f}\to$
$\displaystyle\sqrt{\frac{\bar{\rho}_{*}}{n_{0}}}W\mathbf{x}+\sqrt{\bar{\eta}_{*}-\bar{\zeta}_{*}}\theta$
(S317)

where $\bar{f}\mathrel{\mathop{:}}=\bar{\sigma}(W\mathbf{x}/\sqrt{n_{0}})$ is
the random feature representation of the test point $\mathbf{x}$ and
$y\mathrel{\mathop{:}}=\sqrt{\frac{1}{n_{0}}}\beta^{\top}\mathbf{x}$ is its
corresponding label, which we assume has no additive noise (including noise on
the test labels merely shifts by the total error by an irreducible additive
constant). The centered parameters141414With some abuse of notation, we use
$\bar{\sigma}(\cdot)$ to denote the centering with respect to the training
scale $s$ and $s_{*}$ for unasterisked and asterisked quantities respectively.
are defined as, $\bar{\zeta}$ and $\bar{\zeta}_{*}$ as
$\bar{\zeta}=s\bar{\rho}$ and $\bar{\zeta}_{*}=s_{*}\bar{\rho}_{*}$ with,

$\displaystyle\bar{\eta}$
$\displaystyle\mathrel{\mathop{:}}=\mathbb{E}_{z\sim\mathcal{N}(0,s)}[\bar{\sigma}(z)^{2}]\,,\;\;$
$\displaystyle\bar{\rho}$
$\displaystyle\mathrel{\mathop{:}}=(\frac{1}{s}\mathbb{E}_{z\sim\mathcal{N}(0,s)}[z\bar{\sigma}(z)])^{2}\,,$
(S318) $\displaystyle\bar{\eta}_{*}$
$\displaystyle\mathrel{\mathop{:}}=\mathbb{E}_{z\sim\mathcal{N}(0,s_{*})}[\bar{\sigma}(z)^{2}]\,,\;\;$
$\displaystyle\bar{\rho}_{*}$
$\displaystyle\mathrel{\mathop{:}}=(\frac{1}{s_{*}}\mathbb{E}_{z\sim\mathcal{N}(0,s_{*})}[z\bar{\sigma}(z)])^{2}\,.$
(S319)

The new objects $\Theta$, $\theta$ are matrices of the appropriate shapes with
i.i.d. standard Gaussian entries independent of the other random variables
under consideration. The constants $\bar{\eta}$, $\bar{\zeta}$,
$\bar{\eta}_{*}$, $\bar{\zeta}_{*}$ are chosen so that the mixed moments up to
second order are the same for the original and linearized matrices. The
Gaussian equivalents defined are essentially constructed via a Taylor
expansion of the nonlinearity $\bar{\sigma}(\cdot)$. The explicit calculations
defining these expressions can be found via the Gaussian moment matching
technique in Adlam et al. (2019).

Having appropriately linearized the random feature matrices, we can map back
to the definition of the original activation functions by recalling that the
variables $\zeta$ and $\zeta_{*}$ are defined in Sec. 5.1 as $\zeta=s\rho$ and
$\zeta_{*}=s_{*}\rho_{*}$, and

$\displaystyle\eta$
$\displaystyle\mathrel{\mathop{:}}=\bar{\eta}=\mathbb{E}_{z\sim\mathcal{N}(0,s)}[(\sigma(z)-\mathbb{E}_{z\sim\mathcal{N}(0,s)}[\sigma(z)])^{2}]=\mathbb{V}_{z\sim\mathcal{N}(0,s)}[\sigma(z)]$
(S320) $\displaystyle\rho$
$\displaystyle\mathrel{\mathop{:}}=\bar{\rho}=(\frac{1}{s}\mathbb{E}_{z\sim\mathcal{N}(0,s)}[z(\sigma(z)-\mathbb{E}_{z\sim\mathcal{N}(0,s)}[\sigma(z)])])^{2}=(\frac{1}{s}\mathbb{E}_{z\sim\mathcal{N}(0,s)}[z(\sigma(z)])^{2}$
(S321) $\displaystyle\eta_{*}$
$\displaystyle\mathrel{\mathop{:}}=\bar{\eta_{*}}=\mathbb{E}_{z\sim\mathcal{N}(0,s_{*})}[(\sigma(z)-\mathbb{E}_{z\sim\mathcal{N}(0,s_{*})}[\sigma(z)])^{2}]=\mathbb{V}_{z\sim\mathcal{N}(0,s_{*})}[\sigma(z)]$
(S322) $\displaystyle\rho_{*}$
$\displaystyle\mathrel{\mathop{:}}=\bar{\rho_{*}}=(\frac{1}{s_{*}}\mathbb{E}_{z\sim\mathcal{N}(0,s_{*})}[z(\sigma(z)-\mathbb{E}_{z\sim\mathcal{N}(0,s)}[\sigma(z)])])^{2}=(\frac{1}{s_{*}}\mathbb{E}_{z\sim\mathcal{N}(0,s_{*})}[z(\sigma(z)])^{2}\,.$
(S323)

by using the definition of the centering. Effectively, this is equivalent to
using the linearizations,

$\displaystyle F\to\;\bar{F}$
$\displaystyle\to\sqrt{\frac{\rho}{n_{0}}}WX+\sqrt{\eta-\zeta}\Theta$ (S324)
$\displaystyle f\to\;\bar{f}$
$\displaystyle\to\sqrt{\frac{\rho_{*}}{n_{0}}}W\mathbf{x}+\sqrt{\eta_{*}-\zeta_{*}}\theta$
(S325)

directly to compute the test error.

Simply by definition, the teacher function is exactly linear,

$\displaystyle Y=$
$\displaystyle\sqrt{\frac{1}{n_{0}}}\beta^{\top}X+\varepsilon$ (S326)
$\displaystyle y=$
$\displaystyle\sqrt{\frac{1}{n_{0}}}\beta^{\top}\mathbf{x}\,.$ (S327)

We emphasize that the restriction to linear teacher functions is simply a
matter of convenience and simplicity—nonlinear teacher neural networks can be
studied by means of an analogous linearization process, as discussed in Adlam
and Pennington (2020a), which merely requires introducing additional moment-
matching constants and additional i.i.d. standard Gaussian terms $\theta$, and
$\theta_{y}$.

In the high-dimensional limit the bulk statistics we compute defining the
error, bias, and variance are invariant to the above replacements by
linearized Gaussian equivalents. We remark that further intuition for these
linearized information-plus-noise replacements can be gathered from the
universality results of Banna et al. (2015, 2020). To summarize briefly, the
final expressions we compute are tracial (nonlinear) functions of several
random matrices. A large body of universality results stemming from Banna et
al. (2020); Erdos (2019) show that in the limit it is sufficient to match them
to an equivalent tracial functional of a (linearized) rational function of
random matrices, whose moments match their nonlinear counterparts.

More specifically, two basic trace objects arise in the calculations, which
take the form

$\mathrm{tr}(AB)\,,\quad\text{and}\quad\mathrm{tr}\left({A\frac{1}{B-zI}}\right)\,.$
(S328)

For the first case, $\mathrm{tr}(AB)$, if the random matrices $A$ and $B$ are
independent, its asymptotics can be understood via concentration of measure
arguments. Conditionally on $A$, it is sufficient to replace $B$ with an
equivalent matrix designed to match low-order moments, at the expense of an
error that vanishes asymptotically. The second case,
$\mathrm{tr}\left({A\frac{1}{B-zI}}\right)$, is a bit more involved. Our
arguments proceed by 1) using the linearized matrices in Eqs. S324, LABEL: and
S325 to express the trace object as a rational function of i.i.d. Gaussian
matrices, and then 2) using the linear pencil method151515This is essentially
an iterative application of the Schur complement formula. to express it as the
trace of a large (inverted) block matrix (Helton et al., 2018; Mingo and
Speicher, 2017). The reason the linearization over $B$ preserves the
asymptotic statistics even for general $A$ with correlated entries stems from
the matrix Dyson equation (Erdos, 2019).

###  I.3 Decomposition of the test loss

Recall the expression for the test loss in Eq. S143,

$E_{\Sigma^{*}}=E_{1}+E_{2}+E_{3}$ (S329)

with

$\displaystyle E_{1}$
$\displaystyle=\mathbb{E}_{(\mathbf{x},\beta,\varepsilon)}y(\mathbf{x})^{2}=\mathbb{E}_{(\mathbf{x},\varepsilon)}\mathrm{tr}(y(\mathbf{x})y(\mathbf{x})^{\top})$
(S330) $\displaystyle E_{2}$
$\displaystyle=-2\mathbb{E}_{(\mathbf{x},\beta,\varepsilon)}(K_{\mathbf{x}}^{\top}K^{-1}Y)y(\mathbf{x})=-2\mathbb{E}_{(\mathbf{x},\varepsilon)}\mathrm{tr}(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top}y(\mathbf{x}))$
(S331) $\displaystyle E_{3}$
$\displaystyle=\mathbb{E}_{(\mathbf{x},\beta,\varepsilon)}(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top})^{2}=\mathbb{E}_{(\mathbf{x},\varepsilon)}\mathrm{tr}(K_{\mathbf{x}}^{\top}K^{-1}Y^{\top}YK^{-1}K_{\mathbf{x}})\,,$
(S332)

where the kernels $K=K(X,X)$ and $K_{\mathbf{x}}=K(X,\mathbf{x})$ are given
by,

$K=\frac{F^{\top}\\!F}{n_{1}}+\gamma I_{m}\qquad\text{and}\qquad
K_{\mathbf{x}}=\frac{1}{n_{1}}F^{\top}f\,.$ (S333)

Using the cyclicity and linearity of the trace, the expectation over
$\mathbf{x}$ requires the computation of

$\mathbb{E}_{\mathbf{x}}K_{\mathbf{x}}K_{\mathbf{x}}^{\top}\,,\qquad\mathbb{E}_{\mathbf{x}}y(\mathbf{x})K_{\mathbf{x}}^{\top}\,,\qquad\mathbb{E}_{\mathbf{x}}y(\mathbf{x})y(\mathbf{x})^{\top}\,.$
(S334)

As described in Sec. I.2, without loss of generality we consider the case of a
linear teacher, and Eqs. S327, LABEL: and S325 read

$y=\frac{1}{\sqrt{n_{0}}}\beta^{\top}\mathbf{x}\,\qquad\text{and}\qquad f\to
f^{\text{lin}}=\frac{\sqrt{\rho_{*}}}{\sqrt{n_{0}}}W\mathbf{x}+\sqrt{\eta_{*}-\zeta_{*}}\theta\,.$
(S335)

The expectations over $\mathbf{x}$ are now trivial and we readily find,

$\displaystyle\mathbb{E}_{\mathbf{x}}K_{\mathbf{x}}K_{\mathbf{x}}^{\top}$
$\displaystyle=\frac{1}{n_{1}^{2}}F^{\top}\left({\frac{\rho_{*}}{n_{0}}W\Sigma^{*}W^{\top}+(\eta_{*}-\zeta_{*})I_{n_{1}}}\right)F$
(S336)
$\displaystyle\mathbb{E}_{\mathbf{x}}y(\mathbf{x})K_{\mathbf{x}}^{\top}$
$\displaystyle=\frac{\sqrt{\rho_{*}}}{n_{0}n_{1}}\beta^{\top}\Sigma^{*}W^{\top}F$
(S337) $\displaystyle\mathbb{E}_{\mathbf{x}}y(\mathbf{x})y(\mathbf{x})^{\top}$
$\displaystyle=\frac{1}{n_{0}}\beta\Sigma^{*}\beta^{\top}$ (S338)

One may interpret the substitution in Eq. S335 as a tool to calculate the
expectations above to leading order, which generates terms like Eq. S328.
Next, we recall the definition, $Y=\beta^{\top}X/\sqrt{n_{0}}+\varepsilon$,
and, as above, we consider the leading order behavior with respect to the
random variables $\beta$ to find

$\displaystyle\mathbb{E}_{\beta,\varepsilon}\left[{Y^{\top}Y}\right]$
$\displaystyle=\frac{1}{n_{0}}X^{\top}X+\sigma_{\varepsilon}^{2}I_{m}$ (S339)
$\displaystyle\mathbb{E}_{\beta,\varepsilon}\left[{Y^{\top}\mathbb{E}_{\mathbf{x}}y(\mathbf{x})K_{\mathbf{x}}^{\top}}\right]$
$\displaystyle=\frac{\sqrt{\rho_{*}}}{n_{0}^{3/2}n_{1}}X^{\top}\Sigma^{*}W^{\top}F\,.$
(S340)

Putting these pieces together, we have

$\displaystyle E_{1}$ $\displaystyle=\frac{\mathrm{tr}(\Sigma^{*})}{n_{0}}$
(S341) $\displaystyle E_{2}$ $\displaystyle=E_{21}$ (S342) $\displaystyle
E_{3}$ $\displaystyle=E_{31}+E_{32}\,,$ (S343)

where,

$\displaystyle E_{21}$
$\displaystyle=-2\frac{\sqrt{\rho_{*}}}{n_{0}^{3/2}n_{1}}\mathbb{E}\mathrm{tr}\left(X^{\top}\Sigma^{*}W^{\top}FK^{-1}\right)$
(S344) $\displaystyle E_{31}$
$\displaystyle=\sigma_{\varepsilon}^{2}\mathbb{E}\mathrm{tr}\left(K^{-1}\Sigma_{3}K^{-1}\right)$
(S345) $\displaystyle E_{32}$
$\displaystyle=\frac{1}{n_{0}}\mathbb{E}\mathrm{tr}\left(K^{-1}\Sigma_{3}K^{-1}X^{\top}X\right)$
(S346)

and,

$\Sigma_{3}=\frac{\rho_{*}}{n_{0}n_{1}^{2}}F^{\top}W\Sigma^{*}W^{\top}F+\frac{\eta_{*}-\zeta_{*}}{n_{1}^{2}}F^{\top}F\,.$
(S347)

###  I.4 Decomposition of the bias and variance

Note that it is sufficient to calculate the bias term since, given the total
test loss, the total variance can be obtained as
$V_{\Sigma^{*}}=E_{\Sigma^{*}}-B_{\Sigma^{*}}$. Following the total
multivariate bias-variance decomposition of Adlam and Pennington (2020b), for
each random variable in question we introduce an i.i.d. copy of it denoted by
either the subscript $1$ or $2$. We can then write,

$\displaystyle B_{\Sigma^{*}}$
$\displaystyle=\mathbb{E}_{(\mathbf{x},y)}(y-\mathbb{E}_{(W,X,\varepsilon)}\hat{y}(\mathbf{x};W,X,\varepsilon))^{2}$
(S348)
$\displaystyle=\mathbb{E}_{(\mathbf{x},y)}\mathbb{E}_{(W_{1},X_{1},\varepsilon_{1})}\mathbb{E}_{(W_{2},X_{2},\varepsilon_{2})}(y-\hat{y}(\mathbf{x};W_{1},X_{1},\varepsilon_{1}))(y-\hat{y}(\mathbf{x};W_{2},X_{2},\varepsilon_{2}))$
(S349) $\displaystyle=\frac{\mathrm{tr}(\Sigma^{*})}{n_{0}}+E_{21}+H_{000}\,,$
(S350)

where an expression for $E_{21}$ was given previously and $H_{000}$ satisfies

$\displaystyle H_{000}$
$\displaystyle=\mathbb{E}\hat{y}(\mathbf{x};W_{1},X_{1},\varepsilon_{1})\hat{y}(\mathbf{x};W_{2},X_{2},\varepsilon_{2})\,,$
(S351)

where the expectations are over
$\mathbf{x},W_{1},X_{1},\varepsilon_{1},W_{2},X_{2}$, and $\varepsilon_{2}$.
Recalling the definition of $\hat{y}$,

$\hat{y}(\mathbf{x};W,X,\varepsilon)\mathrel{\mathop{:}}=Y(X,\varepsilon)K(X,X;W)^{-1}K(X,\mathbf{x};W)\,$
(S352)

and the techniques described in the previous section, it is straightforward to
analyze the above term. First note we can write,

$\displaystyle\mathbb{E}_{\mathbf{x}}K(X_{1},\mathbf{x};W_{1})K(\mathbf{x},X_{2};W_{2})$
$\displaystyle=\frac{\rho_{*}}{n_{0}n_{1}^{2}}F_{11}^{\top}W_{1}\Sigma^{*}W_{2}^{\top}F_{22}$
(S353)

since the $f$ linearizations use different sources of auxiliary randomness.
Here we have defined $F_{11}\equiv F(W_{1},X_{1})$ and $F_{22}\equiv
F(W_{2},X_{2})$. Now we proceed to calculate $H_{000}$ as

$\displaystyle H_{000}$
$\displaystyle=\mathbb{E}\hat{y}(\mathbf{x};W_{1},X_{1},\varepsilon)\hat{y}(\mathbf{x};W_{2},X_{2},\varepsilon_{2})$
(S354)
$\displaystyle=\mathbb{E}K(\mathbf{x},X_{2};W_{2})K(X_{2},X_{2};W_{2})^{-1}Y(X_{2},\varepsilon_{2})^{\top}Y(X_{1},\varepsilon_{1})K(X_{1},X_{1};W_{1})^{-1}K(X_{1},\mathbf{x};W)$
(S355)
$\displaystyle=\mathbb{E}\mathrm{tr}\big{(}K(X_{2},X_{2};W_{2})^{-1}X_{2}^{\top}X_{1}K(X_{1},X_{1};W_{1})^{-1}K(X_{1},\mathbf{x};W)K(\mathbf{x},X_{2};W_{2})\big{)}$
(S356)
$\displaystyle=\frac{\rho_{*}}{n_{0}^{2}n_{1}^{2}}\mathbb{E}\mathrm{tr}\left(K_{22}^{-1}X_{2}^{\top}X_{1}K_{11}^{-1}F_{11}^{\top}W_{1}\Sigma^{*}W_{2}^{\top}F_{22}\right)$
(S357) $\displaystyle\equiv E_{4}\,,$ (S358)

where in the second-to-last line we have defined $K_{11}\equiv
K(X_{1},X_{1};W_{1})$ and $K_{22}\equiv K(X_{2},X_{2};W_{2})$.

###  I.5 Summary of linearized trace terms

We now summarize the requisite terms needed to compute the total test error,
bias, and variance after using cyclicity of the trace to rearrange several of
them. In the following, we slightly change notation in order to make explicit
the dependence on the population covariance matrices $\Sigma$ and
$\Sigma^{*}$. To be specific, whereas above we assumed that the columns of
$X_{1}$ and $X_{2}$ were drawn from multivariate Gaussians with covariance
$\Sigma$, below we assume that they are drawn from multivariate Gausssians
with identity covariance. This change is equivalent to replaceing
$X_{1}\to\Sigma^{1/2}X_{1}$ and $X_{2}\to\Sigma^{1/2}X_{2}$ in the above
expressions. We use this definition so that $X_{1}$, $X_{2}$, $W_{1}$,
$W_{2}$, and $\Theta$ all have i.i.d. standard Gaussian entries. From the
previous computations, we can now write the requisite terms as,

$\displaystyle\Sigma_{3}$ $\displaystyle=$
$\displaystyle\frac{\rho_{*}}{n_{0}n_{1}^{2}}F_{11}^{\top}W_{1}\Sigma^{*}W_{1}^{\top}F_{11}+\frac{\eta_{*}-\zeta_{*}}{n_{1}^{2}}F_{11}^{\top}F_{11}$
(S359) $\displaystyle E_{21}$ $\displaystyle=$
$\displaystyle-2\frac{\sqrt{\rho_{*}}}{n_{0}^{3/2}n_{1}}\mathrm{tr}\left(X_{1}^{\top}\Sigma^{1/2}\Sigma^{*}W_{1}^{\top}F_{11}K_{11}^{-1}\right)$
(S360) $\displaystyle E_{31}$ $\displaystyle=$
$\displaystyle\sigma_{\varepsilon}^{2}\mathrm{tr}\left(K_{11}^{-1}\Sigma_{3}K_{11}^{-1}\right)$
(S361) $\displaystyle E_{32}$ $\displaystyle=$
$\displaystyle\frac{1}{n_{0}}\mathrm{tr}\left(K_{11}^{-1}\Sigma_{3}K_{11}^{-1}X_{1}^{\top}\Sigma
X_{1}\right)$ (S362) $\displaystyle E_{4}$ $\displaystyle=$
$\displaystyle\frac{\rho_{*}}{n_{0}^{2}n_{1}^{2}}\mathrm{tr}\left(F_{22}K_{22}^{-1}X_{2}^{\top}\Sigma
X_{1}K_{11}^{-1}F_{11}^{\top}W_{1}\Sigma^{*}W_{2}^{\top}\right)$ (S363)
$\displaystyle E_{\Sigma^{*}}$ $\displaystyle=$
$\displaystyle\frac{\mathrm{tr}(\Sigma^{*})}{n_{0}}+E_{21}+E_{31}+E_{32}$
(S364) $\displaystyle B_{\Sigma^{*}}$ $\displaystyle=$
$\displaystyle\frac{\mathrm{tr}(\Sigma^{*})}{n_{0}}+E_{21}+E_{4}$ (S365)
$\displaystyle V_{\Sigma^{*}}$ $\displaystyle=$ $\displaystyle
E_{\Sigma^{*}}-B_{\Sigma^{*}}$ (S366)

In the remainder of this section we use the machinery of operator-valued free
probability (Mingo and Speicher, 2017) and a series of lengthy algebraic
computations to compute the limiting tracial expressions in
$E_{21},E_{31},E_{32},E_{4}$, from which the total test error, bias, and
variance can be reconstructed.

###  I.6 Calculation of error terms

To compute the test error, bias, and total variance, we need to evaluate the
asymptotic trace objects appearing in the expressions for $E_{21}$, $E_{31}$,
$E_{32}$, and $E_{4}$, defined in the previous section. As these expressions
are essentially rational functions of the random matrices $X$, $W$, $\Theta$,
$\Sigma$, and $\Sigma^{*}$, these computations can be accomplished by
constructing a linear pencil (Far et al., 2006) and using the theory of
operator-valued free probability (Mingo and Speicher, 2017). These techniques
and their application to problems of this type have been well-established
elsewhere (Adlam et al., 2019; Adlam and Pennington, 2020a, b), we only sketch
the mathematical details, referring the reader to the literature for a more
pedagogical overview. Instead, we focus on presenting the details of the
requisite calculations.

Relative to the prior work of Adlam and Pennington (2020a, b), the main
challenge in the current setting is generalizing the calculations to include
an arbitrary training covariance matrix $\Sigma$. This generalization is
facilitated by the general theory of operator-valued free probability, and in
particular through the subordinated form of the operator-valued self-
consistent equations that we first present in Eq. S397. The form of this
equation enables the simple computation of the operator-valued R-transform of
the remaining random matrices, $W$, $X$, and $\Theta$, which are all i.i.d.
Gaussian and can therefore be obtained simply by using the methods of Far et
al. (2006). The remaining complication amounts to performing the trace in Eq.
S397, which becomes an integral over the LJSD $\mu$ in the limit. While this
might in general lead to a complicated coupling of many transcendental
equations, it turns out that the trascendentality can be entirely factored
into a single scalar fixed-point equation, whose solution we denote by $x$
(see Eq. S420), and the remaining equations are purely algebraic given $x$. To
facilitate this particular simplification, it is necessary to first compute
all of the entries in the operator-valued Stieltjes transform of the kernel
matrix $K$, which we do in Sec. I.6.1. Using these results, we compute the
remaining error terms in the subsequent sections.

As a matter of notation, note that throughout this entire section whenever a
matrix $X$, $X_{1}$, or $X_{2}$ appears it is composed of i.i.d. standard
Gaussian entries as in Sec. I.5. This differs from the notation of the main
paper, but we follow this prescription to ease the already cumbersome
presentation. This definition of $X$ allows us to explicitly extract and
represent the training covariance $\Sigma$ in our calculations.

####  I.6.1 $K^{-1}$

Define the block matrix $Q^{K^{-1}}$ as,

$Q^{K^{-1}}=\left(\begin{array}[]{cccccc}I_{m}&\frac{\sqrt{\eta-\zeta}\Theta^{\top}}{\gamma\sqrt{n_{1}}}&\frac{\sqrt{\rho}X^{\top}}{\gamma\sqrt{n_{0}}}&0&0&0\\\
-\frac{\Theta\sqrt{\eta-\zeta}}{\sqrt{n_{1}}}&I_{n_{1}}&0&0&-\frac{\sqrt{\rho}W}{\sqrt{n_{1}}}&0\\\
0&0&I_{n_{0}}&-\Sigma^{1/2}&0&0\\\
0&-\frac{W^{\top}}{\sqrt{n_{1}}}&0&I_{n_{0}}&0&0\\\
0&0&0&0&I_{n_{0}}&-\Sigma^{1/2}\\\
-\frac{X}{\sqrt{n_{0}}}&0&0&0&0&I_{n_{0}}\\\ \end{array}\right)\,.$ (S367)

Then block matrix inversion (i.e. repeated applications of the Schur
complement formula) shows that,

$\displaystyle G^{K^{-1}}_{1,1}$ $\displaystyle=$
$\displaystyle\gamma\,\bar{\mathrm{tr}}(K^{-1})$ (S368) $\displaystyle
G^{K^{-1}}_{2,2}$ $\displaystyle=$
$\displaystyle\gamma\,\bar{\mathrm{tr}}(\hat{K}^{-1})$ (S369) $\displaystyle
G^{K^{-1}}_{3,3}$ $\displaystyle=$ $\displaystyle
G^{K^{-1}}_{6,6}=1-\frac{\sqrt{\rho}\,\bar{\mathrm{tr}}\left(\Sigma^{1/2}W^{\top}FK^{-1}X^{\top}\right)}{\sqrt{n_{0}}n_{1}}$
(S370) $\displaystyle G^{K^{-1}}_{4,3}$ $\displaystyle=$ $\displaystyle
G^{K^{-1}}_{6,5}=\,\bar{\mathrm{tr}}(\Sigma^{1/2})-\frac{\sqrt{\rho}\,\bar{\mathrm{tr}}\left(\Sigma
W^{\top}FK^{-1}X^{\top}\right)}{\sqrt{n_{0}}n_{1}}$ (S371) $\displaystyle
G^{K^{-1}}_{5,3}$ $\displaystyle=$ $\displaystyle
G^{K^{-1}}_{6,4}=\frac{\gamma\sqrt{\rho}\,\bar{\mathrm{tr}}\left(\Sigma^{1/2}W^{\top}\hat{K}^{-1}W\right)}{n_{1}}$
(S372) $\displaystyle G^{K^{-1}}_{6,3}$ $\displaystyle=$
$\displaystyle\frac{\gamma\sqrt{\rho}\,\bar{\mathrm{tr}}\left(\Sigma
W^{\top}\hat{K}^{-1}W\right)}{n_{1}}$ (S373) $\displaystyle G^{K^{-1}}_{3,4}$
$\displaystyle=$ $\displaystyle
G^{K^{-1}}_{5,6}=-\frac{\sqrt{\rho}\,\bar{\mathrm{tr}}\left(FK^{-1}X^{\top}W^{\top}\right)}{\sqrt{n_{0}}n_{1}\psi}$
(S374) $\displaystyle G^{K^{-1}}_{4,4}$ $\displaystyle=$ $\displaystyle
G^{K^{-1}}_{5,5}=1-\frac{\sqrt{\rho}\,\bar{\mathrm{tr}}\left(\Sigma^{1/2}W^{\top}FK^{-1}X^{\top}\right)}{\sqrt{n_{0}}n_{1}}$
(S375) $\displaystyle G^{K^{-1}}_{5,4}$ $\displaystyle=$
$\displaystyle\frac{\gamma\sqrt{\rho}\,\bar{\mathrm{tr}}\left(\hat{K}^{-1}WW^{\top}\right)}{n_{1}\psi}$
(S376) $\displaystyle G^{K^{-1}}_{3,5}$ $\displaystyle=$ $\displaystyle
G^{K^{-1}}_{4,6}=-\frac{\sqrt{\rho}\,\bar{\mathrm{tr}}\left(\Sigma^{1/2}XK^{-1}X^{\top}\right)}{n_{0}}$
(S377) $\displaystyle G^{K^{-1}}_{4,5}$ $\displaystyle=$
$\displaystyle-\frac{\sqrt{\rho}\,\bar{\mathrm{tr}}\left(\Sigma
XK^{-1}X^{\top}\right)}{n_{0}}$ (S378) $\displaystyle G^{K^{-1}}_{3,6}$
$\displaystyle=$
$\displaystyle-\frac{\sqrt{\rho}\,\bar{\mathrm{tr}}\left(K^{-1}X^{\top}X\right)}{n_{0}\phi}\,,$
(S379)

where $G^{K^{-1}}_{i,j}$ denotes the normalized trace of the $(i,j)$-block of
the inverse of $Q^{K^{-1}}$, and we have defined
$\hat{K}=\frac{1}{n_{1}}FF^{\top}+\gamma I_{n_{1}}$.

We aim to compute the limiting values of these trace terms as
$n_{0},n_{1},m\to\infty$, as they will be related to the error terms of
interest. Both here and in the sequel, to ease the already cumbersome
presentation, we use $G$ to denote the limiting values as well as the non-
limiting values and we refrain from explicitly denoting the limit operation
itself, noting its existing can be inferred from context.

To proceed, recall that the asymptotic block-wise traces of the inverse of
$Q^{K^{-1}}$ can be determined from its operator-valued Stieltjes transform
(Mingo and Speicher, 2017). The simplest way to apply the results of Far et
al. (2006); Mingo and Speicher (2017) is to augment $Q^{K^{-1}}$ to form the
the self-adjoint matrix $\bar{Q}^{K^{-1}}$,

$\bar{Q}^{K^{-1}}=\left(\begin{array}[]{cc}0&{[Q^{K^{-1}}]}^{\top}\\\
Q^{K^{-1}}&0\end{array}\right)\,,\\\ $ (S380)

and observe that we can write $\bar{Q}^{K^{-1}}$ as,

$\begin{split}\bar{Q}^{K^{-1}}&=\bar{Z}+\bar{Q}^{K^{-1}}_{W,X,\Theta}+\bar{Q}^{K^{-1}}_{\Sigma}\\\
&=\left(\begin{array}[]{cc}0&Z^{\top}\\\
Z&0\end{array}\right)+\left(\begin{array}[]{cc}0&{[Q^{K^{-1}}_{W,X,\Theta}]}^{\top}\\\
Q^{K^{-1}}_{W,X,\Theta}&0\end{array}\right)+\left(\begin{array}[]{cc}0&{[Q^{K^{-1}}_{\Sigma}]}^{\top}\\\
Q^{K^{-1}}_{\Sigma}&0\end{array}\right)\,,\end{split}$ (S381)

where $Z=I_{m+4n_{0}+n_{1}}$, and,

$\displaystyle Q^{K^{-1}}_{W,X,\Theta}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccccc}0&\frac{\sqrt{\eta-\zeta}\Theta^{\top}}{\gamma\sqrt{n_{1}}}&\frac{\sqrt{\rho}X^{\top}}{\gamma\sqrt{n_{0}}}&0&0&0\\\
-\frac{\Theta\sqrt{\eta-\zeta}}{\sqrt{n_{1}}}&0&0&0&-\frac{\sqrt{\rho}W}{\sqrt{n_{1}}}&0\\\
0&0&0&0&0&0\\\ 0&-\frac{W^{\top}}{\sqrt{n_{1}}}&0&0&0&0\\\ 0&0&0&0&0&0\\\
-\frac{X}{\sqrt{n_{0}}}&0&0&0&0&0\\\ \end{array}\right)$ (S388) $\displaystyle
Q^{K^{-1}}_{\Sigma}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccccc}0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\
0&0&0&-\Sigma^{1/2}&0&0\\\ 0&0&0&0&0&0\\\ 0&0&0&0&0&-\Sigma^{1/2}\\\
0&0&0&0&0&0\\\ \end{array}\right)\,.$ (S395)

Note that we have separated the i.i.d. Gaussian matrices $W,X,\Theta$ from the
constant terms and from the $\Sigma$-dependent terms. Denote by
$\bar{G}^{K^{-1}}$ the block matrix

$\bar{G}^{K^{-1}}=\left(\begin{array}[]{cc}0&{[G^{K^{-1}}]}^{\top}\\\