text:
\section{Introduction}
Optimal transport (OT) \citep[for a detailed mathematical discussion see e.g.][]{villani2008optimal,santambrogio2015optimal} has been a focus of attention in various research fields for a long time. More recently, its powerful geometric features promoted by improved computational tools \citep[see e.g.][]{chizat2018scaling,peyre2019computational,guo2020fast,lin2020fixed} have turned OT into a promising new tool for modern data analysis with applications in machine learning \citep{frogner2015learning,arjovsky2017wasserstein,schmitz2018wasserstein,yang2018low,vacher2021dimension}, computer vision \citep{baumgartner2018visual,kolkin2019style}, genetics \citep{evans2012phylogenetic,schiebinger2019optimal}, cell biology \citep{gellert2019substrate,klatt2020empirical,tameling2021colocalization}, image processing \citep{pitie2007automated,rabin2015convex,bonneel2016wasserstein,tartavel2016wasserstein} and statistical inference \citep{sommerfeld2018inference,lee2018minimax,mena2019statistical,panaretos2019statistical}, among others.\\
However, the wide range of such applications also surfaced some limitations of classical OT. In particular, the assumption of equal total mass intensity of the measures is often inappropriate. A straightforward strategy to overcome this issue is to normalise the measures' total intensities. However, this preprocessing step has an immediate impact on data analysis, in particular on the underlying transport plan. For example, when matching point clouds of different sizes the resulting plan distributes mass among several points, whereas often it is desired to match points one-to-one which is favourable in many applications (see \Cref{fig:uotvsot}). Attempts to circumvent this issue have led to a range of \emph{unbalanced optimal transport} (UOT) proposals \citep[see e.g.][]{figalli2010optimal,liero2018optimal,chizat2018unbalanced,balaji2020robust,mukherjee2021outlier,heinemann2022kantorovich}. These formulations extend optimal transport concepts to general positive measures by either fixing the total amount of mass to be transported in advance or by penalising the hard marginal constraints inherent in OT. All these approaches also give rise to associated barycenters, generalising the popular notion of \emph{OT barycenters} \citep{agueh2011barycenters} to measures of unequal total intensity \citep{chizat2018scaling,friesecke2021barycenters,heinemann2022kantorovich}. \\
The UOT formulation considered here is the $(p,C)$-Kantorovich-Rubinstein distance (KRD) (see \Cref{sec:krd}) whose structural properties and related $(p,C)$-barycenter have recently been studied in detail \citep{heinemann2022kantorovich}. The $(1,1)$-KRD essentially corresponds to the notion of extended Kantorovich norms \citep{kantorovich1958space,hanin1992kantorovich} considered in the context of Lipschitz spaces and signed measures. In this context, \cite{guittet2002extended} first introduced a discrete formulation of the problem, where he established a Linear Program (LP) formulation which carries over to the general $(p,C)$-KRD.
For illustration, a comparison between the $(p,C)$-barycenter and the $p$-Wasserstein barycenter in a simple example is displayed in \Cref{fig:krvswsbary}. From a data analysis point of view we find it particularly appealing that for the $(p,C)$-KRD there is a clear geometrical connection between its penalty $C$ and the structural properties of the corresponding UOT plans and $(p,C)$-barycenter. More precisely, it is shown in \cite{heinemann2022kantorovich} that $C$ controls the largest scale at which mass transport is possible in an optimal plan. This interpretation of the $(p,C)$-KRD allows to easily design it to respect different structural properties of the data and thus makes it a prime candidate for statistical tasks in data analysis. Furthermore, each support point of any $(p,C)$-barycenter is contained in a finite set characterised by the value of $C$. This allows to adapt OT solvers for the unbalanced problems.\\
\begin{figure}
  \centering
  \subfloat[][]{\includegraphics[width=0.4\linewidth]{graphics/otvsuot1v2.png}}%
  \qquad
  \subfloat[][]{\includegraphics[width=0.4\linewidth]{graphics/otvsuot2.png}}%
  \caption{Transport between two measures (blue and brown) with their support points located in $[0,1]^2.$ The respective transport plans between them are displayed by red lines where the thickness of a line is proportional to the transported mass. \textbf{(a)} The measures have been normalised to probability measures (the blue points have mass $1/6$ and the brown points with mass $1/13$). \textbf{(b)} The UOT plan for the $(2,2)$-KRD between the two unnormalised measures (all points have mass $1$).}
  \label{fig:uotvsot}
\end{figure}
Though, due to the unbalanced nature of the problem, the task of sampling from the underlying measures requires alternative sampling schemes and different statistical modelling. While for OT between probability measures there is a canonical sampling model by i.i.d. replications from the measures, this fails for UOT, since the considered measures are not necessarily probability measures. In this work, we address this issue and suggest a framework underpinning UOT based statistical data analysis. To this end, we analyse the $(p,C)$-KRD and its barycenter in three specific statistical models motivated by applications in randomised algorithms and microscopy tasks. Notably, these models also provide a framework which potentially allows to treat the alternative UOT models mentioned above.
\subsection{Kantorovich-Rubinstein Distance}\label{sec:krd}
Let $(\mathcal{X},d)$ be a finite metric space with finite cardinality $M$ and denote by
\begin{align*}
\mathcal{M}_+(\mathcal{X}) \coloneqq \left\lbrace \mu \in \mathbb{R}^{\vert\mathcal{X}\vert} \,\mid\, \mu(x)\geq 0 \ \forall x \in \mathcal{X} \right\rbrace
\end{align*}
the set of non-negative measures\footnote{A non-negative measure on a finite space $\mathcal{X}$ is uniquely characterized by the values it assigns to each singleton $\{ x\}$. To ease notation we write $\mu(x)$ instead of $\mu(\{x\})$. The corresponding $\sigma$-field is always to be understood as the powerset of $\mathcal{X}$.} on $\mathcal{X}$. For a measure $\mu\in\mathcal{M}_+(\mathcal{X})$ its total mass is defined as $\mathbb{M}(\mu)\coloneqq \sum_{x\in\mathcal{X}}\mu(x)$ and the subset $\mathcal{P}(\mathcal{X})\subset \mathcal{M}_+(\mathcal{X})$ of measures with total mass one is the set of probability measures. If $\pi\in \mathcal{M}_+(\mathcal{X}\times\mathcal{X})$ is a measure on the product space $\mathcal{X} \times \mathcal{X}$ its marginals are defined as $\pi(x,\mathcal{X})\coloneqq \sum_{x^\prime} \pi(x,x^\prime)$ and $\pi(\mathcal{X},x^\prime)\coloneqq \sum_{x^\prime \in \mathcal{X}}\pi(x,x^\prime)$, respectively. For two measures $\mu,\nu\in\mathcal{M}_+(\mathcal{X})$ define the set of \emph{non-negative sub-couplings} as
\begin{align}\label{eq:subcouplings}
\begin{split}
    \Pi_{\leq}(\mu,\nu)\coloneqq \lbrace \pi\in \mathcal{M}_+(\mathcal{X}\times\mathcal{X}) \mid \,&\pi(x,\mathcal{X}) \leq  \mu(x),\, \pi(\mathcal{X},x^\prime)  \leq  \nu(x^\prime)\, \forall \, x,x^\prime\in\mathcal{X}  \rbrace.
\end{split}
\end{align}
Following \cite{heinemann2022kantorovich}, for $p\geq 1$ and a parameter $C>0$, the $(p,C)-$\emph{Kanto\-rovich-Rubinstein distance} (KRD) between two measures $\mu,\nu\in\mathcal{M}_+(\mathcal{X})$ is defined as
\begin{align}\label{eq:krdistance}
\begin{split}
\text{KR}_{p,C}(\mu,\nu)\coloneqq \Bigg( \min_{\pi\in \Pi_{\leq}(\mu,\nu)} &\sum\limits_{x,x^\prime \in \mathcal{X}}d^p(x,x^\prime)\pi(x,x^\prime)\\&+C^p\left(\frac{\mathbb{M}(\mu)+\mathbb{M}(\nu)}{2}-\mathbb{M}(\pi)\right)\Bigg)^{\frac{1}{p}}.
\end{split}
\end{align}
For any $p\geq 1$, it defines a distance on the space of non-negative measures $\mathcal{M}_+(\mathcal{X})$ and it naturally extends the well-known \emph{$p$-th order OT distance} defined for measures of equal total mass \citep{heinemann2022kantorovich}.\\
\subsection{Kantorovich-Rubinstein Barycenters}\label{sec:pcbary}
\begin{figure}
  \centering
  \subfloat[][]{\includegraphics[width=0.32\linewidth]{graphics/wsvskrabrydatafull.png}}%
  \
  \subfloat[][]{\includegraphics[width=0.32\linewidth]{graphics/wsvskrbary_ws.png}}%
    \
  \subfloat[][]{\includegraphics[width=0.32\linewidth]{graphics/wsvskrbary_kr.png}}%
  \caption{\textbf{(a)} $J=10$ measures with mass $1$ at each support point and different total mass intensities (each colour corresponds to a different measure) superimposed on top of each other. \textbf{(b)} The OT barycenter (for squared Euclidean cost) of the normalised measures. \textbf{(c)} The $(2,0.3)$-barycenter of the unnormalised measures (see \eqref{eq:KRbarycenter} for a rigorous definition).}%
  \label{fig:krvswsbary}
