Gradient Flow Algorithms for Density Propagation in Stochastic Systems

TL;DR

This paper introduces a novel gradient flow-based computational framework for propagating joint probability density functions in stochastic nonlinear systems, avoiding traditional discretization limitations and enabling fast, non-parametric uncertainty analysis.

Contribution

It develops a time-discretized, geometry-inspired method that solves infinite-dimensional PDEs via proximal recursions on the space of PDFs, with guaranteed convergence and efficiency.

Findings

Achieves fast, non-parametric joint PDF propagation.

Proves contractiveness of the fixed point recursion.

Demonstrates effectiveness through numerical examples.

Abstract

We develop a new computational framework to solve the partial differential equations (PDEs) governing the flow of the joint probability density functions (PDFs) in continuous-time stochastic nonlinear systems. The need for computing the transient joint PDFs subject to prior dynamics arises in uncertainty propagation, nonlinear filtering and stochastic control. Our methodology breaks away from the traditional approach of spatial discretization or function approximation -- both of which, in general, suffer from the "curse-of-dimensionality". In the proposed framework, we discretize time but not the state space. We solve infinite dimensional proximal recursions in the manifold of joint PDFs, which in the small time-step limit, is theoretically equivalent to solving the underlying transport PDEs. The resulting computation has the geometric interpretation of gradient flow of certain free energy functional with respect to the Wasserstein metric arising from the theory of optimal mass transport. We show that dualization along with an entropic regularization, leads to a cone-preserving fixed point recursion that is proved to be contractive in Thompson metric. A block co-ordinate iteration scheme is proposed to solve the resulting nonlinear recursions with guaranteed convergence. This approach enables remarkably fast computation for non-parametric transient joint PDF propagation. Numerical examples and various extensions are provided to illustrate the scope and efficacy of the proposed approach.

Equations

$$\mathrm{d}\mathbf{x}=\mathbf{f}\left(\mathbf{x},t\right)\mathrm{d}t+\mathbf{g}(\mathbf{x},t)\mathrm{d}\mathbf{w},\mathbf{x}(t=0)={\mathbf{x}}_0,$$

$$\frac{\partial \rho }{\partial t}=-\nabla \cdot \left(\rho \mathbf{f}\right)+\frac{1}{2}\sum _{i,j=1}^n\frac{{\partial }^2}{\partial x_i\partial x_j}(\rho \mathbf{g}{\mathbf{g}}^{\top }{)}_{ij},$$

$$\rho (\mathbf{x},t=0)={\rho }_0(\mathbf{x})\text{(given)}.$$

$$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=-\nabla \phi \left(\mathbf{x}\right),\mathbf{x}(0)={\mathbf{x}}_0,$$

$${\mathbf{x}}_k=\underset{\mathbf{x}\in {\mathbb{R}}^n}{\mathrm{argmin}}\frac{1}{2}\parallel \mathbf{x}-{\mathbf{x}}_{k-1}{\parallel }_2^2+h\phi (\mathbf{x})+o(h),k\in \mathbb{N},$$

$${\mathbf{x}}_k-{\mathbf{x}}_{k-1}=-h\nabla \phi ({\mathbf{x}}_{k-1}),$$

$${\mathbf{x}}_k$$

$$=\underset{\mathbf{x}\in {\mathbb{R}}^n}{\mathrm{argmin}}\frac{1}{2}\parallel \mathbf{x}-{\mathbf{x}}_{k-1}{\parallel }_2^2+⟨\mathbf{x}-{\mathbf{x}}_{k-1},h\nabla \phi ({\mathbf{x}}_{k-1})⟩$$

$$+h\phi ({\mathbf{x}}_{k-1}),$$

$${\mathrm{prox}}_{h\phi }^{\parallel \cdot {\parallel }_2}({\mathbf{x}}_{k-1}):=\underset{\mathbf{x}\in {\mathbb{R}}^n}{\mathrm{argmin}}\frac{1}{2}\parallel \mathbf{x}-{\mathbf{x}}_{k-1}{\parallel }_2^2+h\phi (\mathbf{x}),$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\phi =⟨\nabla \phi ,-\nabla \phi ⟩=-\parallel \nabla \phi {\parallel }_2^2$$

$${\mathcal{D}}_2:=\{\rho :{\mathbb{R}}^n\mapsto \mathbb{R}\mid \rho \ge 0,{\int }_{{\mathbb{R}}^n}\rho =1,{\mathbb{E}}_{\rho }[{\mathbf{x}}^{\top }\mathbf{x}]<\infty \}.$$

$${\mathrm{prox}}_{h\Phi }^{d(\cdot )}({\varrho }_{k-1}):=\underset{\varrho \in {\mathcal{D}}_2}{\mathrm{arginf}}\frac{1}{2}d{\left(\varrho ,{\varrho }_{k-1}\right)}^2+h\Phi (\varrho ),$$

$${\varrho }_k={\mathrm{prox}}_{h\Phi }^{d(\cdot )}({\varrho }_{k-1}),k\in \mathbb{N},{\varrho }_0(\mathbf{x}):={\rho }_0(\mathbf{x}).$$

$$W({\pi }_1,{\pi }_2):={\left(\underset{\mathrm{d}\pi \in \Pi \left({\pi }_1,{\pi }_2\right)}{\inf }{\int }_{\mathcal{X}\times \mathcal{Y}}s(\mathbf{x};\mathbf{y})\mathrm{d}\pi \left(\mathbf{x},\mathbf{y}\right)\right)}^{\frac{1}{2}},$$

$$W^2(\nu ,(1-\tau ){\rho }_1+\tau {\rho }_2)\le (1-\tau )W^2(\nu ,{\rho }_1)+\tau W^2(\nu ,{\rho }_2).$$

$$D_{\mathrm{KL}}({\pi }_1||{\pi }_2)$$

$$={\int }_{{\mathbb{R}}^n}{\rho }_1(\mathbf{x})\log \frac{{\rho }_1(\mathbf{x})}{{\rho }_2(\mathbf{x})}\mathrm{d}\mathbf{x},$$

$$\mathrm{d}\mathbf{x}=-\nabla \psi \left(\mathbf{x}\right)\mathrm{d}t+\sqrt{2{\beta }^{-1}}\mathrm{d}\mathbf{w},\mathbf{x}(0)={\mathbf{x}}_0,$$

$$\frac{\partial \rho }{\partial t}=\nabla \cdot \left(\rho \nabla \psi \right)+{\beta }^{-1}\Delta \rho ,\rho (\mathbf{x},0)={\rho }_0(\mathbf{x}).$$

$$F(\rho )$$

$$\frac{\partial \rho }{\partial t}=\nabla \cdot \left(\rho \nabla \zeta \right),\text{where}\zeta :={\beta }^{-1}\left(1+\log \rho \right)+\psi ,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}F=-{\mathbb{E}}_{\rho }\left[\parallel \nabla \zeta {\parallel }_2^2\right],$$

$$\frac{\partial \rho }{\partial t}=\nabla \cdot \left(\rho \nabla \left(\psi +\rho \ast \varphi \right)\right)+{\beta }^{-1}\Delta \rho ,\rho (\mathbf{x},0)={\rho }_0(\mathbf{x}),$$

$$\mathrm{d}\mathbf{x}=-\left(\nabla \psi \left(\mathbf{x}\right)+\nabla \left(\rho \ast \varphi \right)\right)\mathrm{d}t+\sqrt{2{\beta }^{-1}}\mathrm{d}\mathbf{w},\mathbf{x}(0)={\mathbf{x}}_0.$$

$$F(\rho ):={\mathbb{E}}_{\rho }\left[\psi +{\beta }^{-1}\log \rho +\rho \ast \varphi \right],$$

$$\zeta :={\beta }^{-1}\left(1+\log \rho \right)+\psi +\rho \ast \varphi ,$$

$${\varrho }_k$$

$$=\underset{\varrho \in {\mathcal{D}}_2}{\mathrm{arginf}}\frac{1}{2}{W}^2({\varrho }_{k-1},\varrho )+hF(\varrho ),k\in \mathbb{N},$$

