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 Faster-Fast and Free Memory Method ($\fthreem$) to address these scaling issues of KMVM for tall ($10^8\sim 10^9$) and skinny ($D\leq7$) data. Extensive experiments demonstrate that $\fthreem$ has empirical linear time and memory complexity with a relative error of order $10^{-3}$ and can compute a full KMVM for a billion points 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, improving speed 1.5-5.5 times at the cost of $<1\%$ drop in accuracy. We further demonstrate competitive results on 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 <cit.>, <cit.>, numerical analysis <cit.>, fluid dynamics <cit.>, and machine learning <cit.>. For a dataset of size $n$, KMVM using direct computation has complexity and memory footprint $\mathcal{O}(n^2)$, both unfeasible for modern large scale applications where $n\approx 10^9$ is becoming increasingly common. Pioneering contributions presented in the Fast Multipole Method (FMM) <cit.> amend the complexity of these problems to $\mathcal{O}(n \log{(\epsilon^{-1})})$, where $\epsilon$ is the chosen error tolerance, with varying reductions in memory footprint for data restricted to dimension $D=2$. Subsequent developments in <cit.> mainly focused on extending approximations for a broader set of kernels for a fixed dimensionality $D\leq 3$, tailored for
problems in physics with narrow data such as electrostatics, stellar dynamics, Stokes flow, and acoustic problems, amongst others.

In this paper, we introduce Faster-Fast and Free Memory Method (F$^3$M), a novel algorithm built upon the FFM <cit.> framework to perform KMVM on a GPU for tall and skinny ($D\leq 7$) data of order $n\sim10^9$ in under a minute with user-specified error tolerance, providing between $2-8500$ times speed-up over existing methods.

Notations. We use capital and lower case bold letters to represent matrices and vectors, respectively. In this paper, we will work with matrices $\bfX\in\RR^{n_x\times D}$, $\bfY\in\RR^{n_y\times D}$ and vector $\bfb\in\RR^{n_y}$. For a kernel $k$, the goal for KMVM is to compute $\bfv := \bfK\cdot \bfb$, where $\bfK := k(\bfX,\bfY) = \{k(\bfx_i, \bfy_j)\}_{i=1, j=1}^{n_x, n_y}$, and $\bfx_i, \bfy_j$ denote the $i^{th}, j^{th}$ row of $\bfX,\bfY$ respectively.

§ MOTIVATION AND RELATED WORK

Kernel methods are often limited by their $\mathcal{O}(n^2)$ memory footprint and computational complexity for KMVM. These constraints make scaling beyond $n=10^6$ challenging. Many recent developments have been made to improve both of these constraints, ranging from hardware acceleration using GPUs in KeOps <cit.>, to various approximation techniques proposed in <cit.>. In this work, we focus our attention on kernel independent KMVM methods.

KeOps.   Charlier et al. <cit.> proposes a map-reduce scheme to compute kernels using exactly $\mathcal{O}(n)$ memory and $\mathcal{O}(n^2)$ complexity on GPU. This is achieved by computing the full KMVM product on-the-fly by summing $v_i = \sum_{j=1}^{n} k(\mathbf{x}_i,\mathbf{y}_j)b_j$ directly, without ever storing the kernel matrix $\mathbf{K}\in \mathbb{R}^{n\times n}$ explicitly. Extensive experiments show that this method is practical when $n\leq10^6$, as the GPU hardware acceleration allows the KMVM product to be computed in less than a second on a conventional GPU. Moreover, the method places no constraint on the number of features $D$ it can be applied to, making it favourable for KMVM on medium size datasets. In application contexts, KeOps is currently adopted into conjugate gradient solver FALKON <cit.> as part of the default pipeline.

A brief summary of the evolution of the FMM family. Our contributions in $^*$.

Comparison between methods(* indicate ours).

Method     FMM    KeOps     FFM     $\text{F}^{2.5}\text{M}^*$    $\text{F}^3\text{M}^*$
Kernel Independent

Linear Time

Linear Memory

Restriction in $D$    $\leq3$         $\leq3$    $\leq7$     $\leq7$

c]@l@Scales to $n=10^9$

under $1$ hour [1]

under 1 minute! [1]

[KMVM applied to 3D data for $n=10^9$ on a Nvidia V100 GPU.]

The Fast and Free Memory Method (FFM).  While KeOps can theoretically scale to a billion points, it becomes practically infeasible as the $\mathcal{O}(n^2)$ complexity would imply a computational time of $10^6$ seconds, or roughly $11$ days. To overcome this billion points barrier, Aussal et al. <cit.> deploys a geometric space partitioning scheme, and proposed the Fast and Free Memory Method (FFM), a KMVM approach that extends the FMM <cit.> 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 <cit.>. Compared to KeOps, FFM demonstrates both linear memory and time complexity in experiments and scales to compute a billion-points KMVM on a smaller CPU cluster under $4$ hours, outscaling the GPU implementation of FMM <cit.>. While $4$ hours is a 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$ <cit.>, both FMM and FFM can only be applied to $D\leq3$ data, a price to pay for the speed-up of KMVM operations when $n=10^9$. Furthermore, we show in our experiments that a direct FFM port to GPU gives unstable results for $n=10^9, D=3$ for non-trivial data simulations (bottom row in Appendix <ref>).

Our contribution. To surpass the billion point barrier while maintaining high-speed and stable computation, we propose $\fhalfm$ and our main algorithm $\fthreem$, the first pair of KMVM algorithms that can reliably scale to $n=10^9$ on skinny data using a single GPU. We build $\fhalfm$ 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 $\fthreem$, which runs $2.0-33.3$ times quicker and more stable than a direct port of FFM on GPU. See Fig. <ref> and Table <ref> for an overview and comparisons of the methods. We summarise our contribution as follows:

1.  We propose Faster-Fast and Free Memory Method $(\fthreem)$, 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 $n\leq 10^9$ data in under a minute. Codebase is released here <cit.>.

2.  We characterize theoretical time and memory complexity of $\fthreem$.

3.  We run extensive KMVM experiments of $\fthreem$ on a variety of tall and skinny data with $n\leq10^9$, demonstrating empirical linear time and memory scaling, and achieving speedups between 2–8500 times when compared to FFM (GPU and CPU) and KeOps.

4.  We run a practical application of $\fthreem$ as a drop-in replacement for KeOps in conjugate gradient solver FALKON <cit.> for kernel ridge regression and classification (KRR) on 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 <cit.>, SVGP <cit.> and SVGR <cit.> with significant speed-ups.

§ BACKGROUND

The FFM method considers KMVM for a kernel $k$ evaluated on two data matrices $\mathbf{X},\mathbf{Y}$ and $\mathbf{b}$ are weights associated with $\mathbf{Y}$. The KMVM is expressed as $\textbf{v} :=  k(\mathbf{X},\mathbf{Y})\cdot {\mathbf{b}}=\mathbf{K}\cdot \mathbf{b}$. For example, $\mathbf{b}$ could be the weights in a KRR or the strength of electronic charges. As $n_x$ and $n_y$ are taken to be very large, a full computation is unfeasible. In this section, we illustrate and detail the main steps of FFM, before presenting our improvements in Section 4.

Enclosing $\mathbf{X}$ and $\mathbf{Y}$ within a box, "x" marks the center of the box. Numbers under the the boxes denotes edge length. In the right plot, we have enclosed the blue points with the largest box.

For illustration purposes, we first consider a simple 2D KMVM. Our goal is to calculate $k(\mathbf{X},{\mathbf{Y}})\cdot \mathbf{b}$ for $\mathbf{X},\mathbf{Y}$ in <Ref>. The intuition behind FFM is to reduce the complexity of calculating the full KMVM by partitioning $\mathbf{X}$ and
${\mathbf{Y}}$ such that certain calculations can be approximated in a 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 a large enough box that can just enclose $\mathbf{X}$ or ${\mathbf{Y}}$. The edge length of this box is calculated as

$$\mathcal{E}:= \max\left(\max_d(x_{\max}^{(d)}-x_{\min}^{(d)}),\max_d(y_{\max}^{(d)}-y_{\min}^{(d)})\right)$$
where $x_{\max}^{(d)},x_{\min}^{(d)}$ denotes the largest value and the smallest value along the $d$-dimension in $\mathbf{X}$ and similarly for ${\mathbf{Y}}$. <Ref> illustrates this enclosing procedure.

Defining near and far-field  In FMM, an octree <cit.> is applied to recursively partition data into

Recursive partitioning of $\mathbf{X}$ and $\mathbf{Y}$ for 2D data. Far-field interactions are colored green while near-field interactions are colored orange and $\delta$ denotes the euclidean distance between the centers.

