IST Austria
[Klosterneuburg,
Austria]hana.kourimska@ist.ac.athttps://orcid.org/0000-0001-7841-0091 No
<EMAIL_ADDRESS>IST Austria
[Klosterneuburg, Austria]
Inria Sophia Antipolis, Université Côte d’Azur
[Sophia Antipolis, France]
m.h.m.j.wintraecken@gmail.comhttps://orcid.org/0000-0002-7472-2220Supported by
the European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie grant agreement No. 754411, and the Austrian science
fund (FWF) grant No. M-3073 This research has been supported by the European
Research Council (ERC), grant No. 788183, by the Wittgenstein Prize, Austrian
Science Fund (FWF), grant No. Z 342-N31, and by the DFG Collaborative Research
Center TRR 109, Austrian Science Fund (FWF), grant No. I 02979-N35.

###### Acknowledgements.

We are greatly indebted to Fred Chazal for sharing his insights. We further
thank Erin Chambers, Christopher Fillmore, and Elizabeth Stephenson for early
discussions and all members of the Edelsbrunner group (Institute of Science
and Technology Austria) and the Datashape team (Inria) for the atmosphere in
which this research was conducted. Hana Dal Poz Kouřimská, André Lieutier, and
Mathijs Wintraecken Theory of computation $\rightarrow$ Computational geometry

#  The medial axis of closed bounded sets is Lipschitz stable with respect to
the Hausdorff distance under ambient diffeomorphisms

Hana Dal Poz Kouřimská  André Lieutier  Mathijs Wintraecken

###### Abstract

We prove that the medial axis of closed sets is Hausdorff stable in the
following sense: Let ${\mathcal{S}}\subseteq\mathbb{R}^{d}$ be (fixed) closed
set (that contains a bounding sphere). Consider the space of $C^{1,1}$
diffeomorphisms of $\mathbb{R}^{d}$ to itself, which keep the bounding sphere
invariant. The map from this space of diffeomorphisms (endowed with some
Banach norm) to the space of closed subsets of $\mathbb{R}^{d}$ (endowed with
the Hausdorff distance), mapping a diffeomorphism $F$ to the closure of the
medial axis of $F({\mathcal{S}})$, is Lipschitz.

This extends a previous stability result of Chazal and Soufflet on the
stability of the medial axis of $C^{2}$ manifolds under $C^{2}$ ambient
diffeomorphisms.

###### keywords:

Medial axis, Hausdorff distance, Lipschitz continuity

##  1 Introduction

In [8], Federer introduced the _reach_ of a (closed) set
${\mathcal{S}}\subset\mathbb{R}^{d}$ as the infimum over all points in
${\mathcal{S}}$ of the distance from these points to the _medial axis_
$\operatorname{{ax}}({\mathcal{S}})$, the set of points in $\mathbb{R}^{d}$
for which the closest point in ${\mathcal{S}}$ is not unique. Although the
(complement of the) medial axis was studied by Federer, the term ‘medial axis’
was coined by Blum [5]. Federer also introduced the reach at a point
$p\in{\mathcal{S}}$ to be the distance from $p$ to the medial axis of
${\mathcal{S}}$. We now call this quantity the _local feature size_ [1] and
denote it by $\operatorname{lfs}(p)$. We illustrate the three notions in
Figure 1.

Figure 1: A set ${\mathcal{S}}$ of positive reach (in black/gray), the medial
axis $\operatorname{{ax}}({\mathcal{S}})$ of ${\mathcal{S}}$ (in green), and
the local feature size $\operatorname{lfs}(.)$ of a few points in
${\mathcal{S}}$. The reach is attained — for example — at point $p_{2}$. We
always assume that the set ${\mathcal{S}}$ as well as its medial axis are
bounded, and that the bounding sphere is contained in ${\mathcal{S}}$ — we
explain this technicality in Section 2. In light green we depict the points
lying on the closure of $\operatorname{{ax}}({\mathcal{S}})$ that are not part
of the medial axis itself.

Federer proved that the reach is stable under $C^{1,1}$ diffeomorphisms of the
ambient space. Here, a $C^{1,1}$ map is a $C^{1}$ map whose derivative is
Lipschitz, and a $C^{1,1}$ diffeomorphism is a $C^{1,1}$ bijective map whose
inverse is also $C^{1,1}$. Chazal and Soufflet [6] proved that the medial axis
is stable with respect to the Hausdorff distance under ambient
diffeomorphisms, but under stronger assumptions than the work of Federer,
namely assuming that ${\mathcal{S}}$ is a $C^{2}$ manifold and the distortion
is a $C^{2}$ diffeomorphism of the ambient space. Chazal and Soufflet based
their work on earlier results by Blaschke [4], which were not as strong as
Federer’s.

In this paper we extend the stability result of the medial axis. More
concretely, we generalize the result of Chazal and Soufflet [6] to arbitrary
closed sets and $C^{1,1}$ diffeomorphisms of the ambient space, and show that
the Hausdorff distance between the medial axes of the closed set and its image
is bounded in terms of Lipschitz constants associated to the diffeomorphism of
the ambient space. Our result follows from the work of Federer [8] and in fact
shortens the proof in [6] significantly.

We note that our bounds on the Hausdorff distance say nothing about the
topology of the medial axis, which is known to be highly unstable (see for
example [3]), although it preserves the homotopy type (see [11]).

### Motivation

The medial axis has many real world applications — among others, in: robot
motion planning, triangulation algorithms, graphics, shape recognition, and
shape learning, see for example [2, 7, 9, 10, 12, 13, 14, 15]. When proving
properties of algorithms based on the medial axis, authors generally assume
the real RAM model. However, it is not (a priori) clear if using (noisy) real
world data or the output of other computer programs as input for these
algorithms would yield answers that are close to the ground truth. To be able
to prove that the output is correct, we need (numerical) stability of medial
axis. In this paper we provide stability results using almost no assumptions
on the set itself.

### Outline

After a section with preliminaries we state and prove the main stability
result in Section 3. In Section 4 we reformulate this result in terms of norms
on Banach spaces. This also exhibits the fact that the stability of the medial
axis is Lipschitz in the following sense: We think of the set ${\mathcal{S}}$
as fixed and consider the map from the space of diffeomorphisms (endowed by a
norm which makes it a Banach space) to the space of closed subsets of
$\mathbb{R}^{d}$ (endowed with the Hausdorff distance), mapping each
diffeomorphism $F:\mathbb{R}^{d}\to\mathbb{R}^{d}$ to the closure of the
medial axis of $F({\mathcal{S}})$. The Lipschitz constant then only depends on
the diameter of the bounding ball of the set ${\mathcal{S}}$.

##  2 Preliminaries: Sets of positive reach and the closest point projection

In this section we recall some definitions and results concerning the medial
axis and sets of positive reach. Essentially, we need three ingredients from
the literature to prove our main theorem: the notions related to the closest
point projection, the properties of the generalized normal and tangent spaces,
and Federer’s result on the stability of the reach under ambient
diffeomorphisms.

We write $d(\cdot,\cdot)$ for the Euclidean distance between two points, and
the distance between a point and a set. That is, for any closed set
${\mathcal{S}}$ and point $p$,

$d(p,{\mathcal{S}})=\inf_{q\in{\mathcal{S}}}d(p,q).$

We denote the Hausdorff distance between two sets $A,B\subseteq\mathbb{R}^{d}$
by $d_{H}(A,B)$:

$d_{H}(A,B)=\max\left\\{\sup_{a\in A}d(a,B),\sup_{b\in B}d(b,A)\right\\}.$

We write $B(c,r)$, resp. $S(c,r)$, to denote balls, resp. spheres, with centre
$c$ and radius $r$. Lastly, $\left\lvert\cdot\right\rvert$ denotes the
Euclidean norm, and $\left\lVert\cdot\right\rVert$ an operator norm.

### The closest point projection and related notions

The projection of points in the ambient space $\mathbb{R}^{d}$ to the (set of)
closest point(s) of the set ${\mathcal{S}}\subseteq\mathbb{R}^{d}$ is denoted
by $\pi_{{\mathcal{S}}}$, and illustrated in Figure 2.

Figure 2: The closest point projection to the set ${\mathcal{S}}$ of four
points in $\mathbb{R}^{2}$. When a point lies on the medial axis
$\operatorname{{ax}}({\mathcal{S}})$, the closest point projection consists of
more points.

The _medial axis_ of ${\mathcal{S}}$ is the set of all points
$p\in\mathbb{R}^{d}$ where the set $\pi_{{\mathcal{S}}}(p)$ consists of more
than one point:

$\operatorname{{ax}}({\mathcal{S}})=\left\\{p\in\mathbb{R}^{d}\mid\\#\pi_{{\mathcal{S}}}(p)>1\right\\}.$

Here, $\\#\pi_{{\mathcal{S}}}(p)$ denotes the cardinality of the set
$\pi_{{\mathcal{S}}}(p)$.

For a point $p\in{\mathcal{S}}$, the _local feature size of $p$_ is the
distance from $p$ to the medial axis of the set ${\mathcal{S}}$:

$\operatorname{lfs}(p)=d(p,\operatorname{{ax}}({\mathcal{S}})).$

Finally, the _reach_ of the set ${\mathcal{S}}$ is the infimum of the local
feature size over all its points:

$\operatorname{rch}({\mathcal{S}})=\inf_{p\in{\mathcal{S}}}\operatorname{lfs}(p)=\inf_{p\in{\mathcal{S}}}d(p,\operatorname{{ax}}({\mathcal{S}})).$

Throughout this paper we assume that ${\mathcal{S}}\subseteq\mathbb{R}^{d}$ is
a closed set. We shall further assume that the set ${\mathcal{S}}$ as well as
its medial axis are bounded, and that the bounding sphere of ${\mathcal{S}}$
is contained in ${\mathcal{S}}$ itself. More specifically, we assume that
there exists a closed ball $B$ of positive radius such that
${\mathcal{S}}\subseteq B$, and $\partial B\subseteq{\mathcal{S}}$. We call
$\partial B$ the bounding sphere of ${\mathcal{S}}$ (see Figure 1 for
illustration).

The addition of the bounding sphere $\partial B$ to the set ${\mathcal{S}}$ is
necessary to obtain the desired bound on the Hausdorff distance between the
two medial axes of the set ${\mathcal{S}}$ and its image under the ambient
diffeomorphism. Indeed, consider the following example, illustrated in Figure
3.

Let the set ${\mathcal{S}}$ consist of two points in the plane,
${\mathcal{S}}=\\{p,q\\}\subseteq\mathbb{R}^{2}$. The medial axis of
${\mathcal{S}}$ is then the bisector line of $p$ and $q$. After a generic
perturbation $F$ of $p$ and $q$ — that is, not a translation and not a
perturbation in the direction $\pm(p-q)$ — the bisector line
$\operatorname{{ax}}(F({\mathcal{S}}))$ of the perturbed points intersects the
bisector $\operatorname{{ax}}({\mathcal{S}})$ of the original pair. The
Hausdorff distance between these two non-parallel lines is infinite, and thus
unboundable.

Figure 3: In black the set ${\mathcal{S}}$ and its medial axis, in light blue
the perturbed set and its medial axis. The Hausdorff distance between
$\operatorname{{ax}}({\mathcal{S}})$ and
$\operatorname{{ax}}(F({\mathcal{S}}))$ is infinite.

At the same time, the addition of the bounding sphere to the considered set
${\mathcal{S}}$ is not a restriction. Indeed, the medial axes of
${\mathcal{S}}$ and ${\mathcal{S}}\setminus\partial B$ are the same in the
interior of $B$ sufficiently far away from $\partial B$. More precisely:

  * •

Any point $x\in\operatorname{{ax}}({\mathcal{S}})$, such that
$\pi_{\mathcal{S}}(x)\cap\partial B=\emptyset$, lies on the medial axis
$\operatorname{{ax}}({\mathcal{S}}\setminus\partial B)$.

  * •

Conversely, if a point $x$ lies on the medial axis
$\operatorname{{ax}}({\mathcal{S}}\setminus\partial B)$, and any (and thus
every) point $q\in\pi_{{\mathcal{S}}\setminus\partial B}(x)$ satisfies
$d(x,q)<d(x,\partial B)$, then $x\in\operatorname{{ax}}({\mathcal{S}})$.

A recurring strategy in this article is to start at a point $p$ on the set
${\mathcal{S}}$, move away from this point in a ‘normal’ direction, and see if
by projecting using the closest point projection $\pi_{\mathcal{S}}$ we get
back to $p$. To this end, we define the _projection range_.

######  Definition 2.1 (Projection range).

Let $p\in{\mathcal{S}}$ be a point and $v\in\mathbb{R}^{d}$ a vector. The
_projection range_ $d(p,v,\pi_{\mathcal{S}})$ in direction $v$ is the maximal
distance one can travel from $p$ along $v$ such that the closest point
projection yields only the point $p$:

$d(p,v,\pi_{\mathcal{S}})=\sup\\{\lambda\in\mathbb{R}\mid\pi_{\mathcal{S}}(p+\lambda
v)=\\{p\\}\\}.$

Since $\pi_{\mathcal{S}}(p)=\\{p\\}$, the projection range is canonically non-
negative. Furthermore, the directions for which the range is positive are key
to our study, because of the following property:

######  Lemma 2.2 (Theorem 4.8 (6) of [8]).

Consider a point $p\in{\mathcal{S}}$ and a vector $v\in\mathbb{R}^{d}$. If

$\displaystyle 0<d(p,v,\pi_{\mathcal{S}})<\infty,$

then $p+d(p,v,\pi_{\mathcal{S}})\cdot
v\in\overline{\operatorname{{ax}}({\mathcal{S}})}$.

We call these special directions _back projection vectors_ :

######  Definition 2.3 (Unit back projection vectors).

For a point $p\in{\mathcal{S}}$, $\operatorname{UBP}(p,{\mathcal{S}})$ is the
set of unit vectors with a positive projection range:

$\operatorname{UBP}(p,{\mathcal{S}})=\left\\{u\in\mathbb{R}^{d}\mid\left\lvert
u\right\rvert=1\text{ and }0<d(p,u,\pi_{\mathcal{S}})<\infty\right\\}.$

We further define

$\displaystyle\operatorname{UBP}({\mathcal{S}})$
$\displaystyle=\left\\{(p,u)\in{\mathcal{S}}\times\mathbb{R}^{d}\,\middle|\,u\in\operatorname{UBP}(p,{\mathcal{S}})\right\\},$
$\displaystyle\operatorname{BP}({\mathcal{S}})$
$\displaystyle=\left\\{(p,\lambda
u)\in{\mathcal{S}}\times\mathbb{R}^{d}\,\middle|\,(p,u)\in\operatorname{UBP}({\mathcal{S}}),\lambda\geq
0\right\\}.$

Thanks to Lemma 2.2, the following map is well-defined:

$\displaystyle\pi_{\operatorname{{ax}},{\mathcal{S}}}:$
$\displaystyle\operatorname{UBP}({\mathcal{S}})\to\overline{\operatorname{{ax}}({\mathcal{S}})},$
$\displaystyle(p,u)\mapsto p+d(p,u,\pi_{\mathcal{S}})u.$ (1)

### The generalized tangent and normal space

Back projection vectors are intricately related to the generalized tangent and
normal spaces.

######  Definition 2.4 (Definitions 4.3 and 4.4 of [8]).

Let $p\in{\mathcal{S}}$. The _generalized tangent space_
$\operatorname{Tan}(p,{\mathcal{S}})$ is the set of vectors
$u\in\mathbb{R}^{d}$, such that either $u=0$ or, for every $\varepsilon>0$
there exists a point $q\in{\mathcal{S}}$ with

$\displaystyle 0<$ $\displaystyle|q-p|<\varepsilon$ and
$\displaystyle\left|\frac{q-p}{|q-p|}-\frac{u}{|u|}\right|<\varepsilon.$

The _generalized normal space_ $\operatorname{Nor}(p,{\mathcal{S}})$ consists
of vectors $v\in\mathbb{R}^{d}$ such that $\langle v,u\rangle\leq 0$ for all
$u\in\operatorname{Tan}(p,{\mathcal{S}})$. Vectors contained in the
generalized tangent, resp. normal, space are called _tangent_ , resp. _normal_
, _to ${\mathcal{S}}$ at $p$._

The generalized tangent and normal spaces are illustrated in Figure 4.

Figure 4: The (affine) generalized tangent and normal spaces of four points in
the set ${\mathcal{S}}\subset\mathbb{R}^{2}$, in light blue and violet,
respectively.

### Stability of the reach under ambient diffeomorphisms

Our last ingredient is the following result by Federer.

######  Theorem 2.5 (Stability of the reach under ambient diffeomorphisms,
Theorem 4.19 of [8]).

Pick two constants $0<t<\operatorname{rch}({\mathcal{S}})$ and $s>0$. If the
map

$F:\\{x\in\mathbb{R}^{d}\mid d(x,{\mathcal{S}})<s\\}\to\mathbb{R}^{n}$

is injective and continuously differentiable, and the maps $F$, $F^{-1}$, and
$DF$ are Lipschitz continuous with Lipschitz constants
$\operatorname{Lip}(F)$, $\operatorname{Lip}(F^{-1})$,
$\operatorname{Lip}(DF)$, respectively, then the reach
$\operatorname{rch}(F({\mathcal{S}}))$ of the image of the set ${\mathcal{S}}$
under the map $F$ is lower-bounded by

$\operatorname{rch}(F({\mathcal{S}}))\geq\min\left\\{\frac{s}{\operatorname{Lip}(F^{-1})},\frac{1}{\left(\frac{\operatorname{Lip}(F)}{t}+\operatorname{Lip}(DF)\right)\left(\operatorname{Lip}(F^{-1})\right)^{2}}\right\\}.$

##  3 Stability of the medial axis under ambient diffeomorphisms

In this section we prove the main result of this paper, Theorem 3.12. This
theorem extends earlier work by Chazal and Soufflet [6]. Its proof relies on
Federer’s result on the stability of the reach, Theorem 2.5. To give a more
geometrical interpretation we introduce the concept of a weakly tangent sphere
and ball, and a maximal empty weakly tangent ball.

######  Definition 3.1 (Weakly tangent sphere and ball).

Let $p\in{\mathcal{S}}$. A sphere is called _weakly tangent to ${\mathcal{S}}$
at $p$_ if it contains the point $p$ and its centre lies in the generalized
normal space $\operatorname{Nor}(p,{\mathcal{S}})$. In other words, spheres
weakly tangent to ${\mathcal{S}}$ at $p$ are spheres with centres $p+v$ and
radii $\left\lvert v\right\rvert$, for a vector
$v\in\operatorname{Nor}(p,{\mathcal{S}})$.

A ball is called _weakly tangent to ${\mathcal{S}}$ at $p$_ if its boundary
sphere is _weakly tangent to ${\mathcal{S}}$ at $p$_.

######  Remark 3.2.

Using the definition of $\operatorname{Nor}(p,{\mathcal{S}})$, a weakly
tangent ball can also be defined as follows: A ball $B(c,r)$ is weakly tangent
at $p$ if and only if its centre $c$ and radius $r$ satisfy

$\left(p+\operatorname{Tan}({\mathcal{S}},p)\right)\cap B(c,r)=\\{p\\}.$

We also remark:

######  Lemma 3.3.

Let $p\in{\mathcal{S}}$ and $v\in\mathbb{R}^{d}$, and suppose that for some
$\lambda>0$ we have $\pi_{\mathcal{S}}(p+\lambda v)\neq\\{p\\}$. Then, for all
$\lambda^{\prime}\geq\lambda$, we have
$\pi_{\mathcal{S}}(p+\lambda^{\prime}v)\neq\\{p\\}$ and for all
$\lambda^{\prime}>\lambda$, that
$p\notin\pi_{\mathcal{S}}(p+\lambda^{\prime}v)$.

######  Proof 3.4.

We first note that the statement is empty if $v=0$. Let $v\neq 0$, and
consider the nested family of balls
$B(p+\lambda^{\prime}v,\lambda^{\prime}|v|)$, parametrized by
$\lambda^{\prime}>0$.

Because $\pi_{\mathcal{S}}(p+\lambda v)\neq p$, the (closed) ball $B(p+\lambda
v,\lambda|v|)$ contains a point $q\in{\mathcal{S}}$ other than $p$. Since the
balls $B(p+\lambda^{\prime}v,\lambda^{\prime}|v|)$ are nested, the point $q$
lies inside every ball $B(p+\lambda^{\prime}v,\lambda^{\prime}|v|)$ with
$\lambda^{\prime}\geq\lambda$. Moreover, $q$ lies in the interior of
$B(p+\lambda^{\prime}v,\lambda^{\prime}|v|)$ for $\lambda^{\prime}>\lambda$.
Hence, for every $\lambda^{\prime}\geq\lambda$, we have that
$\pi_{\mathcal{S}}(p+\lambda^{\prime}v)\neq\\{p\\}$ and for
$\lambda^{\prime}>\lambda$, that
$p\notin\pi_{\mathcal{S}}(p+\lambda^{\prime}v)$.

Figure 5: Two families of balls weakly tangent to the set
${\mathcal{S}}\subset\mathbb{R}^{2}$ (in blue). Each family contains a unique
maximal empty ball (in purple). Notice that the centre of the maximal empty
ball weakly tangent at the point $p_{1}$ lies at the medial axis
$\operatorname{{ax}}({\mathcal{S}})$, while the centre of the maximal empty
ball weakly tangent at the point $p_{2}$ only lies at its closure,
$\overline{\operatorname{{ax}}({\mathcal{S}})}$.

Lemma 3.3 essentially tells us that a family of weakly tangent balls
$\left\\{B(p+\lambda v,\lambda|v|)\right\\}_{\lambda\geq 0}$ contains at most
one which is maximal with respect to inclusion among those whose interior is
disjoint from the set ${\mathcal{S}}$. Two such families are illustrated in
Figure 5.

We call such balls _maximal empty_. For the purpose of this article, we define
maximal empty balls in terms of unit back projection vectors (Definition 2.3).
To see that each maximal empty ball is indeed weakly tangent, we emphasise:

######  Lemma 3.5.

If $(p,v)\in\operatorname{BP}({\mathcal{S}})$, then
$(p,v)\in\operatorname{Nor}({\mathcal{S}})$. That is,

$\operatorname{BP}({\mathcal{S}})\subseteq\operatorname{Nor}({\mathcal{S}}).$

In particular, for any pair $(p,u)\in{\operatorname{UBP}({\mathcal{S}})}$ and
radius $\lambda\geq 0$, the ball $B(p+\lambda u,\lambda)$ is weakly tangent to
${\mathcal{S}}$.

######  Proof 3.6.

Because $(p,v)\in\operatorname{BP}({\mathcal{S}})$, there exists a ball
$B(p+\lambda v,\lambda|v|)$ whose intersection with ${\mathcal{S}}$ is empty.
To derive a contradiction, suppose that
$(p,v)\notin\operatorname{Nor}({\mathcal{S}})$. Then there exists a vector
$w\in\operatorname{Tan}(p,{\mathcal{S}})$ whose inner product with $v$ is
positive, and therefore the point $p+w$ lies into the interior of $B(p+\lambda
v,\lambda|v|)$, cf. Remark 3.2. This contradicts the fact that $B(p+\lambda
v,\lambda|v|)\cap{\mathcal{S}}=\emptyset$.

######  Remark 3.7.

For general closed sets, the converse of Lemma 3.5, that is,
$\operatorname{Nor}({\mathcal{S}})\subseteq\operatorname{BP}({\mathcal{S}})$,
is not true. One counter-example is the graph of the function
$x\mapsto|x|^{3/2}$ at the origin. However, the inclusion
$\operatorname{Nor}({\mathcal{S}})\subseteq\operatorname{BP}({\mathcal{S}})$
holds for sets of positive reach, thanks to Theorem 4.8 (12) of [8] (Lemma
B.1).

######  Definition 3.8 (Maximal empty weakly tangent ball).

