#  FSP-Laplace: Function-Space Priors for the Laplace Approximation in
Bayesian Deep Learning

Tristan Cinquin
University of Tübingen
<EMAIL_ADDRESS>
&Marvin Pförtner11footnotemark: 1
University of Tübingen
<EMAIL_ADDRESS>
&Vincent Fortuin
Helmholtz AI, TU Munich
<EMAIL_ADDRESS>
&Philipp Hennig
University of Tübingen
<EMAIL_ADDRESS>
&Robert Bamler
University of Tübingen
<EMAIL_ADDRESS>Equal contribution

###### Abstract

Laplace approximations are popular techniques for endowing deep networks with
epistemic uncertainty estimates as they can be applied without altering the
predictions of the neural network, and they scale to large models and
datasets. While the choice of prior strongly affects the resulting posterior
distribution, computational tractability and lack of interpretability of
weight space typically limit the Laplace approximation to isotropic Gaussian
priors, which are known to cause pathological behavior as depth increases. As
a remedy, we directly place a prior on function space. More precisely, since
Lebesgue densities do not exist on infinite-dimensional function spaces, we
have to recast training as finding the so-called weak mode of the posterior
measure under a Gaussian process (GP) prior restricted to the space of
functions representable by the neural network. Through the GP prior, one can
express structured and interpretable inductive biases, such as regularity or
periodicity, directly in function space, while still exploiting the implicit
inductive biases that allow deep networks to generalize. After model
linearization, the training objective induces a negative log-posterior density
to which we apply a Laplace approximation, leveraging highly scalable methods
from matrix-free linear algebra. Our method provides improved results where
prior knowledge is abundant, e.g., in many scientific inference tasks. At the
same time, it stays competitive for black-box regression and classification
tasks where neural networks typically excel.

##  1 Introduction

Neural networks (NN) have become the workhorse for many machine learning
tasks, but they do not quantify the uncertainty arising from data scarcity—the
_epistemic_ uncertainty. They are therefore incapable of providing a measure
of confidence along with their predictions, which is needed for safety-
critical applications, decision making, and scientific modeling [Effenberger
et al., , Indrayan, 2020, Tazi et al., 2023]. As a solution, Bayesian neural
networks (BNNs) cast training as approximating the Bayesian posterior
distribution over the weights (i.e., the distribution of the model parameters
given the training data), thus naturally capturing epistemic uncertainty.
Various methods exist for approximating the (intractable) posterior, either by
a set of samples (MCMC) or by a parametric distribution (e.g., variational
inference or Laplace approximation). While sampling methods are asymptotically
exact for large numbers of samples, they are typically expensive for large
models, whereas variational inference (VI) and the Laplace approximation are
popular because they are scalable. The Laplace approximation is a particularly
appealing method as it does not alter the maximum-a-posteriori (MAP)
predictions. Compared to VI, its differential nature allows the Laplace
approximation to scale to large datasets and models, and has been shown to
provide well-calibrated uncertainty estimates [Immer et al., 2021a, Daxberger
et al., 2022a].

##### The need for function-space priors in BNNs.

While the choice of prior strongly influences the uncertainty estimates
obtained from the posterior [Fortuin, 2022, Li et al., 2024], there are hardly
any methods for elicitation of informative priors in the literature on BNNs,
with the notable exceptions of Sun et al. [2019] and Tran et al. [2022]. This
is not just a conceptual issue. For instance, the default isotropic Gaussian
prior, commonly thought of as uninformative, actually carries incorrect
beliefs about the posterior (uni-modality, independent weights, etc.)
[Fortuin, 2022, Knoblauch et al., 2019] and has known pathologies in deep NNs
[Tran et al., 2022]. As network weights are not interpretable, formulating a
good prior on them is virtually impossible. Current methods work around this
issue by model selection, either through expensive cross-validation, or type-
II maximum likelihood estimation with a Laplace approximation under an
isotropic Gaussian prior, which can only be computed exactly for small
networks, and has known issues with normalization layers [Immer et al., 2021b,
Antorán et al., 2022]. As a solution, Sun et al. [2019] proposed a VI method
to specify priors directly on the function implemented by the BNN, with
promising results. Function-space priors incorporate interpretable knowledge
about the variance, regularity, periodicity, or length scale, building on the
extensive Gaussian Process (GP) literature [Rasmussen and Williams, 2006]. But
it turns out that the method by Sun et al. [2019] is difficult to use in
practice: the Kullback-Leibler (KL) divergence regularizing VI is infinite for
many function-space priors of interest [Burt et al., 2020], and estimating
finite variants of it is challenging [Ma and Hernández-Lobato, 2021].

Here, we address the lack of a procedure to specify informative priors in the
Laplace approximation by proposing a method to use interpretable GP priors.
Motivated by MAP estimation theory, our method first derives an objective
function that regularizes the neural network in function space and whose
minimizer corresponds to the MAP estimator of a GP on the space of functions
represented by the network. We then apply the Laplace approximation to the
log-posterior density induced by the objective, allowing us to effectively
incorporate beliefs from a GP prior into a deep network (see figure 1).
Efficient matrix-free methods from numerical linear algebra allow scaling to
large models and datasets. We demonstrate the effectiveness of our method by
showing improved results in cases where prior knowledge is available and
competitive performance for black-box regression and classification tasks
where neural networks typically excel. We make the following contributions:

  1. 1.

We propose a novel objective regularizing deep nets in function space using a
GP prior.

  2. 2.

We develop an efficient and scalable Laplace approximation that endows the
neural network with epistemic uncertainty reflecting the beliefs specified by
the GP prior.

  3. 3.

A range of experiments shows that our method improves performance on tasks
with prior knowledge in the form of a kernel, while showing competitive
performance for black-box regression, classification, out-of-distribution
detection, and Bayesian optimization tasks.

(a) FSP-Laplace posterior predictive of a BNN under different choices of GP
prior

(b) FSP-Laplace compared to related approximate inference methods and models

Figure 1: FSP-Laplace allows for efficient approximate Bayesian neural network
(BNN) inference under interpretable function space priors. Using our method,
it is possible to encode functional properties like smoothness, lengthscale,
or periodicity through a Gaussian process (GP) prior. The gray data points in
the plots are noisy observations of a periodic function.

##  2 Preliminaries: Laplace approximation in weight space

Consider a dataset ${\mathbb{D}}=({\bm{X}},{\bm{Y}})$ of inputs
${\bm{X}}=({\bm{x}}^{(i)})_{i=1}^{n}\in{\mathbb{X}}^{n}$ and targets
${\bm{Y}}=({\bm{y}}^{(i)})_{i=1}^{n}\in{\mathbb{Y}}^{n}$. We consider a model
of the data that is parameterized by a neural network
${\bm{f}}\colon{\mathbb{X}}\times{\mathbb{W}}\to{\mathbb{R}}^{d^{\prime}}$
with weights ${\bm{w}}\in{\mathbb{W}}\subset{\mathbb{R}}^{p}$, a likelihood
$p({\bm{Y}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{f}}({\bm{X}},{\bm{w}}))\coloneqq\prod_{i=1}^{n}p({\bm{y}}^{(i)}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{f}}({\bm{x}}^{(i)},{\bm{w}}))$,
and a prior $p({\bm{w}})$.

We seek the Bayesian posterior
$p({\bm{w}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\mathbb{D}})\propto
p({\bm{Y}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{f}}({\bm{X}},{\bm{w}}))\,p({\bm{w}})$.
As it is intractable in BNNs, approximate inference methods have been
developed. The Laplace approximation [MacKay, 1992a] approximates the
posterior by a Gaussian distribution
$q({\bm{\mathrm{w}}}={\bm{w}})={\operatorname{\mathcal{N}}\left({\bm{w}};{\bm{w}}^{\star},{\bm{\Lambda}}^{-1}\right)}$,
whose parameters are obtained by finding a MAP estimate
${\bm{w}}^{\star}\in\operatorname*{arg\,min}_{{\bm{w}}\in{\mathbb{W}}}R_{\text{WS}}({\bm{w}})$
of the weights, where

$R_{\text{WS}}({\bm{w}})\coloneqq-\log
p({\bm{w}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\mathbb{D}})=-\log
p({\bm{w}})-\sum_{i=1}^{n}\log
p({\bm{y}}^{(i)}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{f}}({\bm{x}}^{(i)},{\bm{w}}))+\text{const.}$
(2.1)

and computing the Hessian of the negative log-posterior
${\bm{\Lambda}}=-\left.{\bm{\mathrm{H}}}_{{\bm{w}}}\log
p({\bm{w}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\mathbb{D}})\right|_{{\bm{w}}={\bm{w}}^{\star}}\in{\mathbb{R}}^{p\times
p}$.

##### The linearized Laplace approximation

Computing ${\bm{\Lambda}}$ involves the Hessian of the neural network w.r.t.
its weights. This is generally associated with high computational cost. To
make the Laplace approximation scalable, it is common to linearize the network
around ${\bm{w}}^{\star}$ before computing ${\bm{\Lambda}}$ [Immer et al.,
2021c],

${\bm{f}}({\bm{x}},{\bm{w}})\approx{\bm{f}}({\bm{x}},{\bm{w}}^{\star})+{\bm{J}}_{{\bm{w}}^{\star}}({\bm{x}})({\bm{w}}-{\bm{w}}^{\star})\eqqcolon{\bm{f}}^{\text{lin}}({\bm{x}},{\bm{w}})$
(2.2)

with Jacobian
${\bm{J}}_{{\bm{w}}^{\star}}({\bm{x}})=\left.{\bm{\mathrm{D}}}_{{\bm{w}}}\,{\bm{f}}({\bm{x}},{\bm{w}})\right|_{{\bm{w}}={\bm{w}}^{\star}}$.
Thus, the approximate posterior precision matrix ${\bm{\Lambda}}$ is

$\displaystyle{\bm{\Lambda}}=-\left.{\bm{\mathrm{H}}}_{{\bm{w}}}\log
p({\bm{w}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\mathbb{D}})\right|_{{\bm{w}}={\bm{w}}^{\star}}$
$\displaystyle=-\left.{\bm{\mathrm{H}}}_{{\bm{w}}}\log
p({\bm{w}})\right|_{{\bm{w}}={\bm{w}}^{\star}}\underbracket{-\sum_{i=1}^{n}{\bm{J}}_{{\bm{w}}^{\star}}({\bm{x}}^{(i)})^{\top}{\bm{L}}^{(i)}_{{\bm{w}}^{\star}}{\bm{J}}_{{\bm{w}}^{\star}}({\bm{x}}^{(i)})}_{\\!\\!\\!\\!\\!=-\left.{\bm{\mathrm{H}}}_{{\bm{w}}}\log
p({\bm{Y}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{f}}^{\text{lin}}({\bm{X}},{\bm{w}}))\right|_{{\bm{w}}={\bm{w}}^{\star}}},$

where
${\bm{L}}^{(i)}_{{\bm{w}}^{\star}}=\left.{\bm{\mathrm{H}}}_{{\bm{f}}}\,p({\bm{y}}^{(i)}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{f}})\right|_{{\bm{f}}={\bm{f}}({\bm{x}}^{(i)},{\bm{w}}^{\star})}$.
The Hessian of the negative log-likelihood (NLL) under the linearized network
coincides with the _generalized Gauss-Newton matrix_ (GGN) of the NLL under
the non-linearized network. The GGN is equal to the Hessian if the model
perfectly fits the data, or if ${\bm{f}}({\bm{x}}^{(i)},\,\cdot\,)$ is linear
for all $i=1,\dotsc,n$. Crucially, the GGN is positive-(semi)definite even if
$-\left.{\bm{\mathrm{H}}}_{{\bm{w}}}\log
p({\bm{Y}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{f}}({\bm{X}},{\bm{w}}))\right|_{{\bm{w}}={\bm{w}}^{\star}}$
is not.

##  3 FSP-Laplace: Laplace approximation under function-space priors

A conventional Laplace approximation to BNNs requires a prior in weight space,
with the issues discussed in section 1. In the following, we present FSP-
Laplace, a method for computing Laplace approximations under interpretable GP
priors in function space. section 3.1 introduces an objective function that is
a log-density under local linearization. section 3.1 proposes a scalable
algorithm for the linearized Laplace approximation with a function-space prior
using matrix-free linear algebra.

###  3.1 Laplace approximations in function space

We motivate our Laplace approximation in function space through the lens of
MAP estimation in an (infinite-dimensional) function space (section 3.1) under
a GP prior. As Lebesgue densities do not exist in (infinite-dimensional)
function spaces, we cannot minimize equation 2.1 to find the MAP estimate. We
address this issue using a generalized notion of MAP estimation, resulting in
a minimizer of a regularized objective on the reproducing kernel Hilbert space
${\mathbb{H}}_{{\bm{\Sigma}}}$.

##### MAP estimation in neural networks under Gaussian process priors

