1 Introduction↩︎

1.1 Optimizing over Markov processes↩︎

In this paper, we study the problem of finding a Markov process that connects two prescribed distributions while minimizing a given cost. Let \((\Omega,\mathcal{B}(\Omega))\) be a Borel state space. We denote by \(\mathcal{M}(\Omega)\) the set of finite signed Borel measures on \(\Omega\), and by \[\begin{align} \mathcal{P}(\Omega):=\{\rho\in\mathcal{M}(\Omega):\rho\ge 0,\;\rho(\Omega)=1\} \end{align}\] the set of probability measures on \(\Omega\). Given initial and terminal distributions \(\mu,\nu\in\mathcal{P}(\Omega)\), and one-step cost functions \(c^s:\Omega\times\Omega\to\mathbb{R}\) for \(s=0,\ldots,T-1\), we seek intermediate laws \(\rho^s\in\mathcal{P}(\Omega)\) and Markov transition kernels \(K^s(x,dy)\) solving \[\begin{align} \inf_{\left\{ K^s \right\}_{s=0}^{T-1}, \left\{ \rho^s \right\}_{s=0}^T } & \sum_{s=0}^{T-1} \int_\Omega\int_\Omega c^s(x,y)\,K^s(x,\mathop{}\!\mathrm{d}y)\mathop{}\!\mathrm{d}\rho^s(x) \tag{1} \\ \text{s.t.}\quad& \rho^0 = \mu,\quad \rho^T = \nu, \tag{2} \\ & \rho^{s+1} = \left( K^s \right)^* \rho^s, \quad s=0,\ldots,T-1, \tag{3} \\ & K^s \in \mathcal{K}^s, \quad s=0,\ldots,T-1, \tag{4}\\ & \rho^s \ge 0, \quad \rho^s(\Omega) = 1,\quad s=0,\ldots,T. \tag{5} \end{align}\] Here \(\mathcal{K}^s\) denotes the admissible class of Markov kernels at step \(s\). The kernels are allowed to depend on \(s\), so the process is generally time-inhomogeneous. In the Markov evolution constraint 3 , \((K^s)^*\) denotes the adjoint action of the kernel on probability measures: if \(K^s\) acts on test functions by \[\begin{align} K^s f(x)=\int_{\Omega}f(y)\,K^s(x,dy), \end{align}\] then, for every Borel set \(B\in\mathcal{B}(\Omega)\) and probability measure \(\rho\in\mathcal{P}(\Omega)\), \[\begin{align} (K^s)^*\rho(B)=\int_{\Omega}K^s(x,B)\mathop{}\!\mathrm{d}\rho(x). \end{align}\]

A special case of problem 1 is dynamic optimal transport. Classical optimal transport chooses a coupling between two endpoint distributions and minimizes a transportation cost [1], [2]. Dynamic optimal transport instead seeks a curve of probability measures \(\left\{ q(\, \cdot \,, t)\in \mathcal{P}(\Omega): t\in[0,1] \right\}\) with prescribed initial and terminal laws \(q(\, \cdot \,,0)\) and \(q(\, \cdot \,,1)\). When \(\Omega \subset \mathbb{R}^d\), the Benamou–Brenier formula realizes this path through a velocity field [3]: \[\begin{align} \label{problem:BB} \mathcal{W}_p^p(\mu,\nu)=\inf_{q, v}\left\{\int_0^1\int_{\Omega}\|v(x,t)\|^p\mathop{}\!\mathrm{d}q(x,t)\mathop{}\!\mathrm{d}t:\partial_t q +\nabla \cdot(v\,q)=0, \; q(\, \cdot \,, 0)=\mu, \;q(\, \cdot \,, 1)=\nu\right\}. \end{align}\tag{6}\] Here \(\mathcal{W}_p\) denotes the \(p\)-Wasserstein distance and \(\left\{ v(\, \cdot \,, t): \Omega \rightarrow \mathbb{R}^d \mid t \in [0,1] \right\}\) is the velocity field transporting the mass. In 3, we show that, after time discretization, the unconstrained kinetic-cost case of problem 1 exactly recovers the optimal value and the grid-time measures of the Benamou–Brenier problem 6 .

Beyond the unconstrained kinetic-cost case, Problem 1 also describes more general time-discrete distributional dynamics through the admissible kernel classes \(\mathcal{K}^s\). Depending on the application, these classes may encode support, transition-probability, or local-update constraints on the one-step evolution. The single-spin process studied in 4 is one concrete example. We focus on settings in which the corresponding coupling constraints admit tractable convex local outer approximations.

1.2 Sequential-coupling formulation↩︎

The Markov kernel formulation 1 can equivalently be written as a sequential-coupling formulation. Given \(\rho^s\) and \(K^s\), define \[\begin{align} \label{eq:pi323832K} \pi^s(\mathop{}\!\mathrm{d}x,\mathop{}\!\mathrm{d}y) := \rho^s(\mathop{}\!\mathrm{d}x)K^s(x,\mathop{}\!\mathrm{d}y). \end{align}\tag{7}\] Then \(\pi^s\in\mathcal{P}(\Omega\times\Omega)\) is the joint law of two consecutive states \((X^s,X^{s+1})\). Let \({\rm P}_{\rm L},{\rm P}_{\rm R}:\Omega\times\Omega\to\Omega\) be the canonical projections onto the first and second components. The marginals of \(\pi^s\) satisfy \[\begin{align} \label{eq:rho323832pi} \left( {\rm P}_{\rm L} \right)_\sharp\pi^s = \rho^s, \quad \left( {\rm P}_{\rm R} \right)_\sharp\pi^s = \rho^{s+1}, \quad s = 0, \ldots, T-1 \end{align}\tag{8}\] where \(\left( \, \cdot \, \right)_\sharp\) denotes the push-forward operator. Conversely, assuming \(\Omega\) is a standard Borel space, any joint law \(\pi^s\) admits a disintegration with respect to its first marginal, yielding a Markov transition kernel \(K^s\) that satisfies 7 and is unique \(\rho^s\)-almost everywhere. In the finite-state case, this is simply \(K^s(y|x)=\pi^s(x,y)/\rho^s(x)\) whenever \(\rho^s(x)>0\). Therefore, optimizing over the transition kernels \(K^s\) is equivalent to optimizing over the joint couplings \(\pi^s\). Let \(\mathcal{D}^s\) be the sequential-coupling form of the Markov kernel constraints 4 : \[\label{defiD} \mathcal{D}^s:=\left\{ \pi^s\in \mathcal{P}(\Omega\times \Omega):\;\exists\;\rho^s\in \mathcal{P}(\Omega),\;K^s\in \mathcal{K}^s\;{\rm s.t.}\;\pi^s(\mathop{}\!\mathrm{d}x,\mathop{}\!\mathrm{d}y)=\rho^s (\mathop{}\!\mathrm{d}x)K^s(x,\mathop{}\!\mathrm{d}y) \right\},\tag{9}\] which means a nonnegative joint law \(\pi^s\) belongs to \(\mathcal{D}^s\) if and only if its conditional transition kernel belongs to \(\mathcal{K}^s\). Then 1 can equivalently be written as \[\begin{align} \inf_{\{\pi^s\}_{s=0}^{T-1}} \; & \sum_{s=0}^{T-1} \int_{\Omega\times\Omega} c^s(x,y)\mathop{}\!\mathrm{d}\pi^s(x,y) \tag{10}\\ \text{s.t.}\quad& \left( {\rm P}_{\rm L} \right)_\sharp\pi^0 = \mu, \quad \left( {\rm P}_{\rm R} \right)_\sharp\pi^{T-1} = \nu, \tag{11} \\ & \left( {\rm P}_{\rm R} \right)_\sharp\pi^s = \left( {\rm P}_{\rm L} \right)_\sharp\pi^{s+1}, \quad s=0,\ldots,T-2, \tag{12} \\ & \pi^s\in \mathcal{D}^s, \quad s=0,\ldots,T-1, \tag{13} \\ & \pi^s \ge 0, \quad \pi^s(\Omega \times \Omega)=1,\quad s=0,\ldots,T-1.\tag{14} \end{align}\] The intermediate distributions \(\left\{ \rho^s \right\}_{s=0}^{T}\) are absorbed into the couplings: once \(\{\pi^s\}_{s=0}^{T-1}\) is known, the intermediate marginals at the prescribed grid times are recovered as \[\begin{align} \rho^0=\left( {\rm P}_{\rm L} \right)_\sharp\pi^0, \quad \rho^{s+1}=\left( {\rm P}_{\rm R} \right)_\sharp\pi^s, \quad s=0,\ldots,T-1. \end{align}\]

Our central goal is to solve problem 10 in high dimensions. This is challenging because representing the full couplings \(\left\{ \pi^s \right\}_{s=0}^{T-1}\) suffers from the curse of dimensionality [4], [5]. To alleviate this issue, we develop convex relaxation methods that retain only low-order local information of \(\pi^s\), such as a sparse collection of marginals or cluster moments. These relaxations, developed in 2, provide computable lower bounds for problem 10 and recover low-order statistics of the intermediate distributions.

It remains to recover Markov kernels in the original problem 1 . If the full couplings were available, the kernels could be obtained by disintegration as in 7 . After relaxation, however, only partial local information is available. In the dynamic optimal transport case, the Benamou–Brenier dynamics is represented by a velocity field; 3 shows how this velocity can be recovered from dual variables of 10 , yielding an associated Markov kernel through the induced flow. These dual variables can, in turn, be approximated by solutions of the dual convex relaxation. For more general constrained Markov processes, we fit a parametrized family of kernels to the recovered local statistics. 4 develops this fitting procedure and illustrates it with a Markov process between Ising distributions based on Glauber dynamics.

1.3 Our contributions↩︎

We summarize the main contributions.

• We recast the Markov process optimization problem 1 into a sequential-coupling formulation 10 , and propose a convex relaxation framework for solving it in high dimensions without incurring exponential computational costs. The framework can handle both discrete and continuous state spaces, and the relaxed solution provides low-order statistics of the intermediate distributions.

• We demonstrate that problem 1 encompasses dynamic optimal transport as a special case, and we develop a procedure to recover the underlying Benamou–Brenier dynamics from the relaxed solution of 10 .

• We show that our procedure extends beyond dynamic optimal transport to handle other types of dynamics. As a concrete example, we illustrate our method with a Markov process between Ising models governed by Glauber dynamics.

1.4 Related work↩︎

Dynamic optimal transport is a central special case of our framework. The Kantorovich problem and the Benamou–Brenier formula give the classical static and dynamic viewpoints on Wasserstein transport [1]–[3]. After the momentum substitution, the Benamou–Brenier problem is convex. Eulerian augmented-Lagrangian and proximal-splitting methods exploit this structure effectively in low dimensions, but their space–time grids suffer from the curse of dimensionality [3], [6]. High-dimensional alternatives use Lagrangian or sample-based representations and neural parameterizations of velocity fields or flow maps, including neural discretizations of velocity fields and flow maps [7]–[9], TrajectoryNet-type methods [10], and flow-matching variants based on minibatch optimal transport or convex parameterizations [11], [12]. These methods avoid full spatial grids and scale well empirically, but their training problems are generally nonconvex. In finite-sample implementations, endpoint matching and kinetic costs are estimated from samples, minibatch couplings, or penalties; sampling, approximation, optimization, and ODE-discretization errors are therefore intertwined, intermediate laws are usually assessed empirically, and certified lower bounds for the original unregularized action are generally unavailable. Our approach instead optimizes local marginals or cluster moments through a convex program, providing a computable lower bound, low-order intermediate statistics, and approximate dual velocity information. Before spatial relaxation, the sequential-coupling formulation also has the exact Benamou–Brenier value on any prescribed time grid.

Beyond the unconstrained Benamou–Brenier setting, a related literature restricts the admissible distributional dynamics themselves. Constrained dynamic optimal transport on parameterized families and dynamical transport for nonlinear control-affine systems replace the free velocity field by model-specific evolution laws [13], [14]. After time discretization and, when possible, elimination of auxiliary controls, such models may induce one-step kernel or coupling constraints of the type considered here. Their numerical treatment, however, is generally tailored to the prescribed dynamics; in particular, the control-affine method in [14] still relies on spatial discretization.

Schrödinger bridge problems also connect prescribed endpoint distributions through Markovian dynamics, but minimize path-space relative entropy with respect to a prior process [15], [16]. Recent diffusion and discrete variants provide scalable stochastic interpolations [17]–[19]. By contrast, our formulation allows general additive one-step costs and hard constraints on admissible transition kernels, without requiring a reference process or entropic regularization.

At the relaxation level, moment–SOS methods provide convex hierarchies for generalized moment and polynomial optimization [20], [21], and have been applied to static optimal transport [22]. The associated dense moment matrices, however, grow rapidly with dimension. Sparse SOS hierarchies address this issue by exploiting correlative sparsity [23], [24]; more directly, the cluster marginal and moment relaxations in [25] provide the closest methodological antecedent to our approach for high-dimensional static transport. Convex moment relaxations have also been developed for Fokker–Planck dynamics [26]. Our framework combines these ideas by applying local marginal or moment representations to each adjacent coupling and linking them through multiple constraints. The resulting program is therefore a convex relaxation of the full multistage problem, rather than a collection of independent static relaxations.

1.5 Organization↩︎

The rest of the paper is organized as follows. 2 develops the marginal and moment relaxations for the sequential-coupling problem 10 . 3 identifies dynamic optimal transport as a special case of problem 10 and shows how the underlying Benamou–Brenier dynamics can be recovered from the relaxed solution of 10 . 4 develops a kernel-fitting procedure for more general Markov processes and illustrates it with a Markov process between Ising models governed by Glauber dynamics. 5 presents numerical experiments and 6 concludes the paper.

2 Convex relaxation↩︎

In this section, we develop a convex relaxation framework for problem 10 . The construction is motivated by the marginal and cluster moment relaxations for static optimal transport in [25], whose basic idea is to exploit the locality and sparse interaction structure of the problem by representing a coupling \(\pi\) only through its local information on selected coordinate clusters and cluster pairs. In the sequential-coupling formulation, this local representation is applied to every coupling \(\pi^s\), and the constraints in 10 are imposed through the retained local variables. In 2.1, we introduce the basic notation and cluster structure used for this representation. 2.2 presents the marginal relaxation, which is natural for finite or discretized state spaces. 2.3 presents the moment relaxation, which is well suited to continuous state spaces, where direct marginal discretization is computationally costly.

