Gradient Flows in Filtering and Fisher-Rao Geometry

arXiv:1710.00064v2 [math.OC] 29 Oct 2017

Abhishek Halder1 , and Tryphon T. Georgiou2
Abstract— Uncertainty propagation and filtering can be interpreted as gradient flows with respect to suitable metrics in the
infinite dimensional manifold of probability density functions.
Such a viewpoint has been put forth in recent literature,
and a systematic way to formulate and solve the same for
linear Gaussian systems has appeared in our previous work
where the gradient flows were realized via proximal operators
with respect to Wasserstein metric arising in optimal mass
transport. In this paper, we derive the evolution equations as
proximal operators with respect to Fisher-Rao metric arising in
information geometry. We develop the linear Gaussian case in
detail and show that a template two step optimization procedure
proposed earlier by the authors still applies. Our objective is
to provide new geometric interpretations of known equations
in filtering, and to clarify the implication of different choices
of metric.

I. I NTRODUCTION
This paper concerns with an emerging viewpoint that gradient flows on the infinite-dimensional manifold of probability density functions can be constructed to approximate both
the propagation and measurement update steps in filtering.
In the systems-control literature, continuous-time filtering
theory has traditionally been approached from a “transport
viewpoint” where given the noisy process and measurement
models, flow of the posterior probability density function
(PDF) of the state vector is described by the KushnerStratonovich stochastic partial differential equation (PDE)
[1], [2] (or equivalently by the Duncan-Mortensen-Zakai
PDE [3]–[5] for the unnormalized posterior). In the absence
of measurement update, the problem reduces to that of
uncertainty propagation subject to the drift and diffusion in
the process dynamics, and the evolution of the state PDF
is then governed by the Fokker-Planck or Kolmogorov’s
forward PDE [6]. Due to numerical difficulties to solve these
PDEs, in practice one resorts to approximate the evolution
of the state PDFs using sequential Monte Carlo algorithms,
which remain computationally challenging in general. This
motivates us to look for alternative formulations that are
theoretically equivalent to solving the associated transport
PDEs.
One such alternative is the “variational viewpoint” that has
been put forth in recent literature [7], [8]. Specifically, let us
consider the system of Itô stochastic differential equations
1 Department of Applied Mathematics & Statistics, University of California, Santa Cruz, CA 95064, USA, ahalder@ucsc.edu
2 Department of Mechanical and Aerospace Engineering, University of
California, Irvine, CA 92617, USA, tryphon@uci.edu

(SDEs) for the state (x(t)) and measurement (z(t)) vectors
p
(1a)
dx(t) = −∇U (x) dt + 2β −1 dw(t),
dz(t)

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

(1b)

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. Also,
v(t) is assumed to be independent of w(t) and independent
of the initial state x(0). In the absence of (1b), the FokkerPlanck PDE for the state PDF ρ(x, t) corresponding to (1a)
is given by
∂ρ
= ∇ · (∇U (x)ρ) + β −1 △ρ, ρ(x, 0) = ρ0 (x). (2)
∂t
In the absence of (1a), the Kushner-Stratonovich PDE [1],
[2] for the state PDF ρ+ (x, t), conditioned on the history of
measurements till time t, corresponding to (1b) is given by
h
⊤
dρ+ (x(t), t) = c(x(t), t) − Eρ+ {c(x(t), t)} R−1
i
dz(t) − Eρ+ {c(x(t), t)}dt ρ+ (x(t), t). (3)

In the “variational viewpoint” for uncertainty propagation,
one sets up a discrete time-stepping scheme of the form
̺k (x, h) = arg inf
̺

1 2
d (̺, ̺k−1 ) + hΦ(̺), k ∈ N,
2

(4)

with step-size h, a distance functional d(·, ·) between two
PDFs, and a functional Φ(·) that depends on the drift and
diffusion coefficients in the process model. The functionals
d(·, ·) and Φ(·) are to be chosen such that ̺k (x, h) →
ρ(x, t = kh) as h ↓ 0, thus establishing consistency
between the solutions of (2) and (4). Similar recursion for
the measurement update takes the form
−
1 2
̺+
k (x, h) = arg inf 2 d (̺, ̺k ) + hΦ(̺), k ∈ N,

(5)

̺

+
where ̺−
k and ̺k approximates the prior and posterior PDFs,
respectively, and the functional Φ(·) depends on the measurement model and noisy observations. Again, the choices of
+
d(·, ·) and Φ(·) must guarantee ̺+
k (x, h) → ρ (x, t = kh)
as h ↓ 0, to allow consistency between the solutions of (3)
and (5).

