Latent common manifold learning with alternating
diffusion: analysis and applications

arXiv:1602.00078v2 [physics.data-an] 3 Aug 2017

Ronen Talmona , Hau-Tieng Wub,c
a Viterbi

Faculty of Electrical Engineering, Technion - Israel Institute of Technology, Haifa, Israel
of Mathematics and Department of Statistical Science, Duke University, Durham, NC, USA
c Mathematics Division, National Center for Theoretical Sciences, Taipei, Taiwan

b Department

Abstract
The analysis of data sets arising from multiple sensors has drawn significant research attention over the years. Traditional methods, including kernel-based methods,
are typically incapable of capturing nonlinear geometric structures. We introduce a latent common manifold model underlying multiple sensor observations for the purpose
of multimodal data fusion. A method based on alternating diffusion is presented and
analyzed; we provide theoretical analysis of the method under the latent common manifold model. To exemplify the power of the proposed framework, experimental results
in several applications are reported.
Keywords: common manifold, alternating diffusion, sensor fusion, multimodal
sensor, data fusion, seasonality, diffusion maps, manifold learning

1. Introduction
One of the long-standing challenges in signal processing is the fusion of information acquired by multiple, multimodal sensors. The problem of information fusion has
become particularly central in the wake of recent technological advances, which have
led to extensive collection and storage of multimodal data. Nowadays, many devices
and systems, e.g., cell-phones, laptops, and wearable-devices, incorporate more than
one sensor, often of different types. Of particular interest in the context of this paper
are the massive data sets of medical recordings and healthcare-related information, acquired routinely, for example, in operation rooms, intensive care units, and clinics. The
availability of such distinct and complementary information calls for the development
of new theories and methods, leveraging it toward achieving concrete data analysis
objectives, such as filtering and prediction, in a broad range of fields.
Problems in multimodal signal processing has been studied for many years and
has been approached from various research directions [27]. A classic approach for
such problems is Canonical Correlation Analysis (CCA) [22], which recovers highly

Email addresses: ronen@ef.technion.ac.il (Ronen Talmon), hauwu@math.duke.edu
(Hau-Tieng Wu)

Preprint submitted to Elsevier

August 4, 2017

correlated linear projections from two data sets. To extend the linear setting and to
address aspects of nonlinearities, CCA was applied in a kernel space (e.g., [28, 1]).
Recently, ample work based on the optimization criterion of CCA and kernels has
been presented, addressing multi-view problems, and in particular, using multi-kernel
learning (e.g. [2, 29]) and a variety of manipulations and combinations of kernels, e.g.
[5, 11, 10, 35, 36, 51, 25, 3, 23, 50, 48, 33, 34, 17, 38].
Our exposition begins by addressing a particular baseline problem. Consider multiple sensors measuring the same physical phenomenon, such that the properties of the
physical phenomenon are manifested as a hidden manifold (which we would like to extract), while each sensor presents its own deformation and has its own sensor-specific
effects (hidden nuisance variables, which we would like to suppress). We assume that
the relations between the measurements and the nuisance variables are unknown. The
goal is to uncover the common latent manifold and to suppress the sensor-specific variables, thereby extracting the essence of the data and separating the relevant information
from the irrelevant information.
This baseline problem highlights an important aspect in the analysis of multimodal
data sets; that is, the sensor-related variables may not be strictly related to noise and
interferences. Often, such variables exhibit “structures”, such as the position and orientation of the acquiring sensor, environmental effects, and channel characteristics.
To address this, we propose an approach based on manifold learning. The power of
manifold learning can be exploited in this setting, since it is designed to capture nonlinear topological and geometric structures underlying data, and it does not require
prior model knowledge, which can be particularly hard to obtain in the case of multiple modalities. Manifold learning is a class of nonlinear data-driven methods, e.g.
[47, 42, 14, 4], often used to extract the underlying structures in a given data set. Of
particular interest in the context of this paper is diffusion map (DM), [6], in which
discrete diffusion processes are constructed on the given data points; these diffusion
processes are designed to capture the geometry of the underlying variability in a single
data set. Multimodal data present a challenge to such a geometric analysis approach,
since multiple sensors often lead to undesired geometric structures stemming from the
diversity of the different acquisition techniques used in the sensors, making it more
difficult to identify and extract only the “important” variables. Nevertheless, multiple
data sets from various sensors encompass more information, and therefore, enable us
to recover a more reliable description of the measured (physical) phenomenon. Based
on manifold learning, several methods have been proposed to analyze simultaneously
multiple data sets. One approach is to concatenate the vectors representing each data
set into one vector [8]; however, the question of how each data set should be scaled and
concatenated naturally arises, especially if the data sets are acquired by very different
modalities. To address such scaling aspects, it has been proposed in [24] to use DM
to obtain a low-dimensional “standardized” representation of each data set, and then to
concatenate these representations. However, such methods aggregate all the variables
from all the given data sets, and they neither distinguish the important information nor
discard the sensor-specific variables.
Our research direction involving geometric analysis encompasses several significant advantages. First, the method we present is data-driven and “model-free” in the
sense that in addition to the manifold assumption, it does not rely on prior knowledge.
2

In multimodal problems, this is an important advantage, since it circumvents the need
to design an appropriate model for each modality, as well as the “hard wiring” required
for the fusion of different data sets. Second, manifold learning methods are typically
formulated in general settings, and therefore, do not require strong assumptions on the
nature of the data or on the nature of the sensors. As a result, our method is restricted
to neither certain applications nor to multi-view problems, consisting of data acquired
only by a single type of sensors. Third, the combination of geometric analysis, enabling to integrate subtle patterns and structures underlying data, and the availability of
multiple data sets providing complementary information, gives rise to the discovery of
intrinsic structures. Furthermore, the diffusion approach has been shown to be robust
to noise [15, 16], and hence allows us to devise a reliable approach for extracting the
interesting information from highly noisy data.
Recently, a data-driven method to recover the common latent variable underlying
multiple, multimodal sensor data based on alternating products of diffusion operators
was presented [30, 31]; we refer to this method as alternating diffusion (AD). It was
shown both theoretically and in illustrative examples (with real recordings) that this
method extracts useful information about the common source of variability in a multisensor experiment as if it were a typically diffusion operator applied directly to data
sampled from only the interesting/common source of variability. The formulation and
analysis in [30] are based on a setting including only metric spaces and the common
source of variability is identical in the different sensors.
In the current work, we extend [30] and enhance the theoretical results by explicitly
introducing a setting with a common manifold, which can be well accessed by different
sensors, while different types of deformations are introduced due to various effects
of the specific sensors. We focus on the case in which there is a common manifold
underlying all the given data sets, and propose a data-driven method for recovering
the common manifold and constructing its representation. More specifically, in an
analogous way to the classic diffusion geometry approach [6], we show that by using a
product of diffusion operators in DM, we are able to approximate a modified/deformed
Laplace operator on the underlying common manifold.
Our main contribution in this paper can be summarized as follows. (i) While most
current manifold learning methods only address a single data set arising from a single manifold, we extend this basic setting and present a technique addressing several
data sets arising from multiple manifolds. (ii) We introduce a new concept of “nonlinear manifold filtering”, that is, removing the influence of the nuisance variables, with
a rigorous theoretical foundations and analysis. (iii) The ability to extract the common manifold underlying several data sets enables us to propose new methods for the
multimodal sensor fusion problem. These methods consist of a scheme to incorporate
nonlinear geometric priors into a manifold filtering procedure and are demonstrated
on illustrative (real) examples. Specifically, we show application to the analysis of
sleep dynamics as well as to seasonal pattern detection (a problem which is commonly
encountered in time series analysis and will be introduced in the sequel).
The remainder of this paper is organized as follows. In Section 2, the common
manifold model for multiple data sets collected from multimodal sensors is introduced.
With the common manifold model, the AD algorithm is formulated and its analysis is
provided. To further study the behavior of AD in the common manifold model, in Sec3

tion 3, asymptotical analysis results are provided. More details about the AD algorithm
for practical purposes are discussed in Section 4. In Section 5, we exploit the common
manifold model and design an AD algorithm to detect seasonal patterns in time series.
We provide a seasonality index to quantify the seasonal effects. In Section 6, a sleep
data set is studied, demonstrating the power of the common manifold model and the
AD algorithm in the medical field. We will show that among different pairs of sensors,
different physiological information is obtained. In Section 7, the prototypical extension problem of manifold representations is considered and discussed for the common
manifold model and the AD algorithm. Discussion and future directions are outlined
in Section 8. The proofs of the theoretical results from Section 3 are presented in the
Appendix.
2. Common Manifold Model and Alternating Diffusion
The success of the AD algorithm has been shown in different problems, for example, for sleep depth analysis [31]. To analyze the algorithm, in [30], the common
variable as well as the nuisance variables specific for each sensor are modeled by metric measure spaces; specifically, it is shown that the diffusion distance calculated from
the data via AD is equivalent to an “effective alternating-diffusion distance”, which is
defined only on one common variable shared by the two sensors [30, Theorem 5 and
Equation 79]. However, in some problems, the data may exhibit additional structures,
which could be extracted and exploited. In this paper, we consider such a case where
the common variable has a distinct geometric structure, which is modelled by a manifold. Yet, even if we assume that the system we observe remains fixed during the
observation, the data collected from different sensors might depend on the observation
procedure or be contaminated by different irrelevant information. In particular, the
common manifold might be deformed differently by different sensors. As a result, special focus is given to possible dependencies between the “common manifold” and the
sensors. These considerations are manifested in the setting presented in this section,
and our focus in the remainder of the paper is on studying the geometric information
that can be obtained from AD.
2.1. Common manifold model
Consider a setup consisting of multiple sensors observing a system or a phenomenon
of interest simultaneously. To simplify the exposition, we focus on a setup with merely
two sensors, noting that our analysis can be generalized to multiple sensors with slight
modifications.
2.1.1. Geometric model
Suppose there exists a common structure underlying the two sensor observations,
which is modeled by a low dimensional common manifold, and suppose that each sensor introduces various deformations and interferences, which are modeled by other
irrelevant nuisance structures. Mathematically, denote the common structure of interest, M, by a d-dim Riemannian manifold with the metric tensor g (i) , which depends on
the sensor, where i = 1, 2 is the sensor index, and let dg(i) denote the distance function
4

Sensor

S2

Sensor

s(2)

N1

S1

s(1)

M

The Underlying
Space Structure

N2

Figure 1: A diagram illustrating our setting with a common manifold underlying two sensor observations.

on M associated with g (i) . To accommodate possible discrepancies and deformations
between the common manifold observed through the two sensors, we allow the metrics
to be different, i.e., g (1) 6= g (2) . For i = 1, 2, let Ni be a compact metric space with the
distance function dNi , which models the irrelevant nuisance structures. Let Si denote
the observable metric space with the distance function dSi , which models the space of
the collected/accessible data by the i-th sensor. Assume that
s(i) : M × N1 × N2 → Si

(1)

is a smooth isometric embedding of M × Ni into Si modeling how the i-th sensor
collects data, where i = 1, 2. Importantly, note that s(1) ignores N2 (as it represents
the interferences specific to the the second sensor), and that s(2) ignores N1 (as it
represents the interferences specific to the first sensor). In other words, each sensor
typically acquires two structures: a deformation of the common manifold of interest
and an additional nuisance structure. The model described here implies that we do not
have access to the structures underlying each sensor (the deformed common manifold
Mi and the nuisance structure Ni ), nor to the mapping to the observable sensor space
(1). In addition, (1) means that for (x, y, z), (x0 , y 0 , z 0 ) ∈ M × N1 × N2 , where
x, x0 ∈ M, y, y 0 ∈ N1 and z, z 0 ∈ N2 , we have
dS1 (s(1) (x, y, z), s(1) (x0 , y 0 , z 0 ))2 = dg(1) (x, x0 )2 + dN1 (y, y 0 )2

(2)

dS2 (s(2) (x, y, z), s(2) (x0 , y 0 , z 0 ))2 = dg(2) (x, x0 )2 + dN2 (z, z 0 )2 .
In practice, we do not have the full access to the metric structure (2), but only to the
observable distance functions dSi .
In summary, different sensors have access to the information of interest from the

5

common geometric object M. However, the acquired information is deformed by two
different sources – the nuisance variables modeled by Ni , and the deformation induced
by each sensor modeled by the different metric on the manifold M.
A diagram of the geometric model is depicted in Figure 1.
2.1.2. Statistical model
To model the discrete data sets collected by the sensors, let (Ω, F, P ) be a probability space, where Ω is the event space, F is the sigma algebra on Ω, and P is a probability measure defined on F. Consider a random vector S : (Ω, F, P ) → M × N1 × N2 .
The dataset sampled from the i-th sensor, i = 1, 2, is modeled by a random vector
Si = s(i) ◦ S. Note that via the first sensor, only M and N1 are observed in S1 , while
via the second sensor, only M and N2 are observed in S2 .
We further assume that conditional on M, the nuisance variables modeled by N1
and N2 (introduced by the different sensors) are sampled independently. Let ν :=
νM×N1 ×N2 = S∗ P denote the induced probability measure on M × N1 × N2 . By the
assumption that conditional on M the nuisance variables are independent, we have
ν(x, y, z) = νM (x)νN1 |M (y|x)νN2 |M (z|x),

(3)

where νM (x) is the marginal distribution on M, and νNi |M (·|x) is the conditional
distribution on Ni .
To enhance the readability of the paper, we summarize the notations in Table 1.
2.2. Alternating diffusion under the common manifold model
Based on the common manifold model, we apply AD to analyze data collected simultaneously from two sensors, possibly of different modalities. The goal is to extract
the common manifold from the observed data via AD. Throughout the paper, we use
three sets of notations. The first set will be defined on the accessible data, typically
using an observable kernel and an observable diffusion. The second set is intermediate
and designed to describe the relationship between the observation and the hidden common manifold. The third set will be defined on the hidden common manifold using an
(inaccessible) effective kernel and effective diffusion.
Definition 2.1 (Observable Diffusion Kernels). Let P̃ (i) ∈ C([0, ∞)), i = 1, 2, be two
kernels that decay sufficiently fast and are associated with the two sensors. Define
P̃(i) ((x, y, z), (z 0 , y 0 , z 0 )) := P̃ (i)

(i)
(i) 0 0 0 
S (i) (s (x, y, z), s (x , y , z ))

d

√



(4)

(i)

P̃ ((x, y, z), (x0 , y 0 , z 0 ))

P(i) ((x, y, z), (x0 , y 0 , z 0 )) := R

(i)

M×N1 ×N2

P̃ ((x, y, z), (x00 , y 00 , z 00 ))dν(x00 , y 00 , z 00 )

(i)

.

where P is referred to as a diffusion kernel with bandwidth  associated with the
observable metric measure space Si = s(i) (M × N1 × N2 ).

6

Symbol (i = 1, 2)
(M, g (i) )
(Ni , dNi )
(Si , dSi )
s(i)

D(i)
P̃ (i)
D
E
D(ei )
P̃ (e)
D(e)
P̃ (Oi ,y)
ι(i)
K̃ (i)
T
K̃ (e)
(i)
expx
(i)
∇
dV (i)
Ric(i)
s(i)
Π(i)
∆(i)

Meaning
d-dim Riemannian manifold with the Riemannian metric g (i)
Compact metric space with the metric dNi
Metric space with the metric dSi
Sensor collecting data
Bandwidth parameter of the kernel function
Observable diffusion operator on Si
Kernel associated with D(i)
Observable AD operator starting from the first sensor
Marginalization operator
Effective diffusion operator on (M, g (i) )
Kernel associated with D(ei )
Effective AD operator starting from the first sensor
Observable diffusion kernel with the fixed y
Embed M into Rp so that g (i) is the induced metric via ι(i)
Diffusion kernel on (M, g (i) )
Reduced AD operator starting from (M, g (1) )
Kernel associated with T
Exponential map at x
Levi-Civita connection associated with the metric g (i)
Volume form associated with the metric g (i)
Ricci curvature of (M, g (i) )
Scalar curvature of (M, g (i) )
second fundamental form of the embedding ι(i)
Laplace-Beltrami operator of (M, g (i) )
Table 1: Summary of symbols used throughout the paper.

(i)

Note that the diffusion kernel P is in a normalized form, i.e.,
Z
P(i) ((x, y, z), (x0 , y 0 , z 0 ))dν(x0 , y 0 , z 0 ) = 1.
M×N1 ×N2

The main benefit of a diffusion kernel in such a normalized form is that the normalization helps to eliminate unwanted non-intrinsic quantities, such as the terms depending
on the specific kernel in the asymptotical analysis. See, for example, the difference
between Lemma 8 and Proposition 10 in [6]1 . For each kernel,  is referred to as the
bandwidth of the kernel. Further note that while the bandwidth may depend on the
sensor, for simplicity, we assume that the kernels share the same  value.
In this section, we show that when studied under a suitable assumption, the influ1 Note that in [30] the analysis is carried out using forward diffusion, while in this paper we use backward
diffusion for consistency with the standard diffusion maps framework presented in [6].

7

