FOA based DMDに関するノート（アルゴリズム編）

はじめに

個人的に興味深いアルゴリズムだったので、自分の勉強のためにAnantharamu & Mahesh(2019)およびAnantharamu et al.(2019)のFOA based DMDについてまとめる。

要約を先に述べるなら、このアルゴリズムは「Gram-Schmidtの正規直交化法を用いてデータ行列のArnoldi分解とQR分解を並行して行うもの」といえる。この文章が理解できる人なら元論文を読むだけでおおよそ内容を把握できるだろう。残念ながら、自分はよくわからなかったので、以下長々と続くノートを作ってようやく理解した（？）次第である。

なお、以下の内容には間違いが含まれるかもしれないので、大いにご注意ください。

背景

DMDの概要

Dynamic Mode Decomposition（動的モード分解、以下DMD）は、複雑な時系列データから動的モードと呼ばれる特徴的な構造を抽出し、その振る舞いを調べる解析手法の一つである。最初に流体力学の分野で提案されたため（Schmid(2010), Rowley et al.(2009)）、主に流体シミュレーション（CFD）の結果を解釈するためなどに用いられている。一方、データ駆動型の手法であることから様々な現象に適用できるため、現在では分野を問わず広範な領域で応用が研究されている。

DMDに関する詳しい解説などは以下の参考文献などを参照くだされ。

• 動的モード分解による多次元時系列解析：日本語の解説
• Modal Analysis of Fluid Flows: An Overview：網羅的な解説
• Dynamic Mode Decomposition (Overview)：youtubeの解説動画
• PyDMD：Pythonによるライブラリ

基本的な考え方

詳細は上記の参考文献に譲るが、一応ここでも基本的な考え方だけ簡単に説明しておく。DMDでは与えられた$m$個の時系列データを分割して次のような2つのデータ行列を作成する。

\begin{align}
\mathbf{X}_{1}^{m-1} &= \begin{bmatrix} \mathbf{x}_{1} & \mathbf{x}_{2} & \cdots & \mathbf{x}_{m-1} \end{bmatrix} \\
\mathbf{X}_{2}^{m} &= \begin{bmatrix} \mathbf{x}_{2} & \mathbf{x}_{3} & \cdots & \mathbf{x}_{m} \end{bmatrix} \\
\end{align}

ここで、$\mathbf{x}_{k} \in \mathbb{C}^{n}$は時刻$t_{k}$におけるデータベクトル（あるいはスナップショットともいう）であり、CFDでいえばある時刻の全格子点の物理量をまとめたものである。また、$\mathbf{X}_{i}^{j}$は$i$番目から$j$番目までのスナップショットを並べた行列である。このとき、DMDではデータの時間発展が次のような線形力学系で近似できると仮定する。

\mathbf{X}_{2}^{m} \approx \mathbf{A} \mathbf{X}_{1}^{m-1}

ここで、$\mathbf{A}$は系の時間発展を与える行列であり、上式は$\mathbf{x}_{k+1} = \mathbf{A} \mathbf{x}_{k}$という線形力学系をまとめて書いた形になっている。

行列$\mathbf{A}$は上式の残差ノルム$\| \mathbf{X}_{2}^{m} - \mathbf{A} \mathbf{X}_{1}^{m-1} \|_{2}$を最小化する行列として次式で求められる。¹

\mathbf{A} \approx  \mathbf{X}_{2}^{m} \mathbf{X}_{1}^{{m-1}^{+}}

ここで、$\mathbf{X}_{1}^{{m-1}^{+}}$は$\mathbf{X}_{1}^{m-1}$の擬似逆行列である。そして、目的である動的モードは$\mathbf{A}$の固有ベクトルとして求められる（実際のアルゴリズムでは計算コストを抑えるためにもう少し工夫をするが、その辺の話は省略）。

逐次処理の必要性

上述の通り、DMDのアルゴリズムは非常に単純であり、Pythonなどでスクリプトを組むこと自体は容易である。しかし、オリジナルのアルゴリズムにはメモリ消費がとても大きいという欠点がある。これは処理するデータをすべて一つの行列にまとめて計算を行うためである。特に、CFDではデータベクトルの長さは格子点数で決まり、その数は数百万から数千万以上にもなるため、必要なメモリが数百GBを超えることもあり得る。そのため、何らかの工夫をしなければ大規模データに対してDMDを適用することができない。

