commit cb6de7a253628739abb0dc2ddefa6eb5af819b60
parent 0b1eb6b8410b08e8d413f1d5da92cac3b7fd9c4d
Author: Silas Brack <s174433@student.dtu.dk>
Date: Mon, 23 Jan 2023 13:28:37 +0100
More writing.
Diffstat:
5 files changed, 126 insertions(+), 92 deletions(-)
diff --git a/chapters/introduction.tex b/chapters/introduction.tex
@@ -7,8 +7,12 @@
Bayesian methods allow for estimates of uncertainty which enable more
efficient usage of data (i.e., via active learning) and avoid
-overfitting. Furthermore, these uncertainty estimates improve model
+overfitting.
+Furthermore, these uncertainty estimates improve model
interpretability and assessment of model predictive confidence.
+These methods have been successfully applied to a wide range of classical problems in statistics, though attempts to apply them to deep learning have had limited success.
+are particularly useful in the context of deep learning,
+where the models are often black-boxes and the data is often scarce.
% \section{The State of Research in the Area}
\section{Current Methods}
@@ -35,7 +39,7 @@ As such, there is demand for a method which exhibits the same computational cost
The Laplace approximation~\sidecite{daxberger2021laplace,laplace1774memoire,} is a simple, yet principled, posterior approximation suitable for Bayesian modeling.
In the Laplace approximation, to compute the posterior distribution of the model parameters, we need to compute the Hessian of the loss with respect to the model parameters.
-Typically, this is done by storing the Hessian matrix explicitly, and then computing the inverse of this matrix.
+Typically, this is done by storing the Hessian matrix explicitly and computing the inverse of this matrix.
Since the Hessian matrix is a \(D \times D\) square matrix with \(D^2\) elements, where \(D\) is the number of parameters, this is intractable for large models, which can contain millions or billions of parameters.
To overcome this, we can use a crude approximation, such as only storing the Hessian matrix's diagonal~\sidecite{lecun1989optimal,denker1990transforming}, its Kronecker factorisation~\sidecite{heskes2000natural,martens2015optimizing,botev2017practical}, and other low-rank approximations.%~\sidecite{}.
All of these Hessian approximations have successfully been applied to the Laplace approximation~\sidecite{ritter2018scalable,ritter2018online}
diff --git a/chapters/literature.tex b/chapters/literature.tex
@@ -44,12 +44,12 @@ Note that Adam approximates the second moment of the gradient of the loss, that
\section{Laplace Approximation}\label{sec:laplace}
In the Laplace approximation (LA)~\sidecite{laplace1774memoire,bishop1995neural,mackay2003information,daxberger2021laplace}, the posterior is approximated by a Gaussian, similarly to mean-field VI.
-However, instead of finding the optimal Gaussian distribution locations and scales by maximising the ELBO, the LA finds the location by computing the MAP solution \(\bm \theta_{\textsc{map}}\) and the scale by approximating the log posterior with a second degree Taylor expansion around this solution (\(\bm \theta_0 = \bm \theta_{\textsc{map}}\)) and determining the curvature via its Hessian matrix \(\bm \Lambda = - \left.\nabla^2_{\bm \theta} \log p(\bm \theta \given \bm y) \right|_{\bm \theta_\mathrm{MAP}}\).
+However, instead of finding the optimal Gaussian distribution locations and scales by maximising the ELBO, the LA finds the location by computing the MAP solution \(\bm \theta_{\textsc{map}}\) and the scale by approximating the log posterior with a second degree Taylor expansion around this solution (\(\bm \theta_0 = \bm \theta_{\textsc{map}}\)) and determining the curvature via its Hessian matrix \(\bm \Lambda = - \left.\nabla^2_{\bm \theta} \log p(\bm \theta \given \bm y) \right|_{\bm \theta_\textsc{map}}\).
Since the Taylor expansion is performed around the MAP solution, the first order derivative is zero, and the expansion is simply given by
\begin{align}
\ln p(\bm \theta \given \bm y)
\approx{} & \ln p(\bm \theta_{\textsc{map}} \given \bm y) + \nonumber
- \\ & + \frac{1}{2} (\bm \theta - \bm \theta_{\textsc{map}})\T \left( \left.\nabla^2_{\bm \theta} \ln p(\bm \theta \given \bm y) \right|_{\bm \theta_\mathrm{MAP}} \right) (\bm \theta - \bm \theta_{\textsc{map}}) \label{eq:laplace-taylor}
+ \\ & + \frac{1}{2} (\bm \theta - \bm \theta_{\textsc{map}})\T \left( \left.\nabla^2_{\bm \theta} \ln p(\bm \theta \given \bm y) \right|_{\bm \theta_\textsc{map}} \right) (\bm \theta - \bm \theta_{\textsc{map}}) \label{eq:laplace-taylor}
\\ \Rightarrow \tilde{p}(\bm \theta \given \bm y) \approx{}& p(\bm \theta_{\textsc{map}} \given \bm y) \exp\left(- \frac{1}{2} (\bm \theta - \bm \theta_{\textsc{map}})\T \bm \Lambda (\bm \theta - \bm \theta_{\textsc{map}})\right) \propto p(\bm \theta \given \bm y) \nonumber
\end{align}
where \(\tilde{p}(\bm \theta \given \bm y)\) corresponds to the unnormalised posterior.
@@ -57,7 +57,7 @@ Normalising it yields
\begin{align}\label{eq:laplace-normal}
p(\bm\theta \given \bm y) \approx{} & \sqrt{\frac{\det (\bm \Lambda)}{(2 \pi)^D}} \exp\left(- \frac{1}{2} (\bm \theta - \bm \theta_{\textsc{map}})\T \bm \Lambda (\bm \theta - \bm \theta_{\textsc{map}})\right) \nonumber
% \\={}& \sqrt{\frac{\det (- \bm H)}{(2 \pi)^D}} \exp\left(- \frac{1}{2} (\bm \theta - \bm \theta_0)\T (- \bm H) (\bm \theta - \bm \theta_0)\right) \nonumber
- \\={}& \normal(\bm \theta_{\textsc{map}}, \bm \Lambda^{-1}).
+ \\={}& \normal(\bm \theta \given \bm \theta_{\textsc{map}}, \bm \Lambda^{-1}).
\end{align}
\marginnote{\cref{eq:laplace-normal} only holds when the precision matrix \(\bm \Lambda\) is positive definite.}We have now shown that approximating the log-posterior with a second degree Taylor expansion around the \(\bm\theta_{\textsc{map}}\) corresponds to approximating the posterior with a Gaussian distribution given by \(\normal(\bm \theta_{\textsc{map}}, \bm \Lambda^{-1})\) where \(\bm \Lambda = - \left.\nabla^2_{\bm \theta} \log p(\bm \theta \given \bm y) \right|_{\bm \theta_{\textsc{map}}}\).
@@ -136,9 +136,9 @@ minimum eigenvalue determined by the perturbation \(\alpha\) which guarantees po
% \subsection{Maximum a Posteriori Estimation}
% \emph{Maximum a posteriori} estimation finds a point estimate of the posterior given by its maximum.
-% It can therefore be interpreted as a Delta distribution estimate of the posterior \(p(\bm \theta \given \bm y) \approx \delta(\bm \theta_\mathrm{MAP})\), where \(\bm \theta_{\text{MAP}}\) is defined as
+% It can therefore be interpreted as a Delta distribution estimate of the posterior \(p(\bm \theta \given \bm y) \approx \delta(\bm \theta_\textsc{map})\), where \(\bm \theta_{\text{MAP}}\) is defined as
% \begin{align}
-% \bm \theta_\mathrm{MAP}
+% \bm \theta_\textsc{map}
% ={} & \max_{\bm \theta} p(\bm \theta \given \bm y) \nonumber \\
% ={} & \max_{\bm \theta} \log p(\bm \theta \given \bm y) \nonumber \\
% ={} & \max_{\bm \theta} \left[ \log \lik + \log p(\bm \theta) \right] \nonumber \\
diff --git a/chapters/sampling.tex b/chapters/sampling.tex
@@ -21,6 +21,13 @@ In \cref{eq:sample-reparam}, the square root of a matrix yields another matrix w
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.
+Some methods exist which compute the inverse square root of a matrix up to a rotation, such as the Cholesky decomposition.
+However, these methods are not suitable for large matrices, as they require the computation of the full matrix, which is not feasible for large matrices.
+There are also methods which can compute the inverse square root of a low-rank matrix that is constructed as the inner product of two rectangular matrices by performing the inverse square root on the outer product (which is low-dimensional for low-rank problems), such as the Woodbury matrix identity.
+However, the precision matrix we are interested in is not itself low-rank but a sum of low-rank matrices, and thus this method is not suitable for our problem.
+This is discussed in more detail in \cref{sec:laplace,sec:preconditioning}.
+As such, these are not adequate methods for computing the inverse square root using only matrix-vector products.
+
\section[Contour Integral Quadrature]{Contour Integral Quadrature}[Contour Integral Quad]
\sidetextcite{pleiss2020fast} propose one such method.
@@ -63,7 +70,7 @@ it is then also trivial to compute the square root product by multiplying \(\bm
\def\xmax{3.0}
\def\xmin{-0.5}
\def\R{1.3}
-\begin{marginfigure}[-7cm]
+\begin{marginfigure}%[-7cm]
\centering
\input{figures/cauchy_integral.tex}
\caption[Representation of integration circle for the Cauchy integral formula.]{Representation of integration circle \(\Gamma\) contained within the domain of \(f\) as defined according to the Cauchy integral formula.}
@@ -178,34 +185,35 @@ Additionally, since msMINRES is an iterative Krylov method, it is also sensitive
To mitigate the effect of ill conditioning, we can use a preconditioner \(\bm P\) such that \(\bm P^{-1} \bm K\) is well-conditioned.
\section{Preconditioning}
+\label{sec:preconditioning}
One drawback of contour integral quadrature if that the accuracy its of approximation of \(\bm \Lambda^{-1/2} \bm \varepsilon_0\) suffers when the precision matrix \(\bm \Lambda\) is poorly conditioned.\sidenote{A matrix is considered poorly conditioned when the ratio of its largest eigenvalue to its lowest eigenvalue, its \emph{conditioning number} (\(\kappa := \lambda_{\mathrm{max}} / \lambda_{\mathrm{min}}\)), is very large.}
In the case of the Laplace approximation, as per \cref{eq:ggn-laplace}, we have that the precision matrix is given by
\begin{align}\label{eq:ggn-laplace}
\bm \Lambda = \sum_i^N \bm J_i\T \bm H \bm J_i + \alpha \identity,
- \\ \bm \theta \sim \normal(\bm \mu, \bm \Lambda^{-1})
+ \\ \bm \theta \sim \normal(\bm \theta_{\textsc{map}}, \bm \Lambda^{-1})
\end{align}
where \(N\) is the number of observations in your dataset.
Since each \(J_i\T H J_i\) is, for \(i \in 1 \ldots N\), a low-rank positive semi-definite matrix, it has \(o\) non-zero eigenvalues.\sidenote{With rank \(o\), where \(o\) is the number of neural network outputs.}
Let us define \(\lambda_1, \ldots, \lambda_o\) as the eigenvalues of a matrix sorted in decreasing order.
We can then
-Assume that the non-zero spectrum of \(J_i\T H J_i\) is relatively flat, i.e., \(\lambda_1 \not\gg \lambda_o\).
+Assume that the non-zero spectrum of \(\bm J_i\T \bm H \bm J_i\) is relatively flat, i.e., \(\lambda_1 \not\gg \lambda_o\).
In a hand-wavy way, we have that the ``more similar'' the eigenvectors (and assuming that the eigenvectors have the same ordering between the two matrices) are, the more the largest eigenvalue will be to the sum of the eigenvalues.
For preconditioning, we need to find a preconditioner matrix \(P\) and its inverse \(P^{-1}\) such that
\begin{align}
- P^{-1} A \approx{} & \identity \Rightarrow P \approx A
- \\ P^{-1} \Sigma^{-1} ={}& P^{-1} \left( \textstyle{\sum^N_i} J_i\T H J_i + \alpha \identity \right)
+ \bm P^{-1} \bm A \approx{} & \identity \Rightarrow \bm P \approx \bm A
+ \\ \bm P^{-1} \bm \Lambda ={}& \bm P^{-1} \left( \textstyle{\sum^N_i} \bm J_i\T \bm H \bm J_i + \alpha \identity \right)
\end{align}
-where \(\Sigma^{-1}\) is our precision matrix, approximated by a GGN Hessian.
+where \(\bm \Lambda\) is our precision matrix, approximated by a GGN Hessian.
The eigenvalues of our GGN approximation are given by
\begin{align}
- \lambda\left( \textstyle{\sum^N_i} J_i\T H J_i + \alpha \identity \right) = \big\{ \underbrace{s_1, \ldots, s_p}_{\leq N \cdot k}, \underbrace{\alpha, \ldots, \alpha}_{\geq n - N \cdot k} \big\}
+ \lambda\left( \textstyle{\sum^N_i} \bm J_i\T \bm H \bm J_i + \alpha \identity \right) = \big\{ \underbrace{s_1, \ldots, s_p}_{\leq N \cdot k}, \underbrace{\alpha, \ldots, \alpha}_{\geq n - N \cdot k} \big\}
\end{align}
-where \(\{s_1, \ldots, s_p\}, p \leq N \cdot k\) are the non-zero eigenvalues of the rank-\(p\) matrix \(\sum^N_i J_i\T H J_i\).
+where \(\{s_1, \ldots, s_p\}, p \leq N \cdot k\) are the non-zero eigenvalues of the rank-\(p\) matrix \(\sum^N_i \bm J_i\T \bm H_i \bm J_i\).
If a preconditioner approximates all eigenvalues and eigenvectors of a matrix, then the preconditioner and the matrix are equal and the preconditioner is exact.
However, we often do not have a preconditioner which is exact, but rather an approximation of the matrix.
@@ -226,25 +234,25 @@ The more values you add, the steeper the spectrum becomes (even the non-zero one
Typically, when preconditioning, you need a matrix which is similar to the matrix you are preconditioning.
In this way, the preconditioner acts as an ``initial best guess'' of the value to be computed.
-As such, for CIQ, it would be expected that the preconditioner should approximate the value of \(H^{-1/2}\).
-In reality, we only need an approximation of \(H^{-1}\).
-WAIT IS THIS TRUE? WE STILL NEED \(H^{1/2} b\), WE JUST CALCULATE THIS VALUE USING CIQ.
+As such, for CIQ, it would be expected that the preconditioner should approximate the value of \(\bm \Lambda^{-1/2}\).
+In reality, we only need an approximation of \(\bm \Lambda\).
+WAIT IS THIS TRUE? WE STILL NEED \(\bm \Lambda^{1/2} \bm b\), WE JUST CALCULATE THIS VALUE USING CIQ.
-As such, the conditioning number of our preconditioner also becomes an issue, since we have two conflicting objectives: to have a preconditioner which is as similar to \(H\) as possible and to have a preconditioner which we can invert easily.
-The first point lends itself toward the use of preconditioners with a similar conditioning number as \(H\).
+As such, the conditioning number of our preconditioner also becomes an issue, since we have two conflicting objectives: to have a preconditioner which is as similar to \(\bm \Lambda\) as possible and to have a preconditioner which we can invert easily.
+The first point lends itself toward the use of preconditioners with a similar conditioning number as \(\bm \Lambda\).
However, to invert an arbitrary positive semi-definite matrix, we will often resort to using the conjugate gradient method which is sentitive to the conditioning number of the preconditioner.
-Thus, there is an intrinsic balance to be struck between a preconditioner which is very similar to \(H\) (but whose conditioning number is itself an issue) and a preconditioner which is less similar to \(H\) (and will therefore not improve the conditioning of the problem).
+Thus, there is an intrinsic balance to be struck between a preconditioner which is very similar to \(\bm \Lambda\) (but whose conditioning number is itself an issue) and a preconditioner which is less similar to \(\bm \Lambda\) (and will therefore not improve the conditioning of the problem).
-Let us consider the toy example of a preconditioner which we can guarantee to have the same eigenvectors as \(H\) and for which we can control the eigenvalues.
+Let us consider the toy example of a preconditioner which we can guarantee to have the same eigenvectors as \(\bm \Lambda\) and for which we can control the eigenvalues.
One example of such a preconditioner would be???
-Depending on the rank of the preconditioner, we can then eliminate as many of the largest / lowest eigenvalues of \(H\) as the rank of the preconditioner.
-If we have a rank-one preconditioner whose only eigenvector is the largest eigenvector of \(H\), the second-largest eigenvalue becomes the largest.
+Depending on the rank of the preconditioner, we can then eliminate as many of the largest / lowest eigenvalues of \(\bm \Lambda\) as the rank of the preconditioner.
+If we have a rank-one preconditioner whose only eigenvector is the largest eigenvector of \(\bm \Lambda\), the second-largest eigenvalue becomes the largest.
In the best-case scenario, you can then eliminate
-If \(H\) has a very flat spectrum, then the problem cannot be ill conditioned and thus preconditioning is immaterial.
+If \(\bm \Lambda\) has a very flat spectrum, then the problem cannot be ill conditioned and thus preconditioning is immaterial.
% This analogy can be extended in a hand-wavy way to an imperfect preconditioner which doesn't exactly approximate the eigenvalue
Finally, performance of computing
-In CIQ, we will, in each iteration, perform one matrix-vector product \(H v\) and one preconditioner-vector product \textbf{check this}, as seen in \cref{alg:ciq}.\sidenote{Because of this, to avoid the preconditioner becoming a performance bottleneck, the speed of the preconditioner-vector products should be significantly greater than that of the matrix-vector products.}
+In CIQ, we will, in each iteration, perform one matrix-vector product \(\bm \Lambda \bm v\) and one preconditioner-vector product \textbf{check this}, as seen in \cref{alg:ciq}.\sidenote{Because of this, to avoid the preconditioner becoming a performance bottleneck, the speed of the preconditioner-vector products should be significantly greater than that of the matrix-vector products.}
As such, if the preconditioner-vector product is similarly or more expensive to compute compared to the matrix-vector product, then preconditioning will ....
However, an effective preconditioner will significantly improve the convergence speed of CIQ, and this effect can often outweigh the cost of computing the preconditioner.
@@ -252,12 +260,13 @@ However, an effective preconditioner will significantly improve the convergence
\label{sec:adam-preconditioner}
If we calculate the GGN approximation of the Hessian by linearizing over the whole model and loss (as described in \cref{sec:practical-ggn}), we get
\begin{align}
- P ={} & \nabla \nabla\T + \alpha \identity
- \\ 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}.
+ \bm P ={} & \nabla \nabla\T + \alpha \identity
+ \\ \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.
DO I NEED TO PROVE THIS?
-The problem with using this approximation is that the rank of \(P\) is 1.
+The problem with using this approximation is that the rank of \(\bm P\) is 1.
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, it seems that \emph{this is not the case}).
This would then leave the second-highest eigenvalue unaffected.
@@ -300,29 +309,29 @@ Since \(N \gg B\), then the cost of computing the preconditioner will be negligi
\sidecite{harbrecht2012low,bach2013sharp,gardner2018gpytorch}
\begin{align}
- \Sigma^{-1} ={} & \sum_i^N J_i\T J_i + \alpha \identity
- \\ \approx{}& L\T L + \alpha \identity
+ \bm \Lambda ={} & \sum_i^N \bm J_i\T \bm H_i \bm J_i + \alpha \identity
+ \\ \approx{}& \bm L\T \bm L + \alpha \identity
\end{align}
-where \(L \in \reals^{k \times n}\), with \(k\) being the chosen rank of the pivoted Cholesky factorisation and \(n\) being the number of parameters.
+where \(\bm L \in \reals^{k \times n}\), with \(k\) being the chosen rank of the pivoted Cholesky factorisation and \(n\) being the number of parameters.
We can then use this approximation as a preconditioner by computing its inverse via the Woodbury matrix identity (\cref{sec:woodbury}):
\begin{align}
- P ={} & L\T L + \alpha \identity
- \\ P^{-1} ={}& \alpha^{-1} \identity - \alpha^{-2} L\T \left( \identity + \alpha^{-1} L L\T \right)^{-1} L
+ \bm P ={} & \bm L\T \bm L + \alpha \identity
+ \\ \bm P^{-1} ={}& \alpha^{-1} \identity - \alpha^{-2} \bm L\T \left( \identity + \alpha^{-1} \bm L \bm L\T \right)^{-1} \bm L
\end{align}
we then obtain the preconditioner-vector and inverse-preconditioner-vector products as in Equation~\ref{eq:pivoted-cholesky-preconditioner}.
The ``scariest'' computation in the inverse preconditioner is inverting a \(k \times k \) matrix, where \(k\) is the rank of our pivoted Cholesky factorisation (which we choose).
\begin{align}\label{eq:pivoted-cholesky-preconditioner}
- P v ={} & L\T ( L v ) + \alpha v
- \\ P^{-1} v ={}& \alpha^{-1} v - \alpha^{-2} L\T ( \underbrace{\identity + \alpha^{-1} L L\T}_{k \times k} )^{-1} ( L v )
+ \bm P \bm v ={} & \bm L\T ( \bm L \bm v ) + \alpha \bm v
+ \\ \bm P^{-1} \bm v ={}& \alpha^{-1} \bm v - \alpha^{-2} \bm L\T ( \underbrace{\identity + \alpha^{-1} \bm L \bm L\T}_{k \times k} )^{-1} ( \bm L \bm v )
\end{align}
-We can precompute the factorisation-inverse product \(L\T ( \identity + \alpha^{-1} L L\T )^{-1} \in \reals^{n \times k}\), meaning that we only need to perform two matrix multiplications for both the preconditioner- and inverse preconditioner-vector products.
-This would then require storing the factorisation \(L \in \reals^{k \times n}\) and the factorisation-inverse product \(\reals^{n \times k}\) matrix.\sidenote{We can also precompute the inverse \(( \identity + \alpha^{-1} L L\T )^{-1} \in \reals^{k \times k}\), which means we would only need to store the factorisation \(L \in \reals^{k \times n}\) and this inverse \(k \times k\) matrix, in which case we would just need to perform one more matrix multiplication (multiplying \(L\T\) with this \(k \times k\) matrix) in every preconditioner call.}
+We can precompute the factorisation-inverse product \(\bm L\T ( \identity + \alpha^{-1} \bm L \bm L\T )^{-1} \in \reals^{n \times k}\), meaning that we only need to perform two matrix multiplications for both the preconditioner- and inverse preconditioner-vector products.
+This would then require storing the factorisation \(\bm L \in \reals^{k \times n}\) and the factorisation-inverse product \(\reals^{n \times k}\) matrix.\sidenote{We can also precompute the inverse \(( \identity + \alpha^{-1} \bm L \bm L\T )^{-1} \in \reals^{k \times k}\), which means we would only need to store the factorisation \(\bm L \in \reals^{k \times n}\) and this inverse \(k \times k\) matrix, in which case we would just need to perform one more matrix multiplication (multiplying \(\bm L\T\) with this \(k \times k\) matrix) in every preconditioner call.}
ADD SOMEWHERE:
As we sum over multiple low-rank matrices, we will likely find that the minimum non-zero eigenvalue gets progressively smaller compared to the largest.
When we have a single rank \(o\) GGN matrix, we have
-The preconditioning of \(J J\T\) (which should be low-rank) will progressively become both \emph{less} low-rank and its condition number will progressively become higher. MERGE WITH:
+The preconditioning of \(\bm J_i \bm H_i \bm J_i\T\) (which should be low-rank) will progressively become both \emph{less} low-rank and its condition number will progressively become higher. MERGE WITH:
Regarding the eigenvalues and eigenvectors, we consider a best-case scenario.
Assume that the pivoted Cholesky factorisation finds the \(k\) largest eigenvectors.
@@ -336,7 +345,7 @@ Thus, when the batch size \(B\) is large, the condition number will increase and
How much the condition number increases with the batch size depends on the eigenspectrum of the per-observation Hessians (i.e., how close the eigenvectors are to each other).
REWRITE THIS, IT'S KINDA HARD TO EXPLAIN.
-Unfortunately, if we want to use the pivoted Cholesky decomposition to approximate \(\sum^B_i J_i\T J_i\), we are largely trying to approximate the higher eigenvectors of the sum.
+Unfortunately, if we want to use the pivoted Cholesky decomposition to approximate \(\sum^B_i \bm J_i\T \bm H_i \bm J_i\), we are largely trying to approximate the higher eigenvectors of the sum.
In order for this to be the case, we will want to oversample vectors from this matrix which correspond to the largest eigenvalues.
A common approach to sampling the largest eigenvalues is to sample vectors with probability proportional to the diagonal value~\sidecite{harbrecht2012low}.
However, it is not possible to generally and efficiently determine the diagonal of a matrix via matrix-vector products.
@@ -345,7 +354,7 @@ Often, it is necessary to multiply a vector of zeros with a value of one at the
Alternatively, we can approximate the diagonal of the GGN matrix as the element-wise product of the gradient with itself, similarly to the preconditioner described in \cref{sec:adam-preconditioner} and the GGN approximation computed in \cref{eq:fully-linear-ggn}.
In conclusion, to find the optimal parameters for the pivoted Cholesky preconditioner (the number of observations \(B\) and the rank of the pivoted Cholesky factorisation \(k\)), there are some rules of thumb to consider.
-Notably, the quality of the preconditioner depends largely on the ease of inverting the \(\identity + \alpha^{-1} L L\T\) matrix, which depends on its conditioning number.
+Notably, the quality of the preconditioner depends largely on the ease of inverting the \(\identity + \alpha^{-1} \bm L \bm L\T\) matrix, which depends on its conditioning number.
In practice, increasing the number
I've found that it's most effective to keep \(k\) low, such that the \(k \times k\) outer product of the pivoted Cholesky factorisation matrix \(L\) approximates few enough of the largest eigenvectors that the eigenvalue of the largest is close to that of the smallest (meaning this outer product is well-conditioned).
You thus need a large enough \(k\) to approximate enough large eigenvectors and a small enough \(k\) to ensure the conditioning number of the outer product is low.
@@ -370,7 +379,7 @@ However, due to the speed, effectiveness, and theoretical properties of other pr
\subsection[A Recap on Preconditioning]{A Recap on Preconditioning}[A Recap]
-It is worth noting the issues that arise from scaling the term in the \(\sum J\T H J\) preconditioner to account for the summing over a Hessian of a limited batch size.\sidenote{Or, more precisely, the term which is equivalent to the \(J\T H J\) term in the Hessian and is consequently used to approximate the largest eigenvalue.}
+It is worth noting the issues that arise from scaling the term in the \(\sum \bm J_i\T \bm H_i \bm J_i\) preconditioner to account for the summing over a Hessian of a limited batch size.\sidenote{Or, more precisely, the term which is equivalent to the \(\bm J_i\T \bm H_i \bm J_i\) term in the Hessian and is consequently used to approximate the largest eigenvalue.}
\begin{enumerate}
\item Never scale your preconditioner up
@@ -382,16 +391,16 @@ It is worth noting the issues that arise from scaling the term in the \(\sum J\T
\varepsilon_0\T \varepsilon_0 \sim \chi^2(D)
\end{align}
-Samples \(\varepsilon\) which are normally distributed but do not have a mean of zero and a standard deviation of one can be normalised and their sum of square deviations be calculated by computing their Mahalanobis distances.
-
+Samples \(\varepsilon\) which are normally distributed but do not have a mean of zero and a standard deviation of one can be normalised and their sum of square deviations be calculated by computing their Mahalanobis distances, which will be Chi-square distributed with \(D\) degrees of freedom, as per
+%
\begin{align}
- d ={} & \sqrt{(\varepsilon - \mu)\T \Sigma^{-1} (\varepsilon - \mu)}\label{eq:mahalanobis-distance}
- \\ \Rightarrow{}& d\T d \sim \chi^2(D) \nonumber
+ \bm d ={} & \sqrt{(\bm \varepsilon - \bm \mu)\T \bm \Lambda (\bm \varepsilon - \bm \mu)}\label{eq:mahalanobis-distance}
+ \\ \Rightarrow{}& \bm d\T \bm d \sim \chi^2(D) \nonumber
\end{align}
-
-This has the advantage of only requiring computation of the precision matrix \(\Sigma^{-1}\).
-For the Laplace approximation, this can easily be computed as $\Sigma^{-1} = - \nabla^2_\theta \log p(\theta \given x)$.
-Since \(- \nabla^2_\theta \log p(\theta \given x)\) can be computed exactly, this calculation can be performed to evaluate the whether a set of samples \(\varepsilon\) are drawn from a distribution \(\normal(\mu, \Sigma)\) while only accessing the precision matrix in the evaluation.
+%
+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.
As such, this method can be used to evaluate whether the approximate normal samples obtained from CIQ are correctly distributed.
\section{Limitations}
diff --git a/chapters/theory.tex b/chapters/theory.tex
@@ -4,6 +4,7 @@
% 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)\).
Our goal is to find the optimal parameters \(\bm\theta^*\) which best model the observed data.
@@ -12,20 +13,17 @@ From a frequentist perspective, we can define a loss function \(\mathcal{L}(\bm\
Often, we formulate the loss function as the negative log-likelihood of the data under the model, \(\mathcal{L}(\bm\theta) = -\log \lik\).
In this case, the goal of optimisation is to find the parameters \(\bm\theta^*\) which maximise the likelihood of the data.\sidenote[][5cm]{Maximising the likelihood is equivalent to minimising the negative log likelihood, since the log function is strictly monotonically increasing.}
-From a Bayesian perspective, we define a prior distribution \(p(\bm\theta)\) and a likelihood function \(\lik\) such that the goal of optimisation is to find the parameters \(\bm\theta^*\) which maximise the posterior distribution \(p(\bm\theta \given \bm y)\).
-This way, we can quantify the uncertainty in our model parameters \(\bm\theta\) by computing the posterior distribution \(p(\bm\theta \given \bm y)\).
+The likelihood function is problem-specific, and is defined as the probability of the data given the parameters \(\bm\theta\).
+Typically, the likelihood is chosen to be a Gaussian distribution with mean \(f_{\bm\theta}(\bm x)\) and precision \(\rho\) (which is frequently chosen to be equal to one) for regression problems, and a categorical distribution with logits \(f_{\bm\theta}(\bm x)\) for classification problems.
+These two cases correspond to the mean-squared error (MSE) and cross-entropy (CE) negative log-likelihood loss functions, respectively.
-Due to the assumption of independence between observations, we can factorise the negative log-likelihood loss as
+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 x \given \bm \theta) ={} & \prod_i^N p(x_i \given \bm \theta) \nonumber
- \\ \mathcal{L}(\bm\theta) = -\log p(\bm x \given \bm \theta) ={}& - \sum_i^N \log p(x_i \given \bm \theta).
-\end{align}
-
-Similarly, for the negative log-posterior loss, we have
-\begin{align}\label{eq:negative-log-posterior}
- \mathcal{L}(\bm\theta) = -\log p(\bm \theta \given \bm x) ={} & - \log p(\bm x \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}}
+ \\ \mathcal{L}(\bm\theta) \equiv \mathcal{L}(\bm \theta, \bm x) := -\log p(\bm x \given \bm \theta) ={}& - \sum_i^N \log p(x_i \given \bm \theta).
\end{align}
+This allows us to compute the loss as the sum of the losses for each individual observation.
+This loss function is defined in an optimisation problem to find the parameters \(\bm\theta^*\) of the model which minimise this loss, thereby maximising the performance of the model on the data.
\subsection[Stochastic Gradient Descent]{Stochastic Gradient Descent}[Gradient Descent]
@@ -38,27 +36,16 @@ 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 decreasing function of \(t\).
-The gradient \(\nabla_{\bm\theta} \mathcal{L}(\bm\theta_t)\) is computed using a single training example \(\bm x_t, y_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}.
The gradient is computed using the chain rule as
%
\begin{equation*}
\nabla_{\bm\theta} \mathcal{L}(\bm\theta_t) = \nabla_{\bm\theta} \log \lik(\bm x_t, y_t \given \bm\theta_t) \nabla_{\bm\theta} f_{\bm\theta_t}(\bm x_t)
\end{equation*}
%
-where \(\nabla_{\bm\theta} f_{\bm\theta_t}(\bm x_t)\) is the gradient of the neural network with respect to its parameters \(\bm\theta_t\), and \(\eta_t\) can be interpreted as the step size of the gradient descent algorithm.
+where \(\nabla_{\bm\theta} f_{\bm\theta_t}(\bm x_t)\) is the gradient of the neural network output with respect to its parameters at iteration \(t\), \(\bm\theta_t\), and \(\eta_t\) can be interpreted as the step size of the gradient descent algorithm.
It can either be constant throughout training (\(\eta_t = \eta\)) or adaptive (if it varies throughout training).
-\subsection{Adam}
-
-% How does adam work?
-
-\sidecite{sutskever2013importance}
-\sidecite{kingma2014adam}
-
-% Integrate with the previous section
-Adam approximates the second moment of the gradient of the loss, that this corresponds to the variance of this gradient, and that this is the diagonal of the Fisher information matrix, as we defined in \cref{sec:fisher-information} (\cref{eq:gradient-covariance})!
-% This is a common approximation in practice, as the diagonal of the Fisher information matrix is often sufficient to approximate the curvature of the log-likelihood.
-
\subsection[Natural Gradient Descent]{Natural Gradient Descent}[The Natural Gradient]
% https://agustinus.kristia.de/techblog/2018/03/14/natural-gradient/
@@ -94,6 +81,17 @@ Instead, it is often approximated by the diagonal of the Hessian.
For a recent, detailed discussion of the natural gradient, see \sidetextcite{martens2020new}.
\sidecite{wu2019logan}
+\subsection{Adam}
+
+% How does adam work?
+
+\sidecite{sutskever2013importance}
+\sidecite{kingma2014adam}
+
+% Integrate with the previous section
+Adam approximates the second moment of the gradient of the loss, that this corresponds to the variance of this gradient, and that this is the diagonal of the Fisher information matrix, as we define in \cref{sec:fisher-information} (\cref{eq:gradient-covariance}).
+% This is a common approximation in practice, as the diagonal of the Fisher information matrix is often sufficient to approximate the curvature of the log-likelihood.
+
However, these optimisation methods all identify a single set of parameters as the optimal solution.
In the next section, we will discuss how to learn a distribution over parameters instead.
@@ -101,38 +99,61 @@ In the next section, we will discuss how to learn a distribution over parameters
\label{sec:bayesian-deep-learning}
Bayesian deep learning is a framework for deep learning that allows for modelling the parameters as random variables instead of as single optima in the loss landscape.
-It built from Bayes' theorem, which allows for uncertainty quantification by
+From a Bayesian perspective, we define a prior distribution \(p(\bm\theta)\) and a likelihood function \(\lik\) such that the goal of optimisation is to find the parameters \(\bm\theta^*\) which maximise the posterior distribution \(p(\bm\theta \given \bm y)\).
+This way, we can quantify the uncertainty in our model parameters \(\bm\theta\) by computing the posterior distribution \(p(\bm\theta \given \bm y)\).
+It is obtained from Bayes' theorem, which allows for uncertainty quantification by
\begin{align}\label{eq:bayes-theorem}
- p(\bm \theta \given \bm x) = \frac{p(\bm x \given \bm \theta)\,p(\bm \theta)}{p(\bm x)}
+ p(\bm \theta \given \bm x) = \frac{p(\bm x \given \bm \theta)\,p(\bm \theta)}{p(\bm x)},
+\end{align}
+where \(p(\bm \theta \given \bm x)\) is the posterior distribution over the parameters, \(p(\bm \theta)\) is the prior distribution over parameters, \(p(\bm x \given \bm \theta)\) is the likelihood of the data given the parameters, and \(p(\bm x)\) is the marginal likelihood of the data.
+The marginal likelihood is obtained by marginalising over the parameters, i.e., \(p(\bm x) = \int p(\bm x \given \bm \theta)\,p(\bm \theta)\,d\bm \theta\).
+Due to this integral, the marginal likelihood is often intractable.
+However, since the posterior distribution is proportional to the product of the likelihood by the prior, we can obtain the unnormalised posterior from the likelihood and the prior\marginnote{Since the marginal likelihood is a normalisation constant and doesn't depend on the parameters, it does not affect the optimisation problem.}
+\begin{align}
+ p(\bm \theta \given \bm x) \propto p(\bm x \given \bm \theta)\,p(\bm \theta) =: \tilde{p}(\bm \theta \given \bm x) .
+\end{align}
+
+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 x) ={} & - \log p(\bm x \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}}.
+\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.
+Because of this, the prior is often chosen to be a Gaussian distribution with zero mean and a diagonal covariance matrix with some precision \(\alpha\), i.e.,
+\begin{align}
+ p(\bm\theta) = \normal(\bm\theta \given \bm 0, \alpha^{-1}\bm I).
\end{align}
-Instead of learning a single set of parameters, Bayesian deep learning learns a distribution over parameters.
+The likelihood term is defined as in the frequentist setting (see \cref{sec:optimisation-frequentist}).
+However, the posterior distribution itself is typically intractable, so approximations to the posterior are used instead.
-However, the posterior distribution is typically intractable, so approximations to the posterior are used instead.
+% Instead of learning a single set of parameters, Bayesian deep learning learns a distribution over parameters.
-\section{Laplace Approximation}
+\section[The Laplace Approximation]{The Laplace Approximation}[Laplace Approximation]
\label{sec:laplace}
In the Laplace approximation (LA)~\sidecite{laplace1774memoire,bishop1995neural,mackay2003information,daxberger2021laplace}, the posterior is approximated by a Gaussian, similarly to mean-field VI\@.
-However, instead of finding the optimal Gaussian distribution locations and scales by maximising the ELBO, the LA finds the location by computing the MAP solution \(\bm \theta_{\textsc{map}}\) and the scale by approximating the log posterior with a second degree Taylor expansion around this solution (\(\bm \theta_0 = \bm \theta_{\textsc{map}}\)) and determining the curvature via its Hessian matrix \(\bm \Lambda = - \left.\nabla^2_{\bm \theta} \log p(\bm \theta \given \bm y) \right|_{\bm \theta_\mathrm{MAP}}\).
+However, instead of finding the optimal Gaussian distribution locations and scales by maximising the ELBO, the LA finds the location by computing the MAP solution \(\bm \theta_{\textsc{map}}\) and the scale by approximating the log posterior with a second degree Taylor expansion around this solution (\(\bm \theta_0 = \bm \theta_{\textsc{map}}\)) and determining the curvature via its Hessian matrix \(\bm \Lambda = - \left.\nabla^2_{\bm \theta} \log p(\bm \theta \given \bm y) \right|_{\bm \theta_\textsc{map}}\).
Since the Taylor expansion is performed around the MAP solution, the first order derivative is zero, and the expansion is given by
\begin{align}
\ln p(\bm \theta \given \bm y)
\approx{} & \ln p(\bm \theta_{\textsc{map}} \given \bm y) + \nonumber
- \\ & + \frac{1}{2} (\bm \theta - \bm \theta_{\textsc{map}})\T \left( \left.\nabla^2_{\bm \theta} \ln p(\bm \theta \given \bm y) \right|_{\bm \theta_\mathrm{MAP}} \right) (\bm \theta - \bm \theta_{\textsc{map}}) \label{eq:laplace-taylor}
+ \\ & + \frac{1}{2} (\bm \theta - \bm \theta_{\textsc{map}})\T \left( \left.\nabla^2_{\bm \theta} \ln p(\bm \theta \given \bm y) \right|_{\bm \theta_\textsc{map}} \right) (\bm \theta - \bm \theta_{\textsc{map}}) \label{eq:laplace-taylor}
\\ \Rightarrow \tilde{p}(\bm \theta \given \bm y) \approx{}& p(\bm \theta_{\textsc{map}} \given \bm y) \exp\left(- \frac{1}{2} (\bm \theta - \bm \theta_{\textsc{map}})\T \bm \Lambda (\bm \theta - \bm \theta_{\textsc{map}})\right) \propto p(\bm \theta \given \bm y) \nonumber
\end{align}
where \(\tilde{p}(\bm \theta \given \bm y)\) corresponds to the unnormalised posterior.
-Normalising it yields
+Normalising it by the normalisation constant \(p(\bm x)\) gives the normalised posterior\sidenote{This normalisation constant is computed as the marginal likelihood for the Laplace approximation in \cref{sec:maximising-the-evidence}.}
\begin{align}\label{eq:laplace-normal}
p(\bm\theta \given \bm y) \approx{} & \sqrt{\frac{\det (\bm \Lambda)}{(2 \pi)^D}} \exp\left(- \frac{1}{2} (\bm \theta - \bm \theta_{\textsc{map}})\T \bm \Lambda (\bm \theta - \bm \theta_{\textsc{map}})\right) \nonumber
% \\={}& \sqrt{\frac{\det (- \bm H)}{(2 \pi)^D}} \exp\left(- \frac{1}{2} (\bm \theta - \bm \theta_0)\T (- \bm H) (\bm \theta - \bm \theta_0)\right) \nonumber
- \\={}& \normal(\bm \theta_{\textsc{map}}, \bm \Lambda^{-1}).
+ \\={}& \normal(\bm \theta \given \bm \theta_{\textsc{map}}, \bm \Lambda^{-1}).
\end{align}
-\marginnote{The equality in \cref{eq:laplace-normal} only holds when the precision matrix \(\bm \Lambda\) is positive definite.}We have now shown that approximating the log-posterior with a second degree Taylor expansion around the \(\bm\theta_{\textsc{map}}\) corresponds to approximating the posterior with a Gaussian distribution given by \(\normal(\bm \theta_{\textsc{map}}, \bm \Lambda^{-1})\) where \(\bm \Lambda = - \left.\nabla^2_{\bm \theta} \log p(\bm \theta \given \bm y) \right|_{\bm \theta_{\textsc{map}}}\).
+\marginnote{The equality in \cref{eq:laplace-normal} only holds when the precision matrix \(\bm \Lambda\) is positive definite.}
+We have now shown that approximating the log-posterior with a second degree Taylor expansion around the \(\bm\theta_{\textsc{map}}\) corresponds to approximating the posterior with a Gaussian distribution given by \(\normal(\bm \theta \given \bm \theta_{\textsc{map}}, \bm \Lambda^{-1})\) where \(\bm \Lambda = - \left.\nabla^2_{\bm \theta} \log p(\bm \theta \given \bm y) \right|_{\bm \theta_{\textsc{map}}}\).
In practice, we can determine the Hessian of the loss with respect to the parameters, where the loss is given by the negative log-posterior.
-In this case, we get \(\bm \Lambda = \left. - \nabla^2_{\bm \theta} \log p(\bm \theta \given \bm y) \right|_{\bm \theta_{\textsc{map}}} = \left. \nabla^2_{\bm \theta} \mathcal{L} \right|_{\bm \theta_{\textsc{map}}}\) and the posterior is approximated by \(\normal(\bm \theta_{\textsc{map}}, \bm \Lambda^{-1})\).
+In this case, we get \(\bm \Lambda = \left. - \nabla^2_{\bm \theta} \log p(\bm \theta \given \bm y) \right|_{\bm \theta_{\textsc{map}}} = \left. \nabla^2_{\bm \theta} \mathcal{L} \right|_{\bm \theta_{\textsc{map}}}\) and the posterior is approximated by \(\normal(\bm \theta \given \bm \theta_{\textsc{map}}, \bm \Lambda^{-1})\).
\cref{eq:negative-log-likelihood-ggn}
@@ -158,7 +179,7 @@ Furthermore, if all the neural network's second partial derivatives are continuo
In this case, the Hessian is generally dominated by the block-diagonal~\sidecite{martens2015optimizing}.
Typically, the loss function is considered as the negative log-likelihood or negative log-posterior.
-In this framing, if considering, e.g., the negative log-posterior loss from \cref{eq:negative-log-posterior}, we get ... from \cref{eq:negative-log-posterior-hessian}.\sidenote{For the common prior given by a zero-mean Gaussian \(\theta \sim \normal(\bm 0, \alpha^{-1} \identity)\), the log-prior regularisation term is given by \(\nabla^2_{\bm\theta} \log p(\bm \theta) = \alpha \identity\).}
+In this framing, if considering, e.g., the negative log-posterior loss from \cref{eq:negative-log-posterior}, we get ... from \cref{eq:negative-log-posterior-hessian}.\sidenote{For the common prior given by a zero-mean Gaussian \(\bm \theta \sim \normal(\bm 0, \alpha^{-1} \identity)\), the log-prior regularisation term is given by \(\nabla^2_{\bm\theta} \log p(\bm \theta) = \alpha \identity\).}
\begin{align}\label{eq:negative-log-posterior-hessian}
\nabla^2_{\bm\theta} \mathcal{L}(\bm \theta) ={} & \nabla^2_{\bm\theta} \left( \sum_i^N \log p(x_i \given \bm \theta) + \log p(\bm \theta) \right) \nonumber
\\ ={}& \nabla^2_{\bm\theta} \sum_i^N \log p(x_i \given \bm \theta) + \nabla^2_{\bm\theta} \log p(\bm \theta) \nonumber
@@ -188,7 +209,7 @@ This approximation is given by
%
% In this approximation, we have \(\nabla^2_{\bm \theta} f_{\bm \theta_0}(\bm x) = 0\).
In this way, we are linearising our neural network function.
-By combining \cref{eq:ggn-whole-hessian} and \cref{eq:nn-taylor} we then obtain the Gauss-Newton Hessian approximation,
+By combining \cref{eq:ggn-whole-hessian,eq:nn-taylor} we then obtain the Gauss-Newton Hessian approximation,
%
\begin{align*}
\nabla^2_{\bm\theta} \mathcal{L}
diff --git a/chapters/training.tex b/chapters/training.tex
@@ -3,8 +3,9 @@
\labch{training}
\section{Maximising the Evidence}
+\label{sec:maximising-the-evidence}
-According to Bayes' rule, the posterior probability for a given model is equal to
+As seen in \cref{eq:bayes-theorem}, the posterior probability for a given model is equal to
\begin{align}
p(\bm \theta \given \bm x) = \frac{p(\bm x \given \bm \theta) p(\bm \theta)}{p(\bm x)},
\end{align}
@@ -13,8 +14,7 @@ where \(p(\bm x)\) is known as the evidence. It's also known as the \emph{margin
p(\bm x) ={} & \int_{\bm \theta} p(\bm \theta \given \bm x)\,d\bm \theta.
\end{align}
This term is typically intractable, since it involves integrating over the entire parameter space.
-However, certain
-For the Laplace approximation, we approximate the posterior as \(p(\bm \theta \given \bm x) \approx \normal(\bm \theta \given \bm \theta_{\textsc{map}}, \bm \Lambda^{-1})\), where \(\bm \Lambda = -\nabla^2_{\bm\theta} \mathcal L(\bm \theta_{\textsc{map}})\) is the Hessian of the negative log-posterior at the maximum a posteriori (MAP) estimate \(\bm \theta_{\textsc{map}}\).
+For the Laplace approximation, we approximate the posterior as \(p(\bm \theta \given \bm x) \approx \normal(\bm \theta \given \bm \theta_{\textsc{map}}, \bm \Lambda^{-1})\), where the posterior precision \(\bm \Lambda = -\nabla^2_{\bm\theta} \mathcal L(\bm \theta_{\textsc{map}})\) is the Hessian of the negative log-posterior at the maximum a posteriori (MAP) estimate \(\bm \theta_{\textsc{map}}\).
We thus obtain the approximate Laplace posterior
\begin{align}\label{eq:laplace-posterior}
\log p(\bm \theta \given \bm x) \stackrel{\textsc{la}}{\approx}{} & \log p(\bm \theta_{\textsc{map}} \given \bm x) + \nabla_{\bm \theta} \log p(\bm \theta_{\textsc{map}} \given \bm x) \cdot (\bm \theta - \bm \theta_{\textsc{map}})
@@ -69,14 +69,14 @@ While these terms do not affect the optimisation problem, they are useful for co
\\ 3. \(\alpha / 2 \norm{\bm \theta}^2\)
\end{tabular}
-\sidetextcite{fong2020marginal} argue that maximising the marginal (and consequently the log-marginal) is equivalent to performing k-fold cross validation for all values \(k = 1, \ldots, \infty\) and choosing the model with the highest average posterior probability (across each of \(k\) folds and across all values of \(k\)).
+\sidetextcite{fong2020marginal} show that maximising the marginal (and consequently the log-marginal) is equivalent to performing k-fold cross validation for all values \(k = 1, \ldots, \infty\) and choosing the model with the highest average posterior probability (across each of \(k\) folds and across all values of \(k\)).
\section[The Determinant Lower Bound]{The Determinant Lower Bound}[The Lower Bound]
We compute an \emph{upper} bound on the log-determinant of the Hessian of the negative log-posterior \(\bm \Lambda \in \reals^{D \times D}\), i.e.,
\begin{align}
\log \det \bm \Lambda \leq{} & B_D(\mu_1, \mu_2, \beta)
- \\ ={}& \begin{pmatrix} \log \beta & \log t \end{pmatrix}
+ \\ :={}& \begin{pmatrix} \log \beta & \log t \end{pmatrix}
\begin{pmatrix} \beta & t \\ \beta^2 & t^2 \end{pmatrix}^{-1}
\begin{pmatrix} \mu_1 \\ \mu_2 \end{pmatrix},
\end{align}
