Scalable K-Medoids via True Error Bound and Familywise Bandits

TL;DR

This paper introduces a new theoretical framework for true K-Medoid error, proposes a scalable distributed algorithm MCPAM, and demonstrates its efficiency and accuracy on large semi-metric datasets, including one with 1 billion points.

Contribution

It formalizes the true K-Medoid error, provides convergence analysis, and develops MCPAM, a scalable distributed algorithm for large-scale clustering.

Findings

MCPAM achieves significant computational savings.

Error decreases at a rate of Θ(1/n^{2/3}).

Successfully clusters 1 billion points on semi-metric spaces.

Abstract

K-Medoids(KM) is a standard clustering method, used extensively on semi-metric data.Error analyses of KM have traditionally used an in-sample notion of error,which can be far from the true error and suffer from generalization gap. We formalize the true K-Medoid error based on the underlying data distribution.We decompose the true error into fundamental statistical problems of: minimum estimation (ME) and minimum mean estimation (MME). We provide a convergence result for MME. We show $\errMME$ decreases no slower than $\Theta(\frac{1}{n^{\frac{2}{3}}})$, where $n$ is a measure of sample size. Inspired by this bound, we propose a computationally efficient, distributed KM algorithm namely MCPAM. MCPAM has expected runtime $\mathcal{O}(km)$,where $k$ is the number of medoids and $m$ is number of samples. MCPAM provides massive computational savings for a small tradeoff in accuracy. We verify the quality and scaling properties of MCPAM on various datasets. And achieve the hitherto unachieved feat of calculating the KM of 1 billion points on semi-metric spaces.

Tables

Table 1. Table 1: Various K-medoid algorithms as ν ^ ^ ^ ^ 𝜈 \hat{\hat{\nu}}

Estimator	$R$ (iid samples)	$S$	$T_i$
EXHAUSTIVE	$X_1,\dots ,X_n$	$R^K$	$R$
PAM	$X_1,\dots ,X_n$	Walk through $R^K$	$R$
IPAM	$X_1,\dots ,X_{n+mn}$	Walk through $R_1^K$ ($R_1={\left\{X_i\right\}}_{i=1}^n$)	${\left\{X_{ni+j}\right\}}_{j=1}^n$, $i\in [1,m]$
CLARANS	$X_1,\dots ,X_n$	Walk through $R^K$	R
RAND	$X_1,\dots ,X_m$	$R$	random subset of $R$

Table 2. Table 2: Scaling of MCPAM on BillionOne dataset. Due to the large size of the dataset, we only run with 12 workers.

$K$	Num Workers	Num MCPAM Iters	Runtime (s)	Memory (MB)
1	12	2	5800	4030
2	12	2	5870	4030
5	12	3	6310	4030
10	12	3	7300	4030
15	12	6	8510	4030
20	12	8	9830	4030
30	12	10	11800	4030
50	12	5	15700	4030
100	12	1	25400	4030

Equations

$$\mathrm{Ecc}(\bar{x})$$

$$\hat{\mathrm{E}}\mathrm{cc}(\bar{x},T)$$

$$\hat{\hat{\nu }}(S,T_1,\cdots ,T_m)$$

$${\text{err}}_1$$

$${\text{err}}_3$$

$${\text{err}}_3:={\text{err}}_{\text{MME}}+{\text{err}}_{\text{MME}}{\text{err}}_{\text{ME}}+{\text{err}}_{\text{ME}}$$

$$E{\text{err}}_1$$

$${\bar{X}}_i$$

$$X_{ij}$$

$$\hat{\mathrm{E}}\mathrm{cc}({\bar{X}}_i,\{X_{ij}{\}}_{j=1}^m)$$

$$||{\bar{X}}_i-X_{ij}{||}_2^2\sim {\chi }^2(1,X_i^2)$$

$$\kappa (\lambda )=3+12\frac{(1+4\lambda )}{(1+2\lambda {)}^2}$$

$$\kappa (\lambda )\le \kappa (0)=15$$

$$\mathrm{Var}(\hat{\mathrm{E}}\mathrm{cc}(\bar{x},\{y_j\}))$$

$${\hat{\mathrm{E}}\mathrm{cc}}_{\text{lo}}(\bar{x})$$

$${\text{err}}_3:={\text{err}}_1+{\text{err}}_1{\text{err}}_2+{\text{err}}_2$$

$${\text{err}}_3$$

$$=\frac{E\left|\mathrm{Ecc}(\hat{\hat{\nu }})-\mathrm{Ecc}(\hat{\nu })+\mathrm{Ecc}(\hat{\nu })-\mathrm{Ecc}(\nu )\right|}{E(\mathrm{Ecc}(\nu ))}$$

$$=\frac{E\left|\mathrm{Ecc}(\hat{\hat{\nu }})-\mathrm{Ecc}(\hat{\nu })\right|+\left|\mathrm{Ecc}(\hat{\nu })-\mathrm{Ecc}(\nu )\right|}{E(\mathrm{Ecc}(\nu ))}$$

$$={\text{err}}_1\frac{E(\mathrm{Ecc}(\hat{\nu }))}{E(\mathrm{Ecc}(\nu ))}+{\text{err}}_2$$

$${\text{err}}_3$$

$${\text{err}}_3:={\text{err}}_{\text{MME}}+{\text{err}}_{\text{MME}}{\text{err}}_{\text{ME}}+{\text{err}}_{\text{ME}}$$

$${\text{err}}_1$$

$${\text{err}}_2$$

$$E{\text{err}}_1$$

$$E\mathrm{Ecc}(\hat{\hat{\nu }})=\underset{S}{E}\underset{\{T_i{\}}_{i=1}^m}{E}\mathrm{Ecc}(\hat{\hat{\nu }})$$

$$\underset{\{T_i{\}}_{i=1}^m}{E}\mathrm{Ecc}(\hat{\hat{\nu }})=\underset{\{T_i{\}}_{i=1}^m}{E}\mathrm{Ecc}(\hat{\hat{\nu }}(S,T_1,\dots ,T_m))|S$$

$$\underset{\{T_i{\}}_{i=1}^m}{E}\mathrm{Ecc}(\hat{\hat{\nu }})|S\le (1+\frac{p}{100})\mathrm{Ecc}(\hat{\nu })|S$$

$$E|\mathrm{Ecc}(\hat{\hat{\nu }})-\mathrm{Ecc}(\hat{\nu })|$$

$$\le \underset{S}{E}(1+\frac{p}{100})\mathrm{Ecc}(\hat{\nu })-\underset{S}{E}\mathrm{Ecc}(\hat{\nu })$$

$$=\frac{p}{100}E\mathrm{Ecc}(\hat{\nu })$$

$${\hat{i}}_{\min }:=\underset{i}{\mathrm{argmin}}{\hat{\mu }}_i$$

$$ϵ({\mu }_1,\dots ,{\mu }_m):=\underset{E_{ij}}{E}\frac{{\mu }_{{\hat{i}}_{\min }}-{\mu }_{i_{\min }}}{{\mu }_{i_{\min }}}$$

$${\delta }_{i:m}$$

Full text

1 Introduction

K-medoids [26] is an extremely general clustering method. Given a data sample, it finds K data points that are the centers of K clusters. The K points together are called a K-medoid. K-medoids is extensively used on semi-metric data, that doesn’t necessarily respect the triangle inequality. Eg internet RTT [12], RNA-Seq analysis [35], recommender systems [28] where Euclidean metric doesn’t apply. See [26, 4] for good introductions. Other K-medoid algorithms include [34, 25, 42, 9, 37, 33].

Traditionally, K-medoids has used an in-sample notion of error [9, 4]. For example in the $K=1$ case, the most central point of the sample is taken as the reference/true point to calculate error. This is analogous 111 The analogy can be formalized by recognizing the medoid as a generalization of the mean to using the sample mean $\hat{\mu}$ as the reference to calculate error, instead of the true mean $\mu$. However $\hat{\mu}$ is a fluctuating quantity and $|\hat{\mu}-\mu|$ can be quite large. To reiterate, any dataset is a random, limited sampling of an underlying distribution $\mathcal{D}$. All sample quantities are fluctuating, noisy approximations to underlying distributional quantities, and should be treated as such. Again, from the perspective of mixture models, using the true mean $\mu$ as the reference point is standard practice [44, 40, 29, 36, 6]. The use of in-sample error can lead to generalization gap issues [20] Our contributions include:

• •

Formalizing the true K-medoid error

• •

Fundamental insights into the $K$-medoids problem, by showing that it decomposes into two basic statistical problems: minimum estimation (ME) & minimum mean estimation (MME).

• •

A fundamental convergence result for MME. We show $\text{err}_{\text{MME}}$ decreases no slower than $\Theta(\frac{1}{n^{\frac{2}{3}}})$. Where $n$ is a measure of sample size.

• •

Inspired by above analysis, a new extremely scalable K-medoid algorithm (MCPAM). MCPAM has average runtime $\mathcal{O}(km)$, where $k$ is number of medoids and $m$ is number of samples. This makes it the first linear 1-medoid algorithm (expected runtime).

We provide detailed comparisons of these contributions to prior literature in sections 3.4,4.1.

2 Problem Formulation

Let $(\mathbb{X},d,\mathcal{D}_{\mathbb{X}})$ be a semi-metric space equipped with a probability distribution $\mathcal{D}_{\mathbb{X}}$. Let $\mathbb{X}^{K}=\mathbb{X}\times\cdots\times\mathbb{X}$ be the cartesian product of $k$ copies of $\mathbb{X}$. An element $\bar{x}\in\mathbb{X}^{K}$ is a K-tuple and its $l^{\text{th}}$ entry is denoted $\bar{x}[l]$. The distance $\Delta$ from a point $x\in\mathbb{X}$ to a K-tuple $\bar{x}\in\mathbb{X}^{K}$ is the minimum of the $k$ componentwise distances: $\Delta(x,\bar{x})=\min_{l=1}^{k}d(x,\bar{x}[l])$. Intuitively, $\bar{x}$ represents $K$ cluster centers and $\Delta$ is the distance to the nearest center. This is the standard distance for K-medoids [26]. The probabilistic setting has been explored in [33], but only from the perspective of runtime calculation, not for error calculation.

Definition 1** (Eccentricity and Medoid).**

The average distance of a K-tuple $\bar{x}$ to the points in $\mathbb{X}$ is the eccentricity $\mathrm{Ecc}$. A k-medoid $\nu$ is a K-tuple that has minimum eccentricity.

[TABLE]

Eccentricity is inverse centrality. We now develop estimators for $\nu$. Let $\bar{X}$ be a $\mathbb{X}^{K}$ valued random variable. Let $S=\{\bar{X}_{i}\}_{i=1}^{m}$ be $m$ random observations from $\mathbb{X}^{K}$.

Definition 2** (True Sample Medoid).**

The true sample medoid $\hat{\nu}$ is a minimizer of $\mathrm{Ecc}$, over the K-tuples in $S$: $\hat{\nu}(S):=\operatorname*{arg\,min}_{\bar{X}_{i}\in S}\mathrm{Ecc}(\bar{X}_{i})$

It is widely recognized [4] that K-medoids on $\mathbb{R}$ with k = 1 and the $L_{1}$ metric is the median. We will use this in a running example to illustrate various concepts.

Example 1** (1-medoid on $\mathbb{R}$).**

Let $\mathbb{X}=\mathbb{R}$, $d$ be $L_{1}$ metric, $\mathcal{D}_{\mathbb{X}}=\mathcal{N}(\mu=100,\sigma=100)$ & $K=1$. Now $\mathbb{X}^{K}=\mathbb{R}$ & $\mathrm{Ecc}(\bar{x})=$ $\int_{-\infty}^{\infty}|x-\bar{x}|\frac{1}{\sqrt{2*\pi*100^{2}}}\exp{-\frac{(x-100)^{2}}{2*100^{2}}}dx$. This gives $\nu=100$ (the median). Given a sample of $\mathbb{X}^{K}$, $S=\{90,\ 170,\ 60,\ 200,\ 190\}$, $\mathrm{Ecc}(90)\simeq 80.187$. This is the minimum on $S$, so $\hat{\nu}=90$.

