# Convex Clustering through MM: An Efficient Algorithm to Perform Hierarchical
Clustering

Daniel J.W. Touw<EMAIL_ADDRESS>
Patrick J.F. Groenen<EMAIL_ADDRESS>
Econometric Institute
Erasmus University Rotterdam
P.O. Box 1738
3000 DR Rotterdam
The Netherlands Yoshikazu Terada<EMAIL_ADDRESS>
Graduate School of Engineering Science
Osaka University
1-3 Machikaneyama-cho, Toyonaka
Osaka 560-8531, Japan

###### Abstract

Convex clustering is a modern method with both hierarchical and $k$-means
clustering characteristics. Although convex clustering can capture the complex
clustering structure hidden in data, the existing convex clustering algorithms
are not scalable to large data sets with sample sizes greater than ten
thousand. Moreover, it is known that convex clustering sometimes fails to
produce hierarchical clustering structures. This undesirable phenomenon is
called cluster split and makes it difficult to interpret clustering results.
In this paper, we propose convex clustering through majorization-minimization
(CCMM)—an iterative algorithm that uses cluster fusions and sparsity to
enforce a complete cluster hierarchy with reduced memory usage. In the CCMM
algorithm, the diagonal majorization technique makes a highly efficient update
for each iteration. With a current desktop computer, the CCMM algorithm can
solve a single clustering problem featuring over one million objects in seven-
dimensional space within 70 seconds.

Keywords: convex clustering, convex optimization, hierarchical clustering,
majorization-minimization algorithm, unsupervised learning

##  1 Introduction

Clustering is a core method in unsupervised machine learning and is used in
many fields. Applications range from customer segmentation in marketing
analytics to object detection in image processing. Over the past two decades,
a new method called convex clustering has been introduced that combines
aspects of the two popular techniques: $k$-means (MacQueen, 1967) and
hierarchical clustering (see, for example, Gan et al., 2007).

The convex clustering framework was introduced by Pelckmans et al. (2005) who
formulated clustering as shrinkage of centroid distances subject to a convex
loss function where each $p$-dimensional object $\mathbf{x}_{i}$ is
represented by its own centroid $\mathbf{a}_{i}$. Distance shrinkage is
applied to the pairs of $\mathbf{a}_{i}$ and $\mathbf{a}_{j}$ using a penalty
term, and clusters emerge via the occurrence of equal cluster centroids.
Hocking et al. (2011) have explained convex clustering as a convex relaxation
of hierarchical clustering, while Lindsten et al. (2011) have staged it as a
convex relaxation of $k$-means clustering. The common loss function that is
considered in these studies is

$L(\mathbf{A})=\frac{1}{2}\|\mathbf{X-A}\|^{2}+\lambda\sum_{i<j}w_{ij}\|\mathbf{a}_{i}-\mathbf{a}_{j}\|,$
(1)

where $\mathbf{A}$ is a $n\times p$ representation matrix whose rows
$\mathbf{a}_{i}^{\top}$ are centroids corresponding to objects, $w_{ij}$ is a
user-defined nonnegative weight that reflects the importance of clustering
objects $i$ and $j$, $\mathbf{X}$ is the data matrix with rows
$\mathbf{x}_{i}^{\top}$, and $\lambda$ is a penalty-strength parameter that
controls the shrinkage of cluster-centroids. Our notation $\|\cdot\|$ denotes
either the $\ell_{2}$-norm of a vector or the Frobenius norm of a matrix. The
second term in (1) consists of $\|\mathbf{a}_{i}-\mathbf{a}_{j}\|$
representing the Euclidean distance between rows $i$ and $j$ of $\mathbf{A}$.
Each distance is a fusion penalty and behaves like a grouped Lasso (Yuan and
Lin, 2006). The minimization of the loss function causes some of these
distances $\|\mathbf{a}_{i}-\mathbf{a}_{j}\|$ to shrink towards zero, thus
interpreted as grouping objects $i$ and $j$. In related literature, this
penalty is often defined using the more general $\ell_{q}$-norm. In this
study, we focus on the $\ell_{2}$-norm since it guarantees a single
partitioning of clusters for a given $\lambda$.

Without loss of generality, we now assume each column of $\mathbf{X}$ has a
mean zero, and we propose a normalized version of (1), that is,

$L(\mathbf{A})=\kappa_{\epsilon}\|\mathbf{X-A}\|^{2}+\lambda\kappa_{\text{pen}}\sum_{i<j}w_{ij}\|\mathbf{a}_{i}-\mathbf{a}_{j}\|,$
(2)

where $\kappa_{\epsilon}=(2\|\mathbf{X}\|^{2})^{-1}$ and
$\kappa_{\text{pen}}=(\|\mathbf{X}\|\sum_{i<j}w_{ij})^{-1}$. Normalization
constants $\kappa_{\epsilon}$ and $\kappa_{\text{pen}}$ are chosen so that the
clustering of the objects in $\mathbf{X}$, given a particular value of
$\lambda$, proves invariant to $n$ and $p$, as well as to the scales of
$w_{ij}$ and $\mathbf{X}$.

One vital property of convex clustering is its strong convexity of the loss
function in $\mathbf{A}$ (see the proof in Appendix A). In contrast to
$k$-means, the minimum solution of (2) with a given $\lambda\geq 0$ proves
unique, and convex clustering here averts the local minima issue.
Additionally, the minimization of the loss function for a sequence of values
for $\lambda$ yields a series of solutions $\mathbf{A}_{\lambda}$ that
reflects the movement of cluster centroids through the $p$-dimensional space.
These paths are commonly referred to as the clusterpath and allow for dynamic
visualizations of results (Weylandt et al., 2020). Figure 1a shows a
clusterpath that results from minimizing (2) for a data set featuring seven
objects in $\mathbb{R}^{2}$. In Figure 1b, these results are converted into a
dendrogram—a common visual representation in hierarchical clustering.

(a) Clusterpath

(b) Clusterpath dendrogram

Figure 1: In Panel (a), a clusterpath is computed for seven data points in
$\mathbb{R}^{2}$ by minimizing (2) for a sequence of values for
$\lambda\in[0,39]$. The objects in the $7\times 2$ data matrix $\mathbf{X}$
are denoted by the blue dots with the paths traced by their cluster centroids
in $\mathbf{A}$ appearing in grey. For $\lambda=39$, the paths of each cluster
centroid join at the central black dot. In Panel (b), the corresponding
dendrogram appears using the square root of $\lambda$ as the height where
objects are clustered.

Convex clustering incurs several shortcomings. First, Hocking et al. (2011)
have noted for particular choices of $w_{ij}$ in (2) that centroids equal for
some value $\lambda_{1}$ can reseparate for some value
$\lambda_{2}>\lambda_{1}$. This potential violation of cluster hierarchy is
undesirable from an interpretative point of view and should be avoided.
Second, apart from special focus on the weights, the shapes of clusters
recovered by the minimization of (2) tend to be spherical, which is a limiting
factor in the application of the method. Third, existing algorithms scale
poorly for data sets exceeding a few thousand objects. This study tackles
these issues: we develop an algorithm that enforces cluster hierarchy, allows
nonconvex cluster-shape retrieval, and is computationally at least 200 times
more efficient than state-of-the-art implementations.

The rest of this paper is organized as follows. Section 2 details three
different approaches for choosing $w_{ij}$ in (2). We focus on the third
(final) choice, sparse weights, to further the understanding of the effect
they have on the clusterpath. In Section 3, we combine sparsity and cluster
fusions into a majorization-minimization (MM) algorithm to minimize the loss
function. Such an iterative algorithm, where a simple optimization problem
proxies a more complex problem, aptly converges toward the optimum of the
complex problem (Hunter and Lange, 2004). We dub this resulting algorithm
convex clustering through majorization-minimization (CCMM). In Sections 4 and
5, we present results of the simulation study and conclude our paper.

##  2 Hierarchy and Weights in Convex Clustering

One potential pitfall of convex clustering is that increasing $\lambda$ does
not always guarantee a complete cluster hierarchy. In this section, we offer a
solution for this. To make convex clustering scalable for big $n$, sparsity of
weights is necessary where many $w_{ij}$ are set to zero. Sun et al. (2021)
have shown that using such sparsity in weights equips their method to minimize
a single instance of (2) in roughly six minutes for a data set of 200,000
objects. However, under such extreme sparsity, it is likely that some groups
of objects are not connected to the remaining objects via nonzero weights.
Consequently, the original convex clustering problem can be reduced to
separate optimization problems. It is this reducibility that renders a
complete cluster hierarchy impossible. Here, we provide a formal proof of this
property and craft a solution that has minimal impact on the scalability of
our algorithm.

###  2.1 Enforcing Hierarchical Clusters

To address cluster splits, we first note that Hocking et al. (2011)
conjectured that splits do not occur when weights $w_{ij}$ are nonincreasing
in $\|\mathbf{x}_{i}-\mathbf{x}_{j}\|$. This has been empirically confirmed by
Chi and Lange (2015) and Weylandt et al. (2020), as well as in our research.
Therefore, if the weights satisfy this condition, we can fuse the rows of
$\mathbf{A}$ where $\|\mathbf{a}_{i}-\mathbf{a}_{j}\|=0$ and store all unique
centroids in $c\times p$ matrix $\mathbf{M}$. We next define $n\times c$
cluster membership matrix $\mathbf{U}$ as having elements

$u_{ik}=\begin{cases}1&\text{if object $i$ belongs to cluster $k$, and}\\\
0&\text{otherwise},\end{cases}$

and denote row $i$ of $\mathbf{U}$ by row vector $\mathbf{u}_{i}^{\top}$ and
column $k$ by column vector $\mathbf{u}_{k}$. Instead of minimizing (2) with
respect to $\mathbf{A}$, we perform minimization on $\mathbf{M}$ and continue
the fusion process whenever $\|\mathbf{m}_{k}-\mathbf{m}_{l}\|=0$. The matrix
$\mathbf{A}$ can be retrieved simply via $\mathbf{A=UM}$.

We propose preventing any merged clusters from splitting, thereby ensuring a
hierarchy of clusters with increasing $\lambda$. This strategy also yields a
substantial gain in computational efficiency since $\mathbf{M}$ has fewer rows
than $\mathbf{A}$.

###  2.2 Introducing Sparsity

In existing literature, three general approaches for choosing the weights can
be identified. The first type, used by Pelckmans et al. (2005) and Lindsten et
al. (2011), simply sets all $w_{ij}=1$ in a unit-weights approach. However,
using convex clustering in combination with unit weights leads to quite poor
clustering performance. This is illustrated by the clusterpath in Figure 2a
where almost no clustering takes place before each individual path reaches the
center. Hocking et al. (2011) have found that weights that decrease in the
distance between $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ allow for better
cluster recovery. They suggest using weights proportional to the Gaussian
distribution: $w_{ij}=\exp(-\psi\|\mathbf{x}_{i}-\mathbf{x}_{j}\|^{2})$ with
hyperparameter $\psi$ scaling these high-dimensional Euclidean distances.

