# Planar Curve Registration using Bayesian Inversion

Andreas Bock<EMAIL_ADDRESS>Colin J. Cotter<EMAIL_ADDRESS>Robert
C. Kirby<EMAIL_ADDRESS>
###### Abstract

We study parameterisation-independent closed planar curve matching as a
Bayesian inverse problem. The motion of the curve is modelled via a curve on
the diffeomorphism group acting on the ambient space, leading to a _large
deformation diffeomorphic metric mapping_ (LDDMM) functional penalising the
kinetic energy of the deformation. We solve Hamilton’s equations for the curve
matching problem using the Wu-Xu element [S. Wu, J. Xu, Nonconforming finite
element spaces for $2m^{\text{th}}$ order partial differential equations on
$\mathbb{R}^{n}$ simplicial grids when $m=n+1$, Mathematics of Computation 88
(316) (2019) 531–551] which provides mesh-independent Lipschitz constants for
the forward motion of the curve, and solve the inverse problem for the
momentum using Bayesian inversion. Since this element is not affine-equivalent
we provide a _pullback theory_ which expedites the implementation and
efficiency of the forward map. We adopt ensemble Kalman inversion using a
negative Sobolev norm mismatch penalty to measure the discrepancy between the
target and the ensemble mean shape. We provide several numerical examples to
validate the approach.

###### keywords:

Closed curve matching , Nonconforming finite element method , Bayesian inverse
problem

###### PACS:

87.57.N

###### MSC:

65M60 , 65P10 , 65M32

††journal: Computers & Mathematics with Applications

[dtu]organization=Department of Applied Mathematics and Computer Science,
Technical University of Denmark, addressline=Richard Petersens Plads, Building
324, city=Kongens Lyngby, postcode=2800, country=Denmark
[imperial]organization=Department of Mathematics, Imperial College London,
addressline=180 Queen’s Gate, South Kensington, city=London, postcode=SW72RH,
country=United Kingdom [baylor]organization=Department of Mathematics, Baylor
University, addressline=1410 S.4th Street, Sid Richardson Science Building,
city=Waco, postcode=76706, state=Texas, country=United States of America

##  1 Introduction

Closed curve matching is a central problem in shape analysis where the goal is
to bring into alignment two closed curves in
$\text{Emb}(S^{1},\,\mathbb{R}^{d})$ called the _template_ and the _target_
[1]. For unparameterised curves, the shape space for these objects is
$Q=\text{Emb}(S^{1},\,\mathbb{R}^{d})\setminus\text{Diff}_{+}(S^{1})$ [2, 3].
This quotient space disassociates the curve from arbitrary reparameterisation
since they do not affect the range of the curves in question. This gives rise
to studying the commuting left and right actions of two Lie groups,
$G=\text{Diff}_{+}(\mathbb{R}^{2})$ and $H=\text{Diff}_{+}(S^{1})$ as in [4]:

$GQ=\text{Emb}(S^{1},G.\mathbb{R}^{2}),\qquad
HQ=\text{Emb}(H.S^{1},\mathbb{R}^{2}).$ (1)

In the context of developing algorithms for planar curve matching, these group
actions must be explicitly discretised. In this paper we our shape space with
the so-called _outer_ metric inherited by $G$ which acts on the ambient space.
This is in contrast to inner metrics intrinsically defined on the embedded
shape [5], see [2] for a comparison. To treat the parameterisation, one can
parameterise elements of $H$ using its Lie algebra and exploit its vector
space structure. In this paper we consider a mismatch penalty that eliminates
the need to treat $H$ explicitly. Instead we note that two closed curves
$c_{1}$ and $c_{2}$ are similar when the difference between the indicator
function $\mathbb{1}$ evaluated on their interiors is small. For some linear
differential operator $\mathcal{C}$ we therefore we define the mismatch, or
_misfit_ , between them as:

$\mathfrak{E}(c_{1},c_{2})=\|\mathbb{1}_{c_{1}}-\mathbb{1}_{c_{2}}\|_{\mathcal{C}}^{2},$
(2)

where
$\|f\|_{\mathcal{C}}^{2}=\langle\mathcal{C}^{-1}f,\mathcal{C}^{-1}f\rangle_{L^{2}}$
over some computational domain described later. For the outer metric we take
the LDDMM approach [6] and consider a one-parameter family of velocities
$t\mapsto u_{t}$ encoding the motion of the ambient space (and therefore the
shape) which simultaneously provides a distance measure.

We discretise the velocity field using finite elements, specifically the Wu-Xu
element [7]. This element provides a nonconforming discretisation for sixth
order operators; sixth order is necessary for the diffeomorphism to be
sufficiently smooth for the computations that we undertake. The implementation
of this element in Firedrake [8] is made possible by applying the theory of
[9] and techniques for code generation in [10]. Given certain assumptions on
the structure of our problem we can identify this entire family of velocities
with a single initial momentum defined as a function over the template. We
eliminate its evolution equation by using the analytical solution, and
restrict the initial conditions to only generate geodesics in the space of
unparameterised curves. This results in a _forward map_ , taking as input the
momentum and providing the diffeomorphism whose action maps the template to
the target curve. After obtaining a finite element discretisation of this map
we apply massively parallel and derivative-free ensemble Kalman inversion
which we use to invert the forward map for the initial momentum determining
the geodesic motion of the curve.

###  1.1 Previous work

Diffeomorphic registration has enjoyed a rich literature since the seminal
works [11, 12]. For curves specifically, [13, 14] present the first algorithms
for modelling curve matching via gradient descent methods. [6] represents
curves as measures onto which a Hilbert structure is endowed, and computations
of both the outer metric and the curves are done via radial reproducing
kernels producing $C^{\infty}$ velocities. In particular, curves were
represented as geometric currents. [15] studies such a varifold-based loss
function for elastic metrics, see also [16, 17, 18] for numerical frameworks
for $H^{2}$ metrics. [19] contains a review of methods related to elastic
curves.

In this paper we are concerned with higher-order metrics using finite
elements. While there is typically a loss of regularity incurred by these
methods, they offer more computationally efficient methods than e.g. kernel
methods. Finite elements also benefit from spatial adaptivity allowing for
local refinement e.g. close to embedded curves. Closest to our approach in
terms of discretisation are [20, 21] where a _particle-mesh_ method is
employed for curve matching where the curve was discretised into a finite set
of particles, acted on by an outer metric. However, we consider instead an
outer metric _finite element_ discretisation (as opposed to the intrinsic
metric in [5]). [22] presents an adaptive Eulerian FEM discretisation of the
velocity field for LDDMM using $C^{1}$ cubic Hermite elements and compares the
deformations generated using $C^{\infty}$ fields to assess the effect of the
loss of regularity. Smooth mesh deformations are also of interest in shape
optimisation where the aim is to transform a mesh such that some functional is
minimised. Finite element methods are also adopted here, with deformation
fields being discretised using B-splines [23], harmonic polynomials or
Lagrange finite elements depending the desired resolution or order [24]. Using
the finite element space introduced in [7] we can guarantee that the Lipschitz
norm remains bounded under mesh refinement without resorting to spline or
kernel discretisations. As mentioned, we use Firedrake [8] for all our
numerical experiments, see also [25] for an extension of this package for
shape optimisation.

Our formulation eliminates the need to integrate the momentum equation via its
analytical solution thereby improving on the typically larger cost of
Hamiltonian shooting based methods [26] compared to an LDDMM formulation [6].
We only need to solve an elliptic equation to obtain the velocity and use a
simple variational Euler scheme to evolve the diffeomorphism. Traditional
approaches in numerical shape analysis often apply a shooting procedures to
determine the initial momentum transporting the image or landmarks to the
desiderata, see e.g. [27, 28]. Bayesian approaches have been employed before
in the context of shape analysis, see e.g. [29] where function space Markov
Chain Monte Carlo is used to characterise the posterior density of momenta
generating a given shape. Similar to our approach is [30] in which ensemble
Kalman inversion [31, 32] is applied to recover the momentum for landmark
matching.

###  1.2 Organisation

Section 2 contains an introduction to diffeomorphic curve matching and the
associated Hamiltonian systems, We also discuss the application of the finite
element approach using the Wu-Xu element from [7] and the discretisation of
the velocity equation. Section 3 contains the transformation theory for the
Wu-Xu element, and Section 4 contains details of the discretisation of the
Hamiltonian equations. Next, Section 5 discusses the Bayesian inverse problem,
and Section 6 contains numerical results. Section 7 contains a summary.

##  2 Diffeomorphic registration

Let $\Omega$ be a connected convex subset of $\mathbb{R}^{d}$, $d=2$, with
polygonal boundary $\partial\Omega$. We study maps $q\in
Q=H^{1}(S^{1},\mathbb{R}^{d})$ from a template curve
$\Gamma_{0}\in\text{Emb}(S^{1},\,\Omega)$ to a target curve
$\Gamma_{1}\in\text{Emb}(S^{1},\,\Omega)$ whose motion is restricted by the
differential equation:

$\dot{q}_{t}=u_{t}\circ q_{t}\,,$ (3)

where $u_{t}$, $t\in[0,1]$ is a family of time-dependent vector fields on
$\Omega$ with some prescribed spatial smoothness. A _geodesic path_ between
two such parameterised curves $\Gamma_{0}$ and $\Gamma_{1}$ is defined as a
path minimising the associated kinetic energy in $u$:

$\frac{1}{2}\int_{0}^{1}\|u_{t}\|^{2}\,\text{d}t,$ (4)

where $\|\cdot\|$ dominates the Lipschitz norm. In fact, since $u_{t}$ is
supported on $\Omega$ it generates a curve on $\text{Diff}(\Omega)$ [1] of the
entire ambient space via:

$\dot{\varphi}_{t}=u_{t}\circ\varphi_{t},\quad\varphi_{0}=\textbf{id},$ (5)

whose motion restricted to the curve $q_{0}\circ S^{1}$ equals the $q_{t}\circ
S^{1}$ at time $t\in[0,1]$. As the kinetic energy measures distances between
two elements of $\text{Emb}(S^{1},\,\mathbb{R}^{d})$ via velocity defined over
the entire field $\Omega$, we refer to this associated distance measure as an
_outer_ metric on the shape space $\text{Emb}(S^{1},\,\mathbb{R}^{d})$.

###  2.1 Hamiltonian system

Here we take a Hamiltonian approach [33] and introduce the momentum $p_{t}\in
T^{*}Q$ occupying the linear cotangent space, which we assume has enough
regularity so that it has a Fréchet-Riesz representer in $L^{2}(S^{1})$ (also
denoted $p_{t}$, with some abuse of notation). We extremise the following the
functional:

$\displaystyle S=\int_{0}^{1}\frac{1}{2}\|u_{t}\|^{2}+\langle
p_{t},\dot{q}_{t}-u_{t}\circ q_{t}\rangle\,\text{d}t,$

where $\langle h,g\rangle=\int_{S^{1}}h\cdot g\,\text{d}{\theta}$. Taking
variations i.e. $\delta S=0$ leads to Hamilton’s equations for curve matching
for $t\in[0,1]$:

$\displaystyle\int_{0}^{1}\langle\delta p,\dot{q}_{t}-u_{t}\circ
q_{t}\rangle\,\text{d}t=0,$ $\displaystyle\forall\delta p\in L^{2}(S^{1}),$
(6a) $\displaystyle\int_{0}^{1}\langle\dot{p}_{t}-\nabla u_{t}^{\top}\circ
q_{t}p_{t},\delta q\rangle\,\text{d}t=0,$ $\displaystyle\forall\delta q\in Q,$
(6b) $\displaystyle{\frac{1}{2}}\frac{\delta\|u_{t}\|^{2}}{\delta u}-\langle
p_{t},\delta u\circ q_{t}\rangle=0.$ (6c)

where $\delta p$, $\delta u$ and $\delta q$ are space-time test functions. The
following theorem shows that we can solve (6b) analytically:

######  Theorem 1.

The solution $p_{t}$ to (6b) is at all times $t\geq 0$ given by
$p_{t}=\nabla\varphi_{t}^{-\top}\circ q_{0}p_{0}$.

###### Proof.

See A. ∎

To generate _parameterisation-independent_ geodesics as in [4] we replace the
initial condition $q_{0}$ by $q_{0}\circ\eta$, where
$\eta\in\text{Diff}_{+}(S^{1})$ in the case of planar curves is an arbitrary
reparameterisation. As a result of this quotient representation
$\text{Emb}(S^{1},\,\mathbb{R}^{d})\setminus\text{Diff}_{+}(S^{1})$ of curves
we minimise over all $\eta$ leading to the _horizontality condition_ on the
momentum. This means that the momentum $p_{0}$ has no tangential component and
can therefore be described by a one-dimensional signal,
$\tilde{p}_{0}:S^{1}\mapsto\mathbb{R}$:

$\displaystyle p_{0}=\mathbf{n}_{q_{0}}\tilde{p}_{0}$

where $\mathbf{n}_{q_{0}}:S^{1}\rightarrow\mathbb{R}^{2}$ is the outward
normal of the template. Thus, along with Theorem 1 we have the following
characterisation,

$p_{t}=\varphi_{t}^{-\top}\circ q_{0}^{-\top}\mathbf{n}_{q_{0}}\tilde{p}_{0}.$
(7)

This generates trajectories of geodesics between unparameterised curves. The
entire geodesic motion of the curve can therefore be determined by a one-
dimensional signal along the initial curve $q_{0}$. To summarise this section
we are concerned with integration of the following reduced Hamiltonian system
for $t\in[0,1]$:

$\displaystyle{\frac{1}{2}}\frac{\delta\|u_{t}\|^{2}}{\delta
u}=\langle\varphi_{t}^{-\top}\circ q_{0}\mathbf{n}_{q_{0}}\tilde{p}_{0},\delta
u\circ q_{t}\rangle,$ (8a) $\displaystyle\dot{q}_{t}=u_{t}\circ q_{t},$ (8b)

with $q_{0}$ and $\tilde{p}_{0}$ fixed and boundary conditions
$u_{t}|_{\partial\Omega}=0$ for all $t\in[0,1]$. Next we discuss a
discretisation of (8).

###  2.2 Outer metric via finite elements

From Picard-Lindelhöf analysis it is clear that the Banach space ordinary
differential equation (ODE) (8b) require a pointwise Lipschitz condition on
$u_{t}$. As such, $u_{t}$ must occupy at least ${W^{1,\infty}(\Omega)^{d}}$
when $q_{0}\in L^{\infty}(S^{1})$, see [29, Theorem 5] (see also Corollary 7
in this reference for other host spaces). Dupuis [12] establishes sufficient
conditions accomplishing the same in a Hilbertian setting. The Hilbertian
setting is better suited to finite element methods. This is in contrast with
${W^{1,\infty}(\Omega)^{d}}$ which is only a Banach space and, to the best of
the authors’ ability, is not easy to approximate numerically111[34]
approximates by means of a fixed point linearisation solutions to the
nonlinear $\infty$-harmonic equation [35].. We therefore request a norm
$\|\cdot\|$ such a way that a solution to (8a) ensures that this condition is
met, which in turn implies global existence and uniqueness of (8b) by the
references above. For $d=2,\,3$, $H_{0}^{3}(\Omega)$ is contained in
$\mathsf{C}^{1}(\bar{\Omega})$ and so is Lipschitz on the interior [36,
Theorem 2.5.1]. As such, we want to describe a discretisation of (8a) ensuring
a type of $H^{3}$ regularity as the follow theorem shows.

######  Theorem 2.

Let $O$ be a convex bounded Lipschitz domain in $\mathbb{R}^{d}$ with
polygonal boundary and $O_{h}$ a shape-regular, quasi-uniform triangulation
thereof [37] for some mesh size $h>0$. Suppose further that $u$ is continuous
on $\bar{O}$, $u|_{K}\in H^{3}(K)^{d}$ for $K\in O_{h}$ and that there exists
an operator $B$ inducing the norm $\|u\|_{B}^{2}=\sum_{K\in
O_{h}}\|u\|_{B(K)}^{2}$, where we define $\|u\|_{B(K)}^{2}=\int_{K}Bu\cdot
u\,\text{d}x$ such that $\|u\|_{H^{3}(K)^{d}}\lesssim\|u\|_{B(K)}$. Then $u\in
W^{1,\infty}(O)^{d}$.

###### Proof.

The embedding theorem for homogeneous Sobolev spaces (i.e. with zero traces)
into the space $\mathsf{C}^{j}(\bar{O})$ are well-known. However, since the
trace $\gamma_{K}u$ of $u$ on $\partial K$, $K\in O_{h}$ may not be zero. By
[38, Theorem 4.12], $H^{3}(K)\hookrightarrow\mathsf{C}_{B}^{1}(K)$, where:

$\displaystyle\mathsf{C}_{B}^{1}(K)=\\{u\in\mathsf{C}^{1}(K)\;|\;D^{\boldsymbol{\alpha}}u\textnormal{
is bounded on }K,\;|{\boldsymbol{\alpha}}|\leq 1\\}.$

This means any $H^{3}(K)$ function has a continuous representative with almost
everywhere bounded first derivatives on $K$. Since
$u\in\mathsf{C}^{0}(\bar{O})$, $u$ is a continuous function with its first
derivative a.e. bounded, implying a Lipschitz condition. To summarise:

$\displaystyle\|u\|_{W^{1,\infty}(K)^{d}}^{2}$
$\displaystyle\lesssim\|u\|_{H^{3}(K)^{d}}^{2}\lesssim\|u\|_{B(K)}^{2}$

Summing over the elements $K\in\Omega$ and squaring:

$\displaystyle\|u\|_{W^{1,\infty}(O)^{d}}^{2}\lesssim\|u\|_{B}^{2}.$

where we have used that $u$ is a continuous function with essentially bounded
gradient. ∎

In light of this theorem we approximate the space of velocity fields by a
nonconforming finite element space (see e.g. [39, Section 10.3]) This way we
can guarantee the necessary Lipschitz properties of our functions without
having to impose higher-order _global_ continuity of the finite-dimensional
solution spaces.

In Section 4 we use the $H^{3}$-nonconforming finite element space presented
in [7, Section 4] in a discretisation of (8). We choose the operator
$B=(\textbf{id}-\alpha\Delta)^{2m}$ for a given positive constant $\alpha$
leading to the following bilinear form:

$a_{\Omega}(u,v)=\sum_{i=1}^{d}\int_{\Omega}\sum_{j=0}^{m}\alpha^{j}\binom{m}{j}D^{j}u^{i}\cdot
D^{j}v^{i}\,\text{d}x=\int_{\Omega}Bu\cdot v\,\text{d}x,$ (9)

where $x\cdot y$ is the Euclidean inner product, $D^{0}=\textbf{id}$, and

$D^{j}=\begin{cases}\nabla D^{j-1}&j\text{ is odd},\\\ \nabla\cdot
D^{j-1}&j\text{ is even}.\end{cases}$

##  3 A pullback theory for the Wu-Xu element

The Wu-Xu element provides an opportunity to tackle this problem in a
(nonconforming) $H^{3}$ setting, but it presents challenges for
implementation. Although we can construct its basis on a reference element,
say, using the FIAT package [40], the Wu-Xu elements do not form an affine
equivalent family [39] under pullback. Consequently, we apply the theory
developed in [9], which gives a generalization of techniques developed for the
$C^{1}$ conforming Argyris element [41, 42].

To fix ideas, put a reference triangle $\widehat{K}$ with vertices by
$\\{\widehat{\mathbf{v}}_{i}\\}_{i=1}^{3}$. For any nondegenerate triangle $K$
with vertices $\\{\mathbf{v}_{i}\\}_{i=1}^{3}$, we let
$F:T\rightarrow\widehat{K}$ denote the affine mapping sending each
$\mathbf{v}_{i}$ to the corresponding $\widehat{\mathbf{v}}_{i}$ and $J_{T}$
its Jacobian matrix.

$\hat{\mathbf{v}}_{1}$$\hat{\mathbf{v}}_{2}$$\hat{\mathbf{v}}_{3}$$\mathbf{v}_{1}$$\mathbf{v}_{2}$$\mathbf{v}_{3}$$F:K\rightarrow\hat{K}$$\hat{K}$$K$
Figure 1: Affine mapping to a reference cell $\hat{K}$ from a typical cell
$K$. Note that here $F$ maps from the physical cell $K$ to the reference cell
$\hat{K}$ rather than the other way around.

We adopt the ordering convention used in [43], where edge $e_{i}$ of any
triangle connects the vertices other than $i$. We take the unit tangent
$\mathbf{t}_{i}=\begin{bmatrix}t_{i}^{x}&t_{i}^{y}\end{bmatrix}^{T}$ to from
the vertex of lower number to the higher one. The normal to edge $i$ is
defined by counterclockwise rotation of the tangent, so that
$\mathbf{n}_{i}=R\mathbf{t}_{i}$, where $R=\left[\begin{smallmatrix}0&1\\\
-1&0\end{smallmatrix}\right]$. The normals, tangents, and edge midpoints for
the reference element $\widehat{K}$ will include hats:
$\widehat{\mathbf{n}}_{i}$, $\widehat{\mathbf{t}}_{i}$, and
$\widehat{\mathbf{e}}_{i}$. The pull-back of any function $\widehat{f}$
defined on $\widehat{K}$ is given by

$F^{*}(\widehat{f})=\widehat{f}\circ F,$ (10)

and the push-forward of functionals $n$ acting on functions defined over $K$
is

$F_{*}(n)=n\circ F^{*},$ (11)

so that

$F_{*}(n)(\widehat{f})=(n\circ F^{*})(\widehat{f})=n(\widehat{f}\circ F)$ (12)

Finite element implementation requires local shape functions
$\\{\psi_{i}^{K}\\}_{i=1}^{N}$ that are restrictions of the global basis to
cell $K$. These are taken dual to a set of _nodes_ or _degrees of freedom_
$\\{n_{i}^{K}\\}_{i=1}^{N}$ in the sense that

$n_{i}^{K}(\psi_{j}^{K})=\delta_{ij}.$

In practice, one typically computes the basis $\\{\hat{\psi}_{i}\\}_{i=1}^{N}$
dual to some nodes $\\{\hat{n}_{i}\\}_{i=1}^{N}$ over the reference element
$\hat{K}$. For affine equivalent families (like the Lagrange basis), the
physical basis functions are the pullbacks of reference element shape
functions, so that

$\psi_{i}^{K}=F^{*}(\widehat{\psi}_{i}).$

Equivalently, the nodes are preserved under push-forward, with

$F_{*}(n^{K}_{i})=\hat{n}_{i}.$

We may express these relations in a kind of vector-notation. If $\hat{\Psi}$
is a vector whose entries are $\widehat{\psi}_{i}$, then in the affine
equivalent case, $F^{*}(\hat{\Psi})$ contains the basis on cell $K$, and also
$F_{*}(\mathcal{N})=\widehat{\mathcal{N}}$. For non-equivalent families, these
relations fail, but we can hope to construct a matrix $M$ such that

$\Psi=MF^{*}(\hat{\Psi})$ (13)

contains the correct vector of basis functions on $T$. The matrix $M$ will
depend on the particular geometry of each cell, but if it is sparse this
amounts to a considerable savings over directly constructing the basis on each
triangle. Our theory in [9] proceeds by transforming the actions of the
functionals on the finite element space. The finite element functionals are
defined on some infinite-dimensional space (e.g. twice-continuously
differentiable functions), and we let $\pi$ denote the restriction of
functionals to the finite-element space and $\hat{\pi}$ the corresponding
restriction on the reference element. Then, we look for a matrix $V$ such that

$VF_{*}(\pi\mathcal{N})=\widehat{\pi}\widehat{\mathcal{N}},$ (14)

and can prove [9, Theorem 3.1] that

$M=V^{T}.$ (15)

For any triangle $K$ and integer $k\geq 0$, we let $\mathcal{P}^{k}(K)$ denote
the space of polynomials of degree no greater than $k$ over $K$. Letting
$\lambda_{i}$ be the barycentric coordinates for $K$ (equivalently, the
Lagrange basis for $\mathcal{P}^{1}(K)$), we let
$b_{K}=\lambda_{1}\lambda_{2}\lambda_{3}$ be the standard cubic bubble
function over $K$. We also need notation for the linear functionals defining
degrees of freedom. We let $\delta_{\mathbf{x}}$ denote pointwise evaluation
of some (continuous) function:

$\delta_{\mathbf{x}}(p)=p(\mathbf{x}).$ (16)

We let $\delta_{\mathbf{x}}^{\mathbf{x}}$ denote the derivative in some
direction $\mathbf{s}$ at a point $\mathbf{x}$:

$\delta^{\mathbf{s}}_{\mathbf{x}}(p)=\mathbf{s}^{T}\nabla p(\mathbf{x})$ (17)

Repeated superscripts will indicate higher derivatives. We use block notation
will for gradients and sets of second-order derivatives, such as

$\nabla_{\mathbf{x}}=\begin{bmatrix}\delta_{\mathbf{x}}^{\mathbf{x}}&\delta_{\mathbf{x}}^{\mathbf{y}}\end{bmatrix}^{T}$
(18)

for the gradient in Cartesian coordinates at a point $\mathbf{x}$, and

$\bigtriangleup_{\mathbf{x}}=\begin{bmatrix}\delta_{\mathbf{x}}^{\mathbf{x}\mathbf{x}}&\delta_{\mathbf{x}}^{\mathbf{x}\mathbf{y}}&\delta_{\mathbf{x}}^{\mathbf{y}\mathbf{y}}\end{bmatrix}^{T}$
(19)

for the unique components of the Hessian matrix. We will use superscripts in
the block notation to indicate the derivatives taken in other directions than
the Cartesian ones, such as $\nabla^{\mathbf{n}\mathbf{t}}$ containing the
derivatives with respect to a normal vector $\mathbf{n}$ and tangent vector
$\mathbf{t}$ for some part of the boundary. Similarly,
$\bigtriangleup^{\mathbf{n}\mathbf{t}}$ will contain the second partials in
each direction and the mixed partial in both directions.

The Wu-Xu elements also utilise integral moments of normal derivatives, and we
shall also need averages tangential and mixed derivatives over edges to
perform the transformations. Given any directional vector $\mathbf{s}$, we
define the moment of the derivative in the direction $\mathbf{s}$ over edge
$\mathbf{e}$ by:

$\mu^{\mathbf{s}}_{\mathbf{e}}(f)=\int_{\mathbf{e}}\mathbf{s}\cdot\nabla
f\,ds,\\\ $ (20)

Similarly, we let $\mu^{\mathbf{s}_{1}\mathbf{s}_{2}}_{\mathbf{e}}$ to denote
the functionals computing moments of second (possibly mixed) directional
derivatives over an edge. Now, we define the pair of $H^{3}$ nonconforming
triangles considered in [7]. Note that there are two spaces given: a space
compatible with sixth-order problems, and a _robust_ space that is stable for
second, fourth and sixth-order problems. We define function space
$\mathcal{W}(K)$ over some triangle $K$ by

$\mathcal{W}(K)=\mathcal{P}^{3}+b_{K}\mathcal{P}^{1},$ (21)

and the function space for the robust element will be

$\mathcal{W}_{r}(K)=\mathcal{P}^{3}+b_{K}\mathcal{P}^{1}+b_{K}^{2}\mathcal{P}^{1},$
(22)

where $\mathcal{P}^{k}$ is the standard space of polynomials of degree $k$.
Note that we have $\dim\mathcal{W}(K)=12$ and $\dim\mathcal{W}_{r}(K)=15$
since $b_{K}\in\mathcal{P}^{3}\cap b_{K}\mathcal{P}^{1}$. The degrees of
freedom for the two elements are quite similar. We can parametrise
$\mathcal{W}_{r}(K)$ by

$\mathcal{N}=\begin{bmatrix}\delta_{\mathbf{v}_{1}}&\nabla^{T}_{\mathbf{v}_{1}}&\delta_{\mathbf{v}_{2}}&\nabla^{T}_{\mathbf{v}_{2}}&\delta_{\mathbf{v}_{3}}&\nabla^{T}_{\mathbf{v}_{3}}&\mu^{\mathbf{n}_{1}\mathbf{n}_{1}}_{\mathbf{e}_{1}}&\mu^{\mathbf{n}_{2}\mathbf{n}_{2}}_{\mathbf{e}_{2}}&\mu^{\mathbf{n}_{3}\mathbf{n}_{3}}_{\mathbf{e}_{3}}\end{bmatrix}^{T}.$
(23)

That is, the degrees of freedom consist of point values and gradients at each
vertex, together with moments of the second normal derivative along edges. For
the robust element, we also use the moments of the first normal derivatives,
so that

$\mathcal{N}_{r}=\setcounter{MaxMatrixCols}{12}\begin{bmatrix}\delta_{\mathbf{v}_{1}}&\nabla^{T}_{\mathbf{v}_{1}}&\delta_{\mathbf{v}_{2}}&\nabla^{T}_{\mathbf{v}_{2}}&\delta_{\mathbf{v}_{3}}&\nabla^{T}_{\mathbf{v}_{3}}&\mu^{\mathbf{n}_{1}}_{\mathbf{e}_{1}}&\mu^{\mathbf{n}_{2}}_{\mathbf{e}_{2}}&\mu^{\mathbf{n}_{3}}_{\mathbf{e}_{3}}&\mu^{\mathbf{n}_{1}\mathbf{n}_{1}}_{\mathbf{e}_{1}}&\mu^{\mathbf{n}_{2}\mathbf{n}_{2}}_{\mathbf{e}_{2}}&\mu^{\mathbf{n}_{3}\mathbf{n}_{3}}_{\mathbf{e}_{3}}\end{bmatrix}^{T}.$
(24)

Wu and Xu actually define the degrees of freedom as average of these moments
over the relevant facets, although this does not affect unisolvence or other
essential properties. For the reference element, it will be helpful to use
their original definition. For some edge $\mathbf{e}$ of $\hat{K}$, define

$\hat{\mu}^{\hat{\mathbf{s}}}_{\hat{\mathbf{e}}}(f)=\tfrac{1}{|\hat{\mathbf{e}}|}\int_{\hat{\mathbf{e}}}\hat{\mathbf{s}}\cdot\hat{\nabla}f\,d\hat{s},$
(25)

and similarly define moments second directional derivatives over reference
element edges. The reference element nodes for $\mathcal{W}(\hat{K})$ will be
taken as

$\widehat{\mathcal{N}}=\begin{bmatrix}\delta_{\hat{\mathbf{v}}_{1}}&\hat{\nabla}^{T}_{\hat{\mathbf{v}}_{1}}&\delta_{\hat{\mathbf{v}}_{2}}&\hat{\nabla}^{T}_{\hat{\mathbf{v}}_{2}}&\delta_{\hat{\mathbf{v}}_{3}}&\hat{\nabla}^{T}_{\hat{\mathbf{v}}_{3}}&\hat{\mu}^{\hat{\mathbf{n}}_{1}\hat{\mathbf{n}}_{1}}_{\hat{\mathbf{e}}_{1}}&\hat{\mu}^{\hat{\mathbf{n}}_{2}\hat{\mathbf{n}}_{2}}_{\hat{\mathbf{e}}_{2}}&\hat{\mu}^{\hat{\mathbf{n}}_{3}\hat{\mathbf{n}}_{3}}_{\hat{\mathbf{e}}_{3}}\end{bmatrix}^{T},$
(26)

and for $\mathcal{W}_{r}(\hat{K})$ we will use

$\widehat{\mathcal{N}}=\setcounter{MaxMatrixCols}{12}\begin{bmatrix}\delta_{\hat{\mathbf{v}}_{1}}&\hat{\nabla}^{T}_{\hat{\mathbf{v}}_{1}}&\delta_{\hat{\mathbf{v}}_{2}}&\hat{\nabla}^{T}_{\hat{\mathbf{v}}_{2}}&\delta_{\hat{\mathbf{v}}_{3}}&\hat{\nabla}^{T}_{\hat{\mathbf{v}}_{3}}&\hat{\mu}^{\hat{\mathbf{n}}_{1}}_{\hat{\mathbf{e}}_{1}}&\hat{\mu}^{\hat{\mathbf{n}}_{2}}_{\hat{\mathbf{e}}_{2}}&\hat{\mu}^{\hat{\mathbf{n}}_{3}}_{\hat{\mathbf{e}}_{3}}&\hat{\mu}^{\hat{\mathbf{n}}_{1}\hat{\mathbf{n}}_{1}}_{\hat{\mathbf{e}}_{1}}&\hat{\mu}^{\hat{\mathbf{n}}_{2}\hat{\mathbf{n}}_{2}}_{\hat{\mathbf{e}}_{2}}&\hat{\mu}^{\hat{\mathbf{n}}_{3}\hat{\mathbf{n}}_{3}}_{\hat{\mathbf{e}}_{3}}\end{bmatrix}^{T}$
(27)

Note that this redefinition has no effect in the case of an equilateral
reference triangle with unit edge length. For the more common case of a right
isosceles reference triangle, however, this will eliminate the need for logic
indicating to which reference element edges the edges of each triangle
correspond.

(a) $\mathcal{W}(K)$

(b) $\mathcal{W}_{r}(K)$

Figure 2: Degrees of freedom for the Wu-Xu (left) and robust Wu-Xu (right)
elements. Point values are given by dots, gradients by cirles, while averages
of first and second normal derivatives are given by thin and thick arrows,
respectively.

The derivative degrees of freedom in both Wu-Xu elements are not preserved
under push-forward, and since we have only normal derivatives on the edges, we
cannot immediately obtain the correct nodes by taking linear combinations.
Consequently, we must develop a _compatible nodal completion_ [9, Definition
3.4]. For the Wu-Xu elements, this contains all the original degrees of
freedom plus the integrals of tangential and mixed normal/tangential
derivatives. Such a completion is shown for the standard Wu-Xu element in
Figure 3(a). A completion for the robust element includes the first normal
moments and tangential moments as well, as showin in Figure 3(b).

(a) Wu-Xu

(b) Robust Wu-Xu

Figure 3: Compatible nodal completions for the Wu-Xu and robust Wu-Xu elements

We define

$\mathcal{M}_{1,i}=\begin{bmatrix}\mu^{\mathbf{n}_{i}}_{\mathbf{e}_{i}}&\mu^{\mathbf{t}_{i}}_{\mathbf{e}_{i}}\end{bmatrix}^{T}$
(28)

to be the vector of the moments of the normal and tangential derivatives on a
particular edge. We also let $\widehat{\mathcal{M}_{1,i}}$ contain the
corresponding reference element nodes. We only need $\mathcal{M}_{1,i}$ and
$\widehat{\mathcal{M}_{1,i}}$ for the robust element. Both elements require

$\mathcal{M}_{2,i}=\begin{bmatrix}\mu^{\mathbf{n}_{i}\mathbf{n}_{i}}_{\mathbf{e}_{i}}&\mu^{\mathbf{t}_{i}\mathbf{t}_{i}}_{\mathbf{e}_{i}}&\mu^{\mathbf{n}_{i}\mathbf{t}_{i}}_{\mathbf{e}_{i}}\end{bmatrix}^{T}$
(29)

containing the unique second derivative moments on each edge. We similarly
define $\widehat{\mathcal{M}_{2,i}}$ to contain the reference element integral
averages. The compatible nodal completion for $(K,\mathcal{P}(K),\mathcal{N})$
is

$\mathcal{N}^{C}=\begin{bmatrix}\delta_{\mathbf{v}_{1}}&\nabla^{T}_{\mathbf{v}_{1}}&\delta_{\mathbf{v}_{2}}&\nabla^{T}_{\mathbf{v}_{2}}&\delta_{\mathbf{v}_{3}}&\nabla^{T}_{\mathbf{v}_{3}}&\mathcal{M}_{2,1}^{T}&\mathcal{M}_{2,2}^{T}&\mathcal{M}_{2,3}^{T}\end{bmatrix}^{T},$
(30)

with the hatted equivalents comprising $\hat{\mathcal{N}}^{C}$ on the
reference cell. The completed set of nodes for the robust element is

${\scriptsize\mathcal{N}_{r}^{C}=\setcounter{MaxMatrixCols}{12}\begin{bmatrix}\delta_{\mathbf{v}_{1}}&\nabla^{T}_{\mathbf{v}_{1}}&\delta_{\mathbf{v}_{2}}&\nabla^{T}_{\mathbf{v}_{2}}&\delta_{\mathbf{v}_{3}}&\nabla^{T}_{\mathbf{v}_{3}}&\mathcal{M}_{1,1}^{T}&\mathcal{M}_{1,2}^{T}&\mathcal{M}_{1,3}^{T}&\mathcal{M}_{2,1}^{T}&\mathcal{M}_{2,2}^{T}&\mathcal{M}_{2,3}^{T}\end{bmatrix}^{T},}$
(31)

Now, the matrix $V$ from (14) will be obtained in factored form

$V=EV^{c}D,$ (32)

where each matrix plays a particular role. $D$ is a rectangular matrix
expressing the completed nodes in terms of the given physical nodes. $V^{c}$
is a block diagonal matrix relating the push-forward of the reference nodal
completion to the physical nodal completion, and $E$ is a Boolean matrix
selecting actual finite element nodes from the completion. For the Wu-Xu
element, $D$ is $18\times 12$, $V^{c}$ is $18\times 18$, and $E$ is $12\times
18$. For the robust element, $D$ is $24\times 15$, $V^{c}$ is $24\times 24$,
and $E$ is $15\times 24$.

Now, we define the matrix $D$, which expresses the members of
$\mathcal{N}^{C}$ as linear combinations of the members of $\mathcal{N}$.
Clearly, the rows corresponding to members of $\mathcal{N}^{C}$ also appearing
in $\mathcal{N}$ will just have a single nonzero in the appropriate column.
For the Wu-Xu element, the remaining nodes are all integrals of quantities
over edges, and we can use the Fundamental Theorem of Calculus to perform this
task. Let $\mathbf{e}$ be an edge running from vertex $\mathbf{v}_{a}$ to
$\mathbf{v}_{b}$ with unit tangent and normal $\mathbf{t}$ and $\mathbf{n}$,
respectively. We have

$\mu^{\mathbf{t}}_{\mathbf{e}}(f)=\int_{\mathbf{e}}\mathbf{t}^{T}\nabla
fds=f(\mathbf{v}_{b})-f(\mathbf{v}_{a})=\delta_{\mathbf{v}_{b}}(f)-\delta_{\mathbf{v}_{a}}(f).$
(33)

In a similar way, the moments of the second tangential and mixed derivatives
on $\mathbf{e}$ can be expressed as differences between components of the
gradients at endpoints by:

$\begin{split}\mu^{\mathbf{n}\mathbf{t}}_{\mathbf{e}}(f)&=\mathbf{n}^{T}\left(\nabla_{\mathbf{v}_{b}}f-\nabla_{\mathbf{v}_{a}}f\right),\\\
\mu^{\mathbf{t}\mathbf{t}}_{\mathbf{e}}(f)&=\mathbf{t}^{T}\left(\nabla_{\mathbf{v}_{b}}f-\nabla_{\mathbf{v}_{a}}f\right),\end{split}$
(34)

and we have that $\mathcal{N}^{C}=D\mathcal{N}$, or

${\tiny\begin{bmatrix}\delta_{\mathbf{v}_{1}}\\\
\delta^{\mathbf{x}}_{\mathbf{v}_{1}}\\\
\delta^{\mathbf{y}}_{\mathbf{v}_{1}}\\\ \delta_{\mathbf{v}_{2}}\\\
\delta^{\mathbf{x}}_{\mathbf{v}_{2}}\\\
\delta^{\mathbf{y}}_{\mathbf{v}_{2}}\\\ \delta_{\mathbf{v}_{3}}\\\
\delta^{\mathbf{x}}_{\mathbf{v}_{3}}\\\
\delta^{\mathbf{y}}_{\mathbf{v}_{3}}\\\
\hline\cr\mu^{\mathbf{n}_{1}\mathbf{n}_{1}}_{\mathbf{e}_{1}}\\\
\mu^{\mathbf{n}_{1}\mathbf{t}_{1}}_{\mathbf{e}_{1}}\\\
\mu^{\mathbf{t}_{1}\mathbf{t}_{1}}_{\mathbf{e}_{1}}\\\
\mu^{\mathbf{n}_{2}\mathbf{n}_{2}}_{\mathbf{e}_{2}}\\\
\mu^{\mathbf{n}_{2}\mathbf{t}_{2}}_{\mathbf{e}_{2}}\\\
\mu^{\mathbf{t}_{2}\mathbf{t}_{2}}_{\mathbf{e}_{2}}\\\
\mu^{\mathbf{n}_{3}\mathbf{n}_{3}}_{\mathbf{e}_{3}}\\\
\mu^{\mathbf{n}_{3}\mathbf{t}_{3}}_{\mathbf{e}_{3}}\\\
\mu^{\mathbf{t}_{3}\mathbf{t}_{3}}_{\mathbf{e}_{3}}\end{bmatrix}=\setcounter{MaxMatrixCols}{12}\begin{bmatrix}1&0&0&0&0&0&0&0&0&0&0&0\\\
0&1&0&0&0&0&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0&0&0&0&0\\\
0&0&0&1&0&0&0&0&0&0&0&0\\\ 0&0&0&0&1&0&0&0&0&0&0&0\\\
0&0&0&0&0&1&0&0&0&0&0&0\\\ 0&0&0&0&0&0&1&0&0&0&0&0\\\
0&0&0&0&0&0&0&1&0&0&0&0\\\ 0&0&0&0&0&0&0&0&1&0&0&0\\\ \hline\cr
0&0&0&0&0&0&0&0&0&1&0&0\\\
0&0&0&0&-n_{1,x}&-n_{1,y}&0&n_{1,x}&n_{1,y}&0&0&0\\\
0&0&0&0&-t_{1,x}&-t_{1,y}&0&t_{1,x}&t_{1,y}&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&1&0\\\
0&-n_{2,x}&-n_{2,y}&0&0&0&0&n_{2,x}&n_{2,y}&0&0&0\\\
0&-t_{2,x}&-t_{2,y}&0&0&0&0&t_{2,x}&t_{2,y}&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&1\\\
0&-n_{3,x}&-n_{3,y}&0&n_{3,x}&n_{3,y}&0&0&0&0&0&0\\\
0&-t_{3,x}&-t_{3,y}&0&t_{3,x}&t_{3,y}&0&0&0&0&0&0\\\
\end{bmatrix}\begin{bmatrix}\delta_{\mathbf{v}_{1}}\\\
\delta^{\mathbf{x}}_{\mathbf{v}_{1}}\\\
\delta^{\mathbf{y}}_{\mathbf{v}_{1}}\\\ \delta_{\mathbf{v}_{2}}\\\
\delta^{\mathbf{x}}_{\mathbf{v}_{2}}\\\
\delta^{\mathbf{y}}_{\mathbf{v}_{2}}\\\ \delta_{\mathbf{v}_{3}}\\\
\delta^{\mathbf{x}}_{\mathbf{v}_{3}}\\\
\delta^{\mathbf{y}}_{\mathbf{v}_{3}}\\\
\hline\cr\mu^{\mathbf{n}_{1}\mathbf{n}_{1}}_{\mathbf{e}_{1}}\\\
\mu^{\mathbf{n}_{2}\mathbf{n}_{2}}_{\mathbf{e}_{2}}\\\
\mu^{\mathbf{n}_{3}\mathbf{n}_{3}}_{\mathbf{e}_{3}}\\\ \end{bmatrix}.}$ (35)

The matrix $V^{C}$ is obtained by relating the push-forwards of the nodal
completion to their reference counterparts.

We can convert between the Cartesian and other orthogonal coordinate systems
(e.g. normal/tangential) representations as follows. Given a pair of
orthogonal unit vectors $\mathbf{n}$ and $\mathbf{t}$, we can define an
orthogonal matrix $G$ by:

$G=\begin{bmatrix}\mathbf{n}&\mathbf{t}\end{bmatrix}^{T}.$ (36)

In particular, we will use $G_{i}$ to have the normal and tangential vectors
to edge $i$ of triangle $K$ and $\widehat{G}_{i}$ those for triangle
$\widehat{K}$. The multivariate chain rule readily shows that

$\nabla_{x}=G^{T}\nabla^{\mathbf{n}\mathbf{t}}_{\mathbf{x}}.$ (37)

Similarly, letting $\mathbf{n}=\begin{bmatrix}n_{x}&n_{y}\end{bmatrix}^{T}$
and $\mathbf{t}=\begin{bmatrix}t_{x}&t_{y}\end{bmatrix}^{T}$, we define the
matrix $\Gamma$ by

$\Gamma=\begin{bmatrix}n_{x}^{2}&2n_{x}t_{x}&t_{x}^{2}\\\
n_{x}n_{y}&n_{x}t_{y}+n_{y}t_{x}&t_{x}t_{y}\\\
n_{y}^{2}&2n_{y}t_{y}&t_{y}^{2}\end{bmatrix},$ (38)

and the chain rule gives

$\bigtriangleup_{x}=\Gamma\bigtriangleup_{x}^{\mathbf{n}\mathbf{t}}.$ (39)

Although $G$ is an orthogonal matrix, $\Gamma$ is not. A similar calculation
also shows gives that:
$\bigtriangleup^{\mathbf{n}\mathbf{t}}_{\mathbf{x}}=\Gamma^{-1}\bigtriangleup_{\mathbf{x}}$,
where

$\Gamma^{-1}=\begin{bmatrix}n_{x}^{2}&2n_{x}n_{y}&n_{y}^{2}\\\
n_{x}t_{x}&n_{x}t_{y}+n_{y}t_{x}&n_{y}t_{y}\\\
t_{x}^{2}&2t_{x}t_{y}&t_{y}^{2}\end{bmatrix}.$ (40)

We will also need to transform derivatives under pull-back. Using the chain
rule,

$\nabla(\hat{\psi}\circ F)=J^{T}\hat{\nabla}\hat{\psi}\circ F.$ (41)

Combining this with (37) lets us relate the normal and tangential derivatives
in physical space to the normal and tangential derivatives in reference space.

$\nabla_{\mathbf{x}}^{\mathbf{n}\mathbf{t}}=GJ^{T}\widehat{G}^{T}\hat{\nabla}^{\hat{\mathbf{n}}\hat{\mathbf{t}}}_{\hat{\mathbf{x}}}.$
(42)

We can perform a similar calculation for second derivatives. With the entries
of the Jacobian matrix as:

$J=\begin{bmatrix}\tfrac{\partial x}{\partial\hat{x}}&\tfrac{\partial
x}{\partial\hat{y}}\\\ \tfrac{\partial y}{\partial\hat{x}}&\tfrac{\partial
y}{\partial\hat{y}}\end{bmatrix},$ (43)

we define the matrix

$\Theta=\begin{bmatrix}\left(\tfrac{\partial\hat{x}}{\partial
x}\right)^{2}&2\tfrac{\partial\hat{x}}{\partial
x}\tfrac{\partial\hat{y}}{\partial x}&\left(\tfrac{\partial\hat{y}}{\partial
x}\right)^{2}\\\ \tfrac{\partial\hat{x}}{\partial
y}\tfrac{\partial\hat{x}}{\partial x}&\tfrac{\partial{\hat{x}}}{\partial
y}\tfrac{\partial\hat{y}}{\partial x}+\tfrac{\partial\hat{x}}{\partial
x}\tfrac{\partial\hat{y}}{\partial y}&\tfrac{\partial\hat{y}}{\partial
x}\tfrac{\partial\hat{y}}{\partial y}\\\
\left(\tfrac{\partial\hat{x}}{\partial
y}\right)^{2}&2\tfrac{\partial\hat{x}}{\partial
y}\tfrac{\partial\hat{y}}{\partial y}&\left(\tfrac{\partial\hat{y}}{\partial
y}\right)^{2}\end{bmatrix},$ (44)

so that for $\mathbf{x}=F(\hat{\mathbf{x}})$,

$\bigtriangleup_{\mathbf{x}}=\Theta\hat{\bigtriangleup}_{\hat{\mathbf{x}}}.$
(45)

The inverse of $\Theta$ follows by reversing the roles of reference and
physical variables:

$\Theta^{-1}=\begin{bmatrix}\left(\tfrac{\partial
x}{\partial\hat{x}}\right)^{2}&2\tfrac{\partial
x}{\partial\hat{x}}\tfrac{\partial y}{\partial\hat{x}}&\left(\tfrac{\partial
y}{\partial\hat{x}}\right)^{2}\\\ \tfrac{\partial
x}{\partial\hat{y}}\tfrac{\partial
x}{\partial\hat{x}}&\tfrac{\partial{x}}{\partial\hat{y}}\tfrac{\partial
y}{\partial\hat{x}}+\tfrac{\partial x}{\partial\hat{x}}\tfrac{\partial
y}{\partial\hat{y}}&\tfrac{\partial y}{\partial\hat{x}}\tfrac{\partial
y}{\partial\hat{y}}\\\ \left(\tfrac{\partial
x}{\partial\hat{y}}\right)^{2}&2\tfrac{\partial
x}{\partial\hat{y}}\tfrac{\partial y}{\partial\hat{y}}&\left(\tfrac{\partial
y}{\partial\hat{y}}\right)^{2}\end{bmatrix}$ (46)

We can also relate the second-order derivatives in normal/tangential
coordinates under pullback by

$\bigtriangleup^{\mathbf{n}\mathbf{t}}_{\mathbf{x}}=\Gamma\Theta\hat{\Gamma}^{-1}\bigtriangleup_{\hat{\mathbf{x}}}^{\hat{\mathbf{n}}\hat{\mathbf{t}}}.$
(47)

From here, we will let $G_{i}$ and $\hat{G}_{i}$ denote the matrices
containing normal and tangent vectors to edge $\mathbf{e}_{i}$ of a generic
triangle $T$ and the reference triangle $\hat{T}$, respectively, with similar
convention for the other geometric quantities $\Gamma$ and $\Theta$. For any
vector $\mathbf{s}$, edge $\mathbf{e}$, and smooth function $f=f\circ F$, we
have

$\int_{\mathbf{e}}\mathbf{s}^{T}\nabla
fds=\int_{\mathbf{e}}\mathbf{s}^{T}\hat{\nabla}f\circ
Fds=\int_{\hat{\mathbf{e}}}\mathbf{s}^{T}\hat{\nabla}fJ_{\mathbf{e},\hat{\mathbf{e}}}d\hat{s},$
(48)

where the Jacobian $J_{\mathbf{e},\hat{\mathbf{e}}}$ is just the ratio of the
length of $\mathbf{e}$ to that of the corresponding reference element edge
$\hat{\mathbf{e}}$. Applying this to the normal and tangential moments and
using (42), we have that:

$\mathcal{M}_{1,i}=|\mathbf{e}_{i}|G_{i}J^{T}\hat{G}_{i}^{-1}\widehat{\mathcal{M}}_{1,i},$
(49)

where the factor of $|\hat{\mathbf{e}}_{i}|$ in the denominator of the
Jacobian is merged with the reference element moments to produce
$\widehat{\mathcal{M}_{1,i}}$. Hence, the slight modification of reference
element nodes avoids extra data structures or logic in identifying reference
element edge numbers. Then, we can use (47) to express each
$\mathcal{M}_{2,i}$ in terms of the reference element nodes

$\mathcal{M}_{2,i}=|\mathbf{e}|\Gamma_{i}\Theta\hat{\Gamma}_{i}^{-1}\widehat{\mathcal{M}}_{2,i}.$
(50)

We define vectors

$B^{1,i}=\tfrac{1}{|\mathbf{e}_{i}|}\hat{G}_{i}J^{-T}G_{i}^{T},\ \ \
B^{2,i}=\tfrac{1}{|\mathbf{e}_{i}|}\hat{\Gamma}_{i}\Theta^{-1}\Gamma_{i}^{-1},$
(51)

and hence $V^{C}$ is the block-diagonal matrix

$V^{C}=\begin{bmatrix}1&&&&&&&&\\\ &J^{-T}&&&&&&&\\\ &&1&&&&&&\\\
&&&J^{-T}&&&&&\\\ &&&&1&&&&\\\ &&&&&J^{-T}&&\\\ &&&&&&B^{2,1}&&\\\
&&&&&&&B^{2,2}&\\\ &&&&&&&&B^{2,3}\end{bmatrix},$ (52)

with zeros of the appropriate shapes in the off-diagonal blocks. The
extraction matrix $E$ is just the $12\times 18$ Boolean matrix selecting the
members of $N$ from $N^{C}$:

$E=\setcounter{MaxMatrixCols}{18}\begin{bmatrix}1&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&1&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&1&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&1&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&1&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&1&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&1&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&1&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&1&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&1&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&1&0&0\\\ \end{bmatrix}$ (53)

Multiplying $EV^{C}D$ out and defining

$\begin{split}\beta^{i,x}&=n^{i}_{x}B^{2,i}_{12}+t^{i}_{x}B^{2,i}_{13},\\\
\beta^{i,y}&=n^{i}_{y}B^{2,i}_{12}+t^{i}_{y}B^{2,i}_{13},\end{split}$ (54)

we obtain for $V$

${\tiny
V=\setcounter{MaxMatrixCols}{12}\begin{bmatrix}1&0&0&0&0&0&0&0&0&0&0&0\\\
0&\tfrac{\partial x}{\partial\hat{x}}&\tfrac{\partial
y}{\partial\hat{x}}&0&0&0&0&0&0&0&0&0\\\ 0&\tfrac{\partial
x}{\partial\hat{y}}&\tfrac{\partial y}{\partial\hat{y}}&0&0&0&0&0&0&0&0&0\\\
0&0&0&1&0&0&0&0&0&0&0&0\\\ 0&0&0&0&\tfrac{\partial
x}{\partial\hat{x}}&\tfrac{\partial y}{\partial\hat{x}}&0&0&0&0&0&0\\\
0&0&0&0&\tfrac{\partial x}{\partial\hat{y}}&\tfrac{\partial
y}{\partial\hat{y}}&0&0&0&0&0&0\\\ 0&0&0&0&0&0&1&0&0&0&0&0\\\
0&0&0&0&0&0&0&\tfrac{\partial x}{\partial\hat{x}}&\tfrac{\partial
y}{\partial\hat{x}}&0&0&0\\\ 0&0&0&0&0&0&0&\tfrac{\partial
x}{\partial\hat{y}}&\tfrac{\partial y}{\partial\hat{y}}&0&0&0\\\
0&0&0&0&-\beta^{1,x}&-\beta^{1,y}&0&\beta^{1,x}&\beta^{1,y}&B^{2,1}_{11}&0&0\\\
0&-\beta^{2,x}&-\beta^{2,y}&0&0&0&0&\beta^{2,x}&\beta^{2,y}&0&B^{2,2}_{11}&0\\\
0&-\beta^{3,x}&-\beta^{3,y}&0&\beta^{3,x}&\beta^{3,y}&0&0&0&0&0&B^{2,3}_{11}\\\
\end{bmatrix}.}$ (55)

The same considerations lead to a similar derivation of $E$, $V^{c}$, and $D$
for the robust element, resulting in

${\tiny
V=\setcounter{MaxMatrixCols}{15}\begin{bmatrix}1&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&\tfrac{\partial x}{\partial\hat{x}}&\tfrac{\partial
y}{\partial\hat{x}}&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&\tfrac{\partial
x}{\partial\hat{y}}&\tfrac{\partial
y}{\partial\hat{y}}&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&1&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&\tfrac{\partial
x}{\partial\hat{x}}&\tfrac{\partial y}{\partial\hat{x}}&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&\tfrac{\partial x}{\partial\hat{y}}&\tfrac{\partial
y}{\partial\hat{y}}&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&1&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&\tfrac{\partial x}{\partial\hat{x}}&\tfrac{\partial
y}{\partial\hat{x}}&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&\tfrac{\partial
x}{\partial\hat{y}}&\tfrac{\partial y}{\partial\hat{y}}&0&0&0&0&0&0\\\
0&0&0&-B^{1,1}_{12}&0&0&B^{1,1}_{12}&0&0&B^{1,1}_{11}&0&0&0&0&0\\\
-B^{1,2}_{12}&0&0&0&0&0&B^{1,2}_{12}&0&0&0&B^{1,2}_{11}&0&0&0&0\\\
-B^{1,3}_{12}&0&0&B^{1,3}_{12}&0&0&0&0&0&0&0&B^{1,3}_{11}&0&0&0\\\
0&0&0&0&-\beta^{1,x}&-\beta^{1,y}&0&\beta^{1,x}&\beta^{1,y}&0&0&0&B^{2,1}_{11}&0&0\\\
0&-\beta^{2,x}&-\beta^{2,y}&0&0&0&0&\beta^{2,x}&\beta^{2,y}&0&0&0&0&B^{2,2}_{11}&0\\\
0&-\beta^{3,x}&-\beta^{3,y}&0&\beta^{3,x}&\beta^{3,y}&0&0&0&0&0&0&0&0&B^{2,3}_{11}\\\
\end{bmatrix}}$ (56)

for $V$, where $\beta$ is as defined in (54).

##  4 Discretisation

We now describe the discretisations of the Hamiltonian system (8) using a
function space introduced in the previous section. Smooth solutions to (8a)
generates the following curve of diffeomorphisms:

$\dot{\varphi}_{t}=u_{t}\circ\varphi_{t},\quad\varphi_{0}=\textbf{id},$ (57)

where the domain of $\varphi_{t}$ is $\Omega_{0}$. This subsumes the left
action on the curve $q_{0}$ in (8b). Our approach is therefore to solve (8a)
for an _outer metric_ in tandem with integrating the diffeomorphism defined
over the entire domain and moving the mesh, thereby automatically providing a
solution to a discrete analogue of (8b). We denote by $\mathcal{T}_{0}$ denote
a shape-regular, quasi-uniform triangulation of the template domain
$\Omega_{0}$. Let $\mathcal{E}_{h}$ denote the _mesh skeleton_ of
$\mathcal{T}_{0}$ and $\mathring{\mathcal{E}}_{h}$ the subset of
$\mathcal{E}_{h}$ whose elements do not intersect $\partial\Omega_{0}$. We
place the following assumption on the initial triangulation $\mathcal{T}_{0}$:

######  Assumption 1.

$\mathcal{T}_{0}$ is constructed such that the range of $q_{0}$ is described
by a subset of $\mathring{\mathcal{E}}_{h}$.

Using the definition in (21) we define the vector-valued Wu-Xu space defined
over $\Omega_{0}$:

$V(\Omega_{0})=\\{v\in
L^{2}(\Omega)^{2}\,|\,v_{i}|_{K}\in\mathcal{W}(K),\,K\in\mathcal{T}_{0},\,i=1,2\\}.$

Further, let $0=t_{0},t_{\Delta T},\ldots t_{T-1}=1$ denote $T$ uniformly
distributed points and use an Euler discretisation of the time derivate in
(5), where we let $\varphi_{t_{k}}\in V(\Omega_{0})$, $u_{t_{k}}\in
V(\Omega_{0})$:

$\varphi_{t_{k+1}}=\varphi_{t_{k}}+u_{t_{k}}\circ\varphi_{t_{k}}\Delta T.$
(58)

For sufficiently small $\Delta T$, $\varphi_{t_{k}}$ is a diffeomorphism of
$\Omega$ [44]. Using the notation
$\Omega_{t_{k}}=\varphi_{t_{k}}\circ\Omega_{0}$,
$V(\Omega_{t_{k}}):=\\{f\,|\,f\circ\varphi_{t_{k}}^{-1}\in V(\Omega_{0})\\}$
and by noting that $q_{t_{k}}=\varphi_{t_{k}}\circ q_{0}$ we obtain a discrete
analogue of (8) where $\hat{u}_{t_{k}}\in V(\Omega_{t_{k}})$:

$\displaystyle
a_{\Omega_{t_{k}}}(\hat{u}_{t_{k}},\hat{v})=\int_{S^{1}}\nabla\varphi_{t_{k}}^{-\top}\circ
q_{0}\mathbf{n}_{q_{0}}\tilde{p}_{0}\cdot\hat{v}\circ q_{0}\,\text{d}{s},$
$\displaystyle\forall\hat{v}\in V(\Omega_{t_{k}}),$ (59a)
$\displaystyle\varphi_{t_{k+1}}=\varphi_{t_{k}}+\hat{u}_{t_{k}}\Delta T,$
(59b)

for $k=0,\ldots,T-1$, where
$\varphi_{t_{k}}\circ\partial\Omega_{0}=\textbf{id}$ owing to the homogeneous
Dirichlet boundary conditions implied by (59a). At each time step $k$ after
the solution of (59a), the mesh is moved according to (58) upon which the
equation (60):

$q_{t_{k+1}}=q_{t_{k}}+u_{t_{k}}\circ q_{t_{k}}\Delta T,$ (60)

is automatically satisfied. The underlying coordinate field of the mesh itself
is chosen to be a Lagrange subspace of $V(\Omega_{0})$, so that the map
$q_{0}\mapsto\varphi_{t_{k}}\circ q_{0}$ is a diffeomorphism. At $k=0$, the
assembly of the right-hand side is in practice done by integration over
$q_{0}\circ S^{1}$, which means that we can supply an initial “momentum”
signal $\mathbf{p}_{0}\in\tilde{p}_{0}\circ q_{0}^{-1}\in L^{2}(q_{0}\circ
S^{1})$ (now defined over initial curve) to encode the entire geodesic flow of
$\varphi$, and thereby of the embedded curve. Figure 4 show examples of
forward integration of this system for various $\mathbf{p}_{0}$ and
$q_{0}\circ S^{1}$ (the initial meshes were generated using gmsh [45]). Note
that the norm of the velocity present in (8) is confined to certain energy
levels determined by the initial momentum as the system is integrated. In the
fully discrete analogue we can only hope to establish approximate conservation
of the Hamiltonian. The importance of this nebulous since we only integrate
over fixed time intervals, and is subject to future work.

The computational cost of integrating (59) is dominated by the inversion of
the discrete bilinear form. Mesh-based methods readily facilitate parallel
computations (e.g. matrix-vector products in a Krylov subspace method), which
along with preconditioning strategies are competitive with fast multipole
methods. They also offer flexibility in choosing bilinear form (which can be
altered according to an informed modelling choice or application). Finally,
mesh adaptivity is also an option. For the application at hand a _graded_ mesh
with a fine resolution in the vicinity of the curve and coarser elements
closer to the boundary can both increase accuracy and the computational burden
of the method.

Figure 4: Top row: template domains $\Omega_{0}$ with curves highlighted.
Bottom row: $\varphi_{1}\circ\Omega_{1}$ computed by integrating (59) from two
templates with different initial momenta.

##  5 Inverse problem

We now consider the matching problem using the data misfit functional in (2).
We wish to estimate the momentum $\mathbf{p}_{0}:=\tilde{p}_{0}\circ
q_{0}^{-1}\in\mathbf{P}$ that generates the curve $t\mapsto q_{t}$. That is,
$\mathbf{p}_{0}$ is the momentum object defined on the computational domain
$q_{0}\circ S^{1}$. We drop explicit dependence on the template as it remains
fixed during computation as well as the time dimension of the initial
momentum. To ease the notation we use boldface $\boldsymbol{\tau}$ to
represent the smoothed version of the indicator function on the interior of a
curve $q$ the i.e.:

$\boldsymbol{\tau}=\mathcal{C}^{-1}\mathbb{1}_{q}.$

We define the _forward operator_ :

$\mathbf{p}\mapsto\mathcal{F}(\mathbf{p})=\boldsymbol{\tau}:=\mathcal{C}^{-1}\mathbb{1}_{q_{1}},$
(61)

where $q_{1}$ is the solution at $t=1$ given by solving (59) using $q_{0}$ and
$\mathbf{p}$ as initial conditions, i.e. the time-$1$ flow map of the initial
curve. Given a target shape $\boldsymbol{\tau}_{\textnormal{target}}$ the
inverse problem of interest is therefore to recover the momentum
$\mathbf{p}^{*}$ such that:

$\boldsymbol{\tau}_{\textnormal{target}}\approx\mathcal{F}(\mathbf{p}^{*})+\xi,$

where the noise $\xi\in\mathcal{N}(0,\mathcal{C})$ is a Gaussian measure with
mean zero and covariance operator $\mathcal{C}$ defined later on.

To tackle this inverse problem, we use an Ensemble Kalman iteration. We let
$N$ denote the number of ensemble members and let $\mathbf{p}^{j}$,
$j=1,\ldots,N$ denote the momenta corresponding to ensemble member $j$. The
ensemble mean momentum and the mean predicted shape are:

$\displaystyle\bar{\mathbf{p}}:=N^{-1}\sum_{j=1}^{N}\mathbf{p}_{j},\qquad\bar{\boldsymbol{\tau}}:=N^{-1}\sum_{j=1}^{N}\boldsymbol{\tau}^{j},$
(62)

where
$\boldsymbol{\tau}^{j}=\mathcal{C}^{-1}\mathbb{1}_{\mathcal{F}(\mathbf{p}^{j})}$.
The Kalman update operator is defined by:

$\displaystyle\mathfrak{K}_{\mathbf{P}}=\text{Cov}_{\mathbf{P}Q}[\text{Cov}_{QQ}+\xi
I]^{-1},$ (63)

where $\xi$ is a regularisation parameter described later and $I$ is the
identity operator. The actions above are given by:

$\displaystyle\text{Cov}_{QQ}[\cdot]=\frac{1}{N-1}\sum_{j=1}^{N}(\boldsymbol{\tau}^{j}-\bar{\boldsymbol{\tau}})\langle\boldsymbol{\tau}^{j}-\bar{\boldsymbol{\tau}},\cdot\rangle_{L^{2}},$
(64a)
$\displaystyle\text{Cov}_{\mathbf{P}Q}[\cdot]=\frac{1}{N-1}\sum_{j=1}^{N}(\mathbf{p}^{j}-\bar{\mathbf{p}})\langle\boldsymbol{\tau}^{j}-\bar{\boldsymbol{\tau}},\cdot\rangle_{L^{2}}.$
(64b)

The data misfit at iteration $k$ of the EKI is defined as:

$\mathfrak{E}^{k}=\|\boldsymbol{\tau}_{\textnormal{target}}-\bar{\boldsymbol{\tau}}\|_{L^{2}(\Omega)}^{2}.$
(65)

The prediction and analysis steps are summarised below:

  1. 1.

_Prediction_ : For each ensemble member $j$, compute
$\boldsymbol{\tau}^{j}=\mathcal{C}^{-1}\mathbb{1}_{\mathcal{F}(\mathbf{p}^{j})}$
and the average $\bar{\boldsymbol{\tau}}$ using (62).

  2. 2.

_Analysis_ : Update each ensemble momentum:

$\mathbf{p}^{j+1}=\mathbf{p}^{j}+\mathfrak{K}_{\mathbf{P}}(\boldsymbol{\tau}_{\textnormal{target}}-\boldsymbol{\tau}^{j}).$

##  6 Numerical experiments

We present numerical experiments showing that the ensemble Kalman inversion
(EKI) algorithm is able to find suitable approximations of a target given a
random initial ensemble. Section 6.1 describes how we generate the synthetic
target that we will use as matching targets. Section 6.2 summarises the
parameters that we have chosen to in our experiments to match the synthetic
data, and section 6.3 contains the numerical results.

###  6.1 Synthetic data

For simplicity we fix the template curve throughout our experiments and choose
the unit circle. The initial mesh is that shown in the top left vignette of
Figure 4. The computation domain $\Omega_{0}$ is a triangulation of
$[-10,10]^{2}$ with mesh resolution222This is the maximum diameter $h_{K}$ of
any triangle $K$ in the triangulation. $h=1$. We have taken $\alpha=0.5$ in
(9), $T=15$ time steps and have used piecewise constant finite elements on the
mesh skeleton to represent $\mathbf{P}$ (although we compute only with
functions supported over the submanifold $q_{0}\circ
S^{1}\subset\mathcal{E}_{h}$). We use the forward operator described
previously to generate synthetic targets for this set of parameters. Applying
the forward operator $\mathcal{F}$ to the momenta in (66) below we produce the
targets seen in Figure 5.

$\displaystyle\mathbf{p}^{\dagger}_{\text{contract}}=-1.38\pi,$ (66a)
$\displaystyle\mathbf{p}^{\dagger}_{\text{squeeze}}=\begin{cases}0.83\pi
e^{-y^{2}/5}&x<-0.3\\\
\frac{5}{3}\pi\sin(x/5)|y|&\textnormal{otherwise}\end{cases},$ (66b)
$\displaystyle\mathbf{p}^{\dagger}_{\text{star}}=2.6\pi\cos(2\pi x/5),$ (66c)
$\displaystyle\mathbf{p}^{\dagger}_{\text{teardrop}}=\begin{cases}-3\pi\,\textnormal{sign}(y)&y<0\\\
3\pi e^{-x^{2}/5}&\textnormal{otherwise}\end{cases}.$ (66d)

Figure 5: Synthetic targets generated using the momentum in (66). A simple
contraction of the shape, a squeezed figure with a right bias, a star and a
teardrop shape).

With $\mathbf{p}^{\dagger}$ we associate the following relative error at each
iterate $k$:

$\mathcal{R}^{k}=\|\bar{\mathbf{p}}^{k}-\mathbf{p}^{\dagger}\|_{L^{2}(q_{0}\circ
S^{1})}/\|\mathbf{p}^{\dagger}\|_{L^{2}(q_{0}\circ S^{1})}.$ (67)

The _consensus deviation_ $\mathcal{S}^{k}$ of an ensemble at iteration $k$ in
equation (68) is defined below:

$\mathcal{S}^{k}=N^{-1}\sum_{j=1}^{N}\|\mathbf{p}^{j,k}-\bar{\mathbf{p}}^{k}\|_{L^{2}(q_{0}\circ
S^{1})},$ (68)

where $\mathbf{p}^{j,k}$ denotes the momentum of ensemble member $j$ at
iteration $k$. This quantity is a useful diagnostic which measures the
information remaining in the ensemble after iteration $k$. Since EKI relies on
estimates of the forecast covariance, consensus is reached when
$\mathcal{S}^{k}$ approaches zero, at which point the algorithm can be
stopped.

In all our simulations we invert the system in (59a) using the direct solver
MUMPS [46]; investigating a preconditioned iterative solver is subject to
future work. For details on the validation of the implementation of the Wu-Xu
element in Firedrake, see [47, Appendix B] and the Zenodo entry [48].

###  6.2 Experimental setup

We now describe the setup we have used to test the EKI. Firstly, we have taken
$T=10$ and $\alpha=1$ so the parameters differ from those used to generate the
synthetic targets. Recall that EKI requires an initial ensemble, in this case
of momenta. The basis coefficients of the momenta was sampled from a uniform
distribution over the interval $[-25,25]$, with different realisations for
each ensemble member. The parameter $\xi$ in (63) determines the ratio between
the influence of the prediction covariance on the Kalman gain. We set
$\xi=10^{-3}$ in (63), although adaptive tuning of this parameter to avoid
overfitting is possible; an early termination rule is suggested in [32,
Equation 10]. We choose $\mathcal{C}=(\textbf{id}-\kappa\Delta)^{-1}$,
$\kappa=10$ in (65), as this smoothes out the mismatch sufficiently for our
computational domain. The quality of the matching is directly related to the
size and variance of the ensemble as the solution is sought as a linear
combination of its members. We conduct experiments for two ensemble sizes,
$N=20$, $N=40$ and $N=80$. These were chosen since, with the parameter set as
above, $\dim\mathbf{P}=48$ in order to develop an understanding of how EKI
performs when the ensemble size is smaller than, similar to and larger than
the dimension of the state, while still keeping the overall computational cost
such that the experiments can be done in a reasonable amount of time. The case
where $N<<\dim\mathbf{P}$ is the de facto situation for ensemble methods as
the MC approximation allows for a computationally feasible method. In general,
small ensemble sizes can lead to filter inbreeding (the forecast covariance is
underestimated), filter divergence (the gain does not adequately correct the
ensemble), or spurious correlations [49]. We comment on each of these later.

###  6.3 Results

We have run EKI $10$ times for each value of $N$ with different draws for the
initial ensemble to assess the robustness of the algorithm with respect to the
starting point. Figure 6 shows examples of the numerical results we obtain for
curve matching using EKI. Note that only five iterations of EKI were necessary
to produce the results shown in this section to reach a relative tolerance
below $3\%$. Qualitatively a larger ensemble size leads to a better match,
which is to be expected. Ensemble methods such as the EKI offer a practical
advantage to gradient methods given their inherent parallelisability. Indeed,
the _prediction_ step discussed in Section 5 can be done in parallel for each
ensemble member. We therefore start $N$ processes corresponding to each
member, and used a Message Passing Interface (MPI) [50] implementation to
exchange information between them (the MPI _reduce_ operation, specifically).
Thanks to this parallelisation, five iterations of EKI takes less than two
minutes for $N=20$, five minutes for $N=40$ and 14 minutes for $N=80$ on a
2021 MacBook Pro333Apple M1 Pro chip, 16 GB of memory..

Figure 6: Final reconstruction of the targets in figures 5 produced by EKI for
$N=20$ (top row, $N=40$ (middle row) and $N=80$ (bottom row).

Figure 7, 8 and 9 show the relative errors, data misfits and consensus
deviations for our experiments across the selected targets and ensemble sizes.
These all decrease at various rates in the early iterations after which they
stagnate. As the Kalman gain corrects the ensemble members, and therefore the
motion of their respective curves, the data misfit decreases meaning that each
member improves its prediction. This increases consensus in the ensemble,
which explains what is seen in figure 9. We notice from the data misfits and
the momentum consensus that higher values of $N$ provides a more accurate
approximation of the true momentum, which explains the accuracy of the matches
seen in figure 6. Note that the relative error, which is a surrogate for
posterior consistency, also decreases (albeit not with a clear pattern across
shapes and ensemble sizes). Since the forward operator in use is heavily
nonlinear, theoretical convergence results are not readily obtained at this
stage. We highlight that the Kalman gain is very efficient in correcting the
ensemble momenta, with the consensus decreasing exponentially in the early
stages of the algorithm. We comment on this later. The conclusion is therefore
that even at a modest ensemble size, EKI performs well. It is not certain that
the same behaviour that we see above (i.e. few iterations are needed) will
scale with $N$ and the size of the problem, but the results are promising for
research in this direction.

Higher values of $\xi$ were found to slow the convergence of the algorithm
compared to the selected value, which is consistent with the behaviour for
landmark-based EKI [30], and we do not comment on it further. We noticed that
the value of $\kappa$ also influences the convergence of the EKI; for small
values the operator $\textbf{id}-\kappa\Delta$ approaches the identity, and
since the mismatch $X$ is computed from point evaluations of the finite
element function, information can be lost if the grid is not sufficiently
refined. A larger value of $\kappa$ “spreads out” the mismatch which improves
convergence for coarser grids.

Figure 7: Relative error (equation (67)). The rows corresponds to $N=20$,
$N=40$ and $N=80$ and the columns to the targets in figure 5.

Figure 8: Data misfits (equation (65)). The rows corresponds to $N=20$, $N=40$
and $N=80$ and the columns to the targets in figure 5.

Figure 9: Momentum consensus deviation (equation (68)). The rows corresponds
to $N=20$, $N=40$ and $N=80$ and the columns to the targets in figure 5.

##  7 Summary and outlook

In this paper we have presented a parameterisation- and derivative-free method
for matching closed planar curves. A moving mesh discretisation of Hamilton’s
equations for curves was described using the induced diffeomorphism of a
vector field occupying the Wu-Xu finite element space. We also describe a
transformation theory for this element facilitating a computationally
performant forward model for use in the associated inverse problem. Finding
the momentum encoding the forward motion of the template matches a desired
curve was treated as a Bayesian inverse problem in section 5 and EKI was used
to approximate its solution. The numerical results presented in section 6
suggests that the method shows great promise. Not only does is it easy to
implement, the EKI is shown to quickly reach ensemble consensus meaning that
it is efficient in exploring the subspace spanned by the initial ensemble.
This is in part thanks to the momentum being a one-dimensional signal on the
template. Treating the mismatch term in a negative Sobolev norm was shown to
increase both accuracy of our results and robustness to mesh resolution. We
also showed that the method is robust to the choice of initial ensemble even
when the ensemble size is less than half the dimension of the forward problem.
Further, assuming the forward operator is scalable as the mesh is refined
(large-scale PDE solves are common in many areas of scientific computing, and
the inverse needed in the Kalman gain scales cubically in $N$ [51]).

Future work includes proving convergence of the finite element discretisation
for (59) and subsequently using these error estimates to quantify error in a
rigorous treatment of the Bayesian inverse problem [52]. As indicated in [53],
some challenges exist for nonconforming finite element methods with singular
source terms. The template considered in this paper is a piece-wise linear
curve. An obvious extension would be to apply isoparametric methods to cater
for piece-wise higher-order polynomial curves. The effect of this would only
affect the right-hand side and would not affect regularity results for the
velocity. An advantage of the finite element method for curves is also that it
allows for adaptivity e.g. refinement of the mesh only in the vicinity of the
embedded template. We considered problems of modest size to illustrate the
discretisation and the EKI. As the mesh is refined, it is likely the case that
the dimension of the forward operator dwarfs the size of the ensemble and
effects of the MC approximation are more pronounced. This is the case for
ensemble methods for e.g. numerical weather prediction and several techniques
exist to counter these effects [49] (e.g. localisation or covariance
inflation). In particular, localisation methods may be suitable to assume
conditional independence between separated states (i.e. parts of the shape
that are distant in physical space) so as to counter spurious correlations.

##  Appendix A Proof of Theorem 1

The momentum satisfies

$\dot{p}_{t}+\nabla u_{t}^{\top}\circ q_{t}p_{t}=0$

Using the ansatz we verify:

$\displaystyle\dot{p}_{t}+\nabla u_{t}^{\top}\circ q_{t}p_{t}$
$\displaystyle=\dot{J_{t}^{-\top}}p_{0}+\nabla u_{t}^{\top}J_{t}^{-\top}p_{0}$
$\displaystyle=-J_{t}^{-\top}d(J_{t}^{\top})J_{t}^{-\top}p_{0}+\nabla
u_{t}^{\top}\circ q_{t}J_{t}^{-\top}p_{0}$
$\displaystyle=-J_{t}^{-\top}(dJ_{t})^{\top}J_{t}^{-\top}p_{0}+\nabla
u_{t}^{\top}\circ q_{t}J_{t}^{-\top}p_{0}$
$\displaystyle=-J_{t}^{-\top}(\nabla u_{t}\circ
q_{t}J_{t})^{\top}J_{t}^{-\top}p_{0}+\nabla u_{t}^{\top}\circ
q_{t}J_{t}^{-\top}p_{0}$ $\displaystyle=-J_{t}^{-\top}J_{t}^{\top}\nabla
u_{t}^{\top}\circ q_{t}J_{t}^{-\top}p_{0}+\nabla u_{t}^{\top}\circ
q_{t}J_{t}^{-\top}p_{0}$ $\displaystyle=-\nabla
u_{t}^{\top}J_{t}^{-\top}p_{0}+\nabla u_{t}^{\top}\circ
q_{t}J_{t}^{-\top}p_{0}$ $\displaystyle=0\,.$

## References

  * [1] L. Younes, Shapes and diffeomorphisms, Vol. 171, Springer, 2010.
  * [2] P. W. Michor, D. Mumford, An overview of the Riemannian metrics on spaces of curves using the Hamiltonian approach, Applied and Computational Harmonic Analysis 23 (1) (2007) 74–113.
  * [3] M. Bauer, M. Bruveris, P. W. Michor, Overview of the geometries of shape spaces and diffeomorphism groups, Journal of Mathematical Imaging and Vision 50 (1-2) (2014) 60–97.
  * [4] C. J. Cotter, D. D. Holm, Geodesic boundary value problems with symmetry, arXiv preprint arXiv:0911.2205 (2009).
  * [5] M. Bauer, M. Bruveris, A new Riemannian setting for surface registration, arXiv preprint arXiv:1106.0620 (2011).
  * [6] J. Glaunès, A. Qiu, M. I. Miller, L. Younes, Large deformation diffeomorphic metric curve mapping, International journal of computer vision 80 (3) (2008) 317.
  * [7] S. Wu, J. Xu, Nonconforming finite element spaces for $2m^{\text{th}}$ order partial differential equations on $\mathbb{R}^{n}$ simplicial grids when $m=n+1$, Mathematics of Computation 88 (316) (2019) 531–551.
  * [8] F. Rathgeber, D. A. Ham, L. Mitchell, M. Lange, F. Luporini, A. T. McRae, G.-T. Bercea, G. R. Markall, P. H. Kelly, Firedrake: automating the finite element method by composing abstractions, ACM Transactions on Mathematical Software (TOMS) 43 (3) (2016) 1–27.
  * [9] R. C. Kirby, A general approach to transforming finite elements, The SMAI journal of computational mathematics 4 (2018) 197–224.
  * [10] R. C. Kirby, L. Mitchell, Code generation for generally mapped finite elements, ACM Transactions on Mathematical Software (TOMS) 45 (4) (2019) 1–23.
  * [11] U. Grenander, M. I. Miller, Computational anatomy: An emerging discipline, Quarterly of applied mathematics 56 (4) (1998) 617–694.
  * [12] P. Dupuis, U. Grenander, M. I. Miller, Variational problems on flows of diffeomorphisms for image matching, Quarterly of Applied Mathematics (1998) 587–600.
  * [13] J. Glaunes, A. Trouvé, L. Younes, Diffeomorphic matching of distributions: A new approach for unlabelled point-sets and sub-manifolds matching, in: Proceedings of the 2004 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2004. CVPR 2004., Vol. 2, IEEE, 2004, pp. II–II.
  * [14] M. Vaillant, J. Glaunès, Surface matching via currents, in: Biennial International Conference on Information Processing in Medical Imaging, Springer, 2005, pp. 381–392.
  * [15] M. Bauer, M. Bruveris, N. Charon, J. Møller-Andersen, A relaxed approach for curve matching with elastic metrics, ESAIM: Control, Optimisation and Calculus of Variations 25 (2019) 72.
  * [16] M. Bauer, M. Bruveris, P. Harms, J. Møller-Andersen, Curve matching with applications in medical imaging, arXiv preprint arXiv:1506.08840 (2015).
  * [17] M. Bauer, M. Bruveris, P. Harms, J. Møller-Andersen, A numerical framework for Sobolev metrics on the space of curves, SIAM Journal on Imaging Sciences 10 (1) (2017) 47–73.
  * [18] E. Hartman, Y. Sukurdeep, E. Klassen, N. Charon, M. Bauer, Elastic shape analysis of surfaces with second-order Sobolev metrics: a comprehensive numerical framework, arXiv preprint arXiv:2204.04238 (2022).
  * [19] K. Bharath, S. Kurtek, Analysis of shape data: From landmarks to elastic curves, Wiley Interdisciplinary Reviews: Computational Statistics 12 (3) (2020) e1495.
  * [20] C. J. Cotter, The variational particle-mesh method for matching curves, Journal of Physics A: Mathematical and Theoretical 41 (34) (2008) 344003.
  * [21] C. J. Cotter, A. Clark, J. Peiró, A reparameterisation based approach to geodesic constrained solvers for curve matching, International Journal of Computer Vision 99 (1) (2012) 103–121.
  * [22] A. Günther, H. Lamecker, M. Weiser, Flexible shape matching with finite element based LDDMM, International Journal of Computer Vision 105 (2) (2013) 128–143.
  * [23] K. Höllig, Finite element methods with B-splines, Vol. 26, SIAM, 2003.
  * [24] A. Paganini, F. Wechsung, P. E. Farrell, Higher-order moving mesh methods for PDE-constrained shape optimization, SIAM Journal on Scientific Computing 40 (4) (2018) A2356–A2382.
  * [25] A. Paganini, F. Wechsung, Fireshape: a shape optimization toolbox for Firedrake, Structural and Multidisciplinary Optimization 63 (2021) 2553–2569.
  * [26] F. X. Vialard, L. Risser, D. Rueckert, C. J. Cotter, Diffeomorphic 3D image registration via geodesic shooting using an efficient adjoint calculation, International Journal of Computer Vision 97 (2) (2012) 229–241.
  * [27] A. Bock, A. Arnaudon, C. Cotter, Selective metamorphosis for growth modelling with applications to landmarks, in: International Conference on Geometric Science of Information, Springer, 2019, pp. 39–48.
  * [28] M. I. Miller, A. Trouvé, L. Younes, Geodesic shooting for Computational Anatomy, Journal of Mathematical Imaging and Vision 24 (2) (2006) 209–228.
  * [29] C. J. Cotter, S. L. Cotter, F.-X. Vialard, Bayesian data assimilation in shape registration, Inverse Problems 29 (4) (2013) 045011.
  * [30] A. Bock, C. J. Cotter, Learning landmark geodesics using the ensemble Kalman filter, Foundations of Data Science 3 (4) (2021) 701.
  * [31] M. A. Iglesias, K. J. Law, A. M. Stuart, Ensemble Kalman methods for inverse problems, Inverse Problems 29 (4) (2013) 45001.
  * [32] M. A. Iglesias, A regularizing iterative ensemble Kalman method for PDE-constrained inverse problems, Inverse Problems 32 (2) (2016) 25002.
  * [33] J. Glaunès, A. Trouvé, L. Younes, Modeling planar shape variation via Hamiltonian flows of curves, in: Statistics and analysis of shapes, Springer, 2006, pp. 335–361.
  * [34] O. Lakkis, T. Pryer, An adaptive finite element method for the infinity Laplacian, in: Numerical Mathematics and Advanced Applications-ENUMATH 2013, Springer, 2015, pp. 283–291.
  * [35] E. Barron, L. Evans, R. Jensen, The infinity Laplacian, Aronsson’s equation and their generalizations, Transactions of the American Mathematical Society 360 (1) (2008) 77–101.
  * [36] W. P. Ziemer, Weakly differentiable functions: Sobolev spaces and functions of bounded variation, Vol. 120, Springer Science & Business Media, 2012.
  * [37] A. Ern, J. L. Guermond, Theory and practice of finite elements, Vol. 159, Springer Science & Business Media, 2013.
  * [38] R. A. Adams, J. J. Fournier, Sobolev spaces, 2nd Edition, Elsevier, 2003.
  * [39] S. C. Brenner, L. R. Scott, L. R. Scott, The mathematical theory of finite element methods, Vol. 3, Springer, 2008.
  * [40] R. C. Kirby, Algorithm 839: FIAT, a new paradigm for computing finite element basis functions, ACM Transactions on Mathematical Software 30 (4) (2004) 502–516.  doi:10.1145/1039813.1039820.
  * [41] J. H. Argyris, I. Fried, D. W. Scharpf, The TUBA family of plate elements for the matrix displacement method, Aeronautical Journal 72 (1968) 701–709.  doi:10.1017/S000192400008489X.
  * [42] V. Domínguez, F.-J. Sayas, Algorithm 884: A simple Matlab implementation of the Argyris element, ACM Transactions on Mathematical Software 35 (2) (2008) 16\.  doi:10.1145/1377612.1377620.
  * [43] M. E. Rognes, R. C. Kirby, A. Logg, Efficient assembly of $H(\mathrm{div})$ and $H(\mathrm{curl})$ conforming finite elements, SIAM Journal on Scientific Computing 31 (6) (2010) 4130–4151.
  * [44] G. Allaire, M. Schoenauer, Conception optimale de structures, Vol. 58, Springer, 2007.
  * [45] C. Geuzaine, J. F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities, International Journal for Numerical Methods in Engineering 79 (11) (2009) 1309–1331.
  * [46] P. R. Amestoy, I. S. Duff, J.-Y. L’Excellent, J. Koster, MUMPS: a general purpose distributed memory sparse solver, in: International Workshop on Applied Parallel Computing, Springer, 2000, pp. 121–130.
  * [47] A. Bock, On inverse problems and finite elements in shape analysis, Ph.D. thesis, Imperial College London (2020).
  * [48] Software used in ‘On inverse problems and finite elements in shape analysis’ (2020).  doi:10.5281/zenodo.3832141.
URL https://doi.org/10.5281/zenodo.3832141

  * [49] R. Petrie, Localization in the ensemble Kalman filter, MSc Thesis, University of Reading (2008).
  * [50] W. Gropp, W. D. Gropp, E. Lusk, A. Skjellum, A. D. F. E. E. Lusk, Using MPI: portable parallel programming with the message-passing interface, Vol. 1, MIT press, 1999.
  * [51] J. Mandel, Efficient implementation of the ensemble Kalman filter, University of Colorado at Denver and Health Sciences Center, Center for Computational Mathematics, 2006.
  * [52] S. L. Cotter, M. Dashti, A. M. Stuart, Approximation of Bayesian inverse problems for PDEs, SIAM Journal on Numerical Analysis 48 (1) (2010) 322–345.
  * [53] A. Bock, C. Cotter, A note on error analysis for a nonconforming discretisation of the tri-Helmholtz equation with singular data, arXiv preprint arXiv:2101.04181 (2021).