Generalized OT Problems

This chapter changes the optimization problem rather than only the ground distance. Barycenters average several measures, multi-marginal OT couples many measures at once, low-rank and capacity constraints restrict the admissible plans, inverse OT learns the cost from observed transport, and weak or martingale OT acts on conditional laws. These models remain close to Kantorovich optimization, but the unknown can now be a family of couplings, a factored plan, a learned cost, or a coupling subject to nonlinear conditional constraints.

from pathlib import Path
import sys

from IPython.display import Image as DisplayImage
from IPython.display import display

here = Path.cwd()
myst_dir = None
for candidate in [here, here.parent, here / "myst", here.parent / "myst", here.parent.parent / "myst"]:
    if (candidate / "ot4ml_web.py").exists():
        myst_dir = candidate.resolve()
        sys.path.insert(0, str(myst_dir))
        break

if myst_dir is None:
    raise RuntimeError("Could not locate myst/ot4ml_web.py")

repo_root = myst_dir.parent
thumbnails = repo_root / "notebooks-figures" / "thumbnails"

def show_book_figure(name, width=760):
    display(DisplayImage(filename=str(thumbnails / f"{name}.png"), width=width))

OT Barycenters¶

Barycenters ask how to average probability measures rather than points. This section explains the variational definition, the special closed forms in one dimension and for Gaussians, and the entropic algorithms used in practice.

Frechet Means¶

The natural formulation is a Frechet-mean problem on the space of probability measures: the unknown is the barycenter measure itself, and its support is not prescribed. It uses the continuous Kantorovich value $\mathcal L_c$ defined in (40).

Unlike a coupling, the barycenter is a new probability measure on $\Xx$. Since the weights $\lambda_s$ are nonnegative, problem (1) is convex in $\alpha$: Proposition Proposition: Convexity in the Marginals and Concavity in the Cost shows that the continuous Kantorovich value is jointly convex in its two marginals. Agueh and Carlier introduced this problem, following earlier ideas of Carlier and Ekeland Agueh & Carlier, 2011Carlier & Ekeland, 2010. For the quadratic cost on $\Xx=\RR^d$, a barycenter exists under the finite-second- moment assumption. It is unique if at least one positive-weight input is absolutely continuous; more general criteria ensure uniqueness through an essentially unique multi-marginal barycentric map. Discrete existence, consistency, and fixed-point constructions are studied in Anderes et al., 2016Esteban et al., 2016Le Gouic & Loubes, 2016.

Fixed-support discrete barycenters¶

For computation, one often turns the preceding infinite-dimensional problem into a finite one by prescribing possible barycenter locations. Assume the inputs are discrete,

Choose candidate barycenter sites $(y_i)_{i=1}^n$ and restrict the unknown to $\alpha=\sum_i a_i\delta_{y_i}$. For each input $s$, the cost $\mathcal L_c(\alpha,\beta_s)$ then becomes a finite Kantorovich problem with cost matrix

Thus its value is $\mathcal L_{C_s}(a,b_s)$, in the notation of the discrete Kantorovich problem (9).

This construction is a finite-dimensional restriction, not an exact discrete reduction of the general barycenter problem. In the ordinary two-marginal Kantorovich problem, once both marginals are discrete, the two supports are known and the whole problem is exactly the matrix optimization (9) on their product support. For barycenters, the input supports do not determine the support of the unknown barycenter: a minimizer may place mass outside the chosen sites $(y_i)_i$ and outside the union of the input supports. Once the candidate support is fixed, however, the nonnegative weights $\lambda_s$ and Proposition Proposition: Joint Convexity of Discrete OT show that problem (4) is convex in $a$: the discrete Kantorovich value is jointly convex in its two histograms.

For quadratic costs, the multi-marginal formulation of Section Multimarginal OT shows that, for discrete inputs, one may choose a barycenter supported on weighted averages $\sum_s\lambda_s x_{s,i_s}$ of one support point from each input. This exact candidate set can contain $\prod_s n_s$ points, but Corollary Corollary: Sparse Discrete Barycenters shows that there exists a barycenter for which at most $\sum_s n_s-S+1$ of them carry positive mass. Prescribing the support $(y_i)_i$ before solving (4) is nevertheless a numerical approximation, because the active weighted averages are not known in advance.

Figure Div moves beyond this degenerate case and compares barycenter grids obtained from one-dimensional quantile averaging and two-dimensional entropic transport under the same bilinear corner weights.