smaller boxes $B_p^X \subset \mathbf{X}, B_q^Y \subset \mathbf{Y}$, with $p,q$ denoting box indices. Here each box corresponds to a subset of rows in the data matrix. Let us also denote $\mathbf{b}_q$ as the partition of $b_j$'s grouped with the same indices as $B_q^Y$. To calculate the KMVM between two boxes $B_p^X,B_q^Y$ with the grouped vector $\mathbf{b}_q$, for each $\mathbf{x}_i \in B_q^X$, we compute

\begin{align}
v_i^{p,q} = \sum_{\mathbf{y}_j\in B_q^Y, b_j \in \mathbf{b}_q} k(\mathbf{x}_i,\mathbf{y}_j)b_j
\end{align}

with $\mathbf{v}^{p, q}=[v_1^{p,q}\hdots v_{n_x}^{p,q}]$. Now the target $\mathbf{v}$ can be computed as $\mathbf{v}=\sigma([\mathbf{v}^{p=1},\hdots,\mathbf{v}^{p=P}]^\top)$, where $\mathbf{v}^{p}=\sum_{q=1}^Q\mathbf{v}^{p,q}$, $P,Q$ denote the total number of boxes and $\sigma(\cdot)$ a permutation such that $v_i^p$ appear in the same order as $\mathbf{x}_i$ appears in $\mathbf{X}$. <Ref> shows how boxes are recursively partitioned.

Far and near-field interactions  FFM relies on a divide-and-conquer
strategy to effectively compute a 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 approximated using Lagrange interpolation for speed, explained in the paragraph below. The partitioning procedure in FFM is recursive, where the recursion depth $\texttt{tree\_depth}$ controls the size of the edge $l=\frac{\mathcal{E}}{2^{\texttt{tree\_depth}}}$ 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 $2l$, i.e. $\texttt{IsFarField}:=\Vert \mathbf{x}_{\text{center}}-\mathbf{y}_{\text{center}} \Vert \geq 2l$. While $l$ for each box will decrease with the number of divisions, this rule ensures a fixed minimal distance for a given depth for far-field interactions. <Ref> illustrates how far(green) and near(orange)-field interactions arise between $\mathbf{X}$ and ${\mathbf{Y}}$ when $\texttt{tree\_depth}$ increases.

