Cost-augmented Schrödinger bridges on graphs are exactly solvable: a Feynman-Kac tilt replaces learned control

Cost-augmented Schrödinger bridges on graphs are exactly solvable: a Feynman–Kac tilt replaces learned control Akshay Balsubramani akshay@vac.bio

arXiv:2610.02195v1 [cs.LG] 1 Oct 2026

Abstract The generalized Schrödinger bridge on a graph moves mass between two distributions while charging a cost for the states visited. It has been approached by learning the rates of a controlled continuous-time Markov chain, with a temporal-difference penalty that restores the cost. A state cost folds into the reference process as a Feynman–Kac tilt. The cost-augmented bridge is then a plain bridge against the tilted reference, and the penalty is unnecessary. The bridge is computed exactly by alternating two endpoint rescalings, each one sparse matrix-exponential application; nothing is discretized in time or learned. The alternation converges at a rate set by the endpoint coupling alone. For a quadratic congestion cost on time-averaged occupancies, damped best response around the exact bridge is gradient descent on a strongly convex function, and its residual bounds its error. On a protein-folding model, a free-energy cost lowers the expected barrier of the folding paths. On the learned approach’s road network, roll-outs of the exact bridge match the target within sampling error, and on networks with millions of intersections its memory grows linearly.

1

Introduction

Transporting probability mass across a graph under operational constraints is the common form of a large class of problems: routing shipments through a supply network, assigning workers to tasks, steering a molecular simulation from one metastable basin to another. Routing and steering ask for a time-indexed plan, which a static coupling does not supply. For an assignment, posed on a layered graph, the plan reduces to the coupling. In all three problems, each state visited along the way adds to the cost. The Schrödinger bridge [24] is a time-indexed plan. Among path laws with the given endpoint distributions, it minimizes the relative entropy to a reference process. A running cost on the visited states turns it into the generalized bridge. In continuous state spaces, generalized Schrödinger bridge matching learns a control that accounts for task-specific state costs [26]. On a graph, recent work poses the generalized bridge as a learning problem, because it judges a direct solve prohibitive on large graphs [40]. Neural networks parameterize the potentials that define the rates of a controlled continuous-time Markov chain. An iterative-proportional-fitting objective fits them to the endpoint marginals, but the running cost, kept outside the reference, cancels out of that objective (Section 2). A temporal-difference penalty, adapted from a continuous-space solver for mean-field games [25], restores the cost. The penalty adds a weight to be tuned, swept over six values in the supply-chain experiments of [40]. That weight trades the running cost against the terminal constraint: a larger weight lowers the peak occupancy and worsens the terminal match to the target. The penalty’s continuous-space form borrows target networks and replay buffers from deep reinforcement learning to help stabilize training [25]. A running cost f that depends only on the state folds into the reference measure as an exponential tilt. The generalized bridge is then identically a plain bridge against the tilted reference (Section 2). With the cost inside the reference nothing cancels, so no temporal-difference term is needed. The tilt shifts the diagonal: it replaces the reference generator Q by Q − diag(f ). The Schrödinger system then becomes a pair of linear ordinary differential equations, solved by matrix exponentials. The shift leaves every edge rate in place. The optimal control therefore reweights each edge rate by a ratio of Schrödinger potentials, so the controlled chain moves only along the edges of the graph. The bridge follows by alternating two endpoint rescalings, each one sparse matrix-exponential application (Section 3). Both the tilt and the potential-ratio form of the control are classical; Section 2 gives their sources. On a finite graph, the two together let sparse linear algebra solve the problem exactly. Here exact means that nothing is discretized in time and nothing is learned: the solver applies matrix exponentials of the tilted generator to vectors numerically. In continuous space, a bridge with a quadratic state cost is exactly solvable in a stronger sense, since its kernel has a closed form [39].

1

The identity does not cover costs that depend on the marginal; congestion is the canonical example. For a congestion cost proportional to the sum of the squared time-averaged occupancies of the congestible states, we compute the congestion equilibrium by damped best response, with the exact bridge as the inner oracle. Damped best response is gradient descent on a strongly convex function whose curvature has explicit bounds. It therefore converges from any start when the step lies below a threshold set by the upper bound. Near the equilibrium it also converges for larger steps, up to a threshold set by the curvature there. Whatever the step, the residual of each iterate bounds its distance to the equilibrium (Section 4). The exact solver returns the optimal controlled rate in the form targeted by the learned approach, so it can replace that approach directly. We evaluate it on that approach’s data: a Markov state model of chignolin folding, rebuilt from the molecular-dynamics trajectories of its rare-event scenario, and the Washington DC road network of its supply-chain scenario. On both, the exact bridge meets its endpoint marginals. On the road network, roll-outs of its controlled rate match the target within sampling error. On the same network, the tested dampings converge below the threshold set by the curvature at the equilibrium and fail to converge above it. Where the largest eigenvalue of the curvature governs convergence, the residual contracts at the rate predicted from that eigenvalue. On the assignment construction of [40], the exact solution recovers the minimum-cost assignment on every instance at concentrations within the band of that paper’s main table. The learned solver misplaces a pair at the largest size tested, where the exact solution still recovers every instance. On chignolin, the controlled process folds in every roll-out. A free-energy cost lowers the expected barrier of the folding paths. The fitting iteration converges at a rate that depends on the graph only through the second singular value of the endpoint coupling. On road networks of up to 2.07 million intersections, the solver’s memory grows in proportion to the node count.

2

The tilt identity

Let G = (X, E) be a directed graph on a finite node set with |X| = n, and let Q be aPgenerator supported on E: Qxy ≥ 0 is the rate of the transition x → y for x ̸= y, Qxy = 0 whenever (x, y) ∈ / E, and y Qxy = 0. Write R for the law on [0, τ ] of the chain with generator Q started from a distribution of full support, let µ, ν be distributions on X, and denote by qt the time-t marginal of a path measure q. Every q considered below starts at µ, so the initial law of R changes the objective only by a constant. Let f : X → R be a running cost. The generalized Schrödinger bridge is the path measure Z τ  q ⋆ = arg min Eq f (Xt ) dt + KL(q ∥ R) (1) q : q0 =µ, qτ =ν

0

Rτ e by dR/dR e Theorem 2.1 (Tilt identity). Define the tilted path measure R ∝ exp − 0 f (Xt ) dt . Then Z τ  h Rτ i e − log Z, Eq f (Xt ) dt + KL(q ∥ R) = KL(q ∥ R) Z = ER e− 0 f

(2)

0

e for every q absolutely continuous with respect to R. Consequently q ⋆ is the plain Schrödinger bridge between µ and ν against R, tA the normalized path measure of the semigroup e with A = Q − diag(f )

(3)