Both (4) and (5) can be viewed as evaluating suitable
(infinite dimensional) proximal operators [9], [10] proxdhΦ (·)
of functional hΦ(·) with respect to the distance functional
d(·, ·). In other words, recursion (4) can be succinctly written
−
d
as ̺k = proxdhΦ (̺k−1 ). Similarly, ̺+
k = proxhΦ (̺k ) for (5).
In particular, (4) and (5) can be seen as the gradient descent
of the functional Φ(·) with respect to the distance d(·, ·). For

a parallel with finite dimensional gradient descent, see [11,
Section I].
For the proximal recursion (4) associated with uncertainty
propagation, it was shown in [7] that if d2 (·, ·) is chosen to
be the squared Wasserstein-2 distance, which equals the cost
of optimal transport between two probability measures under
consideration, then Φ(·) is the free energy functional, defined
to be the sum of an energy functional and scaled negative
of the differential entropy functional (see [11, Section II] for
details). In this case, the functional Φ(·) depends on both
U (·) and β appearing in the process model (1a).
Laugesen, Mehta, Meyn and Raginsky [8] were first to
consider the proximal recursion (5) associated with measurement update, where it was shown that if 12 d2 (·, ·) is
chosen to be the Kullback-Leibler divergence, then Φ(·) is
the “expected quadratic surprise” functional, given by
Φ(̺) :=

1
E̺ {(yk − c(x))⊤ R−1 (yk − c(x))},
2

(6)

where yk := h1 ∆zk , ∆zk := zk − zk−1 , and {zk−1 }k∈N is
the sequence of samples of z(t) at {tk−1 }k∈N for tk−1 :=
(k − 1)h. We refer this particular variational formulation as
the LMMR scheme.
In [11], a general framework was introduced to bring any
stable controllable linear system to the canonical form (1)
for which the proximal recursions (4) or (5) apply. Also,
a two step optimization strategy was proposed to compute
these recursions, and using these tools, the LMMR scheme
in the linear Gaussian setting was shown to recover the
Kalman-Bucy filter [12]. It was also found [11, Section
IV.B] that the recursion (5) with d2 (·, ·) chosen to be the
squared Wasserstein-2 distance, and Φ(·) as in (6) in the
linear Gaussian case, results a Luenberger-type estimator
with static gain matrix, unlike the Kalman-Bucy filter.
It becomes apparent that the choice of the distance functional d(·, ·) cannot be arbitrary, and to appeal the gradient
descent interpretation, should preferably define a metric on
the manifold of PDFs. While Wasserstein-2 distance is a
metric [13, p. 208], the Kullback-Leibler divergence is not.
This naturally leads to the question whether one can find a
metric d(·, ·) and functional Φ(·) in (5) such that the filtering
equations can be seen as the gradient descent of Φ(·) with
respect to metric d(·, ·). In this paper, by taking d(·, ·) to be
the geodesic distance induced by the Fisher-Rao metric [15]
and Φ(·) as in (6), we answer this in the affirmative for the
linear Gaussian setting. Specifically, a variant of the two-step
optimization template proposed in [11] allows us to perform
explicit computation for (5) with Fisher-Rao metric, and is
shown to recover the Kalman-Bucy filter in the limit h ↓ 0.
It is perhaps not surprising that the functional Φ(·) in this
paper and in the LMMR scheme [8] are the same, the choice
of d(·, ·) are different, and yet both of these two schemes
are shown to recover the same filtering equations. This is
due to the well-known fact [14, Ch. 2, p. 26–28] that the

Kullback-Leibler divergence between a pair of PDFs where
one PDF is infinitesimally perturbed from the other, equals
to half the squared geodesic distance in Fisher-Rao metric
measured between the same.
The rest of this paper is structured as follows. In Section
II, we review some basic aspects of Fisher-Rao geometry
on the manifold of probability density functions. In Section
III, focusing on the linear Gaussian case, we approximate
stochastic estimator as gradient descent with respect to the
geodesic distance induced by the Fisher-Rao metric. Section
IV concludes the paper.

Notations and Preliminaries
As in [11],R we denote the space of PDFs on Rn by D :=
{ρ : ρ ≥ 0, Rn ρ = 1}, and the space Rof PDFs with finite
second moments by D2 := {ρ ∈ D | Rn x⊤ x ρ(x)dx <
∞}. The notation Dµ,P is used for the space of PDFs
which share the
R same mean vector µ and same covariance
matrix P := Rn (x − µ)(x − µ)⊤ ρ(x)dx. The following
inclusion is immediate: Dµ,P ⊂ 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,R when the probability density is
to be specified, Eρ {·} := Rn (·)ρ(x)dx.