The first step of the Laplace approximation is to find a MAP estimate of the
neural network weights ${\bm{w}}$, i.e., the minimizer of the negative log-
density function in Equation 2.1 w.r.t. the Lebesgue measure as a “neutral”
reference measure. In our method, we regularize the neural network in function
space using a $d^{\prime}$-output GP prior
${\bm{\mathrm{f}}}\sim{\operatorname{\mathcal{GP}}\left({\bm{\mu}},{\bm{\Sigma}}\right)}$
with index set ${\mathbb{X}}$. However, since (nonparametric) GPs take values
in infinite-dimensional (Banach) spaces of functions, no such reference
Lebesgue measure exists [Oxtoby, 1946]. Rather, the GP induces a prior
_measure_ $\mathrm{P}_{{\mathbb{B}}}$ on such spaces [Pförtner et al., 2022,
Section B.2]. We are thus interested in finding the “mode” of the posterior
measure $\mathrm{P}_{{\mathbb{B}}}^{{\bm{Y}}}$ under
$\mathrm{P}_{{\mathbb{B}}}$ defined by the Radon-Nikodym derivative
$\mathrm{P}_{{\mathbb{B}}}^{{\bm{Y}}}(\mathrm{d}{\bm{f}})\propto\exp\left(-\Phi^{\bm{Y}}(\mathrm{d}{\bm{f}})\right)\mathrm{P}_{{\mathbb{B}}}(\mathrm{d}{\bm{f}})$
where $\Phi^{\bm{Y}}$ is the _potential_ , in essence the negative log-
likelihood functional of the model. Similar to Bayes’ rule in finite
dimensions, this Radon-Nikodym derivative relates the prior measure to the
posterior measure by reweighting and can thus be thought of as a
generalization of a likelihood. Since there is no Lebesgue measure, the
(standard) mode is undefined, and we therefore use the framework in Lambley
[2023], which works with so-called _weak modes_.

######  Definition 1 (Weak Mode [see e.g., Lambley, 2023, Definition 2.1]).

Let ${\mathbb{B}}$ be a separable Banach space and let $\mathrm{P}$ be a
probability measure on $({\mathbb{B}},\mathcal{B}\left({\mathbb{B}}\right))$.
A _weak mode_ of $\mathrm{P}$ is any point
${\bm{f}}^{\star}\in\operatorname{supp}{\mathbb{B}}$ such that

$\limsup_{r\downarrow
0}\frac{\mathrm{P}(B_{r}({\bm{f}}))}{\mathrm{P}(B_{r}({\bm{f}}^{\star}))}\leq
1\qquad\text{for all ${\bm{f}}\in{\mathbb{B}}$.}$

The intuition for the weak mode is the same as for the finite-dimensional mode
(indeed they coincide when a Lebesgue measure exists), but the weak mode
generalizes this notion to infinite-dimensional (separable) Banach spaces.
Under certain technical assumptions on the GP prior
${\bm{\mathrm{f}}}\sim{\operatorname{\mathcal{GP}}\left({\bm{\mu}},{\bm{\Sigma}}\right)}$
and the potential $\Phi^{\bm{Y}}$ (assumptions A.1 and A.2), Lambley [2023]
shows that any solution to

$\operatorname*{arg\,min}_{{\bm{f}}\in{\mathbb{H}}_{{\bm{\Sigma}}}}\underbracket{\Phi^{\bm{Y}}({\bm{f}})+\frac{1}{2}\lVert{\bm{f}}-{\bm{\mu}}\rVert_{{\mathbb{H}}_{{\bm{\Sigma}}}}^{2}}_{\eqqcolon
R_{\text{FSP}}({\bm{f}})}$ (3.1)

is a weak mode of the posterior probability measure
$\mathrm{P}_{{\mathbb{B}}}^{{\bm{Y}}}$.

Now that we know that the weak mode of the posterior measure
$\mathrm{P}_{{\mathbb{B}}}^{{\bm{Y}}}$ is the solution to an optimization
problem in the RKHS ${\mathbb{H}}_{{\bm{\Sigma}}}$, we can relate it to an
optimization problem in weight space ${\mathbb{W}}$, for which we can apply
tools from deep learning. Informally speaking, we assume that the intersection
of the set of functions represented by the neural network
${\mathbb{F}}\coloneqq\\{{\bm{f}}(\,\cdot\,,{\bm{w}})\colon{\bm{w}}\in{\mathbb{W}}\\}\subset({\mathbb{R}}^{d^{\prime}})^{{\mathbb{X}}}$
and ${\mathbb{H}}_{{\bm{\Sigma}}}$ is non-empty, and constrain the
optimization problem in equation 3.1 such that ${\bm{f}}$ belongs to both sets

$\operatorname*{arg\,min}_{{\bm{f}}\in{\mathbb{H}}_{{\bm{\Sigma}}}\cap{\mathbb{F}}}R_{\text{FSP}}({\bm{f}})$
(3.2)

Unfortunately, the framework by Lambley [2023] cannot give probabilistic
meaning to optimization problems with hard constraints of the form
${\bm{f}}\in{\mathbb{F}}$. To address this, we adopt and elaborate on the
informal strategy from Chen et al. [2021a, Remark 2.4], leaving a full proof
for future work. Denote by
$d_{{\mathbb{B}}}({\bm{f}},{\mathbb{F}})\coloneqq\inf_{{\bm{f}}_{0}\in{\mathbb{F}}}\lVert{\bm{f}}_{0}-{\bm{f}}\rVert_{{\mathbb{B}}}$
the distance of a function ${\bm{f}}\in{\mathbb{B}}$ to the set
${\mathbb{F}}\subset{\mathbb{B}}$. Under assumption A.3,
$d_{{\mathbb{B}}}({\bm{f}},{\mathbb{F}})=0$ if and only if
${\bm{f}}\in{\mathbb{F}}$ by lemma A.2. Hence, we can relax the constraint
${\bm{f}}\in{\mathbb{F}}$ by adding
$\frac{1}{\lambda}d_{{\mathbb{B}}}({\bm{f}},{\mathbb{F}})$ to the objective in
equation 3.1. The resulting optimization problem then has a rigorous
probabilistic interpretation: propositionPropRelaxedObjectiveMAP Let
assumptions A.1, A.2 and A.3 hold. For $\lambda>0$, define
$\Phi^{{\bm{Y}},\lambda}\colon{\mathbb{B}}\to{\mathbb{R}},{\bm{f}}\mapsto\Phi^{{\bm{Y}}}({\bm{f}})+\lambda^{-1}d_{{\mathbb{B}}}({\bm{f}},{\mathbb{F}}).$
A function ${\bm{f}}^{\star}$ is a weak mode of the posterior measure
$\mathrm{P}_{{\mathbb{B}}}(\mathrm{d}{\bm{f}})\propto\exp\left(-\Phi^{{\bm{Y}},\lambda}(\mathrm{d}{\bm{f}})\right)\mathrm{P}_{{\mathbb{B}}}^{{\bm{Y}},\lambda}(\mathrm{d}{\bm{f}})$
if and only if ${\bm{f}}^{\star}\in{\mathbb{H}}_{{\bm{\Sigma}}}$ minimizes

$R_{\text{FSP}}^{\lambda}\colon{\mathbb{H}}_{{\bm{\Sigma}}}\to{\mathbb{R}},{\bm{f}}\mapsto\Phi^{{\bm{Y}},\lambda}({\bm{f}})+\frac{1}{2}\lVert{\bm{f}}-{\bm{\mu}}\rVert_{{\mathbb{H}}_{{\bm{\Sigma}}}}^{2}.$

Intuitively speaking, the term
$\frac{1}{\lambda}d_{{\mathbb{B}}}({\bm{f}},{\mathbb{F}})$ forces the
minimizers of $R_{\text{FSP}}^{\lambda}$ to converge to functions in
${\mathbb{F}}$ as $\lambda\to 0$. This is formalized in definition 1.
propositionPropRelaxedObjectiveMAPConvergence Let
$\\{\lambda_{n}\\}_{n\in{\mathbb{N}}}\subset{\mathbb{R}}_{>0}$ with
$\lambda_{n}\to 0$, and
$\\{{\bm{f}}^{\star}_{n}\\}_{n\in{\mathbb{N}}}\subset{\mathbb{H}}_{{\bm{\Sigma}}}$
such that ${\bm{f}}^{\star}_{n}$ is a minimizer of
$R_{\text{FSP}}^{\lambda_{n}}$. Then
$\\{{\bm{f}}^{\star}_{n}\\}_{n\in{\mathbb{N}}}$ has an
${\mathbb{H}}_{{\bm{\Sigma}}}$-weakly convergent subsequence with limit
${\bm{f}}^{\star}\in{\mathbb{F}}$. We conjecture that
${\bm{f}}^{\star}\in{\mathbb{F}}$ is a minimizer of equation 3.2.

Finally, under certain technical assumptions, Theorem 4.4 from Cockayne et al.
[2019] can be used to show that, as $\lambda\to 0$,
$\mathrm{P}_{{\mathbb{B}}}^{{\bm{Y}},\lambda}$ converges111w.r.t. an integral
probability metric (IPM) to $\mathrm{P}_{\mathbb{F}}^{{\bm{Y}}}$ defined by
the Radon-Nikodym derivative
$\mathrm{P}_{\mathbb{F}}^{{\bm{Y}}}(\mathrm{d}{\bm{f}})\propto\exp\left(-\Phi^{{\bm{Y}}}(\mathrm{d}{\bm{f}})\right)\mathrm{P}_{\mathbb{F}}(\mathrm{d}{\bm{f}})$
where $\mathrm{P}_{\mathbb{F}}$ is the regular conditional probability measure
$\mathrm{P}_{{\mathbb{B}}}(\cdot\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{\mathrm{f}}}\in{\mathbb{F}})$.

##### The FSP-Laplace objective as an unnormalized log-density.

After training the neural network using the objective function
$R_{\text{FSP}}({\bm{w}})\coloneqq
R_{\text{FSP}}({\bm{f}}(\,\cdot\,,{\bm{w}}))$, we intuitively want to use the
same objective function to compute an approximate posterior over the weights
using a Laplace approximation. For this to be well-defined,
$R_{\text{FSP}}({\bm{w}})$ needs to be a valid unnormalized negative log-
density, i.e., ${\bm{w}}\mapsto\exp(-R_{\text{FSP}}({\bm{w}}))$ needs to be
integrable. However, due to continuous symmetries in the weights, this often
fails to be the case.

Our method works around this issue by means of local linearization. More
precisely, after training, we linearize the network around ${\bm{w}}^{\star}$
and apply a Laplace approximation to
$R_{\text{FSP}}^{\text{lin}}({\bm{w}})\coloneqq
R_{\text{FSP}}({\bm{f}}^{\text{lin}}(\,\cdot\,,{\bm{w}}))$. The RKHS
regularizer in $R_{\text{FSP}}^{\text{lin}}$ can be rewritten as

$\displaystyle\frac{1}{2}\lVert{\bm{f}}^{\text{lin}}(\,\cdot\,,{\bm{w}})-{\bm{\mu}}\rVert_{{\mathbb{H}}_{{\bm{\Sigma}}}}^{2}$
$\displaystyle=\frac{1}{2}({\bm{w}}-{\bm{\mu}}_{{\bm{w}}^{\star}})^{\top}{\bm{\Sigma}}_{{\bm{w}}^{\star}}^{\dagger}({\bm{w}}-{\bm{\mu}}_{{\bm{w}}^{\star}})+\text{const}.,$

where
$({\bm{\Sigma}}_{{\bm{w}}^{\star}}^{\dagger})_{ij}\coloneqq\langle({\bm{J}}_{{\bm{w}}^{\star}})_{i},({\bm{J}}_{{\bm{w}}^{\star}})_{j}\rangle_{{\mathbb{H}}_{{\bm{\Sigma}}}}$,
${\bm{v}}_{i}\coloneqq\langle({\bm{J}}_{{\bm{w}}^{\star}})_{i},{\bm{f}}(\,\cdot\,,{\bm{w}}^{\star})-{\bm{\mu}}\rangle_{{\mathbb{H}}_{{\bm{\Sigma}}}}$,
${\bm{\mu}}_{{\bm{w}}^{\star}}\coloneqq{\bm{w}}^{\star}-{\bm{\Sigma}}_{{\bm{w}}^{\star}}{\bm{v}}$.
Crucially, ${\bm{\Sigma}}_{{\bm{w}}^{\star}}$ is positive-(semi)definite,
which means that $R_{\text{FSP}}^{\text{lin}}$ is normalizable so long as
${\mathbb{W}}\subset\operatorname{im}\left({\bm{\Sigma}}_{{\bm{w}}^{\star}}\right)$.
Note that this approximation also induces a (potentially degenerate) Gaussian
“prior”
${\bm{\mathrm{w}}}\sim{\operatorname{\mathcal{N}}\left({\bm{\mu}}_{{\bm{w}}^{\star}},{\bm{\Sigma}}_{{\bm{w}}^{\star}}\right)}$
over the weights.