show_book_figure("barycenters-four-shapes")

Wasserstein barycenter grids for four corner measures. The left panel uses the one-dimensional formula $Q_{u,v}=\sum_{i,j}\lambda_{ij}(u,v)Q_{ij}$ for one Gaussian law and three asymmetric two-Gaussian mixtures, and displays densities reconstructed from the averaged quantiles. The right panel computes entropic Wasserstein barycenters on a common pixel grid for the cat, two-disk, cross and clover silhouettes, using the normalized squared ground cost, $\epsilon=4\cdot10^{-4}$ and a Sinkhorn tolerance of $5\cdot10^{-8}$. The barycenters are rendered as density images with values clamped at their $95\%$ quantile rather than by threshold contours. Colors interpolate between the four corners and encode the same bilinear weights in both panels.

The interactive demo below keeps the exact one-dimensional formula visible: the two coordinates set bilinear weights on the four corner laws, the middle panel averages their quantile functions, and the right panel reconstructs the resulting barycenter density.

One-Dimensional Case¶

On the line, barycenters become linear after the quantile change of variables. This gives the rare case where the barycenter is explicit rather than the solution of a high-dimensional optimization problem.

Gaussian Case¶

Gaussian barycenters show that the same separation as in the Gaussian Wasserstein formula persists: means average linearly, while covariances average according to the Bures--Wasserstein geometry.

If all input covariances are singular, the same contraction argument still gives a Gaussian barycenter, but the uniqueness step can fail and non-Gaussian barycenters may coexist. Thus nondegeneracy is essential when asserting that every barycenter is Gaussian.

Figure Div illustrates the nonlinear covariance interpolation characterized above; increasing anisotropy makes the simultaneous rotation and rescaling of the Bures--Wasserstein barycenter especially visible.

show_book_figure("barycenters-gaussian-covariances")

Bures--Wasserstein barycenters of centered Gaussian covariance matrices. Each panel shows a $5\times5$ grid of barycenter ellipses for four corner covariances, without separate input panels: the corner ellipses are the four input covariances themselves. The right grid uses more anisotropic inputs, making the nonlinear rotation and scaling of covariance barycenters more visible.

The interactive Gaussian demo compares the Bures covariance barycenter with a plain Euclidean covariance average under the same weights. The difference is most visible for rotated, anisotropic covariances: the Euclidean average blends matrix entries, whereas the Bures barycenter follows the geometry induced by quadratic Gaussian transport.

Sliced and Radon Barycenters¶

Slicing gives a scalable surrogate for high-dimensional barycenters by applying the one-dimensional quantile formula in every projection direction. For measures on $\RR^d$, one replaces $\Wass_2$ by the sliced distance $\SW_2$ introduced in Definition Definition: Sliced Wasserstein Distance and interpreted through the Radon transform in Section Sliced Wasserstein Distances:

The constraint that all projected measures come from the same $\alpha$ is the nontrivial part. A cheaper Radon-domain approximation drops this consistency constraint and minimizes directly over one-dimensional projected laws $(\gamma_\theta)_\theta$:

For each $\theta$, this is a one-dimensional barycenter, hence its quantile is the weighted average of the projected quantiles. For two inputs $\beta_0$ and $\beta_1$, define

Thus $Q_t$ is the directionwise quantile field, and when $\gamma_{t,\theta}$ has a density we denote it by $h_t(\theta,\cdot)$. The relaxed value is a lower bound on the sliced-barycenter value. If the minimizing family is Radon-consistent, meaning that $\gamma_\theta=(P_\theta)_\sharp\bar\alpha$ for a common probability measure $\bar\alpha$ and almost every $\theta$, then $\bar\alpha$ is an exact sliced barycenter. In general, independently computed one-dimensional barycenters do not satisfy the range conditions of the Radon transform. One therefore reconstructs a density in a least-squares sense, usually through a regularized Radon pseudoinverse.

Let $h(\theta,t)$ denote a density of $\gamma_\theta$. We use the one-dimensional Fourier transform in $t$ given by

whenever Fourier inversion is valid.

Formula (27) is the filtered back-projection representation of the Radon pseudoinverse used in tomography Herman, 1980; $|\omega|^{d-1}$ is its ramp multiplier. Since this multiplier amplifies high frequencies, one typically chooses a bandwidth $\Omega>0$ and an even low-pass window $\chi$ with $\chi(0)=1$, and replaces the ramp by

