## Gradient Flow Algorithms for Density Propagation in Stochastic Systems

Kenneth F. Caluya, and Abhishek Halder

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.

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]-[7], and in motion planning [8]-[10]. In these applications, the quantity of interest is the joint PDF ρ ( x , t ) and its statistics (e.g., transient moments and marginal PDFs). The second type of applications require ρ ( x , t ) as an intermediate step toward computing other quantities of interest. For example, in nonlinear filtering [11], [12], the joint PDF ρ ( x , t ) serves as the prior in computing the posterior (i.e., conditional state) PDF. In probabilistic model validation [13]-[15] and controller verification [16], computing ρ ( x , t ) helps in quantifying the density-level prediction-observation mismatch. All these applications require fast computation of ρ ( x , t ) in a scalable and unified manner, as opposed to developing algorithms in a case-by-case basis.

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

## I. INTRODUCTION

Consider the continuous-time dynamics of the state vector x ( ) t ∈ R n governed by an Itˆ o stochastic differential equation (SDE)

$$\stackrel { \cdot } { \cdot } \quad \text{d} x = f \left ( x, t \right ) \text{d} t \, + \, g ( x, t ) \text{d} w, \, \ x ( t = 0 ) = x _ { 0 }, \quad \ \left ( 1 \right ) \text{ } \text{ but } \text{c}$$

where the joint probability density function (PDF) for the initial condition x 0 is a known function ρ 0 ; we use the notation x 0 ∼ ρ 0 . The process noise w ( ) t ∈ R m is Wiener and satisfy E [d w i d w j ] = δ ij d t for all i, j = 1 , . . . , n , where δ ij = 1 for i = j , and zero otherwise. Then, the flow of the joint PDF ρ ( x , t ) for the state vector x ( ) t (i.e., x ∼ ρ ) is governed by the partial differential equation (PDE) initial value problem:

$$\frac { \partial \rho } { \partial t } = - \nabla \cdot ( \rho f ) + \frac { 1 } { 2 } \sum _ { i, j = 1 } ^ { n } \frac { \partial ^ { 2 } } { \partial x _ { i } \partial x _ { j } } ( \rho g g ^ { \top } ) _ { i j }, \quad ( 2 a )$$

$$\rho ( x, t = 0 ) = \rho _ { 0 } ( x ) \ \ ( \text{given} ). \quad ( 2 b )$$

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 ρ ( x , t ) that solves a PDE initial value problem of the form (2). From an application standpoint, the need for computing ρ ( 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

K.F. Caluya, and A. Halder are with the Department of Applied Mathematics, University of California, Santa Cruz, CA 95064, USA; kcaluya,ahalder@ucsc.edu

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 ρ ( 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

$$\frac { d x } { d t } = - \nabla \varphi \left ( x \right ), \ \ x ( 0 ) = x _ { 0 },$$

where x x , 0 ∈ R n , and the continuously differentiable function ϕ : R n → R ≥ 0 . The flow x ( ) t generated by (3) can be realized via variational recursion

$$\mathfrak { r } \text{ belief} _ { \text{$m$-PDF} } \quad x _ { k } = \arg \min _ { x \in \mathbb { R } ^ { n } } \frac { 1 } { 2 } \, \| \, x - x _ { k - 1 } \, \| _ { 2 } ^ { 2 } + h \, \varphi ( x ) + o ( h ), \ \ k \in \mathbb { N }, \ \ ( 4 )$$

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

$$x _ { k } - x _ { k - 1 } = - h \nabla \varphi ( x _ { k - 1 } ),$$

Fig. 1: The gradient descent on the manifold of PDFs can be described by successive evaluation of proximal operators to recursively update PDFs from time t = ( k -1) h to t = kh for k ∈ N , and time-step h &gt; 0 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-a5666c51-328a-4564-a8d8-07bbee9ecfa4.jpg)

as

$$x _ { k } & = \arg \min _ { x \in \mathbb { R } ^ { n } } \frac { 1 } { 2 } \left \| x - ( x _ { k - 1 } - h \nabla \varphi ( x _ { k - 1 } ) ) \right \| _ { 2 } ^ { 2 } & \text{appeared} \\ & = \arg \min _ { x \in \mathbb { R } ^ { n } } \frac { 1 } { 2 } \left \| x - x _ { k - 1 } \right \| _ { 2 } ^ { 2 } + \left \langle x - x _ { k - 1 }, h \nabla \varphi ( x _ { k - 1 } ) \right \rangle & \text{prod} \\ & \quad + \ h \varphi ( x _ { k - 1 } ), \ \ ( 5 ) \quad 1 ) \\ \cdot \quad \cdot \cdot \cdot \ e \ \dots \ \cdot \cdot \cdot \cdot \quad \cdot \cdot \cdot \quad \cdot \cdot \quad \text{letters $f$}$$

↦

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 ϕ ( x ) at x = x k . In (4), the mapping x k -1 → x k given by

$$\text{prox} _ { h \varphi } ^ { \| \cdot \| _ { 2 } } ( x _ { k - 1 } ) \coloneqq \underset { x \in \mathbb { R } ^ { n } } { \arg \min } \frac { 1 } { 2 } \left \| x - x _ { k - 1 } \left \| _ { 2 } ^ { 2 } + h \, \varphi ( x ), \quad ( 6 ) \quad \text{and} \, \Delta \, d$$

is called the 'proximal operator' [20, p. 142] of hϕ w.r.t. the standard Euclidean metric ‖ · ‖ 2 . Notice that ϕ ( ) · serves as a Lyapunov function since the quantity

$$\frac { d } { d t } \varphi = \langle \nabla \varphi, - \nabla \varphi \rangle = - \left \| \, \nabla \varphi \, \right \| _ { 2 } ^ { 2 } \, \quad \quad ( 7 ) \, \begin{smallmatrix} i. e., \, \mathbb { E } _ {, } \\ \text{to be} \\ \text{altern} \end{smallmatrix}$$

equals 0 at the stationary point of (3), and &lt; 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 ∗

↦

$$\mathcal { D } _ { 2 } \coloneqq \{ \rho \, \colon \mathbb { R } ^ { n } \mapsto \mathbb { R } \, | \, \rho \geq 0, \int _ { \mathbb { R } ^ { n } } \rho = 1, \, \mathbb { E } _ { \rho } [ x ^ { \top } x ] < \infty \}. \quad \ \ \ t i e _ { \mathbb { W } }$$

↦

Specifically, let d ( · , · ) be a distance metric on the manifold D 2 , and let the functional Φ : D 2 → R ≥ 0 . Then, for some chosen step-size h &gt; 0 , the infinite dimensional proximal operator of h Φ w.r.t. the distance metric d ( · , · ) , given by

$$\text{prox} _ { h \Phi } ^ { d ( \cdot ) } ( \varrho _ { k - 1 } ) \coloneqq \underset { \varrho \in \mathcal { D } _ { 2 } } { \arg \inf } \frac { 1 } { 2 } d ( \varrho, \varrho _ { k - 1 } ) ^ { 2 } + h \, \Phi ( \varrho ), \quad ( 8 ) \quad i t \, i s \, d e f i$$

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

$$\varrho _ { k } = \text{prox} _ { h \Phi } ^ { d ( \cdot ) } ( \varrho _ { k - 1 } ), \ \ k \in \mathbb { N }, \ \varrho _ { 0 } ( x ) \coloneqq \rho _ { 0 } ( x ). \quad \ ( 9 )$$

Just as the proximal recursion (4) approximates the finite dimensional flow (3), similarly it is possible to design d ( · , · ) and Φ( ) · in (8) as function of the drift f and diffusion 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 f is a gradient vector field and g is a scalar multiple of identity matrix, then the distance d ( · , · ) can be taken as the Wasserstein-2 metric [22] with Φ( ) · as the free energy functional. In particular, the solution of (9) was shown to converge to the flow of (2), i.e., glyph[rho1] k ( x ) → ρ ( x , t = kh ) in strong L 1 ( R n ) sense, as h ↓ 0 . The resulting variational recursion

∗ We denote the expectation operator w.r.t. the measure ρ ( x )d x as E ρ [ · ] .

(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 f and 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 Φ( ) · . 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).

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 〈· , ·〉 to denote the Euclidean inner product. In particular, 〈 A B , 〉 := trace( A B glyph[latticetop] ) denotes Frobenius inner product between matrices A and B , and 〈 a b , 〉 := a glyph[latticetop] b denotes the inner product between column vectors a and b . The symbols ∇ and ∆ denote the (Euclidean) gradient and the Laplacian operators, respectively. We use N ( µ, σ 2 ) to denote a univariate Gaussian PDF with mean µ and variance σ 2 . Likewise, N ( µ , Σ ) denotes a multivariate Gaussian PDF with mean vector µ and covariance matrix Σ . The expectation operator w.r.t. the PDF ρ is denoted by E ρ [ ] · , i.e., E ρ [ · ] := ∫ R n ( ) · ρ d x . The operands log( ) · , exp( ) · and ≥ 0 are to be understood as element-wise. The notations glyph[circledot] and glyph[circledivide] denote element-wise (Hadamard) product and division, respectively. We use I n to denote the n × n identity matrix. The symbols 1 and 0 stand for column vectors of appropriate dimension containing all ones, and all zeroes, respectively.

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 d π 1 ( x ) := ρ 1 ( x )d x and d π 2 ( y ) := ρ 2 ( y )d y , supported respectively on X Y , ⊆ R n , is denoted as W π ,π ( 1 2 ) (equivalently, W ρ ,ρ ( 1 2 ) whenever the measures π , π 1 2 are absolutely continuous so that the respective PDFs ρ , ρ 1 2 exist); it is defined as