Figure 2b shows that Gaussian weights indeed outperform unit weights in terms
of clustering, but fail in allowing recovery of the original two interlocking
half moons. Finally, as an extension to Gaussian weights, Chi and Lange (2015)
set $w_{ij}$ to zero when object $i$ is not among the $k$ nearest neighbors of
$j$. This approach improves clustering performance (see Figure 2c) with the
added advantage of greatly reducing the computational complexity of
minimization algorithms. In this paper, a scaled version of the $k$-nn
Gaussian weights is used to ensure that weights remain unaffected by the scale
of the data, yielding the following representation of the sparse weight matrix
$\mathbf{W}$

$w_{ij}=w_{ji}=\begin{cases}\exp\left(-\phi\frac{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|^{2}}{\operatorname*{mean}_{i^{\prime},j^{\prime}}\|\mathbf{x}_{i^{\prime}}-\mathbf{x}_{j^{\prime}}\|^{2}}\right)&\text{if
$j\in\mathcal{S}_{i}^{(k)}\lor i\in\mathcal{S}_{j}^{(k)}$}\\\
0&\text{otherwise},\end{cases}$ (3)

where $\mathcal{S}_{i}^{(k)}$ denotes the set of $k$ nearest neighbors of $i$,
$\phi\geq 0$ is a tuning parameter, and
$\|\mathbf{x}_{i}-\mathbf{x}_{j}\|^{2}$ is scaled by the mean of all squared
distances.

(a) Unit weights

(b) Gaussian weights

(c) $k$-nn Gaussian weights

Figure 2: Comparison of three clusterpaths computed for the same data using
unit weights in the penalty for Panel (a), Gaussian weights with $\phi=4$ for
Panel (b), and $k$-nn Gaussian weights with $k=10$ and $\phi=4$ for Panel (c).
A black dot marks the location of the final cluster in each clusterpath.

###  2.3 Enforcing Irreducibility

Weight sparsity can improve the flexibility of convex clustering and slash
computation time. However, sparsity also incurs a new challenge: the emergence
of different groups of objects, for example $\mathcal{C}_{k}$ and
$\mathcal{C}_{l}$ where $w_{ij}=0$ for all $i\in\mathcal{C}_{k}$ and
$j\in\mathcal{C}_{l}$, that trigger reducibility of the convex clustering
problem where a loss function can split into standalone sub-problems. Such a
$\mathbf{W}$ we call “disconnected” per the corresponding disconnected graph
in graph theory. The theorem below presents irreducibility (connectedness) of
$\mathbf{W}$ as a sufficient condition for a clusterpath ending in a single
partition.

###### Theorem 1

If $\mathbf{W}$ is connected, then for some sufficiently large, finite value
of $\lambda$, a global minimum of (2) exists where all rows of the solution
$\mathbf{A}$ are identical to the mean of the rows in $\mathbf{X}$.

This theorem is a more general version of Proposition 1 in Pelckmans et al.
(2005), and we provide the proof in Appendix B. Using Theorem 1, we can
determine what happens if $\mathbf{W}$ is disconnected as formalized in the
following corollary.

###### Corollary 2

Let $\mathbf{A}_{\mathcal{C}_{k}}$ be the $|\mathcal{C}_{k}|\times p$ matrix
that results from selecting only rows $i\in\mathcal{C}_{k}$ from $\mathbf{A}$,
and let $\mathbf{W}_{\mathcal{C}_{k}}$ be the
$|\mathcal{C}_{k}|\times|\mathcal{C}_{k}|$ matrix that results from selecting
only rows $i\in\mathcal{C}_{k}$ and columns $j\in\mathcal{C}_{k}$ from
$\mathbf{W}$ where $|\mathcal{C}_{k}|=\bm{1}^{\top}\mathbf{u}_{k}$ denotes the
size of cluster $k$. If $\mathbf{W}$ is disconnected, the convex clustering
problem can disaggregate into $K$ separate clustering problems where a
sufficiently large, finite value for $\lambda$ minimizes each problem via
$\mathbf{A}_{\mathcal{C}_{k}}=\bm{1}\bar{\mathbf{x}}_{k}^{\top}$. Here,
$\mathcal{C}_{k}$ comprise the $K$ subsets of objects where
$\mathbf{W}_{\mathcal{C}_{k}}$ is connected and
$\mathbf{W}_{\mathcal{C}_{k}\cup\mathcal{C}_{l}}$ is disconnected for any
$k\neq l$, and $\bar{\mathbf{x}}_{k}$ is the mean of $\mathbf{x}_{i}$ for
$i\in\mathcal{C}_{k}$.

Proof  For disconnected $\mathbf{W}$, the loss function $L(\mathbf{A})$
equates the sum of $K$ distinct loss functions:

$\displaystyle L(\mathbf{A})$
$\displaystyle=\kappa_{\epsilon}\sum_{i=1}^{n}\|\mathbf{x}_{i}-\mathbf{a}_{i}\|^{2}+\lambda\kappa_{\text{pen}}\sum_{i<j}w_{ij}\|\mathbf{a}_{i}-\mathbf{a}_{j}\|$
$\displaystyle=\sum_{k=1}^{K}\left[\kappa_{\epsilon}\sum_{i\in\mathcal{C}_{k}}\|\mathbf{x}_{i}-\mathbf{a}_{i}\|^{2}+\lambda\kappa_{\text{pen}}\sum_{\begin{subarray}{c}i,j\in\mathcal{C}_{k}\\\
i<j\end{subarray}}w_{ij}\|\mathbf{a}_{i}-\mathbf{a}_{j}\|\right].$

The minimum of this sum is obtained by independently minimizing each sub-
problem. Following Theorem 1, the rows of $\mathbf{A}_{\mathcal{C}_{k}}$ can
be computed as
$\mathbf{a}^{(k)}=|\mathcal{C}_{k}|^{-1}\sum_{i\in\mathcal{C}_{k}}\mathbf{x}_{i}=\mathbf{\bar{x}}_{k}$
for sufficiently large $\lambda$.
Since the sparse specification of $\mathbf{W}$ in (3) does not guarantee a
connected weight matrix, Theorem 1 and Corollary 2 indicate that it may be
impossible to obtain a single cluster when performing convex clustering for
large $\lambda$. Still, hierarchical clustering methods require all objects to
eventually merge into a single cluster. To attain a single-cluster result from
convex clustering using a sparse weight matrix, connectedness of $\mathbf{W}$
must be guaranteed.

(a) Weights without symmetric circulant

(b) Weights with symmetric circulant

(c) Clusterpath without symmetric circulant

(d) Clusterpath with symmetric circulant

Figure 3: The effect of the use of a symmetric circulant matrix to ensure
connectedness of the weight matrix on convex clustering. In Panel (a), two
objects are connected by an edge if their weight as determined by (3) is
nonzero and two objects are connected by an edge in Panel (b) if their weight
as determined by (4) is nonzero. In both cases $k=3$ is used. The additional
edges that are added due to the symmetric circulant matrix are indicated by
light grey dashed lines. The clusterpaths that result from convex clustering
with the two different weight matrices are shown in panels (c) and (d). The
colors indicate the true clusters according to the data generating process and
the black dot the location of the final cluster(s) in each clusterpath.

To do so, we propose adding nonzero weights to the sparse weight matrix
$\mathbf{W}$ to ensure its connectedness. To avoid expensive pre-processing
steps, we enlist a symmetric circulant matrix to provide a rule for adding
non-zero entries to $\mathbf{W}$. A symmetric circulant is a square matrix
$\mathbf{H}$ where each row is identical to the row above, but shifted one
position to the right and satisfying $h_{ij}=h_{ji}$ (Gower and Groenen,
1991). One example of the most sparse symmetric circulant not forming the
identity matrix for $n=6$ is

$\mathbf{H}=\begin{bmatrix}0&1&0&0&0&1\\\ 1&0&1&0&0&0\\\ 0&1&0&1&0&0\\\
0&0&1&0&1&0\\\ 0&0&0&1&0&1\\\ 1&0&0&0&1&0\\\ \end{bmatrix}.$

As a weight matrix, $\mathbf{H}$ is connected by definition. Therefore, if
nonzero weights are added to $\mathbf{W}$ according to the nonzero elements of
$\mathbf{H}$, connectedness of $\mathbf{W}$ is guaranteed while adding no more
than $n$ nonsparse elements. In a graphical context, Figure 3a features
nonzero weights based on three nearest neighbors as edges between objects in a
generated data set. In Figure 3b, the additional nonzero weights per the
symmetric circulant matrix appear in light grey. The new formulation of the
weight matrix is

$w_{ij}=w_{ji}=\begin{cases}\exp\left(-\phi\frac{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|^{2}}{\operatorname*{mean}_{i^{\prime},j^{\prime}}(\|\mathbf{x}_{i^{\prime}}-\mathbf{x}_{j^{\prime}}\|^{2})}\right)&\text{if
$j\in\mathcal{S}_{i}^{(k)}\lor i\in\mathcal{S}_{j}^{(k)}\lor h_{ij}=1$},\\\
0&\text{otherwise}.\end{cases}$ (4)

Figure 3 also depicts the effect of $\mathbf{W}$ on the clusterpaths that
result from minimizing (2). In Figure 3c, the disconnected $\mathbf{W}$ from
(3), using only $k$-nn, fails to produce a clusterpath leading to a complete
hierarchy. In contrast, the combination of $k$-nn and a symmetric circulant in
(4) produces a $\mathbf{W}$ which nicely produces a complete hierarchy (see
Figure 3d).

The exact level of sparsity depends on the data. The $k$-nn relation is not
symmetric, and enforcement of symmetry in (3) and (4) can vary the number of
nonzero weights. An upper bound for the number of non-sparse elements
following (3) is $kn$, then increased to $(k+1)n$ when ensuring connectedness
via a symmetric circulant matrix in (4). This contrasts with the dense
specification of $\mathbf{W}$ always demanding the storage of $n(n-1)/2$
weights that proves quadratic in $n$.

##  3 Minimization Algorithm

Cluster fusions play an important role in the CCMM algorithm by ensuring that
no cluster splits occur while reducing the computational burden. Therefore, we
first show in Section 3.1 how fusions are enlisted to formulate the loss
function in terms of the unique centroids in matrix $\mathbf{A}$. In Section
3.2, a surrogate function based on the reformulated loss function is derived
that can be minimized analytically in linear time—making it an excellent
candidate for an iterative minimization algorithm. In Section 3.3, we present
the complete CCMM algorithm and discuss some of its properties.

###  3.1 Cluster Fusions

