\section{Introduction}

The barycenter in the space of probability measures equipped with the Wasserstein distance, first 
introduced by \cite{agueh_carlier_wbary} gained a lot of popularity in recent years.
The most striking property of the Wasserstein space is probably how its geodesics can be interpreted as a displacement of 
particles in $\RR^d$, allowing the Wasserstein barycenter of measures to take into account the geometry of the 
underlying space $\RR^d$.
Thanks to its geometric interpretation, it proved useful in a variety of applications.
For instance, it has been applied in image processing:
\cite{julien_texture_mixing} developed an algorithm for
texture synthesis built around this measure, \cite{barre_atm_wbary} used the
Wasserstein barycenter of images to improve the precision of gas emission sourcing,
and it played a central role in \cite{gramfort_foptand} to compare brain image data. 
In the field of fairness, \cite{barrio_fairness} and \cite{legouic_fairness} leveraged
the distribution to tackle the problem of fairness in regression and classification problems.
Barycenters also found applications in Bayesian inference; \cite{srivastava2015wasp,bayes_inference_example}
proposed a method to accelerate the computation of posterior distributions for Bayesian inference 
using the Wasserstein barycenter.

More formally, let us denote by 
$(P_2(\RR^d),W_2)$ the set of all measures defined on $\RR^d$ with finite second order moment, endowed with the Wasserstein distance
\begin{equation*}
    W_2 : (\mu, \nu) \longmapsto 
        \sqrt{\inf_{\gamma} \int_{\RR^d \times \RR^d} \nabs{x-y}^2 \Di\gamma(x,y)} \,,
\end{equation*}
where the infimum is taken over the set of couplings of $\mu, \nu$.
A Wasserstein barycenter $b$ of $\mu_1, \dots, \mu_n$ with weights $\lambda_1, \dots, \lambda_n$
is a minimizer of the map $\nu \mapsto \sum \lambda_i W_2^2(\mu_i, \nu)$,
\ie the measure $b$ corresponds to a Fréchet mean of $\mu_1, \dots, \mu_n$ in the Wasserstein
space $(P_2(\RR^d), W_2)$.
Equivalently, the multimarginal formulation of the barycenter
problem asserts that $b$ is obtained by pushing forward the minimizer $\gamma^\star$ of 
a functional $G$ (see equation \eqref{eqn:G_functional}) defined over the set of couplings of the 
$\mu_i$'s. Under mild conditions on the $\mu_i$'s,
both problems admit unique solutions, and yield the same measure.
Section \ref{sect:notions} provides more details on the Wasserstein barycenter.

Most methods in the literature propose to estimate the Wasserstein barycenter with a discrete measure \eg
\cite{cuturi2014fast,solomon2015convolutional,benamou2015iterative}.
However, this approach does not scale well with the dimension as noted by \cite{altschuler2021wasserstein}.
Indeed, they showed that computing the infimum of $\nu \mapsto \sum \lambda_i W_2^2(\mu_i, \nu)$ when the
$\mu_i$ are discrete measures is a NP-hard problem in the dimension, even
when only approximate solutions are acceptable.
We can, however, avoid estimating the density altogether
and focus on generating samples distributed according to the barycenter of known measures.
Given the broad applicability of the Wasserstein barycenter and of sampling techniques in general,
we believe that such sampling procedures should be part of the statistician's toolbox.

We are motivated by the success of sampling methods in high dimensional settings, which is
due to the possibility of integrating functions against the target measure without resorting to discretization on a large grid.
Markov chain methods have become popular tools to this end. 
The most well known among them are probably Hamiltonian Monte Carlo methods, variants of the 
Metropolis-Hastings algorithm relying on simulations to generate diverse samples from a distribution
(\cite[Chapter 11]{bishop_prml}), and Langevin diffusion
(\cite[Chapter 4]{pavliotis_stochastic_processes_and_apps_2015}), which transports
points along random trajectories to redistribute them according to the target measure. Such transportation
methods are computationally attractive since one only needs to store minimal information about how to
move the particles at each iteration. 