But $\mathrm{Ecc}$ is not computable in practice. We need to approximate it. Let $T=\{X_{j}\}_{j=1}^{n}$ be a random sample of size $n$ from $\mathbb{X}$. Let us have $m$ such random samples $T_{i}=\{X_{ij}\}_{j=1}^{n}$ $i\in[1,m]$. Let all $X_{ij}$ follow distribution $\mathcal{D}_{\mathbb{X}}$.

Definition 3** (Sample Eccentricity and Sample Medoid).**

The average distance of a K-tuple $\bar{x}$ to a random sample $T=\{X_{j}\}_{j=1}^{n}$ of $\mathbb{X}$ is the sample eccentricity $\mathrm{\widehat{E}cc}$. The sample K-medoid $\hat{\hat{\nu}}$ is a K-tuple from $S$ that has minimum sample eccentricity:

[TABLE]

$\hat{\nu}$ has one level of approximation to $\nu$, namely the use of $S$ as a proxy for $\mathbb{X}^{K}$. Whereas $\hat{\hat{\nu}}$ has two levels of approximation to $\nu$, the additional one being the use of $T$ as a proxy for $\mathbb{X}$.

Example 2** (1-medoid on $\mathbb{R}$ (cont’d)).**

Consider the setting of example 1. Let $T=S$. Then $\mathrm{\widehat{E}cc}(\bar{x})=\frac{1}{5}\sum_{x\in T}|x-\bar{x}|$ and $\hat{\hat{\nu}}=170$. With two more data points $\{90,\ 170,\ 60,\ 200,\ 190,\ -10,\ 150\}$, $\hat{\hat{\nu}}=150$. $\hat{\nu}$ and $\nu$ do not change.

The $S$, $m$, $T_{i}$ and $n$ are central to $\hat{\hat{\nu}}$. We will reuse them throughout the paper, so we reiterate them as a formal definition.

Definition 4** ($\mathbb{X}^{K}$ Sample: $S$).**

$S=\{\bar{X}_{i}\}_{i=1}^{m}$* is a sample of size $m$ from $\mathbb{X}^{K}$. Each $\bar{X}_{i}$ is a $\mathbb{X}^{K}$ valued random variable.*

Definition 5** ($\mathbb{X}$ Sample: $T$).**

For each $\bar{X}_{i}\in S$, we have a sample $T_{i}$ of size $n$ from $\mathbb{X}$. $T_{i}=\{X_{ij}\}_{j=1}^{n}$ will be used to estimate eccentricity of $\bar{X}_{i}$. $X_{ij}$ is a $\mathcal{D}_{\mathbb{X}}$ distributed random variable.

We now express a number of existing K-medoid algorithms as sample medoids $\hat{\hat{\nu}}(S,T_{1},\ldots,T_{m})$ by appropriate choice of $S,T_{i}$. Most existing algorithms derive the $S$ and $T_{i}$ samples from a common iid sample $R$ of $\mathbb{X}$. For instance in PAM [26], $R=\{X_{1},\ldots,X_{n}\}$ is $n$ iid samples from $\mathbb{X}$. The $T_{i}$ are all equal to $R$, i.e. $T_{i}=\{X_{ij}\}_{j=1}^{n}=\{X_{j}\}_{j=1}^{n}=R$. $S$ is a subset of $R^{K}=R\times\cdot\times R$. In more detail, $S$ is constructed by picking a $\bar{X}_{\text{current}}$ from $R^{K}$ at random and then:

Calculating $\mathrm{\widehat{E}cc}$ for all single swap 222 A single swap neighbour of a $\bar{X}$ is got by swapping exactly one entry of $\bar{X}$ with another $X_{i}\in R$

neighbours of $\bar{X}_{\text{current}}$ 2. 2.

Setting $\bar{X}_{\text{current}}$ to neighbour with lowest $\mathrm{\widehat{E}cc}$ 3. 3.

Repeating from step 1 until no further decrease in $\mathrm{\widehat{E}cc}$

$S$ is all the $\bar{X}$ for which $\mathrm{\widehat{E}cc}$ is calculated (and minimized over). $S$ is a function of $R^{K}$. Since $R^{K}$ is a random sample of $\mathbb{X}^{K}$, it follows that $S$ is too. Finally note $m=|S|\leq|R^{K}|=n^{K}$. We also describe IPAM, an iid version of PAM. Here, $S$ is derived from the first $n$ samples of $R$ in the same way. However for each $\mathrm{\widehat{E}cc}$ calculation in step 1, a fresh batch of $n$ samples are taken.

CLARA [25] is essentially PAM, but $n$ is chosen to be quite small. Other algorithms such as CLARANS [34] and RAND [9], are also expressed as $\hat{\hat{\nu}}$ in table 1. Algorithms TOPRANK [37], trimed [33] and meddit [4] are closely related to our framework. They solve ME problem via an exhaustive exploration of $S$. They solve the MME problem via different estimators, for example trimed uses the triangle inequality to limit $n$ for certain $\mathrm{\widehat{E}cc}(\bar{X}_{i})$ evaluation, etc. We do not express iid versions of algorithms besides PAM. But the modification is the same: for each $\mathrm{\widehat{E}cc}$ evaluation, use a fresh sample from $\mathbb{X}$.

2.1 Issues with Existing K-medoid Error; Definition of True k-Medoid Error

The literature has traditionally focused on bounding the error relative to the EXHAUSTIVE $\hat{\hat{\nu}}$ (see [9, 42]), in some cases even to zero (see [37, 33, 4]). This is important and valuable, but uses a sample quantity as reference. Using $\hat{\hat{\nu}}$ as reference is akin to using sample mean (resp. median) as the reference for mean (resp. median) estimation. $\hat{\hat{\nu}}$ is a fluctuating sample quantity and can be quite far from the true center of the distribution. For the K = 1 case, EXHAUSTIVE $\hat{\hat{\nu}}$ is the sample median and the fluctuation of $\hat{\hat{\nu}}$ is seen in examples 1 & 2. The median example is quite pertinent, since the random nature of sample median is widely recognized and various confidence intervals and estimates of error to true median have been developed for the sample median [32, 7]. In the extreme case, $\hat{\hat{\nu}}$ may not be a consistent estimator of $\nu$ and may never converge. Taking an example from mean estimation, we note that the sample mean of the Cauchy distribution never converges.

From a ML perspective, in-sample errors suffer from generalization gap. This has been studied in the context of clustering in [20]. For the same reason it is standard practice in mixture modeling [44, 40, 29, 36, 6]. to use the true centers as the references. Errors using $\nu$ as reference will address all these gaps. We now define error with respect to $\nu$.

Definition 6** (K-Medoid Errors).**

We have two partial notions of errors:

[TABLE]

And one true error

[TABLE]

We study $\text{err}_{1}$ and $\text{err}_{2}$, as they provide a crucial decomposition of $\text{err}_{3}$ into canonical problems.

3 Results on K-Medoids Errors

The chief results of this section:

• •

Decomposition of $\text{err}_{3}$ into two canonical problems

• •

A new fundamental convergence result for one of these, namely the minimum mean estimation (MME) problem.

3.1 $\text{err}_{3}$ Decomposition

We start by defining two fundamental statistical problems 7 8. Then we decompose $\text{err}_{3}$ into these problems.

Definition 7** (Minimum Mean Estimation (MME)).**

Let $m$ distributions $\{\mathcal{E}_{i}\}_{i=1}^{m}$ be given. Let $n$ samples from each distribution be given, $E=\{E_{ij}\},i\in[1,m],j\in[1,n]$. We want to identify the distribution with minimum mean.

Let $\mu_{i}$ be the mean of the $i$th distribution. Let $i_{\min}$ be the index of the minimum mean distribution. Let $\hat{i}_{\min}(E)$ be an estimate. Then define the MME error as $\text{err}_{\text{MME}}:=|\mu_{\hat{i}_{\min}}-\mu_{i_{\min}}|/|\mu_{i_{\min}}|$

Definition 8** (Minimum Estimation).**

Given $m$ samples from a distribution, estimate the minimum of the distribution $x_{\min}$ using the min of the sample $\hat{x}_{\min}$. Then define the ME error as $\text{err}_{\text{ME}}:=|\hat{x}_{\min}-x_{\min}|/|x_{\min}|$

Theorem 1.

[TABLE]

The structuring of the true error $\text{err}_{3}$ (proof at appendix A, theorem 1) into these canonical problems holds with great generality. This is one of our chief contributions, as it reveals the internal workings of K-medoids. Essentially K-medoids has been decomposed into a ME problem followed by a MME problem.

The ME problem is well studied in the Extreme Value Theory literature [5, 41, 27]. A central theorem in EVT is the Fisher-Tippet-Gnedenko (FTG) theorem providing asymptotic distributions for the sample minimum. This reveals a deep connection between the K-medoids problem and the dynamics of the FTG theorem. The various K-medoid strategies are now re-interpreted as variance reduction strategies for sample min. The path is now open to leverage the extensive EVT literature for better min estimation strategies.

The MME problem is a novel continuous variant of the discrete best arm identification problem found in [16, 4]. Formulation of the MME problem as a stand-alone problem will again help craft novel K-medoid strategies by solving MME in isolation. Finally, examining the tradeoff between ME and MME components of the error will be crucial in designing optimal K-medoid algorithms. An initial version of this has been successfully done in section 4 to develop the MCPAM algorithm.

3.2 Error Bounds

Before we develop bounds on errors, we develop a model for distance distributions. Given a R.V $X$ with distribution $\mathcal{D}_{\mathbb{X}}$. We have a family of random variables $\Delta(X,\bar{x})$, indexed by $\bar{x}\in\mathbb{X}^{K}$. This is a family of distance distributions. We have developed a very general model, called the power-variance family to encode such families.

Definition 9**.**

Consider a family $\mathcal{F}$ of distributions $\mathcal{E}$ parametrized by the mean $\mu$: $\mathcal{F}=\{\mathcal{E}(\mu)|\mu\in\mathbb{R}^{>0}\}$ Such that the variance is a power function of the mean: $\sigma^{2}(\mu):=\alpha\mu^{\beta}+k$ The parameters $\alpha,\beta,\gamma,k$ are constants for the entire family $\mathcal{F}$ and satisfy: $\beta\geq 0,\ \alpha>0,\ \mu\geq\gamma>0,\ \alpha\gamma^{\beta}\geq-k$

Definition 10**.**

Consider a triple $(\mathbb{X},d,\mathcal{D}_{\mathbb{X}})$. Let $X$ be a R.V with distribution $\mathcal{D}_{\mathbb{X}}$. We then have the family of random variables $\mathcal{F}=\{\Delta(X,\bar{x})|\bar{x}\in\mathbb{X}^{K}\}$. If $\mathcal{F}$ is a power-variance family as per definition 9, then the triple $(\mathbb{X},d,\mathcal{D}_{\mathbb{X}})$ is termed power-variance compatible.

The inequality $\mu\geq\gamma>0$ models the strict positivity of distances. The inequality $\alpha\gamma^{\beta}\geq-k$ is needed to ensure positivity of $\sigma^{2}$. This model is quite general and permits us to easily encode distance distributions. As a concrete example of definitions 9 & 10, consider example 3.

To calculate $\operatorname*{\mathbb{E}}\text{err}_{\text{MME}}$ requires calculating probabilities of the form:

$\mathbb{P}(\hat{\mu}_{i}<\hat{\mu}_{1}\cdots\hat{\mu}_{(i-1)},\hat{\mu}_{(i+1)}\cdots\hat{\mu}_{m})$. Where $\hat{\mu}_{i}$ is the sample mean of $\mathcal{E}(\mu_{i})\in\mathcal{F}$. Even for $m=3$ and Gaussian $\mathcal{E}$ this calculation requires the integration of a complicated function of a polynomial of $\Phi$ (standard normal cdf) and is intractable. One of our chief theorems is the upper bound of expected MME error (theorem 5 from appendix D.2). In simplified form:

Theorem 5.

Let $\mathcal{F}$ be a power-variance distribution family, having continuous cdfs and uniformly bounded kurtosis333 kurtosis (or Pearson kurtosis) is the fourth standardized moment $\operatorname*{\mathbb{E}}\left[\frac{X-\mu}{\sigma}\right]^{4}$. $\leq\kappa^{\mathrm{UB}}$. Let $\text{err}_{\text{MME}}$ be as in definition 7 and let the samples from the distributions be independent. Let $0<p\leq p_{\max},m>m_{\min},\beta<\frac{4}{3},n=\left\lceil Cm^{1/2}/p^{3/2}\right\rceil$ . Where $p_{\max}$ & $m_{\min}$ are constants defined in defns 22 & 23. And $C$ is a constant for a given $\mathcal{F}$. Then: $\text{err}_{\text{MME}}<\frac{p}{100}$

