\section{Introduction}

The recent discovery of topological insulators has generated great interest in the fields of condensed matter physics and materials science, largely driven by their exotic surface electronic properties.~\cite{Qi2010,Qi2011,Hasan2010}  These materials are insulating in the bulk, but support topologically protected helical edge or surface state in two- or three-dimensional systems.~\cite{Kane2005a,Kane2005b,Qi2010,Qi2011,Hasan2010}  The surface states are predicted to have special properties that could be useful for practical applications in spintronics and quantum computation.~\cite{Moore2010}  To date, topological states have been observed in a number of materials,~\cite{Teo2008,Hsieh2008,Bernevig2006,Konig2007,Brune2011,Fu2007a,Zhang2009,Xia2009,Chen2009,Xiao2010a,Chadov2010,Lin2010a,Feng2011,Feng2010,Al-Sawai2010,Liu2011a,Gofryk2011,Lin2010b,Yan2010a,Chen2010,Sato2010,Kuroda2010,Sun2010,Yan2010b,Kim2010,Jin2011,Zhang2011,Chen2011,Liu2011} opening the door for detailed and systematic investigations of various topological phenomena in laboratories.  While the essential features of topological insulators have been established, our understanding of the fundamental physics in relation to basic material properties is far from complete.  This type of knowledge will be crucial in potential device applications.

In this work, we theoretically demonstrate the possibility of tuning the topological order in cubic semiconductors by applying external strain, using InSb as a prototype system.  Our choice of the material system is motivated by the fact that many technologically important semiconductors crystalize in one of the cubic crystal structures, including the group-IV elements in the diamond structure, and the III-V and II-VI compounds in the zinc-blende structure.  Ternary compound semiconductors with the tetrahedral coordination, such as half-Heuslers~\cite{Xiao2010a,Chadov2010,Lin2010a} and chalcopyrites~\cite{Feng2011}, also belong to this category.  We first establish, based on tight-binding analysis, that lattice \emph{expansion} can induce a topological phase transition of small band-gap cubic semiconductors via a band inversion, and further breaking of the cubic symmetry leads to a topological insulating phase. Using density functional theory calculations, we then demonstrate that InSb enters into a topologically nontrivial state at a reasonable biaxial lattice expansion around $2\%\--3\%$.  The relationship between mechanical strain and topological order has been discussed in layered bismuth compounds~\cite{Young2011} as well as half-Heuslers~\cite{Xiao2010a,Chadov2010,Lin2010a}.  Our analysis provide a physical picture of these findings.

\section{General Considerations}

The band topology of materials with time-reversal symmetry is characterized by the so-called $Z_2$ topological index.  In general, the calculation of the $Z_2$ index requires the knowledge of the occupied Bloch wave functions in the entire Brillouin zone because it is a global property of the energy bands.  However, in most practical situations the $Z_2$ index can be determined by counting the number of band inversions at the time-reversal invariant $\bm k$-points.  In passing from a topologically trivial to nontrivial state, a band inversion, i.e., switching of occupied and unoccupied bands with opposite parity around the Fermi level, must happen.  If an odd number of band inversions occur within the Brillouin zone, then the material could be a topological insulator.

\begin{figure}[b]
\includegraphics[width=\columnwidth]{Fig1}
\caption{\label{fig:order}(Color online) The band order of cubic diamond/zinc-blende semiconductors at the $\Gamma$ point with $\lambda_\text{SO} > 0$.  (a) The normal band order.  (b)-(c) The inverted band order can be obtained from the normal band order by either increasing the spin-orbit coupling $\lambda_\text{SO}$ (b) or by decreasing the coupling potentials $V_{ss}$ and $V_{pp}$ (c).  The band inversion strength is defined as $\Delta E=E_{\Gamma_6}-E_{\Gamma_8}$, which has positive value in normal band order and negative value in inverted band order.}
\end{figure}

We now specialize to zinc-blende semiconductors.  In a typical band structure, the low-energy electronic properties is dominated by bands around the $\Gamma$ point.  Away from the $\Gamma$ point the valence and conduction bands are well separated.  Therefore we only need consider the band order at $\Gamma$.  The relevant states are the $\Gamma_6$ anti-bonding state of the $s$-orbitals, and the $\Gamma_7$ and $\Gamma_8$ bonding states of the $p$-orbitals.  Using the tight-binding model,~\cite{Slater1954,Chadi1977} the band energies of these states are obtained
\begin{align}
E_{\Gamma_6} &= \frac{\varepsilon_{s1} + \varepsilon_{s2}}{2}
+ \sqrt{\Bigl(\frac{\varepsilon_{s1} - \varepsilon_{s2}}{2}\Bigr)^2 + V_{ss}^2} \;, \\
E_{\Gamma_7} &= \frac{\varepsilon_{p1} + \varepsilon_{p2}}{2}
- \sqrt{\Bigl(\frac{\varepsilon_{p1} - \varepsilon_{p2}}{2}\Bigr)^2 + V_{pp}^2}
- \lambda_\text{SO} \;, \\
E_{\Gamma_8} &= \frac{\varepsilon_{p1} + \varepsilon_{p2}}{2}
- \sqrt{\Bigl(\frac{\varepsilon_{p1} - \varepsilon_{p2}}{2}\Bigr)^2 + V_{pp}^2}
+ \frac{1}{2}\lambda_\text{SO} \;,
\end{align}
where $\varepsilon_{s1,2}$ and $\varepsilon_{p1,2}$ are the $s$- and $p$-orbital energies of the cation and anion, respectively, $V_{ss}$ and $V_{pp}$ are the coupling potentials of the $s$-$s$ and $p$-$p$ bonds, and $\lambda_\text{SO}$ is the spin-orbit coupling strength.  Note that $\Gamma_6$ and $\Gamma_{7,8}$ states have opposite parity.

Usually, the $\Gamma_6$ state is located above the $\Gamma_{7,8}$ states, forming the conduction band minimum.  This is the so-called normal band order [see Fig.~\ref{fig:order}(a)].    There are two ways to change the band order.  (i) If the spin-orbit coupling strength is increased, depending on the sign of $\lambda_\text{SO}$, either the $\Gamma_8$ ($\lambda_\text{SO} > 0$) or the $\Gamma_7$ ($\lambda_\text{SO} < 0$) state could rise above the $\Gamma_6$ state, realizing the necessary band inversion [see Fig.~\ref{fig:order}(b)].  (ii) In addition to varying $\lambda_\text{SO}$, which is an intrinsic property of the material, the band order can be also changed by varying the coupling potentials $V_{ss}$ and $V_{pp}$ [see Fig.~\ref{fig:order}(c)].  Increasing the lattice constant leads to a decrease of the coupling potentials, which lowers the $\Gamma_6$ anti-bonding state and raises the $\Gamma_{7,8}$ bonding states.  At sufficiently large lattice constant, the band inversion can be realized.  Thus in the design of topological insulators, both factors should be taken into account.

