\section{Introduction}\label{sec:introduction}}
\else
\section{Introduction}
\label{sec:introduction}
\fi

\IEEEPARstart{T}{his} demo file is intended to serve as a ``starter file''
for IEEE Computer Society journal papers produced under \LaTeX\ using
IEEEtran.cls version 1.8b and later.
I wish you the best of success.

\hfill mds

\hfill August 26, 2015

\subsection{Subsection Heading Here}
Subsection text here.

\subsubsection{Subsubsection Heading Here}
Subsubsection text here.

\section{Conclusion}
The conclusion goes here.

\appendices
\section{Proof of the First Zonklar Equation}
Appendix one text goes here.

\section{}
Appendix two text goes here.

\ifCLASSOPTIONcompsoc

  \section*{Acknowledgments}
\else

  \section*{Acknowledgment}
\fi

The authors would like to thank...

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\section{Introduction}
This demo file is intended to serve as a ``starter file''
for IEEE Computer Society conference papers produced under \LaTeX\ using
IEEEtran.cls version 1.8b and later.
I wish you the best of success.

\hfill mds

\hfill August 26, 2015

\subsection{Subsection Heading Here}
Subsection text here.

\subsubsection{Subsubsection Heading Here}
Subsubsection text here.

\section{Conclusion}
The conclusion goes here.

\ifCLASSOPTIONcompsoc

  \section*{Acknowledgments}
\else

  \section*{Acknowledgment}
\fi

The authors would like to thank...

\section{Introduction}\label{sec:introduction}}
\IEEEPARstart{L}{arge} Scale Visual Geo-Localization (LSVGL) has been thought of as a key factor in automotive driving, augmented reality, robotics navigation, etc, and the past decade has seen the rapid development of LSVGL in computer vision [1, 2] and robotics communities [3, 4]. LSVGL has become a very important technical means of urban computing [5, 6]. Despite its long success, there is an increasing concern in LSVGL due to the varieties in the real world.

LSVGL has always been treated as an instance retrieval task, which proposes to find the most similar image by querying a large trained geo-tagged database, and almost all the LSVGL systems in recent years try to investigate the most satisfying network architecture to extract a more discriminate image representation. Two steps are generally needed to achieve the geometric location: 1). Training a large scale geo-tagged database; 2). Extracting the descriptor and estimating the optimal matching.
\begin{figure}
	\centering
	\includegraphics[scale=0.35]{graph/quintuplet-loss.pdf}
	\caption{a. The classic triplet loss selects a positive sample and a negative sample. During the optimization process, the optimizer makes the anchor as close as possible to the positive sample and away from the negative sample. b. On the basis of triplet loss, the quadruplet loss selects an additional negative sample to make the anchor farther away from all the negative samples. c. We propose a new loss called quintuplet loss specifically for the LSVGL task. Different from the classic triplet loss and quadruplet loss, two positive samples are selected for the quintuplet loss (two are selected in the diagram, we will prove that it is best to choose two positive samples, both theoretically and experimentally) to strengthen the positive-anchor constraints while retaining the quadruplet loss selection strategy to pick negative samples.}
	\label{fig:cmp_with_other_loss}
\end{figure}

The sophisticated hand-crafted local invariant feature extractors contributed significantly to LSVGL in the early stages, such as  Scale-Invariant Feature Transformation (SIFT) and Speed-up Robust Features (SURF), and some global image features, as GIST, bag-of-words approaches. These traditional feature extraction techniques combined with machine learning algorithms such as SVM and AdaBoost have achieved remarkable results in many fields.

In recent years, with the enhancement of GPU-based computation power and data volume, Convolutional Neural Networks (CNNs) have greatly promoted the development of visual recognition, and have achieved state-of-the-art in many areas of recognition tasks such as object detection, scene classification, and image recognition. The principle of CNN has been developed for nearly 40 years since it was proposed in the 1980s. The most urgent research shows that the general features extracted by CNN have strong generalization ability. The features extracted in a certain field can be transferred and well generalized to other visual recognition tasks. Shallow features such as SIFT can easily cause the semantic gap problem in LSVGL, that is, features extracted from different locations may be very similar. Deep feature vectors extracted using CNN networks with deep layers stacking can effectively alleviate this problem.

Triplet loss was initially proposed in \cite{facenetSchroffKP15} to learn robust facial embedding and online triple mining. Different from the traditional softmax loss, triplet loss combines the anchor sample, positive sample, negative sample to form a tuple and helps the network to learn a kind of face embedding, so that minimizing the distance between the anchor face and positive samples, and maximizing the distance between the anchor and negative samples in the embedding space. Owing to the success of triplet loss, the loss derived from triplet loss has also been widely concerned \cite{quadruplet17, trihard17}. TriHard loss \cite{trihard17} selects the negative sample closest to the training sample as hard sample and uses hard sample to replace randomly selected negative sample in triplet loss. Quadruplet loss \cite{quadruplet17} is trained with two negative sample constraints, and the calculation of relative distance in triplet loss is modified to absolute distance. Margin Sample Mining Loss (MSML) \cite{mlms17} takes advantage of TriHard loss and quadruplet loss. After great success in the field of face recognition, researchers began to try to apply triplet loss to LSVGL and also achieved enormous great success \cite{arandjelovic2016netvlad, zhu2018attention, liuliu19}.

Unlike most of the work dedicated to proposing network architectures, we try to use existing network models to analyze the objective situation of the problem and explore whether the problem can be optimized from the perspective of loss. We believe that the biggest drawback of triplet loss and quadruplet loss is that they focus too much on maximizing the distance between anchor samples and negative samples, while ignoring the strong constraints between anchor samples and positive samples. Generally, due to the time difference between the anchor samples and the positive samples and the arbitrariness of the spatial position, although the image at the same location inevitably has a perspective deviation, that is, we can reach a consensus: an anchor sample consists of multiple positive samples in the location partial image composition, which is a significant difference between LSVGL and face recognition. Regardless of triple loss or quadruple loss, only one positive sample is considered, and we believe that this loss is not completely suitable for the LSVGL scenario. We put forward the hypothesis that adding all potential positive samples to the loss to train together will have a gain effect on LSVGL.

In this paper, we propose a new loss called QUInTuplet loss (QUITLoss), which attempts to maintain the advantages of the triple loss and adds all potential positive samples to the loss. This is because there is usually a difference in perspective between the anchor sample and the positive sample, but this problem does not exist in the field of facial recognition. We start with this general common sense and theoretically analyze the rationality of QUITLoss. After a careful comparison, we concluded that our QUITLoss could make up for the shortcomings of previous methods. Fig. \ref{fig:cmp_with_other_loss} shows the main idea of QUITLoss.

We verified the role of QUITLoss in experiments. We corrected the loss in the previous method with QUITLoss and performed ablation experiments. The experimental results verify the effectiveness of our method for the LSVGL task.

In the following, we overview the main contents of our method and summarize the contributions:

1). We embed several metric learning losses into state-of-the-art LSVGL framework and make up for the gap in the current LSVGL solutions that only use triplet loss and do not compare with other metric learning methods.

2). We propose a new metric learning loss strategy to compensate for the precision wastage of the existing methods by considering all potential positive samples to optimize the model which outperforms other metric learning losses on LSVGL.

The paper is organized as bellow: related works with more details are presented in section \ref{sec:related_work}. Then we give a briefly overview of the proposed method and a detailed introduction in section \ref{sec:overview} and \ref{sec:method}, respectively. Datasets and experiments are presented in section \ref{sec:experiments}. Conclusions and outlook are presented in section \ref{sec:conclusion}.

\section{Related Work}\label{sec:related_work}
In this section, we briefly review the methods related to LSVGL. First, we summarize the development of CNNs, and then list the research of metric learning and LSVGL.

\textbf{Deep CNNs.}
Many popular deep networks including AlexNet, VGGNet, GoogleNet, and ResNet networks have been proposed in recent years. Many works have proved that the VGG network pre-trained on the ImageNet image classification dataset is better than other baseline models in the LSVGL task. Therefore, in this paper, we also choose VGG16 as the benchmark network like most other papers. VGGNet is a deep convolutional neural network developed by Oxford University's Visual Geometry Group and researchers from Google DeepMind. Its main contribution is to show that the depth of the network is a key part of the excellent performance of the algorithm. Based on the VGGNet model, many excellent models have been improved, such as the use of new optimization algorithms and multi-model fusion. Furthermore, recognition of salient objects in a scene is helpful to extract the fixed and unchanging structural features at a certain location, the image representation extracted by VGGNet are used in many image recognition tasks\cite{multi-salient, hybrid-graph,contour-transfer18, webly-salient-20}. Yang \cite{multi-salient} proposed a multi-scale cascade network to identify the most visually conspicuous objects in an image and Luo \cite{hybrid-graph} presented a hybrid graph neural network to interweaving the multi-scale features for crowd density as well as its auxiliary task together, to better capture and consolidate multi-scale high-level context information for object skeleton detection, \cite{contour-transfer18} automatically converted an existing deep contour detection model into a salient object detection model without using any manual salient object masks by grafting a new branch onto a well-trained contour detection network. Luo \cite{webly-salient-20} proposed a webly-supervised learning method for salient object detection with no pixel-wise annotations.

\textbf{Deep metric learning.}
Metric learning is a distance function that measures similarity: similar objects are close to each other, and dissimilar ones are far away. While deep learning is still asleep, the Mahalanobis distance in Euclidean space solves the judgment of distance in traditional metric learning, Cross-view Quadratic Discriminant Analysis (XQDA) and Keep It Simple and Straightforward Metric learning (KISSME) were both classic metric learning methods in the part. Deep metric learning currently mainly uses the network to extract embedding, and then uses L2-distance to measure distances in the embedding space.

In recent years, a large amount of work has applied metric learning from person re-identification to LSVGL \cite{arandjelovic2016netvlad, liuliu19, kim2017learned, zhu2018attention}. Two images at the same location are defined as a positive pair, while two images at different locations are defined as a negative pair. Triplet loss has been used with great success in LSVGL \cite{arandjelovic2016netvlad, zhu2018attention}. As an upgraded version of triplet loss, the triplet hard sample mining loss (TriHard) selects the hardest negative sample, even if the anchor is farthest from the nearest (hardest) negative sample. Quadruplet loss extends TriHard loss, retains the idea of 'hardest', and simultaneously selects one more negative sample and adds the absolute distance between the negative samples to the training. Margin Sample Mining Loss (MSML) is a combination of the advantages of the above metric learning methods. It selects the most similar positive samples (the hardest) and the most similar negative samples (the hardest) to optimize the quadruplet loss. Unfortunately, TriHard, quadruplet, MSML did not verify performance in LSVGL, in this paper we make up for this regret.

\textbf{Deep LSVGL related methods.}
LSVGL remains a challenging task in urban environments due to the repetitive structures, illumination condition, and differences in perspective. While there have been many improvements in designing better LSVGL systems \cite{arandjelovic2016netvlad, zhu2018attention, kim2017learned, liuliu19, chen2017deep}. Arandjelovic \cite{arandjelovic2016netvlad} performed learning for LSVGL and proposed the NetVLAD representations, which significantly outperform the previous local-feature-based representations on several benchmarks. Zhu \cite{zhu2018attention} proposed an attention mechanism to solve the LSVGL task. However, none of these works actually training the CNNs for the task considered anchors are almost assigned into multiple images, and these images should be used as positive samples to network learning.

\begin{figure*}
	\centering
	\includegraphics[scale=0.98]{graph/network.pdf}
	\caption{A schematic diagram of metric learning in deep neural networks.}
	\label{fig:metric_learning_diagram}
\end{figure*}

\section{An Overview of Our Method}\label{sec:overview}
We briefly introduce the overview of the proposed method. We maintain the view that the current methods to solve the LSVGL tasks have two main shortcomings: (1). The metric learning methods used for LSVGL is too limited. Almost the current state-of-the-art methods use triplet loss, however many excellent metric learning methods have already been proposed, such as quadruplet loss, TriHard, MSML methods, etc. and (2). We can summarize the LSVGL problem into the task of image retrieval in urban environments. Although it is the same location, it is often taken during retrieval compared with the image used in training the geometric database, the image of this place will appear infection factors such as lighting, perspective, occlusion, and the dynamic scene that everyone knows. Many works are also working on this kind of problem and have achieved quite good results, but there is almost no work to focus on during the training phase, anchor exists in multiple positive samples, and the most commonly used triplet loss only picks the nearest positive sample to learn, which undoubtedly reduces the constraints of positive examples on anchor. Fig. \ref{fig:metric_learning_diagram} shows a schematic diagram of metric learning in deep neural networks.

To address the first shortcoming, we first summarize the current mainstream metric learning methods, and use the current state-of-the-art LSVGL network NetVLAD to compare the impact of multiple metric learning methods on geo-localization performance.

To address the second shortcoming, we considered the characteristics of LVSGL, that is, the problem caused by the arbitrariness of the perspective that anchors may exist in multiple positive samples, we propose a new metric learning method called QUITLoss. QUITLoss is originally designed to expand only one positive sample into multiple positive samples based on triplet loss. We use QUITLoss to correct each metric learning method and compare the performance of other metric learning methods after correction.

\section{Proposed Quintuplet Loss}\label{sec:method}
First, we introduce several metric learning methods, including triplet loss, TriHard loss, quadruplet loss, and MSML loss. And then we formally introduce the proposed QUITLoss.

