Probabilistic smallest enclosing ball in high dimensions via subgradient sampling

TL;DR

This paper introduces a new algorithm for the probabilistic smallest enclosing ball problem in high-dimensional spaces, reducing computational complexity and enabling applications in kernel-based shape fitting methods.

Contribution

It presents a novel combination of sampling and stochastic subgradient descent to improve the efficiency of pSEB algorithms in high dimensions.

Findings

Reduces exponential dependence on dimension to linear in the algorithm.

Extends support vector data description (SVDD) to probabilistic cases.

Applicable to shape fitting in high-dimensional Hilbert spaces.

Abstract

We study a variant of the median problem for a collection of point sets in high dimensions. This generalizes the geometric median as well as the (probabilistic) smallest enclosing ball (pSEB) problems. Our main objective and motivation is to improve the previously best algorithm for the pSEB problem by reducing its exponential dependence on the dimension to linear. This is achieved via a novel combination of sampling techniques for clustering problems in metric spaces with the framework of stochastic subgradient descent. As a result, the algorithm becomes applicable to shape fitting problems in Hilbert spaces of unbounded dimension via kernel functions. We present an exemplary application by extending the support vector data description (SVDD) shape fitting method to the probabilistic case. This is done by simulating the pSEB algorithm implicitly in the feature space induced by the kernel function.

Equations

$$f(c)=\sum _{i=1}^{N}m(c,P_i),$$

