# Elucidating the solution structure of the K-means cost function using energy landscape theory

Cite as: J. Chem. Phys. 156, 054109 (2022); https://doi.org/10.1063/5.0078793

Submitted: 15 November 2021 • Accepted: 16 January 2022 • Accepted Manuscript Online: 17 January

2022 • Published Online: 02 February 2022

<figure><img src="image_1.png" /></figure>

L. Dicks and

<figure><img src="image_2.png" /></figure>

D. J. Wales

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

View Online

<figure><img src="image_5.png" /></figure>

Export Citation

<figure><img src="image_6.png" /></figure>

CrossMark

## ARTICLES YOU MAY BE INTERESTED IN

* Symmetry-specific orientational order parameters for complex structures

The Journal of Chemical Physics 156, 054108 (2022); https://doi.org/10.1063/5.0076915

* Surrogate Hessian accelerated structural optimization for stochastic electronic structure theories

The Journal of Chemical Physics 156, 054104 (2022); https://doi.org/10.1063/5.0079046

* Equation-of-motion cavity quantum electrodynamics coupled-cluster theory for electron attachment

The Journal of Chemical Physics 156, 054105 (2022); https://doi.org/10.1063/5.0078795

Lock-in Amplifiers
up to 600 MHz

<figure><img src="image_7.png" /></figure>

Zurich
Instruments

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>
# Elucidating the solution structure of the K-means cost function using energy landscape theory

Cite as: J. Chem. Phys. 156, 054109 (2022); doi: 10.1063/5.0078793

Submitted: 15 November 2021 • Accepted: 16 January 2022 •

Published Online: 2 February 2022

<figure><img src="image_11.png" /></figure>

View Online

<figure><img src="image_12.png" /></figure>

Export Citation

<figure><img src="image_13.png" /></figure>

CrossMark

L. Dicks ID and D. J. Walesª) ID

## AFFILIATIONS

University Chemical Laboratories, Lensfield Road, Cambridge CB2 1EW, United Kingdom

## ABSTRACT

The K-means algorithm, routinely used in many scientific fields, generates clustering solutions that depend on the initial cluster coordinates. The number of solutions may be large, which can make locating the global minimum challenging. Hence, the topography of the cost function surface is crucial to understanding the performance of the algorithm. Here, we employ the energy landscape approach to elucidate the topography of the K-means cost function surface for Fisher's Iris dataset. For any number of clusters, we find that the solution landscapes have a funneled structure that is usually associated with efficient global optimization. An analysis of the barriers between clustering solutions shows that the funneled structures result from remarkably small barriers between almost all clustering solutions. The funneled structure becomes less well-defined as the number of clusters increases, and we analyze kinetic analogs to quantify the increased difficulty in locating the global minimum for these different landscapes.

© 2022 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0078793

## I. INTRODUCTION

Clustering is an important component of data analysis, where the aim is to partition a given dataset into groups whose members have high similarity. It can help us to reveal the principles that determine the underlying organization, providing insight into the phenomena that govern the data. Many clustering methods exist,¹ and the K-means algorithm² is among the most popular, with previous applications to chemical modeling,³ bioinformatics,⁴,⁵ image compression,⁶ and finance.⁷

For a given number of clusters and initial cluster positions, a K-means solution is obtained by minimizing the sum-of-squares cost function until a local minimum is attained. The number of local minima, which corresponds to clusterings of varying quality, can be large. Consequently, the quality of K-means solutions can be highly dependent on the scheme for generating initial cluster positions, and the location of the optimal solution is an NP-hard (non-deterministic polynomial-time hard) problem.⁸,⁹ Many initialization schemes have been developed,¹⁰⁻¹² with variations for the large datasets that are becoming commonplace.¹³⁻¹⁵

For any initialization scheme, the topography of the cost function surface governs our ability to locate good K-means solutions. The surface topography, described by the K-means minima and the

intermediate behavior of the cost function, depends on both the number of clusters and the dataset. An analysis of the convergence properties of Lloyd's algorithm can be used to describe the behavior of the cost function around local minima,¹⁶ and its geometrical properties.¹⁷ The organization of the minima themselves can be understood through the distinct types of K-means clustering solutions,¹⁸ but existing methods that determine the relation between minima for a given dataset rely on their statistical properties,¹⁹ without directly accessing the topography of the K-means cost function surface. The present contribution provides a new approach to this analysis.

The K-means cost function defines a surface analogous to the potential energy surfaces (PESs) routinely studied in physics and chemistry. The PES is a fundamental object of interest in molecular science, and consequently, many methods have been devised for exploring these high-dimensional landscapes. The computational tools developed in the energy landscape approach enable us to locate stationary points, minima, and the transition states between them. On the PES.²⁰ The application of this framework to the K-means cost function can explain the success, or failure, of K-means on given datasets and aid the location of good clustering solutions. The present contribution introduces the exploration of the K-means landscape and shows how the methodology developed to analyze

<footnote>ª) Author to whom correspondence should be addressed: dw34@cam.ac.uk</footnote>
molecular potential energy surfaces can be applied to gain new insight into this clustering procedure.

## II. METHODS

Fisher's Iris dataset,²¹ a small database containing 150 data points of four features each, was considered first in this study. The dataset contains an equal number of measurements of three Iris species and is routinely used to benchmark clustering algorithms.

In the energy landscape approach, the surface organization is described in terms of minima and transition states, which are stationary points where the gradient vanishes. Minima have only non-negative eigenvalues of the Hessian matrix, and any displacement of internal coordinates increases the cost function. Transition states, here defined as stationary points with a single negative Hessian eigenvalue,²² are a maximum with respect to one eigendirection and a minimum in the orthogonal directions that correspond to internal coordinate displacements.

For any fixed number of clusters, K-means solutions were located. The connectivity of these clustering solutions was revealed by finding the transition states between them, and with transition state information, a network of interconnected minima was built. The resulting networks were analyzed to determine the topography of the cost function surface. These steps are summarized in the following sections.

### A. Locating minima

The K-means algorithm partitions N data points, each with $N_f$ features, into a prespecified number of clusters, K. The aim of K-means clustering is to minimize the cost function,

$$J(\mu) = \sum_{i=1}^{N} \sum_{k=1}^{K} r_{ik} \|\mathbf{x}_i - \mathbf{\mu}_k\|^2, \quad (1)$$

where $\mathbf{x}$, an $(N_f \times N)$-dimensional matrix, contains the coordinates of the N fixed data points, and each $\mathbf{x}_i$, the ith column of $\mathbf{x}$, contains the $N_f$ components of data point i. $\mu$ contains the cluster coordinates, an $(N_f \times K)$-dimensional matrix, with $\mu_k$ (the kth column of $\mu$) being the position of cluster k in the $N_f$-dimensional space. The cluster positions are the optimizable parameters. $\mathbf{r}$ is an $(N \times K)$-dimensional matrix that indicates the assignment of each data point to the K clusters. $r_{ik}$ gives the assignment of data point i to cluster k, and each data point is assigned to the nearest cluster in Euclidean distance, according to

$$r_{ik} = \begin{cases} 1 & \text{if } k = \min_j \|\mathbf{x}_i - \mathbf{\mu}_j\|^2, \\ 0 & \text{if otherwise.} \end{cases} \quad (2)$$

