Geometric diffusions as a tool for harmonic analysis

and structure definition of data: Diffusion maps

R. R. Coifman*†, S. Lafon*, A. B. Lee*, M. Maggioni*, B. Nadler*, F. Warner*, and S. W. Zucker‡
*Department of Mathematics, Program in Applied Mathematics, Yale University, 10 Hillhouse Avenue, New Haven, CT 06510; and ‡Department
of Computer Science, Yale University, 51 Prospect Street, New Haven, CT 06510

Contributed by R. R. Coifman, January 13, 2005

We provide a framework for structural multiscale geometric orga- extend the scope of application of these ideas to various tasks, such
nization of graphs and subsets of ⺢n. We use diffusion semigroups as regression of empirical functions, by adjusting the infinitesimal
to generate multiscale geometries in order to organize and repre- descriptions, and the description of the long-time asymptotics of
sent complex structures. We show that appropriately selected stochastic dynamical systems.
eigenfunctions or scaling functions of Markov matrices, which The simplest way to introduce our approach is to consider a
describe local transitions, lead to macroscopic descriptions at set X of normalized data points. Define the ‘‘quantized’’ corre-
different scales. The process of iterating or diffusing the Markov lation matrix C ⫽ {c ij}, where c ij ⫽ 1 if (x i䡠x j) ⬎ 0.95, and c ij ⫽
matrix is seen as a generalization of some aspects of the Newto- 0 otherwise. We view this matrix as the adjacency matrix of a
nian paradigm, in which local infinitesimal transitions of a system graph on which we define an appropriate Markov process to start
lead to global macroscopic descriptions by integration. We provide our analysis. A more continuous kernel version can be defined
a unified view of ideas from data analysis, machine learning, and 2
as c ij ⫽ e (1⫺(xi䡠xj)/␧) ⫽ e ⫺(储xi⫺xj储 /2␧). The remarkable fact is that the
numerical analysis.
eigenvectors of this ‘‘corrected correlation’’ can be much more
meaningful in the analysis of data than the usual principal
he geometric organization of graphs and data sets in ⺢n is a
T central problem in statistical data analysis. In the continuous
Euclidean setting, tools from harmonic analysis, such as Fourier
components as they relate to diffusion and inference on the data.
As an illustration of the geometric approach, suppose that the
data points are uniformly distributed on a manifold X. Then it is
decompositions, wavelets, and spectral analysis of pseudo- known from spectral graph theory (8) that if W ⫽ {w ij} is any
differential operators, have proven highly successful in many areas symmetric positive semi-definite matrix, with nonnegative en-
such as compression, denoising, and density estimation (1, 2). In this
tries, then the minimization of

