# Diffusion Representation for Asymmetric Kernels

Alvaro Almeida Gomez  Antônio J. Silva Neto  Jorge P. Zubelli  IMPA, Est. D.
Castorina, 110, Jardim Botânico, Rio de Janeiro, 22460-320, Brazil  IPRJ-UERJ,
R. Bonfim 25, Nova Friburgo 28625-570, Brazil  Khalifa University, P.O. Box
127788, Abu Dhabi, UAE

###### Abstract

We extend the diffusion-map formalism to data sets that are induced by
asymmetric kernels. Analytical convergence results of the resulting expansion
are proved, and an algorithm is proposed to perform the dimensional reduction.

A coordinate system connected to the tensor product of Fourier basis is used
to represent the underlying geometric structure obtained by the diffusion-map,
thus reducing the dimensionality of the data set and making use of the speedup
provided by the two-dimensional Fast Fourier Transform algorithm (2-D FFT).

We compare our results with those obtained by other eigenvalue expansions, and
verify the efficiency of the algorithms with synthetic data, as well as with
real data from applications including climate change studies.

###### keywords:

Diffusion maps; Dimensional reduction; FFT; Asymmetric kernels

<EMAIL_ADDRESS>A. Gomez<EMAIL_ADDRESS>J. Silva Neto
<EMAIL_ADDRESS>P. Zubelli

## Highlights

  * 1.

A novel methodology for asymmetric-kernel data dimension reduction is
developed.

  * 2.

Dimension reduction is based on the highly efficient FFT algorithm in higher
dimensions.

  * 3.

Numerical evidence indicates that the tensor product of the FFT basis
performance is faster than eigenvalue-based methods.

  * 4.

Geometric features of complex data sets are revealed in key examples such as
the Möbius strip and the sphere.

  * 5.

The methodology is employed in meteorological applications to identify regions
of largest temperature variation.

##  1 Introduction

Data compression has been studied extensively in many applications. See [1,
2]. Several dimensionality reduction algorithms are based on the spectral
decomposition of symmetric linear operators which induce geometric structures
in the data set. A classical example of such operators is the Laplacian matrix
associated to an undirected graph [3]. In Ref. [4], the eigenvalue
decomposition of the Laplacian matrix is used to reduce the dimensionality of
data sets in such way that the local information is preserved. The diffusion-
map approach is based on using the symmetric normalized Laplacian matrix [5].
Compared with other dimensionality reduction methods diffusion-map assumes
that the data set resides in a lower dimensional manifold, and uses
approximations of the Laplace-Beltrami operator to reveal relevant parameters
in the data set.

The spectral decomposition theorem does not hold for integral operators with
asymmetric kernels. Therefore, we cannot use the diffusion-map framework to
represent diffusion distances induced by asymmetric kernels. Moreover, if we
use the spectral decomposition in some symmetric normalization of an
asymmetric kernel, its performance requires a computational complexity
$O(n^{3})$, for $n\times n$ matrices.

In order to reduce the above $O(n^{3})$ complexity and to deal with more
general kernels, we present a new framework to represent the diffusion
geometry induced by asymmetric kernels. We use the 2-D FFT to compute this
representation. The main advantage of using this representation is that
compared to the eigenvector representation, the computation time decreases. In
fact, for a matrix of dimension $n\times n$, the complexity of the 2-D FFT is
$O(n^{2}\times log(n))$. The choice of the 2-D FFT to represent the data set
structure is based on the fact that the Fourier basis diagonalizes the
Laplacian defined on the Euclidean Torus.

In this paper, we deal with data sets whose structure is induced by an
asymmetric kernel. Based on Refs. [6, 7, 8], we use alternative orthonormal
bases to reduce the dimensionality of the data sets in such way that the
geometric structure is preserved. Here, we work with the representation theory
of diffusion distances in the context of changing data, proposed in Ref. [9].
First, we find a representation form for the classical diffusion geometry, and
then we find a representation form for changing data. To do that, we start
with the case $t=1$ and then extend it for any time $t$.

This paper is organized as follows, in Section 2 we give a brief exposition of
the classical representation theory for diffusion distances proposed in Refs.
[5, 9, 10, 11]. In Section 3 we present our framework to represent diffusion
distances when the structure in the data set is induced by an asymmetric
kernel. Finally, in Sections 4 and 5, we show some experiments with
applications and results, as well as draw some conclusions and directions for
further research.

##  2 Classical Diffusion-Map Theory

In this section, we review some classical results on diffusion-map theory. We
define diffusion distances in a measure space, and we recall some results
related to the representation of such diffusion distances.

###  2.1 Diffusion distance

Assume that our data set $(X,\mu)$ is a measure space, and let $k:X\times
X\to\mathbb{R}_{\geq 0}$ be a non-negative symmetric kernel, which is used to
measure the local connectivity between two points $x$ and $y$. If $X$ is also
a metric space, the most classical example of these kernels is the Gaussian
kernel given by

$e^{-d^{2}(x,y)/2\sigma^{2}},$

where $d$ is the distance function, and $\sigma^{2}$ is the scaling parameter.
We define the associated Markov kernel $\rho(x,y)$ by

$\rho(x,y):=\frac{k(x,y)}{\sqrt{v(x)}\sqrt{v(y)}}\mbox{ ,}$ (2.1)

where $v(x)$ is the volume form defined as $v(x)=\int_{X}k(x,y)\,d\mu(y)$.
Assuming that the volume form never vanishes and that $k\in L^{2}(X\times X)$,
then the operator $A:\,L^{2}(X)\to\,L^{2}(X)$ given by

$A(f)(x):=\int_{X}\rho(x,y)\,f(y)\,d\mu(y),$ (2.2)

is compact and self-adjoint. According to the spectral theorem, we can write

$\rho(x,y)=\sum_{n\in\mathbb{N}}\lambda_{n}\phi_{n}(x)\ \phi_{n}(y)\mbox{ , }$

where $\\{\phi_{n},\lambda_{n}\\}$ is the spectral decomposition of the
operator $A$. For any natural number $t$, we define the diffusion distance at
time $t$ between two points $x$ and $y$ by

$D^{t}(x,y):={\|\,\rho^{t}(x,\cdot)-\rho^{t}(y,\cdot)\|}_{L^{2}(X)}\mbox{ ,}$
(2.3)

where $\rho^{t}$ is the kernel of the integral operator $A^{t}$. Here, $A^{t}$
is the composition of the operator $A$, a total of $t$ times. The kernel
$\rho^{t}(x,y)$ measures the probability that the points $x$ and $y$ are
connected by a path in time $t$. Observe that the distance $D^{t}$ is an
average over all of the paths in time $t$ connecting $x$ to $y$. Therefore,
the diffusion distance is robust to noisy data. We see that the quantity
$D^{t}$ is small when there are many paths of length $t$ connecting $x$ and
$y$. Using the spectral decomposition of the operator $A$, we can write the
diffusion distance as

$D^{t}(x,y)=\sqrt{\sum_{n\in\mathbb{N}}\lambda_{n}^{t}(\phi_{n}(x)-\phi_{n}(y))^{2}}\mbox{
.}$ (2.4)

The above expression allows us to reduce the dimensionality of the diffusion
geometry, namely, we embed our data set in a lower dimensional space using the
diffusion-map $\psi_{k}^{t}:X\,\to\mathbb{R}^{k}$, where

$\psi_{k}^{t}(x)=(\lambda_{i}^{t/2}\,\phi_{i}(x))_{i=1}^{k}\mbox{ .}$

###  2.2 Diffusion-map for changing data