Starting from any given cluster position, the cost function is locally minimized using Lloyd's algorithm,²³ which alternately: (i) calculates the optimal data point assignment, $\mathbf{r}$, for the current cluster positions, $\mu$, and (ii) places each cluster, $\mu_k$, at the center of the data points assigned to it. These two steps are repeated until no change in data point assignment is observed, and the resulting $\mu$ is a valid K-means solution. Each K-means solution is guaranteed

to be a local, but not global, minimum on the cost function surface, and many trials with different initial cluster coordinates may be necessary to find a good clustering solution.

A naïve initialization of cluster coordinates was employed to enumerate K-means minima for the benchmark in question. Uniform random points were drawn in the intervals $\{x_{j,min}, x_{j,max}\}$, where $x_{j,min/max}$ is the minimum/maximum extent in the jth direction. For each specified number of clusters, K, $10^6$ initial cluster coordinates were generated, and minima with empty clusters were discarded. We considered integer values of K between three and eight.

### B. Determining the connectivity of minima: The K-means landscape

#### 1. Cost function surface

The K-means cost function can be written in terms of all possible cluster assignments,

$$J(\mu) = \min_{j \in S} \left( \sum_{i=1}^{N} \sum_{k=1}^{K} r_{ik}^{(j)} \|\mathbf{x}_i - \mathbf{\mu}_k\|^2 \right), \quad (3)$$

where S is the set, of size $K^N$, containing all valid $\mathbf{r}$ matrices, those in which each data point is assigned to a single cluster. Here, $r_{ik}$ must be one for any of the K entries in each row i, with zeros otherwise, giving K possibilities for each row, so the N rows produce $K^N$ possible $\mathbf{r}$. The $r_{ik}^{(j)}$ that gives the minimum value of the cost function, for given $\mu$, is found by assigning each data point to the nearest cluster in the Euclidean distance. All $r_{ik}^{(j)}$ that correspond to permutations of members within the K clusters are considered equivalent.

The double summation in Eq. (3), for any fixed $\mathbf{r}^{(j)}$, describes a quadratic, and hence, there is a single quadratic form associated with each possible cluster assignment. Many of these cluster assignments do not appear on the cost function surface because they are not the optimal cluster assignment at any $\mu$ (unfeasible clusterings). Other cluster assignments can be optimal, but not in the region encompassing the minimum of their corresponding quadratic (non-minimal clusterings). K-means solutions are the minima of those quadratics that correspond to the optimal cluster assignment in the region that contains their minimum (minimal clusterings), as shown in Fig. 1.

#### 2. Transition state analogs

A transition state is defined as a point with a zero gradient and a single negative eigenvalue of the Hessian matrix.²² The spectrum of Hessian eigenvalues specifies that a transition state is a maximum with respect to one eigendirection and a minimum in all orthogonal directions, aside from modes corresponding to continuous symmetries with zero Hessian eigenvalues. There are no such modes for the K-means cost function considered here.

The gradient, for fixed $\mathbf{r}$, is

$$g(\mu) = -2 \sum_{i=1}^{N} \sum_{k=1}^{K} r_{ik}(\mathbf{x}_i - \mathbf{\mu}_k), \quad (4)$$
<figure><img src="image_14.png" /><figcaption>FIG. 1. The K-means cost function surface for two clusters (K = 2) and a one-dimensional dataset (N = 1) containing seven data items (N = 7). Intersections between different clusterings are denoted by red lines, and each clustering is labeled and visualized below. The data points are distributed along the one feature, x₁, and assignment to the two clusters is indicated by red or blue dots. The surface contains two minimal clusterings (r⁽²⁾ and r⁽⁴⁾), where the minima are shown by white dots. The remaining four clusterings on the surface are non-minimal. r⁽⁷⁾ is an example of an unfeasible clustering that does not appear on the K-means surface. The surface contains a plane of symmetry at μ₁ = μ₂, and reflected clusterings correspond to the permutation of the cluster labels. We only consider one of the possible permutations.</figcaption></figure>

and the associated Hessian is

$$H = \begin{pmatrix} 2n_1 \mathbf{I} & 0 & \cdots & 0 \\ 0 & 2n_2 \mathbf{I} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 2n_K \mathbf{I} \end{pmatrix}, \quad (5)$$

where I is an (N<sub>f</sub> × N<sub>f</sub>)-dimensional identity matrix and n<sub>j</sub> is the number of data points assigned to cluster j. The gradient and Hessian above are correct for any given cluster assignment. However, at a change in the cluster assignment, both the gradient and Hessian are undefined.

A change in the cluster assignment corresponds to the intersection of two quadratics on the cost function surface, defined by different cluster assignments. In almost all cases, the change in the cluster assignment corresponds to a single data point exchanged between two clusters. The condition for intersection is that the two clusters are equidistant from the data point they exchange, and all other data points retain the same cluster assignment; the cost function is then continuous across all intersections.

Because the gradient of the cost function is undefined at intersections between quadratics, no intersection seam can contain a transition state. Consequently, the cost function surface contains no transition states in the conventional sense, as the only points on the surface where the gradient vanishes are the minima of minimal clusterings. However, the minimum of an intersection seam, known as a minimum-energy crossing point (MECP), can function as a transition state in the present context. An MECP is minimal in all directions corresponding to motion along the intersection seam and can be a maximum in the remaining direction orthogonal to the intersection seam. We consider an MECP that is a maximum orthogonal to the intersection seam as a transition state analog (Fig. 2). Transition state analogs, despite not having the correct spectrum of Hessian eigenvalues, because the Hessian is undefined, have the appropriate curvature for a transition state. Non-maximal MECPs do not have the appropriate curvature and are not considered further in the analysis of networks.

A transition state analog, when it exists between two quadratics, is the maximum on the minimum energy pathway between the clusterings. Therefore, transition state analogs give the minimum height of the cost function barrier between two different clusterings. An intersection seam does not have to contain an

<figure><img src="image_15.png" /><figcaption>FIG. 2. A portion of the K-means surface given in Fig. 1. r⁽²⁾ and r⁽⁴⁾ are minimal clusterings, with minima (i) and (vi). For all labeled points, the positions of clusters μ₁ and μ₂ (larger black dots) and the corresponding cluster assignments are given below the surface. (iii) is the MECP and a transition state analog between states r⁽²⁾ and r⁽³⁾. (ii) and (iv) are non-MECP points on the intersection seam that corresponds to the motion of only one cluster from minimum (i). (v) is an MECP between r⁽³⁾ and r⁽⁴⁾ that is non-maximal in the direction orthogonal to the intersection seam.</figcaption></figure>
MECP because it is possible that another intersection seam may appear before the minimum is attained. In this case, the minimum height of the cost function barrier will occur at the triple intersection, which would correspond to the exchange of two data points simultaneously, and is, therefore, unlikely to occur in practice.

### 3. Locating transition state analogs

There are numerous methods for locating transition states,²⁴⁻²⁷ but many of them rely on the function being smooth, and with discontinuous first and second derivatives, these approaches cannot be applied to transition state analogs. Therefore, we employ a penalty-constrained MECP optimization algorithm, used for approximately locating MECPs on conical intersections.²⁸ A conical intersection, the intersection of two surfaces that represent electronic states, is analogous to an intersection of quadratics on the K-means cost function surface.