2.1 Notation and cluster graph↩︎

Let \([d]=\{1,\ldots,d\}\) be the coordinate set, and assume that the state space has the product form \(\Omega=\Omega_1\times\cdots\times\Omega_d\). For \(x\in\Omega\), write \(x=(x_1,\ldots,x_d)\) where \(x_i\in\Omega_i\). A cluster is a subset of coordinates. We choose a collection of clusters \[\begin{align} \mathcal{C}= \left\{ A_1, \ldots, A_K \right\}, \quad K = |\mathcal{C}|, \; A_i \subset \left[ d \right], \; 1 \le i \le K \end{align}\] to specify which local groups of coordinates will be represented explicitly. The purpose of the clusters is to retain local marginals while avoiding a full representation on all \(d\) coordinates. If each cluster is small, then distributions on one cluster or on a pair of clusters remain low-dimensional while still capturing local dependence structure. Unless otherwise stated, we assume that \(\mathcal{C}\) is a partition of \(\left[ d \right]\): \[\begin{align} A_a\cap A_b=\varnothing \quad (a\ne b), \quad \bigcup_{a=1}^K A_a=[d]. \end{align}\] For any \(A\subset[d]\), define \[\begin{align} x_A=(x_i)_{i\in A}, \quad \Omega_A=\prod_{i\in A}\Omega_i. \end{align}\] For an adjacent pair \(z = (x,y)\in\Omega\times\Omega\), define \[\begin{align} z_A=(x_A,y_A), \quad z_A \in Z_A=\Omega_A\times\Omega_A. \end{align}\] Thus \(z_A\) contains the coordinates in cluster \(A\) at two consecutive time layers. Throughout the paper, lowercase subscripts such as \(i,j\), together with numerical subscripts such as \(1,2\), refer to individual coordinates, whereas capital subscripts such as \(A,B\) refer to clusters.

The pairwise information between clusters retained by the relaxation is encoded by a cluster graph \[\begin{align} \mathcal{G}=(\mathcal{C},\mathcal{E}), \end{align}\] where the vertex set \(\mathcal{C}\) is the family of all clusters and \(\mathcal{E}\) is the edge set: \[\begin{align} \mathcal{E}\subset \left[ \mathcal{C} \right]_2, \quad \left[ \mathcal{C} \right]_2 = \left\{ AB:A,B\in\mathcal{C}, A \neq B \right\}. \end{align}\] For brevity, we write \(AB\) for the unordered pair \(\left\{ A,B \right\}\). An edge \(AB\in\mathcal{E}\) means that we retain the joint information between clusters \(A\) and \(B\). Sparse choices of \(\mathcal{G}\) lead to smaller convex programs.

For a probability measure \(\eta\) on \(\Omega\), its one-cluster and cluster-pair marginals are denoted by \[\begin{align} \eta_A=({\rm P}_A)_{\sharp} \, \eta, \quad \eta_{AB}=({\rm P}_{AB})_{\sharp} \, \eta, \qquad A \in \mathcal{C}, \quad AB \in [\mathcal{C}]_2, \end{align}\] where \({\rm P}_A:\Omega\to\Omega_A\) and \({\rm P}_{AB}:\Omega\to\Omega_A\times\Omega_B\) are the canonical projections. Similarly, for a coupling \(\pi^s\) on \(\Omega\times\Omega\), we denote its one-cluster and cluster-pair marginals by \[\begin{align} \pi^s_A = ({\rm P}_A )_\sharp \, \pi^s, \quad \pi^s_{AB} = ({\rm P}_{AB})_\sharp \, \pi^s, \qquad A \in \mathcal{C}, \quad AB \in [\mathcal{C}]_2. \end{align}\] Here, by a slight abuse of notation, the same symbols \(\mathrm{P}_A\) and \(\mathrm{P}_{AB}\) are used for the corresponding projections on \(\Omega \times \Omega\).

For an array of functions \(\left\{ f_a(z) \right\}_{a \in \mathcal{A}}\) over the two-layer variable \(z = (x,y) \in \Omega \times \Omega\), we introduce selection operators \({\rm R}_x\) and \({\rm R}_y\), which retain only the entries that depend on the left and right time layers, respectively. \[\begin{align} \label{def:RxRy} \begin{aligned} {\rm R}_x(\left\{ f_a(z) \right\}_{a \in \mathcal{A}}) \mathrel{\vcenter{:}}= \left\{ f_a \mid \exists \, \widetilde{f} \text{ such that } f_a(x,y) = \widetilde{f}(x), a \in \mathcal{A} \right\}, \\ {\rm R}_y(\left\{ f_a(z) \right\}_{a \in \mathcal{A}}) \mathrel{\vcenter{:}}= \left\{ f_a \mid \exists\, \widetilde{f} \text{ such that } f_a(x,y) = \widetilde{f}(y), a \in \mathcal{A} \right\}. \end{aligned} \end{align}\tag{15}\] For example, applying \({\rm R}_x\) and \({\rm R}_y\) to a matrix of basis functions yields \[\begin{align} {\rm R}_x\left( \begin{bmatrix} 1 & y_2 \\ x_1^2 & x_1y_2 \end{bmatrix}\right) = \begin{bmatrix} 1 \\ x_1^2 \end{bmatrix}, \quad {\rm R}_y\left( \begin{bmatrix} 1 & y_2 \\ x_1^2 & x_1y_2 \end{bmatrix}\right) = \begin{bmatrix} 1 \\ y_2 \end{bmatrix}. \end{align}\]

Finally, given a finite signed measure \(\eta\) on a measurable space \(W\) and an integrable scalar-, vector-, or matrix-valued test function \(\Xi:W\to\mathbb{C}^{m\times n}\), we write its moment under \(\eta\) as \[\begin{align} \eta(\Xi)=\int_W \Xi(w)\mathop{}\!\mathrm{d}\eta(w), \end{align}\] with the integral understood entrywise. For example, when \(\Xi_{AB}\) is a function on the set \(Z_A \times Z_B = \left( \Omega_A \times \Omega_A \right) \times \left( \Omega_B \times \Omega_B \right)\), its moment under the cluster-pair marginal \(\pi_{AB}\) is \[\begin{align} \pi_{AB}^s(\Xi_{AB})=\int_{Z_A\times Z_B}\Xi_{AB}(z_A,z_B)\mathop{}\!\mathrm{d}\pi_{AB}^s(z_A,z_B). \end{align}\]

2.2 Marginal relaxation↩︎

In this subsection, we present the marginal relaxation of problem 10 . Given a specified cluster decomposition \(\mathcal{C}\) of the coordinate set \([d]\), which is intended to group strongly correlated coordinates together, the full coupling \(\pi^s(x,y)\) on each time interval \(s = 0, \ldots, T-1\) is replaced by its local marginals on clusters and cluster pairs \[\begin{align} \label{notation:reduced32var32Mar} \{\pi_A^s\}_{A\in\mathcal{C}}, \quad \{\pi_{AB}^s\}_{AB\in \left[ \mathcal{C} \right]_2}. \end{align}\tag{16}\] Here \(\pi_A^s\) is a probability measure on \(Z_A\), and \(\pi_{AB}^s\) is a probability measure on \(Z_{A}\times Z_B\). In practice, the size of each cluster is usually much smaller than \(d\), so these local objects can be represented directly on the chosen finite state space or spatial grid. We also use the cluster graph \(\mathcal{G}= (\mathcal{C}, \mathcal{E})\) introduced above. Its edge set specifies which cluster-pair marginals are retained in the objective and subsequent constraints.

Assume that the one-step cost either admits the local decomposition \[\begin{align} \label{eq:relaxed32cost32Mar} c^s(z) = \sum_{A\in\mathcal{C}} c_{A}^s(z_A) + \sum_{AB\in\mathcal{E}} c_{AB}^s(z_A,z_B), \end{align}\tag{17}\] or is approximated by the right-hand side. If the decomposition is exact and the local marginals are induced by a full coupling \(\pi^s\), then \[\begin{align} \int_{\Omega\times\Omega}c^s(z)\mathop{}\!\mathrm{d}\pi^s(z) = \sum_{A\in\mathcal{C}}\int_{Z_A}c_A^s(z_A)\mathop{}\!\mathrm{d}\pi_A^s(z_A)+\sum_{AB\in\mathcal{E}}\int_{Z_A\times Z_B}c_{AB}^s(z_A,z_B)\mathop{}\!\mathrm{d}\pi_{AB}^s(z_A,z_B). \end{align}\] For example, the quadratic cost \(c^s(x,y)= \|x-y\|^2\) decomposes over a partition \(\mathcal{C}\) as \(\|x-y\|^2 = \sum_{A\in\mathcal{C}} \|x_A-y_A\|^2\). Thus, in this case, only one-cluster cost terms are needed. More general local costs may include the cluster-pair terms \(c_{AB}^s\).

We next impose constraints on the reduced variables.

2.2.0.1 1. Marginal consistency constraints.

For every pair \(AB \in [\mathcal{C}]_2\) and each \(s = 0, \ldots, T-1\), the one-cluster and cluster-pair marginals must be compatible in the sense that \[\label{constraint-local32consis32Mar} ({\rm P}_{A} )_\sharp \, \pi^s_{AB} = \pi^s_A, \quad ({\rm P}_B)_\sharp \, \pi^s_{AB} = \pi^s_B, \quad \forall AB \in [\mathcal{C}]_2.\tag{18}\] These constraints are necessary if the family \(\left\{ \pi_A^s, \pi_{AB}^s \right\}\) is to be interpreted as the collection of one-cluster and cluster-pair marginals of a common global coupling \(\pi^s\).

2.2.0.2 2. Initial and terminal constraints.

The initial and terminal constraints in the full sequential-coupling problem, namely 11 , require the left marginal of the first coupling to be \(\mu\) and the right marginal of the last coupling to be \(\nu\). In the marginal relaxation, these constraints are relaxed to hold only on the retained cluster-pair marginals: \[\begin{align} \label{constraint-endpoint32Mar} \left( {\rm P}_{\rm L} \right)_\sharp\pi_{AB}^0=\mu_{AB}, \quad \left( {\rm P}_{\rm R} \right)_\sharp\pi_{AB}^{T-1}=\nu_{AB}, \quad AB\in\mathcal{E}. \end{align}\tag{19}\] If every cluster is incident to at least one edge in \(\mathcal{E}\), then the corresponding initial and terminal constraints for one-cluster marginals follow from 18 and 19 . If isolated clusters are allowed, their initial and terminal constraints are added explicitly.

2.2.0.3 3. Mass conservation constraints.

The right marginal of \(\pi^s\) can be interpreted as the distribution of mass arriving at time \(t_{s+1}\), whereas the left marginal of \(\pi^{s+1}\) represents the distribution of mass leaving at time \(t_{s+1}\). The mass conservation constraints 12 ensure that these two marginals coincide. In the marginal relaxation, this condition is imposed after projection onto the retained cluster pairs: \[\begin{align} \label{constraint-time32consis32Mar} \left( {\rm P}_{\rm R} \right)_\sharp\pi_{AB}^s=\left( {\rm P}_{\rm L} \right)_\sharp\pi_{AB}^{s+1}, \quad AB\in\mathcal{E}, \quad s = 0,\ldots T-2. \end{align}\tag{20}\] Thus 20 is the projected local form of mass conservation across adjacent time intervals.

2.2.0.4 4. Markov kernel constraints.

The admissible-kernel constraint in the Markov formulation, \(K^s \in \mathcal{K}^s\), is expressed in the sequential-coupling formulation as \(\pi^s \in \mathcal{D}^s\) where \(\mathcal{D}^s\) denotes the admissible class of couplings at step \(s\). At the cluster level, we replace \(\mathcal{D}^s\) by a convex local admissible set \(\mathcal{D}_{\mathrm{Mar}}^s\), and we impose \[\begin{align} \label{constraint-control32Mar} \{\pi_A^s,\pi_{AB}^s\}_{A\in\mathcal{C},\,AB\in[\mathcal{C}]_2}\in\mathcal{D}_{\mathrm{Mar}}^s,\qquad s=0,\ldots,T-1. \end{align}\tag{21}\]

The set \(\mathcal{D}_{\mathrm{Mar}}^s\) is a convex outer approximation that contains the local marginals induced by any admissible coupling \(\pi^s \in \mathcal{D}^s\). Thus, this constraint remains a relaxation of the original Markov kernel constraint.

2.2.0.5 5. Local positivity constraints.

The nonnegativity constraint 14 for the full coupling implies nonnegativity of all its local marginals. We impose this at the cluster-pair level: \[\begin{align} \label{constraint-local32positivity32Mar} \pi_{AB}^s \ge 0, \quad AB \in [\mathcal{C}]_2,\quad s = 0,\cdots T-1. \end{align}\tag{22}\] The one-cluster nonnegativity \(\pi_A^s \ge 0\) follows from 18 and 22 , provided each cluster appears in at least one retained pair. If isolated one-cluster variables are allowed, their nonnegativity is imposed explicitly.

2.2.0.6 6. Global positivity constraints.

The preceding local constraints do not guarantee that the family \(\left\{ \pi_A^s, \pi_{AB}^s \right\}_{A,AB}\) is induced by a genuine global coupling \(\pi^s\). We therefore impose the following necessary global positivity condition. If the local variables are induced by a genuine coupling \(\pi^s\), then for every family of square-integrable test functions \(\left\{ f_A \in L^2(\pi_A^s) \right\}_{A \in \mathcal{C}}\), we have \[\begin{align} & \sum_{A\in\mathcal{C}}\int_{Z_A} f_A(z_A)^2 \mathop{}\!\mathrm{d}\pi_A^s(z_A) + 2\sum_{AB\in[\mathcal{C}]_2}\int_{Z_A\times Z_B} f_A(z_A)f_B(z_B) \mathop{}\!\mathrm{d}\pi_{AB}^s(z_A,z_B) \nonumber \\ & = \pi^s\left( \left( \sum_{A \in \mathcal{C}}f_A(z_A) \right)^2 \right) \ge 0. \end{align}\] In the finite-state or gridded setting, this condition is equivalent to the PSD constraint \[\begin{align} \label{constraint-PSD32Mar} \begin{bmatrix} \operatorname{Diag}(\pi_{A_1}^s) & \pi_{A_1A_2}^s & \cdots & \pi_{A_1A_K}^s\\ (\pi_{A_1A_2}^{s})^\top & \operatorname{Diag}(\pi_{A_2}^s) & \cdots & \pi_{A_2A_K}^s\\ \vdots & \vdots & \ddots & \vdots\\ (\pi_{A_1A_K}^s)^{\top} & (\pi_{A_2A_K}^s)^{\top} & \cdots & \operatorname{Diag}(\pi_{A_K}^s) \end{bmatrix} \succeq 0, \quad s=0,\ldots, T-1, \end{align}\tag{23}\] where \(K = |\mathcal{C}|\) is the number of clusters and \(\operatorname{Diag}(\pi_{A_k}^s)\) denotes the diagonal matrix with diagonal entries given by \(\pi_{A_k}^s\). We write this compactly as \((\pi_A^s,\pi_{AB}^s)_{A \in \mathcal{C}, AB \in [\mathcal{C}]_2}\succeq 0\).