Using model linearization, we can also establish a correspondence between the
Gaussian process prior in function space and the induced prior over the
weights. Namely, the resulting expression is the Lebesgue density of the
${\mathbb{H}}_{{\bm{\Sigma}}}$-orthogonal projection of the GP prior onto the
finite-dimensional subspace spanned by the “feature functions” learned by the
neural network $({\bm{J}}_{{\bm{w}}^{\star}})_{i}$. Hence, our model inherits
the prior structure in function space on the features induced by the Jacobian,
zeroing out the probability mass in the remaining directions.

###  3.2 Algorithmic Considerations

##### Training with the FSP-Laplace objective function.

Evaluating the FSP-Laplace objective proposed in the previous section requires
computing the RKHS norm of the neural network. Unfortunately, this does
generally not admit a closed-form expression. Hence, we use the approximation

$\lVert{\bm{f}}(\,\cdot\,,{\bm{w}})-{\bm{\mu}}\rVert_{{\mathbb{H}}_{{\bm{\Sigma}}}}^{2}\approx\underbrace{({\bm{f}}({\bm{C}},{\bm{w}})-{\bm{\mu}}({\bm{C}}))^{\top}{\bm{\Sigma}}({\bm{C}},{\bm{C}})^{-1}({\bm{f}}({\bm{C}},{\bm{w}})-{\bm{\mu}}({\bm{C}}))}_{=\lVert{\bm{h}}_{\bm{C}}(\,\cdot\,,{\bm{w}})\rVert_{{\mathbb{H}}_{{\bm{\Sigma}}}}^{2}},$

where ${\bm{C}}\in{\mathbb{X}}^{n_{\bm{C}}}$ is a set of $n_{\bm{C}}$ _context
points_ and
${\bm{h}}_{\bm{C}}({\bm{x}},{\bm{w}})\coloneqq{\bm{\Sigma}}({\bm{x}},{\bm{C}}){\bm{\Sigma}}({\bm{C}},{\bm{C}})^{-1}({\bm{f}}({\bm{C}},{\bm{w}})-{\bm{\mu}}({\bm{C}}))\in{\mathbb{H}}_{{\bm{\Sigma}}}$.
The function ${\bm{h}}_{\bm{C}}$ is the minimum-norm interpolant of
${\bm{f}}(\,\cdot\,,{\bm{w}})-{\bm{\mu}}$ at ${\bm{C}}$ in
${\mathbb{H}}_{{\bm{\Sigma}}}$. Hence, the estimator of the RKHS norm provably
underestimates, i.e.,
$\lVert{\bm{h}}_{\bm{C}}(\,\cdot\,,{\bm{w}})\rVert_{{\mathbb{H}}_{{\bm{\Sigma}}}}^{2}\leq\lVert{\bm{f}}(\,\cdot\,,{\bm{w}})-{\bm{\mu}}\rVert_{{\mathbb{H}}_{{\bm{\Sigma}}}}^{2}$.

During training, we need to compute
$\lVert{\bm{h}}_{\bm{C}}(\,\cdot\,,{\bm{w}})\rVert_{{\mathbb{H}}_{{\bm{\Sigma}}}}^{2}$
at every optimizer step. Since this involves solving a linear system in
$n_{\bm{C}}$ unknowns and, more importantly, forwarding the $n_{\bm{C}}$
context points through the neural network, we need to keep $n_{\bm{C}}$ small
for computational efficiency. We find that sampling an i.i.d. set of context
points at every training iteration from a distribution $\mathrm{P}_{\bm{C}}$
is an effective strategy for keeping $n_{\bm{C}}$ low while ensuring that the
neural network is appropriately regularized. The resulting training procedure
is outlined in algorithm 2.

##### Efficient linear Laplace approximations of the FSP-Laplace objective.

Once a minimum ${\bm{w}}^{\star}$ of $R_{\text{FSP}}$ is found, we compute a
linearized Laplace approximation at ${\bm{w}}^{\star}$. The Hessian of
$R_{\text{FSP}}$ is then

${\bm{\Lambda}}={\bm{\Sigma}}_{{\bm{w}}^{\star}}^{\dagger}-\sum_{i=1}^{n}{\bm{J}}_{{\bm{w}}^{\star}}({\bm{x}}^{(i)})^{\top}{\bm{L}}^{(i)}_{{\bm{w}}^{\star}}{\bm{J}}_{{\bm{w}}^{\star}}({\bm{x}}^{(i)}).$

Again, the RKHS inner products in the entries of
${\bm{\Sigma}}_{{\bm{w}}^{\star}}^{\dagger}$ do not admit general closed-form
expressions. Hence, we use the same strategy as before to estimate
${\bm{\Sigma}}_{{\bm{w}}^{\star}}^{\dagger}\approx{\bm{J}}_{{\bm{w}}^{\star}}({\bm{C}})^{\top}{\bm{\Sigma}}({\bm{C}},{\bm{C}})^{-1}{\bm{J}}_{{\bm{w}}^{\star}}({\bm{C}})$
at another set ${\bm{C}}$ of context points. Unlike above, for a Laplace
approximation, it is vital to use a large number of context points to capture
the prior beliefs well. Luckily, ${\bm{\Sigma}}_{{\bm{w}}^{\star}}^{\dagger}$
only needs to be computed once, at the end of training. But for large networks
and large numbers of context points, it is still infeasible to compute or even
represent ${\bm{\Sigma}}_{{\bm{w}}^{\star}}^{\dagger}$ in memory. To address
this problem, we devise an efficient routine for computing (a square root of)
the approximate posterior covariance matrix ${\bm{\Lambda}}^{\dagger}$,
outlined in algorithm 1.

Algorithm 1 Linearized Laplace approximation with Gaussian process priors

FSP-Laplace${\bm{f}}$,
${\operatorname{\mathcal{GP}}\left({\bm{\mu}},{\bm{\Sigma}}\right)}$,
${\bm{C}}$, ${\mathbb{D}}$, ${\bm{w}}^{\star}$

1:${\bm{v}}\leftarrow\nicefrac{{{\bm{J}}_{{\bm{w}}^{\star}}({\bm{C}}){\bm{1}}}}{{\lVert{\bm{J}}_{{\bm{w}}^{\star}}({\bm{C}}){\bm{1}}\rVert_{2}}}$
$\triangleright$ Initial Lanczos vector (using forward-mode autodiff)

2:${\bm{L}}\leftarrow\Call{Lanczos}{{\bm{\Sigma}}({\bm{C}},{\bm{C}}),{\bm{v}}}\in{\mathbb{R}}^{n_{\bm{C}}\times
r}$ $\triangleright$ $r\ll\min(n_{\bm{C}}d^{\prime},p)$ and
${\bm{L}}{\bm{L}}^{\top}\approx{\bm{\Sigma}}({\bm{C}},{\bm{C}})^{\dagger}$

3:${\bm{M}}\leftarrow{\bm{J}}_{{\bm{w}}^{\star}}({\bm{C}})^{\top}{\bm{L}}\in{\mathbb{R}}^{p\times
r}$ $\triangleright$ using backward-mode autodiff

4:$({\bm{U}}_{\bm{M}},{\bm{D}}_{\bm{M}},\dotsc)\leftarrow\operatorname{svd}({\bm{M}})$
$\triangleright$
${\bm{M}}{\bm{M}}^{\top}={\bm{U}}_{\bm{M}}{\bm{D}}_{\bm{M}}^{2}{\bm{U}}_{\bm{M}}^{\top}$
where ${\bm{U}}_{\bm{M}}\in{\mathbb{R}}^{p\times r}$

5:${\bm{A}}\leftarrow{\bm{D}}_{\bm{M}}^{2}-\sum_{i=1}^{n}{\bm{U}}_{\bm{M}}^{\top}{\bm{J}}_{{\bm{w}}^{\star}}({\bm{x}}^{(i)})^{\top}{\bm{L}}^{(i)}_{{\bm{w}}^{\star}}{\bm{J}}_{{\bm{w}}^{\star}}({\bm{x}}^{(i)}){\bm{U}}_{\bm{M}}$
$\triangleright$ using forward-mode autodiff

6:$({\bm{U}}_{\bm{A}},{\bm{D}}_{\bm{A}})\leftarrow\operatorname{eigh}({\bm{A}})$
$\triangleright$ ${\bm{U}}_{\bm{A}},{\bm{D}}_{\bm{A}}\in{\mathbb{R}}^{r\times
r}$

7:${\bm{S}}\leftarrow{\bm{U}}_{\bm{M}}{\bm{U}}_{\bm{A}}{\bm{D}}_{\bm{A}}^{\nicefrac{{\dagger}}{{2}}}\in{\mathbb{R}}^{p\times
r}$ $\triangleright$ ${\bm{S}}{\bm{S}}^{\top}\approx{\bm{\Lambda}}^{\dagger}$

8:Find largest $k$ such that
$\operatorname{tr}({\bm{J}}_{{\bm{w}}^{\star}}({\bm{C}}){\bm{S}}_{:,1:k}({\bm{S}}_{:,1:k})^{\top}{\bm{J}}_{{\bm{w}}^{\star}}({\bm{C}})^{\top})\leq\operatorname{tr}({\bm{\Sigma}}({\bm{C}},{\bm{C}}))$

9:return
${\operatorname{\mathcal{N}}\left({\bm{w}}^{\star},{\bm{S}}_{:,1:k}({\bm{S}}_{:,1:k})^{\top}\right)}$

Our method is _matrix-free_ , i.e. it never explicitly constructs any big
(i.e., $p\times p$, $p\times n_{\bm{C}}$, or $n_{\bm{C}}\times n_{\bm{C}}$)
matrices in memory. This allows the method to scale to large models.

We start by computing a rank $r\ll\min(n_{\bm{C}}d^{\prime},p)$ approximation
of the (pseudo-)inverted kernel Gram matrix
${\bm{\Sigma}}({\bm{C}},{\bm{C}})^{\dagger}\approx{\bm{L}}{\bm{L}}^{\top}$
using fully-reorthogonalized Lanczos iteration [Golub and Van Loan, 2013,
Section 10.1] with an embedded $LDL^{\top}$-factorization [Golub and Van Loan,
2013, Section 11.3.5]. This only needs access to matrix-vector products with
${\bm{\Sigma}}({\bm{C}},{\bm{C}})$, which can be implemented efficiently
without materializing the matrix in memory. Moreover, kernel Gramians
typically exhibit rapid spectral decay, which makes the low-rank approximation
particularly accurate. The low-rank factors ${\bm{L}}$ then yield a rank $r$
approximation of
${\bm{\Sigma}}_{{\bm{w}}^{\star}}^{\dagger}\approx{\bm{M}}{\bm{M}}^{\top}$
with ${\bm{M}}\coloneqq{\bm{J}}_{{\bm{w}}^{\star}}({\bm{C}})^{\top}{\bm{L}}$,
which is embarassingly parallel and can be computed using backward-mode
autodiff to avoid materializing the network Jacobians in memory. An
eigendecomposition of
${\bm{M}}{\bm{M}}^{\top}={\bm{U}}_{\bm{M}}{\bm{D}}_{\bm{M}}^{2}{\bm{U}}_{\bm{M}}^{\top}$
can be computed in $\mathcal{O}(pr^{2})$ time from a thin SVD of
${\bm{M}}={\bm{U}}_{\bm{M}}{\bm{D}}_{\bm{M}}{\bm{V}}_{\bm{M}}^{\top}$. The
eigenvectors ${\bm{U}}_{\bm{M}}$ define an orthogonal projector
${\bm{P}}\coloneqq{\bm{U}}_{\bm{M}}{\bm{U}}_{\bm{M}}^{\top}$ onto the range
and an orthogonal projector
${\bm{P}}_{0}\coloneqq{\bm{I}}-{\bm{U}}_{\bm{M}}{\bm{U}}_{\bm{M}}^{\top}$ onto
the nullspace of ${\bm{M}}{\bm{M}}^{\top}$. This means that we have the
decomposition
${\bm{\Lambda}}={\bm{P}}{\bm{\Lambda}}{\bm{P}}^{\top}+{\bm{P}}_{0}{\bm{\Lambda}}{\bm{P}}_{0}^{\top}={\bm{U}}_{\bm{M}}{\bm{A}}{\bm{U}}_{\bm{M}}^{\top}+{\bm{P}}_{0}{\bm{\Lambda}}{\bm{P}}_{0}^{\top}$
with

${\bm{A}}\coloneqq{\bm{U}}_{\bm{M}}^{\top}{\bm{\Lambda}}{\bm{U}}_{\bm{M}}={\bm{D}}_{\bm{M}}^{2}-\sum_{i=1}^{n}{\bm{U}}_{\bm{M}}^{\top}{\bm{J}}_{{\bm{w}}^{\star}}({\bm{x}}^{(i)})^{\top}{\bm{L}}^{(i)}_{{\bm{w}}^{\star}}{\bm{J}}_{{\bm{w}}^{\star}}({\bm{x}}^{(i)}){\bm{U}}_{\bm{M}}.$