Lagrange interpolation   We review Lagrange interpolation used for far-field approximations in FFM. Given a function $f(x):[-1,1] \to \mathbb{R}$ and $r+1$ unique points $s_i \in [-1,1],\quad i=0,\hdots,r$, there exists a unique polynomial $p_r(x)$ of degree $\leq r$ that interpolates $f$ at $p_r(s_i)=f(s_i)$. The Lagrange polynomial is given by $p_{r}(t)=\sum_{i=0}^{r} f(s_i) \mathcal{L}_{i}(t),$
where $\mathcal{L}_{i}(t)=\frac{\prod_{j=0, j \neq i}^{r}\left(t-s_{j}\right)}{\prod_{j=0, j \neq i}^{r}\left(s_{i}-s_{j}\right)}, \quad i=0, \ldots, r.$
We are free to chose the degree $r$ as well as the points $s_i$ to interpolate through. The choice of $s_i$ is especially important in minimizing large oscillations around the edges of the interpolation interval (Runge's phenomenon <cit.>). For this reason, Chebyshev nodes of the second kind are used <cit.>
$s_{i}=\cos \theta_{i},\text{ where } \theta_{i}=\frac{i \pi}{r}, \quad i=0, \ldots, r$.

Interpolating $k(\mathbf{x},\mathbf{y})$   By noticing that $k(\mathbf{x},\mathbf{y})$ is a bivariate function, we can apply Lagrange interpolation twice, thus interpolating $k(\mathbf{x},\mathbf{y})$ as
$k(\mathbf{x}, \mathbf{y}) \approx \sum_{i=1}^{r_{X}} \mathcal{L}_{i}(\mathbf{x}) \sum_{j=1}^{r_{Y}} k\left(\mathbf{s}^x_{i}, \mathbf{s}^y_{j}\right) \mathcal{L}_{j}(\mathbf{y}).$ Here $r_X,r_Y$ denotes the number of the interpolation nodes and $\mathbf{s}^x_{i},\mathbf{s}^y_{j}\in \mathbb{R}^D$ denotes the grid of interpolation nodes for $B_p^X$ and $B_q^Y$. Note that since $\mathbf{x}\in \mathbb{R}^D$, we take
$\mathcal{L}_{i}(\mathbf{x}):= \prod_{d=1}^D \frac{\prod_{j=0, j \neq i}^{r}\left(x^{(d)}-s^{(d)}_{j}\right)}{\prod_{j=0, j \neq i}^{r}\left(s^{(d)}_{i}-s^{(d)}_{j}\right)}, \quad i=0, \ldots, r.$ These operations can be vectorized and computed sequentially on-the-fly with linear memory footprint
$\mathbf{v} \approx \mathbf{L}_{X}^{T} \cdot \overbrace{(\mathbf{K} \cdot\underbrace{\left(\mathbf{L}_{Y} \cdot \mathbf{b}\right)}_{\mathbf{v}_1})}^{\mathbf{v}_2},$
which is done by first computing $\textbf{v}_1$, then $\textbf{v}_2$ and lastly $\textbf{v}$. A far-field KMVM between two boxes $\textbf{v}_{p,q}=k(B_p^X,B_q^Y)\cdot \textbf{b}_q$ is then approximated by using double Lagrange interpolation according to <Ref>.

Here we approximate a far-field interaction between two boxes. In $\textbf{1.}$ we first normalize the data between $[-1,1]$ and conduct 2D interpolation of $k(\cdot, \mathbf{y})$ while summing over $\mathbf{b}$. In 2., the interpolation for $k(\mathbf{x},\cdot)$ will also be normalized and hence we need to translate the distance between the boxes by calculating $\mathbf{K}\cdot \mathbf{v}_1$. Lastly in 3. we interpolate $k(\mathbf{x},\cdot)$ while we sum $\mathbf{v}_2$.

§ FASTER-FFM ($\FHALFM$ AND $\FTHREEM$)

To fully leverage the port of FFM to GPU, we enhance FFM with novel approximation procedures for improved complexity and memory optimizations to scale to $n=10^9$. We coin this improved version Faster-FFM ($\fthreem$). The capabilities of $\fthreem$ against previous methods are summarized in <Ref>.

§.§ 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:Box-to-threadblock 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-trivial 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 1 consists of points with $i\in [1,\hdots,500]$ and box 2 with  $i\in [501,\hdots,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 <ref>. No native sorting methods – We found that LibTorch <cit.> sorting methods often led to out-of-memory (OOM) due to allocation of large long-type vectors on GPU. When $n=10^9$, this implies

Here we assume $D=2$, hence we only divide each box $2^D=2$ times each time.

allocating 8GB of memory, 25% of the 32GB card used, making it a necessity to avoid native sorting methods. In-place grouping data on boxes – Due to infeasible LibTorch sorting methods, we additionally had to design an algorithm that finds a permutation that would group $\mathbf{X}$ into its corresponding boxes in linear time and memory. We used a count and increment-based strategy that would: 1. Count the number of points in each box during the assignment operation ($\mathcal{O}(n)$) and store the count in a vector $\xi$. Then run a cumulative sum over $\xi$, starting from 0. 2. Initialize a $n$ long permutation vector $\pi$. Using the counting vector $\xi$, we would re-run the assignment operation and arrange a point with index $i$ as following $\pi[\texttt{atomicAdd}(\xi'[\texttt{box\_index}],1)+\xi[\texttt{box\_index}]]=i$, where $\xi'$ is a running count of points in each box. We specifically have to use the function  to increment the count for each box in parallelized GPU environments to avoid thread locks.

We refer to the $\texttt{box\_division\_cum\_hash}$ and $\texttt{box\_division\_assign\_hash}$ function in $\texttt{n\_tree.cu}$ for exact details.Ensuring interactions are sorted – To avoid any unnecessary sorting, we ensure that the matrix containing interactions is always sorted by recursively dividing old interactions. We illustrate the procedure in <Ref>.

KMVM times on several datasets between   FFM(GPU) vs   $\fthreem$. U&N was only run up to $n=5\cdot 10^8$

We refer to the $\texttt{get\_new\_interactions}$ function in $\texttt{n\_tree.cu}$ for the exact implementation.

However, we found that these optimizations and porting alone were not enough to scale to $n=10^9$ on 3D datasets, as <Ref> demonstrates. FFM doesn't remove empty boxes or handle boxes with few points in them and keeps exponentially creating new empty boxes and interactions, thus leading to out-of-memory (OOM) errors on non-uniform data (see Appendix <ref>).

§.§ Scaling to $n=10^9$ on GPU ($\fhalfm$)

In this section, we detail the memory enhancements that allow $\fthreem$ to consistently scale to $n=10^9$.

Removing empty boxes with hash list indexing  To ensure linear memory on GPU, we only keep a reindexing vector $\sigma$ of size $n_x$(resp. $n_y$) in memory during the computation of the algorithm in addition to 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 $\texttt{tree\_depth}$, this method quickly becomes pathological. To amend this, we propose a linear complexity formula to retrieve the box index $\beta_i$ a point $\mathbf{x}\in\mathbf{X}\subset \mathbb{R}^D$ belongs to
$\beta_i = \underbrace{\sum_{d=1}^D 2^{\texttt{tree\_depth}\cdot (d-1)}}_{\text{Summing over $D$ dimensions}}  \cdot \underbrace{\lfloor 2^{\texttt{tree\_depth}} \frac{x_d -  \alpha_d}{\mathcal{E}}   \rfloor}_{ \substack{\in \{0,1\}, \text{ Denotes left or right}\\\text{of center of box edge}}
where $\alpha_d$ denotes the minimum value of $\mathbf{X}$ in dimension $D$ and $x_d$ is the value of $\mathbf{x}$ 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, we use a hash list to store $\beta_i$ and the order $i$. We can then group points $\{\mathbf{x}_i\}_{i=1}$ to their respective ordering $i$ using the hash list in $\mathcal{O}(n)$ time in contrast to $\mathcal{O}(n\cdot2^{D\cdot \texttt{tree\_depth}})$ by direct computation.

Handling boxes with few points with small field  In cases when the number of points in each box can 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_p^X,B_q^X$ if both have a small number of points, i.e. if

Full grid vs sparse grid.

$|B_p^X|+|B_q^X|\leq \rho$,
for some threshold number $\rho$. To minimize the computations needed, $\rho$ can be set to $\rho=r_X+r_Y$. This intuitively allows $\fthreem$ to directly compute interactions that are too small to benefit from interpolation savings (i.e. $|B_p^X|+|B_q^X|\leq r_X+r_Y$), thus limiting memory usage by stopping partitions from dividing further than necessary. In higher dimensions where the division rate is faster, $\rho$ can be set to a higher value to limit memory usage at the expense of more direct computations which are slower.

Sparse grids  As the number of Lagrange polynomials increases exponentially with dimension, we implement sparse grids <cit.> to allow for a finer selection of interpolation nodes. With sparse grids, the number of nodes needed grows slower <cit.>, thus saving memory. We give an example of a sparse grid versus a full grid in 2D in <Ref>.

§.§ Speeding up $\fhalfm$ $\left(\fhalfm \to \fthreem\right)$

Smoothness criteria  FFM speeds up its computations with minimal loss in accuracy by selectively interpolating interactions that are far apart. To improve speed, we introduce the smoothness criterion to widen the selection of interactions that can be interpolated with minimal loss in accuracy. For a Gaussian Kernel $k(\mathbf{x},\mathbf{y})=\exp{\left(\frac{\Vert \mathbf{x}-\mathbf{y}\Vert^2}{2\gamma^2}\right)}$, with lengthscale $\gamma$, the smoothness criteria is defined as $\texttt{is\_smooth}:=\frac{\frac{1}{N_{B_p^X}}\sum_{\mathbf{x}_i\in B_p^X}\sum_{d=1}^D\Vert x_i^{(d)} - \Bar{x}^{(d)} \Vert^2}{2\gamma^2} +\frac{\frac{1}{N_{B_q^Y}}\sum_{\mathbf{y}_j\in B_q^Y}\sum_{d=1}^D\Vert y_j^{(d)} - \Bar{y}^{(d)} \Vert^2}{2\gamma^2}\leq \eta$

Dimension $d=3$

Dimension $d=5$

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 a gaussian kernel. We observe that its not beneficial to interpolate at all when the square distance exceeds 5 and that is is sufficient to only use $3^D$ nodes when square distance is $\leq 0.01$.

between an adjacent interaction of boxes $B_p^X$, $B_q^Y$. The quantity computed can be understood as “Effective Variance” (EV), as it considers total variation in the exponent of the Gaussian kernel. We justify the smoothness criteria with the following proposition.

Consider $\mathbf{x},\mathbf{y} \in \mathcal{X} \subset \mathbb{R}^d$ such that $d(\mathbf{x},\mathbf{y}):=\frac{\Vert \mathbf{x}-\mathbf{y}\Vert^2}{2\gamma^2}=\frac{1}{2\gamma^2}\sum_{i}^d(x^{(i)}- y^{(i)})^2 \leq \eta<1$ for all $\mathbf{x},\mathbf{y}$. When interpolating $k(\mathbf{x},\mathbf{y}) = \exp\left(-d(\mathbf{x},\mathbf{y})\right)$ using bivariate Lagrange interpolation $\mathcal{L}_r(\mathbf{x},\mathbf{y}):=\mathbf{L}_{X}^{T} \cdot \mathbf{K} \cdot \left(\mathbf{L}_{Y} \cdot \mathbf{b}\right)$ with degree $r=2p$, for any $p\in \mathbb{N}_{>0}$ there exist nodes $\mathbf{s}^{\mathbf{x}},\mathbf{s}^{\mathbf{y}}$ for $\mathcal{L}_r(\mathbf{x},\mathbf{y})$ such that the pointwise interpolation error is bounded by $\mathcal{O}(\eta^{p+1})$.

See Appendix <ref> for proof. Hence for small $\eta < 1$, we see that the error becomes small for well specified $\mathcal{L}_r(\mathbf{x},\mathbf{y})$.  To avoid calculating the sample variance during computations which costs $\mathcal{O}(n)$, we exploit that data is partitioned into hypercubes with a known edge $\mathcal{E}$ and take the upper bound of the variance in each cube as $\frac{\mathcal{E}^2}{4}$ along a dimension. A proof for this bound is provided in <Ref>. Adjacent interactions are then classified as smooth when $\sum_{d}^D \frac{\mathcal{E}^2}{4\cdot 2 \gamma} (\text{EV of }B_p^X) + \frac{\mathcal{E}^2}{4\cdot 2 \gamma} (\text{EV of }B_q^Y) =  \frac{D\mathcal{E}^2}{\gamma^2\cdot 4}\leq \eta$ which only costs $\mathcal{O}(1)$ to compute.

Adaptive far-field approximation   To further improve speed we introduce an 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 <cit.>, however, these bounds cannot be directly used to create an adaptive interpolation rule. We thus simulate KMVM errors for $k(\mathbf{X},\mathbf{Y})\cdot \mathbf{b}$ where $\mathbf{X},\mathbf{Y}$ are uniformly distributed and $\mathbf{b}$ is normally distributed. We fix a distance between $\mathbf{X}$ and $\mathbf{Y}$ and vary the squared of this distance between boxes against nodes in <Ref>. We use a Gaussian Kernel with $\gamma=\frac{1}{\sqrt{2}}$.

Based on <Ref>, we use the following rule for selecting the number of interpolation nodes for far-field interactions

\begin{equation*}
\begin{split}
r_{\text{far}}(r)= \left\{\begin{array}{ll}
\min(r,3^D) &\text{ if }  \left (\frac{\mathcal{E}}{2^{c_{\text{depth}}}}\right)^2\cdot \frac{1}{2\gamma^2}\leq 0.01  \\
r &\text{ if } 0.01<\left (\frac{\mathcal{E}}{2^{c_{\text{depth}}}}\right)^2\cdot \frac{1}{2\gamma^2}\leq 5 \\
0 &\text{ if } 5<\left (\frac{\mathcal{E}}{2^{c_{\text{depth}}}}\right)^2\cdot \frac{1}{2\gamma^2} \\
\end{array}\right.
\end{split}
\end{equation*}

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 <cit.> evaluated at the Chebyshev nodes of the second kind. As this is a 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-invariant kernel by recalculating the Taylor expansion for smoothness criteria and rerunning the simulation for adaptive far-field approximation.

§.§ Complexity

The time complexity of FFM is $\mathcal{O}(n\log{(n)})$ <cit.> and we use a similar derivation strategy for $\fthreem$ to obtain a complexity that is dependent on the effective variance limit $\eta$ (chosen parameter) and the box width $\mathcal{E}$ (data). We first present two propositions needed to derive the complexity of $\fthreem$.

A far-field interaction between two boxes containing $n_x$ and $n_y$ points respectively has time complexity $\mathcal{O}(n)$, where $n=\max(n_x,n_y)$.

Given $n$ data points in dimension $D$, the maximum number of divisions $\text{Tree}_{\texttt{max divisions}}$ is given by

\begin{equation}
\text{Tree}_{\texttt{max divisions}} = \log_{2^D}(n).
\end{equation}

With the above results, the complexity of FFM is taken as the maximum number of divisions multiplied by the complexity of far-field interactions at each division which yields $\mathcal{O}(n\log{(n)})$. We remark that near-field interactions between boxes containing only 1 data point have linear time complexity, hence the results hold.

Given a KMVM with edge $\mathcal{E}$ (dependent on data $\mathcal{X},\mathcal{Y}$), lengthscale $\gamma$, effective variance limit $\eta$, $n$ data points and data dimension $D$, $\fthreem$ has time complexity $\mathcal{O}(n\cdot \log_2\left(\frac{D \cdot  \mathcal{E}^2}{\gamma^2 \cdot 4 \cdot \eta} \right))$, which can be taken as $\mathcal{O}(n\cdot \log_2\left(\frac{C}{\eta} \right))$ where $C\propto \frac{D\cdot \mathcal{E}^2}{\gamma^2}$.

Memory footprint  As our implementation uses the same partitioning strategy as FFM, the theoretical memory complexity remains $\mathcal{O}(n)$ for $\fthreem$ (see <cit.> for proof). However, this does not accurately reflect the memory footprint of the actual implementations, whose memory mostly depends on the number of interactions stored. We summarize these memory footprints for FFM and $\fthreem$ in <Ref> below.

The number of interactions $M_i$ against tree depth $i$ of FFM and $\fthreem$ grows as $\mathcal{O}\left(M_{i-1}2^{2\cdot D }- m_{i}^{\text{far}} \right)$ and
$$\mathcal{O}\left(M_{i-1} 2^{2\cdot D}- (m_i^{\text{empty}})^2 - m_{i}^{\text{far}}-m_{i}^{\text{smooth}}-m_{i}^{\text{small}} \right)$$
respectively.  Here $M_{-1} = \frac{1}{2^{2D}}$ and $m_{0}^{\text{far}}=m_{0}^{\text{smooth}}=m_{0}^{\text{small}}=m_0^{\text{empty}}=0$ and $m_{i}^{\text{far}},m_{i}^{\text{smooth}},m_{i}^{\text{small}},m_i^{\text{empty}}$ denotes the number of far-field, smooth field, small field interactions and the number of empty boxes respectively at depth $i>0$.

We see that the additional approximations presented in $\fthreem$ also impacts memory footprint, as the additional $(m_i^{\text{empty}})^2,m_{i}^{\text{smooth}},m_{i}^{\text{small}}$ terms removes a substantial amount of interactions at each $i$, significantly slowing down the growth of interactions, reducing memory growth. The efficacy of $(m_i^{\text{empty}})^2,m_{i}^{\text{smooth}},m_{i}^{\text{small}}$ is widely dependent on data. As an example, data with points very close to each other would significantly benefit $m_{i}^{\text{smooth}}$ more, as the closeness of points would imply more smooth interactions. If points are sparsely spread out, $(m_i^{\text{empty}})^2,m_{i}^{\text{small}}$ 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 <Ref>.

§ EXPERIMENTS

We demonstrate the utility of $\fthreem$ over a variety of experiments using the Gaussian kernel $k(\mathbf{x},\mathbf{y})=\exp{-(\frac{\Vert \mathbf{x}-\mathbf{y}\Vert^2}{2\gamma^2})}$.[$\fthreem$ is kernel agnostic, however we choose the Gaussian kernel for simplicity.] We generate data such that the EV (see section <ref>) varies between $0.1, 1, 10$ for data of sizes $n=10^6,10^7,10^8,10^9$. The parameters used for $\fthreem$ are $\eta=0.1,0.2,0.3,0.5$ and $r=2^D,3^D,4^D$ with a cap at $r=2048$. The error for the approximated KMVM product $\hat{\mathbf{v}}$ is calculated as $\text{Relative error}:= \frac{\Vert \hat{\mathbf{v}}-\mathbf{v}\Vert^2}{\Vert\mathbf{v}\Vert^2}$, where the true KMVM product $\mathbf{v}$ is obtained by calculating the full KMVM on a subset $\mathbf{X}'$ consisting  of the first 5000 points in $\mathbf{X}$ against the entire dataset in double precision, i.e. $\mathbf{v}=k(\mathbf{X}',\mathbf{X})\cdot \mathbf{b}$, where we fix $\mathbf{b} \sim \mathcal{N}(0,I_{n})$.
All experiments were run on NVIDIA 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^9$. It should be noted that $n=10^9$ can only be run up to $D=3$, as $\mathbf{X}$ and $\mathbf{b}$ itself cannot fit in memory for higher dimensions with the GPUs we had available. For details on how $\fthreem$ scales across multiple GPUs, see Appendix <ref>.

$\fthreem$ (GPU) compared to results reported in <cit.> for FFM(CPU). $\fthreem$ achieves a $90\times$ speed up on a billion data points. $\fthreem$ used parameters $r=64$ and $\eta=0.5$

3c|FFM (12 CPU cores)     3c$\fthreem$ (GPU, Ours)
n          Time (s)        Error                      Memory        Time (s)              Error                                     Memory            Speedup
$10^6$     $33.4$          $1.35\cdot 10^{-4}$        100 MB        $0.08\pm0.00$         $3\cdot 10^{-4}\pm 7 \cdot 10^{-5} $      $\sim$ 28 MB      $417\times$

$10^7$     $169$           $1.98\cdot 10^{-4}$        1GB           $1.16\pm0.04$         $3\cdot 10^{-4}\pm 1.2 \cdot 10^{-4}$     $\sim$ 280 MB     $145\times$

$10^8$     $1499$          $1.81\cdot 10^{-4}$        10 GB         $12.45\pm 0.06$       $2\cdot 10^{-4}\pm 5 \cdot 10^{-5}$       $\sim$ 2.8 GB     $120\times$

$10^9$     $11340$         $3.11\cdot 10^{-4}$        100 GB        $125.90 \pm 0.52$     $3\cdot 10^{-4}\pm 1.3 \cdot 10^{-5}$     $\sim$ 28 GB      $90\times$

Run time and relative error of all KMVM experiments for $\fthreem$. The slope is computed by regressing $\log_{10}(\text{Time (s)})$ against $\log_{10}(n)$. A slope of 1 implies $\mathcal{O}(n)$ scaling.

7c|Time(s)                                                                                7cRelative Error

$n / D$                     1                              2                               3                                4                                 5                                 6                                 7                                 1                                   2                                   3                                   4                                   5                                   6                                   7
$10^6$                      $\makecell{0.2\\ \pm0.1}$      $\makecell{0.2\\ \pm0.2}$       $\makecell{0.2\\ \pm0.1}$        $\makecell{0.4\\ \pm0.2}$         $\makecell{0.9\\ \pm0.9}$         $\makecell{3.0\\ \pm2.4}$         $\makecell{3.1\\ \pm2.1}$         $\makecell{0.0018\\ \pm0.0021}$     $\makecell{0.0023\\ \pm0.0028}$     $\makecell{0.0005\\ \pm0.0009}$     $\makecell{0.0014\\ \pm0.0013}$     $\makecell{0.0022\\ \pm0.0022}$     $\makecell{0.0303\\ \pm0.0314}$     $\makecell{0.0294\\ \pm0.0282}$

$10^7$                      $\makecell{0.7\\ \pm0.3}$      $\makecell{1.1\\ \pm0.4}$       $\makecell{1.7\\ \pm0.7}$        $\makecell{4.0\\ \pm2.4}$         $\makecell{12.7\\ \pm10.0}$       $\makecell{79.9\\ \pm88.6}$       $\makecell{88.2\\ \pm90.2}$       $\makecell{0.0016\\ \pm0.0021}$     $\makecell{0.0019\\ \pm0.0023}$     $\makecell{0.0006\\ \pm0.0013}$     $\makecell{0.0017\\ \pm0.0019}$     $\makecell{0.0057\\ \pm0.0032}$     $\makecell{0.0325\\ \pm0.0261}$     $\makecell{0.0273\\ \pm0.0179}$

$10^8$                      $\makecell{3.5\\ \pm1.1}$      $\makecell{7.4\\ \pm2.1}$       $\makecell{17.1\\ \pm7.3}$       $\makecell{40.6\\ \pm28.7}$       $\makecell{76.0\\ \pm83.2}$       $\makecell{525.1\\ \pm410.9}$     $\makecell{512.9\\ \pm471.4}$     $\makecell{0.002\\ \pm0.0026}$      $\makecell{0.0025\\ \pm0.0032}$     $\makecell{0.0007\\ \pm0.0013}$     $\makecell{0.0023\\ \pm0.0021}$     $\makecell{0.0051\\ \pm0.0032}$     $\makecell{0.0376\\ \pm0.0239}$     $\makecell{0.0444\\ \pm0.017}$

$2.5 \cdot 10^8$            N/A                            N/A                             N/A                              $\makecell{76.7\\ \pm25.9}$       $\makecell{340.7\\ \pm247.6}$     OOM                               OOM                               N/A                                 N/A                                 N/A                                 $\makecell{0.0034\\ \pm0.0029}$     $\makecell{0.0076\\ \pm0.0012}$     N/A                                 N/A

$5 \cdot10^8$               N/A                            N/A                             $\makecell{74.9\\ \pm28.5}$      $\makecell{289.4\\ \pm222.8}$     $\makecell{631.5\\ \pm997.0}$     OOM                               OOM                               N/A                                 N/A                                 $\makecell{0.0012\\ \pm0.0013}$     $\makecell{0.0023\\ \pm0.0026}$     $\makecell{0.0041\\ \pm0.003}$      N/A                                 N/A

$10^9$                      $\makecell{29.4\\ \pm8.9}$     $\makecell{73.8\\ \pm23.4}$     $\makecell{174.0\\ \pm76.6}$     OOM                               OOM                               OOM                               OOM                               $\makecell{0.0024\\ \pm0.0024}$     $\makecell{0.0025\\ \pm0.0027}$     $\makecell{0.0009\\ \pm0.0018}$     N/A                                 N/A                                 N/A                                 N/A
Slope                       0.78                  0.85                   0.99                    0.98                     0.99                     1.06                     1.02

$\mathcal{O}(n\log(n))$     7c|1.11

KMVM experiments  We consider a wide variation of generated datasets to simulate different real-world scenarios to test $\fthreem$ on. For the $k(\mathbf{X},\mathbf{X})$-case we consider uniformly and normally distributed data ($D=1,2,3,4,5,6,7$) together with data simulated from Brownian motion, fractional Brownian motion, and Clustered data ($D=1,2,3$). For the $k(\mathbf{X},\mathbf{Y})$-case we consider uniformly distributed $\bfx$ and normal distributed $\bfy$ ($D=1,2,3,4,5,6,7$). See Appendix <ref> for visualizations of data. We have to consider smaller $n$ for the $k(\mathbf{X},\mathbf{Y})$-case when $D\geq3$, as twice the amount of data needs to be stored. For $D=3/(4,5)/(6,7)$ we instead consider at most $n=5\cdot 10^8/2.5\cdot 10^8/10^8$. It should be noted that $D=7$ is a hard limit for geometric partitioning-based methods, since for $D=8$, we would have $2^{8\cdot 2} \cdot 2^{8\cdot 2}\approx 4.3 \cdot 10^9$ interactions after only 2 divisions. This number of interactions cannot even be represented by a 32-bit integer. We summarize the runs in <Ref> and plot the error and time complexity in <Ref> for each dataset when $D=3$. We find that $\fthreem$ maintains sub-linear empirical complexity up to $D=6$, where we have to set small field limit $\rho$ to a larger number to 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=7$ 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.

Time complexity and relative error for each 3D dataset.

We further replicate the data used in the first experiment in <cit.> and compare $\fthreem$ against FFM (CPU) in <Ref>.

Kernel Ridge Regression experiment  We apply $\fthreem$ to FALKON <cit.>, where we replace their KMVM operation with $\fthreem$ and compare performance and speed in solving Kernel Ridge Regression (KRR). The KMVM operation currently used for smaller dimensions is KeOps <cit.>. Given some data $\mathbf{X}\in \mathbb{R}^{N\times d}$ we want to find the solution $\mathbf{\alpha} = \left( k(\mathbf{X},\mathbf{X})+\lambda I \right)^{-1}\mathbf{b}$ where $\lambda$ is the ridge parameter that stabilizes the inverse. FALKON is a Nyström approximation based solver that requires a subsample $\mathbf{X}'\in \mathbb{R}^{M\times d}$ of $\mathbf{X}$ to approximate the inverse computation. We focus the experiments on tall and skinny data and take $n=10^9, d\leq 3$ with $M=10^5$ for all experiments. We consider uniformly and normally sampled data, the Open Street Map (OSM) dataset <cit.> and a classification task on the NYC Taxi dataset <cit.>, where we predict whether the customer will tip based on trip distance, trip time and fare cost. To construct $\mathbf{b}$ on synthetic problems, we first take a subset $\mathcal{D}\in \mathbb{R}^{1000\times d}$ of $\mathbf{X}$ and sample $\alpha \sim \mathcal{N}(0,I_{1000\times 1000})$. We then calculate $\mathbf{b}=k(\mathbf{X},\mathcal{D})\cdot \alpha + \varepsilon$, where $\varepsilon\sim \mathcal{N}(0,0.1)$. We run KRR for $EV=0.1,1,10$ on synthetic data, and report the average $R^2$ (AUC for NYC Taxi) and training time in <Ref>. For the real world datasets OSM and NYC Taxi, we fix the lengthscale using the median heuristic proposed in <cit.> averaged our results over the 3 runs.

FALKON using default KMVM vs FALKON with $\fthreem$.

2c|FALKON with default KMVM     2c|FALKON with $\fthreem$

Dataset      $n$        $D$     $M$        $R^2$                      Time (s)                 $R^2$                      Time (s)               Error diff     Speedup

Uniform      $10^9$     3       $10^5$     $0.975\pm 0.034$           $7631\pm 2$              $0.976\pm 0.038$           $2234\pm 429$          0%            5.31

Normal       $10^9$     3       $10^5$     $0.893\pm 0.118$           $7631\pm 2$              $0.902\pm 0.114$           $2234\pm 429$          1%            3.41

OSM          $10^9$     2       $10^5$     $0.932\pm0.056$            $6752 \pm 13$            $0.943\pm 0.043$           $1670\pm48$            1.2%          4.04

NYC Taxi     $10^9$     3       $10^5$     $0.526\pm 0.029$ (AUC)     $6963\pm 69$             $0.526\pm 0.030$ (AUC)     $4535\pm 7$            0%            1.53

Ablation study between FFM(GPU) and $\text{F}^{3}$M  As much of the improved performance can be attributed to our GPU implementation, we conduct an ablation study of FFM(GPU) against $\fthreem$ and KeOps in <Ref>. 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 is an exact method, the dataset distribution has no effect on computational time. Here, the smoothness criteria and adaptive far-field technique improve computational time. We find that $\fthreem$ achieves a speed-up between $2.0-33.3\times$ against FFM(GPU) and $8.0-8500\times$ speed-up against KeOps.

Comparison between $\fthreem$, FFM(GPU) and KeOps. It should be noted that KeOps is only run up to $n=10^8$ for all experiments (a run for $n=10^9$ would take weeks). The times for $n>10^8$ are extrapolated for KeOps. FFM(GPU) could only run on uniform data for $n=10^9$.

3c|$\text{F}^3$M time (s)                                                                                3c|FFM(GPU) time (s)                                                    1c|KeOps time (s)                    1c|c]@c@Speedup
vs KeOps     3cc]@c@Speedup vs

OSM                               Taxi                              $D=3$                                                                   OSM                           Taxi                            $D=3$                           $D=3$                                                                                OSM              Taxi              $D=3$

n                                                                    ($D=2$)                           ($D=3$)                                                                                ($D=2$)                       ($D=3$)                                                                                ($D=2$)          ($D=3$)

$10^6$                                                               $\makecell{0.1 \pm 0.0}$          $\makecell{0.5 \pm 0.2}$          $\makecell{0.2 \pm 0.1}$                                                $\makecell{0.7 \pm 0.2}$      $\makecell{1.2 \pm 0.2}$        $\makecell{0.4 \pm 0.2}$        1.56                                                   8.0                                                                                7.0     2.4      2.0

$10^7$                                                               $\makecell{1.1 \pm 0.4}$          $\makecell{3.3 \pm 0.6}$          $\makecell{1.7 \pm 0.7}$                                                $\makecell{10.1 \pm 0.8}$     $\makecell{110.0 \pm 1.9}$      $\makecell{3.7 \pm 1.0}$        148.6                                                  87.0                                                                                9.2     33.3     2.2

$10^8$                                                               $\makecell{10.8 \pm 3.9}$         $\makecell{30.4 \pm 6.9}$         $\makecell{17.1 \pm 7.3}$                                               OOM                           $\makecell{769.5 \pm 27.9}$     $\makecell{43.7 \pm 12.2}$      1.49e4                                                 873.0                                                                                N/A              25.3     2.6

$10^9$                                                               $\makecell{101.8 \pm 36.6}$     $\makecell{290.0 \pm 69.4}$     $\makecell{174.0 \pm 76.6}$                                           OOM                           OOM                             $\makecell{517.4 \pm 60.8}$     1.49e6                                                 8583.0                                                                                N/A              N/A               3.0

Error                                                                0.0005                            0.0012                            0.0008                                                                  0.0005                        0.0002                          0.0003                          0.0

Complexity     3c|$\mathcal{O}\left(n\cdot \log_2\left(\frac{D \cdot \mathcal{E}^2}{\gamma^2 \cdot 4 \cdot \eta} \right)\right)$     3c|$\mathcal{O}\left(n \log{(n)}\right)$                                1c|$\mathcal{O}\left(n^2\right)$     1l

Ablation study between $\fhalfm$ and $\text{F}^{3}$M  We provide and additional ablation study between $\fhalfm$ and $\fthreem$ in <Ref>. The results are quite similar to the comparison between $\fthreem$ and FFM(GPU). Here, we see that smooth field and adaptive far-field approximation ($\fthreem$) both improve speed and also memory usage as smooth field helps approximate more interactions. We can thus infer empirically that the $m_i^{\text{smooth}}$ term in <Ref> has a significant impact on reducing the memory footprint of interactions.

Comparison between $\fthreem$, $\fhalfm$ and KeOps. It should be noted that KeOps is only run up to $n=10^8$ for all experiments (a run for $n=10^9$ would take weeks). The times for $n>10^8$ are extrapolated for KeOps.

3c|$\text{F}^3$M time (s)                                                                                3c|$\fhalfm$ time (s)                                          1c|KeOps time (s)           1c|c]@c@Speedup
vs KeOps     3cc]@c@Speedup vs

OSM                             Taxi                            $D=3$                                                                       OSM                          Taxi                            $D=3$                           $D=3$                                                                                OSM               Taxi                      $D=3$

n                                                                    ($D=2$)                         ($D=3$)                                                                                ($D=2$)                      ($D=3$)                                                                                ($D=2$)           ($D=3$)

$10^6$                                                               $\makecell{0.1 \pm 0.0}$        $\makecell{0.5 \pm 0.2}$        $\makecell{0.2 \pm 0.1}$                                                    $\makecell{1.4\pm 0.3}$      $\makecell{2.0\pm 0.3}$         $\makecell{0.3 \pm 0.2}$        1.56                                                   8.0                                                                                14.0     4.0              1.5

$10^7$                                                               $\makecell{1.1 \pm 0.4}$        $\makecell{3.3 \pm 0.6}$        $\makecell{1.7 \pm 0.7}$                                                    $\makecell{11.3\pm 0.4}$     $\makecell{110.8\pm 2.0}$       $\makecell{4.5 \pm 2.1}$        148.6                                                  87.0                                                                                10.3     33.6             2.6

$10^8$                                                               $\makecell{10.8 \pm 3.9}$       $\makecell{30.4 \pm 6.9}$       $\makecell{17.1 \pm 7.3}$                                                   OOM                          $\makecell{777.4 \pm 21.5}$     $\makecell{42.5 \pm 13.1}$      1.49e4                                                 873.0                                                                                N/A               25.6             2.5

$10^9$                                                               $\makecell{101.8 \pm 36.6}$     $\makecell{290.0 \pm 69.4}$     $\makecell{174.0 \pm 76.6}$                                                 OOM                          OOM                             $\makecell{488.4 \pm 55.3}$     1.49e6                                                 8583.0                                                                                N/A               N/A                       2.8

Error                                                                0.0005                          0.0012                          0.0008                                                                      0.0004                       0.0002                          0.0016                          0.0

Complexity     3c|$\mathcal{O}\left(n\cdot \log_2\left(\frac{D \cdot \mathcal{E}^2}{\gamma^2 \cdot 4 \cdot \eta} \right)\right)$     3c|$\mathcal{O}\left(n \log{(n)}\right)$                               1c|$\mathcal{O}\left(n^2\right)$

Gaussian process regression experiments We further compare $\fthreem$ as a drop-in KMVM operation applied to Black-box Matrix Multiplication <cit.> for Gaussian Processes, compared to KISS-GP <cit.>, an approximate Gaussian process using cubic interpolation for kernel approximation. We mimic the setup in <cit.> 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\leq 3$. We demonstrate the results in <Ref>. As exact GP using $\fthreem$ demonstrates competitive results even when compared to SVGP <cit.> and SGPR <cit.>, we hypothesize that many high-dimensional datasets conform to the manifold hypothesis <cit.>, allowing $\fthreem$ to be widely applicable out-of-the-box even in high-dimensional settings.

Gaussian process regression results. Exact GP using $\fthreem$ shows improved results and scaling compared to KISS-GP. SGPR and KISS-GP could not scale to the HouseElectric dataset.

4c|RMSE                                                                                4cTraining time (s)

Dataset     $n$           $d$     1cc]@c@Exact GP
($\text{F}^3$M)     KISS-GP     c]@l@SGPR
($m=512$)     c]@l@SVGP
($m=1024$)     1cc]@c@Exact GP
($\text{F}^3$M)     KISS-GP      c]@l@SGPR
($m=512$)     c]@l@SVGP