\end{figure}
The $(p,C)$-KRD also allows to define a notion of a barycenter for a collection of measures with (potentially) different total masses. Assume $(\mathcal{X},d)$ to be embedded in some connected ambient space\footnote{We assume the metric on $\mathcal{X}\subset \mathcal{Y}$ to be the metric of $\mathcal{Y}$ restricted to $\mathcal{X}$.} $(\mathcal{Y},d)$, e.g., a Euclidean space. Define the (unbalanced) \emph{$(p,C)$-Fr\'echet functional}
\begin{equation}\label{eq:KRbarycenter}
        F_{p,C}(\mu)=\frac{1}{J}\easysum{i=1}{J}KR_{p,C}^p(\mu^i,\mu).
\end{equation}
Any minimiser of this functional in $\mathcal{M}_+(\mathcal{Y})$ is said to be a \emph{$(p,C)$-Kantorovich-Rubinstein barycenter} of $\mu_1,\dots,\mu_J$ or \emph{$(p,C)$-barycenter} for short\footnote{For the sake of readability, the weights in this definition are fixed to $1/J$. Adaptation of all results to arbitrary positive weights $\lambda_1,\dots ,\lambda_J$ summing to one is straightforward.}. The objective functional $F_{p,C}$ is referred as (unbalanced) \emph{$(p,C)$-Fr\'echet functional} and the so called \emph{Borel barycenter application} is defined as
\[
T^{L,p}(x_1,\dots ,x_L) \in  \argmin_{y \in \mathcal{Y}} \sum_{i=1}^L d^p(x_i,y).
\]
Define the \emph{full centroid set}\footnote{There are scenarios where multiple sets fulfil the definition of the centroid set, since there might be multiple points that minimise the barycentric application. In this case, a fixed representative is chosen and there still exists a choice of centroid set which contains the support of the $(p,C)$-barycenter.} of the measures
\begin{align}\label{eq:KRfullcentroid}
\begin{split}
\mathcal{C}_{KR}(J,p)=\Big\{ y \in \mathcal{Y} \ \vert \ \exists L\geq \lceil J/2 \rceil, \ \exists (i_1,\dots,i_L)\subset \{1,\dots,J\},\\ x_1,\dots,x_L: \ x_l\in \text{supp}(\mu_{i_l}) \\ \forall l=1,\dots,L: \
y=T^{L,p}(x_1,\dots,x_L)\Big\},
 \end{split}
\end{align}
and based on it the \emph{restricted centroid set}
\begin{align}
\begin{split}
\label{eq:KRCentroid}
\mathcal{C}_{KR}(J,p,C)=\Big\{ &y=T^{L,p}(x_1,\dots,x_L)\in\mathcal{C}_{KR}(J,p)\, \mid \forall 1\leq l \leq L: \\ &d^p(x_l,y)\leq C^p; \ \es{i=1}{L}d^p(x_l,y) \leq \frac{C^p (2L-J)}{2} \Big\}.
\end{split}
\end{align}
According to \cite{heinemann2022kantorovich}, any $(p,C)$-barycenter is finitely supported and its support is included in the restricted centroid set $\mathcal{C}_{KR}(J,p,C)$. This is critical, as it allows to restrict theoretical analysis as well as computational methods to the scenario all measures involved have finite support.
\subsection{Statistical Models and Contributions}\label{sec:sampling}
In practice, one often does not have access to the population measures $\mu,\nu,\mu^1,\dots ,\mu^J$, respectively. Instead, these measures have to be estimated from data. For probability measures the most common statistical model assumes access to (i.i.d.) data $X_1,\dots,X_N\sim \mu$ (and similar for $\nu,\mu^1,\dots,\mu^J$). A commonly used estimator is the empirical measure $\hat{\mu}_N=(1/N)\sum_{k=1}^{N}\delta_{X_k}$, where $\delta_X$ denotes the (random) point measure at location $X$. This can be used to estimate OT and its plan. In general, this yields a statistically efficient estimator \citep{dudley1969speed,ajtai1984optimal,sommerfeld2019optimal,weed2019sharp,hutter2021minimax,manole2021plugin,hundrieser2022empirical} see, however, \cite{niles2022minimax} for certain improvements. For OT barycenters based on empirical measures significantly less is known, though recently some progress has been made in the context of finitely supported measures \citep{heinemann2022randomized}. \\
Extending such results to general measures is not obvious and requires the need for alternative statistical modelling which is the topic of this paper. Exemplary, we propose three approaches motivated by different applications.
\subsubsection*{Multinomial Model}
For the \emph{multinomial model} the measures are normalised to define probability measures such that sampling as described above is possible. The resulting empirical estimators are rescaled to the original total intensities. More precisely, consider i.i.d. random variables $X_1,\ldots,X_N \stackrel{\text{i.i.d.}}{\sim} \frac{\mu}{\mathbb{M}(\mu)}$, where the total intensity $\mathbb{M}(\mu)$ is assumed to be known. The corresponding unbiased empirical estimator is then defined as
\begin{align}\label{eq:multinomialmeasure}
\hat{\mu}_N\coloneqq \frac{\mathbb{M}(\mu)}{N} \sum_{x \in \mathcal{X}} \lvert\{ k\in \{1,\dots,N\} \ \vert \ X_k=x \}\rvert \ \delta_{x}.
\end{align}
This approach can be understood as extension of the classical sampling approach for probability measures to measures of arbitrary total intensities. A key motivation for this model is resampling for randomised computations of UOT. In real world data analysis it is common to encounter data (e.g. high-resolution images) which are out of reach for current state of the art solvers for OT. One idea in this scenario is to replace each measure by its empirical version and then use these measures as surrogates. This idea has been considered in detail for the $p$-Wasserstein distance \citep{sommerfeld2019optimal} and the $p$-Wasserstein barycenter \citep{heinemann2022randomized}. Statistical deviation bounds allow to balance computational complexity and accuracy of approximation in terms of the sample size $N$. In the context of this model we extend these results to the $(p,C)$-KRD for measures of unequal total intensities.
\subsubsection*{Bernoulli Model}
For the \emph{Bernoulli model} we consider measures $\mu$ with $\mu(x)=1$ for all $x \in \text{supp}(\mu)$. Thus, the measure $\mu$ represents a point cloud in the ambient space $\mathcal{Y}$ with the total mass being the cardinality of the point cloud. We restrict the model to this setting, but we stress that generalisations to arbitrary masses at the individual locations are straightforward. We then assume to observe independent Bernoulli random variables $B_x{\sim} \text{Ber}(s_x)$ with a fixed \emph{success probability} $s_x\in [0,1]$ for each location $x\in\mathcal{X}$. We denote $s_\mathcal{X}:=(s_{x_1},\dots,s_{x_{M}})$ and refer to $s_\mathcal{X}$ as \emph{success vector}. A suitable unbiased estimator for $\mu$ is defined by
\begin{align}\label{eq:bermeasure}
\hat{\mu}_{s_\mathcal{X}}:= \sum_{x\in\mathcal{X}} \frac{B_x}{s_x}\delta_x.
\end{align}
The corresponding Bernoulli field $(B_x)_{x\in \mathcal{X}}$ is a prototypical model for incomplete data, where data is missing at random and the Bernoulli variables serve as labels for this. It further arises, e.g. in the context of generalised linear models where the regressor $X$ is linked to $B_x$ by a link function in a non-parametric fashion. The Bernoulli model also underlies the sampling scheme for the subsequent Poisson intensity model, which occurs e.g. in various imaging devices, such as fluorescence cell microscopy. There, fluorescent markers are, e.g. chemically attached to each protein within a complex protein ensemble and then are excited with a laser beam. The resulting, emitted photons indicate the spatial position of the objects of interest in the proper experimental setup \citep{kulaitis2021resolution}. However, the marker has a limited \emph{labelling efficiency} $s_x \in (0,1]$ at each location $x\in \mathcal{X}$ and we only observe a location which has been labelled by the marker and finally emits photons.
\subsubsection*{Poisson Intensity Model}
For the \emph{Poisson intensity model} we fix a parameter $t>0$ and a success probability $s\in [0,1]$. Consider a collection of $\lvert \mathcal{X} \rvert$ independent Poisson random variables $P_x{\sim} \text{Poi}(t\mu(x))$ with intensity $t\mu(x)$ at each location $x\in\mathcal{X}$ and independently from that Bernoulli random variables $B_x\sim \text{Ber}(s)$ for each $x \in X$. A suitable unbiased estimator for $\mu$ is defined by
\begin{align}\label{eq:poissonmeasure}
\hat{\mu}_{t,s}:= \frac{1}{st}\sum_{x\in\mathcal{X}} B_x P_x\delta_x.
\end{align}
In contrast to the Bernoulli model, in this model the success probability is assumed to be homogeneous (as opposed to the inhomogeneous probabilities in the Bernoulli model, though such generalisations are straightforward) and the values of the population measures at each support point are not necessarily equal to one. Hence, we have two independent layers of randomness in the construction of this empirical measure. First, we draw a location $x$ with a certain probability $s$, then we observe random photon counts driven by a Poisson distribution based on the mass of $\mu$ and the value of $t$.\\
This model is motivated by various tasks in photonic imaging, for example, fluorescence microscopy, X-ray imaging and positron emission tomography (PET), see \cite{munk2020statistical} for a survey. The finite space $\mathcal{X}$ represents the center of bins of a detection interface used to measure the emitted photons. The value $\mu(x)$ corresponds to the integrated underlying photon intensity over its respective bin. This intensity is proportional to an external source, such as a laser duration in fluorescence microscopy and modelled by the parameter $t>0$. The Bernoulli random variable $B_x$ models the possibility that in the bin of $x$ a photon can not be recorded. This might be due to various effects that cause thinning, such as limited labelling efficiency, dead time of cameras or a loss of photons due to sparse detector tubes. The value of $P_x$ corresponds to the number of photons which have been measured at the bin of $x$. Note, that besides $B_x$, there might be also additional effects present which do not disable the whole bin, but just prevent a single photon from being measured. All this causes a thinning of the process and is incorporated in the probability $s^\prime\in (0,1]$  that a single photon at any bin and any point in time can not be measured. In this case, the model can be shown to be equivalent to a Poisson model with parameter $ts^\prime >0$ instead of $t$ \citep{aspelmeier2015modern}. Hence, this kind of thinning corresponds to a reparametrisation of the original model and is thus a special case of this general Poisson intensity model.\\

\subsubsection{Sampling Bounds}
Let $\mu,\nu \in\mathcal{M}_+(\mathcal{X})$ and let $\hat{\mu},\hat{\nu}$ be empirical versions of these measures generated with one of the aforementioned sampling mechanisms (recall \eqref{eq:multinomialmeasure},\eqref{eq:bermeasure} and \eqref{eq:poissonmeasure}). We are interested in tight upper bounds for the quantities
\begin{align*}
    (i) \ \mathbb{E}\left[ KR_{p,C}(\hat{\mu},{\mu})\right],  \quad (ii) \ \mathbb{E}\left[ \lvert KR_{p,C}(\hat{\mu},\hat{\nu})-KR_{p,C}({\mu},{\nu})\rvert \right]. \\
\end{align*}
In particular, we show that there exist constants $\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Mult}}(C)$, $\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Pois}}(C)$, $\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Ber}}(C)$ such that for any $p\geq 1$ and for any measure $\mu$ and its estimator $\hat{\mu}$ in each of the three statistical models it holds
\begin{align*}
\mathbb{E}\left[ \text{KR}_{p,C}\left( \hat{\mu},\mu\right)\right]\leq
\begin{cases}
\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Mult}}(C)^{\frac{1}{p}} \, N^{-\frac{1}{2p}}, &\text{ if }\hat{\mu}=\hat{\mu}_N,\\[3ex]
\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Poi}}(C)^{\frac{1}{p}} \, \phi(t,s)^{\frac{1}{p}}, &\text{ if }\hat{\mu}=\hat{\mu}_{t,s},\\[3ex]
\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Ber}}(C)^{\frac{1}{p}} \, \psi(s_\mathcal{X})^{\frac{1}{p}}, &\text{ if }\hat{\mu}=\hat{\mu}_{s_\mathcal{X}},
\end{cases}
\end{align*}
where
\begin{align*}
    \phi(t,s)=\begin{cases}
   \left(2(1-s)\mathbb{M}(\mu)+\frac{s}{\sqrt{t}}\sum_{x\in \mathcal{X}}\sqrt{\mu(x)} \right), & C\leq \min_{x\neq x^\prime}d(x,x^\prime) \\
  \left( \frac{1}{st}\mathbb{M}(\mu)+\frac{1-s}{s}\sum_{x\in \mathcal{X}}\mu(x)^2 \right)^{\frac{1}{2}}, &\text{else}, \\
    \end{cases}
\end{align*}
and
\begin{align*}
   \psi(s_\mathcal{X})= \begin{cases}
   \left(2\sum_{x\in \mathcal{X}} (1-s_x) \right), & C\leq \min_{x\neq x^\prime}d(x,x^\prime) \\
  \left( \sum_{x\in \mathcal{X}}\frac{1-s_x}{s_x} \right)^{\frac{1}{2}}, &\text{else}. \\
    \end{cases}
\end{align*}
Notably, in the multinomial model for $N\rightarrow\infty$, in the Poisson model for $t\rightarrow\infty$, $s\rightarrow 1$ and in the Bernoulli model for $s_{\mathcal{X}}\rightarrow \mathbf{1}_M$, these upper bounds tend to zero. The proof of these bounds (see \Cref{chp:sampbounds}) is inspired by partitioning strategies used to control the empirical deviation of the OT distance between probability measures \citep{dereich2013constructive,fournier2015rate,weed2019sharp} and based on a tree approximation of the underlying space $\mathcal{X}$ \citep{sommerfeld2019optimal}. Finally, we make use of a recently established closed form solution of the $(p,C)$-KRD on ultra-metric trees by \cite{heinemann2022kantorovich}.\\
The constants in these bounds depend on the structural properties of the measures and space, such as covering numbers and total mass intensity. For a detailed discussion of the constants see \Cref{chp:sampbounds}. Exemplary, for measures supported in $([0,1]^D,d_2)$, $p=2$ and $d_{\min}^{\infty}:={\min_{x\neq x^\prime}d_\infty(x,x^\prime)}$, it holds for $D<4$ (see for details \Cref{sec:explicitB})
\begin{align*}
\mathcal{E}_{2,\mathcal{X},\mu}^{\text{Poi}}(C)\leq D
\begin{cases}
    \frac{C^2}{2\sqrt{D}}-\frac{2^{D/2+3}}{1-2^{D/2-2}},\hspace{0.5cm} &\text{if } C\geq 2\\[0.3cm]
  \frac{2^{5+l(D/2-2)}}{1-2^{D/2-2}},  \hspace{0.5cm}&\text{if }  2^{1-l}\leq C< 2^{2-l},\\[0.3cm]
\frac{C^2}{2}, \hspace{0.5cm}&\text{if } C\leq d_{\min}^{\infty} .
\end{cases}
\end{align*}
For $D=4$ it holds
\begin{align*}
\mathcal{E}_{2,\mathcal{X},\mu}^{\text{Poi}}(C)\leq 4
\begin{cases}
    \frac{C^2}{4}-2\left(2-\lvert \mathcal{X} \rvert^{-\frac{1}{4}} \right)^2+(2+8\log_2\left(\vert \mathcal{X}\vert\right)), \hspace{0.5cm}&\text{if }  C\geq 2,\\[0.3cm]
\left(2+8\log_2\left(\vert \mathcal{X}\vert\right)-l\right),  \hspace{0.5cm}&\text{if }   2^{1-l}\leq C< 2^{2-l},\\[0.3cm]
\frac{C^2}{2} \hspace{0.5cm}&\text{if }  C\leq d_{\min}^{\infty},
\end{cases}
\end{align*}
and for $D>4$, it holds
\begin{align*}
\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Poi}}(C)\leq D
\begin{cases}
     \frac{C^2}{2\sqrt{D}}-2\left(2-\lvert \mathcal{X} \rvert^{-\frac{1}{D}} \right)^2+2\vert \mathcal{X}\vert^{\frac{1}{2}-\frac{2}{D}} \left( 1+\frac{2^{2+D/2}}{2^{D/2-2}-1}\right) \hspace{0.5cm}&\text{if } C\geq 2,\\[0.3cm]
\bigg( \vert \mathcal{X}\vert^{\frac{1}{2}-\frac{2}{D}}+\frac{2^{2+D/2}}{2^{D/2-2}-1}\left( \lvert \mathcal{X}\rvert^{\frac{1}{2}-\frac{2}{D}}-2^{(D/2-2)(l-1)}
\right)\bigg), \hspace{0.5cm}&\text{if }   2^{1-l}\leq C< 2^{2-l},\\[0.3cm]
\frac{C^2}{2}, \hspace{0.5cm}&\text{if } C\leq d_{\min}^{\infty}.
\end{cases}
\end{align*}
In particular, for $D<4$ the constants are independent of the cardinality of the space $\lvert \mathcal{X} \rvert$, for $D=4$ the dependence is logarithmic and for $D>4$ of order $1/2-2/D$. However, it should be stressed that the upper bound in general depends on the measure's intrinsic dimension, e.g. if a measure is supported on a $D^\prime<D$ dimensional subspace of $[0,1]^D$, then these upper bounds could be stated in terms of $D^\prime$.\\
The constants in the Bernoulli model coincide with those for the Poisson model. For the multinomial model it holds for $ C\leq d_{\min}^{\infty}$ and any $D\geq 1$ that
\[
\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Mult}}(C)\leq 2^{-1}D\sqrt{\mathbb{M}(\mu)}\sum_{x\in \mathcal{X}}\sqrt{\mu(x)}.
\]
For $C>2$, it holds
\begin{align*}
\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Mult}}(C)\leq D
\begin{cases}
\frac{2^{D/2+3}}{1-2^{D/2-2}} \hspace{2cm}&\text{if } D<4,\\[0.3cm]
 (2+8\log_2\left(\vert \mathcal{X}\vert\right)) \hspace{2cm}&\text{if }   D=4,\\[0.3cm]
 2\vert \mathcal{X}\vert^{\frac{1}{2}-\frac{2}{D}} \left( 1+\frac{2^{2+D/2}}{2^{D/2-2}-1}\right), \hspace{2cm}&\text{if } D>4.
\end{cases}
\end{align*}
Otherwise, the constants also coincide with those in the Poisson model. Hence, the main difference between the multinomial and the Poisson model is the fact that for large $C$ the latter contains an additional term which originates from having to estimate the population level total mass intensity, while in the multinomial model, this is assumed to be known. \\
To further extend these bounds to the context of the $(p,C)$-Fr\'echet functionals and $(p,C)$-barycenters \eqref{eq:KRbarycenter}, we utilise certain structural properties of the $(p,C)$-barycenter. Let $\mu^*$ be any $(p,C)$-barycenter of the population measures $\mu^1,\dots,\mu^J$ and let $\hat{\mu}^*$ be any barycenter of the empirical measures $\hat{\mu}^1,\dots,\hat{\mu}^J$. From $(i)$, we can derive a bound on the error in the Fr\'echet functional and in conjugation with theory from linear programming, we then extend this control to a more refined statement on the $(p,C)$-barycenters (see \Cref{thm:kr_boundPoi}). The caveat is the fact that neither $\mu^*$ nor $\hat{\mu}^*$ is necessarily unique. Thus, we need to control the error in terms of the respective optimal set $\mathbf{B}^*$ and its empirical counterpart $\hat{\mathbf{B}}^*$. For this, we consider the quantities
\begin{align*}
    (iii) \ \mathbb{E}\left[F_{p,C}(\hat{\mu}^*)-F_{p,C}({\mu}^*)\right], \quad (iv)  \  \mathbb{E}\left[ \sup_{\hat{\mu}^\star \in \hat{\mathbf{B}}^\star} \inf_{\mu^\star \in \mathbf{B}^\star} KR_{p,C}^p(\mu^\star,\hat{\mu}^\star) \right].