The relation between the band gap and lattice constant has been discussed previously in the context of band gap pressure coefficient.~\cite{Wei1999,Li2006}  In a more realistic situation, there are other factors that could affect the band positions, such as the $p$-$d$ coupling in the compounds with active $d$-orbitals in the valence bands, as well as the kinetic energy effect under pressure.  A detailed analysis can be found in Refs. \onlinecite{Wei1999,Li2006}; the general trend we discussed above still holds.

Finally, if $E_{\Gamma_8} > E_{\Gamma_6}$ in the inverted band structure ($\lambda_\text{SO} > 0$), the twofold degenerate ${\Gamma}_6$ state is fully occupied while the fourfold degenerate ${\Gamma}_8$ state is half filled, resulting in a zero band-gap semiconductor [for example, see Fig.~\ref{fig:hydrostatic}(c) and \ref{fig:hydrostatic}(d)]. To realize a truly insulating ground state, a possible way is to apply a biaxial strain so as to break the cubic symmetry and thereby open a band gap at $\Gamma$ point  [see Fig.~\ref{fig:uniaxial}(d)].

\begin{figure}
\includegraphics[width=\columnwidth]{Fig2}
\caption{(Color online) The band structures of InSb at the equilibrium lattice constant (a) and with 1\% (b), 2\% (c), and 3\% (d) hydrostatic lattice expansions. The $\Gamma_{6}$, $\Gamma_{7}$, and $\Gamma_{8}$ state are denoted by red, black, and blue color, respectively. The size of dots is proportional to the probability of $s$-orbit projection.}
\label{fig:hydrostatic}
\end{figure}

\begin{table}
\caption{The relative energy changes of states $E_{\Gamma_6}$, $E_{\Gamma_7}$, and $E_{\Gamma_8}$ under the hydrostatic lattice expansions ($a=b=c$) ranging from $0\%$ to $4\%$.}\label{tab:energy1}
\begin{ruledtabular}
\begin{tabular}{cccc}
\multicolumn{1}{c}{\textrm{a(\%)}} &
\multicolumn{1}{c}{\textrm{V(Bohr$^3$)}} &
\multicolumn{1}{c}{\textrm{$E_{\Gamma_6}-E_{\Gamma_8}$(eV)}} &
\multicolumn{1}{c}{\textrm{$E_{\Gamma_8}-E_{\Gamma_7}$(eV)}} \\
\hline

 0  &  458.92429  &  0.25321  &  0.66261  \\
 1  &  472.83016  &  0.08445  &  0.65357  \\
 2  &  487.01414  & -0.08042  &  0.64499  \\
 3  &  501.47899  & -0.24181  &  0.63672  \\
 4  &  516.22745  & -0.39837  &  0.62886  \\

\end{tabular}
\end{ruledtabular}
\end{table}

\section{Methodology}

Having established the general tuning principle, next we verify whether the anticipated topological phase transition can happen in real materials at a reasonable lattice expansion. We choose InSb as a prototype system because it has a relatively small band gap of 0.235 eV,\cite{Madelung2004} indicating that InSb is on the verge of becoming a topological insulator.

The band structure calculations in this work were performed using full-potential linearized augmented plane-wave method\cite{Singh1994,Blugel2006} implemented in package \textsc{wien2k}.\cite{Blaha2001}  A converged ground state was obtained using $10^4$ $\bm k$-points in the first Brillouin zone with $R_\text{MT} K_\text{max} = 9.0$, where $R_\text{MT}$ represents the smallest muffin-tin radius and $K_\text{max}$ is the maximum size of reciprocal-lattice vectors. The muffin-tin radius was set to $2.5$ Bohr for both  In and Sb atoms. Wave functions and potentials inside the atomic sphere are expanded in spherical harmonics up to $l=10$ and $4$, respectively.  Spin-orbit coupling is included by a second-variational procedure,~\cite{Singh1994} where states up to 9 Ry above Fermi level are included in the basis expansion.  The modified Becke-Johnson exchange potential together with local-density approximation for the correlation potential (MBJLDA)~\cite{Tran2009} was used to obtain the band structures.  The MBJDA potential can effectively mimic the behavior of orbital-dependent potential around the band gap, and it is expected to obtain accurate positions of states near the band edge~\cite{Tran2009}.  In fact,  both standard  local-density approximation (LDA)\cite{Kohn1965,Perdew1992} and generalized gradient approximation (GGA)\cite{Perdew1996} calculations predict that InSb is a zero-band gap semiconductor with an inverted band order, while the MBJLDA calculation produce a band gap about 0.253 eV, which is in excellent agreement with the experimental value of 0.235 eV.~\cite{Madelung2004}  The MBJLAD method has been applied for small band-gap half-Heuslers in the context of topological insulators~\cite{Feng2010,Al-Sawai2010} and yields favorable comparison with more sophisticated GW method.~\cite{Vidal2011}

We consider two types of lattice expansion. The first is the hydrostatic lattice expansion,  equally increasing the lattice constants along all the three axes. This situation is more of academic interest rather than for practical applications. In the second case, we consider the biaxial lattice expansion, where the crystal structure is expanded in the $ab$-plane and relaxed in the $c$-axis. This can be achieved by growing the sample on a substrate with a larger lattice constant. Although the MBJLDA potential gives very accurate band gap and band order, it has no corresponding definition about its energy functional.~\cite{Tran2009}  Hence, both the LDA and GGA for the exchange-correlation potential were used to optimize the lattice constants.

\begin{figure}
\includegraphics[width=\columnwidth]{Fig3}
\caption{(Color online) The evolution of band order of InSb with nonhydrostatic expansion (a) in $ab$-plane and (b) along $c$-axis. The black solid circles are the results without lattice optimization, orange down-diamonds and purple up-diamonds are optimized by LDA and GGA potential receptively. (c) and (d) are normal band order and inverted band order, which correspond to the states marked by dotted circles in (a). The size of dots is proportional to the probability of $s$-orbit projection.}
\label{fig:uniaxial}
\end{figure}

\section{Result and Discussion}

