Year: 2019

Abstract:
Recurrent neural networks (RNN) are powerful time series modeling tools in machine learning. It has been successfully applied in a variety of fields such as natural language processing (Mikolov et al. (2010), Graves et al. (2013), Du et al. (2015)), control (Fei & Lu (2017)) and traffic forecasting (Ma et al. (2015)), etc. In those application scenarios, RNN can be viewed as implicitly modelling a stochastic dynamic system. Another type of popular neural network, deep (feed-forward) neural network has also been successfully applied in different engineering disciplines, whose approximation capability has been well characterized by universal approximation theorem (Hornik et al. (1989), Park & Sandberg (1991), Lu et al. (2017)). However, the underlying approximation capability of RNN has not been fully understood in a quantitative way. In our paper, we consider a stochastic dynamic system with noisy observations and analyze the approximation capability of RNN in synthesizing the optimal state estimator, namely optimal filter. We unify the recurrent neural network into Bayesian filtering framework and show that recurrent neural network is a universal approximator of optimal finite dimensional filters under some mild conditions. That is to say, for any stochastic dynamic systems with noisy sequential observations that satisfy some mild conditions, we show that (informal) ∀ > 0,∃ RNN-based filter, s.t. lim sup k→∞ ∥∥x̂k|k − E[xk|Yk]∥∥ < , where x̂k|k is RNN-based filter’s estimate of state xk at step k conditioned on the observation history and E[xk|Yk] is the conditional mean of xk, known as the optimal estimate of the state in minimum mean square error sense. As an interesting special case, the widely used Kalman filter (KF) can be synthesized by RNN.

text: ∀ > 0,∃ RNN-based filter, s.t. lim sup k→∞ ∥∥x̂k|k − E[xk|Yk]∥∥ < , where x̂k|k is RNN-based filter’s estimate of state xk at step k conditioned on the observation history and E[xk|Yk] is the conditional mean of xk, known as the optimal estimate of the state in minimum mean square error sense. As an interesting special case, the widely used Kalman filter (KF) can be synthesized by RNN.

