A Dynamic Low-Rank Fast Gaussian Transform

## Abstract

The Fast Gaussian Transform (FGT) enables subquadratic-time multiplication of an n × n 𝑛 𝑛 n\times n italic_n × italic_n Gaussian kernel matrix 𝖪 i , j = exp ⁡ ( − ‖ x i − x j ‖ 2 2 ) subscript 𝖪 𝑖 𝑗 superscript subscript norm subscript 𝑥 𝑖 subscript 𝑥 𝑗 2 2 \mathsf{K}_{i,j}=\exp(-\|x_{i}-x_{j}\|_{2}^{2}) sansserif_K start_POSTSUBSCRIPT italic_i , italic_j end_POSTSUBSCRIPT = roman_exp ( - ∥ italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_x start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) with an arbitrary vector h ∈ ℝ n ℎ superscript ℝ 𝑛 h\in\mathbb{R}^{n} italic_h ∈ blackboard_R start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT , where x 1 , … , x n ∈ ℝ d subscript 𝑥 1 … subscript 𝑥 𝑛 superscript ℝ 𝑑 x_{1},\dots,x_{n}\in\mathbb{R}^{d} italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , … , italic_x start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ∈ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT are a set of fixed source points. This kernel plays a central role in machine learning and random feature maps. Nevertheless, in most modern data analysis applications, datasets are dynamically changing (yet often have low rank), and recomputing the FGT from scratch in (kernel-based) algorithms incurs a major computational overhead ( ≳ n greater-than-or-equivalent-to absent 𝑛 \gtrsim n ≳ italic_n time for a single source update ∈ ℝ d absent superscript ℝ 𝑑 \in\mathbb{R}^{d} ∈ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT ). These applications motivate a dynamic FGT algorithm, which maintains a dynamic set of sources under kernel-density estimation (KDE) queries in sublinear time while retaining Mat-Vec multiplication accuracy and speed. Assuming the dynamic data-points x i subscript 𝑥 𝑖 x_{i} italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT lie in a (possibly changing) k 𝑘 k italic_k -dimensional subspace ( k ≤ d 𝑘 𝑑 k\leq d italic_k ≤ italic_d ), our main result is an efficient dynamic FGT algorithm, supporting the following operations in log O ⁢ ( k ) ⁡ ( n / ε ) superscript 𝑂 𝑘 𝑛 𝜀 \log^{O(k)}(n/\varepsilon) roman_log start_POSTSUPERSCRIPT italic_O ( italic_k ) end_POSTSUPERSCRIPT ( italic_n / italic_ε ) time: (1) Adding or deleting a source point, and (2) Estimating the “kernel-density” of a query point with respect to sources with ε 𝜀 \varepsilon italic_ε additive accuracy. The core of the algorithm is a dynamic data structure for maintaining the projected “interaction rank” between source and target boxes, decoupled into finite truncation of Taylor and Hermite expansions.

## 1 Introduction

The fast Multipole method (FMM) was described as one of the top 10 most important algorithms of the 20th century []. It is a numerical technique that was originally developed to speed up calculations of long-range forces for the $n$-body problem in theoretical physics. FMM was first introduced in 1987 by Greengard and Rokhlin [], based on the multipole expansion of the vector Helmholtz equation (see Section A). By treating the interactions between far-away basis functions using the FMM, the underlying matrix entries $M_{ij}\in{\mathbb{R}}^{n\times n}$
(encoding the pairwise “interaction” between $x_{i},x_{j}\in{\mathbb{R}}^{d}$)
need not be explicitly computed nor stored for matrix-vector operations –
This technique allows to improve the naïve $O(n^{2})$ matrix-vector multiplication time to quasi-linear time $\approx n\cdot\log^{O(d)}(n)$, with negligible (polynomial-small) additive error.

Since the discovery of FMM in the late 80s, it had a profound impact on scientific computing and has been extended and applied in many different fields, including physics, mathematics, numerical analysis and computer science [, , , , , , , , , , , , , , ]. To mention just one important example, we note that FMM plays a key role in efficiently maintaining the SVD of a matrix under low-rank perturbations, based on the Cauchy structure of the perturbed eigenvectors [].

In the context of machine learning,
the FMM technique can be extended to the evaluation of matrix-vector products with certain *Kernel matrices*
$\mathsf{K}_{i,j}=f\left(\|x_{i}-x_{j}\|\right)$,
most notably, the *Gaussian Kernel*
$\mathsf{K}_{i,j}=\exp(-\|x_{i}-x_{j}\|_{2}^{2})$
[].
For any query vector $q\in{\mathbb{R}}^{n}$, the *fast Gaussian transform* (FGT) algorithm outputs an arbitrarily-small *pointwise additive* approximation to $\mathsf{K}\cdot q$, i.e.,
a vector $z\in{\mathbb{R}}^{n}$ such that

$$ $\displaystyle\|\mathsf{K}\cdot q-z\|_{\infty}\leq\varepsilon,$ $$

in merely
$n\log^{O(d)}(\|q\|_{1}/\varepsilon)$ time, which is dramatically faster than naïve matrix-vector multiplication ($n^{2}$) for constant dimension $d$. Note that the (poly)logarithmic dependence on $1/\varepsilon$ means that FGT can achieve *polynomially-small* additive error in quasi-linear time, which is as good as exact computation for all practical purposes. The crux of FGT is that the $n\times n$ matrix $\mathsf{K}$ can be stored *implicitly*, using a clever spectral-analytic decomposition of the geometrically-decaying pairwise distances (“interaction rank”, more on this below).

Kernel matrices play a central role in machine learning [, ],
as they allow to extend convex optimization and learning algorithms to nonlinear feature spaces and even to non-convex problems [, , , , , ].
Accordingly, matrix-vector multiplication with kernel matrices is a basic operation in many ML optimization tasks, such as Kernel PCA and ridge regression [, , , ],
Gaussian-process regression (GPR) []),
Kernel linear system solvers (via Conjugate Gradient []), and in fast implementation of the dynamic “state-space model” (SSM) for sequence-correlation modeling (which crucially relies on the Multipole method []),
to mention a few.
The related data-structure problem of *kernel density estimation* of a point [, , , , , ]

$$ $\displaystyle\mathsf{KDE}(X,y)=\frac{1}{n}\sum_{i=1}^{n}\mathsf{K}(x_{i},y)$ $$

has various applications in data analysis and statistics
[, , ], and is the main subroutine in the implementation of transfer learning using kernels (see [, ] and references therein, and the Related Work section below).
As such, speeding up matrix-vector multiplication with kernel matrices, such as FGT,
is an important question in theory and practice.

One drawback of FMM and FGT techniques, however, is that they are *static* algorithms, i.e., they assume a *fixed* set of $n$ data points $x_{i}\in{\mathbb{R}}^{d}$. By contrast, most aforementioned ML and data analysis applications are *dynamic* by nature and need to process rapidly-evolving datasets to maintain prediction and model accuracy.
One example is the renewed interest in *online regression* [, ], motivated by *continual learning* theory [].
Indeed, it is becoming increasingly clear that many static optimization algorithms do not capture the requirements of real-world applications [, , , , , , ].

Notice that changing a single source-point $x_{i}\in{\mathbb{R}}^{d}$ generally affects an *entire row* ($n$ distances $\|x_{i}-x_{j}\|$) of the matrix $\mathsf{K}$. As such, naively re-computing the static FGT on the modified set of distances, incurs a prohibitive computational overhead ($n\gg d$). This raises the natural question of whether it is possible to achieve *sublinear*-time insertion and deletion of source points, as well as “local” *kernel-density estimation* (KDE) queries [, ], while maintaining speed and accuracy of matrix-vector multiplication queries:

Is it possible to ‘dynamize’ the Fast Gaussian Transform, in sublinear time? Can the exponential dependence on $d$ [] be mitigated if the data-points
$x_{i}$ lie in a $k$-dimensional subspace of ${\mathbb{R}}^{d}$?

The last question is motivated by the recent work of [], who observed that kernel-based methods and algorithms typically involve *low-rank* datasets,
(where the “intrinsic” dimension is $w\ll d$), in which
case one could hope to circumvent the exponential dependence on $d$ in the aforementioned (static) FMM algorithm [, ].

### 1.1 Main Result

Our main result is an affirmative answer to the above question. We design a fully-dynamic FGT data structure, supporting *polylogarithmic*-time updates and “density estimation” queries, while retaining quasi-linear time for arbitrary Mat-Vec queries ($\mathsf{K}q$).

More formally, for a set of $N$ “source” points $s_{1},\ldots,s_{N}$, the $j$-th coordinate $(\mathsf{K}q)_{j\in[N]}$ is $G(s_{j})=\sum_{i=1}^{N}q_{i}\cdot e^{-\|s_{j}-s_{i}\|_{2}^{2}/\delta}$, which measures the kernel-density at $s_{j}$ (“interaction” of $s_{j}$ with the rest of the sources). More generally, for any “target” point $t\in{\mathbb{R}}^{d}$, let

$$ $G(t):=\sum_{i=1}^{N}q_{i}\cdot e^{-\|t-s_{i}\|_{2}^{2}/\delta}$ $$

denote the *kernel density* of $t$ with respect to the sources, where each source $s_{i}$ is equipped with a *charge* $q_{i}$. Our data structure supports fully-dynamic source updates and density-estimation queries in *sublinear* time.
Observe that this immediately implies that
entire Mat-Vec queries ($\mathsf{K}\cdot q$) can be computed in quasi-linear time $N^{1+o(1)}$.
The following is our main result:

###### Theorem 1.1 (Dynamic Low-Rank FGT, Informal version of Theorem F.2 ) .

Let $\mathfrak{B}$ denote a $w$-dimensional subspace $\subset{\mathbb{R}}^{d}$. Given a set of source points $s$, and charges $q$, there is a (deterministic) data structure that maintains a fully-dynamic set of $N$ source vectors $s_{1},\cdots,s_{N}\in\mathfrak{B}$
under the following operations:

- •
Insert/Delete$(s_{i}\in{\mathbb{R}}^{d},q_{i}\in{\mathbb{R}})$ Insert or Delete a source point $s_{i}\in{\mathbb{R}}^{d}$ along with its “charge” $q_{i}\in{\mathbb{R}}$, in $\log^{O(w)}(\|q\|_{1}/\varepsilon)$ time. The intrinsic subspace $B$ could change as the source points are updated.
- •
Density-Estimation$(t\in\mathfrak{B})$ For any point $t\in\mathfrak{B}\subset{\mathbb{R}}^{d}$, output the kernel density of $t$ with respect to the sources, i.e., output $\widetilde{G}$ such that $G(t)-\varepsilon\leq\widetilde{G}\leq G(t)+\varepsilon$ in $\log^{O(w)}(\|q\|_{1}/\varepsilon)$ time.

We note that when $w=d$, the costs of our dynamic algorithm match the statistic FGT algorithm.

As one might expect, our data structure applies to a more general subclass of ‘geometrically-decaying’ kernels $\mathsf{K}_{i,j}=f(\|x_{i}-x_{j}\|)$ ($f(tx)\leq(1-\alpha)^{t}f(x)$), see Theorem B.5 for the formal statement of our main result.
It is also noteworthy that our data structure is deterministic, and therefore handles even *adaptive* update sequences [, , ]. This feature is important in adaptive data analysis and in the use of dynamic data structures for accelerating *path-following* iterative optimization algorithms [], where proximity to the original gradient flow (linear) equations is crucial for convergence, hence the data structure needs to ensure the approximation guarantees hold against *any* outcome of previous iterations.

###### Theorem 1.1 (Dynamic Low-Rank FGT, Informal version of Theorem F.2 ) .

#### Remark on Dynamization of
“Decomposable” Problems

A data structure problem $\mathbf{P}(D,q)$ is called *decomposable*, if a query $q$ to the *union* of two separate datasets can be recovered from the two *marginal* answers of the query on each of them separately, i.e.,

$$ $\displaystyle\mathbf{P}(D_{1}\cup D_{2},q)=g\left(\mathbf{P}(D_{1},q),\mathbf{ P}(D_{2},q)\right)$ $$

for some function $g$.
A classic technique in data structures [] asserts that decomposable data structure problems can be (partially) dynamized in a *black-box* fashion – It is possible to convert any *static* DS for $\mathbf{P}$ into a dynamic one supporting incremental updates, with an amortized update time
$t_{u}\sim(T/N)\cdot\log(N)$, where
$T$ is the preprocessing
time of building the static data structure, and $N$ is the input size.
It is not hard to see that Matrix-Vector multiplication over a field (with row-updates to the matrix), is a decomposable problem (since $(A+B)q=Aq+Bq$), and so one might hope that the dynamization of static FMM/FGT methods is an immediate consequence of decomposability.
This reasoning is, unfortunately, incorrect, since changing even a
*single* input point $x_{i}\in{\mathbb{R}}^{d}$, perturbs $n$ distances (i.e., an entire *row* in the kernel matrix $\mathsf{K}$), and so the aforementioned reduction is prohibitively expensive (yields update time at least $n\gg d$ for adding/removing a point).

#### Notations

For a vector $x$, we use $\|x\|_{2}$ to denote its $\ell_{2}$-norm, $\|x\|_{1}$, $\|x\|_{0}$ and $\|x\|_{\infty}$ for its $\ell_{1}$-norm, $\ell_{0}$-norm and $\ell_{\infty}$-norm. We use $\widetilde{O}(f)$ to denote $f\cdot\mathrm{poly}(\log f)$. For a vector $x\in{\mathbb{R}}^{d}$ and a real number $p$, we say $x\leq p$ if $x_{i}\leq p$ for all $i\in[d]$. We say $x\geq p$ if there exists an $i\in[d]$ such that $x_{i}\geq p$. For a positive integer $n$, we use $[n]$ to denote a set $\{1,2,\cdots,n\}$.

### 1.2 Related Work

#### Structured Linear Algebra

Multiplying an $n\times n$ matrix $M$ by an arbitrary vector $q\in{\mathbb{R}}^{n}$ generally requires $\Theta(n^{2})$ time, and this is information-theoretically optimal since merely reading the entries of the matrix requires $\sim n^{2}$ operations. Nevertheless,
if $M$ has some *structure* ($\widetilde{O}(n)$-bit description-size), one could hope for quasi-linear time for computing $M\cdot q$. Kernel matrices ($\mathsf{K}_{ij}=f(\|x_{i}-x_{j}\|)$), which are the subject of this paper, are special cases of such *geometric-analytic* structure, as their $n^{2}$ entries are determined by only $\sim n$ points in ${\mathbb{R}}^{d}$, i.e., $O(nd)$ bits of information.
There is a rich and active body of work in *structured linear algebra*, exploring various “algebraic” structures that allow quasi-linear time matrix-vector multiplication, most of which relies on (novel) extensions of the *Fast Fourier Transform* (see [, , ] and references therein).

A key difference between FMMs and the aforementioned FFT-style line of work is that the latter develops *exact* Mat-Vec algorithms, whereas FMM techniques must inevitably resort to (small) approximation, based on the *analytic* smoothness properties of the underlying function and metric space [, ]. This distinction makes the two lines of work mostly incomparable.

#### Comparison to LSH-based KDEs

A recent line of work due to [, , , , ]
develops fast KDE data structures based on *locality-sensitive hashing* (LSH), which seems possible to be dynamized naturally (as LSH is dynamic by nature).
However, this line of work is incomparable to FGT, as it solves KDE in the *low-accuracy* regime, i.e., the runtime dependence on $\varepsilon$ of these works is $\mathrm{poly}(1/\varepsilon)$ (but polynomial in $d$), as opposed to FGT ($\mathrm{poly}\log(1/\varepsilon)$ but exponential in $d$). Additionally, some work (e.g., []) also needs an upper bound of the ground-truth value $\mu_{\star}=\mathsf{K}\cdot q$, and the efficiency of their data structure depends on $\mu_{\star}^{-O(1)}$, while FGT does not need any prior knowledge of $\mu_{\star}$.