To avoid cluster splits and to reduce the computational burden, rows of
$\mathbf{A}$ that are identical can be merged as proposed by Hocking et al.
(2011). This approach tracks all unique rows of $\mathbf{A}$ and stores them
in the $c\times p$ matrix $\mathbf{M}$ where $c$ denotes the number of unique
rows or, equivalently, the current number of clusters. The decision to fuse
two clusters is based on whether $\mathbf{m}_{k}$ and $\mathbf{m}_{l}$ are
sufficiently close given some threshold $\varepsilon_{f}$. If
$\|\mathbf{m}_{l}-\mathbf{m}_{k}\|\leq\varepsilon_{f}$ for some $k\neq l$,
then rows $\mathbf{m}_{k}$ and $\mathbf{m}_{l}$ in $\mathbf{M}$ are replaced
by the weighted average

$\mathbf{m}_{new}=\frac{|C_{k}|\mathbf{m}_{k}+|C_{l}|\mathbf{m}_{l}}{|C_{k}|+|C_{l}|}.$
(5)

If $\mathbf{UM}$ replaces $\mathbf{A}$, we can rewrite the loss function in
(2) in terms of $\mathbf{M}$. As stated earlier, the advantage here is that
the complexity of computing an update for $\mathbf{M}$ in an iterative
algorithm no longer scales in the number of objects $n$, but in the number of
clusters $c$ ranging from $n$ to 1. The new loss function becomes

$\displaystyle L(\mathbf{M})$
$\displaystyle=\kappa_{\epsilon}\|\mathbf{X-UM}\|^{2}+\lambda\kappa_{\text{pen}}\sum_{i<j}w_{ij}\|\mathbf{u}_{i}^{\top}\mathbf{M}-\mathbf{u}_{j}^{\top}\mathbf{M}\|$
$\displaystyle=\kappa_{\epsilon}\|\mathbf{X-UM}\|^{2}+\lambda\kappa_{\text{pen}}\sum_{k<l}\mathbf{u}_{k}^{\top}\mathbf{Wu}_{l}\|\mathbf{m}_{k}-\mathbf{m}_{l}\|.$
(6)

However, since $\varepsilon_{f}$ is a value near zero, $\mathbf{m}_{k}$ and
$\mathbf{m}_{l}$ may not equate, thus yielding a difference between
$L(\mathbf{M})$ and the $L(\mathbf{M}_{new})$ that results from fusing
$\mathbf{m}_{k}$ and $\mathbf{m}_{l}$. The following theorem states that as
$\varepsilon_{f}$ approaches zero, so does
$|L(\mathbf{M})-L(\mathbf{M}_{new})|$. The proof is in Appendix C.

###### Theorem 3

Let the matrix $\mathbf{M}$ contain the $c$ unique rows of $\mathbf{A}$, and
let $\varepsilon_{f}$ be the threshold used for fusing rows of $\mathbf{M}$.
Assume that $\|\mathbf{m}_{k}-\mathbf{m}_{l}\|\leq\varepsilon_{f}$ for some
$k\neq l$. Let $\mathbf{m}_{new}$ be the weighted average of $\mathbf{m}_{k}$
and $\mathbf{m}_{l}$ computed as in (5), and let $\mathbf{M}_{new}$ be the
matrix that results from setting $\mathbf{m}_{k}$ and $\mathbf{m}_{l}$ to
$\mathbf{m}_{new}$. Then, for fixed $\lambda$ and $\mathbf{W}$, the absolute
error $|L(\mathbf{M})-L(\mathbf{M}_{new})|$ tends toward zero as
$\varepsilon_{f}\rightarrow 0$.

An important consequence of Theorem 3 is that for sufficiently small
$\varepsilon_{f}$ there are no order effects of fusing two rows of
$\mathbf{M}$.

###  3.2 Majorization-Minimization Update

To minimize the loss function in (6), we enlist MM. Here, an objective
function difficult to minimize is replaced by a simpler or surrogate function
with a minimum that can be computed analytically. A similar procedure was
first described by Weiszfeld (1937), but Ortega and Rheinboldt (1970) were the
first to apply it in the context of a line search. Later, it was developed as
a general minimization method by independently by De Leeuw () and Voß and
Eckhardt (1980). In machine learning literature, MM is also known as the
concave-convex procedure (CCCP; Yuille and Rangarajan, 2003). Additional
examples can be found in Lange et al. (2000) and Hunter and Lange (2004).

MM exploits a majorization function that touches the target function at a
“supporting” point and that is never exceeded by the target function over its
domain. Let $g(\mathbf{M},\mathbf{M}_{0})$ denote this majorization function
with $\mathbf{M}_{0}$ as the supporting point. The criteria can be summarized
as

$g(\mathbf{M},\mathbf{M}_{0})\geq L(\mathbf{M})\quad\text{and}\quad
g(\mathbf{M}_{0},\mathbf{M}_{0})=L(\mathbf{M}_{0}).$

The key is to choose a function for $g(\mathbf{M},\mathbf{M}_{0})$ that is
quadratic in $\mathbf{M}$ where its minimum can be computed analytically for
use in the next iteration as the new supporting point, thus converging toward
a minimum of the target function. If the target function is strictly convex
and coercive, this is the unique minimum (see, for example, Boyd and
Vandenberghe, 2004).

The first step in defining $g(\mathbf{M},\mathbf{M}_{0})$ is to majorize the
$\ell_{2}$-norm of the differences between rows of $\mathbf{M}$. If $a,b\geq
0$, then $(ab)^{1/2}\leq(a+b)/2$ since $0\leq(a^{1/2}-b^{1/2})^{2}$.
Substituting $a=\|\bm{\theta}\|^{2}$ and $b=\|\bm{\theta}_{0}\|^{2}$ yields

$\|\bm{\theta}\|\leq\frac{1}{2}\frac{\|\bm{\theta}\|^{2}}{\|\bm{\theta}_{0}\|}+\frac{1}{2}\|\bm{\theta}_{0}\|,$
(7)

for $\|\bm{\theta}_{0}\|>0$. Applying this inequality to
$\|\mathbf{m}_{k}-\mathbf{m}_{l}\|$, using
$\|\mathbf{m}_{k}-\mathbf{m}_{l}\|^{2}=(\mathbf{m}_{k}-\mathbf{m}_{l})^{\top}(\mathbf{m}_{k}-\mathbf{m}_{l})=\operatorname*{tr}(\mathbf{m}_{k}-\mathbf{m}_{l})(\mathbf{m}_{k}-\mathbf{m}_{l})^{\top}=\operatorname*{tr}\mathbf{M}^{\top}\mathbf{E}_{kl}\mathbf{M}$
subject to
$\mathbf{E}_{kl}=(\mathbf{e}_{k}-\mathbf{e}_{l})(\mathbf{e}_{k}-\mathbf{e}_{l})^{\top}$
with $\mathbf{e}_{k}$ being the $k$th column of the $c\times c$ identity
matrix, finally yields

$\displaystyle L(\mathbf{M})\leq$
$\displaystyle\;\kappa_{\epsilon}\|\mathbf{X-UM}\|^{2}+\frac{\lambda\kappa_{\text{pen}}}{2}\sum_{k<l}\mathbf{u}_{k}^{\top}\mathbf{Wu}_{l}\left(\frac{\|\mathbf{m}_{k}-\mathbf{m}_{l}\|^{2}}{\|\mathbf{m}_{0,k}-\mathbf{m}_{0,l}\|}+\|\mathbf{m}_{0,k}-\mathbf{m}_{0,l}\|\right)$
$\displaystyle=$
$\displaystyle\;\kappa_{\epsilon}\|\mathbf{X-UM}\|^{2}+\frac{\lambda\kappa_{\text{pen}}}{2}\operatorname*{tr}\mathbf{M}^{\top}\mathbf{C}_{0}\mathbf{M}+c(\mathbf{M}_{0}),$
(8)

where

$\mathbf{C}_{0}=\sum_{k<l}\frac{\mathbf{u}_{k}^{\top}\mathbf{Wu}_{l}}{\|\mathbf{m}_{0,k}-\mathbf{m}_{0,l}\|}\mathbf{E}_{kl},$

with $c(\mathbf{M}_{0})$ the collection of constant terms not depending on
$\mathbf{M}$. We should note that using cluster fusions in the algorithm
prevents a division by zero in the computation of $\mathbf{C}_{0}$.

At this point, the function in (8) is quadratic in $\mathbf{M}$ and is
continuously differentiable, aptly suiting it for use as
$g(\mathbf{M},\mathbf{M}_{0})$ in the MM algorithm. Here, the update for
$\mathbf{M}$ satisfies

$\displaystyle\frac{\partial
g(\mathbf{M},\mathbf{M}_{0})}{\partial\mathbf{M}}=2\kappa_{\epsilon}\mathbf{U}^{\top}\mathbf{UM}-2\kappa_{\epsilon}\mathbf{U}^{\top}\mathbf{X}+\lambda\kappa_{\text{pen}}\mathbf{C}_{0}\mathbf{M}$
$\displaystyle=\mathbf{O}$
$\displaystyle\left(\mathbf{U}^{\top}\mathbf{U}+\frac{\lambda\kappa_{\text{pen}}}{2\kappa_{\epsilon}}\mathbf{C}_{0}\right)\mathbf{M}$
$\displaystyle=\mathbf{U}^{\top}\mathbf{X},$ (9)

where $\mathbf{O}$ is a matrix of zeroes. However, solving the system of
equations in (9) incurs a complexity of $\mathcal{O}({c^{3}})$—a costly
operation when $c$ is large. As $\mathbf{U}^{\top}\mathbf{U}$ is a diagonal
matrix, the issue lies with the matrix $\mathbf{C}_{0}$. We thus take another
step in the majorization process by eliminating
$\operatorname*{tr}\mathbf{M}^{\top}\mathbf{C}_{0}\mathbf{M}$ from (8).
Groenen et al. (1999) treated a similar problem by carrying out a majorization
step using upper the bounds of the eigenvalues of $\mathbf{C}_{0}$. If
$\mathbf{D}_{0}$ is a diagonal matrix that features these upper bounds as its
elements, then $\mathbf{D}_{0}-\mathbf{C}_{0}$ is positive semi-definite.
Consequently, the following inequality can be derived

$\displaystyle 0$
$\displaystyle\leq\operatorname*{tr}[(\mathbf{M}-\mathbf{M}_{0})^{\top}(\mathbf{D}_{0}-\mathbf{C}_{0})(\mathbf{M}-\mathbf{M}_{0})]$
$\displaystyle\operatorname*{tr}\mathbf{M}^{\top}\mathbf{C}_{0}\mathbf{M}$
$\displaystyle\leq\operatorname*{tr}\mathbf{M}^{\top}\mathbf{D}_{0}\mathbf{M}-2\operatorname*{tr}\mathbf{M}^{\top}(\mathbf{D}_{0}-\mathbf{C}_{0})\mathbf{M}_{0}+\operatorname*{tr}\mathbf{M}_{0}^{\top}(\mathbf{D}_{0}-\mathbf{C}_{0})\mathbf{M}_{0}.$

