Gradient Flows in Uncertainty Propagation
and Filtering of Linear Gaussian Systems

arXiv:1704.00102v4 [math.OC] 29 Sep 2017

Abhishek Halder, and Tryphon T. Georgiou
University of California, Irvine
Abstract— The purpose of this work is mostly expository and
aims to elucidate the Jordan-Kinderlehrer-Otto (JKO) scheme
for uncertainty propagation, and a variant, the LaugesenMehta-Meyn-Raginsky (LMMR) scheme for filtering. We point
out that these variational schemes can be understood as
proximal operators in the space of density functions, realizing
gradient flows. These schemes hold the promise of leading to
efficient ways for solving the Fokker-Planck equation as well as
the equations of non-linear filtering. Our aim in this paper is
to develop in detail the underlying ideas in the setting of linear
stochastic systems with Gaussian noise and recover known
results.

I. I NTRODUCTION
n
Consider the gradient flow dx
dt = −∇ψ(x) in R , where
∇ is the gradient (w.r.t. the Euclidean metric) of a function
ψ(x), and consider the discretization

xk = xk−1 − h∇ψ(xk−1 ), for k ∈ N.

was to show that the minimizer of (2) approximates the
solution ρ(x, t) of the heat equation
∂ρ(x, t)
= ∆ρ(x, t), with ρ(x, 0) = ρ0 (x),
∂t
in the sense that ̺k (x, h) → ρ(x, t = kh), as h ↓ 0. Thus,
(2) establishes the remarkable result that the heat equation is
the gradient descent flow of the (negative) entropy integral
with respect to the Wasserstein metric.
An analogous JKO-like scheme was introduced recently
in Laugesen et al. [9] for the measurement update-step in
continuous-time filtering. More specifically, let us consider
the general system of Itô stochastic differential equations
(SDE’s)
p
(4a)
dx(t) = −∇U (x) dt + 2β −1 dw(t),
dz(t)

= c(x(t), t) dt + dv(t),

(4b)

(1)

where x ∈ Rn , z ∈ Rm , β > 0, U (·) is a potential, the
process and measurement noise processes w(t) and v(t) are
Wiener and satisfy E [dwi dwj ] = Qij dt ∀ i, j = 1, . . . , n
and E [dvi dvj ] = Rij dt ∀ i, j = 1, . . . , m, with Q, R ≻ 0,
respectively. Then x(t) and z(t) represent state and sensor
measurements at time t. Further, as usual, v(t) is assumed to
be independent of w(t) and independent of the initial state
x(0). Given the history of noise corrupted sensor data up to
time t, the filtering problem requires computing the posterior
probability distribution that obeys the Kushner-Stratonovich
stochastic PDE [19]–[21].

the solution, which depends on the choice of the step size
h, satisfies xk (h) → x(t = kh), as h → 0.

For the special case of trivial state dynamics, i.e., dx = 0,
and R the identity, Laugesen et al. [9] introduced

−
̺+
k (x, h) = arg inf {DKL ̺k̺k + hΦ(̺)}, k ∈ N, (5)

As is well known in finite-dimensional optimization,
1
xk = arg min{ kx − (xk−1 − h∇ψ(xk−1 ))k2 }
2
x
1
= arg min kx − xk−1 k2 + hψ(x) + o(h).
2
x
By recursively evaluating the proximal operator [1], [2]
xk

k·k

= proxhψ (xk−1 )
1
= arg min{ kx − xk−1 k2 + hψ(x)},
2
x

The Jordan-Kinderlehrer-Otto (JKO) scheme, introduced
in [3], is a similar recursion in the infinite-dimensional space
of density functions with respect to the Wasserstein geometry
[4], namely,
̺k (x, h) = arg min
̺

1 2
W (̺, ̺k−1 ) + hS(̺), k ∈ N, (2)
2 2

where W2 (·, ·) denotes the Wasserstein-2 distance between
two (probability) density functions,
Z
̺(x) log(̺(x))dx
(3)
S(̺) :=
Rn

is the negative differential entropy functional, and dx is the
volume element. In other words, (2) can be viewed as the
2
proximal operation proxW
hS (̺k−1 ). The main result in [3]

̺∈D2

with

1
E̺ {(yk − c(x))⊤ R−1 (yk − c(x))}, (6)
2
where yk is the noisy measurement in discrete-time defined
via yk := h1 ∆zk , ∆zk := zk − zk−1 , and {zk−1 }k∈N the
sequence of samples of z(t) at {tk−1 }k∈N for tk−1 := (k −
1)h. Laugesen et al. [9] proved that the LMMR equation (5)
approximates the solution of
h
⊤
dρ+ (x(t), t) = c(x(t), t) − Eρ+ {c(x(t), t)} R−1

dz(t) − Eρ+ {c(x(t), t)}dt ρ+ (x(t), t), (7)
Φ(̺) :=

i.e., of the Kushner-Stratonovich PDE corresponding to
+
dx = 0, in the sense that ̺+
k (x, h) ⇀ ρ (x(t), t) over

t ∈ [(k − 1)h, kh), as h ↓ 0. Thus, they showed that
in this special case, the Kushner-Stratonovich PDE is the
gradient descent of functional Φ(·) with respect to DKL , i.e.,
−
KL
computed by proxD
hΦ (̺k ).
The purpose of the present paper is to develop this circle of
ideas, namely, that both uncertainty propagation and filtering
can be viewed as gradient flows in the special case of linear
stochastic systems with Gaussian noise. In fact, we consider
the general case of the linear stochastic system
dx(t) = Ax(t) dt + B dw(t),

(8)

