Score-based Neural Ordinary Differential Equations for Computing Mean Field Control Problems

## Abstract

Abstract Classical neural ordinary differential equations (ODEs) are powerful tools for approximating the log-density functions in high-dimensional spaces along trajectories, where neural networks parameterize the velocity fields. This paper proposes a system of neural differential equations representing first- and second-order score functions along trajectories based on deep neural networks. We reformulate the mean field control (MFC) problem with individual noises into an unconstrained optimization problem framed by the proposed neural ODE system. Additionally, we introduce a novel regularization term to enforce characteristics of viscous Hamilton–Jacobi–Bellman (HJB) equations to be satisfied based on the evolution of the second-order score function. Examples include regularized Wasserstein proximal operators (RWPOs), probability flow matching of Fokker–Planck (FP) equations, and linear quadratic (LQ) MFC problems, which demonstrate the effectiveness and accuracy of the proposed method.

## 1 Introduction

Score functions have been widely used in modern machine learning algorithms, particularly generative models through time-reversible diffusion (Song et al., ). The score function can be viewed as a deterministic representation of diffusion in stochastic trajectories (Carrillo et al., ). In this representation, one reformulates the Brownian motion by the gradient of the logarithm of the density function, after which deterministic trajectories involving score functions can approximate the probability density function. These properties have inspired algorithms for simulating stochastic trajectories or sampling problems that converge to target distributions (Wang et al., ; Lu et al., ). Typical applications include modeling the time evolution of probability densities for stochastic dynamics and solving control problems constrained by such dynamics.

While score functions provide powerful tools for modeling stochastic trajectories, their computations are often inefficient, especially in high-dimensional spaces. Classical methods, such as kernel density estimation (KDE) (Chen, ), tend to perform poorly in such settings due to the curse of dimensionality (Terrell & Scott, ).

Recently, neural ODEs have emerged as efficient ways of estimating densities. In particular, one uses neural networks to parameterize the velocity fields and then approximates the logarithm of density function along trajectories. The time discretizations of neural ODEs can be viewed as normalization flows in generative models.

Several natural questions arise. Can we approximate the trajectory of score functions by constructing a set of neural ODEs (in continuous time) or normalization flows (in discrete time)? Furthermore, can we efficiently solve stochastic control problems with these neural ODEs or normalization flows?

In this paper, we propose a formulation for the first- and second-order score functions using two additional neural ODEs involving high-order derivatives of the velocity fields. We also develop a class of high-order normalizing flows in the discrete-time update of these neural ODEs. As an application, we use this high-order normalizing flow to solve the MFC problem with individual noises. The MFC problem generalizes stochastic control problems by incorporating a running cost that depends on the state distribution, thus capturing the interaction between individual agents and the population density. In our reformulation, the MFC problem is recast into an optimal control problem involving state trajectories, densities, and score functions, which can be efficiently approximated using the proposed neural ODE system. Additionally, a regularization term from the Karush-Kuhn-Tucker (KKT) system of the MFC problem is proposed to enhance the optimization. In this regularization, we approximate the viscous HJB equations using the first- and second-order score dynamics. Numerical examples in various MFC problems demonstrate the effectiveness of the proposed optimization method with high-order neural differential equations.

Related work. The normalizing flow has emerged as a powerful technique for solving inference problems, allowing for the construction of complex probability distributions in a tractable manner. This approach was popularized by works such as Rezende & Mohamed (); Papamakarios et al. (). The advent of Neural ODEs (Chen et al., ) has attracted considerable attention in the community, opening up a new paradigm for continuous-time deep learning models. Following this breakthrough, numerous extensions have been proposed to enhance the expressiveness and flexibility of these models, such as the augmented Neural ODEs by Dupont et al. ().

Regarding score-based models, Song et al. () studies the score-based reversible time diffusion models in generative modeling. In this case, the score function often comes from the OU process. One needs to solve the score-matching problem to learn the time-dependent score function. In addition, the first-order score dynamic from neural ODEs has been introduced by Boffi & Vanden-Eijnden () to work on the simulation of the Fokker-Planck (FP) equation. Rather than using score dynamics directly to compute the probability flow ODEs associated with FP equations, they employed score-matching methods to achieve efficient computations.

The study of MFC problems has become crucial in the last decade Cardaliaguet et al. (); Bensoussan et al. (); Fornasier & Solombrino (). MFC studies strategic decision-making in large populations where individual players interact through specific mean-field quantities. This control formulation is also useful in generative models Zhang & Katsoulakis (). For example, neural network-based methods have been employed to solve MFC problems. Hu & Lauriere () provided a comprehensive review of neural network approaches for MFC. For first-order MFC problems, Ruthotto et al. (); Huang et al. () developed Neural ODE-based methods. In parallel, Ruthotto et al. () introduced numerical algorithms based on the characteristic lines of the Hamilton-Jacobi-Bellman (HJB) equation. For second-order MFC problems, neural networks have also been employed. Lin et al. () utilized generative adversarial networks (GANs) to approximate minimization systems for MFC problems, where one neural network models the population density, and the other parameterizes the value function. Reisinger et al. () proposed a PDE-based iterative algorithm to tackle MFC problems. Along this line, MFC problems are generalizations of dynamical optimal transport problems, known as Benamou-Brenier’s formulas Benamou & Brenier (). A famous example is the Jordan-Kinderlehrer-Otto (JKO) scheme (Jordan et al., ), a variational time discretization in approximating the gradient drift FP equation. The one-step iteration of the JKO scheme can be viewed as an MFC problem. In machine learning computations, the neural JKO scheme (Xu et al., ; Vidal et al., ) was introduced to compute the FP equation, and Lee et al. () extended this to approximate general nonlinear gradient flows through a generalized deep neural JKO scheme. This scheme approximates a deterministic MFC problem in each time interval. Other works that use machine learning methods to solve MFC problems include Dayanikli et al. (); Zhou et al. (); Dayanikli et al. (). Different from previous results, our work designs first- and second-order score dynamics to compute second-order MFC problems.

The organization of this paper is as follows. In section 2, we formulate the neural dynamical system involving first- and second-order score functions, with a specific example based on the Gaussian distribution in terms of linear mapping equations. Section 3 discusses the discrete time evolution of the neural ODE system, which can be viewed as a system of normalization flows for approximating first- and second-order score functions. In section 4, we design an algorithm using a score-based ODE system to solve the second-order MFC problem. Several numerical examples are presented in section 5 to demonstrate the accuracy and effectiveness of the proposed algorithm.

## 2 Score-based neural ODE in continuous time

In this section, we propose a neural dynamical system that evolves along a trajectory of a random variable. The system outputs the logarithm of the density function, as well as its gradient vector and Hessian matrix. We show that these neural dynamics are crucial for efficiently computing the score functions in high-dimensional density estimation problems, which are key in generative models and MFC problems.

We first clarify the notations. Denote a variable $x\in\mathbb{R}^{d}$.
$\nabla_{x}$ denotes the gradient or Jacobian matrix of a function w.r.t. the variable $x$. The gradient is always a column vector. $\nabla_{x}\cdot$ and $\nabla_{x}^{2}$ are divergence and Hessian w.r.t. the variable $x$. $\left|\cdot\right|$ is the absolute value of a scalar, the $l_{2}$ norm of a vector, or the Frobenius norm of a matrix. $\operatorname{Tr}(\cdot)$ denotes the trace of a squared matrix.