To obtain $\mathbf{D}_{0}$, Gershgorin’s Circle Theorem (Bronson, 1989) can be
used where eigenvalues of $\mathbf{C}_{0}$ fall within $[c_{kk}-\sum_{l\neq
k}|c_{kl}|,c_{kk}+\sum_{l\neq k}|c_{kl}|]$. Due to the structure of
$\mathbf{C}_{0}$, the upper bound $c_{kk}+\sum_{l\neq k}|c_{kl}|=2c_{kk}$,
which means that $\mathbf{D}_{0}$ can be computed as twice the diagonal of
$\mathbf{C}_{0}$. Substituting the inequality for
$\operatorname*{tr}\mathbf{M}^{\top}\mathbf{C}_{0}\mathbf{M}$ into the
majorized loss function in (8) yields

$\displaystyle L(\mathbf{M})\leq$
$\displaystyle\;\kappa_{\epsilon}\|\mathbf{X-UM}\|^{2}+\frac{\lambda\kappa_{\text{pen}}}{2}\operatorname*{tr}\mathbf{M}^{\top}\mathbf{D}_{0}\mathbf{M}-\lambda\kappa_{\text{pen}}\operatorname*{tr}\mathbf{M}^{\top}(\mathbf{D}_{0}-\mathbf{C}_{0})\mathbf{M}_{0}+c(\mathbf{M}_{0})$
$\displaystyle=$ $\displaystyle\;g(\mathbf{M},\mathbf{M}_{0}).$

Then, the matrix $\mathbf{M}$ that minimizes $g(\mathbf{M},\mathbf{M}_{0})$
can be computed as

$\displaystyle\frac{\partial
g(\mathbf{M},\mathbf{M}_{0})}{\partial\mathbf{M}}=$
$\displaystyle\;2\kappa_{\epsilon}\mathbf{U}^{\top}\mathbf{UM}-2\kappa_{\epsilon}\mathbf{U}^{\top}\mathbf{X}+\lambda\kappa_{\text{pen}}\bigl{(}\mathbf{D}_{0}\mathbf{M}-(\mathbf{D}_{0}-\mathbf{C}_{0})\mathbf{M}_{0}\bigr{)}=\mathbf{O}$
$\displaystyle\mathbf{M}=$
$\displaystyle\;\left(\mathbf{U}^{\top}\mathbf{U}+\frac{\lambda\kappa_{\text{pen}}}{2\kappa_{\epsilon}}\mathbf{D}_{0}\right)^{-1}\left(\mathbf{U}^{\top}\mathbf{X}+\frac{\lambda\kappa_{\text{pen}}}{2\kappa_{\epsilon}}(\mathbf{D}_{0}-\mathbf{C}_{0})\mathbf{M}_{0}\right).$
(10)

Note that the first term is a diagonal matrix with elements

$\left(\mathbf{U}^{\top}\mathbf{U}+\frac{\lambda\kappa_{\text{pen}}}{2\kappa_{\epsilon}}\mathbf{D}_{0}\right)_{kk}^{-1}=\left(|\mathcal{C}_{k}|+\frac{\lambda\kappa_{\text{pen}}}{2\kappa_{\epsilon}}d_{0,kk}\right)^{-1}$

that will now permit straightforward updates. In Section 3.3, we present the
complete algorithm that uses this update to perform convex clustering.

###  3.3 Implementation

In Algorithm 1, we present the CCMM algorithm used to minimize (6). The
algorithm is divided into six steps, each performing a straightforward task.
Step 1 initializes several variables used in the minimization. Notable are
$\widetilde{\mathbf{X}}$ and $\widetilde{\mathbf{W}}$ used to track
$\mathbf{U}^{\top}\mathbf{X}$ and $\mathbf{U}^{\top}\mathbf{W}\mathbf{U}$ to
avoid unnecessary computations. The main loop of the algorithm starts by
computing the update per (10) in Step 2. To reduce the number of iterations
required to reach convergence, we make use of a relaxed update (De Leeuw and
Heiser, ), which we call “step doubling.” When step doubling, the update for
$\mathbf{M}$ is computed as $2\mathbf{M}_{1}-\mathbf{M}_{0}$, where
$\mathbf{M}_{1}$ is the result from Step 2 and $\mathbf{M}_{0}$ the supporting
point. In practice, this type of update is applied after a number of burn-in
iterations $b$ to avoid “stepping over” the minimum. This accelerated update
is computed in Step 3, and in our experiments we typically set $b$ to 25.

In Step 4, rows in $\mathbf{M}_{0}$ are fused. As Hocking et al. (2011) have
noted, the ideal approach for determining fusions would be to evaluate the
change in loss function per centroid pair. However, this may not be viable
when it incurs a large number of additional computations. Therefore, the
authors suggest using a threshold that triggers a cluster fusion when
$\|\mathbf{m}_{k}-\mathbf{m}_{l}\|\leq\varepsilon_{f}$. We propose a
formulation of $\varepsilon_{f}$ that reflects the scale of the data while
also securing a threshold not too small in cases of very similar (duplicate)
objects. Here, $\varepsilon_{f}$ computes as

$\varepsilon_{f}=\tau\operatorname*{median}_{i,j}\left(\|\mathbf{x}_{i}-\mathbf{x}_{j}\|\right),$
(11)

where $\tau$ denotes a specified fraction of the median of Euclidean distances
between objects. For very large data sets, the median of medians may estimate
this value. In our experiments, we typically set $\tau$ to $10^{-3}$.

Input: $\mathbf{X}$, $\mathbf{W}$, $\lambda$, $\varepsilon_{c}$,
$\varepsilon_{f}$, $b$

Output: $\mathbf{A}$

$\mathbf{M}_{0}:=\mathbf{X}$

$\triangleright$ Step 1

$\mathbf{U}:=\mathbf{I}_{n}$

$L_{1}:=L(\mathbf{M}_{0})$

$L_{0}:=(1+2\varepsilon_{c})L_{1}$

$\widetilde{\mathbf{X}}:=\mathbf{X}$

$\widetilde{\mathbf{W}}:=\mathbf{W}$

$\kappa_{\epsilon}:=\bigl{(}{2\|\mathbf{X}\|^{2}}\bigr{)}^{-1}$

$\kappa_{\text{pen}}:=\bigl{(}{\|\mathbf{X}\|\sum_{i<j}w_{ij}}\bigr{)}^{-1}$

$\gamma:=\lambda\kappa_{\text{pen}}/(2\kappa_{\epsilon})\;$

$t:=1$

while  _$(L_{0}-L_{1})/L_{1} >\varepsilon_{c}$_ do

$\mathbf{C}_{0}:=\sum_{k=2}^{c}\sum_{l=1}^{k-1}(\widetilde{w}_{kl}/\|\mathbf{m}_{0,k}-\mathbf{m}_{0,l}\|)\mathbf{E}_{kl}$

$\triangleright$ Step 2

$\mathbf{D}_{0}:=\operatorname*{diag}(\mathbf{C}_{0})$

$\mathbf{M}_{1}:=\bigl{(}\mathbf{U}^{\top}\mathbf{U}+\gamma\mathbf{D}_{0}\bigr{)}^{-1}\bigl{(}\widetilde{\mathbf{X}}+\gamma(\mathbf{D}_{0}-\mathbf{C}_{0})\mathbf{M}_{0}\bigr{)}$

if  _$t >b$_ then $\triangleright$ Step 3

$\mathbf{M}_{0}:=2\mathbf{M}_{1}-\mathbf{M}_{0}$

else

$\mathbf{M}_{0}:=\mathbf{M}_{1}$

while
_$\exists(k,l):\|\mathbf{m}_{0,k}-\mathbf{m}_{0,l}\|\leq\varepsilon_{f}$_ do
$\triangleright$ Step 4

Fuse rows of $\mathbf{M}_{0}$ and obtain $\mathbf{U}_{new}$

$\widetilde{\mathbf{X}}:=\mathbf{U}_{new}^{\top}\widetilde{\mathbf{X}}$

$\widetilde{\mathbf{W}}:=\mathbf{U}_{new}^{\top}\widetilde{\mathbf{W}}\mathbf{U}_{new}$

$\mathbf{U}:=\mathbf{U}\mathbf{U}_{new}$

$L_{0}:=L_{1}$

$\triangleright$ Step 5

$L_{1}:=L(\mathbf{M}_{0})$

$t:=t+1$

$\mathbf{A}:=\mathbf{U}\mathbf{M}_{0}$

$\triangleright$ Step 6

Algorithm 1 Convex clustering through majorization-minimization (CCMM). Key
parts of the algorithm are Step 2 where the update derived in (10) is
computed, and Step 4 where cluster fusions are performed. Step 3 is executed
only after a number of burn-in iterations to accelerate convergence. Remaining
steps perform initializations and bookkeeping.

The fusion process outputs the $c_{old}\times c_{new}$ membership matrix
$\mathbf{U}_{new}$ that maps the old clusters to the new ones. This matrix can
be used to update $\widetilde{\mathbf{X}}$ and $\widetilde{\mathbf{W}}$ as
$\mathbf{U}_{new}^{\top}\widetilde{\mathbf{X}}$ and
$\mathbf{U}_{new}^{\top}\widetilde{\mathbf{W}}\mathbf{U}_{new}$, respectively.
Additionally, Step 4 repeats if there are rows in the updated matrix
$\mathbf{M}_{0}$ for which
$\|\mathbf{m}_{k}-\mathbf{m}_{l}\|\leq\varepsilon_{f}$.

In Step 5, the loss function value is updated. When the relative decrease in
the loss exceeds some value $\varepsilon_{c}$, the algorithm reverts to Step
2. At convergence, the resulting matrix $\mathbf{A}$ is computed and returned
to the user. For the results presented in this paper, we set the value of
$\varepsilon_{c}$ to $10^{-6}$. Finally, when computing a clusterpath for
multiple values of $\lambda$, warm starts can be used by setting
$\mathbf{M}_{0}$ to the solution of the minimization for the previous value of
$\lambda$ in Step 1.

To analyze algorithmic complexity, we focus on two of the six steps. First, in
Step 2, the update derived in (10) is computed. A key feature of this update
is that it can be done in a single loop over the nonzero elements of the
$c\times c$ matrix $\widetilde{\mathbf{W}}$. Hence, this step has a complexity
of $\mathcal{O}({npk})$ where $n$ is an upper bound on $c$. Second, in Step 4,
the rows of $\mathbf{M}_{0}$ fuse, and several matrices update. Using the
nonzero elements of $\widetilde{\mathbf{W}}$, a single pass through Step 4
also has complexity $\mathcal{O}({npk})$. However, since Step 4 may recur, it
can be shown in the absolute worst case that the number of passes is $n$, and
worst-case complexity is $\mathcal{O}({n^{2}pk})$.

The number of repetitions in Steps (2-5) depends on how far from the minimum
the algorithm is initialized. Near the optimum, the asymptotic convergence
rate of MM algorithms is linear (De Leeuw, ). However, it has been shown that
during the first few iterations the asymptotic linear rate does not apply, and
convergence is faster (Havel, 1991).

