Hierarchical Barycenter Problem Overview

• Hierarchical barycenter problem is a multilevel aggregation framework that combines within-group barycenters with global optimization to preserve underlying structure.
• It leverages optimal transport, conditional simulation with independence penalties, and graph-based Fréchet means to address non-Euclidean and incomplete data challenges.
• The approach ensures well-posedness and computational tractability through entropy regularization, Jensen inequalities, and scalable multiscale algorithms.

The hierarchical barycenter problem denotes a family of multilevel aggregation problems in which barycentric structure is imposed across groups, covariate strata, or graph scales rather than only across a single flat collection of inputs. In the optimal-transport setting, the basic object is the Wasserstein barycenter, defined as a minimizer of a weighted sum of transport costs between a candidate measure and input measures; hierarchical variants first aggregate within groups and then across groups. In recent work, the same phrase also refers to a single-stage conditional simulation problem in which hierarchy is encoded through independence constraints with respect to multiple covariate subsets, including missingness or group indicators, and to a multiscale graph Fréchet-mean problem solved by coarse-to-fine refinement (Han et al., 2024, Tabak et al., 31 Jul 2025, Gavra et al., 2018).

1. Formalizations and scope

For a Polish metric space $\Omega \subset \mathbb{R}^d$ with ground cost $c(x,y)=\|x-y\|^p$ and associated $p$-Wasserstein metric $W_p$, the classical barycenter of measures $\{\mu_i\}_{i=1}^n \subset P_p(\Omega)$ with weights $\lambda_i \ge 0$, $\sum_i \lambda_i=1$, is

$$\mu^{*} = \arg\min_{\mu \in P_p(\Omega)} \sum_{i=1}^{n} \lambda_i W_p^p(\mu, \mu_i).$$

This problem is equivalent to a multi-marginal optimal transport problem, and existence holds generally, while uniqueness requires additional convexity or absolute continuity hypotheses (Tabak et al., 31 Jul 2025).

A grouped or nested hierarchical construction is formulated by first partitioning inputs into groups $G=1,\dots,M$, computing group barycenters

$$\nu_G \in \operatorname*{argmin}_{\nu \in \mathcal{P}_2(X)} \sum_j w_{G,j}\,W_2^2(\nu,\mu_{G,j}),$$

and then computing a global barycenter

$c(x,y)=\|x-y\|^p$0

The corresponding direct barycenter is

$c(x,y)=\|x-y\|^p$1

which aggregates all inputs in one step (Han et al., 2024).

The phrase “hierarchical barycenter” is not uniform across the literature. In one line of work it denotes the nested Wasserstein construction above; in another, it denotes a single optimization over a deformation map $c(x,y)=\|x-y\|^p$2 that removes dependence on multiple covariate subsets $c(x,y)=\|x-y\|^p$3 rather than a nested $c(x,y)=\|x-y\|^p$4 program; in graph settings it denotes a multiscale Fréchet-mean heuristic on clustered graphs (Tabak et al., 31 Jul 2025, Gavra et al., 2018).

Framework	Objective	Hierarchy mechanism
Wasserstein barycenter	$c(x,y)=\|x-y\|^p$5	Groupwise barycenters followed by a global barycenter
Conditional simulation HB	$c(x,y)=\|x-y\|^p$6	Independence from multiple covariate subsets and indicators
Graph multiscale barycenter	$c(x,y)=\|x-y\|^p$7	Clustering, coarse solve, then local refinement

A common misconception is that these formulations are interchangeable. The conditional simulation paper explicitly states that it does not use the nested $c(x,y)=\|x-y\|^p$8 formulation; instead, hierarchy is encoded by independence penalties with respect to selected covariate subsets, including a group or missingness indicator $c(x,y)=\|x-y\|^p$9 (Tabak et al., 31 Jul 2025).

2. Hierarchical Wasserstein barycenters on extended metric measure spaces

