Skip to main content

3 posts tagged with "svd-for-beginners"

View All Tags

特異値分解入門(1)

· 5 min read

互いに相関を持った正規乱数のセット(多変量正規乱数)を生成すること、これが金融工学のモンテカルロシミュレーションにおいて最も重要な技術となります。 NN 系列の正規乱数列を生成する場合を考えます。 これらの乱数列の相関構造を記述するのが相関行列と呼ばれる N×NN\times N の正方行列 C{\boldsymbol C} です。 最終的に

C=QQt{\boldsymbol C}={\boldsymbol Q}{\boldsymbol Q}^\text{t}

という形に分解することが目標となります(ここで t は転置行列です)。分解して得られた行列 Q{\boldsymbol Q} を使って実際に相関乱数を生成する方法については、これも一般的な数値計算の教科書に載っていますのでそちらをご覧ください。

すぐ後で説明するように、この分解にはコレスキー分解を使えと、多くの教科書で教えていますが、実務上コレスキー分解ではうまくいかない場合があります。
NtRand あるいは Numerical Technologies 社の製品では相関行列の分解に特異値分解という方法を採用しています。
これから3回にわたって、なぜコレスキー分解ではダメなのか、そしてなぜ特異値分解だとうまくいくのかを説明します。

相関行列について