冘
paper, we extend multiscale harmonic analysis to discrete graphs
and subsets of ⺢n. We use diffusion semigroups to define and Q共f兲 ⫽ wij共fi ⫺ fj兲2,
generate multiscale geometries of complex structures. This frame-
i, j
work generalizes some aspects of the Newtonian paradigm, in which
local infinitesimal transitions of a system lead to global macroscopic where f is a function on the data set X with the additional
descriptions by integration, the global functions being characterized constraint of unit norm, is equivalent to finding the eigenvec-
by differential equations. We show that appropriately selected tors of D⫺1/2WD1/2, where D ⫽ {dij} is a diagonal matrix with
eigenfunctions of Markov matrices (describing local transitions, or diagonal entry dii equal to the sum of the elements of W along
affinities in the system) lead to macroscopic representations at 2
the ith row. Belkin et al. (5) suggest the choice wij ⫽ e⫺(储xi⫺ xj储 /␧),
different scales. In particular, the top eigenfunctions permit a in which case the distortion Q clearly penalizes pairs of points
low-dimensional geometric embedding of the set into ⺢k, with k ⬍⬍ that are very close, forcing them to be mapped to very close
n, so that the ordinary Euclidean distance in the embedding space values by f. Likewise, pairs of points that are far away from
measures intrinsic diffusion metrics on the data. Many of these each other play no role in this minimization. The first few
ideas appear in a variety of contexts of data analysis, such as spectral
eigenfunctions {␾k} are then used to map the data in a
graph theory, manifold learning, nonlinear principal components,
nonlinear way so that the closeness of points is preserved. We
and kernel methods. We augment these approaches by showing that
will provide a principled geometric approach for the selection
the diffusion distance is a key intrinsic geometric quantity linking
spectral theory of the Markov process, Laplace operators, or of eigenfunction coordinates.
kernels, to the corresponding geometry and density of the data. This general framework based upon diffusion processes leads to
This opens the door to the application of methods from numerical efficient multiscale analysis of data sets for which we have a
analysis and signal processing to the analysis of functions and Heisenberg localization principle relating localization in data to
transformations of the data. localization in spectrum. We also show that spectral properties can
be employed to embed the data into a Euclidean space via a
Diffusion Maps diffusion map. In this space, the data points are reorganized in such
The problem of finding meaningful structures and geometric de- a way that the Euclidean distance corresponds to a diffusion metric.
scriptions of a data set X is often tied to that of dimensionality The case of submanifolds of ⺢n is studied in greater detail, and we
reduction. Among the different techniques developed, particular show how to define different kinds of diffusions to recover the
attention has been paid to kernel methods (3). Their nonlinearity intrinsic geometric structure, separating geometry from statistics.
as well as their locality-preserving property are generally viewed as More details on the topics covered in this section can be found in
a major advantage over classical methods like principal component ref. 9. We also propose an additional diffusion map based on a
analysis and classical multidimensional scaling. Several other meth- specific anisotropic kernel whose eigenfunctions capture the long-
ods to achieve dimensional reduction have also emerged from the time asymptotics of data sampled from a stochastic dynamical
field of manifold learning, e.g., local linear embedding (4), Lapla- system (10).
cian eigenmaps (5), Hessian eigenmaps (6), local tangent space
alignment (7). All these techniques minimize a quadratic distortion
measure of the desired coordinates on the data, naturally leading to †To whom correspondence should be addressed. E-mail: coifman-ronald@[Link].
the eigenfunctions of Laplace-type operators as minimizers. We © 2005 by The National Academy of Sciences of the USA

7426 –7431 兩 PNAS 兩 May 24, 2005 兩 vol. 102 兩 no. 21 [Link]兾cgi兾doi兾10.1073兾pnas.0500334102
Construction of the Diffusion Map. From the above discussion, the 2
Dm 共x, y兲 ⫽ ã共m兲共x, x兲 ⫹ ã共m兲共y, y兲 ⫺ 2ã共m兲共x, y兲.
data points can be thought of as being the nodes of a graph whose
weight function k(x, y) (also referred to as ‘‘kernel’’ or ‘‘affinity The quantity a(x, y), which is related to ã(x, y) according to Eq.
function’’) satisfies the following properties: 2, can be interpreted as the transition probability of a diffusion
process, while a(m)(x, y) represents the probability of transition from
Y k is symmetric: k(x, y) ⫽ k(y, x), x to y in m steps. To this diffusion process corresponds the distance
Y k is positivity preserving: for all x and y in X, k(x, y) ⱖ 0, Dm(x, y), which defines a metric on the data that measures the rate
Y k is positive semi-definite: for all real-valued bounded func- of connectivity of the points x and y by paths of length m in the data,
tions f defined on X, and, in particular, it is small if there are a large number of paths

冕冕
connecting x and y. Note that, unlike the geodesic distance, this
k共x, y兲f共x兲f共y兲d␮共x兲d␮共y兲 ⱖ 0, metric is robust to perturbations on the data.
The dual point of view is that of the analysis of functions
X X
defined on the data. The kernel ã (m)(x, 䡠) can be viewed as a
where ␮ is a probability measure on X. bump function centered at x that becomes wider as m increases.
The distance D 2m(x, y) is also a distance between the two bumps
The construction of a diffusion process on the graph is a ã (m)(x, 䡠) and ã (m)(y, 䡠):

冕
classical topic in spectral graph theory [weighted graph Lapla-
cian normalization (8)], and the procedure consists in renor-
malizing the kernel k(x, y) as follows: for all x 僆 X,
2
D2m 共x, y兲 ⫽ 兩ã共m兲共x, z兲 ⫺ ã共m兲共y, z兲兩2dz.

冕
X