The crystal structure of InSb is zinc-blende with space group $F\bar{4}3m$ (No. 216).
The unit cell of InSb contains two inequivalent atoms, In and Sb, with the fractional coordinates (0,0,0) and (0.25,0.25,0.25), respectively. The initial equilibrium lattice constant of 6.479 $\text{\AA}$ is taken from previous  experimental measurements.\cite{Madelung2004}

Figure~\ref{fig:hydrostatic} shows the band structures of InSb at its equilibrium lattice constant and under hydrostatic lattice expansion ranging from $1\%$ to $3\%$.  The relative energy changes of the states $E_{\Gamma_6}$, $E_{\Gamma_7}$, and $E_{\Gamma_8}$ are listed in Table~\ref{tab:energy1}.  We find that, as anticipated, the band inversion takes place upon lattice expansion, somewhere between 1\% and 2\%.  In addition, we find that the spin-orbit coupling strength, characterized by $E_{\Gamma_8} - E_{\Gamma_7}$, stays almost a constant with a slight decrease as the lattice constant increases.  Hence the evolution of the band order during lattice expansion is dominated by the change of the coupling potentials $V_{ss}$ and $V_{pp}$.

Although a hydrostatic expansion already converts the normal band order into the inverted band order, a nonhydrostatic strain is still needed to creat a band gap at $\Gamma$ point by breaking the cubic symmetry.  As shown in Fig.~\ref{fig:uniaxial}, we investigate the evolution of band order with the nonhydrostatic expansion both in $ab$-plane and along the $c$-axis.  Note that in the tetragonal symmetry the band inversion strength $\Delta E$ is redefined as the energy difference between the $s$-like states and the top of valence bands.  First, we consider the lattice expansion in the $ab$-plane with fixed lattice constant along the $c$-axis (without lattice optimization).  The $\Delta E$ change from positive values to negative values as shown in Fig. \ref{fig:uniaxial}(a) (black solid circles). The band structures of two typical states are presented in Fig.~\ref{fig:uniaxial}(c) and \ref{fig:uniaxial}(d).  One is normal band order with $1\%$ expansion [Fig.~\ref{fig:uniaxial}(c)] and another is inverted band order with $3\%$ expansion [Fig.~\ref{fig:uniaxial}(d)].  Then, we consider the lattice optimization along the $c$-axis by using both LDA and GGA potential.  It is well known that LDA usually underestimates the lattice constant and GGA overestimates it, hence these results give us a lower and upper bound.  The true transition point for the expansion in $ab$-plane is in the range of $2\%\sim3\%$. Finally, if we expand the lattice along $c$-axis, there exist similar topological order transition. The only quantitative difference is that a larger lattice expansion is needed, as shown in \ref{fig:uniaxial}(b).

Finally, we mention the possibility that the tensile biaxial expansion may also be realized by growing the cubic semiconductor on the top of a piezoelectric substrate whose lattice can be changed by applying an electric field. In this way, the topological order of the semiconductor can be tuned with electric field.

\section{Summary}

In summary, we have explored the strain tuning of topological band order in cubic semiconductors by the combination of simple tight-binding analysis and density functional theory calculation. We have predicted that InSb can realize the topological insulating state under the biaxial lattice expansion of $2\%-3\%$. This work provides a generic guiding principle for tuning the topological order of the cubic semiconductors and offers an opportunity for experimentally exploring the properties of the topological surface states on such technologically relevant materials for practical applications.

\acknowledgements

This work is supported by the Laboratory Directed Research and Development Program of Oak Ridge National Laboratory, managed by UT-Battelle, LLC, for the U. S. Department of Energy.


\section{Background}

Let us be given some data represented as a $k$-dimensional manifold $X \subset \mathbb{R}^d$, embedded into a $d$-dimensional Euclidean space. In many applications $d$ is very large while the intrinsic dimension of the data $k$ is small, and one tries to study the structure of the manifold rather than its $d$-dimensional embedding.
Such a structure can be characterized by the means of the {\em Laplace-Beltrami operator}.
In the discrete setting, the manifold is often represented by a weighted graph with vertices $\{\bb{x}_{1},\hdots,\bb{x}_{n}\} \subset X$ and edge weights $w_{ij}=k(\bb{x}_{i},\bb{x}_{j})$ representing local connectivity using e.g. Gaussian kernel (see \cite{Luxburg}).
The Laplace-Beltrami operator can be discretized\footnote{There exist many different constructions of the discrete Laplacian. For the sake of simplicity, we adopt the symmetric Laplacian. Our framework is applicable to other discretization as well.} as $\bb{L}=\bb{D}^{-1/2}(\bb{D} - \bb{W})\bb{D}^{-1/2}$, where $\bb{W} = (w_{ij})$ and $\bb{D} = \mathrm{diag}(\sum_{j}w_{ij})$. Such a discretization is often referred to as {\em symmetric normalized Laplacian} and admits a unitary diagonalization $\bb{L}=\bb{V}\bb{\Lambda}\bb{V}^{\mathrm{T}}$, $\bb{V}\Vv^{\mathrm{T}}=\bb{I}_{n}$
with the eigenvalues $\lambda_{1}=0\leq\lambda_{2}\leq\hdots\leq\lambda_{n}$.
Geometric constructions associated with eigenvectors and eigenvalues of the Laplacian play an important role in machine learning, since several archetypical  problems can be formulated in these terms:

{\bf Eigenmaps. }
Non-linear dimensionality reduction methods try to capture the intrinsic low-dimensional structure of the manifold $X$.
 %
Belkin and Niyogi (\citeyear{Belkin02laplacianeigenmaps}) showed that finding
a neighborhood-preserving $k$-dimensional embedding of $X$ can be posed as the
minimum eigenvalue problem,
\begin{eqnarray}
\min_{\bb{V}\in\mathbb{R}^{n\times k}}\mathrm{tr}\,(\bb{V}^{\mathrm{T}}\bb{L} \bb{V}) \,\,\, \mathrm{s.t.} \,\,\, \bb{V}^{\mathrm{T}}\bb{V}=\bb{I}. \label{eq:belkin}
\end{eqnarray}
This problem is minimized by setting $\bb{V}$ to be the matrix containing the first $k$ eigenvectors of $\bb{L}$, thus effectively embedding the data by means of the eigenfunctions of the Laplace-Beltrami operator (the null eigenvector is usually discarded). Such an embedding is referred to as {\em Laplacian eigenmap}.
More generally, a {\em diffusion map} is given as a mapping of the form $\Psi = (K(\lambda_2)\bb{v}_2,\hdots, K(\lambda_k) \bb{v}_{k})$,
where $K(\lambda)$ is some transfer function acting as a ``low-pass filter'' on eigenvalues $\lambda$ (\cite{Coifman05geometricdiffusions, Coifman}).