Let $(p,u)\in{\operatorname{UBP}({\mathcal{S}})}$. A weakly tangent ball
$B(p+\lambda u,\lambda)$ is called _maximal empty_ to ${\mathcal{S}}$ if
$\lambda=d(p,u,\pi_{\mathcal{S}})$, or, equivalently, if
$\pi_{\operatorname{{ax}},{\mathcal{S}}}(p,u)=p+\lambda u$.

(Maximal empty) weakly tangent balls satisfy the following properties. Let
$(p,u)\in\operatorname{UBP}({\mathcal{S}})$.

  * •

For any radius $0<\lambda\leq d(p,u,\pi_{\mathcal{S}})$, the interior of the
ball $B(p+\lambda u,\lambda)$ is disjoint from the set ${\mathcal{S}}$. This
follows directly from Definition 3.8 and Lemma 3.3.

  * •

The centres of maximal empty weakly tangent balls lie on the closure of the
medial axis of ${\mathcal{S}}$. This is due to Lemma 2.2 and the definition of
the map $\pi_{\operatorname{{ax}},{\mathcal{S}}}$ (equation (1)).

The following lemma moreover tells us, that each point on the medial axis is a
centre of a maximal empty weakly tangent ball.

######  Lemma 3.9 (Surjectivity on $\operatorname{{ax}}({\mathcal{S}})$).

For any point $x\in\operatorname{{ax}}({\mathcal{S}})$ and
$p\in\pi_{{\mathcal{S}}}(x)$, there exists a vector
$u\in\operatorname{UBP}(p,{\mathcal{S}})$ such that
$\pi_{\operatorname{{ax}},{\mathcal{S}}}(p,u)=x$. In other words, $B(x,|x-p|)$
is a maximally empty weakly tangent ball. Moreover, we have that

$\operatorname{{ax}}({\mathcal{S}})\subseteq\pi_{\operatorname{{ax}},{\mathcal{S}}}\left(\operatorname{UBP}({\mathcal{S}})\right)\subseteq\overline{\operatorname{{ax}}({\mathcal{S}})}.$

######  Proof 3.10.

Let $Q=\pi_{{\mathcal{S}}}(x)$ be the subset of ${\mathcal{S}}$ that is
closest to $x$. Because $x\in\operatorname{{ax}}({\mathcal{S}})$, $Q$ contains
at least two points, one of them being $p$. We write $\lambda=\left\lvert
x-p\right\rvert$. Since ${\mathcal{S}}$ and
$\operatorname{{ax}}({\mathcal{S}})$ are disjoint, $\lambda>0$, and thus we
can define $u=\frac{x-p}{\lambda}$.

Since the interior of the ball $B(x,\lambda)$ does not intersect
${\mathcal{S}}$, it in particular does not intersect
$\operatorname{Tan}(p,{\mathcal{S}})$ and thus $B(x,\lambda)$ is weakly
tangent at $p$ by Remark 3.2. Let us now consider the nested family
$B(p+\lambda^{\prime}u,\lambda^{\prime})$ of weakly tangent balls at $p$. By
definition, $\partial B(x,\lambda)\cap{\mathcal{S}}=Q$ and therefore
$B(p+\lambda^{\prime}u,\lambda^{\prime})\cap{\mathcal{S}}=p$ for
$\lambda^{\prime}<\lambda$. At the same time, Lemma 3.3 yields that for
$\lambda^{\prime}>\lambda$, $p\notin\pi_{\mathcal{S}}(p+\lambda^{\prime}u)$.
Hence the projection range in direction $u$ equals
$d(p,u,\pi_{\mathcal{S}})=\lambda$ and we obtain
$\pi_{\operatorname{{ax}},{\mathcal{S}}}(p,u)=p+\lambda u=x$ directly from
Definition 3.8. The fact that
$\pi_{\operatorname{{ax}},{\mathcal{S}}}\left(\operatorname{UBP}({\mathcal{S}})\right)\subseteq\overline{\operatorname{{ax}}({\mathcal{S}})}$
is due to Lemma 2.2, as was already noted above the statement of this Lemma.

We are now almost ready to state our main theorem. Before phrasing the result,
we walk the reader through the assumptions and fix the notation on the way.
The assumptions are illustrated in Figure 6.

Figure 6: The setting of Theorem 3.12.

######  Assumption 3.11.

  * •

We assume that the set ${\mathcal{S}}$ has a bounding sphere of radius $r$,
which we denote by $S(r)$.

  * •

We consider a $C^{1}$ diffeomorphism $F:\mathbb{R}^{d}\to\mathbb{R}^{d}$ such
that the Lipschitz constants of $F$ and $F^{-1}$ are bounded by $L_{F}$, and
the Lipschitz constants of the differentials $DF$ and $DF^{-1}$ are bounded by
$L_{DF}$. We call such a diffeomorphism a $C^{1,1}$ diffeomorphism.

  * •

We further assume that the map $F$ leaves the bounding sphere $S(r)$
invariant, that is, $F(S(r))=S(r)$.

  * •

We pick a point $c\in\operatorname{{ax}}({\mathcal{S}})$, a point
$p\in\pi_{{\mathcal{S}}}(c)$, and write $\rho=|c-p|$. Observe that since
${\mathcal{S}}\cap\operatorname{{ax}}({\mathcal{S}})=\emptyset$, $\rho$ is
positive. By Lemma 3.9, the ball $B(c,\rho)$ is a maximal empty weakly tangent
ball to ${\mathcal{S}}$ at $p$. Moreover, we define $u=\frac{c-p}{\left\lvert
c-p\right\rvert}$ and note that $u\in{\operatorname{UBP}}(p,{\mathcal{S}})$.

  * •

We denote the tangent hyperplane to the boundary sphere of $B(c,\rho)$ at $p$
by $p+T$. The hyperplane $T$ is the orthocomplement of the vector $u$:
$T=u^{\perp}$.

  * •

We work with the unit vector at $F(p)$ that points inside the image of the
ball $B(c,\rho)$ and is orthogonal to the hyperplane $D_{p}F(T)$. We denote
this vector by $u^{\prime}$.

The following properties of $u$ and $u^{\prime}$ play an important rule in the
proof of the theorem:

######  Claim 1.

$\displaystyle
u^{\prime}=\frac{(D_{p}F^{t})^{-1}(u)}{|(D_{p}F^{t})^{-1}(u)|},$ (2)

where $D_{p}F^{t}$ is the transpose matrix (or the adjoint operator) of $DF$
at the point $p$, defined by

$\displaystyle\forall v_{1},v_{2},\qquad\langle
v_{1},D_{p}F(v_{2})\rangle=\langle D_{p}F^{t}(v_{1}),v_{2}\rangle.$

######  Claim 2.

Let $\|D_{p}F-\operatorname{Id}\|\leq\varepsilon<1$. Then the angle $\angle
u,u^{\prime}$ between the vectors $u$ and $u^{\prime}$ satisfies

$\cos\angle u,u^{\prime}\geq\sqrt{1-\varepsilon^{2}}.$

The proof of both of these properties is an exercise in linear algebra, which
we defer to the appendix.

######  Theorem 3.12.

Under the above assumptions, there exists a maximal empty weakly tangent ball
$B({c^{\prime}},\rho^{\prime})$ to the set $F({\mathcal{S}})$, whose boundary
sphere has an internal normal $u^{\prime}$. In particular, the ball
$B({c^{\prime}},\rho^{\prime})$ is tangent to the affine hyperplane
$F(p)+D_{p}F(T)$. Its radius $\rho^{\prime}$ is bounded by

$\rho^{\prime}\in\left[\frac{\rho}{(L_{F})^{3}+\rho
L_{DF}(L_{F})^{2}},\frac{(L_{F})^{3}\rho}{1-\rho L_{DF}(L_{F})^{2}}\right].$

If, moreover, the distortions of both $F$ and $DF$ are bounded, that is, for
all $x\in\mathbb{R}^{d}$,

$\displaystyle|F(x)-x|\leq\varepsilon_{1},\qquad\|DF_{x}-\operatorname{Id}\|\leq\varepsilon_{2}<1,$
(3)

and $r\cdot L_{DF}(L_{F})^{2}\leq 1/2$, then the map
$\pi_{\operatorname{{ax}},{\mathcal{S}}}$ satisfies

$\displaystyle|\pi_{\operatorname{{ax}},{\mathcal{S}}}(p,u)-\pi_{\operatorname{{ax}},F({\mathcal{S}})}(F(p),u^{\prime})|$
$\displaystyle\leq
2r\sqrt{1+(L_{F})^{6}\left(1+4rL_{DF}(L_{F})^{2}\right)^{2}-2(L_{F})^{3}\left(1+4rL_{DF}(L_{F})^{2}\right)\sqrt{1-(\varepsilon_{2})^{2}}}+\varepsilon_{1},$

Thus, the Hausdorff distance between the medial axes of ${\mathcal{S}}$ and
its image $F({\mathcal{S}})$ is bounded by

$\displaystyle
d_{H}(\operatorname{{ax}}({\mathcal{S}}),\operatorname{{ax}}(F({\mathcal{S}})))$
$\displaystyle\leq
2r\sqrt{1+(L_{F})^{6}\left(1+4rL_{DF}(L_{F})^{2}\right)^{2}-2(L_{F})^{3}\left(1+4rL_{DF}(L_{F})^{2}\right)\sqrt{1-(\varepsilon_{2})^{2}}}+\varepsilon_{1}.$
(4)

The bound $|F(x)-x|\leq\varepsilon_{1}$ is really necessary, because we want
our theorem to accommodate for rotations and translations, which rotate and
translate the medial axis without changing distances and hence have Lipschitz
constant 1. We further stress that if the diffeomorphism $F$ is close to the
identity, its Lipschitz constant satisfies $L_{F}\geq 1$, because by
assumption $F$ leaves the bounding sphere $S(r)$ invariant, and $L_{DF}$ is
close to zero.

######  Proof 3.13 (Proof of Theorem 3.12).

We first derive the bounds for the radius $\rho^{\prime}$. As the first step,
we apply Theorem 2.5 to the boundary sphere $S(c,\rho)$ of the maximal empty
weakly tangent ball $B(c,\rho)$. In particular, we can choose the constant $s$
in Theorem 2.5 arbitrarily large, and the constant $t$ arbitrarily close to
the reach $\operatorname{rch}(S(c,\rho))=\rho$, to obtain:

$\operatorname{rch}\left(F(S(c,\rho))\right)\geq\frac{1}{\left(\frac{L_{F}}{\rho}+L_{DF}\right)(L_{F})^{2}}=\frac{\rho}{(L_{F})^{3}+\rho
L_{DF}(L_{F})^{2}}=:\rho_{1}.$

This means that no open ball of radius $\rho_{1}$ tangent to the set
$F(S(c,\rho))$ actually intersects $F(S(c,\rho))$. In addition, since the set
$F(B(c,\rho))$ does not contain any points of $F({\mathcal{S}})$ in its
interior, no ball of radius $\rho_{1}$ that is tangent to $F(S(c,\rho))$ and
whose centre lies inside $F(S(c,\rho))$ contains any point of
$F({\mathcal{S}})$.

The unit vector $u^{\prime}\in DF_{p}(T)^{\perp}$ (defined in (2)) is defined
such that the point $F(p)+\rho_{1}u^{\prime}$ lies inside the distorted ball
$F(B(c,\rho))$. Due to the above observation, the ball
$B(F(p)+\rho_{1}u^{\prime},\rho_{1})$ is weakly tangent to $F({\mathcal{S}})$
at $F(p)$ and contains no points of $F({\mathcal{S}})$ in its interior.

Let us now consider the weakly tangent ball
$B(F(p)+\rho^{\prime\prime}u^{\prime},\rho^{\prime\prime})$, whose radius
$\rho^{\prime\prime}$ satisfies

$\rho^{\prime\prime}>\frac{(L_{F})^{3}\rho}{1-\rho
L_{DF}(L_{F})^{2}}=:\rho_{2}.$

To shorten up the notation, we set

$F(p)+\rho^{\prime\prime}u^{\prime}=:c^{\prime\prime}.$