let v共x兲 ⫽ k共x, y兲d ␮ 共y兲, The eigenfunctions have the classical interpretation of an orthonor-
X
mal basis, and their frequency content can be related to the
spectrum of operator A in what constitutes a generalized Heisenberg

MATHEMATICS
APPLIED
and set principle. The key observation is that, for many practical examples,
k共x, y兲 the numerical rank of the operator A(m) decays rapidly as seen from
a共x, y兲 ⫽ . Eq. 4 or from Fig. 1. More precisely, since 0 ⱕ ␭i ⱕ ␭0 ⫽ 1, the
v共x兲
kernel ã(m)(x, y), and therefore the distance Dm(x, y), can be
Notice that we have the following conservation property: computed to high accuracy with only a few terms in the sum of 4,
that is to say, by only retaining the eigenfunctions ␾i for which ␭i2m

冕 X
a共x, y兲d␮共y兲 ⫽ 1. [1]
exceeds a certain precision threshold. Therefore, the rows (the
so-called bumps) of Am span a space of lower numerical dimension,
and the set of columns can be down-sampled. Furthermore, to
generate this space, one just needs the top eigenfunctions, as
Therefore, the quantity a(x, y) can be viewed as the probability prescribed in Eq. 4. Consequently, by a change of basis, eigenfunc-
for a random walker on X to make a step from x to y. Now we tions corresponding to eigenvalues at the beginning of the spectrum
naturally define the diffusion operator have low frequencies, and the number of oscillations increase as one

冕
moves further down in the spectrum.
The link between diffusion maps and distances can be sum-
Af共x兲 ⫽ a共x, y兲f共y兲d␮共y兲. marized by the spectral identity
X

As is well known in spectral graph theory (8), there is a spectral

储⌽m共x兲 ⫺ ⌽m共y兲储2 ⫽ 冘
jⱖ0
␭j2m共␾j共x兲 ⫺ ␾j共y兲兲2 ⫽ Dm
2
共x, y兲,

theory for this Markov chain, and if à is the integral operator

defined on L 2(X) with the kernel which means that the diffusion map ⌽ m embeds the data into a

冑
Euclidean space in which the Euclidean distance is equal to the
v共x兲 diffusion distance D m. Moreover, the diffusion distance can be
ã共x, y兲 ⫽ a共x, y兲 , [2] accurately approximated by retaining only the terms for which
v共y兲
␭ j2m remains numerically significant: the embedding
then it can be verified that à is a symmetric operator. Conse-
quently, we have the following spectral decomposition x 哫 x̆ ⫽ 共␭m
0 ␾ 0共x兲, ␭ 1 ␾ 1共x兲, . . . , ␭ j0 ␾ j0共x兲兲
m m

ã共x, y兲 ⫽ 冘 ␭i2␾i共x兲␾i共y兲, [3] satisfies

冘
j0⫺1
iⱖ0
2
Dm 共x, y兲 ⫽ ␭j2m共␾j共x兲 ⫺ ␾j共y兲兲2共1 ⫹ O共e⫺␣m兲兲
where ␭ 0 ⫽ 1 ⱖ ␭ 1 ⱖ ␭ 2 ⱖ . . . . Let ã (m)(x, y) be the kernel of j⫽0
à m. Then we have

ã共m兲共x, y兲 ⫽ 冘
iⱖ0
␭i2m␾i共x兲␾i共y兲. [4]
⫽ 储x̆ ⫺ y̆储2共1 ⫹ O共e⫺␣m兲兲.

Therefore, there exists an m 0 such that for all m ⱖ m 0, the

diffusion map with the first j 0 eigenfunctions embeds the data
Lastly, we introduce the family of diffusion maps {⌽ m} by into ⺢ j0 in an approximately isometric fashion, with respect to the

冉 冊
diffusion distance D m.
␭m
0 ␾ 0共x兲
⌽m共x兲 ⫽ ␭1 ␾1共x兲 ,
m
The Heat Diffusion Map on Riemannian Submanifolds. Suppose that
· the data set X is approximately lying along a submanifold M 傺
·
· ⺢ n, with a density p(x) (not necessarily uniform on M). This kind
and the family of diffusion distances {D m} defined by of situation arises in many applications ranging from hyperspec-