The numerical reconstruction below uses the super-Gaussian window $\chi(s)=e^{-|s|^4}$. Choose $\eta_t\geq0$ so that the positive part below has nonzero mass, and define the nonnegative, unit-mass reconstruction

For the endpoints, set $A_i=\rho_i$ when $\beta_i=\rho_i\d x$, $i\in\{0,1\}$. In the figure, the small threshold $\eta_t$ only suppresses finite-angle inversion ghosts. The resulting regularized density is generally only a least-squares approximation to the independently averaged slices. This fast construction was introduced for sliced and Radon Wasserstein barycenters in Bonneel et al., 2015, but it is not the exact constrained sliced barycenter. With all directions, the Radon transform is injective by the Cramér--Wold theorem Cramér & Wold, 1936; inconsistency comes from failure of Radon range conditions such as antipodal symmetry and moment consistency; with finitely many directions, the sampled Radon operator is also non-injective.

Figure Div follows the resulting density, projected-density and quantile fields through a cat-to-heart interpolation.

show_book_figure("sliced-radon-barycenter")

Radon-domain sliced barycentric interpolation between the cat and heart densities. The columns correspond to $t=0,0.2,\ldots,1$. The first row shows the endpoint densities and the intermediate reconstructions $A_t$ defined in (32) from the windowed pseudoinverse (31). The second row shows the projected-density fields $h_t$ (labeled $R_t$ in the figure), obtained by converting the directionwise quantile barycenters back into one-dimensional densities. The third row shows the quantile fields $Q_t$ defined in (23).

Sinkhorn for Barycenters¶

A key difference with the regularized two-marginal OT problem is that there is no canonical reference measure $\alpha\otimes\beta$, because the barycenter $\alpha$ is unknown. To reduce complexity, one usually fixes a candidate support for the barycenter and solves the discrete problem (4); this introduces a discretization error but keeps the number of unknowns manageable.

One can then use the entropy-only convention of (2) and approximate (4) by

for some $\epsilon>0$. This is a smooth convex minimization problem, which can be tackled using gradient descent Cuturi & Doucet, 2014. An alternative is to use a descent method, typically quasi-Newton, on the semi-dual Cuturi & Peyré, 2016; this is useful when adding extra regularization on the barycenter, for instance to impose smoothness.

A simple but effective approach developed in Benamou et al., 2015 observes that (33) has the same minimizers as the weighted KL projection problem

subject to

Here $K_s\eqdef e^{-\C_s/\epsilon}$. The barycenter $a$ is implicitly encoded in the common row marginal

The two objectives differ only by constants depending on $(\C_s,\epsilon)_s$, not on the couplings or barycenter. Assume below that every $\C_s$ is finite and every $b_s$ is positive; zero-weight target atoms can be deleted before the iteration. The optimal couplings then have scaling form

and the generalized Sinkhorn iterations are

The geometric mean enforces the fact that all couplings share the same barycenter marginal.

The scaling cycle has an exact dual-optimization interpretation; it is not merely a sequence of marginal normalizations. Write $u_s=e^{f_s/\epsilon}$ and $v_s=e^{g_s/\epsilon}$, with all operations understood componentwise. The proposition below shows that these potentials maximize the concave dual (41). With $(f_s)_s$ fixed, the objective separates over $s$, and its exact maximizer in each $g_s$ is

which is precisely the $v_s$-update. Conversely, with $(g_s^+)_s$ fixed, exact maximization over the coupled block $(f_s)_s$ under $\sum_s\lambda_s f_s=0$ gives

Thus a complete generalized Sinkhorn cycle is exact two-block coordinate ascent on the dual, or equivalently alternating minimization of its negative. Indeed, the constraint follows from $\log a^+=\sum_s\lambda_s\log q_s$, while the first-order conditions require $e^{f_s^+/\epsilon}q_s=a^+$ for every $s$. In the primal formulation (34), the same cycle alternates weighted KL projections onto the target-column constraints and the common-row-marginal constraint Benamou et al., 2015.

Classical applications include two-dimensional image interpolation, three-dimensional shape interpolation, and barycenters on surfaces where the ground cost is the square of the geodesic distance Solomon et al., 2015.

Wasserstein-Over-Wasserstein and Barycenters¶

