論理の流刑地

地獄の底を、爆笑しながら闊歩する

ダンゴムシでもわかる特異値分解(SVD)

無脊椎動物なめんな。節足動物なめんな。

Introduction

この世に生を受けて30年以上たつけど、もう5万回くらいSVDのやりかたを覚えて忘れてまた覚え直すを繰り返している。
肌身に染み込んで行かないのが文系ブレインの限界ですよね...

ゆえに、SVDの概要・用途・実装方法をここに記す。

[参考文献]
①:千葉大羽石研究室のサイト(授業「医用画像工学」のスライド)
②:小野滋さんのWeb上の公開資料(Greenacre1984の和訳・解説)
③:ボレガラ・岡崎・前原, 2016, 「ウェブデータの機械学習」,pp.158-162.

SVDの概要

基本の定義

任意のm × nの実数行列A (rank = k) があるときに、次元m×mの直交行列Vと、次元n× nの直交行列U
そして非対角成分がゼロで対角成分 \lambda_iが、
  \ \lambda_i > 0  \ (i <= k ), \  \lambda_i = 0 (i > k )

となる行列Sが存在し(この非ゼロの対角成分を 特異値とよぶ)

 A = V S U^T

と変形できる。

便利な別表現

 V = \begin{bmatrix} v_1 & \ldots & v_m   \end{bmatrix}
 U = \begin{bmatrix} u_1 & \ldots & u_n   \end{bmatrix}

とすれば、以下のようにも、表現できる。
 A = V S U^T = \Sigma_{i=1}^ k \lambda_i v_i u_i^T

基本性質

前提知識:フロベニウスノルムによる距離の定義

 || \boldsymbol{A}||_F = \sqrt{\Sigma_{i=1}^m \Sigma_{j=1}^n |a_{ij}|^2  } = \sqrt{ trace(A^T A) }
をフロベニウスノルムという。

ちなみに2つめの等号においては、次元数の等しい行列 A, Bの成分の積の総和は、 A B^Tのトレースに等しいという性質(参考URL)を利用している。
 trace(AB) = trace(BA) なので trace(A^T A ) = trace(A A^T)


このフロベニウスノルムを使うと、二つの行列A, Bの距離を || \boldsymbol{A - B}||_F として定義できる。

性質①:最小二乗近似を得る方法としてのSVD

行列AのSVDが得られているとき、降順に並んだj \ ( j <= k ) 個の特異値と、それに対応する左右の特異ベクトルからなる
  A_{j^*}  = \Sigma_{i = 1}^j \lambda_j u_i  v_i  は、
すべてのランクjの行列X_{j^*}のなかで、距離|| \boldsymbol{A - X_{j^*}}||_F を最小化するような行列となる。


すなわち、この意味で  A_{j^*} は、Aの「ランクjの最小二乗近似」であるといえる。

性質②:行列Aと転置行列A^Tとの積の成分分解

SVDにより、 A = V S U^Tが得られているとき、以下が自明に導かれる。


 A A^T = V S U^T (U^T)^T S^T V^T = V S U^T U S^T V^T = V S S^T V^T
 A^T A  =  (U^T)^T S^T  V^T V S U^T  = U  S^T  V^T V S U^T = U S^T S U^T

ここで、 S S^T = S_{row}^2 ; \ S^T S = S_{column}^2 と表現することとすると
 A A^T =  V S_{row}^2 V^T \iff   AA^TV = V S_{row}^2 \iff AA^Tv_i = \lambda_i^2 v_i
  A^T A  =  U S_{column}^2 U^T \iff   A^TAU = U S_{column}^2  \iff A^TAu_i = \lambda_i^2 u_i


となり、 A A^T  A^T Aが分析対象として意味のある行列のときは、
それぞれを左右の特異ベクトルv_i ,\  u_iと特異値\lambda_iの二乗を用いてランクkで表すことができる。

性質②の具体的な例

 A A^T  A^T Aが分析対象として意味のある行列のとき」とかいきなり言われてもピンとこないと思うので、一つ例をあげよう。

たとえばm人分のケースについてnつの変数の観測値を保存した行列X =\begin{bmatrix} x^1 & \ldots & x^n   \end{bmatrix} があり、
そこから各列について、変数の平均値\bar{x^j}を差し引いて中心化した行列X^*を得る。


すると、ランクn×nの行列(X^*)^TX^* は、変数の分散共分散行列(にスカラmをかけたもの)になるし、
いっぽうで、ランクm×mの行列(X^*)^TX^* は個人間の距離(にスカラnをかけたもの)を表した行列(対角成分は重心からの距離)になる。

したがって、X^*についてSVDを適用しておけば、
変数のばらつき(分散・共分散)や、個人間の距離についての全ての情報を、k個のベクトルの線形和として表現できる。

(厳密にはもっと説明が必要だけど)これが、いわゆる主成分分析の大枠である。
カテゴリカルな変数に対して用いられる「対応分析」や「多重対応分析」と言われる手法も、 A = V S U^TAをどう選ぶかだけの違いに帰着して説明することができる。
(詳しくは別稿で)
その意味で、SVDは多くの分析手法を統合的に理解するうえで非常に強力な方法・枠組みだといえよう。

性質③:2つの特異ベクトル間の変換

もう一度、SVDの別表現を記載する。
 A = \Sigma_{i=1}^ k \lambda_i v_i u_i^T

この式の両辺に、右特異ベクトル u_j \ ( 1 <= j <= k )をかけると、特異ベクトルの正規直交性から、以下のようになる。

 A u_j  = \Sigma_{i=1}^ k \lambda_i v_i u_i^T u_j   = \lambda_j v_j  u_j^T u_j = \lambda_j  v_j