3DRoad               278,319       3       0.297 $\pm$0.036                                                                                0.314 $\pm$0.01      0.661 $\pm$ 0.010                                                              1c|0.481 $\pm$ 0.002                                          27.8 $\pm$18.0                                                                                312.9 $\pm$10.8       720.5 $\pm$ 330.4                                                              2045.1 $\pm$ 191.4

Song                 329,820       90      0.369 $\pm$0.029                                                                                0.57 $\pm$0.298      0.803 $\pm$ 0.002                                                              0.998 $\pm$ 0.000                                                               7.2 $\pm$3.1                                                                                1705.2 $\pm$115.6     473.3 $\pm$ 187.5                                                              2373.3 $\pm$ 184.9

Buzz                 373,280       77      0.967 $\pm$0.002                                                                                0.997 $\pm$0.05      0.300$\pm$ 0.004                                                      0.304 $\pm$ 0.012                                                               33.5 $\pm$9.0                                                                                542.7 $\pm$0.8        1754.8 $\pm$ 1099.6                                                            2780.8 $\pm$ 175.6

HouseEletric         1,311,539     9       0.308 $\pm$0.006                                                                                OOM                  OOM                                                                            0.084 $\pm$ 0.005                                                      79.8 $\pm$23.1                                                                                N/A                   N/A                                                                            22062.6 $\pm$ 282.0