The proof is in Appendix A.1. A generator has zero row sums; those of A are −f , so A is not itself a generator. For f ≥ 0 the semigroup etA is e is that sub-Markov: the tilted process is killed at rate f , so its path measure has total mass Z. The tilted path measure R measure divided by Z. The mass defect does not matter. A Schrödinger bridge is defined against its reference only up to normalization, so the two endpoint rescalings absorb the defect exactly. In (2) it appears only as the q-independent e is supported on the same edge set as R. constant log Z. The tilt is also diagonal, so it changes no off-diagonal entry and R τA In particular, for τ > 0, e has a positive (x, y) entry exactly when y can be reached from x by transitions of positive rate, whatever the cost. We assume that every state in the support of ν can be reached from every state in the support of µ. The block of eτ A on the two supports is then positive, so the bridge exists and its endpoint scalings are positive and unique up to reciprocal constant factors [28]. Every instance in Section 5 satisfies this assumption. The optimal control inherits the edge set.

2

Corollary 2.2 (Optimal control and support preservation). Let (φt , φ̂t ) be the Schrödinger potentials of the bridge against e so that qt⋆ = φt φ̂t . The optimal controlled generator is R, u⋆t (x, y) = Qxy

φt (y) , φt (x)

x ̸= y

(4)

and Qxy = 0 implies u⋆t (x, y) = 0. The proof is in Appendix A.1. Equation (4) is an executable policy: a particle at x consults only its out-neighbors, so the graph’s topology is respected by construction, with no penalty enforcing it. It coincides with the form targeted by the learned approach [40], so the exact solver substitutes for it without reformulating the problem. Both the tilt and the potential-ratio form of the control are classical [24]. In linearly solvable Markov decision problems, whose control cost is a relative entropy, an exponential transform of the value makes the Bellman equation linear in the tilted transition matrix, with optimal transitions that reweight the uncontrolled ones by the transformed value of the next state [42]. Path-integral control linearizes its Bellman equation by the same transform [19]. For a diffusion, a positive solution of a generalized Schrödinger system, a pair of linear equations coupled through their boundary values, yields an optimal control of the bridge with a running cost [9]. In continuous space, iterative Markovian fitting has been extended to the generalized bridge of [26] by recasting it as a plain bridge against a Feynman–Kac transform of Brownian motion [31]. Discrete-time bridges on networks are taken against priors whose rows need not sum to one, with transitions of the same form [8]. On a graph in discrete time, transport of paths regularized by relative entropy, with edge costs folded into the reference transition matrix W , reduces to scaling a block of its fundamental matrix (I − W )−1 [15]. For random walks on directed graphs, a path cost enters as a relative entropy to the prior path law reweighted by the exponential of minus the path cost, a measure that need not be normalized [12]. Concurrent work states the tilt for controlled Markov jump processes, whose endpoint scalings are fitted by iterative proportional fitting when every jump rate is controlled [37], and for discrete-time chains on a finite state space [27]. The absence of a temporal-difference term. The learned approach holds the reference process fixed and adds the running cost to its objective as a separate term [40]. Written against the untilted reference, the potentials obey ∂t φt = −Qφt + f φt and ∂t φ̂t = Q⊤ φ̂t − f φ̂t , in which f appears with opposite signs and cancels from ∂t log qt⋆ = ∂t log φt + ∂t log φ̂t . The two cost terms cancel in the same way in the iterative-proportional-fitting objective of the learned approach, so on its own that objective enforces the endpoint marginals but not consistency with the cost [40]. A temporal-difference penalty restores that consistency. It is the graph form of a loss from a continuous-space solver for mean-field games, whose forward-backward stochastic differential equations, once discretized, resemble temporale the cost is already inside A, so iterative proportional fitting sees it and the difference learning [25]. Written against R, cancellation has nowhere to occur. Both descriptions are of the same problem; only the second leaves the cost visible to the algorithm.

3

An exact continuous-time solver

Theorem 2.1 reduces the bridge to a plain bridge against the path measure of the semigroup etA , normalization aside. The ⊤ Schrödinger system is then linear, ∂t φt = −Aφt and ∂t φ̂t = A⊤ φ̂t , with solutions φs = e(t−s)A φt and φ̂t = etA φ̂0 . Iterative proportional fitting alternates the two endpoint conditions, φ̂0 ← µ ⊘ φ0 ,

⊤ φτ ← ν ⊘ eτ A φ̂0 ,

φ0 ← eτ A φτ

(5)

where ⊘ is entrywise division. Each half-step is one application of a matrix exponential to a vector, computed by a scaled truncated Taylor series that never forms eτ A and touches A only through sparse products, whose number grows at most in proportion to ∥τ A∥1 [1]. Each application therefore takes work at most proportional to |E| ∥τ A∥1 once that norm exceeds a small constant. The whole time axis is handled at once, with no grid: φt at any t is recovered by a single further application. Sinkhorn iterations that apply a heat kernel to vectors, without ever forming the kernel matrix, compute transport on meshes and images [36] and on graphs [17]. There the kernel is symmetric and the output is a static coupling. The fitting iteration (5) is the Sinkhorn iteration on the block of eτ A between the two supports. Only the entries of the scalings on the two supports enter (5), so when a support is small the same iteration can run on the block itself: its columns 3

take one application each, after which every iteration is a product with a small dense matrix. Either way the iteration converges linearly: near its fixed point the error decays by a factor σ22 per iteration, where σ2 is the second singular value of −1/2 −1/2 Dµ ΠDν , Π is the optimal endpoint coupling, and Dµ , Dν are the diagonal matrices of the two marginals on their supports [21, 41]. Started from constant scalings, the number of iterations to a tolerance ε is therefore log(1/ε)/ log(1/σ22 ) up to an offset set by the starting point. On the 11 road-network problems of Section 5.4, the counts are 0.81 to 1.07 times this law (Figure 1c). The singular value σ2 approaches one as the coupling approaches a block-diagonal one, whose blocks do not communicate; the graph enters the rate only through Π. Because Π is the solution itself, σ2 is known only once the problem is solved, so the law explains an iteration count only after the fact. The iteration stops when the logarithm of every scaling on the support of µ changes by less than ε. A change measured against the largest scaling, the rule behind the solves of Sections 5.2 and 5.3, can fall below ε while scalings many orders of magnitude smaller are still moving (Section 5.4). Those solves read their marginal errors out of the potentials, so each reported accuracy is measured directly, whatever the stopping rule.

seconds per iteration

peak allocation (MB)

102 101 104

105

nodes n

(b)

(c)

102

104

fitting iterations

(a)

103

101 100 10 1

106

107

|E| max|Axx|

109

103

102

102

103

log(1012) /log(1/ 22)

104

