# Pattern Matching in Doubling Spaces ††thanks:  This work was supported by
Basic Science Research Program through the National Research Foundation of
Korea (NRF) funded by the Ministry of Education (2017R1D1A1B04036529).

Corentin Allair
École Polytechnique, Paris, France
<EMAIL_ADDRESS>Antoine Vigneron
School of Electrical and Computer Engineering
UNIST, Republic of Korea
<EMAIL_ADDRESS>Corresponding author

###### Abstract

We consider the problem of matching a metric space $(X,d_{X})$ of size $k$
with a subspace of a metric space $(Y,d_{Y})$ of size $n\geqslant k$, assuming
that these two spaces have constant doubling dimension $\delta$. More
precisely, given an input parameter $\rho\geqslant 1$, the $\rho$-distortion
problem is to find a one-to-one mapping from $X$ to $Y$ that distorts
distances by a factor at most $\rho$. We first show by a reduction from
$k$-clique that, in doubling dimension $\log_{2}3$, this problem is NP-hard
and W[1]-hard. Then we provide a near-linear time approximation algorithm for
fixed $k$: Given an approximation ratio $0<\varepsilon\leqslant 1$, and a
positive instance of the $\rho$-distortion problem, our algorithm returns a
solution to the $(1+\varepsilon)\rho$-distortion problem in time
$(\rho/\varepsilon)^{O(1)}n\log n$. We also show how to extend these results
to the minimum distortion problem in doubling spaces: We prove the same
hardness results, and for fixed $k$, we give a $(1+\varepsilon)$-approximation
algorithm running in time
$(\operatorname{dist}(X,Y)/\varepsilon)^{O(1)}n^{2}\log n$, where
$\operatorname{dist}(X,Y)$ denotes the minimum distortion between $X$ and $Y$.

##  1 Introduction

A metric space has doubling dimension $\delta$ if any ball can be covered by
at most $2^{\delta}$ balls of half its radius. When $\delta=O(1)$, we say that
this space is doubling. (See Section 2.) For instance, the Euclidean space
$\mathbb{R}^{d}$ has doubling dimension $O(d)$, hence doubling spaces are
generalizations of fixed-dimensional Euclidean spaces.

In this paper, we study pattern matching problems in doubling spaces. Given
two doubling spaces $(X,d_{X})$ and $(Y,d_{Y})$ of doubling dimension
$\delta$, and sizes $|X|=k$ and $|Y|=n$, where $k\leqslant n$, our goal is to
find a subspace of $Y$ that resembles the pattern $X$. More precisely, we
consider the $\rho$-distortion problem and the minimum distortion problem,
which we describe below.

Given $\rho\geqslant 1$, the $\rho$-distortion problem is to find, if it
exists, a mapping $\sigma:X\to Y$ such that

$(1/\rho)d_{X}(x,x^{\prime})\leqslant
d_{Y}(\sigma(x),\sigma(x^{\prime}))\leqslant\rho d_{X}(x,x^{\prime})$ (1)

for all $x,x^{\prime}\in X$. It follows from this definition that $\sigma$ is
injective.

The $\rho$-distortion problem is analogous to the problem of matching two
point-sets in Euclidean space under rigid transformations, which are
compositions of translations and rotations. If, in addition, we allow scaling,
then an analogous problem in general metric spaces is the minimum distortion
problem. The goal is to minimize the distortion
$\operatorname{dist}(\sigma)=\operatorname{expansion}(\sigma)\times\operatorname{expansion}(\sigma^{-1})$
over all injections $\sigma:X\to Y$, where

$\mathrm{expansion}(\sigma)=\max_{\begin{subarray}{c}x,x^{\prime}\in X\\\
x\neq
x^{\prime}\end{subarray}}\frac{d_{Y}(\sigma(x),\sigma(x^{\prime}))}{d_{X}(x,x^{\prime})}\quad\text{and}\quad\mathrm{expansion}(\sigma^{-1})=\max_{\begin{subarray}{c}x,x^{\prime}\in
X\\\ x\neq
x^{\prime}\end{subarray}}\frac{d_{X}(x,x^{\prime})}{d_{Y}(\sigma(x),\sigma(x^{\prime}))}.$

The minimum of $\operatorname{dist}(\sigma)$ over all injections $\sigma:X\to
Y$ is denoted $\operatorname{dist}(X,Y)$, and it is easy to see that
$\operatorname{dist}(X,Y)\geqslant 1$. The minimum distortion problem was
introduced by Kenyon et al. [21] in the case where $k=n$, and thus $\sigma$ is
a bijection.

Motivated by applications to natural language processing, bioinformatics and
computer vision, Ding and Ye [13] recently proposed a practical algorithm for
a pattern matching problem in doubling spaces. However, this algorithm may
only return an approximation of a local minimum, and its time bound is not
given as a function of the input size. One of our goals is thus to provide a
provably efficient algorithm for pattern matching in doubling spaces. Another
motivation for our work is that, even though the complexity of the minimum
distortion problem has been studied for several types of metrics, it appears
that no result was previously known for two doubling metrics. (See the
comparison with previous work below.)

#### Our results.

We first give a hardness result: We show that for any $\rho\geqslant 1$, the
$k$-clique problem reduces to $\rho$-distortion in doubling dimension
$\log_{2}3$. It implies that the $\rho$-distortion problem is NP-hard, and is
W[1]-hard when parameterized by $k$ (Corollary 10). It also shows that this
problem cannot be solved in time $f(k)\cdot n^{o(k)}$ for any computable
function $f$, unless the exponential time hypothesis (ETH) is false (Corollary
11).

On the positive side, we present a near-linear time approximation algorithm
for small values of $k$. More precisely, if $\rho\geqslant 1$ and
$0<\varepsilon\leqslant 1$, our algorithm returns in $2^{O(k^{2}\log
k)}(\rho^{2}/\varepsilon)^{2k\delta}n\log n$ time a solution to the
$(1+\varepsilon)\rho$-distortion problem whenever a solution to the
$\rho$-distortion problem exists (Theorem 28). In this time bound, it is
reasonable to assume that $\rho$ is a small constant, say $\rho\leqslant 10$,
as a larger value would mean that we allow a relative error of more than 900%
in the quality of the matching, which is probably too much for most
applications.

We also show how to extend these results to the minimum distortion problem. In
particular, we show that the minimum distortion problem cannot be solved in
time $f(k)\cdot n^{o(k)}$ for any computable function $f$, unless ETH is
false, and we give a $(1+\varepsilon)$-approximation algorithm running in
$2^{O(k^{2}\log
k)}(\operatorname{dist}(X,Y)^{2k\delta}/\varepsilon^{2k\delta+O(1)})n^{2}\log
n$ time (Theorem 34). Here again, it is reasonable to assume that
$\operatorname{dist}(X,Y)=O(1)$, and then for any fixed $k$, this algorithm is
an FPTAS with running time $(1/\varepsilon)^{O(1)}n^{2}\log n$.

#### Comparison with previous work.

One of the main differences between our results and previous work on point
pattern matching under rigid transformations, or on the minimum distortion
problem, is that we parameterize the problem by $k$, and hence $k$ is regarded
as a small number. The advantage is that the dependency of our time bounds in
$n$ are low (near-linear or near-quadratic). However, we obtain an exponential
dependency in $k$, which may be unavoidable due to our hardness results.

Geometric point pattern matching problems have been studied extensively. (See
for instance the survey by Alt and Guibas [1].) In the fixed-dimensional
Euclidean space $\mathbb{R}^{d}$, these problems are usually tractable, as the
space of transformations has a constant number of degrees of freedom. For
instance, when $k=n$, we may want to decide whether $X$ and $Y$ are congruent,
which means that there is a rigid transformation $\mu$ such that $\mu(X)=Y$.
Alt et al. [2] showed how to find such a transformation in time $O(n^{d-2}\log
n)$, when it exists. In practice, however, we cannot expect that point
coordinates are known exactly, so it is unlikely that an exact match exists.
We may thus want to find the smallest $\varepsilon>0$ such that each point of
$X$ is brought to distance at most $\varepsilon$ from a point in $Y$. (In
other words, we allow an additive error $\varepsilon$.) Chew et al. gave an
$O(k^{3}n^{2}\log^{2}(nk))$-time algorithm to solve this problem in the plane
under rigid transformations [8].

As mentioned above, to the best of our knowledge, the only work published so
far on pattern matching in doubling spaces presents a practical algorithm for
matching two doubling spaces [13]. However it has not been proven to return a
good approximation of the optimal solution in the worst case. Other problems
studied in doubling spaces include approximate near-neighbor searching [3, 9,
20], spanners [4], routing [7, 16], TSP [26], clustering [15], Steiner forest
[6] …

The minimum distortion problem has been studied under various metrics, when
$k=n$. Kenyon et al. [21] gave a polynomial-time algorithm for line metrics
(1-dimensional point sets) when $\operatorname{dist}(X,Y)<5+2\sqrt{6}$. They
also gave an algorithm that computes $\operatorname{dist}(X,Y)$ when $d_{X}$
is the metric associated with an unweighted graph over $X$ and $d_{Y}$ is the
metric associated with a bounded degree tree over $Y$, with a running time
exponential in the maximum degree and doubly exponential in
$\operatorname{dist}(X,Y)$. For general metrics, the minimum distortion is
hard to approximate within a factor less than $\log^{1/4-\gamma}n$, for any
$\gamma>0$ [23]. Hall and Papadimitriou [18] showed that even for line
metrics, the distortion is hard to approximate when it is large.

When $k\leqslant n$, Fellows et al. [14] showed that the problem of deciding
whether $\operatorname{dist}(X,Y)\leqslant D$ is fixed-parameter tractable
when parameterized by $\Delta$ and $D$, where $d_{X}$ is the metric associated
with an unweighted graph over $X$, and $d_{Y}$ is the metric associated with a
tree of degree at most $\Delta$ over $Y$. For two unweighted graph metrics,
Cygan et al. [11] showed that the problem cannot be solved in time $2^{o(n\log
n)}$ unless ETH is false. When $X$ is an arbitrary finite metric space and $Y$
is a subset of the real line, Nayyeri and Raichel [25] showed that a constant-
factor approximation of $\operatorname{dist}(X,Y)$ can be computed in time
$\Phi(X)^{O((\operatorname{dist}(x,y))^{2})}(kn)^{O(1)}$, where $\Phi(X)$ is
the spread of $X$. (See Section 2.)

In summary, the previously known theoretical results on the minimum distortion
problem are either hardness results, or algorithms for cases where $d_{Y}$ is
a subset of a line metric or a tree metric. The algorithms presented in this
paper, on the other hand, apply when $d_{X}$ and $d_{Y}$ are doubling metrics,
which are generalizations of fixed-dimensional Euclidean metrics.

#### Our approach.

In Section 3, we present hardness results on the $\rho$-distortion problem. We
reduce an instance $G(V,E)$ of $k$-clique to an instance of $\rho$-distortion
consisting of two metric spaces $(X,d_{X})$ and $(Y,d_{Y})$ of sizes $k$ and
$km$, respectively, where $m=|V|$. The pattern $(X,d_{X})$ is an ultrametric,
with exponentially increasing distances. The space $(Y,d_{Y})$ consists of $k$
rings, each ring consisting of $m$ points regularly spaced on a circle of
perimeter 1. Each of these $m$ points is associated with a vertex of the
$k$-clique instance. (See Figure 1a.) The distances between the rings increase
exponentially (Figure 2), and the input graph is encoded by having slightly
longer edges for pairs of vertices lying in different rings that correspond to
edges in the input graph (Figure 1b.) We prove that these two spaces have
doubling dimension $\log_{2}3$, and that this instance of $\rho$-distortion is
equivalent to the $k$-clique instance we started from.

In Section 4, we give a self-contained description of a first approximation
algorithm for the $\rho$-distortion that runs in time $2^{O(k^{2}\log
k)}(\rho^{2}/\varepsilon)^{3k\delta}n+O(kn\log\Phi(Y))$ where $\Phi(Y)$ is the
spread of $Y$. (See Section 2.) We first construct a navigating net over $Y$
[24]. A navigating net is essentially a coordinate-free quadtree that records
a metric space. It represents $Y$ at all resolutions $r$ where $r$ is a powers
of 2. At each scale $r$, our navigating net records an $r$-net $Y_{r}$ of $Y$,
which is a subset of $Y$ whose points are a distance at least $r$ apart, and
such that the radius-$r$ balls centered at $Y_{r}$ cover $Y$. (See Figure 3.)

Let $r_{X}$ be the smallest scale that is at least $\rho$ times the diameter
of $X$. Our algorithm constructs, for each point $y\in Y_{r_{X}}$, a sparse
set of matchings whose images are in the radius-$3r_{X}$ ball centered at $y$.
(By sparse, we mean that any two such matchings send at least one point of $X$
to two points of $Y$ that are a distance at least $\varepsilon
r_{X}/(2\rho^{2})$ apart.) The union of these sets of matchings over all $y\in
Y_{r_{X}}$ is denoted $L(X,\varepsilon,r_{X})$, and we show that any solution
to the $\rho$-distortion problem is close to at least one matching in
$L(X,\varepsilon,r_{X})$.

We compute $L(X,\varepsilon,r_{X})$ recursively from $L(P,\beta,r_{P})$ and
$L(Q,\beta,r_{Q})$ for $\beta=\varepsilon/(8k-8)$, where $P$ and $Q$ form a
partition of $X$. More precisely, we obtain $P$ and $Q$ by running Kruskal’s
algorithm on $X$, and stopping at the second last step. It ensure that $P$ and
$Q$ are well separated, and it follows that any $\rho$-matching
$\bar{\sigma}:X\to Y$ can be approximated by a combination of two
$(1+\beta)\rho$-matchings $\bar{\sigma}^{P}_{\beta}:P\to Y$ and
$\bar{\sigma}^{Q}_{\beta}:Q\to Y$ recorded in $L(P,\beta,r_{P})$ and
$L(Q,\beta,r_{Q})$, respectively. (See Figure 4.) After computing
$L(X,\varepsilon,r_{X})$, we simply return one of the matchings that it
records, if any.

