Numerical Computation of Weil-Peterson Geodesics in the Universal Teichm\"uller Space
Matt Feiszli and Akil Narayan
We propose an optimization algorithm for computing geodesics on the universal Teichm\"uller space $T(1)$ in the Weil-Petersson ($W P$) metric. Another realization for $T(1)$ is the space of planar shapes, modulo translation and scale, and thus our algorithm addresses a fundamental problem in computer vision: compute the distance between two given shapes. The identification of smooth shapes with elements on $T(1)$ allows us to represent a shape as a diffeomorphism on $S^1$. Then given two diffeomorphisms on $S^1$ (i.e., two shapes we want connect with a flow), we formulate a discretized $W P$ energy and the resulting problem is a boundary-value minimization problem. We numerically solve this problem, providing several examples of geodesic flow on the space of shapes, and verifying mathematical properties of $T(1)$. Our algorithm is more general than the application here in the sense that it can be used to compute geodesics on any other Riemannian manifold.
minimization.tex
% Algorithm detail for finding shape-connecting geodesic.
\begin{algorithmic}

\STATE Input: Initial and final particle positions $\V{q}^0$ and $\V{q}^{T+1}$
\STATE Initialize: $\V{v}^t = \V{q}^{T+1} - \V{q}^0$ for all $t$
\WHILE{Not converged}
  \STATE Compute norm, momenta, and kernel coefficients from \eqref{eq:gsvd-solution}
  \STATE Compute standard unconstrained gradient $\nabla E$ from \eqref{eq:standard-gradient-formula} and \eqref{eq:wpnorm-variation-particle}
  \STATE Form projected natural gradient $\widetilde{\nabla E}_n$ from \eqref{eq:operator-Y-definition} and \eqref{eq:discrete-path-metric}.
  %\eqref{eq:standard-gradient-projection}, \eqref{eq:natural-gradient}, and \eqref{eq:natural-gradient-projection}.
  \STATE Find $\varepsilon$ such that $E(V - \varepsilon \widetilde{\nabla E}_n) < E(V)$
  \STATE Update $V \gets V - \varepsilon \widetilde{\nabla E}_n$, update particles using \eqref{eq:particle-update}
\ENDWHILE
\STATE Output: Velocity field $V$
\end{algorithmic}

std-preamble-beamer.tex
\usepackage{amsmath, amsthm}
\usepackage{color}

%\usepackage{pictex}


%\newtheorem{theorem}{Theorem}
%\newtheorem{corollary}[theorem]{Corollary}
%\newtheorem{definition}{Definition}
%\newtheorem{lemma}{Lemma}
%\newtheorem{proposition}{Proposition}
%\newtheorem{example}{Example}

\newcommand{\Disc}{\mathbb{D}}
\newcommand{\Chat}{\hat{\mathbb{C}}}
\newcommand{\Dstar}{\Disc^{*}}
\newcommand{\Hpl}{\mathbb{C}_{+}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\HS}{\mathbb{H}}
\newcommand{\E}{\mathbb{E}}
\newcommand{\dHS}{\partial \HS^{3}}
\newcommand{\B}{\mathbb{B}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\N}{\mathbb{N}}

\newcommand{\Area}[1]{\mathrm{Area}[#1]}
\newcommand{\BigArea}[1]{\mathrm{Area}\Big[#1\Big]}
\newcommand{\Length}[1]{\mathrm{L}[#1]}
\newcommand{\dist}[2]{\mathrm{d}(#1, #2)}

\newcommand{\ud}{\mathrm{d}}
\newcommand{\argu}[1]{\mathrm{arg}(#1)}

\newcommand{\poincare}{Poincar\'e\ }
\newcommand{\mob}{M\"obius\ }
\newcommand{\mobius}{M\"obius\ }
\newcommand{\teichmuller}{Teichm\"uller\ }

\newcommand{\mplus}{M_{r}^{+}}
\newcommand{\mminus}{M_{r}^{-}}

\renewcommand{\Re}{\mathrm{Re}}
\renewcommand{\Im}{\mathrm{Im}}

\def\Circle{\mathbb{T}}
\def\domega{\partial\Omega}
\def\ddisc{\partial\Disc}
\def\Dome{S_{\Omega}}
\def\rhodome{\rho_{\Dome}}
\def\partialz{\frac{\partial}{\partial z}}
\def\partialt{\frac{\partial}{\partial t}}
\def\H{\mathcal{H}}

\def\wo{w_{o}}
\def\zo{z_{o}}
\def\z1{z_{1}}

\def\ddt{\frac{\partial}{\partial t}}

\def\dbar{\bar \partial}

% Geometric groups
\def\Diff{\mathbf{Diff}}
\def\Vec{\mathbf{Vec}}
\def\Homeo{\mathbf{Homeo}}

\newcommand{\myvec}[1]
{
   \begin{pmatrix} #1 \end{pmatrix}
}

%universal teichmuller space
\def\uts{\mathbf{Teich}(\Delta)}
\def\ubb{L^{\infty}(\Delta)_{1}}     % unit ball of beltrami coefficients
\def\bdb{L^{\infty}(\Delta)}            % bdd, meas beltramis on the disk

% differentials:
\newcommand{\mathd}{\mathrm{d}}
\newcommand{\dx}[1]{\mathd #1}
\newcommand{\ddx}[1]{\frac{\mathd}{\mathd #1}}
\newcommand{\ppx}[1]{\frac{\partial}{\partial #1}}
\newcommand{\dfdx}[2]{\frac{\mathd #1}{\mathd #2}}
\newcommand{\pfpx}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\red}[1]{{\bf \textcolor{red}{#1}}}
\newcommand{\blue}[1]{{\bf \textcolor{blue}{#1}}}

std-preamble.tex
\usepackage{amsmath, amsthm}
%\usepackage{pictex}


\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}{Lemma}
\newtheorem{proposition}{Proposition}
\newtheorem{example}{Example}

\theoremstyle{remark}
\newtheorem*{rem}{Remark}

\theoremstyle{definition}
\newtheorem{definition}{Definition}

%\input{/Users/mattfeiszli/Documents/synch/tex/std-preamble-beamer.tex}
%\input{std-preamble-beamer.tex}
wp-geodesics-matt2.tex
\documentclass[11pt]{amsart}
%\usepackage{geometry}                % See geometry.pdf to learn the layout options. There are lots.
\usepackage[left=1.25in,right=1.25in,top=1in]{geometry}
\geometry{letterpaper}                   % ... or a4paper or a5paper or ... 
%\geometry{landscape}                % Activate for for rotated page geometry
%\usepackage[parfill]{parskip}    % Activate to begin paragraphs with an empty line rather than an indent
\usepackage{graphicx,verbatim,array}
\usepackage{graphics,url}
\usepackage{amssymb}
%\usepackage{epstopdf}
%\DeclareGraphicsRule{.tif}{png}{.png}{`convert #1 `dirname #1`/`basename #1 .tif`.png}
\usepackage{amsmath,amsthm}

\input{std-preamble.tex}
\input{std-preamble-beamer.tex}

\usepackage[textwidth=1.25in,textsize=tiny]{todonotes}
\usepackage[footnotesize,bf]{caption}

%\graphicspath{{../ch-bdiff/}{../ch-weld/}{../../preso/resources/}}

\title[Computation of Weil-Peterson Geodesics]{Numerical Computation of Weil-Peterson Geodesics in the Universal
Teichm\"uller Space}
\author[Feiszli and Narayan]{Matt Feiszli and Akil Narayan}
%\date{}                                           % Activate to display a given date or no date


\newcommand{\tN}{\mathcal{N}}
\newcommand{\tNH}{\mathcal{N(\mathcal{H})}}

\newcommand{\vp}{v_{\perp}}
\newcommand{\vr}{v_{r}}
\newcommand{\jpoly}[3]{p_{#1}^{(#2,#3)}}
\newcommand{\wjacobi}[2]{\omega^{(#1, #2)}}
\newcommand{\acos}{\textrm{acos}\,}
\newcommand{\D}{\mathbb{D}}
\newcommand{\proj}{\mathrm{proj}}
\renewcommand{\ker}{\mathrm{ker}}
\newcommand{\ran}{\mathrm{ran}}
\newcommand{\mspan}{\mathrm{span}}
\newcommand{\rank}{\mathrm{rank}}
\newcommand{\pinv}{\dagger}

\newcommand{\dfndx}[3]{\frac{\mathrm{d}^{#2}#1}{\mathrm{d}#3^{#2}}}
\newcommand{\px}[1]{\partial #1}
%\newcommand{\ppx}[1]{\frac{\partial}{\partial #1}}
%\newcommand{\pfpx}[2]{\frac{\partial #1}{\partial #2}}
%\newcommand{\dfdx}[2]{\frac{\mathrm{d} #1}{\mathrm{d} #2}}
%\newcommand{\ddx}[1]{\frac{\mathrm{d}}{\mathrm{d}#1}}

\newcommand{\UTS}{PSL_{2}(\R) \backslash \Diff(S^{1})}

\newcommand{\lmpt}{\mathcal{Q}}

\newcommand{\lm}{\mathcal{L}^{M}}

\newcommand{\V}[1]{\mathbf{#1}}
%\newcommand{\V}[1]{#1}

\newcommand{\aknote}[1]{\red{#1}}
\newcommand{\mfnote}[1]{\blue{#1}}
%\newcommand{\mfnote}[1]{\vskip 5pt \noindent {\color{blue}{\bf Matt}:} {\tt #1}\vskip 5pt}
%\newcommand{\aknote}[1]{\linespread{0.6}\todo[color=red!20]{#1}\linespread{1}}
% uncomment to 'remove' notes:
%\renewcommand{\aknote}[1]{}
\renewcommand{\mfnote}[1]{#1}

\usepackage{ulem}
\usepackage{algorithm,algorithmic}

\begin{document}
\begin{abstract}
We propose an optimization algorithm for computing geodesics on the universal Teichm\"uller space $T(1)$ in the Weil-Petersson ($W P$) metric. Another realization for $T(1)$ is the space of planar shapes, modulo translation and scale, and thus our algorithm addresses a fundamental problem in computer vision: compute the distance between two given shapes. The identification of smooth shapes with elements on $T(1)$ allows us to represent a shape as a diffeomorphism on $S^1$. Then given two diffeomorphisms on $S^1$ (i.e., two shapes we want connect with a flow), we formulate a discretized $W P$ energy and the resulting problem is a boundary-value minimization problem. We numerically solve this problem, providing several examples of geodesic flow on the space of shapes, and verifying mathematical properties of $T(1)$. Our algorithm is more general than the application here in the sense that it can be used to compute geodesics on any other Riemannian manifold.
\end{abstract}
\maketitle
%We begin by summarizing some computational issues in the numerical solution of
%the geodesic equation for the Weil-Peterson metric on the universal
%Teichm\"uller space and introducing the relevant background material.

\section{Introduction and Background}

Representation and comparison of shapes is a central problem in computer vision.  In the past several decades, many approaches to represent, compare, and classify shapes have been presented (See %\cite{veltkamp-shape}, 
\cite{veltkamp_shape_2001}, %\cite{mumford-perception} 
\cite{mumford_pattern_2002} for review and discussion).  The space of 2D shapes is inherently nonlinear; this poses fundamental difficulties in computer vision when attempting object recognition and statistics.  In \cite{sharon_2d-shape_2006}, Mumford and Sharon describe a construction based on conformal mapping which makes the space of simple closed plane curves into a Riemannian metric space.  The space itself is in fact the universal Teichm\"uller space with the Weil-Peterson metric; in this paper, we describe a numerical solver for geodesics in this space. In Sections \ref{sec:intro-welding} and \ref{sec:intro-wp} we describe the optimization problem we wish to solve: minimization of an energy functional with boundary value constraints. We note that a minimization algorithm proposed in \cite{sharon_2d-shape_2006} was applied to only relatively simple shapes because of numerical difficulties. In this work we aim to apply our algorithm to more general, complicated shapes. Our work is competitive with a recent approach based on shooting for the analogous boundary value problem \cite{kushnarev-2012}.

The outline of this paper is as follows: Section \ref{sec:discrete} describes the mathematics of computations for the Weil-Peterson metric using discrete samples of a velocity field. Section \ref{sec:computations} discusses the fully discrete algorithm along with our technique for satisfying the boundary constraints. Finally, Section \ref{sec:results} presents numerical results that illustrate the effectiveness of the algorithm.

%\section{Background}

\subsection{Conformal Welding}\label{sec:intro-welding}

The construction in \cite{sharon_2d-shape_2006} is classical in Teichm\"uller theory.  Given a simply-connected planar region $\Omega$ bounded by a smooth Jordan curve (this is our definition of a ``shape''), one constructs a pair of conformal maps $\Phi_{+}: \Delta \to \Omega_{+}$ and $\Phi_{-}: \Delta \to \Omega_{-}$, which map the exterior and interior of the unit disk to the exterior and interior of $\Omega$, respectively.  The exterior map is normalized to fix $\infty$ and have real derivative there.  The interior map is only defined up to right multiplication by the three-parameter M\"obius group $PSL_{2}(\R)$ of conformal self-maps of the unit disk.  Both $\Phi_{+},\Phi_{-}$ extend continuously to the boundary $S^{1}$, and their composition $\Phi_{-}^{-1} \circ \Phi_{+}$, restricted to the boundary, is a map $\Psi: S^{1} \to S^{1}$.  


A remarkable result (see %\cite{ahlfors-quasiconformal
\cite{ahlfors_lectures_2006}, for example) is that this result is almost an isomorphism: the space of shapes, modulo translation and scale, is isomorphic to the group $\Diff(S^{1})$, modulo conformal self-maps of the disk.  This provides an elegant way of making the space of shapes into a metric space:  we take an element of the coset space $\UTS$, i.e. an equivalence class of diffeomorphisms of the circle, as a representation of a shape.  This space is known as the universal Teichm\"uller space and was initially studied in the context of Riemann surfaces %\cite{ahlfors-quasiconformal, hubbard-book}; 
\cite{ahlfors_lectures_2006,hubbard_teichmuller_2006}
it also arises  in string theory %\cite{bowick-rajeev, nag}
\cite{bowick_string_1987,nag_rm_1990}.  In \cite{sharon_2d-shape_2006} elements of $\UTS$ are called ``fingerprints''; in the mathematical literature, elements of the broader class of quasisymmetric homeomorphisms of $S^{1}$ are commonly known as ``welding maps'' and each has an associated quasicircle.  We illustrate an example of a equivalence class of welding maps in Figure \ref{fig:welding-maps} for a simple shape.

\begin{figure}
  \begin{center}
    \resizebox{1.0\textwidth}{!}{
    \includegraphics{welding-maps}
    }
  \end{center}
  \caption{A triangular region $\Omega$ (left) transformed into an equivalence class of welding maps $\phi$ (right). Four members of the equivalence class are shown. Four segments of $\partial \Omega$ are drawn with differing linestyles to highlight their corresponding segment on each of the welding maps.}
  \label{fig:welding-maps} 
\end{figure}

\subsection{The Weil-Peterson Riemannian metric}\label{sec:intro-wp}

The Lie algebra to $\UTS$ is the space of vector fields on the circle modulo the subspace spanned by $1, \cos$ and $\sin$.
If $v(\theta) = \sum_{n=-\infty}^{\infty} a_{n} e^{in\theta} \frac{\partial}{\partial \theta}$ is a real-valued vector field on the circle, the Weil-Peterson norm is 
\begin{align}\label{eq:wp-norm-fourier}
\| v \|_{WP}^{2} = \sum_{n=2}^{\infty} (n^{3}-n) |a_{n}|^{2} = \frac{1}{2} \left( \left\| v^{(3/2)} \right\|^2 - \left\|v^{(1/2)}\right\|^2 \right)
\end{align}
where $\left\|\cdot\right\|$ denotes the $L^2$ norm.  \footnote{This formula, while explicit, is troublesome in numerical computations; the boundary values of conformal maps have so much high-frequency content, even for curves with real-analytic boundaries, that the number of Fourier coefficients required for for an accurate global representation is prohibitive.}  This may be rewritten as
\begin{align*}
\| v \|_{WP}^{2} = \langle Lv, v \rangle
\end{align*}
for the positive semidefinite, self-adjoint operator
\begin{align}\label{eq:operator-L-definition}
L = -\mathcal{H}(\partial^{3} - \partial)
\end{align}
where $\mathcal{H}$ is the Hilbert transform.  Note that $L$ has a kernel which is exactly the span of $1, \cos$ and $\sin$.  As mentioned above, the corresponding vector fields in the span of $\left\{ \sin\theta \frac{\partial}{\partial \theta}, \cos\theta \frac{\partial}{\partial \theta}, \frac{\partial}{\partial \theta}  \right\}$  are infinitesimal M\"obius maps; i.e. they span the Lie algebra $psl_{2}(\R)$ of the M\"obius group $PSL_{2}(\R)$.  Hence, our norm is indeed a norm on the quotient space.

By right-translations we extend this norm to the entire space; that is, if $\phi(\theta, t) \equiv \phi_{t}(\theta)$ is a curve in $\UTS$, then any any time $t$ we pull back the velocity $\dot \phi_{t}$ under the derivative of right translation to obtain
\begin{align*}
v_{t}(\theta) = \dot \phi_{t} \circ \phi_{t}^{-1}
\end{align*}
The length of a curve $\phi:[0,1] \times S^{1} \to S^{1}$ in $\Diff(S^{1})$ is obtained by integrating
\begin{align*}
L[\phi] = \int_{0}^{1} \| \dot \phi_{t} \circ \phi_{t}^{-1} \|_{WP} dt
\end{align*}
The resulting Riemannian metric is known as the Weil-Peterson metric.  The classical definition was in terms of Beltrami differentials of quasiconformal self-maps of the hyperbolic plane and was originally studied in the context of Teichm\"uller theory.  The formulation we present here occurs naturally in string theory where it was discovered as the unique K\"ahler metric on $\UTS$ %\cite{segal-cft}
\cite{segal_definition_2004}.  See %\cite{nag} 
\cite{nag_rm_1990} for a thorough discussion of how this K\"ahler metric agrees with the classical WP metric.

It is a fact from variational calculus that minimizing length is equivalent to minimizing energy: the path of least energy is a constant-speed parametrization of the path of least length.  Thus given a pair of points $h_{0}, h_{1} \in \Diff(S^{1})$ we find the geodesic by finding the path $\phi$ of minimum energy
\begin{align}\label{eq:energy-definition}
E[\phi] = \int_{0}^{1} \| \dot \phi_{s} \circ \phi_{s}^{-1} \|_{WP}^{2} ds,
\hskip 10pt \phi_0 = h_0, \hskip 10pt \phi_1 = h_1
\end{align}
Finding this minimizing $\phi$ is the problem we are concerned with.  It is a fact that all sectional curvatures of the WP metric are negative %\cite{takhtajan-teo-wp}
\cite{teo_weil-petersson_2006}, which in finite dimensions would imply that minimizing geodesics are unique.  While the existence and uniqueness of minimizing geodesics is a subtle issue in infinite dimensions, both questions were recently answered in the affirmative in %\cite{fgb-wp-unique}
\cite{gay-balmaz_geometry_2009}:  Any two points in $\UTS$ are joined by a unique geodesic, and further, solutions to the geodesic equation exist for all time and the resulting welding maps are Sobolev $H^{s}$ for any $s < 3/2$, where the inequality is sharp.  




It is an essential fact that the WP metric we have described here is, by construction, \textit{right-invariant}.  That is, given any $h_{0}, h_{1}$, define the map $g$ by $h_{0} = g \circ h_{1}$ and let $\phi(\theta, s)$ be the shortest path from $g$ to the identity.  Right-invariance of the metric means that the length of $\phi(\theta, s)$ is the same as length of $\phi(h_{1}(\theta), s)$, which immediately implies that the shortest path from $h_{0}$ to $h_{1}$ is in fact $\phi(h_{1}(\theta), s)$.  Hence, computing geodesics between two welding maps reduces to the case where one of the maps is the identity.   Therefore in what follows we consider only paths ending at the identity.




\section{The Discrete Problem}\label{sec:discrete}

%The reason to $W P$ norms is to make computation of path lengths feasible. Again
%we opt to instead use the energy \eqref{eq:energy-definition}. 

We now discretize the problem \eqref{eq:energy-definition} in the case where $h_{0}: S^{1} \to S^{1}$ is an arbitrary welding map and $h_{1}$ is  the identity map.  This is a boundary-value minimization problem: find $v_s = \dot{\phi}_s \circ \phi_s^{-1}$ subject to initial and
terminal conditions $\phi_0 = h_{0}$ and $\phi_1 = \mathrm{Id}$ such that $\int \|v_s
\|_{WP}^2 \dx{s}$ is minimized. Becuase we eventually adopt a Lagrangian framework we introduce $q \equiv \phi$, and given a velocity field $v$ then 
\begin{align*}
  q(t,x) = q(0,x) + \int_0^t v(\xi, q(\xi,x)) \dx{\xi}
\end{align*}
for all $0 \leq t \leq 1$.

We track $M$ particles on the range of $q$ as they evolve in $t$.
For some fixed $x_m$ on $S^1$, $q(t, x_m)$ are the locations of these particles. The velocity field $v(t, \cdot)$ at the locations $q(t,x_m)$ is likewise stored.
We store these velocity field values and recover the particle positions $q(t,x_m)$ by integrating rather than directly storing the particle positions. 
%(More formally, we choose $M$ values $q_m(1)$ on $S^1$ and track the particle locations
%$q_m(t) = q_t(q_m(1))$ along evolution in $t$.)  The velocity
%field $v_s$ at each point in time will be
%represented by discrete point-evaluations at the locations $q_t(q_m(1))$; 
In the next subsection we will describe how we compute the $WP$ norm at each time. We delay discretizing the $t$ variable until Section \ref{sec:computations}.

%\footnote{The sensitivity of the WP norm to noise is difficult to overstate.  We found that... }. 

\subsection{Computing the WP Norm}




%We now seek to extend our result to all vector fields on $V$. We first make use
%of the following simple fact:
%\begin{lemma}\label{lemma:kernel-totality}
%  For any configuration $\mathcal{L}^3$ of three distinct points on $S^1$, the
%  subspace $\ker_V L = \mspan \{1, \sin, \cos\}$ is total for Dirac interpolation.
%\end{lemma}
%The important consequence of the above is that the coordinate map taking any
%basis of $\ker_V L$ to Dirac evaluations on $\mathcal{L}^M$ can be represented
%by an $M \times 3$ matrix of full rank.
%\begin{lemma}\label{lemma:v-minimizer}
%Let $G$ be an $M\times M$ positive matrix, $B$ an $M \times K$ full-rank matrix
%with $K < M$. Then given some $w \in R^M$, the $G^{-1}$-orthogonal projection of
%$w$ onto $\ker(B^T) = \ran (B)^\perp$ is given by 
%\begin{align*}
%  v = P_G^B w \doteq (I - G^{-1} B (B^T G^{-1} B)^{-1} B^T) w
%\end{align*}
%\end{lemma}
%The proof of the above is an exercise in Lagrange multipliers. We are concerned
%with the case when $G^{-1}$ is a Riemannian metric; however, in our case we only want
%to use the metric for a subspace of $\R^M$ and is not intended for application
%on all $\R^M$. Therefore, we seek to start on the subspace $\ran(B)^\perp$ and
%extend the metric to all $\R^M$, assigning zero norm to any needed extensions.
%In other words we wish to identify $P^L_G$ as the differential of a Riemannian
%submersion.\aknote{Does this sentence make sense?}
%\begin{lemma}
%Assume that $G^{-1}$ is a metric on $\ker(B^T)$, where $B$ is a full-rank $M
%\times K$ matrix with $K < M$. Then the norm-minimizing extension of $G^{-1}$ to
%all $\R^M$ is given by  
%\begin{align*}
%  H = P_G^B G^{-1}.
%\end{align*}
%\end{lemma}
%\begin{proof}
%Given $\V{w}_1, \V{w}_2 \in \R^M$, apply the $G^{-1}$ metric to $P^B_G \V{w}_1$
%and $P^B_G \V{w}_2$.
%\end{proof}
%By construction $H$ has rank $M-K$. This immediately yields the following.
%\begin{corollary}
%\label{corollary-green-gradient}
%In the quotient space $\lm$, the gradient of the WP norm at $\Theta$ is
%\begin{align*}
%\frac{\partial}{\partial v_{j}} \| \widetilde v \|_{WP}^{2} = \left[ 2v^{T}H \right]_{j} 
%\end{align*}
%\begin{proof}
%This is just the usual formula for the derivative of a quadratic form.
%\end{proof}
%\end{corollary}
%Corollary \ref{corollary-green-gradient} is the first step to providing a
%gradient for our optimization routine. 
%%at each timestep yields a
%%vector field on $[0,1] \times S^{1}$ which we can then project onto
%%$T_{\phi}P_{h_{0}}$ to perform our optimization.
%
%There is another way to see that the previous arguments yield the metric we
%desire. If $v$ is a velocity field
%corresponding to velocity field, then the momentum $L v$ is a sum of Dirac
%delta functions. Since our manifold $\mathcal{L}^M$ is defined by a finite
%number of point evaluations, this is equivalent to describing velocity fields by
%momenta that are finite sums of Dirac delta functions. 
%
%The question then is how does a momentum field described by an impluse train of
%Dirac functions behave on this metric?  The Euler-Poincare differential equation
%describes the flow of geodesics on Teichm\"{u}ller space. It is well-known that
%self-similar solutions to this equation, known as {\em teichons}, exist. These
%teichons can take the form of an impulse train. However, in order for the
%Euler-Poincare equation to be well-posed, the coefficients of the impulse train
%must satisfy certain constraints. Namely, the momentum field must be orthogonal
%to $\ker_V L$ -- this is equivalent to asking for orthogonality to $\ran(B)$.
%The natural metric on the momentum field is $G$, so we would like the
%$G$-orthogonal projection onto $\ran(B)$. This turns out to be $P^B_G$.
%Therefore, the metric on the $M$-dimensional space of momenta strengths on
%$\mathcal{L}^M$ is given by $(P^G_B)^T G P^G_B = G P^G_B$, and on the $M$-dimensional space of
%velocity strengths on $\mathcal{L}^M$, it is given by $P^B_G G^{-1}$.
%
%Having computed the horizontal lifts satisfying the desired constraints, we will
%find it useful to present an explicit formula.
%\begin{corollary}\label{cor:v-field-wp-min}
%Fix $\Theta$ and let the matrix $G$ be as in Proposition
%\ref{prop:v-horizontal-lift}. For any $\V{v} \in T_\Theta \mathcal{L}^M$, 
%\begin{align}\label{eq:v-discrete-minimizer}
%  \| v \|^2_{WP(\lmpt)} = \min\{ (\V{v} - \V{w})^T G^{-1} (\V{v}-\V{w}) \}
%\end{align}
%where the minimization is performed over all $\V{w} \in \ran(B)$.
%The horizontal lift $\widetilde{v}$ of $\V{v}$ to $V$ is given by
%\begin{align*}
%  \V{w} &= (B^T G^{-1} B)^{-1} B^T \V{v} \\
%  \widetilde{v} &= \widetilde{w} + \sum_{j=1}^M \left( G^{-1} (\V{v}-G^{-1} B \V{w})
%  \right)_j G(\cdot - \theta_j), \\
%  \widetilde{w} &= \sum_{j=1}^3 \V{w}_j
%  b_j(\cdot), \\
%\end{align*}
%where $\{b_j\}_{j=1}^3$ is the basis for $\ker_V L$ used to construct $B$:
%$B_{m,k} = b_k(\theta_m)$. We will hereafter refer to $\V{p} = H\V{v}$ as the
%momenta coefficients.
%\end{corollary}
%
%\noindent Thus for any $\V{v} \in T_\Theta \mathcal{L}^M$ we can identify a
%corresponding element $\widetilde{v}$ from $T_h\mathbf{Diff}(S^1)_{hor}$; the
%matrix $H = G^{-1} P^B_G$ also allows us to define an appropriate cotangent vector:
%\begin{align*}
%  v \mapsto L \widetilde{v} = \sum_{m=1}^M (H\V{v})_m \delta'_{\theta_m},
%\end{align*}
%where $\delta'_\theta$ is the Dirac measure $L^2$-orthogonalized against
%$\ker_V L$ centered at $\theta \in S^1$.

We first consider the problem of computing the norm of a vector field on $S^1$.  Subsequent sections will consider computation of path energy and its gradient to find geodesics.   

Our discretization provides us with information about the velocity field only at a discrete set of points.  Our first task is to interpolate this velocity field to a field defined on all of $S^{1}$ and then compute the norm of this lift. The computation proceeds in 3 steps.  We first collect some standard results which let us do this when the metric has no kernel.  Next we compute the optimal interpolant and its norm when the metric has a kernel; this is done by projecting the interpolation data onto the orthogonal complement of the kernel and computing the norm of the lift.  Finally, we generalize the problem to allow an arbitrary basis for the interpolating functions; this allows for more flexibility in the implementation and also suggests a more stable method of computation, which we detail later.

\subsubsection{Lifting the velocity field into the orthogonal complement of the kernel}\label{sec:lifting-v}

We need to extend a vector field defined on a finite subset of $S^{1}$ to a vector field on all of $S^{1}$.  Here we describe a method for finding the smoothest possible extension using Green's functions for the operator $L$.  In this section we consider only vector fields orthogonal to the kernel; we extend to more general cases next.

Consider the space of all vector fields on $S^1$, and a subspace on which
$\|\cdot\|_{W P}$ is a proper norm:
\begin{align*}
  V &= \{v: \|v \|_{W P} < \infty\}, & \widetilde{V} = \{v \in V: v \perp 
  \ker_V L\},
\end{align*}
where, explicitly, $\ker_V L = \mathrm{span}\{1, \sin, \cos\}$.
Let $\lm$ be the $M$-manifold of configurations of $M$ points $q_{1} <
q_{2} < ... < q_{M}$ on the unit circle and write $\lmpt = \left\{
q_{1}, q_{2}, ..., q_{M}  \right\}$. The tangent space
$T_{\lmpt}\lm$ is the $M$-dimensional space of vector fields $v = (v_{1},
v_{2}, ..., v_{M})$ supported on the set $\lmpt$.  We consider all possible
extensions of $v$ to vector fields defined on the entire circle, and induce a
norm on $T_{\lmpt}\lm$ by 
\begin{align}\label{eq:v-global-minimizer}
\| v \|_{WP(\lmpt)} \equiv \inf \left\{ \| \widetilde v \|_{WP} : \widetilde{v}(q_{m}) = v_{m}, 1 \le m \le M  \right\}
\end{align}
Conversely, given the basepoint $\lmpt$, consider the evaluation map taking a global vector field $\widetilde v$, defined on all of $S^{1}$, to $T_{\lmpt}\lm$ by
\begin{align*}
\widetilde v \to \left( \widetilde v(q_{1}), \widetilde v(q_{2}), ..., \widetilde v(q_{M})   \right)
\end{align*}
This map induces a splitting of any $\widetilde v$ into two components.  The set of vector fields for which $\widetilde v(q_{m}) = 0$ for all $1 \le m \le M$ is the \textit{vertical subspace} (of the \textit{Lie algebra} of $\Diff(S^{1})$: the set of smooth vector fields on $S^{1}$), and its $WP$-orthogonal complement is the \textit{horizontal subspace}; the minimizing extension in \eqref{eq:v-global-minimizer} of $v \in T_{\lmpt}\lm$ is known as a \textit{horizontal lift}.


As the next proposition shows, the horizontal subspace is an $M$-dimensional subspace spanned by translates of Green's function for the operator $L = -\mathcal{H}(\partial^{3} - \partial)$.  This Green's function is known explicitly, see \cite{kushnarev_teichons_2009}:
\begin{align*}
G(\theta) = 2 \sum_{n=2}^{\infty} \frac{\cos(n\theta)}{(n^{3}-n)} = (1-\cos \theta) \log\left[ 2(1-\cos \theta)  \right] + \frac{3}{2}\cos \theta -1
\end{align*}
The following results are standard facts which may be proven using the reproducing kernel property of Green's function and the fact that the vector fields $G(\theta-q_{m})\frac{\partial}{\partial\theta}$ lie in the horizontal subspace.
\begin{proposition}[Computing Horizontal lifts]
\label{prop:horizontal-lift-nokernel}
\ \newline
\begin{enumerate}
\item The set of vector fields
\begin{align*}
\left\{ G(\theta - q_{m}) \frac{\partial}{\partial\theta} \right\}_{m=1}^{M}
\end{align*}
is a basis for the component of the
horizontal subspace in $\widetilde{V}$.
%\begin{proof}
%By the reproducing kernel property of Green's function, we know that 
%\begin{align*}
%\langle v, G(\cdot - \theta_{j}) \rangle_{WP} = v(\theta_{j})
%\end{align*}
%for all $j = 1,2,...,M$ and any vector $v \in T_{h}\Diff(S^{1}) \cap \widetilde{V}$.
%In particular, if $v$ is vertical we have
%\begin{align*}
%\langle v, G(\cdot - \theta_{j} )\rangle_{WP} = v(\theta_{j}) = 0
%\end{align*}
%so $G\frac{\partial}{\partial\theta}$ is orthogonal to the vertical subspace;
%i.e. $v \in T_{h}\Diff(S^{1})_{hor}\cap \widetilde{V}$.  Further, these vector
%fields are linearly independent.  Indeed, if we have $\{ a_{j} \}_{j=1}^{M}
%\subset \R$ so that
%\begin{align*}
%\sum_{j=1}^{M}a_{j} G(\theta - \theta_{j}) = 0
%\end{align*}
%then for any $v$ we have
%\begin{align*}
%0 = \langle v, \sum_{j=1}^{M}a_{j} G(\theta - \theta_{j}) \rangle_{WP} = \sum_{j=1}^{M}a_{j}v(\theta_{j})
%\end{align*}
%and hence the only possible choice is that all $a_{j} = 0$.
%\end{proof}
\item The horizontal lift on $\widetilde{V}$ of the tangent vector $(v_{1}, ..., v_{M})
\in T_{\lmpt}\lm$ is the vector field $\widetilde v(\theta)
\frac{\partial}{\partial \theta}$ for
\begin{align*}
\widetilde v(\theta) = \sum_{i=1}^{M}G(\theta - q_{i})p_{i} \frac{}{}\quad where \quad p_{i} = \sum_{j=1}^{M}G^{-1}_{ij}v_{j}
\end{align*}
and $G_{ij}$ is the positive definite symmetric matrix
\begin{align}\label{eq:G-matrix}
G_{ij} = \left\langle G(\cdot - q_i), G(\cdot - q_j) \right\rangle_{W P} = G(q_{i} - q_{j} )
\end{align}
Further,
\begin{align*}
\| \widetilde v \|_{WP}^{2} = \sum_{i,j}G^{-1}_{ij} v_{i}v_{j}
\end{align*}
\end{enumerate}
%\begin{proof}
%The basis identified in Proposition \ref{prop:vector-field-basis} allows us to translate
%this into a linear system, and the solution is as presented above.
%Clearly the vector field $\widetilde v$ above is horizontal.  We first verify
%that it does the right thing at the $\theta_{j}$'s:
%\begin{align*}
%\widetilde v(\theta_{k}) &= \sum_{i}G(\theta_{k} - \theta_{i})u_{i} \\
%&= \sum_{i}G_{ki}\sum_{j}G^{-1}_{ij}v_{j} \\
%&= \sum_{i,j} G_{ki}G^{ij}v_{j} \\
%&= \sum_{j} (GG^{-1})_{kj} v_{j} \\
%&= v_{k}
%\end{align*}
%We now compute the WP length:
%\begin{align*}
%\| \widetilde v \|_{WP}^{2} &= \left\langle \sum_{i}G(\cdot-  \theta_{i})u_{i}, \sum_{j}G(\cdot - \theta_{j})u_{j}\right\rangle_{WP} \\
%&= \sum_{i,j}G(\theta_{i} - \theta_{j}) u_{i}u_{j}  \\
%&= \sum_{i,j}G_{ji} \sum_{k}G^{-1}_{ik}v_{k} \sum_{l} G^{-1}_{jl}v_{l} \\
%&= \sum_{j,k,l} I_{jk} G^{-1}_{kl} v_{k} v_{l} \\
%&= \sum_{k,l} G^{-1}_{kl}v_{k}v_{l}
%\end{align*}
%Renaming the dummy variables completes the proof.
%\end{proof}
\end{proposition}

\subsubsection{Lifting the velocity field}
Now we extend these results to compute horizontal lifts which take into account the kernel.  


\begin{lemma}\label{lemma:v-minimizer}
Let $G$ be an $M\times M$ positive matrix, $B$ an $M \times K$ full-rank matrix
with $K < M$. Then given some $v \in \R^M$, the projection of $v$ onto the
$G^{-1}$-orthogonal complement of the kernel of $B^T$ is given by 
\begin{align*}
  P_G^B v \doteq (I - B (B^T G^{-1} B)^{-1} B^TG^{-1}) v
\end{align*}
\end{lemma}
%\begin{proof}
%  We want to solve the minimization problem 
%  \begin{align*}
%    \min_{\V{v}} \V{v}^T G^{-1} \V{v} = \min_{\V{k}} (\V{w} - B \V{k})^T G^{-1} (\V{w} - B \V{k})
%  \end{align*}
%  Taking variations in $\V{k}$, we obtain the equation
%  \begin{align*}
%    (\V{w} - B \V{k})^T G^{-1} B = 0.
%  \end{align*}
%  Solving for $\V{k}$ yields the result.
%\end{proof}
In our case $B$ is an $M \times 3$ matrix containing point evaluations of a basis for the three-dimensinal kernel of the $W P$ norm.  Once we have subtracted out the contribution of the kernel, we may compute the lift and its norm.  The next result is a corollary of section \ref{sec:lifting-v} and the previous lemma.

\begin{corollary}\label{cor:v-field-wp-min}
Fix $\lmpt$ and let the matrix $G$ be as in Proposition
\ref{prop:horizontal-lift-nokernel}. For any $v \in T_\lmpt \lm$, 
\begin{align}\label{eq:v-discrete-minimizer}
  \| v \|^2_{WP(\lmpt)} = v^{T}G^{-1}P_{G}^{B}v  %\min\{ (\V{v} - \V{w})^T G^{-1} (\V{v}-\V{w}) \}
\end{align}
%where the minimization is performed over all $\V{w} \in \ran(B)$.
Further, the horizontal lift $\widetilde{v}$ of $v$ to $V$ is given by
\begin{align*}
  \widetilde{v} &= \sum_{j=1}^3 w_j
  b_j(\cdot)  + \sum_{j=1}^M \left( G^{-1} P_G^Bv
  \right)_j G(\cdot - q_j)
\end{align*}
where
\begin{align*}
  w &= (B^T G^{-1} B)^{-1} B^T G^{-1}v 
\end{align*}
and $\{b_j\}_{j=1}^3$ is the basis for $\ker_V L$ used to construct $B$:
$B_{m,k} = b_k(q_m)$. 
\end{corollary}
\begin{rem}
The vector $p = G^{-1}P_{G}^{B}v$ contains what are known as the momenta coefficients. Similarly,we call $w$ the kernel coefficients.
\end{rem}
The interpolant produced from Corollary \ref{cor:v-field-wp-min} is notably different than the horizontal lift from the quotient space in Proposition \ref{prop:horizontal-lift-nokernel}. This is exemplified in Figure \ref{fig:minimal-interpolation}.
\begin{figure}
  \begin{center}
    \resizebox{1.0\textwidth}{!}{
      \includegraphics{minimal-interpolation}
    }
  \end{center}
  \caption{Interpolation procedures from (a) Proposition \ref{prop:horizontal-lift-nokernel} on the quotient space $\widetilde{V}$, and from (b) Corollary \ref{cor:v-field-wp-min} on $V$. Left: the data is taken from the function $v(\theta) = \sin (2 \theta) + \cos(\theta)$ on 15 randomly distributed $\theta$ locations. Right: the data is taken as Gaussian random variable perturbations of function values from $v(\theta) = \sin (2 \theta)$ on 15 randomly distributed $\theta$ locations.}
  \label{fig:minimal-interpolation}
\end{figure}
%By construction $H$ has rank $M-K$. This immediately yields the following.
%\begin{corollary}
%\label{corollary-green-gradient}
%In the quotient space $\lm$, the gradient of the WP norm at $\Theta$ is
%\begin{align*}
%\frac{\partial}{\partial v_{j}} \| \widetilde v \|_{WP}^{2} = \left[ 2v^{T}H \right]_{j} 
%\end{align*}
%\begin{proof}
%This is just the usual formula for the derivative of a quadratic form.
%\end{proof}
%\end{corollary}
%Corollary \ref{corollary-green-gradient} is the first step to providing a
%gradient for our optimization routine. 
%%at each timestep yields a
%%vector field on $[0,1] \times S^{1}$ which we can then project onto
%%$T_{\phi}P_{h_{0}}$ to perform our optimization.
%
%There is another way to see that the previous arguments yield the metric we
%desire. If $v$ is a velocity field
%corresponding to velocity field, then the momentum $L v$ is a sum of Dirac
%delta functions. Since our manifold $\mathcal{L}^M$ is defined by a finite
%number of point evaluations, this is equivalent to describing velocity fields by
%momenta that are finite sums of Dirac delta functions. 
%
%The question then is how does a momentum field described by an impluse train of
%Dirac functions behave on this metric?  The Euler-Poincare differential equation
%describes the flow of geodesics on Teichm\"{u}ller space. It is well-known that
%self-similar solutions to this equation, known as {\em teichons}, exist. These
%teichons can take the form of an impulse train. However, in order for the
%Euler-Poincare equation to be well-posed, the coefficients of the impulse train
%must satisfy certain constraints. Namely, the momentum field must be orthogonal
%to $\ker_V L$ -- this is equivalent to asking for orthogonality to $\ran(B)$.
%The natural metric on the momentum field is $G$, so we would like the
%$G$-orthogonal projection onto $\ran(B)$. This turns out to be $P^B_G$.
%Therefore, the metric on the $M$-dimensional space of momenta strengths on
%$\mathcal{L}^M$ is given by $(P^G_B)^T G P^G_B = G P^G_B$, and on the $M$-dimensional space of
%velocity strengths on $\mathcal{L}^M$, it is given by $P^B_G G^{-1}$.
%
%Having computed the horizontal lifts satisfying the desired constraints, we will
%find it useful to present an explicit formula.
%\begin{corollary}\label{cor:v-field-wp-min}
%Fix $\Theta$ and let the matrix $G$ be as in Proposition
%\ref{prop:v-horizontal-lift}. For any $\V{v} \in T_\Theta \mathcal{L}^M$, 
%\begin{align}\label{eq:v-discrete-minimizer}
%  \| v \|^2_{WP(\lmpt)} = \min\{ (\V{v} - \V{w})^T G^{-1} (\V{v}-\V{w}) \}
%\end{align}
%where the minimization is performed over all $\V{w} \in \ran(B)$.
%The horizontal lift $\widetilde{v}$ of $\V{v}$ to $V$ is given by
%\begin{align*}
%  \V{w} &= (B^T G^{-1} B)^{-1} B^T \V{v} \\
%  \widetilde{v} &= \widetilde{w} + \sum_{j=1}^M \left( G^{-1} (\V{v}-G^{-1} B \V{w})
%  \right)_j G(\cdot - \theta_j), \\
%  \widetilde{w} &= \sum_{j=1}^3 \V{w}_j
%  b_j(\cdot), \\
%\end{align*}
%where $\{b_j\}_{j=1}^3$ is the basis for $\ker_V L$ used to construct $B$:
%$B_{m,k} = b_k(\theta_m)$. We will hereafter refer to $\V{p} = H\V{v}$ as the
%momenta coefficients.
%\end{corollary}


Up to this point, we have taken Green's functions centered at the interpolation nodes as 
as the basis for our interpolating space.  As discussed, this produces the norm-minimal lift.  While this optimality is nice, and is in fact what our current implementation ultimately uses, one might wish to choose a different space of interpolating functions.  For example, the first $N$ complex exponentials, or a wavelet basis adapted to the problem might be reasonable choices.  We address this now; it is only a small modification to the calculations above and admits a stable solution which we present later.


Let $\widetilde{V}_N$ be any $N$-dimensional subspace of $\widetilde{V}$ such that the interpolation problem for data collocated at $\lmpt$ is unisolvent.  If
$F = \left\{ f_n(\cdot) \right\}$ is a basis of $\widetilde{V}_N$, define the $M \times N$
matrix $\lambda_F$ to have entries 
\begin{align*}
%\lambda(F)_{m,n} = \langle f_n, \delta_{\theta_m}\rangle_{W P} = f_n(\theta_m)
\left(\lambda_F\right)_{m,n} = \langle f_n, G(\cdot-q_m)\rangle_{W P} = f_n(q_m)
\end{align*}
Let $G_F$ be the $N \times N$ Gram matrix for the $f_{n}$:
 \begin{align*}
\left(G_F\right)_{m,n} = \langle f_m, f_n\rangle_{WP}
\end{align*}
 For $v \in \R^M$ interpolation data on $\lmpt$, define $\widetilde{V}_N(v)$ as
the subset of functions from $\widetilde{V}_N \oplus \ker_V L$ that interpolate
to $v$. Then define
\begin{align*}
  \| v \|^2_{WP(\lmpt), \widetilde{V}_N} = %\min_{v \in \widetilde{V}_N (v)} \| v\|^2_{W P}
  \min_{v \in \widetilde{V}_N \oplus \ker_V L} \| v\|^2_{W P}
\end{align*}
The results from the previous section may be rephrased by replacing $G \gets G_{F}$
where $G_{F} = \lambda_F G_F^{-1} \lambda_F^T$.  Since $\widetilde{V}_N$ is total for interpolation, then $G_F$ is invertible. One finds that all formulas carry over essentially up to change-of-basis.  This, then, is our most general result for computing lifts, which we state in terms of an operator $L_{F}$.

\begin{proposition}\label{prop:metric-definition}
In the notation above, set
\begin{align*}
L_{F} = G_{F}^{-1} P^B_{G_{F}}
\end{align*}
Then 
\begin{align}\label{eq:other-subspace-norm}
  \| v \|_{WP(\lmpt), \widetilde{V}_N}^2 %&= \V{v}^T H \V{v} \\\nonumber
                %&=: \V{v}^T \left(Q^{-1}  - Q^{-1} B ( B^T Q^{-1} B)^{-1} B^T
                %Q^{-1}\right)\V{v} \\\nonumber
                &= v^T L_{F} v \\\nonumber
                &=: v^T p.
\end{align}
The reconstruction $\widetilde{v} \in \widetilde{V}_N \oplus \ker_V
L$ is given by 
\begin{align*}
  \widetilde{v} = \sum_{n=1}^N c_n f_n + \sum_{j=1}^3 w_j b_j
\end{align*}
where
\begin{subequations}
\label{eq:other-subspace-coefficients}
\begin{align}
  w &= \left(B^T G_{F}^{-1}  B \right)^{-1} B^T G_{F}^{-1} v \\
  c &= G_F^{-1} \lambda_F^T G_{F}^{-1} (v - B w).
\end{align}
\end{subequations}
\end{proposition}

Note that unlike the result of Corollary \ref{cor:v-field-wp-min}, the momenta
 $p \equiv L_{F}v$ are no longer the reconstruction coefficients for the lift
 unless $\widetilde{V}_N$ is the span of Green's functions centered at $\lmpt$. In any case, $L_F$ is a symmetric, rank $M-3$, positive semidefinite matrix.

Whether using Green's function interpolation or another basis, using
\eqref{eq:other-subspace-norm} and \eqref{eq:other-subspace-coefficients} to
compute either the norm or any coefficients can be problematic. The matrices in
the explicit formula given can easily be ill-conditioned, depending on $F$ and
$\mathcal{Q}$. We do not directly solve this issue here, but instead we reformulate the problem in
a format that we observe is more robust with respect to numerical precision issues than the above formulas. This is done in Section \ref{sec:computations}.

\subsection{Path energy and its gradient}\label{sec:pathlength-notation}

Computing the length and energy of a path is relatively straightforward given the work in the previous section. If $v^t$ is a vector of velocity evaluations at time $g$, we have
\begin{align*}
  E = \int_0^1 \|\widetilde{v}\|^2_{W P} \dx{t} \int_0^1 v^t L_F^t v^t \dx{t},
\end{align*}
where $L^t_F$ is the metric on the particle positions $q^t$.
%\subsection{The Gradient of Energy}
We compute the gradient of energy for our discrete boundary-value problem in two steps.  We first compute the gradient for the unconstrained problem; i.e. we compute an update to the velocity field which simply makes the energy smaller, ignoring the boundary conditions at times $0$ and $1$.  We then project this update onto the space of admissible updates; i.e. those updates which preserve the boundary conditions.  This two-step trick is equivalent to simply taking the first variation of energy for the fixed-endpoint problem, but is considerably simpler to derive, carries  little performance penalty, and makes it possible to alter the space of admissible updates at run-time easily.  

Since we have yet to discretize the temporal variable, we wait to compute the uncontrained gradient until Section \ref{sec:computations}. For now we concentrate on the projection into the space of admissible updates.

Define the energy of a velocity field as above.  Now consider two manifolds of
velocity fields.  $\mathcal{N}$ contains all velocity fields defined on $M$
particles at each time, and $\mathcal{M} \subset \mathcal{N}$ contains those
which respect the boundary-value conditions.
\begin{align*}
\mathcal{N} &= \left\{ \widetilde{v}:\{1,...,M\} \times [0,1] \to \R     \right\} \\
\mathcal{M} &= \left\{ V \in \mathcal{N} \left| q_m(1) = q_m(0) + \int_0^1 \widetilde{v}(\xi,q_m(\xi)) \dx{\xi} \right. \right\}
\end{align*}

%Let $P_{h_{0}}$ be the manifold of all paths $\phi: S^{1} \times [0,1] \to S^{1}$ which join $h_{0}$ to the identity; this is a submanifold of the manifold of all paths in $\Homeo(S^{1})$.   We construct this as an implicit submanifold given by the constraints that for each $\theta$ we must have
%\begin{align}
%\label{condition-submanifold-gets-there}
%\int_{0}^{1} v_{t}(\phi_{t}(\theta)) dt = h_{0}^{-1}(\theta) - \theta
%\end{align}
%along with topological constraints, say $\phi_{t}' > 0$ for all $t$.  
%\mfnote{I was thinking of the $\theta$'s as particles which initially lie in the range of $h_{0}$, which I think is consistent with how the code was working.}
%We now discretize time into $N$ steps and choose $M$ initial particles $\theta_{j}$ in the range of $h_{0}$.  Rather than storing their positions at each time step, we store their velocities $v_{j}(k)$ for $k \in 1,2,3,...,N$.  For each particle $\theta_{j}$, the condition \eqref{condition-submanifold-gets-there} on the velocity field becomes
%\begin{align*}
%F_{j}(v) = \sum_{k=1}^{N} v_{j}(k) = h_{0}^{-1}(\theta_{j}) - \theta_{j}
%\end{align*}
%if we use a simple first-order approximation $\phi_{n}(\theta_{j}) = \theta_{j}(0) + \sum_{k=1}^{n}v_{j}(k)$ to compute particle positions at each timestep.  Now let $F$ be the vector-valued function
%\begin{align}\label{eq:velocity-constraints}
%F(v) = \myvec{F_{1}(v) \\ F_{2}(v)\\ \vdots \\ F_{M}(v)}
%\end{align}
%\mfnote{Yes, it'll follow from basic ODE results applied to the geodesic equation.  I don't know the exact degree of regularity in $t$ here, but every result I've seen on this question just used the usual trick to convert the geodesic equation to a system of first-order ODEs and then cites the usual stuff.}
%
%
%
%
%In coordinate $v$, where $v(\cdot,t)$ is a vector field on $S^{1}$ for any fixed $t$,  the manifold $P_{h_{0}}$ has tangent space given by $\mathrm{null}[DF]$.  We can compute the path of least energy by finding the constrained gradient of energy $\widetilde \nabla E$; that is, the direction in $T_{\phi}P_{h_{0}}$ in which arc length increases most rapidly.  This is simply the projection onto $T_{\phi}P_{h_{0}}$ of the gradient of arc length in the whole space (i.e. the space of all paths regardless of endpoints).  This is trivial to compute; we will need an easy lemma first for computing the gradient of a function on a submanifold $M$ of a manifold $N$ equipped with an inner product.  Given a function $f:N \to \R$, let $Df$ be its differential.  The next lemma relates the gradient $D(\left.f\right|_{M})$, i.e. the differential of the restriction of $f$ to $M$, to the differential $Df$.

We defined the energy $E: \mathcal{N} \to \R$ as a function on $\mathcal{N}$, and in the next section we will describe the computation of its gradient $\nabla E$ on $\mathcal{N}$.  %With this gradient in hand, the next lemma allows us to find the direction of steepest ascent when restricting to admissible updates tangent to $M$.
With this gradient in hand, we obtain the direction of steepest ascent as a projection of $\nabla E$ (a coordinate-dependent quantity) onto the manifold $\mathcal{M}$ by projection in the following metric.

Let $\mathcal{P}$ be the space of smooth paths on $\UTS$.  Let $q(t,x,s) = q_{t;s}(x)$ be (a representative of) a smooth curve in $\mathcal{P}$.  Variable $x$ is position on $S^{1}$, $t$ parametrizes a path for given $s$, and $s$ parametrizes the curve of paths.  We define the velocity $v_{t;s}$ along the path by 
\begin{align*}
q(t,x,s) = q(0,x,s) + \int_{0}^{t} v_{\xi;s}(q(\xi,x,s)) d\xi
\end{align*}
I.e. for fixed $s$ the path of a particle $q(\cdot, x, s)$ is obtained by integrating the Eulerian vector field $v_{t,s}$ in $t$.  To construct a metric we place a norm on the Lie algebra and extend to the whole manifold by right translations.  We write
\begin{align}\label{eq:pathspace-velocity}
w_{t;s} = \ppx{s}\left( \dot{q} \circ q^{-1}\right)_{t;s}
\end{align}
and introduce the norm
\begin{align}\label{eq:pathspace-norm}
\left\| w \right\|_{\mathcal{P}}^{2} \equiv \int_{0}^{1} \| w_{t;s} \|_{WP}^{2} dt
\end{align}

%\aknote{
%  \input{content/pathvariation}
%}

\mfnote{
  We take the metric on path space to be the Riemannian metric induced by the norm $\| \cdot \|_{\mathcal{P}}$.  However, we wish to work on the manifolds $\mathcal{N}$ and $\mathcal{M}$.  With $q = q(t,x,s)$, $\dot{q} = \ppx{t} q$, and $q' = \ppx{x} q$, we have 
  \begin{align*}
    \pfpx{\left(\dot{q} \circ q^{-1}\right)_{t;r}}{r}(x) &= \left( \pfpx{\dot{q}}{s} \circ q^{-1}\right)_{t;r} + (\dot{q} \circ q^{-1})'_{t;r} \left(\pfpx{q^{-1}}{s}\right)_{t;r} \\
    %\pfpx{\left(\dot{q}_{r;t} \circ q_{r;t}^{-1}\right)}{s} + \left(\dot{q}_{t;r} \circ q_{t;r}^{-1}\right)' \ppx{s} q^{-1}_{t;r} \\
    &= \left(\pfpx{v}{s}\right)_{t;r} + v_{t;r}' \left(\pfpx{q^{-1}}{s}\right)_{t;r} \\
    &= \left(\pfpx{v}{s}\right)_{t;r} - v_{t;r}' \left(q^{-1}\right)'_{t;r} \left(\pfpx{q}{s}\circ q^{-1}\right)
  \end{align*}
  Direct computation yields
  \begin{align*}
  \frac{\partial q_{s}}{\partial s}(t,x) &= \frac{\partial}{\partial s} q(0,x,s) + \frac{\partial}{\partial s}\int_{0}^{t} v_{\xi;s}(q(\xi,x,s)) d\xi \\
  &= \int_{0}^{t} \frac{\partial}{\partial s} v_{\xi;s}(q(\xi,x,s)) d\xi 
  %&= \int_{0}^{t} \frac{\partial v^{\ell}_{\xi;s}}{\partial s}(x) d\xi
  \end{align*}
  Therefore
  \begin{align*}
    \pfpx{\left(\dot{q} \circ q^{-1}\right)_{t;r}}{r}(x) &= \mathcal{Y}\left[ \pfpx{v}{s}\right],
  \end{align*}
  where the linear operator $\mathcal{Y}$ is given by 
  \begin{align}\label{eq:operator-Y-definition}
    \mathcal{Y}(t,x,s) = \mathrm{Id} - v' \left(q^{-1}\right)' \int_0^t (\cdot)_{\xi} \dx{\xi}
  \end{align}
}
%where $v^{\ell}_{t;s}$ is the Langragian vector field 
%\begin{align*}
%v_{t;s}^{\ell}(x) \equiv v_{t;s}(q(t,x,s))
%\end{align*}
%This yields that $\ppx{t} \pfpx{q_s}{s} = \ppx{s} v^l_{t;s}(x)$.
This operator $\mathcal{Y}$ gives a mapping from velocity field updates (i.e. tangent vectors to $\mathcal{N}$ or $\mathcal{M}$) to tangent vectors $\frac{\partial q}{\partial s}$ to $\mathcal{P}$.  Hence our metric on path space immediately induces a metric on $\mathcal{N}$ and the submanifold $\mathcal{M}$.  We make our projection of the update to the velocity field with respect to this metric.

%\begin{lemma}\label{lemma:v-projection}
%At any point $V \in M$, let $w^{*} \in T_{V}M$ be the $L^{2}$ projection
%\begin{align*}
%w^{*} \equiv \mathrm{proj}_{\mathrm{null}[A]} \left( \nabla E \right)
%\end{align*}
%%maximizes $ \langle \nabla E, w  \rangle_{L^{2}}$ over all $w\in T_{\V{v}}M$.
%Then
%\begin{align*}
%w^{*} = \arg \max_{w \in T_{V}M} D_{w}E
%\end{align*}
%where $D_{w}$ is the directional derivative in direction $w$.
%\begin{proof}
%Split $\nabla E$ as
%\begin{align*}
%\nabla E = (\nabla E)^{T} + (\nabla E)^{\perp}
%\end{align*}
%where $(\nabla E)^{T}$ and $(\nabla E)^{\perp}$ are $L^{2}$-tangent and $L^{2}$-orthogonal to $M$, respectively.  Let $w \in T_{V}M$ and consider the variation $E(V + \epsilon w)$.  Then  
%\begin{align*}
%\frac{\partial E}{\partial \epsilon} &= \langle \nabla E, w \rangle_{L^{2}} \\
%&= \langle (\nabla E)^{T}, w \rangle_{L^{2}} + \langle (\nabla E)^{\perp}, w \rangle_{L^{2}} \\
%&= \langle (\nabla E)^{T}, w \rangle_{L^{2}}
%\end{align*}
%By the Schwarz inequality, the maximizing direction is $w^{*} =  (\nabla E)^{T} = \mathrm{proj}_{\mathrm{null}[A]}  \left( \nabla E \right)$.
%\end{proof}
%\end{lemma}


%At any point $p$ we may write 
%\begin{align*}
%\nabla E = (\nabla E)^{T} + (\nabla E)^{\perp}
%\end{align*}
%where $(\nabla E)^{T}$ and $(\nabla E)^{\perp}$ are $L^{2}$-tangent and $L^{2}$-orthogonal to $M$, respectively.  We seek a direction $v^{*}$, tangent to $M$, which maximizes
%\begin{align*}
%\langle \nabla E, v \rangle_{L^{2}} = \langle (\nabla E)^{T}, v  \rangle_{L^{2}} 
%\end{align*}
%By the Schwarz inequality, this is clearly given by $v^{*} =  (\nabla E)^{T} = \mathrm{proj}_{\mathrm{null}[A]}  \left( \nabla E \right)$.









\section{Computations}\label{sec:computations}
We now address the more practical issues of implementation of the methods described. We first formulate computation of the $W P$ norm as a classical numerical linear algebra problem. Following this is presentation of the temporal discretization; having discretized both the temporal and $S^1$ variables, we can compute the unconstrained coordinate gradient of the energy. Finally, the submanifold projection for gradient updates is discretized.

\subsection{The WP norm}
%In order to compute the horizontal lift $\widetilde{v}$, there is an alternative
%to direct computation of the matrix $L_{F} \equiv H_{F}^{-1}P_{H_{F}}^{B}$ as defined in Section
%\ref{sec:lifting-v}. 
It is not necessary to explicitly form the metric matrix $L_F$ as defined in Proposition \ref{prop:metric-definition}. We can accomplish the same task in a more stable manner using some established numerical linear algebra results. Given an $M
\times N$ matrix $A$, $A^\pinv$ denotes its Moore-Penrose pseudoinverse. For any
$\V{b} \in \R^N$, the vector $A^\pinv \V{b}$ is the minimum-norm, least squares
solution to $A \V{x} = \V{b}$.  Here, both minimum-norm and least
squares refer to the $\ell^2$ (Euclidean) metric on $\V{x}$ and $\V{b}$,
respectively. 

For computation of the $W P$ norm, we develop the ideas for a general subspace
$V_N$ of $V$ and then specialize to
the minimal-norm lift over all $V$ using Green's functions. Our goal is the interpolation of data
$\V{v} \in \R^M$ collocated at $\lmpt$ with an element from
$\widetilde{V}_N \oplus \ker_V L$; thsu we are essentially trying to find the minimal norm solution to
\begin{align}\label{eq:full-minnorm-system}
  \left[\begin{array}{cc} \lambda_F & B \end{array}\right] \left(
  \begin{array}{c} \V{c} \\ \V{w} \end{array} \right) =:
  A \V{x} = 
  \V{v},
\end{align}
where the norm over the coefficients $\V{c}$ is given by the positive (Gram)
matrix $G_F$ and that over the kernel coefficients $\V{w}$ is zero. Again we
assume that $\widetilde{V}_N$ is total for interpolation on $\lmpt$,
implying that $A \V{x} = \V{v}$ has at least one solution. Let $G_{F_0}$ be a
block diagonal matrix with $G_F$ in the upper-left block and a $3 \times 3$
zero matrix in the lower-right block.  It is easy to show that $\ker (B)$ is
trivial so long as $M \geq 3$.  This implies that $\ker(G_{F_0}) \cap \ker(A)$
is trivial and therefore the minimum norm solution to
\eqref{eq:full-minnorm-system} is unique \cite{eldn_weighted_1982}.  This is an
easy way to show uniqueness of the lift defined in
\eqref{eq:v-global-minimizer}.

Restating \eqref{eq:full-minnorm-system}, we are trying to solve the
minimization problem:
\begin{align}\label{eq:lse-problem}
  \text{minimize} \hspace{0.5cm} \left\| R \V{x} \right\|_2
  \hspace{0.5cm} \text{subject to} \hspace{0.5cm}
  A \V{x} = \V{v},
  %\left[\begin{array}{cc} \lambda(\Phi) & B \end{array}\right] \V{x} = \V{v},
\end{align}
where $R$ is any matrix square root of $G_{F_0}$. (I.e.
$R^T R = G_{F_0}$.) One of the standard tools for solving
least squares equality constrained problems is the generalized singular value
decomposition (GSVD) \cite{loan_generalizingsingular_1976}. Since $A$ and $G(F)_0$ have the same number of
columns, the GSVD matrix decomposition of $(A, R)$ is given by 
\begin{subequations}
\label{eq:gsvd}
\begin{align}
  A &= U C X^{-1} \\
  R &= V S X^{-1},
\end{align}
\end{subequations}
where $U$ and $V$ are orthogonal matrices, $X$ is invertible, and $C$ and $S$
each have only one non-vanishing diagonal with non-negative entries and satisfy
$C^T C + S^T S = I$. (Here, $V$ denotes a matrix in the GSVD decomposition, and not a linear subspace.) The solution to \eqref{eq:lse-problem} is given by
\begin{align}\label{eq:gsvd-solution}
  \V{x} = \left(\begin{array}{c} \V{c} \\ \V{w}\end{array}\right) &= X C^\pinv U^T \V{v},
\end{align}
The above matrix mapping $\V{v}$ to $\V{x}$ is a weighted pseudoinverse of $A$
with respect to $R$.  The momenta are $\V{p} = U C^{\pinv T} X^T G_{F_0}
\V{x} = L_F \V{v}$.  Our experience is that this is a much more stable way to compute
momenta than the explicit matrix relations used to define $L_F$. Note that
since $C^\pinv$ is a matrix with only one non-vanishing diagonal, applications of
$C^{\pinv}$ can be accomplished with a simple vector-vector multiply. The norm is
$\|\V{v}\|_{WP(\lmpt),\widetilde{V}_N} = \V{c}^T G_F\V{c} = \V{v}^T \V{p}$.

For the horizontal lift from the whole space $V$, we need only take the space
$\widetilde{V}_N$ to be the $M$-dimensional space of Green's functions centered
at the $q_m$. Notable simplifications in this case are that (a) $G_F = G$,
with $G$ being the Gram matrix of Green's functions \eqref{eq:G-matrix} and (b) the momentum and basis
coefficients coincide: $\V{p} = \V{c}$.

\subsection{Temporal discretization}
Since flow along the geodesic is reversible, we use a quadrature method that is symmetric with respect to the endpoints. We rewrite the flow of particles in an equivalent form
\begin{align}\label{eq:symmetric-quadrature}
  q_m(t) = \frac{1}{2} \left(q_m(0) + \int_0^t \widetilde{v}(\xi, q_m(\xi)) \dx{\xi}\right) + \frac{1}{2} \left(q_m(1) - \int_t^1 \widetilde{v}(\xi, q_m(\xi)) \dx{\xi}\right)
\end{align}
Now choose $T$ ordered points $s^t$ on $(0,1)$; let $\V{q}^t$ and $\V{v}^t$ denote the $M$ particle positions and velocities, respectively, at those times. The endpoint particle positions $\V{q}^0$ and $\V{q}^{T+1}$ are given.
Let $h^t$ denote the quadrature weight associated with data at time $s^t$ so that $\int_0^1 f(s) \dx{s} \simeq \sum_{t=1}^T h^t f(s^t)$.

Following \eqref{eq:symmetric-quadrature} we use the quadrature rule to
to integrate particle positions forward from $\V{q}^0$ to $\V{q}^t$, and backwards from $\V{q}^{T+1}$ to $\V{q}^t$, and average the result. For example if the representation is piecewise-linear, then the symmetric velocity-to-particle map is 
\begin{align}\label{eq:explicit-particle-update}
  \V{q}^t &= \frac{1}{2}\left(\V{q}^0 + \sum_{r=1}^{t-1} h^r \V{v}^r +
  \frac{h^t}{2} \V{v}^t\right) + \frac{1}{2}\left( \V{q}^{T+1} - \sum_{r=t+1}^T
  h^r \V{v}^r - \frac{h^t}{2} \V{v}^t\right), & t\geq 1
  %\V{q}^t &= \frac{h^1}{2} \V{v}^1 + \sum_{r=2}^t h^r \V{q}^r, & r>1
\end{align}
The scheme described above expresses particle positions $\V{q}^t$ linearly with the velocity field
values $\V{v}^r$. Collect all the velocity
evaluations into a matrix $V$, of size $M \times T$; do the same for the
particle positions in a matrix $Q$.  Then regardless of the choice of linear temporal quadrature rule, there exists a $T \times T$ matrix $Z$
such that 
\begin{align}\label{eq:particle-update}
  Q = \frac{1}{2} \left(Q^0 + Q^{T+1}\right) + V Z,
\end{align}
where $Q^0$ and $Q^{T+1}$ are matrices with the known vectors $\V{q}^0$ and
$\V{q}^{T+1}$ repeated. The entries of the matrix $Z$ depends on the choice of temporal representation and quadrature. For the piecewise-linear choice, it has entries $Z_{r,t} = \frac{1}{2} c_{r,t} h_r$, where 
\begin{align}\label{eq:crt-def}
c_{r,t} = \left\{\begin{array}{rl} 
                   0, & r = t, \\
                   1, & r < t, \\
                   -1, & r > t.
                 \end{array}
          \right.
\end{align}
In order to evaluate the energy $\int_0^1 \| v_s\|_{W P}^2 \dx{s}$, we build the quadrature factor $h^t$ into the norm at each point in time. This can be accomplished by simply replacing $G(F)$ by $h^t G(F)$ at each point in time. Therefore
\begin{align}\label{eq:pathlength-discretization}
  E = \sum_{t=1}^T \| \V{v}^t \|_{WP(\lmpt),\widetilde{V}_N}^2 = \sum_{t=1}^T 
  (\V{v}^t)^T L_{F} \V{v}^t %= \sum_{t=1}^T (\V{v}^t)^T \V{p}^t,
\end{align}
where $L_{F}$ is time-varying, depending explicitly on $\V{q}^t$,
the particle locations at time $s^t$, and also depending proportionally on $h^t$.

To keep the presentation simple we have described use of a piecewise-linear
quadrature. In practice we find it is sometimes more efficient to use a Legendre-Gauss-Lobatto quadrature rule for $t
\in [0,1]$, which enables a high-order polynomial representation for each
particle's velocity field. When to use a high-order representation versus a piecewise linear representation depends on whether one expects large $t$-derivatives in the velocity field.

\subsection{The gradient of the WP norm}
At some time $s^t$, we have seen that if $L_{F} = G_{F}^{-1}P_{G_{F}}^{B}$ then
%\begin{align*}
%\|\V{v}^t\|_{WP(\lmpt),\widetilde{V}_N} = (\V{v}^t)^{T} L_{F} \V{v}^t
%\end{align*}
%$H = Q^{-1} P_Q^B$, then
%$\|\V{v}^t\|_{WP(\lmpt),\widetilde{V}_N} = \V{v}^t H \V{v}^t$. 
formula \eqref{eq:pathlength-discretization} defines the norm, and to minimize
we must compute variations of this with respect to velocity evaluations. Since
the particle positions are influenced by \eqref{eq:particle-update}, the
matrices $L_F$ and therefore the energy change in nontrivial ways when we vary
$\V{v}^t$. To simplify the procedure, we first fix $t$ and consider only
variations in $\|\V{v}^t\|_{WP(\lmpt)}$ with respect to $\V{v}^r$ for $r =
1,\ldots,T$.

%The update \eqref{eq:explicit-particle-update} shows that
%$\V{q}^t$ does not depend on $\V{v}^t$. Therefore, $L_{F}$ at time $t$ is independent
%of $\V{v}^t$. As such, 
Direct variation of the quadratic form yields
\begin{align}\nonumber
  \ppx{\V{v}^r} \|\V{v}^t\|^2_{WP(\lmpt),\widetilde{V}_N} &= 2 \delta_{t,r} (\V{v}^t)^T L_{F} + (\V{v}^t)^T \left(\ppx{\V{v}^r} L_F \right) \V{v}^t \\ \label{eq:wpnorm-variation-velocity}
  &= 2 \delta_{t,r} (\V{p}^t)^T + (\V{v}^t)^T \left(\ppx{\V{v}^r} L_F \right) \V{v}^t
\end{align}
where $\delta_{i,j}$ is the Kronecker delta.
%and for any particle $m = 1, 2, \ldots, M$ we have 
%\begin{align*}
%  \ppx{q^t_m} \|\V{v}^t\|^2_{WP(\lmpt),\widetilde{V}_N} = (\V{v}^t)^T
%  \ppx{q^t_m}\left(L_{F}\right)
%  \V{v}^t,
%\end{align*}
The metric $L_F$ and its components $G_{F}$ and $P_{G_{F}}^B$ are generated from $\V{q}^t$, depending on $\V{v}^t$. Using simple properties of
matrix algebra, a straightforward but messy computation yields a formula for the
variation of the objective with respect to the particle positions.
\begin{lemma}\label{lemma:wpnorm-particle-variation}
Let $\lambda(F')$ be the $M \times N$ matrix with entries
$\lambda(F')_{m,n} = f_n'(q_m)$ and $B'$ be the $M \times 3$ matrix
with entries $B'_{m,j} = b_j'(q_m)$.  For two matrices $C$ and $D$ of the
same size, $C \circ D$ denotes the elementwise product.  Then 
\begin{align}\label{eq:wpnorm-variation-particle}
 \ppx{\V{q}^t} \|\V{v}^t\|^2_{WP(\lmpt),\widetilde{V}_N} = -2 \left[ \V{p}^t \circ
 \left(\lambda(F') G(F)^{-1} \lambda(F)^T \V{p}^t + B'
 \V{w}^t\right)\right]^T.
\end{align}
In particular, if $F$ is a basis of Green's functions centered at $\lmpt$, then 
\begin{align}\label{eq:wpnorm-variation-particle-hlift}
 \ppx{\V{q}^t} \|\V{v}^t\|^2_{WP(\lmpt)} = -2 \left[ \V{p}^t \circ
 \widetilde{v}'(\V{q}^t)\right]^T,
\end{align}
where $\widetilde{v}'$ is the derivative of the horizontal lift from $V$.
\end{lemma}
The salient result of Lemma \ref{lemma:wpnorm-particle-variation} is that
computing variations of the $W P$ norm with respect to particles is quite easy
once we have the momentum coefficients $\V{p}$ and kernel coefficients $\V{w}$
in hand from the GSVD. The only additionaly difficulty could come in computing
entries for $\lambda(F')$. In all the straightforward choices for $F$ we make,
it is not a deterrent. 

\subsection{The gradient of energy}
With the gradient of the $W P$ norm computed, computing the energy is now
straightforward. Combining \eqref{eq:particle-update},
\eqref{eq:wpnorm-variation-velocity}, and \eqref{eq:wpnorm-variation-particle},
we have
\begin{align}\label{eq:standard-gradient-formula}
%  \pfpx{E}{\V{v}^t} = 2 (\V{p}^t)^T - 
%                     \sum_{r \neq t} h_r c_{t,r} \pfpx{\|\V{v}^r\|^2_{WP(\lmpt)}}{\V{q}^r},
%  \pfpx{E}{\V{v}^t} = 2 (\V{p}^t)^T - 
%                     2 \sum_{r}  Z_{t,r} \pfpx{\|\V{v}^r\|^2_{WP(\lmpt)}}{\V{q}^r},
\frac{1}{2} \pfpx{E}{\V{v}^t} =  (\V{p}^t)^T - 
                      \sum_{r}  Z_{t,r} \left[ \V{p}^t \circ
 \left(\lambda(F') G(F)^{-1} \lambda(F)^T \V{p}^t + B'
 \V{w}^t\right)\right]^T
\end{align}
where again we recall that the
momenta coefficients $\V{p}^t$ already have the scaling factor $h^t$ included.

Although the above equation is an exact formula for the gradient, updating
velocity fields with these values will not respect the endpoint constraint
$\V{q}^1 = \V{q}^0 + \sum_t h_t \V{v}^t$, and so we project into the space of admissible updates.

The discrete versions of the manifolds presented in Section \ref{sec:pathlength-notation} can be defined as follows: 
\begin{align*}
N &= \left\{ \V{v}:[1,...,M] \times [1, ..., T] \to \R     \right\} \\
M &= \left\{ \V{v} \in N \left| \V{q}^1 = \V{q}^0 + \sum_t h_t \V{v}^t \right. \right\}
\end{align*}
Since the condition $\V{q}^1 = \V{q}^0 + \sum_t h_t \V{v}^t$ is linear, we may write it as
\begin{align*}
\V{q}^1 - \V{q}^0 = A \V{v}
\end{align*}
for some constraint matrix $A$. We see that $M$ is an implicit submanifold of
$N$: it is the $\V{q}^1 - \V{q}^0$ level set of the function $A \V{v}$.  By the
implicit function theorem, the tangent space to $M$ at any point is equal to
the nullspace of the differential $DA$, and we have simply $[DA] = A$. That is,
the admissible updates to a given velocity field $V$ are exactly those vectors
$\V{v}$ in the nullspace of $A$.

In order to perform the appropriate projection onto the nullspace of $A$, we place our update along a curve on the space of paths $\mathcal{P}$ and apply the norm induced by the metric \eqref{eq:pathspace-velocity} and \eqref{eq:pathspace-norm}. We then seek a discretization of the semidefinite operator $\mathcal{Y}^\dagger L \mathcal{Y}$ given by \eqref{eq:operator-Y-definition} and \eqref{eq:operator-L-definition}, where $\mathcal{Y}^\dagger$ is the $L^2$ adjoint of $\mathcal{Y}$.

At a fixed value of $s^t$, we have already approximated $L$ with $L_F$ depnding
on $\V{q}^t$.  Therefore we concentrate on $\mathcal{Y}$: the integral
$\int_0^{s^t} (\cdot)_\xi \dx{\xi}$ is approximated by the $t$-th column of
$Z$. Let $\widetilde{Z} = Z^T \otimes I_M$, where $I_M$ is the $M\times M$ identity
matrix; then $\V{q} = \widetilde{Z} \V{v}$. 
%We compute the vector
%$(q^{-1})'(\V{q})$ via {\em black magic}. 
The factor $v'(\V{q})$ can be
computed at each $t$ via $v'(\V{q}^t) = \lambda(F') \V{c}^t + B' \V{w}^t$. We
collect these factors into a diagonal $M T \times M T$ matrix $W$ with the
$t$-th block diagonal entry $(W^t)_{j,j} = v'(q^t_j) (q^{-1})'(q^t_j)$.  The
result then is the modified metric on updates $\ppx{\V{v}} E$:
\begin{align}\label{eq:discrete-path-metric}
  \widetilde{L_F} = (I - W \widetilde{Z})^T L_F (I - W \widetilde{Z}).
\end{align}
I.e., the matrix $Y := I - W \widetilde{Z}$ is our approximation to the operator $\mathcal{Y}$.
We use this metric both to form the natural gradient and to orthogonally project into the nullspace of the constraint matrix $A$.

%We can use this gradient as an update vector, but it is known that the natural gradient (cite Amari) is a more robust method to perform optimization. The natural gradient is formed by operating on the standard gradient by the inverse of the metric. In other words, for each column $(\widetilde{\nabla E})^t$ of $\widetilde{\nabla E}$, we form the natural gradient 
%\begin{align}\label{eq:natural-gradient}
%(\nabla E)^t_n = L_F^{+} \nabla E^t
%\end{align} 
%where $L_F$ is the metric on the particles $\V{q}^t$ at time $s^t$. Again, $L_F$ need not be formed explicitly, but instead the relevant matrices from the GSVD may be used to perform the mapping. Finally this natural gradient does not respect the endpoint conditions. In order to form the appropriate projection, we use the discrete versions of \eqref{eq:pathspace-velocity} and \eqref{eq:pathspace-norm}. The metric is the 
%
%, so we perform a final orthogonal projection, but this time it should be performed with respect to the metric $L_F$ rather than the Euclidean metric. We present a lemma that allows us an easier way to perform this final projection.
%\begin{lemma}
%  $(L_F - H_F^{-1}) L_F^+ = \V{0}$.
%\end{lemma}
%The proof is simple: $L_F^+ = P_{H_F}^B H_F$, and the result follows since
%$P_{H_F}^B$ is idempotent. Since $(\nabla E)^t_n \in \ran(L_F^+)$, we can perform the final natural gradient projection with respect to the strictly positive matrices $H_F$ rather than the semipositive $L_F$. Let $H^t_F$ denote the matrix at time $s^t$. Then combining the results of Lemmas \ref{lemma:v-projection} and \ref{lemma:v-minimizer}, our final gradient update matrix takes the form:
%\begin{subequations}\label{eq:natural-gradient-projection}
%\begin{align}
%  \V{u} &= H (\nabla E)_n \V{h}, \\
%  \widetilde{\nabla E}_n^t &= (\nabla E)^t_n - h^t H_F^t \V{u},
%\end{align}
%\end{subequations}
%where $H$ is the weighted Harmonic mean of the unprojected metric matrices $(H_F^t)^{-1}$:
%\begin{align*}
%  H = \left(\sum_{t=1}^T (h^t)^2 H_F^t \right)^{-1}.
%\end{align*}
%Again, if we take the horizontal lift so that $F$ is the span of collocated Green's functions, then considerable simplification occurs: $H_F^t = (h^t)^{-1} G^t$.

\section{Algorithm and Results}\label{sec:results}
With all the necessary derivations complete, a summary of the minimization algorithm is presented in Algorithm Listing \ref{alg:min-algorithm}.
One detail we have not yet addressed is a method for computing welding maps -- i.e. for determining $\phi$ given a discrete collection of ordered samples on $\partial \Omega$. A few standard methods exist for numerically computing conformal welds \cite{driscoll_schwarz-christoffel_2002,li_numerical_1998,symm_integral_1966,symm_numerical_1967}, and we settle upon the Zipper algorithm \cite{marshall_convergence_2007}; specifically we use the (simplest) `geodesic' version. We use the Zipper algorithm to (a) construct a representation of a conformal weld given discrete samples on the shape boundary $\partial \Omega$, (b) `interpolate' the welding map to any starting particle locations $\V{q}^0$,$\V{q}^{T+1}$ of our choosing, and (c) `invert' the weld: compute a representation of a simple closed curve from discrete samples of a homeomorphism of $S^1$. There are other possibilities for accomplishing these conformal welding tasks \cite{sharon_2d-shape_2006}.

\begin{algorithm}
\input{minimization}
\caption{A simple gradient descent algorithm for optimization.}\label{alg:min-algorithm}
\end{algorithm}

Our stopping criterion for the optimization iteration is defined by monitoring the relative objective decrease at each step in tandem with the Euclidean vector 2-norm of the projected natural gradient (normalized with respect to the size of the vector). At convergence, the former is less than $10^{-8}$ and the latter is less than $10^{-6}$.

One necessary measure of convergence is constancy of the $W P$ norm along the path. For all the results we show, the relative variation measure $(\max_t \|\V{v}^t\|^2_{W P(\lmpt)})/(\min_t \|\V{v}^t\|^2_{W P(\lmpt)}) - 1$ is no larger than $10^{-3}$ and is usually $\mathcal{O}(10^{-6})$.

In all the tests we use 150 particles. The full set of particles is obtained by computing 50 uniformly distributed particle samples on the ranges of the two welding maps $\phi_0$ and $\phi_1$, and 50 on the domain of the welding maps, and then taking the union of the sets. This encourages inclusion of resolvable features from both the interior and exterior of the welds in the computation.

For the temporal discretization, we have found that using a high-order Gauss-Lobatto quadrature rule works very well for welding maps whose derivatives are not large. When the welding maps exhibit high frequency content (a result of e.g. protrusions in the shape), then a standard piecewise linear choice works well. In this latter case we use 150 equispaced points in the time variable.

\subsection{Zero-distance welds}
As a first test we compute the energies of the paths connecting welding maps in
the same equivalence class (which are expected to vanish). The welding maps are
the maps $\phi_i$ in Figure \ref{fig:welding-maps}. In Table \ref{fig:equivalence-class-energies} we show the computed energy at convergence. The shown near-vanishing energies confirms validity of the computation.

\begin{figure}
  \begin{tabular*}{\textwidth}{m{0.35\textwidth}c}
    \begin{tabular}{r|c|}
        \vspace{-200pt}& \\  % gddamn tex hack
        \multicolumn{1}{r}{} & \multicolumn{1}{c}{$E$} \\\cline{2-2}
        $\phi_1 \leftrightarrow \phi_2$ & $7.530\times 10^{-12}$ \\\cline{2-2}
        $\phi_1 \leftrightarrow \phi_3$ & $1.249\times 10^{-9} $ \\\cline{2-2}
        $\phi_1 \leftrightarrow \phi_4$ & $2.853\times 10^{-9} $ \\\cline{2-2}
        $\phi_2 \leftrightarrow \phi_3$ & $1.729\times 10^{-8} $ \\\cline{2-2}
        $\phi_2 \leftrightarrow \phi_4$ & $6.633\times 10^{-8} $ \\\cline{2-2}
        $\phi_3 \leftrightarrow \phi_4$ & $2.069\times 10^{-10}$ \\\cline{2-2}
    \end{tabular}
    &
    \resizebox{0.6\textwidth}{!}{
      \includegraphics{zero-distance-welds}
    }
  \end{tabular*}
  \caption{Left: computed path energies between welding maps in the same equivalence
  class shown in Figure \ref{fig:welding-maps}. The near-zero value of the
  energies verifies that the computed paths are accurate geodesics. Right:
  evolution snapshots of $\phi(t,\theta)$ at $t=0, 0.2, 0.4, \ldots, 1$ where
  $\phi(0,\cdot) = \phi_3$ and $\phi(0,\cdot) = \phi_4$.}\label{fig:equivalence-class-energies}
\end{figure}

\subsection{Path length vs aspect ratio}
Let $\psi^{r}_t$ be a geodesic where $\psi^{r}_0$ is the welding map for an ellipse of aspect ratio $r$, and $\psi^{r}_1$ is the identity. In 
Figure \ref{fig:ellipse-aspect-ratio} we show results of simulations for path lengths $L[\psi^r]$ versus aspect ratio $r$ ranging from 1 to 5. The results suggest that the asymptotic relation is approximately linear, and these results match those obtained in \cite{kushnarev-2012} very well. 

\begin{figure}
  \begin{center}
    \resizebox{1.0\textwidth}{!}{
      \includegraphics{ellipse-aspect-ratio}
    }
  \end{center}
  \caption{Left: computed geodesic energy for evolving an ellipse of given aspect ratio to a circle using the minimization algorithm (solid line) and using the shooting procedure from \cite{kushnarev-2012} (dashed line). (The lines overlap.) Right: evolution of the ellipse shown at equidistant points of the path parameter $s$.}
  \label{fig:ellipse-aspect-ratio}
\end{figure}

\subsection{Hyperbolicity test}
We verify the negative curvature of Teichm\"uller space $\UTS$. Consider an ellipse of a fixed aspect ratio corresponding to the weld $\phi_1$. For rotations of the ellipse by angles $\pm 2 \pi/3$ we generate two more welding maps $\phi_2$ and $\phi_3$. We consider these three points in $\UTS$ and form a triangle using $\phi_i$ as the vertices and corresponding geodesic paths on $\UTS$ as the vertex connections. Let the $t=0$ velocity field that pushes $\phi_i$ to $\phi_j$ at $t=1$ be denoted $v_{i,j}$. Then the angle at vertex $\phi_i$ is given by
\begin{align*}
  \alpha_{i} = \frac{\langle v_{i,i\oplus 1}, v_{i, i\ominus 1} \rangle_{W P}}{\left\|v_{i,i\oplus 1}\right\|_{W P} \left\|v_{i,i\ominus 1}\right\|_{W P}},
\end{align*}
where $\oplus$ and $\ominus$ denote modular addition and subtraction, respectively, on the set $\{1, 2, 3\}$. Because all sectional curvatures of the $W P$ metric are negative, we expect that $\sum_{i=1}^3 \alpha_i \leq \pi$. We verify this fact in Figure \ref{fig:negative-curvature}. We again compare the results against those computed in \cite{kushnarev-2012} and obtain similar results.

\begin{figure}
  \begin{center}
    \resizebox{1.0\textwidth}{!}{
      \includegraphics{negative-curvature}
    }
  \end{center}
  \caption{Left: Angle sum versus aspect ratio using the minimization algorithm and the results from \cite{kushnarev-2012}. Right: evolution along geodesics for a triplet of rotated ellipses. Each ellipse-like shape is the true shape along the geodesic, but the translational placement along a hyperbolic triangle is a cartoon-like representation of the flow on $T(1)$, suggesting the negative curvature of the space.}
  \label{fig:negative-curvature}
\end{figure}

\subsection{Shapes with corners}
Non-smooth shapes do not result in diffeomorphic welds, and so the distance from any element of $T(1)$ to these non-smooth welds in the $W P$ metric is infinite. Our algorithm cannot compute a geodesic of infinite length, so we expect to run into limitations in the algorithm as we attempt to compute geodesics to shapes with sharp corners. We generate a triangular shape with rounded corners in the plane, and parameterize the smoothness of the corner by $\alpha \geq 1$. As shown on the right-hand side of Figure \ref{fig:sobolev-smoothness}, $\alpha=1$ corresponds to a sharp corner (a triangle in Euclidean space), and increasing $\alpha$ rounds the corners of the triangle. For many values of $\alpha$, we compute the geodesic from this rounded triangle to the identity and compile how the path length depends on $\alpha$. The results are shown on the left-hand side of Figure \ref{fig:sobolev-smoothness}. 

As we decrease $\alpha$ down to the critical value of $1$, we do observe a sharp growth in the path length. When $\alpha$ gets close to $1 + 10^{-3}$, our algorithm begins to reach its limit: at termination of the algorithm, the relative variation of the $W P$ norm along the path increases to about $0.1\%$. Results from computations with $\alpha -1 $ smaller than $10^{-3}$ no longer appear to be geodesics. (I.e., the $W P$ norm variation along the path becomes larger.) While this illustrates limitations of the algorithm, we are able to verify the expected increase in path length as $\alpha \downarrow 1$.

\begin{figure}
  \begin{center}
    \resizebox{1.0\textwidth}{!}{
      \includegraphics{sobolev-smoothness}
    }
  \end{center}
  \caption{Left: $W P$ distance from a circle (the identity weld) to a triangle with smoothness parameter $\alpha$. Right: triangular shapes corresponding to values of the $\alpha$ parameter.}
  \label{fig:sobolev-smoothness}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Shape images

\subsection{Shapes from the MPEG-7 dataset}
We conclude our investigation with computation of geodesics from shapes in the MPEG-7 CE-shape-1 collection of planar shapes\cite{mpeg7}. Figures \ref{fig:mpeg7-1333}, \ref{fig:mpeg7-785}, and \ref{fig:mpeg7-875} show computed geodesics, and list the path length in each case. Also, contour graphs of the velocity field $v(t,\theta)$ are displayed as well. 

We note that many shapes in this database have welding maps whose derivatives are singular to machine precision, or vanish to machine precision.\footnote{In this case the conformal maps comprising the weld are `crowded'; the max-min ratio of the derivative magnitude is large. A robust strategy for finite-precision computation in this situation is an open problem.}

When this happens, any algorithm implemented in finite precision for computing geodesics from this welding map will fail. We anticipate that any future work aimed at fixing the problem will first need to address the difficulty in computing conformal welds for these problematic shapes.

\begin{figure}
  \begin{center}
    \resizebox{1.0\textwidth}{!}{
      \includegraphics{mpeg7-1333}
    }
  \end{center}
  \caption{Left: shape evolution for {\sc mpeg-7} shape 1333, a rotated tree. Right: contour graph of the velocity field. The path length on $T(1)$ is $10.25$.} 
  \label{fig:mpeg7-1333}
\end{figure}

%\tikzsetnextfilename{mpeg7-682}
%\begin{figure}
%  \begin{center}
%    \resizebox{1.0\textwidth}{!}{
%      \input{tikz/mpeg7-682}
%    }
%  \end{center}
%  \caption{Shape evolution and velocity field for {\sc mpeg7} shape 682. PL=6.802 (great)}
%  \label{fig:mpeg7-682}
%\end{figure}

\begin{figure}
  \begin{center}
    \resizebox{1.0\textwidth}{!}{
      \includegraphics{mpeg7-785}
    }
  \end{center}
  \caption{Left: shape evolution for {\sc mpeg-7} shape 785, a fish. Right: contour graph of the velocity field. The path length on $T(1)$ is $8.889$.}
  \label{fig:mpeg7-785}
\end{figure}

%\tikzsetnextfilename{mpeg7-811}
%\begin{figure}
%  \begin{center}
%    \resizebox{1.0\textwidth}{!}{
%      \input{tikz/mpeg7-811}
%    }
%  \end{center}
%  \caption{Shape evolution and velocity field for {\sc mpeg7} shape 811. PL=4.457 (pretty good)}
%  \label{fig:mpeg7-811}
%\end{figure}

\begin{figure}
  \begin{center}
    \resizebox{1.0\textwidth}{!}{
      \includegraphics{mpeg7-875}
    }
  \end{center}
  \caption{Left: shape evolution for {\sc mpeg-7} shape 875, a rotated mitre. Right: contour graph of the velocity field. The path length on $T(1)$ is $7.486$.}
  \label{fig:mpeg7-875}
\end{figure}

%\tikzsetnextfilename{mpeg7-1374}
%\begin{figure}
%  \begin{center}
%    \resizebox{1.0\textwidth}{!}{
%      \input{tikz/mpeg7-1374}
%    }
%  \end{center}
%  \caption{Shape evolution and velocity field for {\sc mpeg7} shape 1374. PL=10.95 (ok)}
%  \label{fig:mpeg7-1374}
%\end{figure}

%\tikzsetnextfilename{mpeg7-1377}
%\begin{figure}
%  \begin{center}
%    \resizebox{1.0\textwidth}{!}{
%      \input{tikz/mpeg7-1377}
%    }
%  \end{center}
%  \caption{Shape evolution and velocity field for {\sc mpeg7} shape 1377. PL=9.887 (ok)}
%  \label{fig:mpeg7-1377}
%\end{figure}

%\tikzsetnextfilename{mpeg7-631}
%\begin{figure}
%  \begin{center}
%    \resizebox{1.0\textwidth}{!}{
%      \input{tikz/mpeg7-631}
%    }
%  \end{center}
%  \caption{Shape evolution and velocity field for {\sc mpeg7} shape 631., PL=12.34 (ok)}
%  \label{fig:mpeg7-631}
%\end{figure}

%\tikzsetnextfilename{mpeg7-903}
%\begin{figure}
%  \begin{center}
%    \resizebox{1.0\textwidth}{!}{
%      \input{tikz/mpeg7-903}
%    }
%  \end{center}
%  \caption{Shape evolution and velocity field for {\sc mpeg7} shape 903. PL = 4.889 (bad)}
%  \label{fig:mpeg7-903}
%\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Shape images

\section{Conclusion}
In this paper we have proposed and implemented a numerical solver for computing geodesics on the universal Teichm\"uller space $T(1)$ in the $W P$ metric. Given a planar shape, a conformal weld in the coset space $\UTS$ uniquely identifies this shape. We prescribe two given welding maps as the endpoints of the path in $\UTS$, and minimize the $W P$ energy of the path, subject to the endpoint constraints. This provides a constructive way to metrize the space of planar shapes, modulo rigid translation and scaling. 

We demonstrated the applicability of our method by solving for numerous geodesics; our results compare favorably with the results obtained by a shooting method in \cite{kushnarev-2012}. The main difficulty with our procedure stems from the inability to handle welds whose composite conformal maps contain large derivatives; this is in fact a limitation for any procedure of which we are aware. Our future work will be directed at addressing this problem.

Our optimization algorithm can, in effect, be used to minimize energy in any Riemannian metric, and our approach for satisfying endpoint constraints on manifolds can likewise be utilized in other applications.

\newpage
\bibliographystyle{amsplain}
\bibliography{wp-bib}

\end{document}  

wp-geodesics-matt2.bbl
\providecommand{\bysame}{\leavevmode\hbox to3em{\hrulefill}\thinspace}
\providecommand{\MR}{\relax\ifhmode\unskip\space\fi MR }
% \MRhref is called by the amsart/book/proc definition of \MR.
\providecommand{\MRhref}[2]{%
  \href{http://www.ams.org/mathscinet-getitem?mr=#1}{#2}
}
\providecommand{\href}[2]{#2}
\begin{thebibliography}{10}

\bibitem{mpeg7}
\emph{{Shape data for the MPEG-7 core experiment CE-Shape-1}},
  \url{http://www.cis.temple.edu/~latecki/TestData/mpeg7shapeB.tar.gz}.

\bibitem{ahlfors_lectures_2006}
Lars~V. Ahlfors, \emph{Lectures on quasiconformal mappings}, 2 ed., American
  Mathematical Society, August 2006.

\bibitem{bowick_string_1987}
M.~J. Bowick and S.~G. Rajeev, \emph{String theory as the {K\"ahler} geometry
  of loop space}, Physical Review Letters \textbf{58} (1987), no.~11, 1158.

\bibitem{driscoll_schwarz-christoffel_2002}
Tobin~A. Driscoll and Lloyd~N. Trefethen, \emph{{Schwarz-Christoffel} mapping},
  1 ed., Cambridge University Press, June 2002.

\bibitem{eldn_weighted_1982}
Lars Eldén, \emph{A weighted pseudoinverse, generalized singular values, and
  constrained least squares problems}, {BIT} \textbf{22} (1982), no.~4,
  487--502.

\bibitem{gay-balmaz_geometry_2009}
Francois Gay-Balmaz, Jerrold~E. Marsden, and Tudor~S. Ratiu, \emph{The geometry
  of the universal {Teichm\"uller} space and the {Euler-Weil-Petersson}
  equations}, Tech. report, 2009.

\bibitem{hubbard_teichmuller_2006}
John~H Hubbard, \emph{{Teichm\"uller} theory and applications to geometry,
  topology, and dynamics}, Matrix Editions, Ithaca, {NY}, 2006 (English).

\bibitem{kushnarev-2012}
S.~Kushnarev and A.~Narayan, \emph{Approximating the {Weil-Petersson} metric
  geodesics on the universal {Teichm\"uller} space by singular solutions},
  Submitted (2012).

\bibitem{kushnarev_teichons_2009}
Sergey Kushnarev, \emph{Teichons: Solitonlike geodesics on universal
  teichmüller space}, Experimental Mathematics \textbf{18} (2009), no.~3,
  325--336.

\bibitem{li_numerical_1998}
Bao~Cheng Li and Stavros Syngellakis, \emph{Numerical conformal mapping based
  on the generalised conjugation operator}, Mathematics of Computation
  \textbf{67} (1998), no.~222, 619--639, {ArticleType:} research-article / Full
  publication date: Apr., 1998 / Copyright © 1998 American Mathematical
  Society.

\bibitem{loan_generalizingsingular_1976}
Charles F.~Van Loan, \emph{Generalizing the singular value decomposition},
  {SIAM} Journal on Numerical Analysis \textbf{13} (1976), no.~1, 76--83.

\bibitem{marshall_convergence_2007}
Donald~E. Marshall and Steffen Rohde, \emph{Convergence of a variant of the
  zipper algorithm for conformal mapping}, {SIAM} Journal on Numerical Analysis
  \textbf{45} (2007), no.~6, 2577.

\bibitem{mumford_pattern_2002}
David Mumford, \emph{Pattern theory: the mathematics of perception}, in {ICM},
  Higher Education Press, 2002, p.~2002.

\bibitem{nag_rm_1990}
Subhashis Nag and Alberto Verjovsky, \emph{{$\textrm{Diff}(S^1)$} and the
  {Teichm\"uller} spaces}, Communications in Mathematical Physics (1965-1997)
  \textbf{130} (1990), no.~1, 123--138.

\bibitem{segal_definition_2004}
Graeme Segal, \emph{The definition of conformal field theory}, Topology,
  geometry and quantum field theory, London Math. Soc. Lecture Note Ser., vol.
  308, Cambridge Univ. Press, Cambridge, 2004, p.~421–577.

\bibitem{sharon_2d-shape_2006}
E.~Sharon and D.~Mumford, \emph{{2D-Shape} analysis using conformal mapping},
  International Journal of Computer Vision \textbf{70} (2006), no.~1, 55--75.

\bibitem{symm_integral_1966}
George~T. Symm, \emph{An integral equation method in conformal mapping},
  Numerische Mathematik \textbf{9} (1966), no.~3, 250--258.

\bibitem{symm_numerical_1967}
\bysame, \emph{Numerical mapping of exterior domains}, Numerische Mathematik
  \textbf{10} (1967), no.~5, 437--445.

\bibitem{teo_weil-petersson_2006}
Leon A.~Takhtajan Teo and Lee-Peng, \emph{{Weil-Petersson} metric on the
  universal {Teichm\"uller} space}, American Mathematical Society, August 2006.

\bibitem{veltkamp_shape_2001}
Remco~C. Veltkamp, \emph{Shape matching: Similarity measures and algorithms},
  Shape Modeling and Applications, International Conference on (Los Alamitos,
  {CA}, {USA}), {IEEE} Computer Society, 2001, p.~0188.

\end{thebibliography}
