---
abstract: 'We study the regularisation induced in neural networks by Gaussian noise injections (GNIs). Though such injections have been extensively studied when applied to data, there have been few studies on understanding the regularising effect they induce when applied to network activations. Here we derive the explicit regulariser of GNIs, obtained by marginalising out the injected noise, and show that it is a form of Tikhonov regularisation which penalises functions with high-frequency components in the Fourier domain. We show analytically and empirically that such regularisation produces calibrated classifiers with large classification margins and that the explicit regulariser we derive is able to reproduce these effects.'
author:
- |
    \
    \
    \
    \
    \
    ${}^{1}$Department of Statistics, University of Oxford\
    ${}^{2}$Alan Turing Institute, London\
    ${}^{3}$LCTI, Télécom Paris, Institute Polytechnique de Paris\
    ${}^{4}$Oxford-Man Institute, University of Oxford
bibliography:
- 'references.bib'
title: Explicit Regularisation in Gaussian Noise Injections
---

Introduction {#intro}
============

Noise injections are a family of methods that involve adding or multiplying samples from a noise distribution, typically an isotropic Gaussian, to the weights and activations of a neural network during training. The benefits of such methods are well documented. Models trained with noise often generalise better to unseen data and are less prone to overfitting [@Srivastava2014; @Kingma; @Poole2014].

Even though the regularisation conferred by Gaussian noise injections (GNIs) can be observed empirically, and the benefits of noising data are well understood theoretically [@Bishop1995; @Cohen2019; @Webb1994], there have been few studies on understanding the benefits of methods that inject noise *throughout* a network. Here we study the *explicit* regularisation of such injections, which is the added positive term we obtain to the loss function when we marginalise out the noise we have injected. In particular, we show that the regularisers we derive are a form of Tikhonov Regularisation which penalises some function space norm of a neural network [@Girosi1990; @Burger2003; @Bishop1995]. Here we show that noise injections induce a penalisation of a Sobolev space norm and in particular we show that this corresponds to a prior that penalises functions with high-frequency components in the Fourier domain.

Here we show analytically that such regularisation promotes larger classification margins which make models robust to noise, explaining the empirical results of @Rahaman2019 which established that neural networks with lower-frequency spectra are more robust to input perturbations. Finally we show empirically that the explicit regulariser also improves model calibration.

Concretely our contributions are:

-   We derive an analytic form for the explicit regulariser induced by GNIs.

-   We establish a connection between this explicit regulariser and Tikhonov regularisation, showing that it penalises neural networks with high-frequency content in the Fourier domain.

-   Finally, we show analytically and empirically that this regularisation induces larger classification margins and better calibration of models.

Background
==========

Gaussian Noise Injections {#sec:sgd_under_injections}
-------------------------

Training a neural network involves optimising network parameters to maximise the marginal likelihood of a set of labels ${{\boldsymbol{\mathbf{}}}}{y}$ given features ${{\boldsymbol{\mathbf{}}}}{x}$ via gradient descent. With a training dataset composed of $N$ data-label pairs $\mathcal{D} = \{({{\boldsymbol{\mathbf{}}}}{x}_0, {{\boldsymbol{\mathbf{}}}}{y}_0) ... ({{\boldsymbol{\mathbf{}}}}{x}_{N-1}, {{\boldsymbol{\mathbf{}}}}{y}_{N-1})\}$ and a feed-forward neural network with $D$ parameters divided into $L$ layers: ${{\boldsymbol{\mathbf{}}}}{\theta} = \{{{\boldsymbol{\mathbf{}}}}{W}_1,...,{{\boldsymbol{\mathbf{}}}}{W}_L\}$, ${{\boldsymbol{\mathbf{}}}}{\theta} \in \mathbb{R}^D$, our objective is to minimise the expected negative log likelihood and find the optimal set of parameters ${{\boldsymbol{\mathbf{}}}}{\theta}^*$ satisfying: $$\begin{aligned}
    {{\boldsymbol{\mathbf{}}}}{\theta}_* = \operatorname*{arg\,min}_{{{\boldsymbol{\mathbf{}}}}{w}} \mathcal{L}( \mathcal{D}; {{\boldsymbol{\mathbf{}}}}{w}),
    \qquad
    \mathcal{L}( \mathcal{D}; {{\boldsymbol{\mathbf{}}}}{w}):=- \mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{x},{{\boldsymbol{\mathbf{}}}}{y} \sim \mathcal{D}}\left[ \log p_{{{\boldsymbol{\mathbf{}}}}{w}}({{\boldsymbol{\mathbf{}}}}{y}|{{\boldsymbol{\mathbf{}}}}{x})\right]\,.
    \label{eq:marginal_likelihood}\end{aligned}$$

Under stochastic optimisation algorithms, such as Stochastic Gradient Descent (SGD), we estimate equation by sampling a mini-batch of data-label pairs $\mathcal{B} \subset \mathcal{D}$. $$\mathcal{L}_{\mathrm{SGD}}(\mathcal{B}; {{\boldsymbol{\mathbf{}}}}{\theta}) = -\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{x},{{\boldsymbol{\mathbf{}}}}{y} \sim \mathcal{B}} \log p_{{{\boldsymbol{\mathbf{}}}}{\theta}}({{\boldsymbol{\mathbf{}}}}{y}|{{\boldsymbol{\mathbf{}}}}{x}) \approx \mathcal{L}( \mathcal{D}; {{\boldsymbol{\mathbf{}}}}{\theta})$$

Consider an $L$ layer network with no noise injections. We obtain the activations ${{{\boldsymbol{\mathbf{}}}}{h}} = \{{{{\boldsymbol{\mathbf{}}}}{h}}_0, ... , {{{\boldsymbol{\mathbf{}}}}{h}}_{L-1} \}$, where ${{\boldsymbol{\mathbf{}}}}{h}_{0}={{\boldsymbol{\mathbf{}}}}{x}$ is the input data. Let ${{\boldsymbol{\mathbf{}}}}{\epsilon}$ be the set of noise injections at each layer: ${{\boldsymbol{\mathbf{}}}}{\epsilon} = \{{{{\boldsymbol{\mathbf{}}}}{\epsilon}}_0, ... , {{{\boldsymbol{\mathbf{}}}}{\epsilon}}_{L-1} \}$. When performing a noise injection procedure, the value of the next layer’s activations depends on the noised value of the previous layer. We denote the intermediate, soon-to-be-noised value of an activation, for a network consisting of dense layers (a.k.a. a multi-layer perceptron: MLP) as $\hat{{{\boldsymbol{\mathbf{}}}}{h}}_{i}$ and the subsequently noised value as $\tilde{{{\boldsymbol{\mathbf{}}}}{h}}_{i}$: $$\begin{aligned}
\hat{{{\boldsymbol{\mathbf{}}}}{h}}_{i}= \phi\left(\tilde{{{\boldsymbol{\mathbf{}}}}{h}}_{i-1}{{\boldsymbol{\mathbf{}}}}{W}_{i}\right)\,, \qquad  
\tilde{{{\boldsymbol{\mathbf{}}}}{h}}_{i} = \hat{{{\boldsymbol{\mathbf{}}}}{h}}_{i} \circ {{\boldsymbol{\mathbf{}}}}{\epsilon}_{i}\,, 
\label{eq:noise_recursion}\end{aligned}$$ where $\phi$ is some non-linearity, and $\circ$ is some element-wise operation. We can, for example, add or multiply Gaussian noise to each hidden layer unit. The multiplicative case can be rewritten as an activation-scaled addition[^1]. Both these cases give us: $$\begin{aligned}
\tilde{{{\boldsymbol{\mathbf{}}}}{h}}_i &= \hat{{{\boldsymbol{\mathbf{}}}}{h}}_i + {{\boldsymbol{\mathbf{}}}}{\epsilon}_i, 
\qquad {{\boldsymbol{\mathbf{}}}}{\epsilon}_i \sim \mathcal{N}(0,{{\boldsymbol{\mathbf{}}}}{\sigma}_i^2), \label{eq:noise_add} \\
\tilde{{{\boldsymbol{\mathbf{}}}}{h}}_i &= \hat{{{\boldsymbol{\mathbf{}}}}{h}}_i + {{\boldsymbol{\mathbf{}}}}{\epsilon}_i, 
\qquad {{\boldsymbol{\mathbf{}}}}{\epsilon}_i \sim \mathcal{N}\left(0,\hat{{{\boldsymbol{\mathbf{}}}}{h}}_i^2{{\boldsymbol{\mathbf{}}}}{\sigma}_i^2\right).
\label{eq:noise_mult}\end{aligned}$$

Here we focus our analysis on noise *additions*, but through equation we can translate our results to the multiplicative case. The expected cost function to be minimized by SGD is estimated by Jensen’s inequality such that: $$\begin{aligned}
\tilde{\mathcal{L}}_{\mathrm{SGD}}(\mathcal{B}; {{\boldsymbol{\mathbf{}}}}{\theta}) \geq \mathbb{E}_{\substack{\tilde{{{\boldsymbol{\mathbf{}}}}{h}} \sim p(\tilde{{{\boldsymbol{\mathbf{}}}}{h}}|{{\boldsymbol{\mathbf{}}}}{h}) \\
{{\boldsymbol{\mathbf{}}}}{x},{{\boldsymbol{\mathbf{}}}}{y} \sim \mathcal{B}}} \log   p_{{{\boldsymbol{\mathbf{}}}}{\theta}}({{\boldsymbol{\mathbf{}}}}{y}|\tilde{{{\boldsymbol{\mathbf{}}}}{h}}) 
\label{eq:invariant_form_jensen}\end{aligned}$$ with $p(\tilde{{{\boldsymbol{\mathbf{}}}}{h}}|\hat{{{\boldsymbol{\mathbf{}}}}{h}})=\prod_i^L p(\tilde{{{\boldsymbol{\mathbf{}}}}{h}}_i|\hat{{{\boldsymbol{\mathbf{}}}}{h}}_i, {{\boldsymbol{\mathbf{}}}}{\epsilon}_i)$ and $p(\tilde{{{\boldsymbol{\mathbf{}}}}{h}}_i|\hat{{{\boldsymbol{\mathbf{}}}}{h}}_i, {{\boldsymbol{\mathbf{}}}}{\epsilon}_i)=  \mathcal{N}(\tilde{{{\boldsymbol{\mathbf{}}}}{h}_i}|\hat{{{\boldsymbol{\mathbf{}}}}{h}}_i,{{\boldsymbol{\mathbf{}}}}{\sigma}^2)$.