#### Kernel Methods in ML

Kernel methods can be thought of as instance-based learners: rather than learning some fixed set of parameters corresponding to the features of their inputs, they instead “remember” the $i$-th training example $(x_{i},y_{i})$ and learn for it a corresponding weight $w_{i}$. Prediction for unlabeled inputs, i.e., those not in the training set, is treated using an application of a *similarity* function $\mathsf{K}$ (i.e., a kernel) between the unlabeled input $x^{\prime}$ and each of the training-set inputs $x_{i}$. This framework is one of the main motivations for the development of kernel methods in ML and high-dimensional statistics [].
There are two main themes of research on kernel methods in the context of machine learning: The first one is focused on understanding the expressive power and generalization of learning with kernel feature maps [, , , , , , , ]; The second line is focused on the *computational* aspects of
kernel-based algorithms [, , , , , , , , , , ]. We refer the reader to these references for a much more thorough overview of these lines of research and the role of kernels in ML.

## 2 Technical Overview

In this section, we provide a streamlined overview of our data structure (Theorem (B.5)).

In Section 2.1, we review the *offline* FGT algorithm [, ] and analyze the computational costs.
In Section 2.2, we illustrate the technique of estimating $G(t)$ for an arbitrary target vector $t\in{\mathbb{R}}^{d}$.
In Section 2.3, we explain that the data structures support the dynamic setting where the source vectors are allowed to come and leave.
In Section 2.4, we describe how to extend the data structure to a more general kernel function.
In Section 2.5, we show that if the source and target vectors come from a low dimensional subspace, the data structure can bypass the curse of dimension.
In Section 2.6, we modify the data structure to support the scenario where the rank of data points varies across iterations.

Figure: Figure 1: An illustration of the source-target boxing our data structure maintains in high dimensional space, using the “hybrid” of Taylor-Hermite expansions.
Refer to caption: x1.png

### 2.1 Offline FGT Algorithm

We first review []’s offline FGT algorithm.
Consider the following easier problem: given $N$ source vectors $s_{1},\dots,s_{N}\in{\mathbb{R}}^{d}$, and $M$ target vectors $t_{1},\dots,t_{M}\in{\mathbb{R}}^{d}$, estimate

$$ $\displaystyle G(t_{i})=\sum_{j=1}^{N}q_{j}\cdot e^{-\|t_{i}-s_{j}\|_{2}^{2}/\delta}$ $$

for any $i\in[M]$, in quasi-linear time.

Following [, ], our algorithm subdivides $B_{0}=[0,1]^{d}$
into smaller boxes with sides of length $L=r\sqrt{2\delta}$ parallel to the axes, for a fixed $r\leq 1/2$, and then assign each source $s_{j}$ to the box ${\cal B}$ in which it lies and each target $t_{i}$ to the box ${\cal C}$ in which it lies. Note that there are $(1/L)^{d}$ boxes in total. Let $N(B)$ and $N(C)$ denote the number of non-empty source and target boxes, respectively.
For each target box ${\cal C}$, we need to evaluate the total field due to sources in all boxes. Since each box ${\cal B}$ has side length $r\sqrt{2\delta}$, only a fixed number of source boxes ${\cal B}$ can contribute more than $\|q\|_{1}\varepsilon$ to the field in a given target box ${\cal C}$, where $\varepsilon$ is the precision parameter.
Hence, for a target vector in box ${\cal C}$, if we only count the contributions of the source vectors in its $(2k+1)^{d}$ nearest boxes where $k$ is a parameter, it will incur an error that can be upper bounded as follows:

$$ $\displaystyle\sum_{j:\|t-s_{j}\|_{\infty}\geq kr\sqrt{2\delta}}|q_{j}|\cdot e^ {-\|t-s_{j}\|_{2}^{2}/\delta}\leq\|q\|_{1}\cdot e^{-2r^{2}k^{2}}$ (1) $$

When we take $k=\log(\|q\|_{1}/\varepsilon)$, this error becomes $o(\varepsilon)$.

For a single source vector $s_{j}\in{\cal B}$, its field $G_{s_{j}}(t)=q_{j}\cdot e^{-\|t-s_{j}\|^{2}/\delta}$ has the following Taylor expansion at $t_{\cal C}$ (the center of ${\cal C}$):

$$ $\displaystyle G_{s_{j}}(t)=\sum_{\beta\geq 0}{\cal B}_{\beta}(j,{\cal C})\left (\frac{t-t_{\cal C}}{\sqrt{\delta}}\right)^{\beta},$ (2) $$

where $\beta\in\mathbb{N}^{d}$ is a multi-index,

$$ $\displaystyle{\cal B}_{\beta}(j,{\cal C})=q_{j}\cdot\frac{(-1)^{\|\beta\|_{1}} }{\beta!}\cdot H_{\beta}\left(\frac{s_{j}-t_{\cal C}}{\sqrt{\delta}}\right),$ $$

and $H_{\beta}(x)$ is the multi-dimensional Hermite function indexed by $\beta$ (Definition A.7).
We can also control the truncation error of the first $p^{d}$ terms by $\varepsilon$ for $p=\log(\|q\|_{1}/\varepsilon)$ (see Lemma E.6).
Then, for a fixed source box ${\cal B}$, the field can be approximated by

$$ $\displaystyle\sum_{\beta\leq p}C_{\beta}({\cal B},{\cal C})\left(\frac{t-t_{ \cal C}}{\sqrt{\delta}}\right)^{\beta},$ $$

where $C_{\beta}({\cal B},{\cal C}):=\sum_{j\in{\cal B}}{\cal B}_{\beta}(j,{\cal C})$.
Hence, for each query point $t$, we just need to locate its target box ${\cal C}$, and then $G(t)$ can be approximated by:

$$ $\displaystyle\widetilde{G}(t)=$ $\displaystyle\leavevmode\nobreak\ \sum_{{\cal B}\in\mathsf{nb}({\cal C})}\sum_ {\beta\leq p}C_{\beta}({\cal B},{\cal C})\left(\frac{t-t_{\cal C}}{\sqrt{ \delta}}\right)^{\beta}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \sum_{\beta\leq p}C_{\beta}({\cal C})\left( \frac{t-t_{\cal C}}{\sqrt{\delta}}\right)^{\beta},$ $$

where $\mathsf{nb}({\cal C})$ is the set of $(2k+1)^{d}$ nearest-neighbor of ${\cal C}$ and

$$ $\displaystyle C_{\beta}({\cal C}):=\sum_{{\cal B}\in\mathsf{nb}({\cal C})}C_{ \beta}({\cal B},{\cal C}).$ $$

Notice that we can further pre-compute $C_{\beta}({\cal C})$ for each target box ${\cal C}$ and $\beta\leq p$. Then, the running time for each target point becomes $O(p^{d})$.

For the preprocessing time, notice that each $C_{\beta}({\cal B},{\cal C})$ takes $O(N_{\cal B})$-time to compute, where $N_{\cal B}$ is the number of source points in ${\cal B}$. Fix a $\beta\leq p$.
Consider the computational cost of $C_{\beta}({\cal C})$ for all target boxes ${\cal C}$. Note that each source box can interact with at most $(2k+1)^{d}$ target boxes.
Therefore, the total running time for computing $\{C_{\beta}({\cal C}_{\ell})\}_{\ell\in[N(C)]}$ is bounded by

$$ $O\left(N\cdot(2k+1)^{d}+M\right).$ $$

Then, the total cost of the preprocessing is

$$ $\displaystyle O\left(N\cdot(2k+1)^{d}\cdot p^{d}+M\cdot p^{d}\right).$ $$

By taking $p=\log(\|q\|_{1}/\varepsilon)$ and $k\leq\log(\|q\|_{1}/\varepsilon)$, we get an algorithm with $\widetilde{O}_{d}(N+M)$-time for preprocessing and $\widetilde{O}_{d}(1)$-time for each target point.

We note that this algorithm also supports fast computing $\mathsf{K}q$ for any $q\in{\mathbb{R}}^{d}$ and $\mathsf{K}\in{\mathbb{R}}^{n\times n}$ with $\mathsf{K}_{i,j}=e^{-\|s_{i}-s_{j}\|_{2}^{2}/\delta}$. Roughly speaking, for each query vector $q$, we can build this data structure, and then the $i$-th coordinate of $\mathsf{K}q$ is just $G(s_{i})$, which can be computed in poly-logarithmic time. Hence, $\mathsf{K}q$ can be approximately computed in nearly-linear time with $\ell_{\infty}$ error at most $\varepsilon$.

### 2.2 Online Static KDE Data Structure (Query-Only)

Next, we consider the same static setting, except target queries $t\in{\mathbb{R}}^{d}$ arrive online, and the goal is to estimate $G(t)$ for an arbitrary vector in *sublinear* time.
To this end, note that
if $t$ is contained in a non-empty target box ${\cal C}_{\ell}$, then $G(t)$ can be approximated using pre-computed $C_{\beta}({\cal C}_{\ell})$ in poly-logarithmic time. Otherwise, we need to add a new target box ${\cal C}_{N(C)+1}$ for $t$ and compute $C_{\beta}({\cal C}_{N(C)+1})$, which takes time

$$ $\displaystyle\sum_{{\cal B}\in\mathsf{nb}({\cal C}_{N(C)+1})}O(N_{{\cal B}}).$ $$

Note, however, that this linear scan naïvely takes $O(N)$ time in the worst case.
Indeed, looking into the coefficients $C_{\beta}({\cal B},{\cal C})$:

$$ $\displaystyle C_{\beta}({\cal B},{\cal C})=\sum_{j\in{\cal B}}q_{j}\cdot\frac{ (-1)^{\|\beta\|_{1}}}{\beta!}\cdot H_{\beta}\left(\frac{s_{j}-t_{\cal C}}{ \sqrt{\delta}}\right)$ $$

reveals that the source vectors $s_{j}$ are “entangled” with $t_{\cal C}$, so evaluating $C_{\beta}({\cal B},{\cal C})$ brute-forcely for a new target box ${\cal C}$, incurs a linear scan of all source vectors in ${\cal B}$.

To “disentangle” $s_{j}$ and $t_{\cal C}$, we use the Taylor series of Hermite function (Eq. (5)):

$$ $\displaystyle\leavevmode\nobreak\ H_{\beta}\left(\frac{s_{j}-t_{\cal C}}{\sqrt {\delta}}\right)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ H_{\beta}\left(\frac{s_{j}-s_{\cal B}}{\sqrt {\delta}}+\frac{s_{\cal B}-t_{\cal C}}{\sqrt{\delta}}\right)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \sum_{\alpha\geq 0}\frac{(-1)^{\|\alpha\|_{1 }}}{\alpha!}\left(\frac{s_{j}-s_{\cal B}}{\sqrt{\delta}}\right)^{\alpha}H_{ \alpha+\beta}\left(\frac{s_{\cal B}-t_{\cal C}}{\sqrt{\delta}}\right),$ $$

where $s_{\cal B}$ denotes the center of the source box ${\cal B}$.

Hence, $C_{\beta}({\cal B},{\cal C})$ can be re-written as:

$$ $\displaystyle\leavevmode\nobreak\ C_{\beta}({\cal B},{\cal C})$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \sum_{j\in{\cal B}}q_{j}g(\beta)\sum_{\alpha \geq 0}g(\alpha)\left(\frac{s_{j}-s_{\cal B}}{\sqrt{\delta}}\right)^{\alpha}H_ {\alpha+\beta}\left(\frac{s_{\cal B}-t_{\cal C}}{\sqrt{\delta}}\right)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ g(\beta)\sum_{\alpha\geq 0}A_{\alpha}({\cal B })H_{\alpha+\beta}\left(\frac{s_{\cal B}-t_{\cal C}}{\sqrt{\delta}}\right),$ $$

where $g(x)=(-1)^{\|x\|_{1}}/x!$ and

$$ $\displaystyle A_{\alpha}({\cal B}):=\sum_{j\in{\cal B}}q_{j}g(\alpha)\left( \frac{s_{j}-s_{\cal B}}{\sqrt{\delta}}\right)^{\alpha}.$ (3) $$

Now, $A_{\alpha}({\cal B})$ does not rely on the target box and can be pre-computed, hence we can compute $C_{\beta}({\cal B},{\cal C})$ without going over each source vector. However, there is a price for this conversion, namely, that now $C_{\beta}({\cal B},{\cal C})$ involves summing over all $\alpha\geq 0$, so we need to somehow truncate this series while controlling the overall truncation error for $G(t)$, which appears difficult to achieve.

To this end, we observe that this two-step approximation is equivalent to first forming a truncated Hermite series of $e^{\|t-s_{j}\|^{2}/\delta}$ at the center of the source box $s_{\cal B}$, and then transforming all Hermite expansions into Taylor expansions at the center of a *target* box $t_{\cal C}$. More formally, the Hermite approximation of $G(t)$ is

$$ $\displaystyle\leavevmode\nobreak\ G(t)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \sum_{{\cal B}}\sum_{\alpha\leq p}(-1)^{\| \alpha\|_{1}}A_{\alpha}({\cal B})H_{\alpha}\left(\frac{t-s_{\cal B}}{\sqrt{ \delta}}\right)+\mathrm{Err}_{H}(p),$ $$

where $|\mathrm{Err}_{H}(p)|\leq\varepsilon$ (see Lemma E.2). Hence, we can Taylor-expand each $H_{\alpha}$ at $t_{\cal C}$ and get that:

$$ $\displaystyle G(t)=$ $\displaystyle\leavevmode\nobreak\ \sum_{\beta\leq p}C_{\beta}({\cal C})\left( \frac{t-t_{\cal C}}{\sqrt{\delta}}\right)^{\beta}+\mathrm{Err}_{T}(p)+\mathrm{ Err}_{H}(p),$ $$

where

$$ $\displaystyle|\mathrm{Err}_{H}(p)|+|\mathrm{Err}_{T}(p)|\leq\varepsilon$ $$

(for the formal argument, see Lemma E.5).

###### Remark 2.1 .

The original FGT paper contains a flaw in the error estimation, which was partially fixed in [] for the Hermite expansion. Later, [] corrected the error in both Hermite and Taylor expansions. However, their proofs are brief and use different notations that are adapted for their dual-tree algorithm. We provide more detailed and user-friendly proofs for the correct error estimations in Section E. We believe that they are of independent interest to the community.

This means that, at preprocessing time, it suffices to compute $A_{\alpha}({\cal B})$ for all source boxes and all $\alpha\leq p$, which takes

$$ $\displaystyle\sum_{k\in[N(B)]}O\left(p^{d}\cdot N_{{\cal B}_{k}}\right)= \leavevmode\nobreak\ O\left(p^{d}\cdot N\right)=\leavevmode\nobreak\ \widetilde{O}_{d}(N).$ $$

time. Then, at query time, given an arbitrary query vector $t$ in a target box ${\cal C}$, we compute

$$ $\displaystyle C_{\beta}({\cal C})=h(\beta)\sum_{{\cal B}\in\mathsf{nb}({\cal C })}\sum_{\alpha\leq p}A_{\alpha}({\cal B})H_{\alpha+\beta}\left(\frac{s_{{\cal B }}-t_{{\cal C}}}{\sqrt{\delta}}\right),$ $$

which takes
$O\left(d\cdot p^{d}\cdot(2k+1)^{d}\right)=\mathrm{poly}\log(n)$
time, so long as $d=O(1)$ and $\varepsilon=n^{-O(1)}$.

###### Remark 2.1 .

### 2.3 Dynamization