それでは具体的にどうするかというと、データ行列全体を一度に計算せず、いくつかのベクトルや部分行列に分割して逐次的に処理することでメモリを節約する。これにはDMD自体のアルゴリズムはもちろん、内部的な行列計算の実装方法などにも工夫が必要になる。このような逐次的な手法はいくつか提案されており、今回紹介するAnantharamu & Mahesh(2019)のFOA based DMDもそのような手法の1つである。以下、このアルゴリズムについて説明していくが、その準備としてまずはKrylov DMDとArnoldi法について述べる。

Krylov DMD

発案者の論文であるSchmid(2010)で提案されたDMDのアルゴリズムには以下の2つのタイプがある。

1. Krylov部分空間に基づく方法
2. 特異値分解（SVD）に基づく方法

以下、便宜的にそれぞれKrylov DMDとSVD based DMDと呼ぶことにする。上述の擬似逆行列を用いる方法がSVD based DMDであり、数値的安定性が高いという理由から、現在主流になっているのはこちらの方だと思われるが、FOA based DMDはもう一方のKrylov DMDに基づく方法であるため、まずはこれについて説明する。

上述の通り、DMDではデータの時間発展が$\mathbf{x}_{k+1} = \mathbf{A} \mathbf{x}_{k}$という形で表されると考える。このとき、データ行列は次のように書ける。

\begin{align}
\mathbf{X}_{1}^{m-1} &= \begin{bmatrix} \mathbf{x}_{1} & \mathbf{x}_{2} & \mathbf{x}_{3} & \cdots & \mathbf{x}_{m-1} \end{bmatrix} \\
&= \begin{bmatrix} \mathbf{x}_{1} & \mathbf{A}\mathbf{x}_{1} & \mathbf{A}^{2} \mathbf{x}_{1} & \cdots & \mathbf{A}^{m-2} \mathbf{x}_{1} \end{bmatrix} 
\end{align}

これより、$\mathbf{X}_{1}^{m-1}$の列空間

\mathcal{R}\left( \mathbf{X}_{1}^{m-1} \right) = \operatorname*{span} \left\{ \mathbf{x}_{1}, \mathbf{x}_{2}, \mathbf{x}_{3}, \cdots, \mathbf{x}_{m-1} \right\}

は、次のようなKrylov部分空間に等しい。

\mathcal{K}_{m-1} \left( \mathbf{A}, \mathbf{x}_{1} \right) = \operatorname*{span} \left\{ \mathbf{x}_{1}, \mathbf{A}\mathbf{x}_{1}, \mathbf{A}^{2} \mathbf{x}_{1}, \cdots, \mathbf{A}^{m-2} \mathbf{x}_{1} \right\}

したがって、最終時刻のスナップショット$\mathbf{x}_{m}$がこの部分空間に含まれる（$\mathbf{x}_{m} \in \mathcal{K}_{m-1} \left( \mathbf{A}, \mathbf{x}_{1} \right)$）とすれば²、$\mathbf{x}_{m}$は次式のようにそれまでの時刻のスナップショットの線形結合で表すことができる。

\mathbf{x}_{m} = c_{1} \mathbf{x}_{1} + c_{2} \mathbf{x}_{2} + \cdots + c_{m-1} \mathbf{x}_{m-1} + \mathbf{r}

ここで、$\mathbf{c} = \begin{bmatrix} c_{1} \ c_{2} \ \cdots \ c_{m-1} \end{bmatrix}^{\top}$は展開係数、$\mathbf{r}$は残差ベクトルである。行列でまとめて書けば、

\mathbf{x}_{m} = \mathbf{X}_{1}^{m-1} \mathbf{c}+ \mathbf{r}

となる。これより、改めて次式が成り立つ。

\mathbf{A} \mathbf{X}_{1}^{m-1} = \mathbf{X}_{2}^{m} = \mathbf{X}_{1}^{m-1} \mathbf{S} + \mathbf{r} \mathbf{e}_{m-1}^{\top}

ここで、$\mathbf{e}_{m-1} \in \mathbb{R}^{m-1}$は$m-1$番目の要素のみが$1$となる単位ベクトルである。また、$\mathbf{S}$は次のような同伴行列（Companion matrix）と呼ばれる形となる。

\mathbf{S} = 
\begin{bmatrix}
0 & & & & c_{1} \\
1 & 0 & & & c_{2} \\
& \ddots & \ddots & & \vdots \\
& & 1 & 0 & c_{m-2} \\
 & & & 1 & c_{m-1} \\
\end{bmatrix}