\subsection{Background}
To extract image depth embedding, metric learning aims to learn a mapping: $g(x): R^{W \times H} \rightarrow R^D$ which maps semantically similar images from the raw data manifold in $R^{W \times H}$ onto metrically close points in $R^D$ \cite{trihard17}, where $W$ and $H$ are the width and height of input image respectively and $D$ is the dimension of the final embedding. The deep metric learning aims to find the mapping through minimizing the metric loss of training data. Then a metric function $L(x,y): R^D \times R^D \rightarrow R$ is designed to measure similarity in the embedding space.

\textbf{Triplet Loss}. Triplet loss generates more distinctive features than softmax loss \cite{facenetSchroffKP15}. A triplet in triplet loss contains three different images $\{I^{(a)}, I^{(p)}, I^{(n)}\}$, while $I^{(a)}$ is the anchor image and $I^{(p)}$ and $I^{(n)}$ represent the positive image in the same place with $I^{(a)}$ and negative image in another different place with $I^{(a)}$ respectively. A $W$-by-$H$ image generates a $D$-dimensional feature after network mapping, and a triplet of features $\{g^{(a)}, g^{(p)}$, $g^{(n)}\}$ would be embedded into the triplet loss formulation as bellow:
\begin{equation}
L_{tri}=h(\overbrace{||g^{(a)}-g^{(p)}||_2}^{to\ shorten} - \overbrace{||g^{(a)}-g^{(n)}||_2}^{to\ enlarge} + \alpha)
\label{eq:triplet}
\end{equation}

where $h(x)=max(x,0)$ and $\alpha$ is the margin to distinguish between the positive samples and the negative ones. The first term $||g^{(a)}-g^{(p)}||_2$ shortens the distance between $I^{(a)}$ and $I^{(p)}$, and the second term $||g^{(a)}-g^{(n)}||_2$ enlarges the distance between $I^{(a)}$ and $I^{(n)}$.

\textbf{Quadruplet loss.} Compared with the triplet loss, quadruplet loss \cite{quadruplet17} contains one more negative. A quadruplet selects four different images $\{I^{(a)}, I^{(p)}, I^{(n_1)}, I^{(n_2)}\}$, where $I^{(a)}$ and $I^{(p)}$ have the same meaning as triplet loss, $I^{(n_1)}$ and $I^{(n_2)}$ are negative samples. The quadruplet loss is formulates as bellow:
\begin{equation}
\begin{aligned}
L_{quad}&=h(\overbrace{||g^{(a)}-g^{(p)}||_2 - ||g^{(a)}-g^{(n_1)}||_2 + \alpha}^{relative\ distance}) \\
&+ h(\overbrace{||g^{(a)}-g^{(p)}||_2 - ||g^{(n_1)}-g^{(n_2)}||_2 + \beta}^{absolute\ distance})
\end{aligned}
\label{eq:quadruplet}
\end{equation}
where $\alpha$ and $\beta$ are the margin in two terms.  The first term is a triplet loss, which focuses on the relative distance between positive pair and negative pair, the second term adds the absolute distance by considering different negative images.

\textbf{TriHard Loss}. As the possible number of quadruplets grows rapidly as the data volume gets larger, loss Eq. (\ref{eq:quadruplet}) can be further optimized. In order to relieve this,  triplet loss with hard sample mining computes a batch of samples together. For each image, the TriHard loss picks the most similar sample with different place $I^{(h,n)}=argmin_n (||g^{(a)}-g^{(n)}||_2)$ to generate a triplet $\{I^{(a)}, I^{(p)}, I^{(h,n)}\}$, then estimate loss by Eq. (\ref{eq:triplet}).

\textbf{Margin Sample Mining Loss (MSML).} MSML picks the most dissimilar positive pairs and most similar negative pair in the whole batch, as:
\begin{equation}
L_{msml}=h(||g^{(a)}-g^{(h,p)}||_2-||g^{(h,n_{1})}-g^{(h,n_{2})}||_2+\alpha)
\end{equation}

where $g^{(h,n_{1})}$ and $g^{(h,n_{1})}$ can share the same identity with $I^{(a)}$ or not.

\subsection{QUITLoss}

\subsubsection{Definition}
We apply a new positive example mining strategy for LSVGL named QUITLoss. Different from the loss introduced above, QUITLoss selects multiple positive samples to integrate into the final loss calculation:
\begin{equation}
L_{quit}=\sum_{i=1}^k h(||g^{(a)}-g^{(p_{(i)})}||_2 - \mathcal{P} + \alpha)
\end{equation}
where $g^{(p_{(i)})}$ is the feature of the $i$-th positive of the anchor image $I^{(a)}$, and $\mathcal{P}$ represents multiple negative sample selection methods, which will be explained in detail below. Fig. \ref{fig:quit-loss-graph} shows the learning process.

\begin{figure}[h]
	\centering
	\includegraphics[scale=0.45]{graph/quit-loss.pdf}
	\caption{QUITLoss keeps the same learning strategy with triplet loss but selects multiple positives.}
	\label{fig:quit-loss-graph}
\end{figure}

We design QUITLoss as a general concept, not a specific mathematical form. Different negative sample selection strategies can form specific expressions.
When $\mathcal{P}=||g^{(a)}-g^{(h, n)}||_2$, the final QUITLoss can be formed after fusion with the TriHard loss:
\begin{equation}
L_{quit\_trihard}=\sum_{i=1}^k h(||g^{(a)}-g^{(p_{(i)})}||_2 - ||g^{(a)}-g^{(h, n)}||_2 + \alpha)
\label{eq:trihard-quit}
\end{equation}
where $g^{(h, n)}$ is the feature of the nearest negative sample of the anchor image $I^{(a)}$.

We set $\mathcal{P}=||g^{(a)}-g^{(n_1)}||_2$ and rewrite the mathematical expression after the fusion of QUITLoss and quadruplet loss:
\begin{equation}
\begin{aligned}
L_{quit\_quad}&=\sum_{i=1}^k h(||g^{(a)}-g^{(p_{(i)})}||_2 - ||g^{(a)}-g^{(n_1)}||_2 + \alpha) \\&+ \sum_{i=1}^k h(||g^{(a)}-g^{(p_{(i)})}||_2 - ||g^{(n_1)}-g^{(n_2)}||_2 + \beta)
\end{aligned}
\end{equation}

\subsubsection{Gradient}
In addition, we also provide the gradient information of TriHard loss modified by QUITloss. The gradient of the loss function $L_{quit\_trihard}$ for each variable is calculated as follows:
\begin{equation}
\begin{aligned}
\frac{\partial L}{\partial {g^{(a)}}} &= \sum_{i=1}^k 2(g^{(a)}-g^{(p_{(i)})} - (g^{(a)}-g^{(h,n)})) \\&=\sum_{i=1}^k 2(g^{(h,n)} - g^{(p_{(i)})})
\end{aligned}
\end{equation}

\begin{equation}
\begin{aligned}
\frac{\partial L}{\partial g^{(p_{(i)})}} &= 2(g^{(a)}-g^{(p_{(i)})})(-1) =2(g^{(p_{(i)})}-g^{(a)}), \\& i=1, 2, ..., k
\end{aligned}
\end{equation}

\begin{equation}
\frac{\partial L}{\partial g^{(h, n)}} = -2(g^{(a)}-g^{(h, n)})(-1)=2(g^{(a)}-g^{(h, n)})
\end{equation}
Other loss forms can be deduced in the same way.

In summary, compares with other metric learning losses, our QUITLoss has the following advantages. First, QUITLoss makes up for the problem that the metric learning methods in the current LSVGL task do not consider the anchor images assigned to multiple positive samples caused by the difference in perspective. Second, QUITLoss does not completely discard the existing methods, but provide a general idea, which can be extended to the existing metric learning loss. Finally, QUITLoss is easy to implement.

\begin{table*}[t]
\renewcommand\arraystretch{1.5}
  \centering
  \caption{Comparison of state-of-the-art deep metric learning loss. All the results are from the 512-D representations based on VGG-16 architecture.}
  \label{tab:performance_comparison}
  \setlength{\tabcolsep}{4mm}{
    \begin{tabular}{|c|c|c|c|c|c|c|c|}
    \hline

    \multirow{2}{*}{Base model} &
    \multirow{2}{*}{Methods} &
    \multicolumn{3}{c|}{Pitts250k-test} & \multicolumn{3}{c|}{Tokyo 24/7}\cr\cline{3-8}
    & &Recall@1&Recall@5&Recall@10&Recall@1&Recall@5&Recall@10\cr
	\hline
     \multirow{8}{*}{VGG16}& NetVLAD\cite{arandjelovic2016netvlad} &0.8066&0.9088&0.9306 & 0.6000 & 0.7365& 0.7905\cr\cline{2-8}
    & triplet& 0.8169 &  0.9112 &  0.9329 &0.6286 & \textbf{0.7714}& 0.8032\cr
    & triplet + QUITLoss (ours)& \textbf{0.8226} & \textbf{0.9118} & \textbf{0.9341} & \textbf{0.6349}  & 0.7587 & \textbf{0.8159}\cr\cline{2-8}
    & quad& 0.7950 &  0.8966 &  0.9209 & \textbf{0.6508}& 0.7841& 0.8317\cr
    & quad+QUITLoss (ours)& \textbf{0.8030} & \textbf{0.9010} &  \textbf{0.9252} & 0.6476 & \textbf{0.7905} & \textbf{0.8413} \cr\cline{2-8}
    & trihard& 0.8333 &0.9203&0.9313& 0.6730& 0.7905 & 0.8381 \cr
    & trihard + QUITLoss (ours)& \textbf{0.8403} & \textbf{0.9216} & \textbf{0.9401} & \textbf{0.6952} & \textbf{0.8159} & \textbf{0.8476} \cr\cline{2-8}
    & MSML& 0.7431 & 0.8639 & 0.8958 & 0.6032& 0.7492& 0.8032\cr\cline{2-8}
    \hline
    \end{tabular}}
\end{table*}

\section{Experiments}\label{sec:experiments}

\subsection{Datasets and Evaluation Methodology}
We report results on two publicly available datasets. \textbf{Pittsburgh (Pitts250k)} \cite{torii2013visual}  is a popular LSVGL benchmark and contains 250k database images downloaded from Google Street View and 24k queries generated from Street View but taken at different times, years apart.  Similar to \cite{arandjelovic2016netvlad}, we divide Pitts250k into three roughly equal parts for training, validation, and testing, namely Pitts250k-train, Pitts250k-valid, and Pitts250k-test respectively. Each part contains approximately 8k training images and 8k queries. To facilitate faster training, we use a smaller subset (Pitts30k) which contains 10k training images for training.

\textbf{Tokyo 24/7} \cite{torii201524} is another extremely challenging dataset which contains 76k database images and 315 queries.  The database images were only taken at daytime from Google Street View as described above, but the queries were taken under extreme illumination conditions.

\textbf{Evaluation metric.} We follow the standard LSVGL evaluation procedure \cite{retrieval12, torii201524, torii2013visual, arandjelovic2016netvlad}. The localization is deemed correct if at least one of the top N retrieved database images are within $d$=25 meters from the query. $Recall@N$ is the percentage of correctly localized queries for different values of N.

\begin{figure*}[t]
	\centering
	\includegraphics[scale=0.56]{graph/hardest-positive-mining.pdf}
	\caption{Schematic diagram of positive sample selection in QUITLoss. (1). Image (a) is the anchor image, and image (b)-(i) share the same geometric location with the anchor (a). In our experiment, image (b)-(i) are considered as potential positives, and other metric learning losses only choose the nearest image (b) as the positive,  but the brick red church exists in multiple images (b) and (c), we choose two positives to our QUITLoss. (2). If our QUITLoss keeps the hardest positive selection strategy, It is entirely possible to pick out (d)-(i), those are certainly not positive samples from the perspective of human supervision.}
	\label{fig:not_hardest_positive}
\end{figure*}

\subsection{Implementation Details}
Each image is resized into 224 $\times$ 224 pixels with data augmentation. The data augmentation includes randomly flipping, zooming, blurring, and contrast enhancement. We use the VGG16 pre-trained model as our base model, the margin $\alpha$ and $\beta$ are set to 0.3 and 0.2 respectively. During gradient descent, the SGD optimizer is used and the initial learning rate is set to $10^{-4}$ and drops 0.5 times every 5 epochs. The training is terminated if the performance is not improved for 10 consecutive epochs.

\begin{table}[t]
\renewcommand\arraystretch{1.2}

  \caption{Experimental results on the influence of different hyperparameter $k$ on the results of geographical location.}
  \label{tab:learning_k}
  \setlength{\tabcolsep}{6mm}{
    \begin{tabular}{|c|c|c|c|}
    \hline
    \multirow{2}{*}{$k$}&
    \multicolumn{3}{c|}{Pitts250k-test}\cr\cline{2-4}
    &Recall@1&Recall@5&Recall@10\cr
    \hline
    1 & 80.66 &  90.88 & 93.06\cr\hline
    2 & \bf{82.34} & \bf{90.94} & \bf{93.41} \cr\hline
    3 & 81.64 & 90.75 & 93.18 \cr\hline
    4 & 74.79 & 86.28 & 89.93 \cr\hline
    \end{tabular}}
\end{table}

\textbf{Learning the hyperparameter k.} To begin with our experiments, we first study the impacts of the hyperparameter k of our method. We set different parameters to verify the impact of different numbers of positive samples on the Pittks250k-test dataset. The network parameter is the baseline network when $k$=1. The experiment found that the network does not converge when $k$ is greater than 4, and the localization results are best when $k$ = 2, which also meets the theoretical requirements. Table \ref{tab:learning_k} shows the details. The same observation can be seen on the Tokyo 24/7 dataset. The experimental results also verify our cognition of LSVGL, that is, due to the arbitrariness of the two shooting angles, the anchor image is almost allocated to two positive samples. Thus, we empirically select k=2 in all our following experiments.