Figure 1: The number of fitting iterations follows the second singular value of the endpoint coupling, while memory and time per iteration grow with the size of the problem. Road networks of the ROAD-FLOW benchmark from Washington DC (n = 9,559) to Texas (n = 2,073,870); 40 sources and 40 targets within 50 hops of one random center, running cost uniform on [0, 0.5], horizon 100. (a) The solver’s peak allocation grows in proportion to the node count, 346 to 412 bytes per node; dashed, proportional to n. (b) Wall-clock time per fitting iteration against the arc count times τ maxx |Axx |, for the seven networks timed in one container (filled) and for Washington DC with the endpoints of [40] at horizons 25 to 400, timed together in a second container (open). The dashed lines, one proportional fit per container, nearly coincide. Washington DC, the one network timed in both containers, costs 2.0 × 10−8 seconds per arc and unit of τ maxx |Axx | in the first and 4.1 × 10−8 in the second. The networks lie within a factor 2.4 of one constant (Washington DC, timed again at the end, within 1.1% of its first timing); the horizons lie within 3% (horizon 25, timed again at the end, within 1.3% of its first timing). (c) Fitting iterations to a change of 10−12 in the logarithm of the scalings, against log(1012 )/ log(1/σ22 ): the seven networks (filled), four further endpoint sets on New Hampshire (open circles) and the horizons of (b) (open squares); dashed, equality; dotted, a factor 2 either way. At the end of every solve the error contracts per iteration by σ22 to within 0.03%. The counts are 0.81 to 1.76 times the prediction and in its order on the networks (rank correlation 1.00). Every solve in (c) meets both endpoint marginals to a total variation below 3 × 10−13 . Across a horizon, the potential’s overall scale can change by a factor of up to e∥A∥τ , which exceeds the range of double precision once the cost is large relative to 1/τ . Propagation is therefore chunked, with the vector renormalized after each chunk. This is exact here: both the fitting iteration and the marginal φt φ̂t are invariant to a global rescaling of φ. Evaluating (4) across a free-energy barrier likewise requires the ratio in the log domain, since φ can span tens of orders of magnitude between basins and vanishes outside the target support as t approaches τ .

4

Marginal-dependent costs

A congestion cost depends on the marginal, so Theorem 2.1 does not apply: the tilt would depend on the solution. A penalty affine in the logarithm of the number of drivers on a route depends on the marginal and still leaves a linear problem: the 4

equilibrium of the resulting mean-field routing game solves a linearly solvable Markov decision process by one backward recursion [38]. In continuous space, a variational mean-field game with diffusion and a quadratic Hamiltonian is equivalent to minimizing, over path measures with a fixed initial marginal, a relative entropy plus a running cost of the marginals Rτ and a terminal cost [7]. We charge congestion instead on the occupancy averaged over the horizon, p̄ q = τ −1 0 qt dt, P through the quadratic potential F (p) = β2 x∈C p(x)2 on a set C ⊆ X of congestible states, for instance every state other than the endpoints. The objective KL(q∥R) + τ F (p̄ q ) is strictly convex in q over the convex set cut out by the endpoint constraints, since the relative entropy is strictly convex and τ F (p̄ q ) is convex, p̄ q being linear in q. Under the reachability assumption of Section 2 that set is not empty, so the minimizer q ◦ exists and is unique. Its linearization at a given q is the running cost β p̄ q on C, which depends on the state alone, so the best response to it is a bridge of Section 3. We reach q ◦ by damped best response. For a running cost g on C write p̄(g) ∈ RC for the time-averaged marginal of the bridge with cost g, restricted to C, and Φ(g) for the optimal value of (1) with f = g. The time average is computed from the potentials by Gauss–Legendre quadrature, whose error decays faster than any power of the number of points because t 7→ φt φ̂t is entire, analytic on the whole complex plane. From the current occupancy p ∈ RC , the step solves the bridge with cost βp and mixes, p ← (1 − α) p + α p̄(βp). Proposition 4.1 (Damped best response is gradient descent). Assume that every state in the support of ν can be reached from every state in the support of µ. The damped best-response step is a gradient step of size α on h(p) = 12 ∥p∥22 −

Φ(βp) βτ

(6) ◦

which satisfies I ⪯
∇2 h ⪯ (1 + βτ /2)I and has the restriction p◦ of p̄ q to C as its unique minimizer. For every α ∈ 0, 4/(2+βτ ) the iterates satisfy ∥pk −p◦ ∥2 ≤ ρk ∥p0 −p◦ ∥2 with
ρ = max{1−α, α(1+βτ /2)−1}. Let m ≤ L 2 ◦ be the extreme eigenvalues of ∇ h(p ) − I, and for α ∈ 0, 2/(1 + L) let ρ◦ = max{1 − α(1 + m), α(1 + L) − 1} < 1. For every ϵ ∈ (0, 1 − ρ◦ ) there is a ball around p◦ that the step maps into itself, and on it ∥pk+1 − p◦ ∥2 ≤ (ρ◦ + ϵ) ∥pk − p◦ ∥2 . Every p satisfies ∥p − p◦ ∥2 ≤ ∥p − p̄(βp)∥2 . The proof is in Appendix A.2. The last inequality of the proposition, a certificate, holds whatever the damping, so the residual reported with every equilibrium bounds its error, up to the error of the inner bridge solves. The proof bounds ∇2 h − I by the smaller of βτ maxx∈C p̄(x) and βτ /2, which is small when no state holds much of the occupancy. The local rate ρ◦ is smallest at α = 2/(2 + m + L); above 2/(1 + L) the step’s Jacobian at p◦ has the eigenvalue 1 − α(1 + L) < −1, so the equilibrium is unstable. The global range is narrow, 0.016 at βτ = 250 and 0.0016 at βτ = 2500. The converging dampings of Section 5.3 exceed it by factors of 22 to 449 while staying below 2/(1 + L), where the proposition covers them only near p◦ . For βτ from 10 to 250 the global range is nearly the best possible among ranges that depend on βτ alone. On the path a–s–t–b with µ = δs , ν = δt , τ = 1 and C = {a, b}, where s and t jump to each neighbor at rate 104 and a and b return at rate 0.05, the equilibrium occupancy splits evenly between a and b. For each such βτ its curvature L comes within 1.5% of βτ /2, so the local threshold comes within 1.5% of the global range. Section 5.3 therefore measures the convergence of its dampings from where they start, the occupancy of the plain bridge. That section also measures L on a road network and checks the threshold. Started from the occupancy of a bridge, damped best response is the generalized conditional-gradient method applied to the primal objective, which with the step 1/(k + 1) is fictitious play in potential mean-field games [23]. By Fenchel duality it coincides with gradient ascent on the dual, with the primal and dual of [2] exchanged; ours is the case of a quadratic regularizer. In traffic assignment, travelers in logit stochastic user equilibrium choose routes by a logit model of their travel times. That equilibrium solves a strictly convex program in the link travel times. The gradient of the program is the difference between the link flows that produce those times and the flows assigned by the route-choice model at those times [35]. The method of successive averages steps along minus that gradient, so it is a form of steepest descent on that program. In Markovian route choice a route is built by choosing an arc at every node with a discrete choice model. For that model the dual of the equilibrium problem is a smooth concave function of the link costs whose gradient is the assigned link flows minus the inverse of the link cost function [6, 33]. The function βh of Proposition 4.1 is the travel-time program, with Φ/τ in the role of its expected minimum route cost and link costs βp. Damped best response is then the method of successive averages run with a constant step. For entropic optimal transport with a convex functional of the coupling, the continuous-time form of the same construction, a flow toward the best response computed by one Sinkhorn solve against an exponentially tilted reference, converges exponentially [20]. In the mean-field Schrödinger bridge, diffusing particles interact through a nonlocal potential, which makes the problem nonconvex [3, 13]. A fixed-point 5