We use d(·) for the differential, D(·) for the Jacobian, and
vec(·) for the vectorization operator. For taking derivatives of
matrix valued functions, we utilize the Jacobian identification
rules [16, p. 199, Table 2]. Notation In stands for the n × n
identity matrix, and k · kFro for the Frobenius norm. The
set of real n × n matrices is denoted as Mn , the set of
symmetric matrices as Sn ⊂ Mn , and the set of positive
+
definite matrices as S+
n ⊂ Sn . The cone Sn is a smooth
differentiable manifold with tangent space TX S+
n = Sn ,
where X ∈ S+
.
The
symbols
⊗
and
⊕
denote
Kronecker
n
product and Kronecker sum1 , respectively. We will have
multiple occasions to use the following three well-known
facts: (1) vec(·) is a linear operator, (2) (A ⊗ B)(C ⊗
D) = (AC ⊗ BD), and (3) an identity relating vec(·) and
Kronecker product:

vec(ABC) = C ⊤ ⊗ A vec(B).

(7)

All matrix logarithms in this paper are to be understood
as principal logarithms. Furthermore, the following property
(see [17], [18, Theorem 1.13(c)]) will be useful

A log(B)A−1 = log ABA−1 ,

(8)

which also holds when log is replaced by any analytic matrix
function.
1 For an m × m matrix A and an n × n matrix B, the Kronecker sum
is mn × mn matrix A ⊕ B := A ⊗ In + Im ⊗ B.

II. F ISHER -R AO G EOMETRY

and a parametric version of the geodesic distance (10) is

Given PDF ρ ∈ D and tangent vectors u, v ∈ Tρ D,
the Fisher-Rao metric is defined as a Riemannian metric
Gρ (u, v) on D, given by
Z
uv
d x.
(9)
Gρ (u, v) :=
Rn ρ
√
Let ψ :=
ρ, and consider the space of
R “square-root
PDFs”, given by Ψ := {ψ : ψ ≥ 0, Rn ψ 2 = 1}.
Geometrically, Ψ represents the non-negative orthant of the
(infinite dimensional) unit sphere imbedded in the Hilbert
space of square-integrable functions
on Rn equipped with
R
√
the inner product hψ1 , ψ2 i := Rn ψ1 ψ2 . For ψi := ρi ,
i = 1, 2, the (minimal) geodesic distance induced by FisherRao metric, denoted hereafter as dFR (ρ1 , ρ2 ), is then simply
the minor arc-length along the great circle connecting ψ1 and
ψ2 , i.e.,
√ √
(10)
dFR (ρ1 , ρ2 ) = arccosh ρ1 , ρ2 i.
Since any point on the chord joining ψ1 and ψ2 can be
represented as (1 − t)ψ1 + tψ2 , t ∈ [0, 1], hence the
minimal geodesic ρ(x, t) connecting ρ1 and ρ2 , is obtained
by parameterizing the arc ψ(x, t) connecting ψ1 and ψ2 , i.e.,
ψ(x, t) =

t)2

(1 − t)ψ1 + tψ2

1/2