where w(t) is a Wiener process as before, though possibly
not of the same dimension as x. We suppose that the uncertain initial condition x(0) has a known Gaussian PDF, the
matrix A is Hurwitz, and that the diffusion matrix B is such
that (A, B) is a controllable pair. For this, we recover the
well-known propagation equations (see for example [6, Ch.
3.6]) for the mean and covariance of the state x(t) out of the
JKO-scheme via a two-step optimization. The applicability
of the JKO-scheme to (8) is not immediately obvious since
the development in [3] requires the state dynamics to be
in the canonical form (4a) with the drift being a gradient
and the diffusion coefficient being a positive scalar. We
further show that this two-step optimization procedure that
we introduce, can be used to derive the Kalman-Bucy filter
from a generalized version of the LMMR equation (5). We
remark that variational schemes for estimator/observer design
based on gradient flows can also be seen as regularized
dynamic inversion in the spirit of [23].
Notation
Throughout we use bold-faced upper-case letters for matrices, and bold-faced lower case letters for vectors. The
notation I stands for identity matrix of appropriate dimension, we use tr(·) and det(·) to respectively denote the trace
and determinant of a matrix, and the symbols ∇ and △
denote the gradient and Laplacian operators, respectively. We
denote the space of probability
R density functions (PDFs) on
Rn by
D
:=
{ρ
:
ρ
≥
0,
Rn ρ = 1}, by D2 := {ρ ∈
R
D | Rn x⊤ x ρ(x)dx < ∞} the space of PDFs with finite
second moments, by Dµ,P denote the space of PDFs which
share the
R same mean vector µ and same covariance matrix
P := Rn (x − µ)(x − µ)⊤ ρ(x)dx. Likewise, let Dµ,τ
denote the space of PDFs which have the same mean µ and
same trace of covariance τ := tr(P ) > 0. Clearly, Dµ,P ⊂
Dµ,τ ⊂ D2 ⊂ D. We use the symbol N (µ, P ) to denote
a multivariate Gaussian PDF with mean µ, and covariance
P . The notation x ∼ ρ means that the random vector x has
PDF ρ; and E {·} denotes the expectation operator while,
when
the probability density is to be specified, Eρ {·} :=
R
(·)ρ(x)dx.
n
R
II. JKO S CHEME

IN

G ENERAL

We now discuss in some detail the JKO scheme for the
case of the diffusion process in (4a), and the corresponding

Fokker-Planck equation [7]
∂ρ
= ∇ · (∇U (x)ρ) + β −1 △ρ, ρ(x, 0) = ρ0 (x). (9)
∂t
To this end we first introduce the Wasserstein metric, the
free energy, and the Kullback-Leibler divergence.
The Wasserstein-2 distance W2 (ρ1 , ρ2 ) between a pair
of PDFs ρ1 (x), ρ2 (y) ∈ D (or, even between probability
measures, in general), supported on X , Y ⊆ Rn , is

 12
Z
2
(10)
W2 (ρ1 , ρ2 ) :=
inf
kx − yk2 dσ(x, y)
dσ∈Π(ρ1 ,ρ2 )

X ×Y

where Π (ρ1 , ρ2 ) is a probability measure on the product
space X × Y having finite second moments and marginals
ρ1 , ρ2 , respectively. It is well known that W2 : D × D 7→
[0, ∞) is a metric [4, p. 208]. Further, its square W22 (ρ1 , ρ2 )
represents the smallest amount of “work” needed to “morph”
ρ1 into ρ2 [5]. The infimum is achieved over a space of
measures, and under mild assumptions, the minimizing dσ
has support on the graph of the optimal “transportation map”
T : X 7→ Y that pushes ρ1 to ρ2 . Alternatively, one may
view the optimization problem in (10) as seeking the joint
distribution of two random vectors x and y, distributed
according to
 ρ1 and ρ2 respectively, that minimizes the
variance E kx − yk22 .

Another important notion of distance that enters into our
discussion, which however is not a metric, is the KullbackLeibler divergence (also known as relative entropy) between PDFs or positive Rmeasures in general. This is given
dρ1
dρ1
1
by DKL (dρ1 kdρ2 ) := ( dρ
dρ2 ) log( dρ2 )dρ2 where dρ2 denotes the Radon-Nikodym derivative. When1 dρi = ρi (x)dx,
i ∈ {1, 2}, are absolutely continuous with respect to the
Lebesgue measure, then
Z
ρ1 (x)
ρ1 (x) log
DKL (dρ1 kdρ2 ) =
dx.
ρ
2 (x)
Rn

Gradient flow requires
an energy functional, which we
R
denote by E (ρ) := U (x)ρ(x)dx, where U (·) is the potential energy. Then, a stochastically driven gradient flow is
modeled by the Itô SDE (4a) and the Fokker-Planck equation
(9) for the corresponding PDF as before. The stationary
solution of (9) is the Gibbs distribution ρ∞R(x) = Z1 e−βU(x) ,
where the normalization constant Z := Rn e−βU(x) dx is
known as the partition function. The distance to equilibrium
which, in a way, quantifies the amount of work that the
system can deliver, is captured by the so-called free energy
functional F (ρ), defined as the sum of the energy functional
E (ρ) and the negative differential entropy S (ρ) given in (3),
that is,
F (ρ) := E (ρ) + β −1 S (ρ)


= β −1 DKL ρke−βU(x) .

(11a)
(11b)

1 Here we use a slight abuse of notation in that we denote both, the
measure and the density with the same symbol.

For the case of (4a), the JKO scheme becomes