We restrict our analysis to networks with non-linearities that are broadly moment preserving [@Vladimirova2019UnderstandingLevel], such that the accumulated noise at a given layer is still Gaussian. The properties of the accumulated noise on $\tilde{{{\boldsymbol{\mathbf{}}}}{h}}_n$ depends on the properties of the activation function $\phi(\cdot)$ used in the network. For instance, Gaussian noise becomes sub-Gaussian when fed through a bounded activation function, such as a sigmoid, by Hoeffding’s Lemma [@Bentkus2004; @Vladimirova2019UnderstandingLevel]. A class of functions that are broadly moment preserving are those that obey the *extended envelope property*, as they have *asymptotically equivalent moments*, see Appendix A and [@Vladimirova2019UnderstandingLevel]. Activation functions that satisfy this property are among the most widely used in deep learning (Rectified Linear Units ($\mathrm{ReLU}$) and Exponential Linear Units ($\mathrm{ELU}$) among others). Even more complex operations, such as max-pooling in convolutional networks, satisfy this property [@Vladimirova2019UnderstandingLevel].

### Sensitivity and Classification Margins {#sec:sensitivity_background}

Sophisticated deep learning models can be brittle [@Szegedy2014; @Shamir2019; @Goodfellow2015; @Papernot2016; @Moosavi-Dezfooli2016]. Adding small unstructured noise to data, perturbations that would not fool a human can alter neural network predictions. A model’s weakness to perturbations is called the *sensitivity* of the model.

Typically, models with larger classification margins are less sensitive to input perturbations [@Sokolic2017; @Jakubovitz2018; @Cohen2019; @Liu2018; @Li2018]. Such margins are the distance in data-space between a point ${{\boldsymbol{\mathbf{}}}}{x}$ and a classifier’s decision boundary. Larger margins mean that a classifier associates a larger region centered on a point ${{\boldsymbol{\mathbf{}}}}{x}$ to the same class. Intuitively this means that noise added to ${{\boldsymbol{\mathbf{}}}}{x}$ is still likely to fall within this region, leaving the classifier prediction unchanged.

Tikhonov Regularisation
-----------------------

Tikhonov regularisation involves adding some regulariser to the loss function, which encodes a notion of ‘smoothness’ of a function $f$ [@Bishop1995]. As such, by design, regularisers of this form have been shown to have beneficial regularisation properties when used in the training objective of neural networks by smoothing the loss landscape [@Girosi1990; @Burger2003]. If we have a loss of the form $\tilde{\mathcal{L}}_{\mathrm{SGD}}(\mathcal{B}; {{\boldsymbol{\mathbf{}}}}{\theta})$, the Tikhonov regularised loss becomes: $$\tilde{\mathcal{L}}_{\mathrm{SGD}}(\mathcal{B}; {{\boldsymbol{\mathbf{}}}}{\theta}) + \lambda\|f_{\theta}\|^2_{\mathcal{H}}$$ where $f_{\theta}$ is the function with parameters $\theta$ which we are learning and $\|\cdot\|_{\mathcal{H}}$ is the norm in the function space $\mathcal{H}$ and $\lambda$ is a (multidimensional) penalty which penalises elements of $\|f_{\theta}\|^2_{\mathcal{H}}$ unequally, or is data-dependent [@Tikhonov1977; @Bishop1995].

The Explicit Effect of Gaussian Noise Injections
================================================

