commit b9dacadb234f4426c4b6364e5460d54e184988ef
parent eecccef8c324dfb57d3af7010f8863301b8a3b68
Author: Silas Brack <s174433@student.dtu.dk>
Date: Tue, 3 Jan 2023 01:38:41 +0100
Ran ltex.
Diffstat:
8 files changed, 32 insertions(+), 33 deletions(-)
diff --git a/chapters/acknowledgments.tex b/chapters/acknowledgments.tex
@@ -5,7 +5,7 @@
{
To S{\o}ren's group --- in particular, Marco --- for always making me feel like a part of the group, and like my contribution was valued.
- To Marco and Frederik for giving me the chance to dip my toes in ``real'' research in the rawest, realest way possible, with what I perceive to be all of the highs and the lows which constitute it.
+ To Marco and Frederik for giving me the chance to dip my toes in ``real'' research in the rawest, realest way possible, with what I perceive to be all the highs and the lows which constitute it.
To S{\o}ren, for always being up for a great chat, for giving me the opportunity to work on difficult and novel problems (such as preconditioning), and for approaching everything with kindness.
I don't think I could have chosen a better supervisor.
diff --git a/chapters/experiments.tex b/chapters/experiments.tex
@@ -30,7 +30,7 @@ This experiment will then answer two questions:
We then compose a function \((g \circ f) (\bm x) := \sum^N_{i=1} (2 x_i + 0.5)^2\) in terms of
\(f(\bm x) = 2 \bm x + 0.5, g(\bm x) := \sum^N_{i=1} x_i^2\).
-This is equivalent to the sum of square errors loss \(g\) for a simple linear model \(f\).
+This is equivalent to the sum-of-square-error loss \(g\) for a simple linear model \(f\).
Since \(f(x)\) is a linear function, the GGN matrix will exactly equal the actual Hessian.
We then initialise \(\bm x\) to a random vector and compute the matrix-vector product \(\bm K \bm v\) where \(\bm v\) is a random vector using three Hessian-vector products: the JAX efficient HVP and GVP functions and a manual Hessian instantiation and multiplication.
We then measure the total wall-clock time for computing this Hessian and multiplying it by \(\bm v\).
@@ -75,7 +75,7 @@ However, CIQ is typically extremely close to the exact solution for low values o
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.
-Root mean square error of my percentiles vs the real chi-sq percentiles? Feels a bit weird.
+Root-mean-square error of my percentiles vs the real chi-sq percentiles? Feels a bit weird.
As we vary \(q \in [1, \ldots, 20]\), how does performance vary?
What metrics can we use?
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 reparameterisation trick:
+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:
\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/introduction.tex b/chapters/introduction.tex
@@ -53,8 +53,7 @@ in modern optimizers, such as Adam. We use a sampling-based training
procedure, where a sample is first drawn from a Gaussian
weight-posterior, a gradient step is performed on this sampled neural
network, and the variance of the weight-posterior is updated with an
-approximate Hessian. Approximating the Hessian is the most time
-consuming and painful-to-engineer step of this training procedure. In
+approximate Hessian. Approximating the Hessian is the most time-consuming and painful-to-engineer step of this training procedure. In
this project, we tap into the potential of modern machine learning
frameworks to efficiently approximate the Hessian.
diff --git a/chapters/mixture.tex b/chapters/mixture.tex
@@ -17,7 +17,7 @@ where \(\bm \mu_i\) is computed by shifting the MAP solution based on ??? and \(
\bm \mu_i ={} & \bm \mu_{\textsc{map}} - \nabla_{\bm\theta} F_n (\bm \theta) \nabla_{\bm\theta}^2 F_n (\bm \theta)
\\ \bm \Sigma^{-1}_i ={}& - \nabla^2_{\bm \theta} \log p(\bm \theta \given x_i)
\end{align}
-The reparameterisation trick can then be performed using only a single observation to determine the covariance matrix
+The reparametrisation trick can then be performed using only a single observation to determine the covariance matrix
\begin{align}
\bm \varepsilon ={} & \bm \mu_i + \bm \Sigma^{-1/2}_i \bm \varepsilon_0,
\end{align}
diff --git a/chapters/preface.tex b/chapters/preface.tex
@@ -7,7 +7,7 @@
I suppose that, because of this, there is some kind of expectation of anticipation, of the sensation that a culmination/apogee of a.
I, not always being the biggest adherent to the idea of symbol-for-the-sake-of-symbol, wasn't sure what to expect of this experience
- The work in this thesis is the product of a five-month-long research project conducted in the Section for Cognitive Systems of DTU Compute, during which I helped develop some of the ideas in the paper ``Large-Scale Laplace''(cite).
+ The work in this thesis is the product of a five-month-long research project conducted in the Section for Cognitive Systems of DTU Compute, during which I helped develop some ideas in the paper ``Large-Scale Laplace''(cite).
}
\begin{flushright}
diff --git a/chapters/sampling.tex b/chapters/sampling.tex
@@ -5,20 +5,20 @@
\section{Introduction}
In many methods, such as the Laplace approximation (\cref{sec:laplace}), we need to sample from a multivariate normal distribution parameterised by a mean vector \(\bm \mu\) and a covariance matrix \(\bm \Sigma\).
-To sample from a Gaussian distribution, you can apply the \emph{reparameterisation trick} to compute a vector of samples as
+To sample from a Gaussian distribution, you can apply the \emph{reparametrisation trick} to compute a vector of samples as
\begin{align}\label{eq:sample-reparam}
\bm \varepsilon ={} & \bm \mu + \bm \Sigma^{1/2} \bm \varepsilon_0,
\end{align}
where \(\bm \varepsilon_0\) are samples drawn from a standard normal distribution, i.e., \(\bm \varepsilon_0 \sim \normal(\bm 0, \identity)\).
Using the Laplace approximation, we estimate the precision matrix as the Hessian matrix of the loss with respect to the parameters, as discussed in \cref{sec:laplace}.
-Because of this, the covariance matrix is given by the inverse of this Hessian and we thus sample as per
+Because of this, the covariance matrix is given by the inverse of this Hessian, and we thus sample as per
\begin{align}\label{eq:sample-laplace}
\bm \varepsilon = \bm \mu + \bm \Lambda^{-1/2} \bm \varepsilon_0,
\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\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}.
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]
@@ -111,7 +111,7 @@ Since \(\sigma\) is then imaginary (its real part is zero), \(\tau := \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\).
-Since we will be sampling a fixed number of quadrature points from the circle, to improve the closeness of the approximation of the integral, it would be beneficial to choose a closed circle which intersects the real axis \emph{at} the minimum and maximum eigenvalues, since this will avoid sampling quadrature points which lie outside of the spectrum of \(\bm K\).
+Since we will be sampling a fixed number of quadrature points from the circle, to improve the closeness of the approximation of the integral, it would be beneficial to choose a closed circle which intersects the real axis \emph{at} the minimum and maximum eigenvalues, since this will avoid sampling quadrature points which lie outside the spectrum of \(\bm K\).
% ------------------------------------------
% JACOBI ELLIPTICAL FUNCTIONS
@@ -213,11 +213,11 @@ 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\}
\end{align}
-where \(\{s_1, \ldots, s_p\}, p \leq N \cdot k\) are the nonzero 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 J_i\T H J_i\).
Notice that if a preconditioner approximates all eigenvalues and eigenvectors of a
Something something spectrum, shape of spectrum, how steep is it?
-The more values you add, the steeper the spectrum becomes (even the nonzero one)
+The more values you add, the steeper the spectrum becomes (even the non-zero one)
\begin{enumerate}
\item Scaling
@@ -241,7 +241,7 @@ Thus, there is an intrinsic balance to be struck between a preconditioner which
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.
One example of such a preconditioner would be???
Depending on the rank of the preconditioner, we can then eliminate a number of the largest / lowest eigenvalues of \(H\) equal to 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.
+If we have a rank-one preconditioner whose only eigenvector is the largest eigenvector of \(H\), 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 \(H\) has a very
@@ -249,7 +249,7 @@ If \(H\) has a very
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}.
-As such, if the preconditioner-vector product is similarly or more expensive to compute compared to the matrix-vector product, then preconditoineing will ....
+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.
\subsection[Adam-esque fully linear GGN]{Adam-esque fully linear GGN}[Fully Linear GGN]
@@ -263,9 +263,9 @@ If we calculate the GGN approximation of the Hessian by linearizing over the who
DO I NEED TO PROVE THIS?
The problem with using this approximation is that the rank of \(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 highest precision (empirically, it seems that \emph{this is not the case}).
+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.
-If the rank of our precision matrix is greater than one, then this value seems likely to be similarly large to the originally largest eigenvalue.
+If the rank of our precision matrix is greater than one, then this value seems likely to be similarly large to originally the largest eigenvalue.
Therefore, the conditioning number of our problem will not be significantly reduced.
As such, this preconditioner is not worth pursuing.
@@ -279,8 +279,8 @@ This would have the advantage of
% \Span_\lambda\left( J\T J + \alpha \identity \right)
\end{align}
-If the eigenvectors of \(J_i\T J_i\) are all (approximately) the same, then the \(\textstyle{\sum^N_i} J_i\T J_i\) term has \(k\) nonzero eigenvalues (the same as \(J\T J\)).
-If these eigenvectors are all different (orthogonal?), then the aforementioned term will have \(B \cdot k\) nonzero eigenvalues.
+If the eigenvectors of \(J_i\T J_i\) are all (approximately) the same, then the \(\textstyle{\sum^N_i} J_i\T J_i\) term has \(k\) non-zero eigenvalues (the same as \(J\T J\)).
+If these eigenvectors are all different (orthogonal?), then the aforementioned term will have \(B \cdot k\) non-zero eigenvalues.
Alternatively, we can also instead estimate the eigenvectors and eigenvalues of \(\textstyle{N / B \sum^B_i} J_i\T J_i\).
This would \emph{hopefully} have an eigenbasis that is closer to that of \(\textstyle{\sum^N_i} J_i\T J_i\) without needing to calculate the precision matrix over the whole dataset (we may need to guarantee all classes are represented in \(x_i, i \sim \mathcal{B}\)).
@@ -291,7 +291,7 @@ Now, dealing with poor conditioning only in the inversion of our preconditioner
\[\textstyle{\sum^N_i} s_1 v_i / \alpha \leq \kappa \leq N s_1 / \alpha\]
TALK ABOUT PRECONDITIONER WITH ALL OBSERVATIONS BUT DIFFERENT ALPHA
-NOTE THAT THIS POSES ANOTHE ISSUE: THE FACT THAT THE PRECONDITIONER BECOMES WAY MORE EXPENSIVE ON ALL OBSERVATIONS
+NOTE THAT THIS POSES ANOTHER ISSUE: THE FACT THAT THE PRECONDITIONER BECOMES WAY MORE EXPENSIVE ON ALL OBSERVATIONS
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.
@@ -313,7 +313,7 @@ We can then use this approximation as a preconditioner by computing its inverse
\\ P^{-1} ={}& \alpha^{-1} \identity - \alpha^{-2} L\T \left( \identity + \alpha^{-1} L L\T \right)^{-1} L
\end{align}
we then obtain the preconditioner-vector and inverse-preconditioner-vector products as in Equation~\ref{eq:pivoted-cholesky-preconditioner}.
-Notice that the ``scariest'' computation in the inverse preconditioner is inverting a \(k \times k \) matrix, where \(k\) is the rank of our pivoted Cholesky factorization (which we choose).
+Notice that 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 )
@@ -331,11 +331,11 @@ In terms of eigenvalues and eigenvectors,
Assume that the pivoted Cholesky factorisation finds the \(k\) largest eigenvectors.
If the fall-off in the eigenspectrum of this factorisation is steep enough, then the largest eigenvector \(s_1\) will be significantly larger than the smallest eigenvector \(s_k\), especially if the value of \(B\) is large enough that the
-Unfortunately, if 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 J_i\T 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, in general, possible to efficiently determine the diagonal of a matrix via matrix-vector products.
-Often,it is necessary to multiply a vector of zeros with a value of one at the index at which to find the diagonal.
+Often, it is necessary to multiply a vector of zeros with a value of one at the index at which to find the diagonal.
Alternatively, it is possible to approximate the diagonal of the GGN matrix as the element-wise product of the gradient with itself, similarly to the preconditioner described in \cref{ssec:adam-preconditioner} and the GGN approximation computed in \cref{eq:fully-linear-ggn}.
diff --git a/chapters/theory.tex b/chapters/theory.tex
@@ -58,10 +58,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 of the neural network's second partial derivatives are continuous, then the Hessian is symmetric and positive definite.
+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}.
-Typically, the loss functions is considered as the negative log-likelihood or negative log-posterior.
+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\).}
\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
@@ -120,12 +120,12 @@ For log-likelihood loss, we get
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}.
-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 some given model input \(x\).
+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.
Another choice is to linearise the whole likelihood function with respect to the model parameters.
-This yields \(g(y) = -\log(y)\) and thus we obtain
+This yields \(g(y) = -\log(y)\), and thus we obtain
\begin{align}\label{eq:fully-linear-ggn}
\nabla^2_{\bm \theta} \mathcal{L}
={} & \nabla f \T \cdot \nabla^2 g \cdot \nabla f \nonumber
@@ -139,7 +139,7 @@ This yields \(g(y) = -\log(y)\) and thus we obtain
The diagonal of this matrix is given by the element-wise square of the gradient \( \nabla_{\bm\theta} \mathcal{L} \odot \nabla_{\bm\theta} \mathcal{L} \).
Note that the choice of application of the generalised Gauss-Newton approximation is not based on having a reduced computational cost --- both the Hessian and the GGN approximation require the same number of forward and backward passes.
-However, the GGN approximation is more better behaved than the exact Hessian, thereby making it more suitable for practical applications.
+However, the GGN approximation is better behaved than the exact Hessian, thereby making it more suitable for practical applications.
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.
@@ -151,7 +151,7 @@ A GGN approximation is any approximation which fulfils \cref{eq:ggn-whole-hessia
\section{Fisher Information}\label{sec:fisher-information}
% https://agustinus.kristia.de/techblog/2018/03/11/fisher-information/
-Consider the gradients of the log-likelihood optimization objective, given by the expression \(\nabla_{\bm\theta} \log\lik\).
+Consider the gradients of the log-likelihood optimisation objective, given by the expression \(\nabla_{\bm\theta} \log\lik\).
We can show that the expectation of these gradients is zero:
\begin{align*}
\E{\lik}{\nabla_{\bm\theta} \log\lik} = & {} \E{\lik}{\nabla \log \lik}
@@ -163,7 +163,7 @@ We can show that the expectation of these gradients is zero:
\\ =&{} 0
\end{align*}
-Let us now analyze the covariance of these gradients.
+Let us now analyse the covariance of these gradients.
It is given by
%
\begin{align}\label{eq:gradient-covariance}
@@ -184,7 +184,7 @@ It can, however, be approximated empirically by using samples from an empirical
%
Note, however, that this matrix is distinct from the Fisher information matrix~\cite{martens2020new, kunstner2019limitations}, since
-Another nontrivial property of the Fisher information matrix is that it is equivalent to the negative expectation of the Hessian of the log-likelihood.
+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.
%
\begin{align*}
H_{\bm\theta} \left(\log \lik\right) = & {} \nabla^2_{\bm\theta} \log \lik
@@ -211,7 +211,7 @@ LIMITATIONS OF FISHER INFORMATION (MAYBE ESPECIALLY EMPIRICAL?)
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\) parameterized by \(\bm\theta\) by following the direction of steepest descent in the Fisher information metric.
+\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.
It can be interpreted as
% The algorithm can be seen in \cref{natural-gradient}.
%
