# Evolution of $K$-means solution landscapes with the addition of dataset
outliers and a robust clustering comparison measure for their analysis

L. Dicks  Yusuf Hamied Department of Chemistry, Lensfield Road, Cambridge CB2
1EW, United Kingdom  IBM Research Europe, Hartree Centre, Sci-Tech Daresbury,
United Kingdom  D. J. Wales<EMAIL_ADDRESS>Yusuf Hamied Department of
Chemistry, Lensfield Road, Cambridge CB2 1EW, United Kingdom

###### Abstract

The $K$-means algorithm remains one of the most widely-used clustering methods
due to its simplicity and general utility. The performance of $K$-means
depends upon location of minima low in cost function, amongst a potentially
vast number of solutions. Here, we use the energy landscape approach to map
the change in $K$-means solution space as a result of increasing dataset
outliers and show that the cost function surface becomes more funnelled.
Kinetic analysis reveals that in all cases the overall funnel is composed of
shallow locally-funnelled regions, each of which are separated by areas that
do not support any clustering solutions. These shallow regions correspond to
different types of clustering solution and their increasing number with
outliers leads to longer pathways within the funnel and a reduced correlation
between accuracy and cost function. Finally, we propose that the rates
obtained from kinetic analysis provide a novel measure of clustering
similarity that incorporates information about the paths between them. This
measure is robust to outliers and we illustrate the application to datasets
containing multiple outliers.

$K$-means, clustering, outliers, landscapes, topography, cluster comparison

##  I  Introduction

Clustering is ubiquitous in data analysis and, consequently, many diverse
clustering methods have been developed. Amongst the most popular is the hard,
partitional $K$-means algorithm.Lloyd1982  The popularity of $K$-means stems
from its efficiency, algorithmic simplicity, and its ability to produce
clusterings in agreement with more complex and expensive algorithms.Jain2004 ;
deSouto2008  $K$-means clusters a given dataset into a prespecified number of
clusters, $K$. For some given initial cluster positions, $K$-means produces a
valid clustering solution by minimising the sum-of-squares cost function until
a local minimum is attained. The cost function surface is non-convex and can
support many local minima, which correspond to clusterings of varying quality.
The number of local minima for the $K$-means cost function can be vast, and
enumeration is impossible for all but trivial cases.Mahajan2012

Many schemes for initialising cluster centres have been proposed, which are
usually repeatedly applied and, after minimisation of each initial position,
the best clustering solution is retained. Many popular initialisation schemes
use uniform random samples of the data.Steinley2003 ; Bradley1998 ; Faber1994
Such methods have a low computational cost, allowing many uncorrelated initial
cluster positions to be generated. More complex initialisation schemes draw
non-uniform data samples,Arthur2007 ; Mohammad2009 ; Bachem2016  or pair with
other algorithms, such as hierarchical clustering,Milligan1980  or principal
component analysis.Su2004  These schemes generally have a greater
computational cost, but produce local minima of higher quality, allowing fewer
distinct initialisations to be considered. Additionally, it is possible to use
known solutions to guide future searching, rather than generating uncorrelated
trials, and $K$-means has been paired with genetic algorithms,Krishna1999 ;
Franti2000 ; Marghny2014  Monte Carlo-style schemes,Hand2005  and various
other optimisation algorithms.Hatamlou2012 ; Yao2016 ; Das2018 ; Manju2018 ;
Xie2019

The relative performance of all the above methods depends heavily on the
dataset structure and the number of clusters.Steinley2007 ; Filippone2008 ;
Rai2010 ; Celebi2013  The effect of dataset structure on $K$-means performance
is encoded by the cost function topography, as the ability of initialisation
schemes to locate low-valued $K$-means minima depends upon the organisation of
the solution space. The cost function can be coarse-grained in terms of the
local minima and the transition states that connect them, as in the energy
landscape framework.Wales03  The structure, thermodynamics, and kinetic
properties of physical systems can be described from these stationary points,
and we have recently extended this approach to explore the structure of the
$K$-means cost function.Dicks2022 ; Wu2023

The dependence of clustering accuracy on dataset structure means that datasets
are frequently modified prior to clustering to change the structure to favour
higher accuracy.Luai2006  This preprocessing commonly includes removal of
outliers, defined as data points that differ significantly from the rest,
which can badly degrade the quality of clustering that is
achievable.GarciaEscudero1999 ; Zhang2003  Hence care is often taken to
identify, and remove, outliers prior to clustering. However, there are various
methods to make $K$-means more robust to outliers.Mohammad2009 ; Xu2016 ;
Gupta2017 ; Huang2018 ; Brodinova2019 ; Zhang2020  Clustering itself can
reveal outliers, and there are many examples of outlier detection during
clustering.AlZoubi2009 ; Jiang2001 ; Marghny2014 ; Zhang2021  The effect of
outliers on clustering accuracy is well studied, but the effect they have on
the organisation of the solution space is unknown, and this is the question we
address in this contribution.

We first describe the algorithms of the energy landscape approach, and the
extensions required to explore cost function surfaces. Using this framework we
analyse the change in the $K$-means solution space that results from the
presence of dataset outliers, and for two datasets we find that the landscape
becomes more funnelled. Subsequently, we determine the properties of
clustering solutions using kinetic analysis to further examine the nature of
the funnelling. Lastly, we propose the use of metrics based upon rates as a
robust cluster comparison measure that will be especially useful in describing
global optimisation properties.

##  II Methods

###  II.1 Local minima

The $K$-means algorithm aims to minimise the sum-of-squares cost function,

$J(\bm{\mu})=\sum\limits^{N}_{i=1}\sum\limits^{K}_{k=1}r_{ik}\|\mathbf{x}_{i}-\mathbf{\bm{\mu}}_{k}\|^{2}.$
(1)

$\mathbf{x}$ is an ($N_{\mathrm{f}}\times N$)-dimensional matrix that contains
the $N$ data points, each of $N_{\mathrm{f}}$ features. $\bm{\mu}$ contains
the coordinates of all $K$ clusters. Each individual data point and cluster
are denoted by the $N_{\mathrm{f}}$-dimensional vectors $\mathbf{x}_{i}$ and
$\bm{\mu}_{k}$, respectively. For any given $\bm{\mu}$, data points are
assigned to the nearest cluster centre in Euclidean distance according to

$r_{ik}=\begin{cases}1,&\text{if
}k=\displaystyle\min\limits_{j}\|\mathbf{x}_{i}-\bm{\mu}_{j}\|^{2},\\\
0,&\text{if otherwise.}\end{cases}$ (2)

