$\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}=(\boldsymbol{\theta}_{0}+\boldsymbol{u})(\boldsymbol{\theta}_{0}+\boldsymbol{u})^{\top}=\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}+\boldsymbol{\theta}_{0}\boldsymbol{u}^{\top}+\boldsymbol{u}\boldsymbol{\theta}_{0}^{\top}+\boldsymbol{u}\boldsymbol{u}^{\top}$
(98)

To reduce the attack strength, we want to minimise the the attack overlap $A$
which is
$\boldsymbol{\theta}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}$.

To minimise the term
$\boldsymbol{\theta}^{\top}\boldsymbol{u}\boldsymbol{u}^{\top}\boldsymbol{\theta}$,
we would like to make the ERM procedure chose a vector $\boldsymbol{\theta}$
that is perpendicular to $\boldsymbol{u}$. This can be done without affecting
generalisation performances as $\boldsymbol{u}\perp\boldsymbol{\theta}_{0}$. A
way to remove this, is to regularise in the direction $\boldsymbol{u}$. This
can be done with a directional $\ell_{2}$ regularisation or with adversarial
training and choosing
$\boldsymbol{\Sigma}_{\boldsymbol{\delta}}=\boldsymbol{u}\boldsymbol{u}^{\top}$.
With this choice of $\boldsymbol{\Sigma}_{\boldsymbol{\delta}}$, the values
$\boldsymbol{\theta}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{u}^{\top}\boldsymbol{\theta}$
and
$\boldsymbol{\theta}^{\top}\boldsymbol{u}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\theta}$
can be decreased in the same way.

For the last term
$\boldsymbol{\theta}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\theta}$,
this recipe does not work. We can interpret this term as being proportional to
the square of the overlap $m$ for isotropic data. Thus, it is obvious that
protecting such a direction, can only be achieved by reducing the norm of
$\boldsymbol{\theta}$, which reduces the overlap $m$, which in turn cannot be
changed without hurting the angle between teacher and student and thus
generalisation error. This term is thus at the heart of the trade-off between
generalisation and boundary term.

Thus shown that all the components of the eigenvectors of
$\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}$ that are perpendicular to the
teacher can be “regularised away”.

Table 1: We isolate robust features and demonstrate that attacking the non-
robust features contributes most to boundary error. For CIFAR10, we keep the
first 2500 principle components, and evaluate the target 1 vs 9. For
FashionMNIST, we keep the first 500 principle components and evaluate odd
versus even targets.

Cifar 10
---
Attack Direction | $E_{\mathrm{train}}$ | $E_{\mathrm{gen}}$ | $E_{\mathrm{adv}}$
Uniform | 0.0749 | 0.2225 | 0.2645
Robust | 0.0749 | 0.2225 | 0.2230
Non-Robust | 0.0749 | 0.2225 | 0.2645

Fashion MNIST
---
Attack Direction | $E_{\mathrm{train}}$ | $E_{\mathrm{gen}}$ | $E_{\mathrm{adv}}$
Uniform | 0.0326 | 0.0397 | 0.0632
Robust | 0.0326 | 0.0397 | 0.0451
Non-Robust | 0.0326 | 0.0397 | 0.0578

Table 2: We protect the non-robust features and find that identity and inverse
proportionality to the explained variance is good compared to protecting
proportional to the explained variance.

Cifar 10
---
Protection Direction | $E_{\mathrm{train}}$ | $E_{\mathrm{gen}}$ | $E_{\mathrm{adv}}$
Uniform | 0.0801 | 0.1940 | 0.2300
Proportional | 0.0742 | 0.2205 | 0.2630
Inverse | 0.0882 | 0.2040 | 0.2270

Fashion MNIST
---
Protection Direction | $E_{\mathrm{train}}$ | $E_{\mathrm{gen}}$ | $E_{\mathrm{adv}}$
Uniform | 0.0351 | 0.0391 | 0.0515
Proportional | 0.0350 | 0.0388 | 0.0548
Inverse | 0.0342 | 0.0391 | 0.0521

##  Appendix G Experiments on Real Data

In this Appendix we demonstrate the robustness metric and defense-strategies
evaluated the Cifar 10 [28] and FashionMNIST [76] datasets.

###  G.1 Finding Non-Robust Features

The margin-based robustness given in Eq. 9 can be interpreted in terms of the
student estimate $\boldsymbol{\theta}$ and it can be written as

$\mathcal{R}^{\mathrm{emp.}}_{j}=\sum_{i=1}^{n}y_{i}\boldsymbol{\theta}_{j}(\boldsymbol{x}_{i})_{j}-\frac{\varepsilon_{g}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}}$
(99)

which can be evaluated component-wise.

Our aim is to apply this measure on components obtained by PCA. The motivation
for choosing PCA as a method of feature extraction and dimensionality
reduction is twofold, on the one hand it creates uncorrelated variables and on
the other hand, it creates them in order of explained variance. This produces
a data-set similar to the BFM with a power-law distributed spectrum.

We begin our procedure by computing the principal components and projecting
all test and training data. Next, we optimise the weights with
$\varepsilon_{t}=0$. Next, we increase $\varepsilon_{t}$ and evaluate the
component-wise robustness-measure Eq. 99.

By hand, we choose a cutoff-value ($2000$ for Cifar10 and for $4000$
FashionMNIST) for the component-wise robustness-measure to produce a mask for
robust features. All values higher than the cutoff are considered to be
robust. We choose the constant such that the boundary error is explained to a
large degree through the non-robust features. This is evaluated by
constraining the attack to all features, the robust features and the non-
robust features. The results of this experiment is shown in Table 1.

###  G.2 Training the Non-Robust Features Adversarially

A principle finding of ours is the discussion of different defence strategies.
Here, we provide some simple experiments on real data.

In the main text, we suggested to defend the non-robust features either
proportional to their variance, inverse proportional to their variance or
uniformly. We continue the experiment from the previous paragraph and consider
the PCA explained variance as the relevant metric for the proportionality. To
be explicit, we perform three different adversarial training runs with
different directional choices for $\boldsymbol{\Sigma}_{\boldsymbol{\delta}}$.
First, we mask the defence by setting the directions for robust features to
zero. Next, we define the uniform defence by putting all remaining directions
to 1, the proportional defence by setting the directions proportional to their
explained variance and the inverse defence by setting the directions
proportional to the inverse of the explained variance. The results of this
experiment are given in Table 2. Indeed, protecting inversely or uniformly
appears to work well.

##  Appendix H Settings of Figures in Main Text

In this Appendix we provide the parameters for the figures presented in the
main text.

###  H.1 Figure 2

In Fig. 2 we consider we consider a SWFM with a single block. Considering the
figures from left to right, the entries for the parameter specifying the block
are

Table 3: Parameters for Fig. 2
Plot Label | $\psi_{\ell}$ | $t_{\ell}$ | $\Delta_{\ell}$ | $\Upsilon_{\ell}$ | $\varepsilon_{g}$ | $\lambda$ | $\tau$
---|---|---|---|---|---|---|---
Low Robustness/Low Usefulness | 0.5 | 2 | 1 | 1 | 0.2 | $10^{-3}$ | 0.05
Low Robustness/High Usefulness | 0.5 | 8 | 1 | 1 | 0.2 | $10^{-3}$ | 0.05
High Robustness/Low Usefulness | 2 | 0.5 | 1 | 1 | 0.2 | $10^{-3}$ | 0.05
High Robustness/High Usefulness | 2 | 2 | 1 | 1 | 0.2 | $10^{-3}$ | 0.05

The points describe the simulation using ERM routine explained in Appendix J,
where we have used a dimension $d=1000$ and averaged over 20 runs. We show the
mean and the standard deviation on the mean.

###  H.2 Figure 3

In Fig. 3 we build a SWFM model where we fix two different blocks of features.
The percentage for the two features are $(\phi_{1},\phi_{2})=(0.01,0.99)$. The
experiment is performed at a high sample complexity $\alpha=10^{2}$. The other
values for the plots are

Table 4: Parameters for Fig. 3
Plot Label | $(\psi_{1},\psi_{2})$ | $(t_{1},t_{2})$ | $(\Delta_{1},\Delta_{2})$ | $(\Upsilon_{1},\Upsilon_{2})$ | $\varepsilon_{g}$ | $\lambda$ | $\tau$
---|---|---|---|---|---|---|---
Left | (1,1) |  ($10^{3}$,1) |  ($10^{3}$,1) |  ($10^{3}$,1) | 0.006 | $10^{-3}$ | 0.05
Right | (1,1) |  ($10^{3}$,1) | $\operatorname{diag}(\boldsymbol{\Sigma}_{\boldsymbol{\delta}}^{b})$ | $\operatorname{diag}(\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}^{b})$ | 0.006 | $10^{-3}$ | 0.05

In the right plot, we compose a vector orthogonal to the diagonal to the
teacher covariance
$\boldsymbol{v}\perp\operatorname{diag}\boldsymbol{\Sigma}_{\boldsymbol{\theta}}$.
Then, we choose
$\boldsymbol{\Sigma}_{\boldsymbol{\delta}}^{b}=\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}^{b}=a\boldsymbol{v}\boldsymbol{v}^{\top}+N$.
We add a variance 1 Gaussian noise matrix $N$ to make the matrix non-singular
and we increase the importance of the first component through a factor $a$.
All matrices are normalised by their trace.

###  H.3 Figure 4

In Fig. 4 there are three blocks of images, we start from the first block on
the left. Here we show an SWFM consisting of two blocks $k=2$. The values are
filled as follows