Combining the preceding constraints gives the full marginal relaxation: \[\begin{align} \min_{\left\{ \pi^s_A \right\}_{\mathcal{C}}, \left\{ \pi^s_{AB} \right\}_{[\mathcal{C}]_2}} \; & \sum_{s=0}^{T-1} \left[ \sum_{A\in\mathcal{C}}\int_{Z_A}c_A^s\mathop{}\!\mathrm{d}\pi_A^s+\sum_{AB\in\mathcal{E}}\int_{Z_A\times Z_B}c_{AB}^s\mathop{}\!\mathrm{d}\pi_{AB}^s \right] \tag{24} \\ \text{s.t.}\quad\quad & \left\{ \pi_A^s, \pi_B^s, \pi_{AB}^s \right\}_{AB \in [\mathcal{C}]_2} \text{ is consistent \eqref{constraint-local32consis32Mar}}, \tag{25}\\ & \text{initial and terminal constraints \eqref{constraint-endpoint32Mar} for } \left\{ \pi_{AB}^i \right\}_{\mathcal{E}}, \; i=0 \text{ or } T-1, \tag{26}\\ & \text{mass conservation \eqref{constraint-time32consis32Mar} for } \left\{ \pi_{AB}^s \right\}_{\mathcal{E}}, \tag{27} \\ & \text{Markov kernel constraints \eqref{constraint-control32Mar} for } \left\{ \pi_A^s, \pi_{AB}^s \right\}_{\mathcal{C}, [\mathcal{C}]_2}, \tag{28} \\ & \text{local positivity \eqref{constraint-local32positivity32Mar} for } \left\{ \pi_{AB}^s \right\}_{[\mathcal{C}]_2}, \tag{29}\\ & \text{global positivity \eqref{constraint-PSD32Mar}: } (\pi_{A}^s,\pi_{AB}^s)_{\mathcal{C}, \, [\mathcal{C}]_2}\succeq 0. \tag{30} \end{align}\] If \(\mathcal{D}^s_{\mathrm{Mar}}\) is described by affine equalities and inequalities, then 24 is a doubly nonnegative block semidefinite program (SDP): it combines pairwise entrywise nonnegativity with a PSD constraint for each time interval.

Following the approach in [25], we also consider two simpler relaxations of 24 . The first simplification drops the global positivity constraints 30 and retains pair variables only on the edge set \(\mathcal{E}\). Accordingly, the marginal consistency, Markov kernel, and local positivity constraints are imposed only for \(AB \in \mathcal{E}\). The non-edge variables \(\pi_{AB}^s\) with \(AB\notin\mathcal{E}\) are then omitted. This gives \[\begin{align} \min_{\left\{ \pi^s_A \right\}_{\mathcal{C}}, \left\{ \pi^s_{AB} \right\}_{\mathcal{E}}} \; & \sum_{s=0}^{T-1}\left[ \sum_{A\in\mathcal{C}}\int_{Z_A}c_A^s\mathop{}\!\mathrm{d}\pi_A^s+\sum_{AB\in\mathcal{E}}\int_{Z_A\times Z_B}c_{AB}^s\mathop{}\!\mathrm{d}\pi_{AB}^s \right] \tag{31} ^1}\\ \text{s.t.}\quad\quad & \left\{ \pi_A^s, \pi_B^s, \pi_{AB}^s \right\}_{AB \in \mathcal{E}} \text{ is consistent \eqref{constraint-local32consis32Mar}}, \tag{32}^1-a}\\ & \text{initial and terminal constraints \eqref{constraint-endpoint32Mar} for } \left\{ \pi_{AB}^i \right\}_{\mathcal{E}}, \; i = 0 \text{ or } T-1, \tag{33}^1-b}\\ & \text{mass conservation \eqref{constraint-time32consis32Mar} for } \left\{ \pi_{AB}^s \right\}_{\mathcal{E}}, \tag{34} ^1-c}\\ & \text{projected Markov kernel constraints: } \left\{ \pi_A^s, \pi_{AB}^s \right\}_{\mathcal{C}, \, \mathcal{E}} \in \mathcal{D}^s_{\mathrm{Mar}, \mathcal{E}}, \tag{35} ^1-d}\\ & \text{local positivity \eqref{constraint-local32positivity32Mar} for } \left\{ \pi_{AB}^s \right\}_{\mathcal{E}}, \tag{36} ^1-e} \end{align}\] Here \(\mathcal{D}^s_{\mathrm{Mar}, \mathcal{E}}\) in 35 denotes the projection of the constraint set \(\mathcal{D}^s_{\mathrm{Mar}}\) onto the retained variables. When \(\Omega\) is finite and these constraints are affine, 31 reduces to a linear program.

The second simplification drops the local positivity constraints 29 while retaining the global positivity constraints. The resulting relaxation is \[\begin{align} \min_{\left\{ \pi^s_A \right\}_{\mathcal{C}}, \left\{ \pi^s_{AB} \right\}_{[\mathcal{C}]_2}} \; & \sum_{s=0}^{T-1} \left[ \sum_{A\in\mathcal{C}}\int_{Z_A}c_A^s\mathop{}\!\mathrm{d}\pi_A^s+\sum_{AB\in\mathcal{E}}\int_{Z_A\times Z_B}c_{AB}^s\mathop{}\!\mathrm{d}\pi_{AB}^s \right] \tag{37} ^2}\\ \text{s.t.}\quad\quad & \left\{ \pi_A^s, \pi_B^s, \pi_{AB}^s \right\}_{AB \in [\mathcal{C}]_2} \text{ is consistent \eqref{constraint-local32consis32Mar}}, \tag{38}^2-a}\\ & \text{initial and terminal constraints \eqref{constraint-endpoint32Mar} for } \left\{ \pi_{AB}^i \right\}_{ \mathcal{E}}, i = 0 \text{ or } T-1, \tag{39}^2-b}\\ & \text{mass conservation \eqref{constraint-time32consis32Mar} for } \left\{ \pi_{AB}^s \right\}_{ \mathcal{E}}, \tag{40} ^2-c}\\ & \text{Markov kernel constraints \eqref{constraint-control32Mar} for } \left\{ \pi_A^s, \pi_{AB}^s \right\}_{\mathcal{C},\, [\mathcal{C}]_2}, \tag{41} ^2-d}\\ & \text{global positivity \eqref{constraint-PSD32Mar}: } (\pi_{A}^s,\pi_{AB}^s)_{\mathcal{C}, \,[\mathcal{C}]_2}\succeq 0. \tag{42} ^2-e} \end{align}\] The variables \(\pi_{AB}^s\) are now allowed to be signed measures. The global positivity constraints still imply \(\pi_A^s\ge 0\), because each diagonal block in the finite-dimensional PSD matrix is \(\operatorname{Diag}(\pi_A^s)\succeq 0\). However, they do not imply pointwise nonnegativity of the pair variables \(\pi_{AB}^s\).

2.3 Cluster moment relaxation↩︎

We now introduce the cluster moment relaxation. As in the marginal relaxation, we fix a cluster decomposition \(\mathcal{C}\) and the corresponding cluster graph \(\mathcal{G}= (\mathcal{C}, \mathcal{E})\). The moment relaxation replaces the local marginals used in the marginal relaxation by finitely many moments against prescribed cluster bases. We first define the cluster bases.

Fix a degree \(r \in \mathbb{N}^+\). For each coordinate \(i \in [d]\), choose a finite one-dimensional basis of order at most \(r\), denoted by \(\left\{ \phi_{i,j}: \Omega_i \rightarrow \mathbb{C} \right\}_{j=0}^r\), where \(\phi_{i,0} \equiv 1\). Typical choices include monomial bases, orthogonal polynomial bases, Fourier bases, or other bases adapted to the one-dimensional domain \(\Omega_i\). Recall from 2.1 that, for a cluster \(A \in \mathcal{C}\), \(z_A = (x_A, y_A) \in \Omega_A \times \Omega_A\), where \(x_A\) and \(y_A\) denote the variables of the left time layer \(t_s\) and right time layer \(t_{s+1}\), respectively. For a multi-index \[\begin{align} \alpha = (\alpha^x, \alpha^y), \quad \alpha^x = (\alpha_i^x)_{i \in A}, \quad \alpha^y = (\alpha_i^y)_{i \in A} \end{align}\] we define its total degree by \[|\alpha| = \sum_{i \in A} \alpha_i^x + \sum_{i\in A} \alpha_i^y.\] Let \(\mathcal{I}_A^r = \left\{ \alpha = (\alpha^x, \alpha^y) \mid |\alpha| \le r \right\}\). For each \(\alpha \in \mathcal{I}_A^r\), define the cluster basis function on \(Z_A\) by \[\begin{align} \Phi_{A,\alpha}(z_A) = \prod_{i\in A}\phi_{i, \alpha_i^x}(x_i) \phi_{i, \alpha_i^y}(y_i). \end{align}\] We collect these functions into a column vector \(\Phi_A = \left( \Phi_{A, \alpha} \right)_{\alpha \in \mathcal{I}_A^r}\), which defines the two-layer cluster basis associated with cluster \(A\). This basis consists of functions supported only on the variables \(z_A = (x_A, y_A)\) and having total degree at most \(r\). For example, if \(A = \left\{ 1 \right\}\), \(r = 2\) and \(\phi_{1,j}(u_1) = u_1^j\), then, up to ordering, \(\Phi_A = \left( 1, x_1, x_1^2, y_1, y_1^2, x_1y_1 \right)^\top\).

We now define the moment variables. For each time interval \(s = 0, \ldots, T-1\) and each cluster \(A \in \mathcal{C}\), introduce the one-cluster moment matrix: \[\begin{align} \label{eq:MA} M_A^s \mathrel{\vcenter{:}}= \pi_A^s(\Phi_A \Phi_A^*) = \int_{Z_A}\Phi_A(z_A)\Phi_A(z_A)^* \mathop{}\!\mathrm{d}\pi_A^s(z_A), \end{align}\tag{43}\] where the integral is understood entrywise and \(\left( \cdot \right)^*\) denotes the Hermitian transpose. Similarly, for each cluster pair \(AB \in [\mathcal{C}]_2\), the cross-moment matrix is defined as \[\begin{align} \label{eq:MAB} M_{AB}^s \mathrel{\vcenter{:}}= \pi_{AB}^s(\Phi_A \Phi_B^*) = \int_{Z_A \times Z_B}\Phi_A(z_A)\Phi_B(z_B)^* \mathop{}\!\mathrm{d}\pi_{AB}^s(z_A,z_B). \end{align}\tag{44}\] The decision variables in the moment relaxation are the collections \[\begin{align} \label{notation:reduced32var32Mom} \{M_A^s\}_{A\in\mathcal{C}}, \quad \{M_{AB}^s\}_{AB\in \left[ \mathcal{C} \right]_2}, \quad s = 0, \ldots, T-1. \end{align}\tag{45}\] These variables at each time step \(s\) can be assembled into a single Hermitian matrix \[\begin{align} \label{def:Ms} M^s:= \begin{bmatrix} M_{A_1}^s & M_{A_1A_2}^s & \cdots & M_{A_1A_K}^s \\ \left( M_{A_1A_2}^s \right)^* & M_{A_2}^s & \cdots & M_{A_2A_K}^s\\ \vdots & \vdots & \ddots & \vdots\\ \left( M_{A_1A_K}^s \right)^* & \left( M_{A_2A_K}^s \right)^* & \cdots & M_{A_K}^s \end{bmatrix}. \end{align}\tag{46}\] The overall decision variable is the block-diagonal matrix \(M = \mathrm{Diag} \left( M^0, \cdots, M^{T-1} \right)\).

We now describe the objective and constraints in terms of these moment variables. Recall from the previous subsection that we represent the one-step cost \(c^s\) via a local decomposition 17 : \[\begin{align} c^s(z) = \sum_{A\in\mathcal{C}} c_{A}^s(z_A) + \sum_{AB\in\mathcal{E}} c_{AB}^s(z_A,z_B). \end{align}\] Assume that each local cost component lies in the span of the corresponding product basis, or is replaced by its projection onto that span. Let \(C_A^s\) and \(C_{AB}^s\) denote the representations of \(c_A^s\) and \(c_{AB}^s\) in the bases formed by the elements of \(\Phi_A \Phi_A^*\) and \(\Phi_A \Phi_B^*\), respectively. Specifically, \[\begin{align} c_A^s=\langle C_A^s,\Phi_A\Phi_A^*\rangle, \quad c_{AB}^s =\langle C_{AB}^s,\Phi_A\Phi_B^*\rangle. \end{align}\] In the exact case, the local objective at step \(s\) can be written as \[\begin{align} \label{eq:cost32decomp32moment} \int_{\Omega \times \Omega}c^s(x,y) \mathop{}\!\mathrm{d}\pi^s(x,y) = \sum_{A\in\mathcal{C}}\langle C_{A}^s,M_A^s\rangle+\sum_{AB\in\mathcal{E}}\langle C_{AB}^s,M_{AB}^s\rangle. \end{align}\tag{47}\] We next impose constraints on the reduced moment variables. By the definition of selection operators in 15 , \({\rm R}_x(\Phi_A\Phi_B^*)\) and \({\rm R}_y(\Phi_A\Phi_B^*)\) extract the submatrices of entries in \(\Phi_A\Phi_B^*\) that depend exclusively on the left-layer variables \((x_A, x_B)\) and right-layer variables \((y_A,y_B)\), respectively. For each time interval \(s = 0, \ldots, T-1\), cluster \(A \in \mathcal{C}\) and cluster pair \(AB \in [\mathcal{C}]_2\), we define the following shorthand notation: \[\begin{align} (M_A^s)_{\rm L} \mathrel{\vcenter{:}}= \pi_A^s({\rm R}_x(\Phi_A\Phi_A^*)), & \quad (M_A^s)_{\rm R} \mathrel{\vcenter{:}}= \pi_A^s({\rm R}_y(\Phi_A\Phi_A^*)), \\[1ex] (M_{AB}^s)_{\rm L} \mathrel{\vcenter{:}}= \pi_{AB}^s({\rm R}_x(\Phi_A\Phi_B^*)), & \quad (M_{AB}^s)_{\rm R} \mathrel{\vcenter{:}}= \pi_{AB}^s({\rm R}_y(\Phi_A\Phi_B^*)). \end{align}\]

