# Sketching the Heat Kernel: Using Gaussian Processes to Embed Data

Anna Gilbert and Kevin O’Neill

###### Abstract

This paper introduces a novel, non-deterministic method for embedding data in
low-dimensional Euclidean space based on computing realizations of a Gaussian
process depending on the geometry of the data. This type of embedding first
appeared in [1] as a theoretical model for a generic manifold in high
dimensions.

In particular, we take the covariance function of the Gaussian process to be
the heat kernel, and computing the embedding amounts to sketching a matrix
representing the heat kernel. The Karhunen-Loève expansion reveals that the
straight-line distances in the embedding approximate the diffusion distance in
a probabilistic sense, avoiding the need for sharp cutoffs and maintaining
some of the smaller-scale structure.

Our method demonstrates further advantage in its robustness to outliers. We
justify the approach with both theory and experiments.

##  1 Introduction

Recent success in the analysis of high-dimensional data is often attributed,
at least in part, to the tendency of such data to have some sort of underlying
low-dimensional structure. For example, a collection of photographs of the
same object but taken at different angles will compose a small subset of the
space of all digital images. To take advantage of this structure, it is often
desirable to embed the original data in a relatively low-dimensional Euclidean
space, allowing one to better understand this structure or more quickly run
algorithms to analyze it (such as those for clustering).

In this paper, we introduce a novel method for embedding data in low-
dimensional Euclidean spaces. While motivated by theoretical results on
manifolds (the oft-assumed underlying structure of high-dimensional data), the
method ultimately relies on a heat kernel and a notion of distance, suggesting
it may bear fruit when applied to any of the large class of metric spaces on
which a well-behaved heat kernel exists, such as weighted graphs (see [20] and
the related literature, for instance).

The method works by constructing a Gaussian process $f$ on the data. Let
$f_{1},\ldots,f_{k}$ denote independent realizations of $f$. Then, the
embedding of the data into $\mathbb{R}^{k}$ is computed via the formula

$h^{k}(x)=\frac{1}{\sqrt{k}}(f_{1}(x),\ldots,f_{k}(x)).$ (1)

We call an embedding of the form (1) a _Gaussian process embedding_.

Gaussian process embeddings seem to have first appeared in [1]. There, the
authors show that the reach of the normalized Gaussian process embedding of a
fixed manifold $M$ converges almost surely to a constant depending only on $M$
and the choice of Gaussian process as $k\to\infty$. The result may be
interpreted as saying that for most manifolds embedded in high dimensions, the
extrinsic geometry will not be an obstacle for analysis. The purpose of the
Gaussian process embedding is to determine a notion of “most manifolds” in a
way which can be analyzed using previous work on Gaussian processes on
manifolds.

Our motivation lies in the followup work of [11], which established that the
metric induced by such Gaussian process embeddings converges in $C^{2}$ to the
original metric on $M$ almost surely as $k\to\infty$. For practical
applications, we work with finite subsets of $\mathbb{R}^{n}$ in place of
continuous manifolds, so our focus will be on the straight-line distance of
the embeddings: $\|h^{k}(x)-h^{k}(y)\|_{2}$. We show in Propositions 2.4 and
3.2 that this distance in the embedding approximates a diffusion distance on
the original data.

We provide a brief summary of the diffusion maps embedding and related
diffusion distance [12, 7]. Given a smooth, decaying kernel on a subset $X$ of
$\mathbb{R}^{n}$, one may normalize large powers of this kernel to approximate
the heat kernel on $X$. In particular, this process allows one to construct a
heat kernel on a finite subset beginning with only an affinity matrix. The
eigenfunctions of the heat kernel are then used to embed the data into a low-
dimensional Euclidean space, where the straight-line distance between points
is an approximation of the diffusion distance. The diffusion distance, while
perhaps different from the original distance, is itself a valuable quantity,
as it better reflects the connectivity of the data as a graph. Furthermore,
the metric induced by the diffusion distance on a Riemannian manifold agrees
with the original metric, but the diffusion distance is more robust than the
geodesic distance when working with finite sets. See Subsections 2.1 and 4.1
for more.

While embedding with diffusion maps remains a useful tool, Gaussian process
embeddings, in contrast, do not rely on a cutoff of which top eigenvalues to
consider. Rather, as shown by the Karhunen-Loève expansion (Theorem 2.3), they
approximate the diffusion distance by combining all the eigenvectors/functions
in each component, allowing one to preserve small-scale information which is
neglected by diffusion maps. See Subsection 2.2 for the example of a stretched
flat torus. Experimentally, we find that Gaussian process embeddings also
perform well with regards to robustness to outliers. They are also at least as
easy to compute as diffusion maps; constructing a Gaussian process embedding
according to a Gaussian process with covariance matrix equal to the heat
kernel at time $2t$ amounts to sketching the heat kernel at time $t$.
Specifics of computation are discussed more in Subsection 4.5.

We now take a moment to discuss two issues of embeddability in Euclidean space
and their relevance with regards to our method. The Whitney embedding theorem
(as typically first proven in a graduate course) states that any compact
$d$-dimensional manifold may be embedded into $\mathbb{R}^{2d+1}$ (see [16],
for instance). As the authors of [1] remark, 1 is almost surely an embedding
when $k\geq 2d+1$. The Nash embedding theorem requires $k$ on the order of
$d^{2}$ for an embedding into $\mathbb{R}^{k}$ to be isometric; the main
result of [11] shows this occurs for Gaussian process embeddings as
$k\to\infty$.

The paper of Jones, Maggioni, and Schul [9] shows that small neighborhoods of
a $d$-dimensional compact Riemannian manifold may be parameterized by $d$
cleverly chosen eigenfunctions of the Laplace-Beltrami operator on the entire
manifold. An article by Bates [3] proves that the entire manifold may be
embedded in $\mathbb{R}^{N}$ via the lowest $N$ eigenfunctions, where $N$
depends on $d$, a lower bound for injectivity radius, a lower bound for Ricci
curvature, and a volume bound. By varying such quantities, one may require
arbitrarily large $N$ to get an embedding (isometric or not).

An intermediate result in [9] is that a compact, $d$-dimensional Riemannian
manifold may be locally parameterized by its heat kernel evaluated at $d$
distinct and well-chosen points. Computing such a parameterization would
require taking a subsampling of the columns of the heat kernel matrix. Viewing
subsampling as a form of sketching offers parallels with our method.

The second issue regarding embedding involves data which cannot be embedded in
a bi-Lipschitz manner in Euclidean space of any dimension, such as the
Heisenberg group [19]. This may be considered another reason to use the
diffusion distance, which still exists on non-Riemannian manifolds and yet may
be realized in an embedding in Euclidean space.

In Section 2, we formalize the theory of Gaussian process embeddings and use
Dudley’s inequality to determine a rate of convergence of the induced
distances to the diffusion distance. We discuss in detail the Karhunen-Loève
expansion for Gaussian fields and use it to analyze the differences between
Gaussian process embeddings and diffusion maps, which are also formalized in
this section.

Section 3 focuses on the discretization of the tools and results developed in
Section 2, taking us one step closer to constructing Gaussian process
embeddings of finite datasets.

Section 4 states explicitly the algorithms used to compute Gaussian process
embeddings, their runtimes, and how they relate to diffusion maps. We also
develop a version of the algorithm in which sketching by a Gaussian matrix is
replaced with sketching by a matrix of i.i.d. symmetric Bernoulli variables.

Lastly, Section 5 includes a wide variety of experiments addressing
comparisons with diffusion maps, with a particular emphasis on robustness to
outliers and the dimension of Euclidean space required for an embedding.

##  2 Continuous Theory of Embeddings

###  2.1 Diffusion Maps

Let $(M,g_{M})$ be a compact, Riemannian manifold without boundary whose
Laplace-Beltrami operator has eigenvalues
$0=\lambda_{0}\leq\lambda_{1}\leq\lambda_{2}\leq\ldots\nearrow\infty$ and a
corresponding orthonormal basis of eigenfunctions
$\\{\varphi_{i}\\}_{i=0}^{\infty}$. The heat kernel on $M$ at time $t$ may
then be written as

$k_{t}(x,y)=\sum_{i=0}^{\infty}e^{-\lambda_{i}t}\varphi_{i}(x)\varphi_{i}(y).$

Thus, viewed as a linear transformation on $L^{2}(M)$ defined by
$f(x)\mapsto\int_{M}f(x)k_{t}(x,y)dx$, the heat kernel has an orthonormal
basis of eigenfunctions $(e^{-\lambda_{i}t}\varphi_{i}(x))_{i=0}^{\infty}$
with corresponding eigenvalues $(e^{-\lambda_{i}t})_{i=0}^{\infty}$.

Let $\ell^{2}$ denote the Hilbert space of real-valued, square-summable
sequences $(a_{i})_{i=1}^{\infty}$ with metric
$d_{\ell^{2}}\left((a_{i})_{i=1}^{\infty},(b_{i})_{i=1}^{\infty}\right)=\sum_{i=1}^{\infty}(a_{i}-b_{i})^{2}$.
The following result of Berard, Besson, and Gallot may be viewed as motivation
for diffusion maps (analogous to our inspiration in [11]).

######  Theorem 2.1.

[5] Let $(M,g_{M})$, $(\varphi_{i})_{i=0}^{\infty}$, and
$(\lambda_{i})_{i=0}^{\infty}$ be as above. Then,

$\psi_{t}:x\mapsto\sqrt{2}(4\pi)^{n/2}t^{(n+2)/4}\left(e^{-\lambda_{j}t/2}\varphi_{j}(x)\right)_{j=1}^{\infty}$
(2)

is an embedding of $M$ into $\ell^{2}$ for all $t>0$. Furthermore, the
pullback metric $\psi_{t}^{*}d_{\ell^{2}}$ is asymptotic to $g_{M}$ as $t\to
0$.

Observe that the top eigenfunction $\varphi_{0}$ is dropped from (2) as it is
trivial (equal to a nonzero constant function).

One may desire to simply define the diffusion maps embedding to be

$x\mapsto(e^{-\lambda_{1}t}\varphi_{1}(x),\ldots,e^{-\lambda_{k}t}\varphi_{k}(x)),$
(3)

where the constant scaling factor from (2) is dropped and $t/2$ is replaced
with $t$ in the exponent as this is a mere re-scaling of the time variable.
However, in practice one may not have direct access to the heat kernel on a
given set (especially if it is a finite sample from a manifold), so it is
beneficial to derive a version for more arbitrary kernels.

Let $(X,\mu)$ be a measure space with a symmetric, positive semidefinite
kernel $a:X\times X\to[0,\infty)$. The kernel $a$ will be chosen to be similar
to the heat kernel, but for now we maintain generality as we will make two
distinct choices for it.

If

$\int_{X}\int_{X}a(x,y)^{2}d\mu(x)d\mu(y)<\infty,$