The diffusion distance for changing data proposed in Ref. [9] compares data
points between parametric data sets. We define $X_{\alpha}$ as the data set
$X$ endowed with the kernel structure $k_{\alpha}$. As above, for each kernel
$k_{\alpha}$, we consider the associated Markov kernel $\rho_{\alpha}$ defined
in Eq. (2.1), and the operator $A_{\alpha}$ as in Eq. (2.2). To compare the
data structure $X_{\alpha}$ with $X_{\beta}$, we define the dynamic diffusion
distance $D^{t}:X_{\alpha}\times X_{\beta}\to\mathbb{R}_{\geq 0}$ by

$D^{t}(x_{\alpha},y_{\beta})={\|\,\rho_{\alpha}^{t}(x,\cdot)-\rho_{\beta}^{t}(y,\cdot)\|}_{L^{2}(X)}\mbox{
. }$

Furthermore, in Ref. [9], the global diffusion distance at time $t$ is defined
by

$\mathbb{D}^{t}(X_{\alpha},X_{\beta})=\sqrt{\int_{X}(D^{t}(x_{\alpha},x_{\beta}))^{2}\,d\mu(x)}\mbox{
. }$

The global diffusion distance measures the change from the data structure
$X_{\alpha}$ to $X_{\beta}$. Under mild assumptions on the data set $X$ and
the family of kernels $\\{k_{\alpha}\\}$, the dynamic diffusion distance and
the global diffusion distance can be computed using the spectral decomposition
of the operators $A_{\alpha}$, as in Eq. (2.4). See Ref. [9].

##  3 Diffusion Representation for Asymmetric Kernels

Assume that $(X,\mu)$ is a measure space, and consider an asymmetric kernel
$k$, which is any square integrable measurable function $k:X\times
X\to\mathbb{R}_{\geq 0}$. Observe that here we do not require that the kernel
$k$ be symmetric. As example of these kernels, we assume that $(X,d)$ is a
metric space, and consider the weighted Gaussian kernel defined by

$w(x,y)\;e^{-d^{2}(x,y)/2\sigma^{2}}\mbox{ , }$

where $w$ is the weight function. Weighted Gaussian kernels measure how
information is distributed locally from $x$ to $y$. Note that the distribution
of the information may not be uniform. In order to deal with more general
models, we do not use the Markov normalization given by Eq. (2.1) to define
the diffusion distance, instead, we use the diffusion kernel $k$ to define it.
More specifically, we work with the diffusion distance at time $t$ given by

$D^{t}(x,y):={\|\,k^{t}(x,\cdot)-k^{t}(y,\cdot)\|}_{L^{2}(X)}\mbox{ ,}$ (3.1)

where $k^{t}$ is the kernel of the operator $A^{t},$ and $A$ is the integral
operator defined as

$A(f)(x):=\int_{X}k(x,y)\,f(y)\,d\mu(y).$ (3.2)

###  3.1 Diffusion representation for t=1

We now design a representation for the diffusion distance given by Eq. (3.1),
where $k$ is an asymmetric kernel. Suppose that
$\\{W_{m_{1}}\\}_{m_{1}\in\mathbb{Z}}$ and
$\\{W_{m_{2}}\\}_{m_{2}\in\mathbb{Z}}$ are two orthonormal bases of
$L^{2}(X),$ and that there exists a positive constant $M$ such that for any
$m_{1}\in\mathbb{Z},$ and all $x\in X$,

$\|W_{m_{1}}(x)\|\leq M.$

We recall that in such case the tensor product $\\{W_{m_{1}}\otimes
W_{m_{2}}\\}_{(m_{1},m_{2})\in\mathbb{Z}\times\mathbb{Z}}$ defined by

$W_{m_{1}}\otimes W_{m_{2}}(x,y)=W_{m_{1}}(x)\,{W}_{m_{2}}(y),$ (3.3)

is an orthonormal basis of $L^{2}(X\times X)$, (for more details, see [12, 13,
14]). We note, in passing, that our approach also works if the $L^{2}(X)$
space is finite dimensional, in which case $\mathbb{Z}$ should be substituted
by a finite index set. To develop our theory we assume the following
hypothesis on the kernel $k$.

######  Hypothesis 3.1.

Suppose that $k\in L^{2}(X\times X),$ and for a.e $x\in X,$ the kernel
function $k(x,\cdot)$ belongs to the space $L^{2}(X)$.

If we assume the above hypothesis, we can use the basis given by Eq. (3.3) to
write the kernel $k(x,y)$ as

$k(x,y)=\sum_{(m_{1},m_{2})\in\mathbb{Z}\times\mathbb{Z}}a_{(m_{1},m_{2})}\,W_{m_{1}}\otimes
W_{m_{2}}(x,y),$ (3.4)

where $a_{(m_{1},m_{2})}$ are the coefficients. Using this decomposition we
obtain a representation form for the diffusion distance at time $t=1$.

######  Theorem 3.1 (Diffusion representation for $t=1$).

Assume that the kernel $k$ satisfies the Hypothesis 3.1, and that the
representation of $k$ in the coordinate system (3.3) is given by Eq. (3.4). If
the coefficients satisfy the summability condition

$\sum_{(m_{1},m_{2})\in\mathbb{Z}\times\mathbb{Z}}\|a(m_{1},m_{2})\|<\infty,$
(3.5)

then the diffusion distance at time $t=1$ has the representation form

$(D^{1}(x,y))^{2}=\sum_{m_{2}\in\mathbb{Z}}\|\sum_{m_{1}\in\mathbb{Z}}a(m_{1},m_{2})\,(W_{m_{1}}(x)-W_{m_{1}}(y))\|^{2}\mbox{
.}$ (3.6)

###### Proof.

Using the summability condition of Eq. (3.5), we can write the function
$k(x,\cdot)-k(y,\cdot)$ as

$\sum_{m_{2}\in\mathbb{Z}}\left(\sum_{m_{1}\in\mathbb{Z}}a(m_{1},m_{2})\,(W_{m_{1}}(x)-W_{m_{1}}(y))\right)W_{m_{2}}(\cdot).$

Since the set $\\{W_{m_{2}}\\}$ is an orthonormal basis, we conclude that

$(D^{1}(x,y))^{2}=\sum_{m_{2}\in\mathbb{Z}}\|\sum_{m_{1}\in\mathbb{Z}}a(m_{1},m_{2})\,(W_{m_{1}}(x)-W_{m_{1}}(y))\|^{2}.$

∎

For practical purposes, we do not use the representation formula given in
Theorem 3.1 to approximate the diffusion distance, because this representation
includes two sums with many terms. Instead, we use an approximation involving
sums of few terms, this is established in Theorem 3.2. In order to prove this
theorem, we first prove an auxiliary lemma.

######  Lemma 3.1 (Approximation lemma).

Consider the function

$f_{k_{1},k_{2}}:L^{1}(\mathbb{Z}\times\mathbb{Z})\times
L^{1}(\mathbb{Z}\times\mathbb{Z})\times X\times X\to\mathbb{R}_{\geq 0},$

defined by

$f_{k_{1},k_{2}}(a_{1},a_{2},x,y):=\sum_{\|m_{2}\|\leq
k_{2}}\|\sum_{\|m_{1}\|\leq
k_{1}}a_{1}(m_{1},m_{2})\,W_{m_{1}}(x)-a_{2}(m_{1},m_{2})\,W_{m_{1}}(y)\|^{2}\mbox{
.}$

Suppose that $a_{1}$ and $a_{2}$ are sequences in
$L^{1}(\mathbb{Z}\times\mathbb{Z})$, then for each $\delta>0,$ there exist
positive integers $\overline{k_{1}}$ and $\overline{k_{2}}$, such that if
$\overline{k_{1}}\leq k_{1},$ $\overline{k_{2}}\leq k_{2},$ and if $x,y\in X,$
the following inequality holds

$\|f_{k_{1},k_{2}}(a_{1},a_{2},x,y)-f(a_{1},a_{2},x,y)\|<\delta\mbox{ , }$