To derive a contradiction, we assume that
$B(c^{\prime\prime},\rho^{\prime\prime})$ is maximal empty. This is equivalent
to assuming that
$\textrm{int}B(c^{\prime\prime},\rho^{\prime\prime})\cap{F}({\mathcal{S}})=\emptyset$,
and thus $B(c^{\prime\prime},\rho^{\prime\prime})$ is a maximal empty weakly
tangent ball to $F(p)$. Similarly to the beginning of the proof, we now apply
Theorem 2.5 to the map $F^{-1}$ and the boundary sphere
$S(c^{\prime\prime},\rho^{\prime\prime})=\partial
B(c^{\prime\prime},\rho^{\prime\prime})$. As a result, the reach of
$F^{-1}(S(c^{\prime\prime},\rho^{\prime\prime}))$ is at least

$\displaystyle\operatorname{rch}\left(F^{-1}(S(c^{\prime\prime},\rho^{\prime\prime}))\right)\geq$
$\displaystyle\frac{\frac{(L_{F})^{3}\rho}{1-\rho
L_{DF}(L_{F})^{2}}}{(L_{F})^{3}+\frac{(L_{F})^{3}\rho}{1-\rho
L_{DF}(L_{F})^{2}}L_{DF}(L_{F})^{2}}=\rho.$

We conclude that there exists a ball that is tangent to the set
$F^{-1}(S(c^{\prime\prime},\rho^{\prime\prime}))$ at $F^{-1}(F(p))=p$, whose
radius is larger than $\rho$, and that does not contain any points of
${\mathcal{S}}$ in its interior. This contradicts the fact that the ball
$B(c,\rho)$ is maximal empty, and completes the proof of the first part of the
statement.

We now prove the bounds on the distortion of the map
$\pi_{\operatorname{{ax}},{\mathcal{S}}}$. Let
$\rho^{\prime}\in[\rho_{1},\rho_{2}]$ be the radius of the maximal empty
weakly tangent ball at $F(p)$ in the direction $u^{\prime}$, and write
$c^{\prime}:=F(p)+\rho^{\prime}u^{\prime}$ for its centre. We stress that, as
a consequence of Lemma 2.2 (Theorem 4.8 (6) of [8]) ,
$c^{\prime}\in\overline{\operatorname{{ax}}(F({\mathcal{S}}))}$, but it is not
necessarily true that $c^{\prime}\in{\operatorname{{ax}}(F({\mathcal{S}}))}$.

The goal is to estimate the distance between the two centres
$c=\pi_{\operatorname{{ax}},{\mathcal{S}}}(p,u)$ and
$c^{\prime}=\pi_{\operatorname{{ax}},{\mathcal{S}}}(F(p),u^{\prime})$. Indeed,
since $c-p=\rho u$ and $c^{\prime}-F(p)=\rho^{\prime}u^{\prime}$,

$\displaystyle\left\lvert c-c^{\prime}\right\rvert$ $\displaystyle=\left\lvert
c-p+p-F(p)+F(p)-c^{\prime}\right\rvert=\left\lvert\rho
u+p-F(p)-\rho^{\prime}c^{\prime}\right\rvert$
$\displaystyle\leq\left\lvert\rho
u-\rho^{\prime}u^{\prime}\right\rvert+\left\lvert F(p)-p\right\rvert.$

Due to the bound (3), $\left\lvert F(p)-p\right\rvert\leq\varepsilon_{1}$.
Furthermore, thanks to Claim 2,

$\displaystyle\left\lvert\rho
u-\rho^{\prime}u^{\prime}\right\rvert^{2}=\rho^{2}+(\rho^{\prime})^{2}-2\rho\rho^{\prime}\cos\angle
u,u^{\prime}\leq\rho^{2}+(\rho^{\prime})^{2}-2\rho\rho^{\prime}\sqrt{1-(\varepsilon_{2})^{2}}.$

Recalling that $\rho^{\prime}\in[\rho_{1},\rho_{2}]$, we thus obtain

$\displaystyle\left\lvert\rho u-\rho^{\prime}u^{\prime}\right\rvert\leq$
$\displaystyle\max\left(\sqrt{\rho^{2}+(\rho_{1})^{2}-2\rho\,\rho_{1}\cos(\arcsin(\varepsilon_{2}))},\sqrt{\rho^{2}+(\rho_{2})^{2}-2\rho\,\rho_{2}\cos(\arcsin(\varepsilon_{2}))}\right)$
$\displaystyle=\max\left(\sqrt{\rho^{2}+(\rho_{1})^{2}-2\rho\,\rho_{1}\sqrt{1-(\varepsilon_{2})^{2}}},\sqrt{\rho^{2}+(\rho_{2})^{2}-2\rho\,\rho_{2}\sqrt{1-(\varepsilon_{2})^{2}}}\right).$

Hence,

$\displaystyle|c-c^{\prime}|\leq\max\left(\sqrt{\rho^{2}+(\rho_{1})^{2}-2\rho\,\rho_{1}\sqrt{1-(\varepsilon_{2})^{2}}},\sqrt{\rho^{2}+\left(\rho_{2}\right)^{2}-2\rho\,\rho_{2}\sqrt{1-(\varepsilon_{2})^{2}}}\right)+\varepsilon_{1}.$
(5)

As the last step, we simplify the expression (5) (at the cost of weakening the
bounds). For this, we assume that $\rho L_{DF}(L_{F})^{2}\leq 1/2$, so that

$\displaystyle\rho_{1}$ $\displaystyle=\frac{\rho}{(L_{F})^{3}+\rho
L_{DF}(L_{F})^{2}}\geq\frac{\rho}{(L_{F})^{3}}\left(1-\rho\frac{L_{DF}}{L_{F}}\right),$
(6) $\displaystyle\rho_{2}$ $\displaystyle=\frac{(L_{F})^{3}\rho}{1-\rho
L_{DF}(L_{F})^{2}}\leq\rho(L_{F})^{3}\left(1+2\rho L_{DF}(L_{F})^{2}\right),$
(7)

where we used that, for $x\in[0,1/2]$, $\frac{1}{1+x}\geq 1-x$ and
$\frac{1}{1-x}\leq 1+2x$. We note that both $\rho_{1}$ and $\rho_{2}$ tend to
$\rho$ as $L_{F}$ tends to $1$ and $L_{DF}$ tends to $0$. We now consider
$|\rho_{1}-\rho|$ and $|\rho_{2}-\rho|$, and claim that

$|\rho_{1}-\rho|,|\rho_{2}-\rho|\leq\rho(L_{F})^{3}\left(1+2\rho
L_{DF}(L_{F})^{2}\right)-\rho.$

For $|\rho_{2}-\rho|=\rho_{2}-\rho$, the claim holds thanks to (7). To
establish this for $|\rho_{1}-\rho|$ requires a small calculation:

$\displaystyle|\rho_{1}-\rho|=\rho-\rho_{1}$
$\displaystyle\leq\rho-\frac{\rho}{(L_{F})^{3}}\left(1-\rho\frac{L_{DF}}{L_{F}}\right)$
(due to (6)) $\displaystyle\leq\rho(L_{F})^{3}\left(1+2\rho
L_{DF}(L_{F})^{2}\right)-\rho$ (assuming the claim holds) $\displaystyle
2\rho$ $\displaystyle\leq\rho(L_{F})^{3}\left(1+2\rho
L_{DF}(L_{F})^{2}\right)+\frac{\rho}{(L_{F})^{3}}\left(1-\rho\frac{L_{DF}}{L_{F}}\right)$
(reformulating the previous inequality) $\displaystyle 2$
$\displaystyle\leq(L_{F})^{3}+\frac{1}{(L_{F})^{3}}+2\rho
L_{DF}(L_{F})^{5}-\rho\frac{L_{DF}}{(L_{F})^{4}},$

where the final inequality holds because $x^{3}+x^{-3}\geq 2$, and
$2x^{5}-x^{-4}\geq 0$, for $x\geq 1$. We now consider the function

$\displaystyle f(\delta)$
$\displaystyle=\rho^{2}+\rho^{2}(1+\delta)^{2}-2\rho^{2}\,(1+\delta)\sqrt{1-(\varepsilon_{2})^{2}}$
$\displaystyle=\rho^{2}\left(\delta^{2}+2\left(1-\sqrt{1-(\varepsilon_{2})^{2}}\right)\delta+2\left(1-\sqrt{1-(\varepsilon_{2})^{2}}\right)\right).$

The function $f$ is a second order polynomial in $\delta$ and because all
coefficients are positive, the maximum of $f$ on the interval
$[-\delta_{m},\delta_{m}]$ is a attained at $\delta_{m}$, that is,

$\displaystyle\sup_{\delta\in[-\delta_{m},\delta_{m}]}f(\delta)=f(\delta_{m}).$

By combining these results, we see that

$\displaystyle|c-c^{\prime}|$
$\displaystyle\leq\sqrt{f\left((L_{F})^{3}\left(1+2\rho
L_{DF}(L_{F})^{2}\right)-1\right)}+\varepsilon_{1}$
$\displaystyle=\sqrt{\rho^{2}+\left(\rho(L_{F})^{3}\left(1+2\rho
L_{DF}(L_{F})^{2}\right)\right)^{2}-2\rho\left(\rho(L_{F})^{3}\left(1+2\rho
L_{DF}(L_{F})^{2}\right)\right)\sqrt{1-(\varepsilon_{2})^{2}}}$
$\displaystyle\phantom{=}+\varepsilon_{1}$
$\displaystyle=\rho\sqrt{1+(L_{F})^{6}\left(1+2\rho
L_{DF}(L_{F})^{2}\right)^{2}-2(L_{F})^{3}\left(1+2\rho
L_{DF}(L_{F})^{2}\right)\sqrt{1-(\varepsilon_{2})^{2}}}+\varepsilon_{1}.$

Because both $f(\delta)$ and the bound (7) are monotone in $\rho$, and $\rho$
is bounded by the radius $r$ of the bounding sphere $S(r)$, we conclude that

$\displaystyle|c-c^{\prime}|$ $\displaystyle\leq
2r\sqrt{1+(L_{F})^{6}\left(1+4rL_{DF}(L_{F})^{2}\right)^{2}-2(L_{F})^{3}\left(1+4rL_{DF}(L_{F})^{2}\right)\sqrt{1-(\varepsilon_{2})^{2}}}+\varepsilon_{1}.$
(8)

For every point $c$ in $\operatorname{{ax}}({\mathcal{S}})$ we have found a
point $c^{\prime}$ in $\overline{\operatorname{{ax}}(F({\mathcal{S}}))}$ whose
distance is bounded by (8), and therefore the one-sided Hausdorff distance
between the two medial axes $\operatorname{{ax}}({\mathcal{S}})$ and
$\overline{\operatorname{{ax}}(F({\mathcal{S}}))}$ is bounded by the same
quantity. Because the symmetrical formulation of the statement, the same bound
holds for the Hausdorff distance.

##  4 Quantifying $C^{1,1}$ diffeomorphisms as deviations from identity

In this section we reformulate the main result in terms of norms on Banach
spaces. This reformulation offers a more theoretical insight, and we believe
the reformulated bounds may be easier to work with in certain applications.
Indeed, in the context of practical numerical computations, a bound on the
Lipschitz constant of an operator — or, at least, a modulus of continuity —
allows to control the condition number. This control is particularly useful
when we calculate with objects such as the medial axis [3], whose (numerical)
stability is often problematic in practice.

As we will see below, for this reformulation we somewhat strengthen our
assumptions.

We decompose a diffeomorphism $F$ into the identity map
$\mathds{1}_{\mathbb{R}^{d}}$ on $\mathbb{R}^{d}$, and a displacement field
$\varphi$:

$F=\mathds{1}_{\mathbb{R}^{d}}+\varphi.$

For the choice of the displacement field, we restrict ourselves to the vector
space $\mathcal{U}$ of all $C^{1,1}$ maps $\varphi$ from $\mathbb{R}^{d}$ to
$\mathbb{R}^{d}$ whose restriction to the exterior of a certain bounding ball
$B(r)$, $\mathbb{R}^{d}\setminus B(r)$, equals $0$.111This is a bit more
restrictive than assuming that the restriction to the bounding sphere $S(r)$
is $0$, but it simplifies matters in this section.

A natural norm associated to $\mathcal{U}$ is one that makes it a Banach
space. A typical choice, inherited from general Banach spaces of $C^{1,1}$
functions, would be for example, for $\varphi\in\mathcal{U}$,

