時系列データ解析 基礎編 其の八 - Lyapunov Spectrum (2008年01月20)

前回はシステムの初期値依存性の強さを表す Lyapunov Exponent について説明したが、それを成分に分解した(ような) Lyapunov Spectrum というものもある。 2点間の距離がどの方向にどれだけ伸縮されたかを成分に分けて計算しようというものだ。 ところが、アトラクターの伸縮は当然局所的なものであって Lyapunov Exponent のように大局的な数値を計算するのは不可能なように思える。 Lyapunov Spectrum の考え方は、簡単に言うと、局所的な伸縮をある基準点に揃えて 伸縮の様子を分解してみようというものだ。

上記の Lyapunov Spectrum アイデアを式に書くと以下のようになる。

(35)

ただし、JnXn の近傍で求めた ヤコビ行列(Jaconbian)、ui は固有ベクトル(eigenvector)、Λi は固有値(eigenvalue)。 ヤコビ行列(Jaconbian)を J+nJn のように二乗にしているわけは 固有値を実数にするため。

注意しなければいけないのは、Jn の固有ベクトルは各 Xn の近傍毎に異なるため 単純に Jn を掛ける訳にはいかない。 掛ける前に適当な変換を行って 固有値の絶対値の大きい(小さい)順に固有ベクトルの方向を揃える必要がある。

Jn の計算は そう単純ではない。 m 次元の埋め込み空間(embedded space)の場合 Jnm2 の成分を持つ。 Jn を求めるには 各成分が Xn の近傍で最適値を取るようにオプティマイゼイションを行う必要がある。

具体的は、Xn の全ての近傍点 Xk について 下の χ を最小にするように Jn(τ) の各成分を決定する。

(36)

(37)

ここで、U(Xn)Xn の近傍を表し、Xk は その中に含まれる点、τ は遅延時間を表す。

遅延時間 τ を突然入れた理由は、後々 Lyapunov Spectrum の時間発展を観るため、時間経過も含めた Jn(τ) を計算することを明示するためだ。

ところが、(35)にはコンピュータで計算する上で致命的な問題点がある。。。。N が大きくなると ΠJn が発散するのだ。

そこで、今回は以下のように積を和に変えて計算した。 実際これが正しい処理だったのか ちょっと不安が残る。

(38)

やっていることは

  1. ある Xn について Jn(τ) を求める
  2. Jn(τ) から 固有ベクトル(eigenvector) un,i(τ) と 固有値(eigenvalue) Λn,i(τ) を計算する
  3. 固有値(eigenvalue) の大きい順に 固有ベクトル(eigenvector) と 固有値(eigenvalue) を並べる
  4. 全ての Xn について 1-3 を繰り返す
  5. 固有値(eigenvalue) Λn,i(τ) から Lyapunov Spectrum λi(τ) の平均値を計算する

式(13)、(14)で与えられる時系列データ 図-13図-17 について Lyapunov Spectrum を求めた結果が以下の図のようになる。 ただし、埋め込み空間(embedded space) の次元は m=3 とした。 また、参考のために Lyapunov Exponent もプロットした。

図-24

Lyapunov Spectrum は指数であるので、式(13)、(14)で与えられる時系列データは、1方向に伸ばされ、1方向に縮められ、残るもう一つの方向にはほとんど変化が無いということがわかる。

Lyapunov Spectrum の値が Lyapunov Exponent に比べて極端に大きいが、これはXnXkXn+1Xk+1 等はそれぞれ近い距離に位置するのだけど、XnXn+1XkXk+1 が近い位置にあるとは限らないことによる。

下の図は、最大の 固有値(eigenvalue) を持つ 固有ベクトル(eigenvector) un,1(0) を3次元空間内にステレオグラム(平行法)でプロットしたものだ。 固有ベクトル(eigenvector) 自体は単位長さのベクトルであるので、全ての点は半径1の球面上にプロットされている。 アトラクター(図-17)がどの方向に伸ばされているかが分かる。 ただし、基準となる方向は適当とったものであり、図-17の座標系には合っていない。

図-25