Transcript 発表資料

1
Kernel Subspace Method by Stochastic Realization
for Learning Nonlinear Dynamical Systems
（Neural Information Processing Systems 2006）
東京大学大学院航空宇宙工学専攻 東京大学先端科学技術研究センター
河原 吉伸
矢入 健久，町田 和雄
動的システム学習
• 時系列入出力データから，その生成モデルを推定する
– モデルは対象システムの制御や動的特性の解析などに利用
入出力データ
{u(t), y(t), t  0, 1,}
動的システム学習
動的モデル
状態空間モデル
ARMAモデル
x(t  1)  Ax(t )  B u(t )  v

 y(t )  C x(t )  Du(t )  w
y(t )  
y(t  i)  
u(t  i)  v
i1
i0
などなど…
応用 : 対象システムの制御，動的特性の解析，…
2
状態空間モデル
線形の場合
3
非線形の場合
x(t  1)  Ax(t )  B u(t )  v x(t  1)  f ( x(t ), u(t ))  v


 y(t )  C x(t )  Du(t )  w  y(t )  g ( x(t ),u(t ))  w
• 直接観測されない状態ベクトルを用いたモデル表現．
– 過去の全ての観測に関する十分統計量となっているため，ARMAモデルな
どの入出力のみに関するモデルに比べ，正確に動的システムを表現可能
• その他，物理科学などで得られた解析モデルとのアナロジーが得
られやすいなどの特徴があるが，比較的推定は困難．
状態空間学習の分類
入出力データ
{u(t), y(t), t  0, 1,}
モデル・データ間の距離の最小化
部分空間上での幾何学的演算
予測誤差 ⇒ 予測誤差法
尤 度 ⇒ EMアルゴリズム
直交分解，CCAによる確率実現
⇒ 部分空間法
システム行列
A, B, C, D, K
カルマンフィルター
状態ベクトル
4
x(t ), t  0, 1,
○ 高精度（但し初期値依存）
× 局所解の問題，要反復計算
状態ベクトル
x(t ), t  0, 1,
最小二乗法
システム行列
A, B, C, D, K
○ 高速で数値的に安定，一意解
× 推定精度はやや劣る
Existing Subspace Methods for Learning
Nonlinear Dynamical Systems
• A nonlinear algorithm is essential for learning complex systems
which cannot be expressed sufficiently with linear models.
• Existing nonlinear subspace methods
– Nonlinear approaches with neural networks
・ Based on embedding and regression with neural nets [VVS00]
– Nonlinear approaches with reproducing kernels
・ Method for Hammerstein systems with LS-SVM [GPSM05]
・ Method using Kernel CCA on retarded coordinate vectors [VSB+04]
=> ・ Operate directly toward input-output data
・ Need additional nonlinear regression frameworks
– Other approaches
・ Using a conditional expectation algorithm for nonlinear CCA [LB92]
5
6
再生核ヒルベルト空間（RKHS）
…
RKHS  k
正定値カーネル: k : Ω Ω  R

1
:
1
k
(

,
x
)


(
x  Ω)
－
－ 対称性 k ( x, y)  k ( y, x)
－ 再生性 f , ( x)  f ( x)
－ 正定値性 in, j 1ci c j k ( xi , x j )  0
正定値カーネルにより特徴空間内の
・ 関数値
・ 内積
f ( x)  f , ( x)
( x), ( y)  k ( x, y)
暗黙的に，高次元の特徴空間内
でのデータ解析が可能となる．
が計算可能となる．

xj
xi 
xi 
xi
元のデータ空間
Ω
カーネル特徴空間 (RKHS)

k
射影定理（1）
7
• ある時刻 tに関して，過去の入出力 p(t )と未来の入力 u (t )
から，未来の出力 f (t )を予測する問題を考える．
未来の入力
u(t ), u(t  1),, u(t  l  1)
t
t に沿った
 u (t)
f (t )
未来の出力（真）
y(t ), y(t  1),, y(t  l  1)
fˆ (t )
未来の出力（予測）
 t への射影行列