In other words, under general conditions, we can control the MME error to arbitrarily small tolerance $p$ with a $n=\Theta(m^{1/2}/p^{3/2})$. A more general version is theorem 4 in appendix D.2. Under general conditions it gives the bound: $\text{err}_{\text{MME}}<\Theta(m^{1/3}/n^{2/3})$ These are fundamental results for MME and one of our chief contributions. They are the analog of the standard square root rate of convergence of Monte Carlo algorithms. Unlike the asympototic CLT based convergence our result is an exact upper bound when its conditions are satisfied.

Next, we use this bound to control $\text{err}_{1}$ for iid estimators such as IPAM etc. The proof is in theorem 2 in appendix A). A simplified statement:

Theorem 2.

Consider a power-variance compatible triple $(\mathbb{X},d,\mathcal{D}_{\mathbb{X}})$. Consider $\hat{\hat{\nu}}$ with $S\perp T_{i}$ and iid $\{T_{i}\}_{i=1}^{m}$. Let the num samples $m,n$ and $p>0$ be such that the conditions of theorem 5 are met. Then:

[TABLE]

It is straightforward to show monotonic decrease of $\text{err}_{\text{ME}}$ for raising $m$ via Jensens or FTG theorem. Combined with the above result, this controls $\text{err}_{3}$.

3.3 Example and Experimental Verification

We give a concrete example to illustrate the above. And follow it up with experimental verification of the error bound.

Example 3** (Eccentricity of Gaussian).**

Let $(\mathbb{X},d,\mathcal{D}_{\mathbb{X}})=(\mathbb{R},{||}.{||}_{2}^{2},\mathcal{N}(0,1))$ and $K=1$. Then $\mathbb{X}^{K}=\mathbb{R}$. Let $S=\{\bar{X}_{i}\}_{i=1}^{m}$ and $T_{i}=\{X_{ij}\}_{j=1}^{n}$ be as per IPAM:

[TABLE]

Then

[TABLE]

Given the $\bar{X}_{i}$, the distances $\Delta(X_{ij},\bar{X}_{i})$ are distributed as

[TABLE]

Where $\chi^{2}(1,\lambda)$ denotes a non-central Chi-squared random variable with 1 degree of freedom and having mean, variance: $\mu=1+\lambda$, $\sigma^{2}=2(1+2\lambda)$. In our case We can rewrite the variance as $\sigma^{2}=4\mu-2$. This is a power-variance distribution with $\beta=1,\alpha=4,k=-2,\gamma=1$. Clearly all requirements on the parameters are satisfied, including $\alpha\gamma^{\beta}\geq-k$.

We need to uniformly upper bound kurtosis for our upper bounds to hold. The kurtosis $\kappa$ of non-central Chi-squared R.V is given by:

[TABLE]

This is strictly decreasing for all $\lambda>0$. And has a maximum at $\lambda=0$. This gives the upper bound:

[TABLE]

This distribution model is conceptually similar to the well known Tweedie family of distributions, but more suited to model families of distance distributions. In figure 1, we experimentally verify the error bound for the setting of example 3.

3.4 Related Work

MME like problems have been well studied under the rich theory of multi armed bandits. However our setting has significant differences from earlier settings. For instance, in stochastic linear bandits [1, 8], the decision space is $\mathbb{R}^{d}$ and rewards are linear. We don’t make any assumptions about the decision space.

A more direct comparison can be made to [2]. The decision space $D$ is arbitrary. But the rewards over $D$ are restricted to be a discrete set with a gap parameter $\Delta>0$. Analyzing K-Medoids true error requires consideration of a continuous rewards set with $\Delta\in[0,\infty)$. In exchange, our problem has more structure in terms of the relation between the rewards and noise (power variance), which they do not assume. Due to the differing nature of problems the error bounds are also different $\Theta(\log(n)/n\times m/\Delta^{4}_{\max})$ vs $\Theta(m^{1/3}/n^{2/3})$. A direct comparison is not possible, but we have taken their total regret, normalized by total amount of data ($mn$) as a rough analog and assumed all gaps $\Delta_{i\geq 2}=\Delta_{\max}$ (which decreases the bound). There are also links to the best arm identification problems found in [16, 4]. Again the discrete gap parameter $\Delta>0$ in both settings means that the upper bounds found in those papers go to infinity in our continuous setting $\Delta\in[0,\infty)$. We refer to our MME setting as familywise bandits. Since the worst case bound is over arms that are from a continuous family of distributions.

To our knowledge, we are the first to (i) define the true error $\text{err}_{3}$ (ii) decompose $\text{err}_{3}$ into canonical problems (iii) prove a convergence result for familywise bandits / $\text{err}_{1}$. Various K-medoid error analyses include [9, 37, 33, 4] (1-medoid), [43] (PAM) and [42] (PAMAE). However, error is calculated with respect to a sample medoid: EXHAUSTIVE $\hat{\hat{\nu}}$.

4 Faster, Large Scale $K$-Medoids

Lemma 7 (appendix D.1) suggests that for Gaussian and sub-Gaussian eccentricity distributions, $\text{err}_{1}$ is exceedingly fast decreasing in $n$ for a given tolerance $T$. Given that EVT error convergence rates are usually slow [27, 23], this suggests that more computational resources are used to control $\text{err}_{\text{ME}}$ and fewer to control $\text{err}_{\text{MME}}$. This suggests an optimal K-medoids regime of $m>>n$. This is inline with the 1-medoid algorithms RAND, TOPRANK, trimed and meddit And differs from the K-medoid algorithms PAM, CLARA, PAMAE ($m=\Theta(n)$) 444 To understand the statement $m=\Theta(n)$ for PAM, consider the earlier description of PAM. Corresponding to a single iteration of PAM steps 1 - 3, there are $kn$ single swaps. It is well known that there are a small number of PAM iterations, hence $m=\Theta(kn)=\Theta(n)$. Similarly for other algorithms

& CLARANS ($m<<n$). We propose a novel K-medoid algorithm MCPAM in the $m>>n$ regime.

MCPAM (Monte Carlo PAM) estimates the $n$ required for a given error tolerance using a Monte Carlo approach (listing 1). The core idea is to estimate $\mathrm{Ecc}$ of candidate medoids via sequential Monte Carlo sampling. We start with an initial medoid and $n=10^{3}$ and increment by factors of $10$ until a medoid swap having lower $\mathrm{Ecc}$ with high probability, is found. We chose $z_{\alpha}$ to give $1-\alpha$ coverage and calculate symmetric confidence intervals:

[TABLE]

Loop 5 of MCPAM is called MCPAM_INNER.

practical optimizations: $\tau,n_{\max}$ are practical controls for time vs accuracy. In the $m>>n$ regime, the confidence intervals of $\mathrm{\widehat{E}cc}(x^{\text{cur}},\{x_{i}\})$ are quite small compared to the error bars of $\mathrm{\widehat{E}cc}(x_{il},\{y_{i}\})$. As a practical optimization we use $\mathrm{\widehat{E}cc}(x^{\text{cur}},\{x_{i}\})$ instead of $\mathrm{\widehat{E}cc}_{\text{lo}}(x^{\text{cur}},\{x_{i}\})$, $\mathrm{\widehat{E}cc}_{\text{hi}}(x^{\text{cur}},\{x_{i}\})$. Practically, we find a strong correlation between $\mathrm{Ecc}_{\text{lo}}(\bar{x}_{\text{minlo}})$ & $\mathrm{Ecc}_{\text{hi}}(\bar{x}_{\text{minhi}})$. Hence we use $\bar{x}_{\text{minhi}}$ for both lines 9 and 11, and combine by exiting without swapping if $\mathrm{\widehat{E}cc}(\bar{x}_{\text{minhi}},\{y_{j}\})\pm z_{\alpha}\sigma_{\mathrm{\widehat{E}cc}(\bar{x}_{\text{minhi}},\{y_{j}\})}$ is within $\tau$ of $\mathrm{\widehat{E}cc}(x^{\text{cur}},\{x_{i}\})$. Else we either swap or continue the loop depending on whether $\mathrm{\widehat{E}cc}(x^{\text{cur}},\{x_{i}\})$ is higher or lower.

Single Medoid Variant: When $k=1$, the set of swaps is unchanged as $x^{\text{cur}}$ changes (line 4). This allows us to unify the outer and inner loops (lines 3 & 5). Hence we propose a simplified 1-medoid variant of MCPAM. We repeat at line 18 unconditionally, modify line 14 to be a continue of loop 5. and modify line 10 to terminate the program. Then the only way to exit loop 5 and the program is via condition 9.

Distributed MCPAM: The computationally heavy steps are 7,8 $\mathcal{O}(kmn)$ & 1,13 $\mathcal{O}(km)$. These are parallel operations on $\{x_{i}\}$ and each $x_{i}$ needs the full set $\{y_{j}\}$ for its calculations. Assume the $n$ chosen by MCPAM will be $n<<m$ (this is confirmed in theorem 6). Then, given $c$ worker nodes and 1 master node, $\{x_{i}\}$ is distributed in $c$ chunks to workers. Since $n<<m$, we simply sample on the master and broadcast a copy of $\{y_{j}\}$ to each worker. Our costs are given in theorem 6.

4.1 Theoretical Guarantees and Related Work

To simplify computational costs we assume $m,n>k$ in this subsection. We analyze MCPAM_INNER, with $\tau=0,n_{\max}=\infty$ (we will lose some accuracy in exchange for speed if $\tau>0,n_{\max}<\infty$). In the MCPAM setting, we assume that the minimum gap $\Delta$ between the eccentricities of all k-points constructable on $\{x_{i}\}$, is strictly positive and indpendent of $m$. This is similar to the assumption in [4]. Simplified versions of our convergence results are:

Theorem 6.

[MCPAM $K\geq 1$] When $\Delta>0$ and $\Delta=\Theta(1)$. MCPAM_INNER (loop 5) has expected runtime $\mathcal{O}(km)$. Upon exit from MCPAM_INNER, with high probability we have either: (i) a decrease in $\mathrm{Ecc}$ from $\mathrm{Ecc}(x^{\text{cur}})$ or (ii) the swap set constructed around $x^{\text{cur}}$ has no smaller $\mathrm{Ecc}$. MCPAM provides a true confidence interval on the eccentricity for the final medoid estimate. For the distributed version, computational cost is $\mathcal{O}(\frac{km}{c})$ and communication cost is $\Theta(1)$.

Theorem 7.

[MCPAM $K=1$] When $\Delta>0$ and $\Delta=\Theta(1)$. MCPAM has expected runtime $\mathcal{O}(m)$. It finds the sample medoid with high probability and provides a true confidence interval on the eccentricity of the same.

Proofs are in appendix E.1: theorems 7, 6. The MCPAM outer loop almost always runs only a few times (like other K-medoid methods), so $\mathcal{O}(km)$ is the practical runtime. Note, the confidence intervals are not for estimates of $\mathrm{Ecc}(\nu)$, but are true confidence intervals for $\mathrm{Ecc}(x^{\text{cur}})$. To our knowledge, MCPAM is the first fully general (i.e. semi-metric) K-medoids algorithm to have: (i) expected linear runtime in 1-medoid case (ii) computational cost of $\mathcal{O}(km)$ (detailed comparisons in ln 224 - 237) (iii) constant communication cost, providing massive scalability in distributed setting (iv) been scaled to 1 billion points. All this while closely matching PAM’s quality (fig 2).

We compare to the 1-medoid algorithms first. RAND [9] computes $\epsilon$-approximate sample 1-medoid w.h.p in $\mathcal{O}(\frac{m\log(m)}{\epsilon})$ time for all finite datasets. However, $\epsilon$ is measured relative to the network diameter $\Delta$, which can be fragile. The 1-medoid algorithms TOPRANK, trimed, meddit all have various distributional assumptions on the data. TOPRANK [37] finds the sample 1-medoid w.h.p in $\mathcal{O}(m^{\frac{5}{3}}\log^{\frac{4}{3}}m)$. trimed [33] finds the sample 1-medoid in $\mathcal{O}(m^{\frac{3}{2}}2^{\Theta(d)})$ for dimension $d$ while requiring the distances to satisfy the triangle inequality. This and the exponential dependence on $d$ significantly limit the applicability of trimed. meddit [4] finds the sample 1-medoid w.h.p in $\mathcal{O}(m\log(m))$ time, when $\Delta=\Theta(1)>0$. These are all worst case times.