The barycenter formula does not require a finite list of inputs. The Wasserstein-over-Wasserstein viewpoint of Section Wasserstein Over Wasserstein allows one to replace the discrete family $(\beta_s,\lambda_s)_s$ by a law $\mathfrak A\in\mathcal P(\mathcal P(\Omega))$ over probability measures. Such population Wasserstein barycenters were studied for random probability measures by Bigot and Klein, in general geodesic settings by Le Gouic and Loubes, and on Riemannian manifolds by Kim and Pass Bigot & Klein, 2018Le Gouic & Loubes, 2016Kim & Pass, 2017. Assume, for instance, that $c$ is lower semicontinuous and that there exists at least one $\beta_0\in\mathcal P(\Omega)$ with

together with the usual compactness or coercivity hypotheses ensuring existence of minimizers. For instance, these assumptions are automatic when $\Omega$ is compact and $c$ is continuous. The barycenter correspondence is then

When this set is a singleton, we denote its element by

which defines a nonlinear flattening map $\mathfrak A\mapsto\widetilde\alpha_{\mathfrak A}$ from laws over measures back to measures on $\Omega$. When $\mathfrak A=\sum_s\lambda_s\delta_{\beta_s}$, this is exactly the finite barycenter problem above. This map should be contrasted with the linear collapsed, or barycentric, mixture $\bar\alpha_{\mathfrak A}$ of Definition Definition: Collapsed, Or Barycentric, Mixture, which simply averages the input measures themselves. The two operations agree in degenerate linear situations, but in general $\widetilde\alpha_{\mathfrak A}$ is a geometric average in transport space, whereas $\bar\alpha_{\mathfrak A}$ is an ordinary mixture in the ambient linear space of measures.

The next result records the corresponding law of large numbers. It is useful when a dataset is itself made of probability measures, for instance populations of histograms, posterior distributions or shapes. We state it in the compact setting, where no moment or tightness side conditions are needed; non-compact extensions require the usual integrability assumptions. Consistency of Wasserstein barycenters and related statistical constructions is developed in Boissard et al., 2015Le Gouic & Loubes, 2016Zemel & Panaretos, 2019; streaming and large-scale uses of many input measures appear for instance in Staib et al., 2017Srivastava et al., 2015Srivastava et al., 2018.

Thus, (50) is the classical law of large numbers on Wasserstein space, and (51) is its linear image under the collapse map. By contrast, (52) is nonlinear: it recomputes a Wasserstein barycenter from the empirical law over measures. The number $p$ of input measures should not be confused with the number $n$ of samples used to approximate each input measure, studied in Section Sample Complexity. In applications one often observes $p$ empirical measures, each made of roughly $n$ atoms, hence about $np$ points in total. Balancing the error due to finitely many input laws against the error due to finitely sampled input laws is a separate statistical and computational tradeoff.

Toward Central Limit Theorems on Wasserstein Space¶

The same hierarchy suggests a central-limit refinement of the preceding law of large numbers, but the nonlinear geometry makes this substantially more delicate. For the linear collapse $\bar\alpha_{\widehat{\mathfrak A}_p}$, testing against a fixed $f\in C(\Omega)$ reduces the question to the classical scalar central limit theorem for the random variable $\alpha\mapsto\int f\,\mathrm d\alpha$. For the nonlinear barycenter $\widetilde\alpha_{\widehat{\mathfrak A}_p}$, however, there is no canonical vector difference $\widetilde\alpha_{\widehat{\mathfrak A}_p}-\widetilde\alpha_{\mathfrak A}$ inside $\mathcal P(\Omega)$. One has to choose a local linearization. When the population barycenter is sufficiently regular so that the optimal map $T_p$ from $\widetilde\alpha_{\mathfrak A}$ to $\widetilde\alpha_{\widehat{\mathfrak A}_p}$ exists, this amounts to asking whether $\sqrt p\,(T_p-\mathrm{Id})$ converges in a Hilbert space such as $L^2(\widetilde\alpha_{\mathfrak A})$. In nonsmooth settings one must instead work with optimal-plan or logarithmic-map coordinates. Even after such a linearization, an infinite-dimensional CLT requires tightness of the rescaled tangent variables and a genuine Radon Gaussian random element; in a Hilbert space, the associated covariance must be trace class. A cylindrical Gaussian limit alone is therefore not a probability law on the tangent space. This obstruction explains why Wasserstein-space CLTs are more rigid than the weak laws above.

