Dear Reviewer 1bD1,

Thank you for your time and effort spent on reviewing our paper. In the following, we provide responses to your specific concerns and questions.

---

The relationship to related works regarding the handling of missing values seems not to be clearly presented. For example, how is the proposed method related to

• [Smola et al., Kernel methods for missing variables, 2005]
• [Jafrasteh et al., Gaussian processes for missing value imputation, 2023]
• [Imani et al, Nested Gaussian process modeling and imputation of high-dimensional incomplete data under uncertainty, 2019]

Thank you for providing these additional references! They are very interesting and insightful. However, we want to point out that the type of missing values considered in these references differs from our setting.

(1) Missing Input Features

Suppose we are trying to model a latent function with two input dimensions and one output dimension:

(x_1, x_2) --> y

The references which you provided deal with the scenario where, for example, x_1 is missing in some training data examples and x_2 is missing in others. For n examples of d-dimensional training inputs, this corresponds to missing features in the n x d input data matrix.

(2) Missing Output Observations in a Cartesian Product Space

Our work considers the case of input data which can be expressed as a Cartesian product but is missing output values for some of the input data points.

Suppose we are modeling a latent function

(x_1, x_2) --> y

where x_1 can be a, b, or c, and x_2 can be 1, 2, or 3, such that the whole input space is

In our setting, missing values refer to missing output observations in the structured input space, for example, missing the y value corresponding to (b, 1) or (c, 3). However, each input observation itself is complete, i.e., we do not consider missing values such as (b, ?) or (?, 3).

In particular, (1) requires imputation of missing values in individual training examples before such data can be used for GPs (or modeling in general). In contrast, (2) does not obstruct the application of GPs in general, it just prevents the use of scalable methods based on Kronecker products. We focus on improving scalability by proposing a method which facilitates the use of Kronecker products in the presence of (2). Imputing missing values is not the main goal of our contribution. Therefore, our related work discussion and empirical evaluation focuses on other scalable GP methods, but we are happy to include a paragraph in the related work section of our camera-ready version which discusses techniques for (1) and clarifies the differences to (2).

Smola et al. (2005) also mention "incomplete labels" in a setting with multiple outputs. In the context of GPs, it is possible to treat multiple outputs as an additional (discrete) input dimension:

(x, idx) --> y[idx]

In this case, if we consider the Cartesian product {x} x {idx}, "incomplete labels" would correspond to (2). However, Smola et al. (2005) are proposing a method to integrate these missing labels out, whereas, in our setting, "incomplete labels" would just mean that we have less training data.

---

I cannot follow the sentence in pp. 2, " While the above methods scale gracefully to large datasets, they are all are fundamentally limited by performing inference of approximate GP models.". Why are the sparse and variational methods fundamentally limited?

Sparse and variational methods are fundamentally limited because they consider a small set of inducing points to represent a large number of training data points. Although these methods are typically designed to minimize the KL divergence to the true posterior ("get as close as possible to the model which uses all training data points"), in general, they cannot match the true posterior if the number of inducing points is smaller than the number of actual training data points, where the latter is required to achieve computational benefits. Therefore, "fundamentally limited" refers to the fact that sparse and variational GPs consider inducing point approximations of the true model, where true refers to using all training data points. In contrast, iterative methods consider inference in the true model using numerical techniques to solve the associated linear systems, which are capable of returning the true posterior up to numerical precision.

---

We believe that we addressed your concerns. If you have follow-up questions, please do not hesitate to ask them. Otherwise, we kindly ask you to consider increasing your score.

Dear Reviewer UJyJ,

Thank you for your time and effort spent on reviewing our paper. In the following, we provide answers to your specific questions.

---

Can the authors elaborate on why results in Table 1 showcase improved predictive RMSE performance of SVGP, yet not best negative log-likelihood?

• I understand that LKGP might produce the most sensible uncertainty estimates, yet would appreciate if the authors would elaborate on the reasons/hypothesis for the RMSE gap between LKGP and SVGP