2.3.0.1 1. Moment consistency constraints.

For each cluster \(A\in\mathcal{C}\), the one-cluster moment matrix \(M_A^s=\pi_A^s(\Phi_A\Phi_A^*)\) contains entries indexed by products of basis functions in the same variables \(z_A=(x_A,y_A)\). Different pairs of basis functions may therefore generate the same product function. For example, if \(A=\{1\}\) and the one-dimensional basis is the monomial basis \(\phi_{1,j}(u_1)=u_1^j\), then the products \(x_1\cdot x_1^3\) and \(x_1^2\cdot x_1^2\) both represent the same function \(x_1^4\). The corresponding entries of \(M_A^s\) must therefore be equal. These identities give linear equality constraints on \(M_A^s\), for every \(A\in\mathcal{C}\) and \(s=0,\ldots,T-1\). We denote them compactly by saying that \[\begin{align} \label{constraint-consis32Mom} M^s \text{ is consistent, } \quad s = 0, \ldots, T-1. \end{align}\tag{48}\] Note that no analogous constraint is imposed on the cross-moment matrix \(M_{AB}^s=\pi_{AB}^s(\Phi_A\Phi_B^*)\) since the functions in \(\Phi_A\) and \(\Phi_B\) are supported on disjoint coordinate sets.

2.3.0.2 2. Initial and terminal constraints.

The initial and terminal constraints 19 of the marginal relaxation are further relaxed at the moment level to hold only for one-cluster moments and retained cluster-pair moments. Explicitly, we have \[\label{constraint-endpoint32Mom} \begin{align} (M_A^0)_{\rm L}=\mu_A\left({\rm R}_x(\Phi_A\Phi_A^*)\right), & \quad (M_A^{T-1})_{\rm R}=\nu_A\left({\rm R}_y(\Phi_A\Phi_A^*)\right), \quad A\in\mathcal{C}\\[1ex] (M_{AB}^0)_{\rm L}=\mu_{AB}\left({\rm R}_x(\Phi_A\Phi_B^*)\right), & \quad (M_{AB}^{T-1})_{\rm R}=\nu_{AB}\left({\rm R}_y(\Phi_A\Phi_B^*)\right), \quad AB\in\mathcal{E}. \end{align}\tag{49}\]

2.3.0.3 3. Mass conservation constraints.

The mass conservation constraints 20 are relaxed to the moment level for all clusters \(A \in \mathcal{C}\) and retained cluster pairs \(AB \in \mathcal{E}\): \[\label{constraint-time32consis32Mom} \begin{align} (M_A^s)_{\rm R} & =(M_A^{s+1})_{\rm L}, \quad A\in\mathcal{C}, \quad s=0,\ldots,T-2 \\ (M_{AB}^s)_{\rm R} & =(M_{AB}^{s+1})_{\rm L}, \quad AB\in\mathcal{E}, \quad s=0,\ldots,T-2. \end{align}\tag{50}\]

2.3.0.4 4. Markov kernel constraints.

The Markov kernel constraints 21 are imposed at the moment level through a convex set \(\mathcal{D}^s_{\mathrm{Mom}}\): \[\begin{align} \label{constraint-control32Mom} \{M_A^s,M_{AB}^s\}_{A \in \mathcal{C}, \, AB \in [\mathcal{C}]_2}\in\mathcal{D}^s_{\mathrm{Mom}}, \quad s=0,\ldots,T-1. \end{align}\tag{51}\] This is the moment-level relaxation of the admissible transition constraint \(\pi^s\in\mathcal{D}^s\). For instance, if coordinate \(i\) is required to remain frozen at step \(s\), then the support constraint \(x_i-y_i=0\) can be imposed by linear equations of the form \[\begin{align} \int (x_i-y_i)q(z)\mathop{}\!\mathrm{d}\pi^s(z)=0 \end{align}\] for all retained test functions \(q\) such that \((x_i-y_i)q\) lies in the span of the retained product basis.

2.3.0.5 5. Local positivity constraints.

The local positivity constraints 22 require that, for each pair \(AB \in [\mathcal{C}]_2\), the triple \((M_A^s, M_B^s, M_{AB}^s)\) is compatible with a nonnegative measure on \(Z_A \times Z_B\). We impose them through a convex relaxation of local moment representability. For each \(AB\in[\mathcal{C}]_2\), let \(\mathcal{K}_{AB}^{\mathrm{loc}}\) denote the chosen convex set of admissible local moment triples. We impose \[\begin{align} \label{constraint-local32positivity32Mom} (M_A^s,M_B^s,M_{AB}^s)\in\mathcal{K}_{AB}^{\mathrm{loc}}, \quad AB\in[\mathcal{C}]_2, \quad s=0,\ldots,T-1. \end{align}\tag{52}\]

2.3.0.6 6. Global positivity constraints.

The global positivity constraints 23 are relaxed by restricting the square-integrable test functions to the span of the cluster bases. Consider the function set \(\mathcal{F}= \left\{ f= \sum_{A \in \mathcal{C}} v_A^*\Phi_A(z_A) \mid v_A \in \mathbb{C}^{|\Phi_A|} \right\}\). If the moments are induced by a genuine coupling \(\pi^s\), then, for every \(f \in \mathcal{F}\), we have \[\begin{align} \sum_{A\in\mathcal{C}}v_A^* M_A^s v_A+2\mathrm{Re}\sum_{AB\in[\mathcal{C}]^2}v_A^* M_{AB}^s v_B = \pi^s\left( |f|^2 \right) \ge 0. \end{align}\] This is equivalent to requiring the moment matrix defined in 46 to be positive semidefinite: \[\begin{align} \label{constraint-PSD32Mom} M^s\succeq 0, \quad s=0,\ldots,T-1. \end{align}\tag{53}\]

Combining these constraints gives the cluster moment relaxation: \[\begin{align} \min_{\left\{ M^s \right\}_{s=0}^{T-1}} \; & \sum_{s=0}^{T-1} \left[ \sum_{A\in\mathcal{C}}\langle C_{A}^s,M_A^s\rangle+\sum_{AB\in\mathcal{E}}\langle C_{AB}^s,M_{AB}^s\rangle \right] \tag{54} }\\ \text{s.t.}\quad& M^s \text{ is consistent \eqref{constraint-consis32Mom}}, \tag{55}-a} \\ & \text{initial and terminal constraints \eqref{constraint-endpoint32Mom} for } \left\{ M_A^i, M_{AB}^i \right\}_{\mathcal{C}, \, \mathcal{E}}, \; i = 0 \text{ or } T-1, \tag{56}-b} \\ & \text{mass conservation \eqref{constraint-time32consis32Mom} for } \left\{ M_A^s, M_{AB}^s \right\}_{\mathcal{C}, \, \mathcal{E}}, \tag{57} -c}\\ & \text{Markov kernel constraints \eqref{constraint-control32Mom} for } \left\{ M_A^s, M_{AB}^s \right\}_{\mathcal{C}, \, [\mathcal{C}]_2}, \tag{58} -d} \\ & \text{local positivity \eqref{constraint-local32positivity32Mom}: }(M_A^s,M_B^s,M_{AB}^s)\in\mathcal{K}_{AB}^{\mathrm{loc}}, \quad AB\in[\mathcal{C}]_2, \tag{59}-e} \\ & \text{global positivity \eqref{constraint-PSD32Mom}: } M^s\succeq 0. \tag{60} -f} \end{align}\] This is a semidefinite relaxation of the sequential-coupling problem. As in the marginal case, we obtain two further relaxations by dropping either the global positivity constraints \(M^s \succeq 0\) or the local positivity constraints \((M_A^s,M_B^s,M_{AB}^s)\in \mathcal{K}_{AB}^{\mathrm{loc}}\).

\[\begin{align} \min_{\substack{\left\{ M^s_A \right\}_\mathcal{C}\\ \left\{ M_{AB}^s \right\}_\mathcal{E}}} \; & \sum_{s=0}^{T-1} \left[ \sum_{A\in\mathcal{C}}\langle C_{A}^s,M_A^s\rangle+\sum_{AB\in\mathcal{E}}\langle C_{AB}^s,M_{AB}^s\rangle \right] \tag{61} ^1}\\ \text{s.t.}\quad& M^s \text{ is consistent \eqref{constraint-consis32Mom}}, \tag{62} ^1-a}\\ & \text{initial and terminal constraints \eqref{constraint-endpoint32Mom} for } \left\{ M_A^i, M_{AB}^i \right\}_{\mathcal{C}, \, \mathcal{E}} ,\; i = 0 \text{ or } T-1, \tag{63}^1-b}\\ & \text{mass conservation \eqref{constraint-time32consis32Mom} for } \left\{ M_A^s, M_{AB}^s \right\}_{\mathcal{C}, \, \mathcal{E}}, \tag{64} ^1-c}\\ & \text{projected Markov kernel constraints: }\{M_A^s,M_{AB}^s\}_{\mathcal{C}, \, \mathcal{E}} \in \mathcal{D}^s_{\mathrm{Mom},\mathcal{E}}, \tag{65} ^1-d}\\ & \text{local positivity \eqref{constraint-local32positivity32Mom}: }(M_A^s,M_B^s,M_{AB}^s)\in\mathcal{K}_{AB}^{\mathrm{loc}}, \quad AB\in \mathcal{E}. \tag{66} ^1-e} \end{align}\] Problem 61 is obtained from 54 by dropping the global positivity constraints and discarding the cross-moment variables \(M^s_{AB}\) for \(AB \notin \mathcal{E}\). The set \(\mathcal{D}_{\mathrm{Mom},\mathcal{E}}^s\) denotes the projection of the constraint set \(\mathcal{D}_{\mathrm{Mom}}^s\) onto the retained variables. Dropping the local positivity constraints 59 from 54 yields the following relaxation: \[\begin{align} \min_{\left\{ M^s \right\}_{s=0}^{T-1}} \; & \sum_{s=0}^{T-1} \left[ \sum_{A\in\mathcal{C}}\langle C_{A}^s,M_A^s\rangle+\sum_{AB\in\mathcal{E}}\langle C_{AB}^s,M_{AB}^s\rangle \right] \tag{67} ^2}\\ \text{s.t.}\quad& M^s \text{ is consistent \eqref{constraint-consis32Mom}}, \tag{68}^2-a} \\ & \text{initial and terminal constraints \eqref{constraint-endpoint32Mom} for } \left\{ M_A^i, M_{AB}^i \right\}_{\mathcal{C}, \, \mathcal{E}}, \; i=0 \text{ or } T-1, \tag{69}^2-b} \\ & \text{mass conservation \eqref{constraint-time32consis32Mom} for } \left\{ M_A^s, M_{AB}^s \right\}_{\mathcal{C}, \, \mathcal{E}}, \tag{70} ^2-c}\\ & \text{Markov kernel constraints \eqref{constraint-control32Mom} for } \left\{ M_A^s, M_{AB}^s \right\}_{\mathcal{C}, \, [\mathcal{C}]_2}, \tag{71} ^2-d} \\ & \text{global positivity \eqref{constraint-PSD32Mom}: } M^s\succeq 0. \tag{72} ^2-e} \end{align}\]

3 Recovering Benamou–Brenier dynamics from problem 10↩︎

In this section, we show that the sequential-coupling problem 10 contains the Benamou–Brenier problem 6 as a special case. We work on a convex set \(\Omega\subset\mathbb{R}^d\), and let \(\mu, \nu \in \mathcal{P}_p(\Omega)\), where \(\mathcal{P}_p(\Omega)\) denotes the set of probability measures on \(\Omega\) with finite \(p\)-th moment. Fix a time grid \(0 = t_0 < t_1 < \cdots < t_T =1\). For \(1\le p<\infty\), choose the one-step cost \[\begin{align} \label{eq:dot-cost323461} c^s(x,y)=\frac{\|x-y\|^p}{(\Delta t_s)^{p-1}}, \quad \Delta t_s=t_{s+1}-t_s. \end{align}\tag{73}\] We impose no Markov kernel constraints. In the notation of problem 10 , this means \(\mathcal{D}^s = \mathcal{P}(\Omega \times \Omega)\) in 13 . With these choices, problem 10 becomes: \[\begin{align} \inf_{\left\{ \pi^s \right\}_{s=0}^{T-1}} \; & \sum_{s=0}^{T-1} \int_{\Omega\times\Omega} \frac{\|x-y\|^p}{(\Delta t_s)^{p-1}} \mathop{}\!\mathrm{d}\pi^s(x,y) \tag{74} \\ \text{s.t.}\quad& \left( {\rm P}_{\rm L} \right)_\sharp\pi^0 = \mu, \quad \left( {\rm P}_{\rm R} \right)_\sharp\pi^{T-1} = \nu, \tag{75} \\ & \left( {\rm P}_{\rm R} \right)_\sharp\pi^s = \left( {\rm P}_{\rm L} \right)_\sharp\pi^{s+1}, \quad s=0,\ldots,T-2, \tag{76} \\ & \pi^s \ge 0, \quad s=0,\ldots,T-1. \tag{77} \end{align}\]

This section shows how solving 74 recovers an optimal solution of the Benamou–Brenier problem 6 . In 3.1, we prove that these two problems have the same optimal value and that a primal solution of 74 yields the Benamou–Brenier measure curve at the prescribed grid times. In 3.2, under a standard connected-support condition, we show that a dual solution of 74 recovers the Benamou–Brenier velocity field at the prescribed grid times. Finally, in 3.3, we explain how to extract an approximate velocity field from the moment relaxation for 74 .

3.1 Recovering the measure curve in 6 from primal solution of 74↩︎

In this subsection, 1 shows that an optimal solution of 74 recovers the values of a Benamou–Brenier measure curve at the prescribed grid times. We also show that 74 has the same optimal value as the Benamou–Brenier problem 6 .