$$W ( \pi _ { 1 }, \pi _ { 2 } ) \coloneqq \left ( \inf _ { d \pi \in \Pi ( \pi _ { 1 }, \pi _ { 2 } ) } \int _ { \mathcal { X } \times \mathcal { Y } } s ( x ; y ) \, d \pi \left ( x, y \right ) \right ) ^ { \frac { 1 } { 2 } }, \quad ( 1 0 )$$

where s ( x y ; ) := ‖ x -y ‖ 2 2 is the squared Euclidean distance in R n , and Π( π , π 1 2 ) denotes the collection of all joint probability measures on X × Y having finite second moments, with marginals π 1 and π 2 , respectively.

The existence and uniqueness of the minimizer in (10) is guaranteed. It is well-known [22, Ch. 7] that W π ,π ( 1 2 ) defines a metric on the manifold D 2 . This means that W π ,π ( 1 2 ) ≥ 0 with equality if and only if π 1 = π 2 , the symmetry: W π ,π ( 1 2 ) = W π ,π ( 2 1 ) , and that W π ,π ( 1 2 ) satisfies the triangle inequality. Its square, W 2 ( π , π 1 2 ) equals [26] the minimum amount of work required to transport π 1 to π 2 (or equivalently, ρ 1 to ρ 2 ), and vice versa. For

↦

any PDF ν , the function ρ → W 2 ( ρ, ν ) is convex on D 2 , i.e., for any ρ , ρ 1 2 ∈ D 2 , and 0 ≤ τ ≤ 1 , we have

$$W ^ { 2 } ( \nu, ( 1 - \tau ) \rho _ { 1 } + \tau \rho _ { 2 } ) \, \leq \, ( 1 - \tau ) W ^ { 2 } ( \nu, \rho _ { 1 } ) + \tau W ^ { 2 } ( \nu, \rho _ { 2 } ). \ ( 1 1 )$$

Definition 2. ( Kullback-Leibler divergence ) The Kullback-Leibler divergence, also known as relative entropy, between two probability measures d π 1 ( x ) := ρ 1 ( x )d x and d π 2 ( y ) := ρ 2 ( y )d y , denoted as D KL ( π 1 || π 2 ) , is defined as

$$\cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \cdots, \\ D _ { K L } ( \pi _ { 1 } | | \pi _ { 2 } ) & = \int \left ( \frac { d \pi _ { 1 } } { d \pi _ { 2 } } \right ) \log \left ( \frac { d \pi _ { 1 } } { d \pi _ { 2 } } \right ) d \pi _ { 2 } & \quad \begin{array} { c } \text{as the} \\ \text{energy} \\ \text{heat eq} \\ \text{$e$} \end{array} \\ \cdots, \cdots$$

where d π 1 d π 2 denotes the Radon-Nikodym derivative. Henceforth, we use the notational equivalence D KL ( π 1 || π 2 ) ≡ D KL ( ρ 1 || ρ 2 ) .

From Jensen's inequality, D KL ( ρ 1 || ρ 2 ) ≥ 0 ; however, D 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 ρ 1 and ρ 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 x ( ) t and its joint PDF ρ ( 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 ρ ( 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.

## A. FPK Gradient Flow

## Consider the Itˆ o SDE

