Near Linear Time Approximation Schemes for Clustering of Partially Doubling Metrics

A simple idea to solve this special case is to project all points from $X$ to the plane $Y$ and solve the resulting low-dimensional problem. Such an approach will result in a constant approximation as argued in the following. For the analysis we can think of the projection as moving all points to $Y$. By the triangle inequality, this will change the cost of any solution by at most the sum of distances the points have been moved. Since this sum is also a lower bound on the connection cost of any solution to the high-dimensional input, we get that any solution on the projection is a constant approximation. For similar reasons approximating the distances from $x$ to $y\in Y$ by $𝐝(x,{\pi }_Y(x))+𝐝({\pi }_Y(x),y)$, where ${\pi }_Y(x)$ denotes the point from $Y$ closest to $x$, only gives a constant approximation. Instead of replacing $x$ by a single point, our new idea is to replace $x$ by a set of points $N_x$ (which we call proxies), which is an $\epsilon 𝐝(x,{\pi }_Y(x))$-net of a ball $B_Y({\pi }_Y(x),𝐝(x,{\pi }_Y(x))/\epsilon )\subseteq Y$ around ${\pi }_Y(x)$ of radius $𝐝(x,{\pi }_Y(x))/\epsilon$. (Formally, a $\rho$-net of a set $B$ is a subset $N\subseteq B$, such that 1) the interpoint distances in $N$ are at least $\rho$, and 2) every point in $B$ has a nearby point in $N$ within distance $\rho$.) Taking ${\min }_{u\in N_x}𝐝(x,u)+𝐝(u,y)$ for $y\in Y$ we approximate all distances within a factor of $(1+\epsilon )$ as summarized in the following lemma.

For technical reasons, our algorithm does not directly compute the proxy set $N_x$ around the exact nearest neighbor ${\pi }_Y(x)$ of $x$ in $Y$. Instead, we will first compute an $O(1)$-approximate solution $S\subseteq Y$, find in $S$ an $O(1)$-approximate nearest neighbor $s$ of $x$, and construct $N_x$ as an $\epsilon 𝐝(x,s)$-net of the ball $B_Y(s,𝐝(x,s)/\epsilon )$. This replacement introduces an extra $\epsilon 𝐝(x,S)$ additive error in Lemma 6, which will not affect the correctness of our analysis. We note that conceptually similar ideas also appear in [16].

We combine the above dimension reduction idea with the algorithm of [16]. We start from a standard hierarchical decomposition for doubling metrics by Talwar [41], which is an analogue to the randomly-shifted quadtree in Euclidean space. The decomposition $\mathscr{H}$ is constructed on top of $Y$, and has $L=O(\log m)$ levels. For each $0\le \mathrm{\ell }\le L$, level $\mathrm{\ell }$ of $\mathscr{H}$, denoted by ${\mathscr{H}}_{\mathrm{\ell }}$, is a partition of $Y$ into clusters of diameter at most $2^{\mathrm{\ell }}$, and ${\mathscr{H}}_{\mathrm{\ell }-1}$ is a refinement of ${\mathscr{H}}_{\mathrm{\ell }}$. $\mathscr{H}$ can be represented as a tree, where each node corresponds to a cluster and has $2^{O(\mathrm{ddim})}$ children. In [16], a dynamic program is run on $\mathscr{H}$ to compute a $(1+\epsilon )$-approximation for facility location. In our setting, one immediate issue is that clients in $X$ are not directly defined on $\mathscr{H}$, but are replaced by a proxy set $N_x\subseteq Y$. Therefore, it is not immediately clear that one can follow a similar approach as [16].

For the sake of presentation, let us first have a brief review of the argument of [16], for the case $X=Y$. They use the notion of portals, which originates from [2]. Roughly speaking, the portal set $P_C$ for a cluster $C\in {\mathscr{H}}_{\mathrm{\ell }}$ is an $\epsilon 2^{\mathrm{\ell }}$-net of $C$. The level $\mathrm{\ell }$ portals are the union of $P_C$ over $C\in {\mathscr{H}}_{\mathrm{\ell }}$. For a pair of points $x,y\in Y$, the actual distance between $x$ and $y$ will be replaced by the portal-respecting distance ${𝐝}_{\mathrm{port}}^{\mathscr{H}}(x,y)$, which is the length of the portal-respecting path between $x$ and $y$. Specifically, let $\mathrm{\ell }$ be the highest level where $\{x,y\}$ is cut w.r.t. $\mathscr{H}$, i.e., the highest level where $x$ and $y$ fall into different clusters. The portal-respecting path between $x$ and $y$ starts from $x$ (a portal at level $0$), each step connecting the current portal to the closest portal at one level up, until reaching a portal $p$ at level $\mathrm{\ell }$. Portal $p$ is then connected to another portal $q$ at the same level, and the path goes all the way down from $q$ to $y$. It can be shown that the portal-respecting distance ${𝐝}_{\mathrm{port}}^{\mathscr{H}}(x,y)$ is upper bounded by $𝐝(x,y)+O(\epsilon )2^{\mathrm{\ell }}$. A number of techniques are developed in [16] to bound $\mathrm{\ell }$, the level where $\{x,y\}$ is cut.