In Section 5, we show how to improve the time bound to $2^{O(k^{2}\log
k)}(\rho^{2}/\varepsilon)^{2k\delta}n\log n$. We achieve it using the
approximate near-neighbor (ANN) data structure by Cole and Gottlieb [9].
First, this data structure allows us to efficiently prune the sets of
matchings recorded at layer $r_{X}$, by only inserting a new matching if it is
far enough from all previously inserted matching. It can be checked by
performing a constant number of ANN queries in the set of matchings. As the
space of matchings is doubling (Corollary 3), it takes logarithmic time. This
saves a factor $k(\rho^{2}/\varepsilon)^{k\delta}/\log n$ in our time bound.
Second, instead of computing the whole navigating net, which takes time
$O(n\log\Phi(Y))$, we show how to compute any layer in $O(n\log n)$ time using
ANN queries. As our algorithm only requires $k$ layers of the navigating net,
it removes the dependency on $\Phi(Y)$ from the time bound.

Finally, in Section 6, we show how to extend our results on the
$\rho$-distortion problem to the minimum distortion problem. For the hardness
result, it suffices to add an extra point to $X$ and $Y$ that is far enough
from the other points, in order to make the reduction work. For the
algorithms, we use the reduction by Kenyon et al. [21] of the minimum
distortion problem to the $\rho$-distortion problem, which we speed-up using
exponential search, and using a well-separated pairs decomposition, which
allows us to reduce the number of candidate values for
$\operatorname{dist}(X,Y)$.

##  2 Notation and Preliminary

Let $(S,d_{S})$ be a finite metric space. The ball $\operatorname{b}(x,r)$
centered at $x$ with radius $r$ is the set of points $x^{\prime}\in S$ such
that $d_{S}(x,x^{\prime})\leqslant r$. The minimum and maximum interpoint
distances in $S$ are denoted $\operatorname{dmin}(S)$ and
$\operatorname{diam}(S)$, respectively. In other words,
$\operatorname{diam}(S)$ is the diameter of $S$. The spread of $S$ is the
ratio $\Phi(S)=\operatorname{diam}(S)/\operatorname{dmin}(S)$. The distance
from a point $x$ to a subset $T$ of $S$ is $d_{S}(x,T)=\min_{t\in
T}d_{S}(x,t)$. The distance between two sets $T$ and $U$ is
$d_{S}(T,U)=\min_{t\in T,u\in U}d_{S}(t,u)$.

A metric space $(S,d_{S})$ has doubling dimension $\delta$ if any ball of
radius $r$ is contained in the union of at most $2^{\delta}$ balls of radius
$r/2$. When $\delta=O(1)$, we say that this space is doubling. This notion of
dimension generalizes the dimension of a Euclidean space: In particular, the
Euclidean space $\mathbb{R}^{d}$ has doubling dimension $O(d)$ [17]. In this
paper, we will consider spaces of constant doubling dimension, so we assume
that $\delta=O(1)$. We will need the following packing lemma:

######  Lemma 1 ([22]).

If a metric space $(S,d_{S})$ has doubling dimension $\delta$, then
$|S|\leqslant(4\Phi(S))^{\delta}$.

We will also make use of the fact that a product of doubling metrics is
doubling. It was probably known, but we could not find a reference, so we
include a proof.

######  Lemma 2.

Let $(S_{1},d_{1}),\dots,(S_{k},d_{k})$ be metric spaces with doubling
dimensions $\delta_{1},\dots,\delta_{k}$, respectively. Then the product
metric $(S,d_{S})$ where $S=S_{1}\times\dots\times S_{k}$ and

$d_{S}((u_{1},\dots,u_{k}),(v_{1},\dots,v_{k}))=\max(d_{1}(u_{1},v_{1}),\dots,d_{k}(u_{k},v_{k}))$

for all $(u_{1},\dots,u_{k})$, $(v_{1},\dots,v_{k})\in S$ has doubling
dimension $\delta_{1}+\dots+\delta_{k}$.

###### Proof.

We prove this lemma for $k=2$; the general case follows by induction on $k$.
Let $b=\operatorname{b}((u_{1},u_{2}),r)$ be a ball of radius $r$ in $S$. Then
$b=b_{1}\times b_{2}$ where $b_{1}=\operatorname{b}(u_{1},r)$ and
$b_{2}=\operatorname{b}(u_{2},r)$ are balls of radius $r$ in $S_{1}$ and
$S_{2}$, respectively. So $b_{1}$ is covered by $m_{1}\leqslant
2^{\delta_{1}}$ balls $b_{1}^{1},\dots,b_{1}^{m_{1}}$ balls of radius $r/2$.
Similarly, $b_{2}$ is covered by $m_{2}\leqslant 2^{\delta_{2}}$ balls
$b_{2}^{1},\dots,b_{2}^{m_{2}}$ balls of radius $r/2$.

Let $(v_{1},v_{2})\in b$. Then we have $d_{1}(u_{1},v_{1})\leqslant r$ and
$d_{2}(u_{2},v_{2})\leqslant r$, which means that $v_{1}\in b_{1}$ and
$v_{2}\in b_{2}$. So there exist $i_{1}$ and $i_{2}$ such that $v_{1}\in
b_{1}^{i_{1}}$ and $v_{2}\in b_{2}^{i_{2}}$. In other words, $(v_{1},v_{2})\in
b_{1}^{i_{1}}\times b_{2}^{i_{2}}$. We have just proved that $b$ is contained
in the union of the Cartesian products $b_{1}^{i_{1}}\times b_{2}^{i_{2}}$.
There are $m_{1}\times m_{2}\leqslant 2^{\delta_{1}+\delta_{2}}$ such
Cartesian products, and each one of them is a ball of radius $r/2$ in $S$.
Therefore, $(S,d)$ has doubling dimension $\delta_{1}+\delta_{2}$. ∎

We call a mapping $\sigma$ satisfying Equation (1) a $\rho$-matching from $X$
to $Y$. The distance between two matchings $\sigma$ and $\sigma^{\prime}$ from
$X$ to $Y$ is $d_{M}(\sigma,\sigma^{\prime})=\max_{x\in
X}d_{Y}(\sigma(x),\sigma^{\prime}(x))$. We denote by $x_{1},x_{2},\dots,x_{k}$
the $k$ elements of $X$. So a matching $\sigma:X\to Y$ can be identified with
a sequence of $k$ points $(y_{1},\dots,y_{k})$ where
$y_{1}=\sigma(x_{1}),\dots,y_{k}=\sigma(x_{k})$. In other words, the space of
matchings from $X$ to $Y$ can be identified with $(Y^{k},d_{M})$. Then it
follows from Lemma 2 that:

######  Corollary 3.

The space of matchings from $X$ to $Y$ has doubling dimension $k\delta$.

##  3 Reduction from $k$-clique

Given an integer $k\geqslant 1$ and a graph $G(V,E)$, the $k$-clique problem
is to decide whether there exists a subset $C\subseteq V$ of $k$ vertices such
that any two of these vertices are connected by an edge in $E$. This subset
$C$ is called a $k$-clique. In this section, we present a reduction from the
$k$-clique problem to the $\rho$-distortion problem.

###  3.1 Construction

So let $G(V,E)$ be an instance of $k$-clique with $m$ vertices. We denote
$V=\\{v_{1},\dots,v_{m}\\}$, and we assume that $m\geqslant 24$. For any
$\rho\geqslant 1$, we will show how to construct an equivalent instance of the
$\rho$-distortion problem consisting of two metric spaces $(X,d_{X})$ and
$(Y,d_{Y})$ of respective sizes $k$ and $km$, and of doubling dimension
$\log_{2}3$.

Figure 1: (a) An $m$-point ring gadget $R_{i}$. (b) Two rings $R_{i}$ and
$R_{i^{\prime}}$. The edges drawn between $R_{i}$ and $R_{i^{\prime}}$
represent distances equal to ${2}^{\max(i,i^{\prime})}$, and the other
distances between a point in $R_{i}$ and a point in $R_{i^{\prime}}$ are
${2}^{\max(i,i^{\prime})}-(1/m)$

Let us first build $Y$, and its associated metric $d_{Y}$. We define an
$m$-point ring gadget (Figure 1a) as a set $R_{i}=\\{p_{i1},..,p_{im}\\}$ of
$m$ points spaced regularly on a circle of perimeter 1, so that the distance
between any two points is the usual distance along this circle:

$d_{Y}(p_{ij},p_{ij^{\prime}})=(1/m)\cdot\min{(j^{\prime}-j,m+j-j^{\prime})}\text{
whenever }1\leqslant j\leqslant j^{\prime}\leqslant m.$ (2)

Figure 2: Distances between rings. For instance, the distance between any
point in $R_{2}$ and any point in $R_{4}$ is ${2}^{4}$ or ${2}^{4}-(1/m)$.

We define $Y=R_{1}\cup\dots\cup R_{k}$ as the disjoint union of $k$ ring
gadgets. The index $j$ of the point $p_{ij}\in R_{i}$ corresponds to the
vertex $v_{j}$ of $G$. The distances between points in different rings are
defined as follows. (See Figure 1b and Figure 2.) For any
$i,i^{\prime}\in\\{1,\dots,k\\}$ such that $i\neq i^{\prime}$, and for any
$j,j^{\prime}\in\\{1,\dots,m\\}$,

$d_{Y}(p_{ij},p_{i^{\prime}j^{\prime}})=\begin{cases}{2}^{\max(i,i^{\prime})}&\text{if
}(v_{j},v_{j^{\prime}})\in E,\text{ and}\\\
{2}^{\max(i,i^{\prime})}-(1/m)&\text{otherwise}.\end{cases}$ (3)

Thus, the distance between two vertices in different rings is a power of ${2}$
if and only if their associated vertices in $G$ are connected by an edge. The
distance between two points in the same ring is given by Equation (2).

The pattern set $X$ consists of $k$ distinct points $x_{1},\dots,x_{k}$. We
associate with $X$ the following distance function $d_{X}$:

$d_{X}(x_{i},x_{i^{\prime}})=\begin{cases}0&\text{if $i=i^{\prime}$, and}\\\
{2}^{\max(i,i^{\prime})}\rho&\text{otherwise.}\end{cases}$

We now prove that the two spaces $(X,d_{X})$ and $(Y,d_{Y})$ are metric spaces
with doubling dimension $\log_{2}3$.

######  Lemma 4.

The distance function $d_{Y}:Y^{2}\to\mathbb{R}$ defined above is a metric.

###### Proof.

First observe that $d_{Y}$ is non-negative, symmetric, and that
$d_{Y}(p,p^{\prime})=0$ if and only if $p=p^{\prime}$. Within one ring,
$d_{Y}$ respects the triangle inequality because it is the usual distance
along a circle. So we only need to prove that the triangle inequality holds
for points lying in two or three different rings. Let
$i,i^{\prime},i^{\prime\prime}\in\\{1,\dots,k\\}$ and
$j,j^{\prime},j^{\prime\prime}\in\\{1,\dots,m\\}$. If $i$, $i^{\prime}$ and
$i^{\prime\prime}$ are distinct, we have

$d_{Y}(p_{ij},p_{i^{\prime}j^{\prime}})+d_{Y}(p_{i^{\prime}j^{\prime}},p_{i^{\prime\prime}j^{\prime\prime}})\geqslant{2}^{\max(i,i^{\prime})}+{2}^{\max(i^{\prime},i^{\prime\prime})}-\frac{2}{m}\geqslant{2}^{\max(i,i^{\prime\prime})}\geqslant
d_{Y}(p_{ij},p_{i^{\prime\prime}j^{\prime\prime}}).$

If $i\neq i^{\prime}$ and $j^{\prime}\neq j^{\prime\prime}$, then we have

$d_{Y}(p_{ij},p_{i^{\prime}j^{\prime}})+d_{Y}(p_{i^{\prime}j^{\prime}},p_{i^{\prime}j^{\prime\prime}})\geqslant\left({2}^{\max(i,i^{\prime})}-\frac{1}{m}\right)+\frac{1}{m}\geqslant{2}^{\max(i,i^{\prime})}\geqslant
d_{Y}(p_{ij},p_{i^{\prime}j^{\prime\prime}}).$

If $i\neq i^{\prime}$, then we have

$d_{Y}(p_{ij},p_{i^{\prime}j^{\prime}})+d_{Y}(p_{i^{\prime}j^{\prime}},p_{ij^{\prime\prime}})\geqslant
2\left({2}^{\max(i,i^{\prime})}-\frac{1}{m}\right)\geqslant{2}^{i}\geqslant
d_{Y}(p_{ij},p_{ij^{\prime\prime}}).$

This completes the proof that $d_{Y}$ is a metric. ∎

######  Lemma 5.

The metric space $(Y,d_{Y})$ has doubling dimension $\log_{2}3$.

###### Proof.

Let $b_{0}=\operatorname{b}(p,r)$ be an arbitrary ball, where $p=p_{ij}$. We
assume that the index $j^{\prime}$ in $p_{ij^{\prime}}$ is defined modulo $m$,
hence $p_{ij^{\prime}}=p_{i(j^{\prime}+m)}$. Let
$t=\lfloor\log_{2}(r+1/m)\rfloor$.

Suppose that $t<i$. Then $\log_{2}(r+1/m)<i$, and thus $r<{2}^{i}-1/m$. As the
distance from any point in $R_{i}$ to any point in another ring is at least
${2}^{i}-1/m$, it follows that $b_{0}\subseteq R_{i}$.

  * •

If $\lceil mr/2\rceil\geqslant m/4$, then $1+mr/2\geqslant m/4$. As
$m\geqslant 24$, it implies that $r\geqslant 5/12$. Let $c=p_{i1}$,
$c^{\prime}=p_{i(1+\lfloor m/3\rfloor)}$ and $c^{\prime\prime}=p_{i(1-\lfloor
m/3\rfloor)}$. As $m\geqslant 24$, we have $7m/24\leqslant\lfloor
m/3\rfloor\leqslant m/3$, and thus
$d_{Y}(c,c^{\prime})=d_{Y}(c,c^{\prime\prime})\leqslant 1/3$ and
$d_{Y}(c^{\prime},c^{\prime\prime})\leqslant 1-2\cdot(7/24)=5/12$. So we can
cover the whole ring $R_{i}$ (and thus $b_{0}$) with the balls of radius $r/2$
centered at $c$, $c^{\prime}$ and $c^{\prime\prime}$.

  * •