To prove the result, we first introduce the intermediate distributions induced by a feasible sequence of couplings. For any feasible sequence \(\left\{ \pi^s \right\}_{s=0}^{T-1}\), define \[\begin{align} \rho^0 \mathrel{\vcenter{:}}= \left( {\rm P}_{\rm L} \right)_\sharp\, \pi^0, \quad \rho^{s+1} \mathrel{\vcenter{:}}= \left( {\rm P}_{\rm R} \right)_\sharp\, \pi^s, \quad s = 0, \ldots, T-1. \end{align}\] By the initial and terminal constraints 75 and the mass-conservation constraints 76 , we have \[\begin{align} \rho^0 = \mu, \quad \rho^T = \nu, \quad \pi^s \in \Pi(\rho^s, \rho^{s+1}), \quad s = 0, \ldots, T-1, \end{align}\] where \(\Pi(\rho^s, \rho^{s+1})\) denotes the set of all couplings between \(\rho^s\) and \(\rho^{s+1}\). Conversely, any sequence \(\left\{ \rho^s \right\}_{s=0}^T\) with \(\rho^0 = \mu\) and \(\rho^T = \nu\), together with couplings \(\pi^s\in \Pi(\rho^s, \rho^{s+1})\), defines a feasible solution of 74 . Therefore, minimizing 74 over the sequence of couplings is equivalent to first choosing the intermediate marginals at the prescribed grid times \(\left\{ \rho^s \right\}_{s=0}^T\), and then choosing an optimal coupling between each pair of consecutive marginals.

For fixed \(\rho^{s}\) and \(\rho^{s+1}\), minimizing over \(\pi^s \in \Pi(\rho^s, \rho^{s+1})\) gives \[\begin{align} \inf_{\pi^s\in\Pi(\rho^s,\rho^{s+1})}\int_{\Omega\times\Omega}\frac{\|x-y\|^p}{(\Delta t_s)^{p-1}}\mathop{}\!\mathrm{d}\pi^s(x,y)=\frac{1}{(\Delta t_s)^{p-1}} \mathcal{W}_p^p(\rho^s,\rho^{s+1}). \end{align}\] Thus 74 is equivalent to the reduced form \[\begin{align} \label{problem:dot39} \inf_{\left\{ \rho^s \in \mathcal{P}_p(\Omega) \right\}_{s=0}^T}\; \sum_{s=0}^{T-1}\frac{1}{(\Delta t_s)^{p-1}}\mathcal{W}_p^p(\rho^s,\rho^{s+1}),\qquad \rho^0=\mu,\qquad \rho^T=\nu. \end{align}\tag{78}\] We call an optimal measure curve \(\{q(\, \cdot \,, t)\}_{t\in[0,1]}\) in 6 a Benamou–Brenier geodesic between \(\mu\) and \(\nu\). The following proposition explains how to recover the Benamou–Brenier geodesic from the solution of 74 .

Thus, with \(\mathcal{D}^s = \mathcal{P}(\Omega \times \Omega)\) and the cost defined in 73 , problem 10 recovers a Benamou–Brenier geodesic at the prescribed grid times and introduces no time-discretization error in the optimal value. Applying the marginal or moment relaxations from 2 to this formulation then yields low-order statistics of the intermediate distributions and a lower bound for \(\mathcal{W}_p^p\).

3.2 Recovering the velocity in 6 from dual solution of 74↩︎

In this subsection, we derive the dual of problem 74 and show that, under a standard connected-support condition, its dual potentials encode the Benamou–Brenier velocity field at the grid times. We restrict the presentation to the case \(p=2\). Analogous formulas can be derived for \(1 \le p < \infty\). The main statement is as follows.

Thus, in this special case under the above assumptions, we can recover the Benamou–Brenier velocity field at the grid times from the gradients of the dual potentials. When this velocity field generates a flow map, the associated Markov kernel in problem 1 is obtained from the induced transition along the flow.

Before proving 2, we recall the standard Benamou–Brenier duality result.

We now prove 2.

3.3 Approximating the velocity in 6 from moment relaxation for 74↩︎

3.2 shows how to recover the Benamou–Brenier velocity field from the dual potentials of problem 74 in the quadratic case \(p=2\). In this subsection, we explain how to obtain an approximation of the velocity field from a dual solution of the cluster moment relaxation. We focus on the simplified relaxation 67 , where the local positivity constraints are dropped while the global positivity constraints are retained.

The relaxed dual problem may have more than one optimizer. Throughout this subsection, the phrase the dual solution refers to any fixed optimal solution of the relaxed dual problem. The construction below is applied to this optimizer, and its common additive constant has no effect on the recovered velocity field.

We derive the dual of the moment relaxation 67 for problem 74 directly from a measure-level formulation. Let \[\begin{align} \mathcal{Z}=\Omega\times\Omega,\quad z=(x,y)\in \mathcal{Z},\quad s=0,\ldots,T-1. \end{align}\] Instead of optimizing over moment matrices, the relaxation can be viewed as optimizing over signed measures \(\eta^s\in\mathcal{M}(\mathcal{Z})\). Let \(\Phi(z)=(\Phi_A(z_A))_{A\in\mathcal{C}}\) denote the vector of retained cluster basis functions. Define the associated finite-dimensional SOS cone \[\begin{align} \Sigma^s:=\{\langle S^s,\Phi(z)\Phi(z)^*\rangle:S^s\succeq 0\}. \end{align}\] Equivalently, \(\Sigma^s\) consists of functions of the form \(\sigma(z)=\sum_{\ell}|q_{\ell}(z)|^2\), where \(q_{\ell}\in\operatorname{span}\{\Phi\}\). Let \[\begin{align} \mathcal{W}_\mathcal{G}\mathrel{\vcenter{:}}= \operatorname{span}\left\{ {\rm R}_x \left( \left\{ \Phi_A\Phi_A^*: A \in \mathcal{C} \right\} \cup \left\{ \Phi_A \Phi_B^*: AB \in \mathcal{E} \right\} \right) \right\}, \end{align}\] where \({\rm R}_x\) selects the entries that depend only on the \(x\)-copy of \(z = (x,y)\). We can write the measure-level relaxation as \[\begin{align} \inf_{\{\eta^s\}_{s=0}^{T-1}}\quad & \sum_{s=0}^{T-1}\langle c^s,\eta^s\rangle \tag{79} \\ \text{s.t.}\quad& \langle f \circ {\rm P}_{\rm L},\eta^0\rangle = \langle f,\mu\rangle, \quad \forall f\in\mathcal{W}_\mathcal{G}, \\ & \langle f \circ {\rm P}_{\rm R},\eta^{T-1}\rangle=\langle f,\nu\rangle,\quad \forall f\in\mathcal{W}_\mathcal{G}, \\ & \langle f\circ {\rm P}_{\rm R},\eta^s\rangle = \langle f \circ {\rm P}_{\rm L},\eta^{s+1}\rangle, \quad \forall f \in\mathcal{W}_\mathcal{G},\quad s=0,\ldots,T-2, \\ & \langle\sigma,\eta^s\rangle\ge 0,\quad \forall\sigma\in\Sigma^s,\quad s=0,\ldots,T-1. \tag{80}\end{align}\] This formulation is equivalent to the moment-matrix relaxation 67 applied to problem 74 . Indeed, given signed measures \(\eta^s\), define \[M^s=\int_{\mathcal{Z}}\Phi(z)\Phi(z)^*\mathop{}\!\mathrm{d}\eta^s(z).\] Then, for any \(\sigma(z)=\langle S^s,\Phi(z)\Phi(z)^*\rangle\in\Sigma^s\), we have \[\begin{align} \label{eq:measure-matrix} \langle\sigma,\eta^s\rangle=\left\langle S^s,\int_{\mathcal{Z}}\Phi(z)\Phi(z)^*\mathop{}\!\mathrm{d}\eta^s(z)\right\rangle=\langle S^s,M^s\rangle. \end{align}\tag{81}\] Therefore \(\left\langle \sigma,\eta^s \right\rangle \ge 0\) for any \(\sigma\in\Sigma^s\) is equivalent to \(M^s\succeq 0\). Conversely, let \(M^s\) be a feasible moment matrix for relaxation 67 for problem 74 . The consistency constraints imply that \(M^s\) defines a well-defined linear functional \(L^s\) on the finite-dimensional space \(\mathcal{V}^s=\operatorname{span}\{\text{entries of } \Phi(z)\Phi(z)^*\}\) after identifying repeated product functions. Since \(\mathcal{V}^s\) is finite dimensional, any such linear functional can be represented on \(\mathcal{V}^s\) by a signed atomic measure: there exist points \(z_\ell \in \mathcal{Z}\) and weights \(w_\ell \in \mathbb{R}\) such that \(\eta^s=\sum_{\ell}w_{\ell}\delta_{z_{\ell}}\) and \[\begin{align} \int f\mathop{}\!\mathrm{d}\eta^s=L^s(f),\quad \forall f\in\mathcal{V}^s, \quad s = 0, \ldots, T-1. \end{align}\] In particular, \[M^s=\int \Phi(z)\Phi(z)^*\mathop{}\!\mathrm{d}\eta^s(z).\] Thus the measure formulation 79 and the matrix formulation 67 applied to problem 74 have the same feasible retained moments and the same optimal value.

We now derive the dual of 79 . Let \[\widehat{\lambda}^0, \widehat{\lambda}^T \in \mathcal{W}_\mathcal{G}; \qquad \widehat{\lambda}^s \in \mathcal{W}_\mathcal{G}, \quad s = 1, \ldots, T-1\] be the dual potentials associated with the endpoint and mass-conservation constraints, respectively, and let \[\tau^s \in \Sigma^s, \quad s = 0, \ldots, T-1\] be the multipliers associated with the relaxed positivity constraints 80 . The Lagrangian of the problem 79 is \[\begin{align} \mathcal{L} = &\sum_{s=0}^{T-1}\langle c^s,\eta^s\rangle + \langle\widehat{\lambda}^0\circ {\rm P}_{\rm L},\eta^0\rangle-\langle\widehat{\lambda}^0,\mu\rangle +\langle\widehat{\lambda}^T , \nu \rangle-\langle\widehat{\lambda}^T\circ {\rm P}_{\rm R} , \eta^{T-1}\rangle \\ & +\sum_{s=0}^{T-2}\left[\langle\widehat{\lambda}^{s+1}\circ {\rm P}_{\rm L},\eta^{s+1}\rangle-\langle\widehat{\lambda}^{s+1}\circ {\rm P}_{\rm R},\eta^s\rangle\right] -\sum_{s=0}^{T-1}\langle\tau^s,\eta^s\rangle \\ = & \langle\widehat{\lambda}^T,\nu\rangle-\langle\widehat{\lambda}^0,\mu\rangle+\sum_{s=0}^{T-1}\left\langle c^s+\widehat{\lambda}^s\circ {\rm P}_{\rm L}-\widehat{\lambda}^{s+1}\circ {\rm P}_{\rm R}-\tau^s,\eta^s\right\rangle. \end{align}\] Since \(\eta^s\) is a signed measure, the infimum of \(\mathcal{L}\) over \(\eta^s\) is finite only if \[\begin{align} c^s(x,y)+\widehat{\lambda}^s(x)-\widehat{\lambda}^{s+1}(y) = \tau^s(x,y), \quad s=0,\ldots,T-1. \end{align}\] Substituting the quadratic Benamou–Brenier cost \(c^s(x,y) = \|x-y\|^2 / \Delta t_s\), the dual problem is \[\begin{align} \sup_{\widehat{\lambda}^s\in \mathcal{W_\mathcal{G}}, \, \tau^s \in\Sigma^s} \quad & \int_{\Omega}\widehat{\lambda}^T(x)\mathop{}\!\mathrm{d}\nu(x)-\int_{\Omega}\widehat{\lambda}^0(x)\mathop{}\!\mathrm{d}\mu(x) \\ \text{s.t.}\quad& \frac{\|x-y\|^2}{\Delta t_s}-\widehat{\lambda}^{s+1}(y)+\widehat{\lambda}^s(x)=\tau^s(x,y),\quad s=0,\ldots,T-1. \label{conssos} \end{align}\tag{82}\] Using the Gram representation of \(\Sigma^s\), the constraint 82 can be written as \[\begin{align} \label{eq:SOS32identity} \frac{\|x-y\|^2}{\Delta t_s} - \widehat{\lambda}^{s+1}(y)+\widehat{\lambda}^s(x)=\langle S^s,\Phi(z)\Phi(z)^*\rangle,\qquad S^s\succeq 0,\qquad s=0,\ldots,T-1. \end{align}\tag{83}\] This shows that the dual of the moment relaxation 79 for problem 74 is an SOS relaxation of the exact dual problem ?? . Indeed, the pointwise nonnegativity constraint in the exact dual problem ?? is replaced by the stronger SOS certificate in the retained cone \(\Sigma^s\). Consequently, the dual potentials \(\widehat{\lambda}^s\) should be interpreted as approximations of the exact dual potentials \(\lambda^s\) in the prescribed finite-dimensional basis \(\mathcal{W_\mathcal{G}}\).

We now approximate the velocity field. By 2, in the quadratic Benamou–Brenier setting, the optimal velocity field is given at the grid times by \(v(x,t_s)=\frac{1}{2}\nabla\lambda^s(x)\). Therefore, from the moment dual potentials, we define the approximate velocity field associated with the fixed dual solution by \[\begin{align} \label{eq:approximate32v} \widehat{v}(x,t_s)=\frac{1}{2}\nabla\widehat{\lambda}^s(x), \quad s=0,\ldots,T-1. \end{align}\tag{84}\] The SOS identity 83 gives a concrete formula for \(\widehat{\lambda}^s\). Equating the terms in 83 that depend only on \(x\), we obtain, up to an additive constant, \[\begin{align} \label{eq:approximate32lambda} \widehat{\lambda}^s(x)=-\frac{\|x\|^2}{\Delta t_s}+R_x(\langle S^s,\Phi(z)\Phi(z)^*\rangle). \end{align}\tag{85}\] Combining 84 and 85 then yields \[\begin{align} \widehat{v}(x,t_s)=-\frac{x}{\Delta t_s} + \frac{1}{2}\nabla R_x(\langle S^s,\Phi(z)\Phi(z)^*\rangle). \end{align}\]

In summary, the primal solution of the moment relaxation recovers low-order statistics of the intermediate distributions, while its dual potentials approximate the Benamou–Brenier potentials through an SOS relaxation. Differentiating these potentials gives the approximate Benamou–Brenier velocity field, which can then be used to generate the associated approximate flow and hence the corresponding Markov transition kernel in problem 1 .