We can express the effect of the Gaussian noise injection on the cost function as: $$\tilde{\mathcal{L}}_{\mathrm{SGD}}(\mathcal{B};{{\boldsymbol{\mathbf{}}}}{\theta}, {{\boldsymbol{\mathbf{}}}}{\epsilon}) =  \mathcal{L}_{\mathrm{SGD}}(\mathcal{B}; {{\boldsymbol{\mathbf{}}}}{\theta}) + \Delta\mathcal{L}({{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L})$$ where ${{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}$ is the noise accumulated on the final layer $L$ from the noise additions ${{\boldsymbol{\mathbf{}}}}{\epsilon}$ on the previous hidden layer activations. We formulate an estimate for $\Delta \mathcal{L}$ that is architecture independent and is positive semi-definite (PSD), and thus a viable regulariser. It is important that the regularisation term we add to our loss is PSD. Non-PSD regularisers can change sign from batch to batch and do not give a consistent objective to optimise, making them unfit as regularisers [@Botev2017; @Sagun2018; @Wei2020].

Here we consider the case where we noise all layers except the *final* predictive layer which we consider to have no activation function for the purposes of our analysis. First and foremost we need an analytic form of the accumulated noise ${{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}$. We denote each layer’s Jacobian as ${{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x}) \in \mathbb{R}^{N_L \times N_k}, \quad {{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x})_{i,j} = \partial {{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x})_i / \partial {{\boldsymbol{\mathbf{}}}}{h}_{k}({{\boldsymbol{\mathbf{}}}}{x})_j$ and have:

Consider an $L$ layer neural network, experiencing GNIs ${{\boldsymbol{\mathbf{}}}}{\epsilon}_{i}$ at each layer $i \in [0, \dots, L - 1]$. Assuming the Hessians, of the form $\nabla^2 {{\boldsymbol{\mathbf{}}}}{h}_{l}({{\boldsymbol{\mathbf{}}}}{x})\rvert_{{{\boldsymbol{\mathbf{}}}}{h}_{m}({{\boldsymbol{\mathbf{}}}}{x})}, m < l$ where $l,m$ index over layers, are finite, the noise accumulated at the final layer, ${{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}$, is given by: $${{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L} = \sum_{k=0}^{L-1}{{\boldsymbol{\mathbf{}}}}{\epsilon}_{k}{{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x}) + \mathcal{O} ({{\boldsymbol{\mathbf{}}}}\gamma)$$ $\mathcal{O} ({{\boldsymbol{\mathbf{}}}}\gamma)$ represents higher order terms in ${{\boldsymbol{\mathbf{}}}}{\epsilon}$ that tend to zero in the limit of small variance noise. \[prop:accum\_noise\]

See Appendix \[app:accum\_noise\] for a proof. Using this result we can now derive the explicit regulariser induced by GNIs by marginalising out the injected noise ${{\boldsymbol{\mathbf{}}}}{\epsilon}$ from the term $\Delta \mathcal{L}({{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L})$.

Consider an $L$ layer neural network experiencing GNIs ${{\boldsymbol{\mathbf{}}}}{\epsilon}_{i}$ at each layer $i \in [0, \dots, L - 1]$. We assume the third derivative of $\mathcal{L}({{\boldsymbol{\mathbf{}}}}{x})$ w.r.t ${{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x})$ is finite, which is the case for most loss functions (see proof for details). Using the form of ${{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}$ in Proposition \[prop:accum\_noise\]; we can marginalise out the injected noise ${{\boldsymbol{\mathbf{}}}}{\epsilon}$ to obtain an added regulariser: $$\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{\epsilon}} \left[ \Delta \mathcal{L}({{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}) \right] = \frac{1}{2}\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{x} \sim \mathcal{B}} \left[ \sum_{k=0}^{L-1}\sigma^2_k \mathrm{Tr}({{\boldsymbol{\mathbf{}}}}{J}^T_{k}({{\boldsymbol{\mathbf{}}}}{x})
    {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x}))  \right] +   O({{\boldsymbol{\mathbf{}}}}{\kappa}) 
    \label{eqn:deltaL}$$ where $\mathcal{L}({{\boldsymbol{\mathbf{}}}}{x})$ is the loss for a datapoint ${{\boldsymbol{\mathbf{}}}}{x}$, ${{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x})$ is $\nabla^2 \mathcal{L}\rvert_{{{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x})} \in \mathbb{R}^{N_L
\times N_L}$ ($N_L$ is the number of output neurons), and $O({{\boldsymbol{\mathbf{}}}}{\kappa})$ represents higher order terms in ${{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}$ that disappear in the limit of small variance noise. \[prop:explicit\_reg\]

See Appendix \[app:exp\_reg\] for a proof of this. For the remainder of this work, we assume that noise injections are of small variance and that as such, the higher order terms of Proposition \[app:accum\_noise\] and Theorem \[prop:explicit\_reg\] can be ignored. In the case of Proposition \[prop:accum\_noise\], for many of the functions that obey the extended envelope property (eg. $\mathrm{ReLU}$), which are locally linear, this result is locally exact. $\mathcal{O} ({{\boldsymbol{\mathbf{}}}}\gamma)$ is exactly zero if the variance of the noise is smaller than what is required to shift outside the linear region.

In Appendix \[app:exp\_reg\], we also show that the first set of terms in $R$ corresponds to the sum of the traces of the Gauss-Newton approximations of the Hessians ${{\boldsymbol{\mathbf{}}}}{H}_{k}({{\boldsymbol{\mathbf{}}}}{x})$, i.e the second order derivatives of the loss with respect to each layer’s activations. Smaller traces are measures of a smoother loss landscape and conversely larger traces indicate a ‘peakier’ and steeper landscape. Explicitly penalising this Hessian means that we are more likely to land in wider (smoother) minima, which has been shown, although this is a point of contention [@Dinh2017], to induce networks with better generalisation properties [@Keskar2019OnMinima; @Jastrzebski2017]. Though this offers one explanation for the regularisation offered by GNIs, we can glean more insight into the mechanism underpinning this regularisation by considering the cases of regression and classification models separately, which we do in the following section.

**Regularisation in Regression.** In the case of regression one of the most commonly used loss functions is the mean-squared error (MSE), which is defined for a point $n$ as follows: $$\mathcal{L}({{\boldsymbol{\mathbf{}}}}{x}) = ({{\boldsymbol{\mathbf{}}}}{y}_n - {{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x}))^2
\label{eq:mse}$$ For this loss, the Hessians in Theorem \[prop:explicit\_reg\] are simply the identity matrix. Ignoring the higher order terms in Theorem \[prop:explicit\_reg\] for small variance noise injections, the explicit regularisation term, which is guaranteed to be PSD, is: $$R \approx \frac{1}{2} \mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{x} \sim \mathcal{B}} \left[ \sum_{k=0}^{L-1} \sigma^2_k (\|{{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x})\|^2_F) \right]
    \label{eq:mse_reg_term}$$

where $\sigma^2_k$ is the variance of the noise ${{\boldsymbol{\mathbf{}}}}{\epsilon}_{k}$ injected at layer $k$ and $||\cdot||_F$ is the Frobenius norm. See Appendix \[app:exp\_reg\_regression\] for a proof.

**Regularisation in Classification.** In the case of classification, we consider the case of a cross-entropy (CE) loss. Recall that we consider our network outputs ${{\boldsymbol{\mathbf{}}}}{h}_L$ to be the pre-$\mathrm{softmax}$ of logits of the final layer ${{\boldsymbol{\mathbf{}}}}{L}$. We denote ${{\boldsymbol{\mathbf{}}}}{p}({{\boldsymbol{\mathbf{}}}}{x})=\mathrm{softmax}({{\boldsymbol{\mathbf{}}}}{h}_L({{\boldsymbol{\mathbf{}}}}{x}))$. The loss is thus: $$\mathcal{L}({{\boldsymbol{\mathbf{}}}}{x}) = - \sum_{c=0}^M {{\boldsymbol{\mathbf{}}}}{y}_{n,c} \log ({{\boldsymbol{\mathbf{}}}}{p}({{\boldsymbol{\mathbf{}}}}{x}))_c),
\label{eq:cross_entropy}$$ where $c$ indexes over the $M$ possible classes of the classification problem. The hessian ${{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x})$, which is $\nabla^2 \mathcal{L}({{\boldsymbol{\mathbf{}}}}{x})\rvert_{{{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x})} \in \mathbb{R}^{N_L \times N_L}$, is easy to compute and has the form: $${{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x})_{i,j} = 
\begin{cases}
{{\boldsymbol{\mathbf{}}}}{p}({{\boldsymbol{\mathbf{}}}}{x})_i(1 - {{\boldsymbol{\mathbf{}}}}{p}({{\boldsymbol{\mathbf{}}}}{x})_j) & i = j \\
-{{\boldsymbol{\mathbf{}}}}{p}({{\boldsymbol{\mathbf{}}}}{x})_i{{\boldsymbol{\mathbf{}}}}{p}({{\boldsymbol{\mathbf{}}}}{x})_j & i \neq j \\
\end{cases}$$ This Hessian is itself PSD and if we ignore the higher order terms in Theorem \[prop:explicit\_reg\] for small variance noise injections, the explicit regularisation term, which is PSD can be approximated by: $$R \approx \frac{1}{2} \mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{x} \sim \mathcal{B}} \left[ \sum_{k=0}^{L-1} \sigma_k^2 \sum_{i,j}(\mathrm{diag}( {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}))^T{{{\boldsymbol{\mathbf{}}}}{J}}^2_{k}({{\boldsymbol{\mathbf{}}}}{x}) )_{i,j} \right]
\label{eq:ce_reg_term}$$ where as before $\sigma^2_k$ is the variance of the noise ${{\boldsymbol{\mathbf{}}}}{\epsilon}_{k}$ injected at layer $k$. See Appendix \[app:exp\_reg\_classification\] for a detailed demonstration of this.

\

In Figure \[fig:marg\_approx\] we show that models trained with $R$ and GNIs have similar test-set loss and parameter Hessians throughout training, even when using ELU activations which test the first-order approximations for small variance noise we have made. This means that models trained with $R$ and GNIs have almost identical trajectories through the loss landscape and that they have smoother trajectories (smaller Hessian trace) than models trained without regularisation.

In both classification and regression the explicit regulariser penalises the norm of the Jacobians $\|{{\boldsymbol{\mathbf{}}}}{J}_k\|^2_F$. As we now show, this corresponds to a penalisation in the Fourier domain.

Fourier Domain Regularisation
=============================

In this section, we show that the regularisers we have derived are equivalent to Tikhonov regularisation in the Hilbert-Sobolev space with index $(1,2)$, which is formally defined as follows:

The Sobolev space of index (k,p), denoted $W^{k,p}(\Omega), \Omega \subset \mathbb{R}^d$, where $k$ is a non negative integer and $p\geq1$, is the space of locally integrable functions $f: \Omega \to \mathbb{R}$ such that for every $\alpha < k$ the weak derivative $D^{\alpha}f$ exists and $D^{\alpha}f \in L^p(\Omega)$. The norm in such a space is given by $\|v\|_{W^{k,p}(\Omega)} = (\sum_{|\alpha|_1 \leq k} \int_{\Omega} |D^{\alpha}f(x)|^p dx)^{\frac{1}{p}}$.

[@Hornik1991] have shown that neural networks with continuous activations, which have continuous derivatives up to order $k$, such as the sigmoid function, are universal approximators in the Sobolev spaces of order $k$. Recently [@Czarnecki2017] have shown that networks that use $\mathrm{ReLU}$-like activation functions are *also* universal approximators in the Sobolev spaces of order 1. Hence, we can view a neural network to be a parameter that indexes a function in the Sobolev space with index $(1,2)$, i.e. $f_{{{\boldsymbol{\mathbf{}}}}{\theta}} \in {W^{1,2}(\Omega)}$.

For a function, $f_{\theta}$, with a $d$-dimensional input and a single output if $f_{\theta} \in {W^{1,2}(\Omega)}, \Omega \subset \mathbb{R}^d$ holds, then we have that [@Novak]: $$\int_\Omega \|D f_{\theta}({{\boldsymbol{\mathbf{}}}}{x})\|^2_{L_2(\Omega)}d{{\boldsymbol{\mathbf{}}}}{x} = \int_\Omega (\prod_{j=1}^d{{\boldsymbol{\mathbf{}}}}{\omega}_j^2)|\mathcal{F}({{\boldsymbol{\mathbf{}}}}{\omega})|^2 d{{\boldsymbol{\mathbf{}}}}{\omega}, \Omega \subset \mathbb{R}^d
\label{eq:sobolev_norm}$$ where $\mathcal{F}$ is the Fourier transform of $f_{\theta}$ and $j$ indexes over ${{\boldsymbol{\mathbf{}}}}{\omega} = [\omega_1, \dots, \omega_d]$. In the case where the dataset contains finitely many points, the integral is approximated by a summation over datapoints and input dimensions, and in the case of SGD this summation occurs over a sampled batch. Further, in our case, the regularisation is applied to each intermediate set of activations and to each output neuron. Assuming differentiable and continuous activation functions, then the Jacobians in the explicit regulariser (equations and ) are equivalent to the weak derivative, giving: $$\begin{aligned}
 R &= \frac{1}{2}\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{x} \sim \mathcal{B}} \left[ \sum_{k=0}^{L-1} \sum_{i} \lambda_i \|D f^k_{\theta}({{\boldsymbol{\mathbf{}}}}{h}_{k}({{\boldsymbol{\mathbf{}}}}{x}))_i\|^2_{L_2(\Omega)} \right] \\
 &\approx \frac{1}{2}\sum_{k=0}^{L-1} \sum_{i} \int_\Omega \lambda_i (\prod_{j=1}^d{{\boldsymbol{\mathbf{}}}}{\omega}_j^2)|\mathcal{F}^k_i({{\boldsymbol{\mathbf{}}}}{\omega})|^2 d{{\boldsymbol{\mathbf{}}}}{\omega}, \Omega \subset \mathbb{R}^d
 \label{eq:sobolev_reg}\end{aligned}$$ where ${{\boldsymbol{\mathbf{}}}}{h}_0 ={{\boldsymbol{\mathbf{}}}}{x}$, $i$ indexes over output neurons, $f^k_{\theta}(\cdot)_i$ is the function from layer $k$ to the $i^{\mathrm{th}}$ network output, $d_k$ is the dimensionality of layer $k$, and $\mathcal{F}^k_{i}$ is the Fourier transform of $f^k_{\theta}(\cdot)_i$.

In our setting $\lambda_i = \sigma_k^2$ is constant for regression. For classification, $\lambda_i$ is data-dependent $\lambda_i = (\sigma_k^2\mathrm{diag}({{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}))^T)_i$ and becomes a part of the integral as in [@Tikhonov1977] and @Bishop1995. Generally our derived $R$ are equivalent to penalising the term associated with $\alpha = 1$ in the Sobolev norm with coefficient $\lambda_i$, and penalising the term associated with $\alpha = 0$ with coefficient 0 (see definition above for the Sobolev norm). This follows the Tikhonov regularisation definition of [@Tikhonov1977] and @Bishop1995 which allows for different coefficients to weight each term of the sum that defines the Hilbert space norm and allows for data-dependent $\lambda_i$.

As such, the explicit regulariser is the sum of the Fourier regularisers for each *subnetwork* within a network. Training with Gaussian noise is equivalent to a prior in Sobolev space which favours smooth functions with low-frequency components. We demonstrate this empirically in Figure \[fig:network\_fourier\] for models trained with noise and the explicit regulariser. Clearly, these models learn a function with low-frequency components, that doesn’t overfit to data.

![As in @Rahaman2019, we train 6-layer deep 256-unit wide ReLU networks trained to regress the function $\lambda(z) = \sum_i \sin(2\pi k_iz+\phi(i))$ with $k_i  \in(5,10,\dots,45,50)$. We train these networks with no noise (Baseline), with GNIs of variance $\sigma^2=0.1$ injected into each layer except the final layer (Noise), and with the explicit regulariser for regression in (Exp Reg). The first row shows the Fourier spectrum (x-axis) of the networks (calculated using Lemmas 1 and 2 of @Rahaman2019) as training progresses (y-axis) averaged over 10 training runs. Colours show each frequency’s amplitude clipped between 0 and 1. The second row shows samples of randomly generated target functions and the function learnt by the networks. []{data-label="fig:network_fourier"}](figures/spectral_dynamics_baseline10x.png "fig:"){width="28.00000%"} ![As in @Rahaman2019, we train 6-layer deep 256-unit wide ReLU networks trained to regress the function $\lambda(z) = \sum_i \sin(2\pi k_iz+\phi(i))$ with $k_i  \in(5,10,\dots,45,50)$. We train these networks with no noise (Baseline), with GNIs of variance $\sigma^2=0.1$ injected into each layer except the final layer (Noise), and with the explicit regulariser for regression in (Exp Reg). The first row shows the Fourier spectrum (x-axis) of the networks (calculated using Lemmas 1 and 2 of @Rahaman2019) as training progresses (y-axis) averaged over 10 training runs. Colours show each frequency’s amplitude clipped between 0 and 1. The second row shows samples of randomly generated target functions and the function learnt by the networks. []{data-label="fig:network_fourier"}](figures/spectral_dynamics_noise10x.png "fig:"){width="28.00000%"} ![As in @Rahaman2019, we train 6-layer deep 256-unit wide ReLU networks trained to regress the function $\lambda(z) = \sum_i \sin(2\pi k_iz+\phi(i))$ with $k_i  \in(5,10,\dots,45,50)$. We train these networks with no noise (Baseline), with GNIs of variance $\sigma^2=0.1$ injected into each layer except the final layer (Noise), and with the explicit regulariser for regression in (Exp Reg). The first row shows the Fourier spectrum (x-axis) of the networks (calculated using Lemmas 1 and 2 of @Rahaman2019) as training progresses (y-axis) averaged over 10 training runs. Colours show each frequency’s amplitude clipped between 0 and 1. The second row shows samples of randomly generated target functions and the function learnt by the networks. []{data-label="fig:network_fourier"}](figures/spectral_dynamics_exp_reg10x.png "fig:"){width="28.00000%"}\

GNIs and Model Sensitivity {#sec:sensitivity}
==========================

@Rahaman2019 showed empirically that functions biased towards lower frequencies in the Fourier domain are less sensitive to noisy data, and there is ample evidence demonstrating that models trained with noised data are less sensitive to perturbations by inducing larger classification margins [@Cohen2019; @Liu2018; @Li2018].

Given this, our aim is twofold. First, we want to establish a theoretical understanding of the link between a classifier’s Fourier spectrum and its classification margins. Second, because noisy data has already been shown to induce larger classification margins, we want to ascertain whether extending noise to all intermediate layers of a network further decreases network sensitivity.

Classifier Sensitivity
----------------------

As discussed previously, a model’s sensitivity to noise is intrinsically linked to the size of its classification margins. Intuitively a classification margin can be described as the minimal perturbation norm necessary to fool the classification function. @Sokolic2017 and @Jakubovitz2018 define a classification margin $M$ that is the radius of the largest metric ball centered on a point ${{\boldsymbol{\mathbf{}}}}{x}$ to which a classifier assigns ${{\boldsymbol{\mathbf{}}}}{y}$, the true label.

Consider a classifier that outputs a correct prediction for the true class $A$ associated with a point ${{\boldsymbol{\mathbf{}}}}{x}$. Then the first order approximation for the l2-norm of the classification margin $M$, which is the minimal perturbation necessary to fool a classifier, is lower bounded by: $$\begin{aligned}
M({{\boldsymbol{\mathbf{}}}}{x}) \geq \frac{({{\boldsymbol{\mathbf{}}}}{h}^A_{L}({{\boldsymbol{\mathbf{}}}}{x}) - {{\boldsymbol{\mathbf{}}}}{h}^B_{L}({{\boldsymbol{\mathbf{}}}}{x}) )}{\sqrt{2} \|{{\boldsymbol{\mathbf{}}}}{J}_0({{\boldsymbol{\mathbf{}}}}{x}))|_F}. 
\label{eq:radius_robust_frob}\end{aligned}$$ We have ${{\boldsymbol{\mathbf{}}}}{h}^A_{L}({{\boldsymbol{\mathbf{}}}}{x}) \geq {{\boldsymbol{\mathbf{}}}}{h}^B_{L}({{\boldsymbol{\mathbf{}}}}{x})$, where ${{\boldsymbol{\mathbf{}}}}{h}^A_{L}({{\boldsymbol{\mathbf{}}}}{x})$ is the $L^{th}$ layer activation (pre-softmax) associated with the true class $A$, and ${{\boldsymbol{\mathbf{}}}}{h}^B_{L}({{\boldsymbol{\mathbf{}}}}{x})$ is the second largest $L^{th}$ layer activation.

Networks that have lower-frequency spectrums and consequently have smaller norms of Jacobian (as in equation ), will have larger classification margins and will be less sensitive to perturbations. This broadly follows our prior discussion, but expresses it mathematically and explains the empirical observations of @Rahaman2019.

What does this entail for GNIs applied to each layer of a network ? We can view the penalisation of the norms of the Jacobians, induced by GNIs for each layer $k$, as an unpweighted penalisation of $\|{{\boldsymbol{\mathbf{}}}}{J}_0({{\boldsymbol{\mathbf{}}}}{x})\|_F$. By the chain rule ${{\boldsymbol{\mathbf{}}}}{J}_0$ can be expressed in terms of any of the other network Jacobians ${{\boldsymbol{\mathbf{}}}}{J}_0({{\boldsymbol{\mathbf{}}}}{x}) = {{\boldsymbol{\mathbf{}}}}{J}_k({{\boldsymbol{\mathbf{}}}}{x})\frac{\partial {{\boldsymbol{\mathbf{}}}}{h}_k}{{{\boldsymbol{\mathbf{}}}}{x}} \forall k \in [0 \dots L]$. We can write $\|{{\boldsymbol{\mathbf{}}}}{J}_0({{\boldsymbol{\mathbf{}}}}{x})\|_F = \|{{\boldsymbol{\mathbf{}}}}{J}_k({{\boldsymbol{\mathbf{}}}}{x})\frac{\partial {{\boldsymbol{\mathbf{}}}}{h}_k}{{{\boldsymbol{\mathbf{}}}}{x}}\|_F \leq \|{{\boldsymbol{\mathbf{}}}}{J}_k({{\boldsymbol{\mathbf{}}}}{x})\|_F\|\frac{\partial {{\boldsymbol{\mathbf{}}}}{h}_k}{{{\boldsymbol{\mathbf{}}}}{x}}\|_F$. Minimising $\|{{\boldsymbol{\mathbf{}}}}{J}_0({{\boldsymbol{\mathbf{}}}}{x})\|_F$ is equivalent to minimising $\|{{\boldsymbol{\mathbf{}}}}{J}_k({{\boldsymbol{\mathbf{}}}}{x})\|_F$ and $\|\frac{\partial {{\boldsymbol{\mathbf{}}}}{h}_k}{{{\boldsymbol{\mathbf{}}}}{x}}\|_F$, and upweighted penalisations of $\|{{\boldsymbol{\mathbf{}}}}{J}_k({{\boldsymbol{\mathbf{}}}}{x})\|_F$ should translate into a shrinkage of $\|{{\boldsymbol{\mathbf{}}}}{J}_0({{\boldsymbol{\mathbf{}}}}{x})\|_F$.

To illustrate this, consider the artificial situation where the Jacobians of the network are invertible. Recall that by the chain rule ${{\boldsymbol{\mathbf{}}}}{J}_i({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{J}_0^{-1}({{\boldsymbol{\mathbf{}}}}{x}) = (\frac{\partial {{\boldsymbol{\mathbf{}}}}{h}_k}{{{\boldsymbol{\mathbf{}}}}{x}})^{-1} = \nabla {{\boldsymbol{\mathbf{}}}}{h}_k^{-1}\rvert_{{{\boldsymbol{\mathbf{}}}}{x}}$ (for any k indexing over layers). Given the form of the accumulated noise in Proposition \[prop:accum\_noise\], we can state that adding noise to each layer is broadly equivalent to adding the following noise to data: $$\epsilon_0 \sim \mathcal{N}(0,{{\boldsymbol{\mathbf{}}}}{\sigma}_0^2 +  \sum_{i=1}^{L-1}(\nabla {{\boldsymbol{\mathbf{}}}}{h}_k^{-1}\rvert_{{{\boldsymbol{\mathbf{}}}}{x}}{{\boldsymbol{\mathbf{}}}}{\sigma}^T_k)^2))
\label{eq:equiv_variance_input}$$

In this case, training with GNIs is equivalent to training with *anisotropic* noise with variance ${{\boldsymbol{\mathbf{}}}}{\Sigma}({{\boldsymbol{\mathbf{}}}}{x})^2 = {{\boldsymbol{\mathbf{}}}}{\sigma}^2(1 +  \sum_{i=1}^{L-1}(\nabla {{\boldsymbol{\mathbf{}}}}{h}_k^{-1}\rvert_{x})^2)$. For an equivalent model that only adds noise with variance ${{\boldsymbol{\mathbf{}}}}{\sigma}^2$ to data, we can state that ${{\boldsymbol{\mathbf{}}}}{\Sigma}^2({{\boldsymbol{\mathbf{}}}}{x}) \geq {{\boldsymbol{\mathbf{}}}}{\sigma}^2$, element-wise. Substituting ${{\boldsymbol{\mathbf{}}}}{\Sigma}^2({{\boldsymbol{\mathbf{}}}}{x})$ into equation , we can see that this translates into a greater penalisation of $\|{{\boldsymbol{\mathbf{}}}}{J}_0({{\boldsymbol{\mathbf{}}}}{x})\|_F$. As such, noising each layer should induce a smaller $\|{{\boldsymbol{\mathbf{}}}}{J}_0({{\boldsymbol{\mathbf{}}}}{x})\|_F$, and larger classification margins than solely noising data. We support this empirically in Figure \[fig:adv\_attack\] and show that the regularisers $R$ replicate this effect.

Calibration
-----------

\

A neural network classifier gives a prediction $\hat{y}({{\boldsymbol{\mathbf{}}}}{x})$ with confidence $\hat{p}({{\boldsymbol{\mathbf{}}}}{x})$ (the probability attributed to that prediction) for a datapoint ${{\boldsymbol{\mathbf{}}}}{x}$. Perfect calibration consists of being as likely to be correct as you are confident: $p(\hat{y}=y|\hat{p}=r)=r, \quad \forall r\in[0,1]$

In Appendix \[app:fourier\_capacity\], we establish a connection between the Fourier domain and model capacity, which is a measure of model complexity, showing that models with low-frequency components correspond to less ‘complex’ models with lower capacity. @Guo2017 show that models with lower capacity are better calibrated and output a broader range of prediction probabilities. We show that this holds true for models trained with GNIs and with the explicit regulariser, in Figure \[fig:calibration\]. We leave establishing a formal connection between the Fourier domain and calibration for future work.

Related Work
============

Many variants of GNIs have been proposed to regularise neural networks. @Poole2014 extend this process to its logical conclusion and apply noise to all computational steps in a neural network layer. Not only is noise applied to the layer input it is applied to the layer output and to the pre-activation function logits. The authors allude to explicit regularisation but only derive a result for a single layer auto-encoder with a single noise injection. Similarly, @Bishop1995 derive an analytic form for the explicit regulariser induced by noise injections on *data* and show that such injections are equivalent to Tikhonov regularisation in an unspecified function space.

Recently @Wei2020 conducted similar analysis to ours, dividing the effects of Bernoulli dropout into *explicit* and *implicit* effects. Their work is built on that of @Mele1993 and @Helmbold2015 who perform this analysis for linear neural networks. @Arora2020 derive an explicit regulariser for Bernoulli dropout on the final layer of a neural network.

Conclusion
==========

In this work, we derived analytic forms for the explicit regularisation induced by Gaussian noise injections. Having characterised the explicit regulariser as a form of Tikhonov regularisation which penalises networks with high-frequency content in the Fourier space, we leverage this knowledge to show that noise injections produce wider classification margins and better calibration.

Impact Statement {#impact-statement .unnumbered}
================

This paper uncovers a new mechanism by which a widely used regularisation method operates and paves the way for designing new regularisation methods which take advantage of our findings. Regularisation methods produce models that are not only less likely to overfit, but also have better calibrated predictions that are more robust to adversarial attack. As such improving our understanding of such methods is critical as machine learning models become increasingly ubiquitous and embedded in decision making.

Extended Envelope Property {#app:eep}
==========================

As described in @Vladimirova2019UnderstandingLevel, we make us of the Extended Envelope Property.

[Extended Envelope Property:]{} A non-linear function $\phi:\mathbb{R} \to \mathbb{R}$ is said to obey the *extended envelope property* if $\exists~c_1,c_2 \geq 0 ~ d_1, d_2 \geq 0 $ such that:

-   $|\phi(x)| \geq c_1 +d_1|x|, ~ \forall x\in \mathbb{R}^+ ~\mathrm{or} ~ x\in \mathbb{R}^-$

-   $|\phi(x)| \leq c_2 +d_2|x|, \forall x\in \mathbb{R}$

Further we also use the related idea of Asymptotic Equivalence.

[Asymptotic Equivalence]{} Two sequences $a_k$, $b_k$ are asymptotically equivalent ($a_k \asymp b_k$) if $\exists d>0, D>0$ such that: $d \leq \frac{a_k}{b_k} \leq D ~ \forall k \in \mathbb{N}$.

For any non-linearity $\phi: \mathbb{R} \to \mathbb{R} $ that obeys the extended envelope property, then for any *symmetric* random variable ${{\boldsymbol{\mathbf{}}}}{X}$ we can write:

$$\|\phi({{\boldsymbol{\mathbf{}}}}{X})\|_k \asymp \|{{\boldsymbol{\mathbf{}}}}{X}\|_k ~ \forall k \geq 1
\label{eqapp:moment-preservation}$$

where $\|\cdot\|_k$ is the $k^{th}$ norm.

Accumulated Noise Derivation {#app:accum_noise}
============================

For an $L$ layer neural network, experiencing GNIs ${{\boldsymbol{\mathbf{}}}}{\epsilon}_{i}$ at each layer $i \in [0, \dots, L - 1]$, the the noise accumulated at the final layer, ${{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}$, is given by: $${{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L} = \sum_{k=0}^{L-1}{{\boldsymbol{\mathbf{}}}}{\epsilon}_{k}{{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x}) + \mathcal{O} ({{\boldsymbol{\mathbf{}}}}\kappa_L)$$ $\mathcal{O} ({{\boldsymbol{\mathbf{}}}}\kappa_L)$ represents higher order terms that tend to zero in the limit of small variance noise injections. \[appprop:accum\_noise\]

Recall that ${{\boldsymbol{\mathbf{}}}}{h}$ denotes the vanilla activations of the network, those we obtain with no noise injection. If we consider layer 0 to be the data ${{\boldsymbol{\mathbf{}}}}{x}$, we can define the accumulated noise on a layer by recursion. At layer 0 we have we have no accumulated noise: ${{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{0}=0$. Here we can invoke Taylor’s theorem on ${{\boldsymbol{\mathbf{}}}}{h}_1({{\boldsymbol{\mathbf{}}}}{x} + {{\boldsymbol{\mathbf{}}}}{\epsilon}_0)$ around the original input ${{\boldsymbol{\mathbf{}}}}x$. Namely, if we assume that all terms in Hessian of ${{\boldsymbol{\mathbf{}}}}{h}_1({{\boldsymbol{\mathbf{}}}}{x})$ are finite (i.e. $|\partial^2 {{\boldsymbol{\mathbf{}}}}{h}_1({{\boldsymbol{\mathbf{}}}}x)_i / \partial {{\boldsymbol{\mathbf{}}}}x_j {{\boldsymbol{\mathbf{}}}}x_k|<\infty\, \forall i,j,k$)

$${{\boldsymbol{\mathbf{}}}}{h}_1({{\boldsymbol{\mathbf{}}}}{x} + {{\boldsymbol{\mathbf{}}}}{\epsilon}_0) = {{\boldsymbol{\mathbf{}}}}{h}_1({{\boldsymbol{\mathbf{}}}}{x}) + \frac{\partial {{\boldsymbol{\mathbf{}}}}{h}_1}{\partial {{\boldsymbol{\mathbf{}}}}{x}}{{\boldsymbol{\mathbf{}}}}{\epsilon}_0 + \mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa}_0)
\label{eqapp:layer_taylor_expansion}$$

where once again, $\mathcal{O}({{\boldsymbol{\mathbf{}}}}\kappa_0)$ represents asymptotically dominated higher order terms as the variance of ${{\boldsymbol{\mathbf{}}}}{\epsilon}_0$ tends to 0. As such we have ${{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{1}= \frac{\partial {{\boldsymbol{\mathbf{}}}}{h}_1}{\partial {{\boldsymbol{\mathbf{}}}}{x}}{{\boldsymbol{\mathbf{}}}}{\epsilon}_0 + \mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa}_0)$.

Repeating this process for each layer, and assuming that all Hessians of the form $\nabla^2 {{\boldsymbol{\mathbf{}}}}{h}_l\rvert_{{{\boldsymbol{\mathbf{}}}}{h}_m({{\boldsymbol{\mathbf{}}}}{x})} \quad \forall m < l$ are finite, we obtain for a layer $k$: $${{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{k} = \sum_{i=0}^{k-1}\frac{\partial {{\boldsymbol{\mathbf{}}}}{h}_k}{\partial {{\boldsymbol{\mathbf{}}}}{h}_{i}}{{\boldsymbol{\mathbf{}}}}{\epsilon}_{i} + \mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa_k})
\label{eqapp:cumulative_noise}$$

where once again, $\mathcal{O}({{\boldsymbol{\mathbf{}}}}\kappa_k)$ represents asymptotically dominated higher order terms as the variance of all ${{\boldsymbol{\mathbf{}}}}{\epsilon}_k, k\in[0,\dots,k]$ tends to 0.

Let ${{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x}) \in \mathbb{R}^{N_L \times N_k}$ ($N_L$ the number of output neurons and $N_k$ the number of neurons in layer $k$) be defined element-wise as $({{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x}))_{i,j} = \frac{\partial ({{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x}))_i}{\partial ({{\boldsymbol{\mathbf{}}}}{h}_{k}({{\boldsymbol{\mathbf{}}}}{x}))_j}$. The accumulated noise at layer $L$ is: $${{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L} = \sum_{k=0}^{L-1}{{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{\epsilon}_{k} + \mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa}_L)
\label{eqapp:cumulative_noise_app}$$

For simplicity of notation we denote $\mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa}_L)$ as $\mathcal{O}({{\boldsymbol{\mathbf{}}}}{\gamma})$ completing the proof.

Explicit Regularisation Derivation {#app:exp_reg}
==================================

Consider an $L$ layer neural network experiencing GNIs ${{\boldsymbol{\mathbf{}}}}{\epsilon}_{i}$ at each layer $i \in [0, \dots, L - 1]$. We assume the the third derivative of $\mathcal{L}({{\boldsymbol{\mathbf{}}}}{x})$ w.r.t ${{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x})$ is finite, which is the case for most loss functions. Using the form of ${{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}$ in Proposition \[appprop:accum\_noise\]; we can marginalise out the injected noise ${{\boldsymbol{\mathbf{}}}}{\epsilon}$ to obtain an added regulariser: $$\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{\epsilon}} \left[ \Delta \mathcal{L}({{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}) \right] = \frac{1}{2}\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{x} \sim \mathcal{B}} \left[ \sum_{k=0}^{L-1}\sigma^2_k \mathrm{Tr}({{\boldsymbol{\mathbf{}}}}{J}^T_{k}({{\boldsymbol{\mathbf{}}}}{x})
    {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x}))  \right] +   O({{\boldsymbol{\mathbf{}}}}{\kappa}) 
    \label{eqn:deltaL}$$ where $\mathcal{L}({{\boldsymbol{\mathbf{}}}}{x})$ is the loss for a datapoint ${{\boldsymbol{\mathbf{}}}}{x}$, ${{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x})$ is $\nabla^2 \mathcal{L}\rvert_{{{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x})} \in \mathbb{R}^{N_L
\times N_L}$ ($N_L$ is the number of output neurons), and $O({{\boldsymbol{\mathbf{}}}}{\kappa})$ represents higher order terms in ${{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}$ that disappear in the limit of small variance noise. \[appprop:explicit\_reg\]

Recall that ${{\boldsymbol{\mathbf{}}}}{h}$ denotes the vanilla activations of the network, those we obtain with no noise injection. Let us *not* inject noise in the final, predictive, layer of our network such that the noise on this layer is accumulated from the noising of previous layers.

Our network loss can be defined in terms of our final layer activation ${{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x})$ for a datapoint ${{\boldsymbol{\mathbf{}}}}{x}$, using Taylor’s Theorem. Denoting $\nabla^2 \mathcal{L}\rvert_{{{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x})} \in \mathbb{R}^{N_L
\times N_L}$ ($N_L$ is the number of output neurons) as ${{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x})$ we have, and assuming that third derivatives of the loss w.r.t ${{\boldsymbol{\mathbf{}}}}{h}_L$ are finite ($|\partial^3 \mathcal{L}({{\boldsymbol{\mathbf{}}}}x) / \partial {{\boldsymbol{\mathbf{}}}}{h}_L({{\boldsymbol{\mathbf{}}}}x)_i {{\boldsymbol{\mathbf{}}}}{h}_L({{\boldsymbol{\mathbf{}}}}x)_j {{\boldsymbol{\mathbf{}}}}{h}_L({{\boldsymbol{\mathbf{}}}}x)_k|<\infty\, \forall i,j,k$) $$\begin{aligned}
    \mathcal{L}({{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x}) + {{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}({{\boldsymbol{\mathbf{}}}}{x})) &=  \mathcal{L}({{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x})) +  \nabla \mathcal{L}\rvert_{{{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x})}{{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}({{\boldsymbol{\mathbf{}}}}{x}) + \frac{1}{2}{{\boldsymbol{\mathbf{}}}}{\mathcal{E}}^T_{L}({{\boldsymbol{\mathbf{}}}}{x}) {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}) {{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}({{\boldsymbol{\mathbf{}}}}{x}) + \mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa})\nonumber \\
    &=  \mathcal{L}({{\boldsymbol{\mathbf{}}}}{x}) + \nabla \mathcal{L}\rvert_{{{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x})}{{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}({{\boldsymbol{\mathbf{}}}}{x}) + \frac{1}{2}{{\boldsymbol{\mathbf{}}}}{\mathcal{E}}^T_{L}({{\boldsymbol{\mathbf{}}}}{x}) {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}) {{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}({{\boldsymbol{\mathbf{}}}}{x}) + \mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa})
    \label{eqapp:2nd_order_taylor}\end{aligned}$$ Where $\mathcal{O}({{\boldsymbol{\mathbf{}}}}\kappa)$ represents asymptotically dominated higher order terms that go to zero in the limit of small variance noise injections. The second equality comes from the fact that the loss from the vanilla activation ${{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x})$ is equivalent to the loss resulting from the non-noised data ${{\boldsymbol{\mathbf{}}}}{x}$.

Given this definition of the accumulated noise, we can now return to equation . We have: $$\begin{aligned}
    \mathcal{L}({{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x}) + {{\boldsymbol{\mathbf{}}}}{\mathcal{E}}_{L}({{\boldsymbol{\mathbf{}}}}{x})) 
    =&  \mathcal{L}({{\boldsymbol{\mathbf{}}}}{x}) + \sum_{k=0}^{L-1}  \nabla \mathcal{L}\rvert_{{{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x})}{{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{\epsilon}_{k}({{\boldsymbol{\mathbf{}}}}{x}) \\ &+ \frac{1}{2}(\sum_{k=0}^{L-1} {{\boldsymbol{\mathbf{}}}}{\epsilon}^T_{k}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{J}^T_{k}({{\boldsymbol{\mathbf{}}}}{x})) {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}) (\sum_{k=0}^{L-1} {{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{\epsilon}_{k}({{\boldsymbol{\mathbf{}}}}{x}))  + \mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa}) + \mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa}_L) \\ 
    =&  \mathcal{L}({{\boldsymbol{\mathbf{}}}}{x}) + \sum_{k=0}^{L-1}  \nabla \mathcal{L}\rvert_{{{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x})}{{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{\epsilon}_{k}({{\boldsymbol{\mathbf{}}}}{x}) + \frac{1}{2}\sum_{k=0}^{L-1} {{\boldsymbol{\mathbf{}}}}{\epsilon}^T_{k}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{J}^T_{k}({{\boldsymbol{\mathbf{}}}}{x}) {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{\epsilon}_{k}({{\boldsymbol{\mathbf{}}}}{x}) \nonumber \\
    &+ \frac{1}{2}(\sum_{k=0}^{L-1} {{\boldsymbol{\mathbf{}}}}{\epsilon}^T_{k}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{J}^T_{k}({{\boldsymbol{\mathbf{}}}}{x}) {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}) \sum_{j \neq k}^{L-1} {{\boldsymbol{\mathbf{}}}}{J}_{j}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{\epsilon}_{j}({{\boldsymbol{\mathbf{}}}}{x})) + \mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa}) + \mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa}_L)
    \label{eqapp:full_taylor_loss}\end{aligned}$$ Because both $\mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa})$ and $\mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa}_L)$ are higher order terms that tend to 0 in the limit of small variance noise injections, we write their sum for notation simplicity solely as $\mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa})$. Because the injected noise at each layer is isotropic, symmetric centered on 0, and is i.i.d we have many of the multiplicative terms in the first and second order terms of the Taylor series can be ignored: $$\begin{aligned}
    &\int {{\boldsymbol{\mathbf{}}}}{\epsilon}_{k}({{\boldsymbol{\mathbf{}}}}{x})p({{\boldsymbol{\mathbf{}}}}{\epsilon}_k)d_{{{\boldsymbol{\mathbf{}}}}{\epsilon}_k} = 0 
     \label{eqapp:marg}
    \\
     &\int {{\boldsymbol{\mathbf{}}}}{\epsilon}^T_{k}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{J}^T_{k}({{\boldsymbol{\mathbf{}}}}{x}) {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}) 
    {{\boldsymbol{\mathbf{}}}}{J}_{j}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{\epsilon}_{j}({{\boldsymbol{\mathbf{}}}}{x})p({{\boldsymbol{\mathbf{}}}}{\epsilon}_k)p({{\boldsymbol{\mathbf{}}}}{\epsilon}_j)d_{{{\boldsymbol{\mathbf{}}}}{\epsilon}_k}d_{{{\boldsymbol{\mathbf{}}}}{\epsilon}_j} = 0, \ \forall j \neq k
    \label{eqapp:marg_Hess_offdiag}
    \\
    &\int {{\boldsymbol{\mathbf{}}}}{\epsilon}^T_{k}({{\boldsymbol{\mathbf{}}}}{x}) {{\boldsymbol{\mathbf{}}}}{J}^T_{k}({{\boldsymbol{\mathbf{}}}}{x}) {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}) 
    {{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x})
    {{\boldsymbol{\mathbf{}}}}{\epsilon}_{k}({{\boldsymbol{\mathbf{}}}}{x})p({{\boldsymbol{\mathbf{}}}}{\epsilon}_{k})p({{\boldsymbol{\mathbf{}}}}{\epsilon}_k)d_{{{\boldsymbol{\mathbf{}}}}{\epsilon}_k} = \sigma^2_k \mathrm{Tr}({{\boldsymbol{\mathbf{}}}}{J}^T_{k}({{\boldsymbol{\mathbf{}}}}{x})
    {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x}))
    \label{eqapp:marg_Hess_diag}\end{aligned}$$ As such: $$\begin{aligned}