We find that ${\bm{P}}_{0}{\bm{\Lambda}}{\bm{P}}_{0}^{\top}$ is negligible in
practice, and hence, we approximate
${\bm{\Lambda}}\approx{\bm{U}}_{\bm{M}}{\bm{A}}{\bm{U}}_{\bm{M}}^{\top}$.

Finally, by computing an eigendecomposition of
${\bm{A}}={\bm{U}}_{\bm{A}}{\bm{D}}_{\bm{A}}{\bm{U}}_{\bm{A}}^{\top}\in{\mathbb{R}}^{r\times
r}$ in $\mathcal{O}(r^{3})$ time, we obtain an eigendecomposition of
${\bm{\Lambda}}\approx{\bm{U}}_{\bm{M}}{\bm{A}}{\bm{U}}_{\bm{M}}^{\top}=({\bm{U}}_{\bm{M}}{\bm{U}}_{\bm{A}}){\bm{D}}_{\bm{A}}({\bm{U}}_{\bm{M}}{\bm{U}}_{\bm{A}})^{\top}$,
which can be used to obtain a matrix-free representation of (a square root of)
the (pseudo-)inverse of ${\bm{\Lambda}}$. Unfortunately, due to numerical
imprecision, it is difficult to distinguish between zero eigenvalues of
${\bm{A}}$ and those with small positive magnitudes. Since we need to pseudo-
invert the matrix to obtain the covariance matrix, this can explode predictive
variance. As a remedy, we use a heuristic based on the observation that, in a
linear-Gaussian model, the trace of the posterior covariance is always upper-
bounded by the trace of the prior covariance. Hence, we impose that
$\operatorname{tr}\left({\bm{J}}_{{\bm{w}}^{\star}}({\bm{C}}){\bm{\Lambda}}^{\dagger}{\bm{J}}_{{\bm{w}}^{\star}}({\bm{C}})\right)\leq\operatorname{tr}\left({\bm{\Sigma}}({\bm{C}},{\bm{C}})\right)$.
To this end, we successively zero-out the smallest eigenvalues in
${\bm{D}}_{\bm{A}}$ until the condition is fulfilled. This turns out to be an
effective strategy to combat exploding predictive variances in practice.

##  4 Experiments

We evaluate FSP-Laplace on synthetic and real-world data, demonstrating that
our method effectively incorporates beliefs specified through a GP prior;
improves performance on regression tasks for which we have prior knowledge in
the form of a kernel (Mauna Loa and ocean current modeling); and shows
competitive performance on regression, classification, out-of-distribution
detection, and Bayesian optimization tasks compared to baselines.

##### Baselines

We compare FSP-Laplace to a deterministic neural network fitted using maximum
a-posteriori estimation with an isotropic Gaussian prior on the weights (MAP)
and a neural network for which we additionally compute the standard linearized
Laplace approximation (Laplace) [Immer et al., 2021c]. We also compare our
method to a Gaussian process (GP) [Rasmussen and Williams, 2006], when the
size of the dataset allows it, and to a sparse Gaussian process (sparse GP)
[Hensman et al., 2013]. We use the full Laplace approximation when the size of
the neural network allows it, and otherwise consider the K-FAC or diagonal
approximations [Ritter et al., 2018].

##### Qualitative evaluation on synthetic data

We consider two synthetic tasks: a $1$-dimensional regression task with
randomly drawn noisy measurements of the function $y=\sin(2\pi x)$, and the
$2$-dimensional two-moons classification task from the sci-kit learn library
[Pedregosa et al., 2011]. For the regression task, data points are shown as
gray circles, functions drawn from the posterior as green lines, the inferred
mean function as a red line in figures C.1, C.2, 1, C.3 and C.4. For the
classification task, we plot the predictive mean and 2-standard-deviation of
class $1$ (blue circles) in figures C.7 and C.6. We fit the FSP-Laplace to the
data with different GP priors and find that it successfully adapts to the
specified beliefs. For instance, our method makes the functions generated by
the network periodic within the support of the context points (figure C.2),
and controls their smoothness (figure 1) and length scale (figure C.3) without
modifying the network’s architecture, or adding features. Such flexibility is
impossible with the isotropic Gaussian weight-space priors used in the Laplace
and MAP baselines. These results carry over to classification, where our
method shows a smooth decision boundary when equipped with a RBF kernel
(figure C.7) and a rough decision boundary when equipped with a Matern-1/2
kernel (figure C.6). Unlike the Laplace and MAP baselines, whose predictions
remain confident beyond the support of the data, the FSP-Laplace’s posterior
mean reverts to the zero-mean prior and its posterior variance increases.

###  4.1 Quantitative evaluation on real-world data

We now move on to investigate FSP-Laplace on two real-world scientific
modeling tasks: forecasting the concentration of $\text{CO}_{2}$ at the Mauna
Loa observatory and predicting ocean currents in the Gulf of Mexico. We then
assess the performance of our method on standard benchmark regression,
classification, out-of-distribution, and Bayesian optimization tasks. When
reporting results, we bold the value with the highest mean as well as score
whose standard-error bars overlap.

##### Mauna Loa dataset

We consider the task of modeling the monthly average atmospheric
$\text{CO}_{2}$ concentration at the Mauna Loa observatory in Hawaii from 1974
to 2024 using data collected from the NORA global monitoring
laboratory222https://gml.noaa.gov/ccgg/trends/data.html [Keeling and Whorf,
2000]. This 1-dimensional dataset is very accurately modeled by a combination
of multiple kernels proposed in Rasmussen and Williams [2006, Section 5.4.3].
We equip FSP-Laplace with this informative kernel and with additional periodic
features, and compare to a GP with the same prior and to the linearized
Laplace with same additional periodic features. Using periodic features
(partially) reflects the prior knowledge carried by the kernel. Results are
presented in table 1 and figure C.8 in Appendix. We find that incorporating
prior beliefs both via an informative prior and periodic features in FSP-
Laplace results in a strong decrease in mean squared error (MSE) over Laplace
and GP baselines. In terms of expected log-likelihood, both neural networks
under-estimate the likelihood scale, which results in poorer scores than the
GP.

Table 1: Results for the Mauna Loa $\text{CO}_{2}$ prediction and ocean
current modeling tasks. Incorporating knowledge via the GP prior in our FSP-
Laplace improves performance over the standard Laplace.

Dataset |  Expected log-likelihood ($\uparrow$) |  Mean squared error ($\downarrow$)
---|---|---
| FSP-Laplace (ours) | Laplace | GP | FSP-Laplace (ours) | Laplace | GP
Mauna loa |  -6’015.39 $\pm$ 443.13 |  -16’139.85 $\pm$ 2’433.70 | -2’123.90 $\pm$ 0.00 | 8.90 $\pm$ 1.15 |  28.26 $\pm$ 6.87 |  35.66 $\pm$ 0.00
Ocean current | 0.0428 $\pm$ 0.0258 |  -126’845.4063 $\pm$ 14’470.5270 |  -0.5069 $\pm$ 0.0000 |  0.0184 $\pm$ 0.0008 |  0.0320 $\pm$ 0.0016 | 0.0131 $\pm$ 0.0000

##### Ocean current modeling

We further evaluate FSP-Laplace’s ability to take into account prior knowledge
by considering an ocean current modeling task where we are given sparse
$2$-dimensional observations of ocean drifter buoys velocities, and we are
interested in estimating ocean currents further away from the buoys. For this,
we consider the GulfDrifters dataset [Lilly and Pérez-Brunius, 2021] and we
follow the setup by Shalashilin [2024]. We incorporate known physical
properties of ocean currents into the considered models by applying the
Helmholtz decomposition to the GP prior’s kernel [Berlinghieri et al., 2023]
as well as to the neural network. Results are shown in table 1 and in figure
2. We find that FSP-Laplace strongly improves over both Laplace and GP in
terms of expected log-likelihood, and performs competitively in terms of mean
squared error (MSE).

Figure 2: Results for the ocean current modeling experiment. We report the
mean velocity vectors, the norm of their standard-deviation and the squared
errors of compared methods. Unlike the Laplace, we find that FSP-Laplace
accurately captures ocean current dynamics.

##### Image classification

We further evaluate our method on the MNIST [LeCun and Cortes, 2010] and
FashionMNIST [Xiao et al., 2017] image classification datasets using a
convolutional neural network. We compare our model to the Laplace, MAP, and
Sparse GP baselines, as the scale of the datasets forbids exact GP inference.
Results are presented in table 2. While these datasets are particularly
challenging to our method, which is regularized in function space, we find
that FSP-Laplace performs strongly compared to baselines, obtaining the
highest expected log-likelihood and predictive accuracy. FSP-Laplace also
yields well-calibrated models with low expected calibration error (ECE).

Table 2:  Test expected log-likelihood, accuracy, expected calibration error
and OOD detection accuracy on MNIST and FashionMNIST. FSP-Laplace performs
strongly among baselines obtaining the highest expected log-likelihood,
predictive accuracy and out-of-distribution detection accuracy.

Dataset | Metric | FSP-Laplace (ours) | Laplace | MAP | Sparse GP
---|---|---|---|---|---
MNIST |  Log-likelihood ($\uparrow$) | -0.034 $\pm$ 0.001 |  -0.047 $\pm$ 0.001 |  -0.039 $\pm$ 0.001 |  -0.301 $\pm$ 0.002
Accuracy ($\uparrow$) |  -0.991 $\pm$ 0.000 |  -0.989 $\pm$ 0.000 |  -0.987 $\pm$ 0.000 |  -0.930 $\pm$ 0.001
ECE ($\downarrow$) |  -0.005 $\pm$ 0.000 |  -0.014 $\pm$ 0.001 |  -0.003 $\pm$ 0.000 |  -0.035 $\pm$ 0.001
OOD accuracy ($\uparrow$) |  -0.978 $\pm$ 0.003 |  -0.902 $\pm$ 0.009 |  -0.889 $\pm$ 0.009 |  -0.904$\pm$ 0.020
FMNIST |  Log-likelihood ($\uparrow$) | -0.261 $\pm$ 0.003 |  -0.330 $\pm$ 0.003 |  -0.281 $\pm$ 0.002 |  -0.479 $\pm$ 0.002
Accuracy ($\uparrow$) |  -0.909 $\pm$ 0.000 |  -0.896 $\pm$ 0.000 |  -0.900 $\pm$ 0.001 |  -0.848 $\pm$ 0.001
ECE ($\downarrow$) |  -0.021 $\pm$ 0.002 |  -0.005 $\pm$ 0.001 |  -0.016 $\pm$ 0.002 |  -0.027 $\pm$ 0.001
OOD accuracy ($\uparrow$) |  -0.998 $\pm$ 0.000 |  -0.844 $\pm$ 0.007 |  -0.794 $\pm$ 0.010 |  -0.938 $\pm$ 0.007

##### Out-of-distribution detection

We now investigate whether the epistemic uncertainty of FSP-Laplace is
predictive of out-of-distribution data by evaluating it on out-of-distribution
detection. We follow the setup by Osawa et al. [2019] and report the accuracy
of a single threshold to classify OOD from in-distribution (ID) data based on
the predictive uncertainty. Additional details are provided in section B.2.
FSP-Laplace performs strongly, obtaining the highest out-of-distribution
detection accuracy table 2. This can be further observed in figure C.9 in the
Appendix, where the predictive entropy of ID data points is tightly peaked
around $0$, whereas the predictive entropy of OOD data points is highly
concentrated around the maximum entropy of the softmax distribution
($\ln(10)\approx 2.3$).

##### Bayesian optimization

We finally evaluate the epistemic uncertainty of FSP-Laplace as a surrogate
model for Bayesian optimization. We consider a setup derived from Li et al.
[2024], comparing to the linearized Laplace and to a GP. Results are
summarized in figure 3. We find that our method performs particularly well on
lower-dimensional tasks, where it converges both faster and to higher rewards
than Laplace, and noticeably strongly improves over a Gaussian process on PDE.
On higher-dimensional tasks, our method performs comparatively well to Laplace
and GP baselines.

Figure 3: Results using our method (FSP-Laplace) as a surrogate model for
Bayesian optimization. We find that FSP-Laplace performs particularly well on
lower-dimensional problems, where it converges more quickly and to higher
rewards than the Laplace, obtaining comparable scores as the Gaussian process
(GP).

##  5 Related work

##### Laplace approximation in neural networks

First introduced by MacKay [1992b], the Laplace approximation gained strong
traction in the Bayesian deep learning community with the introduction of
scalable log-posterior Hessian approximations [Ritter et al., 2018, Martens,
2020], and the so-called linearized Laplace, which solves the underfitting
issue observed with standard Laplace [Immer et al., 2021c, Foong et al.,
2019]. In addition to epistemic uncertainty estimates, the Laplace
approximation also provides a method to select prior parameters via marginal
likelihood estimation [Immer et al., 2021b, Antorán et al., 2022]. Recent work
has made the linearized Laplace more scalable by restricting inference to a
subset of parameters [Daxberger et al., 2022b], by exploiting its GP
formulation [Immer et al., 2021c] to apply methods from the scalable GP
literature [Deng et al., 2022, Ortega et al., 2024, Scannell et al., 2024], or
by directly sampling from the Laplace approximation without explicitly
computing the covariance matrix [Antorán et al., 2023]. While these approaches
use the GP formulation to make the linearized Laplace more scalable, we are
unaware of any method using GP priors to incorporate interpretable prior
beliefs within the Laplace approximation.