Coifman et al. PNAS 兩 May 24, 2005 兩 vol. 102 兩 no. 21 兩 7427
 ⌬p(x)兾p(x) is a potential term capturing the influence of the
nonuniform density. By writing the nonuniform density in a Bolt-
zmann form, p(x) ⫽ e⫺U(x), the infinitesimal operator can be
expressed as
⌬␾ ⫹ 共储ⵜU储2 ⫺ ⌬U兲␾. [5]

This generator corresponds to the forward diffusion operator

and is the adjoint of the infinitesimal generator of the backward
operator, given by
⌬␾ ⫺ 2ⵜ␾䡠ⵜU. [6]

As is well known from quantum physics, for a double well

potential U, corresponding to two separated clusters, the first
nontrivial eigenfunction of this operator discriminates between
the two wells. This result reinforces the use of the standard graph
Laplacian for computing an approximation to the normalized cut
problem, as described in ref. 13 and more generally for the use
of the first few eigenvectors for spectral clustering, as suggested
by Weiss (14).
To capture the geometry of a given manifold, regardless of the
density, we propose a different normalization that asymptotically
recovers the eigenfunctions of the Laplace–Beltrami (heat) oper-
ator on the manifold. For any rotation-invariant kernel k␧(x, y) ⫽
h(储x ⫺ y储2兾␧), we consider the normalization described in the box
below. The operator A␧ can be used to define a discrete approxi-
mate Laplace operator as
I ⫺ A␧
⌬␧ ⫽ ,
␧

and it can be verified that ⌬ ␧ ⫽ ⌬ 0 ⫹ ␧ 1/2R ␧, where ⌬0 is a

multiple of the Laplace–Beltrami operator ⌬ on M, and R ␧ is
bounded on a fixed space of bandlimited functions. From this, we
can deduce the following result.

Theorem 2.1. Let t ⬎ 0 be a fixed number, then as ␧ 3 0, A ␧ ⫽

t/␧

(I ⫺ ␧⌬ ␧) t/␧ ⫽ (I ⫺ ␧⌬ 0) t/␧ ⫹ O(␧ 1/2) ⫽ e ⫺t⌬0 ⫹ O(␧ 1/2), and the

kernel of A t/␧
␧ is given as

a共␧t/␧兲共x, y兲 ⫽ 冘
jⱖ0
␭j共2t/␧兲␾j共␧兲共x兲␾j共␧兲共y兲

Fig. 1. The spectra of powers of A (a), and the diffusion embedding of a

⫽ 冘
jⱖ0
2
e⫺␯j t␾j共x兲␾j共y兲 ⫹ O共␧1/2兲
mixture of two materials with different heat conductivity (b and c). The
original geometry (b) is mapped as a ‘‘butterfly’’ set, in which the red (higher
conductivity) and blue phases are organized according to the diffusion they
⫽ ht共x, y兲 ⫹ O共␧1/2兲,
generate: the cord length between two points in the diffusion space measures
the quantity of heat that can travel between these points.
where { ␯ j2} and { ␾ j} are the eigenvalues and eigenfunctions of the
limiting Laplace operator, h t(x, y) is the heat diffusion kernel at
time t, and all estimates are relative to any fixed space of band-
tral imagery to image processing to vision. For instance, in the limited functions.
latter field, a model for edges can be generated by considering
pixel neighborhoods whose variability is governed by a few Approximation of the Laplace–Beltrami Diffusion Kernel.
parameters (11, 12). 1. Let p ␧(x) ⫽ 兰 X k ␧(x, y)p(y)dy, and form the new kernel k̂ ␧(x̆,
We consider isotropic kernels, i.e., kernels of the form y) ⫽ k ␧(x, y)兾p ␧(x)p ␧(y).