$$\displaystyle{\bm{\varrho}}_{k}=\underset{{\bm{\varrho}}}{\operatorname*{arg\>min}}\bigg{\{}\underset{\bm{M}\in\Pi({\bm{\varrho}}_{k-1},{\bm{\varrho}})}{\min}\frac{1}{2}\langle\bm{C}_{k},\bm{M}\rangle+h\>\langle\bm{\psi}_{k-1}$$

$$\displaystyle+\beta^{-1}\log{\bm{\varrho}},{\bm{\varrho}}\rangle\bigg{\}},\quad k\in\mathbb{N},$$

$${\mathbf{\psi }}_{k-1}(i):=\psi \left({\mathbf{x}}_{k-1}^i\right),i=1,\dots ,N.$$

$${\mathbf{C}}_k(i,j):=s({\mathbf{x}}_k^i;{\mathbf{x}}_{k-1}^j)=\parallel {\mathbf{x}}_k^i-{\mathbf{x}}_{k-1}^j{\parallel }_2^2,i,j=1,\dots ,N,$$

Full text

Keywords: Proximal operator, Fokker-Planck-Kolmogorov equation, optimal transport, gradient descent, uncertainty propagation.

I Introduction

Consider the continuous-time dynamics of the state vector $\bm{x}(t)\in\mathbb{R}^{n}$ governed by an Itô stochastic differential equation (SDE)

[TABLE]

where the joint probability density function (PDF) for the initial condition $\bm{x}_{0}$ is a known function $\rho_{0}$; we use the notation $\bm{x}_{0}\sim\rho_{0}$. The process noise ${\bm{w}}(t)\in\mathbb{R}^{m}$ is Wiener and satisfy $\mathbb{E}\left[{\rm{d}}w_{i}{\rm{d}}w_{j}\right]=\delta_{ij}{\rm{d}}t$ for all $i,j=1,\ldots,n$, where $\delta_{ij}=1$ for $i=j$, and zero otherwise. Then, the flow of the joint PDF $\rho(\bm{x},t)$ for the state vector $\bm{x}(t)$ (i.e., $\bm{x}\sim\rho$) is governed by the partial differential equation (PDE) initial value problem:

[TABLE]

The second order transport PDE (2a) is known as the Fokker-Planck or Kolmogorov’s forward equation [1]. Hereafter, we will refer it as the Fokker-Planck-Kolmogorov (FPK) PDE.

In this paper, we consider the problem of density or belief propagation, i.e., the problem of computing the transient joint PDF $\rho(\bm{x},t)$ that solves a PDE initial value problem of the form (2). From an application standpoint, the need for computing $\rho(\bm{x},t)$ can be motivated by two types of problems. The first is dispersion analysis, where one is interested in predicting or analyzing the uncertainty evolution over time, e.g., in meteorological forecasting [2], spacecraft entry-descent-landing [3, 4], orientation density evolution for liquid crystals in chemical physics [5, 6, 7], and in motion planning [8, 9, 10]. In these applications, the quantity of interest is the joint PDF $\rho(\bm{x},t)$ and its statistics (e.g., transient moments and marginal PDFs). The second type of applications require $\rho(\bm{x},t)$ as an intermediate step toward computing other quantities of interest. For example, in nonlinear filtering [11, 12], the joint PDF $\rho(\bm{x},t)$ serves as the prior in computing the posterior (i.e., conditional state) PDF. In probabilistic model validation [13, 14, 15] and controller verification [16], computing $\rho(\bm{x},t)$ helps in quantifying the density-level prediction-observation mismatch. All these applications require fast computation of $\rho(\bm{x},t)$ in a scalable and unified manner, as opposed to developing algorithms in a case-by-case basis.

Given its widespread applications, problem (2) has received sustained attention from the scientific computing community where the predominant solution approaches have been spatial discretization and function approximation – both of which, in general, suffer from the “curse of dimensionality” [17]. The purpose of this paper is to pursue a systems-theoretic variational viewpoint for computing $\rho(\bm{x},t)$ that breaks away from the “solve PDE as a PDE” philosophy, and instead solves (2) as a gradient descent on the manifold of joint PDFs. This emerging geometric viewpoint for uncertainty propagation and filtering has been reported in our recent work [18, 19], but it remained unclear whether this viewpoint can offer computational benefit over the standard PDE solvers. It is not at all obvious whether and how an infinite-dimensional gradient descent can numerically obviate function approximation or spatial discretization. The contribution of this paper is to demonstrate that not only this is possible, but also that the same enables fast computation.

To conceptualize the main idea, we appeal to the metric viewpoint of gradient descent, where a continuous-time gradient flow is realized by small time-step recursions of a proximal operator with respect to (w.r.t.) a suitable metric. For example, consider the finite dimensional gradient flow

[TABLE]

where ${\bm{x}},{\bm{x}}_{0}\in\mathbb{R}^{n}$, and the continuously differentiable function $\varphi:\mathbb{R}^{n}\rightarrow\mathbb{R}_{\geq 0}$. The flow ${\bm{x}}(t)$ generated by (3) can be realized via variational recursion

[TABLE]

in the sense that as the step-size $h\downarrow 0$, we have $\bm{x}_{k}\rightarrow\bm{x}(t=kh)$, i.e., the sequence $\{{\bm{x}}_{k}\}$ converges pointwise to the flow ${\bm{x}}(t)$. This can be verified by rewriting the Euler discretization of (3), given by

[TABLE]

as

[TABLE]

where we used the fact that adding and omitting constant terms do not change the arg min. From (5), one can arrive at (4) by invoking first order approximation of $\varphi({\bm{x}})$ at ${\bm{x}}={\bm{x}}_{k}$. In (4), the mapping ${\bm{x}}_{k-1}\mapsto{\bm{x}}_{k}$ given by

[TABLE]

is called the “proximal operator” [20, p. 142] of $h\varphi$ w.r.t. the standard Euclidean metric $\parallel\cdot\parallel_{2}$. Notice that $\varphi(\cdot)$ serves as a Lyapunov function since the quantity

[TABLE]

equals [math] at the stationary point of (3), and $<0$ otherwise.

In the infinite dimensional setting, we are interested in computing the flow generated by (2) via gradient descent on the manifold of joint PDFs with finite second (raw) moments, denoted as***We denote the expectation operator w.r.t. the measure $\rho({\bm{x}}){\rm{d}}{\bm{x}}$ as $\mathbb{E}_{\rho}\left[\cdot\right]$.

[TABLE]

Specifically, let $d(\cdot,\cdot)$ be a distance metric on the manifold $\mathscr{D}_{2}$, and let the functional $\Phi:\mathscr{D}_{2}\mapsto\mathbb{R}_{\geq 0}$. Then, for some chosen step-size $h>0$, the infinite dimensional proximal operator of $h\Phi$ w.r.t. the distance metric $d(\cdot,\cdot)$, given by

[TABLE]

can be used to define a proximal recursion (Fig. 1):

[TABLE]

Just as the proximal recursion (4) approximates the finite dimensional flow (3), similarly it is possible to design $d(\cdot,\cdot)$ and $\Phi(\cdot)$ in (8) as function of the drift $\bm{f}$ and diffusion $\bm{g}$ in (2a) such that the proximal recursion (9) approximates the infinite dimensional flow (2). This idea was first proposed in [21], showing that when $\bm{f}$ is a gradient vector field and $\bm{g}$ is a scalar multiple of identity matrix, then the distance $d(\cdot,\cdot)$ can be taken as the Wasserstein-2 metric [22] with $\Phi(\cdot)$ as the free energy functional. In particular, the solution of (9) was shown to converge to the flow of (2), i.e., $\varrho_{k}({\bm{x}})\rightarrow\rho({\bm{x}},t=kh)$ in strong $L_{1}(\mathbb{R}^{n})$ sense, as $h\downarrow 0$. The resulting variational recursion (9) has since been known as the Jordan-Kinderlehrer-Otto (JKO) scheme [23], and we will refer to the FPK operator (2a) with such assumptions on $\bm{f}$ and $\bm{g}$ to be in “JKO canonical form”. Similar gradient descent schemes have been derived for many other PDEs; see e.g., [24] for a recent survey.