{\bf Diffusion distances. }
Coifman et al. (\citeyear{Coifman05geometricdiffusions, Coifman}) related the eigenmaps to heat diffusion and random processes on manifolds and defined a family of {\em diffusion metrics} that in the most general setting can be written as
\begin{equation}
d^{2}(\bb{x}_{i},\bb{x}_{j})=\sum_{l} K(\lambda_l) (v_{il} - v_{jl})^2 = \| \Psi(\bb{x}_i) - \Psi(\bb{x}_j) \|^2_2.
\label{eq:diffdist}
\end{equation}
Particular choice of $K(\lambda) = e^{-\lambda t}$ gives the {\em heat diffusion distance}, related to the connectivity of points $\bb{x}_i, \bb{x}_j$ on the manifold by means of diffusion process of length $t$.
Such distances are intrinsic and thus invariant to manifold embedding and are robust to topological noise.

 {\bf Spectral clustering. }
 %
\cite{Ng01onspectral} showed a very efficient and robust clustering approach based on the observation that the multiplicity of the null  eigenvalue of $\bb{L}$ is equal to the number of connected components of $X$. The corresponding eigenvectors act as indicator functions of these components.
  %
  Embedding the data using these eigenvectors and then applying some standard clustering algorithm such as K-means was shown to produce significantly better results than clustering the high-dimensional data directly.

\section{Discussion and Conclusions}

We presented a framework for multi-modal data analysis using approximate joint diagonalization of Laplacian matrices, naturally extending the
classical construction of diffusion geometry to the multi-modal setting. This construction allowed an almost straightforward extension of
various diffusion-geometric data analysis tools such as spectral clustering and manifold learning based on diffusion maps.
In follow-up studies, we intend to show multi-modal extensions of other related techniques such as spectral hashing.

We also showed that many previously proposed approaches to multi-modal spectral clustering are nearly equivalent
and try to solve some version of the joint approximate diagonalization problem.
From the numerical perspective, existing methods were tailored for computing the null joint eigenvectors that are sought for in clustering problems.
The underlying optimization problems are poorly suited for broader applications of diffusion geometry such as non-linear dimensionality reduction and manifold learning, where many or all eigenvectors of the Laplacians are of interest.
While approximate joint  diagonalization methods developed in the signal processing community for source separation problems can address the latter case, they were initially developed for full matrices and do not take advantage of the sparse structure of Laplacians.

To the best of our knowledge, there currently exists no efficient tool to compute the joint eigenvectors of very large sparse matrices, akin Matlab's \texttt{eigs}.
We believe that the presented construction makes the need of such a tool central enough to deserve the interest of the entire machine learning community. In future work, we will consider extending standard methods for eigendecomposition of large sparse matrices to the joint diagonalization case.

\section{Introduction}

The Laplacian operator and related constructions play a pivotal role
in a wide range of applications in machine learning, pattern recognition,
and computer vision community. It has been shown that many problems
in these fields boil down to finding some eigenvectors and eigenvalues
of a Laplacian constructed on some high-dimensional data.
Important examples include {\em spectral clustering} (\citet{Ng01onspectral}) where clusters
are determined by the first eigenvectors of the Laplacian; {\em eigenmaps} (\citet{Belkin02laplacianeigenmaps})
and more generally {\em diffusion maps} (\cite{Coifman}), where one tries to find
a low-dimensional manifold structure using the first smallest eigenvectors
of the Laplacian; and {\em diffusion metrics} (\cite{Coifman05geometricdiffusions}) measuring the ``connectivity''
of points on a manifold and expressed through the eigenvalues and
eigenvectors of the Laplacian.
Other applications heavily relying
on the properties of the Laplacian include
{\em spectral
graph partitioning} (\cite{DingHeetal2001}), {\em spectral hashing} (\cite{weiss2008spectral}), spectral correspondence, image segmentation (\cite{Shi97normalizedcuts}),
and
shape analysis (\cite{Levy:2006:LET:1136647.1136965}). Because of
the intimate relation between the Laplacian operator, Riemannian geometry,
and diffusion processes, it is common to encounter the
umbrella term {\em spectral} or {\em diffusion geometry} in relation to the above
problems.

These applications have been considered mostly
in the context of uni-modal data, i.e., a single data space. However,
many applications involve observations and measurements of data done
using different modalities, such as multimedia documents  (\cite{weston2010large,rasiwasia2010new,mcfee2011learning}), audio and video (\cite{kidron2005pixels,alameda2011finding}), or medical imaging modalities like PET and CT (\cite{bronstein2010data}).
Such problems of multimodal (or multi-view)
data analysis have gained increasing interest in the computer vision
and pattern recognition communities, however, there have been only
few attempts extending the powerful spectral methods to such settings.

In this paper, we propose a general framework allowing to extend different
diffusion and spectral methods to the multimodal setting by finding a common eigenbasis of multiple Laplacians. Numerically, this problem is posed as {\em approximate
joint diagonalization} of several matrices. Such methods have received limited attention in the numerical mathematics community (\cite{Bunse-Gerstner:1993}) and have been employed for joint diagonalization of covariance matrices in blind source separation applications by \cite{CardosoBlind1993,Cardoso96jacobiangles,YeredorJD202,ZieglerBlindsource05}. To the best of our knowledge, this is the first time they are applied to spectral embeddings.
Besides providing a principled approach to data fusion, our framework gives a theoretical explanation to existing methods for multimodal data analysis. In particular, we show that many recent works on multi-view clustering by \cite{de2005spectral,MaAnchor2008,tang2009clustering,Cai2011,kumarco} can be considered a particular instance of our framework.

\section{Multimodal diffusion geometry}

Recently, we witness increasing popularity of attempts to analyze different ``views'' or modalities of data.
Such data can be modeled as
$m$ different manifolds $X^1 \subset \mathbb{R}^{d_1}, \hdots, X^m \subset \mathbb{R}^{d_m}$, which can have embeddings of different dimensionality ($d_1,\hdots, d_m$) and sometimes different structure.
We are interested in analyzing these manifolds simultaneously in order to extract their joint intrinsic structure.
We assume that we are given $n$ corresponding samples $\{(\bb{x}^i_1,\hdots,\bb{x}^i_n)\subset \mathbb{R}^{d_i}\}_{i=1}^m$ on the manifolds
and can construct the Laplacian matrices $\bb{L}_1,\hdots, \bb{L}_m$ as described in the previous section.

