Nonlinear Kalman filter based on duality relations between continuous and
discrete-state stochastic processes
Jun Ohkubo

arXiv:1402.2025v2 [] 24 Sep 2015

Graduate School of Science and Engineering, Saitama University,
255 Shimo-Okubo, Sakura, Saitama, 338-8570, Japan
(Dated: January 22, 2018)
A new application of duality relations of stochastic processes is demonstrated. Although conventional usages of the duality relations need analytical solutions for the dual processes, we here employ
numerical solutions of the dual processes and investigate the usefulness. As a demonstration, estimation problems of hidden variables in stochastic differential equations are discussed. Employing
algebraic probability theory, a little complicated birth-death process is derived from the stochastic
differential equations, and an estimation method based on the ensemble Kalman filter is proposed.
As a result, the possibility for making faster computational algorithms based on the duality concepts
is shown.
I.

INTRODUCTION

Duality relations between stochastic processes have
been investigated mainly in mathematical physics, and
it has been shown that the duality concepts are useful to
investigate nonequilibrium systems [1]. Ranging from the
asymmetric simple exclusion processes (for example, see
Ref. [2]) to stochastic differential equations (for example,
see Ref. [3]), the duality concepts have been widely employed. As for the duality concepts, see a recent review
paper in Ref. [4].
The previous works for the duality concepts are based
on analytical solutions for the dual stochastic processes;
when the original process is intractable analytically, its
tractable dual process is solved. However, it is still unclear whether the duality concepts can be combined with
numerical methods and work well or not. Recent developments for duality relations for stochastic processes
enable us to derive dual birth-death processes from various stochastic differential equations [5–8], and the derived dual birth-death processes are sometimes analytically intractable. Is it possible to utilize such analytically
intractable birth-death processes in useful ways? In addition, is there a merit to use the duality concept?
The main aim of the present paper is to seek the possibility for the combination of the duality concept and numerical methods. As a demonstration, we here consider
a filtering problem; there are unobserved hidden states in
stochastic differential equations, and our task is to estimate the hidden states only from partially observed data.
A new estimation method based on the duality concepts
is proposed, which is based on the ensemble Kalman filter
(EnKF) [9]. As a result, we see that the method based on
the duality concepts can work at least for a simple nonlinear system. In addition, there is a possibility that the
duality concepts enable us to make faster computational
algorithms from a novel viewpoint; we should perform
numerical simulations for the dual birth-death processes
in advance, and the numerical results can be used for
the time-evolution for the original stochastic differential
equations with arbitrary initial conditions. Hence, there

is no need to perform Monte Carlo simulations for the
original stochastic differential equations at each measurement time step; we can reuse the numerical results for
the dual birth-death processes repeatedly.
The present paper is constructed as follows. In Sec. II,
the model used in the present paper is explained. Section
III is a brief review of the EnKF. The main proposal in
the present paper is given in Sec. IV; the derivation of
the dual birth-death process and the usage of the duality
relation are explained. In Sec. V, results of a demonstration of the new algorithm and comparisons with the
EnKF are given. Section VI is for concluding remarks.

II.
A.

MODEL

Time-evolution of the state variables