上式はそのまま基底の変換であるから、残差の部分を除いて$\mathbf{S} = \mathbf{X}_{1}^{{m-1}^{-1}} \mathbf{A} \mathbf{X}_{1}^{m-1}$という相似変換になっている。したがって、$\mathbf{S}$の固有値は$\mathbf{A}$の固有値と等しく、結局$\mathbf{S}$の固有値問題を解けばよいことになる。なお、$\mathbf{c}$は$\mathbf{X}_{1}^{m-1} \mathbf{c} = \mathbf{x}_{m}$という線形方程式を解けば求められるが、$\mathbf{X}_{1}^{m-1}$は一般に非正方行列（通常、$n \gg m$）であるため、最小二乗解を求めることになる。これは擬似逆行列を用いて$\mathbf{c} = \mathbf{X}_{1}^{{m-1}^{+}} \mathbf{x}_{m}$とするか、またはQR分解$\mathbf{X}_{1}^{m-1} = \mathbf{Q}\mathbf{R}$を用いて$\mathbf{R} \mathbf{c} = \mathbf{Q}^{*} \mathbf{x}_{m}$を解くことで求められる。

Arnoldi法

次に、上述のKrylov DMDに関連するアルゴリズムであるArnoldi法について説明する。というのも、FOA based DMDの"FOA"は"Full Orthogonalization Arnoldi"の略であり、Arnoldi法を修正したアルゴリズムが核となるためである。
なお、以下のArnoldi法に関する記述は主に数値線形代数の数理とHPCに基づく。

Ritz-Galerkin法とArnoldi分解

DMDで求めたいものは$\mathbf{A}$の固有対である。しかし、データ行列は通常ランク落ちしており、$\mathbf{A}$のサイズに等しいすべての固有対を求める必要はない。また、DMDを適用する動機もデータの振る舞いから支配的な要素のみを抽出することにあるので、その意味でも系を特徴付ける少数の固有対のみが求まれば十分である。

大規模行列から少数の固有対のみを求める固有値問題の近似解法として射影法（反復法）が知られている。これは元の問題より小さな部分空間に行列を射影し、必要な精度に応じて部分空間を拡大しながら近似固有対を計算する方法である。部分空間として何を利用するか、近似解を定めるためにどういった条件を課すかによって様々な方法があるが、ここではKrylov部分空間を用いたRitz-Galerkin法について述べる。

標準固有値問題

\mathbf{A} \mathbf{x} = \lambda \mathbf{x}

を考える。$\mathbf{A} \in \mathbb{C}^{n \times n}$に対し、元の次元より小さいk次元部分空間$\mathcal{V}_{k}$（$k \ll n$）の正規直交基底を$\{v_{i}\}_{i=1}^{k}$とし、これを並べた直交行列を$\mathbf{V}_{k} = \begin{bmatrix} \mathbf{v}_{1} \ \mathbf{v}_{2} \ \cdots \ \mathbf{v}_{k} \end{bmatrix}$とおく。また、固有値と固有ベクトルの近似解をそれぞれ$\theta, \mathbf{u}$とする。

$\mathcal{V}_{k}$の中で近似解を構成すると、$\mathbf{u} \in \mathcal{V}_{k}$より、適当な長さ$k$のベクトル$\mathbf{y}$に対して$\mathbf{u} = \mathbf{V}_{k} \mathbf{y}$と表される。また、近似固有対の残差を$\mathbf{r} = \mathbf{A} \mathbf{u} - \theta \mathbf{u}$として、これがなるべく小さくなるように次のような条件を課す（Galerkin条件）。

\mathbf{r} \perp \mathcal{V}_{k} \Leftrightarrow \mathbf{V}_{k}^{*} \mathbf{r} = \mathbf{0}

すなわち、残差の$\mathcal{V}_{k}$への射影が$\mathbf{0}$になるような解を求める。これより、$\mathbf{V}_{k}^{*}\mathbf{V}_{k} = \mathbf{I}_{k}$なので、

\begin{align}
& \mathbf{V}_{k}^{*} \left( \mathbf{A} \mathbf{u} - \theta \mathbf{u} \right) = \mathbf{V}_{k}^{*} \left( \mathbf{A} \mathbf{V}_{k} \mathbf{y} - \theta \mathbf{V}_{k} \mathbf{y} \right) = \mathbf{0} \\
& \therefore \mathbf{V}_{k}^{*} \mathbf{A} \mathbf{V}_{k} \mathbf{y} = \theta \mathbf{y} 
\end{align}

