\section{Introduction}
Thin elastic objects provide a fascinating range of applications, particularly
in nanotechnology, where small lightweight constructions can undergo large
deformations that are controlled by external stimuli, cf., e.g.,~\cite{SZTDI12}.
An important class of
these objects are elastic ribbons which are thin elastic bands of small width.
A characteristic feature is that these bands are unshearable and hence curvature
can only occur orthogonal to the strip's inplane direction.  A particular
configuration that has been addressed in the literature are M\"obius strips of
small width, cf.~\cite{Wund62,StaHei07,BarHor15}.

To effectively simulate the highly nonlinear behavior of thin elastic objects
dimensionally reduced mathematical descriptions are necessary. In the case of
ribbons a one-dimensional model has been proposed by Sadowsky in~\cite{Sado30a,Sado30b}
which was obtained via formal arguments; further justifications have been given
in~\cite{Wund62,KirFri15}. A complete, mathematically rigorous derivation
from three-dimensional hyperelasticity has been carried out in~\cite{FHMP16}, which
shows that a correction of Sadowsky's functional is necessary. The correction
concerns the important case of vanishing curvature and thereby provides a more
general description. Additionally, this eliminates a singularity in
Sadowsky's functional which is particularly important for the stable
discretization and iterative solution. The nonsmooth character of the
corrected functional underlines the expected feature that the behavior of ribbons is
substantially more complicated than that of elastic rods, cf.~\cite{MorMul03}
for a corresponding rigorously derived mathematical model. Our numerical experiments
confirm that the correction in the one-dimensional functional is relevant.
A comparison of certain instabilities of twisted rods and ribbons is provided in~\cite{AudSef15}.
The presence of such effects supports the importance of developing stable numerical
approximation schemes.

The model identified in~\cite{FHMP16} describes the deformation of an elastic ribbon
of length $L>0$ via a frame $(y,b,d)$ consisting of functions
\[
y,b,d : (0,L) \to \mathbb{R}^3,
\]
defined on the centerline $(0,L)$ of the undeformed ribbons. The frame condition
requires the functions to satisfy
\[
[y',b,d] \in SO(3)
\]
pointwise in $(0,L)$ with the set of proper orthogonal matrices~$SO(3)$.
In particular, $y$ provides an arclength parametrization of
the deformed centerline. The arclength property $|y'|^2=1$ models the physical
feature that thin elastic objects are inextensible in the bending regime.
The directors~$b$ and~$d$ define the deformation of the band orthogonal to the
centerline. Since the band cannot be sheared, the in-plane curvature
vanishes, i.e., the additional condition
\[
y'' \cdot b = 0
\]
occurs, which in view of the relation $y'\cdot b=0$ is equivalent to the
condition $y'\cdot b'=0$. A frame $(y,b,d)$ describing the deformation of an
elastic band for given boundary conditions minimizes the dimensionally reduced
elastic energy
\[
E[y,b,d] = \int_0^L \overline{Q}(y''\cdot d,b'\cdot d) \dv{x}.
\]
Appropriately scaled body forces can be included in the minimization problem.
Because of the frame condition and the constraint $y''\cdot b = 0$,
the functional can be expressed in terms of the curvature $\kappa$ and torsion $\tau$
of the curve parametrized by~$y$ and described by the frame $(y,b,d)$. These
quantities are given by
\[
\kappa^2 = (y''\cdot d)^2 = |y''|^2, \quad \tau^2 = (b'\cdot d)^2 = |b'|^2
\]
and we may consider the equivalent functional
\[
E[y,b] = \int_0^L \overline{Q}(|y''|,|b'|) \dv{x}.
\]
The omitted director~$d$ can be reconstructed via $d=y'\times b$.

Various aspects of the outlined variational problem make the numerical
approximation of minimizers difficult. First, the function $\overline{Q}$ is convex and
continuously differentiable but not twice continuously differentiable.
Hence, to apply iterative solution techniques,
regularizations are required. Second, the nonlinear pointwise constraints
have to be discretized appropriately and incorporated in the solution
strategy. Third, the nonlinear dependence of the integrand on highest
derivatives complicates the analysis of quadrature effects. Fourth, the
iterative minimization has to be carefully done to reliably decrease the
energy and avoid the occurrence of irregular stationary configurations.

We devise a discretization of $E$ that uses piecewise cubic, $C^1$ conforming
finite elements for the approximation of $y$, and piecewise linear,
$C^0$ conforming finite elements for $b$. The variable $d$ is eliminated
via the relation $d=y'\times b$. A twice continuously differentiable
regularization $\overline{Q}_\delta$ of
$\overline{Q}$ is used and the orthogonality relations $y'\cdot b = 0$ and
$y'\cdot b'=0$ are incorporated using penalization, i.e., we consider
\[\begin{split}
E_{\delta,\varepsilon,h}[y_h,b_h] = \int_0^L & \overline{Q}_\delta(|y_h''|,|b_h'|) \dv{x} \\
&+ \frac{1}{2\varepsilon_1} \int_0^L (y_h'\cdot b_h)^2 \dv{x}
+ \frac{1}{2\varepsilon_2} \int_0^L (y_h'\cdot b_h')^2 \dv{x}.
\end{split}\]
The unit length conditions $|y'|^2=1$ and $|b|^2=1$ are imposed at
the nodes of the partitioning defining the finite element spaces, i.e.,
\[
\mathcal{I}_h^{1,0} |y_h'|^2 = \mathcal{I}_h^{1,0} |b_h|^2 = 1
\]
with the nodal interpolation operator $\mathcal{I}_h^{1,0}$ related to piecewise
linear, continuous finite element functions. We derive three convergence
results for the discretized functional and its minimizers. Under an
approximability condition on minimizers of the continuous problem we
show that these are correctly approximated by discrete minimizers.
{Assuming the density of smooth ribbon frames or a regularity property of a
minimizer} we are able to verify
the convergence of the discrete fucntionals in the sense of $\Gamma$-convergence
to the corrected Sadowsky functional.
Most generally, we deduce that the functionals $E_{\delta,\varepsilon,h}$ are
$\Gamma$-convergent to the corrected Sadowsky functional with a penalized
treatment of the orthogonality relation $y'\cdot b' =0$, i.e., for a fixed
parameter $\varepsilon_2$.  A generalization of the results based on Jensen's
inequality includes the use
of quadrature in nonquadratic terms which makes the devised method
fully practical even in the case of minimial regularity properties.
{Justifications of the assumed density and regularity results are currently
under investigation~\cite{FHMP-in-prep-pre}.}

The discretized elastic energies $E_{\delta,\varepsilon,h}$ are particularly suitable
for the minimization using gradient flows, i.e., a discrete version of the
formal evolution equations
\[
\partial_t y = -\nabla_y E[y,b] + (\lambda y')', \quad \partial_t b = - \nabla_b E[y,b] - \mu b,
\]
for given initial and boundary data and suitable definitions of the gradients with
respect to the variables $y$ and $b$. The Lagrange multipliers $\lambda$ and $\mu$
correspond to the pointwise unit length constraints. To define a discrete
variant we use a splitting of the (regularized) energy density $\overline{Q}_\delta$ into
quadratic and nonlinear parts, i.e., we have
\[
\overline{Q}_\delta(\alpha,\beta) = \frac12 \alpha^2 + \frac{5}{2} \beta^2 + \psi_\delta(\alpha^2,\beta^2).
\]
With
backward difference quotient operator~$d_t$ to approximate the time derivatives
and suitable implicit and explicit treatments of variations of the
discrete energy $E_{\delta,\varepsilon,h}$, i.e., the time stepping scheme reads
\[\begin{split}
 (d_t y_k,w)_\star + &(y_k'',w'') +  \varepsilon_1^{-1} (y_k'\cdot b_{k-1},w'\cdot b_{k-1})
+ \varepsilon_2^{-1} (y_k'\cdot b_{k-1}',w'\cdot b_{k-1}') \\
& = - (\psi_{\delta,1}'(|y_{k-1}''|^2,|b_{k-1}'|^2) y_{k-1}'',w''), \\
 (d_t b_k,r)_\dagger + & 5 (b_k',r') +  \varepsilon_1^{-1} (y_k'\cdot b_k,y_k'\cdot r)
+ \varepsilon_2^{-1} (y_k'\cdot b_k',y_k'\cdot r') \\
& = - (\psi_{\delta,2}'(|y_{k-1}''|^2,|b_{k-1}'|^2) b_{k-1}',r'),
\end{split}\]
subject to linearized pointwise constraints
\[
d_t y_k' \cdot y_{k-1}' = 0, \quad d_t b_k \cdot b_{k-1} = 0.
\]
The equations are linear and decoupled in the variables~$y_k$ and~$b_k$.
For this decoupling the penalized approximation of the orthogonality relations
using separately convex functionals is essential. We show that the
resulting discrete time-stepping scheme is energy decreasing and
that the violation of the constraints is controlled uniformly.

An alternative approach to the numerical simulation of elastic ribbons
is obtained by discretizing the two-dimensional nonlinear plate bending model
from~\cite{FrJaMu02} and the finite element method devised in~\cite{Bart13,Bart15-book}.
This has the disadvantage of discretizing a
two-dimensional domain but the advantages that the model leads to
the constrained minimization of a quadratic functional and incorporates
the width parameter which may be relevant in certain applications. However,
the explicit occurrence of the width requires the use of very
fine discretizations.

The outline of this article is as follows. We specify employed notation and
elementary facts about finite element methods in Section~\ref{sec:prelim}.
The continuous variational problem and its regularized approximation
are addressed in Section~\ref{sec:reduced}. In Section~\ref{sec:fem} we
define a discrete approximation of the regularized variational problem
and investigate its convergence to the continuous one. A discrete
gradient flow is devised and analyzed in Section~\ref{sec:iterative}.
Numerical experiments are reported in Section~\ref{sec:num_ex}.

\section{Preliminaries}\label{sec:prelim}
\subsection{Finite element spaces}
We always let $\mathcal{T}_h$ denote a partitioning $x_0<x_1<\dots<x_N$
of the interval $(0,L)$ into elements $T_j = [x_{j-1},x_j]$, $j=1,2,\dots,N$,
with  maximal mesh-size $h=\max_{j=1,\dots,N} h_{T_j}$ with
$h_{T_j}= |x_j-x_{j-1}|$. Sets of piecewise polynomial functions of
degree at most $p$
and which are $k$~times continuously differentiable are denoted
by the sets
\[
\mathcal{S}^{p,k}(\mathcal{T}_h) =
\big\{v_h \in C^k(0,L): v_h|_T \in P_p(T) \text{ for all } T\in \mathcal{T}_h\big\}.
\]
Associated with these spaces is a nodal interpolation operator
\[
\mathcal{I}_h^{p,k}: C^k(0,L)\to \mathcal{S}^{p,k}(\mathcal{T}_h).
\]
We work with the spaces $\mathcal{S}^{1,0}(\mathcal{T}_h)$ of piecewise linear, globally
continuous functions and $\mathcal{S}^{3,1}(\mathcal{T}_h)$ of piecewise cubic, globally
continuously differentiable functions for which the interpolation
operators are entirely determined by the conditions
\[
\mathcal{I}_h^{1,0}v(x_j) = v(x_j)
\]
and
\[
\mathcal{I}_h^{3,1}w (x_j) = w(x_j), \quad \big(\mathcal{I}_h^{3,1}w \big)'(x_j) = w'(x_j)
\]
for $v\in C^0([0,L])$ and $w\in C^1([0,L])$ and
$j=0,1,\dots,N$, respectively. The operators satisfy the
estimates
\[
\| v- \mathcal{I}_h^{1,0} v\|_{H^r(0,L)}  \le c_{1,0} h^{2-r} \|v''\|,
\]
for $r=0,1$, and
\[
\| w- \mathcal{I}_h^{3,1} w\|_{H^r(0,L)}  \le c_{3,1}h^{3-r} \|v^{(3)}\|,
\]
for $r=1,2$, respectively. Moreover,
$\|(\mathcal{I}_h^{3,1} w)'''\|_{L^2(T)}  \le c_{3,1}' \|w^{(3)}\|_{L^2(T)}$
for all $T\in \mathcal{T}_h$ and $w\in H^3(T)$.
We also note the inverse estimate
\begin{equation}\label{eq:inv_est}
\|v_h\|_{L^\infty(T)} \le c_{\rm inv} h_T^{-1/2} \|v_h\|_{L^2(T)}
\end{equation}
for all $T\in \mathcal{T}_h$ and $v_h\in P_p(T)$ with a constant $c_{\rm inv}$
depending on the polynomial degree~$p$. We will always assume that
the partitioning $\mathcal{T}_h$ is quasiuniform so that the minimal mesh-size
is comparable to the maximal mesh-size, i.e., $\min_{T\in\mathcal{T}_h} h_T
\ge c h$.

\subsection{Difference quotients}
Given a step-size $\tau>0$ and an arbitrary sequence $(a^k)_{k=0,\dots,K}$,
the backward
difference quotient $d_t a^k$ is for $k=1,2,\dots,K$ defined via
\[
d_t a^k = \tau^{-1} \big(a^k-a^{k-1}\big).
\]
A binomial formula implies the identity
\[
d_t a^k \cdot a^k = \frac12 d_t |a^k|^2 + \frac{\tau}{2} |d_t a^k|^2.
\]

\subsection{Lebesgue and Sobolev spaces}
Throughout this article we use standard notation for Sobolev and Lebesgue
spaces. We abbreviate the $L^2$ norm via
\[
\|\cdot \| = \|\cdot \|_{L^2(0,L)}.
\]
Throughout
we use standard embedding results, e.g., that $H^1(0,L)$ is compactly
embedded in $C^0(0,L)$.

\section{Reduced model and regularized approximation}\label{sec:reduced}
\subsection{Elastic ribbons}
We consider a dimensionally reduced model for narrow ribbons rigorously derived
in~\cite{FHMP16}. The one-dimensional variational problem is defined via the elastic
energy
\[
E[y,b] = \int_0^L \overline{Q}(|y''|,|b'|) \dv{x}.
\]
The vectors $y'$, $b$, and $d=y'\times b$ are required to constitute an orthonormal
frame, i.e.,
\[
[y',b,d]\in SO(3)
\]
almost everywhere in $(0,L)$. This implies that $y$ and $b$ satisfy the constraints
\[
|y'|^2 = 1, \quad |b|^2 =1, \quad y'\cdot b =0.
\]
Moreover, the derivation of the model leads to the constraint that the inplane
curvature component vanishes, i.e.,
\[
y'' \cdot b = 0,
\]
which in view of $y'\cdot b=0$ is equivalent to $y'\cdot b'=0$.
The function $\overline{Q}$ is given by
\[
\overline{Q}(\alpha,\beta) =
\begin{cases}
(\alpha^2 + \beta^2)^2/\alpha^2 & \mbox{if } |\alpha| \ge |\beta|, \\
4 \beta^2 & \mbox{if } |\alpha| \le |\beta|.
\end{cases}
\]
The following lemma provides structural properties of $\overline{Q}$.

\begin{lemma}[Convexity]\label{la:convex_and_diffb}
The function~$\overline{Q}:\mathbb{R}_{\ge 0}^2 \to \mathbb{R}$ is continuously differentiable and convex so
that the mapping
\[
(z,r) \mapsto \int_0^L \overline{Q}(|z|,|r|) \dv{x}
\]
is weakly lower semicontinuous on $L^2(0,L;\mathbb{R}^3)\times L^2(0,L;\mathbb{R}^3)$.
\end{lemma}

\begin{proof}
For $\alpha\ge \beta$ we have
\[
\nabla \overline{Q}(\alpha,\beta) = \begin{bmatrix} 2 \alpha - 2 \beta^4/\alpha^3 \\ 4 \beta + 4 \beta^3/\alpha^2 \end{bmatrix},
\]
and for $\alpha\le \beta$ that
\[
\nabla \overline{Q}(\alpha,\beta) = \begin{bmatrix} 0 \\ 8 \beta \end{bmatrix},
\]
so that~$\overline{Q}$ is continuously differentiable.
It is obvious that $D^2 \overline{Q}$ is positive semi-definite if $\alpha\le \beta$.
Letting $q=\beta/\alpha$ the Hessian $D^2\overline{Q}$ is for $\alpha\ge \beta$ given by
\[
D^2 \overline{Q}(\alpha,\beta) = \begin{bmatrix} 2 + 6 q^4 & - 8 q^3 \\ - 8 q^3 & 4 + 12 q^2 \end{bmatrix}.
\]
For $x=(x_1,x_2)\in \mathbb{R}^2$ we have, using $q^2 \le 1$,
\[\begin{split}
x^{\sf T} D^2 \overline{Q}(\alpha,\beta) x
&= 2 x_1^2 + 4 x_2^2 + q^2 ( 6 q^2 x_1^2 - 16 q x_1 x_2 + 12 x_2^2)  \\
&\ge q^2 (8 q^2 x_1^2 - 16 q x_1 x_2 + 16 x_2^2) \ge 0,
\end{split}\]
since, e.g., $16 q x_1 x_2 \le 4 q^2 x_1^2 + 16 x_2^2$. A Taylor expansion yields the convexity
property $\overline{Q}(\widetilde{\alpha},\widetilde{\beta}) \ge \overline{Q}(\alpha,\beta) + \nabla \overline{Q}(\alpha,\beta) \cdot ((\widetilde{\alpha},\widetilde{\beta})-(\alpha,\beta))$.
\end{proof}

\subsection{Regularization}
A smooth approximation of the function~$\overline{Q}$ is obtained by using
a regularized modulus function, i.e., we assume that for every $\delta\ge 0$ we
are given a function
\[
|\cdot|_\delta : \mathbb{R} \to \mathbb{R}_{\ge 0}
\]
that satisfies $|\cdot|_\delta \in C^2(\mathbb{R})$, obeys the estimates
\begin{equation}\label{eq:mod_reg}
|\cdot|_\delta \ge c_0 \delta, \quad |\cdot|_\delta' \le c_1, \quad |\cdot|_\delta'' \le c_2 \delta^{-1},
\end{equation}
and approximates the modulus function
uniformly, i.e., for all $x\in \mathbb{R}$ we have
\begin{equation}\label{eq:uni_reg}
\big||x|_\delta -|x|\big| \le c_{{\rm uni}} \delta.
\end{equation}
A possible choice is the function $|x|_\delta = (x^2+\delta^2)^{1/2}$.
Motivated by the relation $2 \max\{x,y\} = x + y + |x-y|$ we define
\[
\overline{Q}_\delta(\alpha,\beta) = 2\beta^2 +  \frac12 \alpha^2 + \frac12 \beta^2 + \frac12 |\alpha^2-\beta^2|_\delta +
\frac{2 \beta^4}{ \alpha^2 + \beta^2 + |\alpha^2-\beta^2|_\delta},
\]
and note that $\overline{Q}_\delta$ coincides with $\overline{Q}$ for $\delta=0$.
To prove the consistency of our discretization, that will be based on the
regularization, we note the following estimate.

\begin{lemma}[Local Lipschitz continuity]\label{la:cont_reg}
Given $\alpha,\beta,\widetilde{\alpha},\widetilde{\beta} \in \mathbb{R}$ we have
\[\begin{split}
\big|\overline{Q}(\alpha,\beta) - \overline{Q}_\delta (\widetilde{\alpha},\widetilde{\beta})\big|
 \le 21 \big(|\beta-\widetilde{\beta}||\beta+\widetilde{\beta}| + |\alpha-\widetilde{\alpha}||\alpha+\widetilde{\alpha}|+ c_{{\rm uni}} \delta\big).
\end{split}\]
\end{lemma}

\begin{proof}
We define
\[
s^2 = \alpha^2+\beta^2+|\alpha^2-\beta^2|, \quad
\widetilde{s}^2 = \widetilde{\alpha}^2+\widetilde{\beta}^2+|\widetilde{\alpha}^2-\widetilde{\beta}^2|_\delta,
\]
and note that
\[
\overline{Q}(\alpha,\beta) = 2 \beta^2 + \frac12 s^2 + \frac{2 \beta^4}{s^2}, \quad
\overline{Q}_\delta(\widetilde{\alpha},\widetilde{\beta}) = 2 \widetilde{\beta}^2 + \frac12 \widetilde{s}^2 + \frac{2 \widetilde{\beta}^4}{\widetilde{s}^2}.
\]
Straightforward calculations lead to
\[
\frac{\beta^4}{s^2} -\frac{\widetilde{\beta}^4}{\widetilde{s}^2}
= \frac{(\widetilde{s}^2-s^2)\beta^4}{s^2\widetilde{s}^2} + \frac{(\beta^2-\widetilde{\beta}^2)(\beta^2+\widetilde{\beta}^2)}{\widetilde{s}^2},
\]
and, by exchanging the roles of the variables, in case $s^2>0$,
\[
\frac{\widetilde{\beta}^4}{\widetilde{s}^2} - \frac{\beta^4}{s^2}
= \frac{(s^2-\widetilde{s}^2) \widetilde{\beta}^4}{\widetilde{s}^2 s^2} + \frac{(\widetilde{\beta}^2-\beta^2)(\widetilde{\beta}^2+\beta^2)}{s^2}.
\]
On combining the two identities, we find that
\[
\Big| \frac{\beta^4}{s^2} -\frac{\widetilde{\beta}^4}{\widetilde{s}^2} \Big|
\le
\big(|\widetilde{s}^2-s^2|+ |\beta^2-\widetilde{\beta}^2| \big)
\min \Big\{ \frac{\beta^4}{s^2\widetilde{s}^2} + \frac{\beta^2+\widetilde{\beta}^2}{\widetilde{s}^2},
\frac{\widetilde{\beta}^4}{s^2\widetilde{s}^2} + \frac{\beta^2+\widetilde{\beta}^2}{s^2}\Big\}.
\]
Noting that $s^2 \ge \beta^2$
and $\widetilde{s}^2\ge \widetilde{\beta}^2 + c_0 \delta \ge \widetilde{\beta}^2$, the estimate
implies that we have
\[
\Big| \frac{\beta^4}{s^2} -\frac{\widetilde{\beta}^4}{\widetilde{s}^2} \Big|
\le \big(|\widetilde{s}^2-s^2|+ |\beta^2-\widetilde{\beta}^2| \big)
\min \Big\{1+2 \frac{\beta^2}{\widetilde{\beta}^2+c_0\delta}, 1+ 2 \frac{\widetilde{\beta}^2}{\beta^2} \Big\}.
\]
A case distinction relating $\beta^2$ and $\widetilde{\beta}^2+c_0\delta$
shows that the minimum on the right-hand side is always bounded by~3.
Using~\eqref{eq:uni_reg}, Lipschitz continuity of the modulus function,
and binomial formulas, show that we have
\[
\big|s^2 -\widetilde{s}^2\big| \le 2 |\alpha-\widetilde{\alpha}||\alpha+\widetilde{\alpha}| + 2 |\beta-\widetilde{\beta}||\beta+\widetilde{\beta}| + c_{{\rm uni}} \delta.
\]
Combining the estimates implies that we have
\[\begin{split}
\big|\overline{Q}(\alpha,\beta) - & \overline{Q}_\delta (\widetilde{\alpha},\widetilde{\beta})\big|
\le 2 |\beta^2-\widetilde{\beta}^2| + \frac12 |s^2-\widetilde{s}^2| + 2 \Big|\frac{\beta^4}{s^2}- \frac{\widetilde{\beta}^4}{\widetilde{s}^2}\Big| \\
&\le (2+6) |\beta^2-\widetilde{\beta}^2| + \big(\frac12 + 6\big) |s^2-\widetilde{s}^2|  \\
& \le |\beta-\widetilde{\beta}||\beta+\widetilde{\beta}|\big(2+6+13 \big) + 13|\alpha-\widetilde{\alpha}||\alpha+\widetilde{\alpha}|
+ \frac{13}{2} c_{{\rm uni}} \delta.
\end{split}\]
This implies the asserted estimate.
\end{proof}

\subsection{Approximation}
Lemma~\ref{la:cont_reg} implies the uniform convergence property
$\overline{Q}_\delta \to \overline{Q}$. Another consequence is the strong approximability
of the variational problem defined by~$E$. More generally, we consider the
approximation by variational problems with energy functionals~$E_{\delta,\varepsilon_2}$
using the regularized function $\overline{Q}_\delta$
and a penalized treatment of the relation $y'\cdot b'=-y''\cdot b =0$,
i.e.,
\[
E_{\delta,\varepsilon_2}[y,b] = \int_0^L \overline{Q}_\delta(|y''|,|b'|) \dv{x}
+ \frac{1}{2\varepsilon_2} \int_0^L (y'\cdot b')^2 \dv{x}.
\]
We incorporate boundary conditions via a bounded linear operator $L_{\rm bc}$,
which is assumed to depend only on the boundary values of $y$, $y'$, and $b$.
It defines the set of admissible frames
\[\begin{split}
\mathcal{A}_{\rm frame} = \big\{  (y,b) &\in H^2(0,L;\mathbb{R}^3)\times H^1(0,L;\mathbb{R}^3) : \\
& [y',b,y'\times b]\in SO(3) \ \text{a.e. in } (0,L), \  L_{\rm bc}[y,b] = \ell_{\rm bc} \big\}.
\end{split}\]
Note that the relation $y'\cdot b'=0$ is not included in the set~$\mathcal{A}_{\rm frame}$.
The variational problem under consideration then reads as follows.
\begin{equation}\label{eq:cont_prob}
\tag{$\textrm{P}_{\delta,\varepsilon_2}$}
\left\{
\begin{array}{l}
\text{Given $\varepsilon_2,\delta\ge 0$ find a minimizing pair} \\
\text{$(y,b) \in \mathcal{A}_{\rm frame}$ for } E_{\delta,\varepsilon_2}[y,b].
\end{array}
\right.
\end{equation}
The original problem formally corresponds to $\delta=0$ and $\varepsilon_2=0$ when the
penalty term is interpreted as a strict constraint. Due to the convexity of
$\overline{Q}$ a solution exists for $\delta=0$ and $\varepsilon_2 \ge 0$.
Our numerical scheme will be a discrete variant of the following approximation result.

\begin{proposition}[Strong approximability]\label{prop:strong_approx}
For every $(y,b) \in \mathcal{A}_{\rm frame}$ with $y'\cdot b' = 0$
and every $\sigma>0$ there exists
\[
(y_\sigma,b_\sigma)\in \big[C^\infty(0,L;\mathbb{R}^3)\big]^2 \cap \mathcal{A}_{\rm frame}
\]
with
\[
\|y-y_\sigma\|_{H^2(0,L)} + \|b-b_\sigma\|_{H^1(0,L)}
+ \|y_\sigma'\cdot b_\sigma'\| \le \sigma,
\]
such that
\[
\big| E[y,b] - E_{\delta,\varepsilon_2}[y_\sigma,b_\sigma]\big| \le
c_{{\rm sa}} \big(\sigma^2\varepsilon_2^{-1} +\delta+\sigma\big).
\]
In particular, the functionals $\big(E_{\delta,\varepsilon_2}\big)_{\delta,\varepsilon_2}$ are
$\Gamma$-convergent to $E$ for $(\delta,\varepsilon_2)\to 0$ and to $E_{0,\varepsilon_2}$ for
fixed $\varepsilon_2>0$ and $\delta \to 0$.
\end{proposition}

\begin{proof}
A regularized frame with the asserted properties
can be constructed using mollification and correction
techniques,~cf.~\cite{BarRei19-in-prep-pre}.
Lemma~\ref{la:cont_reg} and the identity
$y_\sigma'\cdot b_\sigma' = (y_\sigma-y)'\cdot b_\sigma'+ y'\cdot (b_\sigma-b)'$
then imply the approximation result. The $\Gamma$-convergence results
are immediate consequences of the estimates.
\end{proof}

\begin{remark}\label{rem:explode}
We remark that the quantity
\[
[y,b]_\sigma = \|y_\sigma\|_{H^3(0,L)} + \|b_\sigma\|_{H^2(0,L)}
\]
is in general unbounded as $\sigma\to 0$.
If $(y,b) \in H^3(0,L;\mathbb{R}^3)\times H^2(0,L;\mathbb{R}^3)$ then it remains bounded.
\end{remark}

\section{Finite element approximation}\label{sec:fem}
We write the regularized minimization problem with a penalized treatment of
the orthogonality relations as
\[\begin{split}
E_{\delta,\varepsilon}[y,b] &= \int_0^L \overline{Q}_\delta(|y''|,|b'|) \dv{x} + P_\varepsilon[y,b] \\
&= \frac12 \int_0^L |y''|^2 + 5 |b'|^2 + \psi_\delta(|y''|^2,|b'|^2) \dv{x} + P_\varepsilon[y,b],
\end{split}\]
where for $\varepsilon = (\varepsilon_1,\varepsilon_2)$ we define
\[
P_\varepsilon[y,b] = \frac{1}{2\varepsilon_1} \int_0^L (y'\cdot b)^2 \dv{x}
+ \frac{1}{2\varepsilon_2} \int_0^L (y'\cdot b')^2 \dv{x},
\]
and the function $\psi_\delta$ is defined via
\[
\psi_\delta(\alpha^2,\beta^2) = |\alpha^2-\beta^2|_\delta + \frac{4 \beta^4}{\alpha^2+\beta^2 + |\alpha^2-\beta^2|_\delta}.
\]
The minimization of $E_{\delta,\varepsilon}$ is subject to the pointwise unit length conditions
\[
|y'|^2 = 1, \quad |b|^2 = 1.
\]
The corresponding discrete problem uses the conforming subspaces
\[\begin{split}
V_h & = \mathcal{S}^{3,1}(\mathcal{T}_h)^3 \times \mathcal{S}^{1,0}(\mathcal{T}_h)^3 \\
&\subset V = H^2(0,L;\mathbb{R}^3)\times H^1(0,L;\mathbb{R}^3).
\end{split}\]
By imposing the  nonlinear constraints at the nodes of the partitioning $\mathcal{T}_h$ we
arrive at the following discrete problem.
\begin{equation}\label{eq:discr_prob}
\tag{$\textrm{P}_{\delta,\varepsilon,h}$}
\left\{ \begin{array}{l}
\text{Find a minimizing pair $(y_h,b_h) \in V_h$ for} \\[1.5mm]
\displaystyle{
E_{\delta,\varepsilon,h}[y_h,b_h]
= \frac12 \int_0^L |y_h''|^2 + 5 |b_h'|^2 + \psi_\delta(|y''|^2,|b'|^2) \dv{x}} \\
\displaystyle{ \hspace*{3.5cm} + P_\varepsilon[y_h,b_h]} \\[3mm]
\text{subject to}
\displaystyle{\quad \mathcal{I}_h^{1,0} |y_h'|^2 = \mathcal{I}_h^{1,0} |b_h|^2 = 1 \text{ and } L_{\rm bc}[y_h,b_h] = \ell_{\rm bc}.}
\end{array}\right.
\end{equation}
The following proposition implies that discrete minimizers converge to
minimizers of the continuous problem provided that these satisfy an
approximability condition. Because of the assumption that the
operator $L_{\rm bc}$ only depends on the boundary values of $y$, $y'$, and $b$,
these are preserved
under the nodal interopolation operators $\mathcal{I}_h^{3,1}$ and $\mathcal{I}_h^{1,0}$.

\begin{proposition}[Partial $\Gamma$-convergence]\label{prop:gamma_conv}
The discrete
constrained energy functionals $E_{\varepsilon,\delta,h}$ approximate the
continuous constrained energy functional $E$ in the sense of the following
assertions. \\
(i) Let $(y,b)\in V$. If $(\sigma,\varepsilon_1,\varepsilon_2,\delta) \to 0$ as $h\to 0$
with $\sigma>0$ as in Proposition~\ref{prop:strong_approx} such that
\[
h^2 \varepsilon_1^{-1} + \sigma^2 \varepsilon_2^{-1} + h^2 [y,b]_\sigma^2 \varepsilon_2^{-1} \to 0
\]
then there exists a sequence
of discrete admissible pairs $(y_h,b_h) \in V_h$ with $y_h \to y$
in $H^2$ and $b_h\to b$ in $H^1$ as $h\to 0$ and such that
\[
E_{\delta,\varepsilon,h}[y_h,b_h] \to E[y,b].
\]
(ii) For every bounded sequence $(y_h,b_h)\in V_h$ of admissible
pairs with weak limit $(y,b) \in V$ we have that the limitting pair
is admissible and that
\[
E[y,b] \le \liminf_{(\varepsilon,\delta,h)\to 0} E_{\delta,\varepsilon,h}[y_h,b_h].
\]
\end{proposition}

\begin{proof}
(i) We use the result of Proposition~\ref{prop:strong_approx} and consider
for $\sigma>0$ the smooth frame $(y_\sigma,b_\sigma)\in \mathcal{A}_{\rm frame}$ that
approximates the general frame $(y,b) \in \mathcal{A}_{\rm frame}$. We then
consider for $h>0$ the pair $(y_h,b_h)$ defined via
\[
y_h = \mathcal{I}_h^{3,1} y_\sigma, \quad b_h = \mathcal{I}_h^{1,0} b_\sigma.
\]
The pair $(y_h,b_h)$ is admissible in the discrete minimization problem
 i.e., the functions satisfy the boundary conditions and the nodal unit
length constraints. Since $y_\sigma'\cdot b_\sigma= 0$ we find that
\[
y_h' \cdot b_h = (y_h'-y_\sigma')\cdot b_h + y_\sigma' \cdot (b_h-b_\sigma),
\]
and hence, we have
\[
\|y_h' \cdot b_h\| \le c h.
\]
Similarly, with the quantity $c_\sigma= [y,b]_\sigma$ from Remark~\ref{rem:explode},
that bounds the norms
$\|y_\sigma\|_{H^3(0,L)}$ and $\|b_\sigma\|_{H^2(0,L)}$, we have that
\[
\|y_h'\cdot b_h'\| \le \sigma + c h c_\sigma.
\]
{The repeated application of Lemma~\ref{la:cont_reg} provides the estimate
\[\begin{split}
\big| & \overline{Q}_\delta(\alpha,\beta) - \overline{Q}_\delta(\widetilde{a},\widetilde{b})\big| \\
&\le \big|\overline{Q}_\delta(\alpha,\beta) - \overline{Q}(\alpha,\beta)\big|
+ \big|\overline{Q}(\alpha,\beta) - \overline{Q}(\widetilde{a},\widetilde{b})\big|
+\big|\overline{Q}_\delta(\widetilde{a},\widetilde{b}) - \overline{Q}_\delta(\widetilde{a},\widetilde{b})\big|  \\
&\le  c \big(|\beta-\widetilde{\beta}||\beta+\widetilde{\beta}| + |\alpha-\widetilde{\alpha}||\alpha+\widetilde{\alpha}| + \delta\big).
\end{split}\]
In combination with Proposition~\ref{prop:strong_approx} this yields that}
\[\begin{split}
\big|E[y,b]& -E_{\delta,\varepsilon,h}[y_h,b_h] \big| \\
&\le \big|E[y,b] - E_{\delta,\varepsilon_2}[y_\sigma,b_\sigma] \big|
  + \big|E_{\delta,\varepsilon_2}[y_\sigma,b_\sigma]-E_{\delta,\varepsilon,h}[y_h,b_h]  \big| \\
&\le c \big(\sigma^2 \varepsilon_2^{-1} + \delta  + \sigma\big)
+ c \big(h c_\sigma + \delta + h^2 \varepsilon_1^{-1} + (\sigma^2+ h^2 c_\sigma^2)\varepsilon_2^{-1} \big).
\end{split}\]
The conditions of the proposition imply that the right-hand tends
to zero as $h\to 0$.  \\
(ii) Assume that $(y_h,b_h)_{h>0}$ is a sequence of admissible discrete
pairs with weak limit
$(y,b)\in V$. We may assume
that the sequence $E_{\delta,\varepsilon,h}[y_h,b_h]$ is bounded. Compactness properties
of embeddings imply that $|y'|^2= 1$ and $|b|^2 =1$
in $(0,L)$. Since $P_\varepsilon[y_h,b_h]$ is bounded we have
that $y_h'\cdot b_h$ and $y_h'\cdot b_h'$
converge to zero in $L^2(0,L)$. Noting that $y_h'\to y$
in $L^\infty(0,L;\mathbb{R}^3)$ it follows that $y'\cdot b=0$
and $y'\cdot b' = 0$. Since $L_{\rm bc}$ is bounded and linear it is weakly continuous
and we verify
that $L_{\rm bc}[y,b] = \ell_{\rm bc}$. Hence, we conclude that $(y,b) \in \mathcal{A}_{\rm frame}$.
Since the functional~$E$ is weakly lower semicontinuous, since the functionals
$P_\varepsilon$ are
nonnegative, and since the regularization of the integrand $\overline{Q}$ is uniformly
controlled by Lemma~\ref{la:cont_reg}, we find that
\[\begin{split}
E[y,b] & \le \liminf_{h\to 0} E[y_h,b_h] \\
&\le \liminf_{(\varepsilon,h,\delta)\to 0} E_{\delta,\varepsilon,h}[y_h,b_h] + c \delta \\
&= \liminf_{(\varepsilon,h,\delta)\to 0} E_{\delta,\varepsilon,h}[y_h,b_h].
\end{split}\]
This proves the proposition.
\end{proof}

\begin{remark}
Let $\sigma=h^\alpha$, $\varepsilon_j = h^{\beta_j}$, $j=1,2$, and $\delta = h^\gamma$
with $\alpha,\beta_1,\beta_2,\gamma>0$. Then the condition of the proposition reads
\[
h^{2-\beta_1} + h^{2\alpha-\beta_2} + h^{2-\beta_2}  [y,b]_{h^\alpha}^2 \to 0.
\]
Assuming that $[y,b]_\sigma \le c_{{\rm appr}} \sigma^{-p}$ for some $p\ge 0$
then this condition becomes
\[
h^{2-\beta_1} + h^{2\alpha-\beta_2} + h^{2-\beta_2 -2 p \alpha} \to 0.
\]
Hence, for $\beta_1<2$ and $\beta_2$ sufficiently small the condition is
satisfied.  If $y\in H^3(0,L;\mathbb{R}^3)$ and
$b\in H^2(0,L;\mathbb{R}^3)$ then we have $p=0$ and the condition is satisfied
for $\beta_2 < \min\{ 2,2\alpha\}$.
\end{remark}

Proposition~\ref{prop:gamma_conv} implies the convergence of discrete
almost-minimizers under certain conditions.

\begin{corollary}[Approximation of regular minimizers]\label{cor:conv_mins}
Let $(y_h,b_h)_{h>0}$ be a sequence of discrete minimizers for the
problems~\eqref{eq:discr_prob} possibly up to
tolerances $\varepsilon_{{\rm tol}}$ that vanish as $h\to 0$. If
an accumulation point $(y,b)\in V$ of the sequence is such that
there exists a sequence $\sigma\to 0$ such that
\[
h^2 \varepsilon_1^{-1} + \sigma^2 \varepsilon_2^{-1} + h^2 [y,b]_\sigma^2 \varepsilon_2^{-1} \to 0
\]
as $h\to 0$ then every weak accumulation point
$(y,b)\in V$ as $(\delta,\varepsilon,h)\to 0$ minimizes the constrained energy
functional~$E$.
\end{corollary}

\begin{proof}
Since the functionals $E_{\delta,\varepsilon,h}$ are equi-coercive a sequence
of almost-minimizers is bounded and the result is a direct consequence
of Proposition~\ref{prop:gamma_conv}.
\end{proof}

If the regularization result of Proposition~\ref{prop:strong_approx}
can be refined such that $y_\sigma'\cdot b_\sigma'= 0$ for every $\sigma>0$ then
a $\Gamma$-convergence result can be deduced.

\begin{corollary}[$\Gamma$-convergence to $E$]\label{cor:gamma_dense}
Assume that for every $\sigma>0$ there exists a pair $(y_\sigma,b_\sigma)$
as in Proposition~\ref{prop:strong_approx} with the additional property
that $y_\sigma'\cdot b_\sigma' =0$ in $(0,L)$. If $\varepsilon_1>0$ and
$\varepsilon_2>0$ are chosen such that
\[
h^2 \varepsilon_1^{-1} + h^2 \varepsilon_2^{-1} \to 0
\]
as $h\to 0$ then the constrained functionals $(E_{\delta,\varepsilon,h})$
are $\Gamma$-convergent to the constrained functional~$E$.
\end{corollary}

\begin{proof}
The functionals $E_{\delta,\varepsilon,h}$ are equi-coercive and the condition
of Proposition~\ref{prop:gamma_conv} reduces to the condition
$h^2 \varepsilon_1^{-1} + h^2 [y,b]_\sigma^2 \varepsilon_2^{-1} \to 0$ as $h\to 0$.
\end{proof}

A general approximation result can be established if
the penalty parameter $\varepsilon_2>0$ is fixed so that the corrected
Sadowsky functional with penalized treatment of the orthogonality
$y'\cdot b'=0$ is approximated.

\begin{corollary}[$\Gamma$-convergence to a penalized approximation]\label{cor:gamma_pen}
Let $\varepsilon_2>0$ be fixed. Then the constrained functionals
$(E_{\delta,\varepsilon,h})$ are $\Gamma$-convergent to the constrained
functional $E_{0,\varepsilon_2}$ as $(\delta,\varepsilon_1,h)\to 0$ provided that
$h^2 \varepsilon_1^{-1} \to 0$.
\end{corollary}

\begin{proof}
The result is an immediate consequence of the proof of
Proposition~\ref{prop:gamma_conv} and the equi-coercivity
of the functionals $E_{\delta,\varepsilon,h}$.
\end{proof}

The numerical realization of the functionals $E_{\delta,\varepsilon,h}$ requires
the use of quadrature. For an element $T=[x_1,x_2]$ and a
function $r\in C^1(T)$ we use averaging operators
\[\begin{split}
\mathcal{A}_h r'|_T &= |T|^{-1} \int_T r' \dv{x} = h_T^{-1} \big(r(x_2)-r(x_1)\big), \\
\mathcal{M}_h r|_T &= |T|^{-1} \int_T \mathcal{I}_h^{1,0} r \dv{x} =  \frac12 \big( r(x_1) + r(x_2)\big).
\end{split}\]
We note that we have the estimates
\[\begin{split}
\|y_h'- \mathcal{M}_h y_h'\|_{L^\infty(0,L)} &\le c h^{1/2} \|y_h''\|, \\
\|\mathcal{A}_h y_h'' - y_h''\| & \le c h \|y_h'''\|.
\end{split}\]
Moreover, the operator $\mathcal{A}_h:L^2(0,L)\to L^2(0,L)$
is bounded and self-adjoint; as a consequence of Jensen's or H\"older's
inequality we have
\begin{equation}\label{eq:jensen}
\|\mathcal{A}_h y_h'' \| \le \|y_h''\|.
\end{equation}
With the help of the operators $\mathcal{A}_h$ and $\mathcal{M}_h$ we can define a
discrete functional which allows for exact numerical integration and
which $\Gamma$-converges to the continuous variational problem under slightly
more restrictive conditions on the parameters.

\begin{proposition}[Quadrature]\label{prop:quadrature}
The results of Proposition~\ref{prop:gamma_conv} and {Corollaries~\ref{cor:conv_mins},
\ref{cor:gamma_dense}, and~\ref{cor:gamma_pen}}
remain valid if $E_{\delta,\varepsilon,h}$ is replaced by the functional
\[\begin{split}
E_{\delta,\varepsilon,h}^{{\rm quad}}[y_h,b_h] & =
\frac12 \int_0^L |y_h''|^2 + 5 |b_h'|^2 + \psi_\delta(|\mathcal{A}_h y_h''|^2,|b_h'|^2) \dv{x} \\
& \qquad + \frac{1}{2\varepsilon_1} \int_0^L \mathcal{I}_h^{1,0} \big[(y_h'\cdot b_h)^2\big] \dv{x}
 + \frac{1}{2\varepsilon_2} \int_0^L  (\mathcal{M}_h y_h'\cdot b_h')^2 \dv{x},
\end{split}\]
i.e., it approximates $E$ in the sense of Proposition~\ref{prop:gamma_conv}
if the additional condition $h\varepsilon_1^{-1} + h^{1/2} \varepsilon_2^{-1/2}$ is
satisfied.
\end{proposition}

\begin{proof}
We discuss necessary changes in the proof of Proposition~\ref{prop:gamma_conv}. \\
The penalty terms involving the operators $\mathcal{I}_h^{1,0}$ and $\mathcal{M}_h$
can be replaced by the penalty terms without these operators leading to
error terms that vanish as $h\to 0$,
i.e., for a bounded sequence $(y_h,b_h)_{h>0}$
in $H^2(0,L;\mathbb{R}^3)\times H^1(0,L;\mathbb{R}^3)$ we have for the first penalty term that
\[\begin{split}
\Big|\frac{1}{\varepsilon_1}\int_0^L   \mathcal{I}_h^{1,0} \big[(y_h'\cdot b_h)^2\big] -  & (y_h'\cdot b_h)^2 \dv{x} \Big|
\le \varepsilon_1^{-1} \|(y_h'\cdot b_h)^2 - \mathcal{I}_h^{1,0} (y_h' \cdot b_h)^2\|_{L^1(0,L)}  \\
&\le c h \varepsilon_1^{-1} \|(y_h'\cdot b_h) (y_h''\cdot b_h + y_h'\cdot b_h')\|_{L^1(0,L)} \\
&\le c h \varepsilon_1^{-1} \|y_h'\cdot b_h\| \big(\|y_h''\| \|b_h\|_{L^\infty(0,L)}
  + \|y_h'\|_{L^\infty(0,L)} \|b_h'\|\big) \\
&\le c h \varepsilon_1^{-1},
\end{split}\]
while for the second penalty term we have
\[\begin{split}
\Big|\frac{1}{\varepsilon_2}\int_0^L  & (\mathcal{M}_h y_h'\cdot b_h')^2 - (y_h'\cdot b_h')^2 \dv{x} \Big| \\
&\le \varepsilon_2^{-1} \|\mathcal{M}_h y_h' - y_h'\|_{L^\infty(0,L)} \|b_h'\|^2 \|\mathcal{M}_h y_h' + y_h'\|_{L^\infty(0,L)}
\le c h^{1/2} \varepsilon_2^{-1}.
\end{split}\]
(i) For the attainment result we note that
the Lipschitz continuity of $\psi_\delta$ (see Lemma~\ref{la:props_psi} below)
yields the estimate
\[\begin{split}
\int_0^L & \Big|\psi_\delta(|\mathcal{A}_h y_h''|^2,|b_h'|^2) - \psi_\delta(|y_h''|^2,|b_h'|^2)\Big| \dv{x} \\
&\le c_\psi \int_0^L |\mathcal{A}_h y_h'' - y_h''||\mathcal{A}_h y_h''+ y_h''| \dv{x}
\le c h \|y_h'''\| \le c h c_\sigma,
\end{split}\]
where we used that that the interpolation operator $\mathcal{I}_h^{3,1}$ is elementwise
stable in~$H^3$. Hence, we have that
\[
\big|E_{\delta,\varepsilon,h}[y_h,b_h]-E_{\delta,\varepsilon,h}^{{\rm quad}}[y_h,b_h]\big| \to 0.
\]
(ii) To establish the lower bound property we note that if
 $y_h \rightharpoonup y$ in $H^2(0,L;\mathbb{R}^3)$ and $b_h\rightharpoonup b$ in $H^1(0,L;\mathbb{R}^3)$
then~\eqref{eq:jensen} and the weak lower semicontinuity of $E$ identified
in Lemma~\ref{la:convex_and_diffb} show that
\[\begin{split}
\liminf_{h,\delta \to 0} \frac12 & \int_0^L |y_h''|^2 + 5 |b_h'|^2 + \psi_\delta(|\mathcal{A}_h y_h''|^2,|b_h'|^2) \dv{x} \\
&\ge \liminf_{h,\delta \to 0} \int_0^L |\mathcal{A}_h y_h''|^2 + 5 |b_h'|^2 + \psi_\delta(|\mathcal{A}_h y_h''|^2,|b_h'|^2) \dv{x} \\
&\ge \liminf_{h,\delta \to 0} \int_0^L |\mathcal{A}_h y_h''|^2 + 5 |b_h'|^2 + \psi(|\mathcal{A}_h y_h''|^2,|b_h'|^2) \dv{x}  - c \delta \\
&= \liminf_{h \to 0} \int_0^L \overline{Q} (|\mathcal{A}_h y_h''|,|b_h'|) \dv{x} \\
&\ge \int_0^L |y''|^2 + 5 |b'|^2 + \psi(|y''|^2,|b'|^2) \dv{x}.
\end{split}\]
This implies that $E[y,b] \le \liminf_{(\varepsilon,\delta,h)\to 0}
E_{\delta,\varepsilon,h}^{{\rm quad}}[y_h,b_h]$ and proves the assertion.
\end{proof}

\begin{remark}
An inspection of the proof of Proposition~\ref{prop:gamma_conv}
implies that no condition on $\varepsilon_1$ is required if quadrature
is used. For the recovery sequence $(y_h,b_h)$ constructed in~(i)
we have that the nodal values of $y_h'\cdot b_h$ and hence the
penalty term vanish. For a sequence $(y_h,b_h)$ in~(ii) with bounded
energies we have $\|\mathcal{I}_h[(y_h'\cdot b_h)^2]\|_{L^1(0,L)} \le c \varepsilon$
which implies that $y'\cdot b=0$ for accumulation points.
\end{remark}

\section{Iterative minimization}\label{sec:iterative}

\subsection{Discrete gradient flow}
We devise a decoupled gradient scheme with a linearized treatment
of the pointwise constraints to compute stationary configurations.
For ease of readability we often omit the spatial discretization
parameter $h$ in what follows and denote
\[
Y_h = \mathcal{S}^{3,1}(\mathcal{T}_h)^3, \quad B_h = \mathcal{S}^{1,0}(\mathcal{T}_h)^3.
\]
The variational derivatives of the (discretized) penalty
functional $P_{\varepsilon,h}$ with respect to $y$ and $b$
are denoted by $\delta_y P_\varepsilon$ and $\delta_b P_\varepsilon$, respectively.
For the function
$\psi_\delta$ the partial derivatives are denoted by $\psi_{\delta,1}' $ and
$\psi_{\delta,2}'$. We assume that the boundary conditions are separated,
i.e., that we have
\[
L_{\rm bc}[y,b] = \big(L_{\rm bc}^1[y],L_{\rm bc}^2[b]\big),
\]
and $\ell_{\rm bc} = (\ell_{\rm bc}^1,\ell_{\rm bc}^2)$.

\begin{algorithm}[Gradient descent]\label{alg:grad_desc}
Choose $\tau>0$ and $(y_0,b_0)\in Y_h\times B_h$ with
\[
L_{\rm bc}^1[y_0] = \ell_{\rm bc}^1, \quad L_{\rm bc}^2[b_0] = \ell_{\rm bc}^2
\]
and the nodal length conditions
\[
|y_0'|^2=1, \quad |b_0|^2 = 1;
\]
set $k=1$. \\
(1) Compute $y_k \in Y_h$ such that
\[\begin{split}
 (d_t y_k,w)_\star + (y_k'',w'') +  \delta_y P_\varepsilon[y_k&,b_{k-1};w] = \\
& - (\psi_{\delta,1}'(|\mathcal{A} y_{k-1}''|^2,|b_{k-1}'|^2) \mathcal{A} y_{k-1}'',\mathcal{A} w'')
\end{split}\]
for all $w\in Y_h$ subject to the nodal constraints
\[
d_t y_k'\cdot y_{k-1}' = 0, \quad w'\cdot y_{k-1}' = 0
\]
and the boundary conditions
\[
L_{\rm bc}^1[d_t y_k] = 0, \quad L_{\rm bc}^1[w] = 0.
\]
(2) Compute  $b_k \in B_h$ such that
\[\begin{split}
 (d_t b_k,r)_\dagger + 5 (b_k',r') + \delta_b P_\varepsilon[y_k&,b_k;r] = \\
& - (\psi_{\delta,2}'(|\mathcal{A} y_{k-1}''|^2,|b_{k-1}'|^2) b_{k-1}',r')
\end{split}\]
for all $r\in B_h$ subject to the nodal constraints
\[
d_t b_k \cdot b_{k-1} = 0, \quad r \cdot b_{k-1} = 0
\]
and the boundary conditions
\[
L_{\rm bc}^2[d_t b_k] = 0, \quad L_{\rm bc}^2[r] = 0.
\]
(3) Stop the iteration if
$\|d_t y_k\|_\star + \|d_t b_k\|_\dagger \le \varepsilon_{{\rm stop}}$;
set $k\to k+1$ and continue with~(1) otherwise.
\end{algorithm}

We remark that all steps of Algorithm~\ref{alg:grad_desc} define linear
systems of equations that can be solved efficiently.

\subsection{Linearization estimates}
The iterative numerical treatment of the model problem requires the use of
linearizations of nonlinearities occurring in the energy functional.
The following lemmas provide estimates on derivatives of the
nonquadratic part of the regularized function~$\overline{Q}_\delta$.

\begin{lemma}[Derivative bounds]\label{la:props_psi}
The function $\psi_\delta$, for $r,s\in \mathbb{R}_{\ge 0}$ defined via
\[
\psi_\delta(r,s) =  |r-s|_\delta + \frac{4 s^2}{r+s+|r-s|_\delta},
\]
satisfies
\[
|\nabla \psi_\delta| \le c_{\psi,1}, \quad |D^2 \psi_\delta|\le c_{\psi,2} \delta^{-1}.
\]
\end{lemma}

\begin{proof}
Let $\sigma_\delta = |\cdot|_\delta'$ and $D_\delta = \sigma_\delta'= |\cdot|_\delta''$. We have
\[\begin{split}
\partial_1 \psi_\delta (r,s) &= \sigma_\delta(r-s) - \frac{4s^2}{(r+s+|r-s|_\delta)^2}(1+\sigma_\delta(r-s)) \\
\partial_2 \psi_\delta (r,s) &= - \sigma_\delta(r-s) + \frac{8s}{r+s+|r-s|_\delta}
- \frac{4 s^2}{(r+s+|r-s|_\delta)^2}(1-\sigma_\delta(r-s)).
\end{split}\]
The formulas imply the uniform boundedness of the partial derivatives.
For higher order partial derivatives we have, e.g.,
\begin{multline*}
\partial_1^2 \psi_\delta (r,s) = D_\delta(r-s) + \frac{8s^2}{(r+s+ |r-s|_\delta)^3}\big(1+\sigma_\delta(r-s)\big)^2 \\
  - \frac{4s^2}{(r+s+|r-s|_\delta)^2}D_\delta(r-s).
\end{multline*}
The remaining second order derivatives lead
to similar formulas. Incorporating the assumptions on the
second derivative of the regularized modulus function and using
a case distinction relating $r+s$ to $\delta$ imply the bounds
on the Hessian of $\psi_\delta$.
\end{proof}

A Taylor formula is needed to control the explicit treatment of the nonlinear
part~$\psi_\delta$ in Algorithm~\ref{alg:grad_desc}.

\begin{lemma}[Taylor formula]\label{la:taylor_psi}
For all $a,b,\widetilde{a},\widetilde{b} \in \mathbb{R}^3$ there exist $\widehat{a},\widehat{b} \in \mathbb{R}^3$ such that
\[\begin{split}
\frac12 \big(\psi_\delta(|a|^2,|b|^2) - \psi_\delta(|\widetilde{a}|^2,|\widetilde{b}|^2)& \big)
- \nabla \psi_\delta(|\widetilde{a}|^2,|\widetilde{b}|^2) \cdot
\begin{bmatrix} \widetilde{a} \cdot (a-\widetilde{a})  \\ \widetilde{b} \cdot (b-\widetilde{b}) \end{bmatrix} \\
& =  \Gamma(\widehat{a},\widehat{b},a,b,\widetilde{a},\widetilde{b}),
\end{split}\]
with a function $\Gamma$ that satisfies
\[
\big|\Gamma(\widehat{a},\widehat{b},a,b,\widetilde{a},\widetilde{b})\big| \le
c_\Gamma \big[\delta^{-1} \big(|a|^2 + |\widetilde{a}|^2 + |b|^2 + |\widetilde{b}|^2\big) + 1\big]
\big(|a-\widetilde{a}|^2 + |b-\widetilde{b}|^2\big).
\]
\end{lemma}

\begin{proof}
Defining $g(a,b) = (1/2) \psi_\delta(|a|^2,|b|^2)$ we have
\[
g(a,b) -  g(\widetilde{a},\widetilde{b}) - D g(\widetilde{a},\widetilde{b}) \begin{bmatrix} a-\widetilde{a} \\ b-\widetilde{b} \end{bmatrix}
=  \frac12 \begin{bmatrix} a-\widetilde{a} \\ b-\widetilde{b} \end{bmatrix}^{\sf T}
 D^2g(\widehat{a},\widehat{b})  \begin{bmatrix} a-\widetilde{a} \\ b-\widetilde{b} \end{bmatrix},
\]
where $(\widehat{a},\widehat{b})$ belongs to the line segment connecting $(a,b)$ with
$(\widetilde{a},\widetilde{b})$. We have
\[
\partial_1 g(\widehat{a},\widehat{b}) = \partial_1 \psi_\delta(|\widehat{a}|^2,|\widehat{b}|^2) \widehat{a}, \quad
\partial_2 g(\widehat{a},\widehat{b}) = \partial_2 \psi_\delta(|\widehat{a}|^2,|\widehat{b}|^2) \widehat{b},
\]
and consequently,
\[\begin{split}
\partial_1^2 g(\widehat{a},\widehat{b}) &= 2\partial_1^2 \psi_\delta(|\widehat{a}|^2,|\widehat{b}|^2) \widehat{a} \widehat{a}^{\sf T}  + \partial_1 \psi_\delta(|\widehat{a}|^2,|\widehat{b}|^2) I_{3\times 3}, \\
\partial_2^2 g(\widehat{a},\widehat{b}) &= 2 \partial_2^2 \psi_\delta(|\widehat{a}|^2,|\widehat{b}|^2) \widehat{b} \widehat{b}^{\sf T} + \partial_2 \psi_\delta(|\widehat{a}|^2,|\widehat{b}|^2) I_{3\times 3}, \\
\partial_1 \partial_2 g(\widehat{a},\widehat{b}) &= 2\partial_1 \partial_2 \psi_\delta(|\widehat{a}|^2,|\widehat{b}|^2) \widehat{b} \, \widehat{a}^{\sf T}.
\end{split}\]
The bounds of Lemma~\ref{la:props_psi} imply the asserted result.
\end{proof}

\subsection{Energy monotonicity}
Stability of Algorithm~\ref{alg:grad_desc} can be established
if the step size~$\tau$
satisfies a smallness condition defined by the regularization
parameter~$\delta$ and the mesh-size~$h$. The following proposition
implies that the iteration approximates stationary points for
the discrete energy functional.

\begin{proposition}[Energy decay]\label{prop:ener_decay}
Assume that there exist $c_\star,c_\dagger >0$ such that for all
$w\in Y_h$ and $r \in B_h$ we have
\[
\|w''\| \le c_\star \|w\|_\star, \quad \|r'\| \le c_\dagger \|r\|_\dagger.
\]
Then there exists a constant $c_{{\rm ed}}>0$ such that
for the uniquely defined iterates $(y_k,b_k)_{k=0,1,\dots}$ of
Algorithm~\ref{alg:grad_desc} we have
\[
E_{\delta,\varepsilon,h}[y_K,b_K]
+ \big(1-c_{{\rm ed}} \tau \delta^{-1}h^{-1}\big)
\tau \sum_{k=1}^K \big(\|d_t y_k\|_\star^2 + \|d_t b_k\|_\dagger^2 \big)
\le E_{\delta,\varepsilon,h}[y_0,b_0].
\]
Moreover, if $2 c_{{\rm ed}} \tau \le \delta h$ we have that
\[
\|\mathcal{I}_h^{1,0} |b_k|^2-1\|_{L^\infty(0,L)} + \|\mathcal{I}_h^{1,0}|y_k'|^2 -1\|_{L^\infty(0,L)}
\le  c_{{\rm ul}} \tau e_0,
\]
with $e_0 = E_{\delta,\varepsilon,h}[y_0,b_0]$.
\end{proposition}

\begin{proof}
We adopt an inductive argument and first prove an intermediate bound which
will allow us to control the explicit treatment of the nonlinearities related
to $\psi$. The conditional energy decay estimate then implies the control of
the violation of the pointwise constraints. Well-posedness of the iteration
is an immediate consequence of the Lax--Milgram lemma.  \\
{\em Step 1.} We choose $w=d_t y_k$ and $r= d_t b_k$ in Steps~(1) and~(2) of
Algorithm~\ref{alg:grad_desc}. Incorporating~\eqref{eq:jensen} this implies that
\begin{multline}\label{eq:rough_est}
\|d_t y_k\|_\star^2 + \|d_t b_k\|_\dagger^2 +
d_t \Big(\frac12 \|y_k''\|^2 + \frac52 \|b_k'\|^2 + P_{\varepsilon,h}[y_k,b_k] \Big) \\
\le \|\nabla \psi_\delta(|\mathcal{A} y_{k-1}''|^2,|b_{k-1}'|^2)\|_{L^\infty(0,L)}
\big(\|y_{k-1}''\| \| d_t y_k''\| + \|b_{k-1}'\| \|d_t b_k''\| \big).
\end{multline}
Here, we used that the separate convexity of $P_{\varepsilon,h}$ implies that we have
\[
P_{\varepsilon,h}[y_k,b_{k-1}] + \delta_y P_{\varepsilon,h}[y_k,b_{k-1};y_{k-1}-y_k]
\le P_{\varepsilon,h}[y_{k-1},b_{k-1}]
\]
and
\[
P_{\varepsilon,h}[y_k,b_k] + \delta_b P_{\varepsilon,h}[y_k,b_k;b_{k-1}-b_k]
\le P_{\varepsilon,h}[y_k,b_{k-1}]
\]
which by summation and multiplication by $1/\tau$ yields that
\[\begin{split}
d_t P_{\varepsilon,h}[y_k,b_k] &= \frac{1}{\tau}\big(P_{\varepsilon,h}[y_k,b_k]-P_{\varepsilon,h}[y_{k-1},b_{k-1}]\big) \\
&\le \delta_y P_{\varepsilon,h} [y_k,b_{k-1};d_t y_k] + \delta_b P_{\varepsilon,h}[y_k,b_k;d_t b_k].
\end{split}\]
In an inductive argument we assume that $E_{\delta,\varepsilon}[y_{k-1},b_{k-1}] \le e_0$
which yields
\[
\frac12 \|y_{k-1}''\|^2 + \frac52 \|b_{k-1}'\|^2 + P_{\varepsilon,h}[y_{k-1},b_{k-1}] \le e_0.
\]
With~\eqref{eq:rough_est} and Lemma~\ref{la:props_psi} we find that
\[\begin{split}
\tau  \Big(\frac12 \|d_t y_k\|_\star^2 & + \frac12 \|d_t b_k\|_\dagger^2\Big)
+ \Big(\frac12 \|y_k''\|^2 + \frac52 \|b_k'\|^2 + P_{\varepsilon,h}[y_k,b_k] \Big) \\
& \le e_0 + \tau \frac{c_{\psi,1}^2}{2} \max\{c_\star^2,c_\dagger^2\} \big(\|y_{k-1}''\|^2 + \|b_{k-1}'\|^2\big) \\
& \le e_0 + \tau \frac{c_{\psi,1}^2}{2} \max\{c_\star^2,c_\dagger^2\} 2 e_0
= (1+\tau C_0) e_0 \le 2 e_0,
\end{split}\]
where we assumed that $\tau C_0 \le 1$.  The next step shows that $E_{\delta,\varepsilon,h}[y^k,b^k]\le e_0$
and hence implies that the right-hand side can be replaced by~$e_0$. \\
{\em Step 2.} We again choose $w=d_t y_k$ and $r= d_t b_k$ in Steps~(1) and~(2) of
Algorithm~\ref{alg:grad_desc} but use the Taylor formula of Lemma~\ref{la:taylor_psi}
to verify that
\begin{multline*}
d_t \Big(\frac12 \|y_k''\|^2 + \frac52 \|b_k'\|^2
+ \int_0^L \psi_\delta(|\mathcal{A} y_k''|^2,|b_k'|^2) \dv{x}
+ P_{\varepsilon,h}[y_k,b_k] \Big) \\
+ \|d_t y_k\|_\star^2 + \|d_t b_k\|_\dagger^2
\le \int_0^L \Gamma(\widehat{y}_k,\widehat{b}_k,y_{k-1},b_{k-1},y_k,b_k) \dv{x},
\end{multline*}
where, using that $\|\mathcal{A} r\|_{L^\infty(0,L)} \le \|r\|_{L^\infty(0,L)}$,
\begin{multline*}
\int_0^L \Gamma(\widehat{y}_k,\widehat{b}_k,y_{k-1},b_{k-1},y_k,b_k) \dv{x} \\
\le c_\Gamma \tau \delta^{-1} \big(\|y_k''\|_{L^\infty(0,L)}^2 + \|y_{k-1}''\|_{L^\infty(0,L)}^2 +
\|b_k'\|_{L^\infty(0,L)}^2 + \|b_{k-1}'\|_{L^\infty(0,L)}^2+1\big) \\
\times  \big(\|d_t y_k''\|^2 + \|d_t b_k'\|^2\big).
\end{multline*}
Using the available bounds and the inverse estimate~\eqref{eq:inv_est},
we deduce that
\[
\int_0^L \Gamma(\dots) \dv{x} \le c_\Gamma \tau \delta^{-1} \big(4 c_{\rm inv} h^{-1} 8 e_0+1\big)
\max\{c_\star^2,c_\dagger^2\} \big(\|d_t y_k\|_\star^2 + \|d_t b_k\|_\dagger^2\big).
\]
This implies that we have
\[
d_t E_{\delta,\varepsilon}[y_k,b_k] + \big(1- c_{{\rm ed}} \tau \delta^{-1} h^{-1}\big)
\big(\|d_t y_k\|_\star^2 + \|d_t b_k\|_\dagger^2\big) \le 0.
\]
A summation over $k=1,2,\dots,K$ and multiplication by~$\tau$
yield the energy estimate for some~$K$, which implies the induction hypothesis. \\
{\em Step 3.} To derive the bounds that control the violation of the constraints
we note that at the nodes we have
\[
|b_k|^2 = |b_{k-1}|^2 + \tau^2 |d_t b_k|^2
= \dots = |b_0|^2 + \tau^2 \sum_{r=1}^k |d_t b_k|^2.
\]
Incorporating the identity $|b_0|^2 =1$ and taking the $L^\infty$ norm implies the
constraint violation estimate.  The same argument applies to $|y_k'|^2$.
\end{proof}

\section{Numerical Experiments}\label{sec:num_ex}

\subsection{Realization}
With the quadrature introduced in Proposition~\ref{prop:quadrature},
Algorithm~\ref{alg:grad_desc}
can be realized exactly and only requires solving linear systems of
equations. These define certain linear saddle-point problems with
unique solutions. We always use the inner products related to the
norms
\[
\|w\|_\star^2 = \|w\|^2 + \|w''\|^2, \quad \|r\|_\dagger^2 = \|r\|^2 + \|r'\|^2,
\]
and the relations
\[
\tau = h/10, \quad \varepsilon_1 = h, \quad \varepsilon_2 = h^{1/2}, \quad \delta = h^{1/2}.
\]
In view of Proposition~\ref{prop:ener_decay} this choice for~$\tau$ may
be too optimistic in general. However, we always observed energy monotonicity
in our experiments which may be related to an additional regularity property.
We confirm our theoretical findings by showing that our numerical scheme
reliably and efficiently detects
stationary configurations of low energy via discrete evolutions from given
initial states. In what follows we visualize the ribbon by introducing
an artificial positive width. In addition to the discretized dimensionally reduced
eleastic energy $E_{\delta,\varepsilon,h}$ and the penalty term $P_{\varepsilon,h}$
we also investigate the behavior of the bending and twist energies given by
\[
E_{{\rm bend}}[y] = \frac12 \int_0^L |y''|^2 \dv{x}, \quad
E_{{\rm twist}}[b] = \frac12 \int_0^L |b'|^2 \dv{x}.
\]
We always consider clamped boundary conditions for the deformation~$y$ and
Dirichlet boundary conditions for the director~$b$ which are specified
via the initial states, i.e., we impose
\[
y(x) = y_0(x), \quad y'(x) = y_0'(x), \quad b(x) = b_0(x).
\]
at the endpoints $x\in \{0,L\}$. To satisfy the condition $y_0'\cdot b_0 = 0$
the vector field $b_0$ is obtained by rotating a non-tangential vector field~$n_0$
around the tangent vector $y_0'$. The condition $y_0'\cdot b_0' =0$ is not
taken into account in the construction of $b_0$ but will be approximated via
the involved penalization.
The discrete initial states are defined
via nodal interpolation of the continuous initial states, i.e.,
\[
y_h^0 = \mathcal{I}_h^{3,1} y_0, \quad b_h^0 = \mathcal{I}_h^{1,0} b_0.
\]
All computations are carried out on triangulations that are given by
uniform partitions of the intervals $(0,L)$ into $N$ subintervals of
length $h=L/N$. We occasionally refer to the quantities $|y_h''|$ and
$|b_h'|$ as curvature and torsion of a ribbon. We always ran the discrete
gradient flow with the fixed time horizon $T=10$.

\subsection{M\"obius ribbon}
We consider boundary conditions that lead to the formation of a M\"obius
strip of small, vanishing width. The following example defines a twisted
frame for the unit circle which will deform into a stationary configuration
for the reduced elastic energy.

\begin{example}[M\"obius ribbon]\label{ex:moebius}
Set $L=2 \pi$ and define for $x\in (0,L)$
\[
y_0(x) = \big( \cos(x),\sin(x), 0\big)
\]
and with $\alpha= 3/2$
\[
b_0 = \cos(\alpha \cdot) (y_0' \times n_0) \times y_0'
+ \sin(\alpha \cdot) y_0' \times n_0 ,
\]
with $n_0 = [y_0']_2^\perp$ being the rotation of the first two components
of $y_0'$ by $\pi/2$, i.e., $[a]_2^\perp = (-a_2,a_1,0)$ for $a\in \mathbb{R}^3$.
\end{example}

Figure~\ref{fig:band_snaps_exp_1} shows snapshots of the discrete evolution
defined by Algorithm~\ref{alg:grad_desc} in Example~\ref{ex:moebius} for a partitioning
of $(0,L)$ into $N=320$ intervals. We observe an immediate change in shape
of the ribbon towards a M\"obius strip. The initially constant curvature
becomes nonconstant and appears to vanish in neighborhoods of three points
in the nearly stationary configuration. This effect can also be seen in
Figure~\ref{fig:band_finals_exp_1} where the stationary state is visualized
from two other perspectives. In points where curvature vanishes
the original Sadowsky functional is singular and the correction from~\cite{FHMP16}
becomes particularly relevant. The pure bending energy defined by the deformation
of the centerline
increases during the discrete gradient flow evolution, while the twist energy
and the dimensionally reduced energy decay monotonically and quickly become
nearly stationary, cf.~Figure~\ref{fig:ener_exp_1}. The discrete director
fields $b_h^0$ and $b_h^k$ with $k=5092$ are displayed at selected deformed
nodes in Figure~\ref{fig:directors_exp_1}. An experimental convergence
behavior of the stationary energies
for discretizations with $h=L/N$, $N=80, 160, 320, 640$ is provided
in Table~\ref{tab:conv_eners_exp_1}. We observe that the energies increase
with finer resolutions and their differences become smaller. These results also
indicate that the scheme does not tend to get stuck at local minima.

\begin{table}
\begin{tabular}{|l|c|c|c|c|} \hline
$N\sim h^{-1}$ & 80 & 160 & 320 & 640 \\\hline
$E^{{\rm quad}}_{\delta,\varepsilon,h}[y_h,b_h]$ & 14.9255 & 32.0886 & 35.5842 & 37.6865 \\ \hline
\end{tabular}
\vspace*{2mm}
\caption{\label{tab:conv_eners_exp_1} Stationary energy values in
Example~\ref{ex:moebius} for different mesh-sizes $h=L/N$.}
\end{table}

\begin{figure}[p]
\includegraphics[width=6.5cm]{band_exp_1_0.jpg} \hspace*{-8mm}
\includegraphics[width=6.5cm]{band_exp_1_102.jpg} \\[-3mm]
\includegraphics[width=6.5cm]{band_exp_1_204.jpg}  \hspace*{-8mm}
\includegraphics[width=6.5cm]{band_exp_1_306.jpg} \\[-3mm]
\includegraphics[width=6.5cm]{band_exp_1_408.jpg} \hspace*{-8mm}
\includegraphics[width=6.5cm]{band_exp_1_510.jpg} \\[-3mm]
\includegraphics[width=6.5cm]{band_exp_1_2448.jpg} \hspace*{-8mm}
\includegraphics[width=6.5cm]{band_exp_1_5092.jpg}
\caption{\label{fig:band_snaps_exp_1} Snapshots of the discrete gradient flow
in Example~\ref{ex:moebius} after $k=0, 102, 204, 306, 408, 510, 2448, 5092$ iterations
leading to a M\"obius strip, colored by curvature and torsion;
crosses indicate where torsion dominates curvature.}
\end{figure}

\begin{figure}[p]
\includegraphics[width=5.6cm]{final_A_band_exp_1.jpg}  \hspace*{-3mm}
\includegraphics[width=5.6cm]{final_B_band_exp_1.jpg} \\
\caption{\label{fig:band_finals_exp_1} Nearly stationary configuration in
Example~\ref{ex:moebius} from different perspectives colored by curvature and
torsion; crosses indicate where torsion dominates curvature.}
\end{figure}

\begin{figure}
\includegraphics[width=5.8cm]{director_exp_1_0.jpg}
\includegraphics[width=5.8cm]{director_exp_1_5092.jpg} \\
\caption{\label{fig:directors_exp_1} Initial and nearly stationary director
fields on the deformed centerline
in Example~\ref{ex:moebius}; the vectors at every fourth node are displayed.}
\end{figure}

\begin{figure}
\includegraphics[width=7.7cm]{ener_exp_1.pdf}
\caption{\label{fig:ener_exp_1} Energy decay in Example~\ref{ex:moebius}
accompanied by an increase of bending energy; twist energy and penalty functional
decay monotonically.}
\end{figure}

\subsection{Twisted helix}
A helix is a canonical example of a twisted band. We consider here the
evolution from a helical configuration where the torsion is inconsistent
with the winding of the helix.

\begin{example}[Twisted helical ribbon]\label{ex:helix}
Let $L = 2\pi$ and define for $x\in (0,L)$
\[
y_0(x) = \big(c x, d  \cos(\beta x), d \sin(\beta x) \big),
\]
where $c=0.95$ and $d= (1-c^2)^{1/2}/\beta$ with $\beta=2$, and
with $\alpha=1$ let
\[
b_0 = \Pi_{S_1}\big[\cos(\alpha \cdot) (y_0' \times n_0) \times y_0'
+ \sin(\alpha \cdot ) y_0' \times n_0\big] ,
\]
where $n_0 = [y_0']_2^\perp$, cf.~Example~\ref{ex:moebius}, and
$\Pi_{S_1}:\mathbb{R}^3\setminus \{0\}\to \mathbb{R}^3$ denotes the projection
onto the unit sphere.
\end{example}

The application of Algorithm~\ref{alg:grad_desc} to the initial
data defined by Example~\ref{ex:helix} leads to the snapshots
displayed in Figure~\ref{fig:band_snaps_exp_2}. In contrast
to Example~\ref{ex:moebius} we observe in Example~\ref{ex:helix}
that in the nearly stationary configuration curvature is dominated
by torsion along the entire centerline, cf.~Figure~\ref{fig:band_finals_exp_2}.
Interestingly, this is also the case for the inital configuration but
changes during the evolution. The discrete, dimensionally reduced energy,
the twist energy, and the penalty value decrease montonically during
the evolution while the bending energy shows a moderate increase as
can be seen in the plot of Figure~\ref{fig:ener_exp_2}. The initial
and nearly stationary director fields are displayed in
Figure~\ref{fig:directors_exp_2}. The stationary energies for a
sequence of refined partitions are provided in Table~\ref{tab:conv_eners_exp_2}.
In comparison with Example~\ref{ex:moebius} we observe here a non-monotone
behavior and smaller differences.

\begin{table}
\begin{tabular}{|l|c|c|c|c|} \hline
$N\sim h^{-1}$ & 80 & 160 & 320 & 640 \\\hline
$E^{{\rm quad}}_{\delta,\varepsilon,h}[y_h,b_h]$ & 27.7554 & 26.5432 & 26.3554 & 26.4050\\ \hline
\end{tabular}
\vspace*{2mm}
\caption{\label{tab:conv_eners_exp_2} Stationary energy values in
Example~\ref{ex:helix} for different mesh-sizes $h=L/N$.}
\end{table}

\begin{figure}[p]
\includegraphics[width=6.5cm]{band_exp_2_0.jpg} \hspace*{-8mm}
\includegraphics[width=6.5cm]{band_exp_2_408.jpg} \\[-2mm]
\includegraphics[width=6.5cm]{band_exp_2_816.jpg}  \hspace*{-8mm}
\includegraphics[width=6.5cm]{band_exp_2_1224.jpg} \\[-2mm]
\includegraphics[width=6.5cm]{band_exp_2_1632.jpg} \hspace*{-8mm}
\includegraphics[width=6.5cm]{band_exp_2_2040.jpg} \\[-2mm]
\includegraphics[width=6.5cm]{band_exp_2_2448.jpg} \hspace*{-8mm}
\includegraphics[width=6.5cm]{band_exp_2_5092.jpg}
\caption{\label{fig:band_snaps_exp_2} Snapshots of the discrete gradient flow
in Example~\ref{ex:helix} after $k=0, 408, 816, 1224, 1632, 2040, 2448, 5092$ iterations
starting from a helix, colored by curvature and torsion;
crosses indicate where torsion dominates curvature.}
\end{figure}

\begin{figure}[p]
\includegraphics[width=5.6cm]{final_C_band_exp_2.jpg}  \hspace*{-3mm}
\includegraphics[width=5.6cm]{final_B_band_exp_2.jpg} \\
\caption{\label{fig:band_finals_exp_2} Nearly stationary configuration in
Example~\ref{ex:helix} from different perspectives colored by curvature and
torsion; crosses indicate where torsion dominates curvature.}
\end{figure}

\begin{figure}
\includegraphics[width=5.8cm]{director_exp_2_0.jpg}
\includegraphics[width=5.8cm]{director_exp_2_5092.jpg} \\
\caption{\label{fig:directors_exp_2} Initial and nearly stationary director
fields on the deformed centerline
in Example~\ref{ex:helix}; the vectors at every fourth node
are displayed.}
\end{figure}

\begin{figure}
\includegraphics[width=7.7cm]{ener_exp_2.pdf}
\caption{\label{fig:ener_exp_2} Energy decay in Example~\ref{ex:helix}
accompanied by a small increase of bending energy; twist energy and penalty functional
decay monotonically.}
\end{figure}

\medskip
\noindent
{\em Acknowledgments.} The author acknowledges stimulating discussions
with Peter Hornung. The author would like to thank the Isaac Newton
Institute for Mathematical Sciences for support and hospitality
during the programme {\em Geometry, compatibility and structure preservation
in computational differential equations} when work on this paper
was undertaken. This work was supported by EPSRC grant number EP/R014604/1.
The author also acknowledges support by the German Research Foundation (DFG)
via the research unit FOR 3013 {\em Vector- and tensor-valued surface PDEs}.

\bibliographystyle{abbrv}


\section{Introduction}
The success of information retrieval in a given task critically depends on the quality of the underlying data. Another issue is that in many domains knowledge bases are spread across various data sources \cite{nickel2015review} and it is crucial to be able to combine information from different sources.
In this work, we focus on knowledge bases in the form of Knowledge Graphs (KGs), which are particularly suited for information retrieval \cite{singhal2012introducing}.
Joining information from different KGs is non-trivial, as there is no unified schema or vocabulary. The goal of the entity alignment task is to overcome this problem by \emph{learning} a matching between entities in different KGs. In the typical setting some of the alignments are known in advance (seed alignments) and the task is therefore supervised. More formally, we are given graphs $G_L = (V_L, E_L)$ and $G_R = (V_R, E_R)$
with a seed alignment $A = {(l_i, r_i)}_i \subseteq V_L \times V_R$.
It is commonly assumed that an entity $v \in V_L$ can match at most one entity $v' \in V_R$.
Thus the goal is to infer alignments for the remaining nodes only.

Graph Convolutional Networks(GCN) \cite{kipf2016semi,gilmer2017neural}, which have been recently become increasingly popular, are at the core of state-of-the-art methods for entity alignments in KGs \cite{DBLP:conf/emnlp/WangLLZ18,DBLP:conf/acl/CaoLLLLC19,xu2019cross,DBLP:conf/ijcai/ZhuZ0TG19,anonymous2020deep}.
In this paper, we thoroughly analyze one of the first GCN-based entity alignment methods, GCN-Align~\cite{DBLP:conf/emnlp/WangLLZ18}.
Since the other methods we are studying can be considered as extensions of this first paper and have a similar architecture, our goal is to understand the importance of its individual components and architecture choices.
\noindent In summary, our contribution is as follows:
\begin{enumerate}
    \item We investigate the reproducibility of the published results of a recent GCN-based method for entity alignment and uncover differences between the method's description in the paper and the authors' implementation.
    \item We perform an ablation study to demonstrate the individual components' contribution.
    \item We apply the method to numerous additional datasets of different sizes to investigate the consistency of results across datasets.
\end{enumerate}

\section{Related work}
In this section we review previous work for the entity alignment for Knowledge Graphs and revise datasets and current evaluation process. We believe this is useful for practitioners, since we discover some pitfalls, especially when implementing evaluation scores and selecting datasets for comparison. The overview of methods, datasets and metrics is provided in Table~\ref{tab:related}.
\begin{table}
    \centering
    \caption{Overview of related work in the field of entity alignment for knowledge graphs with their used datasets and metrics.}
    \label{tab:related}
    \scriptsize
\begin{tabular*}{\linewidth}{*{4}{@{\extracolsep{\fill}}l}@{\extracolsep{\fill}}}
\toprule
Method & Datasets & Metrics & Code\\ %
\midrule
MTransE \cite{DBLP:conf/ijcai/ChenTYZ17} & WK3l-15K, WK3l-120K, CN3l & H@10(, MR) & yes\\%
IPTransE \cite{DBLP:conf/ijcai/ZhuXLS17} & DFB-\{1,2,3\} & H@\{1,10\}, MR & yes\\%
JAPE \cite{DBLP:conf/semweb/SunHL17} & DBP15K(JAPE) & H@\{1,10,50\}, MR & yes\\%
KDCoE \cite{DBLP:conf/ijcai/ChenTCSZ18} & WK3l-60K & H@\{1,10\}, MR & yes\\%
BootEA \cite{DBLP:conf/ijcai/SunHZQ18} & DBP15K(JAPE), DWY100K & H@\{1,10\}, MRR & yes\\%
SEA \cite{DBLP:conf/www/PeiYHZ19} & WK3l-15K, WK3l-120K & H@\{1,5,10\}, MRR & yes\\%
MultiKE \cite{DBLP:conf/ijcai/ZhangSHCGQ19} & DWY100K & H@\{1,10\}, MR, MRR & yes\\%
AttrE \cite{DBLP:conf/aaai/TrisedyaQZ19} & DBP-LGD,DBP-GEO,DBP-YAGO & H@\{1,10\}, MR & yes\\%
RSN \cite{DBLP:conf/icml/GuoSH19} & custom DBP15K, DWY100K & H@\{1,10\}, MRR &  yes\\%
\textbf{GCN-Align} \cite{DBLP:conf/emnlp/WangLLZ18} & DBP15K(JAPE) & H@\{1,10,50\}& yes\\%
CL-GNN \cite{DBLP:conf/acl/XuWYFSWY19} & DBP15K(JAPE) & H@\{1,10\} & yes\\%
MuGNN \cite{DBLP:conf/acl/CaoLLLLC19} & DBP15K(JAPE), DWY100K & H@\{1,10\}, MRR & yes\\%
NAEA \cite{DBLP:conf/ijcai/ZhuZ0TG19} & DBP15K(JAPE), DWY100K & H@\{1,10\}, MRR & no \\
\bottomrule
\end{tabular*}
 \end{table}
\subsection{Methods}
While the problem of entity alignments in Knowledge Graphs has been tackled historically by researching vocabularies which are as broad as possible, and establish them as a standard, recent approaches take a more data-driven view.
Early methods use classical knowledge graph link prediction models such as TransE~\cite{DBLP:conf/nips/BordesUGWY13} to embed the entities of the individual knowledge graphs using a intra-KG link prediction loss, and differ in what they do with the aligned entities.
For instance MTransE~\cite{DBLP:conf/ijcai/ChenTYZ17} learns a linear transformation between the embedding spaces of the individual graphs using $L_2$-loss.
BootEA~\cite{sun2018bootstrapping} adopts a bootstrapping approach and iteratively labels the most likely alignments and utilizes them for further training. In addition to the alignment loss, embeddings of aligned entities are swapped regularly to calibrate embedding spaces against each other. SEA~\cite{pei2019semi} learns mapping between embedding spaces in both directions and additionaly adds cycle-consistency loss. Therefore the distance between original embedding of an entity and its representation, which was first translated to another space and then back from it, is penalized.
IPTransE~\cite{DBLP:conf/ijcai/ZhuXLS17} embeds both KGs into the same embedding space and uses a margin-based loss to enforce the embeddings of aligned entities to become similar.
RSN~\cite{DBLP:conf/icml/GuoSH19} model generates sequences using different types of random walks which can move between graphs when visiting aligned entities. The generated sequences are feed to adapted recurrent model.
JAPE~\cite{DBLP:conf/semweb/SunHL17}, KDCoE~\cite{DBLP:conf/ijcai/ChenTCSZ18}, MultiKE~\cite{DBLP:conf/ijcai/ZhangSHCGQ19} and AttrE \cite{DBLP:conf/aaai/TrisedyaQZ19} utilize attributes available for some entities and additional information like names of entities and relationships.
Graph Neural Network (GNN) based models
\cite{DBLP:conf/emnlp/WangLLZ18,DBLP:conf/acl/CaoLLLLC19,xu2019cross,DBLP:conf/ijcai/ZhuZ0TG19,anonymous2020deep}\footnote{Please note, that while \cite{DBLP:conf/ijcai/ZhuZ0TG19} does not state explicitly that they use GNNs, their model is very similar to \cite{velivckovic2017graph}.} have in common that they use GNN to create node representations by aggregating node representations together with representations of their neighbors.
Most of GNN approaches do not distinguish between different relations and either consider all neighbors equally \cite{DBLP:conf/emnlp/WangLLZ18,xu2019cross,anonymous2020deep} or use attention \cite{DBLP:conf/acl/CaoLLLLC19} to weight the representations of the neighbors for the aggregation.

\subsection{Datasets}
\begin{table}[ht]
    \centering
    \caption{
    Overview of used datasets with their sizes in the number of triples (edges), entities (nodes), relations (different edge types) and alignments.
    For Wk3l, the alignment is provided as a directed mapping on a entity level.
    However, there are additional triple alignments.
    Following a common practice as e.g. \cite{DBLP:conf/www/PeiYHZ19} we can assume that an alignment should be symmetric, and that we can extract entity alignments from the triple alignments.
    Doing so, we obtain the number of alignments given in brackets.
    }
    \label{tab:my_label}
    \scriptsize
    \begin{tabular*}{\linewidth}{*{2}{l}*{5}{@{\extracolsep{\fill}}r}@{\extracolsep{\fill}}}
    \toprule
    Dataset & Subset & Graph & Triples & Entities & Relations & Alignments \\
    \midrule
    \multirow{6}{*}{DBP15k (full)}
    & \multirow{2}{*}{fr-en} & fr &
    192,191 & 66,858 & 1,379 & \multirow{2}{*}{15,000}\\
    & & en &
    278,590 & 105,889 & 2,209 & \\
    & \multirow{2}{*}{ja-en} & ja &
    164,373 & 65,744 & 2,043 & \multirow{2}{*}{15,000}\\
    & & en &
    233,319 & 95,680 & 2,096 & \\
    & \multirow{2}{*}{zh-en} & zh &
    153,929 & 66,469 & 2,830 & \multirow{2}{*}{15,000}\\
    & & en &
    237,674 & 98,125 & 2,317 & \\
    \midrule
    \multirow{6}{*}{DBP15k (JAPE)}
    & \multirow{2}{*}{fr-en} & fr &
    105,998 & 19,661 & 903 & \multirow{2}{*}{15,000}\\
    & & en &
    115,722 & 19,993 & 1,208 & \\
    & \multirow{2}{*}{ja-en} & ja &
    77,214 & 19,814 & 1,299 & \multirow{2}{*}{15,000}\\
    & & en &
    93,484 & 19,780 & 1,153 & \\
    & \multirow{2}{*}{zh-en} & zh &
    70,414 & 19,388 & 1,701 & \multirow{2}{*}{15,000}\\
    & & en &
    95,142 & 19,572 & 1,323 & \\
    \midrule
    \multirow{4}{*}{WK3l-15k}
    & \multirow{2}{*}{en-de} & en &
    209,041 & 15,127 & 1,841 & 1,289 (10,383) \\
    & & de &
    144,244 & 14,603 & 596 & 1,140 (10,383) \\
    & \multirow{2}{*}{en-fr} & en &
    203,356 & 15,170 & 2,228 & 2,498 (\phantom{1}8,024) \\
    & & fr &
    169,329 & 15,393 & 2,422 & 3,812 (\phantom{1}8,024) \\
    \midrule
    \multirow{4}{*}{WK3l-120k}
    & \multirow{2}{*}{en-de} & en &
    624,659 & 67,650 & 2,393 & 6,173 (50,280) \\
    & & de &
    389,554 & 61,942 & 861 & 4,820 (50,280) \\
    & \multirow{2}{*}{en-fr} & en &
    1,375,406 & 119,749 & 3,109 & 36,749 (87,836) \\
    & & fr &
    760,497 & 118,592 & 2,336 & 36,013 (87,836) \\
    \midrule
    \multirow{4}{*}{DWY-100k}
    & \multirow{2}{*}{dbp-wd} & dbp &
    463,294 & 100,000 & 330 & \multirow{2}{*}{100,000}\\
    & & wd &
    448,774 & 100,000 & 220 & \\
    & \multirow{2}{*}{dbp-yg} & dbp &
    428,952 & 100,000 & 302 & \multirow{2}{*}{100,000}\\
    & & yg &
    502,563 & 100,000 & 31 & \\
    \bottomrule
    \end{tabular*}
\end{table}
The datasets used by entity alignments methods generally derive from large-scale open-source data source such as DBPedia~\cite{DBLP:conf/semweb/AuerBKLCI07}, YAGO~\cite{DBLP:conf/cidr/MahdisoltaniBS15}, or Wikidata~\cite{wikidata}.
While there is the DWY-100k dataset, which comprises 100k aligned entities across the three aforementioned individual knowledge graphs, most of the datasets, such as DBP15k, or WK3l derive from a single multi-lingual database.
There, subsets are formed according to a specific language, and entities which occur in multiple languages and are linked accordingly are used as alignments.

As an interesting observation we found out that all papers which evaluate on DBP15k, do not evaluate on the full DBP15k dataset\footnote{Available at \url{http://ws.nju.edu.cn/jape/}} (which we refer to as \emph{DBP15k (full)}), but rather use a smaller subset provided by the authors of JAPE~\cite{DBLP:conf/semweb/SunHL17} in their GitHub repository\footnote{\url{https://github.com/nju-websoft/JAPE/blob/master/data/dbp15k.tar.gz}}, which we call \emph{DBP15k-JAPE}.
The smaller subsets were created by selecting a portion of entities (around 20k of 100k) which are popular, i.e. appear in many triples as head or tail.
The number of aligned entities stays the same (15k).
As the paper only reports the dataset statistics of the larger dataset, and does not mention the reduction of the dataset, subsequent papers
also report the statistics of the larger dataset, although experiments use the smaller variant \cite{DBLP:conf/semweb/SunHL17,sun2018bootstrapping,DBLP:conf/emnlp/WangLLZ18,DBLP:conf/acl/CaoLLLLC19,DBLP:conf/ijcai/ZhuXLS17}.

\subsection{Scores}
It is common practice to only consider the entities being part of the test alignment as potential matching candidates.
Although we argue that ignoring entities exclusive to a single graph as potential candidates does not reflect well the use-case situation\footnote{In the typical scenario it is not known in advance, which entities have matching and which not. Therefore the resulting score is too optimistic. However, we advocate to investigate this shortcoming further in future work}, we follow this evaluation scheme for our experiments to maintain comparability.

In the following description of evaluation measures we focus only on the case of aligning one node $l_i \in V_L$ with a ground truth alignment $r_i \in V_R$.
The right-to-left alignment is handled analogously.
Let $V_R^* = \{v_r \in V_R \mid \exists v_l \in V_L: (v_l, v_r) \in A_e\}$ denote the set of matching candidates in the right graph.
For a node $l_i$, the entity alignment models generates a score $f(l_i, v_j)$ for each matching candidate $v_j \in V_R^*$.
Afterwards, the candidates are sorted according to their score, and the rank $rank(l_i, r_i)$ is computed as the index of the ground truth match $r_i$ in this sorted list (1-based).
The \emph{mean rank (MR)} is simply the mean over the ranks for all alignments.
$$
MR = \frac{1}{|A_e|}\sum \limits_{i=1}^{|A_e|} rank(l_i, r_i)
$$
The \emph{mean reciprocal rank (MRR)} is the mean over all reciprocal ranks.
$$
MRR = \frac{1}{|A_e|}\sum \limits_{i=1}^{|A_e|} \frac{1}{rank(l_i, r_i)}
$$
It is naturally bounded between 0 and 1, where 1 corresponds to a perfect score.
Moreover, its value is dominated by small ranks, and it is less sensitive to larger ones.
The \emph{hits at $k$ (H@k)} is the percentage of alignments where the rank was at most $k$, i.e. equivalent to the recall at $k$.
$$
H@k = \frac{|\{(l_i, r_i) \in A_e \mid rank(l_i, r_i) \leq k\}|}{|A_e|}
$$ %
\section{Method}
GCN-Align~\cite{DBLP:conf/emnlp/WangLLZ18} is a GCN-based approach to embed all entities from both graphs into a common embedding space.
Each entity $i$ is associated with \emph{structural} features $h_i^s \in \mathbb{R}^d$, which are initialized randomly and updated during training.
The features of all entities in a single graph are combined to the feature matrix $H^s$.
Subsequently, a two-layer GCN is applied.
A single GCN layer is described by
$$
H^{(i+1)} = \sigma\left(\hat{D}^{-\frac{1}{2}}\hat{A}\hat{D}^{-\frac{1}{2}}H^{(i)}W^{(i)}\right)
$$
with $\hat{A} = A + I$, where $A$ is the adjacency matrix, and $\hat{D}_{ii} = \sum \limits_{j=1}^n \hat{A}_{ij}$ is the diagonal node degree matrix.
The input of the first layer is set to $H^{(0)} = H^s$, and $\sigma$ is non-linear activation function.
For the first layer, $\sigma=ReLU$, and the second layer uses the identity.
The output of the second layer is considered as the structural embedding, denoted by $s_i = H^{(2)}_i \in \mathbb{R}^{d}$.
Both graphs are equipped with its own node features, but the convolution weights $W^{(i)}$ are shared across the graphs.

The adjacency matrix is derived from the knowledge graph by first computing a score, called \emph{functionality}, for each relation as the ratio between the number of different entities which occur as head, and the number of triples in which the relation occurs
$$
\alpha_r := \frac{|\{v \in V \mid \exists v': (v, r, v') \in T\}|}{|\{t \in T \mid \exists v, v' \in V: (v, r, v) \in T\}|}
$$
Analogously, the \emph{inverse functionality} $\alpha_r'$ is obtained by replacing the nominator by the number of different tail entities.
The final adjacency matrix is obtained as
$$
A_ij = \sum \limits_{(e_i, r, e_j)} \alpha_r' + \sum \limits_{(e_j, r, e_j)} \alpha_r
$$

In addition, each entity $i$ is also equipped with attributes $h_i^a \in \mathbb{R}^{d'}$ which are combined into a graph attribute matrix $H^a$.
Anagously, the attributes are processed by a GCN for each graph with convolution weights shared across the graphs, resulting in attribute embeddings $a_i \in \mathbb{R}^{d'}$.

The attribute and structure GCNs are optimized separately using SGD.
As loss function, a margin-rank loss is used, exemplary for the structure embedding
$$
L = \sum \limits_{(r_i, l_i) \in A} \sum \limits_{(r_j, l_j) \in A_i^-} \left[\|s_i^L - s_i^R\|_1 + \gamma - \|s_j^L - s_j^R\|\right]_{+}
$$
where $[x]_+ = \max \{0, x\}$, and the margin $\gamma$ is a hyperparameter chosen separately for structure and attribute embeddings.
$A_i^-$ denotes a set of negative samples constructed by either replacing the left or the right entity with a random entity from the same graph.

In order to compare two nodes from both graphs, the $L_1$ distance between their embeddings is used, normalised by the dimensionality.
$$
score(v_i^L, v_j^R) = - \left(\beta \frac{\|s_i^L - s_j^R\|_1}{d} + (1 - \beta) \frac{\|a_i^L - a_j^R\|_1}{d'} \right)
$$
Here, $\beta$ is a hyperparameter for the tradeoff between structural and attribute similarity.

\subsection{Implementation Differences}
The code\footnote{\url{https://github.com/1049451037/GCN-Align}} provided by the authors differs in a few aspects from the method described in the paper.
First, instead of using the full DBP15k dataset having the dataset sizes as reported in the paper, a smaller version is used.
Second, when computing the adjacency matrix, $fun(r)$ and $ifun(f)$ are set to at least 0.3.
Third, the node features are always normalised to unit Euclidean length before passing them into the network.
Finally, there are no convolution weights.
This fact is particularly interesting, as this means that the whole GCN does not contain a single parameter, but is just a fixed function on the learned node embeddings. %

\section{Experiments}
In initial experiments we were able to reproduce the results reported in the paper using the implementation provided by the authors.
Moreover, we are able to reproduce the results using our own implementation, and settings adjusted to the authors' code.
In addition, we replaced the adjacency matrix based on functionality and inverse functionality by a simpler version, where $a_{ij} = \{(h,r,t) \in T \mid h=e_i, t=e_j\}$.
We additionally use $\hat{D}^{-1}\hat{A}$ instead of the symmetric normalization.
In total, we see no difference in performance between our simplified adjacency matrix, and the authors' one.
We identified two aspects which affect the model's performance:
Not using convolutional weights, and normalizing the variance when initializing node embeddings.
We provide empirical evidence for this finding across numerous datasets.
\begin{table}
    \centering
    \caption{Hyperparameter grid used for large-scale hyperparameter search on \emph{DBP15K (JAPE) zh-en}.}
    \begin{tabular*}{\linewidth}{l@{\extracolsep{\fill}}l@{\extracolsep{\fill}}r}
        \toprule
        hyperparameter & abbrev. & value range\\
        \midrule
        optimizer & opt & \{Adam, SGD\}\\
        learning rate & lr & \{0.1, 0.5, 1, 10, 20\}\\
        number of layers & \#layers & \{1, 2, 3\}\\
        number of negative samples & \#neg & \{5, 50, 100\}\\
        number of epochs & \#epochs & \{10, 500, 2000, 3000\}\\
        \bottomrule
    \end{tabular*}
    \label{tab:hyperparameters}
\end{table}

\begin{table}
    \centering
    \caption{Optimal Hyperparameters found for DBP15k (JAPE), zh-en with and without convolution weights, and with two different embedding initialization variances.}
    \begin{tabular*}{\linewidth}{l*{4}{@{\extracolsep{\fill}}r}}
        \toprule
        \textbf{Weights} & \multicolumn{2}{c}{no} & \multicolumn{2}{c}{yes} \\
        \textbf{Variance Emb. Init.} &  1 & $n^{-1/2}$ &  1 & $n^{-1/2}$ \\
        \midrule
        \#epochs    &  2,000 &  3,000 &  2,000 &  2,000 \\
        \#neg &    50 &   100 &    50 &    50 \\
        \#layers      &     2 &     2 &     3 &     2 \\
        lr            &     1 &     1 &     1 &     1 \\
        opt     &  adam &   sgd &  adam &  adam \\
        \bottomrule
    \end{tabular*}
    \label{tab:best_zh_en}
\end{table}

\begin{table}
    \centering
    \caption{
    Ablation study on using convolution weights and different embedding initialisation.
    A detailed description can be found in the text.
    }
    \label{tab:result}
    \scriptsize
\begin{tabular*}{\linewidth}{*{2}{l}*{10}{@{\extracolsep{\fill}}r}@{\extracolsep{\fill}}}
\toprule
        & \textbf{Weights} & \multicolumn{2}{c}{No} & \multicolumn{2}{c}{Yes}\\
         \multicolumn{2}{r}{\textbf{Variance. Emb. Init.}} & $1$ & $n^{-1/2}$ & $1$ &  $n^{-1/2}$ \\
         \midrule
         & {} & \multicolumn{4}{c}{H@1} \\
        \midrule
        \textbf{DBP15K (full)} & \textbf{fr-en} &  31.51 $\pm$ 0.16 &  27.64 $\pm$ 0.22 &  21.82 $\pm$ \phantom{0}0.39 &  16.73 $\pm$ \phantom{0}0.59\\
        & \textbf{ja-en} &  33.26 $\pm$ 0.10 &  29.06 $\pm$ 0.23 &  26.21 $\pm$ \phantom{0}0.33 &  20.78 $\pm$ \phantom{0}0.16 \\
        & \textbf{zh-en} &  31.15 $\pm$ 0.15 &  22.55 $\pm$ 0.27 &  24.96 $\pm$ \phantom{0}0.71 &  18.85 $\pm$ \phantom{0}0.99 \\
        \midrule
\textbf{DBP15K (JAPE)} & \textbf{fr-en} &  45.37 $\pm$ 0.13 &  41.03 $\pm$ 0.13 &  35.36 $\pm$ \phantom{0}0.33 &  30.50 $\pm$ \phantom{0}0.38 \\
        & \textbf{ja-en} &  45.53 $\pm$ 0.18 &  40.29 $\pm$ 0.09 &  35.81 $\pm$ \phantom{0}0.53 &  31.46 $\pm$ \phantom{0}0.15 \\
        & \textbf{zh-en} &  43.30 $\pm$ 0.12 &  39.37 $\pm$ 0.20 &  33.61 $\pm$ \phantom{0}0.49 &  29.94 $\pm$ \phantom{0}0.35 \\
        \midrule
\textbf{DWY100K} & \textbf{wd} &  58.50 $\pm$ 0.05 &  54.07 $\pm$ 0.05 &  50.13 $\pm$ \phantom{0}0.11 &  38.85 $\pm$ \phantom{0}0.31 \\
        & \textbf{yg} &  72.82 $\pm$ 0.06 &  67.06 $\pm$ 0.03 &  67.36 $\pm$ \phantom{0}0.10 &  60.67 $\pm$ \phantom{0}0.30 \\
        \midrule
\textbf{WK3l-120K} & \textbf{en-de} &  10.10 $\pm$ 0.03 &   9.17 $\pm$ 0.05 &   9.02 $\pm$ \phantom{0}0.17 &   6.75 $\pm$ \phantom{0}0.12 \\
        & \textbf{en-fr} &   8.28 $\pm$ 0.03 &   7.38 $\pm$ 0.03 &   7.26 $\pm$ \phantom{0}0.11 &   5.07 $\pm$ \phantom{0}0.16 \\
        \midrule
\textbf{WK3l-15K} & \textbf{en-de} &  16.57 $\pm$ 0.12 &  14.41 $\pm$ 0.23 &  17.43 $\pm$ \phantom{0}0.38 &  12.66 $\pm$ \phantom{0}0.30 \\
        & \textbf{en-fr} &  17.07 $\pm$ 0.15 &  16.16 $\pm$ 0.16 &  15.98 $\pm$ \phantom{0}0.16 &  12.41 $\pm$ \phantom{0}0.18 \\
 \midrule
     & {} & \multicolumn{4}{c}{H@10} \\
\midrule
\textbf{DBP15K (full)} & \textbf{fr-en} &  68.38 $\pm$ 0.32 &  63.41 $\pm$ \phantom{0}0.14 &  59.26 $\pm$ \phantom{0}0.55 &  48.55 $\pm$ \phantom{0}0.92 \\
        & \textbf{ja-en} &  68.22 $\pm$ 0.09 &  61.95 $\pm$ 0.17 &  61.12 $\pm$ \phantom{0}0.51 &  49.56 $\pm$ \phantom{0}0.38 \\
        & \textbf{zh-en} &  67.46 $\pm$ 0.11 &  56.03 $\pm$ 0.21 &  59.07 $\pm$ \phantom{0}1.10 &  50.32 $\pm$ \phantom{0}1.52 \\
        \midrule
\textbf{DBP15K (JAPE)} & \textbf{fr-en} &  82.48 $\pm$ 0.08 &  79.11 $\pm$ 0.07 &  74.71 $\pm$ \phantom{0}0.27 &  69.72 $\pm$ \phantom{0}0.36 \\
        & \textbf{ja-en} &  79.77 $\pm$ 0.14 &  75.13 $\pm$ 0.20 &  73.05 $\pm$ \phantom{0}0.52 &  67.18 $\pm$ \phantom{0}0.28 \\
        & \textbf{zh-en} &  77.63 $\pm$ 0.05 &  73.66 $\pm$ 0.28 &  71.16 $\pm$ \phantom{0}0.17 &  66.22 $\pm$ \phantom{0}0.51 \\
        \midrule
\textbf{DWY100K} & \textbf{wd} &  86.26 $\pm$ 0.05 &  81.30 $\pm$ 0.03 &  79.65 $\pm$ \phantom{0}0.20 &  69.73 $\pm$ \phantom{0}0.25 \\
        & \textbf{yg} &  92.13 $\pm$ 0.04 &  87.57 $\pm$ 0.04 &  88.64 $\pm$ \phantom{0}0.09 &  83.76 $\pm$ \phantom{0}0.27 \\
        \midrule
\textbf{WK3l-120K} & \textbf{en-de} &  27.13 $\pm$ 0.02 &  24.92 $\pm$ 0.03 &  25.49 $\pm$ \phantom{0}0.26 &  20.83 $\pm$ \phantom{0}0.29 \\
        & \textbf{en-fr} &  23.73 $\pm$ 0.04 &  21.57 $\pm$ 0.05 &  22.16 $\pm$ \phantom{0}0.24 &  16.31 $\pm$ \phantom{0}0.35 \\
        \midrule
\textbf{WK3l-15K} & \textbf{en-de} &  42.43 $\pm$ 0.13 &  38.63 $\pm$ 0.17 &  45.24 $\pm$ \phantom{0}0.47 &  37.03 $\pm$ \phantom{0}0.30 \\
        & \textbf{en-fr} &  49.68 $\pm$ 0.15 &  48.18 $\pm$ 0.15 &  47.64 $\pm$ \phantom{0}0.53 &  41.87 $\pm$ \phantom{0}0.36 \\
\midrule
        & {} & \multicolumn{4}{c}{MR}\\
        \midrule
        \textbf{DBP15K (FULL)} & \textbf{fr-en} &   203.90 $\pm$ 3.80 &   262.24 $\pm$ 3.23 &   123.09 $\pm$ 15.43 &   208.00 $\pm$ 12.04 \\
        & \textbf{ja-en} &   206.17 $\pm$ 4.21 &   358.53 $\pm$ 3.65 &   138.80 $\pm$ 12.87 &   238.24 $\pm$ 24.09 \\
        & \textbf{zh-en} &   168.80 $\pm$ 2.59 &   149.08 $\pm$ 2.70 &   279.49 $\pm$ 38.78 &   206.36 $\pm$ 17.60 \\
        \midrule
\textbf{DBP15K (JAPE)} & \textbf{fr-en} &   109.64 $\pm$ 1.56 &   117.59 $\pm$ 2.91 &    130.75 $\pm$ \phantom{0}8.48 &    133.14 $\pm$ \phantom{0}7.09 \\
        & \textbf{ja-en} &   144.81 $\pm$ 1.89 &   195.19 $\pm$ 3.44 &    146.42 $\pm$ \phantom{0}6.45 &   221.92 $\pm$ 12.22 \\
        & \textbf{zh-en} &   181.37 $\pm$ 4.05 &   215.23 $\pm$ 4.53 &   172.05 $\pm$ 12.72 &    236.72 $\pm$ \phantom{0}2.84 \\
         \midrule
\textbf{DWY100K} & \textbf{wd} &   277.08 $\pm$ 8.28 &   460.32 $\pm$ 9.17 &   500.61 $\pm$ 24.10 &   563.29 $\pm$ 28.92 \\

        & \textbf{yg} &    49.32 $\pm$ 2.71 &   102.50 $\pm$ 3.69 &    105.52 $\pm$ \phantom{0}4.63 &     67.71 $\pm$ \phantom{0}3.82 \\
         \midrule
\textbf{WK3l-120K} & \textbf{en-de} &  2753.75 $\pm$ 6.69 &  2280.31 $\pm$ 8.97 &  2843.96 $\pm$ 53.29 &  2289.02 $\pm$ 36.71 \\
        & \textbf{en-fr} &  4438.81 $\pm$ 9.29 &  4110.23 $\pm$ 7.90 &  4551.39 $\pm$ 55.29 &  4007.91 $\pm$ 59.01 \\
         \midrule
\textbf{WK3l-15K} & \textbf{en-de} &   247.74 $\pm$ 1.09 &   233.29 $\pm$ 2.66 &    263.16 $\pm$ \phantom{0}6.75 &    197.40 $\pm$ \phantom{0}5.39 \\
        & \textbf{en-fr} &   196.16 $\pm$ 1.09 &   176.32 $\pm$ 1.03 &    249.77 $\pm$ \phantom{0}7.71 &    184.72 $\pm$ \phantom{0}3.32 \\
\midrule
        & {} & \multicolumn{4}{c}{MRR}\\
        \midrule
        \textbf{DBP15K (full)} & \textbf{fr-en} &  43.59 $\pm$ 0.08 &  39.30 $\pm$ 0.18 &  33.83 $\pm$ \phantom{0}0.46 &  27.03 $\pm$ \phantom{0}0.72 \\
        & \textbf{ja-en} &  44.68 $\pm$ 0.06 &  39.92 $\pm$ 0.20 &  37.59 $\pm$ \phantom{0}0.37 &  30.31 $\pm$ \phantom{0}0.19 \\
        & \textbf{zh-en} &  43.09 $\pm$ 0.10 &  33.55 $\pm$ 0.19 &  36.15 $\pm$ \phantom{0}0.80 &  29.21 $\pm$ \phantom{0}1.16 \\
        \midrule
\textbf{DBP15K (JAPE)} & \textbf{fr-en} &  57.95 $\pm$ 0.10 &  53.78 $\pm$ 0.05 &  48.31 $\pm$ \phantom{0}0.26 &  43.37 $\pm$ \phantom{0}0.33 \\
        & \textbf{ja-en} &  57.14 $\pm$ 0.13 &  51.96 $\pm$ 0.07 &  48.03 $\pm$ \phantom{0}0.55 &  43.30 $\pm$ \phantom{0}0.15 \\
        & \textbf{zh-en} &  54.89 $\pm$ 0.09 &  50.88 $\pm$ 0.15 &  45.97 $\pm$ \phantom{0}0.39 &  41.93 $\pm$ \phantom{0}0.40 \\
        \midrule
\textbf{DWY100K} & \textbf{wd} &  68.33 $\pm$ 0.03 &  63.68 $\pm$ 0.04 &  60.50 $\pm$ \phantom{0}0.14 &  49.56 $\pm$ \phantom{0}0.30 \\
        & \textbf{yg} &  79.74 $\pm$ 0.04 &  74.29 $\pm$ 0.03 &  74.93 $\pm$ \phantom{0}0.09 &  68.76 $\pm$ \phantom{0}0.28 \\
        \midrule
\textbf{WK3l-120K} & \textbf{en-de} &  16.05 $\pm$ 0.03 &  14.73 $\pm$ 0.03 &  14.73 $\pm$ \phantom{0}0.21 &  11.70 $\pm$ \phantom{0}0.18 \\
        & \textbf{en-fr} &  13.65 $\pm$ 0.02 &  12.34 $\pm$ 0.02 &  12.41 $\pm$ \phantom{0}0.16 &   9.04 $\pm$ \phantom{0}0.23 \\
        \midrule
\textbf{WK3l-15K} & \textbf{en-de} &  25.40 $\pm$ 0.10 &  22.81 $\pm$ 0.16 &  26.94 $\pm$ \phantom{0}0.38 &  21.06 $\pm$ \phantom{0}0.25 \\
        & \textbf{en-fr} &  27.98 $\pm$ 0.16 &  26.76 $\pm$ 0.13 &  26.55 $\pm$ \phantom{0}0.23 &  22.20 $\pm$ \phantom{0}0.21 \\
\bottomrule
\end{tabular*}
\end{table}
We fix using convolution weights and the variance for the normal distribution from which the embedding vectors are initialized and optimize the other hyperparameters according to validation H@1 (80/20\% train-validation split) on \emph{DBP15K (JAPE) zh-en} in a large-scale hyperparameter search, comprising 1,440 experiments.
The hyperparameter grid is given in Table~\ref{tab:hyperparameters}, and Table~\ref{tab:best_zh_en} shows the best parameters found for DBP15k (JAPE) zh-en for the four different settings.
For each dataset, we perform a smaller hyperparameter search to fine-tune LR, \#epochs \& \#layers for each dataset (again 80/20 split).
Their optimal parameters are given in the appendix, in Table~\ref{tab:best_hparams}.
We evaluate the best models on the official test set.
Our results regarding Hits@1 (H@1), Hits@10 (H@10), mean rank (MR) and mean reciprocal rank (MRR) are summarised in Table~\ref{tab:result}.

\paragraph{Node Embedding Initialization}
Comparing the columns of Table~\ref{tab:result} we can observe the influence of the node embedding initialization.
Using the settings from the authors' code, i.e. not using weights, a choosing a variance of $n^{-1/2}$ actually results in inferior performance in terms of H@1, as compared to use a standard normal distribution.
These findings are consistent across datasets.

\paragraph{Convolution Weights}
The first column of Table~\ref{tab:result} corresponds to the weight usage and initialization settings used in the code for GCN-Align.
We achieve slightly better results than published in \cite{DBLP:conf/emnlp/WangLLZ18}, which we attribute to a more exhaustive parameter search.
Interestingly, all best configurations use Adam optimizer instead of SGD.
Adding convolution weights degrades the performance across all datasets and subsets thereof but one as witnessed by comparing the first two columns with the last two columns.

\section{Conclusion}
In this work, we reported our experiences when implementing the Knowledge Graph alignment method GCN-Align. We pointed at important differences between the model described in the paper and the actual implementation and quantified their effects in the ablation study.
For future work, we plan to include other methods for entity alignments in our framework.
\section*{Acknowledgements}
This work has been funded by the German Federal Ministry of Education and Research (BMBF) under Grant No. 01IS18036A and by the Bavarian Ministry for Economic Affairs, Infrastructure, Transport and Technology through the Center for Analytics-Data-Applications (ADA-Center) within the framework of “BAYERN DIGITAL II”. The authors of this work take full responsibilities for its content.

\bibliographystyle{splncs04}