Finally, a dynamic program is run on $\mathscr{H}$ w.r.t. portal-respecting distance ${𝐝}_{\mathrm{port}}^{\mathscr{H}}$. Each entry of the DP table is encoded by a cluster $C$ of $\mathscr{H}$, and a configuration ${𝐚}_C$ indicating how $C$ interacts with other clusters as well as the current solution (facility set) $F$ via its portals $P_C$. Roughly speaking, the configuration is a vector ${𝐚}_C\in {ℝ}^{P_C}$, whose $p$-th coordinate ${𝐚}_C(p)$ equals to the minimum portal-respecting distance between portal $p$ and facilities in $F$. The value stored in table entry $(C,{𝐚}_C)$, denoted by $f(C,{𝐚}_C)$, is defined as the minimum facility location cost inside of $C$, given the configuration ${𝐚}_C$, i.e.,

To compute the DP table $f(C,{𝐚}_C)$, the algorithm enumerates all children $D$ of $C$ and all possible configurations. If the configurations $\{{𝐚}_D:D\text{ is a child of }C\}$ are consistent with ${𝐚}_C$, and the sum of $f(D,{𝐚}_D)$ over all children $D$ is less than the current value of $f(C,{𝐚}_C)$, it then updates $f(C,{𝐚}_C)$ to be the sum of $f(D,{𝐚}_D)$.

In our setting, if we would like to combine our approach with the algorithm of [16], the first step is to combine our dimension reduction idea with the portal-respecting distance. This is done in a rather straightforward way, by modifying (1) to

Now, a main challenge is to bound the error incurred by the portal-respecting distance, namely, ${\widehat{𝐝}}_{\mathrm{port}}^{\mathscr{H}}(x,y)-𝐝(x,y)$. In [16], this can be done by bounding the highest cutting level for a single pair of points $\{x,y\}$. However, since $x$ is now represented by $N_x$, we have to bound the cutting level of $\{u,y\}$ for every $u\in N_x$. Nevertheless, we argue that the above error can still be effectively bounded, by only considering the cutting level for two sets $B_Y({\pi }_Y(x),𝐝(x,Y)/\epsilon )$ and $\{{\pi }_Y(x),y\}$. Specifically, if $j$ is the highest level where $B_Y({\pi }_Y(x),𝐝(x,Y)/\epsilon )$ is cut, and $\mathrm{\ell }$ is the highest level where $\{{\pi }_Y(x),y\}$ is cut, then we show that

We further integrate the dimension reduction with the dynamic program framework in [16]. The algorithm is run on $\mathscr{H}$. Whenever we compute the connection cost from a client $x\in X$ to some set of facilities $F\subseteq Y$, we replace $𝐝(x,F)$ with ${\widehat{𝐝}}_{\mathrm{port}}^{\mathscr{H}}(x,F)$. Two challenges come with such a replacement. First, recall that ${\widehat{𝐝}}_{\mathrm{port}}^{\mathscr{H}}(x,F)$ depends on ${𝐝}_{\mathrm{port}}^{\mathscr{H}}(u,F)$ for all proxies $u\in N_x$; therefore, to compute it we have to enumerate all $u\in N_x$, which can be done only when $N_x$ is entirely contained in some cluster $C$. This is different from [16], where the connection cost of $x$ can be trivially computed at the leaf node $\{x\}$. Second, even if we successfully define such a cluster $C$ which contains $N_x$, it is unclear how we access the distance ${𝐝}_{\mathrm{port}}^{\mathscr{H}}(u,F)$.

To resolve the first issue, we find for every $x\in X$ a suitable cluster $C(x)$ that entirely contains $N_x$ in the preprocessing stage, and “defer” the computation of ${\widehat{𝐝}}_{\mathrm{port}}^{\mathscr{H}}(x,F)$ to cluster $C(x)$. More concretely, for a client $x\in X$ and a cluster $C$ on $\mathscr{H}$, we say $x$ is revealed in $C$ if $B_Y({\pi }_Y(x),𝐝(x,Y)/\epsilon )\subseteq C$, and say $x$ is newly revealed in $C$, denoted $C(x)=C$, if $C$ is the lowest level cluster where $x$ is revealed. Then each DP table entry $(C,{𝐚}_C)$ is defined to be the minimum revealed facility location cost of $C$, given the configuration ${𝐚}_C$, i.e.,

Since all clients are revealed in the largest cluster (root node) $Y$, the root node stores the optimal facility location cost.