The CCMM algorithm was coded in C++ while making extensive use of the linear
algebra library Eigen (Guennebaud and Jacob, 2010). We used Rcpp (Eddelbuettel
and François, 2011; Eddelbuettel, 2013; Eddelbuettel and Balamuta, 2017) to
interface the program with the programming language R and pybind11 (Jakob et
al., 2017) to interface it with Python.111Software packages that implement
CCMM for both R and Python are available at https://github.com/djwtouw/CCMMR
and https://github.com/djwtouw/CCMMPy. A key factor in scalability is sparsity
in the weights as defined in (4). To efficiently find the $k$ nearest
neighbors of each object, we enlisted a $k$-d tree (Bentley, 1975) implemented
by Arya et al. (2019) and Pedregosa et al. (2011).

##  4 Numerical Experiments

In this section, we present results of several numerical experiments. The main
contribution of our research is an algorithm that performs convex clustering
more efficiently than prior published algorithms. Therefore, we first present
a comparison between CCMM versus its main competitors. All computational
results in this paper were obtained using a desktop featuring an AMD Ryzen 5
2600X processor at 3.6 GHz using 16 GB of memory running Ubuntu version
20.04.1.

###  4.1 Simulated Data

To be competitive, CCMM should minimize (2) faster than state-of-the-art
convex clustering algorithms including, to our knowledge, today’s fastest
approach: SSNAL (Sun et al., 2021). Beyond comparing CCMM to SSNAL, we
consider AMA (Chi and Lange, 2015), which is available as an R-package.222The
source code for AMA can be found at
https://cran.r-project.org/src/contrib/Archive/cvxclustr/ and the source code
for SSNAL at https://blog.nus.edu.sg/mattohkc/softwares/convexclustering/.
Additionally, the code to generate the results presented in this section is
available at https://github.com/djwtouw/CCMM-paper.

For our contest of three algorithms, we chose an approach similar to that
described in Sun et al. (2021). Data sets ranging from 1,000 to 5,000 objects
were generated per the two interlocking half-moons data generating process
(see Figure 2). For each value of $n$, we used ten different realizations of
the data. Since AMA and SSNAL typically minimize the unscaled loss in (1), we
adjusted the scaling constants in CCMM so that each algorithm minimized the
same problem. Additionally, as the software package of SSNAL does not allow
for custom (connected) weight matrices, we used the weights from (3) without
the scaling component, opting to set $\phi=2$ and $k=15$. The final input
common to the three algorithms was the sequence for $\lambda$ that describes
the clusterpath, set to $\\{0.0,0.2,\ldots,109.8,110.0\\}$.

For the CCMM algorithm, we set $\varepsilon_{c}$ to $10^{-6}$ and $b$ to 25,
computing $\varepsilon_{f}$ according to (11) with $\tau=10^{-3}$. For AMA,
the default tolerance for convergence of $10^{-3}$ was used along with the
default step size parameter $\nu$ (which depends on $\mathbf{W}$ and $n$). For
parameters in the SSNAL algorithm, we applied default values provided by the
authors. The tolerance for terminating the algorithm was set to $10^{-6}$. For
warm starts, 50 ADMM iterations were used for the first value for $\lambda$
and 20 iterations for subsequent values.

(a) Average elapsed time per clusterpath

(b) Average elapsed time per CCMM iteration

Figure 4: In Panel (a): the mean computation times for performing convex
clustering on the two interlocking half moon data for the number of objects
ranging from 1,000 to 5,000. In Panel (b): the mean computation time per
iteration of the CCMM algorithm for the two interlocking half moons ranging
from 6,000 to 20,000 objects.

In Figure 4a, we present mean elapsed times obtained experimentally. Results
show that for $n=$ 1,000, CCMM was roughly 8,610 times faster than AMA and 236
times faster than SSNAL. For $n=$ 5,000, speed factors rose to 10,486 and 478
times, respectively. This convincingly proves that CCMM is the fastest
algorithm to perform convex clustering. For Figure 4b, we continued the
analysis of the two half-moon data sets for $n$ ranging from 6,000 to 20,000
using only CCMM. The figure shows the average elapsed time per minimization
iteration. Even though we derived in Section 3.3 each iteration’s theoretical
worst-case complexity to be $\mathcal{O}({n^{2}pk})$, we observe a trend
linear in $n$.

To be a viable alternative to SSNAL, a new algorithm should not only minimize
the loss function in less time, but also yield a similar solution. To assess
the solution quality of AMA, SSNAL and CCMM, we compared values of the loss
function that each algorithm obtained for data sets comprising 1,000 and 5,000
objects. In Figure 5, we use the SSNAL output as a benchmark for reporting the
relative values obtained by AMA and CCMM. Being the slowest, AMA also incurs
the largest values for minimizations of (1). As Figure 5b shows, the relative
quality of the minimum obtained by AMA is marred by the sample size, while
this is not the case for CCMM. The minimum attained by CCMM deviates at most
0.01% from that attained by SSNAL for both $n=$ 1,000 and $n=$ 5,000.

(a) Relative losses for $n=$ 1,000

(b) Relative losses for $n=$ 5,000

Figure 5: In Panel (a): the mean value of the loss function obtained by CCMM,
SSNAL, and AMA relative to the value obtained by SSNAL. Computed for ten
realizations of the two interlocking half moon data with $n=$ 1,000. In Panel
(b): the same as in Panel (a), but for $n=$ 5,000.

###  4.2 Empirical Data

Besides simulated data, we used several empirical data sets to compare the
three algorithms: the banknotes authentication, musk and MAGIC gamma telescope
data sets (Dua and Graff, 2019), and the MNIST data set (Deng, 2012). In Table
1, the first two columns feature the number of objects ($n$) and variables
($d$) for each data set. To attain reasonable computation times, we used UMAP
(McInnes et al., 2018) to reduce dimensionality to the number of categories in
the data ($p$). Next, the weight matrices were computed using $k=50$ and
$\phi=2$. To determine the sequence for $\lambda$, we first obtained the
smallest value where the number of clusters stopped decreasing. Then, the
sequence was chosen as 200 values along a cubic function from $0$ to
$\lambda_{max}$. Notably, omitting a symmetric circulant matrix to ensure a
connected weight matrix meant that none of the clusterpaths terminated in a
single cluster. The fourth column lists the minimum number of clusters
$c_{min}$.

Table 1: Clusterpaths for real data sets computed by three different algorithms. Dimension reduction was applied to each $n\times d$ data set using UMAP where the target number of dimensions ($p$) was the true number of clusters. Each clusterpath was computed for 200 values of $\lambda$ that terminated at the minimum number of clusters attainable ($c_{min}$). This value exceeds one since the software package for SSNAL forbids the use of custom weight matrices. The three rightmost columns list the number of seconds required by each algorithm to compute the final clusterpath. A hyphen means that the algorithm was unable to complete the minimization. |  |  |  |  |  |  | Elapsed time (s)
---|---|---|---|---|---|---|---
|  | $n$ | $d$ | $p$ | $c_{min}$ |  | AMA | SSNAL | CCMM
Banknotes |  | 1,372 | 4 | 2 | 7 |  | 402.23 | 51.75 | 0.18
Musk |  | 6,598 | 166 | 2 | 27 |  | - | 197.34 | 0.11
Gamma |  | 19,020 | 10 | 2 | 2 |  | - | 6,343.01 | 20.09
MNIST |  | 60,000 | 784 | 10 | 6 |  | - | - | 284.43

The number of seconds passed for each of the algorithms while computing the
clusterpaths is reported in the three rightmost columns of Table 1. A missing
result signals insufficient memory (16 GB system) for that algorithm. Results
show that the speedup of CCMM over the other algorithms is dependent on the
data analyzed. CCMM is about 2,200 times faster than AMA for the banknotes
data, which differs significantly from the results for the generated data
sets. The same holds for the speedup with respect to SSNAl, which ranges from
roughly 285 to almost 1,800.

###  4.3 Analysis of a Large Data Set

In prior sections, we have established CCMM as the most competitive algorithm
to perform large-scale convex clustering. Here, we showcase CCMM analyzing a
data set containing more than one million measurements of the power
consumption by a single household located in Sceaux—a commune part of the
Paris metropolitan area333The individual household electric power consumption
data set is made available in the UCI Machine Learning Repository (Dua and
Graff, 2019)..

Data consist of measurements sampled at one-minute intervals from the 16th of
December 2006 through 16th of December 2008. Metrics include the global active
and reactive power (in kilowatts), voltage (in volts), global intensity (in
amperes), and the activity at three sub-metering stations (in watt-hours). The
first sub-metering station monitored a kitchen, the second a laundry room, and
the third a water-heater and air-conditioner. In total, there are 1,048,571
measurements, and the mean of each variable is reported in the first column of
Table 2. Before performing the cluster analysis, each variable was
standardized to have a mean of zero and variance of one. The weight matrix was
computed using $k=15$ and $\phi=0.5$ in combination with the symmetric
circulant matrix to ensure a complete cluster hierarchy.

Table 2: Variable means for the individual household electric power consumption data for the complete data set (in the leftmost column) and for the two clusters (A and B) that were found using convex clustering. |  | Complete | Cluster A | Cluster B
---|---|---|---|---
Global active power |  | 1.11 | 0.62 | 2.15
Global reactive power |  | 0.12 | 0.11 | 0.13
Voltage |  | 239.97 | 240.65 | 238.52
Global intensity |  | 4.72 | 2.68 | 9.05
Energy sub-metering No. 1 |  | 1.18 | 0.63 | 2.35
Energy sub-metering No. 2 |  | 1.47 | 0.97 | 2.54
Energy sub-metering No. 3 |  | 5.94 | 0.39 | 17.71
Number of measurements |  | 1,048,571 | 712,749 | 335,821

For this analysis, we took an approach distinct from computing the clusterpath
for a predetermined sequence for $\lambda$. While useful for comparing
different algorithms, it demands advance knowledge of which value for
$\lambda$ maps to a specific cluster number output. For this analysis, we thus
initialized $\lambda$ at 0.01 and repeatedly minimized the loss function,
increasing $\lambda$ by 2.5% after solving each minimization. As with
computing the clusterpath, the solution for the minimization problem at the
prior value for $\lambda$ was used as a warm start for the next value for
$\lambda$. After reaching 20 clusters, the procedure entered a refinement
phase. If the number of clusters was reduced by more than one as $\lambda$
increased, a midpoint value was chosen to obtain a cluster hierarchy as
detailed as possible. The procedure was terminated after reaching two clusters
since the final reduction in cluster number is trivial. In total, 984
instances of the convex-clustering loss function were solved. The entire
procedure was completed in 19.1 hours, thus solving each instance in 70
seconds, on average. To our knowledge, this is the largest data set analyzed
to date via convex clustering.

