# StreaMRAK a Streaming Multi-Resolution Adaptive Kernel Algorithm

Andreas Oslandsbotn111University of Oslo, Problemveien 7, 0315 Oslo, Norway
222Simula School of Research and Innovation, Martin Linges Vei 25, 1364
Fornebu, Norway Željko Kereta333University College London, Gower St, WC1E 6BT
London, England Valeriya Naumova444Simula Research Laboratory, Martin Linges
vei 25, 1364 Fornebu, Norway Yoav Freund555University of California San Diego,
9500 Gilman Dr, La Jolla, CA 92093, United States Alexander
Cloninger††footnotemark:

###### Abstract

Kernel ridge regression (KRR) is a popular scheme for non-linear non-
parametric learning. However, existing implementations of KRR require that all
the data is stored in the main memory, which severely limits the use of KRR in
contexts where data size far exceeds the memory size. Such applications are
increasingly common in data mining, bioinformatics, and control. A powerful
paradigm for computing on data sets that are too large for memory is the
streaming model of computation, where we process one data sample at a time,
discarding each sample before moving on to the next one.

In this paper, we propose StreaMRAK \- a streaming version of KRR. StreaMRAK
improves on existing KRR schemes by dividing the problem into several levels
of resolution, which allows continual refinement to the predictions. The
algorithm reduces the memory requirement by continuously and efficiently
integrating new samples into the training model. With a novel sub-sampling
scheme, StreaMRAK reduces memory and computational complexities by creating a
sketch of the original data, where the sub-sampling density is adapted to the
bandwidth of the kernel and the local dimensionality of the data.

We present a showcase study on two synthetic problems and the prediction of
the trajectory of a double pendulum. The results show that the proposed
algorithm is fast and accurate.

##  1 Introduction

Machine learning algorithms based on kernel ridge regression (KRR) [1] is an
active field of research [2, 3, 4, 5, 6], with applications ranging from time
series prediction in finance [7], parameter inference in dynamical systems
[8], to pairwise learning [9], face recognition [10] and drug estimation and
gene analysis in biomedicine [11, 12]. This paper develops a streaming
variation of KRR using a radial kernel, a new sub-sampling scheme, and a
multi-resolution formulation of the learning model.

Many popular data analysis software packages, such as
$\mbox{Matlab}^{\mbox{\tiny TM}}$ require loading the entire dataset into
memory. While the size of computer memory is growing fast, the size of
available data sets is growing much faster, limiting the applicability of in-
memory methods. 666$\mbox{Simulink}^{\mbox{\tiny TM}}$, a companion software
to $\mbox{Matlab}^{\mbox{\tiny TM}}$ supports streaming but has a much more
limited computational model, targeted at signal processing applications.

Streaming [13] is a computational model where the input size is much larger
than the size of memory. Streaming algorithms read one input item at a time,
update their memory, and discard the item. The computer memory is used to
store a model or a sketch of the overall data distribution, which is orders of
magnitude smaller than the data itself. The development of streaming
algorithms is experiencing increased popularity in the face of big data
applications such as data mining [14] and bioinformatics [15], where data sets
are typically too large to be kept in-memory. Many big data applications call
for non-linear and involved models, and thus, the development of non-
parametric and non-linear models is critical for successful learning.

Among the most popular non-parametric learning algorithms are kernel methods,
which include well-known learning schemes such as the support vector machine
(SVM) and KRR, to name a few. The appeal of kernel methods lies in their
strong theoretical foundation [1, 16], as well as their ability to map complex
problems to a linear space without requiring an explicit mapping. A common
class of kernels are radial kernels
$k({\bm{\mathsf{x}}},\tilde{\bm{\mathsf{x}}})=\Phi(\|{\bm{\mathsf{x}}}-\tilde{\bm{\mathsf{x}}}\|/r)$
for ${\bm{\mathsf{x}}},\tilde{\bm{\mathsf{x}}}\in{\cal
X}\subseteq\mathbb{R}^{D}$ and $r>0$ [17]. An example is a Gaussian kernel,
for which the shape parameter $r$ is referred to as the kernel bandwidth. What
is more, radial kernels are universal kernels (with a few exceptions [18]),
meaning that they can approximate any bounded continuous function on ${\cal
X}$ arbitrarily well. However, in high dimensions kernel methods suffer from
the ”curse of dimensionality” and require large amounts of training data to
converge. Furthermore, the computational complexity, memory requirement, and
the number of parameters to learn grow unbounded with the number of training
samples, a drawback known as the ”curse of kernelization” [19]. In the context
of streaming, the prospect of unbounded data streams makes this shortcoming
even more detrimental.

Although kernel-based learning schemes are typically formulated as convex
optimization problems, which do not require tuning hyper-parameters such as
learning rate etc., there is still a need to determine the optimal kernel. For
the Gaussian kernel, this amounts to selecting the bandwidth. Classically, an
optimal kernel is chosen through batch techniques such as leave-one-out and
k-fold cross-validation [20, 21, 22]. However, these approaches are
inefficient as they spend significant time evaluating bad kernel hypotheses
and often use multiple runs over the data, which is impossible in a streaming
setting.

Despite the universality of radial kernels on ${\cal X}\in\mathbb{R}^{D}$,
this only guarantees the convergence of the model in the asymptotic regime and
does not provide finite sample bounds. As a reaction, several works have shown
the benefit of combining multiple kernels from a dictionary of kernel
hypotheses. These strategies include multi-kernel learning (MKL) [4, 23, 24,
25, 26], multi-scale analysis [27, 28], and the Laplacian pyramid (LP) [29,
30]. Combining these strategies with a localized kernel gives a frequency and
location-based discretization similar to multi-resolution analysis, a well-
established concept in signal processing and functional approximation through
concepts such as wavelets [31, 32], diffusion wavelets [33, 34], and graph
wavelets [35, 36, 37].

To meet a need for non-linear non-parametric algorithms for streaming data, we
propose the streaming multi-resolution adaptive kernel algorithm (StreaMRAK) -
a computationally and memory-efficient streaming variation of KRR. StreaMRAK
is a streaming algorithm that combines a streaming sub-sampling scheme with a
multi-resolution kernel selection strategy and adapts the kernel bandwidth $r$
and the sub-sample density to each other over several levels of resolution.
Furthermore, StreaMRAK addresses the curse of dimensionality in a novel way,
through the sub-sampling scheme and multi-resolution formulation.

###  1.1 Setting

We consider a finite sample data-cloud ${\cal X}$, $|{\cal X}|=n$, that is
sampled i.i.d. according to a fixed but unknown distribution ${\cal P}$ over
$\mathbb{R}^{D}$. The target is a bounded and continuous function
$f:\mathbb{R}^{D}\to\mathbb{R}$. We assume that the points in ${\cal X}$ are
placed in a sequence and that their order is random. 777The assumption that
the sequence is randomly ordered allows us to draw statistical conclusions
from prefixes. Each instance ${\bm{\mathsf{x}}}_{i}\in{\cal X}$, for
$i\in[n]$, paired with a label $y_{i}$ where
$y_{i}=f({\bm{\mathsf{x}}}_{i})+\varepsilon_{i}$ and $\varepsilon_{i}\sim{\cal
N}(0,\sigma)$ represents noise. The task of learning is to train a model
$\widehat{f}$ that is a good approximation of the target function $f$.

In this work, we think about the intrinsic dimension of ${\cal X}$ as a local
quantity, meaning it depends on the region ${\cal A}\subseteq{\cal X}$ and the
radius $r$ at which we consider the point cloud. Rooted in this way of
thinking about the intrinsic dimension, StreaMRAK is designed to handle
domains where the local intrinsic dimension changes across different regions
and resolutions.

To estimate the local intrinsic dimension in a ”location and resolution
sensitive” manner, we use the concept of the doubling dimension of a set,
defined in Def. 1.2. We note that our definition of the doubling dimension is
related to the definition used in [38, 39].

######  Definition 1.1 (Covering number).

Consider a set ${\cal A}$ and a ball ${\cal B}({\bm{\mathsf{x}}},r)$, with
$r>0$ and ${\bm{\mathsf{x}}}\in{\cal A}$. We say that a finite set ${\cal
S}\subset{\cal B}(x,r)$ is a covering of ${\cal B}(x,r)$ in ${\cal A}$ if
${\cal A}\cap{\cal B}(x,r)\subset\cup_{{\bm{\mathsf{x}}}_{i}\in{\cal S}}{\cal
B}({\bm{\mathsf{x}}}_{i},r/2)$. We define the covering number $\kappa({\cal
A},{\bm{\mathsf{x}}},r)$ as the minimum cardinality of any covering of ${\cal
B}({\bm{\mathsf{x}}},r)$ in ${\cal A}$.

######  Definition 1.2.

(Doubling dimension) The doubling dimension $\texttt{ddim}({\cal A},r)$ of a
set ${\cal A}$ is defined as $\texttt{ddim}({\cal A},r)=\lceil\log\kappa({\cal
A},{\bm{\mathsf{x}}},r)\rceil$. For an interval ${\cal
I}\subset\mathbb{R}_{>0}$ we define the doubling dimension as the least upper
bound over $r\in{\cal I}$, that is $\texttt{ddim}({\cal A},{\cal
I})=\max_{r\in{\cal I}}\texttt{ddim}({\cal A},r)$.

Using Def. 1.2 we say that the intrinsic dimension of ${\cal X}$ changes with
the location if there exist ${\cal A}_{1},{\cal A}_{2}\subset{\cal X}$ such
that $\texttt{ddim}({\cal A}_{1},r)\neq\texttt{ddim}({\cal A}_{2},r)$ for some
$r>0$. Similarly, we say that the intrinsic dimensionality of ${\cal X}$
changes with the resolution, if there exist $r_{1}\neq r_{2}$ such that the
doubling dimension $\texttt{ddim}({\cal A},r_{1})\neq\texttt{ddim}({\cal
A},r_{2})$ for ${\cal A}\subset{\cal X}$.

In Fig. 1 we consider three examples to provide further insight for the
doubling dimension. In Fig. LABEL:subfig:change_loc we see a domain shaped
like a dumbbell, where the spheres are high dimensional, and the bar
connecting them is lower-dimensional, showing how the dimension can change
with the location. Meanwhile, Fig. LABEL:subfig:change_noise illustrates a
lower-dimensional manifold, embedded in $\mathbb{R}^{3}$, with manifold noise
$\zeta_{m}$. We see that when the resolution is sufficiently small, so that
$r\approx\zeta_{m}$, the doubling dimensionality increases towards the
dimension of the ambient space $\mathbb{R}^{D}$. Furthermore, Fig.
LABEL:subfig:change_scale shows a point cloud that is locally $2$-dimensional,
but is embedded in a $3$-dimensional space. By reducing $r$ we can resolve
this lower dimensionality, but if it is reduced further, we would eventually
resolve the noise level, and the doubling dimension increases again.

We also mention two special cases. First, for large enough $r$ any set in
$\mathbb{R}^{D}$ has a doubling dimension of at most $D$. Second, if ${\cal
A}$ is a finite set of points in $\mathbb{R}^{D}$ and $r$ is smaller than the
minimal distance between two points, then the number of balls of radius
$r^{\prime}<r$ required to cover ${\cal A}$ is at most the number of points.
Therefore the doubling dimension of ${\cal A}$ at the range $(0,r]$ is zero.
In other words, any actual (and therefore finite) training set has dimension
zero for a small enough $r$, as illustrated in Fig. LABEL:subfig:change_scale.

(a)

(b)

(c)

Figure 1: Three examples of variation in the intrinsic dimension. The coloring
of the point clouds illustrates depth.

As an example of how the intrinsic dimension might change with respect to
regions and resolutions, we consider a double pendulum system, a well-known
chaotic system that depends heavily on its initial conditions [40]. Systems
with multiple pendulum elements are well known in engineering applications
such as mechanical and robotic systems with several joints and are studied for
their chaotic properties [41].

Let $\omega_{i}=\dot{\theta_{i}}$, for $i=1,2$, be angular velocities. We
initialize $500$ pendulums with $\theta_{1},\theta_{2}=180$,
$\omega_{1}=10^{-2}$ and $\omega_{2}=10^{-1}$, measured in degrees, and
perturb the angular velocities by $\varepsilon_{i}\sim{\cal
N}(0,10^{-2}|\omega_{i}|)$. We iterate the system for $T=500$ time steps and
let
${\bm{\mathsf{s}}}^{(i)}_{t}=[\theta_{1}(t),\,\theta_{2}(t),\,\omega_{1}(t),\,\omega_{2}(t)]\in\mathbb{R}^{4}$
be the state of pendulum $i$ at $t\in\mathbb{N}$. We think of each state as a
training point in a point cloud ${\cal X}\subset\mathbb{R}^{4}$, where for
instance the target function can be
${\bm{\mathsf{s}}}^{(i)}_{t+\Delta}=f({\bm{\mathsf{s}}}^{(i)}_{t})$, with
$\Delta\in\mathbb{N}$.

In Fig. LABEL:subfig:phaseDiagramIntroduction_2d we visualize the trajectory
of four pendulums $P_{0},P_{1},P_{2},P_{3}$, for which the trajectories are
indistinguishable until a bifurcation occurs around $T=300$ time steps, and
the trajectories start to diverge. In Fig.
LABEL:subfig:phaseDiagramIntroduction_3d_before_bifurcation and Fig.
LABEL:subfig:phaseDiagramIntroduction_3d we zoom in on the trajectory of all
$500$ pendulums in regions before and after the bifurcation. These two
regions, ${\cal A}_{1}$ and ${\cal A}_{2}$, are indicated by a blue and red
circle, respectively, in Fig. LABEL:subfig:phaseDiagramIntroduction_2d. From
the figures, it is clear that learning the trajectory in ${\cal A}_{1}$ is
significantly easier than in ${\cal A}_{2}$, where learning the trajectory is
more affected by the curse of dimensionality.

We also note that the trajectories that remain close after the bifurcation
will remain so until a new bifurcation occurs. The take-home message is that
predicting the trajectory of a double pendulum is a hard problem because of a
few regions where the intrinsic dimension blows up and makes the prediction
hard. However, between these regions, the trajectory is easier to describe.
The spirit of this work aims to reduce the effort in such regions ${\cal
A}_{2}$ where the training data exhibit high intrinsic dimension and focus
more on those regions ${\cal A}_{1}$ where the data has a lower intrinsic
dimension, i.e. is more well behaved.

(a)

(b)

(c)

(d)

Figure 2: (a) Illustration of a double pendulum. Here $l$ and $m$ are the
length and mass of the pendulum rods, and $\theta_{1},\theta_{2}$ are the
angles. Furthermore, $(x_{1},y_{1})$ and $(x_{2},\,y_{2})$ are the positions
of the point masses of the two pendulums. (b) Phase diagram of four double
pendulums $P_{0},P_{1},P_{2},P_{3}$, iterated for $T=500$ time steps. The
bifurcation point at step $T=300$ is indicated with a black solid circle. (c)
and (d) includes the $\omega_{1}$ axis and zoom in on respectively the blue
and red circles in (b).

###  1.2 Contribution and comparison to related work

Contributions of this work can be divided into three components.

  1. (C1)

A multi-resolution variation of the state-of-the-art KRR solver FALKON [2],
using the LP, which refines the predictions at each level of resolution by
regressing on the errors from the previous level.

  2. (C2)

A novel sub-sampling scheme for kernel methods, tailored for use in
combination with the LP, that can handle the curse of dimensionality and does
not require the data to be in-memory.

  3. (C3)

Development of a streaming variation of FALKON, where the time and memory
requirements depend on the doubling dimensionality and the level of
resolution, instead of the number of training points. see Props. 5.1-5.3.

In the following, we give further details on these contributions and compare
them to related work.

The computational backbone of StreaMRAK is based on the state-of-the-art KRR
solver FALKON [2], which among other things combines sub-sampling and
preconditioning to process large data sets efficiently. However, FALKON relies
on selecting an optimal kernel bandwidth, which can be inefficient and hard
within a streaming setting.

Inspired by the success of existing multi-resolution approaches [4, 23, 24,
25, 26, 27], our first contribution (C1) addresses the issue of selecting an
optimal kernel bandwidth by introducing a multi-resolution reformulation of
FALKON using a variation of the LP.

The LP scheme originated in image representation [42] and was introduced to
machine learning by [29] for efficient data representation. The LP refines the
prediction at several levels of resolution, and at each level, reduces the
bandwidth used at the previous level by a constant factor. In doing so, the
selection of a single optimal bandwidth is avoided, and the resulting approach
has greater flexibility. The LP is similar to ideas in wavelet analysis that
have shown great success in numerous applications. However, typical wavelet
architectures [31, 32, 33, 34, 35, 36] require upfront construction of a
wavelet basis, which is not compatible with a data-adaptive kernel. In this
work, we aim to show that the LP is a viable multi-resolution scheme and can
be modified to the streaming setting. Furthermore, we experimentally show that
it significantly improves the estimation accuracy, and inspired by [30] we
provide convergence bounds for the LP in the context of radial kernels and
KRR.

Let us now discuss our second contribution (C2). FALKON addresses the curse of
kernelization by combining Nyström sub-sampling, conjugate gradient, and
preconditioning, and achieves time and memory requirements of
$\mathcal{O}(n\sqrt{n})$ and $\mathcal{O}(n)$ respectively, where $n$ is the
number of samples. In recent years there have been several efforts to address
the curse of kernelization in similar ways through sub-sampling techniques
such as sketching [3, 5], randomized features [43, 44, 45, 46] and Nyström
sub-sampling [47, 48, 49, 50, 51]. However, despite their successes, these
techniques are in principle in-memory type algorithms since they require
access to the training data in advance of the training and are not optimized
for streaming.

Furthermore, FALKON selects the sub-samples uniformly over the input domain
${\cal X}$, and the LP uses the same training set for each level. However,
when learning with a radial kernel, the density of samples should be related
to the bandwidth of the kernel. Otherwise, a too-small bandwidth will lead to
bad interpolation properties, while a too-large bandwidth gives an ill-
conditioned system [27]. Since the LP scheme reduces the kernel bandwidth at
each level of resolution, it would be problematic to use the same sub-sample
density. Furthermore, due to the curse of dimensionality, the covering number
increases exponentially with the doubling dimension. Therefore, if doubling
dimensionality varies across different regions of the domain ${\cal X}$, as
illustrated by Fig. LABEL:subfig:change_loc, then the number of sub-samples
necessary to maintain the density for a given bandwidth will also vary.

Our second contribution (C2) provides an alternative sub-sampling strategy,
which adapts the sub-sampling density to the kernel bandwidth. This strategy
is similar to tuning the kernel bandwidth to the data, which is used in online
algorithms to avoid the use of cross-validation [52, 53, 54]. Although
expensive to calculate, especially in high dimensions, a similar strategy is
used in graph-based methods, where the kernel bandwidth is adapted to the
$k$-nearest neighbor distance of each training point [55].