For the second issue, let us consider how the DP table is updated. The revealed cost of $C$ can be decomposed into two parts – (1) the total revealed cost of $C$’s child clusters, and (2) the newly revealed cost of $C$, i.e., ${\sum }_{x:C(x)=C}{\widehat{𝐝}}_{\mathrm{port}}^{\mathscr{H}}(x,F)$. The first part can be obtained from the corresponding table entries of $C$’s children. For the second part, the key observation is that if $x$ is newly revealed in $C$ ($C(x)=C$), then the proxy set of $x$ is entirely contained in the portal set of $C$, i.e., $N_x\subseteq P_C$. Therefore, for each $u\in N_x$, ${𝐝}_{\mathrm{port}}^{\mathscr{H}}(u,F)$ can be directly obtained from the configuration ${𝐚}_C$ which has the information of the connection of $P_C$, and thus ${\widehat{𝐝}}_{\mathrm{port}}^{\mathscr{H}}(x,F)$ can be effectively computed.

Our final algorithm has a running time of $2^{2^t}\tilde{O}(n+m)$ where $t=O(\mathrm{ddim}\log \frac{\mathrm{ddim}}{\epsilon })$. We remark that this is faster than the algorithm of [16], which has $t^2$ instead of $t$ in the (double) exponent. This improvement comes from the fact that one can improve the analysis of their hierarchical composition [16]. The fact, that such an improvement is possible had been observed in [10] referring to the paper [1]. We provide a self-contained proof of this fact in the full version [19, Lemma 2.5].

For the reason above, we extend Talwar’s hierarchical decomposition to the case that only a subset of the metric space is doubling. We start from constructing Talwar’s decomposition on top of $X$; denote the resulting decomposition by $\mathscr{H}$. It then remains to decide how to add the points in $Y$ to $\mathscr{H}$. Intuitively, we want that every point $y\in Y$ always lies in the same cluster as ${\pi }_X(y)$, the nearest neighbor of $y$ in $X$. More concretely, our plan is to assign to every $y\in Y$ a suitable level $0\le h(y)\le L$, and add $y$ to $\mathscr{H}$ in such a way that

1. 1.

   $\{y\}$ is a leaf node at level $h(y)$;

2. 2.

   $y$ and ${\pi }_X(y)$ are in the same cluster at levels higher than (or equal to) $h(y)+1$.

We call such cluster $\{y\}$ an ornament at level $h(y)$. After adding all $y\in Y$ to $\mathscr{H}$, we obtain a hierarchical decomposition for $X\cup Y$, denoted by $𝒫$.

It could be tricky to define for every ornament $\{y\}$ the level $h(y)$ to which it should be attached as a leaf. On the one hand, we want $h(y)$ to be sufficiently large, so that $y$ can be covered by some level $h(y)$ portal (i.e., $𝐝(y,P_C)\le \epsilon 2^{h(y)}$ for the cluster $C\ni y$). Therefore, we can extend the definition of portal-respecting path to $𝒫$, and define portal-respecting distance ${𝐝}_{\mathrm{port}}^{𝒫}(x,y)$ on $𝒫$ the same way as ${𝐝}_{\mathrm{port}}^{\mathscr{H}}(x,y)$. On the other hand, we want $h(y)$ to be sufficiently small, so that the error incurred by ${𝐝}_{\mathrm{port}}^{𝒫}(x,y)$, namely ${𝐝}_{\mathrm{port}}^{𝒫}(x,y)-𝐝(x,y)=O(\epsilon )2^{h(y)}$, is negligible.

As a first attempt, consider choosing $h(y)=\log (𝐝(y,X)/\epsilon )$. Then we have $𝐝(y,{\pi }_X(y))\le \epsilon 2^{h(y)}$, and we can further argue that $y$ is within distance $\epsilon 2^{h(y)+1}$ to some level $h(y)+1$ portal. Therefore, the portal-respecting distance ${𝐝}_{\mathrm{port}}^{𝒫}(x,y)$ can be defined the same way as ${𝐝}_{\mathrm{port}}^{\mathscr{H}}(x,y)$. However, this choice of $h(y)$ becomes problematic in terms of the error ${𝐝}_{\mathrm{port}}^{𝒫}(x,y)-𝐝(x,y)$. Consider an arbitrary point $x\in X$, and assume $h(y)$ is the highest level where $\{x,y\}$ is cut. Then the error is ${𝐝}_{\mathrm{port}}^{𝒫}(x,y)-𝐝(x,y)=O(\epsilon )2^{h(y)}=O(𝐝(y,X))$, which is too large to afford.

