TITLE: Homomorphism Related Questions
QUESTION [1 upvotes]: Let G and G' be finite groups and $\phi$: G $\rightarrow$ G' a homomorphism.
(a) Suppose that |G'| is prime. Show that $\phi$ must either be trivial or surjective.
(b) Suppose that |G| is prime. Show that $\phi$ must either be trivial or injective.
(c) Suppose that |G| and |G'| are the same prime number. Show that $\phi$ is either trivial or an isomorphism.
(d) Suppose that |G| and |G'| are relatively prime. Show that $\phi$ must be trivial.
My approach for (a) is: Since |G'| is prime. The the only two subgroups of G' must be the id. element and G' itself. Let H be a subgroup of G, then $\phi (H)$ must be a subgroup in G'. Then $\phi (H)$ must be either the id. element or the G'. Therefore, it's trivial or surjective. Is this correct??
for (b), could I do it with the same approach: |G| is prime. Then the subgroups are id. element and G itself. Define H to be the subgroup of G. But I do not how to go from here.
for (c). Both |G| and |G'| are prime. Then the subgroups of G are id. element and G. For G', its {e'} and G'.  If H (H is subgroup of G) is {e}, then $\phi$ is trivial. If H is G, then it maps G to G' which is one to one and onto. I wonder how can I prove its one-to-one and onto here.
for (d). I really don't know what to start with. Should go with the same approach?

REPLY [0 votes]: a. You're on the right track. $\phi(G)$ is a subgroup of $G^{\prime}$, hence is either trivial or all of $G^{\prime}$.
b. Try the same argument as (a), but now applied to the subgroup $\mathrm{ker}(G)$ of $G$.
c. Combine (a) and (b).
d. $\phi(G)$ is a subgroup of $G^{\prime}$, so its order divides $|G^{\prime}|$. But it's also a quotient of $G$, so its order divides $|G|$. Since $|G|$ and $|G^{\prime}|$ are coprime, this means that $\phi(G)$ is trivial.


TITLE: Conversion of the units BeV (US) and GeV (UN)
QUESTION [1 upvotes]: What is the difference in the definition of a billion electron volts in United states (US) and United Nations (UN)?
When the US people say billion, do they mean $10^{12}$ or $10^9$?

REPLY [1 votes]: In the US, 1 billion = $10^9$.
The difference is between the Long and short scales. The US uses the short scale, where a billion is $10^9$. In the long scale, a billion is $10^{12}$.
In the short scale, every term after a million (billion, trillion, etc.) is 1,000 times bigger than the previous one. So, million = $10^6$, billion = $10^9$, trillion = $10^{12}$.
In the long scale, every term after a million is 1,000,000 times bigger. So, million is again $10^6$, but billion = $10^{12}$, and trillion = $10^{18}$.

If you just use SI prefixes, though, those are the same everywhere. $10^9 eV$ is one GeV, whereas $10^{12} eV$ is one TeV.


text:
\begin{document}

\title{Stochastic Saddle-Point Optimization for the Wasserstein Barycenter Problem

\thanks{
The work of A. Gasnikov in Section \ref{sec:saddle_point_repr} was partially supported by the Ministry of Science and Higher Education of the Russian Federation (Goszadaniye) 075-00337-20-03, project no. 0714-2020-0005.
The work of P. Dvurechensly in Section \ref{sec:algorithms} was funded by Russian Science Foundation (project 18-71-10108).

}}

\author{Daniil Tiapkin \\
 HSE University, Moscow, Russia \\
              \texttt{dtyapkin@hse.ru}
\And
        Alexander Gasnikov\\
        Moscow Institute of Physics and Technology, Moscow, Russia \\
              Institute for Information Transmission Problems, Moscow, Russia \\
              HSE University, Moscow, Russia\\
             Weierstrass Institute for Applied Analysis and Stochastics, Berlin, Germany\\
              \texttt{gasnikov@yandex.ru}
              \And
        Pavel Dvurechensky\\
                Weierstrass Institute for Applied Analysis and Stochastics, Berlin, Germany \\
                Institute for Information Transmission Problems, Moscow, Russia \\
                \texttt{pavel.dvurechensky@wias-berlin.de}
}

\maketitle

\begin{abstract}
We consider population Wasserstein barycenter problem for random
probability measures supported on a finite set of points and generated by an online stream of data. This leads to a complicated stochastic optimization problem where the objective is given as an expectation of a function given as a solution to random optimization problem. We employ the structure of the problem and obtain a convex-concave stochastic saddle-point reformulation of this problem.
In the setting when the distribution of random
probability measures is discrete, we propose a stochastic optimization algorithm and estimate its complexity.
The second result, based on kernel methods, extends the previous one to the arbitrary distribution of random
probability measures. Moreover, this new algorithm has a total complexity better than the Stochastic Approximation approach combined with the Sinkhorn algorithm in many cases. We also illustrate our developments by a series of numerical experiments.

\keywords{Stochastic Optimization \and Saddle-Point Optimization \and Computational Optimal Transport \and Wasserstein Barycenter Problem}

\end{abstract}

\section{Introduction}
\label{intro}
In this paper we consider stochastic optimization problems, which arise in computational optimal transport when the goal is to estimate a population Wasserstein barycenter \cite{agueh2011barycenters} (or Fr\'echet mean w.r.t. Wasserstein distance) of a probability distribution on the Wasserstein space of probability measures. Wasserstein barycenter problem has recently attracted a lot of attention from the Machine Learning, statistics and optimization community, see \cite{peyre2018optimaltransport,dvinskikh2020stochastic} and references therein. From the computational point of view, approximating population Wasserstein barycenter is a challenging optimization problem since it contains several layers of complications. The first layer requires to define for $p\geq 1$ the $p$-Wasserstein distance $\W_p(r,c)$ between two probability measures $r,c$, which in this paper we assume to be discrete probability distributions which belong to the standard simplex $\Delta^n$. $\W_p(r,c)$ is defined as an optimal value in an optimal transport problem, where the goal is to find a transport plan between vectors $r$,$c$ such that the total transportation cost is minimal (see \eqref{eq:ot_problem} for a formal definition). The next layer is that the measure $c$ is assumed to be random with some probability distribution $P_c$ on $\Delta^n$ and the population barycenter $r^*$ minimizes the expected Wasserstein distance, i.e. $r^*=\argmin_{r \in \Delta^n} \E_c \W_p^p(r, c)$.
This is a stochastic optimization problem, where the objective is defined as a solution to some other random optimization problem. The third layer of complication is that the dimension $n$ can be very large, e.g. $10^6$. To motivate this setting we refer to \cite{boissard2015distribution}, where the authors consider a template estimation problem from a random sample of its transformations. A particular example of such template can be images, considered in \cite{cuturi2014fast}, where empirically the Wasserstein barycenter gives the best quality of the template image reconstruction. For a $10^3$ by $10^3$ image we obtain that $n=10^6$.

Similarly to other stochastic optimization problems, see e.g. \cite{nemirovski2009robust}, there are two approaches to the population Wasserstein barycenter problem, which we refer to as offline (also known as Sample Average Approximation, SAA) and online (also known as Stochastic Approximation, SA). The offline approach assumes that a number of random probability measures $c_i$, $i=1,...,m$ is sampled in advance and all the measures are stored in memory. Then the population Wasserstein barycenter problem is approximated by the empirical Wasserstein barycenter problem $\hat{r}=\argmin_{r \in \Delta^n} \frac{1}{m} \sum_{i=1}^m\W_p^p(r, c_i)$, i.e. the expectation is approximated by the finite-sum. Funally it is assumed that an optimization algorithm has access to all of the measures $c_i$, $i=1,...,m$ immediately, probably in parallel, and the goal is to minimize this finite-sum.
So far, nearly all algorithms in the literature for approximating Wasserstein barycenter are designed for this finite-sum problem, to name a few  \cite{cuturi2014fast,benamou2015iterative,staib2017parallel,claici2018stochastic,lin2020revisiting,dvinskikh2020improved,heinemann2020randomised}. Moreover, in this setting, it is possible to use a big arsenal of modern decentralized distributed optimization algorithms, see for example, \cite{dvurechensky2018decentralize,dvinskikh2019primal,gorbunov2019optimal,hendrikx2020optimal,gorbunov2020recent}, where such algorithms are developed with best-known complexity bounds for general finite-sum optimization with cheap dual stochastic oracle, which turns out to be available for the Wasserstein distance.

Despite all these developments, the offline approach does not consider clearly the question of how well the solution to the finite-sum problem approximates the target population Wasserstein barycenter which is a solution to stochastic optimization problem. Moreover, in practice it may not be possible to generate all the measures in advance, the storage could be problematic, and the data consisting of random measures can appear as an online stream.
Thus, in this paper, we focus on the less studied online setting with all measures unavailable in advance.
Online estimation techniques have been very recently developed in \cite{mensch2020online} for Wasserstein distance.
Online setting for population Wasserstein barycenter has started recently to attract interest in the statistical community \cite{kroshnin2019statistical,chewi2020gradient} where population Wasserstein barycenters are studied in a different from ours setting of Bures-Wasserstein space of Gaussian probability measures, and in optimization community \cite{dvinskikh2020stochastic}, where online approach is compared to offline.

The approach of the latter paper is based on a widespread idea of entropy regularized optimal transport \cite{cuturi2013sinkhorn,peyre2018optimaltransport}. One of the main negative sides of the entropic regularization is that if the goal is to find the an $\varepsilon$-approximation to the original non-regularized Wasserstein distance or barycenter, the regularization parameter has to be chosen proportional to $\varepsilon$ \cite{dvurechensky2018computational,lin2019efficient,kroshnin2019complexity,lin2020revisiting}. This leads to numerical instability \cite{cuturi2013sinkhorn,stonyakin2019gradient} in the Sinkhorn's algorithm for Wasserstein distance and Sinkhorn-type algorithms for Wasserstein barycenters \cite{benamou2015iterative}.
Alternative methods, e.g. based on accelerated gradient method \cite{dvurechensky2018computational,guminov2019acceleratedAM,lin2019efficient,lin2020revisiting}, similarly to Sinkhorn's algorithm require the regularization parameter to be small.
Moreover, entropic regularization leads to a blurred reconstruction of the original template \cite{cuturi2014fast}.

Thus, in this paper we focus on the following question and give a pisitive answer to it.
{\textit{Is it possible to find a population Wasserstein barycenter of a set of discrete probability measures in the online setting and without the use of entropic regularization and without the calculation of (regularized) Wasserstein distance or its (sub)gradients.}}

To propose such an online approach that does not require to calculate Wasserstein distance and its subgradient, we use a saddle-point representation for the Wasserstein barycenter problem. Since the population Wasserstein barycenter problem is a convex stochastic programming problem, we obtain partially stochastic convex-concave saddle-point representation. Then we adapt Stochastic Mirror Descent \cite{nemirovski2009robust} for this, partially infinite-dimensional, saddle-point problem. Finally, by using a reproducing kernel Hilbert space trick inspired by \cite{genevay2016stochastic}, we demonstrate how to make our approach practical and obtain finite complexity bounds.
We compare the proposed methods with the online approach of \cite{dvinskikh2020stochastic} and obtain a regime in which our complexity bounds are better. We also illustrate numerically the regime in which our methods have better performance than the online algorithm in \cite{dvinskikh2020stochastic}.

\paragraph{Notation.}
We define \(\Mat_{n \times m}(X)\) as a space of all matrices of size \(n \times m\) with entries from the set \(X\).
We denote by \([n]\) an \(n\)-element set \(\{1, 2, \ldots, n\}\), and \(e^v, e^A, \log(v), \log(A), \sqrt{v}, \sqrt{A}\) for \(v \in \mathbb{R}^n, A \in \Mat_{n \times n}(\mathbb{R})\) as a point-wise exponent, logarithm and square root respectively.
Also we define \(\Delta^n\) as \(n\)-dimensional probability simplex \(\Delta^n  = \{ (s_1, \ldots, s_n) \mid \forall i \in [n] : s_i \geq 0, \ \ \sum_{i = 1}^n s_i = 1\}\).
\(\mathbf{1}_n\) is \(n\)-dimensional vector consisting of ones. If the dimension is clear from the context, the subscript is omitted. As \(e_i\) we denote the \(i\)-th coordinate vector. The probability measure induced by a random variable \(\xi\) will be denoted by \(P_\xi\).
For matrix \(M \in \Mat_{n \times m}(X)\) we denote  by \(M_{(i)}\) the \(i\)-th row of this matrix.

\section{Background}
\label{sec:background}

\subsection{Background in optimal transport}
\label{sec:ot_background}
In this section, we recall some basic definitions of the optimal transport, following \cite{peyre2018optimaltransport}. Since we deal only with discrete measures, we consider only the discrete-discrete optimal transport problem.

\paragraph{Optimal transport problem.}
    For a fixed non-negative matrix \(C\) and two discrete probability measures \(r,c \in \Delta^n\)  with \(n\)-element support define \textit{a transportation cost} between measures \(r\) and \(c\) associated with the cost matrix \(C\) as a solution to the following optimization problem
    \begin{equation}\label{eq:ot_problem}
        L_{C}(r,c) = \min_{X \in \mathcal{U}(r,c)} \langle C, X \rangle,
    \end{equation}
    where \(X\) is called \textit{a transport plan} and \(\mathcal{U}(r,c)\) is \textit{a transport polytope}, defined as
    \(
        \mathcal{U}(r,c) = \{ X \in \Mat_{n \times n}(\mathbb{R}_+) \mid X \mathbf{1} = r, \ X^T \mathbf{1} = c \}
    \),
    and $\langle \cdot, \cdot \rangle$ is the Frobenius dot-product of two matrices.

If \(r\) and \(c\) are probability measures onto discrete \(n\)-element metric space \((\mathcal{M}, d)\), we can define as a \textit{\(p\)-Wasserstein distance} between these two measures as the \(p\)-th root of a transportation cost associated with the matrix \(D^p_{i,j} = d(x_i, x_j)^{p}\), where \(x_i\) and \(x_j\) are elements of \(\mathcal{M}\). Formally,
\(
    \W_p(r,c) = \left( L_{D^p}(r,c) \right)^{1/p}.
\)
With a slight abuse of notation, we will refer to \(L_C(r,c)\) for an arbitrary non-negative cost matrix \(C\) as a Wasserstein distance too, and denote it by \(\W(r,c)\) when matrix \(C\) was fixed.

\paragraph{Dual problem.} The linear program in the definition of a transportation cost can be reformulate using so-called Kantorovich duality in two ways, using some reformulation of results from \cite{peyre2018optimaltransport}:
\begin{equation}\label{eq:dual_simple}
    L_{C}(r,c) = \max_{\substack{\lambda, \mu \in \mathbb{R}^n\\-C_{i,j}-\lambda_i-\mu_j \leq 0}} - \langle \lambda, r \rangle - \langle \mu, c \rangle
\end{equation}
and the equivalent one is
\begin{equation}\label{eq:dual_problem}
    L_{C}(r,c) = \max_{\mu \in \mathbb{R}^n} - \langle \lambda^*(\mu, C), r \rangle - \langle \mu, c \rangle,
\end{equation}
where \(\lambda^* \colon \mathbb{R}^n \times \Mat_{n \times n}(\mathbb{R}) \to \mathbb{R}^n\) is defined element-wise:
\begin{equation}\label{eq:lambda_def}
    \lambda_i^*(\mu, C) = \max_{j \in [n]} ( -C_{i,j} - \mu_j).
\end{equation}

\paragraph{Barycenter definition.}
Suppose we have a random variable \(\xi \colon \Omega_\xi \to \Delta^n\) on probability simplex, or, equivalently, on the space of probability measures onto \(\mathcal{M}\). Then we can define a \(p\)-Wasserstein (population) barycenter w.r.t. \(\xi\) as the solution to the following optimization problem:
\[
    r_* = \argmin_{r \in \Delta^n} \E \W_p^p(r, \xi)
\]
However, we are interested in a more general situation that depends not on metric but on an arbitrary non-negative cost matrix \(C\):
\[
    r_* = \argmin_{r \in \Delta^n} \E\left[L_C(r, \xi)\right] = \argmin_{r \in \Delta^n} \E\left[\W(r, \xi)\right].
\]

\subsection{Notation for saddle-point problems}
\label{not_saddle_point}

We refer to \cite{bubeck2014convex,nemirovski2009robust} for more detailed background. In this section we consider only some useful definition.

We are interested in the following convex-concave saddle point problems:
\[
    \min_{x \in \mathcal{X}} \max_{y \in \mathcal{Y}} F(x,y).
\]
The typical way to evaluate the quality of the algorithms that outputs the pair \((\widetilde x,\widetilde y)\) is to bound the so-called duality gap:
\[
    \max_{y \in \mathcal{Y}} F(\widetilde x, y) - \min_{x \in \mathcal{X}} F(x, \widetilde{y}) \leq \varepsilon.
\]

For the stochastic setting, where \((\widetilde x, \widetilde y)\) are random variables, we are the most interested in the confidence region criteria:
\begin{equation}\label{eq:confidence_region_criteria}
    \Prob\left[\max_{y \in \mathcal{Y}} F(\widetilde x, y) - \min_{x \in \mathcal{X}} F(x, \widetilde{y}) \geq \varepsilon\right] \leq \sigma.
\end{equation}
However, there is another criteria that was studied:
\begin{equation}\label{eq:expectation_criteria}
    \E\left[ \max_{y \in \mathcal{Y}} F(\widetilde x, y) - \min_{x \in \mathcal{X}} F(x, \widetilde{y}) \right] \leq \varepsilon.
\end{equation}

\section{Saddle-point representation}
\label{sec:saddle_point_repr}
Using a dual approach and the definition of a Wasserstein barycenter, we have the following problem:
\[
    \min_{r \in \Delta^n} \E \W(r,\xi) = \min_{r \in \Delta^n} \E \left[ \max_{\mu \in \mathbb{R}^n} - \langle \lambda^*(\mu, C), r \rangle - \langle \mu, \xi \rangle \right]
\]

For this problem, we will apply Theorem 14.60 from \cite{rockafellar2009va} to the space of all \(P_\xi\)-measurable functions \(\mathcal{F}\). Clearly, it is decomposable. Also clear that the function under maximum is a normal integrand and finite, since the barycenter is defined. Thus, we have the next equality:
\begin{equation}\label{eq:barycenter_basic}
    \min_{r \in \Delta^n} \E \mathcal{W}(r, \xi) = \min_{r \in \Delta^n} \sup_{f_\mu \in \mathcal{F}} \E\left[-\langle \lambda^*(f_\mu(\xi), C), r \rangle - \langle f_\mu(\xi), \xi \rangle \right].
\end{equation}
The Wasserstein barycenter is a solution to this saddle-point problem.

\subsection{Dual variable bound}
\label{sec:dual_var_bound}

For the theoretical guarantees of the future algorithms, we need to construct a bound on the optimal value of the dual variable in (\ref{eq:dual_problem}). It will be obtained by using the properties of an entropy regularization of the optimal transport problem and some additional assumptions on \(r_*\).

The core idea is to use known bound from \cite{dvurechensky2018computational} on the optimal variable in the regularized case and transfer it into non-regularized case.

\paragraph{Entropy regularization.} Firstly, let us define an entropy-regularized OT problem:
\[
    L_{C}^\gamma(r,c) = \min_{X \in \mathcal{U}(r,c)} \langle C, X \rangle + \gamma H(X),
\]
where \(H(X) = \sum_{i,j=1}^n X_{i,j} \log X_{i,j}\) is an entropy regularizer. After the explicit expanding the Lagrange multipliers, we have the equivalent problem:
\begin{equation}\label{eq:dual_simple_regularized}
    L_{C}^\gamma(r,c) = \max_{\lambda, \mu \in \mathbb{R}^n}  -\langle \lambda, r \rangle - \langle \mu, c \rangle - \gamma \sum_{i,j=1}^n \exp\left(\frac{-C_{i,j} - \lambda_i - \mu_j}{\gamma}  - 1\right).
\end{equation}
Denote by \(\lambda^*_\gamma, \mu^*_\gamma\) the optimal variables for (\ref{eq:dual_simple_regularized}). Our goal is to tie dual variables of regularized and non-regularized problems, and it could be done using the following proposition.

\begin{prop}
    \((\lambda^*_\gamma, \mu^*_\gamma)\) converges to the optimal dual variables of the non-regularized case (\ref{eq:dual_simple}) as \(\gamma \to 0^+\).
\end{prop}
\begin{proof}
    Using Theorem 7.17 and 7.30 from \cite{rockafellar2009va}, it is sufficient to prove that the function under minimization in (\ref{eq:dual_simple_regularized}) converges point-wise to the function in (\ref{eq:dual_simple}).

    Notice that by the properties of the limit it is sufficient to prove that
    \[
        \gamma \exp\left( \frac{\alpha}{\gamma} - 1 \right) \underset{\gamma \to 0^+}{\to} \begin{cases}
            +\infty & \alpha > 0 \\
            0 & \alpha \leq 0
        \end{cases}
    \]
    as a function of \(\alpha\).

    If \(\alpha \leq 0\), it is clear. If \(\alpha > 0\) look at the term under the exponent:
    \[
        \frac{\alpha}{\gamma} - 1 + \log \gamma = \frac{\alpha + \gamma \log \gamma}{\gamma} - 1 \to +\infty
    \]\qed
\end{proof}

The corollary for our purpose is that it is sufficient to make a bound ''with respect to \(\gamma\)'' for \(\lambda^*_\gamma, \mu^*_\gamma\) and use the limit to create a bound on the dual variables of the non-regularized problem.
\paragraph{Bound on dual variables in the non-regularized case.}  First, we make the change of variables: \(u = -\nicefrac{\lambda}{\gamma} - \nicefrac{1}{2}, v = -\nicefrac{\mu}{\gamma} - \nicefrac{1}{2}\). Denote by \(K\) the matrix \(e^{C/\gamma}\). Then we can rewrite (\ref{eq:dual_simple_regularized}) in the next form using the notation of \cite{cuturi2013sinkhorn}:
\begin{equation}\label{eq:dual_regularized}
    L_{C}^\gamma(r,c) = \gamma \max_{v,u} \left[ \langle u, r \rangle + \langle v,c \rangle - \langle \mathbf{1},  \diag{e^u} K \diag{e^v}\rangle  \right] + 2 \gamma
\end{equation}

\textit{Remark.} If \((u^*, v^*)\) are optimal in (\ref{eq:dual_regularized}), then \((u^* + \alpha \mathbf{1}, v^* - \alpha \mathbf{1})\) are optimal too. Hence, we can assume \(\max_i v^* = \max_{i \in [n]}(-\nicefrac{\mu_\gamma^*}{\gamma} - \nicefrac{1}{2}) = -\nicefrac{1}{2} \iff \min_{i \in [n]} (\mu_\gamma^*)_i = 0\).

Now we can use the bound from Lemma 1 of the work \cite{dvurechensky2018computational} in the assumption \(r > 0, c > 0\):
\begin{gather*}
    \max_i (u^*)_i - \min_i (u^*)_i \leq R, \\
    \max_i (v^*)_i - \min_i (v^*)_i \leq R,
\end{gather*}
where \(R := - \log( \nu \min_{i \in [n]}\{r_i, c_i\} ), \nu = e^{-\nicefrac{\Vert C\Vert_\infty}{\gamma}}\).

Use the inverse change the variables, we have that
\[
    \max_{i} (\mu_\gamma^*)_i - \min_i (\mu_\gamma^*)_i \leq \Vert C \Vert_\infty - \gamma \log \min_{i\in[n]} \{r_i, c_i\}.
\]
We know that \(\min_{i \in [n]} (\mu_\gamma^*)_i = 0\). Hence, it could be rewritten as
\[
    \Vert \mu_\gamma^* \Vert_\infty \leq \Vert C \Vert_\infty - \gamma \log \min_{i\in[n]} \{r_i, c_i\}.
\]

By tending \(\gamma\) to \(0^+\) we obtain the following theorem:
\begin{theorem}\label{th:dual_bound}
    There is an optimal solution \(\mu^*\) for a problem (\ref{eq:dual_problem}) such that \( \Vert \mu^* \Vert_\infty \leq \Vert C \Vert_\infty \).
\end{theorem}

Also we have an immediate corollary of the Theorem \ref{th:dual_bound}:
\begin{corollary}\label{cor:variable_bound}
    Assume that \((r^*, f_\mu^*)\) is optimal variables in the saddle-point problem (\ref{eq:barycenter_basic}). Then, if \(r^* > 0\), there exists another optimal \(\hat{f_\mu}\) such that  \(\Vert \hat{f_\mu} \Vert_\infty \leq \Vert C \Vert_\infty\) \(P_\xi\)-almost surely.
\end{corollary}

Using this corollary \ref{cor:variable_bound}, we can reformulate problem \eqref{eq:barycenter_basic} in the following manner:
\begin{equation}\label{eq:barycenter}
    \min_{r \in \Delta^n} \E \mathcal{W}(r, \xi)= \min_{r \in \Delta^n} \sup_{f_\mu \in \mathcal{F}^b} \E\left[-\langle \lambda^*(f_\mu(\xi), C), r \rangle - \langle f_\mu(\xi), \xi \rangle \right],
\end{equation}
where \(\mathcal{F}^b = \{ f \colon \Delta^n \to \mathbb{R}^n \mid f\) is \(P_\xi\)-measurable, \(\Vert f \Vert_\infty \leq \Vert C \Vert_\infty\) \(P_\xi\)-a.s.\(\}\). We will use this form of the problem in the further analysis.

\section{Algorithms}
\label{sec:algorithms}

In this section, we provide algorithms for computation population Wasserstein barycenters using the saddle-point formulation (\ref{eq:barycenter}) on the assumption of existence a positive optimal \(r^*\).

\subsection{Finite support case}
\label{sec:fin_support}

The obvious way to deal with the problem with an infinite-dimensional part is to assume that it is finite-dimensional. Formally, we assume that the random variable \(\xi\) has a finite support \(\Supp(\xi) = \{c_1,\ldots, c_m\}\).

Then any function \(f_\mu \in \mathcal{F}^b\) can be represented as a matrix  of size \(m \times n\) and the value of a function defined \(f_\mu(c_i)\) as an \(i\)-th row of this matrix. The value of a function on points outside of \(\Supp(\xi)\) is zero. The corresponding matrix we denote as \(M(f_\mu)\).

To make the notation simpler, we say that we have a random variable \(\zeta\) over indices \([m]\), such that \(c_{\zeta} = \xi\). Also, we have \(f_\mu(\xi) = M(f_\mu)_{(\zeta)}\). Hence, if we rewrite \eqref{eq:barycenter} in terms of \(\zeta\), we can use only a matrix notation turn our problem to the following:
\begin{equation}\label{eq:dual_finite}
    \min_{r \in \Delta^n} \E \mathcal{W}(r, \xi) = \min_{r \in \Delta^n} \max_{\substack{M \in \Mat_{m \times n}(\mathbb{R}) \\ \Vert M \Vert_\infty \leq \Vert C \Vert_\infty}} \E_\eta \left[-\langle \lambda^*(M_{(\zeta)}, C), r \rangle - \langle M_{(\zeta)}, c_\zeta \rangle \right].
\end{equation}

This non-smooth non-strongly convex saddle-point problem can be solved using the mirror descent algorithm for saddle point problems. We refer to
\cite{bubeck2014convex,nemirovski2009robust} for additional details on this basic algorithm.

\begin{theorem}
    Algorithm \ref{alg:finite_support} outputs a pair \((\widetilde r, \widetilde M)\) that satisfies duality gap \eqref{eq:confidence_region_criteria} criteria with the precision \(\varepsilon\) and the level of confidence \(\sigma\) for the problem \eqref{eq:dual_finite} in
     \[
        N = O\left(\frac{n \max(n\log n, m) \Vert C \Vert_\infty^2}{\varepsilon^2}  \log^2\left(\dfrac{1}{\sigma}\right)  \right)
    \]
    iterations and its total complexity is \(O(nN)\).
\end{theorem}

 \begin{algorithm}\label{alg:finite_support}
        \KwData{\(N\) -- number of iterations;}
        \KwResult{\(\widetilde r\) -- approximation of barycenter;}
        \Begin{
            Set \(r = \left( \nicefrac{1}{n}, \ldots, \nicefrac{1}{n} \right) \in \Delta^n\)\;
            Set \(M = 0_{m \times n} \in \Mat_{m \times n}(\mathbb{R})\)\;
            Set \(\alpha = 2\log n, \beta = 4mn\Vert C \Vert_\infty \)\;
            Set \(\eta = \frac{2}{\Vert C \Vert_\infty \sqrt{8n^2 \log n + 16mn} \cdot \sqrt{5N}}\)\;
            \For{k = 1 to N} {
                Sample \(t\) from \(P_\zeta\)\;
                Sample \(s\) from \(U([n])\)\;
                Sample \(q\) from \(P(r)\)\;
                \(g_s :=  - n \cdot \max_{j \in [n]} (-C_{s,j} - M_{t,j})\)\;
                \(J_q = \argmax_{j \in [n]}(-C_{q,j} - M_{t,j}) \)\;
                \(h_t := -e_{J_q} + c_t,\)\;
                \(r_s := r_s\exp(- \alpha \cdot \eta \cdot g_s)\)\;
                \(r := r / (\sum_{i=1}^n r_i)\)\;
                \(M_{(t)} := M_{(t)} - \beta \cdot \eta \cdot h_t\)\;
                Project \(M\) onto \(\ell_\infty\) box\;
                \(\widetilde r := \frac{1}{k} r + \frac{k-1}{k} \widetilde r\)\;
            }

            Return \(\widetilde r\)\;
        }
        \caption{Finite support version}
\end{algorithm}

\textit{Remark 1.} \(\widetilde r\) will be an \(\varepsilon\)-approximation of the Wasserstein barycenter with probability at least \(1-\sigma\).

\textit{Remark 2.} If we conclude that \(\Vert C \Vert_\infty^2\) and \(\sigma\) are constants and \(m \geq n \log n\), we have a complexity \(O(n^2 m \cdot  \varepsilon^{-2})\). It is the same complexity as the Iterative Bregman Projections (IBP) algorithm \cite{nenna2014ibp,kroshnin2019complexity}.

\begin{proof}
    We have a problem of type
    \[
        \min_{r \in \mathcal{X}} \max_{M \in \mathcal{Y}} \E \Psi(r,M,\zeta),
    \]
    where \(\mathcal{X} := \Delta^n, \mathcal{Y} := \{ M \in \Mat_{m \times n}(\mathbb{R}) \mid \Vert M \Vert_\infty \leq \Vert C \Vert_\infty\}\). In our case
    \[
        \Psi(r,M,\zeta) = -\langle \lambda^*(M_{(\zeta)}, C), r \rangle - \langle M_{(\zeta)}, c_\zeta \rangle,
    \]
    where \(\lambda^*\) is defined in \eqref{eq:lambda_def}.

    To use the mirror descent algorithm, we have to define additional structures on the spaces \(\mathcal{X}\) and \(\mathcal{Y}\).
    \begin{itemize}
        \item For a prox-function for \(\mathcal{X}\), we choose an entropy setup: \(d_\mathcal{X}(r) = \sum_{i=1}^n r_i \log r_i\). We can easily bound a diameter \(R^2_\mathcal{X} = \log n\) and compute an prox-operator as an exponential weighting in \(O(n)\) time. This prox-function respects \(\Vert \cdot \Vert_1\) norm.
        \item For a prox-function for \(\mathcal{Y}\), we choose a default mirror map \(d_\mathcal{Y}(M) = \frac{1}{2} \Vert M \Vert_F^2\). Clearly, we have \(R^2_\mathcal{Y} = m \cdot n \cdot 2 \Vert D_p \Vert_\infty^2\). The projection can be computed in linear time.  This prox-function respects \(\Vert \cdot \Vert_F\) norm.
    \end{itemize}
    Hence, we can endow the space \(\mathcal{Z} : =\mathcal{X} \times \mathcal{Y}\) by the prox-function \(d_{\mathcal{Z}}(r,M) = \frac{1}{2 \log n} d_{\mathcal{X}}(r) + \frac{1}{4 mn \Vert C \Vert_\infty} d_{\mathcal{Y}}(M) \).

    Next, we should define our stochastic (sub)gradient oracles. For this purpose we will use the source of randomness a sample \(t \sim P_\zeta\). We have:
    \begin{align*}
        G(r, M, t) &= \partial_{r} \Psi(r, M, t) = -\lambda^*(M_{(t)}, C), \\
        H(r, M, t) &= \partial_{M_{(t)}} \left[-\Psi(r, M, t) \right] \\
        &= \sum_{i=1}^n r_i \partial_{M_{(t)}} \max_{j \in [n]} \left[- C_{ij} - M_{tj} \right] + c_t.
    \end{align*}
    On each iteration oracle \(H\) gives us \(t\)-th row of the full subgradient matrix (up to the probability of \(c_t\)) and this way can be associated with random coordinate descent.

    However, these subgradients can be calculated in \(O(n^2)\) time and it is too long for us. To fix that issue we apply a pure random coordinate technique to the \(G\) and same technique for \(H\): we will sample index \(i \in [n]\) from distribution associated with \(r\) and calculate only \(i\)-th term of the sum in the definition of \(H\) and the expectation is not changed and stay unbiased.

    Denote the uniform distribution over \([n]\) by \(U([n])\) and the distribution associated with \(r\) by \(P(r)\). Then we can write our final unbiased stochastic oracles using the knowledge of the subgradient of a maximum:
    \begin{align}\label{eq:g_finite}
        g(r, M, t, s) &= - n \cdot \max_{j \in [n]} (-C_{s,j} - M_{t,j}); \\
        \label{eq:h_finite}
        h(r, M, t, q) &= -e_{J_q(M_{(t)}, t)} + c_t,
    \end{align}
    where multiplication by \(n\) in \(g\) is required to unbiased property and \(J_q(\mu, t)\) is an index where \(\lambda^*_q(\mu,t)\) is maximized. These subgradients  can be computed in linear time because they require compute maximum only over \(n\) values.

    Then we should bound their dual norm
     \begin{align*}
        \Vert g(r, M, t, s) \Vert_{\mathcal{X}^*} &= n \left\vert \max_{j \in [n]} (-C_{s,j} - M_{t,j})  \right\vert \leq 2 n \Vert D_p \Vert_\infty, \\
        \Vert h(r, M, t, q) \Vert_{\mathcal{Y}^*} &\leq \Vert c_t \Vert_2 + \Vert e_{J_q(M_{(t)}, c_t)} \Vert_2 \leq 2.
    \end{align*}

    This, for the finite horizon and constant step-size \(\eta = \frac{2}{L\sqrt{5N}}\), where \(L^2 = 8n^2\log n\Vert C \Vert_\infty^2 + 16mn\Vert C \Vert_\infty^2\), we have the required guarenteed by Proposition 3.2 of \cite{nemirovski2009robust} to obtain an \(\varepsilon\) precision with a level of confidence \(\sigma\) if
    \[
        N \geq \frac{5(8n^2\log n\Vert C \Vert_\infty^2 + 16mn\Vert C \Vert_\infty^2)}{\varepsilon^2}\left(8 + 2 \log \frac{2}{\sigma}\right)^2.
    \]
    \qed
\end{proof}

\subsection{Background in the infinite-dimensional optimization}
\label{sec:back_inf_dim_opt}
The previous approach could be viewed as a coordinate descent for the offline problem. To deal with the purely online statement, we have to make less serious assumptions and optimize our function in a <<good>> subspace of \(\mathcal{F}_b\) (see \eqref{eq:barycenter} for the definition), where it can be done effectively. We will follow the idea with a kernel trick from \cite{genevay2016stochastic}. This idea also is widely distributed among the ML community \cite{mohri2018foundations,steinwart2008support}.

At first, we introduce a required background for the infinite-dimensional analysis, follows \cite{bogachev2020real}. The main instrument for us is a Frechet derivative, which defines almost like the derivative for finite dimensions, except for the special requirements on uniformness over the family of sets. For this special type of derivatives, we also have a chain rule, that we will actively use.

The main problem that this derivative is the element of the dual space that typically is not canonically isomorphic to our space, hence, we cannot produce proof of the correctness for gradient descent-like procedures as usual. The solution to this problem is to consider only Hilbert spaces and apply the Riesz representation theorem to produce steps toward gradients direction. In this case, we also have a notion of subgradient in the terms of an inner product (but of infinite dimension space), just as in the finite dimension case:
\[
    f(x) - f(x_0) \geq \langle f'(x_0), x - x_0 \rangle,
\]
that is the most crucial part for the proofs of the correctness of gradient-descent like procedures.

The only question is how to compute this infinite-dimensional object and how to compute the derivatives of linear functionals like \(\varepsilon_x(f) := f(x)\) over \(f\). For this purpose we will use additional assumption: we assume that the space \(\mathcal{H}\) is so-called \textit{a reproducing kernel Hilbert space}.

A Hilbert space \(\mathcal{H}\) of functions \(f \colon \mathcal{X} \to \mathbb{R}\) is called a reproducing kernel Hilbert space (RKHS) if there is a symmetric positive-defined function \(k \colon \mathcal{X} \times \mathcal{X} \to \mathbb{R}\) called a kernel, such that
\begin{enumerate}
    \item \(\forall x \in \mathcal{X} : k(\cdot, x) \in \mathcal{H}\);
    \item \(\forall f \in \Hilb : \langle f, k(\cdot, x) \rangle_\Hilb = f(x)\);
    \item \(\forall x,y \in \mathcal{X} : \langle k(\cdot, x), k(\cdot, y) \rangle_\Hilb = k(x,y)\).
\end{enumerate}

It can be proven \cite{mohri2018foundations} by an almost explicit construction that for each symmetric positive-defined function \(k\) there is such a Hilbert space. However, there is some kernels that are called \textit{universal} \cite{steinwart2008support}: they can uniformly approximate any continuous function if \(X\) is a compact subset of \(\mathbb{R}^d\). The example of a such kernel is a Gaussian RBF kernel \(k(x,x') = \exp(-s \Vert x - x' \Vert_2^2)\). The universality tells us that we do not lose too much by RKHS-assumption.

\subsection{General Kernel Mirror Descent}
\label{sec:general_kmd}

In this subsection we will apply the mirror descent to our paritially infinite-dimensional problem, using the kernel idea. Fix \((\mathcal{H},k)\) is a RKHS of functions from \(\Delta^n\) to \(\mathbb{R}\). Then we will search for an optimal \(f_\mu\) in the space \(\mathcal{H}^n\) of tuples of function from \(\mathcal{H}\), that will form a Hilbert space of functions from \(\Delta^n\) to \(\mathbb{R}^n\) with the inner product \(\langle f, g \rangle_{\Hilb^n} = \sum_{i=1}^n \langle f_i, g_i \rangle_\Hilb\), where \(f_i,g_i\) are \(i\)-th coordinate functions of \(f\) and \(g\) respectively.

Then we will solving to the new optimization problem, deduced from \eqref{eq:barycenter} using RKHS-assumption:
\begin{equation}\label{eq:barycenter_hilb}
    \min_{r \in \Delta^n} \E \mathcal{W}(r, \xi)= \min_{r \in \Delta^n} \sup_{f_\mu \in \mathcal{H}^{n}_b} \E\left[-\langle \lambda^*(f_\mu(\xi), C), r \rangle - \langle f_\mu(\xi), \xi \rangle \right],
\end{equation}
where \(\mathcal{H}^n_b = \{ f \in \mathcal{H}^n \mid \Vert f \Vert_\infty \leq \Vert C \Vert_\infty P_\xi\)-a.s.\(\}\) is an intersection of \(\mathcal{H}^n\) and \(\mathcal{F}^b\) (see \eqref{eq:barycenter}).

The theoretical dependence on the kernel is hidden in two following constants \(R^2 = \sup_{f} \frac{1}{2} \Vert f \Vert^2_\Hilb\) and \(\kappa^2 = \sup_{x} k(x,x)\) that fully depend on kernel. However, in practice, computation of \(R^2\) is a challenging problem.

\begin{theorem}
    For an arbitrary kernel \(k\), algorithm \ref{alg:general_kmd} (Kernel Mirror Descent) outputs the \(\varepsilon\)-approximation of the Wasserstein barycenter with respect to an arbitrary distribution \(P_\xi\) in term of the duality gap in
    \[
        N = O\left(\frac{n^2 \kappa^2 R^2 + \Vert C \Vert^2_\infty \log n}{\varepsilon^2} \log^2\left(\frac{1}{\sigma} \right) \right)
    \] sample iterations and with
    \[
        O\left(\frac{n^5 \kappa^4 R^4 + \Vert C \Vert_\infty^4 n \log^2 n}{\varepsilon^4} \log^4\left(\frac{1}{\sigma}\right)\right)
    \] total complexity.
\end{theorem}

\textit{Remark.}
    This algorithm is from a family of Stochastic Approximation (SA) algorithms and the most correct comparison can be obtained with other SA algorithms.

    From the article of \cite{dvinskikh2020stochastic} we know a total complexity of SA approach based on Sinkhorn algorithm:
    \begin{align*}
        \tilde O\biggl(& \frac{n^3 \Vert C \Vert^2_\infty}{\varepsilon^2} \min\biggl\{\frac{1}{\varepsilon}\sqrt{\frac{n}{\gamma \epsilon}},
        \exp\left(\frac{\Vert C \Vert_\infty \log n}{\varepsilon} \right)\left(\frac{\Vert C \Vert_\infty \log n}{\varepsilon} + \log \left(\frac{\Vert C \Vert_\infty \log n}{\gamma \varepsilon^2} \right) \right)\biggl\} \biggl).
    \end{align*}
    The complexity of this algorithm depends on the parameter \(\gamma\) which true dependence on \(n\) is unknown and in practice this parameter is quite small.
    From the theoretical point of view, if we suggest \(\Vert C \Vert_\infty, R\) and \(\sigma\) as constants, and \(\gamma \leq \nicefrac{\varepsilon}{n^3}\), the complexity of our algorithm will be better in both parameters.

    \begin{algorithm}\label{alg:general_kmd}
        \KwData{\(N\) -- number of iterations, \(r\) -- kernel, \(R^2,\kappa^2\) -- constants associated to the kernel;}
        \KwResult{\(\widetilde r\) -- approximation of barycenter;}
        \Begin{
          Set \(r = \left( \nicefrac{1}{n}, \ldots, \nicefrac{1}{n} \right) \in \Delta^n\)\;
          Set \(f_\mu = 0\)\;
          Set \(\alpha = 2\log n, \beta = 2nR^2\)\;
          Set \(\eta = \frac{2}{\sqrt{8 \log n \Vert C \Vert^2_\infty + 8 n^2 \kappa^2 R^2} \cdot \sqrt{5N}}\)\;
          \For{k = 1 to N} {
                Sample \(c^{(k)}\) from \(P_\xi\)\;
                \For{t = 1 to n} {
                    \((f_\mu(c^{(k)}))_t = \min\{1, \max\{-1, \sum\limits_{i=1}^{k-1} \beta^{(i)}_t k(c^{(k)}, c^{(i)}) \}\}\)\;
                }
                \For {i = 1 to n} {
                    \(J_i = \arg\max_{j \in [n]} (-C_{i,j} - f_\mu(c^{(k)})_{j})\)\;
                    \(g_i = -\max_{j \in [n]} (-C_{i,j} - f_\mu(c^{(k)})_{j}) \)
                }
                \For {t = 1 to n} {
                \[
                    \beta^{(k)}_t = \eta \cdot \beta \cdot \left(-c^{(k)}_t +  \sum_{i=1}^n r^{(k)}_i I\{t = J_{i}(c^{(k)}))\} \right)
                    \]
                }
                \(r := r \cdot \exp(-\eta \cdot \alpha \cdot g)\)\;
                \(r : = r / (\sum_{i=1}^n r_i)\)\;
                \(\bar r := \frac{1}{k} r + \frac{k-1}{k} \bar r\)\;
          }
          Return \(\bar r\)\;
        }
        \caption{Kernel Mirror Descent (KMD)}
    \end{algorithm}

\begin{proof}
As in the previous case, we are dealing with stochastic programming problem of the form
\[
        \min_{r \in \mathcal{X}} \max_{f_\mu \in \mathcal{Y}} \E \Phi(r,f_\mu,\xi),
\]
where \(\mathcal{X} := \Delta^n, \mathcal{Y} := \mathcal{H}^n_b\), and
\[
    \Phi(r,f,c) = -\langle \lambda^*(f(c), C), r \rangle - \langle f(c), c \rangle
\]

For \(\mathcal{X}\) we define entropic setup exactly as in the finite support case: there is no effect of an infinite dimension of \(f_\mu\). For \(\mathcal{Y}\) we define the prox-function \(d_{\mathcal{Y}}(f) = \frac{1}{2} \Vert f \Vert_{H^n}^2 = \frac{1}{2} \sum_{i=1}^n \Vert f_i \Vert_{\mathcal{H}}^2\). The range of the prox-function could be easily calculated: \(\mathcal{R}^2_\mathcal{Y} = n R^2\). To compute the prox-operator, we have to compute the Bregman divergence, using the notion of Frechet derivative:
\[
    B_\mathcal{Y}(f,g) = \frac{1}{2} \Vert g \Vert_{H^n}^2 - \frac{1}{2} \Vert f \Vert_{H^n}^2 - \langle f, g-f\rangle_{\mathcal{H}^n} = \frac{1}{2} \Vert f-g \Vert_{\mathcal{H}^n},
\]
as in the finite-dimensional Euclidean setting. Hence, the prox-operator is a just a projection onto \(\mathcal{H}^n_b\). Overall, for the product space \(\mathcal{Z} := \mathcal{X} \times \mathcal{Y}\) we have \(d_\mathcal{Z}(r, f_\mu) = \frac{1}{2\log n}d_\mathcal{X}(r) + \frac{1}{2n R^2}d_\mathcal{Y}(f_\mu)\).

Our current goal is to compute subgradients and to bound their dual norm. For \(f_\mu\) we compute gradient coordinate-wise, for each fixed coordinate \(t\).
\begin{align*}
    G(r, f, c) &= \partial_{r} \Phi(r, f, t) = -\lambda^*(f(c), C), \\
    H(r, f, c)_t &= \left(\partial_{f}(-\Phi(r,f(c),c))\right)_t \\
    &= \left( \partial_{f}(-\Phi(r,\langle f, k(\cdot,c) \rangle_\Hilb,c))\right)_t \\
    &= k(\cdot, c) \cdot \left( \partial_{\mu}(\langle \lambda^*(\mu, C), r \rangle + \langle \mu, c \rangle)\right)_t \\
    &= k(\cdot, c) \left( c - \sum_{i = 1}^n r_i I\{t = J_i(c)\} \right),
\end{align*}
where \(J_i(c)\) is  one of the indexes in \(\lambda^*_i\) where the maximum value is obtained. Here we can bound norms of both gradients:
\begin{align*}
    \Vert G(r,f,c) \Vert_\infty &\leq 2 \Vert C \Vert_\infty \\
    \Vert H_t(r, f_\mu, c) \Vert_\Hilb &\leq 2 \Vert k(\cdot, c) \Vert_\Hilb \leq 2 \kappa,
\end{align*}
where we used that \( \Vert k(\cdot, c) \Vert_\Hilb = \sqrt {\langle k(\cdot,c), k(\cdot, c) \rangle_\Hilb} = \sqrt{k(c,c)} \leq \kappa\). Hence, for each coordinate we have
\[
    \Vert H(r,f_\mu,c) \Vert_{\Hilb^n}^2 \leq 4 n \kappa^2.
\]

Hence, we can consider a constant step-size \(\eta = \frac{2}{L \sqrt{5N}}\), where \(L^2 = 8 \log n \Vert C \Vert_\infty^2 + 8 n^2 \kappa^2 R^2\) and we will have the required guarantees by \cite{nemirovski2009robust} if
\[
    N \geq \frac{5(8 \log n \Vert C \Vert_\infty^2 + 8 n^2 \kappa^2 R^2)}{\varepsilon^2}\left(8 + 2 \log \frac{2}{\sigma} \right)^2.
\]

The last step is to calculate the complexity of the computation of gradients. Note that we have to use actual value of \(f_\mu^{(k)}\) on the current iteration. For the step-size sequence \(\eta_k\), we have the explicit formula to calculate \(f_\mu^{(k)}\) (if \(f_\mu^{(0)} = 0\)).
\begin{prop}
    If \(f_\mu^{(0)} = 0\), we have a formula to calculate \(f_\mu^{(k)}\):
    \begin{align}\label{eq:beta_calculation}
    &(f_\mu^{(k)})_t = \sum_{i=1}^{k} \beta^{(i)}_t \cdot k(\cdot, c^{(i)}),
    &\beta^{(k)}_t = \eta_k \cdot \left(-c^{(k)}_t +  \sum_{i=1}^n r^{(k)}_i I\{t = J_{i}(c^{(k)}))\} \right),