§ LIMITATIONS AND FURTHER RESEARCH

This work has introduced and implemented $\fthreem$ on GPU, which enables fast KMVM for tall and skinny data up to $n=10^9$. $\fthreem$ has improved complexity which also is controllable through $\eta$, and retains linear memory.  Experiments in higher dimensions also exhibit linear complexity, however requiring more nodes for lower errors. $\fthreem$ 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, $\fthreem$ 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\leq7$. A fruitful direction would be to extend ideas in $\fthreem$ to accommodate higher-dimensional data by considering randomized partitioning <cit.>, decoupling the dependency on $D$ in geometry based partitioning. Further, an exact characterization of how $m_{i}^{\text{far}},m_{i}^{\text{smooth}},m_{i}^{\text{small}},m_i^{\text{empty}}$ grows is left to future work.

The authors sincerely thank Lood van Niekerk and Jean-François Ton for their helpful comments.

§ CHECKLIST

* For all authors...

* Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?

* Did you describe the limitations of your work?

* Did you discuss any potential negative societal impacts of your work?

* Have you read the ethics review guidelines and ensured that your paper conforms to them?

* If you are including theoretical results...

* Did you state the full set of assumptions of all theoretical results?

* Did you include complete proofs of all theoretical results?
See Appendix <ref>, <ref>, <ref>.

