Giga-scale Kernel Matrix-Vector Multiplication on GPU
Authors: Robert Hu
Siu Lun Chau
Dino Sejdinovic
Joan Alexis Glaunès
Venue: NeurIPS 2022 Accept (2022)
Keywords: Kernel
GPU
RKHS
kernel matrix vector multiplication
interpolation
large scale
giga scale
10^9 points
big data

Abstract: Kernel matrix-vector multiplication (KMVM) is a foundational operation in machine learning and scientific computing. However, as KMVM tends to scale quadratically in both memory and time, applications are often limited by these computational constraints. In this paper, we propose a novel approximation procedure coined \textit{Faster-Fast and Free Memory Method} ($\text{F}^3$M) to address these scaling issues of KMVM for tall~($10^8\sim 10^9$) and skinny~($D\leq7$) data. Extensive experiments demonstrate that $\text{F}^3$M has empirical \emph{linear time and memory} complexity with a relative error of order $10^{-3}$ and can compute a full KMVM for a billion points \emph{in under a minute} on a high-end GPU, leading to a significant speed-up in comparison to existing CPU methods. We demonstrate the utility of our procedure by applying it as a drop-in for the state-of-the-art GPU-based linear solver FALKON, \emph{improving speed 1.5-5.5 times} at the cost of $<1\%$ drop in accuracy. We further demonstrate competitive results on \emph{Gaussian Process regression} coupled with significant speedups on a variety of real-world datasets.

Giga-scale Kernel Matrix-
Matrix-Vector Vector Multiplication on GPU
Robert Hu Amazon
Siu Lun Chau
Department of Statistics
robyhu@amazon.cm co. uk
University of Oxford
siu. chau@stats ox ac uk
Dino Sejdinovic
School Computer and Mathematical Sciences
University of Adelaide
dino . sejdinovic@adeeaidde edu. au
Joan Alexis Glaunès MAP5 Université Paris Descartes
alexis. glaunestmi....mm
Abstract
Kernel matrix-vector multiplication (KMVM) is a foundational operation in ma- chine learning and scientific computing. However, as KMVM tends scale quadratically in both memory and time, applications are often limited by these com- putational constraints. In this paper, we propose a novel approximation procedure coined Faster-Fan Fast and Free Memory Method M)) to address these scaling issues of KMVM for tall (108 10") and skinny (D S 7) data. Extensive experiments demonstrate that F'M has empirical linear time and memory complexity with a relative error of order 10-3 and can compute full KMVM for billion points under minute on high-end GPU, leading to significant speed-up in comparison to existing CPU methods. We demonstrate the utility of our procedure by applying it as drop-in for the mate-of-the-ar GPU-based linear solver FALKON, improving speed 1.5-5.5 times at the cost < 1% drop in accuracy. We further demon- strate competitive results Gaussian Process regression coupled with significant speedups on a variety of real-world datasets.
Introduction
Kernel matrix-vector multiplication (KMVM) is one of the most important operations needed in scientific computing with core applications in diffeomorphic registration. geometric learning [1]],
[3I], numerical analysis (28) fluid dynamics [6], and machine learning (27]. For dataset of size n, KMVM using direct computation has complexity and memory footprint (nz?), both unfeasible for modern large scale applications where 10 is becoming increasingly common. Pioneering contributions presented in the Fast Multipole Method (FMM) 10 amend the complexity of these tolerance, with varying reductions in memory
problems to On log (e- (c)),,, where the chosen footprint for data restricted to dimension D 2. Subsequent developments
17) mainly focused
on extending approximations for a broader set of kernels for fixed dimensionality D tailored for problems in physics with narrow data such as electrostatics, stellar dynamics, Stokes flow, and acoustic problems, amongst others.
Work mainly done while the authors were with the Department of Statistics, University of Oxford.
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

persieliaatio oplimisation We imprown the geometric partieeiing reduce memory constraints, leading lo improvement of 5-7
al (2018) introduced kemel apnosdoo interpolation
We further introduce edaptive field and smooth fied approsimation speed from time
(198))988)
2FM()019)
F
Figure 1: A brief summary of the evolution of the FMM family. Our contributions
In this paper, we introduce Faster-Fast and Free Memory Method (F$M), a novel algorithm built upon the FFM [3] framework to perform KMVM on a GPU for tall and skinny (D S 7) data of order 10 in under a minute with user-specified error tolerance, providing between 2 - 8500 times speed-up over existing methods. It should be noted that the constraints on D and n are not inherent formal constraints, but a reflection of practical limits with typical current computational resources.
Notations.
We use capital and lower case bold letters to represent matrices and vectors, respectively. In this paper, will work with matrices X E R Y E RR'' and vector b RRv, For kernel k the goal for KMVM is to compute v := K b, where K == k(X, Y) = (k(x; N,) Ml and Y, denote the ith 5th row of X, Y respectively.
1=1
2 Motivation and Related Work
Kernel methods are often limited by their (n²²) memory footprint and computational complexity for KMVM. These constraints make scaling beyond 10 challenging. Many recent developments have been made to improve both of these constraints, ranging from hardware acceleration using GPUs in KeOps [_] to various approximation techniques proposed in [36 33 3 focus our attention on kernel independent KMVM methods.
9.
In this work, we
KeOps. Charlier et al. 11 proposes map-reduce scheme compute kernels using exactly O(n) memory and Onn²) complexity on GPU. This is achieved by computing the full KMVM product on-the-fly by summing o = E k(x,,,y) yy)b, directly, without ever storing the kernel matrix K E R=*=" explicitly. Extensive experiments show that this method is practical when n S 10 as the GPU hardware acceleration allows the KMVM product to be computed in less than a second conventional GPU. Moreover, the method places no constraint on the number of features D can be applied to, making it favourable for KMVM medium size datasets. In application contexts, KeOps is currently adopted into conjugate gradient solver FALKON 24 26
as part of the default pipeline.
Table I: Comparison between methods( indicate ours).
The Fast and Free Memory Method (FFM). While KeOps can theoretically scale to a bil-
Method
FMM
KeOps
FFM
3M"
F'M
Kemel Independent Linear Time
lion points, becomes practically infeasible as the O(2²) complexity would imply a computa- tional time of 10" seconds, or roughly days. To overcome this billion points barrier, Aussal et al. 3 deploys geometric space partitioning scheme, and proposed the Fast and Free Mem- ory Method (FFM), a KMVM approach that extends the FMM [0] family of algorithms. In contrast to traditional FMM methods, which require specific series expansion of the kernel, FFM deploys Lagrange interpolations to approximate them instead. This allows FFM to be applied to almost any conventional kernel and further enables the user to trade off accuracy with computational efficiencies by controlling the order of the approximating polynomial [19]. Compared to KeOps, FFM demon- strates both linear memory and time complexity in experiments and scales to compute billion-pointe KMVM on a smaller CPU cluster under hours, outscaling the GPU implementation of FMM (21). While 4 hours is significant improvement compared to 11 days from KeOPS, it still renders many machine learning techniques infeasible. Further, as recursive partitioning of the data space scales poorly with D [5]. both FMM and FFM can only be applied to D < 3 data, price to pay for the speed-up of KMVM operations when n = 10 Furthermore, we show in our experiments that a direct FFM port to GPU gives unstable results for n 10", D for non-trivial data simulations (bottom row in Appendix TO).
Linear Memory
Restriction in D GPU
S3
3
S7
ST
Scales under I hour
10"
under 1 minute!
KMVM applied to 3D data for 108 on Nvidia V100 GPU.