iteration on its marginal flow, whose inner steps are Sinkhorn-type solves with the interaction frozen, converges locally under contraction conditions [13]. Linear convergence for a small enough constant step is known for generalized conditional gradient [14] and for its mean-field-game form [23]. For the method of successive averages in logit stochastic user equilibrium, a constant step below a bound set by the most negative eigenvalue of the logit mapping’s Jacobian at the equilibrium converges linearly near it [5]. Near that equilibrium the fixed-point residual bounds the distance to the equilibrium path flows, up to a factor set by the Jacobian of the logit mapping [4]. For the bridge the curvature bound is explicit, which gives the global step range and the certificate above, valid at every p. Entropic fictitious play, on an objective of the same form (a relative entropy plus a convex function of linear moments), converges linearly up to a floor that shrinks with its step [30]. In entropy-regularized mean-field games with finite state and action spaces, the undamped fixed-point iteration is a contraction once the temperature exceeds a threshold set by Lipschitz constants of the game [11]. Here the undamped step α = 1 lies below the local threshold 2/(1 + L) exactly when L < 1. For convex graph-structured tensor problems, generalized Sinkhorn iterations, a different algorithm that ascends the dual block by block, converge globally and, under stricter assumptions, R-linearly (with an error bounded by a geometric sequence) [34].

5

Experiments

We evaluate on the three scenarios of [40] and use its metrics. The rare-event and supply-chain scenarios run on that paper’s data. The first is a Markov state model of the miniprotein chignolin, a finite-state Markov chain estimated from molecular-dynamics trajectories. We rebuild it from the trajectories [18] used in that paper, by the construction in its appendix. The second is the Washington DC road network of the ROAD-FLOW minimum-cost-flow benchmark [22], with that paper’s supply and demand nodes. Because that paper released neither its model nor a complete specification of its reference rates and simulates in discrete time, values across the two papers are comparable in kind, while comparisons within a scenario are like-for-like. The assignment task reconstructs that paper’s graph construction on cost matrices and marginals drawn here, so it is compared at matched concentration.

5.1

Assignment

For an assignment of size K, the construction places an intermediate node Exy on each of the K 2 pairs, with f (Exy ) = Cxy . We draw each cost uniformly from the integers 1, . . . , s − 1, where s is the cost scale swept below, and use uniform endpoint marginals, so the coupling read off the flux is a soft assignment. The layered graph admits exactly two transitions. We charge the cost once per visit to Exy , the discrete-time form of the tilt: the argument of Theorem 2.1 holds with the time integral replaced by a sum, so the bridge is a two-step fitting iteration on the kernel e−C . The concentration of the soft assignment is measured by its row entropy, the mean entropy of the coupling’s rows after each is normalized to sum to one. Concentration is governed by the spread of C relative to the entropic term, so comparing at a single cost scale would not be like-for-like; we sweep the scale and match the row entropies reported in [40]. That paper prints its assignment metrics twice and the two copies disagree at its largest size: the main table gives a row entropy of 0.17 at K = 20 and the appendix table gives 0.22, with the other five columns identical. We match the main table’s band, 0.07 to 0.17, and also report the appendix table’s looser value. Inside that band the exact solution recovers a cost-optimal assignment on every instance at every size tested (Figure 2). Its hard read-out, the permutation of largest total coupling mass, agrees with the optimum on every row. At the appendix table’s 0.22 the recovery is no longer uniform: 19 of the 20 instances at K = 20 are cost-optimal. The exception is the most diffuse of the 20, with a row entropy of 0.54 against a mean of 0.21 over the instances at that scale. At K = 20 the main table’s own value, 0.17, lies between the row entropies of the two nearest scales tested, 0.14 with all 20 instances cost-optimal and 0.21 with 19. The learned solver reports 90% agreement and a cost of 546 against an optimum of 538 at K = 20: the exact solution introduces no learned approximation, so the sweep isolates the error contributed by the entropic regularization alone. At low cost scales (that is, at high entropy) even the exact bridge’s hard read-out departs from the optimum, which locates the failure in the regularized objective itself.

5.2

Rare-event steering

The Markov state model is built on 500 clusters of the backbone dihedral angles of the 10-residue miniprotein chignolin, with transition probabilities estimated at a lag of 1 ns from five trajectories of 1 µs at 300 K. Its free energy is F = − log π in units of kB T , measured from the lowest basin state. A state belongs to the unfolded basin when more than 80% of its 6

K=6

K=8

K = 10

K = 20

1.0 0.9 0.8 0.7 0.6 0.5

0.0

0.2

0.4

0.6

0.8

1.0

row entropy of the coupling (nats)

matched to the band

fraction of instances cost-optimal

mass on the optimal assignment

(a)

1.2

(b)

learned solver's band

1.0 0.9 0.8 0.7 0.6 0.5 0.0

0.2

0.4

0.6

0.8

1.0

row entropy of the coupling (nats)

1.2

Figure 2: Inside the concentration band of the main table of [40], the exact bridge recovers the minimum-cost assignment on every instance. Departures from the optimum appear only at higher entropy. (a) Mass the soft assignment places on the optimal assignment and (b) fraction of instances whose hard read-out is cost-optimal, against the row entropy of the coupling, for K = 6, 8, 10, 20. Each point averages 20 random integer cost matrices at one cost scale; along each curve the scales are 25, 50, 100, 200, 400 and 800, from right to left. The shaded band is the row-entropy range 0.07–0.17 reported for the learned solver in the main table of [40]. Ringed points are the scale matched to it at each size. There the hard read-out agrees with the optimum on every row of every instance. The largest endpoint-marginal violation over the 20 instances is 1.8 × 10−6 , 9.4 × 10−7 , 1.9 × 10−6 and 1.0 × 10−6 for K = 6, 8, 10, 20, set on this layered construction by the budget of the two-step fitting iteration. These violations change no read-out. That paper’s appendix repeats the table with a row entropy of 0.22 at K = 20 in place of 0.17; at the scale matching 0.22, 19 of the 20 instances at that size are cost-optimal. Mass on the optimum falls monotonically as the entropy rises. The departures in (b) scatter at a sample size of 20 instances. Bars are 95% intervals, bootstrap in both coordinates in (a) and Wilson in (b).

7

