Paper deep dive
Nonlinear Non-Gaussian Density Steering with Input and Noise Channel Mismatch: Sinkhorn with Memory for Solving the Control-affine Schrödinger Bridge Problem
Georgiy A. Bondar, Asmaa Eldesoukey, Yongxin Chen, Abhishek Halder
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 98%
Last extracted: 6/21/2026, 6:42:34 AM
Summary
This paper addresses the control-affine Schrödinger bridge problem (CASBP) specifically focusing on the 'channel mismatch' scenario where the input coefficient and noise channel are not proportional. While standard Sinkhorn recursion relies on the Hopf-Cole transform to linearize the problem, this mismatch results in nonlinear PDEs. The authors propose a novel 'Sinkhorn recursion with memory' that leverages the structure of these nonlinear PDEs to solve the CASBP. They provide a proof of local stability for the proposed algorithm and demonstrate its effectiveness through numerical examples, establishing that the CASBP is locally nonlinearly solvable via this method.
Entities (10)
Relation Signals (4)
channel mismatch → causesnonlinearityin → Hopf-Cole-transformed PDEs
confidence 100% · When the channels do not match, the Hopf-Cole-transformed PDEs remain nonlinear
Georgiy A. Bondar → isaffiliatedwith → University of California, Santa Cruz
confidence 100% · Georgiy A. Bondar is with the Department of Applied Mathematics, University of California Santa Cruz
Sinkhorn recursion with memory → solves → control-affine Schrödinger bridge problem
confidence 100% · designing a Sinkhorn recursion with memory that leverages the structure of these nonlinear PDEs, and demonstrate how it solves the control-affine Schrödinger bridge problem
Hopf-Cole transform → isusedin → Sinkhorn recursion
confidence 90% · The mathematical engine behind this approach is the Hopf-Cole transform that recasts the conditions for optimality into a system of boundary-coupled linear PDEs.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Solutions to the Schrödinger bridge problem and its generalizations yield feedback control policies for optimal density steering over a controlled diffusion. To numerically compute the same, the dynamic Sinkhorn recursion has become a standard approach. The mathematical engine behind this approach is the Hopf-Cole transform that recasts the conditions for optimality into a system of boundary-coupled linear PDEs. Recent works pointed out that for the control-affine Schrödinger bridge problem, this exact linearity via Hopf-Cole transform, and thus the standard Sinkhorn recursion, apply only if the control and noise channels are proportional. When the channels do not match, the Hopf-Cole-transformed PDEs remain nonlinear, and no algorithm is available to solve the same. We advance the state-of-the-art by designing a Sinkhorn recursion with memory that leverages the structure of these nonlinear PDEs, and demonstrate how it solves the control-affine Schrödinger bridge problem with input and noise channel mismatch. We prove the local stability of the proposed algorithm.
Tags
Links
- Source: https://arxiv.org/abs/2604.23370v1
- Canonical: https://arxiv.org/abs/2604.23370v1
Trouble viewing inline? Open PDF directly →
Full Text
83,052 characters extracted from source content.
Expand or collapse full text
bstctl:etal, bstctl:nodash, bstctl:simpurl _b:BSTcontrol Nonlinear Non-Gaussian Density Steering with Input and Noise Channel Mismatch: Sinkhorn with Memory for Solving the Control-affine Schrödinger Bridge Problem Georgiy A. Bondar Asmaa Eldesoukey Yongxin Chen Abhishek Halder Georgiy A. Bondar is with the Department of Applied Mathematics, University of California Santa Cruz, CA 95064, USA, gbondar@ucsc.edu.Abhishek Halder and Asmaa Eldesoukey are with the Department of Aerospace Engineering, Iowa State University, Ames, IA 50011, USA, ahalder,asmaae@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 Solutions to the Schrödinger bridge problem and its generalizations yield feedback control policies for optimal density steering over a controlled diffusion. To numerically compute the same, the dynamic Sinkhorn recursion has become a standard approach. The mathematical engine behind this approach is the Hopf-Cole transform that recasts the conditions for optimality into a system of boundary-coupled linear PDEs. Recent works pointed out that for the control-affine Schrödinger bridge problem, this exact linearity via Hopf-Cole transform, and thus the standard Sinkhorn recursion, apply only if the control and noise channels are proportional. When the channels do not match, the Hopf-Cole-transformed PDEs remain nonlinear, and no algorithm is available to solve the same. We advance the state-of-the-art by designing a Sinkhorn recursion with memory that leverages the structure of these nonlinear PDEs, and demonstrate how it solves the control-affine Schrödinger bridge problem with input and noise channel mismatch. We prove the local stability of the proposed algorithm. IEEEkeywords Schrödinger bridge, Sinkhorn algorithm, diffusion, density steering. 1 Introduction The purpose of this work is to investigate the control-affine Schrödinger bridge problem (CASBP) [1, 2] from a computational perspective. The CASBP is a stochastic optimal control problem of the form arginf(ρ,)∈01×∫t0t1ρ[q(t,u)+12‖22]dt (ρ u, u) _01×U _t_0^t_1E_ρ u [q(t, x^u)+ 12\| u\|_2^2 ]\> dt (1a) subject to d=((t,)+(t,))dt+(t,)d, x u= ( f (t, x u )+ g (t, x u ) u )dt+ σ (t, x u )d w, (1b) where 0≤t0<t1<∞0≤ t_0<t_1<∞, the control input ∈ℝm u ^m, the controlled state ∈ℝn x u ^n, the standard Wiener process ∈ℝp w ^p, and for given probability density functions (PDFs) ρ0,ρ1 _0, _1 with finite second moments, the feasible sets 01 \!\!P_01 :=t↦ρ(t,⋅)∈1([t0,t1])∣ρ≥0,∫ρ(t,⋅)=1 :=\t ρ(t,·) ^1([t_0,t_1]) ρ≥ 0, \!ρ(t,·)=1 ∀t∈[t0,t1],ρ(t0,⋅)=ρ0(⋅),ρ(t1,⋅)=ρ1(⋅), ∀ t∈[t_0,t_1],ρ (t_0,· )= _0(·),ρ (t_1,· )= _1(·)\, (2a) :=:[t0,t1]×ℝn↦ℝm∣∥22<∞. := \ u: [t_0,t_1 ]×R^n ^m \| u\|_2^2<∞ \. (2b) The state cost q≥0q≥ 0 is assumed to be bounded and integrable with respect to the measure ρdρ u d x u for all (ρ,)∈01×(ρ u, u) _01×U. So the design objective is to guide the state from the initial PDF ρ0 _0 to the terminal PDF ρ1 _1 over a given finite time horizon [t0,t1][t_0,t_1] via finite energy controls while minimizing an average additive cost (sum of a state cost q≥0q≥ 0 and a control cost 12‖22 12\| u\|_2^2). As standard, the positive semidefinite matrix field :=⊤⪰ := σ 0 is referred to as the diffusion tensor. We make the following standard regularity assumptions on the drift, input, and noise coefficients ,, f, g, σ, respectively. A1 (Non-explosion and Lipschitz coefficients) There exist finite constants c1,c2>0c_1,c_2>0 such that ∀,∈ℝn∀ x, y ^n, ∀t∈[t0,t1]∀ t∈[t_0,t_1], we have ‖(t,)‖2+‖(t,)‖2≤c1(1+‖2)\| f(t, x)\|_2+\| σ (t, x )\|_2≤ c_1 (1+\| x\|_2 ), and ‖(t,)−(t,)‖2≤c2‖−‖2\| f(t, x)- f(t, y)\|_2≤ c_2\| x- y\|_2. A2 (Uniformly lower bounded diffusion tensor) There exists a finite constant c3>0c_3>0 such that ∀∈ℝn∀ x ^n, ∀t∈[t0,t1]∀ t∈[t_0,t_1], we have ⟨,(t,)⟩≥c3‖22 x, (t, x) x ≥ c_3\| x\|_2^2. A3 (Uniformly upper bounded input coefficient) There exists a finite constant c4>0c_4>0 such that ∀∈ℝn∀ x ^n, ∀t∈[t0,t1]∀ t∈[t_0,t_1], we have ‖(t,)‖2<c4\| g(t, x)\|_2<c_4. The CASBP generalizes the classical SBP [3, 4] in that the latter is the following special case of (1): =,==,q=0. f= 0,\, g= σ= I,\,q=0. Given problem data ,,,q,ρ0,ρ1 f, g, σ,q, _0, _1, the CASBP (1)-(2) yields unique Markovian dynamic state feedback policy ∈ u , see e.g., [5]. This control policy guarantees that the controlled state statistics go from ρ0 _0 to ρ1 _1 over the specified time horizon [t0,t1][t_0,t_1], in addition to minimizing the average cost (1a) subject to (1b). Channel mismatch Our particular focus is on solving a generic CASBP with “channel mismatch”, by which we mean that the input coefficient g and the diffusion coefficient σ are not mutually related. In particular, ⊤⪰ g 0 is not proportional to the diffusion tensor . Stated differently, ∄λ>0such thatλ⊤−=. λ>0 that λ g - = 0. (3) This is the situation when the noise and the input do not act on the same subspace [1, Sec. IV-B]. We clarify that the channel mismatch issue is pertinent for n≥2n≥ 2 states, i.e., when the size of the square matrices ⊤, g , are at least 2×22× 2. For n=1n=1, the input and noise channels are trivially matched. Thus, in this work, we assume n≥2n≥ 2. Table 1: Summary of the control-affine Schrödinger bridge problems (CASBPs) and the nature of their solutions CASBP type Solution type Solvable Algorithm Performance guarantee Same input and noise channels linearly dynamic Sinkhorn recursion global convergence Different input and noise channels nonlinearly dynamic Sinkhorn recursion with memory local stability (proposed in this work, Sec. 3.4) (Sec. 4) Related works Several works [5, 6, 7, 8, 9, 10, 11] have studied the stochastic control-theoretic generalizations of the classical Schrödinger bridge, including the control-affine [12, 13, 1, 2] and control-non-affine [14, 15] cases. While works such as [12, 13] proposed solving the CASBP via dynamic Sinkhorn recursions, they made the assumption ⊤∝ g a priori – motivated partly by computational convenience, and partly by that the input and noise channels indeed coincide in practical situations such as when the process noise enters through stochastic actuation or external forcing. More recent works [1, 2] explicitly pointed out that for CASBP with channel mismatch, when the Hopf-Cole transformation [16, 17] is applied to the first order conditions for optimality, a boundary coupled system of linear PDEs are obtained, and that system can be solved by standard Sinkhorn algorithm, only if λ⊤−=λ g - = 0 for some λ>0λ>0. When the channels do not match, the Hopf-Cole-transformed PDE system remains nonlinear, and it is unclear whether a generalized variant of the Sinkhorn algorithm can be designed to solve the same. In the linear Gaussian (i.e., covariance steering) special case of (1), the work in [18] considered different input and noise channels, and derived111In the notation of reference [18], the term λ⊤−λ g - appears there as ⊤−11⊤ B - B_1 B_1 , where ,1 B, B_1 are the state-independent input and noise coefficient matrices, respectively. a system of coupled nonlinear matrix ODE boundary value problem [18, eq. (11)]. That reference approximated the solution of that matrix boundary value problem via time discretization and semidefinite programming. In contrast, the focus of this work is to solve nonlinear non-Gaussian CASBPs with channel mismatch. Motivation The technical motivation behind this work is to remedy the lack of a computational algorithm to solve the generic CASBP with channel mismatch. The channel mismatch is practical and also motivated by cases where the noise may represent modeling uncertainty in addition to disturbances or imperfections in actuation [19, eq. (7.1)], [20, eq. (90)], or when the input is affected by a lag in actuation dynamics [18, Sec. V, Example 2]. In such cases, the proportionality relation ⊤∝ g is violated. More broadly, solving stochastic optimal control problems of the form (1)-(2) is motivated by shaping state distributions via feedback control [21], i.e., direct control of uncertainties. For example, the specification of the initial PDF ρ0 _0 can be interpreted as the initial uncertainty obtained from an estimator, and that of the terminal PDF ρ1 _1 as the desired statistical performance to achieve over the specified time horizon [t0,t1][t_0,t_1]. Connections with linearly solvable stochastic optimal control The proportionality relation ⊤∝ g , i.e., λ⊤−=λ g - = 0 for some λ>0λ>0, has appeared before in the stochastic optimal control literature in a different context. Specifically, it was pointed out in [22, 23, 24] that when solving stochastic optimal control problems of the form (1) with the exception that the endpoint PDF constraints are removed, and a terminal cost of the form ρ[ϕ(u(t1))]E_ρ u [φ( x^u(t_1)) ] for some suitable ϕ(⋅)φ(·) is added in the objective (1a), the associated Hamilton-Jacobi-Bellman PDE can be exactly transformed to a linear PDE provided λ⊤−=λ g - = 0. Such stochastic optimal control problems are called linearly solvable [22, 23, 24] since by standard dynamic programming arguments [25, Ch. 4], the optimal control can then be recovered from the solution of this transformed linear PDE. When the relation ⊤∝ g does not hold, then the problem is not linearly solvable in the sense that computing the optimal control requires solving the second-order nonlinear Hamilton-Jacobi-Bellman PDE. We will explain in Sec. 3 that in the CASBP too, the lack of the same proportionality relation brings nonlinearity in a different way. Contributions Our contributions are the following. • To solve the system of nonlinear PDEs arising in generic CASBP with channel mismatch, we propose a new dynamic Sinkhorn algorithm (Sec. 3.4) that uses–in each epoch–memory from the most recent backward pass to be able to perform a forward pass in time. • We provide a local stability guarantee (Sec. 4) for the proposed Sinkhorn algorithm with memory. • We demonstrate the effectiveness of the proposed generalization of the Sinkhorn algorithm via numerical examples (Sec. 5). At a conceptual level, this work clarifies that when the input and noise channels match, then the CASBP is linearly solvable in the sense the problem reduces to solving a boundary coupled system of linear PDEs in the so-called Schrödinger factors, and standard dynamic Sinkhorn recursions apply with global convergence guarantees. In contrast, when these channels do not match, then the problem reduces to solving a boundary coupled system of nonlinear PDEs in the Schrödinger factors, which is shown to be solvable by dynamic Sinkhorn recursions with memory. However, in this more general case, we only establish local stability guarantees. In this sense, the CASBP then is locally nonlinearly solvable. For the readers’ convenience, this is summarized in Table 1. Organization Sec. 2 gathers some technical preliminaries that are used later in this work. In Sec. 3, we explain the main ideas in the context of state-of-the-art results from existing literature, and then propose the new algorithm. Sec. 4 undertakes an analysis of the proposed Sinkhorn algorithm with memory, and establishes a stability guarantee. Numerical results in Sec. 5 illustrate that the proposed algorithm works well in practice. Sec. 6 concludes the work. 2 Preliminaries We begin our exposition by collecting some concepts and results needed in the sequel. 2.1 Cones, Hilbert’s Projective Metric, and Contraction For D⊆ℝnD ^n, we consider the Banach space L∞(D)L^∞(D) with norm ∥⋅∥∞:=esssup∈D|⋅|\|·\|_∞:=ess _ x∈ D|·|. Throughout, we view all extrema and inequalities involving functions in L∞(D)L^∞(D) in the almost-everywhere sense, and to reduce clutter, we suppress the essess notation. Let K:=f∈L∞(D):f≥0 K:=\f∈ L^∞(D):f≥ 0\ (4) be the closed solid cone of nonnegative functions in L∞(D)L^∞(D), and denote the interior of K by K+K_+. That is, K+:=f∈L∞(D):inf∈Df()>0. K_+:=\f∈ L^∞(D): _ x∈ Df( x)>0\. (5) Given any u,v∈K+u,v∈ K_+, Hilbert’s projective metric dH(⋅,⋅)d_H(·,·) is defined as dH(u,v):=log(infθ>0∣u⪯θvsupθ>0∣θv⪯u) d_H(u,v):= ( \θ>0 u θ v\ \θ>0 θ v u\ ) (6) where ⪯ denotes the conic inequality, i.e., u⪯vu v if and only if v−u∈Kv-u∈ K, see [26, 27, 28]. Equivalently, dH(u,v)=logsupu()v()−loginfu()v()∀u,v∈K+. d_H(u,v)= _ x u( x)v( x)- _ x u( x)v( x) ∀ u,v∈ K_+. (7) The term projective refers to invariance under positive scalings, i.e., dH(u,v)=0d_H(u,v)=0 whenever u=cvu=cv for any c>0c>0. Definition 1 (p-Homogeneity). Consider K,K+K,K_+ as in (4)-(5). Let :K+→K+G:K_+→ K_+ be a map on the interior of K, i.e., a positive map. If (αu)=αp(u)∀u∈K+,α>0, (α u)=α^pG(u) ∀ u∈ K_+,\;α>0, (8) we say G is homogeneous of degree p (or p-homogeneous). Definition 2 (Global Contraction on K+K_+). Consider K,K+K,K_+ as in (4)-(5). The map :K+→K+G:K_+→ K_+ is called globally nonexpansive with respect to dH(⋅,⋅)d_H(·,·) if there exists 0<κ≤10<κ≤ 1 such that dH((u),(v))≤κdH(u,v)∀u,v∈K+. d_H(G(u),G(v))≤κ\,d_H(u,v) ∀ u,v∈ K_+. (9) The map G is said to be globally contractive if (9) holds for 0<κ⪇10<κ 1. The smallest κ for which the map is contractive is called the contraction ratio of G, denoted as κH() _ H(G). For a positive linear map G, a relevant quantity is the projective diameter ProjDiam():=supu,v∈K+dH((u),(v)). ProjDiam(G):= _u,v∈ K_+d_ H(G(u),G(v)). (10) An inequality [29, p. 2380] that will be used later is ProjDiam()≤2supu∈K+dH((u),v), ProjDiam(G)≤ 2 _u∈ K_+d_ H(G(u),v), (11) for any v∈K+v∈ K_+. Theorem 1 (Birkhoff-Bushell Theorem [27, 26]). Consider K,K+K,K_+ as in (4)-(5), the mapping :K+→K+G:K_+→ K_+, and the contraction ratio as in Definition 2. If G is homogeneous of degree p, and is monotone increasing, i.e., u⪯v⟹(u)⪯(v), u v (u) (v), then (9) holds with κH()≤p. _ H(G)≤ p. (12) In the special case G is linear, κH()≤tanh(14ProjDiam()). _ H(G)≤ ( 14ProjDiam(G) ). (13) That is, if :K+→K+G:K_+→ K_+ is linear, and ProjDiam()<∞ProjDiam(G)<∞, then G is globally contractive. 2.2 Comparison Principle for Parabolic PDEs Theorem 2 (Comparison Principle [30, Theorem 9.1]). For D⊂ℝnD ^n, consider a spatio-temporal domain Ω:=[t0,t1]×D :=[t_0,t_1]× D and its parabolic boundary Ω:=(t0×D¯)∪((t0,t1)×∂D), P :=(\t_0\× D)∪((t_0,t_1)×∂ D), (14) where D¯,∂D D,∂ D denote the topological closure and the boundary of D, respectively. Consider the quasilinear parabolic PDE −∂t - _t u(t,)+∑i,j=1naij(t,,u,∇u)∂xixju(t,) u(t, x)\!+\! _i,j=1^na_ij(t, x,u, _ xu) _x_ix_ju(t, x) +b(t,,u,∇u)=0,∀(t,)∈Ω. +b(t, x,u, _ xu)=0, ∀(t, x)∈ . (15) In (15), suppose that aija_ij are independent of u, and b is decreasing in u. Then, for smooth solutions u and v of (15), we have u≤v on Ω⟹u≤v on [t0,t1]×D¯. u≤ v on P u≤ v on [t_0,t_1]× D. (16) 3 Solution of the CASBP and its Computation In Sec. 3.1, we discuss the solution structure for the generic CASBP (1)-(2). In Sec. 3.2, we explain why the case ⊤∝ g is amenable to dynamic Sinkhorn recursion. Then in Sec. 3.3, we motivate our new idea of dynamic Sinkhorn recursion with memory for the case ⊤∝̸ g . The proposed algorithm is summarized in Sec. 3.4. 3.1 Solution Structure for the CASBP Let ⟨⋅,⋅⟩ ·,· denote the Frobenius inner product, and Δρ:=∑i,j=1n∂xixj(ij(t,)ρ(t,)) _ ρ:= _i,j=1^n _x_ix_j ( _ij(t, x)ρ(t, x) ) denote the weighted Laplacian. We use Δ to denote the standard Laplacian, i.e., Δ≡Δ ≡ _ I. It is known [1, Sec. I] that the first-order conditions for optimality for the CASBP (1)-(2) are ∂tρopt+∇⋅(ρopt(+⊤∇S))=12Δρopt, _tρ u_opt+ _ x· (ρ u_opt ( f+ g g _ xS ) )= 12 _ \>ρ u_opt, (primal PDE) (17a) ∂tS+⟨∇S,⟩+12⟨∇S,⊤∇S⟩+12⟨,HessS⟩=q, _tS+ _ xS, f \!+\! 12 _ xS, g _ xS + 12 ,Hess_ xS =q, (dual PDE) (17b) ρopt(t=t0,⋅)=ρ0(⋅),ρopt(t=t1,⋅)=ρ1(⋅). ρ u_opt (t=t_0,· )= _0(·), ρ u_opt (t=t_1,· )= _1(·). (primal boundary conditions) 97.56493pt(primal boundary conditions) (17c) The coupled system of equations (17) is in unknown primal-dual pair (ρopt,S) (ρ u_opt,S ), wherein ρoptρ u_opt is the optimally controlled joint state PDF, and S∈1,2([t0,t1],ℝn;ℝ)S∈ C^1,2 ([t_0,t_1],R^n;R ) is a Lagrange multiplier function222Here 1,2([t0,t1],ℝn;ℝ) C^1,2 ([t_0,t_1],R^n;R ) denotes the class of scalar-valued functions S(t,)S(t, x) which are once continuously differentiable with respect to time t∈[t0,t1]t∈[t_0,t_1] and twice continuously differentiable with respect to state ∈ℝn x ^n.. The optimal control is opt=⊤∇S. u_opt= g _ xS. (18) A direct numerical solution of (17) is challenging since the nonlinear PDEs (17a)-(17b) are coupled in primal-dual pair (ρopt,S) (ρ u_opt,S ), but the boundary conditions (17c) are in terms of ρoptρ u_opt only. To recast (17) in a more amenable form, we fix a parameter λ>0λ>0, and consider the Hopf-Cole transform [16, 17] (ρopt,S)↦(φ^,φ):=(ρoptexp(−S/λ),exp(S/λ)). (ρ u_opt,S ) ( , ):= (ρ u_opt (-S/λ ), (S/λ ) ). Applied to (17), this results in [1, Theorem 2] the transformed system: ∂tφ^+∇⋅(φ^(+φ))−12Δφ^+(qλ+qφ)φ^=0, \!\! _t + _ x· ( ( f+ f_ ) )- 12 _ \> + ( qλ+q_ ) =0, (19a) ∂tφ+⟨∇φ,+φ⟩+12⟨,Hessφ⟩−(qλ+qφ)φ=0, \!\! _t + _ x , f+ f_ + 12 ,Hess_ x - ( qλ+q_ ) =0, (19b) φ^(t0,⋅)φ(t0,⋅)=ρ0(⋅),φ^(t1,⋅)φ(t1,⋅)=ρ1(⋅), (t_0,· ) (t_0,· )= _0(·), (t_1,· ) (t_1,· )= _1(·), (19c) where φ f_ :=(λ⊤−)∇logφ, := (λ g - ) _ x , (20) qφ q_ :=12(∇logφ)⊤(λ⊤−)∇logφ. := 12 ( _ x ) (λ g - ) _ x . (21) The minimizing pair (ρopt,opt) (ρ u_opt, u_opt ) for the CASBP (1)-(2) is then recovered from ρopt(t,⋅) ρ u_opt(t,·) =φ^(t,⋅)φ(t,⋅), = (t,·) (t,·), (22) opt(t,⋅) u_opt(t,·) =λ⊤∇(⋅)logφ(t,⋅), =λ\> g _(·) (t,·), (23) for all t∈[t0,t1]t∈[t_0,t_1]. As standard, the pair (φ^,φ) ( , ) is referred to as the Schrödinger factors since by (22), they provide a factorization of ρoptρ u_opt, the optimally controlled joint state PDF. Remark 1 (Control weight matrix). It is easy to see that replacing 12‖22 12\| u\|_2^2 in (1a) by 12⊤ 12 u Ru with ≻ R 0, generalizes the channel mismatch condition (3) as ∄λ>0 λ>0 such that λ−1⊤−=λ g R^-1 g - = 0. Then the weight matrix in both (20) and (21) becomes (λ−1⊤−)(λ g R^-1 g - ), and the optimal control in (23) becomes opt(t,⋅)=λ−1⊤∇(⋅)logφ(t,⋅) u_opt(t,·)=λ R^-1 g _(·) (t,·). Remark 2 (Value of λ>0λ>0). When the channels match, then ⊤∝ g , and the numerical value of λ>0λ>0, by definition, is the proportionality constant in the relation λ⊤−=λ g - = 0. When the channels do not match, then in all developments starting from Sec. 3.3, we can set λ=1λ=1 without loss of generality. 3.2 The Case ⊤∝ g : Dynamic Sinkhorn Recursion Being a system of nonlinear PDE boundary value problem (BVP), the transformed system (19) at first glance appears to be as difficult as the original system (17). However, a closer inspection reveals an interesting structure in (19). When the input and noise channels match, i.e., λ⊤−=λ g - = 0 for some λ>0λ>0, then both (20) and (21) vanish, and the PDEs (19a)-(19b) become decoupled and linear. The only coupling for the BVP then comes from the bilinear boundary conditions (19c). This paves the way for the dynamic Sinkhorn recursion: φ^0→ℱφ^1→ρ1/φ^1φ1→ℬφ0→ρ0/φ0(φ^0)next, _0~ F~ _1~ _1/ _1~ _1~ B~ _0~ _0/ _0~ ( _0 )_next, (24) where φ^i(⋅):=φ^(t=ti,⋅),φi(⋅):=φ(t=ti,⋅)∀i∈0,1. _i(·):= (t=t_i,·), _i(·):= (t=t_i,·) ∀ i∈\0,1\. In (24), the forward-in-time map ℱF solves a linear PDE initial value problem (IVP) with (19a) from t0t_0 to t1t_1, with initial condition φ^0 _0. The backward-in-time map ℬB solves a linear PDE IVP with (19b) from t1t_1 to t0t_0, with initial condition φ1 _1. In both cases, solution is made possible by λ⊤−=λ g - = 0. It is well-known [29] that such dynamic Sinkhorn recursions are contractive with worst-case linear rate-of-convergence with respect to Hilbert’s projective metric dHd_H defined in (6). Specifically, the recursion (24) is guaranteed to converge to a pair (φ^0,φ1)( _0, _1) that is unique in a projectivized sense, thereby so is the pair (φ^(t,⋅),φ(t,⋅)) ( (t,·), (t,·) ) ∀t∈[t0,t1]∀ t∈[t_0,t_1]. Here “unique in a projectivized sense” means unique up to reciprocal scaling by any constant c>0c>0, i.e., if (φ^0,φ1)( _0, _1) is a solution for (24), then so is (cφ^0,φ1/c)(c _0, _1/c) for any c>0c>0. Thus, the computed pair (cφ^(t,⋅),φ(t,⋅)/c) (c (t,·), (t,·)/c ) is unique up to the choice of c>0c>0. The numerical value of the reciprocal scaling constant c>0c>0 for a converged pair depends on the choice of the initial guess in recursion (24). However, this choice does not affect the unique computation of the original variables of interest, viz. ρoptρ u_opt and opt u_opt, since their recovery from (22)-(23) remains unaffected by reciprocal scaling. Algorithm 1 Sinkhorn with memory 1:Numerical tolerance ε , maximum number of iterations maxiter, positive initial guesses φinit,φ^init _ init, _ init 2:(φ1,φ^0)←(φinit,φ^init) ( _1, _0 )← ( _ init, _ init ) 3:idx←1 idx← 1 ⊳ Initialize recursion index 4:while (err>ε err> ) and (idx<maxiter)( idx< maxiter) do 5: φ1old←φ1 _1 old← _1 6: φ^0old←φ^0 _0 old← _0 7: φ0←ℬ(φ1) _0 ( _1) ⊳ Solve the backward PDE IVP (19b) 8: φ←(φt)t∈[0,1] ←( _t)_t∈[0,1] ⊳ Store history 9: φ^0←ρ0/φ0 _0← _0/ _0 10: φ^1←ℱφ(φ^0) _1 _ ( _0) ⊳ Solve the forward PDE IVP (19a) 11: φ1←ρ1/φ^1 _1← _1/ _1 12: err←maxdH(φ^0,φ^0old),dH(φ1,φ1old) err← \d_H( _0, _0 old),d_H( _1, _1 old)\ 13: idx←idx+1 idx← idx+1 14:end while 15:ρopt←φ^⋅φρ u_opt← · 16:opt←λ⊤∇logφ u_opt←λ g ∇ Result: Optimal solution ρoptρ u_opt, opt u_opt 3.3 The Case ⊤∝̸ g : Dynamic Sinkhorn Recursion with Memory From the Backward Pass When λ⊤−≠λ g - ≠ 0, the terms (20) and (21) are nonzero in general, and the recursion (24) does not apply because then (19a)-(19b) are neither linear nor decoupled. To make algorithmic headway in this more general case, we notice that the BVP (19) still has some structure. In particular, the PDE (19b) is nonlinear yet decoupled from (19a), i.e., the coupling among (19a)-(19b) is one way via the additional terms φ,qφ f_ ,q_ . This motivates our idea of designing a variant of Sinkhorn recursion where starting with an initial guess φ1 _1, one could backward integrate (19b), i.e., evaluate the map ℬB in (24), and store the intermediates (φt)t∈[t0,t1] ( _t )_t∈[t_0,t_1] in a buffer/temporary memory to be accessed for evaluating the map ℱF in (24) that follows. Here and in the sequel, we follow the notational convention that φt(⋅):=φ(t,⋅) _t(·):= (t,·) ∀t∈[t0,t1]∀ t∈[t_0,t_1]. In other words, for solving the IVP with (19a), the most recent (φt)t∈[t0,t1] ( _t )_t∈[t_0,t_1] from the backward pass could be substituted in (19a). Then the resulting linear reaction-advection-diffusion PDE solution can be marched forward in time. This changes the standard Sinkhorn recursion (24) to the new variant: (25) where the box depicts the buffer/temporary memory. In (25), the subscript of the forward-in-time map ℱφF_ emphasizes its dependence on the most recent history of the backward pass. We refer to the recursion (25) as “Sinkhorn with memory”. After each epoch, the buffer is flushed, and the (φt)t∈[t0,t1] ( _t )_t∈[t_0,t_1] is overwritten with the updated evaluation of ℬB. In the following, we concretize this idea into an algorithm. 3.4 Proposed Algorithm Building on the ideas from the previous subsection, we propose Algorithm 1 (Sinkhorn with memory) to solve the generic CASBP (1)-(2) with input and noise channel mismatch. Given the CASBP data q,,,,ρ0,ρ1q, f, g, σ, _0, _1, as in standard Sinkhorn recursion, Algorithm 1 takes three additional inputs: numerical tolerance ε , maximum number of iterations maxiter, and initial guess φinit _ init (everywhere positive). Remark 3 (Symmetry-breaking). Unlike the standard dynamic Sinkhorn recursion (24), there is an asymmetry in the proposed variant (25) in the sense that the proposed recursion must start from the initial guess φ1 _1 (the φinit _ init in Algorithm 1). In contrast, the recursion (24) can either start from the initial guess φ^0 _0 or from φ1 _1. This symmetry-breaking originates from that (19b) does not depend on φ^(t,⋅) (t,·) while (19a) depends on both φ^(t,⋅),φ(t,⋅) (t,·), (t,·). When the channels match, (19a)-(19b) become equation-level decoupled, i.e., then they are only coupled through the boundary condition (19c), thereby allowing the initial guess to be either φ^0 _0 or φ1 _1. In the next Section, we focus on the analysis of Algorithm 1. 4 Analysis of Dynamic Sinkhorn Recursion with Memory from the Backward Pass Given an initial guess, say φinit=φ1(1) _ init= _1^(1), Algorithm 1 is functionally a recursion of the form φ1(j+1)=j(φ1(j)),j=1,…,maxiter, _1^(j+1)=T_j( _1^(j))\>, j=1,…, maxiter, (26) where each mapping j:K+→K+T_j:K_+→ K_+, with K+K_+ as in (5). Unlike the standard dynamic Sinkhorn recursion (24), the map jT_j is not fixed between iterations, which is why we adjoin a subscript j indicating the iteration index. Indeed, ℱφF_ in (25) is dependent on the most recent trajectory φ:=(φt)t∈[0,1] :=( _t)_t∈[0,1], which changes from an iteration to the next. For convenience, we introduce the following notations for the spatial differential operators corresponding to the backward and forward PDEs (19b) and (19a), respectively: ℒ←[φ] \!\!\! L[ ] :=−⟨∇φ,+φ⟩−12⟨,Hessφ⟩+(qλ+qφ)φ, :=- _ x , f+ f_ - 12 ,Hess_ x + ( qλ+q_ ) , (27) ℒ→φ[φ^] \!\!\! L_ [ ] :=−∇⋅(φ^(+φ))+12Δφ^−(qλ+qφ)φ^, :=- _ x· ( ( f+ f_ ) )+ 12 _ \> - ( qλ+q_ ) , (28) where ℒ→φ[φ^] L_ [ ] is meant as the operator for a fixed φ . We also define the division operators i(⋅):=ρi/(⋅),i∈0,1. _i(·):= _i/(·), i∈\0,1\. (29) For a fixed trajectory φ , we define a composite map φT_ as φ(φ1) _ ( _1) :=1(ℱφ(0(ℬ(φ1)))). :=D_1(F_ (D_0(B( _1)))). (30) Thus, at a fixed iteration index j∈1,…,maxiterj∈\1,…, maxiter\, j=φ(j). _j=T_ ^(j). Standard dynamic Sinkhorn recursions, such as (24), are proven [29] to be globally contractive (Definition 2) with respect to Hilbert’s projective metric dH(⋅,⋅)d_H(·,·). Global contractivity leads to global convergence of iterations to a unique fixed point as highlighted earlier in Sec. 3.2. It is then natural to ask what performance guarantee can be established for Algorithm 1. To address this question, we investigate the nature of the constituent maps of jT_j to analyze the proposed algorithm with respect to dH(⋅,⋅)d_H(·,·). For the purpose of analyzing Algorithm 1, we will consider a compact spatial domain D⊂ℝnD ^n. The roadmap for our analysis is as follows. In Sec. 4.1, we discuss non-expansiveness for the maps 0,1,ℬD_0,D_1,B. In Sec. 4.2, we establish sufficient conditions to ensure the map ℱφF_ is locally contractive. In Sec. 4.3, we bring these results together to prove local stability for Algorithm 1 (Theorem 4). 4.1 Nonexpansiveness of 0,1,ℬD_0,D_1,B Using Definition 7, we note that for i∈0,1i∈\0,1\ and for all u,v∈K+u,v∈ K_+, dH(i(u),i(v)) d_ H(D_i(u),D_i(v)) =logsupρi()/u()ρi()/v()−loginfρi()/u()ρi()/v() =\! _ x _i( x)/u( x) _i( x)/v( x)\!-\! _ x _i( x)/u( x) _i( x)/v( x) =dH(u,v). =d_ H(u,v). Hence, the division maps 0D_0 and 1D_1 are in fact isometries with respect to dH(⋅,⋅)d_ H(·,·). The next result establishes that the map ℬB is non-expansive under a sufficient condition. Proposition 1 (Non-expansiveness of ℬB). If qφ+q/λ≥0q_ +q/λ≥ 0, then the backward-in-time map ℬ:K+→K+B:K_+→ K_+ is globally nonexpansive with respect to the Hilbert metric. That is, there exists a constant 0≤κℬ≤10≤ _B≤ 1 such that dH(ℬ(u),ℬ(v))≤κℬdH(u,v)∀u,v∈K+. d_ H(B(u),B(v))≤ _Bd_ H(u,v) ∀ u,v∈ K_+. (31) Proof. Our proof strategy is to verify that the map ℬB is 11-homogeneous (Definition 1) and monotone increasing. We then invoke the Birkhoff-Bushell theorem (Theorem 1 in Sec. 2.1) that establishes nonexpansiveness of ℬB. For α>0α>0, let αφ f_α :=(λ⊤−)∇log(αφ), := (λ g - ) _ x (α ), qαφ q_α :=12(∇log(αφ))⊤(λ⊤−)∇log(αφ). := 12 ( _ x (α ) ) (λ g - ) _ x (α ). Since ∇log(αφ)=∇logφ _ x (α )= _ x , direct computation gives αφ=φ,qαφ=qφ. f_α = f_ , q_α =q_ . Thus, from (27), ℒ←[αφ] L[α ] =−α⟨∇φ,+φ⟩−α2⟨,Hessφ⟩+(qλ+qφ)(αφ) =-α _ x , f+ f_ - α2 ,Hess_ x + ( qλ+q_ )(α ) =αℒ←[φ]. =α L[ ]. That is, the differential operator that corresponds to ℬB, ℒ←[αφ] L[α ], is homogeneous of degree 11. Therefore, if (φt)t∈[t0,t1]( _t)_t∈[t_0,t_1] is the solution to the backward PDE (19b) subject to φ1 _1, then α(φt)t∈[t0,t1]α( _t)_t∈[t_0,t_1] is also the solution to the backward PDE (19b) subject to αφ1α _1. In turn, ℬ(αφ1)=αℬ(φ1) (α _1)= ( _1) as ℬ(φ1)=φ0B( _1)= _0, establishing that ℬB is 11-homogeneous. To prove the monotonicity of ℬB, we make use of the comparison principle for quasilinear parabolic PDEs (Theorem 2). To this end, we define the operator P[η]:=−∂tη−ℒ←[η], P[η]:=- _tη- L[η], for η∈K+η∈ K_+. Using (27), P[η]P[η] can be expressed as P[η] P[η] =−∂tη+∑i,j=1naij(t,,η,∇η)∂xixjη(t,) =- _tη+ _i,j=1^na_ij(t, x,η, _ xη) _x_ix_jη(t, x) +b(t,,η,∇η), +b(t, x,η, _ xη), (32) where aij(t,,η,∇η) a_ij(t, x,η, _ xη) =12Σi,j, = 12 _i,j, b(t,,η,∇η) b(t, x,η, _ xη) =⟨∇η,+η⟩−(qλ+qη)η, = _ xη, f+ f_η - ( qλ+q_η )η, with η f_η :=(λ⊤−)∇logη, := (λ g - ) _ x η, qη q_η :=12(∇logη)⊤(λ⊤−)∇logη. := 12 ( _ x η ) (λ g - ) _ x η. Note that the coefficients aija_ij are independent of η. Further, ∂ηb _ηb =∂η[⟨∇η,+η⟩−(qλ+qη)η] = _η [ _ xη, f+ f_η - ( qλ+q_η )η ] =−12(∇η)⊤(λ⊤−)(∇η)η2−qλ =- 12 ( _ xη) (λ g - )( _ xη)η^2- qλ =−qη−qλ, =-q_η- qλ, (33) so b is decreasing in η if qη+q/λ≥0q_η+q/λ≥ 0. In this case, the comparison principle in Theorem 2 applies to (4.1). Accordingly, we take η(t,⋅)=φt1−tη(t,·)= _t_1-t and ζ(t,⋅)=φt1−t′ζ(t,·)= _t_1-t to be the solutions of the time reversal of (19b), governed by the parabolic operators defined similarly to that in (4.1) and subject to initial conditions η(t0,⋅)=φ1(⋅)η(t_0,·)= _1(·) and ζ(t0,⋅)=φ1′(⋅)ζ(t_0,·)= _1 (·), respectively. If φ1≤φ1′ _1≤ _1 , then by the comparison principle, φt1−t=η(t,⋅)≤ζ(t,⋅)=φt1−t′∀t∈[t0,t1]. _t_1-t=η(t,·)≤ζ(t,·)= _t_1-t ∀ t∈[t_0,t_1]. As ℬ(φ1)=η(t1,⋅)B( _1)=η(t_1,·) and ℬ(φ1′)=ζ(t1,⋅)B( _1 )=ζ(t_1,·), it follows that ℬB is monotone increasing. Having shown that ℬB is 11-homogeneous and monotone increasing, by Theorem 1), we conclude that ℬB is nonexpansive with respect to dH(⋅,⋅)d_H(·,·). ∎ The following corollary offers a simpler sufficient condition in terms of the problem data. Corollary 3 (Positive semi-definiteness of λ⊤−λ g - ). If λ⊤−⪰λ g - 0, then qφ+q/λ≥0q_ +q/λ≥ 0, and Proposition 1 applies. Proof. Since the state cost q≥0q≥ 0, the expression (33) is ≤0≤ 0 whenever λ⊤−⪰λ g - 0. In the absence of state cost (q=0q=0), the conditions qφ+q/λ≥0q_ +q/λ≥ 0 and λ⊤−⪰λ g - 0 are equivalent. ∎ 4.2 The Forward Map ℱφF_ Here, we derive sufficient conditions under which the forward mapping ℱφF_ is locally contractive with respect to dH(⋅,⋅)d_ H(·,·) up to an additive residual. In contrast to the backward-in-time map ℬB, the map ℱφF_ intrinsically depends on the trajectory φ generated from the current iteration. Therefore, when we compare, say ℱφ(⋅)F_ (·) and ℱφ′(⋅)F_ (·), we need to account for the change in inputs of these maps, together with the change in the trajectories φ,φ′ , that underlie their construction. This is formalized in the following triangle inequality: dH(ℱφ(u), d_H(F_ (u), ℱφ′(v))≤dH(ℱφ(u),ℱφ(v)) _ (v))≤ d_H(F_ (u),F_ (v)) +dH(ℱφ(v),ℱφ′(v))∀u,v∈K+. +d_H(F_ (v),F_ (v)) ∀ u,v∈ K_+. (34) We provide two complementary statements (Propositions 2 and 3), each addressing how the forward mapping changes when its two dependencies change separately. Based on them, we point out how the mapping ℱφF_ behaves locally under changing of both dependencies in Proposition 4. Proposition 2 (Contraction of ℱφF_ with φ fixed). For φ fixed, the forward-in-time map ℱφF_ is globally contractive in the Hilbert metric. Proof. With φ fixed, (19a) is a linear second-order parabolic PDE with differential operator ℒ→φ[φ^] L_ [ ] given as in (28). Applying Lemma 1 in [1], we verify that the operator coefficients ∇⋅−( _ x· -( f +φ), + f_ ), 12⟨Hess,⟩ 12 Hess_ x, −∇⋅(+φ)−(qλ+qφ) - _ x·( f+ f_ )- ( qλ+q_ ) are smooth and bounded (due to A1-A3). Then, for fixed φ , (19a) is a linear parabolic PDE with smooth and bounded coefficients, and thereby, it admits [31, 32] a classical fundamental solution Γφ(t,,s,) _ (t, x,s, y). As such, given an initial condition φ^0 _0 supported on a compact D⊂ℝnD ^n, φ^t()=∫DΓφ(t,,t0,)φ^0()d. _t( x)= _D _ (t, x,t_0, y) _0( y) d y. (35) In particular, ℱφ(φ^0)()=φ^1()=∫DΓφ(t1,,t0,)φ^0()d. _ ( _0)( x)= _1( x)= _D _ (t_1, x,t_0, y) _0( y) d y. (36) Theorem 7 in [31] establishes that there exist finite constants C,α1,α2>0C, _1, _2>0 such that C−1ψ1(t1−t0 \!\!C^-1 _1(t_1-t_0 ,−)≤ , x- y)≤ Γφ(t1,,t0,) _ (t_1, x,t_0, y) ≤Cψ2(t1−t0,−)∀,∈D, ≤ C _2(t_1-t_0, x- y) ∀ x, y∈ D, (37) where ψi(t,) _i(t, x) is the fundamental solution of the heat equation ∂tψi=αiΔψi∈1,2. _t _i= _i _i i∈\1,2\. Let m m :=inf,∈DC−1ψ1(t1−t0,−), := _ x, y∈ DC^-1 _1(t_1-t_0, x- y), M M :=sup,∈DCψ2(t1−t0,−). := _ x, y∈ DC _2(t_1-t_0, x- y). Then, m≤Γφ(t1,,t0,)≤M∀,∈D. m≤ _ (t_1, x,t_0, y)≤ M ∀ x, y∈ D. (38) Since 0<ψi<∞0< _i<∞ ∀i∈1,2∀ i∈\1,2\, we have that 0<m≤M<∞0<m≤ M<∞, constituting uniform bounds on the fundamental solution Γφ(t1,,t0,) _ (t_1, x,t_0, y). For any u∈K+u∈ K_+, by (36) and (38), m∫Du()d≤ℱφ(u)()≤M∫Du()d∀∈D. \!\!m _Du( y) d y _ (u)( x)≤ M _Du( y) d y ∀ x∈ D. (39) Letting 1 to be the 1-valued function on D, we bound the projective diameter (see (11)) of ℱφF_ as ProjDiam(ℱφ) ProjDiam(F_ ) ≤2supu∈K+dH(ℱφ(u),) ≤ 2 _u∈ K_+d_ H(F_ (u), 1) =2supu∈K+log(sup∈Dℱφ(u)()inf∈Dℱφ(u)()) =2 _u∈ K_+ ( _ x∈ DF_ (u)( x) _ x∈ DF_ (u)( x) ) ≤2supu∈K+log(M∫Du()dm∫Du()d) ≤ 2 _u∈ K_+ ( M _Du( y) d ym _Du( y) d y ) =2log(Mm)<∞. =2 ( Mm )<∞. Then, by the Birkhoff-Bushell theorem (Theorem 1), we have κH(ℱφ)≤tanh(log(M/m)2)<1. _ H(F_ )≤ ( (M/m)2 )<1. Therefore, ℱφF_ , for a fixed φ , is globally contractive. ∎ Before we proceed to our next statement, we define the operator Γφt→s⋆u _ ^t→ s u for t0≤t≤s≤t1t_0≤ t≤ s≤ t_1 as (Γφt→s⋆u)() ( _ ^t→ s u )( x) :=∫DΓφ(s,,t,)u()d, := _D _ (s, x,t, y)u( y) d y, (40) where Γφ(t,,t,)=δ(−) _ (t, x,t, y)=δ( x- y), the Dirac delta. We remind the reader of the notational convention from Sec. 3 that φ^t,φt _t, _t refer to a pair of snapshots at t∈[t0,t1]t∈[t_0,t_1], while φ^,φ , (without the subscript t) refer to a pair of trajectories. Proposition 3 (Bounding the Hilbert metric between ℱφ(v),ℱφ′(v)F_ (v),F_ (v)). Let φ∗ ^* be a solution of (19b) subject to some final condition φ1∗∈K+ _1^*∈ K_+ and Vφ∗V_ ^* be a neighborhood of φ∗ ^*. Assume the following for all φ,φ′∈Vφ∗ , ∈ V_ ^*. B1 There exist constants mV,MVm_V,M_V where 0<mV≤MV<∞0<m_V≤ M_V<∞ such that mV≤Γφ(t1,,t0,)≤MV∀,∈D. m_V≤ _ (t_1, x,t_0, y)≤ M_V ∀ x, y∈ D. (41) B2 There exists a constant 0≤ℓ<∞0≤ <∞ such that for all t0≤t≤s≤t1t_0≤ t≤ s≤ t_1, ‖Γφt→s⋆u‖∞≤ℓ‖u‖∞∀u∈L∞(D). \| _ ^t→ s u\|_∞≤ \|u\|_∞ ∀ u∈ L^∞(D). (42) B3 There exist constants 0≤ci<∞0≤ c_i<∞, i∈1,2i∈\1,2\ such that supt∈[t0,t1]‖∇ilogφt′−∇ilogφt‖∞ _t∈[t_0,t_1]\|∇^i_ x _t-∇^i_ x _t\|_∞ ≤ci, ≤ c_i, (43) where ∇2:=⟨∇,∇⟩ _ x^2:= _ x, _ x , the standard Laplacian. B4 There exists a constant 0≤c3<∞0≤ c_3<∞ such that supt∈[t0,t1]‖∇logφt‖∞ _t∈[t_0,t_1]\| _ x _t\|_∞ ≤c3. ≤ c_3. (44) B5 For every v∈K+v∈ K_+, there exist constants 0<Ai=Ai(v)<∞0<A_i=A_i(v)<∞, i∈0,1i∈\0,1\, such that for every φ∈Vφ∗ ∈ V_ ^*, if φ is the solution to (19a) that corresponds to φ and subject to φ^0=v _0=v, then supt∈[t0,t1]‖∇iφ^t‖∞ _t∈[t_0,t_1]\|∇^i_ x _t\|_∞ ≤Ai‖v‖∞. ≤ A_i\|v\|_∞. (45) Now let Z:=c1‖∇⋅(λ⊤−)‖∞+(c2+c1c3)‖λ⊤−‖∞Z:=c_1\| _ x· (λ g - )\|_∞+(c_2+c_1c_3)\|λ g - \|_∞ and let αℱ(v):=ℓ(A1c1‖λ⊤−‖∞+A0Z)mV∫Dv()d(t1−t0)‖v‖∞, _F(v):= (A_1c_1\|λ g - \|_∞+A_0Z)m_V _Dv( y) d y(t_1-t_0)\|v\|_∞, (46) for all v∈K+v∈ K_+. Then, for every v∈K+v∈ K_+ and φ,φ′∈Vφ∗ , ∈ V_ ^*, we have dH(ℱφ(v),ℱφ′(v))≤log(1+αℱ(v)1−αℱ(v)),d_H(F_ (v),F_ (v))≤ ( 1+ _F(v)1- _F(v) ), (47) provided that αℱ(v)<1 _F(v)<1. Proof. Let φ,φ′∈Vφ∗ , ∈ V_ ^*. Also let φ^,φ^′ , be solutions of (19a) that correspond to φ and φ′ , respectively, subject to the same initial condition φ^0=φ^0′=v _0= _0 =v, v∈K+v∈ K_+. Letting zt:=φ^t−φ^t′z_t:= _t- _t , and accordingly the trajectories z:=φ^−φ^′z:= - , observe that ∂tz _tz =∂tφ^−∂tφ^′ = _t - _t =ℒ→φ[φ^]−ℒ→φ′[φ^′] = L_ [ ]- L_ [ ] =ℒ→φ[z]+(ℒ→φ[φ^′]−ℒ→φ′[φ^′])⏟=:ℒΔ[φ^′]. = L_ [z]+ ( L_ [ ]- L_ [ ] )_=:L_ [ ]. The difference operator ℒΔ[⋅]L_ [·] reads ℒΔ[φ^′]=⟨∇φ^′,bΔ⟩+φ^′cΔ _ [ ]= _ x ,b_ + c_ (48a) where bΔ b_ :=φ′−φ, := f_ - f_ , (48b) cΔ c_ :=∇⋅(φ′−φ)+(qφ′−qφ). := _ x·( f_ - f_ )+(q_ -q_ ). (48c) By the nonautonomous variation-of-constants formula, also known as Duhamel’s principle [33, Pg. 51][34], it follows that zs=(Γφt0→s⋆zt0)+∫t0s(Γφt→s⋆ℒΔ[φ^t′])dt.z_s= ( _ ^t_0→ s z_t_0 )+ _t_0^s ( _ ^t→ s _ [ _t ] ) dt. In our case zt0=0z_t_0=0, and so zt1=∫t0t1(Γφt→t1⋆ℒΔ[φ^t′])dt.z_t_1= _t_0^t_1 ( _ ^t→ t_1 _ [ _t ] ) dt. Further, ‖zt1‖∞ \|z_t_1\|_∞ ≤∫t0t1‖Γφt→t1⋆ℒΔ[φ^t′]‖∞dt. ≤ _t_0^t_1\| _ ^t→ t_1 _ [ _t ]\|_∞ dt. Due to B3-B5, ℒΔ[φ^t′]∈L∞(D)L_ [ _t]∈ L^∞(D) for all t∈[t0,t1]t∈[t_0,t_1]. Then by B2, we have ‖zt1‖∞ \|z_t_1\|_∞ ≤ℓ∫t0t1‖ℒΔ[φ^t′]‖∞dt ≤ _t_0^t_1\|L_ [ _t ]\|_∞ dt ≤ℓ(t1−t0)supt∈[t0,t1]‖ℒΔ[φ^t′]‖∞. ≤ (t_1-t_0) _t∈[t_0,t_1]\|L_ [ _t ]\|_∞. (49) Now observe that, supt∈[t0,t1]‖bΔ(t)‖∞ _t∈[t_0,t_1]\|b_ (t)\|_∞ =(48b) eq:b-delta= supt∈[t0,t1]‖φ′(t)−φ(t)‖∞ _t∈[t_0,t_1]\| f_ (t)- f_ (t)\|_∞ =(20) def_fphi= supt∈[t0,t1]‖(λ⊤−)(∇logφt′−∇logφt)‖∞ _t∈[t_0,t_1]\| (λ g - )( _ x _t- _ x _t)\|_∞ ≤B3 B2≤ c1‖λ⊤−‖∞. c_1\|λ g - \|_∞. (50) We also have ‖cΔ(t)‖∞=(48c)‖∇⋅(φ′−φ)(t)+(qφ′−qφ)(t)‖∞ \|c_ (t)\|_∞ eq:c-delta=\| _ x·( f_ - f_ )(t)+(q_ -q_ )(t)\|_∞ ≤‖∇⋅(φ′−φ)(t)‖∞+‖(qφ′−qφ)(t)‖∞. ≤\| _ x·( f_ - f_ )(t)\|_∞+\|(q_ -q_ )(t)\|_∞. (51) Note that ∥ \| ∇⋅(φ′−φ)(t)∥∞≤(20) _ x·( f_ - f_ )(t)\|_∞ def_fphi≤ ‖(∇⋅λ⊤−)⋅(∇logφt′−∇logφt)‖∞ \|( _ x·λ g - )·( _ x _t - _ x _t)\|_∞ +‖(λ⊤−)(∇2logφt′−∇2logφt)‖∞, +\| (λ g - )(∇^2_ x _t-∇^2_ x _t)\|_∞, ≤B3c1‖∇⋅λ⊤−‖∞+c2‖λ⊤−‖∞, B2≤c_1\| _ x·λ g - \|_∞+c_2\|λ g - \|_∞, (52) and ∥(qφ′ \|(q_ −qφ)(t)∥∞≤(21)[12∥∇logφt′+∇logφt∥∞× -q_ )(t)\|_∞ def_qphi≤ [ 12\| _ x _t+ _ x _t\|_∞× ∥λ⊤−∥∞∥∇logφt′−∇logφt∥∞], \|λ g - \|_∞\| _ x _t- _ x _t\|_∞ ], ≤B4c1c3‖λ⊤−‖∞. B3≤c_1c_3\|λ g - \|_∞. (53) Therewith, from (51)-(4.2), we have supt∈[t0,t1]‖cΔ(t)‖∞≤ _t∈[t_0,t_1]\|c_ (t)\|_∞≤ c1∥∇⋅λ⊤−∥∞+(c2+c1c3)∥λ⊤−∥∞=:Z. c_1\| _ x·λ g - \|_∞+(c_2+c_1c_3)\|λ g - \|_∞=:Z. (54) From (48), we know that supt∈[t0,t1]‖ℒΔ[φ^t′]‖∞ _t∈[t_0,t_1]\|L_ [ _t ]\|_∞ ≤supt∈[t0,t1](‖∇φ^t′‖∞‖bΔ(t)‖∞) ≤ _t∈[t_0,t_1](\| _ x _t \|_∞\|b_ (t)\|_∞) +supt∈[t0,t1](‖φ^t′‖∞‖cΔ(t)‖∞). + _t∈[t_0,t_1](\| _t \|_∞\|c_ (t)\|_∞). By using (50) and (54), we obtain supt∈[t0,t1]‖ℒΔ[φ^t′]‖∞ \!\!\! _t∈[t_0,t_1]\|L_ [ _t ]\|_∞ ≤supt∈[t0,t1]‖∇φ^t′‖∞c1‖λ⊤−‖∞ ≤ _t∈[t_0,t_1]\| _ x _t \|_∞c_1\|λ g - \|_∞ +supt∈[t0,t1]‖φ^t′‖∞Z + _t∈[t_0,t_1]\| _t \|_∞Z ≤B5(A1c1‖λ⊤−‖∞+A0Z)‖v‖∞. -20.0pt B4≤(A_1c_1\|λ g - \|_∞+A_0Z)\|v\|_∞. (55) By substituting the previous inequality into (49), we have ‖zt1‖∞≤ℓ(t1−t0)(A1c1‖λ⊤−‖∞+A0Z)‖v‖∞. \|z_t_1\|_∞≤ (t_1-t_0)(A_1c_1\|λ g - \|_∞+A_0Z)\|v\|_∞. (56) Now observe that −‖zt1‖∞inf∈Dφ^1′()≤zt1()φ^1′()≤‖zt1‖∞inf∈Dφ^1′(), -\|z_t_1\|_∞ _ x∈ D _1 ( x)≤ z_t_1( x) _1 ( x)≤ \|z_t_1\|_∞ _ x∈ D _1 ( x), (57) for all ∈D x∈ D and further, inf∈Dφ^1′() _ x∈ D _1 ( x) =inf∈D∫Γφ′(t1,,t0,)v()d = _ x∈ D _ (t_1, x,t_0, y)v( y) d y ≥B1mV∫Dv(). B1≥m_V _Dv( y)d y. (58) By substituting (56) and (4.2) into (57), we get −αℱ(v)≤zt1()φ^1′()≤αℱ(v)∀∈D - _F(v)≤ z_t_1( x) _1 ( x)≤ _F(v) ∀ x∈ D (59) where αℱ(v) _F(v) is given in (46). Recall that zt1=φ^1−φ^1′z_t_1= _1- _1 . Consequently, φ^1φ^1′=1+zt1φ^1′. _1 _1 =1+ z_t_1 _1 . (60) By (59) and (60), we obtain 1−αℱ(v)≤φ^1φ^1′=ℱφ(v)ℱφ′(v)≤1+αℱ(v). 1- _F(v)≤ _1 _1 = F_ (v)F_ (v)≤ 1+ _F(v). By definition of the Hilbert metric (Definition 7), dH(ℱφ(v),ℱφ′(v))≤log(1+αℱ)−log(1−αℱ), d_ H(F_ (v),F_ (v))≤ (1+ _F)- (1- _F), (61) resulting in the expression in (47). ∎ Proposition 4 (Bounding the Hilbert metric between ℱφ(u),ℱφ′(v)F_ (u),F_ (v)). Consider the neighborhood Vφ∗V_ ^* and the assumptions in Proposition 3. Then for every φ,φ′∈Vφ∗ , ∈ V_ ^* and u,v∈K+u,v∈ K_+, dH(ℱφ(u),ℱφ′(v)) d_H(F_ (u),F_ (v)) ≤κℱ∗dH(u,v)+ε(v), ≤ _F^*\,d_H(u,v)+ (v), (62) where κℱ∗ _F^* :=supφ∈Vφ∗κH(ℱφ), := _ ∈ V_ ^* _ H(F_ ), (63) ε(v) (v) :=log(1+αℱ(v)1−αℱ(v)), := ( 1+ _F(v)1- _F(v) ), (64) provided that αℱ(v)<1 _F(v)<1. In addition, κℱ∗≤tanh(log(MV/mV)2)<1. _F^*≤ ( (M_V/m_V)2 )<1. Proof. From Proposition 2, we know that for a fixed φ , ℱφF_ is contractive with contraction ratio κH(ℱφ) _H(F_ ) satisfying κH(ℱφ)≤tanh(log(M/m)2)<1. _ H(F_ )≤ ( (M/m)2 )<1. (65) Similar arguments to those in the proof of Proposition 2 together with B1 imply that κH(ℱφ)≤tanh(log(MV/mV)2)<1 _ H(F_ )≤ ( (M_V/m_V)2 )<1 (66) for all φ∈Vφ∗ ∈ V_ ^*. Using the definition in (63), then it holds that κℱ∗≤tanh(log(MV/mV)2)<1. _F^*≤ ( (M_V/m_V)2 )<1. (67) Consequently, for all φ∈Vφ∗ ∈ V_ ^*, it follows that dH(ℱφ(u),ℱφ(v))≤κℱ∗dH(u,v). d_H(F_ (u),F_ (v))≤ _F^*\,d_H(u,v). (68) Also immediately from Proposition 3, we get that for all φ,φ′∈Vφ∗,v∈K+ , ∈ V_ ^*,v∈ K_+, dH(ℱφ(v),ℱφ′(v))≤ε(v), d_H(F_ (v),F_ (v))≤ (v), (69) with ε(v) (v) given by (64). By the triangle inequality dH(ℱφ(u),ℱφ′(v)) d_H(F_ (u),F_ (v)) ≤dH(ℱφ(u),ℱφ(v)) ≤ d_H(F_ (u),F_ (v)) +dH(ℱφ(v),ℱφ′(v)), +d_H(F_ (v),F_ (v)), and the inequalities (68) and (69), we arrive at (62). ∎ 4.3 Local Stability of Algorithm 1 Building on Propositions 1-4, our main result on the stability of Algorithm 1 is as follows. Theorem 4. Consider φ1∗∈K+ _1^*∈ K_+ such that φ1∗=φ∗(φ1∗), _1^*=T_ ^*( _1^*), (70) i.e., φ1∗ _1^* is a fixed point of the update map φ∗T_ ^* defined by the trajectory φ∗ ^*. Let φ^0∗:=0(ℬ(φ1∗)). _0^*:=D_0(B( _1^*)). For r>0r>0, consider a Hilbert ball centered around φ1∗ _1^*, given by Br(φ1∗):=φ1∈K+|dH(φ1,φ1∗)≤r, B_r( _1^*):=\ _1∈ K_+~|~d_ H( _1, _1^*)≤ r\, (71) such that whenever φ1∈Br(φ1∗) _1∈ B_r( _1^*), its corresponding trajectory φ belongs to Vφ∗V_ ^*. Under the assumptions of Propositions 1-4, if r≥ε(φ^0∗)1−κℬκℱ∗, r≥ ( _0^*)1- _B _F^*, (72) and the initial guess of Algorithm 1, φ1(1)=φinit _1^(1)= _ init, belongs to Br(φ1∗)B_r( _1^*), then the updates φ1(j+1)=j(φ1(j)) _1^(j+1)=T_j( _1^(j)) also belong to Br(φ1∗)B_r( _1^*) for all iteration index j∈ℕj . Additionally, limj→∞dH(φ1(j),φ1∗)≤ε(φ^0∗)1−κℬκℱ∗. _j→∞d_ H( _1^(j), _1^*)≤ ( _0^*)1- _B _F^*. (73) Proof. Let φ1∈Br(φ1∗) _1∈ B_r( _1^*) and the corresponding trajectory φ∈Vφ∗ ∈ V_ ^*. Also, let φ^0:=0(ℬ(φ1)) _0:=D_0(B( _1)). Since 0:K+→K+D_0:K_+→ K_+ is an isometry with respect to dH(⋅,⋅)d_ H(·,·), we have dH(φ^0,φ^0∗)=dH(ℬ(φ1),ℬ(φ1∗)). d_ H( _0, _0^*)=d_ H(B( _1),B( _1^*)). (74) Under the assumption of Proposition 1, ℬB is non-expansive, thus satisfying dH(ℬ(φ1),ℬ(φ1∗))≤κℬdH(φ1,φ1∗), d_ H(B( _1),B( _1^*))≤ _Bd_ H( _1, _1^*), (75) with κℬ≤1 _B≤ 1. Then, by (74) and (75), we obtain dH(φ^0,φ^0∗)≤κℬdH(φ1,φ1∗). d_ H( _0, _0^*)≤ _Bd_ H( _1, _1^*). (76) We also note that dH(φ(φ1),φ1∗) d_ H(T_ ( _1), _1^*) =dH(φ(φ1),φ∗(φ1∗)) =d_ H(T_ ( _1),T_ ^*( _1^*)) =dH(ℱφ(φ^0),ℱφ∗(φ^0∗)), =d_ H(F_ ( _0),F_ ^*( _0^*)), (77) because 1:K+→K+D_1:K_+→ K_+ is an isometry with respect to dH(⋅,⋅)d_ H(·,·). Under the assumptions of Propositions 2-4, we have dH(ℱφ(φ^0),ℱφ∗(φ^0∗))≤κℱ∗dH(φ^0,φ^0∗)+ε(φ^0∗). d_ H(F_ ( _0),F_ ^*( _0^*))≤ _F^*\,d_ H( _0, _0^*)+ ( _0^*). (78) Then, by (4.3) and (78), we obtain dH(φ(φ1),φ1∗)≤κℱ∗dH(φ^0,φ^0∗)+ε(φ^0∗). d_ H(T_ ( _1), _1^*)≤ _F^*\,d_ H( _0, _0^*)+ ( _0^*). (79) By (76) and (79), we get dH(φ(φ1),φ1∗)≤κℬκℱ∗dH(φ1,φ1∗)+ε(φ^0∗). d_ H(T_ ( _1), _1^*)≤ _B _F^*\,d_ H( _1, _1^*)+ ( _0^*). (80) The previous relation implies that dH(φ1(2),φ1∗) d_ H( _1^(2), _1^*) =dH(φ(1)(φ1(1)),φ1∗) =d_ H(T_ ^(1)( _1^(1)), _1^*) ≤(80)κℬκℱ∗dH(φ1(1),φ1∗)+ε(φ^0∗). eq:dis-T-2≤ _B _F^*\,d_ H( _1^(1), _1^*)+ ( _0^*). If φ1(1) _1^(1) belongs to Br(φ1∗)B_r( _1^*), then dH(φ1(2),φ1∗)≤rκℬκℱ∗+ε(φ^0∗). d_ H( _1^(2), _1^*)≤ r _B _F^*+ ( _0^*). And as long as the condition (72) holds, φ1(2)∈Br(φ1∗) _1^(2)∈ B_r( _1^*). In a similar way, we can verify that dH(φ1(j), d_ H( _1^(j), φ1∗)≤(κℬκℱ∗)j−1dH(φ1(1),φ1∗)+ _1^*)≤( _B _F^*)^j-1d_ H( _1^(1), _1^*)+ ε(φ^0∗)(1+κℬκℱ∗+⋯+(κℬκℱ∗)j−2), ( _0^*)(1+ _B _F^*+·s+( _B _F^*)^j-2), and φ1(j)∈Br(φ∗) _1^(j)∈ B_r( ^*) if the condition (72) holds. Further, limj→∞dH(φ1(j), _j→∞d_ H( _1^(j), φ1∗)≤rlimj→∞(κℬκℱ∗)j−1 _1^*)≤ r _j→∞( _B _F^*)^j-1 +ε(φ^0∗)limj→∞∑i=0j−2(κℬκℱ∗)i. + ( _0^*) _j→∞ _i=0^j-2( _B _F^*)^i. (81) Since κℬκℱ∗<1 _B _F^*<1, limj→∞(κℬκℱ∗)j=0 _j→∞( _B _F^*)^j=0, and the geometric sum for the second summand in (81) yields (73). ∎ Theorem 4 states that Algorithm 1 is locally stable around a fixed point φ1∗ ^*_1. In that, if the initial guess lies within a Hilbert metric ball centered around φ1∗ _1^*, consecutive iterates remain in this ball, provided that its radius satisfies (72). Further, this Theorem shows that Algorithm 1 converges to a Hilbert metric ball around the fixed point, with radius not exceeding ε(φ^0∗)/(1−κℬκℱ∗) ( _0^*)/(1- _B _F^*), thereby providing an upper bound on the asymptotic error for Algorithm 1. This error is sufficiently small, precisely ε(φ^0∗) ( _0^*) is small, for small channel mismatch and spatial rate of change for such mismatch, quantified respectively by ‖λ⊤−‖∞\|λ g - \|_∞ and ‖∇⋅(λ⊤−)‖∞\| _ x· (λ g - )\|_∞. Remark 4 (Warm start for Algorithm 1). The above discussion suggests a warm start or homotopic approach for practical implementation of Algorithm 1. For instance, if Algorithm 1 fails to converge for large channel mismatch, one may take the initial guess φinit _ init to be a prior output (i.e., converged φ1 _1) of the Algorithm, solved then with small channel mismatch. 5 Numerical Example To illustrate the proposed algorithm, we consider steering the state PDF for a cubic spring-mass-damper with process noise. Specifically, we consider an instance of (1) with m=1m=1 input, n=2n=2 states, p=2p=2 noises, time horizon [t0,t1]=[0,1][t_0,t_1]=[0,1], quadratic state cost q=12()⊤q= 12( x u) Q x u with =diag(1,2) Q=diag(1,2), λ=1λ=1, and (t,x1,x2)=(x2−∂x1V(x1)−γx2),V(⋅)=14(⋅)4,γ=1, \! f(t,x_1,x_2)\!=\!\! pmatrixx_2\\ - ∂ x_1V(x_1)-γ x_2 pmatrix\!\!,\;V(·)= 14(·)^4,\;γ=1, (t,x1,x2)=(01),(t,x1,x2)=. g(t,x_1,x_2)= pmatrix0\\ 1 pmatrix,\; σ(t,x_1,x_2)= I. This is an instance of CASBP with channel mismatch because here, ⊤=[0001] g = bmatrix0&0\\ 0&1 bmatrix and = = I are not proportional. Figure 1: The two pairs of endpoint PDFs used in the numerical example in Sec. 5. We solve the above instance of problem (1) for two pairs of endpoint PDFs (Fig. 1): one pair being ρ0()=((0.25−0.25),120), _0( x)=N ( pmatrix0.25\\ -0.25 pmatrix, 120 I ), ρ1()=12((−0.40.4),140)+12((−0.25−0.25),120), _1( x)= 12N ( pmatrix-0.4\\ 0.4 pmatrix, 140 I )+\! 12N ( pmatrix-0.25\\ -0.25 pmatrix, 120 I ), the other pair being ρ0()=14∑i=14(0i,180),∀0i∈−0.25,0.252, _0( x)= 14 _i=1^4N ( μ_0i, 180 I ),\;∀ μ_0i∈\-0.25,0.25\^2, ρ1()=13((00.4),140)+13((−0.3−0.3),140) _1( x)= 13N ( pmatrix0\\ 0.4 pmatrix, 140 I )+\! 13N ( pmatrix-0.3\\ -0.3 pmatrix, 140 I ) +13((0.3−0.3),140), +\! 13N ( pmatrix0.3\\ -0.3 pmatrix, 140 I ), all re-normalized over [−1,1]2[-1,1]^2. (a) For ρ0,ρ1 _0, _1 as in Fig. 1(a). (b) For ρ0,ρ1 _0, _1 as in Fig. 1(b). Figure 2: Convergence of the proposed Algorithm 1 for the numerical example in Sec. 5. For implementing Algorithm 1 to this problem data, we fix the numerical tolerance ε=10−2 =10^-2 and maxiter=200 maxiter=200. We perform the backward time marching for the nonlinear PDE IVP (19b) via a forward-in-time central-in-space (FTCS) finite difference solver with spatial resolution Δx1=Δx2=2×10−2 x_1= x_2=2× 10^-2 for the domain [−1,1]2[-1,1]^2, and temporal resolution Δt=5×10−5 t=5× 10^-5 for the time horizon [0,1][0,1]. We then forward time march the PDE (19a) by substituting the most recent backward pass solution of (19b), and then applying FTCS finite difference with the same resolution. We then update the boundary condition using (19c), and repeat. The resulting convergence with respect to the Hilbert metric dHd_H is shown in Fig. 2. Figure 3: Snapshots of ρoptρ u_opt (top panel) and opt u_opt (bottom panel) for the numerical example in Sec. 5 with endpoint PDFs as in Fig. 1(a). All subplots are over the state domain [−1,1]2[-1,1]^2. Figure 4: Snapshots of ρoptρ u_opt (top panel) and opt u_opt (bottom panel) for the numerical example in Sec. 5 with endpoint PDFs as in Fig. 1(b). All subplots are over the state domain [−1,1]2[-1,1]^2. The corresponding solution pairs (ρopt,opt) (ρ u_opt, u_opt ) are shown in Fig. 3 and Fig. 4. Figure 5: Snapshots of qφ+q/λq_ +q/λ with λ=1λ=1 for the numerical example in Sec. 5 with endpoint PDFs as in Fig. 1(a). All subplots are over the state domain [−1,1]2[-1,1]^2. Fig. 5 shows the snapshots of qφ+q/λq_ +q/λ (here λ=1λ=1) for the converged solution. This particular quantity appears as a sufficient condition for non-expansiveness of ℬB in Proposition 1, and for our numerical example, has the form qφ+q=12(∇logφ)⊤[−1000]∇logφ+12⊤[1002].q_ +q= 12 ( _ x ) bmatrix-1&0\\ 0&0 bmatrix _ x + 12 x bmatrix1&0\\ 0&2 bmatrix x. From Fig. 5, this quantity, during the backward pass from t=1t=1 to t=0t=0, becomes nonnegative after initial transients, i.e., is nonnegative for most but not for all times. This implies there is room for future improvement for our theoretical guarantees, and that the practical performance of the proposed Algorithm 1 is better than the guarantees proved herein. An improved analysis will require new techniques and will be pursued in our future work. 6 Conclusions This work proposes a computational algorithm for solving the generic control-affine Schrödinger bridge problem with input and noise channel mismatch. This algorithm enables computational synthesis of optimal feedback controller for nonlinear non-Gaussian density steering over a fixed finite horizon. The proposed algorithm, referred to as “Sinkhorn with memory,” can be seen as an extension of the dynamic Sinkhorn recursion – the latter is not applicable in the channel mismatch case. We provide a local stability guarantee for the proposed algorithm. This guarantee takes the form of convergence to a Hilbert metric ball around the fixed point with a radius estimate in terms of the channel mismatch magnitude. Numerical examples highlight that the proposed algorithm performs well in practice. References References [1] A. Teter and A. Halder, “On the Hopf-Cole transform for control-affine Schrödinger bridge,” arXiv preprint arXiv:2503.17640, 2025. [2] A. M. H. Teter, A. Halder, M. D. Schneider, A. S. Perloff, J. Pratt, C. M. Artman, and M. Demireva, “Control-affine Schrödinger bridge and generalized Bohm potential,” IEEE Control Systems Letters, vol. 9, p. 2453–2458, 2025. [3] E. Schrödinger, “Über die Umkehrung der Naturgesetze,” Sitzungsberichte der Preuss Akad. Wissen. Phys. Math. Klasse, Sonderausgabe, vol. IX, p. 144–153, 1931. [4] E. Schrödinger, “Sur la théorie relativiste de l’électron et l’interprétation de la mécanique quantique,” in Annales de L’Institut Henri Poincaré, vol. 2, no. 4. Presses universitaires de France, 1932, p. 269–310. [5] A. Wakolbinger, “Schrödinger bridges from 1931 to 1991,” in Proc. of the 4th Latin American Congress in Probability and Mathematical Statistics, Mexico City, 1990, p. 61–79. [6] A. Blaquiere, “Controllability of a Fokker-Planck equation, the Schrödinger system, and a related stochastic optimal control (revised version),” Dynamics and Control, vol. 2, no. 3, p. 235–253, 1992. [7] 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. [8] 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. [9] A. M. Teter, W. Wang, and A. Halder, “Weyl calculus and exactly solvable Schrödinger bridges with quadratic state cost,” in 2024 60th Annual Allerton Conference on Communication, Control, and Computing. IEEE, 2024, p. 1–8. [10] 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. [11] 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. [12] 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. [13] Y. Chen, T. T. Georgiou, and M. Pavon, “Stochastic control liaisons: Richard Sinkhorn meets Gaspard Monge on a Schrödinger bridge,” Siam Review, vol. 63, no. 2, p. 249–313, 2021. [14] 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. [15] I. Nodozi, J. O’Leary, A. Mesbah, and A. Halder, “A physics-informed deep learning approach for minimum effort stochastic control of colloidal self-assembly,” in 2023 American Control Conference (ACC). IEEE, 2023, p. 609–615. [16] 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. [17] J. D. Cole, “On a quasi-linear parabolic equation occurring in aerodynamics,” Quarterly of Applied Mathematics, vol. 9, no. 3, p. 225–236, 1951. [18] 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. [19] Z. Pan and T. Basar, “Backstepping controller design for nonlinear stochastic systems under a risk-sensitive cost criterion,” SIAM Journal on Control and optimization, vol. 37, no. 3, p. 957–995, 1999. [20] W. Li and M. Krstic, “Stochastic nonlinear prescribed-time stabilization and inverse optimality,” IEEE Transactions on Automatic Control, vol. 67, no. 3, p. 1179–1193, 2021. [21] R. Brockett, “Notes on the control of the Liouville equation,” in Control of Partial Differential Equations: Cetraro, Italy 2010, Editors: Piermarco Cannarsa, Jean-Michel Coron. Springer, 2012, p. 101–129. [22] H. J. Kappen, “Linear theory for control of nonlinear stochastic systems,” Physical Review Letters, vol. 95, no. 20, p. 200201, 2005. [23] E. Todorov, “Efficient computation of optimal actions,” Proceedings of the National Academy of Sciences, vol. 106, no. 28, p. 11 478–11 483, 2009. [24] M. B. Horowitz, A. Damle, and J. W. Burdick, “Linear Hamilton Jacobi Bellman equations in high dimensions,” in 53rd IEEE Conference on Decision and Control. IEEE, 2014, p. 5880–5887. [25] W. H. Fleming and H. M. Soner, Controlled Markov processes and viscosity solutions. Springer, 2006. [26] G. Birkhoff, “Extensions of Jentzsch’s theorem,” Transactions of the American Mathematical Society, vol. 85, no. 1, p. 219–227, 1957. [27] P. J. Bushell, “Hilbert’s metric and positive contraction mappings in a Banach space,” Archive for Rational Mechanics and Analysis, vol. 52, no. 4, p. 330–338, 1973. [28] 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. [29] 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. [30] G. M. Lieberman, Second order parabolic differential equations. World scientific, 1996. [31] D. G. Aronson, “Non-negative solutions of linear parabolic equations,” Annali della Scuola Normale Superiore di Pisa-Scienze Fisiche e Matematiche, vol. 22, no. 4, p. 607–694, 1968. [32] D. Aronson, “Non-negative solutions of linear parabolic equations: An addendum,” Annali della Scuola Normale Superiore di Pisa-Scienze Fisiche e Matematiche, vol. 25, no. 2, p. 221–228, 1971. [33] L. C. Evans, Partial differential equations. American mathematical society, 2022, vol. 19. [34] M. Filali and M. Moussi, “Non-autonomous inhomogeneous boundary Cauchy problems and retarded equations,” Proyecciones (Antofagasta), vol. 22, no. 2, p. 145–159, 2003.