Full text of "Hermite spectral method to 1D forward Kolmogorov equation and its application to nonlinear filtering problems"

See other formats

Hermite spectral method to ID forward Kolmogorov
equation and its application to nonlinear filtering

problems

Xue Luo and Stephen S.-T. Yau, ieee Fellow

Abstract — In this paper, we investigate the Hermite spectral
method (HSM) to numerically solve the forward Kolmogorov
equation (FKE). A useful guideline of choosing the scaling
factor of the generalized Hermite functions is given in this
paper. It greatly improves the resolution of HSM. The conver-
gence rate of HSM to FKE is analyzed in the suitable function
space and has been verified by the numerical simulation.
As an important application and our primary motivation to
study the HSM to FKE, we work on the implementation
of the nonlinear filtering (NLF) problem with a real-time
algorithm developed in 1 17 1. The HSM to FKE is served as
the off-line computation in this algorithm. The translating
factor of the generalized Hermite functions and the moving-
window technique are introduced to deal with the drifting
of the posterior conditional density function of the states in
the on-line experiments. Two numerical experiments of NLF
problems are carried out to illustrate the feasibility of our
algorithm. Moreover, our algorithm surpasses the particle
filter as a real-time solver to NLF.

I. Introduction

The central problem in the field of nonlinear filtering
(NLF) is to give the instantaneous and accurate estimation
of the states based on the noisy observations, if enough
computational resources are provided. In 1960s, Duncan
fTl . Mortensen (T8| and Zakai |23| independently derived
the so-called DMZ equation, which the unnormalized
conditional density function of the states satisfies. Hence
the central problem in NLF is translated into solving the
DMZ equation in the real time and memoryless manner. It
is worthy to point out that the "real time" and "memory-
less" are the most important properties one would like to
maintain in the design of the optimal/suboptimal nonlinear
filters for real applications. More specifically, "memory-
less" refers that one only needs the latest observation to
update the estimation of the states without refering back
to any earlier observation history; "real time" means that
the decision of the states is made on the spot, while the
observation data keep coming in.

It is well known that the exact solution to the DMZ
equation, generally speaking, can not be written in a
closed form. With the well-posedness theory of the DMZ

X. Luo is with the department of Mathematics, Statistics and
Computer Science, University of Illinois at Chicago, Science and
Engineering Offices (M/C 249), 851 S. Morgan Street, Chicago,
IL 60607-7045 xluo6@uic.edu

S. S.-T. Yau is the faculty of Department of mathematical
sciences, Tsinghua University, Beijing, 100084, P.R.China. He
is the Emeritus Distinguished Professor of University of Illinois
at Chicago yau@uic . edu

equation in mind, many mathematicians inake efforts to
seek an efficient algorithm to construct a "good" approxi-
mate solution to the DMZ equation. One of the methods is
the splitting-up method from the Trotter product formula,
which was first described in Besoussan, Glowinski, and
Rascanu , |6|. It has been extensively studied in many
articles later, for instance 1121 . 1141 . 1151 and |19| . In
1990s, Lototsky, Mikulevicius and Rozovskii fT6) devel-
oped a new algorithm (so-called S'^-algorithm) based on
the Cameron-Martin version of Wiener chaos expansion.
However, the above methods require the boundedness of
the drifting term and the observation term (/ and h in
l |l.l[ )), which leaves out even the linear case. To overcome
this restriction, Yau and Yau |22| developed a novel
algorithm to solve the "pathwise-robust" DMZ equation,
where the boundedness of the drift term and observation
term is replaced by some mild growth conditions on /
and h. Nevertheless, they still made the assumption that
the drift term, the observation term and the diffusion term
are "time-invariant". That is to say, /, h and g in ([TTTJ
are not explicitly time-dependent. In 1171 . we generalized
Yau-Yau's algorithm to the most general settings of the
NLF problems - "time-varying" case - where /, h and g
are explicitly time-dependent.

Our study of solving the forward Kolmogorov equation
(FKE) by the Hermite spectral method (HSM) is closely
related to the implementation of the algorithm developed
in 1171 . The detailed formulation of our algorithm could
be found in Appendix A or 1 17|. Briefly speaking, in our
algorithm, we start from the signal based model:

dxt = f{xt,t)dt + g{xt, t)dvt,
dyt = h{xt,t)dt + dwt,

where Xt is a vector of the states of the system at time
t with xo satisfying some initial distribution and yt is
a vector of the observations at time t with yo = 0.
vt and wt are vector Brownian motion processes with

E[dvtdvf] = Q{t)dt and E[dwtdw'^] = S{t)dt, S{t) >
0, respectively. The DMZ equation is derived as

da(x, t) = La{x, t)dt + a[x, t)h^ {x, t)S^^ {t)dyt
a{x, 0) = ao{x),

(1.2)

where a{x) is the unnormalized conditional density funci-

HSM to ID FKE and its application to the NLF problems

ton, (Jo{x) is the density of the initial states xq, and

dxidxj

dxi

(1.3)

To maintain the real-time property, solving the DMZ
equation is translated into solving a FKE off-line and
updating the initial data on-line at the beginning of each
time interval. Let Vk = {0 = ro < ri < ■ ■ ■ < ^ T}
be a partition of [0,r]. The FKE needs to be solved at
each time step is

dui ,

(1.4)

where L is defined as l |1.3[ l. The initial data is updated as
follows

-(x, t) = - ^/i^S ^h^Ui{x,t) on[Ti-i,Ti],

ui{x,0) = ao{x)

Wi(a;,Ti_i) = exp ri_i)5 ^(ri_i)

i > 2,

_ilx,ri_il

(1.5)

where Ui is transformed from o, see the detailed for-
mulation of our algorithm in the Appendix A or 1171 .
From the above description, it is not hard to see that the
FKE \\A\ needs to be solved repeatedly on each time
interval [ri_i,ri]. Thus, it is crucial to obtain a good
approximate solution to \\A) . In this paper, we adopt
HSM to solve FKE for two reasons: on the one hand,
HSM is particularly suitable for functions defined on the
unbounded domain which decays exponentially at infinity;
on the other hand, HSM could be easily patched with
the numerical solution obtained in the previous time step
while the moving-window technique is in use in the on-
line experiments.

The HSM itself is also a field of research, which
could be traced back to 1970s. In |10|, Gottlieb et. al.
gave the example sin a; to illustrate the poor resolution
of Hermite polynomials. To resolve M wavelength of
sin a;, it requires nearly K'P Hermite polynomials. Due
to this fact, they doubted the usefulness of Hermite
polynomials as bases. The Hermite functions inherit the
same deficiency from the polynomials. Moreover, it is lack
of fast Hermite transform (some analogue of fast Fourier
transform). Despite of all these drawbacks, the HSM has
its inherent strength. Many physical models need to solve
a differential equation on an unbounded domain, and
the solution decays exponentially at infinity. From the
computational point of view, it is hard to describe the
rate of decay at infinity numerically or to impose some
artificial boundary condition cleverly on some faraway
"boundary". Therefore the Chebyshev or Fourier spectral
methods are not so useful in this situation. As to the HSM,
how to deal with the behavior at infinity is not necessary.
Recent applications of the HSM can be found in [S], (51 ,
fTTl . (20J, tlJJ, etc.

To overcome the poor resolution, a scaling factor is
necessary to be introduced into the Hermite functions,
refer to fS), (4). It is shown in [4| that the scaling factor
should be chosen according to the truncated modes A'^
and the asymptotical behavior of the function f{x), as
\x\ — >■ oo. Some efforts have been made in seeking
the suitable scaling factor a, see

II, QU, (20), etc. To
optimize the scaling factor is still an open problem, even
in the case that fix) is given explicitly, to say nothing
of the exact solution to a differential equation, which is
generally unknown a-priori. Although some investigations
as far as we know, there is no practical guidelines of
choosing a suitable scaling factor. Nearly all the scaling
factors in the papers with the application of HSM are
obtained by the trial-and-error method. Thus, we believe
it is necessary and useful to give a practical strategy to
pick an appropriate scaling factor and the corresponding
truncated mode for at least the most commonly used
types of functions, i.e. the Gaussian type and the super-
Gaussian type functions. The strategy we are about to give
only depends on the asymptotic behavior of the function.
In the scenario where the solution of some differential
equation needs to be approximated (the exact solution is
unknown), we could use asymptotical analysis to obtain
its asymptotic behavior. Thus, our strategy of picking
the suitable scaling factor is still applicable. A numerical
experiment is also included to verify the feasibility of our
strategy. Although it may not be optimal with respect to
the accuracy, our strategy provides a useful guideline for
the implementations of HSM. In this paper, the precise
convergence rate of the HSM to FKE is obtained in
suitable funciton space by numerical analysis and verified
by a numerical example.

Let us draw our attention back to the implementation
of our algorithm to NLF problems. Through our study
of HSM to FKE, the off-line data could be well pre-
pared. However, when synchronizing the off-line data
with the on-line experiments, to be more specifically,
updating the initial data according to \\.5) on-line, another
difficulty arises due to the drifting of the conditional
density function. The untranslated Hermite functions with
limited truncation modes could only resolve the function
well, if it is concentrated in the neighborhood of the
origin. Let us call this neighborhood as a "window".
Unfortunately, the density function will probably drift
out of the current "window". The numerical evidence is
displayed in Fig. |4.5| To efficiently solve this problem,
we for the first time introduce the translating factor to
the Hermite functions and the moving-window technique
for the on-line implementation. The translating factor
helps the moving-window technique to be implemented
more neatly and easily. Essentially speaking, we shift the
windows back and forth according to the "supporf of the
density function, by tuning the translating factor.

This paper is organized as follows. Section II introduces
the generalized Hermite functions and the guidelines of
choosing suitable scaling factor to improve the resolution;

2

HSM to ID FKE and its application to the NLF problems

section III focuses on the analysis of the convergence
rate of HSM to FKE and a numerical verification is
displayed. Section IV is devoted to the application to the
NLF problems. The translating factor and the moving-
window technique are addressed in detail. Numerical
simulations of two NLF problems solved by our algorithm
are illustrated, compared with the particle filter. For the
readers' convenience, we include the detailed formulation
of our algorithm in Appendix A and the proof of Theorem
|2.1| in Appendix B.

II. Generalized Hermite functions

Let us introduce the generalized Hermite functions
and derive some properties inherited from the physical
Hermite polynomials.

Let L^(R) be the Lebesgue space, which equips with
thenorm||-|| = {Jp\-\'^dx)^ and the scalar product {■,■) .
Let Unix) be the physical Hermite polynomials given by
Unix) = (-l)"e^'9Je~'^', n e Z and n > 0. The
three-term recurrence

■Ho = l, Hi{x) = 2x
and Hn+iix) = 2x'Hn{x) — 2nHn^iix)

(2.1)

is more handy in implementations. One of the well-known
and useful facts of Hermite polynomials is that they are
mutually orthogonal with respect to the weight wix) =

_ 2

e ^ . We define our generalized Hermite functions as

h:^^{x)

=1H„(a(x-/?))e-5"'(-«', (2.2)

for n £ Z and n > 0, where a > 0, /3 G R are some
constants, namely the scaling factor and the translating
factor, respectively. It is readily to derive the following
properties for \2.2y .

I) The {H^'^}'^^o forms an orthogonal basis of
L^{R), i.e.

/ m'^{x)H'^\x)dx^