##### BNNs with function-space priors

Function-space priors in BNNs, in the context of variational inference,
demonstrate improvements in predictive performance compared to weight-space
priors [Sun et al., 2019]. While this idea might seem sound, it turns out that
the KL divergence in the VI objective is infinite for most cases of interest
due to mismatching supports between the function-space prior and BNN’s
predictive posterior [Burt et al., 2020], which therefore requires additional
regularization to be well defined [Sun et al., 2019]. Due to this issue, other
work [Cinquin and Bamler, 2024] considers generalized VI [Knoblauch et al.,
2019] using the regularized KL divergence [Quang, 2019] or abandons
approximating the neural network’s posterior and instead uses deterministic
neural networks as basis functions for Bayesian linear regression [Ma and
Hernández-Lobato, 2021] or the mean of a sparse GP [Wild et al., 2022]. In
contrast, our method does not compute a divergence in function space, but only
the RKHS norm under the prior’s kernel, thus circumventing the issue of
mismatching support. Alternatively, rather than directly placing a prior on
the function generated by a BNN, work has investigated methods to find weight-
space priors whose pushforward approximates a target function-space measure by
minimizing a divergence [Tran et al., 2022, Flam-Shepherd, 2017], using the
Ridgelet transform [Matsubara et al., 2022], or changing the BNN’s
architecture [Pearce et al., 2019].

##### Regularizing neural networks in function space

Arguing that one ultimately only cares about the output function of the neural
network, it has been proposed to regularize neural networks in function space,
both showing that norms could be efficiently estimated and that such
regularization schemes performed well in practice [Benjamin et al., 2019, Chen
et al., 2022, Qiu et al., 2023, Rudner et al., 2023a]. Our method uses the
same RKHS norm estimator as Chen et al. [2022] (however, under a different
kernel) by sampling a new batch of context points at each update step.
Similarly, Rudner et al. [2023a] propose an empirical prior on the weights
that regularizes the neural network in function space, which they use for MAP
estimation or approximate posterior inference. The objective to fit the
empirical prior resembles ours, but unlike our method uses the kernel induced
by the last layer of the neural network and includes an additional Gaussian
prior on the weights.

##  6 Conclusion

We propose a method for performing the Laplace approximation in neural
networks using interpretable Gaussian process priors in function space. This
addresses the issue that conventional applications of approximate Bayesian
inference methods to neural networks require posing a prior in weight space,
which is virtually impossible because weight space is not interpretable. We
address the non-existence of densities in (infinite-dimensional) function
space by generalizing the notion of a MAP estimate to the limit of a sequence
of weak modes of related posterior measures. We further mitigate the
computational cost of calculating high-dimensional curvature matrices using
scalable methods from matrix-free linear algebra. By design, our method works
best in application domains where prior information can be encoded in the
language of Gaussian processes. This is confirmed by experiments on scientific
data and Bayesian optimization. In high dimensional spaces, where explicit
prior knowledge is difficult to state, Gaussian process priors are naturally
at a disadvantage. While we do demonstrate superior performance on image data,
it is yet unclear how to find good function space priors in such high-
dimensional spaces.

## Acknowledgments and Disclosure of Funding

TC, MP, PH, and RB are funded by the DFG Cluster of Excellence “Machine
Learning - New Perspectives for Science”, EXC 2064/1, project number
390727645; the German Federal Ministry of Education and Research (BMBF)
through the Tübingen AI Center (FKZ: 01IS18039A); and funds from the Ministry
of Science, Research and Arts of the State of Baden-Württemberg. MP and PH
gratefully acknowledge financial support by the European Research Council
through ERC StG Action 757275 / PANAMA. VF was supported by a Branco Weiss
Fellowship. RB further acknowledges funding by the German Research Foundation
(DFG) for project 448588364 of the Emmy Noether Programme. The authors thank
the International Max Planck Research School for Intelligent Systems (IMPRS-
IS) for supporting TC and MP.

