# The Square Root Normal Field Distance and Unbalanced Optimal Transport

Martin Bauer, Emmanuel Hartman, Eric Klassen

###### Abstract.

This paper explores a novel connection between two areas: shape analysis of
surfaces and unbalanced optimal transport. Specifically, we characterize the
square root normal field (SRNF) shape distance as the pullback of the
Wasserstein-Fisher-Rao (WFR) unbalanced optimal transport distance. In
addition we propose a new algorithm for computing the WFR distance and present
numerical results that highlight the effectiveness of this algorithm. As a
consequence of our results we obtain a precise method for computing the SRNF
shape distance directly on piecewise linear surfaces and gain new insights
about the degeneracy of this distance.

###### 2020 Mathematics Subject Classification:

49Q10, 49Q22

Martin Bauer and Emmanuel Hartman were partially supported by NSF-grants
1912037 and 1953244.

###### Contents

  1. 1 Introduction
  2. 2 Unbalanced optimal transport and the Wasserstein-Fisher-Rao distance
    1. 2.1 The Wasserstein-Fisher-Rao distance
    2. 2.2 The subspace of all finitely supported measures
    3. 2.3 Flatness of the space of measures with a fixed number of support points
    4. 2.4 A convex optimization problem
    5. 2.5 Implementation and Experiments
  3. 3 The SRNF shape metric as an unbalanced transport problem
    1. 3.1 Shape spaces of surfaces
    2. 3.2 The SRNF framework
    3. 3.3 SRNF shape metric as an unbalanced optimal transport problem.
    4. 3.4 SRNF Computation Experiments
  4. 4 Conclusion

##  1\. Introduction

This paper contributes to two different areas: elastic shape analysis (ESA)
[55] and unbalanced optimal mass transport [13, 35]. The main results of our
article are twofold: first we develop a new algorithm for the numerical
computation of the Wasserstein-Fisher-Rao distance [35, 12, 25] (a form of
unbalanced optimal mass transport), and secondly we establish a connection
between these two areas, which in turn allows for the exact computation of the
Square Root Normal Field distance, which is a widely used similarity measure
in ESA of surfaces; see e.g. [21, 32, 30, 22] and the references therein.
Before we describe the contributions of the present article in more detail, we
will briefly discuss the background of these two fields.

### Background:

In mathematical shape analysis, one is interested in quantifying and
describing the differences between geometric objects, such as point clouds,
geometric curves, or unparametrized surfaces [62, 2, 55, 16]. The main sources
of difficulty in this area are the high (infinite) dimensionality and the non-
linearity of such spaces; e.g., the shape space of surfaces is an (infinite
dimensional) function space modulo several finite and infinite dimensional
group actions. Consequently, even simple operations such as addition or
averaging are not well-defined on such spaces. Riemannian geometry has been
proven to provide a successful framework to tackle this challenging task: in a
Riemannian viewpoint, one considers the space of all shapes of interest
(geometric objects) as an infinite dimensional manifold and equips it with an
(infinite dimensional) Riemannian metric, thereby encoding the invariances of
the objects in the geometry and building a convenient setup for subsequent
statistical analysis. In the context of geometric curves or surfaces this
approach is often referred to as Elastic Shape Analysis (ESA) [55, 61],
despite the analogy to elastic stretching and bending energies being only
loose [46, 47].

In this article we will focus on elastic shape analysis of surfaces, i.e., we
consider Riemannian metrics on the quotient space
$\operatorname{Imm}(M,\mathbb{R}^{3})/\operatorname{Diff}(M)$ of immersions
modulo reparametrizations, where $M$ is a compact two dimensional manifold
(the parameter space), $\operatorname{Imm}(M,\mathbb{R}^{3})$ denotes the
space of immersions of $M$ into $\mathbb{R}^{3}$ and $\operatorname{Diff}(M)$
is the diffeomorphism group of the parameter space. One can define a
Riemannian metric on the quotient spaces, by considering a reparametrization
invariant metric on the space of immersions, such that the projection is a
Riemannian submersion. Over the past years there has been a significant amount
of work dedicated to studying the mathematical properties of such metrics and
in particular sufficient conditions to guarantee non-degeneracy of the
geodesic distance [39] and local wellposedness [4] of the geodesic equations
have been derived. However, the analogs of the global existence and
completeness results, that have been derived in the case of planar curves [9,
8, 33] are still missing. From an application point of view, the most
important task is a fast and robust implementation of the geodesic boundary
value problem, which in the setup of geometric statistics serves as the basis
for any subsequent statistical analysis [41]. In general, the absence of
explicit formulas for geodesics makes this a highly non-trivial task.
Motivated by similar results in the case of geometric curves, several
simplifying transformations have been proposed that locally flatten the
Riemannian metric [26, 27, 22]. The most successful among these is the so-
called Square Root Normal Field (SRNF) transformation [22], which assigns an
invariant (pseudo) distance to the space of immersions by considering the
$L^{2}$-distance between appropriately weighted normal vector fields. The
corresponding shape distance can then be calculated by minimizing over the
action of the reparametrization group, see Section 3.2 for an exact definition
of this framework. Based on the resulting computational ease and convincing
results [32, 3], the SRNF framework has been proven successful in a variety of
applications, see e.g. [28, 23, 37, 31]. In a recent paper [24] certain
degeneracy results for the resulting distance have been characterized, but a
more detailed theoretical study of its properties is still missing; the second
part of this article will contribute towards this aim.

The optimal mass transport (OMT) problem was first formulated as a non-convex
optimization problem on the space of transport maps by Monge in 1781 [40].
Since then a large amount of work has been dedicated to gaining a better
theoretical understanding of this challenging model; we refer to the
monographs [59, 58] for a detailed introduction to the field. Over the past
years OMT has proven successful in a variety of applications, ranging from
Computer Vision to Image Analysis and in particular Statistics and Data
Science, see e.g. [45, 19, 54, 53, 42] and the reference therein. Fueled by
these applications efficient numerical discretizations of OMT have been
developed [42, 5, 38, 10] including in particular the celebrated Sinkhorn
algorithm [15, 52], which efficiently solves an entropic regularized version
of the OMT problem.

The original formulation of OMT is rooted in the assumption that both
densities have the same total mass. Motivated by applications, where this can
be a limiting factor, various formulations of OMT that lift this restriction
have been proposed [43, 36, 25, 35, 12]; such transportation problems are also
called unbalanced transport problems. In particular, a new family of metrics
that interpolates between the Wasserstein and the Fisher-Rao metric has been
introduced in [35, 12, 25]. The theoretical properties of this model, called
Wasserstein-Fisher-Rao distance (WFR) or Hellinger-Kantorovich distance, have
been studied in detail in [34, 14, 35] and, as with traditional optimal
transport, efficient Sinkhorn-type, entropy regularized methods have been
introduced in [13].

### Contributions:

We start our presentation by reviewing the Kantorovich formulation of the WFR
distance, where we will focus on the induced distance on the subspace of all
finitely supported measures on $S^{2}$. In this setting the computation of the
WFR distance reduces to a convex optimization problem on the space of discrete
semi-couplings, i.e., on the space of pairs of constrained matrices, see
Section 2.2 and Lemma 2.4. In Section 2.3, we prove that the set of measures
on $S^{2}$ with a fixed number of support points is a locally flat metric
space with respect to the WFR metric. This result is an extension of a result
in [12], in which the same result is proved for measures on convex Euclidean
domains instead of $S^{2}$. In Section 2.4 we use the formulations of Section
2.2 to develop an efficient coordinate descent algorithm, whose convergence is
ensured by the convexity of the problem. We can find an explicit solution to
the optimization problem when restricted to the space of semi-couplings with a
fixed first (second, resp.) matrix, cf. Lemma 2.4. This gives rise to a
simple, numerically efficient algorithm to compute the optimal semi-coupling.
In Section 2.5 we then present an open source pytorch111Our code is available
at https://github.com/emmanuel-hartman/WassersteinFisherRaoDistance.
implementation of this algorithm and compare it in several experiments to the
entropic regularized Sinkhorn solver of [13].

In the second part of our article, we focus on the SRNF distance. We start by
presenting a unified framework for the SRNF distance that allows us to
incorporate both smooth and piecewise linear surfaces (simplicial complexes).
We then show that this extended framework coincides with the original SRNF
distance when restricted to smooth surfaces, cf. Theorem 3.1. The main
contribution of the second part, which establishes a connection between the
SRNF distance on the space of unparametrized surfaces and the Wasserstein-
Fisher-Rao distance on the space of Borel measures on $S^{2}$, is presented
next; more precisely, we construct a map from the space of piecewise linear
surfaces to the space of finitely supported Borel measures on $S^{2}$, such
that the SRNF shape distance is the pullback of the Wasserstein-Fisher-Rao
distance via this map, see Theorem 3.4. The central building block of this
result is related to the theory of area measures [1], which have a long
history in convex geometry and in particular in Brunn-Minkowski theory, cf.
[49] and the references therein. In the context of shape analysis of curves,
area measures have been recently studied in [11].

Theorem 3.4 highlights the degeneracy of the SRNF distance: in the recent
paper [24] it has been shown that there exist families of non-equivalent
closed surfaces that are indistinguishable by the SRNF distance. Our result
shows that for any closed surface, there exists a unique convex surface such
that the SRNF distance cannot distinguish them. It turns out that this convex
surface is exactly the solution of the well-known Minkowski problem [48, 50],
which allows us to use algorithms of convex geometry to present examples of
such pairs of surfaces that are indistinguishable by the SRNF, cf. Fig. 4.

As a second outcome of Theorem 3.4, we obtain a new algorithm for the precise
SRNF distance computation by reducing it to the solution of the Wasserstein-
Fisher-Rao distance. Thus the algorithm developed in the first part of the
paper directly applies to this situation. In the final section of this
article, we present numerical experiments and, in particular, a comparison to
previous implementations of the SRNF distance.

### Acknowledgements.

We thank Nicolas Charon, Ian Jermyn, Cy Maor, Zhe Su, François-Xavier Vialard,
and the statistical shape analysis group at FSU for helpful discussions during
the preparation of this manuscript.

##  2\. Unbalanced optimal transport and the Wasserstein-Fisher-Rao distance

In this section, we will first recall the Kantorovich formulation of the
recently proposed Wasserstein-Fisher-Rao distance. We will then discuss the
restriction of this distance to the space of finitely supported measures on
$S^{2}$. In our main result of this section, we will construct an efficient
splitting algorithm for the computation of this distance. We will prove the
convergence of our algorithm using the result that the computation of this
distance can be reduced to optimizing a concave function over a finite-
dimensional convex set.

###  2.1. The Wasserstein-Fisher-Rao distance

In recent years, there has been a concerted effort by the optimal transport
community to extend the definition of well-studied classical optimal transport
to unbalanced problems, i.e. to transport problems that allow for expansion
and compression of mass. We will consider a specific example of such a
generalization called the Wasserstein-Fisher-Rao distance that was introduced
independently by [12] and [35]. In the following, we will discuss the
corresponding Kantorovich formulation, as introduced in [14], for the special
case of measures on $S^{2}$.

Therefore we denote by $\mathcal{M}_{+}(S^{2})$ the space of finite Borel
measures on $S^{2}$. To formulate the Kantorovich problem for unbalanced
transport we introduce the notion of a semi-coupling, which is a direct
generalization of the notion of a coupling, which is used in standard OMT:

######  Definition 2.1 (Semi-couplings [14]).

Given $\mu,\nu\in\mathcal{M}_{+}(S^{2})$ the set of all semi-couplings is
given by

$\Gamma(\mu,\nu)=\left\\{(\gamma_{0},\gamma_{1})\in\mathcal{M}_{+}(S^{2}\times
S^{2})^{2}|(\operatorname{Proj}_{0})_{\\#}\gamma_{0}=\mu,(\operatorname{Proj}_{1})_{\\#}\gamma_{1}=\nu\right\\}.$
(2)

The Wasserstein-Fisher-Rao distance from $\mu$ to $\nu$ can be defined as the
infimum of a functional on the space of semi-couplings of $\mu$ and $\nu$.

######  Definition 2.2 (Wasserstein-Fisher-Rao Distance [14, 35]).

The Wasserstein-Fisher-Rao Distance on $\mathcal{M}_{+}(S^{2})$ is given by

$\displaystyle\operatorname{WFR}:\mathcal{M}_{+}(S^{2})\times\mathcal{M}_{+}(S^{2})\to\mathbb{R}^{+}\text{
defined via }$ (3)
$\displaystyle(\mu,\nu)\mapsto\inf\limits_{(\gamma_{0},\gamma_{1})\in\Gamma(\mu,\nu)}\left(\int_{S^{2}\times
S^{2}}\left|\sqrt{\frac{\gamma_{0}}{\gamma}(u,v)}u-\sqrt{\frac{\gamma_{1}}{\gamma}(u,v)}v\right|^{2}d\gamma(u,v)\right)^{1/2}$
(4)

where $\gamma\in\mathcal{M}_{+}(S^{2}\times S^{2})$ such that
$\gamma\ll\gamma_{0},\gamma_{1}$.

The following theorem summarizes the main result about this distance function:

######  Theorem 2.3 (Properties of the Wasserstein-Fisher-Rao Distance [12,
14, 35]).

Given any $\mu,\nu\in\mathcal{M}_{+}(S^{2})$ there exists an optimal
$(\gamma_{0},\gamma_{1})\in\Gamma(\mu,\nu)$ such that

$\operatorname{WFR}(\mu,\nu)^{2}=\int_{S^{2}\times
S^{2}}\left|\sqrt{\frac{\gamma_{0}}{\gamma}(u,v)}u-\sqrt{\frac{\gamma_{1}}{\gamma}(u,v)}v\right|^{2}d\gamma(u,v).$
(5)

Further, $(\mathcal{M}_{+}(S^{2}),\operatorname{WFR})$ is a geodesic length
space.

###  2.2. The subspace of all finitely supported measures

We will now consider the restriction of this metric to the subset of all
finitely supported, finite measures on $S^{2}$

$\mathcal{M}^{<\infty}_{+}(S^{2})=\left\\{\sum_{i=1}^{m}a_{i}\delta_{u_{i}}:m\in\mathbb{N},a_{i}\in\mathbb{R}^{+},\text{
and }u_{i}\in S^{2}\right\\}\subset\mathcal{M}_{+}(S^{2}),$ (6)

where $\delta_{u}$ is the Dirac measure at $u\in S^{2}$. Note that for
$\mu,\nu\in\mathcal{M}^{<\infty}_{+}(S^{2})$ an arbitrary semi-coupling
$(\gamma_{0},\gamma_{1})\in\Gamma(\mu,\nu)$ is not required to be finitely
supported. We will introduce a subset of discrete semi-couplings and show that
optimizing over all valid semi couplings is equivalent to optimizing over our
restricted subset.

######  Definition 2.4 (Discrete semi-couplings).

Let $\mu,\nu\in\mathcal{M}^{<\infty}_{+}(S^{2})$ with
$\mu=\sum_{i=1}^{m}a_{i}\delta_{u_{i}}$ and
$\nu=\sum_{j=1}^{n}b_{j}\delta_{v_{j}}$. A discrete semi-coupling of $\mu$ and
$\nu$ is a pair of $(m+1)\times(n+1)$ matrices $(A,B)$ satisfying the
properties:

  1. (a)

for all $i,j$, $A_{ij}\geq 0$ and $B_{ij}\geq 0$;

  2. (b)

for all $i=1,\dots,m$, $a_{i}=A_{i0}+\dots+A_{in}$;

  3. (c)

for each $j=0,\dots,n$, $A_{0j}=0$;

  4. (d)

for all $j=1,\dots,n$, $b_{j}=B_{0j}+\dots+B_{mj}$;

  5. (e)

for each $i=0,\dots,m$, $B_{i0}=0$.

We denote the set of all discrete semi-couplings of $\mu$ and $\nu$ by
$\mathcal{A}(\mu,\nu)$.

Note that the discrete semi-couplings from $\mu$ to $\nu$ represent a proper
subset of $\Gamma(\mu,\nu)$. Thus, to show that we can compute
$\operatorname{WFR}(\mu,\nu)$ by simply optimizing over $\mathcal{A}(\mu,\nu)$
we need the following lemma.

######  Lemma 2.5.

Let $\mu,\nu\in\mathcal{M}^{<\infty}_{+}(S^{2})$ with
$\mu=\sum_{i=1}^{m}a_{i}\delta_{u_{i}}$ and
$\nu=\sum_{j=1}^{n}b_{j}\delta_{v_{j}}$. Let
$(\gamma_{0},\gamma_{1})\in\Gamma(\mu,\nu)$; then there exists
$(A,B)\in\mathcal{A}(\mu,\nu)$ such that

$\int_{S^{2}\times
S^{2}}\left|\sqrt{\frac{\gamma_{0}}{\gamma}(u,v)}u-\sqrt{\frac{\gamma_{1}}{\gamma}(u,v)}v\right|^{2}d\gamma(u,v)=\sum_{i=0}^{m}\sum_{j=0}^{n}|\sqrt{A_{ij}}u_{i}-\sqrt{B_{ij}}v_{j}|^{2}$
(7)

where $u_{0}=v_{0}=(1,0,0)\in S^{2}$.

###### Proof.

Construct $(A,B)$ as follows:

  * •

For $i\in\\{1,...,m\\}$ and $j\in\\{1,...,n\\}$,
$A_{ij}=\gamma_{0}/\gamma(u_{i},v_{j})$ and
$B_{ij}=\gamma_{1}/\gamma(u_{i},v_{j})$,

  * •