0
 t に沿った t への射影行列
 p(t )
t
*
過去の入出力 w(t  1), w(t  2),
* 結合過程 w(t )  [u(t )', y(t )' ]'
射影定理（2）
• 次の2つの仮定をおく．
[仮定1] 出力 y から入力 へのフィードバックが存在しない


[仮定2]  t が直和分解  t   t   t を持つ（実用的には，入力
が持続的励起条件を満足していれば十分）
u
• このとき，次の射影定理が成り立つ事が導かれる．
過去の入出力と未来の入力に基づいた，未来の出力の最適予測は

fˆ (t )  f (t ) / 



p
(
t
)


u
(t ).
t
t
のように与えられる．このとき，各々の射影行列は次式で表される離
散 Wiener-Hopf 型方程式を満たす．
  pp|u   f p|u ,  uu| p   f u| p
条件付き共分散行列
1
ab|c  a  c, b  c  ab  accc
bc
8
特徴空間における射影定理
• 前褐の2つの仮定が成立する場合，正定値カーネルで定まる特
徴空間においても，同様の条件が成立する事が示せる．
• 同様の手順で，特徴空間における射影定理が導出できる．
特徴空間における過去の入出力と未来の入力に基づいた，特徴空
間における未来の出力の最適予測は


 
 
fˆ  (t )  f  (t ) / 



p
(
t
)


u (t )
t
t
のように与えられる．このとき，各々の射影行列は次式で表される
離散 Wiener-Hopf 型方程式を満たす．
  p p |u   f  p |u          
,
u u |p
f u |p
9
射影定理から分かる事
10
• 最適予測子に関する方程式 fˆ  (t )   p (t )   u (t ) を書き直すと
 y  y(t ) 
 u u(t ) 



 




p
(
t
)







y  y(t  l 1)
u u(t  l 1)


• 一方，状態空間モデルの観測方程式
 y  y(t ) 
 u u(t ) 







x
(
t
)

D





y  y(t  l 1)
u u(t  l 1)


p(t )  x(t ) となるように x(t ) ,  などを選択すればよい
状態ベクトル
• 次のように，（拡張）可観測行列，（拡張）可制御行列，および状
態ベクトルを定義する．
Lf1  f  p |u Lp1  U S V T
: Lf US1/ 2 ,  : S1/ 2V T LTp


x(t)    pp|u
  p(t)  S
1
V T Lp1|u p(t)
1/ 2
– 状態ベクトルが持つべき性質（マルコフ性）を満足する
– 前褐の射影定理と状態空間モデルとの関係が成立する．
• 特徴空間上の出力に関するモデルが得られる．


 



x
(
t

1
)

A
x
(
t
)

B

u
(
t
)

K
e (t )

u









y
(
t
)

C
x
(
t
)

D

u
(
t
)

e
(t )

u
 y
11
元の空間上の状態空間モデル
特徴空間における出力に関する状態空間モデル

 

 

 x (t  1)  A x (t )  B u u(t )  K e (t )

 







y
(
t
)

C
x
(
t
)

D

u
(
t
)

e
(t )

u
 y
t  t  0  と特徴写像の全単射性から
y(t ) / t t  C x x(t )  D u u(t ) が導ける
元の空間における出力に関する状態空間モデル

 

 

 x (t  1)  A x (t )  B u u(t )  K e (t )









y
(
t
)

C

x
(
t
)

D

u
(
t
)

e
(t )
x
u

特徴写像を明示的に含んでいるので計算不可
12
13
カーネルPCAを用いた近似
• 基本的な考え方: 特徴ベクトル z を直接用いる代わりに，カーネ
T
ル主成分で張られる空間上の特徴ベクトル  z Az を用いる
– 係数行列 Az は，例えばグラム行列の固有値分解 Gz
Az  z z1/2 として計算できる
 z  z zT
