# Deep learning of diffeomorphisms for optimal reparametrizations of shapes

Elena Celledoni111Norwegian University of Science and Technology, Trondheim,
Norway<EMAIL_ADDRESS>Helge Glöckner222University of Paderborn,
Germany<EMAIL_ADDRESS>
Jørgen Riseth333Simula Research Laboratory, University of Oslo, Norway
<EMAIL_ADDRESS>and Alexander Schmeding444Nord universitetet i Levanger,
Norway<EMAIL_ADDRESS>
###### Abstract

In shape analysis, one of the fundamental problems is to align curves or
surfaces before computing a (geodesic) distance between these shapes. To find
the optimal reparametrization realizing this alignment is a computationally
demanding task which leads to an optimization problem on the diffeomorphism
group. In this paper, we construct approximations of orientation-preserving
diffeomorphisms by composition of elementary diffeomorphisms to solve the
approximation problem. We propose a practical algorithm implemented in PyTorch
which is applicable both to unparametrized curves and surfaces. We derive
universal approximation results and obtain bounds for the Lipschitz constant
of the obtained compositions of diffeomorphisms.

MSC2020: 65K10 (Primary); 58D05, 46T10 (Secondary)
Keywords: optimal reparametrization, shape analysis, deep learning,
diffeomorphism group

##  1 Introduction

The geometric shape is of major importance in object recognition. The field of
shape analysis gives a formal definition of the shape of an object, and
provides tools to compare these shapes. One way to represent objects
mathematically is to consider them as parametric curves or surfaces defined on
some compact domain $\Omega\subseteq\operatorname{\mathbb{R}}^{n}$. Since the
same shape may be outlined by different parametrizations, shape distances are
typically defined taking the infimum over all possible parametrizations, i.e.
over the group of orientation preserving diffeomorphisms,
$\text{Diff}^{+}(\Omega)$.

In this article, we propose an algorithm for solving optimization problems on
the group of orientation-preserving diffeomorphisms
$\mathrm{Diff}^{+}(\Omega)$ by composition of elementary diffeomorphisms.
Using concepts from deep learning, we restate the reparametrization problem
such that it may be solved by training a residual neural network. This
approach admits an effective, unified framework for optimal reparametrization
of both curves and surfaces. We motivate the approach with results on
infinite-dimensional Lie groups of diffeomorphisms by first describing global
charts for certain diffeomorphism groups, thus implying universal
approximation properties. Furthermore, we establish a priori estimates for the
composition of multiple diffeomorphisms. These estimates are of independent
interest beyond the applications considered in the present paper. Finally, we
provide an implementation of the algorithm based on the PyTorch-framework [32]
and test it for several reparametrization problems. The source code is
available at https://zenodo.org/record/6873960.

###  1.1 Related work

Invertible neural networks are build by the composition of invertible (layer)
maps. One primary example in this class are residual networks [20], where each
layer is of the form $x\mapsto x+f(x).$ Here
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ depends on parameters, e.g. for so
called dense layers $f(x):=\sigma(Ax+b)$, where $A$ is a $d\times d$ matrix,
$b\in\mathbb{R}^{d}$ and $\sigma$ is a scalar nonlinearity applied
componentwise. Residual networks may further be related to approximations of
diffeomorphisms, by interpreting them as flows of ordinary differential
equations [17, 12, 6]. Assuming the parameters to be time-dependent functions
$A(t)$ and $b(t)$, each layer may be interpreted as a forward Euler
discretization of the non-autonomous differential equation

$\dot{x}(t)=f(x,t),\quad x(0)=x,\qquad f(x,t):=\sigma(A(t)x+b(t)),$ (1.1)
$x\mapsto\psi_{L},\quad\psi_{k}=\psi_{k-1}+hf(\psi_{k-1},t_{k}),\quad
k=1,\dots,L,\quad\psi_{0}=x,$

where $h$ is the step-size often taken to be equal to $1$, and with
$t_{k}=kh$. For $f$ Lipschitz continuous in $x$, the Euler method is known to
converge on bounded time intervals to the solution of the differential
equation as $h\rightarrow 0$. A sufficient condition for the invertibility of
the layer map is that $hf$ is Lipschitz continuous, with Lipschitz constant
$\mathrm{Lip}(hf)<1$ (also called $1$-Lipschitz). For activation functions
$\sigma$ with bounded derivative most layer types, e.g. dense layers or
convolutional layers, can be made $1$-Lipschitz by appropriate scaling [9],
see for example [5] for a concrete strategy for choosing such a scaling.

The necessity of invertible layers in neural network is not restricted to the
ODE-interpretation of residual neural networks. For example, normalizing flows
are a class of machine learning models that may be used to artificially
generate data [38]. Given $m$ independent observations of $n$ random variables
one seeks an estimate of their underlying probability density $\rho(x)$. The
approximation of the probability density is achieved by mapping $x$ to a new
set of variables $y=\varphi^{-1}(x)$ with known (for example normal)
distribution $\mu(y)$. Then

$\rho(x)=J_{y}(x)\mu(y(x)),$

where $J_{y}$ is the Jacobian of the map $\varphi^{-1}$. The map
$\varphi^{-1}$ is built by means of a gradient ascent flow increasing the log-
likelihood of the sample data with respect to its density.

While the evolution and growth of the field of deep learning throughout the
last decade have primarily been driven by empirical results, there are
multiple approximation theorems which provide a theoretical basis for the
effectiveness of neural networks. For example, consider the space of one-layer
neural networks of the form

$\psi:[0,1]^{n}\to\operatorname{\mathbb{R}},\qquad\psi(x):=\sum_{j=1}^{N}\alpha_{j}\sigma(w_{j}^{T}x+b_{j}),$
(1.2)

with parameters $\alpha_{j},b_{j}\in\operatorname{\mathbb{R}}$,
$\mathbf{w}_{j}\in\operatorname{\mathbb{R}}^{n}$, and some activation function
$\sigma:\operatorname{\mathbb{R}}\to\operatorname{\mathbb{R}}$. In 1989,
Cybenko [11] proved that such networks are dense in the space of continuous
functions with respect to the supremum norm as $N\to\infty$, provided that
$\sigma$ is continuous and satisfies

$\sigma(t)\rightarrow\left\\{\begin{array}[]{ccc}1&&t\rightarrow+\infty,\\\
0&&t\rightarrow-\infty.\end{array}\right.$

In [1], the authors provide a similar approximation result for diffeomorphisms
isotopic to the identity, defined on a compact manifold $\Omega$. They show
that such diffeomorphisms may be expressed as a finite composition of
exponential maps of vector fields that have been scaled by smooth functions,
as long as the vector fields belong to a bracket-generating family of
functions. In [2], this result is extended to more general maps on $\Omega$.

###  1.2 Structure of the Article

The remainder of this article is structured as follows: Section 2 provides
background on shape analysis and the definition of the proposed deep-learning
approach. Section 3 gives universal approximation results and the theoretical
basis for the deep learning approach. Section 4 describes a priori estimates
of the $C^{k}$-distance between the identity diffeomorphism and the network
$\varphi$, given in terms of elementary diffeomorphisms. Section 5 provides
details on how the proposed algorithm is implemented. Finally, in Section 6 we
test the algorithms for a couple of examples involving both parametric curves,
surfaces and images.

##  2 Background and definition of the proposed numerical method

###  2.1 Shape space and shape metrics

##### Orientation-Preserving Diffeomorphisms

Consider the space of immersions from $\Omega=[0,1]$ into euclidean space
$\operatorname{\mathbb{R}}^{d}$,

$\mathcal{P}:=\text{Imm}(\Omega,\operatorname{\mathbb{R}}^{d})=\left\\{c\in
C^{\infty}(\Omega,\operatorname{\mathbb{R}}^{d})\;:\;c^{\prime}(t)\neq
0,\;\forall t\in\Omega\right\\}.$ (2.1)

To identify two curves representing the same shape, consider the set of
orientation-preserving diffeomorphisms, which consists of monotonically
increasing functions from $\Omega$ onto itself,

$\text{Diff}^{+}(\Omega)=\left\\{\varphi\in
C^{\infty}(\Omega,\Omega)\;:\;\varphi(0)=0,\;\varphi(1)=1,\;\;\varphi^{\prime}(x)>0,\;\forall
x\in\Omega\right\\}.$ (2.2)

The group $\text{Diff}^{+}(\Omega)$ is an infinite-dimensional Lie group. We
shall provide more details on this group in Section 3 below. At the moment, it
is important to notice the canonical right group action of
$\text{Diff}^{+}(\Omega)$ on $\mathcal{P}$ by composition,

$\mathcal{P}\times\text{Diff}^{+}(\Omega)\to\mathcal{P},\qquad(c,\varphi)\mapsto
c\circ\varphi.$

We shall refer to this as reparametrization of $c$ by $\varphi$. With respect
to the reparametrization action we can now describe shape spaces.

##### Shape Space

Elements of shape space are interpreted as unparametrized curves. The shape
space is the set of $\text{Diff}^{+}(\Omega)$-orbits

$\mathcal{S}=\mathcal{P}\;/\;\text{Diff}^{+}(\Omega).$

Our goal is to obtain geometrically sound deformations between elements in
$\mathcal{S}$ and compute distances thereof, by means of a metric
$d_{\mathcal{S}}$. Such a metric is typically constructed by defining a
reparametrization-invariant metric $d_{\mathcal{P}}$ on $\mathcal{P}$. Due to
reparametrization invariance (cf. [8]), $d_{\mathcal{P}}$ descends to a metric

$d_{\mathcal{S}}([c_{1}],[c_{2}])=\inf_{\varphi\in\text{Diff}^{+}(\Omega)}d_{\mathcal{P}}(c_{1},c_{2}\circ\varphi).$
(2.3)

Hence, to compute the distance between two shapes, we need to solve an
optimization problem. The algorithm presented in this article aims to solve
this problem.

##### The SRVT and Q-transform for Curves

To define a suitable distance function $d_{\mathcal{P}}$, one can exploit
certain transformations which transform familiar (Riemannian) metrics to
obtain reparametrization-invariant (Riemannian) metrics. One such
transformation is the square-root velocity transform (SRVT) for curves,
introduced in [36]. It is given by

$\mathcal{R}\colon\mathcal{P}\to
C^{\infty}(\Omega,\operatorname{\mathbb{R}}^{d})\setminus\\{0\\},\quad
c\mapsto\frac{c^{\prime}}{\sqrt{|c^{\prime}|}}.$

It is easy to see that
$\mathcal{R}(c\circ\varphi)(t)=\sqrt{\varphi^{\prime}(t)}\;(\mathcal{R}(c)\circ\varphi)(t)$.
Hence we can pull back the standard $L^{2}$-inner product on function spaces
to obtain a reparametrization-invariant distance:

$d_{\mathcal{P}}(c_{1},c_{2}):=\|\mathcal{R}(c_{1})-\mathcal{R}(c_{2})\|_{L^{2}(\Omega,\operatorname{\mathbb{R}}^{d})}=\left(\int_{\Omega}|\mathcal{R}(c_{1})(t)-\mathcal{R}(c_{2})(t)|^{2}\;dt\right)^{1/2}.$

Another transformation with similar properties as the SRVT, which allows us to
construct a different metric, is the so-called _Q-transform_ [26, 3]:

$Q\colon\mathcal{P}\to C^{\infty}(\Omega,\operatorname{\mathbb{R}}^{d}),\qquad
c(\cdot)\mapsto\sqrt{|c^{\prime}(\cdot)|}c(\cdot).$ (2.4)

The SRVT was generalized to curves on Lie groups and homogeneous spaces [7].

##### Surfaces embedded in 3D space

To extend the framework for computing shape distances to parametric surfaces
embedded in $\operatorname{\mathbb{R}}^{3}$, we will assume for simplicity
that the surfaces are defined on the unit square $\Omega=[0,1]^{2}$. The main
differences between shapes of curves and of surfaces pertain the group of
diffeomorphisms. Further, we need analogues to the SRVT and Q-transform.

Denote by $f_{x},f_{y}$ the partial derivatives of $f$ and by

$a_{f}(\mathbf{x})=\left|f_{x}(\mathbf{x})\times f_{y}(\mathbf{x})\right|$

the area scaling factor of the surface. Now the pre-shape space for surfaces
is

$\mathcal{P}=\left\\{f:C^{\infty}(\Omega,\operatorname{\mathbb{R}})\;:\;a_{f}(x)>0\;\forall
x\in\Omega\right\\}.$ (2.5)

The group of orientation-preserving diffeomorphisms consists of smooth
invertible maps from $\Omega$ onto itself such that the Jacobian determinant
$J_{\varphi}$ is positive,

$\text{Diff}^{+}(\Omega)=\left\\{\varphi\in
C^{\infty}(\Omega,\Omega)\;:\;J_{\varphi}(\mathbf{x})>0\;\forall\mathbf{x}\in\Omega,\;\phi(s)\in\partial\Omega\;\forall
s\in\partial\Omega\right\\}.$

Its tangent space $T_{\varphi}\text{Diff}^{+}(\Omega)$ at $\varphi$ consists
of smooth boundary-respecting vector fields on $\Omega$, see [27].

##### The SRNF and Q-transform of Surfaces

The square root normal field (SRNF), defined in [21], is a generalization of
the SRVT to surfaces. It is defined by

$\mathcal{N}(f)(\mathbf{x})=\sqrt{a_{f}(\mathbf{x})}\mathbf{n}_{f}(\mathbf{x)}$

where $\mathbf{n}_{f}(\mathbf{x})$ is the normal vector field on $f$. Applied
to a reparametrized surface, the transform satisfies the relation

$\mathcal{N}(f\circ\varphi)=\sqrt{J_{\varphi}}(\mathcal{N}(f)\circ\varphi).$

Following the same derivation as was given for curves, we may now define a
reparametrization-invariant pre-shape distance by

$d_{\mathcal{P}}(f_{1},f_{2})=\|\mathcal{N}(f_{1})-\mathcal{N}(f_{2})\|_{L^{2}(\Omega,\operatorname{\mathbb{R}}^{3})}$

since the factor $\sqrt{J_{\varphi}}$ cancels the Jacobian determinant
introduced by the change of variables in the integral. While the SRNF shares
some properties of the SRVT, for example the translation invariance, it is not
invertible. In [23], the authors present several examples of surfaces that
represent different shapes, but have the same image under the SRNF. However,
the SRNF framework has shown some merit in practice, see [37] and the
references therein.

Alternatively, we can use the $Q$-transform for surfaces, which was introduced
in [24]. The $Q$-transform for surfaces is defined by

$Q(f)(\mathbf{x})=\sqrt{a_{f}(\mathbf{x})}f(\mathbf{x})$

and satisfies the same property of reparametrization-invariance with respect
to the $L^{2}$-norm.

###  2.2 The optimization problem on the diffeomorphism group

The metric on shape space is induced by suitable metrics on the pre-shape
space. In the approaches sketched above, this metric is the geodesic distance
of a Riemannian metric arising as pullback of the $L^{2}$-metric via one of
the transformations mentioned. To compute the metric on shape space one needs
to find the elements in the $\operatorname{Diff}^{+}(\Omega)$-orbits which are
closest to each other. However, due to equivariance of the transformations,
this can equivalently be formulated as the following optimization problem for
the $L^{2}$-distance (exploiting that
$(C^{\infty}(\Omega,\operatorname{\mathbb{R}}^{d}),L^{2})$ is a pre-Hilbert
space and the set in which we compute the distance is an open (but non-convex)
subset of this space). Assume that $f_{1},f_{2}\in\mathcal{P}$, and we want to
compute $d_{\mathcal{S}}([f_{1}],[f_{2}])$. Let $\mathcal{T}$ be one of the
transforms discussed above and $q_{i}:=\mathcal{T}(f_{i}),i=1,2$. Then we need
to solve the following optimization problem:

$\displaystyle\varphi^{*}={\arg\inf}_{\varphi\in\mathrm{Diff}^{+}(\Omega)}E(\varphi),\qquad
E(\varphi):=\|q_{1}-\sqrt{\varphi^{\prime}}\,(q_{2}\circ\varphi)\|^{2}_{L^{2}}.$
(2.6)

Note that the optimization problem has been formulated as the squared
$L^{2}$-distance. Solving this optimization problem, we can compute the metric
distance on shape space if both $q_{1}$ and $q_{2}$ are contained in a convex
subset of the image of the transformation. While this is not always the case
(for example, it does not hold for a pair of curves which get mapped by the
SRVT to elements $q_{1}$ and $-q_{1}$), the assumption always holds locally
for mappings which are not too far apart. In any case, solving efficiently
this optimization problem is crucial in shape analysis and the main purpose of
the present paper is to provide tools from machine learning for its solution.

###  2.3 A Riemannian gradient descent method

A gradient descent approach is obtained by representing the gradient of the
functional $E$ by means of an othonormal basis of
$T_{\mathrm{id}}\mathrm{Diff}^{+}(\Omega)$. Suppose $f_{j}\in
T_{\mathrm{id}}\mathrm{Diff}^{+}(\Omega)$, $j=1,2\dots$ is an orthonormal
basis of $T_{\mathrm{id}}\mathrm{Diff}^{+}(\Omega)$ with respect to an inner
product $\langle\cdot,\cdot\rangle_{\mathrm{id}}$ on
$T_{\mathrm{id}}\mathrm{Diff}^{+}(\Omega)$. Consider the left multiplication
$L_{\varphi}:\psi\mapsto\varphi\circ\psi$ on $\mathrm{Diff}^{+}(\Omega)$. We
denote with
$\mathrm{d}L_{\varphi}:T_{\mathrm{id}}\mathrm{Diff}^{+}(\Omega)\rightarrow
T_{\varphi}\mathrm{Diff}^{+}(\Omega)$ the corresponding derivative mapping at
$\mathrm{id}$. The gradient of $E$ with respect to
$\langle\cdot,\cdot\rangle_{\mathrm{id}}$ is defined as the unique vector
field $\mathrm{grad}\,E$ on $\mathrm{Diff}^{+}(\Omega)$ satisfying

$\left.\mathrm{d}E\right|_{\varphi}(w_{\varphi})=\langle\mathrm{d}L_{\varphi}^{-1}(\mathrm{grad}\,E({\varphi})),\mathrm{d}L_{\varphi}^{-1}w_{\varphi}\rangle_{\mathrm{id}},\qquad\forall
w_{\varphi}\in T_{\varphi}\mathrm{Diff}^{+}(\Omega).$

We can represent $\mathrm{grad}\,E({\varphi})$, at least formally, using the
basis of $T_{\mathrm{id}}\mathrm{Diff}^{+}(\Omega)$

$\mathrm{d}L_{\varphi}^{-1}(\mathrm{grad}\,E({\varphi}))=\sum_{j=1}^{\infty}\lambda_{j}f_{j},$
(2.7)

with
$\lambda_{j}=\left.\mathrm{d}E\right|_{\varphi}(\mathrm{d}L_{\varphi}(f_{j}))$,
and $\delta
E_{\mathrm{id}}:=\mathrm{d}L_{\varphi}^{-1}(\mathrm{grad}\,E({\varphi}))$.

In practice, the proposed gradient descent algorithm on
$\mathrm{Diff}^{+}(\Omega)$ is obtained by truncating (2.7) and projecting the
gradient on the subspace of $T_{\mathrm{id}}\mathrm{Diff}^{+}(\Omega)$ spanned
by a finite number $M$ of basis elements. So

$\mathrm{d}L_{\varphi}^{-1}\widehat{\mathrm{grad}\,E}(\varphi)=\sum_{j=1}^{M}\lambda_{j}f_{j},$

and we obtain the following update rule for the gradient descent iteration:

$\varphi^{(n+1)}=\varphi^{(n)}\circ\left(\mathrm{id}-\eta\,\mathrm{d}L_{\varphi^{(n)}}^{-1}\widehat{\mathrm{grad}\,E}\big{(}\varphi^{(n)}\big{)}\right),\quad
n=0,1,\dots,$

with $\varphi_{0}=\mathrm{id}$, where $\eta$ is a scalar parameter optimized
to speed-up convergence and to guarantee the invertibility of
$\varphi^{(n+1)}$. Written in terms of the basis, the gradient descent
iteration becomes

$\varphi^{(n+1)}=\varphi^{(n)}\circ\left(\mathrm{id}-\eta\,\sum_{j=1}^{M}\lambda_{j}^{n}f_{j}\right),\quad
n=0,1,\dots,$ (2.8)

and $\\{\lambda_{j}^{n}\\}_{j=1}^{M}$ are coefficients determined at each
iteration $n$. In [34], this update rule was shown to be equivalent to an
algorithm proposed in [25] and to a similar algorithm created for curves in
[36].

###  2.4 Main method: deep learning of diffeomorphisms

It can be observed that after $L$ iterations of the gradient descent algorithm
(2.8), the diffeomorphism $\varphi^{L}$ is the composition of $L$ elementary
diffeomorphisms each defined as the sum of $M$ basis functions. This
observation may be exploited to create a deep learning approach to the
reparametrization problem.

Let $q_{1}$ and $q_{2}$ be defined as in eq. 2.6. Instead of using the whole
group $\operatorname{Diff}^{+}(\Omega)$, we fix some $L,M\in\mathbb{N}$ and
restrict ourselves to functions of the form

$\varphi=\varphi_{L}\circ\cdots\circ\varphi_{1},\quad\mathrm{where}\quad\varphi_{\ell}:=\mathrm{id}+\sum_{n=1}^{M}w_{n}^{\ell}f_{n},\qquad\ell=1,\dots,L,\qquad$
(2.9)

for a set of basis function $\\{f_{j}\\}^{M}_{j=1}\subset
T_{\operatorname{id}}\text{Diff}^{+}(\Omega)$. Here
$\\{w_{j}^{\ell}\\}_{j=1}^{M}$ are coefficients to be found by the
optimization algorithm. This form is similar to the structure of a deep
residual network [20], and allows us to restate the optimization problem as
the training of the network. To ensure
$\varphi\in\operatorname{Diff}^{+}(\Omega)$, we restrict the set of weights to
some set $\mathcal{W}\subset\operatorname{\mathbb{R}}^{ML}$. Further details
on this restriction are given in section 5.1.

The increased performance of the deep learning approach (2.9) can be seen in
figure 1 which compares the proposed method to the gradient descent algorithm.
This figure indicates that the deep learning approach converges towards a
significantly better solution than gradient descent.

Figure 1: Comparison between gradient descent and multi-layer approaches for
reparametrization of curves. The GD approach used 6 basis functions, whereas
the multi-layer approach used 6 basis functions and 6 layers. (Top): A curve
$c_{2}$ (left) and the same curve with a different parametrization
$c_{1}=c_{2}\circ\varphi$ (right). The goal is to reconstruct $\varphi$.
(Bottom left): The diffeomorphisms found by the different algorithms, compared
to the target diffeomorphism $\varphi$. The multi-layer approach is visually
indistinguishable from the target. (Bottom right): Convergence plots for the
two algorithms. The error is given relative to the initial error.

A preliminary implementation of this approach and proof of concept was
presented in [34].

##  3 Theoretical justification

This section provides the main theoretical motivation for the proposed deep
learning approach for optimization on $\mathrm{Diff}^{+}(\Omega)$. We shall
first recall relevant (infinite-dimensional) Lie group structures on
diffeomorphism groups. Then we show that finite compositions of
diffeomorphisms of the type $\varphi_{\ell}=\operatorname{id}+f_{\ell}$,
$\ell=1,\dots,L$ with $f_{\ell}$ a $1$-Lipschitz vector field can be used to
describe the whole group of diffeomorphisms of a cube $\Omega=[0,1]^{d}$
(Section 3.2). If $\Omega$ is, more generally, a compact convex subset of
${\mathbb{R}}^{d}$ with dense interior (like a disk
$\Omega=\\{(x,y)\in\operatorname{\mathbb{R}}^{2}:x^{2}+y^{2}\leq r^{2}\\}$ in
$\operatorname{\mathbb{R}}^{2}$), analogous diffeomorphisms generate the group
of diffeomorphisms fixing the boundary of $\Omega$ (see Section 3.1). We are
mostly interested in the case of a cube $[0,1]^{d}$; in fact, only the unit
interval $[0,1]$ and square $[0,1]^{2}$ are used in the numerical
calculations. As explained in Section 2.4, eventually our methodology is
implemented first restricting $f_{\ell}$ to a finite-dimensional subspace of
the Lie algebra of vector fields and then constraining $f_{\ell}$ to be
$1$-Lipschitz, see Section 6 for details. This introduces further
approximations which we discuss briefly but do not analyze in detail in the
present work.

###  3.1 Diffeomorphisms fixing the boundary

Lie group structures on diffeomorphism groups of compact manifolds have been
studied in many works (see, e.g., [27, 18, 28, 31]), including the case of
manifolds with boundary or corners. In this section, we consider a compact
convex subset $\Omega\subseteq{\mathbb{R}}^{d}$ with dense interior (having in
mind the main example $\Omega=[0,1]^{d}$). A map on $\Omega$ is smooth if it
extends to a smooth map on an open subset in ${\mathbb{R}}^{d}$; we can
therefore speak about the group $\operatorname{Diff}(\Omega)$ of all
$C^{\infty}$-diffeomorphisms of $\Omega$. Following [15], we consider the
subgroup

$\operatorname{Diff}_{\partial}(\Omega):=\\{\phi\in\operatorname{Diff}(\Omega):\phi(k)=k,\forall
k\in\partial\Omega\\}$

of all smooth diffeomorphisms fixing the topological boundary $\partial\Omega$
pointwise. It has been shown in [15] that if we endow
$\operatorname{Diff}_{\partial}(\Omega)$ with the compact-open
$C^{\infty}$-topology, the group becomes a manifold modeled on the space

$C^{\infty}_{\partial}(\Omega,\operatorname{\mathbb{R}}^{d})=\\{f\in
C^{\infty}(\Omega,\operatorname{\mathbb{R}}^{d}):f(k)=0,\forall
k\in\partial\Omega\\}.$

This space is an infinite-dimensional locally convex space (with respect to
the compact-open $C^{\infty}$-topology) and its elements can be identified
with vector fields which vanish on the boundary. This structure turns
$\operatorname{Diff}_{\partial}(\Omega)$ into an infinite-dimensional Lie
group. A canonical identification yields
$T_{\operatorname{id}}\operatorname{Diff}_{\partial}(\Omega)=C^{\infty}_{\partial}(\Omega,\operatorname{\mathbb{R}}^{d})$.
Further, $\operatorname{Diff}_{\partial}(\Omega)$ admits a global chart

$\displaystyle\kappa\colon\operatorname{Diff}_{\partial}(\Omega)\rightarrow
C^{\infty}_{\partial}(\Omega,\operatorname{\mathbb{R}}^{d}),\quad\phi\mapsto\phi-\operatorname{id}_{\Omega},$
(3.1)

which arises as a restriction of the map $\kappa\colon
C^{\infty}(\Omega,\operatorname{\mathbb{R}}^{d})\rightarrow
C^{\infty}(\Omega,\operatorname{\mathbb{R}}^{d}),f\mapsto
f-\operatorname{id}_{\Omega}$. For vector fields
$f\in\kappa(\operatorname{Diff}_{\partial}(\Omega))$, the inverse
$\kappa^{-1}(f)=f+\operatorname{id}_{\Omega}$ is in
$\operatorname{Diff}_{\partial}(\Omega)$ and we can obtain all elements in
$\operatorname{Diff}_{\partial}(\Omega)$ in this way.

By the preceding, every diffeomorphism fixing the boundary can be expressed as
a vector field (vanishing on the boundary) plus the identity. Thus, one only
needs to approximate vector fields and can ignore the non-linearity of the
infinite-dimensional manifold. However, there are two practical problems for
us to solve:

  1. 1.

Describe the image of the global chart $\kappa$.

  2. 2.

Establish approximation results by sampling from a finite-dimensional subspace
of the Lie algebra.

As the global chart is highly dependent on the geometry of $\Omega$ there is
no easy comprehensive way to describe its image. However, it is easy to
provide sufficient conditions for a vector field to lie inside the image of
$\kappa$.

######  Example 3.0.1.

Let $\Omega=[a,b]$ be a compact interval. Then a vector field $f$ in
$C^{\infty}_{\partial}([a,b],\operatorname{\mathbb{R}})$ is in the image of
$\kappa$ if $f^{\prime}(x)>-1$ for all $x\in[a,b]$ since then
$\operatorname{id}_{[a,b]}+f$ will be monotonically increasing.

The situation is more complicated for higher-dimensional sets. We do not have
any global information concerning the set $\Omega$ then. However, there is a
convenient sufficient condition ensuring that a map is in $\Omega$, which we
describe now.

Fix any norm on $\operatorname{\mathbb{R}}^{d}$ to calculate Lipschitz
constants for functions from $\Omega\subseteq\operatorname{\mathbb{R}}^{d}$ to
$\operatorname{\mathbb{R}}^{d}$, operator norms, and open balls
$B_{r}(x):=\\{y\in\operatorname{\mathbb{R}}^{d}\colon\|y-x\|<r\\}$ for
$x\in\operatorname{\mathbb{R}}^{d}$ and $r>0$. Setting

$\|f\|_{\infty,\operatorname{op}}:=\max_{x\in\Omega}\|f^{\prime}(x)\|_{\operatorname{op}}$
(3.2)

for $f\in C^{\infty}_{\partial}(\Omega,\operatorname{\mathbb{R}}^{d})$, we
obtain a continuous seminorm on
$C^{\infty}_{\partial}(\Omega,\operatorname{\mathbb{R}}^{d})$. Hence

$\mathcal{U}_{1}:=\\{f\in
C^{\infty}_{\partial}(\Omega,\operatorname{\mathbb{R}}^{d})\colon\operatorname{Lip}(f)<1\\}=\\{f\in
C^{\infty}_{\partial}(\Omega,\operatorname{\mathbb{R}}^{d})\colon\|f\|_{\infty,\operatorname{op}}<1\\}$

is an open $0$-neighborhood in
$C^{\infty}_{\partial}(\Omega,\operatorname{\mathbb{R}}^{d})$.

######  Lemma 3.1.

The map $\operatorname{id}_{\Omega}+f$ is an element of
$\operatorname{Diff}_{\partial}(\Omega)$ for all $f\in\mathcal{U}_{1}$. Thus
$\mathcal{U}_{1}\subseteq\kappa(\operatorname{Diff}_{\partial}(\Omega))$ and
it makes sense to consider $\kappa^{-1}(f)=\operatorname{id}_{\Omega}+f$.

###### Proof.

To keep the notation short, we write
$\kappa^{-1}(f):=f+\operatorname{id}_{\Omega}$, though strictly speaking it
makes only sense to use the inverse of the chart once the proof has been
concluded. Given $f\in\mathcal{U}_{1}$, the map
$\kappa^{-1}(f)\colon\Omega\to\operatorname{\mathbb{R}}^{d}$ is injective,
since $\kappa^{-1}(f)(x)=\kappa^{-1}(f)(y)$ with $x\not=y$ in $\Omega$ would
entail that

$\|x-y\|=\|f(y)-f(x)\|\leq\operatorname{Lip}(f)\|x-y\|<\|x-y\|,$

contradiction. Since $\kappa^{-1}(f)(x)=x$ for each $x\in\partial\Omega$, we
have

$\kappa^{-1}(f)(\partial\Omega)=\partial\Omega.$ (3.3)

Using the injectivity of $\kappa^{-1}(f)$, we obtain for the interior
$\Omega^{\circ}$ of $\Omega$

$\kappa^{-1}(f)(\Omega^{\circ})\cap\partial\Omega=\emptyset.$ (3.4)

Since $\|f^{\prime}(x)\|_{\operatorname{op}}<1$, we have
$\kappa^{-1}(f)^{\prime}(x)=\operatorname{id}_{\operatorname{\mathbb{R}}^{d}}+f^{\prime}(x)\in\operatorname{GL}(\operatorname{\mathbb{R}}^{d})$
for each $x\in\Omega^{\circ}$, whence $\kappa^{-1}(f)|_{\Omega^{\circ}}$ is a
local $C^{\infty}$-diffeomorphism at each $x\in\Omega^{\circ}$ (by the Inverse
Function Theorem) and hence an open map.

Next, we pick an $x_{0}$ in the interior $\Omega^{\circ}$ of $\Omega$. Then
the distance

$r:=\min\\{\|x_{0}-y\|\colon y\in\partial\Omega\\}$

is positive. Note that

$B_{r}(x_{0}):=\\{y\in\operatorname{\mathbb{R}}^{d}\colon\|x_{0}-y\|<r\\}\subseteq\Omega^{\circ}.$
(3.5)

In fact, since $B_{r}(x_{0})\cap\partial\Omega=\emptyset$, the disjoint open
sets $B_{r}(x_{0})\cap\Omega^{\circ}\not=\emptyset$ and
$Q:=B_{r}(x_{0})\cap(\operatorname{\mathbb{R}}^{d}\setminus\Omega)$ cover
$B_{r}(x_{0})$. Since $B_{r}(x_{0})$ is connected, $Q=\emptyset$ follows.
There exists $y_{0}\in\partial\Omega$ such that $\|x_{0}-y_{0}\|=r$. Since
$f(y_{0})=0$, we have

$\|f(x_{0})\|=\|f(x_{0})-f(y_{0})\|\leq\operatorname{Lip}(f)\|x_{0}-y_{0}\|<r.$

Hence

$\|\kappa^{-1}(f)(x_{0})-x_{0}\|=\|f(x_{0})\|<r,$

whence $\kappa^{-1}(f)(x_{0})\in\Omega^{\circ}$, by (3.5). Since
$x_{0}\in\Omega^{\circ}$ was arbitrary,

$\kappa^{-1}(f)(\Omega^{\circ})\subseteq\Omega^{\circ}$

is established. Now $\kappa^{-1}(f)(\Omega)$ is compact and hence closed in
$\operatorname{\mathbb{R}}^{d}$. Since
$\kappa^{-1}(f)(\partial\Omega)=\partial\Omega$, we deduce that

$\kappa^{-1}(f)(\Omega^{\circ})=\kappa^{-1}(f)(\Omega^{\circ})\cap\Omega^{\circ}=\kappa^{-1}(f)(\Omega)\cap\Omega^{\circ}$

is closed in $\Omega^{\circ}$. Since $\Omega^{\circ}$ is connected and its
non-empty subset $\kappa^{-1}(f)(\Omega^{\circ})$ is both open and closed in
$\Omega^{\circ}$, we must have
$\kappa^{-1}(f)(\Omega^{\circ})=\Omega^{\circ}$. Being a bijective local
$C^{\infty}$-diffeomorphism, the map
$\kappa^{-1}(f)|_{\Omega^{\circ}}\colon\Omega^{\circ}\to\Omega^{\circ}$ is a
$C^{\infty}$-diffeomorphism. Since $f|_{\partial\Omega}=0$, we have
$\kappa^{-1}(f)|_{\partial\Omega}=\operatorname{id}_{\partial\Omega}$, whence
$\kappa^{-1}(f)(\Omega)=\Omega$. It only remains to show that
$\kappa^{-1}(f)\colon\Omega\to\Omega$ is a $C^{\infty}$-diffeomorphism. By the
preceding, $\kappa^{-1}(f)\colon\Omega\to\Omega$ is a continuous bijective
self-map of the compact topological space $\Omega$ an hence a homeomorphism.
Note that
$\kappa^{-1}(f)^{\prime}(x)=\operatorname{id}_{\operatorname{\mathbb{R}}^{d}}+f^{\prime}(x)\in\operatorname{GL}(\operatorname{\mathbb{R}}^{d})$
for all $x\in\Omega$. If $k\in\mathbb{N}_{0}$ and $\kappa^{-1}(f)$ is a
$C^{k}$-diffeomorphism, then

$(\kappa^{-1}(f)^{-1})^{\prime}(x)=\big{(}\kappa^{-1}(f)^{\prime}(\kappa^{-1}(f)^{-1}(x))\big{)}^{-1}$

is a $\operatorname{GL}(\operatorname{\mathbb{R}}^{d})$-valued
$C^{k}$-function of $x\in\Omega$, entailing that $\kappa^{-1}(f)^{-1}$ is
$C^{k+1}$. Thus $\kappa^{-1}(f)^{-1}$ is smooth and thus
$\kappa^{-1}(f)\in\operatorname{Diff}_{\partial}(\Omega)$. ∎

######  Remark 3.2.

The same reasoning shows that, for each $\ell\in\mathbb{N}$, the set $\\{f\in
C^{\ell}_{\partial}(\Omega,\operatorname{\mathbb{R}}^{d})\colon\operatorname{Lip}(f)<1\\}$
is an open $0$-neighborhood in
$C^{\ell}_{\partial}(\Omega,\operatorname{\mathbb{R}}^{d})$ and
$\operatorname{id}_{\Omega}+f\in\operatorname{Diff}_{\partial}^{C^{\ell}}(\Omega)$
for all $f$ in this $0$-neighborhood.

The main disadvantage of the group $\operatorname{Diff}_{\partial}(\Omega)$ in
the applications we wish to investigate is that its elements fix the boundary
pointwise. Depending on the application, this can be quite unnatural. Hence,
we also consider the full group of $C^{\infty}$-diffeomorphisms for the
special case of a cube $[0,1]^{d}\subseteq\operatorname{\mathbb{R}}^{d}$.

###  3.2 Diffeomorphisms of a cube

In this section, we exclusively consider the cube $\Omega:=[0,1]^{d}$. Then
the group $\operatorname{Diff}(\Omega)$ of smooth diffeomorphisms is an
infinite-dimensional Lie group (cf. [27]). We shall follow the complementary
approach of [13], where diffeomorphism groups of convex polytopes are
discussed. Our goal is to generalize Lemma 3.1 to the diffeomorphism group of
a cube.

For the cube $\Omega$, the diffeomorphism group $\operatorname{Diff}(\Omega)$
is an infinite-dimensional Lie group modeled on the Fréchet space
$C^{\infty}_{\text{str}}(\Omega,\operatorname{\mathbb{R}}^{d})\subseteq
C^{\infty}(\Omega,\operatorname{\mathbb{R}}^{d})$ of smooth vector fields

$f=(f_{1},\ldots,f_{n})\colon\Omega\to{\mathbb{R}}^{d}$

which are _tangent to the boundary_ in the sense that
$f_{j}(x_{1},\ldots,x_{d})=0$ whenever $x_{j}=0$ or $x_{j}=1$ (such vector
fields are also called boundary respecting, or stratified vector fields). For
example, a smooth vector field $f\colon[0,1]\to\operatorname{\mathbb{R}}$ is
tangent to the boundary if it vanishes at the vertices $0$ and $1$. A smooth
vector field $f=(f_{1},f_{2})\colon[0,1]^{2}\to{\mathbb{R}}^{2}$ is tangent to
the boundary if $f_{1}(x,y)$ vanishes whenever $x=0$ or $x=1$, and
$f_{2}(x,y)$ vanishes whenever $y=0$ or $y=1$. The space of stratified vector
fields is a closed vector subspace of
$C^{\infty}(\Omega,\operatorname{\mathbb{R}}^{d})$. As in every Lie group, the
connected component $\operatorname{Diff}(\Omega)_{0}$ of the neutral element
$\operatorname{id}_{\Omega}$ is an open subgroup of
$\operatorname{Diff}(\Omega)$. We recall from [13]:

######  Lemma 3.3.

Let $\Omega=[0,1]^{d}$. Then
$\kappa\colon\operatorname{Diff}(\Omega)_{0}\rightarrow
C^{\infty}_{\mathrm{str}}(\Omega,\operatorname{\mathbb{R}}^{d}),\phi\mapsto\phi-\operatorname{id}_{\Omega}$
has open image and is a chart of $\operatorname{Diff}(\Omega)$, when
considered as a map onto the image. In particular,
$\operatorname{Diff}_{\partial}(\Omega)$ is a Lie subgroup of
$\operatorname{Diff}(\Omega)$.

Again, it is not straightforward to determine the image of the chart $\kappa$.
But we can extend the sufficient criterion encountered for boundary-fixing
diffeomorphisms. To this end, we consider the open subset

$\mathcal{V}_{1}=\\{f\in
C^{\infty}_{\text{str}}(\Omega,\operatorname{\mathbb{R}}^{d}):\text{Lip}(f)<1\\}.$

We prove the following (see [13, Proposition 1.5] for a generalized version):

######  Proposition 3.4.

For every $f\in\mathcal{V}_{1}$, the map $f+\operatorname{id}_{\Omega}$ is
contained in $\operatorname{Diff}(\Omega)$. Hence
$\mathcal{V}_{1}\subseteq\kappa(\operatorname{Diff}(\Omega)_{0})$.

###### Proof.

The proof is by induction on $d\in\mathbb{N}$. If $d=1$, then
$C^{\infty}_{\text{str}}(\Omega,\operatorname{\mathbb{R}})=C^{\infty}_{\partial}(\Omega,\operatorname{\mathbb{R}})$
and the assertion follows from Lemma 3.1. Now assume that $d\geq 2$. Fix
$f\in\mathcal{V}_{1}$. Arguing as in the proof of Lemma 3.1, we see that
$\operatorname{id}_{\Omega}+f$ is injective. Note that each
$(d-1)$-dimensional face $F$ of $\Omega$ is diffeomorphic to $[0,1]^{d-1}$ via
a translation; for example, the upper face $[0,1]\times\\{1\\}$ of the square
$[0,1]^{2}$ can be translated to $[0,1]\times\\{0\\}\cong[0,1]$. Hence
$(\operatorname{id}_{\Omega}+f)(F)=(\operatorname{id}_{F}+f|_{F})(F)=F$ by
induction. As the boundary $\partial\Omega$ is the union of the
$(d-1)$-dimensional faces, we deduce that
$(\operatorname{id}_{\Omega}+f)(\partial\Omega)=\partial\Omega$. Thus (3.3)
holds. Next, we note that $\operatorname{id}_{\Omega}+f$ has invertible
derivative $f^{\prime}(x)$ at each $x\in\Omega^{\circ}$, whence it is an open
map on the interior $\Omega^{\circ}$. To see that
$(\operatorname{id}_{\Omega}+f)(\Omega)\subseteq\Omega$, consider the
projection $\text{pr}_{j}$ onto the $j$th component for $j\in\\{1,\dots,d\\}$.
By compactness of $\Omega$, the minimum
$m_{j}:=\min_{y\in\Omega}(\text{pr}\circ(\operatorname{id}_{\Omega}+f)(y))$
exists. Since $\operatorname{id}_{\Omega^{\circ}}+f|_{\Omega^{\circ}}$ is an
open map, the minimum cannot be attained on the interior of the cube. Thus
$m_{j}$ is attained on the boundary $\partial\Omega$. As the boundary is
mapped to itself and does not contain points with negative coordinates, we
conclude that $m_{j}\geq 0$ for every $j$, whence $x_{j}+f_{j}(x)\geq 0$ for
each $x=(x_{1},\ldots,x_{d})\in\Omega$. Considering the maximum instead, we
conclude that also $x_{j}+f_{j}(x)\leq 1$, whence $x+f(x)\in\Omega$ for all
$x\in\Omega$. Thus $\operatorname{id}_{\Omega}+f$ maps the sets $\Omega$ and
$\partial\Omega$ into themselves, respectively. Arguing as in the proof of
Lemma 3.1, we see that $(\operatorname{id}_{\Omega}+f)(\Omega^{\circ})$ is
both open and closed in the connected set $\Omega^{\circ}$, whence
$(\operatorname{id}_{\Omega}+f)(\Omega^{\circ})=\Omega^{\circ}$. Then
$(\operatorname{id}_{\Omega}+f)(\Omega)=(\operatorname{id}_{\Omega}+f)(\Omega^{\circ})\cup(\operatorname{id}_{\Omega}+f)(\partial\Omega)=\Omega^{\circ}\cup\partial\Omega=\Omega$,
showing that $\operatorname{id}_{\Omega}+f$ is surjective as a map to $\Omega$
and hence bijective. The proof can now be completed as in the proof of Lemma
3.1, and we deduce that
$\operatorname{id}_{\Omega}+f\in\operatorname{Diff}(\Omega)$. ∎

We have the following (strict) inclusions concerning the sets we are working
with

${\mathcal{U}_{1}}$${\mathcal{V}_{1}}$${\kappa(\operatorname{Diff}_{\partial}(\Omega))}$${\kappa(\operatorname{Diff}(\Omega)_{0})}$${C^{\infty}(\Omega,\operatorname{\mathbb{R}}^{d}).}$

In practice, we work with vector fields to approximate diffeomorphisms. As the
algebra of vector fields is infinite dimensional, we have to replace it by a
finite-dimensional subspace. To increase the supply of diffeomorphisms we can
approximate, our method incorporates _compositions_ of diffeomorphisms. Since
$\kappa^{-1}(\mathcal{U}_{1})$ and $\kappa^{-1}(\mathcal{V}_{1})$ are identity
neighborhoods in the respective diffeomorphism groups, they generate the
group,555We are referring to the connected component of the identity of the
respective diffeomorphism group here. i.e. every diffeomorphism can be written
as a finite (but maybe arbitrarily long) composition of elements near the
identity. Passing to the intersection with a finite-dimensional subspace, the
generating property is lost, but compositions will nevertheless yield a much
broader spectrum of diffeomorphisms which can be reached.

######  Remark 3.5.

These observations have been exploited in the _boundaryless_ case for example
in [1]. Since $\operatorname{Diff}(\Omega)_{0}$ is simple in this case, there
are several characterizations available for a subset of vector fields to
generate a sufficiently rich subgroup (see e.g. [1, 2] for a treatment of
flows in the framework of sub-Riemannian geometry on manifolds _without_
boundary). Note that it is not clear how these techniques (e.g. sub-Riemannian
geometry) can be generalized to manifolds with boundary.

Since our method relies on compositions of multiple diffeomorphisms, the next
section provides insight and estimates into such compositions.

##  4 Estimates concerning multiple compositions

If we consider compositions $\phi:=\phi_{L}\circ\cdots\circ\phi_{1}$ of
diffeomorphisms
$\phi_{1},\ldots,\phi_{L}\in\operatorname{Diff}_{\partial}(\Omega)$ for a
large number $L$ of factors, we’d like to be able to control the distance
$\|\phi-\operatorname{id}_{\Omega}\|_{C^{k}}$ of the composite from the
identity $\operatorname{id}_{\Omega}$ in terms of the distances
$\|\phi_{1}-\operatorname{id}_{\Omega}\|_{C^{k}},\ldots,\|\phi_{L}-\operatorname{id}_{\Omega}\|_{C^{k}}$
(where the norms $\|\cdot\|_{C^{k}}$ are as in (4.8)). We show that such an
estimate is possible using a closed formula which applies to all $L$. Let
$\Omega$ be $[0,1]^{d}$, or a compact convex subset of
$\operatorname{\mathbb{R}}^{d}$ with dense interior. For certain positive
integers $M_{0}\leq M_{1}\leq M_{2}\leq\cdots$ which can be calculated
recursively, we have:

######  Proposition 4.1.

For all $k\in\mathbb{N}_{0}$ and all $L\in\mathbb{N}$,

$\displaystyle\|((\operatorname{id}_{\Omega}+f_{L})\circ\cdots\circ(\operatorname{id}_{\Omega}+f_{1}))-\operatorname{id}_{\Omega}\|_{C^{k}}$
(4.1) $\displaystyle\leq$ $\displaystyle
M_{k}\,e^{k\sum_{j=1}^{L}\|f_{j}\|_{C^{k}}}\,(\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}})$

holds for all
$f_{1},\ldots,f_{L}\in\operatorname{Diff}_{\partial}(\Omega)-\operatorname{id}_{\Omega}$
such that $\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}}\leq 1$. Notably,

$\|((\operatorname{id}_{\Omega}+f_{L})\circ\cdots\circ(\operatorname{id}_{\Omega}+f_{1}))-\operatorname{id}_{\Omega}\|_{C^{k}}\leq
M_{k}\,e^{k}(\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}}).$ (4.2)

Here $M_{0}:=M_{1}:=1$. For $k\in\mathbb{N}$ and $j\in\\{1,\ldots,k\\}$, let
$P_{k,j}$ be the set of all partitions $P=\\{I_{1},\ldots,I_{j}\\}$ of
$\\{1,\ldots,k\\}$ into disjoint, non-empty subsets $I_{1},\ldots,I_{j}$ of
$\\{1,\ldots,k\\}$. Then

$M_{k}:=\sum_{j=2}^{k}\sum_{P\in P_{k,j}}M_{|I_{1}|}\cdots
M_{|I_{j}|}\;\;\mbox{for $\,k\geq 2$.}$ (4.3)

If $n_{1}\geq n_{2}\geq\cdots\geq n_{j}\geq 1$ are integers with
$n_{1}+\cdots+n_{j}=k$, write $P(n_{1},\ldots,n_{j})\subseteq P_{j}$ for the
set of all partitions of $\\{1,\ldots,k\\}$ into subsets $I_{1},\ldots,I_{j}$
of cardinality $|I_{a}|=n_{a}$ for $a\in\\{1,\ldots,j\\}$. Then

$M_{k}\,=\,\sum_{j=2}^{k}\,\sum_{n_{1}+\cdots+n_{j}=k}|P(n_{1},\ldots,n_{j})|\,M_{n_{1}}\cdots
M_{n_{j}},$

using a summation over all $(n_{1},\ldots,n_{j})\in\mathbb{N}^{j}$ with
$n_{1}\geq n_{2}\geq\cdots\geq n_{j}$ and $n_{1}+\cdots+n_{j}=k$.
For small $k$, the numbers $M_{k}$ can easily be calculated by hand; their
values are as follows:
$k$ 1 2 3 4 5 6 7 8 9 10 $M_{k}$ 1 1 4 26 236 2752 39208 660032 12818912
282137824
The constants may not be optimal. Note that, for every $k\geq 3$, the set
$\\{1,\ldots,k\\}$ admits $k$ partitions into a singleton and a set with $k-1$
elements. Thus $M_{k}\geq kM_{1}M_{k-1}\geq M_{k-1}$.

######  Remark 4.2.

The numbers $M_{k}$ for $k\in\mathbb{N}$ are well known in combinatorics; they
arise in a classical enumeration problem known as Schröder’s Fourth Problem
(cf. [33, 35]). A longer list can be found in the on-line encyclopedia of
integer sequences (OEIS), sequence A000311.

######  Remark 4.3.

For $m\in\mathbb{N}$, the assertions of Proposition 4.1 remain valid for all
$k\in\\{0,1,\ldots,m\\}$ if
$f_{1},\ldots,f_{L}\in\operatorname{Diff}^{C^{m}}_{\partial}(\Omega)-\operatorname{id}_{\Omega}\subseteq
C^{m}_{\partial}(\Omega,\operatorname{\mathbb{R}}^{d})$ (the proof carries
over).

In the following,
$\mathcal{V}:=\kappa(\operatorname{Diff}_{\partial}(\Omega))=\operatorname{Diff}_{\partial}(\Omega)-\operatorname{id}_{\Omega}$,
which is an open $0$-neighborhood in
$C^{\infty}_{\partial}(\Omega,\operatorname{\mathbb{R}}^{d})$.

######  Remark 4.4.

The proof will show that the condition
$\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}}\leq 1$ in Proposition 4.1 is
unnecessary for the validity of (4.1) in the cases $k=0$ and $k=1$. Actually,

$\|((\operatorname{id}_{\Omega}+f_{L})\circ\cdots\circ(\operatorname{id}_{\Omega}+f_{1}))-\operatorname{id}_{\Omega}\|_{\infty}\leq\|f_{1}\|_{\infty}+\cdots+\|f_{L}\|_{\infty}$
(4.4)

for all $L\in\mathbb{N}$ and $f_{1},\ldots,f_{L}\in\mathcal{V}$. We shall also
see that

$\displaystyle\operatorname{Lip}(((\operatorname{id}_{\Omega}+f_{L})\circ\cdots\circ(\operatorname{id}_{\Omega}+f_{1}))-\operatorname{id}_{\Omega})$
$\displaystyle\leq$
$\displaystyle(1+\operatorname{Lip}(f_{1}))\cdots(1+\operatorname{Lip}(f_{L}))-1$
(4.5) $\displaystyle\leq$ $\displaystyle
e^{\operatorname{Lip}(f_{1})+\cdots+\operatorname{Lip}(f_{L})}-1$ (4.6)
$\displaystyle\leq$ $\displaystyle
e^{\operatorname{Lip}(f_{1})+\cdots+\operatorname{Lip}(f_{L})}(\operatorname{Lip}(f_{1})+\cdots+\operatorname{Lip}(f_{L}))$
(4.7)

for all $L\in\mathbb{N}$ and $f_{1},\ldots,f_{L}\in\mathcal{V}$.

######  Remark 4.5.

In the special case $\Omega=[0,1]^{d}$, the conclusion of Proposition 4.1
remains valid if $\operatorname{Diff}_{\partial}(\Omega)$ is replaced with
$\operatorname{Diff}(\Omega)_{0}$ and vector fields vanishing on the boundary
are replaced with vector fields which are tangent to the boundary. The same
comment applies to Remarks 4.3 and 4.4.

The following notations and facts are useful for the proof of Proposition 4.1.

######  4.1.

If $(E,\|\cdot\|_{E})$ is a normed space, we write
$\overline{B}_{r}^{E}(x):=\\{y\in E\colon\|y-x\|_{E}\leq r\\}$ for the closed
ball of radius $r>0$ around $x\in E$. If also $(F,\|\cdot\|_{F})$ is a normed
space, $k\in\mathbb{N}$ and $\beta\colon E^{k}\to F$ a continuous $k$-linear
map, we define

$\|\beta\|_{\operatorname{op}}:=\sup\\{\|\beta(x_{1},\ldots,x_{k})\|_{F}\colon
x_{1},\ldots,x_{k}\in\overline{B}_{1}^{E}(0)\\}.$

Then
$\|\beta(x_{1},\ldots,x_{k})\|_{F}\leq\|\beta\|_{\operatorname{op}}\|x_{1}\|_{E}\cdots\|x_{k}\|_{E}$
for all $x_{1},\ldots,x_{k}\in E$. If $\Omega$ is a compact topological space
and $g\colon\Omega\to(\mathcal{L}^{k}(E,F),\|\cdot\|_{\operatorname{op}})$ a
continuous map to the space of continuous $k$-linear maps, we let

$\|g\|_{\infty,\operatorname{op}}:=\sup_{x\in\Omega}\|g(x)\|_{\operatorname{op}}\in[0,\infty[\,.$

######  4.2.

Let $E:=\operatorname{\mathbb{R}}^{d}$, endowed with a fixed norm $\|\cdot\|$,
and $\Omega\subseteq E$ be a compact convex subset with non-empty interior.
For $f\in C^{\infty}(\Omega,E)$ and $k\in\mathbb{N}_{0}$, let
$d^{k}f\colon\Omega\times E^{k}\to E$ be the continuous map such that, for all
$x\in\Omega^{\circ}$ and $y_{1},\ldots,y_{k}\in E$,

$d^{k}f(x,y_{1},\ldots,y_{k}):=(D_{y_{k}}\cdots D_{y_{1}}f)(x)$

is the iterated directional derivative at $x$ in the directions
$y_{1},\ldots,y_{k}$. Recall that the topology on
$C^{\infty}_{\partial}(\Omega,E)$ is initial with respect to the mappings
$C^{\infty}_{\partial}(\Omega,E)\to C(\Omega\times E^{k},E)$, $f\mapsto
d^{k}f$ for $k\in\mathbb{N}_{0}$, using the topology of compact convergence on
spaces of continuous functions. Thus, the seminorms

$f\mapsto\sup\\{\|d^{k}f(z)\|_{E}\colon z\in K\\}$

define the locally convex vector topology on
$C^{\infty}_{\partial}(\Omega,E)$, for $k\in\mathbb{N}_{0}$ and compact
subsets $K\subseteq\Omega\times E^{k}$. Abbreviate
$f^{(k)}(x):=d^{k}f(x,\cdot)\in\mathcal{L}^{k}(E,E)$ for $k\in\mathbb{N}$ and
$x\in\Omega$. Then

$f\mapsto\|f^{(k)}\|_{\infty,\operatorname{op}}=\sup\\{\|d^{k}f(z)\|_{E}\colon
z\in\Omega\times\overline{B}^{E}_{1}(0)^{k}\\}$

is a continuous seminorm on $C^{\infty}_{\partial}(\Omega,E)$. Together with
the supremum norm $\|\cdot\|_{\infty}$, the latter seminorms define the
locally convex topological vector space topology of
$C^{\infty}_{\partial}(\Omega,E)$. In fact, for each $k\in\mathbb{N}$ and
compact subset $K\subseteq\Omega\times E^{k}$, we have
$K\subseteq\Omega\times\overline{B}^{E}_{r}(0)^{k}$ for some $r>0$, whence

$\sup\\{\|d^{k}f(x,z)\|_{E}\colon z\in K\\}\leq
r^{k}\|f^{(k)}\|_{\infty,\operatorname{op}}.$

Hence, the vector topology on $C^{\infty}_{\partial}(\Omega,E)$ is also
defined by the sequence of norms $(\|\cdot\|_{C^{k}})_{k\in\mathbb{N}_{0}}$
with $\|\cdot\|_{C^{0}}:=\|\cdot\|_{\infty}$ and

$\|f\|_{C^{k}}:=\max\\{\|f\|_{\infty},\|f^{(1)}\|_{\infty,\operatorname{op}},\ldots,\|f^{(k)}\|_{\infty,\operatorname{op}}\\}.$
(4.8)

The following formula actually holds for locally convex spaces $X$, $Y$, and
$Z$.

######  4.3.

(Faá di Bruno’s Formula) Let $X$, $Y$, and $Z$ be finite-dimensional real
vector spaces, $U\subseteq X$ and $V\subseteq Y$ be locally convex subsets
with dense interior, and $k\in\mathbb{N}$. Let $g\colon U\to Y$ and $f\colon
V\to Z$ be $C^{k}$-maps such that $g(U)\subseteq V$. Then

$d^{k}(f\circ g)(x,y_{1},\ldots,y_{k})=\sum_{j=1}^{k}\sum_{P\in
P_{k,j}}d^{j}f(d^{|I_{1}|}g(x,y_{I_{1}}),\ldots,d^{|I_{j}|}g(x,y_{I_{j}}))$

for all $x\in U$ and $y_{1},\ldots,y_{k}\in E$, where
$P=\\{I_{1},\ldots,I_{j}\\}$ and

$y_{J}:=(y_{i_{1}},\ldots,y_{i_{\ell}})$

for each subset $J\subseteq\\{1,\ldots,k\\}$ of the form
$J=\\{i_{1},\ldots,i_{\ell}\\}$ with $i_{1}<\cdots<i_{\ell}$ (see [16, Theorem
1.3.18 and Remark 1.4.15 (c)]; cf. also [10]).

######  4.4.

We shall use the following fact: If $L\in\mathbb{N}$ and
$x_{1},\ldots,x_{L}\geq 0$, then

$(1+x_{1})(1+x_{2})\cdots(1+x_{L})\leq e^{x_{1}+\cdots+x_{L}}$ (4.9)

as
$\ln((1+x_{1})\cdots(1+x_{L}))=\sum_{j=1}^{L}\ln(1+x_{j})\leq\sum_{j=1}^{L}x_{j}$.
We shall also use that

$e^{x}-1\leq e^{x}\,x\;\;\mbox{for all $x\geq 0$.}$ (4.10)

In fact, $e^{x}-1=\exp(x)-\exp(0)=\exp^{\prime}(\xi)x$ for some $\xi\in[0,x]$
by the Mean Value Theorem, with $\exp^{\prime}(\xi)=\exp(\xi)\leq\exp(x)$ by
monotonicity of $\exp$.

######  4.5.

Let $(E,\|\cdot\|)$ be a normed space and
$(\mathcal{L}(E),\|\cdot\|_{\operatorname{op}})$ be the corresponding algebra
of bounded operators. For all $L\in\mathbb{N}$ and
$\alpha_{1},\ldots,\alpha_{L}\in\mathcal{L}(E)$, we trivially have

$\|(\operatorname{id}_{E}+\alpha_{L})\circ\cdots\circ(\operatorname{id}_{E}+\alpha_{1})\|_{\operatorname{op}}\leq(1+\|\alpha_{L}\|_{\operatorname{op}})\cdots(1+\|\alpha_{1}\|_{\operatorname{op}}).$
(4.11)

Moreover,

$\displaystyle\|((\operatorname{id}_{E}+\alpha_{L})\circ\cdots\circ(\operatorname{id}_{E}+\alpha_{1}))-\operatorname{id}_{E}\|_{\operatorname{op}}$
$\displaystyle\leq$
$\displaystyle((1+\|\alpha_{L}\|_{\operatorname{op}})\cdots(1+\|\alpha_{1}\|_{\operatorname{op}}))-1$
(4.12) $\displaystyle\leq$ $\displaystyle
e^{\|\alpha_{1}\|_{\operatorname{op}}+\cdots+\|\alpha_{L}\|_{\operatorname{op}}}-1$
(4.13) $\displaystyle\leq$ $\displaystyle
e^{\|\alpha_{1}\|_{\operatorname{op}}+\cdots+\|\alpha_{L}\|_{\operatorname{op}}}(\|\alpha_{1}\|_{\operatorname{op}}+\cdots+\|\alpha_{L}\|_{\operatorname{op}})$
(4.14)

In fact,

$\beta:=(\operatorname{id}_{E}+\alpha_{1})\cdots(\operatorname{id}_{E}+\alpha_{L})-\operatorname{id}_{E}=\sum_{j=1}^{m}\sum_{|I|=j}\alpha_{i_{1}}\circ\cdots\circ\alpha_{i_{j}},$

using subsets $I=\\{i_{1},\ldots,i_{j}\\}\subseteq\\{1,\ldots,L\\}$ with
$i_{1}<\cdots<i_{j}$. Thus

$\|\beta\|_{\operatorname{op}}\leq\sum_{j=1}^{L}\sum_{|I|=j}\|\alpha_{i_{1}}\|_{\operatorname{op}}\cdots\|\alpha_{i_{j}}\|_{\operatorname{op}}=((1+\|\alpha_{1}\|_{\operatorname{op}})\cdots(1+\|\alpha_{L}\|_{\operatorname{op}}))-1,$

establishing (4.12). Using (4.9) and (4.10), the inequalities (4.13) and
(4.14) follow.

###### Proof of Proposition 4.1.

We first prove the assertion for $k\in\\{0,1\\}$.
_The case $k=0$_. For $f,g\in\mathcal{V}$, let

$f*g:=((\operatorname{id}_{\Omega}+f)\circ(\operatorname{id}_{\Omega}+g))-\operatorname{id}_{\Omega}=g+f\circ(\operatorname{id}_{\Omega}+g).$

Then $f*g\in\mathcal{V}$ and

$\|f*g\|_{\infty}\leq\|g\|_{\infty}+\|f\|_{\infty}.$

Given $f_{1},\ldots,f_{L}\in\mathcal{V}$, we define $f_{L}*\cdots*f_{1}$
recursively as $f_{L}*(f_{L-1}*\cdots*f_{1})$ and obtain

$\|f_{L}*\cdots*f_{1}\|_{\infty}\leq\|f_{1}\|_{\infty}+\cdots+\|f_{L}\|_{\infty},$
(4.15)

by a straightforward induction. By induction on $2\leq L\in\mathbb{N}$, also

$((\operatorname{id}_{\Omega}+f_{L})\circ\dots\circ(\operatorname{id}_{\Omega}+f_{1}))-\operatorname{id}_{\Omega}=f_{L}*\cdots*f_{1},$

as the left-hand side of the equation equals
$(\operatorname{id}_{\Omega}+f_{L})(\operatorname{id}_{\Omega}+((\operatorname{id}_{\Omega}+f_{L-1})\circ\cdots\circ(\operatorname{id}_{\Omega}+f_{1})-\operatorname{id}_{\Omega}))-\operatorname{id}_{\Omega}$$=f_{L}*(((\operatorname{id}_{\Omega}+f_{L-1})\circ\cdots\circ(\operatorname{id}_{\Omega}+f_{1}))-\operatorname{id}_{\Omega})=f_{L}*(f_{L-1}*\cdots*f_{1})$.
Hence (4.4) is a re-writing of (4.15). Notably, (4.1) holds for $k=0$ with
$M_{0}:=1$.
_The case $k=1$_. Abbreviate $E:=\operatorname{\mathbb{R}}^{d}$. For
$f_{1},\ldots,f_{L}\in\mathcal{V}$, $x_{0}\in\Omega$ and
$x_{j}:=f_{j}(x_{j-1})$ for $j\in\\{1,\ldots,L-1\\}$, we have

$(((\operatorname{id}_{\Omega}+f_{L})\circ\cdots\circ(\operatorname{id}_{\Omega}+f_{1}))-\operatorname{id}_{\Omega})^{\prime}(x_{0})=((\operatorname{id}_{E}+f_{1}^{\prime}(x_{0}))\circ\cdots\circ(\operatorname{id}_{E}+f_{L}^{\prime}(x_{L-1})))-\operatorname{id}_{E}$

with operator norm
$\leq((1+\|f_{1}^{\prime}(x_{0})\|_{\operatorname{op}})\cdots(1+\|f_{L}^{\prime}(x_{L-1})\|_{\operatorname{op}}))-1$,
by (4.12). Passing to the supremum in $x_{0}\in\Omega$, we get

$\|(f_{L}*\cdots*f_{1})^{\prime}\|_{\infty,\operatorname{op}}\leq((1+\|f_{1}^{\prime}\|_{\infty,\operatorname{op}})\cdots(1+\|f_{L}^{\prime}\|_{\infty,\operatorname{op}}))-1,$

which is bounded by $((1+\|f_{1}\|_{C^{1}})\cdots(1+\|f_{L}\|_{C^{1}}))-1$.
Thus (4.5) holds, and (4.6) as well as (4.7) follow using (4.9) and (4.10). As
the latter number is also an upper bound for
$\|f_{L}*\cdots*f_{1}\|_{\infty}$, we deduce that

$\|f_{L}*\cdots*f_{1}\|_{C^{1}}\leq((1+\|f_{1}\|_{C^{1}})\cdots(1+\|f_{L}\|_{C^{1}}))-1$

for all $L\in\mathbb{N}$ and $f_{1},\ldots,f_{L}\in\Omega$. By (4.9) and
(4.10), the right-hand side is bounded by

$e^{\|f_{1}\|_{C^{1}}+\cdots+\|f_{L}\|_{C^{1}}}-1\leq
e^{\|f_{1}\|_{C^{1}}+\cdots+\|f_{L}\|_{C^{1}}}(\|f_{1}\|_{C^{1}}+\cdots+\|f_{L}\|_{C^{1}}),$

for all $L\in\mathbb{N}$ and $f_{1},\ldots,f_{L}\in\mathcal{V}$. Notably,
(4.1) holds for $k=1$ with $M_{1}:=1$.
For $k\geq 2$, the proof is by induction. Details are postponed to Appendix A.
∎

The upper bounds of proposition 4.1 are theoretical. It could therefore be of
interest to see how they relate to the practically achieved $C^{k}$-norms of
diffeomorphisms of the form used in section 6 and how these bounds are
affected by for example the number of compositions. For simplicity, we stick
to the one-dimensional case. In figure 2, the value

$\frac{\|((\operatorname{id}_{\Omega}+f_{L})\circ\cdots\circ(\operatorname{id}_{\Omega}+f_{1}))-\operatorname{id}_{\Omega}\|_{C^{k}}}{\left(\sum_{\ell}^{L}\|f_{\ell}\|_{C^{k}}\right)e^{k\sum_{\ell}^{L}\|f_{\ell}\|_{C^{k}}}}$
(4.16)

is shown for lower-order derivatives with varying number of layers $L$ and
number of basis functions $d$ per layer. For $k=1,\,2$ the estimates are
relatively good with values close to $60\%$ and $15\%$ of $M_{k}$
respectively.

These estimates are largely dependent on the choice of basis functions, so the
results may look different for another set of basis functions than the ones
used here. Further details regarding how the values of figure 2 are computed
can be found in appendix B.

Figure 2: Practical estimates of the $C^{k}$ bounding coefficient in equation
(4.1) for $k=1,2,3$ (left to right). (Top): Varying number of layers $L$.
(Bottom): Varying number of basis functions per layer $M$

##  5 Implementation of the methods

This section describes the implementation of the approach outlined in section
2.4. We first consider the case of curves and then the case of surfaces. In
both cases, we need three parts: a discrete problem formulation, a choice of
basis functions, and a projection operator to ensure that the Lipschitz
constant of each layer is bounded by one.

###  5.1 Curves

##### Problem Formulation

Denote by
$W=\\{\mathbf{w}^{\ell}\\}_{\ell=1}^{L}\subset\operatorname{\mathbb{R}}^{d}$
the collection of weight vectors $\mathbf{w}^{\ell}$, and use the notation

$\varphi(x;W)=\varphi_{L}(x;\mathbf{w}^{L})\circ...\circ\varphi_{1}(x;\mathbf{w}^{1})$
(5.1)

to highlight the dependency of $\varphi$ on $W$, and $\varphi_{\ell}$ on
$\mathbf{w}^{\ell}$, respectively. The goal is to find a set of vectors $W$
such that the function $\varphi(\cdot,W)\in\text{Diff}^{+}(\Omega)$ minimizes
the loss function (2.6). To approximate the shape distance, take a collection
$X=\\{x_{k}\\}_{k=1}^{K}\subset\Omega$ of linearly spaced points, and compute
the mean squared error (MSE) between the $q$-maps evaluated at these points.
In other words, the loss function is approximated by the function

$E(W;X)=\frac{1}{K}\sum_{k=1}^{K}\left|q(x_{k})-\sqrt{\varphi^{\prime}(x_{k};W)}r(\varphi(x_{k};W))\right|^{2}.$
(5.2)

This function should be optimized with respect to the weights $W$ under the
constraints that $\varphi(0,W)=0$, $\varphi(1,W)=1$, and
$\varphi^{\prime}(\cdot:W)>0$. The optimization is handled by a BFGS-
optimizer, with a projection step between the weight updates to ensure that
the layers are invertible.666This approach has performed well for the
experiments in this article, but we have found limited theory on the
convergence of BFGS with a projection step. A large scale application of the
reparametrization algorithm might gain from changing to a constrained
optimization routine, or to another projected quasi-Newton algorithm with a
better theoretical foundation. The projection step is similar to the spectral
normalization approach of [4], and entails a scaling of the weight vectors
such that the Lipschitz constants of each of the residual maps
$x\mapsto\sum_{n=1}^{M}w^{\ell}_{n}f_{n}(x)$ are bounded by one.

##### Basis Functions

The tangent space of the diffeomorphism group on an interval consists of
smooth functions that vanish on the boundary of the interval. A simple basis
which satisfies this relation is a truncated Fourier sine series. However, by
scaling the sine series such that their derivatives are bounded by 1, we
achieve that the size of the feasible set is similar in each direction.
Therefore, choose a basis of functions

$f_{n}(x)=\frac{\sin(n\pi x)}{n\pi},\quad n=1,...,M$ (5.3)

for some $M\in\mathbb{N}$.

##### Projection Operator

To ensure the invertibility of each layer $\varphi_{\ell}$, we require that
the Lipschitz constant of each of the residual maps to be smaller than one.
For this purpose, denote by $\mathcal{W}^{\varepsilon}$ the feasible set of
the weight vectors. To define this set explicitly, consider the upper estimate
of the Lipschitz constants given by the inequalities

$\left|\left(\sum_{n=1}^{M}w_{n}f_{n}(x)\right)-\left(\sum_{n=1}^{M}w_{n}f_{n}(y)\right)\right|\leq\underbrace{\left(\sum_{n=1}^{M}|w_{n}|L_{n}\right)}_{:=L}|x-y|,$
(5.4)

where $L_{n}$ is the Lipschitz constant for the basis function $f_{n}$. For
our choice of basis functions, we have that
$L_{n}=\sup_{x\in\Omega}|f_{n}^{\prime}(x)|=1$ for all $\;n=1,\,...,M$. Thus
$L=\sum_{n=1}^{M}|w_{n}|L_{n}=\|\mathbf{w}\|_{1}$. To ensure that the layers
are diffeomorphisms, we require that $L$ is strictly smaller than one. To
achieve this, we define the feasible set

$\mathcal{W}^{\varepsilon}=\left\\{\mathbf{w}\in\operatorname{\mathbb{R}}^{M}\;:\;\|\mathbf{w}\|_{1}\leq
1-\varepsilon\right\\}.$ (5.5)

for some small $\varepsilon$. If the optimization algorithm makes a weight
update which causes a violation of these constraints, then it needs to be
projected back onto the feasible set. We define this projection through a
simple scaling of the vector, defined by

$\pi\colon\operatorname{\mathbb{R}}^{M}\to\mathcal{W}^{\varepsilon},\quad\pi(\mathbf{w})=\frac{1-\varepsilon}{\max\\{1-\varepsilon,\|\mathbf{w}\|_{1}\\}}\mathbf{w}.$
(5.6)

###  5.2 Surfaces

##### Problem Formulation

The problem formulation for reparametrization of surfaces is analogous to the
case for curves, with the natural extensions to two dimensions: Approximate
the continuous loss function by

$E(W;X)=\frac{1}{K^{2}}\sum_{i,j=1}^{K}\left|q(\mathbf{x}_{i,j})-\sqrt{J_{\varphi}(\mathbf{x}_{i,j};W)}r\left(\varphi(\mathbf{x}_{i,j};W)\right)\right|^{2},$
(5.7)

where
$X=\\{\mathbf{x}_{i,j}\,,\,i,j=1,...,K\\}\subset\Omega(\subset\operatorname{\mathbb{R}}^{2})$,
and $\varphi$ is given on the form (5.1). The loss function should be
optimized with respect to the weights $W$ under the constraints that
$J_{\varphi}(\cdot:W)>0$ and $\varphi(\partial M;W)=\partial M$.

#### Basis Functions

The tangent space of the diffeomorphism group on $\Omega=[0,1]^{2}$ consists
of smooth vector fields such that the vectors are tangential to the boundary
$\partial\Omega$. We construct a basis for this space componentwise, as a
tensor product basis of Fourier sine series in one direction and a full
Fourier basis in the other. This choice was made such as to be as close to the
gradient descent algorithm of [25] which inspired the methods of this article.
In the first component, these functions are given by the three families of
functions,

$\xi_{k}(x,y)=\frac{\sin(\pi kx)}{\pi k},\quad\eta_{k,l}(x,y)=\frac{\sin(\pi
kx)\cos(2\pi ly)}{\pi kl},\quad\phi_{k,l}(x,y)=\frac{\sin(\pi kx)\sin(2\pi
ly)}{\pi kl},$ (5.8)

with corresponding basis functions with the arguments flipped in the second
component. A few example basis vector fields are shown in figure 3. The basis
is truncated such that the coefficients $k,l=1,...,N$ for some maximal
“frequency” $N\in\mathbb{N}$. The total number of basis functions is therefore
$M=2(2N^{2}+N)$, which means that the number of basis functions per layer in
the network will increase quickly when increasing the maximal frequency of the
basis functions.

(a) $\xi_{1}$

(b) $\varphi_{1,1}$

(c) $\tilde{\eta}_{1,1}$

(d) $\tilde{\eta}_{1,3}$

Figure 3: Four examples of basis elements of
$T_{\operatorname{id}}(\text{Diff}^{+}(\Omega))$.

#### Projection Operator

Similarly as in section 5.1, the invertibility of each layer $\varphi_{\ell}$
is ensured by bounding the Lipschitz constants of the residual maps by one.
The projection operator is analogous to the one-dimensional case, and is based
on the approximation (5.4). However, the Lipschitz constant $L_{n}$ of each of
the basis vector fields is no longer simply equal to $1$. Instead, we make use
of the following approximation:

$L_{n}=\sup_{x\in\Omega}\|Df_{n}(x)\|_{2}\leq\sup_{x\in\Omega}\|Df_{n}(x)\|_{F},\quad
n=1,\,...,M,$

where $Df_{n}$ denotes the Jacobian matrix of $f_{n}$, the symbol
$\|\cdot\|_{2}$ denotes the spectral norm, and $\|\cdot\|_{F}$ the Frobenius
norm. According to this estimate, the Lipschitz constant for each of the three
types of basis vector fields are given by

$L^{\xi}_{k}=1,\quad L^{\eta}_{k,l}=\frac{\sqrt{k^{2}+2l^{2}}}{kl},\quad
L^{\phi}_{k,l}=\frac{\sqrt{k^{2}+2l^{2}}}{kl}.$ (5.9)

Similar constants exist for the corresponding vector fields in the
$y$-direction. Adopting the notations from the one-dimensional case, we define

$\pi\colon\operatorname{\mathbb{R}}^{M}\to\mathcal{W}^{\varepsilon},\quad\pi(\mathbf{w})=\frac{1-\varepsilon}{\max\left\\{1-\varepsilon,\left(\sum_{n=1}^{M}|w_{n}|L_{n}\right)\right\\}}\mathbf{w},$
(5.10)

which maps infeasible weight vectors to the feasible set.

###  5.3 Shape Interpolation

An interesting application of the reparametrization algorithm is to create
realistic transitions between different shapes. The shape distance corresponds
to the length of some geodesic between the two parametric shape
representations in the pre-shape space $\mathcal{P}$. Each of the points along
this geodesic corresponds to a parametric representation of some intermediate
shape, providing a natural transformation between the original curves or
surfaces.

In the specific case of reparametrization of curves using the SRVT, we may
compute a geodesic between their shapes. Assume for any $\tau\in(0,1)$ that
the curve $\tau\mathcal{R}(c_{1})(x)+(1-\tau)\mathcal{R}(c_{2})(x)\neq
0,\,\forall\,x\in\Omega$. Then the map

$\tau\mapsto\mathcal{R}^{-1}(\tau\mathcal{R}(c_{1})+(1-\tau)\mathcal{R}(c_{2}))$
(5.11)

with

$\mathcal{R}^{-1}(q)(x)=\int_{0}^{x}q(\xi)|q(\xi)|\,d\xi$ (5.12)

defines a geodesic between the two curves $c_{1}$ and $c_{2}$.

The surface transforms, however, are not invertible. Therefore, we adopt an
alternative approach which seems to work well in many cases. Define the map

$\gamma:\mathcal{P}\times\mathcal{P}\times[0,1]\to
C^{\infty}(\Omega,\operatorname{\mathbb{R}}^{n}),\quad\gamma(f_{1},f_{2},\tau)=\tau
f_{1}+(1-\tau)f_{2}.$

Note that for $f_{1},f_{2}\in\mathcal{P}$, the map
$\tau\mapsto\gamma(f_{1},f_{2},\tau)$ defines a convex combination of two
elements in $C^{\infty}(\Omega,\operatorname{\mathbb{R}}^{n})$ connecting the
two surfaces. Applied directly to the parametrized shapes, the transition
between the two seem “irregular”. To improve upon this irregular shape
transformation, first reparametrize one of the shapes by
$\varphi^{*}\in\operatorname*{arg\,min}_{\varphi^{+}\in\operatorname{Diff}(\Omega)}d_{\mathcal{P}}([f_{1}],[f_{2}\circ\varphi])$
and then apply the shape interpolation according to

$\tau\mapsto\gamma(f_{1},f_{2}\circ\varphi^{*},\tau).$ (5.13)

(a)

(b)

Figure 4: Different interpolation methods between two shapes and how these
curves relate to different shape representations. The dashed red curve
corresponds to the direct linear interpolation of between the two original
parametric shapes. The solid green curve corresponds to a geodesic, whereas
the dotted blue line represents a linear interpolation after reparametrization
given in (5.13).

Figure 4 illustrates the different interpolation approaches, and figure 9
presents an example of the approach (5.13) applied to images.

Figure 5: Comparison of the methods described in section 5.3 interpolation
between a half-circle in white and a circle in black. (Left) Direct linear
interpolation, corresponding to the dashed red curve in 5.3. (Middle) Linear
interpolation after reparametrization, corresponding to the dotted blue curve.
(Right) Geodesic, corresponding to the solid green curve.

###  5.4 Implementation

Since the structure of the diffeomorphisms that we are searching for is
similar to that of a residual neural network, the weight optimization is
analogous to the process of training a neural network. This is typically done
by an iterative line search method such as gradient descent rather than
Newton’s method or other Quasi-Newton algorithms. However, since we are
dealing with a deterministic problem, we have chosen to use the BFGS-algorithm
(see e.g. [30]). The network is implemented using the PyTorch machine learning
framework in Python[32] which takes care of the computation of gradients and
updates the weight vectors. By building our algorithms upon the PyTorch-
framework, the algorithms easily achieve high performance, scalability and
extensibility.

Another noteworthy aspect of this approach is that the final result is a well-
defined diffeomorphism, in the form of a neural network. This network sets no
restrictions on the point collection $X$ which is used in the optimization
step, beyond that each point must belong to the domain. The number of points
in the collection, the point placement, and the ordering may be chosen as the
user prefers. This may be useful if the reparametrized surface will be used in
another application which requires the surface to be represented by a specific
number of points. Moreover, it enables some possibly useful options for the
optimization. For example, replacing or augmenting the fixed grid of points
that have been used within this article with a random sample of points within
the domain might help the optimization algorithm to escape local minima.

##  6 Numerical Experiments

###  6.1 Curves

We start by comparing two different parametric curves representing the same
shape. A simple way to find two such curves is by first selecting a curve $c$,
and reparametrizing the curve with some diffeomorphism $\varphi$. By applying
the reparametrization algorithm to the curves $c_{1}=c\circ\varphi$ and
$c_{2}=c$, we expect to recover a diffeomorphism $\psi\approx\varphi$ such
that the distance $d_{\mathcal{P}}(c_{1},c_{2}\circ\psi)\approx 0$. As a first
example, consider the parametric curve and diffeomorphism given by

$c(t)=[\cos(2\pi t),\sin(4\pi
t)],\qquad\varphi(t)=\frac{\log(20t+1)}{2\log(21)}+\frac{1+\tanh(20(t-0.5))}{4\tanh(10)}.$

The curves $c$ and $c\circ\varphi$ are shown in figure 6. The dots along the
curves represent linearly spaced points on the domain $\Omega=[0,1]$ for the
curve parameter $t$. Figure 6 presents the results of applying the
reparametrization algorithm to these curves using a network with $L=10$ layers
and $M=10$ basis functions per layer. The resulting reparametrization $\psi$
is visually indistinguishable from the true diffeomorphism $\varphi$, and the
components of the reparametrized curve $c\circ\psi$ exactly match the
components of the components of $c\circ\varphi$.

Figure 6: (Left) The curve $c$ to be reparametrized on top, and the target
curve $c\circ\varphi$ below. (Middle) The component functions of the two
curves $c\circ\varphi$ (dotted blue) and $c$ (solid orange) before, and after
reparametrization. (Top right) The true reparametrization $\varphi$ (dotted
black) compared to the reparametrization $\psi$ found by the algorithm.
(Bottom right): Convergence plot for the error as a function of iteration
number, given relative to the initial error.

###  6.2 Surfaces

Similarly as in section 6.1, the performance of the reparametrization
algorithm is investigated by comparing surfaces representing the same shape.
The results are presented in figures 7, 8. The coloring on the surfaces
represents the local area scaling factor $a_{f}$ for the given surface. Once
again, the algorithm does indeed find a diffeomorphism $\psi\approx\varphi$.
Moreover, based on figure 11, the convergence behavior of the
reparametrization algorithm under an increasing number of basis functions and
network layers is similar to the behavior already observed for curves.

The previous examples are using “artificial” surfaces defined by explicit
parametric functions. In most practical applications, however, objects are
represented by discrete observations. For example, images are typically stored
as matrices representing image intensities. Such data must be transformed to
parametrized surfaces before reparametrization. For this purpose, we made use
of a function777torch.nn.functional.grid_sample from the PyTorch library which
implements a bicubic convolution algorithm, [22]. Since we here have chosen to
work with three-dimensional surfaces, it is natural to represent single-
channel images by their graphs.

The fact that the optimal reparametrization is unknown when comparing surfaces
representing different shapes, implies that it is difficult to properly verify
the correctness of the algorithm. However, the shape interpolation method of
section 5.3 provides a good alternative to validate the results; an optimal
reparametrization of one surface relative to the other, should provide a
smoother interpolation between the two. Figure 9 illustrates the interpolation
between two images before and after reparametrization. Before
reparametrization, the intermediate images clearly show elements from each of
the images simultaneously, as one fades away and the other is revealed. After
reparametrization the interpolation represents a smoother transition, with one
of the surfaces being “warped” into the other.

(a) $f\circ\varphi$

(b) $f$

(c) $f\circ\psi$

(d) $\varphi$

(e) $\psi$

Figure 7: Reparametrization of two surfaces representing the same shape

(a) $f\circ\varphi$

(b) $f$

(c) $f\circ\psi$

(d) $\varphi$

(e) $\psi$

Figure 8: Reparametrization of two surfaces representing the same shape Figure
9: Linear interpolation between the two images from the MNIST digits dataset,
before and after reparametrization of the surfaces. The first row corresponds
to direct linear interpolation between the graphs, whereas the second row
corresponds to linear interpolation after reparametrization as given in
(5.13).

###  6.3 Comparing Network Sizes

Figure 10 presents the results of the reparametrization algorithm for a
varying number of layers and basis functions per layer in the network. The
results are also based on the curves in figure 6. These results show that
increasing the number of basis functions, while the number of layers remains
small, does not result in a significant decrease of the shape distance. A
similar result is obtained when keeping the number of basis functions constant
and small, while increasing the number of layers. Moreover, the results
provide evidence that diffeomorphisms may be approximated with high precision
by using relatively few coefficients, when using function compositions.

Figure 11 presents similar results based on the shapes in figure 7. Note that
$N$ represents the “largest frequency” of the basis functions, and not the
number of basis functions. This number is much higher as pointed out in
section 5.2. In other words, a deep network will be able to achieve similar
performance as a wide network, using much fewer coefficients.

Figure 10: Final error after reparametrization of the curves in figure 6, for
varying number of layers $L$ and basis functions per layer $M$.

Figure 11: Final error after reparametrization of the surfaces in figure 7 for
varying number of layers $L$ and basis functions per layer $M$.

##  7 Conclusion

In this paper we proposed a method for the solution of optimization problems
on the group of diffeomorphisms. The method is inspired by deep neural
networks and implemented in PyTorch. An analysis of the universal
approximation properties and a priori estimates of the norms of the obtained
approximations are presented. A number of numerical results illustrate the
merit of the proposed approach in applications of shape analysis.

## Acknowledgements

EC has received funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement No
860124.

##  Appendix A Proof of Proposition 4.1

In this appendix, we provide the details postponed in the proof of 4.1. Let us
recall the statement for the reader’s convenience:

######  A.1Statement of Proposition 4.1.

For all $k\in\mathbb{N}_{0}$ and all $L\in\mathbb{N}$,

$\|((\operatorname{id}_{\Omega}+f_{L})\circ\cdots\circ(\operatorname{id}_{\Omega}+f_{1}))-\operatorname{id}_{\Omega}\|_{C^{k}}\leq
M_{k}\,e^{k\sum_{j=1}^{L}\|f_{j}\|_{C^{k}}}\,(\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}})$

holds for all
$f_{1},\ldots,f_{L}\in\operatorname{Diff}_{\partial}(\Omega)-\operatorname{id}_{\Omega}$
such that $\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}}\leq 1$.

Recall that in Section 4 we have already established the statement for all
$L\in\mathbb{N}$ and $k\in\\{0,1\\}$. Thus all we need to prove is the
statement for all $k\geq 2$.

###### Proof of Proposition 4.1 in the case $k\geq 2$.

The proof is by induction. Assume that $k\geq 2$ and assume that the asserted
estimates already hold for $k-1$ in place of $k$. For $f,g\in\mathcal{V}$,
applying Faá di Bruno’s formula to the first term in
$f*g=(\operatorname{id}_{\Omega}+f)\circ(\operatorname{id}_{\Omega}+g)-\operatorname{id}_{\Omega}$,
we get for all $x\in\Omega$ and $y_{1},\ldots,y_{k}\in\overline{B}^{E}_{1}(0)$

$\displaystyle d^{k}(f*g)(x,y_{1},\ldots,y_{k})$ $\displaystyle=$
$\displaystyle\sum_{j=2}^{k}\sum_{P\in
P_{k,j}}d^{j}f(h(x),d^{|I_{1}|}h(x,y_{I_{1}}),\ldots,d^{|I_{j}|}h(x,y_{I_{j}}))$
(A.2)
$\displaystyle+\,(\operatorname{id}_{E}+f^{\prime}(x))(d^{k}g(x,y_{1},\ldots,y_{k})),$

with $h:=\operatorname{id}_{\Omega}+g$. The norm of the summand in (A.2) is
bounded by

$(1+\|f^{\prime}\|_{\infty,\operatorname{op}})\|g^{(k)}\|_{\infty,\operatorname{op}}\leq
e^{\|f\|_{C^{k}}}\|g\|_{C^{k}};$ (A.3)

this is $\leq M_{k}\,e^{k(\|f\|_{C^{k}}+\|g\|_{C^{k}})}\|g\|_{C^{k}}$. The
norm of a summand in (A.2) is

$\leq\|f^{(j)}\|_{\infty,\operatorname{op}}\|h^{(|I_{1}|)}\|_{\infty,\operatorname{op}}\cdots\|h^{(|I_{j}|)}\|_{\infty,\operatorname{op}}\leq\|f\|_{C^{k}}\|h^{(|I_{1}|)}\|_{\infty,\operatorname{op}}\cdots\|h^{(|I_{j}|)}\|_{\infty,\operatorname{op}}.$

For $a\in\\{1,\ldots,j\\}$, possible cases are that $|I_{a}|=1$; then

$\|h^{(|I_{a}|)}\|_{\infty,\operatorname{op}}=\|\operatorname{id}_{E}+g^{\prime}\|_{\infty,\operatorname{op}}\leq
1+\|g^{\prime}\|_{\infty,\operatorname{op}}\leq
e^{\|g\|_{C^{k}}}=e^{|I_{a}|\,\|g\|_{C^{k}}}.$

If $|I_{a}|\geq 2$, then
$\|h^{(|I_{a}|)}\|_{\infty,\operatorname{op}}=\|g^{(|I_{a}|)}\|_{\infty,\operatorname{op}}\leq\|g\|_{C^{k}}\leq
e^{\|g\|_{C^{k}}}\leq e^{|I_{a}|\,\|g\|_{C^{k}}}$. As
$e^{|I_{1}|\,\|g\|_{C^{k}}}\cdots
e^{|I_{j}|\,\|g\|_{C^{k}}}=e^{k\|g\|_{C^{k}}}\leq M_{|I_{1}|}\cdots
M_{|I_{j}|}e^{k(\|f\|_{C^{k}}+\|g\|_{C^{k}})}$, we get

$\displaystyle\|(f*g)^{(k)}\|_{\infty,\operatorname{op}}$ $\displaystyle\leq$
$\displaystyle\left(\sum_{j=2}^{k}\sum_{P\in P_{k,j}}M_{|I_{1}|}\cdots
M_{|I_{j}|}\right)\|f\|_{C^{k}}e^{k(\|f\|_{C^{k}}+\|g\|_{C^{k}})}$
$\displaystyle+\,M_{k}\,e^{k(\|f\|_{C^{k}}+\|g\|_{C^{k}})}\|g\|_{C^{k}}$
$\displaystyle\leq$ $\displaystyle
M_{k}\,e^{k(\|f\|_{C^{k}}+\|g\|_{C^{k}})}(\|f\|_{C^{k}}+\|g\|_{C^{k}}).$

Thus $\|f*g\|_{C^{k}}\leq
M_{k}\,e^{k(\|f\|_{C^{k}}+\|g\|_{C^{k}})}(\|f\|_{C^{k}}+\|g\|_{C^{k}})$,
whence (4.1) holds for $L=2$. If $f_{1},\ldots,f_{L}\in\mathcal{V}$ are given
with $L\geq 3$ and

$\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}}\leq 1,$ (A.4)

let $f:=f_{L}$ and $g:=f_{L-1}*\cdots*f_{1}$. Using (A.3) and the inductive
hypothesis (the case of $L-1$ factors), we see that the right-hand side of
(A.3) is less or equal

$\begin{split}e^{\|f_{1}\|_{C^{k}}}\|g\|_{C^{k}}&\leq
e^{\|f_{L}\|_{C^{k}}}M_{k}\,e^{k(\|f_{L-1}\|_{C^{k}}+\cdots+\|f_{1}\|_{C^{k}})}(\|f_{L-1}\|_{C^{k}}+\cdots+\|f_{1}\|_{C^{k}})\\\
&\leq
M_{k}\,e^{k\sum_{j=1}^{L}\|f_{j}\|_{C^{k}}}(\|f_{L-1}\|_{C^{k}}+\cdots+\|f_{1}\|_{C^{k}}).\end{split}$

With
$h:=\operatorname{id}_{\Omega}+g=(\operatorname{id}_{\Omega}+f_{L-1})\circ\cdots\circ(\operatorname{id}_{\Omega}+f_{1})$,
the norm of a summand in (A.2) is bounded by

$\|f_{L}\|_{C^{k}}\|h^{(|I_{1}|)}\|_{\infty,\operatorname{op}}\cdots\|h^{(|I_{j}|)}\|_{\infty,\operatorname{op}}.$

For $a\in\\{1,\ldots,j\\}$, possible cases are that $|I_{a}|=1$; then

$\displaystyle\|h^{(|I_{a}|)}\|_{\infty,\operatorname{op}}$ $\displaystyle=$
$\displaystyle\|((\operatorname{id}_{\Omega}+f_{L-1})\circ\cdots\circ(\operatorname{id}_{\Omega}+f_{1}))^{\prime}\|_{\infty,\operatorname{op}}$
$\displaystyle\leq$
$\displaystyle(1+\|f_{L-1}\|_{\infty,\operatorname{op}})\cdots(1+\|f_{1}\|_{\infty,\operatorname{op}})$
$\displaystyle\leq$
$\displaystyle(1+\|f_{L-1}\|_{C^{k}})\cdots(1+\|f_{1}\|_{C^{k}})\leq
e^{|I_{a}|(\|f_{L-1}\|_{C^{k}}+\cdots+\|f_{1}\|_{C^{k}})}$ $\displaystyle\leq$
$\displaystyle
M_{|I_{a}|}e^{|I_{a}|(\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}})}.$

If $|I_{a}|\geq 2$, then

$\displaystyle\|h^{(|I_{a}|)}\|_{\infty,\operatorname{op}}$ $\displaystyle=$
$\displaystyle\|g^{(|I_{a}|)}\|_{\infty,\operatorname{op}}\leq\|g\|_{C^{|I_{a}|}}$
$\displaystyle\leq$ $\displaystyle
M_{|I_{a}|}\,e^{|I_{a}|(\|f_{L-1}\|_{C^{k}}+\cdots+\|f_{1}\|_{C^{k}})}(\|f_{L-1}\|_{C^{k}}+\cdots+\|f_{1}\|_{C^{k}})$
$\displaystyle\leq$ $\displaystyle
M_{|I_{a}|}\,e^{|I_{a}|(\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}})},$

using (A.4) for the final estimate. As

$e^{|I_{1}|(\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}})}\cdots
e^{|I_{j}|(\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}})}=e^{k(\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}})},$

we get

$\displaystyle\|(f_{L}*\cdots*f_{1})^{(k)}\|_{\infty,\operatorname{op}}$
$\displaystyle\leq$ $\displaystyle\left(\sum_{j=2}^{k}\sum_{P\in
P_{k,j}}M_{|I_{1}|}\cdots
M_{|I_{j}|}\right)\|f_{1}\|_{C^{k}}e^{k(\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}})}$
$\displaystyle\;+\,M_{k}\,e^{k(\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}})}(\|f_{2}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}})$
$\displaystyle\leq$ $\displaystyle
M_{k}\,e^{k\sum_{j=1}^{L}\|f_{j}\|_{C^{k}}}(\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}}),$

whence $\|f_{L}*\cdots*f_{1}\|_{C^{k}}\leq
M_{k}\,e^{k\sum_{j=1}^{L}\|f_{j}\|_{C^{k}}}(\|f_{1}\|_{C^{k}}+\cdots+\|f_{L}\|_{C^{k}})$.
This completes the proof. ∎

######  Remark A.1.

Proposition 4.1 was inspired by a recent result from the theory of evolution
equations on infinite-dimensional Lie groups (regularity theory). Consider a
Lie group $G$ modeled on a locally convex space $E$, with neutral element $e$.
Let $TG$ be the tangent bundle, $\mathfrak{g}:=T_{e}G$ and $G\times TG\to TG$,
$(g,v)\mapsto g.v:=TL_{g}(v)$ be the natural left action of $G$ on $TG$, which
uses the tangent map $TL_{g}$ of the left translation $L_{g}\colon G\to G$,
$x\mapsto gx$. The Lie group $G$ is called _$C^{0}$ -regular_ if, for each
continuous path $\gamma\colon[0,1]\to\mathfrak{g}$, there exists a $C^{1}$-map
$\eta\colon[0,1]\to G$ with $\eta(0)=e$ and

$\eta^{\prime}(t)=\eta(t).\gamma(t)\;\;\mbox{for all $\,t\in[0,1]$;}$

and moreover the time $1$-map $C([0,1],\mathfrak{g})\to G$,
$\gamma\mapsto\eta(1)$ is smooth (endowing the domain with the topology of
uniform convergence); see [14, 29]. As shown by Hanusch [19], every
$C^{0}$-regular Lie group is _locally $\mu$-convex_ in the sense of [15].
Thus, for each open identity neighborhood $U\subseteq G$ and
$C^{\infty}$-diffeomorphism $\phi\colon U\to V\subseteq E$ with $\phi(e)=0$,
and each continuous seminorm $q\colon E\to[0,\infty[$, there exists a
continuous seminorm $p\colon E\to[0,\infty[$ such that

$g_{L}\cdots g_{2}g_{1}\in U$

and

$q(\phi(g_{L}\cdots g_{1}))\leq\,p(\phi(g_{L}))+\cdots+p(\phi(g_{1}))$

for all $L\in\mathbb{N}$ and $g_{1},\ldots,g_{L}\in U$ such that
$p(\phi(g_{L}))+\cdots+p(\phi(g_{1}))\leq 1$. As shown in [15], the Lie group
$\operatorname{Diff}_{\partial}(\Omega)$ is $C^{0}$-regular. Applying the
preceding fact to the $C^{\infty}$-diffeomorphism
$\kappa\colon\operatorname{id}_{\Omega}+\mathcal{V}\to\mathcal{V}$, $g\mapsto
g-\operatorname{id}_{\Omega}$ and the seminorm $q:=\|\cdot\|_{C^{k}}$ with a
fixed $k\in\mathbb{N}_{0}$, we infer the existence of a continuous seminorm
$p$ on $C^{\infty}_{\partial}(\Omega,\operatorname{\mathbb{R}}^{d})$ such that

$\|((\operatorname{id}_{\Omega}+f_{L})\circ\cdots\circ(\operatorname{id}_{\Omega}+f_{1}))-\operatorname{id}_{\Omega}\|_{C^{k}}\leq
p(f_{1})+\cdots+p(f_{L})$

for all $L\in\mathbb{N}$ and $f_{1},\ldots,f_{L}\in\mathcal{V}$ with
$p(f_{1})+\cdots+p(f_{L})\leq 1$. We may assume that $p$ is a multiple
$r\|\cdot\|_{C^{\ell}}$ of $\|\cdot\|_{C^{\ell}}$ for some $\ell\geq k$, as
these seminorms define the topology on
$C^{\infty}_{\partial}(\Omega,\operatorname{\mathbb{R}}^{d})$. The point of
Proposition 4.1 is that we can always choose $\ell=k$, and that it provides an
explicit choice for the constant $r>0$.
Previously, $\mu$-regularity has been discussed in quantitative form only for
Banach-Lie groups, where estimates for the Baker-Campbell-Hausdorff
multiplication $*$ on the Banach-Lie algebra $\mathfrak{g}$ enable norms
$\|x_{1}*\cdots*x_{m}\|$ of products of elements
$x_{1},\ldots,x_{m}\in\mathfrak{g}$ to be estimated (see [14]).

##  Appendix B Practical upper bounds for $C^{k}$-norm

This appendix describes the methods that were used for the values in figure 2.
The goal is to estimate how the practically achieved $C^{k}$-norms compare to
the theoretical results of proposition 4.1. For simplicity, we stick to the
one-dimensional case. Since both sides of the inequality 4.1 relies on the
chosen vector fields $f_{\ell}$, we start by reorganizing the inequality
according to

$\frac{\|((\operatorname{id}_{\Omega}+f_{L})\circ\cdots\circ(\operatorname{id}_{\Omega}+f_{1}))-\operatorname{id}_{\Omega}\|_{C^{k}}}{\left(\sum_{\ell}^{L}\|f_{\ell}\|_{C^{k}}\right)e^{k\sum_{\ell}^{L}\|f_{\ell}\|_{C^{k}}}}\leq
M_{k}.$ (B.1)

##### Estimating $C^{k}$-norms

Let $X=\\{x_{i}\\}_{i=1}^{N}\in\Omega$ be a collection of points tightly
distributed throughout the domain $\Omega$. Then we estimate the $C^{k}$-norm
of a function $f\in\Omega$ by
$\|f\|_{C^{k}}=\|D^{k}f\|_{\text{op},\infty}=\max_{x\in\Omega}\|D^{k}f(x)\|_{\text{op}}\approx\max_{\mathbf{x}_{i}\in
X}\|D^{k}f(\mathbf{x_{i}})\|_{\text{op}}$ In the one-dimensional case, this
reduces to

$\max_{x_{i}\in X}\left|\frac{\partial^{k}f}{\partial x^{k}}(x_{i})\right|$
(B.2)

where $X$ is typically chosen as an equidistant grid of points
$X=\left\\{i/N\;\middle|\;i=0,...,N\ \right\\}$.

##### Choosing vector fields

Consider vector fields of the form

$f_{\ell}(x)=\sum_{j=1}^{M}w_{\ell,j}\varphi_{j}(x),$ (B.3)

for basis functions as given in (5.3) and some weight vector
$\mathbf{w}_{\ell}\in\operatorname{\mathbb{R}}^{M}$. Since different weight
vectors might yield different results for the relative norm estimates, we use
two different strategies for initializing the weight vectors:

  1. 1.

Deterministic vector of ones: $w_{\ell,j}=1,\,j=1,...,M,\,\ell=1,...,L$, and

  2. 2.

Normal random sampling: $w_{\ell,j}\sim\mathcal{N}(\mu=0,\sigma^{2}=1)$.

The values of (B.2) presented in figure 2 are chosen as the largest attained
value from strategy one and 500 runs of strategy 2.

##### Ensuring the assumptions hold

The weight vectors cannot be chosen freely, but are assumed to adhere to a
unit norm-sum constraint. To ensure that this condition holds, we will
typically start out with an initial set of weight vectors, which will be
scaled onto the feasible set. Let $W=\\{\mathbf{w}_{\ell}\\}_{\ell=1}^{L}$ be
a collection of weight vectors, from which we create a set of vector fields
$f_{\ell},\,\ell=1,...,L$. Moreover, assume that for this set of vector fields

$\sum_{\ell=1}^{L}\|f_{\ell}\|_{C^{k}}>1.$ (B.4)

The scaling will be performed in a way that is similar to the methods of
section 5.1; Find a scalar $\alpha\in(0,1]$ such that
$\alpha\sum_{\ell=1}^{L}\|f_{\ell}\|_{C^{k}}\leq 1$, i.e.

$\alpha=\min\left\\{1,\frac{1}{\sum_{\ell=1}^{L}\|f_{\ell}\|_{C^{k}}}\right\\}.$

and then choose $\tilde{W}=\\{\alpha\mathbf{w}_{\ell}\\}_{\ell=1}^{L}$.

## References

  * [1] Andrei Agrachev and Marco Caponigro  “Controllability on the Group of Diffeomorphisms”  In _Annales de l’Institut Henri Poincaré C, Analyse non linéaire_ 26.6 European Mathematical Society - EMS - Publishing House GmbH, 2009, pp. 2503–2509  DOI: 10.1016/j.anihpc.2009.07.003
  * [2] Andrei Agrachev and Andrey Sarychev  “Control on the Manifolds of Mappings with a View to the Deep Learning”  In _Journal of Dynamical and Control Systems_ Springer ScienceBusiness Media LLC, 2021  DOI: 10.1007/s10883-021-09561-2
  * [3] Martin Bauer, Martins Bruveris, Stephen Marsland and Peter W. Michor  “Constructing Reparameterization Invariant Metrics on Spaces of Plane Curves”  In _Differential Geometry and its Applications_ 34 Elsevier BV, 2014, pp. 139–165  DOI: 10.1016/j.difgeo.2014.04.008
  * [4] Jens Behrmann et al.  “Invertible Residual Networks”  In _Proceedings of the 36 th International Conference on Machine Learning_ 97, Proceedings of Machine Learning Research  PMLR, 2019, pp. 573–582  URL: https://proceedings.mlr.press/v97/behrmann19a.html
  * [5] Jens Behrmann et al.  “Understanding and Mitigating Exploding Inverses in Invertible Neural Networks”  In _Proceedings of The 24 th International Conference on Artificial Intelligence and Statistics_ 130, Proceedings of Machine Learning Research  PMLR, 2021, pp. 1792–1800  URL: https://proceedings.mlr.press/v130/behrmann21a.html
  * [6] Martin Benning et al.  “Deep learning as optimal control problems: models and numerical methods”  In _Journal of Computational Dynamics_ 6.2, 2019, pp. 171–198  DOI: 10.3934/jcd.2019009
  * [7] Elena Celledoni, Markus Eslitzbichler and Alexander Schmeding  “Shape Analysis on Lie Groups with Applications in Computer Animation”  In _Journal of Geometric Mechanics_ 8.3 American Institute of Mathematical Sciences (AIMS), 2016, pp. 273–304  DOI: 10.3934/jgm.2016008
  * [8] Elena Celledoni, Sølve Eidnes, Markus Eslitzbichler and Alexander Schmeding  “Shape Analysis on Lie Groups and Homogeneous Spaces”  In _Lecture Notes in Computer Science_ Springer International Publishing, 2017, pp. 49–56  DOI: 10.1007/978-3-319-68445-1˙6
  * [9] Elena Celledoni et al.  “Structure-preserving Deep Learning”  In _European Journal of Applied Mathematics_ 32.5 Cambridge University Press (CUP), 2021, pp. 888–936  DOI: 10.1017/s0956792521000139
  * [10] Daniel E. Clark and Jeremie Houssineau  “Faa Di Bruno’s Formula for Gateaux Differentials and Interacting Stochastic Population Processes”  arXiv, 2012  DOI: 10.48550/ARXIV.1202.0264
  * [11] George Cybenko  “Approximation by Superpositions of a Sigmoidal Function”  In _Mathematics of Control, Signals, and Systems_ 2.4 Springer ScienceBusiness Media LLC, 1989, pp. 303–314  DOI: 10.1007/bf02551274
  * [12] Weinan E  “A proposal on machine learning via dynamical systems”  In _Communications in Mathematics and Statistics_ 5.1 Springer, 2017, pp. 1–11
  * [13] Helge Glöckner  “Diffeomorphism Groups of Convex Polytopes”  arXiv, 2022  DOI: 10.48550/ARXIV.2203.09285
  * [14] Helge Glöckner  “Regularity Properties of Infinite-Dimensional Lie Groups, and Semiregularity”  arXiv, 2012  DOI: 10.48550/ARXIV.1208.0715
  * [15] Helge Glöckner and Karl-Hermann Neeb  “Diffeomorphism Groups of Compact Convex Sets”  In _Indagationes Mathematicae_ 28.4 Elsevier BV, 2017, pp. 760–783  DOI: 10.1016/j.indag.2017.04.004
  * [16] Helge Glöckner and Karl-Hermann Neeb  “Infinite-Dimensional Lie Groups” Book in preparation., 2022
  * [17] Eldad Haber and Lars Ruthotto  “Stable architectures for deep neural networks”  In _Inverse Problems_ 34.1 IOP Publishing, 2017, pp. 014004  DOI: 10.1088/1361-6420/aa9a90
  * [18] Richard S. Hamilton  “The inverse function theorem of Nash and Moser”  In _Bulletin of the American Mathematical Society_ 7.1 American Mathematical Society (AMS), 1982, pp. 65–222  DOI: 10.1090/s0273-0979-1982-15004-2
  * [19] Maximilian Hanusch  “Regularity of Lie groups” To appear in Commun. Anal. Geom., cf. arXiv:1711.03508  In _Commun. Anal. Geom._ , 2017
  * [20] Kaiming He, Xiangyu Zhang, Shaoqing Ren and Jian Sun  “Deep Residual Learning for Image Recognition”  In _IEEE Conference on Computer Vision and Pattern Recognition._ , 2016, pp. 770–778
  * [21] Ian H. Jermyn, Sebastian Kurtek, Eric Klassen and Anuj Srivastava  “Elastic Shape Matching of Parameterized Surfaces Using Square Root Normal Fields”  In _Computer Vision – ECCV 2012_ Springer Berlin Heidelberg, 2012, pp. 804–817  DOI: 10.1007/978-3-642-33715-4˙58
  * [22] Robert Keys  “Cubic Convolution Interpolation for Digital Image Processing”  In _IEEE Transactions on Acoustics, Speech, and Signal Processing_ 29.6 Institute of ElectricalElectronics Engineers (IEEE), 1981, pp. 1153–1160  DOI: 10.1109/tassp.1981.1163711
  * [23] Eric Klassen and Peter W. Michor  “Closed Surfaces with Different Shapes That Are Indistinguishable by the Srnf”  arXiv, 2019  DOI: 10.48550/ARXIV.1910.10804
  * [24] 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_ IEEE, 2010  DOI: 10.1109/cvpr.2010.5539778
  * [25] Sebastian Kurtek et al.  “Elastic Geodesic Paths in Shape Space of Parameterized Surfaces”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 34.9 Institute of ElectricalElectronics Engineers (IEEE), 2012, pp. 1717–1730  DOI: 10.1109/tpami.2011.233
  * [26] Meena Mani, Sebastian Kurtek, Christian Barillot and Anuj Srivastava  “A Comprehensive Riemannian Framework for the Analysis of White Matter Fiber Tracts”  In _2010 IEEE International Symposium on Biomedical Imaging: From Nano to Macro_ IEEE, 2010  DOI: 10.1109/isbi.2010.5490185
  * [27] Peter W. Michor  “Manifolds of Differentiable Mappings”  Shiva Pub., 1980, pp. 158
  * [28] John Milnor  “Remarks on infinite-dimensional Lie groups”  In _Relativity, groups and topology. 2_ , 1984
  * [29] Karl-Hermann Neeb  “Towards a Lie Theory of Locally Convex Groups”  In _Japanese Journal of Mathematics_ 1.2 Springer ScienceBusiness Media LLC, 2006, pp. 291–468  DOI: 10.1007/s11537-006-0606-y
  * [30] Jorge Nocedal and Stephen Wright  “Numerical Optimization (Springer Series in Operations Research and Financial Engineering)”  Springer, 2006, pp. 664
  * [31] Hideki Omori  “Infinite Dimensional Lie Transformations Groups”  Springer Berlin Heidelberg, 1974  DOI: 10.1007/bfb0063400
  * [32] Adam Paszke et al.  “PyTorch: An Imperative Style, High-Performance Deep Learning Library”  In _Advances in Neural Information Processing Systems 32_ Curran Associates, Inc., 2019, pp. 8024–8035  URL: http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf
  * [33] John Riordan  “The Blossoming of Schröder’s Fourth Problem”  In _Acta Mathematica_ 137.0 International Press of Boston, 1976, pp. 1–16  DOI: 10.1007/bf02392410
  * [34] Jørgen Nilsen Riseth  “Gradient-based Optimization in Shape Analysis for Reparametrization of Parametric Curves and Surfaces”, 2021  URL: https://hdl.handle.net/11250/2778386
  * [35] Ernst Schröder  “Vier combinatorische Probleme”  In _Z. Math. Phys_ 15, 1870, pp. 361–376
  * [36] Anuj Srivastava, Eric Klassen, Sarang H. Joshi and I.H. Jermyn  “Shape Analysis of Elastic Curves in Euclidean Spaces”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 33.7 Institute of ElectricalElectronics Engineers (IEEE), 2011, pp. 1415–1428  DOI: 10.1109/tpami.2010.184
  * [37] Zhe Su et al.  “Shape Analysis of Surfaces Using General Elastic Metrics”  In _Journal of Mathematical Imaging and Vision_ 62.8 Springer ScienceBusiness Media LLC, 2020, pp. 1087–1106  DOI: 10.1007/s10851-020-00959-4
  * [38] Esteban G. Tabak and Eric Vanden-Eijnden  “Density Estimation by Dual Ascent of the Log-likelihood”  In _Communications in Mathematical Sciences_ 8.1 International Press of Boston, 2010, pp. 217–233