In Figure 6a, a dendrogram in combination with a cluster heatmap shows the
agglomeration of the final 20 clusters. Since these clusters are of different
sizes, the height of each cell row in the heatmap is proportional to the
logarithm of the size of the corresponding cluster. Fifteen of the final
twenty clusters contain two or fewer measurements, and the largest two
combined comprise over 95% of the data. The heatmap distinguishes these final
two clusters by the measurements at the third energy sub-metering station.
Looking at the descriptive statistics of the final two clusters in the second
and third columns of Table 2, we can conclude that cluster A corresponds to
measurements of lower power consumption, B to greater consumption.

(a) Cluster heatmap

(b) Distribution of the final two clusters

Figure 6: In Panel (a): a cluster heatmap for the final 20 clusters in the
individual household electric power consumption data. The height of each row
is proportional to the logarithm of the number of measurements in each of 20
clusters. Cell color is determined by the mean value of the cluster for that
particular variable. In Panel (b): the distribution of timestamps in the final
two clusters. Frequency is determined by the number of times each timestamp
occurred in one of the clusters.

To gain further insight into the final two clusters, Figure 6b shows the
timestamp frequency for both clusters. Energy consumption is generally low
during nighttime when members of the household are likely sleeping. Around
08:00 and 20:00, peaks emerge for cluster B to indicate general times of day
that people are home, activating the heater or air-conditioner, or performing
tasks like cooking or laundering. It should be noted that household members
may stick to different time schedules in their power usage activities during
weekends, holidays, and changing seasons. Still, this figure conveys a dense
level of information about how this household generally spends its days.

##  5 Conclusion

In this paper, we propose a new, efficient algorithm to minimize the convex
clustering loss function and compute the clusterpath. We suggest using a
connected, sparse weight matrix in combination with cluster fusions to
guarantee a complete cluster hierarchy. At the heart of our minimization
algorithm, which we call CCMM, lies an efficient update derived using
majorization-minimization.

Our experiments show that CCMM is up to 10,486 times faster than AMA and up to
478 times faster than SSNAL while losing, at most, 0.01% accuracy.
Furthermore, we analyzed a data set featuring over one million measurements of
the electric power consumption by a single household. On average, solving each
instance of (2) for a particular value for $\lambda$ took 70 seconds. To our
knowledge, this data set is the largest ever clustered using convex
clustering.

One useful future research path may merge CCMM and convex clustering via
algorithmic regularization paths (CARP). Originally developed by Weylandt et
al. (2020), CARP uses small increments in the regularization parameter to
compute a highly detailed clusterpath. For each value of $\lambda$, its prior
solution is used as a warm start, and an ADMM algorithm (Boyd et al., 2011) is
allowed to perform one iteration. In their paper, the authors showed that if
the increments in $\lambda$ are small enough, the difference between the
results of CARP and the exact solution approaches zero. Given the efficiency
of CCMM, we expect that this could suitably replace the ADMM update in CARP.

In sum, CCMM is an algorithm scaling better than state-of-the-art
implementations designed to perform convex clustering. Its excellent
scalability facilitates the analysis of larger data sets in less time than
possible to date.

## A

In this appendix, we prove that the convex clustering loss function (2) is
strongly convex.

Proof  Let $\mathbf{A}$ and $\mathbf{B}$ be points in the domain of the convex
clustering loss function $L$ and $t\in[0,1]$. To prove the strong convexity
(see, e.g., Nesterov, 2004), we need to show that there exists $m>0$ such that

$tL(\mathbf{A})+(1-t)L(\mathbf{B})-L(t\mathbf{A}+(1-t)\mathbf{B})\geq\frac{1}{2}mt(1-t)\|\mathbf{A-B}\|^{2}.$
(12)

The left-hand side (LHS) of (12) can be rewritten as

$\displaystyle\text{LHS}=$
$\displaystyle\;\kappa_{\epsilon}\bigl{(}t\|\mathbf{X-A}\|^{2}+(1-t)\|\mathbf{X-B}\|^{2}-\|\mathbf{X}-t\mathbf{A}-(1-t)\mathbf{B}\|^{2}\bigr{)}$
$\displaystyle\;+\lambda\kappa_{\text{pen}}\sum_{i<j}w_{ij}\bigl{(}t\|\mathbf{a}_{i}-\mathbf{a}_{j}\|+(1-t)\|\mathbf{b}_{i}-\mathbf{b}_{j}\|-\|t(\mathbf{a}_{i}-\mathbf{a}_{j})+(1-t)(\mathbf{b}_{i}-\mathbf{b}_{j})\|\bigr{)}$
$\displaystyle=$
$\displaystyle\;\kappa_{\epsilon}t(1-t)\|\mathbf{A}-\mathbf{B}\|^{2}$
$\displaystyle\;+\lambda\kappa_{\text{pen}}\sum_{i<j}w_{ij}\bigl{(}t\|\mathbf{a}_{i}-\mathbf{a}_{j}\|+(1-t)\|\mathbf{b}_{i}-\mathbf{b}_{j}\|-\|t(\mathbf{a}_{i}-\mathbf{a}_{j})+(1-t)(\mathbf{b}_{i}-\mathbf{b}_{j})\|\bigr{)}.$

By the triangle inequality, the second term is greater or equal to zero.
Therefore, (12) holds with $m>0$ if the normalization constant
$\kappa_{\epsilon}>0$.

## B

This appendix contains the proof of the following theorem from Section 2.3:

Theorem 1 If $\mathbf{W}$ is connected, then for some sufficiently large,
finite value of $\lambda$, a global minimum of (2) exists where all rows of
the solution $\mathbf{A}$ are identical to the mean of the rows in
$\mathbf{X}$.

Proof  The point of interest is the location where the penalty term is equal
to zero. This is the case if $\mathbf{a}_{i}=\mathbf{a}_{j}$ for each
$w_{ij}>0$. In this stage of the proof, the actual value of the rows is not of
interest yet, only the fact that they are identical. Note that the penalty on
the differences between the rows of $\mathbf{A}$ has a nondifferentiable point
if at least one of the distances is zero. To shorten notation, consider the
following definition for the penalty term

$P(\mathbf{A})=\kappa_{\text{pen}}\sum_{i<j}w_{ij}\|\mathbf{a}_{i}-\mathbf{a}_{j}\|.$

For the first stage, let $\mathbf{\hat{A}}$ be any value for $\mathbf{A}$ such
that $P(\mathbf{\hat{A}})=0$. Let the $n\times p$ matrix $\mathbf{\Sigma}$ be
a deviation from $\mathbf{\hat{A}}$ under two restrictions. First,
$\mathbf{\Sigma}$ must satisfy $P(\mathbf{\Sigma})\neq 0$ and second, we only
consider values for $\mathbf{\Sigma}$ such that

$\displaystyle\kappa_{\epsilon}\|\mathbf{X-A}\|^{2}$
$\displaystyle>\kappa_{\epsilon}\|\mathbf{X-(A+\Sigma)}\|^{2}$
$\displaystyle\operatorname*{tr}(\mathbf{X}-\mathbf{\hat{A}})^{\top}(\mathbf{X-\hat{A}})$
$\displaystyle>\operatorname*{tr}(\mathbf{X}-\mathbf{\hat{A}}-\mathbf{\Sigma})^{\top}(\mathbf{X-\hat{A}}-\mathbf{\Sigma})$
$\displaystyle\operatorname*{tr}(\mathbf{X}-\mathbf{\hat{A}})^{\top}(\mathbf{X-\hat{A}})$
$\displaystyle>\operatorname*{tr}(\mathbf{X}-\mathbf{\hat{A}})^{\top}(\mathbf{X-\hat{A}})-2\operatorname*{tr}(\mathbf{X-\hat{A}})^{\top}\mathbf{\Sigma}+\operatorname*{tr}\mathbf{\Sigma}^{\top}\mathbf{\Sigma}$
$\displaystyle\operatorname*{tr}(\mathbf{X-\hat{A}})^{\top}\mathbf{\Sigma}$
$\displaystyle>\frac{1}{2}\operatorname*{tr}\mathbf{\Sigma}^{\top}\mathbf{\Sigma}.$

The motivation for this second restriction is that the deviation must decrease
the first term of the loss function in (2). As the first restriction causes
the penalty term to increase, a choice for $\mathbf{\Sigma}$ that also
increases the first term surely increases the entire loss. Given that we are
investigating the location of the minimum of the loss function, choices for
$\mathbf{\Sigma}$ that do not satisfy the second restriction are not of
interest. The following shows that for any combination of $\mathbf{\hat{A}}$
and $\mathbf{\Sigma}$ there exists a finite value for $\lambda$ such that the
loss corresponding to $\mathbf{\hat{A}}$ is lower than the loss corresponding
to $\mathbf{\hat{A}+\Sigma}$

$\displaystyle L(\mathbf{\hat{A}})<$
$\displaystyle\;L(\mathbf{\hat{A}+\Sigma})$
$\displaystyle\kappa_{\epsilon}\operatorname*{tr}(\mathbf{X-\hat{A}})^{\top}(\mathbf{X-\hat{A}})<$
$\displaystyle\;\kappa_{\epsilon}\operatorname*{tr}(\mathbf{X-\hat{A}-\Sigma})^{\top}(\mathbf{X-\hat{A}-\Sigma})+\lambda
P(\mathbf{\Sigma})$
$\displaystyle\kappa_{\epsilon}\bigl{(}2\operatorname*{tr}\mathbf{\Sigma}^{\top}(\mathbf{X-\hat{A}})-\operatorname*{tr}\mathbf{\Sigma}^{\top}\mathbf{\Sigma}\bigr{)}<$
$\displaystyle\;\lambda P(\mathbf{\Sigma})$ $\displaystyle\lambda>$
$\displaystyle\;\kappa_{\epsilon}\frac{2\operatorname*{tr}\mathbf{\Sigma}^{\top}(\mathbf{X-\hat{A}})-\operatorname*{tr}\mathbf{\Sigma}^{\top}\mathbf{\Sigma}}{P(\mathbf{\Sigma})}.$

The two restrictions imposed on $\mathbf{\Sigma}$ ensure that the resulting
bound for $\lambda$ is positive and that there is no division by zero. This
result shows that for any $\mathbf{\hat{A}}$, there exists a value for
$\lambda$ such that a deviation for which the penalty term of the loss does
not equal zero results in a higher loss.

The next step is to show that the $\mathbf{\hat{A}}$ that minimizes the loss
in (2) is equal to $\bm{1}\mathbf{\bar{x}}^{\top}$. In order to show this, we
first impose the restriction that $\mathbf{W}$ is connected such that
$P(\mathbf{\hat{A}})=0$ if and only if all rows of $\mathbf{\hat{A}}$ are
identical. As the previous part of the proof shows that there exists a finite
value for $\lambda$ such that any deviation $\mathbf{\Sigma}$ from any
$\mathbf{\hat{A}}$ yields a higher value for the loss, the minimum of the loss
function must satisfy $P(\mathbf{A})=0$. Therefore, to find the minimum we now
minimize

