固有直交分解に関する覚書
注意
以下の記述は私的な備忘録であり、内容の正確さに関しては一切保証いたしませんので、悪しからず。
はじめに
流体解析などの分野でデータ解析や縮約モデリングの手法として用いられる固有直交分解の概要について簡単にまとめる。
サンプルとして1次元の線形および非線形の拡散方程式を解き、ROMシミュレーションまで可能か試してみる。
他人が見てもわかりやすい記事にはなっていないと思うけど、何かの足しにでもしてもらえれば僥倖っす。
POD
固有直交分解（Proper Orthogonal Decomposition, POD）とは、時系列データから特徴的な固有モードを抽出する分析手法の一つであり、統計学の分野で主成分分析と呼ばれるものと数学的には等価である。また、求めたモードから支配的な成分のみを抜き出し、元のデータを少数のモードで再構築した縮約モデルを作成し、低コストなシミュレーションを行うこともできる（Reduced Order Model, ROM）。詳しくはこちらを参照のこと。また、主成分分析に関する優れた解説がネットにゴロゴロ落ちているのでそちらも御覧じろ。
ここでは特に流体解析で用いられるSnapshot PODのアルゴリズムについてまとめる。いま、個の時系列データをまとめたデータ行列$\mathbf{X}$を次のように定義する。
ここで、$\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）。
ここで、$\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では代わりに次の固有値問題を解く。
ここで、$\mathbf{X}^T \mathbf{X} \in \mathbb{R}^{m \times m}$、$\mathbf{v}_k \in \mathbb{R}^{m}$であり、元の問題に比べてサイズが小さくなっている点に注意する。そして、求める固有ベクトル＝固有モードは次式から得られる。
これは実際の計算では行列としてまとめて行う。
ここで、$\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$）。対称行列は常に対角化可能であるため、これらを固有ベクトルから定まる直交行列によって対角化し、式変形を行っていくと、両者の零でない固有値が一致すること、固有ベクトルが上式によって互いに変換されることがわかる（この辺は特異値分解に関する解説などを参照）。そこまで厳密な話をしなくても、式に左から$\mathbf{X}$をかけて$\mathbf{u}_k = \mathbf{X} \mathbf{v}_k$と置き換えれば式
になることはすぐに確かめられる（逆も同じ）。
ROMシミュレーション
速度などの物理量を$\mathbf{u}(\mathbf{x},t)$とする。一度その固有モードが求まれば、元の物理量を固有モードの重ね合わせとして次のように近似することができる（Fourier級数展開のイメージ）。
ここで、$\mathbf{\bar{u}}(\mathbf{x})$は時間平均、$a_{j}(t)$は各時刻におけるモードの大きさを表す展開係数である。$r$は採用するモードの数であり、通常、元のデータの次元の数より小さな数で済む場合が多い。元のデータが偏微分方程式に従うとすると、モード展開した式を元の式に代入することで、展開係数に関する常微分方程式が得られ、その時間発展を計算するだけでシミュレーションが可能になる（理想的には）。こうして作ったモデルは縮約モデルや低次元モデルなどと呼ばれる。
これの何が嬉しいかというと、偏微分方程式の計算よりも常微分方程式の方が計算コストが低いので、多数のパラメータについて何度もシミュレーションを行ったり、制御問題などに関連してリアルタイムで計算を行う必要がある場合などに有利になるのだ。
まぁ、実際にはそうそううまくいく例ばかりでもないみたいだが、何とかしてうまくいくようにと頭のいい人たちが日夜頑張ってくれているようだ。
実行例
線形拡散方程式
というわけで、練習として線形の拡散方程式の固有モードを求め、それを使ったROMシミュレーションをやってみようと思う。おなじみ拡散方程式は以下の通り。
ここで、$\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
結果は普通に求まる。
カラーマップにするとこんな感じ。meshgridなデータじゃないのでコンターとかを描くには一工夫要りそうだったので今回はお手軽なimshowだけ（軸ラベルつけ忘れてるが、横軸が時間で縦棒が空間です）。
解が求まったらデータ行列を作り、固有値問題を解く。今回は時間方向の方がデータ点数が大きいため、逆に行列のサイズが大きくなるように計算してしまっているが、まぁ気にしない（修正は簡単だ）。共分散行列は対称行列なので、固有値問題のソルバーは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]
固有値をプロットしてみると、第1モードがかなり大きく、第3モード以降は0に近い値になっている。したがって、第1モードだけでも解の大まかな特徴を表現できていることが期待できる。いくつのモードを採用するのかという問題は、実は難しい問題だったりするらしいのだが、最も単純には、次式の累積寄与率が1に近いかどうかを見ながら判断する。
上式は全体のモードの内、いくつのモードが支配的なのかを見る指標であり、情報量やエネルギーの大きさを表している。これをプロットしてみると、第1モードですでに累積寄与率は95%になっており、第3モードまでで99.9%を超えていることがわかった。そこで、採用するモード数は$r = 3$とする。
# 固有モード r = 3 vn = vec[:,:r]/sqrt(val[:r]) phi = D.dot(vn)
固有モードのプロット。この形は時間平均を差し引いたものである点に注意。時間平均を足せばモード1だけでもGauss分布っぽい形になる。
さて、固有モードがわかれば、今度はROMシミュレーションができる。拡散方程式にモード展開の式を代入し、モードとの内積を取って整理すれば、展開係数に関する常微分方程式が得られる（各モードは互いに直交しているため。Fourier級数の展開係数を求める手続きと同じ）。
ここで、下付き添字の$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の係数が指数関数的に減衰してるのがそれっぽい。
重ね合わせて結果を比較してみると、プロットはほとんど重なっている（計算初期の方が誤差が大きいようだ）。一応誤差、というか差分（$\mathbf{u} - \mathbf{u}_{\mathrm{ROM}}$）のカラーマップも載せておく。
というわけで、拡散方程式の解は3つのモードで表すことができ、ROMシミュレーションも元の結果をよく再現することがわかった。
非線形拡散方程式（多孔質媒体方程式）
固有モード展開が重ね合わせ（線形結合）で表される以上、元の問題が線形ならうまくいくのは当たり前のように思える。そこで、問題を非線形に変えてみて試してみた。例題は次の多孔質媒体方程式と呼ばれる非線形の拡散方程式である（他にもBarenblatt方程式とか色々名前があるらしい）。名前の通り、多孔質媒体中のガスの拡散などを表す方程式のようだ。
ここで、$p$はべき指数であり、非線形の拡散係数の大きさを表す。拡散が濃度のべき乗に比例するので、濃度が大きいところは拡散しやすく、逆に濃度が小さいところはほとんど拡散しない、という感じになるはず。
計算のコードはさっきとほぼ同じなので、solve_ivpに渡す右辺関数の定義だけ載せる（まぁこれも見た目通りに入力すればいいだけだ）。
# パラメータ pw = 3 # PM方程式 def pme(t,u) : ux = np.gradient(u,x) rhs = np.gradient(u**pw*ux,x) return rhs
結果を見てみると、確かにそんな傾向がありそうだ。いわゆる「突っ立ち」のようなものが見えて、自発的に界面ができているように見えるのが面白いかも。
固有値をプロットしてみると、今度はある程度なめらかに分布している模様。当然、累積寄与率もある程度モードを取らないと大きくならない。固有モードは最初の3つだけ示す。
とりあえず、モード数は$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の傾向は拡散方程式と同様だが、今度は高次のモードの正負を色々調整しないとつらい感じなのかな。
結果を見てみると、やはり最初の方が誤差が大きいようだ。時間平均を元に計算するから、拡散方程式だと初期分布を再現しにくいのは当然か。誤差はやはり界面付近で大きくなっている模様。多孔質媒体方程式だと「界面」で勾配がきつくなっているから、そこもモード数を増やさないと再現できないってことなんだろう（矩形波の再現にモード数が必要なのと同じ話？）。
まとめ
1次元の簡単な問題ではあるが、PODとROMをやってみた。本当はモードはそのままでパラメータ（拡散係数とか初期分布とか）を変えたときにはどうなるって話をやるべきなんだろうけど、まとめ作業でちょっと疲れたのでまたの機会にする（お許しを）。
最初はKdV方程式でやろうとしてたんだけど、再構成後の波形がノイジーすぎて断念した。正確には、初期値の再構成はよかったんだが、そこから勾配を計算するとかなりノイズが強くなってしまい（微分だしね）、そのままシミュレーションするのは躊躇われるレベルだったので早々に諦めた（分散項はいうにおよばず）。やっぱり元の結果がまずかったのかなぁ？ ネットで見つかる文献などにはKdV方程式やBurgers方程式でもうまくいったという例があるようなので、自分でもできるかと思ったんだが。ともあれ、POD-ROMは一筋縄ではいかなさそうだと再確認しただけでもよしとしよう（……いいのか？）。