parameterized subset of D2 in such a way that when we
optimize
the functionals 21 W 2 (̺, ̺0 ) and hF (̺) individually
1
̺k (x, h) = arg inf { W22 (̺, ̺k−1 ) + h F (̺)}, k ∈ N, (12) over this chosen subspace, the arginfs (which are achieved)
2
̺∈D2
of the two individual optimization problems match. Hence,
for step-size h > 0, and initialized by a given ̺0 (satisfying the sum of the two has the same arginf over the chosen
F (̺0 ) < ∞). For U (x) ≡ 0, (12) reduces to (2). Solving subspace. In the second step, we optimize over the subspace
(12) results in a sequence of PDFs {̺k (x, h)}k∈N in D2 . parameters. Our choice for the parameterized set of densities
It can be shown following [3] that ̺k (x, h) ⇀ ρ (x(t), t) is Dµ,P ⊂ D2 , i.e., the PDFs with given mean-covariance
weakly in L1 (Rn ) for t ∈ [(k − 1)h, kh), k ∈ N, as h ↓ 0.
pair (µ, P ); the choice of the optimal pair is to be decided
in the second optimization step.
III. JKO S CHEME FOR LINEAR G AUSSIAN SYSTEMS
The development below requires several technical lemmas
that
are collected in the Appendix.
We now develop and solve the JKO scheme for the linear
Gaussian system in (8) with ρ0 = N (µ0 , P0 ) and Q ≡ I,
without loss of generality. Therefore, we are concerned with
the linear Fokker-Planck (Kolmogorov’s forward) PDE
n

∂ρ
∂2
1 X
= −∇ · (ρAx) +
ρBQB ⊤ ij . (13)
∂t
2 i,j=1 ∂xi ∂xj

Under the stated assumptions, it is well-known that (13)
admits a steady-state, which is Gaussian with mean zero
and covariance P∞ ≻ 0 that uniquely solves the algebraic Lyapunov equation AP∞ + P∞ A⊤ + BQB ⊤ = 0.
Also, starting from ρ0 (x) = N (µ0 , P0 ), the transient is
ρ(x(t), t) = N (µ(t), P (t)) where the µ(t) and P (t) satisfy
the following ordinary differential equations (ODEs) [6, Ch.
3.6]
µ̇(t) = Aµ(t), µ(0) = µ0 ,

(14a)

Ṗ (t) = AP (t) + P (t)A⊤ + BQB ⊤ , P (0) = P0 . (14b)
Below, we recover these equations using the JKO scheme.
First, in Section III-A, wep
explain how this is done when
A is symmetric and B ≡ 2β −1 I, β > 0, in which case,
Ax = −∇U (x) for a suitable potential. The general case,
in Section III-B, is more involved and requires to view the
drift as the gradient of a time-varying potential.
p
A. The case where A is symmetric and B ≡ 2β −1 I
p
2β −1 I, the pair (A, B) is controllable.
Since B ≡
Further, since A is Hurwitz and symmetric, Γ := −A ≻ 0,
and utilizing the potential
1 ⊤
x Γx ≥ 0,
2
we can cast (8) in the canonical form (4a). Then,

1 ⊤
E(̺) := E[U (x)] =
µ Γµ + tr (ΓP ) ,
2
where P is the covariance of x. Notice that E(·) depends
on the PDF of x only via its mean and covariance.
U (x) :=

To carry out the optimization (12) over D2 , we adopt
a two-step strategy. Our approach is motivated by the
observation that the objective function in (12) is a sum
of two functionals. In the first step, we choose a suitable

1) Optimizing over Dµ,P : Given ̺0 ≡ ρ0 = N (µ0 , P0 ),
and a µ and P ≻ 0, we first determine
1
̺1 = arg inf { W22 (̺, N (µ0 , P0 )) + h F (̺)}.
(15)
̺∈Dµ,P 2
1
From Lemma 2 we see that arg inf W22 (̺, N (µ0 , P0 ))
2
̺∈Dµ,P
is achieved by ̺ = N (µ, P ) (uniquely). From Lemma 3,
since U (x) = 21 x⊤ Γx, we also know that arg inf h F (̺) is
̺∈Dµ,P

achieved by ̺ = N (µ, P ) (uniquely). Thus, ̺1 = N (µ, P ).
The infimal value in (15) is now the sum of the two infima,
"
 1 !#

2
1
1
h
1
2
2
2
+
k µ − µ0 k2 + tr P + P0 − 2 P0 P P0
2
2β
h
i
−n − n log(2π) − log det(P ) + βµ⊤ Γµ + β tr (ΓP ) . (16)

2) Optimizing over (µ, P ): Equating the gradient of (16)
w.r.t. µ to zero, results µ = φ(µ0 ) := (I + hΓ)−1 µ0 . The
recursion µk = φ(µk−1 ), up to first order in h, becomes
µk = (I − hΓ) µk−1 + O(h2 ).

(17)

We see that this recursion coincides with the solution of
(14a) in the “small h” limit. Specifically, µ(t) = eAt µ0 ⇒
k
µk := µ(t = kh) = eAh µ0 ⇒ µk = eAh µk−1 =
(I + hA) µk−1 + O(h2 ), which is same as (17) since Γ :=
−A. Thus, we have recovered (14a) using discrete timestepping via JKO scheme in the small step-size limit.
Setting the gradient of (16) w.r.t. P to zero (using Lemma
4), we obtain

1 1
1
h
−1
−1 2
P02 − P −1 + hΓ = 0. (18)
I − P02 P0 2 P −1 P0 2
β
− 12

By pre and post multiplying both sides of (18) with P0
1

−1 2
−1
=: Z, we arrive at
and letting P0 2 P −1 P0 2

,

β
β −1
−1
Z − P0 2 (I + hΓ)P0 2 = 0,
h
h
which admits the unique closed-form solution [16, p. 304]

 21 !
h − 12
β
− 12
−I + I + 4 P0 (I + hΓ) P0
. (19)
Z=
2h
β
Z2 +

Expanding (19), we obtain




β
1 h
−I + I + 4 P0 (I + hΓ) P0 +
2h
2 β