The runtimes per iteration for PAM [26], CLARANS [34] are $\mathcal{O}(km^{2})$ (although swaps are subsampled in CLARANS, they are kept proportional to $km$). CLARA [25] is $\mathcal{O}(kn^{2})$, PAMAE [42] is $\mathcal{O}(m+rkn^{2})$ (distributed $\mathcal{O}(m+kn^{2})$), where $n$ is sample size typically set to a multiple of $k$ and $r$ is number of reruns. PAMAE has an accuracy guarantee, bounding the error in the estimate relative to the sample medoid. However, PAMAE is rather restrictive requiring the data come from a normed vector space.

4.2 Results and Comparisions

We perform scaling and quality comparisons on various datasets from literature [13, 30, 15, 45] and 2 newly created datasets (appendix E.2.1). We compare to PAM, as that is the ’gold standard’ in quality for K-medoids. Since PAM uses strictly more data than MCPAM (the $T_{i}$ are much smaller for MCPAM), PAM will have better quality. However, per our motivation, we expect MCPAM to have a small drop in quality with a massive speedup in runtime. This is verified in table 2 and 3. MCPAM was run with K++ initialization [3], 10 times on each dataset. PAM, MCPAM were run with the true number of clusters. We compare to DBSCAN as it is widely used for non-Euclidean distributed clustering [10, 18, 22]. To estimate epsilon parameter for DBSCAN, we (i) visually identify the knee of knn plot, (ii) do extensive grid search around that (iii) pick epsilon giving maximum ARI. The full details of our experimental setup are given in appendix E.2.2. We also compare MCPAM to PAM via clustering cost, again MCPAM is quite close to PAM in quality: figure 5 in appendix E.2.3.

The variance of eccentricity distributions tends to be correlated with the mean. Since $n$ is just required to control the variance of eccentricity distributions relative to their means, it is quite small for a variety of natural and artifical datasets, almost always $n=10^{3}$ in practice. This enables MCPAM to scale to very large dataset sizes. Figure 3 shows MCPAM scaling linearly as $m$ increases keeping underlying distribution same. Figure 4 shows distributed MCPAM results. We are, to the best of our knowledge, the first team to scale a semi-metric K-medoids algorithm to 1 billion points (with 6 attributes) (see BillionOne dataset appendix E.2.1). PAMAE [42] does run on a dataset of around 4 billion, but the data is restricted to Euclidean space, other distributed semi-metric K-medoids algorithms [24, 47, 21, 46, 31] are not run at this scale. To compare to other non-Euclidean algorithms, HPDBSCAN [18] demonstrates runs upto 82 million points with 4 attributes.

5 Conclusions

We have formalized the true K-medoid error $\text{err}_{3}$, and revealed a core dynamic inside this error by decomposing it into the ME and MME problems. We have a new ’familywise bandits’ convergence result to bound MME. Inspired by the tradeoffs in these bounds, we have proposed the first $K(>1)$ medoids algorithm in the $n<<m$ regime: MCPAM. MCPAM has good scaling and quality properties, and is the first K-medoids algorithm to provide true confidence interval on $\mathrm{Ecc}$ of estimated medoid. This work opens some interesting future directions: (i) Identifying conditions under which $\hat{\hat{\nu}}$ diverges. Based on our study of the core $z_{0}$ function in appendix C.1, we suspect there are diverging distributions when the power-variance parameter $\beta>2$. (ii) Our analysis combined with EVT provides a possible, natural approach to analyzing the number of K-medoid outer loop iterations.

Appendix A Medoid Error Analysis: $\text{err}_{1,2,3}$

In this section we analyze $\text{err}_{1},\text{err}_{2},\text{err}_{3}$ and provide an error decomposition and bounds.

Lemma 1** (Decomposition of $\text{err}_{3}$ into $\text{err}_{1}$, $\text{err}_{2}$).**

[TABLE]

Proof.

For a given $S$ and $\{T_{i}\}_{i=1}^{m}$:

[TABLE]

But $\operatorname*{\mathbb{E}}(\mathrm{Ecc}(\hat{\nu}))=(1+\text{err}_{2})\operatorname*{\mathbb{E}}(\mathrm{Ecc}(\nu))$. So:

[TABLE]

∎

Theorem 1** (Decomposition of $\text{err}_{3}$ into Canonical Problems).**

[TABLE]

Proof.

Consider a given $S$ and $\{T_{i}\}_{i=1}^{m}$. For each $\bar{X}_{i}\in S$ & $X_{ij}\in T_{i}$, $\Delta(X_{ij},\bar{X}_{i})$ is distributed with mean $\mathrm{Ecc}(\bar{X}_{i})$. From the definitions it is clear that $\hat{\hat{\nu}}$ is estimating the distribution with minimum mean from these distributions. Further, $\hat{\nu}$ is the distribution with minimum mean. Hence if we consider the minimum mean estimation problem with the randomly chosen means $\mathrm{Ecc}(\bar{X}_{i})$ $\bar{X}_{i}\in S$, we have:

[TABLE]

$S$ is a random sample of $\mathbb{X}^{K}$. Hence $\{\mathrm{Ecc}(\bar{X})|\bar{X}\in S\}$ is a random sample of $\{\mathrm{Ecc}(\bar{x})|\bar{x}\in\mathbb{X}^{K}\}$. The min $\hat{\nu}$ of the first set is estimating the min $\nu$ of the second set. Hence:

[TABLE]

The result is now immediate from lemma 1. ∎

Theorem 2** (Bound on $\text{err}_{1}$).**

Consider a power-variance compatible triple $(\mathbb{X},d,\mathcal{D}_{\mathbb{X}})$. Given samples $S,T_{1},\ldots,T_{m}$, let the following independence assumptions be satisfied:

• •

$S\perp T_{i}\quad\forall i\in[1,m]$**

• •

$X_{ij}\perp X_{i^{\prime}j^{\prime}}$* for $ij\neq i^{\prime}j^{\prime}$*

Let the num samples $m,n$ and $p>0$ be such that the conditions of theorem 5 are met. Then taking expectation over $S,T_{1},\ldots,T_{m}$, we have:

[TABLE]

Proof.

Expectations are over the joint distribution of $S,\{T_{i}\}_{i=1}^{m}$ unless noted otherwise. By independence assumption $S\perp\{T_{i}\}_{i=1}^{m}$ we have that the joint law of $S$ and $\{T_{i}\}_{i=1}^{m}$ is a product measure. By non-negativity of $\mathrm{Ecc}$ we can apply Fubini-Tonelli:

[TABLE]

Consider the inner integral in more detail:

[TABLE]

For a given $\bar{X}_{i}\in S$, $\Delta(X_{ij},\bar{X}_{i})$ is a function of R.V $X_{ij}\in T_{i}$ and is distributed with mean $\mathrm{Ecc}(\bar{X}_{i})$. From the definitions it is clear that $\hat{\hat{\nu}}$ is estimating the distribution with minimum mean from these distributions. Further, $\hat{\nu}$ is the distribution with minimum mean.

By assumption these distributions belong to a power-variance family. Further the assumption $X_{ij}\perp X_{i^{\prime}j^{\prime}}$ gives the required iid structure. Thus we have satisfied the conditions of theorem 5 (upper bound on $n$ for given error percentage $p$). So:

[TABLE]

Now:

[TABLE]

Where we have used equations 1 and 2. The result is now immediate.

∎

Appendix B MME Standard Estimator: Introduction

We define a standard estimator for the MME (Minimum Mean Estimation) problem defined in 7. We are interested in bounding its error. We do this in 3 parts. First we formulate the estimator and its expected error (appendix B). Next, we upper bound the error in the 2-D case, where we restrict to two distributions (appendix C). Finally, we generalize the bound to the M-D case, where we have $m$ distributions (appendix D).

B.1 MME Standard Estimator

Consider the setting of the MME problem (definition 7). We use the sample means $\hat{\mu}_{i}=\frac{1}{n}\sum_{j=1}^{n}E_{ij}$ to estimate $i_{\min}$.

Definition 11** (MME Standard Estimator).**

Let $\hat{i}_{\min}$ be the index of the minimum sample mean:

[TABLE]

Definition 12** (Expected Error of $\hat{i}_{\min}$).**

[TABLE]

Note $\epsilon=\operatorname*{\mathbb{E}}\text{err}_{\text{MME}}$. While dealing with the standard MME estimator (appendices B, C, D), we use a few notational conveniences:

• •

$\hat{i}_{\min}$ will refer specifically to the estimator in definition 11.

• •

$\operatorname*{\mathbb{E}}$ is implicitly over ${E_{ij}}$.

• •

Akin to order statistics we use the notation $\mu_{i:m}$ to denote the $i$th smallest $\mu_{i}$. There is the obvious mapping from the $i:m$ index to the $i$ index, and $\mu_{1:m}=\mu_{i_{\min}}$. Hereafter we use the $i:m$ indexing, for the most part.

We now derive a simple formula for the relative error.

Definition 13** (Relative Exceedances).**

The relative exceedance of the $i:m^{\mathrm{th}}$ mean:

[TABLE]

Definition 14** (Probability of Choosing $i:m^{\mathrm{th}}$ Mean).**

The probability of the $i:m^{\mathrm{th}}$ mean estimate undershooting all other mean estimates: 555For the case where the semi-metric is almost surely restricted to a countable set of values, ties are broken at random with equal probability for each possibility. :

[TABLE]

We now have the following obvious proposition.

Proposition 1** (Formula for Expected Relative Error).**

We can reorder the $\mu_{i}$ without changing $\epsilon$, and in particular:

[TABLE]

Also we have:

[TABLE]

In the next appendix we upper bound the relative error when $m=2$.

Appendix C The Two Mean Case

We want to upper bound $\epsilon$ for $m=2$. We do this in four subsections:

• •

First, we find a more convenient expression for relative error in the 2D case (subsection C.1)

• •

Then, we find an abstract upper bound for $\epsilon$ (subsection C.2). The abstract bound requires the existence of an upper bounding ’$G$’ function, satisfying certain properties.

• •

Subsequently, we show the existence of such a $G$ function for the Gaussian case (subsection C.3)

• •

Finally, we show the existence of a $G$ function for arbitrary distributions (subsection C.4), by extrapolating out of the Gaussian case via the Berry-Esseen theorem.

C.1 Convenient Expression for 2D Relative Error

In this subsection we derive a more convenient expression for relative error in the 2D case. Let the distributions $\mathcal{E}_{i}$ come from a power-variance distribution family $\mathcal{F}$ (see definition 9). The $\mathcal{E}$ are parametized by $\alpha,\beta,\gamma,k$, and hence so is $\epsilon$:

[TABLE]

For notational convenience we suppress the dependence on $\alpha,\beta,\gamma,k,n$, and write:

[TABLE]

Further, we assume the samples $E_{ij}$ (see definition 9) are independent i.e. $E_{i,j}\perp E_{i^{\prime},j^{\prime}}$ when $(i,j)\neq(i^{\prime},j^{\prime})$.

Definition 15** (2D Relative Error (Reparametrized)).**

Given the set $\{\mathcal{E}(\mu_{i:m})|i=1,2\}\subset\mathcal{F}$. Let $\delta_{2:2}$ be as in definition 13. For convenience set $\delta:=\delta_{2:2}$. So $\mu_{2:2}=(1+\delta)\mu_{1:2}$ and we have the obvious reparametrization of $\epsilon$:

[TABLE]

Definition 16** (Difference of Samples).**

Define $D_{j}$ as the difference of the j-th pair of samples from the distributions $\mathcal{E}_{i:m}$:

[TABLE]

Definition 17** (Difference of Sample Means).**

Define $\hat{d}$ as the difference of sample means:

[TABLE]

Definition 18**.**

Define the z-score of zero ($z_{0}$) for the random variable $\hat{d}$:

[TABLE]

Then using $\mu_{2:2}=(1+\delta)\mu_{1:2}$, we get:

[TABLE]

and:

[TABLE]

Proposition 2** (2D Relative Error Formula).**

For the 2D relative error from definition 15, we have:

[TABLE]

