情報論的学習理論テクニカルレポート 2009
Technical Report on Information-Based Induction Sciences 2009 (IBIS2009)

1 次元線形動的システムの特異性とベイズ汎化誤差への影響
Singularities of One-dimensional Linear Dynamical Systems and its
Eﬀect on the Bayesian Generalization Error
内藤 卓人∗
Takuto Naito

山崎 啓介†
Keisuke Yamazaki

Abstract: Linear dynamical systems are widely used in such ﬁelds as system control andtime-dependent data analysis. Such a system can be regarded as a statistical parametricmodel, where the coeﬃcients of the state space equations are unknown and given as parameters. The properties of parameter learning have not yet been established, in spite of a wide
range of applications. Therefore, this paper investigates the system from the viewpoint oflearning theory. It is revealed that the system has singularities in the parameter space.
The generalization error measured by the prediction accuracy for unseen data sequences isreduced, due to the presence of these singularities.
Keywords: Kalman Filter, Bayesian Learning, Time-Series Data Analysis

1Introduction
Linear dynamical systems are widely used for mod-eling practical complex systems with hidden variables
such as object tracking in image processing [4], andposition detection in car navigation systems [6]. Thesystem is described via state space equations containing both observable and hidden variables. The Kalmanﬁlter [5] is an algorithm to estimate the hidden variables from coeﬃcients given preliminarily .
It is important to be able to estimate coeﬃcientsusing the observable data when the coeﬃcients are unknown. The system is regarded as a parametric learning model, in which the coeﬃcients correspond to parameters. As seen in Section 2, the system is expressed
as a generative probability model of the data becausethe process and observation noises are taken into account.
∗ 東 京 工 業 大 学, 知 能 シ ス テ ム 科 学 専 攻, 226-8503 神 奈 川
県 横 浜 市 緑 区 長 津 田 町 4259 R2-5, tel. 045-924-5018, e-mailnaitaku@cs.pi.titech.ac.jp,
Department of Computational Intelligence and Systems Science,
Tokyo Institute of Technology, R2-5, 4259 Nagatsuta, MidoriKu, Yokohama, Kanagawa† 東京工業大学, 精密工学研究所, 226-8503 神奈川県横浜市緑区長津田町 4259 R2-5, tel. 045-924-5018, e-mail k-yam@pi.titech.ac.jp,
Precision and Intelligence Laboratory, Tokyo Institute of Technology, R2-5, 4259 Nagatsuta, Midori-Ku, Yokohama, Kanagawa

Parametric models generally fall into two types, regular and singular. If the relation between the parameter and the expressed probability function is one-toone, a model is referred to as regular. Otherwise, it issingular. Therefore, a singular model has a set of parameters indicating the same function, in which there
are singularities. Because of the singularities, conventional analysis is not applicable; model selection criteria for regular models such as AIC [1] and BIC [7]are inappropriate. An algebraic geometrical method
has been developed for Bayesian learning to reveal theasymptotic generalization error and the marginal likelihood for several singular models [8]. According to its
application to several models, the presence of singularities results in unique properties of the learning process[3, 9].
In spite of a wide range of applications, properties ofa linear dynamical system are still unknown in termsof a learning model. Therefore, the present paper investigates such a system both theoretically and experimentally. We conﬁrm that the system is a singularmodel and analyze the Bayesian generalization errorbased on the algebraic geometrical method. Here, theerror is deﬁned as the prediction accuracy for unseen
time-sequence data. This prediction is diﬀerent fromthat of the conventional Kalman situation in which theprimary concern is the set of hidden variables rather

where I is a unit matrix and Kt is called the Kalman
gain. Firstly, the current state zt is estimated as ẑt|t−1

than the observable sequences. Nevertheless, our analysis can also provide an insight into hidden variable
estimation.

from the estimated state of the previous time t − 1(Eq.3). Then, a more reﬁned value for ẑt|t is calculatedon the basis of ẑt|t−1 after an observation xt is provided

The remainder of the paper is organized as follows.
Section 2 formulates the system. Section 3 introduces(Eq.6).
From the viewpoint of machine learning, a linear dy-Bayesian learning and summarizes the algebraic geometrical method. Section 4 contains our main contributions, deriving a theoretical upper bound of thenamical system can be regarded as a learning model
whose parameters are A, C, D and z1 . The variable z1
indicates the initial state. Let X = (x1 , x2 , . . . , xT ) ∈generalization error and showing experimental results
for the error. Section 5 contains a discussion and ourRp×T be the vector of observations. The probabilityp(X|w), where the parameters w = (A, C, D, z1 ), canconclusions.