#
1 1
2
1
−
1
16h
−1
−
−1
3
2 2
2
2
(I + hΓ) P0 (I + hΓ) P0 + O(h )
P
2!
β2 0


h
−1
−1
= P0 2 I + hΓ − P0−1 P0 2 + O(h2 ).
(20)
β

Z=

−1
2

−1
2

1

−1 2
−1
back into (20), squarSubstituting Z = P0 2 P −1 P0 2
ing, and rearranging, we get that



−1 
−1

1
1
P = I + h Γ − P0−1
P0 I + h Γ − P0−1
+O(h2 )
β
β






1
1
P0 I − h Γ − P0−1
+ O(h2 )
= I − h Γ − P0−1
β
β
= Ψ(P0 ) + O(h2 ),


where Ψ(P0 ) := P0 + h −ΓP0 − P0 Γ + 2β −1 I . Set the
matrix-valued recursion Pk = Ψ (Pk−1 ), where


Ψ (Pk−1 ) := Pk−1 + h −ΓP0 − P0 Γ + 2β −1 I + O(h2 ). (21)

To show that (21) indeed recovers (14b),
p first notice that
substituting A = A⊤ = −Γ and B = 2β −1 I in (14b)
results the Lyapunov differential equation
Ṗ (t) = −ΓP (t) − P (t)Γ + 2β −1 I

subject to P (0) = P0 , which can be solved via the method
of integrating factor as

1
(22)
P (t) = Γ−1 I − e−2Γt + e−Γt P0 e−Γt .
β
Thus, for t = kh, (22) gives

Pk := P (kh) = β −1 Γ−1 I − e−2Γkh + e−Γkh P0 e−Γkh

1) Equipartition of energy coordinate transformation:
Consider
the stationary covariance P∞ associated with
√
(A, 2B) that satisfies
AP∞ + P∞ A⊤ + 2BB ⊤ = 0.

For a system at a stationary distribution, we define the
thermodynamic temperature θ as the average amount of
“energy” per degree of freedom, that is,
1
θ := tr(P∞ ),
n
and, thereby, β := θ−1 the inverse temperature. By pre and
−1
post multiplying (23) with P∞ 2 , and rescaling by θ so as to
preserve the temperature in the new coordinates, we get
√
√
Aep θI + θIA⊤
2θBep ( 2θBep )⊤ = 0,
(24)
ep +

Pk − Pk−1 = 2β −1 hI − hΓP0 − hP0 Γ + O(h2 ),
which is same as (21) derived from JKO scheme. Thus,
we have recovered the covariance evolution through FokkerPlanck dynamics using the time-stepping procedure via JKO
scheme in the small step-size limit.

−1

1

−1

2
, Bep := P∞ 2 B, while
where Aep := P∞ 2 AP∞
the stationary covariance θI reflects equipartition of energy.√The equipartition√of energy co-ordinate transformation
(A, 2B) 7→ (Aep , 2θBep ), corresponds to the state√ −1
transformation x 7→ xep := θP∞ 2 x, leading to
√
dxep (t) = Aep xep (t) dt + 2θBep dw(t).
(25)

This settles how β is to be defined and interpreted in the
context of JKO scheme (11) and (12). On the other hand,
Aep being similar to A, is guaranteed to be Hurwitz but not
symmetric, unless A was symmetric to begin with. Thus, it
remains for us to “symmetrize” Aep and define a suitable
potential energy U (·) as needed in (12). We do this next.
2) Symmetrization transformation: We introduce the
time-varying transformation

= 2β −1 khI + (P0 − khΓP0 − khP0 Γ) + O(h2 ).

Replacing k with k−1 in the latter yields a similar expression
for Pk−1 . Then, subtracting these expressions for Pk−1 from
Pk we obtain that

skew

where

Askew
ep

xep 7→ xsym := e−Aep
:=

1
2 (Aep

skew

t

t

xep

A⊤
ep ).

−
This results in
√
√
(Aep , 2θBep ) 7→ (F (t), 2θG(t)),

with
F (t) := e−Aep

skew

Aep
Asym
ep e

where, similarly, Asym
:=
ep
satisfies

t

skew

, and G(t) := e−Aep

1
2 (Aep

dxsym (t) = F (t)xsym (t) dt +
B. The case of Hurwitz A and controllable (A, B)
√
We scale B into 2B without loss of generality, and
take as initial PDF ̺0 ≡ ρ0 = N (µ0 , P0 ). Since we allow
any Hurwitz√(not necessarily symmetric) A, and any B that
makes (A, 2B) a controllable pair, it is not apparent if
and how one can express (8) in the canonical form (4a).
The main impediment in doing so, is twofold: (1) how to
define the potential energy U (x), and (2) how to interpret and
define the parameter β in the generic case. In the following,
we show that by two successive co-ordinate transformations,
system (8) can indeed be put in the form (4a). Similar
transformations have been mentioned in [17, p. 1464], [18]
in a different context.

(23)

t

Bep ,

+ A⊤
ep ). Thus, xsym (t)
√
2θG(t) dw(t).

(26)

Notice that F (t) is symmetric for all t. Furthermore, observe that the new coordinates xsym is simply obtained by a
(time-varying) orthogonal transformation of the equipartition
of energy coordinates xep . Hence the stationary covariance
of xsym is identical to that of xep , which is θI (from Section
III-B.1). What happens is that the covariance of xsym (t)
tends to the same steady state value as t → ∞ in spite of
the fact that (26) has time varying coefficients. To see this
⊤
in different way, we can rewrite (24) as Bep Bep
= −Asym
ep ,
and deduce that
skew

G(t)G(t)⊤ = e−Aep

skew

⊤ Aep t
Bep Bep
e
= −F (t),
√
√
⇒ F (t)θI + θIF (t) + 2θG(t)( 2θG(t))⊤ = 0.
t