$$m(A,B)=\{\begin{array}{ll}{\max }_{a\in A,b\in B}\Vert a-b\Vert & \text{if }A\ne B\\ 0& \text{if }A=B\end{array}$$

$$g(c_i)=\sum _{j=1}^N\frac{c_i-p_j}{\Vert c_i-p_j\Vert }\cdot \text{mathds}{1}_{c_i\ne p_j}\in \partial f(c_i),$$

$${\mathbb{E}}_{\tilde{c}}\left[f(\tilde{c})-f\left(c^{\ast }\right)\right]\le N\cdot \frac{R^2+(\ell +1)s^2}{2(\ell +1)s},$$

$${\mathbb{E}}_{\tilde{c}}\left[f(\tilde{c})-f\left(c^{\ast }\right)\right]\le \frac{NR}{\sqrt{\ell +1}}\le \frac{Cf(c^{\ast })}{\sqrt{\ell +1}}$$

$$\mathbf{Pr}\left[{\sum }_{P\in \mathcal{P}}m\left(a,P\right)>\left(1+\epsilon \right){\sum }_{P\in \mathcal{P}}m\left(\hat{c},P\right)\right]\le |\mathcal{C}|\cdot e^{-{\epsilon }^2|\mathcal{Q}|/64}.$$

$$c^{\ast }\in {\mathrm{argmin}}_{c\in {\mathbb{R}}^d}{\mathbb{E}}_X\left[m(c,X)\right],$$

$$(1-\epsilon )\cdot {\mathbb{E}}_X\left[\sum _{p\in X}\Vert c-p\Vert \right]\le {\mathbb{E}}_X\left[m(c,X)\right]\le {\mathbb{E}}_X\left[\sum _{p\in X}\Vert c-p\Vert \right],$$

$${\mathbb{E}}_X\left[m(c,X)\right]={\sum }_{P\ne \varnothing }\mathbf{Pr}\left[X=P\right]\cdot m(c,P),$$

$${\mathbb{E}}_X\left[m(\tilde{c},X)\right]\le (1+\epsilon ){\min }_{c\in {\mathbb{R}}^d}{\mathbb{E}}_X\left[m(c,X)\right].$$

$$c^{\ast }\in {\mathrm{argmin}}_{c\in \mathcal{H}}{\max }_{p\in P}\Vert c-\phi (p)\Vert ={\mathrm{argmin}}_{c\in \mathcal{H}}m(c,\phi (P)),$$

$$c^{\ast }\in {\mathrm{argmin}}_{c\in \mathcal{H}}{\mathbb{E}}_X\left[m(c,\phi (X))\right],$$

$$c={\sum }_{q_{u,v}\in Q}{\gamma }_{u,v}\phi (q_{u,v})$$

$$c_{i+1}=c_i-s\cdot \frac{c_i-\phi (q_{j,k})}{\Vert c_i-\phi (q_{j,k})\Vert }=\left(1-\frac{s}{\Vert c_i-\phi (q_{j,k})\Vert }\right)\cdot c_i+\frac{s}{\Vert c_i-\phi (q_{j,k})\Vert }\cdot \phi (q_{j,k}).$$

$${\Vert c_i-\phi (q)\Vert }^2$$

$$={\Vert {\sum }_{w=0}^i{\gamma }_w\phi (q_w)\Vert }^2+{\Vert \phi (q)\Vert }^2-2{\sum }_{w=0}^i{\gamma }_w⟨\phi (q_w),\phi (q)⟩$$

$$={\sum }_{w=0}^i{\sum }_{w^{\prime }=0}^i{\gamma }_w{\gamma }_{w^{\prime }}K(q_w,q_{w^{\prime }})+K(q,q)-2{\sum }_{w=0}^i{\gamma }_{w}K(q_w,q),$$

$${\mathbb{E}}_X\left[m(\tilde{c},\phi (X))\right]\le (1+\epsilon ){\min }_{c\in \mathcal{H}}{\mathbb{E}}_X\left[m(c,\phi (X))\right],$$

$$m(A,C)=\Vert a-c\Vert \le \Vert a-b\Vert +\Vert b-c\Vert \le m(A,B)+m(B,C),$$

$$m(A,C)=0+m(A,C)\le m(A,B)+m(B,C).$$

$$f({\lambda }_{1}x+{\lambda }_{2}y)$$

$$=\sum _{i=1}^N\Vert {\lambda }_1(x-p_i^{\ast })+{\lambda }_2(y-p_i^{\ast })\Vert \le {\lambda }_1\sum _{i=1}^N\Vert x-p_i^{\ast }\Vert +{\lambda }_2\sum _{i=1}^N\Vert y-p_i^{\ast }\Vert$$

$$\le {\lambda }_1\sum _{i=1}^N\underset{p_i\in P_i}{\max }\Vert x-p_i\Vert +{\lambda }_2\sum _{i=1}^N\underset{p_i\in P_i}{\max }\Vert y-p_i\Vert ={\lambda }_{1}f(x)+{\lambda }_{2}f(y)$$

$$|f(x)-f(y)|$$

$$\le \sum _{i=1}^{N}m(x,y)=N\Vert x-y\Vert .$$

$$⟨\frac{c_i-p_j}{\Vert c_i-p_j\Vert }\cdot \text{mathds}{1}_{c_i\ne p_j},c_i-c^{\ast }⟩\ge m(c_i,P_j)-m(c^{\ast },P_j).$$

$$\frac{⟨c_i-p_j,c_i-c^{\ast }⟩}{\Vert c_i-p_j\Vert }=\frac{⟨c_i-p_j,c_i-p_j+p_j-c^{\ast }⟩}{\Vert c_i-p_j\Vert }$$

$$=\frac{⟨c_i-p_j,c_i-p_j-(c^{\ast }-p_j)⟩}{\Vert c_i-p_j\Vert }=\frac{⟨c_i-p_j,c_i-p_j⟩-⟨c_i-p_j,c^{\ast }-p_j⟩}{\Vert c_i-p_j\Vert }$$

$$\stackrel{\text{CSI}}{\ge }\frac{{\Vert c_i-p_j\Vert }^2-\Vert c_i-p_j\Vert \Vert c^{\ast }-p_j\Vert }{\Vert c_i-p_j\Vert }\ge \Vert c_i-p_j\Vert -\Vert c^{\ast }-p_j^{\ast }\Vert =m(c_i,P_j)-m(c^{\ast },P_j),$$

$$⟨\sum _{j=1}^N\frac{c_i-p_j}{\Vert c_i-p_j\Vert }\cdot \text{mathds}{1}_{c_i\ne p_j},c_i-c^{\ast }⟩$$

$$\ge \sum _{j=1}^N\left(m(c_i,P_j)-m(c^{\ast },P_j)\right)=f(c_i)-f(c^{\ast }).$$

$$\mathbb{E}\left[\tilde{g}(c_i)\right]$$

$${\mathbb{E}}_{\tilde{g}}\left[{\Vert c_{i+1}-c^{\ast }\Vert }^2\mid c_i\right]$$

$$={\mathbb{E}}_{\tilde{g}}\left[{\Vert c_i-c^{\ast }\Vert }^2+{\Vert s\tilde{g}(c_i)\Vert }^2-2⟨s\tilde{g}(c_i),c_i-c^{\ast }⟩\mid c_i\right]$$

$$\le {\mathbb{E}}_{\tilde{g}}\left[{\Vert c_i-c^{\ast }\Vert }^2\mid c_i\right]+s^2-2s⟨{\mathbb{E}}_{\tilde{g}}\left[\tilde{g}(c_i)\mid c_i\right],c_i-c^{\ast }⟩$$

Full text

1 Introduction

The (probabilistic) smallest enclosing ball (pSEB) problem in $\mathbb{R}^{d}$ is to find a center that minimizes the (expected) maximum distance to the input points (see Definition 3.1). It occurs often as a building block for complex data analysis and machine learning tasks like estimating the support of high dimensional distributions, outlier detection, novelty detection, classification and robot gathering [9, 27, 28, 30]. It is thus very important to develop highly efficient approximation algorithms for the base problem. This involves reducing the number of points but also keeping the dependence on the dimension as low as possible. Both objectives will be studied in this paper. We will focus on a small dependence on the dimension. This is motivated as follows.

Kernel methods are a common technique in machine learning. These methods implicitly project the $d$-dimensional input data into much larger dimension $D$ where simple linear classifiers or spherical data fitting methods can be applied to obtain a non-linear separation or non-convex shapes in the original $d$-dimensional space. The efficiency of kernel methods is usually not harmed. Despite the large dimension $D\gg d$, most important kernels, and thus inner products and distances in the $D$-dimensional space, can be evaluated in $O(d)$ time [25].

In some cases, however, a proper approximation relying on sampling and discretizing the ambient solution space may require a polynomial or even exponential dependence on $D$. The algorithm of Munteanu et al.[22] is the only fully polynomial time approximation scheme (FPTAS) and the fastest algorithm to date for the pSEB problem in fixed dimension. However, it suffers from the stated problems. In particular, the number of realizations sampled by their algorithm had a linear dependence on $D$ stemming from a ball-cover decomposition of the solution space. The actual algorithm made a brute force evaluation (on the sample) of all centers in a grid of exponential size in $D$. This is prohibitive in the setting of kernel methods since the implicit feature space may have infinite dimension. Even if it is possible to exploit the up to $n$-dimensional subspace spanned by $n$ points in infinite dimensions, we would still have $D=n\gg d$ leading to exponential time algorithms.

To make the probabilistic smallest enclosing ball algorithm viable in the context of kernel methods and generally in high dimensions, it is highly desirable to reduce the dependence on the dimension to a small polynomial occurring only in evaluations of inner products and distances between two (low dimensional) vectors.

1.1 Related work

Probabilistic smallest enclosing ball

The study of probabilistic clustering problems was initiated by Cormode and McGregor [11]. They developed approximation algorithms for the probabilistic settings of $k$-means, $k$-median as well as $k$-center clustering. For the metric $1$-center clustering problem their results are bi-criteria $O(1)$-approximation algorithms with a blow-up on the number of centers. Guha and Munagala [14] improved the previous work by giving $O(1)$-approximations that preserve the number of centers. Munteanu et al.[22] gave the first fully polynomial time $(1+\varepsilon)$-approximation scheme (FPTAS) for the probabilistic Euclidean $1$-center, i.e., the probabilistic smallest enclosing ball problem, in fixed dimensions. The algorithm runs in linear time for sampling a constant number of realizations. Solving the subsampled problem takes only constant time, though exponential in the dimension. This yields a total running time of roughly $O(nd/\varepsilon^{O(1)}+1/\varepsilon^{O(d)})$. Based on $\varepsilon$-kernels for probabilistic data [16], Huang and Li [15] generalized the $(1+\varepsilon)$-approximation of [22] to a polynomial time approximation scheme (PTAS) for Euclidean $k$-center in $\mathbb{R}^{d}$ for fixed constants $k$ and $d$. We note that the running time of the algorithm in [15] grows as a double exponential function of the dimension and it is unclear how to reduce this. We thus base our work on the FPTAS of [22] and reduce its exponential dependence to linear.

Sampling techniques for $1$-median

The work of [22] showed a reduction from the probabilistic smallest enclosing ball to $1$-median problems on (near)-metric spaces. We thus review relevant results on sampling techniques for metric $1$-median. Bădoiu et al.[7] showed that with constant probability, the span of a uniform sample of a constant number of input points contains a $(1+\varepsilon)$-approximation for the $1$-median. This was used to construct a set of candidate solutions of size $O(2^{1/\varepsilon^{O(1)}}\log n)$. A more refined estimation procedure based again on a small uniform sample was used by Kumar et al.[19] to reduce this to $O(2^{1/\varepsilon^{O(1)}})$. Indyk and Thorup [17, 29] showed that a uniform sample of size $O(\log n/\varepsilon^{2})$ is sufficient to approximate the discrete metric $1$-median on $n$ points within a factor of $(1+\varepsilon)$. Ackermann et al.[1] showed how this argument can be adapted to doubling spaces which include the continuous Euclidean space. We adapt these ideas to find a $(1+\varepsilon)$-approximation to the best center in our setting.

Stochastic subgradient descent

One quite popular and often only implicitly used technique in the coreset literature is derived from convex optimization, see [21]. One of the first results of that kind is given in the uniform sampling algorithm of Bădoiu et al.[7] for $1$-median stated above. In each iteration a single point is sampled uniformly. With high probability moving the current center towards that point for a carefully chosen step size improves the solution. Each step can be seen as taking a descent towards a uniformly random direction from the subgradient which equals roughly the sum of directions to all points. Another example is the coreset construction for the smallest enclosing ball problem by Bădoiu and Clarkson [5], where the next point included in the coreset is the one maximizing the distance to the current best center. The direction taken towards that point is again related to the subgradient of the objective function at the current center position. The authors also gave a more explicit application of subgradient descent to the problem with a slightly larger number of iterations. More recently, Cohen et al.[10] developed one of the fastest $(1+\varepsilon)$-approximation algorithms to date for the geometric median problem via stochastic subgradient methods.

Kernel methods in machine learning

Kernel functions simulate an inner product space in large or even unbounded dimensions but can be evaluated via simple low dimensional vector operations in the original dimension of input points [26]. This enables simple spherical shape fitting via a smallest enclosing ball algorithm in the high dimensional feature space, which implicitly defines a more complex and even non-convex shape in the original space. The smallest enclosing ball problem in kernel spaces is well-known as the support vector data description (SVDD) by Tax and Duin [28]. A more subtle connection between (the dual formulations of) several kernel based methods in machine learning and the smallest enclosing ball problem was established by Tsang et al.[30].

1.2 Our contributions and outline

• •

We extend in Section 2 the geometric median in Euclidean space to the more general set median problem. It consists of finding a center $c\in\mathbb{R}^{d}$ that minimizes the sum of maximum distances to sets of points in a given collection of $N$ point sets. We show how to solve this problem via estimation and sampling techniques combined with a stochastic subgradient descent algorithm.

• •

The elements in the collection are sets of up to $n$ points in $\mathbb{R}^{d}$. We discuss in Section 2.1 the possibility of further reducing their size. In the previous work [22] they were summarized via strong coresets of size $1/\varepsilon^{\Theta(d)}$ for constant dimension $d$. This is not an option in high dimensions where, e.g. $d\approx n$. Reviewing the techniques of Agarwal and Sharathkumar [2] we can show that no reduction below $\min\{n,\exp(d^{1/3})\}$ is possible unless one is willing to sacrifice an additional approximation factor of roughly $\sqrt{2}$. However, we discuss the possibility to achieve roughly a factor $(\sqrt{2}+\varepsilon)$-approximation in streaming via the blurred-ball-cover [2] of size $O(1/\varepsilon^{3}\cdot\log 1/\varepsilon)$, and in an off-line setting via weak coresets [5, 6] of size $O(1/\varepsilon)$.

• •

We show in Section 3.1 how this improves the previously best FPTAS for the probabilistic smallest enclosing ball problem from $O(dn/\varepsilon^{3}\cdot\log 1/\varepsilon+1/\varepsilon^{O(d)})$ to $O(dn/\varepsilon^{4}\cdot\log^{2}1/\varepsilon)$. In particular the dependence on the dimension $d$ is reduced from exponential to linear and more notably occurs only in distance evaluations between points in $d$-dimensional Euclidean space, but not in the number of sampled points nor in the number of candidate centers to evaluate.

• •

This enables in Section 3.2 working in very high $D$-dimensional Hilbert spaces whose inner products and distances are given implicitly via positive semidefinite kernel functions. These functions can be evaluated in $O(d)$ time although $D$ is large or even unbounded depending on the kernel function. As an example we extend the well-known support vector data description (SVDD) method to the probabilistic case. SVDD is equivalent to the smallest enclosing ball problem in the implicit high-dimensional feature space.

Please find the missing proofs in Appendix A.

1.3 General notation

We denote the set of positive integers up to $n\in\mathbb{N}$ by $[n]=\{1,\ldots,n\}$. For any vectors $x,y\in\mathbb{R}^{d}$ we denote their inner product by $\langle x,y\rangle=x^{T}y=\sum\nolimits_{i=1}^{d}x_{i}y_{i}$ and the Euclidean norm by $\left\|x\right\|=\sqrt{\langle x,x\rangle}=(\sum\nolimits_{i=1}^{d}x_{i}^{2})^{1/2}$. The Cauchy-Schwarz inequality (CSI) states that $\langle x,y\rangle\leq\left\|x\right\|\left\|y\right\|$. For any convex function $f\colon\mathbb{R}^{d}\rightarrow\mathbb{R}$ we denote by $\partial f(x)=\left\{g\in\mathbb{R}^{d}\mid\forall y\in\mathbb{R}^{d}\colon f(x)-f(y)\leq\langle g,x-y\rangle\right\}$ the subdifferential, i.e., the set of subgradients of $f$ at $x$. For any event $\mathcal{E}$ let the indicator function be $\mathds{1}_{\mathcal{E}}=1$ if $\mathcal{E}$ happens and [math] otherwise. We denote by $B(c,r)=\{x\mid\left\|c-x\right\|\leq r\}$ the Euclidean ball centered at $c\in\mathbb{R}^{d}$ with radius $r\geq 0$. We assume the error parameter satisfies $0<\varepsilon<1/9$. All our results hold with constant probability, say $1/8$, which can be amplified to arbitrary $1-\eta$, $0<\eta<1$, by running ${O(\log 1/\eta)}$ independent repetitions and returning the minimum found.

2 A generalized median problem

The probabilistic smallest enclosing ball (pSEB) problem is to find a center that minimizes the expected maximum distance to points drawn from the input distributions (see Section 3.1 and Definition 3.1), and it can be reduced to two different types of $1$-median problems [22]. One of them is defined on the set of all non-empty locations in $\mathbb{R}^{d}$ where probabilistic points may appear, equipped with the Euclidean distance. The other is defined on the collection of all possible realizations of probabilistic point sets, and the distance measure between a center $c\in\mathbb{R}^{d}$ and a realization $P_{i}\subset\mathbb{R}^{d}$ is the maximum distance between the center and any of the realized points, i.e., $\max\nolimits_{p\in P_{i}}\left\|c-p\right\|$.

We begin our studies with a generalized median problem that we call the set median problem and covers both of these cases.

Definition 2.1** (set median problem).**

Let $\mathcal{P}=\{P_{1},\ldots,P_{N}\}$ be a family of finite non-empty sets where $\forall i\in[N]\colon P_{i}\subset\mathbb{R}^{d}$ and $n=\max\{|P_{i}|\mid i\in[N]\}$. The set median problem on $\mathcal{P}$ consists in finding a center $c\in\mathbb{R}^{d}$ that minimizes the cost function

[TABLE]

where $m(c,P_{i})=\max\nolimits_{p\in P_{i}}\left\|c-p\right\|$.

It was noted in [22] that the distance measure $m$ between any two sets $A,B\subset\mathbb{R}^{d}$ defined by the maximum distance $m(A,B)=\max\nolimits_{a\in A,b\in B}\left\|a-b\right\|$ is not a metric since for any non-singleton set $C\subset\mathbb{R}^{d}$ it holds that $m(C,C)>0$. We stress here that we consider only cases where, as in Definition 2.1, one of the sets $A=\{c\}$ is a singleton, and $B=P$ is an arbitrary non-empty set of points from $\mathbb{R}^{d}$. In order to directly apply results from the theory of metric spaces we thus simply define $m(A,B)=0$ whenever $A=B$.

Lemma 2.2**.**

Let $\mathcal{X}$ be the set of all finite non-empty subsets of $\mathbb{R}^{d}$. We define

[TABLE]

for any $A,B\in\mathcal{X}$. Then $(\mathcal{X},m)$ is a metric space.

Note that in case of singleton sets, the set median problem in Definition 2.1 is equivalent to the well-known Fermat-Weber problem also known as $1$-median or geometric median. Also, for $N=1$ it coincides with the smallest enclosing ball or $1$-center problem.

For both of these problems there are known algorithms based on the subgradient method from convex optimization. Bădoiu and Clarkson [5] gave a simple algorithm for approximating the 1-center within $(1+\varepsilon)$-error. Starting from an initial center, it is iteratively moved a little towards the input point that is furthest away. Note that a suitable subgradient at the current center points exactly into the opposite direction. More precisely if $q\in P$ is a point that is furthest away from the current center $c$ then $(c-q)/\left\|c-q\right\|\in\partial\max_{p\in P}\left\|c-p\right\|$. The algorithm can thus be interpreted as a subgradient descent minimizing $\max_{p\in P}\left\|c-p\right\|$. Weiszfeld’s algorithm [31, 32] was the first that solved the $1$-median problem within additive $O(\varepsilon)$-error in $O(1/\varepsilon)$ subgradient iterations. And indeed one of the fastest algorithms to approximate the geometric median problem to date relies on a stochastic subgradient descent [10]. The crucial ingredients to turn the additive error to a relative error are finding a suitable starting point that achieves a constant approximation and estimating its initial distance to the optimal solution. Another important step is a bound on the Lipschitz constant of the cost function, see also [3]. We will generalize these approaches to minimizing the cost function $f$ of the set median problem from Definition 2.1.

First note that $f$ is a convex function which implies that we can apply the theory of convex analysis and optimization. In particular, it implies that the subdifferential $\partial f(c)$ is non-empty for any center $c\in\mathbb{R}^{d}$ and $c$ is locally optimal if and only if $0\in\partial f(c)$. Moreover, any local optimum is also globally optimal by convexity [4, 23]. This implies that if we find a $(1+\varepsilon)$-approximation to a local minimum of the convex function, the convexity implies that it is a $(1+\varepsilon)$-approximation to the global minimum as well.

To see that $f$ is convex, note that the Euclidean norm is a convex function. Therefore the Euclidean distance to some fixed point is a convex function since every translation of a convex function is convex. The maximum of convex functions is a convex function and finally the sum of convex functions is again convex. We prove this claim for completeness.

Lemma 2.3**.**

The objective function $f$ of the set median problem (see Definition 2.1) is convex.

Next we bound the Lipschitz constant of the function $f$ by $N$. We may get a better bound if we limit the domain of $f$ to a ball of small radius centered at the optimal solution, but the Lipschitz constant cannot be bounded by $o(N)$ in general. We will see later how we can remove the dependence on $N$.

Lemma 2.4**.**

The objective function $f$ of the set median problem (see Definition 2.1) is $N$-Lipschitz continuous, i.e., $|f(x)-f(y)|\leq N\cdot\left\|x-y\right\|$ for all $x,y\in\mathbb{R}^{d}$.

We want to minimize $f$ via the subgradient method, see [23]. For that sake we need to compute a subgradient $g(c_{i})\in\partial f(c_{i})$ at the current center $c_{i}$. To this end we prove the following lemma.

Lemma 2.5**.**

Let $c_{i}\in\mathbb{R}^{d}$ be any center. For each set $P_{j}\in\mathcal{P}$, let $p_{j}\in P_{j}$ be a point with $\left\|c_{i}-p_{j}\right\|=m(c_{i},P_{j})$. We have

[TABLE]

i.e., $g(c_{i})$ is a valid subgradient of $f$ at $c_{i}$.

For brevity of presentation we omit the indicator function in any use of the above Lemma in the remainder of this paper and simply define $\left(c_{i}-p_{j}\right)/\left\|c_{i}-p_{j}\right\|=0$ whenever $c_{i}=p_{j}$.

The subgradient computation takes ${O(d{n}N)}$ time to calculate, since in each of the $N$ terms of the sum we maximize over $|P_{i}|\leq{n}$ distances in $d$ dimensions to find a point in $P_{i}$ that is furthest away from $c$. We are going to discuss the possibility of reducing ${n}$ later. For now we focus on removing the dependence on $N$. To this end we would like to replace the exact subgradient $g(c_{i})$ by a uniform sample of only one nonzero term which points into the right direction in expectation. We formalize this in the following lemma.

Lemma 2.6**.**

Let $c_{i}\in\mathbb{R}^{d}$ be any fixed center. For each set $P_{j}\in\mathcal{P}$, let $p_{j}\in P_{j}$ be a point with $\left\|c_{i}-p_{j}\right\|=m(c_{i},P_{j})$. Let $\tilde{g}(c_{i})$ be a random vector that takes the value $\tilde{g}(c_{i})=\left(c_{i}-p_{j}\right)/\left\|c_{i}-p_{j}\right\|$ for $j\in[N]$ with probability $1/N$ each. Then $\mathbb{E}\left[\left\|\tilde{g}(c_{i})\right\|^{2}\right]\leq 1$ and $\mathbb{E}\left[\tilde{g}(c_{i})\right]={g}(c_{i})/N$, where ${g}(c_{i})\in\partial f(c_{i})$ is the subgradient given in Lemma 2.5.

We can now adapt the deterministic subgradient method from [23] using the random unbiased subgradient of Lemma 2.6 in such a way that the result is in expectation a $(1+\varepsilon)$-approximation to the optimal solution. This method is presented in Algorithm 1. Given an initial center $c_{0}$, a fixed step size $s$, and a number of iterations $\ell$, the Algorithm iteratively picks a set $P_{j}\in\mathcal{P}$ uniformly at random and chooses a point $p_{j}\in P_{j}$ that attains the maximum distance to the current center. This point is used to compute an approximate subgradient via Lemma 2.6. The Algorithm finally outputs the best center found in all iterations.

The following theorem bounds in expectation the quality of the output that our subgradient algorithm returns. It is a probabilistic adaptation of a result in convex optimization [23].

Theorem 2.7**.**

Consider Algorithm 1 on input $\mathcal{P}=\left\{P_{1},\ldots,P_{N}\right\}$ for the set median problem with objective function $f$, see Definition 2.1. Let $c^{*}\in\operatorname{argmin}\limits_{c\in\mathbb{R}^{d}}f(c)$. Let $R=\left\|c_{0}-c^{*}\right\|$. Then

[TABLE]

where the expectation is taken over the random variable $\tilde{c}\in\operatorname{argmin}\limits_{c\in\{c_{i}\mid i=0,\ldots,\ell\}}f\left(c\right)$, i.e., the output of Algorithm 1.

Our aim now is to choose the parameters $\ell,s,$ and $c_{0}$ of Algorithm 1 in such a way that the bound given in Theorem 2.7 becomes at most $\varepsilon f(c^{*})$. To this end we can choose the initial center $c_{0}$ and bound its initial distance $R=\left\|c_{0}-c^{*}\right\|$ proportional to the average cost ${O(f(c^{*})/N)}$ with constant probability using a simple Markov argument.

Lemma 2.8**.**

Choose a set $P$ from $\mathcal{P}=\{P_{1},\ldots,P_{N}\}$ uniformly at random and let $c_{0}$ be an arbitrary point of $P$. Then for any constant $0<\delta_{1}<1$ it holds that $R=\left\|c_{0}-c^{*}\right\|\leq f(c^{*})/(\delta_{1}N)$ with probability at least $1-\delta_{1}$.

Assume that we know the value of $R$, and set the step size to $s=R/\sqrt{\ell+1}$, then Theorem 2.7 and Lemma 2.8 imply that for some constant $C$

[TABLE]

holds with constant probability. We thus only need to run the algorithm for $\ell\in{O(1/\varepsilon^{2})}$ iterations to get within $\varepsilon f(c^{*})$ error. But choosing this particular step size requires to know the optimal center in advance. To get around this, we attempt to estimate the average cost. More formally, we are interested in a constant factor approximation of $f(c^{*})/N$. It turns out that we can do this based on a small sample of the input sets unless our initial center is already a good approximation. But in the latter case we do not care about all the subsequent steps or step sizes, since we are already done after the initialization. The proof technique is originally from [19] and is adapted here to work in our setting with sets of points.

Lemma 2.9**.**

There exists an algorithm that based on a sample $\mathcal{S}\subseteq\mathcal{P}$ of size $|\mathcal{S}|=1/\varepsilon$ returns an estimate $\tilde{R}$ and an initial center $c_{0}$ in time $O(d{n}/\varepsilon)$ such that with constant probability one of the following holds:

1. a)