$$\mathrm d x = - \nabla \psi ( x ) \, \mathrm d t \, + \, \sqrt { 2 \beta ^ { - 1 } } \, \mathrm d w, \ x ( 0 ) = x _ { 0 }, \quad ( 1 3 ) \, \ x - \mathbb { E } _ { \rho } [$$

↦

where the the drift potential ψ : R n → (0 , ∞ ) , the diffusion coefficient β &gt; 0 , and the initial condition x 0 ∼ ρ 0 ( x ) . For the sample path x ( ) t dynamics given by the SDE (13), the flow of the joint PDF ρ ( x , t ) is governed by the FPK PDE initial value problem

$$\frac { \partial \rho } { \partial t } = \nabla \cdot ( \rho \nabla \psi ) \, + \, \beta ^ { - 1 } \Delta \rho, \ \rho ( x, 0 ) = \rho _ { 0 } ( x ). \quad ( 1 4 ) \quad \text{general} i$$

It is easy to verify that the unique stationary solution of (14) is the Gibbs PDF ρ ∞ ( x ) = κ exp( -βψ ( x )) , where the normalizing constant κ is known as the partition function .

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

$$F ( \rho ) \ \coloneqq \ \ & \mathbb { E } _ { \rho } \left [ \psi + \beta ^ { - 1 } \log \rho \right ] \\ & = \ \beta ^ { - 1 } D _ { \text{KL} } \left ( \rho \left \| \exp \left ( - \beta \psi ( \boldsymbol x ) \right ) \right ) \geq 0, \quad ( 1 5 )$$

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

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

and consequently

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

which is &lt; 0 for the transient solution ρ ( x , t ) , and = 0 at the stationary solution ρ ∞ ( x ) = κ exp( -βψ ( 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 ∫ R n ψ ( x ) ρ d x and the internal energy β -1 ∫ R n ρ log ρ d x . If ψ ≡ 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 ≡ W (i.e., the 2-Wasserstein metric in (10)) and the functional Φ ≡ 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 Φ . 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.

## 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

$$\stackrel { \cdot } { \ u t e } \, \text{the} \quad \frac { \partial \rho } { \partial t } = \nabla \cdot \left ( \rho \nabla \left ( \psi + \rho * \phi \right ) \right ) \, + \, \beta ^ { - 1 } \Delta \rho, \quad \rho ( x, 0 ) = \rho _ { 0 } ( x ), \ \ ( 1 8 )$$

↦

where ∗ denotes the convolution in R n , the interaction potential φ : R n → (0 , ∞ ) and is symmetric, i.e., φ ( -v ) = φ ( v ) for v ∈ R n . The associated sample path x ( ) t dynamics has PDF-dependent drift:

$$d x = - \left ( \nabla \psi \left ( x \right ) + \nabla \left ( \rho * \phi \right ) \right ) d t + \sqrt { 2 \beta ^ { - 1 } } \, d w, \ \ x ( 0 ) = x _ { 0 }. \ ( 1 9 )$$

As an example, when φ ( v ) := 1 2 ‖ v ‖ 2 , then ∇ ( ρ ∗ φ )( x ) = x -E ρ [ x ] . Clearly, (18) reduces to (14) in the absence of interaction ( φ ≡ 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

$$F ( \rho ) \, \coloneqq \, \mathbb { E } _ { \rho } \left [ \psi \, + \, \beta ^ { - 1 } \log \rho \, + \, \rho * \phi \right ], \quad \quad ( 2 0 )$$

which is a sum of the potential energy (as before), the internal energy (as before), and the interaction energy 1 2 ∫ R 2 n φ ( x -y ) ρ ( x ) ρ ( y ) d x y d . In this case, (17) holds with

$$\zeta \coloneqq \beta ^ { - 1 } \left ( 1 + \log \rho \right ) + \psi + \rho * \phi, \text{ \quad \quad } \text{(2l)}$$

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

## III. FRAMEWORK

We now describe our computational framework to solve the proximal recursion

$$\varrho _ { k } & = \text{prox} _ { h F ( \cdot ) } ^ { W } ( \varrho _ { k - 1 } ) & & ( 2 2 a ) \\ & = \underset { \varrho \in \mathcal { D } _ { 2 } } { \arg \inf } \, \frac { 1 } { 2 } W ^ { 2 } ( \varrho _ { k - 1 }, \varrho ) + h \, F ( \varrho ), \quad k \in \mathbb { N }, \quad ( 2 2 b ) \\ & & \text{c:dots} \, \cdot$$

with glyph[rho1] 0 ≡ ρ 0 ( x ) (the initial joint PDF) for some small step-size h . For maximal clarity, we develop the framework with the free energy F ( ) · as in (15), i.e., for FPK gradient flow. In Section IV.C, we will show how the same can be generalized when F ( ) · 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 ∈ N .

We discretize time as t = 0 , h, 2 h, . . . , 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 { x i k , glyph[rho1] i k } N i =1 of cardinality N at t k = kh , k ∈ N , where the location of the i th point x i k ∈ R n denotes its state-space coordinate, and the corresponding weight glyph[rho1] i k ∈ R ≥ 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:

$$\varrho _ { k } = \arg \min _ { \varrho } \left \{ \min _ { M \in \Pi ( \varrho _ { k - 1 }, \varrho ) } \frac { 1 } { 2 } \langle C _ { k }, M \rangle + h \, \langle \psi _ { k - 1 } \quad \quad A. \ P r o x \\ + \beta ^ { - 1 } \log \varrho, \varrho \rangle \right \}, \ \ k \in \mathbb { N }, \quad ( 2 3 ) \quad \ \.$$

where the drift potential vector ψ k -1 ∈ R N is given by

$$\psi _ { k - 1 } ( i ) \coloneqq \psi \left ( x _ { k - 1 } ^ { i } \right ), \ \ i = 1, \dots, N. \quad \ \ ( 2 4 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

Here, the probability vectors glyph[rho1] , glyph[rho1] k -1 ∈ S N -1 , the probability simplex in R N . Furthermore, for each k ∈ N , the matrix C k ∈ R N × N is given by

$$C _ { k } ( i, j ) \coloneqq s ( x _ { k } ^ { i } ; x _ { k - 1 } ^ { j } ) = & \| \, x _ { k } ^ { i } - x _ { k - 1 } ^ { j } \, \| _ { 2 } ^ { 2 }, \ i, j = 1, \dots, N, \ ( 2 5 ) ^ { \ \ j o r \, y, z \, \omega }$$

and Π( glyph[rho1] k -1 , glyph[rho1] ) stands for the set of all matrices M ∈ R N × N such that

$$M \geq 0, \ \ M 1 = \varrho _ { k - 1 }, \ \ M ^ { \top } 1 = \varrho. \quad \ ( 2 6 ) \quad a n d \ t$$

Notice that the inner minimization in (23) is a standard linear programming problem if it were to be solved for a given glyph[rho1] ∈ 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 ( M ) := 〈 M , log M 〉 in (23) yields

$$\varrho _ { k } = \arg \min _ { \varrho } \Big \{ \min _ { M \in \Pi ( \varrho _ { k - 1 }, \varrho ) } \frac { 1 } { 2 } \langle C _ { k }, M \rangle + \epsilon H ( M ) \\ + h \langle \psi _ { k - 1 } + \beta ^ { - 1 } \log \varrho, \varrho \rangle \Big \}, \quad ( 2 7 ) \text{ \quad Setting }!$$

where glyph[epsilon1] &gt; 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¨ odinger 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 λ 0 , λ 1 ∈ R N are Lagrange multipliers associated with the equality constraints in (26), then the optimal coupling matrix M opt := [ m opt ( i, j )] N j i, =1 in (27) has the Sinkhorn form

$$m ^ { \text{opt} } ( i, j ) = \exp \left ( \lambda _ { 0 } ( i ) h / \epsilon \right ) \exp \left ( - C _ { k } ( i, j ) / ( 2 \epsilon ) \right ) \\ \exp \left ( \lambda _ { 1 } ( j ) h / \epsilon \right ). \quad ( 2 8 )$$

Since the objective in (27) is proper convex and lower semicontinuous in glyph[rho1] , the strong duality holds, and we consider the Lagrange dual of (27) given by:

$$\lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \lambda \ \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \alpha \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \binom { 1 } { 2 9 } } \\ \mu \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \ \lambda \lambda \ \lambda \lambda \lambda \lambda \lambda \lambda }$$

where

$$F ^ { * } ( \cdot ) \coloneqq \sup _ { \vartheta } \left \{ \langle \cdot, \vartheta \rangle \, - \, F ( \vartheta ) \right \}$$

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.

## A. Proximal Recursions

Given the vectors glyph[rho1] k -1 , ψ k -1 , the matrix C k , and the positive scalars β, h, glyph[epsilon1] in (29), let

$$y \coloneqq \exp ( \lambda _ { 0 } h / \epsilon ), \quad z \coloneqq \exp ( \lambda _ { 1 } h / \epsilon ), \quad \quad ( 3 1 )$$

$$\Gamma _ { k } \coloneqq \exp ( - C _ { k } / 2 \epsilon ), \ \xi _ { k - 1 } \coloneqq \exp ( - \beta \psi _ { k - 1 } - 1 ). \quad ( 3 2 )$$

The following result establishes a system of nonlinear equations for computing λ opt , λ opt in (29), and consequently glyph[rho1] k in (27).

Theorem 1. The vectors λ 0 , λ 1 for y z , ∈ R N from the following system of equations:

0 1 opt opt in (29) can be found by solving

$$y \odot ( \mathbf P _ { k } z ) = \rho _ { k - 1 }, \quad \quad \quad$$

$$z \odot \left ( \Gamma _ { k } ^ { \top } y \right ) = \xi _ { k - 1 } \od$$

and then inverting the maps (31). Let ( y opt , z opt ) be the solution of (33). The vector glyph[rho1] k in (27), i.e., the proximal update (Fig. 1) can then be obtained as

$$\varrho _ { k } = z ^ { \text{opt} } \odot \left ( { \mathbf \Gamma _ { k$$

Proof. From (15), the 'discrete free energy' is

$$F ( \varrho ) = \langle \psi + \beta ^ { - 1 } \log \varrho, \varrho \$$

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

$$F ^ { * } ( \lambda ) = \sup _ { \varrho } \{ \lambda ^ { \top } \varrho - \$$

Setting the gradient of the objective function in (35) w.r.t. glyph[rho1] to zero, and solving for glyph[rho1] yields

$$\varrho _ { \max } = \exp ( \beta ( \lambda - \psi ) - 1 ). \quad \quad ( 3 6 )$$

Substituting (36) back into (35), results

$$F ^ { * } ( \lambda ) = \beta ^ { - 1 } \mathbf 1 ^ { \top } \exp ( \beta ( \lambda - \psi ) - \mathbf 1 ). \quad \quad ( 3 7 )$$

Fixing λ 1 , and taking the gradient of the objective in (29) w.r.t. λ 0 , gives (33a). Likewise, fixing λ 0 , and taking the gradient of the objective in (29) w.r.t. λ 1 gives

$$\nabla _ { \lambda _ { 1 } } F ^ { * } ( - \lambda _ { 1 } ) = z \odot \left ( \mathbf r _ { k } \, ^ { \top } y \right ). \quad \quad ( 3 8 )$$

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

$$\varrho _ { k } = ( M ^ { \text{opt} } ) ^ { \top } 1 = \sum _ { j = 1 } ^ { N } m ^ { \text{opt} } ( j, i ) = z ( i ) \sum _ { j = 1 } ^ { N } \Gamma _ { k } ( j, i ) y ( j ), \\ \text{which is equal to } z \odot \Gamma _ { k } ^ { \top } y, \text{ 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.

## B. Algorithm

1) 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 glyph[rho1] k -1 as input and returns the proximal update glyph[rho1] k as output for k ∈ N . In addition to the data glyph[rho1] k -1 , ψ k -1 , C k , β, h, glyph[epsilon1], N , the Algorithm 1 requires two parameters as user input: numerical tolerance δ , and maximum number of iterations L . The computation in Algorithm 1, as presented, involves making an initial guess for the vector z and then updating y and z until convergence. The iteration over index glyph[lscript] ≤ L is performed while keeping the physical time 'frozen'.

## Algorithm 1 Proximal recursion to compute glyph[rho1] k from glyph[rho1] k - 1

1:

procedure

PROXRECUR(

glyph[rho1]

k

-

1

,

ψ

k

-

1

,

C

k

,

β

,

h

,

glyph[epsilon1]

,

N δ L

,

,

)

2:

Γ

k

←

exp(

-

C

k

/ glyph[epsilon1]

2 )

3:

ξ

←

exp(

-

β

ψ

k

-

1

-

1

)

4:

z

0

←

rand

N

×

1

glyph[triangleright] initialize

5:

z

←

[

z

0

,

0

N

×

(

L

-

1)

]

6:

y

←

[

glyph[rho1]

k

-

1

glyph[circledivide]

(

Γ

k

z

0

)

,

0

N

×

(

L

-

1)

]

7:

glyph[lscript]

= 1

glyph[triangleright] iteration index

$$\stackrel {.. } { 8 } \colon \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$ \sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \sim } { \text{$\sim } - \ell$} \quad & \stackrel { \left | \end{ \end{ \end{pmatrix} }$$

11:

if

‖

y

(:

, glyph[lscript]

+1)

-

y

(:

, glyph[lscript]

)

‖

&lt; δ

&amp;

‖

z

(:

, glyph[lscript]

+1)

-

z

(:

, glyph[lscript]

)

‖

&lt; δ

then

glyph[triangleright] error within tolerance

12:

break

13:

else

$$\iota 4 \colon \quad \ell \leftarrow \ell + 1$$

15:

end if

16: end while

17:

return

glyph[rho1]

k

←

z

(:

, glyph[lscript]

)

glyph[circledot]

(

Γ

glyph[latticetop]

k

y

(:

, glyph[lscript]

)

)

glyph[triangleright] proximal update

18:

end procedure

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

2) Overall algorithm: Samples from the known initial joint PDF ρ 0 are generated as point cloud { x i 0 , glyph[rho1] i 0 } N i =1 . Then for k ∈ N , the point clouds { x i k , glyph[rho1] i k } N i =1 are updated as shown in Fig. 2. Specifically, the state vectors are updated via Euler-Maruyama scheme

Fig. 2: Schematic of the proposed algorithmic setup for propagating the joint state PDF as probability weighted scattered point cloud { x i k , glyph[rho1] i k } N i =1 . The location of the points { x i k } N i =1 are updated via EulerMaruyama scheme; the corresponding probability weights are updated via Algorithm 1. The dashed arrow shown above is present only when the state dynamics is density dependent, as in (19).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-9ed79a12-9891-483b-a980-97ea758c3167.jpg)

applied to the underlying SDE (13) or (19). Explicitly, the EulerMaruyama update for (19) is

$$x _ { k } ^ { i } = x _ { k - 1 } ^ { i } - h \nabla & \left ( \psi ( x _ { k - 1 } ^ { i } ) + \omega ( x _ { k - 1 } ^ { i } ) \right ) \\ & + \sqrt { 2 \beta ^ { - 1 } } \left ( w _ { k } ^ { i } - w _ { k - 1 } ^ { i } \right ), \quad \quad ( 3 9 )$$

where ω ( ) := · ∫ φ ( · -y ) ρ ( y )d y , and w i k := w i ( t = kh ) , k ∈ N . The same for (13) is obtained by setting φ ≡ ω ≡ 0 in (39).

The corresponding probability weights { glyph[rho1] i k } N i =1 are updated via Algorithm 1. Notice that computing C k requires both { x i k -1 } N i =1 and { x i k } N i =1 , and that 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 -∇ ψ in (13), or - ∇ ( ψ + ∇ ( ρ ∗ φ )) 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 EulerMaruyama 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].

## 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 z z , ˜ ∈ K , where K is a non-empty open convex cone. Further, suppose that K is a normal cone, i.e., there exists constant α such that ‖ z ‖≤ α ‖ ˜ z ‖ for

Fig. 3: Comparison of the analytical and proximal solutions of the FPK PDE for (45) with time step h = 10 -3 , and with parameters a = 1 , β = 1 , glyph[epsilon1] = 5 × 10 -2 . Shown above are the time evolution of the ( left ) PDFs, ( middle ) means, and ( right ) variances.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-857b1f17-6c75-4880-bdd0-5448001a02b9.jpg)

z ≤ ˜ z . Thompson [42] proved that K is a complete metric space w.r.t. the so-called Thompson metric given by

$$d _ { T } \, ( z, \widetilde { z } ) \coloneqq \max \{ \log \gamma ( z / \widetilde { z } ), \log \gamma ( \widetilde { z } / z ) \},$$

where γ ( z z / ˜ ) := inf { c &gt; 0 | z ≤ c ˜ z } . In particular, if K ≡ R N &gt; 0 (positive orthant of R N ), then

$$d _ { \text{T} } \left ( z, \widetilde { z } \right ) = \log \max \left \{ \max _ { i = 1, \dots, N } \left ( \frac { z _ { i } } { \widetilde { z } _ { i } } \right ), \max _ { i = 1, \dots, N } \left ( \frac { \widetilde { z } _ { i } } { z _ { i } } \right ) \right \}. \quad ( 4 0 ) \quad \text{homoger}$$

From (25) and (32), C k ( i, j ) ∈ [0 , ∞ ) which implies Γ k ( i, j ) ∈ (0 1] , ; therefore, Γ k is a positive linear map for each k = 1 2 , , . . . . Thus, by (linear) Perron-Frobenius theorem, the map θ 3 is contractive on R N &gt; 0 w.r.t. d T . The map θ 2 is an isometry by Definition 3. As for the map θ 1 , notice that the quantity r := 1 / (1+ βglyph[epsilon1]/h ) ∈ (0 1) , since βglyph[epsilon1]/h &gt; 0 . Therefore, the map θ 1 ( z ) := z r (element-wise exponentiation) is monotone (order preserving) and homogeneous of degree r ∈ (0 1) , on R N &gt; 0 . By Proposition 2, the map θ 1 ( z ) is strictly contractive. Thus, the composition

Proposition 2. [43, Proposition 3.2], [44] Let K be an open, normal, convex cone, and let p : K ↦→ K be an order preserving homogeneous map of degree r ≥ 0 , i.e., p ( c z ) = c r p z ( ) for any c &gt; 0 and z ∈ K . Then, for all z z , ˜ ∈ K , we have

$$d _ { \text{T} } \left ( p ( z ), p ( \widetilde { z } ) \right ) \leq r d _ { \text{T} } \left ( z, \widetilde { z } \right ).$$

In particular, if r ∈ [0 , 1) , then the map p ( ) · is strictly contractive in the Thompson metric d T , and admits unique fixed point in 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

$$z ( \colon, \ell + 1 ) & = \left ( \xi _ { k - 1 } \otimes \left ( \Gamma _ { k } ^ { \top } y ( \colon, \ell ) \right ) \right ) ^ { \frac { 1 } { 1 + \beta \epsilon / h } } \\ & = \left ( \xi _ { k - 1 } \otimes \left ( \Gamma _ { k } ^ { \top } \varrho _ { k - 1 } \otimes ( \Gamma _ { k } z ( \colon, \ell ) ) \right ) \right ) ^ { \frac { 1 } { 1 + \beta \epsilon / h } } \quad ( 4 1 )$$

for glyph[lscript] = 1 2 , , . . . , is strictly contractive in the Thompson metric (40) on R N &gt; 0 , and admits unique fixed point z opt ∈ R N &gt; 0 .

Proof. Rewriting (41) as

$$z ( \colon, \ell + 1 ) = \left ( \left ( \xi _ { k - 1 } \oslash \left ( \Gamma _ { k } ^ { \top } \varrho _ { k - 1 } \right ) \right ) \odot \left ( \Gamma _ { k } z ( \colon, \ell ) \right ) \right ) ^ { \frac { 1 } { 1 + \beta \left ( / h } },$$

↦

and letting η ≡ η k,k +1 := ξ k -1 glyph[circledivide] ( Γ glyph[latticetop] k glyph[rho1] k -1 ) , we notice that iteration (41) can be expressed as a cone preserving composite map θ := θ 1 ◦ θ 2 ◦ θ 3 , where θ : R N &gt; 0 → R N &gt; 0 , given by

$$z ( \colon, \ell + 1 ) = \theta \left ( z ( \colon, \ell ) \right ) = \theta _ { 1 } \circ \theta _ { 2 } \circ \theta _ { 3 } \, \left ( z ( \colon, \ell ) \right ), \quad ( 4 2 ) \, \quad r a n s i e n$$

and θ 1 ( z ) := z 1 1+ βglyph[epsilon1]/h , θ 2 ( z ) := η glyph[circledot] z , θ 3 := Γ k z . Our strategy is to prove that the composite map θ is contractive on R N &gt; 0 w.r.t. the metric d T .

$$\theta = \underset { \text{ strictly contractive} } { \underbrace { \theta _ { 1 } } } \circ \underset { \text{ isometry contractive} } { \underbrace { \theta _ { 2 } } } \circ \underset { \text{ contractive} } { \underbrace { \theta _ { 3 } } }$$

is strictly contractive w.r.t. d T , and (by Banach contraction mapping theorem) admits unique fixed point z opt in R N &gt; 0 .

Corollary 4. The Algorithm 1 converges to unique fixed point ( y opt , z opt ) ∈ R N &gt; 0 × R N &gt; 0 .

Proof. Since y (: , glyph[lscript] + 1) = glyph[rho1] k -1 glyph[circledivide] ( Γ k z (: , glyph[lscript] +1)) , the z iterates converge to unique fixed point z opt ∈ R N &gt; 0 (by Theorem 3), and the linear maps Γ k are contractive (by Perron-Frebenius theory, as before), consequently the y iterates also converge to unique fixed point y opt ∈ R N &gt; 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.

A. Linear Gaussian Systems

For an Itˆ o SDE of the form

$$d x = A x \, d t + B \, d w,$$

it is well known that if x 0 := x ( t = 0) ∼ N ( µ 0 , Σ 0 ) , then the transient joint PDFs ρ ( x , t ) = N ( µ ( ) t , Σ ( )) t where the vectormatrix pair ( µ ( ) t , Σ ( )) t evolve according to the ODEs

$$\overset { \text{gy} } { \Phi } _ { \Phi \wedge } \quad \dot { \mu } ( t ) = A \mu, \quad \mu ( 0 ) = \mu _ { 0 },$$

$$\ u e \quad \dot { \Sigma } ( t ) = A \Sigma ( t ) + A \Sigma ( t ) ^ { \top } + B B ^ { \top }, \ \Sigma ( 0 ) = \Sigma _ { 0 }. \quad ( 4 4 b )$$

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

1) Ornstein-Uhlenbeck Process: We consider the univariate system

$$d x = - a x \, d t + \sqrt { 2 \beta ^ { - 1 } } d w, \ a, \beta > 0, \quad \quad ( 4 5 ) \quad ^ { 2 } \Big | _ { \circ } ^ { * }$$

which is in JKO form (13) with ψ x ( ) = 1 2 ax 2 . We generate N = 400 samples from the initial PDF ρ 0 = N ( µ , σ 0 2 0 ) with µ 0 = 5 and σ 2 0 = 4 × 10 -2 , and apply the proposed proximal recursion for (45) with time step h = 10 -3 , and with parameters a = 1 , β = 1 , glyph[epsilon1] = 5 × 10 -2 . For implementing Algorithm 1, we set the tolerance δ = 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 N ( µ 0 exp( -at , ) ( σ 2 0 -1 aβ ) exp( -2 at ) + 1 aβ ) , and the mean-variance trajectories (computed from the numerical integration of the weighted scattered point cloud data) match with the corresponding analytical solutions.

2) Multivariate LTI: We next consider the multivariate case (43) where the pair ( A B , ) is assumed to be controllable, and the matrix A is Hurwitz (not necessarily symmetric). Under these assumptions, the stationary PDF is N ( 0 Σ , ∞ ) where Σ ∞ 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 A , there does not exist constant symmetric positive definite matrix Ψ such that Ax = -∇ x glyph[latticetop] Ψ 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 changeof-coordinates to (43) with

$$A = \begin{pmatrix} - 1 0 & 5 \\ - 3 0 & 0 \end{pmatrix}, \ \ B = \begin{pmatrix} 2 \\ 2. 5 \end{pmatrix},$$

which satisfy the stated assumptions for the pair ( A B , ) , and implement the proposed proximal recursion on this transformed coordinates with N = 400 samples generated from the initial PDF ρ 0 = N ( µ 0 , Σ 0 ) , where µ 0 = (4 4) , glyph[latticetop] and Σ 0 = 4 I 2 . As before, we set δ = 10 -3 , L = 100 , h = 10 -3 , β = 1 , glyph[epsilon1] = 5 × 10 -2 . Once the proximal updates are done, we transform back the probability weighted scattered point cloud to the original state space coordinates 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 N ( µ ( ) t , Σ ( )) 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 Σ ∞ . As a consequence, the stationary covariance resulting from the proximal recursion oscillates about the true stationary value.

ρ

ρ

x

Fig. 4: Comparison of the analytical ( contour plots ) and proximal ( weighted scattered point cloud ) joint PDFs of the FPK PDE for (43) with time step h = 10 -3 , and with parameters β = 1 , glyph[epsilon1] = 5 × 10 -2 . Simulation details are given in Section IV-A2. The color ( red = high, blue = low ) denotes the joint PDF value obtained via proximal recursion at a point at that time (see colorbar).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-9d09388f-c7ac-41be-9394-4a773db21460.jpg)

4

t

Fig. 5: Comparison of the components of the mean vectors from analytical ( dashed ) and proximal ( solid ) computation of the joint PDFs for (43) with time step h = 10 -3 , and with parameters β = 1 , glyph[epsilon1] = 5 × 10 -2 . Simulation details are given in Section IV-A2.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-4bc01826-d6d2-4071-9c7b-7ea00c1471a7.jpg)

## B. Nonlinear Non-Gaussian System

Next we consider a planar nonlinear system of the form (13) with ψ x , x ( 1 2 ) = 1 4 (1+ x 4 1 )+ 1 2 ( x 2 2 -x 2 1 ) (see Fig. 7). As mentioned in Section II, the stationary PDF is ρ ∞ ( x ) = κ exp( -βψ ( x )) , which for our choice of ψ , 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 ρ 0 = N ( µ 0 , Σ 0 ) with µ 0 = (2 2) , glyph[latticetop] and Σ 0 = 4 I 2 , and set δ = 10 -3 , L = 100 , h = 10 -3 , β = 1 , glyph[epsilon1] = 5 × 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

Fig. 6: Comparison of the components of the covariance matrices from analytical ( dashed ) and proximal ( solid ) computation of the joint PDFs for (43) with time step h = 10 -3 , and with parameters β = 1 , glyph[epsilon1] = 5 × 10 -2 . Simulation details are given in Section IV-A2.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-efa83e48-63af-4595-8328-fc43288d3c4a.jpg)

