Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
Authors: Mehmet Yiğit Balık
Maksim Sinelnikov
Priscilla Ong
Harri Lähdesmäki
Venue: ICML 2025 poster (2025)
Keywords: Gaussian Process
High-dimensional time-series
Variational Autoencoder
Conditional Generation
Reproducing Kernel Hilbert Space
Primary area: deep_learning->sequential_models_time_series

Abstract: High-dimensional time-series datasets are common in domains such as healthcare and economics. Variational autoencoder (VAE) models, where latent variables are modeled with a Gaussian process (GP) prior, have become a prominent model class to analyze such correlated datasets. However, their applications are challenged by the inherent cubic time complexity that requires specific GP approximation techniques, as well as the general challenge of modeling both shared and individual-specific correlations across time. Though inducing points enhance GP prior VAE scalability, optimizing them remains challenging, especially since discrete covariates resist gradient‑based methods. In this work, we propose a scalable basis function approximation technique for GP prior VAEs that mitigates these challenges and results in linear time complexity, with a global parametrization that eliminates the need for amortized variational inference and the associated amortization gap, making it well-suited for conditional generation tasks where accuracy and efficiency are crucial. Empirical evaluations on synthetic and real-world benchmark datasets demonstrate that our approach not only improves scalability and interpretability but also drastically enhances predictive performance.

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
Mehmet Yigit Balik Maksim Sinelnikov
Priscilla Ong
Harri Lähdesmäki
Abstract
healthcare domain. By studying the associations between covariates and samples (i.e., measurements), it is possible to model temporal dynamics, such as disease progression, and impute missing values in the data (Diggle, 2002). However, analysis challenged by factors such as high dimension- ality, non-trivial correlations within and across subjects, time-varying covariates, and missing values and observa- tions (Ramchandran et al., 2021).
High- digh-dimensional time-series datasets are com- mon in domains such as healthcare and economics. Variational autoencoder (VAE) models, where la- tent variables are modeled with a Gaussian pro- (GP) prior, have become a prominent model class analyze such correlated datasets. How- ever, their applications are challenged by the in- herent cubic time complexity that requires spe- cific GP approximation techniques, as well as the general challenge of modeling both shared and individual-specifee correlations across time. Though inducing points enhance GP prior VAE scalability, optimizing them remains challeng- ing, especially since discrete covariates resist gradient-basse methods. In this work, we pro- pose a scalable basis function approximation tech- nique for GP prior VAES that mitigates these chal- lenges and results in linear time complexity, with global parametrization that eliminates the need for amortized variational inference and the as- sociated amortization gap, making it well-suited for conditional generation tasks where accuracy and efficiency are crucial. Empirical evaluations on synthetic and real-world benchmark datasets demonstrate that our approach not only improves scalability and interpretability but also drastically enhances predictive performance.
Several machine learning methods have been developed to address these challenges, with variational autoencoders (VAEs) emerging as an established class of models for representation learning and generative tasks (Kingma Welling, 2014; Rezende et al., 2014). VAEs can learn both expressive neural network parameterized model for high- dimensional data generation well as an efficient amor- tized variational approximation of sample-specific latent variables. However, standard VAEs assume that each data point is independent and identically distributed (iid), which limits their ability to capture correlations between samples in temporal datasets.
VAEs can be enhanced model correlations between sam- ples. Among these enhancements, VAFs equipped with a Gaussian process (GP) prior for latent variables have proven to be particularly effective for capturing time dependencies and correlations between samples and achieve impressive performance improvements (Casale al., 2018). However, this powerful approach comes with cost of computational complexity that scales cubically with the size of the dataset Ceinonero-Candell & Rasmussen, 2005; Bauer et al., 2016) that limits its scalability for large datasets. There have been efforts to tackle these challenges by incorporating inducing points (Titsias, 2009) and stochastic gradient-based opti- mization with mini-batching (Hensman et al., 2013). Opti- mizing inducing point locations becomes non-trivial with discrete covariates their non-continuous kernels preclude gradient-based updates. To our knowledge, no efficient general solution for discrete covariates in GP prior VAEs has been proposed. Finally, their variational posterior com- monly uses a mean-field approximation with parameters tied by a shared encoder, which limits performance and introduces an unavoidable amortization gap (Cremer et al., 2018; Margossian & Blei, 2024).
1. Introduction
High-dimensional time-series datasets are common var- ious domains, including social sciences, economics, and healthcare. These longitudinal datasets contain repeated measurements of subjects over time along with auxiliary covariate information, such as time, age, and gender, in the
'Department of Computer Science, Aalto University, Es- poo, Finland. Correspondence to: Mehmet Yigit Balık waatinettaaiii@aalllooo
Proceedings of the 4rnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
(b) Exact udditive GP
((x Σ (z)) (()) GP(0, K
(d) High-dimensional generation
Py(y
Alliog gratel
Al"ddoooiione divess
(a) Conditional generative model with covarinte intteeiifoonnn information a
3
3
3
(c) Basis function approximation
Σ
20,)
x0)
2
a(3)
14
20
10
Time
Figure 1. Overview of DGBFGP. (a) A conditional generative model incorporates covariate information. (b) Additive GP priors are used but not scalable. (c) We replace exact GP with basis function approximation for scalability. (d) Samples are generated by the decoder.
Contributions. In this work, we propose Deep Genera- tive Basis Function Gaussian Process (DGBFGP), a novel conditional generative model. Inspired by earlier VAE ap- proaches (Ramchandran et al., 2021; Ong et al., 2024), DG- BFGP capable of accurately modeling high-dimensional longitudinal data by capturing complex correlations within and across subjects, specified by extensive continuous and categorical covariates using scalable basis function approx- imation that provides global parameterization and learns parameters using variational inference (VI) (Hoffman et al., 2013). Figure overviews DGBFGP, which uses scalable GP approximations with auxiliary covariates model in- dividual, group, and temporal effects. The model delivers flexible. interpretable representations and scales to large datasets. Our contributions are:
2. Related Work
There has been vast array of suggestions to improve the performance of VAEs. Sohn et al. (2015) proposed the con- ditional VAE (CVAE), which leverages auxiliary covariate information to better guide the latent space representation and generate more context-aware outputs. However, uses standard Gaussian prior, which is not flexible for com- plex tasks. The expressiveness of approximated posteriors can be strengthened by importance sampling (Burda et al., 2016), normalizing flows (Rezende & Mohamed, 2015), RealNVPs (Dinh et al., 2017), and inverse autoregressive flows (Kingma et al., 2016). Yet, these methods still rely on the iid assumption across samples. Hence, they cannot capture sample correlations.
Casale al. (2018) proposed the GP prior VAE (GPPVAE), which incorporates GP priors to include view and object information directly into the latent variables. However, its reliance on restrictive view-object GP product kernel lim- its its ability to capture subject-speciff temporal patterns, making it less suitable for longitudinal studies. Additionally, the pseudo-mini-bant training scheme lacks scalability for large datasets as it ne- cessitates full dataset pass for every training step.
We propose DGBFGP, a conditional deep generative model that leverages basis function approximation for mixed-domain additive GPs to model complex correla- tions within and across subjects, guided by extensive continuous and categorical covariates.
stochastic gradient descent (SGD)
Our global parameterization avoids explicit kemels, runs in linear time, eliminates the amortization gap, overcomes the limitations of categorical inducing point optimization, enhances interpretability in the latent space by quantifying the contributions of different co- variates/effects using Sobol indices, allows for standard mini-batch training, and treats kernel hyperparameters probabilistically.
Building upon the idea of employing GPs within the latent space of VAES, Fortuin et al. (2020) introduced GP-VAE which assumes a separate GP prior for the time-series data of each individual subject. While GP-VAE is tailored for temporal data, its independent GP priors for each subject limit its ability capture shared temporal patterns. Further- more, GP-VAE cannot utilize auxiliary covariate informa- tion beyond time, restricting its flexibility in context-rich
We compare DGBFGP against state-of- of-th--f-t methods in the literature and report highly competitive perfor- mance on synthetic and real- world benchmark datasets.
2

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
Table 1. Features of our proposed method contrasted the related previous methods.
MODEL
SHARED TEMPORAL STRUCTURE
INDIVIDUAL TEMPORAL STRUCTURE
OTHER COVARIATES
MINIBATCHING
GENERATIVE
REFERENCE
COMPLEXITY
CVAE
SOHN ET AL. (2015)
BRITS
ET A1., (2018)
L-NODE
JAKOVLEV ET AL. (2025)
GPPVAE
LIMITED V LIMITED
LIMITED X LIMITED
PSEUDO
ONNH")
CASALE ET AL.
GP-VAE
X
4
O ng)
FORTUIN ET AL. (2020)
VAG
O(NA + M1"
JAZBEC ET AL. (2021)
LVAE
O(NM + -.
(2)
RAMCHANDRAN ET AL. (2021)
SGP-BAE
LIMITED
O(NM 4 M))
TRAN ET AI. (2023)
DGBFGP
O(N E, B
OUR WORK
datasets. To overcome some of these issues, both SVGP- VAE (Jazbec et al., 2021) and LVAE (Ramchandran et al., 2021) were introduced, leveraging inducing points to enable scalable, mini-batch training (Titsias, 2009; Hensman al., 2013). While SVGP- -AE inherits GPPVAE's limitations in incorporating all auxiliary covariates, LVAE is capable of handling an arbitrary number of covariates, allowing greater flexibility in modeling diverse datasets. However, both ap- proaches face challenges in initializing and optimizing the inducing points, which can impact their performance and ap- plicability, particularly in complex scenarios. Finally, Tran al. (2023) proposed SGP-BAE, fully Bayesian autoen- coder using stochastic gradient Hamiltonian Monte Carlo to jointly sample decoder and sparse GP prior parameters. Howver, it still faces the same challenge with categorical covariates due to its reliance on inducing points.
mentary aspects of time-series analysis and offer diverse tools for high-dimensional datasets. Table summarizes the capabilities of the related previous methods.
3. Background
Problem setup. We consider general longitudinal experi- mental designs with P unique instances (i.e. subjects or pa- tients), where each instance has np longitudinal samples. The total number of samples equals the sum of instance- specific samples, - >pp Rp. We denote the high- dimensional samples as E y = RDD. Without loss of generality, we can assume that the samples are ordered according to the instances such that the first nı samples ......... are from instance one, the next n2 samples yn, In.+1:-... tnytng from instance two, etc. Samples across all instances are denoted Y (-1,-... . y v Addition- ally, instance p has corresponding covariate information for each sample, collectively denoted as X [.i,.... 2x NN.. where æ E X' I Xr is Q-dimensional vector includ- ing information such as time, age and gender. We denote the sample-specific L-dimensional (L < D) latent variables as Z [ NN where each z E Z = R'.
Multi-output GPs extend traditional GPs to model corre- lations across multiple output dimensions (Alvarez et al., 2012). By learning dependencies both within and across variables, they offer natural and flexible way to represent coupled latent spaces. This capability makes them particu- larly suitable for tasks where multiple related outputs need to be jointly modeled. However, multi-output GPs do not directly scale to high-dimensional datasets.
3.1. Variational autoencoders
Ong et al. (2024) introduced LMM- which uses linear mixed models (LMMs) as conditional priors in VAEs, en- abling the integration of auxiliary covariates while modeling shared and random effects. By combining the strengths of LMMs and VAEs, it captures nuanced dependencies and is adaptable to diverse applications. Its performance can be further improved with basis functions for greater flcxibility in modeling complex patterns.
Consider a data sample y generated by a joint gener- ative model Pw ( = p=(l)l))p(()) where w (4, 0} denotes parameters. The posterior (=(y)) = Pe (yl ())))=)/pp)))) is generally intractable since the term (w) cannot be computed analytically. A common ap- proach is to employ amortized VI (AVI) (Kingma & Welling, 2014; Rezende et al., 2014) to estimate the true posterior with a parameterized approximation (=(2))). This is done by optimizing a lower bound on the evidence, log (w) w.r.t. all parameters involved
Modeling temporal correlations is also explored through latent neural ordinary differential equations (ODEs) (Rubanova et al., 2019). Iakovlev et al. (2023) proposed L-NODE, latent neural ODE with Bayesian multiple shoot- ing that offers an efficient framework for continuous-time latent dynamics with sparse observations. contrast, clas- sical approaches like BRITS (Cao et al., 2018) use bidi- rectional recurrent networks for fime-series imputation but not generative models. These methods address comple-
log Pa (1) 2 Egg Poogppp y 2)]] (L 3)) || po po))).
A common assumption is that the likelihood P.p (y) =)) the prior (=)),, and the variational posterior (z(2 y)) are all represented using a mean-field approximation. In practice, this implies that distributions factorize over samples.
3

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
Sohn al. (2015) build upon this idea by incorporating auxiliary covariates
generality because the latent variables are further mapped the data space by a neural network parameterized decoder Pe (3 =) that can introduce arbitrary correlations at least as expressive as standard multi-ouiput GPs across dimensions.
X into the encoder and decoder networks. This allows the VI procedure to approximate not just Pw (z v), but Pa (z y. x) for the conditional generative model Pw (y, =/2x) = Po (y) z, a: p)) (z|2z) Although CVAE enables more flexible generative modeling by leveraging known factors, it shares the iid assumption of the VAE.
Motivated by longitudinal datasets, Ramchandran et al. (2021) proposed to use an additive GP prior for each of the latent dimensions:
3.2. Gaussian processes
f((x) = (x +...... 51 (x0((x)( (R) (-())) ~ GP (0. (())) (x!"), (r)' | ))))
p(
GPs extend the concept of the Gaussian probability distri- bution of random variables to stochastic processes, which characterize the properties of entire functions (Rasmussen & Williams, 2006). A univariate GP
where R denotes the number of additive components, each component depends on either single or pair of covari- ates (+() xx)) S X, and 8/t) denotes the kernel hy- perparameters of the ,4th component of the 1 latent dimen- sion. This additive GP model provides interpretability sim- ilar to commonly used LMMs (Laird & Ware, 1982) do and allows for handling mixed -oxedd (e.g. a combina- tion of continuous and categorical) inputs. The additive GP model implies that the factorizable joint distribution p(f(X) 0) from Equation (1) has an additive covariance
f(x) ~ k(z,,"
is fully determined by its mean function jux²) (gen- erally considered as 0) and covariance (i.e. kernel) function 6(z,x'), For a finite collection of inputs X [2 (.1,..., V.J. GP evaluated at X, [(X) = [f(=))....,,( f(xx)] is defined to have joint Gaussian distribution f(X)= ~ E), where E. k(z,,æ,). The covariance function encodes our prior assumptions about the underlying latent functions that generate the data. For example, a kernel for continuous-valuee covariate might assume that the function changes smoothly (ensuring con- tinuity), whereas a kernel for categorical covariate defines properties of the function across the categories. These co- variance functions generally have hyperparameters, which are inferred from data.
structure Σ, E (r) where i,j element of the co- variance matrix V is given by ky (x)" x.)" (
Being able to model correlations comes with chain of trade-offs. We need to keep a kernel matrix with ((2)) memory complexity. In addition, O(N3) time complexity arises from the inversion of the kernel matrix, necessitating reliable approximative GP techniques. Previous GP prior VAES have proposed to use the inducing point approach (Ramchandran et al., 2021; Jazbec et al., 2021) achieving ( + 2) time complexity with M inducing points, or low-rank approximations (Casale et al., 2018) reducing complexity to ONN where H < N is the dimension- ality of the low-rank approximation. Instead, in this work, we propose an efficient method that uses basis function approximation, which results in global parametrization.
3.3. Gaussian process prior variational autoencoders
A GP conditioned auxiliary information can capture non- trivial dependencies and provide a flexible, non-parametric model that can adapt to complex patterns in the data. Pre- viously, researchers have utilized GPs as priors in VAEs to create effective generative models (Casale et al., 2018).
Let f : X Z and f(x)= (z),, (.))).... Sz (T)) and assume that f follows multi-output GP prior (Alvarez et al., 2012):
4. Methods
4.1. Kernel approximation for continuous covariates
f(x) ~ GP(O, K (x x'| 0)),
Stationary GP kernels, such as the squared exponential (SE)
where K(x, x'| |8] a so-called cross-covariance function with hyperparameters 0. Contrary to standard multi-output GPs that typically use a linear model of co-regionalization to define the cross-covariance function (Alvarez et al., 2012), GP prior VAE models assume that each latent dimension has an independent GP prior, i.e., factorizes across dimensions
I===" ), can be writ- kernel (2, x' or E) = 0² exp([=l 2
ten as a function of the distance between the two inputs, k(r) = k(æ, x²) = k(I= - x')), where x E R is a univariate continuous input covariate (e.g. time, age, etc.). According to Bochner's theorem (Akhiezer & Glazman, 2013; Ras- mussen & Williams, 2006), covariance functions can be represented as the Fourier transform of a positive measure since they are positive definite. When a measure has an associated density, it is referred to as the spectral density s(w) of the covariance function. This relationship leads to the Fourier duality between spectral densities and kernels, a
p(f(x) 0) = II NN Ex),
(I)
where Σn denotes the covariance matrix of the th compo- nent of f()). This assumption comes without any loss of
4

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
principle known the Wiener-Khintchin theorem:
s(w) = F(k(r)} = k(r)e
- iwr dr
S
k(r) =())==(s(ww) = 2T I s(w)elurdas dw.
yn
For example, the spectral density corresponding to the SE kernel is given by (u a. = 2n exp( (nm which is Gaussian in the frequency domain.
N
2
Figure 2. Graphical model of DGBFGP. Blank, partially shaded, and shaded circles indicate latent, partially observed, and observed variables, respectively. Non-circled variable S is deterministic.
Solin & Särkkä (2020) demonstrated that an isotropic co-
variance function (k(r) & k(([=l)) has the following eigen- function approximation
obtain the Hilbert space approximation for the y compo- nent of multi-output GP prior by using the approximation from Equation (3) for each latent dimension separately
k(x,x') ~ k(x
Σ s(v (m)()))))) (x'), (2)
where Am and &m are the thth cigenvalue and eigenfunction of the Laplace operator with Dirichlet boundary conditions, including boundary Ω = [] E IR described as
ar1
x(¹) = AA
:
((x(1)),
L
op (m)( = Am@m (z), Ox
x E Ω
where a,l IR-M denotes the linear model parameters for ad- ditive component and latent dimension I, ((x(1)) E RM and a, ~ NO, (p ((r:,,r)) (.))- Here we assume that ker- nel hyperparameters are shared across dimensions although this can be easily generalized, and we explicitly write the prior dependence of art on the corresponding kernel hy- perparameters a, and &, via the spectral density. Since GP prior factorizes across latent dimensions, the prior on A(r) is defined as p(A(rr) I p(a.). Therefore, the full additive GP structure is approximated by
0m (I)) = 0,
x E 82.
In this setup with J > 0, the eigenfunctions and eigen- mm(x+J)) and Am ( amm )2
values are (x) J sin( 2J
respectively. Finally, we can represent Hilbert space approx- imation of GP in a linear parametric form with (xxx)
[dr (z))......, [M(]]] and global parameters R.
f(x) ~ a ((x),
(3)
where a ~ (00 S) and S = diag (s( (<((\s))))) ( (VAM)) Note that the approximation depends on the kernel hyper- parameters via the parameter prior, which for the SE ker- nel can be written as a (0,, E)) and ((a, diag(sse ( (v11 o, .)...... Sae ( (V>as 1)).
R
R A (((())
7=1
which is parameterized by a collection of weights A = (A'),... A()) with prior p(A) = I ((A See Appendix for details. The generative model can then be expressed follows (see Figure 2 for the plate diagram):
For the remainder of this paper, we will use the SE kernel to model continuous covariates. However, our method is applicable to any stationary continuous covariance function.
a, Lognormal 1)
(r = I,...,
4.2. Generative model
E, Lognormal(0,1)
Ong et al. (2024) show that a Fourier basis-based linear model can be utilized to approximate the latent space of GP prior VAEs. Here we show how an additive GP prior VAE model can be defined using the Hilbert space approximation and the associated global parametrization. Let x (r) R de- note one of the continuous covariates from X' as before. To simplify notation, we assume that all additive components use the same number of cigenfunctions, M, for the kernel approximation in Equation (2), but that can be generalized easily. Assuming the latent dimensions are a priori indepen- dent as in previous works (see e.g. Equation (1)), we can
S, 4x) diag( (8eo ( V Am or, er )
R III N(an (,,, (4)
A
r=11=
N
Y 1|A,X
II
Po 1 a,
R=1
where o =
R), 2= ........ RR. and (ym
or, f.) is again the spectral density of the SE kernel that depends on the kernel hyperparameters o, and e, evaluated
5

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
at the square-root of the eigenvalue Am- A log-normal prior is chosen for both the length scale & (Timonen al., 2021) and amplitude or it ensures positivity and accommo- dates wide range of values. The likelihood function Pe is parameterized by a neural network.
are both B dimensional. The GP with interaction kernel given in Equation (5) has the following parametric form f(x) = a J (()), a ~ (0( diag(s)).
We showed in Sections 4.1-4.2 that a,s. the global weights for continuous variables, have multivariate Gaussian prior with diagonal covariance in Equation (4). This structure extends directly to the interaction kernel by incorporating a diagonal covariance with entries that combine the SE ker- nel's 's spectral density and the eigenvalues of the categorical kernel. The basis functions are products of the continuous and categorical kernel bases as shown in Equation (6).
4.3. Interpretable model for mixed-domain covariates
We have so far described a conditional deep generative model with a highly scalable GP prior for continuous covari- ates. But it is also possible to model categorical covariates Merchhn & Hernández-Lobiad 2020)
and kernels utilizing their eigendecompooitioon. Assume that one of the covariates, x E C (1..... C), has C' possible categories, and the corresponding kernel is k C x C - V C R. Let C' denote the C' x C' kemel matrix. Symmetric kernels have the cigendecomposition C' = O where D is a diago- nal matrix containing C' eigenvalues on the diagonal and θ has C' eigenvectors as its columns. Following Timonen & Lähdesmäki (2024), by defining function 9. C - E that outputs the values of the 2th column of θ, the covariance function & can be written as
4.4. Inference
We use VI (Hoffman et al., 2013) to learn the parameters of the generative model by assuming factorized variational posterior, ((4, . E) = ((A))(())))) (C)). Approximate poste- rior ((A) is Gaussian whereas q(or) and q(e) are log-normal. The ELBO is given as (see Appendix C for the derivation)
logp(Y |x) 2 E (A(A)
log Po A, x..)
xx) = Cz,Z' = DD 1 Σ
de Oz
Egtorate) (1( p(A of e)) (7)
Σ dcvc())0,(''..
KL(((o)) |()())
- -KL(g(E) ()))),
The corresponding GP has the parametric form f(z) = ((r)a, where (x) denotes the :t row of θ, and a ~ (O,, D). It is important to note that fx) is the exact GP, not an approximation. Hence, it can be used to model categorical covariates, such as gender and subject IDs. Gen- eralization for new subjects can be achieved by learning only subject-specific parameters (see Appendix B.
where we use Monte Carlo sample and l from the cor- responding distributions q(or) and ((e) in order to estimate (g( (((((A) p(A o E))] These values are then substituted into the spectral density computation. Note that the spectral densities are needed to define p(A A or, L) as shown in Equation (4). The closed form of this KL term is given in Appendix C, Equation (11). which offers simpler solution compared to the more complex formulations involv- ing GPs in previous GP prior VAE models. Additionally, the KL terms for and l are provided in Equation (12). Since the parameterization is global, the ELBO (7) is directly amenable to SGD optimization with mini-batching.
We are typically interested in modeling joint effects of con- tinuous and categorical covariates (Saves et al., 2023). In our setting, we can compute the interactions between the ap- proximation for the stationary kernel and the decomposition of the categorical kernel (Timonen & Lähdesmäki, 2024):
k( x²) (a,, x Rica (x.
We also present a semi-amortized version of our model, denoted as SA-DGBFGP in Appendix E, where AVI is used only for modeling instance-specific characteristics.
M C
(5)
Σ
Σ
s( ( (2)))
1==1
4.5. Computational complexity
where Reee (x. x²) depends only on one continuous variable and hica (ar, .²) depends only on one categorical variable of inputs. Defining B MC., the product kernel can then be expressed as k(xxx²') Σm (2))) where
The Hilbert space approximation admits closed-form eigendecomposition independent of specific kernels, with hyperparameters entering only via spectral densities that also have closed-form solutions, enabling ((1) plug-and- play use of eigenvalues and eigenfunctions. Computing the categorical kernel eigendecomposition once suffices, so overhead is negligible as C < N. Hence, the most com- putationally demanding operation in DGBFGP model is
(z)
s(VAi di s(VA1)d2
(1())((())
<b(x))=
and s
(6)
(2
s( s(VAns)de)
6

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
the multiplication of each data point's basis function rep- resentation, (( (r)) E RRP, with the global linear model parameters, art E R B This computation is performed for
Table 2. Temporal interpolation MSF Rotated MNIST test data.
Method
Latent dimension
MSE
N data points. Extending this operation across all compo- nents while accounting for the possibility that each compo- nent may have different number of basis functions, the total time complexity for single latent dimension becomes O(N Σm=1 rr= (())
BRITS (Cao al., 2018)
0,989 0.0098
SVGP- (Jazhec al., 2021)
16
0.028 + 0.0005
LVAE (Ramchandran et al., 2021)
16
0.026 + 0.0006
SGP-BAE (Tran et al.. 2023)
16
0.023 = 0.0006
L-NODE (lakovlev el al., 2023)
16
0.017 1: 0.0002
DGBPGP (our work)
16 0.009
0.0001
4.6. Sobol indices for quantifying interpretability
5.1. Rotated MNIST
In additive GPs, we can use a sensitivity analysis based on Sobol indices (Sobol, 1993) to quantify the relative con- tribution of additive components (Lu et al., 2022). The Sobol index for 7thh component is defined as SI(e) =
We evaluate our model on the temporal interpolation task using a modified version of the handwritten MNIST dig- its dataset (LeCun al.. 1998). This dataset consists of P 400 unique instances of the digit '3', each rotated across 16 evenly spaced angles covering a full rotation, [0, 2T.. Each data point is characterized by unique instance id and an associated rotation angle. Temporal interpolation in this context refers to predicting intermediate representa- tions of an instance at unseen angles. See Appendix F.1 for additional details.
(a("))
VE [[]]))]]] which is the fraction of the total GP out-
put variance attributable to the 7hh component, thereby of- fering clear metric for assessing its impact. is important to note that in this context, the analysis is focused the contributions of the model components, not directly on in- dividual covariates. A large ST()) signals that component r has a significant contribution, while a small value implies a minor role.
Table 2 provides comparative analysis of our method against existing approaches. Our model achieves signifi- cantly lower MSE on the interpolation task, outperforming previous baselines by a substantial margin in addition to be- ing computationally more efficient. Additionally, Figure 6 in Appendix H presents conditionally generated samples for three test instances, illustrating the ability of our model to capture smooth temporal transitions in the latent space.
5. Experiments
We showcase the effectiveness of our model in temporal interpolation and long-term future prediction through exper- iments on common benchmarks that include synthetic and real-world longitudinal healthcare datasets. To evaluate our approach, we compare it against state-of-the-art GP prior VAEs, including SVGP-VAE (Jazbec et al., 2021), which improves the scalability of GP prior VAEs using inducing points, LVAE (Ramchandran et al., 2021), which uses in- ducing points and is specifically designed for longitudinal datasets, as well as SGP-BAE (Tran et al., 2023), which also uses inducing points but treats all parameters in fully Bayesian manner. We also include L-NODE (Iakovlev et al., 2023) in one experiment with shared governing temporal dynamics across samples, but exclude from experiments with varying group-specific dynamics. Finally, we also included BRITS (Cao et al., 2018), which implements bidirectional RNN. In all experiments, we used similar en- coder and decoder architectures for all models. We assess predictive performance using mean squared error (MSE), reporting the average and standard deviation over five it- erations. In all experiments, we used a fixed number of basis functions, M = 10, approximate continuous ker- nels Details on hyperparameter selection and latent space modeling for each experiment are provided in Appendix F, while neural network architectures are described in Ap- pendix G. The implementation of our model is provided at https: / h///ttttit thu. com/Yigitee Bon/YigiiBalik/// DDGGGGGGG...
In Figure 8, model interpretability is highlighted by exam- ples of latent functions clearly decomposed into two addi- tive components. Normalized Sobol indices show that both components contribute nearly equally, indicating the model effectively captures individual nuances and rotational dy- namics. See Figure 9 (Appendix H) for t-SNE visualization.
5.2. Health MNIST
To better approximate real-world medical data, we introduce modifications to the MNIST digits dataset simulating a high- dimensional longitudinal dataset (Krishnan al., 2015) with missing values at each time point. In this setup, the digits '3' and '6' represent two distinct biological genders, with an equal distribution across both groups A shared aging effect is modeled by gradually shifting all instances toward the bottom right corner over time. Additionally, we assume that half of the subjects in each gender group are healthy, while the other half are unhealthy. The simulated measure- ments remain unchanged for healthy subjects. Whereas, for unhealthy individuals, rotation is applied where the amount of rotation depends on the time to disease diagnosis. To further mimic noisy real-world measurements, we introduce a random rotational jitter to each data point. Finally, 25% of each sample's 's pixels are randomly selected and marked
7

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
(a)
(id), St 0.71
12 (tiee), SI = 0.0G
(b)
0.01
1.0-
0.5
0.03
0.0
MSK
-0.5
0.02
13x
0.01
10
15
0
10
15
1.5
Taxas (fime x gender), SI 0.31
Laxea (diossdd Time x dincsse)) SI
0.15
0.00
DGBFGP
LYAE
SCP-BAE
SVGP- VAE
BRITS
1.0
Ushankte
0.5
(d) 1 (time) faicaatime x gender) fare (dtorasel Time Time disease)
0.0
(P + 12 (2(tiee)) + dineeeiinos Time disesee)
(d) f2 ((ttme)+ fill ftme gender)
-10
f2 + 12 (time)
-15
(20
10
15
15
(Pimme))
Time
Figure
Health MNIST experiment results. (a) Example learned latent functions in the Health MNIST experiment along with normalized Sobol indices averaged over L = 32 dimensions. (b) Predictive MSE of different models with 32 latent dimensions on test samples.
as missing. See Figure 10 for illustration. As result, the dataset includes Q 5 auxiliary covariates, namely age, diseaseTime, gender, and discase along with subject- specific identifier id. The dataset consists of P = 1, 300 unique subjects, each observed over ng = 20 time points.
Table 3. Predictive MSE on test subjects the Physionet dataset.
Latent dimension
Moan prodiction (all training samples)
0.804 0.0000
Mean prodiction (bmining samplas of test subjects)
0.756 + 0.0000
SVGP-VAE (Jazbos et al., 2021)
32
0.802 1.0.009
SGP-RAF (Tran et al, 2023)
32
0.781 +-0.002
For the training set, we include all time points from 1,100 subjects, along with the first five time points from the re- maining subjects. The validation and test sets each contains 100 unique subjects, with their remaining time points used for evaluating the model in the future prediction task.
BRITS (Cao et 2018)
0.732 0.000
(Ramdsandran et 2021)
0.718 0.000
DGBRGP (our work)
0.619 0.009
in the intensive care unit (ICU) over a 48-hour period. Our objective is to model repeated measurements of 36 distinct attributes (we exclude weight from the original set of 37 attributes), including key physiological variables such as body temperature and heart rate. As auxiliary covariates, we incorporate id, time, ICUtype, gender, and mortality.
Figure 3(b) shows performance on the future prediction task, where our approach outperforms others. Figure 3(b) presents an ablation removing additive components to show that the numerical value R alone does not determine expres- siveness: the design of each additive component, especially the time related components, is essential for model perfor- mance even though Figure 3(a) indicates that id is the most dominant Sobol index component. Figure 10 shows visual comparisons for two test subjects, highlighting DGBFGP's ability to capture subject-specific traits and temporal dy- namies accurately. Figure 11 in Appendix depicts the relationship between M and the predictive performance and shows that our approach achieve competitive perfor- mance even with minimal number of basis functions.
In our experiments, we utilize a subset of 3,997 patients provided in "set A" of the dataset. We randomly select 200 subjects for the validation and test sets, allocating 100 sub- jects to each. For these subjects, the first 10 time points are included in the training set to learn individual temporal structure. The remaining time points are reserved for validation and testing. As shown in Table 3, DGBFGP achieves the best performance. While previous models per- form close to the mean prediction, our approach achieves a significantly lower MSE value. These results demonstrate the effectiveness of our method in real-world forecasting tasks. Figure 13 illustrates the learned latent functions and their Sobol indices. The dominant f(ll (id) accounts for heavy individual variation, while the others yield modest but significant capacity gains. Unsurprisingyy the id term prevails by encoding each instance's 's defining traits.
5.3. Physionet
We also utilize healthcare data from the Physionet Challenge 2012 (Silva et al., 2012) to benchmark our model's long- term predictive performance in real-world setting. The dataset comprises approximately 12,000 subjects monitored
8

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
Table Conditional generation MSE on unseen character combi-
nations of the SPRITES dataset.
Ground truth
Method
Latent dimension
MSB
SYGP-VAE
SGP-BAE (Tran al., 2023)
0:0115 0. 0.0021
SVGP-VAE (Jazhec et al., 2021)
64
0. 00008 + 0.0005
LVAE (Ramchandran et 2021)
0.0062 0.0004
LVAE
DGBPGP (our work)
64
0.0024 + 0.0003
Ours
Table MSE comparison for evaluating semi-amortization and generalization across training setups.
Figure 4. Conditionally generated character images from unseen attribute combinations of the SPRITES dataset. An additional example given Appendix Figure 14.
Training sctup
Method
Rotated MNIST
Health MNIST
Prysionet
Standard
SA DGBFGP
DOBEGP"
0.009 1 0.0001
0.010 I 0.0003
0.009 2 0.0000
0.013 1 0,0005
0.619 0.009
0.763 0.007
Independent
test subjects
DOBEGP
EXTRRIPP
SA-GGBKP
0.009 + 0.0003
0,0004
0.0002
0.011 : 0.0005
0.015 0,0006
0.050
0.658 0.026
0.779
0.796 0.007
5.5. Ablation of local latent variable amortization
5.4. SPRITES
In this experiment, we aim to determine whether learning a global encoder for instance-specific characteristics affects overall performance positively or negatively in the standard training setup, where at least some samples of the test sub- jects are included in the training set. Moreover, we assessed the generalization abilities of SA-DGBFGP by completely excluding the training samples of the test subjects from the training set and comparing its performance to both the non-finetuned DGBFGP model (denoted as DGBFGP*) and fine-tuned DGBFGP that has L additional parameters for each new subject (Appendix B.1).
The SPRITES dataset (Yingzhen& & Mandt, 2018) comprises 1, 1,296 unique characters, each defined by four categorical attributes: skin color, hairstyle, top clothing, and bottom clothing, with six possible options per attribute. Addition- ally, each character can perform actions from three cate- gories, each executed in three directions over T 8 time points. This results a total of N = 1, 296 9 8 images, each of size D= 64 64 x 3. The auxiliary covariates used in this experiment are lime, skin, hair, lop, bottom, action, and direction.
In this experiment, we aimed to demonstrate that our model effectively captures shared attributes defining each unique subject across a population without relying on explicit sub- ject IDs. To evaluate this capability, we reserve 296 charac- ters and all their action-time combinations for the test set, while another 296 characters are allocated the validation set. The remaining characters are used for training.
Table 5 displays the MSE values across various datasets. Under our standard training setup, SA-DGBFGP performs similarly to DGBFGP", with only a slight performance decrease observed in the simulated datasets. However, this trend does not hold for the Physionet dataset, where semi- amortization exhibits poor performance. This discrepancy may be attributed to the amortization gap. When it comes to generalization abilities, DGBFGP" cannot perform as well as SA-DGBFGP as expected. However, our proposed fine- tuned DGBFGP outperforms the cemi-amortized model.
Our method achieves noticeably lower MSE than previous approaches, as shown in Table 4. It also generates the most accurate SPRITES characters, as illustrated in Figure 4. These results highlight DGBFGP's ability to capture shared attributes and conditionally generate correlated character trajectories with unseen attribute combinations. Figure 15 shows example latent functions with normalized Sobol in- dices for each component. Unlike previous experiments, action and direction components dominate, highlighting that the latent representation encodes variability tied to those covariates. Hairstyle, clothing, and skin color contribute moderately, while time alone has minimal influence.
6. Conclusions
In this work, we introduce scalable basis-function ap- proximation for GP-prior VAEs that eliminates the usual cubic-time bottleneck. By adopting a global parameteri- zation, our method avoids explicit kernels, runs in linear time, eliminates the amortization gap, overcomes inducing point optimization constraints, supports standard mini-batch training, treats kernel hyperparameters probabilistically, and enhances latent space interpretability via Sobol indices. On various datasets, our approach yields markedly better pre- dictive performance in conditional generation tasks than current state-of-the-an models. We believe this technique will prove invaluable for applications involving correlated samples where accuracy and efficiency are essential.
Table 10 shows that DGBFGP is the fastest model in all experiments among relevant approaches. This demonstrates its improved computational efficiency in addition to more accurate conditional generation. The results were obtained by running models on a single NVIDIA Tesla V100 GPU with 32GB of memory.

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
Acknowledgements
Fortuin, V., Baranchuk, D., Rätsch, G., and Mandt, S. Gp- Deep probabilistic time series imputation. In Interna- tional conference artificial intelligence and statistics, pp. 1651-1661. PMLR, 2020.
We would like to acknowledge the computational resources provided by the Aalto Science-IT. We would like to thank the anonymous reviewers for their thoughtful comments that strengthened this work. This work supported by the Research Council of Finland (decision number: 359135)
Garrido-Mereee E. C. and Hemándeez-oobto D. Dealing with categorical and integer- variables in bayesian optimization with gaussian processes. Newrocomputing, 380:20-35, 2020.
Impact Statement
This paper introduces scalable approach for analyzing high-dimensional time-series data, designed improve computational efficiency, accuracy, and applicability for large-scale correlated datasets. We have not identified any potential negative societal impacts associated with this work.
Hensman, J., Fusi, N., and Lawrence, N. D. Gaussian processes for big data. arXiv preprint arXiv::309.0000 2013.
Hoffman, M. D., Blei, D. M., Wang, C., and Paisley, J. Stochastic variational inference. the Journal of machine Learning research, 14())1300311377, 2013.
References
lakovlev, V., Yildiz, C., Heinonen, M., and Lähdesmäki, H. Latent neural ODEs with sparse bayesian multiple shooting. In The Eleventh International Conference on Learning Representations, 2023.
Akhiezer, N. I. and Glazman, I. M. Theory of linear opera- tors in Hilbert space. Courier Corporation, 2013.
Alvarez, M. A., Rosasco, L., and Lawrence, N. D. Kernels for vector-valued functions: A review. Foundations and Trends@ in Machine Learning, 4(3):195-266, 2012.
Jazbec, M., Ashman, M., Fortuin, V., Pearce, M., Mandt, S.. and Rätsch, G. Scalable gaussian process variational autoencoders. In International Conference on Artificial Intelligence and Statistics, pp. 3511-3519, PMLR, 2021.
Bauer, M., Van der Wilk, M., and Rasmussen, C. E. Under- standing probabilistic sparse gaussian process approxima- tions. Advances neural information processing systems, 29, 2016.
Kingma, D. P. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6988, 2014.
Kingma, D. P. and Welling, M. Auto-encoding variational bayes. In 2nd International Conference on Learning Representations, 2014.
Burda, Y., Grosse, R., and Salakhutdinov, R. Importance weighted autoencoders. In International Conference on Learning Representations, 2016.
Kingma, D. P., Salimans, T., Jozefowicz, R., Chen, X.,
Cao, W., Wang, D., Li, J., Zhou, H., Li, L., and Li, Y. Brits: Bidirectional recurrent imputation for time series. Advances in neural information processing systems, 31, 2018.
Sutskever, I., and Welling, M. Improved variational infer- ence with inverse autoregressive flow. Advances in neural information processing systems, 29, 2016.
Krishnan, R. G., Shalit, U., and Sontag, D. Deep kalman filters. arXiv preprint 1r 2015.
Casale, F. P., Dalea, A., Saglietti, L., Listgarten, J., and Fusi, N. Gaussian process prior variational autoencoders. Advances in neural information processing systems, 31, 2018.
Laird, N. M. and Ware, J. H. Random-effects models for longitudinal data. Biometrics, pp. 934.974, 1982.
Cremer, C., X., and Duvenaud, D. Inference subop- timality in variational autoencoders. In International conference on machine learning, pp. 1078-1086. PMLR, 2018.
LeCun, Y., Bottou, L., Bengio, Y., and Haffner, Gradient- based learning applied document recognition. Proceed- ings of the IEEE, 86(1 12:8227--222 2 1998.
Lu, X., Boukouvalas, A., and Hensman, J. Additive gaus- sian processes revisited. In International conference on machine learning, pp. 14358-14383. PMLR, 2022.
Diggle, P. J. Analysis of longitudinal data. Oxford Univer- sity Press, 2002.
Dinh, L., Sohl- Dicksstttn,, J., and Bengio, S. Density esti- mation using real NVP. In International Conference on Learning Representations, 2017.
Margossian, C. C. and Blei, D. M. Amortized variational inference: When and why? In Uncertainty in Artificial Intelligence, pp. 2434 2449, PMLR, 2024.
10

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
Ong, P., Haussmann, M., and Lähdesmäki, H. Learning high-dimensional mixed models via amortized variational inference. In ICML 2024 Workshop Structured Proba- bilistic Inference & Generative Modeling, 2024.
Sobol, I. M. Sensitivity estimates for nonlinear mathemati- cal models. Math. Model. Comput. Exp., 1(4):407-414, 1993.
Sohn, K., Lee, H., and Yan, X. Learning structured output representation using deep conditional generative models. Advances neural information processing systems, 28, 2015.
Paszke, A., Gross, S., Massa, F,, Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al. Pytorch: An imperative style, performmree deep learning library. Advances in neural information processing systems, 32, 2019.
Solin, A. and Särkkä, S. Hilbert space methods for reduced- rank gaussian process regression. Statistics and Comput- ing, 30(2):419-446, 2020.
Quinonero-Candela, J. and Rasmussen, C. E. A unify- ing view of sparse approximate gaussian process regres- sion. The Journal of Machine Learning Research, 6: 1939-1959, 2005.
Timonen, J. and Lähdesmäki, H. Scalable mixed-domain gaussian process modeling and model reduction for lon- gitudinal data. arXiv preprint arXiv:2111.00 11.02019, 2024.
Ramchandran, S., Tikhonov, G., Kujanpää, K., Koskinen, M., and Lähdesmäki, H. Longitudinal variational autoen- coder. In International Conference Arlificial Intelli- gence and Statistics, pp.
Timonen, J., Mannerström, H., Vehtari, A., and Lähdesmäki, H. Igpr: an interpretable non- con-parametric method for in- ferring covariate effects from longitudinal data, Bioinfor- matics,
PMLR, 2021.
100::186681111 2021.
Rasmussen, C. E. and Williams, C. K. Gaussian processes for machine learning, volume MIT press Cambridge, MA, 2006.
Titsias, M. Variational learning of inducing variables in sparse gaussian processes. In Artificial intelligence and statistics, pp. 567-574. PMLR, 2009.
Rezende, D. and Mohamed, S. Variational inference with normalizing flows. In International conference ma- chine learning, pp. 1530-1538. PMLR, 2015.
Tran, B.-H., Shahbaba, B., Mandt, S., and Filippone, M. Fully bayesian autoencoders with latent sparse gaus- sian processes. In International Conference on Machine Learning, pp. 04409-34430, PMLR, 2023.
Rezende, D. J., Mohamed, S., and Wierstra, D. Stochastic backpropagation and approximate inference in deep gen- erative models. In International conference on machine learning, pp. 1278-1286. PMLR, 2014.
Yingzhen, I.. and Mandt, S. Disentangled sequential autoen- coder. In International Conference on Machine Learning, pp. 5670-5679, PMLR, 2018.
Mayo,,,-Mayol, G., Bürkner, P.-C., Andersen, M. R., Solin, A., and Vehtari, A. Practical hilbert space approximate bayesian gaussian processes for probabilistic program- ming. Statistics and Computing, 33(1):17, 2023.
Rubanova, Y., Chen, R. T., and Duvenaud, D. K. Latent ordinary differential equations for pergularly-samppl time series. Advances in neural information processing systems, 32, 2019.
Saves, P., Diouane, Y., Bartoli, N., Lefebvre, T., and Morlier, A mixed-categorical correlation kemel for gaussian process. Neurocomputing, 550: 150:126440 12472, 2023.
Silva, I., Moody, G., Scott, J., Celi, L. A., and Mark, R. G. Predicting in-hospital mortality of icu patients: The obysionet/computing in cardiology challenge 2012. In 2012 computing in cardiology, pp. 245-248. IEEE, 2012.
Snoek, J., Swersky, K., Zemel, R., and Adams, R. Input warping for bayesian optimization of non-stationary func- tions. In International conference on machine learning, pp. 1674-1682. PMLR, 2014.

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
A. Other relevant background information
Linear Models. For simplicity and interpretability, one can assume that latent variables can be modeled as linear combination of auxiliary covariates æ by standard linear model:
a1 + aQuQ + Aar 6,
where au E RR, A [ai..................... aq] R lxQ and e is noise that comes from a normal distribution with equal variance
across dimensions.
is possible to extend the capabilities of linear models and allow them to model nonlinear relationships by utilizing nonlinear basis functions. basis function extension for latent variables using Fourier features was considered (Ong al., 2024).
B. Model details
A covariate component (r) can correspond to:
A single continuous variable. In that case, the number of basis functions is B = M())
Two variables (one continuous and one categorical with ( many possible values). In that case, the number of basis functions is B(r) = C())
A single categorical variable. In that case, the number of basis functions B C(r), which our work is exclusively used for ID covariate such that C(r) = P.
The associated basis functions are represented as (xx(1)) E BB0 The corresponding weight matrix A(r) is
a,1
Ar)
[anw]] LxB())
a,
where 1 and represent the 1hh latent dimension and 8thh basis function, respectively. The collection of weights of all components can represented as
A =
[4 [40)
A [aī ... To LXQ
Here, ===ER B(r) indicates the size of the auxiliary covariates after applying basis functions.
B.1. Modeling choice for IDs
In our experiments, we did not model interactions between IDs and a continuous variable due to scalability concerns, as the number basis functions would grow to M P, significantly increasing computational complexity. Instead, we modeled IDs separately using diagonal categorical kernel (Appendix Section D.2). Following the model presented Section 4.3, the basis functions then correspond one-hot encodings, and the linear model weights sampled from the standard multivariate normal distribution. This approach simplifies the model while maintaining a manageable number of weights.
A trained model cannot be applied directly to completely new instance with an unseen ID. However, retraining the entire model from scratch is unnecessary. Since covariates other than the ID are shared across all instances, we only need to fine-tune additional L parameters to capture instance-specific characteristics for the new ID while keeping all other parameters fixed.
12

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
C. Derivation of ELBO
Let V = (A, or e) be the generic notation for unobserved random variables and -
We adopted a mean-field approximation and employed the following factorizations for each variational distribution
B N (ant Pribs o
((4)
IIII (an))
IIIIIII r=11=1
(8)
R
( II
Lognormal Or (rro, (ros)
(9)
r=1
(e) II
Lognormal ( Prtg: dite)
(10)
The derivation of ELBO lo
Y X) 2 L X,Y) is given as follows
Ps (Y
V Vp
C XX,Y)
((V)) log
dv
(
El(wy/logge Po (Y
L X]] E (ogp())) Bg g( (
Bo(a)llog
A, X)]
log p(A o. (AA)
+ (B(())(((1(4))
-E(or (ogg)log9(00
- Egce) logg {())
EG(A)
lug AA.)
E (L(g(A) p(A o, R)) (L(g(o) l| - KL(g(() (|)(e))
where, in order to estimate Ey(o)g(e) (g(A) |l p(A [)),, we sample or and & from the corresponding distributions ((rr) and q()) and denote these samples as ā and t. These values are then substituted into spectral density computation, which denote a. Assuming single Monte Carlo sample and using the mean-fiele assumptions from Equations (8)-(10)
as well as closed-form expressions for KI. divergencies between normal and log-normal distributions, the ELBO can be written as follows:
L B
( X,Y) R Eq(A)
1 (vn 1A,
ΣΣΣ
KL(g(arw) (1 plp(anolor....)
n
B Σ
(L((((a)) (()()))
KL((g(e+) (l)())
N
Ey(A)
Σ
log (yn A, 2,)
R
L
(
2 + 2
rrl
log
(11)
28,b
2
V
6=1
R
+
n
+
12g -
log Trag
log Orly
(12)
2
2
13

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
D. Kernels
D.I. Squared exponential kernel
The SE kernel, also known as the Gaussian kernel and the exponentiated quadratic kernel, is a stationary kernel used for continuous variables. It measures similarity between points based on their Euclidean distance, assuming the modeled function is smooth. It is formulated follows:
| r'11
Rse (z, x' 0,, 0² exp(
)
2
where x² the variance determining the magnitude of the variations and is the length scale determining how quickly correlation decays with distance. The spectral density of the SE kernel is given by
lw2
Sec (w [) = 22 2T exp(
2
D.2. Diagonal categorical kernel
The diagonal kernel a simple kernel used for categorical variables. measures similarity by checking whether two inputs belong to the same category. This kernel formulated as:
Rea (x, x')
{
z
otherwise.
This kernel assigns a similarity of to inputs from the same category and 0 otherwise. It is used when categories are independent and no shared structure is assumed between them.
D.3. Zero-sum categorical kernel
The zero-sum kernel also used for categorical variables. It ensures that the sum of similarities across all categories is zero. is defined as:
i x'
kiza (x,, x') -
otherwise,
where the total number of categories. This kernel assigns a similarity of for the same category and distributes the remaining similarity equally among all other categories as -ーー is suitable for modeling categorical variables where some correlation between different categories is expected while maintaining a sum-to-zero property.
D.4. Other kernels
While our experiments were carried out using the standard kernels described above, it is straightforward extend the proposed DGBFGP model incorporate other kernels as well. The only requirements are that the covariance function for each continuous covariate is stationary, such that the eigenfunction approximation is valid and the spectral density exists, and that the covariance function for each categorical covariate is symmetrie, such that the matrix eigendecomposition exists. If dataset contains non-stationarities, it is still possible to apply the proposed DGBFGP model by using e.g. stationary kernels with input warpings (Snoek et al., 2014).
E. Semi-amortized DGBEGP
When new observation arises, we typically do not add additional categories for variables like biological sex. This rule does not apply to covariates that uniquely identify individuals, such ID numbers, each new individual necessarily brings a brand-new category value. In such case, we need to re-train the whole model from scratch preferably, optimize additional L parameters as explained in Section B.I. Another approach to tackle the generalization problem is learning a global function that predicts instance-specific latent variables using AVI while keeping global parameterizations for other covariates. The generative model remains exactly the same model as in DGBFGP. The only change is the variational inference of the instafce-specific part of parameters A. Without loss of generality, we assume that the Rth component
14

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
S
0
yr
Figure 5. Graphical model of SA-DGBFGP.
in DGBFGP is the component responsible of modeling instance- characteristies. Accordingly, this component is removed from the global parameters and denoted as Z = (=1 ..... =²) in SA-DGBFGP to be treated as a local, amortized, instance-specific parameter. With slight misuse of notation, we use A to denote (A(1), ..... ( ((R11) below. BA-DGBFGP whose plate diagram given in Figure 5 is expressed as follows:
~ Lognormal(O, 1)
C, Lognormal 1)
So (x) diag (sse (V Am or, (
R-1 L
A Alo,e-
II
N(an 0, (o,,()),
r=1 1=1
Z II N(2 0.)
N
Y
II (y A, zn, ).,).,
m=
We assume factorizing variational posterior, q( Z, a. e 11) ((A)xs)) Yi)g(o))(()) where Y = pp We adopt the same mean-field assumptions given in Equations (8)-(10) for global parameters, and for local latent variables, we make use of the sample belonging to the first time point of an instance and assume the following factorization:
P II Nix" o3 o
(Z 1Z Y1
P=1
Therefore, the ELBO objective of SA-DGBF-O becomes:
logp(Y x) >
E ((A)((Z vf) [log A,2", ,,2",,/)]]
KL(go(>" 1 11
- Eete [ p(A| ar, e))]
KL(g(o)) |l p(o))
KL(((E) (()).
We note that using the standard normal distribution the prior for local latent variables is analogous predicting global linear parameters that affect specific instance based on its first sample. This is because, when using the diagonal categorical
We also experimented with utilizing all observed time points, but the results were better than those obtained using only the first time point.
15

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
kernel defined in Section D.2, both the D and @ matrices are identity matrices of size C, Ic. which activate only the et entry that is related to the category in global parameter vector, Hence, can plug inferred latent variables into additive structure without any concern.
F. Experimental details
Since we specified a general Lognormal 1) prior for the lengthscale parameters, we standardized all continuous auxiliary covariates have zero mean and unit variance exclusively for training DGBFGP. contrast, other GP prior VAEs do not impose prior kernel hyperparameters, as these optimized jointly with the neural network parameters (Jazbec et al., 2021; Ramchandran al., 2021). Regardless of the approach, samples scaled to the range [0, 1] in the Rotated MNIST, Health MNIST, and SPRITES experiments, while they are standardized in the Physionet experiment.
In all experiments, we used a fixed number of basis fuunction, M = 10, to approximate continuous kernels. Additionally, for all experiments, the boundary condition ! = [-J, J] was defined such that / is approximately c times half the range the corresponding continuous values, following the recommendations in (Riutort-Mayol et al., 2023; Timonen & I Lähdesmäki, 2024). The approximation with a finite number of basis functions is expected to be effective as long as J is not too close the absolute maximum of the continuous covariate (Riutort-Mayol et al., 2023). Based on this criterion, set 1.5.
We implemented our model using PyTorch (Paszke et al., 2019). Across all experiments and methods, we employed the Adam optimizer (Kingma, 2014) with learning rate of 0.001 and a batch size of 256. To ensure stable convergence, the learning rate was dynamically reduced whenever the validation loss did not improve for a specified number of epochs. Specifically, the learning rate was scaled down by a factor of 0.5 no improvement was observed for 10 consecutive epochs. Additionally, early stopping was employed to prevent overfitting, and the final model for each experiment was selected based on the best performance on the respective validation set.
In the remainder of this section, provide detailed specifications for each experiment. We denote the additive components utilizing different covariance functions as follows: f-e(.) for the SE kernel, fsexca(.) for the interaction of SE and CA kernels, and (.a for the CA kernel which is applied exclusively to the id covariate, when applicable.
F.1 Rotated MNIST
For the validation and test sets, we randomly selected 80 distinct instances and randomly sampled four consecutive angles for every instance. The remaining samples were allocated to the training set. As a result, the training set comprises Ntrain = 5, 760 samples, while both the validation and test sets contain Nyal = Neest. = 320 samples each.
For the Rotated MNIST experiment, we modeled the latent space using the following additive kernel structure:
f = f= + ((ngle)))
In contrast, LVAE (Ramchandran al., 2021) employed a more complex setup that incorporates an additional component that models interactions between id and the continuous angle covariate, (angea id), as this configuration was reported to be the most effective in the original work (Ramchandran et al., 2021).
The maximum number of epochs for this experiment was set to 2, 000.
F.2. Health MNIST
The kernel structure for the Health MNIST experiment defined as:
f f (id) + (22) (time) + f(3) (time x gender) + fl (diseaseTime x discase).
Similarly, as in the Rotated MNIST experiment, LVAE used an additional interaction component frexca x id). The maximum number of epochs for this experiment was also set to 2, 000.
16

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
F.3. Physionet
For the Physionet experiment, our model utilized the following kernel structure:
-= / (time) Ssexea (time ICUtype)
+ - (time x gender) - Jex (time x mortality),
whereas LVAE included an additional interaction term, fsexca (time x id).
The maximum number of epochs for this experiment was set to 1, 000.
F.4. SPRITES
For the SPRITES experiment, both our model and LVAE used the same kernel structure:
f = f()) (time) f (time skin) fsexea (time bottom)
(() sexea (time top) + (()) (time x hair)
S (time x action) + sexca (time x direction).
(()
The maximum number of epochs for this experiment was set to 100.
G. Neural network architectures
We present neural network architectures employed across different experiments. Table 6 details the convolutional neural network (CNN) architecture used for the Rotated MNIST experiment; it is similar to the one in (Casale et al., 2018), with minimal modifications. Likewise, Table outlines the CNN architecture used for Health MNIST, closely following (Ramchandran et al., 2021). Table 8 describes the multi-layer perceptron (MLP) architecture employed for the Physionet experiment, and Table 9 shows the CNN-based architecture from (Jazbec et al., 2021) with a minor alteration in the final activation layer. CNN-based architectures, we use a sigmoid activation function in the last layer the generative network, which is better suited to our pre-processed data. By contrast, for the MLP-based architecture, we omit an activation function in the final layer because the processed data are bounded.
Table 6. Neural network architecture used in Rotated MNIST experiment
Hyperparameter
Value
Dimensionality of input
28 x 28
Number of convolution layers
3
Number of filters per convolution layer
72
Kernel size
Inference network
Stride
Number of feedforward layers
Width of feedforward layers
128
Dimensionality of latent space
Activation function of layers
ELU
Dimensionality of input
Number of transposed convolution layers Number of filters per transposed convolution layer
72
Kernel size Stride
3
Generative network
Number of feedforward layers
Width of feedforward layers
128
Activation function of layers
ELU, Sigmoid
17

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
Table 7. Neural network architecture used in Health MNIST experiment
Hyperparameter
Value
Dimensionality of input
36 x 36 2
Number of convolution layers
Number of filters per convolution layer
16, 32
Kernel size
Stride
Pooling
Max pooling
Inference network
Pooling kernel size
2 x 2
Pooling stride
Number feedforward layers Width feedforward layers Dimensionality of latent space
300, 30
Activation function of layers
RELU
Dimensionality of input
Number of transposed convolution layers Number of filters per transposed convolution layer Kernel size Stride
2
16
x
Generative network
Number of feedfforard layers
Width of feedforward layers
30, 300
Activation function layers
RELU, Sigmoid
Table 8, Neural network architecture used in Physionet experiment
Hyperparameter
Value
Dimensionality of input
Number of feedforward layers
2
Inference network
Width of feedforward layers
300, 30
Dimensionality of latent space
L
Activation function layers
RELU
Dimensionality of input Number feedforward layers
L
Generative network
Width of feedforward layers Activation function of layers
30, 300
RELU
Table Neural network architecture used SPRITES experiment
Hyperparameter
Value
Dimensionality of input
3 64 x 64
Number convolution layers
Number filters per convolution layer
Inference network
Kernel size Stride
3 x 3
1,2
Dimensionality of latent space
L
Activation function layers
ELU
Dimensionality of input Number convolution layers Number of filters per convolution layer
L
16
Generative network
Kernel size Stride Activation function layers
3 3
ELU, Sigmoid
18

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
H. Additional results and comparisons
Ground truth
3
3
3
9
3
3
3
3
3
3
3
3
3
3
3
3
3
-VG
3
3
3
3
3
3
3
LVAB
3
3
3
3
3
3
SGP-BAE
3
3
3
3
3
3
3
3
L-NODB
3
3
3
3
3
3
3
Ours
Figure Conditionally generated images three different test instances unseen rotation angles in Rotated MNIST experiment.
0.07
0.06
0.05
MSE
0.04
0.03
0.02
2
10
Figure 7. Test MSE as a function of in the Rotated MNIST experiment.
19

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
iid). SI = 0.48
fole) aaggle), SI = 0.52
1.5
ID: 294
1.0
ID: 107
ID: 399
0.5
0.0
-0.5
-1.0
-1.5
0
2
4
6
0
2
4
Angle
Figure 8. Example learned latent functions in the Rotated MNIST experiment along with normalized Sobal indices averaged J. = 16 dimensions.
Table 10. Run time comparison of GP prior VAEs per epoch in terms of seconds.
Method
Rotated MNIST
Health MNIST
Physionet
SPRITES
LVAE (Ramchandran et al., 2021) SVGP-VAE (Jazbec et al., 2021) DGBFGP (our work)
7.7 + 0.1
34.9 0.4
85.2 0.9
88.3 + 0.5
1.8 + 0.1
13.6 + 0.5
80.7 0.5
35.9 + 0.6
0.5 + 0.0
12.5 + 0.2
13.1 0.1
10.2 + 0.4
20

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
Angle
0.0
0.42
0.84
1.26
1.68
2.09
2.51
2.93
3.35
3.77
4.19
4.61
5.03
5.45
5.86
6.28
Figure t-SNE embedding latent variables the Rotated MNIST data.
21

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
Ground truth
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
SVGP-VAE
3
3
3
6
3
3
3
3
3
3
3
3
SGP-BAE
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
LVAE
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
Ours
3
3
3
3
3
3
3
3
3
3
3
3
Ground
6
6
G
6
6
6
6
SVGP-VAE
6
6
6
6
6
6
6
6
6
6
6
SGP-BAE
6
6
6
6
6
6
6
6
6
6
6
LVAE
6
6
6
6
6
6
6
6
6
6
6
Ours
6
6
6
6
Figure 10. Conditionally generated images of two different test instances for future prediction on the Health MNIST dataset.
0,030-
0.025
MSE
0.020
0.015
0.010
2
M
Figure 11. Test MSE as function M the Health MNIST experiment.
22

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
Gender 0
Gender
All Genders
10.0
75
5.0
-25
-5.
-15
Healthy
Unhealthy
Figure t-SNE embedding the latent variables of the Health MNIST data, color bar representing the diseaseTiine covariate. This covariate is applicable exclusively to unhealthy subjects.
J0a (id), SI 0.85
(time), 0.02
Sexxn (time ICUtype), 0.07
ICU1
100
0.5
10: 1014
2070
1D: 2030
0.5
2015
ID: 2051
1.0
0
20
40
(time x gender), 0.02
Jexxon (time x mortality), 0.04
Male
Douth
0.5
0.0
1..
20
40
0
20
40
Time
Figure 13.. Example learned latent functions in the Physionet experiment along with normalized Sobol indices averaged over L 32 dimensions.
23

---

Bayesian Basis Function Approximation for Scalable Gaussian Process Priors in Deep Generative Models
Ground
truth
SVGP-VAE
LVAE
Ours
Figure 14. Additional example of generated character images from unseen combinations of the SPRITES dataset.
(time), SI 0.01
(time (time skin), SI 0.04
(2)
(33
(time x bottom), SI 0.04
1.5
Jancen
Ainn
Suttomm
buttom
siie
skins
Jolicon
bott.o0m
0.5
-1.0
1.55
0
2
2
4
6
Aexcs (time x top), ST = 0.11
Texeaa (time x hair), SI 0.18
Jc xea (time x action), SI = 0.30
(6)
1.5
halr
halrs
1.0
lnitg
Lnits
spell
t00c 100g
0.5-
1.0-
-155
2
4
6
poxcn (time x direction). SI
(7)
1.5
1.0-
righn
0.5
0.0
0.5
-1.
-1.5
2
4
6
Time
Figure Example learned latent functions in the SPRITES experiment along with normalized Sobol indices averaged over 64 dimensions.
24


Reviews:
[{"review_id": "UOnHeyxxKb", "reviewer": "Reviewer_Fbnr", "summary": "This work presents a scalable basis function-based approximation for Gaussian Process prior Variational Auto-Encoders (GP-VAEs), to overcome the cubic time-complexity (without resorting to inducing-point GP variational inference techniques) and to accomodate shared and individual-specific correlations across time.\n\nTheir method allows for continuous and categorical covariate information to be incorporated for conditional generation, due to the proposed Hilbert space kernel approximation based on the kernels' eigenvalue and eigenfunction decomposition (described in Section 4). More precisely, the authors propose an additive GP prior for VAEs, which is defined using such Hilbert space eigen-decomposition of kernels.\n\nDue to the proposed decomposition, learning can be posed as a variational inference over global parameters (kernel hyperparameters and linear model parameters $A$), with run-time complexity that scales linearly in the size of the dataset (amenable to mini-batching).\n\nResults are reported on synthetic and real-world datasets, showcasing good predictive performance (as measured by Mean Squared Error), at reduced computational complexity.\n\n# After rebuttal\nThe authors provided informative clarifications and additional Sobol Index based analysis of their results for further illustration. I hence lean towards acceptance of the work, with revised final manuscript.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Independence across latent dims and additive GP prior** Correlated GP priors are typically formulated via the linear model of co-regionalization (LMC), which multiplies independent GPs by a factor loading matrix to introduce correlations across latent dimensions. In GP prior VAE models, a neural network–parameterized decoder maps latent variables to likelihood parameters, automatically introducing correlations at least as expressive as those from LMC. This approach is standard in previous GP prior VAE models, so assuming independence across dimensions does not lose generality or cause any limitation.\n\nUsing additive GPs allows the latent effects to be decomposed into additive terms while leveraging a scalable basis function formulation. Although non-additive kernels may be more expressive in theory, our results show that DGBFGP significantly outperforms previous methods—even some that employ non-additive kernels.\n\n**Answer for Q1.** The Hilbert space approximation for GPs has a key property: the eigendecomposition of the Laplace operator with Dirichlet conditions is available in closed form and is independent of the kernel [1]. Similarly, a stationary kernel’s Hilbert space approximation depends on hyperparameters only through its spectral density, which is known for common kernels like SE and Matern. Thus, eigen-decomposition computation is unnecessary; the eigenfunctions and eigenvalues can be used in a plug-and-play fashion with an overall cost of $O(1)$, even when hyperparameters vary. We treat hyperparameters as random variables (Section 4.2) and infer them using VI (Section 4.4).\n\nFor a categorical covariate, we compute the eigendecomposition of a $C$ x $C$ matrix ($C$ is the number of categories), and for the discrete kernels used here, closed-form solution again yields a cost of $O(1)$. Moreover, for arbitrary categorical kernels, with varying or fixed parameters, computational overhead is negligible since $C$ is typically much smaller than $N$.\n\nWe will clarify these details in the revised manuscript.\n\n[1] https://link.springer.com/article/10.1007/s11222-019-09886-w\n\n**Answer for Q2.** We focused on the SE kernel because it is probably the most commonly used and it also worked well in our experiments. In practice, the choice of the kernel depends on the application. As noted at the end of Section 4.1, the Hilbert space approximation applies to any stationary continuous covariance function with a known spectral density, including the common SE, Matern, and many other kernels found in the literature.\n\n**Answer for Q3.** The Hilbert space approximation applies to additive, product, and sums of product kernels (as in our work). For instance, an SE kernel that depends on all input covariates—with a distinct length scale per covariate is equivalent to a product of SE kernels, leading to a product of basis function approximations (as in Section 4.3). Here, the number of terms scales as $M^q$, with $M$ basis functions per covariate and $q$ covariates. This is computationally feasible for a small $q$ but not scalable in general. However, Eq. (4) still factorizes as $A \\mid \\boldsymbol{\\sigma}, \\boldsymbol{\\ell} \\sim \\prod_{l = 1}^L \\mathcal{N}(\\boldsymbol{a}_l \\mid \\boldsymbol{0}, S(\\sigma, \\ell))$, where $S(\\sigma, \\ell)$ is again diagonal, meaning the Gaussian distributions become $M^q$ dimensional. A similar factorization holds even for dependent GP priors before applying LMC (see our response above), in which case modeling a separate factor loading matrix would introduce latent correlations (or would be incorporated into the decoder function as described above).\n\n**Answer for Q4.** Since assessing latent function *quality* is challenging, we propose to quantify the interpretability by measuring contributions of each additive component using the Sobol index, defined for component $r$ as $\\frac{\\mathrm{Var}[f^{(r)}(x^{(r)})]}{\\mathrm{Var}[\\sum_rf^{(r)}(x^{(r)})]}$ . The Sobol index values for the four datasets are shown below. \n\n|Dataset|Component|%|\n|-|-|-|\n|Rotated MNIST|id|48|\n||rotation|52|\n|Health MNIST|age|3|\n||id|71|\n||age x gender|11|\n||diseaseAge x disease|15|\n|Physionet|id|85|\n||time|2|\n||time x ICU|7|\n||time x gender|2|\n||time x mortality|4|\n|Sprites|time|1|\n||time x body|4|\n||time x bottom|4|\n||time x top|11|\n||time x hair|18|\n||time x action|30|\n||time x direction|32|\n\nThe id component’s dominance is expected since it captures each instance’s primary distinguishing features and serves as a natural baseline. Other factors add nuance and boost model capacity.\n\nFor visual interpretability, the latent functions in Figure 1 (from the Health MNIST experiment in Section 5.2) highlight the model’s interpretability. We will include additional visualizations and expand the discussion in the revised Appendix.\n\nFinally, the sensitivity to $M$ is reported in Appendix G (Figures 5 and 7), showing that DGBFGP achieves the best results with $M = 6$ for Rotated MNIST and $M = 4$ for Health MNIST.", "claims_and_evidence": "The claims are generally well-supported.\n\nThe core idea of using kernel eigen-decomposition for an scalable GP approximation is well-grounded in already established theory.\n\nThe paper effectively demonstrates the practical benefits of this approach through empirical evaluation.", "methods_and_evaluation_criteria": "The proposed method, i.e., utilizing additive GP priors and kernel eigen-decomposition, is a sound alternative for GP-VAE inference.\n\nThe evaluation, primarily focused on MSE for predictive performance, is relevant to the task of time-series prediction. \n\nHowever, expanding the evaluation to include metrics assessing the quality of the learned latent spaces would strengthen the analysis.", "theoretical_claims": "The main theoretical results are:\n\n1. The eigen-decomposition of kernels in Sections 4.1-4.3: these results are based on previously known results, so they are correct to the best of my knowledge.\n\n2. The variational ELBO of Equation 7, with details in Appendix C: the presented expression appear correct upon review.", "experimental_designs_or_analyses": "The experimental setup, using modified benchmark datasets (MNIST, Physionet, SPRITES), is appropriate for demonstrating the method's ability to capture time-varying latent forces driving multi-dimensional time-series observations.\n\nHowever, solely relying on the MSE as a metric for evaluation, limits the experimental assessment. Evaluating the quality of the learned latent time series, and its dependence over the number of eigen-functions $M$ used to approximate a kernel, would provide a more comprehensive picture.", "relation_to_broader_scientific_literature": "The key contribution of this work is to combine known results and ideas (eigen-decomposition of kernels with GP priors for VAEs) to propose a model variant that is computation efficient (linear complexity) and performant.\n\nThe authors provide a good overview and description of the main relevant literature (Sections 2, 3, 4.1-4.3), clearly explaining their novelty in combining kernel eigen-decomposition with an additive GP prior for VAEs, which can be defined via a global parametrization that implies increased computational efficiency.", "essential_references_not_discussed": "The main references to related work are well described.", "other_strengths_and_weaknesses": "Strengths:\n\n- The key strength is the effective use of kernel eigen-decomposition to achieve linear-time complexity in GP-VAE inference.\n- The ability to handle both continuous and categorical covariates is a valuable contribution.\n\nWeaknesses:\n\n- The assumption of independence across latent dimensions and the use of additive GP priors per dimension could be limiting.\n\n- The computational cost of computing the eigen-decomposition, especially for varying kernel hyperparameters, is not thoroughly discussed.\n\n- The paper mainly focuses on Squared-Exponential kernels.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "- Please provide a detailed explanation of the computational procedure for calculating the eigen-decomposition of Squared-Exponential kernels, particularly on how this process scales with varying kernel hyperparameters. Quantifying the computational cost would be very helpful.\n\n- Beyond the Squared-Exponential kernel, what other continuous kernels can be efficiently approximated using the proposed eigen-decomposition method?\n\n- Could you elaborate on the challenges and potential solutions for extending the method to non-additive GP priors, specifically addressing the complexity of handling a non-fully factorized matrix A in Equation (4)?\n\n- Would it be feasible to include an evaluation of the quality of the learned latent GP functions, particularly in synthetic experiments? Additionally, please discuss the sensitivity of the results to the choice of M, the number of eigenvalues and eigenfunctions used.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "NWYpwbdw7N", "reviewer": "Reviewer_VgcQ", "summary": "This paper proposes a generative model based on Variational Autoencoders (VAE) where latent variables are assigned a GP prior.\nThe main contribution is to approximate the GP prior with random features so that the model can be optimized through mini-batching and linearly in the number of data.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer's feedback and the opportunity to clarify our contributions. Before addressing specific comments, we emphasize that, unlike traditional methods that use random feature approximations, our approach leverages the Hilbert space approximation of the GP prior.\n\n**Interpretation of the results.** We agree that the interpretability of our parameterization is important. Although quantifying interpretability with neural networks is inherently challenging, our latent function visualizations for the Health MNIST data (see Figure 1) offer clear qualitative insights into the model's behavior. In the final version, we will include additional figures and detailed discussions to further illustrate these insights.\n\nPlease refer to our response to reviewer Fbnr (Q4) for an additional Sobol index metric that quantifies interpretability by measuring the contributions of each additive component.\n\n**Global parameterization.** We thank the reviewer for highlighting our elimination of amortized inference and the related amortization gap [1,2,3]. We will emphasize this contribution in the revised version and add experiments demonstrating the benefits of our approach over an amortized variant. In the DGBFGP model, the instance-specific additive latent components correspond to local latent variables that could be amortized. As detailed in Appendix (Section B.1), these components are modeled using a categorical kernel based on each individual’s \"id\". In the Hilbert space approximation, the instance-specific component becomes an additive offset in the latent space with a standard Gaussian prior, here denoted for simplicity as $\\boldsymbol{a_p}$ for individual $p$. An amortized version then uses $q_{\\phi}(\\boldsymbol{a}_p \\mid Y_p)$ as an amortized variational approximation of the true posterior, where $\\phi$ denotes inference network that is *shared* across all individuals. Our ablation study shows that sharing an inference network across individuals leads to poorer performance compared to our proposed parameterization:\n\n|Method|Rotated MNIST|Health MNIST|\n|-|-|-|\n|Amortized DGBFGP|$0.010 \\pm 0.0004$|$0.012 \\pm 0.0002$|\n|DGBFGP|$0.009 \\pm 0.0001$|$0.009 \\pm 0.0000$|\n\nSimulated datasets using MNIST images allow for easy and accurate amortized inference since the unrotated digit images are available for each individual. However, for more complex datasets, amortized inference would require more complex networks, and the amortization gap is likely to increase. We will include these analyses, along with studies on other datasets, in the revised version.\n\n**Contributions.** GPPVAE [4] was the first GP prior VAE model, with subsequent works building incrementally on it; [5] infers GPs for each trajectory separately, while [6] and [7] use inducing points (with [7] focusing on longitudinal designs using additive kernels). Our contributions are as follows: 1) We argue that the utilization of inducing points in the presence of categorical covariates is problematic and a non-trivial task since their locations cannot be optimized with gradients. To our knowledge, no efficient general solution has been proposed to handle discrete covariates in latent variable models, GP prior VAEs in particular, so far. Our work solves those problems by using the Hilbert space approximation method. 2) All previous works need to explicitly handle the kernels and their approximations one way or the other. Our approach provides a direct way to avoid explicit kernels. 3) We also propose to handle kernel hyperparameters probabilistically using VI (we note that was also done in Tran et al [8]). 4) Instead of performing amortized inference using a shared inference network, we directly optimize the global parameters (provided by our approximation method). This avoids the amortization gap and thereby improves the performance. We now also show and quantify the amortization gap via additional experiments. 5) Overall, we demonstrate that our model is more scalable and outperforms previously proposed methods. 6) Additionally, motivated by reviewers' comments, we now better demonstrate the interpretability of the proposed model by visualizing the latent effects as well as quantifying the contributions of different effects (see our response to Q4 of the reviewer Fbnr). \n\n**Additional competitor.** We thank the reviewer for bringing up SGP-BAE [8]. We run this model on the synthetic datasets and we will include results for the remaining experiments in the final version.\n\n|Method|Rotated MNIST|Health MNIST|\n|-|-|-|\n|SGP-BAE|$0.023 \\pm 0.0006$|$0.024 \\pm 0.0077$|\n|DGBFGP|$0.009 \\pm 0.0001$|$0.009 \\pm 0.0000$|\n\n[1] https://proceedings.mlr.press/v80/cremer18a.html\n\n[2] https://proceedings.mlr.press/v84/krishnan18a.html\n\n[3] https://proceedings.mlr.press/v80/marino18a\n\n[4] https://arxiv.org/abs/1810.11738\n\n[5] https://arxiv.org/abs/1907.04155\n\n[6] https://arxiv.org/abs/2010.13472\n\n[7] https://arxiv.org/abs/2006.09763\n\n[8] https://arxiv.org/abs/2302.04534", "claims_and_evidence": "The claim is that this approach offers better and more interpretable modeling compared to other competitors from the literature.\nThere are some experiments showing better performance but I couldn't find a thorough analysis on interpretability; I believe that this could be a major selling point of the paper and I would encourage the Authors to focus on this aspect more.", "methods_and_evaluation_criteria": "I found the evaluation appropriate and the number of benchmark datasets provided to be sufficient.", "theoretical_claims": "There are no theoretical developments in the paper.", "experimental_designs_or_analyses": "Due to the experimental nature of the paper, I think that the experimental campaign needs to be extensive.\nOverall, I think that the Authors did a good job in selecting a wide range of data sets and reporting some performance measures against competitors.\nAgain, I would encourage the Authors to focus their efforts on the intepretation of the results; this seems to be one of the main selling points of the proposed parameterization fo the GP prior and I believe it needs to be expanded on.", "relation_to_broader_scientific_literature": "I think that the paper does a good job in characterizing the literature, with some exceptions (see below).", "essential_references_not_discussed": "There is an ICML 2023 paper presenting a Bayesian autoencoder where latent variables are given a sparse GP (and Deep GP) prior, which I think could be a good competitor for this work:\n\n[1] B.-H. Tran, B. Shahbaba, S. Mandt, and M. Filippone. Fully Bayesian Autoencoders with Latent Sparse Gaussian Processes. ICML 2023.", "other_strengths_and_weaknesses": "I think that one of the main contributions is not highlighted too well in the paper, that is the one about parameterization.\nThe paper mentioned that thanks to the approximation and inference strategy there is no need for amortized inference; I believe this point could be emphasized in the main text (maybe by dedicating some space to this in Sec 4?) and in the experiments (e.g., by showing the advantage of no amortized inference vs the proposed parameterization).\n\nOverall, one weakness is that the paper risks being just another model in the GP-VAE \"zoo\"; random feature approximations for GPs and Deep GPs are rather common and this might give the impression that the paper is a straightforward combination of known elements.\nI would encourage the Authors to think of ways in which their work can be clearly differentiated from others in this literature and focus any additional experimental efforts in showing the advantages associated with their proposal.", "other_comments_or_suggestions": "Overall the paper is well written, so I don't have specific suggestions for changes in the the writing.", "questions_for_authors": "I've mentioned a few aspects which could improve the paper, and I hope that the rebuttal period will be useful for the Authors to come up with some constructive arguments in favor of acceptance.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "GZZsCLWHcU", "reviewer": "Reviewer_L1D7", "summary": "The authors present an approach for performing efficient variational inference in an additive GP prior VAE (additive GP prior with an MLP parameterised likelihood). The approach uses standard low-rank Hilbert space approximations for the GP kernels, approximating the additive GP prior as an additive linear model. The experimental results demonstrate that the proposed model and approximation outperform a number of GP-VAE baselines on several experiments.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Lacking novelty.** We kindly ask the reviewer to find what we contribute to the literature by referring to our response \"Contributions\" for the reviewer VgcQ.\n\n**Generalization for new instances.** Thank you for giving us the opportunity to clarify this important aspect of our model. We agree with the reviewer in that it is absolutely essential that predictions for new distinct test instances are adapted to each instance. However, we respectfully disagree with the reviewer's claim that our method cannot do that. As detailed in Appendix B.1, our approach generalizes to new distinct test instances by variationally optimizing $L$ instance-specific id parameters using only the data from the test instance from the test instance while keeping all other model parameters fixed at test time, where $L$ is the latent space dimensionality. Given that we need to optimize only a small number of instance-specific variational parameters, it efficiently learns instance-specific characteristics without significant computational overhead. Moreover, our approach to generalize to new distinct test instances is also theoretically more accurate than the amortization-based approach as we can avoid the amortization gap when generalizing to new test instances. \n\nTo clarify this further, we provide additional evidence that our proposed method indeed can generalize to new distinct test instances even better than methods that use amortized VI. In our experiments, the model variant where instance-specific parameters are optimized at test time corresponds to the proposed method and is denoted as DGBFGP, and that is compared to the amortized model variant as well as to a baseline model that does not have any mechanism to generalize to new test instances (denoted by DGBFGP$^*$).  \nOur proposed model outperforms both the baseline and an amortized variant, as shown below:\n\n|Method|Rotated MNIST|Health MNIST|\n|-|-|-|\n|DGBFGP$^*$|$0.071 \\pm 0.0002$|$0.039 \\pm 0.0003$|\nAmortized DGBFGP | $0.013 \\pm 0.0006$ | $0.015 \\pm 0.0005$\nDGBFGP (proposed model) | $0.009 \\pm 0.0003$ | $0.011 \\pm 0.0005$\n\nThis clearly demonstrates that our non-amortized strategy yields highly competitive performance, effectively mitigating the amortization gap. Moreover, our method retains computational efficiency by limiting the optimization to only $L$ parameters per instance.\n\n**Comparison against sparse GP approx.** Thank you for giving us the opportunity to clarify also this important aspect. The LVAE and SVGP-VAE methods both use sparse GP approximation. LVAE method also assumes the same additive GP prior structure as our proposed model, whereas SVGP-VAE is limited to object-view product kernel. In other words, LVAE and SGP-VAE models that are already included in our comparisons implement exactly the kind of comparison that the reviewer is asking. In all experiments that we conducted, we provided a comparison of our method that leverages Hilbert space approximation to those that use sparse GP approximations. We will clarify this in the revised manuscript.", "claims_and_evidence": "The main claims are supported with evidence.", "methods_and_evaluation_criteria": "Yes. The method is evaluated on MNIST-derived, SPRITES, and Physionet based experiments, which it seems are standard for papers in this area.", "theoretical_claims": "No,", "experimental_designs_or_analyses": "The experimental design follow that of previous papers.", "relation_to_broader_scientific_literature": "N/A.", "essential_references_not_discussed": "N/A.", "other_strengths_and_weaknesses": "## Strengths\n\n* The paper is clearly written and easy to follow.\n* The model and approximation method outperform the baselines across the experiments. \n\n## Weaknesses\n\n* Whilst the method is sensible and correct, I can't help but feel the paper is lacking novelty---the methods apply existing GP approximation techniques to an existing model.\n* Although the global variational parameterisation is hailed as a strength of the method; because the variational approximation is not amortised it is unable to generalise to new distinct instances in e.g. the MNIST experiments.\n* I presume that it is possible to use a sparse GP approximation for the additive GP prior VAE model? If so, it would be useful if a comparison was provided.", "other_comments_or_suggestions": "N/A.", "questions_for_authors": "N/A.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ljbg7aorV7", "reviewer": "Reviewer_tjXK", "summary": "This paper addresses the computational challenges of using Gaussian process (GP) priors in Variational Autoencoders (VAEs) for high-dimensional time series analysis. While GP-based VAEs effectively capture temporal dependencies, their cubic time complexity limits scalability. To overcome this limitation, the authors propose a scalable basis function approximation for additive GP priors in VAEs, reducing computational complexity to linear time. The authors claim that their approach is not only computationally efficient but also capable of capturing complex correlations within and across subjects. Empirical evaluations on synthetic and real-world datasets demonstrate that the proposed method enhances computational efficiency while significantly improving predictive performance compared to existing approaches.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments.\n\n**Missing relevant models.** Our experiments include previous GP-VAE methods, as our main goal is to improve these models. We also evaluate an RNN-based method (BRITS) and a latent neural ODE (L-NODE). We also compare against SGP-BAE suggested by reviewer VgcQ. We're happy to run additional comparisons if the reviewer specifies which method they refer to.\n\n**Lacking runtime comparisons.** Initially, we focused on LVAE due to its similarity to our method since only LVAE and DGBFGP can handle an arbitrary number of covariates. However, we agree that adding more baselines gives a bigger picture of DGBFGP's efficiency. We now include SVGP-VAE in our runtime comparison.\n\n|Method|R. MNIST|H. MNIST|Physionet|SPRITES|\n|-|-|-|-|-|\nLVAE|$7.7\\pm0.1$| $34.9\\pm0.4$|$85.2\\pm0.9$|$88.3\\pm0.5$|\nSVGP-VAE|$1.8\\pm0.1$|$13.6\\pm0.5$|$80.7\\pm0.5$|$35.9\\pm0.6$\nDGBFGP|$0.5\\pm0.0$|$12.5\\pm0.2$|$13.1\\pm0.1$|$10.2\\pm0.4$\n\n**VI reference is missing.** We agree. We will update the manuscript with an appropriate citation for SVI, ensuring our methodology is properly contextualized within the broader VI literature.\n\n**Stronger motivation.** Agreed. In the final version, we will expand our discussion to highlight existing models' limitations and how our model improves them. See our response \"Contributions\" for reviewer VgcQ for revised motivation detailing these limitations and our contributions.\n\n**Compare to non-transformer-based alternatives.**\nHigh-dimensional multivariate time series face challenges from complex correlations, time-varying covariates, and missing values. Standard VAEs with an iid normal prior often struggle with these aspects. By incorporating a GP prior into a VAE, our model: 1) Provides reliable uncertainty estimates without extra mechanisms (unlike state-space models or RNNs). 2) Uses tunable kernel functions for transparency, outperforming \"black-box\" approaches. 3) \"filters out\" noise by learning data generating process to reveal the underlying structure.\n\nMoreover, while ODE-based methods struggle with varying dynamics (e.g., between healthy and unhealthy subjects), our GP-based model remains robust. Experiments comparing it to BRITS and L-NODE underscore these advantages, making it especially beneficial for healthcare.\n\n**Answer for Q1.** We agree: our method's performance depends on R, but we view R and its latent components as a modeling choice rather than a hyperparameter. The expressiveness of DGBFGP stems from its additive components that model each covariate. For instance, the Health MNIST experiment in Appendix E.2 uses the following additive model:\n\n$$\nf_{\\mathrm{ca}}^{(1)}(\\mathrm{id})+f_{\\mathrm{se}}^{(2)}(\\mathrm{time})+f_{\\mathrm{se \\times ca}}^{(3)}(\\mathrm{time \\times gender})+f_{\\mathrm{se \\times ca}}^{(4)}(\\mathrm{diseaseTime \\times disease})\n$$\n\nAblation results for the Health MNIST by removing some components from the original implementation:\n\n|Additive model|R|MSE|\n|-|-|-|\n|$f^{(1)}+f^{(2)}+f^{(3)}+f^{(4)}$|4|0.009|\n|$f^{(1)}+f^{(2)}+f^{(4)}$|3|0.010|\n|$f^{(1)}+f^{(2)}+f^{(3)}$|3|0.025|\n|$f^{(1)}+f^{(2)}$|2|0.025|\n|$f^{(1)}$|1|0.031|\n|$f^{(2)}$|1|0.036|\n\nRemoving the 3rd component has little effect since the gender signal is indirectly captured by the digit type. In contrast, 4th component is crucial, as it captures the evolving disease signal; without it, the model cannot discern an instance's health trajectory. Similarly, removing time-related components weakens temporal correlation capture, and eliminating the id term forces all instances to share latent variables, thereby losing individual nuances. This analysis shows that while the numerical value of R doesn't directly dictate expressiveness, the design of each additive component is essential. We will include this analysis into the revised Appendix.\n\n**Answer for Q2.** Since the main focus of this paper is to model latent space using GP priors, for a fair comparison, we followed the architectures and latent dimensionaly choices provided in the previous GP prior VAE works as we denoted in Appendix F.\n\n**Answer for Q3 and Q4.** In models like LVAE and SVGP-VAE, inducing point placement is crucial because they summarize the latent function over the input space. If they don't cover all covariates, key variations can be missed, leading to suboptimal performance, especially with discrete covariates, whose locations can't be optimized with gradients.\n\nFor fair comparisons, we used identical architectures across models, as our goal is to present an alternative GP prior approximation for latent space modeling. The Hilbert space approximation offers global parameterization, eliminating the need for a shared inference network and the associated amortization gap. As [1] shows that the amortization gap is unavoidable in latent variable models, including GPs, our method avoids the suboptimality of a shared encoder network.\n\n[1] Amortized Variational Inference: When and Why?\nhttps://arxiv.org/abs/2307.11018", "claims_and_evidence": "The main claim of the paper on its proposed method include\n\n1) Model Expressiveness: The proposed model, DGBFGP, is claimed to be as expressive as or more expressive than other GP-based VAEs in capturing complex correlations within and across subjects.\n\nThe authors support this claim by evaluating the model on both temporal interpolation and long-term forecasting tasks using standard benchmarks, including Rotated MNIST, Health MNIST, PhysioNet, and SPRITES. The results demonstrate a significant performance improvement over competing models. However, some relevant models are missing from the benchmark comparisons. Including these baselines would provide a more comprehensive and rigorous evaluation of DGBFGP’s expressiveness. \n\n2) Computational efficiency: The paper evaluates computational efficiency by comparing worst-case complexity (Big-O notation), a standard metric for algorithmic runtime. The authors show that DGBFGP achieves a computational complexity of: $O(N\\sum_{p}B^{r})$ which improves upon previous approaches such as Ramchandran et al. (2021) and Jazbec et al. (2021), which have a complexity of $O(NM^2+M^3)$ where M is the number of inducing points.\n\nAdditionally, empirical runtime comparisons are conducted across four datasets (Rotated MNIST, Health MNIST, PhysioNet, and SPRITES), demonstrating that DGBFGP runs significantly faster than LVAE. However, the evaluation is limited to LVAE, and it would be more comprehensive to include runtime comparisons with additional baselines for a more complete assessment of computational efficiency.", "methods_and_evaluation_criteria": "The proposed methods and evaluation criteria are generally appropriate, but there are areas where the evaluation could be more comprehensive.\n\nThe authors evaluate model expressiveness by benchmarking DGBFGP on standard datasets (Rotated MNIST, Health MNIST, PhysioNet, and SPRITES), covering both temporal interpolation and long-term forecasting tasks. While the results show a significant performance improvement, the absence of certain baseline models makes it difficult to fully assess the model’s expressiveness. Including additional comparisons would strengthen the claims.\n\nFor computational efficiency, the paper provides both theoretical analysis (Big-O complexity) and empirical runtime comparisons. The authors demonstrate that DGBFGP reduces computational complexity relative to prior GP-based VAEs and empirically show that it runs significantly faster than LVAE. However, the runtime comparison is only performed against LVAE, and a more thorough evaluation against multiple baselines would provide a clearer understanding of the model’s efficiency across different settings.", "theoretical_claims": "No original theoretical claims are made in the paper.", "experimental_designs_or_analyses": "The experimental design is generally appropriate. Authors evaluate both model expressiveness and computational efficiency across multiple datasets. The performance comparisons demonstrate significant improvements, but the absence of certain baseline models raises concerns about the completeness of the evaluation. Additionally, while the runtime analysis supports the claim of improved efficiency, it is only compared against LVAE, limiting broader conclusions. Including additional baselines for both performance and runtime would strengthen the validity of the results.", "relation_to_broader_scientific_literature": "The paper extends Variational Autoencoders (VAEs) with Gaussian Process (GP) priors for multivariate time-series modeling by introducing a Bayesian basis function approach to approximate mixed-domain additive Gaussian Processes. The authors provide a comprehensive review of existing GP-based VAEs for time-series applications (Sohn et al., 2015; Cao et al., 2018; Iakovlev et al., 2023; Casale et al., 2018; Fortuin et al., 2020; Jazbec et al., 2021; Ramchandran et al., 2021), highlighting the key challenges associated with each method. They then introduce their proposed model, DGBFGP, which addresses these limitations by improving scalability and efficiency while maintaining the expressiveness of GP-based VAEs.", "essential_references_not_discussed": "The paper provides a comprehensive overview of prior work on Gaussian Process (GP) priors in Variational Autoencoders (VAEs) and related applications. However, it does not cite the foundational work on Variational Inference (VI) by Hoffman et al. (2013), despite heavily relying on concepts from Stochastic Variational Inference (SVI) for scalable inference.\n\nHoffman, M. D., Blei, D. M., Wang, C., & Paisley, J. (2013). Stochastic Variational Inference. Journal of Machine Learning Research, 14(1), 1303–1347.", "other_strengths_and_weaknesses": "The paper is well-written with a clear goal and a strong focus on improving GP-based VAEs for time-series forecasting. It effectively highlights its advancements over existing models but could benefit from stronger motivation for its approach. \nThe paper has a highly specific application: time-series forecasting using VAEs with Gaussian Process priors. It would be helpful to contextualize its relevance. Specifically, how does it compare to non-Transformer-based alternatives, and what are its broader implications and use cases? Addressing these questions would strengthen the paper’s impact and clarify its significance in the broader time-series modeling landscape.", "other_comments_or_suggestions": "I have no additional comments or suggestions beyond those already discussed in the other sections.", "questions_for_authors": "Expressiveness and the Role of $R$\n1) The expressiveness of the model likely depends on $R$, the number of latent additive dimensions. How does the model's performance vary as a function of $R$?\n2) I noticed that the latent dimensionality in experiments ranges from 16 to 64. Is there a specific reasoning behind this choice? Was there a trade-off considered between expressiveness, computational efficiency, and overfitting?\n3) I may be misunderstanding aspects of the paper, but given that DGBFGP approximates the additive Gaussian Process kernels used in LVAE, what specifically allows it to achieve significantly better performance?\n4) Is the improvement primarily due to a better approximation method, improved inference procedure, or other architectural differences?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The authors provided a strong rebuttal to the reviewers' critiques. In particular, the authors gave a convincing answer to the novelty of their method compared to other GP random features approximation methods. The experimental results also further strengthened their paper. Therefore, I vote to accept.
Decision: Accept (poster)