then the associated operator

$f\mapsto\int_{X}a(x,y)f(x)dx$

is compact and has the eigendecomposition

$a(x,y)=\sum_{l\geq 0}\lambda_{l}\phi_{l}(x)\phi_{l}(y),$

where $\lambda_{0}\geq\lambda_{1}\geq\ldots\geq 0$. From this, we define

$a^{(t)}(x,y)=\sum_{l\geq 0}\lambda_{l}^{t}\phi_{l}(x)\phi_{l}(y).$

We define the diffusion distance at time $t$, $D^{(t)}:X\times
X\to\mathbb{R}$, by

$D_{a}^{(t)}(x,y)^{2}=\int_{X}(a^{(t)}(x,u)-a^{(t)}(y,u))^{2}d\mu(u)=\sum_{l\geq
0}\lambda_{l}^{2t}(\phi_{l}(x)-\phi_{l}(y))^{2}.$

The diffusion distance is equal to the straight line ($\ell^{2}$) distance
from the embedding

$x\mapsto(\lambda_{l}^{t}\phi_{l})_{l\geq 0}.$ (4)

Thus, we define the diffusion maps embedding by

$x\mapsto(\lambda_{0}^{t}\phi_{0}(x),\ldots,\lambda_{k}^{t}\phi_{k}(x)).$ (5)

By truncating the embedding to lie in $\mathbb{R}^{k}$, we compute the
following approximation of the diffusion distance

$\tilde{D}^{(t)}(x,y)^{2}=\sum_{l=0}^{k}\lambda_{l}^{t}(\phi_{l}(x)-\phi_{l}(y))^{2}.$

Observe that if $a^{(t)}(x,y)$ is the heat kernel at time $t$, then (4) is the
embedding from (2). We discuss how one might obtain the kernel $a$ in a
discretized setting in Subsection 4.1. One often desires to work with a kernel
which is not necessarily symmetric but row-stochastic; however, we will focus
on the symmetric kernels for better comparisons with Gaussian process
embeddings. For full details, see [7, 17, 12].

###  2.2 Challenges of Eigenfunction Embeddings

It is a slight misnomer to refer to (3) as an “embedding” since for any given
$d$-dimensional manifold $M$ it is possible to choose $k$ small enough such
that it is not even injective. For instance, one may take $k$ to be smaller
than the (intrinsic) dimension of $M$ whenever the latter is greater than 1.
However, one may always choose $k$ large enough to obtain an embedding.

That said, for the diffusion maps embedding one may not determine the
appropriate $k$ merely as a function of $d$, as demonstrated by the following
example. Let $M=S^{1}\times rS^{1}$, where $S^{1}$ is the unit circle,
identified with the interval $[0,2\pi)$, and $rS^{1}$ is the circle of radius
$r>1$, identified with the interval $[0,2r\pi)$. A well-known orthonormal
basis of eigenfunctions of the Laplace-Beltrami operator $\Delta_{M}$ of $M$
consists of functions of the form

$\varphi^{cc}_{m_{1},m_{2}}(x,y)=\cos(m_{1}x)\cos(m_{2}y/r),\varphi^{cs}_{m_{1},m_{2}}(x,y)=\cos(m_{1}x)\sin(m_{2}y/r),\\\
\varphi^{sc}_{m_{1},m_{2}}(x,y)=\sin(m_{1}x)\cos(m_{2}y/r),\varphi^{ss}_{m_{1},m_{2}}(x,y)=\sin(m_{1}x)\sin(m_{2}y/r).$
(6)

Here, $m_{1},m_{2}$ are arbitrary nonnegative integers which are allowed to be
zero in the case they lie inside a cosine expression, resulting in a nonzero
constant function. The corresponding eigenvalues are
$m_{1}^{2}+m_{2}^{2}/r^{2}$; thus, the lowest nonzero eigenvalues are (for
non-integer $r$)

$\frac{1}{r^{2}},\frac{1}{r^{2}},\ldots,\frac{\lfloor
r\rfloor^{2}}{r^{2}},\frac{\lfloor r\rfloor^{2}}{r^{2}}.$

The eigenfunctions of $\Delta_{M}$ corresponding to these eigenvalues, and
likewise the top nontrivial eigenfunctions of the heat kernel, are of the form

$\varphi^{cc}_{0,m_{2}}(x,y)=\cos(m_{2}y/r),\varphi^{cs}_{0,m_{2}}(x,y)=\sin(m_{2}y/r),\hskip
18.06749pt1\leq m_{2}\leq\lfloor r\rfloor$

and don’t “see” the $x$-coordinate. As a result, the diffusion maps embedding
in (5) is not truly an embedding for $k\leq 2\lfloor r\rfloor$, and this issue
may be exacerbated by taking $r$ arbitrarily large.

The commonly used spectral embedding [4]

$x\mapsto(\varphi_{1}(x),\ldots,\varphi_{k}(x))$

faces the same obstacle. In contrast, and as observed in [1], Gaussian process
embeddings are almost surely embeddings whenever $k\geq 2d+1$

###  2.3 Gaussian Processes and Riemannian Manifolds

Given a parameter set $T$, a _Gaussian process_ is a collection of Gaussian
random variables $f:=(f_{t})_{t\in T}$. We use the notation $f_{t}$ and $f(t)$
interchangeably. A Gaussian process is _zero-mean_ if $\mathbb{E}f_{t}=0$ for
all $t\in T$. A zero-mean Gaussian process is determined uniquely by its
positive semi-definite covariance function

$C(s,t):=\mathbb{E}[f(s)f(t)].$

Each Gaussian process as above is associated with a _canonical metric_

$D_{C}(s,t):=\left(\mathbb{E}(f_{s}-f_{t})^{2}\right)^{1/2}.$

While actually a pseudometric in the sense that $D_{C}(s,t)=0$ need not imply
$s=t$ (e.g., if $C\equiv 0$), we will keep with the standard terminology.

A _random process_ is a collection of (not necessarily Gaussian) random
variables $(f_{t})_{t\in T}$; we may also refer to its covariance function
(defined the same as above), though two different random processes may have
the same covariance function.

Of great motivation to our work is the following result from [11]. Let
$(M,g_{M})$ be a Riemannian manifold with a Gaussian process $f$ such that

$g_{M}(X,Y)=[Y_{y}X_{x}C(x,y)]|_{x=y},$ (7)

where $X$ and $Y$ are vector fields with values $X_{x},Y_{x}\in T_{x}M$, the
tangent space of $M$. If (7) holds, then we say $g_{M}$ is the _metric induced
by_ $f$. Given a Riemannian manifold $(M,g_{M})$, there always exists a
Gaussian process $f$ inducing $g_{M}$; this is a consequence of the Nash
embedding theorem (as in the proof of Theorem 12.6.1 in [2]). In fact, if
$\iota:M\to\mathbb{R}^{k}$ is an isometric embedding, then given any isotropic
kernel $k(x,y)$ (that is, depending only on $|x-y|$) one may choose $f$ by
setting $C(x,y)=k(\iota(x),\iota(y))$.

Define a (random) map:

$h^{k}(x)=\frac{1}{\sqrt{k}}(f_{1}(x),\ldots,f_{k}(x)),$ (8)

where $f_{1},\ldots,f_{k}$ are independent realizations of $f$. Let
$g_{E}^{k}$ denote the metric induced on $h^{k}(M)$ by the Euclidean metric in
$\mathbb{R}^{k}$.

######  Theorem 2.2.

[11] Let $(M,g_{M})$ be a connected, orientable, compact, $C^{3}$ Riemannian
manifold and $f:M\to\mathbb{R}$ a zero-mean Gaussian process with a.s. $C^{3}$
sample paths inducing the metric $g_{M}$. Then, with probability one, the
pullback metric $(h^{k})^{*}g_{E}^{k}$ converges to $g_{M}$ in the $C^{2}$
topology.

The importance of $C^{2}$-convergence is that it implies the convergence of
intrinsic functionals on $M$, such as volume and Lipschitz-Killing curvatures
whose computations depend on derivatives of $g_{M}$. However, just the
assumption that $M$ be $C^{1}$ would be enough to establish convergence of the
pullback metric. This topic is discussed further in [1].

While motivated by Theorem 2.2, we will not apply it directly or even
discretize it too directly. In particular, we will pick the Gaussian process
with covariance function equal to the heat kernel in order to obtain the
diffusion distance, even if the heat kernel is not isotropic.

###  2.4 Karhunen-Loève Expansion

The following theorem is a fundamental result in the study of Gaussian
processes which can be found in various forms throughout the literature. (See
[14] or Chapter 3 of [2], for instance.)

######  Theorem 2.3 (Karhunen-Loève Expansion).

Let $T\subset\mathbb{R}^{n}$ be compact and $C(x,y)$ be a symmetric, positive
semidefinite function on $T\times T$ such that the integral operator
$g\mapsto\int_{T}C(x,y)g(y)dy$ is compact on $L^{2}(T)$. Let
$\lambda_{1}\geq\lambda_{2}\geq\ldots$ be the eigenvalues of the integral
operator with associated orthonormal eigenfunctions
$\varphi_{1},\varphi_{2},\ldots$.

Then, the Gaussian process $f=(f_{t})_{t\in T}$ with covariance function
$C(x,y)$ is of the form

$f(x)=\sum_{i=1}^{\infty}\xi_{i}\sqrt{\lambda_{i}}\varphi_{i}(x),$ (9)

where $\xi_{1},\xi_{2},\ldots$ are i.i.d. $N(0,1)$ random variables.

The convergence in (9) is in the sense of $L^{2}(T)$ almost surely; however,
one may show the convergence is uniform almost surely if $C(s,t)$ is
Lipschitz. This will be enough to justify our operations with the expression.
(We chose the Lipschitz condition for simplicity’s sake; weaker conditions
would also suffice.)

When dealing with finite data, one may not have a desirable notion of a
geodesic distance or a Riemannian metric at a point; even if the data is
assumed to have been sampled from a manifold, approximation of these
quantities may be considered difficult. However, we may easily derive
estimates on the straight-line distance and an approximation of the Riemannian
metric over distances from the Karhunen-Loève expansion. This approximation
works by replacing $\langle X,Y\rangle$ ($X,Y$ lying in the tangent space of
the manifold at a point $x$) with the Euclidean inner product of secant lines
intersecting at $x$.

######  Proposition 2.4.

Let $T,f,C,\lambda_{i}$, and $\varphi_{i}$ be as in Theorem 2.3. Fix $x,y,z\in
T$ and $k\geq 0$. Then,

$\mathbb{E}\langle
h^{k}(y)-h^{k}(x),h^{k}(z)-h^{k}(x)\rangle=\sum_{i=1}^{\infty}\lambda_{i}(\varphi_{i}(y)-\varphi_{i}(x))(\varphi_{i}(z)-\varphi_{i}(x)).$
(10)