4 Recovering general dynamics from problem 10↩︎

The previous section shows that, when \(\Omega\) is convex, the one-step cost is given by the kinetic cost 73 , and no Markov kernel constraints are imposed (i.e., \(\mathcal{D}^s=\mathcal{P}(\Omega\times\Omega)\) in 13 ), problem 10 exactly recovers the Benamou–Brenier geodesic. In more general settings, such as finite or nonconvex state spaces, the optimal solution of 10 need not correspond to a classical Benamou–Brenier geodesic or to prescribed dynamics. Moreover, the convex relaxations in 2 retain only local statistics of the optimal sequential couplings, rather than the full one-step transition kernels. In this section, we introduce a parametric family of transition kernels \(\{K_{\theta_s}^s\}_{\theta_s}\) and fit the parameters \(\theta_s\) so that the one-step evolution induced by \(K_{\theta_s}^s\) matches the recovered local statistics at the next time step.

We take the marginal relaxation as a representative case. After solving the relaxation, we extract the local marginals of the intermediate distributions \[\begin{align} \{\rho_A^s\}_{A\in\mathcal{C}}, \quad \{\rho_{AB}^s\}_{AB\in\mathcal{E}}, \quad s=0,\ldots,T. \end{align}\] In the kernel-fitting step, we compare the next-time cluster-pair marginals predicted by the parameterized kernel with the recovered ones. The fitting loss may be evaluated on all retained cluster pairs, or only on those local marginals that are affected by the update and are informative for identifying the local transition law. Let \(\mathcal{E}^s_{\mathrm{fit}}\subseteq \mathcal{E}\) denote this selected set of retained cluster pairs. For a cluster pair \(AB \in \mathcal{E}_\mathrm{fit}^s\), define the next-time cluster-pair marginal predicted by the parameterized kernel as \[\begin{align} \label{eq:next-time32marginal} \widehat{\rho}_{AB,\theta_s}^{\, s+1}=({\rm P}_{AB})_\sharp \, \left( (K_{\theta_s}^s)^*\rho^s \right), \end{align}\tag{86}\] where \({\rm P}_{AB}\) denotes the projection onto the coordinates in \(AB\). Formula 86 is written at the full-law level. It can be evaluated locally whenever the \(AB\)-marginal of the one-step kernel depends on the current configuration only through a coordinate neighborhood \(N(AB)\). In a finite state space, this gives \[\begin{align} \widehat{\rho}_{AB,\theta_s}^{\, s+1}(y_{AB})=\sum_{x_{N(AB)}} \rho_{N(AB)}^{\, s} \left( x_{N(AB)} \right) K_{\theta_s,AB}^s \left( y_{AB}\mid x_{N(AB)} \right), \end{align}\] where \(K^s_{\theta_s,AB}\) denotes the output law on cluster \(AB\) induced by \(K^s_{\theta_s}\), and \(\rho^s_{N(AB)}\) is the current marginal on the input coordinates needed to evaluate this local transition law. It remains to obtain the local input marginal \(\rho^s_{N(AB)}\) from the relaxed solution. If \(N(AB)\) is contained in a cluster or a retained cluster pair, then \(\rho^s_{N(AB)}\) is obtained by projection from the corresponding relaxed marginal. Otherwise, we reconstruct it from the available one-cluster and cluster-pair marginals. A standard choice is the Bethe-type reconstruction. For example, on a path graph with single-coordinate clusters, suppose that an update at an interior site \(i\) requires the current marginal on \(\{i-1,i,i+1\}\), but only the pair marginals \(\rho^s_{i-1,i}\) and \(\rho^s_{i,i+1}\) are retained. Then the Bethe-type reconstruction is \[\begin{align} \rho_{i-1,i,i+1}^{s}(a,b,c)=\frac{\rho_{i-1,i}^s(a,b)\rho_{i,i+1}^s(b,c)}{\rho_i^s(b)}, \quad a\in \Omega_{i-1}, \quad b \in \Omega_i, \quad c \in \Omega_{i+1}, \quad \rho_i^s(b) > 0. \end{align}\]

We then fit \(\theta_s\) at each step by minimizing the squared loss \(L^s(\theta_s)\): \[\begin{align} L^s(\theta_s) \mathrel{\vcenter{:}}= \sum_{AB\in\mathcal{E}^s_{\mathrm{fit}}}\left\|\widehat{\rho}_{AB,\theta_s}^{\,s+1}-\rho_{AB}^{s+1}\right\|_2^2+\lambda R(\theta_s), \end{align}\] where \(R\) is a regularization term and \(\lambda\ge 0\) is the penalty strength. Thus the fitting step projects the recovered local statistics, in a least-squares sense, onto the chosen parametric Markov-kernel family.

In what follows, we illustrate this procedure for a constrained Markov process between two Ising models, where only one spin is allowed to change at each step. 4.1 formulates the Ising laws and imposes the single-spin update constraint in the sequential-coupling problem. 4.2 then uses the local marginals of the intermediate distributions obtained from the marginal relaxation to fit a time-inhomogeneous Glauber dynamics, thereby converting the recovered statistics into an explicit Markov kernel.

4.1 Ising models and the single-spin update process↩︎

We consider the finite-state spin configuration space \(\Omega = \left\{ -1,1 \right\}^d\), where each coordinate \(x_i \in \left\{ -1,1 \right\}\) represents the spin at site \(i\). An Ising model is specified by an interaction graph \[\begin{align} \mathcal{G}_{\mathrm{Ising}}=([d],\mathcal{E}_{\mathrm{Ising}}), \quad \mathcal{E}_{\mathrm{Ising}} \subset [d]_2 \mathrel{\vcenter{:}}= \left\{ \{i,j\} : 1 \le i < j \le d \right\}, \end{align}\] pairwise interaction strengths \(\{J_{ij}\}_{ij\in \mathcal{E}_{\mathrm{Ising}}}\), external fields \(h=(h_1,\ldots,h_d)\), and inverse temperature \(\beta>0\). For brevity, we write \(ij\) for the unordered edge \(\left\{ i,j \right\}\). The corresponding probability mass function is \[\begin{align} \label{eq:Ising32model} p_{\theta}(x)=\frac{1}{Z_{\theta}}\exp\left[\beta\left(\sum_{ij\in \mathcal{E}_{\mathrm{Ising}}}J_{ij}x_i x_j+\sum_{i=1}^d h_i x_i\right)\right], \qquad x\in\{-1,1\}^d, \end{align}\tag{87}\] where \(\theta = (\mathcal{G}_{\mathrm{Ising}},J,h,\beta)\) and \(Z_{\theta}\) is the normalizing constant. The source and target laws are taken to be \(\mu=p_{\theta_{\mu}}\) and \(\nu=p_{\theta_{\nu}}\). The initial and terminal laws may have different interaction strengths, external fields, inverse temperatures, or interaction graphs. The graph \(\mathcal{G}_{\mathrm{Ising}}\) records the physical interaction structure of the Ising law and should be distinguished from the cluster graph \(\mathcal{G}= (\mathcal{C}, \mathcal{E})\) used in the relaxation scheme.

We impose a single-spin update constraint on the Markov process. Let \(0 = t_0 < t_1< \cdots < t_T = 1\) be a time grid and let \[\begin{align} i_s\in[d], \qquad s=0,\ldots,T-1, \end{align}\] be the active spin at time step \(s\), during which only the spin \(i_s\) is allowed to change while all other spins remain fixed. Thus the admissible Markov kernels satisfy \[\begin{align} \label{eq:kernel32constraint32Ising} K^s(y\mid x)=0 \quad \text{whenever } (x,y)\notin \Gamma^s, \quad \Gamma^s \mathrel{\vcenter{:}}= \{(x,y)\in\Omega\times\Omega:y_j=x_j \text{ for all } j\ne i_s\}. \end{align}\tag{88}\] In the sequential-coupling formulation, this becomes the linear support constraint \[\begin{align} \pi^s(x,y)=0 \quad \text{whenever } (x,y)\notin \Gamma^s. \end{align}\] The single-spin update process between the initial and terminal laws is therefore the following instance of the sequential-coupling problem 10 : \[\begin{align} \inf_{\{\pi^s\}_{s=0}^{T-1}}\; & \sum_{s=0}^{T-1}\sum_{x,y\in\Omega}c^s(x,y)\pi^s(x,y) \tag{89} \\ \text{s.t.}\quad& \left( {\rm P}_{\rm L} \right)_\sharp\pi^0=\mu,\quad \left( {\rm P}_{\rm R} \right)_\sharp\pi^{T-1}=\nu, \tag{90}\\ & \left( {\rm P}_{\rm R} \right)_\sharp\pi^s = \left( {\rm P}_{\rm L} \right)_\sharp\pi^{s+1},\qquad s=0,\ldots,T-2, \tag{91}\\ & \pi^s(x,y)=0 \quad \text{if } (x,y)\notin \Gamma^s, \tag{92}\\ & \pi^s(x,y)\ge 0. \tag{93} \end{align}\] A simple one-step cost is the quadratic spin-flip cost \(c^s(x,y)=\|x-y\|^2\). Under the single-spin update constraint, this cost is zero if the active spin does not change and equals \(4\) if the active spin flips. Other local costs can also be used.

We solve 89 using the marginal relaxation 31 from 2.2. Since \(\Omega\) is finite and the single-spin support constraints are affine, \(\mathrm{Mar}^1\) becomes a linear program (LP). Let \(\mathcal{C}=\{A_1,\ldots,A_K\}\) be a partition of \([d]\), and let \(\mathcal{G}=(\mathcal{C}, \mathcal{E})\) be the cluster graph used in the relaxation. Solving the resulting LP gives local marginals of the intermediate couplings. We then extract the local marginals of the intermediate distributions: \[\begin{align} \label{eq:cluster32marginals32Ising} \rho_A^s, \quad \rho_{AB}^s, \quad A \in \mathcal{C}, \quad AB \in \mathcal{E}, \quad s = 0, \ldots, T. \end{align}\tag{94}\]

4.2 Fitting a Glauber dynamics↩︎

The marginal relaxation recovers local cluster marginals of the intermediate distributions as in 94 , but it does not directly specify a Markov kernel. We now fit such a kernel from the recovered marginals using the method introduced at the beginning of 4. Because the admissible process has a single-spin update structure, we naturally choose a Glauber kernel. Specifically, over each time interval \([t_s, t_{s+1}]\), the Markov transition from \(\rho^s\) to \(\rho^{s+1}\) is parameterized as a single Glauber update with parameters \(\theta_s\).

Glauber dynamics defines a classical Markov chain for sampling from Ising models. In our parameterization, the inverse temperature is absorbed into the effective external fields and local interactions, so a Glauber kernel is specified by a graph \(\mathcal{G}_{\mathrm{glb}}\), external fields \(h\), and interaction strengths \(J\). The graph \(\mathcal{G}_{\mathrm{glb}} = ([d], \mathcal{E}_{\mathrm{glb}})\) denotes the dependency graph for the fitted Glauber dynamics and is fixed across all time steps. This graph need not coincide with either the initial or the terminal Ising interaction graph, and it is distinct from the cluster graph \(\mathcal{G}= (\mathcal{C}, \mathcal{E})\) used in the marginal relaxation. Let \(N_s\) be the neighborhood of the active spin \(i_s\) at time \(s\) induced by \(\mathcal{G}_{\mathrm{glb}}\): \[\begin{align} N_s \mathrel{\vcenter{:}}= \left\{ j \in [d]: \{i_s,j\} \in \mathcal{E}_\mathrm{glb} \right\}. \end{align}\] At each update step of the Glauber dynamics, the active spin is resampled from its local conditional distribution while all other spins remain fixed: \[\begin{align} \label{eq:Glauber32conditional} q_{\theta_s}^s(x_{i_s} = \sigma\mid x_{N_s})=\frac{\exp(\sigma H^s(x_{N_s}))}{2\cosh(H^s(x_{N_s}))}, \qquad \sigma\in\{-1,1\} \end{align}\tag{95}\] where the local field associated with the neighboring configuration \(x_{N_s}\) is defined as \[\begin{align} \label{def:local32field32H94s} H^s(x_{N_s})=h_{i_s}^s+\sum_{j\in N_s}J_{i_sj}^s x_j. \end{align}\tag{96}\] Consequently, the one-step transition kernel is given by \[\begin{align} \label{Markov32kernel32Ising} K_{\theta_s}^s(y\mid x)=\mathbf{1}_{\{y_j=x_j,\;j\ne i_s\}}q_{\theta_s}^s(y_{i_s}\mid x_{N_s}), \quad x, \;y \in \Omega. \end{align}\tag{97}\] The parameters to be fitted at step \(s\) are \(\theta_s = (h_{i_s}^s,\{J_{i_sj}^s:j\in N_s\})\).

Next, we describe how this kernel predicts the cluster-pair marginals at the subsequent time step. Let \(AB\in\mathcal{E}\) be a retained cluster pair. If the active spin \(i_s \notin AB\), then all coordinates in \(AB\) remain fixed during step \(s\), yielding \[\begin{align} \widehat{\rho}_{AB}^{\,s+1}(y_{AB})=\rho_{AB}^s(y_{AB}), \quad i_s\notin AB. \end{align}\] Such a marginal does not depend on \(\theta_s\). If instead \(i_s\in AB\), then the update probability of \(Y_{i_s}\) depends on the neighboring spins in \(N_s\). The predicted next-time cluster-pair marginal is computed from the current marginal on \((AB\setminus\{i_s\})\cup N_s\): \[\begin{align} \widehat{\rho}_{AB, \theta_s}^{\,s+1}(y_{AB})=\sum_{x_{N_s\setminus AB}}\rho_{(AB\setminus\{i_s\})\cup N_s}^s (y_{AB\setminus\{i_s\}},x_{N_s\setminus AB}) \; q_{\theta_s}^s(y_{i_s}\mid x_{N_s}). \end{align}\] In the conditional law \(q_{\theta_s}^s(y_{i_s}\mid x_{N_s})\), the vector \(x_{N_s}\) is formed by taking \(x_j = y_j\) for \(j \in N_s \cap AB\), and summing over \(x_j\) for \(j \in N_s \setminus AB\). The marginal \(\rho_{(AB\setminus\{i_s\})\cup N_s}^{\,s}\) is extracted from the recovered local marginals: if it is directly available from a retained cluster or cluster-pair marginal, we obtain it by projection; otherwise, it is reconstructed from the available local marginals by the Bethe-type local reconstruction. In particular, if the graph for Glauber dynamics \(\mathcal{G}_\mathrm{glb}\) is chosen so that \(N_s \subset AB\), then the expression simplifies to \[\begin{align} \widehat{\rho}_{AB, \theta_s}^{\,s+1}(y_{AB})=\rho_{AB\setminus\{i_s\}}^s(y_{AB\setminus\{i_s\}}) \; q_{\theta_s}^s(y_{i_s}\mid y_{N_s}). \end{align}\] In summary, the predicted cluster-pair marginals are obtained by keeping the inactive coordinates fixed and resampling the active coordinate according to the fitted Glauber conditional law.