を得る。$\mathbf{V}_{k}^{*} \mathbf{A} \mathbf{V}_{k} \in \mathbb{C}^{k \times k}$は$\mathbf{A}$の部分空間への射影であり、$\{\theta, \mathbf{y}\}$に関する固有値問題は元の問題に比べて規模が非常に小さくなっているため、QR法などによって容易に計算することができる。得られた固有値$\theta$と固有ベクトル$\mathbf{u} = \mathbf{V}_{k} \mathbf{y}$は元の問題に対する残差を最小にする近似解であり、それぞれRitz値およびRitzベクトルと呼ばれる。

部分空間としてKrylov部分空間を用いることを考える。$\mathcal{K}_{k} \left( \mathbf{A}, \mathbf{x}_{1} \right)$から生成されるKrylov列からなる行列を$\mathbf{X}_{1}^{k}$とし、そのQR分解を$\mathbf{X}_{1}^{k} = \mathbf{V}_{1}^{k} \mathbf{R}_{k}$とすると、

\begin{align}
\mathbf{A} \mathbf{X}_{1}^{k} &= \mathbf{A} \mathbf{V}_{1}^{k} \mathbf{R}_{k} \\
\Rightarrow \mathbf{A} \mathbf{V}_{1}^{k} &= \mathbf{A} \mathbf{X}_{1}^{k} \mathbf{R}_{k}^{-1} = \mathbf{X}_{2}^{k+1} \mathbf{R}_{k}^{-1} \\
&= \mathbf{X}_{1}^{k+1} \begin{bmatrix} \mathbf{0}^{\top} \\ \mathbf{I}_{k} \end{bmatrix} \mathbf{R}_{k}^{-1} \\
&= \mathbf{V}_{1}^{k+1} \mathbf{R}_{k+1} \begin{bmatrix} \mathbf{0}^{\top} \\ \mathbf{I}_{k} \end{bmatrix} \mathbf{R}_{k}^{-1}
\end{align}

ここで、

\mathbf{H}_{k+1,k} = \mathbf{R}_{k+1} \begin{bmatrix} \mathbf{0}^{\top} \\ \mathbf{I}_{k} \end{bmatrix} \mathbf{R}_{k}^{-1}

とおくと、$\mathbf{H}_{k+1,k}$は$(k+1) \times k$のHessenberg行列になっており、

\mathbf{H}_{k+1,k} = \begin{bmatrix} \mathbf{H}_{k} \\ h_{k+1,k} \mathbf{e}_{k}^{\top} \end{bmatrix}

と表される。ただし、$h_{k+1,k}$は$\mathbf{H}_{k+1,k}$の$(k+1,k)$成分である。これより、次式のArnoldi分解が成り立つ。

\begin{align}
\mathbf{A} \mathbf{V}_{1}^{k} &= \mathbf{V}_{1}^{k+1} \mathbf{H}_{k+1,k} \\
\mathbf{A} \mathbf{V}_{1}^{k} &= \mathbf{V}_{1}^{k} \mathbf{H}_{k} + h_{k+1,k} \mathbf{v}_{k+1} \mathbf{e}_{k}^{\top}
\end{align}

これより、固有値問題の近似解は

\mathbf{H}_{k} \mathbf{y} = \theta \mathbf{y}

という$\mathbf{H}_{k}$に関する小規模な固有値問題を解くことで計算できる。このとき、残差ベクトルは

\mathbf{r} = h_{k+1,k} y_ {k} \mathbf{v}_{k+1}

となるので³、$\| \mathbf{r} \|_{2}$が十分小さくなるまで部分空間を拡大しながら反復計算を行う。

Gram-Schmidt直交化

Arnoldi法の計算はKrylov列を陽に構成してからそれをQR分解し、さらにHessenberg行列を求める...といったやり方で行う必要はない。いま、$\mathbf{v}_{1}, \mathbf{v}_{2}, \cdots, \mathbf{v}_{j}$（$j \leq k$）までの基底が求まっているとして、新しく$\mathbf{v}_{j+1}$を計算することを考える。Arnoldi分解の式$\mathbf{A} \mathbf{V}_{1}^{k} = \mathbf{V}_{1}^{k+1} \mathbf{H}_{k+1,k}$の$j$列目は以下のように表される。

\mathbf{A}\mathbf{v}_{j} = h_{1,j} \mathbf{v}_{1} + h_{2,j} \mathbf{v}_{2} + \cdots + h_{j,j} \mathbf{v}_{j} +h_{j+1,j} \mathbf{v}_{j+1}