Given our (static) representation of points from the last paragraph, dynamizing the above static KDE data structure now becomes simple.
Suppose we add a source vector $s$ in the source box ${\cal B}$. We first update the intermediate variables $A_{\alpha}(\mathcal{B}),\alpha\leq p$, which takes $O(p^{d})$ time.
So long as the $\ell_{1}$-norm of the updated charge-vector $q$ remains polynomial in the norm of the previously maintained vector, namely

$$ $\displaystyle\sqrt{\log(\|q^{\mathrm{new}}\|_{1})}>\log(\|q\|_{1}),$ $$

we show that one source box can only affect $(2k+1)^{d}$ nearest target box ${\cal C}$; otherwise, when the change is super-polynomial, we rebuild the data structure, but this cost is amortized away. Hence, we only need to update $C_{\beta}({\cal C})$ for those ${\cal C}\in\mathsf{nb}({\cal B})$. Notice that each $C_{\beta}({\cal B},{\cal C})$ can be updated in $O_{d}(1)$ time, so each affected $C_{\beta}({\cal C})$ can also be updated in $O_{d}(1)$ time. Hence, adding a source vector can be done in time
$O((2k+1)^{d}p^{d})=\widetilde{O}_{d}(1)$
as before.
Deleting a source vector follows from a similar procedure.

### 2.4 Generalization to Fast-Decaying Kernels

We briefly explain how the dynamic FGT data structure generalizes to more general kernel functions $\mathsf{K}(s,t)=f(\|s-t\|_{2})$ where $f$ satisfies the 3 properties in Definition B.1 (cf. []).

In the general case,

$$ $\displaystyle G_{f}(t)=\sum_{{\cal B}}\sum_{j\in{\cal B}}q_{j}\mathsf{K}(s_{j} ,t).$ $$

Similar to the Gaussian kernel case, we can first show that only near boxes matter:

$$ $\displaystyle\sum_{j:\|t-s_{j}\|_{\infty}\geq kr}|q_{j}|\cdot f(\|s-t\|_{2})\leq\varepsilon$ $$

by the fast-decreasing property
(P2) in Definition B.1 of $f$ and taking $k=O(\log(\|q\|_{1}/\varepsilon))$.(^1^11Indeed, by property P2 , $f(\Theta(\log(1/\varepsilon^{\prime})))\leq\varepsilon^{\prime}$. Taking $\varepsilon^{\prime}:=\varepsilon/\|q\|_{1}$, we get that $f(\|s-t\|_{2})\leq\varepsilon/\|q\|_{1}$. Hence, the summation is at most $\varepsilon$.)
Then, we can follow the same “decoupling” approach as the Gaussian kernel case to first Hermite expand $G_{f}(t)$ at the center of each source box and then Taylor expands each Hermite function at the center of the target box. In this way, we can show that

$$ $\displaystyle G_{f}(t)\approx\sum_{\beta\leq p}C_{f,\beta}({\cal C})\left( \frac{t-t_{\cal C}}{\sqrt{\delta}}\right)^{\beta}$ $$

where

$$ $\displaystyle C_{f,\beta}({\cal C})=c_{\beta}\sum_{{\cal B}\in\mathsf{nb}({ \cal C})}\sum_{\alpha\leq p}A_{f,\alpha}({\cal B})H_{\alpha+\beta}\left(\frac{ s_{{\cal B}}-t_{{\cal C}}}{\sqrt{\delta}}\right),$ $$

and the approximation error can be bounded since $f$ is truncateable.
$A_{f,\alpha}({\cal B})$ depends on the kernel function $f$ and can be pre-computed in the preprocessing. Then, each $C_{f,\beta}({\cal C})$ can be computed in poly-logarithmic time. Hence, $G(t)$ can be approximately computed in poly-logarithmic time for any target vector $t$.

Figure: Figure 2: An illustration of inserting two source points with corresponding interactions to the data structure.
Refer to caption: x2.png

### 2.5 Handling Points From Low-Dimensional Static Spaces

In many practical problems, the data lies in a low dimensional subspace of ${\mathbb{R}}^{d}$. We can first project the data into this subspace and then perform FGT on ${\mathbb{R}}^{w}$, where $w$ is the rank. The following lemma shows that FGT can be performed on the projections of the data.

###### Lemma 2.2 (Hermite projection lemma in low-dimensional space, informal version of Lemma F.3 ) .

Given $B\in{\mathbb{R}}^{d\times w}$ that defines a $w$-dimensional subspace of ${\mathbb{R}}^{d}$, let

$$ $\displaystyle B^{\top}B=U\Lambda U^{\top}\in{\mathbb{R}}^{w\times w}$ $$

denote the spectral decomposition where $U\in{\mathbb{R}}^{w\times w}$ and a diagonal matrix $\Lambda\in{\mathbb{R}}^{w\times w}$.

We define

$$ $\displaystyle{\mathsf{P}}:=\Lambda^{-1/2}U^{-1}B^{\top}\in{\mathbb{R}}^{w \times d}.$ $$

Then we have for any $t,s\in{\mathbb{R}}^{d}$ from subspace $B$, the following equation holds

$$ $\displaystyle\leavevmode\nobreak\ e^{-\|t-s\|_{2}^{2}/\delta}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \sum_{\alpha\geq 0}\frac{(\sqrt{{1}/{\delta} }{\mathsf{P}}(t-s))^{\alpha}}{\alpha!}h_{\alpha}(\sqrt{{1}/{\delta}}{\mathsf{P }}(t-s)).$ $$

By Lemma 2.2, it suffices to divide ${\mathbb{R}}^{w}$ instead of ${\mathbb{R}}^{d}$ into boxes and conduct Hermite expansion and Taylor expansion on the low-dimensional subspace. More specifically, given the initial source points, we can compute ${\mathsf{P}}$ by SVD or QR decomposition in $N\cdot w^{\omega-1}$-time(^2^22$\omega\approx 2.372$ is the fast matrix multiplication time exponent.), which is of smaller order than the FGT’s preprocessing time(^3^33In practice, we can run numerical algorithms such as randomized SVD that are very fast for low-rank matrices.). Then, we can project each point $s_{i}\in{\mathbb{R}}^{d}$ to $x_{i}:={\mathsf{P}}s_{i}\in{\mathbb{R}}^{w}$ for $i\in[N]$. The remaining procedure in preprocessing is the same as before, but directly working on the low-dimensional sources $\{x_{1},\dots,x_{N}\}$.

In the query phase, consider a target point $t$ in the subspace. We are supposed to compute

$$ $\displaystyle G(t)\approx\sum_{{\cal B}}\sum_{j\in{\cal B}}q_{j}\cdot e^{-\|t- s_{j}\|_{2}^{2}/\delta}.$ $$

By Lemma 2.2, we know that

$$ $\displaystyle G(t)\approx$ $\displaystyle\leavevmode\nobreak\ \sum_{\beta\leq p}C_{\beta}({\cal C})\left( \frac{{\mathsf{P}}(t-t_{\cal C})}{\sqrt{\delta}}\right)^{\beta}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \sum_{\beta\leq p}C_{\beta}({\cal C})\left( \frac{y-y_{\cal C}}{\sqrt{\delta}}\right)^{\beta},$ $$

where ${\cal C}$ is the target box that contains $t$, $y={\mathsf{P}}t$ and $y_{{\cal C}}={\mathsf{P}}t_{\cal C}$ projected points. Moreover, for each $\beta\leq p$ and target box ${\cal C}$, we have

$$ $\displaystyle\leavevmode\nobreak\ C_{\beta}({\cal C})$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{{ \cal B}}\sum_{\alpha\leq p}A_{\alpha}({\cal B})H_{\alpha+\beta}\left(\frac{{ \mathsf{P}}(s_{{\cal B}}-t_{{\cal C}})}{\sqrt{\delta}}\right),$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{{ \cal B}}\sum_{\alpha\leq p}A_{\alpha}({\cal B})H_{\alpha+\beta}\left(\frac{x_{ {\cal B}}-y_{{\cal C}}}{\sqrt{\delta}}\right).$ $$

Similarly, for each $\alpha\leq p$ and source box ${\cal B}$,

$$ $\displaystyle A_{\alpha}({\cal B})=\frac{(-1)^{\|\alpha\|_{1}}}{\alpha!}\sum_{ j\in{\cal B}}q_{j}\cdot\left(\frac{x_{j}-x_{\cal B}}{\sqrt{\delta}}\right)^{ \alpha}.$ $$

Therefore, each query is equivalent to being conducted in a $w$-dimensional space using our data structure, which takes $\log^{O(w)}(\|q\|_{1}/\varepsilon)$-time.

The update can be done in a similar way in the low-dimensional space using the procedure described in Section 2.3. Hence, each update (insertion or deletion) takes $\log^{O(w)}(\|q\|_{1}/\varepsilon)$.

###### Lemma 2.2 (Hermite projection lemma in low-dimensional space, informal version of Lemma F.3 ) .

### 2.6 Handling Points From Low-Dimensional Dynamic Spaces

We note that when we add a new source point to the data structure, the intrinsic rank of the data might change by 1 when the point is not in the subspace. For an inserting source point $s$, consider the rank-increasing case, i.e., $(I-{\mathsf{P}})s\neq 0$. Then, this new source point contributes to one new basis

$$ $\displaystyle u:=\frac{(I-{\mathsf{P}})s}{\|(I-{\mathsf{P}})s\|_{2}}.$ $$

And we can update the projection matrix ${\mathsf{P}}$ by
$\begin{bmatrix}{\mathsf{P}}&u\end{bmatrix}\in{\mathbb{R}}^{(w+1)\times d}$.

However, as the subspace is changed, we need to maintain the intermediate variables $A_{\alpha}({\cal B}),C_{\beta}({\cal C})$. It is easy to observe that for the original projected source and target points or boxes, they can easily be “lifted” to the new subspace by setting zero to the $(w+1)$-th coordinate. We show how to update $A_{\alpha}({\cal B})$ efficiently. For each source box ${\cal B}$ and $\alpha\leq p$, we have

$$ $\displaystyle\leavevmode\nobreak\ A_{(\alpha,0)}^{\mathsf{new}}({\cal B})$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \frac{(-1)^{\|\alpha\|_{1}}\cdot(-1)^{i}}{ \alpha!\cdot i!}\sum_{j\in{\cal B}}q_{j}\cdot\left(\frac{x^{\prime}_{j}-x^{ \prime}_{\cal B}}{\sqrt{\delta}}\right)^{(\alpha,i)}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ A_{\alpha}({\cal B}),$ $$

where $x^{\prime}_{j}$ denotes the lifted point.
And $A_{(\alpha,1)}^{\mathsf{new}}({\cal B})=0$ for all $i>0$. Similarly, for each target box ${\cal C}$,

$$ $\displaystyle\leavevmode\nobreak\ C_{(\beta,i)}^{\mathsf{new}}({\cal C})$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \frac{(-1)^{\|\beta\|_{1}}(-1)^{i}}{\beta!i!}$ $\displaystyle\cdot$ $\displaystyle\leavevmode\nobreak\ \sum_{{\cal B}}\sum_{\alpha\leq p}\sum_{j=0} ^{p}A_{(\alpha,j)}^{\mathsf{new}}({\cal B})H_{(\alpha+\beta,i+j)}\left(\frac{x ^{\prime}_{{\cal B}}-y^{\prime}_{{\cal C}}}{\sqrt{\delta}}\right)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \frac{(-1)^{\|\beta\|_{1}}(-1)^{i}}{\beta!i!}$ $\displaystyle\cdot$ $\displaystyle\leavevmode\nobreak\ \sum_{{\cal B}}\sum_{\alpha\leq p}A_{\alpha} ({\cal B})H_{\alpha+\beta}\left(\frac{x_{{\cal B}}-y_{{\cal C}}}{\sqrt{\delta} }\right)\cdot h_{i}(0)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \frac{(-1)^{i}}{i!}\cdot C_{\beta}({\cal C}).$ $$

Therefore, by enumerating all boxes ${\cal B},{\cal C}$ and indices $\alpha,\beta\leq p$, we can compute $A_{(\alpha,0)}^{\mathsf{new}}({\cal B})$ and $C_{(\beta,i)}^{\mathsf{new}}({\cal C})$ in $\log^{O(w)}(\|q\|_{1}/\varepsilon)$-time.

Then, we just follow the static subspace insertion procedure to insert the new source point $s$. In this way, we obtain a data structure that can handle dynamic low-rank subspaces.

## 3 Conclusion and Future Directions

In this paper, we study the Fast Gaussian Transform (FGT) in a dynamic setting and propose a dynamic data structure to maintain the source vectors that support very fast kernel density estimation, Mat-Vec queries ($\mathsf{K}\cdot q$), as well as updating the source vectors. We further show that the efficiency of our algorithm can be improved when the data points lie in a low-dimensional subspace. Our results are especially valuable when FGT is used in real-world applications with rapidly-evolving datasets, e.g., online regression, federated learning, etc.

One open problem in this direction is, can we compute $\mathsf{K}q$ in

$$ $\displaystyle O(N)+\log^{O(d)}(N/\varepsilon)$ $$

time? Currently, it takes $N\log^{O(d)}(N/\varepsilon)$ time even in the static setting. The lower bounds in [] indicate that this improvement is impossible for some “bad” kernels $\mathsf{K}$ which are very non-smooth. It remains open when $\mathsf{K}$ is a Gaussian-like kernel. It might be helpful to apply more complicated geometric data structures to maintain the interactions between data points.

Another open problem is, can we fast compute Mat-Vec product or KDE for slowly-decaying kernels? The main difficulty is the current FMM techniques cannot achieve high accuracy when the kernel decays slowly. New techniques might be required to resolve this problem.

## Appendix

Roadmap.
In Section A, we provide several notations and definitions about the Fast Multipole Method. In Section B, we present the formal statement of our main result. In Section C, we present our data-structures and algorithms. In Section D, we provide a complete and full for our results. In Section E, we prove several lemmas to control the error. In Section F, we generalize our results to low dimension subspace setting.

## Appendix A Preliminaries

We first give a quick overview of the high-level ideas of FMM in Section A.1. In Section A.2, we provide a complete description and proof of correctness for the fast Gaussian transform, where the kernel function is the Gaussian kernel. Although a number of researchers have used FMM in the past, most of the previous papers about FMM either focus on low-dimensional or low-error cases. We therefore focus on the superconstant-error, high dimensional case, and carefully analyze the joint dependence on $\varepsilon$ and $d$. We believe that our presentation of the original proof in Section A.2 is thus of independent interest to the community.

### A.1 FMM Background

We begin with a description of high-level ideas of the Fast Multipole Method (FMM). Let $\mathsf{K}:{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}_{+}$ denote a kernel function. The inputs to the FMM are $N$ sources $s_{1},s_{2},\cdots,s_{N}\in{\mathbb{R}}^{d}$ and $M$ targets $t_{1},t_{2},\cdots,t_{M}$. For each $i\in[N]$, source $s_{i}$ has associated ‘strength’ $q_{i}$. Suppose all sources are in a ‘box’ ${\cal B}$ and all the targets are in a ‘box’ ${\cal C}$. The goal is to evaluate

$$ $\displaystyle u_{j}=\sum_{i=1}^{N}\mathsf{K}(s_{i},t_{j})q_{i},\leavevmode \nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \forall j\in[M]$ $$

Intuitively, if $\mathsf{K}$ has some nice property (e.g. smooth), we can hope to approximate $\mathsf{K}$ in the following sense:

$$ $\displaystyle\mathsf{K}(s,t)\approx\sum_{p=0}^{P-1}B_{p}(s)\cdot C_{p}(t), \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ s\in{\cal B},t \in{\cal C}$ $$

for some functions $B_{p},C_{p}:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}$,
where $P$ is a small positive integer, usually called the *interaction rank* in the literature [, ].

Now, we can construct $u_{i}$ in two steps:

$$ $\displaystyle v_{p}=\sum_{i\in{\cal B}}B_{p}(s_{i})q_{i},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \forall p=0,1,\cdots,P-1,$ $$