\subsection{Results and Discussion of Different Losses}
We conduct experiments with different losses and provide the results to illustrate the effectiveness of the proposed QUITLoss. The PyTorch version of the open-source code \cite{arandjelovic2016netvlad} is used in our experiment and four metric learning losses, triplet loss, quadruplet loss, TriHard loss, MSML are implemented to compare with the origin triplet loss.

\textbf{Results of Different losses.}
In our experiments, we compared the impact of different metric learning methods on geo-localization performance. The experimental results show that TriHard has advantages over other metric learning methods because it considers using the  nearest (hardest) negative samples to constrain anchors in the metric space. On the Pitts250k-test dataset and under the recall @ 1 indicator, the TriHard method is 2.67\% higher than the state-of-the-art method NetVLAD and is 1.64\% higher than the triplet method we implemented. Quadruplet loss does not significantly help the localization results, because although quadruplet loss uses the absolute distance between two negative samples to constrain the anchor, it still does not have a strong constraint on the anchor with the nearest negative sample. It is worth mentioning that the training of MSML has not been successful, and the performance of MSML on the localization results is greatly reduced. The reason why the hardest (farthest) positive sample cannot be selected will be analyzed later. On the Tokto247 dataset, we also have the same conclusion, see Table \ref{tab:performance_comparison} for details.

\textbf{Discussion of QUITLoss.}
Every metric learning method we use can be corrected with QUITLoss. Experimental results show that the corrected metric learning methods can improve the performance of geo-localization. Similarly, the corrected TriHard loss has better localization results than other metric learning methods. After correction, the triplet loss, quadruplet loss, and TriHard are improved by 0.57\%, 0.80\%, and 0.70\%, respectively, compared to the performance of the methods before correction. The performance reported on Tokyo247 is shown in Table \ref{tab:performance_comparison}.

\textbf{Comparison with state-of-the-art.}
NetVLAD uses Pytorch's built-in functions to calculate triplet loss. For an anchor, multiple negative samples are considered. Taking TriHard as an example, the performance of pure TriHard loss is 2.67\% higher than NetVLAD, and TriHard loss is 3.37\% higher after QUITLoss correction.

\textbf{Why can we not choose the hardest positive?} In LSVGL system, we can't choose the most difficult positive samples, because in panoramic images, we take those images which are close to the anchor image's geographical location as potential positive samples, but there are only two positive samples (we have proved why there are only two) that really belong to the anchor, so in LSVGL only the nearest (not the most difficult) can be selected. As Table \ref{tab:performance_comparison} shows, MSML loss significantly degrades performance, and Fig. \ref{fig:not_hardest_positive} shows that the most difficult is not the anchor's positive sample, and there is no semantic similarity between the most difficult and the anchor. For example, the images that are taken from the opposite perspective proves this conclusion.

\section{Conclusion}\label{sec:conclusion}
In this paper, we propose a new metric learning loss with multiple positive constraints named quintuplet loss for Large Scale Visual Geo-Localization (LSVGL). In our method, we calculate a distance matrix and choose the smallest two distance of positive pairs to form a general strategy. In this way, quintuplet loss used all the positives to train the model. In addition, we also carefully analyze why selecting two positive samples is the best strategy.

We use the pre-trained VGG-16 model as our baseline and the state-of-the-art model NetVLAD to do contrast experiments with different metric learning losses on several benchmarks datasets, including Pitts250K, Tokyo24/7. The results show the metric learning loss corrected by our quintuplet loss gets better performance and learns a finer metric in feature embedding space.

\ifCLASSOPTIONcompsoc

  \section*{Acknowledgments}
\else

\fi
The authors gratefully acknowledge the financial supports by the National Science Foundation of China under grant numbers U1964203, useful suggestions given by Xin Li, Fan Yang, Ao Luo, Hao Guo are also acknowledged.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{unsrt}


\section{Introduction}

The classical theory of Optimal Transport deals with the problem of efficiently transporting mass from a probability distribution into a target one. In the last thirty years, great advances in the understanding of the underlying theory have been achieved \cite{ags, cuturi,sant}. However, only recently these techniques are starting to be applied in order to solve computational problems in a great variety of fields, with logistic problems \cite{bernot, mass, cities, equilibrium}, crowd dynamics \cite{macroscopic, handlingcongestion}, image processing \cite{hmp,mrss,mbmoj,histogram,papadakis, activecontours,colourmapping}, inverse problems \cite{bf,kr} and machine learning \cite{wassersteingan, wassersteinloss, learninggenerative, wassersteintraining, sae, wae} being a few examples.

In this paper we focus on the so-called Benamou--Brenier formula,
which provides an equivalent dynamic formulation of the classical Monge--Kantorovich transport problem \cite{kantorovich}.
Introduced by Benamou and Brenier in \cite{bb}, such formula allows to compute an optimal transport between two probability measures $\rho_0$ and $\rho_1$ on a closed bounded domain $\overline{\Om} \subset \mathbb{R}^d$ through the minimization of the kinetic energy
\begin{equation} \label{intro:bb}
\frac{1}{2}\int_0^1\int_{\overline{\Om}} |v_t(x)|^2 d\rho_t(x)\,,
\end{equation}
among all the pairs $(\rho_t,v_t)$, where $t \mapsto \rho_t$ is a curve of probability measures on $\overline{\Om}$, $v_t \colon \overline{\Om} \to \mathbb{R}^d$ is a time-dependent vector field and the pair $(\rho_t,v_t)$ satisfies distributionally the continuity equation
\begin{equation} \label{intro:cont}
\partial_t \rho_t + \div(\rho_t v_t) = 0 \quad \text{ subjected to } \quad \rho_{t=0}=\rho_0 , \,\, \rho_{t=1}=\rho_1 \,.
\end{equation}
The interest around the Benamou--Brenier formulation is motivated by its remarkable properties.  First, it allows to compute an optimal transport in an efficient way \cite{bb} by means of a convex reformulation of \eqref{intro:bb}, by introducing the momentum $m_t=\rho_tv_t$. More precisely, denoting by  $X:=(0,1) \times \overline{\Om}$ the time-space cylinder, the Benamou--Brenier energy \eqref{intro:bb} can be equivalently defined as a convex functional on the space of bounded Borel measures $\mathcal{M}:=\mathcal{M}(X) \times \mathcal{M}(X;\mathbb{R}^d)$ by setting
\begin{equation} \label{intro:BB_convex}
B (\rho,m):= \frac{1}{2}\int_{X} \left| \frac{dm}{d\rho}(t,x)\right|^2 \, d \rho(t,x) \,,
\end{equation}
whenever $(\rho,m)\in \mathcal{M}$ are such that $\rho \geq 0$, $m \ll \rho$, and $B:=+\infty$ otherwise.
With this change of variables, the continuity equation at \eqref{intro:cont} assumes the form
\begin{equation}\label{eq:contmeasures}
\partial_t \rho + \div m = 0 \,.
\end{equation}
In addition, the  dynamic nature of the Benamou--Brenier reformulation of optimal transport is at the core of many recent developments in the fields of PDEs, optimal transport and inverse problems. Indeed, the dynamic formulation allows to endow the space of probability measures with a differentiable structure \cite{ags,sant}, making possible the characterization of differential equations as gradient flows in spaces of measures \cite{agsinventiones,ottokinderleher2,otto1} or the derivation of sharp inequalities \cite{cordero,maggi,otto_villani}.  Moreover, it motivated recent developments in unbalanced optimal transport theory \cite{chizat2,chizat,kmv,liero}, i.e., when the marginals are arbitrary positive measures. Finally, as the Benamou--Brenier energy provides a description of the optimal flow of the transported mass at each time $t$, which is a valuable information in applications, it was recently employed as a regularizer for variational inverse problems \cite{inprep,inprep2,bf, hmp,mrss,schmitzerwirth2020}.

\smallskip

The goal of this paper is to characterize the extremal points of the unit ball of the Benamou--Brenier energy $B$ at \eqref{intro:BB_convex}, and of a coercive version of it, which is obtained by adding the total variation of $\rho$ to $B$. Both functionals are constrained via the continuity equation \eqref{intro:cont}.  Precisely, we introduce the functional
\begin{equation}\label{intro:J}
J_{\alpha,\beta}(\rho,m):=   \beta B(\rho,m) + \alpha \|\rho\|_{\mathcal{M}(X)}  \quad \text{ subjected to } \quad \partial_t \rho+ \div m=0 \,,
\end{equation}
defined for all $(\rho,m) \in \mathcal{M}$ and $\alpha \geq 0$, $\beta>0$. We then characterize the extremal points of the subset of $\mathcal{M}$ defined by
\[
C_{\alpha,\beta}:=\{(\rho,m) \in \mathcal{M} \, \colon \, J_{\alpha,\beta} (\rho,m)  \leq 1\}\,.
\]
We emphasize that we do not enforce initial conditions to the continuity equation in \eqref{intro:J}. To be more specific, we prove the following result (see Theorem \ref{thm:extremal}).
\begin{thmnonumber}
Let $\alpha \geq 0$, $\beta>0$. The extremal points of the set $C_{\alpha,\beta}$
  are exactly given by the zero measure $(0,0)$ and the pairs of measures $(\rho,m)$ such that
  \[
  \rho=a_{\gamma} \, dt \otimes  \delta_{\gamma(t)},\, \,\,\, \, m=\dot{\gamma}(t) a_\gamma\, dt \otimes  \delta_{\gamma(t)} , \,\,\,\,\, a_\gamma = \left(\frac\beta2 \int_0^1 |\dot\gamma(t)|^2 \, dt + \alpha\right)^{-1}\,,
  \]
  where $\gamma: [0,1] \to \overline{\Om}$ is an absolutely continuous curve with weak derivative $\dot\gamma \in L^2$,
  and such that $a_{\gamma}<+\infty$. If $\alpha = 0$ the condition $a_{\gamma}<+\infty$ is satisfied if and only if $\gamma$ is not constant.
\end{thmnonumber}
We therefore show that the extremal points of the set $C_{\alpha,\beta}$ are pairs of measures $(\rho,m)$, with $\rho$ concentrated on some absolutely continuous curve $\gamma$ in $\overline{\Om}$, and the density of $m$ with respect to $\rho$ is given by $\dot \gamma$ . Notice that such conditions are equivalent to the existence of a measurable field $v \colon X \to \mathbb{R}^d$ such that
\begin{equation}\label{ode}
\dot{\gamma}(t)=v(t,\gamma(t)) \quad \mbox{ for a.e.} \quad t \in (0,1)\,,
\end{equation}
thus showing that $\gamma$ is a characteristic associated to the continuity equation at \eqref{intro:cont} with respect to the field $v$.
We prove the above Theorem in Section \ref{sec:characterizationextremal}, with the aid of a probabilistic version of the superposition principle for positive measure solutions to the continuity equation \eqref{intro:cont} on the domain $(0,1) \times \overline{\Om}$ (see Theorem \ref{thm:lifting}). We mention that the ideas behind such superposition principle are not new, and they were  originally introduced in \cite{ambrosioinventiones} for positive measures on $(0,1) \times \mathbb{R}^d$ (see also \cite{ags, bernard, trevisan}). The result of Theorem \ref{thm:lifting} allows to decompose any measure solution $(\rho,m)$ of the continuity equation \eqref{eq:contmeasures} with bounded Benamou--Brenier energy, as superposition of measures concentrated on absolutely continuous characteristics of \eqref{eq:contmeasures}, that is, curves solving \eqref{ode} with $v=dm/d\rho$. As a consequence, we show any pair of measures that is not of such a form can be written as a proper convex combination of elements of $C_{\alpha,\beta}$ and thus it is not an extremal point.  The opposite inclusion follows from the convexity of the energy and the properties of the continuity equation.

\smallskip

The interest in characterizing extremal points of the Benamou--Brenier energy is not only theoretical. It has been recently shown in \cite{bc} and \cite{chambolle} that in the context of variational inverse problems with finite-dimensional data, the structure of sparse solutions is linked to the extremal points of the unit ball of the regularizer.
In the classical theory of variational inverse problems one aims to solve
\begin{equation}\label{eq:inverseproblemintro}
\min_{u\in \mathcal{U}} \ R(u) \quad \text{ subjected to } \quad  Au = y\,,
\end{equation}
where $\mathcal{U}$ is the target space, $R$ is a convex regularizer, $A$ is a linear observation operator mapping to a finite-dimensional space and $y$ is the observation.
It has been empirically observed that the presence of the regularizer $R$ is promoting the existence of sparse solutions, namely minimizers that can be represented as a finite linear combination of simpler atoms.
While this effect has been well-understood in the case when $\mathcal{U}$ is finite dimensional, the infinite-dimensional case has been only recently addressed \cite{chambolle, bc, sparsetik, flinthweiss, unser3, unser2, unsersplines}. In particular, in \cite{chambolle,bc}, it has been shown that, under suitable assumptions on $R$ and $A$, there exists a minimizer of \eqref{eq:inverseproblemintro} that can be represented as a finite linear combination of extremal points of the unit ball of $R$; namely the atoms forming a sparse solution are the extremal points of the ball of the regularizer.

