とつとつとしてろうとせず

ひまつぶしにどうぞ。

動的モード分解に関する覚書

注意

以下の記述は私的な備忘録であり、内容の正確さに関しては一切保証いたしませんので、悪しからず。

はじめに

流体解析などの分野でデータ解析や縮約モデリングの手法として用いられる動的モード分解の概要についてまとめる。近年、PODと並んで注目されている?手法であり、主に海外で盛んに研究されているようだ(国内での流行り具合はどうなんだろ)。実際には分野や状況に応じてPODと使い分けていくことになるのかな。
関連手法であるKoopmanモード分解については別記事でまとめる予定。

元ネタ

DMDの概要と実行例が載っている素晴らしいサイト。
基本的にはこちらを見てもらえばいいと思うが、PODの場合と違ってネットに日本語の解説があまりないようなので、勉強として自分でもまとめておく。

POD

まずはPOD(固有直交分解)についておさらいすると、時系列データを並べて作ったデータ行列に対して主成分分析を行い、データから特徴的な構造(固有モード)を抽出するのがPODなのだった(名前の由来はPODモードが正規直交基底となるため)。
具体的な計算としてはデータ行列を特異値分解することで行われる。特異値の大きなモードを支配的なモードと解釈し、主要なモードのみを抜き出すことで元のデータからデータ容量を削減した縮約モデル(ROM)を作り、低コストなシミュレーションを行うこともできる。
また、時系列データだけでなく、パラメータを変えて得られたデータに対しても適用できるので、最初に代表的なパラメータに対してデータをサンプリングしておけば、パラメータが動く範囲内に対してはROMによって高速なシミュレーションができる(データを補間するイメージ)。

DMD

PODによる流れ場のモード分解では、時間と空間の相関を同時に計算するため、時空間のモードがうまく分離されず、特に時間的な振る舞いが複数のモードに分散して抽出されてしまう場合がある(系を特徴付ける1つの周波数を再現するために複数の自明でないモードの組み合わせが必要になる)。
動的モード分解(Dynamic Mode Decomposition, DMD)では、時系列データの時間発展に注目してモード分解を行うことで、固有モードとそれに対応する周波数を得ることができるので、あるモードが時間的に成長するのか減衰するのか、または周期的に振動するのかといった、いわゆる安定性に関する情報を得ることもできる。

以下、具体的なアルゴリズムについてまとめる。定式化としてはKrylov部分空間法に基づくものと特異値分解に基づくものがあるが、数学的な解析には前者が使われ、数値計算には後者が用いられるようなので、ここでは特異値分解に基づくアルゴリズムについて述べる。

離散力学系の時間発展

PODと同様、データとしては各時刻におけるスナップショットを1次元のベクトル(状態ベクトル)としてまとめたものを用いる。データ数を m個とすると、このデータベクトル$\mathbf{x}_k$の時間発展は次式のように表される。


\begin{align}
\mathbf{x}_{k+1} = \mathbf{f} \paren{ \mathbf{x}_{k} } \quad \paren{ k = 0, 1, 2, \cdots , m-1 }
\end{align}

ここで、$\mathbf{f}$は系の時間発展を表す関数のベクトルであり、一般に非線形で場合によっては未知である。上式は数学的には離散力学系と呼ばれる定式化になっており、ある時刻の状態が次の時刻にどのような状態に移るのかを表す(連続力学系の場合は上式は時間に関する一階の微分方程式になる)。
DMDはこの時間発展をデータから推定するという発想に基づいており、データ駆動型(data-driven )で方程式によらない(equation free)手法といわれる。つまり、物理モデルが存在しない未知の現象やそもそもモデル化が困難な経済学的な現象(株価の変動とか?)など、$\mathbf{f}$が具体的に分からない現象に対しても、データだけを使ってその振る舞いを分析・予測しようというアプローチを取っている。

さて、上式をまとめて行列で書くと、次式のようになる。


\begin{align}
\mathbf{X}' = \mathbf{F} \paren{ \mathbf{X} }
\end{align}

ここで、データ行列 \mathbf{X} \in \mathbb{R}^{n \times m-1} \mathbf{X}' \in \mathbb{R}^{n \times m-1}は次式で定義される。


\begin{align}
\mathbf{X} = 
\begin{bmatrix}
\vert & \vert  &   & \vert  \\
\mathbf{x}_{0} & \mathbf{x}_{1} & \cdots & \mathbf{x}_{m-2} \\
\vert & \vert  &   & \vert 
\end{bmatrix}
\\
\mathbf{X}' = 
\begin{bmatrix}
\vert & \vert  &   & \vert  \\
\mathbf{x}_{1} & \mathbf{x}_{2} & \cdots & \mathbf{x}_{m-1} \\
\vert & \vert  &   & \vert 
\end{bmatrix}
\end{align}

つまり、上式は元の m個のデータを並べた行列から、時刻を1つずつずらした行列になっており、$\mathbf{X}$と$\mathbf{X}'$の相関関係から時間発展を求めていくことになる。そこで、DMDでは上式を次式で近似する。


\begin{align}
\mathbf{X}' \approx \mathbf{A} \mathbf{X}
\end{align}

すなわち、未知の非線形関数を線形の行列(線形写像)で近似する。元々線形ならともかく、非線形の場合にもこういうことをして上手くいくのかという疑問は当然であって、実際上手くいかないこともあるようだ。ただし、ある状況下では非線形の時間発展を線形近似しても上手くいくことが知られており、この辺りの話はKoopman作用素によるモード分解と関連している(詳細については別記事で触れる)。

力学系の安定性

ちょっと脱線して力学系の話についてまとめる。門外漢のため、俺も詳しくはないのだが、この辺の概念がDMDアルゴリズムや結果を理解するのに必要になるっぽいので、最低限イメージだけでも掴めるようにメモっておく。

連続力学系

まずは連続力学系から考える。一般の状態ベクトルの時間発展は次の微分方程式で表される。


\begin{align}
\dif{}{t} \mathbf{x} = \mathbf{f} \paren{ \mathbf{x}} 
\end{align}

右辺が$t$に陽に依存する場合は非自律系(非自励系)、そうでない場合は自律系(自励系)という。以降、簡単のため自律系を考える。
上述のように \mathbf{f}は一般には非線形であるため、そのまま上式を解くことはできない。そこで、速度が0になる点(平衡点)の周りで上式を線形化し、近似的にその振る舞いを調べるということが行われる。したがって、適当な線形化の下で次式を得る。


\begin{align}
\dif{}{t} \mathbf{x} = \mathbf{A} \mathbf{x} 
\end{align}

上式は線形の微分方程式であり、容易に次のような解を得ることができる。


\begin{align}
\mathbf{x}(t) = e^{\mathbf{A} t} \mathbf{x}(0)
\end{align}

ここで、 e^{\mathbf{A} t}は行列指数関数と呼ばれ、具体的には \mathbf{A}に関する固有分解を \mathbf{A} = \mathbf{U}\boldsymbol{\Lambda}\mathbf{U}^{-1}として、次式で与えられる。

 
\begin{align}
\begin{aligned}
e^{\mathbf{A} t} &= \mathbf{U} e^{\boldsymbol{\Lambda} t} \mathbf{U}^{-1} \\
e^{\boldsymbol{\Lambda} t} &= \operatorname{diag} \paren{ e^{\lambda_{j} t} }
\end{aligned}
\end{align}

つまり、 \mathbf{A}が対角化可能であれば、行列指数関数は \mathbf{A}固有値固有値ベクトルから計算できる(対角化ができなくてもJordan標準形を用いることで同様に計算可能)。

指数関数の性質とEulerの公式から、 \lambda_{j} = \alpha_{j} + i \omega_{j}とすれば、

 
\begin{align}
e^{\lambda_{j} t} = e^{\alpha_{j} t} \paren{ \cos\paren{\omega_{j} t} + i \sin\paren{\omega_{j} t} }
\end{align}

であり、すべての固有値の実部が負( \alpha_{j} < 0)なら系は安定であり、そうでなければ不安定である。したがって、線形系(平衡点)の安定性は行列の固有値から判定できる。

離散力学系

離散力学系に対しても同様に考えることができる。この場合、平衡点に対応するのは \mathbf{f} \paren{ \mathbf{x} } = \mathbf{x}を満たす固定点(あるいは不動点)である。これは差分商(というかEuler陽解法?)