t

=2 Fig. 7: The drift potential ψ x , x ( 1 2 ) = 1 4 (1 + x 4 1 ) + 1 2 ( x 2 2 -x 2 1 ) used in the numerical example given in Section IV-B.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-44e6be13-72cb-4215-9f91-eabd27c73fc2.jpg)

computed via the proximal recursion, tend to the known stationary solution ρ ∞ (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 glyph[lscript] ≤ 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.

## 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 densitydependent sample path dynamics (19). If both the potentials ψ and

10 20 30 40 50 Fig. 8: The proximal ( weighted scattered point cloud ) joint PDFs of the FPK PDE (14) with the drift potential shown in Fig. 7, time step h = 10 -3 , and with parameters β = 1 , glyph[epsilon1] = 5 × 10 -2 . Simulation details are given in Section IV-B. The color ( red = high, blue = low ) denotes the joint PDF value obtained via proximal recursion at a point at that time (see colorbar). In the bottom right plot, the contour lines correspond to the analytical solution for the stationary PDF ρ ∞ .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-c7722097-04f1-41e5-8ed2-8355b2eb3513.jpg)

Fig. 9: The computational times for proximal updates for the simulation in Section IV-B. Here, the physical time-step h = 10 -3 s, and k ∈ N .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-cffe51cb-4fc7-4cdd-b5ef-e250a02869fd.jpg)