for $i\in\\{1,...,m\\}$,
$A_{i0}=a_{i}-\sum_{j=1}^{n}\gamma_{0}/\gamma(u_{i},v_{j})$,

  * •

for $j\in\\{1,...,n\\}$,
$B_{0j}=b_{j}-\sum_{i=1}^{m}\gamma_{1}/\gamma(u_{i},v_{j}),$

  * •

for $i\in\\{0,...,m\\}$, $B_{i0}=0,$ and

  * •

for $j\in\\{0,...,n\\}$, $A_{0j}=0.$

Note that by construction $(A,B)\in\mathcal{A}(\mu,\nu)$. Observe,

$\displaystyle\int\limits_{S^{2}\times
S^{2}}\left|\sqrt{\frac{\gamma_{0}}{\gamma}}u-\sqrt{\frac{\gamma_{1}}{\gamma}}v\right|^{2}d\gamma(u,v)$
$\displaystyle=\int\limits_{S^{2}\times
S^{2}}\left|\sqrt{\frac{\gamma_{0}}{\gamma}}\right|^{2}d\gamma(u,v)+\int\limits_{S^{2}\times
S^{2}}\left|\sqrt{\frac{\gamma_{1}}{\gamma}}\right|^{2}d\gamma(u,v)$ (8)
$\displaystyle\qquad-2\int\limits_{S^{2}\times
S^{2}}\sqrt{\frac{\gamma_{0}}{\gamma}}\sqrt{\frac{\gamma_{1}}{\gamma}}(u\cdot
v)\gamma(u,v)$ (9)
$\displaystyle=\sum_{i=1}^{m}a_{i}+\sum_{j=1}^{n}b_{j}-2\sum_{i=1}^{m}\sum_{j=1}^{n}\sqrt{A_{ij}}\sqrt{B_{ij}}(u_{i}\cdot
v_{j})$ (10)
$\displaystyle=\sum_{i=0}^{m}\sum_{j=0}^{n}|\sqrt{A_{ij}}u_{i}-\sqrt{B_{ij}}v_{j}|^{2}$
(11)

∎

######  Corollary 2.6.

Let $\mu,\nu\in\mathcal{M}^{<\infty}_{+}(S^{2})$ with
$\mu=\sum_{i=1}^{m}a_{i}\delta_{u_{i}}$ and
$\nu=\sum_{j=1}^{n}b_{j}\delta_{v_{j}}$. Then

$\operatorname{WFR}(\mu,\nu)^{2}=\inf\limits_{(A,B)\in\mathcal{A}(\mu,\nu)}\sum_{i=0}^{m}\sum_{j=0}^{n}|\sqrt{A_{ij}}u_{i}-\sqrt{B_{ij}}v_{j}|^{2}$
(12)

######  Remark 2.7.

By writing out the norms inside the summation and excluding terms that sum to
zero, one obtains the following alternative formula for WFR:

$\operatorname{WFR}(\mu,\nu)=\left(\sum_{i=1}^{m}a_{i}+\sum_{j=1}^{n}b_{j}-2\sup_{(A,B)\in\mathcal{A}}\sum_{i=1}^{m}\sum_{j=1}^{n}\sqrt{A_{ij}B_{ij}}(u_{i}\cdot
v_{j})\right)^{1/2}.$ (13)

######  Remark 2.8.

Given a discrete semi-coupling $(A,B)$, the zeroth column of $A$ and zeroth
row of $B$ are included to handle the case where all of the mass at the
corresponding support is destroyed/created rather than transported. This does
not mean, however, that these rows/columns being zero correspond to no
creation/destruction of mass. We will demonstrate this fact in the following
example.

######  Example 2.9.

We consider the example with $\mu=2\delta_{u_{1}}+\delta_{u_{2}}$ and
$\nu=\delta_{v_{1}}+3\delta_{v_{2}}$ with

$u_{1}=v_{1}=\left(\begin{array}[]{c}1\\\ 0\\\ \end{array}\right),\qquad
u_{2}=v_{2}=\left(\begin{array}[]{c}0\\\ 1\\\ \end{array}\right).$

The corresponding weight matrix between the supports is then given by

$\Omega=\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}.$ (14)

Consequently the optimal semi-coupling is given by:

$A=\begin{pmatrix}0&0&0\\\ 0&2&0\\\ 0&0&1\end{pmatrix}\text{ and
}B=\begin{pmatrix}0&0&0\\\ 0&1&0\\\ 0&0&3\end{pmatrix}.$ (15)

Notice that even though the zeroth column of $A$ and zeroth row of $B$ are all
zeros, a unit of mass is destroyed and two units of mass are created. For
contrast we consider a second example with the same masses but different
supports: $\mu=2\delta_{u_{1}}+\delta_{u_{2}}$ and
$\nu=\delta_{v_{1}}+3\delta_{v_{2}}$ where

$u_{1}=v_{1}=\left(\begin{array}[]{c}1\\\ 0\\\ \end{array}\right),\qquad
u_{2}=\left(\begin{array}[]{c}0\\\ 1\\\ \end{array}\right),\qquad
v_{2}=-u_{2}=-\left(\begin{array}[]{c}0\\\ 1\\\ \end{array}\right).$

This time the corresponding weight matrix between the supports is given by

$\Omega=\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix}$ (16)

and an optimal semi-coupling is given by:

$A=\begin{pmatrix}0&0&0\\\ 0&2&0\\\ 1&0&0\end{pmatrix}\text{ and
}B=\begin{pmatrix}0&0&3\\\ 0&1&0\\\ 0&0&0\end{pmatrix}.$ (17)

In this case, two units of mass are destroyed (one from each support of $\mu$)
and three units are created (all three are created at the second support of
$\nu$).

###  2.3. Flatness of the space of measures with a fixed number of support
points

Note that while the main result of this section proves a fundamental fact
about the geometry of the space of finitely supported measures, this result is
not a prerequisite for the rest of the paper. Let ${\mathcal{M}}^{n}_{+}$
denote the set of measures on $S^{2}$ that are supported at precisely $n$
points. It is clear that ${\mathcal{M}}^{n}_{+}$ has a natural structure as a
$3n$-manifold. In the following theorem, we will show that the restriction of
the WFR metric to ${\mathcal{M}}^{n}_{+}$ is locally isometric to the distance
function corresponding to a flat Riemannian metric on this space.

######  Theorem 2.10.

There is a flat Riemannian metric on ${\mathcal{M}}^{n}_{+}$ whose Riemannian
distance function agrees with the Wasserstein-Fisher-Rao metric on a small
neighborhood of every point.

Furthermore, let

${\mathcal{P}}_{n}=\left\\{(x_{1},\dots,x_{n})\in\left(\mathbb{R}^{3}-\\{0\\}\right)^{n}:\text{
for all }i\neq
j,\frac{x_{i}}{|x_{i}|}\neq\frac{x_{j}}{|x_{j}|}\right\\}\subset\left(\mathbb{R}^{3}/\\{0\\}\right)^{n}$

and define

$Q_{n}:{\mathcal{P}}_{n}\to{\mathcal{M}}^{n}_{+}\text{ by
}Q_{n}(x_{1},\dots,x_{n})=\sum_{i=1}^{n}\|x_{i}\|^{2}\delta_{\left(x_{i}/\|x_{i}\|\right)}.$
(18)

Then, for every $x=(x_{1},\dots,x_{n})\in{\mathcal{P}}_{n}$ there exists an
$\epsilon>0$ such that for all $y,z\in{\mathcal{P}}_{n}$ with $|y-x|<\epsilon$
and $|z-x|<\epsilon$ the optimal discrete semi-coupling from $Q_{n}(y)$ to
$Q_{n}(z)$ is a pair of diagonal matrices.

######  Remark 2.11.

Theorem 2.10 is closely related to Theorems 4.1 and 4.2 of [12]. The main
difference is that in [12], the domain of the measure is assumed to be a
convex region in $\mathbb{R}^{d}$, whereas in Theorem 2.10, the domain is
$S^{2}$. Another difference is that in [12], a lower bound on the size of the
neighborhood on which the metric is flat is given. The proof given here is
more elementary (a straightforward application of differential topology) and
thus we hope it is of interest in itself.

Before we are able to prove Theorem 2.10, we need the following technical
lemma.

######  Lemma 2.12.

Let $M$ be an $m$-dimensional manifold, and $C$ a compact subset of $M$. Let
$f:\mathbb{R}^{n}\times M\to\mathbb{R}$ be a $C^{\infty}$ function; let
$x_{0}\in C\subset M$. Define $g:M\to\mathbb{R}$ by $g(x)=f(0,x)$. Assume that
$g$ restricted to $C$ attains an absolute minimum at $x_{0}$, and that $x_{0}$
is the only point of $C$ where $g$ achieves this minimum. Also, assume that
$g$ has a critical point at $x_{0}$ and that the Hessian of $g$ at $x_{0}$ is
positive definite. (Clearly this assumption is independent of which chart
containing $x_{0}$ is used to compute the Hessian.) Now, suppose that there
exists $\zeta>0$ such that for all $p\in\mathbb{R}^{n}$ with $|p|<\zeta$, the
function $M\to\mathbb{R}$ defined by $x\mapsto f(p,x)$ has a critical point at
$x_{0}$. Then there exists $\epsilon>0$ such that for all $p\in\mathbb{R}^{n}$
with $|p|<\epsilon$, the function $C\to\mathbb{R}$ defined by $x\mapsto
f(p,x)$ attains an absolute minimum at $x_{0}$.

###### Proof.

Suppose the lemma is false. Then there exist sequences
$\\{p_{k}\\}\in\mathbb{R}^{n}$ with $p_{k}\to 0$ and $\\{x_{k}\\}\in C$ such
that for all $k\in\mathbb{N}$, $f(p_{k},x_{k})<f(p_{k},x_{0})$. By compactness
of $C$, we may choose a subsequence of $\\{x_{k}\\}$ that converges in $C$.
Continue to denote this subsequence by $\\{x_{k}\\}$ (and denote the
corresponding subsequence of $\\{p_{k}\\}$ by $\\{p_{k}\\}$). Let $c=\lim
x_{k}$. There are two cases to consider:

Case 1. Suppose $c\neq x_{0}$. In that case, by the continuity of $f$,
$f(0,c)=\lim f(p_{k},x_{k})\leq\lim f(p_{k},x_{0})=f(0,x_{0})$, contradicting
the hypothesis $g$ achieves its minimum on $C$ only at the point $x_{0}$.

Case 2. Suppose $c=x_{0}$. By taking subsequences, we may assume that all of
the $x_{k}$ lie in a single chart of $M$; hence, for the rest of this proof we
will replace $M$ by $U\subset\mathbb{R}^{m}$ and $x_{0}$ by
$0\in\mathbb{R}^{m}$. Recall that one assumption in our Lemma is that the
Hessian of $g$ at $0$ is positive definite; denote its smallest eigenvalue by
$\lambda$. (Using the chart, we are now thinking of $g$ as a function
$U\to\mathbb{R}$.) Clearly the condition that the lowest eigenvalue of the
Hessian of the map $y\mapsto f(p,y)$ at $x$ is greater than $\lambda/2$ is an
open condition on $(p,x)$. Hence we can choose $\delta>0$ so that for all
$|p|<\delta$ and $|x|<\delta$, this condition is satisfied. Choose $k_{0}$
such that $|p_{k_{0}}|<\delta$ and $|x_{k_{0}}|<\delta$. Define the map
$\tilde{g}:N_{\delta}(0)\to\mathbb{R}$ by $\tilde{g}(x)=f(p_{k_{0}},x)$.
Define $L:[0,1]\to\mathbb{R}$ by $L(t)=f(p_{k_{0}},tx_{k_{0}})$. Using the
chain rule, we see that
$L^{\prime\prime}(t)=x_{k_{0}}^{T}H(p_{k_{0}},tx_{k_{0}})x_{k_{0}}$, where $H$
denotes the Hessian of $\tilde{g}$. Since we are assuming that the smallest
eigenvalue of $H$ on $N_{\delta}(0)\times N_{\delta}(0)$ is greater than
$\lambda/2$, it follows that $L^{\prime\prime}(t)>|x_{k_{0}}|^{2}\lambda/2$
for all $t\in[0,1]$. Also, recall the assumption that $0$ is a critical point
of $\tilde{g}$; hence $L^{\prime}(0)=0$. From this, it follows that
$f(p_{k_{0}},x_{k_{0}})=L(1)>L(0)=f(p_{k_{0}},0)$. This contradicts our
original construction of the sequences $\\{p_{k}\\}$ and $\\{x_{k}\\}$, which
required that for all $k$, $f(p_{k},x_{k})<f(p_{k},0)$, thereby proving the
lemma. ∎

We are now able to proceed with the proof of Theorem 2.10.

###### Proof of Theorem 2.10.

Clearly, the set ${\mathcal{P}}_{n}$ is a $3n$-manifold, and the symmetric
group $\Sigma_{n}$ acts freely on ${\mathcal{P}}_{n}$ in the obvious way.
Furthermore, if we give ${\mathcal{P}}_{n}$ the standard Euclidean Riemannian
metric restricted from $\mathbb{R}^{3n}$, $\Sigma_{n}$ acts by isometries. It
follows that ${\mathcal{P}}_{n}/\Sigma_{n}$ inherits the structure of a flat
Riemannian manifold.

Let $\mu\in{\mathcal{M}}^{n}_{+}$; so we can write
$\mu=\sum_{i=1}^{n}a_{i}\delta_{u_{i}}$, where $\\{u_{1},\dots,u_{n}\\}$ is a
set of distinct elements of $S^{2}$ and each $a_{i}>0$. Clearly the map
$Q_{n}$, as defined in (18), induces a bijection
${\mathcal{P}}_{n}/\Sigma_{n}\to{\mathcal{M}}^{n}_{+}$.

To complete the proof of the first statement of Theorem 2.10, we just need to
prove that $Q_{n}$ maps a small neighborhood of each point $x$ in
${\mathcal{P}}_{n}$ isometrically to a small neighborhood of $Q_{n}(x)$ in
${\mathcal{M}}^{n}_{+}$. Indeed this will follow directly from the second
statement, i.e., from the fact that if two measures are supported by the same
number of points and have their support points and weights very close to each
other, then the optimal discrete semi-coupling between them is given by the
obvious diagonal matrices coming from the 1-1 correspondence between their
supports and weights.

Let $\mu$ and $\nu$ be two measures on $S^{2}$, both supported at precisely
$n$ points. Write $\mu=\sum_{i=1}^{n}a_{i}\delta_{u_{i}}$ and
$\nu=\sum_{i=1}^{n}b_{i}\delta_{v_{i}}$, where $u_{1},\dots,u_{n}$ are
distinct elements of $S^{2}$ and $a_{1},\dots,a_{n}>0$; similarly,
$v_{1},\dots,v_{n}$ are distinct elements of $S^{2}$ and
$b_{1},\dots,b_{n}>0$. Recall that a discrete semi-coupling from $\mu$ to
$\nu$ is defined to be a pair of $(n+1)\times(n+1)$ matrices $(A,B)$ (where
both indices in each matrix run from $0$ to $n$) satisfying the following
conditions:

  1. (1)

for all $i,j$, $A_{ij}\geq 0$ and $B_{ij}\geq 0$;

  2. (2)

for all $i=1,\dots,n$, $a_{i}=A_{i0}+\dots+A_{in}$;

  3. (3)

for all $j=0,\dots,n$, $A_{0j}=0$;

  4. (4)

for all $j=1,\dots,n$, $b_{j}=B_{0j}+\dots+B_{nj}$;

  5. (5)

for all $i=0,\dots,n$, $B_{i0}=0$.

We define a cost function $C$ on the set of all discrete semi-couplings from
$\mu$ to $\nu$ by

$C(A,B)=\sum_{i=0}^{m}\sum_{j=0}^{n}\left|\sqrt{B_{ij}}v_{j}-\sqrt{A_{ij}}u_{i}\right|^{2}$

An optimal discrete semi-coupling from $\mu$ to $\nu$ is defined to be a
discrete semi-coupling that minimizes $C$. Note that the set of discrete semi-
couplings (i.e., the domain of $C$) varies according to the particular
measures $\mu$ and $\nu$. To remedy this inconvenient fact, we define a
normalized discrete semi-coupling to be pair of matrices $(A,B)$ satisfying
the conditions:

  1. (1)

for all $i,j$, $A_{ij}\geq 0$ and $B_{ij}\geq 0$;

  2. (2)

for all $i=1,\dots,n$, $1=A_{i0}+\dots+A_{in}$;

  3. (3)

for each $j=0,\dots,n$, $A_{0j}=0$;

  4. (4)

for all $j=1,\dots,n$, $1=B_{0j}+\dots+B_{nj}$;

  5. (5)

for each $i=0,\dots,n$, $B_{i0}=0$.

We then define the normalized cost function by

$\bar{C}(A,B)=\sum_{i=0}^{n}\sum_{j=0}^{n}\left|\sqrt{b_{j}B_{ij}}v_{j}-\sqrt{a_{i}A_{ij}}u_{i}\right|^{2},$

where $a_{0}$ and $b_{0}$ should be assigned the value 0.

Note that we are just rescaling each row of $A$ and each column of $B$ to make
their sums 1, and then inserting the relevant scalars into the cost function
so as not to change its behavior. Define an optimal normalized discrete semi-
coupling from $\mu$ to $\nu$ to be a normalized discrete semi-coupling that
minimizes the normalized cost function. Note that the set of normalized
discrete semi-couplings no longer depends on the particular pair of measures
$\mu$ and $\nu$ (as long as they are both supported at precisely $n$ points).
However, the normalized cost function does depend on $\mu$ and $\nu$.
Henceforth, to make this dependence explicit, we write $\bar{C}_{\mu,\nu}$ to
denote the normalized cost function corresponding to the measures $\mu$ and
$\nu$.