((1 −
+ t2 + 2t(1 − t)hψ1 , ψ2 i)

2
p
p
(1 − t) ρ1 (x) + t ρ2 (x)
⇒ρ (x, t) =
√ √ .
(1 − t)2 + t2 + 2t(1 − t)h ρ1 , ρ2 i

(11)

The formula (9), (10) and (11) are non-parametric in the
sense that they do not require the PDFs to have any finite
dimensional parametric co-ordinate representations. If prior
knowledge allows one to consider a known parametric family
of PDFs ρ(x|θ) with r-dimensional parameter vector θ =
(θ1 , . . . , θr )⊤ , then the Fisher-Rao metric (9) admits local
coordinate representation, given by the Fisher information
matrix [19]


∂
∂
ρ(x|θ),
ρ(x|θ)
(12)
gij (θ) = Gρ(x|θ)
∂θi
∂θj


∂
∂
=E
log ρ(x|θ)
log ρ(x|θ)
(13)
∂θi
∂θj

 2
∂
log ρ(x|θ) , i, j = 1, . . . , r. (14)
= −E
∂θi θj
The parametric version of the geodesic (11) becomes
ρ (x, t) ≡ ρ (x|θ(t)), t ∈ [0, 1], where θ(t) solves the EulerLagrange boundary value problem
X
θ̈k (t)+ Γkij (θ(t))θ̇i (t)θ̇j (t) = 0, i, j, k = 1, . . . , r, (15)
i,j

with θ1 := θ(0) and θ2 := θ(1) given, and Γkij are
Christoffel symbols of the second kind defined via gij . The
squared arc-length
X
gij (θ(t)) d θi (t) d θj (t),
(16)
(d s)2 =
i,j

dFR (ρ(x|θ1 ), ρ(x|θ2 )) =
=

Z

t=1

d s(t)

Zt=0
1sX
0

gij (θ(t))θ̇i (t)θ̇j (t) d t. (17)

i,j

Of particular interest to us, is the family of multivariate
Gaussian PDFs. In this case, for µ, µ1 , µ2 ∈ Rn , and
P , P1 , P2 ∈ S+
n , (15) can be written as [21, Theorem 6.1]
µ̈ − Ṗ P −1 µ̇ = 0,

⊤

P̈ + µ̇µ̇ − Ṗ P

−1

(18a)

Ṗ = 0,

(18b)

with t ∈ [0, 1], (µ1 , P1 ) := (µ(0), P (0)) and (µ2 , P2 ) :=
(µ(1), P (1)). Furthermore, (16) becomes
1
trace((P d P )2 ).
(19)
2
The following special cases are well-known [20]–[22]:
(d s)2 = d µ⊤ P −1 d µ +


1/2
,
dFR (N(µ1 , P ) , N(µ2 , P )) = (µ1 − µ2 )⊤P −1 (µ1 − µ2 )
(20)

which is also the Mahalanobis distance [23], and that
dFR (P1 , P2 ) := dFR (N(µ, P1 ) , N(µ, P2 ))
(21)


1
−1/2
−1/2
= √ log P1
P2 P1
2
Fro



2 1/2
1
−1/2
−1/2
trace log P1
P2 P1
(22)
=
2
!1/2
n
1X
(log λi )2
,
=
2 i=1

where λ1 , . . . , λn are the eigenvalues of P1−1 P2 . To the best
of our knowledge, no closed-form expression is known for
dFR (N (µ1 , P1 ), N (µ2 , P2 )).
III. P ROXIMAL R ECURSION WITH RESPECT
F ISHER -R AO M ETRIC

TO

In this Section, we consider approximating continuoustime stochastic estimator via the recursive variational scheme
dFR
−
̺+
k (x, h) = proxhΦ (̺k ) = arg inf
̺∈D2

1 2
d (̺, ̺−
k ) + hΦ(̺), (23)
2 FR

+
where k ∈ N, h is the step-size, and ̺−
k and ̺k are the
approximations of the prior and posterior PDFs, respectively.
We choose the functional Φ(·) as in (6), and focus on the
linear Gaussian case with process and measurement models
√
dx(t) = Ax(t)dt + 2Bdw(t),
(24a)
dz(t) = Cx(t) dt + dv(t),
(24b)

where C ∈ Rm×n and ρ0 = N (µ0 , P0 ). It is wellknown that ρ+ (x(t), t) = N (µ+ (t), P + (t)) and the optimal
estimator is given by the Kalman-Bucy filter [12], consisting
of the following SDE-ODE system:

dµ+ (t) = Aµ+ (t)dt + K(t) dz(t) − Cµ+ (t)dt ,
+

+

+

⊤

⊤

⊤

Ṗ (t) = AP (t)+P (t)A +2BB −K(t)RK(t) ,

(25a)
(25b)

where K(t) := P + (t)C ⊤ R−1 is the Kalman gain.
We next demonstrate that evaluating the proximal recursion (23) for the linear Gaussian case in the h ↓ 0 limit,
indeed recovers the Kalman-Bucy filter (25). To carry out the
optimization over D2 in (23), we adopt a two step strategy
proposed in [11]. In the first step, we optimize the objective
in (23) over the parameterized subspace Dµ,P ⊂ D2 . In the
second step, we optimize over the subspace parameters µ
and P . The main insight behind this strategy comes from the
fact that the objective in (23) is a sum of two functionals.
If we can find a parameterized subspace of D2 , over which
the individual arg infs of these two functionals match, then
it must also be the arg inf of their sum. We detail these
steps below. The ensuing calculations require some technical
results collected in the Appendix.
A. Two Step Optimization
−
−
1) Optimizing over Dµ,P : Let ̺−
k = N (µk , Pk ) be the
prior state PDF at time t = kh. From Lemma 1 in the
Appendix, we know

−
(26)
arg inf d2FR ̺, N (µ−
k , Pk ) = N (µ, P ).
̺∈Dµ,P