φ are convex functions, then (18) admits unique stationary density ρ ∞ ( x ) (see e.g., [31], [45, Section 5.2]). To keep the exposition simple, we consider the univariate case ψ x ( ) = 1 2 ax 2 , φ x ( ) = 1 2 bx 2 , a, b &gt; 0 . Performing the integration appearing in the convolution allows us to rewrite (19) as the mean-reverting process:

$$\mathrm d x = - \left ( ( a + b ) x - b \mu ( t ) \right ) \, \mathrm d t \, + \, \sqrt { 2 \beta ^ { - 1 } } \, \mathrm d w, \ \ x ( 0 ) = x _ { 0 }, \ ( 4 6 )$$

where µ t ( ) is the mean of the transient PDF ρ x, t ( ) . Assuming x 0 ∼ N ( µ , σ 0 2 0 ) , and applying expectation operator to both sides of (46) yields µ t ( ) = µ 0 exp( -at ) . Consequently, the transient PDF for (46)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-13357d79-f714-4df6-b75e-5a04db2b7f67.jpg)

σ

-

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-71c89910-fbd7-44bd-ae27-687facc4d51f.jpg)

t

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-7f9a7189-aec8-49f5-bac7-91d771d032ac.jpg)

t

Fig. 10: Comparison of the analytical and proximal solutions of the McKean-Vlasov flow for (46) with time step h = 10 -3 , ρ 0 = N (5 , 9) , and with parameters a = b = 1 , β = 1 , glyph[epsilon1] = 5 × 10 -2 . Shown above are the time evolution of the transient ( left ) PDFs, ( middle ) means, and ( right ) variances.

at time t is ρ x, t ( ) = N ( µ t , σ ( ) 2 ( ) t ) , with