frames lie beyond 0.5 nm of Cα root-mean-square deviation from the folded structure, and to the folded basin when more than 90% lie within 0.12 nm. At the folded threshold of 0.09 nm stated in [40] no state of the rebuilt model qualifies. The threshold 0.12 nm is the smallest on a grid of 0.01 nm at which the folded basin holds at least a tenth of the frames (32% of them, against 4% at 0.11 nm). At that threshold the trajectories contain five folding events and one unfolding event, so the model’s folding kinetics rest on few transitions, as they must for any model built from these data. The construction removes transitions between clusters whose mean deviations differ by more than 0.1 nm, after which the model keeps the 485 clusters of the largest connected set that remains. The running cost is the free energy scaled by a dimensionless strength λ, so λ = 0 recovers the plain bridge. A path’s barrier is the maximum free energy along the path. The horizon, 58 ns, is calibrated so the uncontrolled process folds 1% of the time, against 0.55% for the uncontrolled chain of [40] over its horizon. Bridge matching has also been used to learn Doob-transformed rates on Markov state models of other systems [43]; on such a model the exact bridge computes those rates directly (Corollary 2.2). At every λ the bridge meets both endpoint marginals to a total variation below 10−16 . The terminal marginal of the bridge is the folded basin’s equilibrium weights, so it folds in every one of 4,000 roll-outs at every λ, against 99.70% and 99.36% reported for the learned bridge without and with the cost. The distribution of the barrier is computed exactly: the probability that a path stays in the sublevel set {F ≤ h} comes from the tilted generator restricted to that set, with the bridge’s endpoint scalings. The expected barrier falls monotonically as the cost strengthens, from 6.20 at λ = 0 to 5.64 at λ = 30. At every strength the barrier is stochastically smaller than under the plain bridge: its distribution function lies on or above the plain bridge’s at every level (Figure 3). At f = 0 the bridge and the uncontrolled process conditioned on folding differ only through their endpoint weights (Corollary 2.2); here those weights raise the expected barrier from 5.78 to 6.20. The running cost brings it back below the conditioned process at a strength between 1 and 3: the expected barrier is 5.83 at λ = 1 and 5.73 at λ = 3. Transition path theory treats a different ensemble, the pieces of an equilibrium trajectory that leave the unfolded set and next enter the folded one [29]. It describes their distribution, probability current and rate through committor functions, the probabilities of reaching each set before the other. On Markov state models of folding, the same theory gives the net flux of folding trajectories [32]. On the rebuilt model no folding path has a barrier below 4.28, the free energy of the lowest unfolded state, so its barriers and the 1.80 and 1.39 reported for the learned bridge measure different models.

(a)

(b)

6.0 5.8

1.0

bridge, exact bridge, roll-outs

uncontrolled, given folding 0

0.1

0.3

1

0.8

Pr(barrier h)

expected barrier (kBT)

6.2

0.6 0.4

uncontrolled, given folding bridge, = 0 bridge, = 1 bridge, = 30

0.2 3

10

strength of the free-energy cost

0.0

30

4

6

8

barrier level h (kBT)

10

Figure 3: On a Markov state model of chignolin folding, a free-energy running cost lowers the expected barrier monotonically. (a) Expected barrier, the maximum free energy along a path, measured from the lowest basin state, against the strength λ of the running cost λF (λ = 0 is the plain bridge; the axis is linear below 0.1 and logarithmic above). Filled markers are exact, from the tilted generator restricted to each sublevel set of F . Open markers average over the roll-outs that fold; they lie within 0.02 of the exact values. The dashed line is the uncontrolled process conditioned on folding, exact; its roll-out estimate from its 42 folding trajectories out of 4,000 is 5.60. (b) The exact distribution of the barrier, Pr(maxt F (Xt ) ≤ h).

8

5.3

Supply-chain routing with congestion

The network has 9,559 intersections and 29,682 arcs whose capacities follow road categories. Reference rates are proportional to capacity and scaled to a unit mean exit rate, so time is measured in mean jump times. The horizon is 100, the number of steps of [40]. Six supply nodes ship 6,360 units in total to six demand nodes. The congestion cost falls on every node other than the supply and demand nodes, as in that paper, which penalizes the occupancy at each step where Section 4 penalizes its average over the horizon. Both endpoint marginals are met to a total variation below 10−12 (Figure 4). Roll-outs of the controlled rate from 5,000 particles end at a total variation of 0.016 from the target, within the range of exact draws of 5,000 particles from the target (0.012 on average, 0.019 at the 95th percentile). On the same network and endpoints the learned solver of [40] ends at 0.03 and its attraction-flow baseline at 0.27. Its graph Schrödinger-bridge baseline, a dynamical formulation whose minimizer satisfies a boundary-value Hamiltonian flow on the probability simplex [10], did not complete on this network because of memory exhaustion. Near the horizon the rate into the target grows like 1/(τ − t), so a uniform time grid saturates there however fine it is. The roll-outs therefore shrink the step geometrically over the last five time units. Congestion spreads the load, weakly at βτ = 250 and more at βτ = 2500, where the mean occupancy of the 100 most occupied nodes falls from 13.4 to 10.4 particles in 5,000 (Figure 4). The largest flow along one arc within one unit of time, 1.09 times its capacity for 5,000 particles (1.38 for the instance’s 6,360 units), runs into the largest demand node during the last unit of time, where the pinned terminal marginal concentrates the arrivals. Because the congestion cost excludes the endpoints, it moves this flow only to 1.08 and 1.05. At βτ = 250 the largest eigenvalue of ∇2 h − I at the equilibrium, measured by Lanczos iteration on finite differences of the best response, is L = 0.1547, a tenth of the bound βτ maxx p◦ (x) = 1.64 from the proof of Proposition 4.1. The critical damping 2/(1 + L) = 1.73 exceeds one. The plain best response converges, contracting by 0.1533 per round over its last three rounds, within 1% of the local rate L; damping 0.35 converges in 42 rounds (Figure 5). At βτ = 2500, L = 1.51 and the critical damping falls to 0.80: the plain best response settles into a two-cycle at a residual of 0.023, while damping 0.35 converges in 42 rounds. At 0.9 times the critical damping, an over-relaxed step of 1.56 at βτ = 250 and a step of 0.72 at βτ = 2500, the residual contracts by 0.800 per round at both strengths, as α(1 + L) − 1 predicts. At 1.1 times it, the residual grows by 1.16 per round at βτ = 250, near the 1.20 the linearization gives; at βτ = 2500 it settles into a two-cycle.

5.4

Scale

The solver’s memory is a few copies of the sparse generator and a constant number of length-n vectors. Each fitting iteration is two sparse matrix-exponential applications. We measure both on the seven road networks of the ROAD-FLOW benchmark [22], from Washington DC (9,559 intersections) to Texas (2,073,870), with 40 sources and 40 targets drawn within 50 hops of one random center, a running cost uniform on [0, 0.5], and horizon 100 (Figure 1). The peak allocation is 346 to 412 bytes per node, 743 MB on Texas. Timed in one container, an iteration costs 2.0 to 4.6 × 10−8 seconds per arc and unit of τ maxx |Axx |, more on the larger networks, so over a 172-fold range of arc counts the time per iteration stays within a factor 2.4 of proportional to |E| τ maxx |Axx |. On Washington DC with the endpoints of [40], timed together in a second container, the time per iteration grows in proportion to the horizon from 25 to 400, to within 3%; horizon 25, timed again at the end, comes within 1.3% of its first timing. The horizon needed to connect distant endpoints therefore sets the cost of a problem as much as the size of the network does. Washington DC, the one network timed in both containers, costs 2.0 × 10−8 seconds per arc and unit of τ maxx |Axx | in the first and 4.1 × 10−8 in the second. The number of iterations follows the endpoint coupling. On 11 problems, the seven networks and four further endpoint sets on New Hampshire, the error contracts per iteration by σ22 to within 5 × 10−5 . The counts are 0.81 to 1.07 times log(1/ε)/ log(1/σ22 ) and in the same order. They range from 561 on Florida to 7,313 on a New Hampshire endpoint set whose coupling nearly splits in two (σ22 = 0.997); Texas takes 1,973. With 40 sources and 40 targets the iteration runs instead on the endpoint block (Section 3), where each iteration is a product with a 40 × 40 matrix: the Texas problem takes 28 minutes on four cores, 27 of them spent on the 40 backward applications that build the block. Where the iteration also ran on full vectors, stopped by a change measured against the largest scaling, on the networks up to Wisconsin (519,157 intersections), both marginals are met below 5 × 10−10 . The terminal scalings of those solves agree to 7 × 10−6 in the logarithm with those of a block iteration stopped by the same rule. At 5 × 104 nodes of its own test graphs, the learned solver of [40] uses 1.4 GB and about one hour; on the 49,109 intersections of Delaware the exact solver’s peak allocation is 18 MB and its block iteration meets both marginals in 26 seconds, on different hardware. The stopping residual matters when the endpoints lie many hops apart relative to the horizon. At horizon 25 the endpoint block of the compared paper’s endpoints spans 100 orders of magnitude. A change measured against the largest 9