In their celebrated work, \cite{jko} showed that the marginals of Langevin
diffusion are distributed according to gradient flows in the Wasserstein space of the Kullback-Leibler divergence with
respect to the stationary measure. This insight, studied rigorously in \cite{ags}, brought to light a
connection between sampling and optimization which added a new perspective to the study of Monte Carlo
algorithms (\eg \cite{vempala_rates_ula} analyze the Unadjusted Langevin Algorithm and 
\cite{chewi2020exponential} analyze the Mirror Langevin diffusion in this framework). 
Inspired by this insight, we aim to minimize the multimarginal
formulation $G$ with gradient descent in the Wasserstein space. This procedure, like
Langevin diffusion, iteratively redistributes randomly initialized points to produce a sample from
the barycenter. 
We therefore recover an approximation of the optimal 
coupling through the samples, which can be desirable in applications. 

However since the barycenter is defined by a constrained optimization problem,
one cannot expect that the constraints will remain satisfied along the gradient flow.
This issue of constraints has been tackled in a variety of ways in the literature. A natural 
solution consists in restricting the domain of possible directions at each iteration,
which hints to the Frank-Wolfe algorithm, \ie to choose the steepest descent direction
in the admissible set.
This approach was studied in \cite{luise2019sinkhorn} within the
context of regularized optimal transport, where the regularized Wasserstein barycenter (also known
as the Sinkhorn barycenter) minimizes a sum of Sinkhorn divergences. 
When the admissible set of directions forms a Reproducing Kernel Hilbert Space (RKHS) with suitable kernel, 
\cite{sinkhorn_bary} propose to generate samples distributed 
according to the Sinkhorn barycenter by iterative pushforward of an initial measure with the map 
$\id_{\RR^d} - h \cdot d$, where $h$ is a step size and $d$ is the direction of steepest descent in the RKHS.
Both methods operate on discrete measures and the authors prove that the
continuous measure is recovered as a (weak) limit when the number of samples increases to ensure consistency.
We refer to \cite{peyre_cot} for details about regularized optimal transport.
We note that in some settings the issue of constraints can be addressed more easily, for example in the Bures-Wasserstein 
manifold (\ie the subspace of Gaussian measures in the Wasserstein space)
where the barycenter problem reduces to a finite dimensional optimization problem.
\cite{chewi_bures} propose a 
gradient descent algorithm for the original formulation of the barycenter problem.
However, we obviously cannot hope that similar properties hold for arbitrary continuous measures.

In this paper, we perform gradient descent on a penalized functional
$F^\al : P_2((\RR^d)^n) \to \RR^+$ obtained from $G$ by adding a penalization term to control the distance 
between the coupling marginals and the $\mu_i$'s.
We weigh the penalization with a coefficient $\al$ and control the induced error with $\al$. Taking advantage of
the differential calculus on the Wasserstein space, we can define a gradient flow for $F^\al$. 
To implement this procedure, we focus on the popular \ssf{SVGD} algorithm introduced in \cite{liu_wang_svgd}
to approximate the gradient of the penalization; technical details are given in Section \ref{sect:simulations}.


\begin{figure}[ht]
   
    \centering
    \vspace*{-1em}
    \includegraphics[ width=\textwidth]{Images/final_final_arctan.png}
    \vspace{-2em}
    \caption{\small Top plot: 200 samples drawn from $\bary{\prN(0,1), \arctan \pf \prN(0,1)}$ using
    \ssf{BARYGD} with a Gaussian kernel and initial step size $h = 10^{-3}$, and the resulting kernel 
    density estimate. Bottom plot: Trajectories of the samples evolving over 220 iterations, from green to
    red samples. We can see the $\arctan$ transport map is approximately reconstructed.}
    \label{fig:norm_arctan}
\end{figure}


\paragraph{Contributions.}

Inspired by the work of Jordan, Kinderlehrer and Otto cited above, we introduce a new sampling 
algorithm, called \ssf{BARYGD}, 
that performs kernelized gradient descent on a well chosen functional $F^\al$ built as a penalized version of 
the now classical multimarginal formulation of the barycenter problem introduced by \cite{agueh_carlier_wbary}.
We show that our method is consistent: with fixed penalization strength $\al$, if gradient descent yields 
a coupling $\geal$ such that $F^\al(\geal) - \inf F^\al \leq \epsilon$, then we get a quantitative bound on 
the Wasserstein distance between the approximate barycenter obtained from $\geal$ and the true barycenter 
obtained from the optimal coupling $\gamma^\star$.
This bound vanishes when $\epsilon \to 0$ and $\al \to \infty$.
As a consequence, we show that the minimizers of $F^\al$ converge to
$\gamma^\star$ when $\al \to \infty$ and quantify the rate of convergence. Furthermore, we 
perform numerical experiments with the algorithm in several settings (see Figures \ref{fig:norm_arctan} and \ref{fig:norm}).

