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

Authors: tristan  cinquin
marvin  pförtner
vincent  fortuin
helmholtz  ai
t u munich
philipp  hennig
robert  bamler
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 08:38:42

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 trained 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 the 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 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 (as is the case in many scientific inference tasks). At the same time, it stays competitive for black-box supervised learning problems, where neural networks typically excel.

Main paper: 
 introduction: 
 Neural networks (NNs) have become the workhorse for many machine learning tasks, but they do not quantify the uncertainty arising from data scarcity-the epistemic uncertainty. NNs therefore cannot estimate their confidence in their predictions, which is needed for safety-critical applications, decision making, and scientific modeling [1] [2] [3] . 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 (e.g., MCMC) or Sparse GP (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.by a parametric distribution (e.g., variational inference or Laplace approximations). While sampling methods can be asymptotically exact for large numbers of samples, they are typically expensive for large models, whereas variational inference (VI) and the Laplace approximation are more scalable. The Laplace approximation is particularly appealing as it does not alter maximum-a-posteriori (MAP) predictions. Compared to VI, its differential nature allows the Laplace approximation to scale to large datasets and models, and it has been shown to provide well-calibrated uncertainty estimates [4, 5] .The need for function-space priors in BNNs. While the choice of prior strongly influences the uncertainty estimates obtained from the posterior [6, 7] , there are hardly any methods for elicitation of informative priors in the literature on BNNs, with the notable exceptions of Sun et al. [8] and Tran et al. [9] . 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.) [6, 10] and has known pathologies in deep NNs [9] . 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 [11, 12] . As a solution, Sun et al. [8] 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 [13] . But it turns out that the method by Sun et al. [8] is difficult to use in practice: the Kullback-Leibler (KL) divergence regularizing VI is infinite for most function-space priors of interest [14] , and estimating finite variants of it is challenging [15] .This paper addresses 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, we first derive an objective function that regularizes the neural network in function space using a GP prior, 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 our method to large models and datasets. We show 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. We propose a novel objective function for deep neural networks that allows incorporating interpretable prior knowledge through a GP prior in function space. 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 in function space. 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. 
preliminaries: laplace approximation in weight space: 
 For a given dataset D = (X, Y ) of inputs X = (x (i) ) n i=1 ∈ X n and targets Y = (y (i) ) n i=1 ∈ Y n , we consider a model of the data that is parameterized by a neural network f : X × W → R d ′ with weights w ∈ W ⊂ R p , a likelihood p(Y | f (X, w)) := n i=1 p(y (i) | f (x (i) , w)), and a prior p(w). We seek the Bayesian posterior p(w | D) ∝ p(Y | f (X, w)) p(w). As it is intractable in BNNs, approximate inference methods have been developed. Among these, the Laplace approximation [16] to the posterior is given by a Gaussian distribution q(w = w) = N w; w ⋆ , Λ -1 , whose parameters are found by MAP estimation w ⋆ ∈ arg min w∈W R WS (w) of the weights, where R WS (w) := -log p(w | D) = -log p(w) -n i=1 log p(y (i) | f (x (i) , w)) + const.(2.1)and computing the Hessian of the negative log-posteriorΛ := -H w log p(w | D)| w=w ⋆ ∈ R p×p .The linearized Laplace approximation. Computing Λ involves the Hessian of the neural network w.r.t. its weights which is generally expensive. To make the Laplace approximation scalable, it is common to linearize the network around w ⋆ before computing Λ [17] ,f (x, w) ≈ f (x, w ⋆ ) + J w ⋆ (x)(w -w ⋆ ) =: f lin (x, w) (2.2)with Jacobian J w ⋆ (x) = D w f (x, w)| w=w ⋆ . Thus, the approximate posterior precision matrix Λ isEQUATIONwhereH (i) w ⋆ = H f log p(y (i) | f ) f =f (x (i) ,w ⋆ ). Thus, under the linearized network, the Hessian of the negative log-likelihood (NLL) coincides with the generalized Gauss-Newton matrix (GGN) of the NLL. Crucially, the GGN is positive-(semi)definite even if -H w log p(Y | f (X, w))| w=w ⋆ is not. 
