SVD | 制御系CAD

制御系CADに関わる数値計算では実行列の特異値分解(Singular-Value Decomposition, SVD)を多用しますが、データ駆動分野の計算でも基盤技術となっています。

●いまデータは次の複素行列にまとめられているとします。

$\displaystyle{(1)\quad X=\left[\begin{array}{ccccc} x_1 & x_2 & \cdots & x_m \end{array}\right]\in{\bf C}^{n\times m} \quad(x_1,x_2,\cdots,x_n\in{\bf C}^n) }$

ここで列ベクトル $x_1,x_2,\cdots,x_n$ はn次元数ベクトルですが、データ駆動分野では、nの値は数個ではなく、数千、数万のオーダーになります。

● $X$ をサイズ $n\times m$ の複素行列（ ${\rm rank}\,X=k$ ）とします。このとき、ユニタリ行列 $U$ と $V$ が存在して

$\displaystyle{(2.1)\quad X=U\Sigma V^H\quad(U\in{\bf C}^{n\times n},\Sigma\in{\bf R}^{n\times m},V\in{\bf C}^{m\times m}) }$

が成り立ちます（ $V^H=V^{*T}$ ）。ここで、実行列 $\Sigma$ は次のように表されます。

$\displaystyle{(2.2)\quad \Sigma= \left[\begin{array}{cc} \widehat{\Sigma} & 0_{k\times m-k}\\ 0_{n-k\times k} & 0_{n-k\times m-k} \end{array}\right]\in{\bf R}^{n\times m} }$

$\displaystyle{(2.3)\quad \widehat{\Sigma} ={\rm diag}\{\sigma_1,\cdots,\sigma_k\}\quad(\sigma_1\ge\cdots\ge\sigma_k>0) }$

MATLABのコマンドでは、次のように求めます。

>> [U,S,V]=svd(X)

●いま $\boxed{n\ge m}$ とします。このとき(2)は次式となります。

$\displaystyle{(3.1)\quad X= \underbrace{\left[\begin{array}{cc} U_1 & U_2 \end{array}\right]}_{U} \underbrace{\left[\begin{array}{cc} \Sigma_1 \\ 0_{n-m\times m} \end{array}\right]}_{\Sigma}V^H =U_1\Sigma_1V^H }$

$\displaystyle{(3.2)\quad \Sigma_1={\rm diag}\{\sigma_1,\cdots,\sigma_k, 0,\cdots,0\}\in{\bf R}^{m\times m} }$

MATLABのコマンドでは、次のように求めます。

>> [U1,S1,V]=svd(X,’econ’)

(2)の場合は $U$ は大規模な行列となりますが、(3)の場合はスリムな行列となりますので、Economic SVDと呼ぶことがあります。

●いま実効階数を $r$ とします。すなわち微小な $\epsilon>0$ に対して

$\displaystyle{(4)\quad \sigma_1\ge\cdots\ge\sigma_r> \epsilon >\sigma_{r+1}\ge\cdots\ge\sigma_k>0 }$

とします。このとき(3)から次式のように近似できます。