In the geometric theory developed for Wasserstein barycenters, the ambient space is an extended metric measure space $p$0, where $p$1 is complete, $p$2 is symmetric, satisfies the triangle inequality, and $p$3 iff $p$4, while $p$5 is a Radon probability measure with full support. The $p$6-Wasserstein space is $p$7, where $p$8 denotes Borel probability measures with finite $p$9-th moment relative to $W_p$0 (Han et al., 2024).

For a probability measure $W_p$1 on $W_p$2 with finite Wasserstein variance

$W_p$3

existence of a Wasserstein barycenter is proved under either of two hypotheses: $W_p$4 is $W_p$5 and $W_p$6 is concentrated on $W_p$7, or $W_p$8 is a probability measure and any $W_p$9 at finite distance to $\{\mu_i\}_{i=1}^n \subset P_p(\Omega)$0 is the starting point of an $\{\mu_i\}_{i=1}^n \subset P_p(\Omega)$1 gradient flow of $\{\mu_i\}_{i=1}^n \subset P_p(\Omega)$2 in Wasserstein space. This removes the usual local compactness constraint and covers non-compact and infinite-dimensional settings, including abstract Wiener spaces and certain $\{\mu_i\}_{i=1}^n \subset P_p(\Omega)$3 spaces (Han et al., 2024).

Uniqueness is obtained in several regimes. If $\{\mu_i\}_{i=1}^n \subset P_p(\Omega)$4 is $\{\mu_i\}_{i=1}^n \subset P_p(\Omega)$5, $\{\mu_i\}_{i=1}^n \subset P_p(\Omega)$6 is supported on $\{\mu_i\}_{i=1}^n \subset P_p(\Omega)$7, and $\{\mu_i\}_{i=1}^n \subset P_p(\Omega)$8, then the barycenter is unique. A second theorem gives uniqueness under a weak Monge property: if optimal transport from $\{\mu_i\}_{i=1}^n \subset P_p(\Omega)$9 to $\lambda_i \ge 0$0 is uniquely solvable whenever $\lambda_i \ge 0$1 and $\lambda_i \ge 0$2, then any $\lambda_i \ge 0$3 with finite mean entropy has a unique barycenter. The mechanism is strict convexity of $\lambda_i \ge 0$4 on the absolutely continuous locus (Han et al., 2024).

Absolute continuity follows from Jensen-type entropy bounds. If $\lambda_i \ge 0$5, then the barycenter $\lambda_i \ge 0$6 satisfies $\lambda_i \ge 0$7, hence $\lambda_i \ge 0$8. For finitely many inputs, multi-marginal optimal transport also yields uniqueness and absolute continuity under hypotheses such as $\lambda_i \ge 0$9 on $\sum_i \lambda_i=1$0 or all $\sum_i \lambda_i=1$1 on $\sum_i \lambda_i=1$2 (Han et al., 2024).

These results transfer directly to nested hierarchical barycenters. Under $\sum_i \lambda_i=1$3, $\sum_i \lambda_i=1$4, or the later $\sum_i \lambda_i=1$5 hypotheses, each group barycenter $\sum_i \lambda_i=1$6 and the global barycenter $\sum_i \lambda_i=1$7 exist; under finite entropy assumptions or the weak Monge property they are unique and absolutely continuous. However, no theorem in this framework asserts associativity, and no general equality

$\sum_i \lambda_i=1$8

is established outside special linear or Hilbert settings (Han et al., 2024).

3. Jensen inequalities, curvature-dimension conditions, and hierarchical consequences

A central structural result is an abstract Jensen inequality on extended metric spaces. If $\sum_i \lambda_i=1$9 is lower semicontinuous on an extended metric space $$\mu^{*} = \arg\min_{\mu \in P_p(\Omega)} \sum_{i=1}^{n} \lambda_i W_p^p(\mu, \mu_i).$$0 and every initial point admits an $$\mu^{*} = \arg\min_{\mu \in P_p(\Omega)} \sum_{i=1}^{n} \lambda_i W_p^p(\mu, \mu_i).$$1 gradient flow satisfying the integral form of $$\mu^{*} = \arg\min_{\mu \in P_p(\Omega)} \sum_{i=1}^{n} \lambda_i W_p^p(\mu, \mu_i).$$2, then for any probability measure $$\mu^{*} = \arg\min_{\mu \in P_p(\Omega)} \sum_{i=1}^{n} \lambda_i W_p^p(\mu, \mu_i).$$3 with finite variance and any barycenter $$\mu^{*} = \arg\min_{\mu \in P_p(\Omega)} \sum_{i=1}^{n} \lambda_i W_p^p(\mu, \mu_i).$$4,

