Squared families are useful conjugate priors
Authors: Russell Tsuchida
Jiawei Liu
Cheng Soon Ong
Dino Sejdinovic
Venue: NeurIPS 2025 poster ()
Keywords: Bayesian inference
conjugate families
Gaussian processes
squared families
Primary area: probabilistic_methods
TL;DR: Squared families of probability distributions are conjugate prior families for many likelihoods with applications including regression.

Abstract: Squared families of probability distributions have been studied and applied in numerous machine learning contexts. Typically, they appear as likelihoods, where their advantageous computational, geometric and statistical properties are exploited for fast estimation algorithms, representational properties and statistical guarantees. Here, we investigate the use of squared families as prior beliefs in Bayesian inference. We find that they can form helpful conjugate families, often allowing for closed-form and tractable Bayesian inference and marginal likelihoods. We apply such conjugate families to Bayesian regression in feature space using end-to-end learnable neural network features. Such a setting allows for a rich multi-modal alternative to Gaussian processes with neural network features, often called deep kernel learning. We demonstrate our method on few shot learning, outperforming existing neural methods based on Gaussian processes and normalising flows.

Squared families are useful conjugate priors
Russell Tsuchida
Jiawei Liu School of Computing
Data Science and Al Group
Monash University
Australian National University 1iu3@anu, edu. au 1
russell. tsuchidadmooaashhh....
Cheng Soon Ong Data61, CSIRO
Dino Sejdinovic
RAIR, AIML
and Australian National University
The University Adelaide
chengsoon. baanggoon.ong@anu... edu.au
dano.sejdhnooriccddelladdd edu. au
Abstract
Squared families of probability distributions have been studied and applied in numerous machine learning contexts. Typically, they appear as likelihoods, where their udvantageous computational, geometric and statistical properties are exploited for fast estimation algorithms, representational properties and statistical guarantees. Here, we investigate the use of squared families prior beliefs in Bayesian inference. We find that they can form helpful conjugate families, often allowing for closed-form and tractable Bayesian inference and marginal likelihoods. We apply such conjugate families to Bayesian regression in feature space using end-to-end learnable neural network features. Such setting allows for rich multi-modal alternative to Gaussian processes with neural network features, often called deep kernel learning. We demonstrate our method on few shot learning, outperforming existing neural methods based on Gaussian processes and normalising flows
Introduction
Squared f'amilies are useful likelihoods Families {g(- feH of probability densities of the form q( {f) x [l where f E H and H is some suitably nice space of vector-valued functions, have recently emerged promising likelihoods in several different communities in machine learning. For example, functions of the form f(.) = Oxp(-) for some θ E mRxx and <6(.) =
(b(:,a)).. .. k.,,, aN)) where (a are data points and is the kernel of reproducing kernel Hilbert space (RKHS). are obtained as minimisers of regularised empirical risk problems [Marteau Ferey al. 2020). When appropriately normalised and applied to divergence minimisation, they can be used to model probability densities [Rudi and Ciliberto| 202T). When m 1, such density models can also be applied to modelling intensities (of Poisson point processes) FFaxman et al., [017]. and Bayesian variants which use Gaussian processes (GPs) instead of RKHS models are also readily available [Walder and Bishop, 2017). Separately to kernel methods, probabilistic circuits [Choi et 2020) can be chosen for f. which can be squared [Sladek] 2023 Loconte et al., 2023a Loconte and Vergari, 2025), or squared and summed [Loconte et al. 2024]. to obtain tractable and expressive models for probability density functions. Finally, neural networks can be used for f. and when they are two-layered and fully-comme conneceddde they often admin tractable and closed-for fooeddfommm normalising constants, marginal distributions and conditional distributions Tsuchida al. 2023
'odee available httpt:://githbb.ccom///eee/////t htttp:///gtttaa.cc////rrrlllslss
git
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