Great question! Our hypothesis is that LKGP might be overfitting while SVGP might be underfitting. The setting is prone to overfitting because values are not missing uniformly at random. Instead, missing values are intentionally accumulated at the end of individual learning curves to simulate early stopping, such that the train (observed) and test (missing) data do not have the same distribution. Thus, fitting the training data distribution will be prone to overfitting. In particular, LKGP uses the whole training dataset and achieves a tight fit due to the exact GP model, potentially leading to overfitting. SVGP is known to be prone to underfitting, which tends to result in higher observation noise and larger length scales. The latter might actually protect against overfitting in this setting.

In Figure 4, LKGP achieves a tight fit with high certainty on the observed (solid black) parts of the learning curves, whereas SVGP does not always fit the observed data well while also estimating a high amount of observation noise and virtually no epistemic uncertainty. The high amount of observation noise can lead to bad NLL despite good mean predictions, while the lack of adequate epistemic uncertainty will penalize the NLL when the mean prediction is inaccurate.

Furthermore, we can compare the performance between train and test sets:

LKGP achieves a higher average rank on train RMSE, which supports the fact that LKGP may be overfitting to the training data distribution, and SVGP ranks worse in terms of train RMSE, suggesting that SVGP might be underfitting. The fact that LKGP still achieves the best test NLL is likely due to superior uncertainty quantification of the exact GP model compared to the variational approximation, which is a commonly observed phenomenon.

We will add a summary of this discussion to the camera-ready version of our manuscript.

---

A key novelty, and critical for the method's efficient use, is the efficient implementation of the latent kernel projections, without explicitly instantiating or multiplying by P, based on zero-padding, slice indexing and lazy evaluations.

• Can the authors elaborate further on these tricks/techniques, to help unfamiliar readers and for the completeness of the manuscript?

Thank you for this suggestion! Below, we provide explanations and basic pseudocode examples, which we will also include in the camera-ready version of the manuscript. In practice, we leverage GPyTorch and its underlying linear operator library, which provide more sophisticated implementations.

Latent Kernel Projections

In our setting, a projection matrix $P$ is a rectangular matrix which consists of an identity matrix where the rows corresponding to missing values are removed.

For example,

$P = \begin{bmatrix} 1 & 0 & 0 & 0 \\\ 0 & 0 & 1 & 0\end{bmatrix}$,

is a submatrix of a 4 x 4 identity matrix without its 2nd and 4th row.

Multiplying $P$ with a vector removes the 2nd and 4th element.

$P \begin{bmatrix} 1 \\\ 2 \\\ 3 \\\ 4 \end{bmatrix} = \begin{bmatrix} 1 \\\ 3 \end{bmatrix}$

Multipliying $P^T$ with a vector inserts zeros as 2nd and 4th element.

$P^T \begin{bmatrix} 1 \\\ 3 \end{bmatrix} = \begin{bmatrix} 1 \\\ 0 \\\ 3 \\\ 0\end{bmatrix}$

To implement $P$ efficiently, we can store the indices of rows which were not removed from the identity matrix. For the example above, idx = [0, 2].

We can avoid explicit matrix multiplications via

x = [1, 2, 3, 4]

Px = x[idx]

and

x = [1, 3]

PTx = zeros(4)

PTx[idx] = x

Lazy Kernel Evaluation

A $n \times n$ kernel matrix $K$ requires $\mathcal{O}(n^2)$ space, which quickly becomes intractable. However, to multiply $K$ with a vector $x$, we do not need to instantiate the full matrix. Instead, given a kernel function which can individually compute each value in $K$, we can iterate over rows of $K$ to perform the calculation using $\mathcal{O}(n)$ space.

for i in range(n): result[i] = k_i @ x

Here, k_i is the i-th row of the kernel matrix, and we rematerialize its values whenever needed to save memory.

---

We believe that we addressed your questions. If you have follow-up questions, please do not hesitate to ask them. Otherwise, we kindly ask you to consider increasing your score.