In view of the above discussion, in Section \ref{sec:application} we apply our characterization of the extremal points of the energy $J_{\alpha,\beta}$ at \eqref{intro:J} to understand the structure of sparse solutions for inverse problems with such transport energy acting as regularizer. We mention that the analysis is carried out for the case $\alpha > 0$, as the functional $J_{0,\beta}$, corresponding to the rescaled Benamou--Brenier energy, lacks of compactness properties (see Remark \ref{rem:compact}). We verify that the assumptions needed to apply the representation theorems in \cite{bc} and \cite{chambolle} are satisfied by $J_{\alpha,\beta}$, and consequently we deduce the existence of a minimizer that is given by a finite linear combination of measures concentrated on absolutely continuous curves in $\overline{\Om}$ (see Theorem \ref{thm:sparserepresentation}). As a specific application of Theorem \ref{thm:sparserepresentation} we consider the setting introduced in \cite{bf}, where the regularizer $J_{\alpha,\beta}$ is coupled with a fidelity term that penalizes the distance between the unknown measure $\rho_t$ computed at $t_1,\ldots,t_N \in (0,1)$, and the observation at such times (see Section \ref{subsec:dynamic}). This setting is relevant for applications, such as variational reconstruction in undersampled dynamic MRI. Employing the previous results we are able to prove the existence of a sparse solution represented with a finite linear combination of measures concentrated on absolutely continuous curves in $\overline{\Om}$ (see Corollary \ref{corollary}).

\smallskip

To conclude, we mention that characterizing the extremal points for a given regularizer has important consequences in devising algorithms able to compute a sparse solution. Notable examples have been proposed for the total variation regularizer in the space of measures \cite{boyd, brediesp} using so-called \emph{generalized conditional gradient methods} (or Frank--Wolfe-type algorithms \cite{frankwolfe}). Inspired by the previous methods, and building on the theoretical results obtained in the present paper, we plan to develop numerical algorithms to compute sparse solutions of dynamic inverse problems with the optimal transport energy $J_{\alpha,\beta}$ as a regularizer \cite{inprep,inprep2}, effectively providing a numerical counterpart to the theoretical framework established in \cite{bf}.
Finally, we remark that similar results to the ones presented in this paper can be obtained for unbalanced optimal transport energies. This has been recently achieved in \cite{superpositioninhomogeneous}, by introducing a novel superposition principle for measure solutions to the inhomogeneous continuity equation.

\section{Mathematical setting and preliminaries}

In this section we give the basic notions about the continuity equation, the Benamou--Brenier energy, and its coercive version $J_{\alpha,\beta}$  anticipated in the introduction. We refer to \cite{ ags, bb, bf} for a more detailed overview. For measure theoretical notions, we refer to the definitions in \cite{afp}.

Given a metric space $Y$ we will denote by $\mathcal{M}(Y)$ (resp.~$\mathcal{M}(Y ; \mathbb{R}^d)$)  the space of bounded Borel measures (resp.~bounded vector Borel measures) on $Y$. Similarly,  $\mathcal{M}^+(Y)$ and $\mathcal{P}(Y)$ denote the set of bounded positive Borel measures and Borel probability measures on $Y$, respectively.
Let \(\Omega \subset \mathbb{R}^d\) be an open, bounded domain with \(d\in \mathbb{N}, d \geq 1\).
Set \(X := (0,1) \times \overline \Omega\),
\[
  \mathcal{M} := \mathcal{M}(X) \times \mathcal{M}( X ; \mathbb{R}^d)\, ,
\]
and
\[
  \mathcal{D} := \{(\rho,m) \in \mathcal{M} \ : \ \partial_t \rho + \div m = 0 \ \ \text{in} \ \
  X \}\,,
\]
where the solutions of the continuity equation are intended in a distributional
sense, that is,
\begin{equation} \label{cont weak}
\int_{X}  \partial_t \varphi \, d\rho +
\int_{X}  \nabla \varphi \cdot  dm   = 0   \quad \text{for all} \quad \varphi \in C^{\infty}_c ( X ) \,.
\end{equation}

We remark that the above weak formulation includes no-flux boundary conditions for the momentum $m$ on $\partial \Omega$. Also, no initial and final data is prescribed in \eqref{cont weak}. Moreover, by standard approximation arguments, we can consider in \eqref{cont weak} test functions in $C^1_c ( X )$ (see \cite[Remark 8.1.1]{ags}).

We now introduce the Benamou--Brenier energy.
For this purpose, define the convex, lower semicontinuous and one-homogeneous map $\Psi \colon \mathbb{R} \times \mathbb{R}^d \to [0,\infty]$ by setting
\
\Psi(t,x):=
\begin{cases}
\frac{|x|^2}{2t} & \text{ if } t>0\,, \\
0    & \text{ if } t=|x|=0 \,, \\
+ \infty& \text{ otherwise}\,.
\end{cases}
\]
The Benamou--Brenier energy $B: \mathcal{M} \to [0,\infty]$ is defined
for every pair $(\rho,m) \in \mathcal{M}$ as
\begin{equation} \label{intro convex}
B (\rho,m):= \int_X \Psi \left( \frac{d\rho}{d\lambda}, \frac{dm}{d\lambda}\right) \, d \lambda \,,
\end{equation}
$\lambda \in \mathcal{M}^+(X)$ is such that $\rho,m \ll \lambda$.
Since $\Psi$ is one-homogeneous, the above representation of $B$ does not depend on $\lambda$.
For some
fixed $\alpha \geq 0$, $\beta > 0$, we consider the following functional
\begin{equation} \label{prel:reg}
  J_{\alpha, \beta}(\rho,m) := \begin{cases}
      \beta B(\rho, m) +
      \alpha \norm{\rho}_{\mathcal{M}(X)} & \,\, \text{ if }(\rho,m) \in \mathcal{D}, \\
      + \infty\qquad  & \,\,   \text{ otherwise},
    \end{cases}
\end{equation}
where $\norm{\cdot}_{\mathcal{M}(X)}$ denotes the total variation norm in  \(\mathcal{M}(X)\).
\begin{rem} \label{rem:compact}
Note that in the definition of $J_{\alpha,\beta}$ we add the total variation of $\rho$ to the Benamou--Brenier energy. If $\alpha > 0$ this choice enforces the balls of the energy $J_{\alpha, \beta}$ to be compact in the weak* topology of $\mathcal{M}$ (see Lemma \ref{lem:prop J}). As a consequence, the functional $J_{\alpha,\beta}$  is a natural regularizer for dynamic inverse problems when the initial and final data are not prescribed \cite{bf}.
We remark that, although in the case $\alpha = 0$ the unit ball of the energy $J_{0, \beta}$ is not compact, we can still characterize its extremal points. However, in this case, due to the lack of coercivity, $J_{0,\beta}$ has limited use as a regularizer for dynamic inverse problems.
\end{rem}

For a measure $\rho \in \mathcal{M}(X)$, we say that $\rho$ disintegrates with respect to time if there exists a Borel family of measures $\{\rho_t\}_{t \in [0,1]}$ in %
$\mathcal{M}(\overline{\Om})$ such that
	\[
	\int_X \varphi(t,x) \, d\rho(t,x) = \int_0^1 \int_{\overline{\Om}} \varphi (t,x) \, d\rho_t(x) \, dt \quad \text{ for all } \quad  \varphi \in L^1_{\rho}(X) \,,
	\]
We denote such disintegration with the symbol $\rho =dt \otimes  \rho_t$. Further, we say that a curve of measures $t \in [0,1] \mapsto \rho_t \in \mathcal{M}(\overline{\Om})$ is narrowly continuous if the map
\[
t \mapsto \int_{\overline{\Om}} \varphi(x) \, d\rho_t(x)
\]
is continuous for each fixed $\varphi \in C(\overline{\Om})$. The family of narrowly continuous curves will be denoted by $C_{\rm w} ([0,1];\M(\overline{\Om}))$. We also introduce $C_{\rm w} ([0,1];\M^+(\overline{\Om}))$, as the family of narrowly continuous curves with values into the positive measures on $\overline{\Om}$.

We now recall several results about $B$, $J_{\alpha,\beta}$ and measure solutions of the continuity equation \eqref{cont weak}, which will be useful in the following analysis. For proofs of such results, we refer the interested reader to Propositions 2.2, 2.4, 2.6 and Lemmas 4.5, 4.6 in \cite{bf}.

\begin{lem}[Properties of $B$] \label{lem:prop B}
	The functional $B$ defined in \eqref{intro convex} is convex, positively one-homogen\-eous and sequentially lower semicontinuous with respect to the weak* topology on $\mathcal{M}$. Moreover it satisfies the following properties:
	\begin{enumerate}
	\item[i)] $B(\rho,m) \geq 0$ for all $(\rho,m)\in \mathcal{M}$,
	\item [ii)] if $B(\rho,m)<+\infty$, then $\rho \geq 0$ and $m \ll \rho$, that is, there exists a measurable map $v \colon X \to \mathbb{R}^d$ such that $m=v\rho$,
	\item [iii)] if $\rho \geq 0$ and $m = v \rho$ for some $v \colon X \to \mathbb{R}^d$ measurable, then
\begin{equation} \label{formula B}
B(\rho,m) =\int_X \Psi (1,v) \, d\rho =  \frac12 \int_X |v|^2 \, d\rho \,.
\end{equation}
	\end{enumerate}
\end{lem}

\begin{lem}[Properties of the continuity equation] \label{lem:prop cont}
	Assume that $(\rho,m) \in \mathcal{M}$ satisfies \eqref{cont weak} and that $\rho \in \mathcal{M}^+(X)$. Then $\rho$ disintegrates with respect to time into $\rho =dt \otimes  \rho_t$,
	%
	%
	%
	%
	where $\rho_t \in \mathcal{M}^+(\overline{\Omega})$ for a.e.~$t$. Moreover $t \mapsto \rho_t(\overline{\Omega})$ is constant, with $\rho_t(\overline{\Omega}) = \rho(X)$ for a.e. $t \in (0,1)$. If in addition $B(\rho,m)< +\infty$, that is,
	\[
	\int_0^1\int_{\overline{\Om}} |v|^2 \, d\rho_t(x) \, dt < +\infty\,,
	\]
where $m=v\rho$ for some $v \colon X \to \mathbb{R}^d$ measurable,
then $t \mapsto \rho_t$ belongs to $C_{\rm w} ([0,1];\M^+(\overline{\Om}))$.
\end{lem}

\begin{lem}[Properties of $J_{\alpha,\beta}$] \label{lem:prop J}
Let $\alpha \geq 0$, $\beta >0$. The functional $J_{\alpha,\beta}$ is non-negative, convex, positively one-homogeneous and sequentially lower semicontinuous with respect to weak* convergence on $\mathcal{M}$. Assume now $\alpha>0$.
		For $(\rho,m) \in \mathcal{M}$ such that $J_{\alpha,\beta}(\rho,m)< +\infty$ we have
	\begin{equation} \label{lem:prop J est}
          \alpha \norm{\rho}_{\mathcal{M}(X)} \leq J_{\alpha,\beta}(\rho,m)
          \,, \quad \min(2\alpha,\beta) \norm{m}_{\mathcal{M}(X;\mathbb{R}^d)} \leq
          J_{\alpha,\beta}(\rho,m) \,.
	\end{equation}
	Moreover, if $\{(\rho^n,m^n)\}_n$ is a sequence in $\mathcal{M}$ such that