We now make a further change of variables. For each normalized discrete semi-
coupling $(A,B)$, we define another pair of $(n+1)\times(n+1)$ matrices
$(P,Q)$ as follows: For each $i,j$, $P_{i,j}=\sqrt{A_{i,j}}$ and
$Q_{i,j}=\sqrt{B_{i,j}}$. If we write the cost function as a function of
$(P,Q)$, the form of the cost function becomes

$G_{\mu,\nu}(P,Q)=\sum_{i=0}^{n}\sum_{j=0}^{n}\left|Q_{i,j}\sqrt{b_{j}}v_{j}-P_{i,j}\sqrt{a_{i}}u_{i}\right|^{2}.$

Note that the corresponding domain of $G_{\mu,\nu}$ is the set of all pairs
$(P,Q)$ of $(n+1)\times(n+1)$ matrices satisfying the following:

  1. (1)

for all $i,j$, $P_{ij}\geq 0$ and $Q_{ij}\geq 0$;

  2. (2)

for all $i=1,\dots,n$, $1=P_{i0}^{2}+\dots+P_{in}^{2}$;

  3. (3)

for each $j=0,\dots,n$, $P_{0j}=0$;

  4. (4)

for all $j=1,\dots,n$, $1=Q_{0j}^{2}+\dots+Q_{nj}^{2}$;

  5. (5)

for each $i=0,\dots,n$, $Q_{i0}=0$.

Let ${\mathcal{T}_{0}}$ denote the domain of $G_{\mu,\nu}$; it is a product of
positive orthants of unit spheres, since each of rows 1 through $n$ of $P$,
and each of columns 1 through $n$ of Q is constrained to be in such an
orthant. (Row 0 of $P$ and column 0 of $Q$ are each required to be the zero
vector.) Also, note that we can remove the first condition on the domain of
$G_{\mu,\nu}$, thereby extending its domain to be a $2n$-fold product of
spheres, with each sphere having dimension $n$. Denote this extended domain by
${\mathcal{T}}$. Clearly, $G_{\mu,\nu}$ is defined (by the same formula) and
smooth on ${\mathcal{T}}$.

Now consider the case $\mu=\nu$. In this case, $G_{\mu,\mu}$ achieves a
minimum value of $0$ on ${\mathcal{T}_{0}}$, and this value is achieved only
at the point where

$P=Q=\begin{pmatrix}0&0\cr 0&I_{n}\end{pmatrix}.$

Denote this particular discrete semi-coupling by $(P_{0},Q_{0})$. In fact,
it’s clear that this value is also a minimum on all of ${\mathcal{T}}$ since
$G_{\mu,\mu}$ is a sum of non-negative terms. (On all of ${\mathcal{T}}$ there
are other points besides $(P_{0},Q_{0})$ where this minimum is achieved.)
$(P_{0},Q_{0})$ is a critical point of $G_{\mu,\mu}$ on ${\mathcal{T}}$, since
it is a point where the minimum is achieved. Furthermore, we will show that
the Hessian of $G_{\mu,\mu}$ at $(P_{0},Q_{0})$ is positive definite.

To see that $G_{\mu,\mu}$ has positive definite Hessian at $(P_{0},Q_{0})$,
reason as follows: Let $(P(t),Q(t))$, for $-a<t<a$, be a path in
${\mathcal{T}}$ such that $(P(0),Q(0))=(P_{0},Q_{0})$, and
$(P^{\prime}(0),Q^{\prime}(0))\neq(0,0)$. We will now show that
$G_{\mu,\mu}(P(t),Q(t))$ has positive second derivative at $t=0$, no matter
which such path is chosen. Note that no summand $G_{\mu,\mu}(P(t),Q(t))$ can
have negative second derivative, since $(P(0),Q(0))$ is a local minimum of
each summand. Hence, it suffices to show that just one summand of
$G_{\mu,\mu}(P(t),Q(t))$ has positive second derivative. Note that the
diagonal entries of $P^{\prime}(0)$ and $Q^{\prime}(0)$ are all zero, since
${\mathcal{T}}$ is a product of spheres. It follows that either
$P^{\prime}(0)$ or $Q^{\prime}(0)$ must have a nonzero entry in an off-
diagonal element. Thus, choose $i_{0}\neq j_{0}$, and assume that
$P^{\prime}_{i_{0},j_{0}}(0)\neq 0$. (The reasoning is the same for the case
$Q^{\prime}_{i_{0},j_{0}}(0)\neq 0$.) First consider the case $j_{0}>0$. The
corresponding summand of $G_{\mu,\mu}(P(t),Q(t))$ is

$\left|Q_{i_{0},j_{0}}(t)\sqrt{b_{j_{0}}}v_{j_{0}}-P_{i_{0},j_{0}}(t)\sqrt{a_{i_{0}}}u_{i_{0}}\right|^{2}.$

A straightforward computation shows that the second derivative of this summand
at $t=0$ is

$2\left|Q^{\prime}_{i_{0},j_{0}}(0)\sqrt{b_{j_{0}}}v_{j_{0}}-P^{\prime}_{i_{0},j_{0}}(0)\sqrt{a_{i_{0}}}u_{i_{0}}\right|^{2}.$

Since $i_{0}\neq j_{0}$, it follows that $u_{i_{0}}$ and $v_{j_{0}}$ are
linearly independent. Hence this second derivative is positive, since we are
assuming that $P^{\prime}_{i_{0},j_{0}}(0)\neq 0$. For the case $j_{0}=0$, we
know that $Q_{i_{0},j_{0}}(t)$ is constant at zero (since it stays in
${\mathcal{T}}$). Therefore $Q^{\prime}_{i_{0},j_{0}}(0)=0$, and it still
follows that the second derivative of this summand is positive. This proves
that $G_{\mu,\mu}(P(t),Q(t))$ has positive second derivative at $t=0$ for
every path $(P(t),Q(t))$ in ${\mathcal{T}}$ with $(P(0),Q(0))=(P_{0},Q_{0})$.
Hence the Hessian of $G_{\mu,\mu}$ at $(P_{0},Q_{0})$ is positive definite.

Claim: For all $\mu,\nu\in{\mathcal{M}}_{n}^{0}$, $(P_{0},Q_{0})$ is a
critical point of $G_{\mu,\nu}$ on ${\mathcal{T}}$.

The proof of this Claim is a straightforward computation. First we construct a
basis for the tangent space $T_{(P_{0},Q_{0})}{\mathcal{T}}$, using the fact
that ${\mathcal{T}}$ is a product of unit spheres. This basis is a union of
two sets. The first set consists of tangent vectors of the form $(E_{i,j},0)$,
where $E_{i,j}$ has the entry 1 in the $(i,j)$-th place and 0’s elsewhere, for
$1\leq i\leq n$, $0\leq j\leq n$, and $i\neq j$. The second set consists of
elements of the form $(0,E_{i,j})$, where $E_{i,j}$ has the entry $1$ in the
$(i,j)$-th space, and 0’s elsewhere, for $0\leq i\leq n$, $1\leq j\leq n$, and
$i\neq j$. It is then a trivial calculation to see that the directional
derivative of $G_{\mu,\nu}$ at $(P_{0},Q_{0})$ in the direction of any of
these tangent vectors vanishes. This proves the claim.

Let $\mu_{0}\in{\mathcal{M}}_{n}^{0}$. We need to prove that there exists an
$\epsilon>0$ such that for all $(\mu,\nu)$ in an $\epsilon$-neighborhood of
$(\mu_{0},\mu_{0})$ in ${\mathcal{M}}_{n}^{0}\times{\mathcal{M}}_{n}^{0}$, the
minimum value of $G_{\mu,\nu}$ on ${\mathcal{T}}_{0}$ is achieved at
$(P_{0},Q_{0})$.

Define a function
$H:({\mathcal{P}}^{n}\times{\mathcal{P}}^{n})\times{\mathcal{T}}\to\mathbb{R}$
by

$H((x,y),(P,Q))=G_{Q_{n}(x),Q_{n}(y)}(P,Q).$

Clearly, $H$ is smooth. We are now in a position to apply Lemma 2.3. The
compact set ${\mathcal{T}}_{0}$ plays the role of $C$ in the Lemma. We have
shown that $(P_{0},Q_{0})$ is a critical point of $(P,Q)\mapsto
H((x,y),(P,Q))$, for all $x,y\in{\mathcal{P}}_{n}$, and that for any
$x_{0}\in{\mathcal{P}}_{n}$, $(P_{0},Q_{0})$ is the location of the unique
minimum of $(P,Q)\mapsto H((x_{0},x_{0}),(P,Q))$ on ${\mathcal{T}_{0}}$.
Having also verified the condition on the Hessian, we can conclude from Lemma
2.3 that there exists an $\epsilon>0$ such that for all $(x,y)$ with
$|x-x_{0}|<\epsilon$ and $|y-x_{0}|<\epsilon$, $H((x,y),(\cdot,\cdot))$
attains its maximum (on the domain ${\mathcal{T}}_{0}$) at the point
$(P_{0},Q_{0})$. Thus the proof of the second statement follows, which implies
at the same time the first statement. ∎

###  2.4. A convex optimization problem

In this section we will study the class of optimization problems, that consist
of maximizing the function

$\displaystyle F:\mathcal{A}$ $\displaystyle\to\mathbb{R}\qquad\text{ defined
via }$ (19) $\displaystyle(A,B)$
$\displaystyle\mapsto\sum_{i=1}^{m}\sum_{j=1}^{n}\sqrt{A_{ij}B_{ij}}\Omega_{ij}.$
(20)

Here $\Omega_{ij}$ is any given weight matrix. Recall that the tuples
$(A,B)\in\mathcal{A}$ are subject to the constraints (a)-(e). In the main
result of this part we will explicitly construct a sequence of semi-couplings
that converges to an optimizer of $F$.

Our motivation for studying this class of optimization problems stems from the
fact that for a particular choice of $\Omega$ it is equivalent to calculating
the WFR-distance between two finitely supported measures. This will lead
directly to an efficient algorithm to numerically calculate this distance.
Note, that most existing methods for estimating the WFR-distance solve an
entropy regularized problem. Our proposed solution instead converges to the
WFR distance by performing an optimization on the polytope of discrete semi-
couplings.

To see this connection between the function $F$ and the WFR-distance let
$\mu=\sum\limits_{i=1}^{m}a_{i}\delta_{u_{i}}$ and
$\nu=\sum\limits_{i=1}^{n}b_{i}\delta_{v_{i}}$ be two finitely supported,
finite measures. Recall from Remark 2.2 that the distance from $\mu$ to $\nu$
can be written as

$\operatorname{WFR}(\mu,\nu)=\left(\sum_{i=1}^{m}a_{i}+\sum_{j=1}^{n}b_{j}-2\sup_{(A,B)\in\mathcal{A}}\sum_{i=1}^{m}\sum_{j=1}^{n}\sqrt{A_{ij}B_{ij}}(u_{i}\cdot
v_{j})\right)^{1/2}.$ (21)

Thus computing $\operatorname{WFR}(\mu,\nu)$ is equivalent to finding
$(A,B)\in\mathcal{A}$ that achieves the supremum

$\sup_{(A,B)\in\mathcal{A}}\sum_{i=1}^{m}\sum_{j=1}^{n}\sqrt{A_{ij}B_{ij}}(u_{i}\cdot
v_{j}),$ (22)

which corresponds to the function $F$ with $\Omega$ being given by
$\Omega_{ij}=(u_{i}\cdot v_{j})$.

######  Remark 2.13 (Generalizations of the Algorithm).

In [14], the Wasserstein Fisher Rao distance is formulated on a domain $M$
with a parameter $\rho\in(0,\infty)$ as

$\frac{1}{\rho^{2}}\operatorname{WFR}^{2}_{\rho}(\mu,\nu)=\mu(M)+\nu(M)-2\sup\limits_{(\gamma_{0},\gamma_{1})\in\Gamma(\mu,\nu)}\int_{d(x,y)<\pi}\cos(d(x,y)/\rho)d(\sqrt{\gamma_{0}\gamma_{1}})(x,y)$
(23)

where
$\sqrt{\gamma_{0}\gamma_{1}}=(\frac{\gamma_{0}}{\gamma}\frac{\gamma_{1}}{\gamma})^{1/2}$
for $\gamma$ such that $\gamma_{0},\gamma_{1}\ll\gamma$. Another popular
distance is the so-called Gaussian Hellinger-Kantorovich distance, which can
be expressed with a parameter $\rho\in(0,\infty)$ as

$\frac{1}{\rho^{2}}\operatorname{GHK}^{2}_{\rho}(\mu,\nu)=\mu(M)+\nu(M)-2\sup\limits_{(\gamma_{0},\gamma_{1})\in\Gamma(\mu,\nu)}\int
e^{\frac{-d^{2}(x,y)}{2\rho}}d(\sqrt{\gamma_{0}\gamma_{1}})(x,y).$ (24)

For finitely supported measures $\mu,\nu$ these distances can be generically
expressed as

$\frac{1}{\rho^{2}}\operatorname{Dist}^{2}_{\rho}(\mu,\nu)=\sum_{i=1}^{m}a_{i}+\sum_{j=1}^{n}b_{j}-2\sup_{(A,B)\in\mathcal{A}}\sum_{i=1}^{m}\sum_{j=1}^{n}\sqrt{A_{ij}B_{ij}}\Omega_{ij}$
(25)

where

$\operatorname{Dist}_{\rho}(\mu,\nu)=\begin{cases}\operatorname{WFR}_{\rho}(\mu,\nu)&\text{
if }\Omega_{ij}=\cos(d(u_{i},v_{j})/\rho)\\\
\operatorname{GHK}_{\rho}(\mu,\nu)&\text{ if
}\Omega_{ij}=e^{\frac{-d^{2}(u_{i},v_{j})}{2\rho}}\end{cases}.$ (26)

As our main result is that the SRNF can be written as the pullback of the WFR
distance with $\rho=1$ for measures on $M=S^{2}$, we developed the algorithm
with this particular case in mind. However, our algorithm does not depend on
the particular choice of $\Omega$ and thus it could be also used for these
more general situations.

First we will show that the optimum of $F$ will be obtained on a set such that
mass is never transported between supports $u_{i}$ and $v_{j}$ where
$\Omega_{ij}$ is negative. Moreover, we show that $F$ obtains its optimum when
all of the mass of a given support is created/destroyed if and only if none of
it can be transported for a suitably small enough cost. In terms of the
discrete semi-couplings this is represented by non-zero entries in the zeroth
column of $A$ and the zeroth row of $B$. Therefore we define the subset
$\tilde{\mathcal{A}}$ of $\mathcal{A}$ of all semi-couplings
$(A,B)\in\mathcal{A}$ satisfying

  1. (1)

$A_{ij}=B_{ij}=0\operatorname{whenever}\Omega_{ij}\leq 0,$

  2. (2)

$B_{0j}=0\text{ if there exists }j\text{ such that }\Omega_{ij}>0$, and

  3. (3)

$A_{i0}=0\text{ if there exists }i\text{ such that }\Omega_{ij}>0$.

We will now show that the value function $F$ obtains its maximum on
$\tilde{\mathcal{A}}$ and that $F$ is concave when restricted to this subset.

######  Lemma 2.14.

The value function $F$ obtains its maximum on $\tilde{\mathcal{A}}$. Moreover,
$\tilde{\mathcal{A}}$ is a convex, compact subset of $\mathbb{R}^{mn+m+n+1}$
and $F$ is concave when restricted to $\tilde{\mathcal{A}}$.

###### Proof.

To show this statement, let $(A,B)\in\mathcal{A}$. We will construct a new
semi-coupling $(\tilde{A},\tilde{B})\in\tilde{\mathcal{A}}$ such that
$F(A,B)\leq F(\tilde{A},\tilde{B})$. Begin by letting
$(\tilde{A},\tilde{B})=(A,B)$. Next for each $i\in\\{1,...,m\\}$ define

$k_{i}=\begin{cases}\operatorname*{arg\,max}\limits_{k\in\\{1,...,n\\}}\Omega_{ik}&\text{if
}\max\limits_{k\in\\{1,...,n\\}}\Omega_{ij}>0\\\
0&\text{otherwise}\end{cases}$ (27)

and for each $j=\\{1,...,n\\}$ define

$l_{j}=\begin{cases}\operatorname*{arg\,max}\limits_{l\in\\{1,...,m\\}}\Omega_{lj}&\text{if
}\max\limits_{l\in\\{1,...,m\\}}\Omega_{lj}>0\\\
0&\text{otherwise}\end{cases}.$ (28)

Step 1: If $\Omega_{ij}<0$, we modify $\tilde{A}_{ij}$ and $\tilde{B}_{ij}$ as
follows:

  * •

Replace $\tilde{A}_{ik_{i}}$ by $A_{ik_{i}}+A_{ij}$ and replace
$\tilde{A}_{ij}$ by $0$.

  * •

Replace $\tilde{B}_{l_{j}j}$ by $B_{l_{j}j}+B_{ij}$ and replace
$\tilde{B}_{ij}$ by $0$.

It is clear that the new point satisfies property (1). In the next stage of
our modification we will ensure that our semi-coupling also satisfies (2) and
(3).

Step 2:

  * •