$$\mu^{*} = \arg\min_{\mu \in P_p(\Omega)} \sum_{i=1}^{n} \lambda_i W_p^p(\mu, \mu_i).$$5

On Wasserstein space, taking $$\mu^{*} = \arg\min_{\mu \in P_p(\Omega)} \sum_{i=1}^{n} \lambda_i W_p^p(\mu, \mu_i).$$6 gives a Wasserstein Jensen inequality under $$\mu^{*} = \arg\min_{\mu \in P_p(\Omega)} \sum_{i=1}^{n} \lambda_i W_p^p(\mu, \mu_i).$$7 and in abstract Wiener or configuration-space settings (Han et al., 2024).

For barycenters $$\mu^{*} = \arg\min_{\mu \in P_p(\Omega)} \sum_{i=1}^{n} \lambda_i W_p^p(\mu, \mu_i).$$8 of $$\mu^{*} = \arg\min_{\mu \in P_p(\Omega)} \sum_{i=1}^{n} \lambda_i W_p^p(\mu, \mu_i).$$9, the entropy inequality takes the form

$G=1,\dots,M$0

On $G=1,\dots,M$1 spaces, a dimensional version is written in terms of

$G=1,\dots,M$2

and standard distortion coefficients $G=1,\dots,M$3 and $G=1,\dots,M$4:

$G=1,\dots,M$5

If $G=1,\dots,M$6, then $G=1,\dots,M$7 and $G=1,\dots,M$8 (Han et al., 2024).

The same work introduces the Barycenter-Curvature-Dimension condition. $G=1,\dots,M$9 requires that for any finitely supported $$\nu_G \in \operatorname*{argmin}_{\nu \in \mathcal{P}_2(X)} \sum_j w_{G,j}\,W_2^2(\nu,\mu_{G,j}),$$0, there exists a barycenter $$\nu_G \in \operatorname*{argmin}_{\nu \in \mathcal{P}_2(X)} \sum_j w_{G,j}\,W_2^2(\nu,\mu_{G,j}),$$1 such that

$$\nu_G \in \operatorname*{argmin}_{\nu \in \mathcal{P}_2(X)} \sum_j w_{G,j}\,W_2^2(\nu,\mu_{G,j}),$$2

If $$\nu_G \in \operatorname*{argmin}_{\nu \in \mathcal{P}_2(X)} \sum_j w_{G,j}\,W_2^2(\nu,\mu_{G,j}),$$3 and $$\nu_G \in \operatorname*{argmin}_{\nu \in \mathcal{P}_2(X)} \sum_j w_{G,j}\,W_2^2(\nu,\mu_{G,j}),$$4 is geodesic, this implies the classical $$\nu_G \in \operatorname*{argmin}_{\nu \in \mathcal{P}_2(X)} \sum_j w_{G,j}\,W_2^2(\nu,\mu_{G,j}),$$5 convexity of entropy along $$\nu_G \in \operatorname*{argmin}_{\nu \in \mathcal{P}_2(X)} \sum_j w_{G,j}\,W_2^2(\nu,\mu_{G,j}),$$6-geodesics. The authors emphasize that $$\nu_G \in \operatorname*{argmin}_{\nu \in \mathcal{P}_2(X)} \sum_j w_{G,j}\,W_2^2(\nu,\mu_{G,j}),$$7 is compatible with, but generally weaker than, full $$\nu_G \in \operatorname*{argmin}_{\nu \in \mathcal{P}_2(X)} \sum_j w_{G,j}\,W_2^2(\nu,\mu_{G,j}),$$8. A finite-dimensional $$\nu_G \in \operatorname*{argmin}_{\nu \in \mathcal{P}_2(X)} \sum_j w_{G,j}\,W_2^2(\nu,\mu_{G,j}),$$9 is defined analogously through the $c(x,y)=\|x-y\|^p$00 inequality (Han et al., 2024).