* If you ran experiments...

* Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)?

* Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)?

* Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)?

* Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)?

* If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...

* If your work uses existing assets, did you cite the creators?

* Did you mention the license of the assets?

* Did you include any new assets either in the supplemental material or as a URL?

* Did you discuss whether and how consent was obtained from people whose data you're using/curating?

* Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content?

* If you used crowdsourcing or conducted research with human subjects...

* Did you include the full text of instructions given to participants and screenshots, if applicable?

* Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable?

* Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation?

§ SYNTHETIC DATA

Note on synthetic datasets   We generated synthetic datasets of different types to measure the ability of $\fthreem$ to deal with dense or sparse data. Dense datasets were generated as independent samples with either uniform or normal distributions. Clustered datasets were generated by sampling cluster centers from a normal distribution, and then recursively sampling sub-cluster centers from a normal distribution with reduced standard deviation and centered at each cluster center, until the desired number of points is attained. Fractional Brownian Motion and Brownian Motion samples were generated as samplings of Fractional Brownian Motion paths with respective Hurst index 0.75 and 0.5. Figure <ref> shows samples of each dataset type in the 2D case.

Uniform dataset

Normal dataset

Uniform and Normal dataset

Clustered dataset

Fractional Brownian Motion

Brownian motion
2D illustrations of the synthetic datasets.