On the other hand, for any ̺ ∈ Dµ,P , we have

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

i
R−1 (yk − Cµ) + trace C ⊤ R−1 CP
= constant. (27)

Consequently
arg inf
̺∈Dµ,P



h
1 2
− 
dFR ̺, N (µ−
E̺ {(yk − Cx)⊤ R−1
k , Pk ) +
2
2

(yk − Cx)} = N (µ, P ), (28)

and the corresponding infimum value is
1 2
hh
−
(yk − Cµ)⊤ R−1
dFR (N (µ, P ), N (µ−
k , Pk )) +
2
2

i
(yk − Cµ) + trace C ⊤ R−1 CP ,

which is not convenient for our next step due to
the lack of availability of closed-form expression for
−
d2FR (N (µ, P ), N (µ−
k , Pk )). This can be circumvented as
follows. Instead of choosing Dµ,P with two free parameters, we choose two different parameterized subspaces, viz.
Dµ,P − and Dµ− ,P , each with one free parameter. We next
k
k
carry out the two step optimization strategy for each of the
two one-parameter subspaces.
Similar to (28), the arg inf for Dµ,P − is N (µ, Pk− ), and
k
using (20), the corresponding infimum is
hh
1
⊤ −1
−
(yk − Cµ)⊤ R−1 (yk − Cµ)
(µ − µ−
k ) Pk (µ − µk ) +
2
2 
i
+ trace C ⊤ R−1 CPk− .
(29)

Likewise, the arg inf for Dµ− ,P is N (µ−
k , P ), and using
k
the notation (21), the corresponding infimum is
1 2
hh
⊤ −1
(yk − Cµ−
(yk − Cµ−
dFR (P , Pk− ) +
k ) R
k)
2
2

i
+ trace C ⊤ R−1 CP . (30)

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


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

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

−1  h ⊤ −1
−1
1
+ C R C = 0,
log Pk+ Pk−
Pk+
2
2


−1
+
−
+ ⊤ −1
⇒Pk Pk
= exp −hPk C R C

= In − hPk+ C ⊤ R−1 C + O(h2 ). (32)
−1
Pre-multiplying both sides of (32) by Pk+
results
−1

−1
In + hPk− C ⊤ R−1 C + O(h2 ),
= Pk−
Pk+
−1 −
⇒ Pk+ = In + hPk− C ⊤ R−1 C
Pk + O(h2 )
= Pk− − hPk− C ⊤ R−1 CPk− + O(h2 ).

(33)

Interestingly, equations (31) and (33) are same as those
obtained (see equations (35) and (36) in [11]) by optimizing
(5) over Dµ,P with same Φ(·) as in (6), but 12 d2 (·, ·) chosen
to be the Kullback-Leibler divergence.
B. Recovering the Kalman-Bucy Filter
We recall two basic relations (equations (29) and (30) in
[11]) for the propagation step obtained by discrete-time stepping of (4) with Wasserstein-2 distance as d(·, ·), and freeenergy functional as Φ(·), for the linear Gaussian process
model (24a):
µk = (In + hA) µk−1 + O(h2 ),

(34)
⊤

⊤

Pk = Pk−1 + h APk−1 +Pk−1 A +2BB



2

+O(h ). (35)

We refer the readers to [11, Section III.B] for details of their
derivations. Intuitively, (34) and (35) can be viewed as the
Euler discretizations of the well-known mean and covariance
(Lyapunov) ODEs for uncertainty propagation.
Using ∆zk = yk h (from Section I), d z(t) = ∆zk +
2
O(h2 ), µ+ (t) d t = µ+
k h+O(h ), (34) and (35), we simplify
(31) as

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

which in the limit h ↓ 0, recovers (25a).

On the other hand, combining (33) with (35) results
+
+
+
Pk+ − Pk−1
= h(APk−1
+ Pk−1
A⊤ + 2BB ⊤ )
+
+
−hPk−1
C ⊤ R−1 CPk−1
+ O(h2 ),

which in the limit h ↓ 0, recovers (25b).
Thus we have demonstrated that the measurement update
step in Kalman-Bucy filter (25) can be interpreted as the
gradient descent of functional (6) with respect to the metric
dFR for c(x) = Cx. This complements our earlier result [11,
Section III] showing the propagation step can be interpreted
as the gradient descent of certain free energy functional with
respect to the Wasserstein-2 metric.