To estimate \(\theta_s\), we match the predicted cluster-pair marginals with the recovered marginals at the subsequent time step. Let \(\mathcal{E}^s_{\mathrm{fit}} = \{AB\in\mathcal{E}:i_s\in AB\}\) denote the set of retained cluster pairs used for fitting at step \(s\). We solve \[\begin{align} \min_{\theta_s}\sum_{AB\in\mathcal{E}^s_{\mathrm{fit}}}\left\|\widehat{\rho}_{AB, \theta_s}^{\,s+1}-\rho_{AB}^{s+1}\right\|_2^2+\lambda R(\theta_s), \end{align}\] where \(\lambda\ge 0\) and \(R(\theta_s)\) is a suitable regularization term.

The Glauber structure also admits a computationally efficient reduced form of this fitting procedure. Let \(B_s=N_s\cup\{i_s\}\). If the recovered marginals are exactly generated by the single-site Glauber update, the log-odds ratio exhibits a linear structure: \[\begin{align} \frac{\rho_{B_s}^{s+1}(x_{N_s}=a,\;x_{i_s} = 1)}{\rho_{B_s}^{s+1}(x_{N_s}=a,\; x_{i_s} = -1)} = \exp(2H^s(a)),\quad a \in \left\{ -1,1 \right\}^{|N_s|}. \end{align}\] Taking logarithms and using the definition of \(H^s\) in 96 yields \[\begin{align} \label{eq:Glauber32fitting32linear32reg} \frac{1}{2}\log\frac{\rho_{B_s}^{s+1}(x_{N_s}=a, \; x_{i_s} = 1)}{\rho_{B_s}^{s+1}(x_{N_s}=a, \; x_{i_s} = -1)} = h_{i_s}^s+\sum_{j\in N_s}J_{i_sj}^s \, a_j. \end{align}\tag{98}\] Consequently, estimating \((h^s_{i_s},\{J^s_{i_sj}:j\in N_s\})\) reduces to a linear regression problem over configurations \(a\in\{-1,1\}^{|N_s|}\). In practice, configurations with very small mass are either omitted or regularized by adding a small positive floor before taking the logarithm.

Repeating this procedure for \(s = 0,\ldots, T-1\) yields a sequence of time-inhomogeneous Glauber kernels \(\left\{ K_{\theta_s}^s \right\}_{s=0}^{T-1}\) as defined in 97 . Given \(X^0\sim\mu\), these kernels define a discrete-time Markov chain by \[\begin{align} X^{s+1}\sim K_{\theta_s}^s(\, \, \cdot \,\mid X^s), \quad s=0,\ldots,T-1. \end{align}\] The resulting process can be used to propagate samples in practice and approximately realize the recovered local intermediate statistics as marginals of an implementable Markov process.

5 Numerical experiments↩︎

In this section, we present numerical experiments illustrating the proposed relaxations and reconstruction methods. In 5.1, we first test our method on Benamou–Brenier dynamics between Gaussian distributions, where the exact geodesic and velocity field are known. This example allows us to evaluate the accuracy of the recovered intermediate statistics and the extracted velocity field, as well as the effect of the cluster graph. We also compare our method with a particle-based back-propagation method to assess their empirical accuracy. In 5.2, we consider a more challenging transport problem from a Gaussian distribution to a Ginzburg–Landau distribution. This experiment illustrates the reconstruction of intermediate low-dimensional marginals and the velocity-induced particle flow beyond the Gaussian setting. Finally, in 5.3, we study a constrained Markov process between Ising distributions, where the marginal relaxation is combined with the fitting procedure in 4 to fit Glauber dynamics.

5.1 Benamou–Brenier dynamics between Gaussian distributions↩︎

We first consider the Benamou–Brenier dynamics between two Gaussian distributions \(\mu = \mathcal{N}(m_0, \Sigma_0)\) and \(\nu = \mathcal{N}(m_1, \Sigma_1)\). In this setting, both the Benamou–Brenier geodesic and the velocity field are available in closed form, so the recovered intermediate moments and velocity coefficients can be compared directly with the ground truth.

The Benamou–Brenier geodesic is the Gaussian curve \(\mathcal{N}(m_t, \Sigma_t)\), where \[\begin{align} \label{eq:Gaussian32mean32cov} m_t=(1-t)m_0+tm_1, \quad \Sigma_t=B_t\Sigma_0B_t^\top, \end{align}\tag{99}\] with \[\begin{align} B_t=(1-t)I+tQ, \quad Q = \Sigma_0^{-1/2}\left(\Sigma_0^{1/2}\Sigma_1\Sigma_0^{1/2}\right)^{1/2}\Sigma_0^{-1/2}. \end{align}\] The optimal velocity field along the Benamou–Brenier geodesic is affine in \(x\). More precisely, \[\begin{align} \label{eq:Gaussian32velocity32field} v(x,t)=A_tx+b_t, \quad A_t=(Q-I)B_t^{-1}, \quad b_t=m_1-m_0-A_tm_t. \end{align}\tag{100}\] These formulas are used to compute the reference mean, covariance, and optimal velocity field at the discrete time grid.

We sample the initial and terminal means \(m_0, m_1 \in \mathbb{R}^d\) independently from the standard normal distribution \(\mathcal{N}(0, I_d)\). The covariance matrices are generated through sparse precision matrices. More precisely, we first generate two symmetric banded matrices \(P_0, P_1\) with bandwidth \(3\). For \(|i-j|\le 3\), \(i\ne j\), the off-diagonal entries are sampled randomly, and the diagonal entries are chosen to make the precision matrices strictly diagonally dominant: \[\begin{align} (P_{\ell})_{ii}=\delta+\sum_{j\ne i}|(P_{\ell})_{ij}|,\qquad \ell=0,1, \end{align}\] where \(\delta>0\) is a diagonal shift. We then set \(\Sigma_0 = P_0^{-1}\) and \(\Sigma_1 = P_1^{-1}\). This construction gives Gaussian distributions with sparse precision matrices and decaying correlations.

We use single-coordinate clusters \(\mathcal{C}= \left\{ \left\{ 1 \right\}, \ldots, \left\{ d \right\} \right\}\). Let \(G_{\mathrm{path}}\) be the one-dimensional path graph on \([d]\). For an integer radius \(r \ge 1\), define the graph power \(G^r_\mathrm{path}\) by connecting \(i\) and \(j\) whenever \(\operatorname{dist}_{G_\mathrm{path}}(i,j)\le r\). In the experiments below, the cluster graph \(\mathcal{G}\) is chosen as \(G_\mathrm{path}^r\). We use the simplified moment relaxation 67 to solve the problem. For each coordinate \(i\), we use the degree-one monomial basis \(\left\{ \phi_{i,0}(x_i) = 1, \phi_{i,1}(x_i) = x_i \right\}\). Thus each moment matrix \(M^s = \pi^s\left( \Phi \Phi^* \right)\) contains the first and second moments of two adjacent time layers, together with their cross-moments. The one-step cost \(c^s(x,y) = \|x-y\|^2/\Delta t\) can be represented exactly by the selected basis.

We set \(d=50\), the precision-matrix diagonal shift \(\delta=1\), the graph radius \(r=16\), and \(T=10\). We use an equispaced time grid on \([0,1]\) with \(\Delta t = 1/T\). The resulting problem is a multi-block SDP with one dense semidefinite block \(M^s\) per time interval, and is solved by MOSEK [29] with a tolerance of \(10^{-10}\). Throughout this subsection, vector and matrix errors are reported as relative Euclidean and relative Frobenius errors, respectively. For sparse cluster graphs, we report a zero-filled covariance error. Specifically, for a recovered covariance matrix \(\widehat{\Sigma}^s\) at time \(t_s\), entries not retained by the cluster graph are set to zero: \[\begin{align} \widehat{\Sigma}^s_{G_\mathrm{path}^r} = \widehat{\Sigma}^s\circ \mathbf{1}_{I\cup E(G_\mathrm{path}^r)}, \quad s = 0, \ldots, T \end{align}\] Here, \(I\) denotes the diagonal entries, \(E(G^r_\mathrm{path})\) denotes the retained off-diagonal entries, \(\circ\) denotes entrywise multiplication, and \(\mathbf{1}_{I\cup E(G^r_\mathrm{path})}\) is the corresponding mask. The reported error is \[\begin{align} \frac{\|\widehat{\Sigma}^s_{G_\mathrm{path}^r}-\Sigma_{t_s}\|_\mathrm{F}}{\|\Sigma_{t_s}\|_\mathrm{F}}, \quad s = 0, \ldots, T, \end{align}\] where the ground truth \(\Sigma_{t_s}\) is computed from 99 . This metric captures both the recovery error on retained entries and the truncation error from omitted covariance entries. When \(r=d-1\), the graph \(G_\mathrm{path}^r\) is complete and the full covariance matrix is recovered from the relaxation.

1 shows that the moment relaxation accurately recovers the intermediate statistics at the prescribed grid times and the affine velocity field.

Next, we study the effect of the graph radius \(r\). Keeping the source and target distributions unchanged from the previous setup, we vary \(r\) from \(6\) to \(20\) in increments of \(2\). For each \(r\), we solve the corresponding relaxation associated with the cluster graph \(G^r_\mathrm{path}\) and evaluate the covariance and velocity coefficient errors at \(t=0.2\) and \(t=0.5\).

2 shows that increasing the graph radius generally improves the recovered covariance and velocity coefficients. This behavior is consistent with the correlation decay of the Gaussian distribution: a larger graph radius retains more of the relevant covariance and cross-covariance structure.

Finally, we compare our approach with a particle-based back-propagation method. In this experiment, we choose dimension \(d = 15\). The initial and terminal Gaussian distributions are generated as above with precision bandwidth \(3\) and diagonal shift \(\delta = 5\). The source and target means are sampled independently as \[\begin{align} m_0 \sim \mathcal{N}(0,I_d), \quad m_1 \sim \mathcal{N}(0.5 \, \mathbf{1}_d, I_d). \end{align}\] To ensure a fair comparison, rather than using the closed-form Gaussian moments from previous experiments, we estimate the initial and terminal moments for the moment relaxation using \(N_{\text{emp}} = 10{,}000\) independent samples drawn from each Gaussian distribution. The relaxation is then solved over the complete graph \(G_\mathrm{path}^{d-1}\).

For the back-propagation method, we draw \(N_{\text{train}}=10,000\) source particles \(\left\{ X_i^0 \right\}_i \sim \mu\). We use the affine form of the Gaussian velocity field and parameterize the velocity field by \[\begin{align} v_\theta(x,t) = A_\theta(t)x + c_\theta(t). \end{align}\] The time-dependent coefficients are expanded in the degree-two monomial basis \(\left\{ 1, t, t^2 \right\}\): \[\begin{align} A_\theta(t)=A_0+tA_1+t^2A_2, \quad c_\theta(t)=c_0+tc_1+t^2c_2, \end{align}\] where \(A_k \in \mathbb{R}^{d\times d}\) and \(c_k \in \mathbb{R}^d\) are trainable parameters. The particles are pushed forward by explicit Euler steps \[\begin{align} X_i^{s+1} = X_i^s + \Delta t\, \bigl(A_\theta(t_s)X_i^s+c_\theta(t_s)\bigr). \end{align}\] The empirical kinetic cost is \[\mathcal{K}(\theta) = \sum_{s=0}^{T-1} \Delta t\, \frac{1}{N_{\text{train}}} \sum_{i=1}^{N_{\text{train}}} \left\| A_\theta(t_s)X_i^s+c_\theta(t_s)\right\|^2 .\] The terminal constraint at \(t=1\) is enforced by matching Hermite moments of degree up to two. Let \(\Psi(x)\) denote the Hermite basis vector centered and normalized with respect to the estimated target Gaussian distribution. The terminal residual is \[r(\theta) = \frac{1}{N_{\text{train}}} \sum_{i=1}^{N_{\text{train}}} \Psi(X_i^T) - \mathbb{E}_{Y\sim\nu}\Psi(Y),\] where the target expectation is estimated by Monte Carlo sampling from the target Gaussian distribution. While various penalties for the terminal constraint can be added to the kinetic cost, the moment-matching penalty presented here yields the best empirical performance. To further refine the solution, we employ an augmented Lagrangian method. The Lagrangian takes the form: \[\mathcal{L}_\beta(\theta,\lambda) = \mathcal{K}(\theta) + \lambda^\top r(\theta) + \frac{\beta}{2} \|r(\theta)\|_2^2 .\] In the reported run, the penalty parameter is initialized as \(\beta=5\) and kept fixed. The method performs \(10\) outer augmented-Lagrangian iterations; in each outer iteration, the velocity parameters are optimized for \(600\) Adam steps with learning rate \(5\times 10^{-3}\). We use a gradient-clipping threshold of \(100\), and a stepwise learning-rate scheduler with decay factor \(0.6\) every \(200\) inner iterations. After each inner solve, the dual coefficients are updated by \(\lambda \leftarrow \lambda+\beta \, r(\theta)\). After training, we push forward a newly generated set of \(N_{\text{test}} = 10,000\) particles using the learned velocity field, and compare the empirical covariance matrices along the path with the exact Gaussian geodesic covariance. We then compare these errors with the covariance error obtained directly from the moment relaxation.

3 shows that the moment relaxation achieves a uniformly smaller covariance error along the trajectory in this experiment. The resulting SDP is also solved within a few seconds, substantially faster than the particle back-propagation training procedure.

5.2 Benamou–Brenier dynamics between more general distributions↩︎

We next consider the Benamou–Brenier dynamics from a Gaussian distribution to a one-dimensional lattice Ginzburg–Landau distribution: \[\begin{align} \nu(y) \propto \exp\left[-\beta\left(\sum_{i=1}^{d+1}\frac{\lambda}{2}\left(\frac{y_i-y_{i-1}}{h}\right)^2+\frac{1}{4\lambda}(1-y_i^2)^2\right)\right],\quad y\in[-L,L]^d, \quad h = \frac{1}{d+1} \end{align}\] with boundary conditions \(y_0 = y_{d+1}=0\). The first term is a nearest-neighbor interaction, while the second is a double-well potential. Thus, the target distribution has local spatial correlation.