これより、Gram-Schmidtの正規直交化法を用いて$\mathbf{A}\mathbf{v}_{j}$を$\mathbf{v}_{1}, \mathbf{v}_{2}, \cdots, \mathbf{v}_{j}$に対して直交化すれば$\mathbf{v}_{j+1}$が求められる。また、Hessenberg行列の要素は直交基底の計算過程で自然に求まる。

まず、

\mathbf{w}_{1} = \mathbf{A} \mathbf{v}_{j}

とおく。上式と$\mathbf{v}_{1}$の内積を取ると、$\mathbf{A} \mathbf{v}_{j} = \sum_{i=1}^{j+1} h_{i,j} \mathbf{v}_{i}$と正規直交性より、対応する展開係数が抽出される。

\mathbf{v}_{1}^{*} \mathbf{w}_{1} = \left( \mathbf{w}_{1} , \mathbf{v}_{1} \right)  = h_{1,j}

ここで、$\left( \cdot , \cdot \right)$はベクトルの内積である。これより、

\mathbf{w}_{2} = \mathbf{w}_{1} - h_{1,j} \mathbf{v}_{1}

とすると、上述の関係から、

\mathbf{w}_{2} = h_{2,j} \mathbf{v}_{2} + \cdots + h_{j,j} \mathbf{v}_{j} +h_{j+1,j} \mathbf{v}_{j+1}

となるので、$\left( \mathbf{w}_{2} , \mathbf{v}_{1} \right) = 0 \Leftrightarrow \mathbf{w}_{2} \perp \mathbf{v}_{1}$が成り立つ。同様に、

\begin{align}
h_{i,j} &= \left( \mathbf{w}_{i} , \mathbf{v}_{i} \right)  \\
\mathbf{w}_{i+1} &= \mathbf{w}_{i} - h_{i,j} \mathbf{v}_{i}
\end{align}

を$i = 1,\cdots,j$まで計算すると、最終的には$j+1$番目の項だけが残り、

\mathbf{w}_{j+1} = h_{j+1,j} \mathbf{v}_{j+1}

が成り立つ。したがって、新しい直交基底は次式で求められる。

\begin{align}
h_{j+1,j} &= \left\| \mathbf{w}_{j+1} \right\|_{2} \\
\mathbf{v}_{i+1} &= \mathbf{w}_{j+1} \left/ \left\| \mathbf{w}_{j+1} \right\|_{2} \right.
\end{align}

実際には、$\mathbf{w}_{i}$は逐次代入で上書きしながら計算すればよく、個別に保存する必要はない。

なお、上記で説明したアルゴリズムでは修正Gram-Schmidt（MGS）法を使用している。これは古典Gram-Schmidt（CGS）法に比べて丸め誤差の伝播を抑えられるためよく用いられるが、そのことが原因で$\mathbf{w}_{i+1}$の更新則は並列効率が悪い。そこで、元論文では並列効率に優れる再直交化処理付きのCGS法を採用している。これは単純にCGS法による直交化を2回繰り返す手法である（ここを参照）。

FOA based DMD

というわけで、前置きが長くなってしまったが、いよいよ本題のFOA based DMDについて述べる。

アルゴリズム

Schmid(2010)でもすでに指摘されているように、Krylov DMDはArnoldi分解の一種と考えることができる。しかし、通常の固有値問題では行列$\mathbf{A}$が既知である一方、DMDでは$\mathbf{A}$が未知であるため、直接Hessenberg行列を求めることができない。⁴

しかし、FOA based DMDはこの「$\mathbf{A}$が既知でなければならない」という問題をアルゴリズムの工夫によって回避し、Arnoldi法に基づく方法によってDMDの計算を可能にしている。以下、この点について説明する。

このアルゴリズムはKrylov DMDの部分で述べた「データ行列の列空間とスナップショットによって生成されるKrylov部分空間が等しい」という仮定に基づいている。データ行列の列空間とスナップショットの正規直交基底の張る空間も等しいので、

\begin{align}
\mathcal{R}\left( \mathbf{X}_{1}^{m-1} \right) &= \mathcal{R}\left( \mathbf{V}_{1}^{m-1} \right) = \mathcal{K}_{m-1} \left( \mathbf{A}, \mathbf{x}_{1} \right) \\
\Leftrightarrow \operatorname*{span} \left\{ \mathbf{x}_{1}, \mathbf{x}_{2}, \cdots,  \mathbf{x}_{m-1} \right\}
&= \operatorname*{span} \left\{ \mathbf{v}_{1}, \mathbf{v}_{2}, \cdots, \mathbf{v}_{m-1} \right\} \\
&= \operatorname*{span} \left\{ \mathbf{x}_{1}, \mathbf{A}\mathbf{x}_{1}, \cdots, \mathbf{A}^{m-2} \mathbf{x}_{1} \right\}
\end{align}