The remaining of this paper is organized as follows. Section II explains the JKO canonical form and the corresponding free energy functional $\Phi(\cdot)$. Our algorithms and convergence results are presented in Section III, followed by numerical simulation results in Section IV. Section V provides various extensions of the basic algorithm showing how the framework proposed here can be applied to systems not in JKO canonical form. Section VI concludes the paper.

We remark here that a preliminary version [25] of this work appeared in the 2019 American Control Conference. This paper significantly expands [25] by incorporating additional results for the so-called McKean-Vlasov flow (Sections II-B and IV-C), and by providing various extensions of the basic algorithm (Section V).

I-1 Notations

Throughout the paper, we use bold-faced capital letters for matrices and bold-faced lower-case letters for column vectors. We use the symbol $\langle\cdot,\cdot\rangle$ to denote the Euclidean inner product. In particular, $\langle\bm{A},\bm{B}\rangle:={\mathrm{trace}}({\bm{A}}^{\top}\bm{B})$ denotes Frobenius inner product between matrices $\bm{A}$ and $\bm{B}$, and $\langle\bm{a},\bm{b}\rangle:=\bm{a}^{\top}\bm{b}$ denotes the inner product between column vectors $\bm{a}$ and $\bm{b}$. The symbols $\nabla$ and $\Delta$ denote the (Euclidean) gradient and the Laplacian operators, respectively. We use $\mathcal{N}(\mu,\sigma^{2})$ to denote a univariate Gaussian PDF with mean $\mu$ and variance $\sigma^{2}$. Likewise, $\mathcal{N}(\bm{\mu},\bm{\Sigma})$ denotes a multivariate Gaussian PDF with mean vector $\bm{\mu}$ and covariance matrix $\bm{\Sigma}$. The expectation operator w.r.t. the PDF $\rho$ is denoted by $\mathbb{E}_{\rho}\left[\cdot\right]$, i.e., $\mathbb{E}_{\rho}\left[\cdot\right]:=\int_{\mathbb{R}^{n}}(\cdot)\rho\>{\rm{d}}{\bm{x}}$. The operands $\log(\cdot)$, $\exp(\cdot)$ and $\geq 0$ are to be understood as element-wise. The notations $\odot$ and $\oslash$ denote element-wise (Hadamard) product and division, respectively. We use $\bm{I}_{n}$ to denote the $n\times n$ identity matrix. The symbols $\bm{1}$ and $\bm{0}$ stand for column vectors of appropriate dimension containing all ones, and all zeroes, respectively.

I-2 Preliminaries

We next collect definitions and some properties of the 2-Wasserstein metric and the Kullback-Leibler divergence, which will be useful in the development below.

Definition 1**.**

(2-Wasserstein metric) The 2-Wasserstein metric between two probability measures ${\rm{d}}\pi_{1}({\bm{x}}):=\rho_{1}({\bm{x}}){\rm{d}}{\bm{x}}$ and ${\rm{d}}\pi_{2}(\bm{y}):=\rho_{2}(\bm{y}){\rm{d}}\bm{y}$, supported respectively on $\mathcal{X},\mathcal{Y}\subseteq\mathbb{R}^{n}$, is denoted as $W(\pi_{1},\pi_{2})$ (equivalently, $W(\rho_{1},\rho_{2})$ whenever the measures $\pi_{1},\pi_{2}$ are absolutely continuous so that the respective PDFs $\rho_{1},\rho_{2}$ exist); it is defined as

[TABLE]

where $s(\bm{x};\bm{y}):=\parallel\bm{x}-\bm{y}\parallel_{2}^{2}$ is the squared Euclidean distance in $\mathbb{R}^{n}$, and $\Pi\left(\pi_{1},\pi_{2}\right)$ denotes the collection of all joint probability measures on $\mathcal{X}\times\mathcal{Y}$ having finite second moments, with marginals $\pi_{1}$ and $\pi_{2}$, respectively.

The existence and uniqueness of the minimizer in (10) is guaranteed. It is well-known [22, Ch. 7] that $W(\pi_{1},\pi_{2})$ defines a metric on the manifold $\mathscr{D}_{2}$. This means that $W(\pi_{1},\pi_{2})\geq 0$ with equality if and only if $\pi_{1}=\pi_{2}$, the symmetry: $W(\pi_{1},\pi_{2})=W(\pi_{2},\pi_{1})$, and that $W(\pi_{1},\pi_{2})$ satisfies the triangle inequality. Its square, $W^{2}(\pi_{1},\pi_{2})$ equals [26] the minimum amount of work required to transport $\pi_{1}$ to $\pi_{2}$ (or equivalently, $\rho_{1}$ to $\rho_{2}$), and vice versa. For any PDF $\nu$, the function $\rho\mapsto W^{2}(\rho,\nu)$ is convex on $\mathscr{D}_{2}$, i.e., for any $\rho_{1},\rho_{2}\in\mathscr{D}_{2}$, and $0\leq\tau\leq 1$, we have

[TABLE]

Definition 2**.**

(Kullback-Leibler divergence) The Kullback-Leibler divergence, also known as relative entropy, between two probability measures ${\rm{d}}\pi_{1}({\bm{x}}):=\rho_{1}({\bm{x}}){\rm{d}}{\bm{x}}$ and ${\rm{d}}\pi_{2}(\bm{y}):=\rho_{2}(\bm{y}){\rm{d}}\bm{y}$, denoted as $D_{\rm{KL}}(\pi_{1}||\pi_{2})$, is defined as

[TABLE]

where $\frac{{\rm{d}}\pi_{1}}{{\rm{d}}\pi_{2}}$ denotes the Radon-Nikodym derivative. Henceforth, we use the notational equivalence $D_{\rm{KL}}(\pi_{1}||\pi_{2})\equiv D_{\rm{KL}}(\rho_{1}||\rho_{2})$.

From Jensen’s inequality, $D_{\rm{KL}}(\rho_{1}||\rho_{2})\geq 0$; however, $D_{\rm{KL}}$ is not a metric since it is neither symmetric, nor does it satisfy the triangle inequality. The Kullback-Leibler divergence (12) is jointly convex in $\rho_{1}$ and $\rho_{2}$.

II JKO Canonical Form

The JKO canonical form refers to a continuous-time stochastic dynamics where the drift vector field is the gradient of a potential function. The potential can be state dependent, or can depend on both the state ${\bm{x}}(t)$ and its joint PDF $\rho({\bm{x}},t)$. In the former case, the associated flow of the joint state PDF is governed by a FPK PDE of the form (2a) whereas in the latter case, the same is governed by the McKean-Vlasov integro-PDE. For both cases, the proposed gradient flow algorithms (Section III) will be able to compute the transient state PDFs $\rho({\bm{x}},t)$. We next describe these canonical forms; generalization of our framework to systems not in JKO canonical form will be given in Sections IV and V.

II-A FPK Gradient Flow

Consider the Itô SDE

[TABLE]

where the the drift potential $\psi:\mathbb{R}^{n}\mapsto(0,\infty)$, the diffusion coefficient $\beta>0$, and the initial condition ${\bm{x}}_{0}\sim\rho_{0}({\bm{x}})$. For the sample path ${\bm{x}}(t)$ dynamics given by the SDE (13), the flow of the joint PDF $\rho\left(\bm{x},t\right)$ is governed by the FPK PDE initial value problem

[TABLE]

It is easy to verify that the unique stationary solution of (14) is the Gibbs PDF $\rho_{\infty}(\bm{x})=\kappa\exp\left(-\beta\psi(\bm{x})\right)$, where the normalizing constant $\kappa$ is known as the partition function.

A Lyapunov functional associated with (14) is the free energy

[TABLE]

that decays [21] along the solution trajectory of (14). This follows from re-writing (14) as