$r_{ik}$ is the element of the cluster assignment matrix $\mathbf{r}$ that
denotes the assignment of data point $i$ to cluster $k$.

The cost function, from any given cluster positions, can be locally minimised
using Lloyd’s algorithm.Lloyd1982  Initial cluster coordinates were generated
using uniform random points drawn within the range of the dataset in each
feature. This inefficient initialisation scheme was chosen to capture both
poor and good clustering solutions, as we aim to explore the complete solution
space. $10^{6}$ initial cluster positions were generated, and subsequently
minimised. Repeated minima were removed, along with those containing empty
clusters.

There are a variety of metrics for evaluating the clusterings generated by
$K$-means. Without known data point assignments, as is common in exploratory
data analysis, internal metrics can assess clustering quality. Internal
metrics use properties of the resulting clusters, such as their compactness
and separation. The central internal metric within $K$-means is the cost
function, which gives the within-cluster variance. However, various other
internal metrics have been proposed.Dunn1974 ; Davies1979 ; Rousseeuw1987

In contrast, external metrics measure the similarity of a clustering to a
given reference. For labelled datasets we refer to the known partitioning as
the ground truth clustering. These ground truth labels are a common reference,
in which case we measure the clustering accuracy. Popular external metrics can
be based upon cluster matchings,Rand1971 ; Hubert1985 ; Fowlkes1983
information theory,Meila2007  or meta-clusters.Cazals2019  We use the adjusted
Rand index (ARI)Rand1971 ; Hubert1985  throughout this work, which is
described in detail in the SI. The adjusted Rand index ranges from $-1$ to 1,
with a value of 1 indicating perfect agreement of two clusterings, and 0
indicating the similarity expected between two random clusterings. Values less
than 0 indicate a lower similarity than expected for random clusterings.

###  II.2 Transition states

For a detailed description of the methodology we employ to explore the
$K$-means cost function topography we refer readers to Ref. (27). A brief
summary is presented here.

Transition states provide the lowest barrier in the intermediate region of a
surface between one minimum and another and, therefore, provide valuable
information about the organisation. There are no transition states, defined as
Hessian index one saddles,Murrell1968  on $K$-means cost function
surfaces.Dicks2022  However, minimum-energy crossing points (MECPs), which lie
at the minimum of an intersection seam between two quadratics defined by
different cluster assignments, can be used as transition states in the present
context. At all points on a seam the Hessian is undefined, but an MECP can
have the correct curvature as a maximum in the direction orthogonal to the
intersection seam, and a minimum with respect to all directions along the
intersection seam.

Most transition state location algorithmsMunro1999 ; Henkelman2000 ;
Peters2004  require a function with smooth derivatives, which is not the case
for $K$-means. Instead, we employ a penalty-constrained MECP optimisation
algorithm.Levine2008  Location of an MECP between a pair of cluster
assignments, $\bm{r}^{(1)}$ and $\bm{r}^{(2)}$, which differ in a single data
point assignment, can be achieved through minimisation of a surrogate
function,

$F_{+}(\bm{\mu},\bm{r}^{(1)},\bm{r}^{(2)},\sigma,\alpha)=\frac{1}{2}\Big{(}J^{(1)}+J^{(2)}\Big{)}+\sigma\left(\frac{\Big{(}J^{(1)}-J^{(2)}\Big{)}^{2}}{\|J^{(1)}-J^{(2)}\|+\alpha}\right).$
(3)

$J^{(j)}=J(\bm{\mu},\bm{r}^{(j)})$ is the $K$-means cost function, evaluated
at cluster position $\bm{\mu}$, assuming a fixed cluster assignment
$\bm{r}^{(j)}$. $\sigma$ and $\alpha$ are parameters that determine the
increase in cost function away from the intersection seam, and smooth the
surface to avoid a discontinuity there, respectively. $\alpha=0.02$ and
$\sigma=30$ were used throughout this work, which produced a difference
between the cost function evaluated under cluster assignments $\bm{r}^{(1)}$
and $\bm{r}^{(2)}$ of $\mathcal{O}(10^{-4})$. This difference was sufficiently
small for our application, and permitted $\mathcal{O}(10^{5})$ transition
state location attempts, but the intersection seam can be found to
progressively higher accuracy by reminimising $F_{+}$ with small increments
applied to $\sigma$. Minimisations were performed using a customised L-BFGS
algorithm.Nocedal1980 ; Liu1989

The minimum of $F_{+}$, our approximation to the MECP, was then evaluated on
$J$ to give the cost function at the transition state between clusterings
$\bm{r}^{(1)}$ and $\bm{r}^{(2)}$. Transition states have two connected
minima, found by following (approximate) steepest-descent paths beginning in
the forwards and backwards direction of the eigenvector corresponding to the
unique negative Hessian eigenvalue. For $K$-means transition states the
required eigenvector must be determined from an alternative function. Here we
use Eq. (3) with subtraction, rather than addition, of the second term. The
function, denoted $F_{-}$, penalises the two quadratics when they are similar
in value and has the desired negative curvature in the direction perpendicular
to the intersection seam. We evaluate the Hessian on the $F_{-}$ surface, at
the minimum of $F_{+}$, to find the negative eigenvalue and its corresponding
eigenvector.

###  II.3 Building a stationary point network

We have summarised the methodology for finding transition states between
neighbouring minima, but to connect a large number of minima we need a scheme
to select appropriate pairs for transition state searches. Pairs are selected
using a distance-based criterion,RoederW18a  which finds the lowest-valued
minimum not currently connected to the global minimum by any sequence of
transition states and connected minima, and selects the nearest clustering
solution in Euclidean distance that is connected to the global minimum.

For each chosen pair of minima, we do not directly attempt a transition state
search, as they may be distant in cluster assignment. Instead, we construct a
linear path composed of discrete, equally-spaced cluster positions, referred
to as images. The number of images is adapted such that the cluster assignment
between any two adjacent images differs by only a single data point, and at
such a change in cluster assignment, a transition state search is attempted.
$F_{+}$ is constructed using the differing cluster assignments of the two
images, and a minimisation is started from the first image. If a valid MECP is
found we compute the Hessian on $F_{-}$ and identify the connected minima by
LBFGS minimisation. The process is repeated for all differing pairs of images,
and the valid transition states are added to a stationary point database,
along with their connectivity.