(2.3)

where Snm is the Kronecker function.

2) H"''^{x) is the n eigenfunction of the following
Strum-Liouville problem

+ \nU{x) = 0,

u{x)))
(2.4)

with the corresponding eigenvalue A„ — 2a^n.
3) By convention, H^-^ = 0, for n < 0. Forn G Z
and n > 0, the three-term recurrence holds:

2a{x - l3)H^'\x) =^/2^H^'f^[x) (2.5)
+ .M^HZ'^.ix);

or

2a\x - I3)m-^{x) ^^H:l\{x)

"+l-^n+^l(^)-

4) The derivative of H"'^ [x) is a linear combination

of H"J,{x) and H^'^,{x):

d.HZ'^ix)

= ^\/A^-H'°fi(3;) - Xn+iH"f^{x)

aH^L^x)-

aH^f.ix). (2.6)

5) Property 1) and 4) yield the "orthogonality" of

{dxm-f'(x)}^^o-

dxH"'^ {x)dxH^'^ {x)dx

4a

(An + A„+i), if m :

I = min{n, m}, if \n — m\ = 2;

0, otherwise.

(2.7)

The generalized Hermite functions form a complete or-
thogonal base in L^(K). That is, any function u G L^(R)
can be written in the form

u{x) = ^ u„H"''^{x)

where {un}^^^ are the Fourier-Hermite coefficients,
given by

a

7^

u{x)HZ'^ {x)dx

(2.8)

Let us denote the subspace spanned by the first + 1
generalized Hermite functions by 7?.jv:

7^iv = span{//o"■''(x),••• , ff°'''(a;)} . (2.9)

In the sequel, we follow the convection in the asymp-
totic analysis that a ~ 6 means that there exists some
constants Ci, C2 > such that C\a < b < C2a; a < b
means that there exists some constant C3 > such that
a < Cs,b. Here, C\,C2 and are generic constants
independent of a, j3 and N.

A. Orthogonal projection and approximation

It is readily shown in )2\\ for a > 0, /3 = that the
difference between an arbitrary function and its orthogonal
projection onto 7?. at in some suitable function space could
be precisely estimated in terms of the scaling factor a and
the truncation mode TV. Let us first introduce the function
space ;9(K), for any integer r > 0,

M^a,;3(K) :={iiG L'(R) : \\u\\r,c,.fi <

fe=0 J

(2.10)

where A^ is in l |2.4[ ( and Uk is the Fourier-Hermite
coefficient in jZ8) . We shall denote W''(R) for short.

3

HSM to ID FKE and its application to the NLF problems

if no confusion will arise. Also, the norms are denoted
briefly as j| • jjr - The larger r is, the smaller space W'^(R)
is, and the smoother the functions in VF'^(R) are. The
index r can be viewed as the indicator of the regularity
of the functions.

Let us define the L^— orthogonal projection P\

,a,/3

L^(R) -i> 7^]v, i.e. given v G -L^(R),

(2.11)

The superscript q, /? will be dropped in P^'*^ in the sequel
if no confusion will arise. More precisely,

JV

Pnv{x) := v„H^''^{x),

n=0

where Vn are the Fourier-Hermite coefficients defined in
\2.S\ . And the truncated error j ju— Pjvuj |r, for any integer
r > 0, has been essentially estimated in Theorem 2.3, 1111
for a = 1, /3 = 0, and in Theorem 2.1 , 1211 for arbitrary
Of > and /3 = 0. For arbitrary q > and /3 7^ 0, the
estimate still holds.

Theorem 2.1: For any it G VF' (R) and any integer <
/i < r, we have

PnuIu <

■2N-

(2.12)

where \u\^ := \ \d^u\\ are the seminorms, if A'^ S> 1.
The proof is extremely similar as those in 1111 and 1211 .
Thus, we omit it here and include it in Appendix B for

B. Guidelines of the scaling factor

From Theorem |2.1[ it is known for sure that any
function in W^'^(R) could be approximated well by the
generalized Hermite functions, provided the truncation
mode A'^ is large enough. However, in practice, "suffi-
ciently" large is impossible. To improve the resolution
of Hermite functions with reasonable large we need
the scaling factor a, as pointed out in |4|. Many efforts
have been made along this direction, refer to (3), l4l .
1131 , etc. However, the optimal choice of a (with respect
to the truncation error) is still an open problem. In this
subsection, we give a practical guideline to choose an
appropriate scaling factor for the Gaussian type and super-
Gaussian type functions.

It is well known that, for smooth functions f(x) =
'^'^^Q fnH"'^ {x), the exponential decay of /„ with
respect to n implies that the infinite sum is dominated
by the first A'^ terms, that is.

o

for A'' 3> 1. Thus, the suitable scaling factor is supposed
to get the Fourier-Hermite coefficients decaying as fast
as possible. Once the coefficient approaching the machine
error (say 10~^®), many other factors such as the roundoff
error will come into play. Hence, it is wise to truncate the
sequence here. Therefore, we need some guidelines of

choosing not only the suitable scaling factor a but also
the corresponding truncation mode A'^.

Suppose the function f{x) concentrates in the neighbor-
hood of the origin and behaves asymptotically as e"'''^'
with some p > and k > 2, as |x| — >■ +00. Our
guidelines are motivated by the following observations:

1) The function / decays exponentially fast, as |a;| — >
00. So, /„ « J_^^ f{x)H"'^ {x)dx, provided L is
large enough, due to \2.S) .

2) For the exact Gaussian function e"''^ , p > 0, the
optimal a is naturally to be \/2p with the truncated
mode A'^ = 1. It follows by the fact that with