\begin{align}
\dif{}{t} \mathbf{x} \approx \frac{ \mathbf{x}_{k+1} - x_{k} }{ \Delta t } = \frac{ \mathbf{f} \paren{ \mathbf{x}_{k} } - x_{k} }{ \Delta t } = 0
\end{align}

を考えると分かりやすいかもしれない。

さて、適当な固定点で線形化すると次式を得る。


\begin{align}
\mathbf{x}_{k+1} = \mathbf{A} \mathbf{x}_{k}
\end{align}

上式は差分方程式であり、時刻 kの解 \mathbf{x}_{k}は初期値から繰り返し行列を作用させることで次のように得られる。

 
\begin{align}
\begin{aligned}
\mathbf{x}_{1} &= \mathbf{A} \mathbf{x}_{0} \\
\mathbf{x}_{2} &= \mathbf{A} \mathbf{x}_{1} = \mathbf{A}^{2} \mathbf{x}_{0} \\
\mathbf{x}_{3} &= \mathbf{A} \mathbf{x}_{2} = \mathbf{A}^{3} \mathbf{x}_{0} \\
& \vdots \\
\therefore \mathbf{x}_{k} &= \mathbf{A}^{k} \mathbf{x}_{0}
\end{aligned}
\end{align}

行列のべき乗は対角化によって簡単に計算することができ、


\begin{align}
\mathbf{A}^{k} = \mathbf{U} e^{\boldsymbol{\Lambda}^{k} } \mathbf{U}^{-1}
\end{align}

となるので、離散力学系の場合、固定点はすべての固有値の絶対値が1より小さければ安定であり( \abs{\lambda_{j}} < 1)、そうでなければ不安定である。

連続系と離散系の関係

連続系と離散系の対応関係は以下のようになる。区別のため、それぞれの時間発展行列を \mathbf{A}_{c}, \ \mathbf{A}_{d}とすると、離散時間が t_{k} = k \Delta tで与えられることを思い出せば、


\begin{align}
\begin{aligned}
\mathbf{x} \paren{(k+1) \Delta t} &= e^{\mathbf{A}_{c} \paren{(k+1)\Delta t}} \mathbf{x} \paren{0} \\
 &= e^{\mathbf{A}_{c} \Delta t } e^{\mathbf{A}_{c} \paren{k \Delta t}} \mathbf{x} \paren{0} \\
 &= e^{\mathbf{A}_{c} \Delta t } \mathbf{x} \paren{k \Delta t} 
\end{aligned}
\end{align}

であるから、


\begin{align}
\mathbf{A}_{d} = e^{\mathbf{A}_{c} \Delta t}
\end{align}

となることが分かる。また、連続系の対角化行列を \mathbf{U}_{c}とすると、


\begin{align}
\begin{aligned}
\mathbf{A}_{d} \mathbf{U}_{c} &= e^{\mathbf{A}_{c} \Delta t} \mathbf{U}_{c} \\
&= \paren{ \mathbf{U}_{c} e^{\boldsymbol{\Lambda}_{c} \Delta t} \mathbf{U}_{c}^{-1} } \mathbf{U}_{c} \\
&= \mathbf{U}_{c} e^{\boldsymbol{\Lambda}_{c} \Delta t} 
\end{aligned}
\end{align}

となることから、連続系と離散系の固有ベクトルは等しく( \mathbf{U}_{d} = \mathbf{U}_{c})、離散系の固有値 \boldsymbol{\Lambda}_{d} = e^{\boldsymbol{\Lambda}_{c} \Delta t} で与えられることが分かる(逆にいえば、離散系の固有値から連続系の固有値 \lambda_{c} = \frac{\ln \lambda_{d}}{\Delta t}で計算できる)。

また、 \abs{e^{\lambda \Delta t}} = \abs{e^{\paren{\alpha + i \omega}\Delta t} } = \abs{e^{\alpha \Delta t}} \abs{e^{i \omega \Delta t}} より、連続系が安定( \alpha < 0)なら対応する離散系も安定であることがわかる( \Delta t > 0 ,\ \abs{e^{i \omega \Delta t}} = 1 なので \abs{e^{\alpha \Delta t} } < 1となる)。

時間発展行列の計算

ともあれ、次なる問題はこの行列$\mathbf{A}$をどう求めるかだ。最初に考えるのは$\mathbf{X}$の逆行列を使って$\mathbf{A} = \mathbf{X}' \mathbf{X}^{-1}$とすればいいじゃんというものだが、そもそもの作り方からしてデータ行列は基本的に非正方行列であるし、かりに正方行列であったとしても正則性(逆行列の存在)は保証されない。そこで、近似的に上式が成り立つような行列を次式の条件から求めるという方針になる(数式の書き方これであってんのかな……?)。


\begin{align}
\mathbf{A} \triangleq \operatorname*{argmin}_{\mathbf{A}} \| \mathbf{X}' - \mathbf{A} \mathbf{X} \|_{F}
\end{align}

ここで、$\| \cdot \|_F$はFrobeniusノルム(行列ノルム)である。上式を満たす行列はMoore-Penroseの擬似逆行列$\mathbf{X}^{\dagger}$を使って次式のように得られる。


\begin{align}
\mathbf{A} \approx \mathbf{X}' \mathbf{X}^{\dagger}
\end{align}

$\mathbf{X}$の特異値分解が求まっていれば、擬似逆行列は次式で計算できる。


\begin{align}
&\mathbf{X} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{*} \\
\Longrightarrow \quad & \mathbf{A} = \mathbf{X}' \mathbf{V}_{r} \boldsymbol{\Sigma}_{r}^{-1} \mathbf{U}_{r}^{*}
\end{align}

ただし、$^*$は複素共役転置である(DMDの計算では固有値固有ベクトル複素数になる)。また、下付きの r特異値分解の時点でランク$r$の低ランク近似を行うことを意味する(特異値の大きな順に r個の要素まで残してそれ以降は打ち切る)。

後は行列$\mathbf{A}$の固有値固有ベクトルを求めてやれば時間発展の周波数と固有モードを得ることができる。

ただし、$\mathbf{A} \in \mathbb{C}^{n \times n}$より、行列$\mathbf{A}$のサイズが非常に大きいため、そのまま固有値問題を解くことは得策ではない。そのため、実際の計算では次式のようにPOD基底  \mathbf{U}_{r}に射影した行列$\tilde{\mathbf{A}}$を求め、こちらに対して固有値問題を解くことで計算コストを小さくする(相似変換であるため、固有値は同じになる)。


\begin{align}
&\tilde{\mathbf{A}} = \mathbf{U}_{r}^{*} \mathbf{A} \mathbf{U}_{r} = \mathbf{U}_{r}^{*} \mathbf{X}' \mathbf{V}_{r} \boldsymbol{\Sigma}_{r}^{-1} \\
\Longrightarrow \quad & \tilde{\mathbf{A}} \mathbf{W} = \mathbf{W} \boldsymbol{\Lambda}
\end{align}

したがって、元の行列 \mathbf{A}の固有モード(DMDモード)が次式のように得られる。


\begin{align}
\begin{aligned}
\mathbf{A} \boldsymbol{\Phi} &= \boldsymbol{\Phi} \boldsymbol{\Lambda} \\
\boldsymbol{\Phi} &= \mathbf{X}' \mathbf{V}_{r} \boldsymbol{\Sigma}_{r}^{-1} \mathbf{W}
\end{aligned}
\end{align}

DMDモードを上式で求める手法をexact DMDと呼ぶ。一方、元々の相似変換の関係からDMDモードを \boldsymbol{\Phi} = \mathbf{U}_{r} \mathbf{W} から求める手法はprojected DMDと呼ぶ(DMDが提案されたオリジナルの論文ではこちらを使っている)。exact DMDは演算に \mathbf{U}_{r} = \mathbf{X} \mathbf{V}_{r} \boldsymbol{\Sigma}_{r}^{-1}ではなく、 \mathbf{U}'_{r} = \mathbf{X}' \mathbf{V}_{r} \boldsymbol{\Sigma}_{r}^{-1}を使っている点に注意。
Tu(2014)によれば、projected DMD  \mathbf{X}のデータしか使っておらず、もはや元の行列  \mathbf{A}固有ベクトルにはなっていないが、exact DMD  \mathbf{X}  \mathbf{X}'の両方のデータを使っており、  \mathbf{A}固有ベクトルになっている。そのため、exact DMDの方が厳密であり、自然な手法ということらしい(実際、今はexact DMDの方が主流っぽい)。

DMDモードを用いた予測