Table 5: Parameters for Fig. 4 (Left)
Plot Label | $(\psi_{1},\psi_{2})$ | $(t_{1},t_{2})$ | $(\Delta_{1},\Delta_{2})$ | $(\Upsilon_{1},\Upsilon_{2})$ | $\varepsilon_{g}$ | $\lambda$ | $\tau$
---|---|---|---|---|---|---|---
Protecting Robust | (5,0.2) | (1,1) | (2,1) | (1,1) | 0.2 | $10^{-3}$ | 0.05
Uniform Defence | (5,0.2) | (1,1) | (1,1) | (1,1) | 0.2 | $10^{-3}$ | 0.05
Protecting Non-Robust | (5,0.2) | (1,1) | (1,2) | (1,1) | 0.2 | $10^{-3}$ | 0.05

All matrices are normalised by their trace. The points describe the simulation
using ERM routine explained in Appendix J, where we have used a dimension
$d=1000$ and averaged over 20 runs. We show the mean and the standard
deviation on the mean.

The figure on the top right, shows a power-law spectrum in $\psi_{\ell}$, we
begin the indexing with 1 and choose a coefficient $\beta=1.5$. Again, all the
matrices are normalised by their trace. The values are constructed as follows,
we consider the block dimension $d_{\ell}=1$, and use the following values

Table 6: Parameters for Fig. 4 (Top, Right) and (Bottom, Left)
$\psi_{\ell}$ | $t_{\ell}$ | $\Delta_{\ell}$ | $\Upsilon_{\ell}$ | $\varepsilon_{g}$ | $\lambda$ | $\tau$
---|---|---|---|---|---|---
$i^{-\beta}$ | 1 | 1 | 1 | 0.2 | $10^{-3}$ | 0.05

In the figure on the bottom right, we show different power-law behaviour when
$\alpha\to\infty$. Again, we normalise each finite size power-law by their
trace, we fix $d=1000$. The parameters are filled exactly as in the previously
described figure, the values can be found in Table 6.

###  H.4 Figure 5

In Fig. 5 for both the left and right plot we considered a BFM with a single
block $k=1$. We fix $\tau=0.05$ and $\varepsilon_{g}=0.2$.

In both figures, we consider $t_{\ell}=1$, $\psi_{\ell}=1$, $\Delta_{\ell}=1$,
$\Upsilon_{i}=1$ and the values for the $\varepsilon_{t}$ are shown in the
legend.

In the left figure, we use our state evolution equations to optimise the
regularisation strength in terms of generalisation error. Each simulation
point hence is created at its respective optimal regularisation.

In the right figure, we compare the learning curves of standard adversarial
training with the learning curves of the approximately equivalent loss, for
this figure we fix the regularisation strength $\lambda=10^{-3}$.

##  Appendix I Replica Computation

In this Section we give a full derivation of the results in Theorem III.1 by
means of the replica approach, a standard method developed in the realm of
statistical physics of disordered systems. Our computational approach can be
found in [32, 16, 3]. For a foundational understanding of this effective yet
heuristic method, we suggest the following books [37, 38].

###  I.1 Gibbs minimisation

The starting point is to define the following Gibbs measure over weights
$\boldsymbol{\theta}\in\mathbb{R}^{d}$. We want the most probable states to be
the ones that minimises the ERM problem in the first place. Then to select
only these states we will consider the zero temperature limit by taking the
parameter $\beta\to\infty$. The measure that we are interested in is

$\mu_{\beta}(\mathrm{d}\boldsymbol{\theta})=\frac{1}{\mathcal{Z}_{\beta}}e^{-\beta\left[\sum_{\mu=1}^{n}g\left(y^{\mu},\boldsymbol{\theta}^{\top}\boldsymbol{x}^{\mu},\boldsymbol{\theta},\boldsymbol{\Sigma}_{\boldsymbol{\delta}},\varepsilon_{t}\right)+\frac{\lambda}{2}\boldsymbol{\theta}^{T}\boldsymbol{\Sigma}_{\boldsymbol{w}}\boldsymbol{\theta}\right]}\mathrm{d}\boldsymbol{\theta}=\frac{1}{\mathcal{Z}_{\beta}}\underbrace{\prod_{\mu=1}^{n}e^{-\beta
g\left(y^{\mu},\boldsymbol{\theta}^{\top}\boldsymbol{x}^{\mu},\boldsymbol{\theta},\boldsymbol{\Sigma}_{\boldsymbol{\delta}},\varepsilon_{t}\right)}}_{P_{g}}\underbrace{e^{-\frac{\beta\lambda}{2}\boldsymbol{\theta}^{T}\boldsymbol{\Sigma}_{\boldsymbol{w}}\boldsymbol{\theta}}\mathrm{~{}d}w_{i}}_{P_{w}}$
(100)

Where $P_{g}$ is the probability distribution associated with the channel and
$P_{w}$ is the prior probability distribution.

Here, $\mathcal{Z}_{\beta}$, is the partition function that normalises the
Gibbs measure and it is given by

$\mathcal{Z}_{\beta}=\int_{\mathbb{R}^{d}}\mathrm{d}\boldsymbol{\theta}e^{-\frac{\beta\lambda}{2}\boldsymbol{\theta}^{T}\boldsymbol{\Sigma}_{\boldsymbol{w}}\boldsymbol{\theta}}\prod_{\mu=1}^{n}e^{-\beta
g\left(y^{\mu},\boldsymbol{\theta}^{\top}\boldsymbol{x}^{\mu},\boldsymbol{\theta},\boldsymbol{\Sigma}_{\boldsymbol{\delta}},\varepsilon_{t}\right)}$
(101)

You do need attention, but the free energy density is sufficient. In the zero
temperature limit, $\beta\to\infty$ the Gibbs measure in Eq. 100 concentrates
around the solutions of the ERM problem. With the replica method, we can
compute the free energy density, it is given by:

$\beta
f_{\beta}=-\lim_{d\rightarrow\infty}\frac{1}{d}\mathbb{E}_{\mathcal{D}}\log\mathcal{Z}_{\beta}$
(102)

To evaluate the quenched average of the free energy is to use the replica
trick

$\lim_{d\rightarrow\infty}\frac{1}{d}\mathbb{E}_{\mathcal{D}}\log\mathcal{Z}_{\beta}=\lim_{r\rightarrow
0}\lim_{d\rightarrow\infty}\frac{1}{d}\frac{\partial_{r}\mathbb{E}_{\mathcal{D}}\mathcal{Z}^{r}}{1}$
(103)

Note that we introduced three limits up to here. The first is the zero
temperature limit ensuring that we find the ground state of our Gibbs measure
which corresponds to the minimum of our ERM problem. The second is the
thermodynamic limit of very large dimension whilst keeping the sampling ratio
fixed. And the third limit stems from the replica trick allowing us to compute
the logarithm of the partition function, it corresponds to setting the number
of replicated systems to zero.

This computation follows for the first part the one in [32]. So we start with
the initial definition of replicated partition function the difference we have
in our case is that we have a dependence on $\varepsilon_{t}$ on the output
probability.

$\displaystyle\mathbb{E}_{\mathcal{D}}\mathcal{Z}_{\beta}^{r}$
$\displaystyle=\prod_{\mu=1}^{n}\mathbb{E}_{\boldsymbol{x}^{\mu}}\prod_{a=1}^{r}\int_{\mathbb{R}^{d}}P_{w}\left(\differential{\boldsymbol{\theta}^{a}}\right)P\left(y^{\mu}\mid\frac{\boldsymbol{x}^{\mu}\cdot\boldsymbol{\theta}^{a}}{\sqrt{d}}\right)$
(104)
$\displaystyle=\prod_{\mu=1}^{n}\int_{\mathbb{R}}\differential{y}^{\mu}\int_{\mathbb{R}^{p}}P_{\boldsymbol{\theta}_{0}}\left(\differential{\boldsymbol{\theta}_{0}}\right)\int_{\mathbb{R}^{d\times
r}}\left(\prod_{a=1}^{r}P_{w}\left(\differential{\boldsymbol{\theta}^{a}}\right)\right)\mathbb{E}_{\boldsymbol{x}^{\mu}}\left[P_{0}\left(y^{\mu}\mid\frac{\boldsymbol{x}^{\mu}\cdot\boldsymbol{\theta}_{0}}{\sqrt{d}}\right)\prod_{a=1}^{r}P_{g}\left(y^{\mu}\mid\frac{\boldsymbol{x}^{\mu}\cdot\boldsymbol{\theta}^{a}}{\sqrt{d}},\boldsymbol{\Sigma}_{\boldsymbol{\delta}},\boldsymbol{\theta}^{a},\varepsilon_{t}\right)\right]$

explicitly we have that the term in $P_{g}$ is

$P_{g}\left(y^{\mu}\mid\frac{\boldsymbol{x}^{\mu}\cdot\boldsymbol{\theta}^{a}}{\sqrt{d}},\boldsymbol{\Sigma}_{\boldsymbol{\delta}},\boldsymbol{\theta}^{a},\varepsilon_{t}\right)=\frac{\sqrt{\beta}}{\sqrt{2\pi}}\exp\quantity(-\beta
g\quantity(y\frac{\boldsymbol{x}^{\mu}\cdot\boldsymbol{\theta}^{a}}{\sqrt{d}}-\frac{\varepsilon_{t}}{\sqrt{d}}\frac{\boldsymbol{\theta}^{a}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}^{a}}{\norm{\boldsymbol{\theta}^{a}}_{2}}))$
(105)

note that $P_{0}$ can be a general noisy channel distribution. the expectation
part is equal to:

$\displaystyle\mathbb{E}_{\boldsymbol{x}^{\mu}}\left[P_{0}\left(y^{\mu}\mid\frac{\boldsymbol{x}^{\mu}\cdot\boldsymbol{\theta}_{0}}{\sqrt{d}}\right)\prod_{a=1}^{r}P_{g}\left(y^{\mu}\mid\frac{\boldsymbol{x}^{\mu}\cdot\boldsymbol{\theta}^{a}}{\sqrt{d}},\boldsymbol{\Sigma}_{\boldsymbol{\delta}},\boldsymbol{\theta}^{a},\varepsilon_{t}\right)\right]$
(106)
$\displaystyle=\int_{\mathbb{R}}\mathrm{d}\nu_{\mu}P_{0}\left(y\mid\nu_{\mu}\right)\int_{\mathbb{R}^{r}}\left(\prod_{a=1}^{r}\mathrm{~{}d}\lambda_{\mu}^{a}P_{g}\left(y^{\mu}\mid\lambda_{\mu}^{a},\boldsymbol{\Sigma}_{\boldsymbol{\delta}},\boldsymbol{\theta}^{a},\varepsilon_{t}\right)\right)\mathbb{E}_{\boldsymbol{x}^{\mu}}\left[\delta\left(\nu_{\mu}-\frac{\boldsymbol{x}^{\mu}\cdot\boldsymbol{\theta}_{0}}{\sqrt{d}}\right)\prod_{a=1}^{r}\delta\left(\lambda_{\mu}^{a}-\frac{\boldsymbol{x}^{\mu}\cdot\boldsymbol{\theta}^{a}}{\sqrt{d}}\right)\right]$

We can still perform the average over the dataset. We have that the new
variables will behave again as Gaussians with the following covariances:

$\rho\equiv\mathbb{E}\left[\nu_{\mu}^{2}\right]=\frac{1}{d}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}_{0},\quad
m^{a}\equiv\mathbb{E}\left[\lambda_{\mu}^{a}\nu_{\mu}\right]=\frac{1}{d}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{a},\quad
Q^{ab}\equiv\mathbb{E}\left[\lambda_{\mu}^{a}\lambda_{\mu}^{b}\right]=\frac{1}{d}\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{b}$
(107)

where one can organise them in a single covariance matrix.

Now we want to perform several change of variables. The first one is the one
in the matrix of overlaps:

$\displaystyle 1\propto$
$\displaystyle\int_{\mathbb{R}}\differential{\rho}\delta\left(d\rho-\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}_{0}\right)\int_{\mathbb{R}^{r}}\prod_{a=1}^{r}\differential{m^{a}}\delta\left(dm^{a}-\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{a}\right)\int_{\mathbb{R}^{r\times
r}}\prod_{1\leq a\leq b\leq
r}\differential{Q^{ab}}\delta\left(dQ^{ab}-\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{b}\right)$
(108) $\displaystyle=$
$\displaystyle\int_{\mathbb{R}}\frac{\differential{\rho}\mathrm{d}\hat{\rho}}{2\pi}e^{-i\hat{\rho}\left(d\rho-\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}_{0}\right)}\int_{\mathbb{R}^{r}}\prod_{a=1}^{r}\frac{\differential{m^{a}}\differential{\hat{m}^{a}}}{2\pi}e^{-i\sum_{a=1}^{r}\hat{m}^{a}\left(dm^{a}-\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{a}\right)}\int_{\mathbb{R}^{r\times
r}}\prod_{1\leq a\leq b\leq
r}\frac{\mathrm{d}Q^{ab}\differential{\hat{Q}^{ab}}}{2\pi}e^{-i\hat{Q}^{ab}\left(dQ^{ab}-\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{b}\right)}$

We also would like to define new overlaps which are

$P=\frac{1}{d}\boldsymbol{\theta}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}\,,\quad
A=\frac{1}{d}\boldsymbol{\theta}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}\,,\quad
F=\frac{1}{d}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}\,,\quad
N=\frac{1}{d}\boldsymbol{\theta}^{\top}\boldsymbol{\theta}\,,$ (109)

which enter into the computation as follows

$\displaystyle 1\propto$
$\displaystyle\int\prod_{a=1}^{r}\differential{P^{a}}\delta\quantity(dP^{a}-\boldsymbol{\theta}^{a}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}^{a})\int\prod_{a=1}^{r}\differential{A^{a}}\delta\quantity(dA^{a}-\boldsymbol{\theta}^{a}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}^{a})$
(110)
$\displaystyle\int_{\mathbb{R}^{r}}\prod_{a=1}^{r}\differential{F^{a}}\delta\left(dF^{a}-\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}^{a}\right)\int\prod_{a=1}^{r}\differential{N^{a}}\delta\quantity(dN^{a}-\boldsymbol{\theta}^{a}\cdot\boldsymbol{\theta}^{a})$
$\displaystyle=$
$\displaystyle\int\prod_{a=1}^{r}\frac{\differential{P^{a}}\differential{\hat{P}^{a}}}{2\pi}e^{-i\hat{P}^{a}(dP^{a}-\boldsymbol{\theta}^{a}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}^{a})}\int\prod_{a=1}^{r}\frac{\differential{A^{a}}\differential{\hat{A}^{a}}}{2\pi}e^{-i\hat{A}^{a}(dA^{a}-\boldsymbol{\theta}^{a}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}^{a})}$
$\displaystyle\int_{\mathbb{R}^{r}}\prod_{a=1}^{r}\frac{\differential{F^{a}}\differential{\hat{F}^{a}}}{2\pi}e^{-i\sum_{a=1}^{r}\hat{F}^{a}\left(dF^{a}-\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}^{a}\right)}\int\prod_{a=1}^{r}\frac{\differential{N^{a}}\differential{\hat{N}^{a}}}{2\pi}e^{-i\hat{N}^{a}(dN^{a}-\boldsymbol{\theta}^{a}\cdot\boldsymbol{\theta}^{a})}$

We finally can write our replicated partition function as the integral of a
functional as follows

$\mathbb{E}_{\mathcal{D}}\mathcal{Z}_{\beta}^{r}=\int\frac{\differential{\rho}\mathrm{d}\hat{\rho}}{2\pi}\prod_{a=1}^{r}\frac{\differential{m^{a}}\differential{\hat{m}^{a}}}{2\pi}\frac{\differential{P^{a}}\differential{\hat{P}^{a}}}{2\pi}\frac{\differential{A^{a}}\differential{\hat{A}^{a}}}{2\pi}\frac{\differential{F^{a}}\differential{\hat{F}^{a}}}{2\pi}\frac{\differential{N^{a}}\differential{\hat{N}^{a}}}{2\pi}\prod_{1\leq
a\leq b\leq
r}\frac{\mathrm{d}Q^{ab}\differential{\hat{Q}^{ab}}}{2\pi}e^{d\Phi^{(r)}}$
(111)

where the $r$ times replicated functional $\Phi^{(r)}$ is

$\displaystyle\Phi^{(r)}$
$\displaystyle=-\rho\hat{\rho}-\sum_{a=1}^{r}m^{a}\hat{m}^{a}-\sum_{1\leq
a\leq b\leq
r}Q^{ab}\hat{Q}^{ab}-\sum_{a=1}^{r}N^{a}\hat{N}^{a}-\sum_{a=1}^{r}A^{a}\hat{A}^{a}-\sum_{a=1}^{r}P^{a}\hat{P}^{a}-\sum_{a=1}^{r}F^{a}\hat{F}^{a}$
(112)
$\displaystyle+\alpha\Psi_{y}^{(r)}\left(\rho,m^{a},Q^{ab},A^{a},N^{a},P^{A},F^{a}\right)+\Psi_{w}^{(r)}\left(\hat{\rho},\hat{m}^{a},\hat{Q}^{ab},\hat{A}^{a},\hat{N}^{a},\hat{P}^{a},\hat{F}^{a}\right)$

we will refer to the elements in the first line of Eq. 112 as the trace term.
Note in Eq. 111 we factored out $d$ such that we can later evaluate the
partition function in the thermodynamic limit using Laplace’s method. We also
have defined the prior part of the free energy $\Psi_{w}$ to be

$\displaystyle\Psi_{w}^{(r)}=$
$\displaystyle\frac{1}{d}\log\left[\int_{\mathbb{R}^{d}}P_{\boldsymbol{\theta}_{0}}\left(\differential{\boldsymbol{\theta}_{0}}\right)e^{\hat{\rho}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}_{0}}\right.$
(113) $\displaystyle\left.\int_{\mathbb{R}^{d\times
r}}\prod_{a=1}^{r}P_{w}\left(\differential{\boldsymbol{\theta}^{a}}\right)e^{\sum_{a=1}^{r}\left(\hat{m}^{a}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{a}+\hat{A}^{a}\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}^{a}+\hat{P}^{a}\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}^{a}+\hat{F}^{a}\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}^{a}+\hat{N}^{a}\boldsymbol{\theta}^{a}\cdot\boldsymbol{\theta}^{a}\right)+\sum_{1\leq
a\leq b\leq
r}\left(\hat{Q}^{ab}\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{b}\right)}\right]$

and the channel part of the free energy $\Psi_{y}$ as

$\Psi_{y}^{(r)}=\log\quantity[\int_{\mathbb{R}}\differential{y}\int_{\mathbb{R}}\mathrm{d}\nu
P_{0}(y\mid\nu)\int\prod_{a=1}^{r}\mathrm{~{}d}\lambda^{a}P_{g}\left(y\mid\lambda^{a},P^{a},N^{a},\varepsilon_{t}\right)\mathcal{N}\left(\nu,\lambda^{a};\mathbf{0},\Sigma^{ab}\right)]$
(114)

where we have used the fact that $(\nu_{\mu},\lambda_{\mu})\>\mu=1,\dots n$
factors over all the data points.

In the thermodynamic limit where $d\to\infty$ with $n/d$ fixed, the integral
in Eq. 111 concentrates around the values of the overlap parameters that
extremize the free entropy $\Phi^{(r)}$ and hence we can get the free energy
density as:

$\beta f_{\beta}=-\lim_{r\to
0^{+}}\frac{1}{r}\operatorname{extr}\Phi^{(r)}=-\lim_{r\to
0^{+}}\partial_{r}\operatorname{extr}\Phi^{(r)}$ (115)