The location of the MECP between any pair of cost function quadratics, defined by cluster assignments $\mathbf{r}^{(1)}$ and $\mathbf{r}^{(2)}$, can be achieved through minimization of the objective function,

$$F_{+}(\mu, \mathbf{r}^{(1)}, \mathbf{r}^{(2)}, \sigma, \alpha) = \frac{1}{2} (J(\mu, \mathbf{r}^{(1)}) + J(\mu, \mathbf{r}^{(2)})) + \sigma \left( \frac{(J(\mu, \mathbf{r}^{(1)}) - J(\mu, \mathbf{r}^{(2)}))^2}{\|J(\mu, \mathbf{r}^{(1)}) - J(\mu, \mathbf{r}^{(2)})\| + \alpha} \right), \quad (6)$$

where $J(\mu, \mathbf{r}^{(j)})$ is the cost function, evaluated under cluster assignment $\mathbf{r}^{(j)}$, at cluster positions $\mu$. $\alpha$ is a smoothing parameter to avoid a discontinuity at the intersection seam. $\sigma$ determines how heavily to penalize the deviation of the current point from the intersection seam.

For finite $\sigma$, the minimum of $F_+$ will not lie at the true MECP. However, the intersection seam can be found to arbitrary accuracy by evaluating the difference between the two quadratics at the minimum of $F_+$ and, if not sufficiently small, minimizing using $\sigma = \sigma + \delta$, where $\delta$ is a small positive value. Here, for each $F_+$, minimization was performed using the L-BFGS (limited-memory Broyden,²⁹ Fletcher,³⁰ Goldfarb,³¹ Shanno³²) algorithm.³³,³⁴ Due to the initial positions being very close to the intersection seam, as described in Secs. II A–II D, we selected a large initial $\sigma$ value ($\alpha = 0.02$ and $\sigma = 30$). This choice of $\sigma$ gave a difference between quadratics of $O(10^{-4})$, which was considered sufficiently small without performing multiple minimizations at increasing values of $\sigma$, and permitted many more transition state analogs to be located on the cost function surface. The MECP, obtained by this minimization of $F_+$, was then evaluated for $J$ to give the cost function at the transition state analog between clusterings $\mathbf{r}^{(1)}$ and $\mathbf{r}^{(2)}$.

For each MECP, we calculate the optimal cluster assignment and check that it is the same as either the $\mathbf{r}^{(1)}$ or $\mathbf{r}^{(2)}$ between which it is a maximum. We only retain MECPs where this condition is true; otherwise, the MECP does not exist on the cost function surface because another intersection has been crossed, and a different lowest-valued quadratic is present at this point.

Each transition state has two connected minima, which are found by following (approximate) steepest-descent paths parallel and antiparallel to the eigenvector corresponding to the single negative Hessian eigenvalue. The Hessian is undefined at transition

state analogs, so to obtain the required eigenvector, an alternative function is constructed, namely,

$$F_{-}(\mu, \mathbf{r}^{(1)}, \mathbf{r}^{(2)}, \sigma, \alpha) = \frac{1}{2} (J(\mu, \mathbf{r}^{(1)}) + J(\mu, \mathbf{r}^{(2)})) - \sigma \left( \frac{(J(\mu, \mathbf{r}^{(1)}) - J(\mu, \mathbf{r}^{(2)}))^2}{\|J(\mu, \mathbf{r}^{(1)}) - J(\mu, \mathbf{r}^{(2)})\| + \alpha} \right), \quad (7)$$

where $\mathbf{r}^{(1)}, \mathbf{r}^{(2)}, \sigma$, and $\alpha$ are given the same values as in $F_+$ and the function only differs in penalizing, rather than favoring, the two quadratics being similar in value. Therefore, $F_-$ has the desired curvature about the intersection seam, with a negative curvature along the direction perpendicular to the intersection. At the minimum of $F_+$, we can evaluate the Hessian for $F_-$ to find the negative eigenvalue and corresponding eigenvector (Fig. 3). Approximate steepest-descent paths from transition state analogs were calculated using the L-BFGS algorithm.

It is worth noting that $F_-$ will not distinguish between transition state analogs and non-maximal MECPs. However, in almost all the cases, the steepest-descent paths from a non-maximal MECP will both end at the same minimum, as the path from the higher-valued clustering will, within a small number of steps, recross the intersection seam. Therefore, non-maximal MECPs provide no information about the connectivity between different minima on the cost function surface.

The steepest-descent paths from a transition state analog do not necessarily end at the minima of the two quadratics for which it is an MECP. The two quadratics between which the transition state analog exists may correspond to non-minimal clusterings, and therefore, the steepest-descent paths will pass through these quadratics before attaining a minimal clustering.

### C. Building the network

Using the methods described above, it is possible to locate K-means minima, transition state analog, and determine their connectivity. Beginning from a database containing all known K-means minima (Sec. II A), we identify the transition states that connect them (Sec. II B). The resulting stationary point database is analogous to a kinetic transition network in molecular science.³⁵⁻³⁸ The database can be represented as a weighted graph, where all minima are nodes and all pairs of minima connected by a transition state analog are joined by an edge weighted by a function of the barrier height.

For connection attempts between minima in the stationary point database, we select pairs using a distance-based criterion.³⁹ This criterion finds the lowest-valued minimum not currently in the subset containing the global minimum, and selects the nearest member of this subset in Euclidean distance. The repeated application of this metric to select minima leads to growth of the subset containing the global minimum until all minima are connected.

An initial linear path was constructed between the chosen pair of minima, composed of a discrete set of equally-spaced cluster positions, referred to as images. This initial path was perturbed through random displacements applied to images separated by a fixed
<figure><img src="image_16.png" /><figcaption>FIG. 3. Transition state analogs for two-dimensional ($N_f = 2, N = 7$) and three-dimensional ($N_f = 3, N = 9$) clustering examples with two clusters ($K = 2$). The assignment of data points to the two clusters, given in black, is indicated by red or blue, and the data point that changes its cluster membership is green. The transition state analogs were located by minimizing $F_+$, and the Hessian is evaluated at $F_^-$. In two dimensions, the intersection seam requires that both the clusters lie on a circle centered at the exchanged data point, and the MECP is the minimum with respect to radius and positions. In three dimensions, the equivalent condition corresponds to clusters lying on a sphere centered at the exchanged data point. In both the cases, a single negative eigenvalue is observed, and the corresponding eigenvector is aligned perpendicular to the intersection seam. The largest components of this eigenvector are indicated by arrows. The other Hessian eigenvectors, with positive eigenvalues, correspond to the motion along the intersection seam.</figcaption></figure>

interval. Each path segment between perturbed images remains linear. The number of images was adapted such that at most one data point changes cluster assignment between any two adjacent images.

Subsequently, the path is analyzed and, at any change in cluster assignment between adjacent images, the first image in the pair is considered a candidate for a transition state analog. For each candidate we construct the relevant $F_+$, using the differing cluster assignments, and minimize to obtain the MECP. Provided the number of images is large, the two images differing in cluster assignment will be close to the intersection seam, so that fewer minimization steps are required. For every valid MECP we subsequently evaluate the Hessian at $F_-$ and find the connected minima, Fig. 4. After evaluating all candidates, the valid transition state analog are added to the stationary point database, along with their connectivity, and the connected minima, if not already present in the database. Minima with unassigned clusters were discarded, apart from the $K = 7$ and $K = 8$ landscapes, where two minima with unassigned clusters were included to allow connection between higher partitionings that proved difficult to connect directly. Two