$L(\mathbf{\hat{A}})=\kappa_{\epsilon}\sum_{i=1}^{n}\|\mathbf{x}_{i}-\mathbf{\hat{a}}\|^{2},$

for which the solution is
$\mathbf{\hat{a}}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{x}_{i}=\mathbf{\bar{x}}$.
This completes the proof that for a connected weight matrix $\mathbf{W}$, the
minimum of (2) for a sufficiently large value for $\lambda$ is attained at
$\mathbf{A}=\bm{1}\mathbf{\bar{x}}^{\top}$.

## C

This appendix contains the proof of the following theorem from Section 3.1:

Theorem 3 Let the matrix $\mathbf{M}$ contain the $c$ unique rows of
$\mathbf{A}$, and let $\varepsilon_{f}$ be the threshold used for fusing rows
of $\mathbf{M}$. Assume that
$\|\mathbf{m}_{k}-\mathbf{m}_{l}\|\leq\varepsilon_{f}$ for some $k\neq l$. Let
$\mathbf{m}_{new}$ be the weighted average of $\mathbf{m}_{k}$ and
$\mathbf{m}_{l}$ computed as in (5), and let $\mathbf{M}_{new}$ be the matrix
that results from setting $\mathbf{m}_{k}$ and $\mathbf{m}_{l}$ to
$\mathbf{m}_{new}$. Then, for fixed $\lambda$ and $\mathbf{W}$, the absolute
error $|L(\mathbf{M})-L(\mathbf{M}_{new})|$ tends toward zero as
$\varepsilon_{f}\rightarrow 0$.

Proof  If $\|\mathbf{m}_{k}-\mathbf{m}_{l}\|\leq\varepsilon_{f}$, the new
cluster centroid is computed as in (5). Let

$\displaystyle\bm{\sigma}_{k}$
$\displaystyle=\mathbf{m}_{k}-\mathbf{m}_{new}=\frac{|C_{l}|}{|C_{k}|+|C_{l}|}(\mathbf{m}_{k}-\mathbf{m}_{l})$
$\displaystyle\bm{\sigma}_{l}$
$\displaystyle=\mathbf{m}_{l}-\mathbf{m}_{new}=\frac{|C_{k}|}{|C_{k}|+|C_{l}|}(\mathbf{m}_{l}-\mathbf{m}_{k}),$

such that

$\|\bm{\sigma}_{k}\|\leq\frac{|C_{l}|}{|C_{k}|+|C_{l}|}\varepsilon_{f}\quad\text{and}\quad\|\bm{\sigma}_{l}\|\leq\frac{|C_{k}|}{|C_{k}|+|C_{l}|}\varepsilon_{f}.$
(13)

Writing out the absolute value of the difference between the loss before and
after merging $\mathbf{m}_{k}$ and $\mathbf{m}_{l}$, we obtain

$\displaystyle|L(\mathbf{M})-L(\mathbf{M}_{new})|=$
$\displaystyle\;\Biggl{|}\kappa_{\epsilon}\Biggl{(}\sum_{i\in
C_{k}}\|\mathbf{x}_{i}-\mathbf{m}_{k}\|^{2}+\sum_{i\in
C_{l}}\|\mathbf{x}_{i}-\mathbf{m}_{l}\|^{2}$ $\displaystyle\;-\sum_{i\in
C_{k}}\|\mathbf{x}_{i}-\mathbf{m}_{new}\|^{2}-\sum_{i\in
C_{l}}\|\mathbf{x}_{i}-\mathbf{m}_{new}\|^{2}\Biggr{)}$
$\displaystyle\;+\lambda\kappa_{\text{pen}}\Biggl{(}\sum_{\begin{subarray}{c}i=0\\\
i\notin\\{k,l\\}\end{subarray}}^{c}\mathbf{u}_{k}^{\top}\mathbf{Wu}_{i}\|\mathbf{m}_{k}-\mathbf{m}_{i}\|$
$\displaystyle\;+\sum_{\begin{subarray}{c}i=0\\\
i\notin\\{k,l\\}\end{subarray}}^{c}\mathbf{u}_{l}^{\top}\mathbf{Wu}_{i}\|\mathbf{m}_{l}-\mathbf{m}_{i}\|+\mathbf{u}_{k}^{\top}\mathbf{Wu}_{l}\|\mathbf{m}_{k}-\mathbf{m}_{l}\|$
$\displaystyle\;-\sum_{\begin{subarray}{c}i=0\\\
i\notin\\{k,l\\}\end{subarray}}^{c}\mathbf{u}_{k}^{\top}\mathbf{Wu}_{i}\|\mathbf{m}_{new}-\mathbf{m}_{i}\|-\sum_{\begin{subarray}{c}i=0\\\
i\notin\\{k,l\\}\end{subarray}}^{c}\mathbf{u}_{l}^{\top}\mathbf{Wu}_{i}\|\mathbf{m}_{new}-\mathbf{m}_{i}\|\Biggr{)}\Biggr{|},$

where the first two lines concern the difference between the first terms of
$L(\mathbf{M})$ and $L(\mathbf{M}_{new})$ and the remaining three lines
concern the difference between the values of the penalty terms. Using the
triangle inequality, this can be written as

$|L(\mathbf{M})-L(\mathbf{M}_{new})|\leq\kappa_{\epsilon}|\Gamma_{1}|+\lambda\kappa_{\text{pen}}|\Gamma_{2}|,$

where $\Gamma_{1}$ is the result of the terms on the first two lines and
$\Gamma_{2}$ is the result of the terms on the remaining four lines. We start
by providing an upper bound for $|\Gamma_{1}|$. Rewriting yields

$\displaystyle|\Gamma_{1}|$ $\displaystyle=\Biggl{|}\sum_{i\in
C_{k}}\bigl{(}\|\mathbf{x}_{i}-\mathbf{m}_{k}\|^{2}-\|\mathbf{x}_{i}-\mathbf{m}_{new}\|^{2}\bigr{)}+\sum_{i\in
C_{l}}\bigl{(}\|\mathbf{x}_{i}-\mathbf{m}_{l}\|^{2}-\|\mathbf{x}_{i}-\mathbf{m}_{new}\|^{2}\bigr{)}\Biggr{|}$
$\displaystyle=\Biggl{|}\sum_{i\in
C_{k}}\bigl{(}\|\mathbf{x}_{i}-\mathbf{m}_{k}\|^{2}-\|\mathbf{x}_{i}-\mathbf{m}_{k}+\bm{\sigma}_{k}\|^{2}\bigr{)}+\sum_{i\in
C_{l}}\bigl{(}\|\mathbf{x}_{i}-\mathbf{m}_{l}\|^{2}-\|\mathbf{x}_{i}-\mathbf{m}_{l}+\bm{\sigma}_{l}\|^{2}\bigr{)}\Biggr{|}$
$\displaystyle=\Biggl{|}\sum_{i\in
C_{k}}\bigl{(}\bm{\sigma}_{k}^{\top}\bm{\sigma}_{k}-2(\mathbf{x}_{i}-\mathbf{m}_{k})^{\top}\bm{\sigma}_{k}\bigr{)}+\sum_{i\in
C_{l}}\bigl{(}\bm{\sigma}_{l}^{\top}\bm{\sigma}_{l}-2(\mathbf{x}_{i}-\mathbf{m}_{l})^{\top}\bm{\sigma}_{l}\bigr{)}\Biggr{|},$

which is bounded by

$\displaystyle|\Gamma_{1}|$ $\displaystyle\leq\sum_{i\in
C_{k}}\bigl{(}\bm{\sigma}_{k}^{\top}\bm{\sigma}_{k}+2|(\mathbf{x}_{i}-\mathbf{m}_{k})^{\top}\bm{\sigma}_{k}|\bigr{)}+\sum_{i\in
C_{l}}\bigl{(}\bm{\sigma}_{l}^{\top}\bm{\sigma}_{l}+2|(\mathbf{x}_{i}-\mathbf{m}_{l})^{\top}\bm{\sigma}_{l}|\bigr{)}$
$\displaystyle\leq\sum_{i\in
C_{k}}\bigl{(}\bm{\sigma}_{k}^{\top}\bm{\sigma}_{k}+2\|\mathbf{x}_{i}-\mathbf{m}_{k}\|\|\bm{\sigma}_{k}\|\bigr{)}+\sum_{i\in
C_{l}}\bigl{(}\bm{\sigma}_{l}^{\top}\bm{\sigma}_{l}+2\|\mathbf{x}_{i}-\mathbf{m}_{l}\|\|\bm{\sigma}_{l}\|\bigl{)}$
$\displaystyle\leq\frac{|C_{k}|\,|C_{l}|}{|C_{k}|+|C_{l}|}\varepsilon_{f}^{2}+\frac{2}{|C_{k}|+|C_{l}|}\Biggl{(}|C_{l}|\sum_{i\in
C_{k}}\|\mathbf{x}_{i}-\mathbf{m}_{k}\|+|C_{k}|\sum_{i\in
C_{l}}\|\mathbf{x}_{i}-\mathbf{m}_{l}\|\Biggr{)}\varepsilon_{f},$

where we made use of the triangle inequality, the Cauchy-Schwarz inequality,
and (13). Each of the terms in the upper bound is now proportional to
$\varepsilon_{f}$. To bound $|\Gamma_{2}|$, consider

$\displaystyle\|\mathbf{m}_{k}-\mathbf{m}_{i}\|-\|\mathbf{m}_{new}-\mathbf{m}_{i}\|$
$\displaystyle=\|\mathbf{m}_{k}-\mathbf{m}_{i}\|-\|\mathbf{m}_{k}-\bm{\sigma}_{k}-\mathbf{m}_{i}\|$
$\displaystyle=\|\bm{\sigma}_{k}+\mathbf{m}_{k}-\bm{\sigma}_{k}-\mathbf{m}_{i}\|-\|\mathbf{m}_{k}-\bm{\sigma}_{k}-\mathbf{m}_{i}\|,$

and by means of the triangle inequality we obtain

$\displaystyle\|\bm{\sigma}_{k}+\mathbf{m}_{k}-\bm{\sigma}_{k}-\mathbf{m}_{i}\|$
$\displaystyle\leq\|\bm{\sigma}_{k}\|+\|\mathbf{m}_{k}-\bm{\sigma}_{k}-\mathbf{m}_{i}\|$
$\displaystyle\|\mathbf{m}_{k}-\mathbf{m}_{i}\|-\|\mathbf{m}_{k}-\bm{\sigma}_{k}-\mathbf{m}_{i}\|$
$\displaystyle\leq\|\bm{\sigma}_{k}\|.$

This can be combined with