The sub-sampling strategy we propose is based on a damped cover-tree (DCT),
which is a modified version of the cover-tree (CT) [39]. The CT is a tree-
based data structure originally intended for nearest neighbor search in metric
spaces. It is closely related to navigation nets [38], but with improved space
requirements: $\mathcal{O}(n)$ in memory and $\mathcal{O}(c^{6}n\log n)$ in
time. In this work, we show how the CT structure can be used to organize the
samples hierarchically with increasing density for each new level in the tree
and how it can adapt the sub-sample density to kernel bandwidths in the LP.

However, the problem with an adaptive sub-sampling strategy is its
vulnerability to the curse of dimensionality. In regions of high doubling
dimensions, the number of samples to achieve a certain density increases
exponentially, as quantified by Def. 1.2. This means that the number of sub-
samples from the CT will quickly grow too large for efficient computing. The
danger is to waste resources on samples from subsets and levels where the
doubling dimension is so large that good interpolation cannot be achieved for
any viable sample sizes. This would only serve to slow down the computation
and not increase the precision.

Due to this, the DCT introduces a damping property, which gradually suppresses
the selection of sub-samples where the doubling dimensionality is large. This
has the additional advantage of allowing to choose more sub-samples from
regions where the doubling dimensionality is small. Thus, the DCT can diminish
the impact of the curse of dimensionality. Furthermore, the DCT can be built
continuously as new samples come in, making it ideal for a streaming
computational model.

Our third contribution (C3), relies on the changes implemented with (C1) and
(C2). In particular, StreaMRAK can operate as a streaming algorithm and
efficiently organizes the sub-samples as it builds the multi-resolution
kernel. Furthermore, the sub-sampling and kernel construction allows for
continuous integration of new training points into the kernel matrix.
Moreover, the DCT, the multi-resolution construction, and the KRR solver can
all be multi-threaded and parallelized.

###  1.3 Organization of the paper

The paper is organized as follows. Section 2 introduces kernel methods and the
FALKON algorithm, as well as the LP. Section 3 introduces the adaptive sub-
sampling scheme and the DCT. StreaMRAK is described in Section 4 and an
analysis of the algorithm is given in Section 5. Finally, Section 6 presents
several numerical experiments and Section 7 gives an outlook for further work.
The Appendix includes further mathematical background and the proofs.

###  1.4 Notation

We denote vectors ${\bm{\mathsf{a}}}\in\mathbb{R}^{D}$ with boldface and
matrices ${\mathbf{A}}\in\mathbb{R}^{n\times m}$ with bold uppercase, and
${\mathbf{A}}^{\top}$ denotes the matrix transpose. We use ${\mathbf{K}}_{nm}$
for kernel matrices, where the dimensionality is indicated by the subscripts.
We reserve $n$ for the number of training samples and $m$ for the number of
sub-samples. The $ij$-th element of a kernel matrix is denoted
$[{\mathbf{K}}_{nm}]_{ij}$, while for other matrices we use
${\mathbf{A}}_{ij}$. The notation $a_{i}$ indicates $i$-th element of a vector
${\bm{\mathsf{a}}}$. Furthermore, we use $f([{\bm{\mathsf{x}}}_{n}])$ to
denote
$(f({\bm{\mathsf{x}}}_{1}),\ldots,f({\bm{\mathsf{x}}}_{n}))^{\top}\in\mathbb{R}^{n}$,
and $[m]$ to denote $\\{i\\}^{m}_{i=1}$. The notation ${\bm{\mathsf{x}}}_{i}$
indicates the $i$-th training example. We use ${\bm{\mathsf{a}}}^{(l)}$ and
${\mathbf{A}}^{(l)}$, where $l$ refers to a specific level in the LP and the
DCT. We take $\|\cdot\|$ to be the $L^{2}$ norm and $\|\cdot\|_{{\cal H}}$ to
be the RKHS norm. We denote the intrinsic dimension of a manifold with $d$ and
the dimension of the embedding with $D$. By $\mathbbm{1}_{{\cal
S}}({\bm{\mathsf{x}}})$ we denote the indicator function, which evaluates to
$1$ if ${\bm{\mathsf{x}}}\in{\cal S}$ and $0$ otherwise, of a set ${\cal
S}\subset\mathbb{R}^{D}$.

##  2 Kernel methods

Consider a positive definite kernel $k:{\cal X}\times{\cal
X}\rightarrow\mathbb{R}$, defined on an input space ${\cal
X}\subset\mathbb{R}^{D}$. Given data
$\\{({\bm{\mathsf{x}}}_{i},y_{i}):i\in[n]\\}$ of samples from ${\cal
X}\times\mathbb{R}^{D}$, kernel ridge regression computes an estimator by
minimising

$\displaystyle\widehat{f}_{n,\lambda}=\operatorname*{argmin}_{f\in\widetilde{\cal
H}}\frac{1}{n}\sum_{i=1}^{n}(f({\bm{\mathsf{x}}}_{i})-y_{i})^{2}+\lambda\left\|{f}\right\|_{\cal
H}^{2},$

where ${\cal H}$ is the Hilbert space induced by the kernel. This allows to
reduce the problem to a linear system

$\displaystyle({\mathbf{K}}_{nn}+\lambda
n{\mathbf{I}}_{n}){\bm{\upalpha}}={\bm{\mathsf{y}}},\text{ for
}[{\mathbf{K}}_{nn}]_{ij}=k({\bm{\mathsf{x}}}_{i},{\bm{\mathsf{x}}}_{j}),\text{
and }{\bm{\mathsf{y}}}=(y_{1},\ldots,y_{n})^{\top}.$ (2.1)

Coefficients ${\bm{\upalpha}}=(\alpha_{1},\ldots,\alpha_{n})^{\top}$ define
the estimator by
$f({\bm{\mathsf{x}}})=\sum_{i=1}^{n}\alpha_{i}k({\bm{\mathsf{x}}},{\bm{\mathsf{x}}}_{i})$.
However, solving (2.1) using traditional methods has a time complexity of
${\cal O}(n^{2})$, which can be costly for large $n$ [2].

FALKON [2] addresses this issue by sub-sampling the columns of
${\mathbf{K}}_{nn}$, which reduces the effective complexity while maintaining
accuracy. Namely, denote
$\Gamma_{n}=\\{{\bm{\mathsf{x}}}_{1},\ldots,{\bm{\mathsf{x}}}_{n}\\}$ and for
$m\ll n$ let
$\widetilde{\Gamma}_{m}=\\{\widetilde{{\bm{\mathsf{x}}}}_{1},\ldots,\widetilde{{\bm{\mathsf{x}}}}_{m}\\}$
be Nyström centers (i.e. a randomly selected subset of $\Gamma_{n}$).
Minimizing

$\widehat{f}_{n,m,\lambda}=\operatorname*{argmin}_{f\in\widetilde{\cal
H}_{M}}\frac{1}{n}\sum_{i=1}^{n}(f({\bm{\mathsf{x}}}_{i})-y_{i})^{2}+\lambda\left\|{f}\right\|_{\cal
H}^{2},$ (2.2)

where $\widetilde{\cal
H}_{m}=\overline{\operatorname{span}}\left\\{{k(\cdot,\widetilde{{\bm{\mathsf{x}}}}_{j}):j\in[m]}\right\\}$,
leads to a linear system

$\displaystyle{\mathbf{H}}\widetilde{{\bm{\upalpha}}}={\bm{\mathsf{z}}},\text{
for }{\mathbf{H}}={\mathbf{K}}_{nm}^{\top}{\mathbf{K}}_{nm}+\lambda
n{\mathbf{K}}_{mm},\text{ and
}{\bm{\mathsf{z}}}={\mathbf{K}}_{nm}{\bm{\mathsf{y}}}.$

Here
$[{\mathbf{K}}_{nm}]_{ij}=k({\bm{\mathsf{x}}}_{i},\widetilde{{\bm{\mathsf{x}}}}_{j})\in\mathbb{R}^{n\times
m}$ is the column-subsampled matrix and the estimator is given by
$\widehat{f}_{n,m,\lambda}({\bm{\mathsf{x}}})=\sum_{j=1}^{m}\widetilde{\alpha}_{j}k({\bm{\mathsf{x}}},\widetilde{{\bm{\mathsf{x}}}}_{j})$.
To further reduce the time complexity FALKON uses a suitable preconditioner to
reduce the condition number. The preconditioner is defined as
${\mathbf{B}}{\mathbf{B}}^{\top}=(n/m{\mathbf{K}}_{mm}^{2}+\lambda
n{\mathbf{K}}_{mm})^{-1}$, which is a natural (lower complexity) approximation
of the ideal preconditioner
${\mathbf{A}}{\mathbf{A}}^{\top}=({\mathbf{K}}_{nm}^{\top}{\mathbf{K}}_{nm}+\lambda
n{\mathbf{K}}_{mm})^{-1}$. We now solve for $\widetilde{{\bm{\upalpha}}}$ from
the system of equations

${\mathbf{B}}^{\top}{\mathbf{H}}{\mathbf{B}}{\bm{\upbeta}}={\mathbf{B}}^{\top}{\bm{\mathsf{z}}},\text{
for }{\mathbf{H}}={\mathbf{K}}_{nm}^{\top}{\mathbf{K}}_{nm}+\lambda
n{\mathbf{K}}_{mm},\text{
}{\bm{\mathsf{z}}}={\mathbf{K}}_{nm}{\bm{\mathsf{y}}},\text{ and
}\widetilde{{\bm{\upalpha}}}={\mathbf{B}}{\bm{\upbeta}}.$ (2.3)

This is solved iteratively, using the conjugate gradients with early stopping.
Choosing $m={\cal O}(\sqrt{n})$ still ensures optimal generalisation (i.e.
same as KRR), while reducing the computational complexity to ${\cal
O}(n\sqrt{n})$.

###  2.1 Streaming adaptation of FALKON

Matrices and vectors involved in the linear system in (2.3) can be separated
into two classes: those that depend only on sub-samples in
$\widetilde{\Gamma}_{m}$; and those
(${\mathbf{K}}_{nm}^{\top}{\mathbf{K}}_{nm}$ and ${\bm{\mathsf{z}}}$) that
also depend on all the training points $\Gamma_{n}$. Critically, terms in both
groups are all of size $m$, which allows to reduce the complexity. Consider
now the set of sub-samples $\widetilde{\Gamma}_{m}$ to be fixed, and assume
new training points, in the form
$\\{({\bm{\mathsf{x}}}_{q},y_{q}):q=n+1,\ldots,n+t)\\}$, are coming in a
stream. We can then update the second class of terms according to

$\displaystyle\big{[}({\mathbf{K}}_{(n+t)m})^{\top}{\mathbf{K}}_{(n+t)m}\big{]}_{ij}$
$\displaystyle=\big{[}({\mathbf{K}}_{nm})^{\top}{\mathbf{K}}_{nm}\big{]}_{ij}+\sum_{q=n+1}^{n+t}k({\bm{\mathsf{x}}}_{q},\widetilde{{\bm{\mathsf{x}}}}_{i})k({\bm{\mathsf{x}}}_{q},\widetilde{{\bm{\mathsf{x}}}}_{j}),$
(2.4)
$\displaystyle\big{[}({\mathbf{K}}_{(n+t)m})^{\top}{\bm{\mathsf{y}}}\big{]}_{i}$
$\displaystyle=z_{i}+\sum_{q=n+1}^{n+t}k({\bm{\mathsf{x}}}_{q},\widetilde{{\bm{\mathsf{x}}}}_{i})y_{q}.$
(2.5)

Thus, only sub-samples $\widetilde{\Gamma}_{m}$, matrices
$\big{(}{\mathbf{K}}_{nm}\big{)}^{\top}{\mathbf{K}}_{nm}$, ${\mathbf{K}}_{mm}$
and ${\bm{\mathsf{z}}}$, need to be stored. However, in order to continuously
incorporate new training points into Eqs. (2.4) and (2.5), sub-samples
$\widetilde{\Gamma}_{m}$ must be determined in advance. Whereas this works if
all the data is provided beforehand, it cannot be done if the data arrives
sequentially. In this work, we address this through a multi-resolution
framework. The overall estimator is composed of a sequence of estimators
defined at different resolution levels of the domain. Correspondingly, the set
of sub-samples $\widetilde{\Gamma}_{m}$ consists of smaller sets
$\widetilde{\Gamma}_{m^{(l)}}^{(l)}$ that correspond to individual levels of
resolution. The sets $\widetilde{\Gamma}_{m^{(l)}}^{(l)}$ are filled as the
data streams in, and once a set for a given level is deemed complete, we
proceed with updating (2.4) and (2.5).

Further details of how the sets $\widetilde{\Gamma}_{m^{(l)}}^{(l)}$ are
constructed, and the corresponding criteria, are provided in Sections 3 and 4.
We begin by describing the multi-resolution framework of estimators.

###  2.2 The Laplacian pyramid

The LP [42, 29] is a multi-resolution regression method for extending a model
$\widehat{f}$ to out-of-sample data points ${\bm{\mathsf{x}}}\in{\cal
X}/\Gamma_{n}$. The LP can be formulated for radial kernels in the form

$k({\bm{\mathsf{x}}}_{i},{\bm{\mathsf{x}}}_{j})=\Phi\bigg{(}\frac{\|{\bm{\mathsf{x}}}_{i}-{\bm{\mathsf{x}}}_{j}\|}{r}\bigg{)},$
(2.6)

where $r>0$ is a shape parameter that determines the decay of $\Phi$ with
respect to $\|{\bm{\mathsf{x}}}_{i}-{\bm{\mathsf{x}}}_{j}\|$, see [17]. The
idea underpinning the LP is to approximate the target function sequentially,
where at each stage we regress on the errors from the previous stage. In other
words, we begin with a rough approximation using a large shape parameter for
which $\Phi$ decays slowly and then improve the approximation by fitting the
resulting error and reducing the shape parameter. In the LP, the estimator at
level $L\in\mathbb{N}$ is defined recursively as

$\widehat{f}^{(L)}({\bm{\mathsf{x}}})=\sum_{l=0}^{L}s^{(l)}({\bm{\mathsf{x}}})=s^{(L)}({\bm{\mathsf{x}}})+\widehat{f}^{(L)}({\bm{\mathsf{x}}}),$
(2.7)

where $\hat{f}^{(0)}=s^{(0)}$, and $s^{(l)}({\bm{\mathsf{x}}})$ is a
correction term defined by

$s^{(l)}({\bm{\mathsf{x}}})=\sum_{i=1}^{n}\alpha^{(l)}_{i}k^{(l)}({\bm{\mathsf{x}}},{\bm{\mathsf{x}}}_{i}).$
(2.8)

The coefficients
${\bm{\upalpha}}^{(l)}=(\alpha^{(l)}_{1},\ldots,\alpha^{(l)}_{n})^{\top}$ are
computed by conducting KRR on the residuals, i.e. errors, from the estimator
at the previous level. Namely,
${\bm{\upalpha}}^{(l)}=({\mathbf{K}}^{(l)}_{nn}+\lambda
n{\mathbf{I}})^{-1}{\bm{\mathsf{d}}}^{(l)}$, where

${\bm{\mathsf{d}}}^{(l)}=\begin{cases}{\bm{\mathsf{y}}},&\text{if}\quad l=0\\\
{\bm{\mathsf{y}}}-\widehat{f}^{(l-1)}([{\bm{\mathsf{x}}}_{n}]),&\text{otherwise}\end{cases}.$
(2.9)

For a FALKON adaption of this scheme, we only need to modify how per-level
coefficients are computed. Following (2.3) we iteratively solve

$({\mathbf{B}}^{(l)})^{\top}{\mathbf{H}}^{(l)}{\mathbf{B}}^{(l)}{\bm{\upbeta}}^{(l)}=({\mathbf{B}}^{(l)})^{\top}\big{(}{\mathbf{K}}^{(l)}_{nm}\big{)}^{\top}{\bm{\mathsf{d}}}^{(l)},$
(2.10)

where ${\mathbf{B}}^{(l)}$ is the corresponding preconditioner, and
${\mathbf{H}}^{(l)}=({\mathbf{K}}_{nm}^{(l)})^{\top}{\mathbf{K}}_{nm}^{(l)}+\lambda
n{\mathbf{K}}_{mm}^{(l)}$, and set
$\widetilde{{\bm{\upalpha}}}^{(l)}={\mathbf{B}}^{(l)}{\bm{\upbeta}}^{(l)}$.

######  Remark 2.1.

In this paper, we construct the kernel matrices ${\mathbf{K}}^{(l)}$ on a
particular class of radial kernels, namely the Gaussian kernel

$k^{(l)}({\bm{\mathsf{x}}},\widetilde{{\bm{\mathsf{x}}}}_{i})=\exp{\bigg{(}-\frac{\|{\bm{\mathsf{x}}}-\widetilde{{\bm{\mathsf{x}}}}_{i}\|^{2}}{2r_{l}^{2}}\bigg{)}},$

where $r_{l}>0$ is the shape parameter (the kernel bandwidth) at level $l$.

##  3 The damped cover tree

This work introduces a data-driven sub-sampling method that we call the damped
cover-tree (DCT). The DCT is a modification of the cover-tree (CT) [39], a
data structure based on partitioning a metric space, initially designed to
facilitate nearest neighbor search. The goal of the DCT is to modify and
simplify the CT to allow a viable sub-sampling scheme.

Let $({\cal X},\|\cdot\|)$ be a normed space where the input domain ${\cal
X}\subset\mathbb{R}^{D}$ is bounded, such that the diameter
$r_{0}=\operatorname{diam}({\cal X})$ is finite. The DCT is a tree structure
where each node $p$ of the tree is associated with a point
${\bm{\mathsf{x}}}_{p}\in{\cal X}$, and which is built sequentially as data
points arrive. Furthermore, let $Q_{l}$ be a set (herein called a cover-set)
containing all the nodes at a level $l\geq 0$ in the given tree. A level is
associated with an integer $l$ and a radius $r_{l}=2^{-l}r_{0}$, where $l=0$
denotes the root level containing only one node and $l$ increases as we
descend deeper into the tree. DCT has three invariants, of which the first two
are also invariants of the CT.

  1. (I1)

(Covering invariant) For all $p\in Q_{l+1}$ there exists $q\in Q_{l}$ such
that $\|{\bm{\mathsf{x}}}_{q}-{\bm{\mathsf{x}}}_{p}\|<r_{l}$.

  2. (I2)

(Separation invariant) For all $q,p\in Q_{l}$ where
${\bm{\mathsf{x}}}_{q}\neq{\bm{\mathsf{x}}}_{p}$, we have
$\|{\bm{\mathsf{x}}}_{q}-{\bm{\mathsf{x}}}_{p}\|>r_{l}$.