Otherwise, the points $p^{\prime}=p_{i(j+\lceil mr/2\rceil)}$ and
$p^{\prime\prime}=p_{i(j-\lceil mr/2\rceil)}$ lie in the open half-circle
centered at $p$, that is, $d_{Y}(p,p^{\prime})<1/4$ and
$d_{Y}(p,p^{\prime\prime})<1/4$. As $mr/2\leqslant\lceil mr/2\rceil<mr/2+1$,
it implies that $r/2\leqslant d_{Y}(p,p^{\prime})<r/2+1/m$ and $r/2\leqslant
d_{Y}(p,p^{\prime\prime})<r/2+1/m$. Therefore
$\operatorname{b}(p^{\prime},r/2)\cup\operatorname{b}(p^{\prime\prime},r/2)$
contains $b_{0}\setminus\\{p\\}$, and thus
$b_{0}\subseteq\operatorname{b}(p,r/2)\cup\operatorname{b}(p^{\prime},r/2)\cup\operatorname{b}(p^{\prime\prime},r/2)$.

Now suppose that $t\geqslant i$ and $t\geqslant 3$. As
$t\leqslant\log_{2}(r+1/m)<t+1$, we have ${2}^{t}-1/m\leqslant
r<{2}^{t+1}-1/m$. For any $i^{\prime}\geqslant t+1$, the distance from $p$ to
any point in $R_{i^{\prime}}$ is at least ${2}^{t+1}-1/m$, so $b_{0}\subseteq
R_{1}\cup\dots\cup R_{t}$. We consider three balls $b_{1}$, $b_{t-1}$ and
$b_{t}$ of radius $r/2$ centered at an arbitrary point $c_{1}\in R_{1}$,
$c_{t-1}\in R_{t-1}$ and $c_{t}\in R_{t}$, respectively. Since
$\frac{r}{2}>\frac{{2}^{t}}{2}-\frac{1}{2m}\geqslant\frac{2}{2}-\frac{1}{48}=\frac{47}{48},$
and $\operatorname{diam}(R_{t})\leqslant 1/2$, we have $R_{t}\subseteq b_{t}$.
Similarly, we have $R_{t-1}\subseteq b_{t-1}$. Let $q\in R_{1}\cup\dots\cup
R_{t-2}$. Then we have $d_{Y}(c_{1},q)\leqslant 2^{t-2}$. As $r/2\geqslant
2^{t-1}-1/(2m)\geqslant 2^{t-2}$, it follows that $q\in b_{1}$. So we just
proved that $R_{1}\cup\dots\cup R_{t-2}\subseteq b_{1}$. Therefore,
$b_{0}\subseteq R_{1}\cup\dots\cup R_{t}\subseteq b_{1}\cup b_{t-1}\cup
b_{t}$.

Similarly, if $t\geqslant i$ and $t=1$, then $t=1=i$ and we have
$b_{0}\subseteq b_{1}$. If $t\geqslant i$ and $t=2$, then $b_{0}\subseteq
b_{1}\cup b_{2}$.

In any case, $b_{0}$ is contained in the union of at most 3 balls of half its
radius, and thus $(Y,d_{y})$ has doubling dimension $\log_{2}3$. ∎

A simpler version of the arguments in the proofs of Lemma 4 and 5 yields the
following.

######  Lemma 6.

$(X,d_{X})$ is a metric space of doubling dimension $\log_{2}3$.

###  3.2 Proof of Correctness

We now prove that our reduction of $k$-clique to the $\rho$-distortion problem
is correct. So given an instance $G=(V,E)$ of the $k$-clique problem, we
construct the metric spaces $(X,d_{X})$ and $(Y,d_{Y})$ as described above.
These two metric spaces form an instance of the $\rho$-distortion problem. We
need to show that these two instances of $k$-clique and $\rho$-distortion are
equivalent.

We first assume that $G$ is a positive instance of $k$-clique. So there is a
clique $\\{v_{j_{1}},\dots,v_{j_{k}}\\}$ of size $k$ in $G$. Let $\sigma:X\to
Y$ be the matching defined by $\sigma(x_{i})=p_{ij_{i}}$ for all
$i\in\\{1,\dots,k\\}$. For any $i\neq i^{\prime}$, there is an edge between
$v_{j_{i}}$ and $v_{j_{i^{\prime}}}$ in $G$. Therefore, we have
$d_{Y}(p_{ij_{i}},p_{i^{\prime}j_{i^{\prime}}})={2}^{\max(i,i^{\prime})}=\frac{1}{\rho}d_{X}(x_{i},x_{i^{\prime}}).$
It means that
$d_{Y}(\sigma(x_{i}),\sigma(x_{i^{\prime}}))=\frac{1}{\rho}d_{X}(x_{i},x_{i^{\prime}})$,
and thus $\sigma$ is a solution to our instance of the $\rho$-distortion
problem.

We now prove the converse. Let $\sigma:X\to Y$ be a solution to our instance
of the $\rho$-distortion problem. We denote $y_{i}=\sigma(x_{i})$ for each
$x_{i}\in X$. Each $y_{i}$ corresponds to a vertex $w_{i}\in V$. More
precisely, we have $y_{i}=p_{i^{\prime}j^{\prime}}$ for some indices
$i^{\prime},j^{\prime}$, and we set $w_{i}=v_{j^{\prime}}$. We want to show
that $w_{1},\dots,w_{k}$ is a $k$-clique in $G$.

######  Lemma 7.

Each ring contains exactly one point $y_{i}$, and $y_{i}\in R_{i}$ whenever
$i\geqslant 3$.

###### Proof.

Let $1\leqslant i_{1}<i_{2}\leqslant k$. Then
$d_{X}(x_{i_{1}},x_{i_{2}})={2}^{i_{2}}\rho$. As $\sigma$ is a solution to the
$\rho$-distortion problem, it implies that
$d_{Y}(y_{i_{1}},y_{i_{2}})=d_{Y}(\sigma(x_{i_{1}}),\sigma(x_{i_{2}}))\geqslant{2}^{i_{2}}\geqslant
4.$ As the diameter of each ring is at most $1/2$, it follows that $y_{i_{1}}$
and $y_{i_{2}}$ lie in different rings. Hence each ring contains exactly one
point $y_{i}$.

We now prove the second part of the lemma, so we may assume that $k\geqslant
3$. The point $x_{k}$ is at distance ${2}^{k}\rho$ from each point
$x_{1},\dots,x_{k-1}$. As $\sigma$ is a solution to the $\rho$-distortion
problem, the distance from $y_{k}$ to any point $y_{1},\dots,y_{k-1}$ is at
least ${2}^{k}$, and since it is the largest interpoint distance in $Y$, it
implies that $d_{Y}(y_{i},y_{k})={2}^{k}$ for any $1\leqslant i<k$. This is
only possible if $y_{k}\in R_{k}$.

Applying the same argument repeatedly to $y_{k-1},\dots,y_{3}$ shows that
$y_{i}\in R_{i}$ for any $i\in\\{3,\dots,k-1\\}$. (As
$d_{X}(x_{1},x_{i})=d_{X}(x_{2},x_{i})$ for any $i\geqslant 3$, the points
$x_{1}$ and $x_{2}$ are indistinguishable, hence we have either $y_{1}\in
R_{1}$ and $y_{2}\in R_{2}$, or $y_{1}\in R_{2}$ and $y_{2}\in R_{1}$.) ∎

We can now prove that:

######  Lemma 8.

The vertices $\\{w_{1},..,w_{k}\\}$ form a $k$-clique in $G$.

###### Proof.

Let $1\leqslant i^{\prime}<i\leqslant k$. By Equation (3), the distance
$d_{Y}(y_{i},y_{i^{\prime}})$ can take only two values: it is equal to
${2}^{i}$ if $(w_{i},w_{i^{\prime}})\in E$, and ${2}^{i}-(1/m)$ otherwise. As
$\sigma$ is a solution to the $\rho$-distortion problem, we also have
$d_{Y}(y_{i},y_{i^{\prime}})=d_{Y}(\sigma(x_{i}),\sigma(x_{i^{\prime}}))\geqslant\frac{1}{\rho}d_{X}(x_{i},x_{i^{\prime}})={2}^{i}.$
It implies that $d_{Y}(y_{i},y_{i^{\prime}})={2}^{i}$, and that
$(w_{i},w_{i^{\prime}})\in E$. ∎

In summary, we obtained the following result.

######  Theorem 9.

The graph $G$ admits a $k$-clique if and only if $(X,d_{X})$, $(Y,d_{Y})$ is a
positive instance of the $\rho$-distortion problem.

###  3.3 Consequences

The construction of $(X,d_{X})$ and $(Y,d_{Y})$ from $G$ is performed in
polynomial time. It is also an FPT-reduction with parameter $k$. As $k$-clique
is NP-complete and $W[1]$-hard, it follows that:

######  Corollary 10.

The $\rho$-distortion problem for doubling spaces of dimension $\log_{2}3$ is
NP-hard, and is $W[1]$-hard when parameterized by $k$.

Unless the exponential time hypothesis (ETH) is false, our reduction also
shows that that the $\rho$-distortion problem cannot be solved in time
$n^{o(k)}$. More precisely, it follows from a known hardness result on
$k$-clique [12, Theorem 14.21] that:

######  Corollary 11.

The $\rho$-distortion problem for doubling spaces of dimension $\log_{2}3$
cannot be solved in time $f(k)\cdot n^{o(k)}$ for any computable function $f$,
unless ETH is false.

##  4 Approximation algorithm for the $\rho$-distortion problem

In this section, we present an approximation algorithm for the
$\rho$-distortion problem. We assume that the doubling dimension is constant,
that is, $\delta=O(1)$. As we saw in Section 3, the $\rho$-distortion problem
is hard, so we relax the problem slightly: Given a parameter
$0<\varepsilon\leqslant 1$, the $(\rho,\varepsilon)$-distortion problem is to
find a $(1+\varepsilon)\rho$-matching whenever a $\rho$-matching exists. If
there is no $\rho$-matching, then our algorithm either returns a
$(1+\varepsilon)\rho$-matching, or it does not return any result.

###  4.1 Navigating nets

Our algorithm records $Y$ in a navigating net, which is a data structure
representing $Y$ at different resolutions. (See Figure 3.) This structure was
introduced by Krauthgamer and Lee [24]. We will use a slightly modified
version of it. Our version of the navigating net has the advantages that each
layer is an $r$-net of $Y$ (see definition below), and that it can easily be
computed in $O(n\log\Phi(Y))$ time. On the other hand, it does not allow
efficient deletions. Several other variations exists. In particular, Cole and
Gottlieb [9] and Har-Peled and Mendel [20] presented more involved data
structures whose time bounds do not depend on the spread.

Figure 3: A navigating net at three consecutive scales $r/2$, $r$ and $2r$.
The white dots represent points of the $r$-nets $Y_{r/2}$, $Y_{r}$ and
$Y_{2r}$. The ellipses represent the balls of radius $r/2$, $r$ and $2r$
centered at the points in $Y_{r/2}$, $Y_{r}$ and $Y_{2r}$, respectively. The
blue edges are the horizontal edges incident to one of the vertices of
$Y_{r/2}$. The red edges are the vertical edges.

For any $r\geqslant 0$, an $r$-net of $Y$ is a subset $Y_{r}\subseteq Y$ such
that $\operatorname{dmin}(Y_{r})\geqslant r$ and $Y\subseteq\bigcup_{y\in
Y_{r}}\operatorname{b}(y,r)$. An $r$-net can be constructed incrementally by
repeatedly adding new points that lie outside of the current union of balls,
until $Y$ is completely covered. Intuitively, an $r$-net represents $Y$ at
resolution $r$.

A scale is a rational number $r=2^{i}$ such that $i\in\mathbb{Z}$. Our
navigating net records a sequence of $r$-nets
$Y_{r_{\min}},Y_{2r_{\min}},\dots,Y_{r_{\max}}$ such that $r_{\min}$ and
$r_{\max}$ are scales satisfying the inequalities
$\operatorname{dmin}(Y)/2<r_{\min}\leqslant\operatorname{dmin}(Y)$ and
$\operatorname{diam}(Y)/2\leqslant r_{\max}<2\operatorname{diam}(Y)$. At the
largest scale $r_{\max}$, the $r_{\max}$-net $Y_{r_{\max}}$ consists of a
single point $y_{\mathrm{root}}$, and thus
$Y=\operatorname{b}(y_{\mathrm{root}},r_{\max})$. At the lowest scale, we set
$Y_{r_{\min}}=Y$. So the navigating net represents $Y$ at all scales $r$ such
that $r_{\min}\leqslant r\leqslant r_{\max}$ by an $r$-net $Y_{r}$. All scales
$r>r_{\max}$ are represented by $Y_{r}=Y_{r_{\max}}=\\{y_{\mathrm{root}}\\}$,
but we do not construct these copies of $Y_{r_{\max}}$ explicitly.