となる。これより、$j$番目のスナップショット$\mathbf{x}_{j}$は正規直交基底$\mathbf{v}_{1}, \mathbf{v}_{2}, \cdots, \mathbf{v}_{j}$の線形結合で表される。展開係数を$b_{i,j}$とすれば、⁵

\begin{align}
\mathbf{x}_{j} &= b_{1,j} \mathbf{v}_{1} + b_{2,j} \mathbf{v}_{2} + \cdots + b_{j,j} \mathbf{v}_{j} = \sum_{i=1}^{j} b_{i,j} \mathbf{v}_{i} \\
&= \mathbf{V}_{1}^{j} \mathbf{b}_{1:j,j}
\end{align}

であり、行列でまとめて書けば、

\mathbf{X}_{1}^{j} = \mathbf{V}_{1}^{j} \mathbf{B}_{j}

となる。ここで、$\mathbf{B}_{j}$は展開係数をまとめた行列であり、$\mathbf{b}_{1:j,j}$は$\mathbf{B}_{j}$の$j$列目の$j$行目までの成分を取り出したベクトルである（以下、下付き添字にコロン$:$が含まれるものはこうしたベクトルである）。上式はQR分解に等しく、$\mathbf{B}_{j}$は上三角行列になる。これより、$\mathbf{x}_{j+1}$に対してArnoldi分解に似た次式のような関係式が成り立つ。

\begin{align}
\mathbf{x}_{j+1} &= \mathbf{A} \mathbf{x}_{j} = \mathbf{A} \mathbf{V}_{1}^{j} \mathbf{b}_{1:j,j} \\
&= \mathbf{A} \mathbf{V}_{1}^{j-1} \mathbf{b}_{1:j-1,j} + \mathbf{A} \mathbf{v}_{j} b_{j,j} \\
&= \mathbf{V}_{1}^{j} \mathbf{H}_{j,j-1} \mathbf{b}_{1:j-1,j} + \mathbf{A} \mathbf{v}_{j} b_{j,j} \\
\end{align}

ここで、$\mathbf{H}_{j,j-1}$は$j \times (j-1)$のHessenberg行列である。上式を$\mathbf{A} \mathbf{v}_{j}$について書き直せば、修正された正規直交基底の生成プロセスが得られる。

\begin{align}
\mathbf{A} \mathbf{v}_{j} 
&= \frac{1}{b_{j,j}} \left( \mathbf{x}_{j+1} - \mathbf{V}_{1}^{j} \mathbf{H}_{j,j-1} \mathbf{b}_{1:j-1,j} \right) \\
&= \frac{1}{b_{j,j}} \left( \mathbf{x}_{j+1} - \mathbf{v}_{1} \sum_{i=1}^{j-1} h_{1,i} b_{i,j} - \sum_{k=2}^{j-1} \mathbf{v}_{k} \sum_{i=k-1}^{j-1} h_{k,i} b_{i,j} \right)
\end{align}

Arnoldi分解の$j$ステップ目においてHessenberg行列の$j$列目が求まっているならば、$\mathbf{A} \mathbf{v}_{j} = \sum_{i=1}^{j+1} h_{i,j} \mathbf{v}_{i}$の関係式より、

\begin{align}
\mathbf{x}_{j+1} = \mathbf{V}_{1}^{j} \mathbf{H}_{j,j-1} \mathbf{b}_{1:j-1,j} + \mathbf{V}_{1}^{j+1} \mathbf{h}_{1:j+1,j} b_{j,j} = \mathbf{V}_{1}^{j+1} \mathbf{H}_{j+1,j} \mathbf{b}_{1:j,j} \\
\end{align}

上式の最右辺はQR分解であるから、改めて$\mathbf{x}_{j+1} = \mathbf{V}_{1}^{j+1} \mathbf{b}_{1:j+1,j+1}$と等しいと置いて両辺の係数を比較すると、

\begin{align}
\mathbf{b}_{1:j+1,j+1} = \mathbf{H}_{j+1,j} \mathbf{b}_{1:j,j}
\end{align}

より、以下のようなHessenberg行列と上三角行列の関係式が成り立つ。

\begin{align}
b_{1,j+1} &= \sum_{i=1}^{j} h_{1,i} b_{i,j} \\
b_{k,j+1} &= \sum_{i=k-1}^{j} h_{k,i} b_{i,j} \quad (k=2,3,\cdots,j+1) \\
\end{align}