and

$$ $\displaystyle\widetilde{u}_{j}=\sum_{p=0}^{P-1}C_{p}(t_{j})v_{p},\leavevmode \nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \forall i\in[M].$ $$

Intuitively, as long as ${\cal B}$ and ${\cal C}$ are well-separated, then $\widetilde{u}_{j}$ is very good estimation to $u_{j}$ even for small $P$, i.e., $|\widetilde{u}_{j}-u_{j}|<\varepsilon$.

Recall that, at the beginning of this section, we assumed that all the sources are in the the same box ${\cal B}$ and ${\cal C}$. This is not true in general. To deal with this, we can discretize the continuous space into a batch of boxes ${\cal B}_{1},{\cal B}_{2},\cdots$ and ${\cal C}_{1},{\cal C}_{2},\cdots$. For a box ${\cal B}_{l_{1}}$ and a box ${\cal C}_{l_{2}}$, if they are very far apart, then the interaction between points within them is small, and we can ignore it. If the two boxes are close, then we deal with them efficiently by truncating the high order expansion terms in $\mathsf{K}$ (only keeping the first $\log^{O(d)}(1/\varepsilon)$ terms). For each box, we will see that the number of nearby relevant boxes is at most $\log^{O(d)}(1/\varepsilon)$.

### A.2 Fast Gaussian Transform

Given $N$ vectors $s_{1},\cdots s_{N}\in{\mathbb{R}}^{d}$, $M$ vectors $t_{1},\cdots,t_{M}\in{\mathbb{R}}^{d}$ and a strength vector $q\in{\mathbb{R}}^{n}$, Greengard and Strain [] provided a fast algorithm for evaluating discrete Gauss transform

$$ $\displaystyle G(t_{i})=\sum_{j=1}^{N}q_{j}e^{-\|t_{i}-s_{j}\|^{2}/\delta}$ $$

for all $i\in[M]$ in $O(M+N)$ time. In this section, we re-prove the algorithm described in [], and determine the exact dependence on $\varepsilon$ and $d$ in the running time.

Without loss of generality, we can assume that all the sources $s_{j}$ and targets are belonging to the unit box ${\cal B}_{0}=[0,1]^{d}$. The reason is, if not, we can shift the origin and rescaling $\delta$.

Let $t$ and $s$ lie in $d$-dimensional Euclidean space ${\mathbb{R}}^{d}$, and consider the Gaussian

$$ $\displaystyle e^{-\|t-s\|_{2}^{2}}=e^{-\sum_{i=1}^{d}(t_{i}-s_{i})^{2}}$ $$

We begin with some definitions. One important tool we use is the Hermite polynomial, which is a well-known class of orthogonal polynomials with respect to Gaussian measure and widely used in analyzing Gaussian kernels.

###### Definition A.1 (One-dimensional Hermite polynomial, [ 48 ] ) .

The Hermite polynomials $\widetilde{h}_{n}:{\mathbb{R}}\rightarrow{\mathbb{R}}$ is defined as follows

$$ $\displaystyle\widetilde{h}_{n}(t)=(-1)^{n}e^{t^{2}}\frac{\mathrm{d}^{n}}{ \mathrm{d}t}e^{-t^{2}}$ $$

The first few Hermite polynomials are:

$$ $\displaystyle\widetilde{h}_{1}(t)=2t,\leavevmode\nobreak\ \widetilde{h}_{2}(t) =4t^{2}-2,\leavevmode\nobreak\ \widetilde{h}_{3}(t)=8t^{3}-12t,\leavevmode \nobreak\ \cdots$ $$

###### Definition A.2 (One-dimensional Hermite function, [ 48 ] ) .

The Hermite functions $h_{n}:{\mathbb{R}}\rightarrow{\mathbb{R}}$ is defined as follows

$$ $\displaystyle h_{n}(t)=e^{-t^{2}}\widetilde{h}_{n}(t)=(-1)^{n}\frac{\mathrm{d} ^{n}}{\mathrm{d}t}e^{-t^{2}}$ $$

We use the following Fact to simplify $e^{-(t-s)^{2}/\delta}$.

###### Fact A.3 .

For $s_{0}\in{\mathbb{R}}$ and $\delta>0$, we have

$$ $\displaystyle e^{-(t-s)^{2}/\delta}=\sum_{n=0}^{\infty}\frac{1}{n!}\cdot\left( \frac{s-s_{0}}{\sqrt{\delta}}\right)^{n}\cdot h_{n}\left(\frac{t-s_{0}}{\sqrt{ \delta}}\right)$ $$

and

$$ $\displaystyle e^{-(t-s)^{2}/\delta}=e^{-(t-s_{0})^{2}/\delta}\sum_{n=0}^{ \infty}\frac{1}{n!}\cdot\left(\frac{s-s_{0}}{\sqrt{\delta}}\right)^{n}\cdot \widetilde{h}_{n}\left(\frac{t-s_{0}}{\sqrt{\delta}}\right).$ $$

###### Lemma A.4 (Cramer’s inequality for one-dimensional, [ 49 ] ) .

For any $K<1.09$,

$$ $\displaystyle|\widetilde{h}_{n}(t)|\leq K2^{n/2}\sqrt{n!}e^{t^{2}/2}.$ $$

Using Cramer’s inequality (Lemma A.4), we have the following standard bound.

###### Lemma A.5 .

For any constant $K<1.09$,
we have

$$ $\displaystyle|h_{n}(t)|\leq K\cdot 2^{n/2}\cdot\sqrt{n!}\cdot e^{-t^{2}/2}.$ $$

Next, we will extend the above definitions and observations to the high dimensional case. To simplify the discussion, we define multi-index notation. A multi-index $\alpha=(\alpha_{1},\alpha_{2},\cdots,\alpha_{d})$ is a $d$-tuple of nonnegative integers, playing the role of a multi-dimensional index. For any multi-index $\alpha\in{\mathbb{R}}^{d}$ and any $t\in{\mathbb{R}}^{t}$, we write

$$ $\displaystyle\alpha!=\leavevmode\nobreak\ \prod_{i=1}^{d}(\alpha_{i}!), \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode \nobreak\ t^{\alpha}=\leavevmode\nobreak\ \prod_{i=1}^{d}t_{i}^{\alpha_{i}}, \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode \nobreak\ D^{\alpha}=\leavevmode\nobreak\ \partial_{1}^{\alpha_{1}}\partial_{2 }^{\alpha_{2}}\cdots\partial_{d}^{\alpha_{d}}.$ $$

where $\partial_{i}$ is the differential operator with respect to the $i$-th coordinate in ${\mathbb{R}}^{d}$. For integer $p$, we say $\alpha\leq p$ if $\alpha_{i}\leq p,\forall i\in[d]$; and we say $\alpha\geq p$ if $\alpha_{i}\geq p,\exists i\in[d]$. We use these definitions to guarantee that $\{\alpha\leq p\}\cup\{\alpha\geq p\}={\mathbb{N}}^{d}$.

We can now define multi-dimensional Hermite polynomial:

###### Definition A.6 (Multi-dimensional Hermite polynomial, [ 48 ] ) .

We define function $\widetilde{H}_{\alpha}:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}$ as follows:

$$ $\displaystyle\widetilde{H}_{\alpha}(t)=\prod_{i=1}^{d}\widetilde{h}_{\alpha_{i }}(t_{i}).$ $$

###### Definition A.7 (Multi-dimensional Hermite function, [ 48 ] ) .

We define function $H_{\alpha}:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}$ as follows:

$$ $\displaystyle H_{\alpha}(t)=\prod_{i=1}^{d}h_{\alpha_{i}}(t_{i}).$ $$

It is easy to see that $H_{\alpha}(t)=e^{-\|t\|_{2}^{2}}\cdot\widetilde{H}_{\alpha}(t)$

The Hermite expansion of a Gaussian in ${\mathbb{R}}^{d}$ is

$$ $\displaystyle e^{-\|t-s\|_{2}^{2}}=\sum_{\alpha\geq 0}\frac{(t-s_{0})^{\alpha} }{\alpha!}h_{\alpha}(s-s_{0}).$ (4) $$

Cramer’s inequality generalizes to

###### Lemma A.8 (Cramer’s inequality for multi-dimensional case, [ 46 , 1 ] ) .

Let $K<(1.09)^{d}$, then

$$ $\displaystyle|\widetilde{H}_{\alpha}(t)|\leq K\cdot e^{\|t\|_{2}^{2}/2}\cdot 2 ^{\|\alpha\|_{1}/2}\cdot\sqrt{\alpha!}$ $$

and

$$ $\displaystyle|H_{\alpha}(t)|\leq K\cdot e^{-\|t\|_{2}^{2}/2}\cdot 2^{\|\alpha \|_{1}/2}\cdot\sqrt{\alpha!}.$ $$

The Taylor series of $H_{\alpha}$ is

$$ $\displaystyle H_{\alpha}(t)=\sum_{\beta\geq 0}\frac{(t-t_{0})^{\beta}}{\beta!} (-1)^{\|\beta\|_{1}}H_{\alpha+\beta}(t_{0}).$ (5) $$

###### Definition A.1 (One-dimensional Hermite polynomial, [ 48 ] ) .

###### Definition A.2 (One-dimensional Hermite function, [ 48 ] ) .

###### Fact A.3 .

###### Lemma A.4 (Cramer’s inequality for one-dimensional, [ 49 ] ) .

###### Lemma A.5 .

###### Definition A.6 (Multi-dimensional Hermite polynomial, [ 48 ] ) .

###### Definition A.7 (Multi-dimensional Hermite function, [ 48 ] ) .

###### Lemma A.8 (Cramer’s inequality for multi-dimensional case, [ 46 , 1 ] ) .

## Appendix B Our Result

Figure: Figure 3: An illustration of deleting a source point from the data structure.
Refer to caption: x3.png

### B.1 Properties of Kernel Function

[] identified the three key properties of kernel functions $\mathsf{K}(s,t)=f(\|s-t\|_{2})$ which allow sub-quadratic matrix-vector multiplication via the fast Multipole method. Our dynamic algorithm will work for any kernel satisfying these properties.

###### Definition B.1 (Properties of general kernel function, [ 1 ] ) .

We define the following properties of the function $f:{\mathbb{R}}\rightarrow{\mathbb{R}}_{+}$:

- •
P1: $f$ is non-increasing, i.e., $f(x)\leq f(y)$ when $x\geq y$.
- •
P2: $f$ is decreasing fast, i.e., $f(\Theta(\log(1/\varepsilon)))\leq\varepsilon$.
- •
P3: $f$’s Hermite expansion and Taylor expansion are truncateable: the truncation error of the first $(\log^{d}(1/\varepsilon))$ terms in the Hermite and Taylor expansion of $\mathsf{K}$ is at most $\varepsilon$.

###### Remark B.2 .

We note that P3 can be replaced with the following more general property:

- •
P4: $\mathsf{K}:{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}$ is *$\{\phi_{\alpha}\}_{\alpha\in{\mathbb{N}}^{d}}$-expansionable*: there exist constants $c_{\alpha}$ that only depend on $\alpha\in{\mathbb{N}}^{d}$ and functions $\phi_{\alpha}:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}$ such that
$\displaystyle\mathsf{K}(s,t)=\sum_{\alpha\in{\mathbb{N}}^{d}}c_{\alpha}\cdot(s
-s_{0})^{\alpha}\cdot\phi_{\alpha}(t-s_{0})$
for any $s_{0}\in{\mathbb{R}}^{d}$ and $s$ close to $s_{0}$. Moreover, the truncation error of the first $(\log^{d}(1/\varepsilon))$ terms is $\leq\varepsilon$.

###### Remark B.3 .

Two examples of kernels that satisfy Properties 1 and 2 are:

- •
$\mathsf{K}(s,t)=e^{-\alpha\|s-t\|^{2}}$ for any $\alpha\in{\mathbb{R}}_{+}$.
- •
$\mathsf{K}(s,t)=e^{-\alpha\|s-t\|^{2p}}$ for any $p\in{\mathbb{N}}_{+}$.

###### Definition B.1 (Properties of general kernel function, [ 1 ] ) .

###### Remark B.2 .

###### Remark B.3 .

### B.2 Dynamic FGT

In this section, we present our main result. We first define the dynamic density-estimation maintenance problem with respect to the Gaussian kernel.

Figure: Algorithm 1 Informal version of Algorithm 2, 3, 4 and 5.

###### Definition B.4 (Dynamic FGT Problem) .

We wish to design a data-structure that efficiently supports any sequence of the following operations:

- •
Init($S\subset{\mathbb{R}}^{d},q\in{\mathbb{R}}^{|S|},\varepsilon\in{\mathbb{R}}$) Let $N=|S|$. The data structure is given $N$ source points $s_{1},\cdots,s_{N}\in{\mathbb{R}}^{d}$ with their charge $q_{1},\cdots,q_{N}\in{\mathbb{R}}$.
- •
Insert$(s\in{\mathbb{R}}^{d},\leavevmode\nobreak\ q_{s}\in{\mathbb{R}})$ Add the source point $s$ with its charge $q_{s}$ to the point set $S$.
- •
Delete$(s\in{\mathbb{R}}^{d})$ Delete $s$ (and its charge $q_{s}$) from the point set $S$.
- •
KDE-Query$(t\in{\mathbb{R}}^{d})$ Output $\widetilde{G}$ such that $G(t)-\varepsilon\leq\widetilde{G}\leq G(t)+\varepsilon$.

The main result of this paper is a fully-dynamic data structure supporting all of the above operations in *polylogarithmic* time:

###### Theorem B.5 (Dynamic FGT Data Structure) .

Given $N$ vectors $S=\{s_{1},\cdots,s_{N}\}\subset{\mathbb{R}}^{d}$,
a number $\delta>0$, and a vector $q\in{\mathbb{R}}^{N}$, let $G:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}$ be defined as $G(t)=\sum_{i=1}^{N}q_{i}\cdot\mathsf{K}(s_{i},t)$ denote the *kernel-density* of $t$ with respect to $S$, where $\mathsf{K}(s_{i},t)=f(\|s_{i}-t\|_{2})$ for $f$ satisfying the properties in Definition B.1 .
There is a dynamic data structure
that supports the following operations:

- •
Init$()$ (Algorithm 2) Preprocess in $N\cdot\log^{O(d)}(\|q\|_{1}/\varepsilon)$ time.
- •
KDE-Query$(t\in{\mathbb{R}}^{d})$ (Algorithm 3) Output $\widetilde{G}$ such that $G(t)-\varepsilon\leq\widetilde{G}\leq G(t)+\varepsilon$ in $\log^{O(d)}(\|q\|_{1}/\varepsilon)$ time.
- •
Insert$(s\in{\mathbb{R}}^{d},q_{s}\in{\mathbb{R}})$ (Algorithm 4) For any source point $s\in{\mathbb{R}}^{d}$ and its charge $q_{s}$, update the data structure by adding this source point in $\log^{O(d)}(\|q\|_{1}/\varepsilon)$ time.
- •
Delete$(s\in{\mathbb{R}}^{d})$ (Algorithm 5) For any source point $s\in{\mathbb{R}}^{d}$ and its charge $q_{s}$, update the data structure by deleting this source point in $\log^{O(d)}(\|q\|_{1}/\varepsilon)$ time.
- •
Query$(q\in{\mathbb{R}}^{N})$ (Algorithm 3) Output $\widetilde{\mathsf{K}q}\in{\mathbb{R}}^{N}$ such that $\|\widetilde{\mathsf{K}q}-\mathsf{K}q\|_{\infty}\leq\varepsilon$, where $\mathsf{K}\in{\mathbb{R}}^{N\times N}$ is defined by $\mathsf{K}_{i,j}=\mathsf{K}(s_{i},s_{j})$ in $N\log^{O(d)}(\|q\|_{1}/\varepsilon)$ time.