###  I.2 Replica Symmetric Ansatz

We propose the following Ansatz for the variables that we have to extremise
over

$\begin{array}[]{rrr}m^{a}=m&\hat{m}^{a}=\hat{m}&\text{ for }a=1,\cdots,r\\\
q^{aa}=Q&\hat{q}^{aa}=-\frac{1}{2}\hat{Q}&\text{ for }a=1,\cdots,r\\\
q^{ab}=q&\hat{q}^{ab}=\hat{q}&\text{ for }1\leq a<b\leq r\\\
P^{a}=P&\hat{P}^{a}=-\frac{1}{2}\hat{P}&\text{ for }a=1,\cdots,r\\\
N^{a}=N&\hat{N}^{a}=-\frac{1}{2}\hat{N}&\text{ for }a=1,\cdots,r\\\
A^{a}=A&\hat{A}^{a}=-\frac{1}{2}\hat{A}&\text{ for }a=1,\cdots,r\\\
F^{a}=F&\hat{F}^{a}=\hat{F}&\text{ for }a=1,\cdots,r\\\ \end{array}$ (116)

Before we take the replica zero limit, let’s check that our ansatz above is
well-defined and does not have an order one term in $\Phi^{(r)}$ that
diverges. For this, we need to ensure that $\lim_{r\rightarrow
0^{+}}\Phi^{(r)}=0$. The trace terms depends on $r$ except for
$\rho\hat{\rho}$, $\lim_{r\rightarrow 0^{+}}\Psi_{y}^{(r)}=\lim_{r\rightarrow
0^{+}}\Psi_{w}^{(r)}=0$ hold. Note that $\lim_{r\rightarrow 0^{+}}\log
C^{r}=0$ but $\lim_{r\rightarrow 0^{+}}\partial_{r}\log C^{r}=\log C$, and
thus all we need to check is the prior part of the free energy in the zero
replica limit.

$\lim_{r\rightarrow 0^{+}}\Phi^{(r)}=-\rho\hat{\rho}$ (117)

For this limit to be zero, we must fix $\hat{\rho}=0$ and note that $\rho$ is
a constant we fixed earlier.

Plugging in the Ansatz, the trace term becomes

$-\rho\hat{\rho}-rm\hat{m}-\frac{r(r-1)}{2}q\hat{q}+\frac{r}{2}Q\hat{Q}+\frac{r}{2}A\hat{A}+\frac{r}{2}N\hat{N}+\frac{r}{2}P\hat{P}-rF\hat{F}$
(118)

Now we take the limit $r\to 0$ after dividing the trace term $T$ (which is no
longer an actual trace as we introduced overlaps beyond the traditional
replica matrix ansatz) by $r$

$\displaystyle T=$
$\displaystyle\frac{1}{2}(q\hat{q}+Q\hat{Q})+\frac{1}{2}P\hat{P}+\frac{1}{2}A\hat{A}+\frac{1}{2}N\hat{N}-m\hat{m}-F\hat{F}$
(119)

We now define $V=Q-q$ and $\hat{V}=\hat{Q}+\hat{q}$ and rewrite the trace term
as follows (by replacing $q$ and $\hat{Q}$).

$\displaystyle T$
$\displaystyle=\frac{1}{2}(q\hat{q}+(V+q)(\hat{V}-\hat{q}))+\frac{1}{2}P\hat{P}+\frac{1}{2}A\hat{A}+\frac{1}{2}N\hat{N}-m\hat{m}-F\hat{F}$
(120)
$\displaystyle=\frac{1}{2}(V\hat{V}+q\hat{V}-V\hat{q})+\frac{1}{2}P\hat{P}+\frac{1}{2}A\hat{A}+\frac{1}{2}N\hat{N}-m\hat{m}-F\hat{F}$

####  I.2.1 Prior Replica Zero Limit

Thus we can proceed plug these ansätze inside Eqs. 114 and 113 we obtain the
following for the prior term

$\displaystyle\Psi_{w}^{(r)}=\frac{1}{d}\log\left[\int_{\mathbb{R}^{d}}P_{\boldsymbol{\theta}_{0}}\left(\differential{\boldsymbol{\theta}_{0}}\right)e^{\hat{\rho}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}_{0}}\right.$
$\displaystyle\int_{\mathbb{R}^{d\times
r}}\prod_{a=1}^{r}P_{w}\quantity(\differential{\boldsymbol{\theta}^{a}})e^{\sum_{a=1}^{r}\quantity(\hat{m}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{a}\hat{F}\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}^{a})+\sum_{1\leq
a<b\leq
r}\quantity(\hat{q}\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{b})}$
(121)
$\displaystyle\left.e^{-\frac{1}{2}\sum_{a=1}^{r}\quantity(\hat{Q}\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{a}+\hat{A}\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}^{a}+\hat{P}\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}^{a}+\hat{N}\boldsymbol{\theta}^{a}\cdot\boldsymbol{\theta}^{a})}\right]$

to perform in the following the $r\to 0^{+}$ limit we can change a bit the
integral by factoring out all the terms.

To perform this simplification we will use the multidimensional Hubbard-
Stratonovic identity which reads

$e^{\frac{1}{2}\sum_{a,b=1}^{r}\boldsymbol{\theta}^{a\top}\left[\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\right]\boldsymbol{\theta}^{b}}=\mathbb{E}_{\boldsymbol{\xi}}\quantity[e^{\boldsymbol{\xi}^{\top}\,\sqrt{\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}}\sum_{a=1}^{r}\boldsymbol{\theta}^{a}}]$
(122)

where $\boldsymbol{\xi}\sim\mathcal{N}(\boldsymbol{0},1_{d})$.

Thus by calling the part inside the $\log$ in Eq. 121 with the letter
$\mathcal{A}$ we have that (putting in $\hat{\rho}=0$)

$\displaystyle\mathcal{A}=$
$\displaystyle\mathbb{E}_{\boldsymbol{\theta}_{0}}\int_{\mathbb{R}^{d\times
r}}\prod_{a=1}^{r}P_{w}\left(\differential{\boldsymbol{\theta}^{a}}\right)e^{-\sum_{a=1}^{r}\left(\hat{m}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{a}-\hat{F}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}^{a}\right)+\frac{1}{2}\sum_{1\leq
a,b\leq
r}\hat{q}\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{b}}$
(123) $\displaystyle\cdot
e^{-\frac{1}{2}\sum_{a=1}^{r}\left(\hat{V}\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{a}+\hat{A}\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}^{a}+\hat{P}\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}^{a}+\hat{N}\boldsymbol{\theta}^{a}\cdot\boldsymbol{\theta}^{a}\right)}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\boldsymbol{\theta}_{0}}\int_{\mathbb{R}^{d\times
r}}\prod_{a=1}^{r}P_{w}\left(\differential{\boldsymbol{\theta}^{a}}\right)e^{-\sum_{a=1}^{r}\quantity(\frac{\hat{V}}{2}\boldsymbol{\theta}^{a}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{a}+\frac{\hat{P}}{2}\boldsymbol{\theta}^{a}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}^{a}+\frac{\hat{A}}{2}\boldsymbol{\theta}^{a}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}^{a}+\frac{\hat{N}}{2}\boldsymbol{\theta}^{a}\boldsymbol{\theta}^{a}+\hat{m}\boldsymbol{\theta}_{0}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{a}+\hat{F}\boldsymbol{\theta}_{0}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}^{a})+\frac{1}{2}\sum_{1\leq
a,b\leq
r}\hat{q}\boldsymbol{\theta}^{a\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{b}}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\boldsymbol{\theta}_{0}}\int_{\mathbb{R}^{d\times
r}}\prod_{a=1}^{r}P_{w}\left(\differential{\boldsymbol{\theta}^{a}}\right)\mathbb{E}_{\boldsymbol{\xi}}\quantity[e^{-\frac{\hat{V}}{2}\boldsymbol{\theta}^{a}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}^{a}-\frac{\hat{P}}{2}\boldsymbol{\theta}^{a}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}^{a}-\frac{\hat{A}}{2}\boldsymbol{\theta}^{a}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}^{a}-\frac{\hat{N}}{2}\boldsymbol{\theta}^{a}\boldsymbol{\theta}^{a}-\boldsymbol{\theta}^{a}\quantity(\hat{m}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}_{0}+\hat{F}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}_{0}-\sqrt{\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}}\boldsymbol{\xi})}]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\boldsymbol{\xi},\boldsymbol{\theta}_{0}}\quantity[\quantity[\int_{\mathbb{R}^{d}}P_{w}\left(\differential{\boldsymbol{\theta}}\right)e^{-\frac{\hat{V}}{2}\boldsymbol{\theta}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}-\frac{\hat{P}}{2}\boldsymbol{\theta}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}-\frac{\hat{A}}{2}\boldsymbol{\theta}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}-\frac{\hat{N}}{2}\boldsymbol{\theta}\boldsymbol{\theta}-\boldsymbol{\theta}\quantity(\hat{m}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}_{0}+\hat{F}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}_{0}-\sqrt{\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}}\boldsymbol{\xi})}]^{r}]$

Then we can take the derivative and limit and obtain

$\Psi_{w}=\lim_{r\to
0^{+}}\partial_{r}\Psi_{w}^{(r)}=\frac{1}{d}\mathbb{E}_{\boldsymbol{\xi},\boldsymbol{\theta}_{0}}\quantity[\log\int_{\mathbb{R}^{d}}P_{w}\left(\differential{\boldsymbol{\theta}}\right)e^{-\frac{\hat{V}}{2}\boldsymbol{\theta}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}-\frac{\hat{P}}{2}\boldsymbol{\theta}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}-\frac{\hat{A}}{2}\boldsymbol{\theta}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}-\frac{\hat{N}}{2}\boldsymbol{\theta}\boldsymbol{\theta}-\boldsymbol{\theta}\quantity(\hat{m}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}_{0}+\hat{F}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}_{0}-\sqrt{\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}}\boldsymbol{\xi})}]$
(124)

