はじめに

リアプノフスペクトル（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)$ の時間発展を考えます。これは、初期時刻 $0$ から時刻 $t$ までの系の線形応答を表す行列 $\boldsymbol{M}(t)$ を用いて次のように書けます¹。

$$ \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$ を用いて、Oseledecの多重エルゴード定理から、次のように収束することが知られています。

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

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

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 \\

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

\end{split} \tag{2}

と変形できます。ここで、上三角行列の積を$\mathbf{R}(t)= \mathbf{R}_t \cdots \mathbf{R}_2\mathbf{R}_1$とおきました。また、行列 $\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}$ です。ゆえに、次のように計算できます。

\begin{split}

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

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

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


\end{split} \tag{2}

ここで、$\boldsymbol{\delta}_i$は$\mu_i$に対応する固有ベクトルです。この基底$\boldsymbol{\delta}_i$に対して、その拡大率をリアプノフスペクトラム$\lambda_i$から計算できます²。（Oseledecの多重エルゴード定理は、リアプノフスペクトラムの収束と存在を保証してくれています）
具体的には、$m$個の基底ベクトル$\boldsymbol{\delta}_i$が張る平行多面体の$m$次元体積 $V_m$ の指数的拡大率は

\begin{split}
    \lim_{t \to \infty} \frac{1}{t} \log V_m = \sum_{i=1}^m \lambda_i
\end{split} \tag{3}

と表せます。
一方、行列 $\mathbf{M}$ によるこの拡大率は、$\mathbf{R}(t)$ の左上の$m \times m$小行列$\mathbf{R}^{mm}(t)$の行列式によって与えられます。（$\mathbf{Q}$は回転や対称移動のみで、体積に影響しない）

\begin{split}
\lim_{t \to \infty} \frac{1}{t} \log V_m 
&= \lim_{t \to \infty} \frac{1}{t} \log| \det(\mathbf{R}^{mm}(t))|\\
&= \lim_{t \to \infty} \frac{1}{t} ( \log|r^{11}(t)|+\log|r^{22}(t)|+\cdots+\log|r^{mm}(t)|)
\end{split} \tag{4}

と計算できます。ここで、$r^{ii}(t)$は$\mathbf{R}(t)$ の対角成分です。

(3,4)式から、
$m=1$のとき、$\lambda_1=\lim_{t \to \infty} \frac{1}{t} \log|r^{11}(t)|$、
$m=2$のとき、$\lambda_1+\lambda_2=\lim_{t \to \infty} \frac{1}{t} (\log|r^{11}(t)|+\log|r^{22}(t)|),~ \therefore\lambda_2=\lim_{t \to \infty} \frac{1}{t} \log|r^{22}(t)|$
$m=\cdots$と逐次的に計算することで、$m$番目のリアプノフ指数 $\lambda_m$ は、

\begin{split}

    \lambda_i
    
    &=\lim_{t \to \infty} \frac{1}{t} \log|r^{ii}(t)| \\
\end{split}

と求められます。
さらに、$\mathbf{R}(t)= \mathbf{R}_t \cdots \mathbf{R}_2\mathbf{R}_1$を戻して、

\begin{split}
    \lambda_i
    
    &=\lim_{t \to \infty} \frac{1}{t} \log|r^{ii}(t)| \\

    &=\lim_{t \to \infty} \frac{1}{t} \log|r^{ii}_t\cdots r^{ii}_1| \\

    &=\lim_{t \to \infty} \frac{1}{t}
    
    \sum_{k=1}^t

    \log{|r_{k}^{ii}|}
\end{split}

と書き直せます。結局のところ各時刻 $k$ の対角成分 $r_{k}^{ii}$ を計算すればよいことが分かります。$\log{|r_{k}^{ii}|}$の平均を取りながら計算することで、オーバーフローも回避できます。

実装（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. ↩
Pikovsky, Arkady, and Antonio Politi. Lyapunov exponents: a tool to explore complex dynamics. Cambridge University Press, 2016. ↩
Sano, Masaki, and Yasuji Sawada. "Measurement of the Lyapunov spectrum from a chaotic time series." Physical review letters 55.10 (1985): 1082. ↩