\paragraph{Organization.}

In Section \ref{sect:notions}, we describe the Wasserstein barycenter in more details.
We present and analyze our method in Section \ref{sect:algorithm}.
Section \ref{sect:simulations} is devoted to implementation details. In Section
\ref{sect:discussion}, we discuss some open questions. 

\paragraph{Notation.}

\emph{All probability measures are assumed to be absolutely
continuous with respect to the Lebesgue measure and, with a slight abuse of notation, we 
will identify measures with their densities.}
We write $X\sim \nu$ when random variable $X$ has distribution $\nu$.
For $\mu, \nu \in P_2(\RR^d)$ and a measurable map $f : \RR^d \to \RR^d$,
we write $\nu = f \pf \mu$ if $f(X) \sim \nu$ when $X \sim \mu$. 
The set of couplings of $\mu$ and $\nu$, \ie probability measures on $\RR^d \times \RR^d$ with marginals
$\mu$ and $\nu$, is denoted by $\Pi(\mu, \nu)$.
We denote marginals by indices so $\gamma_1 = \mu$ and $\gamma_2 = \nu$ for $\gamma \in \Pi(\mu, \nu)$;
the notation extends naturally to couplings of $n$ measures.
We denote Euclidean gradients by $\nabla$
and write $\nabla_W$ for gradients in the Wasserstein space. Finally, we write indifferently $c_t = c(t)$ 
and $c' = \dot{c} = \partial_t c$ when $c : [a,b] \to X$ is a curve in some space $X$.
\section{Wasserstein barycenters}
\label{sect:notions}

This section briefly presents background notions on the Wasserstein space and its barycenters that were first 
introduced in the seminal paper of \cite{agueh_carlier_wbary}.
We refer to \cite{legouic_existence} for questions of existence and stability of the barycenter.

The Wasserstein space over $\RR^d$ is defined as the set $P_2(\RR^d)$ of probability measures over $\RR^d$ 
with finite second order moment, endowed with the distance $W_2$ defined by
\[
    W_2^2(\mu_0,\mu_1) 
        = \inf_{\gamma \in \Pi(\mu_0, \mu_1)} \int \nabs{x-y}^2\Di\gamma(x,y), \quad \forall \mu_0, \mu_1 \in P_2(\RR^d)
\]
where the infimum is taken over the set $\Pi(\mu_0, \mu_1)$ of all couplings $\gamma$ of $\mu_0$ and $\mu_1$.
The Wasserstein space is a geodesic space: for every $\mu_0$ and $\mu_1$ in $P_2(\RR^d)$ there exists 
a path $(\mu_t)_{t\in [0,1]}$ such that 
\[
    W_2(\mu_s,\mu_t) = |s-t| \, W_2(\mu_0,\mu_1), \quad \forall s,t \in [0,1].
\]
Such a path $(\mu_t)_{t\in[0,1]}$ is called a constant-speed \emph{geodesic} between its end points $\mu_0$ and $\mu_1$ 
--- we will often shorten it to \emph{geodesic}.
A functional $G$ defined over $P_2(\RR^d)$ is said to be \emph{geodesically convex} if it is convex along every geodesic.

Let $\mu_1,\dots,\mu_n\in P_2(\RR^d)$ and let $\lambda_1,\dots,\lambda_n$ be strictly positive weights 
such that $\sum \lambda_i = 1$. 
A Wasserstein barycenter of $(\mu_i)$ with weights $(\lambda_i)$ is defined as any measure
\begin{equation}
    \label{eqn:bary_def}
    b \in \argmin_{\nu \in P_2(\RR^d)} \sum \lambda_i W_2^2(\mu_i, \nu),
\end{equation}
\ie a barycenter corresponds to a Fréchet mean of the $\mu_i$ in the metric space $(P_2(\RR^d), W_2)$.
Whenever one of $\mu_1, \dots, \mu_n$ has density with respect to the Lebesgue measure, the barycenter is
unique and always defined. We assume throughout that the $\mu_i$ are absolutely continuous with
respect to the Lebesgue measure on $\RR^d$.