\[
	\sup_n \ J_{\alpha,\beta} (\rho^n,m^n) < + \infty\,,
\]
then $\rho^n = dt \otimes \rho_t^n$ for some $(t \mapsto \rho_t^n) \in C_{\rm w} ([0,1];\M^+(\overline{\Om}))$ and there exists some
	$(\rho,m) \in \mathcal{D}$ with $\rho= dt \otimes \rho_t$, $\rho_t \in C_{\rm w} ([0,1];\M^+(\overline{\Om}))$ such that, up to subsequences,
	\begin{equation} \label{topology} \left\{
	\begin{gathered}
	(\rho^n,m^n) \weakstar (\rho,m) \,\, \text{ weakly* in } \,\, \mathcal{M} \,, \\
	\rho_t^n \weakstar \rho_t \,\, \text{ weakly* in }  \,\, \mathcal{M}(\overline{\Om})\,,\, \text{ for every } \,\, t \in [0,1]\,.
	\end{gathered} \right.
	\end{equation}
\end{lem}

\section{Characterization of extremal points}\label{sec:characterizationextremal}

The aim of this section is to characterize the extremal points of the unit ball of the functional $J_{\alpha,\beta}$ at \eqref{prel:reg} for all $\alpha \geq 0, \beta>0$, namely, of the convex set
\[
C_{\alpha,\beta} := \left\{(\rho,m) \in \mathcal{M} \, \colon \, J_{\alpha,\beta}(\rho,m)  \leq 1 \right\}\,.
\]
To this end, let us first introduce the following set.
\begin{definition}[Characteristics] \label{def:char}
For $\alpha\geq 0$, $\beta > 0$ define the set $\mathcal{C}_{\alpha,\beta}$ of all the pairs $(\rho,m) \in \mathcal{M}$ such that
\[
\rho=a_\gamma \, dt \otimes  \delta_{\gamma(t)}\,, \,\,\,\,\, m=\dot{\gamma}(t)\, a_\gamma \, dt \otimes  \delta_{\gamma(t)} \,, \,\,\,\,\,
a_{\gamma} := \left( \frac{\beta}{2} \int_0^1 |\dot\gamma(t)|^2 \, dt + \alpha \right)^{-1}\,,
\]
where $\gamma \in {\rm AC}^2([0,1];\mathbb{R}^d)$ satisfies $\gamma(t) \in \overline{\Om}$ for each $t \in [0,1]$ and $a_{\gamma}<+\infty$.
\end{definition}
We remind that ${\rm AC}^2([0,1];\mathbb{R}^d)$ denotes the space of absolutely continuous curves having a weak derivative in $L^2$.
We point out that by definition $a_{\gamma} >0$ for all choices of $\alpha\geq 0$, $\beta>0$. Moreover the condition $a_{\gamma}<+\infty$ is always satisfied if $\alpha>0$. When $\alpha=0$ we instead have $a_{\gamma}<+\infty$ if and only if $\int_0^1 |\dot{\gamma}(t)|^2\, dt > 0$, that is, the set $\mathcal{C}_{0,\beta}$ does not contain constant curves.

For the extremal points of $C_{\alpha,\beta}$ we have the following characterization.
\begin{thm} \label{thm:extremal}
Let $\alpha\geq 0$, $\beta>0$ be fixed. Then
	\[
	\ext(C_{\alpha,\beta})=\{(0,0) \} \cup \mathcal{C}_{\alpha,\beta} \,.
	\]
\end{thm}

The proof of Theorem \ref{thm:extremal} is postponed to Section \ref{subsection:proof}. In order to show the inclusion
$\ext(C_{\alpha,\beta}) \subset \{(0,0)\} \cup \mathcal{C}_{\alpha,\beta}$ we will make
use of a superposition principle for measure solutions of the continuity equation \eqref{cont weak}. This result is not new, and it is proved in \cite[Ch 8.2]{ags} for the case $\Omega = \mathbb{R}^d$. In Section \ref{subsec:superposition}  we show that it also holds for bounded closed domains.

\subsection{The superposition principle} \label{subsec:superposition}

 Before stating the superposition principle in $\overline{\Om}$, we introduce the following notation.
Let
\[
\Gamma := \left\{ \gamma \colon [0,1] \to \mathbb{R}^d  \, \colon \, \gamma \text{ continuous} %
\right\}
\]
be equipped with the supremum norm, i.e., $\norm{\gamma}_{\infty}:=\max_{t \in [0,1]} |\gamma(t)|$. For every fixed $t \in [0,1]$ let $e_t \colon \Gamma \to \mathbb{R}^d$ be the evaluation at $t$, that is,  $e_t(\gamma):=\gamma(t)$. Notice that $e_t$ is continuous.  For a measurable vector field $v \colon (0,1) \times \mathbb{R}^d \to \mathbb{R}^d$, we define the following subset of $\Gamma$ consisting of ${\rm AC}^2$ curves solving the ODE \eqref{ode} in the sense of Carath\'eodory:
\[
\Gamma_{v}(\mathbb{R}^d) := \left\{ \gamma \in \Gamma  \, \colon \, \gamma \in {\rm AC}^2([0,1];\mathbb{R}^d),\,\,  \dot{\gamma}(t)=v(t,\gamma(t)) \text{ for a.e. } t \in (0,1) \right\}\,.
\]
Moreover define the set of solutions to the ODE which live inside $\overline{\Om}$ for all times:
\[
\Gamma_{v}(\overline{\Om}) := \left\{ \gamma \in \Gamma_v (\mathbb{R}^d)\, \colon \, \gamma(t) \in \overline{\Om} \,\, \text{ for all } \,\, t \in [0,1]
\right\}  \,.
\]
The superposition principle for probability solutions to \eqref{cont weak} states as follows.

\begin{thm} \label{thm:lifting}
Let $t \in [0,1] \mapsto \rho_t \in \prob (\overline{\Om})$ be a narrowly continuous solution of the continuity equation in the sense of \eqref{cont weak}, for some measurable $v \colon (0,1) \times \overline{\Om} \to \mathbb{R}^d$ such that
\begin{equation} \label{thm:lifting:1}
\int_0^1 \int_{\overline{\Om}} |v(t,x)|^2 \, d\rho_t (x) \, dt < + \infty\,.
\end{equation}
Then there exists a probability measure $\sigma \in \prob(\Gamma)$ concentrated on $\Gamma_{v}( \overline{\Om})$ and such that $\rho_t = (e_t)_\# \sigma$ for every $t \in [0,1]$, that is,
\begin{equation} \label{representation:1}
\int_{\overline{\Om}} \varphi (x) \, d\rho_t (x) = \int_\Gamma \varphi (\gamma(t)) \, d\sigma(\gamma) \quad \text{ for every } \,\, \varphi \in C(\overline{\Om}),\,\, t \in [0,1]\,.
\end{equation}
\end{thm}

\begin{proof}
Let $\bar v :(0,1) \times \mathbb{R}^d \rightarrow \mathbb{R}^d$ be the extension to zero of $v$ to the whole $\mathbb{R}^d$.
Similarly, for each $t \in [0,1]$, let $\bar{\rho}_t \in \prob(\mathbb{R}^d)$ be the extension to zero of $\rho_t$ in $\mathbb{R}^d$. Note that the pair $(\bar \rho, \bar v \, \bar \rho)$ is a solution of the continuity equation in $(0,1) \times \mathbb{R}^d$ in the sense of \eqref{cont weak}. Moreover $\bar{\rho}$ and $\bar{v}$ satisfy \eqref{thm:lifting:1} in $(0,1) \times \mathbb{R}^d$. Therefore we can apply Theorem 8.2.1 in \cite{ags} and obtain a probability measure $\sigma \in \prob (\Gamma)$ concentrated on $\Gamma_{\bar{v}} (\mathbb{R}^d)$ and such that $\bar{\rho}_t = (e_t)_\# \sigma$ for all $t \in [0,1]$, that is,
\begin{equation} \label{representation:2}
\int_{\mathbb{R}^d} \varphi (x) \, d\bar{\rho}_t (x) = \int_\Gamma \varphi (\gamma(t)) \, d\sigma(\gamma) \quad \text{ for every } \,\, \varphi \in C_b(\mathbb{R}^d),\,\, t \in [0,1]\,.
\end{equation}
We claim that $\sigma$ is concentrated on $\Gamma_{v}(\overline{\Om})$. In order to show that, partition $\Gamma_{\bar{v}} (\mathbb{R}^d)$ into
\[
\Gamma_{\bar{v}}(\mathbb{R}^d) = \Gamma_{\bar{v}}(\overline{\Om}) \cup A \,,
\]
where
\[
A:= \left\{ \gamma \in \Gamma_{\bar{v}}(\mathbb{R}^d) \, \colon \, \text{there exists } \,\, \hat{t} \in [0,1] \,\, \text{such that} \,\, \gamma(\hat{t}) \in \overline{\Om}^c \right\} \,.
\]
Notice that, since $\overline{\Om}^c$ is open and $v \equiv  0$ in $\overline{\Om}^c$, the curves in $A$ are constant, so that we can write
\[
A = \left\{ \gamma \in \Gamma_{\bar{v}}(\mathbb{R}^d) \, \colon \,  \gamma(0) \in \overline{\Om}^c \right\} \,.
\]
From this, it follows that $A \subset e_0^{-1}(\overline{\Om}^c)$. Moreover, \eqref{representation:2} implies $\bar{\rho}_0 (\overline{\Om}^c) = \sigma(e_0^{-1}(\overline{\Om}^c))$. Therefore, using that $\bar{\rho}_t$ is concentrated on $\overline{\Om}$, we conclude that $\sigma(A)=0$, showing that $\sigma$ is concentrated on $\Gamma_{\bar{v}}(\overline{\Om})$.
Finally, \eqref{representation:2} implies \eqref{representation:1} since $\bar{\rho}_t$ is supported in $\overline{\Om}$ and it coincides with $\rho_t$ in $\overline{\Om}$. Also $\Gamma_{\bar{v}}(\overline{\Om})=\Gamma_{v}(\overline{\Om})$ by definition of $\bar{v}$, thus concluding the proof.
\end{proof}

\subsection{Proof of Theorem \ref{thm:extremal}} \label{subsection:proof}

Let $\alpha \geq 0$, $\beta >0$. We divide the proof into two parts.

\smallskip

\noindent\textit{Part 1:} $\{ (0,0)\} \cup \mathcal{C}_{\alpha,\beta} \subset \ext(C_{\alpha,\beta})$.
\smallskip

\noindent We start by showing that $\{(0,0)\}\cup \mathcal{C}_{\alpha,\beta} \subset C_{\alpha,\beta}$. The fact that $(0,0) \in C_{\alpha,\beta}$ follows immediately, since $(0,0)$ solves the continuity equation and $J_{\alpha,\beta}(0,0)=0$ (by Lemma \ref{lem:prop B}). Consider now $(\rho,m) \in \mathcal{C}_{\alpha,\beta}$.
Notice that $(\rho,m) \in \mathcal{C}_{\alpha,\beta}$ satisfies the continuity equation in the sense of \eqref{cont weak}: indeed for every $\varphi \in C^1_c ((0,1) \times \overline{\Om})$ we have
\begin{equation} \label{lem:inclusion:13}
\begin{aligned}
\int_{(0,1) \times \overline{\Om}}  \partial_t \varphi \, d \rho +  \nabla \varphi \cdot dm & = a_{\gamma} \int_0^1  \partial_t \varphi (t,\gamma(t)) + \nabla \varphi (t,\gamma(t)) \cdot \dot{\gamma}(t) \, dt
  \\
& = a_{\gamma} \int_0^1 \frac{d}{dt} \varphi (t,\gamma(t)) \, dt =a_{\gamma} ( \varphi (1,\gamma(1)) - \varphi(0,\gamma(0)) ) = 0
\end{aligned}
\end{equation}
since $\varphi$ is compactly supported in $(0,1) \times \overline{\Om}$.
Moreover, thanks to the fact that $\rho \geq 0$ and $m=\dot{\gamma} \rho$, we can invoke \eqref{formula B} to obtain
\begin{equation} \label{inclusion:one}
J_{\alpha,\beta}(\rho,m)
= a_{\gamma} \left( \frac{\beta}{2} \int_0^1  |\dot{\gamma}(t)|^2  \, dt + \alpha\right) = 1\,,
\end{equation}
proving that $(\rho,m) \in C_{\alpha,\beta}$.

We now want to show that any $(\rho,m) \in \{(0,0)\} \cup \mathcal{C}_{\alpha,\beta}$ is an extremal point for $C_{\alpha,\beta}$. Hence assume that $(\rho^1,m^1), (\rho^2,m^2) \in C_{\alpha,\beta}$ are such that
\begin{equation}\label{lem:inclusion:1}
	(\rho,m)=\lambda (\rho^1,m^1) + (1-\lambda) (\rho^2,m^2)
\end{equation}
for some $\lambda \in (0,1)$. We need to show that $(\rho,m)=(\rho^1,m^1) = (\rho^2,m^2)$.
Set $j \in \{1,2\}$.
Since $(\rho^j,m^j)$ is such that $J_{\alpha,\beta}(\rho^j,m^j) \leq 1$, from $ii)$ in Lemma \ref{lem:prop B} we have that $\rho^j \geq 0$ and $m^j=v^j\, \rho^j $ for some Borel field $v^j \colon X \to \mathbb{R}^d$. In particular, if $(\rho,m)=(0,0)$, \eqref{lem:inclusion:1} forces $(\rho^j,m^j)=0$, hence showing that $(0,0)$ is an extremal point of $C_{\alpha,\beta}$.

Let us now consider the case $(\rho,m) \in \mathcal{C}_{\alpha,\beta}$.
By \eqref{inclusion:one} we have $J_{\alpha,\beta}(\rho,m)=1$.
From \eqref{lem:inclusion:1}, convexity of $J_{\alpha,\beta}$, and the fact that $J_{\alpha,\beta}(\rho^j,m^j) \leq 1$, $\lambda \in (0,1)$, we conclude
\begin{equation} \label{lem:inclusion:11}
J_{\alpha,\beta}(\rho^j,m^j)=1\,.
\end{equation}
Since $(\rho^j,m^j)$ solves the continuity equation, $\rho^j \geq 0$ and $J_{\alpha,\beta}(\rho^j,m^j) = 1$, from Lemma \ref{lem:prop cont} we deduce that $\rho^j=dt \otimes \rho_t^j$ for some narrowly continuous curve $t \mapsto \rho^j_t \in \mathcal{M}^+(\overline{\Om})$, with $\rho^j_t(\overline{\Om})$ constant in time. We define $a_j:=\rho^j_0(\overline{\Om})$ and notice that $a_j >0$: Indeed, $a_j=0$ would imply $\rho^j=0$, yielding $J_{\alpha,\beta}(\rho^j,m^j)=J_{\alpha,\beta}(0,0)=0$. This would contradict \eqref{lem:inclusion:11}.
Now, from condition \eqref{lem:inclusion:1}
and uniqueness of the disintegration we deduce
\begin{equation} \label{lem:inclusion:6}
a_{\gamma}\,\delta_{\gamma(t)} = \lambda  \rho_t^1 + (1-\lambda) \rho_t^2  \quad \text{ for every } t  \in [0,1] \,.
\end{equation}
Since $a_j >0$ (and hence $\rho^j_t \neq 0$), the above equality implies that $\supp \rho_t^j = \{\gamma(t)\}$, i.e.,
\begin{equation} \label{lem:inclusion:12}
\rho_t^j =a_j \,\delta_{\gamma(t)} \quad \text{ for every } t  \in [0,1] \,.
\end{equation}
We now show that $v^j=\dot{\gamma}$ on $\supp \rho = \gr(\gamma):=\{ (t,\gamma(t)) \colon t \in (0,1)\}$, that is
\begin{equation} \label{lem:inclusion:field}
	v^j(t,\gamma(t))= \dot\gamma(t) \quad \text{ for a.e. } \,\, t \in (0,1)\,.