$$\mu ( t ) & = \mu _ { 0 } \exp ( - a t ), & ( 4 7 a ) \quad _ { \text{two active} } ^ { \text{s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.s.}$$

$$\sigma ^ { 2 } ( t ) = \left ( \sigma _ { 0 } ^ { 2 } - \frac { 1 } { ( a + b ) \beta } \right ) \exp \left ( - 2 ( a + b ) t \right ) + \frac { 1 } { ( a + b ) \beta }. \ ( 4 7 b ) \ \stackrel { \text{ two case} } { \text{depend} } \\ \text{conserva}$$

Clearly, the stationary PDF is ρ ∞ ( x ) = N (0 1 , / a ( + ) b β ) .

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 ∫ R n × R n φ ( x -y ) glyph[rho1] ( x ) glyph[rho1] ( y )d x y d in (22b), with the linear term ∫ R n φ ( x -y ) glyph[rho1] ( x ) glyph[rho1] k -1 ( y )d x y d , k ∈ N , resulting in a semi-implicit variant of (22), given by

$$\varrho _ { k } = \underset { \varrho \in \mathcal { D } _ { 2 } } { \arg \inf } \, \frac { 1 } { 2 } W ^ { 2 } ( \varrho _ { k - 1 }, \varrho ) + h \, F ( \varrho _ { k - 1 }, \varrho ), \ k \in \mathbb { N }, \quad ( 4 8 ) \quad \text{model } [ 4 \\ d x =$$

with glyph[rho1] 0 ≡ ρ 0 ( x ) (the initial PDF). The proof for the fact that such a semi-implicit scheme guarantees glyph[rho1] k ( x ) → ρ ( x , t = kh ) for h ↓ 0 , where ρ ( 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 ( glyph[rho1] ) = 〈 ψ k -1 + β -1 log glyph[rho1] , glyph[rho1] 〉 . The scheme (48) allows us to write a similar expression in the McKean-Vlasov gradient flow case, as F ( glyph[rho1] ) = 〈 ψ k -1 + D k -1 glyph[rho1] k -1 + β -1 log glyph[rho1] , glyph[rho1] 〉 , where the symmetric matrix D k -1 is given by

$$D _ { k - 1 } ( i, j ) = \phi \left ( x _ { k - 1 } ^ { i } - x _ { k - 1 } ^ { j } \right ), \ \ i, j = 1, \dots, N, \ \ k \in \mathbb { N }.$$

For the particular choice φ x ( ) = 1 2 bx 2 , notice that D k -1 is a (scaled) Euclidean distance matrix on the Euler-Maruyama update. In this case, the associated Euler-Maruyama scheme (39) is

$$x _ { k } ^ { i } \ & = \ x _ { k - 1 } ^ { i } - h \left ( a x _ { k - 1 } ^ { i } + b \left ( x _ { k - 1 } ^ { i } - 1 ^ { \top } \varrho _ { k - 1 } \right ) \right ) \quad \begin{array} { c } \text{stationar} \\ x > 0. \end{array} \\ & + \quad \sqrt { 2 \beta ^ { - 1 } } \left ( w _ { k } ^ { i } - w _ { k - 1 } ^ { i } \right ), \quad k \in \mathbb { N }, \, i = 1, \dots, N, ( 4 9 ) \quad \quad \text{In}$$

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 N ( µ t , σ ( ) 2 ( )) t given by (47). For Fig. 10, the parameter values used in our simulation are h = 10 -3 , β = 1 , glyph[epsilon1] = 5 × 10 -2 , µ 0 = 5 , σ 2 0 = 9 , N = 400 , δ = 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 statedependent 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.

## A. Multiplicative Noise

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

$$\Delta ( 4 8 ) \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \mu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta ^ { 2 } ( 4 8 ) \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nol } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta _ { \nu } \Delta ( 5 0 )$$

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

$$\stackrel { \bullet \text{discrete} } { e \, \text{ scheme} } \\ \stackrel { n - V | \text{asov} } { \log \varrho, \varrho \rangle, } \quad \rho ( x, t ) = \begin{cases} c \exp ( - ( u + v ) ) \left ( \frac { v } { u } \right ) ^ { q / 2 } I _ { q } \left ( 2 \sqrt { u v } \right ), & x > 0, \\ 0, & \text{otherwise}, \end{cases} \quad ( 5 1 )$$

where I q ( ) · is the modified Bessel function of order q , and

$$q & \coloneqq \frac { 2 a \theta } { b ^ { 2 } } - 1, \quad & c \coloneqq \frac { 2 a } { b ^ { 2 } \left ( 1 - \exp ( - a t ) \right ) }, \quad & ( 5 2 a ) \\ \text{is}$$

The transient solutions (51) are non-central chi-squared PDFs. The stationary solution is a Gamma PDF ρ ∞ ( x ) ∝ x q exp ( -2 ax/b 2 ) , x &gt; 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 = ς ( x ) such that the SDE for new variable y has unity diffusion coefficient. From Itˆ o's lemma, it follows that the requisite transformation for (50) is y = ς ( x ) := 2 b √ x , and the resulting SDE in y becomes

$$\mathrm d y = \left \{ \, \left ( \frac { 2 a \theta } { b ^ { 2 } } - \frac { 1 } { 2 } \right ) \frac { 1 } { y } - \frac { a } { 2 } y \right \} \mathrm d t \, + \, \mathrm d w. \quad \quad ( 5 3 )$$

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

$$\psi ( y ) = a y ^ { 2 } / 4 - ( q + 1 / 2 ) \log y.$$

So the proposed proximal algorithm (Algorithm 1) can be applied to (53), and at each time t , the resulting PDF ρ Y ( y, t ) can be

Fig. 11: Comparison of the analytical and proximal transient PDFs of the FPK PDE for (50) with time step h = 10 -3 , and with parameters a = 3 , b = 2 , θ = 2 , x 0 = 5 , glyph[epsilon1] = 5 × 10 -2 . To approximate the analytical PDFs resulting from ρ 0 ( x ) = δ x ( -5) , the proximal recursions were performed with the initial PDF N (5 10 , -4 ) .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-e5da1fec-60c5-40f0-9ce1-6bf15a082b37.jpg)

transformed back to the PDF ρ X ( x, t ) via the change-of-measure formula for the push-forward map x = ς -1 ( y ) = by 2 / 4 . Fig. 11 shows the comparison of the analytical and proximal transient PDFs for (50) resulting from ρ 0 = δ x ( -5) with a = 3 , b = 2 , θ = 2 , h = 10 -3 , glyph[epsilon1] = 5 × 10 -2 , N = 400 , δ = 10 -3 , and L = 100 .

## B. Mixed Conservative-Dissipative Drift

In many engineering applications, one encounters Itˆ o 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

$$\stackrel { \sigma \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \colon \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \, \text{out appl} } { \ddot { q } } = - \frac { \mu q } { \| q \, \| _ { 2 } ^ { 3 } } \, - \, \gamma \dot { q } \, + \, f _ { \text{pert} } ( q ) & \text{conserva} \\ & + \, \sqrt { 2 \beta ^ { - 1 } \gamma } \times \text{ stochastic forcing,} \quad (55) & \varrho _ { k }$$

where q := ( x, y, z ) glyph[latticetop] ∈ R 3 is the relative position vector for the satellite, µ is a constant (product of the Gavitational constant and the mass of Earth), -γ ˙ q models linear drag † , f pert ( 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τ := cos τ , sτ := sin τ , and recalling the relations among spherical coordinates ( r, θ, φ ) and cartesian coordinates ( x, y, z ) , given by r = √ x 2 + y 2 + z 2 , cφ = x/r , sφ = √ 1 -( cφ ) 2 ,

† More generally, for nonlinear drag of the form - ∇ γ ˙ q ˜ V ( ˙) q , the term -γ p in (59) will become - ∇ γ ˜ V ( p ) . This will entail modifying the functional ̂ F in (61) as E ρ [ ˜ V ( p ) + β -1 log ρ ] . The linear drag illustrated here corresponds to the special case ˜ V ( p ) ≡‖ p ‖ 2 2 / 2 .

cθ = z/r , sθ = √ 1 -( cθ ) 2 , we can write f pert ( q ) in spherical coordinates as [51, eqn. (12)-(13)]

$$\begin{pmatrix} & \cdots \dots \dots \dots \dots \dots \dots \dots \dots \end{pmatrix} \\ \Big | \quad & \begin{pmatrix} f _ { r } \\ f _ { \theta } \\ f _ { \phi } \end{pmatrix} _ { \text{pert} } = \begin{pmatrix} \frac { k } { 2 r ^ { 4 } } \left ( 3 ( s \theta ) ^ { 2 } - 1 \right ) \\ - \frac { k } { r ^ { 5 } } s \theta \, c \theta \end{pmatrix} _ { \text{pert} }, \, k \coloneqq 3 J _ { 2 } R _ { E } ^ { 2 } \mu = \, \text{constant, } ( 5 6 ) \\ \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \ dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots 0 \end{pmatrix}$$

and the same in cartesian coordinates as

$$f _ { \text{pert} } ( q ) = \begin{pmatrix} f _ { x } \\ f _ { y } \\ f _ { z } \end{pmatrix} _ { \text{pert} } = \begin{pmatrix} s \theta \, c \phi & c \theta \, c \phi & - s \phi \\ s \theta \, s \phi & c \theta \, s \phi & c \phi \\ c \theta & - s \theta & 0 \end{pmatrix} \begin{pmatrix} f _ { r } \\ f _ { \theta } \\ f _ { \phi } \end{pmatrix} _ { \text{pert} }. \ \ ( 5 7 )$$

In (56), the Earth oblateness coefficient J 2 = 1 082 . × 10 -3 , the radius of Earth R E = 6 3781 . × 10 6 m, and the Earth standard Gravitational parameter µ = 3 9859 . × 10 14 m /s . Modeling the stochastic forcing 3 2 in (55) as standard Gaussian white noise (as in [51]), (55) can then be expressed as an Itˆ o SDE in R 6 :

$$\text{$\quad \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \
\text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$	\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\	\text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \cast \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \cast	\text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad   \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \hat \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \gefrak { \quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad    \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text$\quad \text$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text$\quad \text{$\quad \text{$\quad \text$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\	\text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\	\text{$\quad \text$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text$\quad \text{$\quad \text{$\quad \text{$\quad \text$\quad \text{$\quad \text{$\quad \text$\quad \text{$\quad \text{$\quad \text{$\quad \text$\quad \text$\quad \text{$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text{$\quad \text{$\quad \text{$\quad \text{$\quad \text$\quad \text{$\quad \text{$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text{$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \tilde { \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\	 \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\quad \text$\$$ \, \quad \end{array} }$$

or more succinctly,

$$\binom { \text{d} q } { \text{d} p } = \binom { p } { - \nabla V ( q ) - \gamma p } \, \text{d} t \, + \, \sqrt { 2 \beta ^ { - 1 } \gamma } \, \binom { 0 _ { 3 \times 1 } } { \text{d} w _ { 3 \times 1 } } \,, \ \ ( 5 9 )$$

where p := ( v x , v y , v z ) glyph[latticetop] ∈ R 3 is the velocity vector, and

$$\overbrace { v e s } _ { a n d } \quad V ( q ) \coloneqq V _ { \text{gravitational} } ( q ) + V _ { \text{per} } ( q ),$$

$$\text{and} \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \partial \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odot \quad \odots \quad \text{$\quad $-} \\ \text{ixed} \quad \text{$\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\ \text{$\quad $-} \\ \text{$\quad $f_{z}$} \quad \text{$\quad $f_{z}$} \quad \text{$\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $
\text{$\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad \text{$\quad $-}$$

Introducing a 'Hamiltonian-like function' H ( q p , ) := ‖ p ‖ 2 2 / 2 + V ( q ) , one can verify that the stationary PDF for (59) is ρ ∞ ∝ exp( -βH ( q p , )) , and that the 'total free energy' F ρ ( ) := E ρ [ H + β -1 log ρ ] serves as Lyapunov functional for the associated FPK PDE, i.e., d F d t &lt; 0 . However, the proximal recursion (22b) does not apply as is, instead needs to be modified to account the joint conservative-dissipative effect as

$$5 ) \quad \varrho _ { k } = \underset { \varrho \in \mathcal { D } _ { 2 } } { \arg \inf } \ \frac { 1 } { 2 } \widehat { W } _ { h } ^ { 2 } ( \varrho _ { k - 1 }, \varrho ) + h \gamma \ \widehat { F } ( \varrho ), \ \ k \in \mathbb { N }. \quad ( 6 1 )$$

We implement a recursion from [52, Scheme 2b] where ̂ F glyph[rho1] ( ) := E ρ [ 1 2 ‖ p ‖ 2 2 + β -1 log ρ ] , and ̂ W 2 h ( glyph[rho1] k -1 , glyph[rho1] ) is the optimal mass transport cost (as in squared 2-Wasserstein metric (10)) with modified cost function (modified integrand in (10)), i.e.,

$$\underset { \underline { r }, \, y, \, z } { \cosines } \quad \widehat { W } _ { h } ^ { 2 } ( \rho _ { 1 }, \rho _ { 2 } ) \colon = \inf _ { d \pi \in \Pi ( \pi _ { 1 }, \pi _ { 2 } ) } \int _ { \mathcal { X } \times \mathcal { Y } } \widehat { s } _ { h } \left ( q, p ; \widetilde { q }, \widetilde { p } \right ) \, d \pi \left ( q, p, \widetilde { q }, \widetilde { p } \right ), \\ \underset { \underline { r }, \, y, \, z \right ) } { \longrightarrow } \quad \dots \, \int _ { \mathcal { W } \times \mathcal { Y } \, \dots \, \int \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \$$

where ( q p , ) glyph[latticetop] and ( ˜ ˜ q p , ) glyph[latticetop] are two realizations of the state vector governed by (59)-(60); the integrand in (62) is

$$\hat { s } _ { h } \left ( q, p ; \widetilde { q }, \widetilde { p } \right ) \coloneqq & \ \| \widetilde { p } - p + h \nabla V ( q ) \, \| _ { 2 } ^ { 2 } \\ & + 1 2 \left \| \frac { \widetilde { q } - q } { h } - \frac { \widetilde { p } + p } { 2 } \right \| _ { 2 } ^ { 2 }. \quad \ \ ( 6 3 )$$

That the proximal recursion (61) with the above choices of functionals ̂ F and ̂ W h guarantees glyph[rho1] k ( q p , , h ) → ρ ( q p , , t = kh ) for k ∈ N as h ↓ 0 , where ρ ( q 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

$$\psi _ { k - 1 } ( i ) \ \coloneqq \ \frac { 1 } { 2 } ( p _ { k - 1 } ^ { i } ) ^ { \top } p _ { k - 1 } ^ { i }, \ \ i = 1, \dots, N, \quad ( 6 4 ) \quad \begin{array} { c } 0. 2 0 \\ 0. 0 0 \\ \end{array} \Big | \\ C _ { k } ( i, j ) \ \coloneqq \ \widehat { s } _ { h } ( q _ { k } ^ { i }, p _ { k } ^ { i } ; q _ { k - 1 } ^ { j }, p _ { k - 1 } ^ { j } ), \ \ i, j = 1, \dots, N, ( 6 5 ) \end{array}$$

respectively. Notice that ̂ s h in (63) (and consequently, ̂ W h in (62)) is not a metric (in particular, non-symmetric). Thus, the matrices C k in (65) for k ∈ N are not symmetric.

To apply the proposed framework, we first non-dimensionalize the variables q (m), p (m/s), t (s), w ( √ s) in (59) as

$$q ^ { \prime } = q / R, \ \ p ^ { \prime } = p / ( R / T ), \ \ t ^ { \prime } = t / T, \ \ w ^ { \prime } = w / \sqrt { T }, \ \ ( 6 6 ) \quad \bigoplus _ { \infty } ^ { 1 } \Big | _ { \infty }$$

where R := 4 2164 . × 10 7 mis 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ˆ o's lemma on (59), the non-dimensional SDE in ( q ′ , p ′ ) glyph[latticetop] ∈ R 6 becomes

$$d q ^ { \prime } & = p ^ { \prime } \, d t ^ { \prime }, & ( 6 7 a ) \, \frac { \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \ tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde{ \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \  } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ }
d p ^ { \prime } & = \left \{ - \frac { T \mu } { R ^ { 3 } } \frac { q ^ { \prime } } { \| \tilde { \ } \cdot \, q ^ { \prime } \| \tilde { \ } \cdot \, \tilde { \ } \cdot \, q ^ { \prime } \| \tilde { \ } \cdot \, q ^ { \prime } \| \tilde { \ } \cdot \, q ^ { \prime } \| \tilde { \ } \cdot \, q ^ { \prime } \| \tilde { \ } \cdot \, q ^ { \prime } \| \tilde { \ } \cdot \, q ^ { \prime } \| \tilde { \ } \cdot \, q ^ { \prime } \| \tilde { \ } \cdot \, q ^ { \prime } \| \tilde { \ } \cdot \, q ^ { \prime } \| \tilde { \ } \cdot \, q ^ { \prime } \| \tilde { \ } \cdot \, q ^ { \prime } \| \tilde { \ } \cdot } \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \delta \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \ \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \cdot \, \tilde { \ } \vdots \, \, \,$$

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 ( q ′ , p ′ ) glyph[latticetop] to the same in original variables ( q p , ) glyph[latticetop] 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]).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-841b412d-bf9f-444a-a495-5efa72e5be44.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-dd57fd1b-af9f-4c08-ad90-03d574465bdb.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-c567fc91-0896-40ac-8afe-bb05b5e09f19.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-16331648-8dae-4fa1-9bc9-04d938885418.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-536787e8-945b-4602-a267-2cef4380ca73.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-e7882ef4-77d7-47bb-a422-b1c8c7db4960.jpg)

Fig. 12: Univariate marginal PDFs at t = 0 005 . s for (59)-(60) computed from the joint PDF at that time obtained via the proposed proximal algorithm for recursion (61) with time step h = 10 -5 , and with parameters β = 1 m /s 2 2 , γ = 1 s -1 , glyph[epsilon1] = 5 × 10 -2 , δ = 10 -3 , L = 100 , and N = 400 .

Fig. 13: Univariate marginal PDFs at t = 0 01 . s for (59)-(60) computed from the joint PDF at that time obtained via the proposed proximal algorithm for recursion (61) with time step h = 10 -5 , and with parameters β = 1 m /s 2 2 , γ = 1 s -1 , glyph[epsilon1] = 5 × 10 -2 , δ = 10 -3 , L = 100 , and N = 400 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-4f93e2f4-f629-4fa8-8224-2280f2ad4eaf.jpg)

Fig. 14: The computational times needed for proximal updates in the 6-state numerical example reported in Section V-B. Here, the physical time-step h = 10 -5 s, and k ∈ N .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00533v2-figure-1a094802-6ef6-485d-9e6f-a524572193dc.jpg)

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 R 6 at those times, computed through the proposed proximal algorithm. For this simulation, the initial joint PDF ρ 0 ( q ′ , p ′ ) = N ( µ 0 , Σ 0 ) , where µ 0 = (1 0 0 0 0 0) , , , , , glyph[latticetop] , and Σ 0 = 10 -4 × diag (3 335 . , 6 133 3 933 6 562 9 246 5 761) . , . , . , . , . . The parameter values used in the simulation are: h = 10 -5 , β = 1 m /s , 2 2 γ = 1 s -1 , glyph[epsilon1] = 5 × 10 -2 , δ = 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.

## APPENDIX

- A. Derivation of (28):

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

$$\rightharpoonup _ { k } \rightharpoonup \sum _ { \varrho } \sum _ { \varrho } \left \{ \min _ { M \in \Pi ( \varrho _ { k - 1 }, \varrho ) } \frac { 1 } { 2 h } \langle C _ { k }, M \rangle + \frac { \epsilon } { h } H ( M ) & \frac { \infty } { 1 6 9 } \\ - \frac { \epsilon } { h } 1 ^ { \top } M 1 + F ( \varrho ) \right \}, \quad ( 6 8 ) & \frac { [ 1 4 ] } { \text{met} \\ C o n } )$$

since 1 glyph[latticetop] M 1 = 1 . The Lagrangian L associated with the inner minimization in (68) is given by

$$\mathcal { L } = \frac { 1 } { 2 h } & \langle C _ { k }, M \rangle + \frac { \epsilon } { h } H ( M ) - \frac { \epsilon } { h } \mathbf 1 ^ { \top } M \mathbf 1 + F ( \varrho ) & \cdots _ { A u } \\ & \quad + \langle \lambda _ { 0 }, M \mathbf 1 - \varrho _ { k - 1 } \rangle + \langle \lambda _ { 1 }, M ^ { \top } \mathbf 1 - \varrho \rangle. \quad ( 6 9 ) & \frac { [ 1 6 ] } { \intertext { \text{for} } }$$

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

## REFERENCES

- [1] H. Risken, Fokker-Planck equation: Methods of solution and applications . Springer, 1996.
- [2] M. Ehrendorfer, 'The Liouville equation and its potential usefulness for the prediction of forecast skill. part I: Theory,' Monthly Weather Review , vol. 122, no. 4, pp. 703-713, 1994.
- [3] A. Halder and R. Bhattacharya, 'Beyond Monte Carlo: A computational framework for uncertainty propagation in planetary entry, descent and landing,' in AIAA Guidance, Navigation, and Control Conference , 2010, p. 8029.
- [4] --, 'Dispersion analysis in hypersonic flight during planetary entry using stochastic Liouville equation,' Journal of Guidance, Control, and Dynamics , vol. 34, no. 2, pp. 459-474, 2011.
- [5] S. Hess, 'Fokker-Planck-equation approach to flow alignment in liquid crystals,' Zeitschrift f¨r Naturforschung A u , vol. 31, no. 9, pp. 1034-1037, 1976.
- [6] W. Muschik and B. Su, 'Mesoscopic interpretation of Fokker-Planck equation describing time behavior of liquid crystal orientation,' The Journal of Chemical Physics , vol. 107, no. 2, pp. 580-584, 1997.
- [7] Y. P. Kalmykov and W. T. Coffey, 'Analytical solutions for rotational diffusion in the mean field potential: application to the theory of dielectric relaxation in nematic liquid crystals,' Liquid crystals , vol. 25, no. 3, pp. 329-339, 1998.
- [8] W. Park, J. S. Kim, Y. Zhou, N. J. Cowan, A. M. Okamura, and G. S. Chirikjian, 'Diffusion-based motion planning for a nonholonomic flexible needle model,' in Robotics and Automation, 2005. ICRA 2005. Proceedings of the 2005 IEEE International Conference on . IEEE, 2005, pp. 4600-4605.
- [9] W. Park, Y. Liu, Y. Zhou, M. Moses, and G. S. Chirikjian, 'Kinematic state estimation and motion planning for stochastic nonholonomic systems using the exponential map,' Robotica , vol. 26, no. 4, pp. 419-434, 2008.
- [10] H. Hamann and H. W¨ orn, 'A framework of space-time continuous models for algorithm design in swarm robotics,' Swarm Intelligence , vol. 2, no. 2-4, pp. 209-239, 2008.
- [11] S. Challa and Y. Bar-Shalom, 'Nonlinear filter design using FokkerPlanck-Kolmogorov probability density evolutions,' IEEE Transactions on Aerospace and Electronic Systems , vol. 36, no. 1, pp. 309-315, 2000.
- [12] F. Daum, 'Nonlinear filters: beyond the Kalman filter,' IEEE Aerospace and Electronic Systems Magazine , vol. 20, no. 8, pp. 57-69, 2005.
- [13] A. Halder and R. Bhattacharya, 'Model validation: A probabilistic formulation,' in Decision and Control and European Control Conference (CDC-ECC), 2011 50th IEEE Conference on . IEEE, 2011, pp. 16921697.
- [14] --, 'Further results on probabilistic model validation in Wasserstein metric,' in Decision and Control (CDC), 2012 IEEE 51st Annual Conference on . IEEE, 2012, pp. 5542-5547.
- [15] --, 'Probabilistic model validation for uncertain nonlinear systems,' Automatica , vol. 50, no. 8, pp. 2038-2050, 2014.
- [16] A. Halder, K. Lee, and R. Bhattacharya, 'Optimal transport approach for probabilistic robustness analysis of F-16 controllers,' Journal of Guidance, Control, and Dynamics , vol. 38, no. 10, pp. 1935-1946, 2015.
- [17] R. E. Bellman, Dynamic Programming . Courier Dover Publications, 1957.
- [18] A. Halder and T. T. Georgiou, 'Gradient flows in uncertainty propagation and filtering of linear Gaussian systems,' 2017 IEEE Conference on Decision and Control, arXiv preprint arXiv:1704.00102 , 2017.
- [19] --, 'Gradient flows in filtering and Fisher-Rao geometry,' in 2018 Annual American Control Conference (ACC) . IEEE, 2018, pp. 42814286.
- [20] N. Parikh, S. Boyd et al. , 'Proximal algorithms,' Foundations and Trends R © in Optimization , vol. 1, no. 3, pp. 127-239, 2014.
- [21] R. Jordan, D. Kinderlehrer, and F. Otto, 'The variational formulation of the Fokker-Planck equation,' SIAM Journal on Mathematical Analysis , vol. 29, no. 1, pp. 1-17, 1998.
- [22] C. Villani, Topics in optimal transportation . American Mathematical Soc., 2003, no. 58.

| [23]   | L. Ambrosio, N. Gigli, and G. Savar´, e Gradient flows: in metric spaces and in the space of probability measures . Springer Science & Business Media, 2008.                                                                                           |
|--------|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [24]   | F. Santambrogio, ' { Euclidean, metric, and Wasserstein } gradient flows: an overview,' Bulletin of Mathematical Sciences , vol. 7, no. 1, pp. 87- 154, 2017.                                                                                          |
| [25]   | K. F. Caluya and A. Halder, 'Proximal recursion for solving the Fokker- Planck equation,' in 2019 Annual American Control Conference (ACC) , 2019.                                                                                                     |
| [26]   | J.-D. Benamou and Y. Brenier, 'A computational fluid mechanics so- lution to the monge-kantorovich mass transfer problem,' Numerische Mathematik , vol. 84, no. 3, pp. 375-393, 2000.                                                                  |
| [27]   | R. M. Colombo and M. L´cureux-Mercier, e 'Nonlocal crowd dynamics models for several populations,' Acta Mathematica Scientia , vol. 32, no. 1, pp. 177-196, 2012.                                                                                      |
| [28]   | J. Gomez-Serrano, C. Graham, and J.-Y. Le Boudec, 'The bounded confidence model of opinion dynamics,' Mathematical Models and Methods in Applied Sciences , vol. 22, no. 02, p. 1150007, 2012.                                                         |
| [29]   | Z. Askarzadeh, R. Fu, A. Halder, Y. Chen, and T. T. Georgiou, 'Stability theory of stochastic models in opinion dynamics,' IEEE Transactions on Automatic Control, doi: https://doi.org/10.1109/TAC.2019.2912490 , 2019.                               |
| [30]   | C. Villani, 'Trend to equilibrium for dissipative equations, functional inequalities and mass transportation,' Contemporary Mathematics , vol. 353, p. 95, 2004.                                                                                       |
| [31]   | J. A. Carrillo, R. J. McCann, and C. Villani, 'Kinetic equilibration rates for granular media and related equations: entropy dissipation and mass transportation estimates,' Revista Matem´tica a Iberoamericana , vol. 19, no. 3, pp. 971-1018, 2003. |
| [32]   | J. Karlsson and A. Ringh, 'Generalized Sinkhorn iterations for regular- izing inverse problems using optimal mass transport,' SIAM Journal on Imaging Sciences , vol. 10, no. 4, pp. 1935-1962, 2017.                                                  |
| [33]   | M. Cuturi, 'Sinkhorn distances: Lightspeed computation of optimal transport,' in Advances in neural information processing systems , 2013, pp. 2292-2300.                                                                                              |
| [34]   | J.-D. Benamou, G. Carlier, M. Cuturi, L. Nenna, and G. Peyr´, e 'Iterative Bregman projections for regularized transportation problems,' SIAM Journal on Scientific Computing , vol. 37, no. 2, pp. A1111-A1138, 2015.                                 |
| [35]   | Y. Chen, T. Georgiou, and M. Pavon, 'Entropic and displacement interpolation: a computational approach using the Hilbert metric,' SIAM Journal on Applied Mathematics , vol. 76, no. 6, pp. 2375-2396, 2016.                                           |
| [36]   | Y. Chen, T. T. Georgiou, and M. Pavon, 'On the relation between optimal transport and Schr¨dinger o bridges: A stochastic control viewpoint,' Journal of Optimization Theory and Applications , vol. 169, no. 2, pp. 671-691, 2016.                    |
| [37]   | P. E. Kloeden and E. Platen, Numerical solution of stochastic differential equations . Springer Science & Business Media, 2013, vol. 23.                                                                                                               |
| [38]   | M. Hutzenthaler, A. Jentzen, and P. E. Kloeden, 'Strong and weak diver-                                                                                                                                                                                |

- gence in finite time of Euler's method for stochastic differential equations with non-globally Lipschitz continuous coefficients,' Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences , vol. 467, no. 2130, pp. 1563-1576, 2010.
- [39] M. Hutzenthaler, A. Jentzen, P. E. Kloeden et al. , 'Strong convergence of an explicit numerical method for SDEs with nonglobally Lipschitz continuous coefficients,' The Annals of Applied Probability , vol. 22, no. 4, pp. 1611-1641, 2012.
- [40] M. Hutzenthaler and A. Jentzen, Numerical approximations of stochastic differential equations with non-globally Lipschitz continuous coefficients . American Mathematical Society, 2015, vol. 236, no. 1112.
- [41] D. J. Higham, X. Mao, and A. M. Stuart, 'Strong convergence of Euler-type methods for nonlinear stochastic differential equations,' SIAM Journal on Numerical Analysis , vol. 40, no. 3, pp. 1041-1063, 2002.
- [42] A. C. Thompson, 'On certain contraction mappings in a partially ordered vector space,' Proceedings of the American Mathematical Society , vol. 14, no. 3, pp. 438-443, 1963.
- [43] Y. Lim, 'Nonlinear equations based on jointly homogeneous mappings,' Linear Algebra and Its Applications , vol. 430, no. 1, pp. 279-285, 2009.
- [44] R. D. Nussbaum and D. Hilbert, Hilbert's projective metric and iterated nonlinear maps . American Mathematical Soc., 1988, vol. 391.
- [45] F. Malrieu, 'Uniform propagation of chaos for McKean-Vlasov equations,' https://tinyurl.com/ycxk5g9t, 2011, [Online; accessed 18November-2018].
- [46] J.-D. Benamou, G. Carlier, and M. Laborde, 'An augmented Lagrangian approach to Wasserstein gradient flows and applications,' ESAIM: Proceedings and Surveys , vol. 54, pp. 1-17, 2016.
- [47] M. Laborde, 'On some nonlinear evolution systems which are perturbations of Wasserstein gradient flows,' Topological Optimization and Optimal Transport: In the Applied Sciences , vol. 17, p. 304, 2017.
- [48] J. C. Cox, J. E. Ingersoll Jr, and S. A. Ross, 'A theory of the term structure of interest rates,' Econometrica , vol. 53, no. 2, pp. 385-408, 1985.
- [49] J. K. Møller and H. Madsen, From state dependent diffusion to constant diffusion in stochastic differential equations by the Lamperti transform . DTU Informatics, 2010.
- [50] H. Luschgy and G. Pag` es, 'Functional quantization of a class of Brownian diffusions: a constructive approach,' Stochastic Processes and their Applications , vol. 116, no. 2, pp. 310-336, 2006.
- [51] Y. Sun and M. Kumar, 'Uncertainty forecasting in the perturbed twobody problem via tensor decomposition,' in American Control Conference (ACC), 2016 . IEEE, 2016, pp. 5431-5436.
- [52] M. H. Duong, M. A. Peletier, and J. Zimmer, 'Conservative-dissipative approximation schemes for a generalized Kramers equation,' Mathematical Methods in the Applied Sciences , vol. 37, no. 16, pp. 2517-2540, 2014.
Kenneth F. Caluya Abhishek Halder
Gradient Flow Algorithms for Density Propagation in Stochastic Systems
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.