(a) endpoint total variation 1.0 × 100

uncontrolled

(b) roll-out total (c) top-100 variation occupancy

(d) peak occupancy

not computed

not computed

(e) largest arc flow over capacity

179

not computed

bridge, f = 0

8.7 × 10 13

0.016

13.4

246

1.09

= 250

8.3 × 10 13

0.016

13.0

245

1.08

= 2500 2.5 × 10 13

0.015

10.4

236

1.05

10 12

100

0.00 0.01 0.02 0

10

20 0

200

0

1

Figure 4: On the compared paper’s road network the bridge meets both endpoint marginals, exactly in its read-outs and within sampling error in its roll-outs. Washington DC roads (n = 9,559 intersections, 29,682 arcs) with the six supply and six demand nodes of [40], horizon 100 mean jump times; four processes, one panel per metric, with exact values at the bar ends. (a) The larger total variation of the two exact endpoint read-outs; for the uncontrolled chain, which has no endpoint constraint, the terminal total variation. (b) The terminal total variation of 5,000 roll-outs of the controlled rate (4). The gray band spans exact draws of 5,000 particles from the target up to their 95th percentile (0.019), and the dotted line marks their mean (0.012). (c) to (e) The compared paper’s crowding metrics in expectation, per 5,000 particles: the mean time-averaged occupancy of the 100 most occupied nodes, the largest occupancy of one node at one unit time (both excluding the endpoints), and the largest flow along one arc within one unit time over its capacity. The congestion cost of Section 4 (bottom two bars) spreads the load. For the uncontrolled chain only the terminal total variation and the peak were computed. scaling then falls below ε while scalings tens of orders of magnitude smaller are still moving, leaving a terminal total variation of 0.094. Under the logarithmic rule, as for the counts above, the iteration converges in 837 iterations, 1.76 times the asymptotic count. Across horizons 25 to 400 the counts are 1.09 to 1.76 times the asymptotic count.

6

Discussion

The identity of Section 2 is elementary: a change of reference measure converts a problem posed as learned stochastic control into one solved by sparse linear algebra. The congestion cost of Section 4 adds an outer iteration whose residual bounds its distance to the equilibrium at every step. The exact solver requires the reference generator explicitly. So does the learned approach, whose algorithm takes the reference rates as input [40]. In all three scenarios considered here the reference is known: it is constructed on the assignment’s layered graph, derived from a road network’s capacities, or estimated from trajectories as a Markov state model, so the requirement does not bind. Enumerating the state space is a second requirement. A sampler for an unnormalized target on the product space [N ]D faces 2576 states on a 24 × 24 Ising lattice. The optimal rate there keeps the potential-ratio form of (4); discrete Schrödinger-bridge samplers learn that ratio with a neural network [16].

Acknowledgments We acknowledge the use of large language models in preparing this work.

References [1] Awad H. Al-Mohy and Nicholas J. Higham. Computing the action of the matrix exponential, with an application to exponential integrators. SIAM Journal on Scientific Computing, 33(2):488–511, 2011.

10

= 250, = 1.56 = 250, = 1.91

= 2500, = 0.35 = 2500, = 0.718

= 2500, = 0.878 = 2500, = 1

(b)

1.25

10 2

residual ratio per round

residual p p( p) 2

(a)

= 250, = 0.35 = 250, = 1

10 4 10 6 10 8 10 10 0

20

40

60

round

80

1.00 0.75 0.50 0.25 0.00 0.00

100

0.25

0.50

0.75

1.00

1.25

damping over its critical value 2/(1 + L)

Figure 5: The tested dampings converge below the critical damping 2/(1 + L) and fail to converge above it. The road network of Figure 4. (a) The residual ∥p − p̄(βp)∥2 , which bounds the distance to the equilibrium (Proposition 4.1), against the round, at two congestion strengths. At damping 0.35 the curves for the two strengths start a decade apart and agree to within a factor 1.5 from round 10 on. (b) The ratio of successive residuals over the last three rounds of each run, against the damping in units of its critical value. Lines are the local rate max{1 − α(1 + m), α(1 + L) − 1} of Proposition 4.1. Triangles mark runs that did not converge: the ratio lies above one where the residual grows and at one where it cycles. At βτ = 250 the curvature is L = 0.155 (bound βτ maxx p◦ (x) = 1.64), so the critical damping is 1.732; dampings 0.35, 1 and 1.56 converge and 1.91 does not. At βτ = 2500 the curvature is L = 1.506 (bound βτ maxx p◦ (x) = 12.29), so the critical damping is 0.798; dampings 0.35 and 0.718 converge and 0.878 and 1 do not. On the branch set by L the contraction over the last three rounds is within 1.0% of α(1 + L) − 1 at both strengths. The smallest eigenvalue m of ∇2 h(p◦ ) − I is zero at both strengths, because the equilibrium occupancy is zero at the 37 congestible nodes that lie on no path from a supply node to a demand node. On the branch set by m, the local rate is therefore 1 − α. The runs at damping 0.35 contract faster than that over their last three rounds, by 2.3% at βτ = 250 and 1.0% at βτ = 2500. Above the threshold at βτ = 250 (damping 1.91) the residual grows by 1.16 per round, against α(1 + L) − 1 = 1.20.

11