The $c(x,y)=\|x-y\|^p$01 class is stable under measured Gromov--Hausdorff convergence: if compact $c(x,y)=\|x-y\|^p$02 spaces $c(x,y)=\|x-y\|^p$03 converge in measured Gromov--Hausdorff sense to $c(x,y)=\|x-y\|^p$04, then the limit is again $c(x,y)=\|x-y\|^p$05. Existence of barycenters on $c(x,y)=\|x-y\|^p$06 spaces is then proved under finite variance and finite mean entropy, assuming either exponential volume growth with $c(x,y)=\|x-y\|^p$07 concentrated on $c(x,y)=\|x-y\|^p$08 or that $c(x,y)=\|x-y\|^p$09 is a probability measure (Han et al., 2024).

For hierarchical barycenters, Jensen inequalities propagate through the hierarchy. For any displacement convex functional $c(x,y)=\|x-y\|^p$10, including entropy,

$c(x,y)=\|x-y\|^p$11

This gives hierarchical upper bounds on entropy and related convex functionals. By contrast, direct quantitative bounds on $c(x,y)=\|x-y\|^p$12 are not provided (Han et al., 2024).

The barycentric curvature framework also yields geometric inequalities. On $c(x,y)=\|x-y\|^p$13 spaces, the support $c(x,y)=\|x-y\|^p$14 of the barycenter of the normalized restrictions of $c(x,y)=\|x-y\|^p$15 to bounded measurable sets $c(x,y)=\|x-y\|^p$16 satisfies the multi-marginal Brunn--Minkowski inequality

$c(x,y)=\|x-y\|^p$17

On $c(x,y)=\|x-y\|^p$18 spaces, measurable functions $c(x,y)=\|x-y\|^p$19 with $c(x,y)=\|x-y\|^p$20 integrable and satisfying

$c(x,y)=\|x-y\|^p$21

obey the functional Blaschke--Santaló-type inequality

$c(x,y)=\|x-y\|^p$22

These inequalities are derived from Jensen-type entropy convexity and the multi-marginal barycenter formulation (Han et al., 2024).

4. Independence-driven hierarchical barycenter for conditional simulation

A distinct formulation of the hierarchical barycenter problem arises in conditional probability density simulation with structured and unobserved covariates. Observations are pairs $c(x,y)=\|x-y\|^p$23 sampled from a joint law $c(x,y)=\|x-y\|^p$24, with $c(x,y)=\|x-y\|^p$25 and $c(x,y)=\|x-y\|^p$26, and the goal is to simulate $c(x,y)=\|x-y\|^p$27 for arbitrary $c(x,y)=\|x-y\|^p$28. The starting point is a distributional barycenter construction that seeks a transformation $c(x,y)=\|x-y\|^p$29 such that $c(x,y)=\|x-y\|^p$30 is independent of $c(x,y)=\|x-y\|^p$31 while minimally deforming $c(x,y)=\|x-y\|^p$32, through

$c(x,y)=\|x-y\|^p$33

with $c(x,y)=\|x-y\|^p$34 and typically $c(x,y)=\|x-y\|^p$35 (Tabak et al., 31 Jul 2025).

The hierarchical extension addresses the case in which covariates are defined only on subsets, overlap only partially across groups, or are missing for many samples. Let $c(x,y)=\|x-y\|^p$36 be augmented factors, and let $c(x,y)=\|x-y\|^p$37 be selected subsets that reflect hierarchy depth and availability patterns. The hierarchical barycenter solves

$c(x,y)=\|x-y\|^p$38

where the paper chooses mutual information,