The minimization problem \eqref{eqn:bary_def} is equivalent to the following multimarginal problem.
Let 
$T(x) = \sum \lambda_i x_i$ for any $x \in (\RR^d)^n$ and recall that $\Pi(\mu_1, \dots, \mu_n)$ denotes the set of all couplings of 
the $\mu_i$.
The infimum of the functional
\begin{equation}
    \label{eqn:G_functional}
    G : \gamma \longmapsto \int \sum \lambda_i \nabs{x_i - T(x)}^2 \Di\gamma(x)
\end{equation}
over $\Pi(\mu_1, \dots, \mu_n)$ is achieved for a coupling $\gamma^\star$ such that
\begin{subequations}
    \begin{align}
        & b = \bary{\gamma^\star} \label{eqn:bary_pf}, \\
        & G(\gamma^\star) = \sum \lambda_i W_2^2(\mu_i, \bary{\gamma^\star}) \label{eqn:multmargform}.
    \end{align}
\end{subequations}
Moreover, if one of the $\mu_i$'s is absolutely continuous with respect to the Lebesgue measure, then $\gamma^\star$ is unique.
Note that by introducing the probability measure $P = \sum \lambda_i \delta_{\mu_i}$ on $P_2(\RR^d)$, the right hand side 
of \eqref{eqn:multmargform} corresponds to the variance of $P$.
With this notations, $b$ minimizes
the variance functional $\nu \mapsto \sum \lambda_i W_2^2(\mu_i, \nu)$ of $P$ and
$b$ corresponds to the barycenter of $P$.

In the next section, we leverage the smooth structure of the Wasserstein space in order to develop an optimization scheme 
for $G$ that will be the base of our sampling algorithm.

\section{Sampling with gradient flows}
\label{sect:algorithm}
This section defines the flow gradient of a functional on the Wasserstein space.
Some useful details of this technique introduced in \cite{jko} are gathered
in Appendix~\ref{sect:sosws}.

\subsection{Definition of the problem}

We aim to optimize the functional $G$ defined in Section \ref{sect:notions} 
with a gradient flow on the Wasserstein space. For simplicity, we first define the cost function $c$ by
\begin{equation*}
    c : x \mapsto \sum \lambda_i \nabs{x_i - T(x)}^2
\end{equation*}
then, for any $\gamma \in P_2((\RR^d)^n)$, we have
\begin{equation*}
    G(\gamma) = \int c \Di\gamma.
\end{equation*}

The set of marginal constraints $\Pi(\mu_1, \dots, \mu_n)$ is convex
in $L^2((\RR^d)^n)$ but is not geodesically convex in the Wasserstein space, which makes the multimarginal 
problem a difficult non-convex optimization problem on $P_2(\RR^d)$.
Since a gradient descent scheme on $G$ will inevitably leave the constraint set $\Pi(\mu_1, \dots, \mu_n)$, we modify the problem by enforcing 
the constraints with a penalization term that ensures that the marginals are close to the $\mu_i$'s in $\Penalization$ divergence.
We recall that the $\Penalization$ divergence between two measures $\mu, \nu$ is given by
\[
    \Penalization_\nu(\mu) = \int \bp{\frac{\mu}{\nu} - 1}^2 \Di\nu
        \textrm{  if } \mu \ll \nu \quad\textrm{and}\quad
    \Penalization_\nu(\mu) = +\infty \textrm{  otherwise}.
\]
The problem thus becomes to minimize the functional $F^\al$ defined by
\begin{equation*}
    \gamma \mapsto F^\al(\gamma) := \int c \Di\gamma + \al \sum \lambda_i \Penalization_{\mu_i}(\gamma_i),
\end{equation*}
for $\al\in[0,\infty]$ (we use the convention $0\times\infty=0$).

For $\al$ large enough, the minimizer of $F^\al$ should be close to the minimizer of $G$ with the marginal 
constraints --- which corresponds also to the minimizer of $F^\infty$.
The following proposition ensures that the minimizer exists.
\begin{proposition}
    \label{prop:barygd_falminexist}
    Suppose at least one of the $\mu_i$'s is absolutely continuous.
    Then, for any $\al > 0$, the functional $F^\al$ admits at least one
    minimizer in $P_2((\RR^d)^n)$.
    Moreover, this minimizer is absolutely continuous with respect to the Lebesgue measure.
\end{proposition}