The aim of successive transition state searches is to produce a fully
connected set of minima, where it is possible to move between any pair via a
series of transition states and connected minima (a discrete pathWales2002 ;
Wales2004 ). The repeated application of the distance-based criterion metric
to select minima leads to efficient growth of the global minimum connections
until all minima are connected. The resulting database of minima and
transition states can be represented as a weighted graph, in which minima are
nodes, and all pairs of minima connected by a transition state are joined by
an edge weighted by a function of the barrier height. For physical systems,
this database corresponds to a kinetic transition network,Noe2008 ; Carr2009 ;
Prada2009  which we visualise using disconnectivity graphs.Becker1997 ;
Wales1998  Disconnectivity graphs illustrate the minima and the barriers
separating them for functions in arbitrary dimensions, as described in Sec.
SI.

We construct the analogue of a kinetic transition network for the $K$-means
cost function surface to calculate transition rates between groups of minima.
Such rates have no physical analogue, unlike molecular systems, but they
nevertheless provide a useful estimate of the difficulty in moving between
different regions of cluster solution space. The rates account for all the
intervening barriers in the cost function and all the possible pathways via
any number of intervening local minima. These properties are important in
understanding the $K$-means solution space, and contain important additional
information that is not present in distance calculations of either cluster
positions or assignments.

We compute the rate constant between two minima that are directly connected by
a transition state, using the analogue of unimolecular rate theory.Eyring1935
; Evans1935  The vibrational density of states that appears in this elementary
rate theory is set to unity for our $K$-means applications. The description of
global kinetics for a complete network is based upon the combination of the
rates of these individual steps, and the graph transformation
algorithmTrygubenkoW06b ; Wales2009  was used to calculate all the overall
rates reported in this work. We also extract the fastest path, when
intervening minima are considered in steady state, using Dijkstra’s
algorithmDijkstra1959  with suitable edge weights.Wales2002

The temperature parameter strongly affects the observed pathways and rates.
Higher temperatures lead to faster rates, and pathways that contain larger
barriers. Each cost function landscape was analysed at a fictitious
temperature that produces observed rates in the range between
$\mathcal{O}(10^{-6})-\mathcal{O}(10^{-3})$, which ensures that we favour
pathways with small barriers, and see a clearer separation in rates. The cost
function landscape itself is independent of temperature.

###  II.4 Datasets

We construct $K$-means landscapes for variations of two labelled datasets:
Fisher’s Iris datasetFisher1936  and the glass identification
dataset,Evett1989  both taken from the UCI machine learning data
repository.Dua2019  Fisher’s Iris dataset contains 150 data points of four
features each. There are fifty measurements of each of the three Iris species:
Virginica, Versicolor and Setosa. The Setosa data points form a clearly-
separable cluster, but measurements of the Versicolor and Virginica species
overlap, and cannot be separated by the $K$-means algorithm.

The glass dataset contains 214 data points separated into six clusters of
uneven populations. The dataset has nine features, which exhibit uneven ranges
that vary by two orders of magnitude. The ground truth clusters have
significant overlap, and the accuracy achievable by $K$-means is limited for
this dataset.

Outliers were sequentially added to both datasets, four to the Iris dataset
and three to the glass identification dataset. The addition of each successive
outlier produced a distinct dataset, giving four Iris datasets with an
increasing number of outliers. Outliers were placed distant from both the
original data points and any existing outliers. The distance between each
outlier and the mean of the overall data was not constrained, but is similar
for all outliers. The coordinates were selected at random, and are given in
Tables SI and SII.

##  III Results

Outliers modify $K$-means clustering to produce two alternative types of
solutions:

  1. 1.

The outlier belongs to a cluster containing no data points from the original
dataset.

There remain only $K-O$ clusters to represent the original dataset, where $O$
is the number of outlier clusters. $K$-means performance is degraded due to
poorer partitioning of the original data, and the degradation is greater when
the original dataset is composed of $K$ well-separated clusters with little
overlap.

  2. 2.

The outlier belongs to a cluster containing some original data points.

If an outlier does not belong to its own cluster then it must be accommodated
in a cluster containing original data points. The outlier may skew the cluster
centre a significant distance away from the mean of the original data points,
due to its remote position.

The two outlier solutions distinguished above can be used to understand the
structure of the clustering solutions within $K$-means landscapes. For
datasets with multiple outliers, each can correspond to one of the two cases,
and we provide a classification of the possible structures in Fig. 1.

Figure 1: The structure type scheme used to distinguish clusterings containing
outliers. Each cross represents an outlier in the $K$-means landscape,
irrespective of position in space. Circles surrounding one or more outliers
indicate membership of a distinct cluster with no original data points. All
other outliers are incorporated into clusters with original data points. The
scheme is shown for four outliers, but the ordering remains the same for
datasets with fewer than four.

The structure type, in our specified scheme, is independent of the outlier
performing the function. For example, we do not consider which outlier
corresponds to its own cluster, giving four options for structure type 1 in a
dataset containing four outliers. Each different type 1 solution will have
different values of the cost function depending upon where the outlier lies in
relation to the original data. However, the scheme remains valuable for
understanding the relevant structure of the clusterings.

###  III.1 Cost function topography

The topography of the solution space is central to global optimisation of all
machine learning algorithms, and for $K$-means it determines the ease of
locating the optimal clustering. First, we visualise the global structure of
the solution space, and describe its properties with reference to global
optimisation. Subsequently, we use kinetic analysis to further explore the
intermediate regions of solution space, and finally address the properties of
the clustering solutions.

####  III.1.1 Disconnectivity graphs

We present $K$-means landscapes, visualised using disconnectivity graphs, for
Fisher’s Iris dataset ($K=3$ and $K=6$) and the glass identification dataset
($K=6$) prepared with an increasing number of outliers in Figs. 2, 3 and 4,
respectively. The minima are coloured according to structure type as specified
in the previous section. Despite the large increase in the cost function range
across $K$-means landscapes with additional outliers, there is little change
in the barrier heights between transition states and their connected minima,
which generally remain very small throughout. Furthermore, the barriers to the
global minimum, which can involve pathways containing many minima and
transition states, also remain small for the vast majority of minima. This
structure produces a strong funnelling in the cost function surface. Single-
funnel landscapessocciow98  are characterised by efficient relaxation to a
well-defined global minimum, as for biomolecules that have evolved to perform
a single function.

(a) Original

(b) Outlier 1

(c) Outlier 2

(d) Outlier 3

(e) Outlier 4

Figure 2: $K$-means landscapes for Fisher’s Iris dataset with $K=3$ and an
increasing number of outliers. The scale bar represents the same cost function
range in all plots. Minima are coloured according to the scheme, given in Fig.
1, for distinguishing different types of $K$-means solutions.

(a) Original