$c(x,y)=\|x-y\|^p$39

$c(x,y)=\|x-y\|^p$40, $c(x,y)=\|x-y\|^p$41, and $c(x,y)=\|x-y\|^p$42 is the joint law of $c(x,y)=\|x-y\|^p$43 induced by $c(x,y)=\|x-y\|^p$44 (Tabak et al., 31 Jul 2025).

Hierarchy is encoded by the choice of subsets $c(x,y)=\|x-y\|^p$45. All components defined at a given group or node are included; a categorical indicator $c(x,y)=\|x-y\|^p$46 specifying group membership or missingness pattern is added so that between-group variability is also removed; and low-cardinality subsets such as singleton covariates can optionally be included for robustness when sample sizes are limited. The paper contrasts this construction with a two-level nested $c(x,y)=\|x-y\|^p$47 barycenter,

$c(x,y)=\|x-y\|^p$48

and explicitly states that this nested formulation is not the one used. Instead, hierarchy is represented by simultaneous independence from several covariate subsets and indicators in a single optimization over $c(x,y)=\|x-y\|^p$49 (Tabak et al., 31 Jul 2025).

The empirical problem is built from kernel density estimates. If $c(x,y)=\|x-y\|^p$50 denotes the set of indices on which all covariates in $c(x,y)=\|x-y\|^p$51 are observed, $c(x,y)=\|x-y\|^p$52, and $c(x,y)=\|x-y\|^p$53, $c(x,y)=\|x-y\|^p$54 are kernels, then

$c(x,y)=\|x-y\|^p$55

with

$c(x,y)=\|x-y\|^p$56

The sample objective becomes

$c(x,y)=\|x-y\|^p$57

and the third term in $c(x,y)=\|x-y\|^p$58 is constant in $c(x,y)=\|x-y\|^p$59, so it can be dropped without changing the minimizer (Tabak et al., 31 Jul 2025).

Optimization is performed by regularized gradient descent. If $c(x,y)=\|x-y\|^p$60 and $c(x,y)=\|x-y\|^p$61 is a diagonal curvature approximation with $c(x,y)=\|x-y\|^p$62, then the update is

$c(x,y)=\|x-y\|^p$63

with adaptive $c(x,y)=\|x-y\|^p$64. The gradient is simplified by ignoring derivatives with respect to kernel centers, which the paper describes as yielding a consistent estimator whose bias vanishes with sample size. For $c(x,y)=\|x-y\|^p$65, first-order optimality gives the inverse map

$c(x,y)=\|x-y\|^p$66

This formula is exact at the training samples and extends smoothly to arbitrary $c(x,y)=\|x-y\|^p$67, enabling conditional simulation by drawing barycenter samples $c(x,y)=\|x-y\|^p$68 and mapping them back via $c(x,y)=\|x-y\|^p$69 (Tabak et al., 31 Jul 2025).

Implementation uses empirical measures, Gaussian kernels in $c(x,y)=\|x-y\|^p$70 and diagonal Gaussian bandwidths in $c(x,y)=\|x-y\|^p$71, Silverman-type scaling for the $c(x,y)=\|x-y\|^p$72 bandwidths, and joint cross-validation over $c(x,y)=\|x-y\|^p$73 and the $c(x,y)=\|x-y\|^p$74-kernel bandwidths through held-out conditional log-likelihood. The penalties are parameterized as $c(x,y)=\|x-y\|^p$75 with $c(x,y)=\|x-y\|^p$76. A naive implementation has complexity $c(x,y)=\|x-y\|^p$77 per iteration, and the method is explicitly not Sinkhorn-based or entropic-OT-based (Tabak et al., 31 Jul 2025).

5. Empirical behavior, limitations, and interpretive issues in the conditional formulation