\end{align}
where \(\{c^{(k)}\}\) are samples from \(P_\xi\).
\end{prop}
\begin{proof}
    Prove it by induction on \(k\): basis \(k=0\) is clear, step follows from the formula of a step of a gradient descent:
    \[
        (f_\mu^{(k)})_t = (f_\mu^{(k-1)})_t - \eta_k H_t(r, f_\mu^{(k-1)}, c^{(k)}) = (f_\mu^{(k-1)})_t + \beta^{(k)}_t k(\cdot, c^{(k)}).
    \]\qed
\end{proof}

     The main problem in the kernel approach can be seen there: we need \(O(N)\) iterations to recalculate one coordinate function of \(f_\mu\). For calculation of the subgradient in \(r\) we should calculate \textit{full} vector \(f_\mu\) for a new sampled measure \(c^{(k)}\), it costs \(O(nN)\). Hence, our total complexity is \(O(nN^2)\).
    \qed
\end{proof}

We will use this kernel approach with a Gaussian RBF kernel (\cite{minh2010})
\[
    k(x,x') = \exp\left(-s \Vert x- x'\Vert_2^2\right),
\]
and an information-diffusion kernel for simplex (\cite{lafferty2005diffusion})
\[
    k(x,x') = \exp\left(-\frac{1}{t}\arccos^2\left(\langle \sqrt{x} , \sqrt{x'}\rangle\right)\right).