where we still need to take the limit $d\to\infty$.

####  I.2.2 Channel Replica Zero Limit

Now we can focus on the channel term and rewrite it in a more suitable way for
taking the $r\to 0^{+}$ limit. In a very similar fashion as before we would
like to simplify

$\Psi_{y}^{(r)}=\log\quantity[\int_{\mathbb{R}}\differential{y}\int_{\mathbb{R}}\differential{\nu}P_{0}(y\mid\nu)\int\prod_{a=1}^{r}\differential{\lambda^{a}}P_{g}\left(y\mid\lambda^{a},P,N,\varepsilon_{t}\right)\mathcal{N}\left(\nu,\lambda^{a};\mathbf{0},\Sigma^{ab}\right)]$
(125)

We will indicate the argument of the $\log$ with $\mathcal{B}$. Additionally
we have that the martix of covariances is

$\Sigma=\begin{pmatrix}\rho&m&m&\dots&m\\\ m&Q&q&\dots&q\\\ m&q&Q&\dots&q\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\ m&q&q&\dots&Q\\\ \end{pmatrix}$ (126)

and in addition also the inverse matrix has a Replica Symmetric Structure
which is given from the following elements

$\displaystyle\quantity(\Sigma^{-1})_{00}\equiv\tilde{\rho}$
$\displaystyle=\frac{Q+(r-1)q}{\rho(Q+(r-1)q)-rm^{2}}\,,$
$\displaystyle\quantity(\Sigma^{-1})_{0a}\equiv\tilde{m}$
$\displaystyle=\frac{m}{rm^{2}-\rho(Q+(r-1)q)}\,,$ (127)
$\displaystyle\quantity(\Sigma^{-1})_{aa}\equiv\tilde{Q}$
$\displaystyle=\frac{\rho(Q+(r-2)q)-(r-1)m^{2}}{(Q-q)\left(\rho(Q+(r-1)q)-rm^{2}\right)}\,,$
$\displaystyle\quantity(\Sigma^{-1})_{ab}\equiv\tilde{q}$
$\displaystyle=\frac{m^{2}-\rho q}{(Q-q)\left(\rho(Q+(r-1)q)-rm^{2}\right)}$

and thus there is an implicit dependence on $r$ in the covariance. To check
that the inverse matrix has a RS structure as well one can think of the
formula that is used to evaluate the inverse of a matrix from the cofactors.

Also we look at the determinant of the matrix. There are three different
eigenvalue types

$\displaystyle\lambda_{1}$ $\displaystyle=Q-q\,,\quad$
$\displaystyle\lambda_{2}$
$\displaystyle=\frac{1}{2}\quantity(-Q-q(r-1)-\rho-\tilde{\Delta})\,,\quad$
$\displaystyle\lambda_{3}$
$\displaystyle=\frac{1}{2}\quantity(-Q-q(r-1)-\rho+\tilde{\Delta})\,,$ (128)
$\displaystyle d_{1}$ $\displaystyle=r-1\,,\quad$ $\displaystyle d_{2}$
$\displaystyle=1\,,\quad$ $\displaystyle d_{3}$ $\displaystyle=1\,,$

with $\tilde{\Delta}=\sqrt{4m^{2}r+(Q+q(r-1)-\rho)^{2}}$ and thus one obtains
the determinant. More explicitly we have that

$\displaystyle\det\quantity(2\pi\Sigma)$
$\displaystyle=(2\pi)^{r+1}(Q-q)^{r-1}\frac{1}{4}(-Q-q(r-1)-\rho-\tilde{\Delta})(-Q-q(r-1)-\rho+\tilde{\Delta})$
(129) $\displaystyle=(2\pi)^{r+1}(Q-q)^{r-1}(\rho(Q+(r-1)q)-rm^{2})$

Thus we have that

$\displaystyle\mathcal{B}$
$\displaystyle=\int_{\mathbb{R}}\differential{y}\int_{\mathbb{R}}\differential{\nu}P_{0}(y\mid\nu)e^{-\frac{1}{2}\tilde{\rho}\nu^{2}}\int\prod_{a=1}^{r}\differential{\lambda^{a}}P_{g}\left(y\mid\lambda^{a},P,N,\varepsilon_{t}\right)e^{-\tilde{m}\nu\sum_{a=1}^{r}\lambda^{a}-\frac{1}{2}\tilde{Q}\sum_{a=1}^{r}(\lambda^{a})^{2}-\frac{1}{2}\tilde{q}\sum_{1\leq
a,b\leq r,a\neq b}\lambda^{a}\lambda^{b}-\frac{1}{2}\log\det(2\pi\Sigma)}$
(130)
$\displaystyle=\mathbb{E}_{\xi}\int_{\mathbb{R}}\differential{y}e^{-\frac{1}{2}\log\det(2\pi\Sigma)}\int_{\mathbb{R}}\differential{\nu}P_{0}(y\mid\nu)e^{-\frac{1}{2}\tilde{\rho}\nu^{2}}\quantity[\int\differential{\lambda}P_{g}\left(y\mid\lambda,P,N,\varepsilon_{t}\right)e^{-\frac{\tilde{Q}-\tilde{q}}{2}\lambda^{2}+(\sqrt{-\tilde{q}}\xi-\tilde{m}\nu)\lambda}]^{r}$

Now we can follow a similar procedure as before and define $V=Q-q$ we have
that and the limit is

$\displaystyle\Psi_{y}$ $\displaystyle=\lim_{r\to
0^{+}}\partial_{r}\Psi_{y}^{(r)}=\mathbb{E}_{\xi}\quantity[\int_{\mathbb{R}}\differential{y}\int\frac{\differential{\nu}}{\sqrt{2\pi\rho}}P_{0}(y\mid\nu)e^{-\frac{1}{2\rho}\nu^{2}}\log\quantity[\int\frac{\differential{\lambda}}{\sqrt{2\pi}}P_{g}(y\mid\lambda,P,N,\varepsilon_{t})e^{-\frac{1}{2}\frac{\lambda^{2}}{V}+\left(\frac{\sqrt{q-m^{2}/\rho}}{V}\xi+\frac{m/\rho}{V}\nu\right)\lambda}]]$
(131) $\displaystyle-\frac{1}{2}\log V-\frac{1}{2}\frac{q}{V}$

We would like to rewrite the quantities with the help of the following
definition

$\mathcal{Z}_{0}(y,\omega,V)=\int\frac{\mathrm{d}x}{\sqrt{2\pi
V}}e^{-\frac{1}{2V}(x-\omega)^{2}}P_{0}(y\mid
x)\,,\quad\mathcal{Z}_{y}(y,\omega,V,P,n)=\int\frac{\mathrm{d}x}{\sqrt{2\pi
V}}e^{-\frac{1}{2V}(x-\omega)^{2}}P_{g}(y\mid x,P,N,\varepsilon_{t})\,,$ (132)

The result becomes thus

$\mathbb{E}_{\xi}\quantity[\int_{\mathbb{R}}\differential{y}\mathcal{Z}_{0}\quantity(y,\frac{m}{\sqrt{q}}\xi,\rho-\frac{m^{2}}{q})\log\mathcal{Z}_{y}(y,\sqrt{q}\xi,V,P,N)]$
(133)

Now there are two things that we still need to do : find the form for the
prior term and take the limit $\beta\to\infty$.

###  I.3 Prior term for $\ell_{2}$ regularisation

To be as general as possible we would like to include the case of a possible
non isotropic regularisation. Thus

$P_{w}(\differential{\boldsymbol{\theta}})=\frac{1}{(2\pi)^{d/2}}\exp(-\frac{\beta\lambda}{2}\boldsymbol{\theta}\boldsymbol{\Sigma}_{\boldsymbol{w}}\boldsymbol{\theta})\differential{\boldsymbol{\theta}}$
(134)

We want to calculate the term inside the $\log$ in Eq. 124

$\displaystyle\int_{\mathbb{R}^{d}}P_{w}(\differential{\boldsymbol{\theta}})e^{-\frac{\hat{V}}{2}\boldsymbol{\theta}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}-\frac{\hat{P}}{2}\boldsymbol{\theta}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}-\frac{\hat{A}}{2}\boldsymbol{\theta}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}-\frac{\hat{N}}{2}\boldsymbol{\theta}\boldsymbol{\theta}-\boldsymbol{\theta}\quantity(\hat{m}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\theta}_{0}+\hat{F}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\theta}_{0}-\sqrt{\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}}\boldsymbol{\xi})}$
(135)
$\displaystyle=\frac{\exp\left(\frac{1}{2}\left(-\hat{m}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}-\hat{F}(\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}})^{\top}\boldsymbol{\theta}_{0}+\sqrt{\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}}\boldsymbol{\xi}\right)^{\top}\boldsymbol{\Lambda}^{-1}\left(-\hat{m}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}-\hat{F}(\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}})^{\top}\boldsymbol{\theta}_{0}+\sqrt{\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}}\boldsymbol{\xi}\right)^{\top}\right)}{\sqrt{\det\boldsymbol{\Lambda}}}$

where we defined
$\boldsymbol{\Lambda}=\beta\lambda\boldsymbol{\Sigma}_{\boldsymbol{w}}+\hat{V}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{P}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}+\hat{A}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}+\hat{N}\mathbb{1}$.
Now the prior term becomes after taking the log and using the identity
$\log\operatorname{det}=\tr\log$