$\displaystyle\|\mathbf{m}_{k}-\bm{\sigma}_{k}-\mathbf{m}_{i}\|$
$\displaystyle\leq\|\mathbf{m}_{k}-\mathbf{m}_{i}\|+\|\bm{\sigma}_{k}\|$
$\displaystyle\|\mathbf{m}_{k}-\bm{\sigma}_{k}-\mathbf{m}_{i}\|-\|\mathbf{m}_{k}-\mathbf{m}_{i}\|$
$\displaystyle\leq\|\bm{\sigma}_{k}\|,$

to obtain

$\Bigl{|}\|\mathbf{m}_{k}-\mathbf{m}_{i}\|-\|\mathbf{m}_{new}-\mathbf{m}_{i}\|\Bigr{|}\leq\|\bm{\sigma}_{k}\|.$
(14)

The upper bound for $|\Gamma_{2}|$ is

$\displaystyle|\Gamma_{2}|=$
$\displaystyle\;\Biggl{|}\sum_{\begin{subarray}{c}i=0\\\
i\notin\\{k,l\\}\end{subarray}}^{c}\mathbf{u}_{k}^{\top}\mathbf{Wu}_{i}(\|\mathbf{m}_{k}-\mathbf{m}_{i}\|-\|\mathbf{m}_{new}-\mathbf{m}_{i}\|)$
$\displaystyle\;+\sum_{\begin{subarray}{c}i=0\\\
i\notin\\{k,l\\}\end{subarray}}^{c}\mathbf{u}_{l}^{\top}\mathbf{Wu}_{i}(\|\mathbf{m}_{l}-\mathbf{m}_{i}\|-\|\mathbf{m}_{new}-\mathbf{m}_{i}\|)+\mathbf{u}_{k}^{\top}\mathbf{Wu}_{l}\|\mathbf{m}_{k}-\mathbf{m}_{l}\|\Biggr{|}$
$\displaystyle\leq$ $\displaystyle\;\sum_{\begin{subarray}{c}i=0\\\
i\notin\\{k,l\\}\end{subarray}}^{c}\mathbf{u}_{k}^{\top}\mathbf{Wu}_{i}\|\bm{\sigma}_{k}\|+\sum_{\begin{subarray}{c}i=0\\\
i\notin\\{k,l\\}\end{subarray}}^{c}\mathbf{u}_{l}^{\top}\mathbf{Wu}_{i}\|\bm{\sigma}_{l}\|+\mathbf{u}_{k}^{\top}\mathbf{Wu}_{l}\|\mathbf{m}_{k}-\mathbf{m}_{l}\|$
$\displaystyle\leq$
$\displaystyle\;\Biggl{(}\mathbf{u}_{k}^{\top}\mathbf{Wu}_{l}+\frac{|C_{l}|}{|C_{k}|+|C_{l}|}\sum_{\begin{subarray}{c}i=0\\\
i\notin\\{k,l\\}\end{subarray}}^{c}\mathbf{u}_{k}^{\top}\mathbf{Wu}_{i}+\frac{|C_{k}|}{|C_{k}|+|C_{l}|}\sum_{\begin{subarray}{c}i=0\\\
i\notin\\{k,l\\}\end{subarray}}^{c}\mathbf{u}_{l}^{\top}\mathbf{Wu}_{i}\Biggr{)}\varepsilon_{f},$

where we made use of the triangle inequality, (13), and (14). From these two
upper bounds, we can conclude that if $\varepsilon_{f}\rightarrow 0$, both
$|\Gamma_{1}|\rightarrow 0$ and $|\Gamma_{2}|\rightarrow 0$. Hence,
$|L(\mathbf{M})-L(\mathbf{M}_{new})|\rightarrow 0$ as
$\varepsilon_{f}\rightarrow 0$.

## References

  * Arya et al. (2019) S. Arya, D. Mount, S. E. Kemp, and G. Jefferis.  _RANN: Fast Nearest Neighbour Search (Wraps ANN Library) Using L2 Metric_ , 2019.  URL https://CRAN.R-project.org/package=RANN.  R package version 2.6.1.
  * Bentley (1975) J. L. Bentley.  Multidimensional Binary Search Trees Used for Associative Searching.  _Communications of the ACM_ , 18(9):509–517, 1975\.
  * Boyd and Vandenberghe (2004) S. P. Boyd and L. Vandenberghe.  _Convex Optimization_.  Cambridge University Press, 2004.
  * Boyd et al. (2011) S. P. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein.  Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers.  _Foundations and Trends® in Machine learning_ , 3(1):1–122, 2011.
  * Bronson (1989) R. Bronson.  _Schaum’s Outline of Theory and Problems of Matrix Operations_.  McGraw-Hill, 1989.
  * Chi and Lange (2015) E. C. Chi and K. L. Lange.  Splitting Methods for Convex Clustering.  _Journal of Computational and Graphical Statistics_ , 24(4):994–1013, 2015.
  * Deng (2012) L. Deng.  The MNIST Database of Handwritten Digit Images for Machine Learning Research.  _IEEE Signal Processing Magazine_ , 29(6):141–142, 2012.
  * Dua and Graff (2019) D. Dua and C. Graff.  UCI Machine Learning Repository, 2019.  URL http://archive.ics.uci.edu/ml.
  * Eddelbuettel (2013) D. Eddelbuettel.  _Seamless R and C ++ Integration with Rcpp_.  Springer, 2013.
  * Eddelbuettel and Balamuta (2017) D. Eddelbuettel and J. J. Balamuta.  Extending R with C++: A Brief Introduction to Rcpp.  _PeerJ Preprints_ , 5:e3188v1, 2017.
  * Eddelbuettel and François (2011) D. Eddelbuettel and R. François.  Rcpp: Seamless R and C++ Integration.  _Journal of Statistical Software_ , 40(8):1–18, 2011.
  * Gan et al. (2007) G. Gan, C. Ma, and J. Wu.  _Data Clustering: Theory, Algorithms, and Applications_.  SIAM, 2007.
  * Gower and Groenen (1991) J. C. Gower and P. J. F. Groenen.  Applications of the Modified Leverrier-Faddeev Algorithm for the Construction of Explicit Matrix Spectral Decompositions and Inverses.  _Utilitas Mathematica_ , 40:51–64, 1991.
  * Groenen et al. (1999) P. J. F. Groenen, W. J. Heiser, and J. J. Meulman.  Global Optimization in Least Squares Multidimensional Scaling by Distance Smoothing.  _Journal of Classification_ , 16(2):225–254, 1999\.
  * Guennebaud and Jacob (2010) G. Guennebaud and B. Jacob.  Eigen v3.  http://eigen.tuxfamily.org, 2010.
  * Havel (1991) T. F. Havel.  An Evaluation of Computational Strategies for Use in the Determination of Protein Structure from Distance Constrains Obtained by Nuclear Magnetic Resonance.  _Progress in Biophysics and Molecular Research_ , 56(1):43–78, 1991.
  * Hocking et al. (2011) T. D. Hocking, A. Joulin, F. Bach, and J.-P. Vert.  Clusterpath: An Algorithm for Clustering Using Convex Fusion Penalties.  In _The 28th International Conference on Machine Learning_ , Bellevue, Washington, 2011.
  * Hunter and Lange (2004) D. R. Hunter and K. L. Lange.  A Tutorial on MM Algorithms.  _The American Statistician_ , 58(1):30–37, 2004\.
  * Jakob et al. (2017) W. Jakob, J. Rhinelander, and D. Moldovan.  pybind11—Seamless operability between C++11 and Python, 2017\.  https://github.com/pybind/pybind11.
  * Lange et al. (2000) K. L. Lange, D. R. Hunter, and I. Yang.  Optimization Transfer Using Surrogate Objective Functions.  _Journal of Computational and Graphical Statistics_ , 9(1):1–20, 2000.
  * Leeuw (1977) J. Leeuw, de.  Applications of Convex Analysis to Multidimensional Scaling.  In J.-R. Barra, F. cois Brodeau, G. Romier, and B. van Cutsem, editors, _Recent Developments in Statistics_ , pages 133–146. North Holland Publishing Company, 1977.
  * Leeuw (1994) J. Leeuw, de.  Block-Relaxation Algorithms in Statistics.  In H.-H. Bock, W. Lenski, and M. M. Richter, editors, _Information Systems and Data Analysis_ , pages 308–324. Springer Berlin Heidelberg, 1994.
  * Leeuw and Heiser (1980) J. Leeuw, de and W. J. Heiser.  Multidimensional Scaling with Restrictions on the Configuration.  _Multivariate Analysis_ , 5(1):501–522, 1980\.
  * Lindsten et al. (2011) F. Lindsten, H. Ohlsson, and L. Ljung.  Just Relax and Come Clustering!: A Convexification of K-Means Clustering.  Technical report, Department of Electrical Engineering, Linköping University, Linköping, Sweden, 2011.
  * MacQueen (1967) J. MacQueen.  Some Methods for Classification and Analysis of Multivariate Observations.  In _Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability_ , pages 281–297, 1967.
  * McInnes et al. (2018) L. McInnes, J. Healy, and J. Melville.  UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction, 2018.  URL https://arxiv.org/abs/1802.03426.
  * Nesterov (2004) Y. Nesterov.  _Introductory Lecture on Convex Optimization: A Basic Course_.  Springer New York, 2004.
  * Ortega and Rheinboldt (1970) J. M. Ortega and W. C. Rheinboldt.  _Iterative Solution of Nonlinear Equations in Several Variables_.  Academic Press, 1970.
  * Pedregosa et al. (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine Learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Pelckmans et al. (2005) K. Pelckmans, J. de Brabanter, J. A. K. Suykens, and B. de Moor.  Convex Clustering Shrinkage.  In _PASCAL Workshop on Statistics and Optimization of Clustering Workshop_ , 2005.
  * Sun et al. (2021) D. Sun, K.-C. Toh, and Y. Yuan.  Convex Clustering: Model, Theoretical Guarantee and Efficient Algorithm.  _Journal of Machine Learning Research_ , 22(9):1–32, 2021.
  * Voß and Eckhardt (1980) H. Voß and U. Eckhardt.  Linear Convergence of Generalized Weiszfeld’s Method.  _Computing_ , 25(3):243–251, 1980.
  * Weiszfeld (1937) E. Weiszfeld.  Sur le Point pour lequel la Somme des Distances de $n$ Points Donnés Est Minimum.  _Tohoku Mathemacial Journal_ , 43:355–386, 1937.
  * Weylandt et al. (2020) M. Weylandt, J. Nagorski, and G. I. Allen.  Dynamic Visualization and Fast Computation for Convex Clustering via Algorithmic Regularization.  _Journal of Computational and Graphical Statistics_ , 29(1):87–96, 2020.
  * Yuan and Lin (2006) M. Yuan and Y. Lin.  Model Selection and Estimation in Regression with Grouped Variables.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 68(1):49–67, 2006.
  * Yuille and Rangarajan (2003) A. L. Yuille and A. Rangarajan.  The Concave-Convex Procedure.  _Neural Computation_ , 15(4):915–936, 2003.