したがって、Arnoldi分解の過程で求まるHessenberg行列の要素からQR分解の上三角行列の要素が求められる。

計算の初期過程

アルゴリズムの全体については写すのが面倒なので元論文の方を見てほしいが⁶、そこでの処理で気になった部分があったので書いておく。

DMDでは$\mathbf{A}$が未知なので、Arnoldi分解から計算を始めることはできない。そこで、計算はQR分解の初期値からスタートする。すなわち、

\begin{align}
\mathbf{v}_{1} &= \frac{ \mathbf{x}_{1} }{ \| \mathbf{x}_{1} \|_{2} } \\
b_{1,1} &= \| \mathbf{x}_{1} \|_{2}
\end{align}

として順次計算を進めていくのだが、一番最初の更新式

b_{1,j+1} = \left( \mathbf{h}_{1,1:j-1} , \mathbf{b}_{1:j-1,j} \right)

は、$j=1$のとき明らかにインデックスの範囲がおかしい。この点について元論文では特に何も言っていないので理解に苦しむところなのだが、ひとまずこの処理は$j=1$のとき、

b_{1,2} = h_{1,1} b_{1,1}

になると考えてみる（通常は$j=2$のときにこうなるはず）。そうすると、$h_{1,1}$はまだ未知の量である。すべての行列は$0$で初期化されているため、その値を参照するなら$h_{1,1} = 0$であり、これより、$b_{1,2} = b_{2,2} = 0$となる。

このまま$j=1$のループを計算していくと、

\begin{align}
(1)\ & \mathbf{w}_{1} = \frac{\mathbf{x}_{2}}{b_{1,1}} \\
(2)\ & h_{1,1} = \left( \mathbf{w}_{1} , \mathbf{v}_{1} \right) \\
(3)\ & \mathbf{w}_{2} = \mathbf{w}_{1} - h_{1,1} \mathbf{v}_{1} \\
(4)\ & h_{2,1} = \| \mathbf{w}_{2} \|_{2} ,\ \mathbf{v}_{2} = \frac{ \mathbf{w}_{2} }{ \| \mathbf{w}_{2} \|_{2} } \\
(5)\ & b_{1,2} = h_{1,1} b_{1,1} \\
(6)\ & b_{2,2} = h_{2,1} b_{1,1} \\
\end{align}

ただし、本来代入で上書きする$\mathbf{w}$に番号をつけて区別しており、加えて再直交化処理のための計算は省いた。

ところで、QR分解を$\mathbf{X} = \mathbf{VB}$とすると、Gram-Schmidt直交化によって求められる上三角行列の要素は次のように表される。

\begin{align}
\underset{\mathbf{X}}{\underline{
\begin{bmatrix}
\mathbf{x}_{1} & \mathbf{x}_{2} & \cdots & \mathbf{x}_{k}
\end{bmatrix}
}}
=
\underset{\mathbf{V}}{\underline{
\begin{bmatrix}
\mathbf{v}_{1} & \mathbf{v}_{2} & \cdots & \mathbf{v}_{k}
\end{bmatrix}
}}
\underset{\mathbf{B}}{\underline{
\begin{bmatrix}
\left( \mathbf{x}_{1}, \mathbf{v}_{1} \right) & \left( \mathbf{x}_{2}, \mathbf{v}_{1} \right) & \cdots & \left( \mathbf{x}_{k}, \mathbf{v}_{1} \right) \\
 & \left( \mathbf{x}_{2}, \mathbf{v}_{2} \right) & \cdots & \left( \mathbf{x}_{k}, \mathbf{v}_{2} \right) \\
& & \ddots & \vdots \\
& & & \left( \mathbf{x}_{k}, \mathbf{v}_{k} \right)
\end{bmatrix}
}}
\end{align}

これより、

\begin{align}
b_{1,2} &= h_{1,1} b_{1,1} = \left( \mathbf{w}_{1}, \mathbf{v}_{1} \right) b_{1,1} = \left( \frac{\mathbf{x}_{2}}{b_{1,1}}, \mathbf{v}_{1} \right) b_{1,1} = \left( \mathbf{x}_{2}, \mathbf{v}_{1} \right) \\
b_{2,2} &= h_{2,1} b_{1,1} = \left( \mathbf{w}_{2}, \mathbf{v}_{2} \right) b_{1,1} = \left( \mathbf{w}_{1} - \mathbf{v}_{1} h_{1,1}, \mathbf{v}_{2} \right) b_{1,1} \\
&= \left( \mathbf{w}_{1}, \mathbf{v}_{2} \right) b_{1,1} - \left( \mathbf{v}_{1} h_{1,1}, \mathbf{v}_{2} \right) b_{1,1} = \left( \mathbf{x}_{2}, \mathbf{v}_{2} \right) \\
\end{align}