new minima are selected from the database using the distance-based criterion, and the process is repeated until all minima are connected.

For all $K$, after connection of all minima into a single set, we observe that a small number of transition state analogs are lower in cost function than the higher of their connected minima. A transition state analog, as a maximum orthogonal to the intersection seam, should not be lower in cost function than the minima found by steepest-descent paths. The cause of the reversal in cost function values is the imprecision in the location of MECPs. We expect to locate the intersection seam with an accuracy of $O(10^{-4})$, and the largest gap between a transition state analog and the connected minimum is $1.6 \times 10^{-4}$. Therefore, we add $2 \times 10^{-4}$ to the cost function of each transition state analog to leave the barriers almost unchanged, whilst removing the cases where the connected minimum is higher than the transition state analog, which makes further analysis more convenient. This procedure has no perceptible effect on the visualization of the landscapes presented below.
<figure><img src="image_17.png" /><figcaption>FIG. 4. A scheme for locating transition state analogs and their connected minima, illustrated for one transition state analog on the cost function surface of Fig. 2. The cost function (top left panel) contains one transition state analog, on the intersection seam between the clusterings $r^{(2)}$ and $r^{(3)}$, one non-maximal MECP, and two minima. $F_+$, constructed using $r^{(2)}$ and $r^{(3)}$, is minimized to obtain the white point (top right panel). The equivalent $F_-$ function is constructed (bottom left panel), and the Hessian is evaluated at the minimum of $F_+$ to obtain the eigenvector that we use to define the steepest-descent paths, shown with arrows. The original cost function, $J$, is evaluated at the minimum of $F_+$, and the steepest-descent paths are calculated to obtain the connected minima (bottom right panel). The paths for the forward and backward directions of the eigenvector are shown in orange and blue.</figcaption></figure>

### D. Visualizing the surface

The networks that encode the topography of the cost function landscape were visualized using disconnectivity graphs.40,41 Disconnectivity graphs, which can be constructed for surfaces of arbitrary dimensions, display all minima and faithfully represent the lowest barriers separating them. The barrier between two minima is given by the value of the cost function at the highest transition state analog on a path between them, which may involve a sequence of transition state analogs and connected minima.

The vertical axis of a disconnectivity graph gives the cost function, and a vertical line begins, for each minimum, at the corresponding cost function value. For any given threshold of the cost function, the minima can be separated into disjoint sets (super-basins40), whose members are all connected by barriers below the threshold. At very high values of the threshold, all the minima will belong to a single set, but as the threshold is lowered, the minima

separate into superbasins, containing progressively fewer minima. The superbasin analysis is performed at regular cost function intervals, and the lines corresponding to minima begin on the vertical axis at their respective cost function value and merge at the threshold where they become members of the same set. The horizontal axis is chosen to arrange the minima as clearly as possible.

### III. RESULTS

Fisher's *Iris* dataset has a simple structure, consisting of two clearly separated clusters of differing sizes. The smaller cluster contains data points from the species *Iris setosa*, and the larger cluster contains data from both *Iris virginica* and *Iris versicolor*. Separation of the larger cluster into the constituent species is not possible using the *K*-means algorithm.

*K*-means, at $K \ge 3$, is forced to distribute its clusters between the two underlying clusters of the *Iris* dataset in varying numbers.
These different partitionings of the *K* clusters give rise to distinct groups of clustering solutions. Within a given partitioning, the minima have only small variations in data point assignment relative to those in different partitionings.

The *K*-means landscapes are visualized using disconnectivity graphs in Fig. 5, revealing the presence of a well-defined global minimum for all *K*. The landscapes contain no minima with similar values of the cost function that are separated by large barriers. In the terminology of energy landscapes, this organization is considered unfrustrated.42–44 All global minima are accompanied by a set of low-lying solutions, similar in structure and separated by small barriers, which can be viewed as equivalent clusterings. Any member of this set provides a suitable substitute for the global minimum clustering, as they retain the same partitioning and only differ in cluster boundaries. In all the cases, the global minimum clustering and its associated set of minima correspond to the optimal partitioning of *K*-means clusters between the underlying clusters of the

dataset. However, the gap between partitionings decreases with *K*, as fewer data points are assigned to each cluster, which leads to alternative partitionings becoming more competitive with the global minimum.

The *K*-means landscapes are all quite funneled toward the global minimum, with some interesting structure. Funneled surfaces are associated with self-organizing systems, which exhibit efficient relaxation:45–47 all local minima can reach the global minimum via small barriers. The funneled topography that leads to efficient structure-seeking systems in nature also leads to easier global optimization. Therefore, *K*-means methods that use optimization algorithms to perturb and update local minima are expected to perform well on cost function surfaces such as these.48–53

Subfunnels, which contain minima with alternative morphologies to the global minimum, exhibit the same local topography that characterizes efficient optimization. They also have small downhill barriers from all minima within the subfunnel to the local funnel

<figure><img src="image_18.png" /><figcaption>FIG. 5. Disconnectivity graphs illustrating the changing topography of the *K*-means cost function surface with the cluster number, *K*, for Fisher's *Iris* dataset. The coloring of each minimum indicates its partitioning, given as the number of *K*-means clusters within the smaller of the two underlying clusters of the dataset, as indicated in the scale bar.</figcaption></figure>
bottom. Once within a subfunnel, optimization to its lowest-valued representative is efficient, but transitions to the global minimum funnel involve a much larger barrier. Consequently, subfunnels constitute kinetic traps because they hinder the location of the global minimum.

In the *K*-means landscapes, the subfunnels correspond to different partitionings from the global minimum. The minima within a subfunnel vary locally in data point assignment while retaining the same partitioning. More partitionings are possible for larger *K*, leading to an increase in the number of kinetic traps. The number of minima associated with each partitioning also increases with *K*, due to more minimal clusterings, and the relative number of minima in each subfunnel can vary significantly. For example, at *K* = 6, the number of minima for alternative partitionings of the lowest subfunnel vastly outnumbers the number of minima associated with the global minimum partitioning. In general, there is an increase in the number of minima in lower-lying subfunnels because they correspond to a more even distribution of data points between clusters, which allows them to support more minima, as seen in Fig. 6. Hence, these subfunnels exhibit greater landscape entropy,⁵⁴⁻⁵⁶ which is determined by the number of minima.

In molecular science, various surface topographies have been characterized, notably single-funnel,³⁷ multifunneled,³⁸⁻⁶⁰ and glassy.⁶¹ This classification of landscapes is immediately suggestive of observable properties, such as broken ergodicity and multiple relaxation time scales.⁶²,⁶³ Single-funnel surfaces support efficient optimization to a well-defined global minimum, and correspond to “magic number” clusters, and biomolecules that have evolved to optimize a single function. Multifunneled landscapes arise when

there are two, or more, funnels that are similar in cost function and separated by significant barriers. Multifunnelled surfaces are more challenging for global optimization, and locating the global minimum can be significantly hindered by alternative funnels acting as kinetic traps. In glassy landscapes, there is no clear global minimum if the crystal is excluded, as there are a large number of almost equivalent amorphous minima, corresponding to funnels that are separated by large barriers.⁶⁴,⁶⁵ For these surfaces, global optimization is very challenging, and minima from more than one of these funnels would be needed to fully describe the optimal clustering solutions.

