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.