We add that the standard CT includes a third invariant, the so-called nesting
invariant, which requires $Q_{l}\subseteq Q_{l+1}$, but this is not desired
for our purpose.

To introduce the last invariant of the DCT, we first need the following
definition.

######  Definition 3.1 (The covering fraction).

Let $p\in Q_{l}$ be a node, and ${\bm{\mathsf{x}}}_{p}$ the associated point
in ${\cal X}$. Furthermore, let $\widetilde{C}_{p}=\\{c_{i}\\}_{i=1}^{k}$ be
the children of $p$, and ${\bm{\mathsf{x}}}_{c_{i}}$ the corresponding points
in ${\cal X}$. The covering fraction of a node $p$ is defined as

$\mathfrak{c}\mathfrak{f}(p)=\frac{\operatorname*{Vol}\bigg{(}{\cal
B}({\bm{\mathsf{x}}}_{p},r_{l})\cap\bigcup\limits_{c_{i}\in\widetilde{C}_{p}}{\cal
B}({\bm{\mathsf{x}}}_{c_{i}},r_{l+1})\bigg{)}}{\operatorname*{Vol}{\big{(}{\cal
B}({\bm{\mathsf{x}}}_{p},r_{l})\big{)}}}.$

The covering fraction is the proportion of the volume of ${\cal
B}({\bm{\mathsf{x}}}_{p},r_{l})$ that is covered by balls around its children
of half the radius. This quantity is directly related to (I2), which enforces
the radius $r_{l}$ to reduce by a factor of $2$ for each new level, starting
from an initial radius $r_{0}>0$. The covering fraction allows us to capture
the vulnerability of the standard CT to the curse of dimensionality.

For example, consider two regions ${\cal A}_{1},{\cal A}_{2}\subseteq{\cal
X}$, for which the doubling dimension at radius $r_{l}$ is
$\texttt{ddim}({\cal A}_{1},r_{l})>\texttt{ddim}({\cal A}_{2},r_{l})$. A node
$p\in{\cal A}_{1}$ at level $l$ will then need exponentially more children to
be covered, than a node $q\in{\cal A}_{2}$ at the same level $l$. This
exacerbates the deeper we go into the tree. Therefore, the CT would have
significantly more nodes from regions where the doubling dimension is large.

We recall now that sub-sampling is in kernel methods intended to reduce the
computational complexity. For this purpose, it is desirable to keep the number
of sub-samples from each level within a budget of reasonable size. On the
other hand, a too low sub-sample density will lead to poor interpolation
performance. Due to the exponential growth of the number of nodes with respect
to the doubling dimension, it would be desirable to avoid wasting our budget
on sub-samples from regions and radii with a large doubling dimension, as this
would require dedicating an (exponentially) large number of points to achieve
good interpolation, which is not feasible. Moreover, in high dimensional
regions, we likely cannot learn anything more than a simple function, for
which a lower sampling density would suffice.

To reduce the number of sub-samples from regions of large doubling
dimensionality, we introduce the following damping invariant as the third
invariant of the DCT.

  1. (I3)

(Damping invariant) Let ${\cal D}_{\mathfrak{c}\mathfrak{f}}\in(0,1)$ be some
threshold and let $\widetilde{C}_{p}$ and $\mathfrak{c}\mathfrak{f}(p)$ be as
in Def. 3.1. Then any node $q$ whose parent node $p$ does not satisfy
$\mathfrak{c}\mathfrak{f}(p)\geq{\cal D}_{\mathfrak{c}\mathfrak{f}}$ does not
have children of its own.

The damping invariant forces the tree to devote more resources to regions of
lower doubling dimension by making it harder for nodes in regions with higher
doubling dimensions to have children. In other words, the practical effect of
the damping invariant is to stop the vertical growth of the DCT if the
doubling dimension becomes large. This is because the covering number grows
exponentially with the dimensionality, ensuring
$\mathfrak{c}\mathfrak{f}(p)\geq{\cal D}_{\mathfrak{c}\mathfrak{f}}$ gets
correspondingly harder to achieve.

######  Remark 3.2.

In Section 5.1, we analyze the damping invariant in more detail and show how
the damping suppresses vertical growth of the DCT more for regions of high
doubling dimension than for regions of lower doubling dimensionality.

###  3.1 Construction of the DCT

We now discuss how the DCT is constructed and updated as the data streams in.
First, it is important to restate that we use the DCT to replace the Nyström
sampling, which was in FALKON used to reduce the complexity of the ridge
regressor. Consequently, not all of the streamed data (that is, not every
training point) will be added to the tree, but only those whose inclusion into
the tree would not violate the invariants (I1)-(I3). In other words, the tree
consists of only those training points that help resolve the data space at the
relevant resolution level. Thus, each node $p$ in the DCT is associated with a
unique training sample ${\bm{\mathsf{x}}}_{p}$, but not every training sample
will be represented by a node in the tree. Note that this is different from
the standard CT, which aims to organize all of the training data into a
geometrical leveled data structure.

The construction of the DCT consists of a series of checks which examine
whether adding a given data point to the DCT would, or would not, violate
invariants (I1)\- (I3). When a new point ${\bm{\mathsf{x}}}_{q}$ arrives from
the data stream the goal is to identify the deepest level $l$ for which there
exists a node $p$ such that
$\|{\bm{\mathsf{x}}}_{q}-{\bm{\mathsf{x}}}_{p}\|\leq r_{l}$. This corresponds
to finding the nearest node in the tree that could serve as a parent.

We achieve this in the following way. The first training point is identified
as the root node to which we associate the radius $r_{0}$. For each new point,
we proceed in a top-down manner, starting from the root node888We assume that
all new points ${\bm{\mathsf{x}}}_{q}$ are within a ball of radius $r_{0}$
around this node, which holds for a large enough $r_{0}$. We then check
whether ${\bm{\mathsf{x}}}_{q}$ would violate the separation invariant at the
next level. In other words, if there exists a node $p$ such that
$\|{\bm{\mathsf{x}}}_{q}-{\bm{\mathsf{x}}}_{p}\|<r_{l}$. If such a node does
not exist, then ${\bm{\mathsf{x}}}_{q}$ is added to the set of children of the
root node, and we update the covering fraction estimate for the root node.
Otherwise, if such a node does exist, we repeat the process, checking the
separation invariant among the children of the corresponding node, and proceed
further down the tree.

Assume we arrived to a node $p$ at level $l\geq 1$, and we have
$\|{\bm{\mathsf{x}}}_{q}-{\bm{\mathsf{x}}}_{p}\|\leq r_{l}$. We then check if
$p$ is allowed to have children, that is if the damping invariant is
satisfied. If it is not satisfied, the point ${\bm{\mathsf{x}}}_{q}$ is
dismissed (it is not added to the tree). On the other hand, if $p$ is allowed
to have children, we check whether the separation invariant holds, i.e., if
there exists a child $c$ of the node $p$ such that
$\|{\bm{\mathsf{x}}}_{q}-{\bm{\mathsf{x}}}_{c}\|<r_{l+1}$. If that were the
case, the separation invariant would be violated, and the recursion is applied
again by considering $c$ as the potential parent node. However, if such a
child does not exist, that is, if the separation invariant is not violated,
then ${\bm{\mathsf{x}}}_{q}$ is added to the set of children of the node $p$.
More details are given in Alg. A.1.

Some comments are needed to elucidate how are the steps described above
applied in practice. First, note that the covering fraction from Def. 3.1
cannot be calculated explicitly, since the volume terms require knowing the
intrinsic dimensionality. Therefore, it is necessary to use an estimator
instead. For this purpose, we interpret $\mathfrak{c}\mathfrak{f}(p)$ as the
probability that a sample ${\bm{\mathsf{x}}}\sim\textrm{Uni}({\cal
B}({\bm{\mathsf{x}}}_{p},r))$ will be within ${\cal
B}_{c}\coloneqq\bigcup_{c_{i}\in\widetilde{C}_{p}}{\cal
B}({\bm{\mathsf{x}}}_{c_{i}},r/2)$, where $\widetilde{C}_{p}$ are the children
of $p$. This probability can be estimated by considering the checks of the
separation invariant (I2), conducted on the last $N$ points that were inside
${\cal B}({\bm{\mathsf{x}}}_{p},r)$, as a series of independent random trials.
We use the following running average as an estimator of the covering fraction

$(\mathfrak{c}\mathfrak{f}(p))_{t}=(1-\alpha)(\mathfrak{c}\mathfrak{f}(p))_{t-1}+\alpha\mathbbm{1}_{{\cal
B}_{c}}({\bm{\mathsf{x}}}_{t}),$ (3.1)

where $\mathbbm{1}_{{\cal B}_{c}}({\bm{\mathsf{x}}}_{t})$ is the indicator
function, and $\alpha>0$ is a weighting parameter. This approximates a
weighted average of the outcome of the $N$ last draws (cf. Appendix Appendix
B. Preparatory material). Note that this reduces the memory requirements,
since instead of storing $N$ trial outcomes for each node in the tree, as
required had we used an average of the last $N$ trials, we store only a single
value for each node in the tree.

Second, the separation invariant is in practice too strict since it results in
too few points added to the tree, and thus a worse kernel estimator. Moreover,
checking the separation invariant adds to the computational complexity.
Therefore, we introduce the following relaxation. Assume we have a new point
${\bm{\mathsf{x}}}_{q}$ and arrived at a node $p$ at level $l$. We then first
conduct a random Bernoulli trial, with the failure probability

$q_{\bm{\mathsf{x}}}=\frac{1}{1+\exp{\big{[}\\-h\tan\big{(}\pi(\|{\bm{\mathsf{x}}}_{q}-{\bm{\mathsf{x}}}_{p}\|/r_{l}-\frac{1}{2})\big{)}\big{]}}},$
(3.2)

where $h$ is the hardness of the threshold. In other words, the probability of
failure is proportional to the distance between ${\bm{\mathsf{x}}}_{q}$ and
${\bm{\mathsf{x}}}_{p}$ \- the larger the distance, the more likely the
failure. If the trial’s outcome is a failure, then the check for the
separation invariant is ignored, and the algorithm continues. If it is a
success, we proceed by first checking the separation invariant. This means
that the probability to ignore the separation invariant increases as
${\bm{\mathsf{x}}}_{q}$ gets farther from ${\bm{\mathsf{x}}}_{p}$.

###  3.2 Sub-sampling from the DCT

We now discuss how the DCT is used for sub-sampling the training points. By
organizing the training points into cover-sets $Q_{l}$ the DCT allows a
hierarchical sub-sampling. Even though cover-sets $Q_{l}$ significantly reduce
the number of training points, they are for practical purposes still too large
for efficient sub-sampling. Due to this, we restrict ourselves to a subset
$\widetilde{\Gamma}^{(l)}\subseteq Q_{l}$ of candidate sub-samples called
landmarks.

######  Definition 3.3 (Landmarks).

Let $Q_{l}$ be the cover-set at level $l$ in a DCT. We define the set of
candidate landmarks at level $l$ as
$\widetilde{\Gamma}^{(l)}=\\{{\bm{\mathsf{x}}}_{p}\mid p\in Q_{l}\text{ and
}\mathfrak{c}\mathfrak{f}(p)\geq\mathcal{D}_{\mathfrak{c}\mathfrak{f}}\\}$,
and the set of landmarks (of size $m$) as any subset
$\widetilde{\Gamma}^{(l)}_{m}=\\{\widetilde{{\bm{\mathsf{x}}}}^{(l)}_{1},\ldots,\widetilde{{\bm{\mathsf{x}}}}^{(l)}_{m}\\}\subset\widetilde{\Gamma}^{(l)}$
of size $m$.

Some remarks are in order. First, by Def. 3.3, candidates for landmarks at
level $l$ are only those nodes allowed to have children (according to the
damping invariant (I3)). This design choice implies that the set of candidate
landmarks will contain more points from regions with a lower doubling
dimension than points from regions with a higher doubling dimension. This is
because the larger the doubling dimension is, the more children nodes are
needed to cover a given parent node.

Second, Def. 3.3 suggests using only a subset of candidate landmarks as sub-
samples. We refer to a result from [2] which states that good statistical
accuracy of the estimator is achieved if the number of sub-samples is
proportional to the square root of the number of samples. At level $l$ we
therefore use a set of landmarks which is of size
$m^{(l)}=\delta_{0}\sqrt{\lvert{Q_{l}}\rvert}$, where $\delta_{0}>0$ is a
constant.

The third point that requires attention concerns the question of when the
landmarks should be selected. To that end, we use the covering fraction of a
level, which, with a slight abuse of notation, we denote as
$\mathfrak{c}\mathfrak{f}(Q_{l})$. Moreover, we compute
$\mathfrak{c}\mathfrak{f}(Q_{l})$ as

$(\mathfrak{c}\mathfrak{f}(Q_{l}))_{t}=(1-\alpha)(\mathfrak{c}\mathfrak{f}(Q_{l}))_{t-1}+\alpha\mathbbm{1}_{{\cal
B}_{\text{level}}}({\bm{\mathsf{x}}}_{t}),$ (3.3)

where $\mathcal{B}_{level}=\bigcup\limits_{p\in
Q_{l}}\mathcal{B}({\bm{\mathsf{x}}}_{p}^{(l)},r_{l})$. Moreover, analogously
to the damping invariant, let $\mathcal{D}_{level}\in(0,1)$ be some threshold.
We then say that a level $l$ is sufficiently covered when
$\mathfrak{c}\mathfrak{f}(Q_{l})\geq\mathcal{D}_{level}$.

######  Remark 3.4.

We note that as the level increases, our estimate of
$\mathfrak{c}\mathfrak{f}(Q_{l})$ through Eq. (3.3) will be increasingly more
sensitive to subsets ${\cal A}\subset{\cal X}$ of low doubling dimension than
to subsets of large doubling dimension. This is because the damping invariant
(I3) makes it harder for nodes in high dimensions to have children.
Consequently, we will have fewer points in deeper levels that belong to high
dimensional regions. Because of this, the estimator in Eq. (3.3) is biased
towards using more sub-samples from lower dimensional regions.

Sub-sampling from a level $l$ goes as follows. As training points arrive, we
build the tree and continuously update the covering fraction of a level. Once
that level is sufficiently covered, that is, once
$\mathfrak{c}\mathfrak{f}(Q_{l})\geq\mathcal{D}_{level}$, we extract the set
of landmarks by sub-sampling $m^{(l)}$ points from the pool of candidate
landmarks $\widetilde{\Gamma}^{(l)}$.

##  4 StreaMRAK

In this section, we present StreaMRAK and clarify how it synthesizes concepts
from Sections 2 and 3, and utilizes them in a streaming context. The workflow
of StreaMRAK can be divided into three threads that can run in parallel,
subject to some inter-dependencies. These are the sub-sampling thread, the
training thread, and the prediction thread. Overviews of these threads are
given next, and the reader is referred to Algorithm A.2 in the Appendix for
further details.

###  4.1 Sub-sampling thread

In the sub-sampling thread StreaMRAK collects and organizes the training data
into a DCT. Namely, as new training pairs are collected, the covering (I1) and
separation (I2) are checked, and the covering fraction is updated as described
in Section 3.1. Moreover, the set of landmarks for each level is updated, as
described in Section 3.2. Once the set of landmarks for a given level
$\widetilde{\Gamma}^{(l)}_{m}$ is completed, the landmarks and the estimator
for the corresponding level can be used in the remaining two threads.

###  4.2 Training thread

The model is trained at level $l$ when two conditions are met. First,
coefficients of the previous level $l-1$ in the LP must have been calculated,
i.e. previous training thread must finish. Second, landmarks
$\widetilde{\Gamma}^{(l)}_{m^{(l)}}$ at level $l$ must be ready.

In the first step, we define the kernel matrix on the landmarks by

$[{\mathbf{K}}_{mm}^{(l)}]_{ij}=k^{(l)}(\widetilde{{\bm{\mathsf{x}}}}_{i},\widetilde{{\bm{\mathsf{x}}}}_{j}),\text{
for }\widetilde{{\bm{\mathsf{x}}}}_{i}\in\widetilde{\Gamma}^{(l)}.$ (4.1)

In the second step we consider
$\big{(}{\mathbf{K}}^{(l)}_{nm}\big{)}^{\top}{\mathbf{K}}^{(l)}_{nm}\in\mathbb{R}^{m^{(l)}\times
m^{(l)}}$ and
$\big{(}{\mathbf{K}}^{(l)}_{nm}\big{)}^{\top}{\bm{\mathsf{d}}}^{(l)}_{n}\in\mathbb{R}^{m^{(l)}}$
which in addition to landmarks depend on the training points. They are updated
continuously as new training points come in, according to Eq. (2.4) and Eq.
(2.5). However, they are not updated indefinitely, but only until new training
points do not significantly alter the matrices according to the following
criterion.

######  Definition 4.1.

(Sufficient training points) Let
${\mathbf{A}}_{n}:=({\mathbf{K}}^{(l)}_{nm})^{\top}{\mathbf{K}}^{(l)}_{nm}$,
and
${\bm{\mathsf{b}}}_{n}:=\big{(}{\mathbf{K}}^{(l)}_{nm}\big{)}^{\top}{\bm{\mathsf{d}}}^{(l)}_{n}$.
Let $\delta_{1},\delta_{2},\delta_{3}>0$ be three constants. We consider the
number of training points at a level $l$ sufficient when either
$n\geq\delta_{3}$ or

$\bigg{\|}\frac{{\mathbf{A}}_{n}}{n}-\frac{{\mathbf{A}}_{n+1}}{n+1}\bigg{\|}_{\infty}\leq\delta_{1}\text{
and
}\quad\bigg{\|}\frac{{\bm{\mathsf{b}}}_{n}}{n}-\frac{{\bm{\mathsf{b}}}_{n+1}}{n+1}\bigg{\|}\leq\delta_{2}.$

After enough training samples are collected according to Def. 4.1, the
correction term $s^{(l)}$ is obtained by solving for the coefficients
$\widetilde{\alpha}^{(l)}_{1},\dots,\widetilde{\alpha}^{(l)}_{m^{(l)}}$ using
Eq. (2.10). The new prediction model $\widehat{f}^{(L)}$ is obtained by adding
$s^{(l)}$ to the previous model, according to Eq. (2.7).

###  4.3 Prediction thread

In this thread StreaMRAK makes provides the latest version of the trained LP
model in Eq. (2.7). This means that if $L$ is currently the highest level that
has been trained, the prediction for new points ${\bm{\mathsf{x}}}$ is made
using the model $\widehat{f}^{(L)}({\bm{\mathsf{x}}})$.

##  5 Analysis

In this section, we first analyze the damping invariant of the DCT. We then
offer theoretical results on the convergence properties of the LP in the
context of KRR. Finally, we offer estimates of the time and memory
requirements of StreaMRAK .

###  5.1 Analysis of the DCT