$\|\varphi\|_{C^{1,1}}=\max\left(\|\varphi\|_{\infty},\,\|D\varphi\|_{\infty},\,\operatorname{Lip}(D\varphi)\,\right).$
(9)

Here we used the following notation:

  * •

$\|\varphi\|_{\infty}=\sup_{x\in\mathbb{R}^{d}}|\varphi(x)|$ denotes the sup
norm on $x\mapsto|\varphi(x)|$, where $|\cdot|$ is the Euclidean norm in
$\mathbb{R}^{d}$,

  * •

$\|D\varphi\|_{\infty}=\sup_{x\in\mathbb{R}^{d}}\|D\varphi(x)\|$ denotes the
sup norm on $x\mapsto\|D\varphi(x)\|$ , where $\|D\varphi(x)\|$ is the
operator norm induced by the Euclidean norm on $\mathbb{R}^{d}$.

  * •

We write $\operatorname{Lip}(D\varphi)$ for the Lipschitz semi-norm of
$D\varphi$. The Lipschitz semi-norms of $\varphi$ and $D\varphi$ are defined
as

$\operatorname{Lip}(\varphi)=\sup_{x,y\in\mathbb{R}^{d},\,x\neq
y}\frac{|\varphi(y)-\varphi(x)|}{|y-x|},$

and

$\operatorname{Lip}(D\varphi)=\sup_{x,y\in\mathbb{R}^{d},\,x\neq
y}\frac{\|D\varphi(y)-D\varphi(x)\|}{|y-x|}.$

The norm defined in (9) makes $\mathcal{U}$ into a Banach space, since every
Cauchy sequence in $\mathcal{U}$ has a limit in $\mathcal{U}$. In addition,
any function $\varphi\in\mathcal{U}$ satisfies:

$\displaystyle\operatorname{Lip}(\varphi)$
$\displaystyle=\|D\varphi\|_{\infty},$ (10)
$\displaystyle\|D\varphi\|_{\infty}$ $\displaystyle\leq
r\operatorname{Lip}(D\varphi),$ (11) $\displaystyle\|\varphi\|_{\infty}$
$\displaystyle\leq r\operatorname{Lip}(\varphi)\leq
r^{2}\operatorname{Lip}(D\varphi),$ (12)

since the restriction of $\varphi$ to $\mathbb{R}^{d}\setminus B(r)$ is $0$.
This in turn yields that

$\operatorname{Lip}(D\varphi)\leq\|\varphi\|_{C^{1,1}}\leq\max(1,r,r^{2})\operatorname{Lip}(D\varphi).$

Thus, in $\mathcal{U}$, the norm $\varphi\mapsto\operatorname{Lip}(D\varphi)$
is equivalent to the norm $\varphi\mapsto\|\varphi\|_{C^{1,1}}$.

We can now state slightly less general version of Theorem 3.12 in terms of the
Banach space $(\mathcal{U},\>\varphi\mapsto\operatorname{Lip}(D\varphi)\,)$.

######  Theorem 4.1.

Let ${\mathcal{S}}\subseteq\mathbb{R}^{d}$ be bounded by the ball $B(r)$ of
radius $r>0$, such that $S(r)=\partial B(r)\subseteq{\mathcal{S}}$. Let
further $F$ be a $C^{1,1}$ diffeomorphism from $\mathbb{R}^{d}$ to itself that
leaves the set $\mathbb{R}^{d}\setminus B(r)$ invariant, and define two
displacement fields $\varphi,\tilde{\varphi}\in\mathcal{U}$ such that
$F=\mathds{1}_{\mathbb{R}^{d}}+\varphi$ and
$(\mathds{1}_{\mathbb{R}^{d}}+\tilde{\varphi})\circ(\mathds{1}_{\mathbb{R}^{d}}+\varphi)=\mathds{1}_{\mathbb{R}^{d}}$.

Define
$\varepsilon=\max\left(\operatorname{Lip}(D\varphi),\operatorname{Lip}(D\tilde{\varphi})\right)$.

If $r\varepsilon\leq 1/4$ , the Hausdorff distance between the medial axes of
the set ${\mathcal{S}}$ and its image $F({\mathcal{S}})$ is bounded by

$d_{H}(\operatorname{{ax}}({\mathcal{S}}),\operatorname{{ax}}(F({\mathcal{S}})))\leq\left(1+\sqrt{50}\right)r^{2}\varepsilon+\mathcal{O}\left(r^{3}\varepsilon^{2}\right).$

In particular,

$d_{H}(\operatorname{{ax}}({\mathcal{S}}),\operatorname{{ax}}(F({\mathcal{S}}))){=\mathcal{O}\left(r^{2}\varepsilon\right)}.$

######  Remark 4.2.

Observe that the bound $\mathcal{O}\left(r^{2}\varepsilon\right)$ is
consistent with a scaling by factor $\lambda$:
${\mathcal{S}}\mapsto\lambda\,{\mathcal{S}}$,
$F(\cdot)\mapsto\lambda\,F(\cdot/\lambda)$. Under such a scaling, the radius
$r$ is multiplied by $\lambda$, while the Lipschitz constant
$\operatorname{Lip}(D\varphi)$ — and therefore $\varepsilon$ — is divided by
$\lambda$. Furthermore, the Hausdorff distance
$d_{H}(\operatorname{{ax}}({\mathcal{S}}),\operatorname{{ax}}(F({\mathcal{S}})))$
increases by a factor $\lambda$. By considering a diffeomorphism that
translates the set ${\mathcal{S}}\setminus S(r)$ while keeping the bounding
sphere $S(r)$ fixed, we see that this bound is asymptotically optimal.

######  Proof 4.3 (Proof of Theorem 4.1).

We denote $L_{\varphi}=\operatorname{Lip}(\phi)$ and Expressions (10), (11)
and (12) yield:

$\displaystyle L_{\varphi}$ $\displaystyle\leq r\varepsilon,\qquad
L_{DF}=\varepsilon,\qquad L_{F}\leq 1+L_{\varphi}\leq
1+r\varepsilon,\qquad\varepsilon_{1}\leq
r^{2}\varepsilon,\qquad\varepsilon_{2}\leq r\varepsilon.$ (13)

We deduce

$r\varepsilon\leq{1/4}\Longrightarrow r\varepsilon(1+r\varepsilon)^{2}\leq
1/2\Longrightarrow rL_{DF}(L_{F})^{2}\leq 1/2.$

Thus, the conditions of Theorem 3.12 are satisfied. Next, we reformulate the
inequality (4) of Theorem 3.12. The expression $E$ under the square root at
the right hand side of this inequality is:

$E=1+(L_{F})^{6}\left(1+4rL_{DF}(L_{F})^{2}\right)^{2}-2(L_{F})^{3}\left(1+4rL_{DF}(L_{F})^{2}\right)\sqrt{1-(\varepsilon_{2})^{2}}.$

By replacing $L_{F}$ by $1+L_{\varphi}$ in $E$, the constants, as well as the
degree-one terms in $L_{\varphi}$, ${r}L_{DF}$, and $\varepsilon_{2}$, cancel
out. More precisely,

$E=16r^{2}L_{DF}^{2}+r^{2}\varepsilon_{2}^{2}+24rL_{\varphi}L_{DF}+9L_{\varphi}^{2}+\mathcal{O}(|({r}L_{DF},L_{\varphi},\varepsilon_{2})|^{3}).$
(14)

Finally, by substituting inequalities (13) into (14), we obtain

$E\leq 50r^{2}\varepsilon^{2}+\mathcal{O}\left(r^{3}\varepsilon^{3}\right),$

and

$d_{H}(\operatorname{{ax}}({\mathcal{S}}),\operatorname{{ax}}(F({\mathcal{S}})))\leq\left(1+\sqrt{50}\right)r^{2}\varepsilon+{\mathcal{O}\left(r^{3}\varepsilon^{2}\right).}$

## References

  * [1] N. Amenta and M. Bern.  Surface reconstruction by Voronoi filtering.  Discrete & Computational Geometry, 22(4):481–504, Dec 1999.  doi:10.1007/PL00009475.
  * [2] Nina Amenta, Sunghee Choi, and Ravi Krishna Kolluri.  The power crust.  In Proceedings of the sixth ACM symposium on Solid modeling and applications, pages 249–266, 2001.
  * [3] Dominique Attali, Jean-Daniel Boissonnat, and Herbert Edelsbrunner.  Stability and computation of medial axes - a state-of-the-art report.  In Torsten Möller, Bernd Hamann, and Robert D. Russell, editors, Mathematical Foundations of Scientific Visualization, Computer Graphics, and Massive Data Exploration, pages 109–125, Berlin, Heidelberg, 2009. Springer Berlin Heidelberg.  doi:10.1007/b106657_6.
  * [4] W. Blaschke.  Kreis und Kugel.  Verlag von Veit und Comp., 1916.
  * [5] Harry Blum.  A Transformation for Extracting New Descriptors of Shape.  In Weiant Wathen-Dunn, editor, Models for the Perception of Speech and Visual Form, pages 362–380. MIT Press, Cambridge, 1967.
  * [6] F. Chazal and R. Soufflet.  Stability and finiteness properties of medial axis and skeleton.  Journal of Dynamical and Control Systems, 10(2):149–170, 2004.  URL: http://dx.doi.org/10.1023/B:JODS.0000024119.38784.ff, doi:10.1023/B:JODS.0000024119.38784.ff.
  * [7] Ilke Demir, Camilla Hahn, Kathryn Leonard, Geraldine Morin, Dana Rahbani, Athina Panotopoulou, Amelie Fondevilla, Elena Balashova, Bastien Durix, and Adam Kortylewski.  SkelNetOn 2019: Dataset and challenge on deep learning for geometric shape understanding.  In 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), pages 1143–1151, 2019.  doi:10.1109/CVPRW.2019.00149.
  * [8] H. Federer.  Curvature measures.  Transactions of the America mathematical Society, 93:418–491, 1959\.
  * [9] Seng-Beng Ho and Charles R Dyer.  Shape smoothing using medial axis properties.  IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI-8(4):512–520, 1986.  doi:10.1109/TPAMI.1986.4767815.
  * [10] Jean-Claude Latombe.  Robot motion planning, volume 124.  Springer Science & Business Media, 2012.
  * [11] André Lieutier.  Any open bounded subset of $\mathbb{R}^{n}$ has the same homotopy type as its medial axis.  Computer-Aided Design, 36(11):1029 – 1046, 2004.  Solid Modeling Theory and Applications.  URL: http://www.sciencedirect.com/science/article/pii/S0010448504000065, doi:https://doi.org/10.1016/j.cad.2004.01.011.
  * [12] Punam K Saha, Gunilla Borgefors, and Gabriella Sanniti di Baja.  A survey on skeletonization algorithms and their applications.  Pattern recognition letters, 76:3–12, 2016.
  * [13] Doron Shaked and Alfred M. Bruckstein.  Pruning medial axes.  Computer Vision and Image Understanding, 69(2):156 – 169, 1998\.  URL: http://www.sciencedirect.com/science/article/pii/S1077314297905981, doi:https://doi.org/10.1006/cviu.1997.0598.
  * [14] Andrea Tagliasacchi, Thomas Delame, Michela Spagnuolo, Nina Amenta, and Alexandru Telea.  3d skeletons: A state-of-the-art report.  In Computer Graphics Forum, volume 35, pages 573–597. Wiley Online Library, 2016.
  * [15] Yajie Yan, Kyle Sykes, Erin Chambers, David Letscher, and Tao Ju.  Erosion thickness on medial axes of 3d shapes.  ACM Transactions on Graphics, 35(4):38:1–38:12, July 2016.  URL: http://doi.acm.org/10.1145/2897824.2925938, doi:10.1145/2897824.2925938.

##  Appendix A Proofs of Claims 1 and 2

######  Proof A.1 (Proof of Claim 1).

Indeed, since any invertible matrix $A$ satisfies $(A^{t})^{-1}=(A^{-1})^{t}$,
one has:

$\displaystyle w\in D_{p}F(u^{\perp})$ $\displaystyle\iff\langle
D_{p}F^{-1}(w),u\rangle=0$ $\displaystyle\iff\langle
w,(D_{p}F^{-1})^{t}u\rangle=0$ $\displaystyle\iff\langle
w,u^{\prime}\rangle=0$ $\displaystyle\iff w\in u^{\prime\perp},$

