はじめに

リアプノフスペクトル（Lyapunov Spectrum）を具体的に計算・近似する方法についてのメモです。
途中までの理論・計算方法について、Wikipediaの解説を引用しながら説明します。
その後、本記事の目的であるQR分解によるリアプノフスペクトルの計算方法と実装について解説します。

リアプノフ指数（Lyapunov Exponent）

リアプノフ指数とは、初期値の差が非常に小さい2つの軌道が指数関数的に離れていく度合いを表す量です。
時刻$t$における誤差を$ \boldsymbol{\delta}(t)$として、リアプノフ指数$\lambda$は

\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)$ の解を次の形式で表し、

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

$\boldsymbol{\delta}(0)$ に対する乗数 $\boldsymbol{M}$ の特異値の自然対数をとり、降順に並べたものがリアプノフスペクトルです。
$n$次元離散時間自励力学系は$\boldsymbol{x}(t) \in \mathbb{R}^{n} , \boldsymbol{f}: \mathbb{R}^{n} \xrightarrow{} \mathbb{R}^{n}$を用いて、

    \boldsymbol{x}(t+1) = \boldsymbol{f}[\boldsymbol{x}(t)]

と表されます。$\boldsymbol{f}$ が線形近似可能な場合、$\boldsymbol{x}(t)$におけるヤコビ行列$\mathrm{D}\boldsymbol{f}(\boldsymbol{x}(t))$を用いて、

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

と表されます。したがって、$\boldsymbol{M}$ は $\mathrm{D}_{\boldsymbol{x}}\boldsymbol{f}$ の 0 から t − 1 までの総乗として得ることができるので、

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

と表すことができます。
結局、$\boldsymbol{M}$ は、

    \boldsymbol{M}=\mathrm{D}\boldsymbol{f}(\boldsymbol{x}(t-1))
    \cdots \mathrm{D}\boldsymbol{f}(\boldsymbol{x}(0)) \tag{1}

と表せ、$t$ が十分大きいとき、$\boldsymbol{M}$ の特異値の自然対数をとったものは、リアプノフスペクトルを近似できています。（ここまではWikipediaの解説と同じ）
しかし、データ数 $t$ が有限の値であっても、正のリアプノフ指数が存在すれば、式(1)は数値計算上オーバーフローします。また、式(1)の右辺の積をとる過程で、小さなリアプノフ指数に対応する部分空間は数値的に零空間と区別できなくなってしまいます。
これらの問題を回避するため、QR分解を用いた手法がLyapunovスペクトラムの計算で広く用いられています。

QR分解（復習）

QR分解とは、ある行列$A$を直交行列$Q$と上三角行列$R$に分解することです。

    \mathbf{A} = \mathbf{Q}\mathbf{R}

このとき、行列$A$と上三角行列$R$の特異値は一致します。

Mの特異値をQR分解で求める

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

    \mathrm{D}\boldsymbol{f}(\boldsymbol{x}(0))
    = \mathbf{Q}_1\boldsymbol{R}_1

とQR分解します。
すると，式(1)は

    \boldsymbol{M}
    =\mathrm{D}\boldsymbol{f}(\boldsymbol{x}(t-1))
    \cdots \mathrm{D}\boldsymbol{f}(\boldsymbol{x}(1))
    \mathbf{Q}_1\mathbf{R}_1

となります。
さらに、

    \mathrm{D}\boldsymbol{f}(\boldsymbol{x}(i))\mathbf{Q}_i
    = \mathbf{Q}_{i+1}\mathbf{R}_{i+1}

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

\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_{t}…R_{2}R_{1}$にQR分解したといえます。元の行列の特異値と上三角行列の特異値は一致するので、上三角行列 $R_{t}…R_{2}R_{1}$の特異値を求めればよく、これは上三角行列 $R_{t}…R_{2}R_{1}$の対角成分です。

実際に計算する場合はまず、$\mathbf{R}_k$の対角成分 $[r_1^k, r_2^k, \cdot\cdot\cdot, r_n^k]$ を式(2)の計算を行いながら取り出します。
Lyapunovスペクトラム $[\lambda_1, \lambda_2, \cdot\cdot\cdot, \lambda_n]$ は、$\boldsymbol{M}$ の特異値の自然対数なので、

    [\lambda_1, \lambda_2, \cdot\cdot\cdot, \lambda_n]
    \simeq \frac{1}{t}
    \sum_{k=1}^t
    [\:\log{||r_1^k||}, \log{||r_2^k||}, \cdot\cdot\cdot, \log{||r_t^k||} \:]

と近似することができました。

実装（Python）

上記のリアプノフスペクトルを求めるプロセスをPythonで書いてみました。
Github上に公開しています。
lyapunov_spectrum.pyでは、上記のプロセスが行われます。

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写像でのリアプノフスペクトラムが求まります。
henon.pyの代わりに求めたい力学系があれば、置換してください。

応用

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

連続時間力学系

上記の手法は離散時間力学系の場合で考えていました。時系列と各点におけるヤコビ行列を求めるためには、自励系の場合はscipy.integrate.odeintを逐次的に用いるとよいでしょう。後は、離散時間力学系の時と同様に、lyapunov_spectrum.pyに出力させるだけです。

決定論的法則が未知でアトラクタ上の【時系列データだけ】の場合

アトラクタ（十分な時間が経過したときに軌道が引き付けられる相空間上の有界な集合）上の十分長い時系列データだけある場合にリアプノフスペクトルを求める方法としてSano-Sawadaの方法²というものがあります。
今、十分長い時系列データは既知なので、あとは各点におけるヤコビ行列を求めればよいです。そこで、ヤコビ行列を軌道の時系列から近似することを考えます。

上図のように、ヤコビ行列を求めたい点 $x_k$ を中心とした $\epsilon$-近傍を考えます。この近傍内の点 $x_j, x_{j'}$ について、$\epsilon$ が十分小さいとき、以下に示すTaylor展開が成立します。

    \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}$ が求まる事が分かります。

リザバーコンピューティングの力学的性質の解析

記事を書きましたのでご覧ください。

いかがでしたか？

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

Collet, Pierre, and Jean-Pierre Eckmann. "Concepts and results in chaotic dynamics: a short course." Springer Science & Business Media, 2007. ↩
Sano, Masaki, and Yasuji Sawada. "Measurement of the Lyapunov spectrum from a chaotic time series." Physical review letters 55.10 (1985): 1082. ↩