All the cost function landscapes for the *Iris* dataset exhibit a well-defined funnel for the global minimum. However, the growing number of kinetics traps, and their competition with the global minimum, lead to a more multifunneled character at larger *K*. Therefore, as the number of clusters increases, the cost function landscape changes such that global optimization becomes more challenging, and the cost functions of alternative partitionings approach the value for the global minimum.

The ability of subfunnels to act as kinetic traps can be described in more detail using kinetic analysis,⁶⁶⁻⁶⁸ which provides transition rates or mean first passage times in physical systems. Here, we apply the same framework to the networks of the cost function landscapes. The density of states, which appears in the partition function for each minimum and the transition state theory rate constant formulation,⁶⁹,⁷⁰ was set to unity in the present calculations. Hence, we do not attempt to include the analog of a vibrational density of states, but simply use the appropriate Boltzmann factor to assign relative occupation probabilities for the local minima.

<figure><img src="image_19.png" /><figcaption>FIG. 6. Histograms for each *K*, showing the number of minima at values of the cost function between the highest and the lowest observed minima.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>FIG. 7. Rates calculated between the lowest minimum in each different partitioning, at all *K*. Every node represents a minimum and is labeled by its partitioning, the number of clusters in the smaller of the two underlying clusters of the dataset, and its cost function value. Rates are given in both directions between all different minima.</figcaption></figure>

Overall rates between selected minima in the database, accounting for all possible paths between them, were calculated using the graph transformation method.³⁷ The transition rates between minima in different subfunnels provide a quantitative measure of the difficulty of escaping kinetic traps using optimization algorithms. As for physical systems, we have to specify a temperature at which the rates are calculated. For the cost function surfaces, as the temperature has no physical meaning, we calculate rates at an arbitrary value (*T* = 3) in units consistent with the cost function. For all *K*, we calculated rates for the paths between the lowest-valued minima in each partitioning, which are the global minimum and the lowest minima in all subfunnels.

In almost all the cases, the rate for moving to lower-valued partitionings is higher than the reverse direction, as shown in Fig. 7. These rates provide insight into the efficient global optimization on these funneled landscapes, as they illustrate the preference for escaping kinetic traps and reaching the global minimum compared to the reverse direction. Additionally, we observe that the rates for reaching the global minimum, from all other subfunnels, decrease with increasing *K*, and the reverse rates become more competitive. Therefore, as the number of kinetic traps increases, so does the difficulty in escaping from them, illustrating the likely performance of global optimization algorithms. The lower rates between partitionings at higher *K* are due to the increased number of steps, rather than the

barriers that have to be overcome on the pathways, as shown in Table I. The barriers remain small for all the landscapes, but the number of steps increases due to the greater number of minimal clusterings.

We observe that, as expected, rates for the conversion between partitionings are significantly faster for adjacent partitionings, namely, those that correspond to a difference of one in the number of clusters assigned to the smaller of the underlying clusters of

<table><caption>TABLE I. The average path length, measured in terms of the number of transition state analogs. The average barrier height is tabulated for all transition state analogs, and the largest of these barriers, for pathways downhill between the lowest minima in different partitionings. Averages are only calculated for pathways between partitionings that differ by exchange of one cluster between the two underlying clusters of the dataset.</caption><thead><tr><th>K</th><th>Path length</th><th>Average barrier</th><th>Largest barrier</th></tr></thead><tbody><tr><td>3</td><td>5.0</td><td>0.77</td><td>2.53</td></tr><tr><td>4</td><td>2.5</td><td>1.05</td><td>2.62</td></tr><tr><td>5</td><td>7.0</td><td>1.03</td><td>4.19</td></tr><tr><td>6</td><td>12.0</td><td>0.74</td><td>3.38</td></tr><tr><td>7</td><td>12.5</td><td>0.61</td><td>3.45</td></tr><tr><td>8</td><td>19.8</td><td>0.30</td><td>3.02</td></tr></tbody></table>
<figure><img src="image_21.png" /><figcaption>FIG. 8. Histograms showing the distribution of barrier heights for each cost function landscape. The barrier height is calculated as the difference in cost function between a transition state analog and the higher, or lower, of its two connected minima. These barrier heights were calculated for every transition state analog in the corresponding kinetic transition network for $3 \le K \le 6$, and for the larger $K = 7$ and $K = 8$ databases, we selected a random subset of 50 000 transition state analogs.</figcaption></figure>

the dataset. This rate reduction arises because the fastest pathways between larger changes in partitioning pass through intermediate partitionings and are not connected directly. Consequently, the path length increases to incorporate the additional change in partitioning, and the additional barriers reduce the rate.

For all K, the transition state analogs are very similar in cost function to the higher of their two connected minima (Fig. 8). Hence, there are small barriers for each step on pathways between distant minima, which contribute to the funneling properties of the landscape. The distribution of barrier heights shows that the most probable height, at all K, is small, and transition state analogs with larger barriers become increasingly unlikely. For some surfaces, there are barrier heights at larger values (around 20). These barriers will occur when the transition state analogs are MECPs between two cluster assignments in a higher-valued partitioning, but both the steepest-descent paths end at minima of a lower-valued partitioning. The distributions indicate that most transition state analogs lie within the same partitioning as one of the connected minima.

The small barriers indicate that the vast majority of transition state analogs lie at, or close to, an intersection directly involving the higher of the two connected minima, Table II. Therefore, transition state analogs often correspond to the exchange of a single data point from the higher connected minimum. Furthermore, as a result of the high density of data points, the distance between the transition state analog and the higher connected minimum is small.

The transition state analogs closely follow Hammond's postulate that the transition state resembles the reactant or product most similar in energy,⁷¹ which can be explained using catastrophe theory.⁷²⁻⁷⁴ When the distance between the minimum and the transition state analog is very small, which is expected to occur at high data point density, we expect the corresponding stationary points to merge at a non-Morse point corresponding to a fold catastrophe.⁷²

<table><caption>TABLE II. An analysis of transition state analogs and their relation to their connected minima. Δ**r** and **D** are the average number of data points that vary in assignment and the average difference in Euclidean distance, between the transition state analog and its higher or lower connected minimum, respectively. For the higher connected minima, only transition state analogs with barrier heights in the range 0.0–0.1 were considered, and the proportion of the total transition state analogs these contribute is given by **P**. For the lower connected minima, we give values for the three transition state analogs with the largest barriers.</caption><thead><tr><th rowspan="2">K</th><th colspan="3">Higher</th><th colspan="2">Lower</th></tr><tr><th>Δr</th><th>D</th><th>P</th><th>Δr</th><th>D</th></tr></thead><tbody><tr><td>3</td><td>0.89</td><td>0.012</td><td>0.9</td><td>68.67</td><td>3.047</td></tr><tr><td>4</td><td>0.98</td><td>0.013</td><td>0.96</td><td>99.33</td><td>3.350</td></tr><tr><td>5</td><td>0.96</td><td>0.016</td><td>0.86</td><td>54.00</td><td>2.354</td></tr><tr><td>6</td><td>0.96</td><td>0.017</td><td>0.92</td><td>63.33</td><td>2.69</td></tr><tr><td>7</td><td>0.96</td><td>0.016</td><td>0.95</td><td>58.00</td><td>2.53</td></tr><tr><td>8</td><td>0.96</td><td>0.015</td><td>0.97</td><td>63.67</td><td>1.66</td></tr></tbody></table>
The steepest-descent path to the lower minimum can pass through one or many other clusterings. In most cases, the barrier from the lower connected minimum is only slightly larger than that for the corresponding higher minimum, as reflected in the distributions of barrier heights, which are shifted to slightly higher values. However, for transition state analogs that connect distant minima, the pathway to the lower-valued minimum must pass through multiple non-minimal clusterings and move through a large distance in the configuration space. The difference between the transition state analog and its lower-valued minimum in these cases can be large, as shown in Table II.