ence of the nuisance variables, that is, each metric space Ni representing the sensor
effects and observation specific influences, is erased by the AD procedure, and thus,
can be ignored. Note that in the special case when g (1) = g (2) , i.e., the common manifold is viewed similarly by different sensors, this statement was shown previously in
[30], where the focus is on the diffusion distance.
We begin with the following definition.
Definition 2.2 (Observable Diffusion Operators). For a function f ∈ C(M×N1 ×N2 ),
let D(i) : C(M × N1 × N2 ) → C(M × N1 × N2 ), for i = 1, 2, denote the diffusion
operator on the i-th sensor, defined by
Z
(i)
D f (x, y, z) :=
P(i) ((x, y, z), (x0 , y 0 , z 0 ))f (x0 , y 0 , z 0 )dν(x0 , y 0 , z 0 ).
M×N1 ×N2

(5)
Definition 2.3 (Observable Alternating Diffusion Operator and Kernel). The AD operator starting from the first sensor is defined by
D := D(2) D(1) ,

(6)

while the AD operator starting from the second sensor is defined analogously. Let
P ∈ C((M × N1 × N2 ) × (M × N1 × N2 )) be the AD kernel, defined by
P ((x, y, z), (x00 , y 00 , z 00 ))
(7)
Z
:=
P(2) ((x, y, z), (x0 , y 0 , z 0 ))P(1) ((x0 , y 0 , z 0 ), (x00 , y 00 , z 00 ))dν(x0 , y 0 , z 0 ).
M×N1 ×N2

Note that by definition, we can associate the observable AD kernel and operator:
Z
(Df )(x, y, z) =
P ((x, y, z), (x00 , y 00 , z 00 ))f (x00 , y 00 , z 00 )dν(x00 , y 00 , z 00 ).
M×N1 ×N2

(8)
The second set of notations concerns with the connection between the observations
and the hidden common manifold. Consider the observable AD starting from the first
sensor, and recall that the first sensor only sees M×N1 . Thus, by the definition of P̃ (1)
(1)
in (2) and (4), P ((x, y, z), (x0 , y 0 , z 0 )) only takes (x, y) and (x0 , y 0 ) into account. We
thus introduce the corresponding reduced intermediate diffusion kernels:
P̃(I1 ) ((x0 , y 0 ), (x00 , y 00 )) := P̃(1) ((x0 , y 0 , z 0 ), (x00 , y 00 , z 00 ))
P(I1 ) ((x0 , y 0 ), (x00 , y 00 )) := P(1) ((x0 , y 0 , z 0 ), (x00 , y 00 , z 00 )) ;

(9)

that is, the observable diffusion starting from the first sensor can be simplified by ignoring the contribution of N2 . Similarly, we define
P̃(I2 ) ((x0 , z 0 ), (x00 , z 00 )) := P̃(2) ((x0 , y 0 , z 0 ), (x00 , y 00 , z 00 ))
P(I2 ) ((x0 , z 0 ), (x00 , z 00 )) := P(2) ((x0 , y 0 , z 0 ), (x00 , y 00 , z 00 )) .
8

(10)

Clearly, by defintion, for f ∈ C(M × N1 × N2 ), D(1) f (x, y, ·) is a constant function
for a fixed (x, y). Accordingly, define the reduced intermediate function:
f (I2 ) (x, y) := f (x, y, ·) ,

(11)

which is used to describe the result of the application of D(1) . A similar argument implies that D(2) f (x, ·, z) is a constant function for fixed (x, z), and hence the definition:
f (I1 ) (x, z) := f (x, ·, z)

(12)

Lastly, we define the third set of notations, which consists of the effective counterparts of the observable diffusion operators defined on the hidden common manifold.
These notations are needed to show that in effect the AD erases the nuisance variables.
We need the following auxiliary operator that integrates out the nuisance variables.
Definition 2.4 (Marginalization Operator). Let E : C(M × N1 × N2 ) → C(M) be
the marginalization operator, which is defined by
Z
Ef (x) :=
f (x, y, z)dνN1 |M (y|x)dνN2 |M (z|x).
(13)
N1 ×N2

Clearly, the marginalization operator E is a bounded linear operator, which evaluates the marginal distribution of the collected data. Although a-priori the marginalization operator is unknown and cannot be computed from data (since we do not have
access to the hidden structure of the data, an in particular, to the nuisance variables),
we will show that an equivalent operation is attainable by the observable AD.
With the above definition, when f (x, ·, z) is a constant function for fixed (x, z),
(13) is reduced to
Z
Ef (x) =
f (I1 ) (x, z)dνN2 |M (z|x).
(14)
N2

Similarly, when f (x, y, ·) is a constant function for fixed (x, y), (13) is reduced to
Z
Ef (x) =
f (I2 ) (x, y)dνN1 |M (y|x).
(15)
N1

(e )

Definition 2.5 (Effective Diffusion Kernel and Operator). Let P i (x, x0 ) be the effective diffusion kernel associated with the i-th sensor and defined on M by
Z Z
P(ei ) (x, x0 ) :=
P(Ii ) ((x, y), (x0 , y 0 ))dνNi |M (y|x)dνNi |M (y 0 |x0 ). (16)
Ni

Ni

and let D(ei ) be the corresponding effective diffusion operator, defined by
Z
(ei )
D Ef (x) :=
P(ei ) (x, x0 )Ef (x0 )dνM (x0 ) .
M

9

(17)

(e)

Definition 2.6 (Effective Alternating Diffusion Kernel and Operator). Let P be the
effective AD kernel associated with AD starting from the first sensor, defined by
Z
P(e) (x, x00 ) :=
P(e2 ) (x, x0 )P(e1 ) (x0 , x00 )dνM (x0 ).
(18)
M

and let D(e) be the corresponding effective AD operator, defined by
D(e) := D(e2 ) D(e1 ) .

(19)

By definition we have
(D

(e)

Ef )(x) =

Z

P(e) (x, x00 )Ef (x00 )dνM (x00 ) .

(20)

M

The expansion of the effective diffusion kernel (16) according to the definition of the
diffusion kernel in (4) deserves an additional discussion. Note that


 qd (1) (x,x0 )2 +dN (y,y0 )2 
R
1
g
(1)
0 0
√
dνN1 |M (y |x )
P̃


N1 

.
P(e1 ) (x, x0 ) = E 
 qd (1) (x,x00 )2 +dN (y,y00 )2 

R
1
g
(1)
00
00
00
√
P̃
dν
(x
)dν
(y
|x
)

M
N
|M
1
M×N1

(21)
(1)

(e )

Thus, P 1 (x, x0 ) is a diffusion kernel on M; yet, for a general kernel function P̃ ,
(e )
the kernel P 1 (x, x0 ) cannot be further simplified.
(1)

Remark 2.1. In the special case when P̃ is Gaussian and the nuisance variables are
independent, it is possible to separate it into two terms consisting of the metric on the
common manifold dg(i) and the metric on the nuisance variable dNi . In other words,
in this special case, we have further access to the hidden structure of the data, and in
particular, to the metric defined on each component. Accordingly, the effective kernel
can be simplified into the following normalized form


d (x,x0 )
P̃ (ei ) gi √


P(ei ) (x, x0 ) = R
.
(22)
d (x,x00 )
P̃ (ei ) gi √
dνM (x00 )
M
With the above preparation, we are ready to state the main result of this section,
showing that after marginalization, AD constructed from the two-sensor obserations is
intimately related to a diffusion process defined on the hidden common manifold.
Theorem 2.1. For a fixed i = 1, 2, when f ∈ C(M × N1 × N2 ) is constant on Ni ,
we have
E[D(i) f ] = D(ei ) Ef .

10

(23)

Furthermore, we have
E[Df ] = D(e) Ef.

(24)

Before the proof, note that Theorem 2.1 implies that the marginalization operator
and the AD operator commute. While the effective AD operator D(e) cannot be directly
computed from data, the observable AD can be. However, in practice E is unknown,
and we need to link Df back to the effective AD operator on the common manifold.
We will further address this issue in the sequel.
Proof. Fix i = 1. The proof for i = 2 is analogous. By assumption, f (x, ·, z) is a
constant function for fixed (x, z)2 , so that we have
D(1) f (x, y, z)
Z
=
P(1) ((x, y, z), (x00 , y 00 , z 00 ))f (x00 , y 00 , z 00 )dν(x00 , y 00 , z 00 )
M×N1 ×N2
Z
=
P(I1 ) ((x, y), (x00 , y 00 ))f (I1 ) (x00 , z 00 )dν(x00 , y 00 , z 00 )
M×N1 ×N2
Z
=
P(I1 ) ((x, y), (x00 , y 00 ))Ef (x00 )dνM (x00 )dνN1 |M (y 00 |x00 ).

(25)

M×N1

Note that D(1) f (x, y, ·) is a constant function for a fixed pair (x, y). Thus, we have
E[D(1) f ](x)
Z
=
(D(1) f )(x, y, z)dνN1 |M (y|x)dνN2 |M (z|x)
N1 ×N2
Z
Z
P(I1 ) ((x, y), (x00 , y 00 ))Ef (x00 )
=
N1 ×N2

Z

M×N1

Z

=
N1

Z
=

(26)

M×N1

dνM (x00 )dνN1 |M (y 00 |x00 )dνN1 |M (y|x)dνN2 |M (z|x)
P(I1 ) ((x, y), (x00 , y 00 ))Ef (x00 )dνM (x00 )dνN1 |M (y 00 |x00 )dνN1 |M (y|x)

P(e1 ) (x, x00 )Ef (x00 )dνM (x00 ).

M

Define D(e1 ) Ef (x) := E[D(1) f ](x). Similarly, if f (x, y, ·) is a constant function for
fixed (x, y), the effective diffusion associated with the second sensor is given by
Z
(2)
E[D f ](x) =
P(e2 ) (x, x0 )Ef (x0 )dνM (x0 )
(27)
M

2 Such

a function could be obtained as the result of applying D(2) .

11

and define D(e2 ) Ef (x) := E[D(2) f ](x). To finish the proof, note that
E[Df ] = ED(2) D(1) f = D(e2 ) ED(1) f = D(e2 ) D(e1 ) Ef = D(e) Ef.

(28)

In light of Theorem 2.1, we can remark on the spectral behavior of the observable
and effective AD operators. Note that when  > 0 is finite, although D and D(e)
are compact operators, they are not self-adjoint and hence limited spectral information
can be exploited. In the sequel, we will show that asymptotically, when  approaches
0, these operators approximate a deformed Laplace-Beltrami operator. In particular,
such a limiting operator is self-adjoint, and therefore, at least asymptotically, it gives
the theoretical foundation to consider the spectral analysis of the operators at hand.
Suppose we have Dφ = λφ; that is, φ is the eigenfunction of the observable AD
operator D associated with the eigenvalue λ. Then, by the commutativity shown in
Theorem 2.1 and the linearity of E, we have
D(e) Eφ = EDφ = Eλφ = λEφ ;

(29)

that is, if λ is an eigenvalue of D with the eigenspace Eλ (D), then λ is an eigenvalue
of D(e) with the eigenspace Eλ (D(e) ) containing E(Eλ (D)). While in practice we
can only obtain an approximation of the eigenfunction φ ∈ Eλ (D) rather than Eφ ∈
Eλ (D(e) ), the knowledge of φ carries information about Eφ.
2.3. Approximating the effective AD by the observable AD
Since E is typically unknown, Theorem 2.1 cannot be directly applied to extract
the common manifold information from the given observations. To understand how
the common manifold information can be obtained from the observable AD, in this
subsection, we further explore the relation between the effective AD and the observable
AD, particularly, without taking the marginalization operator E into account.
Before describing the theoretical result, observe that for a function f defined on
M × N1 × N2 , the diffusion on the second sensor, D(2) f , is constant on N1 ; that is,
D(2) integrates information on M and N2 and no information about N1 is embodied in
D(2) f . As a result, intuitively, the additional application of D(1) , namely, D(1) D(2) f ,
not only ignores the information from N2 , but also should not bear information about
N1 as well, since D(2) f is constant on N1 . Yet, rigorous inspection shows that it is not
true in general. Observe for example (25). We have
Z
D(1) f (x, y, z) =
P(I1 ) ((x, y), (x00 , y 00 ))Ef (x00 )dνM (x00 )dνN1 |M (y 00 |x00 )
M×N1
Z hZ
i
=
P(I1 ) ((x, y), (x00 , y 00 ))dνN1 |M (y 00 |x00 ) Ef (x00 )dνM (x00 ) ,
M

N1

which is constant on z, but may still depend on y. To take a closer look, we assume
that the kernel P̃ (1) is a Gaussian kernel so that we could decouple the M and N1 and

12

have
Z
N1

P(I1 ) ((x, y), (x00 , y 00 ))dνN1 |M (y 00 |x00 )
R

P̃
N1

=
R

R

M N1

=R
M

e
e

(1)

P̃ (1)

 qd

g (1)

 qd

g (1)

−dg(1) (x,x00 )2 /

−dg(1)

(x,x0 )2 /

(x,x00 )2 +dN1 (y,y 00 )2 
√
dνN1 |M (y 00 |x00 )


(x,x0 )2 +dN1 (y,y 0 )2 
√
dνN1 |M (y 0 |x0 )dνM (x0 )


R
N1

e−dN1 (y,y

00 2

) /

0 2 /

R
N1

e−dN1 (y,y )

dνN1 |M (y 00 |x00 )

dνN1 |M (y 0 |x0 )dνM (x0 )

,

R
00 2
where N1 e−dN1 (y,y ) / dνN1 |M (y 00 |x00 ) depends on y. Therefore, while Df (x, y, z) =
D(1) D(2) f (x, y, z) is constant in z, it is not necessarily constant in y and it depends on
the geometry of N1 via the kernel integration. To better quantify this dependence, we
introduce the following definition.
Definition 2.7 (Nuisance-dependent diffusion kernel on the common manifold). For a
(N ,y)
fixed y ∈ N1 , let P 1 be the nuisance-dependent diffusion kernel on M, which is
defined in the normalized form by
P(N1 ,y) (x, x0 ) :=

dg(1) (x,x0 ) 
√