fsp-laplace: laplace approximation under function-space priors: 
 A conventional Laplace approximation in neural networks requires a prior in weight space, with the issues discussed in Section 1. We now 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.2 proposes a scalable algorithm for the linearized Laplace approximation with a function-space prior using matrix-free linear algebra.We motivate our Laplace approximation in function space through the lens of MAP estimation in an (infinite-dimensional) function space 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 H Σ . We then constrain the objective to the set of functions representable by the neural network and minimize it using tools from deep learning. Finally, local linearization of the neural network turns this objective into a valid negative log-density, from which we obtain the posterior covariance by computing its Hessian.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 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 ′ -output GP prior f ∼ GP (µ, Σ) with index set X. However, a (nonparametric) GP takes its values in an infinite-dimensional (Banach) space B of functions, where no such reference Lebesgue measure exists [18] . Rather, the GP induces a prior measure P B on B [19, Section B.2] . We are thus interested in the "mode" of the posterior measure P Y B under P B defined by the Radon-Nikodymderivative P Y B (df ) ∝ exp -Φ Y (df ) P B (df )where Φ Y is the potential, in essence the negative log-likelihood functional of the model. In our case, Φ Y (f ) = -n i=1 log p(y (i) | f (x (i) )). Similar to Bayes' rule in finite dimensions, this Radon-Nikodym derivative relates the prior measure to the posterior measure by reweighting. Since there is no Lebesgue measure, the (standard) mode is undefined, and we therefore follow Lambley [20] , using so-called weak modes. Definition 1 (Weak Mode [see e.g., 20, Definition 2.1]). Let B be a separable Banach space and let P be a probability measure on (B, B (B)). A weak mode of P is any point f ⋆ ∈ supp P such thatlim sup r↓0 P(B r (f )) P(B r (f ⋆ )) ≤ 1 for all f ∈ B. (3.1)Above, B (B) denotes the Borel σ-algebra on B and B r (f ) ⊂ B is an open ball with radius r centered at a point f ∈ 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 potential Φ Y (see Assumption A.2 in the appendix), Lambley [20] shows that any solution to arg minf ∈H Σ Φ Y (f ) + 1 2 ∥f -µ∥ 2 H Σ =:RFSP(f ) (3.2)is a weak mode of the posterior probability measure P Y B . Equation (3.2) casts the weak mode of the posterior measure P Y B as the solution of an optimization problem in the RKHS H Σ . We can now relate it to an optimization problem in weight space W, and apply tools from deep learning. Informally speaking, we assume that the intersection of the set of functions represented by the neural networkF := {f ( • , w) : w ∈ W} ⊂ B ⊂ (R d ′ )X and H Σ is non-empty, and constrain the optimization problem in Equation (3.2) such that f belongs to both sets arg minEQUATIONUnfortunately, the framework by Lambley [20] cannot give probabilistic meaning to optimization problems with hard constraints of the form f ∈ F. To address this, we adopt and elaborate on the informal strategy from Chen et al. B on the observation that d B (f , F) + ϵ λ = 0, where ϵ λ is independent centered Gaussian measurement noise with variance λ2 . Formally:Proposition 1. Let Assumptions A.1 to A.3 hold. For λ > 0, define Φ Y ,λ : B → R, f → Φ Y (f ) + 1 2λ 2 d 2 B (f , F).Then the posterior measure P Y ,λ B (df ) ∝ exp -Φ Y ,λ (df ) P B (df ) has at least one weak mode f ⋆ ∈ H Σ , and the weak modes of P Y ,λ B coincide with the minimizers ofR λ FSP : H Σ → R, f → Φ Y ,λ (f ) + 1 2 ∥f -µ∥ 2 H Σ . (3.4) As λ → 0, the term 1 2λ 2 d 2 B (f , F) forces the minimizers of R λ FSP to converge to functions in H Σ ∩ F that minimize R FSP : Proposition 2. Let Assumptions A.1 to A.3 hold. Let {λ n } n∈N ⊂ R >0 with λ n → 0, and {f ⋆ n } n∈N ⊂ H Σ such that f ⋆ n is a minimizer of R λn FSP . Then {f ⋆ n } n∈N has an H Σ -weakly convergent subsequence with limit f ⋆ ∈ H Σ ∩ F. Moreover, f ⋆ is a minimizer of R FSP on H Σ ∩ F.We prove Propositions 1 and 2 in Appendix A.2. To provide some intuition about the mode of convergence in Proposition 2, we point out that H Σ -weak convergence of the subsequence {f ⋆ n k } k∈N implies both (strong/norm) convergence in the path space B of the Gaussian process f (i.e., lim k→∞ ∥f ⋆ n k -f ⋆ ∥ B = 0) and pointwise convergence (i.e., lim k→∞ f ⋆ n k (x) = f ⋆ (x) ∀x ∈ X). Finally, under certain technical assumptions, Theorem 4.4 from Cockayne et al. [22] can be used to show that, as λ → 0, P Y ,λ B converges 2 to P Y F defined by the Radon-Nikodym derivative P Y F (df ) ∝ exp -Φ Y (df ) P F (df ) where P F is the regular conditional probability measureP B (• | f ∈ F).Summarizing informally, we address the nonexistence of a density-based MAP estimator by constructing a family {f ⋆ λ } λ>0 of weak modes of the related "relaxed" posteriors P Y ,λ B (df ) = P Y B (df | d B (f , F)+ϵ λ = 0) (Proposition 1). These weak modes converge to minimizers f ⋆ of the optimization problem in Equation (3.3) as λ → 0 (Proposition 2). Moreover, under certain technical assumptions, the relaxed posteriors P Y ,λ B converge to the "true" posterior P Y F as λ → 0. Given the above, we conjecture that the f ⋆ are weak modes of P Y F , and leave the proof for future work. This motivates using the objective in Equation (3.3) to find the weak mode of the posterior measure P Y F that we wish to Laplace-approximate. The next paragraph shows how this objective becomes a valid log-density.The FSP-LAPLACE objective as an unnormalized log-density. As a first step, we use Algorithm 1, discussed in detail below, to train the neural network using the objective function3 EQUATIONWe then 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 FSP (w) needs to be a valid unnormalized negative log-density, i.e., w → exp(-R FSP (w)) needs to be integrable. However, without weight-space regularization, this often fails to be the case (e.g., due to continuous symmetries in weight space). Our method works around this issue by linearizing the network locally around w ⋆ after training (see Equation (2.2)) and then applying a Laplace approximation to R lin FSP (w) := R FSP (f lin ( • , w)). In this case, the RKHS regularizer in R lin FSP can be rewritten asEQUATIONwhere(Σ † w ⋆ ) ij := ⟨(J w ⋆ ) i , (J w ⋆ ) j ⟩ H Σ (here, Σ † w ⋆ is the Moore-Penrose pseudoinverse), v i := ⟨(J w ⋆ ) i , f ( • , w ⋆ ) -µ⟩ H Σ , µ w ⋆ := w ⋆ -Σ w ⋆ v. Crucially, Σ w ⋆ is positive-(semi)definite. This means that exp(-R lin FSP ( • )) is normalizable over im (Σ w ⋆ ), i.e., we don't integrate over the null space. Note that this approximation also induces a (potentially degenerate) Gaussian "prior" w ∼ N (µ w ⋆ , Σ w ⋆ ) 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 H Σ -orthogonal projection of the GP prior onto the finite-dimensional subspace spanned by the "feature functions" (J w ⋆ ) i learned by the neural network. 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.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 not admit a closed-form expression in general. Hence, we use the approximationEQUATIONwhereC ∈ X n C is a set of n C context points and h C (x, w) := Σ(x, C)Σ(C, C) -1 (f (C, w) - µ(C)) ∈ H Σ . The function h C is the minimum-norm interpolant of f ( • , w)-µ at C in H Σ .Hence, the estimator of the RKHS norm provably underestimates, i.e., ∥hC ( • , w)∥ 2 H Σ ≤ ∥f ( • , w) -µ∥ 2 H Σ . During training, we need to compute ∥h C ( • , w)∥ 2H Σ at every optimizer step. Since this involves solving a linear system in n C unknowns and, more importantly, evaluating the neural network on the n C context points, we need to keep n C small for computational efficiency. We find that sampling an i.i.d. set of context points at every training iteration from a distribution P C is an effective strategy for keeping n C small while ensuring that the neural network is appropriately regularized. The resulting training procedure is outlined in Algorithm 1. EQUATIONEQUATIONEQUATIONFSP , w)7:return w Algorithm 2 Linearized Laplace approximation with Gaussian process priors1: function FSP-LAPLACE(f , GP (µ, Σ), C, D, w ⋆ ) 2: v ← J w ⋆ (C)1 /∥J w ⋆ (C)1∥2▷ Initial Lanczos vector (using forward-mode autodiff)3:L ← LANCZOS(Σ(C, C), v) ∈ R n C d ′ ×r ▷ r ≪ min(n C d ′ , p) and LL ⊤ ≈ Σ(C, C) † 4: M ← J w ⋆ (C) ⊤ L ∈ R p×r ▷ using backward-mode autodiff 5: (U M , D M , . . . ) ← svd(M ) ▷ M M ⊤ = U M D 2 M U ⊤ M where U M ∈ R p×r 6: A ← D 2 M - n i=1 U ⊤ M J w ⋆ (x (i) ) ⊤ H (i) w ⋆ J w ⋆ (x (i) )U M ▷ using forward-mode autodiff 7: (U A , D A ) ← eigh(A) ▷ U A , D A ∈ R r×r 8: S ← U M U A D † /2 A ∈ R p×r ▷ SS ⊤ ≈ Λ † 9:Find smallest k such that diag J w ⋆ (C i )S :,k:r (S :,k:r) ⊤ J w ⋆ (C i ) ⊤ ≤ diag Σ(C i , C i ) ∀i 10:return N w ⋆ , S :,k:r (S :,k:r ) ⊤Efficient linearized Laplace approximations of the FSP-LAPLACE objective. Once a minimum w ⋆ of R FSP is found, we compute a linearized Laplace approximation at w ⋆ . The Hessian of R FSP is thenΛ = Σ † w ⋆ - n i=1 J w ⋆ (x (i) ) ⊤ H (i) w ⋆ J w ⋆ (x (i) ). (3.8)Again, the RKHS inner products in the entries of Σ † w ⋆ do not admit general closed-form expressions. Hence, we use the same strategy as before to estimate Σ † w ⋆ ≈ J w ⋆ (C) ⊤ Σ(C, C) -1 J w ⋆ (C) at another set 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, Σ † w ⋆ 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 Σ † w ⋆ in memory. To address this problem, we devise an efficient routine for computing (a square root of) the approximate posterior covariance matrix Λ † , outlined in Algorithm 2. Our method is matrix-free, i.e. it never explicitly constructs any big (i.e., p × p, p × n C , or n C × n C ) matrices in memory. This allows the method to scale to large models.We start by computing a rank r ≪ min(n C d ′ , p) approximation of the (pseudo-)inverted kernel Gram matrix Σ(C, C) † ≈ LL ⊤ using fully-reorthogonalized Lanczos iteration [23, Section 10.1] with an embedded LDL ⊤ -factorization [23, Section 11.3.5] . This only needs access to matrix-vector products with Σ(C, C), which can be implemented efficiently without materializing the matrix in memory. Moreover, kernel Gramians typically exhibit rapid spectral decay [see e.g., 24], which makes the low-rank approximation particularly accurate. The low-rank factors L then yield a rank r approximation of Σ † w ⋆ ≈ M M ⊤ with M := J w ⋆ (C) ⊤ L, which is embarassingly parallel and can be computed using backward-mode autodiff to avoid materializing the network Jacobians in memory.⊤ = U M D 2 M U ⊤ M can be computed in O(pr 2 ) time from a thin SVD of M = U M D M V ⊤ M .The eigenvectors U M define an orthogonal projector P := U M U ⊤ M onto the range and an orthogonal projector P 0 := I -U M U ⊤ M onto the nullspace of M M ⊤ . This means that we have the decompositionEQUATIONWe find that P 0 ΛP ⊤ 0 is negligible in practice, and hence, we approximateΛ ≈ U M AU ⊤ M (see Table C.2).A D A U ⊤ A ∈ R r×r in O(r 3 ) time, we obtain an eigendecomposition of Λ ≈ U M AU ⊤ M = (U M U A )D A (U M U A ) ⊤, which can be used to obtain a matrix-free representation of (a square root of) the (pseudo-)inverse of Λ. Unfortunately, due to numerical imprecision, it is difficult to distinguish between zero eigenvalues of 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 marginal variance of the posterior is always upper-bounded by the marginal variance of the prior. Hence, we impose that diagJ w ⋆ (C i )Λ † J w ⋆ (C i ) ≤ diag Σ(C i , C i ) for all i = 1, .. . , n C by successively truncating the smallest eigenvalues in D A until the condition is fulfilled. This turns out to be an effective strategy to combat exploding predictive variance in practice.Choice of context points. Methods for regularizing neural networks in function space rely on a set of context points to evaluate the regularizer [8, 9, 25, 26] . The context points should cover the set of input locations where it is plausible that the model might be evaluated when deployed. Popular strategies include uniform sampling from a bounded subset of the input space [8, 9, 25, 26] , from the training data [8] , and from additional (possibly unlabeled) datasets [25, 26] . For MAP estimation, we choose a uniform context point distribution P C or sample from other datasets for high-dimensional inputs. We compute the posterior covariance using samples C from a low-discrepancy sequence, which effectively cover high-dimensional spaces. 