Let $z_{t}\in\mathbb{R}^{d}$ be a flow of random variable for $t\geq 0$ with initialization $z_{0}=x$, which is obtained from a pushforward map $z_{t}=T(t,x)$ in Lagrangian coordinates. For simplicity, we denote $T_{t}:=T(t,\cdot)$, so that $T_{0}$ is an identity map. Throughout this paper, we assume that $T(\cdot,\cdot)$ is smooth and $T_{t}$ is invertible for all $t\geq 0$. We denote the probability distribution of $z_{t}$ by $\rho_{t}$ or $\rho(t,\cdot)$, then $\rho_{t}=T_{t\#}\rho_{0}$, where # is the pushforward operator for distributions. By the change of variable formula, $\rho$ satisfies the well known Monge–Ampère (MA) equation (Gutiérrez & Brezis, )

$$ $\rho(t,T(t,x))\,\det(\nabla_{x}T(t,x))=\rho(0,x)\,.$ (1) $$

Let $f(t,\cdot):\mathbb{R}^{d}\to\mathbb{R}^{d}$ be the vector field of the state dynamic $z_{t}$ in Eulerian coordinates, given by $f(t,T(t,x))=\partial_{t}T(t,x)$ for all $t$ and $x$. Then, the state $z_{t}$ satisfies the ODE dynamic $\partial_{t}z_{t}=f(t,z_{t})$. The density function hence satisfies the continuity equation (also known as the transport equation):

$$ $\partial_{t}\rho(t,x)+\nabla_{x}\cdot(\rho(t,x)f(t,x))=0\,.$ (2) $$

We leave the proofs for this equation and all propositions afterwards in Appendix A.

We denote the logarithm of density along the state trajectory $z_{t}$ by $l_{t}:=\log\rho(t,z_{t})$. We also denote the first- and second-order score functions along the state trajectory by $s_{t}:=\nabla_{z}\log\rho(t,z_{t})$ and $H_{t}:=\nabla_{z}^{2}\log\rho(t,z_{t})$. These score functions are useful in density estimation problems, particularly for applications such as generative models and MFC problems. However, their efficient computations in high-dimensional spaces are challenging problems. In this work, we propose a system of high-order neural ODEs to compute these score functions, formalized in the following proposition.

###### Proposition 1 (Neural ODE system) .

Functions $z_{t}$, $l_{t}$, $s_{t}$, and $H_{t}$ satisfy the following ODE dynamics.

$$ $\displaystyle\partial_{t}z_{t}$ $\displaystyle=f(t,z_{t})\,,$ (3a) $\displaystyle\partial_{t}l_{t}$ $\displaystyle=-\nabla_{z}\cdot f(t,z_{t})\,,$ (3b) $\displaystyle\partial_{t}s_{t}$ $\displaystyle=-\nabla_{z}f(t,z_{t})^{\top}s_{t}-\nabla_{z}(\nabla_{z}\cdot f(t ,z_{t}))\,,$ (3c) $\displaystyle\partial_{t}H_{t}$ $\displaystyle=-\sum_{i=1}^{d}s_{it}\nabla_{z}^{2}f_{i}(t,z_{t})-\nabla_{z}^{2} (\nabla_{z}\cdot f(t,z_{t}))-H_{t}\nabla_{z}f(t,z_{t})-\nabla_{z}f(t,z_{t})^{ \top}H_{t}\,,$ (3d) $$

where $s_{it}$ and $f_{i}$ are the $i$-th component of $s_{t}$ and $f$ respectively.

One corollary of this proposition is that, if we denote the density along the trajectory by $\tilde{l}_{t}:=\rho(t,z_{t})$, then $\tilde{l}_{t}$ satisfies the following ODE:

$$ $\partial_{t}\tilde{l}_{t}=-\nabla_{z}\cdot f(t,z_{t})\,\tilde{l}_{t}\,.$ (4) $$

We note that the first and second order score functions satisfy the following information equality.

###### Proposition 2 (Information equality) .

The following equality holds for all $t\geq 0$,

$$ $\mathbb{E}\left[\operatorname{Tr}(H_{t})\right]=-\mathbb{E}\left[|s_{t}|^{2} \right]\,.$ $$

Example: Centered Gaussian distributions.
To illustrate the ODE dynamics equation 3, we consider a concrete example where the random variable follows a centered Gaussian distribution. Let the pushforward map be linear in $x$, i.e. $T(t,x)=T(t)x$, where $T(t):\mathbb{R}_{+}\to\mathbb{R}^{d\times d}$ is a time-dependent matrix. If we define $A(t):\mathbb{R}_{+}\to\mathbb{R}^{d\times d}$ by $\partial_{t}T(t)=A(t)T(t)$, then the vector field of the state dynamics is $f(t,x)=A(t)x$. Now, let the initial state $z_{0}$ follow a centered Gaussian distribution $N(0,\Sigma(0))$ with covariance matrix $\Sigma(0)$. Under this linear map $T(t)$, the state $z_{t}$ remains a Gaussian distribution, with a time-dependent covariance matrix $\Sigma(t)$. In this case, $z_{t}\sim N(0,\Sigma(t))$, and the time evolution of the covariance matrix follows the matrix ODE (with the derivation in Appendix A.1)

$$ $\partial_{t}\Sigma(t)=A(t)\Sigma(t)+\Sigma(t)A(t)^{\top}\,.$ (5) $$

In this case, $\nabla_{x}f(t,x)=A(t)$ and $\nabla_{x}\cdot f(t,x)=\operatorname{Tr}(A(t))$, and all higher-order spatial derivatives of $f$ vanish. Therefore, the ODE system equation 3 simplifies to

$$ $\displaystyle\partial_{t}z_{t}$ $\displaystyle=A(t)z_{t}\,,\quad\quad\quad\partial_{t}l_{t}=-\operatorname{Tr}( A(t))\,,$ (6) $\displaystyle\partial_{t}s_{t}$ $\displaystyle=-A(t)^{\top}s_{t}\,,\quad\quad\partial_{t}H_{t}=-H_{t}A(t)-A(t)^ {\top}H_{t}\,.$ (7) $$

###### Proposition 1 (Neural ODE system) .

###### Proposition 2 (Information equality) .

## 3 Score-based normalization flows

In this section, we introduce the time discretization of the proposed neural ODE system. This discretized system can be interpreted as a normalization flow. This means that we can efficiently estimate the first- and second-order score functions using a deep neural network function.

### 3.1 A generalization of normalizing flow

Chen et al. () interpreted the deep residual neural network as an ODE with each layer representing one step of the forward Euler scheme. Similar ideas have also been proposed by Haber & Ruthotto (). In this section, we extend these concepts and demonstrate how the first- and second-order score functions can be computed efficiently using the proposed neural ODE system.

We partition the time interval $[0,t_{\text{end}}]$ into $N_{t}$ subintervals with the length $\Delta t=t_{\text{end}}/N_{t}$ and denote the time stamps by $t_{j}=j\,\Delta t$. We apply the forward Euler scheme to discretize the ODE system equation 3. We parametrize the vector field $f(t,z;\theta)$ as a neural network with parameter $\theta$. Here, $\theta=[w_{0},w_{1},w_{2},b_{1},b_{2}]\in\mathbb{R}^{k+k\times d+d\times k+k+d}$ and

$$ $f(t,z;\theta)=w_{2}\,\sigma(w_{0}t+w_{1}z+b_{1})+b_{2}\,,$ (8) $$

where $\sigma\colon\mathbb{R}^{k}\rightarrow\mathbb{R}^{k}$ is a vector function with elementwise activation functions. Here $k$ is the width of the network. This is the typical structure of a neural network with one hidden layer. Then, we discretize the state dynamic $z_{t}$ through

$$ $z_{t_{j+1}}=z_{t_{j}}+\Delta t\,f(t_{j},z_{t_{j}};\theta)\,.$ (9) $$

Similarly, numerical simulations for functions $l_{t}$, $\tilde{l}_{t}$, $s_{t}$, and $H_{t}$ are given by

$$ $\displaystyle l_{t_{j+1}}$ $\displaystyle=l_{t_{j}}-\Delta t\,\nabla_{z}\cdot f(t_{j},z_{t_{j}};\theta)\,,$ (10) $\displaystyle\tilde{l}_{t_{j+1}}$ $\displaystyle=\tilde{l}_{t_{j}}-\Delta t\,\nabla_{z}\cdot f(t_{j},z_{t_{j}}; \theta)\,\tilde{l}_{t_{j}}\,,$ (11) $\displaystyle s_{t_{j+1}}$ $\displaystyle=s_{t_{j}}-\Delta t\left(\nabla_{z}f(t_{j},z_{t_{j}};\theta)^{ \top}s_{t_{j}}+\nabla_{z}(\nabla_{z}\cdot f(t_{j},z_{t_{j}};\theta))\right)\,,$ (12) $\displaystyle H_{t_{j+1}}$ $\displaystyle=H_{t_{j}}-\Delta t\left(\sum_{i=1}^{d}s_{it_{j}}\nabla_{z}^{2}f_ {i}(t_{j},z_{t_{j}};\theta)+\nabla_{z}^{2}(\nabla_{z}\cdot f(t_{j},z_{t_{j}}; \theta))\right.$ (13) $\displaystyle\hskip 72.26999pt+H_{t_{j}}\nabla_{z}f(t_{j},z_{t_{j}};\theta)+ \nabla_{z}f(t_{j},z_{t_{j}};\theta)^{\top}H_{t_{j}}\bigg{)}\,,$ (14) $$

where the derivatives of $f$ are obtained from auto-differentiations of the neural network.

In this way, the map from $z_{0}$, $l_{0}$, $\tilde{l}_{0}$, $s_{0}$, $H_{0}$ to $z_{t_{j}}$, $l_{t_{j}}$, $\tilde{l}_{t_{j}}$, $s_{t_{j}}$, $H_{t_{j}}$, for any $j\geq 2$, can be viewed as deep residual neural networks, with each layer given by a forward Euler time step. For example,

$$ $\displaystyle z_{t_{N_{t}}}$ $\displaystyle=\left(id+\Delta t\,f(t_{N_{t}-1},\cdot;\theta)\right)\circ\cdots \circ\left(id+\Delta t\,f(t_{0},\cdot;\theta)\right)(z_{0})\,,$ (15) $\displaystyle l_{t_{N_{t}}}$ $\displaystyle=\left(id-\Delta t\,\nabla_{z}\cdot f(t_{N_{t}-1},z_{t_{N_{t}-1}} ;\theta)\right)\circ\cdots\circ\left(id-\Delta t\,\nabla_{z}\cdot f(t_{0},z_{0 };\theta)\right)(l_{0})\,,$ $\displaystyle s_{t_{N_{t}}}$ $\displaystyle=\left(id-\Delta t\left(\nabla_{z}f(t_{N_{t-1}},z_{N_{t-1}}; \theta)^{\top}\cdot+\nabla_{z}(\nabla_{z}\cdot f(t_{N_{t-1}},z_{N_{t-1}}; \theta))\right)\right)\circ\cdots\circ$ $\displaystyle\qquad\left(id-\Delta t\left(\nabla_{z}f(t_{0},z_{0};\theta)^{ \top}\cdot+\nabla_{z}(\nabla_{z}\cdot f(t_{0},z_{0};\theta))\right)\right)(s_{ 0})\,,$ $$

where $id$ is the identity mapping function and $\circ$ represents compositions of functions.

### 3.2 Different algorithms for computing the score function

The expression equation 3c or equation 12 provides a way to compute the score function. In this section, we compare several algorithms for computing the score function, and demonstrate the potential efficiency of using formulas equation 3c or equation 12.

One method is to derive the score function from the pushforward map $T(t,x)$ and the MA equation directly. To be more specific, we apply the operator $\nabla_{x}$ and the logarithm on both sides of the MA equation equation 1. We obtain

$$ $\nabla_{x}T(t,x)^{\top}\nabla_{T}\log\rho(t,T(t,x))+\nabla_{x}\log\left(\det( \nabla_{x}T(t,x))\right)=\nabla_{x}\log\rho(0,x)\,.$ (16) $$

Recall that $s_{t}=\nabla_{T}\log\rho(t,T(t,x))$, so the score function can be computed through

$$ $s_{t}=\nabla_{x}T(t,x)^{-\top}\left[\nabla_{x}\log\rho(0,x)-\nabla_{x}\log \left(\det(\nabla_{x}T(t,x))\right)\right]\,.$ (17) $$

This formulation requires computing the inverse of Jacobian or solving related linear systems, resulting in a cost of $\mathcal{O}(d^{3})$ for a single $s_{t}$. As a consequence, the total cost for computing one score trajectory $\{s_{t_{j}}\}_{j=0}^{N_{t}}$ is $\mathcal{O}(N_{t}\,d^{3})$.

As an alternative, we can compute the score $s_{t}$ through the normalizing flow $l_{t}$, i.e., differentiating $l_{t}=\log\rho(t,z_{t})$ w.r.t $z_{t}$. We assume that the width of the neural network equation 8 is $k=\mathcal{O}(d)$. The score function can be reformulated as

$$ $s_{t}=\nabla_{T}\log\rho(t,T(t,x))=\nabla_{x}T(t,x)^{-\top}\nabla_{x}\log\rho( t,T(t,x))=\left(\dfrac{\partial z_{t}}{\partial z_{0}}\right)^{-\top}\dfrac{ \partial l_{t}}{\partial z_{0}}\,,$ (18) $$

where $\dfrac{\partial z_{t}}{\partial z_{0}}$ and $\dfrac{\partial l_{t}}{\partial z_{0}}$ are obtained from auto-differentiations of deep neural network functions in equation 15. Using the chain rule, the total computational computational cost for computing one trajectory of the score function is still $\mathcal{O}(N_{t}\,d^{3})$. Both equation 17 and equation 18 require computing the inverse Jacobian of $T$ or solving related linear systems, resulting in a cubic cost in the spatial dimension. As mentioned by Chen et al. () (section 4), this is the bottleneck of these methods. We leave more detailed discussions of both methods equation 17 and equation 18, and two other methods, to Appendix D.

In contrast, the high-order normalizing flow equation 12 only requires discretizing the ODE dynamics equation 3a and equation 3c, with a total cost of $\mathcal{O}(N_{t}\,d^{2})$ for computing one score trajectory $\{s_{t_{j}}\}_{j=0}^{N_{t}}$. As is shown in table 1, our formulation is more efficient compared with equation 17 and equation 18, especially in high dimensions. Additionally, we are able to compute the second-order score function through 13.

**Table 1: Complexity for different algorithms to compute a score trajectory.**
| pushforward map equation 17 | normalizing flow equation 18 | high-order normalizing flow equation 12 (ours) |
| --- | --- | --- |
| $\mathcal{O}(N_{t}\,d^{3})$ | $\mathcal{O}(N_{t}\,d^{3})$ | $\mathcal{O}(N_{t}\,d^{2})$ |

Example: Gaussian distribution.
Our formulation for the second-order score function also has advantages in the example of Gaussian distributions, where $H_{t}=\Sigma(t)^{-1}$ is exactly the inverse of the covariance matrix. Traditionally, estimating $\Sigma(t)^{-1}$ involves several steps. One needs to sample sufficiently many points $z_{0}^{(n)}$, simulate the state dynamic to obtain $z_{t}^{(n)}$, estimate the covariance matrix from these samples, and finally compute the inverse of this estimation. To achieve an error of $\varepsilon$, at least $\mathcal{O}(\varepsilon^{-2})$ samples are required for accurate covariance estimation. Compared to our second-order score dynamic equation 3d, we only need to pick $\Delta t=\mathcal{O}(\varepsilon)$ and compute $H_{t}$ through equation 13, with a total cost of $\mathcal{O}(\varepsilon^{-1})$. This results in a significantly more efficient algorithm.

## 4 Solving second-order MFC problems

In this section, we apply the score-based normalizing flow to solve the second-order MFC problem. In particular, we demonstrate that the MFC problem can be represented by the optimal control problem of the neural dynamical system in Proposition 1. For readers interested in the connection between MFC problems and generative AI, we refer to the work of Zhang & Katsoulakis (). In the context of machine learning, Neklyudov et al. () also studies similar problems, in which they compute Wasserstein Lagrangian flows. In this work, we can deal with the diffusion term in MFC problems using proposed neural normalization flows.

### 4.1 Formulation of MFC problem

We consider the following MFC problem

$$ $\inf_{v}\int_{0}^{t_{\text{end}}}\int_{\mathbb{R}^{d}}\left[L(t,x,v(t,x))+F(t, x,\rho(t,x))\right]\rho(t,x)\,\mathrm{d}x\,\mathrm{d}t+\int_{\mathbb{R}^{d}}G( x,\rho(t_{\text{end}},x))\rho(t_{\text{end}},x)\,\mathrm{d}x\,,$ (19) $$

where the density $\rho(t,x)$ satisfies the FP equation with a given initialization

$$ $\partial_{t}\rho(t,x)+\nabla_{x}\cdot(\rho(t,x)v(t,x))=\gamma\Delta_{x}\rho(t, x)\,,\quad\rho(0,x)=\rho_{0}(x)\,.$ (20) $$

The corresponding stochastic dynamic is

$$ $\mathrm{d}X_{t}=v(t,X_{t})\,\mathrm{d}t+\sqrt{2\gamma}\,\mathrm{d}W_{t}\,, \quad\quad X_{0}\sim\rho_{0}\,.$ (21) $$

Here $L:\mathbb{R}\times\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$ is the running cost function, and we assume that it is strongly convex in $v$. $F:\mathbb{R}\times\mathbb{R}^{d}\times\mathbb{R}\to\mathbb{R}$ is the cost that involves the density, which distinguishes the MFC problem from the optimal control. $G:\mathbb{R}^{d}\times\mathbb{R}\to\mathbb{R}$ is the terminal cost, which may also involve the density.

Given the density function $\rho(\cdot,\cdot)$, we define the composed velocity field $f$ as

$$ $f(t,x)=v(t,x)-\gamma\nabla_{x}\log\rho(t,x)\,.$ (22) $$

Let $z_{0}=x_{0}$, under the vector field $f$, we obtain the probability flow (Chen et al., ) of the stochastic dynamic $x_{t}$, given by $\partial_{t}z_{t}=f(t,z_{t})$, which coincides with equation 3a. This deterministic dynamic characterizes the probability distribution $\rho(t,x)$. If $z_{0}\sim\rho_{0}$, then the probability distribution for $z_{t}$ is exactly $\rho(t,\cdot)$. With this transformation, the velocity $v(t,x)$ becomes $f(t,x)+\gamma\nabla_{x}\log\rho(t,x)$, and the MFC problem follows

$$ $\displaystyle\inf_{f}$ $\displaystyle\int_{0}^{t_{\text{end}}}\int_{\mathbb{R}^{d}}\left[L\left(t,x,f( t,x)+\gamma\nabla_{x}\log\rho(t,x)\right)+F(t,x,\rho(t,x))\right]\rho(t,x)\, \mathrm{d}x\,\mathrm{d}t$ (23) $\displaystyle+\int_{\mathbb{R}^{d}}G(x,\rho(t_{\text{end}},x))\rho(t_{\text{ end}},x)\,\mathrm{d}x\,,$ $$

subject to the transport equation

$$ $\partial_{t}\rho(t,x)+\nabla_{x}\cdot(\rho(t,x)f(t,x))=0\,,\quad\rho(0,x)=\rho _{0}(x)\,.$ (24) $$

### 4.2 Modified HJB equation for MFC

In this section, we present a system of two PDEs that characterize the optimal solution for the MFC problem. Different from the traditional FP-HJB pair (see (Bensoussan et al., , Chapter 4)), our system consists of a transport equation obtained from equation 24, and a modified HJB equation, tailored for the modified MFC problem equation 23. This characterization could serve as a regularizer to enhance the loss function in numerical algorithms. We define the Hamiltonian $H:\mathbb{R}\times\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$ by

$$ $H(t,x,p)=\sup_{v\in\mathbb{R}^{d}}\left(v^{\top}p-L(t,x,v)\right)\,.$ (25) $$

This definition aligns with classical control theory and is closely related to the maximum principle (Zhou & Lu, ). The solution of the MFC problem is summarized by the following proposition.

###### Proposition 3 .

Let $L$ be strongly convex in $v$, then the solution to the MFC problem equation 23 is as follows. Consider a function $\psi\colon[0,t_{\textrm{end}}]\times\mathbb{R}^{d}\rightarrow\mathbb{R}$, such that

$$ $f(t,x)=\mathrm{D}_{p}H(t,x,\nabla_{x}\psi(t,x)+\gamma\nabla_{x}\log\rho(t,x))- \gamma\nabla_{x}\log\rho(t,x)\,,$ (26) $$

where the density function $\rho(t,x)$ and $\psi\colon[0,t_{\text{end}}]\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ satisfy the following system of equations

$$ $\left\{\begin{aligned} &\partial_{t}\rho(t,x)+\nabla_{x}\cdot\left(\rho(t,x) \mathrm{D}_{p}H\right)=\gamma\Delta_{x}\rho(t,x)\,,\\ &\partial_{t}\psi(t,x)+\nabla_{x}\psi(t,x)^{\top}\mathrm{D}_{p}H-\gamma\nabla_ {x}\cdot\mathrm{D}_{p}H+\gamma\Delta_{x}\psi(t,x)-L(t,x,\mathrm{D}_{p}H)\\ &\quad-\widetilde{F}(t,x,\rho(t,x))+2\gamma^{2}\Delta_{x}\log\rho(t,x)+\gamma^ {2}\left|\nabla_{x}\log\rho(t,x)\right|^{2}=0\,,\\ &\rho(0,x)=\rho_{0}(x)\,,\quad\psi(t_{\text{end}},x)=-\widetilde{G}(x,\rho(t_{ \text{end}},x))-\gamma\log\rho(t_{\text{end}},x)\,,\end{aligned}\right.$ (27) $$

Here, $\mathrm{D}_{p}H$ is short for $\mathrm{D}_{p}H(t,x,\nabla_{x}\psi(t,x)+\gamma\nabla_{x}\log\rho(t,x))$, $\widetilde{F}(t,x,\rho)=\dfrac{\partial}{\partial\rho}(F(t,x,\rho)\rho)=\dfrac
{\partial F}{\partial\rho}(t,x,\rho)\rho+F(t,x,\rho)$, and $\widetilde{G}(x,\rho)=\dfrac{\partial G}{\partial\rho}(x,\rho)\rho+G(x,\rho)$.

In the LQ problem, we let $L(t,x,v)=\frac{1}{2}|v|^{2}$ and $F(t,x,\rho)=0$. The MFC problem becomes

$$ $\inf_{f}\int_{0}^{t_{\text{end}}}\int_{\mathbb{R}^{d}}\frac{1}{2}\left|f(t,x)+ \gamma\nabla_{x}\log\rho(t,x)\right|^{2}\rho(t,x)\,\mathrm{d}x\,\mathrm{d}t+ \int_{\mathbb{R}^{d}}G(x,\rho(t_{\text{end}},x))\rho(t_{\text{end}},x)\, \mathrm{d}x\,,$ (28) $$

subject to equation 24. In this case, the result becomes the Corollary below.

###### Corollary 1 .

In the LQ problem where $L(t,x,v)=\frac{1}{2}|v|^{2}$ and $F(t,x,\rho)=0$, the solution of the MFC problem equation 28 is given by

$$ $f(t,x)=\nabla_{x}\psi(t,x)\,,$ $$

and

$$ $\left\{\begin{aligned} &\partial_{t}\rho(t,x)+\nabla_{x}\cdot\left(\rho(t,x) \nabla_{x}\psi(t,x)\right)=0\,,\\ &\partial_{t}\psi(t,x)+\frac{1}{2}\left|\nabla_{x}\psi(t,x)\right|^{2}+\gamma^ {2}\Delta_{x}\log\rho(t,x)+\frac{1}{2}\gamma^{2}\left|\nabla_{x}\log\rho(t,x) \right|^{2}=0\,,\\ &\rho(0,x)=\rho_{0}(x)\,,\quad\psi(t_{\text{end}},x)=-\widetilde{G}(x,\rho(t_{ \text{end}},x))-\gamma\log\rho(t_{\text{end}},x)\,.\end{aligned}\right.$ (29) $$

If we further define the Fisher information as $I[\rho]:=\int_{\mathbb{R}^{d}}\left|\nabla_{x}\log\rho(x)\right|^{2}\rho(x)\,
\mathrm{d}x$, then the equation for $\psi$ in equation equation 29 becomes

$$ $\partial_{t}\psi(t,x)+\frac{1}{2}\left|\nabla_{x}\psi(t,x)\right|^{2}-\frac{1} {2}\gamma^{2}\dfrac{\delta I[\rho(t,\cdot)]}{\delta\rho(t,\cdot)}(x)=0\,.$ $$

Let $z_{t}$ be the characteristic line of the state trajectory, satisfying $\partial_{t}z_{t}=f(t,z_{t})$ where $f$ is given by equation 26. Then, a residual of the HJB equation along the trajectory $z_{t}$ can be computed to enhance the objective function. For example, in the LQ example in Corollary 1, we know that $\psi$ satisfies

$$ $\partial_{t}\psi(t,z_{t})+\frac{1}{2}\left|\nabla_{z}\psi(t,z_{t})\right|^{2}+ \gamma^{2}\left(\operatorname{Tr}(H_{t})+\frac{1}{2}|s_{t}|^{2}\right)=0\,.$ (30) $$

Also, by equation 22, we can recover the optimal velocity field $v$ through

$$ $v(t,z_{t})=f(t,z_{t})+\gamma s_{t}\,,$ $$

where $f$ is the optimal vector field of the modified MFC problem. This regularization technique could be extended to a more general setting, such as the flow matching problem for overdamped Langevin dynamics. See Corollary 2 in Appendix A for details.

###### Proposition 3 .

###### Corollary 1 .

### 4.3 Numerical algorithms

In this section, we present numerical algorithms to solve the MFC problem. We minimize the objective equation 23 to obtain the optimal vector field $f$. Additionally, we can incorporate the residual of the HJB equation as a regularizer to enhance the loss functional, as discussed after Corollary 1.

The composed velocity field $f$ is parametrized as a neural network, as defined in equation 8. In order to simulate the cost functional numerically, we sample multiple initial points $z_{0}^{(n)}\sim\rho_{0}$ with a batch size $N_{z}$. Then, we can simulate the dynamics of $z_{t}^{(n)}$, $l_{t}^{(n)}$, $\tilde{l}_{t}^{(n)}$, $s_{t}^{(n)}$, and $H_{t}^{(n)}$ numerically for each particle through equation 9, equation 10, equation 11, equation 12, and equation 13 respectively, where the derivatives of $f$ are obtained via auto-differentiation. Then, we simulate the cost functional equation 23 by

$$ $\mathcal{L}_{\text{cost}}=\dfrac{1}{N_{z}}\sum_{n=1}^{N_{z}}\left[\sum_{j=0}^{ N_{t}-1}\left(L\left(t_{j},z_{t_{j}}^{(n)},f(t_{j},z_{t_{j}}^{(n)};\theta)- \gamma s_{t_{j}}^{(n)}\right)+F\left(t_{j},z_{t_{j}}^{(n)},\tilde{l}_{t_{j}}^{ (n)}\right)\right)\Delta t+G(z_{t_{N_{t}}}^{(n)},\tilde{l}_{t_{N_{t}}}^{(n)}) \right]\,.$ (31) $$

With this simulation cost, we minimize this loss $\mathcal{L}_{\text{cost}}$ using the Adam method. The algorithm is summarized in Algorithm 1.

Figure: Algorithm 1 Score-based normalizing flow solver for the MFC problem

To further improve performance, we can add a penalty of the residual of the HJB equation, denoted by $\mathcal{L}_{\text{HJB}}$, to enhance the loss function. This regularization technique is detailed in Appendix B.1. The total loss function is hence given by

$$ $\mathcal{L}_{\text{total}}=\mathcal{L}_{\text{cost}}+\lambda\mathcal{L}_{\text {HJB}}\,,$ (32) $$

where $\lambda\geq 0$ is a weight parameter.

## 5 Numerical results

We present numerical results for several examples in this section, including the regularized Wasserstein proximal operator (RWPO), the flow matching problem, a linear-quadratic (LQ) problem with an entropy potential cost, and an example with double well potential. For examples with exact solutions, we present all errors for the density ($\text{err}_{\rho}$), velocity field ($\text{err}_{f}$), and score function ($\text{err}_{s}$). These errors are calculated as averages over $10$ independent runs, see equation 95 in Appendix B.3 for details. The parameters used for all numerical examples are provided in Appendix B.4. Part of the results are deferred to Appendix C.

### 5.1 Regularized Wasserstein proximal operator

In the RWPO problem, the objective is a regularized Benamou-Brenier formulation for optimal transportation (Benamou & Brenier, ). We minimize the cost functional

$$ $\inf_{v}\int_{0}^{1}\int_{\mathbb{R}^{d}}\frac{1}{2}\left|v(t,x)\right|^{2} \rho(t,x)\,\mathrm{d}x\,\mathrm{d}t+\int_{\mathbb{R}^{d}}G(x)\rho(1,x)\, \mathrm{d}x\,,$ $$

subject to the FP equation $\partial_{t}\rho(t,x)+\nabla_{x}\cdot(\rho(t,x)v(t,x))=\gamma\Delta_{x}\rho(t,x)$. Here we set $G(x)=|x|^{2}/2$ and $\rho_{0}(x)=(8\pi)^{-d/2}\exp\left(-|x|^{2}/8\right)$. We test Algorithm 1 on this problem in $1$, $2$, and $10$ dimensions. These errors are summarized in Table 2. Detailed definitions for $\text{err}_{\rho}$, $\text{err}_{f}$, and $\text{err}_{s}$ are given in Appendix B.3. Additionally, we report the cost gap, which represents the difference between the computed cost and the optimal cost, averaged over 10 independent runs.

**Table 2: Errors for the RWPO problem.**
| errors | $\text{err}_{\rho}$ | $\text{err}_{f}$ | $\text{err}_{s}$ | cost gap |
| --- | --- | --- | --- | --- |
| $1d$ | $3.52\text{\times}{10}^{-3}$ | $3.45\text{\times}{10}^{-2}$ | $8.52\text{\times}{10}^{-3}$ | $1.33\text{\times}{10}^{-2}$ |
| $2d$ | $6.53\text{\times}{10}^{-3}$ | $3.40\text{\times}{10}^{-2}$ | $4.49\text{\times}{10}^{-2}$ | $1.24\text{\times}{10}^{-2}$ |
| $10d$ | $1.68\text{\times}{10}^{-3}$ | $3.15\text{\times}{10}^{-2}$ | $6.42\text{\times}{10}^{-2}$ | $1.69\text{\times}{10}^{-1}$ |

The results are also visualized in Figure 1. The plot on the left shows the cost functional through training in $1$ dimension, which becomes close to the optimal cost in red. The plot in the middle compares the evolution of the density function computed through equation 11 under trained velocity with the true density evolution. Our density dynamic equation 4 accurately captures the density evolution. The plot on the right shows the particle trajectories of $z_{t}$ in $2$ dimensions and compares them with the stochastic dynamics. Our probability flow ODE demonstrates a structured behavior.
We leave further numerical results of this example in Appendix C.1.

Figure: Figure 1: Numerical results for the RWPO problem. Left: cost functional through training with optimal cost in $1d$. Middle: density evolution through equation 4 and comparison with true density in $1d$. Right: trajectories of $z_{t}$ and comparison with the stochastic dynamic. The score dynamic demonstrates a structured behavior compared with the corresponding stochastic trajectory simulated by stochastic differential equations.
Refer to caption: x1.png

Adding HJB regularizer. We also test the regularized Algorithm 2 for this example. We compare all results with weight parameters $\lambda=$1\text{\times}{10}^{-3}$$ and $\lambda=0$ in equation 32 to study the effect of this regularization. All results are presented in table 3. The errors for the regularized algorithm are significantly smaller.

**Table 3: Errors for regularized MFC solver of RWPO problem. The errors for the regularized algorithm are significantly smaller.**
| regularization | $\lambda=$1\text{\times}{10}^{-3}$$ | $\lambda=0$ |  |  |
| --- | --- | --- | --- | --- |
| errors | $\text{err}_{A}$ | $\text{err}_{B}$ | $\text{err}_{A}$ | $\text{err}_{B}$ |
| $1d$ | $2.33\text{\times}{10}^{-3}$ | $3.03\text{\times}{10}^{-3}$ | $1.36\text{\times}{10}^{-2}$ | $1.00\text{\times}{10}^{-2}$ |
| $2d$ | $3.30\text{\times}{10}^{-3}$ | $4.10\text{\times}{10}^{-3}$ | $1.74\text{\times}{10}^{-2}$ | $1.70\text{\times}{10}^{-2}$ |
| $10d$ | $7.52\text{\times}{10}^{-3}$ | $4.45\text{\times}{10}^{-3}$ | $8.34\text{\times}{10}^{-2}$ | $4.02\text{\times}{10}^{-2}$ |

### 5.2 Flow matching for solving FP equations

Flow matching has emerged as an important problem, which is closely related to generative models (Lipman et al., ; Boffi & Vanden-Eijnden, ). The problem is to simulate the probability density function of a stochastic dynamics:

$$ $\mathrm{d}X_{t}=b(t,X_{t})\,\mathrm{d}t+\sqrt{2\gamma}\,\mathrm{d}W_{t},\quad X _{0}\sim\rho_{0}$ (33) $$

where $b\colon\mathbb{R}^{+}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is a known drift vector field, $\rho_{0}\colon\mathbb{R}^{d}\rightarrow\mathbb{R}$ is an initial value probability density function. In this example, we also assume the drift function satisfies $b(t,x)=-\nabla V(x)$, which is the negative gradient of some potential function $V(x)\in C_{\text{loc}}^{4}(\mathbb{R}^{d})$. In this case, SDE equation 33 is the overdamped Langevin dynamic.
To simulate the density function of stochastic process $X_{t}$ in equation 33, we design the following MFC problem. We minimize the objective functional (loss function) as

$$ $\inf_{v}\int_{0}^{t_{\text{end}}}\int_{\mathbb{R}^{d}}\frac{1}{2}\left|v(t,x)- b(t,x)\right|^{2}\rho(t,x)\,\mathrm{d}x\,\mathrm{d}t\,,$ $$

subject to the FP equation equation 20. After the score transformation equation 22, the above second order MFC problem becomes

$$ $\inf_{f}\int_{0}^{t_{\text{end}}}\int_{\mathbb{R}^{d}}\frac{1}{2}\left|f(t,x)- b(t,x)+\gamma\nabla_{x}\log\rho(t,x)\right|^{2}\rho(t,x)\,\mathrm{d}x\,\mathrm {d}t\,,$ $$

subject to the transport equation equation 24.
The state dynamic hence becomes $\partial_{t}z_{t}=-\nabla_{z}V(z_{t})-\gamma s_{t}$. This transformation equation 22 seems to complicate the problem. However, we are able to learn the whole FP equation of the overdamped Langevin dynamic with a given initial distribution $\rho_{0}$ or samples $\{z_{0}^{(n)}\}_{n=1}^{n_{z}}$. More importantly, we can record all intermediate time steps and compute the density function during the entire time domain $[0,t_{\text{end}}]$.

We first present the flow matching problem for an Ornstein–-Uhlenbeck (OU) process, where an explicit solution is available, given in Appendix A.3. The algorithm is described with details in Appendix B.1 and summarized in Algorithm 3. Similar to the previous section, we compare the numerical results with and without regularizations, presented in Table 4. The errors for the regularized algorithm are significantly smaller. We also provide visualized results in Figure 3 and 4 in Appendix C.2, which demonstrates the accuracy of our algorithm.

**Table 4: Errors for the regularized MFC solver of flow matching for OU processes. The errors for the regularized algorithm are significantly smaller.**
| regularization | $\lambda=$1\text{\times}{10}^{-3}$$ | $\lambda=0$ |  |  |
| --- | --- | --- | --- | --- |
| errors | $\text{err}_{A}$ | $\text{err}_{B}$ | $\text{err}_{A}$ | $\text{err}_{B}$ |
| $1d$ | $9.81\text{\times}{10}^{-4}$ | $8.58\text{\times}{10}^{-4}$ | $7.89\text{\times}{10}^{-3}$ | $7.47\text{\times}{10}^{-3}$ |
| $2d$ | $1.50\text{\times}{10}^{-3}$ | $1.37\text{\times}{10}^{-3}$ | $2.02\text{\times}{10}^{-2}$ | $1.40\text{\times}{10}^{-2}$ |
| $10d$ | $4.08\text{\times}{10}^{-3}$ | $6.73\text{\times}{10}^{-3}$ | $7.86\text{\times}{10}^{-2}$ | $3.89\text{\times}{10}^{-2}$ |

In addition, we design an algorithm that partitions a long time interval into several sub-intervals, which could potentially resolve the issue of long time horizon. The model is trained consecutively over each sub-interval, referred to as the multi-stage splicing method. We present an example in $2$ dimensions to demonstrate effectiveness of this splicing method, in which the invariant distribution of the SDE equation 33 is in a double moon shape (cf. equation 99). Our splicing method (Algorithm 4) is able to capture the state dynamic within each interval, and pass the information to the next interval. The detailed numerical implementation and numerical results are presented in Appendix C.2.

We also present two more examples in Appendix C.3 and C.4, where we apply the proposed neural ODE system to approximate the LQ MFC problems.

## 6 Conclusion

In this paper, we propose a neural ODE system to compute evolutions of first- and second-order score functions along trajectories. The forward Euler discretization of neural ODE system satisfies a system of normalization flows along a deep neural network. We then apply the proposed neural ODE system to solve second-order MFC problems. The effectiveness and accuracy of our method are validated through numerical examples of RWPO problems, score-based flow matching problems for FP equations, LQ problems, and the double well potential problems, providing reassurance of its reliability and effectiveness.

In future work, we shall conduct a thorough numerical analysis of the optimal control problem within the neural ODE system. Specifically, for second-order MFC problems, error analysis for the score estimations is necessary. Although we observe some numerical advantages in using the score function to approximate the solution of viscous HJB equation, the underlying error in the current numerical scheme remains unclear. Additionally, exploring neural ODE systems for inverse problems presents an exciting avenue for future research. This direction involves learning and approximating stochastic trajectories from data using the neural ODE system.

## Acknowledgments

M. Zhou is partially supported by AFOSR YIP award No.
FA9550-23-1-0087. M. Zhou and S. Osher are partially funded by AFOSR MURI
FA9550-18-502 and ONR N00014-20-1-2787. Wuchen Li is partially supported
by the AFOSR YIP award No. FA9550-23-1-0087, NSF DMS-2245097, and NSF RTG: 2038080.

## Appendix A Technical details for score-based normalizing flows and MFC problems

We present detailed proofs of all propositions about normalization flows and MFC problems in this section.

### A.1 Proofs for score-based normalizing flow

Continuity equation for density.
Let $T(t,x)$ be smooth and $T(t,\cdot)$ be invertible for all $t\geq 0$. Let $f(t,\cdot):\mathbb{R}^{d}\to\mathbb{R}^{d}$ be the vector field of the state dynamic $z_{t}$, i.e., $f(t,T(t,x))=\partial_{t}T(t,x)$ for all $t$ and $x$. Then, the probability density function satisfies the continuity equation

$$ $\partial_{t}\rho(t,x)+\nabla_{x}\cdot(\rho(t,x)f(t,x))=0\,.$ (34) $$

###### Proof.

Using the chain rule, we have

$$ $\partial_{t}\nabla_{x}T(t,x)=\nabla_{x}\partial_{t}T(t,x)=\nabla_{x}f(t,T(t,x) )=\nabla_{T}f(t,T(t,x))\,\nabla_{x}T(t,x)\,,$ $$

which implies

$$ $\operatorname{Tr}\left[\partial_{t}\nabla_{x}T(t,x)\,\nabla_{x}T(t,x)^{-1} \right]=\operatorname{Tr}\left[\nabla_{T}f(t,T(t,x))\right]=\nabla_{T}\cdot f( t,T(t,x))\,.$ (35) $$

If we compute the derivative of the MA equation equation 1 w.r.t. $t$, we obtain

$$ $\displaystyle 0$ $\displaystyle=\dfrac{\mathrm{d}}{\mathrm{d}t}\left[\rho(t,T(t,x))\det\left( \nabla_{x}T(t,x)\right)\right]$ $\displaystyle=\dfrac{\mathrm{d}}{\mathrm{d}t}\rho(t,T(t,x))\,\det\left(\nabla_ {x}T(t,x)\right)+\rho(t,T(t,x))\,\dfrac{\mathrm{d}}{\mathrm{d}t}\det\left( \nabla_{x}T(t,x)\right)$ $\displaystyle=\left[\partial_{t}\rho(t,T(t,x))+\nabla_{T}\rho(t,T(t,x))\, \partial_{t}T(t,x)\right]\,\det\left(\nabla_{x}T(t,x)\right)$ $\displaystyle\quad+\rho(t,T(t,x))\operatorname{Tr}\left[\partial_{t}\nabla_{x} T(t,x)\,\nabla_{x}T(t,x)^{-1}\right]\,\det\left(\nabla_{x}T(t,x)\right)$ $\displaystyle=\left[\partial_{t}\rho(t,T(t,x))+\nabla_{T}\rho(t,T(t,x))^{\top} f(t,T(t,x))+\rho(t,T(t,x))\nabla_{T}\cdot f(t,T(t,x))\right]$ $\displaystyle\quad\cdot\det\left(\nabla_{x}T(t,x)\right)$ $\displaystyle=\left[\partial_{t}\rho(t,T(t,x))+\nabla_{T}\cdot\left(\rho(t,T(t ,x))f(t,T(t,x))\right)\right]\cdot\det\left(\nabla_{x}T(t,x)\right)\,,$ $$

where we have used equation 35 in the fourth equality. Since $T(t,\cdot)$ is invertible, $\det\left(\nabla_{x}T(t,x)\right)\neq 0$, which implies

$$ $\partial_{t}\rho(t,T(t,x))+\nabla_{T}\cdot\left(\rho(t,T(t,x))f(t,T(t,x)) \right)=0\,,$ $$

for all $t$ and $x$. Therefore,

$$ $\partial_{t}\rho(t,x)+\nabla_{x}\cdot(\rho(t,x)f(t,x))=0\,.$ (36) $$

∎

###### Proposition 1 (Neural ODE system) .

Functions $z_{t}$, $l_{t}$, $s_{t}$, and $H_{t}$ satisfy the following ODE dynamics.

$$ $\displaystyle\partial_{t}z_{t}$ $\displaystyle=f(t,z_{t})\,,$ (37) $\displaystyle\partial_{t}l_{t}$ $\displaystyle=-\nabla_{z}\cdot f(t,z_{t})\,,$ (38) $\displaystyle\partial_{t}s_{t}$ $\displaystyle=-\nabla_{z}f(t,z_{t})^{\top}s_{t}-\nabla_{z}(\nabla_{z}\cdot f(t ,z_{t}))\,,$ (39) $\displaystyle\partial_{t}H_{t}$ $\displaystyle=-\sum_{i=1}^{d}s_{it}\nabla_{z}^{2}f_{i}(t,z_{t})-\nabla_{z}^{2} (\nabla_{z}\cdot f(t,z_{t}))-H_{t}\nabla_{z}f(t,z_{t})-\nabla_{z}f(t,z_{t})^{ \top}H_{t}\,,$ (40) $$

where $s_{it}$ and $f_{i}$ are the $i$-th component of $s_{t}$ and $f$ respectively. Also, $\tilde{l}_{t}=\rho(t,z_{t})$ satisfies

$$ $\partial_{t}\tilde{l}_{t}=-\nabla_{z}\cdot f(t,z_{t})\,\tilde{l}_{t}\,.$ $$

###### Proof.

The dynamic for $z_{t}$ equation 3a comes from the definition directly:

$$ $\partial_{t}z_{t}=\partial_{t}T(t,z_{0})=f(t,T(t,z_{0}))=f(t,z_{t})\,.$ $$

We next show equation 4 and equation 3b.

$$ $\displaystyle\quad\partial_{t}\tilde{l}_{t}=\dfrac{\mathrm{d}}{\mathrm{d}t} \rho(t,z_{t})=\partial_{t}\rho(t,z_{t})+\nabla_{z}\rho(t,z_{t})^{\top}\, \partial_{t}z_{t}=\partial_{t}\rho(t,z_{t})+\nabla_{z}\rho(t,z_{t})^{\top}\,f( t,z_{t})$ $\displaystyle=-\nabla_{z}\cdot[f(t,z_{t})\,\rho(t,z_{t})]+\nabla_{z}\rho(t,z_{ t})^{\top}\,f(t,z_{t})=-\nabla_{z}\cdot f(t,z_{t})\,\rho(t,z_{t})=-\nabla_{z} \cdot f(t,z_{t})\,\tilde{l}_{t}\,.$ $$

where we have used the continuity equation equation 2 in the fourth equality. Consequently,

$$ $\partial_{t}l_{t}=\dfrac{\mathrm{d}}{\mathrm{d}t}\rho(t,z_{t})/\rho(t,z_{t})=- \nabla_{z}\cdot f(t,z_{t})=-\nabla_{z}\cdot f(t,z_{t})\,.$ $$

We next show equation 3c. We will use the fact frequently that $\nabla_{z}\log\rho(t,z_{t})=\nabla_{z}\rho(t,z_{t})/\rho(t,z_{t})$. We first compute $\dfrac{\mathrm{d}}{\mathrm{d}t}(\nabla_{z}\rho(t,z_{t}))$

$$ $\displaystyle\quad\dfrac{\mathrm{d}}{\mathrm{d}t}\left[\nabla_{z}\rho(t,z_{t}) \right]=(\partial_{t}\nabla_{z})\rho(t,z_{t})+\nabla_{z}^{2}\rho(t,z_{t})\, \partial_{t}z_{t}=\nabla_{z}\partial_{t}\rho(t,z_{t})+\nabla_{z}^{2}\rho(t,z_{ t})\,f(t,z_{t})$ (41) $\displaystyle=-\nabla_{z}\left[\nabla_{z}\cdot\left(\rho(t,z_{t})\,f(t,z_{t}) \right)\right]+\nabla_{z}^{2}\rho(t,z_{t})\,f(t,z_{t})$ $\displaystyle=-\nabla_{z}\left[\nabla_{z}\rho(t,z_{t})^{\top}\,f(t,z_{t})+\rho (t,z_{t})\,\nabla_{z}\cdot f(t,z_{t})\right]+\nabla_{z}^{2}\rho(t,z_{t})\,f(t, z_{t})$ $\displaystyle=-\nabla_{z}f(t,z_{t})^{\top}\,\nabla_{z}\rho(t,z_{t})-(\nabla_{z }\cdot f(t,z_{t}))\,\nabla_{z}\rho(t,z_{t})-\rho(t,z_{t})\,\nabla_{z}\left( \nabla_{z}\cdot f(t,z_{t})\right)\,.$ $$

Therefore,

$$ $\displaystyle\quad\partial_{t}s_{t}=\dfrac{\mathrm{d}}{\mathrm{d}t}\left[ \nabla_{z}\log\rho(t,z_{t})\right]=\dfrac{\mathrm{d}}{\mathrm{d}t}\left[\nabla _{z}\rho(t,z_{t})/\rho(t,z_{t})\right]$ $\displaystyle=\dfrac{\mathrm{d}}{\mathrm{d}t}\left(\nabla_{z}\rho(t,z_{t}) \right)/\rho(t,z_{t})-\nabla_{z}\rho(t,z_{t})\,\dfrac{\mathrm{d}}{\mathrm{d}t} \left(\rho(t,z_{t})\right)/\rho(t,z_{t})^{2}$ $\displaystyle=-\nabla_{z}f(t,z_{t})^{\top}\,\nabla_{z}\log\rho(t,z_{t})-\nabla _{z}\cdot f(t,z_{t})\,\nabla_{z}\log\rho(t,z_{t})$ $\displaystyle\quad-\nabla_{z}\left(\nabla_{z}\cdot f(t,z_{t})\right)+\nabla_{z }\log\rho(t,z_{t})\,\nabla_{z}\cdot f(t,z_{t})$ $\displaystyle=-\nabla_{z}f(t,z_{t})^{\top}\,\nabla_{z}\log\rho(t,z_{t})-\nabla _{z}\left(\nabla_{z}\cdot f(t,z_{t})\right)$ $\displaystyle=-\nabla_{z}f(t,z_{t})^{\top}\,s_{t}-\nabla_{z}\left(\nabla_{z} \cdot f(t,z_{t})\right)\,,$ $$

where we have used equation 41 and equation 4 in the third inequality.

Finally we prove equation 3d. In order to simplify notation, we will omit $(t,z_{t})$ when there is no confusion and rewrite $\rho(t,z_{t})$ and $f(t,z_{t})$ as $\rho$ and $f$. We also denote $\partial_{i}$ as the partial derivative w.r.t. the $i$-th variable in $z$. We first compute $\dfrac{\mathrm{d}}{\mathrm{d}t}\nabla_{z}^{2}\rho(t,z_{t})$.

$$ $\displaystyle\quad\dfrac{\mathrm{d}}{\mathrm{d}t}\nabla_{z}^{2}\rho=\nabla_{z} ^{2}\partial_{t}\rho+\nabla_{z}^{3}\rho\cdot f=-\nabla_{z}^{2}[\nabla_{z}\cdot (\rho f)]+\nabla_{z}^{3}\rho\cdot f$ (42) $\displaystyle=-\nabla_{z}^{2}[\nabla_{z}\rho^{\top}f+\rho\,\nabla_{z}\cdot f]+ \nabla_{z}^{3}\rho\cdot f=-\nabla_{z}^{2}\left[\sum_{i=1}^{d}\partial_{i}\rho \,f_{i}+\rho\,\nabla_{z}\cdot f\right]+\nabla_{z}^{3}\rho\cdot f$ $\displaystyle=-\sum_{i=1}^{d}\left(\nabla_{z}^{2}\partial_{i}\rho\,f_{i}+ \partial_{i}\rho\,\nabla_{z}^{2}f_{i}+\nabla_{z}\partial_{i}\rho\,\nabla_{z}f_ {i}^{\top}+\nabla_{z}f_{i}\,\nabla_{z}\partial_{i}\rho^{\top}\right)$ $\displaystyle\quad-\left(\nabla_{z}^{2}\rho\,(\nabla_{z}\cdot f)+\rho\,\nabla_ {z}^{2}(\nabla_{z}\cdot f)+\nabla_{z}\rho\,\nabla_{z}(\nabla_{z}\cdot f)^{\top }+\nabla_{z}(\nabla_{z}\cdot f)\,\nabla_{z}\rho^{\top}\right)+\nabla_{z}^{3} \rho\cdot f$ $\displaystyle=-\sum_{i=1}^{d}\partial_{i}\rho\,\nabla_{z}^{2}f_{i}-\nabla_{z}^ {2}\rho\,\nabla_{z}f-(\nabla_{z}f)^{\top}\nabla_{z}^{2}\rho-\nabla_{z}^{2}\rho \,(\nabla_{z}\cdot f)-\rho\nabla_{z}^{2}(\nabla_{z}\cdot f)$ $\displaystyle\quad-\nabla_{z}\rho\,\nabla_{z}(\nabla_{z}\cdot f)^{\top}-\nabla _{z}(\nabla_{z}\cdot f)\nabla_{z}\rho^{\top}\,.$ $$

We also have

$$ $\nabla_{z}^{2}\log(\rho)=\nabla_{z}\left(\dfrac{\nabla_{z}\rho}{\rho}\right)= \dfrac{\nabla_{z}^{2}\rho}{\rho}-\dfrac{\nabla_{z}\rho\,\nabla_{z}\rho^{\top}} {\rho^{2}}\,,$ (43) $$

which implies

$$ $\dfrac{\nabla_{z}^{2}\rho}{\rho}=H_{t}+s_{t}s_{t}^{\top}\,.$ (44) $$

Dividing equation 42 by $\rho$, we obtain

$$ $\displaystyle\quad\dfrac{\frac{\mathrm{d}}{\mathrm{d}t}\nabla_{z}^{2}\rho}{\rho}$ (45) $\displaystyle=-\sum_{i=1}^{d}\dfrac{\partial_{i}\rho}{\rho}\,\nabla_{z}^{2}f_{ i}-\dfrac{\nabla_{z}^{2}\rho}{\rho}\,\nabla_{z}f-(\nabla_{z}f)^{\top}\dfrac{ \nabla_{z}^{2}\rho}{\rho}-\dfrac{\nabla_{z}^{2}\rho}{\rho}\,(\nabla_{z}\cdot f)$ $\displaystyle\quad-\nabla_{z}^{2}(\nabla_{z}\cdot f)-\dfrac{\nabla_{z}\rho}{ \rho}\,\nabla_{z}(\nabla_{z}\cdot f)^{\top}-\nabla_{z}(\nabla_{z}\cdot f) \dfrac{\nabla_{z}\rho^{\top}}{\rho}$ $\displaystyle=-\sum_{i=1}^{d}s_{it}\nabla_{z}^{2}f_{i}-\left(H_{t}+s_{t}s_{t}^ {\top}\right)\nabla_{z}f-(\nabla_{z}f)^{\top}\left(H_{t}+s_{t}s_{t}^{\top}\right)$ $\displaystyle\quad-(\nabla_{z}\cdot f)\left(H_{t}+s_{t}s_{t}^{\top}\right)- \nabla_{z}^{2}(\nabla_{z}\cdot f)-s_{t}\nabla_{z}(\nabla_{z}\cdot f)^{\top}- \nabla_{z}(\nabla_{z}\cdot f)s_{t}^{\top}\,,$ $$

where we have used equation equation 44 in the second equality.
Finally, we get

$$ $\displaystyle\quad\partial_{t}H_{t}=\dfrac{\mathrm{d}}{\mathrm{d}t}\nabla_{z}^ {2}\log\rho=\dfrac{\mathrm{d}}{\mathrm{d}t}\left(\dfrac{\nabla_{z}^{2}\rho}{ \rho}-\dfrac{\nabla_{z}\rho\,\nabla_{z}\rho^{\top}}{\rho^{2}}\right)$ $\displaystyle=\dfrac{\frac{\mathrm{d}}{\mathrm{d}t}\nabla_{z}^{2}\rho}{\rho}- \dfrac{\nabla_{z}^{2}\rho\frac{\mathrm{d}}{\mathrm{d}t}\rho}{\rho^{2}}-\dfrac{ \left(\frac{\mathrm{d}}{\mathrm{d}t}\nabla_{z}\rho\right)\nabla_{z}\rho^{\top} +\nabla_{z}\rho\left(\frac{\mathrm{d}}{\mathrm{d}t}\nabla_{z}\rho\right)^{\top }}{\rho^{2}}+\dfrac{2\nabla_{z}\rho\,\nabla_{z}\rho^{\top}\,\frac{\mathrm{d}}{ \mathrm{d}t}\rho}{\rho^{3}}$ $\displaystyle=-\sum_{i=1}^{d}s_{it}\nabla_{z}^{2}f_{i}-\left(H_{t}+s_{t}s_{t}^ {\top}\right)\nabla_{z}f-(\nabla_{z}f)^{\top}\left(H_{t}+s_{t}s_{t}^{\top} \right)-(\nabla_{z}\cdot f)\left(H_{t}+s_{t}s_{t}^{\top}\right)$ $\displaystyle\quad-\nabla_{z}^{2}(\nabla_{z}\cdot f)-s_{t}\nabla_{z}(\nabla_{z }\cdot f)^{\top}-\nabla_{z}(\nabla_{z}\cdot f)s_{t}^{\top}+\left(H_{t}+s_{t}s_ {t}^{\top}\right)(\nabla_{z}\cdot f)$ $\displaystyle\quad+\left(\nabla_{z}f^{\top}s_{t}+(\nabla_{z}\cdot f)\,s_{t}+ \nabla_{z}(\nabla_{z}\cdot f)\right)s_{t}^{\top}$ $\displaystyle\quad+s_{t}\left(\nabla_{z}f^{\top}s_{t}+(\nabla_{z}\cdot f)\,s_{ t}+\nabla_{z}(\nabla_{z}\cdot f)\right)^{\top}-2s_{t}s_{t}^{\top}(\nabla_{z} \cdot f)$ $\displaystyle=-\sum_{i=1}^{d}s_{it}\nabla_{z}^{2}f_{i}-\nabla_{z}^{2}(\nabla_{ z}\cdot f)-H_{t}\nabla_{z}f-(\nabla_{z}f)^{\top}H_{t}\,.$ $$

We have used equation 43 in the second equality and used equation 45, equation 44, equation 4, and equation 41 in the fourth equality. This finishes the proof.
∎

###### Proposition 2 (Information equality) .

The following equality holds:

$$ $\mathbb{E}\left[\operatorname{Tr}(H_{t})\right]=-\mathbb{E}\left[|s_{t}|^{2} \right]\,,\quad\textrm{for all $t\geq 0$\,.}$ $$

###### Proof.

For all $t\geq 0$, we have

$$ $\displaystyle\quad\mathbb{E}\left[\operatorname{Tr}(H_{t})+\left|s_{t}\right|^ {2}\right]=\mathbb{E}\left[\nabla_{z}\cdot\left(\nabla_{z}\log\rho(t,z_{t}) \right)+\left|\nabla_{z}\log\rho(t,z_{t})\right|^{2}\right]$ (46) $\displaystyle=\int_{\mathbb{R}^{d}}\left(\nabla_{x}\cdot\left(\nabla_{x}\log \rho(t,x)\right)+\left|\nabla_{x}\log\rho(t,x)\right|^{2}\right)\rho(t,x)\, \mathrm{d}x$ $\displaystyle=\int_{\mathbb{R}^{d}}\left[\nabla_{x}\cdot\left(\dfrac{\nabla_{x }\rho(t,x)}{\rho(t,x)}\right)\rho(t,x)+\dfrac{\left|\nabla_{x}\rho(t,x)\right| ^{2}}{\rho(t,x)^{2}}\rho(t,x)\right]\,\mathrm{d}x$ $\displaystyle=\int_{\mathbb{R}^{d}}\left[-\dfrac{\nabla_{x}\rho(t,x)^{\top}}{ \rho(t,x)}\nabla_{x}\rho(t,x)+\dfrac{\left|\nabla_{x}\rho(t,x)\right|^{2}}{ \rho(t,x)}\right]\,\mathrm{d}x=0\,,$ $$

where we have used the integration by part in the second last equality.
∎

Covariance evolution of centered Gaussian distributions.
Let $z_{0}$ follow a centered Gaussian distribution $N(0,\Sigma(0))$ and $\partial_{t}z_{t}=A(t)z_{t}$. Then $z_{t}$ is also a centered Gaussian distribution $N(0,\Sigma(t))$, where the covariance $\Sigma(t)$ satisfies a matrix ODE:

$$ $\partial_{t}\Sigma(t)=A(t)\Sigma(t)+\Sigma(t)A(t)^{\top}\,.$ (47) $$

###### Proof.

We write down the density function $\rho(t,x)$ for $N(0,\Sigma(t))$

$$ $\rho(t,x)=(2\pi)^{-\frac{d}{2}}\det(\Sigma(t))^{-\frac{1}{2}}\exp\left(-\frac{ 1}{2}x^{\top}\Sigma(t)^{-1}x\right)\,.$ (48) $$

It is sufficient to check that $\rho(t,x)$ satisfies the transport equation

$$ $\partial_{t}\rho(t,x)+\nabla_{x}\cdot\left(A(t)x\,\rho(t,x)\right)=0\,,$ $$

if equation equation 47 holds. We first compute

$$ $\displaystyle\dfrac{\mathrm{d}}{\mathrm{d}t}\det(\Sigma(t))$ $\displaystyle=\det(\Sigma(t))\operatorname{Tr}\left(\Sigma(t)^{-1}\,\partial_{ t}\Sigma(t)\right)=\det(\Sigma(t))\operatorname{Tr}\left(\Sigma(t)^{-1}\,(A(t) \Sigma(t)+\Sigma(t)A(t)^{\top})\right)$ (49) $\displaystyle=\det(\Sigma(t))\operatorname{Tr}\left(\Sigma(t)^{-1}A(t)\Sigma(t )+A(t)^{\top}\right)=2\det(\Sigma(t))\operatorname{Tr}(A(t))\,,$ (50) $$

which implies

$$ $\dfrac{\mathrm{d}}{\mathrm{d}t}\left(\det(\Sigma(t))^{-\frac{1}{2}}\right)=- \frac{1}{2}\det(\Sigma(t))^{-\frac{3}{2}}\,\dfrac{\mathrm{d}}{\mathrm{d}t}\det (\Sigma(t))=-\det(\Sigma(t))^{-\frac{1}{2}}\,\operatorname{Tr}(A(t))\,.$ (51) $$

We also have

$$ $\dfrac{\mathrm{d}}{\mathrm{d}t}\left(\Sigma(t)^{-1}\right)=-\Sigma(t)^{-1}\, \partial_{t}\Sigma(t)\,\Sigma(t)^{-1}=-\Sigma(t)^{-1}A(t)-A(t)^{\top}\Sigma(t) ^{-1}\,,$ (52) $$

which implies

$$ $\displaystyle\quad\dfrac{\mathrm{d}}{\mathrm{d}t}\exp\left(-\frac{1}{2}x^{\top }\Sigma(t)^{-1}x\right)=-\frac{1}{2}x^{\top}\dfrac{\mathrm{d}}{\mathrm{d}t} \left(\Sigma(t)^{-1}\right)\,x\,\exp\left(-\frac{1}{2}x^{\top}\Sigma(t)^{-1}x\right)$ (53) $\displaystyle=\frac{1}{2}x^{\top}\left(\Sigma(t)^{-1}A(t)+A(t)^{\top}\Sigma(t) ^{-1}\right)x\,\exp\left(-\frac{1}{2}x^{\top}\Sigma(t)^{-1}x\right)$ $\displaystyle=x^{\top}A(t)^{\top}\Sigma(t)^{-1}x\,\exp\left(-\frac{1}{2}x^{ \top}\Sigma(t)^{-1}x\right)\,.$ $$

Combining equation 51 and equation 53, we obtain

$$ $\displaystyle\partial_{t}\rho(t,x)$ $\displaystyle=(2\pi)^{-\frac{d}{2}}\,\dfrac{\mathrm{d}}{\mathrm{d}t}\left(\det (\Sigma(t))^{-\frac{1}{2}}\right)\exp\left(-\frac{1}{2}x^{\top}\Sigma(t)^{-1}x\right)$ (54) $\displaystyle\quad+(2\pi)^{-\frac{d}{2}}\det(\Sigma(t))^{-\frac{1}{2}}\,\dfrac {\mathrm{d}}{\mathrm{d}t}\exp\left(-\frac{1}{2}x^{\top}\Sigma(t)^{-1}x\right)$ $\displaystyle=\rho(t,x)\left[-\operatorname{Tr}(A(t))+x^{\top}A(t)^{\top} \Sigma(t)^{-1}x\right]\,.$ $$

Finally, we obtain

$$ $\displaystyle\quad\nabla_{x}\cdot\left(A(t)x\,\rho(t,x)\right)=\operatorname{ Tr}(A(t))\,\rho(t,x)+x^{\top}A(t)^{\top}\,\nabla_{x}\rho(t,x)$ (55) $\displaystyle=\operatorname{Tr}(A(t))\,\rho(t,x)+x^{\top}A(t)^{\top}\left(- \Sigma(t)^{-1}x\right)\rho(t,x)=-\partial_{t}\rho(t,x)\,,$ $$

which finishes the proof.
∎

###### Proof.

###### Proposition 1 (Neural ODE system) .

###### Proof.

###### Proposition 2 (Information equality) .

###### Proof.

###### Proof.

### A.2 Details for the MFC problems

We recall the formulation of the MFC problem

$$ $\inf_{v}\int_{0}^{t_{\text{end}}}\int_{\mathbb{R}^{d}}\left[L(t,x,v(t,x))+F(t, x,\rho(t,x))\right]\rho(t,x)\,\mathrm{d}x\,\mathrm{d}t+\int_{\mathbb{R}^{d}}G( x,\rho(t_{\text{end}},x))\rho(t_{\text{end}},x)\,\mathrm{d}x\,,$ (56) $$

where the density $\rho(t,x)$ satisfies the FP equation with a given initialization

$$ $\partial_{t}\rho(t,x)+\nabla_{x}\cdot(\rho(t,x)v(t,x))=\gamma\Delta_{x}\rho(t, x)\,,\quad\rho(0,x)=\rho_{0}(x)\,.$ (57) $$

This FP equation has a stochastic characterization

$$ $\mathrm{d}x_{t}=v(t,x_{t})\,\mathrm{d}t+\sqrt{2\gamma}\,\mathrm{d}W_{t}\,,$ (58) $$

where $W_{t}$ is a standard Brownian motion. For the well-posedness of the problem, we assume that the running cost $L$ is strongly convex in $v$, then the Hamiltonian, i.e. the Legendre transform (convex conjugate) of $L$, is well defined, given by

$$ $H(t,x,p)=\sup_{v\in\mathbb{R}^{d}}v^{\top}p-L(t,x,v)\,.$ (59) $$

According to standard results in convex analysis, $\mathrm{D}_{v}L$ and $\mathrm{D}_{p}H$ are the inverse functions to each other, in the sense that

$$ $\displaystyle\mathrm{D}_{p}H(t,x,\mathrm{D}_{v}L(t,x,v))$ $\displaystyle=v\,,\quad\forall v\,,$ (60) $\displaystyle\mathrm{D}_{v}L(t,x,\mathrm{D}_{p}H(t,x,p))$ $\displaystyle=p\,,\quad\forall p\,.$ $$

We recall that the composed velocity is

$$ $f(t,x)=v(t,x)-\gamma\nabla_{x}\log\rho(t,x)\,,$ (61) $$

and the probability flow $z_{t}$ is given by

$$ $\partial_{t}z_{t}=f(t,z_{t})\,.$ (62) $$

Note that the stochastic dynamic $x_{t}$ and the probability flow $z_{t}$ share the same probability distribution, while they are different dynamics. This $z_{t}$ serves as a characteristic line for the system. Along this line, a forward Euler scheme gives $\mathcal{O}(\Delta t)$ error, while a direct discretization of the stochastic dynamic has an error term $\mathcal{O}(\sqrt{\Delta t})$.

Next, we present the an analysis for the modified HJB equation of MFC problem, tailored for the probability flow dynamic.

###### Proposition 3 .

Let $L$ be strongly convex in $v$, then solution to the MFC problem equation 23 is as follows. Consider a function $\psi\colon[0,t_{\textrm{end}}]\times\mathbb{R}^{d}\rightarrow\mathbb{R}$, such that

$$ $f(t,x)=\mathrm{D}_{p}H(t,x,\nabla_{x}\psi(t,x)+\gamma\nabla_{x}\log\rho(t,x))- \gamma\nabla_{x}\log\rho(t,x),$ (63) $$

where the density function $\rho(t,x)$ and $\psi\colon[0,t_{\text{end}}]\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ satisfy the following system of equations

$$ $\left\{\begin{aligned} &\partial_{t}\rho(t,x)+\nabla_{x}\cdot\left(\rho(t,x) \mathrm{D}_{p}H\right)=\gamma\Delta_{x}\rho(t,x),\\ &\partial_{t}\psi(t,x)+\nabla_{x}\psi(t,x)^{\top}\mathrm{D}_{p}H-\gamma\nabla_ {x}\cdot\mathrm{D}_{p}H+\gamma\Delta_{x}\psi(t,x)-L(t,x,\mathrm{D}_{p}H)\\ &\quad-\widetilde{F}(t,x,\rho(t,x))+2\gamma^{2}\Delta_{x}\log\rho(t,x)+\gamma^ {2}\left|\nabla_{x}\log\rho(t,x)\right|^{2}=0,\\ &\rho(0,x)=\rho_{0}(x),\quad\psi(t_{\text{end}},x)=-\widetilde{G}(x,\rho(t_{ \text{end}},x))-\gamma\log\rho(t_{\text{end}},x),\end{aligned}\right.$ (64) $$

Here, $\mathrm{D}_{p}H$ is short for $\mathrm{D}_{p}H(t,x,\nabla_{x}\psi(t,x)+\gamma\nabla_{x}\log\rho(t,x))$, $\widetilde{F}(t,x,\rho)=\dfrac{\partial}{\partial\rho}(F(t,x,\rho)\rho)=\dfrac
{\partial F}{\partial\rho}(t,x,\rho)\rho+F(t,x,\rho)$, and $\widetilde{G}(x,\rho)=\dfrac{\partial G}{\partial\rho}(x,\rho)\rho+G(x,\rho)$.

###### Proof.

We start by adding a Lagrange multiplier $\phi(t,x)$ and obtain the augmented objective

$$ $\displaystyle\int_{0}^{t_{\text{end}}}\int_{\mathbb{R}^{d}}\left[\left(L(t,x,f (t,x)+\gamma\nabla_{x}\log\rho(t,x))+F(t,x,\rho(t,x))\right)\rho(t,x)\right.$ (65) $\displaystyle\left.\quad+\phi(t,x)\left(\partial_{t}\rho(t,x)+\nabla_{x}\cdot( f(t,x)\rho(t,x))\right)\right]\,\mathrm{d}x\,\mathrm{d}t+\int_{\mathbb{R}^{d}} G(x,\rho(t_{\text{end}},x))\rho(t_{\text{end}},x)\,\mathrm{d}x$ (66) $\displaystyle=$ $\displaystyle\int_{0}^{t_{\text{end}}}\int_{\mathbb{R}^{d}}\left[\left(L(t,x,f (t,x)+\gamma\nabla_{x}\log\rho(t,x))+F(t,x,\rho(t,x))\right)\rho(t,x)\right.$ (67) $\displaystyle\left.\quad-\partial_{t}\phi(t,x)\rho(t,x)-\nabla_{x}\phi(t,x)^{ \top}f(t,x)\rho(t,x)\right]\,\mathrm{d}x\,\mathrm{d}t$ (68) $\displaystyle\quad+\int_{\mathbb{R}^{d}}(G(x,\rho(t_{\text{end}},x))+\phi(t_{ \text{end}},x))\rho(t_{\text{end}},x)\,\mathrm{d}x-\int_{\mathbb{R}^{d}}\phi(0 ,x)\rho(0,x)\,\mathrm{d}x.$ (69) $$

Taking the variation of equation 68 w.r.t. $\rho(t_{\text{end}},\cdot)$, we get $\phi(t_{\text{end}},x)=-\widetilde{G}(x,\rho(t_{\text{end}},x))$. Taking the variation w.r.t. $f$, we obtain

$$ $\mathrm{D}_{v}L(t,x,f(t,x)+\gamma\nabla_{x}\log\rho(t,x))=\nabla_{x}\phi(t,x)\,,$ $$

which implies

$$ $f(t,x)+\gamma\nabla_{x}\log\rho(t,x)=\mathrm{D}_{p}H(t,x,\nabla_{x}\phi(t,x))\,,$ (70) $$

due to equation 60. Next, we denote $\mathrm{D}_{p}H(t,x,\nabla_{x}\phi(t,x))$ by $\mathrm{D}_{p}H$ for short. We then observe that

$$ $\mathrm{D}_{v}L:=\mathrm{D}_{v}L(t,x,f(t,x)+\gamma\nabla_{x}\log\rho(t,x))= \mathrm{D}_{v}L(t,x,\mathrm{D}_{p}H)=\nabla_{x}\phi(t,x).$ $$

Taking the variation of equation 68 w.r.t. $\rho(\cdot,\cdot)$, we get

$$ $\displaystyle 0$ $\displaystyle=-\gamma\nabla_{x}\cdot\left[\mathrm{D}_{v}L\,\rho(t,x)\right]/ \rho(t,x)+L(t,x,f(t,x)+\gamma\nabla_{x}\log\rho(t,x))$ (71) $\displaystyle\quad+\widetilde{F}(t,x,\rho(t,x))-\partial_{t}\phi(t,x)-\nabla_{ x}\phi(t,x)^{\top}f(t,x)$ (72) $\displaystyle=-\gamma\nabla_{x}\cdot\mathrm{D}_{v}L-\gamma\mathrm{D}_{v}L^{ \top}\nabla_{x}\log\rho(t,x)+L(t,x,\mathrm{D}_{p}H)$ (73) $\displaystyle\quad+\widetilde{F}(t,x,\rho(t,x))-\partial_{t}\phi(t,x)-\nabla_{ x}\phi(t,x)^{\top}\left(\mathrm{D}_{p}H-\gamma\nabla_{x}\log\rho(t,x)\right)$ (74) $\displaystyle=-\gamma\Delta_{x}\phi(t,x)+L(t,x,\mathrm{D}_{p}H)+\widetilde{F}( t,x,\rho(t,x))-\partial_{t}\phi(t,x)-\nabla_{x}\phi(t,x)^{\top}\mathrm{D}_{p}H\,,$ (75) $$

where we have used $\mathrm{D}_{v}L=\nabla_{x}\phi(t,x)$ in the last inequality. Equation equation 75 is the classical HJB equation corresponding to the MFC problem, which coincide with equation (4.6) in Bensoussan et al. () after a sign flip, where our $\mathrm{D}_{p}H$ is their $\hat{v}$. Please note that, throughout our derivation, $\rho(t,x)$ is not a general density, but the density under the optimal control.

We next present the modified HJB equation. We define

$$ $\psi(t,x):=\phi(t,x)-\gamma\log\rho(t,x)\,.$ (76) $$

This definition seems to complicate the system, but eventually gives us an elegant characterization, especially in the case of LQ problem (see Corollary 1). We then observe that

$$ $\displaystyle\quad\partial_{t}\log\rho(t,x)=\dfrac{\partial_{t}\rho(t,x)}{\rho (t,x)}=-\dfrac{\nabla_{x}\cdot(\rho(t,x)f(t,x))}{\rho(t,x)}$ (77) $\displaystyle=-\nabla_{x}\cdot f(t,x)-\nabla_{x}\log\rho(t,x)^{\top}f(t,x)$ $\displaystyle=-\nabla_{x}\cdot\mathrm{D}_{p}H+\gamma\Delta_{x}\log\rho(t,x)- \nabla_{x}\log\rho(t,x)^{\top}\mathrm{D}_{p}H+\gamma\left|\nabla_{x}\log\rho(t ,x)\right|^{2},$ $$

where the last equality is due to equation 70. Therefore, the HJB equation equation 75 is transformed into

$$ $\displaystyle 0$ $\displaystyle=-\gamma\Delta_{x}\psi(t,x)-\gamma^{2}\Delta_{x}\log\rho(t,x)+L(t ,x,\mathrm{D}_{p}H)+\widetilde{F}(t,x,\rho(t,x))-\partial_{t}\psi(t,x)$ $\displaystyle\quad-\gamma\partial_{t}\log\rho(t,x)-(\nabla_{x}\psi(t,x)+\gamma \nabla_{x}\log\rho(t,x))^{\top}\mathrm{D}_{p}H$ $\displaystyle=-\gamma\Delta_{x}\psi(t,x)-\gamma^{2}\Delta_{x}\log\rho(t,x)+L(t ,x,\mathrm{D}_{p}H)+\widetilde{F}(t,x,\rho(t,x))-\partial_{t}\psi(t,x)$ $\displaystyle\quad+\gamma\nabla_{x}\cdot\mathrm{D}_{p}H-\gamma^{2}\Delta_{x} \log\rho(t,x)+\gamma\nabla_{x}\log\rho(t,x)^{\top}\mathrm{D}_{p}H-\gamma^{2} \left|\nabla_{x}\log\rho(t,x)\right|^{2}$ $\displaystyle\quad-(\nabla_{x}\psi(t,x)+\gamma\nabla_{x}\log\rho(t,x))^{\top} \mathrm{D}_{p}H$ $\displaystyle=-\partial_{t}\psi(t,x)-\nabla_{x}\psi(t,x)^{\top}\mathrm{D}_{p}H +\gamma\nabla_{x}\cdot\mathrm{D}_{p}H-\gamma\Delta_{x}\psi(t,x)+L(t,x,\mathrm{ D}_{p}H)$ $\displaystyle\quad+\widetilde{F}(t,x,\rho(t,x))-2\gamma^{2}\Delta_{x}\log\rho( t,x)-\gamma^{2}\left|\nabla_{x}\log\rho(t,x)\right|^{2}.$ $$

Here, $\mathrm{D}_{p}H$ is short for $\mathrm{D}_{p}H(t,x,\nabla_{x}\psi(t,x)+\gamma\nabla_{x}\log\rho(t,x))$. We plug in the shift equation 76 in the first equality; we plug in equation 77 in the second equality. Therefore, we have recovered Proposition 3.
∎

###### Corollary 1 .

In the LQ problem where $L(t,x,v)=\frac{1}{2}|v|^{2}$ and $F(t,x,\rho)=0$, the solution of the MFC problem equation 23 is characterized by

$$ $f(t,x)=\nabla_{x}\psi(t,x)\,,$ $$

and

$$ $\left\{\begin{aligned} &\partial_{t}\rho(t,x)+\nabla_{x}\cdot\left(\rho(t,x) \nabla_{x}\psi(t,x)\right)=0\,,\\ &\partial_{t}\psi(t,x)+\frac{1}{2}\left|\nabla_{x}\psi(t,x)\right|^{2}+\gamma^ {2}\Delta_{x}\log\rho(t,x)+\frac{1}{2}\gamma^{2}\left|\nabla_{x}\log\rho(t,x) \right|^{2}=0\,,\\ &\rho(0,x)=\rho_{0}(x)\,,\quad\psi(t_{\text{end}},x)=-\widetilde{G}(x,\rho(t_{ \text{end}},x))-\gamma\log\rho(t_{\text{end}},x)\,.\end{aligned}\right.$ (78) $$

If we further define the Fisher information as $I[\rho]:=\int_{\mathbb{R}^{d}}\left|\nabla_{x}\log\rho(x)\right|^{2}\rho(x)\,
\mathrm{d}x$, then the equation for $\psi$ becomes

$$ $\partial_{t}\psi(t,x)+\frac{1}{2}\left|\nabla_{x}\psi(t,x)\right|^{2}-\frac{1} {2}\gamma^{2}\dfrac{\delta I[\rho(t,\cdot)]}{\delta\rho(t,\cdot)}(x)=0\,.$ $$

###### Proof.

When $L(t,x,v)=\frac{1}{2}|v|^{2}$, we have $\mathrm{D}_{v}L(t,x,v)=v$ and $\mathrm{D}_{p}H(t,x,p)=p$. So

$$ $\mathrm{D}_{p}H(t,x,\nabla_{x}\psi(t,x)+\gamma\nabla_{x}\log\rho(t,x))=\nabla_ {x}\psi(t,x)+\gamma\nabla_{x}\log\rho(t,x)\,.$ $$

Hence, the equation for $\psi$ becomes

$$ $\displaystyle 0$ $\displaystyle=\partial_{t}\psi(t,x)+\nabla_{x}\psi(t,x)^{\top}\left(\nabla_{x} \psi(t,x)+\gamma\nabla_{x}\log\rho(t,x)\right)-\gamma\Delta_{x}\psi(t,x)$ $\displaystyle\quad-\gamma^{2}\Delta_{x}\log\rho(t,x)+\gamma\Delta_{x}\psi(t,x) -\frac{1}{2}\left|\nabla_{x}\psi(t,x)+\gamma\nabla_{x}\log\rho(t,x)\right|^{2}$ $\displaystyle\quad+2\gamma^{2}\Delta_{x}\log\rho(t,x)+\gamma^{2}\left|\nabla_{ x}\log\rho(t,x)\right|^{2}$ $\displaystyle=\partial_{t}\psi(t,x)+\frac{1}{2}\left|\nabla_{x}\psi(t,x)\right |^{2}+\gamma^{2}\Delta_{x}\log\rho(t,x)+\frac{1}{2}\gamma^{2}\left|\nabla_{x} \log\rho(t,x)\right|^{2}.$ $$

A direct computation gives

$$ $\dfrac{\delta I[\rho]}{\delta\rho}(x)=-2\Delta_{x}\log\rho(x)-\left|\nabla_{x} \log\rho(x)\right|^{2},$ $$

which implies

$$ $\partial_{t}\psi(t,x)+\frac{1}{2}\left|\nabla_{x}\psi(t,x)\right|^{2}-\frac{1} {2}\gamma^{2}\dfrac{\delta I[\rho(t,\cdot)]}{\delta\rho(t,\cdot)}(x)=0\,.$ $$

∎

###### Proposition 3 .

###### Proof.

###### Corollary 1 .

###### Proof.

### A.3 Formulation for flow matching of overdamped Langevin dynamic

In this section, we present the flow matching problem of overdamped Langevin dynamic

$$ $\mathrm{d}X_{t}=b(t,X_{t})\,\mathrm{d}t+\sqrt{2\gamma}\,\mathrm{d}W_{t}\,, \quad\quad X_{0}\sim\rho_{0}\,.$ (79) $$

Here we assume that the drift is the negative gradient of some potential function

$$ $b(t,x)=-\nabla_{x}V(x)\,,$ $$

so that the stationary distribution for the overdamped Langevin dynamic is proportional to $\exp(-V(x)/\gamma)$. The goal for flow matching is to learn a velocity field $v(x)$ that matches $b(t,x)$. One minimizes the objective functional

$$ $\inf_{v}\int_{0}^{t_{\text{end}}}\int_{\mathbb{R}^{d}}\frac{1}{2}\left|v(t,x)- b(t,x)\right|^{2}\rho(t,x)\,\mathrm{d}x\,\mathrm{d}t\,,$ (80) $$

subject to the FP equation

$$ $\partial_{t}\rho(t,x)+\nabla_{x}\cdot(\rho(t,x)v(t,x))=\gamma\Delta_{x}\rho(t, x)\,,\quad\rho(0,x)=\rho_{0}(x)\,.$ (81) $$

The following corollary gives the exact formulation for the solution.

###### Corollary 2 .

For flow matching problem of overdamped Langevin dynamic where $L(t,x,v)=\frac{1}{2}|v-b(t,x)|^{2}$, $F(t,x,\rho)=0$ and $G(x,\rho)=0$, the solution of the MFC problem equation 23 is characterized by

$$ $f(t,x)=\nabla_{x}\psi(t,x)+b(t,x)\,,$ $$

and

$$ $\left\{\begin{aligned} &\partial_{t}\rho(t,x)+\nabla_{x}\cdot\left(\rho(t,x)( \nabla_{x}\psi(t,x)+b(t,x))\right)=0\,,\\ &\partial_{t}\psi(t,x)+\frac{1}{2}\left|\nabla_{x}\psi(t,x)\right|^{2}+\nabla_ {x}\psi(t,x)^{\top}b(t,x)-\gamma\nabla_{x}\cdot b(t,x)\\ &\hskip 43.36243pt+\gamma^{2}\Delta_{x}\log\rho(t,x)+\frac{1}{2}\gamma^{2} \left|\nabla_{x}\log\rho(t,x)\right|^{2}=0\,,\\ &\rho(0,x)=\rho_{0}(x)\,,\quad\psi(t_{\text{end}},x)=-\gamma\log\rho(t_{\text{ end}},x)\,.\end{aligned}\right.$ (82) $$

###### Proof.

Through direct computation, the convex conjugate of $L$ is

$$ $H(t,x,p)=\frac{1}{2}|p|^{2}+p^{\top}b(t,x)\,,$ $$

and $\mathrm{D}_{p}H(t,x,p)=p+b(t,x)$. Plugging these expressions into Proposition 3, we recover the results for this corollary.
∎

###### Corollary 2 .

###### Proof.

## Appendix B Details for numerical implementation

### B.1 HJB regularizer for MFC

In addition to the numerical algorithm in section 4.3, we can add the residual of the HJB equation equation 27 into the loss function. First, let us consider the HJB regularizer in the LQ Gaussian case, as presented in Corollary 1. Here the Gaussian case means that $\rho_{0}$ is also given as a Gaussian distribution, in which the solution of density function in LQ MFC problem stays in a time dependent Gaussian distribution.

In the LQ scenario, parametrizing the vector field $f$ directly as a neural network complicates the system, as terms like $\partial_{t}\psi$ in equation 29 are hard to implement. A more feasible approach is to parametrize $\psi$ as a neural network, allowing the vector field $f$ to be computed as $f(t,x)=\nabla_{x}\psi(t,x)$.

This parametrization, however, requires two additional orders of auto-differentiations compared to the original algorithm. First, if $\psi$ is parametrized as a neural network as in equation 8, then $f(t,x)=\nabla_{x}\psi(t,x)$ must be obtained through auto-differentiation. Second, the term $\Delta_{x}\rho(t,x)$ in equation 29 requires computing the dynamic for $H_{t}$ through equation 13, which requires another order of derivative. This term $H_{t}$ is not required in the original Algorithm 1. These two additional differentiations could potentially complicate the optimization landscape.

Therefore, in order to make a clean and fair comparison in studying the effect of the HJB regularizer, we consider the specific parametrization for $\psi$, given by

$$ $\psi(t,x)=\frac{1}{2}x^{\top}A(t)x+B(t)^{\top}x+C(t)\,,$ (83) $$

where $A(t):[0,t_{\text{end}}]\to\mathbb{R}^{d\times d}$, $B(t):[0,t_{\text{end}}]\to\mathbb{R}^{d}$, and $C(t):[0,t_{\text{end}}]\to\mathbb{R}$. These three functions can be further parametrized into trainable variables $\theta=[\theta^{A},\theta^{B},\theta^{C}]$ with $\theta^{A}\in\mathbb{R}^{(N_{t}+1)\times d\times d}_{\text{sym}}$, $\theta^{B}\in\mathbb{R}^{(N_{t}+1)\times d}$, and $\theta^{C}\in\mathbb{R}^{N_{t}+1}$ after time discretizations. I.e., we only evaluate $\psi$ when $t=t_{j}$ ($j=0,1,\ldots,N_{t}$). Here, $\theta^{A}\in\mathbb{R}^{(N_{t}+1)\times d\times d}_{\text{sym}}$ means that each slice of matrix $\theta^{A}_{j}\in\mathbb{R}^{d\times d}$ is symmetric.
This parametrization is different from equation 8, but the numerical implementation for the ODE discretization equation 9-equation 13 and loss cost equation 31 remain unchanged, with the vector field defined by

$$ $f(t_{j},x;\theta)=\nabla_{x}\psi(t,x;\theta)=\theta^{A}_{j}x+\theta^{B}_{j}\,.$ (84) $$

In the LQ example with Gaussian distribution, the solution is indeed of this form; see Bensoussan et al. () Chapter 6 for details.

With such parametrizations, we can compute the residual of the HJB equation through equation 30 as a regularizer. Here, the term $\partial_{t}\psi(t,z_{t})$ in equation 30 could be approximated using the finite difference scheme. In this work, we use the central difference scheme

$$ $\displaystyle\quad\partial_{t}\psi(t_{j},x;\theta)\approx\dfrac{1}{2\Delta t} \left(\psi(t_{j+1},x;\theta)-\psi(t_{j-1},x;\theta)\right)$ (85) $\displaystyle=\dfrac{1}{2\Delta t}\left(\frac{1}{2}x^{\top}(\theta^{A}_{j+1}- \theta^{A}_{j-1})x+(\theta^{B}_{j+1}-\theta^{B}_{j-1})^{\top}x+(\theta^{C}_{j+ 1}-\theta^{C}_{j-1})\right)\,.$ $$

With this discretization, the numerical residual of the HJB equation at $t_{j}$ (cf. equation 30) is given by

$$ $\displaystyle\mathcal{L}_{\text{HJB}t_{j}}$ $\displaystyle=\dfrac{1}{N_{z}}\sum_{n=1}^{N_{z}}\left[\dfrac{1}{2\Delta t} \left(\psi(t_{j+1},z^{(n)}_{t_{j}};\theta)-\psi(t_{j-1},z^{(n)}_{t_{j}};\theta )\right)+\frac{1}{2}\left|\nabla_{z}\psi(t_{j},z^{(n)}_{t_{j}};\theta)\right|^ {2}\right.$ (86) $\displaystyle\hskip 57.81621pt\left.+\gamma^{2}\left(\operatorname{Tr}\left(H^ {(n)}_{t_{j}}\right)+\frac{1}{2}\left|s^{(n)}_{t_{j}}\right|^{2}\right)\right]\,,$ $$

and the total regularization loss is

$$ $\mathcal{L}_{\text{HJB}}=\sum_{j=1}^{N_{t}-1}\left|\mathcal{L}_{\text{HJB}t_{j }}\right|\Delta t\,.$ (87) $$

Finally, we can minimize $\mathcal{L}_{\text{total}}=\mathcal{L}_{\text{cost}}+\lambda\mathcal{L}_{\text
{HJB}}$ w.r.t. $\theta$ to solve the MFC problem, where $\lambda>0$ is a weight parameter. We summarize the regularized method in Appendix B.2 Algorithm 2

Regularized algorithm for flow matching of overdamped Langevin dynamics
This regularized algorithm could be extended to a more general setting, such as the flow matching problem for an OU process. We recall that the formulation for this problem is presented in Appendix A.3.

We consider the flow matching of OU process with a linear drift function $b(t,x)=-ax$ and an initial Gaussian distribution $N(\mu(0),\Sigma(0))$ where $\mu(0)\in\mathbb{R}^{d}$ and $\Sigma(0)\in\mathbb{R}^{d\times d}_{\text{sym}}$. Then the stochastic state dynamic is given by

$$ $\mathrm{d}x_{t}=-ax_{t}\,\mathrm{d}t+\sqrt{2\gamma}\,\mathrm{d}W_{t}\,,\quad \quad x_{0}\sim N(\mu(0),\Sigma(0))\,.$ (88) $$

Standard calculations inform us that $x_{t}\sim(\mu(t),\Sigma(t))$ remains to be Gaussian, where the mean and variance evolve as

$$ $\mu(t)=\exp(-at)\,\mu(0)\,,$ (89) $$

and

$$ $\Sigma(t)=\dfrac{\gamma}{a}I_{d}+\left(\Sigma(0)-\dfrac{\gamma}{a}I_{d}\right) \exp(-2at)\,.$ (90) $$

This expression provides a reference solution in the numerical test. In this work, we set $a=1$, $\mu(0)$ to be an all-one vector in $\mathbb{R}^{d}$, and $\Sigma(0)=4I_{d}$.

Similar to the algorithm presented above, we parametrize $\psi$ through equation 83. Please note that the solution $\psi$ to the MFC problem is given by Corollary 2. So, the vector field is

$$ $f(t_{j},x;\theta)=\nabla_{x}\psi(t_{j},x;\theta)+b(t_{j},x)=\theta^{A}_{j}x+ \theta^{B}_{j}-ax\,.$ (91) $$

Also, the HJB equation along the state trajectory is written as

$$ $\partial_{t}\psi(t,z_{t})+\frac{1}{2}\left|\nabla_{z}\psi(t,z_{t})\right|^{2}+ \nabla_{z}\psi(t,z_{t})^{\top}b(t,z_{t})-\gamma\nabla_{z}\cdot b(t,_{t})+ \gamma^{2}\left(\operatorname{Tr}(H_{t})+\frac{1}{2}|s_{t}|^{2}\right)=0\,.$ (92) $$

One can view equation 92 as an analog of equation 30 in the flow matching example.

The time discretization for the HJB equation is the same as equation 85, and the residual of the HJB equation at $t_{j}$ is given by

$$ $\displaystyle\mathcal{L}_{\text{HJB}t_{j}}$ $\displaystyle=\dfrac{1}{N_{z}}\sum_{n=1}^{N_{z}}\left[\dfrac{1}{2\Delta t} \left(\psi(t_{j+1},z^{(n)}_{t_{j}};\theta)-\psi(t_{j-1},z^{(n)}_{t_{j}};\theta )\right)+\frac{1}{2}\left|\nabla_{z}\psi(t_{j},z^{(n)}_{t_{j}};\theta)\right|^ {2}\right.$ (93) $\displaystyle\left.+\nabla_{z}\psi(t_{j},z^{(n)}_{t_{j}};\theta)^{\top}b(t_{j} ,z^{(n)}_{t_{j}})-\gamma\nabla_{z}\cdot b(t_{j},z^{(n)}_{t_{j}})+\gamma^{2} \left(\operatorname{Tr}\left(H^{(n)}_{t_{j}}\right)+\frac{1}{2}\left|s^{(n)}_{ t_{j}}\right|^{2}\right)\right]\,,$ $$

and the total regularization loss is still

$$ $\mathcal{L}_{\text{HJB}}=\sum_{j=1}^{N_{t}-1}\left|\mathcal{L}_{\text{HJB}t_{j }}\right|\Delta t\,.$ (94) $$

And the total loss is again $\mathcal{L}_{\text{total}}=\mathcal{L}_{\text{cost}}+\lambda\mathcal{L}_{\text
{HJB}}$. We summarize this method in Algorithm 3 Appendix B.2.

### B.2 Pseudo-code for algorithms

We present all pseudo-codes for all algorithms in this section. Recall that the standard version of the MFC solver Algorithm 1 is presented in the main text. We present the regularized MFC solver for LQ problem Algorithm 2; the regularized MFC solver for the flow matching of OU process Algorithm 3; and the multi-stage splicing method Algorithm 4.

Figure: Algorithm 2 Regularized score-based normalizing flow solver for the MFC problem

Figure: Algorithm 3 Regularized score-based normalizing flow solver for flow matching of OU process

Figure: Algorithm 4 Multi-stage splicing algorithm for MFC problems

### B.3 The errors in the numerical examples

The errors for the standard Algorithm 1. For examples with exact solutions, we present the errors for the density, the velocity field, and the score function. These numerical expressions are given by

$$ $\displaystyle\text{err}_{\rho}$ $\displaystyle=\dfrac{1}{N_{z}}\sum_{n=1}^{N_{z}}\left|\tilde{l}^{(n)}_{t_{N_{t }}}-\rho(t_{N_{t}},z^{(n)}_{t_{N_{t}}})\right|\,,$ (95) $\displaystyle\text{err}_{f}$ $\displaystyle=\dfrac{1}{N_{z}}\dfrac{1}{N_{t}+1}\sum_{n=1}^{N_{z}}\sum_{j=0}^{ N_{t}}\left|f\left(t_{j},z^{(n)}_{t_{j}};\theta\right)-f\left(t_{j},z^{(n)}_{t _{j}}\right)\right|\,,$ $\displaystyle\text{err}_{s}$ $\displaystyle=\dfrac{1}{N_{z}}\sum_{n=1}^{N_{z}}\left|s^{(n)}_{t_{N_{t}}}- \nabla_{z}\log\rho\left(t_{N_{t}},z^{(n)}_{t_{N_{t}}}\right)\right|\,.$ $$

The errors we present in the tables in section 5 and Appendix C are the average errors over $10$ independent runs. Note that there are two sources that contribute to these three errors $\text{err}_{\rho}$, $\text{err}_{f}$, and $\text{err}_{s}$. We take the terminal time $t_{t_{N}}=t_{\text{end}}$ as an example. The first source is the discretization error for the state dynamic, characterized by $z_{t_{\text{end}}}-z_{t_{N_{t}}}$, where $z_{t_{\text{end}}}$ is the end point of the continuous dynamic equation 3a and $z_{t_{N_{t}}}$ is its time discretization. The second source comes from the training error of the neural network $f(\cdot,\cdot;\theta)$. When a function is steep, especially the score function $\nabla\log\rho$ in a non-Gaussian example, a small discretization error could be magnified by a large factor due to this large condition number. This is the reason why the score errors for the double well potential function example in Appendix C.4 are larger compared with other results.

The errors for the double well potential example.
For the double well potential example, it is hard to obtain a reference solution for $t\in(0,t_{\text{end}})$ (cf. Appendix C.4). As an alternative, we compute the error of the velocity field at initial and terminal time $t=0,t_{\text{end}}$ through

$$ $\displaystyle\text{err}_{f_{0}}$ $\displaystyle=\dfrac{1}{N_{z}}\sum_{n=1}^{N_{z}}\left|f\left(0,z^{(n)}_{0}; \theta\right)-f\left(0,z^{(n)}_{0}\right)\right|\,,$ (96) $\displaystyle\text{err}_{f_{end}}$ $\displaystyle=\dfrac{1}{N_{z}}\sum_{n=1}^{N_{z}}\left|f\left(t_{N_{t}},z^{(n)} _{t_{N_{t}}};\theta\right)-f\left(t_{N_{t}},z^{(n)}_{t_{N_{t}}}\right)\right|\,.$ $$

The errors for the regularized Algorithm 2 and 3. For modified algorithms with HJB regularizers, including Algorithm 2 and 3, we compute the error of the parameter $A$ and $B$ directly, given by

$$ $\displaystyle\text{err}_{A}$ $\displaystyle=\dfrac{1}{N_{t}+1}\sum_{j=0}^{N_{t}}\left|\theta^{A}_{j}-A(t_{j} )\right|\,,$ (97) $\displaystyle\text{err}_{B}$ $\displaystyle=\dfrac{1}{N_{t}+1}\sum_{j=0}^{N_{t}}\left|\theta^{B}_{j}-B(t_{j} )\right|\,,$ (98) $$

where $A(\cdot)$ and $B(\cdot)$ denote the true solutions. Here, we do not present the error of $\theta^{C}$ for two reasons. First, according to equation 84 and equation 91, the velocity field does not depend on $\theta^{C}$. Second, we observe from equation 86 and equation 93 that a constant shift of $\theta^{C}$ does not affect the regularization loss. Looking more carefully at equation 85, a constant shift of all the even (or odd) indices of $\theta^{C}$ also does not change the loss equation 86 or equation 93. Therefore, it is not meaningful to compute the error for $\theta^{C}$.

### B.4 Hyperparameters for the numerical examples

We present all the hyperparameters of the numerical tests in this section.

We set the diffusion coefficient $\gamma=1$ for all problems, except for the double well potential example, where $\gamma=0.1$. The reason for this difference is because if we set $\gamma=1$ in the double well example, the terminal distribution $\rho(1,\cdot)$ will be very flat and we will not have a clear bifurcation pattern. In the LQ example in Appendix C.3, we set $\beta=0.1$.

For all the numerical tests, we use a step size $\Delta t=0.01$. The terminal time $t_{\text{end}}=1$ for all the tests except for the double moon example. In the double moon example, we have a multi-stage splicing algorithm, with the first stage being $t\in[0,0.2]$ and the second stage being $t\in[0.2,0.4]$.

For the neural network parametrization equation 8, we use $\tanh$ as the activation function. The widths of the network are $20$, $50$, and $200$ in $1$, $2$, and $10$ dimensions for the RWPO and LQ problems. For the double moon flow matching problem (in $2$ dimensions), the width is $200$. For the double well example, the widths are $50$, $100$, and $200$ in $1$, $2$, and $10$ dimensions. The double moon flow matching problem and the double well problem are harder due to the bifurcation pattern.

For the regularized algorithms (Algorithm 2 and 3), we set the weight parameter for the HJB loss as $\lambda=$1\text{\times}{10}^{-3}$$.

For all the examples, we train the model using the built-in Adam optimizer in Pytorch, with a learning rate $0.01$. The number of training steps is set such that the training loss is roughly stable. The number of training steps are $200$, $200$, and $300$ in $1$, $2$, and $10$ dimensions for both RWPO and LQ problems with Algorithm 1, regularized WPO with Algorithm 2, flow matching for OU process with Algorithm 3. The number of steps for flow matching double moon example and the double well example in $1$, $2$, and $10$ dimensions, the number of training steps is $500$.

## Appendix C Additional numerical results

### C.1 RWPO

For the RWPO problem, after the transformation to probability flow through equation 22, the problem becomes

$$ $\inf_{f}\int_{0}^{1}\int_{\mathbb{R}^{d}}\frac{1}{2}\left|f(t,x)+\gamma\nabla_ {x}\log(\rho(t,x))\right|^{2}\rho(t,x)\,\mathrm{d}x\,\mathrm{d}t+\int_{\mathbb {R}^{d}}G(x)\rho(1,x)\,\mathrm{d}x\,,$ $$

subject to the transport equation

$$ $\partial_{t}\rho(t,x)+\nabla_{x}\cdot(\rho(t,x)f(t,x))=0\,.$ $$

The optimal density evolution is given by

$$ $\rho(t,x)=\left(4\pi(2-t)\gamma\right)^{-\frac{d}{2}}\exp\left(-\dfrac{\left|x \right|^{2}}{4\gamma(2-t)}\right)\,,$ $$

and the optimal velocity field is

$$ $f(t,x)=-\dfrac{x}{2(2-t)}\,.$ $$

In addition to the errors and Figure 1, we present additional numerical results in Figure 2. The first, second, and third rows of Figure 2 show numerical results in $1$, $2$, and $10$ dimensions, respectively. The first column shows the curves of the cost function through training. The light blue shadows plot the standard deviation observed during $10$ independent test runs. The red dashed lines represent the cost under the optimal control field. We observe that our algorithm nearly reaches the optimal cost.

The second column in Figure 2 plots the first dimension of the velocity field at $t=0$. Our neural networks (in blue) accurately captures the true velocity (in orange). The third column presents similar plots for the velocity fields, but at $t=t_{\text{end}}=1$. Our algorithm also computes the velocity fields accurately.

The fourth column in the first row of Figure 2 plots the score function (in blue) at $t=1$ using the data points $\{(z^{(n)}_{t_{N_{t}}},s^{(n)}_{t_{N_{t}}})\}_{n=1}^{N_{z}}$. This curve coincides with the true score function in orange. Note that the score function at $t=0$ is given, so we only plot the score at a terminal time. The second and third figures in the fourth column of Figure 2 show a density plot of the first dimension of the velocity field. These are probability density functions of $f(0,z_{0};\theta)$ and $f(0,z_{0})$, where $z_{0}\sim\rho_{0}$. The density functions are further approximated by the histogram of samples. Such techniques for visualization of high dimensional functions have been applied by Han et al. (). We observe that our neural network accurately captures the velocity fields.

Figure: Figure 2: The first, second and third row: numerical results for RWPO in $1d$, $2d$, and $10d$. First column: cost functional through training. Second and third columns: first dimension of the velocity field at $t=0$ and $t=1$. Fourth column: the first row is the score function at $t=1$, the second and third rows are the density plots for the first dimension of the velocity field. Our numerical results accurately capture the true solutions.
Refer to caption: x4.png

### C.2 Flow matching

In this section, we first present additional numerical results for the flow matching of the OU process. Then, we describe the details of the numerical implementation, such as the double moon potential, and present the numerical results.

Flow matching for OU process. The numerical results of flow matching for OU process are shown in Figure 3 and 4. Figure 3 shows the density evolution of the process in $1$ dimension, which is similar to the second plot in Figure 1. Our algorithm captures both the change of mean value and the shrink of variance accurately. Figure 4 shows the particle dynamic under the trained probability flow in $2$ dimensions and its comparison with the stochastic OU process. We add level sets of the density function with center $\mu(t)$ and radius $\frac{1}{4}\sigma(t)$, $\frac{1}{2}\sigma(t)$, and $\sigma(t)$ for better visualization, where $\sigma(t):=\det(\Sigma(t))^{\frac{1}{2d}}$ denotes the standard deviation. We pick points that are initialized within the circle with center $\mu(0)$ and radius $\sigma(0)$, and record the trajectories of these particles. According to Figure 4, our algorithm accurately captures the change of mean and variance of the OU process. Additionally, compared with the OU process, our probability flow ODE demonstrates structured behaviors.

Figure: Figure 3: Flow matching for OU process in $1$ dimension. Our algorithm accurately captures the density evolution of the state dynamic, including the change of mean and variance.
Refer to caption: x16.png

Figure: Figure 4: Flow matching for OU process in $2$ dimension. First line: particle evolution under the probability flow $z_{t}$ with trained velocity field. Second line: simulation of the stochastic OU process. The circles in the figure are centered at $\mu(t)$, with radius $\frac{1}{4}\sigma(t)$, $\frac{1}{2}\sigma(t)$, and $\sigma(t)$. Our algorithm accurately captures the density evolution of the state dynamic, including the change of mean and variance. The probability flow demonstrate more structured behavior compared with the stochastic OU process.
Refer to caption: x17.png

Double moon example.
In the example with double moon potential, given by

$$ $V(x)=2(|x|-3)^{2}-2\log\left[\exp\left(-2\left(x_{1}-3\right)^{2}\right)+\exp \left(-2\left(x_{1}+3\right)^{2}\right)\right]\,.$ (99) $$

This example has been computed by Wang & Li (); Tan et al. (). We aim to learn the probability flow ODE of the overdamped Langevin dynamic

$$ $\partial_{t}\rho(t,x)-\nabla_{x}\cdot(\rho(t,x)\nabla_{x}V(t,x))=\gamma\Delta_ {x}\rho(t,x)\,,$ $$

where the initial distribution $\rho_{0}$ is $N(0,1)$. There are two moon-shaped patterns from this potential function. Therefore, the state dynamic has a bifurcation phenomenon, which is usually hard to capture.

We set the total time span as $t\in[0,0.4]$. The overdamped Langevin dynamic is already close to its stable distribution at $t=0.4$; see the last scattered plot in Figure 5. Also, our algorithm has demonstrated its ability to learn the dynamic accurately within a longer interval in other examples.

Next, we present the detailed implementation of the multi-stage splicing method for our double moon example. In this toy example, we partition the interval into two stages (sub-intervals) $[0,0.2]$ and $[0.2,0.4]$. We apply Algorithm 1 within the first interval. After training, we save the trained network as a warm-up (initialization) for the training in the next stage. We also record the particles $z_{t}^{(n)}$ at terminal time $t=0.2$ in the first stage, which serves as the initial distribution for training in the next stage.

In the second stage $t\in[0.2,0.4]$, the training process is similar. We inherit the network parameter $\theta$ from the previous stage as initialization. Also, the state particles $z_{t}^{(n)}$ at terminal time from the last stage are utilized as initialization distribution for the new stage. Note that starting form the second stage, there is no longer resampling because we only have finite samples at $t=0.2$. As a consequence, we may encounter the problem of overfitting, which is also known as model collapse. In this work, we apply a $L^{2}$ regularization with weight $0.1$ to avoid this issue. If we only have finite samples for the initial distribution, then regularization should also be added in the first stage of training. We summarize this multi-stage splicing method in Algorithm 4.

The numerical result for the splicing method is shown in Figure 5. The first row shows the particle dynamic of state $z_{t}$ under the trained velocity field within the first stage, which coincides with an overdamped Langevin dynamic in the second row. The third row shows the particle dynamic in the second stage after training, coinciding with the overdamped Langevin dynamic in the fourth row. We also add the level sets for the density function of the stationary distribution for better visualization. These results confirm that our multistage splicing method can capture the Fokker-Planck equation of an overdamped Langevin dynamic in a total time span.

Figure: Figure 5: 2D flow matching double moon. First row: evolution of particles under the trained probability flow $z_{t}$ in the first stage $[0,0.2]$. Second row: evolution of particles under the overdamped Langevin dynamic in the first stage $[0,0.2]$. Third row: evolution of particles under the trained probability flow $z_{t}$ in the second stage $[0.2,0.4]$. Fourth row: evolution of particles under the overdamped Langevin dynamic in the second stage $[0.2,0.4]$. Our multi-stage splicing method captures the overdamped Langevin dynamic correctly.
Refer to caption: x18.png

### C.3 An LQ example with an entropy potential cost

We consider an LQ example in this section. This example was also studied by Lin et al. (). We minimize the cost functional

$$ $\inf_{v}\int_{0}^{t_{\text{end}}}\int_{\mathbb{R}^{d}}\left(\frac{1}{2}\left|v (t,x)\right|^{2}+\frac{1}{2}|x|^{2}+\beta\log(\rho(t,x))\right)\rho(t,x)\, \mathrm{d}x\,\mathrm{d}t+\int_{\mathbb{R}^{d}}G(x)\rho(t_{\text{end}},x)\, \mathrm{d}x\,,$ $$

subject to

$$ $\partial_{t}\rho(t,x)+\nabla_{x}\cdot(\rho(t,x)v(t,x))=\Delta_{x}\rho(t,x)\,, \quad\quad\rho(0,x)=\rho_{0}(x)\,.$ $$

With a score substitution, the problem is equivalent to

$$ $\displaystyle\inf_{f}\int_{0}^{t_{\text{end}}}\int_{\mathbb{R}^{d}}$ $\displaystyle\left(\frac{1}{2}\left|f(t,x)+\nabla_{x}\log(\rho(t,x))\right|^{2 }+\frac{1}{2}|x|^{2}+\beta\log(\rho(t,x))\right)\rho(t,x)\,\mathrm{d}x\, \mathrm{d}t$ $\displaystyle+\int_{\mathbb{R}^{d}}G(x)\rho(t_{\text{end}},x)\,\mathrm{d}x\,,$ $$

subject to

$$ $\partial_{t}\rho(t,x)+\nabla_{x}\cdot(\rho(t,x)f(t,x))=0\,,\quad\quad\rho(0,x) =\rho_{0}(x)\,.$ $$

We define $\alpha:=\left(\sqrt{\beta^{2}+4}-\beta\right)/2$. The initial distribution $\rho_{0}$ is Gaussian $N(0,\frac{1}{\alpha}I_{d})$ and the terminal cost is $G(x)=\frac{\alpha}{2}|x|^{2}$. The optimal density evolution is given by

$$ $\rho(t,x)=\left(\dfrac{\alpha}{2\pi}\right)^{d/2}\exp\left(-\dfrac{\alpha\left |x\right|^{2}}{2}\right)\,.$ $$

We set $\beta=0.1$ and test Algorithm 1 on this example. We remark that we have a term $\beta\log\rho(t,x)$ in the running cost $F$. Therefore, it is better to compute $l_{t}=\log\rho(t,z_{t})$ instead of $\tilde{l}_{t}=\rho(t,z_{t})$ in algorithm 1. Similar to the WPO example, we test our algorithm in $1$, $2$, and $10$ dimensions. The errors are summarized in Table 5, which is similar to the results in Table 2. Our algorithm is able to solve the MFC problem accurately.

The numerical results are also presented in Figure 6. The first, second, and third rows show the results in $1$, $2$, and $10$ dimensions, respectively. The first column is the training curve for the cost functional, which is similar to the first column in Figure 2. Our algorithm gives the correct objective. The second column visualize the density function at terminal time $t_{\text{end}}$ using the data points $\{z^{(n)}_{t_{N_{t}}},\tilde{l}^{(n)}_{t_{N_{t}}}\}_{n=1}^{N_{z}}$. In the $1$ dimension (first row), we plot the density and compare it with true densities directly. In $2$ dimensions, we interpolate the data points and plot the density function $x_{1}\mapsto\rho(t_{\text{end}},(x_{1},0))$ in the figure. In $10$ dimensions, an interpolation is very hard to obtain, so we plot the density in term of the norm of $x$. Our algorithm captures the density function accurately. The third column visualize the score function at terminal time $t_{\text{end}}$ using the data points $\{(z^{(n)}_{t_{N_{t}}},s^{(n)}_{t_{N_{t}}})\}_{n=1}^{N_{z}}$. Again, we plot the score function directly in $1$ dimension. In $2$ and $10$ dimensions, we present a density plot for the first dimension of the score, which gives the probability density function of $\nabla_{z_{1}}\log\rho(t_{\text{end}},z_{t_{\text{end}}})$ (similar to the density plot in Appendix C.1). We observe that our algorithm captures the true density and score functions accurately.

**Table 5: Errors for LQ problem.**
| errors | $\text{err}_{\rho}$ | $\text{err}_{f}$ | $\text{err}_{s}$ | cost gap |
| --- | --- | --- | --- | --- |
| $1d$ | $9.91\text{\times}{10}^{-3}$ | $2.61\text{\times}{10}^{-2}$ | $2.88\text{\times}{10}^{-2}$ | $1.07\text{\times}{10}^{-2}$ |
| $2d$ | $7.39\text{\times}{10}^{-3}$ | $3.74\text{\times}{10}^{-2}$ | $3.58\text{\times}{10}^{-2}$ | $1.97\text{\times}{10}^{-2}$ |
| $10d$ | $1.15\text{\times}{10}^{-3}$ | $2.49\text{\times}{10}^{-2}$ | $2.37\text{\times}{10}^{-2}$ | $1.64\text{\times}{10}^{-1}$ |

Figure: Figure 6: Numerical results for the LQ problem. The first, second, and third rows shows the results in $1$, $2$, and $10$ dimensions respectively. First column: cost functional through training. Second column: visualization for the density function. Third column: visualization for the score function. Our algorithm accurately captures the solution to the problem.
Refer to caption: x20.png

### C.4 Double well potential

In this section we present an MFC example where the terminal cost is a potential function with the double well shape. The formulation is similar to the WPO, with the terminal cost given by

$$ $G(x)=c\left|x-c_{1}\right|^{2}\left|x-c_{2}\right|^{2}\,,$ (100) $$

where $c\in\mathbb{R}^{+}$ and $c_{1},c_{2}\in\mathbb{R}^{d}$ are the minimizers of $G$. This function $G$ does not depend on the density $\rho$. The initial distribution is still a standard Gaussian distribution, i.e., $z_{0}\sim N(0,I_{d})$. As mentioned before, such problem will demonstrate a bifurcation phenomenon, which is hard to capture. We compute this example in $1$, $2$, and $10$ dimensions, and compute all errors in $1$ and $2$ dimensions. We take $c=\frac{1}{4}$ in equation 100 and set two centers as an all $-1$ vector for $c_{1}$ and an all $1$ vector for $c_{2}$.

Computing the reference solution.
Unlike the example in section 5.1, we do not have an explicit solution. But fortunately, we are able to obtain a reference solution using the kernel formula proposed by Li et al. () (equation (10), (11), and (14)). According to their derivations, the optimal density function is given by

$$ $\rho(t,x)=\left(4\pi\gamma\frac{t(t_{\text{end}}-t)}{t_{\text{end}}}\right)^{- \frac{d}{2}}\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}\frac{\exp\left[-\frac{1 }{2\gamma}\left(G(z)+\frac{|x-z|^{2}}{2(t_{\text{end}}-t)}+\frac{|x-y|^{2}}{2t }\right)\right]}{\int_{\mathbb{R}^{d}}\exp\left[-\frac{1}{2\gamma}\left(G( \tilde{y})+\frac{|y-\tilde{y}|^{2}}{2t_{\text{end}}}\right)\right]\,\mathrm{d} \tilde{y}}\rho(0,y)\,\mathrm{d}y\,\mathrm{d}z\,.$ (101) $$

Specifically, the density function as terminal time $t_{\text{end}}$ is

$$ $\rho(t_{\text{end}},x)=\int_{\mathbb{R}^{d}}\frac{\exp\left[-\frac{1}{2\gamma} \left(G(x)+\frac{|x-y|^{2}}{2t_{\text{end}}}\right)\right]}{\int_{\mathbb{R}^{ d}}\exp\left[-\frac{1}{2\gamma}\left(G(z)+\frac{|z-y|^{2}}{2t_{\text{end}}} \right)\right]\,\mathrm{d}z}\rho(0,y)\,\mathrm{d}y\,.$ (102) $$

In addition, the solution to the classic HJB equation (cf. equation 75) is given by

$$ $\phi(t,x)=2\gamma\log\left(\int_{\mathbb{R}^{d}}(4\pi\gamma(t_{\text{end}}-t)) ^{-\frac{d}{2}}\exp\left[-\frac{1}{2\gamma}\left(G(y)+\frac{|x-y|^{2}}{2(t_{ \text{end}}-t)}\right)\right]\,\mathrm{d}y\right)\,,$ (103) $$

with a terminal condition $\phi(T,x)=-G(x)$. With these expressions, we are able to obtain the score function via $\nabla_{x}\log\rho(t,x)=\nabla_{x}\rho(t,x)/\rho(t,x)$. Then, by Corollary 1 and equation 76, the optimal velocity is

$$ $f(t,x)=\nabla_{x}\phi(t,x)-\gamma\nabla_{x}\log\rho(t,x)\,.$ (104) $$

However, when $t\in(0,t_{\text{end}})$, the numerical implementation for $\rho(t,x)$ and $\nabla_{x}\rho(t,x)$ through equation 101 involves three nested integrations in $\mathbb{R}^{d}$, which could potentially result in large errors. Therefore, we instead consider the errors at $t=0$ and $t_{\text{end}}$, where all expressions are relatively easy. At $t=0$, by equation 104, we have

$$ $f(0,x)=\dfrac{\int_{\mathbb{R}^{d}}\frac{x-y}{t_{\text{end}}}\exp\left[-\frac{ 1}{2\gamma}\left(G(y)+\frac{|x-y|^{2}}{2t_{\text{end}}}\right)\right]\,\mathrm {d}y}{\int_{\mathbb{R}^{d}}\exp\left[-\frac{1}{2\gamma}\left(G(z)+\frac{|x-z|^ {2}}{2t_{\text{end}}}\right)\right]\,\mathrm{d}z}-\gamma\dfrac{\nabla_{x}\rho_ {0}(x)}{\rho_{0}(x)}\,.$ (105) $$

At $t=t_{\text{end}}$, we can compute the score function through

$$ $\nabla_{x}\log\rho(t_{\text{end}},x)=\dfrac{-\int_{\mathbb{R}^{d}}\frac{-\frac {1}{2\gamma}\left(\nabla_{x}G(x)+\frac{x-y}{t_{\text{end}}}\right)\exp\left[- \frac{1}{2\gamma}\left(G(x)+\frac{|x-y|^{2}}{2t_{\text{end}}}\right)\right]}{ \int_{\mathbb{R}^{d}}\exp\left[-\frac{1}{2\gamma}\left(G(z)+\frac{|z-y|^{2}}{2 t_{\text{end}}}\right)\right]\,\mathrm{d}z}\rho(0,y)\,\mathrm{d}y}{\int_{ \mathbb{R}^{d}}\frac{\exp\left[-\frac{1}{2\gamma}\left(G(x)+\frac{|x-y|^{2}}{2 t_{\text{end}}}\right)\right]}{\int_{\mathbb{R}^{d}}\exp\left[-\frac{1}{2 \gamma}\left(G(z)+\frac{|z-y|^{2}}{2t_{\text{end}}}\right)\right]\,\mathrm{d}z }\rho(0,y)\,\mathrm{d}y}\,.$ (106) $$

Also, the terminal velocity is given by

$$ $f(t_{\text{end}},x)=-\nabla_{x}G(x)-\gamma\nabla_{x}\log\rho(t_{\text{end}},x)\,.$ (107) $$

Next, we apply the Riemann sum to approximate the expressions equation 105, equation 106, and equation 107 to obtain a reference solution in $1$ and $2$ dimensions. First, we define the following function

$$ $h(y)=\int_{\mathbb{R}^{d}}\exp\left[-\frac{1}{2\gamma}\left(G(z)+\frac{|z-y|^{ 2}}{2t_{\text{end}}}\right)\right]\,\mathrm{d}z\,,$ (108) $$

which has appeared multiple times. We approximated $h(y)$ numerically from the Riemann sum. In $1$ dimension, we use the trapezoid rule with step size $\Delta z=0.01$ to approximate the integrations and truncate the integration within $z\in[-6,6]$ when we compute the reference solution. In $2$ dimensions, similarly, we use a box with size $\Delta z_{1}\times\Delta z_{2}=0.01\times 0.01$ for Riemann sum and truncate the integration within the box $z\in[-6,6]\times[-6,6]$. We compute the value of $h(y)$ on all grid points within $[-4,4]$ in $1$ dimension and $[-4,4]\times[-4,4]$ in $2$ dimensions, and store the values. In this way, integrations w.r.t. $z$ in equation 105 and equation 106 are simulated numerically. Next, we further apply Riemann sum for integrations w.r.t. $y$ with the same step size, truncated for $y\in[-4,4]\,\text{or}\,[-4,4]\times[-4,4]$. In this way, we obtain the function for all $x$ on the grid points. The two bounds $6$ and $4$, and the step size $0.01$, are chosen such that the Riemann sums provide reasonable reference solutions for the problem. Finally, we make an interpolation in the dimensions $1$ or $2$, giving us the desired reference solution.

Results in $1$ dimension.
The errors for $1$ dimensional example are $\text{err}_{\rho}=$1.17\text{\times}{10}^{-2}$$, $\text{err}_{f_{0}}=$8.56\text{\times}{10}^{-3}$$, $\text{err}_{f_{\text{end}}}=$2.39\text{\times}{10}^{-2}$$, and $\text{err}_{s}=$1.65\text{\times}{10}^{-1}$$. The numerical results are illustrated in Figure 7. The two plots in the first row show the velocity fields at $t=0$ and $t=t_{\text{end}}=1$, which nicely approximate the true velocities. In the second row, the figure on the left compares the reference terminal density and the histogram of particles $z_{t_{\text{end}}}$, which match well. The plot in the middle shows the evolution of the density, which is obtained using data points $\{(z_{t_{j}}^{(n)},\tilde{l}_{t_{j}}^{(n)})\}_{j=0,n=1}^{N_{t},\,N_{z}}$. These red lines coincide with the blue dashed lines, which are the true densities. This result demonstrates that our algorithm is able to capture the evaluation from a simple distribution to a complicated double well distribution accurately. The figure on the right plots the score function using data points $\{(z_{t_{N_{t}}}^{(n)},\tilde{l}_{t_{N_{t}}}^{(n)})\}_{n=1}^{N_{z}}$. It successfully captures the S-shape of the true score function.

Figure: Figure 7: double well 1d. Left: true density and particle histogram at $T$; middle: density evolution compared with true density; right: score function at $T$.
Refer to caption: x29.png

Results in $2$ dimensions.
It is hard to compute the nested integrations in Li et al. ()’s paper using Riemann sums in $2$ dimensions. Therefore, instead of computing $\text{err}_{f}$ in equation 95, we compute the errors for the velocity field at $t=0$ and $t=t_{\text{end}}=1$ using the expressions in equation 96. The errors are $\text{err}_{\rho}=$2.62\text{\times}{10}^{-2}$$, $\text{err}_{f_{0}}=$5.11\text{\times}{10}^{-2}$$, $\text{err}_{f_{\text{end}}}=$4.77\text{\times}{10}^{-2}$$, and $\text{err}_{s}=$2.4\text{\times}{10}^{-1}$$.

The results are also shown in Figure 8 and 9. Our learned network successfully captures the velocity field $f$ and the score function which is shown in Figure 8. The first and second rows show each dimension of the velocity field at $t=0$ and $t=t_{\text{end}}=1$ and compare them with the true values. The third row shows the score function. The approximated score function is plotted using the data $\{(z^{(n)}_{t_{N_{t}}},s^{(n)}_{t_{N_{t}}})\}_{n=1}^{N_{z}}$, where we made a non-uniform $2d$ interpolation using the build-in function in Scipy. The reference solution is obtained through differentiation of equations (10) (11) and (14) given by Li et al. () and then apply Riemann sum approximation. We observe that our algorithm correctly captures the velocity field and the score function.

The left plot in Figure 9 shows the scattered points at terminal time $t=1$, which concentrate at the two wells of the potential. We also add the level sets of $\rho(1,\cdot)$ at $\{0.01,0.05,0.1,0.2,0.4,0.6\}$ for better visualization. Note that the expression for $\rho(1,\cdot)$ is given by equation 102 instead of $\frac{1}{Z}\exp(-G(x)/\gamma)$. The plot in the middle shows the particle trajectories under the trained velocity field, which shows a bifurcation phenomenon.

Figure: Figure 8: Image plot for $f(0,\cdot)$ (first row), $f({t_{\text{end}}},\cdot)$ (second row), and score function $\nabla_{x}\log\rho({t_{\text{end}}},\cdot)$. Columns one to four represent the true value in the first dimension, approximated value in the first dimension, true value in the second dimension, and the approximated value in the second dimension. Our algorithm correctly captures the velocity field and score function.
Refer to caption: x34.png

Figure: Figure 9: Left: scattered plot for particles at $t=t_{\text{end}}$ and the level sets of the density for $2$ dimensional example. Middle: particle trajectory plots for $2$ dimensional example, which demonstrates a bifurcation phenomenon. Right: projected scattered plot of particles at $t=t_{\text{end}}$ for $10$ dimensional example.
Refer to caption: x37.png

Results in $10$ dimensions.

Finally, we also test our algorithm in $10$ dimensions. We do not have a reference solution in this case, but we observe that the particles go into two piles, as is shown in the right plot of Figure 9. In this plot, we project the $10$ dimensional particles into $2$ dimensions for visualization.

## Appendix D remarks

We make a few remarks in this section.

The score dynamic equation 3c. Boffi & Vanden-Eijnden () proposed a similar equation (in equation (13)) before. But we have concrete numerical example. Also, we believe that we are the first to propose the ODE dynamic for the second-order score function equation 3d. This equation significantly enhances the accuracy of the regularized MFC solver, as demonstrated in Table 3 and 4.

Parametrization of neural network. In this paper, we parametrize $f$ as a neural network with a single (hidden) layer. As a consequence, the whole map from $z_{0}$ to $z_{t_{\text{end}}}$ is slightly different from the typical structure of a deep neural network, as described by Chen et al. (). The layer structure for $f$ is affine-activation-affine, so the composition of multiple layers becomes affine-activation-affine-affine-activation-affine. In contrast, a typical structure for a deep neural network is affine-activation-affine-activation-affine. The composition of two affine functions is still an affine function, so it is not a big difference.

Related work on theoretical MFC. From a theoretical perspective, convergence properties of MFC problems have been extensively studied by Carmona & Laurière (; ). In addition, Lacker () investigated the limit behavior as the number of agents approaches infinity. The maximum principle, a crucial concept in optimal control, has also been extended to the mean field setting by Meyer-Brandis et al. (), offering an understanding of the theoretical underpinnings of MFC problems. The mean field game (MFG) problem (Cardaliaguet, ) is closely related to MFC. Cardaliaguet et al. () generalize the HJB equation to the master equation, whose properties have been further investigated by Gangbo & Mészáros (); Gangbo et al. (). Bayraktar et al. () study the convergence of particle system on Wasserstein space.

HJB regularization for MFC solver.
For the HJB regularization technique in B.1, we used absolute value instead of a squared residual in equation 87 because it performs slightly better. In fact, both forms of regularization could significantly reduce errors.

As for the reason why it reduces errors, we recall that we minimize the discretized loss function in the original Algorithm 1. As a consequence, there is a trend of overfitting to the discretized loss. Adding a regularizer could resolve this overfitting issue. Such issues for discretizations are also discussed in Zhou & Lu ().

Computational cost for the score function via equation 17 and equation 18.
When we compute the score function through equation 17, the term $\nabla_{x}\log\left(\det(\nabla_{x}T(t,x))\right)$ could be tricky. Note that $\nabla_{x}\log\left(\det(\nabla_{x}T(t,x))\right)=\nabla_{x}\det(\nabla_{x}T(t
,x))/\det(\nabla_{x}T(t,x))$, and the $i$-th component of $\nabla_{x}\det(\nabla_{x}T(t,x))$ can be computed through

$$ $\partial_{x_{i}}\det(\nabla_{x}T(t,x))=\det(\nabla_{x}T(t,x))\operatorname{Tr} (\nabla_{x}T(t,x)^{-1}\,\partial_{x_{i}}\nabla_{x}T(t,x)).$ $$

Therefore, the $i$-th component of $\nabla_{x}\log\left(\det(\nabla_{x}T(t,x))\right)$ in equation 17 is

$$ $\operatorname{Tr}(\nabla_{x}T(t,x)^{-1}\,\partial_{x_{i}}\nabla_{x}T(t,x))\,.$ $$

As a consequence, the cost for computing one single score function $s_{t}$ through equation 17 is either from computing the inverse of Jacobi matrix or solving the related linear system for each dimension. Classical numerical linear algebra methods suggest that the complexity is $\mathcal{O}(d^{3})$, and the cost for one score trajectory $\{s_{t_{j}}\}_{j=0}^{N_{t}}$ is $\mathcal{O}(N_{t}\,d^{3})$.

We can also compute the score function through equation 18, where we assume that the width of the neural network is $\mathcal{O}(d)$, i.e., proportional the dimension. In this way, $\dfrac{\partial z_{t}}{\partial z_{0}}$ and $\dfrac{\partial l_{t}}{\partial z_{0}}$ are computed from auto-differentiations of the deep neural network functions in equation 15. The computational cost for a single score $s_{t_{j}}$ is $\mathcal{O}(j\,d^{3})$, because we need to compute the back propagation of the deep neural network equation 15 with depth $j$ using chain rule. However, the total cost for computing a whole trajectory $\{s_{t_{j}}\}_{j=0}^{N_{t}}$ is $\mathcal{O}(N_{t}\,d^{3})$ rather than $\mathcal{O}(N_{t}^{2}\,d^{3})$ because we can efficiently store intermediate computations. For example, according to equation 9, we can compute

$$ $\dfrac{\partial z_{t_{j+1}}}{\partial z_{0}}=\dfrac{\partial z_{t_{j+1}}}{ \partial z_{t_{j}}}\dfrac{\partial z_{t_{j}}}{\partial z_{0}}=\left(I_{d}+ \Delta t\,\nabla_{z}f(t_{j},z_{t_{j}};\theta)\right)\dfrac{\partial z_{t_{j}}} {\partial z_{0}}$ (109) $$

with additional $\mathcal{O}(d^{2})$ operations if $\dfrac{\partial z_{t_{j}}}{\partial z_{0}}$ is stored. Similarly, by equation 10, we can compute

$$ $\dfrac{\partial l_{t_{j+1}}}{\partial z_{0}}=\dfrac{\partial l_{t_{j}}}{ \partial z_{0}}-\Delta t\left(\dfrac{\partial z_{t_{j}}}{\partial z_{0}}\right )^{\top}\nabla_{z}\left(\nabla_{z}\cdot f(t_{j},z_{t_{j}};\theta)\right)$ (110) $$

with additional $\mathcal{O}(d^{2})$ operations provided that $\dfrac{\partial l_{t_{j}}}{\partial z_{0}}$ and $\dfrac{\partial z_{t_{j}}}{\partial z_{0}}$ are stored. Note that computing the inverse of the Jacobian matrix or solving the related linear systems requires $\mathcal{O}(d^{3})$ operations following classical numerical linear algebra methods. Therefore, computing $s_{t_{j+1}}$ through equation 18 requires additional $\mathcal{O}(d^{3})$ operations if we smartly record the results at $t_{j}$. Consequently, the total cost for computing a whole trajectory $\{s_{t_{j}}\}_{j=0}^{N_{t}}$ is $\mathcal{O}(N_{t}\,d^{3})$. We also remark that Huang et al. () give a derivation similar to equation 18.

Other methods to compute the score function. For the KDE method, an approximation for the density function is obtained through taking the kernel convolution with samples directly. Then, one compute its spatial derivative and obtain the score function via $\nabla_{x}\log\rho(t,x)=\nabla_{x}\rho(t,x)/\rho(t,x)$. However, the choice of kernel and bandwidth could be tricky, and the numerical error could be large, especially in high dimensions (Węglarczyk, ). Score matching (Hyvärinen & Dayan, ) is another commonly used method. With a small trick of integration by part and dropping a constant term that includes the integral with squared score functions, one only requires samples to train a network for the score function. However, training a neural network could be expensive, compared with the KDE method, which evaluates the score function directly.

Regularity for the MFC problem.
Formula equation 3 requires third order derivatives for the velocity field $f$. As a consequence, we need regularity assumption for the MFC problem to guarantee that the ODE system is meaningful. In this work, we assume that $L$ and $F$ are $C^{1,2,2}$ Hölder continuous, and $G$ is $C^{2,2}$ Hölder continuous. According to Schauder estimates (Ladyzhenskaia et al., , Chapter 4), these regularity assumptions guarantee that the velocity field has up to third-order derivatives, ensuring that the ODE system equation 3 is meaningful.

The multi-stage splicing method. For our multi-stage splicing method (cf. Appendix C.2), we partition the total time span evenly and pass the trained neural network from the previous stage into the next stage. When the dynamic becomes complicated, the length of the time interval could also be adapted to suit the problem (Zhou et al., ).

Also, the neural network from the previous stage is probably not an accurate velocity for the new stage, but it does contain more information than a random initialization. Therefore, this nice initialization could serve as a warm-up (Zhou et al., ) for the training in the new stage, which could accelerate the convergence of this algorithm.

Different numerical schemes for time discretizations. In the numerical tests, we have also tried other schemes besides the forward Euler scheme, such as the implicit Euler scheme, the mid-point scheme, and the Runge-Kutta scheme. The numerical results for different schemes are similar. We believe that the time discretization error is not the bottleneck of our algorithm at this moment. It is probably more important to train the neural network properly.

## Peer Reviews

confidence: 4
contribution: 1
initial_rating: 1
presentation: 2
questions: Can the authors address the points raised in the **Weaknesses** above?
soundness: 1
strengths: The technical presentation of the material is clear.
summary: The paper gives a system of ODE equations to solve transport equations for a density $\rho(t,x)$ based on the knowledge of the score $s(t,x) = \nabla \log \rho(t,x)$ and the Hessian $H(t,x) = \nabla \nabla \log \rho(t,x)$. These equations are proposed as a way to solve mean field control problems after learning a drift variationally.
weaknesses: The paper basically contains no new material:

1. Eqs. (3) are well-known and have appeared in many places, including 

https://arxiv.org/abs/2206.00860

https://arxiv.org/abs/2210.04296

https://arxiv.org/abs/2206.04642

All these equations can basically derived by the method of characteristics, assuming that the score $s(t,x) = \nabla \log \rho(t,x)$ and the Hessian $H(t,x) = \nabla \nabla \log \rho(t,x)$ are known. The authors should explained better what is the added values of these equations for the numerical scheme they propose. 

2. The material in Sec. 3 is also standard: Eqs. (7) - (10) are the well-known forward Euler discretization of Eq. (3), and Eq. (11) is just their roll-out version. In addition, the neural architecture proposed in Eq. (5) is just a one-hider layer neural network, which is more more simple that standard approximation use for the score (that use UNet, or DiT, etc). Why use such a simplistic architecture?  This seems limitative as well as unnecessary. The authors should explain better what justified this choice.

3. The algorithm proposed does not solve one of the main (and also well-known) issues with neural ODE, namely that they are not simulation free and as a result are to scale to large problem: in particular, differentiating through Eq. (11) is required in Algorithm 1 and costly. In particular, the authors should explain better why they believe that their algorithm may be scalable to high dimensional problems. 

4. The material in Sec. 4 is again standard. In particular Prop. 3 is a well-known set of forward-backward PDE to solve MFC, and Core. 1 is an immediate generalization of this result for a specific choice of Lagrangian and terminal cost. What is the added value of including these results in main text?

5. The numerical examples are too simple to be convincing. In particular, the example in Sec. 5.1 can be factorized over the dimension (which is why it is analytically solvable), which makes it not very challenging computationally.  Similarly, the problem treated in Sec. 5.2 involves a linear OU process, also factorizable. Including these results as test-cases illustration is okay, but the method should also be tested on more complex examples, like the ones found in the cited papers by Lipman et al. 2022 and Boffi & Vanden-Eijnden 2023.

---

confidence: 2
contribution: 3
initial_rating: 5
presentation: 2
questions: 1. What do you consider to be the main contribution / advance of this method? If someone is only interested in first-order score functions, say, is this method still useful?
2. Are there applications other than second-order MFC where the second-order score function is useful?
3. What are the classical methods for computing (second-order) score functions? Did you compare with any of these?
4. What are the state-of-the-art machine learning methods for estimating (second-order) score functions? Did you compare with any of these?
5. In a number of places, the authors briefly mention the applicability of the method to generative modeling. Could you elaborate on how the proposed method would be useful for generative modeling? What advantage would this offer over existing approaches? 
6. Have you tested the proposed method on higher-dimensional systems?
7. Have you timed the proposed method? How does it depend on the dimension of the system? How does it compare to related methods?

While I generally find the paper quite good, I have sufficient doubts about the contribution and the comparisons with related work that I can't recommend it for acceptance right now. However, I would be happy to increase my scores if these doubts are adequately resolved.
soundness: 2
strengths: 1. The paper is generally clear and well written.
2. The proposed method appears novel and may be impactful (but ultimately a lack of detailed comparisons with the existing literature makes it hard to judge, see below).
3. The novel HJB regularization term appears to significantly improve the results.
summary: This paper proposes a neural ODE system, along with its discretization, for computing first- and second-order score functions. As an application, the authors reformulate second-order mean field control (MFC) problems with individual noises into an unconstrained optimization problem using the proposed neural ODE system, and additionally derive a novel regularization term based on viscous Hamilton-Jacobi-Bellman (HJB) equations. In numerical experiments, including regularized Wasserstein proximal operators, probability flow matching, and linear quadratic MFC problems, the authors demonstrate the accuracy of the proposed method and the benefits of the HJB regularization term.
weaknesses: 1. While the paper is generally well written, there are a couple of issues with the presentation. Firstly, the significance of the paper is not clear to me; what, exactly, is the contribution? Please include a paragraph explicitly stating the main contribution of this paper and how it advances the state-of-the-art.

2. More generally, comparisons with the existing literature are not sufficient. There are two types of work that might be relevant here: methods for computing score functions and methods for solving (second-order) MFC problems. The authors list a large number of methods under related work, but none with sufficient detail to understand exactly how they relate to the proposed method. Please include a more detailed discussion of a handful of key works, so the reader can understand how the proposed method compares to the state-of-the-art in terms of approach.

3. Having identified key related works, please include empirical comparisons with the most relevant methods. so the reader can understand how the proposed method compares to the state-of-the-art in terms of performance.

4. In the second paragraph of the introduction, the proposed method is motivated as follows: 

    “While score functions provide powerful tools for modeling stochastic trajectories, their computations are often inefficient, especially in high-dimensional spaces. Classical methods, such as kernel density estimation (KDE) (Chen, 2017), tend to perform poorly in such settings due to the curse of dimensionality (Terrell & Scott, 1992).” 

    However, the experiments are not very high-dimensional (max. 10). As a result, and also due to the lack of empirical comparisons with other methods, it is not clear whether the proposed method actually solves the problem that was identified at the beginning of the paper.

5. Despite computational efficiency being a motivation, the paper does not include timings of the proposed method for any of the experiments. The closest we get is a discussion of the asymptotic complexity for computing the first-order score function. Please include a table or figure showing empirical runtime measurements for the proposed method across different problem dimensions, comparing these to baseline approaches on the same problems.

---

confidence: 2
contribution: 2
initial_rating: 5
presentation: 3
questions: What are some other methods that could be used in this type of problem setting? For example, can neural networks used for modeling differential equations in this paper be replaced by, e.g., Gaussian process regression? Can you show experimentally that the proposed method is indeed superior when compared to such methods?
soundness: 3
strengths: I am not familiar with mean field control problems, but I believe that the method proposed in this paper is a novel approach as far as understood from the survey and other papers cited in this paper. Numerical experiments have confirmed that the proposed method in fact solves mean field control problems accurately and that the regularization term based on the HJB equation is effective.
summary: In this paper, a method for modelling the trajectory of score functions by using neural ordinary differential equations is proposed. Also, its application to solving mean field control problems is shown, along with a regularization term based on the Hamilton-Jacobi-Bellman equation.
weaknesses: Numerical experiments have been performed only with the proposed method, and no comparison with other methods is shown. Therefore, I am not sure that the proposed method is in fact sperior to existing methods. If there are existing methods that can be applied to the problems used in the experiments, the proposed method should be compared with such methods.

In addition, some theoretical results are presented, but most of them are related to the derivation of equations and so on. No theoretical support for the proposed method, such as the generalization error analysis, is presented.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 2
questions: 1) Is the function L in the MFC objective typically assumed to be strongly convex? 
2) How does this compare with existing, much faster generative models/techniques such as Flow Matching (and variants)? I am aware that those are highly efficient without necessarily being higher-order models.
3) Is there work for control problems beyond LQ problems (which are well-studied)?
4) You mentioned " It is probably more important to train the neural network properly." for the error, can you elaborate a bit more on what is meant by this point?
soundness: 3
strengths: 1. The connection of score-based modeling with mean-field control problems is novel and noteworthy.
2. The theory is rigorous and the paper is largely written in a self-contained manner. 
3. Many cases of regularization/structure are shown including the HJB regularizer, Wasserstein proximal operators, etc.
summary: This work deals with extending the (now) mainstream score-based generative modeling techniques to the case of mean-field control problems. The authors construct high-order normalizing flows (i.e. neural ODE systems for the score) and reformulate the
mean field control (MFC) problem with individual noises into an unconstrained optimization problem framed by the proposed neural ODE system. They estimate the first- and second-order score functions using a deep neural network function.
weaknesses: 1) Clarity of writing could be improved. It is at times hard to follow.
2) This is addressed in the Questions section, but I am not sure/convinced about the efficacy of higher-order ODEs in comparison to methods like Flow Matching or Stochastic Interpolants (in general). It would be nice to see some numerical comparisons (if those can be cast for MFC problems as well)  -- I know there is a theoretical section on flow matching but I would like to see tradeoffs gained/lost with this augmentation.

Meta-review: This paper derives a system of ODE equations to solve density transport based the score function and associated Hessian.  A reformulation of second-order mean field control (MFC) problems using the proposed neural ODEs is explored as an application. A strength of the paper is clarity of presentation. Unfortunately, the biggest weakness is lack of novelty and clear contribution. It is known through the Fokker Plank equations how a PDE describes the time evolution of the density transport process. Substantial parts of the paper are known.  Concerns about novelty, scalability challenges and somewhat unrealistic simplicity of the experiments were not resolved satisfactorily during the rebuttal process. Hence, the decision is that the paper does not meet ICLR acceptance bar.

Decision: Reject