$\varepsilon f(c^{*})/N\leq\tilde{R}\leq\left(2/\varepsilon^{3}\right)\cdot f(c^{*})/N$* and $\left\|c_{0}-c^{*}\right\|\leq 8f(c^{*})/N$;* 2. b)

$\left\|c_{0}-c^{*}\right\|\leq 4\varepsilon f(c^{*})/N$.

Lemma 2.9 has the following consequence. Either the initial center $c_{0}$ is already a $(1+4\varepsilon)$-approximation, in which case we are done. Or we are close enough to an optimal solution and have a good estimate on the step size to find a $(1+4\varepsilon)$-approximation in a constant number of iterations.

Another issue that we need to take care of, is finding the best center in the last line of Algorithm 1 efficiently. We cannot do this exactly since evaluating the cost even for one single center takes time $O(d{n}N)$. However, we can find a point that is a $(1+\varepsilon)$-approximation of the best center in a finite set of candidate centers using a result from the theory of discrete metric spaces.

To this end we can apply our next theorem which is originally due to Indyk and Thorup in [17, 29] and adapted here to work in our setting. The main difference is that in the original work the set of input points and the set of candidate solutions are identical. In our setting, however, we have that the collection of input sets and the set of candidate solutions may be completely distinct and the distance measure is the maximum distance (see Lemma 2.2).