k␧共x, y兲 ⫽ h 冉 ␧
冊
储x ⫺ y储2
.
2. Apply the weighted graph Laplacian normalization to this kernel
by defining v␧(x) ⫽ 兰X k̂␧(x, y)p(y)dy, and by setting a␧(x, y) ⫽ k̂␧(x,
y)兾v␧(x).
2
In ref. 5, Belkin et al. suggest to take k ␧(x, y) ⫽ e ⫺(储x⫺y储 /␧) and Then the operator A ␧f(x) ⫽ 兰 X a ␧(x, y)f(y)p(y)dy is an approx-
to apply the weighted graph Laplacian normalization procedure imation of the Laplace–Beltrami diffusion kernel at time ␧.
described in the previous section. They show that if the density For simplicity, we assume that on the compact manifold M, the
of points is uniform, then as ␧ 3 0, one is able to approximate data points are relatively densely sampled (each ball of radius 公␧
the Laplace–Beltrami operator ⌬ on M. contains enough sample points so that integrals can approximated
However, when the density p is not uniform, as is often the case, by discrete sums). Moreover, if the data only covers a subdomain
the limit operator is conjugate to an elliptic Schrödinger-type of M with nonempty boundary, then ⌬0 needs to be interpreted as
operator having the more general form ⌬ ⫹ Q, where Q(x) ⫽ acting with Neumann boundary conditions. As in the previous

7428 兩 [Link]兾cgi兾doi兾10.1073兾pnas.0500334102 Coifman et al.

 ẋ ⫽ ⫺ⵜU共x兲 ⫹ 冑2ẇ, [7]

where U is the free energy and w(t) is the standard n-

dimensional Brownian motion. Let p(y, t 兩 x, s) denote the
transition probability of finding the system at location y at time
t, given an initial location x at time s. Then, in terms of the
variables {y, t}, p satisfies the forward Fokker–Planck equation
(FPE), for t ⬎ s,

Fig. 2. A dumbbell (a) is embedded by using the first three eigenfunctions ⭸p

⫽ ⵜ䡠共ⵜp ⫹ pⵜU共y兲兲, [8]
(b). Because of the bottleneck, the two lobes are pushed away from each ⭸t
other. Observe also that in the embedding space, point A is closer to the
handle (point B) than any point on the edge (like point C), because there are whereas in terms of the variables {x, s}, the transition probability
many more short paths joining A and B than A and C.
satisfies the backward equation

section, one can compute heat diffusion distances and the corre- ⭸p
⫺ ⫽ ⌬p ⫺ ⵜp䡠ⵜU共x兲. [9]
sponding embedding. Moreover, any closed rectifiable curve can be ⭸s
embedded as a circle on which the density of points is preserved: We
have thus separated the geometry of the set from the distribution As time t 3 ⬁, the solution of the forward FPE converges to the
of the points (see Fig. 3 for an example). steady-state Boltzmann density

Anisotropic Diffusion and Stochastic Differential Equations. So far we e⫺U共x兲

p共x兲 ⫽ , [10]

MATHEMATICS
have considered the analysis of general data sets by diffusion maps, Z

APPLIED
without considering the source of the data. One important case of
interest is when the data x is sampled from a stochastic dynamical where the partition function Z is the appropriate normalization
system. Consider, therefore, data sampled from a system x(t) 僆 ⺢n constant.
whose time evolution is described by the following Langevin The general solution to the FPE can be written in terms of an
equation eigenfunction expansion

Fig. 3. Original spiral curve (a) and the density of points on it (b), embedding obtained from the normalized graph Laplacian (c), and embedding from the
Laplace–Beltrami approximation (d).

Coifman et al. PNAS 兩 May 24, 2005 兩 vol. 102 兩 no. 21 兩 7429
 Fig. 4. The original function f on
the unit square (a), and the first
nontrivial eigenfunction (b). On
this plot, the colors corresponds to
the values of f.