---

Our contribution. To surpass the billion point barrier while maintaining high-speed and stable computation, we propose F2OM and our main algorithm F the first pair of KMVM algorithms that can reliably scale to 10 on skinny data using single GPU. We build F on top of FFM by introducing non-trivial GPU parallelisation and low-level optimisations. We further stabilize and improve the original geometric partitioning scheme in FFM to significantly reduce memory constraints, leading to a relaxation of dimensionality constraints from 3 to 7. At last, we introduce an adaptive far-field and smooth field approximation scheme for kernel interpolation, resulting in our main algorithm M which runs 2.0 33.3 times quicker and more stable than a direct port of FFM on GPU. See Fig. and Table T for overview and comparisons of the methods. We summarise contribution as follows:
1.
We propose Faster-Fast and Free Memory Method (F) a KMVM algorithm building on top of FFM by applying multiple low-level enhancements, GPU parallelisation, and algorithmic computational and memory enhancements, allowing for KMVM operations on S 10 data in under a minute. Codebase is released here 14
2. We characterize theoretical time and memory complexity of F
3. We run extensive KMVM experiments of F'M on variety of tall and skinny data with S 10", demonstrating empirical linear time and memory scaling, and achieving speedups between 2-8500 times when compared to FFM (GPU and CPU) and KeOps.
We run a practical application of F as drop-in replacement for KeOps in conjugate gradient solver FALKON (24 26 for kernel ridge regression and classification (KRR) giga-scale data, obtaining a solution 3.4 times faster with <1% drop in accuracy. We further demonstrate competitive results on Gaussian process regression against KISS-GP [35] SVGP [18] and SVGR (30 with significant speed-ups.
3 Background
The FFM method considers KMVM for kernel k evaluated on data matrices X, Y and are weights associated with Y. The KMVM is expressed as v == k) For example, b could be the weights in a KRR or the strength of electronic charges. As n2 and ny are taken to be very large, a full computation is unfeasible. In this section, we il- lustrate and detail the main steps of FFM, before presenting our improvements Section For illustration purposes, we first consider simple 2D
b.
KMVM. Our goal calculate k(X, Y): b for X, Y in 2 The intuition behind FFM is to reduce the com- plexity of calculating the full KMVM by partitioning X and Y such that certain calculations can be approximated fast manner, based on the pairwise distances between partitions.
Enclosing and partitioning the data. The first step is to partition the data. To begin, we find large enough box that can just enclose or Y. The edge length of this box is calculated as
Figure 2: Enclosing X and Y within box, "x" marks the center of the box. Numbers der the the boxes denotes edge length. In the right plot, we have enclosed the blue points with the largest box.
E max max(<1/0)
x xmau) max
(d)
where 2max, 2 denotes the largest value and the smallest value along the d-dimension in X and similarly for Y. Figure jillustaates this enclosing procedure.
Defining near and far-field
FMM, an octree 23 is applied to recursively partition data into
smaller boxes B. C X, BY Y, with g de- noting box indices. Here each box corresponds to a subset of rows in the data matrix. Let us also denote b, as the partition bj's grouped with the same indices as By To calculate the KMVM between two boxes B.X, BY with the grouped vector by for each E B com-
Figure 3: Recursive partitioning of X and Y for 2D data. Far-field interactions are colored green while near- 'feldd interactions are colored orange and denotes the euclidean distance between the centers.

---