[TABLE]

and consequently

[TABLE]

which is $<0$ for the transient solution $\rho({\bm{x}},t)$, and $=0$ at the stationary solution $\rho_{\infty}({\bm{x}})=\kappa\exp(-\beta\psi(\bm{x}))$. In our context, (17) is an analogue of (7). Notice that the free energy (15) can be seen as the sum of the potential energy $\int_{\mathbb{R}^{n}}\psi({\bm{x}})\rho\>{\rm{d}}{\bm{x}}$ and the internal energy $\beta^{-1}\int_{\mathbb{R}^{n}}\rho\log\rho\>{\rm{d}}{\bm{x}}$. If $\psi\equiv 0$, the PDE (14) reduces to the heat equation, which by (15), can then be interpreted as an entropy maximizing flow.

The seminal result of [21] was that the transient solution of (14) can be computed via the proximal recursion (9) with the distance metric $d\equiv W$ (i.e., the 2-Wasserstein metric in (10)) and the functional $\Phi\equiv F$ (i.e., the free energy (15)). Just as (3) and (4) form a gradient flow-proximal recursion pair, likewise (14) and (9) form the same with the stated choices of $d$ and $\Phi$. From (10), notice that the distance metric $W$ itself is defined as the solution of an optimization problem, hence it is not apparent how to numerically implement the recursion (9) in a scalable manner.

II-B McKean-Vlasov Gradient Flow

In addition to the drift and diffusion, if one has non-local interaction, then the corresponding PDF evolution equation becomes the McKean-Vlasov integro-PDE

[TABLE]

where $*$ denotes the convolution in $\mathbb{R}^{n}$, the interaction potential $\phi:\mathbb{R}^{n}\mapsto(0,\infty)$ and is symmetric, i.e., $\phi(-\bm{v})=\phi(\bm{v})$ for $\bm{v}\in\mathbb{R}^{n}$. The associated sample path $\bm{x}(t)$ dynamics has PDF-dependent drift:

[TABLE]

As an example, when $\phi(\bm{v}):=\frac{1}{2}\parallel\bm{v}\parallel^{2}$, then $\nabla(\rho*\phi)({\bm{x}})=\bm{x}-\mathbb{E}_{\rho}[\bm{x}]$. Clearly, (18) reduces to (14) in the absence of interaction ($\phi\equiv 0$). The McKean-Vlasov equation serves as a model for coupled multi-agent interaction in applications such as crowd movement [27], opinion dynamics [28, 29], population biology and communication systems.

A Lyapunov functional for (18) can be obtained [30] by generalizing the free energy (15) as

[TABLE]

which is a sum of the potential energy (as before), the internal energy (as before), and the interaction energy $\frac{1}{2}\int_{\mathbb{R}^{2n}}\phi({\bm{x}}-{\bm{y}})\rho({\bm{x}})\rho({\bm{y}})\>{\rm{d}}{\bm{x}}{\rm{d}}{\bm{y}}$. In this case, (17) holds with

[TABLE]

which is often referred to as the entropy dissipation functional [31]. As before, the proximal recursion (9) with $d\equiv W$ and $\Phi\equiv F$ (now $F$ given by (20)), approximates the flow (18).

III Framework

We now describe our computational framework to solve the proximal recursion

[TABLE]

with $\varrho_{0}\equiv\rho_{0}({\bm{x}})$ (the initial joint PDF) for some small step-size $h$. For maximal clarity, we develop the framework with the free energy $F(\cdot)$ as in (15), i.e., for FPK gradient flow. In Section IV.C, we will show how the same can be generalized when $F(\cdot)$ is of the form (20). As per the convexity properties mentioned in Section I.2, problem (22) involves (recursively) minimizing sum of two convex functionals, and hence is a convex problem for each $k\in\mathbb{N}$.

We discretize time as $t=0,h,2h,\ldots$, and develop an algorithm to solve (22) without making any spatial discretization. Specifically, we would like to perform the recursion (22) on weighted scattered point cloud $\{{\bm{x}}_{k}^{i},\varrho_{k}^{i}\}_{i=1}^{N}$ of cardinality $N$ at $t_{k}=kh$, $k\in\mathbb{N}$, where the location of the $i$th point ${\bm{x}}_{k}^{i}\in\mathbb{R}^{n}$ denotes its state-space coordinate, and the corresponding weight $\varrho_{k}^{i}\in\mathbb{R}_{\geq 0}$ denotes the value of the joint PDF evaluated at that point at time $t_{k}$. Such weighted point cloud representation of (22) results in the following problem:

[TABLE]

where the drift potential vector $\bm{\psi}_{k-1}\in\mathbb{R}^{N}$ is given by

[TABLE]

Here, the probability vectors ${\bm{\varrho}},{\bm{\varrho}}_{k-1}\in\mathcal{S}_{N-1}$, the probability simplex in $\mathbb{R}^{N}$. Furthermore, for each $k\in\mathbb{N}$, the matrix $\bm{C}_{k}\in\mathbb{R}^{N\times N}$ is given by

[TABLE]

and $\Pi({\bm{\varrho}}_{k-1},{\bm{\varrho}})$ stands for the set of all matrices $\bm{M}\in\mathbb{R}^{N\times N}$ such that

[TABLE]

Notice that the inner minimization in (23) is a standard linear programming problem if it were to be solved for a given ${\bm{\varrho}}\in\mathcal{S}_{N-1}$, as in the Monge-Kantorovich optimal mass transport [22]. However, the outer minimization in (23) precludes a direct numerical approach.

To circumvent the aforesaid issues, following [32], we first regularize and then dualize (23). Specifically, adding an entropic regularization $H(\bm{M}):=\langle\bm{M},\log\bm{M}\rangle$ in (23) yields

[TABLE]

where $\epsilon>0$ is a regularization parameter. The entropic regularization is standard in optimal mass transport literature [33, 34] and leads to efficient Sinkhorn iteration for the inner minimization. Here we point out that there has been parallel work [35, 36] on the connection between optimal mass transport and the so called Schrödinger bridge problem which is a dynamic version of this type of regularization.

In our context, the entropic regularization “algebrizes” the inner minimization in the sense if $\bm{\lambda}_{0},\bm{\lambda}_{1}\in\mathbb{R}^{N}$ are Lagrange multipliers associated with the equality constraints in (26), then the optimal coupling matrix $\bm{M}^{\rm{opt}}:=[m^{\rm{opt}}(i,j)]_{i,j=1}^{N}$ in (27) has the Sinkhorn form

[TABLE]

Since the objective in (27) is proper convex and lower semi-continuous in ${\bm{\varrho}}$, the strong duality holds, and we consider the Lagrange dual of (27) given by:

[TABLE]

where

[TABLE]

is the Legendre-Fenchel conjugate of the free energy $F$ in (15). Next, we derive the first order optimality conditions for (29) resulting in the proximal updates, and then provide an algorithm to solve the same.

III-A Proximal Recursions

Given the vectors ${\bm{\varrho}}_{k-1},\bm{\psi}_{k-1}$, the matrix $\bm{C}_{k}$, and the positive scalars $\beta,h,\epsilon$ in (29), let

[TABLE]

The following result establishes a system of nonlinear equations for computing $\bm{\lambda}_{0}^{\rm{opt}},\bm{\lambda}_{1}^{\rm{opt}}$ in (29), and consequently ${\bm{\varrho}}_{k}$ in (27).

Theorem 1**.**

The vectors $\bm{\lambda}_{0}^{\rm{opt}},\bm{\lambda}_{1}^{\rm{opt}}$ in (29) can be found by solving for $\bm{y},\bm{z}\in\mathbb{R}^{N}$ from the following system of equations:

[TABLE]

and then inverting the maps (31). Let $\left(\bm{y}^{\rm{opt}},\bm{z}^{\rm{opt}}\right)$ be the solution of (33). The vector ${\bm{\varrho}}_{k}$ in (27), i.e., the proximal update (Fig. 1) can then be obtained as

[TABLE]

Proof.