## References

  * [1] Nina Effenberger, Marvin Pförtner, Philipp Hennig, and Nicole Ludwig.  Probabilistic wind speed downscaling for future wind power assessment.  URL https://meetingorganizer.copernicus.org/EGU24/EGU24-2034.html.
  * Indrayan [2020] Abhaya Indrayan.  Aleatory and epistemic uncertainties can completely derail medical research results.  _Journal of postgraduate medicine_ , 66, 03 2020.  doi:10.4103/jpgm.JPGM_585_19.
  * Tazi et al. [2023] Kenza Tazi, Jihao Andreas Lin, Alex S Gardner, ST John, Hong Ge, and Richard E Turner.  Towards more interpretable and robust geospatial modelling with gaussian processes.  _AGU23_ , 2023.
  * Immer et al. [2021a] Alexander Immer, Maciej Korzepa, and Matthias Bauer.  Improving predictions of bayesian neural nets via local linearization, 2021a.
  * Daxberger et al. [2022a] Erik Daxberger, Agustinus Kristiadi, Alexander Immer, Runa Eschenhagen, Matthias Bauer, and Philipp Hennig.  Laplace redux – effortless bayesian deep learning, 2022a.
  * Fortuin [2022] Vincent Fortuin.  Priors in bayesian deep learning: A review, 2022.
  * Li et al. [2024] Yucen Lily Li, Tim G. J. Rudner, and Andrew Gordon Wilson.  A study of bayesian neural network surrogates for bayesian optimization, 2024.
  * Sun et al. [2019] Shengyang Sun, Guodong Zhang, Jiaxin Shi, and Roger Grosse.  Functional variational bayesian neural networks, 2019.
  * Tran et al. [2022] Ba-Hien Tran, Simone Rossi, Dimitrios Milios, and Maurizio Filippone.  All you need is a good functional prior for bayesian deep learning, 2022\.
  * Knoblauch et al. [2019] Jeremias Knoblauch, Jack Jewson, and Theodoros Damoulas.  Generalized variational inference: Three arguments for deriving new posteriors, 2019.
  * Immer et al. [2021b] Alexander Immer, Matthias Bauer, Vincent Fortuin, Gunnar Rätsch, and Mohammad Emtiyaz Khan.  Scalable marginal likelihood estimation for model selection in deep learning, 2021b.
  * Antorán et al. [2022] Javier Antorán, David Janz, James Urquhart Allingham, Erik Daxberger, Riccardo Barbano, Eric Nalisnick, and José Miguel Hernández-Lobato.  Adapting the linearised laplace model evidence for modern deep learning, 2022.
  * Rasmussen and Williams [2006] Carl Edward Rasmussen and Christopher K. I. Williams.  _Gaussian processes for machine learning._ Adaptive computation and machine learning. MIT Press, 2006.  ISBN 026218253X.
  * Burt et al. [2020] David R. Burt, Sebastian W. Ober, Adrià Garriga-Alonso, and Mark van der Wilk.  Understanding variational inference in function-space, 2020.
  * Ma and Hernández-Lobato [2021] Chao Ma and José Miguel Hernández-Lobato.  Functional variational inference based on stochastic process generators.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, _Advances in Neural Information Processing Systems_ , volume 34, pages 21795–21807. Curran Associates, Inc., 2021.  URL https://proceedings.neurips.cc/paper_files/paper/2021/file/b613e70fd9f59310cf0a8d33de3f2800-Paper.pdf.
  * MacKay [1992a] David J. C. MacKay.  The Evidence Framework Applied to Classification Networks.  _Neural Computation_ , 4(5):720–736, 09 1992a.  ISSN 0899-7667.  doi:10.1162/neco.1992.4.5.720.  URL https://doi.org/10.1162/neco.1992.4.5.720.
  * Immer et al. [2021c] Alexander Immer, Maciej Korzepa, and Matthias Bauer.  Improving predictions of bayesian neural nets via local linearization, 2021c.
  * Oxtoby [1946] John C. Oxtoby.  Invariant measures in groups which are not locally compact.  _Transactions of the American Mathematical Society_ , 60:215–237, 1946.  doi:10.1090/S0002-9947-1946-0018188-5.
  * Pförtner et al. [2022] Marvin Pförtner, Ingo Steinwart, Philipp Hennig, and Jonathan Wenger.  Physics-informed Gaussian process regression generalizes linear PDE solvers, 2022.
  * Lambley [2023] Hefin Lambley.  Strong maximum a posteriori estimation in Banach spaces with Gaussian priors.  _Inverse Problems_ , 39(12):125010, 2023.  doi:10.1088/1361-6420/ad07a4.
  * Chen et al. [2021a] Yifan Chen, Bamdad Hosseini, Houman Owhadi, and Andrew M. Stuart.  Solving and learning nonlinear PDEs with Gaussian processes.  _Journal of Computational Physics_ , 447, 2021a.  doi:10.1016/j.jcp.2021.110668.
  * Cockayne et al. [2019] Jon Cockayne, Chris J. Oates, T. J. Sullivan, and Mark Girolami.  Bayesian probabilistic numerical methods.  _SIAM Review_ , 61(4):756–789, 2019.  doi:10.1137/17M1139357.
  * Golub and Van Loan [2013] Gene H. Golub and Charles F. Van Loan.  _Matrix Computations_.  Johns Hopkins Studies in the Mathematical Sciences. The Johns Hopkins University Press, Baltimore, fourth edition, 2013.
  * Hensman et al. [2013] James Hensman, Nicolo Fusi, and Neil D. Lawrence.  Gaussian processes for big data, 2013.
  * Ritter et al. [2018] Hippolyt Ritter, Aleksandar Botev, and David Barber.  A scalable laplace approximation for neural networks.  In _International Conference on Learning Representations_ , 2018.  URL https://openreview.net/forum?id=Skdvd2xAZ.
  * Pedregosa et al. [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Keeling and Whorf [2000] Charles D Keeling and T P Whorf.  Monthly carbon dioxide measurements on Mauna Loa, Hawaii from 1958 to 1998, 2000.  URL https://doi.org/10.1594/PANGAEA.56536.
  * Lilly and Pérez-Brunius [2021] Jonathan M. Lilly and Paula Pérez-Brunius.  GulfDrifters: A consolidated surface drifter dataset for the Gulf of Mexico, January 2021.  URL https://doi.org/10.5281/zenodo.4421585.
  * Shalashilin [2024] Ivan Shalashilin.  Gaussian processes for vector fields and ocean current modelling, March 2024.  URL https://docs.jaxgaussianprocesses.com/examples/oceanmodelling/.
  * Berlinghieri et al. [2023] Renato Berlinghieri, Brian L. Trippe, David R. Burt, Ryan Giordano, Kaushik Srinivasan, Tamay Özgökmen, Junfei Xia, and Tamara Broderick.  Gaussian processes at the helm(holtz): A more fluid model for ocean currents, 2023.
  * LeCun and Cortes [2010] Yann LeCun and Corinna Cortes.  MNIST handwritten digit database.  2010\.  URL http://yann.lecun.com/exdb/mnist/.
  * Xiao et al. [2017] Han Xiao, Kashif Rasul, and Roland Vollgraf.  Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms, 2017.
  * Osawa et al. [2019] Kazuki Osawa, Siddharth Swaroop, Anirudh Jain, Runa Eschenhagen, Richard E. Turner, Rio Yokota, and Mohammad Emtiyaz Khan.  Practical deep learning with bayesian principles, 2019.
  * MacKay [1992b] David John Cameron MacKay.  Bayesian methods for adaptive models.  1992b.  URL https://api.semanticscholar.org/CorpusID:123141880.
  * Martens [2020] James Martens.  New insights and perspectives on the natural gradient method, 2020.
  * Foong et al. [2019] Andrew Y. K. Foong, Yingzhen Li, José Miguel Hernández-Lobato, and Richard E. Turner.  ’in-between’ uncertainty in bayesian neural networks, 2019.
  * Daxberger et al. [2022b] Erik Daxberger, Eric Nalisnick, James Urquhart Allingham, Javier Antorán, and José Miguel Hernández-Lobato.  Bayesian deep learning via subnetwork inference, 2022b.
  * Deng et al. [2022] Zhijie Deng, Feng Zhou, and Jun Zhu.  Accelerated linearized laplace approximation for bayesian deep learning, 2022.
  * Ortega et al. [2024] Luis A. Ortega, Simón Rodríguez Santana, and Daniel Hernández-Lobato.  Variational linearized laplace approximation for bayesian deep learning, 2024.
  * Scannell et al. [2024] Aidan Scannell, Riccardo Mereu, Paul Edmund Chang, Ella Tamir, Joni Pajarinen, and Arno Solin.  Function-space parameterization of neural networks for sequential learning.  In _The Twelfth International Conference on Learning Representations_ , 2024.  URL https://openreview.net/forum?id=2dhxxIKhqz.
  * Antorán et al. [2023] Javier Antorán, Shreyas Padhy, Riccardo Barbano, Eric Nalisnick, David Janz, and José Miguel Hernández-Lobato.  Sampling-based inference for large linear models, with application to linearised laplace, 2023.
  * Cinquin and Bamler [2024] Tristan Cinquin and Robert Bamler.  Regularized kl-divergence for well-defined function-space variational inference in bayesian neural networks, 2024.  URL https://arxiv.org/abs/2406.04317.
  * Quang [2019] Minh Ha Quang.  Regularized divergences between covariance operators and gaussian measures on hilbert spaces, 2019.
  * Wild et al. [2022] Veit D. Wild, Robert Hu, and Dino Sejdinovic.  Generalized variational inference in function spaces: Gaussian measures meet bayesian deep learning, 2022.
  * Flam-Shepherd [2017] Daniel Flam-Shepherd.  Mapping gaussian process priors to bayesian neural networks.  2017\.  URL https://api.semanticscholar.org/CorpusID:160026812.
  * Matsubara et al. [2022] Takuo Matsubara, Chris J. Oates, and François-Xavier Briol.  The ridgelet prior: A covariance function approach to prior specification for bayesian neural networks, 2022.
  * Pearce et al. [2019] Tim Pearce, Russell Tsuchida, Mohamed Zaki, Alexandra Brintrup, and Andy Neely.  Expressive priors in bayesian neural networks: Kernel combinations and periodic functions, 2019.
  * Benjamin et al. [2019] Ari S. Benjamin, David Rolnick, and Konrad Kording.  Measuring and regularizing networks in function space, 2019.
  * Chen et al. [2022] Zonghao Chen, Xupeng Shi, Tim G. J. Rudner, Qixuan Feng, WEIZHONG ZHANG, and Tong Zhang.  A neural tangent kernel perspective on function-space regularization in neural networks.  In _OPT 2022: Optimization for Machine Learning (NeurIPS 2022 Workshop)_ , 2022.  URL https://openreview.net/forum?id=E6MGIXQlKw.
  * Qiu et al. [2023] Shikai Qiu, Tim G. J. Rudner, Sanyam Kapoor, and Andrew Gordon Wilson.  Should we learn most likely functions or parameters?, 2023.
  * Rudner et al. [2023a] Tim G. J. Rudner, Sanyam Kapoor, Shikai Qiu, and Andrew Gordon Wilson.  Function-space regularization in neural networks: A probabilistic perspective, 2023a.
  * Da Costa et al. [2023] Nathaël Da Costa, Marvin Pförtner, Lancelot Da Costa, and Philipp Hennig.  Sample path regularity of Gaussian processes from the covariance kernel, 2023.
  * Kanagawa et al. [2018] Motonobu Kanagawa, Philipp Hennig, Dino Sejdinovic, and Bharath K. Sriperumbudur.  Gaussian processes and kernel methods: A review on connections and equivalences, 2018.
  * Milios et al. [2018] Dimitrios Milios, Raffaello Camoriano, Pietro Michiardi, Lorenzo Rosasco, and Maurizio Filippone.  Dirichlet-based gaussian processes for large-scale calibrated classification.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc., 2018.  URL https://proceedings.neurips.cc/paper_files/paper/2018/file/b6617980ce90f637e68c3ebe8b9be745-Paper.pdf.
  * Clanuwat et al. [2018] Tarin Clanuwat, Mikel Bober-Irizar, Asanobu Kitamoto, Alex Lamb, Yamamoto Kazuaki, and David Ha.  Deep learning for classical japanese literature, 12 2018.
  * HALTON [1960] J.H. HALTON.  On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals.  _Numerische Mathematik_ , 2:84–90, 1960.  URL http://eudml.org/doc/131448.
  * Kingma and Ba [2017] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2017.
  * Biewald [2020] Lukas Biewald.  Experiment tracking with weights and biases, 2020.  URL https://www.wandb.com/.  Software available from wandb.com.
  * Chen et al. [2021b] Hao Chen, Lili Zheng, Raed Al Kontar, and Garvesh Raskutti.  Gaussian process inference using mini-batch stochastic gradient descent: Convergence guarantees and empirical benefits, 2021b.
  * Mckay et al. [1979] M. Mckay, Richard Beckman, and William Conover.  A comparison of three methods for selecting vales of input variables in the analysis of output from a computer code.  _Technometrics_ , 21:239–245, 05 1979.  doi:10.1080/00401706.1979.10489755.
  * Balandat et al. [2020] Maximilian Balandat, Brian Karrer, Daniel R. Jiang, Samuel Daulton, Benjamin Letham, Andrew Gordon Wilson, and Eytan Bakshy.  BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization.  In _Advances in Neural Information Processing Systems 33_ , 2020.  URL http://arxiv.org/abs/1910.06403.
  * Dua and Graff [2017] Dheeru Dua and Casey Graff.  Uci machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Malinin et al. [2021] Andrey Malinin, Liudmila Prokhorenkova, and Aleksei Ustimenko.  Uncertainty in gradient boosting via ensembles, 2021.
  * Bertin-Mahieux et al. [2011] Thierry Bertin-Mahieux, Daniel P.W. Ellis, Brian Whitman, and Paul Lamere.  The million song dataset.  In _Proceedings of the 12th International Conference on Music Information Retrieval (ISMIR 2011)_ , 2011.
  * Bradbury et al. [2018] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang.  JAX: composable transformations of Python+NumPy programs, 2018\.  URL http://github.com/google/jax.
  * Hennigan et al. [2020] Tom Hennigan, Trevor Cai, Tamara Norman, Lena Martens, and Igor Babuschkin.  Haiku: Sonnet for JAX, 2020.  URL http://github.com/deepmind/dm-haiku.
  * Botev and Martens [2022] Aleksandar Botev and James Martens.  KFAC-JAX, 2022.  URL https://github.com/google-deepmind/kfac-jax.
  * Rudner et al. [2023b] Tim G. J. Rudner, Zonghao Chen, Yee Whye Teh, and Yarin Gal.  Tractable function-space variational inference in bayesian neural networks, 2023b.
  * Sensoy et al. [2018] Murat Sensoy, Lance Kaplan, and Melih Kandemir.  Evidential deep learning to quantify classification uncertainty, 2018\.  URL https://arxiv.org/abs/1806.01768.

## Appendix

Algorithm 2 RKHS-regularized model training

FSP-Laplace-Train${\bm{f}}$, ${\bm{w}}^{(0)}$,
${\operatorname{\mathcal{GP}}\left({\bm{\mu}},{\bm{\Sigma}}\right)}$,
$\mathrm{P}_{\bm{C}}$, ${\mathbb{D}}$, $b$

1:$i\leftarrow 0$ minibatch
$\mathcal{B}=({\bm{X}}_{\mathcal{B}},{\bm{Y}}_{\mathcal{B}})\sim{\mathbb{D}}$
of size $b$

2:$R_{\text{FSP}}^{(1)}({\bm{w}}^{(i)})\leftarrow-\frac{n}{b}\sum_{j=1}^{b}\log
p({\bm{y}}^{(j)}_{\mathcal{B}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{f}}({\bm{x}}^{(j)}_{\mathcal{B}},{\bm{w}}^{(i)}))$
$\triangleright$ Negative log-likelihood

3:${\bm{C}}^{(i)}\leftarrow({\bm{c}}^{(i,j)})_{j=1}^{n_{\bm{C}}}\stackrel{{\scriptstyle\textup{\tiny
i.i.d.}}}{{\sim}}\mathrm{P}_{\bm{C}}$ $\triangleright$ Sample $n_{\bm{C}}$
context points

4:$R_{\text{FSP}}^{(2)}({\bm{w}}^{(i)})\leftarrow\frac{1}{2}({\bm{f}}({\bm{C}}^{(i)},{\bm{w}}^{(i)})-{\bm{\mu}}({\bm{C}}^{(i)}))^{\top}{\bm{\Sigma}}({\bm{C}}^{(i)},{\bm{C}}^{(i)})^{-1}({\bm{f}}({\bm{C}}^{(i)},{\bm{w}}^{(i)})-{\bm{\mu}}({\bm{C}}^{(i)}))$

5:${\bm{w}}^{(i+1)}\leftarrow\Call{OptimizerStep}{R_{\text{FSP}}^{(1)}+R_{\text{FSP}}^{(2)},{\bm{w}}^{(i)}}$

6:$i\leftarrow i+1$

7:return ${\bm{w}}^{(i)}$

##  Appendix A Additional theory

###  A.1 Assumptions and their applicability

######  Assumption A.1.

${\bm{\mathrm{f}}}\sim{\operatorname{\mathcal{GP}}\left({\bm{\mu}},{\bm{\Sigma}}\right)}$
is a $d^{\prime}$-output Gaussian process with index set ${\mathbb{X}}$ on
$(\Omega,\mathcal{A},\mathrm{P})$ such that

  1. (i)

the paths of ${\bm{\mathrm{f}}}$ lie ($\mathrm{P}$-almost surely) in a real
separable Banach space ${\mathbb{B}}$ of ${\mathbb{R}}^{d^{\prime}}$-valued
functions on ${\mathbb{X}}$ with continuous point evaluation maps
$\delta_{\bm{x}}\colon{\mathbb{B}}\to{\mathbb{R}}^{d^{\prime}}$, and

  2. (ii)

$\omega\mapsto{\bm{\mathrm{f}}}(\cdot,\omega)$ is a Gaussian random variable
with values in $({\mathbb{B}},\mathcal{B}\left({\mathbb{B}}\right))$.

We denote the law of $\omega\mapsto{\bm{\mathrm{f}}}(\cdot,\omega)$ by
$\mathrm{P}_{{\mathbb{B}}}$.

For this paper, we focus on ${\mathbb{B}}=C({\mathbb{X}})$, with
${\mathbb{X}}$ being a compact metric space. In this case, assumption A.1(i)
can be verified from regularity properties of the prior covariance function
${\bm{\Sigma}}$ [see, e.g., Da Costa et al., 2023]. Moreover, the sufficient
criteria from Pförtner et al. [2022, Section B.2] show that assumption A.1(ii)
also holds in this case.

######  Assumption A.2.

The potential $\Phi^{\bm{Y}}\colon{\mathbb{B}}\to{\mathbb{R}}$ is
(norm-)continuous and, for each $\eta>0$, there is $K(\eta)\in{\mathbb{R}}$
such that $\Phi^{\bm{Y}}({\bm{f}})\geq
K(\eta)-\eta\lVert{\bm{f}}\rVert_{{\mathbb{B}}}^{2}$ for all
${\bm{f}}\in{\mathbb{B}}$.

This holds if the negative log-likelihood functions
$\ell^{(i)}\colon{\mathbb{R}}^{d^{\prime}}\to{\mathbb{R}},{\bm{f}}^{(i)}\mapsto-\log
p({\bm{y}}^{(i)}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{f}}^{(i)})$
are continuous and, for all $\eta>0$, there is $K(\eta)\in{\mathbb{R}}$ such
that
$\ell^{(i)}({\bm{f}}^{(i)})>K({\bm{f}}^{(i)})+\eta\lVert{\bm{f}}^{(i)}\rVert^{2}$
for all ${\bm{f}}^{(i)}\in{\mathbb{R}}^{d^{\prime}}$. For instance, this is
true for a Gaussian likelihood
$\ell^{(i)}({\bm{f}}^{(i)})=\frac{1}{2\lambda_{i}^{2}}\lVert{\bm{y}}^{(i)}-{\bm{f}}^{(i)}\rVert_{2}^{2}$.

######  Assumption A.3.

(i) ${\mathbb{F}}$is a subset of ${\mathbb{H}}_{{\bm{\Sigma}}}$,  (ii)
${\mathbb{F}}$is closed in ${\mathbb{B}}\supset{\mathbb{H}}_{{\bm{\Sigma}}}$,
and  (iii) ${\mathbb{F}}\subset{\mathbb{B}}$has the Heine-Borel property,
i.e., all closed and bounded subsets of ${\mathbb{F}}$ are compact in
${\mathbb{B}}$.

Assumption A.3(i) can be verified using a plethora of results from RKHS
theory. For instance, for Sobolev kernels like the Matérn family used in the
experiments, the RKHS is norm-equivalent to a Sobolev space [see, e.g.,
Kanagawa et al., 2018]. In this case, we only need the NN to be sufficiently
often (weakly) differentiable on the interior of its compact domain
${\mathbb{X}}$. The closure property from assumption A.3(ii) is more difficult
to verify directly. However, if we assume that ${\mathbb{W}}$ is compact and
that the map
${\mathbb{W}}\to{\mathbb{B}},{\bm{w}}\mapsto{\bm{f}}(\cdot,{\bm{w}})$ is
continuous, then ${\mathbb{F}}$ is compact (and hence closed) as the image of
a compact set under a continuous function. This is a reasonable assumption,
since, in practice, the weights of a neural network are represented as machine
numbers with a maximal magnitude, meaning that ${\mathbb{W}}$ is always
contained in an $\ell_{\infty}$ ball of fixed radius. Incidentally,
compactness of ${\mathbb{F}}$ also entails the Heine-Borel property from
assumption A.3(iii). Alternatively, ${\mathbb{F}}$ also has the Heine-Borel
property if it is a topological manifold (e.g., a Banach manifold), since it
is necessarily finite-dimensional.