\end{equation}
By assumption, $\partial_t \rho^j + \div m^j=0$ in the sense of \eqref{cont weak}. %
Therefore, recalling \eqref{lem:inclusion:12} and the fact that $a_j >0$, we get that for each $\varphi \in C^1_c((0,1) \times \overline{\Om})$,
\begin{equation} \label{lem:inclusion:3}
\begin{aligned}
0 & = \int_0^1 \partial_t \varphi (t, \gamma (t)) + \nabla \varphi (t,\gamma(t)) \cdot v^j(t, \gamma(t)) \, dt  \\
& =  \int_0^1 \partial_t \varphi (t, \gamma (t)) + \nabla \varphi (t,\gamma(t)) \cdot \dot{\gamma}(t) \, dt \, +
 \int_0^1 \nabla \varphi (t,\gamma(t)) \cdot ( v^j(t, \gamma(t)) - \dot{\gamma}(t) ) \, dt \\
& = \int_0^1 \nabla \varphi (t,\gamma(t)) \cdot ( v^j(t, \gamma(t)) - \dot{\gamma}(t) ) \, dt\,,
\end{aligned}
\end{equation}
where the last equality follows from \eqref{lem:inclusion:13}, since $a_\gamma >0$. Let $\psi \in C^1_c ((0,1))$ and define $\varphi(t,x):=x_i \psi (t) $, where $x=(x_1,\dots,x_d)$, so that $\varphi$ is a test function for \eqref{lem:inclusion:3}. By plugging $\varphi$ into \eqref{lem:inclusion:3} we obtain
\[
\int_0^1  \psi (t) ( v^j_i(t, \gamma(t)) - \dot{\gamma}_i(t) ) \, dt = 0 \qquad \text{ for all} \,\,\, i \in \{1,\ldots,d\}, \,\,j\in \{1,2\}
\]
where $v_i^j$ and $\dot{\gamma}_i$ are the $i$-th component of $v^j$ and $\dot{\gamma}$, respectively.
This implies that $v^j(t,\gamma(t))=\dot{\gamma}(t)$ for a.e.~$t \in (0,1)$, that is, $v^j=\dot{\gamma}$ a.e.~on $\gr(\gamma)$. With this at hand, by means of \eqref{formula B} we can see that
$J_{\alpha,\beta}(\rho^j,m^j)=a_j /a_\gamma$. Since \eqref{lem:inclusion:11} holds, we obtain $a_j= a_\gamma$, thus proving $(\rho,m)=(\rho^j,m^j)$ and hence extremality for $(\rho,m)$ in $C_{\alpha,\beta}$.

\medskip
\noindent\textit{Part 2:} $\ext(C_{\alpha,\beta}) \subset \{(0,0)\} \cup \mathcal{C}_{\alpha,\beta}$.
\smallskip

\noindent Let $(\rho,m) \in C_{\alpha,\beta}$ be an extremal point. In particular, $J_{\alpha,\beta}(\rho,m) \leq 1$ so that by Lemma \ref{lem:prop B} $ii)$, we obtain $\rho \geq 0$ and $m=v \rho $ for some Borel field $v \colon X \to \mathbb{R}^d$. Notice that by extremality of $(\rho,m)$ and one-homogeneity of $J_{\alpha,\beta}$ we immediately infer
that either
$J_{\alpha,\beta}(\rho,m)=0$ or $J_{\alpha,\beta}(\rho,m)=1$.
If $J_{\alpha,\beta}(\rho,m)=0$, by decomposing $(\rho,m)$ as
\[
(\rho,m) = \frac{1}{2}(2\rho,2m) + \frac{1}{2}(0,0)
\]
and using the extremality of $(\rho,m)$ together with the one-homogeneity of $J_{\alpha,\beta}$ we deduce that $(\rho,m) = (0,0)$.
Thus, we consider the case
\begin{equation}\label{thm:extremal:3}
J_{\alpha,\beta}(\rho,m)=1\,.
\end{equation}
Since by definition, $(\rho,m)$ solves the continuity equation in the sense of \eqref{cont weak} and $J_{\alpha,\beta}(\rho,m) = 1$, we can apply Lemma \ref{lem:prop cont} to obtain that $\rho=a \, dt \otimes \rho_t$ for some narrowly continuous curve $t \mapsto \rho_t \in \prob (\overline{\Om})$, where $a:=\rho(X)>0$.

\medskip
\noindent\textit{Claim:} $\supp \rho_t$ is a singleton for each $t \in [0,1]$.
\medskip

\noindent\textit{Proof of Claim:} The hypotheses of Theorem \ref{thm:lifting} are satisfied, therefore there exists a measure $\sigma \in \prob(\Gamma)$ concentrated on $\Gamma_v (\overline{\Om})$ and such that $\rho_t = (e_t)_\# \sigma$ for every $t \in [0,1]$.
Assume by contradiction that there exists a time $\hat{t} \in [0,1]$ such that $\supp \rho_{\hat{t}}$ is not a singleton.
 Therefore, we can find a Borel set $E \subset \overline{\Om}$ such that
\begin{equation} \label{thm:extremal:2bis}
0<\rho_{\hat{t}} (E)<1, \,\,\,\,\,\, 0< \rho_{\hat{t}} (\overline{\Om} \smallsetminus E )<1 \,.
\end{equation}
Define the Borel set
\[
A:=\{ \gamma \in \Gamma \, \colon \, \gamma(\hat{t}) \in E \} = e^{-1}_{\hat{t}} (E) \,.
\]
By the relation $\rho_t = (e_t)_\# \sigma$ and definition of $A$ we obtain
$\rho_{\hat t}(E) =\sigma(A)$.
Therefore, from \eqref{thm:extremal:2bis}
\begin{equation}\label{thm:extremal:2}
 0<{\sigma} (A)<1,\,\, \, \,\,\, 0 < \sigma (\Gamma \smallsetminus A )<1 \,.
\end{equation}
Define
\[
\begin{aligned}
\lambda_1 & :=a \left( \frac{\beta}{2} \int_0^1 \int_{A} |\dot\gamma(t)|^2 \, d\sigma(\gamma) \, dt + \alpha \, \sigma(A) \right)\, ,\\
\lambda_2 & :=a \left( \frac{\beta}{2} \int_0^1 \int_{A^c} |\dot\gamma(t)|^2 \, d\sigma(\gamma)  \, dt + \alpha \, \sigma(A^c)\right) \,,
\end{aligned}
\]
where $A^c:=\Gamma \smallsetminus A$.
Note that $\lambda_1,\lambda_2$ are well defined (possibly being equal to $+\infty$) as the map
\begin{equation} \label{eq:mapL}
L(\gamma):=\int_0^1 |\dot{\gamma}(t)|^2\, dt  \,\,\,\, \text{ if } \,\, \gamma \in {\rm AC}^2([0,1];\mathbb{R}^d) \,, \,\,\,\, L(\gamma):=+\infty \,\,\,\, \text{ otherwise},
\end{equation}
is lower semicontinuous on $\Gamma$, and hence measurable.
Notice that
\begin{equation} \label{thm:extremal:4}
\lambda_1 + \lambda_2 = a \left(\frac{\beta}{2} \int_0^1 \int_{\Gamma} |\dot\gamma(t)|^2 \, d\sigma(\gamma) \, dt + \alpha \right) =
a \left( \frac{\beta}{2} \int_0^1 \int_{\Gamma} |v(t,\gamma(t))|^2 \, d\sigma(\gamma) \, dt + \alpha \right)\,,
\end{equation}
because $\sigma$ is concentrated on $\Gamma_v (\overline{\Om})$.
Since $v(t,\cdot)$ belongs to $L^2_{\rho_t}(\overline{\Om};\mathbb{R}^d)$ for a.e.~$t \in (0,1)$, by \cite[Theorem 3.6.1]{bogachev} we obtain that the representation formula \eqref{representation:1} holds for $\varphi(x):=v(t,x)$ and a.e.~$t \in (0,1)$, that is,
\begin{equation} \label{thm:extremal:5}
 \int_{\Gamma} |v(t,\gamma(t))|^2 \, d\sigma(\gamma) =
 \int_{\overline{\Om}} |v(t,x)|^2 \, d\rho_t(x) \,\,\, \text{ for a.e.} \,\,
 t \in (0,1)\,.
\end{equation}
Therefore, from \eqref{formula B}, \eqref{thm:extremal:3}, \eqref{thm:extremal:4} and \eqref{thm:extremal:5} we deduce $\lambda_1+\lambda_2=J_{\alpha,\beta}(\rho,m)=1$.

