commit bdb3260af991245e3e6f37b044dbdaddac400edb
parent b9dacadb234f4426c4b6364e5460d54e184988ef
Author: Silas Brack <s174433@student.dtu.dk>
Date: Tue, 3 Jan 2023 02:08:17 +0100
Fixed citations.
Diffstat:
5 files changed, 38 insertions(+), 38 deletions(-)
diff --git a/chapters/experiments.tex b/chapters/experiments.tex
@@ -71,7 +71,7 @@ This is in line with the theory, as noted in \cref{ssec:practical-ggn} (SHOULD I
\section{Sampling Ablation}
The contour integral quadrature method (CIQ) is a Krylov subspace method which converges to the exact solution for \(\lim_{Q \to \infty} \bm s^{\mathrm{CIQ}}_{Q} = \bm K^{-1/2}\).
-However, CIQ is typically extremely close to the exact solution for low values of \(Q\)\sidenote{\textcite{pleiss2020fast} find that \(Q=20\) is approximately sufficient for most problems.}.
+However, CIQ is typically extremely 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.
Need to define a metric which determines if sampling was successful.
diff --git a/chapters/higham.tex b/chapters/higham.tex
@@ -3,7 +3,7 @@ As per the Laplace approximation, we compute its posterior precision as:
\begin{align}
\Sigma^{-1} = J\T J + \alpha \identity \label{eq:ggn-hessian}.
\end{align}
-We can then sample from the Laplace posterior by taking the inverse square root of our precision matrix\cite{fasi2022computing} and using the reparametrisation trick:
+We can then sample from the Laplace posterior by taking the inverse square root of our precision matrix~\sidecite{fasi2022computing} and using the reparametrisation trick:
\begin{align}
\Sigma^{1/2} ={} & \alpha^{-1/2} \identity + J\T \left(J J\T\right)^{-1} \left(\left(J J\T + \alpha \identity\right)^{-1/2} - \alpha^{-1/2} \identity \right) J,\nonumber
\\ \Sigma^{1/2} \epsilon ={}& \alpha^{-1/2} \epsilon + J\T \left(J J\T\right)^{-1} \left(\left(J J\T + \alpha \identity\right)^{-1/2} J \epsilon - \alpha^{-1/2} J \epsilon \right), \label{eq:higham-samples}
diff --git a/chapters/sampling.tex b/chapters/sampling.tex
@@ -18,12 +18,12 @@ Because of this, the covariance matrix is given by the inverse of this Hessian,
\end{align}
where \(\Lambda = \nabla^2_{\bm\theta} \log p(\bm \theta \given \bm x)\).
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\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 A\) (a rotation matrix), \(\bm A\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 A\T \bm \varepsilon \sim \normal(\bm 0, \identity)\).}, in this special case, we simply need to obtain a matrix which returns the original matrix \emph{up to a rotation}.
+However, since a Gaussian distribution is invariant to rotation, in this special case, we simply 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 A\) (a rotation matrix), \(\bm A\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 A\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.
\section[Contour Integral Quadrature]{Contour Integral Quadrature}[Contour Integral Quad]
-\textcite{pleiss2020fast} propose one such algorithm.
+\sidetextcite{pleiss2020fast} propose one such algorithm.
Contour integral quadrature (CIQ) attempts to solve the problem of computing the inverse-square root of a matrix, up to a rotation, by approximating the expression \(f(\bm K) \bm \varepsilon_0 = \bm K^{-1/2} \bm \varepsilon_0\) where \(\bm K\) is a positive semi-definite matrix using Cauchy's integral formula.
To do this, it exploits the Lanczos algorithm, Gaussian quadrature, and the msMINRES algorithm, methods which only require the computation of matrix-vector products \(\bm K \bm v\).
There are then essentially three steps to CIQ, namely:
@@ -75,11 +75,11 @@ From here, it is then trivial to also compute the product of a matrix square roo
% ------------------------------------------
The first step therefore invokes Cauchy's integral formula.
-Cauchy's integral formula is a central theorem in complex analysis which states that, for a holomorphic function \(f\)\marginnote{Note that the square root function is holomorphic for all positive definite real matrices.}, for a closed circular contour \(\Gamma\) in the complex plane (see \cref{fig:cauchy-integral-formula}) which encloses the eigenvalues of \(\bm K\), \(f(\bm K)\) can be approximated as
+Cauchy's integral formula is a central theorem in complex analysis which states that, for a holomorphic function \(f\), for a closed circular contour \(\Gamma\) in the complex plane (see \cref{fig:cauchy-integral-formula}) which encloses the eigenvalues of \(\bm K\), \(f(\bm K)\) can be approximated as\marginnote{Note that the square root function is holomorphic for all positive definite real matrices.}
\begin{align}\label{eq:cauchy-integral-formula}
f(\bm K) ={} & \frac{1}{2 \pi i} \oint_\Gamma f(\tau) \left(\tau \identity - \bm K\right)^{-1}\,d\tau.
\end{align}
-Since, in our case, the eigenvalues of \(\bm K\) will all be real-valued\marginnote{If a real-valued matrix is symmetric, its eigenvalues are real. In our case, \(\bm K\) will be the precision matrix obtained from the Laplace approximation.}, a circle centred on the real axis which encloses the minimum and maximum eigenvalues of \(\bm K\) will suffice.
+Since, in our case, the eigenvalues of \(\bm K\) will all be real-valued, a circle centred on the real axis which encloses the minimum and maximum eigenvalues of \(\bm K\) will suffice.\marginnote{If a real-valued matrix is symmetric, its eigenvalues are real. In our case, \(\bm K\) will be the precision matrix obtained from the Laplace approximation.}
% ------------------------------------------
% QUADRATURE
@@ -106,7 +106,7 @@ We can then apply this to the function \(f(\bm K) = \bm K^{-1/2}\), apply a chan
\end{align}
where \(\sigma\) is obtained from \(\tau\) as a change of variable and \(\Tilde{w}_1, \ldots, \Tilde{w}_Q\) are the \(Q\) quadrature weights.
We can then solve the \(Q\) systems of equations \(\left( \sigma^2_1 \identity - \bm K \right) \bm c_1 = \bm v, \ldots, \left( \sigma^2_Q \identity - \bm K \right) \bm c_Q = \bm v\) to obtain the solutions \(\bm c_q = \left(\sigma^2_q \identity - \bm K \right)^{-1} \bm v\).
-Note that, since the integrand is symmetric with respect to the real axis\sidenote{This is due to the fact that the circle with its centre on the \(x\)-axis is symmetric with respect to said \(x\)-axis.}, only the imaginary component of \(\Gamma_\sigma\) is required.
+Note that, since the integrand is symmetric with respect to the real axis, only the imaginary component of \(\Gamma_\sigma\) is required.\sidenote{This is due to the fact that the circle with its centre on the \(x\)-axis is symmetric with respect to said \(x\)-axis.}
Since \(\sigma\) is then imaginary (its real part is zero), \(\tau := \sigma^2\) is real such that \(\tau = - \Im(\sigma)^2\).
As such, since \(\Im(\sigma)^2\) must be positive, \(\tau\) must then be real-valued and negative.\marginnote{For simplicity, since \(w_q\) happens to also be real-valued and negative, from now on we will redefine the variables \(w_q \gets -\Tilde{w}_q\) and \(t_q \gets -\sigma^2_q\).}
In practice, we incorporate the \(1 / (\pi i) \) term into the quadrature weights and, with the quadrature points \(t_1, \ldots, t_Q\), compute \(\bm K^{-1/2} \bm v \approx \sum^Q_{q=1} w_q \bm c_q\).
@@ -118,7 +118,7 @@ Since we will be sampling a fixed number of quadrature points from the circle, t
% ------------------------------------------
We could sample these quadrature points from the circle uniformly, which corresponds to the regular trapezoid quadrature rule.
-However, the convergence for this quadrature is linear with regard to the condition number of \(\bm K\) \sidecite{pleiss2020fast,hale2008computing}.
+However, the convergence for this quadrature is linear with regard to the condition number of \(\bm K\)~\sidecite{pleiss2020fast,hale2008computing}.
Since the precision matrix computed from the Laplace approximation as the GGN matrix can often have a very low rank and poor conditioning, uniform sampling would then require a very large number of quadrature points \(Q\) and would therefore be inadequate.
Since our GGN precision matrix \(\bm \Lambda = \sum^N_{i=1} \bm J_i\T \bm H_i \bm J_i\) can be low-rank, we would instead prefer to oversample quadrature points which lie close to the minimum eigenvalue.
@@ -179,7 +179,7 @@ We can thereby now compute the quadrature weights and shifts, as seen in \cref{a
% \end{algorithm}
The last step involves solving the \(Q\) linear systems defined in \cref{eq:ciq} to determine the solves \(\bm c_q\).
-As mentioned above, this is done by using msMINRES (\cref{ssec:msminres,alg:msminres}) as \texttt{msMINRES(-K, b, t\_q, P, J)}\sidenote{Note that we negate \(\bm K\), since msMINRES computes \(\left(t_q \identity + \bm K \right)^{-1} \bm v\), while we want to compute \(\left(t_q \identity - \bm K \right)^{-1} \bm v\). Since we can only access the matrix-vector product, we simply use \(- \bm K \bm v\) whenever the matrix-vector product is called.}.
+As mentioned above, this is done by using msMINRES (\cref{ssec:msminres,alg:msminres}) as \texttt{msMINRES(-K, b, t\_q, P, J)}.\sidenote{Note that we negate \(\bm K\), since msMINRES computes \(\left(t_q \identity + \bm K \right)^{-1} \bm v\), while we want to compute \(\left(t_q \identity - \bm K \right)^{-1} \bm v\). Since we can only access the matrix-vector product, we simply use \(- \bm K \bm v\) whenever the matrix-vector product is called.}
This last step is then sensitive to preconditioning
@@ -187,7 +187,7 @@ Mention computational bottleneck being performing the msMINRES for the shifted s
\section{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{Recall that 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.}.
+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{Recall that 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}
@@ -195,7 +195,7 @@ In the case of the Laplace approximation, as per \cref{eq:ggn-laplace}, we have
\\ \bm \theta \sim \normal(\bm \mu, \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\sidenote{With rank \(o\), where \(o\) is the number of neural network outputs.} positive semi-definite matrix, it has \(o\) non-zero eigenvalues.
+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\).
@@ -248,7 +248,7 @@ If \(H\) has a very
% 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\sidenote{Because of this, to avoid the preconditioner becoming a performance bottleneck, the speed of the preconditioner-vector products should ideally be significantly greater than that of the matrix-vector products.} \textbf{check this}, as seen in \cref{alg:ciq}.
+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 ideally 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.
@@ -295,7 +295,7 @@ NOTE THAT THIS POSES ANOTHER ISSUE: THE FACT THAT THE PRECONDITIONER BECOMES WAY
As mentioned before, each msMINRES iteration performs one matrix-vector product and one preconditioner-vector product.
However, even though the preconditioner is constructed in nearly the same way as the GGN matrix, the preconditioner is significantly cheaper to compute.
-For the subsampled GGN preconditioner, the preconditioner-vector product will then cost \(B / N \cdot t_{\mathrm{GGN}}\), where \(t_{\mathrm{GGN}}\) is the computational cost of performing the GGN-vector product\sidenote{Note that this analysis assumes that the computation of the GGN-vector product scales linearly in the number of observations used in the GGN matrix.}.
+For the subsampled GGN preconditioner, the preconditioner-vector product will then cost \(B / N \cdot t_{\mathrm{GGN}}\), where \(t_{\mathrm{GGN}}\) is the computational cost of performing the GGN-vector product.\sidenote{Note that this analysis assumes that the computation of the GGN-vector product scales linearly in the number of observations used in the GGN matrix.}
Since \(N \gg B\), then the cost of computing the preconditioner will be negligible relative to that of the GGN-vector product.
\subsection{Pivoted Cholesky}
@@ -320,7 +320,7 @@ Notice that the ``scariest'' computation in the inverse preconditioner is invert
\end{align}
We can precompute \(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 \(L \in \reals^{k \times n}\) and the aforementioned \(\reals^{n \times k}\) matrix\sidenote{Note that we can also precompute \(( \identity + \alpha^{-1} L L\T )^{-1} \in \reals^{k \times k}\), which means we would only need to store \(L \in \reals^{k \times n}\) and the aforementioned \(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.}.
+This would then require storing \(L \in \reals^{k \times n}\) and the aforementioned \(\reals^{n \times k}\) matrix.\sidenote{Note that we can also precompute \(( \identity + \alpha^{-1} L L\T )^{-1} \in \reals^{k \times k}\), which means we would only need to store \(L \in \reals^{k \times n}\) and the aforementioned \(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.}
ADD SOMEWHERE:
As we sum over multiple low-rank matrices, we will likely find that the minimum non-zero eigenvalue gets progressively smaller relative to the largest.
@@ -364,7 +364,7 @@ However, due to the speed, effectiveness, and theoretical properties of other pr
\subsection[A Recap on Preconditioning]{A Recap on Preconditioning}[A Recap]
-In general, it is worth noting the flaw with scaling the term in the \(\sum J\T H J\) term\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.} preconditioner to account for the summing over a Hessian of a limited batch size
+In general, it is worth noting the flaw with 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.}
\begin{enumerate}
\item Never scale your preconditioner up
diff --git a/chapters/theory.tex b/chapters/theory.tex
@@ -9,7 +9,7 @@ Our goal is to find the optimal parameters \(\bm\theta^*\) which best model the
From a frequentist perspective, we can define a loss function \(\mathcal{L}(\bm\theta): \reals^d \rightarrow \reals\) such that the optimal parameters \(\bm\theta^*\) minimise this loss function.
Often the loss function is defined 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.}.
+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)\).
@@ -59,10 +59,10 @@ Then the Hessian can be written as \(J_{\bm \theta}(\nabla_{\bm\theta} \mathcal{
We are interested in the Hessian of the loss function \(\mathcal{L}(\bm\theta)\) with respect to the parameters \(\bm\theta\).
For a neural network with \(d\) parameters, the Hessian is therefore a \(d \times d\) square matrix.
Furthermore, if all the neural network's second partial derivatives are continuous, then the Hessian is symmetric and positive definite.
-In this case, the Hessian is generally dominated by the block-diagonal \cite{martens2015optimizing}.
+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\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 \(\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
@@ -119,7 +119,7 @@ For log-likelihood loss, we get
\label{ssec:practical-ggn}
Note that the generalised Gauss-Newton approximation requires that the function for which we would like to compute the Hessian, \(h(x)\), can be framed as a composition of two functions \(h(x) = g(f(x))\).
-However, there are many possible choices of functions \(f\) and \(g\) which, when composed, yield \(h\), and these different compositions can yield wildly varying performance for the GGN approximation~\cite{kunstner2019limitations}.
+However, there are many possible choices of functions \(f\) and \(g\) which, when composed, yield \(h\), and these different compositions can yield wildly varying performance for the GGN approximation~\sidecite{kunstner2019limitations}.
Given that, in our case, \(h\) is the loss function with respect to the model parameters of a neural network, a natural choice of functions (and the choice we have assumed above) are \(g: \reals^o \rightarrow \reals\) and \(f: \reals^d \rightarrow \reals^o\) where \(g\) is the loss function with respect to the model output and \(f\) is a function which maps the model parameters to a model output, for a given model input \(x\).
This construction allows us to capture the curvature of the loss with respect to the model output, which can often be easily computed in a simple closed form, without requiring the computation of the curvature of the model function with respect to its parameters; this therefore corresponds to linearising the model function with a Taylor expansion as in \cref{eq:nn-taylor}.
However, this is not the only possible construction.
@@ -143,10 +143,10 @@ However, the GGN approximation is better behaved than the exact Hessian, thereby
This is because the GGN approximation is designed to be positive semi-definite, whereas the exact Hessian can be vulnerable to negative curvature.
This feature can limit its use for optimisation algorithms that require a positive definite Hessian, such as conjugate gradient descent, or the Laplace approximation, which employs the inverse of the Hessian as the covariance matrix of a Gaussian distribution, thereby assuming a positive semi-definite covariance matrix.
-The Levenberg-Marquardt heuristic~\cite{pascanu2013revisiting}
+The Levenberg-Marquardt heuristic~\sidecite{pascanu2013revisiting}
A GGN approximation is any approximation which fulfils \cref{eq:ggn-whole-hessian}.
-\cite{kunstner2019limitations}
+\sidecite{kunstner2019limitations}
\section{Fisher Information}\label{sec:fisher-information}
% https://agustinus.kristia.de/techblog/2018/03/11/fisher-information/
@@ -182,7 +182,7 @@ It can, however, be approximated empirically by using samples from an empirical
\mathcal{F} \approx \frac{1}{N} \sum^N_{n=1} \nabla \log p(y_n \given \bm \theta) \nabla \log p(y_n \given \bm \theta)\T.
\end{equation}
%
-Note, however, that this matrix is distinct from the Fisher information matrix~\cite{martens2020new, kunstner2019limitations}, since
+Note, however, that this matrix is distinct from the Fisher information matrix~\sidecite{martens2020new, kunstner2019limitations}, since
Another non-trivial property of the Fisher information matrix is that it is equivalent to the negative expectation of the Hessian of the log-likelihood.
%
@@ -196,22 +196,22 @@ Another non-trivial property of the Fisher information matrix is that it is equi
The Fisher information matrix can therefore be interpreted as the curvature of the log-likelihood.
RELATIONSHIP WITH GAUSS-NEWTON APPROXIMATION
-for common loss functions, the GGN is equivalent to the Fisher information matrix~\cite{martens2020new}
-\cite{kunstner2019limitations}
+for common loss functions, the GGN is equivalent to the Fisher information matrix~\sidecite{martens2020new}
+\sidecite{kunstner2019limitations}
RELATIONSHIP WITH KL DIVERGENCE
% https://agustinus.kristia.de/techblog/2018/03/14/natural-gradient/
% https://en.wikipedia.org/wiki/Fisher_information#Relation_to_relative_entropy
LIMITATIONS OF FISHER INFORMATION (MAYBE ESPECIALLY EMPIRICAL?)
-\cite{kunstner2019limitations}
+\sidecite{kunstner2019limitations}
\section{The Natural Gradient}
% https://agustinus.kristia.de/techblog/2018/03/14/natural-gradient/
Natural Gradient Descent is an approximate second-order optimisation method. It has an interpretation as optimizing over a Riemannian manifold using an intrinsic distance metric, which implies the updates are invariant to transformations such as whitening. By using the positive semi-definite (PSD) Gauss-Newton matrix to approximate the (possibly negative definite) Hessian, NGD can often work better than exact second-order methods.
-\textcite{amari1998natural} introduced the natural gradient as a way to optimize a function \(f\) parameterised by \(\bm\theta\) by following the direction of the steepest descent in the Fisher information metric.
+side{amari1998natural} introduced the natural gradient as a way to optimize a function \(f\) parameterised by \(\bm\theta\) by following the direction of the steepest descent in the Fisher information metric.
It can be interpreted as
% The algorithm can be seen in \cref{natural-gradient}.
%
@@ -234,12 +234,12 @@ This therefore limits the applicability of the natural gradient.
Instead, it is often approximated by the diagonal of the Hessian.
% This brings us to the Adam optimizer, which is a popular choice for optimizing deep neural networks.
-\cite{kingma2014adam}
-\cite{martens2015optimizing}
-\cite{botev2017practical}
+\sidecite{kingma2014adam}
+\sidecite{martens2015optimizing}
+\sidecite{botev2017practical}
-For a recent, detailed discussion of the natural gradient, see \cite{martens2020new}.
-\cite{wu2019logan}
+For a recent, detailed discussion of the natural gradient, see~\sidetextcite{martens2020new}.
+\sidecite{wu2019logan}
\section[Inverting a low-rank matrix]{Inverting a low-rank matrix}[Woodbury-Morrison]
\label{sec:woodbury}
@@ -281,9 +281,9 @@ Because of this, matrix-vector products alone are sufficient to construct a Kryl
Krylov subspace methods leverage this fact to estimate functions of matrices multiplied with vectors, \(f(\bm K) \bm b\).
These methods therefore identify approximate solutions which are contained within the \(j\)\textsuperscript{th} Krylov subspace, i.e., \(f(\bm K) \bm b \approx \bm c_j \in \mathcal{K}_j( \bm K, \bm b)\).
-These methods are always exact after \(N\) iterations, however are often almost exact after \(J \ll N\) iterations\cite{pleiss2020fast}.
+These methods are always exact after \(N\) iterations, however are often almost exact after \(J \ll N\) iterations~\sidecite{pleiss2020fast}.
One major advantage Krylov subspace methods have is that they only require computation of matrix-vector products \(\bm K \bm b\), thereby not requiring direct instantiation of the matrix \(\bm K\).
-As such, if allocation of memory to store a matrix is infeasible, but the computational graph for this matrix can be determined, then methods such as the conjugate gradient method or the minimum residual method can be used to approximate the inverse of this matrix~\cite{giordano2018covariances, nocedal1999numerical}.
+As such, if allocation of memory to store a matrix is infeasible, but the computational graph for this matrix can be determined, then methods such as the conjugate gradient method or the minimum residual method can be used to approximate the inverse of this matrix~\sidecite{giordano2018covariances, nocedal1999numerical}.
\subsection{The Lanczos Algorithm}
@@ -293,7 +293,7 @@ The power iteration method involves iteratively computing the normalised product
\bm b_{j+1} = \frac{\bm K \bm b_{j}}{\norm{\bm K \bm b_{j}}},
\end{align}
for some number of iterations \(J\), such that \(\lim_{j \to \infty} \bm b_j\) equals the eigenvector corresponding to the largest eigenvalue of \(\bm K\).
-This thereby yields a set of vectors which is proportional to\sidenote{This is because the norm of the Lanczos vectors \(\zeta_j := \norm{\bm K \bm b_{j}}\) is a scalar, which means that \(\bm b_j = \bm K \bm b_{j-1} \zeta_{j-1}^{-1} = \zeta_{j-1}^{-1} \bm K \bm b_{j-1}\), so \(\bm b_j = Z \bm K \bm b_{j-1}\) where \(Z = \prod_{i=0}^{j-1} \zeta_{i}^{-1} \).} \(\{ \bm b_0, \bm K \bm b_0, \bm K^2 \bm b_0, \ldots, \bm K^{j-1} \bm b_0 \}\) which is clearly a spanning set for a Krylov subspace.
+This thereby yields a set of vectors which is proportional to \(\{ \bm b_0, \bm K \bm b_0, \bm K^2 \bm b_0, \ldots, \bm K^{j-1} \bm b_0 \}\) which is clearly a spanning set for a Krylov subspace.\sidenote{This is because the norm of the Lanczos vectors \(\zeta_j := \norm{\bm K \bm b_{j}}\) is a scalar, which means that \(\bm b_j = \bm K \bm b_{j-1} \zeta_{j-1}^{-1} = \zeta_{j-1}^{-1} \bm K \bm b_{j-1}\), so \(\bm b_j = Z \bm K \bm b_{j-1}\) where \(Z = \prod_{i=0}^{j-1} \zeta_{i}^{-1} \).}
However, the power iteration method only uses the last vector \(\bm b_J\), thereby effectively ``wasting'' the other vectors \(\bm b_1, \ldots, \bm b_{J-1}\).
To make use of the entire Krylov subspace \(\mathcal{K}_j(\bm K, \bm b_0)\), the Lanczos algorithm combines the power iteration method with the Gram–Schmidt process, along with the assumption that \(\bm K\) is Hermitian, to produce an orthonormal basis from the Krylov subspace.
This method then converts an eigendecomposition problem for \(\bm K\) into an eigendecomposition problem for a tridiagonal matrix of rank \(J\) whose maximum and minimum eigenvalues will be close to those of \(\bm K\).
@@ -304,7 +304,7 @@ Further adjustments to the Lanczos algorithm yield the conjugate gradient method
\begin{align*}
\bm x = \bm K^{-1} \bm b,
\end{align*}
-where \(\bm K\) is a positive semi-definite\marginnote{Note that a positive semi-definite matrix is square.} matrix.
+where \(\bm K\) is a positive semi-definite matrix.\marginnote{Note that a positive semi-definite matrix is square.}
The algorithm itself can be seen in \cref{alg:cg}.
% EXPLANATION
@@ -403,7 +403,7 @@ Complete elliptic integral of the first kind is given by
\begin{align}
\mathcal{K}(p) = \int_0^{\pi/2} \left( 1 - m \sin(t)^2 \right)^{-1/2}\,dt, \quad p = 1 - m,
\end{align}
-and is defined in the domain\sidenote{Note that if \(p > 1\), then what should be computed is instead the identity \(\mathcal{K}(p) = \mathcal{K}(1/p) / \sqrt{p}\).} \(0 < p \leq 1\).
+and is defined in the domain \(0 < p \leq 1\).\sidenote{Note that if \(p > 1\), then what should be computed is instead the identity \(\mathcal{K}(p) = \mathcal{K}(1/p) / \sqrt{p}\).}
This integral is approximated by
\begin{align}
\mathcal{K}(p) \approx P(p) - \log(p) Q(p)
diff --git a/chapters/training.tex b/chapters/training.tex
@@ -8,9 +8,9 @@ According to Bayes' rule, the posterior probability for a given model is equal t
\begin{align}
p(\theta \given x) = \frac{p(x \given \theta) p(\theta)}{p(x)},
\end{align}
-where \(p(x)\) is known as the evidence\sidenote{It's also known as the \emph{marginal likelihood}, since it's obtained by marginalising out \(\theta\) from the likelihood, as per \(p(x) = \int_\theta p(x \given \theta) p(\theta)\,d\theta\).}.
+where \(p(x)\) is known as the evidence.\sidenote{It's also known as the \emph{marginal likelihood}, since it's obtained by marginalising out \(\theta\) from the likelihood, as per \(p(x) = \int_\theta p(x \given \theta) p(\theta)\,d\theta\).}
-\textcite{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} 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\)).
\section[The Determinant Lower Bound]{The Determinant Lower Bound}[The Lower Bound]