Empirical evaluation of the independence-driven hierarchical barycenter focuses on heterogeneous datasets with missing or partially overlapping covariates. In synthetic missing-covariate experiments with two covariates $c(x,y)=\|x-y\|^p$78 and three subsets of observations—one with only $c(x,y)=\|x-y\|^p$79, one with only $c(x,y)=\|x-y\|^p$80, and one with both—the method is compared against three baselines: a classical barycenter using only the fully observed subset, imputation followed by a classical barycenter, and a classical barycenter with full covariates serving as an oracle upper bound. Reported mean KL divergences are $c(x,y)=\|x-y\|^p$81, $c(x,y)=\|x-y\|^p$82, and $c(x,y)=\|x-y\|^p$83 for the hierarchical barycenter across three tests of increasing complexity, compared with $c(x,y)=\|x-y\|^p$84, $c(x,y)=\|x-y\|^p$85, and $c(x,y)=\|x-y\|^p$86 for the classical barycenter restricted to complete cases, $c(x,y)=\|x-y\|^p$87, $c(x,y)=\|x-y\|^p$88, and $c(x,y)=\|x-y\|^p$89 for imputation plus barycenter, and $c(x,y)=\|x-y\|^p$90, $c(x,y)=\|x-y\|^p$91, and $c(x,y)=\|x-y\|^p$92 for the oracle setting (Tabak et al., 31 Jul 2025).

On a bone mineral density dataset with synthetic missingness, where $c(x,y)=\|x-y\|^p$93 is bone density and $c(x,y)=\|x-y\|^p$94 are gender and age, held-out conditional log-likelihood averaged over $c(x,y)=\|x-y\|^p$95 repetitions is reported as $c(x,y)=\|x-y\|^p$96 for the hierarchical barycenter, $c(x,y)=\|x-y\|^p$97 for the complete-case classical barycenter, $c(x,y)=\|x-y\|^p$98 for the imputation baseline, and $c(x,y)=\|x-y\|^p$99 for the oracle with full covariates. The recovered densities are described as better capturing heteroscedasticity in some subpopulations (Tabak et al., 31 Jul 2025).

In structured-cofactor experiments with two groups and partially overlapping features, the gain depends on how informative the partially overlapping covariates are. For small values of the parameter $p$00, which controls how much one group informs the other, the hierarchical barycenter outperforms the classical barycenter computed only on the fully observed group, while gains diminish as $p$01 grows. In an extrapolation setting where one group does not cover a target region in $p$02, the hierarchical barycenter leverages the second group through the mutual-information constraints and improves extrapolation accuracy over the classical barycenter on the restricted group (Tabak et al., 31 Jul 2025).

The learned barycenter may also act as a residual-variability distribution. In an experiment with bimodal noise driven by an unobserved binary factor, the raw histogram of $p$03 is unimodal because the observed covariates mask the latent structure, whereas the hierarchical barycenter $p$04 reveals clear bimodality more distinctly than a standard barycenter estimated on the smaller fully observed subset. This suggests that removing observed covariate dependence before aggregation can expose latent structure, although that interpretation remains empirical rather than theorem-level (Tabak et al., 31 Jul 2025).

Several limitations are explicit. The paper does not prove new existence, uniqueness, or consistency results for the mutual-information-penalized hierarchical barycenter. The KDE-based mutual-information estimator is subject to the curse of dimensionality, the gradient approximation neglecting center derivatives introduces finite-sample bias, and large-scale computation may require mini-batching, subsampling, or low-rank kernel approximations. Future directions named in the paper include rigorous convergence and stability analysis, scalable solvers such as Nyström or random Fourier features, adaptive selection of $p$05, structured kernels for mixed data, and integration with entropic-OT solvers when appropriate (Tabak et al., 31 Jul 2025).

6. Multiscale graph barycenters as a hierarchical Fréchet-mean problem

On a finite connected weighted graph $p$06 with shortest-path metric $p$07 and node measure $p$08, the barycenter is the $p$09-Fréchet mean

$p$10

Here the object being averaged is a location on the graph, not a probability measure, and the metric is the graph shortest-path distance rather than the Wasserstein distance between distributions. The paper therefore differs fundamentally from Wasserstein barycenter theory even though it adopts hierarchical language (Gavra et al., 2018).