Where $\tilde{d}$ is the standardized (mean zero and unit variance) version of $\hat{d}$ and $F_{\tilde{d}}$ is its cumulative distribution function.

Proof.

The expected relative error is:

[TABLE]

But:

[TABLE]

$\hat{d}$ is standardized (mean zero and variance one) by the linear transformation:

[TABLE]

So

[TABLE]

∎

We recap the quantities used henceforth:

• •

$\delta=\frac{\mu_{2:2}-\mu_{1:2}}{\mu_{1:2}}$

• •

$z_{0}$ is the z-score of zero for R.V $\hat{d}$

• •

$F_{\tilde{d}}$ is the c.d.f of R.V $\tilde{d}$.

These quantities make $\epsilon$ tractable.

C.2 Upper Bounding $\epsilon$ for Two Means: Abstract Bound

In this section we derive an abstract upper bound on $\epsilon(\mu_{1:2},\delta)$. We do this in three stages:

• •

First we upper bound $z_{0}$ for a given $\delta$ over all $\mu_{1:2}$

• •

Second we upper bound $\epsilon$ for a given $\delta$ over all $\mu_{1:2}$

• •

Third we upper bound $\epsilon$ over $\mu_{1:2}$ and $\delta$

C.2.1 Upper Bound on $z_{0}$ for a given $\delta$ over all $\mu_{1:2}$

Define:

[TABLE]

We will upper bound this function in the regime $\beta\in[0,2]$. Consider:

[TABLE]

Where we have defined the constant:

[TABLE]

and in the last step we have used:

[TABLE]

Next, we upper bound the denominator term:

[TABLE]

Define

[TABLE]

When $\forall\delta\geq\delta_{c}$ and $\beta>0$, use $C_{1}\geq(1+\delta)^{\beta}$ and $\delta\geq 1$ to get:

[TABLE]

This also bounds $\mathrm{dr}$ for $\delta\leq\delta_{c}$, since $\mathrm{dr}$ is strictly increasing when $\beta>0$.

We get when $\beta>0$:

[TABLE]

This gives the upper bound on $z_{0}(\delta)$ for $\beta\in(0,2]$:

[TABLE]

When $\beta=0$:

[TABLE]

In both cases the upper bound may be interpreted as a piecewise function, initially linear and then a polynomial. In the $\beta=0$ case, the linear part occupies the whole of $\mathbb{R}^{\geq 0}$ since $\delta_{c}=\infty$ We now unify these two cases.

Definition 19** ($z_{0}^{\mathrm{UB}}$).**

[TABLE]

[TABLE]

[TABLE]

We now have:

[TABLE]

For a given value of $\beta$, $z_{0}^{\mathrm{UB}}$ is continuous on the closed subsets $0\leq\delta\leq\delta_{c}$ and $\delta_{c}\leq\delta<\infty$ (easier to see this by considering the $\beta>0$ case from equation 11 separately). Hence by pasting lemma, $z_{0}^{\mathrm{UB}}$ is continuous in $\delta$.

Next we prove the monotonicity of $z_{0}^{\mathrm{UB}}$. The $\beta=0$ case is trivial. Consider the case $\beta>0$. It is easy to see $z_{0}^{\mathrm{UB}}$ is strictly decreasing on $0\leq\delta\leq\delta_{c}$ (linear part). The polynomial $C\delta^{1-\frac{\beta}{2}}$ is strictly decreasing $\forall\delta>0$, when $\beta<2$. Hence $z_{0}^{\mathrm{UB}}$ is strictly decreasing on $\delta_{c}\leq\delta$. We collect all the above in a proposition.

Proposition 3** (Simple Upper Bound on $z_{0}(\delta)$).**

The function $z_{0}^{\mathrm{UB}}(\delta)$ (definition 19) is continuous, strictly decreasing and invertible on $[0,-\infty)$. Furthermore:

[TABLE]

C.2.2 Upper Bound on $\epsilon$ for a given $\delta$ over all $\mu_{1:2}$

Define:

[TABLE]

Then we have:

Proposition 4** (Upper Bound on $\epsilon(\delta)$).**

If $F_{\tilde{d}}$ is continuous and $\beta\in[0,2]$, then:

[TABLE]

Proof.

Now by proposition 1 (formula for relative error)

[TABLE]

We have used proposition 3 and non-decreasing property of a cdf for the last step. ∎

C.2.3 Maximizing Over $\delta$

We seek an upper bound on $\epsilon(\delta)$ of the following form (essentially a tail bound). For some $T,\delta^{\mathrm{th}}>0$:

[TABLE]

Definition 20** ($\epsilon^{\mathrm{UB}}$).**

Define

[TABLE]

Where $G$ is a function satisfying:

[TABLE]

Then we have the following bound.

Lemma 2** (Abstract Upper Bound on $\epsilon(\delta)$ for Two Means).**

Consider a power-variance distribution family $\mathcal{F}$, with $\beta\in[0,2]$. Let the corresponding $F_{\tilde{d}}$ be continuous. Then, if there exists $G(x)$ such that:

[TABLE]

Then $\epsilon^{\mathrm{UB}}$ (definition 20) satisfies:

[TABLE]

Proof.

We denote the upper bound from proposition 4 by $\epsilon^{\text{UBInt}}$.

[TABLE]

By proposition 3 (properties of $z_{0}^{\mathrm{UB}}$), the inverse of $z_{0}^{\mathrm{UB}}$ exists and we can rewrite $\epsilon^{\text{UBInt}}$ as a function of $z_{0}^{\mathrm{UB}}$. We have $\forall\delta>0$:

[TABLE]

Where we have abused notation using $\epsilon^{\text{UBInt}}(z_{0}^{\mathrm{UB}}(\delta))$ to denote the dependence of $\epsilon^{\text{UBInt}}$ on $\delta$ solely via $z_{0}^{\mathrm{UB}}$.

Then we have a new upper bound

[TABLE]

∎

This is the chief result of this subsection. This is one of our core results, because it is applicable with great generality and gives us an easy way to bound errors for any distribution family $\mathcal{F}$, provided we can find a suitable $G$.

C.3 Upper Bounding $\epsilon$ for Two Means: Gaussian Case

In this subsection we assume that the distribution family ($\mathcal{E}$) is Gaussian ($\mathcal{N}$). As usual the family is further specified by the parameters $\alpha,\beta,\gamma,k,n$. We start by finding an upper bound for $F_{\tilde{d}}=\Phi$. Here $\Phi$ is the cdf of a standard normal. We use an ubiquitous tail bound (proof omitted).

Proposition 5** (A Gaussian Tail Bound).**

This is reproduced from [11]. Let $a>0$, then:

[TABLE]

Then we set

[TABLE]

Observe 666 When $\beta\in[0,2)$, $z_{0}^{\mathrm{UB}}$ is unbounded and so $z_{0}^{\mathrm{UB}}([0,\infty))=[0,-\infty)$. Further $z_{0}^{\mathrm{UB}}(0)=0$ always. that $z_{0}^{\mathrm{UB2}-1}(x)>0$ $\forall x\in(0,-\infty)$ . So by using the tail bound in proposition 5:

[TABLE]

Thus satisfying equation 17 and so $\forall\delta>0$:

[TABLE]

Substituting for $z_{0}^{\mathrm{UB}}$:

[TABLE]

Given a function $f$, for the log (or any strictly increasing function) we have:

[TABLE]

Consider the function:

[TABLE]

Then

[TABLE]

Hence $\frac{\partial f}{\partial x}<0$ iff:

[TABLE]

Where on the last line we used that $g(y)=y^{\frac{1}{2-\beta}}$ is strictly increasing $\forall\beta\in[0,2)$ and $y\geq 0$ We want $f$ to be increasing for all $x\geq\delta_{c}$. So set $x=\delta_{c}$

[TABLE]

Now $\epsilon^{\mathrm{UB3}}_{\mathcal{N}}(\delta)$ is strictly decreasing when the above equation is satisfied. Hereafter we assume this. $\mathcal{R}(f)$ denotes the range of function $f$. Then with $\delta>0$:

[TABLE]

Note: In the $\beta>0$ case:

[TABLE]

For the last step we have used finite numberr of applications of L’Hopital. Hence $T\in\mathcal{R}(\epsilon^{\mathrm{UB3}}_{\mathcal{N}})$ is equivalent to:

[TABLE]

Given a target $T>0$ to bound relative error, we search for a corresponding $\delta^{\mathrm{th}}_{\mathcal{N}}$. If $T\geq\frac{1}{\sqrt{2\pi}C_{2}n^{\frac{1}{2}}}$ ($T$ larger than max-range), then $\epsilon^{\mathrm{UB3}}_{\mathcal{N}}(\delta)<T\quad\forall\delta>0$. Hence $\delta^{\mathrm{th}}_{\mathcal{N}}=0$ in this case.

Next assume $T<\frac{1}{\sqrt{2\pi}C_{2}n^{\frac{1}{2}}}$. Set:

[TABLE]

Now $T\in\mathcal{R}(f_{1})=\mathcal{R}(\epsilon^{\mathrm{UB3}}_{\mathcal{N}})$ and

[TABLE]

When $f_{1}^{-1}(T)\leq\delta_{c}$ i.e. when $n\geq\frac{1}{2\pi C_{2}^{2}T^{2}}\exp\left[-C_{2}^{2}n\delta_{c}^{2}\right]$, we have

[TABLE]

Hence $\delta^{\mathrm{th}}_{\mathcal{N}}=\left(-\frac{2}{C_{2}^{2}n}\log(\sqrt{2\pi n}C_{2}T)\right)^{\frac{1}{2}}$ in this case.

By lemma 2 (abstract upper bound on $\epsilon(\delta)$), we get

[TABLE]

We collect all these in a lemma.

Lemma 3** (Upper Bound on $\epsilon$ for Two Means: Gaussian Case).**

When $\mathcal{E}$ is the Gaussian family and given values for $\alpha,\beta,\gamma,k,n$. Define:

[TABLE]

This is an upper bound:

[TABLE]

Additionally, given a $T>0$ if $n$ satisfies:

[TABLE]

Then $\epsilon^{\mathrm{UB3}}_{\mathcal{N}}$ is strictly decreasing. And we can define a $\delta^{\mathrm{th}}_{\mathcal{N}}$

[TABLE]

Such that:

[TABLE]

Typically we will use this result as follows. Given a target $T$. We will first find the minimum $n$ that satisfies equations 3, 3. Then we can choose $n\geq n_{\min}$ to make $\delta^{\mathrm{th}}_{\mathcal{N}}$ as small as we want.

C.4 Upper Bounding $\epsilon$ for Two Means: General Case

In this subsection we work with an arbitrary distribution family ($\mathcal{E}$). As usual the family is specified by the parameters $\alpha,\beta,\gamma,k,n$. And we want to find a worst case upper bound on $\epsilon$. We split this into two stages. In subsection C.4.1 we find a familywise tail error bound between $F_{\tilde{d}}$ and $\Phi$. In subsection C.4.2 we use this to bound $\epsilon_{\mathrm{gen}}$.

C.4.1 Familywise Tail Bound on $F_{\tilde{d}}$

By the CLT, the Gaussian cdf $\Phi$ is an attractor for the cdf of $\hat{d}_{S}$. Hence we could attempt to generalize the bound in lemma 3 (Gaussian upper bound on $\epsilon$) to non-Gaussian distributions by setting $F_{\tilde{d}}=\Phi$, when ‘$n$ is large enough’. However, there are significant issues with a naive application of CLT in our context.

Varying Distributions

Firstly, we are applying the CLT over a family of distributions. That is, the random variable $\hat{d}$ is the sample mean of one out of a possible family of distributions and not one fixed distribution. Further the family is infinite (indexed by real valued parameters $\mu_{1:2},\mu_{2:2}$). It is easy to construct families such that the CLT requires infinite sample size $n$ to converge to a given tolerance for all members of the family.

Central Convergence vs Tail

The CLT is known to converge fast near the mean. And primarily that is how the CLT is used (to construct confidence intervals around the mean); However we are interested in the tail of the cdf as well. Convergence may require vastly more samples than usual.

To deal with the above challenges we develop a familywise version of the non-uniform Berry-Esseen theorem for the difference of sample averages. Our starting point is the following Berry-Esseen theorem:

Theorem 3** (Non-Uniform Berry-Esseen [38]).**

Let $X_{i}$ be independent random variables such that:

[TABLE]

Then set $W_{n}:=\sum_{i=1}^{n}X_{i}$ so $\mathrm{Var}(W_{n})=1$. And let $F_{W_{n}}$ denote the cdf of $W_{n}$. Then:

[TABLE]

Where $C_{4}<32$.

Then we have a simple corollary.

Proposition 6** (Non-Uniform Berry Esseen for Sample Average).**

Let $Y_{i}$ be $n$ i.i.d random variables such that:

[TABLE]

Let $A_{n}:=\frac{1}{n}\sum_{i}Y_{i}$ be the sample average. Let $W_{n}:=\frac{A_{n}-\mu}{\sigma/\sqrt{n}}$ be the standardized version of $A_{n}$. And let $F_{W_{n}}$ denote the cdf of $W_{n}$. Then:

[TABLE]

Proof.

With the above setup define:

[TABLE]

So

[TABLE]

Thus $X_{i}$ are i.i.d and

[TABLE]

Further

[TABLE]

The conditions of theorem 3 are satisfied and we have:

[TABLE]

∎

For distance distributions it can be hard to calculate $\mathbb{E}|Y_{i}-\mu|^{3}$. We use the common trick of bounding by a higher moment via Jensens. But first a technical statement.

Proposition 7**.**

Let $Y_{0}$ be a random variable and $g_{1}$, $g_{2}$ be functions. Define:

[TABLE]

Then:

[TABLE]

Proof.

Clearly, if $g$ is a function and $Y$ is a random variable, we have:

[TABLE]

If $h$ is a function and define r.v $Z:=g(Y)$, we have:

[TABLE]

By using the above two statements:

[TABLE]

And:

[TABLE]

∎

Proposition 8** (Upper Bound of Third Absolute Moment).**

[TABLE]

Proof.

Recall Jensen’s. When $g$ is concave and finite:

[TABLE]

Let $g(x)=x^{\frac{3}{4}}$. Then $g(x)$ is concave and we have:

[TABLE]

Set $X:=Y^{4}$

[TABLE]

[TABLE]

∎

Proposition 9** (Upper Bound of Third Absolute Moment of a Standardized Difference).**

Let $\kappa_{Y}:=\mathbb{E}\left(\frac{Y-\mu_{Y}}{\sigma_{Y}}\right)^{4}$ denote the kurtosis of a random variable $Y$. Let $Y_{1},Y_{2}$ be independent random variables with $\sigma_{Y_{i}},\kappa_{Y_{i}}<\infty$. and set $Y:=Y_{2}-Y_{1}$. Then:

[TABLE]

Proof.

Define the standardized version of $Y$ as $X:=\frac{Y-\mu_{Y}}{\sigma_{Y}}$. Hence:

[TABLE]

Then by proposition 8 we get:

[TABLE]

Next we derive a formula for $\kappa_{Y}$. The standard formula for kurtosis of sum of 2 random variables is: $X_{1},X_{2}$:

[TABLE]

Where $\kappa(,,,)$ denotes the cokurtosis function. For independent random variables:

[TABLE]

So

[TABLE]

Now set $X_{1}=Y_{2}$ and $X_{2}=-Y_{1}$ to get:

[TABLE]

But $\sigma_{-Y_{1}}=\sigma_{Y_{1}}$, $\kappa_{-Y_{1}}=\kappa_{Y_{1}}$ and $\sigma_{Y_{2}-Y_{1}}=\sqrt{\sigma_{Y_{2}}^{2}+\sigma_{Y_{1}}^{2}}$ so:

[TABLE]

Combining with equation 30 we get:

[TABLE]

∎

Proposition 10** (Familywise Upper Bound on Third Absolute Moment of a Standardized Difference).**

Given a pair $\mu_{1:2},\mu_{2:2}\in[0,\infty)$, Let D be as defined in 16. Let the distribution family $\mathcal{F}$ be such that the kurtosis of any member of the family is upper-bounded (possibly tightly) by a constant $\kappa^{\mathrm{UB}}$. Then:

[TABLE]

Proof.

Define a new familywise upper bound on the kurtosis

[TABLE]

Given an arbitrary pair $\mu_{1:2},\mu_{2:2}$,

Consider an arbitrary $\mu_{1:2},\mu_{2:2}$ pair. Consider $D$ as in definition 16. $D=E_{2:2}-E_{1:2}$ and $E_{2:2}\perp E_{1:2}$. Let $\sigma_{i:2}^{2}=\mathrm{Var}{E_{i:2}}$ and $\kappa_{i:2}=\kappa_{E_{i:2}}$. Since $\sigma_{i:2}<\infty$ (finite function of mean) & $\kappa_{i:2}\leq\kappa^{\mathrm{UB2}}<\infty$, we can apply proposition 9 to get:

[TABLE]

[TABLE]

Clearly this expression is strictly increasing in $\kappa_{i:2}$, keeping everything else fixed. So we have:

[TABLE]

∎

Now we bound the difference between the functions $F_{\tilde{d}}$ and $\Phi$.

Proposition 11** (Non Uniform Berry Esseen for $F_{\tilde{d}}$).**

: Consider $F_{\tilde{d}}$ as defined in proposition 2. Given a distribution $\mathcal{E}(\mu)\in\mathcal{F}$, let $\kappa(\mu)$ denote it’s kurtosis as a function of $\mu$. Further, let:

[TABLE]

And define:

[TABLE]

Then for any $\mu_{1},\mu_{2}\in[\gamma,\infty)$:

[TABLE]

Proof.

Given a pair $\mu_{1},\mu_{2}\in[\gamma,\infty)$, by proposition 10:

[TABLE]

And since $\sigma_{D}<\infty$

[TABLE]

Additionally $D_{i}$ are iid. We can now apply proposition 6 by setting $Y_{i}:=D_{i}$. We get:

[TABLE]

and finally:

[TABLE]

∎

C.4.2 Bounding $\epsilon_{\mathrm{gen}}$

Herafter we assume $\mathcal{F}$ is such that the kurtosis are uniformly bounded:

[TABLE]

Then by the above proposition $\forall x<0$:

[TABLE]

Now $z_{0}^{\mathrm{UB2}-1}(x)>0\ \ \forall x\in(0,-\infty)$. So:

[TABLE]

And we have the required $G$ satisfying equation 17:

[TABLE]

We split this into two terms:

[TABLE]

Where:

[TABLE]

Correspondingly define the relative error components:

[TABLE]

So:

[TABLE]

Easy to see that $\epsilon^{\mathrm{UB3}}_{\mathrm{BE}}$ is strictly decreasing when $\beta<\frac{4}{3}=1.33\ldots$

Now we define:

[TABLE]

$\epsilon^{\mathrm{UB3}}_{\mathrm{gen}}(\delta)$ is strictly decreasing when equation 3 is satisfied and $\beta>\frac{4}{3}$. Henceforth we will assume these conditions are satisfied.

We now seek a formula to invert $\epsilon^{\mathrm{UB3}}_{\mathrm{BE}}$. Let

[TABLE]

Then $\mathcal{R}(f_{1})=(\infty,0)=\mathcal{R}(\epsilon^{\mathrm{UB3}}_{\mathrm{BE}})$. So, given $T>0$ we have the inverse:

[TABLE]

We have:

[TABLE]

The condition $f_{1}^{-1}(T)\leq\delta_{c}$ is equivalent to:

[TABLE]

Now we derive a formula for the inverse of $\epsilon^{\mathrm{UB3}}_{\mathrm{gen}}$. Directly inverting $\epsilon^{\mathrm{UB3}}_{\mathrm{gen}}$ (in the $\delta\leq\delta_{c}$ region) will give a transcendental equation. We get around this with a simple approach. We invert both the component functions ($\epsilon^{\mathrm{UB3}}_{\mathcal{N}}$ and $\epsilon^{\mathrm{UB3}}_{\mathrm{BE}}$) for a target of $\frac{T}{2}$, and take the max of the resulting $\delta$. Since these are strictly decreasing functions, their sum will be $<T$ at that $\delta$. We now work out the details.

Given a desired target $T>0$, let us assume conditions are satisfied ensuring monotonicity and first piecewise invertibility777By first piecewise invertibility we mean that the function inverse is the inverse of the first piecewise segment of $\epsilon^{\mathrm{UB3}}_{\mathcal{N}}$ and $\epsilon^{\mathrm{UB3}}_{\mathrm{BE}}$. Note that the conditions for first piecewise invertibility of $\epsilon^{\mathrm{UB3}}_{\mathcal{N}}$ and $\epsilon^{\mathrm{UB3}}_{\mathrm{BE}}$ are to be applied for a target of $\frac{T}{2}$ and not $T$. Then consider the case:

[TABLE]

Set:

[TABLE]

Next consider the case: $\frac{T}{2}\geq\frac{1}{\sqrt{2\pi}C_{2}n^{\frac{1}{2}}}$ By lemma 3 (Gaussian upper bound on rel error two means):

[TABLE]

So set:

[TABLE]

Then in both cases:

[TABLE]

Hence define $\delta^{\mathrm{th}}_{\mathrm{gen}}(T):=$

[TABLE]

And get:

[TABLE]

Then by lemma 2 (abstract upper bound on relative error):

[TABLE]

We collect all the above in a lemma.

Lemma 4** (Upper Bound on $\epsilon$ for Two Means: General Case).**

Let $\mathcal{E}$ be some distribution family having parameters $\alpha,\beta,\gamma,k,n$. For all distributions $\mathcal{E}_{i:m}$ in this family, let the cdf be continuous and the kurtosis be uniformly bounded:

[TABLE]

Let $\epsilon^{\mathrm{UB3}}_{\mathrm{BE}},\epsilon^{\mathrm{UB3}}_{\mathrm{gen}}$ be defined as:

[TABLE]

This is an upper bound:

[TABLE]

Given a $T$, let $n,T$ satisfy the following conditions:

[TABLE]

and define $\delta^{\mathrm{th}}_{\mathrm{gen}}(T):=$

[TABLE]

Then:

[TABLE]

This $\delta^{\mathrm{th}}_{\mathrm{gen}}$ formula tells us that the $\delta$ threshold required to ensure $\epsilon(\delta)$ less than a specified tolerance is a relatively rapidly decreasing function of the threshold $T$ (either behaving like $\sqrt{-C\log(CT)}$ or $\frac{C}{T^{\frac{1}{2}}}$). For a fixed threshold $T$, the $\delta^{\mathrm{th}}_{\mathrm{gen}}$ is a rapidly decreasing function of $n$ as well (either behaving like $\sqrt{-C\frac{\log(Cn)}{n}}$ or $\frac{C}{n}$). Combined, these indicate small requirements on $n$ to hit a target relative error.

Our usage of the error bound will be to fix an acceptable error threshold $T$ and to fix an acceptable threshold on $\delta^{\mathrm{th}}_{\mathrm{gen}}=\epsilon^{\mathrm{UB3-1}}_{\mathrm{gen}}(T)$. Then we find $n$ as small as possible that guarantees both simultaneously. The bound on $\delta^{\mathrm{th}}_{\mathrm{gen}}$ is required to generalize from the 2-dimensional case (min of two sample means) to the m-dimensional case (min of m sample means).

Appendix D M-Mean Case

In this section we derive a series of results upper bounding $\epsilon$, when we have more than $2$ means. In the first subsection D.1 (Generalization of $2$ Mean Upper Bounds to $m$ means) we prove an important and central result (lemma 5) that generalizes a $2$ mean bound on $\epsilon$ to a $m\geq 2$ mean bound, for all distributions with full generality. In the second subsection D.2 (M-mean Case: Main Results), we state and prove the fundamental result for the minimum mean estimation problem. This is theorem 4

D.1 Generalization of $2$ Mean Upper Bounds to $m$ means

The central result is lemma 5 which provides a generalization mechanism for all distributions provided the $2$ mean case is bounded. Subsequent results generalize the various two mean results of section C.

We start by reducing the $m$ mean $\mathbb{P}_{i:m}$ to its two mean counterpart.

Proposition 12** ($\mathbb{P}_{i:m}$ Upper Bound: Reduction from $m$-Means to $2$-Means).**

:

[TABLE]

Proof.

Using the definition of $\mathbb{P}_{i:m}$ from definition 14:

[TABLE]

∎

Lemma 5** ($\epsilon$ Upper Bound: Reduction from $m$ Means to $2$ Means).**