experiments: 
 We evaluate FSP-LAPLACE on synthetic and real-world data, demonstrating that our method effectively incorporates beliefs specified through a GP prior; that it improves performance on regression tasks for which we have prior knowledge in the form of a kernel (Mauna Loa and ocean current modeling); and that it 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 fit 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) [17] . We further compare our method to FVI [8] , which uses GP priors with VI, to a Gaussian process (GP) [13] when the size of the dataset allows it, and to a sparse Gaussian process (sparse GP) [27] . We use the full Laplace approximation when the size of the neural network allows it, and otherwise consider the K-FAC or diagonal approximations [28] . All neural networks share the same architecture.Qualitative evaluation on synthetic data. We consider two synthetic data tasks: a 1-dimensional regression task with randomly drawn noisy measurements of the function y = sin(2πx), and the 2-dimensional two-moons classification task from the scikit-learn library [29] . For the regression task, data points are shown as gray circles, functions drawn from the posterior as green lines and the inferred mean function as a red line (see Figures 1 and C .1 to C.4). For the classification task, we plot the predictive mean and 2-standard-deviations of the predictions for class 1 (blue circles) in Figures C. 5 and C.6. We apply FSP-LAPLACE to the data with different GP priors and find that it successfully adapts to the specified beliefs. For instance, by varying the GP prior, we can make our method generate functions that are periodic within the support of the context points ( 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. Our method also captures the properties specified by the GP prior better than FVI, especially for rougher Matern-1/2 kernels (Figures 1 and C .5). While a sufficient number of context points is necessary to capture the beliefs specified by the GP prior, we find that, even with a very small number We now move on to investigate FSP-LAPLACE on two real-world scientific modeling tasks: forecasting the concentration of 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 detection, and Bayesian optimization tasks. When reporting results, we bold the score with the highest mean as well as any scores whose standard-error bars overlap.Mauna Loa dataset. We consider the task of modeling the monthly average atmospheric CO 2 concentration at the Mauna Loa observatory in Hawaii from 1974 to 2024 using data collected from the NORA global monitoring laboratory 4 [30]. This 1-dimensional dataset is very accurately modeled by a combination of multiple kernels proposed in Rasmussen and Williams [13, Section 5.4.3] . We equip FSP-LAPLACE with this informative kernel and with additional periodic features, and compare to FVI with the same prior and additional periodic features, to a GP with the same prior, and to the linearized Laplace with the same additional periodic features. Using periodic features (partially) reflects the prior knowledge carried by the kernel. Results are presented by Table 1 and Figure C.9 in the Appendix. We find that incorporating prior beliefs both via an informative prior and periodic features in FSP-LAPLACE significantly reduces the mean squared error (MSE) compared to Laplace and GP baselines, and that our method is also more accurate than FVI, which uses variational inference. In terms of expected log-likelihood, all neural networks under-estimate the likelihood scale, which results in poorer scores than the GP.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 velocities of ocean drifter buoys, and we are interested in estimating ocean currents further away from the buoys. For this, we consider the GulfDrifters dataset [31] and we follow the setup by Shalashilin [32] . We incorporate known physical properties of ocean currents into the considered models by applying the Helmholtz decomposition to the GP prior's kernel [33] as well as to the neural network. Results are shown in Table 1 and in Figure 2 . We find that FSP-LAPLACE 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). Image classification. We further evaluate our method on the MNIST [34] and FashionMNIST [35] image classification datasets using a convolutional neural network. We compare our model to FVI, Laplace, MAP, and Sparse GP baselines, as the scale of the datasets forbids exact GP inference. For FSP-LAPLACE and FVI, we compare using context points drawn from a uniform distribution (RND) and drawn from the Kuzushiji-MNIST dataset (KMNIST) following the setup by Rudner et al. [25] .Results are presented in Table 2 . Although these datasets are particularly challenging to our method, which is regularized in function space, we find that FSP-LAPLACE performs strongly and matches or exceeds the expected log-likelihood and predictive accuracy of best-performing baselines. It also yields well-calibrated models with low expected calibration error (ECE).Out-of-distribution detection. We now investigate whether the epistemic uncertainty of FSP-LAPLACE is predictive for out-of-distribution detection (OOD). We follow the setup by Osawa et al. [36] 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 Appendix B.2. FSP-LAPLACE with context points sampled from KMNIST performs strongly, obtaining the highest out-of-distribution detection accuracy (see Table 2 , note that FSP-LAPLACE makes no assumption on whether context points are in or out of distribution). This can be further observed in Figure C .10 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) ≈ 2.3). With RND context points, our method performs comparably to the Laplace baseline.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. [7] , comparing to FVI, 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. 
related work: 
 Laplace approximation in neural networks. First introduced by MacKay [37] , the Laplace approximation gained strong traction in the Bayesian deep learning community with the introduction of scalable log-posterior Hessian approximations [28, 38] , and the so-called linearized Laplace, which solves the underfitting issue observed with standard Laplace [17, [39] [40] [41] . In addition to epistemic uncertainty estimates, the Laplace approximation also provides a method to select prior parameters via marginal likelihood estimation [11, 12] . Recent work has made the linearized Laplace more scalable by restricting inference to a subset of parameters [42] , by exploiting its GP formulation [17] to apply methods from the scalable GP literature [43] [44] [45] , or by directly sampling from the Laplace approximation without explicitly computing the covariance matrix [46] . While these approaches use the GP formulation to make the linearized Laplace more scalable, we are unaware of any method that uses GP priors to incorporate interpretable prior beliefs within the Laplace approximation.BNNs with function-space priors. In the context of variational inference, function-space priors in BNNs demonstrate improvements in predictive performance compared to weight-space priors [8] . 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 [14] , which therefore requires additional regularization to be well defined [8] . Due to this issue, other work [47] considers generalized VI [10] using the regularized KL divergence [48] or abandons approximating the neural network's posterior and instead uses deterministic neural networks as basis functions for Bayesian linear regression [15] or for the mean of a sparse GP [49] . 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, researchers have investigated methods to find weight-space priors whose pushforward approximates a target function-space measure by minimizing a divergence [9, 50] , using the Ridgelet transform [51] , or changing the BNN's architecture [52] .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 [26, [53] [54] [55] . Unlike FSP-LAPLACE, none of these methods allow to specify informative beliefs via a GP prior. Our method uses the same RKHS norm estimator as Chen et al. [54] (however, under a different kernel) by sampling a new batch of context points at each update step. Similarly, Rudner et al. [26] 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 MAP objective 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. 
conclusion: 
 We propose a method for applying the Laplace approximation to neural networks with 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, leading us to propose a simple objective function. 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.Lemma A.2. Let (X, d) be a metric space and A ⊆ X a closed, nonempty subset with the Heine-Borel property. Then inf a∈A d(x, a) is attained for all x ∈ X.Proof. Let x ∈ X and r > r x := inf a∈A d(x, a). Then A ∩ Br (x) ̸ = ∅ as well as d(x, a) > r for all a ∈ A \ Br (x) and thus inf a∈A d(x, a) = inf a∈A∩ Br(x) d(x, a). Moreover, A ∩ Br (x) is compact by the Heine-Borel property and d(x, • ) is continuous. Hence, the claim follows from the Weierstrass extreme value theorem.Proposition 1. Let Assumptions A.1 to A.3 hold. For λ > 0, define Φ Y ,λ : B → R, f → Φ Y (f ) + 1 2λ 2 d 2 B (f , F).Then the posterior measure P Y ,λ B (df ) ∝ exp -Φ Y ,λ (df ) P B (df ) has at least one weak mode f ⋆ ∈ H Σ , and the weak modes of P Y ,λ B coincide with the minimizers ofEQUATIONProof. d B ( • , F) is (globally) 1-Lipschitz by Lemma A.1 and bounded from below by 0. Hence,Φ Y ,λ = Φ Y + 1 2λ 2 d 2 B ( • , F) is continuous and for all η > 0, there is K(η) ∈ R such that Φ Y ,λ (f ) ≥ K(η) -η∥f ∥ 2 B + 1 2λ 2 d 2 B (f , F) ≥0 ≥ K(η) -η∥f ∥ 2 Bby Assumption A.2. The statement then follows from Theorem 1.1 in Lambley [20] .Proposition 2. Let Assumptions A.1 to A.3 hold. Let {λ n } n∈N ⊂ R >0 with λ n → 0, and {f ⋆ n } n∈N ⊂ H Σ such that f ⋆ n is a minimizer of R λn FSP . Then {f ⋆ n } n∈N has an H Σ -weakly convergent subsequence with limit f ⋆ ∈ H Σ ∩ F. Moreover, f ⋆ is a minimizer of R FSP on H Σ ∩ F.Our proof makes use of ideas from Dashti et al. [58] and Lambley [20] .Proof. Without loss of generality, we assume µ = 0. By Assumption A.2, there are constantsK, α > 0 such that R FSP (f ) ≥ K + α∥f ∥ 2 H Σ for all f ∈ H Σ [20, Section 4.1]. Now fix an arbitrary f ∈ H Σ ∩ F. Then R FSP (f ) = R λn FSP (f ) ≥ R λn FSP (f ⋆ n ) = R FSP (f ⋆ n ) + 1 2λ 2 n d 2 B (f ⋆ n , F) (A.1) ≥ K + α∥f ⋆ n ∥ 2 H Σ + 1 2λ 2 n d 2 B (f ⋆ n , F) (A.2) ≥ K + α∥f ⋆ n ∥ 2 H Σ , and hence ∥f ⋆ n ∥ 2 H Σ ≤ 1 α (R FSP (f ) -K) ,i.e. the sequence {f ⋆ n } n∈N ⊂ H Σ is bounded. By the Banach-Alaoglu theorem [59, Theorems V.3.1 and V.4.2(d)] and the Eberlein-Šmulian theorem [59, Theorem V.13.1], there is a weakly convergent subsequence {f ⋆ n k } k∈N with limit f ⋆ ∈ H Σ . We need to show that f ⋆ ∈ F. From Equation (A.2), it follows that0 ≤ d B (f ⋆ n k , F) ≤ λ n k 2(R FSP (f ) -K -α∥f ⋆ n k ∥ 2 H Σ ) ≤ λ n k 2(R FSP (f ) -K),where the right-hand side converges to 0 as k → ∞. Hence, lim k→∞ d B (f ⋆ n k , F) = 0. The embedding ι : H Σ → B is compact [60, Corollary 3.2.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: See Appendix B.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes]Justification: We report standard error along the mean scores across cross validation folds.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes]Justification: See Appendix Appendix B.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [No] Justification: We believe that there is no direct societal impact of the work perform that needs discussion.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA]Justification: The methods developed in the paper have no risk for misuse.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes]Justification: All external assets (code and data) used in this publication are publicly available, cited in Section 4 and Appendix C.1, and their respective licenses are, to the best of our knowledge, properly respected.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. 
 
 Appendix: A.1, d B ( • , F) : B → R is continuous, which implies that d B (ι[ • ], F) : H Σ → R is sequentially weakly continuous. Hence, d B (f ⋆ , F) = lim k→∞ d B (f ⋆ n k , F) = 0(f ⋆ n k ) ≥ R FSP (f ⋆ ). Now, by Equation (A.1), R FSP (f ) ≥ R FSP (f ⋆ ) + lim sup k→∞ 1 2λ 2 n k d 2 B (f ⋆ n k , F), for all f ∈ H Σ ∩ F. For f = f ⋆ this implies that lim sup k→∞ 1 2λ 2 n k d 2 B (f ⋆ n k , F) = 0. All in all, we arrive at R FSP (f ) ≥ R FSP (f ⋆ ) for all f ∈ H Σ ∩ F, i.e. f ⋆ is a minimizer of R FSP .Regression. We sample points from the corresponding generative modely i = sin(2πx i ) + ϵ with ϵ ∼ N 0, σ 2 n (B.1)using σ n = 0.1 and draw x i ∼ U([-1, -0.5] ∪ [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 (tanh) activation functions. For FSP-LAPLACE, we use 100 context points placed on a regular grid and run a maximum of 500 Lanczos iterations. For FVI, we sample 100 context points drawn from U([-2, 2]) at each update. Except when stated otherwise, we consider a centered GP prior and find the parameters of the covariance function by maximizing the log-marginal likelihood [13] . For the Laplace, we use the full generalized Gauss-Newton matrix and an isotropic Gaussian prior with scale σ 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 [13] .Classification. We sample randomly perturbed data points from the two moons data [29] with noise level σ 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 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 over [-3.75, 3.75] × [-3.75, 3.75] and limit Lanczos to run for at most 500 iterations. For FVI, we sample 100 context points from U([-3.75, 3.75] 2 ) at each update. We consider a centered GP prior and find the parameters of the covariance function by maximizing the log-marginal likelihood [13] using the reparameterization of classifications labels into regression targets from Milios et al. [61] . For the Laplace, we use the full generalized Gauss-Newton matrix and an isotropic Gaussian prior with scale σ 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 [13] . Sparse GP parameters are found by maximizing the ELBO [13] .Mauna Loa. We consider the Mauna Loa dataset which tracks the monthly average atmospheric CO 2 concentration at the Mauna Loa observatory in Hawaii from 1974 to 2024 [30] . 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 (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πt i /T ), cos(2πt i /T )) where t i is the time index and T is the period used in Rasmussen and Williams [13] . For FSP-LAPLACE, we use 100 context points placed on a regular grid and limit Lanczos to run at most 500 iterations. For FVI, we draw uniformly 100 context points at each update. For Laplace, we use the full generalized Gauss-Newton matrix and find the prior scale by maximizing the marginal likelihood [11] .Ocean current modeling. We consider the GulfDrifters dataset [31] and we follow the setup by Shalashilin [32] . 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] × [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 [33] as well as to the neural network f using the following parameterizationf (•, w) = grad Φ(•, w 1 ) + rot Ψ(•, w 2 ) (B.2)where w = {w 1 , w 2 } and Φ(•, w 1 ) and Ψ(•, 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 FVI, we use the same 96 context points. For the linearized Laplace, we use the full generalized Gauss-Newton (GGN) and fix the prior scale to σ p = 1.Image classification. We consider the MNIST [34] and FashionMNIST [35] 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 FVI, 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 × 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 [62] of 28 × 28 gray-scale images during training following Rudner et al. [26] and use 25 ′ 000 points from the Halton low discrepancy sequence [63] to compute the covariance. We also consider sample context points uniformly over the range [p h,w,c min , p h,w,c max ] H,W,C h=1,w=1,c=1 , where H, W and C are respectively the height, width and number of channels of the images, and p h,w,c min = v h,w,c min -0.5 × ∆ h,w,c and p h,w,c max = v h,w,c max + 0.5 × ∆ h,w,c where ∆ h,w,c = v h,w,c max -v h,w,c min is the difference between the minimal (v h,w,c min ) and maximal (v h,w,c max ) values of the data set at pixel index (h, w, c). For FVI, we use the same context point distributions as FSP-LAPLACE. For the Laplace, we use the K-FAC approximation of the generalized Gauss-Newton matrix. We use a Categorical likelihood, the Adam optimizer [64] 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 Weights and Biases [65] 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 [66] using the reparameterization of classifications labels into regression targets from Milios et al. [61] . We optimize over kernel, prior scale, learning-rate, α ϵ (introduced by Milios et al. [61] ) 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. [36] . 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 models and hyper-parameters as for image classification 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. [7] . 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 [13] . We use 400 context points during training and 10 ′ 000 to compute the posterior covariance sampled using latin hypercube sampling [67] . We use the same prior for FVI and the same number of context points during training. The Gaussian process uses a zero mean function and a Matern-5/2 covariance function following Li et al. [7] . For the Laplace approximation, we find the prior scale by maximizing the marginal likelihood [11] . 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 [68] .Regression with UCI datasets. We consider tabular regression datasets from the UCI repository [69] . 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 loglikelihood 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. FSP-LAPLACE uses context points drawn uniformly over the range [p i min , p i max ] d i=1 , where d is the dimension of the feature space, andp i min = v i min -0.5 × ∆ i and p i max = v i max + 0.5 × ∆ i where ∆ i = v i max -v i minis the difference between the minimal (v i min ) and maximal (v i max ) values of the data set at feature index i. For the Laplace, we use the full generalized Gauss-Newton matrix. FVI uses the same context points as FSP-LAPLACE. Neural networks are fit using the Adam optimizer [64] 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 [69] following the setup from Malinin et al. [70] . We aim to separate test data (in-distribution) from a subset of the song dataset [71] (outof-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.We use the JAX [72] and DM-Haiku [73] Python libraries to implement neural networks. The generalized Gauss-Newton matrices used in the Laplace approximations are computed using the KFAC-JAX library [74] . We implemented the GPs and sparse GPs using the GPyTorch library [68] . We conducted experiments the Bayesian optimization experiments using the BOTorch library [68] .All models were fit using a single NVIDIA RTX 2080Ti GPU with 11GB of memory.Regression. We show additional results for the synthetic regression task described in Section 4. We find that 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). -2 -1 0 1 2 -2 0 2 FSP-Laplace (ours) -2 -1 0 1 2 FVI -2 -1 0 1 2 Laplace -2 -1 0 1 2 MAP -2 -1 0 1 2 GP -2 -1 0 1 2 Sparse GP2 -1 0 1 2 -2 0 2 FSP-Laplace (ours) -2 -1 0 1 2 GP -2 -1 0 1 2Figure C.2: Unlike the linearized Laplace, FSP-LAPLACE allows to incorporate periodicity within the support of the data using a periodic prior covariance function and without additional periodic features.- Classification. We show additional results for the two-moons classification task described in Section 4. 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 ( Effect of context points. We provide additional details on the role of the context points in FSP-LAPLACE. The goal of the context points is to regularize the neural network on a finite set of input locations which includes any point where we would like to evaluate the model. During MAP estimation (see Algorithm 1), context points are resampled at each update step to amortize the coverage of the feature space. During the posterior covariance computation (see Algorithm 2), context points are fixed and define where we regularize the model. Context points bare similarity with inducing points in variational GPs [27] in this later step as both define where the model is regularized.2 -1 0 1 2 -2 0 2 λ = 0.15 -2 -1 0 1 2 λ = 0.35 -2 -1 0 1 2 λ = 0.5We show additional results demonstrating the behavior of our model in the low context point regime. prior. We also note that our method requires more context points to capture the beliefs of rougher priors than smooth priors (see Figure C.7).Mauna Loa. We here show the figure associated with the Mauna Loa dataset experiment in Section 4.1. Figure C.9 shows the predictions of FSP-LAPLACE and the baselines on the Mauna Loa dataset. We find that incorporating prior beliefs both via an informative prior and periodic features in FSP-LAPLACE results in an improved fit over FVI, Laplace and GP baselines. Out-of-distribution detection with image data. We present additional results for the out-ofdistribution detection experiment presented in Section 4.1. Figure C.10 shows the distribution of the predictive entropy of in-distribution (ID) and out-of-distribution (OOD) data under the sparse GP, FSP-LAPLACE, FVI 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 ≈ 2.3 nats.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. [25] , Sensoy et al. [75] and track the 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 FVI, 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 which is what we expect from a well calibrated Bayesian model.Regression with UCI datasets. We further evaluate our method on regression datasets from the UCI repository [69] and compare FSP-LAPLACE to FVI, 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 2 . C .1. We find that our method performs well compared to baselines, matching or improving over the mean rank of all Bayesian methods (FVI, Laplace, GP, and Sparse GP) but a slightly lower mean rank than the MAP baseline (1.636 vs. 1.363). 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 matches or improves over the mean rank of Bayesian baselines in terms of expected log-likelihood, and obtains out-of-distribution detection accuracies similar to the best performing BNN. 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. [70] . 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 Appendix B.2. In the context of tabular data, we find that FSP-LAPLACE performs second best among BNNs and is almost as accurate as the Laplace approximation (see Table C .1) which is first. We note that FSP-LAPLACE systematically outperforms FVI in terms of out-of-distribution detection and obtains a higher mean rank (1.818 vs. 3.091).∥P 0 ΛP ⊤ 0 ∥ F is negligible. We provide evidence that the term ∥P 0 ΛP ⊤ 0 ∥ F in Section 3.2 is negligible compared to ∥Λ∥ F in four different configurations. We consider the synthetic regression and classification setups described in Appendix B.1. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.Justification: To the best of our knowledge, we have fully conformed with the code of ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction).w.r.t. an integral probability metric (IPM)In order to simplify notation, we will assume F ⊂ H Σ for the remainder of the main paper.https://gml.noaa.gov/ccgg/trends/data.html 


Decision: The submission considers using a function space prior in place of standard weight priors for the Laplace approximation in neural networks. After several approximations, this boils down to a tractable regularisation that enforces the neural network's prior predictive distribution to match that of a Gaussian process prior at a small number of context points. The empirical results are generally promising, showing the utilities of the new method. The suite of experiments is, however, not as comprehensive as recent Bayesian deep learning / Laplace papers.  

A reviewer raised a substantial concern, resulting in an unsettled debate between the reviewer and authors in the rebuttal. The issue is the theoretical justification of Laplace and its independence from the prior. The AC noted the concern is nitpickingly technically correct. However, what the authors proposed and the approximations to make things tractable, similar to other recent `Laplace` methods, is approximate and no longer strictly Laplace and seems to give good empirical results. Other reviews are generally positive about the submission. The AC thus suggests acceptance.

One thing that the reviewers have surprisingly missed is Rudner et al's tractable functional VI scheme (NeurIPS 2022). The AC encourages the authors to compare FSP-Laplace to this.

Reviews: review: paper_summary: paper_summary: This paper proposes a method to calculate the Laplace approximation of the posterior of neural networks with a prior defined directly in the functional space (Gaussian processes). Due to the absence of Lebesgue densities in infinite-dimensional functional spaces, the notion of weak mode is used to obtain an analogy to MAP estimation in the functional space. Through model linearization, the Laplace approximation is performed at the MAP parameter of the neural network. Experiments on both synthetic and real-world data demonstrate that the proposed method effectively captures the epistemic uncertainty in function fitting and classification.
strength_weakness: strength_weakness: My main concern is the relation between this paper and Sun et al. (2019) [8]. Sun et al. also consider Bayesian inference of function posteriors with a functional prior directly defined in the function spaces (e.g., Gaussian processes). The main difference seems to be that Sun et al. perform fully Bayesian inference, while this work aims to obtain a MAP estimate followed by a Laplace approximation to obtain the posterior. Sun et al.'s work should be included as a baseline, and more discussion on the relation between the two works is needed.

Similar to the difficulty in evaluating the functional KL divergence in Sun et al., when evaluating the "prior probability" (i.e., RKHS norm) of $\\mathbf\{f\}$, a set of context points must be specified to approximate it. My second concern is that this sampling strategy could be inefficient and impractical for high-dimensional datasets. For 1D or 2D datasets, the context points can be sampled uniformly or fixed to grid points. However, for higher-dimensional cases, this approach becomes less effective. This is a critical issue for the proposed method, as in high-dimensional spaces, there is little chance that the context points will adequately cover the support of the data distributions. For example, using KMNIST as the context points for MNIST may not always be feasible in general cases.

Minor:
- L 89, Section 3.1 -> Section 3.2.
- L271, L533, L597: References are missing.
questions: questions: - What will happen if the same context points are used for both training and approximating the Hessian?

Additionally, please also see my questions raised in the Weakness section regarding the relation between this work and Sun et al., 2019, and the challenges of the context point sampling strategy for high-dimensional datasets.
limitations: limitations: Yes.
paper_summary: paper_summary: The paper proposes a new method (FSP-Laplace) to place priors in the function space of deep neural nets and develop scalable Laplace approximations on top of it. The ideas are inspired in MAP estimation theory and the fact that an objective function that is actually regularising the neural network in the function space can be developed. Additionally, to face scalability and mitigate the issues of the computational costs in calculating large curvature matrices, the authors introduce methods from matrix-free linear algebra (i.e. Lanczos). Experimental results show that the FSP-Laplace method works well in such problems were makes sense to encode information in the form of a GP prior, and in general, overcomes the standard Laplace approximation.
strength_weakness: strength_weakness: Even if I (honestly) don't find many weaknesses in the work, I think that the following points make the approach a bit weaker:

- Context points in section 3.2 are kind of an obscure part of the Algorithm to me. Despite it is indicated that the best strategy is to use iid samples and that it is an effective way to compute the FSP-Laplace objective, I think a bit of related work and connections in this point could be needed. In general, it reminds me of inducing points in GPs and kind of pseudo-coresets. I also see that some limitations or issues with this sampling and the number of context points could appear, and somehow it is not really considered in the manuscript (in my opinion).

- For instance, Laplace Redux [Daxberger 2021] states clearly the advantage of using Laplace approximation that it is feasible for both regression and classification problems, once the MAP estimate is obtained. In the case of FSP-Laplace, it is not clear to me if the work is also fully applicable to both classification and regression problems in an easy way, to which likelihood models, and if doing that makes things more difficult on the FSP-Laplace evaluation part or the Laplace approximation side. Perhaps, a bit more of clarification in this direction could be useful.


**References.** 

[Daxberger 2021] -- Laplace Redux – Effortless Bayesian Deep Learning, NeurIPS 2021.
questions: questions: **Q1** -- In L201, it is mentioned that there is an observation on the fact that the trace of the posterior is 'always' upper-bounded by the trace of the prior covariance. Due to this, a little condition is imposed as a remedy. Is this an empirical observation? it is always true? or is there a proof for this?

**After rebuttal comments:** Thanks to the authors for their detailed comments and clarifications to the few questions and points I raised in my first review. I am, in general, excited with the contributions that the paper brings to the community, particularly on Bayesian NNs and topics related to the (linearised) Laplace approximation. I do think that the paper is strong and should be accepted, so I thereby increase my score to 8.
limitations: limitations: Yes
paper_summary: paper_summary: This paper proposes a functional Laplace approximation approach which is able to incorporate function space prior instead of weight space prior used by existing Laplace approximation methods. The function space prior is more meaningful than weight space ones in terms of expressing prior knowledge about the underlying functions. The core idea of this paper is, on top of [20], to resolve the problem of functional probability density issue of naïve extension to functional MAP. The results comparing to weight space LA show the effectiveness of the new functional LA.
strength_weakness: strength_weakness: The comparative experiments are not sufficient. Various efficient LA methods should be compared. 

Ortega, L.A., Santana, S.R. and Hernández-Lobato, D., 2023. Variational linearized Laplace approximation for Bayesian deep learning. arXiv preprint arXiv:2302.12565.

McInerney, J. and Kallus, N., 2024. Hessian-Free Laplace in Bayesian Deep Learning. arXiv preprint arXiv:2403.10671.

Deng, Z., Zhou, F. and Zhu, J., 2022. Accelerated linearized Laplace approximation for Bayesian deep learning. Advances in Neural Information Processing Systems, 35, pp.2695-2708.

It would be better to include some functional BNN as well. 

There are some typos, like Line89 and Line94.
questions: questions: Can you please explain the equation in Line148 in more details?

What is the meaning of 'rigorous probabilistic interpretation' in Line 128? Does that mean the new P_B(df) is a properly defined functional posterior?  

The objective function is similar to the kernel ridge regression. Can you please explain the difference?

Can you please add comparative results with functional BNN?
limitations: limitations: n/a
paper_summary: paper_summary: This work takes a function space view of the posterior distribution, and places a GP on the likelihood with neural network predictors. This generate a posterior distribution, and the authors attempt a Laplace approximation as a which they combine with matrix-free linear algebraic methods to aid in computational tractability. Numerical examples are provided demonstrating their method performs well.
strength_weakness: strength_weakness: Laplace approximations are justified by the Bernstein-von Mises theorem, which states that (under appropriate assumptions) the posterior distribution concentrates around its mode \emph{independent of the prior distribution} in the large-data limit. 
In order to use a Laplace approximation out of the box, it's important for these assumptions to be met, and the model to be in the appropriate pre-limit. 
Singular learning theory is one approach that addresses this for singular models such as neural networks (see https://www.routledge.com/Mathematical-Theory-of-Bayesian-Statistics/Watanabe/p/book/9780367734817). 
However approximating the Hessian, or regularizing it and claiming it's because there's a prior not only does nothing to fix the problem, but perpetuates a line of ill-posed research. 
While taking a function space view opens up the possibility of priors that concentrate with the dataset size and therefore cannot be discarded in the BvM limit, the authors claim that the eigenvalues of their prior term rapidly decay and discard them, killing any possibility that they have considered this concentration. 
other forms of Laplace approximations may be more suitable (see https://arxiv.org/pdf/2110.12922, https://arxiv.org/abs/2307.07785 for related methods using expansions in data-space). 
Having said this, the work does produce decent empirical results as a quadratic approximation of the posterior \emph{inspired by} Laplace approximations, however the current reasoning used to get there is ill-considered.

The dependence of $f$ appearing in Propositions 1 and 2 should be made explicit.

The proofs given in the paper are not sound. In particular, the proof of Proposition 2 is "Analogous to the proof of Theorem 5.2(b) from Lambley [20].". However, Theorem 5.2(b) in [20] is not proved, where it is stated that the proof is similar to another reference. That reference also points to the proof of another corollary, stating it is analogous. A chain this long suggests not only that the proof needs to be written out in the current context to concretely check that the 3 analogies used are sound, but that the authors have not bothered to check the references themselves, let alone write out the proofs internally. The proof of Proposition 1 is also not clear, and needs to be written out, so that it is not burdensome for the reader to confirm the claims.

The claim that the term appearing in line 195 is negligible should be confirmed numerically.
questions: questions: Can you please address the weaknesses above?
limitations: limitations: No, limitations were not adequately addressed. However, there is little view for negative societal impact in work of this nature.

Comments: title: Response
Proof of Proposition 2
comment: Thank you for your response and clarifications. 

> [comparison to deep ensembles]

Thank you for bringing up deep ensembles.
Unfortunately, this concrete criticism comes up too late in the review process for us to run additional experiments a this point.

We first wish to highlight that our contribution is a method to specify informative prior beliefs on the function represented by the neural network in the framework of the Laplace approximation.
Deep ensembles typically use isotropic Gaussian priors on the weights (i.e. $l_2$ regularization) and we are unaware of any method to pose informative function space priors in neural network ensembles.
While interesting, we believe that comparing our method to deep ensembles is not directly relevant for our paper as they are neither related to the Laplace approximation nor to function-space priors.
Just like the other baselines, we expect our method to outperform deep ensembles in cases where a GP accurately reflects prior beliefs.

Outside of the discussion on priors, we further wish to highlight that deep ensembles are also much more expensive to fit than Laplace approximations.
Training multiple neural networks entirely from scratch is unthinkable in large neural networks and in applications where the models are often updated (e.g., Bayesian optimization).
The diversity of predictive ensembles is also known to collapse when the size of the neural network increases ([Abe et al., 2024](https://arxiv.org/abs/2302.00704)).

Finally, in contrast to deep ensembles, Laplace approximations provide a parametric (approximate) posterior probability distribution over the weights of the BNN, which is useful for more than just predictive uncertainty quantification (e.g., weight pruning, model compression, continual learning).

> I would like to state that this is a field largely driven by empirical results. However, the framing of this document is as a theoretically driven and quite technical construction, which as I have outlined is not sound. This would not be an issue, except that the authors seem unwilling to change the tone and framing of the work to address its limitations. 

Thank you for acknowledging that the field is "largely driven by empirical results".
We believe that our results provide strong empirical evidence that our method effectively incorporates beliefs specified by a GP prior and provides sensible uncertainty estimates across many applications.

We respectfully disagree that the paper is "theoretically driven". 
While we include theory to motivate and justify our method, our contribution is a practical algorithm.

Nevertheless, we are willing to take your criticism and we plan to allocate a part of the additional page allowed for the camera ready version to provide more intuition about the method, and we will include a brief summary of empirical results early on in the paper (pointing to the extended results section for details).

> The authors have also not addressed my concerns regarding the decaying eigenvalues of the prior term.
> [Comment referred to above: While taking a function space view opens up the possibility of priors that concentrate with the dataset size and therefore cannot be discarded in the BvM limit, the authors claim that the eigenvalues of their prior term rapidly decay and discard them, killing any possibility that they have considered this concentration.]

As stated in our original rebuttal, we do not consider the BvM theorem to be the justification of our method.
Hence, we do not consider the decaying eigenvalues of the covariance matrices of the finite-dimensional marginals of the prior to be concerning.

> My complaint about the proof of Proposition 2 was not that it was difficult to follow, merely that it did not appear. It is reasonably straightforward linear analysis, but the burden should not be placed upon the reader to chase through multiple resources to find something resembling a proof of a stated proposition . The authors finding it difficult to follow is not good justification for its omission. I am glad the mathematics will appear in the updated version, and Proposition 1 will be clarified.

We are frankly having difficulties understanding this paragraph. In the original review, the reviewer asked for a proof of Proposition 2. We provided this proof in the comment above and promised to include it in the appendix of the camera ready version. We take the reviewer's statement that the proof is "reasonably straightforward linear analysis" as a confirmation that they did not find an error in our proof. We would have expected that this resolves the issue.

> The numerics regarding $P_0 \Lambda P_0^\top$ are reassuring. I hope that the authors intend to include them in their updated work, however no indication has been given that this is the case.

We will naturally include the numerical results supporting that $P_0 \Lambda P_0^\top$ is negligible in the Appendix of the camera ready version of the paper.
Thanks for the authors' response!
Thanks for your response, however it is still unclear if you intend to address any of these concerns in the document.

I think that the analogy between the construction of so-called "Laplace approximations" and the failure of optimization techniques to guarantee convergence outside convex settings is incorrect. First, one can obtain results under relaxations of convexity. However, I believe that the authors were referring to guarantees of convergence to local vs global minima. This is not the same as the situation in the Laplace approximation case, where the construction is known to be incorrect, and a number of alternative constructions have been provided in relevant settings that should not be ignored. The authors do acknowledge that there are issues with the construction. Addressing these limitations up front are in the best interest of both the author, and the field at large. The alternative seems to be acknowledging the flaws in private but trying to sweep things under the rug for the purposes of publication. This is a particular issue in the LA setting, where the empirical results obtained are usually poor when compared to ensembling, which is presumably why the methods are rarely compared to network ensembles. The authors have also not addressed my concerns regarding the decaying eigenvalues of the prior term.

My complaint about the proof of Proposition 2 was not that it was difficult to follow, merely that it did not appear. It is reasonably straightforward linear analysis, but the burden should not be placed upon the reader to chase through multiple resources to find something resembling a proof of a stated proposition . The authors finding it difficult to follow is not good justification for its omission. I am glad the mathematics will appear in the updated version, and Proposition 1 will be clarified.

The numerics regarding $P_0 \Lambda P_0^\top$ are reassuring. I hope that the authors intend to include them in their updated work, however no indication has been given that this is the case.

I would like to state that this is a field largely driven by empirical results. However, the framing of this document is as a theoretically driven and quite technical construction, which as I have outlined is not sound. This would not be an issue, except that the authors seem unwilling to change the tone and framing of the work to address its limitations. Since there has been no indication that they intend to address these issues, I cannot change my score at this stage.
The markdown interpreter does not seem to work correctly in math mode. We apologise for this and can provide a PDF with the proof upon request.

**Proposition 2.**
*Let Assumptions A.1 to A.3 hold.*

*Let $\{\lambda_n\}_{n \in N} \subset R_{> 0}$ with $\lambda_n \to 0$, and $\{ f^\star_n \}_{n \in \N} \subset H_{\Sigma}$ such that $f^\star_n$ is a minimizer of $\R_\text{FSP}^{\lambda_n}$.*

Then 

$\{f_n^\star\}_{n \in \mathbb{N}}$ 

has an $H_{\Sigma}$-weakly convergent subsequence with limit 

$f^\star \in F$.

Our proof makes use of ideas from (Dashti et al., 2013) and (Lambley, 2023).

*Proof.*
Without loss of generality, we assume $\mathbf{\mu} = \mathbf{0}$.
By Assumption A.2, there are constants $K, \alpha > 0$ such that

$R_\text{FSP}(\mathbf{f}) \ge K + \alpha \lVert \mathbf{f} \rVert_{H_{\mathbf{\Sigma}}}^2$

for all $\mathbf{f} \in H_{\mathbf{\Sigma}}$ (Lambley, 2023, Section 4.1).
Now fix an arbitrary $\mathbf{f} \in H_{\mathbf{\Sigma}} \cap F$.
Then

$R_\text{FSP}(\mathbf{f}) = R_\text{FSP}^{\lambda_n}(\mathbf{f})$ 


$\ge R_\text{FSP}^{\lambda_n}(\mathbf{f}^\star_n)$


$= R_\text{FSP}(\mathbf{f}^\star_n) + \lambda_n^{-1} d_{B}(\mathbf{f}^\star_n, F)$


$\ge K + \alpha \lVert \mathbf{f}^\star_n \rVert_{H_{\mathbf{\Sigma}}}^2 + \lambda_n^{-1} d_{B}(\mathbf{f}^\star_n, F)$


$\ge K + \alpha \lVert \mathbf{f}^\star_n \rVert_{H_{\mathbf{\Sigma}}}^2,$

and hence

$\lVert \mathbf{f}^\star_n \rVert_{H_{\mathbf{\Sigma}}}^2
\le \frac{1}{\alpha} \left( R_\text{FSP}(\mathbf{f}) - K \right),$

i.e. the sequence $\{f^\star_n\}_{n \in \mathbb{N}} \subset H_{\Sigma}$ is bounded.

By the Banach-Alaoglu theorem (Conway, 1997, Theorems V.3.1 and V.4.2(d)) and the Eberlein-Šmulian theorem (Conway, 1997, Theorem V.13.1), there is a weakly convergent subsequence $\{ \mathbf{f}^\star_{n_k} \}_{k \in \N}$ with limit $\mathbf{f}^\star \in H_{\mathbf{\Sigma}}$.

It remains to show that $\mathbf{f}^\star \in F$.
From the inequality above, it follows that

$0 \le d_{B}(\mathbf{f}^\star_{n_k}, F) \le \lambda_{n_k} \left( R_\text{FSP}(\mathbf{f}) - K - \alpha \lVert \mathbf{f}^\star_n \rVert_{H_{\mathbf{\Sigma}}}^2 \right),$

where the right-hand side converges to 0 as $k \to \infty$.

Hence, $\lim_{k \to \infty} d_{B}(\mathbf{f}^\star_{n_k}, F) = 0$.

The embedding $\iota \colon H_{\mathbf{\Sigma}} \to B$ is compact (Bogachev, 1998, Corollary 3.2.4) and hence $\{\iota[\mathbf{f}^\star_{n_k}] \}_{k \in \mathbb{N}} \subset B$ converges ($B$-strongly) to $\iota[f^\star] \in B$ (Conway, 1997, Proposition 3.3(a)).

The continuity of $d_{B}(\,\cdot\,, F)$ implies

$d_{B}(\mathbf{f}^\star, F)= d_{B}(\iota[\mathbf{f}^\star], F)$


$= d_{B} \left( \lim_{k \to \infty} \iota[\mathbf{f}^\star_{n_k}], F \right)$


$= \lim_{k \to \infty} d_{B}(\iota[\mathbf{f}^\star_{n_k}], F)$


$= \lim_{k \to \infty} d_{B}(\mathbf{f}^\star_{n_k}, F)$


$= 0,$

and by Assumption A.3 and Lemma A.2 it follows that $\mathbf{f}^\star \in F$.