###### Remark B.6 .

The Query time can be further reduced when the change of the charge vector $q$ is sparsely changed between two consecutive queries. More specifically, let $\Delta:=\|q^{\mathrm{new}}-q\|_{0}$ be the number of changed coordinates of $q$. Then, Query can be done in $\widetilde{O}_{d}(\Delta)$ time.

###### Definition B.4 (Dynamic FGT Problem) .

###### Theorem B.5 (Dynamic FGT Data Structure) .

###### Remark B.6 .

## Appendix C Algorithms

Figure: Algorithm 2 This algorithm are the init part of Theorem B.5.

Figure: Algorithm 3 This algorithm is the query part of Theorem B.5.

Figure: Algorithm 4 This algorithm is the update part of Theorem B.5.

Figure: Algorithm 5 This algorithm is another update part of Theorem B.5.

## Appendix D Analysis

###### Proof of Theorem B.5 .

Correctness of KDE-Query.
Algorithm 2 accumulates all sources into truncated Hermite expansions and transforms all Hermite expansions into Taylor expansions via Lemma E.5, thus it can approximate the function $G(t)$ by

$$ $\displaystyle G(t)=$ $\displaystyle\leavevmode\nobreak\ \sum_{{\cal B}}\sum_{j\in{\cal B}}q_{j}\cdot e ^{-\|t-s_{j}\|_{2}^{2}/\delta}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \sum_{\beta\leq p}C_{\beta}\left(\frac{t-t_{ \cal C}}{\sqrt{\delta}}\right)^{\beta}+\mathrm{Err}_{T}(p)+\mathrm{Err}_{H}(p)$ $$

where $|\mathrm{Err}_{H}(p)|+|\mathrm{Err}_{T}(p)|\leq Q\cdot\varepsilon$ by $p=\log(\|q\|_{1}/\varepsilon)$,

$$ $\displaystyle C_{\beta}=\frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{{\cal B}}\sum _{\alpha\leq p}A_{\alpha}({\cal B})H_{\alpha+\beta}\left(\frac{s_{{\cal B}}-t_ {{\cal C}}}{\sqrt{\delta}}\right)$ $$

and the coefficients $A_{\alpha}({\cal B})$ are defined as Eq. (3).

Running time of KDE-Query.
In line 17, it takes $O(p^{d}N)$ time to compute all the Hermite expansions, i.e., to compute the coefficients $A_{\alpha}({\cal B})$ for all $\alpha\leq p$ and all sources boxes ${\cal B}$.

Making use of the large product in the definition of $H_{\alpha+\beta}$, we see that the time to compute the $p^{d}$ coefficients of $C_{\beta}$ is only $O(dp^{d+1})$ for each box ${\cal B}$ in the range. Thus, we know for each target box ${\cal C}$, the running time is $O((2k+1)^{d}dp^{d+1})$, thus the total time in line 23 is

$$ $\displaystyle O(N(C)\cdot(2k+1)^{d}dp^{d+1}).$ $$

Finally we need to evaluate the appropriate Taylor series for each target $t_{i}$, which can be done in $O(p^{d}M)$ time in line 4. Putting it all together, Algorithm 2 takes time

$$ $\displaystyle\leavevmode\nobreak\ O((2k+1)^{d}dp^{d+1}N(C))+O(p^{d}N)+O(p^{d}M)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ O\left((M+N)\log^{O(d)}(\|q\|_{1}/ \varepsilon)\right).$ $$

Correctness of Update.
Algorithm 4 and Algorithm 5 maintains $C_{\beta}$ as follows,

$$ $\displaystyle C_{\beta}=\frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{{\cal B}}\sum _{\alpha\leq p}A_{\alpha}({\cal B})H_{\alpha+\beta}\left(\frac{s_{{\cal B}}-t_ {{\cal C}}}{\sqrt{\delta}}\right)$ $$

where $A_{\alpha}({\cal B})$ is given by

$$ $\displaystyle A_{\alpha}({\cal B})=\frac{(-1)^{\|\alpha\|_{1}}}{\alpha!}\sum_{ j\in{\cal B}}q_{j}\cdot\left(\frac{s_{j}-s_{\cal B}}{\sqrt{\delta}}\right)^{ \alpha}.$ $$

Therefore, the correctness follows similarly from Algorithm 2.

Running time of Update.
In line 12, it takes $O(p^{d})$ time to update all the Hermite expansions, i.e. to update the coefficients $A_{\alpha}({\cal B})$ for all $\alpha\leq p$ and all sources boxes ${\cal B}$.

Making use of the large product in the definition of $H_{\alpha+\beta}$, we see that the time to compute the $p^{d}$ coefficients of $C_{\beta}$ is only $O(dp^{d+1})$ for each box ${\cal C}_{l}\in\mathsf{nb}(\mathcal{B}_{k})$. Thus, thus the total time in line 17 is

$$ $\displaystyle O((2k+1)^{d}dp^{d+1}).$ $$

Correctness of Query.
To compute $\mathsf{K}q$ for a given $q\in{\mathbb{R}}^{d}$, notice that for any $i\in[N]$,

$$ $\displaystyle(\mathsf{K}q)_{i}=$ $\displaystyle\leavevmode\nobreak\ \sum_{j=1}^{N}q_{j}\cdot e^{-\|s_{i}-s_{j}\| _{2}^{2}/\delta}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ G(s_{i}).$ $$

Hence, this problem reduces to $N$ KDE-Query() calls. And the additive error guarantee for $G(t)$ immediately gives the $\ell_{\infty}$-error guarantee for $\mathsf{K}q$.

Running time of Query.
We first build the data structure with the charge vector $q$ given in the query, which takes $\widetilde{O}_{d}(N)$ time. Then, we perform $N$ KDE-Query, each takes $\widetilde{O}_{d}(1)$. Hence, the total running time is $\widetilde{O}_{d}(N)$.

We note that when the charge vector $q$ is slowly changing, i.e., $\Delta:=\|q^{\mathrm{new}}-q\|_{0}\leq o(N)$, we can Update the source vectors whose charges are changed. Since each Insert or Delete takes $\widetilde{O}_{d}(1)$ time, it will take $\widetilde{O}_{d}(\Delta)$ time to update the data structure.

Then, consider computing $\mathsf{K}q^{\mathrm{new}}$ in this setting. We note that each source box can only affect $\widetilde{O}_{d}(1)$ other target boxes, where the target vectors are just the source vectors in this setting. Hence, there are at most $\widetilde{O}_{d}(\Delta)$ boxes whose $C_{\beta}$ is changed. Let ${\cal S}$ denote the indices of source vectors in these boxes. Since

$$ $\displaystyle G(s_{i})=\sum_{\beta\leq p}C_{\beta}(\mathcal{B}_{k})\cdot\left( \frac{s_{i}-s_{\mathcal{B}_{k}}}{\sqrt{\delta}}\right)^{\beta},$ $$

we get that there are at most $\widetilde{O}_{d}(\Delta)$ coordinates of $\mathsf{K}q^{\mathrm{new}}$ that are significantly changed from $\mathsf{K}q$, and we only need to re-compute $G(s_{i})$ for $i\in{\cal S}$. If we assume that the source vectors are well-separated, i.e., $|{\cal S}|=O(\delta)$, the total computational cost is $\widetilde{O}_{d}(\Delta)$.

Therefore, when the change of the charge vector $q$ is sparse, $\mathsf{K}q$ can be computed in sublinear time.
∎

###### Proof of Theorem B.5 .

## Appendix E Error Estimation

This section provides several technical lemma that are used in Appendix D. We first give a definition.

###### Definition E.1 (Hermite expansion and coefficients) .

Let ${\cal B}$ denote a box with center $s_{\cal B}\in{\mathbb{R}}^{d}$ and side length $r\sqrt{2\delta}$ with $r<1$.
If source $s_{j}$ is in box ${\cal B}$, we will simply denote as $j\in{\cal B}$. Then the Gaussian evaluation from the sources in box ${\cal B}$ is,

$$ $\displaystyle G(t)=\sum_{j\in{\cal B}}q_{j}\cdot e^{-\|t-s_{j}\|_{2}^{2}/ \delta}.$ $$

The Hermite expansion of $G(t)$ is

$$ $\displaystyle G(t)=\sum_{\alpha\geq 0}A_{\alpha}\cdot H_{\alpha}\left(\frac{t- s_{\cal B}}{\sqrt{\delta}}\right),$ (6) $$

where the coefficients $A_{\alpha}$ are defined by

$$ $\displaystyle A_{\alpha}=\frac{1}{\alpha!}\sum_{j\in{\cal B}}q_{j}\cdot\left( \frac{s_{j}-s_{\cal B}}{\sqrt{\delta}}\right)^{\alpha}$ (7) $$

The rest of this section will present a batch of Lemmas that bound the error of the function truncated at certain degree of Taylor and Hermite expansion.

We first upper bound the truncation error of Hermite expansion.

###### Lemma E.2 (Truncated Hermite expansion) .

Let $p$ denote an integer, let $\mathrm{Err}_{H}(p)$ denote the error after truncating the series $G(t)$ (as defined in Eq. (6)) after $p^{d}$ terms, i.e.,

$$ $\displaystyle\mathrm{Err}_{H}(p)=\sum_{\alpha\geq p}A_{\alpha}\cdot H_{\alpha} \left(\frac{t-s_{\cal B}}{\sqrt{\delta}}\right).$ (8) $$

Then we have

$$ $\displaystyle|\mathrm{Err}_{H}(p)|\leq\frac{\sum_{j\in{\cal B}}|q_{j}|}{(1-r)^ {d}}\sum_{k=0}^{d-1}\binom{d}{k}(1-r^{p})^{k}\left(\frac{r^{p}}{\sqrt{p!}} \right)^{d-k}$ $$

where $r\leq\frac{1}{2}$.

###### Proof.

Using Eq. (4) to expand each Gaussian (see Definition E.1) in the

$$ $\displaystyle G(t)=\sum_{j\in{\cal B}}q_{j}\cdot e^{-\|t-s_{j}\|_{2}^{2}/\delta}$ $$

into a Hermite series about $s_{\cal B}$:

$$ $\displaystyle\sum_{j\in{\cal B}}q_{j}\sum_{\alpha\geq 0}\frac{1}{\alpha!}\cdot \left(\frac{s_{j}-s_{\cal B}}{\sqrt{\delta}}\right)^{\alpha}\cdot H_{\alpha} \left(\frac{t-s_{\cal B}}{\sqrt{\delta}}\right)$ $$

and swap the summation over $\alpha$ and $j$ to obtain the desired form:

$$ $\displaystyle\sum_{\alpha\geq 0}\left(\frac{1}{\alpha!}\sum_{j\in{\cal B}}q_{j }\cdot\left(\frac{s_{j}-s_{\cal B}}{\sqrt{\delta}}\right)^{\alpha}\right)H_{ \alpha}\left(\frac{t-s_{\cal B}}{\sqrt{\delta}}\right)=\sum_{\alpha\geq 0}A_{ \alpha}H_{\alpha}\left(\frac{t-s_{\cal B}}{\sqrt{\delta}}\right).$ $$

Here, the truncation error bound is due to Lemma A.8 and the standard equation for the tail of a geometric series.

To formally bound the truncation error, we first rewrite the Hermit expansion as follows

$$ $\displaystyle e^{-\frac{\|t-s_{j}\|_{2}^{2}}{\delta}}=$ $\displaystyle\leavevmode\nobreak\ \prod_{i=1}^{d}\left(\sum_{n_{i}=1}^{p-1} \frac{1}{n_{i}!}\left(\frac{(s_{j})_{i}-(s_{\cal B})_{i}}{\sqrt{\delta}}\right )^{n_{i}}h_{n_{i}}\left(\frac{t_{i}-(s_{\cal B})_{i}}{\sqrt{\delta}}\right)\right.$ $\displaystyle\leavevmode\nobreak\ +\left.\sum_{n_{i}=p}^{\infty}\frac{1}{n_{i} !}\left(\frac{(s_{j})_{i}-(s_{\cal B})_{i}}{\sqrt{\delta}}\right)^{n_{i}}h_{n_ {i}}\left(\frac{t_{i}-(s_{\cal B})_{i}}{\sqrt{\delta}}\right)\right)$ (9) $$

Notice from Cramer’s inequality (Lemma A.5),

$$ $\displaystyle h_{n_{i}}\left(\frac{t_{i}-(s_{\cal B})_{i}}{\sqrt{\delta}} \right)\leq\sqrt{n!}\cdot 2^{n/2}\cdot e^{-(t_{i}-(s_{\cal B})_{i})^{2}/2}.$ $$

Therefore we can use properties of the geometric series (notice $\frac{(s_{j})_{i}-(s_{\cal B})_{i}}{\sqrt{\delta}}\leq r/\sqrt{2}$) to bound each term in the product as follows

$$ $\displaystyle\sum_{n_{i}=1}^{p-1}\frac{1}{n_{i}!}\left(\frac{(s_{j})_{i}-(s_{ \cal B})_{i}}{\sqrt{\delta}}\right)^{n_{i}}h_{n_{i}}\left(\frac{t_{i}-(s_{\cal B })_{i}}{\sqrt{\delta}}\right)\leq\frac{1-r^{p}}{1-r},$ (10) $$

and

$$ $\displaystyle\sum_{n_{i}=p}^{\infty}\frac{1}{n_{i}!}\left(\frac{(s_{j})_{i}-(s _{\cal B})_{i}}{\sqrt{\delta}}\right)^{n_{i}}h_{n_{i}}\left(\frac{t_{i}-(s_{ \cal B})_{i}}{\sqrt{\delta}}\right)\leq\frac{1}{\sqrt{p!}}\cdot\frac{r^{p}}{1- r}.$ (11) $$

Now we come back to bound Eq. (8) as follows

$$ $\displaystyle\mathrm{Err}_{H}(p)=$ $\displaystyle\leavevmode\nobreak\ \sum_{j\in{\cal B}}q_{j}\sum_{\alpha\geq p} \frac{1}{\alpha!}\cdot\left(\frac{s_{j}-s_{\cal B}}{\sqrt{\delta}}\right)^{ \alpha}\cdot H_{\alpha}\left(\frac{t-s_{\cal B}}{\sqrt{\delta}}\right)$ $\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ \left(\sum_{j\in{\cal B}}|q_{j}|\right)\left (e^{-\frac{\|t-s_{j}\|_{2}^{2}}{\delta}}-\prod_{j=1}^{d}\left(\sum_{n_{i}=1}^{ p-1}\frac{1}{n_{i}!}\left(\frac{(s_{j})_{i}-(s_{\cal B})_{i}}{\sqrt{\delta}} \right)^{n_{i}}h_{n_{i}}\left(\frac{t_{i}-(s_{\cal B})_{i}}{\sqrt{\delta}} \right)\right)\right)$ $\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ \frac{\sum_{j\in{\cal B}}|q_{j}|}{(1-r)^{d}} \sum_{k=0}^{d-1}\binom{d}{k}(1-r^{p})^{k}\left(\frac{r^{p}}{\sqrt{p!}}\right)^ {d-k}$ $$

where the first step comes from definition, the second step comes from Eq. (E) and the last step comes from Eq. (10) and Eq. (11) and binomial expansion.

∎

###### Remark E.3 .

By Stirling’s formula, it is easy to see that when we take $p=\log(\|q\|_{1}/\varepsilon)$, this error will be bounded by $\|q\|_{1}\cdot\varepsilon$.

The Lemma E.4 shows how to convert a Hermite expansion at location $s_{\cal B}$ into a Taylor expansion at location $t_{\cal C}$. Intuitively, the Taylor series converges rapidly in the box (that has side length $r\sqrt{2\delta}$ center around $t_{\cal C}$, where $r\in(0,1)$).