## IV. CONCLUSIONS

In this paper, we present a methodology for determining and visualizing the topography of *K*-means cost function surfaces. We applied the methodology to Fisher's *Iris* dataset, varying the number of clusters, to produce these first examples of *K*-means landscapes. The topography of the cost function surfaces underpins the performance of *K*-means clustering, and visualization of *K*-means landscapes is a crucial step to understanding its successes and failures. The landscape classification for different types of dataset and cluster number can be used to predict the performance of *K*-means for novel datasets and which variant of the algorithm may be most successful.

For all cluster numbers, the landscapes of Fisher's *Iris* dataset are broadly funneled toward a well-defined global minimum. The organization of these landscapes arises due to distinct partitionings of the *K* clusters among the two underlying clusters of the *Iris* dataset. Clear partitionings arise in this dataset because it contains well-defined clusters with no outliers. Partitionings, as a method of understanding the organization of *K*-means landscapes, will be present in other datasets, but may be less instructive as the number of outliers grows and the underlying clusters are less clearly separated. The clustered nature of the dataset leads to funneled landscapes. This feature of the *Iris* dataset is a characteristic of many other datasets, and we expect all datasets that have clearly defined underlying clusters to fall into the single-funnel paradigm at low values of *K*.

Single-funnel landscapes are associated with efficient global optimization, and the basin-hopping approach should work well.⁷⁵⁻⁷⁷ Many alternative schemes for *K*-means attempt to locate the global minimum through multiple initialization attempts, with no further perturbations to the located minima, and may also produce low-lying solutions for funneled landscapes. In general, we recommend basin-hopping, especially for more complex landscapes.

As *K* increases, there are more subfunnels in the landscapes that correspond to alternative partitionings. The number of such subfunnels will depend on the exact nature of the dataset, but it is clear that the number increases with *K*. The corresponding reduction in the cost function gap between the global minimum funnel and subfunnels, as *K* increases, leads to a more multifunneled character in the landscapes. It is expected that as *K* continues to increase, the landscapes will progress through multifunneled and perhaps to glassy character.

The growth in subfunnels leads to more challenging global optimization, which we have illustrated by an analog of kinetic analysis

for a molecular landscape. The effect of subfunnels in potentially impeding global optimization was quantified using transition rates from subfunnels to the global minimum funnel. The kinetic analysis shows that it is harder to escape from subfunnels and reach the global minimum as *K* increases, which is largely due to more transition state analogs on the pathways between partitionings, resulting from the increased number of minima. However, each region of the landscape associated with a particular partitioning remains locally funneled, at all *K*.

The distribution of minima among the funnels can change significantly, and we observe an increased proportion of minima in the lower-valued subfunnels as the number of clusters increases. At larger *K*, due to the higher proportion of low-valued minima, it is easier to obtain a clustering that is close in cost function to the global minimum. Therefore, for schemes with multiple initializations, there is an increased probability of obtaining a good solution, judged by the metric of the cost function. However, the growing number of different partitionings and the decreasing number of minima associated with the global minimum partitioning make it more likely that an alternative morphology to the global minimum will be found. For exploratory data analysis, the aim is not only to locate a low-valued minimum but also one with the same structure as the global minimum. Therefore, a much larger number of initializations may be required as *K* increases.

The current study provides another application of the energy landscape framework, which has recently been extended to understand the solution structure of neural networks⁷⁸⁻⁸⁰ and path integral representations.⁸¹ Here, we have used the same approach to analyze the performance of the *K*-means algorithm in terms of *K*-means landscapes. We hope that this topographical perspective can be used to guide future use of *K*-means and allow the understanding and tools developed in molecular science to be applied. Many questions remain about *K*-means landscapes, and we expect that further work addressing complex, real-world datasets will further our understanding of *K*-means clustering.

## SUPPLEMENTARY MATERIAL

See the supplementary material for a detailed description of the software used to generate all results presented in this work (S-I) and application of the methodology to two additional datasets (S-II). Section S-I provides a link to the Gitlab repository from which the software suite can be downloaded, tutorials for use of the software, and a description of its application to novel datasets. Section S-II illustrates *K*-means landscapes for two additional datasets with greater $N_f$ and *N* to provide evidence for the robustness of these algorithms.

## ACKNOWLEDGMENTS