Theorem 2.10**.**

Let $\mathcal{Q}$ be a set of uniform samples with repetition from $\mathcal{P}$. Let $\mathcal{C}$ be a set of candidate solutions. Let $a\in\mathcal{C}$ minimize $\sum_{Q\in\mathcal{Q}}m(a,Q)$ and let $\hat{c}=\operatorname{argmin}_{c\in\mathcal{C}}f(c)$. Then

[TABLE]

Putting all pieces together we have the following Theorem.

Theorem 2.11**.**

Consider an input $\mathcal{P}=\{P_{1},\ldots,P_{N}\}$, where for every $i\in[N]$ we have $P_{i}\subset\mathbb{R}^{d}$ and $n=\max\{|P_{i}|\mid i\in[N]\}$. There exists an algorithm that computes a center $\tilde{c}$ that is with constant probability a $(1+\varepsilon)$-approximation to the optimal solution $c^{*}$ of the set median problem (see Definition 2.1). Its running time is ${O(d{n}/\varepsilon^{4}\cdot\log^{2}1/\varepsilon)}$.

2.1 On reducing the size of the input sets

We would like to remove additionally the linear dependence on $n$ for the maximum distance computations. This is motivated from the streaming extension of [22] where one aims at reading the input once in linear time and all subsequent computations should be sublinear, or preferably independent of $n$. To this end a grid based strong coreset of size $1/\varepsilon^{\Theta(d)}$ was used. However, here we focus on reducing the dependence on $d$, and exponential is not an option if we want to work in high dimensions. It turns out that without introducing an exponential dependence on $d$, we would have to lose a constant approximation factor. Pagh et al.[24] showed that if the coreset is a subset of the input and approximates furthest neighbor queries to within less than a factor of roughly $\sqrt{2}$, then it must consist of $\Omega\left(\min\{n,\exp(d^{1/3})\}\right)$ points. This was shown via a carefully constructed input point set of Agarwal and Sharathkumar [2] who used it to prove lower bounds on streaming algorithms for several extent problems.

In the next theorem, we review the techniques of the latter reference to show a slightly stronger result, namely no small data structure can exist for answering maximum distance queries to within a factor of less than roughly $\sqrt{2}$. In comparison to the previous results [2, 24], it is not limited to the streaming setting, and it is not restricted to subsets of the input.

Theorem 2.12**.**

Any data structure that, with probability at least $2/3$, $\alpha$-approximates maximum distance queries on a set $S\subset\mathbb{R}^{d}$ of size $|S|=n$, for $\alpha<\sqrt{2}\left(1-2/d^{1/3}\right)$, requires $\Omega\left(\min\{n,\exp\left(d^{1/3}\right)\}\right)$ bits of storage.

On the positive side it was shown by Goel et al.[13] that a $\sqrt{2}$-approximate furthest neighbor to the point $c\in\mathbb{R}^{d}$ can always be found on the surface of the smallest enclosing ball of the sets $P_{i}$, using linear preprocessing time $\tilde{O}\left(dn\right)$ and $\tilde{O}\left(d^{2}\right)$ query time. Thus, if we plug in the coresets of Bădoiu and Clarkson [5, 6] of size $O(1/\varepsilon)$ instead of the entire sets $P_{i}\in\mathcal{P}$ to evaluate $m(c,P_{i})$, we would have a sublinear time algorithm (in $n$) after reading the input via a $\sqrt{2}\left(1+\varepsilon\right)$-approximation to any query $m(c,P_{i})$. In a streaming setting the same bound can be achieved via the blurred-ball-cover of Agarwal and Sharathkumar [2] of slightly larger size ${O(1/\varepsilon^{3}\cdot\log 1/\varepsilon)}$. Otherwise, Goel et al.[13] have shown that using ${O(dn^{1+1/(1+\varepsilon)})}$ preprocessing time and $\tilde{O}\left(dn^{1/(1+\varepsilon)}\right)$ query time, one can obtain a $(1+\varepsilon)$-approximation for the furthest neighbor problem. Note that in this case the preprocessing time is already superlinear in $n$, and in particular the exponent is already larger than $1.7$ for $1+\varepsilon<\sqrt{2}$.

3 Applications

3.1 Probabilistic smallest enclosing ball

We apply our result to the probabilistic smallest enclosing ball problem, as given in [22]. In such a setting, the input is a set $\mathcal{D}=\{D_{1},\ldots,D_{n}\}$ of $n$ discrete and independent probability distributions. The $i$-th distribution $D_{i}$ is defined over a set of $z$ possible locations $q_{i,j}\in\mathbb{R}^{d}\cup\{\bot\}$, for $j\in\left[z\right]$, where $\bot$ indicates that the $i$-th point is not present in a sampled set, i.e., $q_{i,j}=\bot\Leftrightarrow\{q_{i,j}\}=\emptyset$. We call these points probabilistic points. Each location $q_{i,j}$ is associated with the probability $p_{i,j}$, such that $\sum_{j=1}^{z}p_{i,j}=1$, for every $i\in\left[n\right]$. Thus the probabilistic points can be considered as independent random variables $X_{i}$.