Given a distribution family $\mathcal{F}$. Suppose the relative error $\epsilon$ in the $2$ mean case is bounded as:

[TABLE]

For some $\delta^{\mathrm{th}},T$. Then define a function:

[TABLE]

This is an upper bound on the relative error in the $m$ mean case:

[TABLE]

Proof.

Let $\mu^{\mathrm{th}}:=\mu_{1:m}(1+\delta^{\mathrm{th}})$ . Now given $\mu_{i:m}$, let $l$ be the index such that

[TABLE]

Then we can split the relative error (using proposition 1) into ‘head’ and ‘tail’ components.

[TABLE]

We can bound the head component in the following manner. Observe $\mu_{i:m}\leq\mu_{1:m}(1+\delta^{\mathrm{th}})\Rightarrow\delta_{i:m}\leq\delta^{\mathrm{th}}$. Hence:

[TABLE]

Where we have used $\sum_{i=1}^{m}\mathbb{P}_{i:m}(\mu_{1:m},\ldots,\mu_{m:m})=1$. Next we bound the tail component.

[TABLE]

Combining:

[TABLE]

∎

Lemma 6** (Abstract Upper Bound on $\epsilon$).**

Given a distribution family $\mathcal{F}$, specified by parameters $\alpha,\beta,\gamma,k,n$. Let the correspoding $F_{\tilde{d}}$ be continuous. Given $m$ distributions from $\mathcal{F}$ specified by:

[TABLE]

Let the conditions of lemma 2 (abstract upper bound on $\epsilon$ for two means) be satisfied. We define:

[TABLE]

for $\delta^{\mathrm{th}},T$ as defined in lemma 2. This is an upper bound on $\epsilon$.

[TABLE]

Proof.

Since the conditions of lemma 2 (abstract upper bound on $\epsilon$ two means) are satisfied we have:

[TABLE]

For $\delta^{\mathrm{th}},T$ as defined in the same lemma. Then by lemma 5 ($\epsilon$ upper bound reduction from $m$ means to $2$ means), we have the upper bound:

[TABLE]

∎

Lemma 7**.**

*(Upper Bound on $\epsilon$: Gaussian Case)

Let $\mathcal{F}$ be the Gaussian family with given values for $\alpha,\beta,\gamma,k,n$. Further, given $m$ distributions from $\mathcal{F}$, specified by*

[TABLE]

Then given a $T>0$, suppose $n$ satisfies equations 3 and 3:

[TABLE]

And $\delta^{\mathrm{th}}_{\mathcal{N}}$ is defined as in equation 27:

[TABLE]

Then define

[TABLE]

Then

[TABLE]

Proof.

Since the conditions of lemma 3 (upper bound on $\epsilon$ with two means: Gaussian case) are satisfied, we have:

[TABLE]

For $T>0$ and $\delta^{\mathrm{th}}_{\mathcal{N}}$ as defined in above lemma. Hence by lemma 5 ($\epsilon$ upper bound reduction from $m$ means to $2$ means):

[TABLE]

∎

Lemma 8**.**

*(Upper Bound on $\epsilon$: General Case)

Let $\mathcal{F}$ be some distribution family, having parameters $\alpha,\beta,\gamma,k,n$. For all distributions $\epsilon(\mu)$ in this family, let the cdf be continuous and the kurtosis be uniformly bounded:*

[TABLE]

Consider $m$ distributions from $\mathcal{F}$, specified by:

[TABLE]

$\mathcal{R}(f)$* denotes the range of function $f$. Given a T, let $n,T$ satisfy the following conditions*

[TABLE]

and define $\delta^{\mathrm{th}}_{\mathrm{gen}}(T;\alpha,\beta,\gamma,k,n):=$

[TABLE]

Then if we define:

[TABLE]

We have:

[TABLE]

Proof.

Since the conditions of lemma 4 (upper bound on $\epsilon$ with two means: General case) are satisfied, we have:

[TABLE]

Hence by lemma 5 ($\epsilon$ upper bound reduction from $m$ means to $2$ means):

[TABLE]

∎

D.2 M-Mean Case: Main Results

In this subsection, we derive a fundamental bound (theorem 4) on the relative error $\epsilon$ for the minimum mean estimation problem. This is a simple closed form bound on the error convergence rate:

[TABLE]

This result is the analog of the standard square root rate of convergence of Monte Carlo algorithms. In theorem 5, we derive conditions on $m,n$ to hit a given tolerance.

Our approach is to simplify lemma 8 (upper bound on $\epsilon$ general case) by expressing the free variable $T$ in terms of $m,n$. We choose a value $T=T^{*}$ to decrease the upper bound.

Theorem 4** (Minimum Mean Estimation: Error Convergence Rate).**

Let $\mathcal{F}$ be some distribution family, having parameters $\alpha,\beta,\gamma,k,n$. For all distributions $\epsilon(\mu)$ in this family, let the cdf be continuous and the kurtosis be uniformly bounded:

[TABLE]

Given $\mu_{i:m}\ \ i=1,\ldots,m$, let the following conditions be satisfied:

[TABLE]

Further if either, the below pair of conditions are both satsified:

[TABLE]

or the below single condition is satisfied:

[TABLE]

then we have a simplified upper bound on the error:

[TABLE]

Where

[TABLE]

Proof.

We can vary the free variable $T$ to decrease the upper bound $\epsilon^{\mathrm{UB4}}_{\mathrm{gen}}(T,m)$. We first provide motivation for our choice of $T$. Let us assume the upper bound is of the following form:

[TABLE]

$f_{2}(T)$ is the second term (BE term) in the $\max$ term in $\delta^{\mathrm{th}}_{\mathrm{gen}}$. Roughly, we are assuming that $\delta^{\mathrm{th}}_{\mathrm{gen}}$ simplifies to $\epsilon^{\mathrm{UB3-1}}_{\mathrm{BE}}$ under widely valid conditions. Then we want to find the minimizer $T^{*}$ of:

[TABLE]

[TABLE]

equating to zero gives

[TABLE]

Further

[TABLE]

Hence $T^{*}$ is the global minimum. The corresponding upper bound:

[TABLE]

Where $C_{7}:=\frac{C_{6}^{\frac{1}{3}}}{C_{2}}\frac{3}{2^{\frac{1}{3}}}$. Now we ask given a $m,n$ when such an upper bound might hold. We answer this by plugging in the value for $T^{*}$ into the conditions of lemma 8 (upper bound on $\epsilon$: general case). Consider the condition:

[TABLE]

Plugging in $T^{*}$:

[TABLE]

Next consider the condition:

[TABLE]

Plugging in $T^{*}$:

[TABLE]

Since conditions 42, 43 are satisfied by assumption and further:

[TABLE]

are also satisfied by assumption, the conditions for lemma 8 (upper bound on $\epsilon$: general case) hold at $T=T^{*}$. Thus we have the upper bound on $\epsilon_{\mathrm{gen}}$

[TABLE]

Where $\delta^{\mathrm{th}}_{\mathrm{gen}}(T)$ is as defined in equation 37. It remains to see when $\delta^{\mathrm{th}}_{\mathrm{gen}}(T^{*})=f_{2}(T^{*})$. We rewrite $\delta^{\mathrm{th}}_{\mathrm{gen}}$. Let

[TABLE]

Where

[TABLE]

Then:

[TABLE]

Now, what are the conditions under which $\delta^{\mathrm{th}}_{\mathrm{gen}}(T^{*})=f_{2}(T^{*})$. We need conditions such that

[TABLE]

This will happen when $T^{*}\geq T_{c}$, i.e. when:

[TABLE]

Next, when $T^{*}<T_{c}$, $\delta^{\mathrm{th}}_{\mathrm{gen}}(T^{*})$ will equal $f_{2}(T^{*})$ if $f_{1}(T^{*})\leq f_{2}(T^{*})$. That is:

[TABLE]

So $\delta^{\mathrm{th}}_{\mathrm{gen}}(T^{*})=f_{2}(T^{*})$ when the following conditions are both met:

[TABLE]

So $\delta^{\mathrm{th}}_{\mathrm{gen}}(T^{*})=f_{2}(T^{*})$ if either condition 44 or the above pair are satisfied. We note that the above pair is easier to satisfy than condition 44. Then:

[TABLE]

and define:

[TABLE]

Completing the proof. ∎

We want to understand the feasible region of $n$ in the above theorem. In practice $n$ will mostly satisfy the pair of conditions:

[TABLE]

and not:

[TABLE]

Hence we will restrict our study to the feasible region when $n$ satisfies the former pair in conjunction with the other conditions on $n$ from theorem 4. We start with an useful definition

Definition 21** (Extended Inverse).**

Given $f:\mathcal{D}\rightarrow\mathcal{R}$ and a target $y^{{}^{\prime}}$. If $|f^{-1}(y^{{}^{\prime}})|\leq 1$, then define the extended inverse $x^{{}^{\prime}}$ of $y^{{}^{\prime}}$ as:

[TABLE]

Note that the extended inverse coincides with the inverse when $y^{{}^{\prime}}\in\mathcal{R}$. Hence, hereafter we use the same notation $f^{-1}(y^{{}^{\prime}})$ for both.

Proposition 13** (Feasible Region for Theorem 4).**

Consider the conditions on $n$ in theorem 4 (minimum mean estimation: error convergence rate). When:

[TABLE]

a feasible region for these conditions is:

[TABLE]

Where $g_{i}$ are the conditions of theorem 4 expressed in functional form, equations 45, 47, 48, 49, 50 respectively. And the inverses are the extended inverses as per definition 21.

Proof.

Consider the conditions on $n$ in theorem 4, specifically where ’conditions on first piecewise component of $\delta^{\mathrm{th}}_{\mathrm{gen}}$’ are being satisfied in conjunction with the rest. These $5$ conditions may be split into four lower bounds and one upper bound. We rewrite these conditions in functional form. The domain for all these functions will be $[1,\infty)$. We start by rewriting the lower bounds. Define:

[TABLE]

Then the lower bound conditions are:

[TABLE]

If for all $n\geq 1$ such that $g_{i}(n)\geq 0$:

[TABLE]

then $g_{i}$ is strictly increasing when it is non-negative valued. Hence the $g_{i}$ have an unique pre-image for [math] (if it exists):

[TABLE]

Again by the above monotonicity property and using definition 21 (extended inverse), we can write the feasible regions as

[TABLE]

We now establish the monotonicity of the $g_{i}$, starting with $g_{2}$.

[TABLE]

If $n\geq 1$ is such that $g_{2}(n)\geq 0$, then:

[TABLE]

But:

[TABLE]

and so:

[TABLE]

then:

[TABLE]

The other three lower bound $g_{i}$ have derivative $1$ everywhere. Next consider the upper bound condition, we rewrite this initially as:

[TABLE]

We will now rewrite this in a more tractable form and establish the same monotonicity property for it as well. Consider the upper bound condition:

[TABLE]

We have $\log\left(\frac{\pi^{\frac{1}{2}}C_{6}^{\frac{1}{3}}}{2^{\frac{5}{6}}m^{\frac{2}{3}}n^{\frac{1}{6}}}\right)>0$ when: $\frac{\pi^{\frac{1}{2}}C_{6}^{\frac{1}{3}}}{2^{\frac{5}{6}}m^{\frac{2}{3}}n^{\frac{1}{6}}}<1$. But:

[TABLE]

This is satisfied when $g_{4}(n)>0$ i.e. when $n>g_{4}^{-1}(0)$. Then for such a $n$ we can rewrite the upper bound condition as:

[TABLE]

then set:

[TABLE]

then the upper bound condition is the same as:

[TABLE]

Now we want to show the strict increase of $g_{5}$. It suffices to show the strict increase of

[TABLE]

But:

[TABLE]

And so it suffices to show the strict increase of $\log(h)$. This in turn is equivalent to:

[TABLE]

When $n>0$, suffices to have:

[TABLE]

Clearly this is satisfied for all $n\geq 1$. Now we have:

[TABLE]

So:

[TABLE]

But:

[TABLE]

So a feasible region defined by the upper bound condition is:

[TABLE]

∎

This is a qualitative result showing that the feasible region has a simple form: an interval. It is easy to see that this interval is non-empty for a very wide range of $n$ and $m$. Since the $g_{5r}^{-1}(0)$ term will grow almost like $\Theta(m^{2})$ whereas the other terms will grow almost like $\Theta(m^{\frac{1}{2}})$. In the following theorem we will formally prove the non-emptiness of the interval by picking the smallest $n$ in the feasible interval (given a $m$). We now give more context on this theorem.