There are nevertheless important settings where such results can be proved. In one dimension, the quantile representation linearizes $\mathcal W_2$, so barycenter fluctuations can be studied through empirical averages of quantile functions. Another finite-dimensional case is the family of non-degenerate Gaussian measures in fixed dimension, where $\mathcal W_2$ reduces to the Bures geometry of means and covariance matrices. Agueh and Carlier Agueh & Carlier, 2017 formulate this Wasserstein-barycenter CLT precisely in tangent coordinates and prove it in a few special cases, including the one-dimensional non-atomic setting and finite laws supported on non-degenerate Gaussian measures. Entropic barycenters give a smoother variant for which central-limit theorems for empirical barycenters are also available Carlier et al., 2021. These results should be read as nonlinear analogues of the statistical limits discussed in Chapter Paragraph, not as a generic Hilbert-space CLT valid on all of $\mathcal P(\Omega)$.

Multimarginal OT¶

Multi-marginal OT couples more than two measures at once. It is the natural language for barycenters, matching with teams and several-body costs, but its tensor dimension is the main computational obstacle.

Definition and Basic Structure¶

The multi-marginal formulation replaces a coupling between two measures by a joint distribution with several prescribed marginals. Given measures $(\alpha_s)_{s=1}^S$ on spaces $(\Xx_s)_{s=1}^S$ and a lower-semicontinuous cost $c:\Xx_1\times\cdots\times\Xx_S\to\RR\cup\{+\infty\}$ bounded from below, the problem reads

where $\Couplings(\alpha_1,\ldots,\alpha_S)$ is the set of probability measures whose $s$-th marginal is $\alpha_s$. This is still a linear program in the discrete setting, but its ambient tensor has size $\prod_s n_s$.

Monge Structure and Splitting-Set Twist¶

As in the two-marginal case, one would like to know when the optimal joint law is induced by deterministic maps from one marginal. The relevant non-degeneracy assumption is stronger than pairwise twist, because the other $S-1$ variables have to be recovered simultaneously. The condition below is the standard multi-marginal analogue used in the Monge-structure theory of Gangbo--Swiech and Pass Gangbo & Swiech, 1998Pass, 2011Pass, 2012Pass, 2015.

Coulomb Cost and Density-Functional Theory¶

A second canonical example, besides barycenters, comes from electronic structure. For $N$ electrons in $\RR^3$, the repulsive Coulomb interaction is the multi-body cost

with the value $+\infty$ on the collision set. Proposition Proposition: Multi-Marginal Monge Structure therefore does not apply verbatim: the Coulomb cost is neither finite nor differentiable on the whole product space. Any finite-energy plan gives zero mass to exact collisions, so the cost is smooth at almost every point charged by the plan, but this removes only the singularity; one must still establish differentiability of the dual potential and twist on the relevant splitting sets. Away from collisions,

For $N=2$, the map from $x_2$ to this vector is injective, so the ordinary two-marginal twist argument can be recovered under the required existence, duality and differentiability hypotheses. For $N\geq3$, however, the displayed total force does not by itself determine the entire tuple $(x_2,\ldots,x_N)$; twist on splitting sets, and hence a Monge representation, is not automatic. The previous proposition thus supplies a mechanism to verify in special Coulomb models, not a general existence theorem for co-motion maps.

If $\rho$ is an electron density with $\int_{\RR^3}\rho(x)\d x=N$ and $\al=\rho/N$ is the associated probability density, the strictly-correlated-electrons relaxation of density-functional theory is the equal-marginal problem

Since the cost and constraints are permutation invariant, symmetrizing any admissible plan does not change its value, so one may equivalently minimize over symmetric plans. This functional gives the smallest possible electron--electron repulsion compatible with the prescribed one-particle density; it appears as the strong-interaction limit in density-functional theory and was connected to optimal transport in Gori-Giorgi et al., 2009Buttazzo et al., 2012Cotar et al., 2013Di Marino et al., 2015. The deterministic ansatz writes a plan through co-motion maps

so that the position of one electron determines the positions of the others. The following cyclic version is the most common structural form of this ansatz in the strictly-correlated-electrons literature.

The Monge-structure proposition above explains the general mechanism that can force graph solutions, while the cyclic co-motion proposition records the additional equal-marginal symmetry used by co-motion maps. For the Coulomb cost, however, the singular repulsion and permutation symmetry make the structure delicate: co-motion maps are optimal in special geometries, but they are not universally optimal, and counterexamples are known Colombo & Stra, 2015Bindini et al., 2020. Thus the DFT problem is both a central application and a warning that multi-marginal OT is richer than a naive deterministic matching problem.

Figure Div shows the same phenomenon in a deliberately small one-dimensional model.