[2] Francis Bach. Duality between subgradient and conditional gradient methods. SIAM Journal on Optimization, 25(1): 115–129, 2015. [3] Julio Backhoff, Giovanni Conforti, Ivan Gentil, and Christian Léonard. The mean field Schrödinger problem: ergodic behavior, entropy estimates and functional inequalities. Probability Theory and Related Fields, 178(1–2):475–530, 2020. doi: 10.1007/s00440-020-00977-8. [4] Debojjal Bagchi. Error bounds for stochastic user equilibrium traffic assignment. M.S.E. thesis, The University of Texas at Austin, 2025. [5] Debojjal Bagchi and Stephen D Boyles. Spectral analysis of the logit mapping and implications for stochastic user equilibrium algorithms. arXiv preprint arXiv:2605.21843, 2026. [6] Jean-Bernard Baillon and Roberto Cominetti. Markovian traffic equilibrium. Mathematical Programming, 111(1–2): 33–56, 2008. doi: 10.1007/s10107-006-0076-2. [7] Jean-David Benamou, Guillaume Carlier, Simone Di Marino, and Luca Nenna. An entropy minimization approach to second-order variational mean-field games. Mathematical Models and Methods in Applied Sciences, 29(08):1553–1583, 2019. doi: 10.1142/S0218202519500283. [8] Yongxin Chen, Tryphon T. Georgiou, Michele Pavon, and Allen Tannenbaum. Robust transport over networks. IEEE Transactions on Automatic Control, 62:4675–4682, 2017. doi: 10.1109/tac.2016.2626796. [9] Yongxin Chen, Tryphon T. Georgiou, and Michele Pavon. Stochastic control liaisons: Richard Sinkhorn meets Gaspard Monge on a Schrödinger bridge. SIAM Review, 63(2):249–313, 2021. doi: 10.1137/20m1339982. [10] Shui-Nee Chow, Wuchen Li, Chenchen Mou, and Haomin Zhou. Dynamical Schrödinger bridge problems on graphs. Journal of Dynamics and Differential Equations, 34(3):2511–2530, 2022. doi: 10.1007/s10884-021-09977-1. [11] Kai Cui and Heinz Koeppl. Approximately solving mean field games via entropy-regularized deep reinforcement learning. In International Conference on Artificial Intelligence and Statistics, pages 1909–1917. PMLR, 2021. [12] Asmaa Eldesoukey and Tryphon T. Georgiou. Schrödinger’s control and estimation paradigm with spatio-temporal distributions on graphs. IEEE Transactions on Automatic Control, 70(4):2466–2478, 2025. doi: 10.1109/TAC.2024. 3485537. [13] Asmaa Eldesoukey, Yongxin Chen, and Abhishek Halder. A generalized Sinkhorn algorithm for mean-field Schrödinger bridge. IEEE Control Systems Letters, 10:1087–1092, 2026. doi: 10.1109/LCSYS.2026.3704530. [14] Saeed Ghadimi. Conditional gradient type methods for composite nonlinear and stochastic optimization. Mathematical Programming, 173(1):431–464, 2019. [15] Guillaume Guex, Ilkka Kivimäki, and Marco Saerens. Randomized optimal transport on a graph: framework and new distance measures. Network Science, 7(1):88–122, 2019. doi: 10.1017/nws.2018.29. [16] Wei Guo, Yuchen Zhu, Xiaochen Du, Juno Nam, Yongxin Chen, Rafael Gómez-Bombarelli, Guan-Horng Liu, Molei Tao, and Jaemoo Choi. Discrete adjoint Schrödinger bridge sampler. arXiv preprint arXiv:2602.08243, 2026. [17] Guillaume Huguet, Alexander Tong, María Ramos Zapatero, Christopher J. Tape, Guy Wolf, and Smita Krishnaswamy. Geodesic Sinkhorn for fast and accurate optimal transport on manifolds. In IEEE International Workshop on Machine Learning for Signal Processing (MLSP), 2023. arXiv:2211.00805. [18] Shinji Iida. Molecular dynamics trajectories of protein folding. Zenodo dataset, version v0.2, 2022. [19] Hilbert J Kappen. Path integrals and symmetry breaking for optimal control theory. Journal of Statistical Mechanics: Theory and Experiment, 2005(11):P11011–P11011, 2005. doi: 10.1088/1742-5468/2005/11/p11011. [20] Anna Kazeykina, Zhenjie Ren, Xiaozhen Wang, and Yufei Zhang. Entropic optimal transport problem with convex functional cost. arXiv preprint arXiv:2503.11843, 2025.

12

[21] Philip A. Knight. The Sinkhorn–Knopp algorithm: Convergence and applications. SIAM Journal on Matrix Analysis and Applications, 30(1):261–275, 2008. doi: 10.1137/060659624. [22] Péter Kovács. Minimum-cost flow algorithms: an experimental evaluation. Optimization Methods and Software, 30(1): 94–127, 2015. [23] Pierre Lavigne and Laurent Pfeiffer. Generalized conditional gradient and learning in potential mean field games. Applied Mathematics & Optimization, 88(3):89, 2023. [24] Christian Léonard. A survey of the Schrödinger problem and some of its connections with optimal transport. Discrete & Continuous Dynamical Systems-A, 34(4):1533–1574, 2014. doi: 10.3934/dcds.2014.34.1533. [25] Guan-Horng Liu, Tianrong Chen, Oswin So, and Evangelos Theodorou. Deep generalized Schrödinger bridge. Advances in Neural Information Processing Systems, 35:9374–9388, 2022. [26] Guan-Horng Liu, Yaron Lipman, Maximilian Nickel, Brian Karrer, Evangelos A. Theodorou, and Ricky T. Q. Chen. Generalized Schrödinger bridge matching. In International Conference on Learning Representations (ICLR), 2024. [27] Tom Lollier and Eyal Neuman. Quantum-enhanced sampling of Schrödinger bridges. arXiv preprint arXiv:2609.27103, 2026. [28] M. V. Menon. Matrix links, an extremization problem, and the reduction of a non-negative matrix to one with prescribed row and column sums. Canadian Journal of Mathematics, 20:225–232, 1968. [29] Philipp Metzner, Christof Schütte, and Eric Vanden-Eijnden. Transition path theory for Markov jump processes. Multiscale Modeling & Simulation, 7(3):1192–1219, 2009. doi: 10.1137/070699500. [30] Atsushi Nitanda, Kazusato Oko, Denny Wu, Nobuhito Takenouchi, and Taiji Suzuki. Primal and dual analysis of entropic fictitious play for finite-sum problems. In International Conference on Machine Learning, pages 26266–26282. PMLR, 2023. [31] Maxence Noble, Marie Scheid, Yazid Janati, Eric Moulines, and Alain Durmus. Twisted Schrödinger bridge matching. arXiv preprint arXiv:2607.16987, 2026. [32] Frank Noé, Christof Schütte, Eric Vanden-Eijnden, Lothar Reich, and Thomas R Weikl. Constructing the equilibrium ensemble of folding pathways from short off-equilibrium simulations. Proceedings of the National Academy of Sciences, 106(45):19011–19016, 2009. doi: 10.1073/pnas.0905466106. [33] Yuki Oyama, Yusuke Hara, and Takashi Akamatsu. Markovian traffic equilibrium assignment based on network generalized extreme value model. Transportation Research Part B: Methodological, 155:135–159, 2022. [34] Axel Ringh, Isabel Haasler, Yongxin Chen, and Johan Karlsson. Graph-structured tensor optimization for nonlinear density control and mean field games. SIAM Journal on Control and Optimization, 62(4):2176–2202, 2024. [35] Yosef Sheffi. Urban Transportation Networks: Equilibrium Analysis with Mathematical Programming Methods. PrenticeHall, Englewood Cliffs, NJ, 1985. ISBN 0-13-939729-9. [36] Justin Solomon, Fernando de Goes, Gabriel Peyré, Marco Cuturi, Adrian Butscher, Andy Nguyen, Tao Du, and Leonidas Guibas. Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains. ACM Transactions on Graphics, 34(4):66:1–66:11, 2015. doi: 10.1145/2766963. [37] Sharan Srinivasan, Berke M. Turkay, and Harsha Honnappa. The BAR-SOT method: Long-term average cost control as stochastic optimal self-transport. arXiv preprint arXiv:2609.17966, 2026. [38] Takashi Tanaka, Ehsan Nekouei, Ali Reza Pedram, and Karl Henrik Johansson. Linearly solvable mean-field traffic routing games. IEEE Transactions on Automatic Control, 66(2):880–887, 2021. [39] Alexis M. H. Teter, Wenqing Wang, and Abhishek Halder. Schrödinger bridge with quadratic state cost is exactly solvable. IEEE Transactions on Automatic Control, 71(5):2903–2917, 2026. doi: 10.1109/TAC.2025.3631521.