(b) Outlier 1

(c) Outlier 2

(d) Outlier 3

(e) Outlier 4

Figure 3: Disconnectivity graphs showing the $K$-means solution landscapes for
variants of Fisher’s Iris dataset with six clusters. The scale bar represents
the same cost function range in each case, and minima are coloured according
to the scheme given for distinguishing different types of $K$-means solutions.

(a) Original

(b) Outlier 1

(c) Outlier 2

(d) Outlier 3

Figure 4: $K$-means solution landscapes for the glass identification dataset,
and variations with additional outliers. Minima are coloured according to the
structure types specified in Fig. 1. The scale bar corresponds to the same
cost function range in each case.

Despite the largely funnelled organisation of minima, there remain regions of
solution space separated from the global minimum by significant barriers in
all the $K$-means landscapes. These regions of the surface are referred to as
kinetic traps, because they slow relaxation to the global minimum in physical
systems, and impede global optimisation algorithms.Bryngelson1995 ;
Onuchic1997  The kinetic trap regions tend to retain a strong local funnelling
to their lowest minimum, and constitute regions of solution space that are
distinct from the global minimum.

The number of kinetic traps is generally comparable in $K$-means landscapes
corresponding to datasets with more outliers, despite the large increase in
the number of different structure types available to the clustering. However,
as the number of outliers increases there is a smaller proportion of minima
associated with these kinetic traps. Therefore, with the small barriers
retained for many minima, and fewer minima in kinetic traps, the landscapes
appear increasingly funnelled for a greater number of outliers.

####  III.1.2 Frustration

One important quantification of landscape funnelling is the frustration, which
we diagnose using the Shannon entropy,Shannon1948  which has previously been
for both molecular energy landscapesdeSouza2017  and $K$-means
landscapes.Wu2023  Increasing frustration values are associated with a
reduction in landscape funnelling due to more kinetic traps, larger barriers,
and more alternative low-valued solutions. Greater frustration highlights an
increased challenge for global optimisation and an increase in the number of
different regions of solution space that need to be considered during
clustering.

Shannon entropy profiles are shown for the Iris and glass landscapes in Fig.
5. The frustration metric is not presented for the Iris dataset with three
clusters, due to the limited number of minima. For both landscapes we see a
smooth decrease in frustration as the number of outliers increases. Therefore,
the stronger funnelling seen in the disconnectivity graphs is further
reflected by the frustration metric, which demonstrates that the addition of
outliers drives funnelling of the landscapes. Furthermore, for the glass
landscapes we observe multiple peaks in the frustration profiles. These
features correspond to strong fluctuations in either the degeneracy of minima,
or the size of their respective basins of attraction, for minima at
progressively higher cost function. Such fluctuations hint at the appearance
of distinct types of clustering becoming available, which would otherwise be
difficult to discern without analysing all the data point assignments.

(a) Iris ($K=6$)

(b) Glass

Figure 5: The Shannon entropy, used as a measure of landscape frustration,
calculated for each solution landscape from the Iris and glass datasets. The
entropy quantifies the complexity of the cost function surface as a function
of the fictitious temperature, where $k_{\mathrm{B}}T$ has the same units as
the cost function. Local minima with higher cost functions are progressively
occupied as the temperature increases.

####  III.1.3 Kinetic analysis

Disconnectivity graphs provide a representation of the overall barriers, but
the application of kinetic analysis produces a more detailed view of the
solution space organisation. Rates and pathways calculated between specified
minima include contributions from both the intermediate barriers and minima
that describe the space between solutions. We separate out two types of
pathways to the global minimum, those starting from minima within the main
funnel and those from kinetic traps.

For the majority of minima not within kinetic traps, the small barriers that
need to be overcome to reach the global minimum permit different structure
types to lie on a funnelled landscape. However, these small downhill barriers
do not imply that the pathways between two different structures and the global
minimum are similar. We compute the pathways from all different structure
types embedded in the funnel to the global minimum to probe the landscape
between the optimal solution and sub-optimal solutions. The cost function
profiles of the fastest pathways between each set of minima are visualised in
Fig. 6 for the Iris datasets, and Fig. S1 for the glass datasets.

(a) Outlier 1

(b) Outlier 2

(c) Outlier 3

(d) Outlier 4

Figure 6: The fastest paths from minima with a given structure type to the
global minimum for the Iris dataset with $K=6$, and its variations with
outliers. All stationary points on the pathway are shown, and changes in
structure type and partition are denoted by triangles and crosses on
transition states, respectively.

The pathways reveal repeated features, such as short plateau regions separated
by sharp decreases in cost function. These features are retained as outliers
are added, and show that the solution space remains composed of many shallow
funnelled regions within a broad funnel that encompasses almost all the space.
The shallow funnels correspond to a set of minima that belong to the same
underlying partitioning or structure type and differ only in the assignment of
a small number of data points, leading to small changes in cost function when
displacing the cluster positions between them. The overall funnel contains all
these regions corresponding to different partitioning and structure types. The
sharp changes in cost function arise due to regions between partitionings and
structure types that do not support any solutions, leading to long
minimisation pathways that result in a significant reduction in $J$.

Outliers permit longer pathways within a funnel, because there are additional
alternative structure types. Similar structure types retain short pathways,
even with four outliers, but a greater proportion of transitions within the
funnel pass through more intermediate shallow funnels to attain the global
minimum. The difference in these short and long pathways are highlighted in
Figs. S2 and S3 for the Iris dataset with four outliers.

For minima in kinetic traps there are much more significant barriers to the
global minimum along the pathways, as shown in Fig. 7 for the Iris landscape,
and Fig. S4 for the glass dataset. These large barriers are associated with
changes in partition and structure type. Intermediate regions of the pathways
still show shallow regions that indicate rearrangement within a partition or
structure type to allow for the larger structural transitions. The glass
pathways are much shorter than those for the Iris dataset, and have fewer
large barriers, which is likely due to the selected kinetic traps having the
same structure type as the global minimum for all but the Outlier 3 landscape.
The pathways on other landscapes still contain large barriers, and these are
likely due to changes in partition, which we cannot diagnose
straightforwardly.

(a) Original

(b) Outlier 1

(c) Outlier 2

(d) Outlier 3

(e) Outlier 4

Figure 7: The fastest paths between minima corresponding to the kinetic traps
given in Fig. S5 and the global minimum for the Iris dataset with $K=6$, and
variants with additional outliers. Paths are represented by minima, at integer
values along the horizontal axis, and the transition states between them.
These minima and transition states are connected by straight lines to guide
the eye. Changes in structure type are denoted by triangles, and changes in
partition by crosses, placed on the transition states that connect the two
differing minima.