Dear Reviewer Zcew,

Thank you for your time and effort spent on reviewing our paper. In the following, we provide answers to your specific questions.

---

(i) It is unclear if the method is still applicable in high-dimensional spaces since we know that for Kronecker-product formulations the number of gridded data (or inducing points) increases exponentially with the number of input dimensions. Does the current method have better scaling with the input dimensionality d?

(ii) Related to the above it is unclear how the authors apply their method to the Sarcos dataset. The Sarcos dataset has input dimensionality d=21, so did you run using this dimensionality or much lower? Also these inputs as far as I know are not gridded. Please state exactly in all datasets which is the dimensionality of the "space" input vector.

Great question! We want to clarify that our method does not impose a grid over each individual input dimension. Instead, we assume that the full input data can be expressed as a Cartesian product of two (or more) factors, where each factor does not have to be structured and can contain multiple dimensions. Please allow us to explain this in the context of our experiments:

Inverse Dynamics Prediction

There are 7 joints with associated position, velocity, acceleration, and torque, and we wish to predict the 7 torques. In our framework, this can be viewed as 22 inputs --> 1 output,

(pos_1, vel_1, acc_1, ..., pos_7, vel_7, acc_7, joint_id) --> torque,

where joint_id ranges from 1 to 7 and specifies which joint's torque the output corresponds to. Here, we partition the input dimensions as

(pos_1, vel_1, acc_1, ..., pos_7, vel_7, acc_7) x (joint_id,),

where the first factor corresponds to an arbitrary configuration of the 7 joints, and joint_id is 1, 2, 3, ..., 7. Therefore, despite joint configurations being 21-dimensional, the full input data is a Cartesian product with two factors of size n_configs and 7 respectively, with n_configs = 5000 in our experiments. Crucially, the joint configurations (pos1, vel1, acc_1, ..., pos_7, vel_7, acc_7) do not need to have any particular structure.

Learning Curve Prediction

We view the problem as

(x_1, ..., x_k, t) --> learning curve value,

and the Cartesian product is

(x_1, ..., x_k) x (t,),

where (x_1, ..., x_k) is a particular hyperparameter configuration and t is the time step of the corresponding training run. Again, (x_1, ..., x_k) do not need to have any structure and the Cartesian product has two factors, with n_configs = 2000, k = 7, and t = 1, ..., 52 in our experiments, such that the full input data is of size 2000 x 52 = 104k.

Climate Data

We model the problem as

(lat, long, t) --> temperature / precipitation,

and partition the input dimensions as

(lat, long) x (t,),

which again leads to a two-factor Cartesian product, consisting of n_locs = 5000 pairs of non-gridded (lat, long) and t = 1, ..., 1000 days of history per location. The dense input data contains 5000 x 1000 = 5M observed temperature / precipitation values. Notably, we could also consider latitudes and longitudes on a Cartesian product space, resulting in an overall three-factor Cartesian product. In general, our method allows an arbitrary number of factors, where each factor considers a distinct set of input dimensions.

---

(iii) Also it will be useful to add a few more details (e.g., in the Appendix) about how do you maximize the log marginal likelihood and particularly how do you deal with the derivatives of the log determinant term.

Thank you for this suggestion! We have included a section in the appendix to explain marginal likelihood optimization and dealing with corresponding derivatives using iterative methods.

In summary, we use conjugate gradients in combination with fast Kronecker MVMs to solve linear systems, and we estimate the trace term in the derivative of the log det using Hutchinson's trace estimator.

To be more explicit, the derivative of the log marginal likelihood is

$\frac{1}{2} (K^{-1}y)^T \frac{\partial K}{\partial \theta} (K^{-1} y) - \frac{1}{2} tr(K^{-1} \frac{\partial K}{\partial \theta})$,

where $K$ is the kernel matrix, $y$ are the training labels, and $\theta$ are the kernel hyperparameters.

Using Hutchinson's trace estimator, the trace term is approximated as