The baseline estimator is an online simulated-annealing method on the graph’s continuous analogue. Events $p$11 are i.i.d. from $p$12, update times are jump times of an inhomogeneous Poisson process with intensity $p$13, and an inverse temperature schedule $p$14 governs annealing. At each update, the process takes a stochastic Brownian-like step on the quantum graph and then moves deterministically toward the new observation along a shortest path, with scale $p$15. Logarithmic schedules $p$16 are described as theoretically safe, whereas linear schedules are faster in practice but can get trapped in local minima (Gavra et al., 2018, Gadat et al., 2016).

The multiscale extension uses a divide et impera strategy. The graph is partitioned into connected clusters $p$17 of nonzero mass. For each cluster, boundary information is recorded for edges crossing cluster boundaries, and a coarse graph $p$18 is built by representing each cluster with a single node $p$19. The representative may be the estimated barycenter of the cluster or a uniformly random node, and the paper reports that random representatives often suffice. An edge between $p$20 and $p$21 is assigned the minimum through-the-boundary path length consistent with the stored boundary distances. The cluster masses define the coarse probability $p$22, and events are projected by cluster membership (Gavra et al., 2018).

A coarse barycenter $p$23 is then estimated on $p$24, and the cluster $p$25 containing $p$26 is “opened up” to full resolution in a multiscale graph $p$27: the coarse node corresponding to $p$28 is replaced by all original nodes in $p$29, internal edges of $p$30 are restored, and boundary nodes of $p$31 are connected to neighboring coarse nodes. The online simulated-annealing estimator is then rerun on $p$32 to obtain the final barycenter estimate. The procedure realizes a hierarchical coarse-to-fine search rather than an exact decomposition theorem (Gavra et al., 2018).

The hierarchical lift is explicitly heuristic. The paper states that there is no theorem asserting that the coarse barycenter’s cluster necessarily contains the true barycenter. The method is motivated by the exact Euclidean decomposition under which the barycenter of a set equals the barycenter of block barycenters weighted by block masses, but this linear identity does not extend exactly to shortest-path graph geometry (Gavra et al., 2018).

Its principal contribution is scalability. The single-scale method requires an all-pairs shortest-path matrix of size $p$33, which becomes infeasible beyond medium-sized graphs. For a New York road graph with $p$34 nodes and $p$35 edges, the baseline would require about $p$36 GB of memory to store the full distance matrix, whereas the multiscale approach uses about $p$37 GB or less depending on clustering. Reported runtimes are about $p$38 h $p$39 min for a partition with roughly $p$40 clusters and about $p$41 h $p$42 min for a Markov-clustering partition with roughly $p$43 clusters; on a YouTube social graph with $p$44 nodes and $p$45 edges, a run takes about $p$46 hours (Gavra et al., 2018).

Empirically, the method matches the single-scale baseline closely on small graphs and remains stable on large graphs. Over $p$47 Monte Carlo runs, success rates are reported as $p$48 for single-scale, $p$49 for multiscale, and $p$50 for multiscale with random representatives on the Paris Metro graph; $p$51, $p$52, and $p$53 on a $p$54-node Facebook subgraph; and $p$55, $p$56, and $p$57 on a $p$58-node Facebook graph. On the New York graph, mean pairwise graph distances between returned centers over repeated runs are approximately $p$59 for one clustering and approximately $p$60 for another, which the paper interprets as good stability. On the YouTube graph, four runs produced two candidate centers at graph distance $p$61, one appearing three times and the other once (Gavra et al., 2018).

A plausible implication is that “hierarchical barycenter problem” serves as an umbrella for three technically distinct research programs: nested Wasserstein barycenters governed by curvature-dimension and entropy methods, independence-driven conditional simulation under partial covariate observability, and coarse-to-fine graph Fréchet-mean estimation. What unifies them is not a single universal objective, but the use of barycentric aggregation under hierarchical structure together with the need to preserve well-posedness, regularity, or computational tractability in spaces that are non-Euclidean, partially observed, or extremely large.