Scalable Gaussian Processes with Latent Kronecker Structure
Authors: Jihao Andreas Lin
Sebastian Ament
Maximilian Balandat
David Eriksson
José Miguel Hernández-Lobato
Eytan Bakshy
Venue: ICML 2025 poster (2025)
Keywords: Gaussian process
Kronecker
Bayesian
scalable
product kernel
gridded data
missing values
Primary area: probabilistic_methods->gaussian_processes
TL;DR: How to leverage efficient product kernels on gridded data with missing values for scalable GP regression

Abstract: Applying Gaussian processes (GPs) to very large datasets remains a challenge due to limited computational scalability. Matrix structures, such as the Kronecker product, can accelerate operations significantly, but their application commonly entails approximations or unrealistic assumptions. In particular, the most common path to creating a Kronecker-structured kernel matrix is by evaluating a product kernel on gridded inputs that can be expressed as a Cartesian product. However, this structure is lost if any observation is missing, breaking the Cartesian product structure, which frequently occurs in real-world data such as time series. To address this limitation, we propose leveraging latent Kronecker structure, by expressing the kernel matrix of observed values as the projection of a latent Kronecker product. In combination with iterative linear system solvers and pathwise conditioning, our method facilitates inference of exact GPs while requiring substantially fewer computational resources than standard iterative methods. We demonstrate that our method outperforms state-of-the-art sparse and variational GPs on real-world datasets with up to five million examples, including robotics, automated machine learning, and climate applications.

Scalable Gaussian Processes with Latent Kronecker Structure
Jihao Andreas Lin 123 Sebastian Ament 1 Maximilian Balandat David Eriksson'
José Miguel Hernandez-Lobato- Eytan Bakshy'
Abstract
GPs are notoriously challenging to scale to large numbers of training examples n. This primarily stems from having to solve an to compute the marginal likelihood and the posterior, which has (nn²) time complexity using direct methods.
Applying Gaussian processes (GPs) to very large datasets remains a challenge due to limited com- putational scalability. Matrix structures, such as the Kronecker product, can accelerate operations significantly, but their application commonly en- tails approximations or unrealistic assumptions. In particular, the most common path to creating Kronecker-structuree kernel matrix is by evalu- ating a product kernel on gridded inputs that can be expressed as a Cartesian product. However, this structure is lost if any observation is missing, breaking the Cartesian product structure, which frequently occurs in real-world data such as time series. To address this limitation, we propose leveraging latent Kronecker structure, by express- ing the kernel matrix of observed values as the projection of a latent Kronecker product. In com- bination with iterative linear system solvers and pathwise conditioning, our method facilitates in- ference of exact GPs while requiring substantially fewer computational resources than standard iter- ative methods. We demonstrate that our method outperforms state-of-thie -fth--fft sparse and variational GPs on real-world datasets with up to five million examples, including robotics, automated machine learning, and climate applications.
linear system involving the kernel matrix
To address these scalability challenges, a plethora of ap- proaches have been proposed, which usually fall into two main cafegories. Sparse GP approaches reduce the size of the underlying linear system by introducing set of inducing points to approximate the full GP: they include conventional (Quinonero- Candela & Rasmussen, 2005) and variational (Titsias, 2009) formulations. Iterative approaches (Gardner et al., 2018a Wang et al., 2019) employ solvers such as the linear conjugate gradient method, often leveraging hardware parallelism and structured kernel matrices.
In particular, Kronecker product structure permits efficient matrix- multiplication (MVM). This structure arises when product kernel is evaluated on data points from a Cartesian product space. A common (though not the only) example of this is spatiotemporal data, as found in climate science for instance, where quantity of interest, such as temperature or humidity, varies across space and is mea- sured at regular time intervals. In particular, assuming ob- servations collected at p locations and q times result in pt data points in total, the kernel matrix generated by evaluat- ing a product kernel on this data can then be expressed as the Kronecker product of two smaller p x p and x q matrices. This allows MVM in Op22q +1²")+pp²) instead of (p²2q²) time and O + q²) instead of O(p²g²) space, such that iter- ative methods require substantially less time and memory. However, in real-world data, this structure frequently breaks due to missing observations, typically requiring reversion to generic sparse or iterative methods that either do not take advantage of special structure or entail approximations to the matrix.
1. Introduction
Gaussian processes (GPs) are probabilistic models prized for their flexibility, data efficiency, and well-calibrated un- certainty estimates. They play a key role in many applica- tions such as Bayesian optimization (Frazier, 2018: Garnett, 2023), reinforcement learning (Deisenroth & Rasmussen, 201 and active learning (Riis et al., 2022). However, exact
In this paper, we propose latent Kronecker structure, which enables efficient inference despite missing values by repre- senting the joint covariance matrix of observed values as lavily-evaluated product between projection matrices and a latent Kronecker product. Combined with iterative methods, which otherwise would not be able to leverage Kronecker structure in the presence of missing values, this reduces the asymptotic time complexity from 0(p2g²) 0pp2g+rq9))
'Meta UUniversity of Cambridge MMax Planck Institute for Intelligent Systems, Tübingen. Correspondence to: Jihao Andreas Lin <jandylin@ meta.com?.
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 20255 by the author(s).

---

