Log-Linear-Time Gaussian Processes Using Binary Tree Kernels
Authors: Michael K. Cohen
Sam Daulton
Michael A Osborne
Venue: NeurIPS 2022 Accept (2022)
Keywords: Probabilistic methods
Gaussian processes
Scalable Gaussian processes

Abstract: Gaussian processes (GPs) produce good probabilistic models of functions, but most GP kernels require $O((n+m)n^2)$ time, where $n$ is the number of data points and $m$ the number of predictive locations. We present a new kernel that allows for Gaussian process regression in $O((n+m)\log(n+m))$ time. Our "binary tree" kernel places all data points on the leaves of a binary tree, with the kernel depending only on the depth of the deepest common ancestor. We can store the resulting kernel matrix in $O(n)$ space in $O(n \log n)$ time, as a sum of sparse rank-one matrices, and approximately invert the kernel matrix in $O(n)$ time. Sparse GP methods also offer linear run time, but they predict less well than higher dimensional kernels. On a classic suite of regression tasks, we compare our kernel against Mat\'ern, sparse, and sparse variational kernels. The binary tree GP assigns the highest likelihood to the test data on a plurality of datasets, usually achieves lower mean squared error than the sparse methods, and often ties or beats the Mat\'ern GP. On large datasets, the binary tree GP is fastest, and much faster than a Mat\'ern GP.

Log-Linear-Time Gaussian Processes Using Binary Tree Kernels
Michael K. Cohen University of Oxford michael. cohendeng. ox. ax.
Samuel Daulton University of Oxford, Meta sdaul tonmetan com
Michael A. Osborne University of Oxford mosb@robots ox, ax. uk
Abstract
Gaussian processes (GPs) produce good probabilistic models of functions, but most GP kernels require ( + m)n²²) time, where n is the number of data points and m the number of predictive locations. We present a new kernel that allows for Gaussian process regression in O((n + m) log(n + m)) time. Our "binary tree" kernel places all data points on the leaves a binary tree, with the kernel depending only on the depth of the deepest common ancestor. We can store the resulting kernel matrix in (n) space in On log time, as a sum of sparse rank-one matrices, and approximately invert the kernel matrix in (n)) time. Sparse GP methods also offer linear run time, but they predict less well than higher dimensional kernels. On classic suite of regression tasks, compare our kernel against Matérn, sparse, and sparse variational kernels. The binary tree GP assigns the highest likelihood to the test data on a plurality of datasets, usually achieves lower mean squared error than the sparse methods, and often ties or beats the Matérn GP. On large datasets, the binary tree GP fastest, and much faster than Matérn GP.
Introduction
Gaussian processes (GPs) can be used to perform regression with high-quality uncertainty estimates. but they are slow. Naïvely, GP regression requires Onn" - n²m) computation time and Onz²) computation space when predicting at locations given n data points [28]. A kernel matrix of size x n must be inverted (or Cholesky decomposed), and then m matrix-vector multiplications must be done with that inverse matrix m linear solves with the Cholesky factors). few methods that we will discuss later achieve (nn²m) time complexity [25, 30].
With special kernels, GP regression can be faster and less space. Inducing point methods, using inducing points, allow regression to be done in (==2(n + m)) time and in O(=2 + zn) space [21, 22, 23, 13]. We will discuss the details of these inducing point kernels later, but they kernels in their own right, not just approximations to other kernels. Unfortunately, these kernels are low dimensional (having z-dimensional Hilbert space), which limits the expressivity of the GP model.
We present a new kernel, the binary tree kernel, that also allows for GP regression O(n + m) space and O((n + m) log(n m)) time (both model fitting and prediction). The time and space complexity of our method is also linear in the depth of the binary tree, which is naïvely linear the can increase the depth sublinearly. Training some kernel parameters takes time quadratic in the depth of the tree. The dimensionality of the binary tree kernel is exponential in the depth of the tree, making it much more expressive than inducing points kerne.. Whereas for an inducing points kernel, the runtime is quadratic in the dimension of the Hilbert space, for the binary tree kernel, it is only logarithmic- exponential speedup.
dimension of the data, although in practice
simple depiction of our kernel shown in Figure 1, which we will define precisely in Section 3. First, we create procedure for placing all data points the leaves of a binary tree. Given the
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