where $f(a_{1},a_{2},x,y)=\lim_{k\to\infty}f_{k,k}(a_{1},a_{2},x,y).$

###### Proof.

Let $b(m_{1},m_{2})$ be defined by

$b(m_{1},m_{2})=a_{1}(m_{1},m_{2})\,W_{m_{1}}(x)-a_{2}(m_{1},m_{2})W_{m_{1}}(y).$

Since the $L^{2}(\mathbb{Z}\times\mathbb{Z})$ norm is smaller than, or equal
to, the $L^{1}(\mathbb{Z}\times\mathbb{Z})$ norm, we obtain that the
expression $\|f_{k_{1},k_{2}}(a_{1},a_{2},x,y)-f(a_{1},a_{2},x,y)\|$ is
bounded from above by

$\left(\sum_{\|m_{2}\|\leq
k_{2}}\sum_{\|m_{1}\|>k_{1}}\|b(m_{1},m_{2})\|+\sum_{\|m_{2}\|>k_{2}}\sum_{m_{1}\in\mathbb{Z}}\|b(m_{1},m_{2})\|\right)^{2}.$

The above expression is dominated by

$\left(M\sum_{(m_{1},m_{2})\notin B(0,k_{1})\times
B(0,k_{2})}\|a_{1}(m_{1},m_{2})\|+\|a_{2}(m_{1},m_{2})\|\right)^{2}.$

Therefore, by Assumption (3.5), we conclude that for a given $\delta>0,$ we
can take $k_{1}$ and $k_{2}$ large, such that the above term is less than or
equal to $\delta.$ ∎

As a consequence of the above lemma, we prove the following theorem, which
states that we can approximate the diffusion distance with finite sums.

######  Theorem 3.2.

Assume that the kernel $k$ satisfies the same hypotheses of Theorem 3.1. Then,
for each $\delta>0,$ there exist positive integers $\overline{k_{1}}$ and
$\overline{k_{2}}$ such that for any $\overline{k_{1}}\leq k_{1},$ and
$\overline{k_{2}}\leq k_{2},$ and all $x,y\in X$, the following inequality
holds

$|f_{k_{1},k_{2}}(a_{1},a_{1},x,y)-(D^{1}(x,y))^{2}|\leq\delta\mbox{ .}$

###### Proof.

The proof is a consequence of Theorem 3.1 together with Lemma 3.1. ∎

###  3.2 Diffusion representation for arbitrary time

Suppose that $t$ is a positive integer denoting an arbitrary time. We now use
the coordinate system of Eq. (3.3) to find the representation form for the
kernel $k^{t+1}$ in terms of the coefficients $a(m_{1},m_{2})$. Let $k$ be an
asymmetric kernel, and suppose that for any $1\leq j\leq t,$ the kernel
$k^{j}$ satisfies all the hypotheses of Theorem 3.1. Under this assumption, we
can use the Fubini’s theorem to write recursively the kernel $k^{t+1}$ of the
operator $A^{t+1}$ as

$k^{t+1}(x,y)=<k(x,\cdot),k^{t}(\cdot,y)>_{L^{2}(X)}\mbox{ .}$ (3.7)

Assuming that the kernel $k^{t}$ has the series representation

$k^{t}(x,y)=\sum_{(m_{1},m_{2})\in\mathbb{Z}\times\mathbb{Z}}a_{(m_{1},m_{2})}^{t}\,W_{m_{1}}\otimes
W_{m_{2}}(x,y)\mbox{. }$

then by Eq. (3.7) we obtain that

$a_{(m_{1},m_{2})}^{t+1}=\sum_{k\in\mathbb{Z}}a_{(m_{1},k)}\,a_{(k,m_{2})}^{t}\mbox{
.}$ (3.8)

Recursively, we obtain that the expression for the coefficients of the kernel
$k^{t+1}$

$a_{(m_{1},m_{2})}^{t+1}=\sum_{n_{1}\in\mathbb{Z}}\sum_{n_{2}\in\mathbb{Z}}\sum_{n_{3}\in\mathbb{Z}}\dots\sum_{n_{t}\in\mathbb{Z}}a_{(m_{1},n_{1})}\,a_{(n_{1},n_{2})}a_{(n_{2},n_{3})}\dots
a_{(n_{t},m_{2})}\mbox{ .}$ (3.9)

Again, the above expression contains infinitely many sums. We now prove that
we can approximate the coefficients of $k^{t}$ using finite sums. The
following lemma establishes this result.

######  Lemma 3.2.

For any $\delta>0$, there exist $k_{0}$, such that for any $k_{0}\leq k$ we
have

$\|a^{t+1}-h^{t+1}_{k}(a)\|_{L^{1}(\mathbb{Z}\times\mathbb{Z})}\leq\delta\mbox{
,}$

where $h^{t+1}_{k}(a)(m_{1},m_{2})$ is the finite sum

$h^{t+1}_{k}(a)(m_{1},m_{2})=\sum_{\|n_{1}\|\leq k}\sum_{\|n_{2}\|\leq
k}\sum_{\|n_{3}\|\leq k}\dots\sum_{\|n_{t}\|\leq
k}a_{(m_{1},n_{1})}\,a_{(n_{1},n_{2})}a_{(n_{2},n_{3})}\dots
a_{(n_{t},m_{2})}\mbox{ . }$

###### Proof.

We prove by induction over $t$, for $t=1$ is clear since

$\displaystyle\sum_{(m_{1},m_{2})\in\mathbb{Z}\times\mathbb{Z}}\lvert\lvert
a_{(m_{1},m_{2})}^{2}-h^{2}_{k}(m_{1},m_{2})\rvert\rvert$
$\displaystyle\leq\sum_{\lvert\lvert n_{1}\rvert\rvert\geq
k}\sum_{m_{1}\in\mathbb{Z}}\sum_{m_{2}\in\mathbb{Z}}\lvert\lvert
a_{(m_{1},n_{1})}a_{(n_{1},m_{2})}\rvert\rvert$ $\displaystyle\leq\lvert\lvert
a\rvert\rvert_{L^{1}}\sum_{\lvert\lvert n_{1}\rvert\rvert\geq
k}\sum_{m_{1}\in\mathbb{Z}}\lvert\lvert a_{(m_{1},n_{1})}\rvert\rvert\mbox{
.}$

Then, by Assumption 3.5, we have that for $k$ large the above inequality is
less than or equal to $\delta.$ We now assume that the claim holds for $t$,
and we prove that also holds for $t+1$. For $k$ large we have that

$\|a^{t}-h^{t}_{k}(a)\|_{L^{1}(\mathbb{Z}\times\mathbb{Z})}\leq\delta.$

The above inequality implies

$\|h^{t}_{k}(a)\|_{L^{1}(\mathbb{Z}\times\mathbb{Z})}\leq(\delta+\|a\|_{L^{1}(\mathbb{Z}\times\mathbb{Z})})\mbox{
.}$

Furthermore, by Eq. (3.8), we have that

$\|a_{(m_{1},m_{2})}^{t+1}-h^{t+1}_{k}(a)(m_{1},m_{2})\|,$

is less or equal to

$\sum_{i\in
I}\|a_{(m_{1},i)}\|\,\|a_{(i,m_{2})}^{t}-h^{t}_{k}(a)(i,m_{2})\|+\sum_{\|i\|\geq
k}\|a_{(m_{1},i)}\|\|h^{t}_{k}(a)(i,m_{2})\|\mbox{ .}$

Using the above inequalities, we obtain the following estimate

$\|a^{t+1}-h^{t+1}_{k}(a)\|_{L^{1}}\leq\delta\|a\|_{L^{1}}+(\delta+\|a\|_{L^{1}})\sum_{\lvert\lvert
i\rvert\rvert\geq k}\sum_{m_{1}\in\mathbb{Z}}\lvert\lvert
a_{(m_{1},i)}\rvert\rvert\mbox{ . }$