expression, using (7), and noticing that d vec(F (X)) =
D F (X) d vec(X), we arrive at D F (X) = X ⊤ ⊕ X as
claimed.
(ii) The directional derivative of G(X) in the direction
Z ∈ TX S+
n = Sn , denoted as d GZ (X), is given by (see
[27, Section 3.1], also [28, Appendix I.B])
log (X + hZ) − log X
d GZ (X) := lim
h→0
h
Z ∞
−1
=
(X + τ In ) Z(X + τ In )−1 d τ. (36)
0

Recall that
vec (d GZ (X)) = D G(X) vec(Z).

IV. C ONCLUSIONS
This paper contributes to an emerging research program
that views the filtering equations as gradient flux or steepest
descent on the manifold of probability density functions. For
gradient descent to be meaningful in infinite dimensions, one
requires a metric with respect to which distance between
probability density functions are to be measured. By using
the geodesic distance induced by the Fisher-Rao metric, we
have shown equivalence between a discrete time-stepping
procedure for gradient descent on the manifold of conditional
probability density functions and the filtering equations in
the linear Gaussian setting. Here, our intent has been to
understand the underlying geometric implications, and to
work out the theoretical details to recover known facts. Our
hope is that the results in this paper will help motivate
developing these ideas in nonlinear filtering setting and to
solve them via proximal algorithms [10].
A PPENDIX
In this Appendix, we collect several technical results that
are used in Section III.
Lemma 1: Given µ, µ0 ∈ Rn and P , P0 ∈ S+
n,

arg inf d2FR (ρ, N (µ0 , P0 )) = N (µ, P ) .

ρ∈Dµ,P

Proof: This is a consequence of Cramér-Rao inequality;
see Theorem 20 in [24, p. 1512], also Section 1 in [25].
Lemma 2 and 3 below provide preparatory steps for
proving Theorem 1 that follows.
Lemma 2: The Jacobians of the matrix valued
F (X) := X 2 for X ∈ Mn , G(X) := log X for
and H(X) := AXB for real A, X, B of sizes
the product is defined, are respectively given by
⊤
(i) D F (X) = X
R ∞ ⊕ X,
(ii) D G(X) = 0 ((X + τ In ) ⊗ (X + τ In ))−1
(iii) D H(X) = B ⊤ ⊗ A.

functions
X ∈ S+
n,
such that

(37)

Applying vec(·) to both sides of (36) yields
Z ∞

vec (d GZ (X)) = vec (X + τ In )−1 Z(X + τ In )−1 d τ
0
Z ∞

(7)
(X + τ In )−⊤ ⊗ (X + τ In )−1 vec(Z) d τ
=
Z0 ∞

−1
=
((X + τ In ) ⊗ (X + τ In )) d τ vec(Z), (38)
0

where the last step is due to the fact that inverse of Kronecker
product equals Kronecker product of inverses (in the same
order), and that X is symmetric. Equating (37) and (38), the
statement follows.
(iii) See [26, Appendix A, Lemma 3].
Z ∞
Lemma 3: For X ∈ S+
,
(X + τ In )−1 (log X) (X +
n
0

τ In )−1 d τ = X −1 log X = (log X)X −1 .

Proof: Consider the spectral decomposition X =
V ΛV −1 , where Λ := diag(λi ) is positive diagonal matrix.
Then log X = V (log Λ) V −1 . Writing In = V V −1 , the
integral of interest equals
Z ∞



log λi
log λi
−1
V diag
dτ V
= V diag
V −1 .
2
λi
0 (λi + τ )
Hence the result.
Theorem 1: For P , P0 ∈ S+
n,


∂ 2
dFR (P , P0 ) = P −1 log P P0−1 = log P0−1 P P −1 .
∂P

Proof: Let F (X) := X 2 , G(X) := log X, H(X) :=
and Ω := 21 F ◦G◦H. From (22), it is clear
∂
that our objective is to compute Y := ∂P
trace (Ω(P )).

−1/2
−1/2
P0
XP0
,

Notice that
trace (d Ω(P ))

d τ,

Proof: (i) Taking the differential operator d(·) to both
sides of F (X) = XX results d F (X) = (d X)X +
X(d X). Applying vec(·) to both sides of this resulting

= (vec(In ))⊤vec (d Ω(P ))
= (vec(In ))⊤D Ω(P ) d vec(P )
1
= (vec(In ))⊤D F (G (H(P ))) D G (H(P )) D H(P )
2
d vec(P ) (39)

where the last step is the chain rule for Jacobians. Since


−1/2
−1/2
−1/2
−1/2
P , P0 , P0
P P0
∈ S+
,
log
P
P
P
∈ Sn ,
n
0
0
{z
}
|
|
{z
}
=H(P )

