API Reference

GraphTransportation

Julia package for discrete transport geometry on graphs, implementing the framework of Erbar, Rumpf, Schmitzer, and Simon — Computation of optimal transport on discrete metric measure spaces.

The package provides:

• Geodesics (discrete_transport, transport_cost): compute discrete transport geodesics between probability measures on a graph via a Galerkin-discretized Chambolle-Pock primal-dual algorithm.
• Barycenters (barycenter): compute discrete transport Fréchet means via gradient descent.
• Coordinate recovery (analysis): recover barycentric coordinates of a measure with respect to a reference family by solving a quadratic programme on the Gram matrix of logarithmic maps.
• Entropic optimal transport barycenters (sinkhorn_barycenter, simplex_regression): entropic optimal transport barycenter computation and coordinate recovery via the Sinkhorn algorithm.

geodesic(G::MarkovGraph, ρA, ρB; method=:socp, kwargs...) -> GeodesicSolution

The discrete transport geodesic between densities ρA and ρB on G, by one of:

• :socp (default): a single second-order-cone program (geodesic_socp). Handles any densities, including boundary-supported ones; time-discretization error O(1/N). Keywords: N, optimizer, silent.
• :shooting: Newton shooting on the Hamiltonian flow (log_map), then the flow is integrated to produce the path. Exact in time, fastest, but requires strictly positive endpoints (PositivityFloorError/assertion otherwise). Keywords: nsteps, tol, maxiters, φ0_init.
• :chambolle_pock: the paper's Galerkin-discretized primal-dual iteration (discrete_transport). Reference implementation; slowest. Keywords: N, tol, maxiters, σ, τ.
• :sinkhorn: the entropic displacement interpolation for a ground cost: the path is the entropic barycenter of the two endpoints at weights (1−t, t) for t = 0, 1/N, …, 1 (the two-reference special case of barycenter(...; method=:sinkhorn)), a different object from the discrete transport geodesic. Requires cost (see ground_cost) and epsilon; keywords N (default 10) and iters (Sinkhorn budget, default 256). W2 is the entropic transport cost ⟨cost, P⟩ of the plan between the endpoints (no entropy term). The path's end columns are the blurred endpoints the Sinkhorn barycenter returns, not ρA/ρB exactly; m, φ0, φ1 are NaN-filled.

The transport metric's mean is the graph's G.mean (see MarkovGraph): :socp and :shooting honor every AdmissibleMean (the SOCP needs QuadLogMean for the logarithmic mean), :chambolle_pock supports only GeometricMean() and errors otherwise, and :sinkhorn's geometry is its ground cost, so it ignores the mean.