####  III.1.4 Clustering solution properties

For a complete picture of the solution space it is important to determine the
properties of clustering solutions, such as accuracy. If minima lower in cost
function do not have higher accuracy then global optimisation may not produce
the best clusterings. The cost function favours solutions that have low
within-cluster variance, which may not necessarily correlate with high
accuracy. Accuracy is here computed as the ARI between the ground truth
cluster labels and the cluster labels of a given minimum. The accuracy
distribution of the minima is given in Fig. 8, and the disconnectivity graphs
coloured by accuracy are given in Figs. S7 and S8.

(a) Iris

(b) Glass

Figure 8: The distribution of clustering solution accuracy for the glass and
Iris dataset ($K=3$). The accuracy was quantified using the adjusted Rand
index (ARI) between each clustering and the ground truth labels.

For the Iris dataset with $K=3$ the cost function correlates strongly with
clustering accuracy, due to the small amount of overlap. The highest
achievable accuracy is 0.730, as the $K$-means algorithm is unable to separate
the Versicolor and Virginica clusters. There is little reduction in the best
accuracy on addition of increasing outliers, because the global minimum
accommodates the outliers by appending them to existing clusters (Fig. S9),
which results in little change in the cluster boundaries from the original
dataset. However, for landscapes of two outliers or more we see low-valued
solutions that correspond to one outlier cluster; the Iris dataset can be well
represented with two clusters due to the overlap of Versicolor and Virginica
data points.

The accuracy distribution of the glass identification dataset exhibits
different behaviour. This dataset contains significant overlap and, although
the accuracy still generally decreases with increasing cost function for the
original glass dataset, the $K$-means cost function is less appropriate for
clustering the dataset. The relationship between the cost function value and
accuracy becomes weaker with a growing number of outliers, and at three
outliers the accuracy has no relation to the cost function. Both the highest
attainable accuracy and the accuracy of the global minimum decrease with
additional outliers. The most accurate clustering moves to a higher cost
function value, relative to the global minimum, and the most accurate
clustering is not the global minimum in any landscape. With additional
outliers the most accurate solutions appear repeatedly at different values of
the cost function.

For the Iris dataset with $K=6$ we cannot evaluate the accuracy, but we can
rationalise the significant kinetic traps in the original Iris landscape by
equating them to different partitionings of the six $K$-means clusters amongst
the three underlying clusters. Partitionings remain useful to understand
kinetic traps, but as the number of outliers increases the structure type
exerts more control over cost function topography. Partitionings become less
clear, as many clusters now incorporate outliers, and a growing number of
kinetic traps can be rationalised by a difference in structure type, and not
partition. Disconnectivity graphs coloured by partitioning are given for Iris
datasets with both cluster numbers in Fig. S10, and the different structure
types of the global minima are shown in Fig. S11.

###  III.2 Rate-based clustering comparison

We propose that the rate computed between clustering solutions can be used as
a measure of their distinctness. All current external metrics for comparing
clusterings use only the initial and final configurations, but the rates
depend on the path between them, including information about the number of
intermediate minima and the barriers for cluster displacements. Hence this
metric encodes the cumulative change in cluster position and assignment,
rather than just the difference between endpoints. It is related to the
difficulty in navigating the landscape between two clusterings, providing
valuable insight into global optimisation performance for $K$-means. Moreover,
rates implicitly weight the importance of data points by the change they
induce in the cost function, which allows a robust measure of cluster
distinctness, even in the presence of dataset outliers.

We examine the correlation between rates, cluster position (distance), and
cluster assignment (ARI) measure in Fig. 9. The clustering similarity measures
were computed between different structure types in all outlier landscapes for
the Iris dataset with six clusters. The structure types differ in the
assignment of outliers, which can be challenging for cluster assignment-based
metrics because a few data point label changes produce a large effect.
Consequently, we observe a weak positive correlation with the ARI. The
distance captures the large distorting effect of the outliers on cluster
positions, and we observe a strong correlation with the rate, which therefore
also captures the importance of individual data points for cluster
distinctness. However, the three metrics probe different features of the
clustering distinctness, and there are many examples where the pathway-based
measure does not simply relate to the distance. We illustrate the comparative
measures alongside the clusterings, in three of the four dimensions, for the
same set of minima in Fig. S12.

Figure 9: Correlation between rates, $k_{12}$, and two other external metrics
for cluster comparison: the distance between cluster centres (left) and the
ARI (right). The pairs of clusterings are different structure types of the
Iris datasets containing outliers. The line of best fit is the solid black
line.

We compute the rates for various transitions that will affect global
optimisation, such as between different structure types and for kinetic traps
to the global minimum. These rates reflect the distinctness of various
clusterings from the global minimum, and hence report on the difficulty of
global optimisation from each starting point. The rates between different
structures embedded in the funnel can vary significantly, even for minima with
comparable cost functions, as shown in Tables 1 and SIV. The rates largely
correlate with the number of different structure types (or more rarely
partitions) that must be visited on the way to reaching the global minimum.
The rates remain similar across all outlier numbers when the structure types
involve a small number of changes, as the similarity of such structures is
retained. However, the rates can be much slower in landscapes with more
outliers, due to the greater difference in arrangements of distant structure
types. The rearrangement of multiple outliers involves more transition states
and intermediate minima, and the rates therefore decrease. The reduction in
slowest rate with outliers is much more pronounced for the glass dataset.
Furthermore, the rates remain remarkably similar between all structure types,
without the large preference for transitions between structures of the same
type seen in the Iris dataset.

Table 1: Transition rates (reduced units) between the set of minima of a given structure type to the global minimum, calculated at an effective temperature of $T=1.0$ (in reduced units consistent with the cost function). Rates are calculated for the Iris dataset with $K=6$, and those datasets with additional outliers. The dot indicates the structure type of the global minimum. Type |  Outlier 1 |  Outlier 2 |  Outlier 3 |  Outlier 4
---|---|---|---|---
0 |  $3.11\times 10^{-3}$ |  $8.82\times 10^{-4}$ |  $2.82\times 10^{-5}$ |  $1.68\times 10^{-7}$
1 |  $\cdot$ |  $1.08\times 10^{-2}$ |  $2.78\times 10^{-3}$ |  $1.94\times 10^{-7}$
2 |  |  $\cdot$ |  $\cdot$ |  $8.81\times 10^{-7}$
3 |  |  |  $8.68\times 10^{-4}$ |  $5.90\times 10^{-7}$
4 |  |  |  $5.11\times 10^{-4}$ |  $5.34\times 10^{-3}$
5 |  |  |  $2.85\times 10^{-3}$ |  $1.10\times 10^{-5}$
6 |  |  |  |  $2.12\times 10^{-7}$
7 |  |  |  |  $1.11\times 10^{-3}$
8 |  |  |  |  $\cdot$
9 |  |  |  |  $1.70\times 10^{-6}$