We set \(d=10\), \(\beta=1/8\), \(\lambda=0.03\), \(T=5\), and \(L=2.5\). We use an equispaced time grid on \([0,1]\) with \(\Delta t=1/T\). We generate \(N_{\text{emp}}=50{,}000\) samples from \(\nu\) by tensor-train conditional sampling [30], [31]. Let \(\widehat{m}_\nu\) and \(\widehat{\Sigma}_\nu\) be the empirical mean and covariance of the training samples from \(\nu\). We choose the source distribution to be a Gaussian distribution \[\mu = \mathcal{N}(\widehat{m}_\nu, 0.7 \,\widehat{\Sigma}_\nu).\] Thus, the source and target have comparable locations and covariance scales. We draw another \(N_{\text{emp}}=50{,}000\) samples from \(\mu\) to estimate the initial moments and solve the Benamou–Brenier path from \(\mu\) to \(\nu\) using the moment relaxation 67 . We use the single-coordinate clusters \(\mathcal{C}= \left\{ \left\{ 1 \right\}, \ldots, \left\{ d \right\} \right\}\) and a path graph \(G_\mathrm{path}\) on \([d]\) as the cluster graph \(\mathcal{G}\). To recover moments for a non-neighboring cluster pair, such as the pair involving \(x_1\) and \(x_5\), we additionally include the edge \(\{1,5\}\) in \(\mathcal{E}\).

We use a monomial cluster basis. For each coordinate, we take \(\phi_{i,j}(x_i)=x_i^j\). Unlike the Gaussian case, where the path is fully determined by the first and second moments, the Ginzburg–Landau target requires higher-order moments to represent its non-Gaussian structure. For a cluster \(A = \left\{ i \right\}\), the monomial cluster basis is \[\begin{align} \Phi_A^{\mathrm{mon}}(z_A) = \left\{ x_i^{\alpha_x}y_i^{\alpha_y}: \alpha_x + \alpha_y \le \deg_{\mathrm{mon}} \right\}. \end{align}\] In the experiments below, we use \(\deg_{\mathrm{mon}} = 11\).

Using the above basis and parameters, we solve 67 . The problem is a multi-block SDP, with one real symmetric PSD block for each coupling \(\pi^s\). To exploit sparsity, we apply chordal decomposition directly to each time-block PSD constraint. The resulting SDPs are solved by MOSEK with a tolerance of \(10^{-10}\). After solving the chordally decomposed SDP, the clique-level primal matrices are reassembled into the original moment matrices by averaging overlapping clique entries and applying PSD completion. These completed matrices are then used to extract moments of intermediate distributions. The resulting objective values and implementation statistics are summarized in 1. The static monomial method is the static cluster moment relaxation for optimal transport from [25] and is included as a reference baseline for the endpoint transport problem. The dynamic monomial row corresponds to the moment relaxation 67 applied to problem 74 .

3pt

We reconstruct the intermediate low-dimensional marginals from the recovered moments by a maximum-entropy step. Let \(M^s_{AB} (\rho^s)\) denote the recovered moments of the \(AB\)-marginal of the distribution \(\rho^s\). These moments are extracted from the left- or right-layer moments of the coupling \(\pi^s\). At the continuous level, the maximum-entropy reconstruction solves \[\begin{align} \max_{\rho^s_{AB}} \quad & -\int \log \frac{\mathop{}\!\mathrm{d}\rho_{AB}^s }{\mathop{}\!\mathrm{d}x_{AB}} \; \mathop{}\!\mathrm{d}\rho^s_{AB} \\ \text{s.t.}\quad& \rho^s_{AB}(R_x(\Phi_A \Phi_B^*)) = M^s_{AB} (\rho^s), \nonumber\label{constraint-moment32ME}\\ & \rho_{AB}^s([-L,L]^{|AB|}) = 1.\end{align}\tag{101}\] Here \(\rho^s_{AB}\) is optimized over all nonnegative measures on \([-L,L]^{|AB|}\). In the computation, we solve a discretized version on a grid and relax the moment-matching constraint 101 by \[\begin{align} \|\rho^s_{AB}(R_x(\Phi_A \Phi_B^*)) - M^s_{AB} (\rho^s)\|_{\infty}\le \varepsilon_{\mathrm{mom}}\max\{1,\|M^s_{AB} (\rho^s)\|_{\infty}\}. \end{align}\] We use \(\varepsilon_{\mathrm{mom}}=10^{-2}\) and \(21\) grid points per coordinate. In 4, the resulting \(2\)-dimensional marginals are visualized by kernel density estimation on \([-L,L]^2\) with \(120\) grid points per coordinate. Here we do not use all moments available from the degree-\(11\) SDP relaxation. Although the relaxation contains monomial moments up to order \(22\), high-degree monomial moments are poorly scaled and can dominate the moment-matching constraints. Therefore we only use recovered monomial moments with total degree \(|\alpha| \le 4\) in the maximum-entropy reconstruction.

4 shows that the maximum-entropy reconstructions from the recovered moments capture the gradual deformation from the Gaussian source to the Ginzburg–Landau target. The neighboring \((1,2)\)-marginal captures local correlations induced by the interaction term, while the non-neighboring \((1,5)\)-marginal illustrates how the relaxation can recover selected longer-range low-dimensional statistics when the corresponding cluster-pair edge is retained.

Finally, we use the dual slack matrices to compute the velocity field \(\left\{ \widehat{v}(\, \cdot \,, t_s) \right\}_{s=0}^{T-1}\) as discussed in 3.3. We then push forward test samples from \(\mu\) by explicit Euler steps. On each interval \([t_s,t_{s+1}]\), we use five Euler substeps: \[\begin{align} X^{s,q+1}=X^{s,q}+\frac{\Delta t}{5}\widehat v(X^{s,q}, t_s),\qquad q=0,\ldots,4. \end{align}\] The final value \(X^{s,5}\) is used as the initial value for the next time interval. To stabilize the visualization, we only transport test samples satisfying \(\|X\|_{\infty}\le 2.2\) and clip the transported samples to \([-L,L]^d\).

In 5, we compare the empirical \((1,2)\)- and \((1,5)\)-marginals across three cases: the source samples, the target samples, and the push-forward samples generated by the recovered velocity fields at \(t = 1\). This demonstrates that the computed dual solution of the moment relaxation provides an implementable approximate flow, not merely intermediate moment statistics.

5.3 Single-spin update process between Ising models↩︎

We finally consider the constrained Markov process between Ising models introduced in 4. The Ising model is specified in 87 . We write \((\beta_\mu, h_\mu, J_\mu)\) and \((\beta_\nu, h_\nu, J_\nu)\) for the parameters of the source law \(\mu\) and target law \(\nu\), respectively. In the experiments below, the initial and terminal laws share the same inverse temperature \(\beta_\mu = \beta_\nu = 0.4\), but have opposite uniform edge interactions \(J_{\mu,ij} = 1\) and \(J_{\nu,ij} = -1\) on the chosen Ising graph. Thus \(\mu\) is ferromagnetic and favors aligned neighboring spins, whereas \(\nu\) is antiferromagnetic and favors alternating neighboring spins. The external fields are set to zero: \(h_\mu = h_\nu = 0\). We generate \(N_{\text{emp}}=20{,}000\) samples \(\left\{ X_r \right\}_{r=1}^{N_{\text{emp}}} \sim \mu\) and \(\left\{ Y_r \right\}_{r=1}^{N_{\text{emp}}} \sim \nu\) using Glauber dynamics, and impose the initial and terminal constraints using the empirical distributions.

We consider two types of Ising graphs \(\mathcal{G}_\mathrm{Ising}\). The first is the one-dimensional path graph on \([d]\), with edge set \(\mathcal{E}_{\mathrm{Ising}} = \left\{ \{i,i+1\}:i = 1, \ldots, d-1 \right\}\). For this \(1\)D Ising model, we choose \(d=30\) and cluster size \(n_c=2\). We decompose the coordinate set into contiguous clusters \[A_k=\{(k-1)n_c+1,\ldots,\min(k\,n_c,d)\}, \quad k = 1, \ldots, \left\lceil{d/n_c}\right\rceil.\] The cluster graph is chosen to be the path graph on these contiguous clusters.

The second example is a \(2\)D Ising model on a \(d_x\times d_y\) grid. In the experiment below, we take \(d_x=d_y=4\), so the total number of spins is \(d=16\). We partition the lattice into four \(2 \times 2\) clusters and choose the cluster graph to be a cycle, as shown in 6.

For both the \(1\)D and \(2\)D models, the source and target Ising parameters are chosen as above. We take the number of time steps \(T=5d\), corresponding to five complete sweeps over the spin system. The active spin at step \(s\) is chosen cyclically as \[i_s = 1 + (s \bmod d), \quad s=0, \ldots, T-1.\] We use the quadratic spin-flip cost \(c^s(x,y) = \|x-y\|^2\), which equals \(4\) if the active spin flips and \(0\) otherwise. We solve the problem using the marginal relaxation 31 , which becomes a linear program. In the implementation, we use a reduced version adapted to the single-spin schedule. At time step \(s\), coupling variables are introduced only for clusters \(A\) with \(i_s \in A\) and retained cluster-pairs \(AB\) with \(i_s \in AB\). Marginals of inactive blocks are propagated from the most recent time at which the corresponding block was updated. This reduced construction gives the same relevant local marginals for the single-spin schedule while substantially reducing the LP size. Solving the LP yields the recovered marginals, denoted by \(\{\rho_A^s\}_{A\in\mathcal{C}}\) and \(\{\rho_{AB}^s\}_{AB\in\mathcal{E}}\) for \(s=0,\ldots,T\).

For the \(1\)D case, after solving the relaxation, we track the two-spin marginal on sites \(\{1,2\}\) after each epoch. Here, a single epoch corresponds to a complete sweep of all spins. 7 shows the evolution of the \((1,2)\)-marginal. The intermediate heatmaps display a gradual transfer of probability mass from aligned configurations to anti-aligned configurations.

After obtaining the local marginals from the relaxation, we fit Glauber dynamics using the linear regression form 98 . For the \(1\)D experiment, we take the Glauber dependency graph to be the path graph on \([d]\). For an interior active spin \(i=i_s\), the parameters \((h_i^s, J_{i,i-1}^s, J_{i,i+1}^s)\) are fitted from the recovered next-time local marginal on \(B_i = \left\{ i-1,i,i+1 \right\}\) by weighted ridge regression. Since \(B_i\) is contained in a retained cluster pair, its marginal can be obtained directly by projection from the relaxed solution. For boundary spins, the same procedure is used with \(B_1 = \left\{ 1,2 \right\}\) and \(B_d = \left\{ d-1,d \right\}\). The fitted parameters define a sequence of Glauber kernels \(\{K_{\widehat{\theta}_s}^s\}_{s=0}^{T-1}\). Starting from the same source samples \(\left\{ X^0_r \right\}_{r=1}^{N_{\mathrm{emp}}}\), we generate particles by \[\begin{align} X^{s+1}_r\sim K_{\widehat{\theta}_s}^s(\cdot\mid X^s_r), \qquad s=0,\ldots,T-1. \end{align}\] From the generated particles, we compute empirical local marginals and compare them with the marginals obtained from the relaxation.

8 compares the recovered marginals from the marginal relaxation and the empirical marginals of the fitted Glauber dynamics on the cluster \(\left\{ 1,2 \right\}\). The bars almost overlap throughout the process, showing that the fitted dynamics closely reproduces the recovered two-spin marginal on this block.

We also visualize a single particle trajectory generated by the fitted dynamics. 9 shows snapshots of one representative particle after each epoch. The trajectory changes by single-spin updates along the cyclic schedule, illustrating one realization of the fitted Markov process whose empirical marginals match the recovered marginals.

For the \(2\)D case, we use the same initial and terminal parameters and the \(2 \times 2\) cluster structure described above. The marginal relaxation is solved using the same single-spin schedule. In the fitting step, we use a dense local graph for Glauber dynamics induced by the retained cluster-pair supports: if \(AB \in \mathcal{E}\) is an edge in the cluster graph, then all spin pairs inside \(AB\) are included in \(\mathcal{G}_\mathrm{glb}\). The local marginal required for fitting the update of spin \(i_s\) is not always directly available from a single retained cluster or cluster-pair marginal, so we approximate it by the Bethe-type reconstruction. 10 compares the recovered marginals from the marginal relaxation and the empirical marginals generated by the fitted Glauber dynamics on the cluster \(A_1 = \left\{ 1,2,5,6 \right\}\). The fitted process reproduces the recovered local marginal on \(A_1\) well across the displayed epochs.

6 Conclusion and discussion↩︎

In this work, we extended the convex-relaxation framework developed for optimal transport in [25] to the optimization of Markov processes via a time-discrete sequential-coupling formulation. The proposed relaxations avoid explicit representation of the full joint law while retaining local statistics of the intermediate distributions, thereby mitigating the curse of dimensionality. In the unconstrained kinetic-cost setting, the framework recovers the Benamou–Brenier geodesic on the prescribed time grid, and the dual variables provide a practical way to reconstruct approximate velocity fields. Beyond the Benamou–Brenier setting, we showed how the recovered local statistics can be converted into implementable dynamics by fitting a parametric transition family, as illustrated by Glauber dynamics for finite spin systems.

Several directions remain open. It would be important to characterize conditions under which the locality of the initial and terminal distributions is preserved along the interpolating process. Such a result would directly justify sparse marginal and cluster-moment relaxations at intermediate times. Another direction is to extend the framework from Markov processes with prescribed transition structures to broader classes of controlled dynamics. Rotation-constrained transport is one representative example: instead of moving mass through an unconstrained velocity field, one seeks to steer a distribution through structured transformations. Incorporating such control variables into the relaxation would broaden the scope of the convex framework for process optimization.

References↩︎

---

1. Department of Statistics, University of Chicago, (hongyi518@uchicago.edu).↩︎

2. Department of Statistics, University of Chicago, (ykhoo@uchicago.edu). The research of this author is partially funded by NSF DMS-2339439, DOE DE-SC0022232, DARPA The Right Space HR0011-25-9-0031, and a Sloan research fellowship.↩︎

3. Department of Statistics, University of Chicago, (ttang@u.nus.edu).↩︎