A probabilistic set $X$ consisting of probabilistic points is also a random variable, where for each random choice of indices $\left(j_{1},\ldots,j_{n}\right)\in\left[z\right]^{n}$ there is a realization $P_{\left(j_{1},\ldots,j_{n}\right)}=X\left(j_{1},\ldots,j_{n}\right)=\left(q_{1,j_{1}},\ldots,q_{n,j_{n}}\right)$. By independence of the distributions $D_{i}$, $i\in\left[n\right]$, it holds that $\mathbf{Pr}\left[X=P_{\left(j_{1},\ldots,j_{n}\right)}\right]=\prod_{i=1}^{n}p_{i,j_{i}}$.

The probabilistic smallest enclosing ball problem is defined as follows. Here we may assume that the distance of any point $c\in\mathbb{R}^{d}$ to the empty set is 0.

Definition 3.1**.**

*([22])

Let $\mathcal{D}$ be a set of $n$ discrete distributions, where each distribution is defined over $z$ locations in $\mathbb{R}^{d}\cup\left\{\bot\right\}$. The probabilistic smallest enclosing ball problem is to find a center $c^{*}\in\mathbb{R}^{d}$ that minimizes the expected smallest enclosing ball cost, i.e.,*

[TABLE]

where the expectation is taken over the randomness of $X\sim\mathcal{D}$.

The authors of [22] showed a reduction of the probabilistic smallest enclosing ball problem to computing a solution for the set median problem of Definition 2.1. Their algorithm distinguishes between two cases. In the first case the probability of obtaining a nonempty realization $P\neq\emptyset$ is small, more formally $\sum_{q_{i,j}\in{Q}}p_{i,j}\leq\varepsilon$, where $Q=\{q_{i,j}\mid q_{i,j}\neq\bot,i\in[n],j\in[z]\}$, and thus we have little chance of gaining information by sampling realizations. However, it was shown that

[TABLE]

where $\mathbb{E}_{X}\left[\sum_{p\in X}\left\|c-p\right\|\right]=\sum_{i,j}p_{i,j}\cdot\left\|c-q_{i,j}\right\|$ is a weighted version of the deterministic 1-median problem, cf. [11], and thus also a weighted instance of the set median problem. In the second case, the probability that a realization contains at least one point is reasonably large. Therefore by definition of the expected value and $m(c,\emptyset)=0$ we have

[TABLE]

which is a weighted version of the set median problem. Depending on these two cases, we sample a number of elements, non-empty locations or non-empty realizations, and solve the resulting set median problem using the samples in Theorem 2.11 for computing the approximate subgradients.

Algorithm 2 adapts this framework. It differs mainly in three points from the previous algorithm of [22]. First, the number of samples had a dependence on $d$ hidden in the $O$-notation. This is not the case any more. Second, the sampled realizations are not sketched via coresets of size $1/\varepsilon^{\Theta(d)}$ any more, as discussed in Section 2.1. Third, the running time of the actual optimization task is reduced via Theorem 2.11 instead of an exhaustive grid search.

Theorem 3.2**.**

Let $\mathcal{D}$ be a set of $n$ discrete distributions, where each distribution is defined over $z$ locations in $\mathbb{R}^{d}\cup\{\bot\}$. Let $\tilde{c}\in\mathbb{R}^{d}$ denote the output of Algorithm 2 on input $\mathcal{D}$, and let the approximation parameter be $\varepsilon<1/9$. Then with constant probability the output is a $(1+\varepsilon)$-approximation for the probabilistic smallest enclosing ball problem. I.e., it holds that

[TABLE]

The running time of Algorithm 2 is ${O(dn\cdot(z/\varepsilon^{3}\cdot\log 1/\varepsilon+1/\varepsilon^{4}\cdot\log^{2}1/\varepsilon))}$.

Comparing to the result of [22], the running time is reduced from $O(dnz/\varepsilon^{O(1)}+1/\varepsilon^{O(d)})$ to $O(dnz/\varepsilon^{O(1)})$, i.e., our dependence on the dimension $d$ is no longer exponential but only linear. Note, in particular, that the factor of $d$ plays a role only in computations of distances between two points in $\mathbb{R}^{d}$. Further the sample size and the number of centers that need to be evaluated do not depend on the dimension $d$ any more. This will be crucial in the next application.

3.2 Probabilistic support vector data description

Now we turn our attention to the support vector data description (SVDD) problem [28] and show how to extend it to its probabilistic version. To this end, let $K\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a positive semidefinite kernel function. It is well known by Mercer’s theorem, cf. [26], that such a function implicitly defines the inner product of a high dimensional Hilbert space $\mathcal{H}$, say $\mathbb{R}^{D}$ where $D\gg d$. This means we have $K(x,y)=\langle\varphi(x),\varphi(y)\rangle$, where $\varphi\colon\mathbb{R}^{d}\rightarrow\mathcal{H}$ is the so called feature mapping associated with the kernel. Examples for such kernel functions include polynomial transformations of the standard inner product in $\mathbb{R}^{d}$ such as the constant, linear or higher order polynomial kernels. In these cases $D$ remains bounded but grows as a function of $d$ raised to the power of the polynomials’ degree. Other examples are the exponential, squared exponential, Matérn, or rational quadratic kernels, which are transformations of the Euclidean distance between the two low dimensional vectors. The dimension $D$ of their implicit feature space is in principle unbounded. Despite the large dimension $D\gg d$, all these kernels can be evaluated in time $O(d)$ [25].

It is known that the SVDD problem is equivalent to the smallest enclosing ball problem in the feature space induced by the kernel function [30], i.e., given the kernel $K$ with implicit feature mapping $\varphi\colon\mathbb{R}^{d}\rightarrow\mathcal{H}$, and an input set $P\subseteq\mathbb{R}^{d}$, the task is to find

[TABLE]

where $\varphi(P)=\{\varphi(p)\mid p\in P\}$.

Now we extend this to the probabilistic setting as we did in Section 3.1. The input is again a set $\mathcal{D}$ of $n$ discrete and independent probability distributions, where $D_{i}\in\mathcal{D}$ is defined over a set of $z$ locations $q_{i,j}\in\mathbb{R}^{d}\cup\{\bot\}$. Note that the mapping $\varphi$ maps the locations $q_{i,j}$ from $\mathbb{R}^{d}$ to $\varphi(q_{i,j})$ in $\mathcal{H}$, and we assume $\varphi(\bot)=\bot$. Then the probabilistic SVDD problem is given by the following adaptation of Definition 3.1.

Definition 3.3**.**

Let $\mathcal{D}$ be a set of $n$ discrete distributions, where each distribution is defined over $z$ locations in $\mathbb{R}^{d}\cup\left\{\bot\right\}$. Let $K\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a kernel function with associated feature map $\varphi\colon\mathbb{R}^{d}\rightarrow\mathcal{H}$. The probabilistic support vector data description (pSVDD) problem is to find a center $c^{*}\in\mathcal{H}$ that minimizes the expected SVDD cost, i.e.,

[TABLE]

where the expectation is taken over the randomness of $X\sim\mathcal{D}$.

Note that the deterministic problem, see Equation (1), is often stated with squared distances [28, 30]. It does not matter whether we minimize the maximum distance or any of its powers or any other monotone transformation. In the probabilistic case this is not true. Consider for instance squared distances. The resulting problem would be similar to a $1$-means rather than a $1$-median problem. Huang et al.[16] have observed that minimizing the expected maximum squared distance corresponds to minimizing the expected area of an enclosing ball in $\mathbb{R}^{2}$. This observation can be generalized to the expected volume of an enclosing ball in $\mathbb{R}^{p}$ when the $p$-th powers of distances are considered. Considering $p=2$ might also have advantages when dealing with Gaussian input distributions due to their strong connection to squared Euclidean distances. In a general setting of the probabilistic smallest enclosing ball problem however, it is natural to minimize in expectation the maximum Euclidean distance, as in Definitions 3.1 and 3.3, since its radius is the primal variable to minimize.

Next we want to show how to find a $(1+\varepsilon)$-approximation for the pSVDD problem. Explicitly computing any center $c\in\mathcal{H}$ takes $\Omega(D)$ time and space which is prohibitive not only when $D=\infty$. Note that for the SEB and SVDD problems, any reasonable center lies in the convex hull of the input points. Since taking the expectation is simply another linear combination over such centers, we can express any center $c\in\mathcal{H}$ as a linear combination of the set of non-empty locations, i.e.,

[TABLE]

The idea is to exploit this characterization to simulate Algorithm 1 and thereby Algorithm 2 to work in the feature space $\mathcal{H}$ by computing the centers and distances only implicitly.

For now, assume that any distance computation can be determined. Note that sampling a set $P_{i}\subset\mathbb{R}^{d}$ is the same as sampling the set $\varphi(P_{i})$ of corresponding points in $\mathcal{H}$ from the same distribution. We assume that we have a set of locations or realizations $\mathcal{P}=\{P_{1},\ldots,P_{N}\}$, with $P_{i}\subset\mathbb{R}^{d}$. The remaining steps are passed to Theorem 2.11 which is based on Algorithm 1. First, we show the invariant that each center $c_{i}$ reached during its calls to Algorithm 1 can be updated such that we maintain a linear combination $c_{i}=\sum_{u,v}\gamma_{u,v}\varphi(q_{u,v})$, where at most $i+1$ terms have $\gamma_{u,v}\neq 0$.