be calculated as follows:T∏

Linear Dynamical Systemsp(X|w) = p(x1 |w)

Linear dynamical systems can be described by state
space models with hidden state variables:

Using the hidden state zt ,2p(xt |x1 , . . . , xt−1 , w).

(8)p(xt |x1 , . . . , xt−1 , w)∫= p(xt |zt , w)p(zt |x1 , . . . , xt−1 , w)dzt .

(9)

t=2

zt+1 = Azt + Dwt ,(1)

xt = Czt + vt ,(2)

Let N (·|µ, Σ) be a multivariate normal distribution

where zt ∈ Rq is the hidden state vector at time t,
xt ∈ Rp is an output vector, wt ∈ Rq and vt ∈ Rp
are process and observation noises, respectively. Thesewith mean µ and covariance matrix Σ. By the definition of a linear dynamical system (Eq.2) and thenoises are assumed follow a standard normal distribution. A ∈ Rq×q is the state matrix , C ∈ Rp×q is thep(xt |zt , w) = N (xt |Czt , I),

derivation of the Kalman ﬁlter,p(zt |x1 , . . . , xt−1 , w) = N (zt |ẑt|t−1 , Pt|t−1 ).

(10)
(11)

output matrix and the elements of D ∈ Rq×q are thecoeﬃcients of the process noise.
The Kalman ﬁlter is known as an eﬃcient recursive

Therefore, p(xt |x1 , . . . , xt−1 , w) is also a normal distribution described byﬁlter that estimates hidden states from a series of outputs. In what follows, the notations ẑn|m and Pn|mp(xt |x1 , . . . , xt−1 , w) = N (xt |C ẑt|t−1 , I + CPt|t−1 C > ).
(12)

represent the estimates of z at time n and its error covariance matrix, respectively, when observations fromt = 1 to t = m are given. The Kalman ﬁlter has twophases: Predict and Update. The algorithms are described as follows:Predictẑt|t−1 = Aẑt−1|t−1(3)>Pt|t−1 = APt−1|t−1 A + DD>

Eq.8 can be expressed asp(X|w) =T∏N (xt |C ẑt|t−1 , I + CPt|t−1 C > ).

(13)

t=1

where we deﬁne ẑ1|0 = z1 and P1|0 = 0.
Let X n = (X1 , X2 , . . . , Xn ) be a set of i.i.d. training
samples. Each Xi is a time sequence deﬁned by Xi =(xi1 , xi2 , . . . , xiT ). The likelihood of the parameter w =(A, C, D, z1 ) can be calculated as(4)L(w) =Updaten∏p(Xi |w)i=1()−1
Kt = Pt|t−1 C I + CPt|t−1 C >()ẑt|t = ẑt|t−1 + Kt xt − C ẑt|t−1(5)

Pt|t = (I − Kt C) Pt|t−1

(7)

>

(6)

=n ∏T
∏

ii
N (xit |C ẑt|t−1, I + CPt|t−1C >)(14)

i=1 t=1ii
where ẑt|t−1
and Pt|t−1
are evaluated using the Kalmanﬁlter.

3

Bayesian Learning and the Gen- m1 , m2 , . . . be the respective orders. The asymptoticform of the generalization error is expressed aseralization Error()
G(n) =This section describes Bayesian learning for time series data and the theoretical analysis of the generalization error.
Let X n = (X1 , X2 , . . . , Xn ) be a set of training samples taken independently and identically from the true
distribution q(X), where n is the number of training
samples. Each Xi (i = 1, . . . , n) is a sequence whose
length is T , i.e. Xi = (xi1 , . . . , xit , . . . , xiT ). Note thatthe sequence data X n are taken as i.i.d. whereas eacha pole z = −λ and its order m have been calculated,
an upper bound is derived as
( 1 )
m−1
λ+oG(n) ≤ −.
(21)
n n log nn log n4Analysis of the Generalization
Error

This section analyzes the Bayesian generalization error for linear dynamical systems. In order to investi-a posteriori probability distribution is deﬁned byn∏
1ϕ(w)p(Xi |w)
Z(X n )
i=1(20)for n → ∞. In many cases, it is not straightforward toﬁnd the largest pole −λ1 and its order m1 [3]. Whensequence Xi is not. Let p(X|w) be a learning model,and ϕ(w) be an a priori probability distribution. The