=G(H(P ))

therefore using Lemma 2, we have
D H(P ) =

−1/2
P0

−1/2
P0

⊗
Z ∞

∈

−1/2

S+
n,

(40)

−1/2

+ τ In ) ⊗
−1
−1/2
−1/2
d τ ∈ S+
(P0
P P0
+ τ In )
n,

D G (H(P )) =

0

(P0

P P0

∂
∂P

trace (Ω(P )), we also


d trace(Ω(P )) = vec Y ⊤ d vec(P ).

(43)

Comparing the right hand sides of (39) and (43), substituting
(40), (41), (42), and using the shorthand H for H(P ) =
−1/2
−1/2
P0
P P0
, we obtain
 Z ∞

1
−1/2
−1/2
P0
⊗ P0
((H + τ In ) ⊗
vec(Y ) =
2
0

−1
(H + τ In )) d τ (log H ⊕ log H) vec(In )
Z ∞


1
1/2
1/2
((H + τ In ) ⊗
vec(Y ) =
⇒ P0 ⊗ P0
2
0

−1
(H + τ In )) d τ (log H ⊕ log H) vec(In ). (44)


1/2
1/2
Using (7), we rewrite the LHS of (44) as vec P0 Y P0
.
On the other hand, observe that (log H ⊕ log H) vec(In ) =
(log H ⊗ In + In ⊗ log H) vec(In ) = vec (2 log H) (using (7) again). Thus (44) simplifies to


1/2
1/2
vec P0 Y P0
Z ∞

=
(H + τ In )−1 ⊗ (H + τ In )−1 d τ vec(log H)
0

=
(2)

Z

∞

Z0

(45)

(H + τ In )−1 ⊗ (H + τ In )−1 vec(log H) d τ

Z

∞

(H + τ In )−1 ⊗ (H + τ In )−1 d τ
Z ∞

= (V ⊗V )
(Λ + τ In )−1 ⊗(Λ + τ In )−1 d τ (V ⊗V )−1 ,
{z
}
| 0
0

2

2

the n × n diagonal matrix L has entries Ljj =
Zwhere
∞
dτ
1
=
for j = (i − 1)n + i, i = 1, . . . , n;
2
(λ
+
τ
)
λ
i
0
Z ∞ i
log λi − log λj
dτ
and Ljj =
for j =
=
(λ
+
τ
)(λ
+
τ
)
λi − λj
i
j
0
[(i − 1)n + iZ+ 1, i(n + 1)], i = 1, . . . , n − 1. Noting that
1
a−b
, we can rewrite
for a, b > 0,
aτ b1−τ d τ =
log a − log b
0
(45) as
Z 1



1/2
1/2
−1
τ
1−τ
(V ⊗V )
Λ ⊗Λ
d τ (V ⊗V ) vec P0 Y P0
{z
}
| 0
=L−1
Z 1



1/2
1/2
=
= vec(log H),
H τ ⊗H 1−τ d τ vec P0 Y P0
0
Z 1
(7)
1/2
1/2
⇒
H τ P0 Y P0 H 1−τ d τ = log H,
(48)
0

providing an integral characterization of (47), i.e., the gradient in Theorem 1 satisfies the matrix equation (48). Indeed,
the solution of (48) is given by [30, equations (5.4.3) and
(5.4.8)]
1/2

P0

1/2

Y P0

=

Z

∞

(H + τ In )−1 (log H) (H + τ In )−1 d τ
0

(Lemma 3)

=

H −1 log H = (log H)H −1 ,

thus establishing the equivalence between (47) and (48).
R EFERENCES
[1] 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.

∞


vec (H + τ In )−1 (log H) (H + τ In )−1 d τ
0 Z

∞
= vec
(H + τ In )−1 (log H) (H + τ In )−1 d τ
=

Remark 1: Starting from (45), one can alternatively derive
(47) via spectral decomposition of H, which reveals an
appealing connection with matrix logarithmic mean [29].
Specifically, letting H = V ΛV −1 , we get

=:L

(41)


−1/2
−1/2
D F (G (H(P ))) = log P0
P P0
⊕


−1/2
−1/2
∈ Sn . (42)
log P0
P P0
By virtue of the definition Y :=
have

Using (8) on (47), the proof is complete.

[2] 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] T.E. Duncan, “Likelihood Functions for Stochastic Signals in White
Noise”. Information and Control, Vol. 16, No. 4, pp. 303–310, 1970.

0

(Lemma 3)

=



vec H −1 log H = vec (log H) H

which yields the matrix equation
1/2