§ NOTE ON MAXIMAL VARIANCE ON AN INTERVAL

Consider a random variable $X \in \mathbb{R}$ with finite  variance with $m=\inf X$ and $M=\sup X$. Then $\text{Var}(X)\leq \frac{(M-m)^2}{4}$.

Define a function $g$ by $g(t)=\mathbb{E}\left[(X-t)^{2}\right]$. Computing the derivative $g^{\prime}$, and solving $g^{\prime}(t)=-2 \mathbb{E}[X]+2 t=0$ yields that $g$ achieves its minimum at $t=\mathbb{E}[X]$ (note that $g^{\prime \prime}>0$ ).
Now, consider the value of the function $g$ at the special point $t=\frac{M+m}{2} .$ It must be the case that
$\operatorname{Var}[X]=g(\mathbb{E}[X]) \leq g\left(\frac{M+m}{2}\right).
Evaluating yields the expression

\begin{equation*}
\begin{split}
&g\left(\frac{M+m}{2}\right)=\\&\mathbb{E}\left[\left(X-\frac{M+m}{2}\right)^{2}\right]=\frac{1}{4} \mathbb{E}\left[((X-m)+(X-M))^{2}\right]
\end{split}
\end{equation*}

Since $X-m \geq 0$ and $X-M \leq 0$, we have

((X-m)+(X-M))^{2} \leq((X-m)-(X-M))^{2}=(M-m)^{2}

implying that

\begin{equation*}
\begin{split}
&\frac{1}{4} \mathbb{E}\left[((X-m)+(X-M))^{2}\right] \leq \\&\frac{1}{4} \mathbb{E}\left[((X-m)-(X-M))^{2}\right]=\frac{(M-m)^{2}}{4}
\end{split}
\end{equation*}

\operatorname{Var}[X] \leq \frac{(M-m)^{2}}{4}

§ DETAILS ON BARYCENTRIC LAGRANGE INTERPOLATION

The barycentric lagrange interpolation is written as

$$L_{i}(t)=\frac{\frac{w_{i}}{t-s_{i}}}{\sum_{i=0}^{r} \frac{w_{i}}{t-s_{i}}},  w_{i}=\frac{1}{\prod_{j=0, j \neq i}^{r}\left(s_{i}-s_{j}\right)}, i=0, \ldots, n

where $w_i$ are known as the barycentric weights. In case of singularities, i.e. when $t=s_j$, we set $L_i(s_j)=\delta_{ij}$. In particular, <cit.> proposes Chebyshev nodes of the second kind $s_{i}=\cos \theta_{i}, \quad \theta_{i}=\frac{i \pi}{r}, \quad i=0, \ldots, r$.
This choice of nodes combined with the scale invariance property of the barycentric form makes the calculation of $w_i$ particularly easy