From (15), the “discrete free energy” is

[TABLE]

Its Legendre-Fenchel conjugate, by (30), is given by

[TABLE]

Setting the gradient of the objective function in $(\ref{DiscreteLF})$ w.r.t. ${\bm{\varrho}}$ to zero, and solving for ${\bm{\varrho}}$ yields

[TABLE]

Substituting (36) back into (35), results

[TABLE]

Fixing $\bm{\lambda}_{1}$, and taking the gradient of the objective in (29) w.r.t. $\bm{\lambda}_{0}$, gives (33a). Likewise, fixing $\bm{\lambda}_{0}$, and taking the gradient of the objective in (29) w.r.t. $\bm{\lambda}_{1}$ gives

[TABLE]

Using (37) to simplify the left-hand-side of (38) results in (33b). To derive (34), notice that combining the last equality constraint in (26) with (28), (31) and (32) gives

[TABLE]

which is equal to $\bm{z}\odot\bm{\Gamma}_{k}^{\top}\bm{y}$, as claimed. ∎

In the following (Section III-B), we propose an algorithm to solve (33) and (34), and then outline the overall implementation of our computational framework. The convergence results for the proposed algorithm are given in Section III-C.

III-B Algorithm

III-B1 Proximal algorithm

We now propose a block co-ordinate iteration scheme to solve (33). The proposed procedure, which we call ProxRecur, and detail in Algortihm 1, takes ${\bm{\varrho}}_{k-1}$ as input and returns the proximal update ${\bm{\varrho}}_{k}$ as output for $k\in\mathbb{N}$. In addition to the data ${\bm{\varrho}}_{k-1},\bm{\psi}_{k-1},\bm{C}_{k},\beta,h,\epsilon,N$, the Algorithm 1 requires two parameters as user input: numerical tolerance $\delta$, and maximum number of iterations $L$. The computation in Algorithm 1, as presented, involves making an initial guess for the vector $\bm{z}$ and then updating $\bm{y}$ and $\bm{z}$ until convergence. The iteration over index $\ell\leq L$ is performed while keeping the physical time “frozen”.

We next outline the overall algorithmic setup to implement the proximal recursion over probability weighted scattered data.

III-B2 Overall algorithm

Samples from the known initial joint PDF $\rho_{0}$ are generated as point cloud $\{{\bm{x}}_{0}^{i},\varrho_{0}^{i}\}_{i=1}^{N}$. Then for $k\in\mathbb{N}$, the point clouds $\{{\bm{x}}_{k}^{i},\varrho_{k}^{i}\}_{i=1}^{N}$ are updated as shown in Fig. 2. Specifically, the state vectors are updated via Euler-Maruyama scheme applied to the underlying SDE (13) or (19). Explicitly, the Euler-Maruyama update for (19) is

[TABLE]

where $\omega(\cdot):=\int\phi(\cdot-\bm{y})\rho(\bm{y}){\rm{d}}\bm{y}$, and $\bm{w}_{k}^{i}:=\bm{w}^{i}(t=kh)$, $k\in\mathbb{N}$. The same for (13) is obtained by setting $\phi\equiv\omega\equiv 0$ in (39).

The corresponding probability weights $\{{\bm{\varrho}}_{k}^{i}\}_{i=1}^{N}$ are updated via Algorithm 1. Notice that computing $\bm{C}_{k}$ requires both $\{\bm{x}_{k-1}^{i}\}_{i=1}^{N}$ and $\{\bm{x}_{k}^{i}\}_{i=1}^{N}$, and that $\bm{C}_{k}$ needs to be passed as input to Algorithm 1. Thus, the execution of Euler-Maruyama scheme precedes that of Algorithm 1.

Remark 1**.**

Our choice of the (explicit) Euler-Maruyama scheme for updating the location of the points in state space was motivated by its simplicity and ease of implementation. Since the diffusion coefficient of (13) or (19) is constant, the Euler-Maruyama scheme is guaranteed to converge strongly to the true solution of the corresponding SDE provided the drift coefficient is globally Lipschitz; see e.g., [37, Ch. 10.2]. If the drift coefficient $-\nabla\psi$ in (13), or $-\left(\nabla\psi+\nabla\left(\rho*\phi\right)\right)$ in (19), is non-globally Lipschitz with superlinear growth, then the Euler-Maruyama scheme is known [38] to diverge in mean-squared error sense. In such cases, one could replace the explicit Euler-Maruyama scheme with the “tamed” Euler-Maruyama scheme [39], or with the partially implicit Euler scheme [40], or with the split-step backward Euler scheme [41].

III-C Convergence

Next, we will prove the convergence properties for Algorithm 1. To this end, the Definition 3 and Proposition 2 given below will be useful in establishing Theorem 3 that follows.

Definition 3**.**

*(Thompson metric)

Consider $\bm{z},\widetilde{\bm{z}}\in\mathcal{K}$, where $\mathcal{K}$ is a non-empty open convex cone. Further, suppose that $\mathcal{K}$ is a normal cone, i.e., there exists constant $\alpha$ such that $\parallel\bm{z}\parallel\leq\alpha\parallel\widetilde{\bm{z}}\parallel$ for $\bm{z}\leq\widetilde{\bm{z}}$. Thompson [42] proved that $\mathcal{K}$ is a complete metric space w.r.t. the so-called Thompson metric given by*

[TABLE]

where $\gamma(\bm{z}/\widetilde{\bm{z}}):=\inf\{c>0\mid\bm{z}\leq c\widetilde{\bm{z}}\}$. In particular, if $\mathcal{K}\equiv\mathbb{R}^{N}_{>0}$ (positive orthant of $\mathbb{R}^{N}$), then

[TABLE]

Proposition 2**.**

[43, Proposition 3.2]**,[44] Let $\mathcal{K}$ be an open, normal, convex cone, and let $\bm{p}:\mathcal{K}\mapsto\mathcal{K}$ be an order preserving homogeneous map of degree $r\geq 0$, i.e., $\bm{p}(c{\bm{z}})=c^{r}\bm{p}(\bm{z})$ for any $c>0$ and $\bm{z}\in\mathcal{K}$. Then, for all $\bm{z},\widetilde{\bm{z}}\in\mathcal{K}$, we have

[TABLE]

In particular, if $r\in[0,1)$, then the map $\bm{p}(\cdot)$ is strictly contractive in the Thompson metric $d_{\rm{T}}$, and admits unique fixed point in $\mathcal{K}$.

Using (40) and Proposition 2, we establish the convergence result below.

Theorem 3**.**

Consider the notations in (31)-(32), and those in Algorithm 1. The iteration

[TABLE]

for $\ell=1,2,\ldots$, is strictly contractive in the Thompson metric (40) on $\mathbb{R}^{N}_{>0}$, and admits unique fixed point $\bm{z}^{\rm{opt}}\in\mathbb{R}^{N}_{>0}$.

Proof.

Rewriting (41) as

[TABLE]

and letting $\bm{\eta}\equiv\bm{\eta}_{k,k+1}:=\bm{\xi}_{k-1}\oslash\left(\bm{\Gamma}_{k}^{\top}\bm{{\bm{\varrho}}}_{k-1}\right)$, we notice that iteration (41) can be expressed as a cone preserving composite map $\bm{\theta}:=\bm{\theta}_{1}\circ\bm{\theta}_{2}\circ\bm{\theta}_{3}$, where $\bm{\theta}:\mathbb{R}^{N}_{>0}\mapsto\mathbb{R}^{N}_{>0}$, given by

[TABLE]

and $\bm{\theta}_{1}(\bm{z}):=\bm{z}^{\frac{1}{1+\beta\epsilon/h}}$, $\bm{\theta}_{2}(\bm{z}):=\bm{\eta}\odot\bm{z}$, $\bm{\theta}_{3}:=\bm{\Gamma}_{k}\bm{z}$. Our strategy is to prove that the composite map $\bm{\theta}$ is contractive on $\mathbb{R}^{N}_{>0}$ w.r.t. the metric $d_{\rm{T}}$.

