Pythonプログラミング(ステップ7・配列・線形代数・行列の分解)

このページでは、NumPyライブラリの線形演算機能を用いる例として、各種の行列の分解についてまとめておく。

行列を行列の積に分解する代表的な操作

LU分解とコレスキー分解

連立一次方程式 $$ A \boldsymbol{x} = \boldsymbol{b} $$ を解く際に、もしも行列$A$が、下三角行列$L$と上三角行列$U$の積 $$ A = L U $$ に分解できるとすると(これはLU分解と呼ばれる)、解くべき方程式は $$ L U \boldsymbol{x} = \boldsymbol{b} $$ となるから、 $$ \boldsymbol{y} = U \boldsymbol{x} $$ と未知変数を置き直せば、まず、 $$ L \boldsymbol{y} = \boldsymbol{b} $$ の解を求め、次に $$ U \boldsymbol{x} = \boldsymbol{y} $$ を解けば、最終的に解 $\boldsymbol{x}$ を求めることができる。

問題を余計に面倒にしただけのようにも思えるが、三角行列で与えられた連立方程式は、掃き出し法でとても簡単に解くことができるので、 一旦LU分解が出来たら、異なる$\boldsymbol{b}$に対して解を得るのは容易である。

また、三角行列の対角成分はその固有値が並んでいるので、その行列式は対角成分の積に等しいことから、$|A| = | LU | = |L| |U|$ によって、LU分解を経由して$A$の行列式を求めることができる。

$A$が正定値対称行列(固有値が全て正であるような対称行列)ならば、下三角行列とその転置のみを使って $$ A = L L^\top $$ と一意に分解できる。これは、特にコレスキー分解と呼ばれている。

QR分解

任意の正方行列$A$は、直交行列$Q$と、上(右)三角行列$R$の積 $$ A = QR $$ に分解することができる(QR分解)。 $A$の基底と、直交基底が、三角行列によって媒介される格好の式になっている。

行列$A$を列ベクトルの並び $$ A = \left( \boldsymbol{a}_1 \boldsymbol{a}_2 \cdots \boldsymbol{a}_n \right) $$ とみなして、これらのベクトルからグラム・シュミットの方法で、正規直交系を構成する手順を考える。

$$ \begin{eqnarray} \boldsymbol{q}_1 & = & \frac{\boldsymbol{a}_1}{\| \boldsymbol{a}_1 \|} \\ \boldsymbol{v}_2 & = & \boldsymbol{a}_2 - (a_2 \cdot q_1) q_1, \;\; \boldsymbol{q}_2 = \frac{\boldsymbol{v}_2}{\| \boldsymbol{v}_2 \|} \\ \vdots & &\\ \boldsymbol{v}_j &=& \boldsymbol{a}_j - (a_j \cdot q_1) q_1 - (a_j \cdot q_2) q_2 - \cdots - (a_j \cdot q_{j-1}) q_{j-1}, \;\; \boldsymbol{q}_j = \frac{\boldsymbol{v}_j}{\| \boldsymbol{v}_j \|} \\ \end{eqnarray} $$

同じ式を並べ替えて $$ \begin{eqnarray} \boldsymbol{a}_1 & = & \| \boldsymbol{a}_1 \| \boldsymbol{q}_1 \\ \boldsymbol{a}_2 & = & (\boldsymbol{a}_2 \cdot \boldsymbol{q}_1) \boldsymbol{q}_1 + \| \boldsymbol{v}_2 \| \boldsymbol{q}_2 \\ \vdots & & \\ \boldsymbol{a}_j & = & (a_j \cdot q_1) q_1 - (a_j \cdot q_2) q_2 - \cdots - (a_j \cdot q_{j-1}) q_{j-1} + \| \boldsymbol{v}_j \| \boldsymbol{q}_j \\ \end{eqnarray} $$ これを行列の形 $$ \left( \boldsymbol{a}_1 \boldsymbol{a}_2 \cdots \boldsymbol{a}_n \right) = \begin{pmatrix} \boldsymbol{q}_1 & \boldsymbol{q}_2 & \cdots & \boldsymbol{q}_n \\ \end{pmatrix} \begin{pmatrix} \| \boldsymbol{a}_1 \| & \boldsymbol{a}_2 \cdot \boldsymbol{q}_1 & \cdots & \boldsymbol{a}_n \cdot \boldsymbol{q}_1 & \\ 0 & \| \boldsymbol{v}_2 \| & & \boldsymbol{a}_n \cdot \boldsymbol{q}_2 & \\ \vdots \\ 0 & 0 & \cdots & \| \boldsymbol{v}_n \| \\ \end{pmatrix} $$ にまとめてみると、グラムシュミットの直交化法によって、$A$が直交行列と上三角行列の積に分解されていることが確認できる。