If $k_{i}\neq 0$, then replace $\tilde{A}_{ik_{i}}$ by $A_{ik_{i}}+A_{i0}$ and
$\tilde{A}_{i0}$ by $0$.

  * •

If $l_{j}\neq 0$, then replace $\tilde{B}_{l_{j}j}$ by $B_{l_{j}j}+B_{0j}$ and
$\tilde{B}_{0j}$ by $0$.

It is clear that $(\tilde{A},\tilde{B})$ is an element of
$\tilde{\mathcal{A}}$ and that $(\tilde{A},\tilde{B})=(A,B)$ if $(A,B)$ was
already in $\tilde{\mathcal{A}}$. Moreover, for each $i\in\\{1,...,m\\}$ and
$j=\\{1,...,n\\}$,

$A_{ij}B_{ij}\Omega_{ij}\leq\tilde{A}_{ij}\tilde{B}_{ij}\Omega_{ij}$ (29)

and thus $F(A,B)\leq F(\tilde{A},\tilde{B})$. In fact, it is clear that the
value of $F$ steadily increases along the line from the original point to the
new point. This contradicts the assumption that $F$ attained a local maximum
at the original point. Thus $F$ obtains its maximum on the subset
$\tilde{\mathcal{A}}$. A parametrized straight line in $\tilde{\mathcal{A}}$
will be of the form $\\{p_{ij}t+q_{ij},y_{ij}t+z_{ij}\\}$. If we restrict $F$
to such a line, it will be a linear combination of functions of the form
$\sqrt{(pt+q)(yt+z)}$, with positive coefficients. It is easy to verify that a
function of this form always has a second derivative $\leq 0$ for all values
of $t$ for which it is defined (i.e., for all values of $t$ for which the
quantity under the square root sign is $\geq 0$). It follows that along any
line in the domain, either $F$ has at most one local maximum, or it is
constant along that line. Therefore, if $F$ has two local maxima on its
domain, then $F$ must be constant on the entire line through these two maxima.
Thus, $F$ is concave when restricted to $\tilde{\mathcal{A}}$. ∎

Thus, we have reformulated our optimization problem as finding the maximum of
a concave function over a convex set (or equivalently the minimum of a convex
function over a convex set). Next, we will introduce two operators on
$\tilde{\mathcal{A}}$: First, let

$T_{1}:\tilde{\mathcal{A}}\to\tilde{\mathcal{A}}\text{ is defined via
}(A,B)\mapsto(E,B)\text{ where }$ (30) $E_{ij}=\begin{cases}A_{ij}&\text{if
$i=0$ or $j=0$}\\\
\frac{a_{i}B_{ij}\Omega_{ij}^{2}}{\sum\limits_{k=1}^{n}B_{ik}\Omega_{ik}^{2}}&\text{if
}\sum\limits_{k=1}^{n}B_{ik}\Omega_{ik}^{2}>0\\\
0&\text{otherwise}\end{cases}.$ (31)

Similarly let

$T_{2}:\tilde{\mathcal{A}}\to\tilde{\mathcal{A}}\text{ be defined via
}(A,B)\mapsto(A,E)\text{ where }$ (32) $E_{ij}=\begin{cases}B_{ij}&\text{if
$i=0$ or $j=0$}\\\
\frac{b_{j}A_{ij}\Omega_{ij}^{2}}{\sum\limits_{k=1}^{m}A_{kj}\Omega_{kj}^{2}}&\text{if
}\sum\limits_{k=1}^{m}A_{kj}\Omega_{kj}^{2}>0\\\ 0&\text{otherwise
}\end{cases}.$ (33)

We will define a sequence of semi-couplings recursively by initializing the
sequence at some $(A,B)_{0}$ in the interior of $\tilde{\mathcal{A}}$ and
iteratively updating the two components of this semi-coupling via $T_{1}$ and
$T_{2}$. The primary goal of this section will be to show that any limit point
of this sequence is a maximizer of $F$, which is formalized in the following
theorem:

######  Theorem 2.15.

Define a sequence of semi-couplings via $(A,B)_{k+1}:=T_{2}(T_{1}((A,B)_{k}))$
where $(A,B)_{0}$ is in arbitrary chosen initialization from the interior of
$\tilde{\mathcal{A}}$. Then the real-valued sequence given by $F((A,B)_{k})$
converges to the optimum of $F$, i.e.,

$\lim\limits_{k\to\infty}F((A,B)_{k})=\operatorname*{arg\,max}\limits_{(A,B)\in\tilde{\mathcal{A}}}F(A,B).$
(34)

The key ingredient for the proof of this statement is the observation that
applying $T_{1}$ maximizes $F$ when we fix the second component and similarly
for $T_{2}$ when fixing the first component. This will allow us to apply
results from coordinate descent analysis to obtain the proof of the above
theorem. This will be the content of the following lemma:

######  Lemma 2.16.

The operators $T_{i}$ restrict to operators on the interior of
$\tilde{\mathcal{A}}$, i.e.,
$T_{i}(\tilde{\mathcal{A}}^{\mathrm{o}})\subset\tilde{\mathcal{A}}^{\mathrm{o}}$.
For a fixed $(A,B)$ in the interior $\tilde{\mathcal{A}}$ let
$\tilde{\mathcal{A}}_{(-,B)}$ be the space of semi-couplings where the second
factor is equal to $B$ and $\tilde{\mathcal{A}}_{(A,-)}$ be the space of semi-
couplings where the first factor is equal to $A$. Thus,

  1. (1)

$T_{1}(A,B)$ uniquely attains
$\sup\limits_{(A^{{}^{\prime}},B)\in\tilde{\mathcal{A}}_{(-,B)}}F(A^{{}^{\prime}},B)$

  2. (2)

$T_{2}(A,B)$ uniquely attains
$\sup\limits_{(A,B^{{}^{\prime}})\in\tilde{\mathcal{A}}_{(A,-)}}F(A,B^{{}^{\prime}})$

###### Proof.

Let$(A,B)$ be in the interior $\tilde{\mathcal{A}}$. Thus, for any $i=1,...,m$
and $j=1,...n$ such that $\Omega_{ij}>0$ we have that $A_{ij}>0$ and
$B_{ij}>0$. We need to show that for any semi-coupling such that the second
matrix is equal to $B$, $F(A^{{}^{\prime}},B)\leq F\left(T_{1}(A,B)\right)$
where equality holds if and only if $(A^{{}^{\prime}},B)=T_{1}(A,B)$.

Therefore, let $(A^{{}^{\prime}},B)$ be an arbitrary semi-coupling such that
the second matrix is equal to $B$ and let $(E,B)=T_{1}(A^{{}^{\prime}},B)$.
For $1\leq i\leq m$ we let

$\rho_{i}=\sum\limits_{j=1}^{n}B_{ij}\Omega_{ij}^{2}\;.$ (35)

Consider the case where $\rho_{i}=0$. Since $(A,B)$ is in
$\tilde{\mathcal{A}}$, $B_{ij}\geq 0$ and $B_{ij}=0$ if and only if
$\Omega_{ij}^{2}=0$. Thus, $B_{ij}=0$ for all $1\leq j\leq n$. Thus,

$\displaystyle\sum\limits_{j=1}^{n}\sqrt{A^{{}^{\prime}}_{ij}B_{ij}}\Omega_{ij}=\sum\limits_{j=1}^{n}\sqrt{E_{ij}B_{ij}}\Omega_{ij}=0,$
(36)

It remains to prove the statement for the case that $\rho_{i}>0$. By Cauchy-
Schwarz we have

$\displaystyle\sum\limits_{j=1}^{n}\sqrt{A^{{}^{\prime}}_{ij}B_{ij}}\Omega_{ij}$
$\displaystyle\leq\sqrt{\left(\sum\limits_{j=1}^{n}A^{{}^{\prime}}_{ij}\right)\left(\sum\limits_{j=1}^{n}B_{ij}\Omega_{ij}^{2}\right)}=\sqrt{a_{i}\sum\limits_{j=1}^{n}B_{ij}\Omega_{ij}^{2}}$
(37)
$\displaystyle=\frac{\sqrt{a_{i}}\sum\limits_{j=1}^{n}B_{ij}\Omega_{ij}^{2}}{\sqrt{\left(\sum\limits_{j=1}^{n}B_{ij}\Omega_{ij}^{2}\right)}}=\sum\limits_{j=1}^{n}\sqrt{\frac{a_{i}B_{ij}\Omega_{ij}^{2}}{\sum\limits_{k=1}^{n}B_{ik}\Omega_{ik}^{2}}B_{ij}}\Omega_{ij}=\sum\limits_{j=1}^{n}\sqrt{E_{ij}B_{ij}}\Omega_{ij}.$
(38)

Therefore,

$\displaystyle F(A^{{}^{\prime}},B)$
$\displaystyle=\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{n}\sqrt{A^{{}^{\prime}}_{ij}B_{ij}}\Omega_{ij}\leq\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{n}\sqrt{E_{ij}B_{ij}}\Omega_{ij}=F\left(T_{1}(A,B)\right).$
(39)

Note that this inequality is strict unless for each $1\leq i\leq m$ such that

$\sum\limits_{k=1}^{n}B_{ik}\Omega_{ik}^{2}>0,$ (40)

$\\{A^{{}^{\prime}}_{ij}\\}$ is a scalar multiple of
$\\{B_{ij}\Omega_{ij}^{2}\\}$. Since we know
$\sum\limits_{j=1}^{n}A^{{}^{\prime}}_{ij}=a_{i}$, equality holds if and only
if

$\\{A^{{}^{\prime}}_{ij}\\}=\left\\{\frac{a_{i}B_{ij}\Omega_{ij}^{2}}{\sum_{k=1}^{n}B_{ik}\Omega_{ik}^{2}}\right\\}=\\{E_{ij}\\}.$
(41)

The proof of (2) follows by a symmetric argument. ∎

Finally, we need to observe the effect of $T_{1}$ and $T_{2}$ on the value of
$F$ at the limit points of our sequence. This is not immediate because our
previous results require that our semi-coupling is in the interior of
$\tilde{\mathcal{A}}$, but the limit point of our sequence could lie on the
boundary of $\tilde{\mathcal{A}}$.

######  Lemma 2.17.

If $(\overline{A},\overline{B})$ is a limit point of $(A,B)_{k}$, then

$F(T_{1}(\overline{A},\overline{B}))=F(T_{2}(T_{1}(\overline{A},\overline{B})))=F(\overline{A},\overline{B}).$
(42)

###### Proof.

Let $(\overline{A},\overline{B})$ be a limit point of $(A,B)_{k}$ and
$(A,B)_{k_{j}}$ be a subsequence that converges to
$(\overline{A},\overline{B})$. Let $\epsilon>0$. Since $F\circ T_{1}$ is
continuous, there exists a $\delta>0$ such that if
$\|(A^{{}^{\prime}},B^{{}^{\prime}})-(\overline{A},\overline{B})\|<\delta$
then
$|F(T_{1}(A^{{}^{\prime}},B^{{}^{\prime}}))-F(T_{1}(\overline{A},\overline{B}))|<\epsilon/3$.

Since $(A,B)_{k_{j}}$ converges to $(\overline{A},\overline{B})$, it follows
from the continuity of $F$ that $F((A,B)_{k_{j}})$ converges to
$F(\overline{A},\overline{B})$. Thus, there exists $k_{j}$ such that

$\|(A,B)_{k_{j}}-(\overline{A},\overline{B})\|<\delta\text{ and
}|F((A,B)_{k_{j}})-F(\overline{A},\overline{B})|<\epsilon/3.$ (43)

Recall, $F((A,B)_{k_{j}})\leq F(T_{1}((A,B)_{k_{j}}))\leq
F(\overline{A},\overline{B})$. Thus,

$|F(T_{1}((A,B)_{k_{j}}\\}))-F((A,B)_{k_{j}})|\leq|F((A,B)_{k_{j}})-F(\overline{A},\overline{B})|<\epsilon/3.$
(44)

Therefore,

$\displaystyle|F(T_{1}(\overline{A},\overline{B}))-F(\overline{A},\overline{B})|\leq|F(T_{1}(\overline{A},\overline{B}))-F(T_{1}((A,B)_{k_{j}}))|$
(45)
$\displaystyle\qquad\qquad+|F(T_{1}((A,B)_{k_{j}}))-F((A,B)_{k_{j}})|+|F((A,B)_{k_{j}})-F(\overline{A},\overline{B})|$
(46) $\displaystyle\qquad<\epsilon/3+\epsilon/3+\epsilon/3=\epsilon$ (47)

The proof that
$F(T_{2}(T_{1}(\overline{A},\overline{B})))=F(\overline{A},\overline{B})$
follows by a similar argument. ∎

We can now proceed with the proof of Theorem 2.15.

###### Proof of Theorem 2.15.

Note that for every $k\in\mathbb{N}$, $F((A,B)_{k+1})\geq F((A,B)_{k})$. Since
$F$ is bounded above. The sequence defined by $F((A,B)_{k})$ converges.

Let $(\overline{A},\overline{B})$ be a limit point of $(A,B)_{k}$. From Lemma
2.4, we have $F(T_{1}((A,B)_{k}))\geq F(A^{{}^{\prime}},B_{k}),$ for each
$(A^{{}^{\prime}},B_{k})\in\tilde{\mathcal{A}}$ where $B_{k}$ is fixed. Taking
the limit of this inequality as $k\to\infty$, we obtain
$F(T_{1}(\overline{A},\overline{B}))\geq F(A^{{}^{\prime}},\overline{B}),$ for
each $(A^{{}^{\prime}},\overline{B})\in\tilde{\mathcal{A}}$ where
$\overline{B}$ is fixed. By Lemma 2.4,
$F(T_{1}(\overline{A},\overline{B}))=F(\overline{A},\overline{B})$.

By the optimality condition (Prop. 2.1.2 in Section 2.1 of [6]) on convex
sets,
$\nabla_{1}F(\overline{A},\overline{B})(\overline{A}-A^{{}^{\prime}})\leq 0$
for all possible $A^{{}^{\prime}}$. Here $\nabla_{1}F$ is the gradient of $F$
with respect to the first block. A similar argument shows that,
$\nabla_{2}F(\overline{A},\overline{B})(\overline{B}-B^{{}^{\prime}})\leq 0$
for all possible $B^{{}^{\prime}}$. Here $\nabla_{2}F$ is the gradient of $F$
with respect to the second block.

Combining the inequalities and using the product structure of
$\tilde{\mathcal{A}}$, we obtain the statement $\nabla
F(\overline{A},\overline{B})((\overline{A},\overline{B})-(A^{{}^{\prime}},B^{{}^{\prime}}))\leq
0$ for all $(A^{{}^{\prime}},B^{{}^{\prime}})$ in $\tilde{\mathcal{A}}$. Thus
by the optimality condition of concave functions on convex sets,
$(\overline{A},\overline{B})=\operatorname*{arg\,max}\limits_{(A,B)\in\tilde{\mathcal{A}}}F(A,B).$
∎

###  2.5. Implementation and Experiments

The theory developed in the previous section directly gives rise to an
algorithm for computing the WFR distance. We outline this procedure in
Algorithm 1. A PyTorch implementation of our algorithm is open source
available at

https://github.com/emmanuel-hartman/WassersteinFisherRaoDistance

Unlike the Sinkhorn-type methods proposed in [13] and implemented in [17], our
method computes the Wasserstein-Fisher-Rao distance without any entropic
regularization.In this implementation, we assume that for all $j$ there exists
$i$ such that $\Omega_{ij}>0$ and for all $i$ there exists $j$ such that
$\Omega_{ij}>0$. As such, the zeroth rows and columns of the optimal discrete
semi-coupling will have all zero entries, so we omit them from our
implementation. Our reason for making this assumption on $\Omega_{ij}$ is that
our main application of this algorithm is to compute SRNF distances between
closed surfaces. For closed surfaces, this assumption will always be valid.
(See Section 3 of this paper for more details.)

We implemented Alg. 1 using PyTorch and perform the computations on the GPU.
The main operations in our algorithm consist of element-wise matrix
multiplication, which leads to a quadratic complexity. This can be also
observed in the computation times in Table 1. As the obtained semi-coupling
matrices are usually sparsely populated, we expect that an implementation
utilizing sparse matrix data types could significantly improve the performance
of the implementation.

Algorithm 1 Calculate $\operatorname{WFR}$

procedure WFR_distance($a,b,u,v,\epsilon$)

$\Omega_{ij}\leftarrow(u_{i}\cdot v_{j})$ if $i,j\geq 1$ and 1 otherwise

$A,B\leftarrow\text{some intitialization in }\tilde{A}^{\mathrm{o}}$,
$err\leftarrow\infty$

$c\leftarrow\sum_{i,j}\sqrt{A_{ij}\cdot B_{ij}}\cdot\Omega_{ij}$

while $\ err\geq\epsilon$

$(A,B)\leftarrow\textsc{WFR\\_iteration}(a,b,\Omega,A,B)$

$c_{\text{prev}}\leftarrow c$

$c\leftarrow\sum_{i,j}\sqrt{A_{ij}\cdot B_{ij}}\cdot\Omega_{ij}$

$err\leftarrow(c-c_{\text{prev}})/c$

return $\left(\sum_{i}a_{i}+\sum_{j}b_{j}-2\sum_{i,j}\sqrt{A_{ij}\cdot
B_{ij}}\cdot\Omega_{ij}\right)^{1/2}$

procedure WFR_iteration($a,b,\Omega,A,B$)

$B^{\prime}_{ij}\leftarrow B_{ij}\cdot\Omega_{ij}\cdot\Omega_{ij}$