\begin{figure}[t!]
\centering
\includegraphics[width=1\linewidth]{figures/swiss_evecs_.pdf}
\caption{\label{fig:evecs}First and second rows: eigenfunctions of the Laplacians of two modalities of the Swiss roll.
Third and fourth rows: joint eigenfunctions of the two Laplacians computed using JADE. Hot colors represent positive values; cold colors represent negative values. }
\end{figure}

Trying to use the eigenvectors $\bb{V}_1,\hdots,\bb{V}_m$ of the Laplacian matrices $\bb{L}_1,\hdots, \bb{L}_m$ is problematic:
for a set of eigenvectors corresponding to an eigenvalue with multiplicity greater than one, we can talk only of eigen sub-space, and any basis spanning it is a valid set of eigenvectors.
As a result, the eigenvectors of the Laplacians in different modalities can be substantially different (Figure~\ref{fig:evecs}, top).

{\bf Joint diagonalization. } A solution is to try to find the eigenbasis of the Laplacians {\em simultaneously}.
This problem is known as {\em joint diagonalization} and consists of finding a set of joint orthogonal eigenvectors $\bar{\bb{V}}$
such that $\bar{\bb{V}}^\mathrm{T} \bb{L}_i \bar{\bb{V}} = \bb{\Lambda}_i$ are diagonal matrices of the eigenvalues of $\bb{L}_i$. Such a common eigenbasis solves the inherent ambiguity in the definition of the eigenvectors and ``couples'' different modalities (Figure~\ref{fig:evecs}, bottom).
However, due to differences between the modalities and the presence of noise, the Laplacian matrices $\bb{L}_1,\hdots, \bb{L}_m$ rarely have a joint eigenbasis (iff they commute).
It is still possible to find an {\em approximate} joint diagonalization by solving
\begin{equation}
\min_{\bar{\bb{V}}}\sum_{i=1}^{m} \mathrm{off}(\bar{\bb{V}}^{\mathrm{T}}\bb{L}_{i}\bar{\bb{V}}), \,\,\, \mathrm{s.t.} \,\,\, \bar{\bb{V}}^\mathrm{T} \bar{\bb{V}} = \bb{I}, \label{eq:costfunction}
\end{equation}
where $\mathrm{off}(\bb{X})$ is some off-diagonality criterion, e.g. the sum of squared off-diagonal elements,  $\mathrm{off}(\bb{X}) = \| \bb{X} - \mathrm{diag}(\bb{X})\|_\mathrm{F}^{2}$.
In this case, $\bar{\bb{V}}^\mathrm{T} \bb{L}_i \bar{\bb{V}}$ are only approximately diagonal; we refer to the average of the diagonal elements $\bar{\bb{\Lambda}} = \frac{1}{m}\sum_{i=1}^m \mathrm{diag}(\bar{\bb{V}}^\mathrm{T} \bb{L}_i \bar{\bb{V}})$ as the {\em joint approximate eigenvalues} of $\bb{L}_1,\hdots, \bb{L}_m$.
This definition allows us to naturally extend the diffusion geometric methods discussed in the previous section (eigenmaps, diffusion distances, spectral clustering, etc.)
to the multimodal setting by simply replacing the eigenvalues and eigenvectors of a single Laplacian $\bb{L}_i$ by the joint eigenvectors  $\bar{\bb{V}}$ and eigenvalues $\bar{\bb{\Lambda}}$ of multiple Laplacians $\bb{L}_1, \hdots, \bb{L}_m$.

{\bf Numerical computation. }
A numerical method for joint diagonalization based on a modified  \emph{Jacobi iteration} traces back to \cite{Bunse-Gerstner:1993}, and it has been used at about the same time by Cardoso and Souloumiac (\citeyear{CardosoBlind1993,Cardoso96jacobiangles}) for joint diagonalization of covariance matrices in the context of blind source separation.
The idea of the standard Jacobi method for eigenvalue calculation is to apply a sequence of plane rotations in order to sequentially minimize the off-diagonal elements of the given matrix. The rotation is applied ``in-place'' and does not require matrix multiplication.
In the modified Jacobi method (referred to as JADE), the rotations are applied to reduce the off-diagonality criterion (\ref{eq:costfunction}) in each step. Let $\bb{R}_{pqcs}$ the (complex) rotation matrix the entries of which are equal to those of the identity matrix except for the elements
\begin{equation}
\left(\begin{array}{cc}
r_{pp} & r_{pq}\\
r_{qp} & r_{qq}
\end{array}\right)=\left(\begin{array}{cc}
c & \bar{s}\\
-s & \bar{c}
\end{array}\right)
\end{equation}
where $\mid c\mid^{2}+\mid s\mid^{2}=1$. \citet{Cardoso96jacobiangles}
show that the problem
\begin{equation}
\min_{\mid c\mid^{2}+\mid s\mid^{2}=1}\sum_{i=1}^{m} \mathrm{off}(\bb{R}_{pqcs}^{\mathrm{T}}\bb{L}_{i}\bb{R}_{pqcs}) \,\,\,
\end{equation}
has a simple explicit solution based on a $3\times3$ eigenvalue problem.
JADE is one of the most common algorithms in the field of joint diagonalization and has complexity comparable to that of the standard Jacobi method.  There are other algorithms, like the  ACDC method of \citet{YeredorJD202}, as well as different versions of the idea of minimizing a suitable cost function on the Stiefel manifold (\citet{Rahbar00geometricoptimization}).

{\bf Analytic computation. }
In the spectral clustering problem, we are looking for the null eigenvectors of the Laplacian.
Assuming that the first $k$ eigenvalues of the Laplacians are zero, we want to find $\bar{\bb{V}} \in \mathbb{R}^{n\times k}$ such that $\bb{L}_i \bar{\bb{V}} = 0$ for all $i=1,\hdots, m$ and $\bar{\bb{V}}^\mathrm{T} \bar{\bb{V}} = \bb{I}$  by reformulating~(\ref{eq:costfunction}) as
\begin{equation}
\min_{\bar{\bb{V}} \in \mathbb{R}^{n\times k}}  \sum_{i=1}^m \| \bb{L}_i \bar{\bb{V}} \|^2_\mathrm{F}, \,\,\, \mathrm{s.t.} \,\,\, \bar{\bb{V}}^\mathrm{T} \bar{\bb{V}} = \bb{I}.
\label{eq:analytic1}
\end{equation}
Since $\sum_{i=1}^m\| \bb{L}_i \bar{\bb{V}} \|^2_\mathrm{F}= \mathrm{tr}\, (\bar{\bb{V}}^\mathrm{T}(\sum_{i=1}^m \bb{L}_i^\mathrm{T} \bb{L}_i )\bar{\bb{V}})$, the problem can be equivalently recast as single-modality clustering with the ``average''  Laplacian matrix
$\bar{\bb{L}} =  \sum_{i=1}^m \bb{L}_i^\mathrm{T} \bb{L}_i$.
We can also consider other averaging operators, e.g. weighted arithmetic mean
$\bar{\bb{L}} =  \sum_{i=1}^m w_i\bb{L}_i$ or harmonic mean
$\bar{\bb{L}} =  (\sum_{i=1}^m \bb{L}^{-1}_i)^{-1}$.
We discuss these methods in the next section.