より，
⇒ （近似）特徴ベクトルが明示的に計算可能 z  z  Tz Az  Gz Az
• 明示的に計算する事が可能な，状態ベクトルおよび非線形状態
空間モデルが得られる:
経験的共分散行列
1/ 2 T ˆ1
k p(t) 特異値分解 Lˆ1 ˆ Lˆ1  Uˆ Sˆ Vˆ T により計算可
・ x(t )  Sˆ Vˆ L
p
f
f p| u p


 
 T
 T


x (t  1)  A x (t )  B Au ku u(t )  K Ae ke e (t )
・
 T

 T




y
(
t
)

C
A
k
x
(
t
)

D
A
k
u
(
t
)

e
(t )
x x
u u




14
アルゴリズム （1）
Step 1. 正則経験的共分散演算子と，その平方根行列を求める:
ˆ f f |u  (GY   I N )2  GY GU (GU   I N )2 GU GY  Lˆ f LˆTf
ˆ p p|u  (GW   I N )2  GW GU (GU   I N )2 GU GW  Lˆ p LˆTp
ˆ f p|u  GY GW  GY GU (GU  I N )2 GU GW
Step 2. 正規化共分散行列の特異値分解を計算する:
Lˆf1 ˆ f p|u Lˆp1  Uˆ Sˆ Vˆ T  U1 S1 V1T
Step 3. 状態ベクトルの推定値を計算する:
 x(l),, x(l  N )   S11/2V1T Lp1 GW
0に近い特異値は無視
アルゴリズム （2）
Step 4. グラム行列の固有値分解により係数行列 Au , Au , Ax を計
算し，次式に最小二乗法を適用してシステム行列を求める．
Xˆ l 1  A  Xˆ l  B  AuT Gˆu,l  w
ˆ x,l  D  AuT G
ˆ u,l  e
Yl|l  C  AxT G
Step 5（カルマンゲインが必要な場合）. 残差の共分散行列を計算
し，代数リッカチ方程式（本文参照,MATLABで一発）を解き，そ
の安定化解を用いてカルマンゲインを求める．
15
16
数値例（１）
• 下記の非線形システムから生成したデータを利用
– 0.05秒毎の4,5次のルンゲ・クッタによるシミュレーション
– 600点を用いて学習し，新たな400点により評価*
– 入力は  0.5～0.5 の間の均一分布から生成
x1 (t )  x2 (t )  0.1cos(x1 (t )) (5x1 (t )  4x13 (t )  x15 (t ))  0.5 cos(x1 (t )) u(t )
x2 (t )  65x1 (t )  50x13 (t ) 15x15 (t )  x2 (t ) 100u(t )
y (t )  x1 (t )
*) 評価は主に次式で表されるシミュレーション誤差を利用 [OM96]

m
s
n
100 y i 1 ( yi )c  ( y i )c
 m ( y ) 
ny c1
i 1
i c

シミュレーション値
観測値
数値例（１）
17
獲得モデルを用いた長期予測の結果
線形部分空間同定法 [KP99]
非線形部分空間同定法 （提案手法）
シミュレーション誤差 : 45.1
〃 : 40.2 （約10％向上）
数値例（２）
• Simulation data of a pH neutralization process in a constant
volume stirring tank.
– Included in DAISY (DAtasets for the Identification of Systems, which
includes several engineering, biological and environmental data).
– Inputs: acid solution flow and base solution flow in litters
Outputs: pH of the solution in the tank
18
数値例（２）
19
獲得モデルを用いた長期予測の結果
線形部分空間同定法 [KP99]
非線形部分空間同定法 （提案手法）
シミュレーション誤差 : 47.0
〃 : 22.7 （約50％向上）
Conclusions and future works
20
• A kernel subspace method based on stochastic realization for
learning nonlinear dynamical systems is proposed.
– The algorithm needs no iterative optimization procedures and can obtain
solutions using fast and reliable numerical schemes (SVD, etc.).
– However, the parameters involve much time and effort for tuning.
• Future works :
– To incorporate a scheme for optimizing, automatically, the parameters into
the proposed method.
– To extend other established subspace methods for learning dynamical
systems to nonlinear frameworks as well.