−1

,

[4] R.E. Mortensen, “Stochastic Optimal Control with Noisy Observations”. International Journal of Control, Vol. 4, No. 5, pp. 455–464,
1966.

1/2

= H −1 log H = (log H) H −1


−1/2
−1/2
−1/2
1/2
,
P0
P P0
⇒Y = P −1 P0 log P0


1/2 −1
−1/2
−1/2
−1/2
P0 P .
= P0
log P0
P P0
P0

(46)

Y P0

[5] M. Zakai, “On the Optimal Filtering of Diffusion Processes”. Probability Theory and Related Fields, Vol. 11, No. 3, pp. 230–243, 1969.

(47)

[6] H. Risken, The Fokker-Planck Equation: Methods of Solution and Applications. Springer Series in Synergetics, Vol. 18, First ed., Springer;
1989.

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

[20] C. Atkinson, and A.F.S. Mitchell, “Rao’s Distance Measure”. Sankhya:
The Indian Journal of Statistics, Series A (1961–2002), Vol. 43, No.
3, pp. 345–365, 1981.

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

[21] L.T. Skovgaard, “A Riemannian Geometry of the Multivariate Normal
Model”. Technical Report, No. 167, Department of Statistics, Stanford
University, April 1981.

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

[22] C.R. Rao, “Differential Metrics in Probability Spaces”. Differential
Geometry in Statistical Inference, Institute of Mathematical Statistics
Lecture Notes and Monographs Series, Vol. 10, pp. 217–240, 1987.

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

[23] P.C. Mahalanobis, “On the Generalized Distance in Statistics”. Proceedings of the National Institute of Sciences (Calcutta), Vol. 2, No.
1, pp. 49–55, 1936.

[11] A. Halder, and T.T. Georgiou, “Gradient Flows in Uncertainty
Propagation and Filtering of Linear Gaussian Systems”. 56th IEEE
Conference on Decision and Control, 2017. arXiv:1704.00102
[math.OC], available online: https://arxiv.org/pdf/1704.
00102.pdf.
[12] 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.
[13] C. Villani, Topics in Optimal Transportation. Graduate Studies in
Mathematics, Vol. 58, First ed., American Mathematical Society; 2003.
[14] S. Kullback, Information Theory and Statistics. First ed., Dover
Publications; 1968.
[15] C.R. Rao, “Information and Accuracy Attainable in the Estimation of
Statistical Parameters”. Bulletin of the Calcutta Mathematical Society,
Vol. 37, pp. 81–91, 1945.
[16] J.R. Magnus, and H. Neudecker, Matrix Differential Calculus with
Applications in Statistics and Econometrics. John Wiley & Sons, Third
edition; 2007.
[17] M.L. Curtis, Matrix Groups. Springer Verlag, First edition; 1979.
[18] N.J. Higham, Functions of Matrices: Theory and Computation. SIAM,
First edition; 2008.
[19] R.A. Fisher, “Theory of Statistical Estimation”. Mathematical Proceedings of the Cambridge Philosophical Society, Vol. 22, No. 5, pp.
700–725, 1925.

[24] A. Dembo, T.M. Cover, and J.A. Thomas, “Information Theoretic
Inequalities”. IEEE Transactions on Information Theory, Vol. 37, No.
6, pp. 1501–1518, 1991.
[25] J.-F. Bercher, “On Generalized Cramér-Rao Inequalities, Generalized
Fisher Information and Characterizations of Generalized q-Gaussian
Distributions”. Journal of Physics A: Mathematical and Theoretical,
Vol. 45, No. 25, pp. 255303, 2012.
[26] 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.
[27] M.B. Ruskai, “Lieb’s Simple Proof of Concavity of (A, B) 7→
TrAp K † B 1−p K and Remarks on Related Inequalities.” International Journal of Quantum Information, Vol. 3, No. 3, pp. 579–590,
2005. arXiv:quant-ph/0404126, available online: https://arxiv.
org/pdf/quant-ph/0404126.pdf.
[28] T.T. Georgiou, “Relative Entropy and the Multivariable Multidimensional Moment Problem”. IEEE Transactions on Information Theory,
Vol. 52, No. 3, pp. 1052–1066, 2006.
[29] R. Bhatia, “The Logarithmic Mean”. Resonance, Vol. 13, No. 6, pp.
583–594, 2008.
[30] F. Hiai, “Matrix Analysis: Matrix Monotone Functions, Matrix Means,
and Majorization”. Interdisciplinary Information Sciences, Vol. 16, No.
2, pp. 139–248, 2010.