this choice, e = H^'°{x). thus, e ''^ is
orthogonal to all the rest of -ff"'", n > 0. That
is, (e-P^^)o ^ and (e-?^^)^ = 0, n > 1. This
suggests that the more matching the asymptotical
behavior of / to e~3" ^ , the faster the Fourier-
Hermite coefficients decays, and the smaller trun-
cation mode A' is.
ture method to compute the Fourier-Hermite coef-
ficients by \2.9i) . The truncation mode A' has to be
chosen such that the roots of Hermite polynomial
Hn+i cover the domain [—aL, aL] where the
integral \2.S\ is contributed most from both / and
n = 0,--- ,iV.
We describe our guidelines for the Gaussian type and
the super-Gaussian type functions separately as follows.
Case I. Gaussian type, i.e. f{x) ~ e

l^l — >■ +00.

2 122
1) e-''^ ~ e--" "

, p > 0, as
00, which yields

2) The integrand in i2.8i is approximately e

Using the machine error 10 to decide the domain

of interest L, i.e. e

10"

yielding that

L ^ V8p-i In 10;

3) Determine the truncation mode A' such that the
roots of Hermite polynomial Hn+i covers approx-
imately {—aL, aL), where aL « 4Vln 10.

Case II. Super-Gaussian type, i.e. f{x) ~ e~^^ , as
^1 — >■ +00 for some fe > 2, p > 0.

1) Notice that e's"^""^ > e"'"''', when x > 1. Thus,

we require that e 2"
that aL ^ \/321nl0;

10 , which implies

2) We match e we 2 ^ near x = ±L yields
that a ^ ySpLt"^ Hence, L x (16p~MnlO)^,
a fcpfc(lnlO)2 (=;

3) Determine the truncation mode A' such that the
roots of Hermite polynomial Hn+i cover approxi-
mately {—aL,aL).

To exam the feasibility of our strategy, we explore the
Gaussian type function f{x) — e~^^ cos(^). Accord-
ing to the strategy in Case I, we choose the scaling factor
a ^ ^ 3.1, L ^ y^^^ « 1.9194 and A^ ^ 24.
As shown in Fig. |2.1[ the truncation error with a = 3.1
decay the most fast with respect to the truncation mode A'^
and approach the machine error near the 20* frequency

4

HSM to ID FKE and its application to the NLF problems

f(x)=cos(>c;i 0)6"™ , |5=0

□ □ □

□ □ □ □

o

+ 0=3.1
□ a=1

10 15
truncation mode N

Fig. 2.1. The truncation error v.s. the truncation mode of

f{x) = cos (^)e-5"=^ is plotted, with ^ = 0, a = 4,3.1
and 1, respectively.

mode. Meanwhile, the decay of the truncation error with
a = 4 and a — 1 are much slower. Moreover, the
truncation mode A'^ = 24 is appropriate in the sense that
the next few coefficients start to grow, due to the roundoff
error.

Remark 2.1: 1) This strategy is very useful. However,
it is not the optimal scaling factor a. For example, if
f{x) — e~2^ , then the optimal scaling factor a — 1
and A'^ = 0, instead of A'' = 24 from our guideline.

2) Although the scaling factor helps to resolve the
function concentrated in the neighborhood of the origin,
it helps little if the function is peaking away from the
origin. The numerical evidence could be found in Table
|4.2| This is the exact reason why we need to introduce
the translating factor to the generalized Hermite functions
when applying to the NLF problems, see the discussion
of translating factor in section IV.B.

III. Hermite spectral method to ID forward

KOLMOGOROV equation
The general ID FKE is in the form

ut{x,t) =p{x,t)uxx{x,t) + q{x,t)ux{x,t)

+ r{x,t)u{x,t), for {x,t) eRxR+
u{x, 0) = cro(x).

(3.1)

The well-posedness of ID FKE has been investigated
in (2|. We state its key result here.

Lemma 3.1: (Besala, (2)) Let p{t,x), q{t,x), r{t,x)
(real valued) together with pa;, p^cx, Qx be locally Holder
continuous in D — (to,^i) x K. Assume that

1) p{t, 2;) > A > 0, V (t, x) G D, for some constant
A;

2) r{t,x) < 0,y{t,x) G V-

3) {r-qx+ Pxx){t, x)<0,\/ (t, x) G V.

Then the Cauchy problem with the initial condition
u(to, x) = uo{x) has a fundamental solution r{t, x; s, z)
which satisfies

< r{t,x;s,z) < c(t- s)~5

for some constant c and

/oo
r(t, x; s, z)dx < 1.
- 00

Moreover, if 110(2;) is continuous and bounded, then

/oo
T{t, x; to, z)uo{z)dz
- 00

is a bounded solution of p.l[ l.
Through the transformation

w^x, t)

□

p2 [s,t)ds,t

(3.2)

where

q{x,t) =p~^{x,t)

<li3:,t) - -p ^Px{x,t)

+ T

p ^pt{s,t)d.s

(3.3)

equation jSTTJ can be simplified to the following FKE
with the diffusion rate equals 1 and without the convection
term.

Wt{x, t) — Wxx{x, t) + V{x, t)w{x, t), for K X
'w{x, 0) = ^10(3;),

(3.4)

where

V{x,t)

1

+ o / qt{s,t)ds + r[x,t)

(3.5)