heading: 1 INTRODUCTION
text: Recurrent neural network (RNN) is a certain type of neural networks characterized by hidden variables that memorize the history of input sequences, and it has been successfully applied and brought amazing results in many different disciplines including computer vision, natural language processing and optimal control, etc. (Mikolov et al. (2010), Graves et al. (2013), Du et al. (2015), Fei & Lu (2017), Ma et al. (2015)). Its huge empirical success in different engineering disciplines is grounded on the expressive power of RNN. However, how to understand the expressive power of RNN in a quantitative way is not fully understood. Even what RNN expresses is not totally clear. Another type of neural network, deep feed forward neural network has been well characterized as a universal function approximator (Hornik et al. (1989), Park & Sandberg (1991)). However, a similar way to characterize the expressive power of RNN is not obvious.
DNN is a mapping from a finite dimensional Euclidean space to another finite dimensional Euclidean space, that is to say it can be regarded as a vector-valued multivariate function. However, RNN is a mapping from a sequence space to another sequence space and the current output depends on both current input and the whole observation history. The input sequence, in principle, can be infinite. The function of RNN is capturing the relationship between input process and output process. For example, in machine translation, the input process (or the observation process) is sentence in one language and the output process (or the state process) is sentence in another language. And in many other RNN’s application scenarios such as traffic forecast and optimal control, the input is a noisy
observation or measurement sequence and the output is an estimate sequence of a certain quantity, e.g., the traffic speed.. We observe that the function of RNN is quite similar to a filter.
In our paper, we propose to characterize the expressive power of RNN in a quantitative way from the perspective of filtering. We consider a discrete filtering system as 1.{
xk = f(xk−1) + g(xk−1)wk−1, yk = h(xk) + vk, (1)
where the state xk at time instant k is an n-dimensional vector, f is an n-dimensional vectorvalued function, g is an n × r matrix-valued function, {wk, k = 0, 1, · · · } is an r-dimensional white Gaussian process and wk ∼ N (0, Q), where Q is the covariance matrix of wk. yk is the m-dimensional observation (measurement) process, h is an m-dimensional vector-valued function, {vk, k = 1, · · · } is an m-dimensional white Gaussian process and vk ∼ N (0, R), where R is the covariance matrix of vk. And we assume that {wk, k = 0, 1, · · · }, {vk, k = 1, · · · } and the initial state x0 are jointly independent. We use Yk to denote the sequence of observations up to time instant k, i.e., Yk := {y1, · · · , yk}. (2) Given the realization of the sequence of observations Yk, the aim of filtering problem is to compute the optimal estimate of xk conditioned on Yk.
Not surprisingly, recurrent neural network has been proposed to do filtering. James Ting Ho Lo showed that Recurrent Multi-layer Perceptron can be used to synthesize optimal filter (Lo (1994)). However, Lo’s approach is based on simply copying and storing the whole observation history in the hidden variables and thus require the time horizon to be finite, which is fundamentally limited. Besides Lo’s work, many efforts have been made to connect RNN and dynamical system. Wilson & Finkel (2009) implemented a neural network based Kalman Filter (KF) but did not provide theoretical analysis on the approximation error. Parlos et al. (2001) proposed an algorithmic approach to do nonlinear filtering using recurrent neural network architecture but did not provide theoretical gurantee. Schäfer & Zimmermann (2006) shows that recurrent neural network is a universal approximator of dynamical system. However, they only consider the deterministic system dynamics and do not analyze the filtering relationship between two stochastic processes. Around the same time our work was under review, a similar paper was also uploaded in arXiv by Lim et al. (2019), The results in the arXiv paper are similar to the June, 2018 graduation thesis of the first author of the current paper in the use of sufficient statistics, and emphasize on the design and implementation of the algorithm. Our current paper has been dramatically improved by introducing the error analysis and asymptotic convergence of the algorithm. And our work is inspired by the finite dimensional filter such as Beneš filter Beneš (1981).
Compared to the existing work, we make the following specific contributions:
• Motivated by the similarity between RNN and filter, we propose to use the ability to approximate optimal filter to characterize the expressive power of RNN. Unlike existing work on expressive power of RNN (Schäfer & Zimmermann (2006)) where only deterministic system is considered, we consider a stochastic dynamic system with noisy observations and analyze the capability of RNN to estimate unknown state.
• We unify RNN-based filter into Bayesian filtering framework. In this framework, the hidden variables of RNN are interpreted as statistics of the observation history. And the evolution of hidden variables are interpreted as the evolution of statistics.
• Based on the Bayesian filtering framework, we derive our main result: Recurrent Neural Networks (RNN) are universal approximators of a large class of optimal finite dimensional filters. That is to say, RNN estimator’s asymptotic estimation error can be as close to minimum mean square error as desired. As an interesting special case, the widely used Kalman Filter can be synthesized by RNN. The consideration of asymptotic error differentiates us from existing work on expressive power of RNN (Schäfer & Zimmermann (2006)).

heading: 2 PRELIMINARY: THE BAYESIAN FRAMEWORK OF FILTERING
text: We first introduce the definition of minimum mean square error estimate.
Definition 1 (Minimum Mean Square Error (MMSE) Estimate (Jazwinski (1970))). Let x̂ be an estimate of random variable x and LMSE := (x− x̂)T (x− x̂). The estimate that minimizes E[LMSE] is called the minimum mean square error estimate. Theorem 1 (Theorem 5.3 in Jazwinski (1970)). Conditioned on the observation history Yk, the minimum mean square error estimate of state xk is the conditional mean E[xk|Yk].
Proof. Proof can be found in our appendix A.1.
The Bayesian filtering consists of recursive prediction and update procedures (Jazwinski (1970)):
Prediction Step p(xk−1|Yk−1) → p(xk|Yk−1): Given the posterior distribution p(xk−1|Yk−1) at instant k − 1, the prior distribution p(xk|Yk−1) of xk satisfies the Chapman-Kolmogorov equation:
p(xk|Yk−1) = ∫ p(xk|xk−1)p(xk−1|Yk−1) dxk−1; (3)
Updating Step p(xk|Yk−1)→ p(xk|Yk): Given prior distribution p(xk|Yk−1), when the observation yk at instant k arrives, the posterior distribution p(xk|Yk) at instant k is given by equation 4,
p(xk|Yk) = p(yk|xk)p(xk|Yk−1)∫ p(yk|xk)p(xk|Yk−1) dxk . (4)
Then we can get the MMSE estimate by simply doing integration: E[xk|Yk] = ∫ p(xk|Yk)xk dxk. (5)
To facilitate subsequent discussion, we make the following definition. Definition 2 (Sufficient Statistic (Beneš (1981))). If the conditional distribution p(xk|Yk) (or p(xk|Yk−1)) can be fully determined by a vector-valued function sk|k (or sk|k−1) of the observation sequence Yk (or Yk−1), then we say sk|k (or sk|k−1) is a sufficient statistic for p(xk|Yk) (or p(xk|Yk−1)).1
Because sufficient statistic sk|k fully determines the posterior distribution p(xk|Yk) and MMSE estimate E(xk|Yk) is a functional of p(xk|Yk), there exists a function γ that maps sk|k to E[xk|Yk],2 i.e.,
E[xk|Yk] = γ(sk|k). (6)

heading: 3 RNN BASED FILTER’S ARCHITECTURE
text: Motivated by the Bayesian framework of filtering, we propose the RNN based filter’s architecture as shown in Fig. 1.
Our RNN based filter’s architecture, Bayesian Filter Net (BFN), consists of three parts: prediction network, update network and estimation network. To mimic the prediction step in Bayesian filtering, we use the prediction network to map the posterior distribution representation vector to a prior distribution representation vector. To mimic the update step in Bayesian filtering, we then use an update network to update the prior state distribution representation vector and the observation to get the posterior state distribution representation vector. Finally, we use an estimation network to map the current posterior state distribution representation vector to current estimation x̂k. We will see in the subsequent discussions, the so-called representation vector or hidden variables indeed can be interpreted as statistics of the underlying conditional distribution.
1 The sufficient statistic sk|k (or sk|k−1) can be any vector-valued function of Yk (or Yk−1) as long as it can fully determine p(xk|Yk) (or p(xk|Yk−1)). Taking linear system as an example, sk|k can be the vector composed of the conditional mean and covariance because these two quantities sufficiently determine the conditional density function of the state in linear filtering system. (See more details in our Appendix A.3.)
2For instance, as for the Kalman Filter in linear system, the sufficient statistic sk|k can be chosen to be the conditional mean and conditional covariance, thus γ function can be the function that copy the conditional mean part of sk|k.

heading: 4 RNN BASED FILTER IS UNIVERSAL
text: We now show that our proposed RNN based filter is universal in that it can approximate a large class of optimal finite dimensional filters to any asymptotic accuracy we desire. We summarize our insight into the diagram 2.
As shown in the diagram 2, we model the conditional prior probability p(xk|Yk−1) ( posterior
probability p(xk|Yk) resp. )
of the state at step k by finite dimensional prior statistics sk|k−1 (posterior statistics sk|k resp.). In finite dimensional filter case (Beneš (1981; 1985); Daum (1987)), there exist finite dimensional statistics that are sufficient, that is to say, the evolution of conditional probability can be fully captured by the evolution of a finite dimensional vector. We denote the evolution function in updating step by ϕ and the evolution function in prediction step by φ. And further, after modelling the the probability distribution as finite dimensional statistics, we use two neural networks to approximate the evolution of them. And to get the final estimate of the state, we use another neural network to approximate the map from the statistics to the optimal estimation.
For the RNN based filter, one can naturally ask:
1. Will the neural network approximation error accumulate and blow up when time goes to infinity and make RNN asymptotically not work at all?
2. How general is this approach? What filter can RNN approximate?
We give the answers in this section.

heading: 4.1 KALMAN FILTER (KF) CAN BE SYNTHESIZED BY RNN
text: When system equation 1 is linear and satisfies Gaussian noise assumption as shown in equation 7, it is well known that the filtering problem can be optimally solved by Kalman Filter (KF) (Kalman (1960)). (See more details in our Appendix A.3.){
xk = Fxk−1 +Gwk−1 yk = Hxk + vk, (7)
where F,G,H are constant matrices with proper dimensions, the initial state x0 is Gaussian, and {wk, k = 0, 1, · · · } and {vk, k = 1, · · · } are two independent white Gaussian sequences that are also independent of the initial state x0.
It can be known that the prior and posterior distributions are Gaussian and fully characterized by the sufficient statistics s composed of mean and covariance matrix. (See more details in our Appendix A.3.) {
sk|k−1 := [mk|k−1, vecT (Pk|k−1)]T ,
sk|k := [mk|k, vecT (Pk|k)]T , (8)
where mk|k−1 and Pk|k−1 are the conditional mean and covariance in step k conditioned on Yk−1 and mk|k and Pk|k are the mean and covariance in step k conditioned on Yk, and vec(◦n1×n2) is the n1n2 × 1 column vector obtained by stacking the columns of the matrix ◦ on top of one another. sk|k (sk|k−1 resp.) is the theoretical statistic that determines the conditional probability distribution of the state xk conditioned on the observation history Yk (Yk−1 resp.) and evolves according to some function ϕ and φ. (See more details in our Appendix A.3.) We also know that there exists some function γ that maps sk|k to MMSE estimate E[xk|Yk]. Thus we have,
sk|k = ϕ(sk|k−1, yk), sk+1|k = φ(sk|k), E[xk|Yk] = γ(sk|k). (9)
We use function ϕ̃ generated by a deep (feedforward) neural network (i.e. update network) to approximate ϕ, and use φ̃ generated by another DNN (i.e. prediction network) to approximate φ. And the numerical statistics computed by RNN are denoted as s̃k|k and s̃k+1|k, i.e.,
s̃k|k = ϕ̃(s̃k|k−1, yk), s̃k+1|k = φ̃(s̃k|k). (10)
And we use function γ̃ generated by the third deep feedforward neural network (i.e. estimation network) to approximate γ in equation 6, i.e.,
x̂k|k = γ̃(s̃k|k). (11)
We consider the probability space (Ω,F ,P) with inner product 〈x, y〉 = E[xT y] and norm ‖x‖ := E1/2[xTx], which is a Hilbert space and denoted as L2(Ω,F ,P). We first state the universal approximation theorem of feedforward neural network before we proceed to show our results. Theorem 2 (Universal Approximation Theorem (Hornik et al. (1989))). For any given compact subset K ⊂ Rn, any given continuous function f defined on K and any given accuracy degree > 0, there exists a function g represented by a single-hidden-layer neural network with non-constant and bounded activation function such that maxx∈K |f(x)− g(x)| < .
Proof. It is a natural corollary of the Thm. 2.1 in Hornik et al. (1989).
Define ek|k := ‖sk|k − s̃k|k‖, which represents the cumulative error caused by the approximation error of φ̃ and ϕ̃. Similarly, we define ek|k−1 := ‖sk|k−1 − s̃k|k−1‖. In the following theorem, we shall give the condition which ensures the cumulative error will not blow up as time k approaches∞. Before we proceed to show our main result, we first establish a key lemma, Lem. 2. We make two assumptions on the system we’ll consider.. Assumption 1. We assume that the linear dynamic system of the state in equation 7 is stable in mean square sense (Samuels (1959)), i.e.,
lim k→∞
‖xk‖ ≤M, (12)
where M is a finite constant.
Assumption 2. The dynamical system equation 7 is uniformly completely observable and uniformly completely controllable.
The definitions of uniformly completely observable and uniformly completely controllable can be found in section 7.5 of Jazwinski (1970). We then state a lemma on the boundedness of conditional covariance. Lemma 1 (Lemma 7.1 in Jazwinski (1970)). If Assumption 2 is satisfied and P0|0 < 03, then Pk|k is uniformly bounded from above for all k ≥ N ,
Pk|k 4
( 1 + αβ
α
) I, k ≥ N, (13)
where N is a positive integer, I is the n× n identity matrix and α, β are positive constants.
Based on Lemma 1, Assumption 1 and Assumption 2, we give the key Lemma 2. Lemma 2. In the discrete linear system equation 7, suppose the Assumption 1 and Assumption 2 are satisfied, then for any given > 0 there exists a compact subsetK ⊂ Rdim(sk|k) such that the statistics computed by KF sk|k−1, sk|k and the statistics s̃k|k−1, s̃k|k computed by RNN based filter with nonconstant and bounded activation function satisfy
∥∥sk|k−11sk|k−1 6∈K∥∥ < , ∥∥sk|k1sk|k 6∈K∥∥ < ,∥∥s̃k|k−11s̃k|k−1 6∈K∥∥ < and ∥∥s̃k|k1s̃k|k 6∈K∥∥ < , where 1A is an indicator function. Proof. Proof can be found in our appendix A.4.
We then derive our main result. Theorem 3. Assume sk|k, k ≥ 0 are the theoretical statistics evolving according to equation 9 and s̃k|k, k ≥ 0 are the real statistics computed by our RNN-based filter evolving according to equation 10. Suppose the Assumption 1 and Assumption 2 are satisfied. Furthermore, we need to assume functions ϕ, φ, γ are Lipschitz, i.e., for any s1, s2,
‖ϕ(s1, y)− ϕ(s2, y)‖ ≤ Cϕ‖s1 − s2‖, ‖φ(s1)− φ(s2)‖ ≤ Cφ‖s1 − s2‖, ‖γ(s1)− γ(s2)‖ ≤ Cγ‖s1 − s2‖,
(14)
where Cϕ and Cφ are Lipschitz constants. If Cϕ and Cφ satisfy |CϕCφ| < 1, then for any > 0, there exists an RNN based filter (with non-constant and bounded activation function) such that
lim sup k→∞ ek|k = lim sup k→∞ ∥∥sk|k − s̃k|k∥∥ < . (15) Furthermore, we have
lim sup k→∞ ∥∥x̂k|k − E[xk|Yk]∥∥ < . (16) Proof. For any δ > 0, we have the following. By Lem. 3, there exists a compact ball K = B(0, r) ⊂ Rdim(sk|k), such that
∥∥sk|k−11sk|k−1 6∈K∥∥ < δ, ∥∥sk|k1sk|k 6∈K∥∥ < δ, ∥∥s̃k|k−11s̃k|k−1 6∈K∥∥ < δ and∥∥s̃k|k1s̃k|k 6∈K∥∥ < δ. By Theorem 2, given any small δϕ, δφ ∈ R+ and δγ ∈ R+, there exist two functions ϕ̃, φ̃ which are represented by the DNN, such that
‖ϕ− ϕ̃‖K∞ ≤ δϕ, ‖φ− φ̃‖K∞ ≤ δφ, ‖γ − γ̃‖K∞ ≤ δγ . (17) where ‖h‖K∞ := maxx∈K |h(x)|. And without loss of generality, we set φ(0) = φ̃(0), ϕ(0) = ϕ̃(0)
and γ(0) = γ̃(0). In the prediction step, based on the evolution equations equation 9 and equation 10, we have ek|k−1 = ‖(sk|k−1 − s̃k|k−1)‖ = ‖φ(sk−1|k−1)− φ̃(s̃k−1|k−1)‖
≤ ‖φ(sk−1|k−1)− φ(s̃k−1|k−1)‖+ ‖φ(s̃k−1|k−1)− φ̃(s̃k−1|k−1)‖ ≤ ‖φ(sk−1|k−1)− φ(s̃k−1|k−1)‖+ ‖(φ(s̃k−1|k−1)− φ̃(s̃k−1|k−1))1s̃k−1|k−1∈K‖
+ ‖(φ(s̃k−1|k−1)− φ̃(s̃k−1|k−1))1s̃k−1|k−1 6∈K‖ ≤ Cφek−1|k−1 + δφ + (Cφ + Cφ̃)δ,
(18)
3Here, X < Y (X 4 Y resp.) if and only if X − Y (Y −X resp.) is positive semi-definite, where X and Y are symmetric matrices.
where the last inequality follows from equation 14 and equation 17 and Cφ̃ is the Lipschitz constant of φ̃. We let δ ′
φ := δφ + (Cφ + Cφ̃)δ.
Similarly, in the updating step, we have ek|k ≤ Cϕek|k−1 + δ ′ ϕ ,where δ ′
ϕ := δϕ + (Cϕ + Cϕ̃)δ. Combining this and equation 18, we obtain
ek|k ≤ Cϕek|k−1 + δ ′ ϕ ≤ (CϕCφ) ek−1|k−1 + ( Cϕδ ′ φ + δ ′ ϕ ) . (19)
Using equation 19 repeatedly, it follows that 4
ek|k ≤ (CϕCφ) k e0|0 + ( Cϕδ ′ φ + δ ′ ϕ ) (CϕCφ)k − 1 CϕCφ − 1 . (20)
Thus lim sup k→+∞
ek|k ≤ ( Cϕ ( δφ + (Cφ + Cφ̃)δ ) + δϕ + (Cϕ + Cϕ̃)δ ) 1 1− CϕCφ as k → ∞
once the condition |CϕCφ| < 1 holds. We choose small enough δφ, δϕ and δ such that( Cϕ ( δφ + (Cφ + Cφ̃)δ ) + δϕ + (Cϕ + Cϕ̃)δ ) 1 1− CϕCφ < . Then we get 15. Now we prove 16.∥∥x̂k|k − E[xk|Yk]∥∥ = ∥∥γ̃(s̃k|k)− γ(sk|k)∥∥ ≤ ∥∥γ(sk|k)− γ(s̃k|k)∥∥+ ∥∥γ(s̃k|k)− γ̃(s̃k|k)∥∥ ≤ ∥∥γ(sk|k)− γ(s̃k|k)∥∥+ ∥∥(γ(s̃k|k)− γ̃(s̃k|k))1s̃k−1|k−1∈K∥∥+ ∥∥(γ(s̃k|k)− γ̃(s̃k|k))1s̃k−1|k−1 6∈K∥∥
?1 ≤Cγek|k + δγ + (Cγ + Cγ̃)δ ?2 ≤ Cγ (CϕCφ)k e0|0 + Cγ ( Cϕδ ′ φ + δ ′ ϕ ) (CϕCφ)k − 1 CϕCφ − 1 + δγ + (Cγ + Cγ̃)δ
(21)
where the inequality ?1 follows equation 14 and equation 17, the inequality ?2 follows equation 59, and Cγ̃ is the Lipschitz constant of γ̃. Thus lim sup
k→+∞ ∥∥x̂k|k − E[xk|Yk]∥∥ ≤ Cγ(Cϕδ′φ + δ ′
ϕ)(−CϕCφ + 1) −1
+ δγ + (Cγ + Cγ̃)δ. Again we can choose small enough δφ, δϕ, δγ and δ such that Cγ(Cϕδ ′ φ + δ ′ ϕ)(−CϕCφ + 1) −1 + δγ + (Cγ +Cγ̃)δ < . Then we obtain the desired 16.
An example satisfying all the assumptions of Thm. 3 can be found in our Appendix A.5, and a general class of systems satisfying assumption |CϕCφ| < 1 can be found in Appendix A.6. We also remark that in our proof, we implicitly require that the Lipschitz constants of φ̃, ϕ̃, γ̃ are uniformly upper bounded. (See more details in our Appendix A.8) Thm. 3 highlights that the optimal filter in linear system with Gaussian noise, Kalman Filter, can be synthesized by RNN. And RNN based filter’s asymptotic error can be as small as wanted under some Lipschitz conditions. That is to say, RNN is an approximator of Kalman filter.

heading: 4.2 RNN BASED FILTER IS A UNIVERSAL APPROXIMATOR OF OPTIMAL FINITE DIMENSIONAL FILTER
text: Thm. 3 shows that Kalman Filter (KF) can be synthesized by RNN. In this section, we’ll try to answer the question "How general is the RNN based filter?" and extend the result into a more general case. We’ll show that any optimal finite dimensional filter can be universally approximated by RNN under some mild conditions. For a general system with noisy observation as shown in equation 1, once the conditional distribution p(xk|Yk) is obtained, the filtering problem is solved. However, we usually need to solve an infinite number of ordinary differential equations (ODE) in order to solve p(xk|Yk). If the distribution p(xk|Yk) admits a finite dimensional sufficient statistics, then we only need to solve a finite number of ODE (Chen (2003)) and we call such filter finite dimensional filter. Finite dimensional filter has been an active research area after the seminal work (Beneš (1981; 1985)) of Beneš. It’s is a large class of filters. Some nontrivial finite dimensional nonlinear filter examples can be found in Daum (1986); Ferrante & Runggaldier (1990); GÜnther (1981); Levine & Marino (1986).
Similarly, we use vector Sk|k to denote the finite dimensional sufficient statistics of the posterior distribution p(xk|Yk) and Sk|k−1 to denote the finite dimensional sufficient statistics of the prior
4See more details in our Appendix A.7.
distribution p(xk|Yk−1). The evolution functions of the statistics are denoted as Φ and Ψ, and the map from Sk|k to conditional mean E(xk|Yk) is denoted as Γ, i.e.,
Sk|k−1 = Φ(Sk−1|k−1), Sk|k = Ψ(Sk|k−1, yk), E(xk|Yk) = Γ(Sk|k). (22)
Similarly, in our proposed neural networks, we use DNN generated function Φ̃ (prediction network) to approximate Φ, use another DNN generated Ψ̃ (update network) to approximate Ψ, and use the third DNN generated function Γ̃ to approximate Γ. And the numerical statistics computed by RNN are denoted as S̃k|k and S̃k|k−1, i.e.,
S̃k|k−1 = Φ̃(S̃k−1|k−1), S̃k|k = Ψ̃(S̃k|k−1, yk), x̂k|k = Γ̃(S̃k|k). (23) We also need the following assumption.
Assumption 3. We assume that for any given > 0 there exists a compact subset K ⊂ Rdim(Sk|k) such that the statistics computed by the optiaml finite dimensional filter Sk|k−1 and Sk|k satisfy∥∥Sk|k−11Sk|k−1 6∈K∥∥ < , and ∥∥Sk|k1Sk|k 6∈K∥∥ < . We can see equation 7 is the special case of system equation 1, and it satisfies the Assumption 3 according to Lemma 2. We further have Lem. 3 and Thm. 4. Lemma 3. In the discrete system equation 1, for any given > 0 there exists a compact subset K ⊂ Rdim(Sk|k) such that the statistics S̃k|k−1, S̃k|k computed by RNN based filter with non-constant and bounded activation function satisfy
∥∥∥S̃k|k−11S̃k|k−1 6∈K∥∥∥ < and ∥∥∥S̃k|k1S̃k|k 6∈K∥∥∥ < , where 1A is indicator function.
Proof. The proof is similar to the step 2 of the proof of Lem. 2.
Theorem 4. Consider a discrete filtering system equation 1 with optimal finite dimensional filter and suppose Sk|k, k ≥ 0 are the theoretical statistics evolving according to equation 22 and S̃k|k, k ≥ 0 are the statistics generated by our RNN based filter and evolving according to equation 23. Suppose the Assumption 3 is satisfied. Furthermore, we need to assume that functions Φ and Ψ are Lipschitz, i.e., for any S1, S2,
‖Ψ(S1, y)−Ψ(S2, y)‖ ≤ CΨ‖S1 − S2‖, ‖Φ(S1)− Φ(S2)‖ ≤ CΦ‖S1 − S2‖, ‖Γ(S1)− Γ(S2)‖ ≤ CΓ‖S1 − S2‖
(24)
where CΨ and CΦ are Lipschitz constants. If CΨ and CΦ satisfy |CΨCΦ| < 1, then for any > 0, there exists a RNN based filter with non-constant and bounded activation function such that
lim sup k→∞ ∥∥∥Sk|k − S̃k|k∥∥∥ < . (25) Furthermore, we have
lim sup k→∞ ∥∥x̂k|k − E[xk|Yk]∥∥ < . (26) Proof. The procedure is similar to the proof of Theorem 3.
Thm. 4 highlights that RNN based filter can not only approximate Kalman filter, but any optimal finite dimensional filter under some Lipschitz conditions. Therefore, RNN’s expressive power is characterized as the universal filtering property.

heading: 5 CONCLUSION
text: In our paper, we try to characterize the expressive power of RNN from the filtering perspective. We unify the recurrent neural network into Bayesian filtering framework and show that recurrent neural network is a universal approximator of optimal finite dimensional filters under some Lipschitz conditions. As an interesting special case, the widely used Kalman filter can be synthesized by RNN. Understanding the expressive power of RNN based filter in more general nonlinear filtering cases (with no finite dimensional sufficient statistics) can be a very interesting future direction.

heading: A APPENDIX

heading: A.1 PROOF OF THEOREM 1
text: Proof. It is known that
E[LMSE ] =E[(xk − x̂k)T (xk − x̂k)] =E [ E [ (xk − x̂k)T (xk − x̂k) ∣∣Yk]] , (27) and let
µk = E[xk|Yk].
Since
E [ (xk − µk)T (µk − x̂k) ∣∣Yk] =(xk − µk)T (µk − µk) = 0,
(28)
and E [ (xk − µk)T (xk − µk) ] (29)
is independent of x̂k, equation 27 becomes E[LMSE ] =E [ E [ (xk − x̂k)T (xk − x̂k) ∣∣Yk]] =E [ E [ (xk − µk + µk − x̂k)T (xk − µk + µk − x̂k)
∣∣Yk]] =E [ E [ (xk − µk)T (xk − µk) + (µk − x̂k)T (µk − x̂k)
+2(xk − µk)T (µk − x̂k) ∣∣Yk]]
=E [ (xk − µk)T (xk − µk) ] + E [ (µk − x̂k)T (µk − x̂k) ] =const + E [ (µk − x̂k)T (µk − x̂k) ] (30)
Therefore E[LMSE ] is clearly minimized by setting x̂k = µk = E[xk|Yk].

heading: A.2 THE EVOLUTION OF DISTRIBUTION IN UPDATING STEP OF BAYESIAN FILTER
text: Let p(xk|Yk) denote the conditional distribution of state xk conditioned on the observation history Yk, then using Bayes’ rule, we have:
p(xk|Yk) Bayes’ rule = p(Yk|xk)p(xk) p(Yk) = p(yk, Yk−1|xk)p(xk) p(yk, Yk−1)
= p(yk|Yk−1, xk)p(Yk−1|xk)p(xk) p(yk|Yk−1)p(Yk−1) Bayes’ rule
= p(yk|Yk−1, xk)p(xk|Yk−1)p(Yk−1)p(xk)
p(yk|Yk−1)p(Yk−1)p(xk)
= p(yk|xk)p(xk|Yk−1)
p(yk|Yk−1)
= p(yk|xk)p(xk|Yk−1)∫ p(yk|xk)p(xk|Yk−1) dxk .
(31)
It can be seen that the posterior distribution p(xk|Yk) is determined by the prior distribution p(xk|Yk−1) and likelihood function p(yk|xk), which depends on the observation model and the known distribution of vk.

heading: A.3 THE EVOLUTION FUNCTION OF STATISTICS FOR KALMAN FILTER
text: The widely used Kalman Filter assumes that the posterior distribution at every instant is Gaussian, and then it can be parameterized by the mean and covariance. We consider the following system:{
xk = Fxk−1 +Gwk−1 yk = Hxk + vk, (32)
where F,G,H are constant matrices with proper dimensions, the initial state x0 is Gaussian, and {wk, k = 0, 1, · · · } and {vk, k = 1, · · · } are two independent white Gaussian sequences that are also independent of the initial state x0 jointly.
Then KF can be viewed as the following recursive relationship:
p(xk−1|Yk−1) = N (xk−1;mk−1|k−1, Pk−1|k−1) p(xk|Yk−1) = N (xk;mk|k−1, Pk|k−1) p(xk|Yk) = N (xk;mk|k, Pk|k),
(33)
where in the prediction step, { mk|k−1 = Fmk−1|k−1
Pk|k−1 = GQG T + FPk−1|k−1F T , (34)
in the updating step,  mk|k = mk|k−1 + Pk|k−1H T (HPk|k−1H T +R)−1 · (yk −Hmk|k−1) Pk|k = Pk|k−1 − Pk|k−1HT (HPk|k−1HT +R)−1
·HPk|k−1,
(35)
and N (x;m,P ) is a Gaussian distribution with mean m, and covariance P . The evolution functions of mean and covariance in the prediction step and updating step can be expressed in the following compact form: {
mk|k−1 = φ1(mk−1|k−1, Pk−1|k−1)
Pk|k−1 = φ2(mk−1|k−1, Pk−1|k−1), (36)
{ mk|k = ϕ1(mk|k−1, Pk|k−1, yk)
Pk|k = ϕ2(mk|k−1, Pk|k−1, yk), (37)
where φ1, φ2, ϕ1 and ϕ2 are defined as{ φ1(m,P ) := Fm
φ2(m,P ) := GQG T + FPFT ,
(38)
and { ϕ1(m,P, y) := m+ PH
T (HPHT +R)−1(y −Hm) ϕ2(m,P, y) := P − PHT (HPHT +R)−1HP,
(39)
Here, mk|k−1 = E[xk|Yk−1], mk|k = E[xk|Yk] (40)
are the conditional means of state conditioned on observation history. Pk|k−1 = E [ (xk −mk|k−1)(xk −mk|k−1)T |Yk−1 ] ,
Pk|k = E [ (xk −mk|k)(xk −mk|k)T |Yk ] (41) are the corresponding conditional covariance. Then we know that the posterior distribution p(xk|Yk) is parameterized by the conditional mean mk|k and covariance Pk|k which are functions of Yk, i.e., p(xk|Yk−1) is parameterized by {mk|k, Pk|k}. And the similar conclusion can be obtained for prior distribution p(xk|Yk−1). Now we vectorize {mk|k, Pk|k} as the statistics
sk|k := [mk|k T , vecT (Pk|k)]T . (42)
Vectorizing the second equation in equation 36 and equation 37, we have[ ϕ1(m,P, y)
vec(ϕ2(m,P, y))
] = [ ϕ1(m, vec−1(vec(P )), y)
vec(ϕ2(m, vec−1(vec(P )), y)) ] , ϕ(m, vec(P ), y) = ϕ(s, y),
(43)
and [ φ1(m,P )
vec(φ2(m,P ))
] , φ(m, vec(P )) = φ(s), (44)
where we omit subscript characters without causing confusion, and vec−1n1×n2 is the inverse operator of vec such that vec−1n1×n2(vec(◦n1×n2)) = ◦n1×n2 . Then we obtain the evolution functions of the statistics: sk|k = ϕ(sk|k−1, yk), sk+1|k = φ(sk|k), (45) where ϕ and φ are defined in equation 43 and equation 44.

heading: A.4 PROOF TO LEM. 2
text: Proof. Step 1: It can be known from the KF that the evolution function for conditional covariance is as follows (see appendix): Pk|k−1 = GQG T + FPk−1|k−1F T Pk|k = Pk|k−1 − Pk|k−1HT (HPk|k−1HT +R)−1
·HPk|k−1, (46)
and the initial value is P0|0. It follows that the conditional covariance P evolves in a deterministic manner and without any randomness according to the evolution equation equation 46, therefore P is independent of the observations. Since
Pk|k = E [( xk −mk|k ) ( xk −mk|k )T ∣∣∣Yk] ,
and Pk|k is independent of Yk, we have Pk|k = E [( xk −mk|k ) ( xk −mk|k )T ] . (47)
According to Lemma 1, we have
Pk|k 4
( 1 + αβ
α
) I, k ≥ N, (48)
where N is a positive integer and α, β are positive constants. It can be easily checked that, there exists a positive constant α0, such that
Pk|k 4 α0I, ∀ k ≥ 0. (49)
According to equation 47, we have
‖xk −mk|k‖2 =E [( xk −mk|k )T ( xk −mk|k )] =tr ( Pk|k ) ,
(50)
where tr(?) denotes the trace of matrix ?. Combining equation 49 and equation 50, we know that
‖xk −mk|k‖ ≤ √ nα0, ∀ k ≥ 0. (51)
Then according to equation 12 and equation 51, we have
‖mk|k‖ ≤ ‖xk −mk|k‖+ ‖xk‖ ≤ √ nα0 +M
≤ C0, ∀ k ≥ 0, (52)
where C0 is a positive constant and can be chosen to be √ nα0 +M .
Now we have the conclusion that mk|k and Pk|k are bounded for all k ≥ 0. Similarly, we know that mk|k−1 and Pk|k−1 are all bounded following the similar procedure as above. According to equation 8, we know that all sufficient statistics sk|k−1, sk|k are bounded.
Step 2: We then prove that ∥∥s̃k|k−11s̃k|k−1 6∈K∥∥ < and ∥∥s̃k|k1s̃k|k 6∈K∥∥ < . Recall that s̃k|k = ϕ̃(s̃k|k−1, yk), s̃k+1|k = φ̃(s̃k|k), (53)
where ϕ̃ and φ̃ are approximating functions generated by neural networks with non-constant and bounded activation functions. Let σ denote the activation functions used by the neural network and assume that |σ(x)| < Cσ,∀x ∈ R. And without loss of generality, we assume ϕ̃(s̃k|k−1, yk) = V1σ(W1p1(s̃k|k−1, yk) + b1) + c1, where p1 is the vector valued function that takes input to the last hidden layer’s output, V1 and W1 are matrices with suitable dimension, b1 is a vector with suitable dimension. Then we’ll have, ∥∥s̃k|k∥∥
= ∥∥ϕ̃(s̃k|k−1, yk)∥∥ ≤ ∥∥V1σ(W1p1(s̃k|k−1, yk) + b1) + c1∥∥ ≤‖V1‖
∥∥σ(W1p1(s̃k|k−1, yk) + b1)∥∥+ ‖c1‖ ≤‖V1‖ √ dim(b1)Cσ + ‖c1‖ ,
(54)
where dim(x) is the dimension of the vector of x. We choose the compact set K to cover the compact ball {x ∈ Rdim(s̃k|k) : ‖x‖E ≤ ‖V1‖ √ dim(b1)Cσ + ‖c1‖} where ‖?‖E denotes Euclidean norm
of ?. Then it can be seen that ∥∥s̃k|k−11s̃k|k−1 6∈K∥∥ < . And ∥∥s̃k|k1s̃k|k 6∈K∥∥ < can be shown in a similar way.

heading: A.5 AN EXAMPLE THAT SATISFIES ALL THE ASSUMPTIONS OF THM. 3
text: Here we verify that the system shown in equation 55 with state’s dimension 1 satisfies all assumptions in Theorem 3. {
xk = (1− α)xk−1 + √ αwk−1 yk = αxk + √ αvk,
(55)
where 0 < α < 1 is a small positive parameter.
• As for Assumption 1, from the state equation of equation 55 we have xk =(1− α)xk−1 + √ αwk−1
=(1− α)2xk−2 + (1− α) √ αwk−2 + √ αwk−1
...
=(1− α)kx0 + k−1∑ i=0 (1− α)k−1−i √ αwi,
(56)
then we have
E[|xk|2] =(1− α)2kE[|x0|2] + k−1∑ i=0 (1− α)2k−2−2iαQ
=(1− α)2kE[|x0|2] + (1− α)2k−2 − (1− α)2
(1− α)2 − 1 αQ.
(57)
It can be easily checked that Assumption 1 is satisfied. • As for Assumption 2, it can be easily checked that system shown in equation 55 satifies
Assumption 2 using the definitions of uniformly completely observable and uniformly completely controllable in section 7.5 of Jazwinski (1970).
A.6 A CLASS OF SYSTEM THAT SATISFIES |CϕCφ| < 1 IN THEOREM 3
Actually this condition |CϕCφ| < 1 can be satisfied by a general class of systems. Let us consider the 1-dimensional linear systems with discretization in continuous observation equation, i.e.,
xtk = axtk−1 + gwtk−1 , ytk = b∆txtk + v ∆ tk , (58)
where a, g, b are constants, ∆t = tk − tk−1 << 1 is the step size of sampling, {wk, k = 1, · · · } is an 1−dimensional white Gaussian process and wk ∼ N (0, Q), and {v∆tk} is 1−dimensional white Gaussian process with distribution N (0, R/∆t). According to Jazwinski (1970), we know that the discrete observation system in equation 58 is the sampled system of continuous observations.
Following equation (34)-(35) in Appendix A.2, we have φ1(m,P ) = am, φ2(m,P ) = gQg + a 2P,
and ϕ1(m,P ) = m+K(y − b∆tm), ϕ2(m,P ) = P −Kb∆tP, where K = Pb(∆t)2((b∆t)2P∆t+R)−1. It can be easily checked that
∇φ = [ a 0 0 a2 ] ,∇ϕ = [ 1−Kb∆t (y −∆tm)∂K∂P
0 1−Kb∆t− b∆tP ∂K∂P
] ,
where ∂K∂P = (b∆t) 2R((b∆t)2P∆t+R)−2. Then we have
Cφ = max{a, a2}, Cϕ = E‖∇ϕ‖2 ≤ E‖∇ϕ‖F = 1 +O(∆t), where ‖ · ‖F is the Frobenius norm. If a < 1, which is a very natural stable condition, then
lim ∆t→0 |CϕCφ| ≤ lim ∆t→0 |(a+O(∆t)) = a < 1.
Apparently, we have |CϕCφ| < 1 with small step size ∆t.

heading: A.7 CALCULATION DETAILS OF INEQUALITY 20
text: ek|k ≤ (CϕCφ) ek−1|k−1 + ( Cϕδ ′ φ + δ ′ ϕ ) ≤ (CϕCφ) ( CϕCφek−2|k−2 + ( Cϕδ ′ φ + δ ′ ϕ )) + ( Cϕδ ′ φ + δ ′ ϕ
) = (CϕCφ) 2 ek−2|k−2 + [CϕCφ + 1] ( Cϕδ ′ φ + δ ′ ϕ
) ...
≤ (CϕCφ)k e0|0 + ( Cϕδ ′ φ + δ ′ ϕ ) k∑ i=0 (CϕCφ) k
= (CϕCφ) k e0|0 + ( Cϕδ ′ φ + δ ′ ϕ ) (CϕCφ)k − 1 CϕCφ − 1 .
(59)
A.8 APPROXIMATE φ AND ϕ USING LIPSCHITZ NEURAL NETWORK
In the proof of Thm. 3, we implicitly require the Lipschitz constants of φ̃ and ϕ̃ are uniformly bounded by another constant C. Note that this assumption won’t affect the approximation capability of Deep Neural Network in approximating Lipschitz continuous function φ and ϕ (Thm. 3 in Anil et al. (2018)). With this requirement, the internal state variables of our proposed RNN based filter is always inside a compact ball independent of the choice of neural
network. And we will have ( Cϕ ( δφ + (Cφ + Cφ̃)δ ) + δϕ + (Cϕ + Cϕ̃)δ ) 1 1− CϕCφ < (Cϕ (δφ + (Cφ + C)δ) + δϕ + (Cϕ + C)δ) 1
1− CϕCφ , thus we can choose small enough δ, δφ and δϕ to make ( Cϕ ( δφ + (Cφ + Cφ̃)δ ) + δϕ + (Cϕ + Cϕ̃)δ ) 1 1− CϕCφ < (Cϕ (δφ + (Cφ + C)δ) + δϕ + (Cϕ + C)δ) 1
1− CϕCφ < .

heading: A.9 EXPERIMENTS

heading: A.9.1 EXPERIMENT SET UP
text: We use PyTorch as the deep learning framework to implement our algorithm and run our simulation on CPU clusters with 5 E5-2630/v3/2.40GHz CPUs, each is equipped with 30GB memory. We use the Mean Squared Error(MSE) as the metric to measure the performance of the filter. For subsequent discussion, we make the following definition. Definition 3 (PMSE). For a set of sampled paths {ω1, ω2, · · · , ωN}, we define the Path Mean Squared Error at time slot k PMSEk to be 1N ∑N i=1 ‖x̂k(ωi)− xk(ωi)‖ 2 2.
Further we define the Temporal Path Mean Squared Error TPMSE. Definition 4 (TPMSE). For a set of sampled paths {ω1, ω2, · · · , ωN} and a horizon of simulation 1 : M , we define the Temporal Path Mean Squared Error (TPMSE) to be 1M ∑M k=1 PMSEk.

heading: A.9.2 TRAINING
text: Alg.1 shows how we train the RNN based filter.

heading: A.9.3 KALMAN FILTER CAN BE SYNTHESIZED BY RNN BASED FILTER EFFICIENTLY
text: We consider a stable linear filtering system as in equation 60.
{ xk = (αAn + In)xk−1 + √ αwk−1
yk = αxk + √ αvk,
(60)
Algorithm 1 Bayesian Filter Net (BFN) training algorithm 1: Setting the sampling path number N , the number of simulated steps M, the optimization horizon
length l and the learning rate λ 2: Sampling x0(ωi) i.i.d from the initial distribution σ0(x) and set y0(ωi) = 0 3: for 0 ≤ m ≤ Ml − 1 do 4: Initialize loss(θ) to be 0 5: for ml ≤ k < (m+ 1)l do 6: for 1 ≤ i ≤ N do 7: xk+1(ωi)←− f(xk(ωi)) + g(xk(ωi))wk−1(ωi) 8: yk+1(ωi)←− h(xk+1(ωi)) + vk(ωi) 9: end for
10: for 1 ≤ i ≤ N do 11: Feed the sampled yk+1(ωi) to the Bayesian-Filter-Net(BFN) and get the corresponding output x̂k+1(ωi) 12: Add the estimation squared error to the loss function loss(θ) ←− loss(θ) + ‖x̂k+1(ωi)− xk+1(ωi)‖22 13: end for 14: end for 15: Doing back propagation and update the network parameters θ ←− θ − λ∇θloss(θ) 16: end for
where α < 1 is a small positive parameter controlling the changing rate of the system state, In is a n × n identity matrix and An is a matrix with elements aij satisfying aij = 1, if i = j + 1, − ( n i−1 ) , if j = n and 0, otherwise.
From Fig. 3a, we can see that as the number of training epochs increase, RNN-based filter’s MSE decrease and finally converge. And from Fig. 3b as the dimension of the hidden variable increases, the converged MSE decreases. The results in Fig. 3 highlight that just with the number of 90 hidden variables, RNN based filter can achieve the performance of Kalman filter. Considering that there are 65 independent statistics in this case, this is a very efficient approximation.

heading: A.9.4 RNN BASED FILTER CAN STILL WORK FOR INFINITE DIMENSIONAL FILTER
text: For finite dimensional filter, our analysis has shown that RNN based filter is a universal approximator of optimal filter. In this subsection, we show that in infinite dimensional filter case (where finite dimensional statistics are not admitted to fully determine the distribution), our RNN based method achieves even better performance to state-of-the-art spectral method Luo & Yau (2013). We consider a special and well-known infinite dimensional case.{
xk = xk−1 + √ αwk−1 yk = αx 3 k + √ αvk
(61)
where α is a positive parameter controlling the change rate of the system. We set the α to be 0.01, the total sampling steps to be 5000, the number of sampled paths to be 1000, the number of training epochs to be 100. We only use two fully connected layers for the prediction network, update network and estimation network. Fig. 4a shows the convergence of PMSE of our Deep-Filter-Net with respect to the time. And fig. 4b shows that our proposed Deep-Filter-Net output x̂k tracks xk well. We also
compare our algorithm’s performance with existing state-of-the-art spectral method. Our proposed DFN’s TPMSE is 0.18, which is much smaller than spectral method’s TPMSE 0.40. One potential reason for its RNN’s applicability in infinite dimensional filter case is that RNN implicitly learns the most relevant finite dimensional statistics.

Review summaries:
This paper shows that RNN (of infinite horizon) can be universal approximators  for any stochastic dynamics system. It does that by using the Bayesian framework of filtering, which shows that if sufficient statistics given the observations so far can be computed at each stage of filtering, then the expected hidden states at each stage can be estimated as well. The paper then follows this framework by constructing deep neural network mapping from “observations so-far” to “sufficient statistics” and by universal approximation theorem this is possible. Then as long as after each stage, the image of such a mapping is bounded into a certain region, and if the mapping is contractive (|C_\psi C_\phi| < 1), this can be applied to arbitrarily long sequence (and thus infinite horizon). 

The paper attempts to establish the asymptotic accuracy of the "RNN" (but not the RNN models that are well-known in the literature - see the below comments) as a universal functional approximator. It considers a general state space model and uses feedforward neural nets (FNN) to learn the filtering and forecast distributions. Based on the well-known universal approximation property of FNNs, the paper shows that their RNN-based filter can approximate arbitrarily well the optimal filter.