# Diffusion and Volume Maximization-Based Clustering of Highly Mixed
Hyperspectral Images

Sam L. Polk  Department of Mathematics, Tufts University  Kangning Cui
Department of Mathematics, City University of Hong Kong  Robert J. Plemmons
Departments of Mathematics and Computer Science, Wake Forest University  James
M. Murphy111Corresponding Author<EMAIL_ADDRESS>Department of
Mathematics, Tufts University

###### Abstract

Hyperspectral images of a scene or object are a rich data source, often
encoding a hundred or more spectral bands of reflectance at each pixel.
Despite being very high-dimensional, these images typically encode latent low-
dimensional structure that can be exploited for material discrimination.
However, due to an inherent trade-off between spectral and spatial resolution,
many hyperspectral images are generated at a coarse spatial scale, and single
pixels may correspond to spatial regions containing multiple materials. This
article introduces the _Diffusion and Volume maximization-based Image
Clustering_ (_D-VIC_) algorithm for unsupervised material discrimination.
D-VIC locates cluster modes—high-density, high-purity pixels in the
hyperspectral image that are far in diffusion distance (a data-dependent
distance metric) from other high-density, high-purity pixels—and assigns these
pixels unique labels, as these points are meant to exemplify underlying
material structure. Non-modal pixels are labeled according to their diffusion
distance nearest neighbor of higher density and purity that is already
labeled. By directly incorporating pixel purity into its modal and non-modal
labeling, D-VIC upweights pixels that correspond to a spatial region
containing just a single material, yielding more interpretable clusterings.
D-VIC is shown to outperform baseline and comparable state-of-the-art methods
in extensive numerical experiments on a range of hyperspectral images,
implying that it is well-equipped for material discrimination and clustering
of these data.

Index Terms: Clustering, Diffusion Geometry, Graphs, Hyperspectral Imaging,
Image Segmentation,
Machine Learning, Spectral Unmixing, Unsupervised Learning.

##  1 Introduction

Hyperspectral images (HSIs) are images of a scene or object, encoding
reflectance at a hundred or more spectral bands per pixel [Eis12, GPC+17,
PMP+11]. HSI remote sensing data, which is generated continuously by airborne
and orbital sensors, has been successfully used for material discrimination in
fields ranging from forensic medicine (e.g., age estimation of forensic traces
[EGVL+12]), conservation (e.g., species mapping in wetlands [AMR10, HMW03]),
and ecology (e.g., estimating water content in vegetation canopies [CKS10]).
The high dimensional characterization of a scene provided in remote sensing
HSI data has motivated its use in the _material discrimination_ problem
[DBG08], wherein machine learning algorithms are meant to separate pixels
based on the constituent materials in the spatial region they correspond to
[GPC+17, PMP+11].

Though hyperspectral imagery has become an essential tool across many
scientific domains, material discrimination using HSI data faces at least two
key difficulties. First, because of an inherent trade-off between spectral and
spatial resolution, HSIs are typically generated at coarse spatial resolution
[JWC+20, Pri97, BDPCV+13, LSR15]. One would prefer an HSI with both a high
spatial resolution (so that individual pixels correspond to spatial regions
containing just one material) and a high spectral resolution (to enable
capacity for material discrimination) [Pri97]. However, an increase in the
spatial resolution of an HSI often comes at the cost of reducing the effective
detection energy entering the recording spectrometer across each spectral band
[JWC+20]. While this effect may be mitigated at least partially by increasing
the aperture of the optical system underlying the spectrometer used to
generate HSI data [JWZ+17], high-aperture instruments generally also have high
volume and weight [JWC+20]. As such, HSI data is generated at a coarse spatial
resolution (often as high as 20 m per pixel). Thus, though some _high-purity_
pixels in an HSI may correspond to spatial regions containing predominantly
just one material, other pixels are mixed: corresponding to spatial regions
with many distinct materials.

Second, expert labeling to generate training data for supervised machine
learning is often impractical due to the large quantities of HSI data
continuously generated by remote sensors [Eis12]. To efficiently analyze
unlabeled HSIs, one may use HSI clustering algorithms, which partition HSI
pixels into groups of points sharing key commonalities [FHT01]. These
algorithms are _unsupervised_ i.e., ground truth labels are not needed to
provide a partition of an HSI [FHT01]. Clustering is usually one of the first
tools available to a user interested in exploring a new dataset. Though
clustering has become an important tool in the field of hyperspectral imagery
[MM18, AG21, ZWL+16, KDP12, WNW+19, CVMZ07, GJC+14, WP17, YYL+18, NMI+20,
GKP14, LQL+21, SWW+20], clustering algorithms that do not directly account for
the fact that HSI pixels are often spectrally mixed may fail to extract
meaningful latent cluster structure. Thus, efficient unsupervised machine
learning algorithms incorporating pixel purity into their labeling are needed
to analyze the rich material structure encoded in HSI data [CP21].

This article introduces the _Diffusion and Volume maximization-based Image
Clustering_ (_D-VIC_) algorithm for unsupervised material discrimination of
HSIs. D-VIC is the first algorithm to simultaneously exploit the high-
dimensional geometry [BAF05, AG21] and abundance structure [CP21, BDPCV+13]
observed in HSIs for the clustering problem. In its first stage, D-VIC locates
cluster modes: high-purity, high-empirical density pixels that are far in
diffusion distance (a data-dependent distance metric [CL06]) from other high-
purity, high-density pixels. These pixels serve as exemplars for all
underlying material structure in the HSI. In its mode selection, D-VIC
downweights high-density pixels that correspond to commonly co-occuring groups
of materials. As such, D-VIC’s exploitation of spectrally mixed structure in
HSI data [JWC+20, Pri97, BDPCV+13, LSR15, BDPCV+13] enables the selection of
modes that better represent the material structure in the scene.

After detecting cluster modes, D-VIC propagates modal labels to non-modal
pixels in order of decreasing density and pixel purity. Since pixel purity is
also incorporated into D-VIC’s non-modal labeling, D-VIC accounts for material
abundance structure in the HSI during its entire labeling procedure. D-VIC is
compared against classical and related state-of-the-art HSI clustering
algorithms on three benchmark real HSI datasets and one synthetic datset built
to demonstrate the key strengths of D-VIC. On each dataset, D-VIC produces a
partition closer to the ground truth than related methods, and moreover enjoys
robustness to parameter selection. Computationally, the proposed algorithm
scales quasilinearly in the size of the HSI, and the empirical runtime of
D-VIC is competitive, suggesting it is well-equipped to cluster large HSIs.

##  2 Background

###  2.1 Background on Unsupervised HSI Clustering

HSI clustering algorithms partition an HSI, denoted $X=\leavevmode\nobreak\
\\{x_{i}\\}_{i=1}^{n}\subset\leavevmode\nobreak\ \mathbb{R}^{D}$ (interpreted
as a point cloud of spectral signatures, with $n$ pixels and $D$ spectral
bands) into $K$ _clusters_ of pixels. The partition, which we call a
_clustering_ of $X$, may be encoded in a labeling vector
$\mathcal{C}\in\\{1,2,\dots,K\\}^{n}$ such that
$\mathcal{C}(x_{i})\leavevmode\nobreak\ =\leavevmode\nobreak\
\mathcal{C}_{i}\in\leavevmode\nobreak\ \\{1,2,\dots,K\\}$ is the label
assigned to the pixel $x_{i}$. Ideally, pixels from any one cluster are in
some sense “related,” and pixels from any two clusters are “unrelated.” The
exact notion of relatedness varies widely in the clustering literature [FHT01,
NJW02, MM19b, GKP14]. Clustering algorithms are _unsupervised_ , meaning that
data points are labeled without the aid of any expert annotations or ground
truth labels. This has motivated the development of clustering algorithms
explicitly built for material discrimination using HSIs [MM18, AG21, ZWL+16,
KDP12, WNW+19, CVMZ07, GJC+14, WP17, YYL+18, NMI+20, GKP14, LQL+21, SWW+20,
CCM14].

Though classical clustering algorithms such as $K$-Means and the Gaussian
Mixture Model (GMM) [FHT01] remain widely used in practice, these algorithms
tend to perform poorly on HSIs for a number of reasons [MM18, AG21]. First,
algorithms that rely on Euclidean distances are prone to the “curse of
dimensionality” on datasets like HSIs that have a high ambient dimension
(i.e., the number of spectral bands is large) [TK06]. Second, HSIs are often
spectrally mixed [JWC+20, Pri97, BDPCV+13, LSR15], and overlap may exist
between clusters in Euclidean space [MM18]. A final complication is that
classical algorithms generally assume that latent clusters in a dataset are
approximately ellipsoidal groups of points that are well-separated in
Euclidean space [FHT01], but clusters in HSIs often exhibit nonlinear
structure. A simple toy dataset, visualized in Figure 1 [NJW02], serves as an
example of dataset with nonlinear structure. This dataset lacks a linear
decision boundary between its $K=2$ latent clusters, and classical algorithms
($K$-Means and GMM [FHT01]) were unable to learn its latent nonlinear cluster
structure. HSIs often contain clusters that can only be separated using a
nonlinear decision boundary [MM18]; thus, algorithms that rely solely on
Euclidean distances are expected to perform poorly at material discrimination
on HSIs.

Figure 1: Example of a toy dataset $(n=1000)$ with nonlinear cluster
structure: two interleaving half-circles. The idealized clustering, visualized
in (a), separates each half-circle. Due to the lack of a linear decision
boundary, however, both $K$-Means ((b)) and GMM ((c)) were unable to extract
latent cluster structure from this simple nonlinear dataset. Both algorithms
were run with 100 replicates.

The limitations outlined above have motivated the application and development
of nonlinear graph-based algorithms for HSI clustering [AG21, ZWL+16, KDP12,
ZCT+17, MM18, WNW+19, CVMZ07, GJC+14, NJW02, CCM14]. Graph-based algorithms
rely on data-generated graphs; pixels are represented as nodes in the graph,
and edges encode pairwise similarity between them. Highly connected regions in
the graph may then be summarized using a nonlinear coordinate transformation
[TdSL00, RS00, BN01, CL06], as is described in more detail in Section 2.2.
Thus, a partition may be obtained by implementing a classical clustering
algorithm on the dimension-reduced dataset. Graph-based clustering algorithms
enjoy a strong theoretical backing that has motivated their development for
HSI clustering [AG21, ZWL+16, KDP12, MM18, WNW+19, CVMZ07, GJC+14, CCM14]. For
example, due to their reliance on a graph representation of an HSI, these
algorithms tend to be robust to small perturbations in the data and noise.
Moreover, theoretical guarantees exist for the successful recovery of latent
cluster structure, even if boundaries between latent clusters are nonlinear
[RCY11, MP22, MM19b]. Despite their exhibited successes, algorithms that rely
solely on graph structure tend to perform poorly on datasets containing
multimodal cluster structure [MM19b, MP22, NG07]; i.e., if a single cluster
has multiple regions of high and low density. Importantly, this includes
spectrally mixed HSIs, the classes of which often contain multiple co-occuring
materials of varying abundances [JWC+20, Pri97, BDPCV+13].