$R_{i}\leftarrow a_{i}/\sum_{j}B^{\prime}_{ij}$

$A_{ij}\leftarrow B^{\prime}_{ij}\cdot R_{i}$

$A^{\prime}_{ij}\leftarrow A_{ij}\cdot\Omega_{ij}\cdot\Omega_{ij}$

$C_{j}\leftarrow b_{j}/\sum_{i}A^{\prime}_{ij}$

$B_{ij}\leftarrow A^{\prime}_{ij}\cdot C_{j}$

return (A,B)

To quantify the performance of our algorithm we perform experiments comparing
our method with entropy regularized methods using a small regularization
parameter. Therefore we construct random pairs of finitely supported measures
with a fixed number of support points. We then calculate the distance using
both our algorithm and the unbalanced Sinkhorn algorithm [13], where we choose
the regularization parameter to be $1\times 10^{-3}$ (this was the smallest
value that led to a stable performance across all experiments). Therefore we
will first discuss the relation between the variables of these two algorithms.
Due to [14], we have

$\operatorname{WFR}^{2}(\mu,\nu)=\inf_{\gamma\in\mathcal{M}(M^{2})}G(\gamma)$
(48)

where

$G(\gamma)=\operatorname{KL}((\operatorname{Proj}_{0})_{\\#}\gamma,\mu)+\operatorname{KL}((\operatorname{Proj}_{1})_{\\#}\gamma,\nu)-\int_{M^{2}}\log(\cos^{2}(d(x,y)\wedge(\pi/2)))d\gamma(x,y).$
(49)

Forthcoming results of by Gallouët, Ghezzi and Vialard [18], show that given a
semi-coupling $(\gamma_{1},\gamma_{2})$ one can produce an optimal transport
plan,

$(\gamma_{1},\gamma_{2})\mapsto\gamma=\sqrt{\gamma_{1}\gamma_{2}}e^{-c/2}\text{
where }c=-\log(\cos^{2}(d(x,y)\wedge(\pi/2)))$ (50)

such that

$G(\gamma)=\int_{M^{2}}\left|\sqrt{\frac{\gamma_{1}}{\gamma}(u,v)}u-\sqrt{\frac{\gamma_{2}}{\gamma}(u,v)}v\right|^{2}d\gamma(u,v)$
(51)

for $\gamma$ such that $\gamma_{1},\gamma_{2}\ll\gamma$. Therefore, we can
compare the transport plans produced by both methods. The Sinkhorn type
algorithm proposed in [13] then solves a regularized optimization problem
given by

$\inf\limits_{\gamma\in\mathcal{M}(M^{2})}G(\gamma)+\lambda\operatorname{Ent}(\gamma)$
(52)

where $\lambda$ is the regularization parameter and
$\operatorname{Ent}(\cdot)$ is an entropic regularization term. To obtain a
fair comparison of the solutions obtained with our algorithm and the solutions
obtained with the Sinkhorn algorithm, we disregarded the final entropy of the
Sinkhorn solution and only compared the corresponding transport costs. The
distances that resulted from our algorithm were consistently smaller (and
consequently more precise) across all experiments, which can be seen in Table
1, where we report the mean errors and variances of the Sinkhorn algorithm as
compared to the obtained distance using our algorithm. For each number of
support points, the relative errors were calculated by repeating the
experiments 100 times. As one can see in this table our method produces
significantly more accurate distances and without having to choose an entropic
regularization parameter.

We also report the corresponding mean computation times for these experiments
using the Sinkhorn algorithm. As the implementation of [17] does not utilize
the GPU, we ported their implementation to PyTorch to be able to have a fair
comparison of the two methods. Both algorithms are run on an Intel 3.2 GHz CPU
with a Gigabyte GeForce GTX 2070 1620 MHz GPU. We used a maximum of 2000
iterations, but we usually observed a much faster convergence. While the GPU
Sinkhorn implementation seems to scale better for larger numbers of support
points, we believe that the above mentioned adaptions would lead to a similar
complexity for Algorithm 1. In addition, as one can see in Figure 1, the
Sinkhorn algorithm has a significantly faster convergence and will thus lead
to a faster computation time. We want to emphasize, however, that our
algorithm solves the exact problem, while the Sinkhorn algorithm only tackles
a regularized problem. This is also mirrored by the fact that in all our
experiments the solutions obtained with our algorithm have a lower distance as
compared to corresponding Sinkhorn solutions, cf. Table 1.

Support | WFR Distance | Sinkhorn Error | Timing in seconds
---|---|---|---
points | Alg. 1 | Sink. | Mean | Variance | Alg. 1 | Sink. (GPU) | Sink. [17]
128 | 24.045 | 25.040 | 3.883% | 1.835% | 0.100 | 1.348 | 0.101
256 | 25.824 | 26.951 | 4.445% | 2.909% | 0.120 | 1.413 | 0.135
512 | 27.526 | 29.532 | 7.303% | 3.071% | 0.179 | 1.470 | 0.219
1024 | 29.724 | 32.302 | 8.673% | 3.490% | 0.398 | 1.523 | 0.509
2048 | 31.728 | 34.533 | 8.929% | 3.976% | 1.372 | 1.901 | 4.506
4096 | 35.061 | 38.027 | 8.382% | 3.699% | 5.569 | 3.152 | 20.323
8192 | 38.385 | 40.426 | 5.250% | 3.154% | 23.472 | 8.732 | 81.905
Table 1. Comparison of the Sinkhorn algorithm of [13] and our Algorithm. For
each number of support points theses results were obtained by calculating the
WFR distance between 100 pairs of randomly chosen measures.

Figure 1. Left: Plot of $||(A,B)-(A^{*},B^{*})||_{L^{2}}$ for 2000 iterations
on pairs of measures with $n$ supports: 256 (Orange), 512 (Green), 1024
(Purple), 2048 (Olive), 4096 (Light Blue), 8192 (Grey). Right: Plot of
$||\gamma-\gamma^{*}||_{2}$ for 2000 iterations of our algorithm (Orange) and
the Sinkhorn algorithm (Blue). The optimal semi-couplings, $(A^{*},B^{*})$,
and optimal plans $\gamma^{*}$ are computed with 10,000 iterations of the
respective algorithms.

##  3\. The SRNF shape metric as an unbalanced transport problem

In this section, we will present the main result of our article: the
interpretation of the SRNF shape distance as an unbalanced OMT problem. Our
result will then allow us to compute the SRNF distance using the algorithm
introduced in Section 2.4 and we will use this to present several numerical
examples at the end of the section.

###  3.1. Shape spaces of surfaces

In all of this section let $M$ be a smooth, connected, compact, oriented
Riemannian 2-dimensional manifold with or without boundary. In addition to the
smooth structure on $M$, we will be also interested in a piecewise linear
structure on it. By [60] any such $M$ indeed admits a Whitehead PL structure.
That is, there exists a polyhedral surface $K$ in $\mathbb{R}^{3}$ and a
homeomorphism called a triangulation $\sigma:K\to M$ such that $\sigma$ is
differentiable with injective differential on each face of $K$.

We denote the space of all Lipschitz immersions of $M$ into $\mathbb{R}^{3}$
by $\operatorname{Imm}_{\operatorname{Lip}}(M,\mathbb{R}^{3})$, i.e.,

$\displaystyle\operatorname{Imm}_{\operatorname{Lip}}(M,\mathbb{R}^{3})=\\{f\in
W^{1,\infty}(M,\mathbb{R}^{3}):Tf\text{ is inj. a.e.}\\}\;.$ (53)

The reason for considering immersions of the Lipschitz class is that this
space has two important subsets: the space of smooth immersions from $M$ to
$\mathbb{R}^{3}$ and the space of functions
$f\in\operatorname{Imm}_{\operatorname{Lip}}(M,\mathbb{R}^{3})$ that are PL
with respect to the given PL structure on $M$.

As we are interested in unparametrized surfaces, we have to factor out the
action of the group of diffeomorphisms. In the context of Lipschitz immersions
the natural group of reparametrizations for us to consider is the group of all
orientation preserving, bi-Lipschitz diffeomorphisms:

$\Gamma_{\operatorname{Lip}}(M)=\\{\gamma\in
W^{1,\infty}(M,M):\;\gamma^{-1}\in W^{1,\infty}(M,M),|D\gamma|>0\text{
a.e.}\\},$ (54)

where $|D\gamma|$ denotes the Jacobian determinant of $\gamma$, which is well-
defined as $D\gamma\in L^{\infty}$. For reasons, that will become clear later
in Section 3.2, we will also consider two subsets of
$\Gamma_{\operatorname{Lip}}(M)$, namely the group of smooth, orientation
preserving diffeomorphisms

$\displaystyle\Gamma_{\infty}(M)=\\{\gamma\in C^{\infty}(M,M):\;\gamma^{-1}\in
C^{\infty}(M,M),|D\gamma|>0\\},$ (55)

and the set of $PL$-homeomorphisms. To define the latter, we recall that a PL-
homeomorphism on $M$ is a homeomorphism $\gamma:M\to M$, such that there
exists some subdivision of $K$ such that $\sigma^{-1}\circ\gamma\circ\sigma$
is linear on each face. We denote the corresponding space of all orientation
preserving homeomorphisms by

$\Gamma_{\operatorname{PL}}(M)=\\{\gamma\in\operatorname{Hom}_{PL}(M,M):\,|D\gamma|>0\text{
a.e.}\\}.$ (56)

Note that any of these reparametrization groups act by composition from the
right on $\operatorname{Imm}_{\operatorname{Lip}}(M,\mathbb{R}^{3})$. In
addition to the action of these reparametrization groups, we also want to
identify surfaces that only differ by a translation. This leads us to consider
the following three quotient spaces:

$\displaystyle\mathcal{S}_{\operatorname{\infty}}:=\operatorname{Imm}_{\operatorname{Lip}}(M,\mathbb{R}^{3})/\Gamma_{\infty}(M)/\operatorname{trans}$
(57)
$\displaystyle\mathcal{S}_{\operatorname{Lip}}:=\operatorname{Imm}_{\operatorname{Lip}}(M,\mathbb{R}^{3})/\Gamma_{\operatorname{Lip}}(M)/\operatorname{trans}$
(58)
$\displaystyle\mathcal{S}_{\operatorname{PL}}:=\operatorname{Imm}_{\operatorname{Lip}}(M,\mathbb{R}^{3})/\Gamma_{\operatorname{PL}}(M)/\operatorname{trans},$
(59)

which will play a central role in the remainder of the article. Note that we
always consider immersions of Lipschitz class and only vary the regularity of
the group acting on this space.

###  3.2. The SRNF framework

The square root normal field (SRNF) map was first introduced by Jermyn et al.
in [21] for the space of smooth immersions. As we will see in the following
this mapping directly extends to all immersions of the Lipschitz class.

For any given $f\in\operatorname{Imm}_{\operatorname{Lip}}(M,\mathbb{R}^{3})$,
the orientation on $M$ allows us to consider the unit normal vector field
$n_{f}:M\to\mathbb{R}^{3}$, which is well-defined as a function in
$L^{\infty}(M,\mathbb{R}^{3})$. Furthermore, let $\\{v,w\\}$ be an orthonormal
basis of $T_{x}M$. Then for any
$f\in\operatorname{Imm}_{\operatorname{Lip}}(M,\mathbb{R}^{3})$ we can also
define the area multiplication factor at $x\in M$ via
$a_{f}(x)=|df_{x}(v)\times df_{x}(w)|$.

The SRNF map is then given by

$\displaystyle\Phi:\operatorname{Imm}_{\operatorname{Lip}}(M,\mathbb{R}^{3})/\operatorname{translations}$
$\displaystyle\to L^{2}(M,\mathbb{R}^{3})$ (60) $\displaystyle f$
$\displaystyle\mapsto q_{f},$ (61)

where

$q_{f}(x)=\sqrt{a_{f}(x)}\;n_{f}(x).$ (62)

We can now use the SRNF to define a pseudometric on
$\operatorname{Imm}_{\operatorname{Lip}}(M,\mathbb{R}^{3})/\operatorname{translations}$
by

$d_{\operatorname{Imm}}(f_{1},f_{2})=\|\Phi(f_{1})-\Phi(f_{2})\|_{L^{2}}.$
(63)

The function $d_{\operatorname{Imm}}$ is only a pseudo-metric due to the non-
injectivity of $\Phi$. Examples of this degeneracy have been discussed
extensively in the recent article [24].

Next we consider a right-action of $\Gamma_{\operatorname{Lip}}(M)$ on
$L^{2}(M,\mathbb{R}^{3})$ that is compatible with the mapping $\Phi$.
Therefore we let

$\displaystyle(g*\gamma)(x)=\sqrt{a_{\gamma}(x)}g(\gamma(x)),$ (64)

where $a_{\gamma}(x)$, the area multiplication factor of $\gamma$ at $x$, is
defined by $a_{\gamma}(x)=\sqrt{|D_{\gamma}(x)|}$. It is easy to check that

$\displaystyle\Phi(f)*\gamma=\Phi(f\circ\gamma)$ (65)

and that the action of $\Gamma_{\operatorname{Lip}}(M)$ on
$L^{2}(M,\mathbb{R}^{3})$ is by linear isometries, if we put the usual $L^{2}$
inner product on $L^{2}(M,\mathbb{R}^{3})$. Thus it follows that the SRNF
pseudometric on $\operatorname{Imm}_{\operatorname{Lip}}(M,\mathbb{R}^{3})$ is
invariant with respect to this action and thus it descends to a pseudometric
on the quotient space $\mathcal{S}_{\operatorname{Lip}}$, which is given by

$d_{\mathcal{S}_{\operatorname{Lip}}}([f_{1}],[f_{2}])=\inf\limits_{\gamma\in\Gamma_{\operatorname{Lip}}(M)}d(f_{1},f_{2}\circ\gamma),\qquad[f_{1}],[f_{2}]\in\mathcal{S}_{\operatorname{Lip}}(M)$
(66)

Since $\Gamma_{\infty}(M)$ and $\Gamma_{\operatorname{PL}}(M)$ are both
subsets of $\Gamma_{\operatorname{Lip}}$ the invariance properties continue to
hold for the actions of the smaller groups and we can also consider the
corresponding distance functions

$\displaystyle
d_{\mathcal{S}_{\infty}}([f_{1}],[f_{2}])=\inf\limits_{\gamma\in\Gamma_{\infty}(M)}d(f_{1},f_{2}\circ\gamma),\qquad[f_{1}],[f_{2}]\in\mathcal{S}_{\infty}(M)$
(67) $\displaystyle
d_{\mathcal{S}_{\operatorname{PL}}}([f_{1}],[f_{2}])=\inf\limits_{\gamma\in\Gamma_{\operatorname{PL}}(M)}d(f_{1},f_{2}\circ\gamma),\qquad[f_{1}],[f_{2}]\in\mathcal{S}_{\operatorname{PL}}(M)$
(68)

The remainder of this section will be devoted to showing that the SRNF
distance for each of these three group actions is equivelent. More precisely,
we aim to prove the following theorem:

######  Theorem 3.1.

Let $f_{1},f_{2}\in\operatorname{Imm}_{\operatorname{Lip}}(M,\mathbb{R}^{3})$.
Then

$d_{\mathcal{S}_{\operatorname{Lip}}}([f_{1}],[f_{2}])=d_{\mathcal{S}_{\infty}}([f_{1}],[f_{2}])=d_{\mathcal{S}_{\operatorname{PL}}}([f_{1}],[f_{2}]).$
(69)

######  Remark 3.2.

Note that this result implies in particular that for smooth immersions
$f_{1},f_{2}$ the SRNF metric as defined in this section is equal to the SRNF
metric considered in [21].

The main ingredient for proving this theorem will be the following lemma
concerning the continuity of the action of $\Gamma_{\operatorname{Lip}}(M)$ on
$L^{2}(M,\mathbb{R}^{3})$. Therefore we first note that
$\Gamma_{\operatorname{Lip}}(M)\subseteq W^{1,\infty}(M,M)\subseteq
W^{1,2}(M,M)$ and thus we can equip $\Gamma_{\operatorname{Lip}}(M)$ with the
$\|\cdot\|_{W^{1,2}}$ norm.

######  Lemma 3.3.

The map

$\displaystyle\left(L^{2}(M,\mathbb{R}^{3}),\|\cdot\|_{L^{2}}\right)\times\left(\Gamma_{\operatorname{Lip}}(M),\|\cdot\|_{W^{1,2}}\right)$
$\displaystyle\to\left(L^{2}(M,\mathbb{R}^{3}),\|\cdot\|_{L^{2}}\right)$ (70)
$\displaystyle(q,\gamma)$ $\displaystyle\mapsto q*\gamma$ (71)

is jointly continuous.

Since the $W^{1,\infty}$ topology dominates the $W^{1,2}$-topology, the lemma
would also hold when equipping $\Gamma_{\operatorname{Lip}}(M)$ with the
Lipschitz topology. The reason for equipping it instead with the
$W^{1,2}$-topology is that the subgroups of smooth and PL diffeomorphisms are
dense with respect to this topology, but not with respect to the
$W^{1,\infty}$-topology. The density of these groups will be a crucial
ingredient for our proof of Theorem 3.1.

###### Proof.

The proof of this result is inspired by [8, Proposition 7], where a related
result for one-dimensional domain space $M$ is shown.
Step 1 (Piecewise constant $q$) Let $q\in L^{2}(M,\mathbb{R}^{3})$ be
piecewise constant, i.e., there exists a disjoint family of sets $M_{j}$ such
that $q=\sum_{j=1}^{N}u_{j}\chi(M_{j})$. Given a sequence
$\gamma_{n}\to\gamma$ in
$\left(\Gamma_{\operatorname{Lip}}(M),\|\cdot\|_{W^{1,2}}\right)$, we need to
show $q*\gamma_{n}\to q*\gamma$ in $L^{2}(M,\mathbb{R}^{3})$. Since $M$ is a
2-manifold, the map

$\displaystyle W^{1,2}(M,M)\to L^{1}(M,M)\text{ defined via }\gamma\mapsto
a_{\gamma}$ (72)

is continuous by the Sobolev multiplication theorem. Thus,

$\int_{M}|a_{\gamma}-a_{\gamma_{n}}|dx\to 0.$ (73)

Further this implies

$\int_{M}|a_{\gamma}-a_{\gamma_{n}}|dx=\int_{M}\left(\sqrt{|a_{\gamma}-a_{\gamma_{n}}|}\right)^{2}dx\geq\int_{M}|\sqrt{a_{\gamma}}-\sqrt{a_{\gamma_{n}}}|^{2}dx$
(74)

and therefore we also have

$\int\limits_{M}\left|\sqrt{a_{\gamma}}-\sqrt{a_{\gamma_{n}}}\right|^{2}dx\to
0.$ (75)

Next we define the sets $M_{i,n}^{j}=\gamma_{n}(M_{j})\cap\gamma(M_{i})$.
Using this we can write the integral via a double sum as:

$\displaystyle\int_{M}|(q\circ\gamma)\sqrt{a_{\gamma}}-(q\circ\gamma_{n})\sqrt{a_{\gamma_{n}}}|^{2}dx=\sum_{i=1}^{N}\sum_{j=1}^{N}\int_{M_{i,n}^{j}}|u_{j}\sqrt{a_{\gamma}}-u_{i}\sqrt{a_{\gamma_{n}}}|^{2}dx$
(76)

Note that $\gamma_{n}\to\gamma$ in $W^{1,\infty}$ and thus for $i\neq j$,
$\mu(M_{i,n}^{j})\to 0$. Meanwhile

$\displaystyle\int_{M_{i,n}^{j}}|u_{j}\sqrt{a_{\gamma}}-u_{i}\sqrt{a_{\gamma_{n}}}|^{2}dx$
$\displaystyle\leq\int_{M_{i,n}^{j}}|u_{j}\sqrt{a_{\gamma}}-u_{i}\sqrt{a_{\gamma}}|^{2}dx+\|q\|^{2}_{\infty}\int_{M_{i,n}^{j}}|\sqrt{a_{\gamma}}-\sqrt{a_{\gamma_{n}}}|^{2}dx$
(77)
$\displaystyle=|u_{j}-u_{i}|^{2}\int_{M_{i,n}^{j}}|a_{\gamma}|dx+\|q\|^{2}_{\infty}\int_{M_{i,n}^{j}}|\sqrt{a_{\gamma}}-\sqrt{a_{\gamma_{n}}}|^{2}dx.$
(78)

Recall, $\gamma\in W^{1,\infty}(M,M)$ thus $|a_{\gamma}|$ is bounded. So

$|u_{j}-u_{i}|^{2}\int\limits_{M_{i,n}^{j}}|a_{\gamma}|dx\to 0$ (79)

as $\mu(M_{i,n}^{j})\to 0$. Additionally,

$\int\limits_{M}\left|\sqrt{a_{\gamma}}-\sqrt{a_{\gamma_{n}}}\right|^{2}dx\to
0\text{ therefore when }i\neq
j,\,\int_{M_{i,n}^{j}}|u_{j}\sqrt{a_{\gamma}}-u_{i}\sqrt{a_{\gamma_{n}}}|^{2}dx\to
0.$ (80)

In the case where $i=j$ we have

$\int_{M_{i,n}^{i}}|u_{i}\sqrt{a_{\gamma}}-u_{i}\sqrt{a_{\gamma_{n}}}|^{2}=\int_{M_{i,n}^{i}}|u_{i}|^{2}\cdot|\sqrt{a_{\gamma}}-\sqrt{a_{\gamma_{n}}}|^{2}dx\leq\|q\|^{2}_{\infty}\int_{M_{i,n}^{i}}|\sqrt{a_{\gamma}}-\sqrt{a_{\gamma_{n}}}|^{2}dx.$
(81)

Thus,

$\int_{M}|(q\circ\gamma)\sqrt{a_{\gamma}}-(q\circ\gamma_{n})\sqrt{a_{\gamma_{n}}}|^{2}dx\to
0.$ (82)

Step 2. (general $q$) Let now $q\in L^{2}(M,\mathbb{R}^{3})$,
$\gamma_{n}\to\gamma$ in $\Gamma_{\operatorname{Lip}}(M)$, and $\epsilon>0$.
Pick $q^{\prime}$ to be piecewise constant such that
$\|q-q^{\prime}\|_{L^{2}}<\epsilon/3$. Then using the fact that
$\Gamma_{\operatorname{Lip}}(M)$ acts by isometries and the triangle
inequality we have

$\displaystyle\|q*\gamma-q*\gamma_{n}\|_{L^{2}}$
$\displaystyle\leq\|q*\gamma-q^{\prime}*\gamma\|_{L^{2}}+\|q^{\prime}*\gamma-q^{\prime}*\gamma_{n}\|_{L^{2}}+\|q^{\prime}*\gamma_{n}-q*\gamma_{n}\|_{L^{2}}$
(83)
$\displaystyle\leq\|q-q^{\prime}\|_{L^{2}}+\|q^{\prime}*\gamma-q^{\prime}*\gamma_{n}\|_{L^{2}}+\|q^{\prime}-q\|_{L^{2}}$
(84)
$\displaystyle\leq\epsilon/3+\|q^{\prime}*\gamma-q^{\prime}*\gamma_{n}\|_{L^{2}}+\epsilon/3.$
(85)

From here we can conclude the result using Step 1.

Step 3. (joint continuity) Let now $q_{n}\to q$ in $L^{2}(M,\mathbb{R}^{3})$,
$\gamma_{n}\to\gamma$ in $\Gamma_{\operatorname{Lip}}(M)$, and $\epsilon>0$.
Using the fact that $\Gamma_{\operatorname{Lip}}(M)$ acts by isometries and
the triangle inequality we have

$\displaystyle\|q*\gamma-q_{n}*\gamma_{n}\|_{L^{2}}$
$\displaystyle\leq\|q*\gamma-q*\gamma_{n}\|_{L^{2}}+\|q*\gamma_{n}-q_{n}*\gamma_{n}\|_{L^{2}}$
(86) $\displaystyle=\|q*\gamma-q*\gamma_{n}\|_{L^{2}}+\|q-q_{n}\|_{L^{2}}$
(87)

From here we can conclude the result using Step 2. ∎

We are now able to prove Theorem 3.1:

###### Proof of Theorem 3.1.

Given $[f_{1}],[f_{2}]\in\mathcal{S}_{\operatorname{Lip}}$, choose
parametrized representations
$f_{1},f_{2}\in\operatorname{Imm}_{\operatorname{Lip}}(M,\mathbb{R}^{3})$.

Part 1.
$(\boldsymbol{d_{\mathcal{S}_{\operatorname{Lip}}}([f_{1}],[f_{2}])=d_{\mathcal{S}_{\infty}}([f_{1}],[f_{2}]))}$
Note that $\Gamma_{\infty}(M)$ is a subset of
$\Gamma_{\operatorname{Lip}}(M)$. Thus,

$\inf\limits_{\gamma\in\Gamma_{\operatorname{Lip}}(M)}\|\Phi(f_{1})-\Phi(f_{2}\circ\gamma)\|_{L^{2}}\leq\inf\limits_{\gamma\in\Gamma_{\infty}(M)}\|\Phi(f_{1})-\Phi(f_{2}\circ\gamma)\|_{L^{2}}.$
(89)

To show the opposite inequality, take an arbitrary
$g\in\Gamma_{\operatorname{Lip}}(M)$ and let $\epsilon>0$. By Lemma 3.2, there
exists some $\delta>0$ such that $\|g-\gamma\|_{W^{1,2}}<\delta$ implies
$\|\Phi(f_{2})*g-\Phi(f_{2})*\gamma\|<\epsilon$. Since $M$ is a 2-manifold,
$C^{\infty}(M,M)$ maps are dense in the Sobolev functions $W^{1,2}(M,M)$ [7].
Since $\Gamma_{\operatorname{Lip}}(M)\subseteq W^{1,2}(M,M)$, there exists
$\gamma\in\Gamma_{\infty}(M)$ such that $\|g-\gamma\|_{W^{1,2}}<\delta$. Thus,

$\displaystyle\|\Phi(f_{1})-\Phi(f_{2})*\gamma\|_{W^{1,2}}$
$\displaystyle\leq\|\Phi(f_{1})-\Phi(f_{2})*g\|_{W^{1,2}}+\|\Phi(f_{2})*g-\Phi(f_{2})*\gamma\|_{W^{1,2}}$
(90) $\displaystyle<\|\Phi(f_{1})-\Phi(f_{2})*g\|_{W^{1,2}}+\epsilon$ (91)

Therefore,

$\inf\limits_{\gamma\in\Gamma_{\operatorname{Lip}}(M)}\|\Phi(f_{1})-\Phi(f_{2}\circ\gamma)\|_{L^{2}}=\inf\limits_{\gamma\in\Gamma_{\infty}(M)}\|\Phi(f_{1})-\Phi(f_{2}\circ\gamma)\|_{L^{2}}.$
(92)

Part 2.
$(\boldsymbol{d_{\mathcal{S}_{\operatorname{Lip}}}([f_{1}],[f_{2}])=d_{\mathcal{S}_{PL}}([f_{1}],[f_{2}]))}$
This follows exactly as in Part 1, using that $\Gamma_{\operatorname{PL}}(M)$
is a dense subset of $\Gamma_{\operatorname{Lip}}(M)$. To prove the density we
use the fact that $\operatorname{Hom}_{PL}(M,M)$ maps are dense in the Sobolev
functions $W^{1,2}(M,M)$ as $\operatorname{dim}(M)=2$, cf. [57]. ∎

###  3.3. SRNF shape metric as an unbalanced optimal transport problem.

First we define a map from $L^{2}(M,\mathbb{R}^{3})$ to the the space of
positive finite Borel measures on $S^{2}$, and then show that computing the
shape distance between two surfaces is equivalent to computing the
Wasserstein-Fisher-Rao distance between the corresponding measures. For $q\in
L^{2}(M,\mathbb{R}^{3})$ and $U\subseteq S^{2}$ open, define

$q^{*}U=\\{x\in M|q(x)\neq 0\text{ and }q(x)/|q(x)|\in U\\}.$ (93)

Then we can define

$\displaystyle\psi:L^{2}(M,\mathbb{R}^{3})\to\mathcal{M}(S^{2})\text{ via
}\psi(q)=\mu_{q}$ (94) $\displaystyle\text{where for }U\subseteq
S^{2},\mu_{q}(U)=\int_{q^{*}U}q(x)\cdot q(x)dx.$ (95)

The main goal of this section is to show that the shape pseudo-distance in the
SRNF framework can be written as a pullback of the Wasserstein-Fisher-Rao
distance via the map $\psi$. We will show this result only for the dense
subset of PL surfaces. We expect that a careful analysis of the continuity of
the SRNF map would allow one to obtain this result also for general Lipschitz
immersions, which we plan to study in future work.

######  Theorem 3.4.

Given two PL surfaces $S_{1}$ and $S_{2}$ parameterized by $f_{1}$ and $f_{2}$
the SRNF shape distance can be computed as an unbalanced transport problem.
More precisely, we have

$d_{\mathcal{S}_{\operatorname{Lip}}}([f_{1}],[f_{2}])=\operatorname{WFR}(\psi\circ\phi(f_{1}),\psi\circ\phi(f_{2})).$
(96)

###### Proof.

Assume that $S_{1}$ and $S_{2}$ are triangulated compact oriented PL surfaces
in $\mathbb{R}^{3}$, and let $f:M\to S_{1}$ be a PL parametrization of
$S_{1}$. Let $\\{\sigma_{i}\\}_{i=1}^{m}$ denote the faces of $S_{1}$, and let
$\\{\tau_{j}\\}_{j=1}^{n}$ denote the faces of $S_{2}$. For each $i$, assume
that $\sigma_{i}$ has area $a_{i}$ and oriented unit normal vector $u_{i}$.
For each $j$, assume that $\tau_{j}$ has area $b_{j}$ and oriented unit normal
vector $v_{j}$. Let $\mathcal{A}$ be the set of all discrete semi-couplings
from $\psi(S_{1})$ to $\psi(S_{2})$.

######  Claim 1.

We may equivalently express the shape distance by

$d_{\mathcal{S}_{\operatorname{Lip}}}([f_{1}],[f_{2}])^{2}=\sum_{i=1}^{m}a_{i}+\sum_{j=1}^{n}b_{j}-2\sup_{g\in\operatorname{Hom}_{PL}^{+}(S_{1},S_{2})}\int_{M}\Phi(f_{1})\cdot\Phi(g\circ
f_{1}).$ (97)

Proof of Claim 1: Recall by Theorem 3.1,

$\displaystyle d_{\mathcal{S}_{\operatorname{Lip}}}([f_{1}],[f_{2}])^{2}$
$\displaystyle=\inf\limits_{\gamma\in\Gamma_{\operatorname{PL}}(M)}\|\Phi(f_{1})-\Phi(f_{2}\circ\gamma)\|^{2}_{L^{2}}$
(98)
$\displaystyle=\inf\limits_{\gamma\in\Gamma_{\operatorname{PL}}(M)}\int\limits_{M}|\Phi(f_{1})-\Phi(f_{2}\circ\gamma)|^{2}dx$
(99)
$\displaystyle=\sum_{i=1}^{m}a_{i}+\sum_{j=1}^{n}b_{j}-2\sup\limits_{\gamma\in\Gamma_{\operatorname{PL}}(M)}\int\limits_{M}|\Phi(f_{1})\cdot\Phi(f_{2}\circ\gamma)|dx$
(100)

Every $\gamma\in\Gamma_{\operatorname{PL}}(M)$, corresponds to
$g=f_{2}\circ\gamma\circ
f_{1}^{-1}\in\operatorname{Hom}_{PL}^{+}(S_{1},S_{2})$. Thus,

$d_{\mathcal{S}_{\operatorname{Lip}}}([f_{1}],[f_{2}])^{2}=\sum_{i=1}^{m}a_{i}+\sum_{j=1}^{n}b_{j}-2\sup_{g\in\operatorname{Hom}_{PL}^{+}(S_{1},S_{2})}\int_{M}\Phi(f_{1})\cdot\Phi(g\circ
f_{1}).$ (101)

This completes the proof of Claim 1. Now recall from Remark 2.2

$\operatorname{WFR}(\psi\circ\phi(f_{1}),\psi\circ\phi(f_{2}))^{2}=\sum_{i=1}^{m}a_{i}+\sum_{j=1}^{n}b_{j}-2\sup_{(A,B)\in\mathcal{A}}\sum_{i=1}^{m}\sum_{j=1}^{n}\sqrt{A_{ij}B_{ij}}(u_{i}\cdot
v_{j})$ (102)

Therefore, showing
$d_{\mathcal{S}_{\operatorname{Lip}}}([f_{1}],[f_{2}])=\operatorname{WFR}(\psi\circ\phi(f_{1}),\psi\circ\phi(f_{2}))$
is equivalent to showing

$\sup_{g\in\operatorname{Hom}_{PL}^{+}(S_{1},S_{2})}\int_{M}\Phi(f_{1})\cdot\Phi(g\circ
f_{1})dx=\sup_{(A,B)\in\mathcal{A}}\sum_{i=1}^{m}\sum_{j=1}^{n}\sqrt{A_{ij}B_{ij}}(u_{i}\cdot
v_{j}).$ (103)

######  Claim 2.

Assume that $(A,B)$ is a discrete semi-coupling from $\psi(S_{1})$ to
$\psi(S_{2})$. Then for all $\epsilon>0$ there is a PL homeomorphism
$g:S_{1}\to S_{2}$ such that

$\left|\int_{M}q_{f}\cdot q_{g\circ
f}-\sum_{i,j}\sqrt{A_{ij}B_{ij}}(u_{i}\cdot v_{j})\right|<\epsilon.$ (104)

Proof of Claim 2: Let $(A,B)$ be a discrete semi-coupling from $\psi(S_{1})$
to $\psi(S_{2})$ such that for each $1\leq i\leq m$ and $1\leq j\leq n$,
$A_{ij},B_{ij}>0$. We will first prove the claim for this restricted case and
extend it to all semi-couplings by continuity.

Let $a:M\to\mathbb{R}$ denote the area multiplication factor of $f$, and let
$u:M\to S^{2}$ denote the unit normal vector function corresponding to $f$.
First we choose a real number $r\in(0,1)$. For each $1\leq i\leq m$, subdivide
$\sigma_{i}$ into $n$ smaller 2-simplexes $\sigma_{ij}$ such that each
$\sigma_{ij}$ has area $A_{ij}$. Similarly, for each $1\leq j\leq n$,
subdivide $\tau_{j}$ into $m$ smaller 2-simplexes $\tau_{ij}$ such that each
$\tau_{ij}$ has area $B_{ij}$. For each $1\leq i\leq m$ and $1\leq j\leq n$,
choose a smaller 2-simplex $\tilde{\sigma}_{ij}$, whose closure is contained
in the interior of $\sigma_{ij}$, such that $\tilde{\sigma}_{ij}$ has area
equal to $rA_{ij}$. Similarly, for each $1\leq i\leq m$ and $1\leq j\leq n$,
choose a smaller 2-simplex $\tilde{\tau}_{ij}$, whose closure is contained in
the interior of $\tau_{ij}$, such that $\tilde{\tau}_{ij}$ has area equal to
$rB_{ij}$.

We now construct an orientation preserving PL homeomorphism $g_{r}:S_{1}\to
S_{2}$. First, for each $1\leq i\leq m$ and $1\leq j\leq n$, define
$g_{r}:\tilde{\sigma}_{ij}\to\tilde{\tau}_{ij}$ to be an arbitrary PL
orientation preserving homeomorphism with constant area multiplication factor.
Note that
$S_{1}-\left(\bigcup\limits_{i=1}^{m}\bigcup\limits_{j=1}^{n}\tilde{\sigma}_{ij}^{\mathrm{o}}\right)$
is homeomorphic to
$S_{2}-\left(\bigcup\limits_{i=1}^{m}\bigcup\limits_{j=1}^{n}\tilde{\tau}_{ij}^{\mathrm{o}}\right)$.
Hence, we can simply extend the homeomorphism $g_{r}$ which we already defined
on the $\tilde{\sigma}_{ij}$’s to a homeomorphism $S_{1}\to S_{2}$ in an
arbitrary manner. Denote the unit normal function coming from the
parametrization $g_{r}\circ f$ of $S_{2}$ by $v_{r}:M\to S^{2}$. Denote the
area multiplication factor of $g_{r}\circ f$ by $m_{r}:M\to\mathbb{R}$.

Write $M=M_{1}\cup M_{2}$, where
$M_{1}=f^{-1}\left(\bigcup\limits_{i=1}^{m}\bigcup\limits_{j=1}^{n}\tilde{\sigma}_{ij}\right)$
and $M_{2}=\overline{M-M_{1}}$. For each $1\leq i\leq m$ and $1\leq j\leq n$,
let $M_{ij}=f^{-1}(\tilde{\sigma}_{ij})$. Note that on each $M_{ij}$,
$m_{r}=\frac{aB_{ij}}{A_{ij}}$. Compute:

$\displaystyle\int_{M_{1}}q_{f}\cdot q_{g_{r}\circ f}$
$\displaystyle=\sum\limits_{i=1}^{m}\sum_{j=1}^{n}\int_{M_{ij}}\sqrt{a}\sqrt{m_{r}}(u_{i}\cdot
v_{j})=\sum\limits_{i=1}^{m}\sum_{j=1}^{n}\int_{M_{ij}}\sqrt{a}\sqrt{\frac{aB_{ij}}{A_{ij}}}(u_{i}\cdot
v_{j})$ (105)
$\displaystyle=\sum\limits_{i=1}^{m}\sum_{j=1}^{n}\int_{M_{ij}}a\sqrt{\frac{B_{ij}}{A_{ij}}}(u_{i}\cdot
v_{j})=\sum\limits_{i=1}^{m}\sum_{j=1}^{n}\operatorname{area}(\tilde{\sigma}_{ij})\sqrt{\frac{B_{ij}}{A_{ij}}}(u_{i}\cdot
v_{j})$ (106)
$\displaystyle=\sum\limits_{i=1}^{m}\sum_{j=1}^{n}\sqrt{\operatorname{area}(\tilde{\sigma}_{ij})}\sqrt{\frac{\operatorname{area}(\tilde{\sigma}_{ij})B_{ij}}{A_{ij}}}(u_{i}\cdot
v_{j})$ (107)
$\displaystyle=\sum\limits_{i=1}^{m}\sum_{j=1}^{n}\sqrt{\operatorname{area}(\tilde{\sigma}_{ij})}\sqrt{\operatorname{area}(\tilde{\tau}_{ij})}(u_{i}\cdot
v_{j})$ (108)
$\displaystyle=\sum\limits_{i=1}^{m}\sum_{j=1}^{n}\sqrt{rA_{ij}}\sqrt{rB_{ij}}(u_{i}\cdot
v_{j}).$ (109)

Meanwhile by the Schwarz inequality,

$\displaystyle\left|\int_{M_{2}}q_{f}\cdot q_{g_{r}\circ f}\right|$
$\displaystyle=\left|\int_{M_{2}}\sqrt{a}\sqrt{m_{r}}(u\cdot
v_{r})\right|\leq\int_{M_{2}}\sqrt{a}\sqrt{m_{r}}\left|u\cdot v_{r}\right|$
(110)
$\displaystyle\leq\int_{M_{2}}\sqrt{a}\sqrt{m_{r}}\leq\sqrt{\int_{M_{2}}a}\sqrt{\int_{M_{2}}m_{r}}$
(111)
$\displaystyle=\sqrt{\mbox{area}\left(S_{1}-\left(\bigcup\limits_{i=1}^{m}\bigcup\limits_{j=1}^{n}\tilde{\sigma}_{ij}^{\mathrm{o}}\right)\right)}\sqrt{\mbox{area}\left(S_{2}-\left(\bigcup\limits_{i=1}^{m}\bigcup\limits_{j=1}^{n}\tilde{\tau}_{ij}^{\mathrm{o}}\right)\right)}$
(112)
$\displaystyle=\sqrt{(1-r)\operatorname{area}(S_{1})}\sqrt{(1-r)\operatorname{area}(S_{2})}.$
(113)

So as we let $r\to 1$,

$\int_{M_{1}}q_{f}\cdot q_{g_{r}\circ
f}\to\sum\limits_{i=1}^{m}\sum_{j=1}^{n}\sqrt{A_{ij}B_{ij}}(u_{i}\cdot
v_{j})\operatorname{and}\int_{M_{2}}q_{f}\cdot q_{g_{r}\circ f}\to 0.$ (114)

Hence,

$\int_{M}q_{f}\cdot q_{g_{r}\circ
f}\to\sum\limits_{i=1}^{m}\sum_{j=1}^{n}\sqrt{A_{ij}B_{ij}}(u_{i}\cdot
v_{j}).$ (115)

Thus Claim 2 follows for the case in which for each $1\leq i\leq m$ and $1\leq
j\leq n$, $A_{ij}>0$ and $B_{ij}>0$. The general case then follows immediately
from the continuity of

$\sum\limits_{i=1}^{m}\sum_{j=1}^{n}\sqrt{A_{ij}B_{ij}}(u_{i}\cdot v_{j})$
(116)

as a function of $(A,B)$. This completes the proof of Claim 2. It follows that

$\sup_{g\in\operatorname{Hom}_{PL}^{+}(S_{1},S_{2})}\int_{M}q_{f}\cdot
q_{g\circ
f}\geq\sup_{\\{A_{ij}\\},\\{B_{ij}\\}\in\mathcal{A}}\sum_{i=1}^{m}\sum_{j=1}^{n}\sqrt{A_{ij}B_{ij}}(u_{i}\cdot
v_{j}).$ (117)

We are left to show the opposite inequality.

######  Claim 3.

Assume $g$ is a PL-homeomorphism from $S_{1}$ to $S_{2}$, then there exists a
discrete semi-coupling $(A,B)$ such that

$\int_{M}q_{f}\cdot q_{g\circ
f}\leq\sum_{i,j}\sqrt{A_{ij}}\sqrt{B_{ij}}(u_{i}\cdot v_{j}).$ (118)

Proof of Claim 3: Let $g:S_{1}\to S_{2}$ be an orientation preserving PL
homeomorphism. For $1\leq i\leq m$ and $1\leq j\leq n$, define
$\sigma_{ij}=g^{-1}(\tau_{j})\cap\sigma_{i}$ and define
$\tau_{ij}=g(\sigma_{ij})$. Now define two $(m+1)\times(n+1)$ matrices $A$ and
$B$ via:

  * •

For $1\leq i\leq m$ and $1\leq j\leq n$,
$A_{ij}=\operatorname{area}(\sigma_{ij})$ and
$B_{ij}=\operatorname{area}(\tau_{ij})$.

  * •

For $0\leq j\leq n$, $A_{j0}=0$.

  * •

For $0\leq i\leq m$, $B_{0i}=0$.

  * •

For $0\leq i\leq m$,

$A_{i0}=a_{i}-\sum\limits_{j=1}^{n}\operatorname{area}(\sigma_{ij}).$ (119)

  * •

For $0\leq j\leq n$,

$B_{0j}=b_{j}-\sum\limits_{i=1}^{m}\operatorname{area}(\tau_{ij}).$ (120)

The pair of matrices $(A,B)$ is a discrete semi-coupling from $\psi(S_{1})$ to
$\psi(S_{2})$ by construction. We say that $(A,B)$ is the semi-coupling
corresponding to the homeomorphism $g$. Let $1\leq i\leq m$ and $1\leq j\leq
n$ and let $M_{ij}=f^{-1}(\sigma_{ij})\subset M$. Denote the area
multiplication factor of $g$ on $\sigma_{ij}$ by $m_{ij}$. Then by the Schwarz
inequality,

$\displaystyle\int_{M_{ij}}q_{f}\cdot q_{g\circ f}$
$\displaystyle=\int_{M_{ij}}\sqrt{a(x)}u_{i}\cdot\sqrt{a(x)m_{ij}(x)}v_{j}\,dx=\int_{M_{ij}}\sqrt{a(x)}\sqrt{a(x)m_{ij}(x)}\,dx(u_{i}\cdot
v_{j})$ (121)
$\displaystyle\leq\sqrt{\int_{M_{ij}}a(x)\,dx\int_{M_{ij}}a(x)m_{ij}(x)\,dx}(u_{i}\cdot
v_{j})=\sqrt{A_{ij}}\sqrt{B_{ij}}(u_{i}\cdot v_{j})$ (122)

Summing over all $i$ and $j$ we obtain:

$\displaystyle\int_{M}q_{f}\cdot q_{g\circ f}$
$\displaystyle=\sum_{i,j}\int_{M_{ij}}q_{f}\cdot q_{g\circ
f}=\sum_{i,j}\int_{M_{ij}}a(x)m(x)(u_{i}\cdot
v_{j})\leq\sum_{i,j}\sqrt{A_{ij}}\sqrt{B_{ij}}(u_{i}\cdot v_{j}).$ (123)

This completes the proof of Claim 3. It follows that,

$\sup_{g\in\operatorname{Hom}_{PL}^{+}(S_{1},S_{2})}\int_{M}q_{f}\cdot
q_{g\circ
f}\leq\sup_{\\{A_{ij}\\},\\{B_{ij}\\}\in\mathcal{A}}\sum_{i=1}^{m}\sum_{j=1}^{n}\sqrt{A_{ij}B_{ij}}(u_{i}\cdot
v_{j}).$ (124)

and thus the theorem is proved. ∎

######  Remark 3.5.

In this section we have defined a mapping from the shape space of PL surfaces
in $\mathbb{R}^{3}$ to the space of finitely supported measures on $S^{2}$. We
have then shown that the SRNF (pseudo-) distance between two surfaces is equal
to the WFR distance between the two corresponding measures. It is shown in
[12, 14, 35] that the space of finitely supported measures on $S^{2}$ is a
geodesic length space. One might hope that geodesics in this space could
somehow be “lifted” to geodesics in the space of PL surfaces. The main problem
with this plan is that there is an infinite-dimensional space of surfaces
corresponding to each measure; see [24] for examples of arbitrarily high
dimensional spaces of surfaces corresponding to a single measure. Hence, there
is no unique way of lifting geodesics in the space of measures to the space of
surfaces. Because of this degeneracy, which is inherent to the SRNF, there is
no direct way to define a “geodesic” in the space of surfaces with respect to
the SRNF distance function. While it is true that existing methods (involving
gradient searches over the reparametrization group; see for example [22]) do
result in plausible-looking deformations from one surface to another, these
deformations are not geodesics in any strict mathematical sense. Thus in order
to define geodesics formally in the space of surfaces, one would have to
resolve the degeneracy of the SRNF distance function by, for example, adding
extra terms to the SRNF metric, see e.g. [21, 56].

The mapping we introduced does, however, restrict to a bijection between the
space of  closed convex PL surfaces in $\mathbb{R}^{3}$ and a subspace of the
finitely supported measures on $S^{2}$. Therefore it would be possible to
define geodesics in the space of closed convex surfaces as lifts of geodesics
in this subspace; this was done recently for the space of convex curves in
[11].

######  Remark 3.6.

Using the mapping that we have defined from the shape space of surfaces in
${\mathbb{R}}^{3}$ to measures on $S^{2}$, we can pull back any distance
function on the space of measures to obtain a pseudo-distance function on the
shape space of surfaces. While the main purpose of this paper is to show that
the SRNF pseudo-distance is obtained in this manner from the WFR distance on
the space of measures, it might be of interest to use other distance functions
on the space of measures to obtain interesting pseudo-distances on the space
of surfaces. A likely candidate for this would be the $\rho$-parametrized WFR
distance function defined in [14]. The corresponding pseudo-distance on
surfaces would have a natural interpretation as assigning different weights to
the direction of the normal vector as opposed to the shrinking or expansion of
area on the surfaces. As explained in Section 2, the algorithm introduced in
Section 2 would also provide a computation of this generalized version of the
SRNF pseudo-distance. Note that any pseudo-distance on surfaces obtained in
this way would have the same degeneracy as the SRNF pseudo-distance.

###  3.4. SRNF Computation Experiments

By Theorem 3.1, we can utilize Algorithm 1 to compute the exact SRNF pseudo-
distance directly between two simplicial meshes. The resulting method is
described in Algorithm 2.

Algorithm 2 Calculate $\operatorname{SRNF}$ shape distance

procedure SRNF_distance($S_{1},S_{2},\epsilon$)

$u_{i}\leftarrow$ the unit normal of the $i$th face of $S_{1}$

$a_{i}\leftarrow$ the area of the $i$th face of $S_{1}$

$v_{j}\leftarrow$ the unit normal of the $i$th face of $S_{2}$

$b_{j}\leftarrow$ the area of the $i$th face of $S_{2}$

return WFR_distance($a,b,u,v,\epsilon$)

To quantify the performance of Algorithm 2 we compare it to the method
introduced in [32]. In their implementation surfaces are assumed to be
represented by a spherical parametrization and the diffeomorphism group is
discretized using spherical harmonics. This in turn allows one to formulate
the SRNF distance computation as a constrained minimization problem over the
coefficients of the reparametrization in the chosen spherical harmonics basis.
In the following, we will refer to this method as the parametrization-based
method.

Most data that one encounters in real applications is, however, not given in
such a parametrized form but rather as a simplicial complex. Thus one first
has to solve the parametrization problem [51], which is of comparable
difficulty to the geodesic boundary problem itself. In our experiments, we
used 24 shapes from the TOSCA dataset (simplicial meshes), for which we also
had access to spherical parametrizations, which have been calculated using the
method of [44] as implemented in [29]. We present the result of this
comparison in Figure 2 which consists of three subplots: Figure 2 (a), which
highlights the distances computed with both methods for 552 pairs of PL
surfaces. The precise distances produced by our method (Orange) are
consistently lower than the parametrization-based distances produced by the
method of [32] (Blue). The mean relative error of the parametrization-based
method compared to our method is $174.582\%$ with a standard deviation of
$90.077\%$. Note that this error consists of both an approximation error of
the spherical parameterization of the simplicial complex and an approximation
error of the optimal reparametrization. In Figure 2 (b), we plot the
correlation between these two methods of computing the distances which have a
Pearson correlation coefficient of 0.793. Note that this is comparable to the
correlation of elastic distances between functions on the line that are either
computed using dynamic programming or computed using an exact algorithm, cf.
[20, 55, 33].

(a) SRNF shape distances produced by Algorithm 2 (Orange) and by the method of
[32] (Blue).

(b) Correlation plot of distances from the method of [32] and Algorithm 2.

Figure 2. Comparison of Algorithm 2 and the method of [32] for 552 pairs of
surfaces from the TOSCA dataset.The shape distances used in this experiment
are sampled from pairs of shapes from the TOSCA data set as displayed in
Figure 5.

The main drawback of our method is that it does not produce an optimal
reparameterization that aligns the two surfaces, i.e., we do not obtain point
correspondences between the two meshes. We can, however, still interpret the
information in the optimal semi-coupling to visualize a “fuzzy” correspondence
between the surfaces. This is visualized in Figure 3: given two PL surfaces
$S_{1}$ and $S_{2}$, we color each face of $S_{2}$ according to its unit
normal vector. Further let $(A,B)$ be the optimal semi-coupling between
$\psi(S_{1})$ and $\psi(S_{2})$. We then color each face of $S_{1}$ with
normal vector $u_{i}$ according to the color of the face of $S_{2}$ with
normal vector $v_{k}$ where $k=\operatorname*{arg\,max}_{j}A_{ij}$, i.e., each
face on $S_{1}$ is colored by the same color as the face where most of its
mass is transported. Examples of such correspondences are presented in Figure
3.

|  | $S_{1}$ |  |  | $S_{2}$
---|---|---|---|---|---
|  |  |  |  |
Figure 3. Fuzzy Point Correspondences: Each face of the source surface
($S_{2}$) on the right is colored according to its normal vector. The faces of
the target surfaces ($S_{1}$) on the left are colored according to the optimal
discrete semi-coupling from $\psi(S_{1})$ to $\psi(S_{2})$.

The phenomenon of distinct closed surfaces that were indistinguishable by the
SRNF shape pseudo-distance was first studied in [24]. As a result of Theorem
3.1, we obtain a full characterization of this phenomenon: Note that the
Wasserstein-Fisher-Rao distance is a true distance and so the SRNF distance
between two surfaces is zero if and only if they are mapped to the same
measure by $\psi$. A useful observation from the study of convex polyhedra,
due to Minkowski, is that every measure on $S^{2}$ satisfying the closure
condition corresponds to a unique (up to translation) closed, convex
polyhedron in $\mathbb{R}^{3}$, see e.g. [48]. Therefore, each closed PL
surface has SRNF distance zero from this unique closed convex polyhedron. In
Figure 4, we give examples of PL surfaces and the convex polyhedron
reconstructed from the corresponding measure. This reconstruction is performed
using the Python package polyhedrec [50] available at
https://github.com/gsellaroli/polyhedrec.

h |  |  |
---|---|---|---
|  |  |
$\wr$ |  $\wr$ |  $\wr$ |  $\wr$
|  |  |
Figure 4. Four pairs of distinct shapes indistinguishable by the SRNF metric:
Each shape on the bottom row is the unique closed, convex polyhedron
reconstructed from the area measure associated with the corresponding shape on
the top row. Thus the SRNF shape distance between each of these pairs is zero.

Despite these known drawbacks, the SRNF pseudo-metric has been demonstrated to
be successful in applications surrounding the classification of surfaces. We
demonstrate this by considering a toy example of shapes from the TOSCA
dataset, that includes 4 cats, 7 dogs, 17 gorillas, 10 horses, and 9
lionesses. We then compute the SRNF pseudo-distance matrix using our Algorithm
2 and visualize the results using 3d multi-dimensional scaling in Figure 5.
One can see that the SRNF distance, despite its degeneracy, produces
meaningful clusters.

Figure 5. Multi-Dimensional Scaling Plot: A subset of the TOSCA dataset is
plotted according to 3d multi-dimensional scaling of the SRNF distance matrix
produced by Algorithm 2.

##  4\. Conclusion

In this article we propose a novel method to _precisely_ compute the SRNF
shape distance between PL surfaces. This method follows from three results
each of which has interesting implications in its own right:

First, we propose a novel method for computing the Wasserstein-Fisher-Rao
distance in unbalanced optimal transport. While fast estimations of this
metric can be achieved by including an entropy regularization term, we propose
a new method that solves exactly for the WFR distance.

Second, we extend the SRNF framework to surfaces parameterized by Lipschitz
immersions. This class of surfaces notably includes both smooth surfaces and
PL surfaces. We then show that this extended framework is consistent with the
original SRNF distance which was formulated in the smooth category.

Finally, we establish an equivalence between the SRNF and the Wasserstein-
Fisher-Rao distance on the space of Borel measures on $S^{2}$. In addition to
establishing a new method to compute the SRNF shape distance, this result
offers insight into theoretical problems that exist surrounding the SRNF shape
distance. For instance, this result gives us tools to analyze the phenomenon
of distinct closed surfaces that are indistinguishable by the SRNF shape
pseudo-distance as highlighted in Figure 4.

Open Questions and Future Work: This project provides several open questions
and opportunities for future work. The first set of open questions concerns
Algorithm 1: in all of Section 2.4, we consider specifically the Wasserstein-
Fisher-Rao distance between measures supported on $S^{2}$. An obvious open
problem is to investigate whether the results of Section 2.4 and the
associated algorithm can be generalized to a larger class of domains and
general unbalanced OMT problems. Extending this algorithm to more general
domains should follow easily from alternate characterizations of the distance
developed in [13, 14, 35]. Generalizing the algorithm to other unbalanced OMT
problems may be more challenging and require restricting to unbalanced OMT
problems that only optimize cost functions satisfying certain properties.
There are also possibilities for future work in making our algorithm more
efficient, by implementing sparse data types and/or absorption thresholds as
entries of the discrete semi-coupling approach zero. Further, as future work,
we would like to update our implementation to work with more general weight
matrices so that our algorithm can be effectively used for other applications
of unbalanced optimal transport.

The main result of this paper gives us new tools for studying the SRNF shape
distance. Currently, we have only proven this result for PL shapes and we
leave the extension of Theorem 3.4 to all surfaces in
$\mathcal{S}_{\operatorname{Lip}}$ as an open problem. We expect that this
result could follow from Theorem 3.4 using the density of PL surfaces in
$\mathcal{S}_{\operatorname{Lip}}$ with respect to the SRNF pseudo-metric and
carefully studying the continuity of all involved operations. Further, there
are opportunities for future work in characterizing the relationships between
shapes that are indistinguishable via the SRNF shape distance. Understanding
this relationship may serve to help in developing meaningful SRNF based
interpolations between shapes based on optimal discrete semi-couplings between
the associated measures.

Another interesting subject for future research is to gain a better
understanding of the set of all shapes of closed surfaces corresponding to a
given measure on $S^{2}$. It seems clear that this set is, in general,
infinite dimensional. See [24] for a construction of arbitrarily high
dimensional sets of this form. Consequently there are several natural
questions that arise: Does the geometry of this set depend on the given
measure on $S^{2}$? Can the shape space of surfaces be thought of as a fiber
bundle over the space of measures on $S^{2}$?

## References

  * [1] AD Alexandrov.  Zur theorie der gemischten volumina von konvexen körpern i.  Mat. Sbornik NS, 1:227–251, 1938.
  * [2] Martin Bauer, Martins Bruveris, and Peter W Michor.  Overview of the geometries of shape spaces and diffeomorphism groups.  Journal of Mathematical Imaging and Vision, 50(1):60–97, 2014.
  * [3] Martin Bauer, Nicolas Charon, Philipp Harms, and Hsi-Wei Hsieh.  A numerical framework for elastic surface matching, comparison, and interpolation.  International Journal of Computer Vision, pages 1–20, 2021.
  * [4] Martin Bauer, Philipp Harms, and Peter W Michor.  Sobolev metrics on shape space of surfaces.  Journal of Geometric Mechanics, 3(4), 2011.
  * [5] Jean-David Benamou and Yann Brenier.  A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem.  Numerische Mathematik, 84(3):375–393, 2000.
  * [6] Dimitri P Bertsekas.  Nonlinear programming.  Journal of the Operational Research Society, 48(3):334–334, 1997\.
  * [7] Fabrice Bethuel and Xiaomin Zheng.  Density of smooth functions between two manifolds in Sobolev spaces.  Journal of functional analysis, 80(1):60–75, 1988.
  * [8] Martins Bruveris.  Optimal reparametrizations in the square root velocity framework.  SIAM Journal on Mathematical Analysis, 48(6):4335–4354, 2016.
  * [9] Martins Bruveris, Peter W Michor, and David Mumford.  Geodesic completeness for Sobolev metrics on the space of immersed plane curves.  In Forum of Mathematics, Sigma, volume 2. Cambridge University Press, 2014.
  * [10] Rainer Burkard, Mauro Dell’Amico, and Silvano Martello.  Assignment problems: revised reprint.  SIAM, 2012.
  * [11] Nicolas Charon and Thomas Pierron.  On length measures of planar closed curves and the comparison of convex shapes.  Annals of Global Analysis and Geometry, 60(4):863–901, 2021.
  * [12] Lenaic Chizat, Gabriel Peyré, Bernhard Schmitzer, and François-Xavier Vialard.  An interpolating distance between optimal transport and Fisher–Rao metrics.  Foundations of Computational Mathematics, 18(1):1–44, 2018.
  * [13] Lenaic Chizat, Gabriel Peyré, Bernhard Schmitzer, and François-Xavier Vialard.  Scaling algorithms for unbalanced optimal transport problems.  Mathematics of Computation, 87(314):2563–2609, 2018.
  * [14] Lenaic Chizat, Gabriel Peyré, Bernhard Schmitzer, and François-Xavier Vialard.  Unbalanced optimal transport: Dynamic and Kantorovich formulations.  Journal of Functional Analysis, 274(11):3090–3123, 2018.
  * [15] Marco Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  In C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 26. Curran Associates, Inc., 2013.
  * [16] Ian L Dryden and Kanti V Mardia.  Statistical shape analysis: with applications in R, volume 995.  John Wiley & Sons, 2016.
  * [17] Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z. Alaya, Aurélie Boisbunon, Stanislas Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, Léo Gautheron, Nathalie T.H. Gayraud, Hicham Janati, Alain Rakotomamonjy, Ievgen Redko, Antoine Rolet, Antony Schutz, Vivien Seguy, Danica J. Sutherland, Romain Tavenard, Alexander Tong, and Titouan Vayer.  Pot: Python optimal transport.  Journal of Machine Learning Research, 22(78):1–8, 2021.
  * [18] Thomas Gallouët, Roberta Ghezzi, and François-Xavier Vialard.  Regularity theory and geometry of unbalanced optimal transport.  arXiv preprint arXiv:2112.11056, 2021.
  * [19] Steven Haker, Lei Zhu, Allen Tannenbaum, and Sigurd Angenent.  Optimal mass transport for registration and warping.  International Journal of computer vision, 60(3):225–240, 2004.
  * [20] Emmanuel Hartman, Yashil Sukurdeep, Nicolas Charon, Eric Klassen, and Martin Bauer.  Supervised deep learning of elastic SRV distances on the shape space of curves.  Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops. 2021., 2021.
  * [21] Ian H Jermyn, Sebastian Kurtek, Eric Klassen, and Anuj Srivastava.  Elastic shape matching of parameterized surfaces using square root normal fields.  In European conference on computer vision, pages 804–817. Springer, 2012.
  * [22] Ian H Jermyn, Sebastian Kurtek, Hamid Laga, and Anuj Srivastava.  Elastic shape analysis of three-dimensional objects.  Synthesis Lectures on Computer Vision, 12(1):1–185, 2017.
  * [23] Shantanu H Joshi, Qian Xie, Sebastian Kurtek, Anuj Srivastava, and Hamid Laga.  Surface shape morphometry for hippocampal modeling in alzheimer’s disease.  In 2016 International Conference on Digital Image Computing: Techniques and Applications (DICTA), pages 1–8. IEEE, 2016.
  * [24] Eric Klassen and Peter W Michor.  Closed surfaces with different shapes that are indistinguishable by the SRNF.  Archivum Mathematicum, 56(2):107–114, 2020.
  * [25] Stanislav Kondratyev, Léonard Monsaingeon, Dmitry Vorotnikov, et al.  A new optimal transport distance on the space of finite radon measures.  Advances in Differential Equations, 21(11/12):1117–1164, 2016.
  * [26] Sebastian Kurtek, Eric Klassen, Zhaohua Ding, and Anuj Srivastava.  A novel Riemannian framework for shape analysis of 3D objects.  In 2010 IEEE computer society conference on computer vision and pattern recognition, pages 1625–1632. IEEE, 2010.
  * [27] Sebastian Kurtek, Eric Klassen, John C Gore, Zhaohua Ding, and Anuj Srivastava.  Elastic geodesic paths in shape space of parameterized surfaces.  IEEE transactions on pattern analysis and machine intelligence, 34(9):1717–1730, 2011.
  * [28] Sebastian Kurtek, Chafik Samir, and Lemlih Ouchchane.  Statistical shape model for simulation of realistic endometrial tissue.  In ICPRAM, pages 421–428, 2014.
  * [29] Sebastian Kurtek, Anuj Srivastava, Eric Klassen, and Hamid Laga.  Landmark-guided elastic shape analysis of spherically-parameterized surfaces.  In Computer graphics forum, volume 32, pages 429–438. Wiley Online Library, 2013.
  * [30] Hamid Laga, Yulan Guo, Hedi Tabia, Robert B Fisher, and Mohammed Bennamoun.  3D Shape analysis: fundamentals, theory, and applications.  John Wiley & Sons, 2018.
  * [31] Hamid Laga, Marcel Padilla, Ian H Jermyn, Sebastian Kurtek, Mohammed Bennamoun, and Anuj Srivastava.  4d atlas: Statistical analysis of the spatiotemporal variability in longitudinal 3D shape data.  arXiv preprint arXiv:2101.09403, 2021.
  * [32] Hamid Laga, Qian Xie, Ian H Jermyn, and Anuj Srivastava.  Numerical inversion of SRNF maps for elastic shape analysis of genus-zero surfaces.  IEEE transactions on pattern analysis and machine intelligence, 39(12):2451–2464, 2017.
  * [33] Sayani Lahiri, Daniel Robinson, and Eric Klassen.  Precise matching of PL curves in $\mathbb{R}^{N}$ in the square root velocity framework.  Geometry, Imaging and Computing, 2(3):133–186, 2015.
  * [34] Matthias Liero, Alexander Mielke, and Giuseppe Savaré.  Optimal transport in competition with reaction: The Hellinger–Kantorovich distance and geodesic curves.  SIAM Journal on Mathematical Analysis, 48(4):2869–2911, 2016.
  * [35] Matthias Liero, Alexander Mielke, and Giuseppe Savaré.  Optimal entropy-transport problems and a new Hellinger–Kantorovich distance between positive measures.  Inventiones mathematicae, 211(3):969–1117, 2018.
  * [36] Jan Maas, Martin Rumpf, Carola Schönlieb, and Stefan Simon.  A generalized model for optimal transport of images including dissipation and density modulation.  ESAIM: Mathematical Modelling and Numerical Analysis, 49(6):1745–1769, 2015.
  * [37] James Matuk, Shariq Mohammed, Sebastian Kurtek, and Karthik Bharath.  Biomedical applications of geometric functional data analysis.  In Handbook of Variational Methods for Nonlinear Geometric Data, pages 675–701. Springer, 2020.
  * [38] Quentin Mérigot.  A multiscale approach to optimal transport.  In Computer Graphics Forum, volume 30, pages 1583–1592. Wiley Online Library, 2011.
  * [39] Peter W Michor and David Mumford.  Vanishing geodesic distance on spaces of submanifolds and diffeomorphisms.  Documenta Mathematica, 10:217–245, 2005.
  * [40] Gaspard Monge.  Mémoire sur la théorie des déblais et des remblais.  Histoire de l’Académie Royale des Sciences de Paris, 1781.
  * [41] Xavier Pennec.  Intrinsic statistics on riemannian manifolds: Basic tools for geometric measurements.  Journal of Mathematical Imaging and Vision, 25(1):127–154, 2006\.
  * [42] Gabriel Peyré, Marco Cuturi, et al.  Computational optimal transport: With applications to data science.  Foundations and Trends® in Machine Learning, 11(5-6):355–607, 2019.
  * [43] Benedetto Piccoli and Francesco Rossi.  Generalized Wasserstein distance and its application to transport equations with source.  Archive for Rational Mechanics and Analysis, 211(1):335–358, 2014\.
  * [44] Emil Praun and Hugues Hoppe.  Spherical parametrization and remeshing.  ACM Transactions on Graphics (TOG), 22(3):340–349, 2003.
  * [45] Yossi Rubner, Carlo Tomasi, and Leonidas J Guibas.  The earth mover’s distance as a metric for image retrieval.  International journal of computer vision, 40(2):99–121, 2000.
  * [46] Martin Rumpf and Max Wardetzky.  Geometry processing from an elastic perspective.  GAMM-Mitteilungen, 37(2):184–216, 2014.
  * [47] Martin Rumpf and Benedikt Wirth.  Variational methods in shape analysis.  Handbook of Mathematical Methods in Imaging, 2:1819–1858, 2015\.
  * [48] Rolf Schneider.  Convex surfaces, curvature and surface area measures.  In Handbook of convex geometry, pages 273–299. Elsevier, 1993.
  * [49] Rolf Schneider.  Convex bodies: the Brunn–Minkowski theory.  Number 151. Cambridge university press, 2014.
  * [50] Giuseppe Sellaroli.  An algorithm to reconstruct convex polyhedra from their face normals and areas.  arXiv preprint arXiv:1712.00825, 2017.
  * [51] Alla Sheffer, Emil Praun, Kenneth Rose, et al.  Mesh parameterization methods and their applications.  Foundations and Trends in Computer Graphics and Vision, 2(2):105–171, 2007.
  * [52] Richard Sinkhorn.  A relationship between arbitrary positive matrices and doubly stochastic matrices.  The annals of mathematical statistics, 35(2):876–879, 1964.
  * [53] Justin Solomon.  Transportation Techniques for Geometric Data Processing.  Stanford University, 2015.
  * [54] Justin Solomon, Raif Rustamov, Leonidas Guibas, and Adrian Butscher.  Wasserstein propagation for semi-supervised learning.  In International Conference on Machine Learning, pages 306–314. PMLR, 2014.
  * [55] Anuj Srivastava and Eric P Klassen.  Functional and shape data analysis, volume 1.  Springer, 2016.
  * [56] Zhe Su, Martin Bauer, Stephen C Preston, Hamid Laga, and Eric Klassen.  Shape analysis of surfaces using general elastic metrics.  Journal of Mathematical Imaging and Vision, 62:1087–1106, 2020\.
  * [57] Jean Van Schaftingen.  Approximation in Sobolev spaces by piecewise affine interpolation.  Journal of Mathematical Analysis and Applications, 420(1):40–47, 2014.
  * [58] Cédric Villani.  Topics in optimal transportation.  Number 58. American Mathematical Soc., 2003.
  * [59] Cédric Villani.  Optimal transport: old and new, volume 338.  Springer Science & Business Media, 2008.
  * [60] John Henry C Whitehead.  On C1-complexes.  Annals of Mathematics, pages 809–824, 1940.
  * [61] Laurent Younes.  Computable elastic distances between shapes.  SIAM Journal on Applied Mathematics, 58(2):565–586, 1998.
  * [62] Laurent Younes.  Shapes and diffeomorphisms, volume 171.  Springer, 2010.