L.D. acknowledges the EPSRC center for doctoral training in computational methods for materials science for funding under Grant No. EP/L015552/1. D.J.W. gratefully acknowledges financial support from the EPSRC. Additionally, D.J.W. gratefully acknowledges an International Chair at the Interdisciplinary Institute for Artificial Intelligence at 3iA Cote d'Azur, supported by the
French government, Reference No. ANR-19-P3IA-0002, which has
provided interactions that furthered the present research project.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available
within this article and its [supplementary material](#).

## REFERENCES

<footnote n="1">¹D. Xu and Y. Tian, "A comprehensive survey of clustering algorithms," *Ann. Data Sci.* **2**, 165–193 (2015).</footnote>

<footnote n="2">²T. F. Gonzalez, "Clustering to minimize the maximum intercluster distance," *Theor. Comput. Sci.* **38**, 293–306 (1985).</footnote>

<footnote n="3">³U. Sengupta, M. Carballo-Pacheco, and B. Strodel, "Automated Markov state models for molecular dynamics simulations of aggregation and self-assembly," *J. Chem. Phys.* **150**, 115101–115113 (2019).</footnote>

<footnote n="4">⁴J. A. Botía, J. Vandrovcova, P. Forabosco, S. Guelfi, K. D'Sa, J. Hardy, C. M. Lewis, M. Ryten, and M. E. Weale, "An additional *k*-means clustering step improves the biological features of WGCNA gene co-expression networks," *BMC Syst. Biol.* **11**, 47 (2017).</footnote>

<footnote n="5">⁵P. B. Frandsen, B. Calcott, C. Mayer, and R. Lanfear, "Automatic selection of partitioning schemes for phylogenetic analyses using iterative *k*-means clustering of site rates," *BMC Evol. Biol.* **15**, 13 (2015).</footnote>

<footnote n="6">⁶J. Paek and J. Ko, "*k*-means clustering-based data compression scheme for wireless imaging sensor networks," *IEEE Syst. J.* **11**, 2652–2662 (2017).</footnote>

<footnote n="7">⁷F. Cai, N.-A. Le-Khac, and M.-T. Kechadi, "Clustering approaches for financial data analysis: A survey," in *International Conference on Data Mining* (Springer-Verlag, 2016), pp. 105–111.</footnote>

<footnote n="8">⁸D. Aloise, A. Deshpande, P. Hansen, and P. Popat, "NP-hardness of Euclidean sum-of-squares clustering," *Mach. Learn.* **75**, 245–248 (2009).</footnote>

<footnote n="9">⁹M. Mahajan, P. Nimbhorkar, and K. Varadarajan, "The planar *k*-means problem is NP-hard," *Theor. Comput. Sci.* **442**, 13–21 (2012).</footnote>

<footnote n="10">¹⁰D. Steinley and M. J. Brusco, "Initializing *K*-means batch clustering: A critical evaluation of several techniques," *J. Classif.* **24**, 99–121 (2007).</footnote>

<footnote n="11">¹¹D. Arthur and S. Vassilvitskii, "*k*-means++: The advantages of careful seeding," in *Proceedings of the 18th Annual ACM-SIAM Symposium on Discrete Algorithms* (Society for Industrial and Applied Mathematics, 2007), pp. 1027–1035.</footnote>

<footnote n="12">¹²O. Bachem, M. Lucic, S. Hamed Hassani, and A. Krause, "Fast and provably good seedings for *k*-means," in *Proceedings of the 30th International Conference on Neural Information Processing Systems* (Curran Associates Inc., 2016), pp. 55–63.</footnote>

<footnote n="13">¹³D. Sculley, "Web-scale *k*-means clustering," in *Proceedings of the 19th International Conference on World Wide Web* (Association for Computing Machinery, 2010), pp. 1177–1178.</footnote>

<footnote n="14">¹⁴M. B. Cohen, S. Elder, C. Musco, C. Musco, and M. Persu, "Dimensionality reduction for *k*-means clustering and low rank approximation," in *Proceedings of the 47th Annual ACM Symposium on Theory of Computing* (Association for Computing Machinery, 2015), pp. 163–172.</footnote>

<footnote n="15">¹⁵M. Capó, A. Pérez, and J. A. Lozano, "An efficient approximation to the *K*-means clustering for massive data," *Knowl.-Based Syst.* **117**, 56–69 (2017).</footnote>

<footnote n="16">¹⁶L. Bottou and Y. Bengio, "Convergence properties of the *K*-means algorithms," in *International Conference on Neural Information Processing Systems* (MIT Press, 1994), pp. 585–592.</footnote>

<footnote n="17">¹⁷S. Bubeck, M. Meilä, and U. von Luxburg, "How the initialization affects the stability of the *k*-means algorithm," *ESAIM: Probab. Stat.* **16**, 436–452 (2012).</footnote>

<footnote n="18">¹⁸W. Qian, Y. Zhang, and Y. Chen, "Structures of spurious local minima in *k*-means," *IEEE Trans. Inf. Theory* **68**, 395 (2021).</footnote>

<footnote n="19">¹⁹P. J. G. Lisboa, T. A. Etchells, I. H. Jarman, and S. J. Chambers, "Finding reproducible cluster partitions for the *k*-means algorithm," *BMC Bioinf.* **14**, S8 (2013).</footnote>

<footnote n="20">²⁰D. J. Wales, *Energy Landscapes* (Cambridge University Press, Cambridge, 2003).</footnote>

<footnote n="21">²¹R. A. Fisher, "The use of multiple measurements in taxonomic problems," *Ann. Eugen.* **7**, 179–188 (1936).</footnote>

<footnote n="22">²²J. N. Murrell and K. J. Laidler, "Symmetries of activated complexes," *Trans. Faraday Soc.* **64**, 371–377 (1968).</footnote>

<footnote n="23">²³S. Lloyd, "Least squares quantization in PCM," *IEEE Trans. Inf. Theory* **28**, 129–137 (1982).</footnote>

<footnote n="24">²⁴L. J. Munro and D. J. Wales, "Defect migration in crystalline silicon," *Phys. Rev. B* **59**, 3969–3980 (1999).</footnote>

<footnote n="25">²⁵G. Henkelman and H. Jónsson, "A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives," *J. Chem. Phys.* **111**, 7010–7022 (1999).</footnote>

<footnote n="26">²⁶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).</footnote>

<footnote n="27">²⁷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).</footnote>

<footnote n="28">²⁸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).</footnote>

<footnote n="29">²⁹C. G. Broyden, "The convergence of a class of double-rank minimization algorithms I. General considerations," *IMA J. Appl. Math.* **6**, 76–90 (1970).</footnote>

<footnote n="30">³⁰R. Fletcher, "A new approach to variable metric algorithms," *Comput. J.* **13**, 317–322 (1970).</footnote>

<footnote n="31">³¹D. Goldfarb, "A family of variable-metric methods derived by variational means," *Math. Comput.* **24**, 23–26 (1970).</footnote>

<footnote n="32">³²D. F. Shanno, "Conditioning of quasi-Newton methods for function minimization," *Math. Comput.* **24**, 647–656 (1970).</footnote>

<footnote n="33">³³J. Nocedal, "Updating quasi-Newton matrices with limited storage," *Math. Comput.* **35**, 773–782 (1980).</footnote>

<footnote n="34">³⁴D. C. Liu and J. Nocedal, "On the limited memory BFGS method for large scale optimization," *Math. Program.* **45**, 503–528 (1989).</footnote>

<footnote n="35">³⁵D. J. Wales, "Energy landscapes: Calculating pathways and rates," *Int. Rev. Phys. Chem.* **25**, 237–282 (2006).</footnote>

<footnote n="36">³⁶F. Noé and S. Fischer, "Transition networks for modeling the kinetics of conformational change in macromolecules," *Curr. Opin. Struct. Biol.* **18**, 154–162 (2008).</footnote>

<footnote n="37">³⁷D. J. Wales, "Calculating rate constants and committor probabilities for transition networks by graph transformation," *J. Chem. Phys.* **130**, 204111 (2009).</footnote>

<footnote n="38">³⁸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).</footnote>

<footnote n="39">³⁹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).</footnote>

<footnote n="40">⁴⁰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).</footnote>

<footnote n="41">⁴¹D. J. Wales, M. A. Miller, and T. R. Walsh, "Archetypal energy landscapes," *Nature* **394**, 758–760 (1998).</footnote>

<footnote n="42">⁴²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).</footnote>

<footnote n="43">⁴³J. N. Onuchic, N. D. Socci, Z. Luthey-Schulten, and P. G. Wolynes, "Protein folding funnels: The nature of the transition state ensemble," *Folding Des.* **1**, 441–450 (1996).</footnote>

<footnote n="44">⁴⁴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).</footnote>
<footnote n="45">⁴⁵K. A. Dill, D. O. V. Alonso, and K. Hutchinson, "Thermal stabilities of globular proteins," *Biochemistry* **28**, 5439–5449 (1989).</footnote>

<footnote n="46">⁴⁶P. G. Wolynes, Z. Luthey-Schulten, and J. N. Onuchic, "Fast-folding experiments and the topography of protein folding energy landscapes," *Chem. Biol.* **3**, 425–432 (1996).</footnote>

<footnote n="47">⁴⁷C. M. Dobson, "Protein folding and misfolding," *Nature* **426**, 884–890 (2003).</footnote>

<footnote n="48">⁴⁸S. Z. Selim and K. Alsultan, "A simulated annealing algorithm for the clustering problem," *Pattern Recognit.* **24**, 1003–1008 (1991).</footnote>