\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{\epsilon}}\mathcal{\tilde{L}}({{\boldsymbol{\mathbf{}}}}{x}) &= 
\mathcal{L}({{\boldsymbol{\mathbf{}}}}{x}) + \frac{1}{2}\sum_{k=0}^{L-1} \sigma^2_k \mathrm{Tr}({{\boldsymbol{\mathbf{}}}}{J}^T_{k}({{\boldsymbol{\mathbf{}}}}{x})
    {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x})) + \mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa}) \\
\label{eqapp:taylor_marginal} 
\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{\epsilon}} \Delta \mathcal{L}({{\boldsymbol{\mathbf{}}}}{x}) &= \frac{1}{2}\sum_{k=0}^L \sigma^2_k \mathrm{Tr}({{\boldsymbol{\mathbf{}}}}{J}^T_{k}({{\boldsymbol{\mathbf{}}}}{x})
    {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x})) + \mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa}) \end{aligned}$$ For ${{\boldsymbol{\mathbf{}}}}{x}$ drawn from some batch $\mathcal{B}$ of size $B$, this becomes: $$\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{\epsilon}} \mathcal{\tilde{L}}_{SGD} = 
\frac{1}{B}\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{x} \sim \mathcal{B}} [\mathcal{L}({{\boldsymbol{\mathbf{}}}}{x}) + \mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{\epsilon}} [\Delta \mathcal{L}({{\boldsymbol{\mathbf{}}}}{x})]] + \mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa})
\label{eqapp:full_eq_hessian_app}$$ We denote the added regulariser, *the explicit regulariser*, and have: $$R = \frac{1}{2B}\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{x} \sim \mathcal{B}} \left[ \sum_{k=0}^L \sigma^2_k \mathrm{Tr}({{\boldsymbol{\mathbf{}}}}{J}^T_{k}({{\boldsymbol{\mathbf{}}}}{x})
    {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x})) \right] + \mathcal{O}({{\boldsymbol{\mathbf{}}}}{\kappa})$$

