Exploiting Hankel-Toeplitz Structures for Fast Computation of Kernel Precision Matrices
Authors: Frida Marie Viset
Anton Kullberg
Frederiek Wesel
Arno Solin
Venue: Accepted by TMLR ()

Abstract: The Hilbert-space Gaussian process (HGP) approach offers a hyperparameter-independent basis function approximation for speeding up Gaussian process (GP) inference by projecting the GP onto $M$ basis functions. These properties result in a favorable data-independent $\mathcal{O}(M^3)$ computational complexity during hyperparameter optimization but require a dominating one-time precomputation of the precision matrix costing $\mathcal{O}(NM^2)$ operations. In this paper, we lower this dominating computational complexity to $\mathcal{O}(NM)$ with no additional approximations. We can do this because we realize that the precision matrix can be split into a sum of Hankel-Toeplitz matrices, each having $\mathcal{O}(M)$ unique entries. Based on this realization we propose computing only these unique entries at $\mathcal{O}(NM)$ costs. Further, we develop two theorems that prescribe sufficient conditions for the complexity reduction to hold generally for a wide range of other approximate GP models, such as the Variational Fourier features approach. The two theorems do this with no assumptions on the data and no additional approximations of the GP models themselves. Thus, our contribution provides a pure speed-up of several existing, widely used, GP approximations, without further approximations

Published in Transactions on Machine Learning Research (MM/2002
Exploiting Hankel-Toeplitz Structures for Fast Computation of Kernel Precision Matrices
Frida Viset i
f. f.m.niseSShueelt..
Delft Center for Systems and Control
Delft University of Technology, The Netherlands
Anton Kullberg f
anton.kullbery Sliu.se
Department of Electrical Engineering Linköping University, Sweden
Frederiek Wesel
Sweesl@Sudelff..tt
Delft Center for Systems and Control Delft University of Technology, The Netherlands
Arno Solin
arno.solin@oalto,
Deportment of Computer Science
Aallo Universily, Finland
Reviewed on OpenReview: https: optp//veee.reeeee..mmeee
id= id= L0DdD
Abstract
The Hilbert-space Gaussian Process (HGP) approach offers hyperparameter-iddeee basis function approximation for speeding up Gaussian Process (GP) inference by projecting the GP onto M basis functions. These properties result in favorable data- mata-independent O compulational complexity during hyperparameter optimization but require dominating one-time precomputation of the precision matrix costing O(NM²) operations. In this paper, we lower this dominating computational complexity to O(NM) with no additional apprerimations. can do this because we realize that the precision matrix can split into a sum of Hankel-Toeplitz matrices, each having (MM) unique entries. Based on this realization we propose computing only these unique entries at O(N costs. Further, we develop two thcorems that prescribe sufficient conditions for the complexity reduction to hold generally for wide range of other approximate GP models, such as the Varialional Fourier Feature (VFF) approach. The two theorems do this with no assumptions on the data and no additional approximations of the GP models themselves. Thus, our contribution provides a pure speed-up of several existing, widely used, GP approximations, without further approzimations.
1 Introduction
Gaussian Processes (GPs, Rasmussen & Williams, 2006) provide flexible formalism for modeling functions which naturally allows for the incorporation of prior knowledge and the production of uncertainty estimates in the form of a predictive distribution. Typically a GP is inst by specifying a prior mean and covariance (kernel) function, which allows for incorporation of prior knowledge. When data becomes available, the GP can then be conditioned on the observations, yiclding a new GP which can be used for predictions and uncertainty quantificatio.
These authors contributed equally to this work. A reference implementation built top of GPJax available https://gtthub..cmm AOKaaaler////gggaaaaaaaeeee

---

Published in Transactions on Machine Learning Research (MM/2024)
10a
f- Wall-clock time (s)
HGP
Ours ((M3) ((M)
40
Number of DPS, M (x10²)
Figure An order of magnitude speed-up without any additional approximations: Wall-clock time to compute the precision matrix for an increasing number M of basis functions.
While all these operations have closed-form expressions for regression, the computations of the mean and covariance of the predictive GP require instantiating and inverting the kernel (Gram) matrix, which encodes pair-wise similaritics between all data. These operations require respectively ((N2) and O(N3) computations, where N is the number of data points. Furthermore, if one wishes to optimize the hyperparameters of the kernel function, which in GPs is typically accomplished by optimizing the log- marginall likelihood, the kernel matrix neods to be instantiated and inverted multiple times, further hindering the applicability of GPS to large-scale datasets.
The ubiquitous strategy reduce this computational complexity consists in approximating the kernel matrix in terms of Basis Functions (BFS), yielding what is essentially a low-rank or "sparse" approximation of the kernel function (see Rasmussen & Williams, 2006; Setthonero-Candde & Rasmussen, 2005; Snelson & Chahramani, 2005). Computational savings can then be achieved by means of the matrix inversion lemma, which requires instantiating and inverting the precision matrir (sum of the outer products of the BFs) instead of the kernel matrix at prediction, thereby lowering the computational complexity of inference to O(N + MM,, where M is the number BFs. If the number of BFs is chosen smaller than the number of data points in the training set (M < N), computational benefits arisc and the computational costs for hyperparameter optimization and inference are dominated by M²).
It is less widely known that this cost can be improved further. A notable BF framework is the HGP (Solin & Särkkä, 2020) which projects the GP onto dense set of orthogonal, This deterministic approximation is particularly altractive as it exhibits fast convergence guarantees to the full GP in terms of the number of BFs for smooth shift-invaria invariant kernels compared to other approximations. Furthermore, the fact that the BFs are hyperparameter-indeee inneeeerrndneeeeeentt speeds up GP hyperparameter optimization considerably, which in the HGP requires only ((M)) operations after one-time precomputation of the precision matrix, costing O(N M2). These favorable properties have ensured relatively wide usage of the IIP (see e.g. Svensson et al., 2016; Berntorp, 2021; Kok & Solin, 2018), and it is available in, e.g,, PyMC (Abril-Pla et al., 2023) and Stan (Riutort-Man Mayol et al., 2022). However, as argued by Lindgren et al. (2022), high number of BFs may be required for faithful approximation of the full model. Thus, the HGP is typically employed in applications where forming the initial O(NM2) projection docs not bocome too heavy.
pnperparantter-indpoodddd
BFs.
In this paper, we reduce this complexity to O(N)) with no additional approzimations (see Fig. 1), by exploiting structural properties of the precision matrix. We remark that these structural properties arise from the BFS that are used to approximate the kernel, not the kernel itself. Our contributions are as follows.
We show that the IGP precision matrix using basis functions defined on a (hyper-)cubical domain can be split in a sum of multilevel Slock-Hanker and a multilevel block-Toepha Toepliz matrix (Fig. 2), where each summand only has ((M) unique entries instead of ( This allows us to determine the elements of the full precision matrix at ((N)) instcad of O(N
This method does not only hold for HGP basis functions. We provide sufficient conditions on the BFS for this reduction in computational complexity to hold. Examples of other BES where the
2

---

Published in Transactions on Machine Learning Research (MM/ 2024)
ID 2D
3D
H
H & H (2)
H (1) @ H
(2)
@ H
(3)
The precision matrix for polynomial basis functions has a nested Hankel structure. The visualization of the matrix is proportionally darker as the logarithm each entry increases. The matrices are computed as the sum of all entries Hn for NI,...,N}. where the expression for H, is given below each matrix.
Figure 2:
conditions hold, and that therefore can be sped up using the same technique are polynomial, (complex) exponential and (coosiiua basis functions. We, therefore, enable the speed-up of other BF-based GP approximations such as variational Fourier features (Hensman et al., 2017).
In the experiments, demonstrate in practice that approach lowers the computational complexity of the HGP by order of magnitude on simulated and real data.
2 Background
A GP is a collection of random variables, any finite number of which have a joint Gaussian distribution (Ras- mussen & Williams, 2006). We denote a zero mean GP by f ~ (P((),((.))) where <(æ, x') RR x RD + R is the kernel, representing the covariance between inputs æ and x'. Given a dataset of input-output pairs (a Wa) GPs used for non-parametric regression and classification by coupling the latent functions I with observations through a likelihood model p(yl f] = II P(y: ((.)) [(=))).
For notational simplicity, we will in the following focus on GP models with a Gaussian (conjugate) likelihood, ys (f) The posterior GP, GP(u, (.), can be written down in closed form by
p. (x,)) = KT [K + 02 y,
(la) (1b)
(x,,, k( x.,) k. (K + o2I)-'k-
where K E NNN and k2 E R are defined elem as Kij == n( 2y) and [k.]: <(<;; x²) for i.j E 1.2,.... N, and k2 is defined similarly. Observations are collected into y RR, Due to the inverse of (K + T)) in the posterior mcan and covariance, the computational cost of a standard GP scales as O(N3), which hinders applicability to large datasets.
2.1 Basis Function Approximations
The prevailing approach in literature to circumvent the O(N5) computational bottleneck is to approximate the GP with a sparse approximation, using a finite number either inducing points or Basis Functions (BFS) (e.g., Rasmuussen & Williams, 2006; Caliionero-Candoe & Rasmussen, 2005; Hensman et al., 2017). The BF reprosentation is commonly motivated by the approximation
<(xxx') xe Φ(x) A d(x').
(2)
We use the notation from Solin & Särkkä (2020) to align with the next section. Here, Φ(-): : R - M are the
BFs, Φ(-) = (0)()))22()).....)))))))))))(()) Ø2(+),... (s)()))).....))))))))) (-) Further, A E MXM are the correspondigg BF weights. Combining
this approximation with the posterior GP, Fq. (1), and applying the Woodbury matrix inversion lemma yields
(x(x)) = ((x.))
( 2 + - & y. ④ + + - A -1) -1)''(((1)
(3a)
E.,,, x()) = o' T
(3b)
3

---

Published in Transactions on Machine Learning Research (MM/2024)
Here, & E RNXM is commonly referred to the regressor matrix and is defined as Pt.: = o(x,)" Further, RMX is the precision matriz which is a central component of the following section. Computing this precision matrix requires (() operations. With this approximation, computing the posterior mean and covariance in Eq. (3) requires instantiating and inverting (TT + + o²A-1) which can be performed with O(NM² M³) operations. If have more samples than BFS (which is the required condition for BFS to give computational savings), i.e., N 2 M, the overall computational complexity is then of O(NM²), Lee, the inversion costs are negligible as compared to the cost of computing the precision matrix.
+ T T E
3 Methods
Our main findings are in the form of two theorems (Theorems and 3.4). These theorems prescribe the necessary conditions that BF expansions need fulfill to be able reduce the computational complexity of computing the precision matrix + 4 from ( O(NM), applicable multiple previous works that rely parametric basis functions (incl. Lázaro-Gredilla Gredilla et al., 2010; Hensman et al., 2017; Solin & Särkkä, 2020; Tompkins & Ramos, 2018; Dutordoir et al., 2020). Since our contribution reduces the cost of computing the precision matrix to MM, it follows that the computational complexity of computing the posterior mean or covariance is also of O(NM). Further, both of the thcorems reduce the memory scaling from (MM2) to O(M). Note that these reductions are without approximations, only relying on the structural properties of the considered modcls.
In the following, we assume that the kernel is lensor product kernel (Rasmussen & Williams, 2006), i.., x(x,2") = 1I. (2((), x(d)), where (.,)) is the kernel along the t dimension. Then, each component of the kernel is approximated with ma BFs such that
(1((),,(0))
(x T Add) (x(4))) (d)
(4)
where (+) R + Rm are the BEs [0(4] 92
IollT along the dth dimension and (d) E maaxa
contains the associated weights. The full kernel can then be approximated as
x(xx I1 (1(())
A(d) (2(0))
(5)
where in this case :(:):Re R. and A E RRM are
Φ(x) = - (1)00 d=1 (d)
(6a)
(6b)
Here, M = -T ma is the total number of BFS. Given this decomposition, the precision matrix can be expressed as
+ T - >. ((x (())))(())))) = EN ( (500
(7)
This decomposition of the precision matrix is key in the following and we will primarily study the individual products (f()) (x."")
provide some intuition, (=)""))" where certain structure may appear that is exploitable to benefit. To consider the precision matrix for polynomial BFs in ID, 2D, and 3D (see Fig, 2). The 1D case (left) has a Hankel structure and the 2D (middle) and 3D (right) cases have 2-level and 3-level block Hankel structure, respectively. It these types of structures that allow a reduction in complexity without approximations. Next, we provide clear technical definitions of the matrix structures and then proceed to state our main findings. All the discussed matrix structures are also visually explained in Table Al in App. E.
3.1 Hankel and Toeplitz Matrices
An m x m matrix H has Hankel structure iff it can be expressod using a vector y (2m-1) containing all the unique entries, according to
72
Ym
73
7m+1
H
(8)
7
7m+1
72m
4

---

Published in Transactions on Machine Learning Research (MM/2024)
Σ HO
Σ TO)
Figure 3: The precision mairix for sinusoidal basis functions in dimension has neither Hankel nor Toeplitz structure. However, it can be decomposed into sum of two matrices, where one has Hankel structure, and one has Toeplitz structure. Here, 49 BFS are placed along one dimension.
In other words, each element of the Hankel matrix on the jth row and the jth column is given by Y----1- Similarly, an
matrix has Toeplitz structure iff it be expressed using vector Y R (22mm1) such that each element on the ith row and the jth column is given by 7-jtm- See Fig. 3 for an example of what a Hankel and Toeplitz matrix visually looks like.
We define a matrix (D) as a D-level block Hankel matrix if it can be expressed as
H (D11)
HDD-1)
Hmp (D-1)
H (DD1)
H
(D11)
H mpti
(D))
(9)
H
(D-1) mp+1
(D-1) 2mp-1
where H (D-) is a block Hankel matrix, and H is simple Hankel matrix. Each submatrix H(! for d E {1..... D} can be indexed by row and column id.ja E mm....mpJ which again yields submatrix H!d-1) which can be indexed in the same manner. Therefore, an equivalent definition of block Hankel matrix is that each individual entry of H given by a set of indices i1, i1.J1,...,, iD. is given by the entry "k.....p of a tensor 7 E K X2xixKxx..xxx KD, where ka = id + ja 1, and K₃ = 2m - 1. Similarly, (D-1) a block Toeplitz matrix each entry f1......... :........ip..................... ip.jo is given by the entry kkp...kpp of tensor
7 E R Kax...X...XXD where ka = ma + la ja, and Ka 2ma 1.
3.2 Block Toeplit Matrices
In addition block Hankel and block Toeplitz structures, we require a slightly more general but highly related structure. We call this structure block Hankel- Tooplitz structure and define it as follows. matrix (DD) has D-level block Hankel Toeplitz structure, if the matrix is defined either
(D-1)
G (D-1)
(D-11 1)
Gmp
G
(D-1)
G
D-11
(D-1)
G(D)
mp+1
(10)
(D-1)
(D-1)
(--1)
Gmp
mp+i
2mp
if level is Hankel, or as
(D-1)
(D-1)
G (D-1)
Gmp
@mptl
G
(-11)
G (D-1)
G(D)
(11)
(D-1)
(D-1) mptt
(DD11)
G2mp
if level D is Toeplitz. Further, G/D-1) are D level block Mankel-Toeplice matrices if D >>> 2, and
simple Hankel Toeplitz matrix if D 1. Each submatrix G(d) for d (1,..., D) be indexed by