p(w|X n ) =m1 − 11λ1−+onn log nn log n(15)where Z(X n ) is a normalizing constant. The Bayesianpredictive distribution is deﬁned by∫p(X|X n ) = p(X|w)p(w|X n )dw.
(16)
The Bayesian generalization error G(n) is deﬁned by[∫]q(X)G(n) = EX nq(X) log
dX,(17)
p(X|X n )gate the eﬀect of redundant hidden states, we study
an essential case, in which the learning model has ahidden variable and the true model generates i.i.d. sequences. This is the simplest setting for singularities toexist in the parameter space because the i.i.d. model
can be regarded as a model with no hidden states. For
simplicity, we assume that the output vector is one
dimensional, where zt , xt , A, C, and D are all scalar.
Moreover, we assume that the ﬁrst hidden state is ﬁxed
as z1 = 0. Formally, the learning model is deﬁned aswhich is the average Kullback information from thetrue distribution to the predictive distribution.
The remainder of this section summarizes the algebraic geometrical method for deriving the asymptotic
form of the error [8]. Let H(w) be the Kullback information from the true distribution q(X) to the learnerp(X|w),
∫q(X)
dX.
(18)
H(w) = q(X) log
p(X|w)
The function ζ(z) of one complex variable z, deﬁnedby∫
ζ(z) = H(w)z ϕ(w)dw,(19)is referred to as the zeta function. It is known that thiszeta function is holomorphic in the region Re(z) > 0,and can be analytically continued to the meromorphic function on the entire complex plane. Then thepoles are all real, negative and rational numbers. Let
0 > −λ1 > −λ2 > . . . be a sequence of poles, andzt+1 = azt + dwt ,(22)xt = czt + vt ,(23)

where zt , xt ∈ R1 and wt and vt are distributed from
N (·|0, 1). The parameter is expressed as w = (a, c, d).
The true model is a one-dimensional normal distribution N (xt |0, 1) for all t, i.e. xt = vt . Following Eq. 13,the true model is given byq(X) =T∏

N (xt |0, 1).

(24)

t=1

4.1

Theoretical analysis

Based on the algebraic geometrical method, the error
has the following bound:
Theorem 4.1 When the true model and a learningmodel are deﬁned by Eq.24 and Eqs 22-23, respectively,
the Bayesian generalization error is bounded above as
follows:( 1 )11−+o,(25)
G(n) ≤2n n log nn log nwhere z1 = 0 and the training sample size n is suﬃciently large.
Sketch of Proof: Because the parameter set {c = 0}
attains p(X|w) = q(X), there is a function fc (w) suchthat H(w) = c2 fc (w). The set {d = 0} ensures thesame property for H(w). Thus, there is a polynomial f (w) such that H(w) = c2 d2 f (w). We can ﬁnd alimited support W of the parameter space, such thatH(w) ≤ Cc2 d2 . Here C is a positive constant. Considering the following zeta function∫{Cc2 d2 }z dw,
ζ1 (z) =(26)W

the pole z = −µ is a lower bound of z = −λ1 [8].
We can ﬁnd a pole µ = 1/2 and its order m = 2.
Combining with Eq. 21, we derive the following leading
terms for the bound,11−,
2n n log n(27)which completes the proof.
End of Proof

Theorem 4.1 as follows.
Corollary 4.1 Under the same setting as Theorem
4.1, the error has an upper bound(1)
1G(n) ≤+o.
(28)
2nn

Experimental results

We experimentally evaluate whether the bound isvalid when ﬁnite training data are given. Samplingfrom the a posteriori distribution, the predictive distribution is given by1 ∑p(X|wj ),M j=1M(29)where (w1 , . . . , wM ) are sampled from p(w|X n ). We
use the Markov chain Monte Carlo (MCMC) method
for the sampling technique [2]. The generalization error is approximated byG(n) ' EX nN[1 ∑N i=1

is a normal distribution for a, c and d.
Figure 1-(a) describes an example of sampling fromthe a posteriori distribution in the parameter space(a, c, d). The vertical and horizontal planes indicate{c = 0} and {d = 0}, respectively. The points arelocated around the subspace {c = 0} ∪ {d = 0}, forwhich the parameters express the true model.
Figure 1-(b) summarizes the error values corresponding to n = 250, 500, 750 and 1, 000. The horizontal andvertical axes describe the number of training data sequences and the error value, respectively. The heavyline depicts experimental values for G(n). The dottedline is the upper bound of Theorem 4.1. The upperbound is valid as seen in the graph.