Scalable Gaussian Processes with Latent Kronecker Structure
Latent Kronecker
Projected Kronecker
Joint Covariance
Kss
KTT
P
P T = P
P
Figure 1. Illustrative example considering input data points (($1, ti), (S1, (2), (S2, tr), ($2, t2); (S2, ta)], consisting of two out of three time steps spatial location s and three out of three time steps at spatial location 8g. Due the missing observation (81, ts), the inputs cannot be expressed as Cartesian product and the joint covariance matrix cannot be directly expressed as a Kronecker product. Latent Kronecker structure assumes that missing observations exist, such that latent Kronecker product can be leveraged for computational benefits. The original joint covariance matrix is exactly recovered by removing the missing rows and columns via projections.
and space complexity from O(p²q") to O(p² + )²). Impor- tantly, unlike sparse approaches, our method does not intro- duce approximations of the GP prior. This avoids common problems of sparse GPs, such as underestimating predic- tive variances (Jankowiak et al., 2020), and overestimating noise (Titsias, 2009; Bauer et al., 2016) due to limited model expressivity with a constant number of inducing points.
Pathwise Conditioning Instead of computing the mean and covariance of the posterior distribution (1), Wilson et al. (2020; 2021) directly express sample from the GP posterior as a random function
(-ly) -
- f(-) + K (xxx(XXx + ooI)- (fx + E)),
where f ~~ GP(,, k) is a sample from the GP prior. fx is its evaluation at the training data X, and e ~ (O, I²I) is random vector. For product kernels on Cartesian product spaces, pathwise conditioning improves the asymptotic com- plexity of drawing samples from the GP posterior (Maddox et al., 2021). We leverage this property in Section 3.
In contrast, our method, Latent Kronecker GP (LKGP), faciliiates highly scalable inference of exact GP models with product kernels. We empirically demonstrate the superior computational scalability, speed, and modeling performance of our approach on various real-world applications including inverse dynamics prediction for robotics, learning curve prediction, and climate modeling, using datasets with up to five million data points.
Sparse and Variational Methods Conventional GP im- plementations scale poorly with the number of training ex- amples n because posterior inference (1) requires linear solve against the kernel matrix, which scales as O(2)) using direct methods. To address this, large number of sparse GP (Snelson & Ghahramani, 2005) approaches that rely on approximating the kernel matrix at a set of inducing points, have been proposed. These methods typically reduce the training cost from O(n²) ) to Onnm + +m)) per gradient step where m is the number of inducing points and m < n.
2. Background and Related Work
Let f X + be stochastic process such that for any finite subset { C CX, the collection of random vari- ables {f(x)) follows multivariate normal distribution. This property defines
Gaussian process (GP), which is uniquely identified by its mean function va()) B[J(-)] and its kernel function (---)) = (())))(()))))) We suppress the dependence on kernel hyperparameters for brevity.
Sparse Gaussian Process Regression (SGPR) an induc- ing point method which introduces variational inference for sparse GPs (Titsias, 2009). The locations of the inducing points are optimized along with the kernel hyperparameters using the evidence lower bound (ELBO), Both the ELBO op- timization and posterior inference require O(mm²), which is significantly faster than exact GP inference when m < n.
To perform GP regression, let the training data consist of inputs }= = X C X' and outputs (y: P= = y E T.... We consider the Bayesian model y; = fxx.) + e,, where f GP(H,, k) and each e, ~N(00 o?). Without loss of generality, we assume μ = 0. The posterior of this model is given by fly ~ GP( k with
Stochastic Variational Gaussian Processes (SVGP) improve the computational complexity of SGPR by decomposing the ELBO into sum of losses over the training labels (Hensman et al., 2013). This makes it possible to leverage stochastic gradient- based optimization for training which reduces the computational complexity to Om3²) ). Consequently, SVGP can be scaled to very large datasets as the asymptotic train- ing complexity no longer depends on n.
Msy ()) (K)K))KK)))) (( (xxx o - K. x (xx + -
where Kxx denotes the covariance matrix [k(x,, x,) [)) and K( Kx and ((,,") denote pairwise evaluations of the same kind. For a comprehensive introduction to GPs, we refer to Rasmussen & Williams (2006).
2

---

Scalable Gaussian Processes with Latent Kronecker Structure
Peak Memory Usage GB
Kernel Evaluation Time in ms
Matrix Multiplication Time in ms
60
150
50
120
12
40
30
80
6
20
10
30
3
0
10k 20k
30k
40%
50k
100M
1B
10k 20k 30k 40k 50k
100M
1B
10k 20k 30k
40k 50k
1M 4M 7M 10M
Number of Training Examples
Number of Training Examples
Number of Training Examples
Without latent Kronecker structure
With latent Kronecker structure
Figure 2. Illustration of computational resources used during kerne evaluation and matrix multiplication on cn-dimensional synthetic datasets of different sizes. Without latent Kronecker structure, memory usage escalates quickly and kernel evaluation time dominates matrix multiplication time asymptotically. With latent Kronecker structure, computations can be scaled to several orders of magnitude larger datasets under similar computational resource usage, and matrix multiplication time dominates kernel evaluation time asymptotically. Reported results are the mean over one hundred repetitions using squared exponential kernel.
Wu al. (2022) proposed the Variational Nearest Neighbor Gaussian Process (VNNGP), which introduces prior that only retains correlations for K nearest neighbors. This results in sparse kernel matrix where the time complexity of estimating the ELBO is O(m,, + m²))3 for a mini- batch of size ny and a mini-batch of m, inducing points.
of b enable the estimation of marginal likelihood gradients for kernel hyperparameter optimization (Lin et al., 2024b). In terms of actual optimization algorithms used as iterative linear system solvers for GPs, Gardner et al. (201 8a) pop- ularized conjugate gradients, Wu et al. (2024) introduced alternating projections, and Lin et al. (2023; 2024a) applied stochastic gradient descent.
Wenger et al. (2024) introduced Computation-Awan Gaus- sian Processes (CaGP), class of GP approximation meth- ods which leverage low-rank approximations while dealing with the overconfidence problems of SVGP, guaranteeing that the posterior variance of CaGP is always greater than the posterior variance of an exact GP. The additional vari- ance is interpreted as computational uncertainty and helps predictions while scaling
In contrast sparse and variational approaches, iterative methods perform inference of the exact GP model up to certain numerical tolerance. We leverage iterative linear system solvers for efficient inference in Section 3.
Structured Kernels In the special case of product ker- nels with multidimensional data on a Cartesian grid, the kernel matrix decomposes into Kronecker product of one- dimensional kernels: K K₁ @ @ Ka- This allows expressing matrix decompositions of K (including Cholesky and eigendecompositions: as individual decompositions of K K,.... Ka. enabling efficient exact GP inference (Bonilla et al., 2007; Stegle et al., 2011; Golub Van Loan, 2013), including for gradient observations (Ament & Gomes, 2022). To complete partially observed grid of observations, Saatçi (2012) and Wilson et al. (2014) propose adding imaginary observations with a large noise variance, giving rise to an approximation which only converges as the artificial noise variance goes to infinity and leads to ill-conditioning.
the method achieve to large datasets.
While the above methods scale to large datasets, they are limited by performing inference of approximate GP models.
Iterative Linear System Solvers Alternatively, costly, or even intractable, direct solves of large linear systems can be avoided by leveraging iterative linear system solvers. To this end, the solution A -bb to the linear system Av b with positive-definite coefficient matrix A is in- stead equivalently expressed the unique optimum of convex quadratic objective,
v= A- +
arg min
1 11
Au T u u'b,
Structured Kernel Approximations In the case of non- gridded inputs,
which facilitates the of iterative optimization algorithms to solve large linear systems in scalable way using MVMs. In the context of GPs, setting A to Kxx + o²I and b to y produces the posterior mean (-) as K (xxv. Setting b to y - (fx + E) generates samples from the GP posterior via pathwise conditioning as f(.) K Other choices
leverage the Structured Kernel In- terpolation (SKI) kernel to approximate the original kernel: Kxx s WKuu T where W is matrix of interpola- tion weights and KUU exhibits structure such as Toeplitz or Kronecker (Wilson & Nickisch, 2015). This algebraic structure can be used for efficient kernel learning using
3

---

Scalable Gaussian Processes with Latent Kronecker Structure
Time
Memory GB
Test RMSE
Test NLL
75
0.15
-1.0
60
60
0.13
-1.1
0.11
-1.2
30
30
15
0.09
-13
0
0
0.07
-14
0.1
0.3
0.5 0.7
0.9
0.1
0.3
0.5
0.7
0.9
0.1
0.3
0.5 0.7
0.9
0.1
0.3 0.5 0.7 0.9
Missing Data Ratio
Missing Data Ratio
Missing Data Ratio
Missing Data Ratio
Iterative GP
Latent Kronecker GP
Asymptotic break-even point
Figure 3. Predicting the inverse dynamics of an anthropomorphic robot amm with seven degrees of freedom. Compared to standard iterative methods, leveraging latent Kronecker structure results in significantly lower runtime and memory requirements while maintaining the same predictive performance. The asymptotic break-even point, at which both methods asymptotically require the same amount of time memory, closely matches the empirical break-even point. Reported results are the mean over 10 different random splits the data.
conjugate gradients for linear solves and Lanczos to approx- imate the log determinant (Dong et al., 2017). SKI is limited low-dimensional inputs as the number inducing points m increases exponentially with the input dimensionality.
space approaches to GP regression are 1) the requirement that the temporal kernel be stationary, and 2) that they scale quadratically or cubically with the spatial dimensionality.
The method proposed herein scales linearly in the spatial dimensionality, does not require the temporal kernel to be stationary, and is not limited to spatiotemporal data. If the method is applied to spatiotemporal data collected in uniform temporal intervals, and the temporal kernel is sta- tionary, the method can be accelerated to be quasi-linear in the number of time steps by leveraging the Toeplitz structure of the temporal kernel matrix (Wilson & Nickisch, 2015).
To address the scalability issues of SKI with respect to the number of input dimensions, Gardner et al. (2018b) pro- posed Structured Kernel Interpolation for Products (SKIP), which approximates the original kernel by Hadamard prod- uet of mnc-dimensional SKI kernels. Every SKI kernel is further approximated by low-rank Lanezos decomposition, which gives rise an algebraic structure that allows for effi- cient MVMs: Kxx R (Q1 T Q1 " @ (Qa Ta QaQa" )) While SKIP addresses the scalability issues of SKI, it is usually not competitive with other methods due to the errors introduced by the approximations (Kapoor et al., 2021).
3. Latent Kronecker Structure
In this paper, we consider a Gaussian process f X' + R defined on a Cartesian product space X' S x T. For illustrative purposes, S could be associated with spatial dimensions and T could refer to time dimension or a task index, commonly used in multi-task GPs (Bonilla et al., 2007), but emphasize that neither of them are limited to one-dimensional spaces. A natural way to model this problem is to define a kernel hx on the product space X. which results in joint covariance
Hierarchical matrix approximation algorithms constitute an- other notable class of methods. Ambikasaran al. (2016) introduced approach which is based on purely algebraic approximations and permits direct matrix inversion. Ryan et al. (2022) proposed an algorithm which uses automatic differentiation and automated symbolic computations which leverage the analytical structure of the underlying kernel. The main limitation of this class methods is also scalabil- ity to higher-dimensional spaces, for which the hierarchical approximations become inefficient.
Cov(J()),[(x)) kx x') kx((s, t), (s', t')),
where s, E S and t, t' E T. However, this results in scalability issues when performing GP regression. the training data consists of n outputs {y, 1 = y R" observed at p spatial locations {s, = C S and q time steps or tasks [tx) T C T. such that n pq. then the joint covariance matrix requires O(p2g²) space and computing its Cholesky factor takes Opp' time.
State-Space Approaches Hartikainen & Särkkä (2010); Sarkka et al. (2013) introduced state-space formulations for (spapio-)temporal GP regression, whose computation has linear complexity with respect to the number of tempo- ral observations for given point in space. See Särkkä & Svensson (2023) for review of the methods. Recently, La- plante et al. (2025) leveraged this formulation to accelerate robust spatiotemporal GP regression. Pförtner et al. (2025) extended the computation-awe aware framework of Wenger et al. (2024) to quantify errors caused by approximate filtering and smoothing methods. The main limitations of the state-
Ordinary Kronecker Structure A common way to re- duce the computational complexity in such settings is to introduce product kernels and Kronecker structure (Bonilla

---

Scalable Gaussian Processes with Latent Kronecker Structure
2
Valirlatinn Loss
6
2
0
6
2
0
-2
20
40
Training Epochs
Training Epochs
Training Epochs
Training Epochs
Observed
Missing
LKGP
SVGP
VNNGP
CaGP
Figure Learning prediction the FaSiion-MNIST MNIST data from the LCBench benchmark (Zimmer al., 2021). Partially observed learning curves are extrapolated into the future. The predictive mean and two standard deviations various GP models are visualized. All models reasonably accurate in the mean, but LKGP produces the most sensible uncertainty estimates, starting with low uncertainty in the observed part of the learning curve and gradually increasing the predicted uncertainty into the missing part. The third row illustrates outlier with significantly different behavior than other learning curves. The sparse methods struggle because the limited number of inducing points is unlikely to cover such case. I.GPP can adapt well because it performs inference the exact GP model.
et al., 2007; Stegle al., 2011; Zhe et al., 2019). Defining xx,xxx)) xx ((s, t), (s', ks(s, (t,
proposing method which facilitates the use of Kronecker products in the presence of missing values.
Dealing with Missing Values In the context of gridded data with missing values, such that n < a key insight is that the joint covariance matrix over observed values is submatrix of the joint covariance matrix over the whole Cartesian product grid. Although the former generally does not have Kronecker structure, the latter does. We leverage this latens Kronecker structure in the latter matrix by pressing the joint covariance matrix over observed values as
where ky only considers spatial locations s, and by solely acts on time steps or tasks tg, allowing the joint covariance matrix Kxx to factorize as
Kxx ks(S,S) & (T.T) Kss @ KTT, nxn
pxp qxq
which can be exploited by expressing decompositions of Kxx in terms of decompositions of Kss and KTT instead. This reduces the asymptotic time complexity to Opp + q²) and space complexity Opp 2²)
Kxx = P (Kss & KTT) PT p4Xn
XXX пxPς
pxp qxq
However, the joint covariance matrix Kxx only exhibits the Kronecker structure if observations y; are available for each time step or task t2 at every spatial location Sj, and this is often not the case. For example. suppose we are con- sidering temperatures y; on different days t& measured by various weather stations at locations Sj. a single weather station does not record the temperature on any given day, then the resulting observations are no longer fully gridded, and thus ordinary Kronecker methods cannot be applied. These missing values refer to missing output observations which correspond to inputs from a Cartesian product space X' = S x T. They do not refer to missing features in the input data, and imputing missing values is not the main goal of our contribution. We focus on improving scalability by
where the projection matrix P is constructed from the iden- tity matrix by removing rows which correspond to missing values. See Figure for an illustration. Crucially, this not an approximation and facilitates inference of the exact GP model. In practice, we implement projections efficiently without explicitly instantiating multiplying by P.
Efficient Inference via Iterative Methods As result of introducing projections, the eigenvalues and eigenvectors of Kxx cannot be expressed in terms of cigenvalues and eigenvectors of Kss and KTT anymore, which prevents the use of Kronecker structure for efficient matrix factorization. However, despite the projections, the Kronecker structure
5

---

Scalable Gaussian Processes with Latent Kronecker Structure
can still be leveraged for fast matrix multiplication. To this end, augment the standard Kronecker product equation, (A @ B) vec(C) vec vec(BCA)), with the aforementioned projections, yielding
in this case, the cost of kernel evaluations is amortized, and the total runtime is dominated by matrix multiplication. But, if kernel matrices do not fit into memory, their values must be rematerialized during each matrix multiplication, leading to many repeated kernel evaluations which may dominate the total runtime. Due to much lower memory requirements, latent Kronecker structure also raises the threshold after which kernel values must be rematerialized. Figure illus- trates memory usage, kernel evaluation times, and MVM times on ten-dimensional synthetic datasets of various sizes, assuming balanced factorization p = Vn.
P(A @ vec(C) = Pveo Bvec '(PT AT)
In practice, vec and vec are implemented as reshaping operations, PT vCc amounts to zero padding, and left- multiplying by P is equivalent slice indexing, all of which typically do not incur significant computational overheads. This facilitates efficient inference in the exact GP model via iterative methods, which only rely MVM to compute so- lutions of linear systems (Gardner et al., 201 8a; Wang et al., 2019). Leveraging the latent Kronecker structure reduces the asymptotic time complexity of MVM from (n²²) to Opp2q + p."), and the asymptotic space complexity from Onn²) to O(p² + q²). Furthermore, if the kernel matrix is not materialized and kernel values are instead calculated demand, also known as lazy kernel evaluation, the asymp- totic space complexity is reduced from (nz) to (p + q).
Efficiency of Latent Kronecker Structure Since latent Kronecker structure can be understood as padding n S pa observations to complete a grid of pg values, a key question is at which point the number of missing values pq - n is large enough such that this padding becomes inefficient and dense matrices without this structure are preferable.
= n/pq be the missing ratio, that is, the relative amount of padding required to complete a grid with pq values. Let be the asymptotic break-even point, that is, the particular missing ratio at which a ker- nel matrix without factorization has the same asymptotic performance as a kernel matrix with latent Kronecker struc- ture. The asymptotic break-even points for MVM time and memory usage are, respectively.
Proposition 3.1.
Posterior Samples via Pathwise Conditioning In the context of product kernels and gridded data, Maddox et al. (2021) proposed to draw posterior samples via pathwise con- ditioning (Wilson et al., 2020; 2021) to exploit Kronecker structure, reducing the asymptotic time complexity from O(p'a") to O(p" + -²)). In our product space notation, the pathwise conditioning equation can be written as
/1 Vhime V p
V
1
1
and m =
q
pa
q²
(fly) ()) = f(((=: ())) + (K)ss i @ ((T)) V.
See Appendix for the derivation. In Section 4, we validate these asymptotic break-even points empirically, and observe that they closely match the empirical break-even points obtained from experiments on real-world data (see Figure 3).
where v is the inverse matrix-vector product
v = (Kss @ KTT + o'1)- (fsXT + t));
and e ~ (O,,02)) To support latent Kronecker structure, we again introduce projections such that
Relationship to High-order GPs Zhe et al. (2019) pro- posed High-order GPs (HOGPs), which model tensor- structured outputs through a kernel defined over low- dimensional latent embeddings of each coordinate, allowing efficient computation via Kronecker algebra. HOGPs and LKGPS both leverage Kronecker products for scalability but differ fundamentally: HOGPs require a complete grid of observations and model dependencies via inferred latent variables, while IKGPss handle missing values via latent projections of a full structured kernel matrix (without infer- ring any latent variables).
v PT (P(KSs
resulting in exact samples from the exact GP posterior. In combination with iterative linear system solvers, latent Kro- necker structure can be exploited to compute using effi- cient MVM. Further, computations can be cached or amor- tized to accelerate Bayesian optimization or marginal likeli- hood optimization (Lin et al., 2024b).
Discussion of Computational Benefits The most appar- ent benefits of using latent Kronecker structure are the im- proved asymptotic time and space complexities of matrix multiplication, as discussed earlier. However, there are fur- ther, more subtle benefits. In particular, the number of kernel evaluations changes from O(n?) evaluations of xx to Opp?) evaluations of ks plus O(g2) evaluations of ky. This may seem irrelevant if kernel matrices fit into memory, because,
4. Experiments
We conduct three distinct experiments to empirically evalu- ate LKGPs: inverse dynamics prediction for robotics, learn- ing curve prediction in an AutoML setting, and prediction of missing values in spatiotemporal climate data. In the first experiment, we compare LKGP against iterative methods
6

---

Scalable Gaussian Processes with Latent Kronecker Structure
Table 1. [eaarning Curve Prediction] Predictive performances and total runtimes of learning curve prediction on every fifth I.C.ench dataset. On average, LKGP produces the best negative log- lileliiiioodd while also requiring the least time. Reported numbers are the mean over 10 random seeds. Results of Appendix B for the results on all datasets.
model are boldfaced and results the second best model are underlined. See Tables to
Model
APSFailure
MiniBooNE
blood
covertype
higgs
kr-vs-kp
segment
Average Rank
LKGP
1.705 + 0.047
0.047 = 0.001
0.278 = 0.017
0.061 = 0.001
0.039 = 0.001
0.008 = 0.000
0.006 + 0.000
1.171 = 0.560
Troin
RMSE
SVGP
3,893 1 0.149
0.201
0.004
0.463 1 0.009
0.240
0.003
0.200
0.003
0.101 1 0.001
0.086 1 0.001
2.771 1 0.680
VNNGP
4,920 + 0.194
0.113
0.001
0.284 + 0.029
0.158
0.003
0.130
0.002
0.067 + 0.001
0.072
0.001
2.314 + 0.708
CaGP
3,972 : 0.147
0.219
0.005
0.511 + 0.011
0.241
0.004
0.208
0.002
0.105 + 0.002
0.091
0.001
3.743 + 0.553
LKGP
2.935 1: 0.150
0.335
1: 0.006
0.747 1 0.016
0.351 = 0.008
0.394 = 0.027
0.261 1 0.005
0.152 1 0.003
2.600 1 0,901
Test RMSH
SVGP
2.716 + 0.145
0.316
0.005
0.540 = 0.013
0.294 = 0.005
0.285
0.005
0.232 + 0.003
0.145 + 0.002
1.657 = 1.068
VNNGP
3,053 +: 0.159
0.705
0.01
0,915
0.017
0.677 = 0.011
0.642 = 0.009
0,591 + 0.006
0.568 + 0.005
3.743 = 0.731
CaGP
2.810 +: 0.154
0.325
0.006
0.593 +: 0.015
0.281
0.004
0.296
0.004
0.239 +: 0.004
0.153 1: 0.002
2.000 + 0.000
LKGP
1.955 1 0.020
-1.573 1 0.022
0.160 1 0.055
-1.311 T 0.017
11.118
0.018
-2,933
0,005
-2.943 1 0,005
1.029 1 0.167
Train NLL
SVGP
3.098 + 0.100
-0.132
+ 0.020
0.652 + 0.021
0.025 = 0.013
-0.150 + 0.012
-0.800 + 0.014
-0.977 + 0.007
2.400 + 0.685
VNNGP
3,005 + 0.041
-0.079
# 0.010
0.639 # 0.045
0.104 = 0.009
-0.056 = 0.009
-0.489
= 0.008
-0.456
+ 0.006
3.286 = 0.881
CaGP
2.794 1: 0.045
-0.042
1: 0.022
0.798
1 0.021
0.034
0.017
-0.103
0.012
-0.718
0.013
-0,900
1 0.009
3.286 1 0.564
LKGP
2.349 + 0.047
0,051
0.026
1.047 0.029
-0.057
0.019
-0.328
0.023
-0.633
0,029
-1,214 + 0.020 0.020
1.257 = 0.602
Test
NLL
SVGP
2.327 + 0.065
0.319
0.781
0.020
0.199
0.016
0.160
0.018
0.065 + 0.032
-0.496 + 0.014
2.543 = 1.024
VNNGP
2.733 1: 0.031
1.046
0.017
1.365 + 0.019
1.011 = 0.016
0.891
= 0.013
0.717
+ 0.013
0,712 + 0.009
3.200 + 1.141
CaGP
2.451 + 0.033
0.321 + 0.021
0.916 = 0.023
0.142 = 0.016
0.171 = 0.013
0.088 + 0.031
-0.463 + 0.013
3.000 = 0.000
LKGP
0.371 +: 0.008
1,458
0.013
0.487
0.003
1.728
0,010
2.447
0.022
1.493
0.017
2.012 11 0.018
1.000 1 0.000
Time min
SVGP
6.473 + 0.032
6,475
0.032
6.479 + 0.031
6.474
0.032
6.492
0.033
6.500
+ 0.030
6.475 + 0.032
2.000 + 0.000
VNNGP
26.34 + 0.443
25.89 = 0.161
25.78 = 0.151
25.85
0.158
25.87
0.181
25,85 = 0.154
26.31 0.464
4.000 = 0.000
CaGP
7,067 1 0.028
7.036
0.028
7.033 1 0.023
7.015
0.020
7.046 0.024
7.024 1 0.019
7,023 1 0.020
3.000 1 0.000
30
15
15
30
100
75
50
25
16
Jan 2020
Apr
Jul
Oct
Jan
Apr
Jul
Oct
Jan 2022
Apr
Jul
Oct
10
30
Temperature in "C
2021
Precipitation in kg/m²
Figure 5. Illustration of daily temperature and precipitation data from the Nordic Gridded Climate Dataset (Tveito et al., 2000; 2005). The heatmaps (left and right) visualize snapshots of single day and subsampled spatial locations. Every spatial location is associated with its own time series (middle). Temperatures (top) exhibit a seasonal periodic trend. Precipitation (bottom) is noisy but locally correlated.
without latent Kronecker structure. In the second and third experiment, compare method to various sparse and variational methods.
periments were conducted on A100 GPUs with total com- pute time of around 2000 hours.
For all our experiments, we start with a gridded dataset and introduce missing values which are withheld during training and used as test data, which provides ground truth informa- tion for missing values. We compute predictive means and variances of LKGP using 64 posterior samples obtained via pathwise conditioning, akin to Lin et al. (2023), All methods are implemented using GPyTorch (Gardner et al., 201 8a), and the source code is available at https: //gtt///////bb//// com/ homn/jandylll inn in http://attetttttrreeeeeeeeeee-ee Krooockernnre-eee-eeee All ex-
Inverse Dynamies Prediction In this experiment, we pre- dict the inverse dynamics of a SARCOS anthropomorphic robot arm with seven degrees of freedom. The SARCOS dataset is publicly available.' In this problem, we wish to learn the mapping of seven joint positions, velocities, and accelerations to their corresponding seven joint torques. The main goal of this experiment is to compare LKGPss to standard iterative methods, because the former effectively
hhttp://gauss: http://gaussaarrrrcees.o.ooo /datt/
7

---

Scalable Gaussian Processes with Latent Kronecker Structure
implements the latter using more efficient matrix algebra. Additionally, we empirically validate the break-even points predicted by Proposition 3.1, at which both methods should require the same amount of time or memory.
particular time step and missing all remaining values. This simulates neural network hyperparameter optimization, where learning curve information becomes available as the neural network continues to train. In this setting, learning curve predictions can be leveraged to save computational resources by pruning runs which are predicted to perform poorly (Elsken et al., 2019). We compare our method to SVGP (Hensman al., 2013), VNNGP (Wu al., 2022), and CaGP (Wenger et al., 2024). See Section 2 for details.
For this experiment, we consider the positions, velocities, and accelerations as S and the seven torques as T. We choose Kg to be squared exponential kernel and ky to be a full-rank ICM kernel (Bonilla al., 2007), demonstrating that LKGP is compatible with discrete kernels. We select subsets of the training data split with p- 5000 joint posi- tions, velocities, and accelerations, and their corresponding q = 7 joint torques, and introduced 10%, 20%, 90% missing values uniformly at random, resulting in n S 35k. The value of p was chosen such that the kernel matrix with- out factorization fits into GPU memory. This results in a fairer comparison since it allows both methods to amortize kernel evaluation costs, such that the comparison is essen- tially between matrix multiplication costs (larger p would force regular iterative methods to recompute kernel values, leading to even higher overall compute times). Details on the model fitting setup are provided in Appendix C.
For this experiment, we consider hyperparameter config- urations as S and learning curve progression steps as T. such that p = 2000, = 52, and thus n I 104000 for each dataset. We choose both kss and by to be squared exponential kernels. Out of the p = 2000 curves per dataset, 10% were provided as fully observed during training and the remaining 90% were partially observed. The stopping point was chosen uniformly at random. See Figure 4 for an illustration and qualitative comparison of the GP methods which we considered for this experiment. Details on the model fitting setup are provided in Appendix C.
Table reports the predictive performance and total runtime for every fifth LCBench dataset (see Appendix B for re- sults on all datasets). On average, SVGP and CaGP achieve better test wwt--mman--uuu-----eeeee duces the best test negative log- ligliiklidood (NLL), providing quantitative evidence for the qualitative observation in Fig- ure which suggests that LKGP produces the most sensible uncertainty estimates.
Figure 3 illustrates the required computational resources and predictive performances for different missing data ra- tios. For small missing ratios, LKGP requires significantly less time and memory. As the missing ratio increases, the iterative methods eventually become slightly more efficient. Notably, the empirical break-even points, where time and memory usage is the same for both methods, matches the asymptotic break-even points predicted by Proposition 3.1. Moreover, the predictive performance of both methods is equivalent across all missing ratios in terms of test root- meean-squaee-ery and test negative log- lieeleilloood, validat- ing that LKGP does not introduce any model approximation.
(RMSE), but LKGP pro-
To explain why LKGP performs worse in terms of test RMSE, we suggest that LKGP might be overfitting. Since missing values are intentionally accumulated at the end of individual learning curves to simulate early stopping, the train (observed) and test (missing) data do not have the same distribution, which makes this setting particularly prone to overfitting. Our hypothesis is empirically supported by observing that LKGP consistently achieves the best per- formance on the training data. The fact that LKGP still achieves the best test NLL is likely due to superior uncer- tainty quantification of the exact GP model compared to variational approximations, which is a commonly observed phenomenon. Furthermore, LKGP requires the least amount of time by a large margin, suggesting that it could potentially be scaled to much larger learning curve datasets.
Learning Curve Prediction In this experiment, we pre- dict the continuation of partially observed learning curves, given some fully observed and other partially observed ex- amples. We use data from LCBench (Zimmer et al., 2021).2
LCBench contains 35 distinct learning curve datasets, each containing 2000 learning curves with 52 steps each, where every step refers to a neural network training epoch. Each learning curve within a particular dataset is obtained by training a neural network on the same dataset but using different hyperparameter configurations. These hyperparam- eters are batch size, learning rate, momentum, weight decay, number of layers, hidden units per layer, and dropout ratio.
Climate Data with Missing Values In this experiment, we predict the daily average temperature and precipitation at different spatial locations. We used the Nordic Gridded Climate Dataset (Tveito et al., 2000; 2005).² This dataset contains observations which are gridded in space and time,
The main goal of this experiment is to evaluate the per- formance under realistic non-uniform pattern of missing values. In particular, learning curves are observed until
httpps: / htt://ittti thub. com/ / amt.mml /LCBench available under the Apache License, Version 2.0
dataseta/insicaa http:/innntttggrideeeee wirridddddddebeeeee wwwpaaiiinnnonnnriiooooo available under the License to Use Copernicus Products
8

---

Scalable Gaussian Processes with Latent Kronecker Structure
Table 2. [Climate Data with Missing Values] Predictive performances and total runtimes of temperature and precipitation prediction across
various missing data ratios. LKGP consistently achieves the best predictive performance and also requires the least time across both datasets and missing ratios. Reported numbers are the mean over random seeds. Results the best model are boldfaced and results of the second best model are underlined.
Temperature Dataset (with Missing Ratio 10%
50%)
Precipitation Dataset (with Missing Ratio 10%
50%)
Model
10%
20%
30%
40%
50%
10%
20%
30%
40%
50%
LKGP
0,06 0.00
0.06
0.00
0.07
0.07 ::
0.07 4: 0.00
0.16 +: 0.00
0.17
:: 0.00
0.17 +: 0.00
0.18
0.00
0.18 0.00
Truin
RMSE
SVGP
0.21 0.00
0.21
+ 0.00
0.22
0.00
0.22 + 0.00
0.23 0.00
0.70 = 0.00
0.71
+ 0.00
0.71 0.00
0.71 = 0.00
0.72 +: 0.00
VNNGP
0.00
0.13 1:3:000 0.00
0.13
0.13 = 0.00
0.13 0.00
0.45 = 0.00
0.47
+: 0.00
0.49 + 0.00
0.52 020:000 0.00
0.58 0.00
CaGP
0.18 0.00
0.19 0.00
0.19
0.19 + 0.00
0.19 +: 0.00
0.60 = 0.00
0.61
+: 0.00
0.61 +: 0.00
0.62 - 0.00
0.62 +: 0.00
LKGP
0.08 0.00
0.08
0.00
0.08 0.08 :: 0.00
0.09
= 0.00
0.09 + 0.00
0.24 = 0.00
0.25 :: 0.00
0.26 - 0.00
0.27 = 0.00
0.28 0.00
Test RMSE
SVGP
0.21 0.00
0:21
0.00
0.00
0.22 0.00
0.23
0.00
0.70 + 0.00
0.71 0.00
0.71
0.72 0.00
0
VNNGP
0.13 0.00
0.13
0.13 0:00
0.13 0.00
0.13
0,46 + 0.00
0.48
0.00
0,50 0.00
0.53 0,00
0.58 0.00
CaGP
0.18 0.00
0.19 + 0.00
0.19 +: 0.00
0.19 + 0.00
0.19 +: 0:00
0.61 + 0.00
0.61
+ 0.00
0.61 + 0.00
0.62 + 0.00
0.63 + 0.00
LKGP
-1.33 +: 0.02
-1.30 = 0.01
11266
- 0.00
-1.21
0.00
-1.18
= 0.00
-0,32 + 0.00
-0.29 + 0.00
-0.26 1: 0.00
-0.23 =: 0.00
+0.18 + 0.00
Train
NLL
SVGP
0..14 0.00
.1.12
+ 0.00
-0.10
+ 0.00
-0.07
+: 0.00
-0.05
+ 0.00
1,07 + 0.00
+ 0.00
+ 0.00
1.09
0.00
1.10 0.00
VNNGP
-0.58 +: 0:00
-0.58
= 0.00
-0.57
+ 0.00
-0.57
+ 0.00
-0.57
+ 0.00
0.64 = 0.00
0.67
+ 0.00
0.71 + 0.00
0.78 + 0.00
0.87 + 0.00
CaGP
-0.22 +: 0.00
-0.21
+: 0.00
-0.20
+ 0.00
-0.18
+ 0.00
-0.17 + 0.00
0,93 + 0.00
0.94
+ 0.00
0.94 + 0.00
0.95 + 0.00
0.96 + 0.00
LKGP
-1.14 +: 0.01
-1.11
= 0.01
-1.07
0.00
-1.03
+: 0.00
-0.99
0.00
-0.02 + 0.00
0.01 + 0.00
0.04
0.00
0.09 = 0.00
0.14 0.00
Test
NLL
SVGP
-0.14 +: 0.00
-0.12
+: 0.00
-0.09
+ 0.00
-0.07
+: 0.00
-0.05
+ 0.00
1.07 + 0.00
1.08 + 0.00
1.08 + 0.00
1.09 + 0.00
1.100 + 0.00
VNNGP
-0.57
0.00
-0.56
0.00
-0.56
0.00
-0,56
0.00
-0.55
0.00
0,66 = 0.00
0.69 0.00
0.73 0.00
0.79
0.00
0.88 0.00
CaGP
-0,22 + 0.00
-0.20
= 0.00
-0.19
+ 0.00
-0.18
+ 0.00
-0.17
0.00
0,94 = 0.00
0.94 + 0.00
0.94 + 0.00
0.95 + 0.00
0.96 +: 0.00
LKGP
28.7 + 0.33
26.9 + 0.41
24.6 + 0.19
23.0 = 0.31
21.2 1: 0.16
15.6 = 0.01
14.2 0.01
12.7 - 0.01
11.3 + 0.00
9.87 + 0.01
Time
mnn
SVGP
150 1 0.08
133
= 0.07
117
+ 0.04
100
+ 0.03
83.3
:: 0.03
150 + 0.02
133 0.13
116 + 0.02
99.9 + 0.05
83.2 + 0.03
VNNGP
158 -: 0.26
140
0 0.17
123 0.33
105 0.16
87.8
4: 0.13
158 0.40
141 + 0.27
123 0.11
105
- 0.13
87.6 - 0.22
CaGP
500 +: 0.15
395
0.06
301 0.12
221 0.01
153
0.01
500 0.02
395 - 0.03
301 + 0.18
221 0.06
153 + 0.07
namely latitude, longitude and days. The main goal of this experiment is to demonstrate the scalability of LKGP large datasets with millions of observations. Additionally, we investigated how increasing the missing data ratio im- pacts predictive performances and runtime requirements. We compare our method to the same sparse and variational methods from the previous experiment.
and CaGP on this problem. We suspect that this is due to the nearest neighbor mechanism working well on these datasets with actual spatial dimensions.
5. Conclusion
We proposed highly scalable exact Gaussian process model for product kernels, which leverages latent Kronecker struc- ture to accelerate computations and reduce memory require- ments for data arranged on partial grid while allowing for missing values. In contrast to existing Gaussian process models with Kronecker structure, our approach deals with missing values by combining projections and iterative linear algebra methods. Empirically, we demonstrated that our method has superior computational scalability compared to the standard iterative methods, and substantially outper- forms sparse variational GPs in terms of prediction quality. Future work could investigate specialized kernels, multi- product generalizations, and heteroskedastic noise models.
For this experiment, we considered latitude and longitude as S and time in days as T. We selected p = 5000 spatial locations uniformly at random together with q = 1000 days of temperature or precipitation measurements, starting on January 1st, 2020, resulting in two datasets of size n S 5M. See Figure 5 for an illustration. We choose ks to be a squared exponential kernel and ky to be the product of squared exponential kernel and periodic kernel to capture seasonal trends. Missing values were selected uniformly at random with ratios of 10%, 20%, .... 50%. Details on the model fitting setup are provided in Appendix C.
Table 2 reports the predictive performances and total run- times of temperature and precipitation prediction. On this large dataset with millions of examples, LKGP clearly out- performs the other sparse and variational methods in both root- rrean-square-ente and negative log-like ligelikki.ood. This is unsurprising, given that LKGP performs inference in the exact GP model while the other methods are subject to a limited number of inducing points, nearest neighbors, or sparse actions. However, due to latent Kronecker structure, LKGP also requires the least amount of time by a large mar- gin. Interestingly, VNNGP consistently outperforms SVGP
Limitations The primary limitation of LKGP is that it employs a product kernel, which assumes that points are only correlated if they are close in both S and T. This is a reasonable assumption in many settings, such as correlated time series data, but will not always be appropriate. The other requirement for our approach is that the data lives on partial grid, which is less restrictive since our model is competitive even if there are a lot of missing values, as our experiments and theory show. Furthermore, if the data does not even live on partial grid, it would be possible to gener- ate an artificial grid, for example via local interpolation.

---

Scalable Gaussian Processes with Latent Kronecker Structure
Impact Statement
Gardner, J., Pleiss, G., Wu, R., Weinberger, K., and Wilson, A. Product Kernel Interpolation for Scalable Gaussian Processes. In International Conference on Artificial Intel- ligence and Statistics, pp. 1407-1416, 2018b.
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
Garnett, R. Bayesian Optimization. Cambridge University Press, 2023.
Acknowledgments
Golub, G. H. and Van Loan, C. F. Matrix Computations. JHU press, 2013.
J. A. Lin was supported by the University of Cambridge Harding Distinguished Postgraduate Scholars Programme, and J. M. Hernandez-Lobato acknowledges support from Turing AI Fellowship under grant EP/VO23756/1. This work was performed using computational resources provided by the Cambridge Tier-2 system operated by the University of Cambridge Research Computing Service.
Hartikainen, J. and Särkkä, S. Kalman filtering and smooth- ing solutions to temporal gaussian process regression models. In 2010 IEEE international workshop on ma- chine learning for signal processing, pp. 379-384. IEEE, 2010.
Hensman, J., Fusi, N., and Lawrence, N. D. Gaussian Processes for Big Data. In Proceedings of the Twenty- Ninth Conference on Uncertainty in Artificial Intellgence, pp. 282-290, 2013.
References
Ambikasaran, S., Foreman- D., Greengard, L., Hogg, D. W., and O'Neil, M. Fast Direct Methods for Gaussian Processes. IEEE Transactions on Pattern Anal- ysis and Machine Intelligence, 38(2):252-25) 265. 2016.
Jankowiak, M.. Pleiss, G., and Gardner, J. Parametric Gaus- sian Process Regressors. In International Conference on Machine Learning, volume 119, pp. 4702-47 12, 2020.
Ament, S. E. and Gomes, C. P. Scalable First-Order Bayesian Optimization via Structured Automatic Dif- ferentiation. In International Conference on Machine Learning, 2022.
Kapoor, S., Finzi, M., Wang, K. A., and Wilson, A. G. G.
SKling on Simplices: Kernel Interpolation on the Per- mutohedral Lattice for Scalable Gaussian Processes. In International Conference on Machine Learning, pp. 5279- 5289, 2021.
Bauer, M., Van der Wilk, M.. and Rasmussen, C. E. Under- standing Probabilistic Sparse Gaussian Process Approxi- mations. In Advances in Neural Information Processing Systems, volume 29, 2016.
Laplante, W., Altamirano, M., Duncan, A. B., Knoblauch, J., and Briol, F.XX Robust and conjugate apatio-temport gaussian processes. In Forty-second International Con- ference on Machine Learning, 2025. URL https: // 1/openreview. net forum? d= 4SWm7 gn.
Bonilla, E. V., Chai, K. M. A., and Williams, K. I. Multi- task Gaussian Process Prediction. In Advances in Neural Information Processing Systems, 2007.
Lin, J. A., Antorán, J., Padhy, S., Janz, D., Hernández- Lobato, J. M., and Terenin, A. Sampling from Gaussian Process Posteriors using Stochastic Gradient Descent. In Advances in Neural Information Processing Systems, 2023.
Deisenroth, M. P. and Rasmussen, C. E. PILCO: A Model- Based and Data-Efficient Approach to Policy Search. In International Conference on Machine Learning, 2011.
Dong, K., Eriksson, D., Nickisch, H., Bindel, D., and Wil- son, A. G. Scalable Log Determinants for Gaussian Pro- cess Kernel Learning. Advances in Neural Information Processing Systems, 30, 2017.
Lin, J. A., Padhy, S., Antorán, J., Tripp, A., Terenin, A.. Szepesvári, C., Hermández-Lobato, J. M., and Janz, D. Stochastic Gradient Descent for Gaussian Processes Done Right. In International Conference on Learning Repre- sentations, 2024a.
Elsken, T., Metzen, J. H., and Hutter, F. Neural Architec- ture Search: A Survey. Journal of Machine Learning Research, 20(55):1-21, 2019.
Lin, J. A., Padhy, S., Mlodozeniec, B., Antorán, J., and Hernández-Lohann J. M. Improving Linear System Solvers for Hyperparameter Optimisation in Iterative Gaussian Processes. In Advances in Neural Information Processing Systems, 2024b.
Frazier, P. A Tutorial on Bayesian Optimization. arXiv e-prints, pp. arXiv:1807.028/ July 2018
Gardner, J., Pleiss, G., Weinberger, K. Q.. Bindel, D., and Wilson, A. G. GPyTorch: Blackbox Matrix-matre Gaus- sian Process Inference with GPU Acceleration. In Ad- vances in Neural Information Processing Systems, 20.
Maddox, W. J., Balandat, M., Wilson, A. G., and Bakshy, E. Bayesian Optimization with High- Digee-iimes Outputs.
10

---

Scalable Gaussian Processes with Latent Kronecker Structure
In Advances in Neural Information Processing Systems, 2021.
Tveito, O. E., Førland, E. J., Heino, R., Hanssen-Bauer, I., Alexandersson, H., Dahlström, B., Drebs, A., Kern- Hansen, C., Jónsson, T., Vaarby-Laursen, E., and West- man, E. Nordic Temperature Maps. DNMI Klima 9/00 KIJMA., 2000.
Pförtner, M., Wenger, J., Cockayne, J., and Hennig, P. Computation-aware kalman filtering and smooth- ing. In Li, Y., Mandt, S., Agrawal, S., and Khan, E. (eds.), Proceedings of The 28th Interna- tional Conference on Artificial Intelligence and Statis- tics, volume 258 of Proceedings of Machine Learn- ing Research, pp. 2071-2079. PMLR, 03-05 May 2025. URL /ttpaaa::////eeeeeee http///prceeee..eeeeeee/oooooese p.lss// v2 v258 2258/pfortne:254 html.
Tveito, O. E., Bjørdal, I., Skjelvág, A. O., and Aune, B. A GIS-based agro-ecoglogical decision system based on gridded climatology. Metoeorl. Appl., 12:57-66 68, 2005.
Wang, K. A., Pleiss, G., Gardner, J. R., Tyree, S., Wein- berger, K. Q.. and Wilson, A. G. Exact Gaussian Pro- cesses a Million Data Points. In Advances in Neural Information Processing Systems, 2019.
J. and Rasmussen, C. E. A unifying view of sparse approximate Gaussian process regression. Journal of Machine Learning Research, 6:1939-1959 2005.
Wenger, J., Wu, K., Hennig. P., Gardner, R., Pleiss, G., and Cunningham, J. P. Computation-Aware Gaussian Processes: Model Selection And Linear-Time Inference. In Advances in Neural Information Processing Systems, 2024.
Rasmussen, C. and Williams, C. I. Gaussian Processes for Machine Learning. MIT Press, 2006.
Wilson, A. and Nickisch, H. Kernel Interpolation for Scal- able Structured Gaussian Processes (KISS-GP). In In- ternational Conference on Machine Learning, pp. 1775- 1784, 2015.
Riis, C., Antunes, F., Hüttel, F., Lima Azevedo, C., and Pereira, F. Bayesian Active Learning with Fully Bayesian Gaussian Processes. In Advances in Neural Information Processing Systems, 2022.
Wilson, A. G., Gilboa, E.,, Nehorai, A., and Cunningham, J. P. Fast Kernel Learning for Multidimensional Pattern Extrapolation. In Advances in Neural Information Pro- cessing Systems, volume 27, 2014.
Ryan, J. P., Ament, S. E., Gomes, C., and Damle, A. The Fast Kernel Transform. In International Conference on Artificial Intelligence and Statistics, 2022.
Saatçi, Y. Scalable Inference Structured Gaussian Pro- cess Models. PhD thesis, University of Cambridge, 2012.
Wilson, J. T., Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. P. Efficiently Sampling Functions from Gaussian Process Posteriors. In International Con- ference on Machine Learning, 2020.
Särkkä, S. and Svensson, L. Bayesian Filtering and Smooth- ing. Cambridge University Press, 2nd edition, 2023. ISBN 9781 1781108926645.
Wilson, J. T., Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. P. Pathwise Conditioning of Gaussian Processes. Journal of Machine Learning Research, 22, 2021.
Sarkka, S.. Solin, A., and Hartikainen, J. Spatiotemporal learning via infnite-dimensional bayesian filtering and smoothing: A look at gaussian process regression through kalman filtering. IEEE Signal Processing Magazine, 30 (4):51-61, 2013.
Wu, K., Wenger, J., Jones, H., Pleiss, G., and Gardner, J. R. Large-Scale Gaussian Processes via Alternating Projection. In International Conference Artificial In- telligence and Statistics, 2024.
doi: 10.1 (tt.///S..013.2442200
Snelson, E. and Ghahramani, Z. Sparse Gaussian Processes using Pseudo-inputs. In Weiss, Y., Schölkopf, B., and Platt, J. (eds.), Advances in Neural Information Process- ing Systems, volume 18. MIT Press, 2005.
Wu, L., Pleiss, G., and Cunningham, Variational Near- est Neighbor Gaussian Process. In International Confer- ence Machine Learning, 2022.
Stegle, O., Lippert, C.,, Mooij, J. M., Lawrence, N., and Borgwardt, K. Efficient inference in matrix-variate Gaus- sian models with iid observation noise. In Advances in Neural Information Processing Systems, 2011.
Zhe, S., Xing, W., and Kirby, R. M. Scalable High-Order Gaussian Process Regression. In International Confer- ence on Artificial Intelligence and Statistics, 2019.
Titsias, M. Variational learning of inducing variables in sparse Gaussian processes. In International Conference on Artificial Intelligence and Statistics, pp. 567-574, 2009.
Zimmer, L.,, Lindauer, M., and Hutter, F. Auto-Py Auto-PyTorch Torch Tabular: Multi-Fidelity MetaLearning for Efficient and Robust AutoDL. IEEE Transactions on Pattern Analysis and Machine Intelligence, 43(9):3079- -0900 2021.

---

Scalable Gaussian Processes with Latent Kronecker Structure
A. Derivation of Asymptotic Break-even Points
To derive the asymptotic break-even points for MVM time and memory usage, as stated in Proposition 3.1, we start by defining the missing ratio 7 n/pq. Alternatively, this can also be written as (1 7)pp. Furthermore, have the properties n. p.q. > 0 and S pq, and thus 7 [0, 1).
The asymptotic time complexity of MVM of
matrix with vector of elements is O For matrix which can be factorized as the Kronecker product of two matrices of sizes p and q respectively, the asymptotic time complexity of MVM is Opp²q + pp²). To find the asymptotic break-even point, we equate these asymptotic time complexities and write n in terms of 7, P. and 9.
n2 = p¹q + pM{; ((1 - )) p2q q- pd
- 7)² p²g + pT: (1 (1-))== 1 1 q P
Since y is the missing ratio, which is between 0 and 1, we conclude that the asymptotic break-even point for time is
Time
For the asymptotic space complexity, we follow the same approach, starting with O and O(p² q²). Again, we set these terms equal to each other and express n in terms of % p. and q.
n² 1² + /² ((1 7))))² 2 1² + q¹,
(1 - 73'2"² p² + 2²
1
(1 (1- q² p
We conclude that the asymptotic break-even point for memory is
7em 1 V
p2
Since the asymptotic time and space complexity kernel evaluations equivalent to the case asymptotic MVM memory usage with (n²)) and O(p² 9²). we conclude that the derivation and result of the second case can be equivalently applied reason about the asymptotic break-even points of kernel evaluation time and memory.
12

---

Scalable Gaussian Processes with Latent Kronecker Structure
B. Additional Experiment Results
Tables to 7 report the predictive performances and total runtimes of our learning curve prediction experiment on all LCBench datasets. See Section 4 for details about the experiment.
Table 3. Predictive performances and total runtimes of learning curve prediction LCBench datasets 1 7. Reported numbers the mean over 10 random seeds. Results of the best model are boldfaced and results of the second best model are underlined.
Model
APSFailure
Amazon
Australian
Fashion
KDDCup09
MiniBooNE
adult
LKGP
1.705 + 0.047 0.273 + 0.005
0.007 + 0.000
0.094 = 0.006
1.609 +: 0.233
0.047 + 0.001
0.053 + 0.001
Train
RMSE
SVGP
3.893 + 0.149
0,402 + 0.007
0.108 = 0.001
0.362 = 0.017
1.775 + 0.247
0.201 + 0.004
0.244 + 0.004
VNNGP
4.920 + 0.194
0.203 + 0.003
0,076 = 0.001
0.341 = 0.174
3.194 = 0.490
0.113 + 0.001
0.131 + 0.002
CaGP
3,972 0.147
0.417 = 0.007
0.119 : 0.001
0.319 = 0.015
1.839 = 0.262
0.219 0.005
0.265 + 0.005
LKGP
2.935 : 0.150
0.386 + 0.007
0.211 = 0.005
0.597 = 0.081
2.080 + 0.360
0.335 + 0.006
0.390 + 0.008
Test RMSE
SVGP
2.716 0.145
0.311 = 0.006
0.195 : 0.001
0.611 = 0.046
2.000
0.330
0.316 0.005
0.351 1 0.005
VNNGP
3.053 1 0.159
0.663 + 0.015
0.593 + 0.005
0.945 = 0.074
2.857 + 0.467
0.705 + 0.013
0.713 1 0.010
CaGP
2.810 + 0.154
0.299
+ 0.006
0.193 + 0.002
0.541 = 0.046
2.058 + 0.331
0.325 + 0.006
0.358 + 0.006
LKGP
1,955 1 0.020
0.122 1 0.019
-2.924 1 0.005
-0,918 1 0.058
1.796 1 0.157
-1.573 1 0.022
1.454 T 0.021
Train
NLI
SVGP
3.098 1 0.100
0.499
1 0.017
-0.746 1 0.010
0.410 = 0,047
1.749 1 0.128
-0.132 1 0.020
0,050 1 0.016
VNNGP
3,005 0.041
0.302
0.005
0.373 = 0.189
2.467
+
0.079
0.033
+ 0.010
CaGP
2.794 0.045
0.524
0.04
0.313 = 0.047
1.927 0.145
0..042
0.141 - 0.018
LKGP
2.349 + 0.047
0.455 + 0.018
-0.798 + 0.013
0.349 = 0.066
1.931
+ 0.195
0.051
+: 0.026
0.146 + 0.009
Test
NL.
SVGP
2.327 : 0.065
0.370
0.018
-0.190 0.012
0.766 = 0.063
1.923
0.172
0.319
0.018
0.401 + 0.016
VNNGP
2.733 + 0.031
0.018
0.753 = 0.007
1.385 = 0.071
2.387
: 0.150
1.046
+ 0.017
1.086 + 0.012
CaGP
2,451 + 0.033
0.355 = 0.017
-0.221 = 0.011
0.723 = 0.123
1.999
0.149
0.321 + 0.021
0.404 + 0.017
LKGP
0.371 + 0.008
0.620 + 0,003
1.680 + 0.014
1.435 = 0.051
0.407
+ 0.013
L,488 0.013
1.539 + 0.011
Time min
SVGP
6,473 + 0.032
6.472 + 0.032
6,475 - 0.033
6.472 = 0.030
6,471
+ 0.032
6,475 + 0.032
6.472 1 0.032
VNNGP
26.34 1 0.443
25.82 1 0.161
25.95 + 0.151
25.84 = 0.171
25.77
+ 0.140
25.89 + 0.161
25.93 1 0.132
CaGP
7.067 + 0.028
7.025
+ 0.021
7.054 + 0.025
7.027 = 0.020
7.084
+ 0.022
7.036 + 0.028
7.070 + 0.026
Table 4. Predictive performances and total runtimes of learning curve prediction LCBench datasets 8 14. Reported numbers are the mcan over 10 random seeds, Results of the best model are boldfaced and results of the second best model are underlined.
Model
airlines
albert
bank
blood
car
christine
cnac-9
LKGP
0.080 0.001
0.055 + 0.002
0.122 = 0.002
0.278 + 0.017
0.017 + 0.000
0.947
0.112
0.005
0.000
Train
RMSF
SVGP
0.293 + 0.004
0.229 + 0.002
0.261 = 0.003
0.463 + 0.009
0.111 + 0.001
1,229 + 0.135
0.075 + 0.001
VNNGP
0.168 + 0.002
0.121 = 0.001
0.135 = 0,001
0.284 : 0.029
0.080 + 0.001
4.379
+ 0,480
0.059 + 0.000
CaGP
0.310 :: 0.004
0.257 : 0.002
0.284 = 0.004
0.511 + 0.011
0.128 + 0.001
1.409
= 0.173
0.074 + 0.001
I.KGP
0.377 + 0.008
0.399 + 0.019
0.371 = 0.008
0.747 + 0.016
0.265 + 0.008
2.829 + 0.487
0.188 + 0.004
Test RMSE
SVGP
0.259 1 0.004
0.308 1 0.005
0.335 0.005
0.540 + 0.013
0.199 + 0.003
2.570
1 0.408
0.166 1 0.002
VNNGP
0.532 1 0.008
0.548 1 0.007
0.712
= 0.009
0.915 1 0.017
0.595
1 0.004
4.600
0.658
0.555 1 0.006
CaGP
0.255 1 0.004
0.321 1 0.004
0.341 - 0.006
0.593 1 0.015
0.204
1 0.002
3.152
1: 0.515
0.176 1 0.003
LKGP
-1.051 1 0.015
-1.122
0.042
-0.664
1 0.016
0.160 1 0.055
-2.528 1. 0.010
1.229
0.200
-3.042 T. 0.003
Train
SVGP
0.181 1 0.012
-0.031
1 0.009
0.093
1 0.013
0.652 1 0.021
-0.713
1 0.006
1.389
1 0.162
-1.113 1 0.010
VNNGP
0.065 1 0.010
-0.116
1 0.005
0.053 - 0.006
0.639
1 0.045
-0.382
+ 0.007
2.788
0.176
-0.620
1: 0.007
CaGP
0.222 0.011
0.095
= 0.007
0.177 = 0.014
0.798
0.021
-0.577
= 0.007
1.709
0.190
-1.084 0.010
LKGP
0.349 + 0.013
0.018 + 0.025
0.231 = 0.014
1.047 + 0,029
-0.537 + 0.018
2.425
+ 0.368
-0.863 + 0.048
Test
NLL
SVGP
0.137 - 0.011
0.215 0.014
0.329 0.015
0.781 - 0.020
-0.173 - 0.017
2.140
::400:00 0.278
-0.349 + 0.025
VNNGP
0.736 0.012
0.772
1 0.007
1.072
= 0.013
1.365
+ 0.019
0.763
0.007
2.852
0.220
0.587
1 0.010
CaGP
0.137 - 0.011
0.266 - 0.010
0.344 = 0.015
0.916
+ 0.023
-0.174
+ 0.013
2.538
- 0.332
-0.242 - 0.034
LKGP
1.182 1 0.008
1.937 1 0.044
0.908 = 0.006
0.487 1 0.003
1.341 1 0.007
0.459
0.063
1.531 0.010
Time
in min
SVGP
6.479 1 0.035
6.472 1 0.032
6.481 = 0.036
6.479
1 0.031
6,474 1 0.030
6.469
0.031
6.495 1 0.033
VNNGP
25.92 1 0.170
25.87
1 0.140
25.81 = 0.156
25.78
1 0.151
25.88
1 0.143
25.83
26.39 1 0.471
CaGP
7.129 + 0.044
7.034 + 0.020
7.019 = 0.028
7.033
+ 0.023
7.042 + 0.023
7.020
7.032 + 0.019
13

---

Scalable Gaussian Processes with Latent Kronecker Structure
Table Predictive performances and total runtimes of learning curve prediction on LCBench datasets 15 21. Reported numbers are the mean over 10 random seeds. Results of the best model are boldfaced and results of the second best model are underlined.
Model
connect-4
covertype
credit-g
dionis
fabert
helena
higgs
LKGP
0.149 1: 0.005
0.061 +: 0.001
0.047 + 0.000
0.008 + 0.000
0.027 + 0.001
0.012 1 0.000
0.039 + 0.001
Train RMSE
SVGP
0.239 1 0.009
0.240 + 0.003
0.133 + 0.001
0.114 + 0.003
0.116 + 0.001
0.108 + 0.001
0.200 - 0.003
VNNGP
0.117
0.003
0.158 1 0.003
0.123 1 0.002
0.088
0.002
0.072 1 0.001
0.088 1 0.001
0.130 0.002
CaGP
0.265 + 0.011
0.241 = 0.004
0.174 0.002
0.123 0.003
0.119 1119 0.001
0.122 - 0.001
0.208 = 0.002
LKGP
0.400
0.021
0.351 1 0.008
0.315 1 0.008
0.249
0.008
0.276 1 0.007
0.257 1 0.030
0.394 0.027
Test RMSE
SVGP
0.363 1: 0.021
0.294 1: 0.005
0.257 1 0.004
0.264 + 0.010
0.230 1 0.004
0.189 1 0.004
0.285 1 0.005
VNNGP
0.722 + 0.034
0.677 = 0.011
0,585
= 0,002
0.646
= 0.015
0.613 = 0.00
0.619 0.010
0.642 = 0.009
CaGP
0.388 + 0.022
0.281 = 0.004
0.270 0.002
0.009 0.254
0.243 + 0.003
0.205 + 0.005
0.296 0.004
LKGP
-0.482 + 0.034
-1.311 + 0.017
-1,611 + 0.010
-2,941 + 0.010
-2.115 + 0.052
-2.727 +: 0.019
-1.718 = 0.018
Train
NLL
SVGP
-0.027 + 0.038
0.025 = 0.013
-0.515 = 0.006
-0.687 = 0.021
-0.672 = 0.010
-0.745 + 0.011
-0.150 = 0.012
VNNGP
-0.101 + 0.024
0.104 + 0.009
-0.132 = 0.010
-0.334 = 0.016
-0,429 0.006
-0.323
0.008
-0.056 = 0.009
CaGP
0.086 - 0.040
0.034 + 0.017
-0,265 = 0.009
-0.557 + 0.022
-0.609 + 0.011
-0.614 + 0.012
-0.103 = 0.012
LKGP
0.132
0.048
-0.057
0.019
-0.202
0.018
-0.551 + 0.033
-0.329
0.060
-0.658
0.054
-0.328 = 0.023
Test
NLL
SVGP
0.267 1: 0.053
0.199 + 0.016
0.018
+ 0.020
0.209 + 0.060
-0.044 1 0.024
-0.215 1: 0.028
0.160 + 0.018
VNNGP
0.958 1: 0.047
1011 + 0.016
0.801
+ 0.005
0.880 + 0.024
0,771 - 0.009
0.826
1 0.014
0.891 + 0.013
CaGP
0.421 + 0.061
0.142 + 0.016
=
0.182 0.062
0.018 + 0.025
-0.155 0.028
0.171 = 0.013
LKGP
0.679 1 0.003
1.728 1 0.010
0.925 1 0.009
1.622 1 0.020
1,432 1 0.018
2.066
0.022
2.447 + 0.022
Time
in min
SVGP
6.480 1 0.032
6.474 1 0.032
6,493 1 0.036
6.476 0.031
6,493 1 0.031
6.476
0.035
6.492 0.033
VNNGP
25.89 0.181
25.85 + 0.158
25.84 0.151
26.34
0.512
25.91 = 0.155
25.86 + 0.167
25.87 0.181
CaGP
7.081 + 0.029
7.015 = 0.020
7.024 = 0.021
7.041 = 0.023
7.015 + 0.026
7.018 + 0.022
7.046 = 0.024
Table Predictive performances and total runtimes of learning curve prediction on LCBench datasets 22 28. Reported numbers are the mean over 10 random seeds. Results of the best model are boldfaced and results of the second best model are underlined.
Model
jannis
jasmine
jungle
kcl
kr-vs-kp
mfcat-
nomao
I.GPP
0.074 + 0.001
0.027 + 0.001
0.044 + 0.001
0.282 + 0.006
0.008 + 0.000
0.006 + 0.000
0.027 + 0.001
Train
RMSE
SVGP
0.224 : 0.004
0.146 = 0.002
0.176 = 0.003
0.408
0.008
0.101 = 0.001
0.099 + 0.001
0.245 + 0.006
VNNGP
0.128 + 0.001
0.088
0.001
0.109 = 0.001
0.185 0.002
0.067 = 0.001
0.079
0.001
0.133 : 0.002
CaGP
0.243 + 0.004
0.149 = 0.001
0,201 = 0.004
0.459 = 0.008
0.105 + 0.002
0.100 = 0.001
0.245 : 0.007
LKGP
0.375 + 0.014
0.290
0.008
0.314 = 0.005
0.497 0.011
0.261 +: 0.005
0.165
0.003
0.324 + 0.015
Test RMSE
SVGP
0.314 1: 0.005
0.264 1 0.004
0.287 = 0.006
0.442 - 0.013
0.232 1 0.003
0.160 1 0.001
0.302 1 0.009
VNNGP
0.675 0.013
0.627 + 0.003
0.683
- 0.012
0.726
+ 0.007
0.591
+ 0.006
0.563 + 0.003
0.702 + 0.018
CaGP
+
0.280 + 0.003
0.312 - 0.005
0.428
0.008
0.239 + 0.004
0.165 + 0.001
0.288 0.008
LKGP
-1.138 1 0.020
-2.094 1 0.031
-1.647
- 0.017
0.163 + 0.023
-2.333 1 0.005
-2.930 1 0.005
-2.035 1 0.034
Trann
NLL
SVGP
-0.055 + 0.016
-0.461
= 0.009
-0.266 = 0.016
0.495
0.021
-0.000
- 0.014
-0.836
+ 0.010
0.048 + 0.025
VNNGP
-0.023 + 0.009
-0.292
- 0.007
-0.127 = 0.010
0.264
0.009
-0.489
+ 0.008
-0.418
+ 0.009
0.032 - 0.012
CaGP
0.038 + 0.016
-0.414
= 0.009
-0.139 = 0.017
0.635 - 0.017
-0.718
- 0.013
-08002
+ 0.009
0.064 - 0.028
LKGP
0.148 + 0.025
-0.382 + 0.018
-0.194 = 0.011
0.638
0.020
-0.633 - 0.029
-1.064 + 0.015
0.104 + 0.034
Test
NLL
SVGP
0.272 + 0.017
0.034
+ 0.012
0.184 = 0.020
0,590
0.009
0.065 : 0.032
.0.3% + 0.013
0.285 = 0.031
VNNGP
0.995 + 0.017
0,851
+ 0.006
0.974 - 0.015
1.105 - 0.009
0.717
+ 0.013
0.711 + 0.006
1.130 = 0.027
CaGP
0.286 + 0.017
0.131 + 0.015
0.259 + 0.017
0.011 0.616
0.088 + 0.031
-0.388 + 0.010
0.229 + 0.030
LKGP
1.327 1 0.008
1.867
+ 0.016
1.881 - 0.011
0.576 + 0.009
1.493 + 0.017
1.942 + 0.018
2.183 1 0.020
in Time
SVGP
6.478 1 0.032
6.495
+ 0.030
6.482 - 0.031
6.491
0.033
6.500
+ 0.030
6.491 + 0.030
6.487 1 0.032
VNNGP
26.40 + 0.449
25.91
25.89 0.162
25.89 0.143
25.85
+ 0.154
25.86 + 0.167
CaGP
7.009 + 0.025
7.019
+ 0.021
7.034 = 0.029
0.024 7.016
7.024 + 0.019
7.034
0.030
7.033 + 0.023
14

---

Scalable Gaussian Processes with Latent Kronecker Structure
Table 7. Predictive performances and total runtimes of learning curve prediction I.Ceench datasets 29 35. Reported numbers the
mean over 10 random seeds. Results of the best model are boldfaced and results of the second best model are underlined.
Model
numerai28.6
phoneme
segment
shuttle
sylvine
vehicle
volkert
LKGP
1.007 + 0.156
0.062 + 0.001
0.006 + 0.000
1.172 = 0.265
0.008 + 0.000
0.017 + 0.000
0.043 + 0.001
Train
RMSF
SVGP
1.669 + 0.188
0.208 + 0.002
0.086 = 0.001
1.316 0.286
0.109 + 0.001
0.077 = 0.001
0.147 + 0.016
VNNGP
5,900 1 0.618
0.132 + 0.001
0.072 + 0.001
2.611 = 0.620
0.075 + 0.001
0.062 1 0.001
0.105 + 0.001
CaGP
1.327 +: 0.157
0.238 +: 0.002
0.091 + 0.001
1.361 = 0.303
0.118 +: 0.002
0.089 + 0.001
0.171 0.003
LKGP
0.831 + 0.099
0.349 +: 0.007
0.152 + 0.003
1.597
= 0.368
0.253 +: 0.007
0.203 + 0.007
0.260 + 0.007
Test RMSF
SVGP
1.104 + 0.120
0.285 + 0.005
0.145 + 0.002
1.568
= 0.371
0.228 + 0.003
0.169 + 0.002
0.200 + 0.021
VNNGP
1.686 + 0.180
0.676 + 0.005
0.568 + 0.005
2.429
= 0.585
0.595 + 0.006
0.576 + 0.006
0.617 + 0.010
CaGP
0.882 + 0.096
0.295 + 0.004
0.153 + 0.002
1.551
0.356
0.241 + 0.003
0.186 + 0.003
0.228 + 0.004
I.KPP
1.393 + 0.140
-1.323 + 0.010
-2.943 + 0.005
1,348
0.219
-2.900 + 0.008
-2.568 + 0.012
-1.668 + 0.029
Train
NLI
SVGP
1.890 + 0.128
-0.112 + 0.007
-0.977 + 0.007
1,435
0.205
-0.728 + 0.011
-10883 + 0.009
-0.321 + 0.036
VNNGP
3.146 1 0.099
0.003 + 0.005
-0.456 + 0.006
2.130
= 0.221
-0,414 0.009
-0.567
0.007
-0.193 1 0.009
CaGP
1.668 1: 0.105
0.015 1 0.007
-0.900
0.009
1.490
= 0.208
-0.624 1 0.014
-0.942
0.008
-0.303 1 0.017
LKGP
1.372 + 0.099
-0.050
0.017
-1.214 = 0.020
1.533 = 0.2247
-0.641
+ 0.026
-0.923
0.034
-0.266 + 0.030
Test
SVGP
1.595 1 0.075
0.168 1 0.014
-0.496
0.014
1.573
= 0.246
-0.021
1 0.026
-0.381
0.019
-0.077 1 0.022
VNNGP
2.734 1 0.100
0.991 1 0.007
0.712
1 0.009
2.079
= 0.221
0.765
1 0.009
0.624
1 0.010
0.883 1 0.016
CaGP
1.519 + 0.099
0.199
+ 0.011
-0,463
+ 0.013
1.577
= 0.219
0.049
+ 0.020
-0.229 1: 0.019
-0.084 1 0.018
LKGP
0.410 + 0.007
1.136 + 0.005
2.012 + 0.018
0.471 = 0.023
1.629
+ 0.018
1.368 + 0.009
1.821 + 0.018
Time mim
SVGP
6.477 + 0.030
6.478 + 0.029
6.475 + 0.032
6.471 = 0.031
6,474 + 0.032
6.471 + 0.031
5.827 + 0.615
VNNGP
25.95 + 0.162
25.82
+ 0.181
26.31
+ 0.464
25.94 = 0.171
25.85
+ 0.168
26.11 + 0.174
26.33 0.442
CaGP
7,044 + 0.026
7.023
+ 0.016
7.023
= 0.020
7.033 = 0.032
7,014 0.024
7.029 + 0.025
7.035 1 0.020
C. Experiment Details
Inverse Dynamics Prediction For both methods, observation noise and kernel hyperparameters were initialized with GPyTorch default values, and inferred by using Adam with a learning rate of 0.1 to maximize the marginal likelihood for 50 iterations. Additionally, both methods use conjugate gradients with a relative residual norm tolerance of 0.01 as iterative linear system solver.
Learning Curve Prediction For methods, observation noise and kernel hyperparameters were initialized with GPyTorch default values, and optimized using Adam. LKGP was trained for 100 iterations using a learning rate of 0.1, conjugate gradients with a relative residual norm tolerance of 0.01 and a pivoted Cholesky preconditioner of rank 100. SVGP was trained for 30 epochs using a learning rate of 0.01, a batch size of 1000, and 10000 inducing points, which were initialized at random training data examples. VNNGP was trained for 1000 epochs using a learning rate of 0.01, a batch size of 1000, inducing points placed at every training example, and 256 nearest neighbors. CaGP was trained for 1000 epochs using a learning rate of 0.1 and 512 actions.
Climate Data with Missing Values For all methods, observation noise and kernel hyperparameters were initialized with GPyTorch default values and optimized using Adam. LKGP was trained for 100 iterations using a learning rate of 0.1, and conjugate gradients with a relative residual norm tolerance of 0.01 and a pivoted Cholesky preconditioner of rank 100. SVGP was trained for 5 epochs using learning rate of 0.001, a batch size of 1000, and 10000 inducing points, which were initialized at random training data examples. VNNGP was trained for 50 epochs using a learning rate of 0.001, a batch size of 1000, 500000 inducing points placed at random training data examples, and 256 nearest neighbors. CaGP was trained for 50 epochs using a learning rate of 0.1 and 256 actions.
15


Reviews:
[{"review_id": "pA0zXvNrc4", "reviewer": "Reviewer_Zcew", "summary": "The paper improves on the scalability of Kronecker product large scale Gaussian processes (GPs). Specifically,  Kronecker product GPs cannot easily deal with missing observation outputs in the gridded data. The paper proposes a trick based on a projection operation that computes the full Kronecker product kernel matrix and projects it to the sub-matrix  associated with the observed outputs. Somehow this operation still allows for faster matrix-vector multiplications and it can be combined with conjugate gradients solvers to achieve \nfast computation for training the GP model or doing inference (e.g., exact sampling from the GP posterior). Three examples mainly in spatio-temporal modeling show that this method can outperform other sparse GP methods in large datasets.", "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": "Dear Reviewer Zcew,\n\nThank you for your time and effort spent on reviewing our paper. In the following, we provide answers to your specific questions.\n\n---\n\n> (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?\n\n> (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.\n\nGreat 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:\n\n**Inverse Dynamics Prediction**\n\nThere 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,\n\n(pos_1, vel_1, acc_1, ..., pos_7, vel_7, acc_7, joint_id) --> torque,\n\nwhere joint_id ranges from 1 to 7 and specifies which joint's torque the output corresponds to. Here, we partition the input dimensions as\n\n(pos_1, vel_1, acc_1, ..., pos_7, vel_7, acc_7) x (joint_id,),\n\nwhere 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**.\n\n**Learning Curve Prediction**\n\nWe view the problem as\n\n(x_1, ..., x_k, t) --> learning curve value,\n\nand the Cartesian product is\n\n(x_1, ..., x_k) x (t,),\n\nwhere (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.\n\n**Climate Data**\n\nWe model the problem as\n\n(lat, long, t) --> temperature / precipitation,\n\nand partition the input dimensions as\n\n(lat, long) x (t,),\n\nwhich 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**.\n\n---\n\n> (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.\n\nThank 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.\n\nIn 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.\n\nTo be more explicit, the derivative of the log marginal likelihood is\n\n$\\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})$,\n\nwhere $K$ is the kernel matrix, $y$ are the training labels, and $\\theta$ are the kernel hyperparameters.\n\nUsing Hutchinson's trace estimator, the trace term is approximated as\n\n$\\frac{1}{s}\\sum_{j=1}^s z_j^T K^{-1} \\frac{\\partial K}{\\partial \\theta} z_j$,\n\nwhere $z_j$ are random probe vectors with $\\mathbb{E} [z_j z_j^T] = I$.\n\nTo circumvent the explicit inverse of $K$, we instead calculate $K^{-1} y$ and $z_j^T K^{-1}$ by solving $s + 1$ linear systems,\n\n$K [v_y, v_1, ..., v_s ] = [y, z_1, ..., z_s]$,\n\nusing conjugate gradients, which only relies on MVMs with $K$ and leverages the Kronecker structure in $K$.\n\n---\n\nWe 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.", "claims_and_evidence": "The paper shows convincing evidence that the proposed method can work in low dimensional datasets. \nIt is unclear if this method is a general purpose GP approximation for high-dimensional inputs.", "methods_and_evaluation_criteria": "Yes, the paper uses some standard real-world regression datasets, such as spatio-temporal datasets, which are very \nappropriate to demonstrate the proposed method. Also Root Mean Square Error and Test Log Likelihood  are the standard \nperformance scores for such regression problems.", "theoretical_claims": "I checked the derivations and they appear to be correct.", "experimental_designs_or_analyses": "The experimental analysis is valid.", "relation_to_broader_scientific_literature": "The paper cites very accurately all the related work in the area.", "essential_references_not_discussed": "No.", "other_strengths_and_weaknesses": "Some possible weakness is that the paper does not discuss much the scalability with respect to the dimensionality of the inputs; \nsee \"Questions for Authors\" below.", "other_comments_or_suggestions": "The paper is very well written. I didn't notice any typos.", "questions_for_authors": "The paper is very well written and the proposed method is based on a nice trick using this projection operation illustrated in Figure 1. \nSome questions I have are the following: \n\n(i) It is unclear if the method is still applicable in high-dimensional spaces since we know that for Kronecker-product formulations  \nthe 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? \n\n(ii) Related to the above it is unclear how the authors apply their method to the Sarcos dataset. The Sarcos dataset (see e.g. https://gaussianprocess.org/gpml/data/) has input dimensionality d=21, so did you run using this dimensionality \nor 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. \n\n(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\nand particularly how do you deal with the derivatives of the log determinant term.", "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": "xKJ6YExwRV", "reviewer": "Reviewer_F5JD", "summary": "In multitask or spatio-temporal Gaussian process regression, scalability issues lead users to adopt a Kronecker, also described as a factored, structure in the covariance between tasks, space, and time. However, this computational gain requires that observations be available for all tasks and inputs, or across all spatial and temporal locations.  \n\nIn this paper, the authors re-express the product-factored kernel in a setting where only some tasks or locations are observed, formulating it as a down-projection of the complete Kronecker product. They implement this down-projection using zero-padding and slice indexing, ensuring low computational overhead when employing matrix-vector product-based forms of Gaussian process inference.", "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": "Dear Reviewer F5JD,\n\nThank you for your time and effort spent on reviewing our paper. In the following, we address your suggestions, concerns, and questions.\n\n---\n\n> I couldn’t find in the paper whether multiple runs were performed with different selections of missing data for the inverse dynamics experiment.\n> - Report the mean and standard deviation over multiple runs of the experiments with different random seeds to ensure more robust results.\n\nThank 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.\n\n**Inverse Dynamics Prediction**\n\nWe 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.\n\nThe relative standard errors (relative to corresponding metric value, averaged across missing ratios) are small, demonstrating stability across replications with different random seeds.\n\n| | Time | Memory | Test RMSE | Test NLL |\n|---|:---:|:---:|:---:|:---:|\n| Iterative GP | 0.85% | 0.00% | 1.27% | 0.33% |\n| LKGP | 0.58% | 0.00% | 1.26% | 0.32% |\n\n**Learning Curve Prediction**\n\nWe 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.\n\nAkin to the above, the relative standard errors (averaged across datasets) are small.\n\n| | Test NLL | Test RMSE | Time |\n|---|:---:|:---:|:---:|\n| LKGP | 1.89% | 1.61% | 0.49% |\n| SVGP | 1.24% | 1.29% | 0.25% |\n| VNNGP | 0.71% | 1.09% | 0.26% |\n| CaGP | 3.18% | 1.19% | 0.11% |\n\n**Climate Data**\n\nDue 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.\n\nAgain, the relative standard errors (averaged across datasets and missing ratios) are small.\n\n| | Test NLL | Test RMSE | Time |\n|---|:---:|:---:|:---:|\n| LKGP | 0.76% | 0.63% | 1.32% |\n| SVGP | 1.11% | 0.29% | 0.09% |\n| VNNGP | 0.13% | 0.27% | 0.40% |\n| CaGP | 1.65% | 0.79% | 0.06% |\n\n---\n\n> 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.\n\nThank 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.\n\nRegarding 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).\n\n---\n\nWe 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.", "claims_and_evidence": "Yes, the authors claim that their method should have the same predictive performance as standard Gaussian process methods while using less memory and compute time, up until the data points are no longer distributed on a grid, according to the asymptotic break-even point. These claims are substantiated with an extensive experimental comparison.", "methods_and_evaluation_criteria": "Yes, the evaluation criteria make sense. The authors explore different types of datasets (kinematics, training curves, geospatial) and evaluate scenarios with and without uniform gridding.", "theoretical_claims": "Yes, the derivation in Appendix A appears correct.", "experimental_designs_or_analyses": "I couldn’t find in the paper whether multiple runs were performed with different selections of missing data for the inverse dynamics experiment.  \n\nAdditionally, 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.  \n\nI recommend that the authors:  \n\n- Clarify the initialization of the variational parameters in SVGP.  \n- Report the mean and standard deviation over multiple runs of the experiments with different random seeds to ensure more robust results.  \n- Randomly initialize the Gaussian processes using sensible hyperparameters for each task.", "relation_to_broader_scientific_literature": "This paper fits well within the broader topic of GP scalability, where the scalability of product-factored kernels has been a common subject of investigation. In particular, this paper contributes to the line of research focusing on GP inference methods that rely solely on matrix multiplications, avoiding Cholesky factorizations and matrix inverses. Within this context, the additional speed-up in matrix multiplication demonstrated by the paper is appreciated.", "essential_references_not_discussed": "I believe all relevant works were cited.", "other_strengths_and_weaknesses": "None.", "other_comments_or_suggestions": "Overall, I found the paper easy to read, and the authors have clearly presented and analyzed their approach.", "questions_for_authors": "No additional questions or than the ones raised in the methodology section.", "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": "jJUYTj3uWq", "reviewer": "Reviewer_UJyJ", "summary": "This work proposes an exact latent Kronecker structure for Gaussian Processes, by expressing the covariance matrix of the observed values (not restricted to a grid structure) as the projection of a latent Kronecker product matrix (defined on a grid structure).\n\nNamely, it allows defining a Kronecker product-based matrix structure for the GP covariance, but accommodates observations missing from the assumed, underlying grid.\n\nThe combination of such latent Kronecker structure with iterative linear system solvers enables inference over the exact GP model (not an approximation to it), at reduced computational and memory complexity.\n\nThe authors demonstrate that the proposed method outperforms state-of-the-art sparse and variational GPs in real-world datasets with many (millions of) examples at reduced complexity.\n\n# After rebuttal\nThe authors have provided an informative rebuttal to all.\n\nHence, I argue for acceptance of this submission, and I would encourage the authors to incorporate key aspects of their rebuttal to the camera-ready version.", "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": "Dear Reviewer UJyJ,\n\nThank you for your time and effort spent on reviewing our paper. In the following, we provide answers to your specific questions.\n\n---\n\n> Can the authors elaborate on why results in Table 1 showcase improved predictive RMSE performance of SVGP, yet not best negative log-likelihood?\n> - 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\n\nGreat 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.\n\nIn 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.\n\nFurthermore, we can compare the performance between train and test sets:\n\n| Average Rank | Train RMSE | Test RMSE | Train NLL | Test NLL |\n|:---:|:---:|:---:|:---:|:---:|\n| LKGP | **1.17** | 2.60 | **1.03** | **1.26** |\n| SVGP | 2.77 | **1.66** | 2.40 | 2.54 |\n\n**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.\n\nWe will add a summary of this discussion to the camera-ready version of our manuscript.\n\n---\n\n> 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.\n> - Can the authors elaborate further on these tricks/techniques, to help unfamiliar readers and for the completeness of the manuscript?\n\nThank 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.\n\n**Latent Kernel Projections**\n\nIn 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.\n\nFor example,\n\n$P = \\begin{bmatrix} 1 & 0 & 0 & 0 \\\\\\ 0 & 0 & 1 & 0\\end{bmatrix}$,\n\nis a submatrix of a 4 x 4 identity matrix without its 2nd and 4th row.\n\nMultiplying $P$ with a vector removes the 2nd and 4th element.\n\n$P \\begin{bmatrix} 1 \\\\\\ 2 \\\\\\ 3 \\\\\\ 4 \\end{bmatrix} = \\begin{bmatrix} 1 \\\\\\ 3 \\end{bmatrix}$\n\nMultipliying $P^T$ with a vector inserts zeros as 2nd and 4th element.\n\n$P^T \\begin{bmatrix} 1 \\\\\\ 3 \\end{bmatrix} = \\begin{bmatrix} 1 \\\\\\ 0 \\\\\\ 3 \\\\\\ 0\\end{bmatrix}$\n\nTo 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].\n\nWe can avoid explicit matrix multiplications via\n\nx = [1, 2, 3, 4]\n\nPx = x[idx]\n\nand\n\nx = [1, 3]\n\nPTx = zeros(4)\n\nPTx[idx] = x\n\n**Lazy Kernel Evaluation**\n\nA $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.\n\nfor i in range(n): result[i] = k_i @ x \n\nHere, k_i is the i-th row of the kernel matrix, and we rematerialize its values whenever needed to save memory.\n\n---\n\nWe 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.", "claims_and_evidence": "All their claims are well supported:\n- The Latent Kronecker Structure model is sound, and its efficiency is both theoretically (analysis in Proposition 3.1) and empirically (results in Figure 3) demonstrated.\n\n- The computational efficiency and inference quality of the proposed framework are clearly evaluated with experiments in challenging real-word experiments, with many samples.", "methods_and_evaluation_criteria": "Absolutely: the experiments are set to evaluate different, yet important criteria of the proposed work (see Experimental Design or Analysis section)", "theoretical_claims": "The proposed latent Kronecker structure that represents the joint covariance matrix of observed values as a product between projection matrices and a latent Kronecker product is correct to the best of my knowledge.", "experimental_designs_or_analyses": "- The Inverse Dynamics Prediction experiment, with results in Figure 3, is well set-up to fairly compare the matrix multiplication costs of LKGP to baseline iterative methods.\n\n- The Learning Curve Prediction experiment, with results in Figure 4 and Table 1, aims at evaluating performance differences (measured via RMSE, log-likelihood and run-time) under realistic, non-uniform missing values of LKGP and sparse GP alternatives.\n\n- The Climate Data with Missing Values experiment with results in Table 2, showcases the scalability of LKGP to real-life, large datasets with millions of observations: it provides improved performance and reduced computational cost.", "relation_to_broader_scientific_literature": "The key contribution of the paper is well defined and presented within the GP inference literature: this work is a highly scalable exact Gaussian process model for product kernels, made possible by leveraging latent Kronecker structure.", "essential_references_not_discussed": "Section 2 of this work presents a very precise and well written overview (plus necessary details) with key references of the Gaussian Process model and the existing different inference frameworks (SparseGP Vs Iterative solvers) needed to situate and understand this work.", "other_strengths_and_weaknesses": "Strengths\n\n- The proposed latent Kronecker structure to enable efficient inference in the context of data on grids with missing values is clever and novel.\n    - Even if efficient matrix factorization is not possible, the authors figure out how to leverage the projected Kronecker structure for fast matrix multiplication.\n\n- The proposed inference methodology is defined over the exact GP (with assumed covariance structure), rather than over an approximation to the assumed GP prior.\n\n- The proposed method is shown to scale to a large number of training examples, at reduced computational cost.\n    - The authors clearly explain and present both the improved asymptotic time and space complexities of matrix\nmultiplication, as well as the reduced number of kernel evaluations.\n\nWeakness\n\n- The main, yet well discussed and acknowledged weakness of this work is its limitation to GP models defined on a product space, via product kernels.", "other_comments_or_suggestions": "None", "questions_for_authors": "- Can the authors elaborate on why results in Table 1 showcase improved predictive RMSE performance of SVGP, yet not best negative log-likelihood?\n    - 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\n\n- 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.\n    - Can the authors elaborate further on these tricks/techniques, to help unfamiliar readers and for the completeness of the manuscript?", "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": "vTlc7BNAdN", "reviewer": "Reviewer_1bD1", "summary": "This paper considers a Gaussian process (GP) regression with data points on a Cartesian grid and proposes a new way of constructing GP gram matrices to deal with missing values. The proposed method represents a gram matrix with missing values as a projection of a latent exact gram matrix, enabling efficient GP inference by combining the Kronecker structure of the gram matrix with iterative methods. The experiments on synthetic and real-world data confirmed the model's validity.", "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": "Dear Reviewer 1bD1,\n\nThank you for your time and effort spent on reviewing our paper. In the following, we provide responses to your specific concerns and questions.\n\n---\n\n> 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\n> - [Smola et al., Kernel methods for missing variables, 2005]\n> - [Jafrasteh et al., Gaussian processes for missing value imputation, 2023]\n> - [Imani et al, Nested Gaussian process modeling and imputation of high-dimensional incomplete data under uncertainty, 2019]\n\nThank 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**.\n\n**(1) Missing Input Features**\n\nSuppose we are trying to model a latent function with two input dimensions and one output dimension:\n\n(x_1, x_2) --> y\n\nThe 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**.\n\n**(2) Missing Output Observations in a Cartesian Product Space**\n\nOur 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.\n\nSuppose we are modeling a latent function\n\n(x_1, x_2) --> y\n\nwhere x_1 can be a, b, or c, and x_2 can be 1, 2, or 3, such that the whole input space is\n\n| (a, 1) | (a, 2) | (a, 3) |\n|---|---|---|\n| (b, 1) | (b, 2) | (b, 3) |\n| (c, 1) | (c, 2) | (c, 3) |\n\nIn 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).\n\nIn 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)**.\n\nSmola 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:\n\n(x, idx) --> y[idx]\n\nIn 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.\n\n---\n\n> 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?\n\nSparse 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.\n\n---\n\nWe 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.", "claims_and_evidence": "The relationship to related works regarding the handling of missing values is not clearly presented. The discussion of related works appears to focus primarily on GP models without missing values.", "methods_and_evaluation_criteria": "The proposed methods and evaluation criteria make sense.", "theoretical_claims": "Theoretical claims seem to be correct.", "experimental_designs_or_analyses": "No issues specified.", "relation_to_broader_scientific_literature": "This paper provides an easy-to-implement algorithm for GP models with missing values, which could be beneficial in many research fields because GP models have been widely utilized.", "essential_references_not_discussed": "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], and [Imani et al, Nested Gaussian process modeling and imputation of high-dimensional incomplete data under uncertainty, 2019]?", "other_strengths_and_weaknesses": "Strength\n- The idea of the paper is clearly presented, and the proposed method is easy-to-implement.\n- The validity of the proposed model was evaluated on a lot of real-world data.\n\nWeakness\n- See above.", "other_comments_or_suggestions": "No other comments.", "questions_for_authors": "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?", "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: Thanks to the authors for their interesting and high quality submission.  Reviewers were uniformly positive, and agree that this contribution fills a gap in the scalable Gaussian Process literature.  Reviewers also strongly encourage the authors to incorporate the additional exposition and experiments from the rebuttal period into the final draft of the manuscript.
Decision: Accept (poster)