$\frac{1}{s}\sum_{j=1}^s z_j^T K^{-1} \frac{\partial K}{\partial \theta} z_j$,

where $z_j$ are random probe vectors with $\mathbb{E} [z_j z_j^T] = I$.

To circumvent the explicit inverse of $K$, we instead calculate $K^{-1} y$ and $z_j^T K^{-1}$ by solving $s + 1$ linear systems,

$K [v_y, v_1, ..., v_s ] = [y, z_1, ..., z_s]$,

using conjugate gradients, which only relies on MVMs with $K$ and leverages the Kronecker structure in $K$.

---

We believe that we addressed your questions and concerns. If you have follow-up questions, please do not hesitate to ask them. Otherwise, we kindly ask you to consider increasing your score.

Dear Reviewer F5JD,

Thank you for your time and effort spent on reviewing our paper. In the following, we address your suggestions, concerns, and questions.

---

I couldn’t find in the paper whether multiple runs were performed with different selections of missing data for the inverse dynamics experiment.

• Report the mean and standard deviation over multiple runs of the experiments with different random seeds to ensure more robust results.

Thank you for this suggestion! All three of our experiments on real-world data were performed using different data splits, and we will include standard errors over multiple runs in the camera-ready version. In the following, we provide details about our random seed and data split procedures.

Inverse Dynamics Prediction

We considered 10 different random splits of the data, as stated in the caption of Figure 3. Each random split considers a different subset of 5000 joint configurations (out of a total of 44484 provided by the full Sarcos dataset) as well as a different pattern of uniformly random missing values.

The relative standard errors (relative to corresponding metric value, averaged across missing ratios) are small, demonstrating stability across replications with different random seeds.

Learning Curve Prediction

We considered 10 different random seeds, as mentioned in the captions of Table 1 and 3. Each random seed uses the same data but a different pattern of missingness, akin to random train / test splits. In particular, out of n learning curves, we first select missing_ratio x n curves uniformly at random. For each of the selected curves, we choose the early stopping point uniformly at random from [0, t_max - 1], ensuring at least 1 missing value. Thus, for each random seed, the split between complete / incomplete curves is randomized, and the missingness patterns within incomplete curves are also randomized.

Akin to the above, the relative standard errors (averaged across datasets) are small.

Climate Data

Due to the large size of the datasets (5M), we only considered 5 different random seeds. The 5M data points themselves are selected once from the NGCD data by drawing 5000 pairs of (latitude, longitude) uniformly at random and taking their corresponding 1000 days of history starting on Jan 1, 2020. For each of the 5 random seeds, the missingness pattern was adapted uniformly at random.

Again, the relative standard errors (averaged across datasets and missing ratios) are small.

---

Additionally, the initialization of hyperparameters and variational parameters is crucial for the predictive performance of both Gaussian process (GP) and sparse GP methods. The fact that GPyTorch’s default initialization was used as-is could negatively affect the results.

Thank you for your comment. Regarding kernel hyperparameters and observation noise, we believe that our experimental design is fair, despite using GPyTorch's default initialization, as, within each task, all GP models use the same kernel paired with the same hyperparameter initialization. While there might be a better initialization, it would arguably benefit all GP models. If certain models require initialization tuning to perform well, this could be considered a downside of that particular model, and similar efforts could be directed towards tuning other models to improve their performance.

Regarding variational parameters / inducing point locations, we initialized them at a uniformly random subset of the training data. This might be suboptimal, and initialization based on, e.g., k-means clustering could potentially lead to better performance. However, given that we consider large datasets of up to 5M data points, performing clustering or similar algorithms as part of the initialization can become expensive. For example, building the kNN table for VNNGP required several hours for each of the climate datasets, despite using the highly optimized faiss-gpu implementation (the time was not included in the reported benchmarks).

---

We believe that we addressed your concerns. If you have follow-up questions, please do not hesitate to ask them. Otherwise, we kindly ask you to consider increasing your score.