$\displaystyle\Psi_{w}$
$\displaystyle=\frac{1}{d}\mathbb{E}_{\boldsymbol{\xi},\boldsymbol{\theta}_{0}}\quantity[\frac{1}{2}\left(-\hat{m}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}-\hat{F}(\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}})^{\top}\boldsymbol{\theta}_{0}+\sqrt{\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}}\boldsymbol{\xi}\right)^{\top}\boldsymbol{\Lambda}^{-1}\left(-\hat{m}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}-\hat{F}(\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}})^{\top}\boldsymbol{\theta}_{0}+\sqrt{\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}}\boldsymbol{\xi}\right)^{\top}]-\frac{1}{2d}\tr\log\boldsymbol{\Lambda}$
(136)
$\displaystyle=\frac{1}{2d}\tr\quantity[\left(\hat{m}^{2}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{m}\hat{F}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}+\hat{m}\hat{F}(\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}})^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{F}^{2}(\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}})^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}+\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\right)\boldsymbol{\Lambda}^{-1}]-\frac{1}{2d}\tr\log\boldsymbol{\Lambda}$

The factor $\frac{1}{d}$ comes from the required scaling on $d$ for the free
entropy and the expectation from our replica zero limit of the prior term.

###  I.4 Zero temperature limit

We now need to take the zero temperature limit for this case. The explicit
scalings of the parameters are

$\begin{array}[]{rrrrrrr}V\rightarrow\beta^{-1}V&q\rightarrow q&m\rightarrow
m&A\rightarrow A&N\rightarrow N&P\rightarrow P&F\rightarrow F\\\
\hat{V}\rightarrow\beta\hat{V}&\hat{q}\rightarrow\beta^{2}\hat{q}&\hat{m}\rightarrow\beta\hat{m}&\hat{A}\rightarrow\beta\hat{A}&\hat{N}\rightarrow\beta\hat{N}&\hat{P}\rightarrow\beta\hat{P}&\hat{F}\rightarrow\beta\hat{F}\\\
\end{array}$ (137)

The limit of the prior term is

$\Psi_{w}=\lim_{\beta\to\infty}\frac{1}{\beta}\Psi_{w}=\frac{1}{2d}\tr\quantity[\left(\hat{m}^{2}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{m}\hat{F}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}+\hat{m}\hat{F}(\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}})^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{F}^{2}(\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}})^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}+\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\right)\boldsymbol{\Lambda}^{-1}]$
(138)

To understand the limit of the channel term, we need to get the following
insight for the limit of the channel partition function

$\displaystyle\mathcal{Z}_{y}(y,\omega,V,P,N)$
$\displaystyle=\int\frac{\mathrm{d}x}{\sqrt{2\pi
V}}e^{-\frac{\beta}{2V}(x-\omega)^{2}}P_{g}(y\mid x,P,N,\varepsilon_{t})$
(139) $\displaystyle=\sqrt{\beta}\int\frac{\mathrm{d}x}{\sqrt{2\pi
V}}e^{-\frac{\beta}{2V}(x-\omega)^{2}}\frac{1}{\sqrt{2\pi}}e^{(-\beta
g(yx-\frac{\varepsilon_{t}}{\sqrt{d}}\frac{P}{\sqrt{N}}))}$
$\displaystyle\underset{\beta\rightarrow\infty}{=}e^{-\beta\mathcal{M}_{Vg(y,\cdot)}(\omega)}$

where we introduced the Moreau envelope defined in Eq. 25.

Then the limit of the channel term becomes

$\Psi_{y}=\lim_{\beta\to\infty}\frac{1}{\beta}\Psi_{y}=-\mathbb{E}_{\xi}\left[\int\differential{y}\mathcal{Z}_{0}\left(y,\frac{m}{\sqrt{q}}\xi,\rho-\frac{m^{2}}{q}\right)\mathcal{M}_{Vg(y,\cdot;P,N,\varepsilon_{t})}(\sqrt{q}\xi)\right]$
(140)

where $\mathcal{M}_{Vg(y,\cdot;A,N,\varepsilon_{t})}$ is the Moreau envelope
of the modified loss function defined in Eq. 7 with the relevant quantities
changed for their overlaps and $\xi\sim\mathcal{N}(0,1)$.

After taking the zero temperature limit, we are left with the following
expression for the free energy density

$\displaystyle\lim_{\beta\rightarrow\infty}f_{\beta}$
$\displaystyle=\mathop{\operatorname{extr}}_{\begin{subarray}{c}V,q,m,A,N,P,F,\\\
\hat{V},\hat{q},\hat{m},\hat{A},\hat{N},\hat{P},\hat{F}\end{subarray}}\left\\{-\frac{1}{2}(q\hat{V}-\hat{q}V)-\frac{1}{2}P\hat{P}-\frac{1}{2}A\hat{A}-\frac{1}{2}N\hat{N}+m\hat{m}+F\hat{F}+\alpha\mathbb{E}_{\xi}\left[\int\differential{y}\mathcal{Z}_{0}\mathcal{M}_{Vg(y,\cdot)}\right]\right.$
(141)
$\displaystyle\left.-\frac{1}{2d}\tr\left(\hat{m}^{2}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{m}\hat{F}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}+\hat{m}\hat{F}(\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}})^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{F}^{2}(\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}})^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}+\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\right)\boldsymbol{\Lambda}^{-1}\right\\}$

###  I.5 Saddle-point equations

The extremisation condition of Eq. 141 can be translated into the overlap
needing to satisfy the following

$\displaystyle\hat{V}$ $\displaystyle=2\alpha\partial_{q}\Psi_{y},$
$\displaystyle q=-2\partial_{\hat{V}}\Psi_{w}$ (142) $\displaystyle\hat{q}$
$\displaystyle=-2\alpha\partial_{V}\Psi_{y},$ $\displaystyle
V=2\partial_{\hat{q}}\Psi_{w},$ $\displaystyle\hat{N}$
$\displaystyle=2\alpha\partial_{N}\Psi_{y},$ $\displaystyle
N=-2\partial_{\hat{N}}\Psi_{w}$ $\displaystyle\hat{P}$
$\displaystyle=2\alpha\partial_{P}\Psi_{y},$ $\displaystyle
P=-2\partial_{\hat{P}}\Psi_{w}$ $\displaystyle\hat{A}$
$\displaystyle=2\alpha\partial_{A}\Psi_{y},$ $\displaystyle
A=-2\partial_{\hat{A}}\Psi_{w}$ $\displaystyle\hat{m}$
$\displaystyle=-\alpha\partial_{m}\Psi_{y},$ $\displaystyle
m=\partial_{\hat{m}}\Psi_{w}$ $\displaystyle\hat{F}$
$\displaystyle=-\alpha\partial_{F}\Psi_{y},$ $\displaystyle
F=\partial_{\hat{F}}\Psi_{w}.$

As we pre-announced we would like to find the stationary values that dominate
the integral and to do so we should derive the exponent with respect to all
the order parameters. The saddle points that depend on $m,q,V,\hat{m},\hat{q}$
and $\hat{V}$ are of a similar form as those found already in [32]. We need
thus to derive with respect to $A,N,P,F,\hat{A},\hat{P},\hat{F}$ and
$\hat{N}$.

####  I.5.1 The Channel Saddle-Point Equations

Let us begin by looking at the derivatives with respect to $P$ and $N$. These
derivatives amount to computing the derivative of the Moreau-envelope with
respect to $P$ and $N$ since we have that

$\displaystyle\partial_{P}\Psi_{y}$
$\displaystyle=\mathbb{E}_{y,\xi}\left[\mathcal{Z}_{0}\quantity(y,\frac{m}{\sqrt{q}}\xi,\rho-\frac{m^{2}}{q})\partial_{P}\mathcal{M}_{Vg(y,\cdot;P,N,\varepsilon_{t})}(\sqrt{q}\xi)\right]\,,$
(143) $\displaystyle\partial_{N}\Psi_{y}$
$\displaystyle=\mathbb{E}_{y,\xi}\left[\mathcal{Z}_{0}\quantity(y,\frac{m}{\sqrt{q}}\xi,\rho-\frac{m^{2}}{q})\partial_{N}\mathcal{M}_{Vg(y,\cdot;P,N,\varepsilon_{t})}(\sqrt{q}\xi)\right]\,,$

In this specific case we have that

$\mathcal{M}_{Vg(y,\cdot,P,N,\varepsilon_{t})}(\omega)=\inf_{x\in\mathbb{R}}\left[\frac{(x-\omega)^{2}}{2V}+\ell\quantity(yx-\varepsilon_{t}\frac{P}{\sqrt{N}})\right]=\mathcal{M}_{V\ell(y,\cdot)}\quantity(\omega-y\frac{\varepsilon_{t}\,P}{\sqrt{N}})$
(144)

where we remind that this specific form is possible since
$y\in\quantity{+1,-1}$. With this we can relate these function the the
derivative of the Moreau envelope with respect to its input as

$\displaystyle\partial_{P}\mathcal{M}_{Vg(y,\cdot,P,N,\varepsilon_{t})}(\omega)$
$\displaystyle=-y\frac{\varepsilon_{t}}{\sqrt{N}}\mathcal{M}^{\prime}_{V\ell(y,\cdot)}\quantity(\omega-y\frac{\varepsilon_{t}P}{\sqrt{N}})$
(145)
$\displaystyle\partial_{N}\mathcal{M}_{Vg(y,\cdot,P,N,\varepsilon_{t})}(\omega)$
$\displaystyle=y\frac{\varepsilon_{t}P}{\sqrt{N^{3}}}\mathcal{M}^{\prime}_{V\ell(y,\cdot)}\quantity(\omega-y\frac{\varepsilon_{t}P}{\sqrt{N}})$

With this, we can write the new equations as