\end{align*}
The term $(iv)$ is the expected deviation of the empirical barycenter which has the largest $(p,C)$-KRD to the set of barycenters of the population measures. The bounds on $(iii)$ and $(iv)$ generalise recent bounds in the context of OT barycenters established in \cite{heinemann2022randomized}. \\
As a major application, we stress that these bounds enable randomised computations of UOT with statistical guarantees. If the size of the population measures is computationally infeasible, then empirical versions of these measures can be used as a proxy for the population level distances and barycenter. The bounds $(i)-(iv)$ allow to tune the problem size against the quality of the approximation. While theoretically all three models allow this approximation approach, the most suitable candidate is clearly the multinomial model. In this resampling scenario, the assumption of known total intensities is natural, as the population measures are known, but solving UOT between them is computationally infeasible. Here, the sample size $N$ provides a strict upper bound on the computational complexity of a given approximation. This is not the case though, for the other two models, where the sampling schemes do not allow to bound the required computational complexity beforehand.\\
One alternative approach, closely related to the multinomial model, is to make use of a subsampling\footnote{Here, subsampling refers to sampling without replacement, while resampling refers to sampling with replacement.} method instead of a resampling one, i.e. to replace the i.i.d. samples $X_1,\dots,X_N$ from $\mu$ by ones drawn without replacement. A natural choice of estimator in this scenario is
\begin{align}\label{eq:subsamp}
    \hat{\mu}_{N}=\frac{1}{\sum_{i=1}^N \mu(X_i)}\sum_{i=1}^N \mu(X_i)\delta_{X_i},
\end{align}
where the mass at each drawn location $x\in \mathcal{X}$ is proportional to the mass of the population measure at $x$ and the total mass intensity is rescaled to the known, true total intensity. This estimator is, due to the sampling without replacement, guaranteed to have $N$ support points, which yields close control on the required runtime for a given approximation. This approach has become popular within the machine learning community where it is referred to as mini-batch OT \citep{fatras2021minibatch,nguyen2021transportation}. An illustration of the potential runtime advantages using the suggested randomised methods is displayed in \Cref{fig:re_vs_sub_time}. For this example, the resampling approach provides an expected relative KRD error of about $5\%$ while achieving a speedup of about a factor of $1000$ compared to the original runtime, while the subsampling approach requires nearly $10\%$ of the original runtime to achieve the same accuracy. A more detailed comparison of the empirical performance of the re- and subsampling model is found in \Cref{sec:realdata}. Though, we note that, in the considered data examples, the resampling approach consistently performed better than the subsampling one. We also study the convergence properties of $(i)-(iv)$ for all three described models in extended simulation studies on a wide range of synthetic datasets. \\
In summary, we provide non-asymptotic deviation bounds for UOT distances and corresponding barycenters. A major finding in this regard is that UOT can be analysed in a variety of statistical models. This extends results for classical OT which has been traditionally only considered in the context of comparing probability distributions. Our results are explicit in their dependency on the cardinality of the measures' support sizes and their total mass intensities. Most importantly, the upper bounds reveal an interesting interaction between the structural properties of the measures and the regularisation parameter $C$. Our results can be used for statistically sound UOT based analysis of data generators with different intensities, as well as for randomised computational schemes with statistical error control.
\begin{figure}
\begin{tabular}{ccc}
      \includegraphics[width=0.32\textwidth]{graphics/IM1.png} &
      \includegraphics[width=0.32\textwidth]{graphics/IM2.png} &
    \includegraphics[width=0.32\textwidth]{graphics/time_comparison_resampling_vs_subsampling.png}
\end{tabular}
\caption{\textbf{Left:} An excerpt of size $300\times 300$ from the STED microscopy data of adult Human Dermal Fibroblasts labelled at MIC60 (a mitochondrial inner membrane complex) in \cite{tameling2021colocalization}. \textbf{Center:} The same type of data labelled at TOM20 (translocase of the outer mitochondrial membrane) in \cite{tameling2021colocalization}. \textbf{Right:} The expected relative $(2,0.1)$-KRD error curves obtained for the two images on the left and in the centre from re- and subsampling for specified runtimes (computed with the CRAN package \emph{WSGeometry}). The orange line indicates an error level of $5\%$. The respective sample sizes are between $100$ and $9000$ for both approaches, while the original images have about $23000$ non-zero pixels.}
  \label{fig:re_vs_sub_time}
\end{figure}

\section{Sampling Bounds for Kantorovich-Rubinstein Distances}\label{chp:sampbounds}
In this section we investigate the Poisson model and its estimator in \eqref{eq:poissonmeasure} in detail. We provide theoretical guarantees for the accuracy of the approximation in terms of their \emph{expected Kantorovich-Rubinstein deviation}. Results for the multinomial and Bernoulli model follow along the same reasoning. Corresponding deviation bounds and the proofs are provided in \Cref{sec:multi} and \Cref{sec:ber}. \\
The proof of the main result and the explicit construction of its constant relies on a tree approximation of the space $\mathcal{X}$. This approximation requires some \emph{depth level} $L\in\mathbb{N}$ and \emph{resolution} $q>1$, based on which we construct minimal $q^{-j}\text{diam}(\mathcal{X})$-coverings\footnote{For a metric space $(\mathcal{X},d)$ an $\epsilon$-cover is a set of points $\left\{ x_1,\ldots,x_m\right\}\subset \mathcal{X}$ such that for each $x\in\mathcal{X}$, there exists some $1\leq i \leq m$ such that $d(x,x_i)\leq \epsilon$. The smallest such set is denoted as $\mathcal{N}(\mathcal{X},d,\epsilon)$.} on $\mathcal{X}$ to obtain our bounds. An illustration of this approximation is given in \Cref{fig:treeapproximation}.
\subsection{Tree Approximation for the Kantorovich-Rubinstein Distance}\label{app:treeapproximation}
Let $\mathcal{T}=(V,E)$ be a rooted, ultrametric tree with height function $h:V\longrightarrow \mathbb{R}_+$ and root $\mathsf{r}$. For two nodes $\textsf{u},\textsf{v}\in V$, denote the unique path between $\mathsf{u}$ and $\mathsf{v}$ in $\mathcal{T}$ by $\mathcal{P}(\mathsf{u},\mathsf{v})$. For a node $\textsf{v}\in V$ its \emph{children} are the elements of the set $\mathcal{C}(\textsf{v})=\left\lbrace \textsf{w}\in V \,\mid\, \textsf{v}\in \mathcal{P}(\textsf{w},\textsf{r})\right\rbrace$. The \emph{parent} $par(\mathsf{v})$ of a node $\textsf{v}$ is the unique node with $(par(\mathsf{v}),\mathsf{v})\in E$ and $h(\mathsf{v})< h(par(\mathsf{v}))$.
For any $C>0$, define the set
\begin{align}\label{eq:subtreeroot}
\mathcal{R}(C)\coloneqq \left\lbrace \textsf{v}\in V\, \mid \, h(\textsf{v})\leq C/2< h(\text{par}(\textsf{v}))\right\rbrace
\end{align}
with the convention that $\mathcal{R}(C)=\lbrace \textsf{r} \rbrace$ if $\frac{C}{2}\geq h(\textsf{r})$. The goal is to control the $(p,C)$-KRD on the finite metric space $(\mathcal{X},d)$ by bounding it from above by a dominating distance $d_\mathcal{T}$ induced\footnote{For two vertices of the tree $\mathcal{T}$, we define their distance $d_{\mathcal{T}}$ as the sum of the weights of the edges included in the unique path between the two vertices. Here, the weight of an edge joining two vertices $v$ and $par(v)$ is given by $h(par(v))-h(v)$.\label{fn:treedist}} from a tree $\mathcal{T}$ with the elements of $\mathcal{X}$ as vertices and a height function $h$ such that $d(x,x^\prime)\leq d_\mathcal{T}(x,x^\prime)$. In this case and by the definition of the Kantorovich-Rubinstein distance it holds for all measures $\mu,\nu\in\mathcal{M}_+(\mathcal{X})$ that
\begin{align}\label{eq:generalKRbound}
\text{KR}_{p,C}(\mu,\nu)\leq \text{KR}_{d_\mathcal{T}^p,C}(\mu,\nu),
\end{align}
where $\text{KR}_{d_\mathcal{T}^p,C}(\mu,\nu)$ denotes the $(p,C)$-KRD w.r.t. the ground space $(\mathcal{X},d_\mathcal{T})$. Moreover, if $\mathcal{T}$ is an ultrametric tree with leaf nodes $L$ and height function $h\colon V\to \mathbb{R}_+$ inducing\cref{fn:treedist} the tree metric $d_\mathcal{T}$ and the two measures $\mu^L,\nu^L\in \mathcal{M}_+(L)$ supported on the leaf nodes of $\mathcal{T}$, then it holds \citep{heinemann2022kantorovich} that
{\begin{align}\label{eq:KRultrametric}
\begin{split}
    &\text{KR}^p_{d_\mathcal{T}^p,C}\left(\mu^L,\nu^L\right)=\\
 &\sum_{\textsf{v}\in \mathcal{R}(C)} \Bigg(2^{p-1}\sum_{\textsf{w}\in \mathcal{C}(\textsf{v})\setminus \lbrace \textsf{v} \rbrace}  \Big(\left( h(par(\textsf{w}))^p-h(\textsf{w})^p \right) \left\vert \mu^L(\mathcal{C}(\textsf{w}))-\nu^L(\mathcal{C}(\textsf{w}))\right\vert \Big) \\[1ex]
&+\left(\frac{C^p}{2}-2^{p-1}h(\textsf{v})^p\right) \left\vert \mu^L(\mathcal{C}(\textsf{v}))-\nu^L(\mathcal{C}(\textsf{v}))\right\vert \Bigg).
\end{split}
\end{align}}
\begin{figure}
  \centering
  \subfloat[][]{\includegraphics[width=0.78\linewidth]{graphics/Treeapproximation1.jpeg}} \hspace{0.2\linewidth}
  \subfloat[][]{\includegraphics[width=0.78\linewidth]{graphics/Treeapproximation2.jpeg}}%
  \caption{\textbf{Ground metric approximation by an ultrametric tree distance:} \textbf{(a)} A finite metric space $(\mathcal{X},d)$ and its covering sets $Q_0$ (black), $Q_1$ (red) and $Q_2$ (green) for $L=2$. \textbf{(b)} Based on the covering sets from (a) an ultrametric tree is constructed. The metric space $\mathcal{X}$ is embedded in level $L+1=3$ and equal to all leaf nodes of that tree.}%
  \label{fig:treeapproximation}
\end{figure}
The construction of $\mathcal{T}$, such that \eqref{eq:generalKRbound} holds, is as follows. Fix some depth level $L\in \mathbb{N}$. For some $q>1$ and level $j=0,\ldots,L$ define the covering set $Q_j\coloneqq \mathcal{N}(\mathcal{X},q^{-j}\text{diam}(\mathcal{X}))\subset \mathcal{X}$ and let $Q_{L+1}\coloneqq\mathcal{X}$. Any point $x\in Q_j$ is considered as a node at level $j$ of a tree $\mathcal{T}$ and denoted as $(x,j)$ to emphasise its level position. For level $j=0$ this yields a single element in $Q_0$ which serves as the root of the tree. For $j=0,\ldots,L$ a node $(x,j)$ at level $j$ is connected to one node $(x^\prime,j+1)$ at level $j+1$ if their distance satisfies $d(x,x^\prime)\leq q^{-j}\text{diam}(\mathcal{X})$ (ties are broken arbitrarily). The edge weight of the corresponding edge is set equal to $q^{-j}\text{diam}(\mathcal{X})$. Consequently, the height of each node only depends on its assigned level $0\leq l\leq L+1$ and is defined as $h_{q,L}\colon \{0,\ldots,L+1\}\to\mathbb{R}$ by
\begin{align}\label{eq:height}
h_{q,L}(l)=\sum_{j=l}^{L} q^{-j}\text{diam}(\mathcal{X})=\frac{q^{1-l}-q^{-L}}{q-1}\text{diam}(\mathcal{X}).
\end{align}
By definition the space $\mathcal{X}$ is embedded in level $L+1$ as the leaf nodes of $\mathcal{T}$ with height $h_{q,L}(L+1)=0$. By a straightforward computation it holds for two points $x,x^\prime \in \mathcal{X}$ considered as embedded in $\mathcal{T}$ as $(x,L+1)$ and $(x^\prime,L+1)$ that
\begin{align}\label{eq:metricupper}
  d^p(x,x^\prime)\leq d_\mathcal{T}^p((x,L+1),(x^\prime,L+1)).
\end{align}
The measures $\mu,\nu$ are embedded into $\mathcal{T}$ as measures $\mu^L,\nu^L$ supported only on leaf nodes of $\mathcal{T}$ and thus it follows from \eqref{eq:metricupper} that
\begin{align*}
\text{KR}_{p,C}(\mu,\nu)\leq \text{KR}_{d_\mathcal{T}^p,C}\left(\mu^L,\nu^L\right).
\end{align*}
In combination with the closed formula from \eqref{eq:KRultrametric} this yields an upper bound on the $(p,C)$-KRD. Whenever clear from the context the notation is alleviated by writing $\textsf{v}\in Q_l$ instead of $(\textsf{v},l) \in Q_l$.