show_book_figure("multimarginal-coulomb-sinkhorn", width=760)

Entropic three-marginal Coulomb transport in one dimension. The three marginals are equal and the pairwise cost is a softened Coulomb repulsion. Each panel shows the $(X_1,X_2)$ marginal of the tensor Sinkhorn solution: small regularization pushes mass away from the collision diagonal, while larger regularization blurs the repulsive structure toward the independent reference.

Multi-Marginal Formulation of Barycenters¶

Wasserstein barycenters are the central example. For the squared Euclidean cost, one can introduce a latent barycenter point and eliminate it explicitly, leading to the multi-marginal cost

This linear support bound, rather than the cardinality of the full product grid, is the standard sparsity estimate for discrete Wasserstein barycenters Anderes et al., 2016.

Entropic Regularization of Multi-Marginal OT¶

As in the two-marginal case, adding an entropic penalty with respect to the product measure $\alpha_1\otimes\cdots\otimes\alpha_S$ leads to scaling algorithms:

The optimizer has the generalized Gibbs form

and generalized Sinkhorn iterations alternately update one potential $f_s$ so that the $s$-th marginal is correct. This formula is direct but, without additional structure, it is mostly a conceptual baseline. In the discrete case, even storing the Gibbs tensor or the coupling requires $\prod_s n_s$ entries, and the unstructured multi-marginal Sinkhorn complexity inherits this exponential dependence on the number of marginals Lin et al., 2019.

Treewidth and Graphical Structure¶

The important exception is when the cost factors over a sparse interaction graph. Let $G=(V,E)$ be a finite undirected graph with $V=\{1,\ldots,S\}$ and suppose, for simplicity, that

The relevant complexity parameter is not merely the number of edges, but the largest intermediate interaction created when variables are summed out.

The last condition is the running-intersection property. A tree decomposition equipped with factors assigned to its bags is also called a junction tree. Equivalently, choose an order in which to eliminate the vertices of $G$. Just before eliminating a vertex, connect all of its remaining neighbors, thereby adding fill-in edges. The induced width of the order is the largest number of remaining neighbors encountered. Treewidth is the minimum induced width over all elimination orders, and the corresponding bags consist of each eliminated vertex together with those neighbors. This equivalent viewpoint explains why treewidth controls exact summation.

The same construction applies to higher-order factors $c_A(x_A)$ indexed by subsets $A\subseteq V$: form the primal interaction graph by connecting every pair of variables occurring in a common factor, then compute the treewidth of that graph.

Junction-Tree Contractions Inside Sinkhorn¶

The treewidth reduction replaces each full tensor contraction in Sinkhorn by exact sum-product messages. In the discrete setting, define the edge Gibbs matrices and current unary factors

Under (86), the current scaled coupling is represented implicitly as

When $G$ is a tree, let $N_G(r)$ denote the neighbors of $r$. The directed sum-product messages satisfy

Cutting the edge $(r,s)$ separates the tree into two components: $m_{r\to s}(i_s)$ is the total contribution of the component containing $r$, conditional on the boundary state $i_s$. Conditioning first on $i_r$ gives the message recursion. A leaf-to-root pass followed by a root-to-leaf pass computes every directed message and hence every current marginal,

For the block $s$ currently selected by cyclic scaling, the exact coordinate update is

with componentwise products and quotients. This changes only the unary factor $h_s$. Updating all blocks at once would instead define a Jacobi scheme. Thus the expensive denominator of one generalized Sinkhorn block update is evaluated by messages rather than by enumerating all multi-indices.

For a general tree decomposition, choose one host bag for each unary factor, assign each edge factor in (88) to a bag containing both endpoints, and denote the product assigned to bag $B_q$ by $\psi_q(i_{B_q})$. For adjacent bags $q,q'\in Q$, set $\mathcal S_{q,q'}=B_q\cap B_{q'}$. Junction-tree messages take the form

After a collect-distribute pass, the calibrated bag belief is

For any $s\in B_q$, summing $\mathfrak b_q$ over $B_q\setminus\{s\}$ gives $\widehat a_s$; the running-intersection property ensures that every bag containing $s$ gives the same result.

Redundant adjacent bags may be contracted, so one can assume that the decomposition is reduced. If its width is $w$, every bag contains at most $w+1$ variables and every separator contains at most $w$. Consequently, if $n_s\leq n$, one collect-distribute inference pass for fixed scalings costs

More precisely, the time is bounded by