As discussed in Section 3, the DCT adds a given training point to the set of
nodes of the tree if conditions (I2) and (I3) are satisfied, and the points
are otherwise discarded. In particular, the damping invariant (I3) makes it
harder for a node to have children. The guiding idea is that damping should
reduce the impact of the curse of dimensionality by making it harder for nodes
in regions of higher doubling dimension to have children, and in doing so it
should effectively stop the vertical growth of the tree in corresponding
regions. Therefore, it is critical to understand how and to what degree the
damping affects high dimensional regions more than low dimensional ones.

In a statistical sense, the damping should treat all nodes in regions of the
same doubling dimension equally. Therefore, to gain insight into the damping,
it suffices to analyze its effects concerning the doubling dimension on a
single node $p$. In this case, the effect of damping can be measured by
analyzing how many training points must pass through $p$, in the sense of Alg.
A.1, before children of $p$ are allowed to have children of their own. This
can be modeled by considering the expected number of training points
${\bm{\mathsf{x}}}_{i}\sim{\rm Uni}\big{(}{\cal
B}({\bm{\mathsf{x}}}_{p},r)\big{)}$ necessary to cover ${\cal
B}({\bm{\mathsf{x}}}_{p},r)$ with balls of radius $r/2$ around points
${\bm{\mathsf{x}}}_{i}$.

Consider ${\bm{\mathsf{x}}}_{i}\sim{\rm Uni}\big{(}{\cal
B}({\bm{\mathsf{x}}}_{p},r)\big{)}$, and let a set ${\cal S}_{p}$ be built in
a succession of trials $i=1,\ldots,N_{t}$ so that

${\bm{\mathsf{x}}}_{i}\in{\cal S}_{p}\text{ if
}\|{\bm{\mathsf{x}}}_{i}-{\bm{\mathsf{x}}}\|\geq\frac{r}{2}\text{ for all
}{\bm{\mathsf{x}}}\in{\cal S}_{p}.$

In other words, a newly sampled point ${\bm{\mathsf{x}}}_{i}$ will only be
added to the set ${\cal S}_{p}$ if it its pairwise distances from all the
points that are already in ${\cal S}_{p}$ are at least $r/2$.

######  Problem 1.

Let $\widetilde{C}_{p}$ denote the set of children of the node $p$,
constructed from the above-described trials. What is the expected number of
trials $N_{t}$ needed to ensure $\mathfrak{c}\mathfrak{f}(p)=1$?

Since there is no unique set ${\cal S}_{p}$ such that the corresponding set of
children $\widetilde{C}_{p}$ ensures $\mathfrak{c}\mathfrak{f}(p)=1$, the
sample space for Problem 1 corresponds to all admissible sets ${\cal S}_{p}$,
which vary in both the number and the location of points they contain.
Characterizing all such sets corresponds to a disordered sphere packing
problem [56], which is an NP-hard combinatorial problem [57]. For a
theoretical analysis of this problem, defining a probability measure over the
sample space is necessary. However, in this level of generality, neither the
sample space nor the probability measure admit a workable definition, with
currently available mathematical tools [56]. Although some theoretical
insights are possible under simplifications on the sample space, this analysis
is restrained to a limited number of spheres and configurations.

Due to these difficulties, we consider a simplified setting where we instead
consider an average case. If the set ${\cal S}_{p}$ is such that ${\cal
B}({\bm{\mathsf{x}}}_{i},r)\subset\bigcup_{{\bm{\mathsf{x}}}_{i}\in{\cal
S}_{p}}{\cal B}({\bm{\mathsf{x}}}_{i},r/2)$, which corresponds to
$\mathfrak{c}\mathfrak{f}(p)=1$, then each of the balls ${\cal
B}({\bm{\mathsf{x}}}_{i},r/2)$ occupies on average $\frac{1}{|{\cal S}_{p}|}$
of the total volume of ${\cal B}({\bm{\mathsf{x}}}_{p},r)$, assuming none of
the balls are covered by a union of other balls. Therefore, as ${\cal S}_{p}$
is being built, adding a point to ${\cal S}_{p}$ will, on average, reduce the
unoccupied volume of ${\cal B}({\bm{\mathsf{x}}}_{p},r)$ by $\frac{1}{|{\cal
S}_{p}|}$. Moreover, it can be shown that the number of elements in such a set
satisfies $2^{\texttt{ddim}-1}\leq|{\cal S}_{p}|\leq 5^{\texttt{ddim}}$, see
Lemma C.1, where $\texttt{ddim}\coloneqq\texttt{ddim}({\cal S}_{p},r)$ is the
doubling dimension of ${\cal S}_{p}$. Based on these considerations we
introduce a simplified setting for the average case of Problem 1.

######  Assumption 1.

Problem 1 can be approximated by dividing the ball ${\cal
B}({\bm{\mathsf{x}}}_{p},r)$ into a union of $c_{d}$ fixed (and known)
disjoint bins ${\cal B}_{i}$ of size $(1/c_{d})\operatorname*{Vol}\big{(}{\cal
B}({\bm{\mathsf{x}}}_{p},r)\big{)}$.

Note that the bins referred to in Assumption 1 correspond to regions around
the children of the node $p$. Assumption 1 reduces the average case of Problem
1 to a form of the classical coupons collector’s problem [58], which considers
$n$ coupons with the same probability of being drawn. Through a series of
randomized trials with replacement, the goal is to obtain a copy of each
coupon. Relevant for Problem 1 is estimating the stopping time $T$, which
counts the number of trials before all coupons are collected, and which
satisfies $\mathbb{E}[T]=nH_{n}$, where $n$ denotes the number of coupons and
$H_{n}$ is the $n$-th harmonic number [58].

In terms of Problem 1, and under Assumption 1, we can therefore identify
$T=N_{t}$, $n=|{\cal S}_{p}|$ and $\mathbb{E}[N_{t}|\text{Node}\,p]=|{\cal
S}_{p}|H_{|{\cal S}_{p}|}$. Combining the bound $\ln(n)+\frac{1}{2}\leq
H_{n}\leq\ln(n)+1$ (from [59]), with the bound on $|{\cal S}_{p}|$ from Lemma
C.1 we have

$2^{\texttt{ddim}-1}((\texttt{ddim}-1)\ln
2+1/2)\leq\mathbb{E}[N_{t}|\text{Node}\,p]\leq
5^{\texttt{ddim}}(\texttt{ddim}\ln 5+1).$ (5.1)

With the same strategy, we can bound the number of trials until the cover-
fraction of a level reaches $1$, as

$2^{l(\texttt{ddim}-1)}({l(\texttt{ddim}-1)}\ln
2+1/2)\leq\mathbb{E}[N_{t}|\text{Level}\,l]\leq
5^{l\texttt{ddim}}(l\texttt{ddim}\ln 5+1).$ (5.2)

From Eq. (5.1) we see that the number of training points
$\mathbb{E}[N_{t}|\text{node}\,p]$ grows exponentially with the doubling
dimensionallity $d$. In other words, significantly more trials are needed to
achieve $\mathfrak{c}\mathfrak{f}(p)={\cal D}_{\mathfrak{c}\mathfrak{f}}$ for
nodes in regions with a large doubling dimension than it is for nodes in
regions with a lower doubling dimension. Consequently, through the damping
invariant, the DCT restricts the vertical growth of the tree comparatively
more the higher the doubling dimension of the local region.

###  5.2 Time and memory requirements

This section analyzes the memory requirements of StreaMRAK , which involve
storing the DCT and the linear system components used to update the
coefficients. Furthermore, we consider the computational requirements, which
consist in solving the coefficient equations. Both the memory and
computational requirements need to be analyzed per level $l$ of the tree due
to the multi-resolution nature of the estimator and the tree organization of
the data.

For the analysis, we consider a simplified setting where we assume that the
doubling dimension is constant for all levels and all subsets of ${\cal X}$,
and that the number of children $c_{d}$ is the same for all nodes. At the end
of the section we describe a more general setting.

In the following, we assume that the growth of the DCT stops at a level $L$.
In other words, level $L$ is the last level at which there are nodes. In
practice, the growth of the DCT slows down exponentially fast with the product
of the doubling dimension $\texttt{ddim}\coloneqq\texttt{ddim}({\cal
X},r_{L})$ and the level $l$. This can be seen from Eq. (5.2), which shows
that the number of training points necessary to fill up a level grows
exponentially with $l\texttt{ddim}$. Therefore, in practice, no new levels
will be added to the DCT when $l\texttt{ddim}$ is large enough, which
effectively makes the last level $L$ independent of the number of training
points. Furthermore, from Lemma C.1 we know that $c_{d}$ is bounded by
$2^{\texttt{ddim}-1}\leq c_{d}\leq 5^{\texttt{ddim}}$, which shows that also
$c_{d}$ is independent of the number of training points.

######  Proposition 5.1.

The memory requirement of StreaMRAK is
$\mathcal{O}\big{(}\sum_{l=0}^{L}c_{d}^{l}\big{)}$.

###### Proof.

The memory requirement of the DCT is determined by the number of nodes in the
tree. Given that the number of children is the same for all nodes. If the
number of children per node is $c_{d}$, then the total number of nodes at
level $l$ is $c_{d}^{l}$. Thus, the memory needed to store the DCT with $L$
levels is ${\cal O}(\sum_{l=0}^{L}c_{d}^{l})$.

To store the linear system on level $l$ we need the matrices
$\big{(}{\mathbf{K}}_{nm^{(l)}}\big{)}^{\top}{\mathbf{K}}_{nm^{(l)}},\,{\mathbf{K}}_{m^{(l)}m^{(l)}}\in\mathbb{R}^{m^{(l)}\times
m^{(l)}}$ and the vector ${\bm{\mathsf{z}}}\in\mathbb{R}^{m^{(l)}}$. The
number of landmarks $m^{(l)}$ at level $l$ is chosen as
$m^{(l)}=\delta_{0}\sqrt{|Q_{l}|}$, where $|Q_{l}|$ is the number of nodes at
level $l$. Since $|Q_{l}|$ is ${\cal O}(c_{d}^{l})$, it follows that
$m^{(l)}\times m^{(l)}$ is also ${\cal O}(c_{d}^{l})$ per level, and the
desired conclusion follows. ∎

Note that with a fixed $L$ and $n$ larger than
$\mathcal{O}\big{(}\sum_{l=0}^{L}c_{d}^{l}\big{)}$, then the memory
requirement is independent of $n$. We also note that if the deepest level
satisfies $L\rightarrow\infty$, then the number of nodes is determined by the
number of training points, and the memory requirement would thus, in the worst
case, become $\mathcal{O}(n)$, the same as for the standard cover-tree.

Next, we discuss the construction of the DCT, where adding a new point to the
set of nodes requires a search through the tree.

######  Proposition 5.2.

Inserting a new point into the DCT, cf. Algorithm A.1, requires ${\cal
O}(c_{d}L)$ operations.

###### Proof.

For a point ${\bm{\mathsf{x}}}_{q}\in{\cal X}$ to be analyzed at level $L$, we
need to have analyzed it at the previous $l<L$ levels. At each level, we must,
in the worst case, check the separation invariant with all children of the
current potential parent $p^{(l)}$, before finding a node $c$ such that
$\|{\bm{\mathsf{x}}}_{q}-{\bm{\mathsf{x}}}_{c}\|\leq 2^{-l}r_{0}$, that would
serve as the next potential parent. This requires at most $c_{d}$ operations
per level, giving $Lc_{d}$ total operations over the $L$ levels. The same
number of operations is necessary if a node is discarded at level $L$. ∎

Lastly, we analyze the computational requirements for solving the linear
system.

######  Proposition 5.3.

The time requirement for solving the linear system in Eq. (2.3) is
$\mathcal{O}\big{(}\delta_{3}m^{(l)}+\big{(}m^{(l)}\big{)}^{3}\big{)}$ per
level, where $\delta_{3}$ is given in Def. 4.1,

###### Proof.

The time requirement of FALKON is ${\cal O}(nmt+m^{3})$ where $n$ is the
number of training points, $m$ the number of landmarks and $t$ the number of
iterations of the conjugate gradient (which has an upper bound). By Def. 4.1,
StreaMRAK uses at most $\delta_{3}$ training samples at each level. Since
$m^{(l)}$ is the number of landmarks at level $l$, the result follows. ∎

Assume that the domain ${\cal X}$ can be divided into disjoint subsets ${\cal
A}_{1},\dots,{\cal A}_{t}\subset{\cal X}$ for which the doubling dimension
$\texttt{ddim}({\cal A}_{i},r_{l})$ differs based on ${\cal A}_{i}$ and radius
$r_{l}$. Let the number of children of a node ${\bm{\mathsf{x}}}_{p}\in{\cal
A}_{i}$ at level $l$ be $c_{d,i,l}$. In this scenario, the growth of the DCT
will stop at different levels $L_{i}$ for different subsets ${\cal A}_{i}$.
The final time and memory requirements would therefore be the sum of the
contribution from each subset ${\cal A}_{i}$. In other words, the memory would
be $\mathcal{O}(\sum_{i=1}^{t}\sum_{l=0}^{L_{i}}c_{d,i,l}^{l})$, and similarly
the time requirement per point insertion would be
$\mathcal{O}(\sum_{i=1}^{t}\sum_{l=0}^{L_{i}}c_{d,i,l})$. We note that
$c_{d,i,l}$ and $L_{i}$ depend on the dimensionality of the data, but are
independent of $n$. Therefore, so are the time and memory requirements.

###  5.3 Convergence of the LP formulation of the KRR

This section analyzes the conditions for which the LP approximates the
training data $y_{i}=f({\bm{\mathsf{x}}}_{i})$, with respect to the number of
levels. A similar analysis was previously done for the LP in the context of
kernel smoothers [30]. However, to the best of our knowledge, this is the
first time the LP formulation of KRR has been analyzed in this way.

Consider the LP estimator $\widehat{f}^{(l)}$ as defined in Eq. (2.7), but
without sub-sampling. From the recurrence relationship for the residuals
${\bm{\mathsf{d}}}^{(l)}$ in Eq. (2.9) by induction it follows

