We thank the reviewer for their comments.

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.

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.

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.

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.

Compare to non-transformer-based alternatives. High-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.

Moreover, 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.

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:

$$f_{\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})$$

Ablation results for the Health MNIST by removing some components from the original implementation:

Removing 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.

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.

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.

For 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.

[1] Amortized Variational Inference: When and Why? https://arxiv.org/abs/2307.11018

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.

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.

To 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$^*$).
Our proposed model outperforms both the baseline and an amortized variant, as shown below:

This 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.

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.

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.

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.

Please 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.

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:

Simulated 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.

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).

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.

[1] https://proceedings.mlr.press/v80/cremer18a.html

[2] https://proceedings.mlr.press/v84/krishnan18a.html

[3] https://proceedings.mlr.press/v80/marino18a

[4] https://arxiv.org/abs/1810.11738

[5] https://arxiv.org/abs/1907.04155

[6] https://arxiv.org/abs/2010.13472

[7] https://arxiv.org/abs/2006.09763

[8] https://arxiv.org/abs/2302.04534

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.

Using 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.

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).

For 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$.

We will clarify these details in the revised manuscript.

[1] https://link.springer.com/article/10.1007/s11222-019-09886-w

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.

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).

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.

The 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.

For 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.

Finally, 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.