From (25) and (32), $\bm{C}_{k}(i,j)\in[0,\infty)$ which implies $\bm{\Gamma}_{k}(i,j)\in(0,1]$; therefore, $\bm{\Gamma}_{k}$ is a positive linear map for each $k=1,2,\ldots$. Thus, by (linear) Perron-Frobenius theorem, the map $\bm{\theta}_{3}$ is contractive on $\mathbb{R}^{N}_{>0}$ w.r.t. $d_{\rm{T}}$. The map $\bm{\theta}_{2}$ is an isometry by Definition 3. As for the map $\bm{\theta}_{1}$, notice that the quantity $r:=1/(1+\beta\epsilon/h)\in(0,1)$ since $\beta\epsilon/h>0$. Therefore, the map $\bm{\theta}_{1}(\bm{z}):=\bm{z}^{r}$ (element-wise exponentiation) is monotone (order preserving) and homogeneous of degree $r\in(0,1)$ on $\mathbb{R}^{N}_{>0}$. By Proposition 2, the map $\bm{\theta}_{1}(\bm{z})$ is strictly contractive. Thus, the composition

[TABLE]

is strictly contractive w.r.t. $d_{\rm{T}}$, and (by Banach contraction mapping theorem) admits unique fixed point $\bm{z}^{\rm{opt}}$ in $\mathbb{R}^{N}_{>0}$. ∎

Corollary 4**.**

The Algorithm 1 converges to unique fixed point $(\bm{y}^{\rm{opt}},\bm{z}^{\rm{opt}})\in\mathbb{R}^{N}_{>0}\times\mathbb{R}^{N}_{>0}$.

Proof.

Since $\bm{y}(:,\ell+1)=\bm{{\bm{\varrho}}}_{k-1}\oslash\left(\bm{\Gamma}_{k}\bm{z}(:,\ell+1)\right)$, the $\bm{z}$ iterates converge to unique fixed point $\bm{z}^{\rm{opt}}\in\mathbb{R}^{N}_{>0}$ (by Theorem 3), and the linear maps $\bm{\Gamma}_{k}$ are contractive (by Perron-Frebenius theory, as before), consequently the $\bm{y}$ iterates also converge to unique fixed point $\bm{y}^{\rm{opt}}\in\mathbb{R}^{N}_{>0}$. Hence the statement. ∎

IV Numerical Results

We now illustrate the computational framework proposed in Section III via numerical examples. Our examples involve systems already in JKO canonical form (Section II), as well as those which can be transformed to such form by non-obvious change of coordinates.

IV-A Linear Gaussian Systems

For an Itô SDE of the form

[TABLE]

it is well known that if ${\bm{x}}_{0}:={\bm{x}}(t=0)\sim\mathcal{N}(\bm{\mu}_{0},\bm{\Sigma}_{0})$, then the transient joint PDFs $\rho({\bm{x}},t)=\mathcal{N}(\bm{\mu}(t),\bm{\Sigma}(t))$ where the vector-matrix pair $\left(\bm{\mu}(t),\bm{\Sigma}(t)\right)$ evolve according to the ODEs

[TABLE]

We benchmark the numerical results produced by the proposed proximal algorithm vis-à-vis the solutions of (44). We consider the following two sub-cases of (43).

IV-A1 Ornstein-Uhlenbeck Process

We consider the univariate system

[TABLE]

which is in JKO form (13) with $\psi(x)=\frac{1}{2}ax^{2}$. We generate $N=400$ samples from the initial PDF $\rho_{0}=\mathcal{N}(\mu_{0},\sigma_{0}^{2})$ with $\mu_{0}=5$ and $\sigma_{0}^{2}=4\times 10^{-2}$, and apply the proposed proximal recursion for (45) with time step $h=10^{-3}$, and with parameters $a=1$, $\beta=1$, $\epsilon=5\times 10^{-2}$. For implementing Algorithm 1, we set the tolerance $\delta=10^{-3}$, and the maximum number of iterations $L=100$. Fig. 3 shows that the PDF point clouds generated by the proximal recursion match with the analytical PDFs $\mathcal{N}\left(\mu_{0}\exp(-at),(\sigma_{0}^{2}-\frac{1}{a\beta})\exp(-2at)+\frac{1}{a\beta}\right)$, and the mean-variance trajectories (computed from the numerical integration of the weighted scattered point cloud data) match with the corresponding analytical solutions.

IV-A2 Multivariate LTI

We next consider the multivariate case (43) where the pair $(\bm{A},\bm{B})$ is assumed to be controllable, and the matrix $\bm{A}$ is Hurwitz (not necessarily symmetric). Under these assumptions, the stationary PDF is $\mathcal{N}(\bm{0},\bm{\Sigma}_{\infty})$ where $\bm{\Sigma}_{\infty}$ is the unique stationary solution of (44b) that is guaranteed to be symmetric positive definite. However, it is not apparent whether (43) can be expressed in the form (13), since for non-symmetric $\bm{A}$, there does not exist constant symmetric positive definite matrix $\bm{\Psi}$ such that $\bm{Ax}=-\nabla{\bm{x}}^{\top}\bm{\Psi}{\bm{x}}$, i.e., the drift vector field does not admit a natural potential. Thus, implementing the JKO scheme for (43) is non-trivial in general.

In our recent work [18], two successive time-varying co-ordinate transformations were proposed which can bring (43) in the form (13), thus making it amenable to the JKO scheme. We apply these change-of-coordinates to (43) with

[TABLE]

which satisfy the stated assumptions for the pair $(\bm{A},\bm{B})$, and implement the proposed proximal recursion on this transformed co-ordinates with $N=400$ samples generated from the initial PDF $\rho_{0}=\mathcal{N}(\bm{\mu}_{0},\bm{\Sigma}_{0})$, where $\bm{\mu}_{0}=(4,4)^{\top}$ and $\bm{\Sigma}_{0}=4\bm{I}_{2}$. As before, we set $\delta=10^{-3},L=100,h=10^{-3},\beta=1,\epsilon=5\times 10^{-2}$. Once the proximal updates are done, we transform back the probability weighted scattered point cloud to the original state space co-ordinates via change-of-measure formula associated with the known co-ordinate transforms [18, Section III.B]. Fig. 4 shows the resulting point clouds superimposed with the contour plots for the analytical solutions $\mathcal{N}(\bm{\mu}(t),\bm{\Sigma}(t))$ given by (44). Figs. 5 and 6 compare the respective mean and covariance evolution. We point out that the change of co-ordinates in [18] requires implementing the JKO scheme in a time-varying rotating frame (defined via exponential of certain time varying skew-symmetric matrix) that depends on the stationary covariance $\bm{\Sigma}_{\infty}$. As a consequence, the stationary covariance resulting from the proximal recursion oscillates about the true stationary value.

IV-B Nonlinear Non-Gaussian System

Next we consider a planar nonlinear system of the form (13) with $\psi(x_{1},x_{2})=\dfrac{1}{4}(1+x_{1}^{4})+\dfrac{1}{2}(x_{2}^{2}-x_{1}^{2})$ (see Fig. 7). As mentioned in Section II, the stationary PDF is $\rho_{\infty}({\bm{x}})=\kappa\exp\left(-\beta\psi({\bm{x}})\right)$, which for our choice of $\psi$, is bimodal. In this case, the transient PDFs have no known analytical solution but can be computed using the proposed proximal recursion. For doing so, we generate $N=400$ samples from the initial PDF $\rho_{0}=\mathcal{N}(\bm{\mu}_{0},\bm{\Sigma}_{0})$ with $\bm{\mu}_{0}=(2,2)^{\top}$ and $\bm{\Sigma}_{0}=4\bm{I}_{2}$, and set $\delta=10^{-3},L=100,h=10^{-3},\beta=1,\epsilon=5\times 10^{-2}$, as before. The resulting weighted point clouds are shown in Fig. 8; it can be seen that as time progresses, the joint PDFs computed via the proximal recursion, tend to the known stationary solution $\rho_{\infty}$ (contour plots in the right bottom sub-figure in Fig. 8).