(27)

The symmetrization xep 7→ xsym leaves the stationary
covariance θI invariant. This guarantees that the definition
of temperature θ stays intact. The coordinate transformations
described above are summarized in Table I.
3) Recovery of the Fokker-Planck solution: We are now
ready to apply the JKO scheme to√ the generic stochastic
linear system dx(t) = Ax(t) dt + 2B dw(t), with initial
PDF ρ(x(0), 0) = N (µ0 , P0 ). To this end, we carry out a
computation akin to the two steps in Section III-A, for the
transformed SDE (26) in the symmetrized coordinate xsym .
From there on, we recover the Fokker-Planck solution in the
original coordinate x.
Since x 7→ xsym is a linear transformation, it follows
that xsym ∼ N (µsym , Psym ) whenever x ∼ N (µ, P ).
Thus, carrying out the first step of the optimization in xsym
coordinate, we get an expression similar to (16) wherein
(µ, P ) is to be replaced by (µsym , Psym ), and (µ0 , P0 ) is
to be replaced by (µsym0 , Psym0 ). To carry out the second
⊤
step of optimization, notice that Asym
= −Bep Bep
 0,
ep
skew
skew
−Aep t sym Aep t
and consequently F (t) = e
Aep e
 0. Thus,
considering the time-varying potential
1
U (xsym ) := − x⊤
F (t)xsym ≥ 0,
2 sym
and setting the partial derivative of the infimal value from
first stage of the optimization w.r.t. µsym to zero, results
the recursion µsymk = (I − hF (kh))−1 µsymk−1 . Recalling
skew √
−1
that xsym = e−Aep t θP∞ 2 x, we arrive at a recursion in
original coordinate:
1

skew

2 Aep
µk = P∞
e

kh

skew

{(I − hF (kh))−1 eAep
e

−Askew
kh
ep

−1
P∞ 2

h

}

µk−1 .

(28)

By series expansion and collecting linear terms in h, one can
verify the following:
2
(I − hF (kh))−1 = I + hAsym
ep + O(h ),
−1 Askew
h
ep

(I − hF (kh))
e

Askew
kh
ep

e

where C ∈ Rm×n , and ρ0 = N (µ0 , P0 ). The conditional
PDF ρ+ (x(t), t) = N (µ+ (t), P + (t)), given measurements
up to time t, is well-known and given by the Kalman-Bucy
filter [22]

dµ+ (t) = Aµ+ (t)dt + K(t) dz(t) − Cµ+ (t)dt , (31a)
Ṗ +(t) = AP +(t)+P +(t)A⊤+2BB ⊤−K(t)RK(t)⊤ (31b)

that specifies a vector SDE and a matrix ODE, respectively,
for the conditional mean µ+ (t) and covariance P + (t). The
initial conditions are µ+ (0) = µ0 , P + (0) = P0 , and
K(t) := P + (t)C ⊤ R−1 is the so-called Kalman gain.
In the sequel, we demonstrate that by applying the twostep optimization strategy we used before in Section III, we
can recover the Kalman-Bucy filter from LMMR-equation
(5) for the linear Gaussian case as the h ↓ 0 limit.
A. LMMR gradient descent scheme
Once again we proceed with carrying out the following
two optimization steps. First, we optimize (5) over Dµ,P ,
and then optimize the minimum value over the choice of
parameters (µ, P ).
−
−
1) Optimizing over Dµ,P : Consider ̺−
k = N (µk , Pk )
to be our prior for the state PDF at time t = kh. Observe
that

inf

̺∈Dµ,P

− 
DKL ̺kN (µ−
k , Pk ) =

inf

̺∈Dµ,P

−E̺ {log

e

F ILTERING

dz(t) = Cx(t) dt + dv(t),

= I + hAep + O(h ),
e

FOR

In this section, we focus on the linear Gaussian filtering
problem, with process model and measurement models
√
dx(t) = Ax(t)dt + 2Bdw(t),

2

−1 Askew
h −Askew
kh
ep
ep

(I − hF (kh))

IV. JKO- LIKE S CHEMES

Z

̺(x) log ̺(x)dx

Rn

−
N (µ−
k , Pk )} ,

(32)

and that

2

= I + hAep + O(h ).

−1
1h
⊤
(µ − µ−
Pk−
(µ − µ−
k)
k )
2
 1

+ tr P (Pk− )−1 − log (2π)n det(Pk− )
2

−
E̺ {log N (µ−
k , Pk )} = −

Hence (28) yields


1
−1
2
Aep P∞ 2 µk−1 + O(h2 )
µk = I + hP∞
= (I + hA) µk−1 + O(h2 ),

− 21

(29)
1

2
where the last equality follows from Aep := P∞ AP∞
.
hA
2
Since µ̇ = Aµ and e = I + hA + O(h ), in the small h
limit, equation (29) thus recovers (14a), as in Section III-A.
A similar straightforward but tedious computation leads to
the matrix recursion

Pk − Pk−1 = h(APk−1 + Pk−1 A⊤ + 2BB ⊤ ) + O(h2 ),
(30)
which in the limit h ↓ 0, is indeed a first-order approximation
of the Lyapunov equation for the original system. We omit
the details for brevity.

remains invariant for all ̺ ∈ Dµ,P . Therefore, the arginf
in (32) is achieved by the Gaussian PDF N (µ, P ) (i.e., the
maximum entropy PDF with given mean-covariance), and
−
the infimal value is precisely DKL (N (µ, P )kN (µ−
k , Pk )).
On the other hand, notice that
1h
1
(yk − Cµ)⊤
E̺ {(yk − Cx)⊤ R−1 (yk − Cx)} =
2
2

i
R−1 (yk − Cµ) + tr C ⊤ R−1 CP
= constant
(33)
inf

̺∈Dµ,P