となるので、インデックスの範囲が$0$となる場合はその部分を$0$とおくことで正しく計算できることが確認できた（実際の実装では処理を分けるなど何らかの工夫が必要と思われる）。⁷

上述の手順によってデータ行列$\mathbf{X}_{1}^{m}$から直交基底$\mathbf{V}_{1}^{m}$、Hessenberg行列$\mathbf{H}_{m,m-1}$、上三角行列$\mathbf{B}_{m}$が求まるので、後は固有値問題$\mathbf{H}_{m-1} \mathbf{y} = \theta \mathbf{y}$を解き、動的モードを$\mathbf{u} = \mathbf{V}_{1}^{m-1} \mathbf{y}$から計算すればよい。このとき、データ行列のランク落ちを考慮して$\mathbf{B}_{m-1}$に対して打ち切り特異値分解（Truncated SVD）を適用することもできる。また、上記で紹介したのはいわゆるバッチ処理のアルゴリズムだが、元論文には逐次処理型（Streaming form）のアルゴリズムも載っているので、詳細については元論文をご参照いただきたい。元々Gram-Schmidt直交化に基づくアルゴリズムなので、逐次処理にすぐ変更できることは容易に想像できると思う。

まとめ

Anantharamu & Mahesh(2019)のFOA based DMDのアルゴリズムについてまとめた。ぶっちゃけ元論文の計算量オーダーとか後方安定性の議論なんかは全然フォローできていないが、一番気になっていた「結局何をやってるのか」については大枠を把握できたと思うので、今回はこれでよしとする。

次回、動作確認を兼ねた実装編に続く。

1. 余談という名の備忘録。Qiitaのインライン数式で下付き文字を書く場合、アンダースコアをエスケープ\_しなければならない（Markdownの強調文字として認識されてしまうため）。また、ノルム記号$\| \cdot \|$を書く場合、\\| \cdot \\|のように二重にエスケープする必要がある（波括弧$\{\}$も同じ）。逆に数式コードブロックの中でエスケープ文字を使うとエラーになる。コピペとかがやりにくいのでこの仕様を改善してくれるとありがたいんだが、TeX記法だけにおもねるわけにもいかないんだろうなぁ。 ↩

2. スナップショットの数が十分大きければ部分空間も十分拡大されるので、この性質がおおむね成り立つと思ってよい。特に（準）周期的な現象が対象であれば、新しい時刻のスナップショットは以前のスナップショットの組み合わせで表現できる。これはモード分解が成り立つために必要な条件ともいえる。そのため、データの数が不足している場合や周期性がない現象から得られたデータに対してはDMDはうまく機能しない。 ↩

3. $\mathbf{r} = \mathbf{A} \mathbf{V}_{1}^{k} \mathbf{y} - \theta \mathbf{V}_{1}^{k} \mathbf{y}$にArnoldi分解の式$\mathbf{A} \mathbf{V}_{1}^{k} = \mathbf{V}_{1}^{k} \mathbf{H}_{k} + h_{k+1,k} \mathbf{v}_{k+1} \mathbf{e}_{k}^{\top}$を代入するだけ。 ↩

4. $\mathbf{X}_{1}^{m-1} = \mathbf{QR}$とすれば、$\mathbf{AQ} \approx \mathbf{QH}$より、同伴行列$\mathbf{S}$を用いて$\mathbf{H} = \mathbf{RSR}^{-1}$となる。もちろん、これは$\mathbf{S}$が既知でなければならず、$\mathbf{H}$を求めることは余分な計算になってしまう。 ↩

5. 元論文では$\beta$だが、便宜上ちょくちょく記号を変更している点に注意。 ↩

6. Qiita（というかMathJax）ではalgorithmパッケージなどが使えないので、数式を使った疑似コードが書けないのが少々残念。コードブロックで実際のコード書けばいいじゃんという考え方もあるが、その場合は依拠する言語の文法に気を使ったり、実装したときに動くようにしなきゃ的なことに気を使ってしまうので、疑似コードで書きたい場合もあると思う。 ↩

7. ここの部分が理解できずに随分時間を取られた...実際に手を動かして計算してみることの重要性を再確認した。 ↩