w_{i}=(-1)^{i} \delta_{i}, \quad \delta_{i}=\left\{\begin{array}{ll}
1 / 2, & i=0 \text { or } i=r \\
1, & i=1, \ldots, r-1
\end{array}\right.

and reduces the complexity of calculating $w_i$ from $\mathcal{O}(r^2)$ to $\mathcal{O}(r)$. The Lagrange interpolation polynomial can then be expressed as $p_{r}(t)=\sum_{i=0}^{r} \frac{\frac{w_{i}}{t-s_{i}}}{\sum_{i=0}^{r} \frac{w_{i}}{t-s_{i}}} f_{i}$.

§ SMOOTH FIELD PROOF

Consider $\mathbf{x},\mathbf{y} \in \mathcal{X} \subset \mathbb{R}^d$ such that $d(\mathbf{x},\mathbf{y}):=\frac{\Vert \mathbf{x}-\mathbf{y}\Vert^2}{2\gamma^2}=\frac{1}{2\gamma^2}\sum_{i}^d(x^{(i)}- y^{(i)})^2 \leq \eta<1$ for all $\mathbf{x},\mathbf{y}$. When interpolating $k(\mathbf{x},\mathbf{y}) = \exp\left(-d(\mathbf{x},\mathbf{y})\right)$ using bivariate Lagrange interpolation $\mathcal{L}_r(\mathbf{x},\mathbf{y}):=\mathbf{L}_{X}^{T} \cdot \mathbf{K} \cdot \left(\mathbf{L}_{Y} \cdot \mathbf{b}\right)$ with degree $r=2p$, for any $p\in \mathbb{N}_{>0}$ there exist nodes $\mathbf{s}^{\mathbf{x}},\mathbf{s}^{\mathbf{y}}$ for $\mathcal{L}_r(\mathbf{x},\mathbf{y})$ such that the pointwise interpolation error is bounded by $\mathcal{O}(\eta^{p+1})$.

Note that $k(\mathbf{x},\mathbf{y})$ is analytic in $d(\mathbf{x},\mathbf{y})$ with Taylor expansion given by
$T_p e^{-d(\mathbf{x},\mathbf{y})}=1-d(\mathbf{x},\mathbf{y})+...+\mathcal{O}\left(d(\mathbf{x},\mathbf{y})^{p+1}\right)$. With $d(\mathbf{x},\mathbf{y})\leq \eta$, it follows that $|T_p k(\mathbf{x},\mathbf{y})-k(\mathbf{x},\mathbf{y})|\leq \mathcal{O}(\eta^{p+1})$. Using triangle inequality, we have $|\mathcal{L}_r(\mathbf{x},\mathbf{y}) - k(\mathbf{x},\mathbf{y})|\leq |T_p k(\mathbf{x},\mathbf{y})-k(\mathbf{x},\mathbf{y})| + |T_p k(\mathbf{x},\mathbf{y})-\mathcal{L}_r(\mathbf{x},\mathbf{y})| \leq \eta^{p+1} + |T_p k(\mathbf{x},\mathbf{y})-\mathcal{L}_r(\mathbf{x},\mathbf{y})|$. We note that $\mathcal{L}_r(\mathbf{x},\mathbf{y})$ contains all the terms of the Taylor expansion, and we can thus choose the nodes $\mathbf{s}^{\mathbf{x}},\mathbf{s}^{\mathbf{y}}$ of $\mathcal{L}_r(\mathbf{x},\mathbf{y})$ such that $|T_p k(\mathbf{x},\mathbf{y})-\mathcal{L}_r(\mathbf{x},\mathbf{y})|=0$ as long as $r=2p$, meaning the polynomial orders are matched.

Note that the same proof strategy can be applied to any kernel $k$ that admits a Taylor expansion.

§ COMPLEXITY

A far-field interaction between two boxes containing $n_x$ and $n_y$ points respectively has time complexity $\mathcal{O}(n)$, where $n=\max(n_x,n_y)$.

Far-field interactions are calculated as

\begin{equation}
\mathbf{v} \approx \underbrace{\mathbf{L}_{X}^{T} \cdot}_{\mathcal{O}(n_x\cdot r_X)} (\underbrace{\mathbf{K} \cdot}_{\mathcal{O}(r_X\cdot r_Y) }\underbrace{\left(\mathbf{L}_{Y} \cdot \mathbf{b}\right)}_{\mathcal{O}(n_y\cdot r_Y)  }).
\end{equation}

As $r_X,r_Y$ are independent of $n_x,n_y$, the complexity becomes $\mathcal{O}(n)$. Further see <cit.> for alternative proof.

Given $n$ data points in dimension $D$, the maximum number of divisions $\text{Tree}_{\texttt{max divisions}}$ is given by

\begin{equation}
\text{Tree}_{\texttt{max divisions}} = \log_{2^D}(n).
\end{equation}

To see this, simply solve for
\frac{n}{2^{D\cdot \text{Tree}_{\texttt{max divisions}}}}=1 \Longrightarrow  \text{Tree}_{\texttt{max divisions}} = \log_{2^D}(n).
Further see <cit.> for alternative proof.

Given a KMVM with edge $\mathcal{E}$ (dependent on data $\mathcal{X},\mathcal{Y}$), lengthscale $\gamma$, effective variance limit $\eta$, $n$ data points and data dimension $D$, $\fthreem$ has time complexity $\mathcal{O}(n\cdot \log_2\left(\frac{D \cdot  \mathcal{E}^2}{\gamma^2 \cdot 4 \cdot \eta} \right))$, which can be taken as $\mathcal{O}(n\cdot \log_2\left(\frac{C}{\eta} \right))$ where $C\propto \frac{D\cdot \mathcal{E}^2}{\gamma^2}$.

Recall that near-field interactions can be smoothly interpolated when $D\cdot \frac{\mathcal{E}^2}{\gamma^2\cdot 4\cdot 2^{\texttt{tree\_depth}}}\leq \eta$. Then all interactions will be interpolated when $\texttt{tree\_depth}\geq \log_2\left(\frac{D \cdot  \mathcal{E}^2}{\gamma^2 \cdot 4 \cdot \eta} \right)$, which implies we can take $\text{Tree}_{\text{max divisions}}=\log_2\left(\frac{D \cdot  \mathcal{E}^2}{\gamma^2 \cdot 4 \cdot \eta} \right)$. Hence the complexity is $\mathcal{O}(n\cdot \log_2\left(\frac{D \cdot  \mathcal{E}^2}{\gamma^2 \cdot 4 \cdot \eta} \right))$.

The number of interactions $M_i$ against tree depth $i$ of FFM and $\fthreem$ grows as $\mathcal{O}\left(M_{i-1}2^{2\cdot D }- m_{i}^{\text{far}} \right)$ and
$$\mathcal{O}\left(M_{i-1} 2^{2\cdot D}- (m_i^{\text{empty}})^2 - m_{i}^{\text{far}}-m_{i}^{\text{smooth}}-m_{i}^{\text{small}}  \right)$$
respectively.  Here $M_{-1} = \frac{1}{2^{2D}}$ and $m_{0}^{\text{far}}=m_{0}^{\text{smooth}}=m_{0}^{\text{small}}=m_0^{\text{empty}}=0$ and $m_{i}^{\text{far}},m_{i}^{\text{smooth}},m_{i}^{\text{small}},m_i^{\text{empty}}$ denotes the number of far-field, smooth field, small field interactions and the number of empty boxes respectively at depth $i>0$. Note that for $i>0$, these are dependent on data.

We prove through induction that the recursion holds for $\fthreem$. We start with the base case $M_{0}=1$, since at depth $i=0$, we only have one box and hence only one interaction. $M_0 = M_{-1}2^{2\cdot D} -(m_0^{\text{empty}})^2 - m_{0}^{\text{far}}-m_{0}^{\text{smooth}}-m_{0}^{\text{small}} = 1$. Clearly at depth 0, there can not be any empty boxes or possible approximations. For the induction step, $M_i= M_{i-1}2^{2\cdot D} -(m_i^{\text{empty}})^2 - m_{i}^{\text{far}}-m_{i}^{\text{smooth}}-m_{i}^{\text{small}}$. To get $M_{i+1}$, each box at depth $i$ is first divided into $2^D$, hence the number of interactions grows by $2^{2\cdot D}$. At depth $i+1$, we can further remove $(m_{i+1}^{\text{empty}})^2$ interactions between empty boxes and further compute $m_{i+1}^{\text{far}}+m_{i+1}^{\text{smooth}}+m_{i+1}^{\text{small}}$ interactions. Then $M_{i+1}=(M_{i-1}2^{2\cdot D} -(m_i^{\text{empty}})^2 - m_{i}^{\text{far}}-m_{i}^{\text{smooth}}-m_{i}^{\text{small}})2^{2\cdot D} -(m_{i+1}^{\text{empty}})^2 - m_{i+1}^{\text{far}}-m_{i+1}^{\text{smooth}}-m_{i+1}^{\text{small}} = M_i 2^{2\cdot D} -(m_{i+1}^{\text{empty}})^2 - m_{i+1}^{\text{far}}-m_{i+1}^{\text{smooth}}-m_{i+1}^{\text{small}}$. Thus the base case and induction step holds which completes our proof. This proof also covers FFM, since FFM can be as a special case for $\fthreem$ without removing empty boxes, smooth field and small field computation.

§ ALGORITHM SUMMARY

We present a summary of the algorithm presented in FFM in <Ref> and the modifications $\fthreem$ does in boldface.

Datasets $\mathbf{X},\mathbf{Y},\mathbf{b}$, kernel $k$, average points threshold $\zeta$
$\textbf{v} =  k(\mathbf{X},\mathbf{Y})\cdot \mathbf{b}$
Initialize treecodes $\tau_x = T(\mathbf{X}),\quad \tau_y = T(\mathbf{Y})$
Initialize near-field interactions as $\mathcal{I}_{\text{near}}=\{0,0\}$

Initialize output $\mathbf{v}=\mathbf{0}$

$\mid\mathcal{I}\mid>0$ and ${\texttt{MaximumBoxSize}}(\tau_x)>\zeta$ and $\texttt{MaximumBoxSize}(\tau_y)>\zeta$
Divide $\tau_x, \tau_y$

Calculate interactions left $\mathcal{I}:=f(\mathcal{I}_{\text{near}})$

Partition $\mathcal{I}$ to $\{\mathcal{I}_{\text{near}},\mathcal{I}_{\text{far}},{\boldsymbol{\mathcal{I}_{\text{smooth}}}}\}$
Throw away interactions that are too far from eachother

Compute far-field interactions $\mathbf{v}\mathrel{+}=\text{FarFieldCompute}(\tau_x,\tau_y,\mathcal{I}_{\text{far}})$
Compute smooth field interactions $\mathbf{v}\mathrel{+}=\text{FarFieldCompute}(\tau_x,\tau_y,\mathcal{I}_{\text{smooth}})$

Compute remaining near-field interactions $\mathbf{v}\mathrel{+}=\text{NearFieldCompute}(\tau_x,\tau_y,\mathcal{I}_{\text{near}})$
FFM ($\fthreem$)

Instead of comparing the average box size to $\zeta$ we compare the maximum box size. When points are non-uniformly distributed, taking the maximum ensures that we don't compute near-field interactions on boxes with many points, since it will be inefficient.

§ SCALABILITY ANALYSIS

We conduct a scalability analysis over $N_{\text{GPU}}=1,2,4,8$. We parallelize the KMVM product by considering the $k(\mathbf{X},\mathbf{X})$-case and divide the work onto multiple GPUs by partitioning each subproduct of the KMVM (see for <Ref> an example when $N_{\text{GPU}}=8$). We take $\mathbf{X}$ to be Uniform and 3 dimensional. We present results in <Ref>.

Partitioning a KMVM product into 8 jobs

Since the V100 cards we use have very high throughput, we only get a performance boost when $n=10^9$.

We also use  to analyze the % of peak throughput of the V100 cards $\fthreem$ can utilize. We run  for 3 dimensional uniform data for $n=10^6,10^7,10^8,10^9$. We present our results in <Ref>.

We used the  metric in  to generate the plot. One GPU was used for this experiment.

§ IMPACT OF $\ETA$ AND $R$ ON PERFORMANCE

The performance of $\fthreem$ is tuned by choosing $\eta$ and $r$ to trade speed against accuracy. In <Ref> we plot how different choices of $\eta$ and $r$ impacts computation time for $\fthreem$ on 3D data.

We see that larger $\eta$, more aggressive  and smaller number of interpolation nodes $r$ improve speed.

§ IMPLEMENTATION OVERVIEW

We provide a skiss of how data is stored and used for $\fthreem$ in <Ref>.

Skiss of how data is stored and used on GPU.

We provide an illustration on near field computations are carried out for $\fthreem$ in <Ref>.

Illustration of how near-field interactions are computed in parallel on GPU. First data is reordered into their corresponding boxes using a permutation vector. Then each box is loaded in parallel into several thread blocks (indexed by $\texttt{blockIdx.x}$), wherein the challenge lies to execute this correctly. The computations are then parallelized across blocks, where only interactions are computed.