Given a minimizer $\gamma^\al$ of $F^\al$, the measure $\bary{\gamma^\al}$, dubbed the associated barycenter to $\gamma^\al$, 
should be close to the barycenter $b$ of $P$ thanks to \eqref{eqn:bary_pf} and \eqref{eqn:multmargform}.
Our next result quantifies this proximity with respect to $\al$ under extra assumptions on $\mu_1, \dots, \mu_n$.

We first assume that $\mu_1, \dots, \mu_n$ satisfy the Poincaré inequality with a strictly positive constant $\Ciso$, \ie that
for any $i \in \eint{1,n}$ and any Lipschitz $f \in L^2(\mu_i)$, we have
\begin{equation}
    \label{def:logsobolev}
    \nabs{f}_{L^2(\mu_i)}^2 \leq \Ciso \nabs{\nabla f}_{L^2(\mu_i)}^2
\end{equation}
where $\nabla f$ is defined Lebesgue-almost everywhere.
Such inequalities are very common in the sampling setting.
We also assume  that $P = \sum \lambda_i \delta_{\mu_i}$ satisfies a variance inequality with constant $k$,
\ie that there exists a strictly positive constant $k$ such that for any 
$\nu \in P_2(\RR^d)$, it holds
\begin{equation*}
    k W_2^2(\nu, b_P) \leq B(\nu) - B(b_P),
\end{equation*}
where $b$ is the barycenter of $\mu_1, \dots, \mu_n$, and $B : \nu \mapsto \sum \lambda_i W_2^2(\mu_i, \nu)$
is the functional appearing in the original formulation of the barycenter problem. 

Such inequalities were introduced by \cite{sturm_varineq} to describe curvature properties of geodesic spaces and have played
a central role in understanding the behavior of the empirical Wasserstein barycenter 
(see \cite{ahidar_coutrix_rates,legouicFastConvergenceEmpirical2019}). We refer to Appendix \ref{sect:proofs} for 
details on variance inequalities. 
We can now state the following result.

\begin{proposition}
    \label{prop:barygd_barstab}
    Let $\al \geq 1$ and $\varepsilon \leq 1$.
    Suppose $\mu_1, \dots, \mu_n$ satisfy a Poincar\'e inequality with constant $\Ciso$ and that 
    $\sum \lambda_i \delta_{\mu_i}$ satisfies a variance
    inequality with constant $\Cvar$. Then there exists a strictly positive constant $C$, only depending
    on the variance $\sigma^2 = \int c \Di\gamma^\star$ and $\Ciso$, such that for any $\geal \in P_2((\RR^d)^n)$
    satisfying $F^\al(\geal) - \inf F^\al \leq \epsilon$, it holds
    \begin{equation}
        \label{ieqn:barygd_stab_wbars_bound}
        \frac{\Cvar}{4} W_2^2(\bary{\geal}, \bary{\gamma^\star}) \leq \epsilon + \frac{C}{\sqrt{\al}}.
    \end{equation}
\end{proposition}

Proposition \ref{prop:barygd_barstab} shows that any approximate minimizer of $F^\al$ produces a barycenter
which distance to the desired barycenter $b$ is of order $O(1 / \sqrt{\al})$ and can therefore be controlled
to arbitrary precision with $\al$.

The fact that the minimization of $F^\al$ is carried out without constraints on the marginals of $\gamma$ allows 
to use classical optimization techniques.
The next section is devoted to the introduction of the gradient flow on the Wasserstein space of $F^\al$ that will be the 
backbone of our algorithm.
\subsection{Minimizing scheme}


We first briefly recall that the Wasserstein gradient of a functional $F : P_2(\RR^d) \to \RR$ at 
a point $\nu \in P_2(\RR^d)$ is a vector field on $\RR^d$ given by the Euclidean gradient of its first variation $f$, \ie
\begin{equation}
    \label{eqn:ngrad_W}
    \nabla_W F(\nu) := x\mapsto\nabla f(x) \in L^2(\RR^d,\RR^d;\nu)
\end{equation}
where $f$ is such that for any signed measure $\delta$ with $\int \!\Di\delta = 0$, we have
\begin{equation*}
    F(\nu + \epsilon \delta) = F(\nu) + \epsilon \int f \Di\delta + o(\epsilon).