• •

The initial center $c_{0}\in\mathcal{H}$ is chosen by sampling uniformly at random a set $P\in\mathcal{P}$ via Lemma 2.8. We take any point $q\in P$, $q\neq\bot$, which maps to $c_{0}=\varphi(q)$. Thus the invariant is satisfied at the beginning, where the corresponding coefficient is $\gamma=1$ and all other coefficients are zero.

• •

In each iteration we randomly sample a set $P\in\mathcal{P}$ to simulate the approximate subgradient $\tilde{g}\left(c_{i}\right)$ at the current point $c_{i}$. The vector $\tilde{g}\left(c_{i}\right)$ is a vector between $c_{i}$ and some point $p_{j,k}=\varphi(q_{j,k})\in\mathcal{H}$, such that $q_{j,k}$ maximizes $\left\|c_{i}-\varphi(q^{\prime})\right\|$ over all $q^{\prime}\in P$ (cf. Lemma 2.6).

To implicitly update to the next center $c_{i+1}$ note that (cf. Algorithm 1)

[TABLE]

Assume the invariant was valid that $c_{i}$ was represented as $c_{i}=\sum_{u,v}\gamma_{u,v}\varphi(q_{u,v})$ with at most $i+1$ non-zero coefficients. Then it also holds for the point $c_{i+1}$ since the previous non-zero coefficients of $\varphi(q_{u,v})$ are multiplied by $1-s/\left\|c_{i}-\varphi(q_{j,k})\right\|$ and the newly added $\varphi(q_{j,k})$ is assigned the coefficient $s/\left\|c_{i}-\varphi(q_{j,k})\right\|$. So there are at most $i+2$ non-zero coefficients.

Therefore, we do not have to store the points $c_{i}$ explicitly while performing Algorithm 1. The implicit representation can be maintained via a list storing points that appear in the approximate subgradients and their corresponding non-zero coefficients.

To actually compute the coefficients, we need to be able to compute Euclidean distances as well as determine $s$. Using Lemma 2.9 we determined the step size $s$ via an estimator $\tilde{R}=\sum_{P\in\mathcal{S}}m(c_{0},P)$ based on a small sample $\mathcal{S}$. In particular this requires distance computations again. To this end, we show how to compute $\left\|c_{i}-\varphi(q)\right\|$ for any location $q\in\mathbb{R}^{d}$. Recall that the kernel function implicitly defines the inner product in $\mathcal{H}$. It thus holds that

[TABLE]

where $w,w^{\prime}\in\{0,\ldots,i\}$ index the locations $q_{w},q_{w^{\prime}}$ with corresponding $\gamma_{w},\gamma_{w^{\prime}}\neq 0$ in iteration $i$. Therefore, we have the following Theorem.

Theorem 3.4**.**

Let $\mathcal{D}$ be a set of $n$ discrete distributions, where each distribution is defined over $z$ locations in $\mathbb{R}^{d}\cup\{\bot\}$. There exists an algorithm that implicitly computes $\tilde{c}\in\mathcal{H}$ that with constant probability is a $(1+\varepsilon)$-approximation for the probabilistic support vector data description problem. I.e., it holds that

[TABLE]

where the expectation is taken over the randomness of $X\sim\mathcal{D}$. The running time of the algorithm is ${O(dn\cdot\left(z/\varepsilon^{3}\cdot\log 1/\varepsilon+1/\varepsilon^{8}\cdot\log^{2}1/\varepsilon\right))}$.

4 Conclusion and open problems

We studied the set median problem in high dimensions that minimizes the sum of maximum distances to the furthest point in each input set. We presented a $(1+\varepsilon)$-approximation algorithm whose running time is linear in $d$ and independent of the number of input sets. We further discussed that in high dimensions the size of the input sets cannot be reduced sublinearly without losing a factor of roughly $\sqrt{2}$. Our work resolves an open problem of [22] and improves the previously best algorithm for the probabilistic smallest enclosing ball problem in high dimensions by reducing the dependence on $d$ from exponential to linear. This enables running the algorithm in high dimensional Hilbert spaces induced by kernel functions which makes it more flexible and viable as a building block in machine learning and data analysis. As an example we transferred the kernel based SVDD problem of [28] to the probabilistic data setting. Our algorithms assume discrete input distributions. It would be interesting to extend them to various continuous distributions. The pSEB problem minimizes the expected maximum distance. When it comes to minimizing volumes of balls or in the context of Gaussian distributions it might be interesting to study higher moments of this variable. This corresponds to a generalization of the set median problem to minimizing the sum of higher powers of maximum distances. Finally we hope that our methods may help to extend more shape fitting and machine learning problems to the probabilistic setting.

Appendix A Omitted proofs

Proof A.1** (Proof of Lemma 2.2).**

The non-negativity and symmetry properties of $m$ follow from the corresponding metric properties in Euclidean space $(\mathbb{R}^{d},\left\|\cdot\right\|)$ and by definition. If $A=B$, then $m(A,B)=0$ holds by definition. Otherwise there exist elements $a\in A$, $b\in B$, $a\neq b$, and thus $m(A,B)\geq\left\|a-b\right\|>0$. This proves the identity of indiscernible elements.

To prove the validity of the triangle inequality, let $A,B,C\in\mathcal{X}$ be distinct. Let $a\in A$, $c\in C$ be points such that $m(A,C)=\left\|a-c\right\|$. For any $b\in B$ it holds that

[TABLE]

using triangle inequality in $(\mathbb{R}^{d},\left\|\cdot\right\|)$ and the definition of $m$. Now consider the cases where at least two sets are equal. In the case that $A=C$ the claim follows from the non-negativity property. In the case that $A=B$, we have

[TABLE]

The case $B=C$ is analogous.

Proof A.2** (Proof of Lemma 2.3).**

Let $x,y\in\mathbb{R}^{d}$, let $\lambda_{1}\in[0,1],\lambda_{2}=1-\lambda_{1}$, and let $p_{i}^{*}$ maximize $\left\|\lambda_{1}x+\lambda_{2}y-p_{i}\right\|$ over all $p_{i}\in P_{i}$. Then

[TABLE]

follows as we have claimed.

Proof A.3** (Proof of Lemma 2.4).**

Fix any $x,y\in\mathbb{R}^{d}$. Let $p_{i}^{*}\in\operatorname{argmax}\nolimits_{p_{i}\in P_{i}}\left\|x-p_{i}\right\|$. By definition of $f$ and applying the triangle inequality to every single term we have

[TABLE]

Proof A.4** (Proof of Lemma 2.5).**

Let $c^{*}\in\operatorname{argmin}_{c\in\mathbb{R}^{d}}f(c)$. We first prove that for each term $j\in[N]$

[TABLE]

Suppose $c_{i}=p_{j}$, then $\left\langle 0,c_{i}-c^{*}\right\rangle=0\geq 0-m(c^{*},P_{j})=m(c_{i},P_{j})-m(c^{*},P_{j})$. Otherwise let $p_{j}^{*}\in P_{j}$ be a point such that $\left\|c^{*}-p_{j}^{*}\right\|=m(c^{*},P_{j})$. We have

[TABLE]

which follows by Cauchy-Schwarz inequality (CSI) and maximality of $p_{j}^{*}$. Now summing Equation (3) over all $j\in[N]$ we have

[TABLE]

Proof A.5** (Proof of Lemma 2.6).**

The vector $\tilde{g}(c_{i})$ is normalized by definition except if an index $j$ is chosen such that $c_{i}=p_{j}$, in which case $\left\|\tilde{g}(c_{i})\right\|=0$. Thus $\mathbb{E}\left[\left\|\tilde{g}(c_{i})\right\|^{2}\right]=\mathbb{E}\left[\left\|\tilde{g}(c_{i})\right\|\right]\leq 1$ holds. Also we have

[TABLE]

cf. Lemma 2.5.

Proof A.6** (Proof of Theorem 2.7).**

Assume that we have reached a center $c_{i}\in\mathbb{R}^{d}$ while running Algorithm 1. Recall from Lemma 2.6 that $\mathbb{E}_{\tilde{g}}\left[\left\|\tilde{g}(c_{i})\right\|^{2}\mid c_{i}\right]\leq 1$ and $\mathbb{E}_{\tilde{g}}\left[{\tilde{g}(c_{i})}\mid c_{i}\right]=g(c_{i})/N$. We have

[TABLE]

The law of total expectation implies, by taking expectations over $c_{i}$ on both sides and rearranging, that

[TABLE]

We sum the equations for $i\in\left\{0,\ldots,\ell\right\}$ such that the terms $\mathbb{E}_{\tilde{g}}\left[\left\|c_{k}-c^{*}\right\|^{2}\right]$ for $k\in\left\{1,\ldots,\ell\right\}$ on both sides cancel and continue this derivation. We have