\begin{lemma}\label{cor:treeapproximation}
Let $(\mathcal{X},d)$ be a finite metric space and let $\mu,\nu\in\mathcal{M}_+(\mathcal{X})$ with total mass $\mathbb{M}(\mu)$ and $\mathbb{M}(\nu)$, respectively. Let $p\geq 1$ and $C>0$. Then for any resolution $q>1$, depth $L\in\mathbb{N}$ and height function \eqref{eq:height} with $$h_{q,L}(k)=\frac{q^{1-k}-q^{-L}}{q-1}\text{diam}(\mathcal{X})$$ it holds that
\begin{align*}
&\text{KR}^p_{p,C}(\mu,\nu)\leq
&\begin{cases}
\left(\frac{C^p}{2}-2^{p-1}h_{q,L}(0)^p\right)\left\vert \mathbb{M}(\mu)-\mathbb{M}(\nu)\right\vert\\[0.1cm]+B_{q,p,L,\mathcal{X}}(1), \\[0.2cm] \hspace{2cm}    \text{if } C\geq 2h_{q,L}(0), \\[0.6cm]
B_{q,p,L,\mathcal{X}}(l), \\[0.2cm] \hspace{2cm} \text{if } 2h_{q,L}(l)\leq  C< 2h_{q,L}(l-1),\\[0.6cm]
\frac{C^p}{2} \text{TV}(\mu,\nu), \\[0.2cm] \hspace{2cm} \text{if } C\leq \left( 2h_{q,L}(L) \vee {\min_{x\neq x^\prime}d(x,x^\prime)}\right), \\
\end{cases}
\end{align*}
where $TV(\mu,\nu)=\sum_{x\in \mathcal{X}} \lvert \mu(x)-\nu(x) \rvert$ is the \emph{total variation distance} and
\begin{align*}
    B_{q,p,L,\mathcal{X}}(l)=2^{p-1}\sum\limits_{j=l}^{L+1}\sum\limits_{x\in Q_j} \left( h_{q,L}(j-1)^p-h_{q,L}(j)^p\right)\left\vert \mu^L(\mathcal{C}(x))-\nu^L(\mathcal{C}(x))\right\vert.
\end{align*}
\end{lemma}
The proof is given in \Cref{app:proofs}.
 \subsection{$\mathbf{(p,C)}$-Kantorovich-Rubinstein Deviation Bound}
Given \Cref{cor:treeapproximation} based on the tree-approximation, we are able to state our main result explicitly.
\begin{theorem}\label{thm:samplingboundKR}
Let $(\mathcal{X},d)$ be a finite metric space and $\mu\in\mathcal{M}_+(\mathcal{X})$. Let $\hat{\mu}_{t,s}$ be the estimator from \eqref{eq:poissonmeasure}. Then, for any $p\geq 1$, resolution $q>1$ and depth $L\in \mathbb{N}$ it holds that
\begin{align*}
       \mathbb{E}\left[ \text{KR}_{p,C}\left( \hat{\mu}_{t,s},\mu\right)\right]\leq
       \mathcal{E}_{p,\mathcal{X},\mu}^{\text{Poi}}(C,q,L)^{1/p}\begin{cases}
   \left(2(1-s)\mathbb{M}(\mu)+\frac{s}{\sqrt{t}}\sum_{x\in \mathcal{X}}\sqrt{\mu(x)} \right)^{\frac{1}{p}},  \\[0.2cm]  \hspace{2cm} \text{if } C\leq \left(2h_{q,L}(L) \vee {\min_{x\neq x^\prime}d(x,x^\prime)}\right) \\[0.6cm]
  \left( \frac{1}{st}\mathbb{M}(\mu)+\frac{1-s}{s}\sum_{x\in \mathcal{X}}\mu(x)^2 \right)^{\frac{1}{2p}}. \\[0.2cm]  \hspace{2cm} \text{else}.
    \end{cases}
\end{align*}
For
\[
A_{q,p,L,\mathcal{X}}(l):=\text{diam}(\mathcal{X})^p2^{p-1}\left(q^{-Lp}\lvert  \mathcal{X} \rvert^{\frac{1}{2}}+\left(\frac{q}{q-1}\right)^p \sum_{j=l}^{L}q^{p-jp}\lvert Q_j \rvert^{\frac{1}{2}}\right),
\]
the constant is equal to
\begin{align*}
\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Poi}}(C,q,L)=\begin{cases}
\left(\frac{C^p}{2}-2^{p-1}\left( \frac{q-q^{-L}}{q-1}\text{diam}(\mathcal{X})\right)^p \right)
+A_{q,p,L,\mathcal{X}}(1),\\[0.2cm] \hspace{2cm}\text{if }  C\geq 2h_{q,L}(0),\\[0.6cm]
A_{q,p,L,\mathcal{X}}(l),\\[0.2cm] \hspace{2cm}\text{if }  2h_{q,L}(l)\leq C< 2h_{q,L}(l-1),\\[0.6cm]
\frac{C^p}{2},\\[0.2cm] \hspace{2cm} \text{if }  C\leq \left(2h_{q,L}(L) \vee {\min_{x\neq x^\prime}d(x,x^\prime)}\right),
\end{cases}
\end{align*}
Furthermore, for $p=1$ the factor $\frac{q}{(q-1)}$ in $A_{q,1,L,\mathcal{X}}(l)$ can be removed for all $l=1,\dots,L$.
\end{theorem}

The constant $\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Poi}}(C,q,L)$ is reminiscent of the constants for similar deviation bounds for optimal transport between finitely supported measures \citep{boissard2014mean,sommerfeld2019optimal}. However, in the case of UOT one finds an interesting case distinction into roughly three cases depending on the relation between the penalty parameter $C$ of the $(p,C)$-KRD and the resolution $q$ and depth $L$ of the tree approximation. The different constants arise from the fact that $C$ controls the maximal range at which transport occurs in an UOT plan. In particular, if $d(x,x^\prime)>C^p$, then for any UOT plan $\pi$ it holds $\pi(x,x^\prime)=0$. If $C$ is sufficiently large, i.e. larger than the diameter of $\mathcal{X}$, then $\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Poi}}(C,q,L)$ coincides with the deviation bounds for usual optimal transport, however, there is an additional summand arising from the necessary estimation of the true total mass $\mathbb{M}(\mu)$. For sufficiently small $C$, e.g. $C<\min_{x\neq x^\prime} d(x,x^\prime)$, the $(p,C)$-KRD is proportional to the TV distance, hence we obtain a constant which is oblivious to the geometry of the ground space (see Lemma~\ref{lem:TVbound}). For an intermediate value of $C$, the UOT problem on the ultra-metric tree $\mathcal{T}$ decomposes into smaller problems on subtrees of $\mathcal{T}$ (for details see the proof of \eqref{eq:KRultrametric} in \cite{heinemann2022kantorovich}) depending on $C$. The expected $(p,C)$-KRD error then depends on the size of these subtrees and the mass estimation error inherent to the total mass on these subtrees.
\begin{remark}\label{rem:infimum}
Since the deviation bound holds for any resolution $q>1$ and depth $L\in\mathbb{N}$ one can optimise and equivalently state upper bounds in terms of the infimum over those parameters. When the dependence on $q$ or $L$ is omitted, it is assumed that the infimum over those parameters has been taken, i.e.
\begin{align*}
    \mathcal{E}_{p,\mathcal{X},\mu}^{\text{Poi}}(C)=\inf_{L\in \mathbb{N},q>1} \mathcal{E}_{p,\mathcal{X},\mu}^{\text{Poi}}(C,q,L).
\end{align*}
\end{remark}
From the reverse triangle inequality we immediately obtain the following corollary from \Cref{thm:samplingboundKR}.
\begin{corollary}
Let $(\mathcal{X},d)$ be a finite metric space and $\mu,\nu\in\mathcal{M}_+(\mathcal{X})$. Let $\hat{\mu}_{t,s},\hat{\nu}_{t,s}$ be the estimator from \eqref{eq:poissonmeasure} for each of these measures, respectively. Then, for any $p\geq 1$, resolution $q>1$ and depth $L\in \mathbb{N}$ it holds that
\begin{align*}
       &\mathbb{E}\left[ \lvert \text{KR}_{p,C}\left( \hat{\mu}_{t,s},\hat{\nu}_{t,s}\right) -  \text{KR}_{p,C}\left( {\mu},{\nu}\right)\rvert \right]\leq\\
       &2\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Poi}}(C,q,L)^{1/p}\begin{cases}
   \left(2(1-s)\mathbb{M}(\mu)+\frac{s}{\sqrt{t}}\sum_{x\in \mathcal{X}}\sqrt{\mu(x)} \right)^{\frac{1}{2p}} , &\text{if } C< \min_{x,x^\prime\in\mathcal{X}}d(x,x^\prime), \\
  \left( \frac{1}{st}\mathbb{M}(\mu)+\frac{1-s}{s}\sum_{x\in \mathcal{X}}\mu(x)^2 \right)^{\frac{1}{2p}}, &\text{else}, \\
    \end{cases}
\end{align*}
where $\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Poi}}(C,q,L)$ is defined as in \Cref{thm:samplingboundKR}.
\end{corollary}
\subsection{Rate Optimality}
In the following, we provide some intuition and consequences of the deviation bound provided in \Cref{thm:samplingboundKR}. The term $\phi(s,t)$ must necessarily contain a sum of terms depending on $s$ and $t$, respectively. In particular, neither choosing $s=1$ nor letting $t$ go to infinity for $s<1$, would yield a zero error. For any fixed $t >0$ and $s=1$ the expected $(p,C)$-KRD error is clearly non-zero as the mass of the measures at each location is in general not estimated correctly. Similarly, for any fixed $s<1$ letting $t\rightarrow \infty$ can not yield an expected $(p,C)$-KRD error of zero as on average $(1-s)\lvert \mathcal{X} \rvert>0$ support points of $\mu$ are not observed. However, for $s=1$ the error vanishes for $t\rightarrow \infty$, as we observe all support points of $\mu$ and then the strong law of large numbers guarantees the convergence of the weights at each location. It remains to verify whether the rate in $t$ is optimal. For this, fix $s=1$ and observe that
\begin{align}\label{eq:naivetv}
  \min\{C,\min_{x\neq x^\prime}d(x,x^\prime)\}^p TV(\mu,\nu)\leq  KR_{p,C}^p(\mu,\nu)\leq \min\{C,\mathrm {diam}(\mathcal{X})\}^p TV(\mu,\nu).
\end{align}
Hence, it suffices to show that $t^{-\frac{1}{2}}$ is the optimal rate for the convergence in total variation distance. It holds
\[
\mathbb{E}\left[ TV(\mu,\hat{\mu}_{t,1})\right]=\sum_{x\in\mathcal{X}}\mathbb{E} \left[ \lvert \mu(x)-\hat{\mu}_{t,1}(x) \vert\right]=\sum_{x\in\mathcal{X}}t^{-1}\mathbb{E} \left[ \lvert \mathbb{E}[P_x] - P_x \vert\right],
\]
where $P_x\sim Poi(\mu(x)t)$. Using the closed form solution of the mean absolute deviation of a Poisson random variable \citep{ramasubban1958mean} and Stirling's formula, we obtain the asymptotic equivalence (in the sense that their ratio converges to one as $t \rightarrow \infty$)
\[
\mathbb{E}\left[ TV(\mu,\hat{\mu}_{t,1})\right]\approx \sum_{x \in \mathcal{X}}(t\mu(x))^{t\mu(x)-\lfloor t\mu (x) \rfloor}t^{-1/2}\sqrt{\mu(x)}\sqrt{2/\pi}.
\]
For $t>\inf_{x \in \text{supp}(\mu)}\mu(x)$, this is bounded from below by
\[
\left(\sqrt{2/\pi}\sum_{x\in \mathcal{X}}\sqrt{\mu(x)} \right)t^{-1/2}.
\]
All combined, the expectation $\mathbb{E}\left[KR_{p,C}^p(\mu,\hat{\mu}_{t,1}) \right]$ has a lower bound which is asymptotically equivalent to $t^{-\frac{1}{2}}$ and hence the rate in $t$ of \Cref{thm:samplingboundKR} is sharp.

\subsection{Explicit Bounds for Euclidean Spaces}\label{sec:explicitB}
While the constants in the previous theorem are valid on arbitrary metric spaces, more explicit bounds can be derived for many practical applications. Thus, assume that $\mathcal{X}\subset \mathbb{R}^D$ and that $d$ is the Euclidean distance $d_2$. We fix $q=2$, since it minimises the constants from \Cref{thm:samplingboundKR} for all choices of $q$ on the interval $[2,\infty)$. Let $d_\infty (x,x^\prime)=\max_{d=1,\dots, D}\lvert x_d-x^\prime_d \rvert$ be the uniform distance and $\text{diam}_\infty(\mathcal{X})$ be the diameter of $\mathcal{X}$. Using the fact that
\begin{align*}
    d_2(x,x^\prime)\leq \sqrt{D}d_\infty (x,x^\prime),
\end{align*}
we bound the $(p,C)$-KRD with respect to the Euclidean distance by the $(p,C)$-KRD with respect to the uniform distance at a price of a factor of $\sqrt{D}$. In particular, we also apply this to the definition of the height function \eqref{eq:height} given for $q=2$ and any $l\in\{0,\ldots,L\}$ by
\begin{align*}
h_{L}(l)=\left(2^{1-l}-2^{-L}\right)\text{diam}_\infty(\mathcal{X}).
\end{align*}
Within this framework, we can compute explicit upper bounds on the constants in \Cref{thm:samplingboundKR}. For $D<2p$ and $L\rightarrow \infty$, it holds
\begin{align*}
\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Poi}}(C)\leq D^{p/2}
\begin{cases}
    \frac{C^p}{2\sqrt{D}}-2^{2p-1}\text{diam}_\infty^p(\mathcal{X})+\text{diam}^p_\infty(\mathcal{X})2^{3p-1} \frac{2^{D/2-p}}{1-2^{D/2-p}}, \\[0.2cm] \hspace{2cm} \text{if } C\geq 2h_{L}(0),\\[0.6cm]
\text{diam}^p_\infty(\mathcal{X})2^{3p-1}  \frac{2^{l(D/2-p)}}{1-2^{D/2-p}}, \\[0.2cm] \hspace{2cm}\text{if }  2h_{L}(l)\leq C< 2h_{L}(l-1),\\[0.6cm]
\frac{C^p}{2}, \\[0.2cm] \hspace{2cm}\text{if } C\leq \left( 2h_{L}(L) \wedge {\min_{x\neq x^\prime}d_\infty(x,x^\prime)}\right).
\end{cases}
\end{align*}
For $D=2p$ and $L=\lfloor\frac{1}{D}\log_2\left(\vert \mathcal{X}\vert\right) \rfloor$, it holds
\begin{align*}
\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Poi}}(C)\leq D^{p/2}
\begin{cases}
    \frac{C^p}{2\sqrt{D}}-2^{p-1}\left(2-\lvert \mathcal{X} \rvert^{-\frac{1}{D}} \right)^p\text{diam}_\infty^p(\mathcal{X})\\[0.1cm]+\text{diam}^p_\infty(\mathcal{X})2^{3p-1}(2^{-2p}+D^{-1}\log_2\left(\vert \mathcal{X}\vert\right)), \\[0.2cm] \hspace{2cm}\text{if }  C\geq 2h_{L}(0),\\[0.6cm]