❤❤❤
→
❤❤❤Coordinate
❤❤❤
Attribute ↓
❤❤

Original

Equipartition of energy

Symmetrization

State vector

x

xep

xsym

(A,

System matrices
Stationary covariance

√

2B)

P∞

(Aep ,

√

2θBep )

θI

(F (t),

√

2θG(t))

θI

TABLE I: Summary of the coordinate transformations for Section III-B.

as well over Dµ,P . Hence


−
arg inf DKL ̺kN (µ−
k , Pk )

B. Alternative JKO-like schemes for filtering

̺∈Dµ,P

+



h
E̺ {(yk − Cx)⊤ R−1 (yk − Cx)} = N (µ, P ),
2

and the corresponding infimum value is

1
− −1
⊤
tr (Pk− )−1 P + (µ−
(µ−
k − µ) (Pk )
k − µ) − n−
2

h
(yk − Cµ)⊤ R−1 (yk − Cµ)
log det (Pk− )−1 P +
2

+ tr C ⊤ R−1 CP .
(34)

2) Optimizing over (µ, P ): Equating the partial derivative
of (34) w.r.t. µ to zero, and setting µ ≡ µ+
k in the resulting
algebraic equation, we get


+
⊤ −1
(Pk− )−1 µ−
yk − Cµ+
k − µk + hC R
k = 0,

−
− ⊤ −1
(35)
yk − Cµ+
⇒µ+
k .
k = µk + hPk C R
On the other hand, equating the partial derivative of (34)
w.r.t. P to zero, and then setting P ≡ Pk+ in the resulting
algebraic equation, we get


(Pk+ )−1 = (Pk− )−1 + hC ⊤ R−1 C ⇒ Pk+ = I + hPk− C ⊤
−1 −
R−1 C
Pk = Pk− − hPk− C ⊤ R−1 CPk− + O(h2 ).
(36)

With ∆zk = yk h, as in Section I,

2
µ+ (t)dt = µ+
k h + O(h ),

and from (29) that
+
2
µ−
k = (I + hA)µk−1 + O(h ).

These, together with (36), allow us to simplify (35) as

2
+
+
+ ⊤ −1
µ+
∆zk − hCµ+
k −µk−1 = hAµk−1 + Pk C R
k + O(h ),

which in the limit h ↓ 0, leads to (31a).

Substituting (30) into (36) we arrive at
+
+
−hPk−1
C ⊤ R−1 CPk−1
+ O(h2 ).

In the limit h ↓ 0, (37) recovers (31b).

̺+
k (x, h) = arg inf
̺∈D2


1 2
W2 ̺, ̺−
+ hΦ(̺),
k
2

k ∈ N, (38)

where the functional Φ(·) is as in (6). The template of the
two-step optimization again applies and, specializing to the
linear Gaussian case, the solution of (38) in the h ↓ 0 limit,
is N (µ+ (t), P + (t)), given by

dµ+ (t) = Aµ+ (t)dt + L dz(t) − Cµ+ (t)dt ,
(39a)

Ṗ + (t)=(A − LC)P + (t)+P + (t)(A − LC)⊤+2BB ⊤(39b)
where L := C ⊤ R−1 , and µ+ (0) = µ0 , P + (0) = P0 . This
follows by noticing from Sections III-A.1 and IV-A.1 that