,
R
dg(1) (x,x0 ) 
(N
,y)
0)
1
√
P̃
dν
(x
M
M


P̃ (N1 ,y)

(30)

where the superscript N1 stands for the dependence on the nuisance variable corresponding to the first sensor and
Z
0 

(N1 ,y) dg (1) (x, x )
√
:=
P̃(I1 ) ((x, y), (x0 , y 0 ))dνN1 |M (y 0 |x0 )
(31)
P̃

N1
q
Z
 dg(1) (x, x0 )2 + dN1 (y, y 0 )2 
√
P̃ (1)
=
dνN1 |M (y 0 |x0 ).

N1
(N2 ,z)

Similarly, for a fixed z ∈ N2 , the nuisance-dependent diffusion kernel on M, P
is defined by
P(N2 ,z) (x, x0 )

:=

dg(2) (x,x0 ) 
√

,
R
dg(2) (x,x0 ) 
√
dνM (x0 )
P̃ (N2 ,z)
M


,

P̃ (N2 ,z)

(32)

where
P̃

(N2 ,z)

0 
g (2) (x, x )

d

√



Z
:=
N2

P̃(I2 ) ((x, z), (x0 , z 0 ))dνN2 |M (z 0 |x0 ).

13

(33)

Note that P̃ (N1 ,y) is a family of kernels defined on M that depend on the nuisance
variable y, and by definition, it is related to the effective diffusion kernel by
EP(N1 ,y) = P(e1 ) .
(N2 ,z)

A similar argument holds for P

(34)

, where we have

EP(N2 ,z) = P(e2 ) .

(35)

With this definition, we summarize the relation between the observable AD and the
effective AD in the following theorem.
Theorem 2.2. Fix (x, y, z) ∈ M × N1 × N2 and a continuous function f that is
constant in N1 . The observable AD on M × N1 × N2 starting from the first sensor
satisfies
Z hZ
i
P(N2 ,z) (x, x0 )P (e1 ) (x0 , x00 )dνM (x0 ) Ef (x00 )dνM (x00 ) .
Df (x, y, z) =
M

M

Proof. By a direct expansion based on the definition of AD, we have
Df (x, y, z)
Z
=

(36)
P(2) ((x, y, z), (x0 , y 0 , z 0 ))

M×N1 ×N2
00

Z

P(1) ((x0 , y 0 , z 0 ), (x00 , y 00 , z 00 ))

M×N1 ×N2
0 0 0

× f (x , y 00 , z 00 )dν(x00 , y 00 , z 00 )dν(x , y , z )
nZ
=
P(I2 ) ((x, z), (x0 , z 0 ))P(I1 ) ((x0 , y 0 ), (x00 , y 00 ))
M×N1 ×N2
M×N1 ×N2
o
× dνN1 |M (y 0 |x0 )dνN2 |M (z 0 |x0 )dνM (x0 ) f (x00 , y 00 , z 00 )dνN1 |M (y 00 |x00 )dνN2 |M (z 00 |x00 )dνM (x00 )
Z
Z
nZ hZ
i
=
P(I2 ) ((x, z), (x0 , z 0 ))dνN2 |M (z 0 |x0 )
P(I1 ) ((x0 , y 0 ), (x00 , y 00 ))dνN1 |M (y 0 |x0 )
M×N1 ×N2
M
N2
N1
o
0
00 00 00
00 00
00 00
× dνM (x ) f (x , y , z )dνN1 |M (y |x )dνN2 |M (z |x )dνM (x00 ) ,
Z

which, by the definition of the nuisance-dependent diffusion kernel on the common
manifold, could be reduced to
Z
nZ
o
00
Df (x, y, z) =
P(N2 ,z) (x, x0 )P(N1 ,y ) (x0 , x00 )dνM (x0 )
M×N1 ×N2

M

× f (x00 , y 00 , z 00 )dνN1 |M (y 00 |x00 )dνN2 |M (z 00 |x00 )dνM (x00 )
Z Z
hZ Z
00
(N2 ,z)
0
=
P
(x, x )
P(N1 ,y ) (x0 , x00 )f (x00 , y 00 , z 00 )
M M
N2 N1
i
00 00
× dνN1 |M (y |x )dνN2 |M (z 00 |x00 ) dνM (x0 )dνM (x00 ) .

14

(37)

To further reduce (37), note that by the assumption that f is constant in y 00 , we have
Z
00 00 00
f (x , y , z ) =
f (x00 , y 000 , z 00 )dνN1 |M (y 000 |x00 ) .
N1

Hence, the integrant quantity inside the bracket in (37) could be simplified by
Z Z
00
P(N1 ,y ) (x0 , x00 )f (x00 , y 00 , z 00 )dνN1 |M (y 00 |x00 )dνN2 |M (z 00 |x00 )
ZN2 ZN1
Z
00
=
P(N1 ,y ) (x0 , x00 )
f (x00 , y 000 , z 00 )dνN1 |M (y 000 |x00 )dνN1 |M (y 00 |x00 )dνN2 |M (z 00 |x00 )
N1
ZN2 N1
Z Z
(N1 ,y 00 ) 0
00
=
P
(x , x )dνN1 |M (y 00 |x00 )
f (x00 , y 000 , z 00 )dνN1 |M (y 000 |x00 )dνN2 |M (z 00 |x00 )
N1
N2 N1
hZ
i
(N1 ,y 00 ) 0
00
00 00
=
P
(x , x )dνN1 |M (y |x ) Ef (x00 )
=P

N1
(e1 )

(x0 , x00 )Ef (x00 ) .

(38)

By plugging (38) into (37), Df (x, y, z) is reduced to
Z hZ
i
P(N2 ,z) (x, x0 )P (e1 ) (x0 , x00 )dνM (x0 ) Ef (x00 )dνM (x00 ).
Df (x, y, z) =
M

M

Theorem 2.2 implies that when the nuisance variables y and z are fixed, the observable AD can be viewed as “an ordinary” diffusion process on the common manifold
R
(N ,z)
with the kernel M P 2 (x, x0 )P (e1 ) (x0 , x00 )dνM (x0 ), which depends on z. Hence,
by Theorem 2.1, the effective AD is related to the observable AD by taking expectation
with respect to the nuisance variable of the second sensor z, when using the observable
AD starting from the first sensor:
Z hZ
i
[D(e) Ef ](x) = E
P(N2 ,z) (x, x0 )P (e1 ) (x0 , x00 )dνM (x0 ) Ef (x00 )dνM (x00 ) .
M

M

(39)
This result indicates that we can gain access to the common manifold information via
the effective AD by viewing the observable AD as a proxy. However, the information
provided by the observable AD depends on the data, and we would only achieve an
extraction of undistorted information on the common manifold through the effective
AD if sufficient amount of data is available to “average out” the nuisance variable.
Based on this result, in the next section, we focus on studying the net result of AD
solely on the common manifold while ignoring the nuisance variables.
Remark 2.2. By the same argument we could obtain a parallel result as (39) for the
forward diffusion operator, which, if applied to the delta measure supported at two
points, recovers the effective alternating-diffusion distance considered in [30, Equation 79]. Notice that in [30], it is the effective alternating-diffusion distance, which
15

embodies the “averaged behavior” of the observable AD kernel, that is discussed,
while in this paper, we focus on the diffusion behavior at each point.
2.4. Alternating diffusion with nuisance variables modeled with manifold structures
Before closing this section, we further study how the observable AD depends on the
data. Under a stronger condition that N2 is a q-dim closed Riemannian manifold (compact and without boundaries) with the metric gN2 and that for each x00 , dνN2 |M (z 00 |x00 )
is absolutely continuous with related to the volume form dVN2 associated with gN2 , we
claim that asymptotically when  is sufficiently small and the Radon-Nikodym derivadν 2 |M (z 00 |x00 )
tive NdV
is constant (uniform sampling conditional on x00 ), Df (x, y, z) is
00
N2 (z )
almost constant in z if we start AD from the first sensor. We mention that while a more
complicated condition could be considered, to simplify the discussion we focus on this
assumption. We have
q
Z
 dg(2) (x, x00 )2 + dN2 (z, z 00 )2 
√
P̃ (2)
dVN2 (z 00 )
(40)

N2
q
Z
 dg(2) (x, x00 )2 + kuk2 
√
P̃ (2)
(1 + Ricz (i, j)ui uj + O(kuk3 ))du
=

Bz
which stems from the change of variables with the normal coordinate at z, Bz :=
z, and 
Ricz is the Ricci curvature
exp−1
z (N2 \Cz ) ⊂ Tz N2 , Cz is the cut locus of
R ∞ (2) √
2
of (N2 , gN2 ) at z. Define ϕ` (a) := 0 P̃
a + s sq−1+` ds for a ≥ 0 and
` = 0, 1, . . .. By changing the Cartesian coordinates to polar coordinates on Tz N2 ,
(40) can be recast as
r
Z
 d (2) (x, x00 )2

g
(q−1)/2
(2)

P̃
+ s2 (sq−1 + Ricz (θ, θ)sq+1 )dsdθ + O((q+2)/2 )

Rq
(41)
= (q−1)/2 |S q−1 |ϕ0 (dg(2) (x, x00 )2 /) + (q+1)/2

|S q−1 |
sy ϕ2 (dg(2) (x, x00 )2 /) + O((q+2)/2 ) ,
q

where we approximate the integration over By by the integration over Rq exploiting the
fast decay assumption of the kernel function, and sz is the scalar curvature of (N2 , gN2 )
at z.
Now, by the uniform sampling assumption, we have
Z
P(I2 ) ((x, z), (x00 , z 00 ))dνN2 |M (z 00 |x00 )
N2

R
=
R

R

M N2

P̃
N2
P̃ (2)

(2)

 qd

 qd

g (2)

g (2)

(x,x00 )2 +dN2 (z,z 00 )2 
√
dVN2 (z 00 )


(x,x00 )2 +dN2 (z,z 00 )2 
√
dVN2 (z 00 )dνM (x00 )


16

,

which by (41) is reduced to
Z
P(I2 ) ((x, z), (x00 , z 00 ))dνN2 |M (z 00 |x00 )
N2

=R

ϕ0 (dg(2) (x, x00 )2 /) +  sqz ϕ2 (dg(2) (x, x00 )2 /) + O(3/2 )
R
ϕ (d (x, x00 )2 /)dνM (x00 ) +  sqz M ϕ2 (dg(2) (x, x00 )2 /)dνM (x00 ) + O(3/2 )
M 0 g (2)

= R

ϕ0 (dg(2) (x, x00 )2 /)
sz
+  G(x, x00 ) + O(3/2 ).
00 )2 /)dν (x00 )
q
ϕ
(d
(x,
x
(2)
M
M 0 g

where
ϕ2 (dg(2) (x, x00 )2 /)
ϕ (d (x, x00 )2 /)dνM (x00 )
M 0 g (2)
R
ϕ0 (dg(2) (x, x00 )2 /) M ϕ2 (dg(2) (x, x00 )2 /)dνM (x00 )
R
−
( M ϕ0 (dg(2) (x, x00 )2 /)dνM (x00 ))2

G(x, x00 ) = R

The above derivation implies that under the assumption that the nuisance variable lies
R
(I )
on a compact manifold and is sampled uniformly, the quantity N2 P 2 ((x, z), (x00 , z 00 ))dνN2 |M (z 00 |x00 )
is almost constant in z, and the dependence on z is of order . Hence, under this strong
condition, we conclude from Theorem 2.2 and from Definition 2.7 that Df (x, y, z) is
almost constant in z if we start AD from the first sensor, and hence by (39), we could
recover the effective AD via the observable AD with a higher order error. It is worthwhile noting that an analogous result can be derived for AD starting from the second
sensor, if N1 is a manifold with the same assumptions.
3. Analysis of the alternating diffusion under the common manifold model
In this section, we study the effective AD under the manifold setup. Based on
the discussion in Section 2, we present the analysis and proof under a reduced setting
which does not contain the nuisance variables. Yet, the results of the analysis under
this setting are transferable to the general setting with the nuisance variables.
Before presenting the analytical results, we take a closer look at the available metrics when we analyze data. First, theoretically, if the metric dSi for the dataset faithfully
reflects the geodesic distance on M, then the analysis becomes simple. However, it is
usually not the case in practice and the best we could expect is that the metric dSi
for the dataset provides a good approximation of the geodesic distance on M; that
is, dg(i) (x, x0 ), for i = 1, 2, in (2) can be approximated from the observations and
their ambient metrics dSi . Denote the approximation of dg(i) (x, x0 ) by d¯g(i) (x, x0 ).
The error introduced by the difference between the metric we have for the dataset
and the geodesic distance might not be easily quantified in general. Thus, for mathematical tractability, we make further assumptions regarding the discrepancy between
dg(i) (x, x0 ) and d¯g(i) (x, x0 ). When the metric is well designed, we can assume that for
sufficiently close x and x0 , dg(i) (x, x0 ) and d¯g(i) (x, x0 ) are close up to a higher order
error. Moreover, in some situations we are able to directly quantify the error, which
helps us to further quantify the information that can be extracted from AD.
17

We now recast the formulation of the problem and the definitions of the diffusion
kernels and operators from Section 2 under the reduced setting, and introduce notations
used in this section for the analysis. Suppose the common manifold M is smoothly
embedded into Rp via ι(i) , i = 1, 2, with the induced metric g (i) and hence the induced
distance function dg(i) . We sample M via two random vectors, Xi := ι(i) ◦ X, where
X : (Ω, F, P ) → M, with the induced measure on M, denoted by µM := X∗ P . This
model entails that under the reduced setting we are not able to access M directly, but
only through the two sensors via Xi ∈ Rp . In this case,
d¯g(i) (x, x0 ) = kι(i) (x) − ι(i) (x0 )k,

(42)

where k · k means the Euclidean distance between any two samples x and x0 from M;
in other words, we only approximate the geodesic distance on M using the Euclidean
distance. In general, other approximations could be used, e.g. based on an embedding
to a manifold or a well-designed metric space, but to simplify the analysis we focus on
the Euclidean space. Clearly, ι(1) (M) and ι(2) (M) are diffeomorphic to each other
−1
via a differomorphism Φ so that Φ = ι(2) ◦ ι(1) . See Figure 2 for an illustration of
this reduced setting.
Reduced Sensor 2
Space

Reduced Sensor 1
Space

M1

M2
◆(1)

◆(2)

M

The Underlying
Common Manifold

Figure 2: A diagram illustrating the reduced setting without nuisance variables.

We study the effective AD starting from the first sensor, defined in (19), by studying
the following diffusion process on M. Note that the effective AD starting from the
second sensor can be analyzed in an analogous way.
Definition 3.1 (Reduced Alternating Diffusion Operator). Take two kernels K̃ (1) and
K̃ (2) , and define the reduced (without the nuisance variables) AD operator T : C(M) →

18

C(M) by
(2)

K̃ (x, x0 )

Z
T f (x) :=
M

R

(2)

K̃ (x, x̄)dµM (x̄)
(1)
hZ
i
K̃ (x0 , x00 )
f (x00 )dµM (x00 ) dµM (x0 )
×
R
(1) 0
M
K̃ (x , x̄)dµM (x̄)
M

(43)

M

where f ∈ C(M), and
K̃(i) (x, x0 ) := K̃ (i)



kι(i) (x) − ι(i) (x0 )k
√



(44)

for i = 1, 2.
(e)

Definition 3.2 (Reduced Alternating Diffusion Kernel). Let K̃ (x, x00 ) be the reduced (without the nuisance variables) AD kernel, defined by
K̃(e) (x, x00 )

(1)

(2)

Z
:=
M

K̃ (x, x0 )K̃ (x0 , x00 )
dµM (x0 ).
R
(1) 0
(x
,
x̄)dµ
(x̄)
K̃

M
M

By definition, we have
Z
Z
K̃(e) (x, x̄)dµM (x̄) =

K̃(2) (x, x̄)dµM (x̄).

(45)

(46)

M

M

and
(e)

R
T f (x) =

M

K̃ (x, x00 )f (x00 )dµM (x00 )
R
(e)
K̃ (x, x̄)dµM (x̄)
M

(47)

The notation used in the remainder of this section is as follows. Note that ι(i) (M),
i = 1, 2, is now a sub-manifold of Rp , and the distance function g (i) is induced from the
canonical metric of Rp . Let dV (i) denote the measure associated with the Riemannian
volume form induced from g (i) . We use the notation ∇(i) , exp(i) , Ric(i) , and Π(i)
to denote the covariant derivative, the exponential map, the Ricci curvature and the
second fundamental form associated with ι(i) , respectively.
We start from the following assumptions.
Assumption 3.1. (A1) The manifold M is d-dim, compact and without a boundary. It
is embedded into Rp via ι ∈ C 4 (M, Rp ) with the metric g induced from the
canonical metric of Rp .
(A2) The kernel functions satisfy the following conditions. For i = 1, 2, K̃ (i) ∈
C 2 ([0, ∞)) are positive, decay exponentially fast and K̃ (i) (0) > 0. Further,
2
2
there exists c1 , c2 > 0 so that K̃ (i) (t) < c1 e−c2 t and |[K̃ (i) ]0 (t)| ≤ c1 e−c2 t .
R
(i)
Denote µl,k := Rd kxkl ∂ k K̃ (i) (kxk)dx < ∞, where l ∈ {0} ∪ N, and ∂ k is the
(i)

k-th derivative, for k = 0, 1, 2. Assume µ0,0 = 1.
19

√
(A3) The bandwidth of the kernel, denoted by , satisfies 0 <  < min{τ, inj(M)},
where τ is the reach of the manifold [40] and inj(M) is the injectivity of M [12,
p 271].
(A4) Assume that dµM is absolutely continuous with related to dV (i) , and we denote
dµM
pi := dV
(i) as the probability density function (p.d.f.) of X on M by the RadonNikodym theorem. Furthermore, assume that pi ∈ C 4 (M) so that 0 < min pi ≤
max pi , where i = 1, 2.
We mention that Assumption (A4) is necessary for the sake of analyzing the asymptotical behavior of AD without the influence of the nuisance variables. In particular it
contributes to the symmetric argument used in the proofs of the theorems.
With the above preparation, we are ready to state our main results. The first theorem states that the effective AD kernel defined under a model ignoring the nuisance
variables behaves essentially like an ordinary diffusion kernel in the normalized form.
In particular, when both kernels are Gaussian, the effective AD kernel is Gaussian as
well.
Theorem 3.2. Suppose Assumptions (A1)-(A4) hold. Take f ∈ C 3 (M), 0 < γ <
(1)
1/2 and x, x00 ∈ M so that x00 = expx v, where v ∈ Tx M and kvkg(1) ≤ 2γ .
(1)

(2)

Denote Rx = [d expx |0 ]−1 [dι(2) ]−1 ∇Φ[dι(1) ][d expx |0 ] : Rd → Rd Then, when 
is sufficiently small, the following holds:
Z
√ 
K̃(e) (x, x00 ) =
K̃ (2) (kRx wk) K̃ (1) kw − v/ k dw + A2, (1, v) + O(3/2 )
Rd

where A2, (1, v) is defined in (A.37), which decays exponentially, A2, (1, 0) is of order
2
O(1), and A2, (1, −v) = A2, (1, v). In particular, if K (1) (t) = K (2) (t) = e−t /π d/2 ,
we have
Z
2 −1/2
Rx vk2 /
e−k(I+Rx )
p
K̃ (2) (kRx wk) K̃ (1) (kw − vk) dw =
.
(48)
det(I + Rx2 )
Rd
On the other hand, when kvkg(1) > 2γ ,
K̃(e) (x, x00 ) = O(3/2 ).

(49)

Note that Rx embodies the difference between the two metrics. AD consists of
two diffusion steps: the first, carried out by K̃ (1) , respects the metric g (1) , and the
second, carried out by K̃ (2) , respects the metric g (2) . In order to study the integrated
behavior of the two consecutive different diffusion steps, we quantify the overall effect
using g (1) via Rx . When ι(1) = ι(2) , that is, when Φ is the identity map, then Rx is
reduced to the identity as well, and the common manifold setup is reduced to the setup
considered in [30].
Based on the behavior of the reduced effective AD kernel K̃ (e) studied in Theorem
3.2, we study the asymptotic behavior of the reduced AD operator T (without the in20

fluence of the nuisance variables). The second theorem states that asymptotically T is
a deformed Laplace-Beltrami operator defined on (M, g (1) ).
Theorem 3.3. Suppose f ∈ C 3 (M). Fix normal coordinates around x associated with
g (1) and g (2) so that {Ei }di=1 ⊂ Tx M is orthonormal associated with g (1) . Set Rx =
(2)
(1)
[d expx |0 ]−1 [dι(2) ]−1 ∇Φ[dι(1) ][d expx |0 ] and by the singular value decomposition
(SVD) Rx = Ux Λx VxT , where Λx = diag[λ1 , . . . , λd ]. Then, when  is sufficiently
small, the AD without the nuisance variables starting from g (1) satisfies
T f (x) = f (x) +

(2) d
(1)
(1)
 (1) 2
µ2,0 X
2∇Ei f (x)∇Ei p1 (x) 
λ
∇
f
(x)
+
i
Ei ,Ei
2d2 i=1
p1 (x)
(1)

+

µ2,0 h (1)
2∇(1) f (x) · ∇(1) p1 (x) i
+ O(3/2 ).
f
(x)
+
∆
2d2
p1 (x)

In particular, when ι(1) = ι(2) , that is, Rx = Id , for every x ∈ M we have
(1)

T f (x) = f (x) + 

(2)

µ2,0 + µ2,0 h (1)
∇(1) p1 (x) · ∇(1) f (x) i
∆ f (x) +
+ O(3/2 )
2
2d
p1 (x)
(50)

This theorem implies that when Φ is not the identity map, then the obtained infinitesimal generator of the AD process is a deformed Laplace-Beltrami operator of M
associated with g (1) . In particular, when K (1) = K (2) , we have
T f (x) = f (x) +

(1) d
(1)
(1)
 (1) 2
µ2,0 X
2∇Ei f (x)∇Ei p1 (x) 
(1
+
λ
)
∇
f
(x)
+
+ O(3/2 ).
i
Ei ,Ei
2d2 i=1
p1 (x)

In addition, when Φ is the identity map, then the infinitesimal generator of the AD
process is precisely the Laplace-Beltrami operator of M associated with g (1) .
The proofs of Theorem 3.2 and Theorem 3.3 appear in Appendix A. The theorems
immediately lead to the following corollary, which describes the asymptotical behavior
of the effective AD operator in a model including the nuisance variables, studied in (19)
and is associated with the AD defined from data, starting from the first sensor (36). We
mention that since in general we are not able to convert (21) into a normalized kernel,
(N ,y)
we need to directly study the observable diffusion kernel P 1 .
Corollary 3.3. Suppose f ∈ C(M × N1 × N2 ) so that Ef ∈ C 3 (M) and νM is
absolutely continuous with respect to the Riemannian measure induced from g (1) so
dνM
4
that p1 = dV
(1) ∈ C (M). Then by the definitions of the AD operators in (6) and

21

(19), and by the commutativity from (24), we obtain
(1)
(1)
d
2∇Ei Ef (x)∇Ei p1 (x) 
C2 X  (1) 2
(51)
λ
∇
Ef
(x)
+
i
Ei ,Ei
2d2 i=1
p1 (x)
C1 h
2∇(1) Ef (x)∇(1) p1 (x) i
+ 2 ∆(1) Ef (x) +
+ O(3/2 ).
2d
p1 (x)

EDf (x) = Ef (x) +

where Ci , i = 1, 2, are constants depending only of the chosen kernel P̃ (i) .
The proof of this corollary appears in Appendix A as well. The corollary states that
after marginalization, the AD operator computed from data is asymptotically (when  is
small) a deformed Laplace-Beltrame operator on the common manifold M associated
with g (1) .
We end this section with two closing remarks. First, we note the importance of
the order of the kernels consisting the AD kernel, as implied by the analytic results
(Theorem 3.2 and Theorem 3.3). Second, these results further show that in order to
compare the outcome of AD starting from the first sensor to those of AD starting from
the second sensor, the deformation of the common manifold in each of the two sensors
has to be taken into account. This issue is illustrated by an example in Section 5.
4. Alternating diffusion algorithm
The AD algorithm is summarized in Algorithm 1, which is a direct discretization
of the observable AD operator.
4.1. Some facts about the spectrum
The discretization of AD is implemented by a direct multiplication of two normalized affinity matrices (kernels), i.e., P := P1 P2 , and it is not obvious that the spectral
theorem can be applied to P. While both P1 and P2 are symmetrizable, the asymmetric matrix P is not symmetrizable or normal in general, so that it does not necessarily
have a real spectrum or a complete eigen-basis. Note that P is a row stochastic matrix,
so that the operator norm of P is bounded by 1, as both P1 and P2 are bounded by 1.
Moreover, by the Perron-Frobenius theory if both P1 and P2 are primitive, P has an
eigenvalue 1, which is simple and is the only eigenvalue with radius 1, corresponding
to the eigenvector [1, . . . , 1]T /n ∈ Rn . Since we use two positive kernels P̃ (1) and
P̃ (2) and since all pairwise affinities are evaluated, these assumptions are satisfied. To
be more precise, define W := W1 D−1
2 W2 , where Di := diag(Wi 1) for i = 1, 2 and
1 = [1, . . . , 1]T ∈ Rn . In general, W is asymmetric but non-negative. We thus have
P = D−1
1 W,

(52)

and note that D1 = diag(W1). In other words, the discretized AD operator P could be
viewed as a stochastic diffusion operator on a directed graph with imbalanced weights
on edges. In general, although W1 and W2 could be made positive definite if the
chosen kernel is Gaussian by the Bochner theorem when there are finite points, the
spectrum of W might not be real, and hence the spectrum of P.

22

Algorithm 1 Alternating Diffusion algorithm
Input: Two data sets Xl := {xl,i }ni=1 ⊂ Rp , where l = 1, 2, are given. x1,i and x2,i are
sampled simultaneously from the sensors, for all i = 1, . . . , n.
Parameters: Pick two positive kernels P̃ (1) and P̃ (2) which decay fast enough. Fix a positive
integer K ≤ n and  > 0.
Output: The first K singular values and singular vectors of the AD starting with the first
or the second sensor.
n×n
1. Build the first affinity matrix W
based on the first data set X1 by
√1 ∈ R
(1)
W1 (i, j) = P̃ (kx1i − x1j k/ ).
2. Build the second affinity matrix W√2 ∈ Rn×n based on the second data set X2
by W2 (i, j) = P̃ (2) (kx2i − x2j k/ ).
3. Build the first diffusion kernel P1 ∈ Rn×n based on the first data set X1 by
(i,j)
P1 (i, j) = PnW1W
.
1 (l,j)
l=1

4. Build the second diffusion kernel P2 ∈ Rn×n based on the second data set X2
(i,j)
.
by P2 (i, j) = PnW2W
2 (l,j)
l=1
5. Run the AD starting with the first sensor with the diffusion operator P2 P1 and
(1)
obtain the first K right singular vectors ui ∈ Rn , i = 1, . . . , K.
6. Run the AD starting with the second sensor with the diffusion operator P1 P2
(2)
and obtain the first K right singular vectors ui ∈ Rn , i = 1, . . . , K.

4.2. Alternating diffusion map algorithm and alternating diffusion distance via SVD
For the asymmetric matrix P  ∈ Rn×n , we
 can always consider the SVD; that
is, P = UΛVT , where U = u1 , . . . , un ∈ O(n) and V = v1 , . . . , vn ∈
O(n) contain the left and right singular vectors, vi and ui , i = 1, . . . , n, and Λ =
diag[σ1 , . . . , σn ] is a n × n diagonal matrix with the corresponding singular values, σi ,
i = 1,P. . . , n, on the diagonal entries. Therefore, the matrix P can be decomposed as
n
P = `=1 σ` u` v`T . By observing that
2

kPei − Pej k =

n
X
`=1

(σ` v` (i) − σ` v` (j))2 ,

(53)

where ei ∈ Rn is the unit vector with the i-th entry 1, consider the following AD map:
Φ : i 7→ [σ1 v1 (i), σ2 v2 (i), . . . , σn vn (i)]T ∈ Rn ,

(54)

and the corresponding AD distance between i and j defined as kΦ(i) − Φ(j)k. Note
that we have
kΦ(i) − Φ(j)k = kPei − Pej k ,
(55)
and kPei −Pej k is analogous to the effective alternating-diffusion distance considered
in [30]. Yet, the AD map and distance here are based on the SVD of P rather than on
23

the eigenvalue decomposition (EVD) of P as presented in [6]3 .
The main benefit of considering kΦ(i) − Φ(j)k as the AD diffusion distance lies
in the capability to handle noise. For example, based on recent advances in analyzing
the SVD in random matrix setup [13], one could balance between the accuracy of the
required AD distance and the noise influence on the result via the truncation scheme.
A systematic study of this direction will be reported in future work.
In general, to realize the idea of “diffusion” [6], consider Pt , where t > 0. Broadly,
if P is viewed as a transition probability matrix of some Markov chain defined on the
samples, Pt consists of the transition probabilities in t steps. In contrast to the standard
diffusion geometry framework [6], which relies on the EVD of P as well as on the tight
connection between the EVD of P and of Pt , the singular values and singular vectors
t
T
of Pt are
 not directly related to those of P. Denote theSVD of P as Ut Λt Vt , where
Ut = ut,1 , . . . , ut,n ∈ O(n), Vt = vt,1 , . . . , vt,n ∈ O(n) contain left and right
singular vectors and Λt = diag[σt,1 , . . . , σt,n ] is a n × n diagonal matrix with the
corresponding singular values on the diagonal entries. Then, by the same argument as
the above, we could consider the following AD map with time t > 0:
Φt : i 7→ [σt,1 vt,1 (i), σt,2 vt,2 (i), . . . , σt,n vt,n (i)]T ∈ Rn ,

(56)

and the corresponding AD distance with diffusion time t > 0 as kΦt (i) − Φt (j)k. As
before, we have kPt ei − Pt ej k = kΦt (i) − Φt (j)k. We remark that a similar approach
is also considered in [37], where the authors consider the time-coupled diffusion maps.
4.3. Non-uniform sampling issue
Typically, the dataset is sampled non-uniformly from the common manifold; that is,
p1 or/and p2 might be non-constant. In the manifold learning society, it has been well
known that the non-uniform sampling effect has possible negative effect [6], whereas
in some cases it is beneficial [39]. One way to reduce the influence of the non-uniform
sampling is the α-normalization proposed in [6]. While the application of this normalization to our AD setup is straightforward, here we summarize the procedure and refer
readers with interest to [6, 44] for details. For 0 ≤ α ≤ 1,  > 0 and a probability
density function p defined on M, we define the following functions for i = 1, 2:
p(i)
 (x) :=

Z

d(i)
,α (x) :=

Z

(i)

K̃(i) (x, y)pi (y)dV (i) (y),

(i)
K̃,α
(x, y) :=

(i)
K̃,α
(x, y)pi (y)dV (i) (y),

(i)
K,α
(x, y) :=

M

K̃ (x, y)
, (57)
α
α
(i)
p  (x)p(i)  (y)
(i)

M
(i)

K̃,α (x, y)
(i)

d,α (x)

.
(i)

Here, p (x) is related to the estimation of the p.d.f. pi at x, denoted by −d/2 p (x).
(i)
The practical meaning of K̃,α (x, y) is a new kernel function at (x, y) adjusted by the
estimated p.d.f. at x and y; that is, the kernel is “normalized” to reduce the influence
3 As described above, in contrast to the standard construction of the kernel, in AD, P does not necessarily
have real eigenvectors.

24

(i)

of the non-uniform p.d.f. pi . The kernel K,α (x, y) is thus another diffusion kernel
(i)
associated with K̃,α (x, y). Since the proof of the α-normalization follows the same
lines as those in [6, 44], we do not present it here.
5. Application to seasonal pattern detection
The ability to extract the common latent manifold underlying multiple manifolds
gives rise to a new approach for detecting latent seasonal patterns in time series [9, 20,
21, 46, 19]. While most existing methods for seasonal pattern detection are based on
trigonometric or Fourier-based analysis, time frequency analysis and parametric estimation, we take a geometric data analysis standpoint. Our approach allows the detection of seasonal patterns hidden in the data and obscured by the observation modality
in addition to noise. Specifically, it is designed to accommodate nonlinearities masking
the information of interest. For example, consider a simple 1-dimensional periodic pattern, represented by the harmonic function cos(2πω0 x) with a single base frequency
ω0 > 0. Suppose this seasonal pattern is distorted by an unknown nonlinear obser√
vation function ι, which takes the form of y when y = cos(2πω0 x) > 0 and y 2
otherwise. Even in this caricature example, neither by trigonometric function matching, Fourier-based analysis, nor parametric estimation, the wrong frequency information might be recovered. In contrast, typical time frequency analysis approaches might
provide redundant information, causing ambiguity.
The primary idea relies on the observation that the geometric manifold representation of a pure 1-dimensional seasonal pattern, which is usually represented by a simple
sinusoidal process, is a 1-dimensional sphere S 1 := {(cos(θ), sin(θ))T ∈ R2 | θ ∈
[0, 2π)}, where θ denotes the intrinsic phase of the observed seasonal oscillation; that
is, the intrinsic manifold of interest, which is to be recovered in order to discover the
seasonal dynamics, is M = S 1 . The time series sampled from the underlying 1dimensional seasonal pattern, denoted by f (t), can converted to a high dimensional
time series, denoted by X , in an observable space S by using a lag map [45]. Clearly,
the geometric structure of the observed (now) high dimensional points might be different from M, as the data might be contaminated by the observation modality and the
embedding process. Using the notation from Section 2, the sampled high dimensional
data X can be modeled as a smooth embedding ι : M × N → S, where N is a metric
space describing various nuisance/interference variables as well as artifacts introduced
by the conversion. We note that while the exposition here focuses on detecting 1dimensional patterns, the extension of the formulation and the detection algorithm to
higher dimensional seasonal patterns is straightforward. For example, d independent
seasonal patterns can be represented by a d-torus: T d = S 1 × · · · × S 1 .
Our approach uses the following rationale: if the data at hand admit a pure 1dimensional seasonal pattern, then, according to Corollary 3.3, we can design a reference manifold to detect its seasonality. We consider a generated reference sample set
from the reference canonical M := S 1 associated with the frequency ω > 0, denoted
by Y, which is the lag map of the time series cos(2πωt). This reference data set from
S 1 represents a pure seasonal pattern. To be more precise, the mapping associated
with the first sensor s(1) : M → Rp is determined by the lag map, where p ∈ N
is the chosen lag step and N1 = ∅; the mapping associated with the second sensor
s(2) : M × N2 → Rp is the composition of the procedure generating the recorded
25

time series from M and the lag map with the lag step p, where N2 represents a space
of possible interferences and nuisance variables introduced during the data acquisition
procedure and the lag embedding process. Here, the given sample set is in X ⊂ Rp
and the goal is to recover M, namely, the underlying seasonal pattern, using AD from
X , the space of the data, and Y, the space of the reference samples. If the given sample
set X embodies a seasonal pattern with the base frequency ω, then it assumed to lie
in a space consisting of an image of S 1 (representing the underlying seasonal pattern
contaminated by the interference/nuisance variables). As a result, by applying AD to
obtain a parameterization of the common manifold, due to the pure seasonal pattern
in the generated reference sample set, the common manifold extracted by the application of AD to the given data set X and the reference set Y is S 1 , i.e., the desired
seasonal pattern underlying the data. Since the given data set does not necessarily exhibit a seasonal pattern with the base frequency ω used to generate the reference set, the
difference between the resulting quantities associated with the AD over the common
manifold model and the diffusion maps over the reference manifold S 1 can be used to
define a seasonality index. Denote the top nontrivial eigenvector of AD over the common manifold model as ψ, and the top nontrivial eigenvector of DM over the reference
manifold associated with the frequency ω > 0 as ψr . Then, define the seasonality index
for the inherent periodicity as
SI(ω) = kψ̂ − ψ̂r k2 ,

(58)

where ψ̂ means the Fourier transform of ψ. Note that this approach is feasible since the
dataset is well-ordered in time. This choice of difference is simple and efficient, yet it
is shown empirically to provide sufficiently accurate results in the tested applications.
Our proposed algorithm for computing the seasonality index for the detection of 1dimensional seasonal patterns based on AD is summarized in Algorithm 2.
To exemplify the algorithm, as well as the theoretical results, we first test our algorithm on a toy problem. Consider a data set of 1000 samples from a 2-dim time series
given by
xn = (cos(2π(n∆t + sin(2πn∆t))), sin(2π(n∆t + sin(2πn∆t))))
1
where ∆t = 1. Clearly, {xn }1000
n=1 can be viewed as a set of samples from S with
nonuniform sampling. Alternatively, we may view the set as a pure 1-dim seasonal
pattern with base frequency ω0 = 1/1000, distorted by a periodic function.
According to Algorithm 2, we create the following 1000 reference samples

rn = (cos(2πω0 n), sin(2πω0 n))
which can be viewed as uniform samples from S 1 .
By Theorem 3.3, where the common manifold is M = S 1 , the metrics are equal
(1)
(ι = ι(2) ), and the sampling density on the reference manifold is uniform, the effective AD operator starting from the reference set Y is given by
(1)

T f (x) = f (x) + 

µ2,0 (1)
∆ f (x) + O(3/2 )
2
26

Algorithm 2 Seasonality Index Computation Based on Alternating Diffusion
Input: a data set {xi }ni=1 in Rd and a tested seasonality frequency ω.
Output: a seasonality index SI(ω) and a baseline seasonality index SIbl (ω).
e ∈ Rn×n based on the given data {xi }n .
1. Build a Gaussian kernel P
i=1
2. Optional, for illustrative purposes:
e and obtain the leading eigenvector ψbl ∈ Rn .
Apply diffusion maps to P
3. Create a synthetic reference manifold:
• Simulate reference samples ri = cos(2πωi), i = 1, . . . , n.

• Create a lag map of samples yi ∈ Rl with arbitrary lag l from ri .
e r ∈ Rn×n from {yi }n .
• Build a Gaussian kernel P
i=1

e r and obtain the leading eigenvector ψr ∈ Rn .
• Apply diffusion maps to P

e r and P
e and obtain the leading eigenvector ψ ∈ Rn .
4. Apply AD to P
5. Compute the seasonality index:
SI(ω) = kψ̂ − ψ̂r k2

(59)

where v̂ denotes the discrete Fourier transform of v.
6. Optional, for illustrative purposes:
Compute the baseline seasonality index:
SIbl (ω) = kψ̂bl − ψ̂r k2

(60)

which is the Laplace-Beltrami operator. As a result, the eigenvectors of the AD operator approximate the eigenfunctions of the Laplace-Beltrami operator, and in particular
in 1-dim, the eigenvector associated with the largest eigenvalue of the AD operator (or
the smallest eigenvalue of the Laplace-Beltrami operator) is ψ(n) ≈ sin(ω0 n), and
hence, ψ̂(ω) ≈ δ(ω − ω0 ), where δ is the Delta function.
In contrast, the effective AD operator, starting from the data set X , where the
nonuniform density is given by p(x) = 1/(2π(1 + 2π cos(2πx)) is given according
to Theorem 3.3 by
(2)

T f (x) = f (x) + 

µ2,0 h (2)
∇(2) p(x) · ∇(2) f (x) i
∆ f (x) +
+ O(3/2 )
2
p(x)

which is a deformed Laplace-Beltrami operator. As a result, the eigenvector associated
with the largest eigenvalue of the AD operator is no longer a sinusoid with a single
base frequency ω0 .
Figure 3 presents the experimental results of such two AD applications, exemplifying our theoretical results, and in particular, the effect of the order of the singleview operators composing the AD operator. In addition, it shows the capability of
Algorithm 2 to accurately extract a seasonal pattern from a data set despite being
27

masked. Note that as discussed in Section 3, since the common manifold components
observed in the data set and in the reference set are not identical, when comparing
the eigenvectors, we need to pull back the resulting eigenvector according to the diffeomorphism. The Matlab code of the numerical implementation is available here:
https://github.com/ronenta2/alternating-diffusion.git.
Next, we test our algorithm on a time series of weekly U.S. finished motor gasoline
products supplied from February 1991 to July 2005 in units of thousands of barrels
per day. Such a time series has an annual seasonal pattern with non-integer period
ω0 = 365.25/7 ≈ 52.179 [9], and it is depicted in Figure 4 (a).
In Fig. 4 (b) we plot the seasonality index, computed according to (59), as a function of the tested frequency ω, ranging from 5 to 60. Figure 4 (c) zooms into the range
of interest around the true base frequency (ω0 = 52.179). As we observe, our index
indeed attains the minimal value at the true (non-integer) based frequency, and thus,
accurately identifies the seasonal pattern underlying this data set. To highlight the contribution of the AD in extracting the seasonality pattern hidden in the data, we compare
the proposed seasonality index with a naive baseline index, computed based on the
comparison between the leading eigenvectors of the reference manifold and the data
obtained by two separate DM applications (without applying AD). This baseline index
is designed to show that the seasonal pattern is hidden in the given data and cannot
be identified simply from the leading eigenvector obtained by, for example, a direct
DM application to the given data set. Figure 4 (d) depicts this naive index, computed
according to (60), as a function of the tested frequency. In contrast to our index, this
index does not capture the underlying annual seasonal pattern.
In Figs. 4 (b)-(c) we observe that the correct seasonality pattern is distinctly identified. In contrast, Fig. 4 (d) shows that the application of AD is critical; it implies
that the seasonality pattern is not the only information hidden in the data and cannot be
simply extracted without the “filtering” of possible nuisance variables and acquisition
procedure deformations attained by the AD procedure.
6. Application in sleep research
The extraction of the common manifold underlying two (or more) data sets from
different sensors can be viewed as, and used for, nonlinear manifold filtering. To
present the main idea and demonstrate its potential, we apply our technique to sleep
data. Sleep is global and systematic physiological activity, which manifests complicated temporal dynamics [32]. Since the sound dynamics associated with sleep is associated with versatile diseases, a full understanding of the sleep dynamics is critical. For
this purpose, different signals from different sensors are typically recorded, aiming to
measure different physiological aspects. For example, in many specialized sleep labs,
various signals such as EEG, ECG, respiration or EMG are recorded.
Recently, we studied the problem of automatic sleep stage identification from a
geometric analysis/manifold learning perspective. In [49], we considered a model in
which there exists a hidden process restricted to a low-dimensional Riemannian manifold governing these measured signals, whereas we only have access to the recorded
time series/signals, which represent physiological processes deformed by the observation procedures. This model, which is referred to as the empirical intrinsic geometry
28

(EIG), is motivated by the assumption that each measurement (e.g., an EEG channel
measuring brain activity or a chest belt measuring respiration) can be affected by numerous factors related to the equipment (e.g., the specific type of sensors and their
exact positions) or to noise, but we have interest in the true intrinsic variable associated
with the sleep cycle.
In [31], we extended this model to multiple sensors and showed that through AD,
we can better capture information on the sleep dynamics from multimodal respiratory
signals, compared with the analysis based only on a single respiratory signal. More
specifically, we applied AD to extract the common source of variability in abdominal
motions and in airflow. We showed that the common variable of these two respiratory
signals (measured by different types of instruments) recovers a true physiological hidden process that is well correlated with the sleep stage. The underlying assumption
is that physiologically there is a dominant controller of the respiratory process, which
cannot be explicitly modeled or accessed in practice. Yet, this controller can be observed by different sensors, for example, from the chest belt movement or the air flow
signal. On the one hand, different observations capture different, complementing aspects of the information about this controller. On the other hand, this information is
deformed in different manners depending on each specific sensor.
Here, we further demonstrate the usefulness of AD, now from a (nonlinear manifold) filtering perspective. Moreover, we attempt to devise a proof-of-concept example
and to associate the common manifold and the output of AD with physiology. We apply
AD (separately) twice to two pairs of channels. In the first application, we use an EEG
channel (O2A1) and a belt sensor located on the chest measuring respiration movements (THO). In the second application, we use the same EEG channel (O2A1) and a
different belt sensor located on the abdomen (ABD), designed to measure respiration
movements as well. See [49, Section III.A] for more details about the experimental
setup and measurements.
The leading eigenvectors resulting from the two AD applications are presented in
Figure 5(a). It is visually evident that the two leading eigenvectors resulting from the
two applications follow the same patterns, a fact that reflects the validity of the assumption that there is a common controller guiding the entire system. In other words, while
the common controller information is masked by each modality (measurement/sensor),
the common controller is extracted by AD. Importantly, note that in both applications,
we can view the AD procedure as a nonlinear filtering of the manifold underlying the
EEG channel from nuisance factors by using the respective manifolds underlying the
two respiratory signals.
To further examine this model, we repeat the procedure and apply AD to a third
pair of measurements. The new pair consists of the same EEG channel (O2A1) and an
EMG channel measuring muscle movements. Figure 5(b) presents the leading eigenvectors of AD applied to the first pair (O2A1 and THO) and the third pair (O2A1 and
EMG). In contrast to Figure 5(a), we now observe that the two leading eigenvectors
resulting from the two applications of AD share less similar patterns. For comparison
purposes, we report the `2 distance between the obtained eigenvectors: the distance
between the two eigenvectors in Figure 5(a) is 0.8 and the distance between the two
eigenvectors in Figure 5(b) 1.2. According to our manifold filtering interpretation, this
result indicates that the manifold filtering of the EEG channel via the respiratory signal
29

and the manifold filtering via the EMG signal are different. This finding can be physiologically explained – the controller common to the EEG signal and the EMG signal
is different from the controller common to the EEG signal and the respiratory signal,
since these sensors are monitoring different parts of our physiological system. In summary, the results indicate that the respiratory signal (THO) and the EMG signal lie on
different manifolds, which contain different physiological information with respect to
the EEG signal, or more generally, with respect to the brain activity.
To validate the consistency of the results and the above statement, in Figure 6 we
repeat the experiment and apply the filtering to a different EEG channel (O1A2), which
is different from O2A1, and as expected we observe similar results. In Figure 6 as well,
the two applications of the manifold filtering of the EEG channel through the two respiratory signals yield a similar representation, whereas the manifold filtering of the
EEG channel through the EMG signal yields a different representation. This again implies the existence of a common manifold hosting the common controller underlying
the two different respiratory signals associated with the respiratory system. In addition, this common controller is different from the common controller of the EEG and
the EMG signals. We report that the same procedure was applied to the same signal
recordings from 10 different subjects/cases and the results were consistent.
The application to sleep data demonstrates the potential of AD in filtering. More
results illustrating the advantage of AD in obtaining signal representation that is well
correlated with the sleep stage are presented in [31]. Further details, discussions, and
the presentation of the full experimental study are beyond the scope of this paper and
will appear in a future publication dedicated to this sleep application.
7. Extending the Common Manifold Representation
One of the key questions in manifold learning is how to extend the manifold representation to new data samples, once the manifold representation has been built from
an initial sample set, without applying the entire, usually computationally demanding
construction procedure [7, 26]. Since many of the manifold learning representations
are given by the eigenvectors of a particular kernel, the Nyström extension is typically
used to achieve this goal [41, 18]. Consider data on a manifold M with probability
measure ν(x), and a pairwise affinity kernel P . Assume that the kernel has an EVD,
Z
1
P (x, x0 )ψk (x0 )dν(x0 ), x ∈ M,
(61)
ψk (x) =
λk M
where ψk and λk are the eigenfunctions and eigenvalues of K, respectively. We can
see from (61) that one of the properties of the EVD is that the eigenfunction at every
point can be written as a (nontrivial) linear combination of the eigenfunction values.
Nyström extension exploits this property; given a new data sample x̃, either on M or
close to M, each eigenfunction at that point is therefore approximated by
Z
1
ψk (x̃) =
P (x̃, x0 )ψk (x0 )dν(x0 ).
(62)
λk M

30

While the Nyström extension is usually applied in order to reduce the number of computations stemming from the EVD, in the context of AD it gives rise to an important
benefit in addition to being more efficient computationally. In the context of AD, we are
given two data sets on two observable manifolds, which are functions of realizations of
triplets of hidden variables (x, y, x) from M × N1 × N2 with joint measure ν(x, y, z).
Recall that we denote by (x, y, z), (x0 , y 0 , z 0 ), (x00 , y 00 , z 00 ) three realizations of the hidden variables, by s1 = s(1) (x, y, z), s01 = s(1) (x0 , y 0 , z 0 ), s001 = s(1) (x00 , y 00 , z 00 ) their
corresponding samples in the observable space S (1) , and by s2 = s(2) (x, y, z), s02 =
s(2) (x0 , y 0 , z 0 ), s002 = s(2) (x00 , y 00 , z 00 ) their corresponding samples in the the observable
space S (2) . By (7), the AD kernel is written as a composition of two kernels
P ((x, y, z), (x00 , y 00 , z 00 ))
(63)
Z
=
P (2) ((x, y, z), (x0 , y 0 , z 0 ))P (1) ((x0 , y 0 , z 0 ), (x00 , y 00 , z 00 ))dν(x0 , y 0 , z 0 ),
M×N1 ×N2

where, by (4), P (1) is computed based on the observations s01 and s001 , and P (2) is
computed based on the observations s2 and s02 .
Substituting (63) into the Nyström extension (62) yields
Z Z
1
P (2) ((x̃, ỹ, z̃), (x0 , y 0 , z 0 ))P (1) ((x0 , y 0 , z 0 ), (x00 , y 00 , z 00 ))
ψk (x̃, ỹ, z̃) =
λk
M×N1 ×N2
(64)
× dν(x0 , y 0 , z 0 )ψk (x00 , y 00 , z 00 )dν(x00 , y 00 , z 00 )
Thus, in order to get an approximation of the eigenfunction of the AD kernel at a new
sample (x̃, ỹ, z̃), we only need to compute the kernel values P (2) ((x̃, ỹ, z̃), (x0 , y 0 , z 0 )),
i.e., the kernel between the new sample (x̃, ỹ, z̃) in the observable space S (2) and all
the existing samples (x0 , y 0 , z 0 ) in S (2) . In other words, the extension only requires
a new sample from only one of the sensors, which is an important benefit in many
applications, in particular, when one sensor is more difficult to obtain than the other.
We demonstrate the extension on the sleep application presented in Section 6. We
randomly take only 75% of the samples from an EEG recording (O2A1) and an EMG
recording. In the initialization stage, we construct the common manifold representation
based on these samples. Then, in an extension stage, we recover the common manifold
representation for the remaining 25% from solely the EMG samples. In Figure 7 we
present both the true manifold representation and the extension we obtain. We observe
that after the initial stage, which requires both measurements, we can accurately recover the common manifold representation from only one of the measurements; that
is, we do not have to evaluate the extension for two kernels, but rather just for one.
Specifically, in the sleep application, this result may have practical implications since
EEG is considered the most informative recording, yet it is also difficult to collect and
usually requires the help of a specialized technician. In addition, the filtering of the
EEG recording via the EMG recording enables us to clean the EEG signal from artifacts and noise and to show high correspondence with the sleep stage. The extension
result presented here implies that after a “calibration” stage, in which both recordings

31

are required, we can get access to the filtering result only based on the EMG recording,
which is easier to collect.
8. Conclusions
In this paper, we introduced a common manifold model underlying sensor observations for the purpose of fusion multimodal sensor data. To study this problem, we
proposed a method based on alternating diffusion and provided theoretical analysis under the common manifold model. Compared with traditional methods, this method is
able to capture the nonlinear manifold structure in the sense that both the topological
and geometrical structures are simultaneously preserved. Several applications were
provided, demonstrating the power of the proposed method and the analytic tools.
One important aspect of the current framework is its extension to more than two
sensors. Note that a straightforward extension of the proposed method to fuse data from
multiple sensors discovers the common manifold underlying all the sensors. Despite
being efficient and appropriate for certain purposes, in some sensor fusion applications
it might erase substantial information. In addition, the notion of nonlinear manifold
filtering, which is only introduced here, along with its filtering capabilities call for further analysis and more substantial mathematical foundations. For example, based on
the results presented here, in a future work, we will present a filtering scheme that
combines the common manifold representations obtained from all the possible pairs of
sensors in a multimodal experiment. This scheme enables us to combine relevant information from multimodal sensors in a pure data-driven manner, while filtering nuisance
phenomena measured only in single sensors, and thus, regarded as less important.
Another future directions will address applications and will focus on the following
three issues. First, we have only shown a proof-of-concept to the potential association
between the common manifold and physiological properties. A large scale study with
statistical analysis is needed to confirm the findings. Second, we may consider a natural generalization by taking the symmetrization into account. Based on the discussion
in Section 2.4, it is clear that the proposed AD scheme could provide a more accurate
information when the AD ends in a sensor with the nuisance that can be well modeled
by a manifold. The comparison of different schemes will be further explored. Third,
following the automatic annotation work in [49], we will investigate the correspondence between the obtained representations of the common manifolds, and the sleep
stage and we will devise improved automatic stage identification methods.
9. Acknowledge
The authors would like to thank Professor Ronald Coifman and Dr. Roy Lederman
for fruitful discussions, and Professor Yu-Lun Lo for sharing the sleep dataset. Hautieng Wu acknowledges the support of Sloan Research Fellow FR-2015-65363. Ronen Talmon was supported by the European Union’s Seventh Framework Programme
(FP7) under Marie Curie Grant 630657 and by the Israel Science Foundation (grant no.
1490/16).

32

References
[1] F. R. Bach and M. I. Jordan. Kernel independent component analysis. The Journal
of Machine Learning Research, 3:1–48, 2003.
[2] F. R. Bach, G. R.G. Lanckriet, and M. I. Jordan. Multiple kernel learning, conic
duality, and the smo algorithm. In Proceedings of the twenty-first international
conference on Machine learning, page 6. ACM, 2004.
[3] X. Bai, B. Wang, C. Yao, W. Liu, and Z. Tu. Co-transduction for shape retrieval.
IEEE Transactions on Image Processing, 21(5):2747–2757, 2012.
[4] M. Belkin and P. Niyogi. Laplacian Eigenmaps for Dimensionality Reduction
and Data Representation. Neural. Comput., 15(6):1373–1396, June 2003.
[5] B. Boots and G. J. Gordon. Two-manifold problems with applications to nonlinear
system identification. In Proc. 29th Intl. Conf. on Machine Learning (ICML),
2012.
[6] R. R. Coifman and S. Lafon. Diffusion maps. Appl. Comput. Harmon. Anal.,
21(1):5–30, 2006.
[7] R. R. Coifman and S. Lafon. Geometric harmonics: a novel tool for multiscale
out-of-sample extension of empirical functions. Applied and Computational Harmonic Analysis, 21(1):31–52, 2006.
[8] M. Davenport, C. Hegde, M. F. Duarte, and R. G. Baraniuk. Joint manifolds for
data fusion. IEEE Transactions on Image Processing, 19(10):2580–2594, 2010.
[9] A. M De Livera, R. J Hyndman, and R. D Snyder. Forecasting time series with
complex seasonal patterns using exponential smoothing. Journal of the American
Statistical Association, 106(496):1513–1527, 2011.
[10] V. R de Sa. Spectral clustering with two views. In ICML workshop on learning
with multiple views, 2005.
[11] V. R. de Sa, P. W. Gallagher, J. M. Lewis, and V. L. Malave. Multi-view kernel
construction. Machine Learning, 79(1-2):47–71, May 2010.
[12] M.P. do Carmo and F. Flaherty. Riemannian Geometry. Birkhauser Boston, 1992.
[13] D. L. Donoho, M. Gavish, and I. M. Johnstone. Optimal Shrinkage of Eigenvalues
in the Spiked Covariance Model. ArXiv e-prints, November 2013.
[14] D. L. Donoho and C. Grimes. Hessian eigenmaps: New locally linear embedding
techniques for high-dimensional data. Proc. Nat. Acad. Sci., 100:5591–5596,
2003.
[15] N. El Karoui. On information plus noise kernel random matrices. Ann. Stat.,
38(5):3191–3216, 2010.

33

[16] N. El Karoui and H.-T. Wu. Connection graph Laplacian methods can be made
robust to noise. Ann. Stat., 2014. Accepted for Publication.
[17] D. Eynard, K. Glashoff, M. M. Bronstein, and A. M. Bronstein. Multimodal diffusion geometry by joint diagonalization of laplacians. arXiv preprint
arXiv:1209.2295, 2012.
[18] C. Fowlkes, S. Belongie, F. Chung, and J. Malik. Spectral grouping using the
Nystrom method. Pattern Analysis and Machine Intelligence, IEEE Transactions
on, 26(2):214–225, 2004.
[19] P. G. Gould, A. B. Koehler, J. K. Ord, R. D. Snyder, R. J. Hyndman, and F. VahidAraghi. Forecasting time series with multiple seasonal patterns. European Journal of Operational Research, 191(1):207–222, 2008.
[20] A. Harvey and S. J. Koopman. Forecasting hourly electricity demand using timevarying splines. Journal of the American Statistical Association, 88(424):1228–
1236, 1993.
[21] A. Harvey, S. J. Koopman, and M. Riani. The modeling and seasonal adjustment
of weekly observations. Journal of Business & Economic Statistics, 15(3):354–
368, 1997.
[22] H. Hotelling. Relations between two sets of variates. Biometrika, 28(3/4):321,
December 1936.
[23] H.-C. Huang, Y.-Y. Chuang, and C.-S. Chen. Affinity aggregation for spectral
clustering. In Computer Vision and Pattern Recognition (CVPR), 2012 IEEE
Conference on, pages 773–780. IEEE, 2012.
[24] Y. Keller, R. Coifman, S. Lafon, and S. W. Zucker. Audio-visual group recognition using diffusion maps. IEEE Transactions on Signal Processing, 58(1):403–
413, January 2010.
[25] A. Kumar and H. Daumé. A co-training approach for multi-view spectral clustering. In Proceedings of the 28th International Conference on Machine Learning
(ICML-11), pages 393–400, 2011.
[26] S. Lafon, Y. Keller, and R. R. Coifman. Data fusion and multicue data matching
by diffusion maps. Pattern Analysis and Machine Intelligence, IEEE Transactions
on, 28(11):1784–1797, 2006.
[27] D. Lahat, T. Adali, and C. Jutten. Multimodal data fusion: An overview of methods, challenges, and prospects. Proceedings of the IEEE, 103(9):1449–1477,
2015.
[28] P. L. Lai and C. Fyfe. Kernel and nonlinear canonical correlation analysis. International Journal of Neural Systems, 10(05):365–377, 2000.

34

[29] G. R.G. Lanckriet, N. Cristianini, P. Bartlett, L. El Ghaoui, and M. I. Jordan.
Learning the kernel matrix with semidefinite programming. The Journal of Machine Learning Research, 5:27–72, 2004.
[30] R. R. Lederman and R. Talmon. Learning the geometry of common latent variables using alternating-diffusion. Appl. Comp. Harmon. Anal., 2015.
[31] R. R. Lederman, R. Talmon, H.-T. Wu, Y.-L. Lo, and R. R. Coifman. Alternating diffusion for common manifold learning with application to sleep stage
assessment. In IEEE International Conference on Acoustics, Speech and Signal
Processing (ICASSP), pages 5758–5762. IEEE, 2015.
[32] T. Lee-Chiong. Sleep Medicine: Essentials and Review. Oxford, 2008.
[33] O. Lindenbaum, A. Yeredor, and M. Salhov. Learning coupled embedding using
multiview diffusion maps. In Proc. Latent Variable Analysis and Signal Separation, pages 127–134. Springer, 2015.
[34] O. Lindenbaum, A. Yeredor, M. Salhov, and A. Averbuch. Multiview diffusion
maps. arXiv preprint arXiv:1508.05550, 2015.
[35] L. Luo, W. Jia, and C. Zhang. Mixed propagation for image retrieval. In IEEE
Fourth International Conference on Multimedia Information Networking and Security (MINES), pages 237–240, 2012.
[36] L. Luo, C. Shen, C. Zhang, and A. J. van den Hengel. Shape similarity analysis by self-tuning locally constrained mixed-diffusion. IEEE Transactions on
Multimedia, 15(5):1174–1183, 2013.
[37] N. F. Marshall and M. J. Hirn.
prints:1608.03628, 2016.

Time Coupled Diffusion Maps.

ArXiv e-

[38] T. Michaeli, W. Wang, and K. Livescu. Nonparametric canonical correlation analysis. arXiv preprint arXiv:1511.04839, 2015.
[39] B. Nadler, S. Lafon, R. R. Coifman, and I. G. Kevrekidis. Diffusion maps, spectral clustering and reaction coordinates of dynamical systems. Appl. Comput.
Harmon. Anal., 21(1):113–127, 2006.
[40] P. Niyogi, S. Smale, and S. Weinberger. Finding the Homology of Submanifolds
with High Confidence from Random Samples. Discrete & Computational Geometry, 23(April 2003):1–23, 2006.
[41] Evert Johannes Nyström. Über die praktische Auflösung von linearen Integralgleichungen mit Anwendungen auf Randwertaufgaben der Potentialtheorie.
Akademische Buchhandlung, 1929.
[42] S. T. Roweis and L. K. Saul. Nonlinear dimensionality reduction by locally linear
embedding. Science, 260:2323–2326, 2000.

35

[43] A. Singer and H.-T. Wu. Vector diffusion maps and the connection Laplacian.
Comm. Pure Appl. Math., 65(8):1067–1144, 2012.
[44] A. Singer and H.-T. Wu. Spectral convergence of the connection laplacian from
random samples. Information and Inference: A Journal of the IMA, 6(1):58–123,
2017.
[45] F. Takens. Detecting strange attractors in turbulence. In David Rand and LaiSang Young, editors, Dynamical Systems and Turbulence, volume 898 of Lecture
Notes in Mathematics, pages 366–381. Springer Berlin Heidelberg, 1981.
[46] J. W. Taylor. Short-term electricity demand forecasting using double seasonal
exponential smoothing. Journal of the Operational Research Society, 54(8):799–
805, 2003.
[47] J. B. Tenenbaum, V. de Silva, and J. C. Langford. A global geometric framework
for nonlinear dimensionality reduction. Science, 260:2319–2323, 2000.
[48] B. Wang, J. Jiang, W. Wang, Z.-H. Zhou, and Z. Tu. Unsupervised metric fusion
by cross diffusion. In Proc. IEEE Conference on Computer Vision and Pattern
Recognition (CVPR), pages 2997–3004, 2012.
[49] H.-T. Wu and Y.-L. Talmon, R.and Lo. Assess sleep stage by modern signal processing techniques. Biomedical Engineering, IEEE Transactions on, 62(4):1159–
1168, 2015.
[50] Z. Wu, Y. Wang, R. Shou, B. Chen, and X. Liu. Unsupervised co-segmentation
of 3d shapes via affinity aggregation spectral clustering. Computers & Graphics,
37(6):628–637, October 2013.
[51] Y. Zhou, X. Bai, W. Liu, and L. J Latecki. Fusion with diffusion for robust visual
tracking. In Advances in Neural Information Processing Systems, pages 2978–
2986, 2012.

36

12
20
10
15
8

6

10

4
5
2

0
0

0
0

0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01

Frequency

0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01

Frequency

(a)

(b)

18
12
16
10

14
12

8

10
6

8
6

4

4
2
2
0
0

0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01

Frequency

0
0

0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01

Frequency

(c)

(d)

Figure 3: The results of the application of AD to the toy problem. (a) The Fourier transform of the eigenvector associated with the nontrivial largest eigenvalue obtained from the diffusion maps application to the
reference data set Y. We observe that indeed the Fourier transform of the eigenvector approximates the function δ(ω−ω0 ) where ω0 = 1/1000. (b) The Fourier transform of the eigenvector associated with the largest
nontrivial eigenvalue obtained from the diffusion maps application to the data set X . We observe that indeed
the eigenvector is deformed and does not only consists of a single frequency ω0 = 1/1000. (c) The Fourier
transform of the eigenvector associated with the largest nontrivial eigenvalue obtained from AD application
starting from the reference set. We observe that indeed, as indicated by Theorem 3.3, the Fourier transform
of the eigenvector approximates the function δ(ω − ω0 ) where ω0 = 1/1000. (d) The Fourier transform of
the eigenvector associated with the largest nontrivial eigenvalue obtained from AD application starting from
the data set. In contrast to (c), and as indicated by Theorem 3.3, here the eigenvector is deformed and does
not solely consists of ω0 = 1/1000.

37

No. of Barrels [units of thousands]

10000
9500
9000
8500
8000
7500
7000
6500
6000

100

200

300

400

500

600

700

Time [days]

(a)
10

10
Seasonality Index
Ground Truth

9

6
5
4

7
6
5
4

3

3

2

2
10

20

30

40

Tested Period

(b)

50

60

8

Seasonality Index

7

1
50

Baseline Seasonality Index
Ground Truth

9

8

Seasonality Index

Seasonality Index

8

1

10
Seasonality Index
Ground Truth

9

7
6
5
4
3
2

51

52

53

Tested Period

(c)

54

55

1

10

20

30

40

50

60

Tested Period

(d)

Figure 4: The results of the application of Algorithm 2 to the time series of weekly U.S. finished motor
gasoline products supplied from February 1991 to July 2005 in units of thousands of barrels per day. (a)
The time series of number of barrels. (b) The seasonality index as a function of the tested frequency ω. (c)
A zoom into the range of interest around the true frequency (ω0 = 52.179). (d) The baseline seasonality
index as a function of the tested frequency.

38

Eigenvectors

Combined O2A1 THO
Combined O2A1 ABD
Only ABD
Only THO
Only O2A1

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Time [sec]

4

x 10

(a)

Eigenvectors

Combined O2A1 THO
Combined O2A1 EMG
Only EMG
Only THO
Only O2A1

0.2

0.4

0.6

0.8

1

Time [sec]

1.2

1.4

1.6

1.8

2
4

x 10

(b)
Figure 5: The leading nontrivial eigenvectors obtained from AD applied to EEG and respiratory signals as
a function of time. (a) The dashed lines represent the leading nontrivial eigenvectors obtained by diffusion
maps applied separately to an EEG recording from a single electrode (O2A1), a signal acquired by a motion
belt located on the chest (THO), and a signal acquired by a motion belt located on the abdomen (ABD). The
solid blue line is the leading nontrivial eigenvector resulting from AD applied to the EEG recording and the
signal from the motion belt on the chest. The solid green line is the leading nontrivial eigenvector resulting
from AD applied to the EEG recording and the signal from the motion belt on the abdomen. (b) The dashed
lines represent the leading nontrivial eigenvectors obtained by diffusion maps applied separately to an EEG
recording from a single electrode (O2A1), a signal acquired by a motion belt located on the chest (THO), and
an EMG recording measuring muscle movements. The solid blue line is the leading nontrivial eigenvector
resulting from AD applied to the EEG recording and the signal from the motion belt on the chest. The solid
green line is the leading nontrivial eigenvector resulting from AD applied to the EEG recording and the EMG
recording.

39

Eigenvectors

Combined O1A2 ABD
Combined O1A2 THO
Only THO
Only ABD
Only O1A2

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Time [sec]

4

x 10

(a)

Eigenvectors

Combined O1A2 ABD
Combined O1A2 EMG
Only EMG
Only ABD
Only O1A2

0.2

0.4

0.6

0.8

1

Time [sec]

1.2

1.4

1.6

1.8

2
4

x 10

(b)
Figure 6: The leading nontrivial eigenvectors obtained from AD applied to EEG and respiratory signals as a
function of time. Same as in Figure 5 only with a different EEG channel (O1A2).

40

0.04
Eigenvector
Extension

0.03

Eigenvector

0.02
0.01
0
−0.01
−0.02
−0.03
−0.04

0.5

1

1.5

2
2.5
Time [sec]

3

3.5

4
4

x 10

(a)
0.04
0.03

Extension

0.02
0.01
0
−0.01
−0.02
−0.03
−0.04
−0.03

−0.02

−0.01

0
0.01
Eigenvector

0.02

0.03

0.04

(b)
Figure 7: The extension of the leading nontrivial eigenvector of alternating diffusion applied to an EEG
recording (O2A1) and the EMG recording. In (a) we plot the true eigenvector and its extension as a function
of time, and in (b) we plot the scatter plot of the true eigenvector samples vs. the extended eigenvector
samples.

41

Appendix A. Proof of Theorems 3.2 and 3.3
Before proceeding, we need a discussion about the diffeomorphism Φ. Denote
Nx ι(i) (M) to be the (p − d)-dim subspace of Rp which contains all vectors perpendicular to dι(i) Tx M, where dι(i) is the total differential of the map ι(i) : M → Rp ,
and dι(i) Tx M := {dι(i) (v)| v ∈ Tx M}; that is, uT v = 0 for all u ∈ Nx ι(i) (M) and
v ∈ dι(i) Tx M. As M is compact and smoothly embedded in Rp , for i = 1, 2, we
could find a tubular neighborhood Tubτ (ι(i) (M)) ⊂ Rp of ι(i) (M), where τ > 0
is chosen so that Tubτ (ι(i) (M)) is homotopic to ι(i) (M). Then Φ could be extended to Tubτ (ι(1) (M)) and becomes a diffeomorphism between Tubτ (ι(1) (M)) and
Tubτ (ι(2) (M)), where we use the same notation to denote the extension. After extension, ∇Φ could be defined on Tubτ (ι(1) (M)). For x ∈ M, ∇Φ|ι(1) (x) is a linear
map from Tι(1) (x) Rp to Tι(2) (x) Rp , which maps dι(1) Tx M to dι(2) Tx M and maps
(i)

Nx ι(1) (M) to Nx ι(2) (M). Note that Πx (u, u) ∈ Nx ι(i) (M) for all u ∈ Tx M.
Recall the following two Lemmas. Below, we will adopt Einstein summation convention to simplify the notation.

Lemma A1. Fix i = 1, 2. In the normal coordinate around x ∈ M, when kvkg(i)  1,
v ∈ Tx M, the Riemannian measure satisfies


1
(i)
k l
3
Ric
(x)v
v
+
O(kvk
)
dv 1 ∧ dv 2 . . . ∧ dv n . (A.1)
dV (i) (exp(i)
v)
=
1
−
x
kl
6
In the polar coordinate v = tθ, where θ ∈ Tx M, kθkg(i) = 1, t > 0, we have


(i)
d−1
d+1
d+2
dV (i) (exp(i)
tθ)
=
t
+
t
Ric
(θ,
θ)
+
O(t
)
dtdθ.
x

(A.2)

Proof. See, for example, [43].
(1)

Lemma A2. Fix x ∈ M and y = expx (v), where v ∈ Tx M with kvkg(1)  1. We
have
(1)

(1)

ι(1) (y) = ι(1) (x) + dι(1) (v) + Q2 (v) + Q3 (v) + O(kvk4g(i) ),
(1)

(A.3)
(1)

where Π(1) is the second fundamental form of ι(1) , Q2 (v) := 21 Π(1) (v, v) and Q3 (v) :=
(1)
1 (1) (1)
(v, v). Further, for z = expx (u), where u ∈ Tx M with kukg(1)  1, we
6 ∇v Π
have
(1)

kι(1) (z) − ι(1) (y)k = ku − vkg(1) +

(1)

+

(1)

kQ2 (u) − Q2 (v)k2
2ku − vkg(1)
(1)

(A.4)

hdι(1) (u − v), Q3 (u) − Q3 (v)i
+ O(kuk4g(1) , kvk4g(1) ).
ku − vkg(1)

42

In particular, suppose v = tθ, where kθkg(1) = 1, we have
(1)

kι(1) (y) − ι(1) (x)k = t −

kQ2 (θ)k2 3
t + O(t4 ).
6

(A.5)

Proof. See, for example, [43] for the proof of (A.3). By (A.3) and a direct expansion,
(1)
we have (A.4). Note that (A.5) comes from (A.3) since kΠ(1) (v, v)k2 +hdι(1) (v), ∇v Π(1) (v, v)i =
0.
An immediate consequence of Lemma A2 is the distance between two points x
(1)
and y = expx (v), where v ∈ Tx M, when measured by the ambient metric; that is,
(1)
(1)
kι(1) (expx (v)) − ι(1) (x)k. We could call the distance kι(1) (expx (v)) − ι(1) (x)k
the ambient distance between x and y. Next, we discuss the difference between two
metrics g (1) and g (2) . Since g (i) is induced from (Rp , can) via ι(i) , and ι(i) , i = 1, 2,
is a smooth embedding of M into Rp , we could find a diffeomorphism Φ : Rp → Rp
so that Φ ◦ ι(1) = ι(2) . Hence, the metrics g (1) and g (2) are related by
kuk2g(2) = h∇Φ|ι(1) (x) dι(1) (u), ∇Φ|ι(1) (x) dι(1) (u)i = k∇Φ|ι(1) (x) dι(1) (u)k2 , (A.6)
where u ∈ Tx M. In the next lemma, we evaluate the ambient distance of two points
on M measured by a different metric.
Lemma A3. Fix x ∈ M. To simplify the notation, we ignore the subscription of
2
∇(1) Φ|ι(1) (x) and use ∇Φ. Similar simplification holds for ∇(1) Φ, dι(1) , Π(1) , etc.
(1)

Suppose y = expx u, where u ∈ Tx M and kukg(1) is small enough. Then we have
(2)

(2)

ι(2) (y) = ι(2) (x) + ∇Φdι(1) u + Q2 (u) + Q3 (u) + O(kuk4g(1) )
(2)

(A.7)

(2)

where Q2 and Q3 are quadratic and cubic polynomials respectively defined by
1
(2)
Q2 (u) := Π(2) (u, u),
2


1
(2)
(1)
2
(1)
Q3 (u) := ∇Φ∇(1)
Π
(u,
u)
+
∇
Φ
u,
Π
(u,
u)
.
u
6

(A.8)

(1)

Further, for z = expx v, where kvkg(1) is also small enough, we have
E
D
(2)
(2)
∇Φdι(1) (u − v), Q2 (u) − Q2 (v)

kι(2) (y) − ι(2) (z)k = k∇Φdι(1) (u − v)k +
k∇Φdι(1) (u − v)k
E
D
(2)
(2)
(2)
(2)
2 ∇Φdι(1) (u − v), Q3 (u) − Q3 (v) + kQ2 (u) − Q2 (v)k2
+
k∇Φdι(1) (u − v)k
+ O(kuk4g(1) , kvk4g(1) ).

43

In particular, suppose v = tθ, where kθkg(1) = 1, we have
D
E
(2)
∇Φdι(1) (θ), Q2 (θ)

(A.9)
kι(2) (y) − ι(2) (x)k = tk∇Φdι(1) (θ)k + t2
k∇Φdι(1) (θ)k
D
E
(2)
(2)
2 ∇Φdι(1) (θ), Q3 (θ) + kQ2 (θ)k2
3
+t
+ O(t4 ).
k∇Φdι(1) (θ)k
Proof. By Lemma A2,
(1)
ι(2) (exp(1)
(exp(1)
(A.10)
x u) = Φ(ι
x u))


1
1
= Φ ι(1) (x) + dι(1) u + Π(1) (u, u) + ∇(1)
Π(1) (u, u) + O(kuk4g(1) )
2
6 u


1 (1)
1 (1) (1)
(2)
(1)
4
= ι (x) + ∇Φ dι u + Π (u, u) + ∇u Π (u, u) + O(kukg(1) )
2
6


1 2
+ ∇ Φ (u, u) + ∇2 Φ u, Π(1) (u, u) + O(kuk4 )
2
(2)
(2)
(2)
= ι (x) + ∇Φdι(1) u + Q2 (u) + Q3 (u) + O(kuk4g(1) ),

where Π(2) (u, u) = ∇ΦΠ(1) (u, u) + ∇2 Φ (u, u) by the chain rule. Thus, we have
(2)
2
kι(2) (exp(1)
(exp(1)
x u) − ι
x v)k

(A.11)

2

= ∇Φdι(1) (u − v)
D
E
(2)
(2)
(2)
(2)
+ 2 ∇Φdι(1) (u − v), Q3 (u) − Q3 (v) + kQ2 (u) − Q2 (v)k2
+ O(kuk4g(1) ku − vkg(1) , kvk4g(1) ku − vkg(1) ).
Note that k∇Φdι(1) (θ)k is uniformly bounded below from zero and above for any
chosen x ∈ M. Thus, by taking the square root, we have (A.9).
Note that in (A.5) the geodesic distance between y and x, dι(1) (u) , is different
from the ambient distance with the error term of order 3, while in (A.9), the error
term is of order 2. This difference comes from the diffeomorphism Φ, in particular its
(2)
Hessian ∇2 Φ showing up in Q2 (u).
Lemma A4. Fix a symmetric matrix S ∈ Rd×d and w ∈ Rd . Then we have the
following integrations:
2
2
2 −1/2
π d/2
Swk2
e−kSuk e−ku−wk du = p
e−k(I+S )
2
det(I + S )
Rd

Z

(A.12)

Proof. Note that these integrations could be viewed as a convolution of two Gaussian
functions, so we could obtain the results by applying the Fourier transform.
We need the following technical lemmas to analyze the effective kernel K (e) . For
44

h ≥ 0, denote
(1)

(1)

B̃h (x) := exp(1)
x (Bh ),

(A.13)

(1)

where Bh = {u ∈ Tx M| kukg(1) ≤ h} ⊂ Tx M is a d-dim disk with the center 0
and the radius h.
Lemma A5. Suppose Assumptions (A1)-(A4) hold, F ∈ L∞ (M) and 0 < γ < 1/2.
(1)
Then, when  is small enough, for all pairs of x, x00 ∈ M so that x00 = expx v and
γ
kvkg(1) > 2 , the following holds:
Z

K̃(e2 ) (x, y)K̃(e1 ) (y, x00 )F (y)dV (1) (y) = O(d/2+3/2 ).

(A.14)

M

Proof. By the assumption that F ∈ L∞ (M), we immediately have
Z
K̃(2) (x, y)K̃(1) (y, x00 )F (x0 )dV (1) (y)

(A.15)

M

≤ kF kL∞

Z
(1)
B̃γ (x)

+ kF kL∞

K̃(2) (x, y)K̃(1) (y, x00 )F (x0 )dV (1) (y)

Z
(1)

K̃(2) (x, y)K̃(1) (y, x00 )F (x0 )dV (1) (y) .

M\B̃γ (x)

To bound the first integration in (A.15), note that the assumption d1 (x, x00 ) > 2γ
(1)
implies that d1 (y, x00 ) ≥ γ for all y ∈ B̃γ (x). Thus, by Lemma A2, we know that
(1)
(1) 00
when  is small enough, kι (y) − ι (x )kRp ≥ (1 + c)γ , where |c| < 1 depends on
the second fundamental form of ι(1) . Thus, when  is small enough, by the exponential
decay assumption of K (1) , we have

 (1)
2γ−1
kι (y) − ι(1) (x00 )k
√
≤ c2 e−c1 (1+c)
= O(d/2+3/2 )
(A.16)
K̃ (1)


 (1)

 (2)
(1)
(2)
(1)
(x00 )k
(y)k
(1) kι (y)−ι
√
√
K̃
is bounded
for all y ∈ B̃γ (x). Since K̃ (2) kι (x)−ι


by 1, we have
Z
(1)
B̃γ (x)

K̃(2) (x, y)K̃(1) (y, x00 )F (x0 )dV (1) (y) = O(d/2+3/2 ).

(A.17)

The second integration
could
 (1)in (A.15)
 be bounded directly by taking Lemma A3 and
(1)
(1)
(x00 )k
(1) kι (y)−ι
√
the fact that K̃
≤ 1 into account. Denote y = expx v so that


kvkg(1) ≥ γ . Note that since Φ is a diffeomorphism, the smallest singular value of
∇Φ|ι(1) (x) is bounded from below, say c0 > 0. Thus, when kvkg(1) ≥ γ and  is small

45

enough, we know
kι(2) (x) − ι(2) (y)k2

(A.18)

2

T

2

= k∇Φ|ι(1) (x) vk + 2(∇Φ|ι(1) (x) v) ∇ Φ|ι(1) (x) (v, v) +

O(kvk4g(1) )

≥ c0 (1 + c0 )kvk2g(1) /2.

Hence, when  is small enough, there exists |c0 | < 1 so that
K̃ (2)



kι(2) (x) − ι(2) (y)k
√




0

≤ c2 e−c1 c0 (1+c )

2γ−1

/2

= O(d/2+3/2 ),

(A.19)

K̃(2) (x, y)K̃(1) (y, x00 )F (x0 )dV (1) (y) = O(d/2+3/2 ),

(A.20)

and hence
Z
M\B̃1,γ (x)

which leads to the proof.
Lemma A6. Suppose Assumptions (A1)-(A4) hold, F ∈ C 3 (M) and 0 < γ < 1/2.
Take a pair of x, x00 ∈ M so that x00 = expx v, where v ∈ Tx M and kvkg(1) ≤ 2γ .
Fix normal coordinates around x associated with g (1) and g (2) , set
(1)
(2)
Rx = [d expx |0 ]−1 [dι(2) ]−1 ∇Φ[dι(1) ][d expx |0 ] : Rd → Rd and by the SVD Rx =
T
Ux Λx Vx , where Ux , Vx ∈ O(d). Then, when  is small enough, the following holds:
Z
K̃(2) (x, y)K̃(1) (y, x00 )F (y)dV (1) (y)
(A.21)
M
h
i
= d/2 F (x)A0 (v) + A2 (F, v) + O(3/2 ) ,
where
Z
A0 (v) :=
Rd

√ 
K̃ (2) (kRx wk) K̃ (1) kw − v/ k dw,

(A.22)

A0 (0) > 0, A0 decays exponentially and for i = 1, 2, A2 are defined in (A.37); both
A0 and A2 decay exponentially as kvk increases. Further, A0 (−v) = A0 (v) and
A2 (F, −v) = A2 (F, v).
2 √
When K (1) and K (2) are both Gaussian, that is, K (1) (t) = K (2) (t) = e−t / π,
we have
A0 (v) = p

π d/2

2 −1/2

det(I + Λ2x )

e−k(I+Λx )

Λx VxT vk2 /

.

(A.23)

This Lemma essentially says that the effective kernel associated with the AD actually still enjoys the exponential decay property of the kernel functions we favor. Further, if the chosen kernels are both Gaussian, then the effective kernel is Gaussian.
Note that due to the diffeomorphism, H (v) is not isotropic, no matter if the kernels

46

are Gaussian or not.
(i)
We remark more about the decomposition Rx = Ux Λx VxT . Note that d expx |0 :
T0 Tx M → Tx M, where i = 1, 2, are unitary and we could view T0 Tx M as Rd , so
we would not distinguish Tx M from Rd under the chosen normal coordinates. Thus,
we could apply the SVD on Rx and obtain Rx = Ux Λx VxT , where Ux , Vx ∈ O(d) and
Λx is a diagonal matrix with non-zero eigenvalues. Clearly, by the diffeomorphism
assumption, we know that the smallest eigenvalue of Λx is uniformly bounded for
x ∈ M away from 0 and from above.
Proof. Since M is compact and K is positive, we have
Z
K̃(2) (x, y)K̃(1) (y, x00 )F (y)dV (1) (y)

(A.24)

M

−
≤ kF kL∞

Z
(1)

K̃(2) (x, y)K̃(1) (y, x00 )F (y)dV (1) (y)

B̃γ (x)

Z
(1)

K̃(2) (x, y)K̃(1) (y, x00 )dV (1) (y) .

M\B̃γ (x)

(1)

Since K (y, x00 ) ≤ 1, by the same bound as that of the second integration of (A.15),
(A.24) is bounded by d/2+3/2 when  is small enough. Next, we handle the main term
we have interest:
Z
K̃(2) (x, y)K̃(1) (y, x00 )F (y)dV (1) (y)
(A.25)
(1)

B̃γ (x)

Z
=

(1)

K̃

(2)

B̃γ (x)



kι(2) (x) − ι(2) (y)k
√



K̃

(1)



kι(1) (y) − ι(1) (x00 )k
√




F (y)dV (1) (y).

(1)

Suppose y = expx u, where u ∈ Tx M with kukg(1)  1 and t ≥ 0. By (A.9) in
Lemma A3, we have
(2)

kι(2) (x) − ι(2) (y)k = k∇Φdι(1) (u)k + Q̃3 (u) + O(kuk4 ),

(A.26)

where
(2)

Q̃3 (u) :=
(2)

E
D
(2)
(2)
2 ∇Φdι(1) (u), Q3 (u) + kQ2 (u)k2
k∇Φdι(1) (u)k

(2)

.

(A.27)

Note that Q̃2 is an odd function and Q̃3 is an even function. Thus, by the Taylor

47

expansion, the first kernel function in (A.25) becomes

 (2)
kι (x) − ι(2) (y)k
√
K̃ (2)
(A.28)





 (2)

k∇Φdι(1) (u)k
k∇Φdι(1) (u)k Q̃3 (u)
kuk5
√
√
√
.
= K̃ (2)
+ [K̃ (2) ]0
+O




(1)

Similarly, by applying (A.4) in Lemma A2, when x00 = expx (v), where v ∈ Tx M
with kvkg(1)  1, we have
(1)

kι(1) (y) − ι(1) (x00 )k = kdι(1) (u − v)k + Q̃3 (u, v) + O(kuk4 , kvk4 ).

(A.29)

where
(1)

(1)

Q̃3 (u, v) =
(1)

(1)

(1)

(1)

kQ2 (u) − Q2 (v)k2 + 2hdι(1) (u − v), Q3 (u) − Q3 (v)i
. (A.30)
2kdι(1) (u − v)k
(1)

Note that Q̃3 (−u, −v) = Q̃3 (u, v). Thus, by the Taylor expansion, the second
kernel function in (A.25) becomes
 (1)

 (1)

kι (y) − ι(1) (x00 )k
kdι (u − v)k
√
√
K̃ (1)
= K̃ (1)
(A.31)



 (1)
 (1)

kuk4 kvk4
kdι (u − v)k Q̃3 (u, v)
(1) 0
√
√
+O √ , √
+ [K̃ ]
.




By the same argument as that of Lemma A5, we could replace the integral domain,
(1)
Bγ ⊂ Tx M, in (A.25) by Tx M with the error of order O(d/2+3/2 ), as  is small
enough. As a result, with the Taylor expansion of F , (A.25) becomes


 (2)



h
kRx uk Q̃3 (u)
kuk5 i
kRx uk
(2) 0
(2)
√
√
√
+ [K̃ ]
+O
K̃
(A.32)




Rd



 (1)


h
ku − vk
ku − vk Q̃3 (u, v)
kuk4 kvk4 i
√
√
√
× K̃ (1)
+ [K̃ (1) ]0
+O √ , √






Z

2

(1)


∇u,u F (x)
× F (x) + ∇(1)
+ O(kuk3 )
u F (x) +
2


(1)
i j
× 1 − Ricij (x)u u + O(kuk3 ) du + O(d/2+3/2 )


= d/2 F (x)A0, (v) + 1/2 A1, (F, v) + A2, (F, v) + O(3/2 ) ,

48

where



kRx uk
ku − vk
(1)
√
√
A0, (v) := 
K̃
K̃
du


Rd




Z
kRx uk
ku − vk
−d/2−1/2
(2)
(1)
√
√
A1, (F, v) := 
K̃
∇(1)
K̃
u F (x)du


Rd
(A.33)
−d/2

Z

(2)



A2, (F, v) := −d/2−1 [A21 (F, v) + A22 (F, v) + A23 (F, v) + A24 (F, v)]
and
 (2)


kRx uk Q̃3 (u) (1) ku − vk
√
√
√
K̃
du



Rd


 (1)

Z
ku − vk Q̃3 (u, v)
kRx uk
(1) 0
(2)
√
√
√
[K̃ ]
du (A.34)
A22 (F, v) := F (x)
K̃



Rd




Z
2
kRx uk
ku − vk
√
√
K̃ (2)
A23 (F, v) :=
K̃ (1)
∇(1)
u,u F (x)du


d
R




Z
kRx uk
ku − vk
(1)
(2)
(1)
√
√
K̃
K̃
A24 (F, v) := F (x)
Ricij (x)ui uj du.


d
R
Z

A21 (F, v) := F (x)

h

[K̃ (2) ]0



Here, we sort the terms according to the order of , and we√claim that Ai (F, v) is of
order O(1), for i = 0, 1, 2. By a change of variable w = u/ , we have
Z
√ 
A0, (v) =
K̃ (2) (kRx wk) K̃ (1) kw − v/ k dw.
(A.35)
Rd

By the assumption of the kernels, H(0) > 0 and H decays exponentially. Similarly,
we have
Z
√ 
A1, (F, v) =
K̃ (2) (kRx uk) K̃ (1) ku − v/ k ∇(1)
(A.36)
u F (x)du
Rd

and
√  (2)
(A.37)
[K̃ (2) ]0 (kRx uk) K̃ (1) ku − v/ k Q̃3 (u)du
Rd
Z
√  (1)
√
+ F (x)
K̃ (2) (kRx uk) [K̃ (1) ]0 ku − v/ k Q̃3 (u, v/ )du
Rd
Z
√ 
2
+
K̃ (2) (kRx uk) K̃ (1) ku − v/ k ∇(1)
u,u F (x)du
d
R
Z
√ 
(1)
+ F (x)
K̃ (2) (kRx uk) K̃ (1) ku − v/ k Ricij (x)ui uj du
Z

A2, (F, v) =F (x)

Rd

Note that since K̃ (1) and K̃ (2) both decay exponentially fast, we know that Ai, (F, v)
also decay exponentially fast as kvk increases.
49

(1)

(2)

By the symmetric properties of Q̃3 (u, v) and Q̃3 (u), A1, is anti-symmetric
associated with v and A2, (F, v) is symmetric associated with v; that is, A1, (F, −v) =
−A1, (F, v) and A2, (F, −v) = A2, (F, v).
2 √
Finally, when both kernels are Gaussian, that is, K̃ (1) (t) = K̃ (2) (t) = e−t / π,
by the fact that Rx = Ux Λx VxT and Lemma A4, the leading order term in (A.32)
Z
√ 2
2
T
A0, (v) =
e−kΛx uk e−ku−Vx v/ k du
(A.38)
Rd

=p

1
det(I +

2 −1/2

Λ2x )

e−k(I+Λx )

Λx VxT vk2 /

.

With Lemma A5 and Lemma A6, we can finish the study of the effective kernel
K (e) . Before proving Theorem 3.2, we mention that in the proof, we focus on how the
AD process behaves on (M, g (1) ) by converting most “measurements” on (M, g (2) )
back to (M, g (1) ). To do so, note that by Assumption 3.1 (A4), dµM = p1 dV (1) =
p2 dV (2) , (43) becomes
(e)

R
T f (x) =

M

(A.39)

(e)

R
=

K̃ (x, x00 )f (x00 )p1 (x00 )dV (1) (x00 )
R
(2)
K̃ (x, x̄)p2 (x̄)dV (2) (x̄)
M

M

K̃ (x, x00 )f (x00 )p1 (x00 )dV (1) (x00 )
R
(2)
K̃ (x, x̄)p1 (x̄)dV (1) (x̄)
M

and (45) becomes
K̃(e) (x, x00 )

(1)

(2)

K̃ (x, x0 )K̃ (x0 , x00 )

Z
=
M

R

(1)

K̃ (x0 , x̄)p1 (x̄)dV (1) (x̄)
M

M

R

(A.40)

p1 (x0 )dV (1) (x0 ).

(A.41)

(1)

(2)

K̃ (x, x0 )K̃ (x0 , x00 )

Z
=

p2 (x0 )dV (2) (x0 )

(1)

K̃ (x0 , x̄)p1 (x̄)dV (1) (x̄)
M

As we will see shortly, although the AD works on two different metrics, most quantities
could be converted to (M, g (1) ), except the estimated distance among two points.
(1)

Proof of Theorem 3.2. To finish the proof, we study (A.41). We focus on x00 = expx v
so that kvkg(1) ≤ 2γ . By the absolute continuity assumption of µM , we have
Z
p̃1, (x0 ) := −d/2
K̃(1) (x0 , x00 )p1 (x00 )dV (1) (x00 ),
(A.42)
M

which by the assumption that p1 ∈ C 4 (M) and a direct expansion becomes
p̃1, (x0 ) = p1 (x0 ) + Q(x0 ) + O(2 ),

50

(A.43)

where
(1)

µ2,0 h (1) 0
Q(x ) :=
s (x ) +
2
0

R
S d−1

i
kΠ(1) (θ, θ)kdθ
0
(1)
0
p
(x
)
−
∆
p
(x
)
,
1
1
24|S d−1 |

(A.44)

where s(1) (x0 ) is the scalar curvature of g (1) at x0 . See, for example, [43, Lemma B.10]
for a proof. Thus, we have
p1 (x0 )
= 1 − Q(x0 ) + O(3/2 ).
p̃1, (x0 )

(A.45)

Since g (1) and p1 are both C 4 , we know that Q ∈ C 2 (M). Hence, by (A.45) we have
 p

1
A1,
, v = A1, (1, v) + O() = O()
(A.46)
p̃1,
since A1, (1, v) = 0, and
A2,

 p

1
, v = A2, (1, v) + O().
p̃1,

(A.47)

By (A.45), the assumptions and Lemma A6, the effective kernel becomes
K̃(e) (x, x00 )
Z
(1)
(2)
−d/2 K̃ (x, x0 )K̃ (x0 , x00 )
=
p1 (x0 )dV (1) (x0 )
R
(1) 0
00 )p (x00 )dV (1) (x00 )
M −d/2
(x
,
x
K̃

1
M
Z
p1 (x0 )
K̃(2) (x, x0 )K̃(1) (x0 , x00 )
= −d/2
dV (1) (x0 )
p̃1, (x0 )
M

(A.48)

= A0, (v) + (A2, (1, v) + Q(x)A0, (v)) + O(3/2 ),

While the proof of Theorem 3.3 becomes standard once the effective kernel is
known, we provide the proof for the sake of self-containedness and to show the net
diffusion outcome on (M, g (1) ). We need the following Lemma, which describes the
asymptotical behavior of a kernel not in the normalization form.
Lemma A7. Suppose Assumption (A1)-(A4) hold. Fix x ∈ M and pick F ∈ C 3 (M).
Fix normal coordinates around x associated with g (1) and g (2) so that {Ei }di=1 ⊂
(2)
(1)
Tx M is orthonormal associated with g (1) . Set Rx = [d expx |0 ]−1 [dι(2) ]−1 ∇Φ[dι(1) ][d expx |0 ]
and by the SVD Rx = Ux Λx VxT , where Λx = diag[λ1 , . . . , λd ]. Then, when  is small
enough, we have
Z
F (x)
−d/2 K̃(e) (x, x00 )F (x00 )dV (1) (x00 ) =
+ C2 (F ) + O(3/2 ) (A.49)
det(Λ
)
x
M
where C2 (F ) is defined in (A.62).
51

(1)

Proof. We focus on x00 = expx v so that kvkg(1) ≤ 2γ . Note that by Lemma (A5),
we will ignore the integration over M\B̃γ (x) with the error of order d/2+3/2 . By
the exponential decay of A0, (v) and A2, (v) and the finite volume assumption of M,
when  is small enough, by Theorem 3.2, (A.49) becomes
Z
−d/2 K̃(e) (x, x00 )F (x00 )dV (1) (x00 )
(A.50)
M
Z
=
−d/2 K̃(e) (x, x00 )F (x00 )dV (1) (x00 ) + O(3/2 )
B̃γ (x)
Z
= −d/2
[A0, (v) + (A2, (1, v) + Q(x)A0, (v))]
Bγ



1 (1) 2
× F (x) + ∇(1)
v F (x) + ∇v,v F (x))
2


(1)
× 1 + Ricij v i v j dv + O(3/2 ).
Expand the multiplication and sort terms in the increasing order, we have to evaluate
the following integrations to obtain (A.50):
Z
−d/2 F (x)
A0, (v)dv
(A.51)
Bγ
Z
−d/2+1/2
A0, (v)∇(1)
(A.52)
v F (x)dv
Bγ
Z
2
1
(1) i j 
−d/2
A0, (v) ∇(1)
(A.53)
v,v F (x) + F (x)Ricij v v dv
2
Bγ
Z
+ −d/2+1 F (x)
(A2, (v) + Q(x)A0, (v))dv.
Bγ

Here, (A.51) is the 0-th order term, (A.52) is the 1-st order term and (A.53) is the 2-nd
order term. By the same argument as the above, we could replace the integral domain
Bγ by Rd in (A.51), (A.52) and (A.53), with a higher order error. By plugging A0 (v)
and changing the integration order, (A.51) becomes
Z
hZ
i
F (x)
K̃ (2) (kRx wk)
K̃ (1) (kw − vk) dv dw + O(3/2 ),
(A.54)
Rd

Rd

which by a direct calculation becomes
F (x)
+ O(3/2 ).
det(Λx )
(1)

(2)

(A.55)

Here, we use the assumption that µ0,0 = µ0,0 = 1, change of variable of u = Rx w
and the SVD Rx = Ux Λx VxT . The first order term (A.52) become 0 due to the anti(1)
symmetry of ∇v F (x) associated with v.

52

To evaluate (A.53), note that by the change of variable we have
Z
2
1
(1) i j 
−d/2
A0, (v) ∇(1)

v,v F (x) + F (x)Ricij v v dv
2
Bγ
Z

1 (1) 2
(1)
=
A0,1 (v) ∇v,v F (x) + F (x)Ricij v i v j dv + O(3/2 )
2
Rd

(A.56)

and


Z

−d/2+1

F (x)
(A2, (v) + Q(x)A0, (v))dv
Bγ
Z
= F (x)
(A2,1 (v) + Q(x)A0,1 (v))dv + O(3/2 ).

(A.57)

Rd
(1) 2

2

By noting that ∇v,v F (x) = v T ∇(1) F (x)v, to simplify (A.56), we denote Sx :=
1 (1) 2
F (x) + F (x)Ric(1) (x), which is a symmetric matrix. By a change of variable
2∇
u = w − v, we have
Z
2
1
(1) i j 
A0,1 (v) ∇(1)
(A.58)
v,v F (x) + F (x)Ricij v v dv
2
d
Z
ZR
i
h
K̃ (1) (kw − vk) v T Sx vdv dw
=
K̃ (2) (kRx wk)
d
Rd
ZR
i
hZ
K̃ (1) (kuk) (w − u)T Sx (w − u)du dw
K̃ (2) (kRx wk)
=
d
Rd
ZR
i
hZ
K̃ (1) (kuk) wT Sx wdu dw
K̃ (2) (kRx wk)
=
Rd
Rd
Z
i
hZ
(2)
K̃ (1) (kuk) uT Sx udu dw
+
K̃ (kRx wk)
Rd
Rd
Z
Z
1
1
=
K̃ (2) (kRx wk) wT Sx wdw +
K̃ (1) (kuk) uT Sx udu
d Rd
d det(Λx ) Rd
where the third equality holds since the crossover term
ZZ
K̃ (2) (kRx wk) K̃ (1) (kuk) uT Swdudw = 0
due to the symmetry of S d−1 , and the last equality holds due to

53

R
Rd

(A.59)

K̃ (2) (kRx wk) dw =

1
det(Λx ) .

Note that
Z

K̃ (2) (kRx wk) wT Swdw

(A.60)

Rd
(2)

(2)

µ2,0
µ2,0
tr(Λx VxT Sx Vx Λx ) =
tr(Vx Λx VxT Sx Vx Λx VxT )
d det(Λx )
d det(Λx )
(2)
d


µ2,0  1 X
2
(1)
λi ∇(1) Ei ,Ei F (x) + Ricii (x)F (x)
=
d det(Λx ) 2 i=1

=

and
Z

(1)

K̃ (1) (kuk) uT Sudu =

Rd


µ2,0 1 (1)
∆ F (x) + s(1) (x)F (x) ,
d 2

(A.61)

Pd
Pd
2
(1)
since ∆(1) F (x) = i=1 ∇(1) Ei ,Ei F (x) and i=1 Ricii (x) = s(1) (x). Hence, we
have obtained a simplification of (A.56).
As a result, by denoting


2
(1)
λi ∇(1) Ei ,Ei F (x) + Ricii (x)F (x)
C2 (F ) :=
d2 det(Λx )

(1)
µ2,0 12 ∆(1) F (x) + s(1) (x)F (x)
+
d2 det(Λx )
Z
+ F (x)
(A2,1 (v) + Q(x)A0,1 (v))dv,
(2) 1
2

µ2,0

Pd

i=1

(A.62)

Rd

we conclude
that (A.53) becomes C2 (F ) + O(3/2 ). Note that we do not simplify
R
F (x) Rd (A2,1 (v) + Q(x)A0,1 (v))dv since this term will be eliminated eventually.
With the above, we could finish the proof of Theorem 3.3. We mention that we
could see from the proof that the net outcome of the AD, while two metrics are involved, eventually could be viewed as a “deformed” diffusion solely with the single
metric g (1) .
Proof of Theorem 3.3. The proof is based on Lemma A7. The numerator of D(e) f (x)
is exactly the same as that stated in Lemma A7; that is,
Z
−d/2 K̃(e) (x, x00 )f (x00 )p1 (x00 )dV (1) (x00 )
(A.63)
M

p1 (x)f (x)
=
+ C2 (f p1 ) + O(3/2 ).
det(Λx )
The denominator of D(e) f (x) is evaluated by plugging f = 1 into Lemma A7; that is,

54

we have
Z
p1 (x)
−d/2 K̃(e) (x, x00 )p1 (x00 )dV (1) (x00 ) =
+ C2 (p1 ) + O(3/2 ). (A.64)
det(Λ
)
x
M
Finally, by the binomial expansion of the denominator, we obtain the result; that is,
(e)

K̃ (x, x00 )f (x00 )p1 (x00 )dV (1) (x00 )
R
(e)
K̃ (x, x00 )p1 (x00 )dV (1) (x00 )
M
i
det(Λx ) h
C2 (f p1 ) − f (x)C2 (p1 ) + O(3/2 ).
= f (x) + 
p1 (x)
R

M

(A.65)

By a direct expansion, we have
(2)

C2 (f p1 ) − f (x)C2 (p1 ) :=

µ2,0

Pd

i=1



2
(1)
(1)
λi p1 (x)∇(1) Ei ,Ei f (x) + 2∇Ei f (x)∇Ei p1 (x)
2d2 det(Λx )

(1)

+

µ2,0 [p1 (x)∆(1) f (x) + 2∇(1) f (x)∇(1) p1 (x)]
(A.66)
2d2 det(Λx )

and hence the conclusion
R
(e)
K̃ (x, x00 )f (x00 )p1 (x00 )dV (1) (x00 )
M
R
(e)
K̃ (x, x00 )p1 (x00 )dV (1) (x00 )
M

(A.67)

= f (x) +

(2) d
(1)
(1)
 (1) 2
µ2,0 X
2∇Ei f (x)∇Ei p1 (x) 
λ
∇
f
(x)
+
i
Ei ,Ei
2d2 i=1
p1 (x)

+

µ2,0 h (1)
2∇(1) f (x)∇(1) p1 (x) i
∆
f
(x)
+
+ O(3/2 ).
2d2
p1 (x)

(1)

Finally, we show the proof of Corollary 3.3.
Proof of Corollary 3.3. Combining (36), Theorem 3.3, and the assumption of νM , for
every fixed (y, z), we have
Z hZ
i
Df (x) =
P(N2 ,z) (x, x0 )P(N1 ,y) (x0 , x00 )dνM (x0 ) Ef (x00 )dνM (x00 )
M

M

(A.68)
= Ef (x) +

(2)
d
µ2,0 (z) h X
λi
2d2
i=1

2

∇(1) Ei ,Ei Ef (x) +

(1)

+

(1)
(1)
2∇Ei Ef (x)∇Ei p1 (x) i

p1 (x)

µ2,0 (z) h (1)
2∇(1) Ef (x)∇(1) p1 (x) i
∆
Ef
(x)
+
+ O(3/2 ),
2d2
p1 (x)

55

where
Z

(i)

µ2,0 (z) =

Rd

kyk2 P̃ (Oi ,z) (kyk)dy.

(A.69)

By the commutativity of E and D, applying the effective operator E bilaterally yields
2∇(2) Ef (x)∇(2) p1 (x) i
C2 (x) h (2)
∆
Ef
(x)
+
(A.70)
2d2
p1 (x)
C1 (x) h (1)
2∇(1) Ef (x)∇(1) p1 (x) i
+
∆ Ef (x) +
+ O(3/2 ),
2
2d
p1 (x)

EDf (x) = Ef (x) +

where Ci (x) =

(i)

R

µ2,0 (z)dνNi |M (z|x). To finish the proof, we show that the depen-

Ni

(1)

(1)

(1)

dence of Ci (x) on x appears on the higher order term. Denote B̃h (x) := expx (Bh ),
(1)
where Bh = {u ∈ Tx M| kukg(1) ≤ h} ⊂ Tx M is a d-dim disk with the center 0
and the radius h > 0. Thus, by the same arguments as above, we have
Z
(i)
µ2,0 (z)dνNi |M (z|x)
(A.71)
Ni
Z Z
=
kyk2 P̃ (Oi ,z) (kyk)dydνNi |M (z|x)
d
ZNi ZR Z
q

d2g(i) (y, x)P̃ (i)
=
d2g(i) (y, x) + d2Ni (z, z 0 ) dV (y)dνNi |M (z 0 |x)dνNi |M (z|x) + O()
B̃ (1/2 )
N
N
Z i i x

i
q
hZ Z
P̃ (i)
d2g(i) (y, x) + d2Ni (z, z 0 ) dνNi |M (z 0 |x)dνNi |M (z|x) dV (y) + O().
=
d2g(i) (y, x)
B̃x (1/4 )

Ni

Ni

By denoting

P̃˜ (i) dg(i) (y, x) :=

Z

Z

Ni

P̃ (i)

Ni

q


d2g(i) (y, x) + d2Ni (z, z 0 ) dνNi |M (z 0 |x)dνNi |M (z|x),

we have
Z

(i)

Ni

µ2,0 (z)dνNi |M (z|x)

Z
=
B̃x (1/4 )

Z
=
Bx (1/4 )

Z
=
Rd

By denoting Ci :=


d2g(i) (y, x)P̃˜ (i) dg(i) (y, x) dV (y) + O()
kyk2 P̃˜ (i) (kyk) dy + O()

kyk2 P̃˜ (i) (kyk) dy + O().

R
Rd

kyk2 P̃˜ (i) (kyk) dy, we hence obtain the conclusion.

56

(A.72)