Therefore, by Assumption (3.5), we conclude that the claim holds for $t+1$. ∎

We now use the above result to design a representation for the diffusion
distance at time $t+1$. This representation is based on finite sums of the
coefficients $a_{(m_{1},m_{2})}$. We establish this result in Theorem 3.3, the
proof involves the following technical lemma.

######  Lemma 3.3 (Continuity).

Consider the function $f_{k_{1},k_{2}}$ as in Lemma 3.1, and suppose that $a$
and $b$ are two sequences in $L^{1}(\mathbb{Z}\times\mathbb{Z}).$ Then, for
any positive number $\epsilon$, there exists a positive number $\delta$, such
that for any pair of sequences $c$ and $d$ satisfying

$\|a-c\|_{L^{1}(\mathbb{Z}\times\mathbb{Z})}\leq\delta\hskip
28.45274pt\mathrm{and}\hskip
28.45274pt\|b-d\|_{L^{1}(\mathbb{Z}\times\mathbb{Z})}\leq\delta\mbox{, }$

then, the following inequality holds for all $x,y\in X,$ and for all positive
integers $k_{1},k_{2},$

$|f_{k_{1},k_{2}}(a,b,x,y)-f_{k_{1},k_{2}}(c,d,x,y)|\leq\epsilon\mbox{. }$

###### Proof.

Let $k_{1},k_{2}$ be positive integers, and define the function $R_{m_{2}}$ as

$R_{m_{2}}(a_{1},a_{2},x,y):=\sum_{\|m_{1}\|\leq
k_{1}}a_{1}(m_{1},m_{2})\,W_{m_{1}}(x)-a_{2}(m_{1},m_{2})\,W_{m_{1}}(y).$

Observe that

$\|R_{m_{2}}(a,b,x,y)-R_{m_{2}}(c,d,x,y)\|\leq M\sum_{\|m_{1}\|\leq
k_{1}}(\|(a-c)(m_{1},m_{2})\|+\|(b-d)(m_{1},m_{2})\|)\mbox{ ,}$ (3.10)

where $M$ is a constant independent of $x$ and $y$. For any real numbers $B$
and $C$, the following inequality holds

$|B^{2}-C^{2}|=|(B-C)^{2}+2BC-2C^{2}|\leq(B-C)^{2}+2|C||B-C|\mbox{ .}$ (3.11)

Applying the above inequality to $B=\|R_{m_{2}}(c,d,x,y)\|,$ and
$C=\|R_{m_{2}}(a,b,x,y)\|$, together with the fact that the $L^{2}$ norm is
smaller than or equal to the $L^{1}$ norm, we obtain that
$|f_{k_{1},k_{2}}(a,b,x,y)-f_{k_{1},k_{2}}(c,d,x,y)|$ is bounded from above by

$M(E^{2}+2(\|a\|_{L^{1}}+\|b\|_{L^{1}})E)\mbox{, }$

where $M$ is a constant which does not depend on $k_{1},k_{2},x,y,$ and

$E:=\|a-c\|_{L^{1}}+\|b-d\|_{L^{1}}\mbox{ .}$

∎

######  Theorem 3.3.

Suppose that for any $1\leq j\leq t+1,$ the kernel $k^{j}$ satisfies the same
hypotheses of Theorem 3.1. Then, for any positive number $\epsilon$, there
exist positive integers $\overline{k_{1}}$, $\overline{k_{2}}$, and
$\overline{k_{3}},$ such that for any natural numbers $k_{1},k_{2},k_{3}$,
satisfying $\overline{k_{1}}\leq k_{1},$ $\overline{k_{2}}\leq k_{2}$ and
$\overline{k_{3}}\leq k_{3}$, the following inequality holds

$|f_{k_{1},k_{2}}(h^{t+1}_{k_{3}}(a),h^{t+1}_{k_{3}}(a),x,y)-(D^{t+1}(x,y))^{2}|\leq\epsilon\mbox{
,}$

where $f_{k_{1},k_{2}}$ is defined as in Lemma 3.1, and $h^{t+1}_{A}$ as in
Lemma 3.2.

###### Proof.

Applying Theorem 3.2 to the kernel $k^{t+1}$, we have that there exist
$\overline{k_{1}}$, $\overline{k_{2}}$, such that for any
$\overline{k_{1}}\leq k_{1},$ and $\overline{k_{2}}\leq k_{2},$ the following
inequality holds

$|f_{k_{1},k_{2}}(a^{t+1},a^{t+1},x,y)-(D^{t}(x,y))^{2}|\leq\epsilon/2\mbox{ ,
}$

for any $x,y\in X$. Moreover, using Lemmas 3.2 and 3.3, there exists a
positive integer $\overline{k_{3}}$ with the property that if
$\overline{k_{3}}\leq k_{3}$, then for all $x,y\in X,$

$|f_{k_{1},k_{2}}(a^{t+1},a^{t+1},x,y)-f_{k_{1},k_{2}}(h^{t+1}_{k_{3}}(a),h^{t+1}_{k_{3}}(a),x,y)|\leq\epsilon/2\mbox{
.}$

Thus, by the triangular inequality we obtain the desired result. ∎

###  3.3 Diffusion representation for changing data

In this section we use the coordinate system of Eq. (3.3) to represent the
dynamic diffusion distance induced by asymmetric kernels. Suppose that
$\\{k_{\gamma}\\}$ is a family of asymmetric kernels defined in the data set
$X$. Again, we consider the diffusion distance without the Markov
normalization, that is, we work with the dynamic diffusion distance given by

$D^{t}(x_{\gamma},y_{\beta})={\|\,k_{\gamma}^{t}(x,\cdot)-k_{\beta}^{t}(y,\cdot)\|}_{L^{2}(X)}.$

We assume that for each parameter $\gamma$, the kernel $k_{\gamma}$ satisfies
all the hypotheses of Theorem 3.1. In this case, we can write the kernel
$k_{\gamma}$ as

$k_{\gamma}(x,y)=\sum_{(m_{1},m_{2})\in\mathbb{Z}\times\mathbb{Z}}a^{\gamma}_{(m_{1},m_{2})}\,W_{m_{1}}\otimes
W_{m_{2}}(x,y).$

The following theorem gives a representation for the dynamic diffusion
distance. The proof of the theorem is similar to that of in Theorem 3.1.

######  Theorem 3.4.

Assume that for all $\gamma$, the kernel $k_{\gamma}$ satisfies all the
hypotheses of Theorem 3.1, then the dynamic diffusion distance
$(D^{t}(x_{\gamma},y_{\beta}))^{2}$ can be written as

$\sum_{m_{2}\in\mathbb{Z}}\|\sum_{m_{1}\in\mathbb{Z}}(a^{\gamma})^{t}(n,m)\,W_{m_{1}}(x)-(a^{\beta})^{t}(n,m)W_{m_{1}}(y)\|^{2}.$

######  Remark 3.1.

Using the same ideas of the proof of Theorem 3.3, we can prove that it is
possible to approximate the dynamic diffusion distance by sums involving few
terms of $a^{\gamma}$ and $a^{\beta}$ . To be more specific, the same
statement of Theorem 3.3 holds if we replace the classical diffusion distance
by the dynamic diffusion distance $D^{t}(x_{\gamma},y_{\beta})$, and the
function $f_{k_{1},k_{2}}(h^{t+1}_{C}(a),h^{t+1}_{C}(a),x,y)$ by the function
$f_{k_{1},k_{2}}(h^{t+1}_{C}(a^{\gamma}),h^{t+1}_{C}(a^{\beta}),x,y).$

We use the previous result to compute the global diffusion distance in terms
of the coefficients $a^{\gamma}$. We recall that the global diffusion distance
between $X_{\alpha},$ and $X_{\beta},$ is defined by