[TABLE]

using the subgradient property and the fact that for any set of positive real valued random variables $X_{i}$ we have $\forall i\colon\min_{j}(X_{j})\leq X_{i}$. This implies $\forall i\colon\mathbb{E}\left[\min_{j}(X_{j})\right]\leq\mathbb{E}\left[X_{i}\right]$ and thus $\mathbb{E}\left[\min_{j}(X_{j})\right]\leq\min_{j}\mathbb{E}\left[X_{j}\right]$. Rearranging and substituting $R=\left\|c_{0}-c^{*}\right\|$ now implies

[TABLE]

Proof A.7** (Proof of Lemma 2.8).**

Define the random variable $X=m(c^{*},P)$. Clearly we have $X\geq 0$. Its expectation equals $\mathbb{E}\left[X\right]=\sum\nolimits_{i=1}^{N}\mathbf{Pr}\left[P=P_{i}\right]\cdot m(c^{*},P_{i})=\sum\nolimits_{i=1}^{N}m(c^{*},P_{i})/N=f(c^{*})/N$. Thus, by Markov’s inequality we have $\mathbf{Pr}\left[X>f(c^{*})/(\delta_{1}N)\right]\leq\delta_{1}$. Now choose an arbitrary $c_{0}\in P$. We have $R=\left\|c_{0}-c^{*}\right\|\leq m(c^{*},P)\leq f(c^{*})/(\delta_{1}N)$ with probability at least $1-\delta_{1}$.

Proof A.8** (Proof of Lemma 2.9).**

Let $\Delta=f(c^{*})/N$. Let $c_{0}$ be the initial center chosen as described in Lemma 2.8 with absolute constant $\delta_{1}=1/8\geq 1/81>\varepsilon^{2}$. Consider the two balls $B_{1}(c_{0},\varepsilon\Delta)$ and $B_{2}(c^{*},\left(1/\varepsilon^{2}\right)\Delta)$. Then $\left\|c_{0}-c^{*}\right\|\leq 8\Delta$ holds with constant probability $1-\delta_{1}$ and clearly $c_{0}\in B_{2}$.

Let $\mathcal{Q}$ consist of all sets of $\mathcal{P}$ that are fully contained in $B_{2}$. We have $|\mathcal{Q}|\geq(1-\varepsilon^{2})N$, since otherwise $f(c^{*})\geq\sum_{P\in\mathcal{P}\setminus\mathcal{Q}}m(c^{*},P)>\left(\varepsilon^{2}N\right)\cdot\left(\Delta/\varepsilon^{2}\right)=f(c^{*})$.

Now, sample a collection $\mathcal{S}$ of $1/\varepsilon$ sets, each uniformly from $\mathcal{P}$. All of our samples are completely contained in the ball $B_{2}$ with constant probability. By a union bound over the elements of $\mathcal{S}$, the probability that this fails is at most $\delta_{2}\leq\varepsilon^{2}\cdot 1/\varepsilon=\varepsilon<1/8$.

Now, let $\tilde{R}=\sum_{P\in\mathcal{S}}m(c_{0},P)$ be our estimate. Clearly we can compute it in time ${O(dn/\varepsilon)}$. We need to show that $\tilde{R}$ is close to the average cost $\Delta$ as we have claimed. To this end, consider the following two cases:

At most $(1-2\varepsilon)|\mathcal{Q}|$ sets of $\mathcal{Q}$ are completely contained in $B_{1}$

In this case we have a constant probability $1-\delta_{3}$, for $\delta_{3}<2\varepsilon<2/8$, that there is a set $Q\in\mathcal{Q}\cap\mathcal{S}$ that contains some point $q\in Q$ that lies outside $B_{1}$. Therefore we have $\tilde{R}\geq m(c_{0},Q)\geq\left\|c_{0}-q\right\|\geq\varepsilon\Delta$. For the upper bound, note that the diameter of $B_{2}$ is $2\Delta/\varepsilon^{2}$. Since all our samples are contained in that ball, we have $\tilde{R}\leq 2\Delta|\mathcal{S}|/\varepsilon^{2}=2\Delta/\varepsilon^{3}$. Thus claim a) holds.

At least $(1-2\varepsilon)|\mathcal{Q}|$ sets of $\mathcal{Q}$ are completely contained in $B_{1}$

Suppose the second item b) does not hold, i.e. we have $R=\left\|c_{0}-c^{*}\right\|>4\varepsilon\Delta$. We can bound the number of sets that are not fully contained in $B_{1}$ by $|\mathcal{P}\setminus\mathcal{Q}|+2\varepsilon|\mathcal{Q}|\leq\varepsilon^{2}N+2\varepsilon N\leq 3\varepsilon N$. Let $\mathcal{B}=\{P\in\mathcal{P}\mid P\subseteq B_{1}\}$ be the remaining family of sets in $\mathcal{P}$ that are fully contained in $B_{1}$. Clearly $|{\mathcal{B}}|\geq(1-3\varepsilon)N$. Now we compare the cost of using the center $c_{0}$ to the optimal cost. For every $P\in\mathcal{P}$ we have $|m(c^{*},P)-m(c_{0},P)|\leq\left\|c_{0}-c^{*}\right\|$ by the triangle inequality. For each $P\in\mathcal{B}$ it holds that $m(c^{*},P)-m(c_{0},P)\geq(\left\|c_{0}-c^{*}\right\|-\varepsilon\Delta)-m(c_{0},P)\geq\left\|c_{0}-c^{*}\right\|-2\varepsilon\Delta$. So using $\left\|c_{0}-c^{*}\right\|>4\varepsilon\Delta$ and $\varepsilon<1/9$ we can deduce

[TABLE]

which contradicts the optimality of $c^{*}$. Thus, claim b) holds in this case.

Proof A.9** (Proof of Theorem 2.10).**

Let $b$ be an arbitrary center in $\mathcal{C}$ with

[TABLE]

If there is no such center then all centers are good approximations, in which case the theorem is trivial. There are at most $|\mathcal{C}|$ choices for $b$. We study the random variable

[TABLE]

where $m(\hat{c},b)=\left\|\hat{c}-b\right\|$, and $h(Q)$ denote the summands of Equation (5). Since $m$ is a metric, by triangle inequality it holds that $X$ is the sum of random variables between [math] and $1$. The bad event is $X\leq|\mathcal{Q}|/2$.

If we denote by $\mathds{1}_{Q\in\mathcal{Q}}$ the indicator function that $Q\in\mathcal{Q}$, then we have $X=\sum_{Q\in\mathcal{Q}}h(Q)=\sum_{Q\in\mathcal{P}}h(Q)\cdot\mathds{1}_{Q\in\mathcal{Q}}$, and it holds that

[TABLE]

Equation (4) and the triangle inequality for any set $Q\in\mathcal{P}$: $m(b,Q)+m(\hat{c},Q)\geq m(\hat{c},b)$, imply that

[TABLE]

Thus it holds that

[TABLE]

It follows that

[TABLE]

and the bad event is bounded by

[TABLE]

where $\gamma=\varepsilon/\left(2+{2}\varepsilon\right)$, and $\gamma\geq\varepsilon/4$ holds since $\varepsilon<1$. Using a Chernoff bound [20] we have that

[TABLE]

since $\gamma\geq\varepsilon/4$ and $\gamma^{2}/(1-\gamma)\geq\varepsilon^{2}/16$. Transforming $X$ in Equation (6) we have

[TABLE]

Defining the set of bad centers $\mathcal{B}=\{b\in\mathcal{C}\mid\sum_{P\in\mathcal{P}}m(b,P)>(1+\varepsilon)\sum_{P\in\mathcal{P}}m(\hat{c},P)\}$, we finally have that

[TABLE]

Then it holds in particular for $a\in\mathcal{C}$ that with probability at least $1-|\mathcal{C}|\exp(-\varepsilon^{2}|\mathcal{Q}|/64)$

[TABLE]

But it holds that $\sum_{Q\in\mathcal{Q}}m(a,Q)\leq\sum_{Q\in\mathcal{Q}}m(\hat{c},Q)$ by optimality of $a$, so the contrapositive yields that

[TABLE]

Proof A.10** (Proof of Theorem 2.11).**

Set $\ell=(68/\varepsilon)^{2}$. Using Lemma 2.9, either our initial center $c_{0}$ satisfies $\left\|c_{0}-c^{*}\right\|\leq 4\varepsilon f(c^{*})/N$. In that case Lemma 2.4 yields $f(c_{0})\leq f(c^{*})+N\left\|c_{0}-c^{*}\right\|\leq(1+4\varepsilon)f(c^{*})$, so the starting point is already a good center.

