---
abstract: 'We present a novel algorithm, called *Links*, designed to perform online clustering on unit vectors in a high-dimensional Euclidean space. The algorithm is appropriate when it is necessary to cluster data efficiently as it streams in, and is to be contrasted with traditional batch clustering algorithms that have access to all data at once. For example, Links has been successfully applied to embedding vectors generated from face images or voice recordings for the purpose of recognizing people, thereby providing real-time identification during video or audio capture.'
address: |
    ^1^Google Inc., USA ^2^Carnegie Mellon University, USA\
    [ ^1^ { [](mailto:memes@google.com), [](mailto:quanw@google.com), [](mailto:liwan@google.com), [](mailto:elnota@google.com) } [@google.com]{} ^2^ [](mailto:cmdowney@cs.cmu.edu) ]{}
bibliography:
- 'supplement\_refs.bib'
title: 'Links: A High-Dimensional Online Clustering Method'
---

Introduction {#sec:intro}
============

Although a wide selection of clustering methods are available [@everitt2011cluster; @hennig2015handbook], most of them assume concurrent access to all data being clustered. Our interest is in efficiently clustering each datum as it becomes available, for applications that require unsupervised learning in real time.

The Links approach is to estimate the probability distribution of each cluster based on its current constituent vectors, to use those estimates to assign new vectors to clusters, and to update estimated distributions with each added vector. The update step includes fixing past cluster assignments where indicated by taking the additional data into account, although this is primarily to improve the internal model over time, since in typical online usage scenarios, each cluster assignment is provided once, at the time a new vector is made available.

Prior work [@straub2015small] addressing online clustering of unit vectors employs a small-variance approximation and is applied to low-dimensional problems such as segmentation of surface normals in 3D. Our approach is complementary in that it uses a high-dimensional approximation, and has been applied to problems with relatively high variance.

Links has been used to cluster CNN-based FaceNet embeddings [@schroff2015facenet] and LSTM-based voice embeddings [@wan2017ge2e]. The results of the latter experiment are presented in a separate paper [@wang2017speaker]. The current paper focuses on the technical details of the algorithm.

Model {#sec:model}
=====

Generative model for a cluster {#sec:generative}
------------------------------

Let $X=\{\mathbf{x}_i\}$ be a set of unit-length vectors in $\mathbb{R}^{N}$. They are confined to the submanifold $S^{N-1}$, and to determine proximity for the purpose of clustering these vectors, we will use the natural metric on this submanifold, which is simply the angle between vectors: $$\angle(\mathbf{x},\mathbf{x'}) = \arccos(\mathbf{x}\cdot\mathbf{x'}).$$

We address the problem of cluster distributions within this submanifold with the following properties:

1.  Each cluster has a *center* vector $\bm{\upmu}$ and its member vectors $\mathbf{x}$ are generated by a probability density that is isotropic in the sense that it only depends on distance from the center, $\mathbf{x}\sim\rho(\mathbf{x};\bm{\upmu})=f\left(\angle(\mathbf{x},\bm{\upmu})\right)$.

2.  The function $f$ is the same for every cluster, so that probability densities for different clusters are related by isometry.

3.  $f(\theta)$ decreases exponentially with $\theta$; for example, as a Gaussian suitably normalized on $S^{N-1}$: $$f(\theta)\propto \mbox{\Large$e^{-\frac{\theta^2}{2\sigma^2}}$}.$$ This ensures that the distribution is reasonably localized, since the exponential decrease compensates for a polynomial factor in the marginal distribution of $\theta$: $$\rho(\angle(\mathbf{x},\bm{\upmu})=\theta)=\mathcal{A}(\sin\theta)^{N-2}f(\theta)$$ where $\mathcal{A}$ is a constant equal to the hypersurface area of $S^{N-2}$, $$\mathcal{A}=\frac{2\pi^{(N-1)/2}}{\Gamma\left(\frac{N-1}{2}\right)}.$$

4.  The prior distribution $\rho(\bm{\upmu})$ for the center of a cluster $\bm{\upmu}$ is constant on $S^{N-1}$ (no unit vector is preferred).

Estimated distribution {#sec:estimated}
----------------------

Given a set $X=\{\mathbf{x}_i\}$ chosen randomly from the same cluster, but without knowledge of the center of the cluster, we would like to estimate the cluster’s probability distribution. The likelihood of the center value $\bm{\upmu}$ is $$\label{eqn:likelihood}
\begin{split}
L(\bm{\upmu};X) & =\prod_i \rho(\mathbf{x}_i;\bm{\upmu}) \\
& \propto \mbox{\Large$e^{-\frac{\sum_i \left(\angle(\mathbf{x}_i,\bm{\upmu})\right)^2}{2\sigma^2}}$}
\end{split}$$ Since the prior $\rho(\bm{\upmu})$ is constant, the posterior $\rho(\bm{\upmu}|X)$ is also proportional to the expression in equation \[eqn:likelihood\]. The maximum likelihood (and maximum a posteriori) center is therefore $$\bm{\hat{\upmu}}=\operatorname*{argmin}_{\bm{\upmu}}\sum_i \big(\angle(\mathbf{x}_i,\bm{\upmu})\big)^2$$ which is the same as the centroid of the vectors $\{\mathbf{x}_i\}$ as defined for a hypersphere according to [@buss2001spherical]. The estimated probability distribution for the cluster is $$\begin{split}
\hat{\rho}(\mathbf{x};X) & \propto \int_{S^{N-1}}  L(\bm{\upmu};X)\rho(\mathbf{x};\bm{\upmu})d\bm{\upmu} \\
& \propto \int_{S^{N-1}}  \mbox{\Large$e^{-\frac{\left(\angle(\mathbf{x},\bm{\upmu})\right)^2+\sum_i \left(\angle(\mathbf{x}_i,\bm{\upmu})\right)^2}{2\sigma^2}}$} d\bm{\upmu}.
\end{split}$$ The probability that a new vector $\mathbf{x}$ belongs to the same cluster can then be estimated as the cumulative amount $$\label{eqn:cumulative}
\int_{\{\mathbf{y}\in S^{N-1}\,\mid\,\hat{\rho}(\mathbf{y};X)\geq\hat{\rho}(\mathbf{x};X)\}} \hat{\rho}(\mathbf{y};X)d\mathbf{y}.$$

High-dimensional approximation {#sec:approximation}
------------------------------

Our primary interest is in problems with relatively large $N$. For example, our typical embedding vectors have $N\ge128$. For large enough $N$, the following are true:

Lemma 1

:   Two randomly chosen vectors $\mathbf{x},\mathbf{x'}$ are almost always almost perpendicular, i.e., $$P(\mathbf{x}\cdot\mathbf{x'}>\delta)<\epsilon$$ for some positive numbers $\delta\ll1$ and $\epsilon\ll1$.

Lemma 2

:   The angle $\theta$ between a cluster center and a random vector from that cluster is almost always almost equal to a global constant $\theta_c$, i.e., $$\label{key}
      P\left(|\theta-\theta_c|>\delta\right)<\epsilon$$ for some positive numbers $\delta\ll\pi$ and $\epsilon\ll1$.

Lemma 3

:   Given two randomly chosen vectors from a cluster with center $\bm{\upmu}$, their components perpendicular to $\bm{\upmu}$ will almost always be almost perpendicular to each other, i.e., $$P\Big(\big(\mathbf{x}-(\mathbf{x}\cdot\bm{\upmu})\bm{\upmu}\big)\cdot\big(\mathbf{x'}-(\mathbf{x'}\cdot\bm{\upmu})\bm{\upmu}\big)>\delta\Big)<\epsilon$$ for some positive numbers $\delta\ll1$ and $\epsilon\ll1$.

To assess whether to add a new vector $\mathbf{x}$ to an existing cluster known to include the $k$ vectors $\left\lbrace\mathbf{x}_i\right\rbrace_{i=1}^k$, we determine a threshold $\mathbf{x}\cdot\bm{\hat{\upmu}}\geq s(k)$ on the cosine similarity between the new vector and the centroid $\bm{\hat{\upmu}}$ of the existing vectors. Using the approximation in lemmas 2 and 3, and assuming $N\gg k$, we can compute vector components in an orthonormal basis including $\bm{\upmu}$, $\frac{1}{\sin\theta_c}(\mathbf{x}-\cos\theta_c\bm{\upmu})$ and $\left\lbrace \frac{1}{\sin\theta_c}\left(\mathbf{x}_i-\cos\theta_c\bm{\upmu}\right)\right\rbrace_{i=1}^k$. This yields

$$\bm{\hat{\upmu}}=\frac{1}{\sqrt{k^2\cos^2\theta_c+k\sin^2\theta_c}}\sum_{i=1}^{k}\mathbf{x}_i$$

and a threshold of $$\label{eqn:single_threshold}
\begin{split}
s(k) & =\frac{k\cos^2\theta_c}{\sqrt{k^2\cos^2\theta_c+k\sin^2\theta_c}} \\
& =\frac{T_c^2}{\sqrt{\frac{1}{k}+\left(1-\frac{1}{k}\right)T_c^2}}
\end{split}$$ where $T_c=\cos\theta_c$, which we call the *cluster similarity threshold*.

Note that $$\lim_{k\rightarrow\infty}\bm{\hat{\upmu}}=\bm{\upmu}$$ and $$\lim_{k\rightarrow\infty}s(k)=T_c,$$ which confirms that as we accumulate more vectors in a given cluster, the center and cosine similarity threshold of the *estimated* distribution approach the center and cosine similarity threshold of the *generative* distribution (i.e., the estimate improves). Since $s(k)$ is a strictly increasing function of $k$, the variance of the estimated distribution decreases with $k$.

Similarly, to assess whether two clusters are the same, we determine a threshold on the cosine similarity between their centroids $\bm{\upmu}_c\cdot\bm{\upmu}_c'\geq s(k,k')$ where, for $N\gg k$ and $N\gg k'$, $$\label{eqn:multi_threshold}
\begin{split}
s(k,k') & =\frac{1}{\sqrt{\left( 1+\frac{1}{k}\tan^2\theta_c\right) \left( 1+\frac{1}{k'}\tan^2\theta_c\right)}} \\
& =\frac{1}{\sqrt{\bigg( 1+\frac{1}{k}\left(\frac{1}{T_c^2}-1\right)\bigg)\bigg(1+\frac{1}{k'}\left(\frac{1}{T_c^2}-1\right)\bigg)}}.
\end{split}$$

Note that equation \[eqn:single\_threshold\] is the special case with $k'=1$, $$s(k,1)=s(k),$$ and $$\label{eqn:multi_limit}
\lim_{k,k'\rightarrow\infty}s(k,k')=1.$$ The latter confirms that the centers estimated from the two sets of cluster points converge.

Algorithm {#sec:algorithm}
=========

Online clustering {#sec:online}
-----------------

Each new input vector is assigned to a cluster as soon as it is produced, with no knowledge of future vectors and no backtracking. A unique ID for that cluster is returned. The clusterer keeps statistical information about the vectors received so far. Although it cannot change a previous answer, it can change the internal representation of cluster statistics, such as improvements to estimated distributions as well as cluster splits and merges when indicated by new information.

Internal representation {#sec:internal}
-----------------------

The Links algorithm’s internal representation is a two-level hierarchy: clusters are collections of subclusters, and subclusters are collections of input vectors. The subclusters are represented as nodes in a graph whose edges join ‘nearby’ nodes (meaning subclusters that likely belong to the same cluster given the data so far), and clusters are defined as connected components of the graph. Whereas subclusters are indivisible, clusters can become split along graph edges in response to changes in subcluster estimated probability distributions as new data is added. Alternatively, subclusters joined by an edge can become merged in response to changes.

The reasons for maintaining this two-level hierarchy (rather than, say, an arbitrary number of levels) are efficiency and practicality. It is efficient because the algorithm scales with number of subclusters rather than number of vectors. It is practical because the key cluster substructure that can affect future cluster IDs is the set of potential split points.

Assessing cluster membership {#sec:assessing}
----------------------------

When a new vector $\mathbf{x}$ is available, compute its cosine similarity to each subcluster centroid $\bm{\hat{\upmu}}_j$, and add it to the most-similar subcluster if the similarity is above a fixed threshold $T_s$. In other words, let $$J = \operatorname*{argmax}_j\{\mathbf{x}\cdot\bm{\hat{\upmu}}_j\}.$$ If $$\label{eqn:subcluster}
\mathbf{x}\cdot\bm{\hat{\upmu}}_J\geq T_s$$ then add $\mathbf{x}$ to subcluster $J$. $T_s$, called the *subcluster similarity threshold*, is a hyperparameter determining the granularity of cluster substructure appropriate for the data.

If inequality \[eqn:subcluster\] does not hold, then start a new subcluster containing just $\mathbf{x}$. Next, use the estimated probability distribution of subcluster $J$ to determine whether to include the new subcluster in the same cluster as $J$, by thresholding the cumulative probability in expression \[eqn:cumulative\]. In the high-dimensional approximation, this means the subcluster is included in the cluster whenever $$\label{eqn:cluster}
\mathbf{x}\cdot\bm{\hat{\upmu}}_J\geq s(k_J)$$ where $k_J$ is the number of vectors in the subcluster $J$. To a first approximation, $s(k)$ is as given in equation \[eqn:single\_threshold\]. This will be further refined in section \[sec:anisotropy\]. If inequality \[eqn:cluster\] does hold, then add an edge to the graph joining the new subcluster to subcluster $J$.

Updating clusters {#sec:updating}
-----------------

When a new vector is added to an existing subcluster, the subcluster’s centroid may change. If this brings it within the subcluster similarity threshold of the centroid of another subcluster currently joined to the first by an edge, then the two are merged. In other words, if $\bm{\hat{\upmu}}_i\cdot\bm{\hat{\upmu}}_j\geq T_s$, then nodes $i$ and $j$ are replaced with a single node containing the vectors of both, and with the edge connections of both. Since the merging process also results in a new subcluster centroid, this check is continued recursively on affected subclusters.

Next, the edges joining affected nodes are checked for validity. The edge joining subclusters $i$ and $j$ is removed if the following does not continue to hold: $$\label{eqn:no_split}
\bm{\hat{\upmu}}_i\cdot\bm{\hat{\upmu}}_j\geq s(k_i,k_j)$$ where $s(k_i,k_j)$ is approximately as given in equation \[eqn:multi\_threshold\], but with improvements to follow in section \[sec:anisotropy\]. After severing a cluster in two by removing an edge, an attempt is made to re-join the two parts by adding an edge from the affected node to a new partner node that does satisfy inequality \[eqn:no\_split\]. If no such partner is found, then the cluster remains permanently split.

Anisotropy {#sec:anisotropy}
----------

Equations \[eqn:single\_threshold\] and \[eqn:multi\_threshold\] were used to determine thresholds for membership in the same cluster as a given subcluster, effectively treating the subcluster’s members as randomly chosen from the cluster and not correlated with each other. If one were to properly take into account intra-subcluster correlations, then one consequence is that the limit in equation \[eqn:multi\_limit\] would be reduced to a positive number $T_p<1$, which we call the *pair similarity maximum*, $$\label{eqn:anisotropic_limit}
\lim_{k,k'\rightarrow\infty}s(k,k')=T_p,$$ whereas the value of $s(1,1)$, which is $T_c^2$, would remain unchanged. Any implicit anisotropy in the cluster distribution, such as an elongation along a preferred axis, will further reduce the value of $T_p$ without changing $s(1,1)$. A simple though approximate way to incorporate these adjustments into the algorithm is to replace $s(k,k')$ and $s(k)$ by the following interpolated versions: $$\label{eqn:anisotropic_multi}
\tilde{s}(k,k')=T_c^2+\frac{T_p-T_c^2}{1-T_c^2}\left(s(k,k')-T_c^2\right)$$ $$\label{eqn:anisotropic_single}
\tilde{s}(k)=\tilde{s}(k,1).$$

Hyperparameter Tuning {#sec:hyperparameter}
---------------------

The similarity thresholds $T_c$, $T_s$ and $T_p$ need to be tuned to best represent the data source. This is done by manually labeling a dataset with cluster IDs, running the clusterer on the data, and adjusting hyperparameters to improve the accuracy of the output cluster IDs. Accuracy is simply fraction of correct IDs. Prior to evaluation, the Hungarian algorithm [@kuhn1955hungarian] is used to map a subset of output cluster IDs bijectively to a subset of ground truth cluster IDs in such a way that produces the best possible accuracy. For some applications an alternate objective has been used; for example, one that gives different weights for conflating IDs vs. fracturing IDs, to reflect the seriousness of each type of error in practise.

Acknowledgements {#sec:ack}
================

The authors would like to thank Dr. Brian Budge and Dr. Navid Shiee for help with APIs and evaluation frameworks used in the implementation of the Links algorithm.