We now proceed with the proof of the claim separately for the cases $\alpha > 0$ and $\alpha = 0$.
Suppose first $\alpha> 0$. Notice that $\lambda_1,\lambda_2 > 0$ thanks to \eqref{thm:extremal:2} and the fact that $a>0$.
Decompose $(\rho,m)$ as
\begin{equation} \label{thm:extremal:6}
(\rho,m)= \lambda_1 (\rho^1, m^1) + \lambda_2 (\rho^2, m^2 ) \,,
\end{equation}
where we defined
\begin{equation} \label{thm:extremal:7}
\rho^j:= \frac{a}{\lambda_j} \, dt \otimes (e_t)_\# \sigma_j \,, \quad m^j:=\rho^j v \,,
\end{equation}
for $j=1,2$, with $\sigma_1 := \sigma \zak A$ and $\sigma_2 := \sigma \zak A^c$. Notice that $\rho^j \in \mathcal{M}^+(X)$, since $\sigma$ is a positive measure concentrated on $\Gamma_v (\overline{\Om})$, and $a,\lambda_j>0$. We now claim that $(\rho^j,m^j) \in C_{\alpha,\beta}$. First, we prove that $\partial_t \rho^j + \div m^j=0$ in the sense of \eqref{cont weak}. Let $j=1$ and fix $\varphi \in C^1_c((0,1)\times \overline{\Om})$. Since $v(t,\cdot)$ belongs to $L^2_{\rho_t}(\overline{\Om};\mathbb{R}^d)$ for a.e. $t \in (0,1)$, by \cite[Theorem 3.6.1]{bogachev},  \eqref{representation:1} and the definition of $\sigma_1$, we get
\begin{equation}\label{eq:continuityeqextremality1}
\begin{aligned}
\int_X \partial_t \varphi \, d\rho^1 +  \nabla \varphi \cdot d m^1 =
& \frac{a}{\lambda_1} \int_0^1 \int_{\overline{\Om}}  \partial_t \varphi (t,x) + \nabla \varphi (t,x) \cdot v (t,x) \, d \left( (e_t)_\# \sigma_1 \right)(x) \, dt \\
 = & \frac{a}{\lambda_1}\int_0^1 \int_{A}  \partial_t \varphi (t,\gamma(t)) + \nabla \varphi (t,\gamma(t)) \cdot v(t,\gamma(t)) \, d \sigma (\gamma) \, dt\,.
\end{aligned}
\end{equation}
Now recall that $\sigma$ is concentrated on $\Gamma_v (\overline{\Om})$ and that $\varphi$ is compactly supported in time, so that
\begin{equation}\label{eq:continuityeqextremality2}
\begin{aligned}
\int_X \partial_t \varphi \, d\rho^1 +  \nabla \varphi \cdot d m^1
 = & \frac{a}{\lambda_1}\int_0^1 \int_{A}  \partial_t \varphi (t,\gamma(t)) + \nabla \varphi (t,\gamma(t)) \cdot \dot{\gamma}(t) \, d \sigma (\gamma) \, dt\\
 = &\frac{a}{\lambda_1}  \int_A \left( \int_0^1 \frac{d}{dt} \varphi (t,\gamma(t)) \, dt \right) \, d\sigma(\gamma) = 0 \,.
\end{aligned}
\end{equation}
The calculation for $j=2$ is similar. Also,
by definition of $(\rho^j,m^j)$ and of $\lambda_j$, one can perform similar calculations to the ones in \eqref{thm:extremal:4}, \eqref{thm:extremal:5}, and prove that $J_{\alpha,\beta}(\rho^j,m^j)=1$. Hence $(\rho^j,m^j) \in C_{\alpha,\beta}$. We now claim that $(\rho^1,m^1) \neq (\rho^2,m^2)$. Suppose by contradiction that $(\rho^1,m^1) = (\rho^2,m^2)$. Then in particular $\rho^1 = \rho^2$, so that by \eqref{thm:extremal:7} we get
\begin{equation}\label{diff}
\frac{(e_t)_\# \sigma_1}{\lambda_1} =\frac{ (e_t)_\# \sigma_2 }{\lambda_2}\quad \text{for a.e. } t\in (0,1)\,.
\end{equation}
As $(\rho^j,m^j)$ are solutions of the continuity equation and $J_{\alpha,\beta}(\rho^j,m^j)=1$, from Lemma \ref{lem:prop cont}  it follows that the maps $ t \mapsto (e_t)_\# \sigma_j$ are narrowly continuous. In particular, \eqref{diff} holds for each $t \in [0,1]$. However, by \eqref{thm:extremal:2} and by definition of $A$, $\sigma_1$, $\sigma_2$,  we have
\[
 [(e_{\hat t})_\# \sigma_1](E) =   \sigma(A)>0 \,, \quad  \quad  [(e_{\hat t})_\# \sigma_2](E) = \sigma(\emptyset) = 0 \,,
\]
which contradicts \eqref{diff}. Therefore $(\rho^1,m^1) \neq (\rho^2,m^2)$, which shows that the decomposition \eqref{thm:extremal:6} is non-trivial. This is a contradiction, since we are assuming that $(\rho,m)$ is an extremal point for $C_{\alpha,\beta}$. Thus the claim follows.

Suppose now that $\alpha = 0$ and define the set
\begin{equation*}
Z := \left\{\gamma \in \Gamma : \int_0^1 |\dot \gamma(t)|^2 \, dt  = 0\right\}\,.
\end{equation*}
Notice that $Z$ is measurable, due to the measurability of the map $L$ at \eqref{eq:mapL}. We claim that $\sigma(Z)=0$. In order to prove that, let $Z^c:= \Gamma \smallsetminus Z$ and define the measures $\sigma_Z :=\sigma \zak Z$, $\sigma_{Z^c} :=\sigma \zak Z^c$, so that $\sigma=\sigma_Z + \sigma_{Z^c}$. Recalling that $\rho_t=(e_t)_{\#}\sigma$ for all $t \in [0,1]$, we can decompose
\begin{equation}\label{eq:convexcomb}
(\rho,m) = \frac{1}{2}(\rho^1,m^1) + \frac{1}{2}(\rho^2,m^2) \,,
\end{equation}
where
\begin{equation*}
\begin{aligned}
\rho^1 & := a \, dt \otimes (e_t)_\# \sigma_{Z^c} \,, \,\,\,\, & m^1:= v \rho^1\, ,\\
\rho^2 & :=  a\, dt  \otimes  (e_t)_\# \sigma_{Z^c}  + 2a \, dt \otimes (e_t)_\# \sigma_{Z}\,, \,\,\,\, & m^2:= v \rho^2\,.
\end{aligned}
\end{equation*}
Let $j=1,2$. Notice that $\rho^j \in \mathcal{M}^+(X)$ since $\sigma$ is a positive measure concentrated on $\Gamma_v (\overline{\Om})$ and $a > 0$. Following similar computation as \eqref{eq:continuityeqextremality1}-\eqref{eq:continuityeqextremality2} we infer that $(\rho^j,m^j)$ solves the continuity equation in the sense of \eqref{cont weak}. Moreover, by definition of $Z$ and the fact that $\sigma$ is concentrated on $\Gamma_v (\overline{\Om})$ we obtain
\begin{equation} \label{eq:fubini}
\begin{aligned}
	\int_0^1 \int_{Z^c} |v(t,\gamma(t))|^2 \, d\sigma(\gamma)\,dt & =
	\int_0^1 \int_{Z^c} |\dot\gamma(t)|^2 \, d\sigma(\gamma)\,dt \\
	& = \int_{Z^c} \int_0^1 |\dot\gamma(t)|^2 \, dt \, d\sigma(\gamma)\\
	 & =  \int_{\Gamma} \int_0^1 |\dot\gamma(t)|^2 \, dt \, d\sigma(\gamma) =  \int_0^1\int_{\Gamma} |v(t,\gamma(t))|^2 \, d\sigma(\gamma) \, dt \,,
\end{aligned}
\end{equation}
where we employed Fubini's Theorem, which holds thanks to the measurability of the map $L$ and the identity \eqref{thm:extremal:5}, the latter implying boundedness of the last term in \eqref{eq:fubini}. By \eqref{eq:fubini} and arguing as in \eqref{thm:extremal:4}-\eqref{thm:extremal:5}, it is immediate to check that $J_{0,\beta}(\rho^j,m^j)=J_{0,\beta}(\rho,m)$. Recalling \eqref{thm:extremal:3} we then obtain $(\rho^j,m^j) \in C_{0,\beta}$.
As $(\rho,m)$ is an extremal point of $C_{0,\beta}$, from \eqref{eq:convexcomb} we deduce that $(\rho^1,m^1) = (\rho^2,m^2)$ and thus $ dt \otimes (e_t)_\# \sigma_{Z} = 0$. In particular, there exists $\hat t \in [0,1]$ such that $(e_{\hat t})_\# \sigma_{Z} = 0$. Hence for every $E \subset \Gamma$ measurable, by the positivity of $\sigma$, we have $\sigma_{Z}(E) \leq  (e_{\hat t})_\# \sigma_{Z}(e_{\hat t}(E))=0$, implying that $\sigma_{Z} = 0$. By \eqref{thm:extremal:2} and the definition of $\lambda_1$, $\lambda_2$, we conclude that $\lambda_1,\lambda_2 > 0$. With this property established, the claim that $\supp \rho_t$ is a singleton for each $t \in [0,1]$ follows by repeating the same arguments of the case $\alpha>0$, employing the decomposition of $(\rho,m)$ as in \eqref{thm:extremal:6}.

\medskip

We have shown that for each $t \in [0,1]$,  $\supp \rho_t$ is a singleton. We now conclude the proof of Theorem \ref{thm:extremal}. Since $\rho_t \in \prob(\overline{\Om})$, the latter implies the existence of a curve $\gamma \colon [0,1] \to \overline{\Om}$ such that $\rho_t=\delta_{\gamma(t)}$ for each $t \in [0,1]$. We will now prove that $\gamma \in {\rm AC}^2([0,1];\mathbb{R}^d)$. By narrow continuity of $t \mapsto \rho_t$, we have that the map $t \mapsto  \, \varphi(\gamma(t))$
is continuous for all $\varphi \in C(\overline{\Om})$. By testing against the coordinate functions $\varphi(x):=x_i$, we obtain continuity for $\gamma$. Consider now $\varphi(t,x):=\xi(t)\eta(x)$ with $\xi \in C^{\infty}_c((0,1))$, $\eta \in C^1(\overline{\Om})$. Notice that the scalar map $t \mapsto \eta(\gamma(t))$ is continuous. Moreover, by testing the continuity equation $\partial_t \rho + \div ( v \rho)=0$ against $\varphi$ we get
\[
\int_0^1 \xi'(t) \, \eta(\gamma(t)) \, dt = - \int_0^1 \xi(t) \, \nabla \eta (\gamma(t)) \cdot v(t,\gamma(t)) \, dt  \,,
\]
which implies that the distributional derivative of the map $t \mapsto \eta(\gamma(t))$ is given by
\[
t \mapsto \nabla \eta (\gamma(t)) \cdot v(t,\gamma(t)) \,.
\]
We now remark that the above map belongs to $L^2((0,1))$, since
\[
\int_0^1 |\nabla \eta (\gamma(t)) \cdot v(t,\gamma(t))|^2 \, dt \leq \norm{\nabla \eta}_{\infty} \int_0^1 |v(t,\gamma(t))|^2 \, dt \leq C\, \norm{\nabla \eta}_{\infty} J_{\alpha,\beta}(\rho,m) < + \infty\,.
\]
Therefore, $t \mapsto \eta(\gamma(t))$ belongs to ${\rm AC}^2([0,1])$ for every fixed $\eta \in C^1(\overline{\Om})$. By choosing $\eta(x):=x_i$, $i = 1 , \dots, d$, we conclude that $\gamma \in {\rm AC}^2([0,1];\mathbb{R}^d)$. Since $\partial_t \rho+ \div (v\rho )=0$, we can repeat the
same argument employed to prove \eqref{lem:inclusion:field}, and infer
\begin{equation}  \label{lem:inclusion:30}
v(t,\gamma(t)) = \dot{\gamma}(t) \,\, \text{ a.e. in } \,\, (0,1) \,.
\end{equation}
From \eqref{lem:inclusion:30} and the fact that $\rho=a \, dt \otimes \delta_{\gamma(t)}$, $m=v\rho$, we then conclude $m=\dot\gamma \rho$. As $a>0$, we can apply \eqref{formula B} to compute
\begin{equation}\label{eq:lastJ}
J_{\alpha,\beta}(\rho,m)=a \left( \frac{\beta}{2} \int_0^1 |\dot{\gamma}(t)|^2 \, dt + \alpha \right) \,.
\end{equation}
Recalling that $J_{\alpha,\beta}(\rho,m)=1$ (see \eqref{thm:extremal:3}), from \eqref{eq:lastJ} we conclude that $a=a_{\gamma}$ with $a_{\gamma}<+\infty$. Therefore $(\rho,m)$ belongs to $\mathcal{C}_{\alpha,\beta}$ according to Definition \ref{def:char}, and the proof of Theorem \ref{thm:extremal} is concluded.

\section{Application to sparse representation for inverse problems with optimal transport regularization} \label{sec:application}

In this section we deal with the problem of reconstructing a family of time-dependent Radon measures given a finite number of observations. To be more specific, let $H$ be a finite dimensional Hilbert space and $A: C_{\rm w} ([0,1];\M(\overline{\Om})) \rightarrow H$ be a linear continuous operator, where continuity is understood in the following sense:
given a sequence $(t \mapsto \rho^n_t)$ in $C_{\rm w} ([0,1];\M(\overline{\Om}))$, we require that
\begin{equation}\label{eq:topologycweak}
\rho_t^n  \weakstar \rho_t \quad \text{ weakly* in } \,\, \mathcal{M}(\overline{\Om})  \,\,\text{ for all }\,\, t\in [0,1] \,\, \text{ implies } \,\, A\rho^n \rightarrow A\rho \ \text{ in } H\,,
\end{equation}
where, with a little abuse of notation, we will denote by $\rho^n$ both the curve $t \mapsto \rho_t^n$, as well as the measure $\rho^n:=dt \otimes \rho^n_t$.

For some given data $y \in H$, we aim to reconstruct a solution $\rho \in C_{\rm w} ([0,1];\M(\overline{\Om}))$ to the dynamic inverse problem
\begin{equation} \label{appl:inverse}
A \rho = y \,.
\end{equation}
For $\alpha >0$ and $\beta >0$ we regularize the above inverse problem by means of the energy $J_{\alpha,\beta}$ defined in \eqref{prel:reg}, following the approach in \cite{bf}. In practice, upon introducing the space
\[
\widetilde{\mathcal{M}}:= C_{\rm w} ([0,1];\M(\overline{\Om})) \times \mathcal{M}(X;\mathbb{R}^d)\,,
\]
we consider the Tikhonov functional $G:\widetilde{\mathcal{M}}\rightarrow \mathbb{R} \cup \{+\infty\}$ defined as
\begin{equation}\label{eq:inversefunctional}
G(\rho,m) =J_{\alpha,\beta}(\rho,m) + F(A\rho)\,,
\end{equation}
where $F: H \rightarrow \mathbb{R} \cup\{+\infty\}$ is a given fidelity functional for the data $y$, which is assumed to be convex, lower semicontinuous and bounded from below. Additionally, we assume that $G$ is proper.
We then replace \eqref{appl:inverse} by
\begin{equation}\label{eq:inverseproblem}
\min_{(\rho,m) \in \widetilde{\mathcal{M}} } \ G(\rho,m) \, .
\end{equation}

\begin{rem}
Two common choices for the fidelity term $F$ in the case $H=\mathbb{R}^k$ are, for example,
\begin{itemize}
\item [i)] $F(x) = I_{\{y\}}(x)$ for a given $y \in \mathbb{R}^k$ that forces the constraint $A\rho = y$,
\item[ii)] $F(x) = \frac{1}{2}\|x-y\|_2^2$ that recovers a classical $l^2$ penalization.
\end{itemize}
\end{rem}

\begin{rem}\label{rem:existence}
Under the above assumptions on $A$ and $F$, problem \eqref{eq:inverseproblem} admits a solution.
Indeed, since $G$ is proper, any minimizing sequence $\{(\rho^n,m^n)\}_n$ is such that $\{G(\rho^n,m^n)\}_n$ is bounded. As $F$ is bounded from below and $J_{\alpha,\beta} \geq 0$, we deduce that $\{J_{\alpha,\beta}(\rho^n,m^n)\}_n$ is bounded. Therefore, Lemma \ref{lem:prop J} implies that $(\rho^n,m^n)$ converges (up to subsequences) to some $(\rho,m) \in \widetilde{\mathcal{M}}$, in the sense of \eqref{topology}. By weak* lower semicontinuity of $J_{\alpha,\beta}$ in $\mathcal{M}$ (see Lemma \ref{lem:prop J}) and by \eqref{eq:topologycweak} together with the lower-semicontinuity of $F$, we infer that $(\rho,m)$ solves \eqref{eq:inverseproblem}.
 %
\end{rem}

It is well-known that the presence of a finite-dimensional constraint in an inverse problem, such as \eqref{appl:inverse}, promotes sparsity in the reconstruction. This observation has been recently made rigorous in \cite{bc} and \cite{chambolle}, where it has been shown that the atoms of a sparse minimizer are the extremal points of the ball of the regularizers. In Theorem \ref{thm:extremal}, we provided a characterization for the extremal points of the ball of $J_{\alpha,\beta}$. Therefore, specializing the above-mentioned results to our setting yields the following characterization theorem for sparse minimizers to \eqref{eq:inverseproblem}.

\begin{thm}\label{thm:sparserepresentation}
Let $\alpha,\beta>0$. There exists a minimizer $(\hat \rho, \hat m) \in \widetilde{\mathcal{M}}$ of \eqref{eq:inverseproblem} that can be represented as
\begin{equation}
(\hat \rho, \hat m) = \sum_{i=1}^p c_i \, (\rho^i,m^i)\,,
\end{equation}
where $p\leq \dim(H)$, $c_i >0$, $\sum_{i=1}^p c_i = J_{\alpha,\beta}(\hat \rho,\hat m)$, and
\[
\rho^i= a_{\gamma_i} \, dt \otimes  \delta_{\gamma_i(t)} \,, \,\,\, m^i = \dot{\gamma}_i \,  \rho^i\,,
\]
where $\gamma_i \in {\rm AC}^2([0,1];\mathbb{R}^d)$ with $\gamma(t) \in \overline{\Om}$ for each $t \in [0,1]$, and $a_{\gamma_i}^{-1}:=\frac{\beta}{2} \int_0^1 |\dot\gamma_i|^2 \, dt + \alpha$.
\end{thm}
In other words, the above theorem ensures the existence of a minimizer of \eqref{eq:inverseproblem} which is a finite linear combination of measures concentrated on the graphs of ${\rm AC}^2$-trajectories contained in $\overline{\Om}$. In Section \ref{subsec:proof2} we give a proof of Theorem \ref{thm:sparserepresentation}, and we conclude the paper with Section \ref{subsec:dynamic}, where we apply the sparsity result of Theorem \ref{thm:sparserepresentation} to dynamic inverse problems with optimal transport regularization, following the approach of \cite{bf}.

\subsection{Proof of Theorem \ref{thm:sparserepresentation}} \label{subsec:proof2}
As already mentioned, the proof is an immediate consequence of Theorem \ref{thm:extremal} and a particular case of \cite[Corollary $2$]{chambolle} (see also \cite[Theorem $1$]{chambolle}). Before proceeding with the proof, for the reader's convenience, we recall Corollary $2$ from \cite{chambolle}. The definitions appearing in the statement below will be briefly recalled in the proof of Theorem \ref{thm:sparserepresentation}.

\begin{thm}[\cite{chambolle}]\label{thm:chambolle}
Let $\mathcal{U}$ be a locally convex space, $H$ be a finite-dimensional Hilbert space, $R:\mathcal{U} \rightarrow [-\infty,+\infty]$, $F: H\rightarrow  [-\infty,+\infty]$ be convex, and $A:\mathcal{U} \rightarrow H$ be linear.
Consider the variational problem
\begin{equation}\label{eq:thmchambolle}
\inf_{u\in \mathcal{U}} \, R(u) + F(Au)\, .
\end{equation}
Suppose that the set of minimizers of \eqref{eq:thmchambolle}, denoted by $S$, is non-empty. Additionally, assume that there exists $\hat u \in \ext(S)$ such that the set
\begin{equation}
\overline{C} = \{u\in \mathcal{U}: R(u) \leq R(\hat u)\}
\end{equation}
is \emph{linearly closed}, the \emph{lineality space} of $\overline{C}$ is $\{(0,0)\}$
and $\inf_{u\in \mathcal{U}} R(u) < R(\hat u)$. Then, exactly one of the following conditions holds:
\begin{itemize}
\item[i)] $\hat{u}$ is a convex combination of at most $\dim (H)$ extremal points of $\overline{C}$,
\item[ii)] $\hat{u}$ is as a convex combination of at most  $\dim (H) -1$ points, which are either extremal points of $\overline{C}$, or belong to an extreme ray of $\overline{C}$.
\end{itemize}
\end{thm}