Explicit Regularisation and Gauss-Newton Hessians {#app:exp_reg_gn}
-------------------------------------------------

Let us now consider the meaning of the terms contained in $R$. First consider the Hessian of the loss with respect to the activations of a given layer $k$: $\nabla^2 \mathcal{L}\rvert_{{{\boldsymbol{\mathbf{}}}}{h}_{k}({{\boldsymbol{\mathbf{}}}}{x})} ={{\boldsymbol{\mathbf{}}}}{H}_{k}({{\boldsymbol{\mathbf{}}}}{x})$. We can decompose this Hessian with respect to the activations of the final layer ${{\boldsymbol{\mathbf{}}}}{h}_L$ as follows:

$${{\boldsymbol{\mathbf{}}}}{H}^{i,j}_{k} = \sum_o \frac{\partial \mathcal{L}}{\partial h_L^o}\frac{\partial^2 h_L^o}{\partial h_k^i\partial h_k^j}  + \sum_{o,l} \frac{\partial h_L^o}{\partial h_k^i}\frac{\partial^2 \mathcal{L}}{\partial h_L^o\partial h_L^l}\frac{\partial h_L^o}{\partial h_k^j}$$

Here we assume that our final layer activations are the logits of our final layer ${{\boldsymbol{\mathbf{}}}}{h}_L$, before any non-linearity is applied (eg. $\mathrm{softmax}$ in the case of classification). Due to the extended envelope property, we assume that ${{\boldsymbol{\mathbf{}}}}{h}_L$ is approximately piece-wise linear wrt to any of the other previous activations, an assumption we made in our prior derivations. As such the second order derivatives in the first sum will be approximately 0. Ignoring the sums in that first term which are not guaranteed to be PSD, results in the the Gauss-Newton approximation, guaranteed to be PSD if ${{\boldsymbol{\mathbf{}}}}{H}_{L}$ is PSD [@Botev2017; @Sagun2018; @Wei2020]:

$${{\boldsymbol{\mathbf{}}}}{H}_{k} \approx \operatorname{\mathbb{E}}[{{\boldsymbol{\mathbf{}}}}{J}^T_{k}{{\boldsymbol{\mathbf{}}}}{H}_{L} {{\boldsymbol{\mathbf{}}}}{J}_{k}]_n
\label{eqapp:gauss-newton-hessian}$$

Regularisation in Regression Models and Autoencoders {#app:exp_reg_regression}
----------------------------------------------------

In the case of regression the most commonly used loss is the mean-square error. For a point $n$ we have:

$$\mathcal{L}({{\boldsymbol{\mathbf{}}}}{x}) = ({{\boldsymbol{\mathbf{}}}}{y}_n - {{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x}))^2
\label{eqapp:mse}$$

In this case, ${{\boldsymbol{\mathbf{}}}}{H}_{L, n}$ is identity. For ${{\boldsymbol{\mathbf{}}}}{x}$ drawn from some batch $\mathcal{B}$ of size $B$, and assuming low variance noise such that we ignore higher order terms in the explicit regulariser of Theorem \[prop:explicit\_reg\], our loss is approximately equivalent to:

$$\begin{aligned}
\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{\epsilon}} \mathcal{\tilde{L}}_{SGD} &\approx  \frac{1}{B}\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{x} \sim \mathcal{B}} [\mathcal{L}({{\boldsymbol{\mathbf{}}}}{x}) + \frac{1}{2}\sum_{k=0}^L \sigma_k^2(\mathrm{Tr}({{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x})^T{{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x})))] \\
&=\frac{1}{B}\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{x} \sim \mathcal{B}} [\mathcal{L}({{\boldsymbol{\mathbf{}}}}{x}) + \frac{1}{2}\sum_{k=0}^L \sigma_k^2(\|{{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x})\|^2_F)]
\label{eqapp:mse_reg}\end{aligned}$$