and thus

$D_{p}F(u^{\perp})=u^{\prime\perp}.$ (15)

In other words, we have shown that $u^{\prime}$ is orthogonal to
$D_{p}F(u^{\perp})=D_{p}F(T)$.

Because

$\langle D_{p}F(u),(D_{p}F^{-1})^{t}(u)\rangle=\langle
D_{p}F^{-1}(D_{p}F(u)),u\rangle=\langle u,u\rangle>0,$

we deduce that $\langle D_{p}F(u),u^{\prime}\rangle>0$. This is in turn
equivalent to $u^{\prime}$ pointing towards the interior of $F(B(c,\rho))$.

######  Proof A.2 (Proof of Claim 2).

We first show that $\angle u,u^{\prime}<\pi/2$. Indeed, define the vector $w$
as

$w=(D_{p}F^{t})^{-1}(u),$

that is, the vector satisfying $u=D_{p}F^{t}(w)$. Then
$u^{\prime}=\frac{w}{\left\lvert w\right\rvert}$ (see equation (2)), and

$\displaystyle|w|\langle u,u^{\prime}\rangle$ $\displaystyle=\langle
u,w\rangle=\langle D_{p}F^{t}w,w\rangle$ $\displaystyle=\langle
w,D_{p}Fw\rangle=\left\lvert w\right\rvert^{2}+\langle
w,(D_{p}F-\operatorname{Id})w\rangle$ $\displaystyle\geq\left\lvert
w\right\rvert^{2}-\left\lvert w\right\rvert^{2}\|D_{p}F-\operatorname{Id}\|$
$\displaystyle>0.$ (because, by assumption, $\|DF_{p}-\operatorname{Id}\|<1$)

Thus, $\langle u,u^{\prime}\rangle>0$, and therefore $\angle
u,u^{\prime}<\pi/2$.

Figure 7: Since $\|v-D_{p}F(v)\|\leq\varepsilon\left\lvert v\right\rvert$, the
vector $D_{p}F(v)$ lies in the green ball $B(v,\varepsilon\left\lvert
v\right\rvert)$. Since $\varepsilon<1$, the angle between $v$ and $D_{p}F(v)$
is upper-bounded by $\arcsin\varepsilon<\pi/2$.

Furthermore, consider a vector $v\in u^{\perp}$. Since $\|v-D_{p}F(v)\| \leq
\|D_{p}F-\operatorname{Id}\||v|\leq\varepsilon|v|$, the angle between $v$ and
$D_{p}F(v)$ is upper-bounded by $\arcsin\varepsilon<\pi/2$, as illustrated in
Figure 7. This yields a bound on the angle between the tangent spaces
$u^{\perp}$ and $D_{p}F(u^{\perp})$:

$\sin\angle u^{\perp},D_{p}F(u^{\perp})=\sin\sup_{v\in u^{\perp},w\in
D_{p}F(u^{\perp})}\angle v,w\>\leq\varepsilon.$ (16)

Using (15) and (16) we deduce that:

$\sin\angle u,u^{\prime}=\sin\angle u^{\perp},u^{\prime\perp}\leq\varepsilon.$

Finally, since $\angle u,u^{\prime}<\pi/2$, $\cos\angle
u,u^{\prime}\geq\sqrt{1-\varepsilon^{2}}$. This concludes the proof.

##  Appendix B Federer’s tubular neighbourhood lemma

######  Lemma B.1 (Federer’s tubular neighbourhood lemma, Theorem 4.8 (12) of
[8]).

Let $p\in{\mathcal{S}}$ and $\operatorname{lfs}(p)>0$. The generalized normal
space to ${\mathcal{S}}$ at $p$ is characterized by the following property:
For any $\rho\in\mathbb{R}$ satisfying $0<\rho<\operatorname{lfs}(p)$,

$\displaystyle\operatorname{Nor}(p,{\mathcal{S}})=\\{\lambda
v\in\mathbb{R}^{d}\mid\lambda\geq
0,|v|=\rho,\pi_{\mathcal{S}}(p+v)=\\{p\\}\\}.$

In particular, $\operatorname{Nor}(p,{\mathcal{S}})$ is a convex cone. The
generalized tangent space $\operatorname{Tan}(p,{\mathcal{S}})$ is the convex
cone dual to $\operatorname{Nor}(p,{\mathcal{S}})$.


# Sequential Analysis of a finite number of Coherent states

Esteban Martínez Vargas<EMAIL_ADDRESS>Física Teòrica: Informació i
Fenòmens Quàntics, Departament de Física, Universitat Autònoma de Barcelona,
08193 Bellatera (Barcelona) Spain

###### Abstract

We investigate an advantage for information processing of ordering a set of
states over making a global quantum processing with a fixed number of copies
of coherent states. Suppose Alice has $N$ copies of one of two quantum states
$\sigma_{0}$ or $\sigma_{1}$ and she gives these states to Bob. Using the
optimal sequential test, the SPRT, we ask if processing the states in batches
of size $l$ is advantageous to optimally distinguish the two hypotheses. We
find that for the symmetric case
$\\{\left|{\gamma}\right\rangle,\left|{-\gamma}\right\rangle\\}$ there is no
advantage of taking any batch size $l$. We give an expression for the optimal
batch size $l_{\text{opt}}$ in the assymetric case. We give bounds
$l_{\text{min}}$ and $l_{\text{max}}$ for when $P_{S}\approx 1$.

##  I Introduction

The efficient detection of quantum phenomena is a matter of fundamental and
practical importance, useful to test a fundamental theory or create a precise
detector for a technological application, for example. This kind of problem
can be framed within the study of hypothesis testing [1]. This topic can be
generalized into quantum hypothesis testing [2, 3]. If the information is
stored in one of several quantum states the problem is usually called quantum
state discrimination [4] as the task concerns with differentiating these
quantum states. Optical systems are very relevant to quantum technologies [5,
6] and therefore the discrimination of optical quantum states is an important
topic of study [7].

Efficient detection implies the best use of the available resources for the
discovery of an event in a signal. A usual approach to analyze the efficiency
of a protocol is to fix a number $N$ of resources and find the apparatus that
minimizes the errors [8]. However, in practice, it is useful to consider
online, on-the-fly detectors such as change-point detection [9, 10, 11, 12].

In its simplest form, quantum state discrimination consists in being given a
state $\rho$ with the promise of being one of two possible states:
$\rho=\sigma_{0}$ or $\rho=\sigma_{1}$ (hypothesis 0 and 1 respectively) and
construct a quantum measurement that distinguishes them with the lowest
possible average error [13, 4, 14]. We call Type-I error for guessing
hypothesis 1 as true while it is false and Type-II for guessing hypothesis 0
as true while being false. Such measurement is described by a Positive
Operator Valued Measure (POVM). One can consider $N$ copies of states and form
tensor states $\rho^{\otimes N}$. The probability of success will be higher
with more copies, as more resources are available [15]. Then, the problem
changes to distinguish between two hypotheses with the lowest possible average
error using the given number of copies. However, the POVM might imply highly
entangled operators which can be hard to build. It is relevant and not trivial
to know how well different strategies behave with respect to the total number
of resources $N$.

Sequential analysis is a statistical framework that addresses the issue of
optimal resource handling [16]. In this framework, the desired error bounds on
the Type-I and Type-II errors are fixed beforehand and the number of average
samples needed to decide within these bounds. The protocol that minimizes this
average number of resources is called the Sequential Probability Ratio Test
(SPRT) [17].

Recently the framework of Sequential Analysis has been introduced to quantum
theory [18]. It considers the problem of having access to quantum measurements
of states. The bounds given in Ref. [18] give us the minimum number of
resources needed using quantum measurements. It was found afterward in [19]
that the general bounds are attainable with adaptive measurements. The present
work can be regarded as an extension of the sequential analysis program when
considering coherent states, which imply an infinite dimensional Hilbert
space. However, here we consider a fixed (or non-adaptive) protocol. The
problem we treat here uses the SPRT and asks for the probability that the
protocol stops with $N$ copies or less with the probabilities of Type-I and
Type-II errors being less than or equal to given probabilities $\alpha$ and
$\beta$ respectively.

Three different strategies are relevant to this work. First, we have the
_general_ case when all the $N$ copies are available at once. This case
includes possibly entangled operators for measurement. Then when the states
are available one by one, we have the _online_ scenario, which implies that
the protocol ignores if there is a horizon in the number of copies and
therefore is optimal at each step of the process (this protocol is well
described in [20]). Finally, there is the _sequential_ scenario, that uses the
SPRT and is closely related to the online one. A relevant difference is that
the SPRT is a test that minimizes the average number of resources needed.

In this article, we explore the freedom of using collective quantum strategies
on subsets of copies of coherent states. The collective strategy involves an
accumulation of information into one mode [21]. We have a setting as in Fig.
(1), Alice gives a state $\rho^{\otimes N}$ to Bob and we investigate if
slicing this set into $N/l$ batches of the state $\rho^{\otimes l}$ and
measuring them in an ordering given by a function $f$ is beneficial for Bob in
terms of distinguishing which state he was given: $\sigma_{0}^{\otimes N}$ or
$\sigma_{1}^{\otimes N}$. The function $f$ only represents the fact that we
are using a statistical method: the SPRT.

Figure 1: Bob slices the set $\rho^{\otimes N}$ into $N/l$ batches
$\rho^{\otimes l}$ and processes each batch in an ordering given by a function
$f$.

Intuition indicates that there should be a trade-off, as measurements with
more copies yield less error, however, if we make the batches too large we
will run out of copies for the SPRT, as Alice handles a finite number of
copies $N$. Therefore, given $N$ there must be an optimal batch size $l$ in
terms of the probability of successfully identifying the given state. We find
that this is not always the case, as there are relevant instances where all
values of $l$ are equivalent.

We first revisit the pure qubit case with unambiguous from Ref. [18] in Sec.
II. Afterward, we treat the problem with coherent states. In Sec. III.1 we
introduce basic notions of the SPRT. In section III.2 we explore the problem
of testing Gaussian distributions and calculate the probabilities for the SPRT
to stopping with $N_{0}$ copies or less depending on which hypothesis is true,
given bounds on the Type-I and Type-II errors. Then we introduce in Secs. IV.1
and IV.2 the problem of measuring coherent states and the quantum strategy of
accumulating the information of several copies into one mode. This leads to
the results of Sec. IV.3 where we explore the optimality of $l$ in several
cases. We end the article in Sec. V with the conclusions.

##  II Sequential Unambiguous POVM

In some cases, nonorthogonal states can be exactly distinguished if we allow
the possibility of outcomes that don’t give information. Such discrimination
protocols are called unambiguous [4]. Here we study an unambiguous protocol
for distinguishing two pure finite-dimensional states. Let us denote without
loss of generality, the two possible states as $\left|{\psi_{0}}\right\rangle$
and $\left|{\psi_{1}}\right\rangle$ as [22]

$\left|{\psi_{a}}\right\rangle=\cos\theta\left|{x}\right\rangle+(-1)^{a}\sin\theta\left|{y}\right\rangle,$
(1)

where we have written them in terms of an orthonormal basis
$\left|{x}\right\rangle$ and $\left|{y}\right\rangle$ of a two-dimensional
Hilbert space and an angle $\theta\in[0,\pi/4]$ between them. Let us denote
the overlap between them as $\langle{\psi_{0}}|{\psi_{1}}\rangle=\cos
2\theta=:c$. We use a three-outcome POVM because the protocol considered here
is unambiguous [4]. Following [18] we have the sequential probability of
success for unambiguously discriminating 2 hypotheses when $N$ copies are
available goes as $P_{S}^{UA}=1-c^{N}$. Remarkably, this is a result that
applies to a global strategy as well as for online strategies. This
equivalence implies that all batch sizes are equivalent. To see this last
statement imagine that Bob makes batches of size $l$ from the original set of
$N$ states. We would therefore have the states
$\left|{\Psi_{a}^{l}}\right\rangle=\left|{\psi_{a}}\right\rangle^{\otimes l}$.
The effective overlap between the redefined copies is
$C=\langle{\Psi_{0}^{l}}|{\Psi_{1}^{l}}\rangle=\langle{\psi_{0}}|{\psi_{1}}\rangle^{l}=c^{l}$.
We would therefore have $N/l$ batches. As we have batches of size $l$ then we
can see this fact as a redefinition of a copy. We have therefore the
probability of success for unambiguous discrimination of these batches as