さて、固有モードが求まれば解の再構成が可能になる。離散力学系 \mathbf{x}_{k+1} = \mathbf{A} \mathbf{x}_{k}の解は固有分解 \mathbf{A}  = \boldsymbol{\Phi} \boldsymbol{\Lambda} \boldsymbol{\Phi}^{\dagger}を用いて次式で与えられるのだった( \boldsymbol{\Phi}が正則とは限らないので擬似逆行列になる)。


\begin{align}
\mathbf{x}_{k} = \boldsymbol{\Phi} e^{\boldsymbol{\Lambda}^{k}} \boldsymbol{\Phi}^{\dagger} \mathbf{x}_{0}
\end{align}

離散系と連続系の対応関係から、 \boldsymbol{\Omega} = \ln \boldsymbol{\Lambda} \left/ \Delta t \right.とおけば、上式を任意の時間に対する解として表現することができる(離散的なデータから連続な時間の解が求まる!)。


\begin{align}
\mathbf{x}(t) = \boldsymbol{\Phi} e^{\boldsymbol{\Omega t}} \mathbf{b}
\end{align}

ただし、 \mathbf{b} = \boldsymbol{\Phi}^{\dagger} \mathbf{x}_{0}とおいた。これは次のようなDMDモードによるモード展開を表している点に注意する( \boldsymbol{\xi}は適当な空間座標)。


\begin{align}
\mathbf{x}(t) = \sum_{j=0}^{r} b_{j}(t) \varphi_{j}\paren{\boldsymbol{\xi}}  = \sum_{j=0}^{r} b_{j}(0) e^{\omega_{j} t} \varphi_{j}\paren{\boldsymbol{\xi}}
\end{align}

PODと異なり、こちらは時間発展を計算するだけなら展開係数の計算は要らないのが特徴だ。原理的には周期的な現象のデータを1周期分だけ取っておけば、それ以降の振る舞いを何周期分でも計算で予測することができる(とはいえ、実際には長時間予測には課題が多いようだ)。
また、DMDモードは直交基底ではないということもPODと異なる点だ。これは一見悪いことのように思えるが、直交性の条件に縛られないことによって柔軟な表現力があるともいえる(主成分分析と独立成分分析の関係に似ているかも)。ROMシミュレーションを考える場合は双対基底(左固有ベクトル)を使うことで展開係数も簡易に計算できる(Wang&Wei(2017))。

計算手順