$\hat{P}=\alpha\frac{\varepsilon_{t}}{\sqrt{N}}\mathbb{E}_{\xi}\left[\int_{\mathbb{R}}\differential{y}y\mathcal{Z}_{0}f_{g}\right]\,,\quad\hat{N}=-\alpha\frac{\varepsilon_{t}\,P}{\sqrt{N^{3}}}\mathbb{E}_{\xi}\left[\int_{\mathbb{R}}\differential{y}y\mathcal{Z}_{0}f_{g}\right]\,,$
(146)

where we have defined

$f_{g}(y,\omega,V,P,N,\varepsilon_{t})=-\mathcal{M}^{\prime}_{Vg(y,\cdot,P,N,\varepsilon_{t})}(\omega)$
(147)

As the channel term does not depend on the overlaps $F$ and $A$ the hat
equations are trivially zero and

$\hat{A}=0\,,\quad\hat{F}=0\,.$ (148)

The remaining three equations can be found as in [32], where the only
difference lies in the dependence of $f_{g}$ on $\hat{P}$ and $\hat{N}$. Note
that here and above we denote by $z^{*}$ the value of the proximal at any
given point of integration.

$\hat{V}=-\alpha\mathbb{E}_{\xi}\left[\int_{\mathbb{R}}\differential{y}\mathcal{Z}_{0}\partial_{\omega}f_{g}\right]\,,\quad\hat{q}=\alpha\mathbb{E}_{\xi}\left[\int_{\mathbb{R}}\differential{y}\mathcal{Z}_{0}f_{g}^{2}\right]\,,\quad\hat{m}=\alpha\mathbb{E}_{\xi}\left[\int_{\mathbb{R}}\differential{y}\partial_{\omega}\mathcal{Z}_{0}f_{g}\right]\,,$
(149)

####  I.5.2 The Prior Saddle-Point Equations

For the prior saddle-point equations our starting point is

$\displaystyle\Psi_{w}$
$\displaystyle=\frac{1}{2d}\tr\quantity[\left(\hat{m}^{2}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{m}\hat{F}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}+\hat{m}\hat{F}(\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}})^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{F}^{2}(\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}})^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}+\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\right)\boldsymbol{\Lambda}^{-1}]$
(150)

where for simplicity of notation we define
$\boldsymbol{\Lambda}=\beta\lambda\boldsymbol{\Sigma}_{\boldsymbol{w}}+\hat{V}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{P}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}+\hat{A}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}+\hat{N}\mathbb{1}$
and we will use
$\boldsymbol{H}=\hat{m}^{2}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{m}\hat{F}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}+\hat{m}\hat{F}(\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}})^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{F}^{2}(\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}})^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}+\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}$.

As the channel equations for $\hat{A}$ and $\hat{F}$ are trivially zero, we
want to start with these derivatives as the following expressions will
simplify considerably.

$\displaystyle A$
$\displaystyle=\partial_{\hat{A}}\Psi_{w}=\frac{1}{d}\tr\quantity[\boldsymbol{H}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\Lambda}^{-2}]$
(151) $\displaystyle F$
$\displaystyle=\partial_{\hat{F}}\Psi_{w}=\frac{1}{d}\tr\quantity[\hat{m}\left(\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\right)\boldsymbol{\Lambda}^{-1}]$

We want to compute a few derivatives of the term $\Psi_{w}$ to obtain
equations for the overlaps $P,N$.

We begin with the hat-variables $\hat{P}$ and $\hat{N}$

$P=\partial_{\hat{P}}\Psi_{w}=\frac{1}{d}\tr\quantity[\boldsymbol{H}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\Lambda}^{-2}]$
(152)

and

$N=\partial_{\hat{N}}\Psi_{w}=\frac{1}{d}\tr\quantity[\boldsymbol{H}\mathbb{1}\boldsymbol{\Lambda}^{-2}]$
(153)

As before the derivative w.r.t. $\hat{q},\hat{a},\hat{n}$ follow from previous
literature as

$V=\partial_{\hat{q}}\Psi_{w}=\frac{1}{d}\tr\quantity[\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\Lambda}^{-1}]\,,\quad
q=\partial_{\hat{V}}\Psi_{w}=\frac{1}{d}\tr\quantity[\boldsymbol{H}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\Lambda}^{-2}]\,,\quad
m=\partial_{\hat{m}}\Psi_{w}=\frac{1}{d}\tr\quantity[\hat{m}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\Lambda}^{-1}]\,.$
(154)

Note that for the numerical evaluation $\frac{1}{d}\tr$ is just the mean of
the eigenspectrum.

###  I.6 Final set of saddle point equations for $\ell_{2}$ regularisation

We state here our final set of saddle point equations for reference

$\displaystyle\begin{cases}\hat{m}=\alpha\mathbb{E}_{\xi}\left[\int_{\mathbb{R}}\differential{y}\partial_{\omega}\mathcal{Z}_{0}f_{g}(\sqrt{q}\xi,P,N,\varepsilon_{t})\right]\\\
\hat{q}=\alpha\mathbb{E}_{\xi}\left[\int_{\mathbb{R}}\differential{y}\mathcal{Z}_{0}f_{g}^{2}(\sqrt{q}\xi,P,N,\varepsilon_{t})\right]\\\
\hat{V}=-\alpha\mathbb{E}_{\xi}\left[\int_{\mathbb{R}}\differential{y}\mathcal{Z}_{0}\partial_{\omega}f_{g}(\sqrt{q}\xi,P,N,\varepsilon_{t})\right]\\\
\hat{P}=\varepsilon_{t}\frac{1}{\sqrt{N}}\alpha\mathbb{E}_{\xi}\quantity[\int_{\mathbb{R}}\differential{y}\mathcal{Z}_{0}yf_{g}(\sqrt{q}\xi,P,N,\varepsilon_{t})]\\\
\hat{N}=-\varepsilon_{t}\frac{P}{N^{3/2}}\alpha\mathbb{E}_{\xi}\quantity[\int_{\mathbb{R}}\differential{y}\mathcal{Z}_{0}yf_{g}(\sqrt{q}\xi,P,N,\varepsilon_{t})]\end{cases}$
(155)
$\displaystyle\begin{cases}m=\frac{1}{d}\tr\quantity[\hat{m}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\Lambda}^{-1}]\\\
q=\frac{1}{d}\tr\quantity[\quantity(\hat{m}^{2}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}})\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\Lambda}^{-2}]\\\
V=\frac{1}{d}\tr\quantity[\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}\boldsymbol{\Lambda}^{-1}]\\\
P=\frac{1}{d}\tr\quantity[\quantity(\hat{m}^{2}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}})\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\Lambda}^{-2}]\\\
N=\frac{1}{d}\tr\quantity[\quantity(\hat{m}^{2}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}})\boldsymbol{\Lambda}^{-2}]\end{cases}$

where we remember the definitions of
$\boldsymbol{\Lambda}=\lambda\boldsymbol{\Sigma}_{\boldsymbol{w}}+\hat{V}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{P}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}+\hat{N}\mathbb{1}$.
Note $\hat{A},\hat{F}$ are exactly zero because the channel is not dependent
on $A,F$. The values for $A$ and $F$ becomes

$A=\frac{1}{d}\tr\quantity[\quantity(\hat{m}^{2}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}+\hat{q}\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}})\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\boldsymbol{\Lambda}^{-2}]\,,\quad
F=\frac{1}{d}\tr\quantity[\hat{m}\left(\boldsymbol{\Sigma}_{\boldsymbol{\boldsymbol{x}}}^{\top}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\upsilon}}\right)\boldsymbol{\Lambda}^{-1}]\,.$
(156)

where we remind that $f_{g}$ is defined in Eq. 25.

##  Appendix J Numerical Empirical Risk Minimisation

Computing the loss, gradient and hessian of the logistic loss, let alone the
adversarially perturbed logistic loss pose well-known numerical challenges
[51]. In this section we explain in detail how to numerically evaluate the
ERM-estimator. To begin, here is the loss function over a data set
${(\boldsymbol{x}_{\mu},y_{\mu})}\>\mu=1,...,N$.

$\mathcal{L}=\sum_{\mu=1}^{N}\log\left(1+\exp\quantity(-y_{\mu}\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}+\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}})\right)$
(157)

To compute the gradient, we switch from a -1/1 labelling to a 0/1 labelling.
We can write the following:

$\displaystyle\mathcal{L}$
$\displaystyle=\sum_{\mu=1}^{N}\log\left(1+\exp\quantity(-y_{\mu}\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}+\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}})\right)$
(158)
$\displaystyle=\sum_{\mu=1}^{N}\delta_{y_{\mu},+1}\log\left(1+\exp\quantity(-\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}+\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}})\right)+\delta_{y_{\mu},-1}\log\left(1+\exp\quantity(\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}+\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}})\right)$
$\displaystyle=\sum_{\mu=1}^{N}y_{\mu}\log\left(1+\exp\quantity(-\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}+\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}})\right)+(1-y_{\mu})\log\left(1+\exp\quantity(\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}+\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}})\right)$
$\displaystyle=\sum_{\mu=1}^{N}-y_{\mu}\log\left(\frac{\exp\quantity(\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}-\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}})}{1+\exp\quantity(\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}-\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}})}\right)+(1-y_{\mu})\log\left(1+\exp\quantity(\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}+\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}})\right)$
$\displaystyle=\sum_{\mu=1}^{N}-y_{\mu}\quantity(\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}-\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}})+y_{\mu}\log\left(1+\exp\quantity(\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}-\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}})\right)$
$\displaystyle\qquad+(1-y_{\mu})\log\left(1+\exp\quantity(\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}+\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}})\right)$

In the second step we are free to choose the labelling to be 0/1 as there is
no more explicit dependence on the $y_{\mu}$ instance.

###  J.1 Computing the Loss

To accurately compute the loss we base our computation on the work of [41] who
showed how to accurately compute $\log(1+\exp(x))$. The idea is to choose the
most suitable approximation depending on the argument of the function.

Finally, we can write the loss with 0/1 labels as