冘
⬁ eigenvectors corresponding to the anisotropic kernel (Eq. 12)
p共x, t兲 ⫽ a je ⫺ ␭ jt
␾j共x兲, [11] capture the long-time asymptotic behavior of the stochastic system
(Eq. 7). Therefore, the diffusion map can be seen as an empirical
j⫽0
method for homogenization (see ref. 10 for more details). These
where ␭j are the eigenvalues of the Fokker–Planck operator, with variables are the right observables with which to implement the
␭0 ⫽ 1 ⬎ ␭1 ⱖ ␭2 ⱖ . . . , and with ␾j(x) the corresponding equation-free complex兾multiscale computations of Kevrekidis et al.
eigenfunctions. The coefficients aj depend on the initial conditions. (see refs. 15 and 16).
A similar expansion exists for the backward equation, with the
eigenfunctions of the backward operator given by ␺j(x) ⫽ eU(x)␾j(x). One-Parameter Family of Diffusion Maps. In the previous sections, we
As can be seen from Eq. 11, the long time asymptotics of the showed three different constructions of Markov chains on a discrete
solution is governed only by the first few eigenfunctions of the data set that asymptotically recover either the Laplace–Beltrami
Fokker–Planck operator. Whereas in low dimensions, e.g., n ⱕ 3, operator on the manifold, the backward Fokker–Planck operator
approximations to these eigenfunctions can be computed via nu- with potential 2U(x) for the normalized graph Laplacian, or U(x) for
merical solutions of the partial differential equation, in general, this the anisotropic diffusion kernel.
is infeasible in high dimensions. On the other hand, simulations of In fact, these three normalizations can be seen as specific cases
trajectories according to the Langevin Eq. 7 are easily performed. of a one-parameter family of different diffusion maps, based on the
An interesting question, then, is whether it is possible to obtain kernel
approximations to these first few eigenfunctions from (large
enough) data sampled from these trajectories. k␧共x, y兲
k共␧␣兲共x, y兲 ⫽ [13]
In the previous section we saw that the infinitesimal generator of p␣␧ 共x兲p␣␧ 共y兲
the normalized graph Laplacian construction corresponds to a
Fokker–Planck operator with a potential 2U(x) (see Eq. 6). There- for some ␣ ⬎ 0.
fore, in general, there is no direct connection between the eigen- It can be shown (9) that the forward infinitesimal operator
values and eigenfunctions of the normalized graph Laplacian and generated by this diffusion is
those of the underlying Fokker–Planck operator 8. However, it is
possible to construct a different normalization that yields infinites- H f共␣兲␾ ⫽ ⌬␾ ⫺ 共e共1⫺␣兲U⌬e⫺共1⫺␣兲U兲 ␾ . [14]
imal generators corresponding to the potential U(x) without the
additional factor of two. One can easily see that the interesting cases are (i) ␣ ⫽ 0,
Consider the following anisotropic kernel. corresponding to the classical normalized graph Laplacian; (ii)
␣ ⫽ 1, yielding the Laplace–Beltrami operator; and (iii) ␣ ⫽ 1兾2
k␧共x, y兲 yielding the backward Fokker–Planck operator.
k̃␧共x, y兲 ⫽
冑p␧共x兲p␧共y兲 [12] Therefore, while the graph Laplacian based on a kernel with
␣ ⫽ 1 captures the geometry of the data, with the density e ⫺U
A similar analysis to that of the previous section shows that the playing absolutely no role, the other normalizations take into
normalized graph Laplacian construction that corresponds to account also the density of the points on the manifold.
this kernel gives in the asymptotic limit the correct Fokker–
Planck operator, e.g., with the potential U(x). Directed Diffusion and Learning by Diffusion
Since the Euclidean distance in the diffusion map space corre- It follows from the previous section that the embedding that one
sponds to diffusion distance in the feature space, the first few obtains depends heavily on the choice of a diffusion kernel. In some

Fig. 5. Pathology slice with partially labeled data (a),

tissue classification from spectra by using 1-nearest neigh-
bors (b), and tissue classification from spectra by using
geometric diffusion (c). The three tissue classes are marked
with blue, green, and pink.

7430 兩 [Link]兾cgi兾doi兾10.1073兾pnas.0500334102 Coifman et al.

