\section{Introduction}
\label{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 or 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 \citep{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 \citep{Bishop1995, Cohen2019, Webb1994}, there have been few studies on understanding the benefits of methods that inject noise \textit{throughout} a network.
Here we study the \textit{explicit} regularisation of such injections, which is a positive term added to the loss function obtained when we marginalise out the noise we have injected.

Concretely our contributions are: 

\begin{itemize}
    \item We derive an analytic form for an explicit regulariser that explains most of GNIs' regularising effect.
    \item We show that this regulariser penalises networks that learn functions with high-frequency content in the Fourier domain and most heavily regularises neural network layers that are closer to the output. See Figure \ref{fig:illustration} for an illustration.
    \item Finally, we show analytically and empirically that this regularisation induces larger classification margins and better calibration of models. 
\end{itemize}

   \begin{figure}[t!]

    \centering
    \includegraphics[width=0.9\textwidth]{figures/frequencydiag2.pdf}
    \caption{
    Here we illustrate the effect of GNIs injected throughout a network's activations. Each coloured dot represents a neuron's activations.
    We add GNIs, represented as circles, to each layer's activations bar the output layer. 
    GNIs induce a network for which each layer learns a progressively lower frequency function, represented as a sinusoid matching in colour to its corresponding layer. 
    }
    \vspace{-1em}
    \label{fig:illustration}
\end{figure}







\section{Background}

\subsection{Gaussian Noise Injections}
\label{sec:sgd_under_injections}

Training a neural network involves optimising network parameters to maximise the marginal likelihood of a set of labels given features via gradient descent.
With a training dataset $\mathcal{D}$ composed of $N$ data-label pairs of the form $(\v{x}, \v{y}) \ \v{x} \in \mathbb{R}^d, \v{y} \in \mathbb{R}^m$ and a feed-forward neural network with $M$ parameters divided into $L$ layers: $\v{\theta} = \{\v{W}_1,...,\v{W}_L\}$, $\v{\theta} \in \mathbb{R}^M$,
our objective is to minimise the expected negative log likelihood of labels $\v{y}$  given data $\v{x}$,  $- \log p_{\v{\theta}}(\v{y}|\v{x})$ , and find the optimal set of parameters  $\v{\theta}^*$ satisfying: 
\begin{align}
    \v{\theta}_* = \argmin_{\v{\theta}} \mathcal{L}( \mathcal{D}; \v{\theta}),
    \qquad
    \mathcal{L}( \mathcal{D}; \v{\theta}):=- \mathbb{E}_{\v{x},\v{y} \sim \mathcal{D}}\left[ \log p_{\v{\theta}}(\v{y}|\v{x})\right]\,.
    \label{eq:marginal_likelihood}
\end{align}
Under stochastic optimisation algorithms, such as Stochastic Gradient Descent (SGD), we  estimate $\mathcal{L}$ by sampling a mini-batch of data-label pairs $\mathcal{B} \subset \mathcal{D}$. 
\begin{equation}
\mathcal{L}(\mathcal{B}; \v{\theta}) = -\mathbb{E}_{\v{x},\v{y} \sim \mathcal{B}} \log p_{\v{\theta}}(\v{y}|\v{x}) \approx \mathcal{L}( \mathcal{D}; \v{\theta}).
\end{equation}
Consider an $L$ layer network with no noise injections and a non-linearity $\phi$ at each layer.
We obtain the activations ${\v{h}} = \{{\v{h}}_0, ... , {\v{h}}_{L} \}$, where $\v{h}_{0}=\v{x}$ is the input data \textit{before} any noise is injected. For a network consisting of dense layers (a.k.a.\ a multi-layer perceptron: MLP) we have that: 
\begin{align}
\v{h}_{k}(\v{x})=
\phi(\v{W}_k \v{h}_{k-1}(\v{x})).
\label{eq:nonoise_acts}
\end{align}
What happens to these activations when we inject noise? 
First, let $\v{\epsilon}$ be the set of noise injections at each layer:  $\v{\epsilon} = \{{\v{\epsilon}}_0, ... , {\v{\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 as $\widehat{\v{h}}_{k}$ and the subsequently noised value as $\widetilde{\v{h}}_{k}$: 
\begin{align}
\widehat{\v{h}}_{k}(\v{x})= \phi\left(\v{W}_{k}\widetilde{\v{h}}_{k-1}(\v{x})\right)\,, \qquad  
\widetilde{\v{h}}_{k}(\v{x}) = \widehat{\v{h}}_{k}(\v{x}) \circ \v{\epsilon}_{k}\,, 
\label{eq:noise_recursion}
\end{align}
where $\circ$ is some element-wise operation. 
We can, for example, add or multiply Gaussian noise to each hidden layer unit. 
In the additive case, we obtain: 
\begin{align}
\widetilde{\v{h}}_k(\v{x}) &= \widehat{\v{h}}_k(\v{x}) + \v{\epsilon}_k, 
\qquad \v{\epsilon}_k \sim \mathcal{N}(0,\sigma_k^2\v{I}). \label{eq:noise_add}
\end{align}
The multiplicative case can be rewritten as an activation-scaled addition:
\begin{align}
\widetilde{\v{h}}_k(\v{x}) &= \widehat{\v{h}}_k(\v{x}) + \v{\epsilon}_k, 
\qquad \v{\epsilon}_k \sim \mathcal{N}\left(0,\widehat{\v{h}}_k^2(\v{x})\sigma_k^2\v{I}\right).
\label{eq:noise_mult}
\end{align}
Here we focus our analysis on noise \textit{additions}, but through equation \eqref{eq:noise_mult} we can translate our results to the multiplicative case.


















\subsection{Sobolev Spaces}
\label{sec:sobolev}
To define a Sobolev Space we use the generalisation of the derivative for multivariate functions of the form $g: \mathbb{R}^d \to \mathbb{R}$. 
We use a multi-index notation $\alpha \in \mathbb{R}^d$ which defines mixed partial derivatives. 
We denote the $\alpha^{\mathrm{th}}$ derivative of $g$ with respect to its input $\v{x}$ as $D^\alpha g(\v{x})$.
\[
D^{\alpha}g = \frac{\partial^{|\alpha|}g}{\partial x^{\alpha_1}_1\dots\partial x^{\alpha_d}_d}
\]
where $|\alpha| = \sum_{i=1}^d|\alpha_i|$. Note that $\v{x}^\alpha = [x_1^{\alpha_1}, \dots, x_d^{\alpha_d}]$ and $\alpha! = \alpha_1!\dots\cdot\dots\alpha_d!$. 
\vspace{1em}

\begin{definition}[\textnormal{\cite{Cucker2002}}]
Sobolev spaces are denoted $W^{l,p}(\Omega), \Omega \subset \mathbb{R}^d$,  where $l$, the order of the space, is a non-negative integer and $p\geq1$.  The Sobolev space of index $(l,p)$ is the space of locally integrable functions $f: \Omega \to \mathbb{R}$ such that for every multi-index $\alpha$ where $|\alpha| < l$ the derivative $D^{\alpha}f$ exists and $D^{\alpha}f \in L^p(\Omega)$.
The norm in such a space is given by 
$\|f\|_{W^{l,p}(\Omega)} = \left(\sum_{|\alpha| \leq l} \int_{\Omega} |D^{\alpha}f(\v{x})|^p d\v{x}\right)^{\frac{1}{p}}$.
\label{def:sobolev}
\end{definition}
For $p=2$ these spaces are Hilbert spaces, with a dot product that defines the $L_2$ norm of a function's derivatives. Further these Sobolev spaces can be defined in a measure space with \textit{finite} measure $\mu$. 
We call such spaces finite measure spaces of the form $W^{l,p}_{\mu}(\mathbb{R}^d)$ and these are the spaces of locally integrable functions such that for every $\alpha$ where $|\alpha| < l, \ D^{\alpha}f \in L^p_{\mu}(\mathbb{R}^d)$, the $L^p$ space equipped with the measure $\mu$. 
The norm in such a space is given by \citep{Hornik1991}: 
\begin{equation}
\|f\|_{W^{l,p}_{\mu}(\mathbb{R}^d)} = \left(\sum_{|\alpha| \leq l} \int_{\mathbb{R}^d} |D^{\alpha}f(\v{x})|^p d\mu(\v{x})\right)^{\frac{1}{p}}, f \in W^{l,p}_{\mu}(\mathbb{R}^d), |\mu(\v{x})| < \infty \ \forall \v{x} \in \mathbb{R}^d
\label{eq:weighted_sobolev_norm}
\end{equation}

Generally a Sobolev space over a compact subset $\Omega$ of $\mathbb{R}^d$ can be expressed as a weighted Sobolev space with a measure $\mu$ which has compact support on $\Omega$ \citep{Hornik1991}. 

\cite{Hornik1991} have shown that neural networks with continuous activations, which have continuous and bounded derivatives up to order $l$, such as the sigmoid function, are universal approximators in the \textit{weighted} Sobolev spaces of order $l$, meaning that they form a dense subset of Sobolev spaces.
Further, \cite{Czarnecki2017} have shown that networks that use piecewise linear activation functions (such as $
\mathrm{ReLU}$ and its extensions) are \textit{also} universal approximators in the Sobolev spaces of order 1 where the domain $\Omega$ is some compact subset of $\mathbb{R}^d$. 
As mentioned above, this is equivalent to being dense in a weighted Sobolev space on $\mathbb{R}^d$  where the measure $\mu$ has compact support. 
Hence, we can view a neural network, with sigmoid or piecewise linear activations to be a parameter that indexes a function in a weighted Sobolev space with index $(1,2)$, i.e. $f_{\v{\theta}} \in W^{1,2}_{\mu}(\mathbb{R}^d)$.







\section{The Explicit Effect of Gaussian Noise Injections}




Here we consider the case where we noise all layers with \textit{isotropic} noise, except the final predictive layer which we also consider to have no activation function. 
We can express the effect of the Gaussian noise injection on the cost function as an added term $\Delta\mathcal{L}$, which is dependent on $\v{\mathcal{E}}_{L}$, the noise accumulated on the final layer $L$ from the noise additions $\v{\epsilon}$ on the previous hidden layer activations. 
\begin{equation}
\widetilde{\mathcal{L}}(\mathcal{B};\v{\theta}, \v{\epsilon}) =  \mathcal{L}(\mathcal{B}; \v{\theta}) + \Delta\mathcal{L}(\mathcal{B};\v{\theta},\v{\mathcal{E}}_{L})
\end{equation}


To understand the regularisation induced by GNIs, we want to study the regularisation that these injections induce \textit{consistently} from batch to batch. 
To do so, we want to remove the stochastic component of the GNI regularisation and extract a regulariser that is of consistent sign. 
Regularisers that change sign from batch-to-batch do not give a constant objective to optimise, making them unfit as regularisers \citep{Botev2017, Sagun2018, Wei2020}.

As such, we study the explicit regularisation these injections induce by way of the expected regulariser, $\expect_{\v{\epsilon} \sim p(\v{\epsilon})} \left[ \Delta\mathcal{L}(\cdot) \right]$ that marginalises out the injected noise $\v{\epsilon}$. 
To lighten notation, we denote this as  $\expect_{\v{\epsilon}} \left[\Delta\mathcal{L}(\cdot)\right]$.
We extract $R$, a constituent term of the expected regulariser that dominates the remainder terms in norm, and is \textit{consistently positive}.

Because of these properties, $R$ provides a lens through which to study the effect of GNIs. 
As we show, this term has a connection to the Sobolev norm and the Fourier transform of the function parameterised by the neural network. Using these connections we make inroads into better understanding the regularising effect of noise injections on neural networks.

To begin deriving this term, we first need to define the accumulated noise $\v{\mathcal{E}}_{L}$. 
We do so by applying a Taylor expansion to each noised layer.
As in Section \ref{sec:sobolev} we use the generalisation of the derivative for multivariate functions using a multi-index $\alpha$.
For example $D^\alpha h_{k,i}(\v{h}_{k-1}(\v{x}))$ denotes the $\alpha^{\mathrm{th}}$ derivative of the $i^{\mathrm{th}}$ activation of the $k^{\mathrm{th}}$ layer ($h_{k,i}$) with respect to the preceding layer's activations  $\v{h}_{k-1}(\v{x})$ and $D^\alpha \mathcal{L}(\v{h}_{k}(\v{x}), \v{y})$ denotes the $\alpha^{\mathrm{th}}$ derivative of the loss with respect to the non-noised activations $\v{h}_{k}(\v{x})$.



\begin{prop}
Consider an $L$ layer neural network experiencing isotropic GNIs at each layer $k \in [0, \dots, L - 1]$ of dimensionality $d_k$. 
We denote this added noise as $\v{\epsilon} = \{{\v{\epsilon}}_0, ... , {\v{\epsilon}}_{L-1} \}$. 
We assume $\v{h}_L$ is in $C^\infty$ the class of infinitely differentiable functions. 
We can define the accumulated noise at layer each layer $k$ using a multi-index $\alpha_k \in \mathbb{N}^{d_{k-1}}$:
\begin{align*}
    \mathcal{E}_{L,i} &=
    \sum_{|\alpha_L|=1}^{\infty}
    \frac{1}{\alpha_L!}
     \left(D^{\alpha_L} h_{L,i}(\v{h}_{L-1}(\v{x}))\right)
     \v{\mathcal{E}}_{L-1}^{\alpha_L}, \ i = 1, \dots, d_L \\
     \mathcal{E}_{k,i} &= \epsilon_{k,i} +  \sum_{|\alpha_k|=1}^{\infty}
    \frac{1}{\alpha_k!}
    \left(D^{\alpha_k} h_{k,i}(\v{h}_{k-1}(\v{x}))\right) 
    \v{\mathcal{E}}_{k-1}^{\alpha_k} , \ i = 1, \dots, d_k,\ k = 1 \dots L-1  \\  \v{\mathcal{E}}_0 &= \v{\epsilon}_0 \ 
\end{align*}
\label{prop:accum_noise}
where $\v{x}$ is drawn from the dataset $\mathcal{D}$, $\v{h}_{k}$ are the activations before any noise is added, as defined in Equation \eqref{eq:nonoise_acts}. 
\end{prop}

See Appendix \ref{app:accum_noise} for the proof. 
Given this form for the accumulated noise, we can now define the expected regulariser. 
For compactness of notation, we denote each layer's Jacobian as $\v{J}_k \in \mathbb{R}^{d_L \times d_k}$ and the Hessian of the loss with respect to the final layer as $\v{H}_L \in \mathbb{R}^{d_L \times d_L}$. 
Each entry of $\v{J}_{k}$ is a partial derivative  of $f^{\theta}_{k,i}$, the function from layer $k$ to the $i^{\mathrm{th}}$ network output, $i = 1 ... d_L$. 
\begin{align*}
    \v{J}_{k}(\v{x}) = \begin{bmatrix} 
    \frac{f^{\theta}_{k,1}}{\partial h_{k,1}} & \frac{f^{\theta}_{k,1}}{\partial h_{k,2}} & \dots \\
    \vdots & \ddots & \\
    \frac{f^{\theta}_{k,d_L}}{\partial h_{k,1}} &        & \frac{f^{\theta}_{k,d_L}}{\partial h_{k,d_k}} 
    \end{bmatrix} , \ 
    \v{H}_{L}(\v{x}, \v{y}) = \begin{bmatrix} 
    \frac{\partial^2 \mathcal{L}}{\partial h^2_{L,1}} & \frac{\partial^2 \mathcal{L}}{\partial h_{L,1}\partial h_{L,2}} & \dots \\
    \vdots & \ddots & \\
    \frac{\partial^2 \mathcal{L}}{\partial h_{L,d_L}\partial h_{L,1}} &        & \frac{\partial^2 \mathcal{L}}{\partial h^2_{L,d_L}} 
    \end{bmatrix} 
\end{align*}

Using these notations we can now define the explicit regularisation induced by GNIs. 

\begin{theorem}
Consider an $L$ layer neural network experiencing isotropic GNIs at each layer $k \in [0, \dots, L - 1]$ of dimensionality $d_k$. 
We denote this added noise as $\v{\epsilon} = \{{\v{\epsilon}}_0, ... , {\v{\epsilon}}_{L-1} \}$. 
We assume  $\mathcal{L}$ is in $C^\infty$ the class of infinitely differentiable functions. 
We can marginalise out the injected noise $\v{\epsilon}$ to obtain an added regulariser: 
\begin{align*}
\expect_{\v{\epsilon}} \left[ \Delta\mathcal{L}(\mathcal{B};\v{\theta},\v{\mathcal{E}}_{L}) \right] 
    &= \mathbb{E}_{(\v{x},\v{y}) \sim \mathcal{B}} \left[\frac{1}{2}\sum_{k=0}^{L-1}\left[\sigma_k^2\mathrm{Tr}\left(\v{J}^\intercal_{k}(\v{x})
    \v{H}_{L}(\v{x}, \v{y})\v{J}_{k}(\v{x})\right)\right] \right]  + \expect_{\v{\epsilon} } \left[\mathcal{C}(\mathcal{B}, \v{\epsilon} )\right]
    \nonumber 
\end{align*}
where $\v{h}_{k}$ are the activations before any noise is added, as in equation \eqref{eq:nonoise_acts}. $\expect_{\v{\epsilon} } \left[\mathcal{C}(\cdot)\right]$ is a remainder term in higher order derivatives.
\label{prop:explicit_reg}
\end{theorem}
See Appendix \ref{app:exp_reg} for the proof and for the exact form of the remainder. 
We denote the first term in Theorem~\ref{prop:explicit_reg} as:
\begin{equation}
     R(\mathcal{B}; \v{\theta}) = \mathbb{E}_{(\v{x},\v{y}) \sim \mathcal{B}} \left[\frac{1}{2}\sum_{k=0}^{L-1}\left[\sigma_k^2\mathrm{Tr}\left(\v{J}^\intercal_{k}(\v{x})
    \v{H}_{L}(\v{x}, \v{y})\v{J}_{k}(\v{x})\right)\right] \right]
\end{equation}

To understand the main contributors behind the regularising effect of GNIs, we first want to establish the relative importance of the two terms that constitute the explicit effect. 
We know that $R$ is the added regulariser for the linearised version of a neural network, defined by its Jacobian.
This linearisation well approximates neural network behaviour for \textit{sufficiently wide} networks \citep{ntk, chizat:hal-01945578, NEURIPS2019_dbc4d84b} in early stages of training \citep{chen2020generalized}, and we can expect $R$ to dominate the remainder term in norm, which consists of higher order derivatives. 
In Figure \ref{fig:kappadom} we show that this is the case for a range of GNI variances, datasets, and activation functions for networks with 256 neurons per layer; where the remainder is estimated as: 
\[\expect_{\v{\epsilon} } \left[\mathcal{C}(\mathcal{B}, \v{\epsilon})\right] \approx \frac{1}{1000}\sum_{i=0}^{1000}\widetilde{\mathcal{L}}(\mathcal{B};\v{\theta}, \v{\epsilon}) - R(\mathcal{B}; \v{\theta}) - \mathcal{L}(\mathcal{B};\v{\theta}).\] 
These results show that $R$ is a significant component of the regularising effect of GNIs. 
It dominates the remainder $\expect_{\v{\epsilon} } \left[\mathcal{C}(\cdot)\right]$ in norm and is always positive, as we will show, thus offering a consistent objective for SGD to minimise.
Given that $R$ is a likely candidate for understanding the effect of GNIs; we further study this term separately in regression and classification settings.

\begin{figure}[t!]
    \centering
     \subfigure[][BHP Sigmoid]{\includegraphics[width=0.25\textwidth]{figures/BHPSigmoid6layersRvsOkappa.pdf}}
     \subfigure[][CIFAR10 ELU]{\includegraphics[width=0.24\textwidth]{figures/CIFAR10ELU6layersRvsOkappa.pdf}}
     \subfigure[][BHP Sigmoid]{\includegraphics[width=0.25\textwidth]{figures/BHPSigmoid6layersratioovertime.pdf}}
    \subfigure[][CIFAR10 ELU]{\includegraphics[width=0.24\textwidth]{figures/CIFAR10ELU6layersratioovertime.pdf}}
   
    \caption{
    In (a,b) we plot $R(\cdot)$ vs $\expect_{\v{\epsilon} } \left[\mathcal{C}(\cdot)\right]$ at initialisation for 6-layer-MLPs with GNIs at each 256-neuron layer with the same variance $\sigma^2 \in [0.1,0.25,1.0,4.0]$ at each layer. Each point corresponds to one of 250 different network initialisation acting on a batch of size 32 for the classification dataset CIFAR10 and regression dataset Boston House Prices (BHP) datasets. The dotted red line corresponds to $y=x$ and demonstrates that for all batches and GNI variances $R$ is greater than $\expect_{\v{\epsilon} } \left[\mathcal{C}(\cdot)\right]$. In (c,d) we plot $\mathrm{ratio}=|\expect_{\v{\epsilon} } \left[\mathcal{C}(\cdot)\right]|/R(\cdot)$ in the first 100 training iteration for 10 randomly initialised networks. Shading corresponds to the standard deviation of values over the 10 networks. $R(\cdot)$ remains dominant in early stages of training as the ratio is less than 1 for all steps. 
    }
    \label{fig:kappadom}
\end{figure}








\paragraph{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 data label pair $(\v{x}, \v{y})$ as:
\begin{equation}
\mathcal{L}(\v{x}, \v{y}) = \frac{1}{2}(\v{y} - \v{h}_{L}(\v{x}))^2.
\label{eq:mse}
\end{equation}
For this loss, the Hessians in Theorem \ref{prop:explicit_reg} are simply the identity matrix. 
The explicit regularisation term, guaranteed to be positive is: 
\begin{equation}
    R(\mathcal{B}; \v{\theta}) = \frac{1}{2} \mathbb{E}_{\v{x} \sim \mathcal{B}} \left[ \sum_{k=0}^{L-1} \sigma^2_k (\|\v{J}_{k}(\v{x})\|^2_F) \right].
    \label{eq:mse_reg_term}
\end{equation}

where $\sigma^2_k$ is the variance of the noise $\v{\epsilon}_{k}$ injected at layer $k$ and $\|\cdot\|_F$ is the Frobenius norm. 
See Appendix \ref{app:exp_reg_regression} for a proof. 



\paragraph{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 $\v{h}_L$ to be the pre-$\mathrm{softmax}$ of the logits of the final layer.
We denote $\v{p}(\v{x})=\mathrm{softmax}(\v{h}_L(\v{x}))$. 
For a pair $(\v{x}, \v{y})$ we have:
\begin{equation}
\mathcal{L}(\v{x}, \v{y}) = - \sum_{c=0}^C \v{y}_{c} \log (\v{p}(\v{x}))_c),
\label{eq:cross_entropy}
\end{equation}
where $c$ indexes over $C$ possible classes.
The hessian $\v{H}_{L}(\cdot)$ no longer depends on $\v{y}$: 
\begin{equation}
\v{H}_{L}(\v{x})_{i,j} = 
\begin{cases}
\v{p}(\v{x})_i(1 - \v{p}(\v{x})_j) & i = j \\
-\v{p}(\v{x})_i\v{p}(\v{x})_j & i \neq j \\
\end{cases}
\label{eq:hessian_classification}
\end{equation} 
This Hessian is positive-semi-definite  and $R(\cdot)$, guaranteed to be positive, can be written as:
\begin{equation}
    R(\mathcal{B}; \v{\theta}) = \frac{1}{2} \mathbb{E}_{\v{x} \sim \mathcal{B}} \left[ \sum_{k=0}^{L-1} \sigma_k^2 \sum_{i,j}(\mathrm{diag}( \v{H}_{L}(\v{x}))^\intercal{\v{J}}^2_{k}(\v{x}) )_{i,j} \right],
\label{eq:ce_reg_term}
\end{equation}
where $\sigma^2_k$ is the variance of the noise $\v{\epsilon}_{k}$ injected at layer $k$.
See Appendix \ref{app:exp_reg_classification} for the proof. 



\begin{figure}[t!]
    \centering
    \subfigure[][SVHN ELU]{\includegraphics[width=0.54\textwidth]{figures/loss_approximations_add_elu_new_flat_svhn.pdf}}
    \hspace{1cm}
    \subfigure[][SVHN MLP $\mathrm{Tr}(\v{H})$ ]{\includegraphics[width=0.3\textwidth]{figures/hessian_approximations_svhn_add_new_runs.pdf}} \\
    \caption{
    Figure (a) shows the test set loss for convolutional models (CONV) and 4 layer MLPs trained on SVHN with $R(\cdot)$ and GNIs for $\sigma^2=0.1$, and no noise (Baseline).
    Figure (b) shows the trace of the network parameter Hessian for a 2-layer, 32-unit-per-layer MLP where $\v{H}_{i,j} = \frac{\partial \mathcal{L}}{\partial w_i\partial w_j}$, which is a proxy for the parameters' location in the loss landscape.  
    All networks use ELU activations. 
    See Appendix \ref{app:exp_reg_results} for more such results on other datasets and network architectures.
    }
    \label{fig:marg_approx}
\end{figure}

To test our derived regularisers, in Figure~\ref{fig:marg_approx} we show that models trained with $R$ and GNIs have similar training profiles, whereby they have similar test-set loss and parameter Hessians throughout training, meaning that they have almost identical trajectories through the loss landscape. This implies that $R$ is a good descriptor of the effect of GNIs and that we can use this term to understand the mechanism underpinning the regularising effect of GNIs. 
As we now show, it penalises neural networks that parameterize functions with higher frequencies in the Fourier domain; offering a novel lens under which to study GNIs. 

\section{Fourier Domain Regularisation}
\label{sec:Fourier}

To link our derived regularisers to the Fourier domain, we use the connection between neural networks and Sobolev Spaces mentioned above. 
Recall that by \cite{Hornik1991}, we can only assume a sigmoid or piecewise linear neural network parameterises a function in a weighted Sobolev space with measure $\mu$, if we assume that the measure $\mu$ has compact support on a subset $\Omega \in \mathbb{R}^d$.
As such, we equip our space with the \textit{probability} measure $\mu(\v{x})$, which we assume has compact support on some subset $\Omega \subset \mathbb{R}^d$ where $\mu(\Omega)=1$. 
We define it such that $d\mu(\v{x}) = p(\v{x})d\v{x}$ where $d\v{x}$ is the Lebesgue measure and $p(\v{x})$ is the data density function. 
Given this measure, we can connect the derivative of functions that are in the Hilbert-Sobolev space $W^{1,2}_{\mu}(\mathbb{R}^d)$ to the Fourier domain. 

\begin{theorem}
Consider a function, $f: \mathbb{R}^d \to \mathbb{R}$, with a $d$-dimensional input and a single output with $f \in W^{1,2}_{\mu}(\mathbb{R}^d)$ where $\mu$ is a probability measure which we assume has compact support on some subset $\Omega \subset \mathbb{R}^d$ such that $\mu(\Omega)=1$.
Let us assume the derivative of $f$, $D^{\alpha} f$, is in  $L^2(\mathbb{R}^d)$ for some multi-index $\alpha$ where $|\alpha|=1$.  
Then we can write that: 
\begin{align}
 \sum_{|\alpha|=1}\|D^\alpha f\|^2_{L^2_{\mu}(\mathbb{R}^d)} &= 
   \int_{\mathbb{R}^d} \sum^d_{j=1} \Bigr|\mathcal{G}(\v{\omega},j)\overline{ \mathcal{G}(\v{\omega},j) * \mathcal{P}(\v{\omega})} \Bigr|  d\v{\omega}  \\ 
 \mathcal{G}(\v{\omega},j) &=  \v{\omega}_j\mathcal{F}(\v{\omega})  \nonumber
\end{align}
where $\mathcal{F}$ is the Fourier transform of $f$, $\mathcal{P}$ is the Fourier transform or the `characteristic function' of the probability measure $\mu$, $j$ indexes over $\v{\omega} = [\omega_1, \dots, \omega_d]$, $*$ is the convolution operator, and $\overline{(\cdot)}$ is the complex conjugate. 
\label{th:plancherel-measure-space}
\end{theorem}

See Appendix \ref{app:plancherel_th} for the proof. 
Note that in the case where the dataset contains finitely many points, the integrals for the norms of the form $ \|D^\alpha f_{\theta}\|^2_{L^2_{\mu}(\mathbb{R}^d)}$ are approximated by sampling a batch from the dataset which is distributed according to the presumed probability measure $\mu(\v{x})$. 
Expectations over a batch thus approximate integration over $\mathbb{R}^d$ with the measure $\mu(\v{x})$ and this approximation improves as the batch size grows. 
We can now use Theorem \ref{th:plancherel-measure-space} to link $R$ to the Fourier domain.

\begin{figure}[t!]
    \centering
     \hspace{0.4cm}
    \includegraphics[width=0.3\textwidth]{figures/spectral_dynamics_baseline.pdf}
    \includegraphics[width=0.3\textwidth]{figures/spectral_dynamics_noise.pdf}
    \includegraphics[width=0.3\textwidth]{figures/spectral_dynamics_exp_reg.pdf} \\
    \vspace{-0.1cm}
    \hspace{-0.4cm}
    \subfigure[Baseline]{\includegraphics[width=0.3\textwidth]{figures/learnt_function_baseline.pdf}}
    \subfigure[Noise]{\includegraphics[width=0.3\textwidth]{figures/learnt_function_noise.pdf}}
    \subfigure[$R$]{\includegraphics[width=0.3\textwidth]{figures/learnt_function_exp_reg_1.pdf}}
    
    \caption{As in \citet{Rahaman2019}, we train 6-layer deep 256-unit wide ReLU networks trained to regress the function $\lambda(z) = \sum_i \sin(2\pi r_iz+\phi(i))$ with $r_i  \in(5,10,\dots,45,50)$. 
    We train these networks with no noise (Baseline), with GNIs of variance $0.1$ injected into each layer except the final layer (Noise), and with the $R(\cdot)$ for regression in \eqref{eq:mse_reg_term}.
    The first row shows the Fourier spectrum (x-axis) of the networks 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. 
    }
    \label{fig:network_fourier}
\end{figure}

\paragraph{Regression} 
Let us begin with the case of regression.
Assuming differentiable and continuous activation functions, then the Jacobians within $R$ are equivalent to the derivatives in Definition \ref{def:sobolev}.
Theorem \ref{th:plancherel-measure-space} only holds for functions that have 1-D outputs, but we can decompose the Jacobians $\v{J}_{k}$ as the derivatives of multiple 1-D output functions.
Recall, that the $i^\mathrm{th}$ row of the matrix $\v{J}_{k}$ is the set of partial derivatives of $f^{\theta}_{k,i}$, the function from layer $k$ to the $i^{\mathrm{th}}$ network output, $i = 1 ... d_L$, with respect to the $k^\mathrm{th}$ layer activations.
Using this perspective, and the fact that each $f^{\theta}_{k,i} \in W^{1,2}_{\mu}(\mathbb{R}^{d_k})$ ($d_k$ is the dimensionality of the $k^{\mathrm{th}}$ layer), if we assume that the probability measure of our space $\mu(\v{x})$ has compact support, we use Theorem \ref{th:plancherel-measure-space} to write:
\begin{align}
 R(\mathcal{B}; \v{\theta}) &= \frac{1}{2}\mathbb{E}_{\v{x} \sim \mathcal{B}} \left[ \sum_{k=0}^{L-1} \sigma_k^2 \sum_{i}\|\v{J}_{k,i}(\v{x})\|^2_{F} \right] \nonumber \\ &= \frac{1}{2}  \sum_{k=0}^{L-1} \sigma_k^2 \sum_{i} \mathbb{E}_{\v{x} \sim \mathcal{B}} \left[\|\v{J}_{k,i}(\v{x})\|^2_{F} \right] \nonumber \\ 
 &\approx \frac{1}{2}  \sum_{k=0}^{L-1} \sigma_k^2 \sum_{i}  \sum_{|\alpha|=1}\|D^\alpha f^{\theta}_{k,i}\|^2_{L^2_{\mu}(\mathbb{R}^{d_k})} \nonumber \\
 &= \frac{1}{2}\sum_{k=0}^{L-1} \sigma_k^2 \sum_{i} \int_{\mathbb{R}^{d_k}} \sum^{d_k}_{j=1} \Bigr|\mathcal{G}^{\theta}_{k,i}(\v{\omega},j)\overline{ \mathcal{G}^{\theta}_{k,i}(\v{\omega},j) * \mathcal{P}(\v{\omega})}\Bigr|   d\v{\omega}
 \label{eq:sobolev_reg_mse}
\end{align}
where $\v{h}_0 =\v{x}$, $i$ indexes over output neurons, and $\mathcal{G}^{\theta}_{k,i}(\v{\omega}, j) = \v{\omega}_j \mathcal{F}^{\theta}_{k,i}(\v{\omega})$, where $\mathcal{F}^{\theta}_{k,i}$ is the Fourier transform of the function $f^{\theta}_{k,i}$.
The approximation comes from the fact that in SGD, as mentioned above, integration over the dataset is approximated by sampling mini-batches $\mathcal{B}$.



If we take the data density function to be the empirical data density, meaning that it is supported on the $N$ points of the dataset $\mathcal{D}$ (i.e it is a set of $\delta$-functions centered on each point), then as the size $B$ of a batch $\mathcal{B}$ tends to $N$ we can write that:
\begin{align}
 \lim_{B\to N} R(\mathcal{B}; \v{\theta}) &=  \frac{1}{2}\sum_{k=0}^{L-1} \sigma_k^2 \sum_{i} \int_{\mathbb{R}^{d_k}} \sum^{d_k}_{j=1} \Bigr|\mathcal{G}^{\theta}_{k,i}(\v{\omega},j)\overline{ \mathcal{G}^{\theta}_{k,i}(\v{\omega},j) * \mathcal{P}(\v{\omega})}\Bigr|   d\v{\omega}.
\end{align}




\paragraph{Classification}
The classification setting requires a bit more work. 
Recall that our Jacobians are weighted by $\mathrm{diag}(\v{H}_{L}(\v{x}))^\intercal$, which has positive entries that are less than 1 by Equation \eqref{eq:hessian_classification}. We can define a new set of measures such that  $d\mu_i(\v{x}) = \mathrm{diag}(\v{H}_{L}(\v{x}))^\intercal_i p(\v{x}) d\v{x}, \ i = 1 \dots d_L$.
Because this new measure is positive, finite and still has compact support, Theorem \ref{th:plancherel-measure-space} still holds for the spaces indexed by $i$: $W^{1,2}_{\mu_i}(\mathbb{R}^d)$. 

Using these new measures, and the fact that each $f^{\theta}_{k,i} \in W^{1,2}_{\mu_i}(\mathbb{R}^{d_k})$, we can use Theorem \ref{th:plancherel-measure-space} to write that for classification models: 
\begin{align}
 R(\mathcal{B}; \v{\theta}) &=  \frac{1}{2}  \sum_{k=0}^{L-1} \sigma_k^2 \sum_{i} \mathbb{E}_{\v{x} \sim \mathcal{B}} \left[\mathrm{diag}(\v{H}_{L}(\v{x}))^\intercal_i \|\v{J}_{k,i}(\v{x})\|^2_{2} \right] 
 \nonumber\\ &\approx \frac{1}{2}  \sum_{k=0}^{L-1} \sigma_k^2 \sum_{i} \sum_{|\alpha|=1}\|D^\alpha f^{\theta}_{k,i}\|^2_{L^2_{\mu_i}(\mathbb{R}^{d_k})} \nonumber \\
 &=  \frac{1}{2}\sum_{k=0}^{L-1} \sigma_k^2 \sum_{i} \int_{\mathbb{R}^{d_k}}  \sum^{d_k}_{j=1} \Bigr|\mathcal{G}^{\theta}_{k,i}(\v{\omega}, j)\overline{ \mathcal{G}^{\theta}_{k,i}(\v{\omega}, j) * \mathcal{P}_i(\v{\omega})}\Bigr|   d\v{\omega}
 \label{eq:sobolev_reg_ce}
\end{align}



Here $\mathcal{P}_i$ is the Fourier transform of the $i^\mathrm{th}$ measure $\mu_i$ and as before $\mathcal{G}^{\theta}_{k,i}(\v{\omega}, j) = \v{\omega}_j \mathcal{F}^{\theta}_{k,i}(\v{\omega})$, where $\mathcal{F}^{\theta}_{k,i}$ is the Fourier transform of the function $f^{\theta}_{k,i}$.
Again as the batch size increases to the size of the dataset, this approximation becomes exact.

For both regression and classification, GNIs, by way of $R$, induce a prior which favours smooth functions with low-frequency components.
This prior is enforced by the terms $\mathcal{G}^{\theta}_{k,i}(\v{\omega}, j)$ which become large in magnitude when functions have high-frequency components, penalising neural networks that learn such functions. In Appendix \ref{app:tikhonov} we also show that this regularisation in the Fourier domain corresponds to a form of \textit{Tikhonov regularisation}.  

In Figure \ref{fig:network_fourier}, we demonstrate empirically that networks trained with GNIs learn functions that don't overfit; with lower-frequency components relative to their non-noised counterparts. 


\paragraph{A layer-wise regularisation} Note that there is a recursive structure to the penalisation induced by $R$. 
Consider the layer-to-layer functions which map from a layer $k-1$ to $k$, $\v{h}_{k}(\v{h}_{k-1}(\v{x}))$. 
With a slight abuse of notation, $\nabla_{\v{h}_{k-1}}  \v{h}_{k}(\v{x})$ is the Jacobian defined element-wise as: 
\[\left(\nabla_{\v{h}_{k-1}}  \v{h}_{k}(\v{x})\right)_{i,j} = \frac{\partial h_{k,i}}{\partial h_{k-1,j}(\v{x})},\]
where as before $h_{k,i}$ is the $i^\mathrm{th}$ activation of layer $k$. 

$\|\nabla_{\v{h}_{k-1}}  \v{h}_{k}(\v{x})\|_2^2$ is penalised $k$ times in $R$ as this derivative appears in $\v{J}_0, \v{J}_1 \dots \v{J}_{k-1}$ due to the chain rule.
As such, when training with GNIs, we can expect the norm of $\|\nabla_{\v{h}_{k-1}}  \v{h}_{k}(\v{x})\|_2^2$ to decrease as the layer index $k$ increases (i.e the closer we are to the network output).
By Theorem \ref{th:plancherel-measure-space}, and Equations \eqref{eq:sobolev_reg_mse}, and \eqref{eq:sobolev_reg_ce}, larger  $\|\nabla_{\v{h}_{k-1}}  \v{h}_{k}(\v{x})\|_2^2$ correspond to functions with higher frequency components.
Consequently, when training with GNIs the layer to layer function $\v{h}_{k}(\v{h}_{k-1}(\v{x}))$ will have higher frequency components than the next layer's  function $\v{h}_{k+1}(\v{h}_{k}(\v{x}))$.

\begin{figure}[t!]

    \centering
    \subfigure[][Baseline  ]{\includegraphics[width=0.245\textwidth]{figures/layer_layer_grad_baseline.pdf}} 
    \subfigure[][GNI  ]{\includegraphics[width=0.245\textwidth]{figures/layer_layer_grad_noise.pdf}}
    \subfigure[][Baseline  ]{\includegraphics[width=0.245\textwidth]{figures/tr_w_baseline.pdf}} 
    \subfigure[][GNI  ]{\includegraphics[width=0.245\textwidth]{figures/tr_w_noise.pdf}} 

    \caption{We use 6-layer deep 256-unit wide ReLU networks on the same dataset as in Figure \ref{fig:network_fourier} trained with (GNI) and without GNI (Baseline). 
    In (a,b), for layers with square weight matrices, we plot 
    $\|\widetilde{\v{W}}_k\|_2^2$.
    In (c,d) we plot the trace of these layers' weight matrices $\mathrm{Tr}(\v{W}_k)$. 
    For GNI models, as the layer index $k$ increases, $\mathrm{Tr}(\v{W}_k)$ and  $\|\widetilde{\v{W}}_k\|_2^2$ decrease, indicating that each successive layer in these networks learns a function with lower frequency components than the past layer.  
    }
    \label{fig:per_layer_reg}
\end{figure}

We measure this layer-wise regularisation in $\mathrm{ReLU}$ networks, using $\nabla_{\v{h}_{k-1}}  \v{h}_{k}(\v{x}) = \widetilde{\v{W}}_k$.
$\widetilde{\v{W}}_k$ is obtained from the original weight matrix $\v{W}_k$ by setting its $i^{\mathrm{th}}$ column to zero whenever the neuron $i$ of the $k^{\mathrm{th}}$ layer is inactive.
Also note that the inputs of $\mathrm{ReLU}$ network hidden layers, which are the outputs of another $\mathrm{ReLU}$-layer, will be positive.
Negative weights are likely to `deactivate' a $\mathrm{ReLU}$-neuron, inducing smaller $\|\widetilde{\v{W}}_k\|_2^2$, and thus parameterising a lower frequency function.
We use the trace of a weight matrix as an indicator for the `number' of negative components.

In Figure \ref{fig:per_layer_reg} we demonstrate that $\|\widetilde{\v{W}}_k\|_2^2$ \textit{and} $\mathrm{Tr}(\v{W}_k)$  decrease as $k$ increases for $\mathrm{ReLU}$-networks trained with GNIs, indicating that each successive layer in these networks learns a function with lower frequency components than the past layer. 
This striation and ordering is clearly absent in the baselines trained without GNIs. 









\paragraph{The Benefits of GNIs}

What does regularisation in the Fourier domain accomplish? 
The terms in $R$ are the traces of the Gauss-Newton decompositions of the second order derivatives of the loss. 
By penalising this we are more likely to land in wider (smoother) minima (see Figure \ref{fig:marg_approx}), which has been shown,  contentiously \citep{Dinh2017}, to induce networks with better generalisation properties \citep{Keskar2019OnMinima, Jastrzebski2017}. 
GNIs however, confer other benefits too.

\textit{Sensitivity.}  A model's weakness to input perturbations is termed the \textit{sensitivity}. 
\citet{Rahaman2019} have shown empirically that classifiers biased towards lower frequencies in the Fourier domain are less sensitive,
and there is ample evidence demonstrating that models trained with noised data are less sensitive \citep{Liu2018, Li2018}.
The Fourier domain - sensitivity connection can be established by studying the \textit{classification margins} of a model (see Appendix \ref{app:classification_margins}). 


\textit{Calibration.}
Given a network's prediction $\hat{y}(\v{x})$ with confidence $\hat{p}(\v{x})$ for a point $\v{x}$, perfect calibration consists of being as likely to be correct as you are confident: $p(\hat{y}=y|\hat{p}=r)=r, \,\, \forall r\in[0,1]$ \citep{Dawid1982, Degroot1982TheForecasters}. 
In Appendix \ref{app:capacity} we show that models that are biased toward lower frequency spectra have lower `capacity measures', which measure model complexity and lower values of which have been shown empirically to induce better calibrated models \citep{Guo2017}. 
In Figure \ref{fig:calibration_app}
we show that this holds true for models trained with GNIs.


\section{Related Work}

Many variants of GNIs have been proposed to regularise neural networks.
\citet{Poole2014} extend this process 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, \citet{Bishop1995} derive an analytic form for the explicit regulariser induced by noise injections on \textit{data} and show that such injections are equivalent to Tikhonov regularisation in an unspecified function space. 

Recently \citet{Wei2020} conducted similar analysis to ours, dividing the effects of Bernoulli dropout into \textit{explicit} and \textit{implicit} effects. 
Their work is built on that of \citet{Mele1993}, \citet{Helmbold2015}, and \citet{wager2013dropout} who perform this analysis for linear neural networks. 
\citet{Arora2020} derive an explicit regulariser for Bernoulli dropout on the final layer of a neural network.
Further, recent work by \citet{dieng2018noisin} shows that noise additions on recurrent network hidden states outperform Bernoulli dropout in terms of performance and bias. 


\section{Conclusion}

In this work, we derived analytic forms for the explicit regularisation induced by Gaussian noise injections, demonstrating
that the explicit regulariser penalises networks with high-frequency content in Fourier space.
Further we show that this regularisation is not distributed evenly within a network, as it disproportionately penalises high-frequency content in layers closer to the network output.
Finally we demonstrate that this regularisation in the Fourier domain has a number of beneficial effects. 
It induces training dynamics that preferentially land in wider minima, it reduces model sensitivity to noise, and induces better calibration.






\newpage

\section*{Acknowledgments}

This research was directly funded by the Alan Turing Institute under Engineering and Physical Sciences Research Council (EPSRC) grant EP/N510129/1.
AC was supported by an EPSRC Studentship.
MW was supported by EPSRC grant EP/G03706X/1.
U\c{S} was supported by the French National Research Agency (ANR) as a part of the FBIMATRIX (ANR-16-CE23-0014) project.
SR gratefully acknowledges support from the UK Royal Academy of Engineering and the Oxford-Man Institute.
CH was supported by the Medical Research Council, the Engineering and Physical Sciences Research Council, Health Data Research UK, and the Li Ka Shing Foundation

\section*{Impact Statement}

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 distribution shifts. 
As such improving our understanding of such methods is critical as machine learning models become increasingly ubiquitous and embedded in decision making. 


\bibliographystyle{plainnat}