DMDアルゴリズムをまとめると、以下のようになる。

  1. 時系列データからデータ行列$\mathbf{X}$と$\mathbf{X}'$を作成
  2. $\mathbf{X}$の特異値分解 \mathbf{X} = \mathbf{U}_{r} \boldsymbol{\Sigma}_{r} \mathbf{V}_{r}^{*}を計算し、低ランク近似を行う
  3. POD基底に射影した時間発展行列を \tilde{\mathbf{A}} = \mathbf{U}_{r}^{*} \mathbf{X}' \mathbf{V}_{r} \boldsymbol{\Sigma}_{r}^{-1} で求める
  4.  \tilde{\mathbf{A}}固有値問題 \tilde{\mathbf{A}} \mathbf{W} = \mathbf{W} \boldsymbol{\Lambda}を解く
  5. 元の行列の固有モード(DMDモード)を \boldsymbol{\Phi} = \mathbf{X}' \mathbf{V}_{r} \boldsymbol{\Sigma}_{r}^{-1} \mathbf{W}で計算
  6. 任意の時刻に対する解は \mathbf{x}(t) = \boldsymbol{\Phi} e^{\boldsymbol{\Omega t}} \boldsymbol{\Phi}^{\dagger} \mathbf{x}_{0}で計算できる( \boldsymbol{\Omega} = \ln \boldsymbol{\Lambda} \left/ \Delta t \right.

まとめ

DMDの概要についてまとめた。本当ならpythonでテスト計算もやるべきなんだけど、力尽きたので別の記事でやることにする(実行例は上述のブログを見てください)。
個人的には離散系と連続系の対応関係について確認できてよかった(論文とか読んでて固有値のlogを取ってるのを見て訳がわからなかったので)。
やっぱり力学系とか現代制御とかをもっと勉強しないとかなぁ。

特異値分解を用いたPODの定式化

注意

以下の記述は私的な備忘録であり、内容の正確さに関しては一切保証いたしませんので、悪しからず。

概要

前回、固有直交分解(POD)をやるときに共分散行列を作って計算していたが、よくよく考えるとそれは元のデータ行列を特異値分解(SVD)するだけでよくね? このことは色んなところに書いてある内容だし、なぜ気づかなかったのか不思議なくらい当たり前のことではないか。
というわけで、特異値分解を用いたPODについてまとめる。また、特異値分解の計算方法についても簡単に触れる。
参考サイト様のように主成分分析に関する優れた解説がネットにあふれているので、分からない向きはそちらを見てつかあさい。

話の流れ

最初に要点だけ書き出すと、PODは固有値分解じゃなくて特異値分解でいいよねというのは以下のような話である。

  1. PODモード$\boldsymbol{\Phi}$はデータ行列$\mathbf{X}$の共分散行列$\mathbf{X}\mathbf{X}^{T}$の正規直交基底(固有ベクトル)である
  2. $\mathbf{X}\mathbf{X}^{T}$は巨大な行列なので、実際には$\mathbf{X}^{T}\mathbf{X}$の基底$\mathbf{V}$を求め、それを元の基底に射影する($\boldsymbol{\Phi} = \mathbf{X} \mathbf{V} \boldsymbol{\Lambda}^{-1/2}$)
  3. 一方、$\mathbf{X}$の特異値分解 \mathbf{X} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^{T}の内、$\mathbf{U}$と$\mathbf{V}$はそれぞれ$\mathbf{X}\mathbf{X}^{T}$の$\mathbf{X}^{T}\mathbf{X}$の基底である
  4. したがって、求めるPODモードは \boldsymbol{\Phi} = \mathbf{U} = \mathbf{X} \mathbf{V} \boldsymbol{\Sigma}^{-1}である。これはPODにおける射影演算そのものである(\boldsymbol{\Lambda} = \boldsymbol{\Sigma}^{2}

以下、もう少し詳しく式展開の流れを追っていこう(ただし、数学的な厳密さはなしで)。

正方行列の対角化

まずは正方行列の対角化の話からおさらいする。
正方行列$\mathbf{A} \in \mathbb{R}^{n \times n}$が半単純(semisimple)であれば、$\mathbf{A}$を次のような行列の積に分解できる。


\begin{align}
\label{eq:rev}
\mathbf{A} = \mathbf{U} \boldsymbol{\Lambda} \mathbf{U}^{-1}
\end{align}

ここで、$\mathbf{U}$は$\mathbf{A}$の固有ベクトルから定まる行列であり、$\boldsymbol{\Lambda}$は$\mathbf{A}$の固有値を対角項に持つ対角行列である。
行列が半単純であるとは線形独立な固有ベクトルが$n$本取れるという意味(単純行列は固有値がすべて異なる。半単純行列は固有値の重複があったとしても固有ベクトルはすべて異なる)であり、対角化可能と同値だ。変換行列$\mathbf{U}$は$\mathbf{A}$の固有ベクトルを並べて作られるが、細かいことをいえば、これは右固有ベクトル$\mathbf{u}$のことだ。


\begin{align}
\begin{aligned}
& \mathbf{A} \mathbf{u} = \lambda \mathbf{u} \\
\Longrightarrow \quad & \mathbf{A} \mathbf{U} = \mathbf{U} \boldsymbol{\Lambda}
\end{aligned}
\end{align}

入門的な線形代数の教科書だと触れてない印象があるのだが、右固有ベクトルと対を成す左固有ベクトル$\mathbf{v}$もある。


\begin{align}
\begin{aligned}
& \mathbf{v}^{T} \mathbf{A} = \lambda \mathbf{v}^{T} \\
\Longleftrightarrow \quad & \mathbf{A}^{T} \mathbf{v} = \lambda \mathbf{v}
\end{aligned}
\end{align}

すなわち、$\mathbf{A}$の左固有ベクトルは$\mathbf{A}^{T}$の右固有ベクトルである。最初の式で転置しているのは行ベクトルを列ベクトルとして統一的に表記するためである。また、$\mathbf{A}$と$\mathbf{A}^{T}$の固有値は同じなので、右固有ベクトルと左固有ベクトル固有値は一致する。

さて、右固有ベクトルで対角化ができるんだから左固有ベクトルでも同じことができる。


\begin{align}
\begin{aligned}
&  \mathbf{A}^{T} \mathbf{V} = \mathbf{V} \boldsymbol{\Lambda} \\
\Longleftrightarrow \quad & \mathbf{V}^{T} \mathbf{A} = \boldsymbol{\Lambda} \mathbf{V}^{T} \\
\Longleftrightarrow \quad & \mathbf{A} = \left( \mathbf{V}^{T} \right)^{-1} \boldsymbol{\Lambda} \mathbf{V}^{T}
\end{aligned}
\label{eq:lev}
\end{align}

 \eqref{eq:rev} \eqref{eq:lev}を見比べれば、次式が導かれる。


\begin{align}
\begin{aligned}
& \mathbf{A} = \mathbf{U} \boldsymbol{\Lambda} \mathbf{U}^{-1} = \left( \mathbf{V}^{T} \right)^{-1} \boldsymbol{\Lambda} \mathbf{V}^{T}  \\
\Longleftrightarrow \quad & \mathbf{U} = \left( \mathbf{V}^{T} \right)^{-1} \land \mathbf{V}^{T} = \mathbf{U}^{-1} \\
\Longleftrightarrow \quad & \mathbf{U}\mathbf{V}^{T} = \mathbf{V}^{T} \mathbf{U} = \mathbf{I}
\end{aligned}
\end{align}

ここで、$\mathbf{I}$は単位行列である。つまり、右固有ベクトルと左固有ベクトルを並べて作った行列は互いに逆行列の関係にあり、取りも直さず、これは右固有ベクトルと左固有ベクトルが直交することを意味する。こういうのを二重直交基底(bi-orthogonal basis)または双対基底(dual basis)というらしい。

というわけで、正方行列の対角化は、右固有ベクトルと左固有ベクトルを使って次式のように表される。

\begin{align}
\label{eq:lrd}
& \mathbf{A} = \mathbf{U} \boldsymbol{\Lambda} \mathbf{V}^{T}
\end{align}

普通、わざわざ式 \eqref{eq:lrd}のような書き方はせず、対角化の表記は右固有ベクトルに基づいたもので統一されている。通常、固有ベクトルといったときには右固有ベクトルを指すので、双対基底の話を持ち出してまで上の表記を使うのは話がややこしくなるからだろう。しかし、行列は積が非可換になるように、右からの作用と左からの作用が区別される分野である(その方が多い?)。そういった意味で、色々学んだ後に対角化の定式化を見直すと、式 \eqref{eq:lrd}のような書き方の方が自然なようにも思われる。

なぜこういう書き方をしたかというと、以下の特異値分解の話と対応するからだ。特異値分解の方が一般的な話なので、固有値分解の表記をそれに合わせておくと話がすっきりすると思う。

最後に特別な場合として、$\mathbf{A}$が対称行列、すなわち$\mathbf{A} = \mathbf{A}^{T}$である場合を考える。応用上重要であり、かつ数学的に話が簡単になるので、一般の教科書にはこっちのことしか載ってなかったりするかもしれない。
このとき、左固有ベクトルは元の行列の転置行列の右固有ベクトルであったことを思い出せば、対称行列の場合は右固有ベクトルと左固有ベクトルが一致する($\mathbf{U} = \mathbf{V} \equiv \mathbf{P}$)。したがって、逆行列の関係式は$\mathbf{P}\mathbf{P}^{T} = \mathbf{P}^{T} \mathbf{P} = \mathbf{I}$となり、これは$\mathbf{P}$が直交行列であることを意味する。これにより、対称行列の場合は次式のように直交行列によって対角化できる(対称行列は半単純であり、常に対角化可能)。


\begin{align}
& \mathbf{A} = \mathbf{P} \boldsymbol{\Lambda} \mathbf{P}^{T}
\end{align}

非正方行列の特異値分解

特異値分解(Singular Value Decomposition, SVD)は正方行列の対角化の一般化であり、非正方行列に対しても固有値分解と同様の分解が可能である。
詳しい証明などはネット上の適当な文献を当たってほしい。

非正方行列$\mathbf{X} \in \mathbb{R}^{n \times m}$を考える。ここでは$ n > m $の場合を考えるが、$ n < m $の場合でも同様である。当然、$ n = m $なら正方行列となる。
さて、非正方行列はそのまま対角化を行うことはできない。というのも、$\mathbf{X}$と積を取ることができるベクトルは\mathbf{w} \in \mathbb{R}^{m}だが、 \mathbf{X} \mathbf{w} \in \mathbb{R}^{n}であり、左辺と右辺でベクトルの次元が合わないのだ。非正方行列では正方行列の固有値問題を素直に適用するはできないのである。

そこで、非正方行列から$\mathbf{X}^{T} \mathbf{X} \in \mathbb{R}^{m \times m}$という対称行列を作る(これが対称行列であることは転置行列の性質からすぐに確かめられる)。すると、対称行列は半正定値であるため、固有値は非負実数であり、かつ常に対角化可能である。そこで、直交行列$\mathbf{V}\in \mathbb{R}^{m \times m}$を使って次式のように書ける。


\begin{align}
\begin{aligned}
%& \mathbf{X}\mathbf{X}^{T} = \mathbf{U} \begin{pmatrix} \boldsymbol{ \Lambda} & \\ & \mathbf{O} \end{pmatrix} \mathbf{U}^{T} \\
& \paren{ \mathbf{X}^{T} \mathbf{X} } \mathbf{V} = \mathbf{V} \boldsymbol{\Lambda} \\
\Longleftrightarrow \quad & \paren{ \mathbf{X} \mathbf{X}^{T} } \mathbf{X}  \mathbf{V} = \paren{ \mathbf{X} \mathbf{V} } \boldsymbol{\Lambda} 
\end{aligned}
\end{align}

すなわち、$\mathbf{X}$から作られるもう1つの対称行列$\mathbf{X} \mathbf{X}^{T}$の固有ベクトルは$\mathbf{X} \mathbf{V}$であり、両者の0でない固有値は一致する。ここで、この固有ベクトルのノルムを計算すると、直交行列のノルムは1($\| \mathbf{V} \| = 1$)なので、


\begin{align}
\begin{aligned}
\paren{ \mathbf{X} \mathbf{V} }^{T} \paren{ \mathbf{X} \mathbf{V} } = \mathbf{V}^{T} \paren{ \mathbf{X}^{T} \mathbf{X} } \mathbf{V} =  \mathbf{V}^{T} \paren{ \mathbf{V} \boldsymbol{\Lambda} \mathbf{V}^{T} } \mathbf{V} = \boldsymbol{\Lambda}
\end{aligned}
\end{align}

となる。ベクトルノルムは$\| x \| = \mathbf{x}^{T} \mathbf{x} = |\mathbf{x}|^{2}$であるから、新しい固有ベクトルの大きさは \sqrt{\lambda}_k \equiv \sigma_kである($\lambda_k \geq 0$なので平方根が取れる)。 \sigma_kを$\mathbf{X}$の特異値と呼ぶ。
よって、$\mathbf{X} \mathbf{X}^{T}$の正規化された固有ベクトルを次式で定める。


\begin{align}
\mathbf{U} = \mathbf{X} \mathbf{V} \boldsymbol{\Sigma}^{-1}
\end{align}

ここで、 \boldsymbol{\Sigma}は特異値の対角行列である。上式を変形すれば、$\mathbf{X}$の特異値分解が得られる。


\begin{align}
\mathbf{X} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^{T}
\end{align}

ここで、$\mathbf{U}$を左特異ベクトル、$\mathbf{V}$を右特異ベクトルと呼ぶ。上式と式 \eqref{eq:lrd}は記号の定義が異なるため同じ式ではないが、正方対称行列の場合は特異値と特異ベクトルはそれぞれ固有値固有ベクトルに一致する。

というわけで、PODモードはデータ行列を特異値分解したときの左特異ベクトルである(この場合、累積寄与率は特異値で考える)。

特異値分解の計算方法

上述の通り、特異値分解は元の行列から対称行列を作り、その固有値固有ベクトルから求められる。しかし、この対称行列は元の行列の行数と列数の大きい方のサイズを持つので、どちらか一方が大きければ固有値問題のサイズが大きくなり、計算コストが大きい(精度もよくないらしい)。そのため、実際の特異値分解の計算ではこうした対称行列を作らないで計算を行うようだ。

基本的な流れは次のような感じ(参考文献の要約)。

  1. 非正方行列$\mathbf{X}$を$\mathbf{X} = \mathbf{Q}\mathbf{R}$とQR分解する。$\mathbf{Q}$が非正方列直交行列、$\mathbf{R}$が正方上三角行列になるため、$\mathbf{R}$の特異値分解が求まれば、 \mathbf{X} = \paren{\mathbf{Q}\mathbf{U}} \boldsymbol{\Sigma} \mathbf{V}^{T}から計算できる。そのため、正方行列の特異値分解が求められればよい。
  2. 対称行列を作らないで特異値分解を求めるには、正方行列を上二重対角行列に変換し、その特異値分解を求めてから逆変換すればよい。
  3. 二重対角化は適当な直交行列を用いたHouseholder変換で計算する。
  4. 二重対角行列の特異値分解は行列の対角化に帰着され、対角化に用いた直交行列と得られた対角行列が特異値分解の近似値となる。手法としてはQR法と分割統治法の2つがよく用いられる。QR法は計算コストが高く、分割統治法に比べて計算速度が遅いが、入力に寄らず得られる特異ベクトルの直交性が高い。分割統治法は並列化向きだがメモリ容量が大きい。
  5. 最後に二重対角化に用いた直交行列と対角化に用いた直交行列の積から求める特異ベクトルに変換する。

たとえば、scipyでSVDを計算する関数(scipy.linalg.svd)は、デフォルトでLAPACKのgesddルーチンを使う(分割統治法)。一方、オプションでgesvdルーチン(QR法)も使えるが、こちらはMatlabOctaveなどで使われているようだ。

その他の行列分解手法

SVDは広く応用されている重要な手法だが、計算量が比較的大きく($ n > m $のとき$O(nm^2)$らしい)、大規模な行列になってくると途端に計算時間がかかってしまう。そのため、大規模問題に対してはSVDを使わないように工夫することもなども行われているようだ。
詳しいアルゴリズムまでは調べられていないのだが、行列の低ランク近似の文脈において、演算量の低減やスパース性の保持などを目的としてSVDに代わる行列分解の方法が研究されているらしい。中でもCX分解(またはCUR分解?)は元の行列からいくつかの列をサンプリングして近似行列を構成することで相対的な誤差が小さい近似行列を生成する確率的手法らしく、興味深い。
POD/DMDの文脈でこういった手法がどの程度必要になってくるかは不明だが、データ行列が巨大になってくるとSVDを適用するのが難しくなってくるため、こういった新しい手法も研究されていくのかもしれない。データ行列そのものを低ランク近似してからPOD/DMDにかけるのか、それともPOD/DMDに含まれるSVDのプロセスで計算量を低減化させるのに用いられるのかは分からないが、データ量というのものは増えこそすれ減ることはないものだから、こういった視点でのアルゴリズムの改良は後々効いてくると思う。

まとめ

PODはSVDでやればいいよ、という話であった。当初考えていたより長くなってしまったが、復習は大事なのでよしとしよう。

固有直交分解に関する覚書

注意

以下の記述は私的な備忘録であり、内容の正確さに関しては一切保証いたしませんので、悪しからず。

はじめに

流体解析などの分野でデータ解析や縮約モデリングの手法として用いられる固有直交分解の概要について簡単にまとめる。
サンプルとして1次元の線形および非線形の拡散方程式を解き、ROMシミュレーションまで可能か試してみる。
他人が見てもわかりやすい記事にはなっていないと思うけど、何かの足しにでもしてもらえれば僥倖っす。

POD

固有直交分解(Proper Orthogonal Decomposition, POD)とは、時系列データから特徴的な固有モードを抽出する分析手法の一つであり、統計学の分野で主成分分析と呼ばれるものと数学的には等価である。また、求めたモードから支配的な成分のみを抜き出し、元のデータを少数のモードで再構築した縮約モデルを作成し、低コストなシミュレーションを行うこともできる(Reduced Order Model, ROM)。詳しくはこちらを参照のこと。また、主成分分析に関する優れた解説がネットにゴロゴロ落ちているのでそちらも御覧じろ。

ここでは特に流体解析で用いられるSnapshot PODのアルゴリズムについてまとめる。いま、 m個の時系列データをまとめたデータ行列$\mathbf{X}$を次のように定義する。


\begin{align}
\mathbf{X} = 
\begin{bmatrix}
\vert & \vert  &   & \vert  \\
\mathbf{x}_{0} & \mathbf{x}_{1} & \cdots & \mathbf{x}_{m-1} \\
\vert & \vert  &   & \vert 
\end{bmatrix}
\end{align}

ここで、$\mathbf{x}_j$は時刻$t_j$での$n$個の格子点における物理量(速度や圧力など)を縦に並べた列ベクトルである。一般に、流体解析における格子点の数は数百万以上であるのに対して、時系列データ(スナップショット)の数は数千から数万程度であり、$ n \gg m $である。すなわち、$\mathbf{X}$は$ n \times m $の縦長の行列になる($\mathbf{X} \in \mathbb{R}^{n \times m}$)。また、データ行列$\mathbf{X}$は元のデータから時間平均値$ \frac{1}{m} \sum_{m} \mathbf{x}_m $を事前に差し引いておく。

PODではデータの特徴的な構造(固有モード)を抽出するため、次の共分散行列$\mathbf{X} \mathbf{X}^T$の固有値問題を解く。これはデータの分散を最大化する方向を求める条件から導かれる(詳しくは(ry)。


\begin{align}
\label{eq:XXT}
\mathbf{X} \mathbf{X}^T \mathbf{u}_k = \lambda_k \mathbf{u}_k
\end{align}

ここで、$\lambda_k$と$\mathbf{u}_k$はそれぞれ$k$番目の固有値固有ベクトルである($\mathbf{u}_k \in \mathbb{R}^{n}$)。
ところで、この問題には$\mathbf{X} \mathbf{X}^T \in \mathbb{R}^{n \times n}$という巨大な密行列が現れ、これをそのまま解くことは計算コストが非常に大きい。そこで、Snapshot PODでは代わりに次の固有値問題を解く。


\begin{align}
\label{eq:XTX}
\mathbf{X}^T \mathbf{X} \mathbf{v}_k = \lambda_k \mathbf{v}_k
\end{align}

ここで、$\mathbf{X}^T \mathbf{X} \in \mathbb{R}^{m \times m}$、$\mathbf{v}_k \in \mathbb{R}^{m}$であり、元の問題に比べてサイズが小さくなっている点に注意する。そして、求める固有ベクトル=固有モードは次式から得られる。


\begin{align}
\boldsymbol{\phi}_k = \mathbf{X} \mathbf{v}_k \left/ \sqrt{\lambda_{k}} \right.
\end{align}

これは実際の計算では行列としてまとめて行う。


\begin{align}
\mathbf{\Phi} = \mathbf{X} \mathbf{V} \mathbf{\Lambda}^{-1/2}
\end{align}

ここで、$\mathbf{\Phi}$と$\mathbf{V}$はそれぞれ固有モードと固有ベクトルを並べて作った行列であり、$\mathbf{\Lambda}$は固有値の対角行列である。

以前引っかかったので、なぜ$\mathbf{X}^T \mathbf{X}$でも求める固有モードが得られるのかについてメモしておく。構成の仕方から明らかなように、$\mathbf{X} \mathbf{X}^T$と$\mathbf{X}^T \mathbf{X}$はどちらも半正定値対称行列であるため、固有値はすべて非負である($\lambda_{k} \geq 0$)。対称行列は常に対角化可能であるため、これらを固有ベクトルから定まる直交行列によって対角化し、式変形を行っていくと、両者の零でない固有値が一致すること、固有ベクトルが上式によって互いに変換されることがわかる(この辺は特異値分解に関する解説などを参照)。そこまで厳密な話をしなくても、式 \eqref{eq:XTX}に左から$\mathbf{X}$をかけて$\mathbf{u}_k = \mathbf{X} \mathbf{v}_k$と置き換えれば式 \eqref{eq:XXT}になることはすぐに確かめられる(逆も同じ)。

ROMシミュレーション

速度などの物理量を$\mathbf{u}(\mathbf{x},t)$とする。一度その固有モードが求まれば、元の物理量を固有モードの重ね合わせとして次のように近似することができる(Fourier級数展開のイメージ)。


\begin{align}
\mathbf{u}(\mathbf{x},t) \approx \mathbf{\bar{u}}(\mathbf{x}) + \sum_{j=1}^{r} a_{j}(t) \boldsymbol{\phi}_{j}(\mathbf{x}) 
\end{align}

ここで、$\mathbf{\bar{u}}(\mathbf{x})$は時間平均、$a_{j}(t)$は各時刻におけるモードの大きさを表す展開係数である。$r$は採用するモードの数であり、通常、元のデータの次元の数より小さな数で済む場合が多い。元のデータが偏微分方程式に従うとすると、モード展開した式を元の式に代入することで、展開係数に関する常微分方程式が得られ、その時間発展を計算するだけでシミュレーションが可能になる(理想的には)。こうして作ったモデルは縮約モデルや低次元モデルなどと呼ばれる。

これの何が嬉しいかというと、偏微分方程式の計算よりも常微分方程式の方が計算コストが低いので、多数のパラメータについて何度もシミュレーションを行ったり、制御問題などに関連してリアルタイムで計算を行う必要がある場合などに有利になるのだ。
まぁ、実際にはそうそううまくいく例ばかりでもないみたいだが、何とかしてうまくいくようにと頭のいい人たちが日夜頑張ってくれているようだ。

実行例

線形拡散方程式

というわけで、練習として線形の拡散方程式の固有モードを求め、それを使ったROMシミュレーションをやってみようと思う。おなじみ拡散方程式は以下の通り。


\begin{align}
\pdif{u}{t} = \nu \pddif{u}{x}
\end{align}

ここで、$\nu = 0.01$は拡散係数である。これを解くプログラムはこんな感じ。今回は微分の計算には擬スペクトル法ではなくてnumpy.gradientを使っている。初期条件はみんな大好きGauss分布です。

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from numpy import dot
from numpy import pi, cos, sqrt
from scipy.linalg import eigh
from scipy.integrate import solve_ivp

# データの生成
L = 2.0 # 空間の長さ
T = 10 # 計算時間
N = 256 # 空間分割数
dt =1.0e-2 # 時間刻み
x = np.linspace(0, L, N) #空間
t = np.arange(0, T, dt) # 時間

# パラメータと初期条件
dd = 0.01
mu = 1.0
sigma = 0.2
u0 = exp(-(x-mu)**2 / sigma**2)

# 線形拡散方程式
def lde(t,u) : 
  ux = np.gradient(u,x)
  rhs = np.gradient(ux,x)
  return dd*rhs

# 数値積分の実行
sol = solve_ivp(lde,[0,T], u0, method='Radau', t_eval=t)
u = sol.y

結果は普通に求まる。

f:id:iqujack-lequina:20180505172545p:plain
拡散方程式の時間発展

カラーマップにするとこんな感じ。meshgridなデータじゃないのでコンターとかを描くには一工夫要りそうだったので今回はお手軽なimshowだけ(軸ラベルつけ忘れてるが、横軸が時間で縦棒が空間です)。

f:id:iqujack-lequina:20180505172533p:plain
拡散方程式のカラーマップ

解が求まったらデータ行列を作り、固有値問題を解く。今回は時間方向の方がデータ点数が大きいため、逆に行列のサイズが大きくなるように計算してしまっているが、まぁ気にしない(修正は簡単だ)。共分散行列は対称行列なので、固有値問題のソルバーはscipy.linalg.eighを使う。また、得られた固有値固有ベクトルは降順になるようにしておく。

# データ行列
u_ave = np.average(u, axis=1) # 列方向の平均
D = u - u_ave.reshape(len(u_ave), 1) # 時間平均を差し引く

# 固有値問題
R = (D.T).dot(D)
val, vec = eigh(R) # R is symmetric

# eighの戻り値は昇順なので逆順にして降順にする
val = val[::-1] 
vec = vec[:, ::-1]

f:id:iqujack-lequina:20180505172527p:plain
共分散行列の固有値

固有値をプロットしてみると、第1モードがかなり大きく、第3モード以降は0に近い値になっている。したがって、第1モードだけでも解の大まかな特徴を表現できていることが期待できる。いくつのモードを採用するのかという問題は、実は難しい問題だったりするらしいのだが、最も単純には、次式の累積寄与率が1に近いかどうかを見ながら判断する。


\begin{align}
\sum_{k=1}^{r} \lambda_{k} \left/ \sum_{k}^{m} \lambda_{k} \right. \approx 1
\end{align}

上式は全体のモードの内、いくつのモードが支配的なのかを見る指標であり、情報量やエネルギーの大きさを表している。これをプロットしてみると、第1モードですでに累積寄与率は95%になっており、第3モードまでで99.9%を超えていることがわかった。そこで、採用するモード数は$r = 3$とする。

f:id:iqujack-lequina:20180505172520p:plain
累積寄与率

# 固有モード
r = 3
vn = vec[:,:r]/sqrt(val[:r])
phi = D.dot(vn)

固有モードのプロット。この形は時間平均を差し引いたものである点に注意。時間平均を足せばモード1だけでもGauss分布っぽい形になる。

f:id:iqujack-lequina:20180505172540p:plain
固有モード

さて、固有モードがわかれば、今度はROMシミュレーションができる。拡散方程式にモード展開の式を代入し、モードとの内積を取って整理すれば、展開係数に関する常微分方程式が得られる(各モードは互いに直交しているため。Fourier級数の展開係数を求める手続きと同じ)。


\begin{align}
\dif{}{t} \mathbf{a} = \nu \paren{  \mathbf{\Phi}^T \mathbf{\bar{u}}_{xx} +  \mathbf{\Phi}^T \mathbf{\Phi}_{xx} \mathbf{a} }
\end{align}

ここで、下付き添字の$xx$は$x$の2回微分を表し、上式はベクトル表記を用いている。また、展開係数の初期値は$\mathbf{a}(0) = \mathbf{\Phi}^T \paren{ \mathbf{u}(\mathbf{x}, 0) - \mathbf{\bar{u}} }$から求められる。これを計算するのが以下のコードである。

# ROMシミュレーション

# 初期値
a0 = (u0 - u_ave).dot(phi)

# 平均値の勾配と分散
uax = np.gradient(u_ave,x)
uaxx = np.gradient(uax,x)

# 固有モードの勾配と分散
phix = np.gradient(phi,x, axis=0)
phixx = np.gradient(phix,x, axis=0)

def lde_rom(t,a) : 
  rhs1 = dot(phi.T, uaxx)
  rhs2 = dot((phi.T).dot(phixx), a)
  return dd*(rhs1 + rhs2)

sol_a = solve_ivp(lde_rom,[0,T], a0, method='Radau', t_eval=t)
a =sol_a.y
u_rom = u_ave.reshape(len(u_ave),1) + dot(phi, a)

展開係数の時間発展はこうなっている。モード1の係数が指数関数的に減衰してるのがそれっぽい。

f:id:iqujack-lequina:20180505180430p:plain
展開係数の時間発展

重ね合わせて結果を比較してみると、プロットはほとんど重なっている(計算初期の方が誤差が大きいようだ)。一応誤差、というか差分($\mathbf{u} - \mathbf{u}_{\mathrm{ROM}}$)のカラーマップも載せておく。

f:id:iqujack-lequina:20180505172543p:plain
元の結果とROMの結果の比較

f:id:iqujack-lequina:20180505172523p:plain
誤差のカラーマップ

というわけで、拡散方程式の解は3つのモードで表すことができ、ROMシミュレーションも元の結果をよく再現することがわかった。

非線形拡散方程式(多孔質媒体方程式)

固有モード展開が重ね合わせ(線形結合)で表される以上、元の問題が線形ならうまくいくのは当たり前のように思える。そこで、問題を非線形に変えてみて試してみた。例題は次の多孔質媒体方程式と呼ばれる非線形の拡散方程式である(他にもBarenblatt方程式とか色々名前があるらしい)。名前の通り、多孔質媒体中のガスの拡散などを表す方程式のようだ。


\begin{align}
\pdif{u}{t} = \pdif{}{x} \paren{ u^p \pdif{u}{x} }
\end{align}

ここで、$p$はべき指数であり、非線形の拡散係数の大きさを表す。拡散が濃度のべき乗に比例するので、濃度が大きいところは拡散しやすく、逆に濃度が小さいところはほとんど拡散しない、という感じになるはず。

計算のコードはさっきとほぼ同じなので、solve_ivpに渡す右辺関数の定義だけ載せる(まぁこれも見た目通りに入力すればいいだけだ)。

# パラメータ
pw = 3

# PM方程式
def pme(t,u) : 
  ux = np.gradient(u,x)
  rhs = np.gradient(u**pw*ux,x)
  return rhs

結果を見てみると、確かにそんな傾向がありそうだ。いわゆる「突っ立ち」のようなものが見えて、自発的に界面ができているように見えるのが面白いかも。

f:id:iqujack-lequina:20180505180315p:plain
多孔質媒体方程式の時間発展

f:id:iqujack-lequina:20180505180305p:plain
多孔質媒体方程式のカラーマップ

固有値をプロットしてみると、今度はある程度なめらかに分布している模様。当然、累積寄与率もある程度モードを取らないと大きくならない。固有モードは最初の3つだけ示す。

f:id:iqujack-lequina:20180505180302p:plain
共分散行列の固有値

f:id:iqujack-lequina:20180505182308p:plain
累積寄与率

f:id:iqujack-lequina:20180505180308p:plain
固有モード

とりあえず、モード数は$r = 20$としてみる。ROMシミュレーションのコードも変更点だけ載せる。この場合、非線形項はどうにもならないので、そのまま計算する(内積と要素積の使い分けに注意)。

# 平均値の勾配
uax = np.gradient(u_ave,x)

# 固有モードの勾配
phix = np.gradient(phi,x, axis=0)

def pme_rom(t,a) : 
  ux = uax + phix.dot(a)
  upw = (u_ave + phi.dot(a))**pw
  nl = upw * ux
  nlx = np.gradient(nl,x, axis=0)
  return (phi.T).dot(nlx)

展開係数の時間発展を5個のモードまでプロットした。モード1の傾向は拡散方程式と同様だが、今度は高次のモードの正負を色々調整しないとつらい感じなのかな。

f:id:iqujack-lequina:20180505180434p:plain
展開係数の時間発展

結果を見てみると、やはり最初の方が誤差が大きいようだ。時間平均を元に計算するから、拡散方程式だと初期分布を再現しにくいのは当然か。誤差はやはり界面付近で大きくなっている模様。多孔質媒体方程式だと「界面」で勾配がきつくなっているから、そこもモード数を増やさないと再現できないってことなんだろう(矩形波の再現にモード数が必要なのと同じ話?)。

f:id:iqujack-lequina:20180505180313p:plain
元の結果とROMの結果の比較

f:id:iqujack-lequina:20180505180301p:plain
誤差のカラーマップ

まとめ

1次元の簡単な問題ではあるが、PODとROMをやってみた。本当はモードはそのままでパラメータ(拡散係数とか初期分布とか)を変えたときにはどうなるって話をやるべきなんだろうけど、まとめ作業でちょっと疲れたのでまたの機会にする(お許しを)。

最初はKdV方程式でやろうとしてたんだけど、再構成後の波形がノイジーすぎて断念した。正確には、初期値の再構成はよかったんだが、そこから勾配を計算するとかなりノイズが強くなってしまい(微分だしね)、そのままシミュレーションするのは躊躇われるレベルだったので早々に諦めた(分散項はいうにおよばず)。やっぱり元の結果がまずかったのかなぁ? ネットで見つかる文献などにはKdV方程式やBurgers方程式でもうまくいったという例があるようなので、自分でもできるかと思ったんだが。ともあれ、POD-ROMは一筋縄ではいかなさそうだと再確認しただけでもよしとしよう(……いいのか?)。

擬スペクトル法によるKdV方程式の数値計算

注意

以下の記述は私的な備忘録であり、内容の正確さに関しては一切保証いたしませんので、悪しからず。

概要

PODのテスト用に時系列データのサンプルが欲しかったので、擬スペクトル法を用いてKdV方程式を解いてみた(PODのサンプルとして適切な選択なのかはおいおい確認することにしよう……)。
実行環境は例によってColaboratoryでござい。
「擬スペクトル法ってみんな使ってるし簡単だよね!」と舐めてかかっていたら、何の報いか色々うまくいかずに引っかかってしまったので、つまずいたところなどをメモしておく。

KdV方程式

Korteweg-de Vries方程式は次式で与えられる非線形波動方程式である(係数の取り方によって異なる定義が複数ある)。

\begin{align}
\frac{\partial u}{\partial t} + \alpha u \frac{\partial u}{\partial x} + \beta \frac{\partial^3 u}{\partial x^3} = 0
\end{align}
ここで、$\alpha, \beta > 0$はパラメータである。今回は特に$\alpha = 1$、$\beta = 0.022^2$とおいて有名な論文の結果が得られるか確認する(Zabusky & Kruskal (1965))。また、領域$[0,L]$は$L=2$の周期領域とする。

擬スペクトル法

擬スペクトル法(pseudo spectral method)は偏微分方程式に対する数値計算法の一つで、以下のように離散化を行う。

  1. 空間について離散フーリエ変換(DFT)を計算し、波数空間で空間微分などを計算(波数をかけるだけ)
  2. 時間についてはそのまま離散化する

これにより、偏微分方程式が波数空間における常微分方程式に変換されるため、Runge-Kutta法などによって時間積分を行い、得られた結果を逆変換すれば求める解が得られる。差分化を行わないので数値拡散がなく、高精度の解を得られるが、非線形項の取扱が難しかったり(波数空間だと計算コストが高いので一度実空間に戻す)、任意の境界形状は表現できなかったりするらしい。
とはいえ、1次元問題であれば、DFTも使えるpythonでやるにはお手軽な手法なので、今回は擬スペクトル法でやってみる。

実行例

コードはこんな感じ。ライブラリのインポート部分は後述のように実行時には使ってない関数も含まれるので注意。

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from numpy import pi, sin, cos, exp, real, imag
from scipy.integrate import solve_ivp
from numpy.fft import fft, ifft, fftfreq, rfft, irfft
from scipy.fftpack import diff as psdiff

# 空間、時間、波数の設定
L = 2.0 # 空間の長さ
T = 10 # 計算時間
N = 256 # 空間分割数
m = N/2+1 # 波数空間のデータ数
dt =1.0e-2 # 時間刻み
x = np.linspace(0, L, N) #空間
t = np.arange(0, T, dt) # 時間
k = fftfreq(m)*m*2*pi/L #波数

# パラメータと初期条件
a = 1.0
b = 0.022**2
u0 = cos(pi*x)

# KdV方程式
def kdv(t,u) : 
  v = rfft(u)
  ux = irfft(1.0j*k*v) # gradient term
  conv = u*ux # convection term
  disp = irfft( (1.0j*k)**3*v ) # dispersion term
  return -a*conv -b*disp

# 数値積分の実行
sol = solve_ivp(kdv, [0, T], u0, method='Radau', t_eval=t)
u = sol.y # 解を代入

コードのコアの部分は右辺の関数の定義部分なわけだが、擬スペクトル法だと非常に簡潔に書けてしまうのが魅力。真面目に差分法のコードを書くとさすがにもう少し長くなるはず。また、今回は勾配の計算のみを波数空間で行う形にしているので、最終結果は実空間での値である。波数空間で時間積分して最後に逆変換という形でも問題ないはず。どの道、非線形項の計算のために一度実空間に戻さないといけないので、変換と逆変換は毎ステップ行うことになる(線形方程式ならこの手の処理は要らないけど)。

結果のアニメーション

せっかくなのでアニメーションを作ろうと思ったのだが、Colaboratory上では一個のセル辺りのデータ容量?の上限が決まっているようで(9M?)、フレーム数が増えると実行できないようだったので、アニメーションだけはローカルの環境で作成した。gifにしようかとも思ったんだが、容量オーバーなのかアップロードに失敗したので、無難にyoutube経由でリンクを貼る方法にした。


最初の1フレーム目にグラフを重ね書きしてしまう不具合が直せなかったのだけど、全体の挙動は見えるので今回はこれでよしとしよう。参考にしたサイト様やWikipediaの動画などと比べると時間進行の速さにずれがあるような気がするが……結果はあまり信用できないのかもしれない。波形の様子なんかは似ているので、大間違いとまではいかないと思うけど。

大雑把なコードの説明

今回使った、あるいは使用を検討した関数などについて、コメントがてらつまずいたところなどをメモしておく。

fftとrfft

numpyで離散フーリエ変換をやる関数は大きく分けて"fft"と"rfft"の2種類がある。コマンドリファレンスなどを見ると、両者の違いは波数空間での負の周波数成分を考慮するかどうからしい。そのため、係数に複素数が含まれるような場合には"fft"を使う。一方、入力が実数のみであるなら負の周波数を考慮しないrfftの方が便利、ということらしい。ただし、ここに一つ落とし穴があった。

擬スペクトル法では関数の勾配を計算するために、関数をフーリエ変換し、波数空間で$i k$を掛けて逆変換を行う。この波数は次のように計算できる。

N = 256 # 空間分割数
m = N/2+1 # 波数空間のデータ数
k = fftfreq(m)*m*2*pi/L #波数

"fftfreq"は離散フーリエ変換時のサンプリング周波数を返す関数であり、定義式上データ点の個数で割られているため、データ点の個数と$2\pi / L$をかけることで波数の配列が得られる。このとき、フーリエ変換後のデータ点の個数には注意が必要だ。fft関数は負の周波数も考慮するため、波数空間でのデータ点の個数は元の空間のデータ点の個数と等しい。一方、rfftは正の周波数のみを考慮するため、データ点の個数は元々の半分になる(正確には、元のデータ点数$n$が偶数なら$(n/2)+1$、奇数なら$(n+1)/2$になる)。したがって、波数のデータ点数とフーリエ変換後のデータ点数が合うように定義していないと勾配の計算時にエラーになる(後で気づいたが、事前に適当なフーリエ変換をしておいて、その長さを"len(v)"のように取得する方がスマートだな)。

solve_ivp

scioyのv1.0から?入ったらしいODEを解く新しい関数。最初は上述のfftを使わないといけないと思っており、計算途中に複素数が入ってくるため、複素数でも同じような使用感で使えるらしいこの関数を使ってみた。結果としてrfftを使うならodeintでええやんというオチになってしまったが、せっかくなので簡単に使い方をメモしておく。

solve_ivpの使い方はodeintとほとんど同じだが、細かいところが色々と変わっている。

from scipy.integrate import solve_ivp

sol = solve_ivp(kdv, [0, T], u0, method='Radau', t_eval=t)
u = sol.y

まず、地味にややこしいのが右辺関数の定義がodeintと逆で"f(t,y)"の順になる。また、関数の引数としてパラメータを追加するような使い方は基本できないらしい(ラムダ式を使えば云々という情報がネットにあったが、俺の場合は特に困らなかったのでスルー)。
積分時間の指定も変わっていて、基本的にはソルバー側が時間刻みを自動で決める方針なので、"[t_start, t_end]"みたいな形で積分する区間を始点と終点のタプルで与える。積分点を指定したい場合は"t_eval"オプションで与える。
初期値とスキームの指定はご覧の通り。また、時間刻みをソルバーが決める関係上、戻り値は解と時間の配列がセットで返ってくる("sol.y"、"sol.t"のようにしてアクセス)。

係数に複素数が含まれる場合でも使い方は同じだが(odeだと書き方が変わるっぽい?)、その場合は初期値を複素数型のデータとして与える必要がある(初期値が実数であったとしても、"0.0j"を足すみたいな感じで複素数型にしておく)。また、スキームによって複素数に対応しているものとそうでないものがあるので、場合によってはスキームを変更する必要が出てくる。

今回は使ってないが、"dense_output"をTrueにすると、結果を時間方向に補間した関数で返してくれるようだ(戻り値は"sol.sol"みたいな感じ)。補間関数として何が使われるかもスキームに依存しているので、使う際はその辺も確認しないといけないだろう。また、補間関数が個別の評価点を使うためなのか、評価点("t_eval")をこちらで指定するとエラーになる。

ちなみに、デフォルトのRunge-Kutta法('RK45')では、今回使った時間刻みの値だと結果にものすごい高周波ノイズが入ってしまったので、陰的Runge-Kutta法('Radau')にスキームを変えている("BDF"でも可。odeintと同じスキームにしたいなら"LSODA"とする)。時間刻みをかなり小さくとればRK45でも計算できるのかもしれないが、問題がstiffならそこは素直に陰的スキームを使う法が適切だろう。

psdiff

最初はscipyのcookbookをもとにpsdiff("scipy.fftpack.diff")を使おうとしていた。これは以下のような感じで擬スペクトル法をもとに勾配を簡単に計算してくれる関数だ。

from scipy.fftpack import diff as psdiff

u = cos(pi*x)
ux = psdiff(u, period=L)
uxxx = psdiff(u, period=L, order=3)

擬スペクトル法を使っているので、関数の周期を入力する必要がある(指定しない場合は$2\pi$)。微分階数はデフォルトで1、それ以上の場合は"order"オプションで指定する。

この関数、使いやすいしコードも見やすくなって擬スペクトル法の計算のためにあるような便利なものなんだが、試しに計算してみると特に高階の計算でなぜか区間の両端にノイズが入ることがあり、使い物にならなかった。正確には、これを使ってもKdV方程式の計算そのものはうまくいったんだけど、それはノイズが出る分散項の係数が小さく設定されているからじゃないのかという疑いが晴れず、一般のケースに適用するには不安が残るため、結局は使わなかったのだ。原因自体はよくわからず、何か使い方を間違えているのかもしれないが、定義通りに計算するとノイズも出ないので、今回は愚直な方法を選んだ次第だ。

まとめ

紆余曲折あったが、それらしい結果が得られたので今回はこれでよしとする。もちろん、真面目にやるなら厳密解とか文献ときちんと合っているのか比較したり、保存性のチェックなんかをしなきゃいけないんだろうけど、面倒だしそれが目的でもないのでやらない。

追記

確認したらpsdiffを使わない方法(定義通りの計算)でも両端にノイズ入ってましたわ……結果が同じなら簡単に書けるpsdiffの方がいいね。確認不足でパチこいてすんません。
でもそうだとすると、余計にシミュレーション結果が怪しくなってくるな。擬スペクトル法ならではの処置を施さないとやっぱりうまいかないんだろうか。だとすると、素人にはちょっと難しい問題だ。
試しにnumpy.gradientで計算したところ、よっぽどきれいな結果(勾配と分散)が得られた。ただし、やはり境界で値がおかしくなるけど(当然っちゃ当然か)。一々境界値修正するのも面倒だしなぁ。すぐにはどうにもなりそうにないので、不完全燃焼だが、ここいらで撤退することにする。

上記の結果はあまり鵜呑みにしないでくださいますよう。

Colaboratoryのscipyのバージョンを上げる

注意

以下の記述は私的な備忘録であり、内容の正確さに関しては一切保証いたしませんので、悪しからず。

問題

以前やったodeintによる数値計算だが、scipy 1.0のドキュメントいわく、"ode"系のコマンドは古いAPIだそうで、今は"solve_ivp"などの新しいAPIに置き換わっているようだ。
別にodeintでも大抵の場合は問題ないと思うが、方程式が複素数の係数を含む場合にはodeintが使えず、少し面倒なやり方をしなければならなそうだったので、新しいコマンドを使ってみるためにColaboratoryのscipyのバージョンを上げることにした。

手順

まずは現行バージョンの確認。

import scipy
scipy.version.full_version

ちなみにnumpyとかもこれで同じように確認できる。結果は

'0.19.1'

と返ってきた。

バージョンを上げるのは次のようにする。おなじみpipコマンドは新しいツールを追加する場合などにも使う。

!pip install --upgrade scipy

実行すると、以下のように旧バージョンがアンインストールされて新しいバージョンがインストールされる。

Looking in indexes: https://pypi.org/simple, https://legacy.pypi.org/simple
Collecting scipy
  Downloading https://files.pythonhosted.org/packages/9c/0b/5deb712a9ea5bb0a1de837d04ef7625c5f3ee44efe7ed0765ceda681d7f1/scipy-1.0.1-cp27-cp27mu-manylinux1_x86_64.whl (46.7MB)
    100% |████████████████████████████████| 46.7MB 599kB/s 
Requirement not upgraded as not directly required: numpy>=1.8.2 in /usr/local/lib/python2.7/dist-packages (from scipy) (1.14.2)
Installing collected packages: scipy
  Found existing installation: scipy 0.19.1
    Uninstalling scipy-0.19.1:
      Successfully uninstalled scipy-0.19.1
Successfully installed scipy-1.0.1

最後にランタイムを再起動する。これをやらないと読み込まれるバージョンが変わらない?ので注意。

再起動後にもう一度バージョン確認を行い、

'1.0.1'

となっていればおk。

試しに先程の新しいコマンドを読み込んでみる。

from scipy.integrate import solve_ivp

エラーが出なければ作業完了(元のバージョンでは「そんなコマンドねーよ」と怒られる)。

さーて、後はやりたい計算をやるだけなんだが……うまくいくかな(汗)。