\text{diam}^p_\infty(\mathcal{X})2^{3p-1}\left(2^{-2p}+D^{-1}\log_2\left(\vert \mathcal{X}\vert\right)-l\right), \\[0.2cm] \hspace{2cm}\text{if }  2h_{L}(l)\leq  C< 2h_{L}(l-1),\\[0.6cm]
\frac{C^p}{2} ,\\[0.2cm] \hspace{2cm}\text{if }  C\leq \left( 2h_{L}(L) \wedge {\min_{x\neq x^\prime}d_\infty(x,x^\prime)}\right).
\end{cases}
\end{align*}
For $D>2p$ and $L=\lfloor\frac{1}{D}\log_2\left(\vert \mathcal{X}\vert\right) \rfloor$, it holds
\begin{align*}
\mathcal{E}_{p,\mathcal{X},\mu}^{\text{Poi}}(C)\leq D^{p/2}
\begin{cases}
     \frac{C^p}{2\sqrt{D}}-2^{p-1}\left(2-\lvert \mathcal{X} \rvert^{-\frac{1}{D}} \right)^p\text{diam}_\infty^p(\mathcal{X})\\[0.1cm]+\text{diam}^p_\infty(\mathcal{X})2^{p-1}\vert \mathcal{X}\vert^{\frac{1}{2}-\frac{p}{D}} \left( 1+\frac{2^{p+D/2}}{2^{D/2-p}-1}\right),\\[0.2cm] \hspace{2cm}\text{if } C\geq 2h_{L}(0),\\[0.6cm]
\text{diam}^p_\infty(\mathcal{X})\bigg( \vert \mathcal{X}\vert^{\frac{1}{2}-\frac{p}{D}}\\
+\frac{2^{p+D/2}}{2^{D/2-p}-1}\left( \lvert \mathcal{X}\rvert^{\frac{1}{2}-\frac{p}{D}}-2^{(D/2-p)(l-1)}
\right)\bigg),\\[0.2cm] \hspace{2cm}\text{if } 2h_{L}(l)\leq C< 2h_{L}(l-1),\\[0.6cm]
\frac{C^p}{2},\\[0.2cm] \hspace{2cm}\text{if } C\leq \left( 2h_{L}(L) \wedge {\min_{x\neq x^\prime}d_\infty(x,x^\prime)}\right).
\end{cases}
\end{align*}
These bounds distinguish cases, based on the value of the penalty $C$, as well as the dimension $D$. The most critical part of these these bounds is the impact of the cardinality of $\mathcal{X}$. If $D<2p$, then there is no dependence on $\lvert \mathcal{X} \rvert$ and the convergence of approximation error of the empirical measure is independent of its support size. If $D=2p$, then $\lvert \mathcal{X} \rvert$ enters through a logarithmic term. If $D>2p$, then the dependence becomes polynomial in $\lvert \mathcal{X} \rvert$. These phase transitions for the dependence on the support size match those for empirical optimal transport \citep{sommerfeld2019optimal}. A novelty for the UOT setting is the additional dependency on the different scales of $C$. This is explained by the previously discussed control of $C$ on the maximal distance at which transport occurs in an optimal plan. The height function is again used to specify the scale induced by a particular choice of the parameter $C$. Notably, the dependence on $\lvert \mathcal{X} \rvert$ does not change on most scales of $C$. There is an exception, however, for sufficiently small values of $C$, where the $(p,C)$-KRD is equal to a scaled total variation distance. Thus, these bounds are completely oblivious to the geometry of $\mathcal{X}$ in $\mathcal{Y}$, though they scale as $\lvert \mathcal{X} \rvert^{\frac{1}{2}}$ which is the same rate we obtain for the TV bound (recall \cref{lem:TVbound}). As a final observation, we note that for $C>2h_L(0)$, these bounds essentially recover analogue bounds for empirical optimal transport \citep{sommerfeld2019optimal}. However, for the $(p,C)$-KRD the bounds include an additional summand based on the estimation error for the measure's total mass intensity.
\subsection{Proofs}
As a preparatory step to prove \Cref{thm:samplingboundKR}, we treat the significantly simpler case of an empirical deviation bound with respect to the total variation distance.
\begin{lemma}[Total Variation Bound]\label{lem:TVbound}
Let $(\mathcal{X},d)$ be a finite metric space and $\mu\in\mathcal{M}_+(\mathcal{X})$ with total mass $\mathbb{M}(\mu)$. Let $\hat{\mu}_{t,s}$ be the estimator from \eqref{eq:poissonmeasure}. Then, for any $p\geq 1$ it holds that
\begin{align*}
\mathbb{E}\left[ \text{KR}^p_{p,C}\left(\hat{\mu}_{t,s},\mu\right)\right]\leq  \frac{C^p}{2}\left(2(1-s)\mathbb{M}(\mu)+\frac{s}{\sqrt{t}}\sum_{x\in \mathcal{X}}\sqrt{\mu(x)} \right).
\end{align*}
\end{lemma}
\begin{proof}
Let $\mu,\nu\in \mathcal{M}_+(\mathcal{X})$. By retaining all common mass between $\mu$ and $\nu$ at place and delete (resp. create) excess mass (resp. deficient mass) we obtain a feasible solution for \eqref{eq:krdistance} with objective value in terms of a total variation distance between $\mu$ and $\nu$. Thus, it holds
\begin{align*}
\text{KR}^p_{p,C}\left(\mu,\nu\right)\leq \frac{C^p}{2} \text{TV}\left(\mu,\nu\right).
\end{align*}
In particular, this holds for $\nu=\hat{\mu}_{t,s}$. Taking expectations yields
\begin{align*}
    \mathbb{E}\left[\text{TV}\left(\hat{\mu}_{t,s},\mu\right)\right]&=\frac{1}{st}\sum_{x\in \mathcal{X}}\mathbb{E}\left[ \lvert P_xB_x - st\mu(x) \rvert \right] \\
    &=\frac{1}{st}\sum_{x\in \mathcal{X}}s\mathbb{E}\left[ \lvert P_x - st\mu(x) \rvert \right]+(1-s)st\mu(x)\\
    &\leq \frac{1}{st}\sum_{x\in \mathcal{X}}s(1-s)\mathbb{E}\left[ P_x \right] +s^2\mathbb{E}\left[  \lvert P_x-t\mu(x) \rvert \right]+(1-s)st\mu(x)\\
    &\leq \frac{1}{st}\sum_{x\in \mathcal{X}} 2s(1-s)t\mu_x +s^2\sqrt{t}\sqrt{\mu(x)}\\
    &=2(1-s)\mathbb{M}(\mu)+\frac{s}{\sqrt{t}}\sum_{x\in \mathcal{X}} \sqrt{\mu(x)}.
\end{align*}
\end{proof}

With Lemma~\ref{cor:treeapproximation} and Lemma~\ref{lem:TVbound} at our disposal we are able to prove \Cref{thm:samplingboundKR}.

\begin{proof}[Proof of \Cref{thm:samplingboundKR}]
Let $\hat{\mu}_{t,s}$ be the estimator from \eqref{eq:poissonmeasure}. We fix $p=1$ and detail the case $p>1$ at the end of the proof. Suppose first that $C\leq \min_{x\neq x^\prime} d(x,x^\prime)$. According to \cite[Theorem 2.2 (ii)]{heinemann2022kantorovich} it holds that
\begin{align*}
\mathbb{E}\left[\text{KR}_{1,C}(\hat{\mu}_{t,s},\mu)\right]=\frac{C}{2}\mathbb{E}\left[\sum_{x\in \mathcal{X}} \left\vert \hat{\mu}_{t,s}-\mu(x) \right\vert\right]=\frac{C}{2}\mathbb{E}\left[\text{TV}(\hat{\mu}_{t,s},\mu)\right].
\end{align*}
This yields the total variation bounds (see Lemma~\ref{lem:TVbound}). Next, consider the tree approximation as outlined in \Cref{app:treeapproximation} and construct an ultrametric tree $\mathcal{T}$ such that $\text{KR}_{1,C}(\hat{\mu}_{t,s},\mu)\leq \text{KR}_{d_\mathcal{T},C}\left(\hat{\mu}_{t,s}^L,\mu^L\right).$
Applying Lemma~\ref{cor:treeapproximation} for $p=1$ where by definition the difference of height function is equal to
\begin{align*}
h_{q,L}(j-1)-h_{q,L}(j)=\frac{\text{diam}(\mathcal{X})}{q-1} \left( q^{2-j}-q^{1-j}\right)=\text{diam}(\mathcal{X})q^{1-j}
\end{align*}
and yields the upper bound
\begin{align*}
&\mathbb{E}\left[\text{KR}_{1,C}(\hat{\mu}_{t,s},\mu)\right]\leq \mathbb{E}\left[\text{KR}_{d_\mathcal{T},C}\left(\hat{\mu}_{t,s}^L,\mu^L\right)\right]\\[3ex]
&=
\begin{cases}
\left( \frac{C}{2}-h_{q,L}(0)\right)\mathbb{E}\left[\lvert \mathbb{M}(\hat{\mu}_{t,s})-\mathbb{M}(\mu)\rvert\right]\\
+\text{diam}(\mathcal{X})\sum\limits_{j=1}^{L+1}q^{1-j}\sum\limits_{x\in Q_j}\mathbb{E}\left[\left\vert \hat{\mu}_{t,s}^L(\mathcal{C}(x))-\mu^L(\mathcal{C}(x))\right\vert\right], &C\geq 2h_{q,L}(0) \\[3ex]
\text{diam}(\mathcal{X})\sum\limits_{j=l}^{L+1}q^{1-j}\sum\limits_{x^\prime\in Q_j} \mathbb{E}\left[\left\vert \hat{\mu}_{t,s}^L(\mathcal{C}(x^\prime))-\mu^L(\mathcal{C}(x^\prime))\right\vert\right], \,  &2h_{q,L}(l)\leq C< 2h_{q,L}(l-1),\\[3ex]
\frac{C^p}{2} \mathbb{E}\left[ \text{TV}(\hat{\mu}_{t,s},\mu)\right],\,  &C\leq \left( 2h_{q,L}(L) \vee {\min_{x\neq x^\prime}d(x,x^\prime)}\right). \\
\end{cases}
\end{align*}
For the estimator from \eqref{eq:poissonmeasure} with $B_x\sim \text{Ber}(s)$ and $P_x\sim \text{Poi}(t\mu(x))$ for all $x \in \mathcal{X}$, it holds
\begin{align*}
&\sum_{x\in Q_l}\mathbb{E}\left[\left\vert \hat{\mu}_{t,s}^L(\mathcal{C}(x))-\mu^L(\mathcal{C}(x))\right\vert\right]\\&= \sum_{x\in Q_l}\
    \frac{1}{st}\mathbb{E}\left[ \left\lvert \sum_{y\in \mathcal{C}(x)} P_yB_y - st\sum_{y\in \mathcal{C}(x)} \mu(y) \right\rvert \right] \\
    &\leq \sum_{x\in Q_l}\ \frac{1}{st} \sqrt{\text{Var}\left(\sum_{y\in \mathcal{C}(x)} P_yB_y\right)}
    \end{align*}
    \begin{align*}
    &= \sum_{x\in Q_l}\ \frac{1}{st} \sqrt{\sum_{y\in \mathcal{C}(x)}\text{Var}(P_yB_y)}\\
    &=\sum_{x\in Q_l}\ \frac{1}{st} \sqrt{\sum_{y\in \mathcal{C}(x)} s(1-s)t\mu(y)+s(1-s)\mu(y)^2+t\mu(y)s^2}\\
    &= \sum_{x\in Q_l}\ \sqrt{ \frac{1-s}{st}\sum_{y\in \mathcal{C}(x)}\mu(y)+ \frac{1-s}{s}\sum_{y\in \mathcal{C}(x)}\mu(y)^2+\frac{1}{t}\sum_{y\in \mathcal{C}(x)}\mu(y)}\\
        &=\sum_{x\in Q_l} \sqrt{ \frac{1}{st}\mu\left(\mathcal{C}(x)\right)+ \frac{1-s}{s}\sum_{y\in \mathcal{C}(x)}\mu(y)^2}\\
        &\leq \sqrt{\lvert Q_l \rvert}\sqrt{\frac{1}{st} \sum_{x\in Q_l}\mu(\mathcal{C}(x))+\frac{1-s}{s} \sum_{x\in Q_l} \sum_{y \in \mathcal{C}(x)}\mu(y)^2   }\\
        &= \sqrt{\lvert Q_l \rvert}\sqrt{\frac{1}{st} \mathbb{M}(\mu)+\frac{1-s}{s} \sum_{x\in \mathcal{X}} \mu(x)^2   }
\end{align*}
Following an analogous computation one bounds the estimation error for the total mass intensity as
\begin{align*}
    \mathbb{E}\left[ \lvert \mathbb{M}(\hat{\mu}_{t,s})-\mathbb{M}(\mu) \rvert \right] \leq \sqrt{\text{Var}(\mathbb{M}(\hat{\mu}_{t,s}))}\leq \sqrt{\frac{1}{st} \mathbb{M}(\mu)+\frac{1-s}{s} \sum_{x\in \mathcal{X}} \mu(x)^2   }.
\end{align*}
Applying both of these bounds to the previous upper bound on the $(p,C)$-KRD in  Lemma~\ref{cor:treeapproximation} yields the claim. \\
For $p>1$, we first observe again that if $C\leq \min_{x\neq x^\prime} d(x,x^\prime)$ then according to \cite{heinemann2022kantorovich} it holds that
\begin{align*}
\mathbb{E}\left[\text{KR}^p_{p,C}(\hat{\mu}_{t,s},\mu)\right]=\frac{C^p}{2}\mathbb{E}\left[\text{TV}(\hat{\mu}_{t,s},\mu)\right]
\end{align*}
which yields the total variation bounds. For more general $C$, we simple repeat the previous calculations with the upper bounds on the difference of height function $h_{q,L}(j-1)^p-h_{q,L}(j)^p\leq \text{diam}(\mathcal{X})^p \left(\frac{q}{q-1}\right)^p q^{p-jp}$.
Since $h_{q,L}(L+1)=0$ we also have $h_{q,L}(L)^p-h_{q,L}(L+1)^p= \text{diam}(\mathcal{X})^pq^{-Lp}$. The expectations are bounded identically as before. Finally, using Jensen's inequality to bound
\[
\mathbb{E}\left[ KR_{p,C}(\mu,\hat{\mu}_{t,s}) \right]\leq \left( \mathbb{E}\left[ KR_{p,C}^p(\mu,\hat{\mu}_{t,s}) \right]\right)^{\frac{1}{p}}
\]
finishes the proof.
\end{proof}

