Introduction
この世に生を受けて30年以上たつけど、もう5万回くらいSVDのやりかたを覚えて忘れてまた覚え直すを繰り返している。
肌身に染み込んで行かないのが文系ブレインの限界ですよね...
ゆえに、SVDの概要・用途・実装方法をここに記す。
[参考文献]
①：千葉大羽石研究室のサイト（授業「医用画像工学」のスライド）
②：小野滋さんのWeb上の公開資料（Greenacre1984の和訳・解説）
③：ボレガラ・岡崎・前原, 2016, 「ウェブデータの機械学習」,pp.158-162.
SVDの概要
基本の定義
任意のm × nの実数行列があるときに、次元m×mの直交行列
と、次元n× nの直交行列
、
そして非対角成分がゼロで対角成分が、
となる行列が存在し（この非ゼロの対角成分を 特異値とよぶ）
と変形できる。
便利な別表現
とすれば、以下のようにも、表現できる。
基本性質
前提知識：フロベニウスノルムによる距離の定義
をフロベニウスノルムという。
ちなみに2つめの等号においては、次元数の等しい行列の成分の積の総和は、
のトレースに等しいという性質（参考URL）を利用している。
（ なので
）
このフロベニウスノルムを使うと、二つの行列の距離を
として定義できる。
性質①：最小二乗近似を得る方法としてのSVD
行列のSVDが得られているとき、降順に並んだ
個の特異値と、それに対応する左右の特異ベクトルからなる
は、
すべてのランクの行列
のなかで、距離
を最小化するような行列となる。
すなわち、この意味では、
の「ランク
の最小二乗近似」であるといえる。
性質②：行列
と転置行列
との積の成分分解
SVDにより、が得られているとき、以下が自明に導かれる。
ここで、と表現することとすると
となり、や
が分析対象として意味のある行列のときは、
それぞれを左右の特異ベクトルと特異値
の二乗を用いてランクkで表すことができる。
性質②の具体的な例
「や
が分析対象として意味のある行列のとき」とかいきなり言われてもピンとこないと思うので、一つ例をあげよう。
たとえばm人分のケースについてnつの変数の観測値を保存した行列があり、
そこから各列について、変数の平均値を差し引いて中心化した行列
を得る。
すると、ランクn×nの行列は、変数の分散共分散行列（にスカラmをかけたもの）になるし、
いっぽうで、ランクm×mの行列は個人間の距離（にスカラnをかけたもの）を表した行列（対角成分は重心からの距離）になる。
したがって、についてSVDを適用しておけば、
変数のばらつき（分散・共分散）や、個人間の距離についての全ての情報を、k個のベクトルの線形和として表現できる。
（厳密にはもっと説明が必要だけど）これが、いわゆる主成分分析の大枠である。
カテゴリカルな変数に対して用いられる「対応分析」や「多重対応分析」と言われる手法も、の
をどう選ぶかだけの違いに帰着して説明することができる。
（詳しくは別稿で）
その意味で、SVDは多くの分析手法を統合的に理解するうえで非常に強力な方法・枠組みだといえよう。
性質③：2つの特異ベクトル間の変換
もう一度、SVDの別表現を記載する。
この式の両辺に、右特異ベクトルをかけると、特異ベクトルの正規直交性から、以下のようになる。
ゆえに、
となり、左特異ベクトルを右特異ベクトルとSVDの適用対象行列、そして特異値によって表現することができる。
また、上の別表現を転置したのちに右からをかけることで、同様に
を得られる。
実装方法
以下の手順を実装すればいい。Simple!
あるいは
（通常は計算量が少ない方、すなわちm < nだったら
と右あるいは左特異行列
を得る。
- 上記の性質③の変換式を利用して、
ちなみにあるいは
は対称行列なので、いくつかの効率的な固有値・固有値分解計算の効率的な方法（ヤコビ法、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!
