commit b6ef1c9c9a8da2729f18917d1ba916723b0d7994
parent 4bc2e66b1d4db5e99a63e5bcf53db33d3878a9a1
Author: Silas Brack <s174433@student.dtu.dk>
Date: Thu, 2 Feb 2023 11:08:13 +0100
Cleaning up.
Diffstat:
4 files changed, 32 insertions(+), 21 deletions(-)
diff --git a/chapters/experiments.tex b/chapters/experiments.tex
@@ -2,9 +2,11 @@
\chapter{Experiments}
\labch{experiments}
-WARNING: THIS SECTION IS BY NO MEANS FINISHED
+WARNING: THIS CHAPTER IS BY NO MEANS FINISHED
-\section[Hessian-Vector Products]{Hessian-Vector Products}[Hessian-Vector Products]
+\section{Hessian-Vector Products}[Hessian-Vector Products]
+
+THINKING OF MOVING THIS SECTION TO \cref{ch:intro} OR \cref{ch:laplace}
In this project, we propose performing the Laplace approximation by only storing the computational graph for the implicit Hessian-vector product, thereby avoiding the explicit computation and storage of the quadratically-scaling Hessian matrix.
How significant of an effect does this have, though?
@@ -72,6 +74,8 @@ This is in line with the theory, as noted in \cref{sec:practical-ggn} (SHOULD IT
\section{Sampling Ablation}
+THINKING OF MOVING THIS SECTION TO \cref{ch:sampling}
+
\begin{marginfigure}
\centering
\includegraphics{ablation/mnist_ablation.pdf}
@@ -83,8 +87,6 @@ Contour integral quadrature (CIQ) converges to the exact solution for \(\lim_{Q
However, CIQ is typically close to the exact solution for low values of \(Q\).\sidenote{\textcite{pleiss2020fast} find that \(Q=20\) is approximately sufficient for most problems.}
We then investigated the number of quadrature points \(Q\) required to determine an adequate approximate solution for the Laplace approximation.
-
-
We vary \(Q \in \{1, \ldots, 20\}\) to determine the number of quadrature points required for the Laplace approximation.
We then compute the Mahalanobis distance between the samples and the true posterior mean, and compare this to the chi-squared distribution with \(D\) degrees of freedom.
We then compute the p-value of the Kolmogorov-Smirnov test for each value of \(Q\) by comparing the empirical CDF of the Mahalanobis distances to the theoretical \(\chi^2_D\) CDF.
@@ -120,7 +122,7 @@ In particular, for models with a very large number of parameters, the time requi
Even if the sampling time scales sublinearly with the number of parameters, it would still prevent methods such as online Laplace, which require the Laplace approximation to be computed at each iteration, from being used in practice.
However, this sampling time constraint is also significant for MCMC methods, which are the most common method for sampling from the posterior distribution of a BNN.
-\section[Sampling for the Sine Function]{Sampling for the Sine Function}[The Sine Function]
+\section{Sampling for the Sine Function}[The Sine Function]
The choice of \(\alpha\) is a non-trivial one
@@ -134,15 +136,24 @@ The choice of \(\alpha\) is a non-trivial one
\begin{figure}
\centering
\includegraphics{chisq_sin.pdf}
- \caption[Chi-squared sample plot for sine curve.]{Chi-squared sample plot for visualising normality of 1000 posterior samples from the full Laplace approximation for a fully connected neural network modelling \(f(x) = \sin(x)\). Since the sample square Mahalanobis distances appear chi-squared distributed, the samples themselves are multivariate normal with correct (fixed) covariance. As can be seen by looking at the histogram and the quantile--quantile plot, samples are multivariate normal.}
+ \caption[Chi-squared sample plot for sine curve.]{Chi-squared sample plot for visualising normality of 1000 posterior samples from the full Laplace approximation for a fully connected neural network modelling \(f(x) = \sin(x)\). Since the sample square Mahalanobis distances appear chi-squared distributed, the samples themselves are multivariate normal with the correct (known) covariance.}
\label{fig:chisq-sin}
\end{figure}
\section[Sampling for MNIST]{Sampling for MNIST}[MNIST]
+Sampling was successful for values of alpha as low as 0.1,
+
\begin{marginfigure}
\centering
\includegraphics{eigenvalues_mnist_55k_0.1_old.pdf}
\caption[Spectrum of the posterior precision on MNIST.]{Hyperexponential distribution of eigenvalues of Laplace posterior precision on MNIST with a convolutional NN with 17k parameters and a 3 by 3 filter. There are two distinct clusters in the spectrum, with the first cluster being over an order of magnitude larger than the second. This first cluster consists of nine eigenvalues, and so we can conclude that the posterior covariance has an effective rank close to nine. FIX THIS PLOT.}
\label{fig:spectrum-mnist}
\end{marginfigure}
+
+\section{Discussion}
+
+The choice of prior precision \(\alpha\) is non-trivial.
+Furthermore, since the prior precision is the lowest eigenvalue of the posterior precision matrix, we
+While we can ``choose'' \(\alpha\) to be a value large enough to reduce conditioning issues, this somewhat invalidates the Bayesian interpretation of the prior, which should theoretically be chosen based on prior knowledge of the problem, can also be learned from the data.
+Choosing the prior in order to reduce ill-conditioning of our posterior precision means we sacrifice some of the benefits of the prior.
diff --git a/chapters/laplace.tex b/chapters/laplace.tex
@@ -1,12 +1,12 @@
\setchapterpreamble[u]{\margintoc}
\chapter{The Laplace Approximation}
-\labch{laplace}
+\label{ch:laplace}
% http://saul.cpsc.ucalgary.ca/pmwiki.php/GradTips/GradTipsLiteratureReview
\section[Optimisation in Deep Learning]{Optimisation in Deep Learning}[Optimisation]
\label{sec:optimisation-frequentist}
-Suppose a neural network is a real-valued function \(f: \reals^N \times \reals^D \rightarrow \reals^O\) parameterised in \(\bm\theta\) which maps an input \(\bm x\) to an output \(f(\bm x; \bm\theta) \equiv f_{\bm\theta}(\bm x)\).
+Suppose a neural network is a real-valued function \(f: \reals^N \times \reals^D \rightarrow \reals^O\) parameterised in \(\bm\theta\) which maps an input \(\bm x\) to an output \(f(\bm x, \bm\theta) \equiv f_{\bm\theta}(\bm x)\).
Our goal is to find the parameters \(\bm\theta^*\) which best model the observed data.
From a frequentist perspective, we can define a loss function on the model \(\mathcal{L}(\bm\theta): \reals^D \rightarrow \reals\) such that the optimal parameters \(\bm\theta^*\) minimise this loss function for the given model.
@@ -19,7 +19,7 @@ These two cases correspond to the mean-squared error (MSE) and cross-entropy (CE
Due to the assumption of independence between observations, we can factorise the likelihood and the negative log-likelihood loss as
\begin{align}\label{eq:negative-log-likelihood}
- p(\bm y \given \bm \theta) ={} & \prod_i^N p(x_i \given \bm \theta),
+ p(\bm y \given \bm \theta) ={} & \prod_i^N p(y_i \given \bm \theta),
\\ \mathcal{L}(\bm\theta) \equiv \mathcal{L}(\bm \theta, \bm y) \coloneqq -\log p(\bm y \given \bm \theta) ={}& - \sum_i^N \log p(y_i \given \bm \theta).
\end{align}
This allows us to compute the loss as the sum of the losses for each individual observation.
@@ -38,7 +38,7 @@ The update rule is given by
where \(\eta_t\) is the learning rate at iteration \(t\).
The learning rate \(\eta_t\) can be constant or adaptive, and is often chosen to be a non-increasing function of \(t\).
The gradient \(\nabla_{\bm\theta} \mathcal{L}(\bm\theta_t, \bm x_t, y_t)\) can be computed separately for each single training example \(\bm x_t, y_t\), as seen in \cref{eq:negative-log-likelihood}.
-Using the chain rule on the composition of the loss into the loss as a function of the model output \(\mathcal{L}\) and the model output as a function of the parameters \(f\) such that \(\mathcal{L}(\bm \theta_t, \bm x_t, y_t) = \mathcal{L}(y_t, f(\bm x_t; \bm \theta_t))\),we compute the gradient as
+Using the chain rule on the composition of the loss into the loss as a function of the model output \(\mathcal{L}\) and the model output as a function of the parameters \(f\) such that \(\mathcal{L}(\bm \theta_t, \bm x_t, y_t) = \mathcal{L}(y_t, f(\bm x_t, \bm \theta_t))\),we compute the gradient as
%
\begin{align}\label{eq:sgd-gradient}
\nabla_{\bm\theta} \mathcal{L}(\bm\theta_t, \bm x_t, y_t) ={} & \nabla_{\bm\theta} \log p(y_t \given \bm \theta) \nabla_{\bm\theta} f_{\bm\theta_t}(\bm x_t),
@@ -102,7 +102,7 @@ The update rule for Adam is thus given by
\\ \bm m_t ={}& (1 - \beta_1^t) \left[ \beta_1 \bm m_{t-1} + (1 - \beta_1) \nabla_{\bm\theta} \mathcal{L}(\bm\theta_t) \right],
\\ \bm v_t ={}& (1 - \beta_2^t) \left[ \beta_2 \bm v_{t-1} + (1 - \beta_2) \left( \nabla_{\bm\theta} \mathcal{L}(\bm\theta_t) \right)^2 \right],
\end{align}
-where \(\odot\) denotes element-wise multiplication, \(\varepsilon\) is a small constant to prevent division by zero, and \(\bm m_t\) and \(\bm v_t\) are the first and second moments of the gradient respectively.
+where \(\odot\) denotes the Hadamard product, \(\varepsilon\) is a small constant to prevent division by zero, and \(\bm m_t\) and \(\bm v_t\) are the first and second moments of the gradient respectively.
The momentum term~\cite{polyak1964some} is inspired by physics, where the momentum of an object is the product of its mass and velocity.
In Adam, the first moment of the gradient is used to compute the momentum of the parameters, increasing the convergence speed of the optimisation~\cite{sutskever2013importance}.
@@ -137,7 +137,7 @@ Furthermore, since the unnormalised posterior from \cref{eq:unnormalised posteri
Similarly to the negative log-likelihood, we can factorise the negative log-posterior loss as per
\begin{align}\label{eq:negative-log-posterior}
\mathcal{L}(\bm\theta) = -\log p(\bm \theta \given \bm y) ={} & - \log p(\bm y \given \bm \theta) - \log p(\bm\theta) \nonumber
- \\ ={}& \underbrace{-\sum_i^N \log p(x_i \given \bm \theta)}_{\text{Negative log-likelihood}} - \underbrace{\log p(\bm\theta)}_{\text{Regularisation}}.
+ \\ ={}& \underbrace{-\sum_i^N \log p(y_i \given \bm \theta)}_{\text{Negative log-likelihood}} - \underbrace{\log p(\bm\theta)}_{\text{Regularisation}}.
\end{align}
Here the regularisation term is given by the prior distribution \(p(\bm\theta)\).
Due to the black-box nature of neural networks, it is often difficult to define a prior distribution over the parameters using previous knowledge.
diff --git a/chapters/sampling.tex b/chapters/sampling.tex
@@ -1,6 +1,6 @@
\setchapterpreamble[u]{\margintoc}
\chapter{Sampling}
-\labch{sampling}
+\label{ch:sampling}
\section{Introduction}
@@ -16,7 +16,7 @@ Because of this, the covariance matrix is given by the inverse of this Hessian,
\begin{align}\label{eq:sample-laplace}
\bm \varepsilon = \bm \mu + \bm \Lambda^{-1/2} \bm \varepsilon_0,
\end{align}
-where \(\bm \Lambda = \nabla^2_{\bm\theta} \log p(\bm \theta \given \bm x)\).
+where \(\bm \Lambda = \nabla^2_{\bm\theta} \log p(\bm \theta \given \bm y)\).
In \cref{eq:sample-reparam}, the square root of a matrix yields another matrix which, when multiplied by itself, is equal to the original matrix.
However, since a Gaussian distribution is invariant to rotation, in this special case, we just need to obtain a matrix which returns the original matrix \emph{up to a rotation}.\sidenote{That is, if we have a vector of standard normal samples \(\bm \varepsilon \sim \normal(\bm 0, \identity)\), then, for an orthogonal matrix \(\bm R\) (i.e., a rotation matrix), \(\bm R\T \bm \varepsilon \sim \normal(\bm 0, \identity)\). We then have that we can sample \emph{non}-standard normal samples as \(\bm \Lambda^{-1/2} \bm R\T \bm \varepsilon \sim \normal(\bm 0, \identity)\).}
We thus need to find a method which computes the inverse square root of a high-dimensional matrix up to a rotation.
@@ -305,7 +305,7 @@ If we calculate the GGN approximation of the Hessian by linearizing over the who
\\ \bm P^{-1} ={}& \beta \cdot \nabla \nabla\T + \alpha^{-1} \identity
\\ \beta ={}& \left(\alpha \norm{\nabla}^2 + \norm{\nabla}^4\right)^{-1} - \left(\alpha \norm{\nabla}^2\right)^{-1},
\end{align}
-where \(\nabla\) is the gradient of the loss with regard to the model parameters.\sidenote{DO I NEED TO PROVE THIS?}
+where \(\nabla\) is the gradient of the loss with regard to the model parameters.\sidenote{DO I NEED TO PROVE THIS? BLAH BLAH BLAH EIGENDECOMPOSITION}
The problem with using this approximation is that the rank of \(\bm P\) is 1.
This means that, best-case scenario, we can reduce the highest eigenvalue of the Hessian.
This occurs when the gradient corresponds to the direction of the highest precision (empirically, this may not be the case).
@@ -348,7 +348,7 @@ This identity can be used to compute the inverse of the GGN matrix with the prio
\begin{align}
\left( \alpha \identity + \bm J\T \bm H \bm J \right)^{-1} \bm v = \alpha^{-1} \identity - \alpha^{-2} \bm J\T \left( \bm H^{-1} + \alpha^{-1} \bm J \bm J\T \right)^{-1} \bm J,
\end{align}
-where \(\alpha\) is the prior precision.
+where \(\alpha\) is the prior precision and \(\bm H\) and \(\bm J\) are the Hessian and the Jacobian as per \cref{eq:laplace-precision-ggn} evaluated at a single observation.
Here, to exploit the low-rank structure of the preconditioner, we only compute the GGN matrix over a single observation.
As such, for this preconditioner to be effect, we assume that the eigenvectors of the GGN matrix are constant across the dataset.
@@ -441,7 +441,7 @@ Early tests of the randomly pivoted Cholesky preconditioner were not successful,
Additionally, we considered potential preconditioners to be computed as traditional approximations of the posterior precision, like the diagonal approximation or the Kronecker factorisation.
These methods have been shown to be relatively effective approximations of the posterior precision in the context of Bayesian neural networks~\cite{daxberger2021laplace}, and so we considered them as potential preconditioners, as they may contain enough information about the full posterior precision to be effective preconditioners.
-However, due to the speed, effectiveness, and theoretical properties of other preconditioners proposed in this section, these methods were not attempted.
+However, due to the speed, effectiveness, and theoretical properties of other preconditioners proposed in this section, these methods were not used.
Overall, both the sub-sampled preconditioner and the pivoted Cholesky preconditioners were found to be effective under the correct settings.
However, the pivoted Cholesky preconditioner is not easy to implement in JAX and was more difficult to tune than the sub-sampled preconditioner.
The sub-sampled preconditioner is thus the preferred preconditioner in this project.
@@ -472,8 +472,8 @@ Samples \(\bm \varepsilon\) which are normally distributed but do not have a mea
\end{align}
%
This has the advantage of only requiring computation of the precision matrix \(\bm \Lambda\).
-For the Laplace approximation, this can easily be computed as $\bm \Lambda = - \nabla^2_{\bm \theta} \log p(\bm \theta \given \bm x)$.
-Since this precision matrix can be computed exactly, this calculation can be performed to evaluate the whether a set of samples \(\bm \varepsilon\) are drawn from a distribution \(\normal(\bm \mu, \bm \Sigma)\) while only accessing the precision matrix in the evaluation.
+For the Laplace approximation, this can easily be computed as $\bm \Lambda = - \nabla^2_{\bm \theta} \log p(\bm \theta \given \bm y)$.
+Since this precision matrix can be computed exactly, this calculation can be performed to evaluate the whether a set of samples \(\bm \varepsilon\) are drawn from a distribution \(\normal(\bm \mu, \bm \Lambda^{-1})\) while only accessing this precision matrix in the evaluation.
To do so, we can compute the squared Mahalanobis distance of each sample \(\bm \varepsilon_i\) from the mean \(\bm \mu\) and compare it to the chi-squared distribution with \(D\) degrees of freedom.
We do this by comparing the histogram of these squared Mahalanobis distances to the appropriate PDF and by visualising the sample empirical percentiles against the theoretical percentiles in a quantile-quantile plot.
As such, this method can be used to evaluate whether the approximate normal samples obtained from CIQ are correctly distributed.
diff --git a/chapters/training.tex b/chapters/training.tex
@@ -1,13 +1,13 @@
\setchapterpreamble[u]{\margintoc}
\chapter{Training}
-\labch{training}
+\label{ch:training}
\section{Maximising the Evidence}
\label{sec:maximising-the-evidence}
As seen in \cref{eq:bayes-theorem}, the posterior probability for a given model is equal to
\begin{align}
- p(\bm \theta \given \bm y) = \frac{p(\bm y \given \bm \theta) p(\bm \theta)}{p(\bm y)},
+ p(\bm \theta \given \bm y) = \frac{p(\bm y \given \bm \theta)\,p(\bm \theta)}{p(\bm y)},
\end{align}
where \(p(\bm y)\) is known as the evidence. It's also known as the \emph{marginal likelihood}, since it's obtained by marginalising out \(\bm \theta\) from the likelihood, as per
\begin{align}