This added term corresponds to the trace of the covariance matrix of the outputs ${{\boldsymbol{\mathbf{}}}}{h}_L$ given an input ${{\boldsymbol{\mathbf{}}}}{h}_{k}$. As such we are penalising the sum of output variances of the approximator; we are penalising the sensitivity of outputs to perturbations in layer $k$ [@Webb1994; @Bishop1995]. The greater the variance of the injected noise, the greater the penalisation. As such, the added regularisation term, which is guaranteed to be PSD, is:

$$R \approx \frac{1}{2B}\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{x} \sim \mathcal{B}} \sum_{k=0}^L \sigma^2_k (\|{{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x})\|^2_F)
    \label{eqapp:mse_reg_term}$$

Given the extended envelope property cited above, because our functions are at *most* linear, we can bound our regularisers using the Jacobian of an equivalent linear network:

$$\sum_{k=0}^L \sigma^2_k(\|{{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x})\|^2) < \sum_{k=0}^L \sigma^2_k(\|{{\boldsymbol{\mathbf{}}}}{J}^{\mathrm{linear}}_{k}({{\boldsymbol{\mathbf{}}}}{x})\|^2) = \sum_{k=0}^L\sigma^2_k(\|{{\boldsymbol{\mathbf{}}}}{W}_L \dots {{\boldsymbol{\mathbf{}}}}{W}_k\|^2)
\label{eqapp:mse_reg_bound}$$