\end{equation*}
For instance, the Wasserstein gradient of $\mu \mapsto \Penalization_{\nu}(\mu)$ is given by 
\[
    \nabla_W \Penalization_{\nu}(\mu) = 2 \nabla \frac{\mu}{\nu}.
\]
Intuitively, when a probability measure is seen as a large number of particles, the Wasserstein 
gradient of a function defines a vector field along which each particle should be moved in order to locally maximize the function.
More details about the differential calculus on the Wasserstein space are provided in Appendix \ref{sect:sosws}.

To sample from the barycenter of $P=\sum \lambda_i\, \delta_{\mu_i}$, we consider the flow described by
\begin{equation}
    \label{eq:algo_cont}
    \begin{cases}
        \partial_t X_t = - \nabla_W F^\al(\gamma(t))(X_t) \\
        X(0) = X_0 \sim \gamma(0),
    \end{cases}
\end{equation}
where $\gamma(t)$ denotes the distribution of $X_t$ at time $t > 0$.
Recall that for a measure $\gamma \in P_2((\RR^d)^n)$, we write $\gamma_i$ for its $i$-th marginal.
Using \eqref{eqn:ngrad_W}, the Wasserstein gradient of $F^\al$ is given by
\begin{equation*}
    \nabla_W F^\al(\gamma)(x)
        = \nabla c(x) + 2 \al \sum \lambda_i \nabla \frac{\gamma_i}{\mu_i}(x),
\end{equation*}
for any $\gamma\in P_2((\RR^d)^n)$ and $x\in\RR^d$.

The explicit Euler scheme for \eqref{eq:algo_cont} provides an approximation of $(X_{t_m})_m$ for $t_{m+1}-t_m=h_m$, with
\begin{equation}
    \label{eqn:bad_euler}
    X_{m+1} = X_m - h_m \left[\nabla c(X_m) + 2 \al \sum \lambda_i \nabla \frac{\gamma_i(t_m)}{\mu_i}(X_m)\right],
\end{equation}
where $\gamma(t_m)$ is the distribution of $X_m$ and $h_m$ is a given step size.

However, note that since $\gamma(t)$ is unknown, this discretization scheme is not implementable as it is.
In order to implement it, we use a kernel to approximate the Wasserstein gradient
of the penalization term in $F^\al$ as is done in \cite{liu_wang_svgd, chewi_legouic_lu_maunu_rigollet_lawgd}.
The next section provides more details.


\section{Implementation}
\label{sect:simulations}

As noted in Section~\ref{sect:algorithm}, Algorithm \eqref{eqn:bad_euler} cannot be implemented directly 
since there is no canonical way to recover the unknown distribution $\gamma(t_m)$ of the current 
state $X_{t_m}$ from the mere knowledge of $X_m$. In this section, we present a concrete implementation of
the scheme \eqref{eqn:bad_euler}.


\paragraph{The \ssf{SVGD} algorithm.}

We consider the popular \ssf{SVGD}
algorithm introduced in \cite{liu_wang_svgd} which uses an explicit kernel $\ksvgd$ independent
of the target measure.
The kernel $\ksvgd$ can be chosen to have a convenient
analytical form (\eg a Gaussian kernel) amenable to direct evaluation. Let
$\pi \propto e^{-V}$, with $V \in \prC^1(\RR^d)$, be the target measure.
\ssf{SVGD} aims to minimize the $\Penalization_\pi$ divergence by transporting randomly initialized points
along the flow described by
\begin{equation}
    \label{eqn:svgd_flow}
    \partial_t X_t = - \iosvgd_{\pi} \nabla \Penalization_{\pi}(\mu_t)(X_t),
\end{equation}
where $X_t \sim \mu_t$ and $\iosvgd_{\pi} : f \mapsto \int \ksvgd(\scdot, x) f(x) \Di\pi(x)$ is the integral operator
associated with $\ksvgd$. Note that integration by parts yields
\[
    \iosvgd_{\pi} \nabla \Penalization_{\pi}(\mu_t)(X_t)
        = \int \bb{\ksvgd(X_t, x) \, \nabla V(x)  - \nabla_2 \ksvgd(X_t, x)} \Di\mu_t(x).
\]
Approximating the integrals by averaging with samples $(X^i_t)_{i=1,\dots,N}$ distributed according to $\mu_t$, 
one gets the \ssf{SVGD} algorithm
\begin{equation}
    \label{eqn:svgd_iter}
    X_{t+1}^i = X_t^i +
        \frac{h_t}{N} \sum_{j=1}^N  \Big\{ -\ksvgd(X_t^i, X_t^j)\, \nabla V(X_t^j)
            + \nabla_2 \ksvgd(X_t^i, X_t^j) \Big\}