Otherwise we have $\varepsilon f(c^{*})/N\leq\tilde{R}\leq\left(2/\varepsilon^{3}\right)\cdot f(c^{*})/N$ and $R=\left\|c_{0}-c^{*}\right\|\leq 8f(c^{*})/N$ by Lemma 2.9. Thus $\varepsilon^{3}\tilde{R}/2\leq f(c^{*})/N\leq\tilde{R}/\varepsilon$. To improve this we run the main loop of Algorithm 2.1 for the step sizes $s=\tilde{R}_{j}/\sqrt{\ell+1}$, where $\tilde{R}_{j}=2^{j-1}\cdot\varepsilon^{3}\tilde{R}$ for all values of $0\leq j\leq\lceil\log(2/\varepsilon^{4})\rceil$. For some particular value of $j$ we have a $2$-approximation given by $f(c^{*})/N\leq\tilde{R}_{j}\leq 2f(c^{*})/N$. In this particular run, setting the step size $s=\tilde{R}_{j}/\sqrt{\ell+1}$ and plugging this into the bound given in Theorem 2.7 we have that

[TABLE]

Using Markov’s inequality we have that

[TABLE]

The best center collected in all repetitions cannot be worse than this particular $\tilde{c}$ or $c_{0}$, see the cases of Lemma 2.9.

Finally we have a collection $\mathcal{C}$ of $|\mathcal{C}|\in O(1/\varepsilon^{2}\cdot\log 1/\varepsilon)$ centers and want to find one of them that is a $(1+\varepsilon)$ approximation for the best center in $\mathcal{C}$ via Theorem 2.10. We sample a collection of $64/\varepsilon^{2}\cdot\ln(8|\mathcal{C}|)\in O(1/\varepsilon^{2}\cdot\log 1/\varepsilon)$ point sets from $\mathcal{P}$ and find the best center for this subset of points which is the final output of our algorithm. By Theorem 2.10 this center is within another factor of $(1+\varepsilon)$ to the best in $\mathcal{C}$ with failure probability at most $\delta_{5}\leq 1/8$. The total approximation factor is thus at most $(1+4\varepsilon)(1+\varepsilon)\leq 1+9\varepsilon$. Rescaling $\varepsilon$ yields the correctness. The total failure probability is at most $\delta=\sum_{i=1}^{5}\delta_{i}\leq 6/8$ by a union bound over all bad events in Lemma 2.9 and this theorem.

We continue with the running time. The initial center $c_{0}$ and the estimate $\tilde{R}$ can be computed in $O(d{n}/\varepsilon)$ time, see Lemma 2.9. The main loop of Algorithm 1 takes $O(d{n})$ in each iteration and runs for $\ell\in O(1/\varepsilon^{2})$ iterations for a fixed step size. But we try $O(\log 1/\varepsilon)$ different step sizes. This makes up a running time of $O(d{n}/\varepsilon^{2}\cdot\log 1/\varepsilon)$. Finally we evaluate the objective function for $O(1/\varepsilon^{2}\log 1/\varepsilon)$ centers for the sample of $O(1/\varepsilon^{2}\cdot\log 1/\varepsilon)$ sets taken via Theorem 2.10. This can be done in time $O(d{n}/\varepsilon^{4}\cdot\log^{2}1/\varepsilon)$ which dominates the running time as we have claimed.

Proof A.11** (Proof of Theorem 2.12).**

We reduce from the indexing problem which is known to have $\Omega\left(n\right)$ one-way randomized communication complexity [18]. Alice is given a vector $a\in\{0,1\}^{n}$. Bob has an index $i\in[n]$, and has to guess the value of the $i$-th bit of $a$, denoted $a_{i}$, with probability $2/3$.

It is known from [2] that there is a centrally symmetric point set $K$ of size $\Omega\left(\exp\left(d^{1/3}\right)\right)$ on the unit hypersphere in $\mathbb{R}^{d}$ centered at the origin, such that for any pair of distinct points $p,q\in K$ it holds that

[TABLE]

unless $p\neq-q$, in which case $\left\|p-q\right\|=2$.

Let $d$ be the smallest integer such that $d\geq 8$ and $n\leq\exp\left(d^{1/3}\right)$. We choose a set of $\exp\left(d^{1/3}\right)$ pairs of centrally symmetric points of $K$. We may assume that there is a lexicographic order of these pairs, so there is a mapping between the indices of $a$ and the pairs of points of $K$ that is known to both, Alice and Bob. Alice constructs the set $S$ by including the first point of the $i$-th pair, denoted $p_{i}$, if and only if $a_{i}=1$. She builds a data structure $\Sigma_{S}$ which she sends to Bob.

Let the data structure be such that for any $x\in\mathbb{R}^{d}$ the answer to a query $\Sigma_{S}(x)$ satisfies

[TABLE]

for some constant $1<\alpha<\sqrt{2}\left(1-2/d^{1/3}\right)$. We consider two cases:

• •

If $a_{i}=1$ then $p_{i}$ is included in $S$. Thus $\Sigma_{S}\left(-p_{i}\right)\geq m\left(-p_{i},S\right)/\alpha=2/\alpha$. Since $\alpha<\sqrt{2}\left(1-2/d^{1/3}\right)$ it holds that $\Sigma_{S}\left(-p_{i}\right)>\sqrt{2}\left(1+2/d^{1/3}\right)$.

• •

If $a_{i}=0$ then $p_{i}\notin S$ and thus $\Sigma_{S}\left(-p_{i}\right)\leq m\left(-p_{i},S\right)\leq\sqrt{2}\left(1+2/d^{1/3}\right).$

Thus if $\alpha<\sqrt{2}\left(1-2/d^{1/3}\right)$ Bob could based on $\Sigma_{S}$ solve the indexing problem by querying $q_{i}=-p_{i}$. Consequently any encoding of $\Sigma_{S}$ uses $\Omega\left(\min\{n,\exp\left(d^{1/3}\right)\}\right)$ bits of space.

Proof A.12** (Proof of Theorem 3.2).**

The correctness of the algorithm follows from [22] and Theorem 2.11. It remains to analyze the running time. In the first case we go through all input distributions and use $k$ independent copies of a weighted reservoir sampler [8, 12] to get the $k$ samples. This takes $O(dnzs)\subseteq O(dnz/\varepsilon^{2}\cdot\log 1/\varepsilon)$ time. The subsampled problem is then solved via Theorem 2.11 with ${n}=1$ in time $O(d/\varepsilon^{4}\cdot\log^{2}1/\varepsilon)$ with failure probability at most $\delta=\sum_{i=1}^{5}\delta_{i}\leq 6/8$.

In the second case, each realization can be sampled similarly in time $O(dnz)$ but the probability that a realization is non-empty can only be lower bounded by $\varepsilon$. However this means that the expected number of samples that we need to take in order to have $k$ non-empty realizations is bounded by at most $k/\varepsilon$. Thus, by an application of Markov’s inequality the probability that we need more than $8k/\varepsilon\in O(1/\varepsilon^{3}\cdot\log 1/\varepsilon)$ trials is bounded by at most $\delta_{6}=1/8$. So we can assume with constant probability that this step succeeds in time $O(dnz/\varepsilon^{3}\cdot\log 1/\varepsilon)$. The subsampled problem is then solved via Theorem 2.11 in time $O(dn/\varepsilon^{4}\cdot\log^{2}1/\varepsilon)$. The final failure probability is at most $\delta=\sum_{i=1}^{6}\delta_{i}\leq 7/8$.

Proof A.13** (Proof of Theorem 3.4).**

With the described adaptations the correctness of the algorithm follows from Theorem 3.2.

The running time increases by a factor that is imposed by the simulation of the distance computations within Algorithm 1. Note that by the invariant there are at most $i+1$ non-zero coefficients in the $i$-th step. Expression (3.2) can thus be evaluated in time $O(i^{2}d)$ assuming $K$ can be evaluated in time $O(d)$. We conclude:

• •

The sampling part of Algorithm 2 does not change and thus runs in time ${O(dnz/\varepsilon^{3}\log 1/\varepsilon)}$.

• •

Estimating $\tilde{R}=\sum_{P\in\mathcal{S}}m(c_{0},P)$ takes time $O(dn/\varepsilon)$ since $i=0$, $|\mathcal{S}|=O(1/\varepsilon)$, and for each $P$ we have $|P|\leq n$.

• •

The subgradient computation in the $i$-th iteration of the main loop takes $O(dni^{2})$ time since it needs to maximize over $n$ distances. This means that for $\ell\in{O(1/\varepsilon^{2})}$ iterations we need ${O(dn\sum_{i=1}^{\ell}i^{2})}={O(dn\ell^{3})}={O(dn/\varepsilon^{6})}$ time. This is repeated ${O(\log 1/\varepsilon)}$ times, which implies a running time of ${O(dn/\varepsilon^{6}\cdot\log 1/\varepsilon)}$.

• •

The evaluation of the minimum at the end reaches as before over ${O(1/\varepsilon^{2}\cdot\log 1/\varepsilon)}$ centers, each defined by $O(\ell)$ non-zero coefficients. Each is evaluated with respect to a sample of ${O(1/\varepsilon^{2}\cdot\log 1/\varepsilon)}$ sets from the input. Since each maximum distance evaluation takes ${O(dn\ell^{2})}={O(dn/\varepsilon^{4})}$ time, we have that the total time for evaluation is ${O(dn/\varepsilon^{8}\cdot\log^{2}1/\varepsilon)}$.

Summing the running times yields the claim.