13

[40] Panagiotis Theodoropoulos, Juno Nam, Evangelos Theodorou, and Jaemoo Choi. Generalized Schrödinger bridge on graphs. In International Conference on Machine Learning (ICML), 2026. [41] Alexis Thibault, Lénaïc Chizat, Charles Dossal, and Nicolas Papadakis. Overrelaxed Sinkhorn–Knopp algorithm for regularized optimal transport. Algorithms, 14(5):143, 2021. [42] Emanuel Todorov. Linearly-solvable Markov decision problems. In Advances in Neural Information Processing Systems, volume 19, 2006. [43] Yinuo Zhang, Sophia Tang, and Pranam Chatterjee. ScooBDoob: Schrödinger bridge with Doob’s h-transform for molecular dynamics. In NeurIPS 2025 Workshop on Structured Probabilistic Inference & Generative Modeling, 2025.

14

Appendices to “Cost-augmented Schrödinger bridges on graphs are exactly solvable: a Feynman–Kac tilt replaces learned control”

A

Proofs

A.1

Proofs for Section 2

e differ by a constant Proof of Theorem 2.1. The generalized-bridge objective, which is the left-hand side of (2), and KL(q∥R) Rτ e independent of q: expanding the latter, KL(q∥R) = Eq [log q − log R + 0 f + log Z], which is the left-hand side plus e is the normalized path measure of etA : By the Feynman–Kac representation, R log Z. The minimizers therefore coincide. R − f multiplying the path measure by e subtracts f from the diagonal, leaving the off-diagonal rates unchanged. e whose h-transform representation gives Proof of Corollary 2.2. By Theorem 2.1 the optimizer is the plain bridge against R, the controlled generator Axy φt (y)/φt (x) off the diagonal, where A is the tilted operator (3). Since A = Q − diag(f ) agrees with Q off the diagonal, Axy = Qxy for x ̸= y, which is (4). The support claim is immediate: φt (x) > 0 wherever qt⋆ (x) > 0, so the ratio is finite on the support of qt⋆ . Moreover u⋆t (x, y) = 0 whenever Qxy = 0.

A.2

Proofs for Section 4

Proof R τ of Proposition 4.1. The bridge with cost g has path law dqg /dR ∝ exp{a(X0 ) + b(Xτ ) − ⟨g, ℓ⟩}, where ℓ(x) = 1{Xt = x} dt is the occupation time of x ∈ C. Here a and b are the logarithms of its endpoint scalings on the supports 0 of µ and ν, fixed by the two marginal constraints up to adding a constant to a and subtracting it from b. By Lagrangian duality Φ(g) = maxa,b {⟨µ, a⟩+⟨ν, b⟩−log ER ea(X0 )+b(Xτ )−⟨g,ℓ⟩ }, attained at those a, b, so ∇Φ(g) = Eqg [ℓ] = τ p̄(g). The scalings depend smoothly on g because the endpoint block is positive, so Φ is twice continuously differentiable. Hence ∇h(p) = p − p̄(βp), so the step is p − α∇h(p). Differentiating Eqg [ℓ] while a and
b move to hold the marginals fixed gives, for v ∈ RC , −v ⊤ ∇2 Φ(g) v = minu0 ,uτ Varqg ⟨v, ℓ⟩ − u0 (X0 ) − uτ (Xτ ) , the variance of ⟨v, ℓ⟩ left after regressing it on additive functions of the endpoints. This is nonnegative and at most Varqg ⟨v, ℓ⟩, with v extended by zero off C. With U uniform on [0, τ ] and independent of the path, ⟨v, ℓ⟩ = τ E[v(XU ) | X]. The law of XU is the time-averaged marginal of qg , whose restriction to C is p̄(g).
Since v vanishes off C, the law of total variance gives Varqg ⟨v, ℓ⟩ ≤ τ 2 Var v(XU ) = τ 2 v ⊤ diag p̄(g) − p̄(g)p̄(g)⊤ v. The last form is at most τ 2 maxx p̄(g)(x) ∥v∥22 . It is 2 also at most τ2 ∥v∥22 , because v(XU ) takes its values among those of v and 0, which for a unit vector lie in an interval of √ length at most 2; a variable confined to such an interval has variance at most 1/2. So ∇2 h(p) = I − βτ ∇2 Φ(βp) lies between I and (1 + βτ /2)I. Moreover ∇h(p◦ ) = 0, because the minimizer of a convex objective is the best response to its own linearization. Since h is strongly convex, it has no other critical point, so p◦ is its unique minimizer and a fixed point
R1 of the step. For any p, p′ the two steps differ by M (p − p′ ), where M = 0 I − α∇2 h(p′ + s(p − p′ )) ds is symmetric with spectrum in [1 − α(1 + βτ /2), 1 − α], so the step contracts by ρ. With p′ = p◦ , M tends to I − α∇2 h(p◦ ) as p → p◦ because ∇2 h is continuous. That limit is symmetric with spectrum in [1 − α(1 + L), 1 − α(1 + m)], so its norm is ρ◦ . Hence ∥M ∥ ≤ ρ◦ + ϵ < 1 for every p in a small enough ball around p◦ , so the step maps the ball into itself. Strong convexity gives ∥p − p◦ ∥2 ≤ ∥∇h(p) − ∇h(p◦ )∥2 = ∥∇h(p)∥2 .

15

Related documents

• Cost-augmented Schrödinger bridges on graphs are exactly solvable: a Feynman-Kac tilt replaces learned control#1157543
• Learning Latent Graph Geometry via Fixed-Point Schrödinger-Type Activation: A Theoretical Study#984058
• Spectral Prioritized Sweeping in Nonstationary Reinforcement Learning#973786
• Reinforcement Learning Assisted Quantum Simulation of Many-Body Excited States and Real-Time Dynamics#1084215
• ΨLM: Coupling Frozen Language and Physics Models through Trainable Latent Bridges#687435