pute
M(x,,,))
(1)
with p0 = [0p ... o. Now the target can be computed as v = ((v v PTT, where v" Σ v P. Q denote the total number of boxes and a (-) a permutation such that of appear in the same order as x; appears in X. Figure shows how boxes are recursively partitioned.
Far and near-field interactions FFM relies on a divide-and strategy to effectively com- pute KMVM product; data is partitioned into boxes and then separated into far-field and near-field interactions, where near-field interactions are computed exactly and far-field interactions are approxi- mated using Lagrange interpolation for speed, explained in the paragraph below. The partitioning procedure in FFM is recursive, where the recursion depth tree_ _ree_depth controls the size of the edge 1 2ixoo,daput of the box. An interaction is defined to be in the far-field if the distance between the two center points of the boxes exceeds 21, i.e. IsFarField == |xcenter Ycemer| 2 2/, While for each box will decrease with the number divisions, this rule ensures fixed minimal distance for a given depth for far-field interactions. Figure 3 illlusraaee how far(green) and pear(orange)--ier interactions arise between X and Y when tree_depth increases.
Lagrange interpolation We review Lagrange interpolation used for far-field approximations in FFM. Given function ((x) [-1,1] IR r-+1 unique points 8; [-1,1], = 0. ......77 there exists a unique polynomial pr (x) of degree I r that interpolates at pr (8,) = f(:)). The Lagrange polynomial is given by (t) = Ei (t), where Ct(t) = i. (t-b,)
O....,
12
well the points s; interpolate through. The choice of 8; is especially important in minimizing large oscillations around the edges of the interpolation interval (Runge's phenomenon (12). For this reason, Chebyshev nodes of the second kind are used = cos 0,, where 0, 1M i= O....,
We free chose the degree
Interpolating y(x,y) By noticing that k(x,y) is bivariate function, interpolation twice, thus interpolating y(x,y) as k(x,y) ss Σ Li(x) E k sy) (y). Here r' denotes the number of the interpolation nodes and s s E R.D denotes the grid of interpolation nodes for B and BY Note that since x E R.D, we take (,(x) =
apply Lagrange
Il2=1 I1 (200)) (2(0))a(())
IT: (1.0 (1 0 r. These operations can be vectorized and computed se-
quentially on-the-fly with linear memory footprint v so L (K (Ly b)), which is done by first
V1
computing V1, then v2 and lastly Here Lx denotes matrix with entries (xyx), where indexes the rows and j the columns, with Ly following the same definition for y instead. A far-field
KMVM between two boxes Vp.g = kBBx By b, is then approximated by using double Lagrange interpolation according Figure
1. Interpolate over Y' and sum over
Translate over interpolation nodes af
3. Interpolate over and sum over
6 Meam
IC 1 -aj Tyaa
K
Miy' symll R
Were (0t -a) IElaso (e Theee (o -as Tifcom (n -a))
Guid Jost
SUNDIAA
Y1 ERY
E
Normalize data
Interpolate and sum
Normalize data
Interpolate
Approsimution complise!
Figure 4: Here we approximate a far-field interaction between two boxes. In 1. we first normalize the data between 1-1, 1] and conduct 2D interpolation of k(- y) while summing over b. In 2. the interpolation for k(xx,) -) will also be normalized and hence need to translate the distance between the boxes by calculating K Vi. Lastly in 3. interpolate k(x,,) while we V2-

---

4 Faster-FFM (F M and M*))
To fully leverage the port of FFM to GPU, we enhance FFM with novel approximation procedures for improved complexity and memory optimizations scale to 10 We coin this improved version Faster-FFM (F M). The capabilities of F against previous methods are summarized in Table |
4.1
CPU to GPU optimizations
In FFM, every computation is serial and on CPU. When moving to GPU, we have parallelized all major computations. These parallelizations are non-trivial and require low-level algorithmic optimizations, with challenges such as:
tor-errbbca alignment- A major challenge in the implementation of both the parallel far-field and near-field computations was correctly aligning thread blocks to boxes. This aligning requirement imposed non-trive anniiiiiial boundary conditions on data indexing when using shared memory. To minimize memory usage of box and block indicators for our implementation, we represented the box belonging of each point as index intervals (i.e. box consists of points with E [1...., 500] and box 2 with i E [501,... 1337], etc.) and modulo arithmetic to infer the block belonging. This clearly requires that the points are sorted or grouped according to their box belonging. However, as we detail in the next paragraph, arranging the points could not be done straightforwardly with native sorting methods. We further illustrate how parallelization is done for calculating near-field interactions in Appendix
No native sorting methods We found that LibTorch 25 sorting methods often led to out-of- memory (OOM) due to allocation of large long-type vectors on GPU. When n 100, this implies
allocating 8GB of memory, 25% of the 32GB
card used, making it a necessity to avoid native sorting methods.
Input: A sorted matrix of interactions Io
Step
Step 1
Stop 2
In-place grouping data on boxes Due to in- feasible LibTorch sorting methods, we addition- ally had to design algorithm that finds a per- mutation that would group into its correspond- ing boxes in linear time and memory. We used count and increment- baased strategy that would: 1. Count the number of points in each box dur- ing the assignment operation (O(n) (1)) and store cumulative
Expand
Reppes right indox
Repeett left Indices.
Expand
right
new arckoes in
dbvie box
Indicos nto
order
indicos
within
new
according
to right
0
D
3
2
ndices, by
Index
shiciing
1
0
inckaes
1
in order
Figure 5: Here we assume D = 1, hence we only divide each box 2n times each time.
the count in vector ç. Then sum over ç, starting from O.
2. Initialize a long permutation vector Us- ing the counting vector ç, we would re-run the assignment operation and arrange a point with index i as following π rlatomii a/ (copd(""""[iiidia points in each box. We specifically have to use the function atomi cAdd to increment the count for each box in parallelized GPU environments avoid thread locks.
: 1) + [ioo__]nexx]]]__ where g is a running count of
We refer to the box_ divisiioncum_ppee hash and box_division_assie hash function in n_ treer. cu for exact details.
Ensuring interactions are sorted sary sorting,
To that the matrix
avoid any unneces- containing
ensure
interactions
is always sorted by recursively dividing old in- teractions. We illustrate the procedure in Fig- 3 We refer to the get_ new interactions function in trreee. cu for the exact implemen- tation.
n 10
10"d
003
Time
203
However, we found that these optimizations and porting alone were not enough to scale to = 10 on 3D datasets, as Figure 6 demonstrates. FFM doesn't 't remove empty boxes or handle boxes with few points in them and keeps expo- nentially creating new empty boxes and inter- actions, thus leading to out-of-memory (OOM) errors non-uniform data (see Appendix 10
1231 Chaterdd F151
Non
Uee that
UM1 Dhakent HMM
Serr
U.S
Datasets
Datascta
Figure 6: KMVM times on several datasets between FFM(GPU) vs 10
F'M. U&N was only run up to
5

---

4.2 Scaling to n = 10 on GPU (F2.0M)
In this section, we detail the memory enhancements that allow F'M to consistently scale to n 10
Removing empty boxes with hash list indexing To ensure linear memory on GPU, only keep a reindexing vector a of size By (resp. ny) in memory during the computation of the algorithm in addition a list of interactions and box centers. This reindexing vector rearranges the data points so they appear in the order of the box they belong to. We optimize both the computation and the memory footprint of these objects by avoiding recursive formulas and hash lists.
Naively, points can be assigned to boxes by direct comparison to all existing box centers. As the number of centers grows exponentially with depth tree _eethh this method quickly becomes pathological. To amend this, we propose a linear complexity formula to retrieve the box index (3:
D
point x E X C R. belongs to B, = 2 (gree_dephh(d-1] [ 12tree_dppth Dd Odd where ag denotes
E
E(0,1), Denotes left or right center of box edge
Summing over dimensions
the minimum value of X in dimension D and Fd is the value of in dimension D. To prevent the number of boxes from growing exponentially, we remove empty boxes with each division. To assign points to the corresponding boxes, group points {x: to their respective ordering using the hash list in (n) time in contrast to On 2D-tree_depppee by direct computation.
hash list to store 3; and the order i. We can then
Handling boxes with few points with small field In cases when the number of points in each box vary greatly, we separately consider the interactions where the number of points in boxes is small. Hence, we say that there is a small field interaction between boxes B. B if both have a small number of points, i.e. if
p
1B I p. for some threshold number p. To minimize the computations needed, p can be set to rx + ry. This intuitively allows F to directly compute interactions that too small to benefit from interpolation savings (i.e. | + Bxx < rx + ry), thus limiting memory usage by stopping partitions from dividing further than necessary. In higher dimensions where the division rate is faster, P can be set to a higher value to limit memory usage the expense of more direct computations which are slower.
Figure 7: Full grid vs sparse grid.
Sparse grids As the number of Lagrange polynomials increases exponentially with dimension, we implement sparse grids (29] to allow for finer selection of interpolation nodes. With sparse grids, the number of nodes needed grows slower [20), thus saving memory. We give an example of a sparse
grid versus full grid in 2D in Figure 7
4.3 Speeding up 2 M (F2- F
Smoothness criteria FFM speeds up its computations with minimal loss
in accuracy by selectively interpolating To improve speed, we introduce the smoothness criterion to widen
interactions that are far
apart.
the selection of interactions that can interpolated with minimal loss in accuracy. For a Gaussian |x-y|"
15 27
213
1004
Kernel y(x,y) = exp 27
Rel Enur
Rd Buns
Eiee
with lengthscale 7, the smoothness criteria is defined as is is_smooth
TaS00 Sqaar ditener
-3
Sqaae dotem
"ax Ex,ED Ed, 11
(a) Dimension d
(b) Dimension d' d= 5
2,2
Figure 8: Plotting relative and absolute error against squared distance between boxes. "0" nodes mean we use the zero vector as an approximation to the KMVM for gaussian kemel. We observe that its not beneficial to interpolate at all when the square distance exceeds 5 and that is sufficient only use 3" nodes when square distance is < 0.01.
Σ Ey,eng
Ed Iy (4)
< n
2-2
between an adjacent interaction of boxes B,,, By The quantity com-
puted can be understood as "Effective of the Gaussian kernel. We justify the
Variance" (EV), as considers total variation in the exponent
smoothness criteria with the following proposition.

---

S " < for all x, y. When interpolating k(x,y) exp (-d(x,)) using bivariate Lagrange interpo- lation Le (x(x.) y)= ;= L K (Ly b) with degree = 2p, for any p E N50 there exist nodes s s for (x y) such that the pointwise interpolation error is bounded by (.
Proposition 1. Consider x,y E X' C RRd such that y)x,,)
xx xx 2T 2 E
See Appendix D for proof. Hence for small < 1, specified C, y). To avoid calculating the sample variance during computations which costs (n)). exploit that data is partitioned into hypercubes with known edge & and take the upper bound of
that the error becomes small for well
the variance each cube as & along a dimension. A proof for this bound is provided in Appendix B
Adjacent interactions are then classified smooth when Ea D E² (EV of B ) 4-27 DE 7 S 7) which only costs (1) to compute.
g (EV of B
4-27
Adaptive far-field approximation To further improve speed introduce adaptive rule to select the number of interpolation nodes used when calculating far-field interactions. Error bounds for multidimensional Lagrange interpolation have been proposed in 22 however, these bounds cannot be directly used to create an adaptive interpolation rule. We thus simulate KMVM errors for kXX, Y) b where X, Y are uniformly distributed and b is normally distributed. We fix distance between X and Y and vary the squared of this distance between boxes against nodes in Figure g We use Gaussian Kernel with 7 V2'
Based on Figure g we use the following rule for selecting the number of interpolation nodes for far-field interactions
(ii(r)3))
if E
S
0.01
29
Pfar (r)
if 0.01
S 5
2
27
if 2 2
24
where r is the number of nodes chosen to interpolate with in the general case.
Barycentric lagrange interpolation We slightly improve the complexity further by implementing barycentric Lagrange interpolation 7 evaluated at the Chebyshev nodes of the second kind. As this is well-known technique, we refer to the appendix for more details. It should be noted that the above methods can straightforwardly be extended to any translation- -inttiiannn kernel by recalculating the Taylor expansion for smoothness criteria and rerunning the simulation for adaptive far-field approximation.
4.4 Complexity
The time complexity of FFM is log (n)) and we use a similar derivation strategy for F'M to obtain a complexity that is dependent on the effective variance limit 7) (chosen parameter) and the box width & (data). We first present two propositions needed to derive the complexity of F'M.
Proposition 2. A far-field interaction between two boxes containing n, and ny points respectively has time complexity Om)) where n max (Mx, Nly).
Proposition 3. Given n. data points in dimension D, the maximum number of divisions
Tree dswiasons is given by
Tree divisions = log (n).
With the above results, the complexity of FFM is taken as the maximum number of divisions multiplied by the complexity of fir-field interactions at each division which yields O(n log (n)). We remark that near-field interactions between boxes containing only I data point have linear time complexity, hence the results hold.
Theorem Given a KMVM with edge (dependent on data X, Y,)), lengthscale y, effective variance
) which limit n, data points and data dimension D. F has time complexity O(n loga (2.4.1)
D.E2
can be taken as O(n logg
( where C P.
Memory footprint As our implementation uses the same partitioning strategy as FFM, the the- oretical memory complexity remains (n) for F'M (see [3] for proof). However, this does not accurately reflect the memory footprint of the actual implementations, whose memory mostly depends
7

---

on the number of interactions stored. We summarize these memory footprints for FFM and F in Theorem Theorem[below.
Theorem 2. The number of interactions M, against tree depth i of FFM and F grows as O M
22 - m, and
0 (M.--122-5 (m, m, for
smooth
amail
m, m74
respectively. Here M-1
0 and cmpty denotes the number far-fiel, smooth field, small field interactions and Tor .mmmooth, mimail m;
22 and m mn
smooth
mi
momall
mi
empty
the number of empty boxes respectively at depth >>
We see that the additional approximations presented in F also impacts memory footprint, as the additional (m. )², mm:moodh, m.!ull terms removes a substantial amount of interactions at each i, significantly slowing down the growth of interactions, reducing memory growth. The efficacy of mmp
(m, y, m. Smmmndh mi mnmmall is widely dependent data. As example, data with points very close to each other would significantly benefit m memmooth more, as the closeness of points would imply more smooth interactions. If points sparsely spread out, (m, m mmmall would provide the most benefit as they remove empty boxes and stops boxes with few points to divide unnecessarily. All proofs can be found in Appendix
Experiments
We demonstrate the utility of F'M over a variety of experiments using the Gaussian kernel k(x, y) =
exp- (( -xxxx=xl
We generate data such that the EV (see section 4.3 varies between 0.1, 1, 10 for data of sizes n = 105, 10", 100 10 The parameters used for FM are 0.1, 0.2, 0.3, 0.5 and 2048. The error for the approximated KMVM product v is r 2D, 30, 10 with a cap ||0-v| where the true KMVM product v is obtained by calculating
2-²
calculated as Relative error == ||
the full KMVM on a subset X' consisting of the first 5000 points in X against the entire dataset in double precision, i.e. v = k(X', X) b, where we fix b ~ (0, In). All experiments were run on NVIDIA V V100-32GB cards, where the data is fitted entirely on the GPU. These cards were chosen since the extra graphic memory is necessary to fit the data on one card when n 10 It should be noted that n 10 can only be run up to D 3, as X and b itself cannot fit in memory for higher dimensions with the GPUs we had available. For details on how F'M scales across multiple GPUs, see Appendix G
Table 2: F' (GPU) compared to results reported in B for FFM(CPU). FM achieves a 90x speed up on a billion data points. F used parameters 64 and 0.5
FFM (12 CPU cores)
F'M (GPU, Ours)
Time (8)
Error
Memory
Time (s)
Error
Memory
Speedup
10"
33.4
1.35 10 4
100 MB
0.08 + 0.00
3- 10 4 4
28 MB
417x
107
169
1.98 10
1GB
1.16 0.04
3- 4 +.104411..-10 10
~ 280 MB
145x
108
1499
1.81 10
10 GB
12.45 + 0.06
2.10 1-110445.00
~.8 GB
120x
109
11340
3.11 10
100 GB
125.90 + 0.52
3- 3-10-1 1.3 10-5
~ 28 GB
90x
KMVM experiments We consider a wide variation of generated datasets to simulate different real- world scenarios to test F'M on. For the k(X, X X)ccsse we consider uniformly and normally distributed data (D 1, 2, 4, 5, 6, 7) together with data simulated from Brownian motion, fractional Brownian motion, and Clustered data (D- 1, 3). For the x(X, Y) ccase we consider uniformly distributed x and normal distributed y (D
5, 6, 7). See Appendix [or visualizations of data. We have to consider smaller n for the k(X, -case when D 2 3, as twice the amount of data needs to be stored. For D= = 3/(4, 5)/(6, 7) we instead consider at most n = 5 10 /2.5 10 /10". It should be noted that D = 7 is hard limit for geometric partitioning- based methods, since for D = 8, we would have 28-2. 28-2 4.3 10 interactions after only divisions. This number of interactions cannot even be represented by a 32-bit integer. We summarize the runs in Table and plot the error and time complexity in Figure 9 for each dataset when D 3. We find that F'M maintains sub-linear empirical complexity up to D- 6, where we have to set small field limit p to a larger number to
2,
3 is kernel agnostic, however choose the Gaussian kernel for simplicity.
8

---

Table 3: Run time and relative error of all KMVM experiments for F.. The slope is computed by regressing
logio (Time (s)) against (n). A slope of I implies (()) scaling.
Himees))
Relative Error
n/D
10f
0.2
0.2
0.2
0.3
6.8
5.0
5.1
00011
000225
00005
01001
00012
OUARS
0X20E
10.1
-022
200.2
20.0
42.4
22.3
2000211
20.0028
2000m0
20.0013
20,0222
200314
2002222
107
0.7
1.1
1.7
4.0
12.7
73.9
88.2
0.0115
000019
0.005
00017
0.0157
0.825
QU2T3
+0.3
-0.4
+0.7
+2.4
+10.0
+58.6
+902
+00021
-0.0225
+00015
+0.0019
=00032
+00261
-0.0199
10
3.5
7.4
17.1
40.6
76.0
525.1
5129
0.02
0.0055
0.0077
0.025
0.0511
0.076
0.0SSE
=2.1
47.3
281.2
+410.
4.71.4
20,0211
20,002
2002
25- 10
NiA
NIA
$40.7
COM
00M
125.9
1247.6
5-10"
N/A
N/A
#28.3
=507.0
COM
00M
N/A
200013
#0.0055
20.03
10P
18.9
20.4
223.4
73.8
1778.
174.0
00M
00M
00M
00M
10.0024
0.004
200027
0.025
10.018
0.0000
N/A
N/A
NAA
N/A
Siige
0.78
0.85
1.99
0,54
639
1.66
1,05
1,11
Table 4: FALKON using default KMVM vs FALKON with F..
FALKON with default KMYM Time (8)
FALKON with F'
Dataset
n
D
M
R²
R²
Time (6)
Emor diff'
Speedup
Uniform
10"
10
0.975 0.034
7631
0.976 0.038
2234
5.31
Normal
10P
3
10a
0.893 + 0.118
7631 +
0.902 + 0.114
2231 = 429
1%
3.41
OSM
10a
2
10a
0.932 1: 0.056
6752 1: 13
0.943 + 0.043
1670 + 48
1.2%
4,04
NYC Taxi
10P
10a
0.526 + 0.029 (AUC)
0963 69
0.526 + 0.030 (AUC)
1535 = 7
1.53
not run out of memory. Further, the error increases in the higher dimensions since we use fewer nodes per dimension when interpolating, owing to the sparse grid technique. We note that D has faster run times than D = 6 which is explained by that for some values of EV, D = 7 doesn' 't run with acceptable errors which skew the run time to datasets where a larger portion of the data can be interpolated.
Rolative lngala)
1 socmmss
Undian, Sope-12
Noncal, Arpe-tr55
Bresias Mitics, alpmr155
Clusenl, skpel D4
Fruckod Boinin Motie, skgeed 18
Ohi Oankgg Aspewilll
Iogul
log
Figure 9: Relative error and time complexity for each 3D dataset.
We further replicate the data used in the first experiment in [3] and compare F'M against FFM (CPU)
in Table
Kernel Ridge Regression experiment We apply F'M to FALKON [24] where we replace their
KMVM operation with F'M and compare performance and speed in solving Kernel Ridge Regression (KRR). The KMVM operation currently used for smaller dimensions is KeOps [I Given some data
want to find the solution a (A(X, X) + AL) b where À is the ridge parameter that stabilizes the inverse. FALKON is a Nyström approximation based solver that requires a subsample X' E MR xd of X to approximate the inverse computation. We focus the experiments on tall and skinny data and take n 100, < with M- = - for all experiments. We consider uniformly and normally sampled data, the Open Street Map (OSM) dataset [[] and classification task the NYC Taxi dataset (2). where we predict whether the customer will tip based on trip distance, trip time and fare cost. To construct b on synthetic problems,
X E N
first take subset D R 1000xd of X and sample a ~ (O0 11000 1000). We then calculate b kXX, D) a + E, where e ~ (0, 0.1). We run KRR for EV 0. 1, 10 synthetic data, and report the average R2 (AUC for NYC Taxi) and training time in Table 4 For the real world datasets OSM and NYC Taxi, we fix the lengthscale using the median heuristic proposed
[6] averaged our results over the runs.
Ablation study between FFM(GPU) and F'M As much of the improved performance can be attributed to our GPU implementation, we conduct an ablation study FFM(GPU) against F'M and KeOps in Table We first present KMVM run times averaged over D = 3 and real-world datasets OSM and NYC Taxi. For KeOps, we only computed the KMVM on uniform data. Since KeOps an exact method, the dataset distribution has no effect computational time. Here, the smoothness

---

criteria and adaptive far-field technique improve computational time. We find that F'M achieves a speed-up between 2.0- 33.3x against FFM()PU) and 8.0 8500x speed-up against KeOps.
Table 5: Comparison between F FFM(GPU) and KeOps. It should be noted that KeOps is only run up to n 10 for all experiments (a run for nt = 10a would take weeks). The times for n > 10 are extrapolated for KeOps. FFM(GPU) could only rum uniform data for n 10
I'M time (s)
FFMGGPU) time 64)
KeOps time 69)
Speedup KeOps
Specdop vs (GPL)
D6M
Taxi
D =
OSM
Tasi (D St
D ==
CSM
Taxi
==3
17 10=2) 2)
(D 5)
(D 21
(D = 3)
(1) = 3)
10"
0.1:000 00
132
02101
0.7
12202
044102
1.56
K0
X.0
2.4
28
107
1.100.4
33 0.6
0.7
10.1 +.8
1000 1.9
3.7 1.0
168.6
57.0
33.3
2.2
30"
3.0
204 8.
17.1 1.3
00M
7095 27.0
13.7.1 12.2
1,4%4
573.0
25.3
26
10"
290.0 08.4
174.0 76.6
0UM
OUM
1,496
$583.0
N6A
N/A
30
D0003
02012
000038
05005
02002
0008
1500
Thooretccal Complexity
(1x logo (
(nnaa
(sF)
Ablation study between F2.5M and F'M We provide and additional ablation study between F2.5M and F in 6 The results are quite similar to the comparison between M and FFM(GPU). Here, we see that smooth field and adaptive far-field approximation (F*M) both improve speed and also memory usage as smooth field helps approximate more interactions. We can thus infer empirically that the mimuoali term in Theorem 2 has a significant impact on reducing the memory footprint of interactions.
Table 6: Comparison between FM, F and KcOps. It should be noted that KcOps is only run up to n 10 for all experiments (a run for n = 100 would take weeks). The times for n > 10 are extrapolated for KeOps.
F'M time (6)
1a time
KcOpa time
Speedup KeOps
Speedup
Taxi
13-3
QSM
13
OSM
1 2)
(2-
(1)
3)
10*
0.1
0.2
0.2
0.3+
0.3
0.3
1,56
8.0
14.0
4.0
1.5
107
0.4
1.7407
1133.0.4
1108120
4.5.1..11
148.6
67.0
10.3
33.6
2.6
102
10.8 38
804=69
1.3
00M
77744 21.5
12.5+18.3
1,494
$73.0
N/A
25.6
2.5
10
101. 36.6
200.
09.4
17600120 26.6
OOM
OOM
466.4.130 15.3
1,0066
$583.0
N/A
N/A
288
0200
3000
Thereziccl Camplesity
Inz. )
(n)
(
Gaussian process regression experiments We further compare F'M as a drop-in KMVM operation applied to Black-box Matrix Multiplication [15] for Gaussian Processes, compared to KISS-GP (35), an approximate Gaussian process using cubic interpolation for kernel approximation. We mimic the setup in 32 and consider the datasets 3DRoad, Song, Buzz and House Electric, where we apply PCA to the last three datasets and take the 3 first principal components for a fair comparison against KISS-GP. which is limited by D < 3. We demonstrate the results in Table 7 As exact GP using F'M demonstrates competitive results even when compared to SVGP [18] and SGPR 30), we hypothesize that many dimensiinnni datasets conform to the manifold hypothesis
13 allowing F M to be
widely applicable out-of-the-bon even in high-dimensional settings.
Table 7: Gaussian process regression results. Exact GP using F M shows improved results and scaling compared to KISS GP. SGPR and KISS- GP could not scale to the Hlustrc dataset.
RMSR
Training SGPR
Danser
(I"ND)
KISS-GP
SGPR
(n
SVGP
(n 1024
Rsact
(F (F'M)
KISS-GP
SVGP
(m 1024)
3DRoad
379.119
0.297 HUB36
0:14- 00.0
0.60 0.000
0481
0.000
27.8 18.0
+1200+108
730.5 + 130
2451 191.4
Song
820,200
0L369 +00029
0.57 +0.298
0801 0.00
0.98 0.000
7.2
1705
473.3 187.5
2373.1 184.9
Buze
375,280
0.67 00002
0.977 -0.05
0.300 0,004
0.304 0.002
335 -9.0
542.7 +08
1754. 1099.5
2780.8 175.6
HouseEletri
1,311,509
0.308 10.00
0OM
00M
0.084 0.008
79.8 23.1
NIA
NIA
22062.6 1. 282.0
6 Limitations and Further Research
This work has introduced and implemented F'M on GPU, which enables fast KMVM for tall and skinny data up to n 10 F'M has improved complexity which also is controllable through 77, and retains linear memory. Experiments in higher dimensions also exhibit linear complexity, however requiring more nodes for lower errors. F can further be directly used as a drop-in KMVM operation, as demonstrated with FALKON and Gaussian process regression, achieving significant speedups and competitive performance on both tasks. As an interpolation based approximation method, F'M is still limited by the exponential growth of interpolation nodes with respect to D, although removing empty boxes, small field and sparse grids allow KMVM for D < A fruitful direction would be to extend ideas in F'M to accommodate higher-dimensional data by considering randomized partitioning 4., decoupling the dependency on D in geometry based partitioning. Further, an exact characterization of how mar m smooth m. small m.; mmmply grows is left future work.
A Acknowledgments
The authors sincerely thank Lood van Nickerk and Jean- Cean-François Ton for their helpful comments.
10

---

References
[1] osm-1billion html Accessed: 2021-05-19.
osm/
[2]
htt.p::/www...nneee gov/
Accessed: 2022-05-14.
[3] M. Aussal and M. Bakry. The fast and free memory method for the efficient computation of convolution kernels, 2019.
[4] A. Backurs, P. Indyk, C. Musco, and T. Wagner. Faster kernel matrix algebra via density estimation. In M. Meila and T. Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 500-510, PMLR, 18-24 Jul 2021.
[5] J. H. Barnes and P. Hut. hierarchical o(n log n) force-calculation algorithm. Nature, 324:446 449, 1986.
[6] J.-M. Belley, P. Belley, E. Colin, and R. Egli. Non-smooth kernels for meshfree methods fluid dynamics. Computers l& Mathematies with Applications, (1:0):1253-12222 2009.
[7] J.-P. Berrut and L. N. Trefethen. Barycentric lagrange interpolation. SIAM Review, 46()):501-- 517, 2004.
[8] S. Börm, M. Fernaneer, and S. Sauter. Variable order, directional h2-matrices for helmholtz problems with complex frequency, 2019.
[9] D. Cai, H. Chow, Y. Saad, and Y. Xi. Smash: Structured matrix approximation by separation and hierarchy, 2017.
[10] Carrier, L. Greengard, and V. Rokhlin. A fast adaptive multipole algorithm for particle simulations. SIAM J. Sci. Stat. Comput., 9(4):669-686 July 1988.
[11] B. Charlier. J. Feydy, J. A. Glaunès, F.-. Collin, and G. Durif. Kerel operations on the GPU, with autodiff, without memory overflows. arXiv preprint arXiv:2004. 2020.
[12] F. Epperson. On the runge example. Am. Math. Monthly, 94(4):329-341, Apr. 1987.
[13] C. Feffermann S. Mitter, and H. Narayanan. Testing the manifold hypothesis. Journal of the American Mathematical Society, 29, 10 2013.
[14] C. for F3M. https//gittbb.oom///tt//////
[15] R. Gardner, Pleiss, D. Bindel, K. Q. Weinberger, and Wilson. Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration, 2018.
[16] D. Garreau, W. Jitkrittum, and M. Kanagawa. Large sample analysis of median heuristic, 2017.
[17] L. Greengard, O'Neil, M. Rachh, and F. Vico. Fast multipole methods for evaluation of layer potentials with locally-corected quadratures, 2020.
[18] Hensman, Fusi, and N. D. Lawrence. Gaussian processes for big data. In Proceedings of the Twenty-Ninth Conference Uncertainty in Artificial Intelligence, UAT' 13, page 282-290, Arlington, Virginia, USA, 2013. AUAI Press.
[19] G. W. Howell. Derivative error bounds for lagrange interpolation: An extension of cauchy's bound for the error of lagrange interpolation. Journal of Approximation Theory, 1991.
[20] W. Kang and C. Wilcox. Mitigating the curse of dimensionality: Sparse grid characteristics method for optimal feedback control and hjb equations, 2015.

---

[21] B. Kohnke, C. Kutzner, A. Beckmann, G. Lube, I. Kabadshow. H. Dachsel, and H. Grubmüller. cuda fast multipole method with highly efficient m21 far field evaluation. The International Journal of High Performance Computing Applications, 35:1049442009964485, 10 2020.
[22] G. K. Leaf and H. G. Kaper. - bounds for multivariate lagrange approximation. SIAM Journal on Numerical Analysis, 1((2):63-338) 3811, 1974.
[23] D. Meagher. Octree encoding: A new technique for the representation, manipulation and display of arbitrary 3-d objects by computer. 10 1980.
[24] G. Mcanti, L. Carratino, L. Rosasco, and A. Rudi, Kernel methods through the roof: Handling billions of points efficiently. In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 14410-14422. Curran Associates, Inc., 2020.
[25] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, Kopf, E. Yang. Z. DeVito, M. Raison, A. Te- jani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala. Pytorch: An imperative style, high-performance deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 8024 8035. Curran Associates, Inc., 2019.
[26] A. Rudi, L. Carratino, and L.. Rosasco. Falkon: An optimal large scale kemnel method, 2018.
[27] B. Scholkopf and J. Smola. Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond. MIT Press, Cambridge, MA, USA, 2001.
[28] C. Schwab and W. L. Wendland. Kernel properties and representations of boundary integral operators. Mathematische Nachrichlen, 1 156(1): 1992.
[29] S. A. Smolyak. Quadrature and interpolation formulas for tensor products of certain class of functions. Dokl. Akad. Nauk SSSR, 142 1963. Transl.: Soviet Math. Dokl. 4:240-243, 1963.
[30] M. Titsias. Variational learning of inducing variables in sparse gaussian processes. In D. van Dyk and M. Welling, editors, Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics, volume 5 of Proceedings of Machine Learning Research, pages 567-574, Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA, 16-18 Apr 2009. PMLR.
[31] D. Tward, T. Brown, Y. Kageyama, J. Patel, Z. Hou, S. Mori, M. Albert, J. Troncoso, and M. Miller, Diffeomorphic registration with intensity transformation and missing data: Applica- tion 3d digital pathology of alzheimer's s disease. Frontiers Neuroscience, 14, 2020.
[32] K. A. Wang, G. Pleiss, J. R. Gardner, S. Tyree, K. Q. Weinberger, and A. G. Wilson. Exact gaussian processes on million data points, 2019.
[33] L. Wang, R. Krasny, and S. Tlupova. A kernel-independen treecode based on barycentric lagrange interpolation, 2019.
[34] A. Wilson and H. Nickisch. Kernel interpolation for scalable structured gaussian processes (kiss-gp). In F. Bach and D. Blei, editors, Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pages 1775-1784, Lille, France, 07-09 Jul 2015. PMLR.
[35] A. G. Wilson and H. Nickisch. Kernel interpolation for scalable structured gaussian processes (kiss-gp). In Proceedings of the 32nd International Conference on International Conference on Machine Learning- Volume 37, ICMLI 15, page 1775-1784. JMLR.org, 2015.
[36] Yang, Duraiswami, Gumerov, and Davis. Improved fast gauss transform and efficient kernel density estimation. In Proceedings Ninth IEEE International Conference on Computer Vision, pages 664-671 vol.1, 2003.
12

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes]
(c) Did you discuss any potential negative societal impacts of your work? [Yes]
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes]
(b) Did you include complete proofs of all theoretical results? [Yes] See Appendix B
D
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material URL)? [Yes]
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes]
(c) Did you report error bars (e.g., with respect to the random seed after running experi- ments multiple times)? [Yes]
(d) Did you include the total amount of compute and the type of resources used (c.g., type of GPUs, internal cluster, cloud provider)? [Yes]
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [N/A]
(b) Did you mention the license of the assets? [N/A]
(c) Did you include any assets either in the supplemental material as URL? [N/A]
(d) Did you discuss whether and how consent was obtained from people whose data you're 're using/curating? [N/A]
(e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
13


Reviews:
[{"review_id": "ztBDNqKmeE3", "reviewer": "Reviewer_TBZE", "summary": "The paper proposes a novel algorithm to approximate kernel matrix-vector multiplication (KMVM), for large kernels and data of dimension less than or equal to 7. Experiments and theoretical analysis suggest the algorithm has linear time and memory complexity. Since KMVM is a key component of other algorithms like conjugate gradients, the authors show their algorithm can be combined with existing large-scale kernel methods and demonstrate significant speed ups.", "questions": "- Maybe mention that $D$ is the dimension of the data in the first paragraph where it is introduced.\n- Are there any assumptions on the type of kernel F$^3$M is useful for? It seems like most experiments were performed on RBF kernels.\n- Line 143 should $u_1$ and $u_2$ be $v_1$ and $v_2$?\n- Can the authors include more detail in Section 4.1 on the low-level GPU optimizations performed in the main text?", "limitations": "No concerns about negative societal impact.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n- The limitations of the method are clearly discussed.\n- The paper makes a strong contribution with a nice algorithmic advance and a technically challenging implementation on GPU.\n- All the main claims of the paper are thoroughly supported by experiments, including an ablation study. Overall a high quality paper.\n\nWeaknesses\n- The algorithm is limited to low-dimensional data, which might hamper the paper's significance in the NeurIPS community. However, this limitation is clearly acknowledged and discussed.\n- I found some sections of the paper hard to follow, particularly Section 4.2, but this may be because I'm not familiar with some of the prior work the paper builds on.\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": "cdPF9frpE5P", "reviewer": "Reviewer_sNpd", "summary": "This paper presents a fast matrix vector multiplication for skinny matrices.\n", "questions": "1. Is it possible to use tensor cores in the proposed algorithm?\n2. How to improve the FLOP efficiency?", "limitations": "As the authors claim, the major limitation of the proposed algorithm is the restriction of the dimension (D <= 7).\nAlthough it may limit the usage, many relational tabular data fit the constraint.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strong Points\n----\n1. This paper targets to solve a well-motivated and important problem.\n2. Related work is clearly discussed.\n3. The empirical results substantially outperform baselines.\n\nWeak Points\n---\n1. The V100 used in the experiments supports tensor cores. I suggest the authors include a discussion of whether tensor cores may be utilized in F^3M.\n2. The GPU utilization of the proposed algorithm is relatively low to the full capacity of V100. I wonder whether the authors could add a discussion of the reason and how to improve the utilization.\n\nPresentation\n----\nThe upper part of Figure 5 is incomplete. And the text is too small for readers.\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": "KEWdZIQ51f", "reviewer": "Reviewer_eghV", "summary": "This paper tries to accelerate the fundamental problem of kernel matrix-vector multiplication (KMVM) operation. The authors propose F^{2.5}M and which is built on classical FFM and specially designed for GPU parallelization. Then they introduce the main algorithm F^{3}M to further accelerate the computation speed. The experiments show that compared with CPU and GPU baselines the new proposed algorithm uses less memory and has a faster speed while controlling the error in an acceptable range.", "questions": "I would like to to see some ablation studies and experiments on F^{2.5}M.", "limitations": "N/A", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "This paper tries to solve the fundamental problem of high memory usage and low speed in large-scale KMVM operations. The writing is clear and the motivation is well-driven. The algorithm is straightforward but seems effective according to the experiment results. The figures in the attachment materials show clearly how the algorithm utilizes the GPU parallelization attribute.\nMinor:\nPictures need to be replaced with vector graphics to improve clarity.", "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": "5gIBTWBF9Ks", "reviewer": "Reviewer_MjYH", "summary": "The paper considers speeding up the matrix-vector multiplication operation in kernel methods for a specific structure (tall and skin) of the kernel matrix. The approximation method builds on the FFM algorithm, which separates data points into far-field and near-field interactions and Lagrange interpolation. The proposed schemes are further implemented on GPUs to exploit their parallel processing power; various speed-ups can be observed compared to the FFM method.", "questions": "The last two items in Table 1 are quite strange as runtime depends on the hardware. What does that mean by running under a minute or an hour?\n\nIt seems that the F3 method is specific to the Gaussian kernel in Section 4.3? It is not clear how the methods depend on hyperparameters when the data set is changed. For example, the suggested parameters in Line 250 are based solely on Figure 7, which may not generalize to new data sets.\n\nTheorem 2 is hard to understand, it needs more discussion. It is not clear how large the different terms (with negative signs) are, relative to the first term in the complexity bound. ", "limitations": "see above", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 1, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The methods contain several levels of approximation and heuristics, with many hyperparameters to be tuned. It is not clear how sensitive these methods are to these parameters. How to best select them for a given dataset?\n\nThe paper is densely packed with approximations and notations, which make it quite difficult to follow at some points. \nIn Section 4.1, the authors repeatedly mentioned **we used non-trivial method...** without any guidance. Please consider updating it. Some parameters are used for a different meaning, e.g., the vector **v** in Line 143 and Line 144.", "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