In particular,

$D_{C}(x,y)^{2}=\mathbb{E}(h^{1}(x)-h^{1}(y))^{2}=\mathbb{E}(h^{k}(x)-h^{k}(y))^{2}=\sum_{i=1}^{\infty}\lambda_{i}(\varphi_{i}(x)-\varphi_{i}(y))^{2}=D_{C}^{(1/2)}(x,y)^{2}.$
(11)

More generally, $D_{A}^{(t)}=D_{A^{2t}}$.

###### Proof.

Observe that

$\mathbb{E}\langle
h^{k}(y)-h^{k}(x),h^{k}(z)-h^{k}(x)\rangle=\frac{1}{k}\mathbb{E}\sum_{i=1}^{k}(f_{i}(y)-f_{i}(x))(f_{i}(z)-f_{i}(x)).$

Thus, it suffices to show that

$\mathbb{E}(f(y)-f(x))(f(z)-f(x))=\sum_{i=1}^{\infty}\lambda_{i}(\varphi_{i}(x)-\varphi_{i}(y))(\varphi_{i}(x)-\varphi_{i}(z)).$

By (9) and the independence of the $\xi_{i}$,

$\displaystyle\mathbb{E}(f(y)-f(x))(f(z)-f(x))$
$\displaystyle=\mathbb{E}\left[\left(\sum_{i=1}^{\infty}\xi_{i}\sqrt{\lambda_{i}}(\varphi_{i}(y)-\varphi_{i}(x))\right)\left(\sum_{i=1}^{\infty}\xi_{i}\sqrt{\lambda_{i}}(\varphi_{i}(z)-\varphi_{i}(x))\right)\right]$
$\displaystyle=\mathbb{E}\left[\sum_{i,j}\xi_{i}\xi_{j}\sqrt{\lambda_{i}\lambda_{j}}(\varphi_{i}(y)-\varphi_{i}(x))(\varphi_{j}(z)-\varphi_{j}(x))\right]$
$\displaystyle=\mathbb{E}\sum_{i=1}^{\infty}\xi_{i}^{2}\lambda_{i}(\varphi_{i}(y)-\varphi_{i}(x))(\varphi_{i}(z)-\varphi_{i}(x)),$

which is equal to the right-hand-side of (10) by the fact that
$\mathbb{E}\xi_{i}^{2}=1$ for $\xi_{i}\sim N(0,1)$.

The second equality in (11) follows immediately as the special case of (10)
where $y=z$. The first and fourth follow trivially from the definitions of the
distances.

∎

###  2.5 Quantitative Convergence

Given a pseudometric space $(T,d)$, let $\operatorname{diam}(T,d)$ denote its
diameter, that is, $\sup_{s,t\in T}d(s,t)$. If $\epsilon>0$, then
$\mathcal{N}(T,d;\epsilon)$ is the covering number of $T$, that is, the
smallest number $N$ such that $T$ is contained in the union of $N$ balls of
radius $\epsilon$.

######  Theorem 2.5.

Let $f$ be a mean-zero Gaussian process on a compact set
$T\subset\mathbb{R}^{n}$ with Lipschitz covariance function $C(x,y)$. Then,
there exists a constant $C>0$ such that

$\mathbb{E}\sup_{x,y\in
T}\left|\|h^{k}(x)-h^{k}(y)\|^{2}-D_{C}(x,y)^{2}\right|\leq\frac{C\operatorname{diam}(T,D_{C})}{\sqrt{k}}\int_{0}^{\operatorname{diam}(T,D_{C})}\log(\mathcal{N}(T,D_{C};\epsilon)+1)d\epsilon.$

We note that the factor of $\operatorname{diam}(T,D_{C})$ on the right-hand-
side provides the correct scaling with respect to the distance in that
replacing the covariance function $C(x,y)$ by a constant multiple of itself
$rC(x,y)$ results in both sides scaling by a factor of $r^{2}$.

To prove Theorem 2.5, we will rely on a couple basic properties of covering
numbers and diameters stated below.

######  Proposition 2.6.

Let $(T,d)$ be a pseudometric space of finite diameter.

  1. 1.

If $d^{\prime}$ is a pseudometric on $T$ such that $d^{\prime}(x,y)=Cd(x,y)$
for some fixed constant $C>0$ and all $x,y\in T$, then

$\operatorname{diam}(T,d^{\prime})=C\operatorname{diam}(T,d)$ (12)

and

$\mathcal{N}(T,d^{\prime},\epsilon)=\mathcal{N}(T,d,\epsilon/C)$ (13)

for all $\epsilon>0$.

  2. 2.

Define a pseudometric $d^{\prime\prime}$ on $T\times T$ by
$d^{\prime\prime}((x_{1},y_{1}),(x_{2},y_{2}))=\sqrt{d(x_{1},x_{2})^{2}+d(y_{1},y_{2})^{2}}$.
Then,

$\operatorname{diam}(T\times
T,d^{\prime\prime})=\sqrt{2}\operatorname{diam}(T,d)$ (14)

and

$\mathcal{N}(T\times
T,d^{\prime\prime},\epsilon\sqrt{2})\leq\mathcal{N}(T,d,\epsilon)^{2}.$ (15)

###### Proof.

We may prove and (12) and (13) by simply observing that a ball of radius
$\epsilon$ in the $d^{\prime}$ pseudometric is a ball of radius $\epsilon/C$
in the original $d$ pseudometric.

Choose $x,y\in T$ such that $d(x,y)\geq\operatorname{diam}(T,d)-\delta$, where
$\delta>0$. Then,

$d^{\prime\prime}((x,x),(y,y))=\sqrt{2}d(x,y)\geq\sqrt{2}\operatorname{diam}(T,d)-\sqrt{2}\delta.$

As $\delta>0$ is arbitrary, $\operatorname{diam}(T\times
T,d^{\prime\prime})\geq\sqrt{2}\operatorname{diam}(T,d)$.

Similarly, if $(x_{1},y_{1}),(x_{2},y_{2})\in T\times T$ are arbitrary, then
$d(x_{1},x_{2}),d(y_{1},y_{2})\leq\operatorname{diam}(T,d)$ and

$d^{\prime\prime}((x_{1},y_{1}),(x_{2},y_{2}))\leq\sqrt{\operatorname{diam}(T,d)^{2}+\operatorname{diam}(T,d)^{2}}=\sqrt{2}\operatorname{diam}(T,d).$

Thus, $\operatorname{diam}(T\times
T,d^{\prime\prime})\leq\sqrt{2}\operatorname{diam}(T,d)$ and
$\operatorname{diam}(T\times
T,d^{\prime\prime})=\sqrt{2}\operatorname{diam}(T,d)$.

Let $B_{1},\ldots,B_{N}$ be a collection of balls of radius $\epsilon$ and
centers $x_{1},\ldots,x_{N}$ whose union is $T$, where
$N=\mathcal{N}(T,d,\epsilon)$. We claim

$\bigcup_{1\leq i,j,\leq N}B_{ij}=T\times T,$

where $B_{ij}$ is the ball of radius $\epsilon\sqrt{2}$ and center
$(x_{i},x_{j})\in T\times T$ under the pseudometric $d^{\prime\prime}$.

Let $(x,x^{\prime})\in T\times T$. Then, there exists $1\leq i,j\leq N$ such
that $x\in B_{i}$ and $x^{\prime}\in B_{j}$. By definition of
$d^{\prime\prime}$, $(x,x^{\prime})\in B_{ij}$. There are $N^{2}$ such
$B_{ij}$, proving our claim. ∎

We say a random variable $X$ is sub-gaussian if there exists $K>0$ such that

$\mathbb{E}\exp(X^{2}/K^{2})\leq 2,$

in which case $\|X\|_{\psi_{2}}$ is the smallest such $K$.

Similarly, $X$ is sub-exponential if there exists $K^{\prime}>0$ such that

$\mathbb{E}\exp(|X|/K^{\prime})\leq 2,$

in which case $\|X\|_{\psi_{1}}$ is the smallest such $K^{\prime}$.

We will use two basic results on sub-gaussian and sub-exponential random
variables, as well as a corollary of Bernstein’s inequality.

######  Lemma 2.7.

If $X$ and $Y$ are sub-gaussian random variables, then

$\|XY\|_{\psi_{1}}\leq\|X\|_{\psi_{2}}\|Y\|_{\psi_{2}}.$ (16)

######  Lemma 2.8.

There exists $C>0$ such that if $X$ is a sub-exponential random variable, then

$\|X-\mathbb{E}X\|_{\psi_{1}}\leq C\|X\|_{\psi_{1}}.$ (17)

Lemmas 2.7 and 2.8 are Lemma 2.7.7 and Exercise 2.7.10 in [21], respectively.

In order to handle sums of subexponential random variables, we will rely on
the following form of Bernstein’s inequality (see Corollary 2.8.3 in [21]).

######  Theorem 2.9 (Bernstein’s Inequality).

Let $X_{1},\ldots,X_{N}$ be independent, mean-zero, sub-exponential random
variables. Then, there exists an absolute constant $c>0$ such that for all
$t\geq 0$,

$P\left\\{\frac{1}{N}\left|\sum_{i=1}^{N}X_{i}\right|\geq t\right\\}\leq
2\exp\left[-c\min\left(\frac{t^{2}}{K^{2}},\frac{t}{K}\right)N\right].$

######  Corollary 2.10.

Let $X_{1},\ldots,X_{N}$ be i.i.d., mean-zero, sub-exponential random
variables. Then,

$\left\|\sum_{i=1}^{N}X_{i}\right\|_{\psi_{1}}\leq\frac{CK}{\sqrt{N}},$ (18)

where $K=\|X_{1}\|_{\psi_{1}}=\cdots=\|X_{N}\|_{\psi_{1}}$.

###### Proof.