Remark 3.2: From the computational point of view, the
form {i.4\ is superior to the original form l |3.1[ ) in general,
when implementing with the HSM.

(i) If both the potential V{x,t) and the initial data

'w{x, 0) are even functions in x, so is the solution to 1 3.4 1

With the fact that the odd modes of the Fourier-Hermite
coefficients of the even functions are identically zeros, it
requires half amount of computations to resolve the even
functions.

(ii) Even when V{x, t) and 'w{x, 0) are not even, it is
still wise to get rid of the convection term, since this term
will drive the states to left and right, and probably out of
the current "window". Shifting of the windows frequently
by the moving-window technique will definitely affect the
computational efficiency.

A. Formulation and convergence analysis

In this subsection, we shall investigate the convergence
rate of the HSM of solving the FKE. Let us consider
the FKE j3.4^ with some source term F{x,t). The weak

5

HSM to ID FKE and its application to the NLF problems

formulation of HSM is to find un{x, t) G TZn such that

{dtUN{x,t),(p)

+ {V{x,t)uN{x,t),ip) + {F{x,t),^),
. un{x, 0) =Pnuo{x),

(3.6)

for all if £ TZn- The convergence rate is stated below:
Theorem 3.2: Assume

-(i + la;]')^ < V{x,t) < c,

for all {x,t) G K X (0,r), for some 7 > and
some constant C. If uo G VK'^(K) and it is the so-
lution to 1 3.4 1 with source term F(x,t), then for u G
L°°{0,T;W''(R)) n i^(0,r;'tV(K)) with r > 27 and

r 47-21- + 2 . 1

iV>max<^ci 2^-1 max{(Q/3)

Proof of Theorem \3.2\ Denote Un ~ Pnu for

a -T max{(a^)*^,l}

it yields that

\\u - UN fit) < ca-^''-^ ma.yi{(apf'' , l}N^''-\

(3.7)

where c* depends only on T, |!itjjioo(o,T;W"-(E)) and

11^111-2(0, T;VK '-(«))•

Before we prove Theorem |3.2[ we need some estimate
on ||x'^i9J!2j((^j;)||2^ fQj. jijjy integer ri,r2 > 0:

Lemma 3.2: For any function u G W^^^^^ (M.), with
some integer ri , r2 > 0, we have

\x''^ d^'^ u\f < a ^'■^ ^ max{(Q^)^''\

(3.8)

Proof: For any integer ri , r2 > 0,

r-i I I 2

X Oj; lt|| =

y^^u„x''^d^'^Hn''^{x)

^ 00 r2+ri

by ([23} and ( |2.6[ l, where for each n fixed, a„,fe is a
product of 2(ri + r2) factors of or ^A„+j, with
— r2 — ri < j < r2 + ri. Let n* > such that

a /? ~ And notice that A„+_, ~ A^+i for

, T- 1 O Tr\r n _L

n + j > and W^fAx) = for n + j < 0. Hence,

we have

00

Jl — Tl * + 1

-2ri-l o\2ri

<Q max{(a/3) ''^ , 1}| juj 1^^+^

simplicity. By i3.4i with source term F{x,t) and the

definition of Un, we obtain that

= {dt{u - Un),<p) = - {ux,'fix} + {V{x,t)u,ifi)
+ {F{x,t),(f) - {dtUN,'p)

^ {dtUN,'p) = - {Ux,'Px} + {V{x,t)u,if)

+ {F{x,t),^), (3.9)
for all (/p G 72. AT. Combine with l |3.6| l, it yields that

{dt{uN - Un), f) ^ - (dxiuN - u), ^x)

+ {V{x,t){uN — u),ip),

for all </? G 7?.]v. Set qn = un — Un- Choose the function
ip — 2qn, then we have

dtWoN^ = - 2\\dxQNf - 2{dx(UN - u), BxQn)
+ 2{V{x,t)QN, Qn)

+ 2{V{x,t)[UN ~u),qn). (3.10)
By Young's inequality,

\{dx{UN ~ u),dxQN)\

< l\\dx{UN -u)f + \\dxgNf.

(3.11)

The assumption V{x,t) < C for {x,t) G M x (0,T)
yields that

{V{x,t)gN,gN) <C\\eNf, (3.12)
for (x,t) G R X (0,T). Moreover, we have

\{V{x,t){UN-u),QN)\

< ^\\V(UN^n)f + l\\gNf,
(3.13)

by Cauchy-Schwartz's inequality. Substitute l |3.11[ l-( |J. 13[ l
into l |3.10| l, we obtain that

dt\\gN\\' -{C+l)\\gNf

<\\V{Un -u)f + ^\\dx{UN - n)f. (3.14)

Notice that V > — (l+|a;p)^, for some 7 > 0. Essentially
by the estimate in Lemma [^2] we can estimate

\\V{UN-U)f

<||(1 + \x\'')-'{UN-u)f < Wix^-' + 1){Un - u)\\''

00

<a"'*^"^max{(a^)''^,l} ^ ^t'iiii-l

n = JV + l

+ \\Un — uf
<a""-'-' miix{iaP)*\l}N^^-''\\u\\l

I — 2r — Iat — '"ll l|2

+ a N kt r-

(3.15)

for any integer r\,r2 > 0, by l|2.3

The estimate of the second term on the right hand side of
l|3.15^ is followed by Theorem |2.1| Again, by Theorem
|2.1| we obtain

\\dx{UN - u)f = \Un - u\i < a-^'-+'N'-'-\\u\\l

(3.16)

6

HSM to ID FKE and its application to the NLF problems

Substitute ([STTSj, \3.16\ into l |3.14^ , we obtain

dtWQNlf -iC+l)\\gN\\''

provided that

Af>max<^ci ^^-i max {(q/3) 1} ^-^^ ,
a^-T max{(a/3)*^,l}"5^

/'

^0

Therefore, it yields that

\\QN\\\t) <a-^^-Vax{(a/3)^^l}iV2^-''

'0

By the triangular inequality and Theorem |2.1|

||u — Mjv||^(t)

< ||£)Ar||' + llM-[/ivl|'

<a-^->-^N^'^-^ [Ml

+ max{(a/?)'''',l} /* e-<^+i'<*-"'||u||?,(s)ds

< c'a-"^'^ max{(a/?)*^, IjA^^^"",

where c* is a constant depending on | |loo(o,t;W(m))>

11,2(0, T;W"-(E))

and T.

S. Numerical verification of the convergence rate

To verify the convergence rate of HSM, we explore an
ID FKE with some source F{x,t). The exact solution
could be found explicitly as our benchmark. The error
v.s. the truncation mode is plotted.

We consider the ID FKE

ut = u-xx — X u + (sin t + cos t + 3x)e

u{x, 0) = xe

(3.17)

for (a;, t) £Rx [0, T]. It is easy to verify that u{x, t) =
[x + sint)e~ 2^ is the exact solution.

Notice that the initial data, the potential and the source
in l |3.17^ are all concentrated around the origin. So, we
set /3 = 0. For notational convenience, we drop /3 in this
example. As to the suitable scaling factor q, from our
strategy in section II.B, we know that it is better to let
a = 1. However, if we do so, the first two modes will
give us extremely good approximation. Hence, the error
v.s. the truncation mode won't be perceivable. Due to this
consideration, we pick a = 1.4 (a little bit away from 1,
but not too far away so that it won't affect the resolution
too much). The formulation l|3.6^ yields

{dtUN,(p) = — {dxUN,dx(p) — {xUN,X(p)

+ {F{x,t),^), (3.18)

for all ((9 G TZn- Take the test functions tp = H"(x), n =
0, 1, • • • , A'^, in l|3.18|l. The numerical solution ujv £ TZn

20 25 30 35

truncation mode N

Fig. 3.2. The L^.g^-ors of the HSM to FKE pi?) v.s. the
truncation mode N = 5, 15, 25, 35 and 45 is plotted, with a =
1.4, /3 = and the time step dt = 10~^.

can be written in the form

UN

The matrix form of ( |3.18^ follows from l |2.5[ > and \2J\ :

dta{t) = Aa{t) + f{t), (3.19)

where a(t) = (ao(t), ai(t), ■ • ■ , ajv(t))^, /(t) =

^/o(t), fi{t), ■ • ■ , ,/]v(t)^ are column vectors with A'^ +

1 entries, fi{t), i = 0, 1, • • • , A^, are the Fourier-Hermite
coefficients and ^4 is a penta-diagonal (A'' + 1) x (A'^ + 1)
constant matrix, where A = —Ai — A2,

2

fc = min{j,j}, |i-j|=2,

, 0, otherwise,

and

{ V(fc + l)(fc + 2)
2q2

fc = min{j,j}, \i-j\=2,

(2i + 1) . .

* =3,

2q2
I, otherwise.

The errors v.s. the truncation mode A'^ at time T =
0.1 is plotted in Fig. |3.2| The ODE l|3TT9j is numerically
solved by central difference scheme in time with the time
step dt = 10"^. It indeed shows the spectral accuracy of
HSM.

IV. Application to nonlinear filtering

PROBLEMS

Recall the brief description of our algorithm in the
introduction (and more details in Appendix A), the off-

7

HSM to ID FKE and its application to the NLF problems

line computation is to numerically solve the FKE ( |1.4| l
repeatedly on each interval [ri,ri+i]. Equation l |1.4| is in
the form of ( |3.1| l with

r(x,t)

\h^/S + Q{gl +ggxx) - fx

where Q, S, /, g and h are in {OJ.

A. Existence and Uniqueness of the solution to the FKE

In this subsection, we interpret the well-posedness
theorem (Lemma 1 3. 1^ for general ID FKE in section III
in the framework of the NLF problems.

Proposition 4.1 (Existence): Let /, g, h in are
Holder continuous functions m T) := R x {ta,ti). gx,
gxx and fx exist and are also Holder continuous in T>.
Assume further that

1)

> A > 0, for some A > 0;

2) S > 0;

3) -\h^/S - fx + Q{gl + ggxx) < C, for some
constant C,

for (x, t) G T). Then there exists a bounded solution
u{x, t) to l |3.1^ , if the initial condition uo{x) is continuous
and bounded.

Proof: Conditions l)-3) in Lemma [3TT] are directly
translated into conditions l)-3) in this proposition with
C < 0. For C > 0, let v{x,t) = e-'^'*~*»)'u(a;, f), then
V satisfies

Vt (a;, t) =p{x, t)v^^ {x,t)+q{x, t)v^ {x, t)

+ {r{x,t)~C)v{x,t), (4.20)

for {x,t) G T), with the initial condition v{x,to) =
110(2;). The coefficients of l |4.20^ satisfy the conditions in
Lemma [3T| Thus, we apply Lemma [3TT] directly to l |4.20^ .
The existence of the solution to \2>A\ follows immediately.

■

Remark 4.3: In practice, the initial data of the condi-
tional density function is either compactly supported or
exponentially decays as \x\ — >■ +00. So, the assumption
on the initial data in Proposition |4.1| is valid.

For concise of notations, we give the uniqueness for
the equation instead of l |3.1| >. It can be easily trans-
formed into each other, due to the bijective transformation

Proposition 4.2 (Uniqueness): There exists
a unique solution to l |3.4| ( in the class that
{u : lim\x\-foo uux = 0} if V{x,t) is bounded
from above in V.

Proof: Case I: Assume V{x,t) < in 2?. Suppose
there exist two distinct solutions to l |3.4^ , say ui and U2.
Denote rj :— ui — U2, and 77 satisfies

Vt = Vxx + V{x,t)r],

(4.21)

in 2? with the initial condition rj{x,to) = 0. Use the
standard energy estimate, i.e. multiplying ( |4.21^ with rj
and integrating with respect to a; in R:

1,

/

T^xW + / V(x,t)r]''dx < -\\Vx\\ < 0,

TABLE 4.2

Trunction error v.s. the "peaking" po of the
Gaussian function f{x)

PO

erroro

errors

-1

3.3 X 10""

1.1 X IQ-^

8.2 X 10-1''

7.7 X IQ-*'

1

1.6 X 10"^^

1.8 X 10"^

2

1.8 X IQ-^

3.3 X IQ-'-'-^

3

7.7 X IQ-"

8.2 X 10"-^^

4

1.1 X IQ-^

1.6 X IQ-'-'-^

The scaling factor is chosen to be 1 according to the guideline in
section II. B and the truncation error is N = 24. The truncation
en'ors with different translating factor /3 is denoted as error^,
which is defined as ||/ - Yln^o UH^'^^W-

by the
lim

that

integration by parts, and the facts that
> rjrjx = and V{x,t) < in D. This yields

ll'7ll'W<

for t £ {to,ti). With the fact that r){x,tQ) = 0, we
conclude that 77 = in D, i.e. ui =112.

Case II: Assume V{x,t) < C, for some C > 0.
We use the strategy in the proof of Proposition |4.1[ Let

v{x,t) = e~^^ °'u{x,t), then v satisfies i3.4i with the
potential V{x, t) — C < in V. By case I, we conclude
the uniqueness of v, so does u. ■
Remark 4.4: The similar conditions as in Proposition
|4. 1| are derived to guarantee the well-posedness of the
"pathwise-robusf DMZ equation in |22| and to estab-
lish the convergence of our algorithm in II7I . They
essentially require that h has to grow relatively faster
then /. They are not restrictive in the sense that most
of the polynomial sensors are included. For example,
f(x) = fox\ g{x) = go(l + a;^)''' and h{x) = hox'- ,
with 5, Q > 0, /o,<7o and fto are constants, j,k,l G N,
provided I > max { , 2fc — 1 } .

B. Translating factor P and moving-window technique

As we mentioned in the introduction, the untranslated
Hermite functions with suitable scaling factor could re-
solve functions concentrated in the neighborhood of the
origin accurately and effectively. However, the states of
the NLF problems could be driven to left and right during
the on-line experiments. It is not hard to imagine that the
"peaking" area of the density function escapes from the
current "window".

The translating factor /3 is introduced under the cir-
cumstance that the function is peaking far away from the
"window" covered by the current Hermite functions. We
translate the current Hermite functions to the "support" of
the function, by letting the translating factor /3 near the
"peaking" area of the function.

In Table 14.21 we list the truncation error of the Gaus-

sian function f{x) = e 2

with various po =

— 1, ■ ■ ■ ,4 and different translating factors /? = or 3.
According to the guidelines in section II.B, the scaling

8

HSM to ID FKE and its application to the NLF problems

u/=o i^L^, + I3j,l3j + L„). For simplicity and clarity,
let us assume further that

1) The operator (i — ^h'^S~^h) is not explicitly
time-dependent;

2) The time steps are the same, i.e. r^+i — Ti — At.
For the storage of the off-line data, on each interval
{—L^+Pjjfij+Lw), it requires to store (TV + 1)^
floating point numbers. Hence, the total (J + 1) intervals
requires to store (J + 1){N -\- Vf' floating point numbers.
As to the number of the flops in the on-line computations,
if no moving-window technique is adopted during the
experiment, for each time step, it requires 0{{N + 1)^)
flops. The number of the flops to complete the experiment
during [0,T] = u'^-^[Ti,Ti+i] is 0{k{N + lf). Suppose
the number of shifting the windows during [0, T] is P,
then the total number of flops is O {{k + P){N + 1)^).

Remark 4.5: If either assumption 1) or 2) is not satis-
fied, then the real time manner won't be affected. That is,
the number of the flops in the on-line experiment remains
the same. But the off-line data will take more storage
as the trade-off. To be more specific, on each interval
(— L„ + /3j, Pj + Lw), it requires to store fc x (A'' + 1)^
floating point numbers, where k is the total number of
time steps. Therefore, the total storage is k{J+l){N+l)'^
floating point numbers.

C. Numerical simulations

In this subsection, we shall solve two NLF problems
by our algorithm: the almost linear sensor and the cubic
sensor. Since the drift term could always be absorbed
into the potential V{x) by the transformation \3.2\ , for
simplicity, in our examples, we set / = 0. Our algorithm
is compared with the particle filter in both examples.
The particle filter is implemented based on the algorithm
described in 1 1 1. And the systematic resampling is adopted
if the effective sample size drops below 50% of the
total number of particles. As we shall see, our algorithm
surpasses the particle filter in the real time manner.

1) Almost linear filter: We start from the signal obser-
vation model

dxt = dvt

{Pj}j=o, for some J > 0, need to be pT^ared before- I '^V* = + 0.25co8Xt)dt + dwt,

hand, such that [—L,L] C u/=o {—L^ + j3j,l3j where Xt, yt € R, vt, wt are scalar Brownian motion

factor is a = 1 and the truncation mode = 24. As
shown in the table, the further the function is peaking
away from the origin, the larger the error is with untrans-
lated Hermite functions. But with appropriate translating
factor, the function could be resolved very well with the
same scaling factor, for example, error,-? ~ 10"" for
f{x) =e-3{— 3)'.

Indeed this fact motivates the idea of moving-window
technique. The suitable width of the window could be pre-
determined if the trunction error of the density function
v.s. various "peaking" po is investigated beforehand. To be
more precise, suppose we know the asymptotical behavior
of the density function of the NLF problem from the
asymptotical analysis, say ~ e"''^ , with some p > 0,
k > 2. According to the guideline in section II.B, the
suitable scaling factor a and the truncation mode with
/3 = could be chosen. With these parameters, the similar
table as Table |4.2| could be obtained, i.e. the truncation
error (erroro) of the function e~*''^~''°' v.s. various po - If
the error tolenrance is given, then the appropriate width of
the window is obtained according to the table. Take Table
|4.2| as a concrete example. If the asymptotical behavior
of the density function is e~2^ , then the scaling factor
a — 1 and the truncation mode — 24. Suppose we
set the error tolerance to be 10~^, then the suitable width
of the window would be 3 + 3 = 6, from the first two
column of Table [4. 2 1 The window covers the origin would
be [-3,3].

Our algorithm with moving-window technique is il-
lustrated in the flowchart Fig. |2.6| It reads as follows.
Without loss of generality, assume that the expectation
of the initial distribution of the state is near 0. During
the experimental time, say [0,r], the state remains in-
side some bounded interval [— L,L], for some L > 0.
We first cover the neighborhood of by the untrans-
lated Hermite functions {H"'^^}^^q, where a, N can
be chosen according to the guidelines in section II.B.
With the given error tolerance, the suitable width of
the window could be pre-defined, denoted as L„. If
[—L,L] C [— Lu,,Lu,], then no moving- window tech-
nique is needed. Hence, the on-line experiment runs
always within the left half loop in Fig. |2.6| Otherwise,

The off-line data corresponding to different intervals
{—Lw+I3j,l3j+Lw) have to be pre-computed and
stored ahead of time. During the on-line experiment, if the
expectation of the state E[xt] moves accross the boundary
of the current "window" (the condition in the rhombic box
in Fig. |2.6| is satisfied), the current "window" is shifted
to the nearby window where E[xt] falls into. That is, the
right half loop in Fig. |2.6| is performed once.

Let us analyze the computational cost of our algorithm.
Notice that only the storage capacity of the off-line data
and the number of the flops for on-line performance
need to be taken into consideration in our algorithm.
Without loss of generality, let us assume as before
E[a;](0) is near and our state is inside [~L,L] C

processes with E[dVf dvt

1, E[dwt dwt

Suppose the signal at the beginning is somewhere near
the origin.

The corresponding FKE in this case is

Ut =

1 _1 2

^u^. ^x

(1 +

COS a; u

(4.22)

Assume further that the initial distribution of xq is
uq{x) — e 2 . This assumption is not crucial at all. The

non-Gaussian ones, for example ito(a;) — e 2 , vvill give
the similar results as the Gaussian one.

It is easv to se e that the asymptotical behavior of the

solution to 1 4.22 1 is e 2 . With the guidelines in section
II.B, we choose a = 1, /3 = and TV = 25 for the

9

HSM to ID FKE and its application to the NLF problems

Almost linear sensor

14

-2 • • • • • • • • • 1

5 10 15 20 25 30 35 40 45 50
time

Fig. 4.3. Almost linear filter is investigated with our algorithm
and the particle filter with 10 and 50 particles. The total experi-
mental time is T = 50s. And the update time is At = 0.01.

starting interval. We shall run the experiment for the total
time T — 50s. Thus, we expect the density function
probably will move out of the starting interval. Table
|4.2| suggests that the appropriate width of the window
should be 3, if the error tolerance is set to be lO^"".
We shall overlap the adjacent windows a little bit to
prevent frequent shifting of windows. Let us take the
width of the overlaped region to be 0.5. Therefore, as the
preparation for the moving-window technique, we shall
prepare the off-line data for [-19.5,-13.5], [-14,-8],
[-8.5,-2.5], [-3,3], [2.5,8.5], [8,14] and [13.5,19.5].
The correpsonding /3s are -16.5,-11,-5.5,0,5.5,11
and 16.5. The barrier in the rhombic box in the flowchart
Fig |2.6| should be 3 (the width of the "window").

Our algorithm is compared with the particle filter with
10 or 50 particles in Fig. |4.3| for the total experimental
time T = 50s. The time step is At = 0.01s. All three
filters show acceptable experimental results. It is clear
(between time 10 to 30) that the particle filter with 50
particles gives closer estimation to our algorithm than that
with 10 particles. But as to the efficiency, our algorithm
is superior to the particle filter, since the CPU times of
particle filter with 10 and 50 particles are 5.00s and
35.75s respectively, while that of our algorithm is only
2.62s. As to the storage, the size of the binary file to keep
the off-line data is only 35.5kB. During this particular on-
line experiment, the window has been shifted for 13 times,
which can't be perceivable from the figure at all. And it
seems that the moving-window technique doesn't affect
the efficiency of our algorithm.

2) Cubic sensor in the channel: We consider cubic
sensor in the channel xt G [—3,3]:

dxt ~ dvt

3 (4.23)
dyt — Xfdt + dwt,

where xt, yt G K, vt, Wt are scalar Brownian motion
processes with Eldvfdvt] = 1, Elduufdiut] — 1.
Assume the initial state is somewhere near 0.

Cubic sensor

5 10 15 20 25 30 35 40 45 50
time

Fig. 4.4. Cubic sensor in the channel is experimented for T =
50, with the time step At = 0.01s, by both particle filter and
our algorithm.

50

X

Fig. 4.5. The normalized density functions are plotted every
other Is for the cubic sensor in the channel.

The FKE is

1 1 6

Ut = -Uxj: — -X u. (4.24)

Furthermore, we assume the initial distribution is uq [x) =
g-i /4 gjjj^-g {jjg gf^fg inside the channel, we set
our translating factor /3 = and the moving-window
technique won't be used. According to section II.B, we
choose the scaling factor a ^ 23 (iB^) ' ^ 2.4637, and
the truncated mode TV ~ 45.

In Fig. |4.4| we compare our algorithm with the particle
filter with 50 particles for T — 50s. The observation data
come in every 0.01s. Fig. |4.4| reads that both filters work
very well. The result of our algorithm nearly overlaps with
that of the particle filter all the time. However, the CPU
time of our algorithm is 4.90s, while that of particle filter
is 37.17s. With our algorithm, the on-line computational
time for every estimation of the state is around 0.001s,
which is 10 times less than the update time 0.01s. This
indicates that our algorithm is indeed a real-time solver.
The normalized density functions, which is defined as
"^^''l — TT, have been plotted every other Is in Fig

max^gj u(x.t] f J o

El

10

HSM to ID FKE and its application to the NLF problems

V. Conclusions

In this paper, we first investigate the HSM applied to the
ID FKE. It is well-known that the choice of the scaling
factor Q is crucial to the resolution of HSM. We give a
practical guidelines to help choosing the suitable one. The
convergence rate of the HSM has been shown rigorously
and has been verified by a numerical experiment. As an
important application, we solve the NLF problem, by
using the algorithm in fTTl . in the last section, where
solving ID FKE serves as the off-line computation. To
capture the state even if it drifts out of the "window",
translating factor of Hermite functions and the moving-
window technique are introduced. The translating factors
help the switch of the windows back and forth easier,
according to the "support" of the density function of the
state. We analyzed the computational complexity of our
algorithm in detail, with respect to the storage capacity
of off-line data and the number of flops of the on-line
computations. Finally, two online experiments - almost
linear filtering and cubic sensor in the channel - are
reported. The feasibility and efficiency of our algorithm
are verified numerically, which surpasses the particle filter
as a real time solver.

Appendix A
The detailed formulation of our algorithm

Starting from the signal model ( |l.l[ l, the DMZ equation
is derived for the unnormalized density function
a{x, t) of the states Xt conditioned on the observation
history Yt = {j/s : < s < f}. In real applications, one
is more interested in the robust state estimators. Hence,
for given observation path yt, let us make an invertible
exponential transformation

a{x,t) ^ explh^ {x,t)S~'^ {t)yt]p{x,t). (A.l)

The "pathwise-robust" DMZ equation is obtained:

f^ix,t) + ^^{h^s-YytPix,t)

the robust DMZ equation with yt

Ti-i < t < Ti, i — 1,2, - ■ ■ ,k

on the interval

= exp(-/i^S ^yt)

■ [exp {h^ S''^yt)p{x,t)]
, pix,0) ^ ao{x).

2

(A.2)

The exact solution to l |A.2| l, generally speaking, doesn't
have a closed form. Hence, we developed an efficient
algorithm to construct a good approximation in 1 17|.

Let us assume that we know the observation time
sequence = ro < ri < ■ • ■ < Tk = T apriorily.
But the observation data {yr;} at each sampling time r,;,
i = 0, ■ ■ ■ ,k are unknown until the on-line experiment
runs. We call the computation "off-line", if it can be per-
formed without any on-line experimental data (or say pre-
computed); otherwise, it is called "on-line" computations.
One only concerns the computational complexity of the
on-line computations, since this hinges the success of "real
time" application.

Denote the observation time sequence sls Vk ~ {0 =
To < Ti < ■ ■ ■ < Tk = T}. Let Pi be the solution of

2

^{x,t) + ^^{h^S-^y yr,_,p,{x,t)

= exp (^-h^ S~^yT,_-^

■ \eyip(h^ S~^yr^_^^^ p,{x,t)^
pi(a;,0) = o-o(a;),

or

_ Pi{x,Ti-i) = pi_i(a;,Ti_i), for i = 2, 3, •

Define the norm of Vk by \'Pk\ — supi<;<j.(ri
Intuitively, as \'Pk\ — > 0, we have

k

in some sense, for all < t < T, where p{x,t) is the
exact solution of ( |A.2| (. To maintain the real time manner,
our algorithm resorts to the following proposition.

Proposition 1.3: For each r^-i < t < Ti, i =

■ ,k.

(A3)

Ti-l).

,■■ ■ ,k, pi{x, t) satisfies i A.3 i if and only if

,{x,t) = exp \h^ {x,t)S'^{t)yr^_,'^p,{x,t), (A.4)

satisfies the FKE ( fL4| l.

The initial data need to be updated as \i.5\ , followed from
JA31.

With the observation time sequence known {ri}*Li,
we obtain a sequence of two-parameter semigroup
{U{t,Ti-i)}*^^i, for Ti-i < t < Ti, generated by
the family of operators {L — |/i"^5'^^/i}t>o. The off-
line computation in our algorithm is to pre-compute
the solutions of \\A) at time t = r^+i, denoted as
{U{Ti+i,T,)(l)i}fl^, where {<;/>; (a::)}gi (chosen as the
initial data at i = Ti) is a set of complete orthonormal base
in L^(R"). These data should be stored in preparation of
the on-line computations.

The on-line computation in our algorithm is consisted
of two parts at each time step r^-i, i = 1, • • • ,k.

• Project the initial condition Ui{x,Ti-i) G L^(R")
at t = Ti-\ onto the base {(j)i(x)}^-^, i.e.,
Ui{x,Ti-i) = '^^^Ui^i4>i{x). Hence, the solution
to (|1.4|l at t = Ti can be expressed as

(A.5)

Ui{x,Ti) ^U{Ti,Ti-l)Ui{x,Ti-l)

OO

= [U{Ti,Ti-l)(jll{x)

computed off-line.

Update the initial condition of ^IA\ at n with the
new observation y^. Let us specify the observation
updates (the initial condition of ^IA\ ) for each
time step. For < t < ri, the initial condition
is ■ui(x,0) = a-o{x). At time t = ti, when the

11

HSM to ID FKE and its application to the NLF problems

observation is available,

U2{x, Ti)

^ exp [h'^(x, Tl)S~^{Tl)yrj^]p2{x, n)

exp [h''^{x,Ti)S^^{Ti)yri]ui{x, n),

with the fact yo = 0. Here, ui{x,ti) =
[W(Ti,O)0i(a;)], where {uij}fZi is com-
puted in the previous step, and {U{Ti,0)(j)i{x)}'j^i
are prepared by off-line computations. Hence, we
obtain the initial condition U2{x,ti) of for the
next time interval ri < t < T2. Recursively, the
initial condition of l|1.4b for r,;_i < t < is

and

Ui{x,Ti^l) =ex.p[h'^ {x,Ti-l)S ^(Ti_l)

(yr,_i - 2/r,_2)]Wi-l(2:,n-l),

(A.6)

for i — 2,3, ■•■ where Ui_i (a;, Ti_i) =

Ui-2,l pl{Ti^l,Ti^2)<t)l{x)].

The approximation of p[x,t), denoted as p{x,t), is ob-
tained

k

P{^,t) = ^XlTi_i,Ti]it)pt{x,t), (A.7)

i = l

where pi{x,t) is obtained from Ui{x,t) by l|A.4[. And

a{x,t) could be recovered by i A.l

Appendix B
The proof of Theorem I2.1I

Proof of Theorem \2.1\ By induction, we first show
that for p = Q. For any integer r > 0,

y OC

II D l|2 V'^ .2

Q ^ — '

n = JV + l

y OC

yT^ >-r -.r ~2

n = N + l

(B.l)

Suppose for 1 < < r, ( |2.12^ holds for ^ — 1. We need
to show that j2.12| ( is also valid for p. It is clear that

\U — Pnu\^ <\dxU — PnOxU]^-!

+ \PNdxU - dxPNu\^^x. (B.2)

On the one hand, due to the assumption for /i — 1, we
apply j2.12[ l to dxU and replace p and r with p ~ 1 and
r — 1, respectively:

\dxu - PNdxu\^-i <a''-''-^N'^\\dxu\\r-i

<a""'-5iV^||ti||^, (B.3)

where the last inequality holds with the observation that

°° 2

a

: / udxH^"
Jm

20F
2^

(a;)cia;

uHl'J^^{x)dx, by (12.6 1

n A

1tn + l 7, — fin-l-

2 ^ 2

On the other hand, by the virtue of \2.()\

PnQxU — OxPnu

oo

^PN'^UndxH^''^ {x) - UndxH"''^ {x)

n=0 n=0

n=0
^ iV + 1

n =

-2 5Z \/An + lUnH"+i(s)
n=0

1 ^

+ 2

= ^\/Ajv+i |^'Ujv-H'^'+i(3;) + ujv+i//^'''(x)j
This yields that

IPnOxU — dxPNu\l,-l

<Xm+i {ul\H^iA^)\l.i+ul+,\H%-^{x)\

(B.4)

due to the property of seminorms. Moreover, we estimate

ul and \H^-'^{x)\l_-^, for k = N,N + l:

OO

ul! < ill < - Pn-iu\\^ < a~^''N'''\\u\\l,

(B.5)

by ||rT}. Similarly, -u^^+i < Q-^'^iy'llull^. And

<a'^\\H'^''^{x)\\l_i, by Lemma [T2|

— l-i/i — 1 ^- — — 1

=Q! Ajv < a Ajv+i,
since {H'^'^)^ = 5fejv, for fc £ 2,'^. Simi larly,

(B.6)

-H'^+iIm-i ^ «"^A^ ,\. Substitute (|B.5l and (|B.6|( into

B.4^ , we get

IPnOxU - dxPNu\l-l <a~^''~^A*'"''A'i^_|_il|M||r

<^2M-2r-l^M-r||^||2^ (B.7)

by the fact that \n = 2Na^. Combine jB.2| , |B3
(|B.7^, we arrive the conclusion.

and

12

HSM to ID FKE and its application to the NLF problems

References

 M. S. Arulampalam, S. Maskell, N. Gordon and T. Clapp,
A tutorial on particle filters for online nonlinear/non-
gaussian bayesian tracking, IEEE Trans. Signal Process.,
50(2):174-188, 2002.

 P. Besala, On the existence of a fundamental solution for
a parabolic differential equation with unbounded coeffi-
cients, Ann. Polonici Math., 29:403-409, 1975.

 J. Boyd, The rate of convergence of Hermite function
series. Math. Comp., 35:1039-1316, 1980.

 J. Boyd, Asymptotic coefficients of Hermite function se-
ries, J. Comput. Phys., 54:382-410, 1984.

 A. Bensoussan, R. Glowinski and A. Rascanu, Approx-
imation of the Zakai equation by the splitting up method,
SIAM J. Control Optim., 28: 1420-1431, 1990.

 A. Bensoussan, R. Glowinski and A. Rascanu, Approx-
imation of some stochastic differential equations by the
splitting up methods, Appl. Math. Optim., 25:81-106,
1992.

 T. Duncan, Probability density for diffusion processes with
applications to nonlinear filtering theory. Ph. D. thesis,
Stanford University, Stanford, CA, 1967.

 J. Fok, B. Guo and T. Tang, Combined hermite spectral-
finite difference method for the fokker-planck equation.
Math. Comp., 71(240): 1497-1528, 2001.

 D. Funaro and O. Kavian, Approximation of some diffu-
sion evolution equation in unbounded domains by Hermite
function. Math. Comp., 37:597-619, 1991.
 D. Gottlieb and S. Orszag, Numerical analysis of spectral
methods: theory and applications, Soc. In. and Appl.
 B.-Y. Guo, J. Shen and C.-L. Xu, Spectral and pseu-
dospectral approximation using Hermite function: Appli-
cation in Dirac equation. Adv. Comput. Math., 19:35-55,
2003.

 I. Gyongy and N. Kiylov, On the splitting-up method
and stochastic partial differential equation, Ann. Probab.,
31:564-591, 2003.

 W. E. Hopkins, Jr, Nonlinear filtering of nondegnerate
diffusions with unbounded coefficients. Ph. D. disserta-
tion. Dep. Elec. Eng., Univ. Maryland, College Park, Nov.
1982.

 K. Ito, Approximation of the Zakai equation for nonlinear
filtering, SIAM J. Control Optim., 34:620-634, 1996.

 K. Ito and B. Rozovskii, Approximation of the Kushner
equation for nonlinear filtering, SIAM J. control Optim.,
38:893-915, 2000.

 S. Lototsky, R. Mikulevicius and B. L. Rozovskii, Nonlin-
ear filtering revisited: aspectral approach, SIAM J. Control
Optim., 35(2):435-461, 1997.

 X. Luo and S. S.-T. Yau, Complete Real Time Solution of
the General Nonlinear Filtering Problem without Memory,
preprint. arXiv: 1208.0962

 R. Mortensen, Optimal control of continuous time stochas-
tic systems. Ph. D. thesis. University of California, Berke-
ley, CA, 1966.

 N. Nagase, Remarks on nonlinear stochastic partial differ-
ential equations: An application of the splitting-up method,
SIAM J. Control Optim. 33:1716-1730, 1995.

 J. W. Schumer and J. P. HoUoway, V/a.sov simulations
using velocity-scaled Hermite representations, J. Comp.
Phys, 144:626-661, 1998.

 X.-M. Xiang and Z.-Q. Wang, Generalized Hermite spec-
tral method and its applications to problems in unbounded
domains, SIAM J. Numer Anal., 48(4): 123 1-1253, 2010.

 S. T. Yau and S. S.-T. Yau, Real time solution of nonlinear
filtering problem without memory 11, SIAM J. Control
Optim., 47(l):230-243, 2008.

 M. Zakai, On the optimal filtering of diffusion processes,
Z. Wahrsch. Verw. Gebiete, 11:230-243, 1969.

13

HSM to ID FKE and its application to the NLF problems

Start

Compute the Her-
mite coefficient

{Mi,n(Tj)}^=0

the base {H-'^{x)}^^^

Update, i.e. evaluate
Uiix,Ti) at +/3}^^o

Input Q, 0, and
initial data uo

Compute E[,x]{0)
and

Evaluate Ui{x,Ti+i)

Evaluate Ui(x, Tj)

using {Ui,niTi)}^^o

Compute Hermite
coefficient {ui,n{Ti)}'^^o
under the shifted

bases _^

Let/3 =

Synchronized with
the corresponding
part of the off-line
data, i.e. the values of

{U(Tt+un)H^-P(x)}';^^^

Fig. 2.6. The flowchart of our algorithm, where /3' e {0j}j=o-

14 