\end{equation}
where the initial points $X_0^1, \dots, X_0^N$ are chosen at random.


\paragraph{Proposed algorithm.}

To minimize $F^\al$, we implement \eqref{eqn:bad_euler}, replacing each gradient of 
the penalization sum, $\sum_{i=1}^n\lambda_i\Penalization_{\mu_i}$, with the \ssf{SVGD} iteration for each marginal.
Thus, we couple $n$ batches of $N$ points and generate $N$ samples approximately distributed according to the
barycenter $\bary{\gamma^\al}$ of a coupling $\gamma^\al$ minimizing $F^\al$. Our iteration takes the form
\begin{equation*}
   
    X_{t+1}^{i,j} 
        = X_t^{i,j} - h_t \nabla_j c(X_t) 
           
            + \frac{\al h_t \lambda_j}{N} \sum_{\ell=1}^N 
                \bb{ \ksvgd(X_t^{i,j}, X_t^{\ell,j}) \, \nabla \log \mu_j(X_t^{\ell,j}) + \nabla_2 \ksvgd(X_t^{i,j}, X_t^{\ell,j}) }
\end{equation*}
where $j \in \eint{1,n}$ is the index of marginal $\mu_j$ and $X_t^{i,j}$, $i \in \eint{1,N}$, stands for the 
$i$-th sample associated with the $j$-th marginal.

\paragraph{Numerical details.}

All the simulations are carried out with a Gaussian kernel 
\begin{equation*}
    \ksvgd(x,y) = \exp\bp{-\nabs{x-y}^2} \quad \forall x, y \in \RR.
\end{equation*}
The step size $h_t$ as well as $\al = \al_t$ vary over iterations, starting from $h_0=0.1$ 
and $\al_0=1000$. We first keep the parameters fixed and let gradient descent concentrate the particles on 
the graph of an increasing function before doubling the penalization strength while keeping $h_t \cdot \al_t$
constant. Thus, we enforce the marginal constraints once the coupling corresponds to an optimal transport.
In Figure \ref{fig:norm}, we the step-size evolves according to AdaGrad (see \cite{duchi_adagrad}).
We give further experimental results in Appendix \ref{sect:app_simuls}.

\section{Conclusion}
\label{sect:discussion}

We have presented a new algorithm to sample from the Wasserstein barycenter of a finite set of measures.
Our numerical experiments illustrate its convergence in several examples. 
We also proved a theoretical bound for samples drawn from almost minimizers of $F^\al$, the gradient 
flow of which inspired our algorithm.
It remains to formally prove that this algorithm, a kernelized gradient descent of $F^\al$, converges to a 
minimum of $F^\al$ and bound the rate of convergence.
On another note, it would be interesting to import other common tools of optimization to the Wasserstein setting
in order to derive a sampling algorithm by optimizing $F^\al$.





\begin{figure}[ht!]
   
    \centering
    \hspace*{-2em}
    \includegraphics[width=0.8\textwidth]{Images/svgd_norms.png}
    \includegraphics[width=0.8\textwidth]{Images/distances_norms.png}
   
    \caption{\small \textit{Top figure}: Top plot: 150 samples drawn from the barycenter of two unit variance normal
        distributions (shown in light gray) using our algorithm with Gaussian kernel with $\al = 10^3$ and $h_m = 10^{-4}$, 
        and the resulting kernel density estimator. 
        Bottom plot: Trajectories of the samples evolving over 1000 iterations, starting from green to red sample.
        The transport map is almost perfectly reconstructed.
        \textit{Bottom figure}: the Wasserstein distance to the true barycenter over iterations plotted in $\log-\log$ scale. 
        From the graph, we read that \ssf{BARYGD} performs as $O(1/\sqrt{m})$ in the number of iterations $m$.
    }
    \label{fig:norm}
\end{figure}

\titleformat{\section}[block]{\sffamily\textbf{\MakeUppercase{#1}}}{}{1em}{}

\section*{Acknowledgments}

We thank Philippe Rigollet for useful discussions.
Thibaut Le Gouic was supported by NSF award IIS-1838071.

\clearpage
\printbibliography

\newpage
\onecolumn