\begin{remark} \label{rem:krdp}
Omitting Jensen's inequality in the last step of the proof of \Cref{thm:samplingboundKR} implies the slightly stronger result
\begin{align*}
       \mathbb{E}\left[ \text{KR}^p_{p,C}\left( \hat{\mu}_{t,s},\mu\right)\right]\leq
       \mathcal{E}_{p,\mathcal{X},\mu}^{\text{Poi}}(C,q,L)\begin{cases}
   \left(2(1-s)\mathbb{M}(\mu)+\frac{s}{\sqrt{t}}\sum_{x\in \mathcal{X}}\sqrt{\mu(x)} \right) , \\[0.2cm] \hspace{2cm} \text{if } C\leq \left(2h_{q,L}(L) \vee {\min_{x\neq x^\prime}d(x,x^\prime)}\right), \\[0.6cm]
  \left( \frac{1}{st}\mathbb{M}(\mu)+\frac{1-s}{s}\sum_{x\in \mathcal{X}}\mu(x)^2 \right)^{\frac{1}{2}}, \\[0.2cm] \hspace{2cm} \text{else}. \\
    \end{cases}
\end{align*}
\end{remark}

\section{Empirical Kantorovich-Rubinstein Barycenters}
 Consider measures $\mu^1,\ldots,\mu^J\in\mathcal{M}_+(\mathcal{X})$ which we replace by $\hat{\mu}^1_{t_1,s_1},\ldots,\hat{\mu}^J_{t_J,s_J}\in\mathcal{M}_+(\mathcal{X})$ as defined in \eqref{eq:poissonmeasure}. We again focus on the Poisson model and treat the remaining two models in the appendix. The previous upper bound on the Kantorovich-Rubinstein distance in \Cref{thm:samplingboundKR} between a measure and its empirical version is used to achieve control on the mean absolute deviation of $(p,C)$-barycenters in terms of their $p$-Fr\'echet functional $F_{p,C}(\mu)=\frac{1}{J}\sum_{i=1}^J \text{KR}_{p,C}^p(\mu^i,\mu)$ from \eqref{eq:KRbarycenter}. We denote
\begin{align*}
\mu^\star\in \argmin_{\mu\in\mathcal{M}_+(\mathcal{Y})} \frac{1}{J} \sum_{i=1}^J \text{KR}_{p,C}^p(\mu^i,\mu),\quad \hat{\mu}^\star\in \argmin_{\mu\in\mathcal{M}_+(\mathcal{Y})} \frac{1}{J} \sum_{i=1}^J \text{KR}_{p,C}^p(\hat{\mu}_{t_i,s_i}^i,\mu)
\end{align*}
and measure the accuracy of approximation of $\mu^\star$ by $\hat{\mu}^\star$ in terms of their mean absolute $p$-Fr\'echet deviation.

\begin{theorem}\label{thm:frechetboundPoi}
Let $\mu^1,\ldots,\mu^J\in\mathcal{M}_+(\mathcal{X})$ and denote $\mathcal{X}_i=\text{supp}(\mu^i)$ for $i=1,\dots ,J$. Consider random estimators $\hat{\mu}_{t_1,s_1}^1,\ldots,\hat{\mu}_{t_J,s_J}^J\in\mathcal{M}_+(\mathcal{X})$ derived from \eqref{eq:poissonmeasure}. Then,
\begin{align*}
\mathbb{E}\left[\lvert F_{p,C}(\mu^\star)-F_{p,C}(\hat{\mu}^\star) \rvert \right] \leq \frac{2p \min \{\mathrm {diam}(\mathcal{Y}),C\}^{p-1}}{J}\easysum{i=1}{J} \mathcal{E}_{1,\mathcal{X}_i,\mu^i}^{\text{Poi}}(C)\phi(t_i,s_i),
\end{align*}
where $\phi$ is given by
\begin{align*}
    \phi(t,s)=\begin{cases}
   \left(2(1-s)\mathbb{M}(\mu)+\frac{s}{\sqrt{t}}\sum_{x\in \mathcal{X}}\sqrt{\mu(x)} \right), & C\leq \min_{x\neq x^\prime}d(x,x^\prime) \\
  \left( \frac{1}{st}\mathbb{M}(\mu)+\frac{1-s}{s}\sum_{x\in \mathcal{X}}\mu(x)^2 \right)^{\frac{1}{2}}, &\text{else}. \\
    \end{cases}
\end{align*}
\end{theorem}

A more elaborate statement gives control over the set of empirical $(p,C)$-barycenters itself. This involves a related linear program that is presented in detail in \Cref{app:sec_LP}.
\begin{theorem}\label{thm:kr_boundPoi}
Let $\mu^1,\ldots,\mu^J\in\mathcal{M}_+(\mathcal{X})$ and denote $\mathcal{X}_i=\text{supp}(\mu^i)$ for $i=1,\dots ,J$. Consider random estimators $\hat{\mu}_{t_1,s_1}^1,\ldots,\hat{\mu}_{t_J,s_J}^J\in\mathcal{M}_+(\mathcal{X})$ derived from \eqref{eq:poissonmeasure}. Let $\mathbf{B}^\star$ be the set of $(p,C)$-barycenters of $\mu^1,\dots,\mu^J$ and $\hat{\mathbf{B}}^\star$ the set of $(p,C)$-barycenters of $\hat{\mu}^1_{t_1,s_1},\ldots,\hat{\mu}^J_{t_J,s_J}$. Then, for $p\geq 1$ it holds that
\begin{align*}
\mathbb{E}\left[ \sup_{\hat{\mu}^\star \in \hat{\mathbf{B}}^\star} \inf_{\mu^\star \in \mathbf{B}^\star} KR_{p,C}^p(\mu^\star,\hat{\mu}^\star) \right]
\leq\frac{p \min \{\mathrm {diam}(\mathcal{Y}),C\}^{p-1} }{V_{P}J}\easysum{i=1}{J} \mathcal{E}_{1,\mathcal{X}_i,\mu^i}^{\text{Poi}}(C)\phi(t_i,s_i),
\end{align*}
where $\phi$ is defined as in \Cref{thm:frechetboundPoi}. The constant $V_P$ is strictly positive and given by
\begin{align*}
V_{P}\coloneqq V_{P}(\mu^1,\dots,\mu^J):=(J+1)\mathrm {diam}(\mathcal{X})^{-p}\underset{v \in V\backslash V^\star}{\min} \ \frac{c^Tv-f^*}{d_1(v,\mathcal{M})},
\end{align*}
where $V$ is the set of feasible vertices from the linear program in \Cref{app:sec_LP}, $V^\star$ is the subset of optimal vertices, $c$ is the cost vector of the program, $f^\star$ is the optimal value, $\mathcal{M}$ is the set of minimisers of the linear program \eqref{eq:lpbary} and $d_1(x,\mathcal{M})=\inf_{y\in \mathcal{M}} \lVert x-y \rVert_1$.
\end{theorem}
The proofs of \Cref{thm:frechetboundPoi} and \Cref{thm:kr_boundPoi} are deferred to \Cref{app:proofs}.
\begin{remark}
For $J=1$ and any $p\geq 1, C>0$ the $(p,C)$-barycenter of $\mu^1$ is just $\mu^1$. Thus, the optimal value of the Fr\'echet functional is zero and it holds
\[
F(\hat{\mu}^*)-F(\mu^*)=KR_{p,C}^p(\mu^1,\hat{\mu}^1_{t,s}).
\]
Consequently, it also holds
\[
\sup_{\hat{\mu}^\star \in \hat{\mathbf{B}}^\star} \inf_{\mu^\star \in \mathbf{B}^\star} KR_{p,C}^p(\mu^*,\hat{\mu}^*)=KR_{p,C}^p(\mu^1,\hat{\mu}^1_{t,s}).
\]
Thus, the rate for the convergence of the $(p,C)$-barycenter of the empirical measures, can in general not be faster than the convergence rate of a single estimator. In particular, the rates in $t$ in \Cref{thm:frechetboundPoi} and \Cref{thm:kr_boundPoi} are sharp.
\end{remark}

\section{Simulations}\label{sec:sims}
In this section we investigate empirically the decay in the expected error for the Poisson model for measures within $\mathcal{X}\subset [0,1]^2$. For the $(p,C)$-KRD we consider two measures $\mu,\nu \in \mathcal{M}_+(\mathcal{X})$ and the \emph{relative $(p,C)$-KRD error}\footnote{We define $0/0:=0$.\label{fn0}}
\begin{align}\label{eq:relKR}
    \mathbb{E}\left[\left\lvert \frac{KR_{p,C}(\hat{\mu}_{t,s},\hat{\nu}_{t,s})-KR_{p,C}(\mu,\nu)}{KR_{p,C}(\mu,\nu)} \right\rvert\right].