logq(Xi ) ].
p(Xi |X n )Discussions and Conclusions
First, let us discuss the upper bound of the gener-alization error. In the regular case, the error has thefollowing asymptotic form,G(n) =( 1 )dim w+o,
2nn log n(31)which means that λ1 = dim w/2 and m1 = 1. Notethat even a singular model has this asymptotic form
if the true and learning models have the same dimension of the hidden state vector. The asymptotic form
indicates that the cost to ﬁt all parameters determines

We omit the proof for lack of space.

p(X|X n ) 'sequences is N = 1, 000. The number in the MCMC
sample is M = 500. We obtain the expectation EX n [·]over 100 sets of training data. The a priori distribution

5If the initial state is unknown and is regarded asa parameter such as w = (a, c, d, z1 ), we can extend

4.2

The experimental settings are as follows. The length of
the time sequence is T = 10. The number of test data(30)

the error as the dimension dim w appears. Comparing
Theorem 4.1 with the regular case, we can derive theresult that the error is much smaller, i.e.
( 1 )11−+o2n n log nn log n( 1 )3<+o,2nn log nG(n) ≤(32)

which conﬁrms that the ﬁtting cost for redundant parameters is not strongly reﬂected in the error.
Thus far, we have focused on prediction of the unseen observable data sequence X. Next, we consider
estimation of the hidden states zt . According to thea posteriori distribution, there are two regions for theoptimal parameters; one is around c = 0 and the other−33

1.5x 10

Result G(n)1/(2n) − 1/(n log n)211d 0−10.5−2−3−2a0

2−3−2−1c01230250(a)500750

1000(b)Figure 1: An example of the a posteriori distribution and the generalization error.
is around d = 0. They imply completely diﬀerent behaviors of the hidden state. The former, c = 0, indicates that a and d can take any value, by whichq(X) = p(X|w). Thus, there are no constraints on
the movement of the hidden state. By taking into account z1 = 0, the latter, d = 0, contrarily implies that
there is no movement because zt = 0 for all times t. If
several hidden variables in the true model stop moving due to disorder in a practical situation, the desired
estimation is d = 0. However, c = 0 can also be an
estimated result; these variables move on the basis ofarbitrarily-estimated a and d. This adverse estimation
can occur along any dimension of the hidden state vector. Therefore, detection of hidden variable size is an
essential problem to solve.

Vol. 19, pp. 716–723, 1974.
[2] Christophe Andrieu, Nando de Freitas, Arnaud
Doucet, and Michael I. Jordan. An introduction
to MCMC for machine learning. Machine Learning, Vol. 50, No. 1-2, pp. 5–43, 2003.
[3] Miki Aoyagi and Sumio Watanabe. Stochastic complexities of reduced rank regression in Bayesian estimation. Neural Networks, Vol. 18, pp. 924–933,2005.
[4] N. Funk. A study of the Kalman ﬁlter applied tovisual tracking. Technical Report Project for CMPUT 652, University of Alberta, 2003.
[5] R. E. Kalman. A new approach to linear ﬁlter-Finally, we state our conclusions. The present paper
establishes that linear dynamical systems are singular models. The singularities ensure that the uppering and prediction problems. J. Basic Engineering,Vol. 82, pp. 35–45, 1960.

bound of the Bayesian generalization error is small.
The experimental results indicate that the bound is[6] D. Obradovic, H. Lenz, and M. Schupfner. Sensor
fusion in siemens car navigation system. Proc. ofMLSP 2004, pp. 655–664, 2004.

valid. Moreover, the a posteriori distribution implies
that estimation of hidden states cannot be appropriate
if there are redundant hidden variables.

Acknowledgment
This research was partially supported by the Ministry of Education, Science, Sports and Culture in Japan,
Grant-in-Aid for Scientiﬁc Research 18079007.

References[1] H. Akaike. A new look at the statistical model
identiﬁcation. IEEE Trans. on Automatic Control,[7] G. Schwarz. Estimating the dimension of a model.
Annals of Statistics, Vol. 6 (2), pp. 461–464, 1978.
[8] Sumio Watanabe. Algebraic analysis for nonidentiﬁable learning machines. Neural Computation,
Vol. 13, No. 4, pp. 899–933, 2001.
[9] Keisuke Yamazaki and Sumio Watanabe. Algebraic geometry and stochastic complexity of hidden
Markov models. Neurocomputing, Vol. 69, No. 1-3,pp. 62–84, dec 2005.