In the present paper, the following Van der Pol-type
model, which was used for a test of the filtering problem
in Ref. [10], is considered:
(
d
dt x1 (t) = x2 (t) + w1 (t),
(1)
d
2
dt x2 (t) = ǫ(1 − x1 (t) )x2 (t) − x1 (t) + w2 (t),
where wi (t) ∈ R is a zero-mean white Gaussian noise with
a covariance matrix Q ∈ R2×2 . Different from the original Van der Pol model, the model Eq. (1) contains the
noise terms. Here, we assume that the noises in Eq. (1)
are not correlated with each other, and then the covariance matrix is a diagonal matrix; Q = diag[Q11 , Q22 ]. In
the following, the vector x(t) = [x1 (t) x2 (t)]T is sometimes used for notational brevity.

B.

Measurements

The time-evolution of the state variable x obeys
Eq. (1). Here, suppose that only one of the state variable x can be observed, and that the measurement is
performed with certain time intervals. That is, although

2

4

time interval for the measurements is ∆τ obs = 0.2. Finally, we obtain the measurement data depicted in Fig. 1.
The aim of the problem here is to estimate the state
variable x(t) (i.e., not only x2 (t), but also x1 (t)) from
the partially measured data in Fig. 1. Although it is
necessary to estimate the covariance matrix Q for the
noise in the model and the variance R of the measurement
noise in practice, in the present paper, we assume that
these parameters are previously known for simplicity.

y (measurements)
x2

3
2
1
0
-1
-2

III.

BRIEF REVIEW OF ENKF

-3
0

2

4

6

8

10

12

14

t
FIG. 1. Measurement data. Although a trajectory of x(t) is
simulated using Eq. (1), the measurement is only performed
for x2 with the discrete time interval ∆τ obs = 0.2. In addition,
each measurement includes measurement noise (the variance
R = 0.04). We also depict the original trajectory for x2 as
the dashed line.

the time-evolution of the model Eq. (1) is continuous,
the measurement results are obtained only for the discrete times {τ1 , τ2 , · · · }. For simplicity, in the present
paper, we basically consider that the time interval of
the measurements, ∆τkobs ≡ τk − τk−1 , is fixed, i.e.,
∆τkobs = ∆τ obs for all k. Note that it is easy to extend
the methods discussed in the present paper to variable
time interval cases.
In summary, the following measurement procedure at
time τk is employed:
y(τk ) = Hx(τk ) + v(τk ),

(2)

where H = [0 1] and v(τk ) is a zero-mean white Gaussian
noise with variance R. Hence, only some parts of the
state variable x2 are observed with the addition of the
measurement noise.
C.

Data used in the present paper

The discrete version of the model Eq. (1) has been used
in Ref. [10], and hence we here employ the following parameters, which are similar to the previous work: ǫ = 1.0,
Q11 = 0.0262, Q22 = 0.008, and R = 0.04. (Compared
with the work in Ref. [10], we use a little larger measurement noise.) Using these parameters, the data for the
estimation problem is created as follows.
Firstly, the time-evolution in Eq. (1) is simulated using
the first-order Euler-Maruyama scheme [11, 12]; the time
interval for the simulation is 10−4 , and the initial conditions are x1 (0) = 0.2 and x2 (0) = 0.1. Secondly, after
the simulation of the state variable x(t), the measurement procedure is performed; only the state variable x2
is extracted, and the measurement noises are added. The

The most famous method for adaptive estimation of
the hidden states is the Kalman filter (as for this topic,
for example, see Refs. [13, 14].) The original Kalman
filter was formulated for linear systems under Gaussian
noise, and its nonlinear extensions have been studied
well. The extensions include the extended Kalman filter, the particle filter, and the EnKF; for the details of
these filters, for example, see Ref. [14].
The aim of the present paper is to demonstrate the
applicability of the duality concepts. Hence, a new numerical method for the filtering problem, which is based
on the duality concepts, will be proposed and discussed
later. As a method for comparison, we here employ the
EnKF [9, 15, 16] because the new method, which will be
introduced in Sec.IV, is based on the EnKF. Hence, we
here briefly explain the EnKF.
The EnKF uses the Monte Carlo simulations; using an
ensemble of many particles, statistical quantities such as
means and covariances are evaluated, and these quantities are employed to calculate an important quantity,
so-called Kalman gain. Note that we here use the problem settings in Sec. II, and hence the measurements are
performed only for the discrete times {τ1 , τ2 , · · · }.
Algorithm for the EnKF
1. Initialization:
Make the initial ensemble. Here, we choose n samples from a Gaussian distribution with mean x(0)
and covariance matrix P (0), where x(0) and P (0)
are chosen arbitrarily. We denote each sample at
time τ0 = 0 as xi (τ0 ) for i = 1, 2, · · · , n.
2. Forecast step (xi (τk−1 ) → xfi (τk )):
Using the time-evolution of Eq. (1), simulate the
path of the state variable for each sample starting
from xi (τk−1 ). For the simulation, for example,
the first-order Euler-Maruyama scheme is available.
The simulated path for sample i is denoted as xfi (t)
(τk−1 < t ≤ τk ).
3. Assimilation step (xfi (τk ) → xi (τk )):
(a) Make realizations of random variables
{vi (τk )}ni=1 as the measurement noises. Each
realization is obtained from the zero-mean
white Gaussian noise with the variance R.

3
(b) Calculate following quantities:
(mean evaluated from the ensemble)

IV.

ALGORITHM BASED ON DUALITY
RELATIONS

n

1X f
x (τk ) =
x (τk ),
n i=1 i
f

(error matrix E f (τk ) ∈ R2×n )
E f (τk )
= [xf1 (τk ) − xf (τk )

···

xfn (τk ) − xf (τk )],

(4)

(unbiased covariance evaluated from the ensemble)
Pb f (τk ) =

1
E f (τk )(E f (τk ))T ,
n−1

(5)

(mean of the measurement noises)
n

v(τk ) =

1X
vi (τk ),
n i=1

A.

(3)

(6)

Here, a simple explanation for the basic concept of the
duality relation is given; we consider here the duality relation between a stochastic differential equation and a
birth-death process. For simplicity, stochastic processes
with only one random variable are considered in this subsection.
Suppose that (xt )t≥0 ∈ R is a sample trajectory of the
stochastic differential equation, and p(x, t) is the probability density at time t. It has been known that some
stochastic differential equations have the corresponding
dual birth-death processes. Denote a dual birth-death
process as (nt )t≥0 ∈ N, whose probability distribution
at time t is P (n, t), and then the following equality is
satisfied [1, 5–7]:
Ex [xnt 0 ] = En [xn0 t ] ,

(unbiased variance of the measurement noise)
n

b k) =
R(τ

1 X
(vi (τk ) − v(τk ))2 ,
n − 1 i=1

(Kalman gain)

−1
b k)
b k ) = Pb f (τk )H T H Pb f (τk )H T + R(τ
K(τ
.

(7)

(8)

4. Modify the forecasted state variables
{xfi (τk )}ni=1 using the measurement at
time τk , i.e., y(τk ), as follows:

b k ) y(τk ) + vi (τk ) − Hxfi (τk ) .
xi (τk ) = xfi (τk ) + K(τ
(9)

Steps 2 and 3 in the above algorithm are performed for
each measurement time step.
In the EnKF, the time evolution in the forecast step
is performed as the nonlinear systems, which gives the
non-Gaussian distribution for {xi (t)}ni=1 even if we start
from a Gaussian distribution. After the time evolution,
at each assimilation step, the conventional Kalman filter
is employed, which means a filter at least up to the second
moment. The final filtered value of the state variable x(t)
is obtained as the mean of the ensemble {xi (t)}ni=1 . In
addition, the estimated error could be obtained from the
(co)variance of the ensemble {xi (t)}ni=1 .
One of the problem in the EnKF is as follows: in order
to obtain more accurate estimation results, we need large
ensemble size n. That is, a small ensemble size gives an
inaccurate Kalman gain, and hence the estimation results
would not be accurate. In general, it is easy to imagine
that the large ensemble size needs high computational
costs. Although it has been clarified that a small ensemble size is often enough for the EnKF in practical cases
[10], it would be preferable if we could avoid the numerical simulations of the time evolution of the ensemble at
each measurement time step.

Basic concept

(10)

where Ex and En are the expectations in the stochastic differential equation (xt )t≥0 starting from x0 and in
the birth-death process (nt )t≥0 starting from n0 , respectively. More explicitly, we can rewrite Eq. (10) as
Z

∞

−∞

n0

p(x, t)x dx =

∞
X

P (n, t)xn0 ,

(11)

n=0

where p(x, 0) = δ(x − x0 ) and P (n, 0) = δn,n0 .
Equation (10) shows that the information about the
stochastic differential equation can be obtained from the
solution of the birth-death process. That is, when we obtain the probability distribution of the birth-death process, P (n, t), with the initial condition n0 = 1, it is possible to evaluate the first order moment, i.e., the mean
value of xt , of the stochastic differential equation, without solving the stochastic differential equation.
There are several advantages of the usage of the duality
relations. It is sometimes easier to treat the birth-death
process, compared with the stochastic processes. For
some specific cases, the analytical solution of the birthdeath process has been obtained. In addition, there are
numerical algorithms to simulate the birth-death process
efficiently. As for numerical evaluations of the stochastic
differential equations, we need some approximation; for
example, the time-discretization is needed in the EulerMaruyama scheme. On the other hand, for example, the
Gillespie algorithm for the birth-death process does not
need the time-discretization [17]. In this sense, the birthdeath process would be more tractable than the stochastic differential equations.
Furthermore, only ‘a’ solution of the birth-death process can be used to obtain the information about the
stochastic differential equations with ‘arbitrary’ initial
conditions.
That is, if n0 = 1, Ex [xt ] = En [xn0 t ] =
P∞
n
n=0 P (n, t)x0 , and P (n, t) is independent of the value

4
of x0 . Hence, using only a solution of the birth-death
process, P (n, t), it is possible to estimate the average
of xt , Ex [xt ], for the stochastic differential equations for
‘arbitrary’ initial conditions x0 . Of course, if we want
to know the first and second moments of xt in Eq. (10),
two solutions of the birth-death process nt with the different initial conditions, n0 = 1 and n0 = 2, are needed.
However, once we have these two solutions of the birthdeath process, the first and second moments, Ex [xt ] and
Ex [x2t ], with arbitrary initial conditions can be evaluated
by using the duality relation in Eq. (10). In contrast,
when these moments are evaluated from the direct simulation of the stochastic differential equation, we need
many sample trajectories with an initial condition x0 ;
if the initial condition x0 is changed, we must perform
many other numerical simulations.
The basic idea for a new algorithm is as follows; the
Monte Carlo simulation in the forecast step in the EnKF
is replaced with the simple numerical evaluation based
on the duality relation.
The remaining problem is as follows: How should we
derive the dual birth-death process from a given stochastic differential equation? In the successive subsections,
we will show the method to obtain the dual birth-death
process, employing a mathematical formalism called the
Doi-Peliti formalism [18–21].
B.

Doi-Peliti formalism

We here briefly review the Doi-Peliti formalism, which
is useful to obtain the duality relations. The Doi-Peliti
formalism is the method similar to the second quantization method in quantum mechanics. Up to now, the
Doi-Peliti formalism has been used in various contexts,
mainly in order to investigate discrete systems such as
chemical reactions, and it has been shown that the algebraic probability theory [22] gives the mathematical basis
of the Doi-Peliti formalism [23].
In the Doi-Peliti formalism, creation operator a† and
annihilation operator a are introduced, which satisfy the
following commutation relation:
[a, a† ] ≡ aa† − a† a = 1,

[a, a] = [a† , a† ] = 0.

(12)

These operators act on a vector in the Fock space, |ni,
as follows:
a† |ni = |n + 1i,

a|ni = n|n − 1i,

(13)

and the vacuum state |0i is characterized by a|0i = 0.
Additionally, vectors {hm|}∞
m=0 satisfy the following orthogonal relation to the vectors {|ni}∞
n=0 :
hm|ni = δm,n n!.

(14)

It has been shown that the Doi-Peliti formalism is
deeply related to the conventional generating function
method [24], and the following correspondences can be
useful to understand the usage of the Doi-Peliti formalism in the duality problem:
|ni ↔ xn ,

a† ↔ x,

a↔

∂
.
∂x

(16)

That is, the differential operator is connected to the annihilation operator, which is used in partial differential
equations derived from stochastic differential equations.
In addition, the annihilation operator acts on the discrete states |ni, which is available to construct the birthdeath process, as shown later. The important point here
is that the Doi-Peliti formalism can bridge continuous
states with discrete states, which corresponds to the connection between a stochastic differential equation and a
birth-death process.

C.

Derivation of the dual birth-death process

The derivation of the simple duality relation, such as
Eq. (10), has been discussed in Ref. [5], and the derivation
based on the Doi-Peliti formalism has also been proposed
[6]. However, in order to treat the estimation problem,
only the simple duality relation in Eq. (10) is not enough;
the simple duality relation in Eq. (10) can deal with only
a very restricted class of stochastic differential equations.
Recently, extended duality relations have been proposed
[7], which is necessary to construct the new algorithm
based on the duality relation. Here, we only show, as
an example, the derivation of a dual birth-death process
from the stochastic differential equations in Eq. (1). For
the mathematical details, see the original paper [7].
First of all, it is needed to construct the corresponding Fokker-Planck equation of the stochastic differential
equations in Eq. (1). (For the derivation of the FokkerPlanck equation from the stochastic differential equations, see, for example, Ref. [12].) The corresponding
Fokker-Planck equation is as follows:
∂
p(x1 , x2 , t)
∂t 


∂
∂ 
2
= −
x2 −
ǫ(1 − x1 )x2 − x1 p(x1 , x2 , t)
∂x1
∂x2


∂2
∂2
1
Q11 2 + Q22 2 p(x1 , x2 , t),
(17)
+
2
∂x1
∂x2

Note that a† a corresponds to the number operator and
a† a|ni = n|ni.

(15)

where p(x1 , x2 , t) is the probability density at time t for
the stochastic differential equations. For notational con-

5
venience, we introduce the following linear operator L∗ :


∂
∂
L∗ x1 ,
, x2 ,
∂x1
∂x2



∂ 
∂
x2 −
ǫ(1 − x21 )x2 − x1
= −
∂x1
∂x2


2
∂2
∂
1
Q11 2 + Q22 2 ,
+
2
∂x1
∂x2

(18)

and hence the Fokker-Planck equation (17) is rewritten
as


∂
∂
∂
∗
p(x1 , x2 , t).
, x2 ,
p(x1 , x2 , t) = L x1 ,
∂t
∂x1
∂x2
(19)
The adjoint operator of Eq. (18), L, is as follows:


∂
∂
, x2 ,
L x1 ,
∂x1
∂x2
∂
∂
∂
= x2
− ǫ(1 − x21 )x2
− x1
∂x1
∂x2
∂x2


1
∂2
∂2
+
Q11 2 + Q22 2 ,
2
∂x1
∂x2

(20)

(21)

where the following correspondences are used:
a1 ↔

∂
,
∂x1

a†2 ↔ x2 ,

a2 ↔

∂
.
∂x2

(22)

∞ X
∞ X
∞
X

P (n0 , n1 , n2 , t̃)|n0 , n1 , n2 i,

where P (n0 , n1 , n2 , t̃) is a probability distribution of a
birth-death process, and considering the time-evolution
equation
∂
|ψ(t̃)i = L(a†0 , a†1 , a1 , a†2 , a2 )|ψ(t̃)i,
∂ t̃
we may have a time-evolution equation for
P (n0 , n1 , n2 , t̃) by comparing the coefficient of a
state vector |n0 , n1 , n2 i on the right and left hand sides.
However, as discussed in the previous work [7], it is
impossible to simply interpret the linear operator in
Eq. (25) as a time-evolution operator for a birth-death
process; the linear operator does not satisfy the probability conservation law, and, in addition, some terms
seem to correspond to ‘negative’ transition rates.
In order to construct an adequate birth-death process,
we need an additional operator b, which satisfies the following relations:
b|+i = |−i,

In addition, as discussed later, it is convenient to introduce the time scaling t = rts t̃; due to the time scaling,
the original Fokker-Planck equation (19) is rewritten as
∂
p(x1 , x2 , rts t̃)
∂ t̃


∂
∂
∗
p(x1 , x2 , rts t̃).
, x2 ,
= rts L x1 ,
∂x1
∂x2

Note that all creation operators in each term in Eq. (25)
are placed on the left side of the annihilation operators;
if not, we must replace the term using the commutation
relation in Eq. (12).
Since the linear operator in Eq. (25) acts on the Fock
space, i.e., the discrete state space {|ni}∞
n=0 , one may
expect that the linear operator simply gives the timeevolution for a birth-death process. That is, defining the
state vector |ψ(t̃)i as

n0 =0 n1 =0 n2 =0



L a†1 , a1 , a†2 , a2

a†1 ↔ x1 ,

= a†0 a†2 a1 + ǫa†0 a†2 a2 − ǫa†0 (a†1 )2 a†2 a2 − a†0 a†1 a2
i
1h
+
(25)
Q11 a†0 a1 a1 + Q22 a†0 a2 a2 .
2

|ψ(t̃)i ≡

and using the correspondence between operators in the
Doi-Peliti formalism and the differential operators in
Eq. (16), we have

= a†2 a1 + ǫ(1 − a†1 a†1 )a†2 a2 − a†1 a2
1
+ [Q11 a1 a1 + Q22 a2 a2 ] ,
2

and therefore the following linear operator is obtained:


L a†0 , a†1 , a1 , a†2 , a2


≡ a†0 L a†1 , a1 , a†2 , a2

b|−i = |+i,

(26)

where |+i and |−i are orthonormal state vectors satisfying h+|+i = h−|−i = 1 and h+|−i = h−|+i = 0.
The following procedure is needed to construct an adequate time-evolution operator for a birth-death process
(although the following procedure might seem complicated, a concrete example will be given soon):

(23)

For each term in the linear operator L, apply the
following procedures.

We here focus on the fact that continuous variables, x1
and x2 , are replaced with creation operators in Eq. (22).
Hence, if we reinterpret a constant as a creation operator,
the following replacement is available:

1. If the coefficient of the term has a negative sign,
replace the negative sign ‘−’ with ‘+b’ using the
operator b.

rts ↔ a†0 ,

(24)

2. Act the operators on |ni, and evaluate the coefficient including the effects of the number operators
a†i ai .

6
3. Subtract a term in order to guarantee the probability conservation law; the term gives the same
coefficient with Step 2, and additionally, the term
consists of the same number of creation and annihilation operators for the same sub-index.
4. In order to compensate the subtracted term, add
the same term with Step 3.
For example, the first term in Eq. (25) is a†0 a†2 a1 , and
then we do not need Step 1. Since
a†0 a†2 a1 |n0 , n1 , n2 i = n1 |n0 + 1, n1 − 1, n2 + 1i,

(27)

the coefficient is n1 , and therefore the term a†1 a1 must
be subtracted in Step 2. On the other hand, the third
term in Eq. (25), −ǫa†0 (a†1 )2 a†2 a2 , has the negative sign,
and hence Step 1 gives +bǫa†0 (a†1 )2 a†2 a2 . The subtracted
term in Step 3 is bǫa†2 a2 because
a†0 (a†1 )2 a†2 a2 |n0 , n1 , n2 i = n2 |n0 + 1, n1 + 2, n2 i,

(28)

which gives the coefficient n2 .
Using the above procedure, we obtain the following linear operator, instead of the linear operator L in Eq. (25):


L a†0 , a0 , a†1 , a1 , a†2 , a2 , b




= L′ a†0 , a0 , a†1 , a1 , a†2 , a2 , b + V a†1 a1 , a†2 a2 ,

(29)

where



L′ a†0 , a0 , a†1 , a1 , a†2 , a2 , b

= (a†0 a†2 a1 − a†1 a1 ) + (ǫa†0 a†2 a2 − ǫa†2 a2 )
+ (ǫba†0 (a†1 )2 a†2 a2 − ǫa†2 a2 ) + (ba†0 a†1 a2 − a†2 a2 )


1
1
†
† †
Q11 a0 a1 a1 − Q11 a1 a1 a1 a1
+
2
2


1
1
†
† †
(30)
Q22 a0 a2 a2 − Q22 a2 a2 a2 a2
+
2
2

and


V a†1 a1 , a†2 a2



1
= a†1 a1 + 2ǫa†2 a2 + a†2 a2 + Q11 (a†1 a1 )(a†1 a1 ) − a†1 a1
2


1
†
†
+ Q22 (a2 a2 )(a2 a2 ) − a†2 a2 .
(31)
2

Note that we used a† a† aa = (a† a)(a† a)−aa, which stems
from [a, a† ] = aa† − a† a = 1.
As shown in the previous work [7], it is possible to
interpret the linear operator L′ in Eq. (30) as a timeevolution operator for a birth-death process. Because

of the operator b, we must consider an additional state
variable, which takes only two states (+ or −), in addition to the state variables n0 , n1 and n2 . See the first
term in Eq. (30); the action of the first term on |ni is
Eq. (27), and hence we can interpret this term as an elementary birth-death process X1 → X2 +X0 with rate n1 ,
where n0 , n1 and n2 is the number of particles X0 , X1 ,
and X2 , respectively. On the other hand, the third term
in Eq. (30) gives X2 → 2X1 + X2 + X0 and the state
change with + → − or − → + with rate ǫn2 ; the state
(+ or −) is changed due to an event corresponding to the
third term. Repeating the similar discussions, finally the
following birth-death process is obtained:

n1

X2 + X0 ,
X1 −→
(i)


ǫn
2

(ii) X2 −−→ X2 + X0 ,



ǫn2

(iii) X2 −
−→ 2X1 + X2 + X0 with S.C.,
n2
(32)
−
(iv)
X
2 → X1 + X0 with S.C.,



n1 (n1 −1)


(v) 2X1 −−−−−−→ X0 ,




n2 (n2 −1)
(vi) 2X −
→ X0 ,
2 −−−−−

where ‘S.C.’ means the state change + → − or − → +.
The term V (a†1 a1 , a†2 a2 ) is called a Feynman-Kac term,
and the important point is that this term is written only
in terms of the number operators a†i ai . Since the number
operators a†i ai does not affect the state vectors, we can
simply replace the number operators with the random
variables in the birth-death process; i.e., V (a†1 a1 , a†2 a2 ) =
V (n1 , n2 ).
The intuitive understanding of the duality relation is
written as follows: we here abbreviate a state related to
the Fokker-Plank equation at time t as FP(t) and that to
the birth-death process as BD(t). In addition, for simplicity, set rts = 1 and hence t = t̃; furthermore, we neglect the Feynman-Kac term here. Then, formally, timeevolution of BD(t) is given by BD(t) = eLt BD(0). In
contrast, that of FP(t) is written as FP(t) = FP(0)eLt ,
when we consider the left-action of L. Hence, the adjoint operator L∗ gives the actual time-evolution of the
Fokker-Planck equation. In addition, we have formally
FP(t)BD(0) = FP(0)eLt BD(0) = FP(0)BD(t); this corresponds to the duality relation. The linear operator L
can be written in terms of both the differential operators
and the creation-annihilation operators, and hence the
stochastic differential equation in Eq. (1) and the birthdeath process in Eq. (32) are connected naturally. Of
course, the above discussion is just an intuitive one, and
for the mathematical explanation of the duality relations,
see the previous work [7].
Here we consider a time-evolution from time 0 to τ .
Because we consider the time-scaling factor rts , the time
interval, τ , in the stochastic differential equations correspond to the time interval τ̃ = τ /rts in the dual birthdeath process. Using the duality relation, we can obtain
the following identity:

7

h
i
Ex x1 (τ )n1 (0) x2 (τ )n2 (0)
"
(Z
τ̃

= En,+ exp
"

− En,− exp

V (n1 (s), n2 (s))ds

0

(Z

)

τ̃

V (n1 (s), n2 (s))ds

0

where the abbreviations x = (x1 , x2 ) and n =
(n0 , n1 , n2 ) are used; we set n0 (0) = 0; En,+ and En,−
are the expectations related to the states + and − in
the birth-death process in Eq. (32), respectively. Using
the duality relation in Eq. (33), we can evaluate various information for the stochastic differential equation
in Eq. (1) from the solution of the birth-death process in
Eq. (32). In order to evaluate the probability distribution and the corresponding Feynman-Kac terms for the
birth-death process in Eq. (32), Monte Carlo simulations
with the Gillespie algorithm are available.
As noted in Sec. IV.A, the initial conditions for the
birth-death process correspond to the order of moments
in the stochastic differential equation. When we use
n1 (0) = 2 and n2 (0) = 0, the time-evolution of the second moment for x1 , i.e., Ex [x1 (τ )2 ], is evaluated.

n (τ̃ )
rts0 x1 (0)n1 (τ̃ ) x2 (0)n2 (τ̃ )

)

n (τ̃ )
rts0 x1 (0)n1 (τ̃ ) x2 (0)n2 (τ̃ )

At the end of this subsection, we give one comment for
the operator b. In the above construction, we replaced
the minus sign ‘−’ with the operator ‘+b’ and avoided the
‘negative’ transition problem. Instead of that, we can use
the following trick: Interpret constants as random variables. That is, we interpret ‘−1’ as a random variable x3 ,
and the (stochastic) differential equation for x3 is set as
dx3 /dt = 0. Then, employing the same discussion above,
x3 is replaced with a†3 in the Doi-Peliti formulation, and

#

,

(33)

finally we obtain the following birth-death process

n1

X2 + X0 ,
(i)
X1 −→



ǫn2

−
−
→
X2 + X0 ,
(ii)
X

2


ǫn2

(iii) X2 −
−→ 2X1 + X2 + X0 + X3 ,
n2
(34)
−
(iv)
X
2 → X1 + X0 + X3 ,



n1 (n1 −1)


(v) 2X1 −−−−−−→ X0 ,



n2 (n2 −1)
(vi) 2X −
−−−−−→ X ,
2

0

and the following duality relation
i
h
Ex x1 (τ )n1 (0) x2 (τ )n2 (0) x3 (τ )n3 (0)
"
(Z
)
τ̃
= En exp
V (n1 (s), n2 (s))ds
0

×
Using the above procedures, the time evolution of the
stochastic differential equation for arbitrary initial conditions can be replaced with that of the dual birth-death
process for specific initial conditions. In addition, if
changing the time-scaling variable rts , we can evaluate
information for various time interval in the stochastic
differential equations from only a result for a single fixed
time interval τ̃ in the birth-death process. For example,
assume that τ̃ = 1; if we want to know the information
of the stochastic differential equation at time t = 1, we
should set rts = 1. On the other hand, the information
at time t = 0.9 can be evaluated from the same results of
the birth-death process with τ̃ = 1 and setting rts = 0.9.
Hence, there is no need to perform the Monte Carlo simulations for various different time intervals; this is also
one of the advantages of the new method.

#

n (τ̃ )
rts0 x1 (0)n1 (τ̃ ) x2 (0)n2 (τ̃ ) x3 (0)n2 (τ̃ )

#

,

(35)

where the abbreviations x = (x1 , x2 , x3 ) and n =
(n0 , n1 , n2 , n3 ) are used. Note that the Feynman-Kac
term does not depend on n3 (t). Because the ‘random’
variable x3 (t) is a time-independent constant and x3 (t) ≡
−1, we can see that the initial condition, n0 (0) = 0
and n3 (0) = 0, recovers the original duality relation in
Eq. (33). This technique, in which constants are interpreted as random variables, would be sometimes useful
to consider more complicated stochastic differential equations. However, in the example considered here, only the
‘negative’ transition rate should be avoided, and the expression in Eq. (33) is simpler in computational viewpoint; although n3 (t) takes any natural numbers, the additional states |+i and |−i takes only two states, and
hence the computational memory in Monte Carlo simulations is largely reduced. From this reason, the additional
operator b is introduced and used in the present paper.
D.

Moment evaluation

The duality relation in Eq. (33) can be used for evaluating the moments in stochastic differential equation in
Eq. (1) starting from the Dirac-delta-type initial conditions, as discussed in Sec. IV.A. However, in the EnKF,
we use an ensemble of samples. In the EnKF, only the

8
mean and covariance matrix of the ensemble are needed,
so we can consider that the ensemble is essentially characterized by a Gaussian distribution. Hence, the Diracdelta-type initial conditions are not enough.
We here write the expectation value, which is taken
by using a Gaussian distribution with mean x(t) and
covariance matrix V , as h· · · i. Hence, the first term in
the r.h.s. in Eq. (33) should be replaced with
iE
h R τ̃
D
n (τ̃ )
En,+ e 0 V ds rts0 x1 (0)n1 (τ̃ ) x2 (0)n2 (τ̃ )
h R τ̃
D
Ei
n (τ̃ )
= En,+ e 0 V ds rts0
x1 (0)n1 (τ̃ ) x2 (0)n2 (τ̃ ) .
We need a little additional calculations in order to evaluate the expectation with the Gaussian distribution; introducing the following notations,
µ1 = hx1 i, µ2 = hx2 i,

(36)

the following recursion formula is known [25]:
hxn1 1 xn2 2 i =µ1 xn1 1 −1 xn2 2 + V11 (n1 − 1) xn1 1 −2 xn2 2
+ V12 n2 xn1 1 −1 xn2 2 −1

(37)

or
hxn1 1 xn2 2 i =µ2 xn1 1 xn2 2 −1 + V21 n1 xn1 1 −1 xn2 2 −1
+ V22 (n2 − 1) xn1 1 xn2 2 −2 .

(38)

Using the above recursion formula, once we have the
probability distribution P (n0 , n1 , n2 , t) for the birthdeath process in Eq. (32) with adequate initial values, it
is possible to evaluate various moments for the stochastic differential equation in Eq. (1) with a Gaussian initial
distribution, using the following duality relation:

i
h
Ex,Gaussian initial x1 (τ )n1 (0) x2 (τ )n2 (0)
"
(Z
#
)
τ̃
D
E
n0 (τ̃ )
n1 (τ̃ )
n2 (τ̃ )
= En,+ exp
V (n1 (s), n2 (s))ds rts
x1 (0)
x2 (0)
"

− En,− exp

E.

0

(Z

τ̃

V (n1 (s), n2 (s))ds
0

Algorithm

We here consider general cases with variable time intervals, and denote the maximum of the time intervals,
τk−1 − τk , as ∆τ . The new Kalman filter based on the
duality relation, called the DuKF here, is as follows:
Algorithm for DuKF
1. Preparation for the duality relations:
Simulate the dual birth-death process in Eq. (32).
We need simulations with five different initial conditions;
(c1) n0 = 0, n1 = 1, n2 = 0,
(c2) n0 = 0, n1 = 0, n2 = 1,
(c3) n0 = 0, n1 = 2, n2 = 0,
(c4) n0 = 0, n1 = 0, n2 = 2,
(c5) n0 = 0, n1 = 1, n2 = 1,
which are necessary to evaluate Ex [x1 (t)],
Ex [x2 (t)], Ex [x1 (t)2 ], Ex [x2 (t)2 ], Ex [x1 (t)x2 (t)],
respectively. For all cases, the additional state
variable is set to ‘+’ initially. The Monte Carlo
simulations are performed from t = 0 to t = ∆τ ,
and the integral of the Feynman-Kac term

)

n (τ̃ )
rts0

#
E
D
n1 (τ̃ )
n2 (τ̃ )
.
x1 (0)
x2 (0)

(39)

V (n1 , n2 ) and the final probability distribution
P (n0 , n1 , n2 , ∆τ ) are evaluated numerically.
2. Initialization:
Set an initial Gaussian distribution with mean
xa (0) and covariance matrix P a (0).
3. Forecast step at τk−1 :
Using the duality relation in Eq. (39), evaluate various moments, which are necessary to characterize
a Gaussian distribution. The ensemble average in
Eq. (39), h· · · i, is taken for the Gaussian distribution with mean xa (τk−1 ) and covariance matrix
P a (τk−1 ). The above procedure gives the mean
xf (τk ) and the covariance Pb f (τk ) of the nonlinear
systems at time τk . Note that the time-scaling variable rts must be selected adequately according to
the ratio between τk − τk−1 and ∆τ .
4. Assimilation step at τk :
Calculate the Kalman gain

−1
b k ) = Pbf (τk )H T H Pb f (τk )H T + R
,
K(τ

(40)

where H is the matrix assigned for the measurement and R the variance of the measurement;

9

4

(a)

x(0) = (0.1, 0.1), [P (0)]11 = 0.1, [P (0)]12 = [P (0)]21 = 0,
and [P (0)]22 = 0.1. Firstly, for the DuKF, we performed
the Monte Carlo simulations for the dual birth-death process using the Gillespie algorithm [17]. For each initial
condition, 1012 sample paths were generated, and the
Feynman-Kac term and the probability distribution were
evaluated.

true
EnKF
DuKF

3

x1

2
1
0
-1
-2
0

2

4

6

8

10

12

14

12

14

t

(b)

true
EnKF
DuKF

4
3

x2

2
1
0
-1
-2
-3
0

2

4

6

8

10

t
FIG. 2. State estimates x1 and x2 for the stochastic differential equations (1); (a) for the unobserved state variable x1 ,
and (b) for the state variable x2 , which is observed with noise.
Note that x2 is measured only at discrete times as shown in
Fig. 1. For the EnKF, we set n = 10 and ∆t = 10−4 in the
Euler-Maruyama scheme.

see Eq. (2). Update the following quantities:

b k ) y(τk ) − Hxf (τk )
xa (τk ) = xf (τk ) + K(τ
(41)
and

b k )H)Pbf (τk ).
P a (τk ) = (1 − K(τ

(42)

Steps 3 and 4 are performed for each measurement time
step.
V.

RESULTS

In order to demonstrate the DuKF, we employ the
algorithm in Sec. IV.E to the problem in Sec. II. For
simplicity, we assume here that the parameters in the
stochastic differential equations, Q and R, are previously known, as explained in Sec. II. We use the following parameters for the initial Gaussian distribution;

Figure 2 shows the results of the state estimation. Note
that only the state variable x2 is observed at discrete
times, as depicted in Fig. 1. Although the observation
contains the measurement noises, as seen in Fig. 1, both
the EnKF and DuKF give adequate estimations. Especially, the non-observed state x1 can be estimated reasonably as shown in Fig. 2. Note that the initial guesses
for x are not far from the true values, and hence it is
difficult to see the initial transient behavior for both the
EnKF and DuKF. From Fig. 2, it may be difficult to
judge whether the DuKF gives better results than the
EnKF or not, but we confirmed that the DuKF gives a
slightly better mean squared error. If we want to obtain
the similar mean squared error, we need larger ensemble
for the EnKF, and the larger ensemble needs more computational time. On the other hands, the forecast step
in the DuKF does not need any Monte Carlo simulation,
and hence the DuKF works rapidly. Actually, in order to
deal with the current data, the EnKF with ensemble size
n = 10 and original discretized time step for simulations
(10−4 ) needs about 0.8 seconds in a standard computer
with Intel Core i5 processor (2.2GHz). In contrast, the
DuKF needs less than 0.1 seconds. As shown below, the
ensemble size n = 10 is not enough large, and if we use
n = 1000 ensembles, 1000 times costs of simulations are
needed for the EnKF. (Of course, the DuKF needs precalculations for the dual birth-death processes. For 1012
sample paths, about one month calculations were needed,
but we can reuse the pre-calculation results repeatedly in
the actual filtering steps.)

In order to see the difference between the EnKF and
DuKF more explicitly, we show the covariance [P (t)]12
calculated in the EnKF and DuKF in Fig. 3. As shown
in Fig. 3, the larger ensemble size in the EnKF gives
the similar covariance with that of the DuKF. The covariance is used to evaluate the Kalman gain, and hence
more accurate covariance is necessary to have better estimations. The ensemble size of order 103 is needed to
obtain the similar covariance with that of the DuKF; it is
very time-consuming. In addition, it would be needed to
choose an adequate time-interval ∆t for the simulation in
the Euler-Maruyama scheme for the stochastic differential equations; the computational time and the precision
of the estimations in EnKF largely depend on ∆t. On
the other hand, we do not need such time discretization
for the DuKF when the Gillespie algorithm is employed,
as discussed before.

covx

1 x2

10

0.06
0.05
0.04
0.03
0.02
0.01
0.00
-0.01
-0.02
-0.03

EnKF, n = 10
EnKF, n = 1000
DuKF

0

2

4

6

8

10

12

14

t
FIG. 3. Covariance between x1 and x2 , [P (t)]12 , calculated
in the EnKF and DuKF. The evaluated covariance is used
to calculate the Kalman gain. The evaluation in the EnKF
approaches to that in the DuKF when larger ensembles are
used.

VI.

CONCLUDING REMARKS

Of course, we do not claim that the proposed DuKF
is always superior to the EnKF for any estimation problem. Actually, the implementation of the EnKF is simpler than that of the DuKF. In addition, there could be
more appropriate methods for specific tasks (see, for example, [26].) In order to make the DuKF more practical,
it is needed to develop more efficient numerical methods
for the dual birth-death processes; very high accuracy
is necessary. Actually, in preliminary works, the famous
Lorenz system was investigated, but the DuKF works
well only for very short observation time intervals, and
then it was not practical. If we have more accurate numerical results, it would be possible to make the DuKF
for the chaotic systems. In order to overcome this problem, it may be possible to employ simulations based on
the important sampling methods [27, 28]. Note that it
has not yet been clarified what factors determine the precalculation time; the preliminary works for the Lorenz
systems could suggest that the chaotic behavior is related to the demands for the high-accuracy of the precalculation, but more detailed studies should be done
in future from both experimental and theoretical view
points for the practical usage of the duality relations.
We hope that the present work open up a new way of
the numerical applications of the duality concepts, and
further efficient numerical methods will be proposed in
future.

The duality relation between stochastic processes has
been still developing, and hence it is expected that new
applications of the duality relations give completely novel
algorithms for various research fields. We demonstrate
one of the applications by using the estimation problem. As a result, pre-calculations for the dual birth-death
processes enable us to avoid the time-consuming Monte
Carlo simulations for each forecast step.

This work was supported in part by grant-in-aid for
scientific research (No. 25870339) from Japan Society
for the Promotion of Science (JSPS).

[1] T.M. Liggett, Interacting particle systems (Classics in
Mathematics) (Springer, Berlin, 2005).
[2] G.M. Schütz, J. Stat. Phys. 86, 1265 (1997).
[3] T. Shiga and K. Uchiyama, Probab. Th. Rel. Fields 73,
87 (1986).
[4] S. Jansen and N. Kurt, Probab. Surv. 11, 59 (2014).
[5] C. Giardinaá, J. Kurchan, F. Redig, and K. Vafayi, J.
Stat. Phys. 135, 25 (2009).
[6] J. Ohkubo, J. Stat. Phys. 139, 454 (2010).
[7] J. Ohkubo, J. Phys. A: Math. Theor. 46, 375004 (2013).
[8] G. Carinci, C. Giardinà, C. Gilberti and F. Redig, Stoch.
Proc. and Appl. 125, 941 (2015).
[9] G. Evensen, J. Geophys. Res. 99, 10143 (1994).
[10] S. Lakshmivarahan and D.J. Stensrud, IEEE Control
Syst. 29, 34 (2009).
[11] P.E. Kloeden and E. Platen, Numerical solution of
stochastic differential equations (Springer, Berlin, 1992).
[12] C. Gardinar, Stochastic methods (Springer, Heidelberg,
2009) 4th ed.
[13] S. Haykin, Adaptive filter theory (Prentice-Hall, New Jersey, 2002) 4th ed.

[14] J.V. Candy, Bayesian signal processing: classical, modern, and particle filtering methods (Wiley, Hoboken, N.J.,
2009).
[15] S. Gillijns, O.B. Mendoza, J. Chandrasekar, B.L.R.
De Moor, D.S. Bernstein, and A. Ridley, What is
the ensemble Kalman filter and how well does it
work?
Proc. 2006 American Control Conference.
doi:10.1109/ACC.2006.1657419.
[16] G. Evensen, IEEE Control Syst. 29, 83 (2009).
[17] D.T. Gillespie, J. Phys. Chem. 81, 2340 (1977).
[18] M. Doi, J. Phys. A: Math. Gen. 9, 1465 (1976).
[19] M. Doi, J. Phys. A: Math. Gen. 9, 1479 (1976).
[20] L. Peliti, J. Physique 46, 1469 (1985).
[21] U.C. Täuber, M. Howard, and B.P. Vollmayr-Lee, J.
Phys. A: Math. Gen. 38, R79 (2005).
[22] A. Hora and N. Obata, Quantum probability and spectral
analysis of graphs (Springer, Berlin, 2007).
[23] J. Ohkubo, J. Phys. Soc. Jpn. 82, 084001 (2013).
[24] M. Droz and A. McKane, J. Phys. A: Math. Gen. 27,
L467 (1994).
[25] R. Willink, Stat. Prob. Lett. 73, 271 (2005).

ACKNOWLEDGEMENT

11
[26] L. Palatella, A. Trevisan and S. Rambaldi, Phys. Rev. E
88, 022901 (2013).
[27] S. Asmussen and P.W. Glynn, Stochastic simulation
(Springer, New York, 2007).

[28] D.P. Landau and K. Binder, A guide to Monte Carlo
simulations in statistical physics (Cambridge University
Press, Cambridge, 2009) 3rd ed.