\end{align}
For the setting of barycenters we consider the \emph{relative $(p,C)$-Fr\'echet error}\cref{fn0}
\begin{align}\label{eq:relFre}
    \mathbb{E}\left[ \frac{F_{p,C}(\hat{\mu}^*)-F_{p,C}(\mu^*)}{F_{p,C}(\mu^*)} \right].
\end{align}
In both cases, the relative error allows for easier comparisons between models than the absolute error. Additionally, for the $(p,C)$-barycenter \eqref{eq:relFre} is readily available from simulations, while numerically considering the quantity in \Cref{thm:kr_boundPoi} is difficult, as it requires all optimal solutions instead of a single one. All computations of the KRD and the $(p,C)$-barycenter in this section are performed using the methods available in the CRAN package \emph{WSGeometry}.
\subsection{Synthetic Datasets}
We consider eight types of measures for our simulations. Let us fix some notation. Let $J\in \mathbb{N}$ be the number of measures generated. Let $U[0,1]^2$ denote the uniform distribution and let $\text{Poi}(\lambda)$ denote a Poisson distribution with intensity $\lambda$. In all eight settings considered below, the measures are of the form
\[
\mu^i=\sum_{k=1}^{K_i}w^i_k\delta_{l^i_k}
\]
for some weights $w^i_k$, locations $l^i_k$ and $K_i\in \mathbb{N}$. If $K_i=K_j$ for all $i,j=1,\dots ,J$, then we omit the index and denote the number of points by $K$. Note, that all measures have been constructed to have their support included in $[0,1]^2$.

\subsubsection*{Poisson Intensities on Uniform Positions (PI), see \Cref{fig:poissondata} (a)}
Let $w^i_1,\dots ,w^i_K\sim \text{Poi}(\lambda)$ for some intensity $\lambda>0$ and $l^i_1,\dots,l^i_K \sim U[0,1]^2$ for $1\leq i\leq J$.
\subsubsection*{Poisson Intensities on a Grid (PIG), see \Cref{fig:poissondata} (b)}
Set $K=M^2$ for $M\in \mathbb{N}$ and let $w^i_1,\dots ,w^i_{M^2}\sim \text{Poi}(\lambda)$ and $l^i_1,\dots,l^i_{M^2}$ be the location of an equidistant $M\times M$ grid in $[0,1]^2$ for $1\leq i\leq J$.
\begin{figure}[H]
  \centering
  \subfloat[][]{\includegraphics[width=0.45\linewidth]{graphics/PI_data.png}}%
  \qquad
  \subfloat[][]{\includegraphics[width=0.45\linewidth]{graphics/PIG_data.png}}%
  \caption{\textbf{(a)} An example of $J=8$ measures from the \emph{PI} dataset with $M=500$ and $\lambda=5$. \textbf{(b)} An example of $J=8$ measures from the \emph{PIG} dataset with $M=23^2$ and $\lambda=5$.}%
  \label{fig:poissondata}
\end{figure}
\subsubsection*{Norm-based Intensities on Uniform Positions (NI), see \Cref{fig:normeddata} (a)}
Fix $J$ locations $l^1_0,\dots,l^J_0 \in [0,1]^2$. Let $l^i_1,\dots,l^i_K \sim U[0,1]^2$ and let $w^i_k=\lVert l^i_k-l_0^i \rVert_2$ for $1\leq i\leq J$.
\subsubsection*{Norm-based Intensities on a Grid (NIG), see \Cref{fig:normeddata} (b)}
Let $K=M^2$ for $M\in \mathbb{N}$. Fix $J$ locations $l^1_0,\dots,l^J_0 \in [0,1]^2$ and let $l^i_1,\dots,l^i_{M^2}$ be the location of an equidistant $M\times M$ grid in $[0,1]^2$ for each $1\leq i\leq J$. Set $w^i_k=\lVert l^i_k-l_0^i \rVert_2$ for $1\leq i\leq J$.
\begin{figure}[H]
  \centering
  \subfloat[][]{\includegraphics[width=0.45\linewidth]{graphics/NI_data.png}}%
  \qquad
  \subfloat[][]{\includegraphics[width=0.45\linewidth]{graphics/NIG_data.png}}%
  \caption{\textbf{(a)} An example of $J=8$ measures from the \emph{NI} dataset with $M=500$. \textbf{(b)} An example of $J=8$ measures from the \emph{NIG} dataset with $M=23^2$.}%
  \label{fig:normeddata}
\end{figure}
\subsubsection*{Nested Ellipses (NE), see \Cref{fig:ellipsedata} (a)}
Let $G_1,\dots ,G_J\sim U\{1,2,3,4,5\}$ and let $K_i=MG_i$ for $M\in \mathbb{N}$. Set all $w_k^i$ equal to $1$ for each $1\leq k \leq K_i$ for $i=1,\dots,J$. Let $t_1,\dots ,t_M$ be a discretisation of $[0,2\pi]$. Let $U^i_1,\dots,U^i_K,V^i_1,\dots,V^i_K\sim U[0.2,1]$. For $1\leq i\leq J$, set
\[
l^i_{M(j_i-1)+k}=0.5(1+3^{-j}(U_{M(j_i-1)+k}\sin(t_k),V_{M(j_i-1)+k}\cos(t_k))^T), \quad j_i=0,\dots ,G_i.
\]
\subsubsection*{Clustered Nested Ellipses (NEC), see \Cref{fig:ellipsedata} (b)}
Let $G^c_1,\dots ,G^c_J\sim \text{Poi}(\lambda_c)$ for $c=1,\dots ,5$. Let $\lambda_3=2$ and set $\lambda_c=1$ else. Let $K_i=M\sum_{c=1}^{5}G^c_i$. Set $w_k^i$ equal to $1$ for $1\leq k\leq K_i$ for $i=1,\dots,J$. Let $t_1,\dots ,t_M$ be a discretisation of $[0,2\pi]$. Let $U^i_1,\dots,U^i_K,V^i_1,\dots,V^i_K\sim U[0.2,1]$. Let $\alpha=(2,12,12,22,12)^T$ and $\beta=(12,2,12,12,22)^T$. Set for $c=1,\dots,5$ and $j_i=0,\dots ,G_i^c$
\[
l^i_{M\left(\sum_{r=1}^{c-1}G^r_i\right)+M(j_i-1)+k}=\frac{1}{24}((3^{-j}U_{M(j_i-1)+k}\sin(t_k)+\alpha_c,3^{-j}V_{M(j_i-1)+k}\cos(t_k)+\beta_c))^T,
\]
where we use the convention that a sum is zero if its last index is smaller than its first one.
\begin{figure}[H]
  \centering
  \subfloat[][]{\includegraphics[width=0.45\linewidth]{graphics/ne_data.png}}%
  \qquad
  \subfloat[][]{\includegraphics[width=0.45\linewidth]{graphics/NEC_Data.png}}%
  \caption{\textbf{(a)} An example of $J=8$ measures from the \emph{NE} dataset with $M=200$. \textbf{(b)} An example of $J=8$ measures from the \emph{NEC} dataset with $M=95$.}%
  \label{fig:ellipsedata}
\end{figure}
\subsubsection*{Spirals of varying Length (SPI), \Cref{fig:spiraldata} (a)}
Let $a_i\sim U[2,4]$ and $b_i\sim U[3,6]$ for $i=1,\dots,J$. Let $K_i=\lceil b_iM \rceil$ and let $t_1,\dots,t_K$ be a discretisation of $[0,b\pi]$. Set $w^i_k=1$ for $k=1,\dots,K_i$ and $i=1,\dots,J$. Set
\[
l^i_k=a_i((t_k\sin(t_k)+64)/140,(t_k\cos(t_k)+70)/130)^T.
\]
\subsubsection*{Clustered Spirals (SPIC), see \Cref{fig:spiraldata} (b)}
Let $a^c_i\sim U[2,4]$ and $b^c_i\sim U[3,6]$ for $i=1,\dots,J$, $c=1,\dots, 5$. Let $K^c_i=\lceil b_i^cM \rceil$ and let $t_1,\dots,t_K$ be a discretisation of $[0,b\pi]$. Set $w^i_k=1$ for $k=1,\dots,K_i$ and $i=1,\dots,J$ and let $\alpha=(0,3,3,6,3)^T$ and $\beta=(3,0,3,3,6)^T$. Set
\[
l^i_{\sum_{r=1}^{c-1}K^r_i+k}=(1/7)(((a_i^ct_k\sin(t_k)+64)/140)+\alpha_c,(a_i^ct_k\cos(t_k)+70)/130+\beta_c)^T,
\]
where we again use the convention that a sum is zero if its last index is smaller than its first one.
\begin{figure}[H]
  \centering
  \subfloat[][]{\includegraphics[width=0.45\linewidth]{graphics/SPI_data.png}}%
  \qquad
  \subfloat[][]{\includegraphics[width=0.45\linewidth]{graphics/SPIC_data.png}}%
  \caption{\textbf{(a)} An example of $J=8$ measures from the \emph{SPI} dataset with $M=110$. \textbf{(b)} An example of $J=8$ measures from the \emph{SPIC} dataset with $M=22$.}%
  \label{fig:spiraldata}
\end{figure}
\subsection{Simulation Results for the $\mathbf{(2,C)}$-Kantorovich-Rubinstein Distance}
\begin{figure}
\begin{tabular}{cc}
      \includegraphics[width=0.45\textwidth]{graphics/comPoi_NE_C0.01.png} &       \includegraphics[width=0.45\textwidth]{graphics/comPoi_NE_C0.1.png} \\
    \includegraphics[width=0.45\textwidth]{graphics/comPoi_NE_C1.png} &       \includegraphics[width=0.45\textwidth]{graphics/comPoi_NE_C10.png} \\
\end{tabular}
\caption{Expected relative $(2,C)$-KRD error for two measures in the Poisson sampling model for the NE class and different success probabilities $s$. For each pair of success probability $s$ and observation time $t$ the expectation is estimated from $1000$ independent runs. Set $M=100$. From top-left to bottom-right we have $C=0.01,0.1,1,10$, respectively.}
  \label{fig:dist_poiNE}
\end{figure}

\begin{figure}
\begin{tabular}{cc}
      \includegraphics[width=0.45\textwidth]{graphics/comPoi_NEC_C0.01.png} &       \includegraphics[width=0.45\textwidth]{graphics/comPoi_NEC_C0.1.png} \\
    \includegraphics[width=0.45\textwidth]{graphics/comPoi_NEC_C1.png} &       \includegraphics[width=0.45\textwidth]{graphics/comPoi_NEC_C10.png} \\
\end{tabular}
\caption{As in \Cref{fig:dist_poiNE}, but for the NEC class and $M=75$.}
  \label{fig:dist_poiNEC}
\end{figure}
In the following, we discuss the results from our simulation studies for the Poisson model for the $(2,C)$-KRD between two measures within one of the eight classes of measures introduced above. Note, that for brevity the plots for some of the classes have been omitted to \Cref{sec:addfigpoi} and exemplary the classes PI, NE and NEC are discussed. The remaining plots in the Figures \ref{fig:dist_poiPIG}, \ref{fig:dist_poiNI}, \ref{fig:dist_poiNIG}, \ref{fig:dist_poiSPI} and \ref{fig:dist_poiSPIC} are found in \Cref{sec:addfigpoi}.\\
For the error of the NE class in \Cref{fig:dist_poiNE} it can be seen that the error is decreasing in $s$ and $t$, but increasing in $C$. Both of these behaviours are in line with the bound in \Cref{thm:samplingboundKR}. The decrease of the error for increasing $s$ and $t$ is immediately clear from our theoretical results. The increase of error for increasing $C$ is based on the fact that in the Poisson model the population total intensities of $\mu$ and $\nu$ are unknown and have to be estimated from the data. The $(p,C)$-KRD penalises mass deviation with a factor scaling with $C$, so naturally for increasing $C$, the errors in the estimation of the true difference of masses yields an increase in the expected relative $(p,C)$-KRD error.
\begin{figure}
\begin{tabular}{cc}
      \includegraphics[width=0.45\textwidth]{graphics/comPoi_PI_C0.01.png} &       \includegraphics[width=0.45\textwidth]{graphics/comPoi_PI_C0.1.png} \\
    \includegraphics[width=0.45\textwidth]{graphics/comPoi_PI_C1.png} &       \includegraphics[width=0.45\textwidth]{graphics/comPoi_PI_C10.png} \\
\end{tabular}
\caption{As in \Cref{fig:dist_poiNE}, but for the PI class and $M=450$.}
  \label{fig:dist_poiPI}
\end{figure}
Notably, while the decrease in $s$ and $t$ is similar for the error of the NEC class in \Cref{fig:dist_poiNEC}, the error is no longer increasing in $C$. Instead the errors increase from $C=0.01$ to $C=0.1$, but then decrease going to $C=1$. Afterwards they increase again at $C=10$. This difference in behaviour is explained by the cluster structure of the measure in NEC. There is still the general trend of increasing error for increasing $C$, as present in the NE class, but now there is an additional change in behaviour based on the fact if transport occurs within clusters or between clusters. From $C=0.1$ to $C=1$, we pass the size of the clusters and the distance between the clusters. Thus, $C=1$ is the first value in our simulation for which inter-cluster transports can occur. This causes a decrease in error, as the impact of the estimation of the total mass intensity of a measure within one cluster is decreased. After this point the usual increase in error for increasing $C$ due to the estimation of the total mass intensity occurs again. For the $(p,C)$-KRD error in the PI class in \Cref{fig:dist_poiPI} this effect is particularly strong. The error increases on average  about two orders of magnitude from $C=0.01$ to $C=10$. This is explained by the fact that the total mass intensity in this class is significantly larger than for the classes NE and NEC, where each location in the support of the measures has mass one. This also causes an increase of the variance of the mass of the empirical measures at each location, which causes a faster increase of error for increasing $C$ at all scales of $C$.

\subsection{Simulation Results for the $\mathbf{(2,C)}$-Barycenter}

\begin{figure}
\begin{tabular}{ccc}
  \includegraphics[width=0.31\textwidth]{graphics/comPoi_NE_bary_C0.1.png} &
    \includegraphics[width=0.31\textwidth]{graphics/comPoi_NE_bary_C1.png} &       \includegraphics[width=0.31\textwidth]{graphics/comPoi_NE_bary_C10.png}
\end{tabular}
\caption{Expected relative Fr\'echet error for the $(2,C)$-barycenter for $J=5$ measures from the NE class in the Poisson sampling model with different success probabilities $s$. For each pair of success probability $s$ and observation time $t$ the expectation is estimated from $100$ independent runs. Set $M=100$. From left to right we have $C=0.1,1,10$, respectively.}
  \label{fig:bary_poiNE}
\end{figure}

It remains to discuss the results from our simulation studies for the Poisson model for the $(2,C)$-barycenter between sets of measures within one of the eight classes of measures introduced above. We again exemplary discuss the classes NE, NEC and PI. The remaining plots in the Figures \ref{fig:bary_poiPIG}, \ref{fig:bary_poiNI}, \ref{fig:bary_poiNIG}, \ref{fig:bary_poiSPI} and \ref{fig:bary_poiSPIC} are found in \Cref{sec:addfigpoi}. We also restrict our analysis to the values of $C=0.1,1,10$, since for $C=0.01$ the $(p,C)$-KRD is close to the TV distance. In particular, for all classes except PIG and NIG, where all measures share the same support grid, the $(2,C)$-barycenter will be close or identical to the zero measure, since the measures in the other classes are almost surely disjoint. Additionally, if the barycenter of the population measures is the zero measure, any empirical barycenter has mass zero as well. Thus, there is little merit in simulating the barycenters in these cases. For the classes NIG and PIG the barycenters are essentially TV-barycenters for small $C$ which removes any geometrically interesting features from the barycenter. Finally, for extremely small values of $C$ the $(p,C)$-barycenter computations tend to become numerically unstable due to either involving values close to machine accuracy and UOT plans for this values of $C$ often being close to the zero measure. Hence, empirical simulations of the expected relative Fr\'echet error would also be less reliable in this regime of values for $C$. In summary, empirical analysis of the properties of the $(p,C)$-barycenter for values of $C$ which are several orders of magnitude smaller than the diameter of $\mathcal{Y}$ is inadvisable. \\
\begin{figure}
\begin{tabular}{ccc}
    \includegraphics[width=0.31\textwidth]{graphics/comPoi_NEC_bary_C0.1.png} &
    \includegraphics[width=0.31\textwidth]{graphics/comPoi_NEC_bary_C1.png} &       \includegraphics[width=0.31\textwidth]{graphics/comPoi_NEC_bary_C10.png}
\end{tabular}
\caption{As in \Cref{fig:bary_poiNE}, but for the NEC class and $M=75$.}
  \label{fig:bary_poiNEC}
\end{figure}
For the relative Fr\'echet errors we observe significant changes in behaviour compared to the relative $(p,C)$-KRD before. Considering the error for the NE class in \Cref{fig:bary_poiNE}, we note that for $C=0.1$ the behaviour in $s$ and $t$ is different than for the $(p,C)$-KRD. Namely, for fixed $s$, the error is in general not strictly decreasing in $t$ and vice versa for fixed $t$, the error is not always strictly decreasing in $s$. This is an interesting effect arising for small values of the product $st$. A point $y\in \mathcal{Y}$ can only be a support point of a $(p,C)$-barycenter if it is in the intersection of at least $J/2$ balls of size $C^p/2$ around support points of different measures (compare the construction of the centroid set in \eqref{eq:KRCentroid}). Now, for small $s$ and $t$ many support points of the population barycenter are not included in the support of the empirical one, since centroid set of the empirical measures is significantly smaller than the population level one. In particular, this can create situations where an increase in $s$ or $t$ on average adds support points to empirical barycenter, which cause the relative error to increase, since placing mass zero at this location, for small $C$, is actually better than placing a potentially larger mass (since we assumed $(ts)^{-1}$ to be relatively small) at this location. Thus, while asymptotically, the rate in \Cref{thm:frechetboundPoi} is optimal, for certain, sufficiently small, values of $s$, $t$ and $C$, the behaviour of the relative Fr\'echet error might be counter-intuitive. For $C=0.1$ and $C=1$, the errors behave quite similarly to the $(p,C)$-KRD setting, though there is essentially no increase in error going from $C=1$ to $C=10$. This is explained by two points. First, the location of the $(p,C)$-barycenter tends to be more centred within the support of the measures (all measures are support on subsets of the unit square), so little transport between the barycenter and the $\mu^i$ occurs at a distance larger than one. Second, the key factor for the increasing error for increasing $C$ in the $(p,C)$-KRD case is the estimation error for the total mass intensities. However, for sufficiently large $C$ the mass of the $(p,C)$-barycenter is the median of the total masses of the $\mu^i$. Since, this quantity is significantly more stable under estimation than the individual total mass intensities, it is to be expected that the mass estimation has little effect on the relative Fr\'echet error. For the error of the NEC class in \Cref{fig:bary_poiNEC} the results are similar to the NE case. We observe similar effects on the dependence of $s$ and $t$ for $C=0.1$ and for $C=1$ and $C=10$, the errors look extremely similar. One notable distinction is the fact the from $C=0.1$ to $C=1$ the errors decrease on average. As before for the NEC class in the $(p,C)$-KRD setting, this can be explained by its cluster structure and $C=1$ being the first value for which inter-cluster transport becomes possible in an UOT plan. This is therefore also the first value of $C$ which allows the $(p,C)$-barycenter to have mass between clusters. Finally, for the error of the PI class in \Cref{fig:bary_poiPI} the value of $C$ only has a minimal effect on the resulting errors. Notably and contrary to the two prior classes, we do not encounter any additional effects for $C=0.1$. This is explained by the in general higher mass intensities of the measures in the PI class, which make the previously described effects due to low values of $s$ and $t$ less likely. Additionally, these measures do not possess any geometrical structures in their support, which could impact the behaviour on different scales. There is again little increase in error for increasing $C$, which is in stark contrast to the PI class in the $(p,C)$-KRD setting, where the error increased by multiple orders of magnitude. This is another strong indicator, that the Fr\'echet error is significantly more stable under $C$, due to the stability of total mass intensity of the empirical barycenter opposed to the total mass intensity of the individual measures.
\begin{figure}
\begin{tabular}{ccc}
     \includegraphics[width=0.31\textwidth]{graphics/comPoi_PI_bary_C0.1.png} &
    \includegraphics[width=0.31\textwidth]{graphics/comPoi_PI_bary_C1.png} &       \includegraphics[width=0.31\textwidth]{graphics/comPoi_PI_bary_C10.png}
\end{tabular}
\caption{As in \Cref{fig:bary_poiNE}, but for the NEC class and $M=450$.}
  \label{fig:bary_poiPI}
\end{figure}
\subsection{Real Data Example}\label{sec:realdata}
In \Cref{fig:realdata} we consider the $(2,0.1)$-KRD between images which are an excerpt from STED microscopy of adult human dermal fibroblast cells (for the full dataset see \cite{tameling2021colocalization}). The images in the Figures \ref{fig:realdata}(a),(b) and Figures \ref{fig:realdata}(c),(d) are visually similar, as they correspond to measurements taken based on two different markers (one at the inner mitochondrial membrane and one at the outer) in the same cells. The $(2,0.1)$-KRD captures this fact in the sense, that the pairwise distance between the measures are smallest for these pairs of images. Utilising UOT on this type of datasets is a potential way of quantifying dissimilarity between the respective measures and extending OT based dissimilarity analysis to measures of unequal total intensity (the total mass intensities in this examples lies roughly between $6200$ and $9500$. \\
We further want to use this dataset to illustrate the performance of the randomised computational approach for the $(p,C)$-KRD based on the multinomial model (recall \eqref{eq:multinomialmeasure}). The $300\times 300$ images here are specifically chosen such that the true distances can still be computed which allows to compare the expected error of the empirical $(p,C)$-KRD for given sample sizes on this data set. We compare the results obtained from the resampling approach (i.e. the estimator from \eqref{eq:multinomialmeasure}) considered in the multinomial model to the subsampling approach (i.e. the estimator from \eqref{eq:subsamp}) obtained by sampling without replacement from the measures instead. In these simulations the maximum sample size is about $1/5$ of the support sizes. This corresponds to a runtime of about $2.5\%$ of the original problem size. While it is clear by construction that for sufficiently large sample sizes, subsampling yields a smaller error than the resampling (as the error approaches zero if the sample size approaches the support size), for smaller sample sizes the resampling can have a better performance. It yields a relative error below $5\%$ at less than $10\%$ of the original support size in all considered instances. This approximation can be achieved in around $0.5\%$ of the original runtime. The subsampling approach does not reach this level of accuracy for the considered sample sizes. Thus, these simulations suggest that randomised computations based on the multinomial model allow for high accuracy approximations of the $(p,C)$-KRD in real data applications at a significantly lower computational cost than the original problem and that for small sample sizes there are scenarios where the resampling approach yields significantly better performance than the subsampling one.
\begin{figure}
  \centering
  \subfloat[][]{\includegraphics[width=0.25\linewidth]{graphics/IM2.png}}
    \subfloat[][]{\includegraphics[width=0.25\linewidth]{graphics/IM1.png}}
      \subfloat[][]{\includegraphics[width=0.25\linewidth]{graphics/IM5.png}}
  \subfloat[][]{\includegraphics[width=0.25\linewidth]{graphics/IM6.png}} \hspace{0.01\linewidth}
      \subfloat[][]{\includegraphics[width=0.33\linewidth]{graphics/image_uot_samples5000_C01_1_2resampvssubsamp.png}}
    \subfloat[][]{\includegraphics[width=0.33\linewidth]{graphics/image_uot_samples5000_C01_1_3resampvssubsamp.png}}
      \subfloat[][]{\includegraphics[width=0.33\linewidth]{graphics/image_uot_samples5000_C01_1_4resampvssubsamp.png}} \hspace{0.01\linewidth}
        \subfloat[][]{\includegraphics[width=0.33\linewidth]{graphics/image_uot_samples5000_C01_2_3resampvssubsamp.png}}
    \subfloat[][]{\includegraphics[width=0.33\linewidth]{graphics/image_uot_samples5000_C01_2_4resampvssubsamp.png}}
      \subfloat[][]{\includegraphics[width=0.33\linewidth]{graphics/image_uot_samples5000_C01_3_4resampvssubsamp.png}}
  \caption{\textbf{(a)-(d):} Excerpts of size $300\times 300$ from the STED microscopy data of adult Human Dermal Fibroblasts in \cite{tameling2021colocalization}. The images have on average about $25000$ non-zero pixels. (a) and (c) have been labelled at MIC60 (a mitochondrial inner membrane complex); (b) and (d) have been labelled at TOM20 (translocase of the outer mitochondrial membrane). \textbf{(e)-(j):} The relative error for the empirical $(2,0.1)$-KRD (obtained from resampling and subsampling) between the four filament structures in (a)-(d) considered as measures in $[0,1]^2$. \textbf{(e)} Between (a) and (b). \textbf{(f)} Between (a) and (c). \textbf{(g)} Between (a) and (d). \textbf{(h)} Between (b) and (c). \textbf{(i)} Between (b) and (d). \textbf{(j)} Between (c) and (d).}%
  \label{fig:realdata}
\end{figure}
\section*{Acknowledgments}
F.\  Heinemann and M.\ Klatt gratefully acknowledge support from the DFG Research Training Group 2088 \textit{Discovering structure in complex data: Statistics meets optimization and inverse problems}. A.\ Munk gratefully acknowledges support from the DFG CRC 1456 \textit{Mathematics of the Experiment A04, C06} and the DFG Cluster of Excellence 2067 MBExC \textit{Multiscale bioimaging--from molecular machines to networks of excitable cells}.

\bibliographystyle{abbrvnat}

redpajama_set_name: RedPajamaArXiv


Your Professional Option for Search Engine Optimization and Marketing | Website Design Services, Maintenance and Hosting Solutions in Livonia, Michigan!
Xtreme Internet Marketing™ – Creating Websites That Get You Notice! Providing Livonia, Michigan Clients Custom Responsive Websites, Local and National Search Engine Optimization Services.
Xtreme Internet Marketing™ works with clients from across Livonia, Michigan, helping with; new site designs, redesigns, updates, modifications, SEO and web hosting services. Our talented team of dedicated professionals will execute together your vision and ideas into reality with creativity, technology, and marketing expertise. We welcome you to experience for yourself our extraordinary customer service, and creative, affordable website designs.
Livonia, Michigan Website Design, Search Engine Optimization and Marketing Services.
We specialize in web design and online marketing. Working with clients across Livonia, Michigan, big and small, creating powerful and effective websites as well as online marketing campaigns. All sites built by Xtreme Internet Marketing™ are custom, fully responsive, Google compliant, SEO-Optimized and include a FREE 1-year SSL certificate to increase search rankings and secure customer data.
Our, "Listed First Local and National Organic Listings", are the Most Popular SEO Service – Offering Our Livonia, Michigan Clients the Best Prices in the Industry with Proven Results!
Develop Your Company Proficiently with the Best Choice (SEO) Search Engine Optimization Services Company in Livonia, Michigan.
Trust the certified experts with your online marketing strategy techniques at an affordable price that gets your web site found on Google, Yahoo!, Bing, and YouTube with local search marketing from Xtreme Internet Marketing™ serving Livonia, Michigan.
Xtreme Internet Marketing™ is Your Affordable Choice for Search Engine Optimization and Marketing Services throughout Livonia, Michigan!
Local and Nationwide – Performance, Results, Rankings, Guaranteed! We can easily say that ninty five percent of our clients are on the First Page of Google and Yahoo Local! Can your existing SEO guy say that?
This is not PPC! This is a month to month paid sevice. NO CONTRACTS, you can stop anytime. We can start your service for as little as $200 per month! Call us to discover more about this amazing service – With Undoubted Results Guaranteed!
We get your website on the First Page of Google and Yahoo Local and Nationwide with the best price in the industry – With Undoubted Results for those across the Livonia, Michigan!
Having your website ranking #1 on search engines brings you more traffic from new customers searching for your services or products. At Xtreme Internet Marketing™, with our (SEO) search engine optimization and Marketing services, it is our mission to increase your calls, leads and sales through organic search engine optimization. Our team of (SEO) search engine optimization experts have the right combination of knowledge, experience and resources to attain this goal for you, with verified results.
It's Why We've Been Named The Top Choice for Search Engine Optimization (SEO) in Livonia, Michigan.
Our proprietary approach to the optimization process is highly customized and tailored to your needs, from beginning to end. There is no one solution that works for everyone. Knowing what strategies to execute at what times, and how to execute them for every business is critical. Apart from that, careful data collection is the key to understanding how to enhance performance.
(SEO) search engine optimization is the process, of improving the volume or best quality traffic to a web site from search engines via natural (organic) or un-paid search results as opposed to search engine marketing which deals with paid inclusion. Most always, the higher a site appears in the search results list, the more visitors it will receive from the search engine.
With that said, it refers to the procedure that permits your web page rank higher in organic search results, making your web page more visible to people who are looking for your brand, product, or service via search engines like Google, Yahoo, and Bing.
As part of the Xtreme Internet Marketing™ (SEO) search engine optimization agreement, you submit a one-time initial setup fee for 30 days of intensive optimization – after which your front month begins at the monthly rate determined by our listing certified experts. This fee covers your business' category in one city. Prices for other cities and categories vary according to factors; such as size of city, number of competitors, number of locations, size of business. Volume discounts are available.
Allowing for optimization of Google's maximum 150 characters, five keyword phrases (that a lot of keywords), to millions of likely consumers within a record setting 24 hrs.. The software forms 28,000 plus, United States city landing pages, optimized for up to five keyword phrases, within one day.
Our proprietary software moderates everything, allowing us to be the topmost local organic (SEO) search engine optimization certified experts in the world. We can now explicitly target any city or state in the United States. Search engines rapidly identify and rank the listings due to our dynamically created local (SEO) search engine optimization landing pages. There is nothing that works better or faster. Nothing!
Our unique (SEO) search engine optimization local landing pages are also optimized for voice search from smart phones, now becoming an industry standard. Local keywords are now and will be the trend for small business online marketing for many years to come. Xtreme Internet Marketing™'s "Listed First" software empowers us to tackle the (SEO) search engine optimization market with targeted and absurdly fast results. It permits us to optimize your listings for 100's to 1000's of keyword phrases. We will maximize your marketing dollar by optimizing your pages for the service or delivery areas your work in or targeting, while still maintaining great profit margins.
As part of the Xtreme Internet Marketing™ (SEO) search engine optimization agreement, you submit a one-time initial setup fee for 30 days of intensive optimization — after which your first month begins at the monthly rate determined by our listing certified experts. This fee covers your business' category in one city. Prices for other cities and categories vary according to factors; such as size of city, number of competitors, number of locations, size of business. Volume discounts are available.
Even if a process was done perfectly in a given time and this will bring the site in front position, nobody and nothing can guarantee that it will remain there forever. Increasing rivals for highest rank spots periodically veer, changes in the indexing algorithms are a non-stop threat to your position in the top rankings.
SEO and Marketing Services in Livonia, Michigan.
Xtreme Internet Marketing™ of Livonia, Michigan has presented (SEO) search engine optimization for highest rank placement; local and national search, organic listings, PPC, social media marketing, professional website design services, maintenance and hosting services for 1000s of folks in Livonia, Michigan and all across the country! Xtreme Internet Marketing™ takes genuine pleasure in helping in your businesses success. Making sure you're getting expert services within (SEO) search engine optimization for highest rank placement in your local and national search, organic listings, PPC, Social Media Marketing (SSM), professional website design services, maintenance and hosting solutions! Our (SEO) search engine optimization, website designers and hosting solutions specialists work their hardest to make sure your website is listed on the first page every day!
Xtreme Internet Marketing™ certainly is the premiere option for (SEO) search engine optimization for highest rank placement; local and national search, organic listings, PPC, social media marketing (SSM), professional website design services, maintenance and hosting solutions in Livonia, Michigan.
Xtreme Internet Marketing™ offers the best possible search engine optimization (SEO), maintenance and hosting solutions when your website stops or your website no longer comes up on the first page when searching on Google, Yahoo!, Bing and MSN in Livonia, Michigan. Permit our SEO, website designers and hosting solutions certified professionals of Livonia, Michigan get your business back online with our first class service. It is critical to us the success of your business by making sure your website, local and national search campaigns are working and back on line right now, continuing to making you money.
Xtreme Internet Marketing™ of Livonia, Michigan has expert and dedicated (SEO) search engine optimization for highest rank placement; local and national search, organic listings, PPC, Social Media Marketing (SSM) consultants waiting to aid you with your website design, (SEO) search engine optimization marketing services or hosting needs!
Social Media Marketing (SSM) is more important than ever. Using Xtreme Internet Marketing™ Social Media Marketing (SSM) professionals, we can help you create a social campaign that will increase traffic, promote your brand and increase your presence in Livonia, Michigan and saturating the marketplace using the most widely used networking platforms such as Pinterest, Google+, LinkedIn, Facebook, and Twitter provide a great way to stay involved with clients, build a community of users, increase brand awareness and upsell your products or services. Social Media Marketing (SSM) is also a great way to offer special deals and promotions in real time to invite new customers and more traffic to your website.
Since no two companies are alike, we provide customized social media packages for each of our clients based on their individual needs. Call us and will go over your business goals and current social media ambitions. Next, we produce a strategy to get your social media endeavours off the ground. We can also leverage your existing social media campaigns and take them to the next level.