The rates between kinetic traps and the global minimum are given in Table 2
and Table SV for the Iris and glass datasets, respectively. These rates
generally decrease as the number of outliers increases. The escape rates from
different kinetic traps within the same landscape exhibit large fluctuations
as a result of the specific features. However, there remains an overall
reduction in escape rates and, consequently, escaping kinetic traps becomes
more difficult as the number of outliers increases.

Table 2: Escape rates (reduced units) from kinetic traps on the solution landscape for the Iris dataset with $K=6$, and variations with additional outliers. The rates are calculated from the set of minima at the bottom of the kinetic trap to the global minimum, at an effective temperature of $T=1.0$ (reduced units). The selected kinetic traps are given in Fig. S2. The geometric mean rate for each landscape is given in the final row. Kinetic trap |  Original |  Outlier 1 |  Outlier 2 |  Outlier 3 |  Outlier 4
---|---|---|---|---|---
1 |  $3.65\times 10^{-4}$ |  $2.25\times 10^{-5}$ |  $2.64\times 10^{-4}$ |  $3.23\times 10^{-4}$ |  $7.25\times 10^{-7}$
2 |  $1.74\times 10^{-5}$ |  $1.02\times 10^{-3}$ |  $6.55\times 10^{-8}$ |  $3.26\times 10^{-6}$ |  $1.12\times 10^{-6}$
3 |  - |  $2.66\times 10^{-5}$ |  $2.20\times 10^{-5}$ |  $2.26\times 10^{-9}$ |  $8.57\times 10^{-7}$
Mean |  $7.97\times 10^{-5}$ |  $8.48\times 10^{-5}$ |  $7.25\times 10^{-6}$ |  $1.34\times 10^{-6}$ |  $8.86\times 10^{-7}$

##  IV Conclusions

Successful $K$-means clustering depends on two factors: efficient location of
the global minimum (or low-valued minima), and the appropriateness of the cost
function, so that low-valued minima are the most accurate. The efficiency of
global optimisation is encoded in the organisation of the solution space. We
find that these landscapes have a funnelled structure for the datasets we
consider with any number of outliers. The degree of funnelling, as quantified
using a frustration metric, increases with the number of outliers.

Kinetic analysis reveals that the overall funnel seen in each landscape is
composed of many different shallow regions that correspond to different types
of structure. These shallow regions of space involve small barriers between
different structure types, and they are separated by cluster solution space
that does not contain any minima. Therefore, there can be large changes in
cost function between these regions associated with long pathways. The number
of such structure types increases with the number of outliers, which can
result in much longer pathways for navigation of the funnel.

We have introduced kinetic analysis of the $K$-means landscapes as a novel
clustering similarity measure. The rate extends external clustering metrics by
incorporating information about pathways between clustering solutions, rather
than just the endpoints. Such a measure can quantify the cumulative change in
cluster position and assignment, which may be more relevant for rationalising
global optimisation performance. We compute rates within the Iris and glass
datasets to analyse the similarity between different regions of the solution
space. We observe that high similarity is conserved between adjacent structure
types across all outlier landscapes, but the lowest similarity decreases as
the number of outliers increases.

A funnelled surface is usually favourable for global optimisation in physical
systems. However, the rates illustrate that it becomes more challenging to
both escape kinetic traps, and relax within the main funnel, to locate the
global minimum with increasing outliers. There is a sharp increase in the
number of ways in which outliers can be accommodated into clusterings as the
number of outliers increases, and this increase leads to longer pathways
between different regions of space. Therefore, although movement between
adjacent structure types remains fast, it becomes more challenging to move
about the main funnel between a greater number of structure types. For
successful global optimisation runs from distinct starting configurations may
be helpful, and an efficient moveset would require both local perturbations,
which optimise cluster boundaries, and larger changes between different
structure types by exchanging outliers. However, with an appropriate moveset
the single-funnel solution space should allow the global minimum to be located
quite efficiently.

The correlation between cost function and accuracy is degraded by the
increased number of outliers. The outliers perturb the cluster centre
positions, leading to reduced accuracy of the low-valued solutions.
Additionally, the increase in the number of structure types produces greater
variation in the accuracy throughout the cost function range. Therefore, for
datasets containing multiple outliers, obtaining the global minimum does not
necessarily give the most accurate clustering solution that $K$-means can
support. Without a clear way to distinguish accurate solutions in cost
function it is important to use an internal metric, such as the silhouette
coefficient,Rousseeuw1987  to independently rank them. However, except for the
case of three outliers in the highly-overlapping glass dataset, we observe
that the low-valued region does contain minima of comparable accuracy to the
most accurate clustering. Hence locating low-valued minima is useful, but it
becomes essential to generate an ensemble of clustering solutions, which can
capture the many possible arrangements of the various outliers.Alqurashi2019

##  V Supplementary material

See the supplementary material for additional methodological details (S-I and
S-III), a description of the datasets used in this work (S-II), and additional
figures referenced in this work (S-IV).

##  VI Acknowledgements

This work was supported by an EPSRC knowledge transfer fellowship. DJW
gratefully acknowledges an International Chair at the Interdisciplinary
Institute for Artificial Intelligence at 3iA Cote d’Azur, supported by the
French government [grant number ANR-19-P3IA-0002], which has provided
interactions that furthered the present research project.