Deep neural networks have recently become popular for material discrimination
in HSIs because of their capacity for prediction on complex data sources like
HSI [WP17, YYL+18, NMI+20, LQL+21, SWW+20]. While these algorithms tend to be
highly accurate on real HSI data, many state-of-the-art deep models for HSI
segmentation still rely on at least a small number of ground truth labels
[WP17, YYL+18] and/or pseudo-labels [LQL+21, WP17]. Moreover, even fully
unsupervised “deep clustering” algorithms [NMI+20, SWW+20] rely on deep neural
networks, which are prone to error from perturbations and noise [NYC15,
SZS+14], and whose strong performance in some settings may be incidental to
the use of deep features [HYV20].

###  2.2 Background on Spectral Graph Theory

As overviewed in Section 2.1, graph-based clustering algorithms learn latent,
possibly nonlinear cluster structure from HSIs by treating pixels as nodes in
an undirected, weighted graph, where connections between pixels are encoded in
a weight matrix $\mathbf{W}\in\mathbb{R}^{n\times n}$ [NJW02, MM19b, MP22,
PM21]. In large datasets like HSIs, edges can be restricted to the first $N\ll
n$ $\ell^{2}$-nearest neighbors (i.e. Euclidean nearest neighbors) and given
unit weight. In other words, $\mathbf{W}_{ij}=1$ if $x_{i}$ is one of the $N$
nearest neighbors of $x_{j}$ or vice versa, and $\mathbf{W}_{ij}=0$ otherwise.
Let $\mathbf{P}=\mathbf{D}^{-1}\mathbf{W}$, where
$\mathbf{D}\in\mathbb{R}^{n\times n}$ is the diagonal degree matrix defined by
$\mathbf{D}_{ii}=\sum_{j=1}^{n}\mathbf{W}_{ij}$. The matrix
$\mathbf{P}\in\mathbb{R}^{n\times n}$ may be interpreted as the transition
matrix for a Markov diffusion process on $X$ and has a unique stationary
distribution $\pi\in\mathbb{R}^{1\times n}$ satisfying $\pi\mathbf{P}=\pi$
[CL06, MM19b]. Define $\\{(\lambda_{i},\psi_{i})\\}_{i=1}^{n}$ to be the
(right) eigenvalue-eigenvector pairs of $\mathbf{P}$, sorted in non-increasing
order so that $1=\lambda_{1}>|\lambda_{2}|>\dots>|\lambda_{n}|\geq 0$. The
first $K$ eigenvectors of $\mathbf{P}$ often concentrate on the $K$ most
coherent subgraphs in the graph underlying $\mathbf{P}$, making these vectors
useful for clustering [NJW02].

####  2.2.1 Background on Spectral Clustering

The Spectral Clustering (SC) algorithm is the implementation of $K$-Means on
the first $K$ eigenvectors of $\mathbf{P}$, usually after a normalization step
[NJW02]. Though SC can often learn latent cluster structure from datasets with
nonlinear and elongated clusters, it tends to perform poorly on datasets like
HSIs that contain multimodal cluster structure [NG07, MM18, MM19b].

####  2.2.2 Background on Diffusion Geometry

Diffusion distances are a family of data-dependent distance metrics that
enable comparisons between points in the context of the Markov diffusion
process encoded in $\mathbf{P}$ [CL06]. Diffusion distances have been
successfully used in a number of applications (e.g., in gene expression
profiling [HBW+16, VDSN+18], data visualization [ZS14, MvDW+19], and molecular
dynamics analysis [RZMC11, ZRMC11, CF18]). Moreover, diffusion distances have
been shown to efficiently capture low-dimensional structure in HSI data,
resulting in excellent clustering performance [MM18, MP22].

Define the _diffusion distance_ at time $t\geq 0$ between pixels
$x_{i},x_{j}\in X$ [CLL+05, CL06, NLCK06] by

$D_{t}(x_{i},x_{j})=\sqrt{\sum_{k=1}^{n}\frac{[(\mathbf{P}^{t})_{ik}-(\mathbf{P}^{t})_{jk}]^{2}}{\pi_{k}}}.$

Diffusion distances are a nonlinear data-dependent distance metric that have a
natural connection to the clustering problem [MM19b, MP22]. To see this, note
that $D_{t}(x_{i},x_{j})$ may be interpreted as the Euclidean distance between
the $i$th and $j$th rows of $\mathbf{P}^{t}$, weighted according to $1/\pi$.
If pixels from the same cluster share many high-weight paths of length $t$,
but paths of length $t$ between any two pixels from different clusters are
relatively low weight, then the $i$th and $j$th rows of $\mathbf{P}^{t}$ are
expected to be nearly equal for pixels $x_{i}$ and $x_{j}$ from the same
cluster and very different if these pixels come from different clusters.
Indeed, in datasets with coherent and well-separated latent cluster structure,
the maximum within-cluster diffusion distance is bounded away from the minimum
between-cluster diffusion distance across a large range of $t$ [MM19b, MP22].

Diffusion distances can be efficiently computed using the eigendecomposition
of $\mathbf{P}$ [CL06, NLCK06, CLL+05]:

$D_{t}(x_{i},x_{j})=\sqrt{\sum_{k=1}^{n}\lambda_{k}^{2t}[(\psi_{i})_{k}-(\psi_{j})_{k}]^{2}}.$

For $t$ sufficiently large so that $|\lambda_{k}|^{2t}\approx 0$ for $k>\ell$,
the sum in diffusion distances can be truncated past the $\ell$th term,
yielding an accurate and efficient approximation of diffusion distances.
Importantly, the relationship between diffusion distances and the
eigendecomposition of $\mathbf{P}$ indicates that diffusion distances may be
interpreted as Euclidean distances after nonlinear dimensionality reduction
via the following dimension-reduced mapping of the ambient space into
$\mathbb{R}^{\ell}$:
$x_{i}\rightarrow[\lambda_{1}^{t}(\psi_{1})_{i}\;\lambda_{1}^{t}(\psi_{2})_{i}\;\dots,\lambda_{1}^{t}(\psi_{\ell})_{i}]$
[CL06, CLL+05, NLCK06]. As such, diffusion distances are related to nonlinear
dimensionality reduction; e.g., locally linear embeddings [RS00], isomap
[TdSL00], and Laplacian eigenmaps [BN01].

HSIs often encode well-defined latent multiscale cluster structure that can be
learned by diffusion-based HSI clustering algorithms by varying the time
parameter $t$ in diffusion distances [PM21, MP22]. Indeed, smaller $t$ values
generally enable the detection of fine-scale local cluster structure, while
larger $t$ values enable the detection of coarse-scale global cluster
structure. However, for algorithms that require $K$ as an input, $t$ must be
tuned to correspond to the desired number of clusters. Thus, the choice of $t$
must be carefully considered when clustering a dataset using an algorithm that
relies on diffusion distances [MP22, PM21].

###  2.3 Background on the LUND Algorithm

The Learning by Unsupervised Nonlinear Diffusion (LUND) clustering algorithm
(Algorithm 1) leverages the natural theoretical relationship between diffusion
geometry and the clustering problem to efficiently partition datasets [CL06,
MM19b, MM18]. LUND relies on two main functions to calculate cluster modes:
points acting as exemplars for all underlying cluster structure. The first is
a _kernel density estimate_ (_KDE_) defined by

$p(x)=\frac{1}{Z}\sum_{y\in NN_{N}(x)}\exp(-\|x-y\|^{2}_{2}/\sigma_{0}^{2}),$

where $NN_{N}(x)$ is the set of $N$ $\ell^{2}$-nearest neighbors of $x$ in
$X$, $\sigma_{0}>0$ is a _KDE scale_ controlling the interaction radius
between points, and $Z$ is a constant normalizing $p(x)$ so that $\sum_{y\in
X}p(y)=1$. By construction, $p(x)$ will be large if the pixel $x$ is close to
its $N$ $\ell^{2}$-nearest neighbors in $X$ and small otherwise. The second
function used by LUND for mode detection is

$\displaystyle\rho_{t}(x)=\begin{cases}\max\limits_{y\in
X}D_{t}(x,y)&x=\operatorname*{argmax}\limits_{y\in X}p(y),\\\
\min\limits_{y\in X}\\{D_{t}(x,y)\ |\ p(y)\geq
p(x)\\}&\text{otherwise}.\end{cases}$

For points that are not empirical density maximizers, $\rho_{t}(x)$ is the
diffusion distance at time $t$ between the pixel $x$ and its $D_{t}$-nearest
neighbor that is higher density.