For zero eigenvalues, (\ref{eq:analytic1}) is akin to (\ref{eq:costfunction}), which justifies the successful use of such ``averaging'' methods in problems of multimodal spectral clustering (\cite{MaAnchor2008,Cai2011}).
However, iterative methods such as JADE explicitly minimizing the off-diagonality criterion (\ref{eq:costfunction}) are more generic and applicable to settings where one has to find all or many joint eigenvectors, e.g., for computing eigenmaps or diffusion distances.

\section{Relation to previous works}

There have been numerous recent works on multimodal spectral-type clustering proposing different ways of fusing multiple modalities based on different principles.
Considering these methods through the prism of joint diagonalization, we show many commonalities and equivalences between algorithms stemming from different motivations and coming from various communities.
\cite{MaAnchor2008} considered detection of shots in video sequences using fusion of video and audio information, employing for this purpose spectral clustering of a Laplacian created as a weighted arithmetic mean of each modality Laplacian.
\cite{tang2009clustering}
used low-rank factorization of the weight matrix, trying to find a common factor $\bb{U}$ such that $\bb{W}_i \approx \bb{U} \bb{\Lambda}_i \bb{U}^\mathrm{T}$
by solving
\begin{equation}
\min_{\bb{U} \in \mathbb{R}^{n\times k}, \Lambda_i \in \mathbb{R}^{n\times n}} \sum_{i=1}^m \| \bb{W}_i - \bb{U} \bb{\Lambda}_i \bb{U}^\mathrm{T} \|_\mathrm{F}^2,
\label{eq:tang1}
\end{equation}
using the quasi-Newton method.
Besides the fact that the factorization is applied to the weight matrix (it can be equivalently applied to the Laplacian), we see here a (non-orthogonal) joint diagonalization problem with an off-diagonality criterion considered by \cite{YeredorJD202}.

\cite{Cai2011} proposed a method for {\em multiview spectral clustering} (MVSC) by solving\footnote{\cite{Cai2011} also impose a non-negativity constraint on the matrix $\bb{V}$ in order to obtain cluster indicators directly and bypass the K-means clustering stage. We ignore this additional constraint for the simplicity of discussion; such a constraint can be added to all the problems discussed in this paper. }
\begin{equation}
\min_{\bb{V}_i, \bb{V} \in \mathbb{R}^{n\times k} } \sum_{i=1}^{m}\mathrm{tr}\,(\bb{V}_{i}^\mathrm{T}\bb{L}_{i}\bb{V}_{i})+\alpha \|\bb{V}_{i}-\bb{V}\|_\mathrm{F}^{2}\,\,\, \mathrm{s.t.} \,\,\, \bb{V}^\mathrm{T} \bb{V} = \bb{I}
\label{eq:Caiminimization}
\end{equation}
The authors show that this problem can be equivalently posed as
\begin{equation}
\max_{\bb{V} \in \mathbb{R}^{n\times k} } \mathrm{tr}\, \left( \bb{V}^\mathrm{T}  \textstyle \sum_{i=1}^{m}\left(\bb{L}_{i} + \alpha \bb{I} \right)^{-1} \bb{V} \right)\,\,\, \mathrm{s.t.} \,\,\, \bb{V}^\mathrm{T} \bb{V} = \bb{I},
\label{eq:Caiminimization1}
\end{equation}
and then employ an iterative algorithm to find the solution $\bb{V}$.
First, we observe that problem (\ref{eq:Caiminimization}) consists of $m$ minimum-eigenvalue problems w.r.t. bases $\bb{V}_i$, with the addition of a coupling term,
encouraging $\bb{V}_i$ as close as possible to some common basis $\bb{V}$ (note that the authors do not impose orthogonality constraints $\bb{V}_i^\mathrm{T} \bb{V} = \bb{I}$, but for $\alpha \gg 0$, the proximity to orthogonal $\bb{V}$ makes $\bb{V}_i$ approximately orthogonal).
Thus, it is possible to interpret~(\ref{eq:Caiminimization})  as a kind of joint diagonalization criterion.
Second, problem (\ref{eq:Caiminimization1}) can be rewritten as a minimum eigenvalue problem
\begin{equation}
\min_{\bb{V} \in \mathbb{R}^{n\times k} } \mathrm{tr}\, \left( \bb{V}^\mathrm{T}  \left( \textstyle \sum_{i=1}^{m}\left(\bb{L}_{i} + \alpha \bb{I} \right)^{-1} \right)^{-1} \bb{V} \right)\,\,\, \mathrm{s.t.} \,\,\, \bb{V}^\mathrm{T} \bb{V} = \bb{I},
\label{eq:Caiminimization2}
\end{equation}
whose solution is given by the matrix composed of the first $k$ eigenvectors of the matrix
$\left( \textstyle \sum_{i=1}^{m}\left(\bb{L}_{i} + \alpha \bb{I} \right)^{-1} \right)^{-1}$. For $\alpha>0,$ this a regularized version of the \emph{harmonic mean} of the Laplacian matrices. We can thus regard the method of \cite{Cai2011} as a particular instance of our joint diagonalization approach discussed in the previous section.