2.0
Mode
Modo
1,5
Mode
%i
A/
Hyperparams
1.0
Piior
Posterior Nm2
10
8.5
Fosterior N-4
A ~ Gamma GSF
$ 1,2,3
Puis(A)
0.0
6
10
X
Figure 1: (Left) Empirical Bayesian estimation of the rate parameter À of a Poisson likelihood in hierarchical model. Here the mode index = 2, is observed, but not the mode rate Aj. A classical empirical Bayesian approach might be to place shared conjugate prior over X (in this case, the Gamma density, see Appendix A), and maximise the marginal likelihood with respect the prior hyperparameters. Classical conjugate family priors are only capable of expressing a unimodal prior belief, unless they are mixture models. (Right) At inference time, we are presented N datapoints from a single mode- - Mode 1, GSF priors (solid curves) are multimodal, and are able to adapt their trimodal prior to the observations quickly, whereas the unimodal Gamma prior (dashed curves) have to slowly adapt their density to pick out the correct mode. GSFs significantly generalise mixture models of simpler base conjugate priors.
Are squared families useful priors? In this work, we examine the use of such families as priors, and more specifically, families of conjugate priors. To the best of our knowledge, we the first to consider such families as priors. Note that this question is distinct from a Bayesian treatment of squared family likelihoods [Walder and Bishop 2017 for example], where squared family is used for likelihood and some prior belief over I (for example, f ~ GP) is employed. Rather, this question is about using some other (not necessarily squared family likelihood) and a squared family prior. We first construct hierarchy of families, called generalised squared families (GSF). We then show when a prior belongs the GSF, the marginal likelihoods, posterior parameter and posterior predictive distributions all admit closed-form expressions. Furthermore, the posterior distributions remain within the GSF. This allows for rich and multimodal expression of prior beliefs within conjugate family (see Figure | for example). We apply our conjugate families to the problem of Bayesian regression feature space, which offers a model called GSF processes (GSFP) with rich multi-modal uncertainty and capabilities well-suited for few-task learning.
2
Background
Conjugate priors Let [1) be family of probability measures each supported on subset of ft. We call IT] a conjugate class for a likelihood p(u u w) if for all π E IT the posterior probability measure n' defined by <'dw | u) p(u ( Jn p(u u | w w')rrddw") is also element of n. An immediate but unhelpful conjugate class is the class @ of all probability measures, since for any likelihood prior in Q results in posterior in Q [Robert] 2007 page 98]. Often when authors speak of conjugate families, they implicitly and loosely also require that such conjugate families have closed-form or tractable normalising constants or parameters. In order facilitate such closed-form updates, conjugate families are therefore required to be small and analytically tractable. Exponential families of priors matched to exponential family likelihoods are arguably the most examples of conjugate priors (see Robert [2007) $3.3.3] and Appendix A). We explicitly differentiate between conjugate families and conjugate families with closed-form normalising constants. We will consider two distinct types of conjugate families: one for classes of squared families, and one for their corresponding base measures. Both conjugacies are with respect to the same likelihood function.
Gaussian processes Gaussian processes (GPs) are flexible nonparametric models for probabilistic inference over functions. Typically one specifies that some function of interest priori follows a GP, takes some observations from the likelihood given the value of the GP, and then computes the posterior conditioned on the observations under the likelihood. If the likelihood is Gaussian, then the posterior process is also GP because the Gaussian distribution is a conjugate prior for the Gaussian likelihood. There are multiple ways to mathematically construct GPs, two of which are the function space view and the weight space view [Rasmussen and Williams 2006 $ 2.2 and $ 2.1].

---

Qu.A
M
Qt- M.P.j
M
qC M0-M,,"'
Qu M.
P 1-3
QM.A
Figure 2: Given base measure μ and feature sb, parameter M M indexes prior Q( M, n4 1) from generalised squared family (GSF) QM. (see Definition D. Given some observed data U and likelihood p(U 1)), the same parameter M E MI indexes the posterior Q(- M, v, v)). where p(U -)((-) (see Proposition 3). The posterior and prior both belong to the conjugate class QMA... a GSF (see Definition DT In certain cases of linear regression in feature space, where the prior is placed over coefficients and the likelihood is Gaussian with expectation described by the regression function, the same parameter M E MI indexes the posterior predictive distribution Q( M, 1", ) for some appropriate v" and rp² (see Corollary 6). The posterior predictive, posterior parameter, and prior distributions all belong to the conjugate class @M, also a GSF. Closed-form expressions for all densities in the GSF QM, as well as the marginal likelihood (see Proposition allow for joint learning of M, Jl and ap of the prior distribution via type maximum likelihood, empirical Bayes (see $ and for closed-form posterior updates.
2
Function-space view A (vector valued) Gaussian process (f(x)) simply is collection of (possibly infinitely many) random dz-dimensional vectors such that every finite subset {f(x,)) follows an Nddzvariate Gaussian distribution. GP is characterised by its vector-valued mean function m : X + dR and positive semidefinite PSD matrix-valued covariance (or kemel) function K X x X + 54, where $ is the set of de x d2 PSD matrices. We write ~ GP(m, K) to mean that f is a GP with mean and covariance functions m and K respectively, meaning that for any finite set of indices X = [x, J vec (f(X)) E R Ndg has mean matrix vec (m(X)) E R Ndz and block partitioned covariance matrix with ijth block K (x;, m.y) sd..
Weight-space view Fintte-feature GPs can be constructed by specifying feature mapping Y- defining a function f(x: (2=) and choosing a Gaussian distribution over the parameters 12 Let Ω E Rd xda follow a matrix Gaussian distribution such that vec(S2 has mean zero and covariance B 3 C' with B E s and C E sit with the notation & denoting the Kronecker product. Let Y X + gdd be a finite feature mapping and define I [==((}) 7(X). Since Gaussian random vectors are closed under addition, one may verify that vec (f(x)) vec (T22T) follows Ndg-variate Gaussian distribution. The covariance matrix is given by
E vee vee ) B (TCTT).
We note that this is a special-case way of handling vector- vector-valued GPs, corresponding with the so- called coregionalisation model [Alvarez et al. 2012]. Our later discussions likely extend more involved vector-valued setups, with more involved notations.
Conditionals and marginals of GPs One attractive property of GPs is that marginal, conditional and evaluations of linear operators of GPs are also GPs. This in particular implies that if one uses a GP function prior or Gaussian parameter prior and one makes observations from Gaussian likelihood, the posterior parameter distribution (if applicable) and the posterior predictive distribution are also Gaussian and available in closed-form Furthermore, the marginal likelihood is available in dooed--orrm and Gaussian. Access to the marginal likelihood allows for empirical Bayes, or type II maximum likelihood for estimation of hyperparameters of the kernel and mean functions. When combined with deep learning, this approach is known as deep kernel learning [Wilson et al.) 2016). This allows for deep neural network features to be used end-to-end in probabilistic regression models. Unfortunately, while expressive deep learning features are used, which leads to improved point estimates, the uncertainty itself is limited to unimodal Gaussian distribution.
3

---

Generalised squared families Let v/ f + IR be some feature mapping, and let μ be measure supported on D2 We construct densities g(w e. M. 4) (with respect measure p)) which are proportional to the squared Euclidean norm of (w)) That
g(w O. (.)
( =
(wb(w) dulww).
(
The normalising constant z(0) satisfies convenient parameter-integral factorottaaa factorisation
Kp
= Tr ( @ ( (ww))))(()) a(dws)
which follows from the cyclic property of the trace and linearity of the integral. The squared family kernel K p.pp is a PSD matrix, as is O O, and the normalising constant may be viewed as an inner product of PSD matrices. The squared family kernel is available in closed-form for many combinations of f2, and ab, with examples from the neural network Gaussian process (Neal 1995 and neural tangent kernel literature [Han et al. 2022 table 1]. random feature literature [Rahimi and Recht 2007]. and classical exponential families Nielsen and Garcia 20224] (see Appendix B) The
parameter factor + + involves no integration, and the integral factor K K does not depend on any
parameters. This factorisation extends beyond the normalising constant to other integrals (such as
those involved in the marginal likelihood, posterior parameter, and posterior predictive distributions)
and is very helpful, allowing factorisation into + + and parameter- innpeeeddenn integral.
The numerator |( ||2 Tr (e* (bb()) may also be understood as inner product of PSD parameter matrix and rank 1 feature matrix. Hence may also parameterise g(w O μ4, 1b)
and its corresponding probability measure Q(dw M, 14, 1b) in terms of M Q O
Tr (M (W))) p(dw),
Q(dw M, p. 1b)
(M) = Tr (MK,,.) w.v).
(1)
x(M)
Distributions of the form generalise mixture models, which are obtained when M is diagonal.
3 Conjugacy of Generalised Squared Families
In order to better describe convenient properties of Bayesian updates, introduce various families of probability densities. We later describe precisely how in some sense these families are closed under Bayesian updating, and how this closure is computationally attractive.
Definition 1. Let M. be set of x PSD matrices, & be a set of nonnegative measures, and be sei of feuture mappings of the form sb : O + RD. A generalised squared family (GSF) is set of probability measures of the form ①
QM. [Q(dw M,
p1, MEM.NERNEE
indexed by M E M., j 2, E V. We allow any of the indexing sets ML, B. v to be singleton, and with an abuse of notation, write the singleton's 's element in place of the set. For example,
QM (Q(dw M. (.) MEE and QM. (QQ M,,, +b) PEA. (.EAAE)
We note that by Von Neumann's 's trace inequality, M. may be taken to be the set of all non-zero PSD matrices long as the squared family kernel K p. is PD, since Tr (MK, vn 2(M > 0. The inclusions QM, & S QM. S QMAA., together with a summary of some of our later results, are visualised Figure 2 Previous work (Tsuchida et al., 2025) called families of the form QM where M. is the set of rank 1 PSD matrices, squared families. Throughout this paper, we use Q for probability measures belonging to GSFs, and and p for arbitrary probability measures and densities (which may also belong to GSFs, where stated).

---

M. u. 4
r(t,)
M.
Wij
w
M.
1- L. N
w Q(dw M, p(u w)
y,
w
Qidwr M.
| M.P.. p(u w)
p(y (Ra(s))
Figure 3: (Left) 3.1 we consider an abstract model where random variables following GSFs w are latent under an observation model we expand the observed node u so that observations take the form of (æz, y) pairs, and assume an observation model of the form y ~ p(y [2y())) Here w = vec is parameter and y is some deterministic feature mapping (for example, from deep neural network). Finally in $3.3 on the setting where the likelihood pyy (2)(z)) and base measure are both Gaussian, allowing us generalise featurf Gaussian process regression. (Right) Our results extend to the case of hierachical model with shared prior over J datasets. We apply this to few-shot learning in $ 4
p(u w). (Middle) In & speciallse setting where
focus
3.1 Marginal likelihood and posterior distributions
We first consider abstract setting where there two variables of interest variable which observe and a latent variable which we do not observe but would like to infer. The latent variable w
priori is distributed w ~ ((w M, ju, tb) according an element of a GSF, and is mapped to the observed variable U = {u.),- implicitly through some arbitrary likelihood function p(u w).
This model is depicted in the left part of Figure
The main trick in performing posterior updates GSFs that the base measure and the likelihood combine to give a new base measure. We derive the marginal likelihood in the main text to expose the simplicity of this idea. Here we make use of the parameter-integral factorisation. The proofs of later results follow a similar manipulation (see Appendices C
Proposition 2. Consider the left model of Figure] 3 Define a measure v(dw) = p(U ( as the product of the likelihood and base measure. The marginal likelihood is the ratio of normalising constants, p(U) Tr (M (M.Kw))/ Tr (M (MK, K w.b)) (where K K depends U through v).
Proof. The marginal likelihood p(U) is given by
v(dw)
1. p(U [w) M,
Tr M In (w)((())) P(U w))(dw)
Tr (M K v. Tr (MK m.
Tr (MK m.
The marginal likelihood does not belong to GSF The posterior distribution does belong to GSF, even for finite N. We may form a GSF which includes both the prior and the posterior, with an updated base measure and unchanged parameter.
Proposition 3. Consider the model in Figure (Left). Define measure v(dw) p(U ((ww. The posterior P(dw U) = Q(dw M,v, v, v)) belongs to a GSF with base measure v. feature mapping 2a and parameter M.
2 Informally, while the marginal likelihood itself does not belong to GSF, it asymptotically belongs to GSF, since the measure is proportional posterior under prior μ. Therefore under mild conditions, by the Bernstein-von Mises theorem. y concentrates a Dirac delta distribution centered at MLE wv. Hence p(U)/ ((dw) z Ir (M (w n (w>)w)))) ) / 'r [M K ). ). This would allow previous results on maximum likelihood estimation with GSFs Tsuchida et al, 2025 to extend to maximum marginal likelihood estimation. We leave rigorous investigation of this observation for future work.
5

---

Proposition implies that the GSF QM... is a conjugate class for likelihood p(u w), where the of all measures that can be expressed in the form (I p(u; ((w)) for some integer h 2 0 and belonging to the support of the likelihood. We note that the posterior predictive p(U = U) is also available, since is essentially of the same form as the marginal likelihood, replacing the GSF prior with a GSF posterior and the training likelihood with the test likelihood.
Proposition Consider the model in Figure 3 (Left). Define a measure v(dw) = p(U ())u(dw) and (dw) = p(U, w)v(dw)) where, abusing notation, the train likelihood p(U w) and test likelihwod pl U. w) are not necessarily the same. The posterior predictive is the ratio of normalising constants, pU,, U)= Trr (MKv v' w/) Tr (M..)))
The following observation illustrates that it is easy to obtain closed-form expressions for the marginal likelihood, posterior parameter and posterior predictive distributions, whenever the base measure is itself conjugate to the likelihood (for exponential family examples, see Appendix A
Remark 5. If A a conjugate class for p(u | w), then v E (up to a normalising constant), so if the squared family kernel K &. is available in closed-form all A, then so is v.
Note that here refer to two different conjugate updates: One if the prior were base measure belonging to A and another for the true prior belonging to a GSF QM. As described in $ 2 and Appendix B closed-fome expressions for the squared family kernel
available for many of
combinations of 0, Jr and 1a [Neal] 1995
Williams
1996 Rahimi and Recht [007 Cho and Saul
2009
Nielsen and Garcia 2024 Han et
2022 Tsuchida al. 2023 for example].
3.2
Regression models
We now consider a more specific instance of the previously described model, as shown in the middle of Figure 3 We consider Bayesian learning of vector-valued functions of the form
(2) where _2T E RRdxxds are readout parameters and y : X - Rd is a feature mapping (from e.g. a deep neural network). We assume our observations are independent samples (X, Y) {(.,, )
f(x:) = 22y(ar,
w ~ Q(dw M, 14, 3b), w = vec(S2) E Rdd where d= did2 d2,
from ply (ty((=)) and some distribution over Let r 7X)). The marginal likelihood,
p
rJ) = P(Y Io )((ww)
M, p, v,,)
serves two purposes. First, as a normalisation constant for the posterior, and second, as an objective function for point estimation of M, μt, and sb. The posterior parameter distribution is given by
P(d Y. Y.T) p(Y rn (dw (1,5))))
P(Y |T)
which in turn induces a posterior over noisy evaluations Y, of f new test points (features) Xx (T+). This posterior is called the (noisy) posterior predictive distribution,
P VY..T.)...).
p(Y. I...T)P(ddS) Y,T),
and is used to form predictions on test data X.. Note that the test likelihood p(Y, |T may be different to the train likelihood p(Y [T e.g. with heteroskedastic noise. The noiseless posterior predictive p(f(X.) Y.I,T.) is obtained whenever the test likelihood has zero variance.
Marginal likelihood, posterior parameter, and posterior predictive distributions The results of Propositions 2 and
4 apply directly to the marginal likelihood, posterior parameter and posterior p(Y rΩ )u((w) and v' v'(ww) = p(Y. T.Ω )d(dw), the marginal likelihood is p(Y Y T, MM Tr (MK,, / / Tr (MK p. the poste- rior parameter distribution P(a Y. T) Q(dw M. v, tb) belongs to a GSF and the posterior predictive distribution is p(Y. T, mark 5 still applies; if μ., and v' remain inside the same conjugate class and K K. is available within that family, then closed-form normalising constants are available. We utilise this observation in the following subsection, focusing conjugacy of the Gaussian family with Gaussian likelihood.
predictive distributions. That is, defining v(dw)
Tr (MK,,.) Tr (MKr*). Note also that Re-

---

3.3 Generalised squared family processes
We now specialise our choice of base measure and likelihood. Take the Radon-Nikodym derivative of the base measure j with respect to Lebesgue measure to be matrix Gaussian,
dp (2) - p(2)) - NV(w ve &
(3)
Take the train likelihood to be multivariate Gaussian with homoskedastic noise with variance v²,
p(Y rQ (Vvooo(V) (Id, @ T)w,"?INa
Take the test likelihood to be multivariate Gaussian with homoskedastic noise with variance y?? PYY r. NV (vec(Y (1 T.)w, I.),,,"IINd) IN.IN.
(5)
With GSF prior, we refer to such model as GSF process (GSFP).
Completing the square and parameter-integral factorisation By Proposition 4 and Gaussian conjugacy, in order to compute the posterior predictive distribution, we need only reverse the factorisation z/ vd(w)
p(Y, r.Ω )p(Y ru (p(dw) as v' (dw) sulw Y, VYY.))p (YY, YY)((Y, which is Gaussian, and same for v(dw)) p(Y In )u(dw) The mean and covariance parameters of these Gaussian measures can be found using the same techniques as in vector-valued Bayesian linear regression. We show only the noise free posterior predictive distribution here, with more details in Appendix [D]
Corollary 6. Consider the special model
3 and 5 applied to Figure
(Right). Define
v' vdww) = pYY, r.n. )p(Y rΩ (2) and
myy vec(M) + (B CT (B % ICT
2² INd.) (oNd))=(vveee)) (vec(Y) (Ids @ I) vec(M))
B C (B s CT (B @ rCTT 2 I so rc).
Suppose that C' and F, are full-rank and N. 2 d1. Then the noise free predictive distribution ((df(X.) Y, r.). mapping and modified base measure,
(d) M. ,b)) belongs to GSF with modified feature
vf(((.)) & 2/ (Loo & (r)))((.)) xb vec (oof(f(((.)))
v((f((.)) - N vec (f(x.)) (Idz & T.)mm, (Ta, @ C @ (T.) (6)
Corollary 6 implies that under the Gaussian likelihood and base measure setting, the GSF QM.A. is a conjugate class, where @ is the class of transformations obtained by composition of vp with linear transformations and A involves deformation via Gaussian conjugate updates, as per (6).
3.4 Related work
To the best of our knowledge, we are the first to consider squared families as priors and posteriors. The nature of the results concerning likelihoods is fundamentally different to that of priors and posteriors, because for example, one does not attempt to estimate parameters in forming a posterior, one merely attempts to form posterior belief. Nevertheless, existing work special cases of squared family likelihoods is relevant in that expressive power, flexibility and tractability of densities intuitively relevant to both likelihoods and posteriors, and uncertainty quantification more generally.
Kernel methods Building off the general results of modelling non-negative functions using squared elements of RKHS (Martcau-Ferey et al. 2020]. Rudi and Ciliberto [22T] model probability distributions, with tractable marginalisation and closure under multiplication. They fit likelihood models within this class by minimising a regularised L2 distance between the target likelihood and the model. While these models are nonparametrie in construction, estimates resulting from minimising the objective satisfy representer theorem type instantiation in finite dimensions. The authors show that the model admits both favourable representational capability (in being able to represent f-times differentiable target density with error less than with number of parameters scaling like ((c-r/a (logr/((og(( //)))//2, where is the dimension of the domain of the distribution), as well as statistical estimation properties (as likelihoods, the proposed estimate converges at rate of N 2в47 (log //2, where N is the number of datapoints).
7

---

Probabilistic circuits Probabilistic circuits [Choi et al. 2020) provide a computational framework for tractable probabilistic modelling. They are constructed as graphs of connected units, and by constraining the graph, allow for tractable marginalisation. However, they need to impose some structure on the units or the connections circuits [Sladek] 2023 Loconte et al., 2023a b Loconte and Vergari 2025] bypass this constraint by squaring the output of the circuit. Sums of squared circuits can also be instantiated [Loconte et al.) indeed quadratic time/memory in the number of units) computation of normalising constants or marginalisation, as well as the inclusions of the function spaces imposed by different classes of probabilistic circuits.
order to ensure nonnegativity. Squared probabilistic
2024). The focus such works is typically in showing tractable (polynomial
Neural nefworks and finite feature models Previous work (ssuchida et al.) 2025 has modelled likelihoods as proportional to the squared Euclidean norm of a function of the form 2 (z), where O a matrix and & vector-valued funetion. Provided the features s/ are rich enough, such models admit a universal approximation property, typically approximating the squared L2 distance between a square root of a target density and the model at a rate of (11/n), where n is parameter count. As likelihoods, they learn target densities with KL divergence to the target density decreasing at a rate of O(1 ((1//N), where N is the number of samples. Such families of models have tractable normalising constants, Fisher information and statistical divergences. Furthermore, in some special cases, such families have not only tractable but indeed closed-form normalising constants, Fisher information and statistical divergences. An example studied previously was squared neural families (Tsuchida et al. 2023]. where is a single hidden layer neural network.
Applications Poisson point processes (PPPs) are helpful variations of density modelling, where one models an intensity instead of density. Whereas in density modelling, realisations are i.i.d. draws from a target density. in intensity modelling, realisations are conditionally i.i.d. draws from target density given the number of realisations. The number of realisations follows a Poisson distribution with an expected number of points equal to the integral of the intensity over the domain. Flaxmman et al.. [2017) used (frequentist) squared elements of RKHSs to model intensity functions, whereas Walder land Bishop [2077] used a (Bayesian) squared Gaussian process prior to model intensity functions. We note that placing prior over a function which when squared gives intensity/density (likelihood) is completely disjoint to the problem we are solving here, which is to use a squared function as a prior density. Probabilistic circuits have been applied in converting knowledge graph embeddings into generative models [Loconte et al., 2023a] Squared neural family likelihoods have been applied to label distribution learning [Zhang al. 2025], which is kind of variation classical regression where instead of the target label being a single class, the target variable is a discrete distribution (i.e. element of the simplex). Hence, here squared neural family models model distribution over distributions with a closed-form normalising constant, expectation, variance, and covariance conditioned on input samples. Probabilistic modelling is leveraged to provide confidence intervals, conformal predictions, active learning, and model ensembling.
4 Experiments on few-shot learning
Few shot learning setting The hyperparameters of nonparametric probabilistic model (such as GP or GSFP) can adapted to a single training dataset (X, Y) with N examples, by maximising the marginal likelihood p( Y X) under the model with respect to the model hyperparameters, also referred to as type II maximum marginal likelihood [Rasmussen and Williams) 2006 5.44. It may also happen that meta dataset, consisting of J datasets {(X,, j, Y J is assumed to arise as J samples from a model sharing the same hyperparameters, in which case the marginal likelihood (x,)²= 1) still serves as a natural objective for tuning model hyperparameters (Ras mussen and Williams 2006 p. 115]. In this case, each of the (X, j, YY;) consists of N, examples, where N, may be constant N across all J datasets, but not necessarily, At inference time, when a new support dataset (X; j. Y and new set of query inputs X are given, one may form the posterior predictive p(Y; X; Y1) using the hyperparameters obtained by maximising the marginal likelihood of the tasks. This classical setting arises from point estimation of the hyperparameters un- der the assumption of shared point hyperparameters (see right side of Figure 33. and has traditionally been referred under the umbrella of multi-task learning [Caruana| 1997 Minka and Picard. 1997].
P(DY))
This setting has recently received renewed attention under the name of few-shot learning, and GP models have been empirically shown to outperform other more recent models on modern benchmarks,
8

---

Figure 4: When can multimodal prior be helpful? (Top) GP using deep learning features (deep kernel learning [Wilson et al. 2016)). (Bottom) GSFP using deep learning features. (Leff) Prior predictive, trained using points (white) drawn from a multimodal distribution over functions (with mean given by red and blue). The GFSP is bimodal, whereas the GP is always unimodal. (Middle) Posteriors for both models given some ambiguous points, which may belong to either red or blue modes, The GP incorrectly places all of its belief over one of the modes, whereas the GFSP reserves belief for both modes. (Right) Given additional data from x 0, both identify the correct mode.
both in the regression and classification setting PPatacchiola et al., 2020 (Deep kernel transfer, or DKT). In the regression setting, further enhancements to GP models, which use continuous normalising flows to transform GP marginal predictive posterior distributions into more complicated distributions, have further been used to improve upon GPs [Sendera et al. 2021] (Non-Gaussian Gaussian processes, or NGGP). Such enhancements address GP prior limitations; namely lack of an ability probabilistically model skewed, non -ann-GGssiitaaan tailed multimodal beliefs. Our GSFP also offer a means to address such limitations (see Figure 4), and are therefore an excellent candidate for few-shot learning. We consider the same scalar-valued (d2 1) regression benchmarks as Senderal et al.
2021), who in turn add to the regression benchmarks of Patacchiola et al. 2020).
Implementation details We learn prior predictive distribution, which is further broken into prior parameter distribution belonging GSF QM.A, + and (deterministic) feature extractor y, by attempting to maximise the marginal likelihood (see Proposition with respect to some tuneable decision variables. The decision variables are an unconstrained O, with O Q M and feature mapping & (which together define the prior over S in (2), as well as the tuneable parameters of the neural network 7. We further parameterise the feature mapping &b as the hidden layer of learnable neural network, (b(w) = a [Ww w b). with decision variables W and Updates within the set A happen automatically as part of the Bayesian inference pipeline and are not part of the decision variables the marginal likelihood (see Proposition 2). We reimplemented competitor models DKT [Patacchiola et al. 2020] and NGGP [Sendera et al. 2021]. as we were not able to get the originally supplied code to run, and through extensive hyperparameter tuning were able to obtain results better than what they originally reported in all but two datasets (NDX 100 and EEG). Full hyperparameter details and run times are given in Appendix IF] Results are shown in Table 4 We observe that GFSP obtains the lowest NLL on 11 out of 12 settings, although sometimes is well within one standard deviation of the best competing method. While our method is designed for situations in which multimodal belief is helpful, e.g. few-shot learning, also apply it to vanilla regression. See Appendix(C G
Discussion, limitations and conclusion
Limitations The main computational hurdle in working with GSFs computing the (exact) normalising constant via the squared family kernel, which has a time complexity of O(n?d), where n the size of the hidden layer and d the dimension of the random variable. practice, the training time of GSFP is slightly longer than that of DKT, but substantially less than that of NGGP (see Table F... The inference time is much lower than NGGP, but slightly more than DKT. As with competitors DKT and NGGP, which are also based on (potentially finite-fean (niterfeature) GP methods. our

---

Table 1: Benchmark results showing test NLLs for regression tasks across 6 meta datasets. For each meta dataset, we train 5 meta models with different random seed via maximum marginal likelihood over the J datasets {(x,, j. Y })} At testing time, we condition (train) on a support (few shot) dataset (Xi, 1: Yj) and evaluate the test NLL on the query dataset (X 22). We evaluate the test NLL on these few shot datasets for 500 random shuffles of the test meta dataset. In total, 5 x 500 200 models are trained. Each entry shows the mean + standard deviation over 2, 500 model evaluations.
Each testing scenario includes in-range and out-of-range, where the support and query set is either similar to or different to the J training datasets. Best results in bold, second best underlined.
Methods
Paaeecaaacchiooaacom al.
2020
Mixed-Noise Sines Sendera al. 2021
QMUL Gong al. 1996
Name
Kernel
NLLb. 4
NLLoa 4
NLL
NLLoa
NLL1, 1
NLLou.
RAF
0.71 4-0.0
-0.62
0.76
0.06
4 0.13
-0733 0.19
-0.00 0.21
DKT
Spectral
-0.82 1 0.05
-0.79
0.06
0.38
0.17
1.78 1 0.49
-0.67 1 0.24
-0.64 0.20
Linear
0..76 + 0.00
0.30
1.07
0.47 0.25
2.31 0.772
-0.66 0.49
0.87 1.18
RBF
-0.74
-0.56
0.14
0.13
0.87
-0.10
0.13=0.70
NGGP
Spectral
-0.84 0.00
-0.73
0.07
0.38=0.17
1.78 0.49
-0.99 + 0.31
-0.77 0.33
NN Linear
-0.80 + 0.06
0,27 + 1.26
0.21 0.89
2.14 1.32
0.01 = 1.31
1.41 1.27
GSFP
NN Lineur
0.85 1. 0.06
0.83
0.07
0.46 0.10
0.98 1. 0.21
-1.99 1. 0.12
-1.15 0.14
Methods
NDX10U
Qin et al.
2017
EBG
Femandez-Fraga et al. 20199
Power (Hebrail and Berard] 2006
Name
Kernel
NLLn J
NLLoa 1
NLL
NLLoa
NLL., 1
NLLov
RBF
1.16 1: 0.03
-1.17 a 0.03
-1.10 0.54
-1.16 1: 1,09
-0.77 1 0.62
-0.34 1 0.70
DKT
Spectral
-0.46 0-43
-0.63 0.30
-1.32 +.33
0.84 + 0.67
-0.47 + 0.62
-0.35++0.70
NN Lineur
-2.16 1.61
-2.41 a 0.96
-1.10
-1.02 1. 1.06
-0.42 1. 0.68
-0.33 0.773
RBF
1.16
0.15
-0,96 0..17
0.05 0.332
1.15
115=
NGGP
Spectral
-1.32
0.84+00.67
+..30
0.33
1.35+00.01
1.35 0.01
NN Lincar
-2.07 1.20
-2.32 = 0.02
-1.05 0.559
-1.18 0.89
-0.27 0.40
0.:10
GSFP
Lineur
-2.49 1.28
22.5 0.53
-1.49 0.44
-1.38 0.81
-0.71 0.339
-0.65 0.45
method has (N2)) and O(N)) memory and time complexity for inference (or (nN) and (n²)
for fenite-feature models), where is the number of training points.
Computational complexity of Bayesian inference In the general case, Bayesian inference is computationally difficult, unless some heavy restriction is placed on the densities (e.g. log concave). This is especially difficult when the dimension is large, due to the curse of dimensionality and the exponential growth of volume. In practice, if one desires an exact computation, one needs to restrict the class of probability distributions. When we restrict the class to GSFs, we able compute exact posteriors. The complexity of this calculation is governed by the inner product of the parameter matrix M and the kernel matrix K u. (as in I.. Assuming the kernel matrix is known, this inner product has complexity O(z²), where n is the number of parameters in the GSF. Crucially, dependence of the complexity on d only appears through the calculation of the kernel matrix K (L In practice, the kernel matrices K p.pp can be computed exactly in linear time in d. This leads to a combined complexity of Odn)²), i.e. exact inference linear in dimension d, which is huge improvement over inexact inference exponential in dimension d. This improvement comes with the restriction of the class of densities to GSFs, however this restriction is not too severe because GSFs are rich density approximators. More precisely navigating the trade-off between universal approximation (making n large) and computational efficiency (making small) is important direction for future work.
Conclusion GSF likelihoods appear in disguise in a wide range of applications, including Poisson point processes [Flaxman et al. 2017 Walder and Bishop 2017 Tsuchida et al. 2024 generative models [Loconte et al., 2023a], and most fundamentally, density estimation [Rudi and Ciliberto 2021). Here we considered the orthogonal setting of GSF priors. Using the parameter-integral factorisation, we found that GSFs form conjugate priors and that for many likelihoods of interest, admit closed-form Bayesian updates. Such closed-form update allows one generalise Gaussian parameter prior and Gaussian likelihood regression in feature space models to GSF priors, allowing for multimodal and rich expressions of uncertainty. We empirically demonstrated that GSFPs perform better than or as well as other models on a range of benchmark few-shot regression problems. In future. should be possible to show asymptotic normality of maximum marginal likelihood estimates, as well as provide generalisation bounds for the problem of density estimation under a hierarchical model, by using the fact (see Footnote 2 that the marginal likelihood asymptotically belongs to a GSF and using known results for GSF likelihoods.
10

---

Acknowledgements
DS was partially supported by the Responsible Al Research Centre (RAIR).
References
Mauricio A Alvarez, Lorenzo Rosasco, Neil D Lawrence, et al. Kernels for vector-valued functions: A review. Foundations and Trends@ in Machine Learning, 4(3):195 266, 2012.
Christopher M Bishop. Pattern recognition and machine learning, volume 4. Springer, 2006.
Lawrence D Brown. Fundamentals of statistical exponential families: with applications in statistical decision theory. 1986.
Rich Caruana. Multitask learning. Machine learning, 28:41-75, 1997.
Youngmin Cho and Lawrence Saul. Kernel methods for deep learning. Advances in Neural Information Processing Systems, volume 22, 2009.
YooJung Choi, Antonio Vergari, and Guy Van den Broeck. Probabilistic circuits: A unifying framework for tractable probabilistic models. 2020.
Andrea Coraddu, Luca Oneto, Alessandro Ghio, Stefano Savio, Davide Anguita, and Massimo Figari. Condition Based Maintenance of Naval Propulsion Plants. UCI Machine Learning Repository, 2014. DOL: htpps//bbi.og/11222445////3
Peter Corke. A robotics toolbox for matlab. IEEE Robotics & Automation Magazine, 3(1):24-32, 2002.
Paulo Cortez, António Cerdeira, Fernando Almeida, Telmo Matos, and José Reis. Modeling wine preferences by data mining from physicochemical properties. Decision support systems, 44)) 547-553, 2009.
Persi Diaconis and Donald YIvisaker. Conjugate priors for exponential families. The Annals of statistics, pages 269-281, 1979.
S.M. Fernandez-Fraga, M.A. Sevves-Fernndde, and J.C. Pedraza-Ortega. Eeg data collection using visual evoked, steady state visual evoked and motor image task, designed brain computer interfaces (bci) development. Data in Brief, 25:103871, 2019, ISSN 2352-3409. doi: cocioree/arriaaa//////5555349999
URL https://wwrt wciencedirect............. com/
Chelsea Finn, Pieter Abbeel, and Sergey Levine. Model-agnostic meta-learning for fast adaptation of deep networks. In Doina Precup and Yee Whye Teh, editors. Proceedings of the 34th International Conference on Machine Learning. volume 70 of Proceedings of Machine Learning Research. pages 1126-1135. 115.. PMLR, 06-11 Aug 2017. URL https://procedaag..ttt paw..pprsss///// finn17a. html
Seth Flaxman, Yee Whye Tch, and Dino Sejdinovic. Poisson intensity estimation with reproducing kernels. In Artificial Intelligence and Statisties, pages 270-279. PMLR, 2017.
Gerritsma, R. Onnink, and A. Versluis. Yacht Hydrodynamics. UCI Machine Learning Repository, 1981. DOE:
S. Gong, S. McKenna, and J.J. Collins. An investigation into face pose distributions. In Proceedings of the Second Intermational Conference on Automatic Face and Gesture Recognition, pages 265-270, 1996. doi: 10.1109/AFGR. 1996.557777
Richard HR Hahnloser, Rahul Sarpeshkar, Misha A Mahowald, Rodney J Douglas, and H Sebastian Seung. Digital selection and analogue amplification coexist in cortex-inspires innspired silicon circuit. Nature, 405(6789):947-951. 2000.
Insu Han, Amir Zandieh, Jaehoon Lee, Roman Novak, Lechao Xiao, and Amin Karbasi. Fast neural kernel embeddings for general activations. In Koyejo. S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 35, pages 35657-35671, 2022.
David Harrison Jr and Daniel I. Rubinfeld. Hedonic housing prices and the demand for clean air. Journal of environmental economics and management, 5(1):11-002, 1978.
11

---

Georges Hebrail and Alice Berard. Individual Household Electric Power Consumption. UCI Machine Learning Repository, 2006. DOI:
José Miguel Hernández-Loban and Ryan Adams. Probabilistic backpropagation for scalable learning of bayesian neural networks. In International conference on machine learning, pages 1861-1869.
PMLR, 2015.
Diederik Kingma. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.05 2014.
Lorenzo Loconte and Antonio Vergari. On faster marginalization with squared circuits via orthonor- malization. In AAAI 2025 Workshop on Connecting Low-rank Representations in AI, 2025.
Lorenzo Loconte, Nicola Di Mauro, Robert Peharz, and Antonio Vergari. How turn your knowledge graph embeddings into generative models. In Thirty-seventh Conference on Neural Information Processing Systems, 2023a. URL hRttsttps::
Lorenzo Loconte, Stefan Mengel, Nicolas Gillis, and Antonio Vergari. Negative mixture models via squaring: Representation and learning. In The 6th Workshop on Tractable Probabilistic Modeling, 2023b.
Lorenzo Loconte, Stefan Mengel, and Antonio Vergari. Sum of squares circuits. arXiv preprint 1rXiv:2408.11778, 2024.
Ulysse Marteau-Ferry, Francis Bach, and Alessandro Rudi. Non-parametric models for non-negative functions. Advances in neural information processing systems, 33:12816-12888 2020.
Thomas Minka and Rosalind W Picard. Learning how to learn is learning with point sets. Unpub- lished manuscript. Available at htlps://tminka. htttimiiikniiittuuuuo. httpaubb.oooppprs/piitssesss/s/e
1997.
Radford M Neal. Bayesian learning for neural networks. PhD thesis, University of Toronto, 1995.
Frank Nielsen and Vincent Garcia. Statistical exponential families: A digest with flash cards. arXiv
preprint arXiv:091 1,4863, 2024.
Massimiliano Patacchiola, Jack Turner, Elliot J Crowley, Michael O'Boyle, and Amos J Storkey. Bayesian meta-leaming for the fho-shot setting via deep kernels. Advances in Neural Information Processing Systems, 33:16108-161 18, 2020.
Tim Pearce, Russell Tsuchida, Mohamed Zaki, Alexandra Brintrup, and Andy Necly. Expressive priors in bayesian neural networks: Kernel combinations and periodic functions. In Uncertainty in artificial intelligence, pages 134-144. PMLR, 2020.
Yao Qin, Dongjin Song, Haifeng Cheng, Wei Cheng. Guofei Jiang, and Garrison W Cottrell. A dual-stage attention-based recurrent neural network for time series prediction. In Proceedings of the 26th International Joint Conference on Artificial Intelligence, pages 2627-2633, 2017.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. Advances in neural information processing systems, 20, 2007.
Prashant Rana. Physicochemical Properties of Protein Tertiary Structure. UCI Machine Learning Repository, 2013. DOI: http:://6ol.org11024433///
Carl Edward Rasmussen and Christopher KI Williams. Gaussian processes for machine learning, volume 2. MIT press Cambridge, MA, 2006.
Christian Robert. The Bayesian choice: from decision-theoretie foundatiins to computational implementation, volume 2. Springer, 2007.
Alessandro Rudi and Carlo Ciliberto. PSD representations for effective probability models. In Advances in Neural Information Processing Systems, volume 34, pages 19411-19422, 2021.
Marcin Sendera, Jacek Tabor, Aleksandra Nowak, Andrzej Bedychaj, Massimiliano Patacchiola, Tomasz Trzcinski, Przemyslaw Spurek, and Maciej Zieba. Non-gaussian gaussian processes for few-shot regression. Advances Neural Information Processing Systems, 34:10885-10298, 2021.
Aleksanteri Sladek. Positive Semi-Definite Probabilistic Circuits. Master's thesis, Aalto University.
School of Science, 2023.
Athanasios Tsanas and Angeliki Xifara. Accurate quantitative estimation of energy performance of residential buildings using statistical machine learning tools. Energy and buildings, 49:560-567, 2012.
12

---

Russell Tsuchida, Tim Pearce, Chris van der Heide, Fred Roosta, and Marcus Gallagher. Avoiding kernel fixed points: Computing with ELU and GELU infinite networks. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 9967-9977, 2021.
Russell Tsuchida, Cheng Soon Ong, and Dino Sejdinovic. Squared neural families: A new class of tractable density models. In Neural Information Processing Systems, 2023.
Russell Tsuchida, Cheng Soon Ong, and Dino Sejdinovic. Exact, fast and expressive poisson point processes via squared neural families. In The 38th Annual AAAI Conference on Artificial Intelligence, 2024.
Russell Tsuchida, Jiawei Liu, Cheng Soon Ong, and Dino Sejdinovic. Squared families: Searching beyond regular probability models, 2025. /attpss//arxiv..rrrreeetrs abs 21128
Pinar Tüfekci. Prediction of full load electrical power output of base load operated combined cycle power plant using machine learning methods. International Journal of Electrical Power & Energy Systems, 60:126-140, 2014.
Christian J Walder and Adrian N Bishop. Fast Bayesian intensity estimation for the permanental process. In International Conference on Machine Learning, pages 3579-3588. PMLR, 2017.
Christopher Williams. Computing with infinite networks. In Advances in Neural Information Processing Systems, volume 9, 1996.
Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P. Xing. Deep kernel learning. In Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, volume 51 of Proceedings of Machine Learning Research, pages 370-378, PMLR, 2016.
I-C Yeh. Modeling of strength of high-performance concrete using artificial neural networks. Cement and Concrete research, 28 1797-1788 1998.
Daokun Zhang, Russell Tsuchida, and Dino Sejdinovic. Label distribution learning using the squared neural family on the probability simplex. In The 41st Conference on Uncertainty Artificial Intelligence, 2025.
Liu Ziyin, Tilman Hartwig, and Masahito Ueda. Neural networks fail to learn periodic functions and how to it. In Advances Neural Information Processing Systems, volume 33, pages 1583-1594, 2020.
13

---

A Conjugacy of exponential families
A.1 General formulation
Here we review the well-known conjugacy exponential families. Let t -> 12 with C Rai and let H be nonnegative measure. Let A be such that
P(dæ 17) exp (1 t(x) A(7)) H (ddx)
integrates to 1. The set {P(dæ mn 10< exp(rt <d)) called an exponential family, conjugate family for an exponential family likelihood has densities the form
exp (x 'n r.((n))
p(n x.r) J ex (x rA())) dn
x exp (-p(-r(())) exp(x Im)
That is, conjugate family is {p(n X, xxr))xeR=> form of exponential family with canonical parameters n = (x,r) and sufficient statistic t(n) (n. A(n)). The conjugate prior has the exponential family form
Note that the right expression is in the
p(n m exp (TT(()) A())
À ( log / exp ( (m)) dn.
Given iid observations X E R.Nxe xax from the likelihood, the posterior is obtained by updating the parameters of the prior with the sufficient statistics,
p(n X, X. = n Σ t(x,), +
A.2 Examples
Poisson-Gamma Take X = {0. 1. I....). h(x) and t(x) = x. This results in a Poisson
likelihood, with rate X exp.) and A(r) exp(1)). The conjugate family is set of distributions of the form
pln r) exp (A(-rn)) exp(xn) exp c-p()))))()) exp(xn)-
Taking À exp v, and accounting for the Jacobian of the transformation, we identify the family of densities of the form
p( X, Xx-1 exp(-rA),
which are Gamma densities with rate and shape X- These have a known (divisive) normalising constant given by T(x)/rx, where here I denotes the Gamma function. The posterior density is then
pA X x,r) E ((m:))1
exp(-(r ( N)A),
with known normalising constant. The corresponding prior distribution over 7) = log X, as an exponential family, is the esponential-gamma distribution, with
x-
p(n X X,
exp (x1) exp((rexp(1)).
T(x)
An example of Poisson-Gammi updates is shown in Figure
Caausian-Gaussian Take X R. h(x) exp (- x² (0-x=/(25))) ((2s)) (2 and t(x) x/s for some known standard deviation s 0. This results in Gaussian likelihood, with known variance s?. mean 87) and A(r) x/2 /2. The conjugate family is set of distributions of the form
p(n x,r) & exp(-rx² /2) exp(xn),
which are Gaussian densities with variance 1/r and mean x/r. These have known (divisive) normalising constant, and
-2 exp (- 2r-1 (n x//)²)
p(n x.
(-2r-1
14

---

Truta
Truth
Prior
Posturior N-20
P(d X
Posterior N=20
Posterior N-10
Preterioe N-40
Posterior N-00
Preterioe N-60
Preterior N-S0
Posterior N
Poeterior N-100
Posterior N- 100
22
0
10
7
Figure 5: Bayesian inference in a well-specified Poisson model, using a Gamma conjugate family for rate (or, in canonical parameterisation, an exponential -roonnnial-Gaamaa conjugate family for log X). Shown are the posterior distributions for the parameter given N 0, 20..... 100 samples from the likelihood.
Truta
Truth
Price
Peior
Preterior N-20
Posterior N-20
Posterior N=40
Postecior N=
Poeterior N=
Postecior N=60
Posterior N-80
Posterior N-80
Posterior N-100
Pusterior N-100
1.0
0.5
135
1.00
1.5
23
m
Figure 6: Bayesian inference in a well-specified Gaussian model (with known variance), using Gaussian conjugate family for mean (in canonical parameterisation, the conjugate family is also Gaussian). Shown are the posterior distributions for the parameter given = 0, 20...., 100 samples from the likelihood.
Under the change of variables m = for mean m,
p(m x.r) /Nxxxxx and
pmm X, x. = =N N((x+(1))) I [x)ss((+ m
N
An example of Gaussian- Gaussian updates is shown Figure
15

---

B Examples of squared family kernels
B.1 Exponential family base measures and exponential features
Take the base measure μ belong to any exponential family with sufficient statistic è. Take the feature ap to be (b(w) = exp (Wi(w)/2) Then the ijth entry of the squared family kernel K n.3b is
(K nw J. exp (ww: w,) ((w)/22 T u(dw),
where un is the ith row of W. This is simply the moment generating function of the sufficient statistic t under the exponential family, which is available in closed-form in terms of the log normalising constant (7) of the exponential family,
K p. exp (A + (wn; + w2))2) - A(r)
where 77 is the canonical parameter of the exponential family. Recall Remark that if yu itself belongs conjugate prior family for the likelihood p(U [-), then v(-) = (p(())p( conjugate family, up to multiplication by the marginal likelihood, and therefore the squared family kernel K n. is available in closed-form for all posterior updates. By Appendix A I does indeed belong to a conjugate family. Log normalising constants are available for many exponential families which are known be conjugate priors for exponential family likelihoods (Diaconis and Yivisaker 1979] (see e.g. Appendix A and closed-form normalising constants in Nielsen and Garcia (2224]
remains within the
B.2 Normalised base measures and random Fourier features
Bochner's theorem says that every stationary kernel is the Fourier transform of nonnegative probability measure, and vice versa [Rahimi and Recht [207]] for example]. Hence for real-valued stationary squared family kernels and probability base measures μ4,
(K k,w)sss exp (i(w; w;)² Tw) u(dw)
= L wi (W))))((W))((())
where (w) cos( w.TTT. w) + isin(w.) ( w) and i is the imaginary unit. This allows the use of base measures with cosine and sine features, provided the base measure has closed-form Fourier transform. For example, the Gaussian, Laplace and Cauchy base measures all have closed-fon kernels (corresponding with Gaussian, Cauchy and Laplacian kernels respectively). Real-valued cosine and sine transforms are also available for many densities, including the Gaussian density.
B.3
Gaussian family base measures and neural network features
Integrals of the form of the squared family kernel arise as so-called neural network Gaussian process kernels, when the features are hidden layers of neural networks and the base measure is Gaussian, That is, the features take the form (w) o (WW + b) for some weight and bias parameters W and b, and some activation function and
(K (K,))) )ij = E [aww] w + b.)a(w)
by)]
Example activation functions with known closed form include the error function [Williams 1996]. the ReLU (and other rectified monomial) function [Cho and Saul 2009 and sine and cosines |Paaree et al.] 2020). Further examples are given in Table of Tsuchida et al. (023) and Table Han et al. 2022 J, including monomials, GeLU, snake and approximations to well-behaved functions.
B.4 Table of some closed-form squared families
16

---

Domain
Base measure
Feature (b(w)
Neefeeence/Notes
Rd
(0, Cov)
ef(WWw)
Williams
1996
Rd
(0,, Cov)
ReLU(WW w)
Cho and Saul
2009
Rd
(00 Cov)
GELU(Wa w)
Tsuchiia et al.
2021
gd
(mean, Cov)
Snake Ww + b)
Tsuchida et al.
2023
Rd
Ammann, Cov)
exp( Ww w + b)
Nielsen and Garcia
2024
{0, .)
(w!)-
exp( Ww b)
Nielsen and Garcia
2024
Sd-1
Uniform
<xp(w=- b)
Brown [1986]
Table 2: All of these settings admit dlosed-form squared family kernels K p.v.- In each case, the
closed-form kernel is computed linear time in the dimension d. Here erf denotes the error function,
gd-1 denotes the unit hypersphere, Snake denotes the snake activation function [Ziyin et al. 2020)
Even further examples are given in Han al. (202 table 1] and Tsuchida al. 2023 table 1].
17

---

C Proofs
Proposition Consider the model Figure 3 (Left). Define measure v(dw) p(U w)
The posterior P(dw U)= Q(dw M, 1, 4) belongs to GSF with base measure feature mapping ap and parameter M.
Proof. The posterior is the product of the likelihood and prior divided by the marginal likelihood,
Tr (MK n. P(U) (M (W)))((U w)) (M (W))))) w))(dw) x Q(dw M.
Proposition Consider the model in Figure 3 (Left). Define a measure v(dw) = p(U (سسارر and v' (dw) = p(U ())(dw), where, abusing notation, the train likelihood p(U w) and test likelihood pU,, w) are not necessarily the same. The posterior predictive is the ratio of normalising constants, p PUU U) Tr (MK, v1 // Tr (MK..).
Proof. This follows by applying Proposition Catttiioontiicee twice, directly. Directly, the posterior predictive obtained by marginalising out the predictive distribution given the parameters with respect to the parameter posterior. That is,
PUU p(U. P(dw u))
/ p(U= w)Q(dw M,
Tr (T((b(W)))(()) p(U. w))(dw)
(MK
Tr (MKu..))
TF (MKu.b)
D Squared probability process calculations
The derivation of the following corollaries may appear long and symbol-laden but is essentially just mixing the parameter-integral factorisation (see $ $(2) with the standard "completing the square"
technique Gaussian linear regression [Bishop] 2006 $2.3.1] Rasmussen and Williams 2006 $2.1.1].
Corollary 7. Consider the special model 3 and 5 applied to Figure (Right). Define vdww) = p(Y. I )p(Y IΩ ()()). The Gadon-Nikodyy derivatives of v and v/ are the same in the finite-feature Gaussian process regression setting. That is,
du (dw) = NNW m'y, C')NN, dx
where
m'y vec (B CT (B & T
(vec(Y) (Lo, @ I') vec(M))
C' = C (B CT
CT T + rC),
and
dv' (dw) = -N(vec(YY) (Ia & u IN + (La, @ T.)C @r.))
d
N(W m Nr,,,,, where
18

---

y (Id, 8 IT. & (T.)C' ST.)) + IN,da) ( - (Id @ m.)mmy)
+ +myy
CC' @ ((T,, & (T,) @ T.) IN.de) (Idg @T.)C'
Here NY hides an inconsequential closed-form factor which depends only on Y. Hence the noisy posterior predictive distribution is available as p(Y, Y.T.T.) = Tr (MKr.v) /)) Tr (M (M,)) whenever K sp can be computed in closed-form for Gaussian base measures.
Noise-free posterior predictive In general, the noisy posterior predictive distribution p( Y Y. Y.T.T. is not a GSF, However, if the posterior predictive distribution evaluated on a large enough test set of size N. 2 d1, and the resulting feature matrix P. is full-rank, then the noiseless posterior predictive distribution is itself a GSF.
Corollary 8. Suppose additionally that C' and T. are full-rank and N. 2 di Then
lim +²+
m (Idz @ r.)T ve (Ids & (T, Tr) vec(Y.)
lim C" 2
===0, 0,
and the noise free predictive distribution (x)f(X.) |Y,T.T.) Q(df(x.) | M, belongs to a GSF with modified f feature mapping and modified base measure,
4s (f(x))) - +b ((La, &x.))((x.)) = = (vece ( (T!f(X))))
v" (df((x.)) & N( vec (f(X.)) (Idz @ m'y (Ia, @ IIa @ of.)) (7)
Proof. v2.
Take the train likelihood to be multivariate Gaussian with homoskedastic noise with variance
P
rΩ = (2mo2)dida/2
exp (- 2v
Tr T( ((-- -rn-)rr r22)))
(2To2)did2/22 exp (- 20
(vec(()) (Idz @ T)w)' T (vec(Y)) (Id, Tiw)
N(( (voc(Y) (I @ T)w, v² Idd,)
Take the Radon-Nikonn Nikoddm derivative of the base measure y with respect to Lebesgue measure to be the density of matrix Gaussian distribution,
exp 1T - -M -M))
du
12 p(52) dA
(2 Bl /C
= =N(w vec(M), B C)
We convert the product p(Y
rΩ ((()) into product of Gaussian probability density
functions. In order to do so, we write p(Y Tu )p() = p())p(1 Y) = >Y, S2) and note that and Ω are jointly Gaussian. Since LΩ E. the parameters of the marginal distribution of can be simply read off as
E veol Y)] = vec(TM),
B[ vec(Y) vec -2II Nd, + E vec (
-I2INa + (Ia, & F) voc(S2 voc(S2 T (Ia, T
I + (Ld & T) 8 Cl @ @T
IIIdd + B (TCT").
so that
p(Y) -N(eve())) Ve
2 INd + & (TCT
(T, vea(M)
19

---

For p( Y), we compute the unnormalised posterior Gaussian probability density function by completing the square of the quadratic, following typical calculation. This yields
p Y) = =NWw m'y,C'r)
where
my = v + @ X & - (Id, @ vec(
Cy B C B CT (B @ Pcr
rC).
Therefore,
dp (2) = ((vee((())))))) (Is, @ F) vec( Ind, + -Bo @ (TCr*) (TCT"))N(u my, Cy),
P(Y
rQ
(8)
Applying the same result again for p(Y.
I. (voo(Y) ( & &T.)., IN.ds) we
may form the probability density function corresponding with the base measure 1/,,
du dA
pYY. r.s ) rQ
(S2)
=N(vec(YY) (Id, @ F) vec (M (M) No, + B @ (TCT)
((vec((y,) (Iz @ D. m. IN.d, + (Ia, & (Ia, & | m
where
m Y. =myy + Clss 8 (U,,, (r 8 r.. v IN.ds) - @ P.)C'y.
8 IN (vec(Y)) (Idz @ T.rmmy
Cr.Y. = Cyy Clo Tr (1Ta, @ &T.jc (L @ @r.)
LL be decomposition in terms of L such that L Rdiddx and suppose Cýy is invertible. We may use the limit definition of the pseudo-inverse to find
Limiting base measure // Let Cy
lim Cy @ ((La, @ (4a @ r.) + v IN.da)
+0+
lim LL (L 8 (Ul, &T. (La @ @r x2 I
-²-
((la, (.),))r))))
The matrix (Id, @ p.) E N.d2xdidd has linearly independent columns, assuming N. is large enough. L has linearly independent rows (since C² is invertible). Therefore
a( 8 &T)L)
(I T..
-Id, ② T
We may then obtain the base measure lim v/ (dw) which has parameters
by.-- A lim m Y. (Id. & &T.) vec(Y.)
-0
lim Ch,y. =. +==
We therefore have that
lim v'(d v'(dw) NV( vec(Y) (Ld, vec M) v2 Nd. +B +B (TCTT)) N(voc(( (Id, @ r, my (Id, @ @r.
+0+
& @ T. lx.y.).
20

---

Integration against this measure can be understood as an evaluation functional, and since ly Y- is linear in Y., a, the evaluation is also GSF. More concretely.
2()) S. [|0b(W)) .v.Y-) dw
((f(x.) (1 T.)my, (1 @ @I,)Cy ( r.)))
= 100 x- x-10xr,)))))11((( veo(f(XX)) (Idz @ T.)my,, (Ta, @ T..)Cy @ r)))
= - T.) veel vec(f((.))))
N( ((vee(f(X.)) (Tdg @ r. r.) m, (Iaz @ (Ta + ).))
Hence the posterior predictive bleongs to GSF with base measure vddf(x.))
N( vec (f(X)) | (1dz @ T.)miy (1s @ @T,) o df()).
21

---

E
Algorithms
In response to reviewer feedback, include a description of the algorithms used in our experiments. During the experiments, we run Algorithm 3 on the meta datset then run Algorithm 2 on the support/query dataset.
Algorithm 1 Marginal likelihood (i.e., Proposition 2)
Input: Trainable final layer parameter O trainable hidden network TRd -> IR", trainable base
measure 14, frozen likelihood p(U w)
Output: Marginal likelihood p(U) [0, 00)
1: Set M = θ θ.
2: Compute K p.
D see Appendix
3: Set v(dw) p(U w))u(d)) then compute K v
>> see Appendix B
4: Set (MM = Tr(MK ,.v)
D> as per
5: return (UU) Tr(MK (/2M/
> as per Proposition
Algorithm GSFP posterior predictive update (i.e., Corollary
Input:
Frozen Gaussian likelihood p(Y 7 X ( support set. Query X" Prior element of GSF parameterised by: final layer parameter O [gmx*, hidden GSF network ?b IR R" Gaussian base measure μu, deep feature extractor 7 X R..
Output: Posterior element of GSF, parameterised by: final layer parameter, hidden network, and base measure.
1: Set M = θ.
2: Set p(= x²)) be the Gaussian posterior predictive.
per as per Corollary as per
3: Set up"
4: Set xMM) Tr(M ". :
5: return Posterior element of GSF parameterised by M, 1/ and sb.. That is, P(df
M, v", ""'""
D as per ① and Proposition
(Here shorthand for f(X)))
Algorithm 3 GSFP Few-shot learning pretraining (i.e. paragraph of Section
Input: Meta-datasets Xx, Yjj for j 1,...,J J. Batch size B {1, (1,.... J). Trainable final layer parameter + R mmxx Trainable hidden network &b Rdd R". Trainable Gaussian base measure μu. Deep feature extractor 7 + Rd.
Output: Trained GSFP prior distribution
for each batch b 1..... [J/b] do
2: Set U, = (XX,)(X)) (). for each index in the batch. Compute Gaussian likelihood p(Y;
7 X (w;) as per (4) for each index in the batch
3: Compute gradients of the marginal likelihood (Algorithm I input (A, 1h, p(Y;
(X))w)) for each in the batch, with respect to parameters O and parameters of 4b, / and 7,
using automatic differentation.
4: Update O and parameters of wb, a and 7, using the computed gradients (e.g. using SGD update rule).
5: Optionally store gradient statistics (e.g. for Adam or other optimisers).
6: end for
7: return Trained GSFP prior distribution described by O wb, p and 7.
22

---

F Experiment details few-shot regression
Here describe the six datasets (two which come with licenses, the rest without), the implemen- tation details, and the run times for each method.
F1
Sines
Sines dataset |Finn al. 2017 is comprised of input from the interval [-5, with corresponding label y = A sin(x + p) + 6, where amplitude A U 1, 5.0] and phase ~ [0, π] are uniformly
distributed on their respective intervals, and (O, 0.1) is a Gaussian noise with mean and standard deviation of 0.1.
In training, 10 input-label pairs sampled with five each for support and query sets. Evaluation performed on 500 inference iterations on 200 data points, with a 5/195 split ratio for support and query sets. Testing inputs are sampled from [-5, 5] for in-distribution evaluation and [-5, 10] for out-o--disttibution evaluation.
F.2 Mixed-Noise Sines
Mixed-Noise Sines [Sendera al. 2021 is variant of sines experiment by utilising input-dependent noise: y = A sin(x + p) + + P where | an absolute value function.
F33 NDX100
NASDAQI00 Small Dataset [Qin ct al. 2017] contains stock prices of 81 major corporations and index values of NASDAQ 100 from July 26 to December 22 in 2016. The data are recorded at frequency of one point per minute, resulting in 390 data points collected in standard trade day, with exceptions of 210 and 180 points on November 25 and December 22 separately.
Training sets and in-distribution evaluation sets are obtained by partitioning NASDAQ 100 index with 70/30 split. nut-of-diitribution evaluation performed the entire data series of Yahoo stock prices. During both training and evaluation, 10 input-label pairs are sampled with a fixed frequency of I point per 10 minutes. Both in-distribution and out-of -iitt--diittibttion evaluations report average performances over 500 sampled data sequences.
E4 EEG
EEG Steady-State Visual Evoked Potential Signals Dataset [Fernandez-Fraga et al. 2019 contains the electroencephalogram (EEG) data from 29 subjects performing different visual tests. The time-series data signals of various lengths (owing different test durations) recorded a frequency of 128 Hz. This dataset available under a Creative Commons Attribution 4.0 International (CC-BY-4.0) licence.
Electrode AF4 data obtained from the first Five Box Visual Test I on Subject 1 from Group A (file name AB01SB1 l_c__vvcv is used for training and in-distribution evaluation following partition ratio of 70/30. Out-of-distributtin evaluation is performed on the entire electrode AF4 data obtained from the first Five Box Visual Test I Subject from Group A (file name A003SB1 1.csv). During both training and evaluation, 10 input-label pairs are sampled with a frequency of 12.8 Hz. Both in-distribution and out-of-distribttion evaluations report average performances over 500 sampled data sequences.
F.5
QMUL
Queen Mary University of London Multiview Face Dataset (QMUL) [Gong et al. 1996 contains normalised facial images of 48 people. There are 133 facial images associated with each person. rep- resenting slice of view sphere where yaw and pitch span {100 x 999, 8 8,9) and B E {10e x j E
2....,2.3. { (-3,- 22. ...... 2.3) respectively.
We use the subset comprising of only 37 people, whose facial images are grayscale, for training and evaluation following a split ratio of 32/5. Randomly sampled head trajectory (a sequence of yaw-pitch pairs) is obtained as: {(a, 10 |A sin(a/10 p) + 3 30)}, where amplitude A (((0.1,5.0)
23

---

and phase ~ UO0 π| are uniformly distributed. Images and tilts along the trajectory are used as input-label pairs for training/evaluation. Note in out-of-distrbuttion experiments, yaw is limited to subset (100 x i (-9, -8.- ..... -8..... }) during training, but has the entire set in the evaluation. Both in-distribution and out-of out-of-distribution evaluations report average performances over 500 sampled data sequences.
F.6
Power
Power (Hebrail and Berard, 2006] dataset contains 2,075,259 electricity consumption readings in
house in Sceaux (France), spanning 47 months between December 2006 and November 2010. The time-series measurements of electricity consumption (watt per hour) recorded at frequency of one reading per minute on individual metres. This dataset licenced under Creative Commons Attribution 4.0 International (CC-BY-4.0) licence.
Readings from Sub metering are used for training and evaluation of the models. Following the setting (enndera et al. [202]] 70% of the data, starting from the beginning of the time series, is used for both training and in-distribution evaluation, while the remaining 30% is used for out-of- distribution evaluation. During both training and evaluation, 10 input-label pairs are sampled with frequency of 1 reading per 10 minutes. Both in-distribution and -utto--distributtin evaluations report average performances over 500 sampled data sequences.
F.7 Implementation details
GSFP implemented using Pytorch and GPytorch frameworks. Features ((x) of inputs æ are obtained using Multilayer Perceptron (MLP) in Sines, Mixed-Noise-SSie -Sines NDX 100 and EEG experiments, and Convolutional Neural Network (CNN) for image data in QMUL experiment. The MLP has two fully-con cull-cconnected lavers. with the first layer mapping X + R followed by ReLU activation function [Hahnloser et al. 2000 and the second one mapping IR Rd. The CNN has three convolutional layers with fixed kernel size of 3, stride of 2 and dilation of 2. Input channel number of the first convolutional layer 3, and upsized to 36 the two subsequent convolutional layers. Output channel number is 36 across all convolutional layers. ReI.U activation function is appended to each of the first two convolutional layers. Output of the last convolutional layer is flattened into a d-dimensional vector. scaled cosine function cos( (-)())) (-)/vd is appended as the final activation.
The squared neural network consists of hidden layer (1b f E where R and readout parameters ( E R.**). The hidden layer is fully connected layer (w(w) (Ww + b) with a weight matrix W E Ra initialised from a standard normal distribution (O I). and bias vector b E IR" initialised be all ones. For the activation function o,, we use the Snake activation function Ziyin et al., 2020].
1
1
1
Snake, (2) = z sin sin²(az) 2
cos(2az) +
a
2a
2az
The dimensions of data, latent feature and output used across the experiments are summarised in Table 3
Table Summary dimensions used across the experiments. Recall that is the dimension of the feature y in the regression model [=x w ((x) (see (2)), and and m are the dimensions of the readout parameter O mmxn in the GSF model (see $
Datasets
Dimensions
Sines
Mixed-Noise Sines
NDX 100
EFG
Power
QMUT.
d
40
40
2,916
10
10
500
m 10
10
500
Both DKT (Patacchiola et al.) 2020 and NGGP [Sendera et al. 2021 are integrated into our code to facilitate their evaluations under the same experimental settings as those of GSFP. Their implementations follow the details specified in their original works, with only exceptions of increasing the dimension the range the feature extractor or feature mapping function (r) from to the 1DX100, EEG, and Power datasets.
24

---

F.8 Training details
All models, including DKT Patacchiola et al. 2020] NGGP [Sendera et al., 2021 and GFSP, are trained with Adam optimiser [Kingma 2014 using fixed learning rate of 0.001 and default beta coefficients of 31 0.9 and 52 0.999 across all experiments. All models are trained on J = 10, 000 training datasets that are randomly sampled from the meta dataset as detailed in & F. F.6 The evaluation results of each meta model is averaged over 500 test datasets, with each entry representing the mean + the standard deviation. Optimal performance of each meta model is obtained with grid search at frequency of 500 training datasets. Both training and evaluation are conducted with single RTX 2080TI GPU and an Intel(R) Xeon(R) Gold 6242 CPU @ 2.80GHz.
F.9 Run times
Table 4: Benchmark results showing training time (mins) and inference speed (tests / meta dataset across regression tasks. The results are evaluated on meta models, and each entry shows the mean standard deviation over these meta models.
Methods
Training Time (mins)
Name
Kernel
Sines
Mixed-Noise Sines
NDX100
EEG
QMUL
Power
RBF
1.04 + 0.02
1.02 + 0.05
1.05 + 0.02
0.98 + 0.01
3.06 0.11
1.01 + 0.05
Spectral
NN Lincar
1.02 + 0.01
0.98 = 0.03
0.97 = 0.01
0.95
0.01
3.12 0.03
0.98 + 0.01
RBF
41.59 +: 0.60
43.57 = 0.58
53.41 +: 2.39
55.97
+ 5.54
57.68
0.88
75.30 2.89
NGGP
Spectral
37.40 1.41
35.19 = 0.66
53.57 + 0.87
61.33
7.22
57.13 + 1.41
70.52 3.48
NN
42.53
5.21
0.73 65.93
NN
3.45
Methods
Inference Speed (tests f
Name
Kernel
Sines
Mixed-Noise Sines
NDX
EEG
QMUL
Power
RBF
160.73
737331 de 1.17
167.84 84 4: 1.98
170.96
1.28
174.98
: 1.88
59.92 d: 4.26
169.38
d: 0.44
DKT
Spectaal
166.22
0.52
163.67 + 0.64
172.98
0.75
175.61
+ 1.44
60.06
1.46
174.97
NN Linear
184.35 d: 1.46
185.00 1: 1.31
195.49
0.68
196.61
1: 0.59
62:25 4: 1.97
195.17
1.12
RBF
5.49 4: 0.07
5.90 0.05
4.33 4: 0.87
4.11 0.46
3.84 0.24
3.85
0.28
NGGP
Special
7.00 4: 0.35
6.80 0.85
4.40 +: 0.58
4.33 : 0.13
3.87 1: 0.09
3.78
0.16
NN Lineur
0.94
GSFP
NN Linear
99/.08 0.97
150.09 1.86
149.61 :
7.14
149.41
25

---

G Experimental details vanilla regression
Here we describe the nine datasets selected by [Hernández-Lobato and Adams 2015] for regression experiments, as well the implementation details, experimental performances in terms of NLL, and run times for each method.
G.1 Boston Housing
The Boston Housing dataset (Harrison Jr and Rubinield 1978 consists of 506 instances with 14 variables describing the socioeconomic characteristics of neighbourhoods, in Boston, Massachusetts. We use 13 features as imput variables and the remaining variable, MEDV (median value of owner- occupied homes in $1000's), serves as the regression target. This dataset is licensed under an Apache open source license.
G.2 Concrete Compressive Strength
The Concrete Compressive Strength dataset Yeh 1998
consists of 1,030 instances with variables describing the compositional structure, the curation duration, and the compressive strength of concrete samples. We use 7 features input variables and the remaining variable, concrete compressive strength (MPa), serves as the regression target. This dataset licensed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license.
G.3 Energy Efficiency
The Energy Efficiency dataset [Tsanas and Xifara 2012 consists of 768 samples with 8 features describing the design details the buildings. We use these features as imput variables and the variable, cooling load, serves as the regression target. This dataset is licensed under Creative Commons Attribution 4.0 International (CC BY 4.0) license.
G.4 Kin8nm
The Kin8nm dataset part of the Kin family datasets [Corke 2002 that consists of 8,192 instances with 8 features simulating the forward kinetics of robotic arm. We use these 8 features as input variables and the remaining variable (y) serves the regression target.
G.5 Condition Based Maintenance of Naval Propulsion Plants
The Condition Based Maintenance of Naval Propulsion Plant dataset [Coraddu et al. [014]] consists of 11,934 samples generated from a numerical simulation naval vessel with a gas turbine propulsion plant. We use the 16 steady-state sensor and operational measurements input variables, and the compressor degradation serves as the regression target. This dataset is licensed under Creative Commons Attribution 4.0 International (CC BY 4.0) license.
G.6 Combined Cycle Power Plant
The Combined Cycle Power Plant dataset Tuückki 2014
consists of 9,568 data points collected from a fully operational power plant over a period of 2006-2011. I. We use the hourly average readings from the four sensors input variables and the remaining variable, net hourly electrical energy output (EP), serves as the regression target. This dataset is licensed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license.
G.7 Physicochemical Properties Protein Tertiary Structure
The Physicochemical Properties of Protein Tertiary Structure dataset [Rna 2013
consists of 45,730 samples of decoy structures. Each sample has nine physicochemical descriptor features, which we use as input variables. The remaining variable, root-mean-squee deviation (RMSD), is the regression target. This dataset licensed under Creative Commons Attribution 4.0 International (CC BY 4.0) license. We use the first 15,000 samples for regression experiment.
26

---

G.8 Wine Quality
The Wine Quality dataset Cortezz al. 2009] has two sub-datasets regarding red wine and white wine from the Portuguese "Vinho Verde region. Following Sermandez-Lobatee and Adams 2013]. we choose the red wine sub-dataset for our regression experiment. The red wine sub-dataset has 1,599 samples, each of which is described with 11 physicochemical features. 'hee quality variable, taking ordinal integers ranging from 10, serves as the regression target. This dataset is licensed under Creative Commons Attribution 4.0 International (CC BY 4.0) license.
G.9 Yacht Hydrodynamics
The Yacht Hydrodynamics dataset [Gerritsma et al. 1981 consists of 308 samples with features derived from hull geometry and hydrodynamic context. The output variable, resistance, serves as the regression target. This dataset licensed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license.
G.10 Implementation details
Most of the implementation details for GSFP, DKT, and NGGP in the regression experiments follow those in Appendix F.7 with differences as follows. (a) The feature extractor is the two-layer MLP; (b)
2. and We originally considered the same few-shot learning benchmark (all datasets) as in Sendera et al. [202T] During the rebuttal period we applied DKT, NGGP and GSFP vanilla regression datasets Boston, Concrete, Energy. Kin8nm, Naval, Power Plant, Protein, Wine and Yacht. In applying our multimodal, non-Gaussian model vanilla regression, we repeat the message Sendera et al. [202T]: "the main goal of [NNGP [Sendera et al. 2021]] was to show improvement of NGGP over standard GPs in the case of few-shot regression Intuii.n that NGGP may be superior to standard GPs in simple regression setting for datasets with non-Gaussian characteristics, but do not expect any improvement otherwise". We stress that we did not perform hyperparameter tuning during the rebuttal period, and the results which follow are merely to demonstrate the performance of the models previously tuned for few-shot regression in vanilla regression settings. As such, we caution against general takeaways.
We fix d =
Both and (b) are applied across the nine regression experiments.
G.11 Training details
All models, including DKT PPatccchoola et al., 2020), NGGP [Sendera al. a [021] and GFSP, are trained with Adam optimiser [ using fixed learning rate of 0.001 and default beta coefficients of 31 0.9 and 32 = 0.999 across all experiments. The train/test split ratio of 80/20 is applied across the datasets described in & G. on each dataset using random seeds ranging from to All models are trained on the train set for 10,000 iterations. The evaluation result of each method is averaged over 5 models, with each entry representing the mean + standard deviation. Optimal performance of each meta model is obtained with grid search at frequency of 100 training iterations. Both training and evaluation are conducted with single RTX 2080TI GPU and an Intel(R) Xeon(R) Gold 6242 CPU 2.80GHz.
G.9 Five models are trained for each method
G.12 Results
As shown in Table 5) we observe that GSFP obtains the lower NLL on 7 out of 9 datasets, and the second lowest NL. on the Protein dataset. For Kin8nm, Naval, Power Plant and Protein datasets, DKT and NGGP with spectral kernels cannot be trained due to the large number of trainable parameters associated with the spectral kernel, which cannot in our available GPU memory.
G.13 Run times
27

---

Table 5: Benchmark results showing test NLLs for regression tasks across datasets deseribed in $ G.G.G G.9 For each dataset, we train 5 models with different random seeds (1-5) via maximum marginal likelihood over the train At testing time, condition the train set and evaluate the NLL on the test set for each model. Each entry shows the mean + standard deviation over model evaluations. Best results in bold, second best underlined.
Methods
Patacchiola
Sendera
RBP
Spectral
NN Lineur
RHP
Spectral
NN Linear
NN Linear
Baston
0.333
0.198
0.491 -0 0.179
0.478 0. 01122
0.271 + 0 0.220
0.439 0.228
0.288 + 0.084
0.345
0.272
Concrete
0.300
0.125
0.302 0.135
0.442 0.082
0.847 0.255
0.447 0.247
0.411 0.407
0.072
0.175
Energy
00.92
0.175
0.277 0.084
-0474
0.118
0.481
0.507 0.027
0.364
1.:44
0.393
Kinkkm
0.153
0.011
0.388 + 0.1055
0.310
0.150
0.655 + 0.113
0.021
0.080
Naval
1,070 0.011
0.994 # 0.156
0.405 +: 0.474
0,830 =: 0.135
22.019 :: 0.219
Power
0.022
03:336
0.137
1.144
Wine
1.162 0.034
1,192 0.050
1.170
0.001
-1.881 264
1.444 0.226
1.696 =: 0.664
01735 2: 0.133
Yacht
-1850 0.251
-12778 1.219
-1.496 0.256
-2,350 à: 0.534
0,542 0.9
0.106 0.591
-2.873 1: 0.238
Table 6: Benchmark results showing training time (hours) and inference speed (tests / s) for each model across regression tasks. The results evaluated models, and each entry shows the mean + standard deviation over these 5 models.
Training Time (hours)
Methods
Datasets
DKT
Patacchiola
2020
NGGP
Sendern
2021
GSFP
RBF
Spectral
NN Linear
RBF
Spectral
NN Linear
NN Linear
Boston
0.034
0.001
0.037
0.001
0.032 0.002
5.042
0.619
4.024
0.198
4.828 + 0.353
0.100
0.003
Concrete
0.001
0.007
0.112
0.005
0.115 0.011
10.000
1.209
10.592
1.400
9:229
1,258
0.176
0.003
Energy
0.048 + 0.000
0.053 0.001
0.045 + 0.001
7.585 + 0.705
7.741 0.645
6.886 + 0.450
0.121
0.004
KinSum
62
0.002
Naval
3.149 0.0+000.0000 0.075
1.367 + 0.052
127 623 + 11.349
108.597 + 5.025
2,025 + 0.008
Power Plant
1.611 0.271
0.539 0.061
76 76.962 + 4.722
68.406 + 5.127
1.519 0. 0.008
Protein
0.509
2.303
144.444
111.314
12.836
2,808 0.020
Wine
0.099 = 0.008
0.175 = 0.017
0.101 2 0.008
16. 680 + 2.621
18.150
3.6 623
18.791 + 4.215
0,206 = 01.0000000 0.00L
Yacht
0.000
0.035
0.032
0.108
4.432
0.206
0.001
(seconda)
Methods
Patacchiola
2020
Seadera
GSFP
RBF
Spectral
Lineur
RBF
Spectral
NN Linear
NN Linear
0.30
0.30 0.00
0.27
7.36 0.07
0:01
Concrete
1.10 1100::00 0.00
1.43 0.05
1.55
0.17
14.11 1 2.33
12.95
1.35
15.74 1: 1.49
1.04
0.00
Energy
0.70
0.79 0.01
0.55
01.00
8.20
0.82
8.38 0.93
8.61 + 0.03
1.02
0.02
Kinšum
Naval
67. 4:261
12.08 - 3.60
178.25 + 3.59
155.69 + 2.97
15.03
0.119
Power Plant
42.51 0.56
27.32 2.16
144.27
4.71
134.66
10.57
9.67
0.02
Protein
118.11 0.12
70.83
2.72
297.26 1 14.70
211.55 4: 4.89
15.10
0.03
Wine
1.70 0.00
2.98 0.10
1.62 1: 0.00
17.08 1 0.72
25.
1.95
21.76 4 5.18
2.10
0.01
28

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: All of the claims in the abstract and introduction are proved theoretically or demonstrated empirically.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: A limitations section is provided in the last section of the paper, discussing computational complexity and lack formal result for maximum marginal likelihoods.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g.. independence assumptions, noiseless settings, model sell-speciiiatiinn asymptotic approximations only holding locally). ). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
29

---

Answer: [Yes]
Justification; All the results in this paper are derived either in the main paper or in the Appendix, as mentioned in the main paper. All are numbered and cross-referenced. The proof of the first Proposition 2 serves as a short proof sketch for the remainder of the results.
Guidelines:
The answer NA means that the paper does not include theoretical results, All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: The broad experimental setup is described in The details of the datasets and all hyperparameters are provided Appendix F
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
30

---

5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: Code is provided in the submission, and will be open access after acceptance. Instructions as well as exact commands are provided to reproduce experimental results.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.co
www.ibimmssiooooooiiiioc.. for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
7 //nips. wwtpccppppppcppu..ccccc gttp://wiiess/ee/eee www.eoewwisssoonnooooooooooo
for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification: All of the training and test split details are described in Appendix F Details on how hyperparameters were chosen is provided in Appendix
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 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, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: Table 4 shows performance metries with mean and standard deviation 500 random shuffles of the test meta dataset and 5 model training from different random seeds (i.e. in total, 2, 500 models are trained and tested in a few-shot setting, per entry the table).
Guidelines:
The answer NA means that the paper does not include experiments.
31

---

The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the enror 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.
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 96% CI, if the hypothesis of Normality of errors not verified.
For asymmetric distributions, the authors should be careful to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Yes, as mentioned in Appendix both training and evaluation are conducted with single RTX 2080TI GPU and an Intel(R) Xeon(R) Gold 6242 CPU 2.80GHz. The compute time is given in Table F.9 No additional compute was required than the experiments reported in the paper. 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 compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (c.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the http:///wrrpppcccpppppppe/elttta
NeurlPS Code of Ethics
Answer: [Yes]
Justification: Code of Ethics adhered to (including considerations for potential harms caused by the research process, societal impact and potential harmful consequences, impact mitigation measures) and there are no special circumstances which apply.
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 deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative socictal impacts of the work performed?
32

---

Answer: [NA]
Justification: Our work is similar to the example pointed out in the guidelines below: "it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster." Our model serves as generic probabilistic inference machine, and in principle could be used enable harmful applications more quickly. However, this applies generically to all research in technology.
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 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 partieular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate 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 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 attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Similar to above, our contributions are sufficiently well-separated from direct applications with high risk for misuse (e.g. language models, image/video generators, problematic datasets,...)
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have 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 best faith effort.
12. Licenses for existing assets
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?
Question:
Answer: [Yes]
33

---

All datasets and competitor models are cited in $ Licenses are discussed in Appendix F and are respected.
Justification:
Guidelines:
The answer NA means that the paper does 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 used and, possible, include 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.
assets are released, the license, copyright information, and terms use in package should be provided. For popular datasets, paperssii thoode com/ com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re- e-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: Yes, well-documented README.md contains instructions to reproduce the results and use the provided software. This will be provided in the full release. license is provided.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/mmoel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethies, workers involved data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
Institutional review board (IRB) approvals or equivalent for research with human subjects
34

---

Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: As above.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: The core method development in this research does involve LLMs any important, original, or non-standard components (in fact, any whatsoever).
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non non-standard components.
Please refer to our LLM policy //neur cc. for what should or should not be described.
35


Reviews:
[{"review_id": "dCvw6FSRZa", "reviewer": "Reviewer_Z38J", "summary": "This work considers the use of Generalised Squared Families (GSF) as **priors** for Bayesian linear feature regression problems.\n\nGeneralised squared densities $q(\\omega | \\theta, \\mu, \\psi)$ are defined with respect to measure $\\mu$, proportionally to the squared Euclidean norm of the inner product $\\theta \\psi(\\omega)$. Due to its analytical properties, the density can be written in terms of a ratio between two analytically tractable traces ---as in Equation (1) ---, generalizing mixture models.\n\nThe authors demonstrate in Section 3 that closed-form expressions for the marginal likelihood, posterior parameter and posterior predictive distributions are attainable for GSFs, whenever the base measure of the GSF is itself conjugate to the likelihood of the model.\n\nBased on these analytical properties, the authors tackle Bayesian (linear in features) regression models, and showcase that for Gaussian likelihoods, Bayesian GSF priors over these models (called GSF processes) result in closed form posterior and marginals of interest.\n\nExperiments showcase the benefits of considering a multi-modal prior/posterior family in few-shot learning examples, where better model-fit is achieved when compared to reasonable baselines.", "questions": "- Can the authors elaborate and more clearly explain the computational complexity of GSFPs? What are the key bottlenecks and how can one implement or overcome them?\n\n- I would encourage the authors to discuss (if any) the potential of GSF priors in classical GP regression problems\n    - what are the benefits/drawbacks of this family of priors?\n    - can the authors showcase empirically how these priors would impact regression performance?", "limitations": "Yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 4, "strengths_and_weaknesses": "**The main strength** and significant contribution of this submission is to extend Gaussian Process regression with multimodal priors and posteriors via GSFs, and to show that these have closed form solutions.\n    \n- Namely, this work can be seen as a generalization of Gaussian prior with Gaussian likelihood regression models in feature space, via GSF priors, that allow for multimodal expressions of uncertainty.\n\nAdditionally,\n- The background section is a thorough and informative summary of the concepts needed to understand this work\n- Section 3, where the Bayesian linear regression in feature space with GSF priors is introduced, is a math-heavy section that clearly presents all the key assumptions and properties needed to derive the closed form expressions for GSFPs.\n\n**The main weakness**, as acknowledged by the authors themselves, is the computational (and memory) cost associated with training a GSFP.", "was_revised": true, "final_justification": "I recommend acceptance of this work based on my original assessment and the informative rebuttal provided by the authors.\n\nThe rebuttal addresses most of our concerns and provides new illustrative experiments to support their claims.\n\nIn addition, the authors have clarified certain aspects of the work and suggested edits to the manuscript that will improve its presentation.", "rebuttal": "Thanks very much for your thoughtful review. See below for our responses to your comments and suggestions. If anything is unclear or needs further clarification, please let us know.\n\n---\n\n## 1. Computational complexity\n\n(Same as response to Reviewer sKy9).\nWe should clarify the distinction between (more general) GSFs applied to Bayesian inference and (more specific) GSFPs  applied to Bayesian inference in regression models.\n\nFirst, in the general case, Bayesian inference is (in)famously computationally difficult. With arbitrary likelihoods and priors, tuning computational algorithms can be practically difficult, and even fail to give good approximations unless some heavy restriction is placed on the densities (e.g. log concave). This is especially difficult when the dimension of the random variable $d$ is large, due to the curse of dimensionality, due to the exponential growth of volume. In practice, if one desires an exact computation, we need to restrict our class of probability distributions. When we restrict the class to squared families, we are able to compute exact posteriors, a thing which is usually associated with very classical models, not neural network based models. In GSFs, the complexity of this calculation is governed by the inner product of the parameter matrix $\\mathbf{M}$ and the kernel matrix $K_{\\mu, \\mathbf{\\psi}}$ (as in equation 1). Assuming the kernel matrix is known, this inner product has complexity $\\mathcal{O}(n^2)$, where $n$ is the number of parameters in the GSF. Crucially, dependence of the complexity on $d$ only appears through the calculation of the kernel matrix $K_{\\mu, \\mathbf{\\psi}}$. In practice, the kernel matrices $K_{\\mu, \\mathbf{\\psi}}$ can be computed exactly in *linear* time in $d$. (All examples in our experiments admit linear complexity in $d$. We will enumerate more explicitly such examples in the main text --- see response to point 4 of Reviewer 7C3F). This leads to a combined complexity of $\\mathcal{O}(d n^2)$, i.e. exact inference linear in dimension $d$, which is a huge improvement over inexact inference exponential in dimension $d$, which one would typically expect from numerical algorithms in general classes of probability distributions. Of course, this improvement comes with the restriction of the class of densities to GSFs, however this restriction is not too severe because GSFs are rich density approximators (as demonstrated empirically in the experiments section, as well as in previous theoretical and empirical investigations in using special cases of them as likelihood models (Rudi and Ciliberto 2021, Flaxman 2017, Walder and Bishop 2017, Loconte et al. 2023a,b, Loconte and Vergari, 2025, Loconte et al., 2024, Tsuchida et al., 2023, Tsuchida et al., 2025)). More precisely navigating the trade-off between universal approximation (making $n$ large) and computational efficiency (making $n$ small) is an important direction for future work, however we believe having $n$ as a tuneable way of navigating this trade-off is a useful perspective.\n\nSecond, in the more specific case of Bayesian inference in regression models, we are only aware of two tractable (polynomial time) methods which allow for exact inference. These are GPs (Gaussian prior with Gaussian likelihood) and NGGPs (combining a GP and a normalising flow). We would also claim that these have desirable complexities when placed in the space of all possible function inference schemes. E.g. doing inference with an exponential prior on parameter weights and a Gaussian likelihood would be computationally very difficult in practice, becuase one would require numerical methods. DKT and NGGP, as well as ours, can indeed be understood as kernel methods. The complexities are similar (although NGGP has some further complexities due to the added normalising flow component). Practically, our method is slightly slower than DKT (i.e. GPs), because of the added GSF, but faster than NGGPs. See Appendix E.9. This is why we chose to highlight it as a limitation. Our method is faster than NGGPs, the only other neural method for tractable inference in function space which we are aware of.\n\nWe will add this discussion to the main paper, with the extra space afforded by extra one page as well as your suggestions  and those of other reviewers. \n\n---\n\n## 2. Classical regression problems \n**What are the benefits/drawbacks of this family of priors (in a classical regression setting)?**\nThe second point above is pertinent. Relative to inference in general function spaces, our method is tractable, with complexity matching GPs (DKT) and NGGP. Practically, our training and inference times sit between these two methods. Empirically, due to the added flexibility of our multimodal prior/posterior, we are able to outperform GPs. Unlike NGGPs, we are also not constrained to invertible transformations of GPs. The same applies to the few-shot learning setting studied in our original submission.\n\nIn terms of an intuitive justification of expressivity and therefore predictive performance in the context of **vanilla** regression, we offer the following hypothesis. We expect in general that the performance follows the ordering from worst to best\n\nGP (classical kernel) $\\leq$ GP (deep kernel learning kernel) $\\approx$ NGGP $\\approx$ GSFP. \n\n(and recall that GP (deep kernel learning kernel) is called DKT when applied to few-shot learning rather than classical regression).\n\nThe reason for this ordering is that:\n\n1. the seminal work of DKL (Wilson et al., 2016) shows empirically that DKL outperforms classical kernels in classical regression. This leverages the expressive power of deep learning to learn more relevant priors for the task at hand than classical covariance functions. In order for a DKL model to represent a classical GP, the deep neural network representing the feature mapping just has to be the identity function. But more generally, DKL can represent non-trivial covariance functions.\n\n2. as discussed in Appendix C of Sendera et al., the \"main goal of [the original NGGP work] was to show improvement of NGGP over standard GPs in the case of a **few-shot regression task**...Intuition is that NGGP may be superior to standard GPs in a simple regression setting for datasets with non-Gaussian characteristics, but do not expect any improvement otherwise\". The empirical results the authors show agreement with this intuition, showing mixed results on real data, but vastly improved results on a synthetic non-Gaussian and multimodal example. In order for the NGGP model to represent a DKT or classical GP model, the normalising flow model just has to be the identity function. But this is only valuable if the data benefits from non-Gaussian modelling. In few-shot learning, where multimodal beliefs are required (depending on the support set), this can be advantageous.\n\n3. NGGPs push GPs through a normalising flow in order to build expressivity. The practical representational limitations of normalising flows are well-known (Stimper et al., 2022). In particular, without special tricks which do not easily apply to NGGPs, normalising flows are only capable of representing *invertible transformations* of the original base distribution (which is in this case, Gaussian) (Papamakarios et al., 2021). Therefore, if the target function of interest is not invertible, it cannot be represented by an NGGP. In contrast, squared families have no such restriction (see e.g. Tsuchida et al. (2023) for a discussion in the context of frequentist likelihoods or Rudi and Ciliberto (2021) for universal approximation but without reference to invertible transformations), and can represent distributions that normalising flows cannot. This is only valuable if the data benefits from non-Gaussian modelling. \n\n**Can the authors showcase empirically how these priors would impact regression performance (in classical regression)?**\n\nDuring the rebuttal period, we applied DKL, NGGP and GSFP to vanilla regression UCI datasets Boston, Concrete, Energy, Kin8nm, Naval, Power Plant, Protein, Wine and Yacht. Please see response to Reviewer sKy9 for results and discussion. These results will appear in the updated submission.\n\n\n---\n\n\n## Additional references\n\nStimper, V., Schölkopf, B. &amp; Miguel Hernandez-Lobato, J.. (2022).  Resampling Base Distributions of Normalizing Flows. AISTATS\n\nPapamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S., & Lakshminarayanan, B. (2021). Normalizing flows for probabilistic modeling and inference. JMLR"}, {"review_id": "4Ta7C6Ku3a", "reviewer": "Reviewer_2rpL", "summary": "This paper introduces the use of squared families as priors. The authors show that when the base measure and the likelihood are conjugate, the posterior, marginal likelihood, and predictive distribution all admit closed-form expressions. The squared prior is capable of modeling multi-modal distributions, and its parameters can be parameterized by neural networks. This ability to capture multi-modality gives the method an advantage over Gaussian processes, as demonstrated in few-shot learning regression experiments.", "questions": "- In GSFP, the base measure is chosen to be a Gaussian with mean $\\text{vec}(\\mathbf{M})$, sharing the parameter $\\mathbf{M}$ with the squared prior. Is this specific parameterization necessary for achieving good performance on regression tasks?\n- For regression tasks, is it possible to adopt other base measures, feature maps, or likelihoods beyond the Gaussian base—for example, those mentioned in Appendix B? How robust is the method to different choices of these components?", "limitations": "See weaknesses and questions above.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 4, "strengths_and_weaknesses": "Strengths\n- The use of squared families as priors is quite novel. The conditions under which closed-form computations are possible are fairly general, making the method potentially applicable to a wide range of modern probabilistic machine learning problems.\n- The paper is well written and easy to follow.\n\nWeaknesses\n- The paper could benefit from a more detailed literature review of recent related work on squared families, including a clearer discussion of their respective settings and how they differ from the current approach. At present, this comparison appears primarily in the early part of the introduction and is somewhat general, which may make it difficult for readers to fully appreciate the distinctions and significance of this work.", "was_revised": true, "final_justification": "This is a good paper, I think it can find application in many problems of probabilistic machine learning. My concerns regarding a more detailed literature review are also addressed in the rebuttal. So, I fully acknowledge the significance of this paper and strongly supports its acceptance.", "rebuttal": "Thanks for your review of our paper on squared family priors, and for your helpful comments for improving the paper. See our response below. Please let us know if you have any additional questions or believe something needs further clarification or elaboration.\n\n---\n\n## 1. Discussion and placement of related work.\n\nWe will add an extended discussion of related works at the end of section 3, as follows. This builds upon the existing discussion in the early part of the introduction, but benefits from being after the technical part of the paper, so more technical details can be explained to the reader. It also further clarifies the distinction between priors/posteriors and likelihoods. Please let us know if this meets your expectation.\n\n\"\"\"\n\nHere we place our contributions with respect to existing work. We reiterate that to the best of our knowledge, we are the first to consider (even special cases of) squared families as priors and posteriors. The precise nature of the results concerning likelihoods is fundamentally different to that of priors and posteriors, because for example, one does not attempt to estimate parameters in forming a posterior, one merely attempts to form a posterior belief. Nevertheless, existing work on special cases of squared family likelihoods is relevant in that expressive power, flexibility and tractability of densities is intuitively relevant to both likelihoods and posteriors, and uncertainty quantification more generally. \n\n**Kernel methods** Building off the general results of modelling non-negative functions using squared elements of RKHS [Marteau-Ferey et al., 2020], Rudi and Ciliberto [2021] model probability distributions, with tractable marginalisation and closure under multiplication. They fit likelihood models within this class by minimising a regularised L2 distance between the target likelihood and the model. While these models are nonparametric in construction, estimates resulting from minimising the objective satisfy a representer theorem type instantiation. The authors show that the model admits both favourable representational capability (in being able to represent a $\\beta$-times differentiable target density with error less than $\\epsilon$ with number of parameters scaling like $\\mathcal{O}(\\epsilon^{-r/\\beta} (\\log(1/\\epsilon))^{r/2})$, where $r$ is the dimension of the domain of the distribution), as well as statistical estimation properties (as likelihoods, the proposed estimate converges at a rate of $N^{-\\frac{\\beta}{2\\beta + r}} (\\log N)^{r/2}$, where $N$ is the number of datapoints).\n\n**Probabilistic circuits** Probabilistic circuits [Choi et al., 2020] provide a computational framework for tractable probabilistic modelling. They are constructed as graphs of connected units, and by constraining the graph, allow for tractable marginalisation. However, they need to impose some structure on the units or the connections in order to ensure nonnegativity. Squared probabilistic circuits [Sladek, 2023, Loconte et al., 2023a,b, Loconte and Vergari, 2025] bypass this constraint by squaring the output of the circuit. Sums of squared circuits can also be instantiated [Loconte et al., 2024]. The focus on such works is typically in showing tractable (polynomial or indeed quadratic time/memory in the number of units) computation of normalising constants or marginalisation, as well as the inclusions of the function spaces imposed by different classes of probabilistic circuits. \n\n**Neural networks and finite feature models** Previous work [Tsuchida et al., 2025] has modelled likelihoods as proportional to the squared Euclidean norm of a function of the form $\\Theta \\psi(x)$, where $\\Theta$ is a matrix and $\\psi$ is a vector-valued function. Provided the features $\\psi$ are rich enough, such models admit a universal approximation property, typically approximating the squared L2 distance between a square root of a target density and the model at a rate of $\\mathcal{O}(1/n)$, where $n$ is a parameter count. As likelihoods, they learn target densities with a KL divergence to the target density decreasing at a rate of $\\mathcal{O}(1/\\sqrt{N})$, where $N$ is the number of samples. Such families of models have tractable normalising constants, Fisher information and statistical divergences. Furthermore, in some special cases, such families have not only tractable but indeed closed-form normalising constants, Fisher information and statistical divergences. An example studied previously was squared neural families [Tsuchida et al., 2023], where $\\psi$ is a single hidden layer neural network. (see Appendix B and new Table in response to Reviewer 7C3F).\n\n**Applications** Poisson point processes (PPPs) are helpful variations on density modelling, where one models an intensity instead of a density. Whereas in density modelling, realisations are i.i.d. draws from a target density, in intensity modelling, realisations are conditionally i.i.d. draws from a target density given the number of realisations. The number of realisations is itself a random variable following a Poisson distribution with an expected number of points equal to the integral of the intensity over the domain.  Flaxman et al. (2017) used squared elements of RKHSs to model intensity functions, whereas Walder and Bishop (2017) used a squared Gaussian process prior to model intensity functions. Both are similar, the latter being a Bayesian extension of the former. We note that placing a prior over a function which when squared gives a intensity/density (likelihood) is completely disjoint to the problem we are solving here, which is to use a squared function as a prior density. Probabilistic circuits have been applied in converting knowledge graph embeddings into generative models [Loconte et al. 2023a]. Squared neural family likelihoods have been applied to label distribution learning [Zhang et al., 2025], which is a kind of variation on classical regression where instead of the target label being a single class, the target variable is a discrete distribution (i.e. an element of the simplex). Hence, here squared neural family models model a distribution over distributions with a closed-form normalising constant, expectation, variance, and covariance conditioned on input samples. Probabilistic modelling is leveraged to provide confidence intervals, conformal predictions, active learning, and model ensembling. \n\n\"\"\"\n\n---\n\n## 2.  Notation overload of $\\mathbf{M}$.\n\nThanks very much for identifying this confusing typo! This was a notational overloading which we did not pick up on when preparing the draft. We do not need to tie the Gaussian base mean parameter to the parameter $\\mathbf{M}$ in the squared prior. We did not tie them when considering regression tasks. We will change the symbol for the Gaussian mean parameter in the revision.\n\n\n---\n\n## 3. Non-Gaussian base measures and likelihoods. \n\nWe will clarify in the updated manuscript.  Our claim is that the GSF priors/posteriors are rich/expressive, i.e. they can model non-Gaussian and multimodal distributions. This is not a statement about likelihoods -- indeed, we do not claim that GSFs are conjugate to all likelihoods. In the case of regression, we rely on a Gaussian likelihood and our results follow from its conjugacy with the Gaussian base measure -- but note that the priors and posteriors have a complex neural-network based non-Gaussian density with respect to that base measure (this is what makes them expressive while preserving conjugacy). We also analyse other likelihoods (e.g. Poisson in Figure 1, others in Appendix B).\n\n(Note that we did not use Gaussian feature maps $\\gamma$ or $\\psi$, and the feature maps do not need to be Gaussian. See table in response to 7C3F. For GSFP regression, only base measure $\\mu$ and likelihood are Gaussian. More generally for GSFs, non-Gaussian base measure $\\mu$ and likelihoods can be used).\n\n---\n\n## Additional references\n\nZhang, D. Label Distribution Learning using the Squared Neural Family on the Probability Simplex. In The 41st Conference on Uncertainty in Artificial Intelligence. 2025"}, {"review_id": "mBsgPuVnws", "reviewer": "Reviewer_sKy9", "summary": "This paper explores *Squared Families* as prior distributions and the conjugacy of prior distributions belonging to *Generalized Squared Families*. GSFs are indexed by a parameter, $\\textbf{M}$, a base measure, $\\mu$, and a feature, $\\Phi$. The authors show GSF priors can admit closed form expressions for commonly used objects (e.g. marginal likelihoods, posterior predictive distributions) when, for example, the base measure is conjugate to the likelihood. The paper explores particular settings (base measure + likelihood pairs) and works through results for those settings. The paper concludes with an experiments section comparing against two baselines on a collection of datasets.", "questions": "1. The computational benefits of GSPs were presented as one of their sources of value, but section 5's limitations makes me question if this is actually true since the complexity presented on lines 269-271 seem comparable to naive GPs. Can the authors please explain what in particular makes the reported complexities desirable? These seem as restrictive as any kernel based method. \n\n2. The paper presents a substantial amount of material, but it is not as accessible as it could be for the work to have broader impact. For instance, Section 3 could be shortened, as the key point about conjugacy for exponential families is already made immediately after Proposition 4. The content in Sections 3.2 and 3.3 could be streamlined, with technical details moved to the appendix. Currently, both regression and generalized squared family processes are presented with equal detail, but shifting the finer points to the appendix and focusing more on the differences between the two settings would make the exposition more effective\n\n3. In the experiments, the methods are evaluated solely on predictive performance, making it difficult to assess whether the multimodal capabilities of squared families are actually leveraged. Aside from Figure 4, do the authors have additional results that demonstrate the claimed benefits? The paper does not convincingly show that the unimodality of Gaussian processes is a practical limitation or that the proposed method meaningfully addresses it. A clearer argument supporting this claim would help strengthen the paper.\n\n\n4. The experiments are limited in both scope, given the narrow range of datasets and baselines used, and in depth, as only NLL is reported. Additionally, it is difficult to determine whether the observed differences on these small-scale datasets reflect meaningful methodological differences or are the result of implementation issues (e.g., line 258). Do the authors have results on additional baselines or datasets, such as standard UCI regression tasks or pose estimation? Including Gaussian processes as a baseline could help clarify whether the multimodal approach is genuinely beneficial. It would also be helpful to report additional metrics, such as RMSE, to better assess performance.\n\nI think question (1) above could be addressed if I am simply misunderstanding the paper and the authors could correct my thinking. Question (2) is certainly something that could be addressed with simple changes to the writing. Question (3/4) are also addressable but I don't know if this is something that can be done in during the short response period. Addressing 1/2 would likely me bump my score but still lean towards reject, while addressing 3/4 would lead me to lean towards weak accept. Please let me know if my criteria is not clear and I will clarify.", "limitations": "Yes.", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "These bullet points present the thoughts I had while reading the paper and more detailed questions / suggestions are contained in \"Questions\". \n\nStrengths:\n- The core idea is unique and as far as I know this particular argument of using squared families as priors is new.\n- The graphics of the paper are well done and help communicate ideas. For example, Figure 4 helps give intuition about the benefit of GSPs' multimodal behavior. \n- The paper's relation to previous work on squared families is present. \n- The presentation contains relevant details and I believe the analytic results are accurate. \n\nWeaknesses:\n- The paper is hard to follow as things are presented without context or implicitly defined (e.g. line 14). \n- The details are sometimes too much and distract from the actual point being made, see suggestion on reworking section 3. \n- The numerical results are limited and make it hard to say there is any benefit to using GSFPs, see suggestion on the experiments section.", "was_revised": true, "final_justification": "I believe the authors directly addressed my concerns for the writing changes that I suggested. However, I believe the authors did not fully address my concerns about the computational complexity since this seems a fundamental issue with the method itself. The numerical results also did not convince me of the benefits in performance or the accuracy of the comparisons. The comment I made about implementation issues (comment from line 258) was largely ignored. While I believe the paper is not yet sufficient in its comparisons I will not object if others think the paper should be accepted.", "rebuttal": "Thanks for your detailed comments and actionable suggestions for improving the paper. See our response below. Please let us know if anything needs to be further clarified or if you have any additional questions.\n\n---\n\n## 1. Computational benefits of GSFs/GSPs. \n\n(same as response to Reviewer Z38J).\nWe should clarify the distinction between (more general) GSFs applied to Bayesian inference and (more specific) GSFPs  applied to Bayesian inference in regression models.\n\nFirst, in the general case, Bayesian inference is (in)famously computationally difficult. With arbitrary likelihoods and priors, tuning computational algorithms can be practically difficult, and even fail to give good approximations unless some heavy restriction is placed on the densities (e.g. log concave). This is especially difficult when the dimension of the random variable $d$ is large, due to the curse of dimensionality, due to the exponential growth of volume. In practice, if one desires an exact computation, we need to restrict our class of probability distributions. When we restrict the class to squared families, we are able to compute exact posteriors, a thing which is usually associated with very classical models, not neural network based models. In GSFs, the complexity of this calculation is governed by the inner product of the parameter matrix $\\mathbf{M}$ and the kernel matrix $K_{\\mu, \\mathbf{\\psi}}$ (as in equation 1). Assuming the kernel matrix is known, this inner product has complexity $\\mathcal{O}(n^2)$, where $n$ is the number of parameters in the GSF. Crucially, dependence of the complexity on $d$ only appears through the calculation of the kernel matrix $K_{\\mu, \\mathbf{\\psi}}$. In practice, the kernel matrices $K_{\\mu, \\mathbf{\\psi}}$ can be computed exactly in *linear* time in $d$. (All examples in our experiments admit linear complexity in $d$. We will enumerate more explicitly such examples in the main text --- see response to point 4 of Reviewer 7C3F). This leads to a combined complexity of $\\mathcal{O}(d n^2)$, i.e. exact inference linear in dimension $d$, which is a huge improvement over inexact inference exponential in dimension $d$, which one would typically expect from numerical algorithms in general classes of probability distributions. Of course, this improvement comes with the restriction of the class of densities to GSFs, however this restriction is not too severe because GSFs are rich density approximators (as demonstrated empirically in the experiments section, as well as in previous theoretical and empirical investigations in using special cases of them as likelihood models (Rudi and Ciliberto 2021, Flaxman 2017, Walder and Bishop 2017, Loconte et al. 2023a,b, Loconte and Vergari, 2025, Loconte et al., 2024, Tsuchida et al., 2023, Tsuchida et al., 2025)). More precisely navigating the trade-off between universal approximation (making $n$ large) and computational efficiency (making $n$ small) is an important direction for future work, however we believe having $n$ as a tuneable way of navigating this trade-off is a useful perspective.\n\nSecond, in the more specific case of Bayesian inference in regression models, we are only aware of two tractable (polynomial time) methods which allow for exact inference. These are GPs (Gaussian prior with Gaussian likelihood) and NGGPs (combining a GP and a normalising flow). We would also claim that these have desirable complexities when placed in the space of all possible function inference schemes. E.g. doing inference with an exponential prior on parameter weights and a Gaussian likelihood would be computationally very difficult in practice, becuase one would require numerical methods. DKT and NGGP, as well as ours, can indeed be understood as kernel methods. The complexities are similar (although NGGP has some further complexities due to the added normalising flow component). Practically, our method is slightly slower than DKT (i.e. GPs), because of the added GSF, but faster than NGGPs. See Appendix E.9. This is why we chose to highlight it as a limitation. \n\nWe will add this discussion to the main paper, with the extra space afforded by extra one page as well as your suggestions from point 2 below. Please let us know if you need any clarifications.\n\n---\n\n## 2. Rearranging of section 3. \n\nWe agree with your suggestion. In the revision, we propose the following changes. Please let us know if you think these are sufficient.\na) Combine section 3.2 and section 3.3 \nb) Most of 3.3 does not need to appear in the main paper, and will be moved to the appendix. All that is important is using a Gaussian base measure and Gausssian likelihood.\nc) Similarly, Corollary 6 can be moved to the Appendix. Corollary 7 is what matters most, as well as a claim that the marginal likelihood is also available in closed-form (with a reference to a more precise derivation in the appendix).\n\n---\n\n## 3. Multimodality is useful.\n\nWe note that the unimodal limitation of Gaussian processes is discussed in prior work, which we cite and use as a competitor model, NGGP (Sendera et al., 2019, see second paragraph, Fig. 1, Fig. 5). Here visualisations of the failure of unimodal GPs on synthetic datasets (similar to our Fig. 4) are given. The practical limitation of GP/DKT is that it can only represent Gaussian beliefs. Any data which requires multimodal beliefs (such as our Fig. 4, or alternatively Sendera et al.'s Fig. 1, Fig. 5) will necessarily require non-Gaussian beliefs, because Gaussians are unimodal. See also our Fig. 1. We will update the manuscript to reflect this. In summary, our argument is that GPs can only represent Gaussian and therefore unimodal beliefs. Thus if a prediction may benefit from a multi-modal belief, a GP cannot realise this benefit. In contrast, a GSF is capable of representing multimodal beliefs (see Fig. 1 and Fig. 4).\n\nIt is hard to visualise such multi-modal functional beliefs on real datasets (or even the datasets themselves), because we are not aware of any with 1 input dimensions. This is why we opted to add a new synthetic fig to the existing literature to demonstrate such a belief in 1 dimension. If you have any additional suggestions, please let us know.\n\n---\n\n## 4. Experiments. \n\nWe originally considered the same few-shot learning benchmark (all datasets) as in Sendera et al., 2021. Note also that our existing QMUL dataset in section 4 is pose estimation.\n\nDuring the rebuttal period we applied DKT, NGGP and GSFP to vanilla regression datasets Boston, Concrete, Energy, Kin8nm, Naval, Power Plant, Protein, Wine and Yacht. \n\nPlease appreciate that the rebuttal period does not allow exhaustive hyperparameter searches (kin8nm took 30 hours for NNGP, 1.3 hours for GSFP), and we will update these experiments with more time beyond the rebuttal period, or caution against general takeaways. In an attempt to simplify hyperparameter search in a way which is fair to all methods, we used the same feature extractor $\\gamma$ for all methods, and shared hyperparameters across all methods. All methods share the same backbone architecture $\\gamma$, a two layer MLP. The only remaining architectural hyperparameters (ignoring batch size, learning rates, etc.) are the parameter sizes $n,m$ for the squared family and the normalising flow architecture for the NGGP. We used the original NGGP paper's architecture from Sendera et al. 2021 for the normalising flow. We used $n=3$ and $m=10$ for the squared family. All datasets use an 80-20 data split. We show mean $\\pm$ std of NLL\n\n| Method | Kernel | Dataset | | |  |  |  | | | |\n|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|\n| |  |    Boston House    |      Concrete      |       Energy       |       Kin8nm      |        Naval       |     Power Plant    |      Protein      |        Wine | Yacht  |\n|   DKL  | RBF | 0.220 ± 0.881 | 0.346 ± 0.717 | -0.298 ± 0.909 | 0.257 ± 0.693 | -1.280 ± 0.120 |   -0.041 ± 0.639   |   1.169 ± 0.765 | 1.097 ± 0.935   |   -1.784 ± 0.515   |\n| |  Spectral |    0.387 ± 1.162   |   _0.085 ± 1.565_  |  _-0.355 ± 1.039_  |   0.370 ± 0.796   |   -1.022 ± 0.198   |   -0.033 ± 0.681   |   1.178 ± 0.752   |    1.130 ± 0.723   |   -1.605 ± 0.976   |\n| | NN Linear |    0.565 ± 1.929   |    0.324 ± 1.452   |   -0.305 ± 1.036   |   0.444 ± 0.775   |  _-1.452 ± 0.391_  |   -0.043 ± 0.651   |   1.168 ± 0.711   | 1.085 ± 0.707   |   -1.317 ± 1.328   |\n| NGGP | RBF |   _0.123 ± 1.060_  |    0.743 ± 0.849   |    0.516 ± 0.516   |   0.259 ± 0.834   |    0.767 ± 0.388   |   -0.130 ± 0.312   | **0.474 ± 1.163** |   -0.825 ± 3.238   | **-3.088 ± 1.805** |\n| | Spectral |    0.361 ± 1.802   |    0.869 ± 0.468   |    0.557 ± 0.326   |   0.531 ± 0.525   |    0.562 ± 0.436   |  _-0.257 ± 0.262_  |   0.671 ± 1.236   | **-1.211 ± 4.056** |    0.525 ± 1.719   |\n| | NN Linear |    0.358 ± 1.220   |    0.543 ± 1.113   |    0.366 ± 0.450   |  _0.184 ± 0.635_  |    0.600 ± 0.428   |   -0.247 ± 0.538   |   1.278 ± 0.543   |  _-0.885 ± 4.858_  |   -0.011 ± 4.132   |\n|   GSFP   |   NN Linear  | **-0.128 ± 0.712** | **-0.150 ± 1.272** | **-1.426 ± 1.169** | **0.138 ± 0.873** | **-1.766 ± 0.180** | **-0.980 ± 1.062** |  _0.577 ± 1.030_  | 0.494 ± 1.744   | _-1.942 ± 0.941_ |\n\nNo room for details here, but we also analyse signal to noise ratio of datasets. NGGP tends to perform better with higher signal to noise (protein, wine, yacht), otherwise it may overfit. GSFP still performs well on these datasets compared with DKL.\n\nMSE is exactly the Gaussian NLL (up to an additive constant) using posterior mean of $f(\\mathbf{x})$ in the Gaussian likelihood. This is meaningful if the posterior can be easily summarised by its mean (e.g. a GP model). If the posterior is not easily summarised by a point predictor like in the case of GSFP or NGGP (e.g. the posterior mean may be in a low-density region between 2 modes), there is no single meaningful MSE value. Note that we report NLL under the full (possibly multimodal) posterior predictive distribution."}, {"review_id": "VSgEAwUoYn", "reviewer": "Reviewer_7C3F", "summary": "The paper discusses the use of squared families as priors in Bayesian inference.\nThis work considers the generalised square family (GSF) of distributions,\nshowing that it is an especially convenient family for Bayesian priors, because\n(under some constraints) (1) the marginal likelihood and posterior predictive\nbecome tractable, (2) the posterior is also GSF. These properties make GSF a good choice\nfor Bayesian regression, and in the experiment section, the authors demonstrate\nthat GSF priors can be used to outperform existing methods (DKT and NGGP) on\nfew-shot learning benchmarks.", "questions": "- You say that GSF is a rich family of distributions, but then (correct me if\n  I'm wrong) experiments are only conducted using Gaussian\n  measures/likelihoods. Are you able to try with any alternatives? Specifically,\n  the categorical likelihood would be very useful.\n- Are you able to combine GSF with e.g. deeper feature extractors? 1 layer is a\n  bit artificial. Results would be more convincing if they included a more\n  thorough hyperparameter search.\n- How exactly (algorithmically) are you performing your Bayesian updates and NN\n  updates? An algorithmic description in the experimental section would help\n  ellucidate this.", "limitations": "Limitations are adequately discussed", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 1, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strenghts:\n\n- The paper has strong results, and comprehensively demonstrates the benefits of\n  GSF priors in Section 3.\n- The experiments cleanly show benefits of GSF over alternative few-shot\n  learning methods in Table 1 and Figure 4.\n\nWeaknesses:\n\n- The introduction, especially the first subsection of squared families, reads\n  more like a background section, lacking a clear problem statement and\n  motivation for the use of GSF priors.\n- The contributions are not clear in the introduction. Two examples: first, the few-shot setting used\n  in the experiments is not disucssed at all; second, authors say 'we first\n  construct a hierarchy of families, ... (GSF)', which on first reading sounds\n  like a contribution, but then confusingly GSFs appear in the background section.\n- Claims that GSF is a rich and 'convenient to work with' family of\n  distributions (e.g. line 100/101, line 164) is not easily verifable from the\n  main text. Concrete examples would be helpful for readers.\n- The experimental results could be more comprehensive, e.g. by using\n  additional, more modern datasets/methods for comparison\n- Experiments are conducted only with a single layer of MLP. I would imagine DKT\n  and NGGP could be improved with deeper networks, potentially leading the reader to doubt the baselines.\n- Figure 4 could be much clearer, e.g. by labelling rows/columns and axes and/or\n  adding a legend.\n- There are no functioning links in the paper (e.g. figure links, table links, refs, etc.)", "was_revised": true, "final_justification": "During rebuttals, authors addressed my concerns, providing new experiments, comprehensive examples of when GSFs are tractable, improved related work (in another review), and a plan for improving the introduction/background. Due to this, I will increase my score.", "rebuttal": "Thanks for the detailed review. Please let us know if you have additional queries.\n\n## 1. Intro reads like a background. Few-shot learning is not discussed\n\nWe will add further details to the few-shot learning setting to the end of the introduction. We will move most of the first paragraphs of section 4 (adapting the text where required) to the end of section 1 in a new few-shot learning subsection. We will add motivation in terms of difficulty of Bayesian inference (see response 1. to Reviewer Z38J). Please read this response - we clarify that our main contributions concern Bayesian inference generally. Our problem statement is then to do fast and expressive closed-form Bayesian inference, and we use the tool of GSFs. As an application, we consider few-shot learning.\n\n---\n\n## 2. Definition 1 (GSFs) in background\n\nWe agree with your suggestion that this belongs in the same section as our contributions. We will move definition 1 to the start of section 3, in a new subsection. This just requires moving the heading of section 3 back a couple of paragraphs. \n\n---\n\n## 3. Concrete examples (rich and convenient to work with)\n\nUsing extra space afforded by Reviewer sKy9's suggestions as well as the extra page, we will add a short table in the main paper as well as an extended paper in the appendix which shows concrete examples of squared families with closed form normalising constant (referring to it in line 100/101 and 164).\n\n| Domain $\\Omega$ |   Base measure $\\mu$ | Feature $\\psi(\\omega) $| Reference and/or Notes | \n| ---- | ---- | ---- | ---- |\n|  |  | NN hidden layer | As per Appendix B.3. **For example:** | \n| | | | |\n| $\\mathbb{R}^d$  | $\\mathcal{N}(\\text{mean},\\text{Cov})$ | $\\text{erf}(\\mathbf{W} \\omega)$ | Williams [1996] | \n| $\\mathbb{R}^d$  | $\\mathcal{N}(\\text{mean},\\text{Cov})$ | $\\text{ReLU}(\\mathbf{W} \\omega)$ | Cho and Saul [2009] | \n| $\\mathbb{R}^d$  | $\\mathcal{N}(\\text{mean},\\text{Cov})$ | $\\text{cos}(\\mathbf{W} \\omega + \\mathbf{b})$ | Pearce et al. [2020] | \n| $\\mathbb{R}^d$  | $\\mathcal{N}(\\text{mean},\\text{Cov})$ | $\\text{GELU}(\\mathbf{W} \\omega)$ | Tsuchida et al. [2021] | \n| $\\mathbb{R}^d$  | $\\mathcal{N}(\\text{mean},\\text{Cov})$ | $\\text{Snake}(\\mathbf{W} \\omega)$ | Tsuchida et al. [2023] | \n| | | | |\n| | | Exponential | As per Appendix B.1. **For example:** |\n| | | | |\n| $\\mathbb{R}^d$  | $\\mathcal{N}(\\mathbf{0},\\mathbf{I})$    | $\\exp(\\mathbf{W} \\omega+\\mathbf{b})$ | e.g. Nielsen and Garcia [2024] | \n| $\\{ 0, 1, 2,\\ldots \\}$  | $(\\omega !)^{-1} \\nu$    | $\\exp(\\mathbf{W} \\omega+\\mathbf{b})$ | e.g.  Nielsen and Garcia [2024] | \n| $\\mathbb{S}^{d-1}$  | Uniform    | $\\exp(\\mathbf{W} \\omega+\\mathbf{b})$ | e.g.  Brown [1986] | \n\n(**Caption:** All of these settings admit closed-form squared family kernels $\\matrix{K}_{\\mu,\\psi}$ in linear time in the dimension $d$. Here erf denotes the error function, $\\mathbb{S}^{d-1}$ denotes the unit hypersphere, Snake denotes the snake activation function and $\\nu$ denotes the counting measure. Further examples are given in Han et al. (Table 1, 2022) and Tsuchida et al. (Table 1, 2023).)\n\n(Note: $\\Omega$ is the domain of the distribution. Openreview does not allow typesetting as in paper.)\n\nWe are also happy to type some of the equations of the closed-form kernels in the text. Please let us know if you would like this as well. To give you an idea of what this would look like, for the Cho and Saul [2009] row in the table above, the $ij$th entry $ \\mathbf{K}_{\\mu,\\psi}(i,j)$ is\n$$ \\mathbf{K}_{\\mu,\\psi}(i,j) = \\frac{\\Vert \\mathbf{w}_i \\Vert \\Vert \\mathbf{w}_j \\Vert }{2\\pi}\\left( \\sin  \\rho - (\\pi - \\rho) \\cos\\rho  \\right), $$\nwhere $\\rho$ is the angle between $\\mathbf{w}_i $ and $\\mathbf{w}_j $, i.e. $\\rho = \\cos^{-1} \\frac{\\mathbf{w}_i^\\top \\mathbf{w}_j  }{\\Vert \\mathbf{w}_i \\Vert \\Vert \\mathbf{w}_j \\Vert}$. The complexity of this calculation is only linear in the dimension $d$ of the distribution (which is the same as the number of entries in $\\mathbf{w}_i$ or $\\mathbf{w}_j$), because we just need to compute inner products and norms. The same linear complexity holds for all examples in the table.\n\n---\n\n## 4. Other methods \n\nWe compared with probabilistic methods for function inference (GPs/DKTs and NGGPs) with exact normalising constant. These are the only competitor methods we are aware of with these properties. Sendera et al. only compared DKT/NGGP. If you can suggest other methods with these properties, we may be able to add if we have time, or at the very least discuss similarities/differences/strengths/weaknesses.\n\n---\n\n## 5. Other datasets \n\nWe used all of the same datasets as Sendera et al. for few shot learning. See our response to Reviewer sKy9 for **9 additional** experiments on vanilla regression, with standard UCI datasets.\n\n---\n\n## 6. Number of layers\n\nNote that experiments are **not** conducted with only a single layer MLP. There is a distinction between the feature mapping $\\mathbf{\\gamma}$ (which is a deep network, common to all methods, and may be a 3 hidden layer conv net or 2 hidden layer MLP depending on the dataset) and the function $\\psi$ which only appears in the GSF, and is a single hidden layer network. See Appendix E.7. Note that we use the same feature mappings for DKT and NGPP as in respective original papers (Patacchiola et al., 2020 and Sendera et al., 2021].\n\n---\n\n## 7. Links\n\nApologies, for this oversight. The final paper will have functioning links.\n\n---\n\n## 8. GSF are rich claim\n\nOur claim is that the GSF priors/posteriors are rich, i.e. they can model non-Gaussian and multimodal distributions. This is not a statement about likelihoods - indeed, we do not claim that GSFs are conjugate to all likelihoods. In the case of regression, we rely on a Gaussian likelihood and our results follow from its conjugacy with the Gaussian base measure -- but note that the priors and posteriors have a complex neural-network based non-Gaussian density with respect to that base measure (this is what makes them expressive while preserving conjugacy). We also analyse other likelihoods (e.g. Poisson in Figure 1, others in Appendix B).\n\n---\n\n## 9. Algorithms\n\nWe will add the following algorithms to the paper (Formatted better in Latex Openreview doesn't allow some things. We skip some  algorithms here due to space, but these are all that are required for GSFP few shot learning). In section 4, we run algorithm 3 on the meta dataset then algorithm 2 on the support/query dataset.\n \n**Alg. 1:** ``Marginal likelihood`` (i.e. Proposition 2). **Inputs:** Trainable final layer parameter $\\mathbf{\\Theta} \\in \\mathbb{R}^{m \\times n}$, Trainable hidden network $\\psi:\\mathbb{R}^d \\to \\mathbb{R}^n$, Trainable base measure $\\mu$, Frozen likelihood $p(\\mathbf{U} \\mid \\omega)$. **Output:** Marginal likelihood $p(\\mathbf{U}) \\in [0,\\infty)$.\n\n0. Set $\\mathbf{M} = \\mathbf{\\Theta}^\\top\\mathbf{\\Theta}$.\n1. Compute $\\mathbf{K}_{\\mu,\\psi}$ (see Appendix B, enlarged with response to Question 4 above).\n2. Set $\\nu(d\\omega) = p(\\mathbf{U} \\mid \\omega) \\mu(d \\omega)$. Compute $\\mathbf{K}_{\\nu,\\psi}$ (see Appendix B, enlarged with response to Question 4 above).\n3. Set $z(\\mathbf{M}) = \\text{Tr}(\\mathbf{M} \\mathbf{K}_{\\mu,\\psi})$, as per (1).\n4. Return $p(\\mathbf{U})=\\text{Tr}(\\mathbf{M} \\mathbf{K}_{\\nu,\\psi})/z(\\mathbf{M}),$ as per Proposition 2.\n\n**Alg. 2:** ``GSFP posterior predictive update`` (i.e. Corollary 7). **Inputs:** Frozen Gaussian likelihood $p(\\mathbf{Y} \\mid  \\gamma(\\mathbf{X}) \\omega)$ on support set. Query $\\mathbf{X}_\\ast$. Prior element of a GSF, parameterised by: Final layer parameter $\\mathbf{\\Theta} \\in \\mathbb{R}^{m \\times n}$, Hidden GSF network $\\mathbf{\\psi}:\\mathbb{R}^d \\to \\mathbb{R}^n$, Gaussian Base measure $\\mu$, Deep feature extractor $\\gamma:\\mathbb{X} \\to \\mathbb{R}^d$ **Output:** Posterior element of a GSF, parameterised by: Final layer parameter, Hidden network and Base measure.\n\n0. Set $\\mathbf{M} = \\mathbf{\\Theta}^\\top\\mathbf{\\Theta}$.\n1. Set $\\nu''(d\\omega)$ to be Gaussian posterior predictive, as per (6).\n2. Set $\\mathbf{\\psi}'$ as per Corollary 7.\n2. Set $z(\\mathbf{M}) = \\text{Tr}(\\mathbf{M} \\mathbf{K}_{\\nu'',\\mathbf{\\psi}'})$, as per (1).\n3. Return Posterior element of a GSF, parameterised by $\\mathbf{M},\\nu''$ and $\\mathbf{\\psi}'$. That is, \n$P(d \\mathbf{f} \\mid \\mathbf{M},\\nu'',\\mathbf{\\psi}') =\\Big( \\text{Tr}( \\mathbf{M} \\psi'( {\\mathbf{f}} ) \\psi'( {\\mathbf{f}} )^\\top )/z(\\mathbf{M}) \\Big)\\nu''(d\\mathbf{f})$, as per (1) and Proposition 2.\n\n($\\mathbf{f}$ is shorthand for $\\mathbf{f}(\\mathbf{X}_\\ast)$).\n\n**Alg. 3:** ``GSFP Few-shot learning pretraining`` (i.e. paragraph 1 of section 4). **Inputs:** Meta datasets $(\\mathbf{X}_j,\\mathbf{Y}_j)$ for $j = 1,\\ldots, J$, Batch size $B \\in \\{1,\\ldots,J\\}$, Trainable final layer parameter $\\mathbf{\\Theta} \\in \\mathbb{R}^{m \\times n}$, Trainable hidden network $\\mathbf{\\psi}:\\mathbb{R}^d \\to \\mathbb{R}^n$, Trainable Gaussian base measure $\\mu$, Deep feature extractor $\\gamma:\\mathbb{X} \\to \\mathbb{R}^d$ **Output:** Trained GSFP prior distribution.\n\n0. For batch $b = 1,\\ldots,\\lfloor J/B \\rfloor + 1$:\n\n    0.a Set $\\mathbf{U}_i = (\\mathbf{Y}_i,\\gamma(\\mathbf{X}_i) )$, for each index $i$ in the batch. Compute Gaussian likelihood $p(\\mathbf{Y}_i \\mid  \\gamma(\\mathbf{X}_i) \\omega)$, as per (4), for each index $i$ in the batch.\n\n    0.b Compute gradients of ``Marginal likelihood``($\\mathbf{\\Theta},\\mathbf{\\psi},\\mu, p(\\mathbf{Y}_i \\mid  \\gamma(\\mathbf{X}_i) \\omega)$) for each $i$ in the batch with respect to $\\mathbf{\\Theta}$ and parameters of $\\mathbf{\\psi},\\mu$ and $\\gamma$ using automatic differentiation.\n\n    0.c Update $\\mathbf{\\Theta}$ and parameters of $\\mathbf{\\psi},\\mu$ and $\\gamma$ using gradients (e.g. using SGD update rule). \n\n    0.d Store gradient information if required (e.g. for Adam etc.).\n\n1. Return $\\mathbf{\\Theta},\\mathbf{\\psi},\\mu,\\gamma$.\n\n---\n\n## References\n\nBrown (1986). Fundamentals of statistical exponential families: with applications in statistical decision theory."}]

Author rebuttal:

Meta-review: Reviewers find this paper to be strong, especially the fact that its results are novel (Reviewer 2rpL, Reviewer sKy9), generalize fundamental classes of models such as mixture models (Reviewer Z38J), and that the paper is well-written with good visuals (Reviewer sKy9), and comprehensive experimental results (Reviewer 7C3F). The main weakness raised is by Reviewer 7C3F, who notes that the paper lacks a clear problem statement: I think one can view this as more of a property of the work's nature, than a weakness. Specifically, this paper is not about solving a specific problem, but rather about introducing and developing a novel general-purpose technique that can be used in many problems. In such a paper, it is important to develop the components of the technique to an adequate standard, which is exactly what the paper does in the form of deriving the key formulas for marginal likelihoods and other quantities of interest, and noting the key properties that a reader might be interested in using. There is great precedent for papers like this ending up useful and impactful over the long run, potentially in ways the authors could not have initially considered, as fundamental research in new model classes creates new kinds of tools for people to use in various ways. As consequence, I am convinced that while Reviewer 7C3F presents this as a weakness, it could simultaneously have been presented as a strength under different framing.

Given these points, and on basis of reviewer consensus, I recommend acceptance.
Decision: Accept (poster)