###### Lemma E.4 (Hermite expansion with truncated Taylor expansion) .

Suppose the Hermite expansion of $G(t)$ is given by Eq. (6), i.e.,

$$ $\displaystyle G(t)=\sum_{\alpha\geq 0}A_{\alpha}\cdot H_{\alpha}\left(\frac{t- s_{\cal B}}{\sqrt{\delta}}\right).$ (12) $$

Then, the Taylor expansion of $G(t)$ at an arbitrary point $t_{0}$ can be written as:

$$ $\displaystyle G(t)=\sum_{\beta\geq 0}B_{\beta}\left(\frac{t-t_{0}}{\sqrt{ \delta}}\right)^{\beta}.$ (13) $$

where the coefficients $B_{\beta}$ are defined as

$$ $\displaystyle B_{\beta}=\frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{\alpha\geq 0} (-1)^{\|\alpha\|_{1}}A_{\alpha}\cdot H_{\alpha+\beta}\left(\frac{s_{\cal B}-t_ {0}}{\sqrt{\delta}}\right).$ (14) $$

Let $\mathrm{Err}_{T}(p)$ denote the error by truncating the Taylor expansion after $p^{d}$ terms, in the box ${\cal C}$ (that has center at $t_{\cal C}$ and side length $r\sqrt{2\delta}$), i.e.,

$$ $\displaystyle\mathrm{Err}_{T}(p)=\sum_{\beta\geq p}B_{\beta}\left(\frac{t-t_{ \cal C}}{\sqrt{\delta}}\right)^{\beta}$ $$

Then

$$ $\displaystyle|\mathrm{Err}_{T}(p)|\leq\frac{\sum_{j\in{\cal B}}|q_{j}|}{(1-r)^ {d}}\sum_{k=0}^{d-1}{d\choose k}(1-r^{p})^{k}\left(\frac{r^{p}}{\sqrt{p!}} \right)^{d-k}$ $$

where $r\leq 1/2$.

###### Proof.

Each Hermite function in Eq. (12) can be expanded into a Taylor series by means of Eq. (5). The expansion in Eq. (13) is due to swapping the order of summation.

Next, we will bound the truncation error. Using Eq. (7) for $A_{\alpha}$, we can rewrite $B_{\beta}$:

$$ $\displaystyle B_{\beta}=$ $\displaystyle\leavevmode\nobreak\ \frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{ \alpha\geq 0}(-1)^{\|\alpha\|_{1}}A_{\alpha}H_{\alpha+\beta}\left(\frac{s_{ \cal B}-t_{\cal C}}{\sqrt{\delta}}\right)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{ \alpha\geq 0}\left(\frac{(-1)^{\|\alpha\|_{1}}}{\alpha!}\sum_{j\in{\cal B}}q_{ j}\left(\frac{s_{j}-s_{\cal B}}{\sqrt{\delta}}\right)^{\alpha}\right)H_{\alpha +\beta}\left(\frac{s_{\cal B}-t_{\cal C}}{\sqrt{\delta}}\right)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{j \in{\cal B}}q_{j}\sum_{\alpha\geq 0}\frac{(-1)^{\|\alpha\|_{1}}}{\alpha!}\left (\frac{s_{j}-s_{\cal B}}{\sqrt{\delta}}\right)^{\alpha}\cdot H_{\alpha+\beta} \left(\frac{s_{\cal B}-t_{\cal C}}{\sqrt{\delta}}\right)$ $$

By Eq. (5), the inner sum is the Taylor expansion of $H_{\beta}((s_{j}-t_{\cal C})/\sqrt{\delta})$. Thus

$$ $\displaystyle B_{\beta}=\frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{j\in{\cal B}} q_{j}\cdot H_{\beta}\left(\frac{s_{j}-t_{\cal C}}{\sqrt{\delta}}\right)$ $$

and Cramer’s inequality implies

$$ $\displaystyle|B_{\beta}|\leq\frac{1}{\beta!}K\cdot Q_{B}2^{\|\beta\|_{1}/2} \sqrt{\beta!}=KQ_{B}\frac{2^{\|\beta\|_{1}/2}}{\sqrt{\beta!}}.$ $$

To formally bound the truncation error, we have

$$ $\displaystyle\mathrm{Err}_{T}(p)=$ $\displaystyle\leavevmode\nobreak\ \sum_{\beta\geq p}B_{\beta}\left(\frac{t-t_{ \cal C}}{\sqrt{\delta}}\right)^{\beta}$ $\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ KQ_{B}\left(\prod_{i=1}^{d}\left(\sum_{n_{i} =0}^{\infty}\frac{1}{\sqrt{n_{i}!}}2^{n_{i}/2}\left(\frac{t-t_{\mathcal{C}}}{ \delta}\right)^{n_{i}}\right)-\prod_{i=1}^{d}\left(\sum_{n_{i}=0}^{p-1}\frac{1 }{\sqrt{n_{i}!}}2^{n_{i}/2}\left(\frac{t-t_{\mathcal{C}}}{\delta}\right)^{n_{i }}\right)\right)$ $\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ \frac{\sum_{j\in{\cal B}}|q_{j}|}{(1-r)^{d}} \sum_{k=0}^{d-1}{d\choose k}(1-r^{p})^{k}\left(\frac{r^{p}}{\sqrt{p!}}\right)^ {d-k}$ $$

where the second step uses $|B_{\beta}|\leq KQ_{B}\frac{2^{\|\beta\|_{1}/2}}{\sqrt{\beta!}}$ and the rest are similar to those in Lemma E.2.
∎

For designing our algorithm, we would like to make a variant of Lemma E.4 that combines the truncations of Hermite expansion and Taylor expansion. More specifically, we first truncate the Taylor expansion of $G_{p}(t)$, and then truncate the Hermite expansion in Eq. (14) for the coefficients.

###### Lemma E.5 (Truncated Hermite expansion with truncated Taylor expansion) .

Let $G(t)$ be defined as Def E.1. For an integer $p$, let $G_{p}(t)$ denote the Hermite expansion of $G(t)$ truncated at $p$, i.e.,

$$ $\displaystyle G_{p}(t)=\sum_{\alpha\leq p}A_{\alpha}H_{\alpha}\left(\frac{t-s_ {\cal B}}{\sqrt{\delta}}\right).$ $$

The Taylor expansion of function $G_{p}(t)$ at an arbitrary point $t_{0}$ can be written as:

$$ $\displaystyle G_{p}(t)=\sum_{\beta\geq 0}C_{\beta}\cdot\left(\frac{t-t_{0}}{ \sqrt{\delta}}\right)^{\beta},$ $$

where the coefficients $C_{\beta}$ are defined as

$$ $\displaystyle C_{\beta}=\frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{\alpha\leq p} (-1)^{\|\alpha\|_{1}}A_{\alpha}\cdot H_{\alpha+\beta}\left(\frac{s_{\cal B}-t_ {\cal C}}{\sqrt{\delta}}\right).$ (15) $$

Let $\mathrm{Err}_{T}(p)$ denote the error in truncating the Taylor series after $p^{d}$ terms, in the box ${\cal C}$ (that has center $t_{\cal C}$ and side length $r\sqrt{2\delta}$), i.e.,

$$ $\displaystyle\mathrm{Err}_{T}(p)=\sum_{\beta\geq p}C_{\beta}\left(\frac{t-t_{ \cal C}}{\sqrt{\delta}}\right)^{\beta}.$ $$

Then, we have

$$ $\displaystyle|\mathrm{Err}_{T}(p)|\leq\frac{2\sum_{j\in{\cal B}}|q_{j}|}{(1-r) ^{d}}\sum_{k=0}^{d-1}{d\choose k}(1-r^{p})^{k}\left(\frac{r^{p}}{\sqrt{p!}} \right)^{d-k}$ $$

where $r\leq 1/2$.

###### Proof.

We can write $C_{\beta}$ in the following way:

$$ $\displaystyle C_{\beta}=$ $\displaystyle\leavevmode\nobreak\ \frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{j \in{\cal B}}q_{j}\sum_{\alpha\leq p}\frac{(-1)^{\|\alpha\|_{1}}}{\alpha!}\left (\frac{s_{j}-s_{\cal B}}{\sqrt{\delta}}\right)^{\alpha}\cdot H_{\alpha+\beta} \left(\frac{s_{\cal B}-t_{\cal C}}{\sqrt{\delta}}\right)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{j \in{\cal B}}q_{j}\left(\sum_{\alpha\geq 0}-\sum_{\alpha>p}\right)\frac{(-1)^{ \|\alpha\|_{1}}}{\alpha!}\left(\frac{s_{j}-s_{\cal B}}{\sqrt{\delta}}\right)^{ \alpha}\cdot H_{\alpha+\beta}\left(\frac{s_{\cal B}-t_{\cal C}}{\sqrt{\delta}}\right)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ B_{\beta}-\frac{(-1)^{\|\beta\|_{1}}}{\beta! }\sum_{j\in{\cal B}}q_{j}\sum_{\alpha>p}\frac{(-1)^{\|\alpha\|_{1}}}{\alpha!} \left(\frac{s_{j}-s_{\cal B}}{\sqrt{\delta}}\right)^{\alpha}\cdot H_{\alpha+ \beta}\left(\frac{s_{\cal B}-t_{\cal C}}{\sqrt{\delta}}\right)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ B_{\beta}+(C_{\beta}-B_{\beta})$ $$

Next, we have

$$ $\displaystyle|\mathrm{Err}_{T}(p)|\leq\left|\sum_{\beta\geq p}B_{\beta}\left( \frac{t-t_{\cal C}}{\sqrt{\delta}}\right)^{\beta}\right|+\left|\sum_{\beta\geq p }(C_{\beta}-B_{\beta})\cdot\left(\frac{t-t_{\cal C}}{\sqrt{\delta}}\right)^{ \beta}\right|$ (16) $$

Using Lemma E.4, we can upper bound the first term in the Eq. (16) by,

$$ $\displaystyle\left|\sum_{\beta\geq p}B_{\beta}\left(\frac{t-t_{\cal C}}{\sqrt{ \delta}}\right)^{\beta}\right|\leq$ $\displaystyle\leavevmode\nobreak\ \frac{\sum_{j\in{\cal B}}|q_{j}|}{(1-r)^{d}} \sum_{k=0}^{d-1}{d\choose k}(1-r^{p})^{k}\left(\frac{r^{p}}{\sqrt{p!}}\right)^ {d-k}.$ $$

Since we can similarly bound $C_{\beta}-B_{\beta}$ as follows

$$ $\displaystyle|C_{\beta}-B_{\beta}|\leq\frac{1}{\beta!}K\cdot Q_{B}2^{\|\beta\| _{1}/2}\sqrt{\beta!}\leq KQ_{B}\frac{2^{\|\beta\|_{1}/2}}{\sqrt{\beta!}},$ $$

we have the same bound for the second term

$$ $\displaystyle\leavevmode\nobreak\ \left|\sum_{\beta\geq p}(C_{\beta}-B_{\beta} )\left(\frac{t-t_{\cal C}}{\sqrt{\delta}}\right)^{\beta}\right|\leq\leavevmode \nobreak\ \frac{\sum_{j\in{\cal B}}|q_{j}|}{(1-r)^{d}}\sum_{k=0}^{d-1}{d \choose k}(1-r^{p})^{k}\left(\frac{r^{p}}{\sqrt{p!}}\right)^{d-k}.$ $$

∎

The proof of the following Lemma is almost identical, but it directly bounds the truncation error of Taylor expansion of the Gaussian kernel. We omit the proof here.

###### Lemma E.6 (Truncated Taylor expansion) .

Let $G_{s_{j}}(t):{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}$ be defined as

$$ $\displaystyle G_{s_{j}}(t)=q_{j}\cdot e^{-\|t-s_{j}\|_{2}^{2}/\delta}.$ $$

The Taylor expansion of $G_{s_{j}}(t)$ at $t_{\cal C}\in{\mathbb{R}}^{d}$ is:

$$ $\displaystyle G_{s_{j}}(t)=\sum_{\beta\geq 0}{\cal B}_{\beta}\left(\frac{t-t_{ \cal C}}{\sqrt{\delta}}\right)^{\beta},$ $$

where the coefficients $B_{\beta}$ is defined as

$$ $\displaystyle B_{\beta}=q_{j}\cdot\frac{(-1)^{\|\beta\|_{1}}}{\beta!}\cdot H_{ \beta}\left(\frac{s_{j}-t_{\cal C}}{\sqrt{\delta}}\right)$ $$

and the absolute value of the error (truncation after $p^{d}$ terms) can be upper bounded as

$$ $\displaystyle|\mathrm{Err}_{T}(p)|\leq\frac{\sum_{j\in{\cal B}}|q_{j}|}{(1-r)^ {d}}\sum_{k=0}^{d-1}{d\choose k}(1-r^{p})^{k}\left(\frac{r^{p}}{\sqrt{p!}} \right)^{d-k}$ $$

where $r\leq 1/2$.

###### Definition E.1 (Hermite expansion and coefficients) .

###### Lemma E.2 (Truncated Hermite expansion) .

###### Proof.

###### Remark E.3 .

###### Lemma E.4 (Hermite expansion with truncated Taylor expansion) .

###### Proof.

###### Lemma E.5 (Truncated Hermite expansion with truncated Taylor expansion) .

###### Proof.

###### Lemma E.6 (Truncated Taylor expansion) .

## Appendix F Low Dimension Subspace FGT

In this section, we consider FGT for data in a lower dimensional subspace of ${\mathbb{R}}^{d}$. The problem is formally defined below:

###### Problem F.1 (Dynamic FGT on a low dimensional set) .

Let $W$ be a subspace of ${\mathbb{R}}^{d}$ with dimension $\dim(S)=w\ll d$. Given $N$ source points $s_{1},\dots,s_{N}\in W$ with charges $q_{1},\dots,q_{N}$, and $M$ target points $t_{1},\dots,t_{M}\in W$, find a dynamic data structure that supports the following operations:

- •
Insert/Delete$(s_{i},q_{i})$ Insert or Delete a source point $s_{i}\in{\mathbb{R}}^{d}$ along with its “charge” $q_{i}\in{\mathbb{R}}$, in $\log^{O(w)}(\|q\|_{1}/\varepsilon)$ time.
- •
Density-Estimation$(t\in{\mathbb{R}}^{d})$ For any point $t\in{\mathbb{R}}^{d}$, output the kernel density of $t$ with respect to the sources, i.e., output $\widetilde{G}$ such that $G(t)-\varepsilon\leq\widetilde{G}\leq G(t)+\varepsilon$ in $\log^{O(w)}(\|q\|_{1}/\varepsilon)$ time.
- •
Query$(q\in{\mathbb{R}}^{N})$ Given an arbitrary query vector $q\in{\mathbb{R}}^{N}$, output $\widetilde{\mathsf{K}q}$ in $N\cdot\log^{O(w)}(\|q\|/\varepsilon)$ time.

Figure: Algorithm 6 Initialization of low-dim FGT.

Figure: Algorithm 7 This algorithm is the query part of Theorem F.2.

Figure: Algorithm 8 This algorithm is the update part of Theorem F.2.

Figure: Algorithm 9 This algorithm is another update part of Theorem F.2.

We generalize our dynamic data structure to solve Problem F.1, which is stated in the following theorem. The computational cost of each update or query depends on the intrinsic dimension $w$ instead of $d$.

###### Theorem F.2 (Low Rank Dynamic FGT Data Structure, formal version of Theorem 1.1 ) .