$(\mathbb{D}^{t}(X_{\gamma},X_{\beta}))^{2}=\int_{X}({D}^{t}(x_{\gamma},x_{\beta}))^{2}d\mu(x).$

######  Theorem 3.5.

Under the same assumptions of Theorem 3.4, the global diffusion distance at
time $t$ can be written as

$(\mathbb{D}^{t}(X_{\gamma},X_{\beta}))^{2}=\sum_{(n,m)\in\mathbb{Z}\times\mathbb{Z}}\|((a^{\gamma})^{t}-(a^{\beta})^{t})(n,m)\|^{2}\mbox{
.}$

###### Proof.

Theorem 3.4 implies that

$(\mathbb{D}^{t}(X_{\gamma},X_{\beta}))^{2}=\int_{X}\sum_{m\in\mathbb{Z}}\|\sum_{n\in\mathbb{Z}}((a^{\gamma})^{t}-(a^{\beta})^{t})(n,m)W_{n}(x)\|^{2}d\mu(x)\mbox{
.}$

Expanding the quadratic form, we obtain that

$\mathbb{D}^{t}(X_{\gamma},X_{\beta})^{2}=\int_{X}\sum_{m\in\mathbb{Z}}\sum_{(i,j)\in\mathbb{Z}\times\mathbb{Z}}((a^{\gamma})^{t}-(a^{\beta})^{t})(i,m)\,\overline{((a^{\gamma})^{t}-(a^{\beta})^{t})(j,m)}\,W_{i}(x)\,\overline{W_{j}}(x)d\mu(x)\mbox{.}$

Using Hölder’s inequality and the fact that $\|W_{n}\|_{L^{2}(X)}=1$, we
conclude that the expression

$S=\sum_{m\in\mathbb{Z}}\sum_{(i,j)\in\mathbb{Z}\times\mathbb{Z}}\int_{X}\|((a^{\gamma})^{t}-(a^{\beta})^{t})(i,m)\,\overline{((a^{\gamma})^{t}-(a^{\beta})^{t})(j,m)}\,W_{i}(x)\,\overline{W_{j}}(x)\|d\mu(x).$

is bounded from above by

$\sum_{m\in\mathbb{Z}}\left(\sum_{n\in\mathbb{Z}}\|((a^{\gamma})^{t}-(a^{\beta})^{t})(n,m)\|\right)^{2}\leq\left(\sum_{m\in\mathbb{Z}}\sum_{n\in\mathbb{Z}}\|((a^{\gamma})^{t}-(a^{\beta})^{t})(n,m)\|)\right)^{2}<\infty\mbox{,
}$

where we used the fact that the $L^{2}$ norm is smaller than or equal to the
$L^{1}$ norm. By the dominated convergence theorem, we can change the order
between the integral and the sums. Moreover, using the fact that $W_{n}$ is an
orthonormal basis of $L^{2}(X),$ we conclude that

$(\mathbb{D}^{t}(X_{\gamma},X_{\beta}))^{2}=\sum_{m\>\in\>\mathbb{Z}}\sum_{n\>\in\>\mathbb{Z}}\|((a^{\gamma})^{t}-(a^{\beta})^{t})(n,m)\|^{2}\mbox{
. }$

∎

###  3.4 Weak representation

The framework developed up to this point uses the absolute summability
condition (3.5). We now relax this assumption, and under an a priori
assumption on the diffusion distance, we design a representation for the
diffusion distance of points lying in a set of large measure. Without loss of
generality, we work with the changing data framework.

######  Theorem 3.6.

Let $t$ be a nonnegative integer, assume that the kernel $k^{t+1}_{\gamma}$
satisfies the Hypothesis 3.1, and also that for any integer $j$, with $1\leq
j\leq t,$ the kernel $k^{j}_{\gamma}$ satisfies all the hypotheses of Theorem
3.1. If we suppose that the dynamic diffusion distance at time $t+1$ is
bounded from above, that is, there exists a positive constant $C$ such that
for any $x,y\in X$, and all indices $\gamma,\beta,$

$D^{t+1}(x_{\gamma},y_{\beta})\leq C.$

Then, for any positive real number $\epsilon$, there exists a positive integer
$k_{0},$ such that for any integer $k\geq k_{0}$, there exists a measurable
set $E_{k}$, with the property that the Lebesgue measure of $X\setminus E_{k}$
satisfies $\mathcal{L}(X\setminus E_{k})<\epsilon$, and such that for any
$x,y\in E_{k}$ the following inequality holds

$\|f_{k,k}(h^{t+1}_{k}(a^{\gamma}_{k}),h^{t+1}_{k}(a^{\beta}_{k}),x,y)-(D^{t+1}(x_{\gamma},y_{\beta}))^{2}\|\leq\overline{C}(\sqrt{\epsilon}+\epsilon)\mbox{,
}$

where $\overline{C}$ is a constant and the functions $f_{k,k},$ $h^{t+1}_{k}$
are defined in Lemmas 3.1, and 3.2, and $a^{\gamma}_{k}$ is the truncated
sequence defined by

