リアプノフスペクトル（Lyapunov Spectrum）のQR分解を用いた求め方

# はじめに
リアプノフスペクトル（Lyapunov Spectrum）を具体的に計算・近似する方法についてのメモです。
途中までの理論・計算方法について、[Wikipediaの解説](https://ja.wikipedia.org/wiki/%E3%83%AA%E3%82%A2%E3%83%97%E3%83%8E%E3%83%95%E6%8C%87%E6%95%B0)を引用しながら説明します。
その後、本記事の目的であるQR分解によるリアプノフスペクトルの計算方法と実装について解説します。


## リアプノフ指数（Lyapunov Exponent）
リアプノフ指数とは、初期値の差が非常に小さい2つの軌道が指数関数的に離れていく度合いを表す量です。
時刻$t$における誤差を$ \boldsymbol{\delta}(t)$として、リアプノフ指数$\lambda$は
```math
\begin{equation}
     ||\boldsymbol{\delta}(t)|| \approx e^{\lambda} ||\boldsymbol{\delta}(0)||
\end{equation}
```
と近似的に表されます。
相空間が $n$ 次元の場合、リアプノフ指数は $n$ 個存在することになります。
これらを降順に並べた${ \lambda_1, \lambda_2, \cdot\cdot\cdot, \lambda_n }$をリアプノフスペクトル（Lyapunov Spectrum）といいます。


## リアプノフスペクトルを求める

まず線形化された誤差ベクトル $\boldsymbol{\delta}(t)$ の時間発展を考えます。これは、初期時刻 $0$ から時刻 $t$ までの系の線形応答を表す行列 $\boldsymbol{M}(t)$ を用いて次のように書けます[^1]。

$$ \boldsymbol{\delta}(t) = \boldsymbol{M}(t) \boldsymbol{\delta}(0) $$

$n$次元離散時間力学系が、状態ベクトル $\boldsymbol{x}(t) \in \mathbb{R}^{n}$ と写像 $\boldsymbol{f}: \mathbb{R}^{n} \to \mathbb{R}^{n}$ を用いて、
$$ \boldsymbol{x}(t+1) = \boldsymbol{f}(\boldsymbol{x}(t)) $$
と表される場合を考えます。基準となる軌道 $\boldsymbol{x}(0), \boldsymbol{x}(1), \dots$ からの微小なずれ $\boldsymbol{\delta}(t)$ の時間発展は、$\boldsymbol{f}$ が滑らかであれば、各点 $\boldsymbol{x}(t)$ におけるヤコビ行列 $Df(\boldsymbol{x}(t))$ を用いて線形近似できます。

$$ \boldsymbol{\delta}(t+1) \approx Df(\boldsymbol{x}(t)) \boldsymbol{\delta}(t) $$

これを時刻 $0$ から $t$ まで繰り返し適用すると、$\boldsymbol{\delta}(t)$ は $\boldsymbol{\delta}(0)$ を用いて次のように表されます。

$$ \boldsymbol{\delta}(t) \approx Df(\boldsymbol{x}(t-1)) \dots Df(\boldsymbol{x}(1)) Df(\boldsymbol{x}(0)) \boldsymbol{\delta}(0)$$

したがって、時刻 $t$ における線形応答を表す行列 $\boldsymbol{M}(t)$ は、$t$ ステップ分のヤコビ行列の積として具体的に与えられます。

$$ \boldsymbol{M}(t) = Df(\boldsymbol{x}(t-1)) \dots Df(\boldsymbol{x}(1)) Df(\boldsymbol{x}(0)) \tag{1} $$

リアプノフスペクトル $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_n$ は、この行列 $\boldsymbol{M}$ の特異値 $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_n \ge 0$ を用いて、以下の極限で定義されます。

$$ \lambda_i = \lim_{t \to \infty} \frac{1}{t} \ln \sigma_i \quad (i=1, \dots, n) $$
(ここで $\sigma_i$ は $t$ ステップ後の行列 $\boldsymbol{M}$ の $i$ 番目の特異値です。)


（ここまでは[Wikipediaの解説](https://ja.wikipedia.org/wiki/%E3%83%AA%E3%82%A2%E3%83%97%E3%83%8E%E3%83%95%E6%8C%87%E6%95%B0)と同じ）
しかし、データ数 $t$ が有限の値であっても、正のリアプノフ指数が存在すれば、式(1)は数値計算上オーバーフローします。また、式(\ref{1})の右辺の積をとる過程で、小さなリアプノフ指数に対応する部分空間は数値的に零空間と区別できなくなってしまいます。
これらの問題を回避するため、QR分解を用いた手法がLyapunovスペクトラムの計算で広く用いられています。


## QR分解で求める

まず，$\mathrm{D}\boldsymbol{f}(\boldsymbol{x}(0))$を

```math

    \mathrm{D}\boldsymbol{f}(\boldsymbol{x}(0))

    = \mathbf{Q}_1\boldsymbol{R}_1

```

とQR分解します。

すると，式(1)は

```math

    \boldsymbol{M}

    =\mathrm{D}\boldsymbol{f}(\boldsymbol{x}(t-1))

    \cdots \mathrm{D}\boldsymbol{f}(\boldsymbol{x}(1))

    \mathbf{Q}_1\mathbf{R}_1

```

となります。

さらに、

```math

    \mathrm{D}\boldsymbol{f}(\boldsymbol{x}(i))\mathbf{Q}_i

    = \mathbf{Q}_{i+1}\mathbf{R}_{i+1}

```

と逐次QR分解を行うことで

```math

\begin{split}

    \boldsymbol{M}

    &=\mathrm{D}\boldsymbol{f}(\boldsymbol{x}(t-1))

    \cdots \mathrm{D}\boldsymbol{f}(\boldsymbol{x}(1))

    \mathbf{Q}_1

    \boldsymbol{R}_1 \\

    &=\mathrm{D}\boldsymbol{f}(\boldsymbol{x}(t-1))

    \cdots 

    \mathbf{Q}_2\mathbf{R}_2\mathbf{R}_1 \\

    & \vdots \\

    &=\mathrm{D}\boldsymbol{f}(\boldsymbol{x}(t-1))\mathbf{Q}_{t-1}

    \cdots

    \mathbf{R}_2\mathbf{R}_1 \\

    &=\mathbf{Q}_{t}\mathbf{R}_{t}

    \cdots 

    \mathbf{R}_2\mathbf{R}_1 \\

\end{split} \tag{2}

```

と変形できます。

このとき、上三角行列同士の積は上三角行列なので、$\boldsymbol{M}$を直交行列$Q_{t}$と上三角行列 ${R} = R_{t}…R_{2}R_{1}$にQR分解したといえます。また、行列 $\boldsymbol{M}$ の特異値 $\sigma_i$ は、行列 $\boldsymbol{M}^\top \boldsymbol{M}$ の固有値を $\mu_i$ ($\mu_1 \ge \dots \ge \mu_n \ge 0$) としたとき、その非負の平方根 $\sigma_i = \sqrt{\mu_i}$ です。ゆえに、次のように計算できます。
```math
\begin{split}

    \sigma_i
    &= (\mu_i)^{1/2}\\
    &= (\boldsymbol{e}_i^\top \boldsymbol{M}^\top \boldsymbol{M}\boldsymbol{e}_i )^{1/2}\\

    &= (\boldsymbol{e}_i^\top 
    \mathbf{R}^\top \mathbf{Q}_{t}^\top
    \mathbf{Q}_{t}\mathbf{R}\boldsymbol{e}_i )^{1/2}\\

    &= |\mathbf{R}\boldsymbol{e}_i| \\

    &= |r^{ii}|
    


\end{split} \tag{2}

```
ここで、$\boldsymbol{e}_i$は$\mu_i$に対応する固有ベクトルで、$r^{ii}$ は $\mathbf{R}$ の $i$ 番目の対角成分を表します。Lyapunovスペクトラム $[\lambda_1, \lambda_2, \cdot\cdot\cdot, \lambda_n]$ は、$\boldsymbol{M}$ の特異値の自然対数なので、



```math
\begin{split}

    \lambda_i

    &=\lim_{t \to \infty} \frac{1}{t} \ln \sigma_i \\

    &\simeq \frac{1}{t}

    \sum_{k=1}^t

    \log{||r_{k}^{ii}||}
\end{split}
```
と近似することができました。
この計算から、結局のところ各時刻 $k$ の対角成分 $r_{k}^{ii}$ を計算すればよいことが分かります。

## 実装（Python）
[上記](#Mの特異値をQR分解で求める)のリアプノフスペクトルを求めるプロセスをPythonで書いてみました。
[Github](https://github.com/OEC-ket-bra/Lyapunov-spectrum)上に公開しています。
`lyapunov_spectrum.py`では、[上記](#Mの特異値をQR分解で求める)のプロセスが行われます。
```python:lyapunov_spectrum.py
    def __call__(self, time_series, ls_range):
         """
        Args:
            time_series (array): 時系列データ（N_x * データ長）
            ls_range (int): 何番目までのリアプノフスペクトルを求めるか
        Returns:
            array: リアプノフスペクトル
        """
        Q = np.eye(self.N_x) # 単位行列
        le_bank = np.zeros(ls_range) #零配列 

        data_length = len(time_series) # データ長

        for x in time_series: # (2)式のプロセス
            Q, R =  np.linalg.qr(self.jacobian(x) @ Q) # QR分解
            R_diag = np.diag(R) # 対角成分の抽出
            le_bank = le_bank + np.log(np.abs(R_diag[:ls_range])) # 対角成分のスタック
            
        return le_bank / data_length
```
`main.ipynb`を上から実行すれば、例として[Henon写像](https://ja.wikipedia.org/wiki/%E3%82%A8%E3%83%8E%E3%83%B3%E5%86%99%E5%83%8F)でのリアプノフスペクトラムが求まります。
`henon.py`の代わりに求めたい力学系があれば、置換してください。

## 応用
結局この手法では、**十分長い時系列データ**と**各点におけるヤコビ行列**があればリアプノフスペクトルを求められることが分かります。したがって、以下のような応用があります。

#### 連続時間力学系
上記の手法は離散時間力学系の場合で考えていました。**時系列**と**各点におけるヤコビ行列**を求めるためには、自励系の場合は[scipy.integrate.odeint](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html)を逐次的に用いるとよいでしょう。後は、離散時間力学系の時と同様に、`lyapunov_spectrum.py`に出力させるだけです。

#### 決定論的法則が未知でアトラクタ上の【時系列データだけ】の場合
[アトラクタ](https://ja.wikipedia.org/wiki/%E3%82%A2%E3%83%88%E3%83%A9%E3%82%AF%E3%82%BF%E3%83%BC)（十分な時間が経過したときに軌道が引き付けられる相空間上の有界な集合）上の**十分長い時系列データ**だけある場合にリアプノフスペクトルを求める方法としてSano-Sawadaの方法[^2]というものがあります。
今、**十分長い時系列データ**は既知なので、あとは**各点におけるヤコビ行列**を求めればよいです。そこで、ヤコビ行列を軌道の時系列から近似することを考えます。

![名称未設定のノート (2)-13 (1).jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3340883/3b371c39-4b18-6baa-8910-310d65f63f51.jpeg)

上図のように、ヤコビ行列を求めたい点 $x_k$ を中心とした $\epsilon$-近傍を考えます。この近傍内の点 $x_j, x_{j'}$ について、$\epsilon$ が十分小さいとき、以下に示すTaylor展開が成立します。
```math
    \boldsymbol{x}_{k+1} = \boldsymbol{f}[\boldsymbol{x}_k] = \boldsymbol{f}[\boldsymbol{x}_j] + \mathrm{D}_{\boldsymbol{x}_j}\boldsymbol{f}(\boldsymbol{x}_j-\boldsymbol{x}_k) + O(||\boldsymbol{x}_j-\boldsymbol{x}_k||^2) \tag{3}
```
式(3)は近傍内の点の個数だけ取れるので、相空間の次元の数よりも多ければ、最小二乗法などで $\mathrm{D}_{\boldsymbol{x}_j}\boldsymbol{f}$ が求まる事が分かります。

#### リザバーコンピューティングの力学的性質の解析
[記事](https://qiita.com/OEC_ket_bra/items/174b72f1708309453b41)を書きましたのでご覧ください。

# いかがでしたか？
この記事が気に入ったなら是非、記事の上下にある**いいね**💓を押して下さいね！それでは皆さんもぜひリアプノフスペクトルをごりごり計算してください！
（[応用](#応用)の内容の実装はそのうちやります。）

[^1]: Collet, Pierre, and Jean-Pierre Eckmann. "Concepts and results in chaotic dynamics: a short course." Springer Science & Business Media, 2007.
[^2]: [Sano, Masaki, and Yasuji Sawada. "Measurement of the Lyapunov spectrum from a chaotic time series." Physical review letters 55.10 (1985): 1082.](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.55.1082)