1 2
−
−
W ̺, N (µk , Pk + hΦ(̺) = N (µ, P ),
arg inf
2 2
̺∈Dµ,P

where the infimum value is


 1 

1
− 21 2
−
− 12
2
P
(P
)
k µ−µ−
k
+
tr
P
+
P
−2
(P
)
k
k 2
k
k
2

h
+ (yk − Cµ)⊤ R−1 (yk − Cµ)+tr C ⊤ R−1 CP .(40)
2

dz(t) = ∆zk + O(h2 ),

+
+
+
Pk+ − Pk−1
= h(APk−1
+ Pk−1
A⊤ + 2BB ⊤ )

The ideas in the LMMR-scheme suggest the possibility
of alternative variational schemes to approximate stochastic
estimators. Such a viewpoint has been put forth in [23],
promoting the notion of regularized dynamic inversion. As an
example, one may consider a gradient descent with respect
to the Wasserstein distance 12 W22 , instead of KL-divergence
DKL in (5). In that case, the posterior may be constructed
according to

(37)

Equating the partial derivative of (40) w.r.t. µ to zero, then
+
+
setting µ ≡ µ+
k , and using (29), we find (µk −µk−1 ) equals

⊤ −1
2
∆zk − hCµ+
hAµ+
k + O(h ),
k−1 + C R

(41)

which in the limit h ↓ 0, results the SDE (39a). Similarly,
using Lemma 4, we equate the partial derivative of (40) w.r.t.
P to zero, and then setting P ≡ Pk+ , we get


(Pk+ )−1 = I + hC ⊤ R−1 C (Pk− )−1 I + hC ⊤ R−1 C ⇒

Pk+ = Pk− − h Pk− C ⊤ R−1 C + C ⊤ R−1 CPk− + O(h2 ),

+
+
which combined with the recursion Pk− = Pk−1
+h(APk−1
+
⊤
⊤
2
+ Pk−1 A + 2BB ) + O(h ) from Section III-B.3, yields

 +
+
+
Pk+ = Pk−1
+h A − C ⊤ R−1 C Pk−1
+ Pk−1
(A−
i

⊤
2
⊤ −1
⊤
C R C + 2BB + O(h ). (42)

In the limit h ↓ 0, recursion (42) gives Lyapunov ODE (39b).
It is instructive to compare the SDE-ODE system (39)
with that in (31). In the case of (39), the estimator is of
a Luenberger type with a static gain matrix L which is
decoupled from the covariance, unlike (31). The estimator
(39) is obviously not optimal in the minimum mean-square
error sense. It is only presented here as a guideline to explore
other variational schemes with desirable properties.
V. C ONCLUDING

REMARKS

Reformulating uncertainty propagation and the filtering
equations as gradient flows [8] is potentially transformative
[3] [9]. The full power of this viewpoint is yet to be
uncovered. Moreover, casting the iterative approximation
steps in the language of proximal operators on the space
of density functions may provide theoretical insights and
computational benefits. A specific direction of future work
would be developing proximal algorithms [2] to numerically
solve the nonlinear filtering problem by recursively solving
convex optimization problems, and to quantify computational
performance of the same with respect to existing sequential
Monte Carlo algorithms like the particle filter. The purpose
of the present paper has been to highlight and elucidate
the ideas in [3] and [9] in the context of linear Gaussian
systems. We hope that this study will help to motivate further
exploration of this topic.
A PPENDIX
In this Appendix, we collect some lemmas that are used in
Sections III and IV. In addition, we will show in Corollary
1 below that applying Lemma 1 and 2 together enables us
to provide an alternative proof of a Theorem in [14], which
might be of independent interest.
Lemma 1: If X and Y are symmetric positive definite
 1
 21
p
1
matrices, then tr X 2 Y X 2
≤ tr (X) tr (Y ).

Proof: From Uhlmann’s variational formula (see [10],
also Theorem 6.1 in [11]) , given any G ≻ 0, we have

 21  p
1
1
≤ tr (XG) tr (Y G−1 ),
(43)
tr
X2Y X2

where the equality in (43) is achieved for the specific choice
 1
− 21 1 1
1
1
1
Gopt = Y 2 X 2 Y X 2
X 2 Y 2 X − 2 . Specializing
(43) for G = Y , and noting that tr (XY ) ≤ tr (X) tr (Y ),
the statement follows.

Lemma 2: Given a PDF ̺0 (x) ∈ D2 with mean
µ0 ∈ Rn , and n × n covariance matrix P0 ≻ 0. Then
inf W22 (̺, ̺0 ) equals

̺∈Dµ,P

kµ−

µ0 k22


 1
 21 
1
2
2
,
+ tr P + P0 − 2 P0 P P0

(44)

and is achieved by push-forward of ̺0 (x) via an affine trans 1
− 12 1
1
1
port map M x + m, where M := P 2 P 2 P0 P 2
P 2,
and m := µ− µ0 , that is, the arginf for (44) is ̺(x)
 =
q
 1
 12
1
1
1
det(P0 )
−2
P 2 P0 P 2
P − 2 (x − µ) + µ0 . In
det(P ) ̺0 P
particular, if ̺0 = N (µ0 , P0 ), then ̺ = N (µ, P ).

Proof: Let ̺0 be as given, and choose any ̺ ∈ Dµ,P .
Let ̺0 and ̺ be obtained by translating ̺0 and ̺ respectively,
such that both ̺0 and ̺ have zero mean. Using (10), we can
directly verify [12, p. 236] that W22 (̺, ̺0 ) =k µ − µ0 k22
+W22 (̺, ̺0 ). On the other hand, it is known [13, p. 11,
Proposition 1.1.6] that

 1
 21 
1
2
2
2
W2 (̺, ̺0 ) ≥ tr P + P0 − 2 P0 P P0
⇒ W22 (̺, ̺0 ) ≥ right hand side of (44).

(45)

Now consider a candidate transport map M x + m where
M and m are functions of P , P0 , µ, µ0 as in the statement.
It suffices to prove that our candidate transport map indeed
achieves the equality in (45). To this end, directly substituting
the expressions for M and m, notice that the push-forward
has mean M µ0 + m = µ, and covariance M P0 M ⊤ =
P . So our candidate transport map (M , m) is feasible.
To show optimality, from (10) notice that W22 (̺, ̺0 ) =
inf tr(P + P0 − 2C), where C := M P0 solves P0 −
C∈Rd×d

CP −1 C ⊤  0, which has known optimal solution Copt :=
 1
− 21 1
1
1
Mopt P0 = P0 P 2 P 2 P0 P 2
P 2 . Since our candi1


−
2
1
1
1
1
P 2 satisfies tr (M P0 ) =
date M := P 2 P 2 P0 P 2


 12
1
1
, the statement follows.
tr (Mopt P0 ) = tr
P02 P P02
In the Corollary below, combining Lemma 1 and 2, we
recover a result in [14, Theorem 3.1].

Corollary 1: Given d-dimensional joint PDF ̺0 with
mean µ0 , covariance P0 ≻ 0, suppose tr(P0 ) = τ0 . For
fixed µ and τ > 0,
√ 2
√
inf W22 (̺, ̺0 ) =
τ − τ 0 + k µ − µ0 k22 , (46)
̺∈Dµ,τ

and is achieved by ̺(x) =

τ0
τ

 d2

̺0

τ0
τ


(x − µ) + µ0 .

Proof: Let us choose P := ττ0 P0 , and from (44)
√ 2
√
observe that inf W22 (ξ, ̺0 ) = ( τ − τ 0 ) + k µ −
ξ∈Dµ,P

µ0 k22 . On the other hand, for any ̺ ∈ Dµ,τ , we know from

(45) that


1

 12
1

W22 (̺, ̺0 ) ≥ τ + τ0 − 2 tr P02 SP02

+ k µ − µ0 k22 ,

where S is the covariance of ̺. Using Lemma 1, we get
 21
 1
1
√
≤ τ τ0 ⇒ W22 (̺, ̺0 ) ≥
tr P02 SP02
√
√ 2
inf W22 (ξ, ̺0 ) =
τ − τ 0 + k µ − µ0 k22 , (47)

[7] H. Risken, The Fokker-Planck Equation: Methods of Solution and Applications. Springer Series in Synergetics, Vol. 18, First ed., Springer;
1989.
[8] L. Ambrosio, N. Gigli, and G. Savaré, Gradient Flows: in Metric
Spaces and in the Space of Probability Measures. Lectures in Mathematics, ETH Zürich, Second ed., Birkhäuser; 2008.
[9] R.S. Laugesen, P.G. Mehta, S.P. Meyn, and M. Raginsky, “Poisson’s
Equation in Nonlinear Filtering”. SIAM Journal on Control and
Optimization. Vol. 53, No. 1, pp. 501–525, 2015.

ξ∈Dµ,P

and that the equality is achieved when S = P = ττ0 P0 . In
 d
that case, det(S) = ττ0 det(P0 ), and hence Lemma 2
yields the arg inf ̺(x) for (46) as
s


 1
 12
1
det(P0 )
− 12
− 12
2
2
̺0 S
S P0 S
S (x − µ) + µ0
det(S)

 τ  d2  τ
0
0
̺0
(x − µ) + µ0 .
=
τ
τ

Lemma 3: If E(·) depends on ̺ only via the mean and
covariance of ̺, then inf F (̺) is achieved by N (µ, P ).
̺∈Dµ,P

Proof: As E(̺) ≡ E(µ, P ), henceR from (11) we get
̺ log ̺ dx. Since
inf F (̺) = E(µ, P ) + β −1 inf

̺∈Dµ,P

̺∈Dµ,P

[10] A. Uhlmann, “The “Transition Probability” in the State Space of a
*-Algebra”. Reports on Mathematical Physics, Vol. 9, No. 2, pp. 273–
279, 1976.
[11] D. Petz, Quantum Information Theory and Quantum Statistics. Theoretical and Mathematical Physics, First ed., Springer; 2008.
[12] C.R. Givens, and R.M. Shortt, “A Class of Wasserstein Metrics for
Probability Distributions”. The Michigan Mathematical Journal, Vol.
31, No. 2, pp. 231–240, 1984.
[13] S.T. Rachev, and L. Rüschendorf, Mass Transportation Problems.
Volume I: Theory. First ed., Springer; 1998.
[14] E.A. Carlen, and W. Gangbo, “Constrained Steepest Descent in the
2-Wasserstein Metric”. Annals of Mathematics, pp. 807–846, 2003.
[15] A. Halder, and E.D.B. Wendel, “Finite Horizon Linear Quadratic
Gaussian Density Regulator with Wasserstein Terminal Cost”. Proceedings of the 2016 American Control Conference, pp. 7249–7254,
2016.

N (µ, P ) is the maximum entropy PDF under prescribed
mean µ and covariance P , hence the statement.

[16] N.J. Higham, and H.M. Kim, “Solving A Quadratic Matrix Equation
by Newton’s Method with Exact Line Searches”. SIAM Journal on
Matrix Analysis and Applications, Vol. 23, No. 2, pp. 303–316, 2001.

Lemma 4: For P , P0 ≻ 0,
 21
1 1
 1
1
∂
1 1  −1
−1 2
tr P02 P P02
= P02 P0 2 P −1 P0 2
P02 .
∂P
2
Proof: We refer the readers to Appendix B in [15].

[17] D. Liberzon, and R.W. Brockett, “Spectral Analysis of Fokker–Planck
and Related Operators Arising from Linear Stochastic Differential
Equations”. SIAM Journal on Control and Optimization, Vol. 38, No.
5, pp. 1453–1467, 2000.

R EFERENCES
[1] H.H. Bauschke, and P.L. Combettes, Convex Analysis and Monotone
Operator Theory in Hilbert Spaces. CMS Books in Mathematics,
Springer; 2011.

[18] R.W. Brockett, and J.C. Willems, “Stochastic Control and the Second
Law of Thermodynamics”. Proceedings of the 1978 IEEE Conference
on Decision and Control including the 17th Symposium on Adaptive
Processes, pp. 1007–1011, 1978.
[19] R.L. Stratonovich, “Application of the Theory of Markov Processes for
Optimum Filtration of Signals”. Radio Eng. Electron. Phys. (USSR),
Vol. 1, pp. 1–19, 1960.

[2] N. Parikh, and S. Boyd, “Proximal Algorithms”. Foundations and
Trends in Optimization. Vol. 1, No. 3, pp. 127–239, 2014.

[20] H.J. Kushner, “On the Differential Equations Satisfied by Conditional
Densities of Markov Processes, with Applications”. Journal of the
SIAM Series A Control, Vol 2, No. 1, pp. 106–119, 1964.

[3] 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.

[21] M. Fujisaki, G. Kallianpur, and H. Kunita, “Stochastic Differential
Equations for the Non Linear Filtering Problem”. Osaka Journal of
Mathematics, Vol. 9, No. 1, pp. 19–40, 1972.

[4] C. Villani, Topics in Optimal Transportation. Graduate Studies in
Mathematics, Vol. 58, First ed., American Mathematical Society; 2003.

[22] R.E. Kalman, and R.S. Bucy, “New Results in Linear Filtering and
Prediction Theory”. Journal of Basic Engineering, Vol. 83, No. 3, pp.
95–108, 1961.

[5] J.-D. Benamou, and Y. Brenier, “A Computational Fluid Mechanics
Solution to the Monge–Kantorovich Mass Transfer Problem”. Numerische Mathematik. Vol. 84, No. 3, pp. 375–393, 2000.
[6] K.J. Åström, Introduction to Stochastic Control Theory. Academic
Press; 1970.

[23] A. Yezzi, and E.I. Verriest, “Nonlinear Observers via Regularized
Dynamic Inversion”. Proceedings of the 2007 American Control
Conference, pp. 1693–1698, 2007.