###  A.2 Proofs

######  Lemma A.1.

Let $({\mathbb{X}},d)$ be a metric space, ${\mathbb{A}}\subseteq{\mathbb{X}}$
nonempty, and

$d(\cdot,{\mathbb{A}})\colon{\mathbb{X}}\to{\mathbb{R}}_{\geq
0},x\mapsto\inf_{a\in{\mathbb{A}}}d(x,a).$

Then $d(\cdot,{\mathbb{A}})$ is 1-Lipschitz.

###### Proof.

For all $x_{1},x_{2}\in{\mathbb{X}}$ we have

$\displaystyle d(x_{2},{\mathbb{A}})$
$\displaystyle=\inf_{a\in{\mathbb{A}}}d(x_{2},a)$ $\displaystyle\leq
d(x_{2},x_{1})+\inf_{a\in{\mathbb{A}}}d(x_{1},a)$
$\displaystyle=d(x_{1},x_{2})+d(x_{1},{\mathbb{A}})$

by the triangle inequality and hence
$d(x_{2},{\mathbb{A}})-d(x_{1},{\mathbb{A}})\leq d(x_{1},x_{2})$. Since this
argument is symmetric in $x_{1}$ and $x_{2}$, this also shows that

$-(d(x_{2},{\mathbb{A}})-d(x_{1},{\mathbb{A}}))=d(x_{1},{\mathbb{A}})-d(x_{2},{\mathbb{A}})\leq
d(x_{2},x_{1})=d(x_{1},x_{2}).$

All in all, we obtain $\lvert
d(x_{2},{\mathbb{A}})-d(x_{1},{\mathbb{A}})\rvert\leq d(x_{1},x_{2})$. ∎

######  Lemma A.2.

Let $({\mathbb{X}},d)$ be a metric space and
${\mathbb{A}}\subseteq{\mathbb{X}}$ a closed, nonempty subset with the Heine-
Borel property. Then $\inf_{a\in{\mathbb{A}}}d(x,a)$ is attained for all
$x\in{\mathbb{X}}$.

###### Proof.

Let $x\in{\mathbb{X}}$ and $r>r_{x}\coloneqq\inf_{a\in{\mathbb{A}}}d(x,a)$.
Then ${\mathbb{A}}\cap\bar{B}_{r}(x)\neq\emptyset$ as well as $d(x,a)>r$ for
all $a\in{\mathbb{A}}\setminus\bar{B}_{r}(x)$ and thus
$\inf_{a\in{\mathbb{A}}}d(x,a)=\inf_{a\in{\mathbb{A}}\cap\bar{B}_{r}(x)}d(x,a)$.
Moreover, ${\mathbb{A}}\cap\bar{B}_{r}(x)$ is compact by the Heine-Borel
property and $d(x,\,\cdot\,)$ is continuous. Hence, the claim follows from the
Weierstrass extreme value theorem. ∎

*

###### Proof.

$d_{{\mathbb{B}}}(\,\cdot\,,{\mathbb{F}})$ is (globally) 1-Lipschitz by lemma
A.1 and bounded from below by 0. Hence,
$\Phi^{{\bm{Y}},\lambda}=\Phi^{{\bm{Y}}}+\frac{1}{\lambda}d_{{\mathbb{B}}}(\,\cdot\,,{\mathbb{F}})$
is continuous and for all $\eta>0$, there is $K(\eta)\in{\mathbb{R}}$ such
that

$\Phi^{{\bm{Y}},\lambda}({\bm{f}})\geq
K(\eta)-\eta\lVert{\bm{f}}\rVert_{{\mathbb{B}}}^{2}+\frac{1}{\lambda}d_{{\mathbb{B}}}({\bm{f}},{\mathbb{F}})\geq
K(\eta)-\eta\lVert{\bm{f}}\rVert_{{\mathbb{B}}}^{2}$

by assumption A.2. The statement then follows from Theorem 1.1 in Lambley
[2023]. ∎

*

###### Proof.

Analogous to the proof of Theorem 5.2(b) from Lambley [2023]. ∎

##  Appendix B Experimental setup

###  B.1 Qualitative experiments with synthetic data

##### Regression

We sample points from the corresponding generative model

$y_{i}=\operatorname{sin}(2\pi
x_{i})+\epsilon\quad\text{with}\quad\epsilon\sim{\operatorname{\mathcal{N}}\left(0,\sigma_{n}^{2}\right)}$
(B.1)

using $\sigma_{n}=0.1$ and draw $x_{i}\sim\mathcal{U}([-1,-0.5]\cup[0.5,1])$.
We plot data points as gray circles, functions sampled from the approximate
posterior as green lines, the empirical mean function as a red line and its
empirical 2-standard-deviation interval around the mean as a green surface.
All neural networks have the same two hidden-layer architecture with $50$
neurons per layer and hyperbolic tangent activation (Tanh) functions. For FSP-
Laplace, we use $100$ context points placed on a regular grid and run a
maxmimum of $500$ Lanczos iterations. Except when stated otherwise, we
consider a GP prior with zero-valued mean function and find the parameters of
the covariance function by maximizing the log-marginal likelihood [Rasmussen
and Williams, 2006]. For the Laplace, we use the full generalized Gauss-Newton
matrix and an isotropic Gaussian prior with scale $\sigma_{p}=1$. The MAP
estimate uses the same prior. We find the parameters of the Gaussian process
priors by maximizing the log-marginal likelihood and the parameters of the
sparse GP by maximizing the evidence lower bound [Rasmussen and Williams,
2006].

##### Classification

We sample randomly perturbed data points from the two moons data [Pedregosa et
al., 2011] with noise level $\sigma_{n}=0.1$. We plot the the data points from
class 0 as red dots and those from class 1 as blue dots. We show the mean
(upper row) and 2-standard-deviation (bottom row) of the probability that a
sample ${\bm{x}}$ belongs to class 1 under the approximate posterior, which we
estimate using $K=100$ samples. We consider a two hidden-layer neural network
with $100$ neurons per layer and hyperbolic tangent activation functions. For
FSP-Laplace, we use $100$ context points placed on a regular grid and limit
Lanczos to run for at most $500$ iterations. We consider a GP prior with zero-
valued mean function and find the parameters of the covariance function by
maximizing the log-marginal likelihood [Rasmussen and Williams, 2006] using
the reparameterization of classifications labels into regression targets from
Milios et al. [2018]. For the Laplace, we use the full generalized Gauss-
Newton matrix and an isotropic Gaussian prior with scale $\sigma_{p}=1$. The
MAP estimate uses the same prior. For the Gaussian process, we Laplace-
approximate the intractable GP posterior and find the prior parameters by
maximizing the log-marginal likelihood [Rasmussen and Williams, 2006]. The
parameters of the sparse GP are found by maximizing the ELBO [Rasmussen and
Williams, 2006].

###  B.2 Quantitative experiments with real-world data

##### Mauna Loa

We consider the Mauna Loa dataset which tracks the monthly average atmospheric
$\text{CO}_{2}$ concentration at the Mauna Loa observatory in Hawaii from 1974
to 2024 [Keeling and Whorf, 2000]. We consider the first $70\%$ of the data in
chronological order as the train set (from 1974 to 2005) and the last $30\%$
as the test set (from 2009 to 2024). We standardize the features (time) and
regression targets ($\text{CO}_{2}$ concentration). We use two hidden-layer
neural networks with hyperbolic tangent activations and $50$ units each. We
augment the input of the neural networks with an additional sine and cosine
transformation of the features i.e., we use the feature vectors
$(t_{i},\sin(2\pi t_{i}/T),\cos(2\pi t_{i}/T))$ where $t_{i}$ is the time
index and $T$ is the period used in Rasmussen and Williams [2006]. For FSP-
Laplace, we use $100$ context points placed on a regular grid and limit
Lanczos to run at most $500$ iterations. For the Laplace, we use the full
generalized Gauss-Newton matrix and find the prior scale by maximizing the
marginal likelihood [Immer et al., 2021b].

##### Ocean current modeling

We consider the GulfDrifters dataset [Lilly and Pérez-Brunius, 2021] and we
follow the setup by Shalashilin [2024]. We use as training data $20$ simulated
velocity measurements (red arrows in figure 2), and consider as testing data
$544$ average velocity measurements (blue arrows in figure 2) computed over a
regular grid on the $[-90.8,-83.8]\times[24.0,27.5]$ longitude-latitude
interval. We standardize both features and regression targets. We incorporate
physical properties of ocean currents into the models by applying the
Helmholtz decomposition to the GP prior’s covariance function [Berlinghieri et
al., 2023] as well as to the neural network $f$ using the following
parameterization

$f(\cdot,{\bm{w}})=\operatorname{grad}\Phi(\cdot,{\bm{w}}_{1})+\operatorname{rot}\Psi(\cdot,{\bm{w}}_{2})$
(B.2)

where ${\bm{w}}=\\{{\bm{w}}_{1},{\bm{w}}_{2}\\}$ and
$\Phi(\cdot,{\bm{w}}_{1})$ and $\Psi(\cdot,{\bm{w}}_{2})$ are two hidden-layer
fully-connected neural networks with hyperbolic tangent activation functions
and $100$ hidden units per layer. We use $96$ context points placed on a
regular grid for the FSP-Laplace and limit Lanczos to run at most $500$
iterations. For the linearized Laplace, we use the full generalized Gauss-
Newton (GGN) and fix the prior scale to $\sigma_{p}=1$.

##### Image classification

We consider the MNIST [LeCun and Cortes, 2010] and FashionMNIST [Xiao et al.,
2017] image classification datasets. We standardizing the images, fit the
models on a random partition of $90\%$ of the provided train splits, keeping
the remaining $10\%$ as validation data, and evaluate the models on the test
data. We report mean and standard-errors across $5$ such random partitions of
the train data with different random seeds. We compare our model to the
Laplace, MAP, and Sparse GP baselines, as the scale of the datasets forbids
exact GP inference. The expected log-likelihood and expected calibration error
are estimated by Monte Carlo integration with $10$ posterior samples.

All neural networks have the same convolutional neural network architecture
with three convolutional layers ($3\times 3$ kernels and output channels $16$,
$32$ and $64$) interleaved with a max-pooling layer, before two fully
connected layers (with output size $128$ and $10$). For FSP-Laplace, we sample
context points from the Kuzushiji-MNIST (KMNIST) dataset [Clanuwat et al.,
2018] of $28\times 28$ gray-scale images during training following Rudner et
al. [2023a] and use $25^{\prime}000$ points from the Halton low discrepancy
sequence [HALTON, 1960] to compute the covariance. We use a Categorical
likelihood, the Adam optimizer [Kingma and Ba, 2017] with a batch size of
$100$ and stop training early when the validation loss stops decreasing. We
optimize hyper-parameters using the Bayesian optimization tool provided by
Biewald [2020] and select the parameters which maximize the average validation
expected log-likelihood across $1$ random partitioning of the provided
training split into training and validation data. We find covariance function
parameters by maximizing the log-marginal likelihood from batches [Chen et
al., 2021b] using the reparameterization of classifications labels into
regression targets from Milios et al. [2018]. We optimize over kernel, prior
scale, learning-rate, $\alpha_{\epsilon}$ (introduced by Milios et al. [2018])
and activation function and select covariance functions among the RBF,
Matern-1/2, Matern-3/2, Matern-5/2 and Rational Quadratic.

##### Out-of-distribution detection with image data

We consider out-of-distribution detection with image data and a Categorical
likelihood following the setup by Osawa et al. [2019]. We aim to partition in-
distribution (ID) data from out-of-distribution (OOD) based on the mean of
entropy of the predictive distribution with respect to $10$ posterior samples
using a single threshold found by fitting a decision stump. We evaluate a
model fit on MNIST using its test data set as in-distribution data (ID) and
the test data set of FashionMNIST as out-of-distribution (OOD), and vice-versa
for a model fit on FashionMNIST. We use the same hyper-parameters as the
models used for the image classification experiments and report mean and
standard-error of our scores across the same $5$ random partitions of the
data.

##### Bayesian optimization