In the following theorem, we use the above result (theorem 4) to answer a practical question that has strong implications for algorithmic design. Given a $m$, we want to find a $n$ (as small as possible) such that $\epsilon$ is below a tolerance:

[TABLE]

Where $p>0$ is a percentage that we want to upper bound the error with. A good answer to this question will resolve some crucial algorithmic design questions. Our strategy is simple. Given a target $\frac{p}{100}$, we will solve our upper bound from equation 39 for $n$, keeping everything else fixed.

Definition 22** ($p_{\max}$).**

Let $\mathcal{F}$ be some distribution family having parameters $\alpha,\beta,\gamma,k,n$. For all distributions $\epsilon(\mu)$ in this family, let the cdf be continuous and the kurtosis be uniformly bounded:

[TABLE]

Given $\mu_{i:m}\in[\gamma,\infty)$ $i=1,\ldots,m$. Then define:

[TABLE]

Definition 23** ($m_{\min}$).**

Let $\mathcal{F}$ be some distribution family having parameters $\alpha,\beta,\gamma,k,n$. For all distributions $\epsilon(\mu)$ in this family, let the cdf be continuous and the kurtosis be uniformly bounded:

[TABLE]

Given $\mu_{i:m}\in[\gamma,\infty)$ $i=1,\ldots,m$ and $0<p\leq p_{\max}$ (required for first piecewise invertibility of $\epsilon^{\mathrm{UB3}}_{\mathrm{BE}}$ at $\frac{T^{*}}{2}$), define

[TABLE]

and define $m_{\min}(p;\alpha,\beta,\gamma,k)$ as the infinimum $m$ that satisfies the equations

[TABLE]

Theorem 5** (Upper Bound on $n$ for Given Tolerance).**

Let $\mathcal{F}$ be some distribution family having parameters $\alpha,\beta,\gamma,k,n$. For all distributions $\epsilon(\mu)$ in this family, let the cdf be continuous and the kurtosis be uniformly bounded:

[TABLE]

Given $\mu_{i:m}\in[\gamma,\infty)$ $i=1,\ldots,m$ and $0<p\leq p_{\max}(;\alpha,\beta,\gamma,k)$ and let $m\in\mathbb{Z}^{+}$ be such that:

[TABLE]

And let:

[TABLE]

And let:

[TABLE]

We have

[TABLE]

Proof.

We consider $n$ of the form:

[TABLE]

Where $C>0$. Then we derive equivalent conditions for the conditions of theorem 4 (minimum mean estimation: error convergence rate). Consider the condition:

[TABLE]

Next consider the condition

[TABLE]

If $\frac{2^{\frac{5}{3}}C^{\frac{1}{3}}}{\pi C_{6}^{\frac{2}{3}}}m^{\frac{3}{2}}\leq 1$ the inequality is trivially satisfied. If $\frac{2^{\frac{5}{3}}C^{\frac{1}{3}}}{\pi C_{6}^{\frac{2}{3}}}m^{\frac{3}{2}}>1$ then the inequality is equivalent to:

[TABLE]

So we replace the original condition with the above stronger condition. Next, consider the condition

[TABLE]

Next consider the condition

[TABLE]

Finally consider

[TABLE]

Because $m^{\frac{1}{2}}>0$. And $-\log\left(\frac{\pi^{\frac{1}{2}}C_{6}^{\frac{1}{3}}}{2^{\frac{5}{6}}C^{\frac{1}{6}}m^{\frac{3}{4}}}\right)>0$ when $m^{\frac{3}{4}}>\frac{\pi^{\frac{1}{2}}C_{6}^{\frac{1}{3}}}{2^{\frac{5}{6}}C^{\frac{1}{6}}}$ (which is equivalent to condition 54)

[TABLE]

Hence if conditions 51, 52, 53, 54 and 55 are satisfied and also:

[TABLE]

then all the requirements for theorem 4 (minimum mean estimation: error convergence rate) are satisfied with $n=Cm^{\frac{1}{2}}$. Thus:

[TABLE]

is a constant for all $m,n$. We want to express this constant as a percentage $p>0$:

[TABLE]

We will now plug this requirement into the previously derived set of conditions. Consider condition 53:

[TABLE]

Consider condition 51

[TABLE]

Next consider condition 52.

[TABLE]

Next consider equation 54:

[TABLE]

Finally consider condition 55.

[TABLE]

Since conditions 57, 58, 59, 61 and 62 are satisfied by assumption and also:

[TABLE]

We have:

[TABLE]

∎

Appendix E MCPAM Proofs and Experimental Details

E.1 Theory

Theorem 6** (MCPAM Guarantees $k\geq 1$).**

We analyze MCPAM inner loop (line 5), with $\tau=0,n_{\max}=\infty$. Let $\Delta=\min|\mathrm{Ecc}(\bar{x}_{i_{1}})-\mathrm{Ecc}(\bar{x}_{i_{2}})|$, with $\bar{x}_{i_{1}},\bar{x}_{i_{2}}$ taken over all possible k-points generatable from $\{x_{i}\}$. Let $\Delta>0$ and be independent of $m$, let the distributions $\mathrm{Ecc}(\bar{x})$ be supported on $\mathbb{R}$. Let $C$ be the runtime cost of loop 5. Then loop $5$ is guaranteed to terminate and (expectation over the $\{T_{i}\}_{i=1}^{m}$)

[TABLE]

Upon exit, with high probability, either one of the following holds:

• •

we have found a smaller $\mathrm{Ecc}$ than $\mathrm{Ecc}(x^{\text{cur}})$

• •

there are no smaller $\mathrm{Ecc}$ than $\mathrm{Ecc}(x^{\text{cur}})$ in the swap set constructed from $x^{\text{cur}}$

And the confidence interval $[\mathrm{\widehat{E}cc}_{\text{lo}}(x^{\text{new}}),\mathrm{\widehat{E}cc}_{\text{hi}}(x^{\text{new}})]$ has true coverage of $\geq\alpha$.

For the distributed version, computational cost is $\mathcal{O}(\frac{km}{c})$ and communication cost is $\Theta(1)$.

Proof.

It is easy to show via SLLN that sample mean and sample variance converge almost surely to the mean and variance respectively, for r.v with finite variance. This gives that the width of the confidence intervals of the $\mathrm{\widehat{E}cc}$ go to zero almost surely as $n$ increases. Since $\Delta>0$, the intervals become non-overlapping and either line 9 or 11 of MCPAM must be satisfied a.s . Hence we will exit from the loop of line 5.

Let $p(n)$ be the probability of exiting from loop 5. We have shown $\lim_{n\rightarrow\infty}p(n)=1$. We are increasing $n$ in steps of 10. Let $F$ be number of iterations of loop 5 before exit. This is a non-homogenous geometric variable. Let $C(F)$ denote the runtime cost. We have per iteration cost $kmn=km10^{2+f}$. By summing the geometric progression $km10^{3},km10^{4},\ldots$ we get:

[TABLE]

The $\operatorname*{\mathbb{E}}(F)$ sum is dependent on the data distribution and independent of $m$. It is easy to see the convergence of this sum, by applying the ratio test. Let $T_{i}$ be the $i^{\text{th}}$ term.

[TABLE]

Hence the average runtime is $\operatorname*{\mathbb{E}}(C(F))=\mathcal{O}(km)$. We believe the sum $\operatorname*{\mathbb{E}}(F)$ will be upper bounded by $\frac{1}{p(10^{3})}$.

The procedure we follow in loop 5 is termed sequential Monte Carlo. It is well studied and has a convergence result [17] similar to the CLT. The confidence interval $[\mathrm{\widehat{E}cc}_{\text{lo}}(\bar{x}),\mathrm{\widehat{E}cc}_{\text{hi}}(\bar{x})]$ contains $\mathrm{Ecc}(\bar{x})$ with $1-\alpha$ probability as $m$ increases. Finally, note that condition 9 is checking (across $il$) $\mathrm{Ecc}(x^{\text{cur}})<\mathrm{Ecc}(x_{il})$ for the box confidence region of $[\mathrm{\widehat{E}cc}_{\text{lo}}(x^{\text{cur}}),\mathrm{\widehat{E}cc}_{\text{hi}}(x^{\text{cur}})]\times[\mathrm{\widehat{E}cc}_{\text{lo}}(x_{il}),\mathrm{\widehat{E}cc}_{\text{hi}}(x_{il})]$. The box confidence region is a superset of the actual $1-\alpha$ ellipsoidal confidence region, even in the correlated means case. The result on the distributed version is immediate from the above. ∎

Theorem 7** (MCPAM Guarantees $k=1$).**

We analyze 1-medoid MCPAM, with $\tau=0,n_{\max}=\infty$. If the conditions of 6 hold with $k=1$, Then MCPAM is guaranteed to terminate and $\operatorname*{\mathbb{E}}C=\mathcal{O}(m)$ Upon exit, we have found the sample medoid with high probability and the confidence interval around the estimate has true coverage of $\geq\alpha$.

Proof.

In the 1-medoid case, we finally exit loop 5 only via condition 9. The arguments of theorem 6 apply essentially unchanged. When we satisfy line 9, $\mathrm{\widehat{E}cc}(x^{\text{cur}})$ now has the lowest $\mathrm{Ecc}$ with probability $1-\alpha$. ∎

E.2 Experiments

E.2.1 Datasets

We use the collection of datasets provided in [14] for most of our evaluation. S1, S2, S3, S4 are datasets of size $5000$ in two dimensions with increasing overlap among a cluster. For ex- S4 will have significantly higher overlap among the clusters compared to S1, S2 and S3. They all have 15 clusters. Leaves is taken from [30] it contains 1600 rows with 64 attributes. There are 100 clusters. letter1 to letter4 is borrowed from [15] each of them have around 4600 rows with 16 attributes with increasing overlap among the classes. They all have 26 clusters.

For large scale run we used Foursquare checkin dataset [45] which contains around 32 million rows.

In addition we have generated two non-Euclidean datasets.

M1,M2,M3,M4: These are synthetic datasets with mixed (numeric and categorical) attributes. Each dataset consists of 2 numeric attributes and 1 categorical attribute with 2 levels. There are 3200 points in total and 32 clusters. Each cluster is a hybrid distribution, the numeric attributes are drawn from uncorrelated multivariate gaussian. The categorical attribute follows a Bernoulli distributon. The datasets M1,M2,M3,M4 are in increasing order of overlap. We use Gowers distance [19] with equal weights.

BillionOne: This is a synthetic, mixed (numeric and categorical attributes) dataset with a billion points. We have 5 real valued attributes, 1 categorical attribute with 24 levels. Each cluster is a hybrid distribution, the numeric attributes are drawn from an uncorrelated multivariate gaussian. The categorical attribute is drawn from a categorical distribution, aka generalized binomial distribution. There are 24 such clusters. We use Gowers distance [19] with equal weights.

E.2.2 Experimental Setup

Software Setup: For PAM we used R’s cluster package: cluster_2.0.3, R version 3.2.3 . For dbscan we used R’s dbscan package: dbscan_1.1-2 . The OS was Ubuntu 16.04 .

Hardware Setup: For small scale time and memory comparison we used commodity 16 GB RAM laptop, with a 6th Generation Intel i7 processor. For Foursquare dataset we used a local cluster of 4 commodity machines. Each with 32 GB RAM and Intel Core i7 CPUs, connected by a 1 Gigabit network. Each machine was running multiple workers, but the workers were isolated in different userspaces, so no two workers were affecting each other despite running on the same machine. For 1 Billion data point run we used 4 C5.4x large each with four workers and a C5.2x large as master.

Distribution Setup: For distribution, we implemented a master worker topology. For which we use Flask to create REST API endpoints and Redis as Message Broker for making asynchronous requests. Given $c$ workers, the data is partitioned into $c$ chunks of $m/c$ points. In our implementation the Master does out of core random sampling. The data is on the hard disk, but is never loaded into memory. Alternatively, it is also possible for master to not have access to any data just knowledge of how many data points each worker has is sufficient.

E.2.3 Experimental Results (Contd.)

Table 5 compares clustering cost between PAM and MCPAM.

Table 2 details the scaling of MCPAM on BillionOne data set.