Paper deep dive
A Generalized Sinkhorn Algorithm for Mean-Field Schrödinger Bridge
Asmaa Eldesoukey, Yongxin Chen, Abhishek Halder
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 93%
Last extracted: 4/10/2026, 1:59:38 AM
Summary
The paper introduces a generalized Sinkhorn algorithm to solve the mean-field Schrödinger bridge (MFSB) problem, which involves controlling a diffusion process with nonlocal interactions. By proposing a generalized Hopf-Cole transform, the authors convert the nonconvex MFSB problem into a system of coupled integro-PDEs, enabling a recursive numerical solution with convergence guarantees.
Entities (5)
Relation Signals (3)
Sinkhorn algorithm → solves → Mean-Field Schrödinger Bridge
confidence 95% · design a Sinkhorn-type recursive algorithm to solve the associated system of integro-PDEs
Hopf-Cole transform → enables → Sinkhorn algorithm
confidence 90% · We propose a generalization of the Hopf-Cole transform for MFSB and, building on it, design a Sinkhorn-type recursive algorithm
Mean-Field Schrödinger Bridge → modeledby → McKean-Vlasov integro-PDE
confidence 90% · the dynamics of pt is described by the McKean-Vlasov integro-PDE
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:The mean-field Schrödinger bridge (MFSB) problem concerns designing a minimum-effort controller that guides a diffusion process with nonlocal interaction to reach a given distribution from another by a fixed deadline. Unlike the standard Schrödinger bridge, the dynamical constraint for MFSB is the mean-field limit of a population of interacting agents with controls. It serves as a natural model for large-scale multi-agent systems. The MFSB is computationally challenging because the nonlocal interaction makes the problem nonconvex. We propose a generalization of the Hopf-Cole transform for MFSB and, building on it, design a Sinkhorn-type recursive algorithm to solve the associated system of integro-PDEs. Under mild assumptions on the interaction potential, we discuss convergence guarantees for the proposed algorithm. We present numerical examples with repulsive and attractive interactions to illustrate the theoretical contributions.
Tags
Links
- Source: https://arxiv.org/abs/2604.06531v2
- Canonical: https://arxiv.org/abs/2604.06531v2
Trouble viewing inline? Open PDF directly →
Full Text
57,049 characters extracted from source content.
Expand or collapse full text
A Generalized Sinkhorn Algorithm for Mean-Field Schrödinger Bridge Asmaa Eldesoukey, Yongxin Chen, Abhishek Halder Asmaa Eldesoukey and Abhishek Halder are with the Department of Aerospace Engineering, Iowa State University, Ames, IA 50011, USA, asmaae,ahalder@iastate.edu.Yongxin Chen is with the School of Aerospace Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA, ychen3148@gatech.edu.This research was partially supported by NSF awards 2111688, 2450377, 2450378. Abstract The mean-field Schrödinger bridge (MFSB) problem concerns designing a minimum-effort controller that guides a diffusion process with nonlocal interaction to reach a given distribution from another by a fixed deadline. Unlike the standard Schrödinger bridge, the dynamical constraint for MFSB is the mean-field limit of a population of interacting agents with controls. It serves as a natural model for large-scale multi-agent systems. The MFSB is computationally challenging because the nonlocal interaction makes the problem nonconvex. We propose a generalization of the Hopf-Cole transform for MFSB and, building on it, design a Sinkhorn-type recursive algorithm to solve the associated system of integro-PDEs. Under mild assumptions on the interaction potential, we discuss convergence guarantees for the proposed algorithm. We present numerical examples with repulsive and attractive interactions to illustrate the theoretical contributions. I Introduction The problem. For some large positive integer N, consider a collection of N interacting agents with noisy controlled dynamics: dXti=σati−1N∑j=1N∇W(Xti−Xtj)dt+σdBti, X_t^i= \\!σ a_t^i- 1N _j=1^N∇ W\! (X_t^i-X_t^j )\!\! \\> t+σ\> B_t^i, (1) where i∈⟦N⟧:=1,…,Ni∈ N :=\1,…,N\, time t∈[0,1]t∈[0,1], and Xti,ati,Bti∈ℝdX_t^i,a_t^i,B_t^i ^d denote the state, control action and (standard Brownian) process noise for the iith agent, respectively. In (1), the constant σ>0σ>0 is the input/noise strength, ∇ denotes the spatial gradient operator, and W(⋅)W(·) is a known interaction potential satisfying the standard assumptions A1-A3: A1 W∈2(ℝd;ℝ)W∈ C^2(R^d;R), A2 W is symmetric, i.e., W(x)=W(−x)∀x∈ℝdW(x)=W(-x)\>∀ x ^d, A3 the Hessian Hess(W)Hess(W) is uniformly upper bounded. In the large population limit (N→∞N→∞), we are interested in designing a minimum-effort feedback control strategy ati=:ut(Xti),i∈⟦N⟧, a_t^i=:u_t (X_t^i ), i∈ N , that guides the initial stochastic state X0i∼pinX_0^i p_in to a final stochastic state X1i∼pfinX_1^i p_fin over the fixed time horizon [0,1][0,1]. Here, pin,pfinp_in,p_fin are the given initial and final probability density functions (PDFs) with finite second moments. The minimum effort objective is understood as that of minimizing the average quadratic control action ∫01∑i=1N|ati|2dtE _0^1 _i=1^N|a_t^i|^2 t where E denotes the expectation operator induced by the underlying probability distribution, and |⋅||·| denotes the Euclidean norm. We refer to N→∞N→∞ as the mean-field limit since the collective dynamics then are governed by the PDF pt≈1N∑i=1NδXti,whereδXtidenotes Dirac delta atXti.p_t≈ 1N _i=1^N _X_t^i,\;where\; _X_t^i\;denotes Dirac delta at\;X_t^\i\. In that limit, −1N∑j=1N∇W(Xti−Xtj)≈−∇W∗pt- 1N _j=1^N∇ W(X_t^i-X_t^j)≈-∇ W*p_t, where ∗* denotes convolution, and the dynamics of ptp_t is described by the McKean-Vlasov integro-PDE [1]: ∂tpt(x)+∇⋅(pt(x)(σut(x)−(∇W∗pt)(x)))=σ22Δpt(x), _tp_t(x)\!+\!∇\!·\! (p_t(x) (σ u_t(x)\!-\!(∇ W*p_t)(x)\! )\! )\!\!=\! σ^22 p_t(x), where ∇⋅∇· denotes the divergence, and Δ is the standard Laplacian operator. Notice that the empirical measure ∑i=1NδXti/N _i=1^N _X_t^i/N is random for finite N. However, in the mean field limit, ptp_t is a deterministic PDF for any t∈[0,1]t∈[0,1]. Thus, our stochastic optimal control problem of interest is the following. Given two PDFs pin,pfinp_in,p_fin supported over subsets of ℝdR^d with finite second moments, determine finite-energy control policy (ut)t∈[0,1](u_t)_t∈[0,1] and the PDF-valued 1 C^1-in-time curve (pt)t∈[0,1](p_t)_t∈[0,1] that solve the variational problem: minimizeut,pt12∫01∫ℝd|ut(x)|2pt(x)dxdt, u_t,p_tminimize 12 _0^1 _R^d|u_t(x)|^2\>p_t(x)\> x\> t, (2a) ∂tpt(x)+∇⋅(pt(x)(σut(x)−(∇W∗pt)(x)))=σ22Δpt(x), _tp_t(x)\!+\!∇\!·\! (p_t(x) (σ u_t(x)\!-\!(∇ W*p_t)(x)\! )\! )\!\!=\! σ^22 p_t(x), (2b) p0=pin,p1=pfin. p_0=p_in, p_1=p_fin. (2c) We refer to (2) as the mean-field Schrödinger bridge (MFSB) problem. Related works. The MFSB differs from the classical Schrödinger bridge (SB) in that the latter is the W≡0W≡ 0 special case of (2). Note in particular that for W≠0W≠ 0, the dynamical constraint (2b) is nonlinear in ptp_t. In the absence of interaction, several control-theoretic works have studied the SB with more general drift and diffusion coefficients [2, 3, 4, 5, 6, 7], with additional state cost [8, 9, 10, 11] and constraints [12, 13, 14, 15]. In comparison, problem (2) is much less explored. In the probability literature, [16] established the large deviation principle for the MFSB problem (2) under the stated assumptions A1-A3 and proved the existence of solutions. That work does not address numerical solutions. In the control literature, [17] proposed to numerically solve (2) by first transforming the dynamic problem to the static large deviations problem, then solving the discrete version of the same as a non-standard multi-marginal SB. The reformulated problem therein remains nonconvex, and [17] solved it using a proximal gradient algorithm that is guaranteed to converge locally at a sublinear rate. Motivation. By suitably modeling or interpreting the interaction potential W, problem (2) has natural applications in steering multi-agent populations as in crowd dynamics [18], population biology [19, 20], and in swarm robotics [21]. For instance, W may model limited communication or information sharing due to trust, privacy, or power constraints. Contributions. Our technical contributions are threefold. • We propose a generalized Hopf-Cole transform for (2) to derive a new Schrödinger system (Sec. I-A). • Using the above, we design a generalized Sinkhorn algorithm (Sec. I-B) to solve (2) numerically. Beyond novelty, the proposed algorithm is a significant structural generalization of existing Sinkhorn algorithms. • We prove local convergence guarantee (Sec. I-B) for the proposed algorithm, and illustrate its effectiveness via numerical examples (Sec. IV). I Conditions for Optimality For computational convenience, we scale the interaction potential by the input/noise strength σ as W↦σ2W σ^2W without loss of generality. Thereby, we consider the following scaled version of (2b): ∂tpt+∇⋅(pt(σut−σ2(∇W∗pt)))=σ22Δpt. _tp_t\!+\!∇\!·\! (p_t (σ u_t\!-\!σ^2(∇ W*p_t)\! )\! )\!\!=\! σ^22 p_t. (2b′) We derive the first-order necessary conditions of optimality for the MFSB via the equivalent saddle point problem: minut,ptmaxψt:=∫01∫ℝd12|ut(x)|2pt(x)+ψt(x)[∂tpt(x)+ _u_t,p_t\; _ _tJ\!\!:=\!\ _0^1\!\!\! _R^d \ 12|u_t(x)|^2p_t(x)+ _t(x) [ _tp_t(x)+ ∇⋅(pt(x)(σut(x)−σ2(∇W∗pt)(x)))−σ22Δpt(x)]dxdt, \!∇\!·\! (p_t(x) (σ u_t(x)\!-\!σ^2(∇ W*p_t)(x)\! )\! )\!\!-\! σ^22 p_t(x) ] \\, dx\, dt, subject to (2c), where the Lagrange multiplier ψt(x)∈1,2([0,1]×ℝd;ℝ) _t(x)∈ C^1,2([0,1]×R^d;R). We use integration by parts, assuming a sufficiently fast decay at infinity so that boundary terms vanish111Boundary terms also vanish in the case of compactly-supported densities., to express J as =∫01∫ℝdpt(x)12|ut(x)|2−[∂tψt(x)+σ⟨ut(x),∇ψt(x)⟩ = _0^1\!\!\! _R^dp_t(x) \ 12|u_t(x)|^2- [ _t _t(x)+σ u_t(x),∇ _t(x) +σ22Δψt(x)]dxdt+σ2∫01∫ℝd×ℝdSt(x,y)dxdydt + σ^22 _t(x) ] \\, dx\, dt+σ^2 _0^1\!\! _R^d×R^dS_t(x,y)\, dx\, dy\, dt +∫ℝd(p1(x)ψ1(x)−p0(x)ψ0(x))dx, + _R^d (p_1(x) _1(x)-p_0(x) _0(x) )\, dx, (3) with222Throughout, we use ⟨⋅,⋅⟩ ·,· to denote the Euclidean inner product. St(x,y):=12⟨∇ψt(x)−∇ψt(y),∇W(x−y)⟩pt(x)pt(y).S_t(x,y):= 12 ∇ _t(x)-∇ _t(y),∇ W(x-y) \,p_t(x)p_t(y). In obtaining (3), we used the symmetry of W, namely W(−x)=W(x)W(-x)=W(x) together with its implication ∇W(−x)=−∇W(x).∇ W(-x)=-∇ W(x). As a result, it holds that ∫ℝd×ℝd⟨∇ψt(x),∇W(x−y)⟩pt(x)pt(y)dxdy _R^d×R^d ∇ _t(x),∇ W(x-y) p_t(x)\,p_t(y)\> dx\> dy =−∫ℝd×ℝd⟨∇ψt(y),∇W(x−y)⟩pt(x)pt(y)dxdy =- _R^d×R^d ∇ _t(y),\!∇ W(x-y) p_t(x)p_t(y)\, dx\, dy =∫ℝd×ℝdSt(x,y)dxdy. = _R^d×R^dS_t(x,y)\, dx\, dy. Minimizing (3) over utu_t gives the optimal control utopt(x)=σ∇ψt(x). u_t opt(x)=σ∇ _t(x). (4) Substituting (4) back into J and minimizing over ptp_t yields a forward Hamilton-Jacobi-Bellman (HJB) equation ∂tψt(x)+σ22|∇ψt(x)|2+σ22Δψt(x)= _t _t(x)+ σ^22|∇ _t(x)|^2+ σ^22 _t(x)= σ2∫ℝd⟨∇ψt(x)−∇ψt(y),∇W(x−y)⟩ptopt(y)dy. σ^2\! _R^d\!\! ∇ _t(x)-∇ _t(y),∇ W(x-y) p_t opt(y)\, dy. (5) By minimizing over ψt _t, we recover (2b′), whereas the last term in (3) is fixed by the constraints (2c). That is, (4), (I), (2b′), and (2c) form the first-order necessary conditions for the MFSB, or equivalently, for the saddle-point problem. It turns out [16] that the optimal PDFs ptoptp_t^opt solving Problem (2) factor into products that are determined by a forward-backward HJB system of coupled nonlinear integro-PDEs. This is summarized in Proposition 1 next. We outline its proof for completeness since it was omitted in [16]. Proposition 1. [16, Corollary 1.2] The solution of Problem (2) can be characterized by ptopt(x)=exp(−2(W∗ptopt)(x)+ψt(x)+ϕt(x)), p_t opt(x)= (-2(W*p_t opt)(x)+ _t(x)+ _t(x)), (6) where ψt _t satisfies the forward HJB equation in (I), and ϕt _t satisfies the backward HJB equation −∂tϕt(x)+σ22|∇ϕt(x)|2+σ22Δϕt(x)= - _t _t(x)+ σ^22|∇ _t(x)|^2+ σ^22 _t(x)= σ2∫ℝd⟨∇ϕt(x)−∇ϕt(y),∇W(x−y)⟩ptopt(y)dy. σ^2\! _R^d\!\! ∇ _t(x)-∇ _t(y),∇ W(x-y) p_t opt(y)\, dy. (7) Equations (I), (6), (1) are subject to the boundary conditions pin(x) p_ in(x) =exp(−2(W∗pin)(x)+ψ0(x)+ϕ0(x)), = (-2(W*p_ in)(x)+ _0(x)+ _0(x)), (8a) pfin(x) p_ fin(x) =exp(−2(W∗pfin)(x)+ψ1(x)+ϕ1(x)). = (-2(W*p_ fin)(x)+ _1(x)+ _1(x)). (8b) Proof. We let btopt(x):=−(∇W∗ptopt)(x)b_t opt(x):=-(∇ W*p_t opt)(x). Then, from (6), −∂tϕt - _t _t =−∂t(logptopt)−2(W∗∂tptopt)+∂tψt, =- _t( p_t opt)-2(W* _tp_t opt)+ _t _t, (9a) σ22|∇ϕt|2 σ^22|∇ _t|^2 =σ22|∇(logptopt)−2btopt−∇ψt|2, = σ^22|∇( p_t opt)-2b_t opt-∇ _t|^2, (9b) σ22Δϕt σ^22 _t =σ22(Δ(logptopt)−2(∇⋅btopt)−Δψt). = σ^22 ( ( p_t opt)-2(∇· b_t opt)- _t ). (9c) By the identity |∇(logptopt)|2+Δ(logptopt)=Δptoptptopt, |∇( p_t opt)|^2+ ( p_t opt)= p_t optp_t opt, (10) the evolution of logptopt p_t opt is governed by ∂t(logptopt)=σ22|∇(logptopt)|2+σ22Δ(logptopt) _t( p_t opt)= σ^22|∇( p_t opt)|^2+ σ^22 ( p_t opt) −σ2(Δψt+∇⋅btopt+⟨∇ψt+btopt,∇(logptopt)⟩). -σ^2 ( _t+∇· b_t opt+ ∇ _t\!+\!b_t opt,∇( p_t opt) ). (11) By substituting (I) into (9a), summing the left-hand sides of (9a)-(9c), and expanding the quadratic term on the right-hand side of (9b), we obtain −∂tϕt(x) - _t _t(x) +σ22|∇ϕt(x)|2+σ22Δϕt(x)= + σ^22|∇ _t(x)|^2+ σ^22 _t(x)= +2W∗(∇⋅(σ2ptopt(∇ψt+btopt)))(x) +2W* (∇·(σ^2p_t opt(∇ _t+b_t opt)) )(x) −σ2⟨∇ϕt(x),btopt(x)⟩−σ2W∗Δptopt(x) -σ^2 ∇ _t(x),b_t opt(x) -σ^2W* p_t opt(x) −σ2∫ℝd⟨∇ψt(y),∇W(x−y)⟩ptopt(y)dy. -\!σ^2\!\! _R^d\! ∇ _t(y),∇ W(x-y) \,p_t opt(y)\, dy. (12) Utilizing the symmetry of W together with integration by parts, we have that W∗(∇⋅(σ2ptopt(∇ψt+btopt)))(x) W* (∇·(σ^2p_t opt(∇ _t+b_t opt)) )(x) = = ∫ℝd⟨∇W(x−y),(σ2ptopt(∇ψt+btopt))(y)⟩dy. _R^d ∇ W(x-y), (σ^2p_t opt (∇ _t+b_t opt ) )(y) \> dy. In addition, W∗Δptopt(x)=∫ℝd⟨∇W(x−y),∇ptopt(y)⟩dy.W* p_t opt(x)= _R^d ∇ W(x-y),∇ p_t opt(y) \, dy. Substituting the previous two equations in (I) yields (1). Combining (2c) and (6) gives (8). ∎ I Schrödinger System and its Solution I-A Hopf-Cole Transform Our idea now is to transcribe the coupled system of equations (I), (6), (1) and (8) into a Schrödinger system, for which we will design a generalized Sinkhorn recursion enabling numerical computation. To do so, we use the Hopf-Cole transform [22, 23] (ψt,ϕt)↦(φt,φ^t):=(expψt,expϕt). ( _t, _t) ( _t, _t):=( _t, _t). (13) Combining (6) and (13), we get ptopt(x)=e−2(W∗ptopt)(x)φt(x)φ^t(x). p_t opt(x)=e^-2(W*p_t opt)(x) _t(x) _t(x). (14a) Furthermore, ∂tφt(x) _t _t(x) +σ22Δφt(x)=σ2φt(x)∫ℝd⟨∇logφt(x) + σ^22 _t(x)=σ^2 _t(x) _R^d\!\! \ ∇ _t(x) −∇logφt(y),∇W(x−y)⟩ptopt(y)dy, -∇ _t(y),∇ W(x-y) \,p_t opt(y) \\, dy, where we have used the identity akin to (10). If we let Qt(ξt,ptopt)(x):=−∫ℝd⟨∇logξt(y),∇W(x−y)⟩ptopt(y)dy, \!\!\!\!\!Q_t( _t,p_t opt)(x)\!:=-\!\! _R^d\!\! ∇ _t(y),\!∇ W(x-y) \,p_t opt(y) dy, then φt _t satisfies ∂tφt+σ2⟨btopt,∇φt⟩+σ22Δφt=σ2Qt(φt,ptopt)φt, _t _t\!+\!σ^2 b_t opt,∇ _t \!+\! σ^22 _t\!=\!σ^2Q_t( _t,p_t opt) _t, (14b) with btopt:=−(∇W∗ptopt)b_t opt:=-(∇ W*p_t opt). Analogously, φ^t _t satisfies ∂tφ^t−σ2⟨btopt,∇φ^t⟩−σ22Δφ^t=−σ2Qt(φ^t,ptopt)φ^t. \!\! _t _t\!-\!σ^2 b_t opt,∇ _t \!-\! σ^22 _t\!=\!-σ^2Q_t( _t,p_t opt) _t. (14c) Equations (14a)-(14c) are subject to the boundary conditions e−2(W∗pin)(x)φ0(x)φ^0(x)=pin(x), e^-2(W*p_ in)(x) _0(x) _0(x)=p_ in(x), (14d) e−2(W∗pfin)(x)φ1(x)φ^1(x)=pfin(x). e^-2(W*p_ fin)(x) _1(x) _1(x)=p_ fin(x). (14e) The system of equations (14) constitutes a nonlinear coupled-through-time Schrödinger system with mean-field interactions. As expected, for W=0W=0, (14) reduces to the classical Schrödinger system. To numerically solve (14), we next propose a generalized Sinkhorn algorithm that simultaneously evades the nonlinearities in the drift and source/sink (reaction) terms in (14b)-(14c). Algorithm 1 Mean-Field Sinkhorn 0: pinp_ in, pfinp_ fin, W σ, p(0),p^(0), such that p0(0)=pin,p1(0)=pfinp_0^(0)=p_ in,p_1^(0)=p_ fin, φ(0),φ^(0) ^(0), ^(0), N1N_1, N2,N3N_2,N_3, tol, θ∈(0,1]θ∈(0,1]. 0: p,φ,φ^p, , . 1: for k=0,1,⋯,N1k=0,1,·s,N_1 do 2: Compute b(k)=−∇W∗p(k)b^(k)=-∇ W*p^(k) 3: Set φj←φ(k) ^j← ^(k) and φ^j←φ^(k) ^j← ^(k) 4: for j=0,1,⋯,N2j=0,1,\,·s,N_2 do 5: Set φ10←φ1j _1^\0\← _1^j and φ^00←φ^0j _0^\0\← _0^j 6: for i=0,1,⋯,N3i=0,1,\,·s,N_3 do 7: Integrate backward in time ∂tφti+1 _t _t^\i+1\ +σ2⟨bt(k),∇φti+1⟩+σ22Δφti+1 +σ^2 b_t^(k),∇ _t^\i+1\ + σ^22 _t^\i+1\ =σ2Qt(φtj,pt(k))φti+1 =σ^2Q_t ( _t^j,p_t^(k) ) _t^\i+1\ (15a) with φt=1i+1(x)=φ1i(x) _t=1^\i+1\(x)= _1^\i\(x) 8: Compute φ^0i+1(x)=pin(x)φ0i+1(x)e2(W∗pin)(x) _0^\i+1\(x)= p_ in(x) _0^\i+1\(x)e^2(W*p_ in)(x) (15b) 9: Integrate forward in time ∂tφ^ti+1 _t _t^\i+1\ −σ2⟨bt(k),∇φ^ti+1⟩−σ22Δφ^ti+1 -σ^2 b_t^(k),∇ _t^\i+1\ - σ^22 _t^\i+1\ =−σ2Qt(φ^tj,pt(k))φ^ti+1 =-σ^2Q_t ( _t^j,p_t^(k) ) _t^\i+1\ (15c) with φ^t=0i+1(x)=φ^0i+1(x) _t=0^\i+1\(x)= _0^\i+1\(x) 10: Compute φ1i+1(x)=pfin(x)φ^1i+1(x)e2(W∗pfin)(x) _1^\i+1\(x)= p_ fin(x) _1^\i+1\(x)e^2(W*p_ fin)(x) (15d) 11: if maxdH(φ1i+1,φ1i),dH(φ^0i+1,φ^0i)<tol \d_H( _1^\i+1\, _1^\i\),d_H( _0^\i+1\, _0^\i\)\< tol then 12: break 13: end if 14: end for 15: Set φj+1←φi+1 ^j+1← ^\i+1\ and φ^j+1←φ^i+1 ^j+1← ^\i+1\ 16: if dH⊕((φj+1,φ^j+1),(φj,φ^j))<told_H (( ^j+1, ^j+1),( ^j, ^j) )\!<\! tol then 17: break 18: end if 19: end for 20: φ(k+1)←φj+1 ^(k+1)← ^j+1 and φ^(k+1)←φ^j+1 ^(k+1)← ^j+1 21: Compute (p(k)=e−2(W∗p(k))φ(k+1)φ^(k+1) (p^(k)=e^-2(W*p^(k)) ^(k+1) ^(k+1) (16) 22: (p(k))←(p(k))/∫(p(k))xC(p^(k)) (p^(k))/ (p^(k))\,dx p(k+1)=θ(p(k))+(1−θ)p(k) \!\!\!\!\!\!\!p^(k+1)= (p^(k))\!+\!(1-θ)p^(k) (17) 23: if dH(p(k+1),p(k))<told_H^S(p^(k+1),p^(k))< tol then 24: break 25: end if 26: end for I-B Generalized Sinkhorn Algorithm We introduce Algorithm 1 that consists of nested fixed-point iterations. Inputs to this algorithm comprise of the problem data pin,pfin,W,σp_in,p_fin,W,σ, the initial guess (p(0),φ(0),φ^(0)) (p^(0), ^(0), ^(0) ), 333Throughout, we adopt notation such as p=(pt)t∈[0,1]p=(p_t)_t∈[0,1] to succinctly represent trajectories/curves over time.and internal parameters, viz. number of iterations N1,N2,N3N_1,N_2,N_3 for the for loops, numerical tolerance tol, and a damping constant θ∈(0,1]θ∈(0,1]. The initial guess (p(0),φ(0),φ^(0))(p^(0), ^(0), ^(0)) of Algorithm 1 can be taken as the solution of the classical (non-interacting) SB problem. The innermost Sinkhorn-type iteration, with index i, solves (14b)-(14e) with fixed ptp_t and reaction terms QtQ_t. As a result, this step consists of time-marching the solution of the backward-forward Kolmogorov equations with reaction terms. Upon completion of the innermost iteration, we obtain trajectories (φ,φ^)( , ) that are used to update the reaction terms and the boundary condition φ1i _1^\i\. The outer iteration in Algorithm 1 updates ptp_t using a damped version of (14a) with a fixed damping parameter θ∈(0,1]θ∈(0,1]. As is well-known [24, 25], damped fixed-point iterates improve numerical stability by mitigating aggressive updates. We defer explaining the termination criteria involving the distances dH(⋅,⋅),dH⊕(⋅,⋅)d_H(·,·),d_H (·,·) and dH(⋅,⋅)d_H^S(·,·) to Section I-C. I-C Convergence Guarantees We now present the convergence guarantees for the proposed Sinkhorn iteration. We consider a compact set Ω⊂ℝd ^d, and the Banach space L∞(Ω)L^∞( ) with norm ∥⋅∥∞:=supx∈Ω|⋅|\|·\|_∞:= _x∈ |·|. 444Throughout, sup or inf are understood as esssup ess and essinf ess .Let :=f∈L∞(Ω):f≥0 a.e.K:=\f∈ L^∞( ):f≥ 0 a.e.\ be the cone of non-negative functions in L∞(Ω)L^∞( ), and let 0K_0 be the interior of K. The Hilbert metric dHd_H [26, 27, 28] between f,g∈0f,g _0 is dH(f,g) \!\!\!\!d_H(f,g) :=logsupx∈Ωf(x)g(x)−loginfx∈Ωf(x)g(x). := _x∈ f(x)g(x)- _x∈ f(x)g(x). (18) The Hilbert metric is projective (pseudo-metric on 0K_0) because dH(f,g)=0d_H(f,g)=0 when f=cgf=cg for any c≥0c≥ 0. We say a positive map :0→0C:K_0 _0 is locally contractive with respect to the Hilbert metric near f∈0f _0 if there exists a neighborhood D⊂0D _0 containing f and a constant λ∈(0,1)λ∈(0,1) such that dH((g),(h))≤λdH(g,h)∀g,h∈D. d_H(C(g),C(h))≤λ\,d_H(g,h) ∀ g,h∈ D. Whenever the previous inequality holds for 0K_0, the map C is globally contractive. Notably, Chen et al. [28] established that a Sinkhorn iteration akin to the innermost iteration in Algorithm 1 is globally contractive in dHd_H, and therefore, by the Banach contraction mapping theorem, converges to a unique fixed point. We rely on this result to conclude that, the innermost iteration in Algorithm 1 also converges to a unique fixed point. Thus, it remains to establish local convergence of the iterations indexed by j,kj,k, which entail the reaction-term update and the fixed-point iteration in p. To this end, we consider the space of positive trajectories =f=(ft)t∈[0,1]:ft∈0∀t∈[0,1],A=\f=(f_t)_t∈[0,1]:f_t _0\ ∀ t∈[0,1]\, and the space of trajectories of probability densities :=h=(ht)t∈[0,1]:∫Ωht(x)=1∀t∈[0,1],S:=\h=(h_t)_t∈[0,1]: _ h_t(x)=1\ ∀ t∈[0,1]\, with ⊂S . We let dH(p,q):=supt∈[0,1]dH(pt,qt)∀p,q∈, d_H^S(p,q):= _t∈[0,1]d_H(p_t,q_t) ∀ p,q , (19) defining a distance between trajectories of probability densities. Let =(φ,φ^),=(μ,μ^)∈× =( , ), μ=(μ, μ) ×A, then we define dH⊕(,):=maxsupt∈[0,1]dH(φt,μt),supt∈[0,1]dH(φ^t,μ^t), \!\!\!d_H ( , μ)\!:=\! \\!\! _t∈[0,1]\!d_H( _t, _t), _t∈[0,1]\!d_H( _t, μ_t)\! \, (20) to measure distances between pairs of positive trajectories. We construct the map :→C:S as (p) (p) :=e−2(W∗p)φ(p)φ^(p), :=e^-2(W*p)\, (p)\, (p), (21) where φ(p),φ^(p) (p), (p) are the limits of the intermediate iterations for a fixed p. For a fixed p∈p , we define :×→×G:A×A ×A, for all j by p(j(p))=j+1(p)j(p)=(φj(p),φ^j(p)). _p( ^j(p))= ^j+1(p) ^j(p)=( ^j(p), ^j(p)). (22) In what follows, we establish the local convergence of the iterations induced by C (Theorem 1) with respect to dH(⋅,⋅)d_H^S(·,·) and by pG_p with respect to dH⊕(⋅,⋅)d_H (·,·) (Theorem 2). Together, Theorems 1 and 2 establish the local convergence for Algorithm 1 (Theorem 3). Theorem 1. Let p∗p^* be a fixed point of the map C defined in (21). For r>0r>0, consider a ball of radius r centered at p∗p^*, given by Br(p∗):=p∈:dH(p,p∗)≤r.B_r(p^*):=\p :d_H^S(p,p^*)≤ r\. In addition to A1-A3, we make the following assumptions B1-B5. B1 W,∇W,ΔW∈L∞(Ω)W,∇ W, W∈ L^∞( ). B2 suptsupp∈Br(p∗)‖∇ψt(p)‖∞≤a1 _t _p∈ B_r(p^*)\|∇ _t(p)\|_∞≤ a_1 for some a1<∞a_1<∞. B3 suptsupp∈Br(p∗)‖∇ϕt(p)‖∞≤a2 _t _p∈ B_r(p^*)\|∇ _t(p)\|_∞≤ a_2 for some a2<∞a_2<∞. B4 ‖ψ1(p)−ψ1(q)‖∞≤c1supt‖pt−qt‖1\| _1(p)- _1(q)\|_∞≤ c_1 _t\|p_t-q_t\|_1 and ‖ϕ0(p)−ϕ0(q)‖∞≤c2supt‖pt−qt‖1\| _0(p)- _0(q)\|_∞≤ c_2 _t\|p_t-q_t\|_1 for some c1,c2<∞c_1,c_2<∞ and ∀p,q∈Br(p∗),p≠q∀ p,q∈ B_r(p^*),p≠ q. B5 suptsupp∈Br(p∗)‖∇pt‖1<a3 _t _p∈ B_r(p^*)\|∇ p_t\|_1<a_3, for some a3<∞a_3<∞. Let W=βW¯W=β W for some scaling β>0β>0, and let the constant λ:=2e2r+1[2βe‖W¯‖∞+2σ2(a1+a2)β‖∇W¯‖∞+c1+c21−eσ2βZ], λ:=2e^2r+1 [ 2βe\| W\|_∞\!+\! 2σ^2(a_1\!+\!a_2)β\|∇ W\|_∞\!+\!c_1\!+\!c_21-eσ^2β Z ], (23) for Z:=‖ΔW¯‖∞+a3‖∇W¯‖∞Z:=\| W\|_∞+a_3\|∇ W\|_∞ and eσ2βZ<1eσ^2β Z<1. If λ<1λ<1, then the map C is locally contractive near p∗p^* with respect to the metric dHd_H^S, i.e., dH((p),(q))≤λdH(p,q)∀p,q∈Br(p∗). d_H^S (C(p),C(q) )≤λ d_H^S(p,q) ∀ p,q∈ B_r(p^*). Additionally, for every p(0)∈Br(p∗)p^(0)∈ B_r(p^*), the iteration p(k+1)=(p(k))p^(k+1)=C(p^(k)) is guaranteed to converge to p∗p^* as k→∞k→∞. Before we proceed with proving Theorem 1, we remark that Assumptions B1-B5 are not forced, and are in fact common for parabolic PDEs such as those considered herein. Assumption B1 holds under bounded drift and reaction terms in the PDEs we consider. Assumptions B2 and B3 hold when the coefficients and the boundary conditions of the HJB equations are sufficiently regular, while Assumption B4 involves the stability of these solutions under perturbations of the HJB coefficients. Lastly, ‖∇pt‖1\|∇ p_t\|_1 is bounded when ptp_t, the solution to McKean-Vlasov dynamics, is of the Sobolev class W1,1W^1,1, which is also satisfied for sufficiently regular coefficients and boundary conditions. Proof of Theorem 1. Letting ℰt(p):=e−2(W∗pt)E_t(p):=e^-2(W*p_t) ∀t∈[0,1]∀ t∈[0,1], we write ((p))t=ℰt(p)φt(p)φ^t(p).(C(p))_t=E_t(p)\, _t(p)\, _t(p). Here, φt(p),φ^t(p) _t(p), _t(p) are viewed as the images of the maps φt,φ^t:→0 _t, _t:S _0, respectively. From the definition (18), one can verify that dH(p1p2,q1q2)≤dH(p1,q1)+dH(p2,q2)d_H(p_1p_2,q_1q_2)≤ d_H(p_1,q_1)+d_H(p_2,q_2) for pi,qi∈0,i∈1,2p_i,q_i _0,i∈\1,2\. Accordingly, dH((p),(q))≤suptdH(ℰt(p),ℰt(q)) d_H^S (C(p),C(q) )≤ _td_H (E_t(p),E_t(q) ) +suptdH(φt(p),φt(q))+suptdH(φ^t(p),φ^t(q)). + _td_H ( _t(p), _t(q) )\!+\! _td_H ( _t(p), _t(q) ). (24) We now investigate the terms on the right-hand side of (I-C) individually, obtaining their respective upper bounds in terms of dH(p,q)d_H^S(p,q) for p,qp,q in the neighborhood Br(p∗)B_r(p^*). A pertinent inequality for our purposes is dH(f,g) d_H(f,g) ≤‖log(f/g)‖∞+‖log(g/f)‖∞ ≤\| (f/g)\|_∞+\| (g/f)\|_∞ =2‖logf−logg‖∞∀f,g∈0. =2\| f- g\|_∞ ∀ f,g _0. (25) Another inequality of interest is ‖f−g‖1≤edH(f,g)−1, \|f-g\|_1≤ e^d_H(f,g)-1, (26) for all f,g∈0f,g _0 satisfying ‖f‖1=‖g‖1\|f\|_1=\|g\|_1 (cf. Corollary 1 [26] and Lemma 2.2 [29]). Using (25), we have dH(ℰt(p),ℰt(q)) d_H (E_t(p),E_t(q) ) ≤4‖β(W¯∗(pt−qt))‖∞ ≤ 4\|β ( W*(p_t-q_t) )\|_∞ ≤4β‖W¯‖∞‖pt−qt‖1, ≤ 4β\| W\|_∞\|p_t-q_t\|_1, (27) where we have used Young’s convolution inequality [30, Proposition 2.33] and B1. Then, by (26), we get dH(ℰt(p),ℰt(q))≤4β‖W¯‖∞(edH(pt,qt)−1). d_H (E_t(p),E_t(q) )≤ 4β\| W\|_∞ (e^d_H(p_t,q_t)-1 ). If p,q∈Br(p∗)p,q∈ B_r(p^*), then dH(pt,qt)≤2rd_H(p_t,q_t)≤ 2r, and it follows that edH(pt,qt)−1≤e2rdH(pt,qt)e^d_H(p_t,q_t)-1≤ e^2rd_H(p_t,q_t). 555ex−1≤xexe^x-1≤ xe^x for x≥0x≥ 0. Thus, for all p,q∈Br(p∗)p,q∈ B_r(p^*), suptdH(ℰt(p),ℰt(q)) \!\! _td_H (E_t(p),E_t(q) ) ≤4β‖W¯‖∞e2rdH(p,q). ≤ 4β\| W\|_∞e^2rd_H^S(p,q). (28) Let γt:=logφt(p)−logφt(q)≡ψt(p)−ψt(q) _t:= _t(p)- _t(q)≡ _t(p)- _t(q). Using (25), we have dH(φt(p),φt(q))≤2‖γt‖∞d_H ( _t(p), _t(q) )≤ 2\| _t\|_∞, and so suptdH(φt(p),φt(q))≤2supt‖γt‖∞. _td_H ( _t(p), _t(q) )≤ 2 _t\| _t\|_∞. (29) Next, we derive an upper bound for supt‖γt‖∞ _t\| _t\|_∞, which by (29), will give an upper bound for suptdH(φt(p),φt(q)) _td_H ( _t(p), _t(q) ). From (I) and by direct calculations, we get ∂tγt _t _t +σ22⟨∇ψt(q)+∇ψt(p)−2W∗pt,∇γt⟩ + σ^22 ∇ _t(q)+∇ _t(p)-2W*p_t,∇ _t +σ22Δγt=Lt+Mt+Rt, + σ^22 _t=L_t+M_t+R_t, (30) where Lt:=σ2⟨∇ψt(q),∇W∗(pt−qt)⟩,L_t:=σ^2 ∇ _t(q),∇ W*(p_t-q_t) , Mt:=σ2∫Ω⟨∇ψt(p)(y),∇W(x−y)⟩(qt−pt)(y)dy,M_t:=σ^2 _ ∇ _t(p)(y),∇ W(x-y) (q_t-p_t)(y)\, dy, and Rt=σ2∫Ωγt(y)∇y⋅(∇W(x−y)qt(y))dy.R_t=σ^2 _ _t(y) _y· (∇ W(x-y)\,q_t(y) )\, dy. We let ηs:=γ1−s _s:= _1-s for s∈[0,1]s∈[0,1], i.e., η is the time reversal of γ. Then, ηs _s solves the forward parabolic PDE −∂sηs - _s _s +σ22⟨∇ψ1−s(q)+∇ψ1−s(p)−2W∗p1−s,∇ηs⟩ + σ^22 ∇ _1-s(q)+∇ _1-s(p)-2W*p_1-s,∇ _s +σ22Δηs=L1−s+M1−s+R1−s. + σ^22 _s=L_1-s\!+\!M_1-s+R_1-s. (31) The previous PDE, whenever its coefficients are bounded, satisfies the conditions of [31, Theorem 2.10], which establishes sups∈[0,1]‖ηs‖∞≤e(supf|η|+sups∈[0,1]‖L1−s+M1−s+R1−s‖∞), _s∈[0,1]\| _s\|_∞≤ e ( _ P_f|η|+\!\! _s∈[0,1]\!\!\|L_1-s\!+\!M_1-s\!+\!R_1-s\|_∞ ), where e≈2.71828e≈ 2.71828 is the Euler number, and f:=(Ω×0)⋃(∂Ω×(0,1)) P_f:=( ×\0\) (∂ ×(0,1)), representing a union set of spatial and temporal boundaries. Equivalently, supt∈[0,1]‖γt‖∞≤e(supb|γ|+supt∈[0,1]‖Lt+Mt+Rt‖∞), _t∈[0,1]\| _t\|_∞≤ e (\! _ P_b|γ|\!+\! _t∈[0,1]\!\!\|L_t\!+\!M_t\!+\!R_t\|_∞ ), (32) where b:=(Ω×1)⋃(∂Ω×(0,1)) P_b:=( ×\1\) (∂ ×(0,1)). Since the same Dirichelt spatial boundary conditions can be enforced on the ψ(p),ψ(q)ψ(p),ψ(q), then supx∈∂Ω|ψt(p)−ψt(q)|=0 _x∈∂ | _t(p)- _t(q)|=0 for all t∈[0,1]t∈[0,1]. Accordingly, only bounds on the difference of ψ at the temporal boundary, namely, ‖ψ1(p)−ψ1(q)‖∞\| _1(p)- _1(q)\|_∞, remain. By Assumption B4, sup∂(Ω×[0,1])|γ|≤c1e2rdH(p,q)∀p,q∈Br(p∗). _∂( ×[0,1])|γ|≤ c_1e^2rd_H^S(p,q) ∀ p,q∈ B_r(p^*). (33) Since ‖Lt‖∞≤σ2‖∇ψt(q)‖∞‖∇W∗(pt−qt)‖∞\|L_t\|_∞≤σ^2\|∇ _t(q)\|_∞\|∇ W*(p_t-q_t)\|_∞, then by B2 and Young’s convolution inequality, it holds that ‖Lt‖∞≤σ2a1β‖∇W¯‖∞‖pt−qt‖1.\|L_t\|_∞≤σ^2a_1β\|∇ W\|_∞\|p_t-q_t\|_1. By virtue of (26), supt‖Lt‖∞≤σ2a1β‖∇W¯‖∞e2rdH(p,q), _t\|L_t\|_∞≤σ^2a_1β\|∇ W\|_∞e^2rd_H^S(p,q), (34) for all p,q∈Br(p∗)p,q∈ B_r(p^*). Likewise, supt‖Mt‖∞≤σ2a1β‖∇W¯‖∞e2rdH(p,q). _t\|M_t\|_∞≤σ^2a_1β\|∇ W\|_∞e^2rd_H^S(p,q). (35) Further, by B5, supt‖Rt‖∞≤σ2βZsupt∈[0,1]‖γt‖∞, _t\|R_t\|_∞≤σ^2β Z _t∈[0,1]\| _t\|_∞, (36) for Z=‖ΔW¯‖∞+a3‖∇W¯‖∞Z=\| W\|_∞+a_3\|∇ W\|_∞. By the triangle inequality, ‖Lt+Mt+Rt‖∞≤‖Lt‖∞+‖Mt‖∞+‖Rt‖∞\|L_t\!+\!M_t+R_t\|_∞≤\|L_t\|_∞+\|M_t\|_∞+\|R_t\|_∞. Hence, by summing (33)-(36), the inequality (32) gives the sought-after bound on supt∈[0,1]‖γt‖∞ _t∈[0,1]\| _t\|_∞. Then, by (29), we establish that suptdH(φt(p),φt(q))≤ _td_H ( _t(p), _t(q) )≤ 2e2r+11−eσ2βZ(2σ2a1β‖∇W¯‖∞+c1)dH(p,q), 2e^2r+11-eσ^2β Z (2σ^2a_1β\|∇ W\|_∞+c_1 )d_H^S(p,q), (37) for all p,q∈Br(p∗)p,q∈ B_r(p^*). All steps used to obtain the previous relation apply verbatim for φ except for the need to revert the time to invoke [31, Theorem 2.10]. The result is suptdH(φ^t(p),φ^t(q))≤ _td_H ( _t(p), _t(q) )≤ 2e2r+11−eσ2βZ(2σ2a2β‖∇W¯‖∞+c2)dH(p,q), 2e^2r+11-eσ^2β Z (2σ^2a_2β\|∇ W\|_∞+c_2 )d_H^S(p,q), (38) for all p,q∈Br(p∗)p,q∈ B_r(p^*) and eσ2βZ<1eσ^2β Z<1. Plugging (28), (I-C) and (I-C) into (I-C) yields the constant λ>0λ>0 in (23). As long as p(0)∈Br(p∗)p^(0)∈ B_r(p^*) and λ<1λ<1, dH(p(1),p∗)=dH((p(0)),(p∗))≤λdH(p(0),p∗),d_H^S (p^(1),p^* )=d_H^S (C(p^(0)),C(p^*) )≤λ\,d_H^S(p^(0),p^*), and p(1)p^(1) stays in Br(p∗)B_r(p^*). Hence, dH(p(k),p∗)≤λkdH(p(0),p∗)=λkr.d_H^S (p^(k),p^* )≤λ^kd_H^S(p^(0),p^*)=λ^kr. Since λk→0λ^k→ 0 as k→∞k→∞, the fixed point iteration p(k+1)=(p(k))p^(k+1)=C(p^(k)) converges to p∗p^*. This concludes the proof. ∎ Theorem 2. Let ∗(p)=(φ∗(p),φ^∗(p)) ^*(p)=( ^*(p), ^*(p)) be a fixed point of the map pG_p defined in (22) for a fixed p. For ρ>0ρ>0, define the neighborhood Vρp:=(p)∈×:dH⊕((p),∗(p))≤ρ.V_ρ^p:=\ (p) ×A:d_H ( (p), ^*(p) )≤ρ\. Assume that B1 holds, and W=βW¯W=β W for some β>0β>0. Assume also that there exists a neighborhood Uρp⊆×U_ρ^p ×A such that Uρp⊇Vρp∪p(Vρp)U_ρ^p V_ρ^p _p(V_ρ^p). Suppose there exist mi=mi(ρ,p)>0,i=1,⋯,4m_i=m_i(ρ,p)>0,i=1,·s,4 such that for all (p)=(φ(p),φ^(p)),(p)=(μ(p),μ^(p))∈Uρp (p)=( (p), (p)), μ(p)=(μ(p), μ(p))∈ U_ρ^p and (p)≠(p) (p)≠ μ(p), ‖logφ1−logμ1‖∞≤m2supt∈[0,1−δ]‖logφt−logμt‖∞, \!\!\!\| _1\!-\! _1\|_∞≤ m_2 _t∈[0,1-δ]\| _t\!-\! _t\|_∞, (39) supt‖∇logφt−∇logμt‖∞≤m1suptdH(φt,μt), \!\!\! _t\|∇ _t\!-\!∇ _t\|_∞≤ m_1 _td_H( _t, _t), (40) for some δ>0δ>0, and analogously for φ^,μ , μ and their initial conditions with m3m_3 and m4m_4. Let the constant Λp=max2em1σ2β‖∇W¯‖∞1−em2,2em3σ2β‖∇W¯‖∞1−em4, \!\!\!\! _p= \ 2em_1σ^2β\|∇ W\|_∞1-em_2, 2em_3σ^2β\|∇ W\|_∞1-em_4 \, (41) with em2<1em_2<1 and em4<1em_4<1. If Λp<1 _p<1, then the map pG_p is locally contractive near ∗(p) ^*(p) with respect to the metric dH⊕d_H , i.e., dH⊕(p((p)),p((p)))≤ΛpdH⊕((p),(p)), d_H (G_p( (p)),G_p( μ(p)) )≤ _pd_H ( (p), μ(p)), for all (p),(p)∈Vρp (p), μ(p)∈ V_ρ^p. Additionally, for any 0(p)∈Vρp ^0(p)∈ V_ρ^p, the sequence j+1(p)=p(j(p)) ^j+1(p)=G_p( ^j(p)) converges to ∗ ^* as j→∞j→∞, up to multiplication of φ∗(p) ^*(p) by a positive constant and division of φ^∗(p) ^*(p) by the same constant. Proof. Let j(p),j(p)∈Vρp ^j(p), μ^j(p)∈ V_ρ^p. We note that j+1=p(j) ^j+1=G_p( ^j) is a limit of the innermost (Sinkhorn) iteration where φtj,φ^tj _t^j, _t^j appear in the reaction terms, and therefore, φtj+1 _t^j+1 satisfies (15a) subject to a final boundary condition φ1j+1 _1^j+1. Analogously, if j+1=p(j) μ^j+1=G_p( μ^j), then j+1 μ^j+1 is a solution to (15a) where μtj,μ^tj _t^j, μ_t^j appear in the reaction terms with the terminal condition μ1j+1 _1^j+1. If we set ζtj+1:=logφtj+1−logμtj+1 _t^j+1:= _t^j+1- _t^j+1, then we have ∂tζtj+1−σ22⟨2(∇W∗pt) _t _t^j+1\!-\! σ^22 2(∇ W*p_t) −∇logφtj+1−∇logμtj+1, -\!∇ _t^j+1-\!∇ _t^j+1, ∇ζtj+1⟩+σ22Δζtj+1=Ftj, ∇ _t^j+1 + σ^22 _t^j+1=F_t^j, where Ftj(x):=−σ2∫Ω⟨∇ζtj(y),∇W(x−y)⟩pt(y)dy,F_t^j(x):=-σ^2 _ ∇ _t^j(y),∇ W(x-y) p_t(y)\, dy, with ζtj:=logφtj−logμtj _t^j:= _t^j- _t^j. Provided that coefficients in the previous PDE are bounded, we can apply again on Theorem 2.10 [31] to obtain supt‖ζtj+1‖∞≤e‖ζ1j+1‖∞+esupt‖Ftj‖∞. _t\| _t^j+1\|_∞≤ e\| _1^j+1\|_∞+e _t\|F_t^j\|_∞. (42) By (39), we have ‖ζ1j+1‖∞≤m2supt∈[0,1−δ]‖ζtj+1‖∞≤m2supt∈[0,1]‖ζtj+1‖∞. \!\!\!\!\| _1^j+1\|_∞≤ m_2\!\! _t∈[0,1-δ]\!\| _t^j+1\|_∞≤ m_2\!\! _t∈[0,1]\!\| _t^j+1\|_∞. (43) Using (40) and Young’s convolution inequality, we get ‖Ftj‖∞≤σ2β‖∇ζtj‖∞‖∇W¯‖∞. \|F_t^j\|_∞≤σ^2β\|∇ _t^j\|_∞\|∇ W\|_∞. (44) Thus, through (42)-(44), we know that supt‖ζtj+1‖∞≤em1σ2β‖∇W¯‖∞1−em2suptdH(φtj,μtj), _t\| _t^j+1\|_∞≤ em_1σ^2β\|∇ W\|_∞1-em_2 _td_H( _t^j, _t^j), (45) for em2<1em_2<1. Using (25), we arrive at suptdH(φtj+1,μtj+1)≤2em1σ2β‖∇W¯‖∞1−em2suptdH(φtj,μtj). _td_H( _t^j+1, _t^j+1)≤ 2em_1σ^2β\|∇ W\|_∞1-em_2 _td_H( _t^j, _t^j). (46) By applying similar arguments applied to φ^j+1(p) ^j+1(p) and μ^j(p) μ^j(p), one can get suptdH(φ^tj+1,μ^tj+1)≤2em3σ2β‖∇W¯‖∞1−em4suptdH(φ^tj,μ^tj). _td_H( _t^j+1, μ_t^j+1)≤ 2em_3σ^2β\|∇ W\|_∞1-em_4 _td_H( _t^j, μ_t^j). (47) By the definition of dH⊕d_H in (20), the estimates in (46) and (47) translate into a contraction estimate in dH⊕d_H , giving dH⊕(p(),p())≤ΛpdH⊕(,), d_H (G_p( ),G_p( μ))≤ _pd_H ( , μ), for all , , μ in VρpV_ρ^p and Λp _p is as given in (41). We set (p)=∗(p) μ(p)= ^*(p) and (p)=0(p) (p)= ^0(p). Hence, as long as 0(p)∈Vρp ^0(p)∈ V_ρ^p and Λp<1 _p<1, then it holds that dH⊕(1,∗)=dH⊕(p(0),p(∗))≤ΛpdH⊕(0,∗),d_H ( ^1, ^* )=d_H (G_p( ^0),G_p( ^*) )≤ _pd_H ( ^0, ^* ), and 1(p) ^1(p) belongs to Vρp⊆UρpV_ρ^p U_ρ^p. Hence, dH⊕(j,∗)≤ΛpjdH⊕(0,∗). d_H ( ^j, ^* )≤ _p^jd_H ( ^0, ^* ). Since Λpj→0 _p^j→ 0 as j→∞j→∞, then limj→∞dH⊕(j(p),∗(p))=0 _j→∞d_H ( ^j(p), ^*(p))=0. Specifically, φj(p)→κ1φ∗(p) ^j(p)→ _1 ^*(p) and φ^j(p)→κ2φ^∗(p) ^j(p)→ _2 ^*(p) for some κ1,κ2>0 _1, _2>0 due to the projective property of dHd_H and inherently dH⊕d_H . However, from (15b), κ1κ2exp(−2W∗pin) _1 _2 (-2W*p_ in) φ0∗(p)φ^0∗(p)=pin _0^*(p) _0^*(p)=p_ in =exp(−2W∗pin)φ0∗(p)φ^0∗(p). = (-2W*p_ in) _0^*(p) _0^*(p). Therefore, κ1=1/κ2 _1=1/ _2, which completes the proof. ∎ Theorem 3. Under the assumptions of Theorems 1 and 2, if p(0)∈Br(p∗)p^(0)∈ B_r(p^*) and λ<1λ<1 and if 0(p(k))∈Vρp(k) ^0(p^(k))∈ V_ρ^p^(k) for all k together with supp∈Br(p∗)Λp<1 _p∈ B_r^(p^*) _p<1. Then, Algorithm 1 converges locally to the fixed point (p∗,φ∗(p∗),φ^∗(p∗))(p^*, ^*(p^*), ^*(p^*)) up to constant scaling of φ∗(p∗),φ^∗(p∗) ^*(p^*), ^*(p^*). This fixed point satisfies the Schrödinger system (14), constituting a solution to the MFSB problem666The existence of a solution follows from the large deviations arguments and its equivalence to optimal control formulation in (2) [16, Proposition 1.1, Lemma 3.6] and by extension to the equivalent formalism here with Cole-Hopf transform.. Proof. Under the assumptions of Theorem 1, if p(0)∈Br(p∗)p^(0)∈ B_r(p^*) and λ<1λ<1, then p(k)∈Br(p∗)p^(k)∈ B_r(p^*) for all k, and the iteration p(k+1)=(p(k))p^(k+1)=C(p^(k)) converges to p∗p^* as k→∞k→∞. For a fixed p=p(k)p=p^(k), if (p(k))∈Vρp(k) ^0(p^(k))∈ V_ρ^p^(k), and supp∈Br(p∗)Λp<1 _p∈ B_r(p^*) _p<1. Then, by Theorem 2, we know that the iteration +(p(k))=p(k)((p(k))) ^j+1(p^(k))=G_p^(k) ( ^j(p^(k)) ) converges to ∗(p(k)) ^*(p^(k)), up to the constant scaling indicated in Theorem 2, as j→∞j→∞. Further, by (I-C), taking φ(p)=φ∗(p(k)),φ(q)=φ∗(p∗) (p)= ^*(p^(k)), (q)= ^*(p^*) limk→∞suptdH(φt∗(p(k)),φt∗(p∗))=0. _k→∞ _td_H ( _t^*(p^(k)), _t^*(p^*) )=0. Analogously, by (I-C), we get limk→∞suptdH(φ^t∗(p(k)),φ^t∗(p∗))=0. _k→∞ _td_H ( _t^*(p^(k)), _t^*(p^*) )=0. Consequently, limk→∞(p(k),φ∗(p(k)),φ∗(p(k)))=(p∗,φ∗(p∗),φ^∗(p∗)), _k→∞(p^(k), ^*(p^(k)), ^*(p^(k)) )=(p^*, ^*(p^*), ^*(p^*) ), up to the constant positive scaling of φ∗(p∗),φ^∗(p∗) ^*(p^*), ^*(p^*) This establishes local convergence to the fixed point satisfying (14). ∎ Figure 1: Main plot: evolution of the controlled state PDFs ptp_t. Left inset: repulsive interaction potential: 5/(2|x|0.2)5/(2|x|^0.2). Right inset: convergence with respect to dHd_H^S. Figure 2: Main plot: evolution of the controlled state PDFs ptp_t. Left inset: attractive interaction potential: −exp(−|x|2/0.3)- (-|x|^2/0.3). Right inset: convergence with respect to dHd_H^S. From a vantage point, Theorems 1-3 are best viewed in a perturbative sense. In that, we see that the λ,Λp<1λ, _p<1 can be achieved for sufficiently weak interaction (β is small) and/or close initial guess to (p∗,φ∗(p∗),φ^∗(p∗))(p^*, ^*(p^*), ^*(p^*)). Heuristically, the latter can be obtained by successively applying Algorithm 1 to obtain a solution for a small β, which can then be used to initialize the algorithm for relatively larger β values. We also note that whenever C is contractive, the damping in (17) delays updating (p(k))C(p^(k)), effectively preventing the density p(k+1)p^(k+1) from escaping the contraction neighborhood due to possible numerical overshooting or oscillations. IV Numerical Results We report two 1D examples for the MFSB problem (2). In both, we use pin(x)∝0.5(exp(−(x−0.5)2/0.08)+exp(−(x+0.4)2/0.08))p_ in(x) 0.5 ( (-(x-0.5)^2/0.08 .)+ (-(x+0.4)^2/0.08 .) ), pfin(x)∝exp(−(x−0.4)2/0.08)p_ fin(x) (-(x-0.4)^2/0.08 .), tol=10−6tol=10^-6. We use the respective non-interacting SB solutions to initialize Algorithm 1. In the first example, we consider the repulsive interaction potential 5/(2|x|0.2)5/(2|x|^0.2), σ2=0.2σ^2=0.2, and θ=0.7θ=0.7. This example is comparable to that in [17, Sec. VI-B]. Our results are depicted in Fig. 1. For the second example, we consider the attractive potential W=−exp(−|x|2/0.3)W=- (-|x|^2/0.3), σ2=1σ^2=1, and θ=1θ=1. The corresponding numerical results are shown in Fig. 2. These figures delineate that the solutions obtained from Algorithm 1 match the given marginals with rapid convergence and high accuracy for both potentials. V Conclusions In this work, we propose a generalized Sinkhorn algorithm to solve the mean-field Schrödinger bridge problem. The solution enables distribution steering via feedback control for interacting multi-agent systems in the large population limit. For the proposed algorithm, we provide a guarantee of local convergence and illustrate it with numerical examples. References [1] H. P. McKean Jr, “A class of Markov processes associated with nonlinear parabolic equations,” Proceedings of the National Academy of Sciences, vol. 56, no. 6, p. 1907–1911, 1966. [2] Y. Chen, T. T. Georgiou, and M. Pavon, “Optimal steering of a linear stochastic system to a final probability distribution, Part I,” IEEE Transactions on Automatic Control, vol. 61, no. 5, p. 1170–1180, 2015. [3] Y. Chen, T. T. Georgiou, and M. Pavon, “Stochastic control liaisons: Richard Sinkhorn meets Gaspard Monge on a Schrodinger Bridge,” SIAM Review, vol. 63, no. 2, p. 249–313, 2021. [4] K. F. Caluya and A. Halder, “Finite horizon density steering for multi-input state feedback linearizable systems,” in 2020 American Control Conference (ACC). IEEE, 2020, p. 3577–3582. [5] K. F. Caluya and A. Halder, “Wasserstein proximal algorithms for the Schrödinger bridge problem: Density control with nonlinear drift,” IEEE Transactions on Automatic Control, vol. 67, no. 3, p. 1163–1178, 2021. [6] I. Nodozi, C. Yan, M. Khare, A. Halder, and A. Mesbah, “Neural Schrödinger bridge with Sinkhorn losses: Application to data-driven minimum effort control of colloidal self-assembly,” IEEE Transactions on Control Systems Technology, vol. 32, no. 3, p. 960–973, 2023. [7] A. Teter and A. Halder, “On the Hopf-Cole transform for control-affine Schrödinger bridge,” arXiv preprint arXiv:2503.17640, 2025. [8] Y. Chen, T. T. Georgiou, and M. Pavon, “Optimal steering of a linear stochastic system to a final probability distribution—Part I,” IEEE Transactions on Automatic Control, vol. 63, no. 9, p. 3112–3118, 2018. [9] A. M. H. Teter, W. Wang, and A. Halder, “Schrödinger bridge with quadratic state cost is exactly solvable,” IEEE Transactions on Automatic Control, p. 1–15, 2025. [10] A. M. Teter, W. Wang, S. Shivakumar, and A. Halder, “Markov kernels, distances and optimal control: A parable of linear quadratic non-Gaussian distribution steering,” arXiv preprint arXiv:2504.15753, 2025. [11] A. M. Teter, I. Nodozi, and A. Halder, “Probabilistic Lambert problem: Connections with optimal mass transport, Schrödinger bridge, and reaction-diffusion PDEs,” SIAM Journal on Applied Dynamical Systems, vol. 24, no. 1, p. 16–43, 2025. [12] K. F. Caluya and A. Halder, “Reflected Schrödinger bridge: Density control with path constraints,” in 2021 American Control Conference (ACC). IEEE, 2021, p. 1137–1142. [13] A. Eldesoukey and T. T. Georgiou, “Schrödinger’s control and estimation paradigm with spatio-temporal distributions on graphs,” IEEE Transactions on Automatic Control, vol. 70, no. 4, p. 2466–2478, 2024. [14] A. Eldesoukey, O. M. Miangolarra, and T. T. Georgiou, “An excursion onto Schrödinger’s bridges: Stochastic flows with spatio-temporal marginals,” IEEE Control Systems Letters, vol. 8, p. 1138–1143, 2024. [15] O. Movilla Miangolarra, A. Eldesoukey, and T. T. Georgiou, “Inferring potential landscapes: A Schrödinger bridge approach to maximum caliber,” Physical Review Research, vol. 6, no. 3, p. 033070, 2024. [16] J. Backhoff, G. Conforti, I. Gentil, and C. Léonard, “The mean field Schrödinger problem: ergodic behavior, entropy estimates and functional inequalities,” Probability Theory and Related Fields, vol. 178, no. 1, p. 475–530, 2020. [17] Y. Chen, “Density control of interacting agent systems,” IEEE Transactions on Automatic Control, vol. 69, no. 1, p. 246–260, 2023. [18] R. M. Colombo and M. Lécureux-Mercier, “Nonlocal crowd dynamics models for several populations,” Acta Mathematica Scientia, vol. 32, no. 1, p. 177–196, 2012. [19] C. M. Topaz, A. L. Bertozzi, and M. A. Lewis, “A nonlocal continuum model for biological aggregation,” Bulletin of mathematical biology, vol. 68, no. 7, p. 1601–1623, 2006. [20] Z. Zhang, Z. Wang, Y. Sun, J. Shen, Q. Peng, T. Li, and P. Zhou, “Deciphering cell-fate trajectories using spatiotemporal single-cell transcriptomic data,” npj Systems Biology and Applications, vol. 12, no. 1, 2025. [21] K. Elamvazhuthi and S. Berman, “Mean-field models in swarm robotics: A survey,” Bioinspiration & Biomimetics, vol. 15, no. 1, p. 015001, 2020. [22] E. Hopf, “The partial differential equation ut+uux=μxxu_t+u_x= _x,” Communications on Pure and Applied Mathematics, vol. 3, no. 3, p. 201–230, 1950. [23] J. D. Cole, “On a quasi-linear parabolic equation occurring in aerodynamics,” Quarterly of Applied Mathematics, vol. 9, no. 3, p. 225–236, 1951. [24] M. A. Krasnosel’skii, “Two remarks on the method of successive approximations,” Uspekhi matematicheskikh nauk, vol. 10, no. 1, p. 123–127, 1955. [25] W. R. Mann, “Mean value methods in iteration,” Proceedings of the American Mathematical Society, vol. 4, no. 3, p. 506–510, 1953. [26] G. Birkhoff, “Extensions of Jentzsch’s theorem,” Transactions of the American Mathematical Society, vol. 85, no. 1, p. 219–227, 1957. [27] T. T. Georgiou and M. Pavon, “Positive contraction mappings for classical and quantum Schrödinger systems,” Journal of Mathematical Physics, vol. 56, no. 3, 2015. [28] Y. Chen, T. Georgiou, and M. Pavon, “Entropic and displacement interpolation: a computational approach using the Hilbert metric,” SIAM Journal on Applied Mathematics, vol. 76, no. 6, p. 2375–2396, 2016. [29] S. Eckstein, “Hilbert’s projective metric for functions of bounded growth and exponential convergence of Sinkhorn’s algorithm,” Probability Theory and Related Fields, vol. 193, no. 1, p. 585–621, 2025. [30] J. van Neerven, Functional analysis. Cambridge University Press, 2022, vol. 201. [31] G. M. Lieberman, Second order parabolic differential equations. World scientific, 1996.