$\widehat{f}^{(l+1)}([{\bm{\mathsf{x}}}_{n}])-f([{\bm{\mathsf{x}}}_{n}])=({\mathbf{I}}-{\mathbf{P}}_{nn}^{(l)})(\widehat{f}^{(l)}([{\bm{\mathsf{x}}}_{n}])-f([{\bm{\mathsf{x}}}_{n}]),$
(5.3)

where
${\mathbf{P}}_{nn}^{(l)}\coloneqq{\mathbf{K}}^{(l)}_{nn}({\mathbf{K}}^{(l)}_{nn}+\lambda
n{\mathbf{I}})^{-1}$, cf. Lemma C.2.

######  Theorem 5.4.

Let $\widehat{f}^{(l)}$ be the LP estimator defined in Eq. (2.7) and let
$\lambda$ be a regularization parameter. Furthermore, let
$0<\sigma_{l,n}\leq\dots\leq\sigma_{l,1}$ be the eigenvalues of
${\mathbf{K}}^{(l)}_{nn}$. For $L>0$ we then have

$\|\widehat{f}^{(L+1)}([{\bm{\mathsf{x}}}_{n}])-f([{\bm{\mathsf{x}}}_{n}])\|\leq\prod_{l=0}^{L}(1-\varepsilon(l))\|\widehat{f}^{(0)}([{\bm{\mathsf{x}}}_{n}])-f([{\bm{\mathsf{x}}}_{n}])\|,\quad\text{where}\quad\varepsilon(l)=\frac{\sigma_{l,n}}{n\lambda+\sigma_{l,n}}.$

From Thm. 5.4 it follows that the LP estimator will converge as
$l\rightarrow\infty$, since $\sigma_{l,n}>0$ and therefore
$1-\varepsilon(l)\in(0,1)$ for all $l$. In Thm. 5.5 we characterise how
$\varepsilon(l)$ depends on the level $l$ to give insight on the nature of
this convergence.

######  Theorem 5.5.

The LP estimator $\widehat{f}^{(l)}$ from Eq. (2.7) converges with increasing
level $L$ to the training data $f({\bm{\mathsf{x}}}_{i})$, cf. Thm. 5.4, with
the rate $\prod_{l=0}^{L}(1-\varepsilon(l))$, where

$1-\varepsilon(l)\leq\big{(}1+C_{1,D}2^{-Dl}\exp\big{(}-C_{2,D}4^{-l}\big{)}/n\lambda\big{)}^{-1},$
(5.4)

for

$C_{1,D}=\frac{1}{2}(6\sqrt{2})^{D}\Gamma(D/2+1)^{\frac{D-1}{D+1}}\bigg{(}\frac{\pi}{9}\bigg{)}^{\frac{D}{D+1}}\bigg{(}\frac{r_{0}}{\delta}\bigg{)}^{D}\quad\text{and}\quad
C_{2,D}=1152\bigg{(}\frac{\pi\Gamma^{2}(D/2+1)}{9}\bigg{)}^{\frac{2}{D+1}}\bigg{(}\frac{r_{0}}{\delta}\bigg{)}^{2},$

where $\Gamma$ is the gamma function.

Furthermore, for $l>\log_{2}(\sqrt{D/2}(r_{0}/\delta))$ we have the tighter
bound

$1-\varepsilon(l)<\bigg{(}1+\big{(}1-2^{1+\frac{1}{\ln{2}}(C_{3}D-g(l))}\big{)}/n\lambda\bigg{)}^{-1},$
(5.5)

where $g(l)=4^{l-\log_{2}{r_{0}/\delta}}$ and $C_{3}=(\ln{(1+1/4)}+2\ln{2})$.

We note that the bound in Eq. (5.4) underestimates the rate of convergence for
lower levels but improves as the levels increase. Furthermore, Thm. 5.5 shows
that the convergence rate increases with the level $l$. In fact, the bound in
Eq. (5.4) can be simplified with an a fortiori bound of the same form, where
$C_{1,D}=\frac{1}{2}\big{(}\frac{12.76}{2^{3/2}}\big{)}^{D}\big{(}\frac{D^{D}}{\Gamma(D/2+1)}\big{)}\big{(}\frac{r_{0}}{\delta}\big{)}^{D}$
and $C_{2,D}=(12.76\sqrt{2}D)^{2}(r_{0}/\delta)^{2}$, which ensures that
$1-\varepsilon(l)$ decreases monotonically for
$l<\log_{2}(\sqrt{D/2}(r_{0}/\delta))+\log_{2}(25.52\sqrt{2})$. see Remark C.3
and Corollary C.4.

On the other hand, when $l>\log_{2}(\sqrt{D/2}(r_{0}/\delta))$ the tighter
bound from Eq. (5.5) ensures that $1-\varepsilon(l)$ continues to decreases
monotonically. Moreover, as $l\rightarrow\infty$ each new level reduces the
residual error by $(1+1/n\lambda)^{-1}$. We can also observe that the
convergence rate is reduced by the number of training points $n$, but this
effect can be mitigated by reducing the regularization parameter $\lambda$. We
also note that Thm. 5.4 and Thm. 5.5 are derived for a vector of numbers on
the training data $\Gamma_{n}\subset{\cal X}$, without assumptions on the
target function. In other words, the LP estimator can approximate the training
data for any function $f:\Gamma_{n}\rightarrow\mathbb{R}$, to arbitrary
precision, by including sufficiently many levels.

######  Corollary 5.6.

If the residual
${\bm{\mathsf{d}}}^{(l)}=(\widehat{f}^{(l)}([{\bm{\mathsf{x}}}_{n}])-f([{\bm{\mathsf{x}}}_{n}]))$
at level $l$ only projects non-trivially onto the eigenvectors with eigenvalue
$\sigma_{l,n}\geq\sigma_{\text{cut-off}}$, then we say the residual is
spectrally band-limited with respect to the kernel. If the residual
${\bm{\mathsf{d}}}^{(l)}$ is spectrally band-limited, then
$1-\epsilon(l)<n\lambda/(n\lambda+\sigma_{\text{cut-off}})$.

##  6 Experiments

This section presents comparative numerical experiments of the proposed
estimator on three problems. In Section 6.1 we consider a one-dimensional
regression problem, and in Section 6.2 we consider a dumbbell-shaped domain
that consists of two 5-dimensional spheres connected by a 2-dimensional plane.
Lastly, in Section 6.3, we forecast the trajectory of a double pendulum, which
is a well-known chaotic system [40].

We compare StreaMRAK with FALKON [2] and an LP modification of KRR (LP-KRR).
Both FALKON and LP-KRR rely on the standard Nyström sub-sampling [47, 48].
Furthermore, FALKON does not rely on a multi-resolution scheme but uses
instead a single bandwidth, found by cross-validation.

Throughout the experiments, we set the threshold for the number of sub-samples
(landmarks) in StreaMRAK to be $10\sqrt{\lvert{Q_{l}}\rvert}$, where $Q_{l}$
is the set of nodes at level $l$ in the DCT. We note that to choose the sub-
sample size, FALKON and LP-KRR require $n$ to be known beforehand. For FALKON
we let the number of Nyström landmarks be $10\sqrt{n}$, where $n$ is the
number of training samples. Meanwhile, for LP-KRR  we sub-sample $\sqrt{n}$
Nyström landmarks, which are then used for all levels.

We also need to pre-select the number of training points for LP-KRR and
FALKON. For FALKON  we use the entire training set, as in [2]. Similarly, it
is also common for the LP to use the entire training set at each level [29,
30]. However, for large data sets, it might be better to include fewer data
points. Therefore, we also use a version of the LP-KRR where we divide the
total training data equally between the levels.

Throughout the experiments, we measure the performance of StreaMRAK, FALKON,
and LP-KRR by estimating the mean square error

$MSE(y,y_{\textit{pred}})=\frac{1}{\Upsilon\Lambda}\sum_{k=1}^{\Upsilon}\frac{1}{n_{k}}||{\bm{\mathsf{y}}}_{k}-{\bm{\mathsf{y}}}_{k}^{\textit{pred}}||^{2},\text{
with }\Lambda=\max_{\begin{subarray}{c}k\in[\Upsilon]\\\
i\in[n_{k}]\end{subarray}}[{\bm{\mathsf{y}}}_{k}]_{i}-\min_{\begin{subarray}{c}k\in[\Upsilon]\\\
i\in[n_{k}]\end{subarray}}[{\bm{\mathsf{y}}}_{k}]_{i},$ (6.1)

where $\Upsilon$ is the number of test runs we average over, $n_{k}$ is the
number of test points at test run $k$, and
${\bm{\mathsf{y}}}_{k},{\bm{\mathsf{y}}}_{k}^{\textit{pred}}\in\mathbb{R}^{n_{t}}$
are the target values and predictions respectively, and $\Lambda$ is the
normalisation factor.

###  6.1 Multi-resolution benchmark

We consider the function,

$f(x)=\sin{\bigg{(}\frac{1}{x+0.01}\bigg{)}},\text{ for
}x\in\left[0,\frac{\pi}{4}\right].$ (6.2)

In the experiment we use a training set of $n=2.2\times 10^{6}$ samples and a
test set of $1.3\times 10^{5}$ samples. We use the non-uniform gamma
distribution $\Gamma(\alpha,\beta)$ with $\alpha=1,\,\beta=2$ to sample the
training data.

The number of training points used at each level in StreaMRAK is determined by
setting $\delta_{1}$ and $\delta_{2}$ from Def. 4.1 to $10^{-3}$. With this
choice, StreaMRAK selects between $30\,244\text{\,}$ and $40\,100\text{\,}$
training points for each level. For comparison, FALKON uses all the
$2.2\text{\times}{10}^{6}\text{\,}$ training points. Furthermore, for LP-KRR
we run two experiments: LP-KRR (1) using $1.1\text{\times}{10}^{5}\text{\,}$
training points at each level and LP-KRR (2) using
$2.2\text{\times}{10}^{6}\text{\,}$ training points at each level.

Table 1: Comparison of StreaMRAK, LP-KRR and FALKON for the target in Eq. (6.2). For each level $l$ we show the number of landmarks, the mean square error (MSE), and the accumulated time to train the prediction model (Time). In parenthesis, in the time column of the FALKON row, is the time to find the optimal bandwidth through cross-validation. | Level | # Landmarks | MSE | Time
---|---|---|---|---
StreaMRAK | 5 | $47\text{\,}$ | $2.55\text{\times}{10}^{-1}\text{\,}$ | $77\text{\,}\mathrm{s}$
10 | $392\text{\,}$ | $3.69\text{\times}{10}^{-2}\text{\,}$ | $116\text{\,}\mathrm{s}$
15 | $1525\text{\,}$ | $8.63\text{\times}{10}^{-6}\text{\,}$ | $497\text{\,}\mathrm{s}$
16 | $2302\text{\,}$ | $6.18\text{\times}{10}^{-6}\text{\,}$ | $1194\text{\,}\mathrm{s}$
LP-KRR (1) $n_{l}=$1.1\text{\times}{10}^{5}\text{\,}$$ | 5 | $1483\text{\,}$ | $2.56\text{\times}{10}^{-1}\text{\,}$ | $143\text{\,}\mathrm{s}$
10 | $1483\text{\,}$ | $3.65\text{\times}{10}^{-2}\text{\,}$ | $413\text{\,}\mathrm{s}$
15 | $1483\text{\,}$ | $8.72\text{\times}{10}^{-6}\text{\,}$ | $825\text{\,}\mathrm{s}$
16 | $1483\text{\,}$ | $6.85\text{\times}{10}^{-6}\text{\,}$ | $922\text{\,}\mathrm{s}$
| 18 | $1483\text{\,}$ | $6.55\text{\times}{10}^{-6}\text{\,}$ | $1136\text{\,}\mathrm{s}$
LP-KRR (2) $n_{l}=$2.2\text{\times}{10}^{6}\text{\,}$$ | 5 | $1483\text{\,}$ | $2.56\text{\times}{10}^{-1}\text{\,}$ | $2963\text{\,}\mathrm{s}$
10 | $1483\text{\,}$ | $3.64\text{\times}{10}^{-2}\text{\,}$ | $8704\text{\,}\mathrm{s}$
18 | $1483\text{\,}$ | $8.91\text{\times}{10}^{-6}\text{\,}$ | $23\,113\text{\,}\mathrm{s}$
StreaMRAK | – | $14\,830\text{\,}$ | $5.7\text{\times}{10}^{-3}\text{\,}$ |  $4642\text{\,}\mathrm{s}$+($27\,930\text{\,}\mathrm{s}$)

(a)

(b)

(c)

(d)

Figure 3: (a)-(d) shows the target function $f(x)$ from Eq. (6.2) as a grey
dotted line. The light-blue circles indicates the predicted values made by
StreaMRAK. Similarly the red triangles indicates the predictions made by LP-
KRR and the dark blue squares the predictions made by FALKON.

Results are presented in Table 1, and the prediction results are illustrated
in Fig.
LABEL:subfig:varsinus_streamrak_pred_lvl3-LABEL:subfig:varsinus_falkon_pred.
The results show that StreaMRAK and both LP-KRR schemes perform much better
than FALKON. The reason is that FALKON uses only one bandwidth $r$, while the
multi-resolution schemes StreaMRAK and LP-KRR, utilize a bandwidth regime
$r_{l}=2^{-l}r_{0}$ that varies with the level $l$. The consequence is that
StreaMRAK and LP-KRR approximate the low-frequency components of $f$ when the
bandwidth is large, and then target the high-frequency components of $f(x)$
gradually as the bandwidth decreases. These results illustrate the benefits of
a multi-resolution scheme over a single bandwidth scheme.

From Table 1, we also observe that LP-KRR (2) is significantly slower than
StreaMRAK and LP-KRR (1). This is because it uses the entire training set at
each level. Therefore, since LP-KRR (1) and LP-KRR (2) achieve comparable
precision, we see that including all training points at each level is not
always necessary.

A closer comparison of StreaMRAK and LP-KRR is given in Fig. 4. In particular,
in Fig. LABEL:subfig:varsin_mse_vs_level we see that the two algorithms
achieve very similar precision. However, comparing the training times in Fig.
LABEL:subfig:varsin_mse_vs_time, we see that StreaMRAK trains each level
faster and therefore achieves better precision earlier than LP-KRR (1).

(a)

(b)

Figure 4: Comparison of StreaMRAK and LP-KRR. (a) shows the mean square error
calculated according to Eq. (6.1) with the target function from Eq. (6.2).
Along the x-axis is the number of levels included in the model. (b) The x-axis
shows the accumulated training time until a level in the LP is completed. The
y-axis shows the MSE of the prediction using the currently available model.
The blue circles indicate the prediction error of StreaMRAK and the red
triangles indicate the prediction error of LP-KRR (1).

(a)

(b)

(c)

Figure 5: (a)-(c) shows the landmarks with their position along the $x_{1}$
axis and the average distance to their 2 nearest neighbors along the y-axis.
Here the red triangles are the Nyström landmarks of LP-KRR and the light-blue
circles the landmarks of StreaMRAK . The grey dotted line is the bandwidth at
the given level.

In Fig. 5 we show the average distance of each landmark to their 2 nearest
neighbors (2-NN distance). Two aspects of the selection require attention. As
opposed to LP-KRR, StreaMRAK selects landmarks such that the 2-NN distance is
comparable to the bandwidth used at a specific level. In addition, StreaMRAK
saves computational power by not choosing landmarks in regions where the 2-NN
distance is too low compared to the bandwidth. In Fig.
LABEL:subfig:avg_nn_lvl_16_streamrak_mrfalkon this can be observed for level
$l=16$ for landmarks with $x\geq 0.2$. Due to the non-uniform sample
distribution with a higher density around $x=0$, the adaptive sub-sampling is
able to select more landmarks in the region close to $x=0$, where $f$
oscillates with high frequency. Furthermore, StreaMRAK stops predicting at
level $16$ because level $17$ is not yet covered with a high enough density of
landmarks. Meanwhile, LP-KRR continues, but as seen from Fig.
LABEL:subfig:varsin_mse_vs_level the improvements after level $15$ are not
significant because the density of Nyström samples is too low compared to the
bandwidth.

###  6.2 Adaptive sub-sampling benchmark

We consider a dumbbell-shaped domain embedded in $\mathbb{R}^{5}$, consisting
of two 5-dimensional spheres connected by a 2-dimensional plane. A projection
of the input domain in $\mathbb{R}^{3}$ is shown in Fig. 7 (a)-(c).
Furthermore, as target we consider the following function,

$f({\bm{\mathsf{x}}})=\begin{cases}A\sin(Bx_{1}+\phi)+(x_{1}+2),&1<x_{1}<3\\\
1,&\text{ otherwise }\\\ \end{cases},\text{ for
}{\bm{\mathsf{x}}}\in[-1,5]\times[-1,1]^{4},$ (6.3)

where $A,B$ and $\phi$ are chosen so that $f\in{\cal
C}^{1}([-1,5]\times[-1,1]^{4},\mathbb{R}^{5})$. For the experiments, we
consider a training set of $1.9\text{\times}{10}^{6}\text{\,}$ samples and a
test set of $6\text{\times}{10}^{5}\text{\,}$ samples, all sampled uniformly
at random from the input domain. We note that we purposefully chose a simple
function in the high dimensional regions because complicated functions in high
dimensions require far too many points to be satisfactorily learned.

To determine the number of training points for StreaMRAK, we let
$\delta_{1}=$1\text{\times}{10}^{-3}\text{\,}$$ and
$\delta_{2}=$1\text{\times}{10}^{-4}\text{\,}$$, cf. Def. 4.1. With this
choice StreaMRAK selects between $30\,100\text{\,}$ and $40\,100\text{\,}$
training points for each level. FALKON again uses all the
$1.9\text{\times}{10}^{6}\text{\,}$ training points and for LP-KRR we consider
two settings: LP-KRR (1) using $1.8\text{\times}{10}^{5}\text{\,}$ training
points at each level, and LP-KRR (2) using $1.9\text{\times}{10}^{6}\text{\,}$
training points at each level.

The results for StreaMRAK, LP-KRR, and FALKON are presented in Table 2. We
observe that StreaMRAK achieves a better prediction than both FALKON and LP-
KRRbecause it adapts the sub-sampling density to the level of resolution.

Table 2: Comparison of StreaMRAK, LP-KRR, and FALKON predictions of the target function in Eq. (6.3). For each level $l$ we show the number of landmarks, the mean square error (MSE), and the accumulated time to train the prediction model (Time). In parenthesis, in the time column of the FALKON row, is the time to find the optimal bandwidth through cross-validation. | Level | # Landmarks | MSE | Time
---|---|---|---|---
StreaMRAK | 4 | $352\text{\,}$ | $1.29\text{\times}{10}^{-3}\text{\,}$ | $64\text{\,}\mathrm{s}$
5 | $2667\text{\,}$ | $1.27\text{\times}{10}^{-3}\text{\,}$ | $1398\text{\,}\mathrm{s}$
6 | $1858\text{\,}$ | $8.31\text{\times}{10}^{-4}\text{\,}$ | $1462\text{\,}\mathrm{s}$
8 | $1329\text{\,}$ | $2.75\text{\times}{10}^{-5}\text{\,}$ | $2307\text{\,}\mathrm{s}$
LP-KRR (1) $n_{l}=$1.8\text{\times}{10}^{5}\text{\,}$$ | 4 | $1375\text{\,}$ | $1.28\text{\times}{10}^{-3}\text{\,}$ | $386\text{\,}\mathrm{s}$
5 | $1375\text{\,}$ | $1.26\text{\times}{10}^{-3}\text{\,}$ | $520\text{\,}\mathrm{s}$
6 | $1375\text{\,}$ | $9.10\text{\times}{10}^{-4}\text{\,}$ | $671\text{\,}\mathrm{s}$
8 | $1375\text{\,}$ | $3.30\text{\times}{10}^{-4}\text{\,}$ | $1064\text{\,}\mathrm{s}$
| 9 | $1375\text{\,}$ | $3.10\text{\times}{10}^{-4}\text{\,}$ | $1287\text{\,}\mathrm{s}$
LP-KRR (2) $n_{l}=$1.9\text{\times}{10}^{6}\text{\,}$$ | 4 | $1375\text{\,}$ | $1.34\text{\times}{10}^{-3}\text{\,}$ | $4160\text{\,}\mathrm{s}$
5 | $1375\text{\,}$ | $1.30\text{\times}{10}^{-3}\text{\,}$ | $5570\text{\,}\mathrm{s}$
6 | $1375\text{\,}$ | $9.44\text{\times}{10}^{-4}\text{\,}$ | $7168\text{\,}\mathrm{s}$
8 | $1375\text{\,}$ | $3.16\text{\times}{10}^{-4}\text{\,}$ | $11\,125\text{\,}\mathrm{s}$
| 9 | $1375\text{\,}$ | $3.01\text{\times}{10}^{-4}\text{\,}$ | $13\,334\text{\,}\mathrm{s}$
FALKON | – | $14\,830\text{\,}$ | $6.8\text{\times}{10}^{-4}\text{\,}$ |  $6590\text{\,}\mathrm{s}$+($37\,561\text{\,}\mathrm{s}$)

To understand the improvement in prediction accuracy, we need to discuss the
effects of landmark selection. In Fig.
LABEL:subfig:dumbell_lm_dist_lvl4-LABEL:subfig:dumbell_lm_dist_lvl6 we show
the projections of landmarks for StreaMRAK and LP-KRR on $\mathbb{R}^{3}$, and
in Fig.
LABEL:subfig:dumbell_lm_avg_distance_lvl4-LABEL:subfig:dumbell_lm_avg_distance_lvl6
the average distance of each landmark to its $7$ nearest neighbors. These
distances are compared with the bandwidth $r_{l}$ selected for the given level
$l$. We see that StreaMRAK selects landmarks in regions where the average
distance to nearest neighbors is comparable to the bandwidth. This means that
in high dimensional regions, which correspond to $x_{1}\in[-1,1]\cup[3,5]$,
the algorithm effectively stops collecting landmarks since it cannot maintain
high enough density. On the other hand, LP-KRR uses Nystrom sub-sampling,
which imposes a uniform selection of landmarks. Consequently, a significant
number of landmarks come from high-dimensional regions.

(a)

(b)

Figure 6: Comparison of StreaMRAK and LP-KRR (1) in the 2-dim and 5-dim
regions of the Dumbbell domain. The solid blue line is StreaMRAK for dimension
$d=2$ while the solid red line is LP-KRR (1) for dimension $d=2$. The grey
dotted line is StreaMRAK for dimension $d=5$ and the dark-grey dashed line is
LP-KRR (1) for dimension $d=5$ (a) shows the mean square error calculated
according to Eq. (6.1). (b) shows the number of landmarks in the 2-dimensional
and the 5- dimensional regions.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7: In the figure, red triangles correspond to LP-KRR and light-blue
circles to StreaMRAK. (a)-(c) shows the landmark distributions projected on
$\mathbb{R}^{3}$ at level $l=4,5,6$ respectively. (d)-(f) shows the average
distance between the 7 nearest neighbors; bandwidth $r_{l}$ is indicated with
a dotted line.

Moreover, Fig. 7 shows that in the case of LP-KRR, the average distance
between the landmarks in high dimensional regions is larger than the bandwidth
$r_{l}$ when $l\geq 5$. As a knock-on effect, LP-KRR makes only small
improvements in high dimensional regions for $l\geq 5$, as seen from Fig.
LABEL:fig:mse_low_and_high_dimensions. Analogous behavior can be observed for
StreaMRAK. However, since StreaMRAK devotes fewer resources to high
dimensional regions, it sub-samples more from the low dimensional region, as
illustrated in Fig. LABEL:fig:comparison_low_and_high_dim_landmarks. The
consequence is that StreaMRAK makes bigger improvements in the low dimensional
region than LP-KRR, as seen from Fig. LABEL:fig:mse_low_and_high_dimensions.
Note that this was not the case in Section 6.1, where the two methods had
similar behavior, but unlike here, the input domain in Section 6.1 did not
consist of regions with different dimensionalities.

###  6.3 Forecasting the trajectory of a double pendulum

We consider the double pendulum, illustrated in Fig.
LABEL:fig:DoublePendIllustration, which we model by the Lagrangian system

${\cal
L}=ml^{2}(\omega_{1}^{2}+\frac{1}{2}\omega_{2}^{2})+ml^{2}\omega_{1}\omega_{2}\cos{(\theta_{1}-\theta_{2})}+mgl(2\cos{\theta_{1}}+\cos{\theta_{2}}),$
(6.4)

under the assumption that the pendulums are massless rods of length
$l_{1}=l_{2}=l$ with masses $m_{1}=m_{2}=m$ centered at the end of each rod.
Here $g$ is the standard gravity, $\omega_{1}\coloneqq\dot{\theta}_{1}$,
$\omega_{2}\coloneqq\dot{\theta}_{2}$ are the angular velocities, and the
angles $\theta_{1}$, $\theta_{2}$ are as indicated in Fig.
LABEL:fig:DoublePendIllustration. For the experiments we let $m=1$, $l=1$ and
$g=10$.

The learning task is to forecast the trajectory of the pendulum, given only
its initial conditions. We let
${\bm{\mathsf{s}}}_{t}=[\theta_{1}(t),\,\theta_{2}(t),\,\omega_{1}(t),\,\omega_{2}(t)]\in\mathbb{R}^{4}$
be the state of the system at step $t\in\mathbb{N}$ and train StreaMRAK, LP-
KRR and FALKON to learn how ${\bm{\mathsf{s}}}_{t}$ maps to a later state
${\bm{\mathsf{s}}}_{t+\Delta}$, for $\Delta\in\mathbb{N}$. The trained model
$\widehat{f}$ is used to forecast the state ${\bm{\mathsf{s}}}_{T}$ for $T>>0$
by recursively predicting
${\bm{\mathsf{s}}}_{t+\Delta}=\widehat{f}({\bm{\mathsf{s}}}_{t})$ from the
initial state ${\bm{\mathsf{s}}}_{0}$ until $t=T$.

For the experiments we consider two settings: a low energy system
${\bm{\mathsf{s}}}_{0}^{\textit{low}}=[-20\degree,\,-20\degree,\,0\degree,\,0\degree]$
and a high energy system
${\bm{\mathsf{s}}}_{0}^{\textit{high}}=[-120\degree,\,-20\degree,\,-7.57\degree,\,7.68\degree]$.
For these systems, we initialize $8000$ pendulums as
${\bm{\mathsf{s}}}_{0}\sim{\cal
N}({\bm{\mathsf{s}}},\sigma({\bm{\mathsf{s}}}))$ for
${\bm{\mathsf{s}}}={\bm{\mathsf{s}}}_{0}^{\textit{low}},\,{\bm{\mathsf{s}}}_{0}^{\textit{high}}$
respectively, where
$\sigma({\bm{\mathsf{s}}})=[0.025|\theta_{1}|,\,0.15|\theta_{2}|,\,0.3|\omega_{1}|,\,0.3|\omega_{2}|]$.
Each pendulum is iterated for $500$ steps, which results in $5\times 10^{6}$
training points distributed in $\mathbb{R}^{4}$. Furthermore, for the test
data we consider $100$ pendulums ${\bm{\mathsf{s}}}_{0}\sim{\cal
N}({\bm{\mathsf{s}}},0.01|{\bm{\mathsf{s}}}|)$ for
${\bm{\mathsf{s}}}={\bm{\mathsf{s}}}_{0}^{\textit{low}},\,{\bm{\mathsf{s}}}_{0}^{\textit{high}}$,
iterated for 500 steps.

To determine the number of training points for StreaMRAK, we let
$\delta_{1},\delta_{2}=10^{-4}$, cf. Def. 4.1. With this choice StreaMRAK
selects between $30\,219\text{\,}$ and $70\,282\text{\,}$ training points for
each level for the low energy system, and between $36\,300\text{\,}$ and
$130\,200\text{\,}$ for the high energy system. Meanwhile, FALKON uses all
$5.0\text{\times}{10}^{6}\text{\,}$ training points and LP-KRR use
$3.9\text{\times}{10}^{5}\text{\,}$ training points at each level.

Results are presented in Table 3 and 4. Furthermore, to illustrate the
prediction results we consider the center of mass
$\overline{M}_{x}({\bm{\mathsf{s}}}_{t})=\frac{1}{2}(x_{1}(t)+x_{2}(t))\in\mathbb{R}$
at state ${\bm{\mathsf{s}}}_{t}$, where $x_{1},\,x_{2}\in\mathbb{R}$ are the
positions of the two pendulum masses as seen in Fig.
LABEL:fig:DoublePendIllustration. The prediction results are illustrated in
Fig. 8 and 9 for the low and high energy pendulums respectively. We calculate
the MSE at each step $t$ separately, such that for a given $t$ we use Eq. 6.1
with ${\bm{\mathsf{y}}}_{k}=\overline{M}_{x}({\bm{\mathsf{s}}}_{t})$,
${\bm{\mathsf{y}}}^{\textit{pred}}_{k}=\overline{M}_{x}({\bm{\mathsf{s}}}^{\textit{pred}}_{t})$
and $\Upsilon=100$.

Table 3: Comparison of StreaMRAK, LP-KRR, and FALKON for the low energy system. For each level $l$ we show the number of landmarks, the MSE at step T=50, and the accumulated time to train the prediction model (Time). In parenthesis, in the time column of the FALKON row, is the time to find the optimal bandwidth through cross-validation. | Level | # Landmarks | MSE(T=50) | Time
---|---|---|---|---
StreaMRAK | 2 | $1\text{\,}$ | $1.12\text{\times}{10}^{-1}\text{\,}$ | $31\text{\,}\mathrm{s}$
5 | $66\text{\,}$ | $3.39\text{\times}{10}^{-5}\text{\,}$ | $498\text{\,}\mathrm{s}$
7 | $659\text{\,}$ | $1.44\text{\times}{10}^{-6}\text{\,}$ | $534\text{\,}\mathrm{s}$
9 | $4085\text{\,}$ | $3.01\text{\times}{10}^{-7}\text{\,}$ | $812\text{\,}\mathrm{s}$
LP-KRR | 2 | $1979\text{\,}$ | $5.93\text{\times}{10}^{-2}\text{\,}$ | $490\text{\,}\mathrm{s}$
5 | $1979\text{\,}$ | $2.73\text{\times}{10}^{-5}\text{\,}$ | $1463\text{\,}\mathrm{s}$
7 | $1979\text{\,}$ | $2.16\text{\times}{10}^{-7}\text{\,}$ | $2395\text{\,}\mathrm{s}$
9 | $1979\text{\,}$ | $1.16\text{\times}{10}^{-8}\text{\,}$ | $3550\text{\,}\mathrm{s}$
FALKON | – | $19\,790\text{\,}$ | $5\text{\times}{10}^{6}\text{\,}$ |  $2934\text{\,}\mathrm{s}$+($1498\text{\,}\mathrm{s}$)
Table 4: Comparison of the StreaMRAK, LP-KRR, and FALKON for the high energy system. For each level $l$ we show the number of landmarks, the MSE at step T=50, and the accumulated time to train the prediction model (Time). In parenthesis, in the time column of the FALKON row, is the time to find the optimal bandwidth through cross-validation. | Level | # Landmarks | MSE(T=50) | Time
---|---|---|---|---
StreaMRAK | 2 | $1\text{\,}$ | $2.70\text{\times}{10}^{-1}\text{\,}$ | $49\text{\,}\mathrm{s}$
5 | $1106\text{\,}$ | $8.53\text{\times}{10}^{-3}\text{\,}$ | $915\text{\,}\mathrm{s}$
7 | $6376\text{\,}$ | $2.16\text{\times}{10}^{-4}\text{\,}$ | $1999\text{\,}\mathrm{s}$
LP-KRR | 2 | $1979\text{\,}$ | $1.72\text{\times}{10}^{-2}\text{\,}$ | $522\text{\,}\mathrm{s}$
5 | $1979\text{\,}$ | $5.09\text{\times}{10}^{-3}\text{\,}$ | $1474\text{\,}\mathrm{s}$
7 | $1979\text{\,}$ | $1.39\text{\times}{10}^{-4}\text{\,}$ | $2431\text{\,}\mathrm{s}$
FALKON | – | $19\,790\text{\,}$ | $5\text{\times}{10}^{6}\text{\,}$ |  $23\,830\text{\,}\mathrm{s}$+($11\,050\text{\,}\mathrm{s}$)

For the low energy system, we see from Fig. LABEL:subfig:mse_vs_time_dp_lowE
how StreaMRAK is trained significantly faster than LP-KRR, although at a cost
of reduced precision. The reduced training time of StreaMRAK is a consequence
of the low doubling dimension of the training data, which allows the selection
of far fewer landmarks for StreaMRAK than what is used at each level in LP-
KRR.

For the high-energy pendulum, we see from Fig.
LABEL:subfig:dp_mse_vs_time_highE that StreaMRAK is again able to achieve good
precision faster than LP-KRR . Furthermore, we see that the number of
landmarks selected for StreaMRAK increases abruptly with the levels,
reflecting the high doubling dimension of the training data. Due to this
StreaMRAK stops the training after level $7$, as the next levels require too
many landmarks. By continuing for $2$ more levels LP-KRR is able to achieve
marginally better precision but at an increased computational cost.

(a)

(b)

(c)

Figure 8: Comparison of StreaMRAK (light blue lines and circles), LP-KRR (red
lines and triangles), and FALKON (dark blue dotted lines and squares) for the
low energy pendulum. (a) Shows the mean square error of the center of mass
$\overline{M}_{x}({\bm{\mathsf{s}}}_{t})$ for the level $7$ prediction, with
step $T$ along the x-axis. (b) shows the true center of mass trajectory as a
grey dotted line and the predictions of StreaMRAK, LP-KRR, and FALKON at level
$9$. (c) The x-axis shows the accumulated training time until a level in the
LP is completed. The y-axis shows the MSE of the predicted system state after
$T=50$ steps. We note that StreaMRAK includes $7$ levels, while LP-KRR
includes $9$.

(a)

(b)

(c)

Figure 9: Comparison of StreaMRAK (light blue lines and circles), LP-KRR (red
lines and triangles), and FALKON (dark blue dotted lines and squares) for the
high-energy pendulum. (a) Shows the mean square error of the center of mass
$\overline{M}_{x}({\bm{\mathsf{s}}}_{t})$ for the level $7$ prediction, with
step $T$ along the x-axis. (b) shows the true center of mass trajectory as a
grey dotted line and the predictions of StreaMRAK, LP-KRR, and FALKON at level
$9$. (c) The x-axis shows the accumulated training time until a level in the
LP is completed. The y-axis shows the MSE of the predicted system state after
$T=50$ steps.

As seen in Fig. LABEL:subfig:cm_pred_doublePend_highE, the forecasting of
StreaMRAK and LP-KRR breaks down after $T\approx 200$ steps. In Fig.
LABEL:subfig:phase_diagram_with_bifurcation we observe the trajectory of a
pendulum with initial condition ${\bm{\mathsf{s}}}_{0}^{\textit{high}}$, as
well as four pendulums with a $0.5\%$ perturbation on the angles $\theta_{1}$
and $\theta_{2}$ in ${\bm{\mathsf{s}}}_{0}^{\textit{high}}$. We observe that
after roughly $T=205$ time steps the trajectory of the five pendulums diverge
significantly from each other. Therefore, it seems that a bifurcation point
occurs around this time, which may explain why all the algorithms are unable
to make good forecasting beyond this point.

(a)

(b)

(c)

Figure 10: (a) Pendulum positions at $T=200$ and $(c)$ The positions at
$T=210$. In (a) and (c), $P_{M}$ is the main pendulum with initial conditions
${\bm{\mathsf{s}}}_{0}^{\textit{high}}$, while $P_{F}$ is the
StreaMRAKforecast of the pendulum position. Similarly, $P_{0}$-$P_{3}$ are
four training pendulums with a perturbation of $0.5\%$ on the initial angles
$\theta_{1}$ and $\theta_{2}$ of the main pendulum. (b) Projection of the
training data on the $\theta_{1}\theta_{2}$-plane. The thick red line is the
main pendulum corresponding to $P_{M}$ and the four grey dotted lines are the
test pendulums $P_{0}$-$P_{3}$, where the X indicates the time $T=205$.

##  7 Outlook

Further development of StreaMRAK is intended with focus on four objectives.

  1. (O1)

Augmentation of the DCT to track the error at each node

  2. (O2)

Improve the estimator in Def. 4.1 and Eq. 3.1.

  3. (O3)

Refinement of previously fitted levels in the LP as new data arrives.

  4. (O4)

Further theoretical analysis of the LP.

Considering objective (O1) we intend to develop the DCT to track the error at
each node. This way the growth can be restricted in regions where the error is
small, which allows for more focus on regions where the error is large. The
intention is that this will reduce the problem complexity even further, while
also increasing the precision. Regarding objective (O2), a drawback with the
estimator in Eq. 3.1 was already mentioned in Remark B.1. Furthermore, for the
estimator in Def. 4.1, we intend to implement and evaluate alternative ways to
estimate the convergence of the matrices. Another focus area will be objective
(O3), as we believe new information may be revealed as new training data
arrive, and refinement of previously fitted levels can therefore be
beneficial. Finally, the theoretical analysis in objective (O4) will focus on
analyzing the generalization error for the LP, particularly in combination
with the adaptive sub-sampling scheme.

##  8 Acknowledgement

We especially would like to thank Prof. Pieter Abeel at UC Berkeley and Asst.
Prof. Sicun Gao at UC San Diego for their input on the double pendulum system,
and for providing a code example for this system. We would also like to thank
Sami Ortoleva at UC San Diego for his discussion on the analysis of the damped
cover-tree. AO is part of the Simula-UCSD-UiO Research and Ph.D. training
program (SUURPh), an international collaboration in computational biology and
medicine funded by the Norwegian Ministry of Education and Research, ŽK is
funded by UK EPSRC grant EP/T000864/1, AC is funded by NSF DMS 1819222,
2012266, and Russell Sage Foundation grant 2196 and YF is funded by the NIH
grant NINDS (PHS) U19NS107466 Reverse Engineering the Brain Stem Circuits that
Govern Exploratory Behavior.

## References

  * [1] Bernhard Schölkopf and Alexander J Smola.  Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond.  MIT press, 1 edition, 2002.
  * [2] Alessandro Rudi, Luigi Carratino, and Lorenzo Rosasco.  FALKON: An optimal large scale kernel method.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Proc. 31th Int. Conf. Neural Inf. Process. Syst., volume 30, pages 3889–3899, 2017.
  * [3] Ahmed El Alaoui and Michael W. Mahoney.  Fast randomized kernel ridge regression with statistical guarantees.  In Proc. 28th Int. Conf. Neural Inf. Process. Syst., volume 1, pages 775–783, 2015.
  * [4] Yuchen Zhang, John Duchi, and Martin Wainwright.  Divide and conquer kernel ridge regression: A distributed algorithm with minimax optimal rates.  J. Mach. Learn. Res., 16:3299–3340, 2015.
  * [5] Haim Avron, Kenneth L. Clarkson, and David P. Woodruff.  Faster kernel ridge regression using sketching and preconditioning.  J. Matrix. Anal. Appl., 38(4):1116–1138, 2017.
  * [6] Evgeny Burnaev and Ivan Nazarov.  Conformalized kernel ridge regression.  In Proc. 15th Int. Conf. Mach. Learn. Appl., pages 45–52, 2017\.
  * [7] Peter Exterkate, Patrick J.F. Groenen, Christiaan Heij, and Dick van Dijk.  Nonlinear forecasting with many predictors using kernel ridge regression.  Int. J. Forecas, 32(3):736–753, 2016.
  * [8] Mu Niu, Simon Rogers, Maurizio Filippone, and Dirk Husmeier.  Fast parameter inference in nonlinear dynamical systems using iterative gradient matching.  In Proc. 33rd Int. Conf. Mach. Learn. Res., pages 1699–1707, 2016\.
  * [9] Michiel Stock, Tapio Pahikkala, Antti Airola, Bernard De Baets, and Willem Waegeman.  A comparative study of pairwise learning methods based on kernel ridge regression.  Neural Comput., 30(8):2245–2283, 2018.
  * [10] Senjian An, Wanquan Liu, and Svetha Venkatesh.  Face recognition using kernel ridge regression.  In Proc. Conf. Comput. Vis. Recognit., pages 1–7, 2007.
  * [11] Benjamin Yee Shing Li, Lam Fat Yeung, and King Tim Ko.  Indefinite kernel ridge regression and its application on QSAR modelling.  Neurocomputing, 158:127–133, 2015.
  * [12] P. Mohapatra, S. Chakravarty, and P. K. Dash.  Microarray medical data classification using kernel ridge regression and modified cat swarm optimization based gene selection system.  Swarm. Evol. Comput., 28:144–160, 2016.
  * [13] S. Muthukrishnan.  Data streams: Algorithms and applications.  Found. Trends Theor. Comput. Sci., 1(2):117–236, 2005.
  * [14] Wei Fan and Albert Bifet.  Mining big data.  ACM SIGKDD Explor. Newsl., 14(2):1–5, 2013.
  * [15] Kun Lan, Dan-Tong Wang, Simon Fong, Lian-Sheng Liu, Kelvin K L Wong, and Nilanjan Dey.  A survey of data mining and deep learning in bioinformatics.  J. Med. Syst., 42(8), 2018.
  * [16] Jyrki Kivinen, Alex J. Smola, and Robert C. Williamson.  Online learning with kernels.  In Proc. 14th Int. Conf. Neural Inf. Process. Syst., page 785–792, 2001.
  * [17] Clint Scovel, Don Hush, Ingo Steinwart, and James Theiler.  Radial kernels and their reproducing kernel hilbert spaces.  J. Complex., 26(6):641–660, 2010.
  * [18] Charles A. Micchelli, Yuesheng Xu, and Haizhang Zhang.  Universal kernels.  J. Mach. Learn. Res., 7:2651–2667, 2006.
  * [19] Zhuang Wang, Koby Crammer, and Slobodan Vucetic.  Breaking the curse of kernelization: Budgeted stochastic gradient descent for large-scale SVM training.  J. Mach. Learn. Res., 13:3103–3131, 2012.
  * [20] Clive R Loader.  Bandwidth selection: Classical or plug-in?  Ann. Stat., 27(2):415–438, 1999.
  * [21] Gavin C. Cawley and Nicola L.C. Talbot.  Fast exact leave-one-out cross-validation of sparse least-squares support vector machines.  Neural Netw., 17(10):1467–1475, 2004.
  * [22] Sylvain Arlot and Alain Celisse.  A survey of cross-validation procedures for model selection.  Stat. Surv., 4:40–79, 2010.
  * [23] Gert R.G. Lanckriet, Nello Cristianini, Peter Bartlett, Laurent El Ghaoui, and Michael I. Jordan.  Learning the kernel matrix with semidefinite programming.  J. Mach. Learn. Res., 5:27–72, 2004.
  * [24] Francis R. Bach, Gert R.G. Lanckriet, and Michael I. Jordan.  Multiple kernel learning, conic duality, and the SMO algorithm.  In Proc. 21th Int. Conf. Mach. Learn., pages 41–48, 2004.
  * [25] Sören Sonnenburg, Gunnar Rätsch, Christin Schäfer, and Bernhard Schölkopf.  Large scale multiple kernel learning.  J. Mach. Learn. Res., 7:1531–1565, 2006.
  * [26] Eduard Gabriel Băzăvan, Fuxin Li, and Cristian Sminchisescu.  Fourier kernel learning.  In Eur. Conf. Comput. Vis., pages 459–473, 2012.
  * [27] Amit Bermanis, Amir Averbuch, and Ronald R. Coifman.  Multiscale data sampling and function extension.  Appl. Comput. Harmon. Anal., 34(1):15–29, 2013.
  * [28] N. Rabin and D. Fishelov.  Multi-scale kernels for Nyström based extension schemes.  Appl. Math. Comput., 319:165–177, 2018.
  * [29] Neta Rabin and Ronald R. Coifman.  Heterogeneous datasets representation and learning using diffusion maps and Laplacian pyramids.  In Proc. 12th Int. Conf. Data Min., pages 189–199, 2012.
  * [30] William Leeb.  Properties of Laplacian pyramids for extension and denoising.  arXiv preprint, 2019.
  * [31] Amara Graps.  An introduction to wavelets.  IEEE Comput. Sci. Eng., 2(2):50–61, 1995.
  * [32] Ali N. Akansu, Wouter A. Serdijn, and Ivan W. Selesnick.  Emerging applications of wavelets: A review.  Phys. Commun., 3(1):1–18, 2010.
  * [33] Ronald R. Coifman and Mauro Maggioni.  Diffusion wavelets.  Appl. Comput. Harmon. Anal., 21(1):53–94, 2006.
  * [34] M. Maggioni and H. N. Mhaskar.  Diffusion polynomial frames on metric measure spaces.  Appl. Comput. Harmon. Anal., 24(3):329–353, 2008.
  * [35] David K. Hammond, Pierre Vandergheynst, and Rémi Gribonval.  Wavelets on graphs via spectral graph theory.  Appl. Comput. Harmon. Anal., 30(2):129–150, 2011.
  * [36] Alexander Cloninger, Haotian Li, and Naoki Saito.  Natural graph wavelet packet dictionaries.  J. Fourier Anal. Appl., 27(3):1–33, 2021.
  * [37] Ernesto De Vito, Zeljko Kereta, Valeriya Naumova, Lorenzo Rosasco, and Stefano Vigogna.  Wavelet frames generated by a reproducing kernel.  J. Fourier Anal. Appl., 27(2):1–39, 2021.
  * [38] Robert Krauthgamer and James R. Lee.  Navigating nets: Simple algorithms for proximity search.  In Proc. 15th Annu. ACM-SIAM Symp. Discrete Algorithms, pages 798–807, 2004.
  * [39] Alina Beygelzimer, Sham Kakade, and John Langford.  Cover trees for nearest neighbor.  In Proc. 23th Int. Conf. Mach. Learn., pages 97–104, 2006.
  * [40] Troy Shinbrot, Celso Grebogi, Jack Wisdom, and James A Yorke.  Chaos in a double pendulum.  Am. J. Phys., 60(6):491–499, 2016.
  * [41] Angelo Marcelo Tusset, Vinícius Piccirillo, Atila Madureira Bueno, José Manoel Balthazar, Danuta Sado, Jorge Luis Palacios Felix, and Reyolando Manoel Lopes Rebello Da Fonseca Brasil.  Chaos control and sensitivity analysis of a double pendulum arm excited by an RLC circuit based nonlinear shaker.  J. Vib. Control, 22(17):3621–3637, 2016.
  * [42] Peter J Burt and Edward H Adelson.  The Laplacian pyramid as a compact image code.  IEEE Trans. commun., 31(4):532–540, 1983.
  * [43] Ali Rahimi and Benjamin Recht.  Random features for large-scale kernel machines.  In Adv. Neural Inf. Process Syst., volume 20, pages 1177–1184, 2008\.
  * [44] Quoc V Le, Tamas Sarlos, and Alex Smola.  Fastfood-computing hilbert space expansions in loglinear time.  In Proc. 30th Int. Conf. Mach. Learn., volume 28, page 244–252, 2013.
  * [45] Zichao Yang, Alexander J. Smola, Le Song, and Andrew Gordon Wilson.  A la carte — learning fast kernels.  In Proc. 18th Int. Conf. Artif. Intell. Stat., volume 38, pages 1098–1106, 2015.
  * [46] Siyuan Ma and Mikhail Belkin.  Sigma-delta and distributed noise-shaping quantization methods for random fourier features.  arXiv preprint, 2021.
  * [47] Christopher Williams and Matthias Seeger.  Using the Nyström method to speed up kernel machines.  In Proc. 14th Annu. Conf. Neural Inf. Process Syst., volume 13, pages 682–688, 2001.
  * [48] Alex J. Smola and Bernhard Schölkopf.  Sparse greedy matrix approximation for machine learning.  In Proc. 1th Int. Conf. Mach. Learn., pages 911 – 918, 2000.
  * [49] Alexander Cloninger.  Prediction models for graph-linked data with localized regression.  In Proc. SPIE Int. Soc. Opt. Eng., volume 10394, 2017.
  * [50] Siyuan Ma, Raef Bassily, and Mikhail Belkin.  The power of interpolation : Understanding the effectiveness of SGD.  In Proc. 35th Int. Conf. Mach. Learn., page 3331–3340, 2018.
  * [51] Siyuan Ma and Mikhail Belkin.  Kernel machines that adapt to GPUs for effective large batch training.  arXiv preprint, 2018.
  * [52] Jiaming Zhang, Hanwen Ning, Xingjian Jing, and Tianhai Tian.  Online kernel learning with adaptive bandwidth by optimal control approach.  IEEE Trans. Neural Netw. Learn. Syst., 32(5):1920–1934, 2021.
  * [53] Badong Chen, Junli Liang, Nanning Zheng, and José C. Príncipe.  Kernel least mean square with adaptive kernel size.  Neurocomputing, 191:95–106, 2016.
  * [54] Haijin Fan, Qing Song, and Sumit B. Shrestha.  Kernel online learning with adaptive kernel width.  Neurocomputing, 175:233–242, 2015.
  * [55] Xiuyuan Cheng and Hau-Tieng Wu.  Convergence of graph Laplacian with kNN self-tuned kernels.  preprint arXiv, 2020.
  * [56] Jeffrey Picka.  Statistical inference for disordered sphere packings.  Stat. Surv., 6:74–112, 2012.
  * [57] M. Hifi and R. M’Hallah.  A literature review on circle and sphere packing problems: Models and methodologies.  Adv. Oper. Res., 2009.
  * [58] Philippe Flajolet, Danièle Gardy, and Loÿs Thimonier.  Birthday paradox, coupon collectors, caching algorithms and self-organizing search.  Discrete Appl. Math., 39(3):207–229, 1992.
  * [59] Gabriel Klambauer.  Problems and Propositions in Analysis.  Marcel Dekker, New York, 1979.
  * [60] Thomas Hofmann, Bernhard Schölkopf, and Alexander J. Smola.  Kernel methods in machine learning.  Ann. Stat., 36(3):1171–1220, 2008.
  * [61] MA Aiserman, Emmanuil M Braverman, and Lev I Rozonoer.  Theoretical foundations of the potential function method in pattern recognition.  Avtomat. i Telemeh., 25(6):917–936, 1964.
  * [62] Bernhard E. Boser, Isabelle M. Guyon, and Vladimir N. Vapnik.  A training algorithm for optimal margin classifiers.  In Proc. 5th Annu. Workshop. Comput. Learn. Theory, pages 144–152, 1992.
  * [63] George S. Kimeldorf and Grace Wahba.  A correspondence between bayesian estimation on stochastic processes and smoothing by splines.  Ann. Math. Stat., 41(2):495–502, 2011.
  * [64] Bernhard Schölkopf, Ralf Herbrich, and Alex J. Smola.  A generalized representer theorem.  In Int. Conf. Comput. Learn. Theory, pages 416–426, 2001.
  * [65] Holger Wendland.  Scattered Data Approximation.  Cambridge University Press, 2004.
  * [66] Danny Gómez.  A more direct proof of Gerschgorin´s theorem.  Mat: Enseñanza Univ., 14(2):119–122, 2006.

## Appendix A. Algorithms

We here denote nodes by $p,q,c$ and
${\bm{\mathsf{x}}}_{p},\,{\bm{\mathsf{x}}}_{q},{\bm{\mathsf{x}}}_{c}\in{\cal
X}\subset\mathbb{R}^{D}$ are the corresponding points.

Algorithm A.1 Insert(point $q$, node $p$, level $l$)

1: We assume $q$ already satisfies
$\|{\bm{\mathsf{x}}}_{q}-{\bm{\mathsf{x}}}_{p}\|\leq 2^{-l}r_{0}$.

2: if $\|{\bm{\mathsf{x}}}_{q}-{\bm{\mathsf{x}}}_{c}\|>2^{-(l+1)}r_{0}$ for
all $c\in\textit{Children}(p)$ then

3: Insert $q$ into $\textit{Children}(c)$.

4: Update_CoverFraction(Parent($Q_{l}$), ”No parent found”)

5: Break

6: else if $\|{\bm{\mathsf{x}}}_{q}-{\bm{\mathsf{x}}}_{c}\|<2^{-(l+1)}r_{0}$
for some $c\in\textit{Children}(p)$ then

7: Consider all children of $c$, namely $\textit{Children}(c)$

8: if $\textit{Children}(c)$ is empty then

9: if Covering fraction of $p$, Def. 3.1, satisfy
$\mathfrak{c}\mathfrak{f}(p)\geq{\cal D}_{\mathfrak{c}\mathfrak{f}}$ for some
threshold ${\cal D}_{\mathfrak{c}\mathfrak{f}}$ then

10: Insert q into $\textit{Children}(c)$

11: Break

12: else

13: Update_CoverFraction(p, ”parent found”) {c is found to be a potential
parent. However, since $\mathfrak{c}\mathfrak{f}(p)<{\cal
D}_{\mathfrak{c}\mathfrak{f}}$ we can not add $q$ to Children(c)}

14: end if

15: else

16: Insert($q$, $c$, $l+1$)

17: end if

18: end if

Algorithm A.2 StreaMRAK(point ${\bm{\mathsf{x}}}$, target $y$)

1: Let $l$ be the level. Let $p_{(0)}$ be the root node, $r_{0}$ the radius of
the root node.

2: Sub-sampling thread

3: Insert ${\bm{\mathsf{x}}}$ into the cover tree with
$\textsc{Insert}({\bm{\mathsf{x}}},p_{(0)},l=0)$. {See Alg. A.1}

4: if a new level has
$\mathfrak{c}\mathfrak{f}(Q_{l})\geq\mathcal{D}_{level}$. then

5: Extract the landmarks at level $l$ as sub-samples, namely
$\Gamma^{(l)}_{m^{(l)}}$.

6: end if

7: Training thread

8: Consider level $l$ and assume that the landmarks $\Gamma^{(l)}_{m^{(l)}}$
are extracted.

9: while $l$ is not sufficiently covered with training points according to
Def. 4.1. do

10: Update
$\big{[}({\mathbf{K}}^{(l)}_{nm})^{\top}{\mathbf{K}}^{(l)}_{nm}\big{]}_{ij}$
and ${\bm{\mathsf{z}}}^{(l)}_{i}$ according to Eq. (2.4) and Eq. (2.5) as new
samples $({\bm{\mathsf{x}}},y)$ arrive, using the landmarks in
$\widetilde{\Gamma}^{(l)}_{m}$ from Def. 3.3.

11: Continuously check if matrices have converged.

12: if Matrices converge according to Def. 4.1 then

13: Update the StreaMRAK regression model $\widetilde{f}^{(L)}$, by including
the correction term $s^{(l)}$ into the Laplacian pyramid, as described in
Section 2.2. Let $L=l$ and update $l=l+1$.

14: end if

15: end while

Algorithm A.3 Update_CoverFraction(node $p$, string s)

1: if s=”No parent found” then

2: Update covering fraction of $p$ with
$\mathfrak{c}\mathfrak{f}(p)=(1-\alpha)\mathfrak{c}\mathfrak{f}(p)$

3: else if s= ”parent found” then

4: Update covering fraction of $p$ with
$\mathfrak{c}\mathfrak{f}(p)=(1-\alpha)\mathfrak{c}\mathfrak{f}(p)+\alpha$

5: end if

## Appendix B. Preparatory material

We offer preparatory material on the damped cover-tree and kernel methods.

### B.1 Preparatory material on the damped cover-tree

This section shows how the recursive formula in Eq. 3.1 approximates the
weighted average of the outcome of the last $N$ random trails. Where the
trails are as described in Section 3.1. By expanding Eq. 3.1 we have
$(\mathfrak{c}\mathfrak{f}(p))_{t}=(1-\alpha)^{t}(\mathfrak{c}\mathfrak{f}(p))_{1}+\alpha\sum_{i=1}^{t-1}(1-\alpha)^{i}\mathbbm{1}_{\mathcal{B}_{c}}({\bm{\mathsf{x}}}_{t-i})$.
Since $(1-\frac{1}{N})^{N}\approx 1/e$, the first term becomes negligible when
$t\gg N$. Similarly, all terms $i>N$ in the sum becomes negligible. This
leaves,

$(\mathfrak{c}\mathfrak{f}(p))_{t}\approx\frac{1}{N}\sum_{i=1}^{N}\bigg{(}1-\frac{1}{N}\bigg{)}^{i}\mathbbm{1}_{\mathcal{B}_{c}}({\bm{\mathsf{x}}}_{t-i})$

which is a weighted average of the outcome of the $N$ last draws as claimed.

######  Remark B.1.

We mention a weakness of the estimator in Eq. (3.1). As follows from Algorithm
A.1, every time a new point ${\bm{\mathsf{x}}}$ is not covered by the existing
children, a new child is added. This consequently updates ${\cal B}_{c}$,
leading to the posterior distribution $\text{Prob}(\mathbbm{1}_{{\cal
B}_{c}}({\bm{\mathsf{x}}})=0|{\bm{\mathsf{x}}})$ to changed every time
$\mathbbm{1}_{{\cal B}_{c}}({\bm{\mathsf{x}}})=0$.

### B.2 Preparatory material on Kernel methods

Kernel methods in the context of reproducing kernel Hilbert spaces (RKHS)
offer a powerful approach to machine learning with a well-established
mathematical foundation [60, 1]. In this paper we consider an input space
${\cal X}\subset\mathbb{R}^{D}$, a corresponding target space ${\cal
Y}\subset\mathbb{R}$ and let $\rho$ be the probability distribution on ${\cal
X}\times{\cal Y}$. Furthermore, we assume an RKHS ${\cal H}_{k}$ generated by
a positive definite kernel $k:{\cal X}\times{\cal X}\rightarrow\mathbb{R}$. In
other words, the eigenvalues $\sigma_{i},\dots,\sigma_{n}$ of the
corresponding kernel matrix
${\mathbf{K}}_{nn}=(k({\bm{\mathsf{x}}}_{i},{\bm{\mathsf{x}}}_{j}))_{i,j=1}^{n}$
satisfies $\sigma_{i}>0$ for all $i\in n$. In this setting the inner product
between two feature vectors
$\phi({\bm{\mathsf{x}}}),\phi({\bm{\mathsf{x}}}^{\prime})\in{\cal H}_{k}$
satisfies the property that
$\left<{\phi({\bm{\mathsf{x}}})},{\phi({\bm{\mathsf{x}}}^{\prime})}\right>_{{\cal
H}_{k}}=k({\bm{\mathsf{x}}},{\bm{\mathsf{x}}}^{\prime})$. This relation, known
as the ”kernel trick” [61, 62], effectively circumvents the need for explicit
construction of non-linear mappings $\phi$.

Given a training set $\\{({\bm{\mathsf{x}}}_{i},y_{i}):i\in[n]\\}$ sampled
according to $\rho$ with $\Gamma_{n}=\\{{\bm{\mathsf{x}}}_{i}:i\in[n]\\}$, we
formulate the kernel ridge regression (KRR) problem as

$\widehat{f}_{n,\lambda}=\operatorname*{argmin}_{f\in\widehat{\cal
H}_{n}}\frac{1}{n}\sum_{i=1}^{n}(f({\bm{\mathsf{x}}}_{i})-y_{i})^{2}+\lambda\left\|{f}\right\|_{\cal
H}^{2},$ (B.1)

where $\lambda>0$ is a regularisation parameter and $\widehat{\cal
H}_{n}=\overline{\textrm{span}}\\{k(\cdot,{\bm{\mathsf{x}}}_{i}):i\in[n]\\}$
is a finite-dimensional subspace of ${\cal H}_{k}$. What is more, for all
$f\in\widehat{\cal H}_{n}$ the Representer theorem [63, 64] guarantees that
there exists coefficients $\alpha_{1},\ldots,\alpha_{n}$ such that the
solution to Eq. (B.1) is on the form

$f({\bm{\mathsf{x}}})=\sum_{i=1}^{n}\alpha_{i}k({\bm{\mathsf{x}}},{\bm{\mathsf{x}}}_{i}).$

Computing the KRR estimator is therefore reduced to solving the linear system

$({\mathbf{K}}_{nn}+\lambda{\mathbf{I}}_{n}){\bm{\upalpha}}={\bm{\mathsf{y}}},$

where ${\bm{\mathsf{y}}}=(y_{1},\ldots,y_{n})^{\top}$,
${\bm{\upalpha}}=(\alpha_{1},\ldots,\alpha_{n})^{\top}$, and
$[{\mathbf{K}}_{nn}]_{ij}=k({\bm{\mathsf{x}}}_{i},{\bm{\mathsf{x}}}_{j})$.

## Appendix C. Proofs and definitions

######  Lemma C.1.

Consider a domain ${\cal X}\in\mathbb{R}^{D}$, a ball ${\cal
B}({\bm{\mathsf{x}}}_{p},r)\subset{\cal X}$ and let
$S=\\{{\bm{\mathsf{x}}}_{i},\,{\bm{\mathsf{x}}}_{j}\in{\cal
B}({\bm{\mathsf{x}}}_{p},r)|\|{\bm{\mathsf{x}}}_{i}-{\bm{\mathsf{x}}}_{j}\|\geq\delta\,\text{for}\,\,i\neq
j\\}$. Furthermore, let the doubling dimension of the set $S$ be
$\texttt{ddim}\coloneqq\texttt{ddim}(S,r)$. We let $c_{d}\coloneqq|S|$ when
$\mathfrak{c}\mathfrak{f}(p)=1$. We then have $2^{\texttt{ddim}-1}\leq
c_{d}\leq 5^{\texttt{ddim}}$.

###### Proof.

The upper bound on $c_{d}$ follows from Lemma C.5 with $r=r_{0}$ and
$\delta=r_{0}/2$. The lower bound follows from the definition of the doubling
dimension 1.2. ∎

######  Lemma C.2.

Let ${\bm{\mathsf{d}}}^{(l)}$ be the residual at level $l$ as defined in Eq.
(2.9). We then have,

${\bm{\mathsf{d}}}^{(l+1)}=({\mathbf{I}}-{\mathbf{K}}_{nn}^{(l)}({\mathbf{K}}_{nn}^{(l)}+\lambda
n{\mathbf{I}})^{-1}){\bm{\mathsf{d}}}^{(l)}$

###### Proof.

Denote ${\bm{\mathsf{s}}}^{(l)}=s^{(l)}([{\bm{\mathsf{x}}}_{n}])$, and note
that
${\bm{\mathsf{s}}}^{(l)}={\mathbf{K}}_{nn}^{(l)}({\mathbf{K}}_{nn}^{(l)}+\lambda
n{\mathbf{I}})^{-1}{\bm{\mathsf{d}}}^{(l)}$. For $l=1$, we have

$\displaystyle{\bm{\mathsf{d}}}^{(1)}$
$\displaystyle={\bm{\mathsf{y}}}-{\bm{\mathsf{s}}}^{(0)}={\bm{\mathsf{y}}}-{\mathbf{K}}_{nn}^{(l)}{\bm{\upalpha}}^{(0)}=({\mathbf{I}}-{\mathbf{K}}_{nn}^{(0)})({\mathbf{K}}_{nn}^{(0)}+\lambda
n{\mathbf{I}})^{-1}{\bm{\mathsf{y}}}.$

We proceed by induction. Assume the statement holds for an $l\geq 2$. We now
have

$\displaystyle{\bm{\mathsf{d}}}^{(l+1)}={\bm{\mathsf{y}}}-\sum_{j=0}^{l}{\bm{\mathsf{s}}}^{(j)}={\bm{\mathsf{d}}}^{(l)}-{\bm{\mathsf{s}}}^{(l)}={\bm{\mathsf{d}}}^{(l)}-{\mathbf{K}}_{nn}^{(l)}({\mathbf{K}}_{nn}^{(l)}+\lambda
n{\mathbf{I}})^{-1}{\bm{\mathsf{d}}}^{(l)}=({\mathbf{I}}-{\mathbf{K}}_{nn}^{(l)}({\mathbf{K}}_{nn}^{(l)}+\lambda
n{\mathbf{I}})^{-1}){\bm{\mathsf{d}}}^{(l)}.$

∎

### C.1 Proof of Thm. 5.4

Let
${\mathbf{P}}_{nn}^{(l)}\coloneqq{\mathbf{K}}^{(l)}_{nn}({\mathbf{K}}^{(l)}_{nn}+\lambda
n{\mathbf{I}})^{-1}$. By definition of ${\bm{\mathsf{d}}}^{(l)}$ and Lemma C.2
it follows

$f([{\bm{\mathsf{x}}}_{n}])-\widehat{f}^{(l+1)}([{\bm{\mathsf{x}}}_{n}])={\bm{\mathsf{d}}}^{(l+1)}=({\mathbf{I}}-{\mathbf{P}}_{nn}^{(l)}){\bm{\mathsf{d}}}^{(l)}=({\mathbf{I}}-{\mathbf{P}}_{nn}^{(l)})(f([{\bm{\mathsf{x}}}_{n}])-\widehat{f}^{(l)}([{\bm{\mathsf{x}}}_{n}])).$
(C.1)

We then have

$\|\widehat{f}^{(l+1)}([{\bm{\mathsf{x}}}_{n}])-f([{\bm{\mathsf{x}}}_{n}])\|\leq\|{\mathbf{I}}-{\mathbf{P}}_{nn}^{(l)}\|\|\widehat{f}^{(l)}([{\bm{\mathsf{x}}}_{n}])-f([{\bm{\mathsf{x}}}_{n}])\|.$
(C.2)

Consider the SVD
${\mathbf{K}}^{(l)}_{nn}={\mathbf{U}}\mathbf{\Sigma}{\mathbf{U}}^{\top}$ where
$\mathbf{\Sigma}=\operatorname{diag}{(\sigma_{l,i})}$ and
$\sigma_{l,n}\leq\dots\leq\sigma_{l,1}$. We then have

$\displaystyle\begin{split}\|{\mathbf{I}}-{\mathbf{P}}^{(l)}_{nn}\|&=\bigg{\|}{\mathbf{U}}\operatorname{diag}{\Big{(}\frac{n\lambda}{n\lambda+\sigma_{l,i}}\Big{)}}{\mathbf{U}}^{\top}\bigg{\|}=\bigg{\|}\operatorname{diag}{\Big{(}\frac{n\lambda}{n\lambda+\sigma_{l,i}}\Big{)}}\bigg{\|}\\\
&=\frac{n\lambda}{n\lambda+\sigma_{l,n}}\coloneqq
1-\varepsilon(l),\end{split}$ (C.3)

and Thm. 5.4 follows recursively from Eq. (C.2) and Eq. (C.3). ∎

### C.2 Proof of Thm. 5.5

To bound the smallest eigenvalue of the kernel matrix
$[{\mathbf{K}}^{(l)}_{nn}]_{ij}=\Phi(\|{\bm{\mathsf{x}}}_{i}-{\bm{\mathsf{x}}}_{j}\|)$,
namely $\sigma_{l,n}$, we will assume that there exists a lower bound on the
minimal distance between any two points
${\bm{\mathsf{x}}}_{i},{\bm{\mathsf{x}}}_{j}\in{\cal X}$, namely
$\delta\coloneqq\min\limits_{i\neq j\in{\cal
X}}\left\|{{\bm{\mathsf{x}}}_{i}-{\bm{\mathsf{x}}}_{j}}\right\|>0$. Consider
the Gaussian
$\Phi({\bm{\mathsf{x}}})=\exp(-\beta\|{\bm{\mathsf{x}}}\|_{2}^{2})$,
$\beta>0$, with the Fourier transform
$\widehat{\Phi}(\omega)=(\pi/\beta)^{D/2}\exp(-\|\omega\|_{2}^{2}/4\beta)$.
From [65, Corollary 12.4] we have the bound

$\sigma_{l,n}\geq
C_{D}2^{D}(2\beta)^{-D/2}\delta^{-D}\exp(-4M_{D}^{2}/(\delta^{2}\beta)),$

where

$M_{D}=12\bigg{(}\frac{\pi\Gamma^{2}(D/2+1)}{9}\bigg{)}^{1/(D+1)}\quad\text{and}\quad
C_{D}=\frac{1}{2\Gamma(D/2+1)}\bigg{(}\frac{M_{D}}{2^{3/2}}\bigg{)}^{D}.$

With $\beta=(\sqrt{2}2^{-l}r_{0})^{-2}$ we then have

$\displaystyle\begin{split}\sigma_{l,n}&\geq
C_{D}2^{D}2^{-Dl}\bigg{(}\frac{r_{0}}{\delta}\bigg{)}^{D}\exp\big{(}-(2\sqrt{2}M_{D})^{2}(r_{0}/\delta)^{2}4^{-l}\big{)}\\\
&=C_{1,D}2^{-Dl}\exp\big{(}-C_{2,D}4^{-l}\big{)}\coloneqq B(l),\end{split}$

where we define

$C_{1,D}=\frac{1}{2}(6\sqrt{2})^{D}\Gamma(D/2+1)^{\frac{D-1}{D+1}}\bigg{(}\frac{\pi}{9}\bigg{)}^{\frac{D}{D+1}}\bigg{(}\frac{r_{0}}{\delta}\bigg{)}^{D}\quad\text{and}\quad
C_{2,D}=1152\bigg{(}\frac{\pi\Gamma^{2}(D/2+1)}{9}\bigg{)}^{\frac{2}{D+1}}\bigg{(}\frac{r_{0}}{\delta}\bigg{)}^{2}.$

The first bound in Thm. 5.5 follows from this result. ∎

######  Remark C.3.

In [65, Thm. 12.3] they also offer an a fortiori bound corresponding to
$M_{D}=6.38D$,
$C_{1,D}=\frac{1}{2}\big{(}\frac{12.76}{2^{3/2}}\big{)}^{D}\big{(}\frac{D^{D}}{\Gamma(D/2+1)}\big{)}\big{(}\frac{r_{0}}{\delta}\big{)}^{D}$
and $C_{2,D}=(12.76\sqrt{2}D)^{2}(r_{0}/\delta)^{2}$.

######  Corollary C.4.

We note that $B(l)$ has a maximum at

$\displaystyle\begin{split}l^{*}&=\frac{1}{2}\log_{2}\bigg{(}\frac{C_{2,D}\log
4}{D\log
2}\bigg{)}=\log_{2}\bigg{(}\sqrt{\frac{D}{2}}\bigg{(}\frac{r_{0}}{\delta}\bigg{)}\bigg{)}+\log_{2}\bigg{(}\frac{4M_{D}}{D}\sqrt{2}\bigg{)}\end{split}$

and is monotonically increasing with $l$ on the interval $l\in(0,l^{*})$.
Furthermore, with the a fortiori expression for $M_{D}$ from Remark C.3 we
have

$l^{*}=\log_{2}\bigg{(}\sqrt{\frac{D}{2}}\bigg{(}\frac{r_{0}}{\delta}\bigg{)}\bigg{)}+\log_{2}\bigg{(}25.52\sqrt{2}\bigg{)}.$

When the level $l$ becomes sufficiently large, the kernel matrix
${\mathbf{K}}^{(l)}_{nn}$ becomes diagonally dominant, and we can therefore
bound the eigenvalues using Garschgorins Theorem [66, Thm. 1.1], which gives

$|\sigma_{l,i}-[{\mathbf{K}}^{(l)}_{nn}]_{jj}|=|\sigma_{l,i}-1|<\sum_{\begin{subarray}{c}q=1,\\\
q\neq j\end{subarray}}^{n}|[{\mathbf{K}}^{(l)}_{nn}]_{jq}|\quad\text{for}\quad
i,j\in[n].$ (C.4)

To find a more explicit bound, we analyze the sum on the right-hand side using
Lemma C.5.

######  Lemma C.5.

Consider a ball ${\cal B}({\bm{\mathsf{x}}},r)\in\mathbb{R}^{D}$ and let
$\delta>0$. The number of points in any (discrete) set of points within ${\cal
B}({\bm{\mathsf{x}}},r)$ that are at least $\delta$ apart,
$S=\\{{\bm{\mathsf{x}}}_{i}\in{\cal
B}({\bm{\mathsf{x}}},r)|d({\bm{\mathsf{x}}}_{i},{\bm{\mathsf{x}}}_{j})\geq\delta\,\text{for}\,\,i\neq
j\\}$, is bounded by $|S|\leq\bigg{(}\frac{2r}{\delta}+1\bigg{)}^{D}$.

###### Proof.

Since the points in $S$ are at least $\delta$ apart, it follows that the balls
${\cal B}({\bm{\mathsf{x}}}_{i},\delta/2)$ are disjoint. Consider now the ball
${\cal B}({\bm{\mathsf{x}}},r+\delta/2)$. All of the balls ${\cal
B}({\bm{\mathsf{x}}}_{i},\delta/2)$ are entirely contained within ${\cal
B}({\bm{\mathsf{x}}},r+\delta/2)$. Since the balls ${\cal
B}({\bm{\mathsf{x}}}_{i},\delta/2)$ are disjoint, it follows that

$|S|\leq\frac{\operatorname*{Vol}\Big{(}{\cal
B}({\bm{\mathsf{x}}},r+\delta/2)\Big{)}}{\operatorname*{Vol}\Big{(}{\cal
B}({\bm{\mathsf{x}}}_{i},\delta/2)\Big{)}}=\bigg{(}\frac{2r}{\delta}+1\bigg{)}^{D}.$

∎

Consider a family of annuli $\\{R_{t}\\}_{t=0}^{\infty}$ where $R_{t}={\cal
B}({\bm{\mathsf{x}}}_{j},2^{t+1}\delta)\backslash{\cal
B}({\bm{\mathsf{x}}}_{j},2^{t}\delta)$. Inspired by [30], we can interpret the
right hand side of Eq. (C.4) as a sum over $\\{R_{t}\\}_{t=0}^{\infty}$. The
entries of ${\mathbf{K}}^{(l)}_{nn}$ are defined as

$[{\mathbf{K}}^{(l)}_{nn}]_{ij}=\exp{\bigg{(}-\frac{\left\|{{\bm{\mathsf{x}}}_{i}-{\bm{\mathsf{x}}}_{j}}\right\|^{2}}{2r_{l}^{2}}\bigg{)}},\quad\forall
i,j\in[n],$

where $r_{l}=2^{-l}r_{0}$ for $r_{0}>0$. It follows

$\displaystyle\begin{split}\sum_{\begin{subarray}{c}q=1,\\\ q\neq
j\end{subarray}}^{n}|[{\mathbf{K}}^{(l)}_{nn}]_{jq}|=\sum_{t=0}^{\infty}\sum_{{\bm{\mathsf{x}}}_{q}\in
R_{t}}k^{(l)}({\bm{\mathsf{x}}}_{j},{\bm{\mathsf{x}}}_{q})\leq\sum_{t=0}^{\infty}\bigg{(}\frac{2^{t+2}\delta}{\delta}+1\bigg{)}^{D}\exp\big{(}-(2^{t}\delta
2^{-1/2}r_{l}^{-1})^{2}\big{)}\end{split}$

where in the first term on the right-hand side we bound the number of summands
using Lemma C.5, and in the second we use
$\|{\bm{\mathsf{x}}}_{q}-{\bm{\mathsf{x}}}_{j}\|\geq 2^{t}\delta$ for
${\bm{\mathsf{x}}}_{q}\in R_{t}$. Note now that for all $T\geq 1$ there exists
$C_{T}>0$ such that $\exp(-r^{2})\leq C_{T}r^{-T}$ holds for all $r>0$. Such a
constant is given by the Lambert W function and satisfies
$C_{T}=\left(\frac{T}{2{\rm e}}\right)^{T/2}$. Moreover, $2^{t+2}+1\leq
2^{t+2+\alpha}$, for $\alpha\geq\ln(1+1/4)/\ln(2)$. Thus,

$\displaystyle\sum_{\begin{subarray}{c}q=1,\\\ q\neq
j\end{subarray}}^{n}|[{\mathbf{K}}^{(l)}_{nn}]_{jq}|\leq
C_{T}\left(\frac{r_{l}}{\delta}\right)^{T}2^{(2+\alpha)D+T/2}\sum_{t=0}^{\infty}2^{t(D-T)}\leq
C_{DT}\left(\frac{r_{l}}{\delta}\right)^{T},$

where using $\sum_{t=0}^{\infty}2^{t(D-T)}\leq 2$, which holds for $D-T<0$, we
let

$2\cdot 2^{D(2+\alpha)+T/2}C_{T}\leq 2\cdot
2^{D(2+\alpha)-T/2(1+1/\ln(2))}T^{T/2}=:C_{DT},$

where we used $\exp(1)\geq 2^{1+1/\ln(2)}$.

We now consider the function

$\displaystyle\begin{split}F(T)&\coloneqq
2^{-\frac{T}{2}(1+1/\ln{2})}2^{-lT}T^{\frac{T}{2}}\bigg{(}\frac{r_{0}}{\delta}\bigg{)}^{T}\\\
&=2^{-\frac{T}{2}(1+1/\ln{2})}2^{-lT}2^{T/2\log_{2}T}2^{T\log_{2}(r_{0}/\delta)}\\\
&=2^{-T/2(B-\log_{2}T)}=2^{f(T)},\end{split}$

Where $B=1+\frac{1}{\ln 2}+2l-2\log_{2}\bigg{(}\frac{r_{0}}{\delta}\bigg{)}$.
$F$ is minimized by

$T^{*}=2^{B-1/\ln 2}=2^{1+1/\ln 2+2l-2\log_{2}(r_{0}/\delta)-1/\ln 2}=2\cdot
4^{l-\log_{2}(r_{0}/\delta)},$

such that

$F(T^{*})=2^{-\frac{T}{2}(B-\log_{2}2^{B}-\log_{2}2^{-1/\ln
2})}=2^{\frac{-T^{*}}{2\ln 2}}=2^{-\frac{1}{\ln
2}4^{l-\log_{2}(r_{0}/\delta)}}.$

Inserting this back and with $\alpha=\ln(1+1/4)/\ln 2$, we have

$\sigma_{l,n}>1-\sum_{\begin{subarray}{c}q=1,\\\ q\neq
j\end{subarray}}^{n}|[{\mathbf{K}}^{(l)}_{nn}]_{jq}|\geq
1-2^{1+\frac{1}{\ln{2}}((\ln{(1+1/4)}+2\ln{2})D-g(l))},\quad
g(l)=4^{l-\log_{2}{r_{0}/\delta}}$

With $C_{3}=(\ln{(1+1/4)}+2\ln{2})$ this leads to

$0<1-\varepsilon(l)<\bigg{(}1+\big{(}1-2^{1+\frac{1}{\ln{2}}(C_{3}D-g(l))}\big{)}/n\lambda\bigg{)}^{-1},$
(C.5)

We note that the bound in Eq. (C.5) holds for $T^{*}>D$ which means that
$l>\log_{2}(\sqrt{D/2}r_{0}/\delta)$. ∎

### C.3 Proof of Corollary 5.6

Follows from Eq. (C.1)-(C.3) with
$P^{(l)}_{nn}=P^{(l)}_{nn,k}+\big{(}P^{(l)}_{nn,k}\big{)}^{\perp}$, where
$P^{(l)}_{nn,k}$ is the projection on the eigenvectors associated with the $k$
largest eigenvalues and
$\big{(}P^{(l)}_{nn,k}\big{)}^{\perp}(\widehat{f}^{(l)}([{\bm{\mathsf{x}}}_{n}])-f([{\bm{\mathsf{x}}}_{n}]))=0$.
∎