$P_{S}^{UA}=1-C^{N/l}=1-c^{N}.$ (2)

The reason for the simple substitution on Eq. (2) is that the global
performance of the unambiguous protocol is achieved by an online strategy
[18]. The online strategy is to apply an unambiguous POVM for each available
copy. Being an unambiguous measurement then the probability of success with
$N$ copies coincide with the probability of stopping at step $N$ because this
measurement yields a zero error answer. Only if we get an inconclusive outcome
we would have to keep on measuring. However, we can wait to have all the $N$
copies and make a global unambiguous measurement and have a result with the
same success probability, therefore there is no gain in the ordering strategy
by Bob in the unambiguous protocol.

A drawback of using an unambiguous protocol is that despite that it yields a
no-error answer, the whole protocol has, in general, a lower probability of
success than a two-outcome POVM. The reason for this is that is a very
restrictive protocol. Also, for mixed states, unambiguous discrimination is
possible only in very restrictive cases.

##  III Probability that the SPRT stops with $N_{0}$ samples or less

###  III.1 Classical SPRT

We first review some basic notions of the SPRT theory by Wald [16]. Consider
that we have $N_{0}<\infty$ independent and identically distributed (i.i.d.)
samples $x_{i}$ of a random variable $X$ that follows the probability
distributions $p(x|0)$ or $p(x|1)$. The $\\{0,1\\}$ index denotes the
hypotheses $0$ or $1$ respectively. Observe that $N_{0}\neq N$, $N$ will
return afterward. We can define a useful variable

$z(x)=\log\frac{p(x|0)}{p(x|1)}$ (3)

where $\log$ denotes the natural logarithm. Thus, with a set of outputs
$\\{x_{1},\ldots,x_{n}\\}$, we have a set of values
$\\{z(x_{1}),\ldots,z(x_{n})\\}$ that we will denote as
$\\{z_{1},\ldots,z_{n}\\}$ for simplicity. At step $n$ we define

$Z_{n}:=\sum_{i=1}^{n}z_{i}.$ (4)

$Z_{n}$ is an example of what is known in the literature as Martingale [23],
which is a stochastic process whose mean value for step $n+1$ is the value of
step $n$. The SPRT consists in observing the value of $Z_{n}$ when a new
sample is available. If $Z_{n}\geq C_{0}$ we will accept hypothesis 0 as true.
If $Z_{n}\leq C_{1}$ we will accept hypothesis 1 as true. If
$C_{1}<Z_{n}<C_{0}$ continue sampling. It can be shown [16] that the bounds
$C_{0}$ and $C_{1}$ can be chosen such that the type I error probability is
$\leq\alpha$ and analogously, such that the type II error probability is
$\leq\beta$ for given $\alpha,\leavevmode\nobreak\ \beta\in[0,1]$. Defining

$A:=\frac{1-\beta}{\alpha}\quad\text{and}\quad B:=\frac{\beta}{1-\alpha},$ (5)

we have that in a very good approximation [16],

$C_{0}\approx\log A\quad\text{and}\quad C_{1}\approx\log B.$ (6)

The SPRT is the sequential test that requires fewer samples on average [17].

We are given $N_{0}$ samples and we restrict to the SPRT. The relevant
probabilities to calculate correspond

$P_{0}(Z_{N_{0}}\geq\log A)\quad\text{and}\quad P_{1}(Z_{N_{0}}\leq\log B)$
(7)

where $P_{i}$ correspond to the probability when hypothesis $i$ is true. Let
us suppose that we are given the hypothesis 0 and 1 with equal priors
therefore the total probability of success is

$P_{S}=\frac{1}{2}P_{0}(Z_{N_{0}}\geq\log
A)+\frac{1}{2}P_{1}(Z_{N_{0}}\leq\log B).$ (8)

###  III.2 Testing Gaussians

Suppose now that $X$ is normally distributed so that the probability
distribution when the hypothesis $i$ is true corresponds to

$p(x|i)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\theta_{i})^{2}}{2\sigma^{2}}}.$
(9)

It is straightforward to show that

$z(x)=\frac{1}{2\sigma^{2}}(2(\theta_{0}-\theta_{1})x+\theta_{1}^{2}-\theta_{0}^{2}).$
(10)

Recalling Eq. (4) we have that

$\sum_{i=1}^{N_{0}}x_{i}=\frac{Z_{N_{0}}2\sigma^{2}-(\theta_{1}^{2}-\theta_{0}^{2})}{2(\theta_{0}-\theta_{1})}.$
(11)

Suppose that each $x_{i}$ has mean $\theta$ and variance $\sigma^{2}$. Observe
that $\sum_{i=1}^{N_{0}}x_{i}$ is a sum of normally distributed random
variables, therefore it is a normally distributed variable with mean
$N_{0}\theta$ and variance $N_{0}\sigma^{2}$ [24].

The stopping condition for the SPRT $Z_{N_{0}}\geq\log A$ translates to

$\sum_{i=1}^{N_{0}}x_{i}\geq\frac{\log
A2\sigma^{2}-(\theta_{1}^{2}-\theta_{0}^{2})}{2(\theta_{0}-\theta_{1})}.$ (12)

Observe that the probability that a normally distributed variable $X$ to take
a value less than or equal $\lambda$ is given by the cumulative probability
$G(\lambda)$. As we want the probability that a variable takes a value less
than or equal to some lambda we need $1-G(\lambda)$. In terms of the Error
function [25] defined as

$\text{Erf}(y)=\frac{2}{\sqrt{\pi}}\int_{0}^{y}e^{-t^{2}}dt,$ (13)

we thus have the probability

$\displaystyle P_{0}($ $\displaystyle Z_{N_{0}}\geq\log A)=$
$\displaystyle\frac{1}{2}\left(1-\text{Erf}\left(\frac{2\sigma^{2}\log
A-(\theta_{1}^{2}-\theta_{0}^{2})-2N_{0}\theta_{0}(\theta_{0}-\theta_{1})}{2(\theta_{0}-\theta_{1})\sqrt{2N_{0}}\sigma}\right)\right).$
(14)

Analogously, we can calculate

$\displaystyle P_{1}($ $\displaystyle Z_{N_{0}}\leq\log B)=$
$\displaystyle\frac{1}{2}\left(1+\text{Erf}\left(\frac{2\sigma^{2}\log
B-(\theta_{1}^{2}-\theta_{0}^{2})-2N_{0}\theta_{1}(\theta_{0}-\theta_{1})}{2(\theta_{0}-\theta_{1})\sqrt{2N_{0}}\sigma}\right)\right).$
(15)

##  IV Coherent states

###  IV.1 Wigner function

A coherent state $\left|{\gamma}\right\rangle$ is described by a complex
number $\gamma$. In the phase space, we can write
$\gamma=q_{\gamma}+ip_{\gamma}$ with $q$ and $p$ denoting quadratures of the
electromagnetic field. The Wigner function of such a state is given by a
Gaussian [26]

$W_{\gamma}=\frac{2}{\pi}\text{exp}[-2(q-q_{\gamma})^{2}-2(p-p_{\gamma})^{2}].$
(16)

To detect a quadrature of the electromagnetic field one normally uses homodyne
detection, which allows us to detect intensity discrepancies in an
electromagnetic field. Explicitly we can detect [26]

$\Delta
I=\sqrt{2}|\gamma|\left\langle\frac{e^{-i\xi}\hat{a}+e^{i\xi}\hat{a}^{\dagger}}{\sqrt{2}}\right\rangle,$
(17)

for the angle $\xi$. Suppose we measure the quadrature $q$ with the coherent
state $\left|{\gamma}\right\rangle$, which corresponds to $\xi=0$ the
probability distribution is a Gaussian with mean $\theta=q_{\gamma}$ and
variance $\sigma^{2}=1/4$.

###  IV.2 Multiple copies

If multiple copies of coherent states are available we can accumulate the
information into one mode [21]. Consider a beam splitter of transmissivity $T$
and reflexivity $R$, if the coherent states $\left|{\gamma}\right\rangle$ and
$\left|{\delta}\right\rangle$ incide into the beam splitter it transforms to

$\left|{\gamma}\right\rangle\otimes\left|{\delta}\right\rangle\rightarrow\left|{\sqrt{T}\gamma+\sqrt{R}\delta}\right\rangle\otimes\left|{-\sqrt{R}\gamma+\sqrt{T}\delta}\right\rangle.$
(18)

Therefore, if $\gamma=\delta$ and we have a 50:50 beam splitter we get
$\left|{\sqrt{2}\gamma}\right\rangle\otimes\left|{0}\right\rangle$. In
general, $l$ copies can be concentrated into one mode. Suppose several beam
splitters are put one after another such that they perform the unitary
transformation $\left|{\gamma}\right\rangle^{\otimes
l}\rightarrow\left|{\sqrt{l}\gamma}\right\rangle\otimes\left|{0}\right\rangle^{l-1}$
[21]. To achieve this, the beam splitters must have transmissivities and
reflectivities given by

$T_{j}=\frac{j}{j+1},\quad R_{j}=\frac{1}{j+1}.$ (19)

###  IV.3 Optimal $l$

We return to the scenario of Fig. (1). Suppose that we are given $N_{0}=N/l$
batches of copies of coherent states. For each batch of $l$ states, we
implement the process of accumulation from section IV.2. Therefore, the
probability distributions we are comparing are given by

$p(x|i)=\sqrt{\frac{2}{\pi}}e^{-2(x-\sqrt{l}\theta_{i})^{2}},$ (20)

where $\theta_{i}$ is given by the real part of the coherent state
$\theta_{i}=\Re(\gamma_{i})$. Notice that
$\theta_{i}\rightarrow\sqrt{l}\theta_{i}$ with respect to section III.2.
Therefore, following Eqs. (14) and (15) we have

$\displaystyle P_{0}($ $\displaystyle Z_{N/l}\geq\log A)=$
$\displaystyle\frac{1}{2}\left(1-\text{Erf}\left(\frac{\frac{1}{2}\log
A-l(\theta_{1}^{2}-\theta_{0}^{2})-2N\theta_{0}(\theta_{0}-\theta_{1})}{\sqrt{2}(\theta_{0}-\theta_{1})\sqrt{N}}\right)\right)$
(21)

and

$\displaystyle P_{1}($ $\displaystyle Z_{N/l}\leq\log B)=$
$\displaystyle\frac{1}{2}\left(1+\text{Erf}\left(\frac{\frac{1}{2}\log
B-l(\theta_{1}^{2}-\theta_{0}^{2})-2N\theta_{1}(\theta_{0}-\theta_{1})}{\sqrt{2}(\theta_{0}-\theta_{1})\sqrt{N}}\right)\right).$
(22)

Observe that $P_{0}(Z_{N/l}\geq\log A)$ and $P_{1}(Z_{N/l}\leq\log B)$ depend
on $l$.

In Fig. (2) the SPRT is illustrated for several values of $l$. The Gaussian
distribution is a numerical approximation truncated in $\\{-10,10\\}$. In that
figure, we observe random realizations, some of which surpass the bound
corresponding to $\log A$, which correspond to the success instances. The mean
value of the sampling distribution corresponds to $\theta_{0}$ and thus we see
that the martingales tend to go upwards.

Figure 2: Examples of martingales $Z_{n}$ for different values of $l$
indicated above each plot. We take $1000$ samples in each figure. In each
figure $N=100$, $\theta_{0}=0.1$, $\theta_{1}=-0.1$, $\alpha=0.01$ and
$\beta=0.05$. The mean value of the sampling distribution is $\theta_{0}$. The
pale blue paths correspond to random realizations of $Z_{n}$. The dark blue is
the mean path over the $1000$ samples.

The cost function that needs to be optimized is the total probability from Eq.
(8). It remains to optimize it over $l$. To this end, we need to investigate
the sum of Error functions. Let us then define