$\mathcal{L}=\sum_{\mu=1}^{N}-y_{\mu}\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}+y_{\mu}\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}}+(1-y_{\mu})\log\left(1+\exp\left(\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}+\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}}\right)\right)+y_{\mu}\log\left(1+\exp\left(\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}-\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}}\right)\right)$
(159)

In [41] the authors showed how to accurately compute $\log(1+\exp(x))$,
additionally [50] inspire us to achieve faster computation by extending the
case distinction by [41] by the $\leq 2$ and using directly the log function.

$\log(1+\exp(x)):=\begin{cases}\exp(x)&x\leq-37\\\
\operatorname{log1p}(\exp(x))&-37<x\leq x_{0}:=-2\\\
\log(1+\exp(x))&x_{0}<x\leq x_{1}:=18\\\ x+\exp(-x)&x_{1}<x\leq x_{2}:=33.3\\\
x&x>x_{2}\end{cases}$ (160)

###  J.2 Computing the Gradient

First we notice that

$\derivative{z}\quantity[\log\quantity(1+e^{z})]=\frac{1}{1+e^{-z}}\,,\quad\derivative{z}\quantity[\frac{1}{1+e^{-z}}]=\frac{e^{z}}{(1+e^{z})^{2}}\,,\quad\derivative{z}\quantity[\frac{1}{1+e^{z}}]=-\frac{1}{2\cosh(z)+2}\,.$
(161)

also we remember the following derivatives

$\gradient_{\boldsymbol{\theta}}\quantity[\boldsymbol{\theta}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}]=2\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}\,,\quad\operatorname{Hess}\quantity(\boldsymbol{\theta}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta})=2\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\,,\quad\gradient_{\boldsymbol{\theta}}\quantity[\norm{\boldsymbol{\theta}}_{2}]=\frac{\boldsymbol{\theta}}{\norm{\boldsymbol{\theta}}_{2}}\,,\quad\operatorname{Hess}\quantity(\norm{\boldsymbol{\theta}}_{2})=\frac{\mathbf{1}}{\norm{\boldsymbol{\theta}}_{2}}-\frac{\boldsymbol{\theta}\boldsymbol{\theta}^{\top}}{\norm{\boldsymbol{\theta}}_{2}^{3}}\,,$
(162)

since we suppose that $\boldsymbol{\Sigma}_{\boldsymbol{\delta}}$ is a
symmetric matrix.

Call the derivative of the optimal attack $\boldsymbol{h}$

$\displaystyle\boldsymbol{h}=\gradient_{\boldsymbol{\theta}}B=\gradient_{\boldsymbol{\theta}}\quantity[\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}}]$
$\displaystyle=\frac{2\varepsilon_{t}}{\sqrt{d}\norm{\boldsymbol{\theta}}_{2}}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}-\frac{\varepsilon_{t}}{\sqrt{d}\norm{\boldsymbol{\theta}}_{2}^{3}}\boldsymbol{\theta}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta}\>\boldsymbol{\theta}$
(163)

We define the arguments of the sigmoid activation $C$ and $\overline{C}$

$C_{\mu}:=\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}+\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}}\,,\qquad\overline{C}_{\mu}=\frac{\boldsymbol{x}_{\mu}^{\top}\boldsymbol{w}}{\sqrt{d}}-\frac{\varepsilon_{t}\boldsymbol{w}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{w}}{\sqrt{d}\sqrt{\boldsymbol{w}^{\top}\boldsymbol{w}}}\,,$
(164)

Then we can write the derivative of the loss as

$\displaystyle\gradient_{\boldsymbol{\theta}}\mathcal{L}$
$\displaystyle=\sum_{\mu=1}^{N}\left[-y_{\mu}\frac{\boldsymbol{x}_{\mu}}{\sqrt{d}}+y_{\mu}\boldsymbol{h}+\frac{(1-y_{\mu})}{1+\exp\left(-C_{\mu}\right)}\quantity(\frac{\boldsymbol{x}_{\mu}}{\sqrt{d}}+\boldsymbol{h})+\frac{y_{\mu}}{1+\exp\left(-\overline{C}_{\mu}\right)}\quantity(\frac{\boldsymbol{x}_{\mu}}{\sqrt{d}}-\boldsymbol{h})\right]$
(165)
$\displaystyle=\sum_{\mu=1}^{N}\left[\boldsymbol{h}\left(\frac{(1-y_{\mu})}{1+\exp\left(-C_{\mu}\right)}+\frac{y_{\mu}}{1+\exp\quantity(\overline{C}_{\mu})}\right)+\frac{\boldsymbol{x}_{\mu}}{\sqrt{d}}\left(\frac{(1-y_{\mu})}{1+\exp\left(-C_{\mu}\right)}-\frac{y_{\mu}}{1+\exp\quantity(\overline{C}_{\mu})}\right)\right]$

With this it is easy to see what the loss per sample is when factoring out the
data $\boldsymbol{x}$ and the derivative of the optimal attack
$\boldsymbol{h}$. We see that their respective contributions are similar. One
can avoid overflows by placing the exponential part of the sigmoid carefully
depending on the argument.

###  J.3 Computing the Hessian

To produce this derivative it is easier to derive in index notation with
respect to $\boldsymbol{\theta}_{i}$. We start from

$\displaystyle\partialderivative{B}{\boldsymbol{\theta}_{j}}{\boldsymbol{\theta}_{i}}$
$\displaystyle=\partialderivative{\boldsymbol{h}_{i}}{\boldsymbol{\theta}_{j}}=\frac{\varepsilon_{t}}{\sqrt{d}}\quantity[\frac{2\quantity(\boldsymbol{\Sigma}_{\boldsymbol{\delta}})_{ij}}{\norm{\boldsymbol{\theta}}_{2}}-\frac{2}{\norm{\boldsymbol{\theta}}_{2}^{3}}\quantity[2\quantity(\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta})_{i}\boldsymbol{\theta}_{j}+\frac{1}{2}\quantity(\boldsymbol{\theta}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta})\,\delta_{ij}]+\frac{3}{\norm{\boldsymbol{\theta}}_{2}^{5}}\quantity(\boldsymbol{\theta}^{\top}\boldsymbol{\Sigma}_{\boldsymbol{\delta}}\boldsymbol{\theta})\,\boldsymbol{\theta}_{i}\boldsymbol{\theta}_{j}]$
(166)

where we assumed $\boldsymbol{\Sigma}_{\boldsymbol{\delta}}$ to be symmetric.

Now we can continue by considering the second term in Eq. 165 and
differentiating it

$\frac{\boldsymbol{x}^{\mu}_{i}}{\sqrt{d}}\quantity[\frac{1-y_{\mu}}{2(\cosh(C_{\mu})+1)}\quantity(\frac{\boldsymbol{x}^{\mu}_{j}}{\sqrt{d}}+\boldsymbol{h}_{j})+\frac{y_{\mu}}{2(\cosh(\overline{C_{\mu}})+1)}\quantity(\frac{\boldsymbol{x}^{\mu}_{j}}{\sqrt{d}}-\boldsymbol{h}_{j})]$
(167)

then we get from the derivative of the first piece a term similar to the
previous one with $\boldsymbol{h}_{i}$ instead of
$\boldsymbol{x}^{\mu}_{i}/\sqrt{d}$ which is

$\boldsymbol{h}_{i}\quantity[\frac{1-y_{\mu}}{2(\cosh(C_{\mu})+1)}\quantity(\frac{\boldsymbol{x}^{\mu}_{j}}{\sqrt{d}}+\boldsymbol{h}_{j})-\frac{y_{\mu}}{2(\cosh(\overline{C_{\mu}})+1)}\quantity(\frac{\boldsymbol{x}^{\mu}_{j}}{\sqrt{d}}-\boldsymbol{h}_{j})]$
(168)

and the last term

$\partialderivative{B}{\boldsymbol{\theta}_{j}}{\boldsymbol{\theta}_{i}}\quantity(\frac{(1-y_{\mu})}{1+\exp\left(-C_{\mu}\right)}+\frac{y_{\mu}}{1+\exp\quantity(\overline{C}_{\mu})})$
(169)

where the value in front is given by Eq. 166 and to complete it one should sum
over all the data points to obtain the hessian. Reorganising the terms a bit
gives us the hessian as

$\displaystyle\partial_{\boldsymbol{\theta}^{2}}\mathcal{L}$
$\displaystyle=\sum_{\mu=1}^{N}\partialderivative{B}{\boldsymbol{\theta}_{j}}{\boldsymbol{\theta}_{i}}\quantity(\frac{(1-y_{\mu})}{1+\exp\left(-C_{\mu}\right)}+\frac{y_{\mu}}{1+\exp\quantity(\overline{C}_{\mu})})$
(170)
$\displaystyle+\boldsymbol{h}_{i}\quantity[\frac{1-y_{\mu}}{2(\cosh(C_{\mu})-1)}\quantity(\frac{\boldsymbol{x}^{\mu}_{j}}{\sqrt{d}}+\boldsymbol{h}_{j})-\frac{y_{\mu}}{2(\cosh(\overline{C_{\mu}})+1)}\quantity(\frac{\boldsymbol{x}^{\mu}_{j}}{\sqrt{d}}-\boldsymbol{h}_{j})]$
$\displaystyle+\frac{\boldsymbol{x}^{\mu}_{i}}{\sqrt{d}}\quantity[\frac{1-y_{\mu}}{2(\cosh(C_{\mu})-1)}\quantity(\frac{\boldsymbol{x}^{\mu}_{j}}{\sqrt{d}}+\boldsymbol{h}_{j})+\frac{y_{\mu}}{2(\cosh(\overline{C_{\mu}})+1)}\quantity(\frac{\boldsymbol{x}^{\mu}_{j}}{\sqrt{d}}-\boldsymbol{h}_{j})]$