\cite{kumarco} proposed
the {\em centroid co-regularization} approach for multimodal clustering based on the minimization of
\begin{equation}
\min_{\bb{V}, \bb{V}_i \in \mathbb{R}^{n\times k}  } \sum_{i=1}^m \mathrm{tr}\, ( \bb{V}_i^\mathrm{T}  \bb{L}_i \bb{V}_i ) - \alpha  \mathrm{tr}\, ( \bb{V}_i \bb{V}_i^\mathrm{T} \bb{V} \bb{V}^\mathrm{T} ) \,\,\, \mathrm{s.t.} \,\,\, \bb{V}_i^\mathrm{T} \bb{V}_i = \bb{I}; \,\,\, \bb{V}^\mathrm{T} \bb{V} = \bb{I}.
\label{eq:kumar2}
\end{equation}
This function is alternatingly minimized, first with respect to the $\bb{V}_i$, then with respect to $\bb{V}$.
Problems~(\ref{eq:kumar2}) and~(\ref{eq:Caiminimization}) are similar in their spirit (the first one uses dissimilarity $\| \bb{V}_i - \bb{V} \|_\mathrm{F}^2$ as coupling term, while the second one the similarity $\mathrm{tr}\, ( \bb{V}_i \bb{V}_i^\mathrm{T} \bb{V} \bb{V}^\mathrm{T} ) = \| \bb{V}_i^\mathrm{T} \bb{V} \|_\mathrm{F}^2$)), and  fall under our joint diagonalization framework.

We must stress that these methods were developed for clustering problems where one has to find the null eigenvectors, and do not adapt easily to other applications of diffusion geometry where one has to find many or all joint eigenvectors of the Laplacians (e.g., computation of diffusion distances).
In particular, iterative solvers used in \cite{tang2009clustering,kumarco,Cai2011} do not scale up to such cases.
On the other hand, algorithms such as modified Jacobi iteration (JADE) are made for finding a full set of joint eigenvectors and have the complexity akin to standard Jacobi iteration. Further speed-up might be achieved by making explicit use of the sparse structure of the Laplacian matrices, which is not taken advantage of in JADE.

\section{Results}

We tested the proposed approach on three applications: dimensionality reduction, diffusion distance, and spectral clustering.
All the datasets and code generating the results in this section are available from {\tt anonymous.com}.
Additional results are shown in the supplementary material.

{\bf Swiss rolls. }
In the first experiment, we used two Swiss roll surfaces with slightly different embedding as two different data modalities. The rolls were constructed in such a way that in each modality there is topological noise (connectivity ``across'' the roll loops) at different points. Laplacians were constructed as in \cite{Belkin02laplacianeigenmaps} using $5$-neighbor
connectivity and Gaussian weights with scale parameter $t$.
Figure~\ref{fig:evecs} shows the first few eigenvectors computed using each Laplacian individually and jointly.
Figure~\ref{fig:eigenmap} shows two-dimensional embeddings of the same surfaces using the first non-trivial eigenvectors.
When using joint eigenvectors, we are able to correctly capture the intrinsic structure of the data.
Figure~\ref{fig:dist} shows the diffusion distance on the Swiss roll surfaces, computed using the first 100 eigenvectors and heat diffusion kernel $K(\lambda) = e^{-1000\lambda}$. Topological noise is clearly visible especially in the first modality, resulting in the distance between two loops to be small. This phenomenon does not occur when using joint eigenvectors.

{\bf Synthetic data clustering. } In the second experiment, we performed clustering on several synthetic multimodal datasets.
Laplacians were constructed using 15 nearest neighbors (10 for the {\em circles}), and Gaussian weight selected using the self-tuning approach of \cite{perona2004self}.
We compare spectral clustering based on single modalities (SC-1 and SC-2) and joint diagonalization obtained using the JADE method of \cite{Cardoso96jacobiangles}; harmonic mean (JD-HM) of Laplacians (\cite{Cai2011}); and a non-spectral Comraf clustering algorithm (\cite{BekkermanJ07}). Quality was measured using the clustering accuracy criterion as defined in \cite{BekkermanJ07}. For {\em Blobs}, accuracy is averaged over 100 experiments ran on randomly generated datasets.

The results are summarized in Figure~\ref{fig:clustering1} and Table~\ref{tab:clustering1}.
Surprisingly, the simple-minded averaging approach performs extremely well; this is consistent with the previously reported results and the success of the methods of \cite{Cai2011} (essentially harmonic mean) and \cite{MaAnchor2008} (arithmetic mean).

\begin{table}[htdp]
\begin{center}
\begin{tabular}{rccccccc}
                        & \small{\bf Clus.} & \small{\bf SC-1} & \small{\bf SC-2} & \small{\bf JADE} & \small{\bf JD-HM} & \small{\bf Comraf} \\
\hline
\small{\em Blobs}	& \small6 & \small91.0$\pm$7.2\% & \small90.8$\pm$7.2\% & \small97.3$\pm$4.2\% & \small98.3$\pm$3.0\% & \small86.9$\pm$8.6\%   \\
\small{\em Circles}	& \small4 & \small65.9\% & \small63.4\% & \small{100.0\%} & \small99.8\% & \small31.4\% \\
\small{\em NIPS}	& \small4 & \small63.3\% & \small75.1\% & \small99.9\% & \small99.9\% & \small51.8\% \\
\hline
\small{\em NUS}	& \small7 & \small83.5\% & \small71.0\% & \small92.4\% & \small80.7\% & \small82.1\%  \\
\hline
 \multirow{2}{*}{\small{\em Caltech}}
	& \small7 & \small73.3\% & \small76.2\% & \small86.7\% & \small84.8\% &  -- \\
	& \small20 & \small66.3\% & \small70.7\% & \small73.3\% & \small76.0\% & -- \\
\hline
\end{tabular}
\end{center}\vspace{-3mm}
\caption{\label{tab:clustering1} Accuracy of different clustering methods.}
\end{table}

\begin{figure}[ht]