By Bernstein’s inequality (Theorem 2.9, there exists $c>0$ such that

$P\left\\{\left|\sum_{i=1}^{N}X_{i}\right|\geq t\right\\}\leq
2\exp\left[-c\min\left(\frac{t^{2}}{K^{2}},\frac{t}{K}\right)N\right].$ (19)

By Proposition 2.7.1 of [21], an equivalent (up to a constant multiple)
definition of $\|\sum_{i=1}^{N}X_{i}\|_{\psi_{1}}$ is the smallest number
$K^{\prime}$ such that

$P\left(\left|\sum_{i=1}^{N}X_{i}\right|\geq t\right)\leq
2\exp(-t/K^{\prime})$ (20)

for all $t>0$. We claim (20) holds with

$K^{\prime}=\frac{1}{\sqrt{c\log(2)}}\frac{K}{\sqrt{N}}.$ (21)

If $t\leq K^{\prime}\log(2)$, then $2\exp(-t/K^{\prime})\geq 1$, so (20) holds
trivially. So suppose going forward that

$t\geq K^{\prime}\log(2).$ (22)

If $t^{2}/K^{2}\geq t/K$, i.e., $t\geq K$, then by (19)

$P\left(\left|\sum_{i=1}^{N}X_{i}\right|\geq t\right)\leq 2\exp(-cNt/K)\leq
2\exp(-t/K^{\prime})$ (23)

since $K^{\prime}\geq\frac{K}{N}\frac{1}{\sqrt{c}\log(2)}\geq c\frac{K}{N}$,
as $c^{3/2}\leq\frac{1}{\sqrt{\log(2)}}$. (One may check in the proof of
Bernstein’s inequality that $c<1$, or simply observe that one may take $c<1$
regardless.)

If $t^{2}/K^{2}\leq t/K$, i.e., $t\leq K$, then by (19)

$P\left(\left|\sum_{i=1}^{N}X_{i}\right|\geq t\right)\leq
2\exp(-cNt^{2}/K^{2}).$ (24)

We want to show that $cNt^{2}/K^{2}\geq t/K^{\prime}$,or

$cNt/K^{2}\geq 1/K^{\prime}.$

By (22),

$cNt/K^{2}\geq cNK^{\prime}\log(2)/K^{2},$

so it suffices to show that

$cNK^{\prime}\log(2)/K^{2}\geq 1/K^{\prime},$

or

$(K^{\prime})^{2}\geq\frac{K^{2}}{cN\log(2)}.$ (25)

Equality holds in (25) by the definition of $K^{\prime}$ in 21, so we are
done. One obtains the arbitrary constant in (18) by combining
$\frac{1}{\sqrt{c\log(2)}}$ with the constant implicit in Proposition 2.7.1 of
[21]. ∎

######  Theorem 2.11 (Dudley Inequality).

Let $(f_{t})_{t\in T}$ be a mean-zero random process on a pseudometric space
$(T,d)$ such that

$\|f_{s}-f_{t}\|_{\Psi_{1}}\leq d(s,t)\text{ for all }s,t,\in T.$

Then,

$\mathbb{E}\sup_{t\in T}|X_{t}|\leq
8\int_{0}^{\operatorname{diam}(T,d)}\log(\mathcal{N}(T,d;\epsilon)+1)d\epsilon.$

Theorem 2.11 follows immediately from the proof of Theorem 11.1 in [15] by
taking the random process to be mean-zero, and restricting to the special case
of subexponential norms.

###### Proof of Theorem 2.5.

Recall that

$h^{k}(x)=\frac{1}{\sqrt{k}}(f_{1}(x),\ldots,f_{k}(x)),$

where each $f_{i}$ is an independent realization of the given Gaussian process
$f$.

Define a random process $(X_{xy})_{(x,y)\in T\times T}$ by

$\displaystyle X_{xy}$
$\displaystyle=\|h^{k}(x)-h^{k}(y)\|^{2}-d_{C}(x,y)^{2}$
$\displaystyle=\frac{(f_{1}(x)-f_{1}(y))^{2}+\ldots+(f_{k}(x)-f_{k}(y))^{2}}{k}-D_{C}(x,y)^{2}.$

By (9) and the independence of the $f_{i}$,

$f_{i}(x)=\sum_{j\geq 0}\xi_{ij}\sqrt{\lambda_{j}}\varphi(x),$

where the $\xi_{ij}$ ($1\leq i\leq k,j\geq 0$) are i.i.d. $N(0,1)$ random
variables.

Thus,

$X_{st}=\frac{1}{k}\sum_{i=1}^{k}\left(\sum_{j\geq
0}\xi_{ij}\sqrt{\lambda_{j}}(\varphi(x)-\varphi(y))\right)^{2}-D_{C}(x,y)^{2}$

and, factoring the difference of squares,

$\displaystyle X_{xy}-X_{wz}$
$\displaystyle=\frac{1}{k}\sum_{i=1}^{k}\left(\sum_{j\geq
0}\xi_{ij}\sqrt{\lambda_{j}}(\varphi(x)-\varphi(y))\right)^{2}-\left(\sum_{j\geq
0}\xi_{ij}\sqrt{\lambda_{j}}(\varphi(w)-\varphi(z))\right)^{2}-(D_{C}(x,y)^{2}-D_{C}(w,z)^{2})$
$\displaystyle=\frac{1}{k}\sum_{i=1}^{k}\left(\sum_{j\geq
0}\xi_{ij}\sqrt{\lambda_{j}}(\varphi(x)-\varphi(y)-\varphi(w)+\varphi(z))\right)\left(\sum_{j\geq
0}\xi_{ij}\sqrt{\lambda_{j}}(\varphi(x)-\varphi(y)+\varphi(w)-\varphi(z))\right)$
$\displaystyle-(D_{C}(x,y)^{2}-D_{C}(w,z)^{2}).$

For each $1\leq i\leq k$,

$Y_{1i}:=\left(\sum_{j\geq
0}\xi_{ij}\sqrt{\lambda_{j}}(\varphi(x)-\varphi(y)-\varphi(w)+\varphi(z))\right)\sim
N\left(0,\sum_{j\geq
0}\lambda_{j}(\varphi(x)-\varphi(y)-\varphi(w)+\varphi(z))^{2}\right)$

and

$Y_{2i}:=\left(\sum_{j\geq
0}\xi_{ij}\sqrt{\lambda_{j}}(\varphi(x)-\varphi(y)+\varphi(w)-\varphi(z))\right)\sim
N\left(0,\sum_{j\geq
0}\lambda_{j}(\varphi(x)-\varphi(y)+\varphi(w)-\varphi(z))^{2}\right).$

By the inequality $(a-b)^{2}\leq 2(a^{2}+b^{2})$, we have

$\displaystyle\sum_{j\geq
0}\lambda_{j}(\varphi(x)-\varphi(y)-\varphi(w)+\varphi(z))^{2}$
$\displaystyle=\sum_{j\geq
0}\lambda_{j}[(\varphi(x)-\varphi(w))-(\varphi(y)-\varphi(z))]^{2}$
$\displaystyle\leq 2\sum_{j\geq
0}\lambda_{j}[(\varphi(x)-\varphi(w))^{2}+(\varphi(y)-\varphi(z))^{2}]$
$\displaystyle=2(D_{C}(x,w)^{2}+D_{C}(y,z)^{2}).$

Thus,

$\|Y_{1i}\|_{\psi_{2}}\leq\sqrt{2(D_{C}(x,w)^{2}+D_{C}(y,z)^{2})}$

and similarly,

$\|Y_{2i}\|_{\psi_{2}}\leq\sqrt{2(D_{C}(x,z)^{2}+D_{C}(y,w)^{2})}.$

As a result, by (16) and (17),

$\displaystyle\|Y_{1i}Y_{2i}-\mathbb{E}Y_{1i}Y_{2i}\|_{\psi_{1}}$
$\displaystyle\leq C\|Y_{1i}Y_{2i}\|_{\psi_{1}}$ $\displaystyle\leq
C\|Y_{1i}\|_{\psi_{2}}\|Y_{2i}\|_{\psi_{2}}$
$\displaystyle=C\sqrt{(D_{C}(x,w)^{2}+D_{C}(y,z)^{2}}\sqrt{(D_{C}(x,z)^{2}+D_{C}(y,w)^{2}}$
$\displaystyle\leq C\operatorname{diam}(T,D_{C})d_{T\times T}((x,y),(w,z)),$

where

$d_{T\times T}((x,y),(w,z)):=\sqrt{(D_{C}(x,w)^{2}+D_{C}(y,z)^{2}}.$

By definition,

$X_{xy}-X_{wz}=\frac{1}{k}\sum_{i=1}^{k}Y_{1i}Y_{2i}-\mathbb{E}Y_{1i}Y_{2i},$

the average of $k$ independent, mean-zero, sub-exponential random variables.
By Corollary 2.10 and the above,

$\|X_{xy}-X_{wz}\|_{\psi_{1}}\leq\frac{C\operatorname{diam}(T,D_{C})}{\sqrt{k}}d_{T\times
T}((x,y),(w,z)).$

Defining a pseudometric on $T\times T$ by

$d^{*}((x,y),(w,z))=\frac{C\operatorname{diam}(T,D_{C})}{\sqrt{k}}d_{T\times
T}((x,y),(w,z))$

we may apply Dudley’s inequality and obtain

$\mathbb{E}\sup_{(x,y)\in T\times T}|X_{xy}|\leq
8\int_{0}^{\operatorname{diam}(T\times T,d^{*})}\log(\mathcal{N}(T\times
T,d^{*};\epsilon)+1)d\epsilon.$

For ease of notation, set $r=\frac{C\operatorname{diam}(T,D_{C})}{\sqrt{k}}$.

By (12), (13), and the change of variables $\epsilon\to\epsilon/r$,

$\displaystyle\mathbb{E}\sup_{(x,y)\in T\times T}|X_{xy}|$ $\displaystyle\leq
8\int_{0}^{r\operatorname{diam}(T\times T,d_{T\times
T})}\log(\mathcal{N}(T\times T,d_{T\times T};\epsilon/r)+1)d\epsilon$
$\displaystyle=8r\int_{0}^{\operatorname{diam}(T\times T,d_{T\times
T})}\log(\mathcal{N}(T\times T,d_{T\times T};\epsilon)+1)d\epsilon.$

By (14), (15), and the change of variables $\epsilon\to\epsilon/\sqrt{2}$,

$\displaystyle\mathbb{E}\sup_{(x,y)\in T\times T}|X_{xy}|$ $\displaystyle\leq
8r\int_{0}^{\sqrt{2}\operatorname{diam}(T,D_{C})}\log\left((\mathcal{N}(T,D_{C};\epsilon/\sqrt{2})+1)^{2}\right)d\epsilon$
$\displaystyle=16\frac{C\operatorname{diam}(T,D_{C})}{\sqrt{k}}\int_{0}^{\sqrt{2}\operatorname{diam}(T,D_{C})}\log(\mathcal{N}(T,D_{C};\epsilon/\sqrt{2})+1)d\epsilon$
$\displaystyle=16\sqrt{2}\frac{C\operatorname{diam}(T,D_{C})}{\sqrt{k}}\int_{0}^{\operatorname{diam}(T,D_{C})}\log(\mathcal{N}(T,D_{C};\epsilon)+1)d\epsilon.$

∎

##  3 Discretization of Gaussian Process Embeddings

Let $T=\\{x_{1},\ldots,x_{N}\\}\subset\mathbb{R}^{n}$. Fix a symmetric,
positive semidefinite, $N\times N$ matrix $A$ and let
$\lambda_{1},\ldots,\lambda_{N}$ denote the eigenvalues of $A$ with
corresponding orthonormal basis of eigenvectors $v_{1},\ldots,v_{N}$. $A$ has
a unique symmetric, positive semidefinite square root we denote $\sqrt{A}$.

If $g\sim N(0,I_{N})$, where $I_{N}$ is the $N\times N$ identity matrix, then

$f:=\sqrt{A}g$

defines a Gaussian process $f:=(f_{i})_{1\leq i\leq N}$ on $T$ with covariance
matrix $A$, where the $i$-th entry of $f$ is associated with $x_{i}$. In
practice, $A$ will depend on the $x_{i}$ explicitly.

We denote the canonical metric of the Gaussian process by

$d_{A}(x_{i},x_{j})=\mathbb{E}[(f_{i}-f_{j})^{2}]^{1/2}.$

To generate a Gaussian process embedding of $X$ into $\mathbb{R}^{k}$, we
compute the $N\times k$ matrix

$Y=\frac{1}{\sqrt{k}}\sqrt{A}G,$

where $G$ is an $N\times k$ matrix of i.i.d. $N(0,1)$ entries, and the rows of
$Y$, denoted $y_{1},\ldots,y_{N}\in\mathbb{R}^{k}$ are the images of
$x_{1},\ldots,x_{N}$, respectively.

The diffusion distance with respect to the matrix $A$ at time $t$ is given by

$d_{A}^{(t)}(x_{i},x_{j})^{2}=\sum_{i=k}^{N}\lambda_{k}^{2t}[(v_{k})_{i}-(v_{k})_{j}]^{2}=\|a_{i}-a_{j}\|^{2},$
(26)

where $v_{i}$ denotes the $i$-th entry of the vector $v$ and
$a_{1},\ldots,a_{N}$ are the rows of $A$.

######  Theorem 3.1 (Discrete Karhunen-Loève Expansion).

Let $f=(f_{i})_{1\leq i\leq N}$, $\lambda_{1},\ldots,\lambda_{N}$, and
$v_{1},\ldots,v_{N}$ be as above. Then,

$f=\sum_{i=1}^{\infty}\xi_{i}\sqrt{\lambda_{i}}v_{i},$ (27)

where $\xi_{1},\ldots,\xi_{N}$ are i.i.d. $N(0,1)$ random variables.

While Theorem 3.1 follows immediately from Theorem 2.3, we note that one may
prove it directly with a simple application of the diagonalization properties
of symmetric matrices and the rotational invariance of the standard
multivariate normal $g$.

######  Proposition 3.2.

Let $T$, $A$, $Y$, $f=(f_{i})_{1\leq i\leq N}$,
$\lambda_{1},\ldots,\lambda_{N}$, and $v_{1},\ldots,v_{N}$ be as above. Fix
$x_{i},x_{j},x_{l}\in T$. Then,

$\mathbb{E}\langle
y_{i}-y_{j},y_{i}-y_{k}\rangle=\sum_{p=1}^{N}\lambda_{p}[(v_{p})_{i}-(v_{p})_{j})((v_{p})_{i}-(v_{p})_{l}].$
(28)

In particular,

$d_{A}(x_{i},x_{j})^{2}=\mathbb{E}\|y_{i}-y_{j}\|^{2}=\sum_{p=1}^{N}\lambda_{p}[(v_{p})_{i}-(v_{p})_{j})]^{2}=d_{A}^{1/2}(x_{i},x_{j})^{2},$
(29)

More generally, $d_{A}^{(t)}=d_{A^{2t}}$ for all $t>0$.

While Proposition 3.2 follows immediately from Proposition 2.4, we note that
it may also be proven easily from the elementary property that
$\mathbb{E}\langle v,g\rangle\langle w,g\rangle=\langle v,w\rangle$ for
$v,w\in\mathbb{R}^{N}$ and $g\sim N(0,I_{N})$.

A discrete version of Theorem 2.5 is as follows:

######  Theorem 3.3.

Let $T$, $A$, $f=(f_{i})_{1\leq i\leq N}$, be as above. Then,

$\mathbb{E}\sup_{1\leq i,j,\leq
N}\left|\|y_{i}-y_{j}\|^{2}-d_{A}(x_{i},x_{j})^{2}\right|\leq\frac{C\operatorname{diam}(T,d_{A})}{\sqrt{k}}\int_{0}^{\operatorname{diam}(T,d_{A})}\log(\mathcal{N}(T,d_{A};\epsilon)+1)d\epsilon.$

As data can be noisy, it would be nice to know that small perturbations in
measurement result in small differences in Gaussian process embeddings.
However, if one repeats the same Gaussian process embedding with the same
covariance matrix (i.e., zero noise), then one may expect entirely different
results. While the two embeddings may have similar behavior, it may be
difficult to measure their similarities or differences.

Thus, we choose to compare two realizations of the Gaussian process embeddings
which are computed by sketching with the same Gaussian matrix.

######  Proposition 3.4 (Robustness of Gaussian process embeddings).

Let $A^{2},B^{2}$ be symmetric positive semidefinite $N\times N$ matrices
which serve as covariance matrices for Gaussian process embeddings in
$\mathbb{R}^{k}$ $Y=\frac{1}{\sqrt{k}}AG$ and $Z=\frac{1}{\sqrt{k}}BG$,
respectively.

Let $\mu_{1},\ldots,\mu_{N}$ denote the eigenvalues of $A-B$. Then,

$\sum_{i=1}^{N}\|y_{i}-z_{i}\|_{2}^{2}\sim\sum_{i=1}^{N}\mu_{i}^{2}\frac{\chi_{i}^{2}(k)}{k},$

where $\chi^{2}(k)$ is the chi-squared distribution with $k$ degrees of
freedom and $\chi^{2}_{1}(k),\ldots,\chi^{2}_{N}(k)$ are independent copies of
$\chi^{2}(k)$. In particular,

$\mathbb{E}\sum_{i=1}^{N}\|y_{i}-z_{i}\|_{2}^{2}=\|A-B\|_{F}^{2}$

and

$\mathbb{V}\sum_{i=1}^{N}\|y_{i}-z_{i}\|_{2}^{2}=2\|A-B\|_{F}^{2}/k,$

where $\|M\|_{F}$ refers to the Frobenius norm of the matrix $M$ and
$\mathbb{V}X$ is the variance of a random variable $X$.

###### Proof.

Denote the columns of $G$ by $g_{1},\ldots,g_{k}\in\mathbb{R}^{N}$ and let
$a_{i}$ and $b_{i}$ be the $i$-th columns of $A$ and $B$, respectively. Then,

$\displaystyle\sum_{i=1}^{N}\|y_{i}-z_{i}\|_{2}^{2}$
$\displaystyle=\frac{1}{k}\sum_{i=1}^{N}\sum_{j=1}^{N}(a_{i}\cdot
g_{j}-b_{i}\cdot g_{j})^{2}$
$\displaystyle=\frac{1}{k}\sum_{j=1}^{N}\sum_{i=1}^{N}[(a_{i}-b_{i})\cdot
g_{j}]^{2}$ $\displaystyle=\frac{1}{k}\sum_{j=1}^{N}\|(A-B)g_{j}\|^{2}$
$\displaystyle=\frac{1}{k}\sum_{j=1}^{N}\sum_{i=1}^{N}\mu_{i}^{2}\xi_{ij}^{2},$

where $\xi_{ij}$ are i.i.d. $N(0,1)$ random variables. Again swapping the
order of summation,

$\displaystyle\sum_{i=1}^{N}\|y_{i}-z_{i}\|_{2}^{2}$
$\displaystyle=\frac{1}{k}\sum_{i=1}^{N}\sum_{j=1}^{N}\mu_{i}^{2}\xi_{ij}^{2}$
$\displaystyle=\sum_{i=1}^{N}\mu_{i}^{2}\frac{\chi_{i}^{2}(k)}{k}.$

∎

We are particularly interested in the case where $A$ and $B$ are
approximations of the heat kernel on two finite sets, one of which may be a
perturbation of the other. Towards this end, it is notable that there are
results in the literature bounding $\|A-B\|_{F}$ for two normalized affinity
matrices $A$ and $B$ formed from close sets of points. In particular, see [13]
for robustness of the bistochastically normalized affinity matrix with respect
to heteroskedastic noise.

We also observe that the quantity $\|A-B\|_{F}^{2}$ is a discretization of
global diffusion distance between the pairs of points associated with $A$ and
$B$ in [6].

##  4 Algorithm

###  4.1 Approximating the Heat Kernel

Here, we attempt to approximate the heat kernel on a manifold from a finite
sample. This approximate heat kernel will take the role of the kernel $a$ in
Subsection 2.1 in determining the diffusion maps embedding and diffusion
distance.

We choose the heat kernel so that the straight-line distances in our embedding
coincide with the commonly-used diffusion distances associated with said heat
kernel. We could attempt to mimic Theorem 2.2 and choose an arbitrary
isotropic kernel. However, in working with a finite set it is difficult to
work with geodesic distances so we restrict our motivation to the results on
straight-line distances.

Let $X=\\{x_{1},\ldots,x_{N}\\}\subset\mathbb{R}^{n}$. Fix $\epsilon>0$ and
define a kernel matrix $K_{\epsilon}$ with $(i,j)$-th entry

$k_{\epsilon}(x_{i},x_{j})=h(|x_{i}-x_{j}|/\epsilon),$

where $h:\mathbb{R}\to\mathbb{R}$ is a smooth function with exponential decay.
A common choice is $h(x)=e^{-|x|^{2}/\epsilon}$. Define a vector

$q_{\epsilon}(x_{i})=\sum_{j=1}^{N}k_{\epsilon}(x_{i},x_{j})$

and replace the original kernel $k_{\epsilon}$ with

$\tilde{k}_{\epsilon}(x_{i},x_{j})=\frac{k_{\epsilon}(x_{i},x_{j})}{q_{\epsilon}(x_{i})q_{\epsilon}(x_{j})}.$

Now let

$v_{\epsilon}^{2}(x_{i})=\sum_{j=i}\tilde{k}_{\epsilon}(x_{i},x_{j})$

and define

$A_{\epsilon}(x_{i},x_{j})=\frac{1}{v_{\epsilon}(x_{i})^{2}}\tilde{k}_{\epsilon}(x_{i},x_{j}).$

The matrix $A_{\epsilon}$ is row-stochastic and, when $X$ is sampled from a
manifold $M\subset\mathbb{R}^{D}$, has been shown to approximate the heat
kernel on $M$. This holds even when the sampling is nonuniform [7]. However,
unlike the heat kernel on a manifold, it is generally not symmetric.

We would prefer to use a symmetric kernel for the purpose of Gaussian process
embeddings. (Sketching $A_{\epsilon}$ will give us a version of (27) with its
left singular vectors, not its eigenvectors, which may not even be
orthogonal.) Thus, we define the symmetrized kernel

$\tilde{A}_{\epsilon}(x_{i},x_{j})=\frac{\tilde{k}_{\epsilon}(x_{i},x_{j})}{v_{\epsilon}(x_{i})v_{\epsilon}(x_{j})}.$

The construction of $\tilde{A}_{\epsilon}$ is summarized in Algorithm 1.

Algorithm 1 Symmetric Normalized Affinity Matrix

Input: Finite subset $X=\\{x_{1},\ldots,x_{N}\\}\subset\mathbb{R}^{n}$, scale
$\epsilon>0$

Form affinity matrix $K$ with entries
$K_{ij}=e^{-|x_{1}-x_{j}|^{2}/\epsilon}$.

Compute vector $q=(q_{1},\ldots,q_{N})$ of row sums
$q_{i}=\sum_{j=1}^{N}K_{ij}$.

Compute $\tilde{K}=\operatorname{diag}(q^{-1})K\operatorname{diag}(q^{-1})$,
where $q^{-1}=(1/q_{1},\ldots,1/q_{N})$.

Form vector $v=(v_{1},\ldots,v_{N})$ of row sums of $\tilde{K}$:
$v_{i}=\sum_{j=1}^{N}\tilde{K}_{ij}$, where $\tilde{K}_{ij}$ is the $(i,j)$-th
entry of $\tilde{K}$.

Return
$\tilde{A}=\operatorname{diag}(1/\sqrt{v})\tilde{K}\operatorname{diag}(1/\sqrt{v})$,
where $1/\sqrt{v}=(1/\sqrt{v_{1}},\ldots,1/\sqrt{v_{N}})$.

The matrix $\tilde{A}_{\epsilon}$ does not have the same provable
approximations to the heat kernel as $A_{\epsilon}$ when $X$ is sampled with
nonuniform density from a manifold. To obtain the best of both worlds– a
kernel which is symmetric and approximates the heat kernel– we use the
bistochastic normalization of $k_{\epsilon}$ [17].

Our goal is to find $d(x)$ such that

$b_{\epsilon}(x_{i},x_{j})=\frac{k_{\epsilon}(x_{i},x_{j})}{d(x_{i})d(x_{j})}$

is bistochastic.

Such a $d(x)$ may be found as follows. (See [17].) Given
$v=(v_{1},\ldots,v_{n})\in\mathbb{R}^{N}$, let $diag(v)$ be the $N\times N$
diagonal matrix with diagonal entries $v_{1},\ldots,v_{N}$.

We want to find a diagonal matrix $D$ such that

$D^{-1}K_{\epsilon}D^{-1}\mathbf{1}=\mathbf{1},$ (30)

where $\mathbf{1}$ is the all-ones vector in $\mathbb{R}^{N}$.

Set $D_{0}=I_{N}$ and inductively define

$D_{i+1}=\operatorname{diag}(K_{\epsilon}D^{-1}\mathbf{1}).$

The limit

$D=\lim_{i\to\infty}D_{i}^{1/2}D_{i+1}^{1/2}$

solves (30). Computationally, we terminate the construction of the $D_{i}$
when $D_{i}^{1/2}D_{i+1}^{1/2}$ solves (30) within some given tolerance.
Following [13], we use Algorithm 2. (In the pseudocode, vector operations such
as division and square roots are componentwise.)

In [17], is it shown that choosing $b_{\epsilon}$ to be bistochastic with
respect to $1/q_{\epsilon}$ leads to a provable approximation of the heat
kernel when the $x_{i}$ are sampled from a manifold via a nonuniform
distribution. However, this sacrifices our ability to compare with Gaussian
process embeddings, so we avoid it, instead choosing to sample from uniform
distributions in our experiments for fairest comparisons.

Algorithm 2 Bistochastic Normalized Affinity Matrix

Input: Finite subset $X=\\{x_{1},\ldots,x_{N}\\}\subset\mathbb{R}^{n}$, scale
$\epsilon>0$, tolerance $\delta>0$ (usually set to $10^{-8}$)

Form affinity matrix $K$ with entries $K_{ij}=e^{-|x_{1}-x_{j}|^{2}/\epsilon}$

Let $d_{0}=\operatorname{diag}(\mathbf{1})$,
$d_{1}=K\operatorname{diag}(d_{0})^{-1}\mathbf{1}$, and
$d_{2}=K\operatorname{diag}(d_{1})^{-1}\mathbf{1}$.

Set $i=0$.

while $\|d_{i+2}/d_{i}-\mathbf{1}\|_{\infty}>\delta$ do

$i\to i+1$.

Compute $d_{i+2}=K\operatorname{diag}(d_{i+1})^{-1}\mathbf{1}$.

end while

Let $d=d_{i+2}^{1/2}d_{i}^{1/2}$.

Return $B=\operatorname{diag}(d)^{-1}K\operatorname{diag}(d)^{-1}$.

Despite the advantages of the bistochastic normalization, we will use the
symmetric normalized kernel $\tilde{A}_{\epsilon}$ in our experiments as well,
as it is simpler and often faster to compute.

###  4.2 Diffusion Maps Algorithm

Now that we have a couple of methods for approximating the heat kernel on a
discrete subset of $\mathbb{R}^{n}$, we use them to formally compute the
diffusion maps embedding found in (5). Observe that in place of taking the
powers of the affinity matrix, we may instead take powers of its eigenvalues.

We make one modification to the formula in (5), dropping the use of the top
eigenfunction, as for the heat kernel this is a constant function which does
nothing to separate points. Upon discretizing, this holds for the bistochastic
normalization (the all-ones vector is trivially the top eigenvector). While
the top eigenvector for the symmetric normalization is $v_{\epsilon}$,
differing entries reflect differences in density more than differences in
location. While these differences do contribute to the diffusion distance, we
run experiments where points are sampled from manifolds via the uniform
distribution to minimize this contribution.

Alternatively, one could include the top eigenvector in the diffusion maps
embedding with the symmetric normalized kernel. However, this would require
going up to 3 dimensions to properly embed a circle rather than 2, so it
appears more fair discard the top eigenvector.

Algorithm 3 Diffusion Maps Embedding

Input: Finite subset $X=\\{x_{1},\ldots,x_{N}\\}\subset\mathbb{R}^{n}$, scale
$\epsilon>0$, $k\in\mathbb{N}$, $t>0$, normalization, (Optional: tolerance
$\delta>0$)

if normalization=symmetric then

Let $A$ be the output of Algorithm 1 with inputs $X$ and $\epsilon$.

else if normalization=bistochastic then

Let $A$ be the output of Algorithm 2 with inputs $X$, $\epsilon$, and
$\delta$.

end if

Compute $k+1$ highest eigenvalues
$\lambda_{0}\geq\lambda_{1}\geq\ldots\geq\lambda_{k}$ of $A$ and corresponding
eigenvectors $v_{0},v_{1},\ldots,v_{k}$.

for $j=1,\ldots,N$ do

Compute
$y_{j}=(\lambda_{1}^{t}(v_{1})_{j},\ldots,\lambda_{k}^{t}(v_{k})_{j}).$

end for

Return $\\{y_{1},\ldots,y_{N}\\}$.

###  4.3 Gaussian process embeddings Algorithm

We have set up Gaussian process embeddings for discrete sets in full
generality in Section 3. At this point, we choose to compute them using the
approximations of the heat kernel for finite subsets of $\mathbb{R}^{n}$.

Specifically, if $A$ is an $N\times N$ matrix representing the heat kernel at
time $t$ on a set $X=\\{x_{1},\ldots,x_{N}\\}\subset\mathbb{R}^{n}$, then we
compute the Gaussian process embedding using $A^{2}$ (the heat kernel at time
$2t$) as the covariance matrix. This allows us to simply sketch the matrix $A$
and not have to compute any square roots.

Here, unlike as in Algorithm 3, we must directly take powers of the heat
kernel matrix in order to approximate the heat kernel at large times.

On the surface, diffusion maps does not have this issue. At any time $t$, one
takes the same top eigenvectors. The individual coordinates may be reweighed,
but the basic topology of the embedding will remain the same. However,
computing diffusion maps efficiently (with respect to time) may still require
taking large powers. (See Subsection 4.5.)

Algorithm 4 Gaussian process embedding

Input: Finite subset $X=\\{x_{1},\ldots,x_{N}\\}\subset\mathbb{R}^{n}$, scale
$\epsilon>0$, $k\in\mathbb{N}$, integer $p>0$, normalization, (Optional:
tolerance $\delta>0$)

if normalization=symmetric then

Let $A$ be the output of Algorithm 1 with inputs $X$ and $\epsilon$.

else if normalization=bistochastic then

Let $A$ be the output of Algorithm 2 with inputs $X$, $\epsilon$, and
$\delta$.

end if

Compute $Y=\frac{1}{\sqrt{k}}A^{p}G$, where $G$ is an $N\times k$ matrix with
i.i.d. $N(0,1)$ entries.

for $j=1,...,N$ do

Let $y_{j}$ equal the $j$-th row of $Y$.

end for

Return $\\{y_{1},\ldots,y_{N}\\}$.

###  4.4 Alternate Forms of Sketching and Non-Gaussian Processes

Motivated by the use of sketching in Algorithm 4, we modify our algorithm by
sketching by a matrix of symmetric Bernoulli random variables in place of
$N(0,1)$ random variables. This generates a (non-Gaussian) random process with
the same covariance matrix ($A^{2}$) as before, provided the symmetric
Bernoulli random variables are pairwise independent.

The pairwise independence, which shows up as a hypothesis in the following
proposition, allows one to store random matrices more easily.

######  Proposition 4.1.

Let $A$ be an $N\times N$, symmetric, positive semidefinite matrix and $X$ be
a mean-zero random variable satisfying $\mathbb{E}X^{2}=1$. Let $G$ be an
$N\times k$ random matrix whose entries are pairwise independent and each
distributed according to $X$.

Define a random matrix $Y=\frac{1}{\sqrt{k}}AG$. Then, for all $1\leq
i,j,l\leq N$,

$\mathbb{E}(y_{i}-y_{j})\cdot(y_{i}-y_{l})=(a_{i}-a_{j})\cdot(a_{i}-a_{k}),$
(31)

where $y_{1},\ldots,y_{N}$ are the rows of $Y$ and $a_{1},\ldots,a_{N}$ are
the rows of $A$.

In other words, Proposition 3.2 holds with random processes which are not
necessarily Gaussian.

###### Proof.

Here, we use the notation $M_{ij}$ to denote the $(i,j)$-th entry of the
matrix $M$.

Expanding the left-hand-side of (31), we have

$\displaystyle\mathbb{E}(y_{i}-y_{j})\cdot(y_{i}-y_{l})$
$\displaystyle=\frac{1}{k}\mathbb{E}\sum_{p=1}^{k}(Y_{ip}-Y_{jp})(Y_{ip}-Y_{jp})$
$\displaystyle=\frac{1}{k}\mathbb{E}\sum_{p=1}^{k}\left(\sum_{q=1}^{N}A_{iq}G_{qp}-\sum_{q=1}^{N}A_{jq}G_{qp}\right)\left(\sum_{q=1}^{N}A_{iq}G_{qp}-\sum_{q=1}^{N}A_{lq}G_{qp}\right)$
$\displaystyle=\frac{1}{k}\sum_{p=1}^{k}\mathbb{E}\left(\sum_{q=1}^{N}(A_{iq}-A_{jq})G_{qp}\right)\left(\sum_{r=1}^{N}(A_{ir}-A_{lr})G_{rp}\right)$
$\displaystyle=\frac{1}{k}\sum_{p=1}^{k}\sum_{q=1}^{N}(A_{iq}-A_{jq})(A_{iq}-A_{lq})$
$\displaystyle=\sum_{q=1}^{N}(A_{iq}-A_{jq})(A_{iq}-A_{lq}),$

where in the penultimate step we apply the pairwise independence of the
$G_{lk}$ and the fact that each of them has variance 1. ∎

We apply Proposition 4.1 by computing another form of Gaussian process
embedding in Algorithm 5, in which the sketching of Algorithm 4 is done with a
matrix of i.i.d. symmetric Bernoulli random variables.

Algorithm 5 Gaussian process embedding II (Symmetric Bernoulli Version)

Input: Finite subset $X=\\{x_{1},\ldots,x_{N}\\}\subset\mathbb{R}^{n}$, scale
$\epsilon>0$, $k,p\in\mathbb{N}$, normalization, (Optional: tolerance
$\delta>0$)

if normalization=symmetric then

Let $A$ be the output of Algorithm 1 with inputs $X$ and $\epsilon$.

else if normalization=bistochastic then

Let $A$ be the output of Algorithm 2 with inputs $X$, $\epsilon$, and
$\delta$.

end if

Compute $Y=\frac{1}{\sqrt{k}}A^{p}G$, where $G$ is an $N\times k$ matrix with
symmetric Bernoulli entries.

for $j=1,...,N$ do

Let $y_{j}$ equal the $j$-th row of $Y$.

end for

Return $\\{y_{1},\ldots,y_{N}\\}$.

###  4.5 Computational Time

One may check by inspection that Algorithm 1 runs in time $O(N^{2})$, where
$N$ is the number of points. Algorithm 2 runs in time $O(IN^{2})$, where $I$
is the total number of iterations in the while loop. While there are various
results in the literature on the convergence rate of similar algorithms (for
instance, [10] and [8]) we include Algorithm 1 in large part for cases in
which these bounds may be large.

To discuss the runtime of Algorithms 3 and 4, it may help to review some
numerical linear algebra. Classical techniques may be used to compute the top
$k$ eigenvectors and eigenvalues of an $N\times N$ matrix $A$ in time
$O(N^{3})$. More recent algorithms in randomized numerical linear algebra can
do the same in time $O(kN^{2})$. The basic idea is to first compute $AG$,
where $G$ is an $N\times(k+p)$ matrix of i.i.d. $N(0,1)$ entries and apply
classical techniques on the resulting $N\times(k+p)$ matrix which is smaller
and requires less runtime. Here, $p$ is a small, fixed integer parameter.
Strong precision is guaranteed with high probability, provided $A$ has a
rapidly decaying spectrum.

In the case $A$ has a slowly decaying spectrum, one uses the power iteration
method. That is, by observing that $A^{m}$ has a rapidly decaying spectrum for
some large $m$, one can apply the above method on $A^{m}$, computing $A^{m}G$
one step at a time (first $AG$, then $A^{2}G=A(AG)$, etc.). This step runs in
time $O(kmN^{2})$, which is still faster than $O(N^{3})$ for large $N$.

See [18] for more on randomized methods in numerical linear algebra.

Now let $A$ refer to the output of either Algorithm 1 or Algorithm 2.

Let us first consider the case where one decides not to take large powers of
$A$. If $A$ has a rapidly decaying spectrum, then the sketch in Algorithm 4
requires an additional $O(kN^{2})$ amount of work. In Algorithm 3, one may
compute the top eigenvalues and eigenvectors in $O(kN^{2})$ steps using the
randomized techniques described above.

Now suppose $A$ has a slowly decaying spectrum. For Algorithm 3, one may
either compute the top eigenvectors in $O(N^{3})$ steps (via the slower
classical methods), or apply the power iteration method and take $O(kmN^{2})$
steps. In theory, one could run Algorithm 4 in $O(kN^{2})$ steps. However, the
slowly decaying spectrum, combined with the Karhunen-Loève expansion,
encourages one to take large powers of $A$ anyway so one would not do this in
practice. (We will see in a moment this takes the same $O(kmN^{2})$ steps as
diffusion in this case.)

If instead, one already desires to compute powers of $A$, the above reasoning
shows that both Algorithm 3 and 4 will take $O(kmN^{2})$ steps, where $m$ is
the desired power. Here, we assume $A^{m}$ will have a rapidly decaying
spectrum, even if $A$ does not.

We note that if one wants to conduct a multiscale analysis using either
embedding, then many of the steps may be done simultaneously for all scales at
once. For Algorithm 3, this is somewhat trivial, as a different scale means a
different rescaling of each individual coordinate. For Algorithm 4, this
follows from the computation of $A^{m}G$ as $A(A(\cdots(AG)))$ as in the power
iteration method.

##  5 Experiments

Except where otherwise specified, our experiments were run as follows. Fix a
probability distribution $q(x)$ on a submanifold $M\subset\mathbb{R}^{D}$.

Fix parameters $N,n,p,k_{\min},k_{\max}\in\mathbb{N}$ and $\epsilon>0$. For
each of $N$ trials, we sample $n$ points from $M$ via the distribution $q$ to
obtain a set $X$. From this same sample of points, we run a subset of the
following algorithms:

  * •

Diffusion maps embedding (Algorithm 3) with the symmetric normalization (DMS)

  * •

Diffusion maps embedding (Algorithm 3) with the bistochastic normalization
(DMB)

  * •

Gaussian process embedding (Algorithm 4) with the symmetric normalization
(GPS)

  * •

Gaussian process embedding (Algorithm 4) with the bistochastic normalization
(GPB)

  * •

Gaussian process embedding with sketching by a symmetric Bernoulli random
matrix (Algorithm 5) with the symmetric normalization (GPSBS)

  * •

Gaussian process embedding with sketching by a symmetric Bernoulli random
matrix (Algorithm 5) with the bistochastic normalization (GPSBB)

Each of the embeddings was computed with target dimensions
$k_{\min},k_{\min}+1,\ldots,k_{\max}$. The set $X$ and the parameters remained
the same across the choice of embedding and target dimension for fair
comparison. Furthermore, whenever a Gaussian matrix was generated for the
purpose of sketching, it was used for both instances of Algorithm 4. The
analogous statement holds for symmetric Bernoulli matrices and Algorithm 5.

From each embedding $f:X\to\mathbb{R}^{k}$, we computed the smallest $L$ such
that

$cd(x,y)\leq\|f(x)-f(y)\|_{2}\leq cLd(x,y)$ (32)

for all $x,y\in X$ and some $c>0$. This $L$ is like the usual biLipschitz
constant of the map $f$; however, it adjusts for rescaling. That is, if there
exists $\lambda>0$ such that $\|f(x)-f(y)\|=\lambda d(x,y)$ for all $x,y\in
X$, then we would like to recover a value of 1, rather than $\lambda$. The
optimal $L$ may be computed by dividing the maximum dilation by the minimum
dilation, that is,

$L=\frac{\max_{x\neq y}[\|f(x)-f(y)\|/d(x,y)]}{\min_{x\neq
y}[\|f(x)-f(y)\|/d(x,y)]}.$

Note that the distance $d$ above is unspecified. It will refer to either the
diffusion distance appropriate for the particular embedding or the Euclidean
distance on $\mathbb{R}^{D}$. In most cases, we will stick to the diffusion
distance as in (26), as it is provably approximated by both diffusion maps and
Gaussian process embeddings; our purpose is to measure how well. However, we
note that since, for embedded Riemannian manifolds, the metric induced by the
diffusion distance is the same as that induced by the Euclidean metric. Hence,
the difference is often mote.

The purpose of computing $L$ is to capture the effects of self-intersection
for Gaussian process embeddings (see Figure 1(b)) and the challenges described
in Subsection 2.2 for diffusion maps (see Figures 2(a) and 2(b)). Both
phenomena involve two points which are far away being sent very close to each
other in the embedding, resulting in a large $L$.

For each embedding, we compute $\log(L)$ (where $L$ is the optimal constant in
(32) and the logarithm is the natural logarithm) and determine the mean and
standard deviation of this value across all trials, fixing the method of
embedding and the target dimension $k$. The mean is plotted against $k$ with
errorbar equal to the standard deviation. Multiple methods are plotted on the
same axes for comparison.

For most experiments, the symmetric normalization and bistochastic
normalization for the heat kernel lead to very similar results. To avoid
clutter, we plot both only when they are graphically distinguishable. When
they are not, we plot the methods using the more commonly used symmetric
normalized kernels.

###  5.1 Manifolds

To begin with, we embed $S^{1}=\\{(x,y):x^{2}+y^{2}=1\\}$, with parameters
$N=200,n=300,p=8,k_{\min}=2,k_{\max}=8,$ and $\epsilon=0.25$. The distribution
on $S^{1}$ is taken to be the uniform measure. Results are plotted in Figure
1(a).

We see that in the case of $S^{1}$, DMS is significantly more effective at
approximating the diffusion distance than GPS. This phenomenon may be
explained in large part by looking at the case $k=2$. The top 2 nontrivial
eigenfunctions of the heat kernel on $S^{1}$ are $\sin(\theta)$ and
$\cos(\theta)$, which are enough to recover the basic structure of $S^{1}$,
even if more eigenfunctions are required to fully obtain the diffusion
distance.

However, a random curve in the plane may be self-intersecting, collapsing two
points together. In the discrete case simulated above, this corresponds to a
very high value of $L$ (see Figure 1(b)). A random curve in $\mathbb{R}^{3}$
may not self-intersect, but in higher dimensions Gaussian process embeddings
still lack the optimization of diffusion maps.

(a) BiLipschitz constants for embeddings of $S^{1}$

(b) A sample Gaussian process embedding of $S^{1}$ into $\mathbb{R}^{2}$

Figure 1: Embeddings of $S^{1}$

Clearly, Gaussian process embeddings are not effective when the target
dimension is $k=2$. However, a value of $L\approx 3$ in higher dimensions may
be suitable for some practical applications.

To demonstrate the advantages of Gaussian process embeddings, we turn to our
example from Subsection 2.2: $S^{1}\times rS^{1}$. Here, we take
$N=100,n=500,p=10,k_{\min}=2,k_{\max}=12$, and $\epsilon=0.3$.

We consider $S^{1}\times rS^{1}$ as a stretched version of the flat torus,
originally living in $\mathbb{R}^{4}$. It may be described by the
parameterized equations

$(u,v)\mapsto(\cos(u),\sin(u),r\cos(v),r\sin(v)).$

We also take $r=3.5$; here, theory predicts that the diffusion maps embedding
will significantly struggle for $k\leq 7$. Sampling on $S^{1}\times 3.5S^{1}$
is with respect to the uniform measure. Results are plotted in Figure 2(a).

(a) Embeddings of $S^{1}\times 3.5S^{1}$

(b) Embeddings of the Klein bottle

Figure 2: Comparison of methods for two manifolds

In dimension $k=2$, the torus is not embeddable and all methods struggle
somewhat equally. Beginning with dimension $k=3$ (when the torus becomes
topologically embeddable), the Gaussian process embeddings start to
distinguish themselves more strongly from the diffusion maps embeddings.
Diffusion maps catches up in higher dimensions when the eigenfunctions which
distinguish points of $S^{1}$ begin to be used.

An oft advantage of diffusion maps is that it denoises data by removing
smaller-scale phenomena; however, the smaller-scale phenomenon here is the
copy of $S^{1}$ in $S^{1}\times rS^{1}$, which should most likely not be
considered noise due to its exact structure.

Lastly, we note that the diffusion maps embedding with the bistochastically
normalized heat kernel fares worse than the one with the symmetric normalized
kernel. While this would normally be considered a disadvantage, in this case
it may be the reflection of a better estimator of the eigenfunctions of the
torus.

Similar, yet more extreme phenomena occur for the Klein bottle, parameterized
as

$(u,v)\mapsto\left((a+b\cos(v))\cos(u),(a+b\cos(v))\sin(u),b\sin(v)\cos(u/2),b\sin(v)\sin(u/2)\right).$

Here, we take $a=10$ and $b=5$ and sample points with respect to the uniform
measure on $(u,v)\in[0,2\pi)^{2}$.

For our parameters, we take $N=100,n=500,p=4,k_{\min}=3,k_{\max}=20$, and
$\epsilon=2$. Results are plotted in Figure 2(b).

###  5.2 Outliers

For our next example, we sample 198 points from the uniform distribution on
$S^{1}$ and adjoin the points $(0,3)$ and $(3,0)$ as outliers. Our parameters
are $N=100,n=200,p=4,k_{\min}=2,k_{\max}=5$, and $\epsilon=0.5$. Results are
depicted in Figure 3(a).

(a) Embeddings of $S^{1}$ plus two outliers

(b) Sample embedding of $S^{1}$ plus two outliers with DMS

(c) Sample embedding of $S^{1}$ plus two outliers with DMS, outliers removed

Figure 3: Embeddings with outliers

Here, the Gaussian process embedding methods significantly outperform the
diffusion maps embeddings in dimensions $k=2$ and $k=3$, before the latter
surpasses them for $k=4,5$. We include a sample embedding when $k=2$ for DMS
in Figure 3(b). The two original outliers remain outliers; however, the 198
points from $S^{1}$ are all concentrated at the origin. Removing the outliers
and zooming in, we obtain the image in Figure 3(c). While some of the
structure of $S^{1}$ remains, there are still issues with regard to self-
intersection and the scale relative to the distance from the outliers.

DMS has to “waste” its top 2 eigenvectors on separating the outliers from the
circle, leaving the circle to collapse on itself. Moving up to $k=3$ nearly
removes the self-intersection problem, but one additional eigenvector is not
enough to fully describe the circle.

Sample embeddings with $k=2$ for GPS may be found in Figures 4, 4, and 4.
While adding outliers does not remove the original self-intersection issue
around embedding $S^{1}$ in $\mathbb{R}^{2}$, the images demonstrate a
reasonable treatment of the outliers in the following sense. The scales of
distances between points of $S^{1}$ and between $S^{1}$ and the outliers may
be viewed simultaneously.

Figure 4: Three sample embeddings of $S^{1}$ plus two outliers with GPS

###  5.3 Sketching with Symmetric Bernoulli Matrices

We replot earlier results on the Klein bottle (Figure 5(a)) and Gaussian point
clouds (Figure 5(b)), this time for the GPS and GPSBS methods. Here, it
appears that both methods, i.e., both forms of sketching, perform similarly.
There does not appear to be any significant disadvantage to saving time or
memory by replacing i.i.d. Gaussians with i.i.d symmetric Bernoulli random
variables. In fact, when using the bistochastically normalized heat kernel
there appears to be a slight advantage to using GPSBB (hidden in the figure by
GPSBS) over GPB. However, more study is needed to determine if this difference
extends to other examples.

(a) Results for the Klein bottle with all 4 types of Gaussian process
embeddings

(b) Results for point clouds with all 4 types of Gaussian process embeddings,
comparison with Euclidean distance

Figure 5: Comparison of Gaussian process embeddings with i.i.d. Gaussian
matrices and with symmetric Bernoulli matrices

###  5.4 Multiscale Analysis

While diffusion maps has a natural multiscale nature, taking powers of the
heat kernel corresponds to taking powers of its eigenvalues while leaving the
eigenvectors the same, resulting in a simple rescaling of the embedding.

For Gaussian process embeddings, the effect is more complicated. For extremely
low powers of the heat matrix, many eigenvalues will be close to 1; by the
Karhunen-Loève expansion, one will expect a random vector with distribution
close to $N(0,I_{N})$. For extremely large powers, the lower eigenvalues will
disappear and one will expect a random vector nearly proportional to $rv_{1}$,
where $v$ is the top eigenvector of the heat kernel and $r\sim N(0,1)$.

To test this effect we ran experiments the same as before except the target
dimension $k$ was fixed while the power $p$ of the heat kernel took on the
values $p=2,4,8,\ldots,2^{P}$ for some natural number $P$. Furthermore, the
diffusion distance was replaced by the Euclidean distance in determining $L$.

Figure 6(a) shows the results for $S^{1}\times 3.5S^{1}$ with
$N=100,n=500,k=8,P=10$, and $\epsilon=.3$. As one can see, for low powers $p$,
GPS struggles nearly as much as DMS before being a clear improvement for
powers $p=16,32$. However, as $p$ increases, GPS gets progressively worse,
eventually surpassing DMS.

Figure 6(b) shows the results for $S^{1}$ with $N=100,n=300,k=2,P=8$, and
$\epsilon=.25$. Similar behavior occurs as with the prior example, though DMS
consistently functions better than GPS, as expected.

(a) Embeddings of $S^{1}\times 3.5S^{1}$ into $\mathbb{R}^{8}$, $p$ varying,
$L$ measured with Euclidean distance

(b) Embeddings of $S^{1}$ into $\mathbb{R}^{2}$, $p$ varying, $L$ measured
with Euclidean distance

Figure 6: A multiscale analysis of the embeddings

Acknowledgments. The authors would like to thank Ronald Coifman and Hanwen
Zheng for productive discussions related to this work. We also thank Jeremy
Hoskins for the title.

## References

  * [1] Robert J. Adler, Sunder Ram Krishnan, Jonathan E. Taylor, and Shmuel Weinberger.  Convergence of the reach for a sequence of Gaussian-embedded manifolds.  Probab. Theory Related Fields, 171(3-4):1045–1091, 2018.
  * [2] Robert J Adler, Jonathan E Taylor, et al.  Random fields and geometry, volume 80.  Springer, 2007.
  * [3] Jonathan Bates.  The embedding dimension of Laplacian eigenfunction maps.  Appl. Comput. Harmon. Anal., 37(3):516–530, 2014.
  * [4] Mikhail Belkin and Partha Niyogi.  Laplacian eigenmaps for dimensionality reduction and data representation.  Neural computation, 15(6):1373–1396, 2003.
  * [5] P. Bérard, G. Besson, and S. Gallot.  Embedding Riemannian manifolds by their heat kernel.  Geom. Funct. Anal., 4(4):373–398, 1994.
  * [6] Ronald R. Coifman and Matthew J. Hirn.  Diffusion maps for changing data.  Appl. Comput. Harmon. Anal., 36(1):79–107, 2014.
  * [7] Ronald R. Coifman and Stéphane Lafon.  Diffusion maps.  Appl. Comput. Harmon. Anal., 21(1):5–30, 2006.
  * [8] Promit Ghosal and Marcel Nutz.  On the convergence rate of sinkhorn’s algorithm, 2022.
  * [9] Peter W. Jones, Mauro Maggioni, and Raanan Schul.  Universal local parametrizations via heat kernels and eigenfunctions of the Laplacian.  Ann. Acad. Sci. Fenn. Math., 35(1):131–174, 2010.
  * [10] Philip A. Knight.  The Sinkhorn-Knopp algorithm: convergence and applications.  SIAM J. Matrix Anal. Appl., 30(1):261–275, 2008.
  * [11] Sunder Ram Krishnan, Jonathan E. Taylor, and Robert J. Adler.  The intrinsic geometry of some random manifolds.  Electron. Commun. Probab., 22:Paper No. 1, 12, 2017.
  * [12] Stephane S. Lafon.  Diffusion maps and geometric harmonics.  ProQuest LLC, Ann Arbor, MI, 2004.  Thesis (Ph.D.)–Yale University.
  * [13] Boris Landa, Ronald R. Coifman, and Yuval Kluger.  Doubly stochastic normalization of the Gaussian kernel is robust to heteroskedastic noise.  SIAM J. Math. Data Sci., 3(1):388–413, 2021.
  * [14] O. P. Le Maître and O. M. Knio.  Spectral methods for uncertainty quantification.  Scientific Computation. Springer, New York, 2010.  With applications to computational fluid dynamics.
  * [15] Michel Ledoux and Michel Talagrand.  Probability in Banach spaces.  Classics in Mathematics. Springer-Verlag, Berlin, 2011.  Isoperimetry and processes, Reprint of the 1991 edition.
  * [16] John M. Lee.  Introduction to smooth manifolds, volume 218 of Graduate Texts in Mathematics.  Springer, New York, second edition, 2013.
  * [17] Nicholas F. Marshall and Ronald R. Coifman.  Manifold learning with bi-stochastic kernels.  IMA J. Appl. Math., 84(3):455–482, 2019.
  * [18] Per-Gunnar Martinsson.  Randomized methods for matrix computations, 2019.
  * [19] Stephen Semmes.  On the nonexistence of bi-Lipschitz parameterizations and geometric problems about $A_{\infty}$-weights.  Rev. Mat. Iberoamericana, 12(2):337–410, 1996.
  * [20] K. T. Sturm.  Diffusion processes and heat kernels on metric spaces.  Ann. Probab., 26(1):1–55, 1998.
  * [21] Roman Vershynin.  High-dimensional probability, volume 47 of Cambridge Series in Statistical and Probabilistic Mathematics.  Cambridge University Press, Cambridge, 2018.  An introduction with applications in data science, With a foreword by Sara van de Geer.