cases, one is interested in constructing diffusion kernels that are different tissue classes (marked with blue, green, and pink in Fig.
data- or task-driven. As an example, consider an empirical function 5a) and are asked to classify all pixels in the image using only
F(x) on the data. We would like to find a coordinate system in which spectral, as opposed to spatial, information. Both Gaussian classi-
the first coordinate has the same level lines as the empirical function fiers and nearest-neighbor classifiers (see Fig. 5b) perform poorly
F. For that purpose, we replace the Euclidean distance in the in this case as there is a gradual change in both shading and chemical
Gaussian kernel by the anisotropic distance. composition in the vertical direction of the tissue sample.
The diffusion framework, however, provides an alternative clas-
D2␧共x, y兲 ⫽ d2共x, y兲兾␧ ⫹ 兩F共x兲 ⫺ F共y兲兩2兾␧2 sification scheme that links points together by a Markov random
walk (see ref. 17 for a discussion): let ␹i be the L1-normalized
The corresponding limit of At/␧ ␧ is a diffusion along the level surfaces characteristic function of the initially labelled set from class i. At a
of F from which it follows that the first nonconstant eigenfunction given time t, we can interpret the diffused label functions (At␹i)i as
of A␧ has to be constant on level surfaces. This is illustrated in Fig. the posterior probabilities of the points belonging to class i. Choose
4, where the graph represents the function F and the colors a time ␶ when the margin between the classes is maximized, and
correspond to the values of the first nontrivial eigenfunction. In then define the label of a point x 僆 X as the maximum a posteriori
particular, observe that the level lines of this eigenfunction are the estimate L(x; ␶) ⫽ argmaxi A␶␹i. Fig. 5c shows the classification of
integral curves of the field orthogonal to the gradient of F. This is the pathology sample using the above scheme. The latter result
clear since we forced the diffusion to follow this field at a much agrees significantly better with a specialist’s view of correct tissue
faster rate, in effect integrating that field. It also follows that any classification.
differential equation can be integrated numerically by a nonisotro- In many practical situations, the user may want to refine the
pic diffusion in which the direction of propagation is faster along the classification of points that occur near the boundaries between
field specified by the equation. classes in state space. One option is to use an iterative scheme,
We now apply this approach to the construction of empirical where the user provides new labelled data where needed and then
models for statistical learning. Assume that a data set has been restarts the diffusion with the new enlarged training set. However,
generated by a process whose local statistical properties vary from if the total data set X is very large, an alternative, more efficient,

MATHEMATICS
location to location. Around each point x, we view all neighboring scheme is to define a modified kernel that incorporates both

APPLIED
data points as having been generated by a local diffusion whose previous classification results and new information provided by the
probability density is estimated by px(y) ⫽ cxexp(⫺qx(x ⫺ y)), where user: for example, assign to each point a score si(x) 僆 [0, 1] that
qx is a quadratic form obtained empirically by PCA from the data reflects the probability that a point x belongs to class i. Then use
in a small neighborhood of x. We then use the kernel a(x, z) ⫽ these scores to warp the diffusion so that we have a set of
兰 px(y)pz(y)dy to model the diffusion. Note that the distance defined class-specific diffusion kernels {Ãi}i that slow down diffusion be-
by this kernel is (兰兩px(y) ⫺ pz(y)兩2dy)1/2, which can be viewed as the tween points with different label probabilities. Choose, for example,
natural distance on the ‘‘statistical tangent space’’ at every point in in each new iteration, weights according to k̃i(x, y) ⫽ k(x, y)si(x)si(y),
the data. If labels are available, the information about the labels can where si ⫽ A␶␹i are the label posteriors from the previous diffusion,
be incorporated by, for example, locally warping the metric so that and renormalize the kernel to be a Markov matrix. If the user
the diffusion starting in one class stays in the class without leaking provides a series of consistent labelled examples, the classification
to other classes. This could be obtained by using local discriminant will speed up in each new iteration and the diffusion will eventually
analysis (e.g., linear, quadratic, or Fisher discriminant analysis) to occur only within disjoint sets of samples with the same labels.
build a local metric whose fast directions are parallel to the Summary
boundary between classes and whose slow directions are transversal
to the classes (see, e.g., ref. 1). In this article, we presented a general framework for structural
In data classification, geometric diffusion provides a powerful multiscale geometric organization of graphs and subsets of ⺢ n.
We introduced a family of diffusion maps that allow the explo-
tool to identify arbitrarily shaped clusters with partially labelled
ration of both the geometry, the statistics and functions of the
data. Suppose, for example, we are given a data set X with N points
data. Diffusion maps provide a natural low-dimensional embed-
from C different classes. Assume our task is to learn a function L :
ding of high-dimensional data that is suited for subsequent tasks
X 3 {1, . . . , C} for every point in X but we are given the labels of
such as visualization, clustering, and regression. In the compan-
only s ⬍⬍ N points in X. If we cannot infer the geometry of the data
ion article (18), we introduce multiscale methods that allow fast
from the label points only, many parametric methods (e.g., Gauss- computation of functions of diffusion operators on the data and
ian classifiers) and nonparametric techniques (e.g., nearest neigh- also present a scheme for extending empirical functions.
bors) lead to poor results. In Fig. 5, we illustrate this with an
example. Here, we have a hyperspectral image of pathology tissue. We thank Ioannis Kevrekidis for very helpful suggestions and discussions
Each pixel (x, y) in the image is associated with a vector {I(x, y)}␭ and Naoki Saito for useful comments during the preparation of the
that reflects the material’s spectral characteristics at different manuscript. This work was partially supported by the Defense Advanced
wavelengths ␭. We are given a partially labelled set for three Research Planning Agency and the Air Force Office of Scientific Research.