\]
This kernels are very similar and their properites interested for us are similar too: \(\kappa^2 = \sup_{x} k(x,x) = 1\), and the value of \(R^2\) is unknown for both kernels. We will use it as a parameter for our algorithm. Note that the tuning of values \(s\) and \(t\) are crucial for both kernels.

\paragraph{Dynamic stepsize}

If we consider the convergence of duality gap in expectation \eqref{eq:expectation_criteria} instead of criteria in term of confidence region, we may use varying stepsize approach from  \cite{nemirovski2009robust} and obtain essentially the same theoretical complexity bounds. Thus, we can use stepsizes
\[
    \eta_t = \frac{\sqrt{2}\left( \sup_{r,r'}\frac{B_{\mathcal{X}}(r,r')}{2\log n} + \sup_{f,f'}\frac{B_{Y}(f,f')}{2nR^2} \right)^{1/2}}{L\sqrt{t}} = \frac{\sqrt{3}}{L\sqrt{t}}
\]
and use the sliding average \(\widetilde r = (\sum_{t=1}^N \eta_t r^{(t)})/(\sum_{t=1}^N \eta_t)\) as the output of our algorithm and obtain fully online algorithm for an infinte horizon.

\subsection{Linear Kernel Mirror Descent}
\label{sec:linear_kmd}
However, we have a kernel that allows easy computation of \(f_\mu\) and all related constants. The cost of this is a poor family of optimized functions: it is a linear kernel \(k(x,x') = \langle x, x' \rangle\) and a family of linear operators \(\Delta^n \mapsto \mathbb{R}^n\). For this kernel we can compute both constants:
\[
    \kappa^2 = \sup_{x \in \Delta^n} \langle x,x \rangle = 1, \ \ \ R^2 = \sup_{f_\mu \in \mathcal{H}^n_b} \Vert f_\mu \Vert_{2}^2 \leq 2 n^2 \Vert C \Vert_\infty^2,
\]
where we used a representation of a linear operator as a matrix of size \(n \times n\). Overall, we have the following complexity in terms of number of gradient computations:
\[
    N = O\left( \frac{n^4 \Vert C \Vert^2_\infty}{\varepsilon^2} \log^2\left(\frac{1}{\sigma}\right)\right).
\]
Moreover, each computation of the gradient can be done on \(O(n^2)\) time: we say that our function \(f_\mu\) is presented by square matrix \(\Theta\), and each gradient step could be represented in the space of matrices.

Overall, we have a total complexity
\[
  N = O\left( \frac{n^6 \Vert C \Vert^2_\infty}{\varepsilon^2} \log^2\left(\frac{1}{\sigma}\right)\right).
\]

It could be noticed that the theoretical complexity of this algorithm is very high, but the main advantage is a lack of parameters to tuning, and a high performance of each step in comparison as to KMD with a general kernel, as to methods based on gradients of the optimal transport distance \cite{dvinskikh2020stochastic}.

\textit{Remark.} All the complexities bounds in sections~\ref{sec:fin_support},~\ref{sec:general_kmd},~\ref{sec:linear_kmd} are based on typically rough enough bounds for the distance to the solution and dual norm of stochastic subgradients. We may expect that in practice we can use significantly smaller constants, that allows to converging faster due to the larger stepsizes. From this side, it can be useful to use adaptive variants of stochastic Mirror Descent (stochastic Universal Mirror Prox): \cite{bach2019universal} or \cite{dvurechensky2018generalized} with adaptive batching idea from \cite{dvinskikh2019adaptive}.

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

\subsection{Finite Support Algorithm.}
\label{sec:finite_support_exp}
The algorithm \ref{alg:finite_support} for measures with a finite support have a good theoretical properties, but it is very impractical because of a high variance of combinations of coordinate descent techniques. So we do not study it empirically, however, we may expect that it could be the only way to compute barycenters in a case of \(n \geq 10^6\) because of very cheap steps.

\subsection{Kernel Mirror Descent.}
\label{sec:kmd_exp}

We will study Kernel Mirror Descent algorithm \ref{alg:general_kmd} in three cases:
\begin{itemize}
    \item Gaussian RBF kernel \(k(x,x') = \exp\left( -s \Vert x - x' \Vert_2^2\right)\);
    \item Information-Diffusion kernel \(k(x,x') = \exp\left( -\frac{1}{t} \arccos^2\left( \langle \sqrt{x}, \sqrt{x'} \rangle\right) \right)\);
    \item Linear kernel \(k(x,x') = \langle x, x' \rangle\).
\end{itemize}

For the last one we implement more computational effective algorithm that works directly with matrices as briefly discussed in the section \ref{sec:linear_kmd}. For the first two algorithms we choose optimal parameter \(R^2\) empirically, and for the parameter of the kernel we perform a grid-search.

We will compare our approach with two algorithms, described in \cite{dvinskikh2020stochastic}:
\begin{itemize}
    \item Gradient descent with gradients of regularized Wasserstein distance, computed by Sinkhorn algorithm;
    \item Mirror descent with gradients of non-regularized Wasserstein distance, computed by solving a dual OT problem via linear programming;
\end{itemize}
We will refer to the first approach as a Sinkhorn-based and to the second as a direct LP-based.

\paragraph{Gaussian measures.}
To compare convergence we consider the auxiliary dataset of randomly generated Gaussian measures. The main reason is an ability to compute the true barycenter over the distribution \cite{delon2020wassersteintype}. We generate 10000 1-d Gaussian measures \(\mathcal{N}(\mu,\sigma^2)\) with \(\mu\) sampled from \(\mathcal{N}(1,4)\) and \(\sigma\) sampled from \(\operatorname{Exp}(1/2)\). The true population barycenter for this distribution over Gaussian measures is \(\mathcal{N}(1,4)\), since the expected value of \(\mu\) is 1 and the expected value of \(\sigma\) is 2.

After it, for each measure we consider it discretization over the segment \([-10;10]\) with \(n=300\) point masses. On the first step we consider distributions over mean and variance with much more general support to generate measures that could be called outliers, and make this model closer to applications.

As a measure of performance we choose \(\W_2\)-distance to true barycenter. We notice that we do not have convergence guarantees in this metric for our algorithms, however, it is one of the very limited number of ways to evaluate barycenters in the online approach where the expectation can not be expressed analytically.

Firstly, we found an optimal parameters for our general kernel methods. For both considered types of general kernels the optimal value was found empirically and is equal to \(R^2 = 25\). To find an optimal kernel parameters, we perform grid search and found out that the RBF kernel's optimal parameter \(s\) should be very small, and the information diffusion kernel's optimal parameter \(t\) should be very large. We choose \(s = 10^{-3}\) and \(t = 10^{3}\) by symmetry.

Secondly, we compare the measures of our two kernel methods with the Sinkhorn-based approach with two values of the regularization parameter \(\gamma = 10^{-2}\) and \(\gamma = 5 \cdot 10^{-5}\). For Gaussian and information diffusion kernels we obtain exactly the same barycenters, so, we illustrate only one of them on the figure \ref{fig:kmd_vs_sinkhorn}.

\begin{figure}
    \centering
    \includegraphics[scale=0.33]{images/kmd_vs_sinkhorn.png}
    \caption{Comparison between the KMD and the Sinkhorn-based method.}
    \label{fig:kmd_vs_sinkhorn}
\end{figure}

For Sinkhorn-based methods we faced numerical instabilities in the computation of gradients that may lead to numerical errors. Also note that the gradient descent with Sinkhorn gradients converges not to the true but to the regularized barycenter.

Then we compare our method with a Gaussian kernel to the direct method. Despite the fact of theoretical effectiveness, the running time of the direct method was so long that we compute the approximate barycenter using only 3000 measures instead of 10000. The comparison is presented on the figure \ref{fig:kmd_vs_lp}. Note that in this two methods we see stochastic noise and rather different <<structure>> of approximate barycenters.

\begin{figure}
    \centering
    \includegraphics[scale=0.33]{images/kmd_vs_lp.png}
    \caption{Comparison between the KMD and the direct LP-based method.}
    \label{fig:kmd_vs_lp}
\end{figure}

Next we compare barycenters computed with the linear KMD and the Sinkhorn-based method on the figure \ref{fig:lin_kmd_vs_sinkhorn}. Note that this algorithm clearly outperforms an algorithm with a relatively big value of \(\gamma = 0.01\), but not with a very small \(\gamma = 5 \cdot 10^{-5}.\) However, our algorithm is parameter-free, hence, we do not required to tuning \(\gamma\) and a number of Sinkhorn iterations to have a relatively good precision.

\begin{figure}
    \centering
    \includegraphics[scale=0.33]{images/lin_kmd_vs_sinkhorn.png}
    \caption{Comparison between the linear KMD and the Sinkhorn-based method.}
    \label{fig:lin_kmd_vs_sinkhorn}
\end{figure}

Finally, we compare all methods using proposed quality metric on the figure \ref{fig:w2_dist_measures}. Gaussian and information diffusion kernels shows almost equal performance, so we present only one of them. The best method in \(\W_2\) distance to the true barycenter is the direct LP-based method. However, the computation of gradients for this method spends too much time. Also we might observe that Sinkhorn-based method with small \(\gamma = 5 \cdot 10^{-5}\) stabilised on the one level. It could be connected to the convergence to the regularized barycenter or to numerical instabilities in gradients caused by a small value of \(\gamma\). Also we could see that approximate barycenters produced by the linear KMD outperforms Sinkhorm barycenters with \(\gamma = 0.01\) not only visually but in the \(\W_2\) distance, and barycenters produced by the KMD with Gaussian kernel become better than the Sinkhorn-based method with a large number of iterations.

\begin{figure}
    \centering
    \includegraphics[scale=0.33]{images/w2_distance_measures.png}
    \caption{Quality of the estimation in \(\W_2\)-distance to the true barycenter.}
    \label{fig:w2_dist_measures}
\end{figure}

Also we compare running time of each compared algorithm in the figure \ref{fig:processing_of_measures} except the direct LP-based method because it proceed 3000 measures in more then one hour (\(\approx 70\) minutes) when all other methods finished in less than 3 minutes. Here we could see one expected effect: the curves corresponding to the general KMD are parabolas whereas all other curves are lines. Another effect on this graph is connected to the numerical instabilities of a Sinkhorn algorithm with a small value of \(\gamma\): the computation of a gradient was stopped when there was a NaN value in the potential vector and the last computed value returned. Also we can see that the information diffusion kernel is slightly more computationally effective, than Gaussian.

\begin{figure}
    \centering
    \includegraphics[scale=0.33]{images/measure_processing_speed.png}
    \caption{Speed of processing of new measures.}
    \label{fig:processing_of_measures}
\end{figure}

\paragraph{MNIST.}
We will study the influence of the kernel parameters and the parameter \(R^2\) for the Kernel Mirror Descent algorithm with an information diffusion kernel by computing the barycenter of 1000 images of hand-written digits 3's from MNIST dataset.

We fix \(R^2 = 25\) and start varying the parameter \(t\) of the information diffusion kernel. The results could be seen in the figure \ref{fig:mnist_different_t}. For the values of \(t\) less than \(0.004\) the image disappears. We could say that the value of \(t \geq 1\) could be considered as a <<safe>> value where we obtain a good barycenter. However, there is a possibility to vary the <<thickness>> of the barycenter by tuning this parameter.

\begin{figure}
    \centering
    \includegraphics[scale=0.35]{images/mnist/kmd_inf_t=0004.png}
    \includegraphics[scale=0.35]{images/mnist/kmd_inf_t=0078.png}
    \includegraphics[scale=0.35]{images/mnist/kmd_inf_t=0336.png}
    \includegraphics[scale=0.35]{images/mnist/kmd_inf_t=1000.png}
    \caption{Comparison of the barycenter computed by the KMD with an information-diffusion kernel with different values of \(t\).}
    \label{fig:mnist_different_t}
\end{figure}

Then we fix \(t = 10\) and start varying parameter \(R\). The result of experiment is presented in \ref{fig:mnist_different_R}. We could describe that this parameter relies on the <<sharpness>> of the barycenter.

\begin{figure}
    \centering
    \includegraphics[scale=0.35]{images/mnist/kmd_inf_r=0001.png}
    \includegraphics[scale=0.35]{images/mnist/kmd_inf_r=0018.png}
    \includegraphics[scale=0.35]{images/mnist/kmd_inf_r=0336.png}
    \includegraphics[scale=0.35]{images/mnist/kmd_inf_r=1000.png}
    \caption{Comparison of the barycenter computed by the KMD with an information-diffusion kernel with different values of \(R^2\).}
    \label{fig:mnist_different_R}
\end{figure}

\section{Conclusion}
\label{sec:conclusion}

In this work, we provided two algorithms for a stochastic approximation approach for the Wasserstein barycenter problem.

The first algorithm uses the additional assumption that we take an expectation only over a finite number of measures. This assumption is very similar to usage of SAA-approach, when we replace our expectation with an empirical one. This algorithm meet the theoretical complexity of well-known Iterative Bregman Projections algorithm, however, from a practical side it suffers from high variance of the combination of coordinate descent-like techniques.

The second algorithm does not makes such huge assumptions and uses an optimization in reproducing kernel Hilbert spaces as a main instrument. It is the first algorithm for the online Wasserstein barycenter problem that does not uses a typically expensive or numerically unstable computation of gradients of the Wasserstein distance or its regularized versions. Moreover, in the experiment section we showed that this approach is practical and comparable to both main online approaches. Also usage of the simplest possible kernel give us a practical parameter-free algorithm.

\bibliographystyle{unsrt}
\bibliography{arxiv.bib}

config: arxiv


TITLE: Continuous probability distribution with no first moment but the characteristic function is differentiable
QUESTION [10 upvotes]: I am looking for an example of a continuous distribution function where the first moment does not exist but the characteristic function is differentiable everywhere.
Cauchy distributions do not fulfill this, as their characteristic functions are not differentiable at $0$. Does anybody here have an example?

REPLY [10 votes]: Below, there is a discrete example because I misread the word continuous.
However, if $X$ has the density $f$ given by
$$f(x)=\begin{cases}0&\mbox{ if }|x|\leqslant 2,\\
\frac C{x^2\log |x|}&\mbox{ if }|x|>2,\end{cases}$$
then $X$ is not integrable and the characteristic function is given by
$$\varphi(t)=2C\int_2^{+\infty}\frac{\cos(tx)}{x^2\log x}\mathrm{d}x.$$

Let us show that $\varphi$ is differentiable at $0$ and $\varphi'(0)=0$. We start from
$$\frac{\varphi(t)-\varphi(0)}{2tC}=\int_2^{1/t}\frac{\cos(tx)-1}{x^2\log x}\frac{dx}t+\int_{1/t}^{+\infty}\frac{\cos(tx)-1}{x^2\log x}\frac{dx}t=:A(t)+B(t).$$
Then using $|\cos u-1|\leqslant u^2$
$$|A(t)|\leqslant \int_2^{1/t}\frac{t^2x^2}{x^2\log x}\frac{dx}t=t\int_2^{1/t}\frac{dx}{\log x}.$$
Let $y:=\log x$, then $x=e^y$ and $dx=e^ydy$:
$$|A(t)|\leqslant t\int_{\log 2}^{-\log t}\frac{e^y}y\mathrm dy,$$
and splitting the interval at $R$, we get
$$|A(t)|\leqslant te^R/\log 2+R^{-1},$$
hence $\lim_{to\to 0}A(t)=0$.
Since
$$|B(t)|\leqslant \frac 1t\int_{1/t}^\infty\frac{2}{x^2\log(x)}dx,$$
we obtain after the substitution $s=1/x$ that
$$|B(t)|\leqslant -\frac 2t\int_0^t\frac{ds}{\log s}.$$
An other transformation $tu=s$ yields
$$|B(t)|\leqslant -\int_0^1\frac{du}{\log u+\log t},$$
which converges to $0$ as $t$ goes to $0$.
A similar argument can be used to show the differentiability at each point.

This example is taken from Stoyanov's book Counter-examples in probability.

Take $X$ such that $$\mathbb P(X=(-1)^kk)=\frac C{k^2\log k},\quad k\geqslant 2.$$
Then $\mathbb E|X|$ is infinite, while the characteristic function
$$\varphi_X(t)=\sum_{k\geqslant 2}\frac C{k^2\log k}e^{it(-1)^kk}$$
is differentiable everywhere: compute $\varphi_X(t+h)-\varphi_X(t)$ approximating $e^{ih(-1)^kk}-1$ by $h(-1)^kk+o(h)$.


TITLE: Domain of a composition of log functions
QUESTION [2 upvotes]: Domain of $$\log_3(\log_{1/3}(\log_4(\log_{1/4} x)))$$
Please guide me to solve this problem. Since it is a composition of logs, I am confused how to start.

REPLY [1 votes]: (not gonna discuss complex arguments that might give a real result, because I don't think that's intended in the question)
For the $log_3$ giving something real, the argument has got to be real and $> 0$. Thus
$$\log_{1/3}(\log_{4}(\log_{1/4}(x))) \ \in  \ (0,\infty)$$
which means
$$\log_{4}(\log_{1/4}(x))  \ \in  \ ((1/3)^\infty,(1/3)^0) = (0,1)$$
which is again
$$\log_{1/4}(x)  \ \in  \ (4^0,4^1) = (1,4)$$
and this means
$$x  \ \in  \ ((1/4)^4, (1/4)^1) = (1/256, 1/4).$$