\begin{proof}[Proof of Theorem \ref{thm:sparserepresentation}]
We just need to verify that we can apply Theorem \ref{thm:chambolle} to the variational problem \eqref{eq:inverseproblem}. So, we choose $\mathcal{U}= \widetilde M$, $R=J_{\alpha,\beta}$ and $F$ and $A$ satisfying the assumptions stated above. Let $S$ be the set of solutions to \eqref{eq:inverseproblem}.

First, notice that in Remark \ref{rem:existence} we have already shown that  the set of minimizers for \eqref{eq:inverseproblem} is non-empty, so that $S \neq \emptyset$. Moreover $S$ is compact with respect to the weak* topology. Indeed, given a sequence $(\rho^n,m^n)$ in $S$ we can use Lemma \ref{lem:prop J} to extract a subsequence (not relabelled) such that $(\rho^n,m^n) \weakstar (\rho,m)$ in $\mathcal{M}$ and $\rho_t^n \weakstar \rho_t$ in $\mathcal{M}(\overline{\Om})$ for every $t\in [0,1]$. Using the sequential lower semicontinuity of $J_{\alpha,\beta}$ with respect to weak* convergence combined with the continuity of $A$ (according to \eqref{eq:topologycweak}) and the lower semicontinuity of $F$, we obtain $(\rho,m) \in S$. We conclude that $S$ is sequentially weakly* compact and hence weakly* compact, thanks to the metrizability of the weak* convergence on bounded sets. Finally note that $S$ is convex thanks to the convexity of $F$ and $J_{\alpha,\beta}$ (Lemma \ref{lem:prop J}). By Krein--Milman's Theorem, we then infer the existence of a $(\hat{\rho},\hat m) \in \ext(S)$.

The lineality space of $\overline{C}$ is defined as $\text{lin}(\overline{C}) = \text{rec}(\overline{C}) \cap (-\text{rec}(\overline{C}))$, where $\text{rec}(\overline{C})$ is the recession cone of $\overline{C}$ defined as the set of all $(\rho,m)\in \mathcal{U}$ such that $\overline C +\mathbb{R}_+(\rho,m)  \subset \overline C$. Hence, from the coercivity of $J_{\alpha,\beta}$ in Lemma \ref{lem:prop J} it is immediate to conclude that $\text{lin}(\overline{C}) = \{(0,0)\}$.
Moreover, $\overline{C}$ is linearly closed if the intersection of $\overline{C}$ with every line is closed. It is easy to verify that, as $\overline{C}$ is weakly* closed (Remark \ref{rem:existence}), it is also linearly closed. Finally, the assumption  $\inf_{(\rho,m)\in \widetilde M} J_{\alpha,\beta}(\rho,m) < J_{\alpha,\beta}(\hat \rho, \hat m)$ is satisfied whenever $(\hat \rho, \hat m) \neq 0$, as in this case $J_{\alpha,\beta}( \hat \rho,  \hat m) >0$, while $\inf_{(\rho,m)\in \widetilde M} J_{\alpha,\beta}(\rho,m)=0$.
Hence, the hypotheses of Theorem \ref{thm:chambolle} for the functional \eqref{eq:inversefunctional} are verified. Notice also that $\overline C$ does not contain extreme rays. In order to prove that, we first recall that a ray of $\overline C$ is any set of the form $r_{p,v} = \{p+ tv : t >0\}$ for $p, v\in \overline C$, $v\neq 0$. An extreme ray of $\overline C$ is a ray $r_{p,v}$ such that for every segment intersecting $r_{p,v}$, the whole segment is contained in $r_{p,v}$.
Thanks to the coercivity of $J_{\alpha,\beta}$ in Lemma \ref{lem:prop J}, it is immediate to see that $\overline C$ contains no rays and thus no extreme rays.
Hence, from either of the conclusions $i),ii)$ in Theorem \ref{thm:chambolle}, we deduce that there exists a minimizer $(\hat \rho, \hat m) \in \mathcal{M}$ of \eqref{eq:inverseproblem} that can be represented as
\begin{equation}\label{eq:1representation}
(\hat \rho, \hat m) = \sum_{i=1}^p c_i (\rho_i,m_i)\, ,
\end{equation}
where $(\rho_i,m_i)\in \ext(C_{\alpha,\beta})$, $p\leq \dim(H)$, $c_i>0$ and $\sum_{i=1}^p c_i = J_{\alpha,\beta}(\hat \rho,\hat m)$.
We remark that if $(\hat \rho, \hat m) = 0$, the assumption $\inf_{(\rho,m)\in \widetilde M} J_{\alpha,\beta}(\rho,m) < J_{\alpha,\beta}(\hat \rho, \hat m)$ in Theorem \ref{thm:chambolle} is not satisfied, but the representation \eqref{eq:1representation} holds trivially.
Using the characterization of extremal points in Theorem \ref{thm:extremal} and \eqref{eq:1representation}, we obtain an explicit sparse representation for solutions of \eqref{eq:inverseproblem} and the proof is achieved.
\end{proof}

\subsection{Dynamic inverse problems} \label{subsec:dynamic}

Theorem \ref{thm:sparserepresentation} provides a representation formula for sparse solutions of \eqref{eq:inverseproblem} that holds for every $A$ and $F$ satisfying the above-stated hypotheses.
A relevant choice for $A$ and $F$ is proposed in \cite{bf} as a model for dynamic inverse problems: in particular, the authors apply their framework to variational reconstruction in undersampled dynamic MRI. In what follows we make an explicit choice of $F$ and $A$ in order to apply Theorem \ref{thm:sparserepresentation} to a special case of the framework in \cite{bf}, namely the case of discrete time sampling, and finite-dimensionality of the data for each sampled time. %

To be more specific, consider a discretization of the interval $[0,1]$ in $N$ points $t_1<t_2<\ldots<t_N$ and assume that we want to reconstruct an element of $C_{\rm w} ([0,1];\M(\overline{\Om}))$, by only making observations at the time instants $t_1,\ldots,t_N$. To this aim, let $H_{t_i}$ be a family of finite-dimensional Hilbert spaces and introduce the product space $\mathcal{H}:=\bigtimes_{i=1}^N H_{t_i}$, normed by $\|y\|_{\mathcal{H}}^2 := \sum_{i=1}^N \|y_i\|^2_{H_{t_i}}$. Let
$A_{t_i} : \mathcal{M}(\overline{\Om}) \to H_{t_i}$ be linear operators, which are assumed to be weak*  continuous for each $i=1,\ldots,N$.
For a given observation $(y_{t_1},\ldots,y_{t_N}) \in  \mathcal{H}$, consider the problem of finding $\rho \in C_{\rm w} ([0,1];\M(\overline{\Om}))$ such that
\[
A_{t_i} \rho_{t_i} = y_{t_i} \,\, \text{ for each } \,\,
i = 1 ,\ldots,N\,.
\]
Following \cite{bf}, we regularize the above problem by
\begin{equation}\label{eq:mri}
\min_{(\rho,m)\in \widetilde{\mathcal{M}}} \, J_{\alpha,\beta}(\rho,m) + \frac{1}{2}\sum_{i=1}^N \|A_{t_i}\rho_{t_i} - y_{t_i}\|^2_{H_{t_i}} \,.
\end{equation}
In order to recast the above problem into the form \eqref{eq:inverseproblem}, let $A : C_{\rm w} ([0,1];\M(\overline{\Om})) \rightarrow \mathcal{H}$ be the linear operator defined by
\[
A\rho := (A_{t_1}\rho_{t_1}, \ldots,A_{t_N}\rho_{t_N})\,.
\]
Notice that $A$ is continuous in the sense of \eqref{eq:topologycweak}, thanks to the assumptions on $A_{t_i}$. We can then equivalently rewrite \eqref{eq:mri} as
\begin{equation}\label{eq:mri2}
\min_{(\rho,m)\in \widetilde{\mathcal{M}}} \, J_{\alpha,\beta}(\rho,m) + \frac{1}{2}\|A\rho - y\|^2_{\mathcal{H}} \,.
\end{equation}
 In this way, we recover a problem of the type of \eqref{eq:inverseproblem}, where $F(x) := \frac{1}{2}\|x - y\|^2_{\mathcal{H}}$. Notice that $F$ is convex, lower semicontinuous and bounded from below. Moreover, the functional in  \eqref{eq:mri2} is proper, since $J_{\alpha,\beta}(0,0)=0$. Hence, we can apply Theorem \ref{thm:sparserepresentation} to conclude the following result.

 \begin{cor} \label{corollary}
Let $\alpha,\beta>0$.
 There exists a minimizer $(\hat \rho, \hat m) \in \widetilde{\mathcal{M}}$ of \eqref{eq:mri} that can be represented as
\begin{equation}\label{eq:final}
(\hat \rho, \hat m) =  \sum_{i=1}^p c_i \, (\rho^i,m^i)\,,
\end{equation}
where $p\leq \dim(\mathcal{H}) = \sum_{i=1}^N \dim(H_i)$, $c_i >0$, $\sum_{i=1}^p c_i = J_{\alpha,\beta}(\hat \rho,\hat m)$, and
\[
\rho^i= a_{\gamma_i} \, dt \otimes  \delta_{\gamma_i(t)} \,, \,\,\, m^i = \dot{\gamma}_i \,  \rho^i\,,
\]
where $\gamma_i \in {\rm AC}^2([0,1];\mathbb{R}^d)$ with $\gamma(t) \in \overline{\Om}$ for each $t \in [0,1]$, and $a_{\gamma_i}^{-1}:=\frac{\beta}{2} \int_0^1 |\dot\gamma_i|^2 \, dt + \alpha$.
 \end{cor}

\begin{rem}
The upper bound $p \leq \sum_{i=1}^N \dim (H_i)$ in the representation formula \eqref{eq:final} might not be optimal. However, a careful analysis of the faces of the ball of the Benamou--Brenier energy, possibly under additional assumptions on the operator $A$ and fidelity term $F$, could be needed to substantiate such conjecture. We leave this question open for future research.
\end{rem}

\section*{Acknowledgements}
We thank the referee for the useful suggestions provided, particularly for encouraging us to include the case $\alpha=0$ in the characterization of Theorem \ref{thm:extremal}, which was previously missing.
KB and SF gratefully acknowledge support by the Christian Doppler Research Association (CDG) and Austrian Science Fund (FWF) through the Partnership in Research
project PIR-27 ``Mathematical methods for motion-aware medical imaging'' and project P 29192 ``Regularization graphs for variational imaging''.
MC is supported by the Royal Society (Newton International Fellowship NIF\textbackslash R1\textbackslash 192048 Minimal partitions as a robustness boost for neural network classifiers). The Institute of Mathematics and Scientific Computing, to which KB, SF and FR are affiliated, is a member of NAWI Graz (\texttt{http://www.nawigraz.at/en/}). The authors  KB, SF and FR are further members of/associated with BioTechMed Graz (\texttt{https://biotechmedgraz.at/en/}).

	%