We may assume that $\log(\Phi(Y))\geqslant 1$ as otherwise, $|Y|=O(1)$ by
Lemma 1, and the $\rho$-distortion problem can be solved in $O(1)$ time by
brute force. So we only construct $Y_{r}$ at $O(\log(\Phi(Y))$ different
scales $r$ such that $\operatorname{dmin}(Y)/2<r_{\min}\leqslant r\leqslant
r_{\max}<2\operatorname{diam}(Y)$.

We construct a graph over these $r$-nets. First, at each scale $r$ such that
$r_{\min}\leqslant r\leqslant r_{\max}$, we connect any two nodes $y$,
$y^{\prime}\in Y_{r}$ such that $d_{Y}(y,y^{\prime})\leqslant 6r$ by an edge,
that we call a horizontal edge. At each scale $r$ such that
$2r_{\min}\leqslant r\leqslant r_{\max}$, we also connect with a vertical edge
each $y\in Y_{r/2}$ to a node $p(y)\in Y_{r}$, called the parent of $y$, such
that $d_{Y}(y,p(y))\leqslant r$. At least one such node $p(y)$ exists since
$Y_{r}$ is an $r$-net. We call $y$ a child of $p(y)$. More generally, we say
that $y$ is a descendant of $y^{\prime}$ if $y$ is a child of $y^{\prime}$, or
$y$ is a child of a descendant of $y^{\prime}$. Conversely, we say that
$y^{\prime}$ is an ancestor of $y$ if $y$ is a descendant of $y^{\prime}$.
When $r^{\prime}>r_{\max}$, the node $y_{\mathrm{root}}\in Y_{r^{\prime}}$ is
an ancestor of any node at a lower level.

These $r$-nets, together with the horizontal and vertical edges, form our
navigating net. The lemma below shows that it has bounded degree.

######  Lemma 12.

Each node of the navigating net has degree $O(1)$.

###### Proof.

The children of $y\in Y_{r}$ are at distance at most $r$ from $y$, and thus at
distance at most $2r$ from each other. As they are in $Y_{r/2}$, they are at
distance at least $r/2$ from each other. So the spread of the set of children
of $y$ is at most 4, thus $y$ has at most $16^{\delta}$ children by Lemma 1.
It follows that any node is incident to at most $1+16^{\delta}=O(1)$ vertical
edges. The nodes adjacent to $y$ via horizontal edges are at distance at most
$6r$ from $y$ and at least $r$ from each other. Hence their spread is at most
$12$, so there are at most $48^{\delta}$ of them. ∎

The lemma below shows that our navigating net can be computed efficiently.

######  Lemma 13.

The navigating net of $Y$ can be computed in $O(n\log\Phi(Y))$ time.

###### Proof.

We construct the navigating net incrementally, adding the nodes one-by-one in
a top-down manner. We first pick an arbitrary $y_{\mathrm{root}}\in Y$, and we
let $r_{\max}$ be the smallest scale that is at least $\max_{y\in
Y}d_{Y}(y_{\mathrm{root}},y)$. We initialize our navigating net with the
single node $y_{\mathrm{root}}$ at scale $r_{\max}$. Initially, we set
$r_{\min}=r_{\max}$, and we will decrease $r_{\min}$ when we insert new
points, if needed. We will maintain the invariants that any two inserted
points are at distance at least $r_{\min}$ from each other, and that
$Y_{r_{\min}}$ is the set of all inserted points.

We repeatedly insert an arbitrary point $y\in Y$, until all points have been
inserted. We now show how to update the data structure while inserting $y$.

We start from the root and move down to scale $r_{\min}$. During this process,
at the current scale $r$, we maintain the set $N(y,r)$ of points in $Y_{r}$
that are at distance at most $6r$ from $y$. We have
$\operatorname{dmin}(N(y,r))\geqslant r$ and
$\operatorname{diam}(N(y,r))\leqslant 12r$, hence $|N(y,r)|\leqslant
48^{\delta}=O(1)$ by Lemma 1.

We now show how to maintain $N(y,r)$ while traversing the data structure from
top to bottom. So we assume that we know $N(y,2r)$, and we show how to find
$N(y,r)$. Suppose that $y^{\prime}\in N(y,r)$. Then we have
$d_{Y}(y^{\prime},y)\leqslant 6r$ and
$d_{Y}(y^{\prime},p(y^{\prime}))\leqslant 2r$. Therefore,
$d_{Y}(y,p(y^{\prime}))\leqslant 8r$, which implies that $p(y^{\prime})\in
N(y,2r)$. So in order to find the points in $N(y,r)$, we can simply go through
all the children of the nodes in $N(y,2r)$, and check for each of them if it
is at distance at most $6r$ from $y$. It takes $O(1)$ time as the navigating
net has constant degree. So we can compute all the sets $N(y,r)$ such that
$r_{\min}\leqslant r\leqslant r_{\max}$ in $O(\log\Phi(Y))$ time as there are
$O(\log\Phi(Y))$ different scales.

If $N(y,r_{\min})$ is empty, then the closest inserted point $y_{c}$ to $y$ is
at distance more than $r_{\min}$, so we don’t need to update $r_{\min}$.
Otherwise, since $y_{c}\in N(y,r_{\,min})$, we can find $y_{c}$ in $O(1)$ time
by brute force. If $d_{Y}(y,y_{c})<r_{\min}$, we set $r_{\mathrm{temp}}$ to be
the largest scale that is at most $d_{Y}(y,y_{c})$, and we create new levels
in the navigating net at all scales $r$ such that $r_{\mathrm{temp}}\leqslant
r<r_{\min}$. As $Y_{r_{\min}}$ consists of all the inserted points, the same
is true for $Y_{r}$ at the new levels, and we can compute the horizontal edges
in $O(1)$ time per node as they are a subset of the horizontal edges at scale
$2r$. We also connect each node by a vertical edge to its copy at scale $2r$.
After updating all levels, we set $r_{\min}=r_{\mathrm{temp}}$.

We now show how to insert $y$ in the navigating net. When inserting $y$ at
level $r$, we check whether all the points in $N(y,r)$ are at distance at
distance at least $r$ from $y$. If it is the case, we insert $y$ into $Y_{r}$,
and we connect $y$ to all the points in $N(y,r)$ by a horizontal edge. Then
let $y^{\prime\prime}$ be the point of $N(y,2r)$ that is closest to $y$. As
$Y_{2r}$ is a $2r$-net, we have $d_{Y}(y,y^{\prime\prime})\leqslant 2r$, so we
set $y^{\prime\prime}=p(y)$.

This algorithm constructs the navigating net incrementally in $O(\log\Phi(Y))$
time per point in $Y$, so the overall running time is $O(n\log\Phi(Y))$. Our
construction ensures that
$\operatorname{dmin}(Y)/2<r_{\min}\leqslant\operatorname{dmin}(Y)$ and
$\operatorname{diam}(Y)/2\leqslant r_{\max}<2\operatorname{diam}(Y)$.

∎

We associate the ball $\operatorname{b}(y,2r)$ with each node $y$ in an
$r$-net $Y_{r}$. These balls have the following properties.

######  Lemma 14.

At any scale $r$, and for any subset $S\subseteq Y$ such that
$\operatorname{diam}(S)\leqslant r$, there exists a node $y\in Y_{r}$ such
that $S\subseteq\operatorname{b}(y,2r)$.

###### Proof.

Let $z$ be a point in $S$. Since $Y_{r}$ is an $r$-net, there exists $y\in
Y_{r}$ such that $d_{Y}(y,z)\leqslant r$. For any $z^{\prime}\in S$, we have
$d_{Y}(z,z^{\prime})\leqslant\operatorname{diam}(S)\leqslant r$ so by the
triangle inequality, $d_{Y}(y,z^{\prime})\leqslant 2r$. In other words,
$S\subseteq\operatorname{b}(y,2r)$. ∎

The lemma below shows that any ball other than the root is contained in the
ball associated with any of its ancestors.

######  Lemma 15.

Let $r$ and $r^{\prime}$ be two scales such that $r<r^{\prime}$. For any $y\in
Y_{r}$, and for any ancestor $y^{\prime}\in Y_{r^{\prime}}$ of $y$, we have
$b(y,2r)\subseteq\operatorname{b}(y^{\prime},2r^{\prime})$.

###### Proof.

When $r^{\prime}=2r$, we have $y^{\prime}=p(y)$, so
$d_{Y}(y,y^{\prime})\leqslant r^{\prime}$ by definition of $p(y)$. It follows
that
$\operatorname{b}(y,2r)\subseteq\operatorname{b}(y^{\prime},2r+r^{\prime})=\operatorname{b}(y^{\prime},2r^{\prime})$.
When $r^{\prime}\geqslant 4r$, the result is obtained by induction. ∎

In summary, the balls $\operatorname{b}(y,2r)$, connected by the vertical
edges, form a tree such that each ball is contained in each of its ancestors.
The horizontal edges will help us traverse this tree within a given level.

###  4.2 Splitting the pattern

Our algorithm proceeds recursively, by partitioning $X$ into two well-
separated subsets $P$ and $Q$ at each stage. More precisely, we will split $X$
as follows. (Remember that $k=|X|$.)

######  Lemma 16.

If $k\geqslant 2$, we can partition $X$ into two non-empty subsets $P$ and $Q$
such that $\operatorname{diam}(X)\leqslant(k-1)\cdot d_{X}(P,Q)$.

###### Proof.

We obtain $P$ and $Q$ by running Kruskal’s algorithm [10] for computing a
minimum spanning tree of $X$, and stopping at the second-last step. So
starting from the forest $(X,\emptyset)$, we repeatedly insert the shortest
edge that connects any two trees of the current forest, until we are left with
exactly two trees $P$ and $Q$. At the last step, $P$ and $Q$ are then
connected with an edge of length $\ell=d_{X}(P,Q)$, which is not shorter than
any edge in the spanning trees we constructed for $P$ and $Q$. Thus,
$\operatorname{diam}(P)\leqslant(|P|-1)\ell$ and
$\operatorname{diam}(Q)\leqslant(|Q|-1)\ell$. It follows that
$\operatorname{diam}(X)\leqslant(|P|-1)\ell+\ell+(|Q|-1)\ell=(k-1)\ell$. ∎

###  4.3 Recording approximate matchings

Our algorithm for the $(\rho,\varepsilon)$-distortion problem records a
collection of approximate matchings for some layers of the navigating net.
More precisely, for some non-empty subset $W$ of $X$, for some
$0<\beta\leqslant 1$ and for some scales
$r\geqslant\rho\operatorname{diam}(W)$, we will construct a data structure
$L(W,\beta,r)$ that records at least one $(1+\beta)\rho$-matching from $W$ to
$Y$ if a $\rho$-matching from $W$ to $Y$ exists. In particular, we will
compute $L(X,\varepsilon,r_{X})$ at an appropriate scale $r_{X}$, which
records a solution to the $(\rho,\varepsilon)$-distortion problem if there is
one. We first give three invariants of $L(W,\beta,r)$, and in the next two
sections we will show how to compute $L(X,\varepsilon,r_{X})$ recursively.

The data structure $L(W,\beta,r)$ records a set $M(y,W,\beta,r)$ of matchings
at each node $y\in Y_{r}$. These sets satisfy the following properties.

######  Property 17.

Let $r$ be a scale such that $r\geqslant\rho\operatorname{diam}(W)$.

  1. (a)

For any $y\in Y_{r}$, each matching $\sigma_{\beta}\in M(y,W,\beta,r)$ is a
$(1+\beta)\rho$-matching from $W$ to $Y$ such that
$\sigma_{\beta}(W)\subseteq\operatorname{b}(y,3r)$.

  2. (b)

For any $y\in Y_{r}$ and any two distinct $\sigma_{\beta}$,
$\sigma_{\beta}^{\prime}\in M(y,W,\beta,r)$, we have
$d_{M}(\sigma_{\beta},\sigma_{\beta}^{\prime})\geqslant\beta r/(2\rho^{2})$.

  3. (c)

For any $\rho$-matching $\sigma:W\to Y$, there exist $y\in Y_{r}$ and
$\sigma_{\beta}\in M(y,W,\beta,r)$ such that
$d_{M}(\sigma,\sigma_{\beta})\leqslant\beta r/\rho^{2}$.

These properties imply the following bound on the sizes of these sets.

######  Lemma 18.

For any $y\in Y_{r}$, we have
$|M(y,W,\beta,r)|=O((\rho^{2}/\beta)^{k\delta})$.

###### Proof.

Property 17a implies that any two matchings in $M(y,W,\beta,r)$ are at
distance at most $6r$ from each other, that is,
$\operatorname{diam}(M(y,W,\beta,r))\leqslant 6r$. Property 17b means that
$\operatorname{dmin}(M(y,W,\beta,r))\geqslant\beta r/(2\rho^{2})$. Therefore,
$\Phi(M(y,W,\beta,r))\leqslant 12\rho^{2}/\beta$. Then by Lemma 1 and
Corollary 3, we have $|M(y,W,\beta,r)|\leqslant(48\rho^{2}/\beta)^{k\delta}$.
∎

###  4.4 Recursive construction

Our algorithm constructs $L(W,\beta,r)$ recursively, for some subsets $W$ of
$X$ and some values $\beta$ and $r$. We start with the base case, then we show
how to compute $L(W,\beta,r^{\prime})$ from $L(W,\beta,r)$ when
$\rho\operatorname{diam}(W)\leqslant r<r^{\prime}$. Finally, we show how to
recursively compute $L(X,\varepsilon,r_{X})$ at an appropriate scale $r_{X}$
by splitting $X$ into $P$ and $Q$ according to Lemma 16. (See Figure 4.)

Figure 4: Recursive construction of $L(X,\varepsilon,r_{X})$. The pattern $X$
is split into $P$ and $Q$. (a) Bottom-up phase: We compute $L(P,\beta,r_{X})$
and $L(Q,\beta,r_{X})$ from $L(P,\beta,r_{P})$ and $L(P,\beta,r_{Q})$. (b) An
approximate matching $\bar{\sigma}_{\varepsilon}:X\to Y$ in
$L(X,\varepsilon,r_{X})$ is obtained by combining two matchings
$\bar{\sigma}_{\beta}^{P}:P\to Y$ and $\bar{\sigma}_{\beta}^{Q}:Q\to Y$.

#### Base case where $|W|=1$ and $r=r_{\min}$.

We have $W=\\{x_{i}\\}$, so a $(1+\beta)\rho$-matching simply maps $x_{i}$ to
any element of $Y$. Therefore, at scale $r=r_{\min}$, we have $Y_{r}=Y$, so
for each $y\in Y$, we record the matching that sends $x_{i}$ to $y$ in
$M(y,W,\beta,r_{\min})$.

#### Bottom-up construction.

Suppose that $L(W,\beta,r)$ has been computed for some scale
$r\geqslant\rho\operatorname{diam}(W)$. Given a larger scale $r^{\prime}>r$,
we now show how to construct $L(W,\beta,r^{\prime})$. (See Algorithm 1.)

Algorithm 1 Computing $L(W,\beta,r^{\prime})$ from $L(W,\beta,r)$

1:procedure BottomUpConstruction

2: for $y^{\prime}\in Y_{r^{\prime}}$ do

3: $M(y^{\prime},W,\beta,r^{\prime})\leftarrow\emptyset$

4: for each descendant $y$ of $y^{\prime}$ at scale $r$ do

5: for each matching $\sigma_{\beta}\in M(y,W,\beta,r)$ do

6: if $d_{M}(\sigma_{\beta},\sigma^{\prime}_{\beta})\geqslant\beta
r^{\prime}/(2\rho^{2})$ for all $\sigma^{\prime}_{\beta}\in
M(y^{\prime},W,\beta,r^{\prime})$ then

7: insert $\sigma_{\beta}$ into $M(y^{\prime},W,\beta,r^{\prime})$

For each $y^{\prime}\in Y_{r^{\prime}}$, we proceed as follows. Initially, we
set $M(y^{\prime},W,\beta,r^{\prime})=\emptyset$. Then for each descendant
$y\in Y_{r}$ of $y^{\prime}$, and for each matching $\sigma_{\beta}\in
M(y,W,\beta,r)$, we check by brute force whether
$d_{M}(\sigma_{\beta},\sigma_{\beta}^{\prime})\geqslant\beta
r^{\prime}/(2\rho^{2})$ for each matching $\sigma_{\beta}^{\prime}$ that was
previously inserted into $M(y^{\prime},W,\beta,r^{\prime})$. If it is the
case, we insert $\sigma_{\beta}$ into $M(y^{\prime},W,\beta,r^{\prime})$, and
otherwise we discard $\sigma_{\beta}$. It ensures that Property 17b holds for
$M(y^{\prime},W,\beta,r^{\prime})$.

We now prove that Property 17a holds. Let $\sigma_{\beta}\in
M(y^{\prime},W,\beta,r^{\prime})$. By construction, $\sigma_{\beta}\in
M(y,W,\beta,r)$ for some descendant $y$ of $y^{\prime}$. As Property 17a holds
at scale $r$, it follows that $\sigma_{\beta}$ is a $(1+\beta)\rho$-matching,
and that $\sigma_{\beta}(W)\subseteq\operatorname{b}(y,3r)$. Since $y$ is a
descendant of $y^{\prime}$, by Lemma 15, we have
$\operatorname{b}(y,2r)\subseteq\operatorname{b}(y^{\prime},2r^{\prime})$, and
thus
$\operatorname{b}(y,3r)\subseteq\operatorname{b}(y^{\prime},2r^{\prime}+r)\subseteq\operatorname{b}(y^{\prime},3r^{\prime})$.
It follows that
$\sigma_{\beta}(W)\subseteq\operatorname{b}(y^{\prime},3r^{\prime})$, and thus
Property 17a holds for $M(y^{\prime},W,\beta,r^{\prime})$.

Finally, we prove that Property 17c holds as well. Let $\sigma:W\to Y$ be a
$\rho$-matching. As Property 17c holds for $L(W,\beta,r)$, there must be a
node $y\in Y_{r}$ and $\sigma_{\beta}\in M(y,W,\beta,r)$ such that
$d_{M}(\sigma,\sigma_{\beta})\leqslant\beta r/\rho^{2}$. Let $y^{\prime}\in
Y_{r^{\prime}}$ be the ancestor of $y$ at scale $r^{\prime}$. If
$\sigma_{\beta}$ was inserted into $M(y^{\prime},W,\beta,r^{\prime})$, then we
are done because $d_{M}(\sigma,\sigma_{\beta})\leqslant\beta
r/\rho^{2}\leqslant\beta r^{\prime}/\rho^{2}$. Otherwise, we must have
inserted a matching $\sigma_{\beta}^{\prime}$ such that
$d_{M}(\sigma_{\beta},\sigma_{\beta}^{\prime})<\beta r^{\prime}/(2\rho^{2})$.
It follows that $d_{M}(\sigma,\sigma_{\beta}^{\prime})\leqslant
d_{M}(\sigma,\sigma_{\beta})+d_{M}(\sigma_{\beta},\sigma_{\beta}^{\prime})\leqslant\beta(r+r^{\prime}/2)/\rho^{2}\leqslant\beta
r^{\prime}/\rho^{2}$, which completes the proof that Property 17 holds for
$L(W,\beta,r^{\prime})$.

We now analyze this algorithm. First we need to find all the descendants of
each node $y^{\prime}\in Y_{r^{\prime}}$. We can do this by traversing the
navigating net, which takes time $O(n\log\Phi(Y))$ as there are
$O(\log\Phi(Y))$ levels in the navigating net. By Lemma 18, we have
$|M(y^{\prime},W,\beta,r^{\prime})|=O((\rho^{2}/\beta)^{k\delta})$. Therefore,
each time we attempt to insert a matching $\sigma_{\beta}$ into
$M(y^{\prime},W,\beta,r^{\prime})$, we compare it with
$O((\rho^{2}/\beta)^{k\delta})$ previously inserted matchings, so it takes
$O(k(\rho^{2}/\beta)^{k\delta})$ time as the distance between two matchings
can be computed in $O(k)$ time. Since $|Y_{r}|\leqslant n$, and each set
$M(y,W,\beta,r)$ has cardinality $O((\rho^{2}/\beta)^{k\delta})$, we spend
$O(nk(\rho^{2}/\beta)^{2k\delta})$ time for computing $L(W,\beta,r^{\prime})$.
So we just proved the following.

######  Lemma 19.

Let $r$ and $r^{\prime}$ be two scales such that
$\rho\operatorname{diam}(W)\leqslant r<r^{\prime}$. Given $L(W,\beta,r)$, we
can compute $L(W,\beta,r^{\prime})$ in
$O(nk(\rho^{2}/\beta)^{2k\delta}+n\log\Phi(Y))$ time.

#### Computing $L(X,\varepsilon,r_{X})$ by splitting $X$.

Suppose that $k\geqslant 2$. Let $r_{X}$ be the smallest scale that is at
least as large as $\rho\operatorname{diam}(X)$, hence
$r_{X}=2^{\lceil\log_{2}(\rho\operatorname{diam}(X))\rceil}$. In particular,
we have $r_{X}/2<\rho\operatorname{diam}(X)\leqslant r_{X}$. If
$r_{X}<r_{\min}$, then we have
$\rho\operatorname{diam}(X)<r_{\min}\leqslant\operatorname{dmin}(Y)$, and
there cannot be any $\rho$-matching, so our algorithm does not return any
matching. Therefore, from now on, we may assume that $r_{X}\geqslant
r_{\min}$.

Let $P$ and $Q$ be the sets obtained by splitting $X$ as described in Lemma
16, so $\operatorname{diam}(X)\leqslant(k-1)\cdot d_{X}(P,Q)$. Let
$\beta=\varepsilon/(8k-8)$, and suppose that $L(P,\beta,r_{X})$ and
$L(Q,\beta,r_{X})$ have been computed earlier. We now show how to compute
$L(X,\varepsilon,r_{X})$.

For any two matchings $\sigma^{P}:P\to Y$ and $\sigma^{Q}:Q\to Y$, we denote
by $\sigma^{P}\cdot\sigma^{Q}$ the matching from $X$ to $Y$ whose restrictions
to $P$ and $Q$ are $\sigma^{P}$ and $\sigma^{Q}$, respectively. In other
words, if $\sigma=\sigma^{P}\cdot\sigma^{Q}$, then we have
$\sigma(x)=\sigma^{P}(x)$ for all $x\in P$ and $\sigma(x)=\sigma^{Q}(x)$ for
all $x\in Q$.

We compute $L(X,\varepsilon,r_{X})$ as follows. (See Algorithm 2.) For each
node $y\in Y_{r_{X}}$, we consider all the pairs of matchings consisting of a
matching $\sigma^{P}_{\beta}\in M(y^{P},P,\beta,r_{X})$ and a matching
$\sigma^{Q}_{\beta}\in M(y^{Q},Q,\beta,r_{X})$, where $y^{P}$ and $y^{Q}$ are
in $Y_{r_{X}}$ and are at distance at most $6r_{X}$ from $y$. Then we consider
$\sigma_{\varepsilon}=\sigma^{P}_{\beta}\cdot\sigma^{Q}_{\beta}$ as a
candidate for being inserted into $M(y,X,\varepsilon,r_{X})$. We first check
whether $\sigma_{\varepsilon}$ is a $(1+\varepsilon)\rho$-matching and
$\sigma_{\varepsilon}(X)\subseteq\operatorname{b}(y,3r_{X})$. If it is the
case, and if $\sigma_{\varepsilon}$ is at distance at least $\varepsilon
r_{X}/(2\rho^{2})$ from any matching previously inserted into
$M(y,X,\varepsilon,r_{X})$, we insert $\sigma_{\varepsilon}$ into
$M(y,X,\varepsilon,r_{X})$.

Algorithm 2 Computing $L(X,\varepsilon,r_{X})$ from $L(P,\beta,r_{X})$ and
$L(Q,\beta,r_{X})$

1:procedure ComputeBySplitting

2: for $y\in Y_{r_{X}}$ do

3: $M(y,X,\varepsilon,r_{X})\leftarrow\emptyset$

4: for each pair of horizontal edges $(y,y^{P})$, $(y,y^{Q})$ do

5: for each pair $\sigma^{P}_{\beta}\in M(y^{P},P,\beta,r_{X})$,
$\sigma^{Q}_{\beta}\in M(y^{Q},Q,\beta,r_{X})$ do

6: $\sigma_{\varepsilon}\leftarrow\sigma_{\beta}^{P}\cdot\sigma_{\beta}^{Q}$

7: if $\sigma_{\varepsilon}$ is a $(1+\varepsilon)\rho$-matching and
$\sigma_{\varepsilon}(X)\subseteq\operatorname{b}(y,3r_{X})$ then

8: if
$d_{M}(\sigma_{\varepsilon},\sigma^{\prime}_{\varepsilon})\geqslant\varepsilon
r_{X}/(2\rho^{2})$ for all $\sigma^{\prime}_{\varepsilon}\in
M(y,X,\varepsilon,r_{X})$ then

9: insert $\sigma_{\varepsilon}$ into $M(y,X,\varepsilon,r_{X})$

We first prove that this algorithm is correct. So we must prove that Property
17 holds for each set $M(y,X,\varepsilon,r_{X})$. Property 17a follows from
the fact that we only insert $\sigma_{\varepsilon}$ if it is a
$(1+\varepsilon)\rho$-matching and if
$\sigma_{\varepsilon}(X)\subseteq\operatorname{b}(y,3r_{X})$. Property 17b
follows from fact that we only insert $\sigma_{\varepsilon}$ if it is at
distance at least $\varepsilon r_{X}/(2\rho^{2})$ from all the previously
inserted matchings.

We now prove that Property 17c holds. So given a $\rho$-matching
$\bar{\sigma}:X\to Y$, we want to prove that there exist $\bar{y}\in
Y_{r_{X}}$ and $\bar{\sigma}_{\varepsilon}\in M(\bar{y},X,\varepsilon,r_{X})$
such that $d_{M}(\bar{\sigma},\bar{\sigma}_{\varepsilon})\leqslant\varepsilon
r_{X}/\rho^{2}$. Let $\bar{\sigma}^{P}$ be the restriction of $\bar{\sigma}$
to $P$. In other words, $\bar{\sigma}^{P}:P\to Y$ is defined by
$\bar{\sigma}^{P}(x)=\bar{\sigma}(x)$ for all $x\in P$. Similarly, let
$\bar{\sigma}^{Q}$ be the restriction of $\bar{\sigma}$ to $Q$. Then
$\bar{\sigma}^{P}$ and $\bar{\sigma}^{Q}$ are $\rho$-matchings. As Property
17c holds for $L(P,\beta,r_{X})$, there exist $\bar{y}^{P}\in Y_{r_{X}}$ and
$\bar{\sigma}^{P}_{\beta}\in M(\bar{y}^{P},P,\beta,r_{X})$ such that
$d_{Y}(\bar{\sigma}^{P},\bar{\sigma}^{P}_{\beta})\leqslant\beta
r_{X}/\rho^{2}$. Similarly, there exist $\bar{y}^{Q}\in Y_{r_{X}}$ and
$\bar{\sigma}^{Q}_{\beta}\in M(\bar{y}^{Q},Q,\beta,r_{X})$ such that
$d_{Y}(\bar{\sigma}^{Q},\bar{\sigma}^{Q}_{\beta})\leqslant\beta
r_{X}/\rho^{2}$.

The matching
$\bar{\sigma}_{\varepsilon}=\bar{\sigma}^{P}_{\beta}\cdot\bar{\sigma}^{Q}_{\beta}$
satisfies $d_{M}(\bar{\sigma},\bar{\sigma}_{\varepsilon})\leqslant\beta
r_{X}/\rho^{2}<\varepsilon r_{X}/\rho^{2}$. The lemma below shows that it is a
$(1+\varepsilon)\rho$-matching.

######  Lemma 20.

The matching $\bar{\sigma}_{\varepsilon}$ is a $(1+\varepsilon)\rho$-matching
from $X$ to $Y$.

###### Proof.

Let $x,x^{\prime}\in X$. If $(x,x^{\prime})\in P^{2}$ or $(x,x^{\prime})\in
Q^{2}$, then

$\frac{1}{(1+\beta)\rho}d_{X}(x,x^{\prime})\leqslant
d_{Y}(\bar{\sigma}_{\varepsilon}(x),\bar{\sigma}_{\varepsilon}(x^{\prime}))\leqslant(1+\beta)\rho
d_{X}(x,x^{\prime})$

since $\bar{\sigma}^{P}_{\beta}$ and $\bar{\sigma}^{Q}_{\beta}$ are
$(1+\beta)\rho$-matchings by Property 17a. As $\beta<\varepsilon$, it follows
that

$\frac{1}{(1+\varepsilon)\rho}d_{X}(x,x^{\prime})\leqslant
d_{Y}(\bar{\sigma}_{\varepsilon}(x),\bar{\sigma}_{\varepsilon}(x^{\prime}))\leqslant(1+\varepsilon)\rho
d_{X}(x,x^{\prime})$

which is the desired inequality.

So we may now assume, without loss of generality, that $x\in P$ and
$x^{\prime}\in Q$, and thus $d_{X}(x,x^{\prime})\geqslant d_{X}(P,Q)$. As $X$
was partitioned into $P$ and $Q$ using Lemma 16, we have
$\operatorname{diam}(X)\leqslant(k-1)\cdot d_{X}(P,Q)$. By the definition of
$r_{X}$, we also have $\rho\operatorname{diam}(X)>r_{X}/2$, and thus
$(k-1)d_{X}(x,x^{\prime})>r_{X}/(2\rho)$. By the triangle inequality, we have
$d_{Y}(\bar{\sigma}_{\varepsilon}(x),\bar{\sigma}_{\varepsilon}(x^{\prime}))\leqslant
d_{Y}(\bar{\sigma}_{\varepsilon}(x),\bar{\sigma}(x))+d_{Y}(\bar{\sigma}(x),\bar{\sigma}(x^{\prime}))+d_{Y}(\bar{\sigma}(x^{\prime}),\bar{\sigma}_{\varepsilon}(x^{\prime})),$
and thus

$\displaystyle
d_{Y}(\bar{\sigma}_{\varepsilon}(x),\bar{\sigma}_{\varepsilon}(x^{\prime}))$
$\displaystyle\leqslant
d_{Y}(\bar{\sigma}(x),\bar{\sigma}(x^{\prime}))+\frac{2\beta r_{X}}{\rho^{2}}$
because $d_{M}(\bar{\sigma},\bar{\sigma}_{\varepsilon})\leqslant\frac{\beta
r_{X}}{\rho^{2}}$ $\displaystyle\leqslant\rho d_{X}(x,x^{\prime})+\frac{2\beta
r_{X}}{\rho^{2}}$ since $\bar{\sigma}$ is a $\rho$-matching
$\displaystyle\leqslant\rho
d_{X}(x,x^{\prime})+\frac{4\beta(k-1)}{\rho}d_{X}(x,x^{\prime})$
$\displaystyle\leqslant\rho
d_{X}(x,x^{\prime})+\frac{\varepsilon}{2\rho}d_{X}(x,x^{\prime})$ because
$\beta=\frac{\varepsilon}{8k-8}$
$\displaystyle\leqslant\left(1+\frac{\varepsilon}{2}\right)\rho
d_{X}(x,x^{\prime}).$

We just proved the right-hand side of Inequality (1). For the left-hand side,
we use the same arguments as above, and obtain

$\displaystyle
d_{Y}(\bar{\sigma}_{\varepsilon}(x),\bar{\sigma}_{\varepsilon}(x^{\prime}))$
$\displaystyle\geqslant
d_{Y}(\bar{\sigma}(x),\bar{\sigma}(x^{\prime}))-\frac{2\beta
r_{X}}{\rho^{2}}\geqslant\frac{1}{\rho}d_{X}(x,x^{\prime})-\frac{2\beta
r_{X}}{\rho^{2}}$
$\displaystyle\geqslant\frac{1}{\rho}d_{X}(x,x^{\prime})-\frac{4\beta(k-1)}{\rho}d_{X}(x,x^{\prime})$
$\displaystyle=\frac{1}{\rho}d_{X}(x,x^{\prime})-\frac{\varepsilon}{2\rho}d_{X}(x,x^{\prime})=\left(1-\frac{\varepsilon}{2}\right)\frac{1}{\rho}d_{X}(x,x^{\prime}).$

So in order to complete the proof, we only need to argue that
$1-\varepsilon/2\geqslant 1/(1+\varepsilon).$ This is equivalent to
$1+\varepsilon/2-\varepsilon^{2}/2\geqslant 1$, which follows from our
assumptions that $0<\varepsilon\leqslant 1$. ∎

As $\bar{\sigma}$ is a $\rho$-matching, we have
$\operatorname{diam}(\bar{\sigma}(X))\leqslant\rho\operatorname{diam}(X)$ and
thus $\operatorname{diam}(\bar{\sigma}(X))\leqslant r_{X}$. By Lemma 14, it
implies that there exists $\bar{y}\in Y_{r_{X}}$ such that
$\bar{\sigma}(X)\subseteq\operatorname{b}(\bar{y},2r_{X})$, and thus
$\bar{\sigma}^{P}(P)\subseteq\operatorname{b}(\bar{y},2r_{X})$. As
$d_{M}(\bar{\sigma}^{P},\bar{\sigma}^{P}_{\beta})\leqslant\beta
r_{X}/\rho^{2}\leqslant r_{X}/8$, we have
$\bar{\sigma}^{P}_{\beta}(P)\subseteq\operatorname{b}(\bar{y},17r_{X}/8)$. By
Property 17a, we also have
$\bar{\sigma}^{P}_{\beta}(P)\subseteq\operatorname{b}(\bar{y}^{P},3r_{X})$, so
the balls $\operatorname{b}(\bar{y},17r_{X}/8)$ and
$\operatorname{b}(\bar{y}^{P},3r_{X})$ intersect, which implies that
$d_{Y}(\bar{y},\bar{y}^{P})<6r_{X}$. The same proof shows that
$d_{Y}(\bar{y},\bar{y}^{Q})<6r_{X}$.

Therefore, our algorithm considers
$\sigma_{\varepsilon}=\sigma_{\beta}^{P}\cdot\sigma_{\beta}^{Q}$ as a
candidate solution for each $\sigma_{\beta}^{P}\in
M(\bar{y}^{P},P,\beta,r_{X})$ and each $\sigma_{\beta}^{Q}\in
M(\bar{y}^{Q},Q,\beta,r_{X})$. In particular, we must have considered the
matching
$\bar{\sigma}_{\varepsilon}=\bar{\sigma}^{P}_{\beta}\cdot\bar{\sigma}^{Q}_{\beta}$.
As $d_{M}(\bar{\sigma},\bar{\sigma}_{\varepsilon})<\varepsilon r_{X}/\rho^{2}$
and $\bar{\sigma}(X)\subseteq b(\bar{y},2r_{X})$, we have
$\bar{\sigma}_{\varepsilon}(X)\subseteq b(\bar{y},3r_{X})$, and thus we must
have attempted to insert $\bar{\sigma}_{\varepsilon}$ into
$M(\bar{y},X,\varepsilon,r_{X})$. If $\bar{\sigma}_{\varepsilon}$ was inserted
into $M(\bar{y},X,\varepsilon,r_{X})$, then we are done. Otherwise, it means
that there exists $\sigma_{\varepsilon}^{\prime}\in
M(\bar{y},X,\varepsilon,r_{X})$ such that
$d_{M}(\bar{\sigma}_{\varepsilon},\sigma_{\varepsilon}^{\prime})<\varepsilon
r_{X}/(2\rho^{2})$. Since
$d_{M}(\bar{\sigma},\bar{\sigma}_{\varepsilon})\leqslant\beta
r_{X}/\rho^{2}<\varepsilon r_{X}/(2\rho^{2})$, it follows that
$d_{M}(\bar{\sigma},\sigma_{\varepsilon}^{\prime})<\varepsilon
r_{X}/\rho^{2}$. In any case, it shows that Property 17c holds. So we obtain
the following result.

######  Lemma 21.

Suppose that $k\geqslant 2$ and $\beta=\varepsilon/(8k-8)$. Then we can
compute $L(X,\varepsilon,r_{X})$ from $L(P,\beta,r_{X})$ and
$L(Q,\beta,r_{X})$ in $O(nk(\rho^{2}/\beta)^{3k\delta})$ time.

###### Proof.

The discussion above shows that our algorithm is correct. We still need to
analyze its running time. Let $y\in Y_{r_{X}}$. As $d_{Y}(y,y^{P})\leqslant
6r_{X}$ and $d_{Y}(y,y^{Q})\leqslant 6r_{X}$, the nodes $y^{P}$ and $y^{Q}$
are connected to $y$ by horizontal edges. As the nodes of the navigating net
have constant degree, it implies that there are $O(1)$ pairs $(y^{P},y^{Q})$
to consider. By Lemma 18, there are $O((\rho^{2}/\beta)^{k\delta})$ matchings
in $M(y^{P},P,\beta,r_{X})$ and $M(y^{Q},Q,\beta,r_{X})$. Therefore, we
consider $O((\rho^{2}/\beta)^{2k\delta})$ matchings $\sigma_{\varepsilon}$
when constructing $M(y,X,\varepsilon,r_{X})$. Each of these matchings is then
compared with previously inserted matchings. By Lemma 18, we have
$|M(y,X,\varepsilon,r_{X})|=O((\rho^{2}/\varepsilon)^{k\delta})=O((\rho^{2}/\beta)^{k\delta})$.
We can compute the distance between two matchings in $O(k)$ time, so it takes
$O(k(\rho^{2}/\beta)^{3k\delta})$ time to compute $M(y,X,\varepsilon,r_{X})$.
As there are at most $n$ nodes $y\in Y_{r_{X}}$, the overall time bound is
$O(nk(\rho^{2}/\beta)^{3k\delta})$. ∎

###  4.5 Putting everything together

We can now describe our algorithm for the $(\rho,\varepsilon)$-distortion
problem. We first compute the navigating net in $O(n\log\Phi(Y))$ time. If
$r_{X}<r_{\min}$, then $\rho\operatorname{diam}X<\operatorname{dmin}(Y)$, and
thus there is no $\rho$-matching. Otherwise, we recursively compute
$L(X,\varepsilon,r_{X})$, as described below. If $L(X,\varepsilon,r_{X})$
contains at least one matching, then we return one of them, which by Property
17a is a $(1+\varepsilon)\rho$-matching. By Property 17c, if there is a
solution to the $\rho$-distortion problem, then at least one
$(1+\varepsilon)\rho$-matching must be recorded in $L(X,\varepsilon,r_{X})$,
which shows that this algorithm indeed solves the
$(\rho,\varepsilon)$-distortion problem.

We now explain how to compute $L(X,\varepsilon,r_{X})$ recursively. The base
case is $k=1$ and $r_{X}=r_{\min}$. As explained above, it can be done in
$O(n)$ time by recording a trivial matching at each leaf node.

When $k\geqslant 2$, we split $X$ into $P$ and $Q$ as described in Lemma 16.
Then we compute recursively $L(P,\beta,r_{P})$ and $L(Q,\beta,r_{Q})$ where
$\beta=\varepsilon/(8k-8)$ and $r_{P}$ and $r_{Q}$ are the smallest scales at
least as large as $\rho\operatorname{diam}(P)$ and
$\rho\operatorname{diam}(Q)$, respectively. We compute $L(P,\beta,r_{X})$ and
$L(Q,\beta,r_{X})$ using Lemma 19, which takes time
$O(nk(\rho^{2}/\beta)^{2k\delta}+n\log\Phi(Y))$. Then we obtain
$L(X,\varepsilon,r_{X})$ by Lemma 21 in $O(nk(\rho^{2}/\beta)^{3k\delta})$
time.

So the running time $T(n,k,\rho,\varepsilon)$ of our algorithm satisfies the
relation

$T(n,k,\rho,\varepsilon)=T(n,k_{1},\rho,\beta)+T(n,k_{2},\rho,\beta)+O(nk(\rho^{2}/\beta)^{3k\delta}+n\log\Phi(Y)),$

where $k_{1}+k_{2}=k$, $1\leqslant k_{1}\leqslant k-1$ and
$\beta=\varepsilon/(8k-8)$. This expression expands into a sum of $k-1$ terms
$O(nk(\rho^{2}(8k)^{k}/\varepsilon)^{3k\delta}+n\log\Phi(Y))$ and $k$ terms
$O(n)$ for the base cases. Hence we have
$T(n,k,\rho,\varepsilon)=2^{O(k^{2}\log
k)}n(\rho^{2}/\varepsilon)^{3k\delta}+O(kn\log\Phi(Y))$. We just proved the
following:

######  Theorem 22.

The $(\rho,\varepsilon)$-distortion problem can be solved in $2^{O(k^{2}\log
k)}(\rho^{2}/\varepsilon)^{3k\delta}n+O(kn\log\Phi(Y))$ time.

Our algorithm may also allow us to find more than one matching. If we return
all the matchings stored in $L(X,\varepsilon,r_{X})$, we obtain a collection
of $(1+\varepsilon)\rho$-matching that approximates all the solutions: By
Property 17c, for each $\rho$-matching $\sigma$, one of the
$(1+\varepsilon)\rho$-matching $\sigma_{\varepsilon}$ that we return
approximates $\sigma$ in the sense that
$d_{M}(\sigma,\sigma_{\varepsilon})\leqslant\varepsilon
r_{X}/\rho^{2}<2\varepsilon\operatorname{diam}(X)/\rho$.

##  5 Improved Algorithm

In this section, we show how the approximation algorithm presented in Section
4 can be improved using an efficient approximate nearest-neighbor (ANN) data
structure. More precisely, we will use this data structure to answer fixed-
radius near neighbor queries efficiently.

First we briefly introduce $(3/2)$-approximate nearest neighbor searching. Let
$(S,d_{S})$ be a doubling space, with doubling dimension $\delta_{S}$. Let $T$
be an $n$-point subset of $S$. Given a query point $q\in S$, its nearest
neighbor in $T$ is a point $t_{q}\in T$ such that $d_{S}(q,t_{q})$ is minimum.
A $(3/2)$-approximate nearest neighbor (ANN) is a point $\hat{q}\in T$ such
that $d_{S}(q,\hat{q})\leqslant(3/2)d_{S}(q,t_{q})$. The data structure by
Cole and Gottlieb [9] returns a $(3/2)$-ANN in time $2^{O(\delta_{S})}\log n$.
This data structure is dynamic: It allows the insertion or deletion of a point
in $2^{O(\delta_{S})}\log n$ time.

###  5.1 Fixed-radius near neighbor searching

We now consider a variation on the problem above. Given a radius $R$ and a
query point $q$, a fixed-radius near neighbor query returns all the points in
$T$ that are at distance at most $R$ from $q$. We now show how to efficiently
answer these queries in a doubling space.

######  Lemma 23.

Suppose that $\operatorname{dmin}(S)\geqslant R/\lambda$ for some constant
$\lambda>0$. Then Gottlieb and Cole’s ANN data structure [9] allows us to
answer fixed-radius near neighbor queries in $2^{O(\delta_{S})}\log n$ time.
Each such query returns $2^{O(\delta_{S})}$ points. This data structure also
supports insertions and deletions in $2^{O(\delta_{S})}\log n$ time.

We now show how to answer a query. Let $q$ be a query point, so we want to
return all the points in $T\cap\operatorname{b}(q,R)$. We repeatedly perform
(3/2)-ANN queries, as long as the distance $d_{S}(q,\hat{q})$ from $q$ to its
approximate near neighbor $\hat{q}$ is at most $3R/2$. If
$d_{S}(q,\hat{q})\leqslant R$, then we report $\hat{q}$. In any case, we
delete $\hat{q}$ from $T$ and insert it into a set $U$. In the end, we find a
point $\hat{q}$ such that $d_{S}(q,\hat{q})>3R/2$, and we complete the
procedure by inserting all the points of $U$ back into $T$.

As $U\subseteq\operatorname{b}(q,3R/2)$ we have
$\operatorname{diam}(U)\leqslant 3R$. Since $U\subseteq S$ we also have
$\operatorname{dmin}(S)\geqslant R/\lambda$. Therefore, $\Phi(U)\leqslant
3\lambda$, and since $U$ has doubling dimension $\delta_{S}$, we have
$|U|=2^{O(\delta_{S})}$ by Lemma 1. Therefore, we answer one fixed-radius near
neighbor query by performing $2^{O(\delta_{S})}$ ANN queries, insertions and
deletions. So a fixed-radius near neighbor query is answered in
$2^{O(\delta_{S})}\log n$ time.

###  5.2 Faster pruning

The bottleneck of our approximation algorithm is when we check whether a
matching $\sigma_{\varepsilon}$ is at distance at least $\varepsilon
r_{X}/(2\rho^{2})$ from each matching $\sigma_{\varepsilon}^{\prime}\in
M(y,X,\varepsilon,r_{X})$. (See Algorithm 2.) As
$|M(y,X,\varepsilon,r_{X})|=O((\rho^{2}/\varepsilon)^{k\delta})$, it takes
time $O(k(\rho^{2}/\varepsilon)^{k\delta})$ for each matching
$\sigma_{\varepsilon}$ using brute force.

Instead of checking by brute force, we record each set
$M(y,X,\varepsilon,r_{X})$ in a $(3/2)$-ANN data structure. Remember that we
identify the space of matchings from $X$ to $Y$ with $Y^{k}$, and thus this
space has doubling dimension $k\delta$ by Corollary 3. By Property 17b, we
have $\operatorname{dmin}(M(y,X,\varepsilon,r_{X}))\geqslant\varepsilon
r_{X}/(2\rho^{2})$. By Lemma 23, this data structure allows us to perform
insertions and fixed-radius near neighbor queries in $2^{O(k\delta)}\log n$
time with radius $R=\varepsilon r_{X}/(2\rho^{2})$. As we assume that
$\delta=O(1)$, this is $2^{O(k)}\log n$ time. We only need one such query to
determine whether $\sigma_{\varepsilon}$ is at distance at least $\varepsilon
r_{X}/(2\rho^{2})$ from each matching $\sigma_{\varepsilon}^{\prime}\in
M(y,X,\varepsilon,r_{X})$.

In summary, we replace the factor $O(k(\rho^{2}/\varepsilon)^{k\delta})$ in
the time bound from Lemma 21 with a factor $2^{O(k)}\log n$. So the running
time of our approximation algorithm is given by the recurrence relation (see
Theorem 22):

$T(n,k,\rho,\varepsilon)=T(n,k_{1},\rho,\beta)+T(n,k_{2},\rho,\beta)+2^{O(k)}(\rho^{2}/\beta)^{2k\delta}n\log
n+O(n\log\Phi(Y)),$

We thus obtain the following time bound.

######  Lemma 24.

The $(\rho,\varepsilon)$-distortion problem can be solved in $2^{O(k^{2}\log
k)}(\rho^{2}/\varepsilon)^{2k\delta}n\log n+O(kn\log\Phi(Y))$ time.

###  5.3 Removing the dependency on the spread

We now show how to remove the $O(kn\log\Phi(Y))$ term in the above time bound.
This term comes from the construction and traversal of the navigating net. But
we do not need to compute more than $k$ levels of the navigating nets. More
precisely, we need to compute it at scale $r_{X}$, and at each scale
$r_{X^{\prime}}$ where $X^{\prime}$ is obtained by recursively splitting $X$
into two subsets according to Lemma 16. In order to perform our bottom-up
construction, we will also need to be able to find, for two such scales
$r<r^{\prime}$, an ancestor in $Y_{r^{\prime}}$ for any node $y\in Y_{r}$. We
now explain how to achieve this using the ANN data structure.

We first show how to compute $Y_{r}$ in $O(n\log n)$ time. We construct
$Y_{r}$ incrementally: For each point $y\in Y$, we check whether it is at
distance at least $r$ from each previously inserted point. It can be checked
in $O(\log n)$ time by performing a fixed-radius near neighbor query. If $y$
is indeed at distance at least $r$ from each previously inserted point, we
insert $y$ into $Y_{r}$. In order to perform this test efficiently, we record
$Y_{r}$ in an ANN data structure. By Lemma 23, as
$\operatorname{dmin}(Y_{r})\geqslant r$, it allows us to answer fixed-radius
near neighbor queries in $O(\log n)$ time with radius $R=r$. Insertions are
also performed in $O(\log n)$ time, so we obtain the following:

######  Lemma 25.

Let $Y$ be an $n$-point metric space of constant doubling dimension. We can
compute an $r$-net of $Y$ in $O(n\log n)$ time.

Our algorithm for the $(\rho,\varepsilon)$-distortion problem requires us to
know the horizontal edges of each $r$-net we consider. They can also be
computed in $O(n\log n)$ time using the same approach as the lemma above,
except that the radius $R$ is now $6r$.

######  Lemma 26.

We can compute all the horizontal edges of $Y_{r}$ in $O(n\log n)$ time.

Finally, at scales $r<r^{\prime}$, we need to be able to find the ancestor in
$Y_{r^{\prime}}$ for each node $y\in Y_{r}$, in order to run our bottom-up
construction of $L(W,\beta,r^{\prime})$ from $L(W,\beta,r)$. As we no longer
have vertical edges, we relax the definition of ancestors as follows: The
ancestor of $y\in Y_{r}$ is a point $y^{\prime}\in Y_{r^{\prime}}$ such that
$d_{Y}(y,y^{\prime})\leqslant r^{\prime}$. (There could be several such
points; we can pick any one of them to be the ancestor of $y$.) This property
suffices for our algorithm to work as described in Section 4.

In order to compute the ancestors, we construct an ANN data structure for
$Y_{r^{\prime}}$. By Lemma 23, as
$\operatorname{dmin}(Y_{r^{\prime}})\geqslant r^{\prime}$, we can construct
the ancestor of any point in $O(\log n)$ time by performing a fixed radius
near-neighbor query with $R=r^{\prime}$.

######  Lemma 27.

Given $Y_{r}$ and $Y_{r^{\prime}}$ such that $r^{\prime}>r$, we can compute
$O(n\log n)$ time the ancestors in $Y_{r^{\prime}}$ of all the points in
$Y_{r}$.

This shows that we can run the algorithm from Section 4 without computing the
whole navigating net, hence the $O(nk\log\Phi(Y))$ term in Lemma 24 is
replaced with $O(kn\log n)$, as we only compute $Y_{r}$ for $k$ different
values of $r$. In summary, we obtain the following result.

######  Theorem 28.

The $(\rho,\varepsilon)$-distortion problem can be solved in $2^{O(k^{2}\log
k)}(\rho^{2}/\varepsilon)^{2k\delta}n\log n$ time.

##  6 Minimum distortion

In this section, we show how to extend our results on the $\rho$-distortion
problem (Sections 3, 4 and 5) to the minimum distortion problem.

###  6.1 Hardness result

We show how to modify the construction from Section 3 so as to obtain a
reduction to the minimum distortion problem. Let $\lambda$ be a large enough
number, say $\lambda=5m\rho^{2}2^{k}$. We add an extra point $x_{k+1}$ to $X$
such that $d_{X}(x_{i},x_{k+1})=\lambda$ for all $i\leqslant k$. We also add a
point $p_{k+1}$ to $Y$ that is at distance $\lambda$ from every other point.
If there is a $k$-clique in $G$, then the mapping with minimum distortion is
the same as in Section 3, where we add $\sigma_{x_{k+1}}=y_{k+1}$. It
satisfies $\operatorname{expansion}(\sigma)=1$ and
$\operatorname{expansion}(\sigma^{-1})=\rho$, hence
$\operatorname{dist}(X,Y)=\operatorname{dist}(\sigma)=\rho$.

Conversely, suppose that $\operatorname{dist}(X,Y)=\rho$. We first argue that
we must have $\sigma(x_{k+1})=p_{k+1}$. If it were not the case, then
$d_{Y}(\sigma(x_{1}),\sigma(x_{k+1}))\leqslant 2^{k}$ so
$\operatorname{expansion}(\sigma^{-1})\geqslant 5m\rho^{2}$. As
$d_{X}(x_{1},x_{2})=4\rho$ and $d_{Y}(\sigma(x_{1}),\sigma(x_{2}))\geqslant
1/m$, we also have $\operatorname{expansion}(\sigma)\geqslant 1/(4\rho m)$ and
thus $\operatorname{dist}(x,y)\geqslant 5\rho/4$, a contradiction.

Therefore, we must have $\sigma(x_{k+1})=y_{k+1}$. It implies that
$\operatorname{expansion}(\sigma)\geqslant 1$. Then the only way to obtain
$\operatorname{dist}(\sigma)\leqslant\rho$ is to use the same construction as
in Section 3, where we match each point $x_{i}$ with a point in $R_{i}$ for
$3\leqslant i\leqslant k$. The corresponding vertices in the $k$-clique
instance form a clique.

######  Corollary 29.

Suppose that $X$ and $Y$ have doubling dimension $\log_{2}3$. For any
$\rho\geqslant 1$, the problem of deciding whether the minimum distortion
$\operatorname{dist}(X,Y)$ is equal to $\rho$ is NP-hard, and is $W[1]$-hard
when parameterized by $k$. It cannot be solved in time $f(k)\cdot n^{o(k)}$
for any computable function $f$, unless ETH is false.

###  6.2 A first algorithm

We now present a decision algorithm. Given two positive reals $e$,
$e^{\prime}$, we would like to know whether there is a mapping $\sigma:X\to Y$
such that $\operatorname{expansion}(\sigma)\leqslant e$ and
$\operatorname{expansion}(\sigma^{-1})\leqslant e^{\prime}$. As observed by
Kenyon et al. [21], this is true if and only if there is mapping
$\sigma^{\prime}:X\to Y$ with
$\operatorname{expansion}(\sigma^{\prime})\leqslant\sqrt{ee^{\prime}}$ and
$\operatorname{expansion}(\sigma^{\prime-1})\leqslant\sqrt{ee^{\prime}}$,
where we use the metric $d^{\prime}_{Y}=\sqrt{e^{\prime}/e}\cdot d_{Y}$ on
$Y$, and we still use $d_{X}$ on $X$. In other words, we are solving the
$\rho$-distortion problem with $\rho=\sqrt{ee^{\prime}}$ on $(X,d_{X})$ and
$(Y,d^{\prime}_{Y})$. So we obtain an approximate answer using our algorithm
from Theorem 28:

######  Lemma 30.

Given $e,e^{\prime}>0$ and $0<\varepsilon\leqslant 1$, our decision algorithm
returns a positive answer if there exists a mapping $\sigma:X\to Y$ such that
$\operatorname{expansion}(\sigma)\leqslant e$ and
$\operatorname{expansion}(\sigma^{-1})\leqslant e^{\prime}$, and returns a
negative answer if there is no mapping $\sigma:X\to Y$ such that
$\operatorname{expansion}(\sigma)\leqslant(1+\varepsilon)e$ and
$\operatorname{expansion}(\sigma^{-1})\leqslant(1+\varepsilon)e^{\prime}$. Its
running time is $2^{O(k^{2}\log k)}(ee^{\prime}/\varepsilon)^{2k\delta}n\log
n.$

Still following the approach by Kenyon et al., we can turn this decision
algorithm into an optimization algorithm by observing that
$\operatorname{expansion}(\sigma)$ and $\operatorname{expansion}(\sigma^{-1})$
are ratios of distances between points of $X$ and $Y$, and thus there are
$O(k^{2}n^{2})$ candidate values for each. So for any distinct
$x,x^{\prime}\in X$ and $y,y^{\prime}\in Y$, we run the algorithm from Lemma
30 using $e=d_{Y}(y,y^{\prime})/d_{X}(x,x^{\prime})$ and
$e^{\prime}=d_{X}(x,x^{\prime})/d_{Y}(y,y^{\prime})$ and with a relative error
ratio $\varepsilon^{\prime}=\varepsilon/3$. More precisely, we consider these
pairs $(e,e^{\prime})$ by increasing value of the product $ee^{\prime}$, and
we stop as soon as we obtain a positive answer. At this point, there must be a
matching $\sigma:X\to Y$ such that
$\operatorname{expansion}(\sigma)\leqslant(1+\varepsilon/3)e$ and
$\operatorname{expansion}(\sigma^{-1})\leqslant(1+\varepsilon/3)e^{\prime}$,
and since $0<\varepsilon\leqslant 1$, we have
$\operatorname{dist}(X,Y)\leqslant(1+\varepsilon/3)^{2}\leqslant(1+\varepsilon)ee^{\prime}$.
On the other hand, we must have $ee^{\prime}\leqslant\operatorname{dist}(X,Y)$
as otherwise, the algorithm would have halted earlier. So taking
$\Delta_{\varepsilon}=ee^{\prime}$, we obtain the following:

######  Lemma 31.

Given $0<\varepsilon\leqslant 1$, we can compute a number
$\Delta_{\varepsilon}$ such that
$\operatorname{dist}(X,Y)\leqslant\Delta_{\varepsilon}\leqslant(1+\varepsilon)\operatorname{dist}(X,Y)$
in time $2^{O(k^{2}\log
k)}(\operatorname{dist}(X,Y)/\varepsilon)^{2k\delta}k^{4}n^{5}\log n$.

###  6.3 Improved algorithm

We now describe a faster version of this optimization algorithm. To this end,
we need to introduce well-separated pairs decompositions (WSPD) [19]. Given a
metric space $(S,d_{S})$ and two subsets $A,B\subseteq S$, we say that $A$ and
$B$ are $(1/\varepsilon)$-separated if
$\operatorname{diam}(A)\leqslant\varepsilon d_{S}(A,B)$ and
$\operatorname{diam}(B)\leqslant\varepsilon d_{S}(A,B)$. As a result, for any
$a,a^{\prime}\in A$ and $b,b^{\prime}\in B$, we have

$d_{S}(a,b)/(1+2\varepsilon)\leqslant
d_{S}(a^{\prime},b^{\prime})\leqslant(1+2\varepsilon)d_{S}(a,b).$ (4)

A WSPD of $S$ of size $m$ is a collection of pairs of subsets
$\\{A_{1},B_{1}\\}$, $\\{A_{2},B_{2}\\}$, $\dots,\\{A_{m},B_{m}\\}$ of $S$
such that $A_{i}$ and $B_{i}$ are $(1/\varepsilon)$-separated for all $i$, and
for any two distinct points $a,b\in S$, there is a pair $A_{i},B_{i}$ such
that $(a,b)\in A_{i}\times B_{i}$ or $(a,b)\in B_{i}\times A_{i}$.

Callahan and Kosaraju [5] showed that a WSPD of size $n/\varepsilon^{O(1)}$
can be implicitly computed in $O(n\log n)+n/\varepsilon^{O(1)}$ time when $S$
is a set of $n$ points in a fixed-dimensional Euclidean space. In particular,
we can obtain a pair of points $(a_{i},b_{i})\in A_{i}\times B_{i}$ for all
$i$ within this time bound. This result was generalized to metrics of fixed
doubling dimension by Talwar [26], but with size and construction time bounds
larger by a factor $O(\log\Phi(S))$. Still for doubling metrics, Har-Peled and
Mendel gave a randomized algorithm that constructs in $O(n\log
n)+n/\varepsilon^{O(1)}$ time a WSPD of size $n/\varepsilon^{O(1)}$ [20],
which can be derandomized using Cole and Gottlieb’s data structure [9]. So if
we let $\ell_{i}=d_{S}(a_{i},b_{i})$, we obtain the following.

######  Lemma 32.

We can compute in time $O(n\log n)+n/\varepsilon^{O(1)}$ a set of lengths
$\\{\ell_{1},\dots,\ell_{m}\\}$ such that $m=n/\varepsilon^{O(1)}$ and for
every pair of distinct points $a,b\in S$, there exists $i$ such that
$\ell_{i}/(1+\varepsilon)\leqslant
d_{S}(a,b)\leqslant(1+\varepsilon)\ell_{i}$.

We can now present an approximate decision algorithm for the minimum
distortion problem. Given a distortion $\Delta$, it returns a positive answer
if $\operatorname{dist}(X,Y)\leqslant\Delta$, and it returns a negative answer
if $\operatorname{dist}(X,Y)>(1+\varepsilon)\Delta$. We first compute using
Lemma 32 a set $\\{\ell_{1},\dots,\ell_{m}\\}$ such that
$m=n/\varepsilon^{O(1)}$ and for all $y,y^{\prime}\in Y$, $y\neq y^{\prime}$,
we have $\ell_{i}/(1-\varepsilon/6)\leqslant
d_{Y}(y,y^{\prime})\leqslant(1+\varepsilon/6)\ell_{i}$ for some $i$. Then for
each $i\in\\{1,\dots,m\\}$, and for each pair of distinct points
$x,x^{\prime}\in X$, we run the decision algorithm of Lemma 30 with expansions
$e=(1+\varepsilon/6)\ell_{i}/d_{X}(x,x^{\prime})$ and
$e^{\prime}=(1+\varepsilon/6)\Delta d_{X}(x,x^{\prime})/\ell_{i}$, and with
relative error ratio $\varepsilon/6$. If one of the answers to a calls to the
algorithm of Lemma 30 is positive, we return a positive answer, and otherwise
we return a negative answer.

We now prove that this algorithm is correct. Suppose that
$\operatorname{dist}(X,Y)\leqslant\Delta$. Then there exists a mapping
$\sigma:X\to Y$ such that
$\operatorname{expansion}(\sigma)\operatorname{expansion}(\sigma^{-1})\leqslant\Delta$.
We must have
$\operatorname{expansion}(\sigma)=d_{Y}(y,y^{\prime})/d_{X}(x,x^{\prime})$ for
some $x,x^{\prime}\in X$ and $y,y^{\prime}\in Y$. So there is $i$ such that

$\frac{1}{1+\varepsilon/6}\cdot\frac{\ell_{i}}{d_{X}(x,x^{\prime})}\leqslant\operatorname{expansion}(\sigma)\leqslant(1+\varepsilon/6)\frac{\ell_{i}}{d_{X}(x,x^{\prime})}$

and thus

$\operatorname{expansion}(\sigma^{-1})\leqslant\frac{\Delta}{\operatorname{expansion}(\sigma)}\leqslant(1+\varepsilon/6)\Delta\frac{d_{X}(x,x^{\prime})}{\ell_{i}}.$

So we have $\operatorname{expansion}(\sigma)\leqslant e$ and
$\operatorname{expansion}(\sigma^{-1})\leqslant e^{\prime}$, and thus the
corresponding call to the decision algorithm of Lemma 30 must have returned a
positive answer.

Now suppose that our algorithm returns a positive answer. So there must be
$x,x^{\prime}\in X$, $i\in\\{1,\dots,m\\}$ and a mapping $\sigma:X\to Y$ such
that $\operatorname{expansion}(\sigma)\leqslant(1+\varepsilon/6)e$ and
$\operatorname{expansion}(\sigma^{-1})\leqslant(1+\varepsilon/6)e$, where
$e=(1+\varepsilon/6)\ell_{i}/d_{X}(x,x^{\prime})$ and
$e^{\prime}=(1+\varepsilon/6)\Delta d_{X}(x,x^{\prime})/\ell_{i}$. It follows
that $\operatorname{dist}(X,Y)\leqslant(1+\varepsilon/6)^{2}ee^{\prime}$, and
thus $\operatorname{dist}(X,Y)\leqslant(1+\varepsilon/6)^{4}\Delta$. As
$0<\varepsilon\leqslant 1$, it implies that
$\operatorname{dist}(X,Y)<(1+\varepsilon)\Delta$.

So we have proved that this decision algorithm is correct. As we run the
algorithm of Lemma 30 $k^{2}n/\varepsilon^{O(1)}$ times, we obtain the
following.

######  Lemma 33.

Given $0<\varepsilon\leqslant 1<\Delta$, the approximate decision algorithm
above returns a positive answer if $\operatorname{dist}(X,Y)\leqslant\Delta$,
and a negative answer if
$\operatorname{dist}(X,Y)\geqslant(1+\varepsilon)\Delta$. Its running time is
$2^{O(k^{2}\log k)}(\Delta^{2k\delta}/\varepsilon^{2k\delta+O(1)})n^{2}\log
n.$

We can turn this decision algorithm into an approximation algorithm for the
minimum distortion using exponential search. So we run this algorithm with
$\Delta=1,2,4,\dots$ until it returns a positive answer for $\Delta=2^{i}$. At
this point we know that $2^{i-1}<\operatorname{dist}(X,Y)<2^{i+1}$. Then we
perform binary search for $\Delta$ in the interval $[2^{i-1},2^{i+1}]$, still
using the decision algorithm of Lemma 33. After $O(\log(1/\varepsilon))$
steps, it yields a $(1+\varepsilon)$-approximation of
$\operatorname{dist}(X,Y)$. So we obtain the following result.

######  Theorem 34.

Given $0<\varepsilon\leqslant 1$, we can compute a number
$\Delta_{\varepsilon}$ such that
$\operatorname{dist}(X,Y)\leqslant\Delta_{\varepsilon}\leqslant(1+\varepsilon)\operatorname{dist}(X,Y)$
in time $2^{O(k^{2}\log
k)}(\operatorname{dist}(X,Y)^{2k\delta}/\varepsilon^{2k\delta+O(1)})n^{2}\log
n$.

## References

  * [1] Helmut Alt and Leonidas J. Guibas.  Discrete geometric shapes: Matching, interpolation, and approximation.  In Handbook of Computational Geometry, pages 121–153. Elsevier, 2000.
  * [2] Helmut Alt, Kurt Mehlhorn, Hubert Wagener, and Emo Welzl.  Congruence, similarity, and symmetries of geometric objects.  Discrete & Computational Geometry, 3:237–256, 1988.
  * [3] Sunil Arya, David M. Mount, Antoine Vigneron, and Jian Xia.  Space-time tradeoffs for proximity searching in doubling spaces.  In Proc. European Symposium on Algorithms, pages 112–123, 2008\.
  * [4] Glencora Borradaile, Hung Le, and Christian Wulff-Nilsen.  Greedy spanners are optimal in doubling metrics.  In Proc. Symposium on Discrete Algorithms, pages 2371–2379, 2019\.
  * [5] Paul B. Callahan and S. Rao Kosaraju.  A decomposition of multidimensional point sets with applications to k-nearest-neighbors and n-body potential fields.  J. ACM, 42(1):67–90, 1995.
  * [6] T.-H. Hubert Chan, Shuguang Hu, and Shaofeng H.-C. Jiang.  A PTAS for the steiner forest problem in doubling metrics.  SIAM J. Comput., 47(4):1705–1734, 2018.
  * [7] T.-H. Hubert Chan, Mingfei Li, Li Ning, and Shay Solomon.  New doubling spanners: Better and simpler.  SIAM J. Comput., 44(1):37–53, 2015.
  * [8] L. Paul Chew, Michael T. Goodrich, Daniel P. Huttenlocher, Klara Kedem, Jon M. Kleinberg, and Dina Kravets.  Geometric pattern matching under euclidean motion.  Comput. Geom., 7:113–124, 1997.
  * [9] Richard Cole and Lee-Ad Gottlieb.  Searching dynamic point sets in spaces with bounded doubling dimension.  In Proc. ACM Symposium on Theory of Computing, pages 574–583, 2006.
  * [10] Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, and Clifford Stein.  Introduction to Algorithms, 3rd Edition.  MIT Press, 2009.
  * [11] Marek Cygan, Fedor V. Fomin, Alexander Golovnev, Alexander S. Kulikov, Ivan Mihajlin, Jakub Pachocki, and Arkadiusz Socala.  Tight lower bounds on graph embedding problems.  J. ACM, 64(3):18:1–18:22, 2017.
  * [12] Marek Cygan, Fedor V. Fomin, Lukasz Kowalik, Daniel Lokshtanov, Dániel Marx, Marcin Pilipczuk, Michal Pilipczuk, and Saket Saurabh.  Parameterized Algorithms.  Springer, 2015.
  * [13] Hu Ding and Mingquan Ye.  On geometric alignment in low doubling dimension.  In Proc. AAAI Conference on Artificial Intelligence, pages 1460–1467, 2019.
  * [14] Michael R. Fellows, Fedor V. Fomin, Daniel Lokshtanov, Elena Losievskaja, Frances A. Rosamond, and Saket Saurabh.  Distortion is fixed parameter tractable.  In Proc. International Colloquium on Automata, Languages and Programming, pages 463–474, 2009.
  * [15] Zachary Friggstad, Mohsen Rezapour, and Mohammad R. Salavatipour.  Local search yields a PTAS for k-means in doubling metrics.  SIAM J. Comput., 48(2):452–480, 2019.
  * [16] Lee-Ad Gottlieb and Liam Roditty.  Improved algorithms for fully dynamic geometric spanners and geometric routing.  In Proc. ACM-SIAM Symposium on Discrete algorithms, pages 591–600, 2008.
  * [17] Anupam Gupta, Robert Krauthgamer, and James R. Lee.  Bounded geometries, fractals, and low-distortion embeddings.  In Proc. IEEE Symposium on Foundations of Computer Science, pages 534–543, 2003.
  * [18] A. Hall and C. Papadimitriou.  Approximating the distortion.  In Proc. APPROX-RANDOM, pages 111–122, 2005.
  * [19] Sariel Har-peled.  Geometric Approximation Algorithms.  American Mathematical Society, 2011.
  * [20] Sariel Har-Peled and Manor Mendel.  Fast construction of nets in low-dimensional metrics and their applications.  SIAM J. Comput., 35(5):1148–1184, 2006.
  * [21] Claire Kenyon, Yuval Rabani, and Alistair Sinclair.  Low distortion maps between point sets.  SIAM J. Comput., 39(4):1617–1636, 2009.
  * [22] Michael Kerber and Arnur Nigmetov.  Metric spaces with expensive distances.  CoRR, abs/1901.08805, 2019.
  * [23] Subhash Khot and Rishi Saket.  Hardness of embedding metric spaces of equal size.  In Proc. APPROX-RANDOM, pages 218–227, 2007.
  * [24] Robert Krauthgamer and James R. Lee.  Navigating nets: Simple algorithms for proximity search.  In Proc. ACM-SIAM Symposium on Discrete Algorithms, pages 798–807, 2004.
  * [25] Amir Nayyeri and Benjamin Raichel.  Reality distortion: Exact and approximate algorithms for embedding into the line.  In Proc. IEEE Symposium on Foundations of Computer Science, pages 729–747, 2015.
  * [26] Kunal Talwar.  Bypassing the embedding: algorithms for low dimensional metrics.  In Proc. ACM Symposium on Theory of Computing, pages 281–290, 2004.