ゆえに、

 v_j = \frac{1}{\lambda_j}  A u_k

となり、左特異ベクトルを右特異ベクトルとSVDの適用対象行列A、そして特異値によって表現することができる。

また、上の別表現を転置したのちに右から v_j \ ( 1 <= j <= k )をかけることで、同様に

 u_j = \frac{1}{\lambda_j}  A^T v_k

を得られる。

固有値分解との関係性

上述の性質②より、
 A A^T =  V S S^T V^T = V  S_{row}^2  V^T
 A^T A = U S^T S U^T =  U  S_{column}^2 U^T

が成立する。よく見ればこれは A A^Tおよび A^T A固有値分解である。
よって特異値\lambda_iというのは、 A A^T A^T A固有値分解したときに得られる固有値平方根である。

実装方法

以下の手順を実装すればいい。Simple!

  1.  A A^Tあるいは A^T A(通常は計算量が少ない方、すなわちm < nだったら A A^T、m > n だったら A^T A)を固有値分解し、その結果から特異値の対角行列Sと右あるいは左特異行列SもしくはUを得る。
  2. 上記の性質③の変換式を利用して、SUの得られていない方を得る。

ちなみに A A^Tあるいは A^T Aは対称行列なので、いくつかの効率的な固有値固有値分解計算の効率的な方法(ヤコビ法、3重対角行列を用いた逆反復法、ユーザーである我々は詳しく知る必要はない..)が使えるよ。

Rで確認

Rにはsvd()というSVDのための関数が備え付けであるが、一応上の方法で計算されていることを確認するために検算してみる

# 適用データ生成
# 元データは、Clausen,1998=2015「対応分析入門」, p.189より
tbl <- matrix( c( 395 , 2456 , 1758 , 147 , 153 , 916 , 694 , 327 , 1347)  , byrow=3  , ncol = 3)
colnames(tbl) <- c("強盗" , "詐欺", "破壊")
rownames(tbl) <- c( "オスロ" , "中部地域", "北部地域")

#分布に変換
tbl_p <-  tbl / sum(tbl)
mgn_c <- apply( tbl_p , 2, sum)#列周辺分布
mgn_r <- apply( tbl_p , 1, sum) #行周辺分布

D_clm <- diag(1/sqrt(mgn_c)) #列周辺分布の-1/2乗を対角成分とした行列
D_row <- diag(1/sqrt(mgn_r)) #行周辺分布の-1/2乗を対角成分とした行列

# SVDを適用する対象の行列をつくる。
mat_A <- D_clm %*% (tbl_p - mgn_r %o% mgn_c) %*% D_row


## SVDを適用する
svd_res <- svd( mat_A) 

#--結果--#
# ちなみに上記の解説とはuとvが逆になってるので、注意が必要(uが左, vが右)
#$d  (特異値)
#[1] 7.927067e-01 8.484083e-02 1.339025e-17
#
#$u(左特異ベクトル)
#          [,1]        [,2]      [,3]
#[1,] -0.9210524  0.02831857 0.3884076
#[2,]  0.2288391 -0.76764442 0.5986273
#[3,]  0.3151112  0.64025003 0.7005604
#$v (右特異ベクトル)
#          [,1]       [,2]      [,3]
#[1,]  0.1744326  0.6379811 0.7500356
#[2,] -0.8992856 -0.2070410 0.3852524
#[3,]  0.4010719 -0.7416968 0.5376125


## 実際に固有値を利用する方法でやってみるver
#今回は行列Aが正方行列なので計算量的にはどっちでもいいが、
#とりあえず、左特異ベクトルを求めてから変換公式で右特異ベクトルを求める形にする。

mat_B <- mat_A %*% t(mat_A) #固有値分解の対象行列を生成
evd_res <- eigen(mat_B, TRUE) #固有値分解する

sval <- sqrt(evd_res$values ) #平方根をとることで特異値をえる
u_vec <- evd_res$vectors #左特異ベクトルを得る

v_vec <-    t(mat_A) %*% u_vec %*% diag(  1/ sval)  #変換公式 v =  1/λ * A^t * u を適用
#補足: 1/λjを対角成分とする対角行列を右からかけることによってj列を1/λj倍している

### 固有値分解をアプローチによる結果 ###
#> sval
#[1] 7.927067e-01 8.484083e-02 8.575436e-09
#> u_vec
#           [,1]        [,2]       [,3]
#[1,]  0.9210524  0.02831857 -0.3884076
#[2,] -0.2288391 -0.76764442 -0.5986273
#[3,] -0.3151112  0.64025003 -0.7005604
#> v_vec
#           [,1]       [,2]          [,3]
#[1,] -0.1744326  0.6379811 -1.941982e-08
#[2,]  0.8992856 -0.2070410 -6.473275e-09
#[3,] -0.4010719 -0.7416968  2.427478e-08
 
#--この場合、二軸で分散を説明しきっているので、3個目/3列目の特異値/特異ベクトルは誤差というか無視して良い
#--2つ目までのの特異値と特異ベクトルは先ほどのsvd()関数の戻り値と完全に一致しており、求める結果が得られた
#--なおひとつめの特異ベクトルは符号が違うが、左右一緒に同一固有値に対応する特異ベクトルの符号を変えさえすれば-1をかけてもよい(正規直交基底の性質と変換公式を考えれば明らか)

Conclusion

さいきんは、ずっと、これをきいている...
www.youtube.com


SVDがとても便利だってわかって、さらにその仕組みがわかったので、
次はSVDの枠組みでMCA(多重対応分析)とCA(対応分析)を説明するYO!
SVDがわかると両手法の仕組みや違い・共通点を超エレガントに説明できるYO!

Enjoy!