Where ${{\boldsymbol{\mathbf{}}}}{J}^{\mathrm{linear}}_{k}({{\boldsymbol{\mathbf{}}}}{x})$ is the gradient evaluated with no non-linearities in our network. This upper bound is reminiscent of $rank-k$ ridge regression, but here we penalise each sub-network in our network [@Kunin1999].

Also note that the regression setting is directly translatable to Auto-Encoders, where the labels are the input data ${{\boldsymbol{\mathbf{}}}}{x}$. We have:

$$\mathcal{L}({{\boldsymbol{\mathbf{}}}}{x}) = ({{\boldsymbol{\mathbf{}}}}{x} - {{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x}))^2
\label{eqapp:mse_ae}$$

Regularisation in Classifiers {#app:exp_reg_classification}
-----------------------------

In the case of classification, we consider the cross-entropy loss. Recall that we consider our network outputs ${{\boldsymbol{\mathbf{}}}}{h}_L$ to be the pre-$\mathrm{softmax}$ of logits of the final layer ${{\boldsymbol{\mathbf{}}}}{L}$. We denote ${{\boldsymbol{\mathbf{}}}}{p}({{\boldsymbol{\mathbf{}}}}{x})=\mathrm{softmax}({{\boldsymbol{\mathbf{}}}}{h}_L({{\boldsymbol{\mathbf{}}}}{x}))$. The loss is thus:

$$\mathcal{L}({{\boldsymbol{\mathbf{}}}}{x}) = - \sum_{c=0}^M {{\boldsymbol{\mathbf{}}}}{y}_{n,c} \log (\mathrm{softmax}({{\boldsymbol{\mathbf{}}}}{h}_{L}({{\boldsymbol{\mathbf{}}}}{x}))_c)
\label{eqapp:ce}$$

where $c$ indexes over the $M$ possible classes of the classification problem. The hessian ${{\boldsymbol{\mathbf{}}}}{H}_L$ in this case is easy to compute and has the form:

$${{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x})_{i,j} = 
\begin{cases}
{{\boldsymbol{\mathbf{}}}}{p}({{\boldsymbol{\mathbf{}}}}{x})_i(1 - {{\boldsymbol{\mathbf{}}}}{p}({{\boldsymbol{\mathbf{}}}}{x})_j) & i = j \\
-{{\boldsymbol{\mathbf{}}}}{p}({{\boldsymbol{\mathbf{}}}}{x})_i{{\boldsymbol{\mathbf{}}}}{p}({{\boldsymbol{\mathbf{}}}}{x})_j & i \neq j \\
\end{cases}$$

As @Wei2020, @Sagun2018, and [@LeCun1998] show, this Hessian is PSD, meaning that $\mathrm{Tr}({{\boldsymbol{\mathbf{}}}}{J}_{k}{{\boldsymbol{\mathbf{}}}}{H}_{L} {{\boldsymbol{\mathbf{}}}}{J}^T_{k})$ will be positive, fulfilling the criteria for a valid regulariser. For ${{\boldsymbol{\mathbf{}}}}{x}$ drawn from some batch $\mathcal{B}$ of size $B$, the cross-entropy loss under small variance noise injections such that we ignore the higher order terms in the explicit regulariser of Theorem \[prop:explicit\_reg\], is approximately equal to:

$$\begin{aligned}
\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{\epsilon}_0 \dots {{\boldsymbol{\mathbf{}}}}{\epsilon}_{L-1}} \mathcal{\tilde{L}}_{SGD} &\approx    \frac{1}{B}\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{x} \sim \mathcal{B}} [\mathcal{L}({{\boldsymbol{\mathbf{}}}}{x}) + \frac{1}{2}\sum_{k=0}^L \sigma_k^2 \sum_{i,j}(  {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}) \circ {{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{{{\boldsymbol{\mathbf{}}}}{J}}^T_{k}({{\boldsymbol{\mathbf{}}}}{x}))_{i,j}]\\
&= \frac{1}{B}\mathbb{E}_{{{\boldsymbol{\mathbf{}}}}{x} \sim \mathcal{B}} [\mathcal{L}({{\boldsymbol{\mathbf{}}}}{x}) + \frac{1}{2}\sum_{k=0}^L \sigma_k^2 \sum_{i,j}(\mathrm{diag}({{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}))^T{{\boldsymbol{\mathbf{}}}}{J}^2_{k}({{\boldsymbol{\mathbf{}}}}{x}) )_{i,j} \\ &+ \frac{1}{2}\sum_{k=0}^L \sigma_k^2 \sum_{\forall i,j \ i \neq j}({{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}) \circ {{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x}){{\boldsymbol{\mathbf{}}}}{{{\boldsymbol{\mathbf{}}}}{J}}^T_{k}({{\boldsymbol{\mathbf{}}}}{x}))_{i,j}]
\label{eqapp:ce_reg}\end{aligned}$$

$\mathrm{diag}( {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}))^T$ is the row vector of the diagonal of ${{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x})$. The first equality us due to the fact that ${{\boldsymbol{\mathbf{}}}}{H}_L$ is symmetric and is due to the commutative properties of the trace operator. The final equality is simply the decomposition of the sum of the matrix product into diagonal and off-diagonal elements. Components of the form $ \sum_{i,j}(\mathrm{diag}( {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}))^T){{\boldsymbol{\mathbf{}}}}{J}^2_{k}({{\boldsymbol{\mathbf{}}}}{x}))_{i,j} $ are minimimised when $\|{{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x})\|^2_F$ is minimised. This corresponds to a penalisation of each of the subnetworks starting at layer $k$ and ending at layer $L$, weighted by the diagonal of the cross-entropy of the predicted class probability, ${{\boldsymbol{\mathbf{}}}}{H}_L$. See Appendix \[app:covariance\] for an empirical explanation as to why the off-diagonal elements of ${{\boldsymbol{\mathbf{}}}}{J}_k{{\boldsymbol{\mathbf{}}}}{J}^T_k$ are approximately 0.

Ignoring these off-diagonal terms, the added PSD regularisation term we obtain is:

$$R \approx \frac{1}{2B}\sum_{k=0}^L \sigma_k^2 \sum_{i,j}(\mathrm{diag}( {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}))^T{{{\boldsymbol{\mathbf{}}}}{J}}^2_{k}({{\boldsymbol{\mathbf{}}}}{x}))_{i,j}
\label{eqapp:ce_reg_term}$$

Given the extended envelope property cited above, because our functions are at *most* linear, we can bound our regularisers using the Jacobian of an equivalent linear network:

$$\sum_{k=0}^L \sigma_k^2 \sum_{i,j}(\mathrm{diag}( {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}))^T{{\boldsymbol{\mathbf{}}}}{J}_{k}({{\boldsymbol{\mathbf{}}}}{x})^2 )_{i,j} <   \sum_{k=0}^L\sigma^2_k\sum_{i,j}(\mathrm{diag}( {{\boldsymbol{\mathbf{}}}}{H}_{L}({{\boldsymbol{\mathbf{}}}}{x}))^T({{\boldsymbol{\mathbf{}}}}{W}_L \dots {{\boldsymbol{\mathbf{}}}}{W}_k)^2 ))_{i,j}
\label{eqapp:ce_reg_bound}$$

