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

ひまつぶしにどうぞ。

特異値分解を用いた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でやればいいよ、という話であった。当初考えていたより長くなってしまったが、復習は大事なのでよしとしよう。