Fig. 9 shows the computational times for the proposed proximal recursions applied to the above nonlinear non-Gaussian system. Since the proposed algorithm involves sub-iterations (“while loop” in Algorithm 1 over index $\ell\leq L$) while keeping the physical time “frozen”, the convergence reported in Section III-C must be achieved at “sub-physical time step” level, i.e., must incur smaller than $h$ (here, $h=10^{-3}$ s) computational time. Indeed, Fig. 9 shows that each proximal update takes approx. $10^{-6}$ s, or $10^{-3}h$ computational time, which demonstrates the efficacy of the proposed framework.

IV-C Non-local Interactions

We now consider a numerical example demonstrating our gradient flow framework for computing the transient PDFs generated by the McKean-Vlasov integro-PDE (18) associated with the density-dependent sample path dynamics (19). If both the potentials $\psi$ and $\phi$ are convex functions, then (18) admits unique stationary density $\rho_{\infty}({\bm{x}})$ (see e.g., [31],[45, Section 5.2]). To keep the exposition simple, we consider the univariate case $\psi(x)=\frac{1}{2}ax^{2}$, $\phi(x)=\frac{1}{2}bx^{2}$, $a,b>0$. Performing the integration appearing in the convolution allows us to rewrite (19) as the mean-reverting process:

[TABLE]

where $\mu(t)$ is the mean of the transient PDF $\rho(x,t)$. Assuming $x_{0}\sim\mathcal{N}(\mu_{0},\sigma_{0}^{2})$, and applying expectation operator to both sides of (46) yields $\mu(t)=\mu_{0}\exp(-at)$. Consequently, the transient PDF for (46) at time $t$ is $\rho(x,t)=\mathcal{N}\left(\mu(t),\sigma^{2}(t)\right)$, with

[TABLE]

Clearly, the stationary PDF is $\rho_{\infty}(x)=\mathcal{N}\left(0,1/(a+b)\beta\right)$.

To benchmark our algorithm with the analytical solution (47), we implement the proximal recursion for (46) with free energy (20). Following [46, Section 4], we replace the non-convex bilinear term $\int_{\mathbb{R}^{n}\times\mathbb{R}^{n}}\phi({\bm{x}}-\bm{y})\varrho({\bm{x}})\varrho(\bm{y}){\rm{d}}{\bm{x}}{\rm{d}}\bm{y}$ in (22b), with the linear term $\int_{\mathbb{R}^{n}}\phi({\bm{x}}-\bm{y})\varrho({\bm{x}})\varrho_{k-1}(\bm{y}){\rm{d}}{\bm{x}}{\rm{d}}\bm{y}$, $k\in\mathbb{N}$, resulting in a semi-implicit variant of (22), given by

[TABLE]

with $\varrho_{0}\equiv\rho_{0}({\bm{x}})$ (the initial PDF). The proof for the fact that such a semi-implicit scheme guarantees $\varrho_{k}({\bm{x}})\rightarrow\rho(\bm{x},t=kh)$ for $h\downarrow 0$, where $\rho({\bm{x}},t)$ is the flow generated by (18), can be found in [47, Section 12.3]. Notice that for the FPK gradient flow, the “discrete free energy” in (23) was $F({\bm{\varrho}})=\langle\bm{\psi}_{k-1}+\beta^{-1}\log{\bm{\varrho}},{\bm{\varrho}}\rangle$. The scheme (48) allows us to write a similar expression in the McKean-Vlasov gradient flow case, as $F({\bm{\varrho}})=\langle\bm{\psi}_{k-1}+\bm{D}_{k-1}{\bm{\varrho}}_{k-1}+\beta^{-1}\log{\bm{\varrho}},{\bm{\varrho}}\rangle$, where the symmetric matrix $\bm{D}_{k-1}$ is given by

[TABLE]

For the particular choice $\phi(x)=\frac{1}{2}bx^{2}$, notice that $\bm{D}_{k-1}$ is a (scaled) Euclidean distance matrix on the Euler-Maruyama update. In this case, the associated Euler-Maruyama scheme (39) is

[TABLE]

i.e., the dashed arrow in Fig. 2 becomes active. Fig. 10 shows that the weighted scattered point cloud solutions for (46) with $a=b=1$, computed through our proximal algorithm match with the analytical solutions $\mathcal{N}(\mu(t),\sigma^{2}(t))$ given by (47). For Fig. 10, the parameter values used in our simulation are $h=10^{-3}$, $\beta=1$, $\epsilon=5\times 10^{-2}$, $\mu_{0}=5$, $\sigma_{0}^{2}=9$, $N=400$, $\delta=10^{-3}$, $L=100$.

V Extensions

In Section IV-A2, we have already seen that systems not in JKO canonical form may be transformed to the same via suitable change-of-coordinates, thus making density propagation for such systems amenable via our framework. In this Section, we provide two extensions along these lines. First, we consider a case of state-dependent diffusion; thereafter, we consider a system with mixed conservative-dissipative drift. In both cases, we use specific examples (instead of general remarks) to help illustrate the extensions of the basic framework. These examples point out the broad scope of the algorithms proposed herein.

V-A Multiplicative Noise

We consider the Itô SDE for the Cox-Ingersoll-Ross (CIR) model [48], given by

[TABLE]

Due to multiplicative noise, (50) is not in JKO canonical form (13). If the initial PDF $\rho_{0}$ is Dirac delta at $x_{0}$, then the FPK PDE for (50) admits closed form solution:

[TABLE]

where $I_{q}(\cdot)$ is the modified Bessel function of order $q$, and

[TABLE]

The transient solutions (51) are non-central chi-squared PDFs. The stationary solution is a Gamma PDF $\rho_{\infty}(x)\propto x^{q}\exp\left(-2ax/b^{2}\right)$, $x>0$. We will benchmark our proximal algorithm against (51).

In order to transcribe (50) in the JKO canonical form, we employ the Lamperti transform [49, 50], where the idea is to find a change of variable $y=\varsigma(x)$ such that the SDE for new variable $y$ has unity diffusion coefficient. From Itô’s lemma, it follows that the requisite transformation for (50) is $y=\varsigma(x):=\frac{2}{b}\sqrt{x}$, and the resulting SDE in $y$ becomes

[TABLE]

Clearly, (53) is in the JKO canonical form (13) with $\beta=2$, and

[TABLE]

So the proposed proximal algorithm (Algorithm 1) can be applied to (53), and at each time $t$, the resulting PDF $\rho_{Y}(y,t)$ can be transformed back to the PDF $\rho_{X}(x,t)$ via the change-of-measure formula for the push-forward map $x=\varsigma^{-1}(y)=by^{2}/4$. Fig. 11 shows the comparison of the analytical and proximal transient PDFs for (50) resulting from $\rho_{0}=\delta(x-5)$ with $a=3$, $b=2$, $\theta=2$, $h=10^{-3}$, $\epsilon=5\times 10^{-2}$, $N=400$, $\delta=10^{-3}$, and $L=100$.

V-B Mixed Conservative-Dissipative Drift

In many engineering applications, one encounters Itô SDEs where the drift vector fields have both dissipative (gradient) and conservative (Hamiltonian) components. For example, stochastic systems arising from Newton’s law in mechanics often have mixed conservative-dissipative structure. Our intent here is to illustrate that such SDEs are amenable to the proposed proximal recursion framework. As an example, we will work out the details for a perturbed two-body problem in celestial mechanics similar to the one treated in [51].

We consider the relative motion of a satellite in geocentric orbit, given by the second order Langevin equation

[TABLE]