\begin{minipage}[b]{0.23\linewidth}
\centering
\includegraphics[width=0.9\linewidth]{figures/figure_parts/figureRandom_1_L1.pdf}
\end{minipage}\hspace{1mm}
\begin{minipage}[b]{0.23\linewidth}
\centering
\includegraphics[width=0.9\linewidth]{figures/figure_parts/figureRandom_2_L2.pdf}
\end{minipage}\hspace{6mm}
\begin{minipage}[b]{0.23\linewidth}
\centering
\includegraphics[width=0.9\linewidth]{figures/figure_parts/figureRandom_1_cardoso.pdf}
\end{minipage}\hspace{1mm}
\begin{minipage}[b]{0.23\linewidth}
\centering
\includegraphics[width=0.9\linewidth]{figures/figure_parts/figureRandom_2_cardoso.pdf}
\end{minipage}\vspace{2mm}
\begin{minipage}[b]{0.23\linewidth}
\centering
\includegraphics[width=0.8\linewidth]{figures/figure_parts/figureCircles_1_L1.pdf}
\end{minipage}\hspace{1mm}
\begin{minipage}[b]{0.23\linewidth}
\centering
\includegraphics[width=0.8\linewidth]{figures/figure_parts/figureCircles_2_L2.pdf}
\end{minipage}\hspace{6mm}
\begin{minipage}[b]{0.23\linewidth}
\centering
\includegraphics[width=0.8\linewidth]{figures/figure_parts/figureCircles_1_wa.pdf}
\end{minipage}\hspace{1mm}
\begin{minipage}[b]{0.23\linewidth}
\centering
\includegraphics[width=0.8\linewidth]{figures/figure_parts/figureCircles_2_wa.pdf}
\end{minipage}\vspace{1mm}
\begin{minipage}[b]{0.23\linewidth}
\centering
\includegraphics[width=0.9\linewidth]{figures/figure_parts/figureNIPS_1_L1.pdf}
\end{minipage}\hspace{1mm}
\begin{minipage}[b]{0.23\linewidth}
\centering
\includegraphics[width=0.9\linewidth]{figures/figure_parts/figureNIPS_2_L2.pdf}
\end{minipage}\hspace{6mm}
\begin{minipage}[b]{0.23\linewidth}
\centering
\includegraphics[width=0.9\linewidth]{figures/figure_parts/figureNIPS_1_wa.pdf}
\end{minipage}\hspace{1mm}
\begin{minipage}[b]{0.23\linewidth}
\centering
\includegraphics[width=0.9\linewidth]{figures/figure_parts/figureNIPS_2_wa.pdf}
\end{minipage}\vspace{-1mm}
\begin{minipage}[b]{0.23\linewidth}
\centering
\small Modality 1
\end{minipage}\hspace{3mm}
\begin{minipage}[b]{0.23\linewidth}
\centering
\small Modality 2
\end{minipage}\hspace{1mm}
\begin{minipage}[b]{0.46\linewidth}
\centering
\small Joint
\end{minipage}\vspace{-3mm}
\caption{\label{fig:clustering1} Clustering synthetics datasets.
Marker size represents ground truth; marker color represents segmentation results (ideally, markers of each type should have a single color).}
\end{figure}

{\bf NUS dataset. }
In the third experiment, we used a subset of the NUS-WIDE dataset \cite{nus-wide-civr09} containing annotated images. The images were selected on purpose to have ambiguous content and annotations (e.g., swimming tigers are also tagged as ``water'' making them confuse e.g. with whales).
As two different modalities, we used the 64-dimensional color histograms and 1000-dimensional bags of words.
Laplacians were constructed using 10 nearest neighbors and Gaussian weight was selected using self-tuning.
Table~\ref{tab:clustering1} shows the performance of different clustering methods, and Figure~\ref{fig:nus} exemplifies the clustered images.

Using JADE joint diagonalization, we produced all the joint eigenvectors of the two modalities Laplacians.
Figure~\ref{fig:roc} (top) shows the distance matrices between the objects in the NUS dataset obtained using uni- and multi-modal diffusion distances (computed with the first $100$ eigenvectors according to~(\ref{eq:diffdist}) using heat diffusion kernel $K(\lambda) = e^{-5\lambda}$).
Ideally, the distance matrix should contain zero blocks on the diagonal (objects of the same class) and non-zero elsewhere (objects from different classes). Thresholding these distances at a set of levels and measuring the false positives/true positive rates (FPR/TPR), we produce the ROC curves that clearly indicate the advantage of using multiple modalities (Figure~\ref{fig:roc}).

In Figure~\ref{fig:fps} (top), we used the diffusion distance to progressively sample the NUS dataset using the farthest point sampling strategy: starting with some point, pick up the second one as most distant from the first; then the third as the most distant from the first and second, and so on. Such sampling is almost-optimal (\cite{hochbaum1985best}) and is known to produce a progressively refined $r$-covering of the set.
In fact, the first $7$ samples produced in this way cover all the classes present in the dataset, which is an indication of the meaningfulness of such a sampling.

{\bf Caltech dataset. }
In the fourth experiment, we repeated the third experiment on a subset of the Caltech-101 dataset with 7 and 20 image classes as in \cite{Cai2011}.
For each image, kernels arising from different visual descriptors were given.
For the 7-clusters experiment, we used the bio-inspired features and 4x4 pyramid histogram of visual words (PHOW); for the 20-clusters experiment, we used geometric blur and 4x4 PHOW descriptors as different modalities, respectively.
Laplacians were constructed from these kernels using
Gaussian weight selected with self-tuning.
Diffusion distances were computed with the first $100$ eigenvectors using the kernel $K(\lambda) = e^{-5\lambda}$.
The results are shown in Figures~\ref{fig:cal}--\ref{fig:roc}.

\begin{figure}[t!]
\centering
\includegraphics[width=1\linewidth]{figures/NUS/nus_clustering_cut_full.pdf} \vspace{-7mm}
\caption{\label{fig:nus}
Spectral clustering of NUS dataset. Shown are a few images and corresponding tags belonging to the same cluster obtained using Tags (top row), Color histogram (second row), and joint modalities (third to fifth row). Groundtruth clusters are shown in different colors. }
\end{figure}

\begin{figure}[t!]
\centering
\includegraphics[width=1\linewidth]{figures/caltech/caltech_clustering_full.pdf} \vspace{-7mm}
\caption{\label{fig:cal}
Spectral clustering of Caltech101 dataset. Shown are a few images and corresponding tags belonging to the same cluster obtained using ht\_bio\_105034 bio-inspired features (top row), 4x4 PHOW (second row), and joint modalities (third and fourth row). Groundtruth clusters are shown in different colors.  }
\end{figure}

\begin{figure}[t!]
\centering
\includegraphics[width=1\linewidth]{figures/diffdist/fps-nus-10.pdf}\vspace{1mm} \\
\includegraphics[width=1\linewidth]{figures/diffdist/fps-caltech-10.pdf} \vspace{-7mm}
\caption{\label{fig:fps}
Farthest point sampling of NUS (top) and Caltech (bottom) datasets using joint diffusion distance.  First point is on the left. Numbers indicate the sampling radius. Note that in both cases, the first 7 samples cover all the image classes. }
\end{figure}

\begin{figure}[t!]
\centering
\includegraphics[width=1\linewidth]{figures/diffdist/roc-nus-1.pdf}\vspace{0mm}\\
\includegraphics[width=1\linewidth]{figures/diffdist/roc-caltech-1.pdf} \vspace{-8mm}
\caption{\label{fig:roc}
Columns one to three: distance matrices, column four: ROC curves, computed on NUS (top) and Caltech (bottom) datasets using joint diffusion distance.  Ambiguities are shown in white. }
\end{figure}