$\displaystyle{(5.1)\quad X\approx \widetilde{X}= \underbrace{\left[\begin{array}{cc} \widetilde{U} & \cdot \end{array}\right]}_{U} \underbrace{\left[\begin{array}{cc} \widetilde{\Sigma} & 0_{r\times m-r} \\ 0_{n-r\times r} & 0_{n-r\times m-r} \end{array}\right]}_{\Sigma'} \underbrace{\left[\begin{array}{c} \widetilde{V}^H\\ \cdot \end{array}\right]}_{V^H} =\widetilde{U}\widetilde{\Sigma}\widetilde{V}^H }$

$\displaystyle{(5.2)\quad \widetilde{\Sigma} ={\rm diag}\{\sigma_1,\cdots,\sigma_r\}\quad(\sigma_1\ge\cdots\ge\sigma_r>\epsilon) }$

●(2.1)、(3.1)、(5.1)の関係はダイアディック形式を用いて次のように表されます。

$\displaystyle{(6)\quad \underbrace{\sum_{i=1}^{m}\sigma_iu_iv_i^H}_{U\Sigma V^H} =\underbrace{\sum_{i=1}^{k}\sigma_iu_iv_i^H}_{U_1\Sigma_1 V^H} \approx \underbrace{\sum_{i=1}^{r}\sigma_iu_iv_i^H}_{\widetilde{U}\widetilde{\Sigma}\widetilde{V}^H} }$

●相関行列 $XX^H$ または $X^HX$ は、(3.1)より、次式のように表されます。

$\displaystyle{(7.1)\quad \begin{array}{l} XX^H=U\left[\begin{array}{cc} \Sigma_1 \\ 0_{n-m\times m} \end{array}\right] \underbrace{V^HV}_{I_m}\left[\begin{array}{cc} \Sigma_1 & 0_{m\times n-m} \end{array}\right] U^H\\=U\left[\begin{array}{cc} \Sigma_1^2 & 0_{m\times n-m} \\ 0_{n-m\times m} & 0_{n-m\times n-m} \end{array}\right] U^H \end{array} }$

$\displaystyle{(7.2)\quad X^HX=V\left[\begin{array}{cc} \Sigma_1 & 0_{m\times n-m} \end{array}\right] \underbrace{U^HU}_{I_n}\left[\begin{array}{cc} \Sigma_1 \\ 0_{n-m\times m} \end{array}\right] V^H=V\Sigma_1^2 V^H }$

これらより

$\displaystyle{(8.1)\quad \begin{array}{l} XX^HU=U\left[\begin{array}{cc} \Sigma_1^2 & 0_{m\times n-m} \\ 0_{n-m\times m} & 0_{n-m\times n-m} \end{array}\right] \end{array} }$

$\displaystyle{(8.2)\quad X^HXV=V\Sigma_1^2 }$

すなわち、非零の特異値は相関行列 $XX^H$ または $X^HX$ の固有値の正の平方根です。もし $X=X^H$ の場合は $X$ の特異値は $X$ の固有値の絶対値となります。これはスカラの場合、次式を意味しています。

$\displaystyle{(9)\quad \sigma=\sqrt{xx^*}=\sqrt{(x_R+jx_I)(x_R-jx_I)}=\sqrt{x_R^2+x_I^2}=|x| }$

データ駆動分野では、nの値は数個ではなく、数千、数万のオーダーになります。一方、ｍの値はそのように大きな値とはなりません。したがって、(8.2)の小規模な固有値問題を解いて $\Sigma_1$ と $V$ を求め、(3.1)に基づいて、(8.1)の大規模な固有値問題の固有ベクトルの一部 $U_1$ を、次式により求めることが行われます。

$\displaystyle{(10)\quad U_1=XV\Sigma_1^{-1} }$

(5.1)に基づく場合は次式となります。

$\displaystyle{(11)\quad \widetilde{U}=X\widetilde{V}\widetilde{\Sigma}^{-1} }$

このような手法はスナップショットと呼ばれています。

●観測データを行ベクトルとしてもつ $X$ を考えます。

$\displaystyle{(12)\quad X=\left[\begin{array}{c} X_{11},\cdots,X_{1m}\\ \vdots\\ %X_{i1},\cdots,X_{im}\\ %\vdots\\ X_{n1},\cdots,X_{nm}\\ \end{array}\right] }$

ここで、 $n$ 回の観測が行われ、各観測では $m$ 個のデータを得るものとしています。これから

$\displaystyle{(13)\quad \bar{X}=\left[\begin{array}{c} 1\\ \vdots\\ 1 \end{array}\right] \left[\begin{array}{c} \displaystyle{\frac{1}{n}\sum_{i=1}^nX_{i1},\cdots,\frac{1}{n}\sum_{i=1}^nX_{im}} \end{array}\right] }$

を引いた行列 $B$ の特異値分解を

$\displaystyle{(14)\quad B=X-\bar{X}=U\Sigma V^H }$

とします。このとき、共分散行列は

$\displaystyle{(15)\quad C=\frac{1}{n-1}BB^H=\frac{1}{n-1}V\Sigma^2 V^H\Rightarrow D=\frac{1}{n-1}\Sigma^2 }$

のように表されます。すなわち共分散行列の主成分ベクトルと、その方向に沿う不変分散値が、行列 $B$ の特異値分解から得られることに留意します。