where $\bm{q}:=(x,y,z)^{\top}\in\mathbb{R}^{3}$ is the relative position vector for the satellite, $\mu$ is a constant (product of the Gavitational constant and the mass of Earth), $-\gamma\dot{\bm{q}}$ models linear drag222More generally, for nonlinear drag of the form $-\gamma\nabla_{\dot{\bm{q}}}\widetilde{V}(\dot{\bm{q}})$, the term $-\gamma\bm{p}$ in (59) will become $-\gamma\nabla\widetilde{V}(\bm{p})$. This will entail modifying the functional $\widehat{F}$ in (61) as $\mathbb{E}_{\rho}\left[\widetilde{V}(\bm{p})\>+\>\beta^{-1}\log\rho\right]$. The linear drag illustrated here corresponds to the special case $\widetilde{V}(\bm{p})\equiv\parallel\bm{p}\parallel_{2}^{2}/2$., $\bm{f}_{\text{pert}}(\bm{q})$ models (deterministic) perturbative force due to the oblateness of Earth, and the stochastic forcing is due to solar radiation pressure, free-molecular aerodynamic forcing etc. Using the shorthands for sines and cosines as $c\tau:=\cos\tau$, $s\tau:=\sin\tau$, and recalling the relations among spherical coordinates $(r,\theta,\phi)$ and cartesian coordinates $(x,y,z)$, given by $r=\sqrt{x^{2}+y^{2}+z^{2}}$, $c\phi=x/r$, $s\phi=\sqrt{1-(c\phi)^{2}}$, $c\theta=z/r$, $s\theta=\sqrt{1-(c\theta)^{2}}$, we can write $\bm{f}_{\text{pert}}(\bm{q})$ in spherical coordinates as [51, eqn. (12)-(13)]

[TABLE]

and the same in cartesian coordinates as

[TABLE]

In (56), the Earth oblateness coefficient $J_{2}=1.082\times 10^{-3}$, the radius of Earth $R_{\rm{E}}=6.3781\times 10^{6}$ m, and the Earth standard Gravitational parameter $\mu=3.9859\times 10^{14}$ m3/s2. Modeling the stochastic forcing in (55) as standard Gaussian white noise (as in [51]), (55) can then be expressed as an Itô SDE in $\mathbb{R}^{6}$:

[TABLE]

or more succinctly,

[TABLE]

where $\bm{p}:=(v_{x},v_{y},v_{z})^{\top}\in\mathbb{R}^{3}$ is the velocity vector, and

[TABLE]

Introducing a “Hamiltonian-like function” $H(\bm{q},\bm{p}):=\parallel\bm{p}\parallel_{2}^{2}/2+V(\bm{q})$, one can verify that the stationary PDF for (59) is $\rho_{\infty}\propto\exp(-\beta H(\bm{q},\bm{p}))$, and that the “total free energy” $F(\rho):=\mathbb{E}_{\rho}\left[H+\beta^{-1}\log\rho\right]$ serves as Lyapunov functional for the associated FPK PDE, i.e., $\frac{{\rm{d}}F}{{\rm{d}}t}<0$. However, the proximal recursion (22b) does not apply as is, instead needs to be modified to account the joint conservative-dissipative effect as

[TABLE]

We implement a recursion from [52, Scheme 2b] where $\widehat{F}(\varrho):=\mathbb{E}_{\rho}\left[\frac{1}{2}\parallel\bm{p}\parallel_{2}^{2}\>+\>\beta^{-1}\log\rho\right]$, and $\widehat{W}_{h}^{2}(\varrho_{k-1},\varrho)$ is the optimal mass transport cost (as in squared 2-Wasserstein metric (10)) with modified cost function (modified integrand in (10)), i.e.,

[TABLE]

where $(\bm{q},\bm{p})^{\top}$ and $(\widetilde{\bm{q}},\widetilde{\bm{p}})^{\top}$ are two realizations of the state vector governed by (59)-(60); the integrand in (62) is

[TABLE]

That the proximal recursion (61) with the above choices of functionals $\widehat{F}$ and $\widehat{W}_{h}$ guarantees $\varrho_{k}(\bm{q},\bm{p},h)\rightarrow\rho(\bm{q},\bm{p},t=kh)$ for $k\in\mathbb{N}$ as $h\downarrow 0$, where $\rho(\bm{q},\bm{p},t)$ is the joint PDF generated by the FPK flow for (59) at time $t$, was proved in [52]. Our proximal algorithm in Section III-B applies by simply modifying (24)-(25) as

[TABLE]

respectively. Notice that $\widehat{s}_{h}$ in (63) (and consequently, $\widehat{W}_{h}$ in (62)) is not a metric (in particular, non-symmetric). Thus, the matrices $\bm{C}_{k}$ in (65) for $k\in\mathbb{N}$ are not symmetric.

To apply the proposed framework, we first non-dimensionalize the variables $\bm{q}$ (m), $\bm{p}$ (m/s), $t$ (s), $\bm{w}$ ($\sqrt{\text{s}}$) in (59) as

[TABLE]

where $R:=4.2164\times 10^{7}$ m is the radius of the nominal geostationary orbit, and $T=86164$ s is its period. In (66), the primed variables are non-dimensionalized. Using (66) and Itô’s lemma on (59), the non-dimensional SDE in $(\bm{q}^{\prime},\bm{p}^{\prime})^{\top}\in\mathbb{R}^{6}$ becomes

[TABLE]

To avoid numerical conditioning issues, we will apply our algorithm for recursion (61) associated with the non-dimensional SDE (67), and then transform the PDF in $(\bm{q}^{\prime},\bm{p}^{\prime})^{\top}$ to the same in original variables $(\bm{q},\bm{p})^{\top}$ via change-of-measure formula. Such a computational pipeline, i.e., density propagation in non-dimensional SDE and then transforming the density back in dimensional variables, is standard in celestial mechanics due to different orders of magnitude in different state variables (see e.g., [51, Sec. III.B]).

Figs. 12 and 13 show the univariate marginal PDFs at $t=0.005$s and $t=0.010$s respectively, associated with the joint PDFs supported on $\mathbb{R}^{6}$ at those times, computed through the proposed proximal algorithm. For this simulation, the initial joint PDF $\rho_{0}(\bm{q}^{\prime},\bm{p}^{\prime})=\mathcal{N}\left(\bm{\mu}_{0},\bm{\Sigma}_{0}\right)$, where $\bm{\mu}_{0}=\left(1,0,0,0,0,0\right)^{\top}$, and $\bm{\Sigma}_{0}=10^{-4}\times\text{diag}\left(3.335,6.133,3.933,6.562,9.246,5.761\right)$. The parameter values used in the simulation are: $h=10^{-5}$, $\beta=1$ m2/s2, $\gamma=1$ s*-1*, $\epsilon=5\times 10^{-2}$, $\delta=10^{-3}$, $L=100$, and $N=400$. The computational times for this 6-state example reported in Fig. 14 reveal that the proximal recursions run remarkably fast (indeed, faster than the runtime for the 2-state example in Section IV-B, cf. Fig. 9).

VI Conclusions

In this paper, novel uncertainty propagation algorithms are presented for computing the flow of the joint PDFs associated with continuous-time stochastic nonlinear systems. By interpreting the PDF flow as gradient descent on the manifold of joint PDFs w.r.t. a suitable metric, the proposed computational framework implements proximal algorithms which are proved to be convergent due to certain contraction properties established herein. Numerical examples are provided to demonstrate the practical use of the proposed algorithm and its efficiency in terms of computational runtime. In contrast to the conventional function approximation algorithms for this problem, the proposed non-parametric framework does not make any spatial discretization, instead performs finite sample probability-weighted scattered data computation in the form of temporal recursion. The location of the atomic measures is delegated to a Euler-Maruyama scheme, thus avoiding spatial discretization. The results of this paper provide computational teeth to the emerging systems-theoretic viewpoint [18, 19] that the PDF flows in uncertainty propagation can be seen as gradient flow.

-A Derivation of (28):

Noting that scaling and translation by constants do not alter the outer argmin in (27), we rewrite the same as

[TABLE]

since $\bm{1}^{\top}\bm{M}\bm{1}=1$. The Lagrangian $\mathcal{L}$ associated with the inner minimization in (68) is given by

[TABLE]

Setting the derivative of (69) w.r.t. the $(i,j)$-th element of $\bm{M}$ equal to zero, followed by algebraic simplification yields (28).