Let $W$ be a subspace of ${\mathbb{R}}^{d}$ with dimension $\dim(S)=w\ll d$. Given $N$ source points $s_{1},\dots,s_{N}\in W$ with charges $q_{1},\dots,q_{N}$, and $M$ target points $t_{1},\dots,t_{M}\in W$, a number $\delta>0$, and a vector $q\in{\mathbb{R}}^{N}$,
let $G:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}$ be defined as $G(t)=\sum_{i=1}^{N}q_{i}\cdot\mathsf{K}(s_{i},t)$ denote the *kernel-density* of $t$ with respect to $S$, where $\mathsf{K}(s_{i},t)=f(\|s_{i}-t\|_{2})$ for $f$ satisfying the properties in Definition B.1 .
There is a dynamic data structure
that supports the following operations:

- •
Init$()$ (Algorithm 6) Preprocess in $N\cdot\log^{O(w)}(\|q\|_{1}/\varepsilon)$ time.
- •
KDE-Query$(t\in{\mathbb{R}}^{d})$ (Algorithm 7) Output $\widetilde{G}$ such that $G(t)-\varepsilon\leq\widetilde{G}\leq G(t)+\varepsilon$ in $\log^{O(w)}(\|q\|_{1}/\varepsilon)$ time.
- •
Insert$(s\in{\mathbb{R}}^{d},q_{s}\in{\mathbb{R}})$ (Algorithm 8) For any source point $s\in{\mathbb{R}}^{d}$ and its charge $q_{s}$, update the data structure by adding this source point in $\log^{O(w)}(\|q\|_{1}/\varepsilon)$ time.
- •
Delete$(s\in{\mathbb{R}}^{d})$ (Algorithm 9) For any source point $s\in{\mathbb{R}}^{d}$ and its charge $q_{s}$, update the data structure by deleting this source point in $\log^{O(w)}(\|q\|_{1}/\varepsilon)$ time.
- •
Query$(q\in{\mathbb{R}}^{N})$ (Algorithm 7) Output $\widetilde{\mathsf{K}q}\in{\mathbb{R}}^{N}$ such that $\|\widetilde{\mathsf{K}q}-\mathsf{K}q\|_{\infty}\leq\varepsilon$, where $\mathsf{K}\in{\mathbb{R}}^{N\times N}$ is defined by $\mathsf{K}_{i,j}=\mathsf{K}(s_{i},s_{j})$ in $N\log^{O(w)}(\|q\|_{1}/\varepsilon)$ time.

###### Problem F.1 (Dynamic FGT on a low dimensional set) .

###### Theorem F.2 (Low Rank Dynamic FGT Data Structure, formal version of Theorem 1.1 ) .

### F.1 Projection Lemma

###### Lemma F.3 (Hermite projection lemma in low-dimensional space, formal version of Lemma 2.2 ) .

Given a subspace $B\in{\mathbb{R}}^{d\times w}$. Let $B^{\top}B=U\Lambda U^{\top}\in{\mathbb{R}}^{w\times w}$ denote the spectral decomposition where $U\in{\mathbb{R}}^{w\times w}$ and a diagonal matrix $\Lambda\in{\mathbb{R}}^{w\times w}$.

We define ${\mathsf{P}}=\Lambda^{-1/2}U^{-1}B^{\top}\in{\mathbb{R}}^{w\times d}$. Then we have for any $t,s\in{\mathbb{R}}^{d}$ from subspace $B$, the following equation holds

$$ $\displaystyle e^{-\|t-s\|_{2}^{2}/\delta}=\sum_{\alpha\geq 0}\frac{(\sqrt{{1}/ {\delta}}{\mathsf{P}}(t-s))^{\alpha}}{\alpha!}h_{\alpha}(\sqrt{{1}/{\delta}}{ \mathsf{P}}(t-s)).$ $$

###### Proof.

First, we know that

$$ $\displaystyle{\mathsf{P}}t=$ $\displaystyle\leavevmode\nobreak\ \Lambda^{-1/2}U^{-1}B^{\top}t$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \Lambda^{-1/2}U^{-1}B^{\top}Bx$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \Lambda^{-1/2}U^{-1}U\Lambda U^{\top}x$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \Lambda^{-1/2}\Lambda U^{\top}x$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \Lambda^{1/2}U^{\top}x$ (17) $$

where the first step follows from ${\mathsf{P}}=\Lambda^{-1/2}U^{-1}B^{\top}$, the second step follows from $t=Bx$ (since $t$ is from low dimension, then there is always a vector $x$), the third step follows $B^{\top}B=U\Lambda U^{\top}$, the forth step follows $U^{-1}U=I$, and the last step follows from $\Lambda^{-1/2}\Lambda=\Lambda^{1/2}$.

Compute the spectral decomposition $B^{\top}B=U\Lambda U^{\top}$, $U\in{\mathbb{R}}^{w\times w}$ is the orthonormal basis, $\Lambda=\mathrm{diag}(\lambda_{1},\dots,\lambda_{k})\in{\mathbb{R}}^{w\times w}$. Let $u_{i}\in{\mathbb{R}}^{w}$ denote the vector that is the transpose of $i$-th row $U\in{\mathbb{R}}^{w\times w}$.
Then we have

$$ $\displaystyle e^{-\|t-s\|_{2}^{2}/\delta}=$ $\displaystyle\leavevmode\nobreak\ e^{-(x-y)^{\top}B^{\top}B(x-y)/\delta}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ e^{-(x-y)^{\top}U\Lambda U^{\top}(x-y)/\delta}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \prod_{i=1}^{w}\Big{(}\sum_{n=1}^{\infty} \frac{1}{n!}(\sqrt{{\lambda_{i}}/{\delta}}\cdot u_{i}^{\top}(x-y))^{n}\cdot h_ {n}(\sqrt{{\lambda_{i}}/{\delta}}\cdot u_{i}^{\top}(x-y))\Big{)}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \sum_{\alpha\geq 0}\frac{\left(\sqrt{1/ \delta}\Lambda^{1/2}U^{\top}(x-y)\right)^{\alpha}}{\alpha!}\cdot h_{\alpha} \left(\sqrt{1/\delta}\Lambda^{1/2}U^{\top}(x-y)\right)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \sum_{\alpha\geq 0}\frac{\left(\sqrt{1/ \delta}\cdot{\mathsf{P}}(t-s)\right)^{\alpha}}{\alpha!}\cdot h_{\alpha}\left( \sqrt{1/\delta}\cdot{\mathsf{P}}(t-s)\right)$ $$

where the first step follows from changing the basis preserves the $\ell_{2}$-distance, the second step follows from $B^{\top}B=U\Lambda U^{\top}$, and the fifth step follows from Eq. (F.1).

∎

###### Lemma F.3 (Hermite projection lemma in low-dimensional space, formal version of Lemma 2.2 ) .

###### Proof.

### F.2 Proof of Main Result in Low-Dimensional Case

###### Proof of Theorem F.2 .

Correctness of KDE-Query.
Algorithm 6 accumulates all sources into truncated Hermite expansions and transforms all Hermite expansions into Taylor expansions via Lemma F.4. Thus it can approximate the function $G(t)$ by

$$ $\displaystyle G(t)=$ $\displaystyle\leavevmode\nobreak\ \sum_{{\cal B}}\sum_{j\in{\cal B}}q_{j}\cdot e ^{-\|t-s_{j}\|_{2}^{2}/\delta}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \sum_{\beta\leq p}C_{\beta}\left(\frac{{ \mathsf{P}}(t-t_{\cal C})}{\sqrt{\delta}}\right)^{\beta}+\mathrm{Err}_{T}(p)+ \mathrm{Err}_{H}(p)$ $$

where $|\mathrm{Err}_{H}(p)|+|\mathrm{Err}_{T}(p)|\leq Q\cdot\varepsilon$ by $p=\log(\|q\|_{1}/\varepsilon)$,

$$ $\displaystyle C_{\beta}=\frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{{\cal B}}\sum _{\alpha\leq p}A_{\alpha}({\cal B})H_{\alpha+\beta}\left(\frac{{\mathsf{P}}(s_ {{\cal B}}-t_{{\cal C}})}{\sqrt{\delta}}\right)$ $$

and the coefficients $A_{\alpha}({\cal B})$ are defined as Line 24.

Running time of KDE-Query.
In line 24, it takes $O(p^{w}N)$ time to compute all the Hermite expansions, i.e., to compute the coefficients $A_{\alpha}({\cal B})$ for all $\alpha\leq p$ and all source boxes ${\cal B}$.

Making use of the large product in the definition of $H_{\alpha+\beta}$, we see that the time to compute the $p^{w}$ coefficients of $C_{\beta}$ is only $O(dp^{w+1})$ for each box ${\cal B}$ in the range. Thus, we know for each target box ${\cal C}$, the running time is $O((2k+1)^{w}dp^{w+1})$, thus the total time in line 30 is

$$ $\displaystyle O(N(C)\cdot(2k+1)^{w}dp^{w+1}).$ $$

Finally, we need to evaluate the appropriate Taylor series for each target $t_{i}$, which can be done in $O(p^{w}M)$ time in line 4. Putting it all together, Algorithm 6 takes time

$$ $\displaystyle\leavevmode\nobreak\ O((2k+1)^{w}dp^{w+1}N(C))+O(p^{w}N)+O(p^{w}M)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ O\left((M+N)\log^{O(w)}(\|q\|_{1}/ \varepsilon)\right).$ $$

Correctness of Update.
Algorithm 8 and Algorithm 9 maintains $C_{\beta}$ as follows,

$$ $\displaystyle C_{\beta}=\frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{{\cal B}}\sum _{\alpha\leq p}A_{\alpha}({\cal B})H_{\alpha+\beta}\left(\frac{{\mathsf{P}}(s_ {{\cal B}}-t_{{\cal C}})}{\sqrt{\delta}}\right)$ $$

where $A_{\alpha}({\cal B})$ is given by

$$ $\displaystyle A_{\alpha}({\cal B})=\frac{(-1)^{\|\alpha\|_{1}}}{\alpha!}\sum_{ j\in{\cal B}}q_{j}\cdot\left(\frac{{\mathsf{P}}(s_{j}-s_{\cal B})}{\sqrt{ \delta}}\right)^{\alpha}.$ $$

Therefore, the correctness follows similarly from Algorithm 6.

Running time of Update.
In line 12, it takes $O(p^{w})$ time to update all the Hermite expansions, i.e. to update the coefficients $A_{\alpha}({\cal B})$ for all $\alpha\leq p$ and all sources boxes ${\cal B}$.

Making use of the large product in the definition of $H_{\alpha+\beta}$, we see that the time to compute the $p^{w}$ coefficients of $C_{\beta}$ is only $O(dp^{w+1})$ for each box ${\cal C}_{l}\in\mathsf{nb}(\mathcal{B})$. Thus, thus the total time in line 17 is

$$ $\displaystyle O((2k+1)^{w}dp^{w+1}).$ $$

Correctness of Query.
To compute $\mathsf{K}q$ for a given $q\in{\mathbb{R}}^{w}$, notice that for any $i\in[N]$,

$$ $\displaystyle(\mathsf{K}q)_{i}=$ $\displaystyle\leavevmode\nobreak\ \sum_{j=1}^{N}q_{j}\cdot e^{-\|s_{i}-s_{j}\| _{2}^{2}/\delta}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ G(s_{i}).$ $$

Hence, this problem reduces to $N$ KDE-Query() calls. And the additive error guarantee for $G(t)$ immediately gives the $\ell_{\infty}$-error guarantee for $\mathsf{K}q$.

Running time of Query.
We first build the data structure with the charge vector $q$ given in the query, which takes $\widetilde{O}_{d}(N)$ time. Then, we perform $N$ KDE-Query, each takes $\widetilde{O}_{d}(1)$. Hence, the total running time is $\widetilde{O}_{d}(N)$.

We note that when the charge vector $q$ is slowly changing, i.e., $\Delta:=\|q^{\mathrm{new}}-q\|_{0}\leq o(N)$, we can Update the source vectors whose charges are changed. Since each Insert or Delete takes $\widetilde{O}_{d}(1)$ time, it will take $\widetilde{O}_{d}(\Delta)$ time to update the data structure.

Then, consider computing $\mathsf{K}q^{\mathrm{new}}$ in this setting. We note that each source box can only affect $\widetilde{O}_{d}(1)$ other target boxes, where the target vectors are just the source vectors in this setting. Hence, there are at most $\widetilde{O}_{d}(\Delta)$ boxes whose $C_{\beta}$ is changed. Let ${\cal S}$ denote the indices of source vectors in these boxes. Since

$$ $\displaystyle G(s_{i})=\sum_{\beta\leq p}C_{\beta}(\mathcal{B}_{k})\cdot\left( \frac{{\mathsf{P}}(s_{i}-s_{\mathcal{B}_{k}})}{\sqrt{\delta}}\right)^{\beta},$ $$

we get that there are at most $\widetilde{O}_{d}(\Delta)$ coordinates of $\mathsf{K}q^{\mathrm{new}}$ that are significantly changed from $\mathsf{K}q$, and we only need to re-compute $G(s_{i})$ for $i\in{\cal S}$. If we assume that the source vectors are well-separated, i.e., $|{\cal S}|=O(\delta)$, the total computational cost is $\widetilde{O}_{d}(\Delta)$.

Therefore, when the change of the charge vector $q$ is sparse, $\mathsf{K}q$ can be computed in sublinear time.
∎

###### Lemma F.4 (Truncated Hermite expansion with truncated Taylor expansion (low dimension version of Lemma E.5 )) .

Let $G(t)$ be defined as Def E.1. For an integer $p$, let $G_{p}(t)$ denote the Hermite expansion of $G(t)$ truncated at $p$, i.e.,

$$ $\displaystyle G_{p}(t)=\sum_{\alpha\leq p}A_{\alpha}H_{\alpha}\left(\frac{{ \mathsf{P}}(t-s_{\cal B})}{\sqrt{\delta}}\right).$ $$

The Taylor expansion of function $G_{p}(t)$ at an arbitrary point $t_{0}$ can be written as:

$$ $\displaystyle G_{p}(t)=\sum_{\beta\geq 0}C_{\beta}\cdot\left(\frac{{\mathsf{P} }(t-t_{0})}{\sqrt{\delta}}\right)^{\beta},$ $$

where the coefficients $C_{\beta}$ are defined as

$$ $\displaystyle C_{\beta}=\frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{\alpha\leq p} (-1)^{\|\alpha\|_{1}}A_{\alpha}\cdot H_{\alpha+\beta}\left({\mathsf{P}}\frac{( s_{\cal B}-t_{\cal C})}{\sqrt{\delta}}\right).$ (18) $$

Let $\mathrm{Err}_{T}(p)$ denote the error in truncating the Taylor series after $p^{w}$ terms, in the box ${\cal C}$ (that has center $t_{\cal C}$ and side length $r\sqrt{2\delta}$), i.e.,

$$ $\displaystyle\mathrm{Err}_{T}(p)=\sum_{\beta\geq p}C_{\beta}\left(\frac{{ \mathsf{P}}(t-t_{\cal C})}{\sqrt{\delta}}\right)^{\beta}.$ $$

Then, we have

$$ $\displaystyle|\mathrm{Err}_{T}(p)|\leq\frac{2\sum_{j\in{\cal B}}|q_{j}|}{(1-r) ^{w}}\sum_{l=0}^{w-1}{w\choose l}(1-r^{p})^{l}\left(\frac{r^{p}}{\sqrt{p!}} \right)^{w-l}$ $$

where $r\leq 1/2$.

###### Proof.

We can write $C_{\beta}$ in the following way:

$$ $\displaystyle C_{\beta}=$ $\displaystyle\leavevmode\nobreak\ \frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{j \in{\cal B}}q_{j}\sum_{\alpha\leq p}\frac{(-1)^{\|\alpha\|_{1}}}{\alpha!}\left (\frac{{\mathsf{P}}(s_{j}-s_{\cal B})}{\sqrt{\delta}}\right)^{\alpha}\cdot H_{ \alpha+\beta}\left(\frac{{\mathsf{P}}(s_{\cal B}-t_{\cal C})}{\sqrt{\delta}}\right)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ \frac{(-1)^{\|\beta\|_{1}}}{\beta!}\sum_{j \in{\cal B}}q_{j}\left(\sum_{\alpha\geq 0}-\sum_{\alpha>p}\right)\frac{(-1)^{ \|\alpha\|_{1}}}{\alpha!}\left(\frac{{\mathsf{P}}(s_{j}-s_{\cal B})}{\sqrt{ \delta}}\right)^{\alpha}\cdot H_{\alpha+\beta}\left(\frac{{\mathsf{P}}(s_{\cal B }-t_{\cal C})}{\sqrt{\delta}}\right)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ B_{\beta}-\frac{(-1)^{\|\beta\|_{1}}}{\beta! }\sum_{j\in{\cal B}}q_{j}\sum_{\alpha>p}\frac{(-1)^{\|\alpha\|_{1}}}{\alpha!} \left(\frac{{\mathsf{P}}(s_{j}-s_{\cal B})}{\sqrt{\delta}}\right)^{\alpha} \cdot H_{\alpha+\beta}\left(\frac{{\mathsf{P}}(s_{\cal B}-t_{\cal C})}{\sqrt{ \delta}}\right)$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\ B_{\beta}+(C_{\beta}-B_{\beta})$ $$

Next, we have

$$ $\displaystyle|\mathrm{Err}_{T}(p)|\leq\left|\sum_{\beta\geq p}B_{\beta}\left( \frac{{\mathsf{P}}(t-t_{\cal C})}{\sqrt{\delta}}\right)^{\beta}\right|+\left| \sum_{\beta\geq p}(C_{\beta}-B_{\beta})\cdot\left(\frac{{\mathsf{P}}(t-t_{\cal C })}{\sqrt{\delta}}\right)^{\beta}\right|$ (19) $$

Using Lemma E.4, we can upper bound the first term in the Eq. (19) by,

$$ $\displaystyle\left|\sum_{\beta\geq p}B_{\beta}\left(\frac{{\mathsf{P}}(t-t_{ \cal C})}{\sqrt{\delta}}\right)^{\beta}\right|\leq$ $\displaystyle\leavevmode\nobreak\ \frac{\sum_{j\in{\cal B}}|q_{j}|}{(1-r)^{w}} \sum_{l=0}^{w-1}{w\choose l}(1-r^{p})^{l}\left(\frac{r^{p}}{\sqrt{p!}}\right)^ {w-l}.$ $$

Since we can similarly bound $C_{\beta}-B_{\beta}$ as follows

$$ $\displaystyle|C_{\beta}-B_{\beta}|\leq\frac{1}{\beta!}K\cdot Q_{B}2^{\|\beta\| _{1}/2}\sqrt{\beta!}\leq KQ_{B}\frac{2^{\|\beta\|_{1}/2}}{\sqrt{\beta!}},$ $$

we have the same bound for the second term

$$ $\displaystyle\leavevmode\nobreak\ \left|\sum_{\beta\geq p}(C_{\beta}-B_{\beta} )\left(\frac{{\mathsf{P}}(t-t_{\cal C})}{\sqrt{\delta}}\right)^{\beta}\right| \leq\leavevmode\nobreak\ \frac{\sum_{j\in{\cal B}}|q_{j}|}{(1-r)^{w}}\sum_{l=0 }^{w-1}{w\choose l}(1-r^{p})^{l}\left(\frac{r^{p}}{\sqrt{p!}}\right)^{w-l}.$ $$

∎

###### Proof of Theorem F.2 .

###### Lemma F.4 (Truncated Hermite expansion with truncated Taylor expansion (low dimension version of Lemma E.5 )) .

###### Proof.

### F.3 Dynamic Low-Rank FGT with Increasing Rank

Figure: Algorithm 10 This algorithm is the update part of Theorem F.5.

Figure: Algorithm 11 This algorithm is another part of Theorem F.5.

We further give an algorithm for FGT when the low-dimensional subspace is dynamic, i.e., the rank may increase with data insertions.

###### Theorem F.5 (Low Rank Dynamic FGT Data Structure) .

Let $W$ be a subspace of ${\mathbb{R}}^{d}$ with dimension $\dim(S)=w\ll d$. Given $N$ source points $s_{1},\dots,s_{N}\in W$ with charges $q_{1},\dots,q_{N}$, and $M$ target points $t_{1},\dots,t_{M}\in W$, a number $\delta>0$, and a vector $q\in{\mathbb{R}}^{N}$,
let $G:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}$ be defined as $G(t)=\sum_{i=1}^{N}q_{i}\cdot\mathsf{K}(s_{i},t)$ denote the *kernel-density* of $t$ with respect to $S$, where $\mathsf{K}(s_{i},t)=f(\|s_{i}-t\|_{2})$ for $f$ satisfying the properties in Definition B.1 .
There is a dynamic data structure
that supports the following operations:

- •
Init$()$ (Algorithm 6) Preprocess in $N\cdot\log^{O(w)}(\|q\|_{1}/\varepsilon)$ time.
- •
KDE-Query$(t\in{\mathbb{R}}^{d})$ (Algorithm 7) Output $\widetilde{G}$ such that $G(t)-\varepsilon\leq\widetilde{G}\leq G(t)+\varepsilon$ in $\log^{O(w)}(\|q\|_{1}/\varepsilon)$ time.
- •
Insert$(s\in{\mathbb{R}}^{d},q_{s}\in{\mathbb{R}})$ (Algorithm 10) For any source point $s\in{\mathbb{R}}^{d}$ and its charge $q_{s}$, update the data structure by adding this source point in $\log^{O(w)}(\|q\|_{1}/\varepsilon)$ time. The subspace dimension $w$ may be increased by 1 if $s$ is not in the original subspace.
- •
Query$(q\in{\mathbb{R}}^{N})$ (Algorithm 7) Output $\widetilde{\mathsf{K}q}\in{\mathbb{R}}^{N}$ such that $\|\widetilde{\mathsf{K}q}-\mathsf{K}q\|_{\infty}\leq\varepsilon$, where $\mathsf{K}\in{\mathbb{R}}^{N\times N}$ is defined by $\mathsf{K}_{i,j}=\mathsf{K}(s_{i},s_{j})$ in $N\log^{O(w)}(\|q\|_{1}/\varepsilon)$ time.

###### Proof.

Since Algorithm 10 updates $A_{\alpha},C_{\beta}$ in the same way as Algorithm 8, the correctness of Procedures KDE-Query and Query follows similarly from Theorem B.5.

Furthermore, Scale takes $O(wd+(N(B)+N(C))\cdot p^{w})$ time. For the correctness, we know that the rows of ${\mathsf{P}}$ form an orthonormal basis for the subspace. For a newly inserted point $s$, if it is not lie in the subspace, $(I-{\mathsf{P}})s$ gives a new basis direction. Therefore, we can easily update ${\mathsf{P}}$ by attaching this vector (after normalization) as a column.
Then, we show that the intermediate variables $A_{\alpha}$ and $C_{\beta}$ can be correctly updated for the new subspace.
For each source box ${\cal B}$ and each $w$-tuple $\alpha\leq p$, we have

$$ $\displaystyle A_{(\alpha,0)}^{\mathsf{new}}({\cal B})=\frac{(-1)^{\|\alpha\|_{ 1}}\cdot(-1)^{i}}{\alpha!\cdot i!}\sum_{j\in{\cal B}}q_{j}\cdot\left(\frac{x^{ \prime}_{j}-x^{\prime}_{\cal B}}{\sqrt{\delta}}\right)^{(\alpha,i)}=A_{\alpha} ({\cal B}),$ $$

where $x^{\prime}_{j}$ denotes the “lifted” point in the new subspace.
And $A_{(\alpha,i)}^{\mathsf{new}}({\cal B})=0$ for all $i>0$, since $(x^{\prime}_{j}-x^{\prime}_{\cal B})_{k+1}=0$.
Similarly, for each target box ${\cal C}$,

$$ $\displaystyle C_{(\beta,i)}^{\mathsf{new}}({\cal C})=$ $\displaystyle\frac{(-1)^{\|\beta\|_{1}}(-1)^{i}}{\beta!i!}\sum_{{\cal B}}\sum_ {\alpha\leq p}\sum_{j=0}^{p}A_{(\alpha,j)}^{\mathsf{new}}({\cal B})H_{(\alpha+ \beta,i+j)}\left(\frac{x^{\prime}_{{\cal B}}-y^{\prime}_{{\cal C}}}{\sqrt{ \delta}}\right)$ $\displaystyle=$ $\displaystyle\frac{(-1)^{\|\beta\|_{1}}(-1)^{i}}{\beta!i!}\sum_{{\cal B}}\sum_ {\alpha\leq p}A_{\alpha}({\cal B})H_{\alpha+\beta}\left(\frac{x_{{\cal B}}-y_{ {\cal C}}}{\sqrt{\delta}}\right)\cdot h_{i}(0)$ $\displaystyle=$ $\displaystyle\frac{(-1)^{i}}{i!}h_{i}(0)\cdot C_{\beta}({\cal C}),$ $$

where the second step follows from $A_{(\alpha,i)}^{\mathsf{new}}({\cal B})=A_{\alpha}({\cal B})\cdot{\bf 1}_{i=0}$.
Therefore, by enumerating all boxes ${\cal B},{\cal C}$ and indices $\alpha,\beta\leq p$, we can correctly compute $A_{(\alpha,0)}^{\mathsf{new}}({\cal B})$ and $C_{(\beta,i)}^{\mathsf{new}}({\cal C})$. Thus, we complete the proof of the correctness of Algorithm 11.
∎

###### Theorem F.5 (Low Rank Dynamic FGT Data Structure) .

###### Proof.

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 2
questions: Please refer to the weakness.
soundness: 3
strengths: 1. This paper has a strong theoretical contribution by solving an open problem of dynamizing the FGT algorithm, which has applications in machine learning and data analysis.
2. This paper includes a novel technical approach combining Taylor and Hermite expansions to enable efficient updates, with careful error analysis and proofs.
3. This paper improves complexity bounds when data lies in low-dimensional subspaces, addressing the curse of dimensionality issue in previous FGT work.
summary: This paper presents a dynamic data structure for maintaining the Fast Gaussian Transform (FGT) under insertions and deletions of source points. The main contribution is an algorithm that supports polylogarithmic-time updates and kernel density estimation queries when the data points lie in a low-dimensional subspace while retaining quasi-linear time for matrix-vector multiplication. The key technical innovation is a hybrid approach combining Taylor and Hermite expansions to "disentangle" source and target points, enabling efficient dynamic updates.
weaknesses: 1. This paper has limited experimental evaluations, the paper focuses primarily on theoretical analysis without empirical validation of real datasets.
2. The assumptions of data points lie in a low-dimensional subspace, which may not hold in practice for many applications.
3. The paper has limited discussion of numerical stability issues that may arise in practice when computing Taylor/Hermite expansions.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 4
questions: See weakness section on how to extend the results to other types of kernels? how to handle tuning parameters as mentioned in the weakness section?
soundness: 3
strengths: The key contribution is the new methodology proposed by this paper to handle efficiently evolving datasets in kernel density estimation (KDE) and matrix-vector multiplication in a dynamic setting. This question itself is very important in practice and the new method has good time complexity that supports fast updates and KDE queries in sublinear time, with quasi-linear time complexity for full matrix-vector multiplication. Overall, the paper’s strengths lie in its innovative adaptation of FGT to a dynamic setting, its computational efficiency, and its broad relevance to various machine learning and statistical applications. The work addresses a crucial need for scalable kernel-based methods in dynamic, real-time data environments.
summary: The paper introduces a Dynamic Low-Rank Fast Gaussian Transform (FGT), designed to handle efficiently evolving datasets in kernel density estimation (KDE) and matrix-vector multiplication. The Fast Gaussian Transform (FGT) traditionally enables subquadratic-time multiplication of an n×n Gaussian kernel matrix with a vector, but it struggles with dynamically changing data. The proposed dynamic FGT algorithm supports fast updates and KDE queries in sublinear time, with quasi-linear time complexity for full matrix-vector multiplication.

The paper proposes a new data structure that dynamically maintains low-rank subspaces, supporting efficient insertion and deletion of points, along with fast KDE queries. The method builds on classical fast multipole methods (FMM) and uses Taylor and Hermite expansions to manage “interaction ranks” between source and target boxes. By focusing on data in low-dimensional subspaces, the method bypasses the curse of dimensionality, achieving sublinear update and query times under certain conditions. This is particularly useful for large, evolving datasets common in fields like machine learning. Importantly, they answered a question in an affirmative way:

Is it possible to ‘dynamize’ the Fast Gaussian Transform, in sublinear time? Can the exponentialdependence on d (Greengard & Strain, 1991) be mitigated if the data-points $x_i$ lie in a k-dimensional subspace of R^d?
weaknesses: I think one problem or actually a question that is worth thinking and needs to be clarified is that the method is mainly tailored specifically for Gaussian kernels. Although in Section 3.4, it is generalized to a broader class of "fast-decaying kernels". It is still unclear if it allows other important kernels such as polynomial kernels, Moreover, as I see the kernel is defined using the Euclidean distance, maybe in manifold settings, it is restrictive and not the proper kernel choice.

The other issue is that the paper does not provide detailed guidance on choosing parameters, for example, in Section 3.1, the kernel bandwidth $\delta$ and the length of the smaller boxes $L$, as it is mainly a theoretical paper without real experiments.

While the paper lists theoretical contributions, since the paper proposes a Fast Gaussian Transform (FGT), which is an algorithm on the issue of fast computation, I think it would be good to do some, at least simple simulations on a dataset that is dynamically changing (as claimed by the paper) to see its empirical performance.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
questions: See above.
soundness: 3
strengths: I think this paper is well-written. The introduction clearly outlines the history of FMM and FGT, and states their drawback clearly: they are static and it is hard to add/delete points. I also think the main result section is also well-written, and state the contributions clearly.

I also think the technical contributions are solid: the problem of efficiently updating FGT for changing datasets is important, and this work seems to be one of the first to address this problem.
summary: The fast gaussian transform (FGT) can be used to efficiently apply a gaussian kernal matrix to a vector. However traditionally the number of points in the gaussian kernal matrix is already fixed, so even adding one single point can require us to recompute the FGT from scratch, which makes it very hard to deal with datasets that are dynamically changing. In this work, the authors propose a dynamic FGT algorithm that supports adding or deleting a source point, and also estimating the kernel-density of a query point, under the assumption that the added data-points lie in a low dimensional subspace of the original vector space.
weaknesses: I do not see any major weaknesses, but I think it would be better if the authors could include some numerical experiments to showcase their results. In particular, how much faster is the dynamic FGT compared with recomputing from scratch?

---

confidence: 5
contribution: 2
initial_rating: 3
presentation: 2
questions: .
soundness: 2
strengths: The problem seems well-motivated.
summary: The paper studies the following problem.
Given $n$ points in $\mathbb{R}^d$ lying in a $k$ dimensional subspace, we would like to construct a data structure that supports insertion, deletion and query.
More precisely, we can insert a new point into the data structure or delete an existing point from the data structure.
Also, we can compute the kernel desnity estimation of the points in the data structure at a certain point within $\epsilon$ error efficiently.
The result relies on fast Gauss transform.
The authors constructure a data structure such that the running time of these operations is at most $\log^{O(k)}(1/\epsilon)$ time.
weaknesses: The result relies heavily on the previous results. 
I am not sure if there is a new idea introduced.

Meta-review: The paper address the problem of kernel computation in sub quadratic time via fast Gaussian transofm on Matrix-Vec computation. The proposed method incorporate dynamization where the data source can be added and deleted, while the dimension of the query points can also vary. Despite the interesting setting, from the discussion, the novelty of the paper is limited to appear in the top conference such as ICLR.

Decision: Reject