## References

  * [1] S. Lloyd.  Least squares quantization in PCM.  IEEE Trans. Inf. Theory, 28:129–137, 1982.
  * [2] A. K. Jain, A. Topchy, M. H. C. Law, and J. M. Buhmann.  Landscape of clustering algorithms.  In Proc. of the IAPR Int. Conf. Pattern Recognit., pages 260–263, 2004.
  * [3] M. C. P. de Souto, I. G. Costa, D. S. A. de Araujo, T. B. Ludemir, and A. Schliep.  Clustering cancer gene expression data: a comparative study.  BMC Bioinform., 9:497, 2008.
  * [4] M. Mahajan, P. Nimbhorkar, and K. Varadarajan.  The planar $K$-means problem is NP-hard.  Theor. Comput. Sci., 442:13–21, 2012.
  * [5] D. Steinley.  Local optima in $K$-means clustering: what you don’t know may hurt you.  Psychol. Methods, 8:294–304, 2003.
  * [6] P. S. Bradley and U. M. Fayyad.  Refining initial points for $K$-means clustering.  Int. Conf. Mach. Learn., 1:91–99, 1998.
  * [7] V. Faber.  Clustering and the continuous $K$-means algorithm.  Los Alamos Science, 22:138–144, 1994.
  * [8] D. Arthur and S. Vassilvitskii.  $k$-means++: the advantages of careful seeding.  In Proc. of the 18th Ann. ACM-SIAM Symp. on Discrete Algorithms, pages 1027–1035, 2007.
  * [9] A. H. Mohammad, C. Vineed, S. Saeed, and J. Z. Mohammed.  Robust partitional clustering by outlier and density insensitive seeding.  Pattern Recognit. Lett., 30:994–1002, 2009.
  * [10] O. Bachem, M. Lucic, S. Hamed Hassani, and A. Krause.  Fast and provably good seedings for $k$-means.  In Proc. of the 30th Int. Conf. on Neural Information Processing Systems, pages 55–63, 2016.
  * [11] G. W. Milligan and P. D. Isaac.  The validation of four ultrametric clustering algorithms.  Pattern Recognit., 12:41–50, 1980.
  * [12] T. Su and J. G. Dy.  Another look at non-random methods for initializing $K$-means clustering.  In Proc. of the 16th IEEE Int. Conf. Tools Art. Intell., pages 784–786, 2004.
  * [13] K. Krishna and M. N. Murty.  Genetic $K$-means algorithm.  IEEE Trans. Syst. Man Cybern., 29:433–439, 1999.
  * [14] P. Fränti.  Genetic algorithm with deterministic crossover for vector quantization.  Pattern Recognit. Lett., 21:61–68, 2000.
  * [15] M. H. Marghny and A. I. Taloba.  Outlier detection using improved genetic $k$-means.  arXiv, 2014.
  * [16] D. J. Hand and W. J. Krzanowski.  Optimising $K$-means clustering results with standard software packages.  Comput. Stat. Data Anal., 49:969–973, 2005.
  * [17] A. Hatamlou, S. Abdullah, and H. Nezamabadi-pour.  A combined approach for clustering based on $k$-means and gravitational search algorithms.  Swarm Evol. Comput., 6:47–52, 2012.
  * [18] M. Yao, Q. Wu, J. Li, and T. Huang.  $k$-walks: clustering gene-expression data using a $k$-means clustering algorithm optimised by random walks.  Int. J. Data Mining Bioinf., 16:121–140, 2016.
  * [19] P. Das, D. K. Das, and S. Dey.  A modified bee colony optimization (MBCO) and its hybridization with $k$-means for an application to data clustering.  Appl. Soft Comput., 70:590–603, 2018.
  * [20] V. N. Manju and A. L. Fred.  AC coefficient and $K$-means cuckoo optimisation algorithm-based segmentation and compression of compound images.  IET Image Process., 12:218–225, 2018.
  * [21] H. Xie, L. Zhang, C. P. Lim, Y. Yu, C. Liu, H. Liu, and J. Walters.  Improving $K$-means clustering with enhanced firefly algorithms.  Appl. Soft Comput., 84:105763–105785, 2019.
  * [22] D. Steinley and M. J. Brusco.  Initializing $K$-means batch clustering: a critical evaluation of several techniques.  J. Classif., 24:99–121, 2007.
  * [23] M. Filippone, F. Camastra, F. Masulli, and S. Rovetta.  A survey of kernel and spectral methods for clustering.  Pattern Recognit., 41:176–190, 2008.
  * [24] P. Rai and S. Singh.  A survey of clustering techniques.  Int. J. Comput. Appl., 7:1–5, 2010.
  * [25] M. E. Celebi, H. A. Kingravi, and P. A. Vela.  A comparative study of efficient initialization methods for the $K$-means clustering algorithm.  Expert Syst. Appl., 40:200–210, 2013.
  * [26] D. J. Wales.  Energy Landscapes.  Cambridge University Press, Cambridge, 2003.
  * [27] L. Dicks and D. J. Wales.  Elucidating the solution structure of the $K$-means cost function using energy landscape theory.  J. Chem. Phys., 156:054109, 2022.
  * [28] Y. Wu, L. Dicks, and D. J. Wales.  Archetypal solution spaces for clustering gene expression datasets in identifying cancer subtypes.  arXiv, 2023.
  * [29] L. A. Shalabi, Z. Shaaban, and B. Kasasbeh.  Data mining: a preprocessing engine.  J. Comput. Sci., 2:735–739, 2006.
  * [30] L. A. García-Escudero and A. Gordaliza.  Robustness properties of $k$-means and trimmed $k$-means.  J. Am. Stat. Assoc., 94:956–969, 1999.
  * [31] J-S. Zhang and Y-W. Leung.  Robust clustering by pruning outliers.  IEEE Trans. Sys. Man Cybern., 33:983–999, 2003.
  * [32] J. Xu, J. Han, K. Xiong, and F. Nie.  Robust and sparse fuzzy $K$-means clustering.  In Proc. of the $25^{\mathrm{th}}$ Int. Joint Conf. Art. Intell., pages 2224–2230, 2016.
  * [33] S. Gupta, R. Kumar, K. Lu, B. Moseley, and S. Vassilvitskii.  Local search methods for $k$-means with outliers.  In Proc. of the VLDB Endowment, pages 757–768, 2017.
  * [34] S. Huang, Y. Ren, and Z. Xu.  Robust multi-view data clustering with multi-view capped-norm $k$-means.  Neurocomputing, 311:197–208, 2018.
  * [35] Š. Brodinová, P. Filzmoser, T. Ortner, C. Breiteneder, and M. Rohm.  Robust and sparse $k$-means clustering for high-dimensional data.  Adv. Data Anal. Classif., 13:905–932, 2019.
  * [36] X. Zhang, Y. He, Y. Jin, H. Qin, M. Azhar, and J. Z. Huang.  A robust $k$-means clustering algorithm based on observation point mechanism.  Complexity, 2020.
  * [37] M. B. Al-Zoubi.  An effective clustering-based approach for outlier detection.  Eur. J. Sci. Res., 28:310–316, 2009.
  * [38] M. F. Jiang, S. S. Tseng, and C. M. Su.  Two-phase clustering process for outliers detection.  Pattern Recognit. Lett., 22:691–700, 2001.
  * [39] Z. Zhang, Q. Feng, J. Huang, Y. Guo, J. Xu, and J. Wang.  A local search algorithm for $k$-means with outliers.  Neurocomputing, 450:230–241, 2021.
  * [40] J. C. Dunn.  Well-separated clusters and optimal fuzzy partitions.  J. Cybern., 4:95–104, 1974.
  * [41] D. L. Davies and D. W. Bouldin.  A cluster separation measure.  IEEE Trans. Pattern Anal. Mach. Intell., 1:224–227, 1979.
  * [42] P. J. Rousseeuw.  Silhouettes: a graphical aid to the interpretation and validation of cluster analysis.  Comput. Appl. Math., 20:53–65, 1987.
  * [43] W. M. Rand.  Objective criteria for the evaluation of clustering methods.  J. Am. Stat. Assoc., 66:846–850, 1971.
  * [44] L. Hubert and P. Arabie.  Comparing partitions.  J. Classif., 2:193–218, 1985.
  * [45] E. B. Fowlkes and C. L. Mallows.  A method for comparing two hierarchical clusterings.  J. Am. Stat. Assoc., 78:553–569, 1983.
  * [46] M. Meilă.  Comparing clusterings – an information based distance.  J. Multivar. Anal., 98:873–895, 2007.
  * [47] F. Cazals, D. Mazauric, R. Tetley, and R. Watrigant.  Comparing two clusterings using matchings between clusters of clusters.  J. Exp. Algorithmics, 24:1–41, 2019.
  * [48] J. N. Murrell and K. J. Laidler.  Symmetries of activated complexes.  Trans. Faraday Soc., 64:371–377, 1968.
  * [49] L. J. Munro and D. J. Wales.  Defect migration in crystalline silicon.  Phys. Rev. B, 59:3969–3980, 1999.
  * [50] G. Henkelman, B. P. Uberuaga, and H. Jónsson.  A climbing image nudged elastic band method for finding saddle points and minimum energy paths.  J. Chem. Phys., 113:9901–9904, 2000.
  * [51] B. Peters, A. Heyden, A. T. Bell, and A. Chakraborty.  A growing string method for determining transition states: comparison to the nudged elastic band and string methods.  J. Chem. Phys., 120:7877–7886, 2004.
  * [52] B. G. Levine, J. D. Coe, and T. J. Martínez.  Optimizing conical intersections without derivative coupling vectors application to multistate multireference second-order perturbation theory (MS-CASPT2).  J. Phys. Chem. B, 112:405–413, 2008.
  * [53] J. Nocedal.  Updating quasi-Newton matrices with limited storage.  Math. Comput., 35:773–782, 1980.
  * [54] D. C. Liu and J. Nocedal.  On the limited memory BFGS method for large scale optimization.  Math. Program., 45:503–528, 1989.
  * [55] K. Röder and D. J. Wales.  Energy landscapes for the aggregation of A$\beta_{17-42}$.  J. Am. Chem. Soc., 140:4018–4027, 2018.
  * [56] D. J. Wales.  Discrete path sampling.  Mol. Phys., 100:3285–3305, 2002.
  * [57] D. J. Wales.  Some further applications of discrete path sampling to cluster isomerization.  Mol. Phys., 102:891–908, 2004.
  * [58] F. Noé and S. Fischer.  Transition networks for modelling the kinetics of conformational change in macromolecules.  Curr. Opin. Struct. Biol., 18:154–162, 2008.
  * [59] J. M. Carr and D. J. Wales.  Refined kinetic transition networks for the GB1 hairpin peptide.  Phys. Chem. Chem. Phys., 11:3341–3354, 2009.
  * [60] D. Prada-Gracia, J. Gómez-Gardeñes, P. Echenique, and F. Falo.  Exploring the free energy landscape: from dynamics to networks and back.  PLoS Comput. Biol., 5:e1000415, 2009.
  * [61] O. M. Becker and M. Karplus.  The topology of multidimensional potential energy surfaces: theory and application to peptide structure and kinetics.  J. Chem. Phys., 106:1495–1517, 1997.
  * [62] D. J. Wales, M. A. Miller, and T. R. Walsh.  Archetypal energy landscapes.  Nature, 394:758–760, 1998.
  * [63] H. Eyring.  The activated complex and the absolute rate of chemical reactions.  Chem. Rev., 17:65–77, 1935.
  * [64] M. G. Evans and M. Polanyi.  Some applications of the transition state method to the calculation of reaction velocities, especially in solution.  Trans. Faraday Soc., 31:875–894, 1935.
  * [65] S. A. Trygubenko and D. J. Wales.  Graph transformation method for calculating waiting times in Markov chains.  J. Chem. Phys., 124:234110, 2006.
  * [66] D. J. Wales.  Calculating rate constants and committor probabilities for transition networks by graph transformation.  J. Chem. Phys., 130:204111–204118, 2009.
  * [67] E. W. Dijkstra.  A note on two problems in connexion with graphs.  Numerische Math., 1:269–271, 1959.
  * [68] R. A. Fisher.  The use of multiple measurements in taxonomic problems.  Ann. Eugen., 7:179–188, 1936.
  * [69] I. W. Evett and E. J. Spiehler.  Rule induction in forensic science.  Knowl. Based Sys., pages 152–160, 1989.
  * [70] D. Dua and C. Graff.  UCI machine learning repository, 2017.
  * [71] N. D. Socci, J. N. Onuchic, and P. G. Wolynes.  Protein folding mechanisms and the multidimensional folding funnel.  Proteins, 32:136–158, 1998.
  * [72] J. D. Bryngelson, J. N. Onuchic, N. D. Socci, and P. G. Wolynes.  Funnels, pathways, and the energy landscape of protein folding: a synthesis.  Proteins, 21:167–195, 1995.
  * [73] J. N. Onuchic, Z. Luthey-Schulten, and P. G. Wolynes.  Theory of protein folding: the energy landscape perspective.  Ann. Rev. Phys. Chem., 48:545–600, 1997.
  * [74] C. E. Shannon.  A mathematical theory of communication.  Bell Syst. Tech. J., 27:379–423, 1948.
  * [75] V. K. de Souza, J. D. Stevenson, S. P. Niblett, J. D. Farrell, and D. J. Wales.  Defining and quantifying frustration in the energy landscape: applications to atomic and molecular clusters, biomolecules, jammed and glassy systems.  J. Chem. Phys., 146:124103, 2017.
  * [76] T. Alqurashi and W. Wang.  Clustering ensemble method.  Int. J. Mach. Learn. & Cyber., 10:1227–1246, 2019.