---

Published in Transactions on Machine Learning Research (MM/2024)
row and column dd {I,..., ma) which again yields submatrix defined either Gld-1) or G
matia-ja
(depending on whether the devel has Hankel structure as in Eq. (10) or Toeplitz structure as in Eq. (11)).
Each entry in a block Hankel Toeplitz matrix can also be expressed by the entry Aa of a tensor y E R RKKi Xxx.xxxXXX where Ka 2md 1, and
id ja 1, if level d is Hankel ma + ia ja: if level d is Toeplitz.
ka
(12)
A crucial property of the block Hankel Toeplitz structure is the preservation of structure under addition. Assume that A and B are two block -Tcpplitz matrices and that they are structurally identical, in the sense that they have the same number of levels, the same number of entries each block, and each level shares either Toeplitz or Hankel properties. Then, let each entry of A and B be given by Ag...k..... for tensor Rm....m and K2..k for 8 RM,..., respectively, with ka defined in Eq. (15). Each entry in A B is then given by the sum of the entries akk....ppp + fk.......p Thus, the sum of two block Hankel-Toeplnn matrices with identical structure is also a block Hlankel-Tooplee matrix. By associativity of matrix addition, a sum Ex= Gn of N Hankel Toeplitz matrices {G....., Gr3 with identical structure is therefore itself a Hankel Toeplitz matrix.
3.3 Kronecker Products of Hankel-Toeplitz Matrices and Block Hankel-Toeplitz Matrices
A special case of a class of matrices G which has block Hankel Toeplitz structure are Kronecker products of D Hankel and Tocplitz matrices (G()),... ...., GDD), i.c.
G
G(d) (
G
2 D(D)
(13)
d=11
An equivalent definition of the Kronecker product gives each entry on the 7th row and jth column in the block Toeplitz matrix G as an expression of the entries on the is" row and jath column of each matrix ((d) according to
Gzj= = -I2m (((a)
(14)
Note that there is one-to-one map between each index i, and the index sets (+i,... "ii..-,;ip) in} and (f1...., in As each matrix G(d) has Hankel Toeplitz structure, the entries can equivalently be defined by vector 7 with 2m - 1 entries. Each entry Gig is therefore given by
G = - The -on
(d)
(15)
where ka is defined in Eq. (12), and whore 7 = =02 d- ((()) is a rank-1 tensor with I[2 (2ma - 1) elements.
3.4 Main Results
We are now ready to state our main findings. The following two theorems rely on the product decomposition Eq. (5) and study each dimension separately, as is evidently possible from Eq. (7). Our first thcorem generalizes result by Greengard et al. (2023) regarding complex exponential basis functions. Our first theorem establishes that if the product (5((((()))))) (x(h)) has Hankel or Toeplitz structure, the resulting precision matrix only has I12,,2ma, 1) unique entries, reducing the computational complexity of instantiating it from O to O(NM). We formalize this in the following theorem.
Theorem 3.1. If the matrir
G = (1)))))))0000000001
(16)
is a Hankel or Toeplitz matriz for all ((d) E R along each dimension d, the information matrir T will be multi-levee block Hankel Toeplitz mabrix, and therefore have l (2ma - 1) unique entries.

---

Published in Transactions on Machine Learning Research (MM/2024)
T
E HO
E SHO
Σ (H
ΣTO
Figure 4: The precision matrix for sinusoidal BFs two dimensions has neither Hankel nor Toeplitz structure. However, it can be decomposed into 2D 4 matrices, which each have block Mankel-Toophee structure. Here, 7 BPS are placed along each of the two dimensions, giving total of 49 BFs.
Proof. Assume that the matrix (I1(") (xn ( (2²") is Hankel Toeplitz. The precision matrix can then be expressed as
T + = En oZ (xn (210) (c)) = =
(17)
where the matrix _2_ ((d)( (In (d) is multi-level Hankel or Toeplitz by definition, see Sec. 3.3. Further, the sum
of several D-level block Hankel Hanpel-Toeplizz matrices itself a D-level block Hankel Hankel-Toeplitz matrix, see Sec. 3.2. Since each matrix G(d) (d)) has at most (2ma- - 1) unique entries, the matrix T T >n Gldd (zn therefore has most M I (2ma - 1) unique entries.
The preceding theorem holds true for, for instance, polynomial and complex exponential BFs, which we establish in Corollaries 3.2 and 3.3.
Corollary 3.2. The precision matris for polynomial BFS defined by
1 (x)c)) = (2/1))80-11
(18)
can be represented by tensor with I 2m 1 entries.
Proof. See App. A for proof.
0
Corollary 3.3. The precision matrir for compler exponential BFs defined by (x(x)= = exp(inj] x) -T exp(irjjaaa,
(19)
can be represented by tensor with IT 2md 1 entries.
Proof. See App. B for a proof.
For some BFs, the structure the product ((d) (xy (()) (b( (x(h (zh")) is more intricate, but still of a favorable, exploitable nature. This is clearly evident from Fig. 4, where the precision matrix for sinusoidal BPS two dimensions is visualized. In particular, for some BFs, the product is the sum of a Hankel and Toeplitz matrix, such that the precision matrix only has I 3ma unique entries, again reducing the computational cost of computing it to O(NM). We formalize this in the following theorem.
Thcorem 3.4. If the product (idd) (d)) (1)")) (x.)))) (IT is the sum of Hankel matrix denoted (( (d) and Toeplitz matrir (z(), ((()), and there exists a function g(d)
(ka, 2il)) such that (1)))
G (2(
g(āa + ja. x()) and Giads (a())) (())) -g ja. x 1.l)), all entries in the precision matria can be represented (dd).(-1)
by tensor (k1, k2,....kpp) kp) with [b 3mg entries.
7

---

Published in Transactions on Machine Learning Research (MM/2024)
Proof. The precision matrix can in this case be expressed as
4 En of 8z (al) (21) (
Σ &2-1 (( (1)) ((a,)-1((h1))) (z.h))
Epe (Tw= (3)) En 8 ((d))e)))) (1))
(20)
where ep = {ep p...,pp dD) E 5 and S = -1 is set containing 2 elements. Each of the 2D matrices [Tl=1 0 Em= -2_ G(d)(( (.x) ( ) is now the Kronecker product between D Hankel or Toeplitz matrices.
The entries of C can be expressed dementtwwse as
20 c=== = Ep=
(T ( Em
E Ile glia + ep Ja: (a.
(21)
If we define a tensor y as
Yk.... .h..kk = >n II (kk, on
(22)
.22mm 1, 2m,, each entry the precision matrix Ct can be expressed as for indices ki = - md. 2 - My,
2n (IT2 c
(d) Cp Visteg ie (
(23)
As each sum kn = e
(d) ja is an integer between 1-m and 2ma, the tensor y will have I 3m, entries.
The preceding theorem applies to, for instance, the BFs in an HGP (Solin & Särkkä, 2020) defined roctangular domain, formalized in Corollary 3.5. It further holds for multiple other works using similar BFs, formalized in Corollary 3.6. We remark that for D = 1, we require ma > 3 for any savings to take effect, whereas for D > 1, ma 2 2 suffices.
Corollary 3.5. The precision matriz in an HGP defined on rectangular domain [- L1) x ... x LDL, Lpl can be represented by tensor with I 3ms entries.
Proof. See App. C for a proof.
Corollary 3.6. The precision matriz in a GP approcimated by sinusoidal and cosine BFs with frequencies on grid (such as the regulor Powrier features described in Hensman el al. (2017) and Wahls el al. (2014), the Fourier approzimations
periodic kernels described in Tompkins & Ramos (2018), the quadrature Fourier features described in Mutng & Krause (2018), the puuspaceddverssin of sparse spectrum BES described in Lázaro- Gredila et al. (2010), or the one-dimensional special case of Dutordoir et al. (2020)) can be represented by a tensor with I1 3ma entries.
Proof. See App. D for a proof.
3.5 Outlook and Practical Use
Both Theorems 3.1 and 3.4 reduce the computational complexity of calculating the entries of the precision matrix Eq. (7) from O(N to O(NM). This enables us to scale the number of BFs significantly more than previously, before running into computational, or storage related, problems. For clarity, the standard IIP is given in Alg, with the original approach in red and our proposed approach in blue. As compared to the standard (offline) HGP, the only change we make is the computation of the precision matrix.
We remark that computing the posterior mean and variance of the BF expansion now costs O(NM M3), now possibly dominated by the M term. This cost also appears in the hyperparameter optimization, choose to optimize the MIL. in Solin & Särkkä (2020). Remedies for this cost out of scope of this
8

---

Published in Transactions on Machine Learning Research (MM/2024)
Algorithm 1 Sketch of an algorithı for Hilbert GP learning and inference. The original approach by Solin & Särkkä (2020) in red, our proposed approach in blue.
Input: Data as input-output pairs {( (e)
test inputs æx, number of basis functions M
Compute + + at cost. O(NMM)
Eq. (7) Eq. (22)
Compute y at cost M(N
Construct + + using Y at cost ((M²)
repeat
Optimize Marginal Log- Likeeikedd (MLL) w.r.t. hyperparameters at cost O(M3)
until Convergence
Perform GP inference using the pre-calculated matrices. This entails computing the posterior mean and covariance, at computational cost of ((M3).
Number of BFs (M)
HGP
Ours
6,400
3,500
900
73.46
Longitude
Wall-clock
73.42
10a
7.45 7.65 7.0
7.45 7.6 7.55 7.8
7.45
7.6 7.56 7.0
Latitude (*)
Latitude (?)
Latitude (?)
-0.1
0.1
0.2
20
40
80
NLPD
(a) Predictive means indicated by color intensity with transparency
(b) Wall-clock time to compute precision matrix
proportional to the predictive variance. Increasing number of nFs by sequentially including each data point from left to right.
NLPD.
Figure 5: Our proposed computational scheme reduces the computation time for datasets with tighffrequency variations, these require many BFs to achieve accurate reconstruction. This underwater magnetic field has lower Negotive Log Predictive Density (NLPD) with a large amount (6400) compared to a smaller amount (400) of BFs. For 6400 BFs, our computational scheme reduced the required time to compute the precision matrix from 2.7 hours to 1.7 minutes.
paper, but could potentially be reduced through the use of efficient approximate matrix inverses and trace estimators (see, e.g., Davies, 2015). Another consequence of these theorems is that multi-agent systems that collaborate to learn the precision matrix such as Jang et al. (2020), Visct et al. (2023) or Pillonetto et al. (2019) can do this by communicating O(M) bits instead ((M2) bits.
4 Experiments
We demonstrate the storage and computational savings of our structure exploiting scheme by means of three numerical experiments. The experiments demonstrate the practical efficiency our scheme using the HGP. We reitcrate that the savings are without additional approximations and the posterior therefore enactly equal to that of the standard HIG.. Further, HIPSS adhere to Theorem 3.4, this demonstration is a representative example of the speedups that can be expected using, e.g., regular Fourier features (Hensman et al., 2017), BF expansions of periodic kernels (Dutordoir et al., 2020). Since Theorem 3.1 requires computing and storing only 2D M components whereas Theorem 3.4 requires 3D M, our experiments demonstrate a practical upper bound on the storage and computational savings, a "worst case".
Our first experiment demonstrates the computational scaling of our scheme on simulated 3D datasct. Secondly, we consider a magnetic field mapping example with data collected by an undervater vessel, as an application where the high-frequency content. of the data requires large amount of BFS to reconstruct the
9

---

Published in Transactions on Machine Learning Research (MM 2024)
HGP
Ours ((M²)
((M)
Wall-rlock time
0.1
3
(a) Full GP
(b) HGP / Ours
Number of BFs, M (x (x100)
(c) Time to compute the precision matrix for increasing M
1500
3000 [mm]]
Figure These experiments recover the results from Solin & Särkkä (2020) exactly for predicting yearly precipitation levels across the US, and measure the wall-clock time needed by our proposed computational
scheme. The HGP efficiently approximates the full GP solution using my m2 = 45, totaling M 2025 BFs.
field. Thirdly, a precipitation dataset is used, mirroring an example in Solin & Särkkä (2020), improving the computational scaling in that particular application even further than the standard HGP. Our freely available IIP reference implementation along with hyperparameter optimization is written for GPJax (Pinder & Dodd, 2022). All timing experiments are run on an HP Elitebook 840 G5 laptop (Intel i7-8550U CPU, 16GB RAM). For fair comparison, we naively loop over data points, to avoid any possible low-level optimization skewing the results.
Computational Scaling We compare the time necessary for computing the precision matrix for N 500 data points for the standard HGP well as for our structure exploiting scheme. The results are presented in Fig. 1 for an increasing number of BFs. After M = 14000, the HGP is no longer fcasible due to memory constraints, whereas our formulation scales well above that, but stop al M = 64000 for clarity in the illustration. Clearly, the structure exploitation gives a significantly lower computational cost than the standard HGP, even for small quantities of BFs. Further, it drastically reduces the memory requirements, where for M 14000, the HGP requires roughly 1.5 GB for storing the precision matrix, while our formulation requires roughly 2.8 MB using 64-bit. floats. This makes it possible for us to scale the number of BFs significantly more before running into computational storage restrictions. It is noteworthy that even though the implementation is in high-level language, we still see significant compulational savings.
Magnetic Field Mapping As HCPS are commonly used to model and estimate the magnetic field for, e.g., mapping purposes (Solin et al., 2015; Kok & Solin, 2018), we consider a magnctic field mapping example and demonstrate the ability of our computational scheme to scale HGPs to spatially vast datasets. The data was gathered sequentially in lawn-mower path in region approximately x km large by an underwater vessel outside the coast of Norway (d=2 with roughly 50/50 split, deterministically split by grid pattern, to ensure reasonable predictions in the entire data domain, see App. F.2 for more details. We vary the amount of BFS and compare the time required to sequentially include each new data point in the precision matrix as well as the Negative Log Predictive Density (NLPD). As the underwater magnetic field covers a large area, a large number of BFs are required to accurately represent the field, see Fig. 5a where the details of the predicted magnetic field captured more accurately for an incrcasing number of BFs. This also apparent from the decrcasing NLPD as the number of BFS increases, see Fig. 5b. At 6400 BES, the necessary computation time is several orders of magnitude lower for our approach compared to the standard HGP.
= 2 and N 1.39 million). The data was split into training set and test set
U.S. Precipitation Data We consider a standard precipitation data set containing US annual precipitation sunmaries for year 1995 (d = 2 and N = 5776) (Vanhatalo & Vehtari, 2008). We exactly mimic the set in Solin & Särkkä (2020) and primarily focus our evaluation on the calculation of the precision matrix. The time for computing the precision matrix is visualized in Fig. 6c, where our approach clearly outperforms the standard HGP. The predictions on a dense grid over the continental US can be found in Figs. 6a and 6b, where the HGP manages to capture both the large-scale as well as the small-scale variations well.
10

---

Published in Transactions on Machine Learning Research SMMM/024]
5
Related Work
The computational scheme that detail here can be used to speed up range of approximate GP and kernel methods with BFs that satisfy the Hankel Toeplitz structure we use (see e.g. Tompkins & Ramos, 2018; Solin & Särkkä, 2020; Hensman et al., 2017). It is also applicable to Lávaro-Gredille et al. (2010) in the special case where the considered frequencies of the BFs are equidistant, as well as Dutordoir et al. (2020) when the input space is onc-dimensional. However, there is also a wide range of GP approximations that do not have the Hankel Toeplitz structure required for Theorem 3.1 or Theorem 3.4 to apply. While the structure we exploit in Theorem 3.4 is apparent in the quadrature Fourier feature approach of Mutný & Krause (2018) when the frequencies are on a structured Cartesian grid, other quadrature- methods are not possible to speed-up in similar ways. The most well-known method is the random Fourier feature approach (Rahimi & Recht, 2007), where speed-up is not possible as the frequencies are sampled at random. Similarly, the Gaussian quadrature approaches of Dao et al. (2017); Shustin & Avron (2022) and the random quadrature approach of Munkhoeva et al. (2018) do not constrain frequencies to a regular grid and therefore do not have the Hankel-Toepae structure required by Theorem 3.4.
As was pointed out by Muinonoco-Candela & Rasmussen (2005), inducing point approaches can also be viewod as BF approximations. The inducing variable approaches essentially summarize the data in a set of inducing variables (Snelson & Ghahramani, 2005; Seeger et al., 2003; Csató & Opper, 2002), which most commonly represent function values at a certain set of inputs, even though other choices are possible (Lázaro-Gredilla & Tigueiras-Vidal. 2009). Another closely related approach which also uses inducing points, with similar computational complexity of O(N is the variational GP pioneered by Titsias (2009). In light of the BF viewpoint, all the aforementioned approaches therefore also involve computing a precision matrix, but even for ordered inducing points on a grid, the particular Hankel-Toeplitz structure we exploit does not exist in general.
A range of related work stacks additional approximations top of the sparse basis function cennttioninducing /inducing point approximations (Yadav et al., 2021; Izmailov et al., 2018; Hensman et al., 2013). The purpose of these additional approximations is typically to either reduce the computational complexity of computing the precision matrix for example by using BFs with compact support, or by implementing an approximate algorithm for inverting the precision matrix. A well-known structure exploiting method is Structured Kernel Interpolation (SKI) (Wilson & Nickisch, 2015; Yadav et al., 2021; Izmailov et al., 2018), which approximates the precision matrix through cubic interpolation between inducing points on a regular grid. An alternative approach that improves the computational complexity of the variational GP, is the Stochastic Variational GP (SVGP) (Hensman et al., 2013). The SVGP is the de facto standard approach for large-scalo GPS, due to the possibility of utilizing mini-batching for training, greatly speeding up hyperparameter (and variational parameter) learning (Hensman et al., 2013; 2015). with implementation available in, e.g., GPyTorch and GPflow (Gardner et al., 2018; de G. Matthows et al., 2017).
While this paper focuses on BF approximations to GPs, there is a lot of work on approximating exact GP regression through means of, e.g., Conjugate Gradient (CG) descent (Gibbs & MacKay, 1996; Artemev et al., 2021; Gardner et al., 2018). These methods view GP regression the solution to lincar system of equations and seek to solve this approximately. This typically reduces the cost of exact GP regression from O(N3) to O(N2),, where I is the number of CG iterations. Davies (2015) uses CG as the driving scalability mechanism, but further reduces the computational complexity to O(N through use of M-efficient kernels, which essentially constitute the approaches from Muinonero-Candell & Rasmussen (2005) as well as compact kernels (see, e-g., Kullberg et al., 2021; Buhmann, 2003; Wu, 1995). CG has also been combined with SKI in the KISS-GP framework of Wilson & Nickisch (2015) to reduce the computational complexity of (mean) inference to ONN + M log M). Pleiss et al. (2018) extended KISS-GP with the Lanczos algorithm to reduce the complexity of computing the predictive variance to O(k) after pre-computation, where k is the number of Lanczos iterations. Recently, Stochastic Gradient Descent (SGD) was introduced as an alternative to CG (Lin et al., 2023a;b), potentially with better performance for ill-conditioned datasets. Technically, any of these ideas are straightforward to include in the HGP as well and could palentially be used to reduce the ((3)) complexity of inference and hyperparameter learning.
Other previous work has discovered that functions of difference matrices (a special case of Toeplitz matrices) are also sums of Hankel and Toeplitz matrices in the one-dimensional case, and Kronecker products of these
11

---

Published in Transactions on Machine Learning Research (MM/2024)
in the two-dimensional case (Strang & MacNamara, 2014). Since our matrices are also Kronecker products of Hankel and Toeplitz matrices, we end up exploiting similar structures as Strang & MacNamara. (2014) arising in a different situation.
6 Conclusion
Our contribution delails a computational approach for exploiting Hankel and Toeplitz structures that appear in multiple BF approximation schemes to kernels for GPs. These structures allow us reduce the computational complexity of computing the corresponding precision matrix from O(N to O(N without further approximations. Further, approach reduces the storage requirement for containing all necessary information about the posterior to make predictions from ((M2) to ((M). The Hankel and Toeplitz structures appear because of the particular BFs that are used to approximate the kernel, not the kernel itsolf. The reduced computational and storage requirements are particularly beneficial in the HGP where more BES allow us to capture higher frequencies of the kernel spectrum, otherwise unattainable without significant computational resources. We foresce that our contribution will allow HGPS to tackle larger problems without the need for extensive specialized hardware, opening up approximate audience. Future work could investigate if the results can be generalized to wider ranges of BFS, which is easily verified through Theorems 3.1 and 3.4. Another potential is to investigate ways of approximately decomposing an already known precision matrix into Tonpel-Topplitz matrices which does not admit an analytical such decomposition. This would yield the same computational benefits but with some potential loss of accuracy.
learning and inference for wider
A reference implementation built on top of GPJax available at http:///ghhabbcoom/AAkkkllllgggp structure.
Impact Statement
This work develops numerical methods for the wide field of machine learning, where the goal is to make existing methods more compute-efficient and open new avenues for developing future methods. There are many potential uses and therefore socictal consequences of such methods, none of which we see the need to specifically highlight here.
Acknowledgements
We would like to express our gratitude to the supervisors (Gustaf Hendeby, Isaac Skog, Kim Batselier, Manon Kok and Rudy Helmons) of the three PhD candidates (Anton Kullberg, Frederick Wesel and Frida Viset) involved in this project. They secured the funding that made this research possible. Their support in providing the necessary resources and their encouragement for our development as independent researchers have been invaluable. Their contributions have thus indirectly shaped this work, and we are grateful for their continued guidance and support. We would like to thank the anonymous reviewers for their numerous suggestions which have greatly improved the quality this paper. Frederiek Wesel, and thereby this work, is supported by the Delft University of Technology AI Labs program. Arno Solin acknowledges funding from the Research Council of Finland (grant id 339730). The underwater magnetic field data used were collected by MARMINE/NTNU research cruise funded by the Research Council of Norway (Norges Forskningsråd, NFR) Project No. 247626/ 770 O30 and associated industrial partners. Occan Floor Geophysics provided the magnetometer that was used for magnetic data acquisition and pre-processed the magnetic data. The authors declare no competing interests.
References
Oriol Abril-Pla, Virgile Andreani, Colin Carroll, Larry Dong, Christopher J Fonnesbock, Maxim Kochurov, Ravin Kumar, Junpeng Lao, Christian C Luhmann, Osvaldo A Martin, et al. PyMC: modern, and comprehensive probabilistic programming framework in Python. PeerJ Computer Science, 9:e1516, 2023.
12

---

Published in Transactions on Machine Learning Research (MM/2024)
Artem Artemev, David R. Burt, and Mark van der Wilk. Tighter bounds the log marginal likelihood of
Gaussian process regression using conjugate gradients. In Proceedings of the 38th Internationa Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pp. 362 372. PMLR, 2021.
Karl Berntorp. Online Bayesian inference and learning of Gaussian-process state space models. Automatica, 129:109613, July 2021.
Martin D. Buhmann. Radial Basis Functions: Theory and Implementations. Cambridge University Press, 2003.
Zhongming Chen, Kim Batselier, Johan A. K. Suykens, and Ngai Wong. Parallelized tensor train learning of polynomial classifiers. IEFF Transactions on Neural Networks and Learning Systems, 29(10):4621 2018.
Lchel Csató and Manfred Opper. Sparse on-line Gaussian processos. Neural Computation, 14(3):641-668, March 2002.
Tri Dao, Christopher M De Sa, and Christopher Ré. Gaussian quadrature for kernel features. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
Alexander James Davies. Effective Implementation of Gaussian Process Regression for Machine Learning. PhD thesis, University of Cambridge, 2015.
Alexander G. de G. Matthews, Mark van der Wilk, Tom Nickson, Keisuke Fujii, Alexis Boukouvalas, Pablo Kcón-Villageee Zoubin Ghahramani, and James Hensman. GPflow: A Gaussian process library using TensorFlow. Journal of Machine Learning Research, 18(40):1-6, 2017.
Vincent Dutordoir, Nicolas Durrande, and James Hensman. Sparse Gaussian processes with spherical harmonic features. In Proceedings of the 37th International Conference on Mochine Learning, volume 119 of Proccedings of Machine Learning Rescorch, pp. 2793-2802. PMLR, 2020.
Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. GPyTorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration. In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
Mark N. Gibbs and David J. C. MacKay. Efficient implementation of Gaussian processes for interpolation. Technical Report, University of Cambridge, 1996.
Philip Greengard, Manas Rachh, and Alex Barnett. Equispaced fourier representations for efficient gaussian process regression from billion data points. arXiv, 2023.
James Hensman, Nicolò Fusi, and Neil D. Lawrence. Gaussian processes for big data. In Proceedings of the 29th Conference on Uncertainty in Artificial Intelligence, UAI, pp. 282-290. AUAI Press, 2013.
James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational Gaussian process classification. In Proceedings of the 18th International Conference on Artificial Intelligence and Statistics, volume 38 of Proceedings of Machine Learning Reseorch, pp. 351-360. PMLR, 2015.
James Hensman, Nicolas Durrande, and Arno Solin. Variational Fourier features for Gaussian processes. The Journal of Machine Learning Research, 18(1) (11))55537-5588 January 2017.
Pavel Izmailov, Alexander Novikov, and Dmitry Kropotov. Scalable Gaussian processes with billions of inducing imputs via tensor train decomposition. In Proceedings of the 21st International Conference Artificial Intelligence and Statislics, volume 84 of Proceedings of Machine. Learning Research, pp. 726-735. PMLR, 2018
Dohyun Jang, Jachyun Yoo, Clark Youngdong Son, Dabin Kim, and H. Jin Kim. Multi-robot active sensing and environmental model learning with distributed Gaussian process. IEEE Robotics and Automation Letters, 5(4):5905-50) 5912, 2020.
13

---

Published in Transactions on Machine Learning Research (MM/2024)
Diederik P. Kingma and Jimmy Ba. Adam. Adam: A method for stochastic optimization. In Proceedings of the 3rd International Conference on Learning Representations (ICLR), 2015.
Manon Kok and Arno Solin. Scalable magnetic field SLAM in 3D using Gaussian process maps. In Proceedings of the 21st International Conference. on Information. Fusion (FUSTON), pp- 1353 1353-1360, 2018.
Anton Kullberg, Isaac Skog, and Gustal Hendeby. Online joint state inference and learning of partially unknown state space models. IEEE Transactions on Signal Processing, 69:4149-411 44611, 2021.
Miguel Lázaro-Gredile and Aníbal Figueiras- Vidal. Inter-domain Gaussian processes for sparse inference using inducing features. In Advances in Neural Information Processing Systems., volume 22, 2009.
Miguel Lázaro-Gree Gredilla, Joaquin Quiñonero- Gandela Carl Edward Rasmussen, and Aníbal R. Figueiras- -idal.. Sparse spectrum Gaussian process regression. Journal of Machine Learning Research, (633)111851188, 2010.
Jihao Andreas Lin, Javier Antorán, Shreyas Padhy, David Janz, José Miguel Hernández-Lobato, and Alexander Terenin. Sampling from Gaussian process postcriors using stochastic gradient descent. In Advances in Neural Information Processing Systems, volume 37, 2023a.
Jihao Andreas Lin, Shreyas Padhy, Javier Antorán, Austin Tripp, Alexander Terenin, Csaba Szepesvári, José Mignel Lonnataez-obato, and David Janz. Stochastic gradient descent for Gaussian processes done right. MrXXiv2330020055 2023b.
Finn Lindgren, David Bolin, and Håvard Rue. The SPDE approach for Gaussian and Gaussann fields: 10 Years and still rumning. Spatial Statistics, 50:100599, August 2022.
Marina Munkhoeva, Yermek Kapushev, Evgeny Burnaev, and Ivan Oseledets. Quadrature- based features for kernel approximation. In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
Mojmír Mutný and Andreas Krause. Efficient high dimensional Bayesian optimization with additivity and quadrature Fourier features. In Advances in Neural Information Processing Systems, volume 32, pp. 9019- 9030, 2018.
A. Novikov, M. Trofimov, and I. Oscledcts. Exponential machines. Bulletin of the Polish Academy of Sciences Technical Sciences, 66(6 (Special Section on Deep Learning: Theory and Practice)):7899-99 2018.
Gianluigi Pillonetto, Luca Schenato, and Damiano Varagnolo. Distributed multi agent Gaussian regression via dimensionnn approximations. IEEE Tronsactions on Pattern Analysis and Machine Intelligence, 41 (9):20099-211 2019.
Thomas Pinder and Daniel Dodd. GPJax: A Gaussian process framework in JAX. Journal of Open, Source Software, 7(75):4455, 2022.
Geoff Pleiss, Jacob Gardner, Kilian Weinberger, and Andrew Gordon Wilson. Constant-tine predictive distributions for Gaussian processes. In Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 4114 4123. PMLR, 2018.
Joaquin QuifioneroCannella and Carl Edward Rasmussen. A unifying view of sparse approximate Gaussian process regression. Journal of Machine Learning Research, 6(65):1939-1959, 2005.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems, volume 20. Curran Associates, Inc., 2007.
C.E. Rasmussen and C.K.I. Williams. Gaussion Processes for Machine Learning. Adaptive Computation and Machine Learning series. MIT Press, 2006.
Gabricl Riutort-Mayol, Paul-Christian Bürkner, Michael R. Andersen, Arno Solin, and Aki Vehtari. Practical Hilbert space approximate Bayesian Gaussian processes for probabilistic programming. Statistics and Computing, 33(1) 1777, December 2022.
14

---

Published in Transactions on Machine Learning Research (MM/2024)
Matthias W. Seeger, Christopher K. I. Williams, and Neil D. Lawrence. Fast forward selection to speed up sparse Gaussian process regression. In International Workshop on Artificial Intelligence and Stalistics, volume 9 of Proceedings of Machine Learning Research, pp. 254 261. PMLR, 2003.
Paz Fink Shustin and Haim Avron. Gauss- Regss-Tegendre features for Gaussian process regression. Journal of Machine Learning Research, 23(92):1-47, 2022.
Edward Snelson and Zoubin Ghahramani. Sparse Gaussian processes using pseudo-inputs. In Advances in Neural Information Processing Systems, volume 18. MIT Press, 2005.
Arno Solin and Simo Särkkä. Hilbert space methods for reduced-rank Gaussian process regression. Statistics and Compuling, 30(2):419 446, March 2020.
Arno Solin, Manon Kok, Niklas Wahlström, Thomas B Schön, and Simo Särkkä. Modeling and interpolation of the ambient magnetic field by Gaussian processes. IEEE Transactions on Robolics, 34(4):1112-113 1127, September 2015.
Gilbert Strang and Shev MacNamara, Functions of difference matrices are Toeplitz plus Hankel. SIAM Review, 56(3):525-546, 2014.
Andreas Svensson, Arno Solin, Simo Särkkä, and Thomas Schön. Computationally efficient Bayesian learning of Gaussian process state space models. In Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, volume 51, pp. 213 221. PMLR, 2016.
Michalis Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Proceedings of the 12th International Conference on Artificial Intelligence and Statistics, volume 5 of Proceedings of Machine Learning Research, pp. 567 574. PMLR, 2009.
Anthony Tompkins and Fabio Ramos. Fourier feature approximations for periodic kernels in time-scries modelling, Proceedings of the AAAI Conference on Artificial Intelligence, 32(1), 2018.
Jarno Vanhatalo and Aki Vehtari. Modelling local and global phenomena with sparse Gaussian processos. In Proceedings of the 24th Conference Uncertainty in Artificial Intelligence (UAI), pp. 571 - 578. AUAI Press, 2008.
Frida Viset, Rudy Helmons, and Manon Kok. Distributed multi-agent magnetic field norm slam with gaussian processes. In 2023 26th International Conference Information Fusion (FUSTON), pp. 1-8, 2023.
Sander Wahls, Visa Koivumen, H. Vincent Poor, and Michel Verhaegen. Learning multidimensional Fourier series with tensor trains. In 2014 IEEE Global Conference on Signal and Information Processing (GlobalSIP), pp. 394-398. IEEE, 2014.
Andrew Wilson and Hannes Nickisch. Kernel interpolation for scalable structured Gaussian processes (KISS GP). In Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pp. 7775-1784 PMLR, 2015.
Zongmin Wu.. Compactly supported positive definite radial functions. Advances in Compulational Mathemalics, 4(1):283-22) 292-,2299222222 1995.
Mohit Yadav, Daniel Sheldon, and Cameron Musco. Faster kernel interpolation for Gaussian processes, In Proceedings of The 24th International Conference Artificial Intelligence and Statistics, volume 130 of Proceedings of Machine Learning Research, pp. 2971 2979. PMLR, 2021.
15

---

Published in Transactions on Machine Learning Research (MM/2024)
Appendix
The majoriy of the appendix covers proofs of the corollaries following from Theorems 3.1 and 3.4. The rest is dedicated to further details on our empirical experiments as well as visual explanations of the structures that are exploited in the main body of the paper. The appendix is organized as follows. App. proves Theorem for polynomial Basis Functions (BFS). App. B proves Theorem 3.1 for complex exponential BFS. App. C proves Theorem 3.4 for Hilbert space BFs defined on rectangular domain. App. D proves Theorem 3.4 for ordinary Fourier features. App. E contains visual representations of the structures explained in the main body of the paper. Lastly, App. F provides full description of all the included experiments and the data used in them, with additional plots and results.
A Proof of Corollary 3.2 (use of Theorem 3.1 for Polynomial Basis Functions)
Polynomial BESS (as used in for example Chen et al. (2018)) are defined along each dimension as
d 0 (x( (a (4) )
(24)
By selecting Pittsa (x((0) piettaat the product of two component BFS becomes
(z No (x(4)) (xl (zld))iatir-) Piatiar
(25)
which shows that the conditions for Theorem 3.1 is satisfied.
B Proof of Corollary 3.3 (use of Theorem 3.1 for Complex Exponential Basis Functions)
Complex exponential BFS (as used in for example Novikov et al. (2018)) defined along each dimension
of (x((0) exp(injaa!!))
(26)
By selecting (x(())(<(() exp(it 1)) the product of two component BFs becomes (1) (z) = exp(itian!)) (d) exp(injaz( = exp(ir (d) = Patio
( (x. (d) 10
(27)
which shows that the conditions for Theorem 3.1 is satisfied.
C Proof of Corollary 3.5 (use of Theorem 3.4 for Hilbert Space Basis Functions a Rectangular Domain)
Hilbert space BFs on a rectangular domain are defined according to (Solin Särkkä, 2020)
rialaa + La)
1
Tiaaaa La)
(x(x)
II
Tnd
sin
2La
cos
2La
2
(28)
where the indices 1, 1,.... M has one-to-one mapping with all possible combinations of the indices i1, in, given that each index i, = (I....,)) D}.
This corresponds to defining the BFs according to Eq. (6a), with each entry of (bd) defined as
(x) + sin Tia(x + Ld)
Tid(x + Lu)
VLd cos
2L
2La
where ia E f1,... fl,...,m. m). Define a lincar function 0i,, : IR + RR as
(x) Tiala + Tra)
2Ld 2
16

---

Published in Transactions on Machine Learning Research (MM/200 2024)
Hence, the BFs can be written
cosl 01)
2 (x)
Ls
cos (0
The product (xn (zh)))) is then given by
cos(O1) cos(f1)
cos(01) cos(
(1 1 (I:(e) = Ld Σ
cos(O cos(Ons)
cos(O cos(OM))os((MM)
Then, note that
cos(u + v) + cos(u - 2)
cos(u) cos(r)
(29)
2
When we apply this to each entry in the matrix, we get that
Φ (d) (x(())
(x)a))
cos(0; +0,
cos(0; + +0a)
cos((0 01)
cos(0; - -0a))
N
N
2T
cos(OM + 01) cos(f. + +xx)
cos(OM 01) cos(0,r 0m1)
AG400
(g1(-1)
where (d.(1) and Gd,(-1) Hankel and Toeplitz matrix, respectively. Then, since cos(O2, + 0z.) + sin (ja g( ja) g(ka), there exists a function g(ka) that fulfills conditions for Theorem 3.4.
D Proof of Corollary 3.6 (use of Theorem 3.4 for Ordinary Fourier Features)
Ordinary Fourier features for separable kernels are defined slightly differently from Hilbert Space BFS (Hensman et al., 2017), in that they consider both sine and cosine function for each considered frequency. The set of BFs are given by
(x) = (c) (z)
Good (x))
(30)
= =(sin(Ar)
sin (2Ax)
(in))a cos(Az)
cos()Ar)
cos(m&A)].
where A determines the spacing of the Fourier features in the frequency domain. The precision matrix therefore be expressed
T
S
T $
sin
sin
cos
stn
(31)
+ ain
cos
cas
and we can apply Thcorem 3.4 directly entrics + o sin block Mankel-Tooplie strueture.
and + diin T
Poss prove that each of these have a
For the matrix + I Pain the product Φann (x( Joiin (x!(d)) ( can be expanded as
(dhen (p)mm))) Jo. (an(iAn(())sinnAAri = (G(O)--))) + (c(d).())
(32)
where
(G)) = = + (d)
(33)
and
(G(d)(-1)) cos( (d( - (r(
(34)
The entry row and column of (l)) can therefore be defined by the function g(i + cos(i j)Aa(") and the entry at row and column j of (-).(-1) is given by -g(i j), satisfying the requirements for Theorem 3.4.
+
17

---

Published in Transactions on Machine Learning Research (MM/ 2024)
For the matrix o Joe 4 sin; the product Pens (un (4) D ) can be expanded as
(s D (x) ( cos(iior sin(jAx(() - (G + (G
(35)
where
(G(d).))) 1is = sin(i + ),
(36)
and
(G())-)) Jis sin(iAx(d)) jA
(37)
The entry row i and column j of ()d).()) can therefore be defined by the function g + - cos( i) and the entry at row i and column of G(d)(-1) is given by -g(i - j), satisfying the requirements for Theorem 3.4.
+
For the matrix + cos 4 coss the product Deos (xn ) (xth)) can be expanded as
[rees )
(20 (2! T
cos()oo (G(d)--)) (G(O Jags
(38)
where
(G()) foy= cos + )Axx(
(39)
and
(G(4)(-1) Ja- cos(iAz(()) - jAz(d)).
(40)
The entry at row i and column j of (d))(() can therefore be defined by the function g(i+j) = cos((it))Arl
and the entry at row i and column j of (-d),(11) is given by g(i-j). An important notion for this matrix which makes it different from Pain & si and & Pcos sin is that this does not exactly satisfy the criteria for Theorem 3.4. The difference is that for the criteria be exactly satisfied, entry (G(d).(-1)) Ja.j would have be equal
-g(i- j) rather than g(i -1)) j). However, by applying the proof of Theorem 3.4, but now noticing that the entries of C $ T cas can be expressed dement-wise
2
D
c,== == glia + ja),
(41)
1-1n=ld=11
which allows us to use the tensor 7 as defined in Eq. (22) to express each entry of C according to
2
Oij = Σ No . in
ep)) (P))
(42)
E Overview of Hankel, Toeplitz, and D-level Block Hankel-Toeplitz Matrix Structures
An overview of Hankel, Tocplitz and D-level Block Hankel Tonkel-Topplitz matrix structures are given in Table A1.
F
Experiment Details
More specific details of the numerical experiments are provided here with additional visualizations and explanations.
F.1 U.S. Precipitation Data
The precipitation data bwo-dimensional with N 5776 data points first considered in (Vanhatalo & Vehtari, 2008). We perform regression in the lat/lon domain and first center the data (but do not perform scaling) and use simple squared-eexpoentiall kernel. We optimized the MLL using GPJax (Pinder & Dodd, 2022) for both the HGP as well as the standard GP for 100 itcrations using Adam (Kingma & Adam, 2015). The hyperparameters of the kernel and likelihood were initialized as / = 1, 0 10 and a, = 1, where is the lengthscale, o' the kernel variance and oe the noise standard deviation. Purely for visualization purposes, the inputs were projected to a local coordinate system given by CRS 5070 which are used for all of the plots.
The original data. is plotted in Fig. A7. The timing experiments were run using m. = 5. 10...., 65 BFS along each dimension, totaling between M 25 and M 4225 BES.
18

---

Published in Transactions on Machine Learning Research (MM/ (MM42024)
Table A1: An overview of the matrix structure and tensor representation for Hankel, Tocplitz, block Hankel, block Toeplitz and block Hankel matrices. The illustrations are examples of matrices with the property described in each row. The illustrations contain one square for each matrix entry, where the color of the square corresponds to the value.
STRUCTURE
DEFINITION
VISUAL
DOMAIN
722
Ym
72
+m+1
HANKEL
H
Y R
7
72m-1
Ym
"/1
TOEPLITZ
T
7m+11
73
7/2
7 R
Y2m-1
7m+1
H (D-1)
H (D-1)
Hmp
(D-1)
H D-1)
H D-1)
H (D11)
D-LEVEL BLOCK HAN-
HD))
mpti
E Kix..XX
KEL
H
H (D-1)
H (ー-1)
mp+1
2mp-1
T(P-1)
T (-11)
T (D11)
T(D-1)
T
D-1)
T D-1))
D-LEVEL TOEPLITZ
BLOCK (T()
7 KiX...XK
T (
T(D-1
T
(D-1)
mp
G (n-1)
D-11)
(D-1)
D-LEVEL BLOCK HANKEL HANKEL-TOUPLIVZ
G
(-11)
D-1)
(D-1)
G(D)
mp+1
E R X
LEVEL IS HANKEL
Gmp
(D-1)
(D-1)
mp#l
2mp
(D-1)
G
G DDD11
Gmp
D-LEVEL BLOCK HANKEI-TOEPLLTZ,
GDD1)
D-1
D-1)
G
Y RKIX..KKK
LEVEL D is TOEPLITZ
G2mp-1
Gmptl
Gmp
(D-1)
19

---

Published in Transactions on Machine Learning Research (MM/2024)
3000 [m]]]
Figure A7: Raw data from precipitation data set. Each data point is visualized cross with color indicating the precipitation. The data set contains mostly low frequency content with some high frequency content apparent in the west coast as well as the south eastern parts.
F.2 Magnetic Field Mapping
The magnetic field data has N ss 1.4 million data points and is collected with an underwater vehicle outside the Norwegian coast. The data used were collected by MARMINE/NINU research cruise funded by the Research Council of Norway (Norges Forskningsrád, NFR) Project No. 247626/O34 and associated industrial partners. Occan Floor Geophysics provided magnetometer that was used for magnetic data acquisition and pre-processed the magnetic data. The data was later split into training set and test set, at roughly 50% split. The nature of the data split is visualized in Fig. A8b. However, in practice, we selected the width of the test squares and the training squares smaller than the one displayed in the illustration and they are merely that big for visualization purposes. The illustration displays squares that are 0.01 latitudinal degrees wide and 0.03 longitudinal degrees tall, corresponding to approximately km Cartesian coordinates this area. The split we actually used was squares which were 0.001 latitudinal degrees wide and 0.003 longitudinal degrees tall, corresponding to approximately 100 m in both directions in Cartesian coordinates in that area. GP regression with squared exponential kernel is able to extrapolate for approximately one lengthscale, but will not necessarily give a very informative prediction one or two lengthscales away from the nearest. measurement. Although we do not know the lengthscale that would optimize the likelihood of the data before using training data to find it, we see from zoomed-in version of Fig. 5b approximately how fast the magnetic field is changing across the spatial dimension and use this to make a reasonable guess at the distance we expect a well-tuned to be able to extrapolate the learned magnetic field. We then project the data into local coordinate system using WGS84 and perform regression in Cartesian coordinates. We center and standardize the data with the mean and standard deviation of the training data. A squared-exponentee kernel was used and we optimize the MLL in GPJax (Pinder & Dodd, 2022) using Adam (Kingma & Adam, 2015) for 100 iterations to find hyperparameters. The hyperparameters were initialized as I 1= 200 m, oe. and o? = 1. The resulting hyperparameters were IsE = 190, of = 0.0675, oge 7.15.
20

---

Published in Transactions on Machine Learning Research (MM/2024)
Training data
Test data
73.46
T3.48
Longitude (*)
Longitnde
73.42
7.45
7.5
7:55
7.6
7.65
7.45
7.5
7:35
7.0
7.65
Latitude (?)
Latitude (?)
(a) Raw magnetic field measurements for the underwater magnetic field data. The plotted data is subsampled to th of the data for visualization purposes.
(b) Data divided into training and test set. Determin- istically split to ensure ensure that the lengthscale is captured properly in the training data.
Figure AS: Magnetic field training data and test data. The data is roughly 50/50 split between training and test set. Both training and test data are normalized only by the mean and standard deviation the training data.
21


Reviews:
[{"review_id": "w6ygz03OVD", "reviewer": "Reviewer_TQev", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The results shown are of interest as the scheme proposed does not introduce additional approximations but is simply a more efficient way to perform the computation. The main contribution is in showing the sufficient conditions for the choice of the basis functions to exploit the prosed speed-up.\n\nThe weaker element of the proposed work are:\n1. a clear illustration of the wider applicability and limitations of the results: \n- the approach works for polynomial, complex exponential and sinusoidal basis functions with frequencies on a grid, for Hilbert–space Gaussian Process defined on a rectangular domain \n- while the efficiency of other quadrature-like methods that are not defined on regular grids cannot be improved, e.g. approaches based on random frequency sampling or based on inducing points\n- it would be useful for the reader to present a list of the use cases and applications that are most commonly addressed with the type of kernels that can be improved and of the cases that need to be modelled with kernels that are not subject to the speed up improvements\n\n2.  It is not immediately clear how and if these results can be generalised, i.e. how one could approach the problem of approximately (but automatically) decomposing a matrix in a composition of Hankel-Toeplitz matrices \n\n3. It would be useful to empirically show how the reduction of the computational and memory costs can allow the solution of problems of practical interest that were previously unfeasible (on a given computational budget), e.g. showing that a class of climate/weather predictions is now achievable on common hardware in reasonable timeframes.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors propose a way to speed up Gaussian Process inference when the precision matrix exhibits a special structure. Specifically they show that if the precision matrix can be decomposed in a combination of Hankel-Toeplitz matrices then they need to compute a number of parameters that scales linearly rather than quadratically with the number of basis functions used for the approximation of the precision matrix. This results in practice in speed ups of one or more orders of magnitude in the computation of the precision matrix, which they confirm in empirical studies.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- a list of the use cases and applications that are most commonly addressed with the type of kernels that can be improved and of the cases that need to be modelled with kernels that are not subject to the speed up improvements\n- a discussion on how these results could be generalised, i.e. how one could approach the problem of approximately (but automatically) decomposing a matrix in a composition of Hankel-Toeplitz matrices \n- a comparison on the relative quality vs speed-up between the proposed approach and alternative approaches (e.g. using inducing point approximations)", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "bZE2BYTrnQ", "reviewer": "Reviewer_2CuJ", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "*Strengths*\nThe speed-up presented by the paper can be beneficial in a number of settings that satisfy the hypothesis of the main theorems in the publication. The authors further presented some numerical experiments on real-world datasets, showing the benefit of their structured approach.\n\n*Weaknesses*\nThe paper focuses on a specific, yet important, structure of the basis function approximation scheme. \nMoreover, the writing of the paper can benefit from careful re-reading, as some parts are not, in my opinion, very clear. \nIn some aspects regarding the writing, the manuscript is unfortunately sloppy: I added below a list of changes, some of which I think are necessary for a publication.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The submitted contribution focuses on the speed of the computation of kernel precision matrices by the exploitation of their possible Hankel-Toeplitz structure. In particular, the analysis is developed in the context of the Hilbert-space Gaussian Process approach, in which a covariance matrix is estimated from a dataset of $N$ data points by projecting a Gaussian process on a suitable set of $M$ basis functions. The standard technique requires $\\mathcal O(NM^2)$ steps for the estimation of the precision matrix: as the needed number of basis functions might be large, the computational cost might be demanding. The authors show that a block Hankel-Toeplitz structure assumption on the precision matrix allows for a speed-up of the analysis, reducing the complexity to $\\mathcal O(NM)$.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- In the introduction, it is claimed that the paper ```We show that each increment of the HGP precision matrix can be split in a multilevel block-Hankel and a multilevel block-Toeplitz matrix```: I am not sure if this sentence is accurate. From what I understand, the specific Hankel-Toeplitz structure considered in the paper is not universal in the HGP approach. If this is the case, I would suggest to modify the sentence. \n- In Eq. 5, is $\\boldsymbol \\phi^{(d)}(\\boldsymbol x)$ to be replaced with $\\boldsymbol \\phi^{(d)}(x^{(d)})$?\n- There are some notation imprecisions after Eq. 9 (e.g., $\\boldsymbol H_k^{D-1}$ should be $\\boldsymbol H_k^{(D-1)}$). Also, the notation in the labels in Fig. 3 does not look consistent with the text: please check.\n- Section 3.3 seems to deal with block Hankel matrices only: if the authors meant to refer to block Hankel-Toeplitz matrices (as the context, the title, and the appearance of $\\boldsymbol G^{(d)}$ suggest), they should specify it. This is actually quite a crucial point to clarify.\n- In Eq. 19 and in Appendix B, the notation $\\pi_{j_d}$ appears: do the authors mean $\\pi j_d$? Please fix.\n- In Eq. 20 the index $n$ in  the summand disappear.\n- In Section 4, at the end of the first paragraph the authors refer to the computation of $2M$ and $3M$ components: were the experiments performed using $D=1$? (From Appendix F.1, it does not seem so). Please clarify.\n- The writing of Section 5 might be improved (just to mention one aspect, the expression ```structure we exploit``` appears five times in two paragraphs...)\n- In the main text *Corollaries* are introduced for specific bases, but in the Appendix, titles of corresponding proofs refer to the main Theorems instead of the corresponding Corollaries. Please fix.\n- What does the subscript $t$ in Eq. 28 refer to? Why does the dependence on $d$ disappear from $L$ in the subsequent formulas?\n- Appendix E is **empty**! If it was supposed to contain Table A1 (as it seems) please add one line to reference the content.", "broader_impact_concerns": "There are no concerns regarding this publication, which is of theoretical/algorithmic nature.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "vMMN7agvJl", "reviewer": "Reviewer_Kfc8", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths ###\n\n- The paper is well-written. The mathematical reasoning is presented clearly.\n- Experiments are done on real world datasets that confirm the theoretical claims.\n- The proposed algorithm appears to be mathematically sound and applicable to a variety of kernels.\n\n### Weaknesses ###\nThe exposition feels incomplete to me.\n- The proposed method decreases the computational complexity of calculating the precision matrix. How does it affect the entire posterior mean/covariance computation (if any)?\n- It would be helpful to readers to include more detailed pseudocode, if only for one type of kernel.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper provides a method to reduce the computational complexity of calculating the precision matrix in Hilbert-space Gaussian Processes without introducing approximations. For many kernels, the precision matrix is the sum of Hankel-Toeplitz matrices, that have $O(M)$ unique entries. Therefore, the computational complexity can be reduced from $O(NM^2)$ to $O(NM)$, where $N$ is the number of data points and $M$ is the number of basis functions. Experiments are done on synthetic and real datasets confirming their theoretical claims.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Critical:\n- Improving the exposition as suggested in the previous question.", "broader_impact_concerns": "The authors state in an impact statement that since their work proposes a new numerical method, it has many potential uses and there is no need to specifically highlight certain societal impacts. I agree; since the proposed method is general-purpose, it is difficult to assign either positive or negative impacts.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This manuscript considers the problem of Gaussian process (GP) inference. A naive implementation of GP inference requires computational effort that grows cubically with the number of training observations, which limits scaling without further effort. This manuscript considers one approach to speed up the inferential process, namely the Hilbert-space Gaussian process (HGP), which operates via projection onto a set of basis functions. The authors show how to exploit convenient structure in HGP precision matrices to accelerate the computation required for the HGP approximation -- without any further approximation -- from scaling quadratically to linearly in the number of basis functions. 

The reviewers were universally positive in their initial assessment of this work, although the initial reviews did identify a number of minor issues. These issues were adequately addressed during the author-reviewer discussion period, and the manuscript was strengthened as a result. Ultimately the reviewers agreed that the paper satisfies both TMLR acceptance criteria and were unified in their recommendation of acceptance. I agree with their assessment.

I encourage the authors to reflect on the reviewer feedback and author-reviewer discussion in preparing the camera-ready version of their manuscript.
Decision: Accept (poster)