We consider Bayesian optimization (BO) problems derived from Li et al. [2024].
More specifically, we use the same setup but change the dimension of the
feature space of the tasks. We report mean and standard error of 5 repetitions
of the tasks across different random seeds. We use two hidden layer neural
networks with hyperbolic tangent activations and $50$ hidden units each. FSP-
Laplace uses a Matern-5/2 covariance function with constant zero mean function
whose parameters are found by maximizing the marginal likelihood [Rasmussen
and Williams, 2006]. We use $400$ context points during training and
$10^{\prime}000$ to compute the posterior covariance sampled using latin
hypercube sampling [Mckay et al., 1979]. The Gaussian process uses a zero mean
function and a Matern-5/2 covariance function following Li et al. [2024]. For
the Laplace approximation, we find the prior scale by maximizing the marginal
likelihood [Immer et al., 2021b]. Unlike FSP-Laplace which uses low-rank
factors to parameterize the posterior covariance, we found that repeatedly
computing the covariance and predictive posterior of the linearized Laplace
with the full and K-FAC generalized Gauss-Newton (GGN) matrices was often
prohibitively slow in the BO setup. We therefore use the K-FAC approximation
to the GGN where possible (Branin and PDE) and the diagonal approximation
otherwise (Ackley, Hartmann, Polynomial and BNN). We implemented this
experiment using the BO Torch library [Balandat et al., 2020].

##### Regression with UCI datasets

We consider tabular regression datasets from the UCI repository [Dua and
Graff, 2017]. Specifically, we perform leave-one-out 5-fold cross validation,
considering $10\%$ of the training folds as validation data, and we report the
mean and standard-error of the average expected log-likelihood on the test
fold. We report the mean rank of the methods across all datasets by assigning
rank 1 to the best scoring method as well as any method who’s standard error
overlaps with the highest score’s error bars, and recursively apply this
procedure to the methods not having yet been assigned a rank. We estimate the
expected log-likelihood using $10$ posterior samples. We encoding categorical
features as one-hot vectors and standardizing the features and labels. We
consider two hidden-layer neural networks with $50$ hidden units each and
hyperbolic tangent activations. All models have a homoskedastic noise model
with a learned scale parameter. Neural networks are fit using the Adam
optimizer [Kingma and Ba, 2017] and we stop training early when the validation
loss stops decreasing. Hyper-parameters are found just as for the image
classification experiment.

##### Out-of-distribution detection with regression data

We consider out-of-distribution (OOD) detection with tabular regression data
from the UCI datasets [Dua and Graff, 2017] following the setup from Malinin
et al. [2021]. We aim to separate test data (in-distribution) from a subset of
the song dataset [Bertin-Mahieux et al., 2011] (out-of-distribution) with the
same number of features based on the variance of the predictive posterior
estimated from $10$ posterior samples using a single threshold obtained via a
decision stump. We process the data just like in the regression experiments,
use the same model hyper-parameters and report mean and standard-error of the
scores across the same $5$ random partitions of the data.

####  B.2.1 Software

We use the JAX [Bradbury et al., 2018] and DM-Haiku [Hennigan et al., 2020]
Python libraries to implement neural networks. The generalized Gauss-Newton
matrices used in the Laplace approximations are computed using the KFAC-JAX
library [Botev and Martens, 2022]. We implemented the GPs and sparse GPs using
the GPyTorch library [Balandat et al., 2020]. We conducted experiments the
Bayesian optimization experiments using the BOTorch library [Balandat et al.,
2020].

####  B.2.2 Hardware

All models were fit using a single NVIDIA RTX 2080Ti GPU with 11GB of memory.

##  Appendix C Additional experimental results

###  C.1 Additional qualitative results

##### Regression

FSP-Laplace successfully adapts to the beliefs specified by different Gaussian
process priors in terms of periodicity (figure C.2), smoothness (figure 1) and
length scale (figure C.3) without modifying the neural network’s architecture
or adding features. We also find that our method effectively regularizes the
model when the data is very noisy (figure C.4). Finally, we show that our
methods needs more context points to capture the behavior of rougher priors
(figure C.5).

Figure C.1: Just like the Gaussian process (GP) and sparse GP, FSP-Laplace
captures the smoothness behavior specified by the RBF covariance function of
the Gaussian process prior.

Figure C.2: Unlike the linearized Laplace, FSP-Laplace allows to incorporate
periodicity within the support of the data using a periodic covariance
function and without additional periodic features.

Figure C.3: FSP-Laplace adapts to the length scale provided by the RBF
covariance function of the Gaussian process prior.

Figure C.4: FSP-Laplace is effectively regularized under strong label noise.

Figure C.5: FSP-Laplace requires more context points to accurately capture the
beliefs carried by rougher priors.

##### Classification

Similar to the GP and sparse GP baselines, we find that our method captures
the behavior of the prior, showing a smooth decision boundary when equipped
with a RBF covariance function (figure C.7) and a rough decision boundary when
equipped with a Matern-1/2 covariance function (figure C.6).

Figure C.6: FSP-Laplace with a Matern-1/2 covariance function against
baselines in the two-moons classification task. Similar to the Gaussian
process (GP) and sparse GP, our method shows a rough decision boundary.

Figure C.7: FSP-Laplace with a RBF covariance function against baselines in
the two-moons classification task. Similar to the Gaussian process (GP) and
sparse GP, our method shows a smooth decision boundary.

###  C.2 Additional quantitative results

##### Mauna Loa

figure C.8 shows the predictions of FSP-Laplace and the baselines on the Mauna
Loa dataset described in section 4.1 and in section C.2. Incorporating prior
beliefs both via an informative prior and periodic features in FSP-Laplace
results in an improved fit over the Laplace and GP baselines.

Figure C.8: Regression on the Mauna Loa dataset. Incorporating prior knowledge
via a kernel tailored specifically to the dataset, our method (FSP-Laplace)
results in a strong decrease in mean square error over baselines.

##### Out-of-distribution detection with image data

We present additional results for the out-of-distribution detection experiment
presented in section C.2. figure C.9 shows the distribution of the predictive
entropy of in-distribution (ID) and out-of-distribution (OOD) data under the
sparse GP, FSP-Laplace and Laplace models. For both MNIST and FashionMNIST, we
find that the predictive entropy of ID data under FSP-Laplace is tightly
peaked around $0$ nats and that the predictive entropy of OOD data strongly
concentrates around its maximum $\ln 10\approx 2.3$ nats.

Figure C.9: Distribution of in-distribution (ID) and out-of-distribution (OOD)
samples for the MNIST and Fashion MNIST image datasets. We observe that the
predictive entropy produced by FSP-Laplace nearly perfectly partition ID and
OOD data. This is reflected in OOD accuracy in table 2.

##### Rotated MNIST and FashionMNIST

We provide an additional experiment studying the behavior of FSP-Laplace under
out-of-distribution data. We consider the setup by Rudner et al. [2023b],
Sensoy et al. [2018] and track the expected log-likelihood, accuracy and
predictive entropy of models trained on MNIST and FashionMNIST for increasing
angles of rotation of the test images. We expect the predictive entropy of a
well-calibrated neural network to grow as the inputs become increasingly
dissimilar to the training data with higher angles of rotation. Similar to the
sparse GP and the linearized Laplace baselines, we find that FSP-Laplace
yields low predictive entropy for small rotation angles and that the
predictive entropy increases with the angle on both MNIST and FashionMNIST.

Figure C.10: Rotated Fashion MNIST experiment.

Figure C.11: Rotated MNIST experiment.

##### Regression with UCI datasets

We further evaluate our method on regression datasets from the UCI repository
[Dua and Graff, 2017] and compare FSP-Laplace to the Laplace, MAP, GP, and
Sparse GP baselines. We perform leave-one-out 5-fold cross-validation, keeping
$20\%$ of the remaining train folds as validation data. We report the mean and
standard-error of the expected log-likelihood with respect to samples from the
posterior across the 5-folds. Additional details can be found in section B.2.
Results are presented in table C.1. We find that our method performs well
compared to baselines, obtaining a higher mean rank than all Bayesian methods
(Laplace, GP, and Sparse GP) but a slightly lower mean rank than the MAP
baseline ($1.364$ vs. $1.182$). In particular, our method is noticeably more
accurate than the Laplace.

Table C.1: Test log-likelihood and out-of-distribution accuracy of evaluated
methods on regression datasets from the UCI repository. We find that our
method obtains a higher mean rank than Bayesian baselines in terms of expected
log-likelihood, and out-of-distribution detection accuracy of the same order
of magnitude than baselines.

Dataset |  Expected log-likelihood ($\uparrow$) |  Out-of-distribution detection accuracy ($\uparrow$)
---|---|---
| FSP-Laplace (ours) | Laplace | MAP | GP | Sparse GP | FSP-Laplace (ours) | Laplace | GP | Sparse GP
Boston | -0.528 $\pm$ 0.028 |  -0.733 $\pm$ 0.019 |  -0.906 $\pm$ 0.274 |  -1.136 $\pm$ 0.035 |  -0.762 $\pm$ 0.077 |  0.905 $\pm$ 0.005 |  0.913 $\pm$ 0.008 |  0.948 $\pm$ 0.004 | 0.969 $\pm$ 0.003
Concrete | -0.297 $\pm$ 0.093 |  -0.538 $\pm$ 0.026 | -0.199 $\pm$ 0.056 |  -0.766 $\pm$ 0.119 |  -0.609 $\pm$ 0.050 |  0.802 $\pm$ 0.024 |  0.859 $\pm$ 0.002 | 0.878 $\pm$ 0.005 |  0.854 $\pm$ 0.007
Energy |  -1.518 $\pm$ -0.043 |  1.423 $\pm$ 0.062 |  -1.754 $\pm$ 0.058 |  -1.370 $\pm$ 0.060 |  -1.518 $\pm$ 0.032 | 1.000 $\pm$ 0.000 | 1.000 $\pm$ 0.000 | 1.000 $\pm$ 0.000 | 1.000 $\pm$ 0.000
kin8nm |  -0.153 $\pm$ 0.011 |  -0.199 $\pm$ 0.009 | -0.105 $\pm$ 0.011 | - |  -0.206 $\pm$ 0.003 |  0.628 $\pm$ 0.003 |  0.609 $\pm$ 0.015 | - | 0.645 $\pm$ 0.004
Naval |  -3.415 $\pm$ 0.052 |  -3.213 $\pm$ 0.028 |  -3.460 $\pm$ 0.048 | - |  -3.474 $\pm$ 0.019 | 1.000 $\pm$ 0.000 | 1.000 $\pm$ 0.000 | - | 1.000 $\pm$ 0.000
Power |  -0.028 $\pm$ 0.012 |  -0.016 $\pm$ 0.011 |  -0.047 $\pm$ 0.013 | - |  -0.026 $\pm$ 0.009 |  0.798 $\pm$ 0.003 | 0.820 $\pm$ 0.003 | - |  0.777 $\pm$ 0.004
Protein | -1.000 $\pm$ 0.004 |  -1.043 $\pm$ 0.006 | -1.003 $\pm$ 0.003 | - |  -1.035 $\pm$ 0.002 |  0.938 $\pm$ 0.005 | 0.992 $\pm$ 0.001 | - |  0.967 $\pm$ 0.001
Wine | -1.192 $\pm$ 0.032 |  -1.403 $\pm$ 0.014 | -1.174 $\pm$ 0.022 |  -1.313 $\pm$ 0.032 |  -1.194 $\pm$ 0.020 |  0.758 $\pm$ 0.010 |  0.763 $\pm$ 0.012 | 0.812 $\pm$ 0.008 |  0.680 $\pm$ 0.014
Yacht |  -0.556 $\pm$ 0.129 |  -0.136 $\pm$ 0.762 |  -1.644 $\pm$ 0.310 |  -1.336 $\pm$ 0.060 |  -1.532 $\pm$ 0.139 | 0.865 $\pm$ 0.020 | 0.851 $\pm$ 0.029 | 0.842 $\pm$ 0.023 |  0.651 $\pm$ 0.011
Wave |  -6.853 $\pm$ 0.046 |  -6.460 $\pm$ 0.025 |  -7.146 $\pm$ 0.195 | - |  -4.909 $\pm$ 0.001 |  0.501 $\pm$ 0.000 | 0.810 $\pm$ 0.023 | - |  0.513 $\pm$ 0.001
Denmark | -0.459 $\pm$ 0.011 |  -0.810 $\pm$ 0.010 |  -0.693 $\pm$ 0.014 | - |  -0.768 $\pm$ 0.001 | 0.702 $\pm$ 0.003 | 0.714 $\pm$ 0.006 | - |  0.626 $\pm$ 0.002
Mean rank | 1.364 | 2.727 | 1.182 | 2.800 | 2.182 | 2.182 | 1.545 | 1.200 | 1.818

##### Out-of-distribution detection on tabular data

We also investigate whether the epistemic uncertainty of FSP-Laplace is
predictive of out-of-distribution data in the regression setting by evaluating
it on out-of-distribution detection following Malinin et al. [2021]. We report
the accuracy of a single threshold to classify OOD from in-distribution (ID)
data based on the predictive uncertainty. More details are presented in
section B.2. In the context of tabular data, we find that FSP-Laplace is
generally weaker compared to the baselines and produces scores similar in
magnitude to the ones obtained with the Laplace approximation (see table C.1).