$a^{\gamma}_{k}(n,m)=\left\\{\begin{array}[]{@{}l@{\thinspace}l}\
a^{\gamma}(n,m)\mbox{, }&\quad\text{if}\quad\|(n,m)\|\leq k\\\ \ 0\mbox{,
}&\quad\text{otherwise}\mbox{ }\\\ \end{array}\right.\mbox{. }$

###### Proof.

Define the truncated kernel $k_{\gamma,k}$ by

$k_{\gamma,k}(x,y)=\sum_{(m_{1},m_{2})\in\mathbb{Z}\times\mathbb{Z}}a^{\gamma}_{k}{(m_{1},m_{2})}\,W_{m_{1}}\otimes
W_{m_{2}}(x,y)\mbox{ . }$

Using the recursive formula in Eq. (3.7), and the fact that
$\lim_{k\to\infty}k_{\gamma,k}=k_{\gamma}$, we see inductively that
$\lim_{k\to\infty}k^{t+1}_{\gamma,k}=k^{t+1}_{\gamma}$ for any non negative
integer $t$, where the convergence is in the $L^{2}(X\times X)$ norm.
Therefore, for any positive number $\epsilon$, there exists $k_{0}$ such that
for any integer $k$ satisfying $k_{0}\leq k,$ we have
$\|k^{t+1}_{\gamma}-k^{t+1}_{\gamma,k}\|_{L^{2}(X\times X)}\leq\epsilon.$
Define the set

$E_{k}=\\{x\in
X|\,\,\,\|k^{t+1}_{\gamma,k}(x,\cdot)-k^{t+1}_{\gamma}(x,\cdot),\|_{L^{2}(X)}\leq\sqrt{\epsilon}\\}\mbox{
.}$

Then, by Chebyshev’s inequality and Fubini’s theorem we obtain that

$\mathcal{L}(X\setminus
E_{k})\leq\frac{1}{\epsilon}\int_{X}\|k^{t+1}_{\gamma,k}(x,\cdot)-k^{t+1}_{\gamma}(x,\cdot),\|^{2}_{L^{2}(X)}dx\leq\epsilon\mbox{
.}$

On the other hand, if $x,y\in E_{k}$, then by Minkowski inequality

$|D^{t+1}(x_{\gamma,k},y_{\beta,k})-D^{t+1}(x_{\gamma},y_{\beta})|\leq
2\sqrt{\epsilon}\mbox{ , }$

where $x_{\gamma,k}$ is the data point $x$ endowed with the kernel structure
$k^{t+1}_{\gamma,k}$. We apply Inequality (3.11) of Lemma 3.3 with
$B=D^{t+1}(x_{\gamma,k},y_{\beta,k}),$ and $C=D^{t+1}(x_{\gamma},y_{\beta})$,
to obtain

$|D^{t+1}(x_{\gamma,k},y_{\beta,k})^{2}-D^{t+1}(x_{\gamma},y_{\beta})^{2}|\leq(1+2M)(\epsilon+\sqrt{\epsilon})\mbox{
.}$

Since the set of non-zero coefficients of the kernel $k_{\gamma,k}$ is finite,
then, by Eq. (3.9), we conclude that $a^{t+1}_{\gamma,k}(n,m)=0$ whenever
$\|(n,m)\|>k,$ and also that for any integers $n,m$ we have
$a^{t+1}_{\gamma,k}(n,m)=h^{t+1}_{k}(a)(n,m)$. Using these facts we conclude
that

$D^{t+1}(x_{\gamma,k},y_{\beta,k})^{2}=f_{k,k}(h^{t+1}_{k}(a^{\gamma}_{k}),h^{t+1}_{k}(a^{\beta}_{k}),x,y)\mbox{
.}$

∎

##  4 Applications, Experiments, and Results

###  4.1 Dimensionality reduction

In this section, we use our representation framework to reduce the dimension
of data sets, in such way that the Euclidean norm of the reduced data set
approximates the diffusion distance. Assume that the data set $X$ is endowed
with the kernel structure $k_{\alpha}.$ Then, by Remark 3.1 the map
$\phi_{k_{1},k_{2}}:X\to\mathbb{C}^{2k_{2}+1}$ defined by

$\phi_{k_{1},k_{2}}(x)(i)=\sum_{\|m_{1}\|\leq
k_{1}}(a^{\gamma})^{t}(m_{1},i-k_{2}-1)\,W_{m_{1}}(x),$ (4.1)

approximates the diffusion distance at time $t$. We summarize the above in
Algorithm 1.

In order to work with a small parameter $k_{2}$, and thus embed our data set
in a low dimensional Euclidean space, we need to use a proper orthonormal
basis of $L^{2}(X)$. Proper bases are those for which the information of the
kernel $k$ is concentrated on the coefficients $a_{(m_{1},m_{2})},$ whose pair
of integers $(m_{1},m_{2})$ are near to the origin.

Algorithm 1 Dimensionality reduction algorithm.

  1. 1.

Take the data in the form of an $M\times M$ (possibly asymmetric) matrix $L$
representing the kernel structure on the data set

  2. 2.

Run the Markov process $t$-times, and use the matrix $L^{t}$ instead of $L$.

  3. 3.

Compute the first components of the matrix $L^{t}$ in a proper orthonormal
tensor basis.

  4. 4.

Embed the data set using the function $\phi_{k_{1},k_{2}}$ of Eq. (4.1)

We note that depending on the data and the mathematical model, we may want to
use a proper normalization for the matrix $L$ in Algorithm 1.

In practical situations, our measure space is a finite data set
$X=(x_{i})^{n}_{i=1}$, endowed with the counting measure. If the kernel matrix
A is symmetric, then one can use the eigenvector basis. More generally, if the
matrix A is asymmetric, one can use the singular value decomposition (SVD) to
write

$\textbf{A}=\sum_{i=1}^{n}c_{i}\,L_{i}\otimes R_{i}\mbox{ , }$

where $L_{i}$ and $R_{i}$ are the left singular vectors and right singular
vectors of A, respectively, and where $c_{i}$ are the singular values of A.
The computational complexity of the SVD of an $n\times n$ matrix is
$O(n^{3})$. In the following experiments, we use the Fourier basis to improve
the computational complexity. In fact, the computational complexity of the 2-D
FFT is $O(n^{2}\log n).$ We recall that the Fourier basis
$(E_{k})^{n-1}_{k=0}$ in the complex vector space $\mathbb{C}^{n}$ is given by

$E_{k}=\frac{1}{\sqrt{n}}(1,e^{2\pi ik/n},e^{2\pi i2k/n},...,e^{2\pi
i(n-1)k/n})\mbox{ .}$ (4.2)

In all experiments we compute the coefficients using the 2D-FFT algorithm. All
experiments were run in Matlab software, using a desktop computer with the
following configuration: Intel core i7-2600 3.4 GHz processor, and 16 GB RAM.

###  4.2 Synthetic data using a symmetric kernel

In this experiment, our data set $X$ consists of $n$ random points in the
sphere $S^{2}$ (Figure 4.1). Here, we use the parametrization

$\displaystyle x(u,v)$ $\displaystyle{}={}$ $\displaystyle\cos{u}$
$\displaystyle\sin{v}\mbox{, }$ $\displaystyle y(u,v)$ $\displaystyle{}={}$
$\displaystyle\sin{u}$ $\displaystyle\sin{v}\mbox{, }$ $\displaystyle z(u,v)$
$\displaystyle{}={}$ $\displaystyle\cos{v}\mbox{, }$

for $0\leq u\leq\pi$ and $0\leq v\leq 2\pi$. We endow the data set $X$ with
the Markov normalization $\rho$ defined in Eq. (2.1), of the Gaussian kernel

$k(x,y)=e^{-||x-y||^{2}}\mbox{ .}$

Figure 4.1: Data set $X$ with $512$ random points in the sphere $S^{2}$.

Our goal is to compare the efficiency between the representation given in Eq.
(2.4) using the eigenvector basis and the representation of Theorem 3.1, using
the Fourier basis of Eq. (4.2). In Figure 4.2, we show the first two
coordinates for each representation with a data set of 512 points. Observe
that the first two coordinates of both representations are similar.

(a)  Eigenvector basis

(b) Real and imaginary part of the Fourier basis

Figure 4.2: Plot of the two dimensional embedding for the data set $X$ using
the eigenvector basis coefficients (a), and the Fourier basis coefficients
(b). Note the scale.

In Figure 4.3, we plot the error, and computational time (in seconds) of the
first two coordinates for several values of $n$. We also remark that by using
the Fourier basis, the performance of the representation is faster than using
the eigenvector basis, and also that the Fourier basis gives an acceptable
error when compared to the eigenvector method.

(a) Error of the approximation

(b) Computational time

Figure 4.3: Plot of the $L^{2}$ error and computational time of the first two
coordinates for different $n\times n$ kernel-sizes, for the data set of random
points in the sphere.

###  4.3 Synthetic data on the Möbius strip

Here, we assume that our data set $X$ is a set of $300$ data points
distributed along of the Möbius strip (Figure 4.4). We use the parametrization

$\displaystyle x(u,v)$
$\displaystyle{}={}(1+\frac{v}{2}\cos{\frac{u}{2}})\cos{u}\mbox{, }$
$\displaystyle y(u,v)$
$\displaystyle{}={}(1+\frac{v}{2}\cos{\frac{u}{2}})\sin{u}\mbox{, }$
$\displaystyle z(u,v)$ $\displaystyle{}={}\frac{v}{2}\sin{\frac{u}{2}}\mbox{,
}$

for $0\leq u\leq 2\pi$ and $-\frac{1}{2}\leq v\leq\frac{1}{2}$. We endow the
data set $X$ with the Markov normalization $\rho$ defined in Eq. (2.1), of the
weight Gaussian kernel

Figure 4.4: Data set $X$ with $300$ random points in the Möbius strip $M^{2}$.

$k(x,y)=(S(x-y)+1)e^{-||x-y||^{2}}\mbox{, }$

where $S(z)$ is the sign function of the angle (in cylindrical coordinates) of
the vector $z$. This kernel measures local information taking into account if
the first two components are rotating clockwise. In this experiment, we
compared the performance of the representation using the SVD, and the
representation using the Fourier basis. In Figure 4.5, we plot the first two
coordinates for each representation. Note that the real part of the
representation given by the Fourier basis allows us to see in more detail the
distribution of the data set. In fact, the representation that uses the
Fourier basis recognizes the rotation of the data set. However, the
presentation that uses SVD does not allow recognizing this feature of $X$.
This is due to the fact that the representation using the SVD approximates the
kernel $k^{T}k$, instead of the kernel $k$ (Figure 4.6). Therefore, the
representation using the SVD does not distinguish some geometric properties of
the data set $X$.

(a)  svd

(b) Real and imaginary part of the Fourier basis

Figure 4.5: Dimensionality reduction using the Eigenvector basis, and Fourier
basis

(a)  Kernel matrix $k$

(b)  Kernel matrix $k^{T}k$

Figure 4.6: Kernel matrix $k$ (a), and normalization $k^{T}k$ (b)

###  4.4 Synthetic data using an asymmetric kernel

Here, we assume that our data set $X$ is a random set of $256$ data points. We
endow this data set with the kernel structure $k$ given by the Tom Jobim
picture of Figure 4.7 whose dimensions are $256\times 256$ pixels. That is,
$k(x,y)$ is defined to be the gray scale value of the pixel coordinates
$(x,y)$. As in the previous experiment, we use the Markov normalization of the
kernel $k$. In this experiment, we compared the performance of the
representation using the SVD and the representation using the Fourier basis.

We stress that our objective with this example is not to try to do image
processing, but rather to use a picture so that we are able to visually assess
the quality of the approximation. This point will be furthered in the sequel.

In Figure 4.8, we plot the approximation of the kernel using the SVD and the
Fourier basis, both using the parameters $k_{1}=256$ and $k_{2}=64$. One may
notice that we see a horizontal modulation both under the SVD and the Fourier
basis methods. However, this modulation is stronger in the Fourier method.
This is due to the fact that we have only used high frequencies to approximate
this image.

In this case we observe that despite using a smaller number of parameters, it
is possible to obtain a good approximation of the original image. In Figure
4.9, we plot the $L^{2}$ error, and computational time (in seconds), in a
logarithmic scale, of the embedding data set using several approximation
orders. As in the previous experiment, we see that using the Fourier basis,
the performance is faster, and provides an acceptable error when compared to
the SVD method. Furthermore, we point out that in the two previous experiments
we did not obtain a better performance with respect to computational time when
we used the truncated SVD instead of the SVD.

Figure 4.7: Synthetic asymmetric kernel structure chosen from a picture of Tom
Jobim in Ipanema Beach (Rio de Janeiro, Brazil). Source [15].

(a) Approximation using the SVD

(b) Approximation using the Fourier basis

Figure 4.8: Plot of the approximation using the SVD and the Fourier basis,
both using the parameters $k_{1}=256$ and $k_{2}=64$

(a) Approximation error

(b) Computational time

Figure 4.9: Plot of the $L^{2}$ error and computational time for the data set
representing the picture of Tom Jobim.

###  4.5 Temperature changes in Brazil

In the last decades the world temperature distribution has presented drastic
changes, in part this is likely due to human activities [16, 17, 18]. Here, we
use the diffusion distance for changing data to detect the regions of Brazil
in which the local temperature variation was the highest in the years 2010 and
2018 (Figure 4.10), compared with the year 2000 (Figure 4.10). In fact, if a
certain region has a great diffusion distance, it means that this particular
region has presented significant changes in its temperature. This experiment
is based on the change detection on hyperspectral imagery data proposed in
Ref. [9]. Our data set consists of $N=13,974$ points, and each point
represents a pixel coordinate of the Brazilian map. These points are a subset
of a picture of size $170\times 170$ pixels. Here we do not take into account
the blank pixels, which correspond to places outside the Brazilian soil. For
each year we endow the data set with the un-normalized kernel $K$ which is
defined on $X\,\times\,X$ by

$K(x_{i,j},x_{k,l})=\frac{1}{\sqrt{N}}T_{\alpha}(x_{k,l})\,e^{-\|(i,j)-(k,l)\|^{2}/2\sigma^{2}}\mbox{
, }$

where $T_{\alpha}(x_{k,l})$ is the temperature in the rectangular pixel
$x_{k,l}$ in the year $\alpha\in$ {2000, 2010, 2018}, and where $\|\cdot\|$ is
the Euclidean distance. In this experiment we use the scaling parameter
$2\sigma^{2}=650$. We obtained similar results with a parameter in a range of
$600\leq 2\sigma^{2}\leq 700$. We use this kernel without normalization in
order to avoid its high computational cost. This data set was taken from the
Brazilian National Institute of Meteorology website [19]. This asymmetric
kernel represents the distribution of the local temperature around the
rectangular pixel $x_{k,l}$.

We use Theorem 3.4 to approximate the dynamic diffusion distance
$(D^{1}(x_{2000},x_{\gamma}))^{2}$, for $\gamma\in\\{2010,2018\\}$. Due to the
high dimensionality of the kernel matrix, the SVD algorithm did not run in the
computer whose configuration is given in Section 4.1. Therefore, we cannot use
the singular vector basis to represent the diffusion distance. Here, we use
the Fourier basis defined in Eq. $(\ref{fourier})$ to represent this diffusion
distance. We approximate the dynamic diffusion distance using the parameters
$k_{1}=13,974$, and $k_{2}\in\\{5,100\\}$. See Figures 4.11 and 4.13 for
$2010$ and $2018$, respectively. The green-yellow scale represents the
intensity of the dynamic diffusion distance, in which the yellow regions have
a greater diffusion distance as compared to the green regions. To detect which
regions have the greatest positive variation, that is, regions in which the
temperature has increased, we use contour plots of the diffusion distance,
taking into account regions where the temperature increased. See Figures 4.12
and 4.14 for 2010 and 2018, respectively.

In Table 1, we show the global diffusion distance for each year. Observe that
the distance is greater in 2018 than in 2010. This suggests that during 2018,
there were more changes in temperature when compared to 2010. In Figures 4.15
and 4.16, we plot the error and computational time of the performance for
several approximation orders $k_{2}$, for 2010 and 2018, respectively. We
evaluate the performance of the orders using the metric

(a) Year $2000$

(b) Year $2010$

(c) Year $2018$

Figure 4.10: Plot of the average Brazilian temperature distribution in years
$2000$, $2010$, $2018$.

$M_{B}=\frac{1}{N}E\times t_{cpu}\mbox{ ,}$ (4.3)

where $E$ the absolute error between $(D^{1}(x_{2000},x_{\gamma}))^{2}$ and
its approximation, and $t_{cpu}$ is the computational time (in seconds) to
compute the approximation. We see that even using small orders, it is possible
to obtain a good performance compared to larger orders.

Year | Global diffusion distance $(\;^{\circ}\mathrm{C})$
---|---
2010 | 25.5506
2018 | 37.5000
Table 1: Global diffusion distance in the years 2010 and 2018, with respect to
year 2000.

(a) $k_{2}=5$

(b) $k_{2}=100$

(c)  $D^{1}(x_{2000},x_{2010})^{2}$

Figure 4.11: Plot of the approximations for the dynamic diffusion distance
$(D^{1}(x_{2000},x_{2010}))^{2}$ with different values of $k_{2}$ (a), (b),
and plot of the dynamic diffusion distance (c).

(a) $k_{2}=5$

(b) $k_{2}=100$

(c)  $D^{1}(x_{2000},x_{2010})^{2}$

Figure 4.12: Contour plot of the representation approximation for the
diffusion distance $(D^{1}(x_{2000},x_{2010}))^{2}$ taking into account the
temperature increase.

(a) $k_{2}=5$

(b) $k_{2}=100$

(c)  $D^{1}(x_{2000},x_{2010})^{2}$

Figure 4.13: Plot of the approximations for the dynamic diffusion distance
$(D^{1}(x_{2000},x_{2018}))^{2}$ with different values of $k_{2}$ (a), (b),
and plot of the dynamic diffusion distance (c).

(a) $k_{2}=5$

(b) $k_{2}=100$

(c)  $D^{1}(x_{2000},x_{2010})^{2}$

Figure 4.14: Contour plot of the representation approximation for the
diffusion distance $(D^{1}(x_{2000},x_{2018}))^{2}$ taking into account the
temperature increase.

(a) Error of the approximation

(b) Computational time

(c) Metric $M_{B}$

Figure 4.15: Plot of the absolute error, computational time, and the metric
given by Eq. (4.3) with different orders to approximate
$(D^{1}(x_{2000},x_{2010}))^{2}$, for the data set with the temperature
distribution in Brazil.

(a) Error of the approximation

(b) Computational time

(c) Metric $M_{B}$

Figure 4.16: Plot of the absolute error, computational time, and the metric
given by Eq. (4.3) with different orders to approximate
$(D^{1}(x_{2000},x_{2018}))^{2}$, for the data set with the temperature
distribution in Brazil.

##  5 Conclusions

In this paper, we treat the problem of dimensionality reduction of data sets
whose structure is given by asymmetric kernels. Our methodology generalizes
the diffusion-map framework to asymmetric kernels, and computes a diffusion
representation based on the kernel coordinates in a proper orthonormal basis.
Our representation depends on two parameters, the first parameter defines the
approximation error and the second one the dimensionality.

In our experiments, we used the Fourier basis to represent the structure of
the data set. This choice is based on the fact that the Fourier basis
diagonalizes the Laplacian operator which is the main example of a diffusive
process. From the numerical viewpoint, the main advantage of using the Fourier
basis is that the 2d-FFT allows us a reduction from linear growth to
logarithmic growth of one of the factors. The latter contributes to the
computational complexity reduction when compared to traditional eigenvalue
methods. In fact, if we consider that our kernel is represented by an $n\times
n$ matrix, then the SVD takes $O(n^{3})$ of operations to be performed,
whereas the 2d-FFT decomposition is $O(n^{2}\log(n))$. This fact was confirmed
in a set of experiments with randomly generated kernels. Observe that the SVD
representation gives a better approximation, i.e., smaller errors. However, if
we use the Fourier basis we can obtain a good approximation of the data set
for a much lower computational cost. Additionally, the use of the Fourier
basis allows to see in more detail some geometric properties of the data set.
This suggests that it is possible to use the Fourier basis as an alternative
to the classic representation by eigenvalues, especially in computers with low
performance.

We perform a few applied experiments to test the theory. In particular, we
apply it to identify which regions of Brazil have presented a greater
variation in the temperature vis a vis other ones during the last decades. In
this experiment, we see that the Amazon region has presented more variations
in its temperature as compared to other places. This observation indicates
that further studies should be performed to investigate the possible reasons
for such variations. In order to avoid increasing the computational cost, we
did not use the Markov normalization in this kernel. Due to the high
dimensionality of the kernel matrix, the SVD algorithm did not run in our
computer for this experiment. However, we managed to execute the algorithm
using the 2d-FFT.

We performed also some experiments with synthetic data using a wavelet basis.
However, we did not obtain an improvement in the computation time or in the
error of the approximation when compared to the Fourier basis and the singular
vector basis.

Asymmetric kernels are present in a number of mathematical models, for
instance in weighted directed graphs. In such graphs, the transition from one
node to another is measured by an asymmetric kernel. In general, asymmetric
kernels are useful to represent a gain or loss of information when we move
from one point to another. Weighted directed graphs are thus used to model
real world problems such as, the traffic in a city, electrical network
systems, water flow in hydrological basins, and commodity trading between
economies. These are natural follow-up avenues to the present work. Another
natural follow-up would be the use of other orthonormal bases. One possibility
would be the diagonalizers of the Laplace-Beltrami operator on certain
manifolds. An example of such bases is given by spherical harmonics.

## Acknowledgements

The authors acknowledge the financial support provided by CAPES, Coordenação
de Aperfeiçoamento de Pessoal de Nível Superior (Finance code 001), CNPq,
Conselho Nacional de Desenvolvimento Científico e Tecnológico, and FAPERJ,
Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro.

## References

## References

  * [1] Gerard L.G. Sleijpen, Peter Sonneveld, and Martin B. van Gijzen.  Bi-cgstab as an induced dimension reduction method.  Applied Numerical Mathematics, 60(11):1100–1114, 2010.  Special Issue: 9th IMACS International Symposium on Iterative Methods in Scientific Computing (IISIMSC 2008).
  * [2] Gang Wu and Fei Li.  A randomized exponential canonical correlation analysis method for data analysis and dimensionality reduction.  Applied Numerical Mathematics, 164:101–124, 2021.  Special Issue on The Seventh International Conference on Numerical Algebra and Scientific Computing.
  * [3] K. Ch. Das.  The Laplacian spectrum of a graph.  Computers $\&$ Mathematics with Applications, 48(5):715 – 724, 2004\.
  * [4] M. Belkin and P. Niyogi.  Laplacian eigenmaps for dimensionality reduction and data representation.  Neural Computation, 15(6):1373–1396, June 2003.
  * [5] R. R. Coifman and S. Lafon.  Diffusion maps.  Applied and Computational Harmonic Analysis, 21(1):5 – 30, 2006\.  Special Issue: Diffusion Maps and Wavelets.
  * [6] M.A Hajji, S Melkonian, and R Vaillancourt.  Representation of differential operators in wavelet basis.  Computers $\&$ Mathematics with Applications, 47(6):1011 – 1033, 2004.
  * [7] M. Salhov, A. Bermanis, G. Wolf, and A. Averbuch.  Diffusion representations.  Applied and Computational Harmonic Analysis, 45(2):324 – 340, 2018\.
  * [8] M. Salhov, A. Bermanis, G. Wolf, and A. Averbuch.  Approximately-isometric diffusion maps.  Applied and Computational Harmonic Analysis, 38(3):399 – 419, 2015\.
  * [9] R. R. Coifman and M. J. Hirn.  Diffusion maps for changing data.  Applied and Computational Harmonic Analysis, 36(1):79 – 107, 2014\.
  * [10] R. R. Coifman, S. Lafon, A. B. Lee, M. Maggioni, B. Nadler, F. Warner, and S. W. Zucker.  Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps.  Proceedings of the National Academy of Sciences, 102(21):7426–7431, 2005.
  * [11] N. F. Marshall and M. J. Hirn.  Time coupled diffusion maps.  Applied and Computational Harmonic Analysis, 45(3):709 – 728, 2018\.
  * [12] M. Pedersen.  Functional Analysis in Applied Mathematics and Engineering.  CRC Press, 1999.
  * [13] L. Grafakos.  Classical Fourier Analysis.  Graduate Texts in Mathematics. Springer New York, 2014.
  * [14] E. M. Stein.  Singular Integrals and Differentiability Properties of Functions.  Princeton University Press, 1970.
  * [15] Diario de cultura.  Tom Jobim in Ipanema Beach picture.  https://www.diariodecultura.com.ar/columnas/crucigrama-antonio-brasileiro/.  Accessed: 2020-06-01.
  * [16] J. Hansen, R. Ruedy, M. Sato, and K. Lo.  Global surface temperature change.  Reviews of Geophysics, 48(4), 2010.
  * [17] Earth observatory.  World of Change: Global Temperatures.  https://earthobservatory.nasa.gov/world-of-change/DecadalTemp.  Accessed: 2020-06-01.
  * [18] United Nations Framework Convention on Climate Change.  The Paris Agreement.  https://unfccc.int/process-and-meetings/the-paris-agreement/d2hhdC1pcy.  Accessed: 2020-06-01.
  * [19] INMET.  Brazilian temperature dataset.  http://www.inmet.gov.br/portal/.  Accessed: 2020-06-01.