Input: $X$ (HSI), $N$ (# nearest neighbors),

$\sigma_{0}$ (KDE scale), $t$ (diffusion time), $K$ (# clusters)

Output: $\mathcal{C}$ (clustering)

1 Build the KDE $p(x)$ using $N$ $\ell^{2}$-nearest neighbors and KDE scale
$\sigma_{0}$;

2 Build $\rho_{t}(x)$, where diffusion distances are computed from a KNN graph
with $N$ edges per pixel;

3 Assign $\mathcal{C}(x_{m_{k}})=k$ for $1\leq k\leq K$, where
$\\{x_{m_{k}}\\}_{k=1}^{K}$ are the $K$ pixels maximizing
$\mathcal{D}_{t}(x)=p(x)\rho_{t}(x)$;

In order of non-increasing $p$-value, for each unlabeled $x\in X$, assign the
label $\mathcal{C}(x)=\mathcal{C}(x^{*})$, where
$x^{*}=\operatorname*{argmin}\limits_{y\in X}\\{D_{t}(x,y)|\;p(y)\geq
p(x)\text{ and }\mathcal{C}(y)>0\\}$.

Algorithm 1 Learning by Unsupervised Nonlinear Diffusion (LUND) [MM19b]

Hence, the maximizers of $\mathcal{D}_{t}(x)=p(x)\rho_{t}(x)$, being high-
density data points far in diffusion distance at time $t$ from other high-
density data points [MM19b], are reasonable choices for cluster exemplars.
LUND designates the $K$ points maximizing $\mathcal{D}_{t}(x)$ as _cluster
modes_ and assigns these points unique labels. LUND then assigns each non-
modal point (in order of non-increasing density) the label of its
$D_{t}$-nearest neighbor that is higher density and already labeled. Under lax
assumptions on density and the graph underlying $\mathbf{P}$, LUND is
guaranteed to learn latent cluster structure, even in datasets with nonlinear,
elongated, and/or multimodal structure. Due to its success on real and
synthetic datasets, many modifications of LUND exist [MM19b, MP22, MM18, PM21,
MM19c, MM19a, Mur20], some built specifically for HSI clustering [MP22, PM21,
MM19c, MM18].

Though LUND has strong performance guarantees for recovering latent cluster
structure on general datasets [MM19b, MP22], it may be subject to error on
HSIs with highly mixed pixels. For such datasets, the spatial region
corresponding to a single pixel may contain multiple materials. As such,
though LUND may be equipped to segment an HSI into clusters of pixels that
contain the same groups of multiple materials, it may be unable to locate
individual materials within a mixed HSI.

###  2.4 Background on Spectral Unmixing

Real-world HSI data is often generated at a coarse spatial resolution; thus,
pixels may correspond to spatial regions containing multiple materials
[CMAC11, Win99, MSS01]. To learn latent material structure from HSIs,
_spectral unmixing_ algorithms decompose each pixel’s spectrum into a linear
combination of _endmembers_ that encode the spectral signature of materials in
the scene. The coefficients of this linear combination indicate the relative
_abundance_ of materials within the spatial region corresponding to that
pixel. Mathematically, a spectral unmixing algorithm learns
$\mathbf{U}=(u_{1}\;u_{2}\;\dots\;u_{m})^{T}\in\mathbb{R}^{m\times D}$ (with
rows encoding the spectral signatures of endmembers) and
$\mathbf{A}\in\mathbb{R}^{n\times m}$ (with rows encoding abundances) such
that $x_{i}\approx\sum_{j=1}^{m}\mathbf{A}_{ij}u_{i}$ for each $x_{i}\in X$
[MSS01]. Usually, $\mathbf{A}$ is normalized so that
$\sum_{j=1}^{m}\mathbf{A}_{ij}=1$ for each $i$; hence, abundances are data-
dependent features storing estimates for the relative frequency of materials
in pixels. The _purity_ of $x_{i}\in X$, defined by $\eta(x_{i})=\max_{1\leq
i\leq m}\mathbf{A}_{ij}$ [CP21], will be large if the spatial region
corresponding to $x_{i}$ is highly homogeneous (i.e., containing predominantly
just one material) and small otherwise. As such, pixel purity and spectral
unmixing may be used to aid in the unsupervised clustering of HSIs [CP21,
GKP14, KDP12].

Spectral unmixing has become an important tool within hyperspectral imagery,
prompting its usage in a number of applications; e.g., image reconstruction
[WFG+18], noise reduction [CMR13], and anomaly detection [QWG+18, MYW18]. The
importance of spectral unmixing in remote sensing has motivated the
development of many algorithms for this task [BDPCV+13]. Geometric methods for
spectral unmixing estimate endmembers by searching for endmembers that form a
simplex of minimal volume, subject to a constraint that at least some nearly
pure pixels exist within the observed HSI pixels [CWLO06, Nev99, BKG95, Boa93,
CCHM09, AG21, BDPCV+13, CMAC11, ND05]. For highly mixed HSIs that lack pure
pixels, statistical methods may be used [BKL+04, ZG07, DMC+09, MBMDC06,
TRK11]. These methods typically treat spectral unmixing as a blind source
separation problem, and though they are often successful at this task,
statistical algorithms are usually more computationally expensive [BDPCV+13].
Additionally, while linear spectral unmixing is well-developed and widely used
in practice, some nonlinear unmixing algorithms (including relying on neural
networks [GAC01, LDF11, COR11, SLP+19, ZSZ+18]) have been developed to account
for nonlinear interactions between endmembers [GAC01, LDF11, SLP+19, ZSZ+18,
YCI13, HADT11, CRH12, OKA18]. Nevertheless, these algorithms typically require
training data or hyperparameter inputs unlike many of the linear mixing models
reviewed above [BDPCV+13].

####  2.4.1 Background on the HySime Algorithm

Hyperspectral Signal Subspace Identification by Minimum Error (_HySime_) is a
standard algorithm for estimating the number of materials $m$ in $X$ [BDN08].
HySime assumes that each $x_{i}\in X$ is of the form $x_{i}=y_{i}+\zeta_{i}$,
where $y_{i}\in\mathbb{R}^{D}$ and $\zeta_{i}\in\mathbb{R}^{D}$ model the
signal and noise associated with $x_{i}$, respectively. If signal vectors are
linear mixtures of $m$ ground truth endmembers (i.e.,
$y_{i}=\sum_{j=1}^{m}\mathbf{A}_{ij}u_{j}$ for $1\leq i\leq m$), then the set
$\\{y_{i}\\}_{i=1}^{m}$ lies on a $m$-dimensional subspace of
$\mathbb{R}^{D}$. With this motivation, HySime estimates the subspace
dimension by balancing the error of projecting signal vectors
$\\{y_{i}\\}_{i=1}^{n}$ onto their first $m$ principal components and the
amount of noise captured by those vectors’ orthogonal space. Though other
algorithms exist for estimating the number of materials in a scene using HSI
data, many of these alternatives rely on hyperparameter inputs to estimate $m$
or have large computational complexity [Cha18]. In contrast, HySime can
compute a high-quality estimate using only $X$ in just $O(D^{2}n)$ operations.

####  2.4.2 Background on the Alternating Volume Maximization

Alternating Volume Maximization (_AVMAX_) is a spectral unmixing algorithm,
requiring $m$ as a parameter, that searches for vectors
$\\{u_{i}\\}_{i=1}^{m}\subset\mathbb{R}^{D}$ that produce an $m$-simplex of
maximal volume, subject to the constraint that each $u_{i}$ lies in the convex
hull of the dataset after PCA dimensionality reduction. Endmembers are
optimized through multiple partial maximization procedures (i.e., keeping
$m-1$ endmembers constant and optimizing for volume while varying the
$m\textsuperscript{th}$ endmember) until convergence [CMAC11]. AVMAX has
recently become popular for spectral unmixing because of its strong
performance guarantees and the rigor behind its optimization framework
[CMAC11, BDPCV+13]. In a noiseless, linearly mixed dataset containing the
optimal endmember set, if each partial maximization problem in AVMAX converges
to a unique solution, AVMAX is guaranteed to converge to the optimal endmember
set [CMAC11]. Moreover, AVMAX can easily be modified to make it robust to
random initialization; one can run multiple replicates of AVMAX in parallel
and choose the endmember set with largest volume. Once endmembers are learned,
abundances can be computed using a nonnegative least squares solver
$(\mathbf{A}_{i1}\;\mathbf{A}_{i2}\;\dots\;\mathbf{A}_{im})=\operatorname*{argmin}_{a\in\mathbb{R}^{m}}\|x_{i}-\sum_{j=1}^{m}a_{j}u_{j}\|_{2}^{2}$
for each $x_{i}\in X$ [BDJ97].

##  3 Diffusion and Volume maximization-based Image Clustering

In highly mixed HSIs, any one pixel may correspond to a spatial region that
contains many materials. Thus, even state-of-the-art algorithms for
unsupervised material discrimination may perform poorly on mixed HSIs, failing
to recover clusterings that can be linked to materials within the scene.
Algorithms that do not directly incorporate a spectral unmixing step into
their labeling may assign clusters that correspond to groups of materials
rather than clusters corresponding to individual materials. Thus, additional
improvements are needed to improve HSI clustering algorithms for material
discrimination on highly mixed images.

Input: $X$ (HSI), $N$ (# nearest neighbors),

$\sigma_{0}$ (KDE scale), $t$ (diffusion time), $K$ (# clusters)

Output: $\mathcal{C}$ (clustering)

1 Estimate $m$, the number of latent endmembers in $X$, using HySime [BDN08];

2 Learn endmembers $\mathbf{U}\in\mathbb{R}^{m\times D}$ using AVMAX [CMAC11]
and abundances $\mathbf{A}\in\mathbb{R}^{n\times m}$ using a nonnegative least
squares solver [BDJ97];

3 For each $x\in X$, calculate pixel purity $\eta(x)$ and empirical density
$p(x)$. Store
$\zeta(x)=\frac{2\bar{p}(x)\bar{\eta}(x)}{\bar{p}(x)+\bar{\eta}(x)}$, where
$\bar{p}(x)=\frac{p(x)}{\max_{y\in X}p(y)}$ and
$\bar{\eta}(x)=\frac{\eta(x)}{\max_{y\in X}\eta(y)}$;

4 Build $d_{t}(x)$, where diffusion distances are computed from a KNN graph
with $N$ edges per pixel;

5 Assign $\mathcal{C}(x_{m_{k}})=k$ for $1\leq k\leq K$, where
$\\{x_{m_{k}}\\}_{k=1}^{K}$ are the $K$ pixels maximizing
$\mathcal{D}_{t}(x)=\zeta(x)d_{t}(x)$;

In order of non-increasing $\zeta(x)$, for each unlabeled $x\in X$, assign $x$
the label $\mathcal{C}(x^{*})$, where
$x^{*}=\operatorname*{argmin}\limits_{y\in
X}\\{D_{t}(x,y)|\;\zeta(y)\geq\zeta(x)\text{ and }\mathcal{C}(y)>0\\}$.

Algorithm 2 Diffusion and Volume maximization-based Image Clustering (D-VIC)

This section introduces the _Diffusion and Volume maximization-based Image
Clustering_ (_D-VIC_) algorithm (Algorithm 2) for unsupervised material
discrimination of HSIs. To learn material abundances, D-VIC first performs a
spectral unmixing step: decomposing the HSI by learning the number of
endmembers $m$ using HySime [BDN08], implementing AVMAX with that $m$ value to
learn endmembers [CMAC11], and calculating abundances and purity through a
nonnegative least squares solver [BDJ97]. D-VIC then estimates empirical
density using a KDE [RL14, MM18, MM19b].

To locate pixels that are both high-density and indicative of an underlying
material, D-VIC calculates

$\zeta(x)=\leavevmode\nobreak\
\frac{2\bar{p}(x)\bar{\eta}(x)}{\bar{p}(x)+\bar{\eta}(x)},$

where $\bar{p}(x)=\frac{p(x)}{\max_{y\in X}p(y)}$ and
$\bar{\eta}(x)=\leavevmode\nobreak\ \frac{\eta(x)}{\max_{y\in X}\eta(y)}$.
Thus, $\zeta(x)$ returns the harmonic mean of $p(x)$ and $\eta(x)$, which are
normalized so that density and purity are approximately at the same scale. By
construction, $\zeta(x)\approx 1$ only at high-density, high-purity pixels
$x$. In contrast, if a pixel $x$ is either low-density or low-purity, then
$\zeta(x)$ will be small. Importantly, $\zeta(x)$ downweights mixed pixels
that, though high-density, correspond to a spatial region containing many
materials. Thus, points with large $\zeta$-values will correspond to pixels
that are modal (due to their high empirical density) and representative of
just one material in the scene (due to their high pixel purity).

D-VIC uses the following function to incorporate diffusion geometry into its
procedure for selecting cluster modes:

$\displaystyle d_{t}(x)=\begin{cases}\max\limits_{y\in
X}D_{t}(x,y)&x=\operatorname*{argmax}\limits_{y\in X}\zeta(y),\\\
\min\limits_{y\in
X}\\{D_{t}(x,y)|\zeta(y)\geq\zeta(x)\\}&\text{otherwise}.\end{cases}$

Thus, a pixel will have a large $d_{t}$-value if it is far in diffusion
distance at time $t$ from its $D_{t}$-nearest neighbor of higher density and
pixel purity. D-VIC assigns modal labels to the $K$ points maximizing of
$\mathcal{D}_{t}(x)=d_{t}(x)\zeta(x)$, which are high-density, high-purity
pixels far in diffusion distance at time $t$ from other high-density, high-
purity pixels.

Figure 2:  Diagram of the D-VIC clustering algorithm. The computational
complexity of each step is colored in red. Note that all steps are quasilinear
with respect to $n$, implying that D-VIC scales well to large HSI datasets.

After labeling cluster modes, D-VIC labels non-modal points according to their
$D_{t}$-nearest neighbor of higher $\zeta(x)$-value that is already labeled.
Importantly, though this non-modal labeling scheme is similar to the one used
in LUND, D-VIC downweights low-purity pixels through $\zeta(x)$ in its non-
modal labeling. Thus, pixel purity is incorporated in all stages of the D-VIC
algorithm through $\zeta(x)$. D-VIC is provided in Algorithm 2, and a
schematic is provided in Figure 2.

###  3.1 Computational Complexity

This section considers the computational complexity of D-VIC. The
computational complexity of the HySime algorithm is $O(D^{2}n)$ operations
[BDN08], whereas the computational complexity of spectral unmixing using AVMAX
and a standard nonnegative least squares solver [BDJ97] is
$O((D^{2}+m^{4}+m^{2}I)n)$ operations, where $I$ is the number of AVMAX
partial maximizations [CMAC11]. We assume that nearest neighbor searches are
performed using _cover trees_ : an indexing data structure that enables
logarithmic nearest neighbor searches [BKL06]. To see this, define the
_doubling dimension_ of $X$ by $d=\log_{2}(c)$, where $c>0$ is the smallest
value for which any ball $B_{p}(p,r)=\\{q\in\leavevmode\nobreak\
X\leavevmode\nobreak\ |\leavevmode\nobreak\ \|p-\leavevmode\nobreak\
q\|\leavevmode\nobreak\ \leq\leavevmode\nobreak\ r\\}$ can be covered by $c$
balls of radius $r/2$. If the spectral signatures of pixels in
$X\subset\mathbb{R}^{D}$ have doubling dimension $d$, a search for the $N$
$\ell^{2}$-nearest neighbors of each HSI pixel using cover trees has
computational complexity $O(NDC^{d}n\log(n))$, where $C$ is a constant
independent of $n$, $D$, $N$, and $d$. Thus, if $\mathbf{W}$ is constructed
using cover trees [BKL06] with $N$ nearest neighbors, and $O(1)$ eigenvectors
of $\mathbf{P}$ are used to calculate diffusion distances, then the
computational complexity of D-VIC is
$O((D^{2}+m^{4}+m^{2}I)n+\leavevmode\nobreak\ NDC^{d}n\log(n))$ [MM19b,
BKL06].

So long as the spatial dimensions of the scene captured by an HSI are not
changed, we expect that $m$ (the expected number of materials within the
scene) will be constant with respect to the number of samples $n$. Similarly,
numerical simulations have shown that, if $m$ remains constant as the number
of samples increases, then $I$ tends to grow only slightly [CMAC11]. If
$m=O(1)$ and $I=O(\log(n))$ with respect to $n$, then the complexity of D-VIC
reduces to $O(NDC^{d}n\log(n))$.

###  3.2 Comparison to Related HSI Clustering Algorithms

In this section, D-VIC is compared theoretically against related graph-based
HSI clustering algorithms; empirical comparisons appear in Section 4.

####  3.2.1 Symmetric Nonnegative Matrix Factorization

Nonnegative matrix factorization (NMF)-based clustering algorithms typically
assume that each cluster in a latent clustering can be represented using just
a single basis vector [GKP14, CHNP04, KDP12]. Though this assumption holds for
datasets sampled from a union of affine spaces, it can fail on HSIs, which
often contain nonlinear structure [ZWL+16, KDP12]. Rather than directly
decomposing the pixels in $X$ as is done in other NMF-based clustering
algorithms [GKP14], the symmetric NMF (SymNMF) algorithm decomposes the weight
matrix $\mathbf{W}$ using the following optimization:
$\mathbf{H}^{*}=\leavevmode\nobreak\
\operatorname*{argmin}_{\mathbf{H}\in\mathbb{R}^{n\times
K}}\|\mathbf{W}-\leavevmode\nobreak\ \mathbf{H}\mathbf{H}^{T}\|_{F}^{2}$
subject to $\mathbf{H}_{ij}\geq 0$. Assuming data points from the same cluster
are weighted highly in the graph, but points from different clusters have low
edge weights, the matrix $\mathbf{H}^{*}$ will naturally encode latent
(possibly nonlinear) cluster structure in $X$. After performing its NMF step,
SymNMF labels points according to the largest entries in $\mathbf{H}^{*}$:
$\mathcal{C}(x_{i})=\operatorname*{argmax}_{1\leq j\leq
K}\mathbf{H}^{*}_{ij}$.

D-VIC and SymNMF are similar in that they both rely on the decomposition of a
matrix encoding a graph representation of $X$ to learn underlying (possibly
nonlinear) latent cluster structure. Nevertheless, there are key differences
that may lead SymNMF to erroneous clusterings. SymNMF does not directly
incorporate empirical density or pixel purity into its labeling, leaving it
prone to error on datasets containing multimodal cluster structure or HSIs
with mixed pixels. In contrast, D-VIC upweights high-density, high-purity
points within its mode selection, making it robust to error on these
problematic classes of data. Moreover, though SymNMF has been shown to perform
well on real HSI data [KDP12], its cluster assignments tend to be sensitive to
the initialization of $\mathbf{H}$.

####  3.2.2 KNN-SSC

The _$K$ -Nearest Neighbors-based Sparse Subspace Clustering_ (_KNN-SSC_)
algorithm modifies SC to learn the best locality-preserving low-rank
representation of the data [LLY+12, ZWL+16]. More precisely, KNN-SSC finds
$(\mathbf{Z}^{*},\mathbf{E}^{*})\in\mathbb{R}^{n\times
n}\times\mathbb{R}^{n\times D}$ minimizing
$\|\mathbf{Z}\|_{*}+\lambda\|\mathbf{E}\|_{2,1}$ ($\lambda>0$ is a
regularization parameter) subject to the following constraints: (i)
$\mathbf{X}=\mathbf{Z}\mathbf{X}+\mathbf{E}$ so that data points are expressed
as linear combinations of other data points; (ii)
$\mathbf{1}\mathbf{Z}=\mathbf{1}$ so that coefficients in these linear
combinations sum to 1; and (iii) $\mathbf{Z}_{ij}>0$ if and only if $x_{j}$ is
one of the $N$ $\ell^{2}$-nearest neighbors of $x_{i}$. SC is then implemented
using the weight matrix
$\mathbf{W}=\frac{1}{2}(|\mathbf{Z}^{*}|+|\mathbf{Z}^{*}|^{T})$ to produce the
KNN-SSC clustering of $X$ [NJW02, ZWL+16].

The KNN-SSC and D-VIC algorithms are similar in that they both rely on a
linear unmixing step to extract additional information from the dataset
[LLY+12, ZWL+16, CP21]. However, KNN-SSC does not incorporate data density
into its clustering procedure, leaving it prone to error on data with
multimodal cluster structure [NG07, MM19b, MP22]. Moreover, while KNN-SSC
relies on the minimization of a convex objective function subject to linear
constraints, there are at least two main difficulties in the KNN-SSC
optimization [ZWL+16]. First, though the objective function
$\|\mathbf{Z}\|_{*}+\lambda\|\mathbf{E}\|_{2,1}$ is convex, it is not smooth.
Second, the constraint that $\mathbf{Z}_{ij}>0$ if and only if $x_{i}$ and
$x_{j}$ are $\ell^{2}$-nearest neighbors significantly restricts the domain.
These complications may yield a much slower convergence, even if the graph
provided is superior. A final concern is that KNN-SSC scales quadratically
with $n$ in storage costs. As such, KNN-SSC may not be suitable for clustering
large HSI datasets.

####  3.2.3 Learning by Unsupervised Nonlinear Diffusion

We remark that, for any one cluster, there may be multiple reasonable choices
for cluster modes: pixels that are exemplary of underlying cluster structure.
In LUND, cluster modes are selected to be high-density pixels that are far in
diffusion distance from other high-density pixels [MM19b, MM18]. However, not
all high-density pixels necessarily correspond to underlying material
structure. A maximizer of $p(x)$ could, for example, correspond to a spatial
region containing a group of commonly co-occurring materials (rather than a
single material). By weighting pixel purity and density equally, D-VIC avoids
selecting such a pixel as cluster mode; thus, D-VIC modes will be both
indicative of underlying material structure and modal, making these points
better exemplars of underlying material structure than the modes selected by
LUND.

Figure 3: Ground truth labels and pixel purity of synthetic dataset. Notice
that empirical density maximizers near the origin are also lowest-purity data
points.

This section illustrates the key difference between LUND and D-VIC by
implementing both algorithms on a simple dataset (visualized in Figure 3)
built to illustrate the idealized scenario where D-VIC outperforms LUND due to
its incorporation of pixel purity. This dataset was generated by sampling
$n=5000$ points from an equilateral triangle in $\mathbb{R}^{2}$ centered at
the origin with edge length $\sqrt{2}$; the $K=3$ vertices of this triangle
served as ground truth endmembers. We sampled 1000 data points from a Gaussian
distribution with a standard deviation of 0.175 centered at each endmember,
keeping only the samples lying within in the convex hull of the ground truth
endmember set. In addition, 2000 data points were sampled from a Gaussian
distribution with zero-mean and a smaller standard deviation of 0.0175. As
such, high-purity points indicative of latent material structure were also
relatively low-density, and density maximizers were engineered so to not be
indicative of latent material structure. Each point was assigned a ground
truth label corresponding to its highest-abundances ground truth endmember.

For both LUND and D-VIC, overall accuracy (OA), defined to be the fraction of
correctly labeled pixels, was optimized for across the same grid of relevant
hyperparameter values (see Appendix A). The optimal clusterings and their
corresponding OA values are provided in Figure 4. These numerical results
illustrate a fundamental limitation of relying solely on empirical density to
select cluster modes in spectrally mixed HSI data. Because empirical density
maximizers are not representative of underlying material structure in this
synthetic dataset, LUND is unable to discriminate data points within the high-
density, low-purity region near the origin, resulting in poor performance and
an OA of 0.739. In contrast, D-VIC downweights high-density points that are
not also high-purity and therefore selects points that are more representative
of the dataset’s underlying material structure as cluster modes. As a result,
D-VIC correctly separates the high-density, low-purity region into three
segments, yielding a substantially higher OA of 0.905, a difference of 0.166
when compared to LUND.

Figure 4: Optimal LUND (OA = 0.739) and D-VIC (OA = 0.905) clusterings of the
synthetic dataset (Figure 3).

##  4 Numerical Experiments

This section contains numerical results indicating the efficacy of D-VIC.
Classical and state-of-the-art clustering algorithms were implemented on three
real HSI datasets described in Section 4.1. Classical algorithms D-VIC was
compared against include $K$-Means [FHT01], $K$-Means applied to the first
principal components of the HSI ($K$-Means+PCA) [FHT01], and GMM applied to
first principal components of the HSI (GMM+PCA) [FHT01]. D-VIC was also
compared against several state-of-the-art HSI clustering algorithms: SC
[NJW02], SymNMF [KDP12], KNN-SSC [AG21, ZWL+16], and LUND [MM19b, MM18] (see
Sections 2.2.1, 2.3, and 3.2).

In all experiments, the number of clusters was set equal to the ground truth
$K$. We clustered the entire image but discarded unlabeled pixels when
comparing clusterings to the ground truth labels. Comparisons were made using
OA and Cohen’s $\kappa$ coefficient: $\kappa=\frac{\text{OA}-p_{e}}{1-p_{e}}$,
where $p_{e}$ is the probability that a clustering agrees with the ground
truth labels by chance. OA was optimized for across hyperparameters ranging a
grid of relevant values for each algorithm (see Appendix A). We report the
median OA across 100 trials for K-Means, GMM, SymNMF, and D-VIC to account for
the stochasticity associated with random initial conditions. Diffusion
distances were computed using only the first 10 eigenvectors of $\mathbf{P}$
in LUND and D-VIC. For D-VIC, AVMAX was implemented 100 times in parallel, and
the endmember set that formed the largest-volume simplex was selected for
later cluster analysis. Numerical experiments can be replicated at
https://github.com/sampolk/D-VIC.

###  4.1 Description of Benchmark HSI Datasets

To illustrate the efficacy of D-VIC, we analyzed three publicly available,
real HSIs often used as benchmarks for new HSI clustering algorithms; see
Figure 5 and the descriptions below. Water absorption bands were discarded,
and pixel spectra were standardized before analysis.

(a)

(b)

(c)

Figure 5: Ground truth labels and first principal component scores for the
real benchmark HSIs analyzed in this article.

  1. 1.

Salinas A (Figure 5(a)) was recorded by the Airborne Visible/Infrared Imaging
Spectrometer (AVIRIS) sensor over farmland in Salinas Valley, California, USA
in 1998 at a spatial resolution of 1.3 m. Spectral signatures, ranging in
recorded wavelength from 380 nm to 2500 nm across 224 spectral bands, were
recorded across $83\times 86$ pixels ($n=7138$). Gaussian noise (with mean 0
and standard deviation $=10^{-7}$) was added to each pixel to differentiate
two pixels with identical spectral signatures. The Salinas A scene contains
$K=6$ ground truth classes corresponding to crop types.

  2. 2.

Jasper Ridge (Figure 5(b)) was recorded by the AVIRIS sensor over the Jasper
Ridge Biological Preserve, California, USA in 1989 at a spatial resolution of
5 m. Spectral signatures, ranging in recorded wavelength from 380 nm to 2500
nm across 224 spectral bands, were recorded across spatial dimensions of
$100\times 100$ pixels ($n=10000$). The Jasper Ridge scene contains $K=4$
ground truth endmembers: road, soil, water, and trees. Ground truth labels
were recovered by selecting the material of highest ground truth abundance for
each pixel.

  3. 3.

Indian Pines (Figure 5(c)) was recorded by the AVIRIS sensor over farmland in
northwest Indiana, USA in 1992 at a low spatial resolution of 20 m. Spectral
signatures, ranging in recorded wavelength from 400 nm to 2500 nm across 224
spectral bands, were recorded across spatial dimensions of $145\times 145$
pixels ($n=21025$). The Indian Pines scene contains $K=16$ ground truth
classes (e.g., crop types and manufactured structures) as well as many
unlabeled pixels.

###  4.2 Performance Comparisons

This section compares clusterings produced by D-VIC against those of related
algorithms (Table 1). On each of the three benchmark HSIs analyzed, D-VIC
produces a clustering that is closer to the ground truth labels than those
produced by related algorithms. The Indian Pines dataset (Fig. 5(c)) is
expected to be challenging because pixels from the same class exist in
multiple segments of the image, and the size of ground truth clusters varies
substantially across the $K=16$ classes. Nevertheless, D-VIC achieves higher
performance than all other algorithms on this challenging dataset. Notably,
though all other algorithms (including state-of-the-art algorithms, such as
LUND) achieve $\kappa$-statistics in the same narrow range of 0.298 to 0.316,
D-VIC achieves a substantially higher $\kappa$-statistic of 0.350. As such,
the incorporation of pixel purity in D-VIC enables superior detection even in
this difficult setting.

Figure 6: Comparison of clusterings produced by D-VIC and related algorithms on the Salinas A HSI. Unlike any comparison method, D-VIS correctly groups pixels corresponding to 8-week maturity romaine (indicated in yellow), resulting in near-perfect recovery of the ground truth labels.  | Salinas A | Jasper Ridge | Indian Pines
---|---|---|---
| OA | $\kappa$ | OA | $\kappa$ | OA | $\kappa$
$K$-Means | 0.764 | 0.703 | 0.784 | 0.703 | 0.383 | 0.315
$K$-Means+PCA | 0.764 | 0.703 | 0.785 | 0.703 | 0.382 | 0.316
GMM+PCA | 0.611 | 0.512 | 0.789 | 0.701 | 0.364 | 0.292
SC | 0.834 | 0.797 | 0.760 | 0.670 | 0.382 | 0.314
SymNMF | 0.828 | 0.791 | 0.662 | 0.542 | 0.365 | 0.304
KNN-SSC | 0.844 | 0.809 | 0.726 | 0.629 | 0.371 | 0.308
LUND | 0.887 | 0.860 | 0.815 | 0.737 | 0.404 | 0.312
D-VIC | 0.976 | 0.970 | 0.865 | 0.805 | 0.445 | 0.350
Table 1:  Performances of D-VIC and related algorithms on benchmark HSIs.
Highest and second-highest performances are bolded and italicized
respectively. The median OA across 100 trials is reported for $K$-Means,
$K$-Means+PCA, GMM+PCA, SC, SymNMF, and D-VIC. D-VIC offers substantially
higher performance on all datasets evaluated.

As visualized in Figure 6, D-VIC achieved nearly perfect recovery of the
ground truth labels for Salinas A. Most notably, though all comparison methods
erroneously separate the ground truth cluster indicated in yellow in Figure
5(a) (corresponding to 8-week maturity romaine), D-VIC correctly groups the
pixels in this cluster, resulting in a performance that was 0.089 higher in OA
and 0.110 in $\kappa$ than the that of LUND, its closest competitor in Table
1. As such, downweighting high-density points that are not also exemplary of
the latent material structure improves not only modal, but also non-modal
labeling. Moreover, what error does exist in the D-VIC clustering of Salinas A
could likely be remedied through spatial regularization or post-processing
smoothing [MM19c, Mur20].

D-VIC similarly achieved much higher performance than related state-of-the-art
graph-based clustering algorithms on Jasper Ridge (as visualized in Figure 7).
This difference in performance was substantially driven by superior separation
of the classes indicated in dark blue (corresponding to tree cover) and green
(corresponding to soil) in Figure 5(b). Indeed, though LUND groups most tree
cover pixels with soil pixels in Figure 7, D-VIC correctly separates much of
the latent structure for this class. The difference between LUND’S and D-VIC’s
clusterings indicates that the pixels corresponding to the tree cover class,
though lower density than pixels corresponding to the soil class, have
relatively high pixel purity.

f

Figure 7: Comparison of clusterings produced by D-VIC and related algorithms
on the Jasper Ridge HSI. D-VIC outperforms all other algorithms, largely due
to superior discrimination of the pixels corresponding to tree and soil
classes (indicated in dark blue and green respectively).

Figure 8: Visualization of D-VIC’s median OA across 50 trials as
hyperparameters $N$ and $\sigma_{0}$ are varied. D-VIC achieves high
performance across a large set of hyperparameters.

###  4.3 Robustness to Hyperparameter Selection

This section analyzes the robustness of D-VIC’s performance to hyperparameter
selection. For each node in a grid of $(N,\sigma_{0})$, D-VIC was implemented
50 times, and the median OA value across these 50 trials was stored.
Performance degraded as $N$ increased substantially past 100, and such a
choice is not advised. In Figure 8, we visualize how the performance of D-VIC
varies with $N$ and $\sigma_{0}$. The relative nominal values in the range of
$\sigma_{0}$ reflect that pixels from the HSIs analyzed in this article are
relatively close to their $\ell^{2}$-nearest neighbors on average. As is
described in Appendix A, the range of $\sigma_{0}$ used for each grid search
is data-dependent, ranging the distribution of $\ell^{2}$-distances between
pixel spectra and their 1000 $\ell^{2}$-nearest neighbors. Thus, the ranges
indicated in Figure 8 cover the relevant values of $\sigma_{0}$.

Figure 9: Analysis of D-VIC’s performance as $t$ varies across $[10,200]$. Values are the median OA across 100 implementations of D-VIC with the optimal $N$ and $\sigma_{0}$ values. Generally, $t$ appears to have little impact on the OA of D-VIC, and D-VIC achieves performances comparable to those reported in Table 1 across $t\in[90,200]$, uniformly across all data sets.  | Salinas A | Jasper Ridge | Indian Pines
---|---|---|---
$K$-Means | 0.04 | 0.10 | 1.04
$K$-Means+PCA | 0.10 | 0.14 | 0.58
GMM+PCA | 0.13 | 0.23 | 2.19
SC | 1.82 | 3.15 | 14.54
SymNMF | 3.50 | 4.42 | 48.29
KNN-SSC | 4.11 | 7.91 | 103.05
LUND | 2.35 | 4.14 | 14.74
D-VIC | 4.95 | 7.64 | 23.70
Table 2: Runtimes (seconds) of D-VIC and related algorithms. D-VIC achieves
runtimes comparable to state-of-the-art algorithms and scales well to the
larger Indian Pines dataset.

It is clear from Figure 8 that D-VIC is capable of achieving high performance
across a broad range of hyperparameters on each HSI. Thus, given little
hyperparameter tuning, D-VIC is likely to output a partition that is
competitive with clusterings reported in Table 1. Figure 8 also motivates
recommendations for hyperparameter selection to optimize the OA of D-VIC.
Larger datasets (e.g., Indian Pines) tend to require larger values of $N$ for
D-VIC to achieve high OA, corresponding with recommendations in the literature
that $N$ should grow logarithmically with $n$ [MM19b]. Additionally, for
datasets with high-purity material classes (e.g., Salinas A), D-VIC achieves
highest OA using large $\sigma_{0}$. This reflects that, as $\sigma_{0}$
increases, the KDE $p(x)$ becomes more constant across the HSI and
$\zeta(x)\approx\eta(x)$. Since purity is an excellent indicator of material
class structure for Salinas A, D-VIC becomes better able to recover the latent
material structure with larger $\sigma_{0}$.

We also analyze the robustness of D-VIC’s performance to the selection of the
diffusion time parameter $t$. Using the optimal values of $N$ and
$\sigma_{0}$, D-VIC was evaluated 100 times at $t$-values ranging
$\\{10,20,\dots,200\\}$. Figure 9, which visualizes the results of this
analysis, indicates D-VIC achieves high OA values across a broad range of $t$;
for each $t\in[90,200]$, D-VIC outputs a clustering with OA equal to or very
close to those reported in Table 1. These results indicate that D-VIC is well-
equipped to provide high-quality clusterings given little or no tuning of $t$.
Indeed, a simple choice of $t=100$ works extremely well across all datasets.

###  4.4 Runtime Analysis

This section compares runtimes of the algorithms implemented in Section 4.2,
where hyperparameters were set to be those which produced the results in Table
1. All experiments were run in MATLAB on the same environment: a macOS Big Sur
system with an 8-core Apple® M1™ Processor and 8 GB of RAM. Each core had a
processor base frequency of 3.20 GHz. Runtimes are provided in Table 2. All
classical algorithms have smaller runtimes than D-VIC, but the performances
reported in Table 1 for these algorithms are substantially less than those
reported for D-VIC. On the other hand, though KNN-SSC and SymNMF achieve
performances competitive to D-VIC, unlike D-VIC, these algorithms appear to
scale poorly to large datasets. Finally, D-VIC outperforms LUND at the cost of
only a small increase in runtime (associated with an additional spectral
unmixing step).

##  5 Conclusion

This article introduces the _D-VIC clustering algorithm_ for unsupervised
material discrimination in HSIs. D-VIC assigns modal labels to high-density,
high-purity pixels within the HSI that are far in diffusion distance from
other high-density, high-purity pixels [CMAC11, MM19b, CP21]. We have argued
that these _cluster modes_ are highly indicative of underlying material
structure, leading to more interpretable and accurate clusterings than those
produced by related algorithms [MM18, MM19b]. Indeed, numerical experiments
presented in Section 4 show that incorporating pixel purity into D-VIC results
in clusterings closer to the ground truth labels on three benchmark real HSI
datasets of varying sizes and complexities. As such, D-VIC is equipped to
perform efficient material discrimination on broad ranges of mixed HSI
datasets.

Future work includes modifying the spectral unmixing step in D-VIC to account
for nonlinear mixing of endmembers. Linear endmember extraction is
computationally inexpensive [BDN08, CMAC11, BDJ97], and it results in strong
performance in D-VIC, but recent years have brought significant advances in
algorithms for the nonlinear spectral unmixing of HSIs [CRH12, HADT11, YCI13].
Modifying the spectral unmixing step in D-VIC may improve performance,
especially for HSIs in which assumptions on linear mixing do not hold [CRH12,
HADT11, YCI13].

Additionally, much of the error in D-VIC’s clusterings may be corrected
through by incorporating spatial information into its labeling. Such a
modification of D-VIC may improve performance on datasets with spatially
homogeneous clusters [MM19c, PM21, FTB+12, GBU13, FLD+15, TBC09, MSB07].
Moreover, it is likely that varying the diffusion time parameter $t$ in D-VIC
may enable the detection of multiple scales of latent cluster structure: a
problem we would like to consider further in future work [MP22, PM21, GKP14].
Finally, we expect D-VIC may be modified for _active learning_ , wherein
ground truth labels for a small number of carefully selected pixels are
queried and propagated across the image [MM19a, Mur20].

## Acknowledgment

The US National Science Foundation partially supported this research through
grants NSF-DMS 1912737, NSF-DMS 1924513, and NSF-CCF 1934553. We thank C.
Schönlieb, M. S. Kotzagiannidis, and D. Coomes for conversations that aided in
the development of D-VIC. In addition, we thank N. Gillis, D. Kuang, H. Park,
C. Ding, M. Abdolali, D. Kun, and I. Gerg for making code for HySime [BDN08],
AVMAX [GK22, CMAC11], SymNMF [KDP12], and KNN-SSC [AG21] publicly available.

##  Appendix A Hyperparameter Optimization

This section describes the hyperparameter optimization performed to generate
Table 1. The parameter grids used for each algorithm are summarized in Table
3. For $K$-Means+PCA and GMM+PCA, we clustered the first $z$ principal
components of the HSI, where $z$ was chosen so that 99% of the variation in
the HSI was maintained after PCA dimensionality reduction. Thus, $K$-Means,
$K$-Means+PCA, and GMM+PCA required no hyperparameter inputs. For stochastic
algorithms with hyperparameter inputs (SC, SymNMF, and D-VIC), we optimized
for the median OA across 10 trials at each node in the hyperparameter grids
described below.

Algorithm | Param. 1 Grid | Param. 2 Grid | Param. 3 Grid
---|---|---|---
$K$-Means | — | — | —
$K$-Means+PCA | — | — | —
GMM+PCA | — | — | —
SC [NJW02] | $N\in\mathcal{N}$ | — | —
SymNMF [KDP12] | $N\in\mathcal{N}$ | — | —
KNN-SSC [AG21, ZWL+16] | $N\in\mathcal{N}$ | $\lambda=10$ | —
LUND [MM19b] | $N\in\mathcal{N}$ | $\sigma_{0}\in\mathscr{D}$ | $t\in\mathscr{T}$
D-VIC | $N\in\mathcal{N}$ | $\sigma_{0}\in\mathscr{D}$ | $t\in\mathscr{T}$
Table 3: Hyperparameter grids for algorithms. The number of nearest neighbors
$N$ took values in $\mathcal{N}$: an exponential sampling of the set
$[10,900]$. The set $\mathscr{D}$ contains $\ell^{2}$-distances between data
points and their 1000 $\ell^{2}$-nearest neighbors. The set $\mathscr{T}$ is
an exponential sampling of the diffusion process:
$\mathscr{T}=\\{0,1,2,\dots,2^{2},\dots,2^{T}\\}$. A “—” indicates a lack of a
hyperparameter input.

All graph-based algorithms relied on adjacency matrices built from sparse KNN
graphs. The number of nearest neighbors was optimized for each algorithm
across $\mathcal{N}$: an exponential sampling of the set $[10,900]$. KNN-SSC’s
regularization parameter was set to $\lambda=10$, motivated by prior work with
this parameter [AG21]. In addition to $N$, LUND and D-VIC clustering
algorithms also require the KDE scale $\sigma_{0}>0$ and diffusion time $t$ as
input. In our grid searches, $\sigma_{0}$ ranged $\mathscr{D}$: a sampling of
the distribution of $\ell^{2}$-distances between data points and their 1000
$\ell^{2}$-nearest neighbors. Similarly, both LUND and D-VIC were implemented
at time steps $t\in\mathscr{T}=\\{0,1,2,2^{2},\dots,2^{T}\\}$, where
$T=\lceil\log_{2}[\log_{\lambda_{2}(\mathbf{P})}(\frac{2\times
10^{-5}}{\min(\pi)})]\rceil$. Searches end at time $t=2^{T}$ because
$\max_{x,y\in X}D_{t}(x,y)\leq 10^{-5}$ for $t\geq\leavevmode\nobreak\ 2^{T}$.
For each dataset, we chose the $t\in\\{0,1,2,2^{2},\dots,2^{T}\\}$ resulting
in maximal OA. We selected $t=2^{5}$ for Salinas A, $t=2^{7}$ for Jasper
Ridge, and $t=0$ for Indian Pines. As noted in Section 4.3, D-VIC is quite
robust to this choice of parameter.

## References

  * [AG21] Maryam Abdolali and Nicolas Gillis.  Beyond linear subspace clustering: A comparative study of nonlinear manifold clustering algorithms.  Comput Sci Rev, 42:1574–0137, 2021.
  * [AMR10] E Adam, O Mutanga, and D Rugege.  Multispectral and hyperspectral remote sensing for identification and mapping of wetland vegetation: a review.  Wetl Ecol Manag, 18(3):281–296, 2010.
  * [BAF05] Charles M Bachmann, Thomas L Ainsworth, and Robert A Fusina.  Exploiting manifold geometry in hyperspectral imagery.  IEEE Trans Geosci Remote Sens, 43(3):441–454, 2005.
  * [BDJ97] Rasmus Bro and Sijmen De Jong.  A fast non-negativity-constrained least squares algorithm.  J Chemom, 11(5):393–401, 1997.
  * [BDN08] José M Bioucas-Dias and José MP Nascimento.  Hyperspectral subspace identification.  IEEE Trans Geosci Remote Sens, 46(8):2435–2445, 2008.
  * [BDPCV+13] José M Bioucas-Dias, Antonio Plaza, Gustavo Camps-Valls, Paul Scheunders, Nasser Nasrabadi, and Jocelyn Chanussot.  Hyperspectral remote sensing data analysis and future challenges.  IEEE Geosci Remote Sens Mag, 1(2):6–36, 2013.
  * [BKG95] J W Boardman, F A Kruse, and R O Green.  Mapping target signatures via partial unmixing of AVIRIS data.  1995\.
  * [BKL+04] M Berman, H Kiiveri, R Lagerstrom, A Ernst, R Dunne, and J F Huntington.  ICE: A statistical approach to identifying endmembers in hyperspectral images.  IEEE Trans Signal Process, 42(10):2085–2095, 2004.
  * [BKL06] Alina Beygelzimer, Sham Kakade, and John Langford.  Cover trees for nearest neighbor.  In Proc ICML, pages 97–104, 2006.
  * [BN01] M Belkin and Partha Niyogi.  Laplacian eigenmaps and spectral techniques for embedding and clustering.  In Proc NeurIPS, volume 14, 2001.
  * [Boa93] Joseph W Boardman.  Automating spectral unmixing of AVIRIS data using convex geometry concepts.  In Proc Annu JPL Airborne Geosci Workshop, volume 1, pages 11–14, 1993.
  * [CCHM09] T Chan, C Chi, Y Huang, and W Ma.  A convex analysis-based minimum-volume enclosing simplex algorithm for hyperspectral unmixing.  IEEE Trans Signal Process, 57(11):4418–4432, 2009.
  * [CCM14] N D Cahill, W Czaja, and D W Messinger.  Schroedinger eigenmaps with nondiagonal potentials for spatial-spectral clustering of hyperspectral imagery.  In Proc SPIE XX, volume 9088, pages 27–39. SPIE, 2014.
  * [CF18] W Chen and A L Ferguson.  Molecular enhanced sampling with autoencoders: On-the-fly collective variable discovery and accelerated free energy landscape exploration.  J of Comput Chem, 39(25):2079–2102, 2018.
  * [Cha18] Chein Chang.  A review of virtual dimensionality for hyperspectral imagery.  IEEE J Sel Top Appl Earth Obs Remote Sens, 11(4):1285–1305, 2018\.
  * [CHNP04] M Catral, Lixing Han, Michael Neumann, and Robert J Plemmons.  On reduced rank nonnegative matrix factorization for symmetric nonnegative matrices.  Linear Algebra Appl, 393:107–126, 2004.
  * [CKS10] J G P W Clevers, L Kooistra, and M E Schaepman.  Estimating canopy water content using hyperspectral remote sensing data.  Int J Appl Earth Obs Geoinf, 12(2):119–125, 2010.
  * [CL06] Ronald R Coifman and Stéphane Lafon.  Diffusion maps.  Appl and Comput Harm Anal, 21(1):5–30, 2006.
  * [CLL+05] Ronald R Coifman, Stephane Lafon, Ann B Lee, Mauro Maggioni, Boaz Nadler, Frederick Warner, and Steven W Zucker.  Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps.  Proc Natl Acad Sci USA, 102(21):7426–7431, 2005.
  * [CMAC11] Tsung Chan, Wing Ma, ArulMurugan Ambikapathi, and Chong Chi.  A simplex volume maximization framework for hyperspectral endmember extraction.  IEEE Trans Geosci Remote Sens, 49(11):4177–4193, 2011.
  * [CMR13] Daniele Cerra, Rupert Müller, and Peter Reinartz.  Noise reduction in hyperspectral images through spectral unmixing.  IEEE Geosci Remote Sens Lett, 11(1):109–113, 2013.
  * [COR11] Adam S Charles, Bruno A Olshausen, and Christopher J Rozell.  Learning sparse codes for hyperspectral imagery.  IEEE Journal of Selected Topics in Signal Processing, 5(5):963–978, 2011.
  * [CP21] Kangning Cui and Robert J Plemmons.  Unsupervised classification of AVIRIS-NG hyperspectral images.  In Proc WHISPERS, pages 1–5. IEEE, 2021.
  * [CRH12] J Chen, C Richard, and Paul Honeine.  Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model.  IEEE Trans Signal Process, 61(2):480–492, 2012.
  * [CVMZ07] G Camps-Valls, T V B Marsheva, and D Zhou.  Semi-supervised graph-based hyperspectral image classification.  IEEE Trans Geosci Remote Sens, 45(10):3044–3054, 2007.
  * [CWLO06] C Chang, C Wu, W Liu, and Y Ouyang.  A new growing method for simplex-based endmember extraction algorithm.  IEEE Trans Geosci Remote Sens, 44(10):2804–2819, 2006.
  * [DBG08] M Dalponte, L Bruzzone, and D Gianelle.  Fusion of hyperspectral and lidar remote sensing data for classification of complex forest areas.  IEEE Trans Geosci Remote Sens, 46(5):1416–1427, 2008.
  * [DMC+09] N Dobigeon, S Moussaoui, M Coulon, J Tourneret, and A O Hero.  Joint Bayesian endmember extraction and linear unmixing for hyperspectral imagery.  IEEE Trans Signal Process, 57(11):4355–4368, 2009.
  * [EGVL+12] Gerda J Edelman, Edurne Gaston, Ton G Van Leeuwen, PJ Cullen, and Maurice CG Aalders.  Hyperspectral imaging for non-contact analysis of forensic traces.  Forensic Sci Int, 223(1-3):28–39, 2012.
  * [Eis12] Michael Theodore Eismann.  Hyperspectral remote sensing.  SPIE, 2012.
  * [FHT01] Jerome Friedman, Trevor Hastie, and Robert Tibshirani.  The elements of statistical learning, volume 1.  Springer Series in Statistics, 2001.
  * [FLD+15] Leyuan Fang, Shutao Li, Wuhui Duan, Jinchang Ren, and Jón Atli Benediktsson.  Classification of hyperspectral images by exploiting spectral–spatial information of superpixel via multiple kernels.  IEEE Trans Geosci Remote Sens, 53(12):6663–6674, 2015.
  * [FTB+12] Mathieu Fauvel, Yuliya Tarabalka, Jon Atli Benediktsson, Jocelyn Chanussot, and James C Tilton.  Advances in spectral-spatial classification of hyperspectral images.  Proceedings of the IEEE, 101(3):652–675, 2012.
  * [GAC01] K J Guilfoyle, M L Althouse, and C Chang.  A quantitative and comparative analysis of linear and nonlinear spectral mixture models using radial basis function neural networks.  IEEE Trans Geosci Remote Sens, 39(10):2314–2318, 2001.
  * [GBU13] Pedram Ghamisi, Jon Atli Benediktsson, and Magnus Orn Ulfarsson.  Spectral–spatial classification of hyperspectral images based on hidden markov random fields.  IEEE Trans Geosci Remote Sens, 52(5):2565–2574, 2013.
  * [GJC+14] Y Gao, R Ji, P Cui, Q Dai, and G Hua.  Hyperspectral image classification through bilayer graph-based learning.  IEEE Trans Image Process, 23(7):2769–2778, 2014.
  * [GK22] Isaac Gerg and David Kun.  Hyperspectral toolbox.  https://github.com/davidkun/HyperSpectralToolbox, 2022.
  * [GKP14] Nicolas Gillis, Da Kuang, and Haesun Park.  Hierarchical clustering of hyperspectral images using rank-two nonnegative matrix factorization.  IEEE Trans Geosci Remote Sens, 53(4):2066–2078, 2014.
  * [GPC+17] Pedram Ghamisi, Javier Plaza, Yushi Chen, Jun Li, and Antonio J Plaza.  Advanced spectral classifiers for hyperspectral images: A review.  IEEE Geosci Remote Sens Mag, 5(1):8–32, 2017.
  * [HADT11] Abderrahim Halimi, Yoann Altmann, Nicolas Dobigeon, and Jean-Yves Tourneret.  Nonlinear unmixing of hyperspectral images using a generalized bilinear model.  IEEE Trans Geosci Remote Sens, 49(11):4153–4162, 2011.
  * [HBW+16] L Haghverdi, M Büttner, F A Wolf, F Buettner, and F J Theis.  Diffusion pseudotime robustly reconstructs lineage branching.  Nat methods, 13(10):845–848, 2016.
  * [HMW03] A Hirano, M Madden, and R Welch.  Hyperspectral image data for mapping wetland vegetation.  Wetl, 23(2):436–448, 2003.
  * [HYV20] Benjamin David Haeffele, Chong You, and Rene Vidal.  A critique of self-expressive deep subspace clustering.  In International Conference on Learning Representations, 2020.
  * [JWC+20] Jianxin Jia, Yueming Wang, Jinsong Chen, Ran Guo, Rong Shu, and Jianyu Wang.  Status and application of advanced airborne hyperspectral imaging technology: A review.  Infr Phys Technol, 104:103115, 2020.
  * [JWZ+17] J Jia, Y Wang, X Zhuang, Y Yao, S Wang, D Zhao, R Shu, and J Wang.  High spatial resolution shortwave infrared imaging technology based on time delay and digital accumulation method.  Infr Phys Technol, 81:305–312, 2017.
  * [KDP12] Da Kuang, Chris Ding, and Haesun Park.  Symmetric nonnegative matrix factorization for graph clustering.  In Proc SDM12, pages 106–117. SIAM, 2012.
  * [LDF11] G A Licciardi and F Del Frate.  Pixel unmixing in hyperspectral data by means of neural networks.  IEEE Trans Geosci Remote Sens, 49(11):4163–4172, 2011.
  * [LLY+12] Guangcan Liu, Zhouchen Lin, Shuicheng Yan, Ju Sun, Yong Yu, and Yi Ma.  Robust recovery of subspace structures by low-rank representation.  IEEE Trans Pattern Anal Mach Intell, 35(1):171–184, 2012.
  * [LQL+21] K Li, Y Qin, Q Ling, Y Wang, Z Lin, and W An.  Self-supervised deep subspace clustering for hyperspectral images with adaptive self-expressive coefficient matrix initialization.  IEEE J Sel Top Appl Earth Obs Remote Sens, 14:3215–3227, 2021.
  * [LSR15] Valero Laparrcr and Raúl Santos-Rodriguez.  Spatial/spectral information trade-off in hyperspectral images.  In Proc IGARSS, pages 1124–1127. IEEE, 2015.
  * [MBMDC06] S Moussaoui, D Brie, A Mohammad-Djafari, and C Carteret.  Separation of non-negative mixture of non-negative sources using a bayesian approach and MCMC sampling.  IEEE Trans Signal Process, 54(11):4133–4145, 2006.
  * [MM18] James M Murphy and Mauro Maggioni.  Unsupervised clustering and active learning of hyperspectral images with nonlinear diffusion.  IEEE Trans Geosci Remote Sens, 57(3):1829–1845, 2018.
  * [MM19a] Mauro Maggioni and James M Murphy.  Learning by active nonlinear diffusion.  Found Data Sci, 1(3):271, 2019.
  * [MM19b] Mauro Maggioni and James M Murphy.  Learning by unsupervised nonlinear diffusion.  J Mach Learn Res, 20(160):1–56, 2019.
  * [MM19c] James M Murphy and Mauro Maggioni.  Spectral–spatial diffusion geometry for hyperspectral image clustering.  IEEE Geosci Remote Sens Lett, 17(7):1243–1247, 2019.
  * [MP22] James M. Murphy and Sam L. Polk.  A multiscale environment for learning by diffusion.  Appl Comput Harm Anal, 57:58–100, 2022.
  * [MSB07] A Mohan, G Sapiro, and E Bosch.  Spatially coherent nonlinear dimensionality reduction and segmentation of hyperspectral images.  IEEE Geosci Remote Sens Lett, 4(2):206–210, 2007.
  * [MSS01] Dimitris Manolakis, Christina Siracusa, and Gary Shaw.  Hyperspectral subpixel target detection using the linear mixing model.  IEEE Trans Geosci Remote Sens, 39(7):1392–1409, 2001.
  * [Mur20] James M Murphy.  Spatially regularized active diffusion learning for high-dimensional images.  Pattern Recognit Lett, 135:213–220, 2020.
  * [MvDW+19] K R Moon, D van Dijk, Z Wang, S Gigante, D B Burkhardt, W S Chen, K Yim, A van den Elzen, M J Hirn, and R R Coifman.  Visualizing structure and transitions in high-dimensional biological data.  Nat Biotechnol, 37(12):1482–1492, 2019.
  * [MYW18] Dandan Ma, Yuan Yuan, and Qi Wang.  Hyperspectral anomaly detection via discriminative feature learning with multiple-dictionary sparse representation.  Remote Sens, 10(5):745, 2018.
  * [ND05] José MP Nascimento and José MB Dias.  Vertex component analysis: A fast algorithm to unmix hyperspectral data.  IEEE Trans Geosci Remote Sens, 43(4):898–910, 2005.
  * [Nev99] R Neville.  Automatic endmember extraction from hyperspectral data for mineral exploration.  In Proc CSRS, 1999.
  * [NG07] Boaz Nadler and Meirav Galun.  Fundamental limitations of spectral clustering.  In Proc NeurIPS, pages 1017–1024, 2007.
  * [NJW02] Andrew Y Ng, Michael I Jordan, and Yair Weiss.  On spectral clustering: Analysis and an algorithm.  In Proc NeurIPS, pages 849–856, 2002.
  * [NLCK06] Boaz Nadler, Stéphane Lafon, Ronald R Coifman, and Ioannis G Kevrekidis.  Diffusion maps, spectral clustering and reaction coordinates of dynamical systems.  Appl Comput Harmon Anal, 21(1):113–127, 2006.
  * [NMI+20] Jakub Nalepa, Michal Myller, Yasuteru Imai, Ken-ichi Honda, Tomomi Takeda, and Marek Antoniak.  Unsupervised segmentation of hyperspectral images using 3-d convolutional autoencoders.  IEEE Geosci Remote Sens Lett, 17(11):1948–1952, 2020.
  * [NYC15] A Nguyen, J Yosinski, and J Clune.  Deep neural networks are easily fooled: High confidence predictions for unrecognizable images.  In Proc IEEE CVPR, pages 427–436, 2015.
  * [OKA18] Savas Ozkan, Berk Kaya, and Gozde Bozdagi Akar.  Endnet: Sparse autoencoder network for endmember extraction and hyperspectral unmixing.  IEEE Trans Geosci Remote Sens, 57(1):482–496, 2018.
  * [PM21] Sam L. Polk and James M. Murphy.  Multiscale clustering of hyperspectral images through spectral-spatial diffusion geometry.  In Proc IGARSS, pages 4688–4691, 2021.
  * [PMP+11] Antonio Plaza, Gabriel Martín, Javier Plaza, Maciel Zortea, and Sergio Sánchez.  Recent developments in endmember extraction and spectral unmixing.  In Opt Remote Sens, pages 235–267. Springer, 2011.
  * [Pri97] John C Price.  Spectral band selection for visible-near infrared remote sensing: spectral-spatial resolution tradeoffs.  IEEE Trans Geosci Remote Sens, 35(5):1277–1285, 1997.
  * [QWG+18] Ying Qu, Wei Wang, Rui Guo, Bulent Ayhan, Chiman Kwan, Steven Vance, and Hairong Qi.  Hyperspectral anomaly detection through spectral unmixing and dictionary-based low-rank decomposition.  IEEE Trans Geosci Remote Sens, 56(8):4391–4405, 2018.
  * [RCY11] K Rohe, S Chatterjee, and B Yu.  Spectral clustering and the high-dimensional stochastic blockmodel.  Ann Stat, 39(4):1878–1915, 2011.
  * [RL14] Alex Rodriguez and Alessandro Laio.  Clustering by fast search and find of density peaks.  Science, 344(6191):1492–1496, 2014.
  * [RS00] S T Roweis and L K Saul.  Nonlinear dimensionality reduction by locally linear embedding.  Science, 290(5500):2323–2326, 2000.
  * [RZMC11] M A Rohrdanz, W Zheng, M Maggioni, and C Clementi.  Determination of reaction coordinates via locally scaled diffusion map.  J Chem Phys, 134(12):03B624, 2011.
  * [SLP+19] Y Su, J Li, A Plaza, A Marinoni, P Gamba, and S Chakravortty.  DAEN: Deep autoencoder networks for hyperspectral unmixing.  IEEE Trans Geosci Remote Sens, 57(7):4309–4321, 2019.
  * [SWW+20] J Sun, W Wang, X Wei, L Fang, X Tang, Y Xu, H Yu, and W Yao.  Deep clustering with intraclass distance constraint for hyperspectral images.  IEEE Trans Geosci Remote Sens, 59(5):4135–4149, 2020.
  * [SZS+14] C Szegedy, W Zaremba, I Sutskever, J Bruna, D Erhan, I Goodfellow, and R Fergus.  Intriguing properties of neural networks.  In Proc ICLR, 2014.
  * [TBC09] Yuliya Tarabalka, Jón Atli Benediktsson, and Jocelyn Chanussot.  Spectral–spatial classification of hyperspectral imagery based on partitional clustering techniques.  IEEE Trans Geosci Remote Sens, 47(8):2973–2987, 2009.
  * [TdSL00] J B Tenenbaum, V de Silva, and J C Langford.  A global geometric framework for nonlinear dimensionality reduction.  Science, 290(5500):2319–2323, 2000.
  * [TK06] S Theodoridis and K Koutroumbas.  Pattern recognition.  Elsevier, 2006.
  * [TRK11] Konstantinos E Themelis, Athanasios A Rontogiannis, and Konstantinos D Koutroumbas.  A novel hierarchical Bayesian approach for sparse semisupervised hyperspectral unmixing.  IEEE Trans Signal Process, 60(2):585–599, 2011.
  * [VDSN+18] D Van Dijk, R Sharma, J Nainys, K Yim, P Kathail, A J Carr, C Burdziak, K R Moon, C L Chaffer, and D Pattabiraman.  Recovering gene interactions from single-cell data using data diffusion.  Cell, 174(3):716–729, 2018.
  * [WFG+18] Li Wang, Yan Feng, Yanlong Gao, Zhongliang Wang, and Mingyi He.  Compressed sensing reconstruction of hyperspectral images based on spectral unmixing.  IEEE J Sel Top Appl Earth Obs Remote Sens, 11(4):1266–1284, 2018\.
  * [Win99] Michael E Winter.  N-FINDR: An algorithm for fast autonomous spectral end-member determination in hyperspectral data.  In Proc Imaging Spectr V, volume 3753, pages 266–275. SPIE, 1999\.
  * [WNW+19] R Wang, N Nie, Z Wang, F He, and X Li.  Scalable graph-based clustering with nonnegative relaxation for large hyperspectral image.  IEEE Trans Geosci Remote Sens, 57(10):7352–7364, 2019.
  * [WP17] H Wu and Saurabh Prasad.  Semi-supervised deep learning using pseudo labels for hyperspectral image classification.  IEEE Trans Image Process, 27(3):1259–1270, 2017.
  * [YCI13] Naoto Yokoya, Jocelyn Chanussot, and Akira Iwasaki.  Nonlinear unmixing of hyperspectral data using semi-nonnegative matrix factorization.  IEEE Trans Geosci Remote Sens, 52(2):1430–1437, 2013.
  * [YYL+18] Xiaofei Yang, Yunming Ye, Xutao Li, Raymond YK Lau, Xiaofeng Zhang, and Xiaohui Huang.  Hyperspectral image classification with deep learning models.  IEEE Trans Geosci Remote Sens, 56(9):5408–5423, 2018.
  * [ZCT+17] Wei Zhu, Victoria Chayes, Alexandre Tiard, Stephanie Sanchez, Devin Dahlberg, Andrea L Bertozzi, Stanley Osher, Dominique Zosso, and Da Kuang.  Unsupervised classification in hyperspectral imagery with nonlocal total variation and primal-dual hybrid gradient algorithm.  IEEE Trans Geosci Remote Sens, 55(5):2786–2798, 2017.
  * [ZG07] A Zare and P Gader.  Sparsity promoting iterated constrained endmember detection in hyperspectral imagery.  IEEE Geosci Remote Sens Lett, 4(3):446–450, 2007.
  * [ZRMC11] W Zheng, M A Rohrdanz, M Maggioni, and C Clementi.  Polymer reversal rate calculated via locally scaled diffusion map.  J Chem Phys, 134(14):144109, 2011.
  * [ZS14] Z Zhao and A Singer.  Rotationally invariant image representation for viewing direction classification in cryo-EM.  J Struct Biol, 186(1):153–166, 2014.
  * [ZSZ+18] X Zhang, Y Sun, J Zhang, P Wu, and L Jiao.  Hyperspectral unmixing via deep convolutional neural networks.  IEEE Geosci Remote Sens Lett, 15(11):1755–1759, 2018.
  * [ZWL+16] Liansheng Zhuang, Jingjing Wang, Zhouchen Lin, Allen Y Yang, Yi Ma, and Nenghai Yu.  Locality-preserving low-rank representation for graph construction from nonlinear manifolds.  Neurocomputing, 175:715–722, 2016.