対角行列を含む分解(対角化)

データ処理にしばしば登場する相関行列や、電気回路、制御、化学プロセス、物理現象等に関係してしばしば登場する行列の多くは対称行列$S$である。 実対称行列の固有値は実数で、その固有値を対角要素に並べた対角行列$\Lambda$と、直行行列$P$を使って $$ S = P \Lambda P^\top $$ のように分解することができる。 ここで、$P^{-1} = P^\top$である。

実対称行列を複素数に拡張したエルミート行列$U$は $$ U = P \Lambda P^* $$ のように分解できる。 ここで、$P^*$は$P$のエルミート共役(転置を行い、要素の複素共役を取った行列)である。

$n$次の正方行列$A$を構成する$n$個の固有ベクトル全てが互いに独立ならば、 固有ベクトル(列ベクトル)を並べて構成した行列$P$を使って、 $$ A = P \Lambda P^{-1} $$ と分解できる。 ここでも$\Lambda$は、対角要素に固有値を並べた対角行列である。

特異値分解

任意の$m$行$n$列の行列$A$は、$m$行$m$列のユニタリー行列 $U$ と $n$行$n$列のユニタリー行列$V$のエルミート共役$V^*$を使って $$ A = U \Sigma V^* $$ のように分解できる。 ここで、行列$A$から構成される列(あるいは行)ベクトルの独立な成分の数(ランク)を$r$とすれば、 $\Sigma$は、サイズ$r$の対角ブロックと、0要素で構成されるブロックで構成される$m$行$n$列の行列である: $$ \Sigma = \left( \begin{array}{@{}c|c@{}} \begin{matrix} \lambda_1 & 0 & & 0 & \\ 0 & \lambda_2 & & 0 & \\ & & \ddots & & \\ 0 & 0 & & \lambda_r & \\ \end{matrix} & O \\ \hline O & O \\ \end{array} \right) $$

ユニタリー行列(直交行列)$U$と$V$はベクトルの回転や反転、対角行列$\Sigma$は伸長や短縮を表すことを思い出すと、任意の線形変換は、回転・反転の後に、軸に沿った伸縮操作を施し、再び回転・反転を行うという3ステップに分解できる、とも解釈できる。

NumPyの機能

NumPyの線形代数ライブラリ(numpy.linalg)には、行列の分解がいくつか実装されている。コードの冒頭で

import numpy as np

が宣言されていることを前提に、表にまとめておく。

関数 機能
L = np.linalg.cholesky(A)
対称行列Aをコレスキー分解した行列を返す
Q,R = np.linalg.qr(A)
AのQR分解
w,P = np.linalg.eigh(U) 
対称行列またはエルミート行列の固有値の一次元配列wと
固有ベクトルを並べた行列(すなわち変換行列)Pを返す。
†固有ベクトルは縦(列)ベクトルとして返される
w,P = np.linalg.eig(A)
正方行列の固有値wと対応する固有ベクトルを並べた行列Pを返す
U,s,V = np.linalg.svd(A,full_matrices=False)
行列を特異値分解した結果を返す。
sは対角要素の一次元配列。
ここで返されるVは$V^*$。

対角成分が一次元配列sで与えられているとき

np.diag(s)

で(非対角要素の0が補われた)対角行列を得ることができる。