$\displaystyle y_{A}$ $\displaystyle:=\frac{\frac{1}{2}\log
A-l(\theta_{1}^{2}-\theta_{0}^{2})-2N\theta_{0}(\theta_{0}-\theta_{1})}{\sqrt{2}(\theta_{0}-\theta_{1})\sqrt{N}}.$
(23) $\displaystyle y_{B}$ $\displaystyle:=\frac{\frac{1}{2}\log
B-l(\theta_{1}^{2}-\theta_{0}^{2})-2N\theta_{1}(\theta_{0}-\theta_{1})}{\sqrt{2}(\theta_{0}-\theta_{1})\sqrt{N}}.$
(24)

Therefore,

$P_{S}=\max_{l}\\{\frac{1}{4}\left(2-\text{Erf}(y_{A})+\text{Erf}(y_{B})\right)\\}.$
(25)

####  IV.3.1 Symmetric case

The frequently used Dolinar receiver [27] normally works with a symmetric pair
of coherent states
$\\{\left|{\gamma}\right\rangle,\left|{-\gamma}\right\rangle\\}$. If we are in
this symmetric case then we have that $\theta_{0}=-\theta_{1}$. This implies

$\displaystyle y_{A}$ $\displaystyle:=\frac{\frac{1}{2}\log
A-4N\theta_{0}^{2}}{2\sqrt{2}\theta_{0}\sqrt{N}}.$ (26) $\displaystyle y_{B}$
$\displaystyle:=\frac{\frac{1}{2}\log
B+4N\theta_{0}^{2}}{2\sqrt{2}\theta_{0}\sqrt{N}}.$ (27)

We see that there is no dependence on $l$, therefore any batch size is equally
good.

####  IV.3.2 Non-symmetric case

Suppose now that $\theta_{0}\neq-\theta_{1}$. In general, $P_{S}$ can have
three behaviors as shown in Fig. (3). We can change the value of $l$ such that
we move in the $x$ axis of the figures in question. The optimization over $l$
depends on the case we have at hand.

Figure 3: We plot the three possible cases for the behavior of the sum of
Error functions that are relevant for the probability of success $P_{S}$.

If we are in case I there is nothing to do, we have that $P_{S}\leq 1/2$. In
this case the best guess for the hypothesis at hand is random.

If we are in case II observe that there is a point where the sum of Error
functions attain a maximum. This maximum can be approximated with the Taylor
expansion of the exponential around 0. Using the Eq. (13) we obtain a Taylor
expansion for the Error function around 0

$\text{Erf}(y)=\frac{2}{\sqrt{\pi}}\sum_{n=0}^{\infty}\frac{(-1)^{n}y^{2n+1}}{n!(2n+1)}.$
(28)

At order zero, we see that

$\text{Erf}(y)\approx\frac{2}{\sqrt{\pi}}y.$ (29)

Using this, by symmetry, we can obtain an approximation to the optimal value
of $l$. Notice that the point where the zero-order approximation in case II in
Fig. (3) cross each other marks the optimal value of the sum of Error
functions. Therefore, the maximum of $P_{S}$ is found when

$-y_{A}\approx y_{B}.$ (30)

We thus approximate value for the optimal $l$, we define

$l_{\text{opt}}:=N+\frac{\log A+\log B}{4(\theta_{1}^{2}-\theta_{0}^{2})}.$
(31)

This value only makes sense when

$0\leq l_{\text{opt}}\leq N.$ (32)

In Fig. (4) we have a graph of the total probability of success $P_{S}$
dependent on $l$ and see that it attains its maximum at $l_{\text{opt}}$ given
by Eq. (31).

Figure 4: For this example we have $\alpha=0.00005$, $\beta=0.2$,
$\theta_{0}=0.2$, $\theta_{1}=-0.1$ and $N=100$.

If we are in case III there are limits for $x$ were $P_{S}\approx 1$ in Fig.
(3) as $x_{\text{min}}$ and $x_{\text{max}}$. This implies bounds for $l$ that
are defined as follows

$\displaystyle-\frac{2}{\sqrt{\pi}}y_{A}$ $\displaystyle=1$
$\displaystyle\frac{2}{\sqrt{\pi}}y_{B}$ $\displaystyle=1.$ (33)

These equations give the limits

$\displaystyle l_{\text{min}}$
$\displaystyle:=\frac{1}{2(\theta_{1}+\theta_{0})}\left(\frac{\log
B}{\theta_{1}-\theta_{0}}+4N\theta_{1}+\sqrt{2N\pi}\right)$ $\displaystyle
l_{\text{max}}$
$\displaystyle:=\frac{1}{2(\theta_{1}+\theta_{0})}\left(\frac{\log
A}{\theta_{1}-\theta_{0}}+4N\theta_{0}-\sqrt{2N\pi}\right).$ (34)

These bounds are only defined for

$\displaystyle 0$ $\displaystyle\leq l_{\text{min}}\leq N,$ $\displaystyle 0$
$\displaystyle\leq l_{\text{max}}\leq N.$ (35)

##  V Conclusions

We extend the study of sequential analysis protocols for coherent states.
Specifically, we study the probability that a specific statistical test, the
SPRT accepts one of two possible hypotheses with $N$, a given number of
resources. In so doing we investigate the duality of collective measurements
with many copies and the necessity of having to process the measurements
optimally with the SPRT.

We find that in the symmetric case,
$\\{\left|{\gamma}\right\rangle,\left|{-\gamma}\right\rangle\\}$ there is no
advantage of taking batches of any size. In contraposition with the adaptive
protocol used in the Dolinar receiver [27] the protocol we consider here is
non-adaptive. The independence with respect to $l$ seems to come from the fact
that we are considering optimal sequential processing. For non-symmetric
cases, two cases are relevant to us. In the first one, there is a unique $l$
that achieves the maximum labeled $l_{\text{opt}}$, which is approximated
using the Taylor expansion of the Error function. The second relevant case
implies a range of values of $l$ for which, using the expansion of the Error
function we define a lower bound $l_{\text{min}}$ and an upper bound
$l_{\text{max}}$ for the range of values of $l$ that attain the optimal
$P_{S}$.

Operationally speaking, the SPRT shows an advantage when considering small
type-I and type-II error probabilities. However, notice that the protocol we
are considering is more general than only making a collective, entangled
measurement. The reason for this is that the batch could be of size $N$
_always_ i.e. $l=N$. We show that in general, this is not the case and that
there is an advantage when taking into account the statistical process.

The treatment here was with the most simple quantum strategy that involves
only pure states and fixed measurements. Perhaps an adaptive strategy in the
measurement apparatus gives more insight into when sequential information
processing is necessary [19].

The results from this work could be generalized to the mixed-state case for
finite-dimension states. It would be necessary a calculation of the
probabilities of Eqs. (14) and (15).

##  VI Acknowledgements

I want to acknowledge useful discussions and suggestions from Ramon Muñoz-
Tapia, Gael Sentís and John Calsamiglia.

## References

  * [1] Thomas M. Cover and Joy A. Thomas.  Elements of Information Theory (Wiley Series in Telecommunications and Signal Processing).  Wiley-Interscience, USA, 2006.
  * [2] Tomohiro Ogawa and Hiroshi Nagaoka.  Strong Converse and Stein’s Lemma in Quantum Hypothesis Testing, pages 28–42.  WORLD SCIENTIFIC, 2 2005.
  * [3] Sarah Brandsen, Kevin D. Stubbs, and Henry D. Pfister.  Reinforcement learning with neural networks for quantum multiple hypothesis testing.  In 2020 IEEE International Symposium on Information Theory (ISIT). IEEE, 6 2020.
  * [4] Stephen M. Barnett and Sarah Croke.  Quantum state discrimination.  Adv. Opt. Photon., 1(2):238–278, Apr 2009.
  * [5] V. Giovannetti, S. Guha, S. Lloyd, L. Maccone, J. H. Shapiro, and H. P. Yuen.  Classical capacity of the lossy bosonic channel: The exact solution.  Physical Review Letters, 92, 1 2004.
  * [6] Christian Weedbrook, Stefano Pirandola, Raúl García-Patrón, Nicolas J. Cerf, Timothy C. Ralph, Jeffrey H. Shapiro, and Seth Lloyd.  Gaussian quantum information.  Reviews of Modern Physics, 84:621–669, 5 2012.
  * [7] J. Calsamiglia, R. Muñoz Tapia, Ll. Masanes, A. Acin, and E. Bagan.  Quantum chernoff bound as a measure of distinguishability between density matrices: Application to qubit and gaussian states.  Physical Review A, 77, 3 2008.
  * [8] C.W. Helstrom.  Quantum Detection and Estimation Theory.  Mathematics in Science and Engineering : a series of monographs and textbooks. Academic Press, 1976.
  * [9] A. Tartakovsky, I. Nikiforov, and M. Basseville.  Sequential Analysis: Hypothesis Testing and Changepoint Detection.  Chapman & Hall/CRC Monographs on Statistics & Applied Probability. Taylor & Francis, 2014.
  * [10] Gael Sentís, Emilio Bagan, John Calsamiglia, Giulio Chiribella, and Ramon Muñoz Tapia.  Quantum change point.  Physical Review Letters, 117, 10 2016.
  * [11] Gael Sentís, John Calsamiglia, and Ramon Muñoz Tapia.  Exact identification of a quantum change point.  Physical Review Letters, 119, 10 2017.
  * [12] Esteban Martínez Vargas and Ramon Muñoz Tapia.  Certified answers for ordered quantum discrimination problems.  Physical Review A, 100, 10 2019.
  * [13] Joonwoo Bae and Leong-Chuan Kwek.  Quantum state discrimination and its applications.  Jul 2017.  J. Phys. A: Math. Theor. 48 083001 (2015).
  * [14] Anthony Chefles.  Quantum state discrimination.  Contemporary Physics, 41(6):401–424, 2000.
  * [15] K. M. R. Audenaert, J. Calsamiglia, R. Muñoz Tapia, E. Bagan, Ll. Masanes, A. Acin, and F. Verstraete.  Discriminating states: The quantum chernoff bound.  Physical Review Letters, 98, 4 2007.
  * [16] A. Wald.  Sequential Analysis.  Dover books on advanced mathematics. Dover Publications, 1973.
  * [17] A. Wald and J. Wolfowitz.  Optimum character of the sequential probability ratio test.  The Annals of Mathematical Statistics, 19(3):326–339, 1948.
  * [18] Esteban Martínez Vargas, Christoph Hirche, Gael Sentís, Michalis Skotiniotis, Marta Carrizo, Ramon Muñoz Tapia, and John Calsamiglia.  Quantum sequential hypothesis testing.  Physical Review Letters, 126, 5 2021.
  * [19] Yonglong Li, Vincent Y. F. Tan, and Marco Tomamichel.  Optimal adaptive strategies for sequential quantum hypothesis testing.  Communications in Mathematical Physics, 392:993–1027, 6 2022.
  * [20] Gael Sentís, Esteban Martínez-Vargas, and Ramon Muñoz-Tapia.  Online identification of symmetric pure states.  Quantum, 6:658, February 2022.
  * [21] Michal Sedlák, Mário Ziman, Vladimír Bužek, and Mark Hillery.  Unambiguous comparison of ensembles of quantum states.  Physical Review A, 77, 4 2008.
  * [22] A. Acín, E. Bagan, M. Baig, Ll. Masanes, and R. Muñoz Tapia.  Multiple-copy two-state discrimination with individual measurements.  Physical Review A, 71, 3 2005.
  * [23] Michael Mitzenmacher and Eli Upfal.  Probability and Computing: Randomized Algorithms and Probabilistic Analysis.  Cambridge University Press, 2005.
  * [24] D.S. Lemons, P. Langevin, and A. Gythiel.  An Introduction to Stochastic Processes in Physics.  Johns Hopkins Paperback. Johns Hopkins University Press, 2002.
  * [25] G.B. Arfken and H.J. Weber.  Mathematical Methods for Physicists.  Elsevier, 2005.
  * [26] Pieter Kok and Brendon W. Lovett.  Introduction to Optical Quantum Information Processing.  Cambridge University Press, 2010.
  * [27] S. J. Jr. Dolinar.  An optimum receiver for the binary coherent state quantum channel.  Quarterly Progress Report, 111, 1973.