1. Hastie, T., Tibshirani, R. & Friedman, J. H. (2001) The Elements of Statistical Learning 11. Pedersen, K. S. & Lee, A. B. (2002) Proceedings of the 7th European Conference on Computer
(Springer, Berlin), pp. 144–155. Vision, Copenhagen, Denmark, eds. Nielsen, M., Heyden, A., Sparr, G. & Johansen, P.
2. Coifman, R. R. & Saito, N. (1994) C. R. Acad. Sci. 319, 191–196. (Springer, Berlin), pp. 328–342.
3. Ham, J., Lee, D. D., Mika, S. & Schölkopf, B. (2003) A Kernel View of the Dimensionality 12. Huggins, P. S. & Zucker, S. W. (2002) Proceedings of the 7th European Conference on
Reduction of Manifolds (Max-Planck-Institut für Biologische Kybernetik, Tübingen, Ger- Computer Vision, Copenhagen, Denmark, eds. Nielsen, M., Heyden, A., Sparr, G. &
many), Tech. Rep. TR-110, pp. 1–9. Johansen, P. (Springer, Berlin), pp. 384–398.
4. Roweis, S. T. & Saul, L. K. (2000) Science 290, 2323–2326. 13. Shi, J. & Malik, J. (2000) IEEE Trans. Pattern Anal. Machine Intell. 22, 888–905.
5. Belkin, M. & Niyogi, P. (2003) Neural Comput. 15, 1373–1396. 14. Weiss, Y. (1999) Segmentation Using Eigenvectors: a Unifying View, Proceedings of the
6. Donoho, D. L. & Grimes, C. (2003) Proc. Natl. Acad. Sci. USA 100, 5591–5596.
Institute of Electrical and Electronics Engineers International Conference on Computer Vision,
7. Zhang, Z. & Zha, H. (2002) Principal Manifolds and Nonlinear Dimension Reduction via
pp. 975–982.
Local Tangent Space Alignment, Tech. Rep. CSE-02-019 (Dept. of Computer Science and
Engineering, Pennsylvania State University), pp. 1–22. 15. Gar, C. W., Kevrekidis, I. G. & Theodoropoulos, C. (2002) Comput. Chem. Eng. 26,
8. Chung, F. R. K. (1997) Spectral Graph Theory (Conference Board of the Mathematical 941–963.
Sciences, Am. Math. Soc., Providence, RI). 16. Kevrekidis, I. G., Gear, C. W., Hyman, J. M., Kevrekidis, P. G., Runborg, O. &
9. Coifman, R. R. & Lafon, S. (2004) Diffusion Maps: Applied and Computational Harmonic Theodoropoulos, K. (2003) Commun. Math. Sci. 1, 715–762.
Analysis (Elsevier, New York), in press. 17. Szummer, M. & Jaakkola, T. (2001) Adv. Neural Inf. Process. Syst. 14, 945–952.
10. Nadler, B., Lafon, S., Coifman, R. R. & Kevrekidis, I. (2004) Applied and Computational 18. Coifman, R. R., Lafon, S., Lee, A. B., Maggioni, M., Nadler, B., Warner, F. & Zucker,
Harmonic Analysis (Elsevier, New York), in press. S. W. (2005) Proc. Natl. Acad. Sci. USA 102, 7432–7437.

Coifman et al. PNAS 兩 May 24, 2005 兩 vol. 102 兩 no. 21 兩 7431