To resolve the issue, we choose $h(y)$ as a slightly smaller value $\log (𝐝(y,X)/\sqrt{\epsilon })$. For this choice, the error becomes ${𝐝}_{\mathrm{port}}^{𝒫}(x,y)-𝐝(x,y)=O(\epsilon )2^{h(y)}=O(\sqrt{\epsilon })𝐝(y,X)$, which is at most $O(\sqrt{\epsilon })𝐝(x,y)$ and thus can be charged to $𝐝(x,y)$. The tradeoff is that $y$ is no longer $\epsilon 2^{h(y)+1}$-covered by the level $h(y)+1$ portal set. Instead, the covering radius becomes $\sqrt{\epsilon }{2}^{h(y)+1}$. Nonetheless, we can still define ${𝐝}_{\mathrm{port}}^{𝒫}(x,y)$ similarly under this weaker covering property, and it does not change the $O(\sqrt{\epsilon })𝐝(x,y)$ error bound above when the cutting level is exactly $h(y)$. We obtain a weaker error bound of ${𝐝}_{\mathrm{port}}^{𝒫}(x,y)-𝐝(x,y)=O(\sqrt{\epsilon })2^{\mathrm{\ell }}$ only if the cutting level $\mathrm{\ell }$ is strictly greater than $h(y)$. In this case, the question reduces to finding the highest cutting level of $\{x,{\pi }_X(y)\}$, which can be answered fairly well using techniques in [16] and the previous section, since both points are in $X$.

Our new hierarchical decomposition $𝒫$ and the portal-respecting distance ${𝐝}_{\mathrm{port}}^{𝒫}(x,y)$ can then be combined with the dynamic program of Cohen-Addad et al. [16] to get our result in the case when the clients are low-dimensional. This seems difficult at first, because each cluster on $𝒫$ now has an unbounded number of child clusters, mainly due to the newly added ornaments. A naive enumeration of the configurations of these ornament children would blow up the time complexity of the DP.

Perhaps surprisingly, we show that we can avoid doing the enumeration for ornaments, and it suffices to only enumerate the configurations of non-ornament children, the number of which is bounded by $2^{O(\mathrm{ddim})}$. Our key observation is that ornaments must be candidate facilities, and thus only the opening cost needs to be computed. Therefore, we can obtain from portals of non-ornament children the information which ornaments are potentially required to be opened as a facility. Once we have this information, we simply select the smallest set of ornaments that serve all the unserved portals. This can be done via solving a set cover problem with a bounded universe.

For the case of the ambient space being one-dimensional we develop a new complexity reduction for the discrete Fréchet distance that can be summarized as follows. This provides an alternative approach that allows to apply the work of Cohen-Addad et al. [16] in a more direct way for this special case.

We believe that this complexity reduction is of independent interest. For example, it allows us to get an improved coreset construction for clustering under the discrete Fréchet distance.

Our idea for the complexity reduction can be described as follows. As a first (simple) step, we show that one can reduce the number of distinct values appearing in a time series of complexity $z$ to $O(\mathrm{\ell }/\epsilon )$ while maintaining the distance to any time series of complexity $\mathrm{\ell }$ up to a factor of $(1+\epsilon )$.

Then we observe that the discrete Fréchet distance between two time series $x$ and $y$, each of fixed complexity, can be written as a minimum over all traversals. Our goal is to describe the function minimizing over all traversals with a function minimizing over a much smaller set. During each traversal, every value $y_i$ of $y$ is matched to a subsequence of $x$. In order to determine the Fréchet distance, it suffices to know the minimum and maximum value of $x$ matched to $y_i$. As it turns out, each traversal can equivalently (but not uniquely!) be described by remembering a sequence of constraints that consist of the minimum and maximum value matched to each $y_i$. Furthermore, the function minimizing over the set of all possible traversals can likewise be described by minimizing over all possible ordered constraint sets. Since we have reduced the number of different values of the time series to $O(\mathrm{\ell }/\epsilon )$, the number of different constraint sets is small. The set of constraints will be called an $\mathrm{\ell }$-profile. It is important to note that the set of all $\mathrm{\ell }$-profiles of a time series $x$ with values from a fixed set $X$ completely determines the Fréchet distance to any time series of complexity $\mathrm{\ell }$. Since there are only a constant number of different sets of $\mathrm{\ell }$-profiles (where the constant depends on $\mathrm{\ell }$ and $\epsilon$) for any time series $x$ with $O(\mathrm{\ell }/\epsilon )$ distinct values, we can replace $x$ by the shortest time series $x^{\prime }$ over the same set of values that has the same set of $\mathrm{\ell }$-profiles. The length of the shortest such time series is a constant that depends on the set of profiles and the number of distinct values of the time series. Since the number of profiles is also a constant depending on $\epsilon$ and $\mathrm{\ell }$, the maximum length of these shortest time series is constant as well. We can compute such a time series using a dynamic programming approach. The resulting time series have length $O(2^{O{(\mathrm{\ell }/\epsilon )}^{(2\mathrm{\ell }+2)}})$.