binary tree, the kernel between two points depends only on the depth of the deepest common ancestor. Because very different tree structures are possible for the data, we can easily form an ensemble of diverse GP regression models. Figure 2 depicts a schematic sample from a binary tree kernel. Note how the posterior mean piecewise flat, but the pieces can be small.
On a standard suite of benchmark regres- sion tasks [25], we show that our kernel usu- ally achieves better negative log likelihood (NLL) than state-of-the
0
sparse methods and conjugate trrgggatt-rrdient---e gradient-based "exact" methods, at lower computational cost in the big-data regime.
0.3
0.8
There are not many limitations to using our ker- nel. The main limitation is that other kernels sometimes capture the relationships in the data better. We do have good procedure for un- derstanding when data has more Matérn charac- ter or more binary tree character (except through running both and comparing training NLL). But given that the binary tree kernel usually outperforms the Matérn, we'll tentatively say the best first guess is that a new dataset has more binary tree character. One concrete limitation for some appli- cations, like Bayesian Optimization, is that the posterior mean is piecewise- flatt so gradient-based heuristics for finding extrema would not work.
1
aa
24
Figure 1: A binary tree kernel with four data points. In this example, (1) 1, k( z2) k(#1, = 0.8, and k(r1, 244) = 0.3.
In contexts where piecewise- flat posterior mean is suitable, we struggle prefer a sparse or sparse variational GP to a binary tree kernel. The most thorough approach would be run both and see which has a better training NLL, but if you had pick one, the binary tree GP seems to be better performing and comparably fast. If minimizing mean-squared error is the objective, the Matérn kernel seems to do slightly better than the binary tree. If the dataset is small, and one needs very fast prediction, Matérn kernel may be the best option. But otherwise, if one cares about well-calibrated predictions, these initial results we present tentatively suggest using a binary tree kernel over the widely-used Matérn kernel.
when one would
The log-linear time and linear space complexity of the binary tree GP, with performance exceeding "normal" kernel, could profoundly expand the viability of GP regression to larger datasets.
2
Preliminaries
Our problem setting is regression. Given a function f -+ R. for some arbitrary set X, we would like to predict f(z) for various
X. What we have are observations of f(x) (zz) for various (other) E X. Let X E X be an n-tuple of elements of X. and let y E K" be an n-tuple of real numbers. such that yı f(,)) + (O, A), for IR X and y comprise our training data.
With an m-tuple of test locations X' E Xmm let y E R", with y/ = f({)). y' is the ground truth for the target locations. Given training data, would like produce distribution over IR for each target location X/, such that it assigns high marginal probability to the unknown y/c Alternatively, we sometimes would like to produce point estimates i/ in order to minimize the squared error (ji v/)².
GP prior over functions is defined by a mean function m - IR, and a kernel le : X' - IR. The expected function value point is defined to be m(x), and the covariance of the function values at two points x and 22 is defined to be ((2),))). Let kxx E Rxx" be the matrix of kernel
Figure 2: A schematic diagram of function sampled from a binary tree kernel. The function is over the interval [0, 1], and points on the interval are placed onto the leaves of depth-4 binary tree according to the first bits of their binary expansion. The sampled function is in black. Purple represents the sample if the tree had depth 3, green depth 2, orange depth 1, and red depth 0.

---

values (Kxx); = k(x,, X2), and let mx E R" be the vector of mean values (mx); = (X,).. For a GP to be well-defined, the kernel must be such that x is positive semidefinite for any X E X For point E X, Let Kx. E I be the vector of kernel values: (Kxx);
a)), and let Ke = xx Let X 2 0 be the variance of observation noise. Let Me and o be the mean and variance of our posterior predictive distribution at I. Then, with K Amy = (Kxx + X1)-1
Ug == (y mx) K X + m(xx
(1)
o : == Ka + A.
(2)
See Williams and Rasmussen [28] for derivation. We compute Equations 1 and 2 for all x E X
3 Binary tree kernel
We now introduce the binary tree kernel. First, we encode our data points as binary strings. So we have X' = 189, where B {00, 1}. and q E N.
If X' = Rd we must map I2 + B". First, we rescale all points (training points and test points) to lie within the box |0, 1]s. (If we have a stream of test points, and one lands outside of the box [0, 1]6, we can either set K, to O for that point or we rescale and retrain in O log n.) time.) Then, for each x E [0, 1] for each dimension, we take the binary expansion up to some precision and for those d x p bits, we permute them using some fixed permutation. We call this permutation the bit order, and it is the same for all x E [0, 1] Note that now q dp. See Figure 3 for an example. We optimize the bit order during training. and we can also form an ensemble of GPs using different bit orders.
For æ E B", let x be the first i bits of x. [expression] evaluates to 1, if expression is true, otherwise O. We now define the kernel:
I1 (0
Definition 1 (Binary Tree Kernel), Given weight vector E IR, with w غ 0 and ||w||1 =1 1,
22
(0.
AS
1001 1001
(T1, Σ Wi
21
22
32 111000
Figure
3: Function from
So the more leading bits shared by xi and Ø2, the larger the covari- ance between the function values. Consider, for example, points 21
[0, 1|2 B
and xz from Figure 1, where 21 is (left, left, right), and x4 is (Ieft, right, right); they share only the first leading "bit". We train the weight vector to maximize the likelihood of the training data.
Proposition 1 (Positive Semidefiniteness) For X E X", for k = Kw Kxx 2 0.
Proof. Let s E U B' be a binary string, and let |s| be the length of Let [1s] E I² with (X(s))s = [x554 = 8 [(a) X is clearly positive semidefinite. Finally, Kxx ΣT Eeer w. wp X and recall wi 2 0, so Kxx 0.
4 Sparse rank one sum representation
In order to do GP regression in (n) space and On log n) time, we develop "Sparse Rank One Sum" representation of linear operators (SROS). This was developed separately from the very similar Hierarchical matrices [1], which we discuss below. In SROS form, linear transformation of a vector can be done in (nn) time instead of On²)) We will store our kernel matrix and inverse kernel matrix in SROS form. The proof of Proposition exemplifies representing a matrix as the sum of sparse rank one matrices. Note that each X1s is sparse-if is large, most Xis) 's are the zero vector.
We now show how to interpret an SROS representation of an n X n matrix. Let [n] = {1, n). For E N. let T. [r]m x [r]m x IR" x R -> " construct linear operator from four vectors. Definition 2 (Linear Operator from Simple SROS Representation). Let p.p' E [r]", and let u, u' e R.. For I E [r), let u- E RR" be the vector where u.-l u, [p, I]] likewise for ul and p/. Then: [p,p'' u, w') + Em wn (u(''''
3

---

We depict Definition 2 in Figure 4. p and p' represent par- titions over elements: all elements with the same integer value in the vector belong to the same partition. Note that r, the number of parts in the partition, need not exceed n, the number of elements being partitioned. If p p' (which is almost always the case for us) and the elements of p. u, and u' were shuffled that all elements in the same partition were next to each other, then L(p, p': u, u') would be block diag- onal. Note that L(p, p', u. ut') is not necessarily low rank. If p is the finest possible partition, and p = pp' L(p. p', u, u') is diagonal. SROS matrices can be thought of as generalization of two types of matrix that are famously amenable to fast com- putation: rank matrices (all points in the same partition) and diagonal matrices (each point in its own partition).
p'
0
1
1
U
1.5
1.5-0.5
-1
0.5
p u
u
0
0
1.5
1.5
-0.5
-0.5
0
-1
-1
0
0.5
0.5
Figure 4: A matrix standard form constructed from matrix in SROS form. The large square depicts the matrix Lpp.p' u, u') E R5x5 with elements colored by value. See u== for a color legend.
We now extend the definition of L to allow for multiple p, p', u, and u' vectors.
Definition 3 (Linear Operator from SROS Representation). Let L: [r]m*a "p]nxg x x R - Let P,
[r]mxg, and let U, U Rnxg Let P, Un etc. be the ith columns of the respective arrays. Then: L(P, (U,U)) E L(P.,, P. (.
Algorithm performs linear transformation of vector using SROS representation (ng)) time.
Algorithm 1 Linear Transformation with SROS Linear Operator. This can be vectorized on Graphical Processing Unit (GPU), using c.g. torch Tensor index_nade add__ for Line and non-slice indexing for Line 6 [19]. Slight restructuring allows vectorization over [q] well.
Require: P. E [r] U.U" E ER* x ER"
Ensure: y == L(P, P', L'U,PP,UUUU)
I: t 0 R"
2: for E ql do
O(ng) time
3: P.P', uu <- P.j P! U.i U
4:
z 2 E R
> 22 will store the dot product ((u')''=I)T T
5: for j E [n] do 2pl, 2//, + u/
((n) time
6: for E [n] do yj Vj 2p;
(n) time
return
We now discuss how to approximately invert a certain kind of symmetric SROS matrix, but our methods could be extended to asymmetric matrices. First, define a partial ordering over partitions. For two partitions p', we say p' < if p' is finer than or equal to p. that is, P/J =py = pj Py, Using that partial ordering, symmetric SROS matrix can be approximately inverted efficiently if for all S i, i' < P., < P..y or P., < P...- As the reader may have recognized, our kernel matrix Kxx can be written as an SROS matrix with this property.
We will write symmetric SROS matrices in a slightly more convenient form. All (w')== must be a constant times up=(, We will store these constants in an array C. Let L(P, C, U) shorthand for L(P, P. U, C @ U), where @ denotes element-wise multiplication. For L(P, C, C,U) be symmetric, it must be the case that Pf:- = Pys = Cj Cys.. Then, all elements of U corresponding to a given wi=l are multiplied by the same constant. We now present an algorithm for calculating (L(P,C,U) A1)-1. for À * 0, which is approximate inversion of L(P. C,). We have not yet analyzed numerical sensitivity for > 0, but we conjecture that all floating point numbers involved need be stored to at least 1og,(1// (L(P.C,U) + X x (L(P, x C. A.PP,,-CCUU) +)-
(1/X) bits. Without loss of generality, let X = 1, and note
By assumption, all columns of P are comparable with respect to the partial ordering above, can reorder the columns of P such that P.; P.. for j. The key identity that we use develop our fast inversion algorithm is the Sherman-Morrison Formula:

---

(A + cou A A
un T A-1
c- A
Starting with A = 1, we add the sparse rank one matrices iteratively, from the finest partition to the coarsest one, updating 1 as we go. We represent (L(P, C, U+++++++ + +1) in the form I+ L(P. C², U*). so we write an algorithm that returns C" and U". We can also quickly calculate log [ C,U)+ at the same time, using the matrix determinant lemma: |A + cur" = (1 + cu c T AA)
(is coorser thon
Theorem 1 (Fast Inversion). For P E [r]maq and C,U E Raxq Pf. 2
for < j.
then there exists C², U' E R=*0, such that (L(P, C, + - L(P.C',U) There exists an algorithm for computing C' and U' that takes (nqq²) time.
Proof. For X E R let X. E E R=x((-1) be columns + through q of matrix X (inclusive).
Let A, = I 1(Pp, C. Ug and A, I + I( 1, U (1 for some C' and U For the base case of = this holds trivially. We show it also holds for i- and we can compute Ch,saq' Ui in O(n(g - i)) time. Let p P.s- u = U.. and c= C..-- Consider u -- where each element is zero unless the corresponding element of p equals I. What do we know about the product A, u== (as seen in Equation 3y?
Now suppose A, can be written as
Because the columns of P from coarser partitions to finer ones, all of the vectors generating the sparse rank one components of T.(p,,+1:q: Ch U ('(+1:;g) are from partitions that are equal to or finer than p. Thus, they are either zero everywhere w² is zero, or zero everywhere wz is Thus, when multiplying
nonzero. Vectors of the second kind can be ignored, as cvv u²² L(P.,+1:q. 1lgg U (1g)) by u²1 the only relevant vectors are filled with zeros except where the corresponding element of equals 1. So we can get rid of those rows of P.,++1:p. Ca and U 11 Suppose there are nı elements of p that equal 1. Then L(Pr, C1 1111 U 1441:1)) ) involves nı rows, and can be computed in (n - -)) time. Moreover, this product, which we'll call (u')== is only nonzero when the corresponding element of p equals 1, so it has the same sparsity pattern as wp=i_ The other component of A, is the identity matrix, and Tu clearly has the same sparsity as u²=!. Thus, returning to Equation 3, when we add w ) to A,, we update A, with an outer product of vectors whose sparsity pattern is the same as that of u==!
For each 4 A; need not be updated with each wp=l one at a time. For I + I', w== and w==l" are nonzero at separate indices, so wn and (u')==l' are nonzero at separate indices, so the extra component of A, that appears after the w²=!" update is irrelevant to the wn= update, because (w==l)T(
O. Since the w²=l update takes 0(n(qq -)) time, all of them together take O mg - i)) time, which equals (n - i)) time. Calculating an element of c' only involves computing the denominator in Equation 3, using a matrix-vector product already computed. So we can write Cr,bag and UAAig by adding a preceding column to C2,++.g and UK, +iggg using the same partition p, and it takes O(n(q i)) time.
Following the induction down to i 0, we have time of (nq²).
= I U'), and total
Algorithm 2 also performs approximate inversion, which we prove in Appendix A. It differs slightly from the algorithm in the proof, but can take full advantage of a GPU speedup. In the setting where all columns of U are identical, observe that in Lines 10 and 11, the same computation is repeated for all l: E [i]. Indeed, in this setting, this block of code can be modified to run in O(n) time rather than Onni), making the whole algorithm in O(nq) time, as shown in Proposition 4 in Appendix A.
matrix which is either represented as a low-rank matrix or as 2 x block matrix of Hierarchical matrices [1]. In our SROS format, many of the sparse rank one matrices overlap, whereas in a Hierarchical matrix, the low-rank matrices do not overlap, and converting an SROS matrix into a Hierarchical matrix would typically be inefficient. Hierarchical matrices admit approximate inversion in O(na² log" n)) time, where a is the maximum rank of the component submatrices [11]. However, this is not an approximation a technical sense, as there is no error bound. At many successive steps in the algorithm, a rank 2a matrix is approximated by rank
Hierarchical matrix
5

---

Algorithm Inverse and determinant of /+ SROS Linear Operator. Lines 5 through 11 can all be easily vectorized on GPU. Lines and 10 require torch. Tensor index_add, or equivalent, and lines 6 and 11 require non-slice indexing, which are not quite as fast as some GPU operations.
Require: P [r]nxg C.U E Px4 = (I (I+ L(P,C,U)) x= log T+ L(P,C,U))
Ensure: L(P,
1: z. C",U" <- 0,0 R
2: for E (q.9 1,.., (.....): 1)do do
> (nq")) time
3:
P.C, u, u +- P. C..., U.ú U
4:
2 2 R
2 will store ((u')=i)) umm where c
Ck if pk: = I
5:
for j E [n] do 2p, + 2p, + cyuj
(n) time
6:
for j E [m] do j + -2// (1 + 2p,)
b (n) time
7:
for / [r] do x t x + log(1 + 21)
> O(n) time because r S n
8:
i > 0 then
9:
y f 0 Rixi
> ( time
10:
for j. k E [n] x [i - 1] do Upyk + Upyk + U/UUu
D O(ni) time
11:
for j. E [n] x [ (6-11 do Ujx - Un + CjulUops
D. O time
return C² U
matrix [10]; to our knowledge there no analysis of how resulting errors might cascade. After converting an SROS matrix to hierarchical form, this rough inversion would take O(ng² log" n) time.
5 Binary tree Gaussian process
We now show that our kernel matrix Kxx can be written in SROS form, with P containing succes- sively finer partitions. Thus, Kxx x can be approximately inverted quickly, for use in Equations and 2. Next, wl'll show that we can efficiently optimize the log likelihood of the training data by tuning the weight vector along with the bit order. The log likelihood can be calculated in O(ng logn) time and then the gradient w.r.t. w in O(nq²) time.
Recall from the proof of Proposition 1: Kxx = EiER w;X
XT XT 1o/= where X(a) E TR
with (Xi))
X
+1 s]]. So we will set P.,.. C... and U.i that L(P.,,, C, U..)
EVER' W;X Let p.. partition the of points X so that points are in the same partition if the first i bits match. Now, requiring the first + 1 bits to match is stricter criterion than requiring the first i bits to match, so the P.,, grow successively finer. For any piece of the partition where the first i bits of the constituent points equals the bitstring s, the corresponding sparse rank one component of Kxxx is w. X [e) XT- So let U. = 1", and let C. = w;1".
[s]-
Proposition 2 (SROS Form Kernel). Kxx L(P, C, C), defined above.
This follows immediately from the definitions. To compute these partitions P.s,, we sort X, which set of bit strings. And then we can easily compute which points have the same first i takes (ng log n.) time. Now note that U. = U. for all i, i', so (Kxx + A1)-1 and [xxx + A can be computed in ((ng) time, rather than (nq²)).
bits. This all
The training negative log likelihood of GP is that of the corresponding multivariate Gaussian on the training data. So: NLL(ww
(y T (Kxx(w) + X + log KK (xx(w)) A7| + n log(2r) )) This can be computed in O(nq) time, since matrix-vector multiplication takes O(ng) time for matrix in SROS form. So if the bit order is unchanged, an optimization step can be done in O(nq) time, and if the data needs to be resorted, then in Onn log n.) time. On the largest dataset we tested (House Electric), with n R 1.3 million and = 88, sorting the data and computing takes about 0.96 seconds on GPU, and then calculating the negative log likelihood takes about another 1.08 seconds. We show in Appendix B how to compute V., NLL in (ng²²) time.
To optimize the bit order and weight vector at the same time, we represent both with single parameter vector 0 E TRy, with |00| 1. To get the bit order from 0, we start with a default bit order and permute the bit order according to permutation that would sort 0 in descending order. To get the

---

Algorithm 3 GP Regression with a binary tree kernel.
Require: X E _=xq y R" X' e B"Xg R À E R+
Ensure: μxy and xx are the predictive means and variances at X', and nll the training negative log likelihood.
I: X + X X"
2: X perm t- Sort(X)
The rows of X are sorted lexically from leading bit to trailing bit. (( log(n + m)) time.
3: for j,ú [n m] x [a] do Pj #of unique rows in X 1:j.1si
4:
Mujia is the first rows and i columns of M. O =m)g) time.
5: P perm-1(Pt)
This "unsorts" the input. O - -m]g) time.
6: P.P' + Pim Patlntm
7: U,U', - - 1mxa, 1xxg (
8: C,C 10C+1"mm 1"
9: C. CT. loggeee - Invert PC,-)
Uses Algorithm Speedup O(nq) time because columns of U are identical.
10: C-1 logdet + - logdet, + (Ag
11: <- LinTransform P.U C ② U
. + +
D Uses Algorithm 1 to compute the Woodbury vector, O(nq) time.
12: px +- --inirannformm (P', U" C U.,))
D O( m)g) time.
13: nll + (y + logdet + n log (2
14: Cuose ÜTpose <- Invert( P, C- 1,,,,-100 1-)
D O(1 + m))g time.
15: Chre Upree +- Exprec n+1m+mm U Uppec 1+1 16: C U Invert( Pr, Cree U
Ommq²) time; extra factor of q because columns Upree are not identical. O(mg) time.
17: ox f X(1 + SSmmachhow((( (Co Jan C U+r)
18: return px', nll
weight vector, we sort 0 in descending order, add a 0 at the end, and compute the differences between adjacent elements. When there are ties in the elements of the choice bit order does not affect the negative log likelihood (or the kernel at all) because the relevant associated weight is O. The negative log likelihood is continuous with respect to 0, and when all values of are unique, differentiable e* /|
with respect to 0. Letting
minimize loss w.r.t. 6 using BFGS [8].
To calculate the predictive mean at a list of predictive locations X', we first multiply y by (Kxx- + A and then we multiply that vector by xxx'. We obtain both Xxx and Xxxx in SROS form as follows. Let X X o X' be the concatenation of the two tplles, now an (n + m))uuple. Writing K x L(P, Ć, U), the arrays on the r.h.s. can be computed in O + log(n + m)) time. Then, with P. C, and U being the first n rows of P. Č2 Ü, Kxx L(P, U). And letting p" and U" be the last m rows, Kxxx L(P, P" C C@U,U" U"). Thus, the predictive mean Hs from Equation 1 can be computed m locations in O((n + m)q log(n + m)) time.
The predictive covariance matrix, which extends the predictive variance from Equation 2, is calculated Σx x + XIn Kx X - the Schur complement. From a property of block matrix inversion, the last m columns of the last m rows of (K x + AT) equals ((K x + Amm+n)/ ( xxx)=1x___ So get the predictive precision matrix in O((n + m)g log(nu m)) time by inverting K x + Al and taking the bottom right m x m block. Then, we get the predictive covariance matrix by inverting that. This takes ((mq²) time, since it does not have the property of all the columns of U being equal. If only want the diagonal elements of an SROS matrix (the independent predictive variances in this case), we can simply sum the rows of Co U @Uin in ((m)) time. Thus, in total, computing the independent predictive variances requires O(n + m)g log(1z n + m) + mq²) time. See Algorithm 3.
x + Al Kxx where denotes
6 Related Work
All existing kernels of which are aware for linear time GP regression unstructured data involve
7

---

inducing points (related to the Nyström approximation [27]) or inducing frequencies. For given set of inducing points Z, for some base kernel k, the inducing point kernel (in its most basic form) the following, although subtle variants exist: k² (I, Kzz K Zzz [21].
Sparse Gaussian Process Regression (SGPR) involves selecting Z, and then using k2 (or variant). Notably, K xx = XxzKzz is low rank, providing computational efficiency. The predictive mean and covariance have compact form, with observational noise A: Pz(Z) = KaZ(AZzz + kzxkxz) > Xzxy and o2 (Z) = Kzz KKzz + 1-1xxx - Kzz'.
Titsias's s [24] sparse variational kernel is also low rank and uses inducing points. The sparse variational
GP (SVGP) is constructed as the solution to a variational inference problem. It depends on inducing points Z. data points X, and observed function values y. We have focused on Gaussian processes with 0 mean, but the SVGP method uses a nonzero prior mean along with a kernel: m mvsu ) = Mu (Z) and J-SGPP (2r x') = k(x,,²)) x2 K Zz" + zz (Z). Given the dependence on X and y. this is not true probability distribution over function space. The variational problem underlying this kernel also provides guidance in how to select the inducing points Z. For further discussion of the kernel underlying the SVGP method, see Wild et al. [26].
An inducing point kernel with 2 inducing points produces z 2-dimensional reproducing kernel Hilbert space (RKHS). The dimensionality of the RKHS relates to the expressivity of the kernel. Whereas an inducing point method buys a z-dimensional RKHS for the price of O(22n) n) time and O(2n) space, the binary tree kernel produces 29 dimdimensional RKHS in (gq) time and space- an exponential improvement. (Observe that we can find 2 linearly independent functions of the form k(,, ((-,a))-nne one for each of the 29 leaves might belong to.) Wilson and Nickisch [29] develop method for speeding up inducing point methods significantly, especially in low-dimensional settings.
given set of m inducing vectors 8t. k(x,, x²) = 1/m >n c2ss (x - z')). Dutordoir et al. [6] propose an inducing frequencies kernel for low dimensional data, in which k(2,, x²) is a special function of x x'.
Lázaro-Gredill, al. [16] propose an inducing frequencies kernel:
On one- one-dimensional data, filtering/smoothnng methods perform Bayesian inference over functions in (nn) time [14, 12]. few non-O(n) methods bear mentioning. We are not the first to consider a kernel over points on the leaves of tree [18, 17] on the leaves of multiple trees [7]. but their methods take (nn²) time. On certain kinds of structured data, Toeplitz solvers achieve ((n)) time complexity [30]. Cutajar et al.'s 's [3], Wang et al.'s 's [25], and others' use of a conjugate gradients solver to replace inversion/factooizzaionn has unclear time complexity between On²) and O(n²), depending on the kernel matrix spectrum.
7 Experiments
In Table 1, we compare our binary tree kernel and a binary tree ensemble (see Appendix C for details on the ensemble) against three baseline methods: exact GP regression using a Matérn kernel, sparse Gaussian process regression (SGPR) [23], and stochastic variational Gaussian process (SVGP) [13]. We evaluate our method on the same open-access UCI datasets [4] as Wang et al. [25], using their same training, validation, and test partitions, and we compare against the baseline results they report. For the binary tree (BT) kernels, we use p = min(S, [150/d] d] + 1), and recall q = pd. We set x 1/n. We train the bit order and weights to minimize training NLL. For the binary tree ensemble (BTE), we use 20 kernels. For the Matérn kernel, we use Blackbox Matrix-Matrix multiplication (BBMM) [9], which uses the conjugate gradients method to calculate matrix-vector products with (Kxx + A7)-1. SGPR uses 512 data points and SVGP uses 1,024 inducing points. We report the mean and two standard errors across 3 replications with different dataset splits. For further experimental details, see Appendix C. BTE achieves the best NLL on 6/12 datasets, and best RMSE on 5/12 datasets (including some ties). Out of the 4 largest datasets, BT/BTE is fastest on 3. The run
'The conjugate gradients method takes On (()) time, where is the condition number of the kernel matrix. Poggio al. [20] say "claims about the condition number of random matrix A should also apply to kernel matrices with random data." If they mean a Wishart random matrix (which should be if, e.g. k(x,, x²) = x x'), 2²), that would be the square of the condition number of the corresponding Gaussian random matrix, which grows as (nn) [2]. Putting it all together, we get O(n") for conjugate gradients. We don't 't know how quickly preconditioning can reduce the condition number.
8

---

10
and Matém (BBMM); m = 1.57#0.16
SVGP; m 1.1:0.04
SGPR; m 0.7520.07
BT/BTE; m 0.58
10a
Run time (m)
10a
104
10
10c
data points
Figure 5: Run times given dataset size. For BT, the trendline is calculated controlling for log(4) with affine regression, and then setting g 150. The slope w.r.t. log(g) is 2.82 + 1.06. Theoretically, all slopes are too low except for that of SVGP, presumably because of overhead in the small-data regime.
times are plotted in Figure 5. The code is available at
BT performs noticeably worse than BTE for test NLL on CTSlice due over-fitting. There enough degrees of freedom when optimizing the bit order (d = 385) that BT kernel can over-fit to the training data. The ensemble over multiple bit orders is much more robust.
DATASET
d
BTE
BT
MATERN (BBMM)
SGPR
SVGP
POLITELE
9,600
26
-0.625 = 0.035
-0.400 : 0.040
-0.180 0.056
-0.094 + 0.008
-D.001 0.008
ELEVAIURS
10,623
18
0.649 0.032
0.646 + 0.023
0.619 0.005
0.130 0.060
0.609
11,122
-0.708 + 0.433
-0.806
0.11920.044
0.291
0,272 0.0019
KIN40x
25,600
0.004
0.008
-0.558
0.236
PROTIIN
29,267
0.781 + 0.023
0.845 + 0.026
1.015+0.0500
0,970 0.010
1,035 + 0.006
KeGGDIR
31,248
20
-1.031 0.020
-1.029
0.021
-0.199 0.381
-1.123 0.016
-0.940 + 0.020
CTstarm:
-..073
0.005
KEGGU
27
D.667
0.007
0.667 0.007
-0.419
0.027
0.984
0.66G 0.007
3DROAD
238,319
3
-0.251
= 0.009
0.252
0.006
0.909 +
0.943 0.002
0.697 0 002
SONG
329,820
90
1,330
0.003
1,331 # 0.003
1.206 0.024
1,213 0.003
1,417 000
373,280
0.003
0.267
0.106
0,224 0,050
HouseRune
1,311,530
2.569 = 0.006
2.402
0.012
1.....12 0.001
11.100 0.039
POLITTLE
+ 0.006
0.004
+ 0.012
BLEVATURS
10,623
18
0.478
0.021
D.476 + 9.018
0.394 0.006
0.437 0.011
0.899 0.009
BIKE
11.122
17
0.118 + 0.057
0.103
0.029
0.220 1 0.002
0,362
0,303 1 0.004
KIN40s
0.006
0.099
0.623
0.656 0.010
KROODUR
31,248
20
0.086 1 0.003
0.086
0.003
0.086 = 0.005
0.104 0.003
0.096 : 0.001
CTsLaEE
34,240
385
0.116 + 0.009
0.132 + 0.009
0,262 + 0.448
0.011
1,003 + 005
KEGGU
3DROAD
278,319
0.002
0.186 0.001
0.101
SDNG
329,820
90
0.914 0.003
0.916 - 0.003
0.807 = 0.024
0.803 0.002
0.99 + 0.000
BUZZ
373,280
77
0.801 0.002
0.801 4 0.0U2
0.288 0.018
0.300 0.004
0.304 0.012
HOUsEELEC
1,311,530
0.001
0.001
0.055
POLITELE:
26
5.16
1.162
1.15 0.068
ELEVATORS
1.16
0.092
BIKE
11.122
17
2.68 0.15
0.69 0.015
1.33
1.28 0.093
KiNs0x
25,600
1.44 0.028
0.71 11.45
1.62
3.26 0.23
2,27
3.31
CTSLaUE:
34,240
385
62.01 + 0.92
3.32 5.0
2.16
3.87 0.334
KEGGU
40,708
27
7.45 0.54
6,32" 1..41
2.22
4.78 0.4
3DROAD
278.319
1.83 0.12
126.37 20.92*
12.01:5.51
34.09 3.119
90
31.87
33,79
7.89
Buzz
373,250
77
20.18
+ 6.86
571.13 66.34"
29.25
18.33
46.35 2.99
HOUsHELOC
1,311,539
11
118.41 + 3.93
575.64 6.94
367.71 4.7
Table 1: NLL (top), RMSE (middle). and run time in minutes (bottom) regression datasets, using a single GPU (Tesla V100-SXM2- 16GB for BT and BTE and Tesla V100-SXMM-32GB for the other methods). The asterisk indicates an estimate of the time from the reported training time on 8 GPUS, assuming linear speedup number of GPUs and independent noise training times per GPU. All columns except BT and BTE come from Wang et al. [25].

---

8
Discussion
We have proven that the binary tree kernel GP is scalable. Our empirical results suggest that it often outpredicts not just other scalable methods, but even the popular Matérn GP. If the results in this paper replicate in other domains, it could obviate wide usage of classic GP kernels like the Matérn kernel, well as inducing point kernels. Sometimes, our kernel fails to capture patterns in the data; some functions' values simply do not covary this way. But other kernels we tested seemed to fail like that even more.
Our contributions to linear algebra and kernel design may significantly increase the size of data sets on which GPs can do state-of-the-ar modelling.
References
[1] Mario Bebendorf. Hierarchical matrices. Springer, 2008.
[2] Zizhong Chen and Jack Dongarra. Condition numbers Gaussian random matrices. SIAM Journal on Matrix Analysis and Applications, 27(3):603-620, 2005.
[3] Kurt Cutajar, Michael Osborne, John Cunningham, and Maurizio Filippone. Preconditioning kernel matrices. In International conference on machine learning, pages 2529-2538. PMLR,, 2016.
[4] Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. URL http:////rr ics. uci odu/ml.
[5] Philippe Duchon, Philippe Flajolet, Guy Louchard, and Gilles Schaeffer. Boltzmann samplers for the random generation of combinatorial structures. Comb. Probab. Comput., 13(4-5): 577-625, jul 2004. ISSN 0963-5483. doi:
[6] Vincent Dutordoir, Nicolas Durrande, and James Hensman. Sparse Gaussian processes with spherical harmonic features. In International Conference on Machine Learning, pages 2793- 2802. PMLR, 2020.
[7] Dai Feng and Richard Baumgartner. Random forest (RF) kernel for regression, classification and survival. CoRR, 0bs/20009000033 2020. URL https: //arxiv...vv./// org/ htg//bs/20000 00089.
[8] Roger Fletcher. Practical methods of optimization. John Wiley & Sons, New York, 2nd edition, 2013.
[9] Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. Gpytoreh: Blackbox matrix matriix Gaussian process inference with gpu acceleration, Advances in neurul information processing systems, 31, 2018.
[10] Wolfgang Hackbusch. A sparse matrix arithmetic based on H-matrices. Part I: Introduction to H-matrices. Computing, 62 1999.
[11] Wolfgang Hackbusch, Boris N Khoromskij. and Ronald Kriemann. Hierarchical matrices based on a weak admissibility criterion. Computing, 73(3)2207-243, 2004.
[12] Jouni Hartikainen and Simo Särkkä. Kalman filtering and smoothing solutions to temporal Gaussian process regression models. In 2010 JEEE international workshop on machine learning for signal processing, pages 379-384. IEEE, 2010.
[13] James Hensman, Nicolò Fusi, and Neil D. Lawrence. Gaussian processes for big data. In Proceedings of the Tweth Conference Uncertainty in Artificial Intelligence, UAI' 13, page 282-290, Arlington, Virginia, USA, 2013. AUAI Press.
[14] RE Kalman. new approach linear filtering and prediction problems. Journal of Basic Engineering, 83(1):95-108, 1960.
10

---

[15] Sanyam Kapoor, Mare Finzi, Ke Alexander Wang, and Andrew Gordon Gordon Wilson. Skiing on simplices: Kernel interpolation on the permutohedral lattice for scalable gaussian processes. In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 5279- 5289. PMLR, 18 24 Jul 2021. URL htl.https:///e https::/pproceeddngggeee/// mlr press/ https///19//ppppoo/ html.
[16] Miguel Lázaro-Gredie Greaill,,, Joaquin Quinonero-Caddlla, Carl Edward Rasmussen, and Aníbal R Figueiras-Vidde Sparse spectrum Gaussian process regression. The Journal of Machine Learning Research, 11:1865-188 1881, 2010.
[17] Julien-Charles Lévesque, Audrey Durand, Christian Gagné, and Robert Sabourin. Bayesian optimization for conditional hyperparameter spaces. In 2017 International Joint Conference on Neural Networks (IJCNN), pages 286-293, 2017.
[18] Xingchen Ma and Matthew Blaschko. Additive tree-structured covariance function for con- ditional parameter spaces in Bayesian optimization. In Silvia Chiappa and Roberto Calandra, editors, Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pages 1015-1025. PMLR, 26-28 Aug 2020.
[19] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang. Zachary VeVito. Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high- performance deep learning library. In Advances in Neural Information Processing Systems 32, pages 8024-803 2034-0035.00 Curran Associates, Inc., 2019. URL http://paperrs neurips cc/paper/ 9015-pry 1155--yytrch--- an- imperatiie- style- www.hhhihhhhhheeeeeee pwwformance- -eabeeineerringg -iinggry.bbr.... pdf.
[20] Tomaso Poggio, Gil Kur, and Andrzej Banburski. Double descent in the condition number. arXiv preprint arXiv: 19/2.06/90, 2019.
[21] Joaquin Quinonero-Candela and Carl Edward Rasmussen. A unifying view of sparse approxi- mate Gaussian process regression. The Journal of Machine Learning Research, 6:199--1959 2005.
[22] Edward Snelson and Zoubin Ghahramani. Sparse Gaussian processes using pseudo-inpure In Weiss, B. Schölkopf, and J. Platt, editors, Advances in Neural Information Processing Systems, volume 18. MIT Press, 2005.
[23] Michalis Titsias. Variational learning of inducing variables in sparse Gaussian processes. In David van Dyk and Max Welling, editors, Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics,, volume 5 of Proceedings of Muchine Learning Research, pages 567-574, Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA, 16-18 18 Apr 2009. PMLR. URL https::/prooeedding mlr
html.
[24] Michalis Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Artificial intelligence and statisties, pages 567-574. PMIR, 2009.
(25) Ke Wang, Geoff Pleiss, Jacob Gardner, Stephen Tyree, Kilian Q Weinberger, and Andrew Gor- don Wilson. Exact Gaussian processes million data points. Advances in Neural Information Processing Systems, 32, 2019.
[26] Veit Wild, Motonobu Kanagawa, and Dino Sejdinovic. Connections and equivalences be- tween the Nyström method and sparse variational Gaussian processes. arXiv preprint arXiv:2106.01121. 2021.
[27] Christopher Williams and Matthias Seeger. Using the Nyström method to speed up kernel machines. In T. Leen, T. Dietterich, and V. Tresp, editors, Advances in Neural Information Processing Systems, volume 13. MIT Press, 2000. URL https: wttp:///dn/oo.o..../e...oo neurips. ottpae/rrooooziii//eeeeaaaaaaaae
Paper pdf.

---

[28] Christopher K Williams and Carl Edward Rasmussen. Gaussian processes for machine learning. volume 2. MIT press Cambridge, MA, 2006.
[29] Andrew Wilson and Hannes Nickisch. Kernel interpolation for scalable structured gaussian processes (kiss-gp). In International conference on machine learning, pages 1775-1784. PMLR, 2015.
[30] Yunong Zhang, William E Leithead, and Douglas Leith. Time-series Gaussian process regression based on Toeplitz computation of O operations and Onn)-level storage. In Proceedings of the 44th IEEE Conference on Decision and Control, pages 3711-3716. IEEE, 2005.
12


Reviews:
[{"review_id": "wf5VjLeG5Vp", "reviewer": "Reviewer_GqxJ", "summary": "The paper develops a novel binary tree kernel function for fast Gaussian process computation. The kernel is based on the similarity of permuted bit strings of the two input features, which induces a kernel matrix with a sparse rank one sum representation. The computational gain is achieved by efficient algorithms for matrix inversion and matrix-vector multiplication. ", "questions": "Major comments:\n\n1. As far as I understand, for any given bit permutation, $||x_1^{raw} - x_2^{raw}||_2 \\rightarrow 0$ does not imply $k_w(x_1, x_2) \\rightarrow ||w||_1$, where $x_1^{raw}, x_2^{raw} \\in \\mathbb{R}^d$ are the raw features, and $x_1, x_2 \\in \\mathbb{R}^q$ are their bit representation. However, in some low-dimensional problems such as spatial applications, it makes sense that data with two closer features (in terms of Euclidean distance) have stronger correlation, which cannot be captured by the proposed kernel. I am curious how the proposed method would perform in such scenarios. \n\n2. It appears to me that after the bit permutation is updated during model training, we need to re-sort the bit feature matrix $X$ to obtain a new $P$. Does this mean the actual training time complexity is $O(n\\log{n})$ instead of $O(n)$?\n\n3. The proposed kernel seems to be highly nonstationary in nature, which could be one of the reasons why it outperforms the standard stationary Mat\\'ern model in real applications where the stationarity assumption does not make sense. It would be interesting to see how it compares with other nonstationary GP models such as Zhang \\& Williamson (2019) and Trapp et al. (2020).\n\n4. I think the high-dimensional kernel parameter $\\theta$ could be somewhat hard to learn. How sensitive is the performance to the initial values of $\\theta$? Also, the performance of BT is based on multiple initializations, which may not be fair to other competing methods.\n\n\nMinor comments:\n\n* Could you please elaborate more on why the matrix inversion algorithm is called \"approximate\" inversion?\n* Figure 2: Are those piecewise constant functions the posterior means of binary tree kernel GPs? I think a realization from the proposed GP should be piecewise smooth rather than piecewise constant.\n\n\nReferences:\n\nTrapp, M., Peharz, R., Pernkopf, F., \\& Rasmussen, C. E. (2020). Deep structured mixtures of gaussian processes. In International Conference on Artificial Intelligence and Statistics (pp. 2251-2261). PMLR.\n\nZhang, M. M., \\& Williamson, S. A. (2019). Embarrassingly parallel inference for Gaussian processes. Journal of Machine Learning Research.", "limitations": "Please see the major comments above.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n* Instead of approximating a kernel matrix (or its inverse), the paper presents a novel idea of constructing a kernel with a computationally favorable structure.\n* The proposed sparse rank one sum representation and algorithms could be applicable to other settings.\n* The method is nicely illustrated.\n\nWeakness (see major comments below for details):\n* Unclear performance in low-dimensional settings.\n* Training time complexity appears to be $O(n\\log{n})$ instead of $O(n)$.\n* The paper could benefit from additional experiment results.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "IpDRLpsrslO", "reviewer": "Reviewer_Lnbh", "summary": "The manuscript proposes a GP covariance function to model multivariate piecewise constant functions so that computationally efficient inference and prediction becomes possible. The authors propose numerical primitives to compute matrix-vector multiplications, the inverse and log-determinant of the system matrix that scale linearly in the number of datapoints, the input dimension of the data and the inverse stepsize of the grid underlying the piecewise constant function both in terms of runtime and space. ", "questions": "* The piecewise constant nature of the modeled function should be emphasized more clearly.\n* The space/time footprint of the algorithms depends linearly (derivative of hyperparameters quadratically) on the data dimension and the stepsize of the underlying grid. This information should be presented more prominently.\n* Also, the fact that the covariance is degenerate (low-dimensional) could be made more accessible.\n* Given that in line 29, the manuscript puts a lot of emphasis on the variance estimate of inducing point methods, I would have expected a detailed comparison of the predictive variance of the proposed (degenerate) covariance with an exact GP with non-degenerate covariance and an inducing point method approximation.\n* Can the numerical algorithms be cast such that also non-Gaussian likelihoods can be treated via expectation propagation, variational Bayes or Laplace approximation?\n* An analysis, why \"better [..] likelihoods\" were achieved could be instructive. In particular, how does the likelihood depend on q? How are hyperparameters optimised\n* There is no comparison to the exact model in terms of modeling accuracy.", "limitations": "Yes", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "1) Clarity\nThe paper is mostly well written and the technical content is rather well accessible. The Algorithms 1+2 and would become simpler if consistent vector algebra was used instead of loops over indices.\n\n2) Originality\nThe proposed covariance function and the derived algorithmical approaches seem to not have been explored before.\n\n3) Significance\nThe method is (currently) specific to GP regression and could extend the range of GP models to high-dimensional data whenever piecewise constant functions and rank-deficient covariances are appropriate.\n\n4) Reproducibility\nAs the code to run the experiments is included in the submission and the datasets are public, it should be rather simple to reproduce the results.\n\n5) Empirical analysis\nThe manuscript contains a large number of numerical results gathered in a big table regarding marginal likelihood, prediction error and runtime. Some more results on more subtle aspects of the covariance function itsel would have been much appreciated.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "DR1qxlWQe7u", "reviewer": "Reviewer_9UL1", "summary": "The paper discusses a new type of kernel for Gaussian process regression based on binary-tree partitions of the input space. The structure of the kernel allows for a fast (linear time) inversion of the kernel matrix which is benefitable for large datasets. The Gaussian process with this kernel gives piecewise-constant predictions.", "questions": "Since the authors claim the superioirity of their kernel over other methods of scaling up GP regression, it is interesting how the experiments were carried out.\nIn specific, it is known that a Matern kernel performance depends on the choice of the smoothness parameter. How did you choose this parameter in your experiments?\n\nIt is also interesting how did you choose the $\\lambda$ parameter of your kernel (which is basically noise level)? Is it possible to train this parameter along with the order and weights?\n\nAlso an interesting thing is how does the choice of the precision $p$ affect the performance? In particular, since the train time scales quadratically with the precision. How did you come up with the precision in the experiments and is it possible to automatically train the precision? Does it makes sense to ensemble over not only different bit orders but also over different precisions?\n\nIt would be nice to see some comparision ageints other partition-based methods like Random Forest Kernel.\n", "limitations": "Authors have clearly stated the limitations of their work. The predictions provided with a GP with binary-tree kernel are piecewise-constant. That makes the kernel unusable in Bayesian Optimization. It would be nice to discuss limitations of the choise of the precision level $p$ and noise level $\\lambda$. ", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 1, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "*Originality*. The work is original. It is somewhat similar in nature to works related to Random Forest kernels but I am not aware of other work based on bit representations of the input data.\n\n*Clarity*. One area where the paper could use improvement is the clarity, especially notation-wise. It is quite hard to read the heavy-weight superscripts throughout the proof of the theorem 1 in particular. It seems like there is a typo on the line 178: missing \\prime? There is also some notation inconsistencies like using the same letters $i$, $j$, $k$ with a different meaning in similar contexts. Also it seems like the proof of the Theorem 1 in the main paper is in between a high-level sketch and a precise mathematical formulation which makes it harder to read. Could be benefitial if the proof outlined a high-level sketch with the details hidden in the supplement.\n\n*Quality*. The quality of research is high. The authors provided proofs of the mathematical statements they make. The authors have also made empirical performance comparison against popular GP kernels, such as Matern or SVGP method. The authors have provided a review of a related and competing methods of scaling up GP regression. An interesting thing would be to compare their method to other partition-based methods, such as a Random Forest Kernel.\n\nThe authors claim that their kernel widely outperforms other methods such as exact regression with a Matern kernel or a sparse GP regression. However, as presented in the Table 1, this is perhaps an overstatement. In particular, Matern kernel beats the BT(E) for moderately-sized datasets both in test RMSE and computational speed, and in terms of the test NLL, it is not a clear cut that a BT kernel is ultimately a better choice.\n\nIt would be also interesting to see some ablation studies regarding the choice of the bit precision and the noise variance $\\lambda$. \n\n*Significance*. The work provides a new type of kernel that could be used in setttings where a piecewise-constant predictions are suitable. The authors provide an algorithm for a computationally-easy inversion of the (noised) covariance matrix. Since GP regression is notoriously heavy on the computational side, the new method to speed it up is a significant contribution.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "2jI-Y9SPKUG", "reviewer": "Reviewer_a9uh", "summary": "\nIn this paper, the authors propose a binary tree-based kernel for Gaussian Process acceleration.  The proposed kernel enables a fast computation of the negative log-likelihood in $O(nq)$ time.", "questions": "\n1. Could the authors explain more about the time complexity regarding the above concern?\n\n2. It is better to provide a detailed comparison with a strong inducing pionts-based GP (Kapoor et al. 2021) for better evaluation.\n\n\nKapoor et al. SKIing on Simplices: Kernel Interpolation on the Permutohedral Lattice for Scalable Gaussian Processes. ICML 2021.", "limitations": "NA", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "\nPros.\n1.  The proposed binary tree-based kernel enables a fast (approximate) compuation of  $(K_{X,X} + \\lambda I)^{-1}$ and $|K_{X,X} + \\lambda I |$, which is important for GP accelration.\n\n2.  The fast computation algorithm (Algorithm 1) may be potential for acceleration in other tasks relying on the binary tree-based features.\n\n\n\nCons.\n1.  The time complexity of Algorithm 2 may be up to $O(2^q)$.\n From Eq.(3),  the number of $s \\in \\mathbb{B}^r$ is $2^r$, and it can be up to $2^q$. As a result, the number of partitions can be as \nlarge as $2^q$. The time complexity of the  for loop in step 7 in Algorithm 2  is then up to $O(2^q)$. This can be much larger than $O(nq)$. \nNote that $q=dp$,  $O(2^q)$ is quite large.\n\n2. The authors argue that inducing points-based GP acceleration is less competitive. However,  in the experiments,  a detailed comparison with strong inducing points-based methods is missing.   It is unconvinced to support the claim.\n\n3. Some symbols are confusing. \n(1). the index $i$ in $(X_{[s]})_i$ confuse with the index in $w_i$ in  Eq.(3).\n(2). The $\\boldsymbol{i}$ in step 8 in Algorithm 2 is not defined or used in Algorithm 2. \n(3).  $m$ in Line 237-251 is confused with the number of partitions in Definition 2 in Line 120.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept