# Dynamical Programming for off-the-grid dynamic Inverse Problems

Vincent Duval  Robert Tovey

(CEREMADE, CNRS, UMR 7534, Université Paris-Dauphine, PSL University 75016
Paris, France
INRIA-Paris, MOKAPLAN, 75012 Paris, France
)

###### Abstract

In this work we consider algorithms for reconstructing time-varying data into
a finite sum of discrete trajectories, alternatively, an off-the-grid sparse-
spikes decomposition which is continuous in time. Recent work showed that this
decomposition was possible by minimising a convex variational model which
combined a quadratic data fidelity with dynamical Optimal Transport. We
generalise this framework and propose new numerical methods which leverage
efficient classical algorithms for computing shortest paths on directed
acyclic graphs. Our theoretical analysis confirms that these methods converge
to globally optimal reconstructions which represent a finite number of
discrete trajectories. Numerically, we show new examples for unbalanced
Optimal Transport penalties, and for balanced examples we are 100 times faster
in comparison to the previously known method.

##  1 Introduction

The signal-processing task of dynamical super-resolution involves retrieving
fine-scale features, in space and time, from a signal which evolves over time.
A convex variational model was recently proposed for such tasks using Optimal
Transport (OT) to regularise the associated inverse problem [Bredies et al.,
2021, Bredies and Fanzon, 2020]. This new approach allows the decomposition of
a signal into a finite sum of smooth curves, for example to track the centers
of multiple particles in time with smooth trajectories. Similar ideas were
explored for a specific example in Alberti et al. [2019] where the shape of
curves is built into the model, and without appealing to OT.

In this work we focus on dynamical super-resolution problems regularised by
OT. We can consider potential models to be partitioned into two classes,
depending on whether the particles have constant mass/brightness in time, or
if mass is allowed to vary. These classes are referred to as balanced or
unbalanced problems respectively, mathematically encoded in the choice of OT
cost. Current literature provides analysis for the balanced Benamou-Brenier
(BB) [Bredies et al., 2021] and the unbalanced Wasserstein-Fisher-Rao (WFR)
[Bredies et al., 2020a] energies, both are shown to reconstruct data into a
finite number of smooth curves with constant or smoothly varying mass. Initial
numerical experiments for the Benamou-Brenier model have also been carried out
showing great promise [Bredies et al., 2020b], however current methods are too
slow for large-scale applications.

Let $\Omega$ be an open bounded spatial domain. At the heart of the analysis
of Bredies et. al. is the interplay between measures $\rho(t,x)$ defined on
the time-space cylinder $\left[0,1\right]\times\overline{\Omega}$, and
measures $\sigma(\gamma)$ defined on the space of curves $\gamma=(h,\xi)$ with
mass $h\in\mathscr{C}(\left[0,1\right])$ and trajectory
$\xi\in\mathscr{C}(\left[0,1\right];\overline{\Omega})$. We can think of
$\rho$ as representing the evolving physical volume that can be observed with
a microscope, whereas $\sigma$ more efficiently represents a collection of
(trajectories of) particles, that we wish to reconstruct. In the present work,
building upon the above-mentioned interplay, we focus on variational problems
in the space of measures $\sigma$ on paths. Such a point of view is fruitful,
as the problem more closely resembles the Beurling-LASSO which is now well-
understood [Bredies and Pikkarainen, 2013, Azaïs et al., 2015, Duval and
Peyré, 2015, Poon et al., 2021]. But its main advantage is that it paves the
way for “off-the-grid” numerical methods when solving such dynamical inverse
problems. Recent works in the field of sparse spike recovery [Bredies and
Pikkarainen, 2013, Boyd et al., 2017, Denoyelle et al., 2019] have
demonstrated that it is possible to design efficient numerical solvers without
reconstructing the unknown on a grid, by exploiting a conditional gradient
descent / Frank-Wolfe approach together with a good knowledge of the
regularising term. Indeed, the Frank-Wolfe minimisation algorithm and its
variants (see the review Jaggi [2013]) build iterates that are convex
combinations of the extreme points of level sets of the regulariser; being
able to easily encode and handle such extreme points makes it possible to
solve variational problems in a continuous (or up to floating point) setting.
Moreover, having iterates that are convex combinations of a few extreme points
of the level sets of the regulariser is particularly relevant, as it is known
in inverse problems that some solutions have precisely that structure when the
observation consists of a finite number of linear measurements [Unser et al.,
2017, Bredies and Carioni, 2019, Boyer et al., 2019].

In the particular case of dynamical super-resolution using optimal transport
regularisation, those extreme points are described by the key results of
Bredies et al. [2021, 2020a] (see also Ambrosio et al. [2008]) as measures
supported on the above-mentioned curves
$\gamma=(h,\xi)\colon\left[0,1\right]\to\mathds{R}\times\overline{\Omega}$.
Even though only a small subset of
$L^{\infty}(\left[0,1\right];\mathds{R}\times\overline{\Omega})$ is involved,
optimising over such a set is a challenging task that is addressed in Bredies
et al. [2020b] in the case of balanced optimal transport. Yet, the computation
times are hardly compatible with a practical applications.

The aim of the present paper is to propose an alternative optimisation method
which keeps the key properties of Bredies et al. [2020b] (i.e. iterates made
of a combination of atoms in a continuous setting) while drastically improving
the computation time.

###  1.1 Outline of the paper

In Section 2 we gather several results on the probabilistic representation of
the solutions of the continuity equation that are useful to understand the
perspective of measures on paths. Section 3 studies the case of the Benamou-
Brenier regularisation: we explain how the original problem is equivalent to a
variational problem on the space of paths, and we describe a new algorithm to
solve it, which is a particular instance of the framework introduced in
Section 5. Our main numerical contribution is to reformulate the main Frank-
Wolfe step as a shortest-path problem on an ordered, weighted, directed,
acyclic graph. This is demonstrated in detail in Section 3.4. Section 4
describes the unbalanced case, corresponding to the Wasserstein-Fisher-Rao
regularisation: for clarity, we have decided to present both cases separately.

The main contributions of this work are presented in Section 5 where we
introduce our general variational framework, numerical algorithm, and
convergence theorem. The proofs of convergence for inexact and stochastic
Frank-Wolfe are given in Section 6 where we confirm that our proposed
algorithm converges almost-surely to an optimal reconstruction of our
variational model. The stochastic algorithm is similar to that considered by
Silveti-Falls et al. [2020], although we compute the step-size adaptively
rather than fixing it a priori. Our novel reformulation to the shortest-path
problem is explained in Section 7, the dynamical programming approach is
standard but presented in the current context for clarity.

Finally, in Section 8 we present thorough numerical results demonstrating the
properties of stochastic and deterministic variants of our proposed algorithm.
Examples with the BB energy are compared with the algorithm of Bredies et al.
[2020b], and new unbalanced examples are reconstructed with a WFR penalty.

###  1.2 Notation

##### Convex sets and extreme points.

Let $V$ be a linear space. For all $\sigma_{0},\sigma_{1}\in V$, we define the
closed line segment between $\sigma_{0}$ and $\sigma_{1}$ as
$\left[\sigma_{0},\sigma_{1}\right]\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left\\{\,\lambda\sigma_{0}+(1-\lambda)\sigma_{1}\mid
0\leq\lambda\leq 1\,\right\\}$. Similarly, we define the open line segment
$\left]\sigma_{0},\sigma_{1}\right[\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left[\sigma_{0},\sigma_{1}\right]\setminus\\{\sigma_{0},\sigma_{1}\\}$.

######  Definition 1.1.

Let $D\subseteq V$. $D$ is called _convex_ if
$\left]\sigma_{0},\sigma_{1}\right[\subset D$ for all
$\sigma_{0},\sigma_{1}\in D$. We say
$\sigma\in\operatorname{Ext}\left(D\right)$ is an _extreme point_ of $D$ if
$\sigma\in D$ and there are no points $\sigma_{0},\sigma_{1}\in D$ such that
$\sigma\in\left]\sigma_{0},\sigma_{1}\right[$. In other words,

$\forall\lambda\in\left]0,1\right[,\forall\sigma_{0},\sigma_{1}\in
D,\quad\left(\sigma=\lambda\sigma_{0}+(1-\lambda)\sigma_{1}\right)\quad\implies\quad\left(\sigma_{0}=\sigma_{1}=\sigma\right).$
(1.1)

Furthermore, it is possible to define the notion of face (and elementary face)
of a convex set, which extends the notion of extreme point to higher-
dimensional sets. We refere to Dubins [1962] for more detail.

##### Measure spaces.

For a separable metric space $\Gamma$ and Banach space $X$, we define
$\mathscr{C}_{b}(\Gamma;X)$ to be the set of continuous bounded functions from
$\Gamma$ to $X$. When $X=\mathds{R}$, we simply write
$\mathscr{C}_{b}(\Gamma)$. Recall that for any Borel measure $\sigma$ on
$\Gamma$, we can define the unique non-negative Borel measure
${\left|\sigma\right|}\in\mathcal{M}^{+}(\Gamma)$ by

${\left|\sigma\right|}(A)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\sup\left\\{\,\sum_{i=1}^{n}{\left|\sigma(A_{i})\right|}\mid
n\in\mathds{N}^{*},\\{A_{1},\ldots,A_{n}\\}\ \mbox{Borel partition of
$A$}\,\right\\}$ (1.2)

for all Borel measurable sets $A\subset\Gamma$. We denote by
$\mathcal{M}(\Gamma)$ the space of signed Borel measures $\sigma$ with finite
total variation, i.e.
${\left\lVert\sigma\right\rVert}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\int_{\Gamma}\mathop{}\\!\mathrm{d}{\left|\sigma\right|}<+\infty$.
The total variation ${\left\lVert\cdot\right\rVert}$ defines a norm on
$\mathcal{M}(\Gamma)$, but it is sometimes more convenient to use the _narrow
topology_ , i.e. the weakest topology on $\mathcal{M}(\Gamma)$ which makes the
integration against continuous bounded functions a continuous linear form. The
narrow topology is equivalent to the weak-* topology on
$(\mathscr{C}_{b}(\Gamma))^{\prime}$, in particular, a sequence
$\\{\sigma^{n}\\}_{n\in\mathds{N}}\subset\mathcal{M}(\Gamma)$ converges to
$\sigma^{*}\in\mathcal{M}(\Gamma)$ in the narrow topology (denoted
$\sigma^{n}\stackrel{{\scriptstyle*}}{{\rightharpoonup}}\sigma$) if

$\forall\phi\in\mathscr{C}_{b}(\Gamma),\quad\lim_{n\to+\infty}\int_{\Gamma}\phi\mathop{}\\!\mathrm{d}\sigma^{n}=\int_{\Gamma}\phi\mathop{}\\!\mathrm{d}\sigma^{*}.$
(1.3)

The support of $\sigma\in\mathcal{M}^{+}(\Gamma)$ is defined as

$\operatorname{supp}(\sigma)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left(\bigcup\left\\{\,U\mid\sigma(U)=0,\ U\text{ is
open}\,\right\\}\right)^{c}.$ (1.4)

This is a closed set satisfying
$\sigma(\operatorname{supp}(\sigma))=\sigma(\Gamma)$.

##### Function domains.

In this work we consider two measure domains, the time-space cylinder

${\left[0,1\right]\times\overline{\Omega}}\quad\text{for an open, bounded,
convex domain }\quad\Omega\subseteq\mathds{R}^{d},\ d\geq 1$ (1.5)

and a closed set $\Gamma$ of continuous (weighted) curves which are viewed as
pairs $\gamma=(h,\xi)$ where $h(t)$ is the mass at time $t$ and $\xi(t)$ is
the location. A formal definition will be given in Lemma 2.1. In the time-
space cylinder we will denote measures
$\rho\in\mathcal{M}({\left[0,1\right]\times\overline{\Omega}})$ with test
functions $\psi\in\mathscr{C}({\left[0,1\right]\times\overline{\Omega}})$, and
similarly on the space of curves $\sigma\in\mathcal{M}(\Gamma)$ and
$\phi\in\mathscr{C}_{b}(\Gamma)$.

##### Narrowly continuous measures.

An important subspace of
$\mathcal{M}({\left[0,1\right]\times\overline{\Omega}})$ is the space of
narrowly continuous measures. With a slight abuse of the standard notation, we
will say
$\rho\in\mathscr{C}_{w}(\left[0,1\right];\mathcal{M}(\overline{\Omega}))\subset\mathcal{M}({\left[0,1\right]\times\overline{\Omega}})$
if there exists a map $t\mapsto\rho_{t}\in\mathcal{M}(\overline{\Omega})$
(informally, $\rho_{t}=\rho(t,\cdot)$ is a “time slice” of $\rho$) such that

$\forall\psi\in\mathscr{C}(\overline{\Omega}),\qquad\left[t\mapsto\int_{\overline{\Omega}}\psi(x)\mathop{}\\!\mathrm{d}\rho_{t}(x)\right]\in\mathscr{C}(\left[0,1\right])$
(1.6)

and

$\forall\psi\in
L^{1}_{\rho}({\left[0,1\right]\times\overline{\Omega}}),\qquad\int_{\left[0,1\right]\times\overline{\Omega}}\psi(t,x)\mathop{}\\!\mathrm{d}\rho(t,x)=\int_{0}^{1}\left(\int_{\overline{\Omega}}\psi(t,x)\mathop{}\\!\mathrm{d}\rho_{t}(x)\right)\mathop{}\\!\mathrm{d}t.$
(1.7)

Given a measure on paths $\sigma\in\mathcal{M}(\Gamma)$ such that
$\int_{\Gamma}{\left\lVert
h\right\rVert}_{\infty}\mathop{}\\!\mathrm{d}\sigma(h,\xi)<+\infty$, one may
define the family of measures
$(e_{t})_{\sharp}{\sigma}\in\mathcal{M}(\overline{\Omega})$, for
$t\in\left[0,1\right]$, by

$\forall\psi\in\mathscr{C}(\overline{\Omega}),\qquad\int_{\overline{\Omega}}\psi(x)\mathop{}\\!\mathrm{d}{[(e_{t})_{\sharp}{\sigma}]}(x)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\int_{\Gamma}h(t)\psi(\xi(t))\mathop{}\\!\mathrm{d}\sigma(h,\xi).$
(1.8)

Formally, $(e_{t})_{\sharp}{\sigma}$ is the image measure of $\sigma$ by the
evaluation at time $t$. That family is narrowly continuous, and as we explain
in Theorem 2.2 below, it is the evaluation (disintegration) of some measure
$\Theta(\sigma)\in\mathscr{C}_{w}(\left[0,1\right];\mathcal{M}(\overline{\Omega}))$.

##### The continuity equation.

######  Definition 1.2.

Let $\rho\in\mathcal{M}({\left[0,1\right]\times\overline{\Omega}})$ be a
measure. We say that $\rho$ satisfies the _continuity equation_ if there
exists $v\in
L^{1}_{{\left|\rho\right|}}({\left[0,1\right]\times\overline{\Omega}};\mathds{R}^{d})$,
$g\in L^{1}_{{\left|\rho\right|}}({\left[0,1\right]\times\overline{\Omega}})$
such that

$\forall\psi\in\mathscr{C}_{c}^{1}(\left]0,1\right[\times\overline{\Omega}),\qquad\int[\partial_{t}\psi+\nabla\psi\cdot
v+\psi g]\mathop{}\\!\mathrm{d}\rho=0.$ (1.9)

##### Lower semi-continuity.

A function $w\colon\Gamma\to\mathds{R}\cup\\{+\infty\\}$ is called _lower
semi-continuous_ if the sets $\left\\{\,\gamma\in\Gamma\mid
w(\gamma)\leq\lambda\,\right\\}$ are closed for each $\lambda\in\mathds{R}$.
As $\Gamma$ is a metric space, it is equivalent to say that for all
$\gamma\in\Gamma$ and sequences $\gamma_{n}\to\gamma$,
$w(\gamma)\leq\liminf_{n\to+\infty}w(\gamma_{n})$.

##  2 Preliminary results

In this work we define ${\Gamma_{0}}$ to be the set of continuous weighted
paths in $\overline{\Omega}$, modelling particles with (varying) mass $h(t)$
at location $\xi(t)$, defined by

${\Gamma_{0}}=\left\\{\,\gamma=(h,\xi)\mid h\in\mathscr{C}(\left[0,1\right]),\
\xi\colon\left[0,1\right]\rightarrow\overline{\Omega},\ \xi|_{\\{h\neq
0\\}}\text{ is continuous}\,\right\\}.$ (2.1)

For technical reasons we permit curves $\xi$ which may not be continuous at
points $t$ where $h(t)=0$. Heuristically, the location $\xi(t)$ is not
necessarily important if the particle has no mass and cannot be observed. The
other technical requirements for the set ${\Gamma_{0}}$ is that it must be a
complete separable metric space, as explored in the following lemma. We follow
the suggestion of [Bredies et al., 2020a, Prop. 3.6] who used the flat metric
on the space
$\left\\{\,h\delta_{\xi}\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})\mid(h,\xi)\in\Gamma,\
h\geq 0\,\right\\}$.

######  Lemma 2.1.

Define
$\operatorname{d}_{\Gamma}\colon{\Gamma_{0}}\times{\Gamma_{0}}\to\left[0,+\infty\right[$
by

$\displaystyle\operatorname{d}_{\Gamma}((h_{1},\xi_{1}),(h_{2},\xi_{2}))$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\sup_{t\in\left[0,1\right]}\operatorname{d}_{F}((h_{1}(t),\xi_{1}(t)),(h_{2}(t),\xi_{2}(t)))\qquad\text{where}$
(2.2) $\displaystyle\operatorname{d}_{F}((r_{1},x_{1}),(r_{2},x_{2}))$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\begin{cases}|r_{1}|+|r_{2}|&r_{1}r_{2}\leq 0\text{ or
}|x_{1}-x_{2}|\geq 2\\\
|r_{1}-r_{2}|+\min(|r_{1}|,|r_{2}|)|x_{1}-x_{2}|&\qquad\text{else,}\end{cases}$
(2.3)

then $(\nicefrac{{{\Gamma_{0}}}}{{\sim}},\operatorname{d}_{\Gamma})$ is a
complete separable metric space where

$(h_{1},\xi_{1})\sim(h_{2},\xi_{2})\iff h_{1}=h_{2}\quad\text{and}\quad\forall
t\in\\{h_{1}\neq 0\\},\quad\xi_{1}(t)=\xi_{2}(t).$ (2.4)

Convergence of a sequence $\gamma_{n}=(h_{n},\xi_{n})\in{\Gamma_{0}}$ in the
metric $\operatorname{d}_{\Gamma}$ can equivalently be stated as:

$\left[\gamma_{n}\stackrel{{\scriptstyle\operatorname{d}_{\Gamma}}}{{\to}}(h,\xi)\right]\iff\Big{[}h_{n}\to
h\text{ in }\mathscr{C}(\left[0,1\right])\text{ and for all $\varepsilon>0$,
}\xi_{n}\to\xi\text{ in }\mathscr{C}(\\{|h|\geq\varepsilon\\})\Big{]}$ (2.5)

Furthermore, for any
$\psi\in\mathscr{C}({\left[0,1\right]\times\overline{\Omega}})$, we have
$\Psi\in\mathscr{C}(\left[0,1\right]\times{\Gamma_{0}})$ where

$\forall t\in\left[0,1\right],\
(h,\xi)\in{\Gamma_{0}},\qquad\Psi(t,h,\xi)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}h(t)\psi(t,\xi(t)).$ (2.6)

###### Proof.

_Metric properties._ To see that
$(\nicefrac{{{\Gamma_{0}}}}{{\sim}},\operatorname{d}_{\Gamma})$ is a metric
space, it is sufficient to show that
$(\nicefrac{{\mathds{R}\times\overline{\Omega}}}{{\sim}},\operatorname{d}_{F})$
is a metric space. Following the argument of [Bredies et al., 2020a, Lemma
3.3] in the case $h\geq 0$, observe that for all $r_{1},r_{2}\in\mathds{R}$,
$x_{1},x_{2}\in\overline{\Omega}$

$\operatorname{d}_{F}((r_{1},x_{1}),(r_{2},x_{2}))=\sup_{c_{1},c_{2}\in\mathds{R}}\left\\{\,r_{1}c_{1}-r_{2}c_{2}\mid|c_{1}|,|c_{2}|\leq
1,\ |c_{1}-c_{2}|\leq|x_{1}-x_{2}|\,\right\\}.$ (2.7)

For $\min(r_{1},r_{2})\geq 0$, this has already been shown in Bredies et al.
[2020a]. The case $\max(r_{1},r_{2})\leq 0$ is equivalent, therefore it is
true for all $r_{1}r_{2}\geq 0$. For the final case, without loss of
generality $r_{1}\geq 0\geq r_{2}$, so the supremum is achieved by
$c_{1}=c_{2}=1$.

It is clear that $\operatorname{d}_{\Gamma}\geq 0$. To verify the triangle
inequality, let
$(r_{1},x_{1}),(r_{2},x_{2}),(r_{3},x_{3})\in\mathds{R}\times\overline{\Omega}$.
Let $c_{1},c_{3}\in\left[-1,1\right]$ such that
${\left|c_{1}-c_{3}\right|}\leq{\left|x_{1}-x_{3}\right|}$. By symmetry of
$\operatorname{d}_{F}$, we may assume without loss of generality that
$c_{1}\leq c_{3}$. Then

$\displaystyle
c_{3}-c_{1}\leq{\left|x_{1}-x_{3}\right|}\leq{\left|x_{2}-x_{1}\right|}+{\left|x_{3}-x_{2}\right|},$
(2.8)

hence there exists $c_{2}\in\left[c_{1},c_{3}\right]$ such that
${\left|c_{2}-c_{1}\right|}\leq{\left|x_{2}-x_{1}\right|}$ and
${\left|c_{2}-c_{3}\right|}\leq{\left|x_{2}-x_{3}\right|}$. Thus,

$\displaystyle
r_{1}c_{1}-r_{3}c_{3}=(r_{1}c_{1}-r_{2}c_{2})+(r_{2}c_{2}-r_{3}c_{3})\leq\operatorname{d}_{F}((r_{1},x_{1}),(r_{2},x_{2}))+\operatorname{d}_{F}((r_{2},x_{2}),(r_{3},x_{3})).$
(2.9)

Taking the supremum over $c_{1},c_{3}$ yields

$\displaystyle\operatorname{d}_{F}((r_{1},x_{1}),(r_{3},x_{3}))\leq\operatorname{d}_{F}((r_{1},x_{1}),(r_{2},x_{2}))+\operatorname{d}_{F}((r_{2},x_{2}),(r_{3},x_{3}).$
(2.10)

Hence $\operatorname{d}_{F}$ satisfies the triangle inequality and so does
$\operatorname{d}_{\Gamma}$.

Finally we must check
$\operatorname{d}_{\Gamma}(\gamma_{1},\gamma_{2})=0\iff\gamma_{1}\sim\gamma_{2}$.
To see $h_{1}=h_{2}$, observe for any $t\in\left[0,1\right]$ either
$|h_{1}(t)|+|h_{2}(t)|=0$, or $|h_{1}(t)-h_{2}(t)|=0$, but in any case,
$h_{1}(t)=h_{2}(t)$. Now whenever $h_{1}(t)\neq 0$, we must be in the case
$|\xi_{1}(t)-\xi_{2}(t)|\leq 2$, therefore
$|h_{1}(t)||\xi_{1}(t)-\xi_{2}(t)|=0$, hence $\xi_{1}(t)=\xi_{2}(t)$ as
required.

_Convergence._ Fix a sequence $\gamma_{n}=(h_{n},\xi_{n})\in{\Gamma_{0}}$ and
point $\gamma=(h,\xi)\in{\Gamma_{0}}$. Note that for all $r_{i}\in\mathds{R}$,
$x_{i}\in\overline{\Omega}$,

$|r_{1}-r_{2}|\leq\operatorname{d}_{F}((r_{1},x_{1}),(r_{2},x_{2}))\leq|r_{1}|+|r_{2}|.$
(2.11)

Suppose $\operatorname{d}_{\Gamma}(\gamma_{n},\gamma)\to 0$, then we
immediately have ${\left\lVert h_{n}-h\right\rVert}_{\infty}\to 0$. Also, for
any $\varepsilon>0$ choose $N_{\varepsilon}\in\mathds{N}$ such that

$\forall n\geq N_{\varepsilon},\qquad{\left\lVert
h_{n}-h\right\rVert}_{\infty}\leq\frac{\varepsilon}{2}\quad\text{and}\quad\operatorname{d}_{\Gamma}(\gamma_{n},\gamma)\leq\frac{\varepsilon}{2}.$
(2.12)

Observe that for all $t\in\\{|h|\geq\varepsilon\\}$ and $n\geq
N_{\varepsilon}$, if $|\xi_{n}(t)-\xi(t)|\geq 2$, then
$\operatorname{d}_{\Gamma}(\gamma_{n},\gamma)\geq|h_{n}(t)|+|h(t)|>\varepsilon$,
contradicting the choice of $N_{\varepsilon}$. Therefore we have the uniform
bound

$\sup_{|h(t)|\geq\varepsilon}|\xi_{n}(t)-\xi(t)|\leq\sup_{|h(t)|\geq\varepsilon}\frac{2}{\varepsilon}\min(|h_{n}(t)|,|h(t)|)|\xi_{n}(t)-\xi(t)|\leq\frac{2}{\varepsilon}\operatorname{d}_{\Gamma}(\gamma_{n},\gamma).$
(2.13)

This concludes the “$\implies$” direction of (2.5). Conversely, suppose
${\left\lVert h_{n}-h\right\rVert}_{\infty}\to 0$ and for any $\varepsilon>0$,
${\left\lVert\xi_{n}-\xi\right\rVert}_{L^{\infty}(\\{|h|\geq\varepsilon\\})}\to
0$. Fix $\varepsilon>0$ and $N_{\varepsilon}\in\mathds{N}$ such that for all

$\forall n\geq N_{\varepsilon},\qquad{\left\lVert
h_{n}-h\right\rVert}_{\infty}\leq\frac{\varepsilon}{2}\quad\text{and}\quad{\left\lVert\xi_{n}-\xi\right\rVert}_{L^{\infty}(\\{|h|\geq\varepsilon\\})}\leq\min(2,\varepsilon).$
(2.14)

Then, from (2.11), for all $n\geq N_{\varepsilon}$ we have

$\displaystyle\operatorname{d}_{\Gamma}(\gamma_{n},\gamma)$
$\displaystyle\leq\sup_{t\in\left[0,1\right]}\left\\{\begin{array}[]{l}|h_{n}(t)|+|h(t)|\\\
\operatorname{d}_{F}((h_{n}(t),\xi_{n}(t)),(h(t),\xi(t)))\end{array}\right.$
$\displaystyle\begin{array}[]{c}\text{if }|h(t)|,\\\ \text{else}\end{array}$
(2.19)
$\displaystyle\leq\sup_{t\in\left[0,1\right]}\left\\{\begin{array}[]{l}5\varepsilon/2\\\
\varepsilon/2+{\left\lVert
h\right\rVert}_{\infty}\varepsilon\end{array}\right.$
$\displaystyle\begin{array}[]{c}\text{if }|h(t)|,\\\ \text{else.}\end{array}$
(2.24)

In either case we have
$\limsup_{n\to+\infty}\operatorname{d}_{\Gamma}(\gamma_{n},\gamma)\leq
O(\varepsilon)$ for all $\varepsilon>0$, therefore
$\operatorname{d}_{\Gamma}(\gamma_{n},\gamma)\to 0$ as required.

_Separability._ Recall by the Weierstrass approximation theorem that
$\mathscr{C}(\left[0,1\right];{\left[0,1\right]\times\overline{\Omega}})$ is
separable, let $\\{\gamma_{n}=(h_{n},\xi_{n})\\}_{n\in\mathds{N}}$ be a dense
set (in the $\mathscr{C}$ sense). Fix $\gamma=(h,\xi)\in{\Gamma_{0}}$. Denote
the closed domains $U_{k}=\\{k|h|\geq 1\\}$. By the Urysohn-Brouwer lemma, for
each $k\in\mathds{N}$ there exists an extension of $\xi$, say
$\widehat{\xi}_{k}\in\mathscr{C}(\left[0,1\right];\overline{\Omega})$, such
that $\xi=\widehat{\xi}_{k}$ on $U_{k}$. Therefore there exists a sequence
$n_{k}$ such that

$\forall k\in\mathds{N},\qquad{\left\lVert
h_{n_{k}}-h\right\rVert}_{\infty}\leq\frac{1}{k}\quad\text{and}\quad{\left\lVert\xi_{n_{k}}-\xi\right\rVert}_{L^{\infty}(\\{k|h|\geq
1\\})}\leq\frac{1}{k}.$ (2.25)

From (2.5), it is clear that
$\gamma_{n_{k}}\stackrel{{\scriptstyle\operatorname{d}_{\Gamma}}}{{\to}}\gamma$,
therefore $\\{\gamma_{n}\\}_{n\in\mathds{N}}$ is also countable and dense in
${\Gamma_{0}}$.

_Completeness._ Let $\gamma_{n}=(h_{n},\xi_{n})$ be an arbitrary
$\operatorname{d}_{\Gamma}$-Cauchy sequence, we are required to show that
there exists $\gamma\in{\Gamma_{0}}$ such that
$\operatorname{d}_{\Gamma}(\gamma_{n},\gamma)\to 0$. By the lower bound in
(2.11), we have

${\left\lVert
h_{n}-h_{m}\right\rVert}_{\infty}\leq\operatorname{d}_{\Gamma}(\gamma_{n},\gamma_{m}),$
(2.26)

so $h_{n}$ is Cauchy in $\mathscr{C}(\left[0,1\right])$. As
$\mathscr{C}(\left[0,1\right])$ is complete, there must be a limit point
$h\in\mathscr{C}(\left[0,1\right])$ such that
$h_{n}\stackrel{{\scriptstyle\mathscr{C}}}{{\to}}h$. Repeating the argument in
(2.13) confirms that $\xi_{n}$ is also Cauchy in
$\mathscr{C}(\\{|h|\geq\varepsilon\\};\overline{\Omega})$ for each
$\varepsilon>0$, therefore there is a (unique) curve
$\xi\in\mathscr{C}(\\{h\neq 0\\};\overline{\Omega})$ such that $\xi_{n}\to\xi$
in each $\mathscr{C}(\\{|h|\geq\varepsilon\\};\overline{\Omega})$,
$\varepsilon>0$. We may extend $\xi$ by a constant $x_{0}\in\Omega$ on
$\\{h=0\\}$. Due to (2.5), we obtain that $\gamma_{n}\to(h,\xi)$.

_Continuity._ Fix
$\psi\in\mathscr{C}({\left[0,1\right]\times\overline{\Omega}})$ and define
$\Psi$ as in (2.6). As we have now confirmed that
$(\nicefrac{{{\Gamma_{0}}}}{{\sim}},\operatorname{d}_{\Gamma})$ is a metric
space, we can use the sequential definitions of continuity. Suppose
$\tau_{n}\in\left[0,1\right]$, $\gamma_{n}=(h_{n},\xi_{n})\in{\Gamma_{0}}$ and
$\tau_{n}\to t$,
$\gamma_{n}\stackrel{{\scriptstyle\operatorname{d}_{\Gamma}}}{{\to}}\gamma$,
observe for each $n$ we have

$\displaystyle|\Psi(\tau_{n},\gamma_{n})-\Psi(t,\gamma)|$
$\displaystyle=|h_{n}(\tau_{n})\psi(\tau_{n},\xi_{n}(\tau_{n}))-h(t)\psi(t,\xi(t))|$
(2.27)
$\displaystyle=|(h_{n}(\tau_{n})-h(t)+h(t))\psi(\tau_{n},\xi_{n}(\tau_{n}))-h(t)\psi(t,\xi(t))|$
(2.28)
$\displaystyle\leq|h_{n}(\tau_{n})-h(t)|{\left\lVert\psi\right\rVert}_{\infty}+|h(t)||\psi(\tau_{n},\xi_{n}(\tau_{n}))-\psi(t,\xi(t))|.$
(2.29) $\displaystyle\leq[{\left\lVert
h_{n}-h\right\rVert}_{\infty}+|h(\tau_{n})-h(t)|]{\left\lVert\psi\right\rVert}_{\infty}+|h(t)||\psi(\tau_{n},\xi_{n}(\tau_{n}))-\psi(t,\xi(t))|.$
(2.30)

Now consider the limit of $n\to+\infty$. As
$h\in\mathscr{C}(\left[0,1\right])$, $h(\tau_{n})\to h(t)$. The
characterisation of convergence in $\operatorname{d}_{\Gamma}$ in (2.5) also
confirms ${\left\lVert h_{n}-h\right\rVert}_{\infty}\to 0$ and either
$h(t)=0$, or

$|\xi_{n}(\tau_{n})-\xi(t)|\leq|\xi_{n}(\tau_{n})-\xi(\tau_{n})|+|\xi(\tau_{n})-\xi(t)|\to
0.$ (2.31)

In either case, we see that

$\limsup_{n\to+\infty}|\Psi(\tau_{n},\gamma_{n})-\Psi(t,\gamma)|\leq\limsup_{n\to+\infty}|h(t)||\psi(\tau_{n},\xi_{n}(\tau_{n}))-\psi(t,\xi(t))|=0,$
(2.32)

therefore $\Psi$ is continuous at $(t,\gamma)$. ∎

A key analytical tool in related prior works (c.f. Bredies et al. [2021],
Bredies and Fanzon [2020], Bredies et al. [2020a, b]) is a mapping between
measures $\rho\in\mathcal{M}({\left[0,1\right]\times\overline{\Omega}})$ which
satisfy the continuity equation, and measures on (weighted) paths
$\sigma\in\mathcal{M}({\Gamma_{0}})$, making some structures become more
apparent. We will now recap and expand on previous results. The first theorem
collects results from Ambrosio et al. [2008], Bredies et al. [2020a] and
provides minor extensions for the scope of the current work.

######  Theorem 2.2.

Let $\sigma\in\mathcal{M}({\Gamma_{0}})$. If $\int_{\Gamma_{0}}{\left\lVert
h\right\rVert}_{1}\mathop{}\\!\mathrm{d}|\sigma|(h,\xi)<+\infty$, then there
is a unique finite Borel measure
$\Theta(\sigma)\in\mathcal{M}({\left[0,1\right]\times\overline{\Omega}})$ such
that

$\forall\psi\in\mathscr{C}({\left[0,1\right]\times\overline{\Omega}}),\qquad\int_{\left[0,1\right]\times\overline{\Omega}}\psi(t,x)\mathop{}\\!\mathrm{d}\Theta(\sigma)(t,x)=\int_{\Gamma_{0}}\left(\int_{0}^{1}h(t)\psi(t,\xi(t))\mathop{}\\!\mathrm{d}t\right)\mathop{}\\!\mathrm{d}\sigma(h,\xi).$
(2.33)

Moreover,

  1. 1.

The mapping
$\Theta\colon\left\\{\,\sigma\in\mathcal{M}({\Gamma_{0}})\mid\int_{\Gamma_{0}}{\left\lVert
h\right\rVert}_{1}\mathop{}\\!\mathrm{d}|\sigma|<+\infty\,\right\\}\to\mathcal{M}({\left[0,1\right]\times\overline{\Omega}})$
is linear.

  2. 2.

Equality (2.33) holds for all $\psi\in
L^{1}_{|\Theta(\sigma)|}({\left[0,1\right]\times\overline{\Omega}})$.

  3. 3.

If $\int_{\Gamma_{0}}{\left\lVert
h\right\rVert}_{\infty}\mathop{}\\!\mathrm{d}{\left|\sigma\right|}<+\infty$,
then
$\Theta(\sigma)\in\mathscr{C}_{w}(\left[0,1\right];\mathcal{M}(\overline{\Omega}))$.

  4. 4.

If there exist Borel measurable functions
$v\colon{\left[0,1\right]\times\overline{\Omega}}\rightarrow\mathds{R}^{d}$
and $g:{\left[0,1\right]\times\overline{\Omega}}\rightarrow\mathds{R}$ such
that

$\displaystyle h^{\prime}(t)=g(t,\xi(t))h(t)\text{ for $\sigma$-a.e. $(h,\xi)$
and a.e. $t\in\left]0,1\right[$},$ (2.34)
$\displaystyle\xi^{\prime}(t)=v(t,\xi(t))\text{ for $\sigma$-a.e. $(h,\xi)$
and a.e. $t$ such that $h(t)>0$},$ (2.35) and
$\displaystyle\int_{\Gamma_{0}}\int_{0}^{1}(1+|v(t,\xi(t))|+|g(t,\xi(t))|)\,|h(t)|\,\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}{\left|\sigma\right|}(h,\xi)<+\infty,$
(2.36)

then $\int_{\Gamma_{0}}{\left\lVert
h\right\rVert}_{\infty}\mathop{}\\!\mathrm{d}{\left|\sigma\right|}<+\infty$
and $\Theta(\sigma)$ satisfies the continuity equation (1.9).

Conversely, given
$\rho\in\mathcal{M}({\left[0,1\right]\times\overline{\Omega}})$, if $\rho\geq
0$ satisfies the continuity equation (1.9) and

$\int_{\left[0,1\right]\times\overline{\Omega}}(1+|v(t,x)|^{2}+|g(t,x)|^{2})\mathop{}\\!\mathrm{d}\rho(t,x)<+\infty,$
(2.37)

then $\rho=\Theta(\sigma)$ for some $\sigma\in\mathcal{M}^{+}({\Gamma_{0}})$
such that (2.34)-(2.36) hold and $\int_{\Gamma_{0}}{\left\lVert
h\right\rVert}_{\infty}\mathop{}\\!\mathrm{d}\sigma<+\infty$.

###### Proof.

By Lemma 2.1, for any
$\psi\in\mathscr{C}({\left[0,1\right]\times\overline{\Omega}})$, the map

$(h,\xi)\mapsto\int_{0}^{1}h(t)\psi(t,\xi(t))\mathop{}\\!\mathrm{d}t=\int_{0}^{1}\Psi(t,h,\xi)\mathop{}\\!\mathrm{d}t$
(2.38)

is continuous (hence Borel) in ${\Gamma_{0}}$ and dominated by ${\left\lVert
h\right\rVert}_{1}{\left\lVert\psi\right\rVert}_{\infty}$, therefore the
right-hand side of (2.33) is well-defined. This induces a linear form on
$\mathscr{C}({\left[0,1\right]\times\overline{\Omega}})$ which is moreover
bounded, since

$\left|\int_{\Gamma_{0}}\int_{0}^{1}h(t)\psi(t,\xi(t))\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}\sigma(h,\xi)\right|\leq{\left\lVert\psi\right\rVert}_{\infty}\int_{\Gamma_{0}}{\left\lVert
h\right\rVert}_{1}\mathop{}\\!\mathrm{d}{\left|\sigma\right|}.$ (2.39)

By the Riesz representation theorem, that linear form is represented by a
unique Radon measure
$\Theta(\sigma)\in\mathcal{M}({\left[0,1\right]\times\overline{\Omega}})$.
This confirms that $\Theta$ satisfies the required properties for point (1).

Points (2)-(4) have been proved in Bredies et al. [2020a] under the additional
assumptions $\sigma\geq 0$ and $\inf_{t\in\left[0,1\right]}h(t)\geq 0$. To
apply these results, we need the following modified Hahn-Jordan decomposition.

> Claim: For any $\sigma\in\mathcal{M}({\Gamma_{0}})$ with
> $\int_{\Gamma_{0}}{\left\lVert
> h\right\rVert}_{1}\mathop{}\\!\mathrm{d}{\left|\sigma\right|}<+\infty$ there
> exists $\sigma^{+},\sigma^{-}\in\mathcal{M}^{+}({\Gamma_{0}})$ such that
>
>
> $\sigma^{\pm}\left(\left\\{\,(h,\xi)\in{\Gamma_{0}}\mid\inf_{t\in\left[0,1\right]}h(t)<0\,\right\\}\right)=0$
> (2.40)
>
> with $\int_{{\Gamma_{0}}}{\left\lVert
> h\right\rVert}_{1}\mathop{}\\!\mathrm{d}\sigma^{\pm}<+\infty$ and
> $\Theta(\sigma)=\Theta(\sigma^{+})-\Theta(\sigma^{-})$.
>
> ###### Proof of claim..
>
> The Hahn-Jordan decomposition gives $\sigma=\max(0,\sigma)-\max(0,-\sigma)$,
> therefore it is sufficient to consider the case $\sigma\geq 0$. Define the
> maps $T^{\pm}\colon{\Gamma_{0}}\rightarrow{\Gamma_{0}}$ by
>
> $\forall(h,\xi)\in{\Gamma_{0}},\qquad
> T^{\pm}(h,\xi)\stackrel{{\scriptstyle\mbox{\tiny def.}}}{{=}}(\max(0,\pm
> h),\xi).$ (2.41)
>
> Because
> $\operatorname{d}_{\Gamma}(T^{\pm}(\gamma_{1}),T^{\pm}(\gamma_{2}))\leq\operatorname{d}_{\Gamma}(\gamma_{1},\gamma_{2})$,
> $T^{\pm}$ are continuous (therefore Borel), and we can define the image
> measures $\sigma^{\pm}\stackrel{{\scriptstyle\mbox{\tiny
> def.}}}{{=}}(T^{\pm})_{\sharp}\sigma$. Since the push-forward operation does
> not increase the total variation, we have
> $\sigma^{+},\sigma^{-}\in\mathcal{M}^{+}({\Gamma_{0}})$. Moreover, (2.40)
> holds, and by construction of the image measure
>
>
> $\int_{{\Gamma_{0}}}\phi(h,\xi)\mathop{}\\!\mathrm{d}\sigma^{\pm}=\int_{\Gamma_{0}}\phi(\max(0,\pm
> h),\xi)\mathop{}\\!\mathrm{d}\sigma$ (2.42)
>
> for all $\phi\in\mathscr{C}_{b}({\Gamma_{0}})$, and by monotone convergence
>
> $\int_{{\Gamma_{0}}}{\left\lVert
> h\right\rVert}_{1}\mathop{}\\!\mathrm{d}\sigma^{\pm}=\int_{\Gamma_{0}}{\left\lVert\max(0,\pm
> h)\right\rVert}_{1}\mathop{}\\!\mathrm{d}\sigma\leq\int_{\Gamma_{0}}{\left\lVert
> h\right\rVert}_{1}\mathop{}\\!\mathrm{d}\sigma.$ (2.43)
>
> Finally, we confirm that for all
> $\psi\in\mathscr{C}({\left[0,1\right]\times\overline{\Omega}})$,
>
>
> $\displaystyle\int_{\left[0,1\right]\times\overline{\Omega}}\psi(t,x)\mathop{}\\!\mathrm{d}\Theta(\sigma)$
> $\displaystyle=\int_{\Gamma_{0}}\int_{0}^{1}h(t)\psi(t,\xi(t))\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}\sigma$
> (2.44)
> $\displaystyle=\int_{\Gamma_{0}}\int_{0}^{1}[\max(0,h(t))-\max(0,-h(t))]\psi(t,\xi(t))\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}\sigma$
> (2.45)
> $\displaystyle=\int_{{\Gamma_{0}}}\int_{0}^{1}h(t)\psi(t,\xi(t))\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}(\sigma^{+}-\sigma^{-})$
> (2.46)
> $\displaystyle=\int_{\left[0,1\right]\times\overline{\Omega}}\psi(t,x)\mathop{}\\!\mathrm{d}(\Theta(\sigma^{+})-\Theta(\sigma^{-}))$
> (2.47)
>
> as required. ∎

With this decomposition, we can apply the results of Bredies et al. [2020a] to
each $\sigma^{\pm}$. Point (2) is a direct consequence of [Bredies et al.,
2020a, Lemma 4.4]. Under the assumptions of point (3), the same lemma also
guarantees the existence of a disintegration $\rho_{t}^{\pm}$ of
$\Theta(\sigma^{\pm})$, in the sense of (1.7). In particular,

$\forall\psi\in\mathscr{C}(\overline{\Omega}),\
t\in\left[0,1\right],\qquad\int_{\overline{\Omega}}\psi(x)\mathop{}\\!\mathrm{d}\rho_{t}^{\pm}(x)=\int_{{\Gamma_{0}}}h(t)\psi(\xi(t))\mathop{}\\!\mathrm{d}\sigma^{\pm}(h,\xi).$
(2.48)

The same property holds for $\Theta(\sigma)$ by linearity. Point (3) requires
$t\mapsto\int_{\overline{\Omega}}\psi\mathop{}\\!\mathrm{d}(\rho_{t}^{+}-\rho^{-}_{t})$
to be continuous for all $\psi\in\mathscr{C}(\overline{\Omega})$. By Lemma
2.1, for each $t\in\left[0,1\right]$ the function
$\Psi(t,h,\xi)\stackrel{{\scriptstyle\mbox{\tiny def.}}}{{=}}h(t)\psi(\xi(t))$
is continuous on $\left[0,1\right]\times{\Gamma_{0}}$ and

$\left|\int_{\overline{\Omega}}\psi\mathop{}\\!\mathrm{d}(\rho_{\tau}-\rho_{t})\right|\leq\int_{\Gamma_{0}}|h(\tau)\psi(\xi(\tau))-h(t)\psi(\xi(t))|\mathop{}\\!\mathrm{d}{\left|\sigma\right|}(h,\xi)=\int_{\Gamma_{0}}|\Psi(\tau,\gamma)-\Psi(t,\gamma)|\mathop{}\\!\mathrm{d}{\left|\sigma\right|}(\gamma).$
(2.49)

The integrand is pointwise bounded by $2{\left\lVert
h\right\rVert}_{\infty}{\left\lVert\psi\right\rVert}_{\infty}$, therefore we
conclude that the limit of the integral as $\tau\to t$ is 0 by dominated
convergence.

Point (4) and its converse are addressed by [Bredies et al., 2020a, Thm. 4.2].
For the forward direction, we again consider the modified Hahn-Jordan
decomposition. The assumptions (2.34)-(2.36) are only assumed to hold for
almost every $t$ and $h$, therefore they also hold for the choice of measures
$\sigma^{\pm}$ given in the claim replacing $h$ with $\max(0,\pm h)$. We can
then apply [Thm. 4.2 Bredies et al., 2020a] to each component and sum for the
result. The converse is exactly the statement of [Thm. 4.2 Bredies et al.,
2020a] because we assume $\rho\geq 0$. ∎

The previous theorem makes precise statements about when $\Theta(\sigma)$
satisfies particular smoothness properties, the assumption
${\left\lVert\sigma\right\rVert}<+\infty$ is too weak for any useful
properties. In this work we would like $\Theta$ to be narrowly continuous and
only to consider measures $\sigma$ such that
$\Theta(\sigma)\in\mathscr{C}_{w}(\left[0,1\right];\mathcal{M}(\overline{\Omega}))$.
The following lemma shows that we can choose a subset of ${\Gamma_{0}}$ which
only removes un-desired elements from the image of $\Theta$. Recall
${\left\lVert\sigma\right\rVert}<+\infty$ for all
$\sigma\in\mathcal{M}({\Gamma_{0}})$, if $\Gamma\subset{\Gamma_{0}}$ is
sufficiently “small”, then this implicit assumption is sufficient for
$\Theta|_{\mathcal{M}(\Gamma)}$ to have the desired properties.

######  Lemma 2.3.

Define the set

$\quad\text{for some}\quad
p\in\left[1,+\infty\right],\qquad\Gamma_{p}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left\\{\,\gamma=(h,\xi)\in{\Gamma_{0}}\mid{\left\lVert
h\right\rVert}_{p}\leq 1\,\right\\},$ (2.50)

then

$\left\\{\,\Theta(\sigma)\mid\sigma\in\mathcal{M}({\Gamma_{0}}),\
\int_{{\Gamma_{0}}}{\left\lVert
h\right\rVert}_{p}\mathop{}\\!\mathrm{d}{\left|\sigma\right|}<+\infty\,\right\\}=\left\\{\,\Theta(\widehat{\sigma})\mid\widehat{\sigma}\in\mathcal{M}(\Gamma_{p})\,\right\\}$
(2.51)

and
$\Theta\colon\mathcal{M}(\Gamma_{p})\to\mathcal{M}({\left[0,1\right]\times\overline{\Omega}})$
is narrowly continuous.

Furthermore, if $p=+\infty$, then $\forall t\in\left[0,1\right]$,
$(e_{t})_{\sharp}{\colon}\mathcal{M}(\Gamma_{\infty})\to\mathcal{M}(\overline{\Omega})$
is also narrowly continuous.

In particular, the sequentially we have that, for any sequence
$\sigma^{n}\xrightharpoonup{*}\sigma$ narrowly in $\mathcal{M}(\Gamma_{p})$:

$\displaystyle\text{for all }p\in\left[1,+\infty\right],$
$\displaystyle\qquad\Theta(\sigma^{n})\stackrel{{\scriptstyle*}}{{\rightharpoonup}}\Theta(\sigma)\text{
narrowly in }\mathcal{M}({\left[0,1\right]\times\overline{\Omega}}),$ (2.52)
$\displaystyle\text{if }p=+\infty,\ \forall t\in\left[0,1\right],$
$\displaystyle\qquad(e_{t})_{\sharp}{\sigma^{n}}\stackrel{{\scriptstyle*}}{{\rightharpoonup}}(e_{t})_{\sharp}{\sigma}\text{
narrowly in }\mathcal{M}(\overline{\Omega}).$ (2.53)

###### Proof.

The “$\supset$” inclusion is clear: if
$\widehat{\sigma}\in\mathcal{M}(\Gamma_{p})$, extend $\widehat{\sigma}$ by 0
such that $\widehat{\sigma}\in\mathcal{M}({\Gamma_{0}})$, then

$\int_{{\Gamma_{0}}}{\left\lVert
h\right\rVert}_{p}\mathop{}\\!\mathrm{d}{\left|\widehat{\sigma}\right|}\leq\int_{\Gamma_{p}}1\mathop{}\\!\mathrm{d}{\left|\widehat{\sigma}\right|}={\left\lVert\widehat{\sigma}\right\rVert}<+\infty$
(2.54)

as required. Conversely, suppose $\int_{{\Gamma_{0}}}{\left\lVert
h\right\rVert}_{p}\mathop{}\\!\mathrm{d}{\left|\sigma\right|}<+\infty$. Note
as $h\mapsto\max(1,{\left\lVert h\right\rVert}_{p})$ and
$h\mapsto\frac{h}{\max(1,{\left\lVert h\right\rVert}_{p})}$ are continuous in
$\mathscr{C}(\left[0,1\right])$ and
$\operatorname{d}_{\Gamma}(\gamma_{1},\gamma_{2})\geq{\left\lVert
h_{1}-h_{2}\right\rVert}_{\infty}$, both $(h,\xi)\mapsto\max(1,{\left\lVert
h\right\rVert}_{p})$ and $(h,\xi)\mapsto(\frac{h}{\max(1,{\left\lVert
h\right\rVert}_{p})},\xi)$ are continuous w.r.t. $\operatorname{d}_{\Gamma}$.
Define $\widehat{\sigma}$ through the (rescaled) push-forward
$\widehat{\sigma}(h,\xi)=\max(1,{\left\lVert h\right\rVert}_{p})\cdot
T_{\sharp}\sigma(h,\xi)$ where $T\colon{\Gamma_{0}}\to\Gamma_{p}$ is defined
by the map $T(h,\xi)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left(\frac{h}{\max(1,{\left\lVert
h\right\rVert}_{p})},\xi\right)$. $T$ is Borel measurable, therefore
$\widehat{\sigma}$ is a Borel measure which satisfies

$\forall\phi\in\mathscr{C}_{b}(\Gamma_{p}),\qquad\int_{\Gamma_{p}}\phi\mathop{}\\!\mathrm{d}\widehat{\sigma}=\int_{{\Gamma_{0}}}\max(1,{\left\lVert
h\right\rVert}_{p})\phi\left(\frac{h}{\max(1,{\left\lVert
h\right\rVert}_{p})},\xi\right)\mathop{}\\!\mathrm{d}\sigma(h,\xi).$ (2.55)

Furthermore, $\widehat{\sigma}\in\mathcal{M}(\Gamma_{p})$ as

$\left|\int_{\Gamma_{p}}\phi\mathop{}\\!\mathrm{d}\widehat{\sigma}\right|\leq{\left\lVert\phi\right\rVert}_{\infty}\int_{{\Gamma_{0}}}(1+{\left\lVert
h\right\rVert}_{p})\mathop{}\\!\mathrm{d}{\left|\sigma\right|}(h,\xi),$ (2.56)

so ${\left\lVert\widehat{\sigma}\right\rVert}<+\infty$. The last step is to
confirm $\Theta(\widehat{\sigma})=\Theta(\sigma)$. For any
$\psi\in\mathscr{C}({\left[0,1\right]\times\overline{\Omega}})$ observe

$\displaystyle\int_{\left[0,1\right]\times\overline{\Omega}}\psi(t,x)\mathop{}\\!\mathrm{d}\Theta(\widehat{\sigma})(t,x)$
$\displaystyle=\int_{\Gamma_{p}}\left(\int_{0}^{1}h(t)\psi(t,\xi(t))\mathop{}\\!\mathrm{d}t\right)\mathop{}\\!\mathrm{d}\widehat{\sigma}(h,\xi)$
(2.57) $\displaystyle=\int_{{\Gamma_{0}}}\max(1,{\left\lVert
h\right\rVert}_{p})\left(\int_{0}^{1}\frac{h(t)}{\max(1,{\left\lVert
h\right\rVert}_{p})}\psi(t,\xi(t))\mathop{}\\!\mathrm{d}t\right)\mathop{}\\!\mathrm{d}\sigma(h,\xi)$
(2.58)
$\displaystyle=\int_{\left[0,1\right]\times\overline{\Omega}}\psi(t,x)\mathop{}\\!\mathrm{d}\Theta(\sigma)(t,x).$
(2.59)

This shows the “$\subset$” direction, therefore we conclude the equality in
(2.51). For continuity of $\Theta$, in the $p\geq 1$ context [Thm. 1, Sec. 1.6
Yosida, 1980] states that $\Theta$ is continuous in
$(\mathscr{C}({\left[0,1\right]\times\overline{\Omega}}))^{\prime}$ if and
only if for every
$\psi\in\mathscr{C}({\left[0,1\right]\times\overline{\Omega}})$ there exists a
finite collection $\Phi\subset\mathscr{C}_{b}(\Gamma_{p})$ such that

$\forall\sigma\in\mathcal{M}(\Gamma_{p}),\qquad\left|\int_{\left[0,1\right]\times\overline{\Omega}}\psi(t,x)\mathop{}\\!\mathrm{d}\Theta(\sigma)\right|\leq\max_{\phi\in\Phi}\left|\int_{\Gamma_{p}}\phi(h,\xi)\mathop{}\\!\mathrm{d}\sigma(h,\xi)\right|.$
(2.60)

From Lemma 2.1, define $\Psi(t,h,\xi)=h(t)\psi(t,\xi(t))$, then
$\phi(h,\xi)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\int_{0}^{1}\Psi(t,h,\xi)\mathop{}\\!\mathrm{d}t$ is continuous
and bounded as
$|\phi(h,\xi)|\leq{\left\lVert\psi\right\rVert}_{\infty}{\left\lVert
h\right\rVert}_{1}\leq{\left\lVert\psi\right\rVert}_{\infty}$ for each
$(h,\xi)\in\Gamma_{p}$, $p\geq 1$. This confirms the continuity of $\Theta$
because

$\int_{\left[0,1\right]\times\overline{\Omega}}\psi(t,x)\mathop{}\\!\mathrm{d}\Theta(\sigma)=\int_{\Gamma_{p}}\left(\int_{0}^{1}h(t)\psi(t,\xi(t))\mathop{}\\!\mathrm{d}t\right)\mathop{}\\!\mathrm{d}\sigma(h,\xi)=\int_{\Gamma_{p}}\phi(h,\xi)\mathop{}\\!\mathrm{d}\sigma(h,\xi)$
(2.61)

for each $\sigma\in\mathcal{M}(\Gamma_{p})$, therefore $\Phi=\\{\phi\\}$ is
sufficient. The $p=+\infty$ case is special because for each
$t\in\left[0,1\right]$, the map $\phi(h,\xi)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\Psi(t,h,\xi)$ is continuous and bounded by
${\left\lVert\psi\right\rVert}_{\infty}$. This leads to

$\int_{\overline{\Omega}}\psi(t,x)\mathop{}\\!\mathrm{d}(e_{t})_{\sharp}{\sigma}=\int_{\Gamma_{\infty}}\left(h(t)\psi(t,\xi(t))\right)\mathop{}\\!\mathrm{d}\sigma(h,\xi)=\int_{\Gamma_{\infty}}\phi(h,\xi)\mathop{}\\!\mathrm{d}\sigma(h,\xi),$
(2.62)

so we conclude similarly that $(e_{t})_{\sharp}{}$ is continuous in
$(\mathscr{C}(\overline{\Omega}))^{\prime}$ as required. ∎

In the rest of this work we want to work on a domain
$D\subset\mathcal{M}^{+}({\Gamma_{0}})$ such that $\Theta|_{D}$ has nice
properties with respect to the space
$\mathscr{C}_{w}(\left[0,1\right];\mathcal{M}(\overline{\Omega}))$. In
particular, Theorem 2.2 shows that we may choose
$D\subset\left\\{\,\sigma\in\mathcal{M}^{+}({\Gamma_{0}})\mid\int_{{\Gamma_{0}}}{\left\lVert
h\right\rVert}_{\infty}\mathop{}\\!\mathrm{d}\sigma<\infty\,\right\\}$. Lemma
2.3 shows that it is equivalent to consider
$D\subset\mathcal{M}^{+}(\Gamma_{\infty})$, from a numerical perspective there
is very little difference. Analytically we will always choose
$D\subset\mathcal{M}^{+}(\Gamma)$ for $\Gamma$ a closed subset of
$\Gamma_{\infty}$ in order to have the additional continuity properties.

######  Remark 2.1.

Note that the definition of $\Gamma_{p}$ is consistent with the relation
$\sim$, so $(\nicefrac{{\Gamma_{p}}}{{\sim}},\operatorname{d}_{\Gamma})$ is
also a metric space. In the rest of the paper, we require that $\Gamma$ si a
closed subset of $\nicefrac{{\Gamma_{\infty}}}{{\sim}}$ in order for it to be
a complete separable metric space.

##  3 The balanced Benamou-Brenier example

The purpose of this section is to motivate our framework using a specific
example which has been studied in the seminal work of Bredies et al. [2021,
2020b]. Following Bredies et al. [2021, 2020b], we define for $\alpha,\beta>0$
the Benamou-Brenier penalty
$\mathcal{W}\colon\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})\to\mathds{R}\cup\\{+\infty\\}$
by

$\mathcal{W}(\rho)\stackrel{{\scriptstyle\mbox{\tiny def.}}}{{=}}\inf_{v\in
L^{2}_{\rho}(\left[0,1\right]\times\overline{\Omega};\mathds{R}^{d})}\Bigg{\\{}\int_{\left[0,1\right]\times\overline{\Omega}}\left[\alpha+\frac{\beta}{2}|v|^{2}\right]\mathop{}\\!\mathrm{d}\rho\
\operatorname{\;s.t.\;}\
\partial_{t}\rho+\operatorname{div}(v\rho)=0\Bigg{\\}},$ (3.1)

where the continuity equation is satisfied in the sense of (1.9). The function
with $\alpha=0$, $\beta=1$ is referred to as the Benamou-Brenier energy whose
main properties can be found in [Santambrogio, 2015, Sec. 5.3.1]. By Theorem
2.2, if $\mathcal{W}(\rho)<+\infty$, then
$\rho\in\mathscr{C}_{w}(\left[0,1\right];\mathcal{M}(\overline{\Omega}))$. In
particular, there is a curve $t\mapsto\rho_{t}$ in
$\mathcal{M}^{+}(\overline{\Omega})$, $0\leq t\leq 1$, which is continuous in
the narrow topology and such that

$\forall\psi\in
L^{1}(\left[0,1\right]\times\overline{\Omega};\rho),\quad\int_{\left[0,1\right]\times\overline{\Omega}}\psi(t,x)\mathop{}\\!\mathrm{d}\rho(t,x)=\int_{0}^{1}\left(\int_{\overline{\Omega}}\psi(t,x)\mathop{}\\!\mathrm{d}\rho_{t}(x)\right)\mathop{}\\!\mathrm{d}t.$
(3.2)

Thanks to this property, given times $t_{0}=0<\ldots<t_{T}=1$, it is justified
to consider “slices” $\rho_{t_{j}}$ and to assume that we are given data
$b_{j}\in\mathds{R}^{m}$ (for simplicity assume $m$ does not change with $j$)
at time $t_{j}$, corresponding to linear observations of $\rho_{t_{j}}$, given
by linear operators
$A_{j}\colon\mathcal{M}(\overline{\Omega})\to\mathds{R}^{m}$. To solve the
corresponding inverse problem, [Bredies et al., 2021, (43)] proposes the
minimisation of

$\mathcal{E}(\rho)=\frac{1}{2}\sum_{j=0}^{T}{\left\lVert
A_{j}\rho_{t_{j}}-b_{j}\right\rVert}_{2}^{2}+\mathcal{W}(\rho).$ (3.3)

We assume that each $A_{j}$ is narrowly continuous (that is, in the topology
induced by $\mathscr{C}_{b}(\overline{\Omega})$). The more involved case of
continuous-time observations is handled in Bredies and Fanzon [2020].

Minimising such a model comes with very clear reconstruction properties.

######  Theorem 3.1 ([Bredies et al., 2021, Thms. 7, 10]).

Let $\alpha,\beta>0$ and

$\Gamma\stackrel{{\scriptstyle\mbox{\tiny def.}}}{{=}}\left\\{\,(h,\xi)\mid
h\equiv 1,\
\xi\in\operatorname{AC}^{2}(\left[0,1\right];\overline{\Omega})\,\right\\}.$
(3.4)

Then for all
$\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})$ with
$\mathcal{W}(\rho)<+\infty$, there exists $\sigma\in\mathcal{M}^{+}(\Gamma)$
such that $\rho=\Theta(\sigma)$, as defined in Theorem 2.2, and $v\in
L^{2}_{\rho}({\left[0,1\right]\times\overline{\Omega}};\mathds{R}^{d})$ such
that

$\mathcal{W}(\rho)=\int_{\left[0,1\right]\times\overline{\Omega}}\left[\alpha+\frac{\beta}{2}|v|^{2}\right]\mathop{}\\!\mathrm{d}\rho,\qquad\partial_{t}\rho+\operatorname{div}(v\rho)=0,$
(3.5)

and

$\text{for a.e. }t\in\left[0,1\right]\text{ and }\sigma\text{-a.e.
}\xi,\qquad\xi^{\prime}(t)=v(t,\xi(t)).$ (3.6)

Moreover, there exists a minimiser
$\displaystyle\rho^{*}\in\operatorname*{argmin}_{\rho\in\mathcal{M}^{+}(\left[0,1\right]\times\overline{\Omega})}\mathcal{E}(\rho)$
such that

$\text{for some }a_{i}\geq 0,\
\xi^{i}\in\operatorname{AC}^{2}(\left[0,1\right];\overline{\Omega}),\qquad\forall
t\in\left[0,1\right],\qquad\rho^{*}_{t}=\sum_{i=1}^{m(T+1)}a_{i}\delta_{\xi^{i}(t)}.$
(3.7)

######  Remark 3.1.

The formal definition of $\operatorname{AC}^{2}$ is given in [Ambrosio et al.,
2008, Sec. 1.1], but here it can be thought of as $H^{1}$, the space of
$L^{2}$ functions with weak derivative in $L^{2}$. See Lemma A.4 for more
details.

###### Proof of Theorem 3.1..

Fix $\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})$ with
$\mathcal{W}(\rho)<+\infty$. To apply [Thm. 7 Bredies et al., 2021], we need
to show that there exists $v\in
L^{2}_{\rho}({\left[0,1\right]\times\overline{\Omega}};\mathds{R}^{d})$ which
achieves the minimum of (3.1). To do this, note that the set

$\left\\{\,v\in
L^{2}_{\rho}({\left[0,1\right]\times\overline{\Omega}};\mathds{R}^{d})\mid\int_{\left[0,1\right]\times\overline{\Omega}}\left[\alpha+\frac{\beta}{2}|v|^{2}\right]\mathop{}\\!\mathrm{d}\rho\leq\mathcal{W}(\rho)+1,\
\partial_{t}\rho+\operatorname{div}(v\rho)=0\,\right\\}$ (3.8)

is bounded, therefore there is a weak limit-point $v\in
L^{2}_{\rho}({\left[0,1\right]\times\overline{\Omega}};\mathds{R}^{d})$ with
$\mathcal{W}(\rho)=\int[\alpha+\frac{\beta}{2}|v|^{2}]\mathop{}\\!\mathrm{d}\rho$.
The weak form of the continuity equation is preserved under weak limits in
$v$, therefore the triplet $(\rho,v,0)$ also satisfies (1.9). Then, [Bredies
et al., 2021, Thm. 7] (or Theorem 2.2(4)) confirms the existence of
$\sigma\in\mathcal{M}^{+}(\Gamma)$ with the desired properties. The second
statement, concerning the minimiser $\rho^{*}$, is exactly the result of
[Bredies et al., 2021, Thm. 10]. ∎

These results highlight a very close relationship between the representations
$\rho$ and $\sigma$, i.e. dynamical measures and measures on paths. In the
following subsections we first show that minimisation of (3.3) over
$\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})$ can be
expressed as an equivalent minimisation over
$\sigma\in\mathcal{M}^{+}(\Gamma)$. This reformulation will then allow us to
propose new and more efficient numerical schemes to compute reconstructions.

###  3.1 Reformulation in space of measures on paths

Here we propose a new energy
$\operatorname{E}\colon\mathcal{M}^{+}(\Gamma)\to\mathds{R}$ whose minimisers
coincide with those of
$\mathcal{E}\colon\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})\to\mathds{R}$.
The main problem is non-injectivity of $\Theta$, so we do not necessarily
expect to have $\operatorname{E}(\sigma)=\mathcal{E}(\Theta(\sigma))$ for all
$\sigma$. On the other hand, the value of the data-fidelity is unaffected by
the non-uniqueness. For any
$\rho\in\mathscr{C}_{w}(\left[0,1\right];\mathcal{M}(\overline{\Omega}))$ and
$\sigma\in\mathcal{M}^{+}(\Gamma)$ such that $\rho=\Theta(\sigma)$, we have

$\forall
t\in\left[0,1\right],\qquad\rho_{t}=\Theta(\sigma)_{t}=(e_{t})_{\sharp}{\sigma},$
(3.9)

(see (1.8) and Theorem 2.2). Motivated by this equivalence of “time slice”
measures, we propose the energy

$\displaystyle\forall\sigma\in\mathcal{M}^{+}(\Gamma),\qquad\operatorname{E}(\sigma)$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\frac{1}{2}\sum_{j=0}^{T}{\left\lVert
A_{j}(e_{t_{j}})_{\sharp}{\sigma}-b_{j}\right\rVert}_{2}^{2}+\operatorname{W}(\sigma),$
(3.10) $\displaystyle\operatorname{W}(\sigma)$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\int_{\Gamma}\int_{0}^{1}\left[\alpha+\frac{\beta}{2}|\xi^{\prime}|^{2}\right]\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}\sigma(1,\xi).$
(3.11)

Since the mass $h$ is constant in this example, we will write $\xi$ or
$(1,\xi)\in\Gamma$ interchangeably for the remainder of this section. Key
results from Bredies et al. [2021] show that
$\operatorname{W}\colon\mathcal{M}^{+}(\Gamma)\to\mathds{R}$ is equivalent to
$\mathcal{W}\colon\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})\to\mathds{R}$
in the following way.

######  Lemma 3.2.

The energy $\operatorname{E}$ is convex and lower semi-continuous. Also, for
any $\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})$,

$\mathcal{W}(\rho)=\min_{\sigma\in\mathcal{M}^{+}(\Gamma)}\left\\{\,\operatorname{W}(\sigma)\mid\Theta(\sigma)=\rho\,\right\\}.$
(3.12)

We conclude that

$\min\left\\{\,\operatorname{E}(\sigma)\mid\sigma\in\mathcal{M}^{+}(\Gamma)\,\right\\}=\min\left\\{\,\mathcal{E}(\rho)\mid\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})\,\right\\},$
(3.13)

and there is a minimising
$\sigma^{*}\in\operatorname*{argmin}_{\sigma\in\mathcal{M}^{+}(\Gamma)}\operatorname{E}(\sigma)$
of the form

$\text{for some }a_{i}\geq 0,\
(1,\xi^{i})\in\Gamma,\qquad\sigma^{*}=\sum_{i=1}^{m(T+1)}a_{i}\delta_{(1,\xi^{i})}.$
(3.14)

###### Proof.

As $\Gamma\subset\Gamma_{\infty}$, Lemma 2.3 confirms that
$A_{j}(e_{t_{j}})_{\sharp}{}\colon\mathcal{M}^{+}(\Gamma)\to\mathds{R}^{m}$ is
linear and narrowly continuous for each $j=0,\ldots,T$, so $\operatorname{E}$
is still convex and lower semi-continuous.

Secondly, fix
$\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})$. If
$\mathcal{W}(\rho)=+\infty$, then
$\mathcal{W}(\rho)\geq\min\left\\{\,\operatorname{W}(\sigma)\mid\Theta(\sigma)=\rho\,\right\\}$
is clear. Otherwise, let $\widehat{\sigma}\in\mathcal{M}^{+}(\Gamma)$ be the
measure given by Theorem 3.1 satisfying $\rho=\Theta(\widehat{\sigma})$, then

$\displaystyle\mathcal{W}(\rho)$
$\displaystyle=\int_{\left[0,1\right]\times\overline{\Omega}}\left[\alpha+\frac{\beta}{2}|v|^{2}\right]\mathop{}\\!\mathrm{d}\rho$
Theorem 3.1 (3.15)
$\displaystyle=\int_{\Gamma}\int_{0}^{1}\left[\alpha+\frac{\beta}{2}|v(t,\xi(t))|^{2}\right]\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}\widehat{\sigma}(\xi)$
Theorem 2.2(2) (3.16)
$\displaystyle=\int_{\Gamma}\int_{0}^{1}\left[\alpha+\frac{\beta}{2}|\xi^{\prime}(t)|^{2}\right]\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}\widehat{\sigma}(\xi)=\operatorname{W}(\widehat{\sigma}).$
$\displaystyle\eqref{eq: BB v = x'}$ (3.17)

This confirms
$\mathcal{W}(\rho)\geq\min\left\\{\,\operatorname{W}(\sigma)\mid\Theta(\sigma)=\rho\,\right\\}$
for all $\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})$.

The converse holds by Jensen’s inequality. In particular, because
$\mathcal{W}$ is convex, proper, and lower semi-continuous ([Bredies et al.,
2021, Lemma 2]), by [Thm. 5 Rockafellar, 1974] there exists a collection
$\\{(c_{i},\psi_{i})\\}_{i\in
I}\subset\mathds{R}\times\mathscr{C}({\left[0,1\right]\times\overline{\Omega}})$
such that

$\forall\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}}),\qquad\mathcal{W}(\rho)=\sup_{i\in
I}\left(c_{i}+\int_{\left[0,1\right]\times\overline{\Omega}}\psi_{i}\mathop{}\\!\mathrm{d}\rho\right).$
(3.18)

Because $\mathcal{W}$ is positively homogeneous, this simplifies with
$c_{i}=0$. For any $\widehat{\sigma}\in\mathcal{M}^{+}(\Gamma)$, $i\in I$, we
have
$\Theta(\widehat{\sigma})\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})$
and

$\int_{\left[0,1\right]\times\overline{\Omega}}\psi_{i}\mathop{}\\!\mathrm{d}\Theta(\widehat{\sigma})=\int_{\Gamma}\int_{0}^{1}\psi_{i}(t,\xi(t))\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}\widehat{\sigma}\leq\int_{\Gamma}\sup_{j\in
I}\int_{0}^{1}\psi_{j}(t,\xi(t))\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}\widehat{\sigma}\\\
=\int_{\Gamma}\left[\sup_{j\in
I}\int_{\left[0,1\right]\times\overline{\Omega}}\psi_{j}\mathop{}\\!\mathrm{d}\Theta(\delta_{\xi})\right]\mathop{}\\!\mathrm{d}\widehat{\sigma}=\int_{\Gamma}\mathcal{W}(\Theta(\delta_{\xi}))\mathop{}\\!\mathrm{d}\widehat{\sigma}.$
(3.19)

It is shown in [Bredies et al., 2021, (19)] that
$\mathcal{W}(\Theta(\delta_{\xi}))=\int_{0}^{1}\alpha+\frac{\beta}{2}|\xi^{\prime}(t)|^{2}\mathop{}\\!\mathrm{d}t$.
Now taking a supremum over $i\in I$ gives

$\mathcal{W}(\Theta(\widehat{\sigma}))\leq\int_{\Gamma}\mathcal{W}(\Theta(\delta_{\xi}))\mathop{}\\!\mathrm{d}\widehat{\sigma}=\int_{\Gamma}\int_{0}^{1}\alpha+\frac{\beta}{2}|\xi^{\prime}(t)|^{2}\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}\widehat{\sigma}=\operatorname{W}(\widehat{\sigma}).$
(3.20)

This shows
$\mathcal{W}(\rho)\leq\min\left\\{\,\operatorname{W}(\sigma)\mid\Theta(\sigma)=\rho\,\right\\}$
for all $\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})$,
which becomes equality when combined with the first result.

Finally we consider the equivalence of minimums and minimisers. With $a_{i}$
and $\xi^{i}$ given by $\rho^{*}$ in (3.7), define $\sigma^{*}$ as in (3.14).
One can check that $\Theta(\sigma^{*})=\rho^{*}$. Observe each measure has
equal energy:

$\displaystyle\mathcal{W}(\rho^{*})$
$\displaystyle=\sum_{i=1}^{m(T+1)}a_{i}\left[\int_{0}^{1}\alpha+\frac{\beta}{2}|(\xi^{i})^{\prime}(t)|^{2}\mathop{}\\!\mathrm{d}t\right]$
[Bredies et al., 2021, Thm. 10] (3.21)
$\displaystyle=\operatorname{W}\left(\sum_{i=1}^{m(T+1)}a_{i}\delta_{\xi^{i}}\right)=\operatorname{W}(\sigma^{*}),$
(linearity of $\operatorname{W}$) (3.22)

therefore

$\min_{\sigma\in\mathcal{M}^{+}(\Gamma)}\operatorname{E}(\sigma)\leq\frac{1}{2}\sum_{j=0}^{T}{\left\lVert
A_{j}(e_{t_{j}})_{\sharp}{\sigma^{*}}-b_{j}\right\rVert}_{2}^{2}+\operatorname{W}(\sigma^{*})=\mathcal{E}(\rho^{*})=\min_{\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})}\mathcal{E}(\rho).$
(3.23)

Conversely, for any $\sigma\in\mathcal{M}^{+}(\Gamma)$ we have

$\min_{\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})}\mathcal{E}(\rho)\leq\mathcal{E}(\Theta(\sigma))\leq\frac{1}{2}\sum_{j=0}^{T}{\left\lVert
A_{j}(e_{t_{j}})_{\sharp}{\sigma}-b_{j}\right\rVert}_{2}^{2}+\operatorname{W}(\sigma)=\operatorname{E}(\sigma).$
(3.24)

We conclude that
$\min_{\sigma\in\mathcal{M}^{+}(\Gamma)}\operatorname{E}(\sigma)=\operatorname{E}(\sigma^{*})=\mathcal{E}(\rho^{*})=\min_{\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})}\mathcal{E}(\rho)$
as required. ∎

This theorem justifies that $\Theta$ can be thought of as a reparametrisation;
the two optimisation problems are equivalent. Although different measures
$\sigma$ might represent the same mass $\rho$ (see for instance [Bredies et
al., 2020b, Rem. 4.2]), the point of view of measures on paths $\sigma$ is
still interesting as it provides plausible trajectories for physical
particles.

###  3.2 Minimisation algorithm

A major advantage of Frank-Wolfe / Generalised Conditional Gradient methods is
their ability to minimise energies like (3.3) while constructing iterates
which share the same sparse structure (3.7) as the minimisers (see for
instance Bredies and Pikkarainen [2013], Boyd et al. [2017], Denoyelle et al.
[2019] in the case of point sources recovery). The method proposed in Bredies
et al. [2020b] follows that approach, and the core of the algorithm is the
minimisation of a linear functional on the level sets of $\mathcal{W}$ (see
[Bredies et al., 2020b, Prop. 3.6]). The description of the extreme points of
those level sets in Bredies et al. [2021] makes it possible to reformulate
that problem as a variational problem on curves.

In the present section, we point out that, while the algorithm derived in
Bredies et al. [2020b] is quite natural, there are many possible ways to
define a Frank-Wolfe–based algorithm minimising (3.3).

Relying on the linear map $\Theta$, we note that the minimisation of (3.10) is
equivalent to

$\displaystyle\min_{\sigma\in\mathcal{M}^{+}(\Gamma)}\frac{1}{2}\sum_{j=0}^{T}{\left\lVert
A_{j}(e_{t_{j}})_{\sharp}{\sigma}-b_{j}\right\rVert}_{2}^{2}+\int_{\Gamma}w\mathop{}\\!\mathrm{d}\sigma$
(3.25) $\displaystyle\mbox{where}\quad\forall\xi\in\Gamma,\quad
w(\xi)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\alpha+\frac{\beta}{2}\int_{0}^{1}|\xi^{\prime}(t)|^{2}\mathop{}\\!\mathrm{d}t.$
(3.26)

The standard Frank-Wolfe procedure consists in iteratively minimising the
linearisation of (3.25) on a compact set (for some suitable topology) to yield
descent directions. At iteration $n$, let
$\eta_{j}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}A_{j}^{*}\left(A_{j}(e_{t_{j}})_{\sharp}{\sigma^{n}}-b_{j}\right)$
denote the linearisation of each summand of the fidelity term. We are led to
minimise

$\min_{\sigma\in
C}\int_{\Gamma}\left[\sum_{j=0}^{T}\eta_{j}(\xi(t_{j}))+w(\xi)\right]\mathop{}\\!\mathrm{d}\sigma,$
(3.27)

where $C\subset\mathcal{M}^{+}(\Gamma)$ is a suitably chosen compact convex
set which contains the minimisers of $\operatorname{E}$ on
$\mathcal{M}^{+}(\Gamma)$. A first choice is to consider the set

$C_{1}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left\\{\,\sigma\in\mathcal{M}^{+}(\Gamma)\mid\int_{\Gamma}w\mathop{}\\!\mathrm{d}\sigma\leq\operatorname{E}(0)\,\right\\},$
(3.28)

which does contain the minimisers of (3.25) since they must satisfy
$\int_{\Gamma}w\mathop{}\\!\mathrm{d}\sigma\leq\operatorname{E}(\sigma)\leq\operatorname{E}(0)$.
By Lemma A.5, the set of extreme points of $C_{1}$ is
$\\{0\\}\cup\left\\{\,\frac{\operatorname{E}(0)}{w(\xi)}\delta_{\xi}\mid\xi\in\Gamma\,\right\\}$,
which yields the descent direction

$\xi^{n+1}\in\operatorname*{argmin}_{\xi\in\Gamma}\frac{\sum_{j=0}^{T}\eta_{j}(\xi_{j})+w(\xi)}{w(\xi)}=\operatorname*{argmin}_{\xi\in\Gamma}\frac{\sum_{j=0}^{T}\eta_{j}(\xi_{j})}{w(\xi)}.$
(3.29)

That is precisely the descent direction used in [Bredies et al., 2020b,
(4.30)], although in that context (3.27) is written as a minimisation over
$\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})$.

However, alternative choices for $C$ are possible. We propose to use

$C_{2}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left\\{\,\sigma\in\mathcal{M}^{+}(\Gamma)\mid\int_{\Gamma}\alpha\mathop{}\\!\mathrm{d}\sigma\leq\operatorname{E}(0)\,\right\\},$
(3.30)

which contains the minimisers of (3.25) too, since they must satisfy
${\int_{\Gamma}\alpha\mathop{}\\!\mathrm{d}\sigma\leq\int_{\Gamma}w\mathop{}\\!\mathrm{d}\sigma\leq\operatorname{E}(\sigma)\leq\operatorname{E}(0)}$.
Applying Lemma A.5, we see that the set of extreme points of $C_{2}$ is
$\\{0\\}\cup\left\\{\,\frac{\operatorname{E}(0)}{\alpha}\delta_{\xi}\mid\xi\in\Gamma\,\right\\}$,
which yields the linear minimisation step

$\xi^{n+1}\in\operatorname*{argmin}_{\xi\in\Gamma}\frac{\sum_{j=0}^{T}\eta_{j}(\xi_{j})+w(\xi)}{1}=\operatorname*{argmin}_{\xi\in\Gamma}\sum_{j=0}^{T}\eta_{j}(\xi_{j})+w(\xi).$
(3.31)

We found (3.31) more convenient than (3.29) as it is amenable to dynamic
programming techniques (see Section 3.4). The whole minimisation algorithm is
summarised in Algorithm 1.

######  Remark 3.2.

In measure spaces, compactness is equivalent to boundedness and tightness. In
both cases, boundedness comes from $\inf_{\gamma\in\Gamma}\varphi>0$, and
tightness comes from the fact that $w$ has compact sub-levelsets. The set
$C_{1}$ is compact by Lemma A.3. On the other hand, $C_{2}$ is only bounded,
but sub-levelsets of (3.27) are still compact, as argued in Lemma 5.4.

1:Choose $\sigma^{0}=0$, $n\leftarrow 0$

2:repeat

3: $\eta_{j}\leftarrow
A_{j}^{*}(A_{j}(e_{t_{j}})_{\sharp}{\sigma^{n}}-b_{j})\in\mathscr{C}(\overline{\Omega})$
$\triangleright$ Gradient of data term

4: Compute $\xi^{n+1}$ according to (3.29) or (3.31)$\triangleright$ Linear
oracle step

5: Choose $\mu^{n+1}\propto\delta_{\xi^{n+1}}$ such that
$\mu^{n+1}\in\operatorname{Ext}\left(C\right)$ $\triangleright$ Choice of
scaling

6: $\sigma^{n+1}\leftarrow(1-\lambda_{n})\sigma^{n}+\lambda_{n}\mu^{n+1}$
$\triangleright$ Some choice of $\lambda_{n}\in\left[0,1\right]$

7: $n\leftarrow n+1$

8:until converged

Algorithm 1 Frank-Wolfe algorithm for the Benamou-Brenier example

###  3.3 Discrete-time formulation

In the previous two sections (Lemma 3.2 and Section 3.2) we have shown that
minimisers $\sigma^{*}$ of (3.10) can be approximated whilst only considering
spatially discrete measures. Because of the discrete-time data, we can further
simplify the numerical computations by considering discrete-time measures.
Observe

$\displaystyle\operatorname{E}(\sigma^{*})$
$\displaystyle=\min_{\sigma\in\mathcal{M}^{+}(\Gamma)}\Bigg{\\{}\frac{1}{2}\sum_{j=0}^{T}{\left\lVert
A_{j}(e_{t_{j}})_{\sharp}{\sigma}-b_{j}\right\rVert}_{2}^{2}+\operatorname{W}(\sigma)\Bigg{\\}}$
(3.32)
$\displaystyle=\min_{\begin{subarray}{c}\rho_{j}\in\mathcal{M}^{+}(\overline{\Omega})\\\
j=0,\ldots,T\end{subarray}}\
\min_{\begin{subarray}{c}\sigma\in\mathcal{M}^{+}(\Gamma)\\\
(e_{t_{j}})_{\sharp}{\sigma}=\rho_{j}\end{subarray}}\Bigg{\\{}\frac{1}{2}\sum_{j=0}^{T}{\left\lVert
A_{j}\rho_{j}-b_{j}\right\rVert}_{2}^{2}+\operatorname{W}(\sigma)\Bigg{\\}}$
(3.33)
$\displaystyle=\min_{\begin{subarray}{c}\rho_{j}\in\mathcal{M}^{+}(\overline{\Omega})\\\
j=0,\ldots,T\end{subarray}}\Bigg{\\{}\frac{1}{2}\sum_{j=0}^{T}{\left\lVert
A_{j}\rho_{j}-b_{j}\right\rVert}_{2}^{2}+\min_{\begin{subarray}{c}\sigma\in\mathcal{M}^{+}(\Gamma)\\\
(e_{t_{j}})_{\sharp}{\sigma}=\rho_{j}\end{subarray}}\operatorname{W}(\sigma)\Bigg{\\}}.$
(3.34)

The minimal $\rho_{j}^{*}$ will inherit the same sparse structure from
$\rho^{*}$. This highlights that the minimal energy can be computed exactly
whilst only considering measures which are discrete in both time and space,
then the continuous-time minimiser is recovered as a
$\operatorname{W}$-minimising interpolation. In the $\sigma$ measures-on-paths
representation, this continuous-time interpolation can be written path-wise.
We simply replace any discrete path $(\xi^{i}(t_{0}),\ldots,\xi^{i}(t_{T}))$
with the Benamou-Brenier geodesic curve interpolating those points.

######  Lemma 3.3.

Assume that $\Omega$ is convex and
$\sigma^{*}\in\operatorname*{argmin}_{\sigma\in\mathcal{M}^{+}(\Gamma)}\operatorname{E}(\sigma)$,
then

$\text{for }\sigma^{*}\text{-a.e.
}\xi,\qquad\xi\in\operatorname*{argmin}_{\widetilde{\xi}\in\Gamma}\left\\{\int_{0}^{1}\left[\alpha+\frac{\beta}{2}|\widetilde{\xi}^{\prime}|^{2}\right]\mathop{}\\!\mathrm{d}t\quad\operatorname{\;s.t.\;}\quad\widetilde{\xi}(t_{j})=\xi(t_{j})\
\text{for all }j\right\\}.$ (3.35)

This is an example of Lemma B.1. In this specific case, as commented in
[Bredies et al., 2020b, Rem. 4.10], paths are piecewise linear. If we restrict
$\Gamma$ to this set of geodesics

$\Gamma\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left\\{\,(h,\xi)\in\\{1\\}\times\operatorname{AC}^{2}(\left[0,1\right];\overline{\Omega})\mid\mbox{$\xi^{\prime}$
is constant on $\left]t_{j-1},t_{j}\right[$, $1\leq j\leq T$}\,\right\\},$
(3.36)

then we can write $\operatorname{W}$ in its explicit discrete-time form:

$\operatorname{W}(\sigma)=\int_{\Gamma}\left[\alpha+\frac{\beta}{2}\sum_{j=1}^{T}\frac{|\xi(t_{j})-\xi(t_{j-1})|^{2}}{t_{j}-t_{j-1}}\right]\mathop{}\\!\mathrm{d}\sigma(\xi).$
(3.37)

###  3.4 Linear oracle computation

The most challenging step of Algorithm 1 for our choice of $C$ is the
minimisation (3.31), the linear oracle computation. Fix
$\sigma^{n}\in\mathcal{M}^{+}(\Gamma)$. Utilising the discrete-time form of
$\operatorname{W}$, (3.31) becomes

$\displaystyle\xi^{n+1}$
$\displaystyle\in\operatorname*{argmin}_{\xi\in\overline{\Omega}^{T+1}}\widetilde{\operatorname{E}}_{T}(\xi),$
(3.38) $\displaystyle\text{where }\forall J=0,\ldots,T,\
\xi\in\overline{\Omega}^{T+1},$
$\displaystyle\widetilde{\operatorname{E}}_{J}(\xi)$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\sum_{j=0}^{J}\eta_{j}(\xi_{j})+\alpha+\frac{\beta}{2}\sum_{j=1}^{J}\frac{|\xi_{j}-\xi_{j-1}|^{2}}{t_{j}-t_{j-1}}$
(3.39) $\displaystyle\text{and }\forall j=0,\ldots,T,$ $\displaystyle\eta_{j}$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}A_{j}^{*}(A_{j}(e_{t_{j}})_{\sharp}{\sigma^{n}}-b_{j}).$ (3.40)

Unless $\eta_{j}$ has some special structure, computing $\xi^{n+1}$ is a
challenging non-convex optimisation. As $\Omega$ is a non-discrete set, it
cannot be computed exactly in finite time. Even if $\Omega$ is discretised,
the search space is still of size $|\Omega|^{T+1}$. In Bredies et al. [2020b]
they use a stochastic search of $\overline{\Omega}^{T+1}$ to approximate
$\xi^{n+1}$. We propose taking a discrete subset, and then use the discrete
minimiser for $\xi^{n+1}$.

Assume that we are given a family of finite sets
$\Lambda_{j}\subseteq\overline{\Omega}$ (“finite grids” which discretise
$\overline{\Omega}$), for $j=0,\ldots,T$. For each $J=0,\ldots,T$ and
$y\in\Lambda_{J}$ define the optimal curves

$Y^{*}[y,J]\in\operatorname*{argmin}_{Y\in\overline{\Omega}^{J+1}}\left\\{\,\widetilde{\operatorname{E}}_{J}(Y)\mid
Y_{J}=y,\ Y_{j}\in\Lambda_{j}\text{ for }j=0,\ldots,J-1\,\right\\}.$ (3.41)

Observe that for all $J=1,\ldots,T$ and $y\in\Lambda_{J}$,

$\displaystyle\min_{\begin{subarray}{c}Y\in\prod_{j=0}^{J}\Lambda_{j}\\\
Y_{J}=y\end{subarray}}\widetilde{\operatorname{E}}_{J}(Y)$
$\displaystyle=\min_{\begin{subarray}{c}Y\in\prod_{j=0}^{J}\Lambda_{j}\end{subarray}}\left\\{\widetilde{\operatorname{E}}_{J-1}(Y)+\eta_{J}(y)+\alpha+\frac{\beta}{2}\frac{|y-Y_{J-1}|^{2}}{t_{J}-t_{J-1}}\right\\}$
(3.42)
$\displaystyle=\eta_{J}(y)+\alpha+\min_{\begin{subarray}{c}y^{\prime}\in\Lambda_{J-1}\end{subarray}}\Bigg{\\{}\frac{\beta}{2}\frac{|y-y^{\prime}|^{2}}{t_{J}-t_{J-1}}+\min_{\begin{subarray}{c}Y\in\prod_{j=0}^{J-1}\Lambda_{j}\\\
Y_{J-1}=y^{\prime}\end{subarray}}\widetilde{\operatorname{E}}_{J-1}(Y)\Bigg{\\}}.$
(3.43)

In particular, we can choose $Y^{*}[y,J]$ inductively to be

$\displaystyle
Y^{*}[y,J]=\begin{pmatrix}Y^{*}[y^{\prime},J-1]&y\end{pmatrix}\quad\text{for}\quad
y^{\prime}\in\operatorname*{argmin}_{y^{\prime}\in\Lambda_{J-1}}\left[\frac{\beta}{2}\frac{|y-y^{\prime}|^{2}}{t_{J}-t_{J-1}}+\widetilde{\operatorname{E}}_{J-1}(Y^{*}[y^{\prime},J-1])\right].$
(3.44)

The key complexity result is that we can compute this minimiser in linear time
in $T$, rather than the naive exponential scaling.

##  4 The unbalanced Wasserstein-Fisher-Rao example

In the last section we showed that optimisation of the Benamou-Brenier energy
can be reformulated as an optimisation over the space of measures on paths. In
this section we will justify that this is also true for the Wasserstein-
Fisher-Rao (WFR) example considered in Bredies et al. [2020a], where the
transportation cost allows for changes of mass. The two arguments are very
similar and so we will only emphasise the differences and reference the
different published results.

The regularisation penalty considered in [Bredies et al., 2020a, (33) and
(152)] is
$\mathcal{W}\colon\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})\to\mathds{R}\cup\\{+\infty\\}$
such that for all
$\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})$,

$\mathcal{W}(\rho)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\inf_{\begin{subarray}{c}v\in
L^{2}_{\rho}(\left[0,1\right]\times\overline{\Omega};\mathds{R}^{d})\\\ g\in
L^{2}_{\rho}(\left[0,1\right]\times\overline{\Omega})\end{subarray}}\Bigg{\\{}\int_{\left[0,1\right]\times\overline{\Omega}}\left[\alpha+\frac{\beta}{2}|v|^{2}+\frac{\beta\delta^{2}}{2}g^{2}\right]\mathop{}\\!\mathrm{d}\rho\
\operatorname{\;s.t.\;}\
\partial_{t}\rho+\operatorname{div}(v\rho)=g\rho\Bigg{\\}}$ (4.1)

for some $\alpha,\beta,\delta>0$, and the continuity equation is satisfied in
the sense of (1.9). As in (3.3), we choose the energy

$\mathcal{E}(\rho)=\frac{1}{2}\sum_{j=0}^{T}{\left\lVert
A_{j}\rho_{t_{j}}-b_{j}\right\rVert}_{2}^{2}+\mathcal{W}(\rho).$ (4.2)

Much is already known about the structure of minimisers of this energy.

######  Theorem 4.1 ([Bredies et al., 2020a, Thms. 4.2,6.4]).

Let $\alpha,\beta,\delta>0$ and

$\Gamma\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left\\{\,(h,\xi)\colon\left[0,1\right]\to\left[0,1\right]\times\overline{\Omega}\mid\sqrt{h}\in\operatorname{AC}^{2}(\left[0,1\right]),\
\sqrt{h}\xi\in\operatorname{AC}^{2}(\left[0,1\right];\mathds{R}^{d})\,\right\\}.$
(4.3)

Then for all
$\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})$ with
$\mathcal{W}(\rho)<+\infty$, there exists $\sigma\in\mathcal{M}^{+}(\Gamma)$,
$v\in L^{2}_{\rho}({\left[0,1\right]\times\overline{\Omega}};\mathds{R}^{d})$,
$g\in L^{2}_{\rho}({\left[0,1\right]\times\overline{\Omega}})$ such that
$\rho=\Theta(\sigma)$,

$\mathcal{W}(\rho)=\int_{\left[0,1\right]\times\overline{\Omega}}\left[\alpha+\frac{\beta}{2}|v|^{2}+\frac{\beta\delta^{2}}{2}g^{2}\right]\mathop{}\\!\mathrm{d}\rho,\qquad\partial_{t}\rho+\operatorname{div}(v\rho)=g\rho,$
(4.4)

and

$\displaystyle\xi^{\prime}(t)$
$\displaystyle=v(t,\xi(t))\quad\text{for}\quad\text{a.e. }t\in\\{h>0\\}\text{
and }\sigma\text{-a.e. }(h,\xi),$ (4.5) $\displaystyle h^{\prime}(t)$
$\displaystyle=g(t,\xi(t))h(t)\quad\text{for}\quad\text{a.e.
}t\in\left[0,1\right]\text{ and }\sigma\text{-a.e. }(h,\xi).$

Moreover, there exists a minimiser
$\displaystyle\rho^{*}\in\operatorname*{argmin}_{\rho\in\mathcal{M}^{+}(\left[0,1\right]\times\overline{\Omega})}\mathcal{E}(\rho)$
such that

$\text{for some }a_{i}\geq 0,\ (h^{i},\xi^{i})\in\Gamma,\qquad\forall
t\in\left[0,1\right],\qquad\rho^{*}_{t}=\sum_{i=1}^{m(T+1)}a_{i}h^{i}(t)\delta_{\xi^{i}(t)}.$
(4.6)

###### Sketch proof..

If $\mathcal{W}(\rho)<+\infty$, the existence of optimal $v$ and $g$ for (4.1)
is again given by weak compactness in
$L^{2}_{\rho}({\left[0,1\right]\times\overline{\Omega}})$. The result of
[Bredies et al., 2020a, Thm. 4.2] (or Theorem 2.2) shows that
$\rho=\Theta(\sigma)$ for some $\sigma\in\mathcal{M}^{+}(\Gamma)$ satisfying
(4.5). Finally, the structure of $\rho^{*}$ is given by [Bredies et al.,
2020a, Thm. 6.4]. ∎

###  4.1 Reformulation in space of measures on paths

The function $\mathcal{E}$ can be reformulated into an equivalent energy
$\operatorname{E}\colon\mathcal{M}^{+}(\Gamma)\to\mathds{R}\cup\\{+\infty\\}$
defined by

$\displaystyle\forall\sigma\in\mathcal{M}^{+}(\Gamma),\qquad\operatorname{E}(\sigma)$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\frac{1}{2}\sum_{j=0}^{T}{\left\lVert
A_{j}(e_{t_{j}})_{\sharp}{\sigma}-b_{j}\right\rVert}_{2}^{2}+\operatorname{W}(\sigma),$
(4.7) $\displaystyle\operatorname{W}(\sigma)$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\int_{\Gamma}\int_{0}^{1}\left[\alpha+\frac{\beta}{2}|\xi^{\prime}|^{2}+\frac{\beta\delta^{2}}{2}\left(\frac{h^{\prime}}{h}\right)^{2}\right]h\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}\sigma(h,\xi).$
(4.8)

######  Remark 4.1.

Since for all $(h,\xi)\in\Gamma$, $\sqrt{h}$ and $\sqrt{h}\xi$ are absolutely
continuous, they are differentiable almost everywhere in $\left[0,1\right]$,
hence $\xi$ is differentiable at a.e. $t$ such that $h(t)>0$. Hence, the
integrand in (4.8) makes sense when regarded as

$\displaystyle\int_{\\{h>0\\}}\left[\alpha+\frac{\beta}{2}{\left|\left(\frac{\sqrt{h}\xi}{\sqrt{h}}\right)^{\prime}(t)\right|}^{2}+{2\beta\delta^{2}}{\left|(\sqrt{h})^{\prime}(t)\right|}^{2}\right]h\mathop{}\\!\mathrm{d}t.$
(4.9)

We can now replicate the extension performed in Lemma 3.2

######  Lemma 4.2.

The energy $\operatorname{E}$ is convex and lower semi-continuous. Also, for
any $\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})$,

$\mathcal{W}(\rho)=\min_{\sigma\in\mathcal{M}^{+}(\Gamma)}\left\\{\,\operatorname{W}(\sigma)\mid\Theta(\sigma)=\rho\,\right\\}.$
(4.10)

We conclude that

$\min\left\\{\,\operatorname{E}(\sigma)\mid\sigma\in\mathcal{M}^{+}(\Gamma)\,\right\\}=\min\left\\{\,\mathcal{E}(\rho)\mid\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})\,\right\\},$
(4.11)

and there is a minimising
$\sigma^{*}\in\operatorname*{argmin}_{\sigma\in\mathcal{M}^{+}(\Gamma)}\operatorname{E}(\sigma)$
of the form

$\text{for some }a_{i}\geq 0,\
(h^{i},\xi^{i})\in\Gamma,\qquad\sigma^{*}=\sum_{i=1}^{m(T+1)}a_{i}\delta_{(h^{i},\xi^{i})}.$
(4.12)

###### Sketch proof..

As before, Lemma 2.3 confirms that $\operatorname{E}$ is convex and lower
semi-continuous and the existence results of Theorem 4.1 show that
$\mathcal{W}(\rho)\geq\min_{\sigma\in\mathcal{M}^{+}(\Gamma)}\left\\{\,\operatorname{W}(\sigma)\mid\Theta(\sigma)=\rho\,\right\\}$
for all $\rho\in\mathcal{M}^{+}({\left[0,1\right]\times\overline{\Omega}})$.

Conversely, [Bredies et al., 2020a, Lemma A.6] shows that $\mathcal{W}$ is
convex, proper, lower semi-continuous, and positively homogeneous. Therefore
the argument of Lemma 3.2 by Jensen’s inequality holds, giving us
$\mathcal{W}(\rho)\leq\int_{\Gamma}\mathcal{W}(\delta_{\gamma})\mathop{}\\!\mathrm{d}\sigma$
for any $\sigma\in\mathcal{M}^{+}(\Gamma)$ such that $\Theta(\sigma)=\rho$. It
is the result of [Bredies et al., 2020a, Prop. 3.9] that confirms
$\int_{\Gamma}\mathcal{W}(\delta_{\gamma})\mathop{}\\!\mathrm{d}\sigma=\operatorname{W}(\sigma)$.

Finally we consider the equivalence of minimums and minimisers. Define
$\rho^{*}$ as in (4.6) and $\sigma^{*}$ as in (4.12), again we can confirm
$\Theta(\sigma^{*})=\rho^{*}$. The required equality
$\mathcal{W}(\rho^{*})=\sum_{i=1}^{m(T+1)}a_{i}\operatorname{W}(\delta_{\gamma^{i}})$
is given by [Bredies et al., 2020a, Thm. 6.4]. The rest of the proof follows
as before. ∎

###  4.2 Discrete-time formulation

Problems of the form (4.7) also have a discrete-time structure inherited from
discrete-time data. The same argument from Section 3.3 (namely the use of
(3.34) and Lemma B.1) is valid for unbalanced examples. In the context of the
WFR penalty, this is stated as follows.

######  Lemma 4.3.

Assume that $\Omega$ is convex and
$\sigma^{*}\in\operatorname*{argmin}_{\sigma\in\mathcal{M}^{+}(\Gamma)}\operatorname{E}(\sigma)$,
define

$d_{\alpha,\beta,\delta}(\gamma_{0},t_{0},\gamma_{1},t_{1})=\inf_{\gamma\in\Gamma}\Bigg{\\{}\int_{t_{0}}^{t_{1}}\left[\alpha+\frac{\beta}{2}|\xi^{\prime}(t)|^{2}+\frac{\beta\delta^{2}}{2}\left(\frac{h^{\prime}(t)}{h(t)}\right)^{2}\right]h(t)\mathop{}\\!\mathrm{d}t\\\
\text{such that }\gamma=(h,\xi),\ \gamma(t_{0})=\gamma_{0},\
\gamma(t_{1})=\gamma_{1}\Bigg{\\}}.$ (4.13)

Then $\sigma^{*}$-a.e. $\gamma=(h,\xi)$ achieves minimal energy. In particular

$\int_{0}^{1}\left[\alpha+\frac{\beta}{2}|\xi^{\prime}|^{2}+\frac{\beta\delta^{2}}{2}\left(\frac{h^{\prime}}{h}\right)^{2}\right]h\mathop{}\\!\mathrm{d}t=\sum_{j=1}^{T}d_{\alpha,\beta,\delta}(\gamma(t_{j-1}),t_{j-1},\gamma(t_{j}),t_{j}).$
(4.14)

The analytical form of these WFR-geodesics is not so convenient as for the
Benamou-Brenier penalty, although for $\alpha=0$ it is known ([Chizat et al.,
2018a, Thm. 5.6])

$d_{0,\beta,\delta}((h_{0},\xi_{0}),t_{0},(h_{1},\xi_{1}),t_{1})=\frac{4\beta\delta^{2}}{t_{1}-t_{0}}\left[\tfrac{h_{1}+h_{0}}{2}-\sqrt{h_{1}h_{0}}\cos\left(\min\left(\tfrac{|\xi_{1}-\xi_{0}|}{2\delta},\pi\right)\right)\right].$
(4.15)

###  4.3 Linear oracle computation

Similarly to Section 3.4, the linear oracle computation is to find

$\gamma^{n+1}\in\operatorname*{argmin}_{\gamma\in(\left[0,1\right]\times\overline{\Omega})^{T+1}}\widetilde{\operatorname{E}}_{T}(\gamma),$
(4.16)

where, for all $\gamma=(h,\xi)$, $j,J=0,\ldots,T$ we define

$\widetilde{\operatorname{E}}_{J}(\gamma)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\sum_{j=0}^{J}h_{j}\eta_{j}(\xi_{j})+\sum_{j=1}^{J}d_{\alpha,\beta,\delta}(\gamma_{j-1},t_{j-1},\gamma_{j},t_{j})\quad\text{and}\quad\eta_{j}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}A_{j}^{*}(A_{j}(e_{t_{j}})_{\sharp}{\sigma^{n}}-b_{j}).$ (4.17)

This formulation is much less explicit than in Section 3.4, but in Section 5.3
we will see that this is sufficient to be covered by our general framework. In
particular, the numerical method explained in Section 7 is also applicable to
the WFR example.

##  5 Proposed framework and analytical results

In this section we will introduce the framework which we use to generalise the
examples of Sections 3 and 4, and state the main results which will be proven
in later sections.

###  5.1 Variational formulation

In this work we consider dynamical inverse problems with a discrete-time
structure. In particular, there exist observation times
$t_{j}\in\left[0,1\right]$, $j=0,\ldots,T$ and narrowly continuous linear
operators $A_{j}\colon\mathcal{M}(\overline{\Omega})\to\mathds{R}^{m}$. The
operators $A_{j}$ are described by
$a_{i}^{j}\in\mathscr{C}(\overline{\Omega})$ such that

$\forall\rho\in\mathcal{M}(\overline{\Omega}),\ i=1,\ldots,m,\
j=0,\ldots,T,\qquad(A_{j}\rho)_{i}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\int_{\overline{\Omega}}a_{i}^{j}(x)\mathop{}\\!\mathrm{d}\rho(x).$
(5.1)

As stated at the end of Section 2, we work with a closed set of curves

$\Gamma\subset\Gamma_{\infty}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left\\{\,\gamma=(h,\xi)\mid
h\in\mathscr{C}(\left[0,1\right];\left[-1,1\right]),\
\xi\colon\left[0,1\right]\to\overline{\Omega},\ \xi|_{\\{h\neq 0\\}}\text{ is
continuous}\,\right\\},$ (5.2)

so that $\forall t\in\left[0,1\right]$, the map
$(e_{t})_{\sharp}{}\colon\mathcal{M}(\Gamma)\to\mathcal{M}(\overline{\Omega})$
is narrowly continuous where

$\forall\sigma\in\mathcal{M}(\Gamma),\
\psi\in\mathscr{C}(\overline{\Omega}),\qquad\int_{\overline{\Omega}}\psi(x)\mathop{}\\!\mathrm{d}[(e_{t})_{\sharp}{\sigma}](x)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\int_{\Gamma}h(t)\psi(\xi(t))\mathop{}\\!\mathrm{d}\sigma(h,\xi).$
(5.3)

We can therefore define
$\operatorname{F}\colon\mathcal{M}(\Gamma)\to\mathds{R}$ by

$\text{for some convex }\operatorname{F}_{j}\in
C^{1}(\mathds{R}^{m};\mathds{R}),\qquad\operatorname{F}(\sigma)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\sum_{j=0}^{T}\operatorname{F}_{j}(A_{j}\left[(e_{t_{j}})_{\sharp}{\sigma}\right]).$
(5.4)

For lower semi-continuous $w,\varphi\colon\Gamma\to\left[0,+\infty\right]$
define
$\operatorname{W}\colon\mathcal{M}^{+}(\Gamma)\to\left]-\infty,+\infty\right]$,
$D\subset\mathcal{M}^{+}(\Gamma)$ by

$\displaystyle\forall\sigma\in D,$ $\displaystyle\operatorname{W}(\sigma)$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\int_{\Gamma}w(\gamma)\mathop{}\\!\mathrm{d}\sigma(\gamma)$ (5.5)
where $\displaystyle D$ $\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left\\{\,\sigma\in\mathcal{M}(\Gamma)\quad\mid\quad\sigma\geq
0,\quad\int_{\Gamma}\varphi(\gamma)\mathop{}\\!\mathrm{d}\sigma\leq
1\,\right\\}.$ (5.6)

Our general energy is of the form $\operatorname{E}\colon
D\to\left]-\infty,+\infty\right]$ defined by

$\forall\sigma\in
D,\qquad\operatorname{E}(\sigma)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\operatorname{F}(\sigma)+\operatorname{W}(\sigma).$ (5.7)

Similarly to the examples in Sections 3 and 4, we could use the map
$\Theta\colon\mathcal{M}(\Gamma)\to\mathscr{C}_{w}(\left[0,1\right];\mathcal{M}(\overline{\Omega}))$
to represent an energy
$\mathcal{E}\colon\mathscr{C}_{w}(\left[0,1\right];\mathcal{M}(\overline{\Omega}))\to\left]-\infty,+\infty\right]$
via the formula

$\mathcal{E}(\rho)=\sum_{j=0}^{T}\operatorname{F}_{j}(A_{j}\rho_{j})+\inf\left\\{\,\operatorname{W}(\sigma)\mid\sigma\in
D,\ \Theta(\sigma)=\rho\,\right\\}.$ (5.8)

For instance, the example of Section 3 is recovered by choosing
$\Gamma=\left\\{\,(h,\xi)\in{\Gamma_{0}}\mid h\equiv 1\,\right\\}$,

$\operatorname{F}_{j}(A_{j}\rho)=\frac{1}{2}{\left\lVert A_{j}\rho-
b_{j}\right\rVert}_{2}^{2},\quad\text{and}\quad
w(h,\xi)=\int_{0}^{1}\alpha+\frac{\beta}{2}|\xi^{\prime}(t)|^{2}\mathop{}\\!\mathrm{d}t.$
(5.9)

The choice of $\varphi$ in Bredies et al. [2020b] was equivalent to
$\varphi(\gamma)=\frac{w(\gamma)}{\operatorname{E}(0)}$, whereas we suggest
the default of $\varphi(\gamma)=\frac{\alpha}{\operatorname{E}(0)}$. Both
functions are strictly positive but sufficiently small not to change the
minimisers of $\operatorname{E}$ or $\mathcal{E}$ (i.e.
$\operatorname*{argmin}_{\sigma\in
D}\operatorname{E}(\sigma)=\operatorname*{argmin}_{\sigma\in\mathcal{M}^{+}(\Gamma)}\operatorname{E}(\sigma)$,
see Section 3.2).

The generality of $\Gamma$ allows us to treat different applications with the
same analysis, for example:

$\bf{\Gamma\subset\left\\{\,(h,\xi)\in{\Gamma_{0}}\mid h\equiv 1\,\right\\}}$

This accounts for the example of Section 3. If
$\sigma\in\mathcal{M}^{+}(\Gamma)$, then $\Theta(\sigma)\geq 0$ and mass is
preserved on paths (e.g.
$t\mapsto\int_{\Omega}\mathop{}\\!\mathrm{d}\Theta(\sigma)_{t}$ is constant).

$\bf{\Gamma\subset\left\\{\,(h,\xi)\in{\Gamma_{0}}\mid
h\in\mathscr{C}(\left[0,1\right];\left[0,1\right])\,\right\\}}$

This accounts for the example of Section 4. We still have $\Theta(\sigma)\geq
0$, but $t\mapsto\int_{\Omega}\mathop{}\\!\mathrm{d}\Theta(\sigma)_{t}$ is not
(necessarily) constant. In particular, mass can be created or destroyed
(continuously) at any time.

$\bf{\Gamma\subset\Gamma_{\infty}=\left\\{\,(h,\xi)\in{\Gamma_{0}}\mid{\left\lVert
h\right\rVert}_{\infty}\leq 1\,\right\\}}$

In the general case $\Theta(\sigma)$ is a general signed measure, even when
$\sigma\geq 0$. In words, $\sigma$ can give positive weight to negative mass
curves. The only constraint is that
$\Theta(\sigma)\in\mathscr{C}_{w}(\left[0,1\right];\mathcal{M}(\overline{\Omega}))$
is continuous in time.

Any choice of energy in this framework leads to a sparse reconstruction in the
space of curves. The following theorem is equivalent to that of Theorems 3.1
and 4.1.

######  Theorem 5.1.

If $A_{j}$ are given by (5.1) and
$\varphi,w\colon\Gamma\to\left[0,+\infty\right]$ are lower semi-continuous,
then $\operatorname{F}$ and $\operatorname{E}$ are lower semi-continuous.
Furthermore, if $\operatorname{F}$ is bounded below, $w$ or $\varphi$ have
compact sub-levelsets, and $\inf_{\gamma\in\Gamma}\varphi(\gamma)>0$, then
$\operatorname{E}|_{D}$ has compact sub-levelsets. There exists a choice of
minimiser $\sigma^{*}\in\operatorname*{argmin}_{\sigma\in
D}\operatorname{E}(\sigma)$ such that

$\text{for some }a_{i}\geq 0,\
\gamma^{i}\in\Gamma,\qquad\sigma^{*}=\sum_{i=1}^{s}a_{i}\delta_{\gamma^{i}}$
(5.10)

for some $s\leq m(T+1)+1$. If in addition
$\int_{\Gamma}\varphi\mathop{}\\!\mathrm{d}\sigma^{*}<1$, then $s\leq m(T+1)$.

The proof is given in Appendix A.

###  5.2 Inexact sliding Frank-Wolfe algorithm

As previously stated, we will use a variant of the Frank-Wolfe algorithm to
take advantage of the sparse structure of reconstructions. We will also allow
for inexact computation of the linear oracle, as is necessary in this
application. With a slight abuse of notation, define
$\operatorname{F}^{\prime}\colon D\to\mathscr{C}_{b}(\Gamma)$ so

$\forall\sigma\in
D,\qquad\operatorname{F}^{\prime}(\sigma)=\left[(h,\xi)\mapsto\sum_{j=0}^{T}h(t_{j})\eta_{j}(\xi(t_{j}))\right]\quad\text{where}\quad\eta_{j}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}A_{j}^{*}\nabla\operatorname{F}_{j}(A_{j}(e_{t_{j}})_{\sharp}{\sigma}).$
(5.11)

It can be verified that $\operatorname{F}^{\prime}$ is a representation of the
Gateaux derivative of $\operatorname{F}$. For all $\sigma\in D$ define the
exact linear oracle $\Phi\colon D\to D$ by a choice

$\Phi(\sigma)\in\operatorname*{argmin}_{\widetilde{\sigma}\in
D}\sum_{j=0}^{T}\int_{\Gamma}\left(\operatorname{F}^{\prime}(\sigma)[\gamma]+w(\gamma)\right)\mathop{}\\!\mathrm{d}\widetilde{\sigma}(\gamma).$
(5.12)

Lemma 5.4 will justify that this minimiser is well-defined. The following
result is our main convergence proof for Algorithm 2, which clarifies the
sense of approximation $\mu^{n+1}\approx\Phi(\sigma^{n})$.

######  Theorem 5.2.

Suppose $\\{\sigma^{n},\mu^{n},\lambda_{n}\\}_{n\in\mathds{N}}$ are generated
by Algorithm 2 and there exists
$\sigma^{*}\in\operatorname*{argmin}_{\sigma\in D}\operatorname{E}(\sigma)$.
If $\operatorname{F}$ has bounded curvature, i.e.

$C\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}2\sup_{\begin{subarray}{c}\sigma,\widetilde{\sigma}\in D\\\
\lambda\in\left]0,1\right[\end{subarray}}\frac{\operatorname{F}((1-\lambda)\sigma+\lambda\widetilde{\sigma})-\operatorname{F}(\sigma)-\lambda\int_{\Gamma}\operatorname{F}^{\prime}(\sigma)\mathop{}\\!\mathrm{d}[\widetilde{\sigma}-\sigma]}{\lambda^{2}}<+\infty,$
(5.13)

and

$\liminf_{n\to+\infty}\operatorname{gap}_{n}\leq
0\quad\text{where}\quad\operatorname{gap}_{n}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\int_{\Gamma}\left(\operatorname{F}^{\prime}(\sigma^{n})+w\right)\mathop{}\\!\mathrm{d}[\mu^{n+1}-\sigma^{*}],$
(5.14)

then
$\lim_{n\to+\infty}\operatorname{E}(\sigma^{n})=\operatorname{E}(\sigma^{*})$.

The proof of this theorem is in Section 6.1. Note that
$\operatorname{gap}_{n}$ is equivalent to the standard linear duality gap,
c.f. Jaggi [2013], informally written as
$\operatorname{gap}_{n}=\left\langle\nabla\operatorname{E}(\sigma^{n}),\Phi(\sigma^{n})-\sigma^{*}\right\rangle$.
The motivation for this formulation is to enable the proof of almost-sure
convergence for a stochastic variant.

######  Theorem 5.3.

Let $\sigma^{*}$ and $\\{\sigma^{n},\mu^{n},\lambda_{n}\\}_{n\in\mathds{N}}$
be as in Theorem 5.2 where $\mu^{n}$ is chosen stochastically. Let
$\mathcal{S}^{N}$ be the filtration generated by
$\\{\sigma^{n},\mu^{n}\\}_{0\leq n\leq N}$, and let
$\delta\colon\left[0,+\infty\right[\to\left[0,1\right]$ such that

$\forall n\in\mathds{N},\ \varepsilon\geq
0,\qquad\mathds{P}(\operatorname{gap}_{n}\leq\varepsilon|\mathcal{S}^{n})\geq\delta(\varepsilon)\quad\mathcal{S}^{n}\text{-almost
surely}.$ (5.15)

If $\operatorname{F}$ has bounded curvature and $\delta(\varepsilon)>0$ for
all $\varepsilon>0$, then
$\lim_{n\to+\infty}\operatorname{E}(\sigma^{n})=\operatorname{E}(\sigma^{*})$
almost-surely.

The proof of this theorem is in Section 6.2. The key idea is that for every
$n$ there is uniformly positive chance of finding a good approximation of
$\sigma^{*}$ or $\Phi(\sigma^{n-1})$.

Silveti-Falls et al. [2020] proves a similar convergence result in the
stochastic setting, our random $\operatorname{gap}_{n}$ is an equivalent error
estimate to their “$\lambda^{s}_{k}$”. Their proof requires (assumption P.8)
something similar to
$\sum_{n=0}^{\infty}\lambda_{n}\operatorname{gap}_{n}<+\infty$ where
$\lambda_{n}$ is chosen deterministically, which necessitates the random error
decaying at an a priori known rate. This is unrealistic in our application,
and so we allow $\lambda_{n}$ to be adaptive to stabilise the convergence.

1:Choose $\sigma^{0}\in D$ with $\operatorname{E}(\sigma^{0})<+\infty$,
$n\leftarrow 0$,

2:repeat

3: $\mu^{n+1}\approx\Phi(\sigma^{n})$ $\triangleright$ inexact linear oracle

4:
$\lambda_{n}\leftarrow\operatorname*{argmin}_{\lambda\in\left[0,1\right]}\operatorname{E}((1-\lambda)\sigma^{n}+\lambda\mu^{n+1})$
$\triangleright$ exact line search

5: Choose $\sigma^{n+1}$ such that
$\operatorname{E}(\sigma^{n+1})\leq\operatorname{E}((1-\lambda_{n})\sigma^{n}+\lambda_{n}\mu^{n+1})$
$\triangleright$ Sliding step

6: $n\leftarrow n+1$

7:until converged

Algorithm 2 Inexact sliding Frank-Wolfe algorithm

###  5.3 Inexact linear oracle computation

As in Section 3.4, the most challenging computational component of Algorithm 2
is evaluating the (inexact) linear oracle, i.e. (5.12). Similarly to Section
3.4, we can show that it is equivalent to optimise directly over $\Gamma$.

######  Lemma 5.4.

Suppose $\varphi,w\colon\Gamma\to\left[0,+\infty\right]$ are lower semi-
continuous. If $\inf_{\gamma\in\Gamma}\varphi(\gamma)>0$ and either $\varphi$
or $w$ have compact sub-levelsets, then there exists a choice of
$\Phi(\widetilde{\sigma})$ such that

$\Phi(\widetilde{\sigma})\in\\{0\\}\cup\left\\{\,\varphi(\gamma^{*})^{-1}\delta_{\gamma^{*}}\mid\gamma^{*}\in\operatorname*{argmin}_{\gamma\in\Gamma}\frac{\eta(\gamma)+w(\gamma)}{\varphi(\gamma)}\,\right\\}.$
(5.16)

where $\eta(h,\xi)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\sum_{j=0}^{T}h(t_{j})\eta_{j}(\xi(t_{j}))\in\mathscr{C}_{b}(\Gamma)$.

###### Proof.

Recall the definition of $\Phi$,

$\Phi(\widetilde{\sigma})\in\operatorname*{argmin}_{\sigma\in
D}\widetilde{\operatorname{E}}(\sigma)\quad\text{where}\quad\widetilde{\operatorname{E}}(\sigma)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\int_{\Gamma}\left(\eta+w\right)\mathop{}\\!\mathrm{d}\sigma.$
(5.17)

The properties of $D$ come from Lemma A.5, in particular $D$ is convex, closed
and bounded, so $\inf_{\sigma\in
D}\int_{\Gamma}\eta\mathop{}\\!\mathrm{d}\sigma>-\infty$ and
$\widetilde{\operatorname{E}}$ is lower semi-continuous (Lemma A.1). It is
therefore well-posed to consider minimisers of $\widetilde{\operatorname{E}}$.
We show that there is a choice
$\Phi(\widetilde{\sigma})=\sigma^{*}\in\operatorname{Ext}\left(D\right)$.

Case $D$ is compact:

If $\varphi$ has compact sub-levelsets, then $D$ is compact by Lemma A.5.
Bauer’s principle therefore states that there exists a point

$\sigma^{*}\in\operatorname{Ext}\left(D\right)\quad\text{such
that}\quad\widetilde{\operatorname{E}}(\sigma^{*})=\inf_{\sigma\in
D}\widetilde{\operatorname{E}}(\sigma).$ (5.18)

Else:

Otherwise, $w$ has compact sub-levelsets, so the sub-levelset

$U\stackrel{{\scriptstyle\mbox{\tiny def.}}}{{=}}\left\\{\,\sigma\in
D\quad\mid\quad\widetilde{\operatorname{E}}(\sigma)\leq 1\,\right\\}$ (5.19)

is compact by Theorem A.6, convex, and non-empty because $0\in D$,
$\widetilde{\operatorname{E}}(0)=0$. Application of Bauer’s principle now
gives a point $\sigma^{*}\in\operatorname{Ext}\left(U\right)$, to complete the
proof we will confirm $\sigma^{*}\in\operatorname{Ext}\left(D\right)$.

Suppose for contradiction that there exists $\sigma_{0},\sigma_{1}\in
D\setminus\\{\sigma^{*}\\}$ with
$\sigma^{*}\in\left]\sigma_{0},\sigma_{1}\right[$. Without loss of generality,
we assume that $\sigma^{*}=\frac{1}{2}\sigma_{0}+\frac{1}{2}\sigma_{1}$ and we
define

$\forall\lambda\in\left[0,1\right],\qquad\sigma_{\lambda}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\lambda\sigma_{1}+(1-\lambda)\sigma_{0}\in
D,\qquad\sigma_{\frac{1}{2}}=\sigma^{*}.$ (5.20)

Consider the function
$f(\lambda)=\widetilde{\operatorname{E}}(\sigma_{\lambda})$, which is linear
and
$f(\frac{1}{2})=\widetilde{\operatorname{E}}(\sigma^{*})\leq\widetilde{\operatorname{E}}(0)=0$.
As $f$ is continuous, there exists $\varepsilon>0$ with
$f(\frac{1}{2}\pm\varepsilon)\leq\frac{1}{2}$. We conclude that
$\sigma_{\frac{1}{2}\pm\varepsilon}\in U$ and
$\sigma^{*}\in\left]\sigma_{\frac{1}{2}-\varepsilon},\sigma_{\frac{1}{2}+\varepsilon}\right[$,
contradicting the assumption that
$\sigma^{*}\in\operatorname{Ext}\left(U\right)$.

In both cases we see there is a choice
$\Phi(\widetilde{\sigma})=\sigma^{*}\in\operatorname{Ext}\left(D\right)$.
Finally, by Lemma A.5 we have

$\operatorname{Ext}\left(D\right)=\\{0\\}\cup\left\\{\,\varphi(\gamma)^{-1}\delta_{\gamma}\mid\varphi(\gamma)<+\infty\,\right\\},$
(5.21)

so (5.16) is just the required statement

$\sigma^{*}\in\operatorname*{argmin}_{\sigma\in\operatorname{Ext}\left(D\right)}\widetilde{\operatorname{E}}(\sigma)\subset\\{0\\}\cup\left\\{\,\varphi(\gamma^{*})^{-1}\delta_{\gamma^{*}}\mid\gamma^{*}\in\operatorname*{argmin}_{\gamma\in\Gamma}\widetilde{\operatorname{E}}\left(\varphi(\gamma)^{-1}\delta_{\gamma}\right)\,\right\\}.$
(5.22)

∎

We use this formulation numerically because it is easier to optimise over a
set of curves than a set of measures on curves. Because of this, we define our
discrete minimiser

$\Phi_{d}(\sigma)\in\\{0\\}\cup\left\\{\,\varphi(\gamma_{d}^{*})^{-1}\delta_{\gamma_{d}^{*}}\mid\gamma_{d}^{*}\in\operatorname*{argmin}_{\gamma\in\prod_{j=0}^{T}\Lambda_{j}}\frac{\eta(\gamma)+w(\gamma)}{\varphi(\gamma)}\,\right\\}$
(5.23)

for some discrete sets $\Lambda_{j}\subset\mathds{R}\times\overline{\Omega}$.
For us to compute this discrete minimiser efficiently, we must further assume

$\displaystyle\varphi(\gamma)$ $\displaystyle=\varphi_{0}$
$\displaystyle\text{ for some constant }\varphi_{0}>0,$ (A1) $\displaystyle
w(\gamma)$
$\displaystyle=\sum_{j=1}^{T}\operatorname{step}_{j}(\gamma(t_{j-1}),\gamma(t_{j}))$
$\displaystyle\text{ for some
}\operatorname{step}_{j}\colon(\mathds{R}\times\overline{\Omega})^{2}\to\mathds{R}.$
(A2)

These properties were clear in Section 3.4 where for all
$\gamma=(h,\xi)\in\Gamma$

$\varphi_{0}=\frac{\alpha}{\operatorname{E}(0)},\qquad\operatorname{step}_{j}(\gamma(t_{j-1}),\gamma(t_{j}))=\alpha(t_{j}-t_{j-1})+\frac{\beta}{2}\frac{|\xi(t_{j})-\xi(t_{j-1})|^{2}}{t_{j}-t_{j-1}}.$
(5.24)

In particular, the step function of the Benamou-Brenier penalty is Wasserstein
optimal transport, and the step function of the dynamic Wasserstein-Fisher-Rao
penalty is the static Wasserstein-Fisher-Rao penalty.

As before, it is possible to recover the continuous curve $\gamma$ from the
vector $(\gamma(t_{0}),\ldots,\gamma(t_{T}))$ using geodesics of $w$. If
$\sigma^{*}\in\operatorname*{argmin}_{\sigma\in\mathcal{M}^{+}(\Gamma)}\operatorname{E}(\sigma)$
and $\int_{\Gamma}\varphi\mathop{}\\!\mathrm{d}\sigma^{*}<1$, then
$\sigma^{*}$ is supported on geodesics of $w$ (Lemma B.1). Otherwise, if
$\int_{\Gamma}\varphi\mathop{}\\!\mathrm{d}\sigma^{*}=1$, then the continuous-
time interpolation of $\gamma$ must balance the costs of $w(\gamma)$ and
$\varphi(\gamma)$.

Under assumptions (A1) and (A2), the problem of evaluating $\Phi_{d}$ can be
formulated as finding the minimal path on an ordered, weighted, directed
acyclic graph. The vertices are the points $(t_{j},y)$ for $y\in\Lambda_{j}$,
the edge weights are given by $\eta_{j}$ and $\operatorname{step}_{j}$, and
the time index provides an ordering and prevents cycles in the graph. In
Section 7 we will write out the algorithm fully for clarity, but we state here
the known complexity bound.

######  Theorem 5.5 ([Cormen et al., 2009, Sec. 24.2]).

Suppose assumptions (A1) and (A2) hold and $|\Lambda_{j}|\leq N$ for each $j$,
then $\gamma_{d}^{*}$ can be computed with complexity

$\text{total
time}=O\left(N\sum_{j=0}^{T}\operatorname{cost}(\eta_{j})+N^{2}\sum_{j=1}^{T}\operatorname{cost}(\operatorname{step}_{j})\right)$
(5.25)

where $\operatorname{cost}(\eta_{j})$ is the cost of evaluating $\eta_{j}$
once etc.

######  Remark 5.1.

In graph terminology, the two terms of (5.25) represent the number of vertices
plus the number of edges, asymptotically $O(NT)$ and $O(N^{2}T)$ respectively.
If, for example, the paths have a maximum velocity $V$, then the number of
edges per vertex is reduced from $N$ to $O(N(VT^{-1})^{d})$. This results in a
reduced total complexity of $O(NT+N^{2}T^{1-d}V^{d})$.

##  6 Frank-Wolfe convergence

In this section we analyse Algorithm 2, first proving Theorem 5.2 in Section
6.1, then the stochastic variant Theorem 5.3 is in Section 6.2.

###  6.1 Deterministic Frank-Wolfe analysis

The proof of Theorem 5.2 closely follows the argument of [Jaggi, 2013, Thm. 1]
where $\operatorname{gap}_{n}=O(n^{-1})$, whereas we only assume
$\operatorname{gap}_{n_{k}}=O(k^{-1})$ on some sub-sequence $n_{k}$. This
relaxed assumption on $\operatorname{gap}_{n}$ is accounted for by the
linesearch in $\lambda$. This guarantees both that small steps are made for
bad search directions, and large steps are made for good ones. For example,
for all $n\in\mathds{N}$ and $\mu^{n+1}\in D$, we have a monotone descent,

$\operatorname{E}(\sigma^{n+1})\leq\operatorname{E}((1-\lambda_{n})\sigma^{n}+\lambda_{n}\mu^{n+1})=\inf_{\lambda\in\left[0,1\right]}\operatorname{E}((1-\lambda)\sigma^{n}+\lambda\mu^{n+1})\leq\operatorname{E}(\sigma^{n}).$
(6.1)

The remaining assumption in Theorem 5.2 is the smoothness of
$\operatorname{F}$, i.e. $C<+\infty$ where

$C\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\max\left(\frac{\operatorname{E}(\sigma^{0})-\operatorname{E}(\sigma^{*})}{2},\
2\sup_{\begin{subarray}{c}\sigma,\widetilde{\sigma}\in D\\\
\lambda\in\left]0,1\right[\end{subarray}}\frac{\operatorname{F}((1-\lambda)\sigma+\lambda\widetilde{\sigma})-\operatorname{F}(\sigma)-\lambda\int_{\Gamma}\operatorname{F}^{\prime}(\sigma)\mathop{}\\!\mathrm{d}[\widetilde{\sigma}-\sigma]}{\lambda^{2}}\right).$
(6.2)

For all $n\in\mathds{N}$ and $\lambda\in\left[0,1\right]$, by the convexity of
$\operatorname{F}$ we have

$\displaystyle\operatorname{F}((1-\lambda)\sigma^{n}+$
$\displaystyle\lambda\mu^{n+1})\leq\operatorname{F}(\sigma^{n})+\lambda\int_{\Gamma}\operatorname{F}^{\prime}(\sigma^{n})\mathop{}\\!\mathrm{d}[\mu^{n+1}-\sigma^{n}]+\frac{C\lambda^{2}}{2}$
(6.3)
$\displaystyle=(1-\lambda)\operatorname{F}(\sigma^{n})+\lambda\left[\operatorname{F}(\sigma^{n})+\int_{\Gamma}\operatorname{F}^{\prime}(\sigma^{n})\mathop{}\\!\mathrm{d}[\mu^{n+1}-\sigma^{n}]\right]+\frac{C\lambda^{2}}{2}$
(6.4)
$\displaystyle=(1-\lambda)\operatorname{F}(\sigma^{n})+\lambda\left[\operatorname{F}(\sigma^{n})+\int_{\Gamma}\operatorname{F}^{\prime}(\sigma^{n})\mathop{}\\!\mathrm{d}[\sigma^{*}-\sigma^{n}+\mu^{n+1}-\sigma^{*}]\right]+\frac{C\lambda^{2}}{2}$
(6.5)
$\displaystyle\leq(1-\lambda)\operatorname{F}(\sigma^{n})+\lambda\left[\operatorname{F}(\sigma^{*})+\int_{\Gamma}\operatorname{F}^{\prime}(\sigma^{n})\mathop{}\\!\mathrm{d}[\mu^{n+1}-\sigma^{*}]\right]+\frac{C\lambda^{2}}{2}.$
(6.6)

The linearity of $\operatorname{W}$ gives the same approximation,

$\displaystyle\operatorname{W}((1-\lambda)\sigma^{n}+\lambda\mu^{n+1})$
$\displaystyle=(1-\lambda)\operatorname{W}(\sigma^{n})+\lambda\operatorname{W}(\mu^{n+1})$
(6.7)
$\displaystyle=(1-\lambda)\operatorname{W}(\sigma^{n})+\lambda\left[\operatorname{W}(\sigma^{*})+\operatorname{W}(\mu^{n+1}-\sigma^{*})\right].$
(6.8)

Recalling
$\operatorname{gap}_{n}=\int_{\Gamma}\left(\operatorname{F}^{\prime}(\sigma^{n})+w\right)\mathop{}\\!\mathrm{d}[\mu^{n+1}-\sigma^{*}]$,
we summarise

$\operatorname{E}((1-\lambda)\sigma^{n}+\lambda\mu^{n+1})\leq(1-\lambda)\operatorname{E}(\sigma^{n})+\lambda[\operatorname{E}(\sigma^{*})+\operatorname{gap}_{n}]+\frac{C}{2}\lambda^{2}.$
(6.9)

We can now combine these approximations to complete the proof.

###### Proof of Theorem 5.2..

Since $\liminf_{n\to\infty}\operatorname{gap}_{n}\leq 0$ there exist a
subsequence $(n_{k})_{k\in\mathds{N}}$ such that
$\operatorname{gap}_{n_{k}}\leq\frac{C}{k}$ for all $k\geq 1$. For all
$\lambda\in\left[0,1\right]$,

$\displaystyle\operatorname{E}(\sigma^{n_{k+1}})-\operatorname{E}(\sigma^{*})$
$\displaystyle\leq\operatorname{E}(\sigma^{n_{k}+1})-\operatorname{E}(\sigma^{*})$
(6.10)
$\displaystyle\leq\operatorname{E}((1-\lambda)\sigma^{n_{k}}+\lambda\mu^{n_{k}+1})-\operatorname{E}(\sigma^{*})$
(6.11)
$\displaystyle=(1-\lambda)[\operatorname{E}(\sigma^{n_{k}})-\operatorname{E}(\sigma^{*})]+\lambda\operatorname{gap}_{n_{k}}+\frac{C}{2}\lambda^{2}.$
(6.12)

Now we proceed by induction. Note that
$\operatorname{E}(\sigma^{n_{2}})-\operatorname{E}(\sigma^{*})\leq\operatorname{E}(\sigma^{0})-\operatorname{E}(\sigma^{*})\leq
2C$ by (6.1) and (6.2). Next, if for some $k\geq 2$,
$\operatorname{E}(\sigma^{n_{k}})-\operatorname{E}(\sigma^{*})\leq\frac{4C}{k}$,
we choose $\lambda=\frac{2}{k}\in\left[0,1\right]$, and we observe

$\displaystyle\operatorname{E}(\sigma^{n_{k+1}})-\operatorname{E}(\sigma^{*})$
$\displaystyle\leq\left(1-\frac{2}{k}\right)\frac{4C}{k}+\frac{2}{k}\frac{C}{k}+\frac{C}{2}\frac{4}{k^{2}}$
(6.13) $\displaystyle=\frac{4C}{k}\left(1-\frac{2}{k}+\frac{1}{k}\right)$
(6.14) $\displaystyle\leq\frac{4C}{k}\left(1-\frac{1}{k+1}\right)$ (6.15)
$\displaystyle=\frac{4C}{k+1}.$ (6.16)

We conclude by induction that the inequality holds for all $k\geq 2$. By
monotonicity,
$\lim_{n\to+\infty}\operatorname{E}(\sigma^{n})=\lim_{k\to+\infty}\operatorname{E}(\sigma^{n_{k}})=\operatorname{E}(\sigma^{*})$
as required. ∎

###  6.2 Stochastic Frank-Wolfe analysis

The previous proof shows that Algorithm 2 converges whenever
$\liminf_{n\to+\infty}\operatorname{gap}_{n}=0$. The main aim in this section
is to prove Theorem 5.3 which shows this is satisfied almost surely. We will
then explore two explicit and practical examples in which it applies.

####  6.2.1 Convergence proof

In Theorem 5.3 we assume
$\mathds{P}(\operatorname{gap}_{n}\leq\varepsilon|\mathcal{S}^{n})\geq\delta(\varepsilon)>0$
for all $n\in\mathds{N}$ and $\varepsilon>0$. The first task is to translate
this assumption into a statement about
$\liminf_{n\to+\infty}\operatorname{gap}_{n}$. For fixed $\varepsilon>0$,
observe

$\displaystyle\mathds{P}\left[\liminf_{n\to+\infty}\operatorname{gap}_{n}>\varepsilon\right]$
$\displaystyle=\mathds{P}\left[\bigcup_{N\in\mathds{N}}\left\\{\inf_{n\geq
N}\operatorname{gap}_{n}>\varepsilon\right\\}\right]$ (6.17)
$\displaystyle\leq\sum_{N\in\mathds{N}}\mathds{P}\left[\inf_{n\geq
N}\operatorname{gap}_{n}>\varepsilon\right]$ (sub-additivity) (6.18)
$\displaystyle\leq\sum_{N\in\mathds{N}}\mathds{P}\left[\bigcap_{n\geq
N}\left\\{\operatorname{gap}_{n}>\varepsilon\right\\}\right]$ (6.19)
$\displaystyle=\sum_{N\in\mathds{N}}\lim_{M\to+\infty}\mathds{P}[G_{N}^{M}]$
(monotone convergence) (6.20)

where

$G_{N}^{M}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\bigcap_{n=N}^{M}\left\\{\operatorname{gap}_{n}>\varepsilon\right\\}.$
(6.21)

The probability of these events can be bounded inductively for $M\geq N$,

$\displaystyle\mathds{P}(G_{N}^{M+1})$
$\displaystyle=\mathds{E}\left[\mathds{E}(\mathds{1}_{G_{N}^{M}\cap\\{\operatorname{gap}_{M+1}>\varepsilon\\}}|\mathcal{S}^{M})\right]$
(6.22)
$\displaystyle=\mathds{E}\left[\mathds{E}(\mathds{1}_{G_{N}^{M}}\cdot\mathds{1}_{\\{\operatorname{gap}_{M+1}>\varepsilon\\}}|\mathcal{S}^{M})\right]$
(6.23)
$\displaystyle=\mathds{E}\left[\mathds{1}_{G_{N}^{M}}\cdot\mathds{E}(\mathds{1}_{\\{\operatorname{gap}_{M+1}>\varepsilon\\}}|\mathcal{S}^{M})\right]$
$\displaystyle(G_{N}^{M}\in\mathcal{S}^{M})$ (6.24)
$\displaystyle=\mathds{E}\left[\mathds{1}_{G_{N}^{M}}\cdot\mathds{P}(\operatorname{gap}_{M+1}>\varepsilon|\mathcal{S}^{M})\right]$
(def. of $\mathds{P}$) (6.25)
$\displaystyle\leq\mathds{E}\left[\mathds{1}_{G_{N}^{M}}\cdot(1-\delta(\varepsilon))\right]$
(def. of $\delta$) (6.26)
$\displaystyle=(1-\delta(\varepsilon))\mathds{P}(G_{N}^{M})\ \ldots$ (6.27)
$\displaystyle\leq(1-\delta(\varepsilon))^{M+1-N}.$ (6.28)

We can now complete the proof.

###### Proof of Theorem 5.3..

To see
$\lim_{n\to+\infty}\operatorname{E}(\sigma^{n})=\operatorname{E}(\sigma^{*})$
almost-surely, observe

$\displaystyle\mathds{P}(\lim_{n\to+\infty}\operatorname{E}(\sigma^{n})=\operatorname{E}(\sigma^{*}))$
$\displaystyle\geq\mathds{P}(\liminf_{n\to+\infty}\operatorname{gap}_{n}\leq
0)$ (Theorem 5.2) (6.29)
$\displaystyle=1-\mathds{P}(\liminf_{n\to+\infty}\operatorname{gap}_{n}>0)$
(6.30) $\displaystyle=1-\lim_{\varepsilon\to
0}\mathds{P}(\liminf_{n\to+\infty}\operatorname{gap}_{n}>\varepsilon)$
$\displaystyle\begin{subarray}{c}\text{(monotone}\\\
\text{convergence)}\end{subarray}$ (6.31) $\displaystyle\geq
1-\lim_{\varepsilon\to
0}\sum_{N\in\mathds{N}}\lim_{M\to+\infty}\mathds{P}\left[G_{N}^{M}\right]$
$\displaystyle\eqref{eq: subadditive prob}$ (6.32) $\displaystyle\geq
1-\lim_{\varepsilon\to
0}\sum_{N\in\mathds{N}}\lim_{M\to+\infty}(1-\delta(\varepsilon))^{M-N}$
$\displaystyle\eqref{eq: G_N^M to infty}$ (6.33) $\displaystyle=1$ (6.34)

as required. ∎

####  6.2.2 Independent random gaps variant

In this section we focus on the case

$\mu^{n}=\widetilde{\Phi}(\sigma^{n-1},z^{n})$ (6.35)

for some approximate linear oracle $\widetilde{\Phi}$, every $n\in\mathds{N}$,
where $z^{n}\sim Z$ is a random instance of an ambient probability density. In
the context of Section 5.3, $z^{n}$ could be a random discretisation of
$\Gamma$ and $\mu^{n}$ is the discrete minimiser. The following corollary
shows that, if $z^{n}$ are independent and have a uniformly positive
probability of being close to $\sigma^{*}$, then the conditions of Theorem 5.3
are satisfied.

######  Corollary 6.1.

Suppose $\mu^{n}=\widetilde{\Phi}(\sigma^{n-1},z^{n})$ for some $z^{n}\sim Z$
and define $\widetilde{\delta}\colon\left[0,+\infty\right[\to\left[0,1\right]$
by

$\widetilde{\delta}(\varepsilon)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\inf_{\sigma\in D}\mathds{P}_{z\sim
Z}\left(\int_{\Gamma}\operatorname{F}^{\prime}(\sigma)+w\mathop{}\\!\mathrm{d}[\widetilde{\Phi}(\sigma,z)-\sigma^{*}]\leq\varepsilon\right).$
(6.36)

If $z^{n}$ is independent of the filtration $\widetilde{\mathcal{S}}^{n-1}$
generated by $\\{z^{k},\sigma^{k}\\}_{k=0}^{n-1}$ for all
$n\in\mathds{N}^{*}$, then $\widetilde{\delta}$ is a valid choice for $\delta$
in Theorem 5.3.

######  Lemma 6.2.

Suppose $X$ is an $\mathcal{S}$-measurable random variable and $Z\perp X$
(i.e. $Z$ is independent of $X$) is another random variable. Then, for all
integrable Borel measurable functions $g$

$\mathds{E}(g(X,Z)|\mathcal{S})=\mathds{E}\left[\mathds{E}_{Z}(g(X,Z))|\mathcal{S}\right].$
(6.37)

This is a standard result in probability theory for which a proof is provided
in Appendix C.

###### Proof of Corollary 6.1..

We need to link the values of $\delta$ from Theorem 5.3 with
$\widetilde{\delta}$. Recall that
$\operatorname{gap}_{n}=\int_{\Gamma}\left(\operatorname{F}^{\prime}(\sigma^{n})+w\right)\mathop{}\\!\mathrm{d}[\widetilde{\Phi}(\sigma^{n},z^{n+1})-\sigma^{*}]$
is measurable w.r.t. the filtration generated by
$\\{\widetilde{\mathcal{S}}^{n},z^{n+1}\\}$, with
$\widetilde{\mathcal{S}}^{n}\perp z^{n+1}$. We therefore have

$\displaystyle\mathds{P}(\operatorname{gap}_{n}\leq\varepsilon|\mathcal{S}^{n})$
$\displaystyle\equiv\mathds{E}\left[\mathds{1}_{\operatorname{gap}_{n}\leq\varepsilon}\,|\mathcal{S}^{n}\right]$
(6.38)
$\displaystyle=\mathds{E}\left[\mathds{E}(\mathds{1}_{\operatorname{gap}_{n}\leq\varepsilon}|\widetilde{\mathcal{S}}^{n})\,|\mathcal{S}^{n}\right]$
$\displaystyle(\widetilde{\mathcal{S}}^{n}\supset\mathcal{S}^{n})$ (6.39)
$\displaystyle=\mathds{E}\left[\mathds{E}(\mathds{E}_{z^{n+1}\sim
Z}(\mathds{1}_{\operatorname{gap}_{n}\leq\varepsilon})|\widetilde{\mathcal{S}}^{n})\,|\mathcal{S}^{n}\right]$
(Lemma 6.2) (6.40)
$\displaystyle\geq\mathds{E}\left[\mathds{E}(\widetilde{\delta}(\varepsilon)|\widetilde{\mathcal{S}}^{n})\,|\mathcal{S}^{n}\right]$
(def. of $\widetilde{\delta}$) (6.41)
$\displaystyle=\widetilde{\delta}(\varepsilon).$ (6.42)

This bound is uniform in $n$, therefore $\widetilde{\delta}$ satisfies the
definition of $\delta$. ∎

A natural choice of $\widetilde{\Phi}$ might be where $z^{n}$ is a random
discretisation and $\widetilde{\Phi}$ is the discrete minimiser from (5.23).
In the context of this work, we expect $\operatorname{supp}(\sigma^{*})$ to be
a discrete set and can choose
$z^{n}\subset(\left[-1,1\right]\times\overline{\Omega})^{T+1}$ (which is
compact). If $|z^{n}|\geq|\operatorname{supp}(\sigma^{*})|$, then it should be
straight-forward to show $\mathds{P}_{z^{n}\sim
Z}\left(\int_{\Gamma}\operatorname{F}^{\prime}(\sigma^{n})+w\mathop{}\\!\mathrm{d}[\widetilde{\Phi}(\sigma^{n-1},z^{n})-\sigma^{*}]\leq\varepsilon\right)>0$.
To show $\widetilde{\delta}(\varepsilon)>0$, the final challenge is to show
uniformity in $\sigma^{n-1}$. For example, this follows if
$\operatorname{F}^{\prime}(\sigma)+w\in\mathscr{C}((\left[-1,1\right]\times\overline{\Omega})^{T+1})$
is Lipschitz with Lipschitz constant independent of $\sigma$.

Another natural approach would be, for known $\sigma^{n-1}$, to choose for
$z^{n}$ with guaranteed approximation properties for $\Phi(\sigma^{n-1})$,
rather than the fixed but unknown target of $\sigma^{*}$. This approach is
addressed by the following corollary.

######  Corollary 6.3.

Suppose $\mu^{n}=\Phi(\sigma^{n-1},z^{n})$ for some $z^{n}\sim Z$ and define
$\widetilde{\delta}\colon\left[0,+\infty\right[\to\left[0,1\right]$ by

$\widetilde{\delta}(\varepsilon)\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\inf_{\sigma\in D}\mathds{P}_{z\sim
Z}\left(\int_{\Gamma}\left(\operatorname{F}^{\prime}(\sigma)+w\right)\mathop{}\\!\mathrm{d}[\Phi(\sigma,z)-\Phi(\sigma)]\leq\varepsilon\right).$
(6.43)

If $z^{n}$ is independent of the filtration $\widetilde{\mathcal{S}}^{n-1}$
generated by $\\{z^{k},\sigma^{k}\\}_{k=0}^{n-1}$ for all
$n\in\mathds{N}^{*}$, then $\widetilde{\delta}$ is a valid choice for $\delta$
in Theorem 5.3.

###### Proof.

This corollary is a strictly weaker statement than the previous. Fix
$\sigma\in D$ and define $\eta\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\operatorname{F}^{\prime}(\sigma)+w$, then

$\displaystyle\int_{\Gamma}\eta\mathop{}\\!\mathrm{d}[\Phi(\sigma,z)-\Phi(\sigma)]\leq\varepsilon$
$\displaystyle\iff\int_{\Gamma}\eta\mathop{}\\!\mathrm{d}[\Phi(\sigma,z)-\sigma^{*}]\leq\int_{\Gamma}\eta\mathop{}\\!\mathrm{d}[\Phi(\sigma)-\sigma^{*}]+\varepsilon$
(6.44)
$\displaystyle\iff\int_{\Gamma}\eta\mathop{}\\!\mathrm{d}[\Phi(\sigma,z)-\sigma^{*}]\leq\inf_{\widetilde{\sigma}\in
D}\int_{\Gamma}\eta\mathop{}\\!\mathrm{d}[\widetilde{\sigma}-\sigma^{*}]+\varepsilon$
(6.45)
$\displaystyle\implies\int_{\Gamma}\eta\mathop{}\\!\mathrm{d}[\Phi(\sigma,z)-\sigma^{*}]\leq\varepsilon.$
(6.46)

We conclude that the assumptions of Corollary 6.1 hold, and so the conclusion
is the same. ∎

In this case, we see that $\widetilde{\Phi}$ only needs to approximate the
exact linear oracle $\Phi(\sigma^{n-1})$, so we can have convergence by
choosing random discretisations
$z^{n}\subset(\left[-1,1\right]\times\overline{\Omega})^{T+1}$ with
$|z^{n}|\geq|\operatorname{supp}(\Phi(\sigma^{n-1}))|=1$. In our setting, the
criterion of this second corollary is much easier to verify.

##  7 General dynamical programming method

As stated at the end of Section 5.3, we compute optimal curves of (5.23) using
the standard algorithm for minimal paths on ordered, weighted, directed,
acyclic graphs described in [Cormen et al., 2009, Sec. 24.2]. In this section
we expand on the precise computations for a fixed $\sigma\in D$. Recall for
all $\gamma\in\Gamma\subset\Gamma_{\infty}$ we have
$\gamma(t)\in\left[-1,1\right]\times\overline{\Omega}$ for all
$t\in\left[0,1\right]$. Under assumptions (A1) and (A2), the optimisation
problem can therefore be written as

$\displaystyle\gamma_{d}^{*}$
$\displaystyle\in\operatorname*{argmin}_{\gamma}\left\\{\,\sum_{j=0}^{T}\eta_{j}(\gamma_{j})+\sum_{j=1}^{T}\operatorname{step}_{j}(\gamma_{j-1},\gamma_{j})\mid\gamma\in\prod_{j=0}^{T}\Lambda_{j}\,\right\\}$
(7.1)
$\displaystyle\quad\text{where}\quad\eta_{j}(h_{j},\xi_{j})\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}h_{j}[A_{j}^{*}\nabla\operatorname{F}_{j}(A_{j}(e_{t_{j}})_{\sharp}{\sigma})](\xi_{j})$
(7.2)

and $\Lambda_{j}\subseteq\left[-1,1\right]\times\overline{\Omega}$ are finite
sets (“grids” in the mass-location space), $j=0,\ldots,T$. Similarly to
Section 3.4, we define the functions

$\displaystyle\widetilde{\operatorname{E}}_{J}(\gamma)$
$\displaystyle\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\sum_{j=0}^{J}\eta_{j}(\gamma_{j})+\sum_{j=1}^{J}\operatorname{step}_{j}(\gamma_{j-1},\gamma_{j}),$
(7.3) $\displaystyle Y^{*}[y,J]$
$\displaystyle\in\operatorname*{argmin}_{Y}\left\\{\,\widetilde{\operatorname{E}}_{J}(Y)\mid
Y\in\prod_{j=0}^{J}\Lambda_{j},\ Y_{J}=y\,\right\\}$ (7.4)

for each $J=0,\ldots,T$ and $y\in\Lambda_{J}$. This generates $\gamma_{d}^{*}$
through the computation

$\gamma_{d}^{*}=Y^{*}[y^{*},T],\qquad
y^{*}\in\operatorname*{argmin}_{y\in\Lambda_{T}}\widetilde{\operatorname{E}}_{T}(Y^{*}[y,T]).$
(7.5)

Observe that for all $J=1,\ldots,T$ and $y\in\Lambda_{J}$,

$\displaystyle\min_{\begin{subarray}{c}Y\in\prod_{j=0}^{J}\Lambda_{j}\\\
Y_{J}=y\end{subarray}}\widetilde{\operatorname{E}}_{J}(Y)$
$\displaystyle=\min_{\begin{subarray}{c}Y\in\prod_{j=0}^{J}\Lambda_{j}\end{subarray}}\left\\{\widetilde{\operatorname{E}}_{J-1}(Y)+\eta_{J}(y)+\operatorname{step}_{J}(Y_{J-1},y)\right\\}$
(7.6)
$\displaystyle=\eta_{J}(y)+\min_{\begin{subarray}{c}y^{\prime}\in\Lambda_{J-1}\end{subarray}}\Bigg{\\{}\operatorname{step}_{J}(y^{\prime},y)+\min_{\begin{subarray}{c}Y\in\prod_{j=0}^{J-1}\Lambda_{j}\\\
Y_{J-1}=y^{\prime}\end{subarray}}\widetilde{\operatorname{E}}_{J-1}(Y)\Bigg{\\}}.$
(7.7)

In particular, we can choose $Y^{*}[y,J]$ inductively to be

$Y^{*}[y,J]=\begin{pmatrix}Y^{*}[y^{\prime},J-1]&y\end{pmatrix}\quad\text{for}\quad
y^{\prime}\in\operatorname*{argmin}_{y^{\prime}\in\Lambda_{J-1}}\left[\operatorname{step}_{J}(y^{\prime},y)+\widetilde{\operatorname{E}}_{J-1}(Y^{*}[y^{\prime},J-1])\right].$
(7.8)

##  8 Numerical results

For numerical experiments we implement variants of Algorithm 2 using the
linear-oracle strategy discussed in Section 7. The expanded form of this
algorithm is given in Algorithm 3. The choice of $A_{j}$,
$\operatorname{F}_{j}$, $t_{j}$, and $\operatorname{step}_{j}$ define the
optimisation problem, then the final choices of $\Lambda_{j}$ and $k$ dictate
the variant of the algorithm. We always choose the sliding step to select a
local minimum of $\operatorname{E}$ (or $\widetilde{\operatorname{E}}$) as
computed by an implementation of the L-BFGS algorithm. All code to reproduce
the results and figures in this work can be found
online111https://gitlab.inria.fr/rtovey/DP-for-dynamic-IPs.

The classical sliding Frank-Wolfe algorithm (Algorithm 1) can be recovered by
choosing $\Lambda_{j}=\left[-1,1\right]\times\overline{\Omega}$ and $k=1$. The
potential for $k>1$ was considered in [Bredies et al., 2020b, Sec 5.1.5] as a
“multistart” parameter. The idea is that the approximation of optimal curves
(e.g. $Y^{*}[\cdot,T]$) is very expensive, so some of the other near-optimal
curves should also be used to improve efficiency of the algorithm. Suppose
$\Lambda_{j}$ is chosen to approximate
$\left[-1,1\right]\times\overline{\Omega}$ with a grid of $M$ masses in
$\left[-1,1\right]$ and $N^{d}$ points in $\overline{\Omega}$. From Theorem
5.5 we know that the computational complexity of the linear oracle at every
iteration is $O((N^{d}M)^{2})$. If the “true curves’ are easy to find, then we
can hope to reduce this to $O(\frac{N^{2d}M^{2}}{k})$ per curve. An
alternative motivation is provided by Ding et al. [2020] who propose the
“$k$FW” algorithm. They argue that the inclusion of multiple extreme points at
each iteration allows the algorithm to move directly towards an extreme face,
rather than zig-zagging between extreme points. Note that a similar
observation is made in Flinth et al. [2021] in the context of sparse spike
recovery, where a collection of Dirac masses are added at each iteration. The
sliding step in line 13 of Algorithm 3 further addresses this property,
allowing non–zig-zagging movement.

We consider two families of algorithms implementing Algorithm 3, one
stochastic and the other deterministic. Throughout this section we fix
$\Omega\stackrel{{\scriptstyle\mbox{\tiny def.}}}{{=}}\left]0,1\right[^{2}$
and only consider non-negative curves $h\geq 0$, so the algorithms can be
stated as:

$(k,N,M)$-random mesh

The “multistart” parameter is $k$. For each $n$ and $j$ we generate new
independent uniformly random points
$H=\\{h_{i}^{j}\sim\mathcal{U}[\left[0,1\right]]\\}_{i=1}^{M}$,
$X=\\{x_{i}^{j}\sim\mathcal{U}[\left[0,1\right]^{2}]\\}_{i=1}^{N^{2}}$ and

$\Lambda_{j}=\left\\{\,(h,x)\mid h\in H,\ x\in X\,\right\\}.$ (8.1)

When $M=0$ we take $H=\\{1\\}$ (i.e. a balanced mesh).

$(k,N,M)$-uniform mesh

The “multistart” parameter is $k$. For each $n$ we choose
$\Lambda_{0}=\ldots=\Lambda_{T}$ to be of the form

$\Lambda_{j}=\Lambda_{j}(\widetilde{N})\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left\\{\,(h,(x_{1},x_{2}))\mid h\in\\{0,M^{-1},\ldots,1\\},\
x_{1},x_{2}\in\\{0,\widetilde{N}^{-1},\ldots,1\\}\,\right\\}$ (8.2)

for some $\widetilde{N}\geq 1$. We start with $\widetilde{N}=16$ at $n=0$ then
increment $\widetilde{N}\leftarrow 2\widetilde{N}$ whenever
$\operatorname{E}(\sigma^{n+1})=\operatorname{E}(\sigma^{n})$ (tested after
line 13 of Algorithm 3). This process is terminated once $\widetilde{N}>N$.
Again, if $M=0$ then we only allow the balanced mass $h=1$.

For all problems, $(1,+\infty,+\infty)$-random and
$(1,+\infty,+\infty)$-uniform mesh algorithms are equivalent to the exact
Algorithm 1. For balanced problems, such as in Section 3, it is sufficient to
use the triplet $(1,\infty,0)$. The random algorithms with $N,M\geq 1$ are
guaranteed to converge (eventually) to the exact minimiser by Theorem 5.3. On
the other hand, the uniform path algorithms have nice computational properties
which allows $Y^{*}[\cdot,T]$ to be computed more efficiently at larger $N$.
Whilst all results for random meshes are asymptotic, the uniform path
algorithms also provide a clear stopping criterion. The algorithm stops at
iteration $n$ when
$\operatorname{E}(\sigma^{n+1})=\operatorname{E}(\sigma^{n})$, often this will
mean $\sigma^{n+1}=\sigma^{n}$. Let $D_{d}\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left\\{\,\sigma\in
D\mid\operatorname{supp}(\sigma)\subset\prod_{j=0}^{T}\Lambda_{j}(N)\,\right\\}$
be the discretisation of $D$. Expanding on the linesearch computation, due to
the convexity and smoothness of $\operatorname{E}$,

$\forall\lambda\in\left[0,1\right],\qquad\operatorname{E}((1-\lambda)\sigma^{n}+\lambda\mu^{n+1})-\operatorname{E}(\sigma^{n})=\lambda\int_{\Gamma}(\operatorname{F}^{\prime}(\sigma^{n})+w)\mathop{}\\!\mathrm{d}[\mu^{n+1}-\sigma^{n}]+O(\lambda^{2}).$
(8.3)

If the linesearch terminates with $\lambda=0$, then clearly
$\int_{\Gamma}(\operatorname{F}^{\prime}(\sigma^{n})+w)\mathop{}\\!\mathrm{d}[\mu^{n+1}-\sigma^{n}]\geq
0$. The optimality given by the dual gap (see Jaggi [2013]), can therefore be
stated

$\operatorname{E}(\sigma^{n})-\min_{\sigma\in
D_{d}}\operatorname{E}(\sigma)\leq\sup_{\sigma\in
D_{d}}\int_{\Gamma}(\operatorname{F}^{\prime}(\sigma^{n})+w)\mathop{}\\!\mathrm{d}[\sigma^{n}-\sigma]=\int_{\Gamma}(\operatorname{F}^{\prime}(\sigma^{n})+w)\mathop{}\\!\mathrm{d}[\sigma^{n}-\mu^{n+1}]\leq
0.$ (8.4)

We conclude that $\sigma^{n}$ is at least optimal up to a spatial resolution
of $\frac{1}{N}$.

1:Set $\sigma^{0}\leftarrow 0\in\mathcal{M}^{+}(\Gamma)$, $n\leftarrow 0$, fix
$k\in\mathds{N}$

2:repeat

3: Let
$\eta_{j}=A_{j}^{*}\nabla\operatorname{F}_{j}(A_{j}(e_{t_{j}})_{\sharp}{\sigma^{n}})\in\mathscr{C}(\overline{\Omega})$
for $j=0,\ldots,T$,
$\widetilde{\operatorname{E}}((h,\xi))=\sum_{j=0}^{T}\eta_{j}(\xi_{j})h_{j}+\sum_{j=1}^{T}\operatorname{step}_{j}((h_{j-1},\xi_{j-1}),(h_{j},\xi_{j}))$

4: Choose $\Lambda_{j}\subset\left[0,1\right]\times\overline{\Omega}$,
$j=0,\ldots T$ $\triangleright$ the discrete mesh

5: Compute $Y^{*}[y,T]\in\prod_{j=0}^{T}\Lambda_{j}$ for all
$y\in\Lambda_{T}$, as in (7.8) $\triangleright$ discrete optimal paths

6: Choose
$\widetilde{\gamma}^{1},\ldots,\widetilde{\gamma}^{k}\in\operatorname{Image}(Y^{*}[\cdot,T])$
with least energy in $\widetilde{\operatorname{E}}$ $\triangleright$ select
best $k$ endpoints

7: Find
$\gamma^{1},\ldots\gamma^{k}\in\left[0,1\right]\times\overline{\Omega}$ with
$\widetilde{\operatorname{E}}(\gamma^{i})\leq\widetilde{\operatorname{E}}(\widetilde{\gamma}^{i})$
$\triangleright$ sliding step on linearised problem

8: Set $\sigma\leftarrow\sigma^{n}$, re-order index $i$ such that
$\widetilde{\operatorname{E}}(\gamma^{1})\leq\widetilde{\operatorname{E}}(\gamma^{2}),\ldots$

9: for $i=1,\ldots,k$ do

10:
$\lambda\leftarrow\operatorname*{argmin}_{\lambda\in\left[0,1\right]}\operatorname{E}((1-\lambda)\sigma+\lambda\varphi(\gamma^{i})^{-1}\delta_{\gamma^{i}})$
$\triangleright$ exact line search

11:
$\sigma\leftarrow(1-\lambda)\sigma+\lambda\varphi(\gamma^{i})^{-1}\delta_{\gamma^{i}}$

12: end for

13: Choose $\sigma^{n+1}$ such that
$\operatorname{E}(\sigma^{n+1})\leq\operatorname{E}(\sigma)$ $\triangleright$
sliding step on exact problem

14: $n\leftarrow n+1$

15:until converged

Algorithm 3 Inexact sliding Frank-Wolfe algorithm

###  8.1 Benamou-Brenier example

First we compare directly with the numerical results presented in Bredies et
al. [2020b] for the model discussed in Section 3. In particular, the energy we
seek to minimise is
$\operatorname{E}\colon\mathcal{M}^{+}(\Gamma)\to\left]-\infty,+\infty\right]$
defined by

$\operatorname{E}(\sigma)=\frac{1}{2}\sum_{j=0}^{T}{\left\lVert
A_{j}(e_{t_{j}})_{\sharp}{\sigma}-b_{j}\right\rVert}_{2}^{2}+\int_{\Gamma}\int_{0}^{1}\left(\alpha+\frac{\beta}{2}|\xi^{\prime}(t)|^{2}\right)\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}\sigma(h,\xi)$
(8.5)

where $t_{j}=\frac{j}{T}$ for each $j=0,\ldots,T$,

$\Gamma\stackrel{{\scriptstyle\mbox{\tiny
def.}}}{{=}}\left\\{\,(h,\xi)\in\\{1\\}\times\operatorname{AC}^{2}(\left[0,1\right];\overline{\Omega})\mid\mbox{$\xi^{\prime}$
is constant on $\left]t_{j-1},t_{j}\right[$, $1\leq j\leq T$}\,\right\\},$
(8.6)

and $A_{j}\colon\mathcal{M}(\left[0,1\right]^{2})\to\mathds{R}^{m}$ represents
a finite number of smoothed Fourier samples. The precise details are given in
[Bredies et al., 2020b, Sec. 6]. In view of the convexity of $\Omega$, (8.5)
reformulates as

$\operatorname{E}(\sigma)=\frac{1}{2}\sum_{j=0}^{T}{\left\lVert
A_{j}(e_{t_{j}})_{\sharp}{\sigma}-b_{j}\right\rVert}_{2}^{2}+\int_{\Gamma}\left(\sum_{j=1}^{T}(t_{j}-t_{j-1})\left(\alpha+\frac{\beta}{2}\frac{|\xi(t_{j})-\xi(t_{j-1})|^{2}}{(t_{j}-t_{j-1})^{2}}\right)\right)\mathop{}\\!\mathrm{d}\sigma(h,\xi).$
(8.7)

The two phantoms are also from [Bredies et al., 2020b, Sec. 6]. In the
notation of this work, we would say that, for example, phantom 1 is
represented by

$\displaystyle\sigma_{1}^{B}=\delta_{(h^{1},\xi^{1})}+\delta_{(h^{2},\xi^{2})}\quad\text{where}\quad
h^{1}=h^{2}=1,\ \xi^{1}(t)$ $\displaystyle=(0.2+0.6t,\ 0.2+0.6t),$ (8.8)
$\displaystyle\xi^{2}(t)$ $\displaystyle=(0.8-0.6t,\ 0.2+0.6t).$ (8.9)

Similarly phantom 2 is the sum of 3 Dirac masses, both phantoms are shown here
in Figure 1. In the setting of Section 5.3, we choose $\varphi_{0}=0.1$ (as 10
is much larger than 2 or 3 which is the mass of phantoms 1 and 2
respectively), and

$\operatorname{step}_{j}((h_{0},\xi_{0}),(h_{1},\xi_{1}))=\alpha(t_{j}-t_{j-1})+\frac{\beta}{2}\frac{|\xi_{1}-\xi_{0}|^{2}}{t_{j}-t_{j-1}}=\frac{\alpha}{T}+\frac{\beta
T}{2}|\xi_{1}-\xi_{0}|^{2}.$ (8.10)

For phantom 1 we have $\alpha=\beta=0.5$, $T=21$, and $\alpha=\beta=0.1$,
$T=51$ for phantom 2. The algorithm of Bredies et al. [2020b] found
online222https://github.com/panchoop/DGCG_algorithm/commit/553a564fd8641abcfac6067ebf51a900a6a91d0f
is run with original parameters as a baseline, although with a time limit of 5
days when necessary. This is compared with multiple variants of the random and
uniform algorithms described at the beginning of the section. The uniform
algorithms are run to convergence and the random variants are run for 100 and
10,000 iterations for phantoms 1 and 2 respectively. Figure 2(a) shows the
reconstructions of each algorithm with default parameters, each image is
visually equivalent. The convergence behaviour is shown in Figure 2(b). The
energy plots confirm that each reconstruction has approximately the same
energy, although we find that the random mesh algorithm finds the lowest
energy, closely followed by the uniform mesh. Similarly, the sparsity of each
final reconstruction is equal. The greatest difference between algorithms is
run-time. The random and uniform mesh algorithms are over 100 times faster
than that of Bredies et al. [2020b] in both examples.

We also replicate the noise-scenarios for phantom 2 as tested in Bredies et
al. [2020b]. Our only modification of the baseline algorithm is to remove the
early-stopping routine and run the algorithm for the minimum of 21 iterations
or 5 days (c.f. [Bredies et al., 2020b, Table 1]). We again use the
$(1,25,0)$-random mesh and $(1,256,0)$-uniform mesh algorithms for comparison
in each example. In the three noisy scenarios we add 20%, 60%, and 60%
Gaussian white noise to the data respectively. The first two scenarios use
$\alpha=\beta=0.1$ while the third scenario uses $\alpha=\beta=0.3$. Seen in
Figure 3, both the random and uniform mesh algorithms converge with similar
rates, while the algorithm of Bredies et. al. is still at least 100 times
slower. We see the expected behaviour that the uniform variant converges to a
(possibly non-optimal) energy. As predicted by Theorem 5.3, the random variant
often finds an even lower energy, despite having a much smaller value of $N$.

Figure 1: Synthetic phantoms of the form
$\sigma=\sum_{i}\delta_{(h^{i},\xi^{i})}$ for balanced (top row) and
unbalanced (bottom row) examples. Colour indicates the time $t$, the solid
line indicates the positions $\xi^{i}(t)$, and the width of the overlayed band
is proportional to $h^{i}(t)$.

(a) Final reconstructions visualised equivalently to those in Figure 1.

(b) Convergence of energy and sparsity of reconstructions. For each phantom,
every energy is translated by the smallest energy found by any method.

Figure 2: Comparison of algorithm from Bredies et al. [2020b], the
$(1,25,0)$-random mesh algorithm, and $(1,256,0)$-uniform mesh algorithm