相関行列 C{\boldsymbol C} は、その ij 成分が確率変数 XiX_{i} と XjX_{j} との相関係数 ρij\rho_{ij} である行列です。 この定義から対角成分 ρii\rho_{ii} は確率変数 XiX_{i} と自分自身との相関、つまり 1 となります。また、ρij=ρji\rho_{ij}=\rho_{ji} も定義から明らかです。
つまりこの行列は、

  • N×NN\times N正方行列
  • 成分は実数
  • 対角成分は 1
  • 対角成分以外は 1ρij1-1\leq\rho_{ij}\leq 1
  • 対称(ρij=ρji\rho_{ij}=\rho_{ji}

という性質を持っています。
更に相関係数の定義から、この行列は半正定値であることが示されます(半正定値って何だ?それは行列の固有値が全て0以上であること。固有値って何かって?それは勉強してください)。 つまり相関行列はその定義から半正定値実対称行列であるということです。

コレスキー分解 (Cholesky decomposition)

C=QQt{\boldsymbol C}={\boldsymbol Q}{\boldsymbol Q}^\text{t}

という分解を行う技として、どの教科書にも書いてあるのが**コレスキー分解という方法です。 どの教科書にも載っているので、ここでは説明しません。
この手法は正定値実対称行列を目的の形に分解してくれるので、正定値実対称行列である相関行列に対して使えば目標達成!メデタシ、メデタシとなる筋書きです。
しかし...
実はコレスキー分解は容易に失敗します。その代表的な状況が
観測データが少ない場合です。 確率変数の数が100個で、その観測データの数が100個未満の場合は、定義に従った正しい相関行列であってもコレスキー分解は失敗します。 なぜなら、この状況で作られる相関行列は必然的に正定値行列になるからです。コレスキー分解は正定値行列専用**です。
それでは、少なくとも正しく作られた相関行列では間違いなく

C=QQt{\boldsymbol C}={\boldsymbol Q}{\boldsymbol Q}^\text{t}

と分解してくれるより強力な方法はないでしょうか?

特異値分解入門(2)

· 7 min read

ここで相関行列 C{\boldsymbol C} が満たしているはずの性質を今一度列挙してみましょう。

  • 成分 ρij\rho_{ij} は実数
  • 対称行列 ρij=ρji\rho_{ij}=\rho_{ji}
  • 半正定値
  • 対角成分 ρii=1\rho_{ii}=1
  • 対角成分以外 1ρij1-1\leq\rho_{ij}\leq 1

この相関行列 C{\boldsymbol C} を、

C=QQt{\boldsymbol C}={\boldsymbol Q}{\boldsymbol Q}^\text{t}

という形に分解するのが目的です(覚えてましたか?)
前回その方法としてコレスキー分解を紹介しましたが、どうも安心して使えなさそうだとわかりました。
次に紹介するのは固有値分解です。これは半正定値実対称行列を分解してくれます。
ではここからは線形代数の勉強です(証明などの詳細は省きます。教科書に必ず載っていますので調べてみてください)。

事実1

任意の実正方行列 A{\boldsymbol A} は、ある正方行列 U{\boldsymbol U} とこれまた正方行列、ただし対角成分のみしか値がない対角行列 W{\boldsymbol W} を使って、

A=UWU1{\boldsymbol A}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^{-1}

という形に分解される(これを行列のスペクトル分解といいます)。
ここで行列 U{\boldsymbol U} は行列 A{\boldsymbol A} の固有ベクトルを横に並べたもので、対角行列 W{\boldsymbol W} の成分は行列 A{\boldsymbol A} の固有値です。

事実2

更に行列 A{\boldsymbol A} が対称行列であった場合、行列の固有ベクトルが全て直交する。すると、

U1=Ut{\boldsymbol U}^{-1}={\boldsymbol U}^\text{t}

となる。というわけで

A=UWU1=UWUt{\boldsymbol A}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^{-1}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^\text{t}

事実3

更に更に行列 A{\boldsymbol A} が半正定値だとすれば、行列の固有値は全て0以上。つまり行列 W{\boldsymbol W} は、

(λ1000λ20vdots00λN)=(λ1000λ20vdots00λN)(λ1000λ20vdots00λN)=WW{\small \begin{pmatrix}\lambda_1&0&\cdots&0\\0&\lambda_2&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\lambda_{N}\end{pmatrix}=\begin{pmatrix}\sqrt{\lambda_1}&0&\cdots&0\\0&\sqrt{\lambda_2}&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{\lambda_{N}}\end{pmatrix}\begin{pmatrix}\sqrt{\lambda_1}&0&\cdots&0\\0&\sqrt{\lambda_2}&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{\lambda_{N}}\end{pmatrix}}={\boldsymbol W}^\prime {\boldsymbol W}^\prime

と出来るではないか!すると、

A=UWU1=UWUt=UWWUt{\boldsymbol A}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^{-1}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^\text{t}={\boldsymbol U}{\boldsymbol W}^\prime {\boldsymbol W}^\prime {\boldsymbol U}^\text{t}

となります。ここで (AB)t=BtAt({\boldsymbol A}{\boldsymbol B})^\text{t}={\boldsymbol B}^\text{t}{\boldsymbol A}^\text{t} という事実(積の順番が入れ替わる)と、対角行列は転置しても変化しないこと(この場合 (W)t=W({\boldsymbol W}^\prime)^\text{t}={\boldsymbol W}^\prime)に注意すると、

A=UWU1=UWUt=UWWUt=UW(UW)t{\boldsymbol A}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^{-1}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^\text{t}={\boldsymbol U}{\boldsymbol W}^\prime {\boldsymbol W}^\prime {\boldsymbol U}^\text{t}={\boldsymbol U}{\boldsymbol W}^\prime ({\boldsymbol U}{\boldsymbol W}^\prime)^\text{t}

やった!! UW=Q{\boldsymbol U}{\boldsymbol W}^\prime = {\boldsymbol Q} とおけば目標達成じゃないか!

A=QQt{\boldsymbol A}={\boldsymbol Q}{\boldsymbol Q}^\text{t}

しかし...
実務上(ここがポイント)、相関行列がいつもいつもタチの良い形であるとは限らないのです!もちろん実対称行列であることはどんな場合も間違いないですが(複素数の相関や、100×234100\times 234の相関行列を作れるものなら作ってみろ)、 半正定値を満たさない場合は往々にしてありえるのです。例えば...

  • データに欠損があって、欠損データを適当に埋める
  • 乱数列のうち、いくつかのペアについては相関が(観測とは無関係に)与えられている
  • 測定誤差

などの場合、相関係数に矛盾が出ることがあるのです。相関係数が矛盾するとはどういうことか、ひとつ最も大げさな例を示してみましょう。
以下の3変数の相関行列を見てください。

C=(111111111){\boldsymbol C}=\begin{pmatrix}1&1&-1\\1&1&1\\-1&1&1\end{pmatrix}

実数であるとか、対称であるとか相関行列の性質は確かに満たしているんですが...何がおかしいか分かりますか?

  • 第1変数と第2変数の相関係数 ρ12\rho_{12} は 1。第1変数と第2変数は完全相関で全く同じ値になる
  • 第2変数と第3変数の相関係数 ρ23\rho_{23} は 1。第2変数と第3変数は完全相関で全く同じ値になる

ってことは必然的に第1変数と第3変数は同じ値になるはずです。しかし与えられた相関係数 ρ13\rho_{13} は -1 になっている!つまり矛盾しているわけです。
要するに相関係数がありえない組み合わせになっている場合がありえるのです(さすがに例のような大げさなのはないでしょうけど)。このとき、相関行列は半正定値ではなくなってコレスキー分解が機能しなくなります。

実務上、相関行列が与えられたら(それに矛盾があろうとなかろうと)計算がストップしてしまうことは望ましくありません。とにかく安定して何らかの答えを出すことが求められます。 もちろん正しい相関行列(半正定値実対称行列)は正しく分解されることが前提です。

そこで NtRand では相関行列の分解に特異値分解(Singular Value Decomposition,SVD)という手法を採用しています。

特異値分解入門(3)

· 7 min read

不正な相関行列 C{\boldsymbol C} が与えられると何が起きるが考えてみましょう。
前回の解説の "事実3" で悲劇が起きます。つまり、

(λ1000λ20vdots00λN)=(λ1000λ20vdots00λN)(λ1000λ20vdots00λN){\small \begin{pmatrix}\lambda_1&0&\cdots&0\\0&\lambda_2&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\lambda_{N}\end{pmatrix}=\begin{pmatrix}\sqrt{\lambda_1}&0&\cdots&0\\0&\sqrt{\lambda_2}&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{\lambda_{N}}\end{pmatrix}\begin{pmatrix}\sqrt{\lambda_1}&0&\cdots&0\\0&\sqrt{\lambda_2}&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{\lambda_{N}}\end{pmatrix}}

とする時になって、ハッと気が付くのです...「固有値が負だ!」と。
そしてコンピューターは無情にも「負の平方根は計算できません」というメッセージを残して息絶えます。
ちなみに第1回に載せた不正な相関行列、

C=(111111111){\boldsymbol C}=\begin{pmatrix}1&1&-1\\1&1&1\\-1&1&1\end{pmatrix}

の固有値は 1,2,2-1,2,2 となります。

この困難を打破してとにかく答えを出したい!そこでいよいよ特異値分解の登場です。

先ずは特異値分解で任意の行列(なんと正方行列でなくてもいい!)A{\boldsymbol A} はどのように分解されるか見てみましょう。

A=UWVt{\boldsymbol A}={\boldsymbol U}{\boldsymbol W}{\boldsymbol V}^\text{t}

となります。とりあえずここでは A{\boldsymbol A} として実対称正方行列であるところの相関行列 C{\boldsymbol C} だけを対象にして順に説明していきましょう。

行列 C{\boldsymbol C} から M=CCt{\boldsymbol M}={\boldsymbol C}{\boldsymbol C}^\text{t} なる行列を計算します。この行列は自動的に、

  • 正方行列
  • 対称行列
  • 半正定値行列

となるんです!(もちろんもとの行列が実数の行列ならば M{\boldsymbol M} も実数)。最後の特徴に注目してください。簡単に言うと、元の行列(C{\boldsymbol C})2乗することで負だった固有値があったとしてもそれが正になるということです。
というわけで行列 M{\boldsymbol M} は互いに直交する固有ベクトルを持ち、固有値は全て0以上となります。この固有ベクトルを並べた行列が U{\boldsymbol U} です。
次に、N=CtC{\boldsymbol N}={\boldsymbol C}^\text{t}{\boldsymbol C} なる行列を作ります。この行列は M{\boldsymbol M} の場合と全く同様に

  • 正方行列
  • 対称行列
  • 半正定値行列

となります。N{\boldsymbol N} の固有ベクトルを並べた行列が V{\boldsymbol V} です。
最後に W{\boldsymbol W} は行列 C{\boldsymbol C} の固有値の絶対値を対角に並べた対角行列になります。

ここまでの手順で分かるように、各行列 U{\boldsymbol U}V{\boldsymbol V}W{\boldsymbol W} は元の行列 C{\boldsymbol C} がどんなものであれ安定的に得ることができます(計算は止まらない!)。 前にも述べましたがこれが最も重要な点です。

とにもかくにも行列 C{\boldsymbol C} は、

C=UWVt{\boldsymbol C}={\boldsymbol U}{\boldsymbol W}{\boldsymbol V}^\text{t}

という形に分解されました。次に

W=(λ1000λ20vdots00λN)=(λ1000λ20vdots00λN)(λ1000λ20vdots00λN)=WW{\small {\boldsymbol W}=\begin{pmatrix}|\lambda_1|&0&\cdots&0\\0&|\lambda_2|&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&|\lambda_{N}|\end{pmatrix}=\begin{pmatrix}\sqrt{|\lambda_1|}&0&\cdots&0\\0&\sqrt{|\lambda_2|}&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{|\lambda_{N}|}\end{pmatrix}\begin{pmatrix}\sqrt{|\lambda_1|}&0&\cdots&0\\0&\sqrt{|\lambda_2|}&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{|\lambda_{N}|}\end{pmatrix}}={\boldsymbol W}^\prime {\boldsymbol W}^\prime

とします(対角成分は固有値の絶対値なので安心して平方根がとれる)。すると、

C=UWVt=UWWVt=UW(VW)t{\boldsymbol C}={\boldsymbol U}{\boldsymbol W}{\boldsymbol V}^\text{t}={\boldsymbol U}{\boldsymbol W}^\prime {\boldsymbol W}^\prime {\boldsymbol V}^\text{t}={\boldsymbol U}{\boldsymbol W}^\prime({\boldsymbol V}{\boldsymbol W}^\prime)^\text{t}

となります。

あと一息です

さてここで行列 C{\boldsymbol C} が正しい相関行列、つまり半正定値だとどうなるか。その場合は U=V{\boldsymbol U}={\boldsymbol V} となって特異値分解は前回の解説の結果と同じになります。 つまり Q=VW{\boldsymbol Q}={\boldsymbol V}{\boldsymbol W}^\prime とすることで、相関行列は C=QQt{\boldsymbol C}={\boldsymbol Q}{\boldsymbol Q}^\text{t} となって目標達成です。
しかしそうでない場合はどうでしょう。その場合もやはり同じく Q=VW{\boldsymbol Q}={\boldsymbol V}{\boldsymbol W}^\prime とすることで手をうちませんか?もちろんこの場合は CQQt{\boldsymbol C}\neq {\boldsymbol Q}{\boldsymbol Q}^\text{t} ではありますが(UV{\boldsymbol U}\neq{\boldsymbol V}だから)...でも、

  • 安定的に答えが出る
  • 正しい相関行列の場合は正しい答えが出る

という望ましい性質を持ってます。
因みに正しい相関行列ではない場合にも、QQt{\boldsymbol Q}{\boldsymbol Q}^\text{t} の各成分の符号は元の行列 C{\boldsymbol C} のそれに等しくなります(正の相関が負になることはない)。

再び

C=(111111111){\boldsymbol C}=\begin{pmatrix}1&1&-1\\1&1&1\\-1&1&1\end{pmatrix}

を例に計算してみましょう。特異値分解の結果は、

C=(0.5770.7390.3480.5770.6710.4660.5770.0680.814)U(100020002)W(0.5770.5770.5770.7390.6710.0680.3480.4660.814)Vt{\boldsymbol C}=\underbrace{\begin{pmatrix}0.577&-0.739&-0.348\\-0.577&-0.671&0.466\\0.577&0.068&0.814\end{pmatrix}}_{{\boldsymbol U}}\underbrace{\begin{pmatrix}1&0&0\\0&2&0\\0&0&2\end{pmatrix}}_{{\boldsymbol W}^\prime}\underbrace{\begin{pmatrix}-0.577&0.577&-0.577\\-0.739&-0.671&0.068\\-0.348&0.466&0.814\end{pmatrix}}_{{\boldsymbol V}^\text{t}}

となります(UV{\boldsymbol U}\neq {\boldsymbol V}であることが分かりますか?第1固有ベクトルの符号が逆転しています!こうすることで、固有値の方を正に揃えているという訳です)。そして、

Q=VW=(0.5770.5770.5770.7390.6710.0680.3480.4660.814)(100020002){\boldsymbol Q}={\boldsymbol V}{\boldsymbol W}^\prime=\begin{pmatrix}-0.577&0.577&-0.577\\-0.739&-0.671&0.068\\-0.348&0.466&0.814\end{pmatrix}\begin{pmatrix}1&0&0\\0&\sqrt{2}&0\\0&0&\sqrt{2}\end{pmatrix}

となります。ここから QQt{\boldsymbol Q}{\boldsymbol Q}^\text{t} を計算すると、

QQt=(5/31/31/31/35/31/31/31/35/3){\boldsymbol Q}{\boldsymbol Q}^\text{t}=\begin{pmatrix}5/3&1/3&-1/3\\1/3&5/3&1/3\\-1/3&1/3&5/3\end{pmatrix}

となります(うーん、確かにC{\boldsymbol C}とはちょっと違うなぁ。残念)。

以上で相関行列の分解に特異値分解を用いる理由とその方法の解説を終了します。