\bm{\hL{i}}^{(t)} - \eta \nabla \EnergyL{5}[\bm{\hL{i}}].$$ For completeness, we give the second-order derivative in Appendix. Its computation, however, involves integrating over the Voronoi facets and thus is intractable in general.

To solve for $\NU$, we rewrite the objective of the WB  as $$\label{eq:wb3}
\begin{split}
    \min_{\NU \in \MUS(\Y)} \EnergyL{6}[\NU]
    &\eqdef \frac{1}{\N} \sum_{i = 1}^{\N} \int_{\XL{i}}\|\x - \TLU{i}{*}(\x)\|_{2}^{2}d\MUL{i}(\x) \\
    & = \frac{1}{\N} \sum_{i = 1}^{\N} \sum_{\IdxCentroid = 1}^{\K} \int_{\RegionL{i, \IdxCentroid}}\|\x - \yL{\IdxCentroid}\|_{2}^{2}d\MUL{i}(\x),
\end{split}$$ $s.t.\ \yL{\IdxCentroid} = \TLU{i}{*}(x)\ \forall \x \in \XL{i}$. The critical point of this quadratic energy w.r.t. each $\yL{\IdxCentroid}$ has a closed form: $$\label{eq:vwb_dy}
\begin{split}
    \yLU{\IdxCentroid}{*}
    = \frac{\sum_{i = 1}^{\N} \int_{\RegionL{i, \IdxCentroid}} \x d\MUL{i}(\x)}{N\sum_{i = 1}^{N} \int_{\RegionL{i, \IdxCentroid}} d\MUL{i}(\x)} 
    \approx \frac{\sum_{i = 1}^{\N} \sum_{\x \in \RegionL{i, \IdxCentroid}} x \MUL{i}(\x)}{N\sum_{i = 1}^{\N} \sum_{\x \in \RegionL{i, \IdxCentroid}} \MUL{i}(\x)}, \nonumber
\end{split}$$ which is the center of mass of its correspondence across all measures. The latter expression is the “stochastic” version.

![Ten random nested ellipses (top) averaged according to the Euclidean distance (left) and the Wasserstein distance (right) as implemented by VWB. For a better visual, we use the Euclidean sum instead. Middle is the Euclidean sum after re-centered. The VWB preserves the topology (rainbow colors) of the ellipses.[]{data-label="fig:vwb"}](ed_vs_wd_color.pdf){width="0.9\columnwidth"}

The last step is to derive 