<footnote n="49">⁴⁹G. P. Babu and M. N. Murty, "A near-optimal initial seed value selection in K-means algorithm using a genetic algorithm," *Pattern Recognit. Lett.* **14**, 763–769 (1993).</footnote>

<footnote n="50">⁵⁰G. P. Babu and M. N. Murty, "Simulated annealing for selecting optimal initial seeds in the K-means algorithm," *Ind. J. Pure Appl. Math.* **25**, 85–94 (1994).</footnote>

<footnote n="51">⁵¹K. Krishna and M. N. Murty, "Genetic K-means algorithm," *IEEE Trans. Syst. Man Cybern.* **29**, 433–439 (1999).</footnote>

<footnote n="52">⁵²P. Fränti, "Genetic algorithm with deterministic crossover for vector quantization," *Pattern Recognit. Lett.* **21**, 61–68 (2000).</footnote>

<footnote n="53">⁵³T. Niknam and B. Amiri, "An efficient hybrid approach based on PSO, ACO and k-means for cluster analysis," *Appl. Soft Comput.* **10**, 183–197 (2010).</footnote>

<footnote n="54">⁵⁴F. Sciortino, W. Kob, and P. Tartaglia, "Thermodynamics of supercooled liquids in the inherent-structure formalism: A case study," *J. Phys.: Condens. Matter* **12**, 6525–6534 (2000).</footnote>

<footnote n="55">⁵⁵T. V. Bogdan, D. J. Wales, and F. Calvo, "Equilibrium thermodynamics from basin-sampling," *J. Chem. Phys.* **124**, 044102 (2006).</footnote>

<footnote n="56">⁵⁶D. J. Wales, "Surveying a complex potential energy landscape: Overcoming broken ergodicity using basin-sampling," *Chem. Phys. Lett.* **584**, 1–9 (2013).</footnote>

<footnote n="57">⁵⁷N. D. Socci, J. N. Onuchic, and P. G. Wolynes, "Protein folding mechanisms and the multidimensional folding funnel," *Proteins: Struct., Funct., Bioinf.* **32**, 136–158 (1998).</footnote>

<footnote n="58">⁵⁸J. P. K. Doye, M. A. Miller, and D. J. Wales, "The double-funnel energy landscape of the 38-atom Lennard-Jones cluster," *J. Chem. Phys.* **110**, 6896–6906 (1999).</footnote>

<footnote n="59">⁵⁹Y. Chebaro, A. J. Ballard, D. Chakraborty, and D. J. Wales, "Intrinsically disordered energy landscapes," *Sci. Rep.* **5**, 10386 (2015).</footnote>

<footnote n="60">⁶⁰K. Röder, J. A. Joseph, B. E. Husic, and D. J. Wales, "Energy landscapes for proteins: From single funnels to multifunctional systems," *Adv. Theory Simul.* **2**, 1800175–1800199 (2019).</footnote>

<footnote n="61">⁶¹V. K. de Souza and D. J. Wales, "Energy landscapes for diffusion: Analysis of cage-breaking processes," *J. Chem. Phys.* **129**, 164507 (2008).</footnote>

<footnote n="62">⁶²J. P. K. Doye and D. J. Wales, "On potential energy surfaces and relaxation to the global minimum," *J. Chem. Phys.* **105**, 8428–8445 (1996).</footnote>

<footnote n="63">⁶³D. J. Wales, "Decoding the energy landscape: Extracting structure, dynamics and thermodynamics," *Philos. Trans. R. Soc. London, Ser. A* **370**, 2877–2899 (2012).</footnote>

<footnote n="64">⁶⁴S. P. Niblett, V. K. de Souza, J. D. Stevenson, and D. J. Wales, "Dynamics of a molecular glass former: Energy landscapes for diffusion in ortho-terphenyl," *J. Chem. Phys.* **145**, 024505 (2016).</footnote>

<footnote n="65">⁶⁵V. K. de Souza and D. J. Wales, "The potential energy landscape for crystallisation of a Lennard-Jones fluid," *J. Stat. Mech.: Theory Exp.* **2016**, 074001.</footnote>

<footnote n="66">⁶⁶E. W. Dijkstra, "A note on two problems in connexion with graphs," *Numer. Math.* **1**, 269–271 (1959).</footnote>

<footnote n="67">⁶⁷S. A. Trygubenko and D. J. Wales, "Kinetic analysis of discrete path sampling stationary point databases," *Mol. Phys.* **104**, 1497–1507 (2006).</footnote>

<footnote n="68">⁶⁸S. A. Trygubenko and D. J. Wales, "Graph transformation method for calculating waiting times in Markov chains," *J. Chem. Phys.* **124**, 234110 (2006).</footnote>

<footnote n="69">⁶⁹W. Forst, *Theory of Unimolecular Reactions* (Academic Press, New York, 1973).</footnote>

<footnote n="70">⁷⁰K. J. Laidler. *Chemical Kinetics*. Harper & Row, New York, 1987.</footnote>

<footnote n="71">⁷¹G. S. Hammond, "A correlation of reaction rates," *J. Am. Chem. Soc.* **77**, 334–338 (1955).</footnote>

<footnote n="72">⁷²D. J. Wales, "A microscopic basis for the global appearance of energy landscapes," *Science* **293**, 2067–2069 (2001).</footnote>

<footnote n="73">⁷³R. Gilmore, *Catastrophe Theory for Scientists and Engineers* (Wiley, New York, 1981).</footnote>

<footnote n="74">⁷⁴T. V. Bogdan and D. J. Wales, "New results for phase transitions from catastrophe theory," *J. Chem. Phys.* **120**, 11090–11099 (2004).</footnote>

<footnote n="75">⁷⁵Z. Li and H. A. Scheraga, "Monte Carlo-minimization approach to the multiple-minima problem in protein folding," *Proc. Natl. Acad. Sci. U. S. A.* **84**, 6611–6615 (1987).</footnote>

<footnote n="76">⁷⁶D. J. Wales and J. P. K. Doye, "Global optimization by basin-hopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms," *J. Phys. Chem. A* **101**, 5111–5116 (1997).</footnote>

<footnote n="77">⁷⁷D. J. Wales and H. A. Scheraga, "Global optimization of clusters, crystals and biomolecules," *Science* **285**, 1368–1372 (1999).</footnote>

<footnote n="78">⁷⁸A. J. Ballard, J. D. Stevenson, R. Das, and D. J. Wales, "Energy landscapes for a machine learning application to series data," *J. Chem. Phys.* **144**, 124119 (2016).</footnote>

<footnote n="79">⁷⁹R. Das and D. J. Wales, "Machine learning prediction for classification of outcomes in local minimisation," *Chem. Phys. Lett.* **667**, 158–164 (2017).</footnote>

<footnote n="80">⁸⁰A. J. Ballard, R. Das, S. Martiniani, D. Mehta, L. Sagun, J. D. Stevenson, and D. J. Wales, "Energy landscapes for machine learning," *Phys. Chem. Chem. Phys.* **19**, 12585–12603 (2017).</footnote>

<footnote n="81">⁸¹C. L. Vaillant, S. C. Althorpe, and D. J. Wales, "Path integral energy landscapes for water clusters," *J. Chem. Theory Comput.* **15**, 33–42 (2018).</footnote>