${{\boldsymbol{\mathbf{}}}}{J}_k{{\boldsymbol{\mathbf{}}}}{J}^T_k$ as a Covariance Matrix {#app:covariance}
-----------------------------------------------------------------------------------------

${{\boldsymbol{\mathbf{}}}}{J}_k{{\boldsymbol{\mathbf{}}}}{J}^T_k$, can be interpreted as the covariance of the network outputs given noise in layer $k$ [@Bishop1995; @Webb1994]. For relatively shallow networks, the off-diagonal elements of this metric *within* the network are likely to be small and ${{\boldsymbol{\mathbf{}}}}{J}_k{{\boldsymbol{\mathbf{}}}}{J}^T_k$ can be approximated by ${{\boldsymbol{\mathbf{}}}}{J}^2_k$ [@Poole2016; @Hauser2017; @Farquhar2020; @Aleksziev]. This is also true, though to a lesser extent, for the data layer $k$, i.e the extraction of independent features occurs in the very first layer of the network [@Poole2016]. See Figure \[fig:covariance\_examples\] for a demonstration that the off-diagonal elements of ${{\boldsymbol{\mathbf{}}}}{J}_k^T{{\boldsymbol{\mathbf{}}}}{J}_k$, are negligible for smaller networks.

In Figure \[fig:output\_uncertainty\] we show that models trained with Gaussian noise and the derived explicit regularisers induce networks with much smaller norms of the Jacobian.

\

Fourier Domain and Model Capacity {#app:fourier_capacity}
---------------------------------

A model’s capacity is a measure of a model’s ‘complexity’. Formally if we have a model class $\mathcal{H}$, then the capacity assigns a non-negative number to each hypothesis in the model class $\mathcal{M} : \{\mathcal{H},\mathcal{D}_{train}\} \to \mathbb{R}^+$, where $\mathcal{D}_{train}$ is the training set and a lower capacity is an indicator of better model generalisation [@Neyshabur2017]. Generally, deeper and narrower networks induce large capacity models that are likely to overfit and generalise poorly [@Zhang2019]. The network Jacobian’s spectral norm, Frobenius norm, and the spectral norm of the product of weights ($\prod_{{{\boldsymbol{\mathbf{}}}}{W}_k \in \theta} \|{{\boldsymbol{\mathbf{}}}}{W}_k\|_F$) are good approximators of model capacity and are clearly linked to the explicit regularisers we have derived in equations and [@Guo2017; @Neyshabur2017; @Neyshabur].

As we have shown, the Frobenius norm of the network Jacobian corresponds to a norm in Sobolev space which is a measure of a network’s high-frequency components in the Fourier domain. From this we offer the first theoretical results on why norms of the Jacobian are a good measure of model capacity: as low-frequency functions correspond to smoother functions that are less prone to overfitting, a smaller norm of the Jacobian is thus a measure of a smoother ‘less complex’ model.

Classification Margins {#app:margin}
======================

Our aim is to establish a classification margin $M$ that is radius of the largest metric ball centered on a point ${{\boldsymbol{\mathbf{}}}}{x}$ to which a classifier assigns the same class. As mentioned before @Sokolic2017 and @Jakubovitz2018 establish a stricter margin where $R$ is radius of the largest metric ball centered on a point ${{\boldsymbol{\mathbf{}}}}{x}$ to which a classifier assigns ${{\boldsymbol{\mathbf{}}}}{y}$, the true label. Though the definitions are different, their proofs still hold for our definition.

By Lemma 1 of @Jakubovitz2018, the first order approximation for the distance between an input ${{\boldsymbol{\mathbf{}}}}{x}$, assigned class $A$ by the network, and a perturbed input classified to the boundary hyper-surface separating the classes $A$ and $k$ for an $L_2$ distance metric is given by:

$$d = \frac{({{\boldsymbol{\mathbf{}}}}{h}^A_{L}({{\boldsymbol{\mathbf{}}}}{x}) - {{\boldsymbol{\mathbf{}}}}{h}^{k}_{L}({{\boldsymbol{\mathbf{}}}}{x}) )}{\|{{\boldsymbol{\mathbf{}}}}{J}^A_{0}({{\boldsymbol{\mathbf{}}}}{x}) - {{\boldsymbol{\mathbf{}}}}{J}^k_{0}({{\boldsymbol{\mathbf{}}}}{x})\|_2}$$

we have ${{\boldsymbol{\mathbf{}}}}{h}^A_{L}({{\boldsymbol{\mathbf{}}}}{x}) \geq {{\boldsymbol{\mathbf{}}}}{h}^{k}_{L}({{\boldsymbol{\mathbf{}}}}{x})$, where ${{\boldsymbol{\mathbf{}}}}{h}^A_{L}({{\boldsymbol{\mathbf{}}}}{x})$ is the $L^{th}$ layer activation (pre-softmax) associated with the most probable class $A$ assigned to ${{\boldsymbol{\mathbf{}}}}{x}$, and ${{\boldsymbol{\mathbf{}}}}{h}^{k}_{L}({{\boldsymbol{\mathbf{}}}}{x})$ is the $L^{th}$ layer activation associated with some other class $k$. ${{\boldsymbol{\mathbf{}}}}{J}^k_{0}({{\boldsymbol{\mathbf{}}}}{x})$ is then the Jacobian between neuron $k$ of the $L^{th}$ layer and the input ${{\boldsymbol{\mathbf{}}}}{x}$.

The $L_2$ norm of the minimal perturbation required to alter the network’s classification is then given by:

$$M({{\boldsymbol{\mathbf{}}}}{x}) = \frac{({{\boldsymbol{\mathbf{}}}}{h}^A_{L}({{\boldsymbol{\mathbf{}}}}{x}) - {{\boldsymbol{\mathbf{}}}}{h}^{B}_{L}({{\boldsymbol{\mathbf{}}}}{x}) )}{\|{{\boldsymbol{\mathbf{}}}}{J}^A_{0}({{\boldsymbol{\mathbf{}}}}{x}) - {{\boldsymbol{\mathbf{}}}}{J}^B_{0}({{\boldsymbol{\mathbf{}}}}{x})\|_2}$$

where $B$ is the second most probable class according to network outputs.

We know that:

$$\begin{aligned}
     M^2({{\boldsymbol{\mathbf{}}}}{x}) &= \frac{({{\boldsymbol{\mathbf{}}}}{h}^A_{L}({{\boldsymbol{\mathbf{}}}}{x}) - {{\boldsymbol{\mathbf{}}}}{h}^{B}_{L}({{\boldsymbol{\mathbf{}}}}{x}))^2}{\|{{\boldsymbol{\mathbf{}}}}{J}^A_{0}({{\boldsymbol{\mathbf{}}}}{x}) - \nabla_{{{\boldsymbol{\mathbf{}}}}{x}} {{\boldsymbol{\mathbf{}}}}{J}^B_{0}({{\boldsymbol{\mathbf{}}}}{x})\|^2_2} \\
    &= \frac{({{\boldsymbol{\mathbf{}}}}{h}^A_{L}({{\boldsymbol{\mathbf{}}}}{x}) - {{\boldsymbol{\mathbf{}}}}{h}^{B}_{L}({{\boldsymbol{\mathbf{}}}}{x}))^2}{\|{{\boldsymbol{\mathbf{}}}}{J}^A_{0}({{\boldsymbol{\mathbf{}}}}{x})\|^2_2 - 2{{\boldsymbol{\mathbf{}}}}{J}^A_{0}({{\boldsymbol{\mathbf{}}}}{x})({{\boldsymbol{\mathbf{}}}}{J}^B_{0}({{\boldsymbol{\mathbf{}}}}{x}))^T + \|{{\boldsymbol{\mathbf{}}}}{J}^B_{0}({{\boldsymbol{\mathbf{}}}}{x})\|^2_2} \\
    &\geq \frac{({{\boldsymbol{\mathbf{}}}}{h}^A_{L}({{\boldsymbol{\mathbf{}}}}{x}) - {{\boldsymbol{\mathbf{}}}}{h}^{B}_{L}({{\boldsymbol{\mathbf{}}}}{x}))^2}{2(\|{{\boldsymbol{\mathbf{}}}}{J}^A_{0}({{\boldsymbol{\mathbf{}}}}{x})\|^2_2 +  \|{{\boldsymbol{\mathbf{}}}}{J}^B_{0}({{\boldsymbol{\mathbf{}}}}{x})\|^2_2)}\end{aligned}$$

Since $\sum_{k} \|{{\boldsymbol{\mathbf{}}}}{J}^k_{0}({{\boldsymbol{\mathbf{}}}}{x})\|^2_2 =  \|{{\boldsymbol{\mathbf{}}}}{J}_{0}({{\boldsymbol{\mathbf{}}}}{x})\|_F$, where $k$ indexes over classes, we have the upper bound:

$$(\|{{\boldsymbol{\mathbf{}}}}{J}^A_{0}({{\boldsymbol{\mathbf{}}}}{x})\|^2_2 +  \|{{\boldsymbol{\mathbf{}}}}{J}^B_{0}({{\boldsymbol{\mathbf{}}}}{x})\|^2_2) < \|{{\boldsymbol{\mathbf{}}}}{J}_{0}({{\boldsymbol{\mathbf{}}}}{x})\|^2_F$$

As such:

$$M^2({{\boldsymbol{\mathbf{}}}}{x}) \geq \frac{({{\boldsymbol{\mathbf{}}}}{h}^k_{L}({{\boldsymbol{\mathbf{}}}}{x}) - {{\boldsymbol{\mathbf{}}}}{h}^{B}_{L}({{\boldsymbol{\mathbf{}}}}{x}))^2}{2\|{{\boldsymbol{\mathbf{}}}}{J}_{0}({{\boldsymbol{\mathbf{}}}}{x})\|^2_F}$$

From which we have:

$$M({{\boldsymbol{\mathbf{}}}}{x}) \geq \frac{({{\boldsymbol{\mathbf{}}}}{h}^A_{L}({{\boldsymbol{\mathbf{}}}}{x}) - {{\boldsymbol{\mathbf{}}}}{h}^{B}_{L}({{\boldsymbol{\mathbf{}}}}{x}))}{\sqrt{2}\|{{\boldsymbol{\mathbf{}}}}{J}_{0}({{\boldsymbol{\mathbf{}}}}{x})\|_F}$$

As such:

$$M({{\boldsymbol{\mathbf{}}}}{x}) = \frac{({{\boldsymbol{\mathbf{}}}}{h}^A_{L}({{\boldsymbol{\mathbf{}}}}{x}) - {{\boldsymbol{\mathbf{}}}}{h}^{B}_{L}({{\boldsymbol{\mathbf{}}}}{x}))}{\sqrt{2}\|{{\boldsymbol{\mathbf{}}}}{J}_{0}({{\boldsymbol{\mathbf{}}}}{x})\|_F}$$

is a valid classification margin.

Finally we can state that for networks that use activation functions that obey the extended envelope property:

$$M({{\boldsymbol{\mathbf{}}}}{x}) \geq \frac{({{\boldsymbol{\mathbf{}}}}{h}^A_{L}({{\boldsymbol{\mathbf{}}}}{x}) - {{\boldsymbol{\mathbf{}}}}{h}^{B}_{L}({{\boldsymbol{\mathbf{}}}}{x}))}{\sqrt{2}\|{{\boldsymbol{\mathbf{}}}}{J}_{0}({{\boldsymbol{\mathbf{}}}}{x})\|_F}$$

The second bound comes from the fact that our functions are at *most* linear (see Appendix \[app:eep\]), and that we can bound our regularisers using the Jacobian of an equivalent linear network, which is the product of weights.

Measuring Calibration {#app:calibration}
=====================

A neural network classifier gives a prediction $\hat{y}({{\boldsymbol{\mathbf{}}}}{x})$ with confidence $\hat{p}({{\boldsymbol{\mathbf{}}}}{x})$ (the probability attributed to that prediction) for a datapoint ${{\boldsymbol{\mathbf{}}}}{x}$. Perfect calibration consists of being as likely to be correct as you are confident: $$p(\hat{y}=y|\hat{p}=r)=r, \quad \forall r\in[0,1]$$ To see how closely a model approaches perfect calibration, we plot reliability diagrams [@Guo2017; @Niculescu-Mizil2005], which show the accuracy of a model as a function of its confidence over $M$ bins $B_m$.

$$\begin{aligned}
\mathrm{acc}(B_m) &=\frac{1}{|B_m|}\sum_{i\in B_m}\mathbf{1}(\hat{y}_i = y_i) \\
\mathrm{conf}(B_m) &=\frac{1}{|B_m|}\sum_{i\in B_m}\hat{p}_i\end{aligned}$$

We also calculate the Expected Calibration Error (ECE) [@Naeini2015], the mean difference between the confidence and accuracy over bins: $$\mathrm{ECE} = \sum_{m=1}^M \frac{|B_m|}{N}|\mathrm{acc}(B_m) - \mathrm{conf}(B_m)|
        \label{eqapp:ece}$$ However, note that ECE only measures calibration, not refinement. For example, if we have a balanced test set one can trivially obtain ECE $\approx 0$ by sampling predictions from a uniform distribution over classes while having very low accuracy.

Results Explicit Regularisation {#app:exp_reg_results}
===============================

\

\
\

Classification Margins and Calibration Results {#app:exp_reg_results_margins}
==============================================

\

\
\
\

Network Hyperparameters {#app:hyperparams}
=======================

All networks were trained using stochastic gradient descent with a learning rate of 0.001 and a batch size of 512.

All MLP networks, unless specified otherwise, are 2 hidden layer networks with 512 units per layer.

All convolutional (CONV) networks are 2 hidden layer networks. The first layer has 32 filters, a kernel size of 4, and a stride length of 2. The second layer has 128 filters, a kernel size of 4, and a stride length of 2. The final output layer is a dense layer.

[^1]: In our notation ${{\boldsymbol{\mathbf{}}}}{v}^2$ is a vector composed of the elementwise square components of ${{\boldsymbol{\mathbf{}}}}{v}$. Similarly, for matrices ${{\boldsymbol{\mathbf{}}}}{M}^2$ is elementwise square.