All methods return a GeodesicSolution. Its W2 is the squared distance; ρ is the n × (steps+1) density path and m the |E| × steps momentum path. The endpoint potentials φ0, φ1 are the gradients of W2 with respect to each endpoint (the SOCP's duals); for :shooting they are derived from the flow's velocity potentials (φ0 = -2φ(0), φ1 = 2φ(1)), and for :chambolle_pock they are not available (NaN-filled). Caveat for :chambolle_pock: its momentum path agrees with the other methods away from the endpoints (to O(1/N)), but the first-interval momentum m0 is a boundary-cell artifact of the Galerkin discretization, off by roughly 50% and not converging with N; it is what the paper's momentum-based analysis uses, and is self-consistent there, but it is not the initial momentum in the sense the other two methods return. status is the JuMP termination status for :socp and :converged for the other methods; solvetime is wall-clock seconds.

transport_cost(Q, μ, ν; N=64, σ=0.5, τ=0.5, maxiters=2^16, tol=1e-10,
               progress=false) -> Float64

Return the discrete transport distance W(μ, ν) on the graph defined by Q, computed as √(action(discrete_transport(Q, μ, ν; ...))).

Arguments

• Q: row-stochastic Markov transition matrix defining the graph
• μ, ν: probability densities w.r.t. the stationary distribution of Q
• N: number of time steps for the geodesic discretization (default 64)
• σ, τ: Chambolle-Pock step sizes; must satisfy σ·τ < 1 (default 0.5)
• maxiters: maximum Chambolle-Pock iterations (default 65536)
• tol: convergence tolerance on the density change between iterates
• progress: show a progress spinner if true

transport_cost(G::MarkovGraph, ρA, ρB; method=:socp, kwargs...) -> Float64

The discrete transport distance 𝒲(ρA, ρB) (not squared): sqrt(geodesic(...).W2). See geodesic for the methods and keywords.

barycenter(M, weights, Q; h=1.0, maxiters=100, tol=1e-8,
           geodesic_tol=1e-10, geodesic_steps=100,
           return_stats=false, initialization_index=1,
           initialization=nothing, uniform_init=false,
           geodesic_warmstart=true, verbose=true,
           stagnation_window=100, stagnation_tol=0.01,
           h_min=1e-8) -> ν  (or (ν, norm_diffs, variances) if return_stats=true)

Compute the discrete transport barycenter of the reference measures in M with respect to weights via gradient descent on the graph defined by Q.

Convergence is assessed on the gradient norm ‖∇J‖; a stagnation detector halves the step size when relative improvement over the last stagnation_window iterates falls below stagnation_tol.

Arguments

• M: V × p matrix whose columns are reference probability densities
• weights: non-negative length-p vector summing to 1
• Q: row-stochastic Markov transition matrix defining the graph
• h: initial WGD step size (default 1.0)
• maxiters: maximum number of WGD iterations
• tol: gradient-norm convergence threshold
• geodesic_tol: convergence tolerance passed to each discrete_transport call
• geodesic_steps: number of time steps for each geodesic
• return_stats: if true, return (ν, norm_diffs, variances) instead of ν
• initialization: explicit starting measure; takes precedence over index-based init
• initialization_index: column of M used as the starting point (default 1)
• uniform_init: if true, start from the uniform density (density ≡ 1)
• geodesic_warmstart: reuse geodesics from the previous iteration as warm starts
• verbose: show a ProgressMeter display with per-iteration statistics
• stagnation_window: number of iterations over which to measure relative improvement
• stagnation_tol: minimum relative improvement before halving h
• h_min: minimum allowed step size; stagnation halving stops here

barycenter(G::MarkovGraph, refs, λ; method=:socp, kwargs...) -> (ν, J, info)

The discrete transport barycenter of the reference densities refs (a vector of densities on G) with weights λ, i.e. the minimizer of J(ν) = Σᵢ λᵢ 𝒲²(refsᵢ, ν), by:

• :socp (default): one joint second-order-cone program (barycenter_socp), solved to its global optimum. info = (; geodesics) holds one GeodesicSolution per reference with λᵢ > 0. Keywords: N, optimizer, silent.
• :shooting: intrinsic (Riemannian) gradient descent with exact-in-time geodesics: each iteration log-maps ν to every reference (log_map, warm-started from the previous iteration), forms the descent direction Σᵢ λᵢ φ0ᵢ as a potential, and moves along the geodesic with exp_map. Requires strictly positive references and works with every AdmissibleMean, including the exact LogarithmicMean, which the SOCP cannot represent. info = (; iters, status, J_hist, grad_hist, h): the objective and the Riemannian gradient norm per iteration (including the final point), the final step size, and status — :converged (gradient norm below tol), :stalled (no step decreases the objective by more than ftol relative: the descent has reached the precision of the log maps), or :maxiters (warns). Keywords: h (step, default 1, halved on a positivity-floor hit or an objective increase, doubled back toward its initial value after an accepted step), maxiters (200), tol (gradient-norm threshold, 1e-5), ftol (relative objective decrease below which the descent stops, 1e-12; on the 49-node USA graph the objective's noise floor is about 1e-13 relative and the gradient's about 1e-7), log_tol (residual tolerance of the log maps, 1e-12; the objective's noise is roughly the square of the potentials' error), nsteps (integrator steps, 150), init (starting density, default the λ-weighted average of the references), verbose. Steps are accepted only on a strict objective decrease, so J_hist is monotone; a step's decrease is about h‖g‖², so asking for tol below √(ftol·J) typically ends as :stalled rather than :converged. First-order, so it converges linearly; the SOCP remains the certificate.
• :chambolle_pock: the paper's intrinsic gradient descent with Chambolle-Pock geodesics (the barycenter(M, weights, Q) method). info = (; norm_diffs, variances) are the descent's per-iteration statistics. Keywords: h, maxiters, tol, geodesic_tol, geodesic_steps, verbose, and the rest of that method's options. J is evaluated afterwards with Chambolle-Pock geodesics at the same settings.
• :sinkhorn: the entropically regularized Wasserstein barycenter for a ground cost (Benamou et al. 2015; Bonneel, Peyré & Cuturi 2016; sinkhorn_barycenter). This is a different object from the discrete transport barycenter of the other methods: it depends on the choice of cost (see ground_cost) and on epsilon, and even for a single reference it returns a blurred copy of that reference. Both cost and epsilon are required keywords; iters (default 256) is the Sinkhorn budget. J = Σᵢ λᵢ ⟨cost, Pᵢ⟩ over the entropic plans Pᵢ (no entropy term) and is not comparable with the J of the other methods; info = (; cost, epsilon, iters, marginal_errors). refs are densities with respect to π as elsewhere (converted to probability vectors internally).

Returns the barycenter ν, the objective value J at ν, and the method-specific info named tuple.

analysis(ν, M, Q; N=100, tol=1e-10, compute_condition=false,
         return_system=false) -> weights

Recover the barycentric coordinates of ν with respect to the reference measures in M (a V × p matrix whose columns are probability densities w.r.t. the stationary distribution of Q).

The method proceeds in three steps:

1. For each reference measure M[:,i], compute the geodesic from ν to M[:,i] via discrete_transport and extract the initial tangent vector m[1,:,:] (the logarithmic map at ν).
2. Assemble the p × p Gram matrix A[i,j] = ⟨log_ν(M[:,i]), log_ν(M[:,j])⟩_{g(ν)} under the metric tensor g(ν).
3. Solve the simplex-constrained quadratic programme min_{w≥0, Σw=1} w'Aw exactly (solve_simplex_qp), which recovers the weights at the discrete transport barycenter.

Optional keyword arguments:

• N: number of time steps for each geodesic computation (default 100).
• tol: convergence tolerance for each discrete_transport call.
• compute_condition: if true, prints the condition number of A.
• return_system: if true, returns (weights, A) instead of weights alone.

analysis(G::MarkovGraph, target, refs; method=:socp, kwargs...) -> λ̂ (or (λ̂, A))

Recover the barycentric coordinates of target with respect to the reference densities refs: compute the geodesic from target to each reference, build the Gram matrix of their tangent vectors at target, and solve min_{λ∈Δ} λᵀAλ. The geodesic solver is:

• :socp (default): analyze_socp — endpoint potentials from the geodesic SOCP. Keywords: N, optimizer, convention.
• :shooting: analyze_shooting — potentials from log_map; interior data only. Keywords: nsteps, tol, φ0_inits.
• :chambolle_pock: the paper's analysis(ν, M, Q) — initial momenta from discrete_transport. Keywords: N, tol.
• :sinkhorn: Wasserstein barycentric coordinates for a ground cost (Bonneel, Peyré & Cuturi 2016; simplex_regression): L-BFGS over the simplex on ½‖P(λ) − target‖², with P(λ) the entropic barycenter of refs, differentiated through the Sinkhorn iterations. Not a Gram-matrix method, so return_system and compute_condition are not supported (an error, not silently ignored). Requires cost and epsilon; keywords iters (Sinkhorn budget, default 256) and α0 (initial pre-softmax point, default 0, i.e. uniform weights). The result depends on cost and epsilon, and recovers exactly the weights of a barycenter synthesized by barycenter(...; method=:sinkhorn) with the same cost, epsilon and iters.

Common keywords (Gram-matrix methods): compute_condition, return_system (also return the Gram matrix A).

A barycenter is recovered to solver tolerance only by the method (and time resolution) that synthesized it; other methods recover it to their discretization error, because each checks stationarity in its own discrete convention.

GeodesicSolution

Result of solving the geodesic SOCP (geodesic_socp). Note the naming convention: W2 is the squared discrete transport distance; the metric itself is sqrt(W2).

Fields:

• W2::Float64: squared discrete transport distance, i.e. the SOCP objective value.
• ρ::Matrix{Float64}: the n × (N+1) density path, ρ[:, i+1] at time i/N.
• m::Matrix{Float64}: the |E| × N momentum path.
• m0::Vector{Float64}: m[:, 1], the initial momentum (convenience).
• φ0::Vector{Float64}, φ1::Vector{Float64}: the discrete potentials at the two endpoints, read off the duals of the endpoint constraints ρ[:,1] == ρA and ρ[:,N+1] == ρB and divided by π, so that they are exactly the gradient of W2 with respect to each endpoint density under the π-weighted pairing: δW2 = ⟨φ0, δρA⟩_π + ⟨φ1, δρB⟩_π. This is the discrete Hamilton–Jacobi potential (m ≈ -½ θ ∇φ in the continuum limit; see _geodesic_block!), and is the quantity the potential-based analysis Gram matrix (analyze_socp) is built from. Each is defined only up to an additive constant; JuMP/Clarabel return the representative with ⟨φ, 1⟩ fixed by the solver, which is harmless since only ∇φ is ever used.
• status: JuMP termination status of the solve.
• solvetime::Float64: wall-clock solve time in seconds.

MarkovGraph

Graph primitive for the SOCP formulation.

Stores an undirected graph together with a reversible Markov chain on it, using a per-edge (rather than per-node-pair-matrix) representation: momenta, potential gradients, and other edge quantities are each a single Vector{Float64} of length |E|, with an arbitrary but fixed orientation per edge. Antisymmetry of edge fields (Lemma 2.4 of Erbar et al.) is therefore hard-coded into the representation.

Fields:

• n::Int: number of vertices.
• E::Vector{Tuple{Int,Int}}: oriented edge list, E[e] = (x, y) meaning the edge is oriented from x to y.
• π::Vector{Float64}: stationary distribution.
• Q::SparseMatrixCSC{Float64,Int}: transition rate matrix.
• κ::Vector{Float64}: edge weights κ[e] = Q[x,y] * π[x] for e = (x,y), which by reversibility equals Q[y,x] * π[y].
• D::SparseMatrixCSC{Float64,Int}: n × |E| incidence matrix with D*m == graph_divergence(G, m) for any edge field m, cached at construction so repeated divergences (e.g. one per SOCP time step) don't rebuild it. D[x,e] = -Q[x,y], D[y,e] = Q[y,x] for edge e = (x,y).
• mean::AdmissibleMean: the mobility θ(s, t) of the transport metric on this graph (GeometricMean() by default; see AdmissibleMean). The mean is part of the geometry, like the graph itself, so it lives here rather than as a per-call option: every geodesic, barycenter and analysis computed from the same MarkovGraph uses the same metric, and a target synthesized on one graph is analysed with the mean it was made with. Chambolle-Pock supports only the geometric mean; the SOCP needs QuadLogMean rather than LogarithmicMean; the shooting maps accept all of them.

geodesic_socp(G::MarkovGraph, ρA, ρB; N=10, optimizer=Clarabel.Optimizer, silent=true, check=true) -> GeodesicSolution

Compute the discrete transport geodesic between densities ρA and ρB on G as a single second-order-cone program, rather than via the Chambolle-Pock primal-dual iteration (discrete_transport). Returns the squared distance W2 = ‖ρA - ρB‖_𝒲²; the metric distance is sqrt(W2).

The mobility θ is the graph's G.mean: GeometricMean() (default), ArithmeticMean(), HarmonicMean() or QuadLogMean(K) (the logarithmic mean by Gauss–Legendre quadrature, K power cones per edge and time step; a graph built with LogarithmicMean() has no conic form and errors here, use MarkovGraph(G; mean=QuadLogMean(8))). See _mean_cone! for the conic representations and AdmissibleMean for the theory.

N is the number of time-discretization intervals (h = 1/N); the returned ρ has N+1 columns and m has N columns. optimizer is any solver JuMP can dispatch to that supports rotated second-order cone constraints (Clarabel by default). With check=true (default) a solve that does not reach OPTIMAL/ALMOST_OPTIMAL throws rather than returning the solver's last iterate (see _check_solved).

barycenter_socp(G, refs, λ; N=10, optimizer=Clarabel.Optimizer, silent=true, check=true)
    -> (ν::Vector{Float64}, J::Float64, geodesics::Vector{GeodesicSolution})

Compute the discrete transport barycenter of reference measures refs (a vector of probability densities on G, one per reference) with weights λ as a single joint second-order-cone program, rather than by intrinsic gradient descent (barycenter). Replicates the geodesic SOCP block (_geodesic_block!) once per reference with λ[i] > 0, sharing a single free right endpoint ν across all of them.

Returns:

• ν: the barycenter, i.e. the shared optimal right endpoint.
• J: the optimal objective value, Σᵢ λᵢ 𝒲_h²(refs[i], ν).
• geodesics: one GeodesicSolution per reference with λ[i] > 0 (in the same relative order as refs), giving the geodesic from that reference to ν. Each carries its endpoint potentials (φ0 at the reference, φ1 at ν, both for the unweighted geodesic — the λ[i] factor is divided out of the duals). KKT stationarity of the joint program with respect to ν is exactly Σᵢ λᵢ φ1ᵢ = const on the support of ν, which is the stationarity condition the potential-based analysis (analyze_socp) checks.

References with λ[i] == 0 are dropped entirely rather than solved with a zero weight. The mobility is G.mean, and check behaves, as in geodesic_socp. λ must be a probability vector (λ .>= 0, sum(λ) ≈ 1).

analyze_socp(G, target, refs; N=10, optimizer=Clarabel.Optimizer,
             convention=:potential, compute_condition=false, return_system=false)
    -> λ̂ (or (λ̂, A))

The :socp analysis backend: recover the barycentric coordinates of target with respect to refs by solving, for each reference, the geodesic SOCP from target to it (geodesic_socp), building a p × p Gram matrix of the resulting tangent vectors at target, and solving the simplex QP min_{λ∈Δ} λᵀAλ shared with the existing analysis (solve_barycentric_coordinates_qp).

convention selects which tangent-vector representation the Gram matrix is built from:

• :potential (default): the Riemannian Gram matrix at target, A_ij = Σ_e κ_e θ(target)_e (∇φ_i)_e (∇φ_j)_e, with φ_i the endpoint potential GeodesicSolution.φ0 — the exact gradient of 𝒲_h²(target, ref_i) with respect to target. Since a barycenter_socp solution is, by its KKT conditions, exactly a point where Σᵢ λᵢ φ_i = const on its support, this recovers the synthesis weights of an SOCP barycenter to solver tolerance at any N, including N=2.
• :momentum: the earlier implementation, which used the initial momentum m0 as the tangent vector with the dense metric_tensor(target) weighting, exactly as the Chambolle-Pock-based analysis does. m0 lives on the first time interval (m0 = -(1/2h) θ(ρ̄_{1/2}) ∇ψ_{1/2}, with the midpoint density and the half-step potential), so it is only an O(h) proxy for the endpoint potential, and Σᵢ λᵢ m0ᵢ ≈ 0 is not the stationarity condition the SOCP barycenter satisfies. At coarse N this mismatch can be large (tens of percent recovery error at N=2 on irregular graphs). Kept for comparison experiments only.

The mobility for both the geodesic solves and the Gram-matrix weighting is G.mean; recovery is exact for a target synthesized on the same graph (same mean).

refs is a vector of reference probability densities on G. Returns the recovered weight vector λ̂ (exact simplex QP on the small p × p Gram matrix), or (λ̂, A) if return_system=true.

hamiltonian(G::MarkovGraph, ρ, φ) -> Float64

H(ρ,φ) = ½ Σ_e θ(ρ_x,ρ_y) (∇φ)_e² κ_e. Requires ρ strictly positive (see ρ_floor).

hamiltonian_flow(G::MarkovGraph, ρ, φ) -> (ρ̇, φ̇)

Equations of motion:

ρ̇(x) = Σ_y θ(ρ(x),ρ(y)) (φ(x)-φ(y)) Q(x,y)  =  -div(θ(ρ)∘∇φ)(x)
φ̇(x) = -½ Σ_y ∂₁θ(ρ(x),ρ(y)) (φ(x)-φ(y))² Q(x,y)

with θ = G.mean and ∂₁θ = partial_s(G.mean, ·, ·) (for the geometric mean ½√(t/s), which divides by ρ(x); hence the positivity floor). Requires ρ strictly positive.

integrate_hamiltonian(G::MarkovGraph, ρ0, φ0; nsteps=150, T=1.0, floor_rtol=1e-6)
    -> (ρ_path, φ_path)

Integrate the Hamiltonian flow forward from (ρ0, φ0) over [0, T] using fixed-step classical RK4 rather than a symplectic integrator: over a unit time interval at 100-200 steps, symplecticity is a nicety, not a requirement. Returns the full paths as n × (nsteps+1) matrices.

Positivity guard: if any density in the proposed next step would fall below ρ_floor(G; rtol=floor_rtol), the step is halved (up to 4 times) before giving up and throwing PositivityFloorError, which callers should catch and act on (fall back to geodesic_socp or mollify).

ρ_floor(G::MarkovGraph; rtol=1e-6) -> Float64

Minimum density the shooting ODE machinery will tolerate: rtol relative to min(π) (every entry point asserts minimum(ρ) > ρ_floor). Below this, θ_geo and its partial derivative ∂₁θ_geo(s,t) = ½√(t/s) are numerically unsafe (division by a near-zero s), and the caller should fall back to geodesic_socp or a mollified approximation (log_map_mollified) instead.

PositivityFloorError <: Exception

Thrown by integrate_hamiltonian when a trajectory would drive some density below ρ_floor even after repeated step bisection. log_map catches exactly this type (and nothing broader) to damp its initial guess and shorten line-search steps; any other exception from the flow propagates.

weighted_laplacian(G::MarkovGraph, ν) -> SparseMatrixCSC

The n × n weighted graph Laplacian L_θ(ν) = ∇ᵀ Diag(κ ∘ θ(ν)) ∇, i.e. (L φ)(x) = Σ_y κ_{xy} θ(ν_x, ν_y) (φ(x) − φ(y)). Symmetric positive semidefinite with kernel spanned by the constants, and related to the Hamiltonian flow by π ∘ ρ̇ = L_θ(ν) φ (equivalently ρ̇ = −div(θ(ν)∘∇φ), see hamiltonian_flow).

solve_weighted_laplacian(G::MarkovGraph, ν, b) -> φ

Solve L_θ(ν) φ = b for the unique solution in the gauge ⟨φ, 1⟩_π = 0. b must be orthogonal to the constants (sum(b) ≈ 0), which every right-hand side arising here is (π ∘ (target − ν) for two probability densities, or ∇ᵀ(κ ∘ m) for any momentum). Implemented via the rank-one regularization (L + π πᵀ) φ = b: for such b this has the same solution as L φ = b and enforces the gauge automatically.

momentum_to_potential(G::MarkovGraph, ν, m) -> φ

Recover the potential φ with m = θ(ν) ∘ ∇φ (gauge ⟨φ, 1⟩_π = 0) from a momentum edge vector m, by solving L_θ(ν) φ = ∇ᵀ(κ ∘ m). If m is not exactly a gradient field this returns the θ(ν)-weighted least-squares projection, i.e. the potential of the gradient part of m in the Hodge sense.

exp_map(G::MarkovGraph, ν, tangent; t=1.0, nsteps=150, kind=:auto, floor_rtol=1e-6) -> ρ_end

The Riemannian exponential map at ν. Integrates the Hamiltonian flow (integrate_hamiltonian) from (ν, φ0) for time t and returns the endpoint density.

tangent is either a potential φ0 ∈ Rⁿ (kind=:potential) or a momentum m0 ∈ R^{|E|} (kind=:momentum), in which case φ0 is recovered first with momentum_to_potential. With kind=:auto (default) the kind is inferred from the length, which is ambiguous when n == |E| (e.g. a cycle); then kind must be given. The potential is gauge-fixed to ⟨φ0, 1⟩_π = 0 before integrating (this does not change the flow's ρ).

Requires ν strictly positive (see ρ_floor); errors, rather than returning garbage, if the flow hits the positivity floor before time t.

log_map(G::MarkovGraph, ν, target; φ0_init=nothing, tol=1e-9, maxiters=50, nsteps=150,
        floor_rtol=1e-6, verbose=false) -> (; φ0, m0, W2, iters, residual)

The Riemannian logarithm at ν, by single shooting. Solves F(φ0) := ρ(1; ν, φ0) − target = 0 with a damped Newton iteration over the mean-zero potentials (n − 1 unknowns, the gauge ⟨φ0, 1⟩_π = 0 and the mass constraint each removing one dimension), with the Jacobian computed by ForwardDiff through integrate_hamiltonian and a backtracking line search on ‖F‖_π.

Initialization is the linearized geodesic L_θ(ν) φ0 = π ∘ (target − ν), exact to first order in target − ν, unless φ0_init (a potential) is given, e.g. from a previous solve at a nearby base point. This is the warm-start mechanism; no (base, target)-keyed cache is kept, callers hold their own φ0. Cold starts on interior data take 2–4 Newton iterations; the initial guess is damped by halving if it overshoots the positivity floor (far-apart concentrated endpoints, ~9 iterations).

Near the positivity floor the integrator's step bisection makes the residual piecewise-smooth in φ0 (ForwardDiff differentiates whichever branch the trajectory took), so Newton can stall at the scale of those jumps; this is the failure mode for strongly mollified boundary data (ε ≲ 1e-4) and is reported as a line-search failure. Multiple shooting is not implemented.

Returns the potential φ0, the momentum m0 = θ(ν) ∘ ∇φ0, the squared distance W2 = 2H(ν, φ0), the Newton iteration count, and the final residual ‖F‖_π. Errors if the iteration has not reached tol after maxiters steps, or if the shooting trajectory persistently hits the positivity floor (fall back to geodesic_socp, or mollify with log_map_mollified). Requires ν and target strictly positive (see ρ_floor).

log_map_mollified(G::MarkovGraph, ν, target; εs=(1e-2, 1e-3, 1e-4), tol=1e-7, kwargs...)
    -> (; W2, W, φ0, m0, fit, εs, Ws, approximate=true)

log_map for endpoints with zero or near-zero entries, where shooting cannot be run directly. Both endpoints are mollified toward the uniform density, ρ_ε = (1−ε)ρ + ε·𝟙 (still a probability density), log_map is run for each ε in εs (warm-starting each from the previous, coarser ε), and the distance is extrapolated to ε → 0 by a least-squares fit W(ε) ≈ W₀ + a·√ε (the mollification perturbs the endpoints by O(√ε) in 𝒲). Returns the extrapolated W and W2 = W², plus the potential/momentum from the smallest solved ε (at the mollified base point ν_ε, so only approximately a tangent at ν), the fit coefficients (W₀, a), and the raw per-ε distances. Levels whose shooting fails (very small ε makes the near-boundary geodesic stiff) are skipped with a warning as long as two remain. Results are flagged approximate=true; for boundary-supported data geodesic_socp needs no such fallback and is the reference. Empirically the raw W at the smallest ε is often already as accurate as the extrapolation: on 5x5-grid probes the mollification error decayed faster than √ε (successive differences shrank by ~7x per decade of ε, vs. ~3.2x predicted), and the fit over-corrected by 0.3–0.5% while the raw smallest-ε value was within 0.05–0.15% of the SOCP reference. The √ε model is kept as specified; both values are returned so the numerics can be compared with the theory. The shooting tolerance tol defaults to a looser 1e-7 here (the result is approximate anyway, and the near-boundary trajectories make the last digits of the residual hard to reach). Remaining kwargs go to log_map.

analyze_shooting(G::MarkovGraph, target, refs; nsteps=150, tol=1e-9, φ0_inits=nothing,
                 compute_condition=false, return_system=false) -> λ̂ (or (λ̂, A))

The :shooting analysis backend: like analyze_socp, but each reference's potential comes from log_map(G, target, ref) (the Hamiltonian velocity potential φ0 at target) instead of the geodesic SOCP's endpoint dual. The Gram matrix and simplex QP are shared with analyze_socp (internal helper potential_gram_qp). Requires strictly positive target and refs; a failure on any one reference propagates (there is no per-reference fallback to the SOCP), and the same weighted Laplacian is refactored once per reference, which is negligible next to the Newton solves.

φ0_inits, if given, is a vector of warm-start potentials, one per reference.

Note that this checks stationarity in the Hamiltonian flow's (RK4, nsteps) discretization, which differs from both barycenter_socp's and the Chambolle-Pock descent's. A barycenter synthesized by either of those is therefore recovered only to O(h) in that method's time step, not to solver tolerance — see the analyze_socp docstring for why synthesis and analysis conventions must match.

discrete_transport(Q, μ, ν; N=64, σ=0.5, τ=0.5, maxiters=2^16, tol=1e-10,
                   progress=false, initialization=nothing) -> ErbarBundle

Compute the discrete transport geodesic from μ to ν on the graph defined by the Markov transition matrix Q, using the Chambolle-Pock primal-dual algorithm of Erbar et al. 2020.

Returns an ErbarBundle encoding the full geodesic; its vector.m[1,:,:] field is the initial tangent vector (logarithmic map at μ), and action(result) gives the squared discrete transport distance.

Arguments

• Q: row-stochastic Markov transition matrix defining the graph
• μ, ν: probability densities w.r.t. the stationary distribution of Q
• N: number of time steps for the geodesic discretization (default 64)
• σ, τ: Chambolle-Pock step sizes; must satisfy σ·τ < 1 (default 0.5)
• maxiters: maximum Chambolle-Pock iterations (default 65536)
• tol: convergence tolerance on the density change between iterates
• progress: show a progress spinner if true
• initialization: an ErbarBundle from a prior solve to warm-start from

action(u::ErbarBundle) -> Float64

Compute the discrete action functional of the curve encoded in u, as defined in equation (18) of Erbar et al. 2020:

A(ρ, m) = (h/2) · Σ_t Σ_{x,y} (m[t,x,y]² / θ[t,x,y]) · Q[x,y] · π[x]

where h = 1/N is the time step. The integrand uses the convention m²/θ = 0 when m = θ = 0, and +Inf when m ≠ 0 and θ = 0.

The square root of the action is the discrete Wasserstein distance between the boundary measures u.cache.μ and u.cache.ν.

sinkhorn_barycenter(coords, measures, target, cost, epsilon; iters=256)

Compute the Wasserstein barycenter of measures with weights coords using the Sinkhorn algorithm with regularization epsilon.

simplex_regression(measures, target, cost, epsilon; iters=256, α0=zeros(S), optim_options=Optim.Options())

Wasserstein barycentric coordinates of target with respect to the columns of measures (Bonneel, Peyré & Cuturi 2016, §4.3): minimize E_L(λ) = ½‖P^(L)(λ) − target‖² over the simplex by L-BFGS on α with λ = softmax(α), using the gradient from sinkhorn_differentiate with the softmax Jacobian applied. α0 = 0 is the paper's λ0 = 1/S. iters is the Sinkhorn iteration budget used for both the objective and the gradient. Returns λ̂ ∈ Δ.

ground_cost(G::MarkovGraph, rule::Symbol; t=graph_diameter(G), normalize=true) -> Matrix{Float64}

A ground cost matrix on the nodes of G for the entropic (Sinkhorn) barycenter, barycenter(...; method=:sinkhorn):

• :shortest_path: squared BFS hop count (edges of G unweighted).
• :diffusion: squared diffusion distance at time t under the chain G.Q, D_t(x,y)² = Σ_z (Qᵗ[x,z] − Qᵗ[y,z])² / π(z). t defaults to the graph diameter so no pair has zero distance. This uses the chain G.Q itself, which coincides with the random walk on the adjacency only for unweighted chains.

With normalize=true (default) the matrix is divided by its maximum so entries lie in [0, 1], which keeps the kernel exp(-cost/epsilon) from underflowing at the usual epsilon. Throws if the graph is disconnected.

graph_diameter(G::MarkovGraph) -> Int

Largest BFS hop count between any two nodes, treating the edges of G as unweighted.

markov_chain_from_edge_list(E) -> (Q, π)

Construct the uniform random walk Markov chain on the graph defined by the edge list E. Each element of E is a tuple (i, j) representing an undirected edge. All edges are given equal weight (unweighted graph), so the transition probability from node x to each neighbor is 1 / degree(x). The stationary distribution satisfies π[x] ∝ degree(x).

Returns (Q, π) where Q is the row-stochastic transition matrix and π is the stationary distribution vector.

markov_chain_from_adjacency_matrix(A) -> (Q, π)

Construct the uniform random walk Markov chain from a binary adjacency matrix A. Equivalent to markov_chain_from_weight_matrix(A). The transition matrix is Q[x,y] = A[x,y] / degree(x) and the stationary distribution is π[x] ∝ degree(x).

Returns (Q, π).

markov_chain_from_weight_matrix(W) -> (Q, π)

Construct a weighted random walk Markov chain from a non-negative symmetric weight matrix W. The transition probability from node x to node y is Q[x,y] = W[x,y] / sum(W[x,:]), and the stationary distribution is π[x] = sum(W[x,:]) / sum(W).

Returns (Q, π).

stationary_from_transition(Q) -> π

Compute the stationary distribution of a row-stochastic transition matrix Q by solving the linear system (Q' - I)π = 0 with the normalization constraint sum(π) = 1. The overdetermined system is assembled and solved via least squares (\).

triangle_markov_chain() -> (Q, π)

Construct the uniform random walk Markov chain on the 3-cycle (triangle).

square_markov_chain() -> (Q, π)

Construct the uniform random walk Markov chain on the 4-cycle (square).

cube_markov_chain() -> (Q, π)

Construct the uniform random walk Markov chain on the 3-cube (8 nodes, 12 edges). Returns (Q, π) where Q is the transition matrix and π the stationary distribution.

hypercube_markov_chain() -> (Q, π)

Construct the uniform random walk Markov chain on the 4-cube (16 nodes, 32 edges). All edges have equal weight, so the stationary distribution is uniform.

weighted_hypercube_markov_chain() -> (Q, π)

Construct a randomly weighted random walk Markov chain on the 4-cube (16 nodes, 32 edges). Edge weights are drawn uniformly from {1, …, 10} and symmetrized, so the stationary distribution is non-uniform. Results vary across calls.

grid_markov_chain() -> (Q, π)

Construct the uniform random walk Markov chain on a 9-node wheel-like graph: an 8-cycle (nodes 1–8) with a central hub (node 9) connected to the even nodes.

grid_markov_chain(n) -> (Q, π)

Construct the uniform random walk Markov chain on the n × n grid graph (n² nodes with nearest-neighbor edges).

triangular_prism_markov_chain() -> (Q, π)

Construct the uniform random walk Markov chain on the triangular prism graph (6 nodes: two triangular faces connected by three edges).

triangle_with_tail_markov_chain() -> (Q, π)

Construct the uniform random walk Markov chain on a 3-cycle with one pendant edge (node 4 attached to node 3).

T_markov_chain() -> (Q, π)

Construct the uniform random walk Markov chain on the T-shaped graph (path 1–2–3 with a branch 2–4).

double_T_markov_chain() -> (Q, π)

Construct the uniform random walk Markov chain on the double-T graph (8 nodes: two T-shaped graphs — nodes 1–4 and nodes 5–8 — connected at their corresponding vertices: 1↔5, 2↔6, 3↔7, 4↔8).

ma_house_markov_chain() -> (Q, π)

Construct the uniform random walk Markov chain on the Massachusetts state house district adjacency graph. Nodes are legislative districts; edges connect districts that share a non-trivial boundary (point adjacencies are excluded).

Reads district geometries from the bundled shapefile at src/experiments/HOUSE2021/HOUSE2021_POLY.shp. Requires LibGEOS and Shapefile to be loaded.

graph_gradient(Q, f) -> Matrix

Compute the graph gradient ∇_G(f), where G is the graph induced by the transition matrix Q and f is a function defined on the nodes.

Entry (i, j) of the result is f[i] - f[j] when the edge (i,j) exists (i.e. Q[i,j] ≠ 0) and zero otherwise. The matrix is antisymmetric. Asserts that Q and f have compatible dimensions.

graph_gradient(G::MarkovGraph, φ::AbstractVector) -> Vector{Float64}

Compact-edge-vector counterpart of the dense graph_gradient: (∇φ)[e] = φ[x] - φ[y] for the oriented edge e = (x, y).

add_graph_gradient!(dest, Q, f)

Add the graph gradient ∇G(f) directly into dest in-place. Equivalent to `dest .+= graphgradient(Q, f)` but with no intermediate array allocation: the loop body writes only to edge positions and skips non-edges via the sparsity check on Q.

graph_divergence(Q, m) -> Vector

Compute the graph divergence of the V×V edge field m with respect to the graph induced by Q. Entry i of the result is

(div m)[i] = (1/2) · Σ_j Q[i,j] · (m[j,i] - m[i,j])

graph_divergence(G::MarkovGraph, m::AbstractVector)

Compact-edge-vector counterpart of the dense graph_divergence, specialized to a single scalar per undirected edge (rather than a full antisymmetric V×V matrix). Agrees with the dense version under m_dense[x,y] = m[e], m_dense[y,x] = -m[e], and therefore satisfies the same adjoint identity ⟨φ, div m⟩_π = -⟨∇φ, m⟩_Q (with ⟨m,w⟩_Q := Σ_e κ[e] m[e] w[e]).

Implemented as G.D * m (see the D field), so m may hold Float64s or, e.g., JuMP variables/affine expressions when building an SOCP model.

graph_divergence!(out, Q, m_t)

Compute the graph divergence of the V×V edge field m_t into the pre-allocated length-V vector out, with no allocations.

laplacian_from_transition(Q) -> Matrix

Return the graph Laplacian L = Q - D where D is the diagonal degree matrix D[i,i] = Σ_j Q[i,j].

metric_tensor(ρ, mean=geomean)

Return the V × V metric tensor matrix for the probability density ρ.

Entry (i, j) is mean(ρ[i], ρ[j]) for i ≠ j and zero on the diagonal. The default mean is the geometric mean, giving the Riemannian metric associated to the graph Wasserstein geometry of Erbar et al. 2020. Passing mean=logmean gives the metric for the logarithmic-mean variant.

metric_tensor(G::MarkovGraph, ρ::AbstractVector, mean=G.mean) -> Vector{Float64}

Compact-edge-vector counterpart of the dense metric_tensor: θ[e] = mean(ρ[x], ρ[y]) for the oriented edge e = (x, y), with mean defaulting to the graph's own G.mean (any callable (s, t) -> θ can be passed explicitly). Does not include the edge weight κ; the Riemannian inner product of two potential gradients at ρ is ⟨∇φ, ∇ψ⟩_ρ = Σ_e κ[e] θ[e] (∇φ)[e] (∇ψ)[e].

avg_operator(N)

Return the sparse time-averaging operator of size (N-1) × N.

Row i contains 0.5 at columns i and i+1, so multiplying by a length-N vector computes the average of each consecutive pair. In the Erbar et al. 2020 discretization, this maps the node density curve ρ ∈ V{n,h}^1 (N+1 time points) to its time-averaged version ρavg ∈ V{n,h}^0 (N time points) via `avgoperator(N+1) * ρ`.

finite_difference_operator(N)

Return the sparse forward finite-difference operator of size (N-1) × N.

Multiplying by a length-N vector f returns the vector of forward differences [f[2]-f[1], f[3]-f[2], …, f[N]-f[N-1]].

AdmissibleMean

Abstract supertype of the means θ(s, t) (mobilities) that define a discrete transport metric: ‖∇φ‖²_ρ = Σ_e κ_e θ(ρ_x, ρ_y) (∇φ)_e². An admissible mean is continuous, symmetric, positively 1-homogeneous, concave, positive on (0,∞)² and normalized so θ(s, s) = s (Maas 2011). Concavity is what makes the action m²/θ jointly convex, hence the geodesic/barycenter programs convex.

Every concrete mean θ is callable, θ(s, t), and has partial_s(θ, s, t) (the partial derivative in the first argument; partial_t follows by symmetry). Both are generic in the element type so ForwardDiff can differentiate through them.

For all s, t > 0: HarmonicMean() ≤ GeometricMean() ≤ LogarithmicMean() ≤ ArithmeticMean(). The arithmetic mean is the only one with θ(0, t) ≠ 0, so it is the only one under which mass can flow out of an empty node; see the package documentation on boundary behavior.

θ(s,t) = √(st). The package's default, following Erbar et al. (2020).

θ(s,t) = (s+t)/2. Does not vanish at an empty node (θ(0,t) = t/2).

θ(s,t) = 2st/(s+t). The smallest admissible mean of the four; bounded ∂₁θ near zero.

θ(s,t) = (s−t)/(ln s − ln t) (Λ). The mean for which the heat flow is the gradient flow of the relative entropy (Maas 2011). Evaluated by a series near s = t so that the value and its derivative are smooth there (needed for ForwardDiff through the shooting maps). Not finitely conic-representable; the SOCP uses QuadLogMean instead.

QuadLogMean(K::Int=8)

Gauss–Legendre approximation of the logarithmic mean, Λ_K(s,t) = Σ_k w_k s^{α_k} t^{1−α_k} ≈ ∫₀¹ s^α t^{1−α} dα = Λ(s,t), with K nodes. Λ_K is itself an admissible mean (symmetric, homogeneous, concave, positive), so a transport metric built on it is exact, not approximate; and each term is a power-cone constraint, which is how the SOCP represents the logarithmic mean. Worst relative error against Λ over density ratios up to 1000: 8.5e-4 (K=4), 6e-7 (K=6), 1.2e-10 (K=8), 2e-15 (K=12).

partial_s(θ::AdmissibleMean, s, t)

∂θ/∂s(s, t). By symmetry partial_t(θ, s, t) = partial_s(θ, t, s).

partial_t(θ::AdmissibleMean, s, t)

∂θ/∂t(s, t), i.e. partial_s(θ, t, s) by the symmetry of admissible means.

geomean(x, y) -> Real

Geometric mean of real numbers x and y. Returns -Inf if either argument is negative (convention used by the graph Wasserstein metric tensor).

logmean(s, t; tol=1e-5) -> Real

Logarithmic mean of real numbers s and t:

L(s, t) = (s - t) / (log s - log t)

Returns -Inf if either argument is negative, 0 if either is zero and the other non-negative, and (s + t) / 2 (first-order Taylor approximation) when |s - t| ≤ tol.

logmean_partial_s(s, t; tol=1e-5)

Partial derivative of the logarithmic mean with respect to s.

For |s - t| < tol (near the diagonal), returns the limiting value 0.5. Otherwise returns the exact derivative:

∂/∂s[(s - t)/(log s - log t)] = (-s + t + s·log(s/t)) / (s·(log s - log t)²)

logmean_partial_t(s, t; tol=1e-5)

Partial derivative of the logarithmic mean with respect to t.

For |s - t| < tol (near the diagonal), returns the limiting value 0.5. Otherwise returns the exact derivative:

∂/∂t[(s - t)/(log s - log t)] = (s - t - t·log(s/t)) / (t·(log s - log t)²)

ErbarVector{M, T}

An element of the Hilbert space H described in equation (27) of Erbar et al. 2020, representing a discretized curve (ρ, m) together with the auxiliary variables (θ, ρ_minus, ρ_plus, ρ_avg, q) used by the Chambolle-Pock solver.

Type parameters:

• M: concrete 2-D array type for node fields (ρ, ρ_avg, q); typically Matrix{Float64}
• T: concrete 3-D array type for edge fields (m, θ, ρ_minus, ρ_plus); typically Array{Float64,3}

Parameterising on the concrete types allows Julia to fully specialise all operations on ErbarVector fields in the Chambolle-Pock inner loop.

Fields

• ρ: density curve, shape (N+1) × V
• m: edge flux, shape N × V × V
• θ: dual edge variable, shape N × V × V
• ρ_minus, ρ_plus: edge density copies, shape N × V × V
• ρ_avg: time-averaged density, shape N × V
• q: dual node variable, shape N × V

ErbarCache

Immutable container for the graph OT problem metadata: everything that needs to be computed only once given (Q, μ, ν, N).

Constructors:

• ErbarCache(Q, π, μ, ν, N) — use a pre-computed stationary distribution π
• ErbarCache(Q, μ, ν, N) — computes π internally
• ErbarCache(existing, μ_new, ν_new) — cheap warm-start copy reusing the cached linear-system factorisations (ceh_sys, avg_sys depend only on Q and N)

Fields

• Q: Markov transition matrix
• μ, ν: boundary probability densities
• π: stationary distribution of Q
• N: number of geodesic time steps
• ceh_sys: factorised continuity-equation linear system
• avg_sys: factorised time-averaging linear system

ErbarBundle{M, T}

A pair (cache::ErbarCache, vector::ErbarVector) representing the complete state of a graph OT problem: the problem metadata and the current iterate.

Type parameters {M, T} mirror those of ErbarVector so that accessing bundle.vector in hot-path code returns a concretely-typed value, enabling full specialisation of the Chambolle-Pock inner loop without dynamic dispatch.

Constructors:

• ErbarBundle(Q, μ, ν, N) — compute π internally
• ErbarBundle(Q, π, μ, ν, N) — use pre-computed stationary distribution
• ErbarBundle(cache, vector) — assemble from existing components

combine!(c::ErbarVector, a::ErbarVector, b::ErbarVector, α::Number, β::Number)

Compute the linear combination c = α·a + β·b in-place across all seven component arrays of the ErbarVector (ρ, m, θ, ρminus, ρplus, ρ_avg, q). No intermediate arrays are allocated; all operations use broadcasting with .=.

combine!(c::ErbarBundle, a::ErbarBundle, b::ErbarBundle, α::Number, β::Number)

Bundle-level wrapper for combine! on ErbarVectors. Computes c.vector = α·a.vector + β·b.vector in-place; the cache fields of the bundles are not modified.

assign!(c::ErbarVector, a::ErbarVector)

Copy all component arrays of a into c in-place (deep value copy, not pointer aliasing). Equivalent to c = copy(a) but reuses c's existing allocations. Covers all seven fields: ρ, m, θ, ρminus, ρplus, ρ_avg, q.

assign!(c::ErbarBundle, a::ErbarBundle)

Bundle-level wrapper for assign! on ErbarVectors. Copies a.vector into c.vector in-place; the cache fields are not modified.