SVD -

表紙(FrontPage) | 編集(管理者用) | 差分 | 新規作成 | 一覧 | RSS | 検索 | 更新履歴

SVD 特異値分解

データを表す行列を A とする。A はほとんどの場合正方行列ではない。また対称行列でもない。特異値分解では A をもとにして $AA^T$ 、 $A^TA$ 　という行列を作ることが分解の肝になっている。そうすると、 $AA^T$ 、 $A^TA$ 　は A の性質を引き継いでいて、正方行列で、しかも対称行列にもなる。こうすると扱いやすくなる。

具体例を構成するために、R 言語で行列の計算をしてみる。

 library(MASS)
 gyou <- 2
 retu <- 3
 mat1 <- matrix(0, gyou, retu)
 mat1[1, 1:3] <- c(1, 4, 7)
 mat1[2, 1:3] <- c(2, 6, 3)

こうすると mat1 の値は

 　     [,1] [,2] [,3]
 [1,]    1    4    7
 [2,]    2    6    3

この行列は「3変数（3次元）のデータが2つある」と解釈できる。

$AA^T$ に相当するものは　mat1 %*% t(mat1)　この値は

      [,1] [,2]
 [1,]   66   47
 [2,]   47   49

確かに対称な正方行列ができている。できた正方行列の行、列の数は、元の行列の行の数になる。

$AA^T$ 、 $A^TA$ 　は、複数の行列を組み合わせて表す（分解する）ことができる。これらの行列は正方行列なので「固有値分解」を用いて分解ができる。それらの結果を用いると A 自体の分解もできて、それが特異値分解になる。

固有値分解は正方行列を固有値と固有ベクトルを組み合わせて表す（分解する）方法である。そもそも固有値、固有ベクトルとはなにか。単純に考えると、これらは、対象とする行列をうまく分解することを可能にする、「数値」と「ベクトル」の組み合わせというように考えられる。固有値は必ず固有ベクトルとセットになって計算される。「値」というから数値が一つあるだけなのかというとそうではなく、元の行列の列の数と同じだけの個数の、固有値と固有ベクトルのセットが生じる。

行列 A があるとする。固有ベクトルを B、固有値を λ とすると、

$A \small %*% \Large B = \lambda \small %*% \Large B$ 　という関係が成り立つように、それらの値が決められる (%*% は行列、ベクトルのかけ算）。手で計算するのはとても大変だが、R 言語などを使うと簡単に計算できる。左辺では行列 A だったところが、右辺では固有値 λ に置き換わっている。このことから、固有値は元の行列の性質を保持しながら縮約した、大切な部分を取り出したようなものと考えることもできる。実際に様々な分野で、固有値を求めることがとても重要な意味を持つ数値を得ることにつながることがよく知られている。

上で書いた具体例で計算してみる。固有値、固有ベクトルを R で計算するには eigen() という関数を使う。上の例の、 $AA^T$ に相当する　mat1 %*% t(mat1)　を計算すると　eigen((mat1 %*% t(mat1))　で計算できて、

 $values
 [1] 105.262433   9.737567
 
 $vectors
           [,1]       [,2]
 [1,] -0.7674517  0.6411068
 [2,] -0.6411068 -0.7674517

になる。

固有値分解では、固有値 λ、固有ベクトル B を用いて、元の行列 A を以下のように分解する。

固有ベクトルは固有値の数と同じだけの個数があり、1番、2番・・と順番がある。それらの縦のベクトルを左から並べて行列を作る。それを P とする。固有値も、1番、2番・・と順番がある。それらを対角成分にした行列を作る。それを Λ とする。

そうすると、元の行列 A は、 $A = P \Lambda P^{-1}$ 　という、三個の行列をかけ算した形に分解できる。 $P^{-1}$ 　は、行列 P の逆行列である。これを固有値分解という。A は正方行列でないといけない。

これも上の例で行ってみると、固有値を並べた対角行列を作るのに diag(), 固有ベクトルの逆行列をつくるのに solve() を使い、

 i <- eigen((mat1 %*% t(mat1))
 i$vectors %*% diag(i$values) %*% solve(i$vectors)
     [,1] [,2]
 [1,]   66   47
 [2,]   47   49

確かに、元の mat1 %*% t(mat1) と同じ値になっている。

この「固有値分解」については、「物理数学の直観的方法第 2 版」(長沼伸一郎著）という本で解説されている。

上の例では2行3列の横に長い行列にしたが、3行2列のものでもやってみる。

 library(MASS)
 gyou <- 3
 retu <- 2
 mat1 <- matrix(0, gyou, retu)
 mat1[1, 1:2] <- c(1, 4)
 mat1[2, 1:2] <- c(2, 6)
 mat1[3, 1:2] <- c(9, 3)

上の例の、 $AA^T$ に相当する　mat1 %*% t(mat1)　を計算すると　

      [,1] [,2] [,3]
 [1,]   17   26   21
 [2,]   26   40   36
 [3,]   21   36   90

eigen((mat1 %*% t(mat1))　は

 eigen() decomposition
 $values
 [1]  1.182800e+02  2.871998e+01 -3.250408e-15
 
 $vectors
            [,1]       [,2]        [,3]
 [1,] -0.2943417 -0.4848891  0.82355662
 [2,] -0.4782796 -0.6713208 -0.56619518
 [3,] -0.8274126  0.5605452  0.03431486

固有値を並べた対角行列を作るのに diag(), 固有ベクトルの逆行列をつくるのに solve() を使い、

 i <- eigen((mat1 %*% t(mat1))
 i$vectors %*% diag(i$values) %*% solve(i$vectors)
 
      [,1] [,2] [,3]
 [1,]   17   26   21
 [2,]   26   40   36
 [3,]   21   36   90

この場合は行の数が3なので、3行3列の行列になり、確かに、元の mat1 %*% t(mat1) と同じ値になっている。

$A$ 　を特異値分解する場合も、行列を組み合わせて表す。

まず A は　 $A = U \Sigma V^{T}$ 　と分解できるということにする。形は固有値分解と似ている。

U と V は、直交行列と呼ばれる性質を持つということにする（後でそれが有効に使える）。

$\Sigma$ 　は、固有値分解と同じく対角行列ということにする（後でそれが有効に使える）。

これらの仮定が成り立つような U, V, $\Sigma$ 　が常に存在することは、数学的に示されているのだろう。そうでないと困る。

$A = U \Sigma V^{T}$ 　だから、 $AA^T = [ U \Sigma V^{T} ] [ U \Sigma V^{T} ]^{T}$ 　と、書くことができる。

右側の転置の部分は、転置の性質から、順番が入れ替わって、 $V \Sigma^{T} U^{T}$ 　になる。

上に書いたように、U と V は、直交行列と呼ばれる性質を持つということにしていた。この場合逆行列と転置行列が同じものになる。　 $U^T = U^{-1}$ 　 $V^T = V^{-1}$

だから　 $AA^T = [ U \Sigma V^{T} ] [ V \Sigma^{T} U^{T} ] = U \Sigma \Sigma^{T} U^{-1}$ 　　

上に書いたように、 $\Sigma$ 　は、固有値分解と同じく対角行列ということにしていたので、転置しても変化がない。だから $\Sigma \Sigma^{T} = \Sigma^{2}$ 　になる。

だから、 $AA^T = U \Sigma^{2} U^{-1}$ 　になる。これは　 $AA^T$ 　を固有値分解した形　 $P \Lambda P^{-1}$ 　に相当する。このことから、U を決定するには、　 $AA^T$ 　を固有値分解すればよいことになる。つまり、A をもとにして $AA^T$ 　という行列を作ると、 $AA^T$ 　は A の性質を引き継いでいて、正方行列になる。正方行列なので固有値分解できて、U は、 $AA^T$ 　の固有ベクトルを並べて作った行列になる。U を使って A 自体の分解（特異値分解）ができる。

$\Sigma^2$ 　は、 $AA^T$ の固有値を対角成分に並べたものになる。自乗がそういう値になるので、 $\Sigma$ 　は、それらの固有値の平方根の絶対値（正の値）ということにして、それを特異値という。これも A 自体の分解（特異値分解）に使われる。

同様に　 $A^TA$ 　という行列も作る。 V は、 $A^TA$ 　の固有ベクトルを並べて作った行列になる。この場合、V の列の数を U の行の数と一致させないといけない。そうしないと　 $V{^T}$ 　を右からかけ算するときに、行列のかけ算ができない。

上で「 $\Sigma^2$ 　は、 $AA^T$ の固有値を対角成分に並べたものになる。」と書いた。それだけでなく、 $A^TA$ 　の固有値を対角成分に並べた物も　 $\Sigma^2$ 　になる。つまり　 $\Sigma^2$ 　を計算する方法は二通りあることになる。

実際に計算してみると、片方の　 $\Sigma^2$ 　と、もう一方の　 $\Sigma^2$ 　は、計算（行列とのかけ算）に必要な部分は全く同じ値が並んでいて、どちらか一方は行列とのかけ算に使えない余計な固有値がさらに付いている。だから二通りの　 $\Sigma^2$ 　で、そのまま行列とのかけ算に使える方を採用すればよいことになる。

これで　 $U, \Sigma, V$ 　をそろえることができ、組み合わせることで A の特異値分解ができる。

分解で生じた U は、元の行列 A と同じ行数を持っていて、元の行列が含む情報を縮約している。例えば、元の行列 A の１行目のデータと３行目のデータに相関があるなら、行列 U でも１行目のデータと３行目のデータに相関がある。

分解で生じた V は、この文章でも後で出てくるが、元のデータ行列の分散共分散行列（データを規準化すれば相関係数）を計算することに使うことができる。

上の2行3列の mat1 の例で計算してみると、

i <- eigen( mat1 %*% t(mat1) )　として、i は

 $values
 [1] 105.262433   9.737567
 $vectors
           [,1]       [,2]
 [1,] -0.7674517  0.6411068
 [2,] -0.6411068 -0.7674517

になる。

これから　 $U$ は i$vectors になる。 $\Sigma$ は diag(sqrt(i$values)) になる。

つぎに j <- eigen( t(mat1) %*% mat1 )　として、j は

 $values
 [1]  1.052624e+02  9.737567e+00 -1.388973e-16
 $vectors
           [,1]       [,2]       [,3]
 [1,] -0.1997773 -0.2864267  0.9370426
 [2,] -0.6741342 -0.6538304 -0.3435823
 [3,] -0.7110781  0.7003324  0.0624695

になる。

こちらで出てくる固有値は、最初の二つは　eigen( mat1 %*% t(mat1) )　の方の固有値と一致している。もう一つの固有値は行列 U とのかけ算に使えない。だから　eigen( mat1 %*% t(mat1) )　の方の固有値を　 $\Sigma$ 　に使う。

行列とのかけ算ができるように、 $V$ は j$vectors （3行3列）の 1, 2 列目（U の行数）を取って　j$vectors[ , 1:2] にする。

元の mat1 の値は　i$vectors %*% diag( sqrt(i$values) ) %*% t( j$vectors[ , 1:2]) )

になって、計算すると

     [,1] [,2] [,3]
 [1,]    1    4    7
 [2,]    2    6    3

で、確かに元の値と一致している。

上の3行2列の mat1 の例で計算してみると、

i <- eigen( mat1 %*% t(mat1) )　として、i は

 $values
 [1]  1.182800e+02  2.871998e+01 -3.250408e-15
 
 $vectors
            [,1]       [,2]        [,3]
 [1,] -0.2943417 -0.4848891  0.82355662
 [2,] -0.4782796 -0.6713208 -0.56619518
 [3,] -0.8274126  0.5605452  0.03431486

これから　 $U$ は i$vectors になるが、今度はこちらが列が多すぎて行列のかけ算ができない。だから　i$vectors[,1:2]　にする。

つぎに j <- eigen( t(mat1) %*% mat1 )　として、j は

 $values
 [1] 118.28002  28.71998
 
 $vectors
            [,1]       [,2]
 [1,] -0.7997319  0.6003572
 [2,] -0.6003572 -0.7997319

になる。

この場合　 $V$ は j$vectors （2行2列）でそのまま行列とのかけ算ができる。

また j$values の方の固有値を　 $\Sigma$ 　に使うことで行列とのかけ算ができる。

元の mat1 の値は　i$vectors[,1:2] %*% diag( sqrt(j$values) ) %*% t( j$vectors )

になって、計算すると、

      [,1] [,2]
 [1,]    1    4
 [2,]    2    6
 [3,]    9    3

で、確かに元の値と一致している。

一つ問題がある。2行3列の mat1 の場合、この計算では U は i$vectors で

 $vectors
           [,1]       [,2]
 [1,] -0.7674517  0.6411068
 [2,] -0.6411068 -0.7674517

V は j$vectors[, 1:2] で

 $vectors
           [,1]       [,2]      
 [1,] -0.1997773 -0.2864267  
 [2,] -0.6741342 -0.6538304 
 [3,] -0.7110781  0.7003324

しかし R言語の svd() 関数では、svd(mat1) は

 $d
 [1] 10.259748  3.120508
 $u
           [,1]       [,2]
 [1,] -0.7674517 -0.6411068
 [2,] -0.6411068  0.7674517
 $v
           [,1]       [,2]
 [1,] -0.1997773  0.2864267
 [2,] -0.6741342  0.6538304
 [3,] -0.7110781 -0.7003324

で、正負が入れ違っているところがある。どこかに計算法の違いがあるのだろう。しかし元の mat1 が再生できるように分解されていることは変わらない。

特異値分解と主成分分析 PCA

特異値分解は主成分分析と関係があり、上の文章で出てきた $\Sigma^2$ 　や $V$ 　を使って計算することができる。

元の行列として3行2列のものを使ってみる。これは「二変数（二次元）のデータが、3点存在する」と解釈することができる。

 library(MASS)
 gyou <- 3
 retu <- 2
 mat1 <- matrix(0, gyou, retu)
 mat1[1, 1:2] <- c(1, 4)
 mat1[2, 1:2] <- c(2, 6)
 mat1[3, 1:2] <- c(9, 3)

まず、この行列を「列ごとの平均 = 0」「列ごとの分散 = 1」になるように変換しておく。これは、R 言語では scale(mat1) で行うことができる。この変換した行列を　 $A$ 　と書くことにする。

この変換した行列の、分散共分散行列を計算する。なぜそうするかというと、主成分分析では処理対象になるデータの間でなるべく分散が大きくなるように、変数をうまく混ぜ合わせることによってデータを変換（射影）する。そのことから分散がとても大切な値になるので、その計算をすることが必須になる。

変換した行列　 $A$ 　の分散共分散行列は、 $\frac{1}{[N - 1]} A^T A$ 　という数式で表現できる。 $A^T A$ 　は、上に書いた SVD でも出てきた。このことから主成分分析と SVD に関係があると言うこともできる。

この場合の　 $N$ 　は行列の行の数になる。 N から 1 を引くのは、不偏分散を求める（標本分散では N で割る。しかし母集団からサンプルとして採った標本から、元の母集団の分散を求めるには、N-1 で割らないといけない)からで、R 言語の VAR() も N-1 で割り算した不偏分散の値になっている。そこでここでもそれに合わせる。

 なぜ N-1 が適切かと言うことに関しては、「高校数学の美しい物語」というすばらしい
 ページに証明が記載されている。

この式だけを見てもどうしてこれが分散共分散になるのかすぐにはわからないが、3行2列なら手で紙に書いてみる、R 言語で確かめることがすぐにできるのでやってみると、「確かにそうだ」ということがわかる。

一つの列に対する不偏分散は　 $\frac{\Sigma [X_i \ - \ \hat{ X } ]^2 }{[N - 1]}$ 　で計算できる。 $\hat{ X }$ 　は、その列の平均値で、N は行の数（＝データの個数）　　不偏分散＝その標本の元になっている、想定される母集団が示す分散の値（の期待値）

最初にそれぞれの列の平均値を 0 にするように scale をしてあるので、その列の平均値　 $\hat{ X }$ 　は省略できる。そうすると、例えば一列目のデータの分散　 $\sigma_{11}$ 　は　 $\frac{X_{11}^2 + X_{21}^2 + X_{31}^2 }{[3 - 1]}$ 　になる。

共分散の場合は2乗ではなく二つの X の積になり、この場合も平均値を 0 にしてあるので、一列目と二列目の共分散　 $\sigma_{12}$ 　は　 $\frac{X_{11} X_{12} + X_{21} X_{22} + X_{31} X_{32} }{[3 - 1]}$ 　となる。これらは結局行列のかけ算で表せて、 $\frac{1}{[N - 1]} A^T A$ 　で表現できる。N は、A の行の数　　

SVD でも同じような行列のかけ算がでてきた。SVD では、行列 A は　 $A = U \Sigma V^{T}$ 　と分解できるということにした。この場合、V を導くために　 $A^T A$ 　を考えた。 $A^T A = [ U \Sigma V^{T} ]^{T} [ U \Sigma V^{T} ]$ 　と、書くことができる。

右側の転置の部分は、転置の性質から、順番が入れ替わって、 $V \Sigma^{T} U^{T}$ 　になる。上に書いたように、U と V は、直交行列と呼ばれる性質を持つということにしていた。この場合逆行列と転置行列が同じものになる。　 $U^T = U^{-1}$ 　 $V^T = V^{-1}$

また上に書いたように、 $\Sigma$ 　は、固有値分解と同じく対角行列ということにしていたので、転置しても変化がない。だから $\Sigma \Sigma^{T} = \Sigma^{2}$ 　になる。

だから、 $A^T A = V \Sigma^{2} V^{-1}$ 　になる。

そして、 $A$ 　の分散共分散行列 VAR(A) は　 $\frac{1}{[N - 1]} A^T A$ 　　N は、A の行の数

だったから、 $VAR[A]$ 　は　 $\frac{1}{[N - 1]} V \Sigma^{2} V^{-1}$ 　で計算できて、V と　 $\Sigma$ 　は SVD の場合と同じ方法で、A から計算することができる。ここにおいても、主成分分析（で重要な分散共分散行列）と SVD との関連が出てくる。

次に、主成分分析でとても重要な部分である「データの変換（射影）」について考える。

射影とはどういうことか？　この場合、「複数の変数（次元）をもつデータ点を、それぞれの変数をうまく混ぜ合わせる（例えば変数が x, y なら kx + ly という値を作る）ことで、数直線上の特定の位置（値）に変換する」ということになる。変数が 3 つ x, y, z なら kx + ly + mz というように、もっとたくさん変数がある場合も同じように変換する。

今の例では 3行2列で、2変数（二次元）のデータが3つある。それらについて、kx + ly と言う値を作る。これによって、二次元の3つのデータが、一次元の3つのデータに変換されることになる。これは変数がもっと多い場合も同じで、主成分分析が変数が多い高次元のデータを低次元に圧縮する機能を実現する元になっている。

化学の分野には「混成軌道」という考え方がある。s 軌道と p 軌道を足し合わせることで混ぜ合わせて、分子を構成する新しい軌道を作る。それと似ている。

二次元の場合 (x, y) が kx + ly に変換されることが射影になるが、これは (k, l) というベクトルを考えると、(x, y) と (k, l) の内積を計算するのと同じである。もっと変数が増えても同じように内積を計算すればよくて便利なので、そういうことにする。

今分析している行列（3データ2変数）　 $A$ 　を射影した値を S、射影変換に使うベクトルを a とすると、

$S = A \cdot a$ 　ということになる。　 $\cdot$ 　は内積

実際には、a はベクトルが元のデータ行列の変数の数だけ並んだ行列になり、変換もそれらのベクトル一つずつを用いて行うことになる。それによって第一主成分、第二主成分、･･･というように値が出てくる。

主成分分析では、変換に使うベクトル（を並べた行列） a を決めるために、「射影した結果 S の分散が、最も大きくなるように a を決める」という基準を用いる。

S の分散は、上で A の分散を求めたときと同じように考えることができるので、

$VAR(S) = \frac{1}{[N - 1]} S^T S$ 　で、これに　 $S = A \cdot a$ 　を代入すると、

$VAR(S) = \frac{1}{[N - 1]} S^T S = \frac{1}{[N - 1]} [A \cdot a]^T [A \cdot a] = \frac{1}{[N - 1]} a^T A^T A a$ 　

うまい具合に　 $\frac{1}{[N - 1]} A^T A$ 　の部分は VAR(A) と一致するので、書き直すと

$VAR(S) = a^T VAR(A) a$ 　になる。

この VAR(S) が最大になるように a を決めればよい。それには、未定乗数法という、こういう計算によく使われるテクニックを使う。このことについては優れた解説がたくさん公開されているので、それらを今回の場合に合わせて自分にとってわかりやすいように長ったらしく書き直してみる。

未定乗数法

未定乗数法を使うには、まず最大・最小にしたい値を決める。しかしそれだけだと、「その値が無限に大きくなる」とか「無限に小さくなる」条件が答えになってしまうことが多くなる。そこで、「ベクトル a の長さ＝１」とか、「(x, y) は x^2 + y^2 = 1 を満たさないといけない」のような条件、制限をつける必要がある。そういった条件を拘束条件という。

今回の例では、「a はベクトルなので長さがあるが、それは常に 1 」という制限を付ける。だから、まずベクトル a の長さを求めないといけない。

まず簡単な例として、二次元の場合を考える。二次元の平面で原点から (a1, a2) へ向かうベクトルを考えると、その長さは　 $\sqrt{ a1^2 + a2^2 }$ 　で計算できる。しかし平方根がつくと多次元を考えた際に計算がやりにくくなる。今の例では、ベクトル a の長さを 1 に制限する。だから a の長さの二乗も 1 になる。これを条件にした方がやりやすい。だからそうすると　 $a1^2 + a2^2 = 1$ 　が常に成り立つという制限にしてもよい。今回はそうする。

a が 3 次元になると、2 次元の三角形の (a1, a2) の頂点から、a1, a2 の軸と直交する方向へ一つ軸を伸ばすことになる。そうすると、一つの辺の長さは SQRT(a1^2 + a2^2)、もう一つの辺の長さは a3 という新しい三角形ができる。その新しい三角形の、もう一つの辺の長さが必要な距離で、SQRT(a1^2 + a2^2 + a3^2)　になり、2 乗すると a1^2 + a2^2 + a3^2 になる。もっと次元が増えた多次元なら、同じことを繰り返すことになり a1^2 + a2^2 + … = 1　ということになる。都合のよいことに、これは行列のかけ算で表現できて、

$a^T a = 1$ 　と書くことができる（a は縦に数値が並んだベクトル）。この式の左辺を関数と考えると、引数は a になり、a で微分することもできる。 a は a1, a2, … と値が縦に並んだベクトル（またはベクトルがさらに並んだ行列）なので、それらの項一つずつで微分する偏微分をすることができる。

a が二次元のベクトルの場合で微分してみると、a1^2 + a2^2 = 1 で、左辺に二つの項がある。左辺を a1 で微分（a2 は定数とするので偏微分）すると 2 * a1 になる。a2 で偏微分すると 2 * a2 になる。関数が最大・最小になる条件を求める際は「微分した値が０」という式を使うことが多いので、微分・偏微分できることはとても重要・必須になる。

今回　 $VAR(S) = a^T VAR(A) a$ 　の VAR(S) を最大にしたい。上に書いたように、関数の最大、最小になる条件を決めるには、微分してその値が 0 になるという条件を使う。行列･ベクトルの状態で考えるのは私の頭ではやりにくいので、ここでも一番簡単な、a が二次元のベクトルの場合を例題にしてみる。

$a = \left( \begin{array}{GC+23} a_1 \\ a_2 \end{array} \right)$ 　とすると、VAR(A) は 2 行 2 列の分散共分散行列になる（そうでないとかけ算ができない）。その値を　 $VAR(A) = \left( \begin{array}{GC+23} \large S_{11} \ S_{12} \\ \large S_{12} \ S_{22} \end{array} \right)$ 　とする。分散共分散行列は対称行列なので、S12 = S21 で S11, S12, S22 で表せる。

だから　 $VAR(S) = a^T VAR(A) a$ 　は　 $\left( \begin{array}{GC+23} a_1 \ a_2 \end{array} \right) \left( \begin{array}{GC+23} \large S_{11} \ S_{12} \\ \large S_{12} \ S_{22} \end{array} \right) \Large \left( \begin{array}{GC+23} a_1 \\ a_2 \end{array} \right)$ 　　と書くことができる。

これを計算すると

$\left( \begin{array}{GC+23} \large a_1 S_{11} + a_2 S_{12} , \ a_1 S_{12} + a_2 S_{22} \end{array} \right) \Large \left( \begin{array}{GC+23} a_1 \\ a_2 \end{array} \right)$ 　 $= \ S_{11} a_1^2 + 2 a_1 a_2 S_{12} + S_{22} a_2^2$ 　

という、a1 と a2 の関数になる。これらを a1, a2 でそれぞれ偏微分すると

a1 で：　 $2 S_{11} a_1 + 2 S_{12} a_2$

a2 で：　 $2 S_{22} a_2 + 2 S_{12} a_1$

ここまでが準備で、ここで未定乗数を導入した式を作る。

今回　 $VAR(S) = a^T VAR(A) a$ 　を最大にしたい。また拘束条件は　「　 $a_1^2 + a_2^2 = a^T a = 1$ 　が常に成り立つ」　だった。それらのことを利用して、以下のような式を作る。

$L = a^T VAR(A) a \ - \ \lambda[a^T a - 1]$ 　　　　　　 $\lambda$ 　は未定乗数と呼ばれる数で、拘束条件をうまく式に導入するために用いられる。

拘束条件から、 $a^T a - 1$ 　は、つねに 0 になる。だから上の式の L が最大の値を取る条件では、 $VAR(S) = a^T VAR(A) a$ 　も最大の値を取る。そこで L が最大の値を取る条件を求めるために L を a1, a2 でそれぞれ偏微分してその値が 0 とすると、

$\frac{\partial L}{\partial a_1} = 2 S_{11} a_1 + 2 S_{12} a2 - \lambda [2 a_1] = 0$

$\frac{\partial L}{\partial a_2} = 2 S_{22} a_2 + 2 S_{12} a1 - \lambda [2 a_2] = 0$ 　

と計算することができる。 $\lambda$ 　のおかげで、式に拘束条件がうまく導入できている。

うまい具合に、上の式の　 $\lambda$ 　の部分を右辺に移すと、行列を用いた書き方　 $A \small %*% \Large B \ = \ \lambda \small %*% \Large B$ 　（%*% は行列のかけ算）に書き換えることができる。書いてみると、

$\left( \begin{array}{GC+23} \large S_{11} \ S_{12} \\ \large S_{12} \ S_{22} \end{array} \right) \Large \left( \begin{array}{GC+23} a_1 \\ a_2 \end{array} \right) \ = \ \lambda \Large \left( \begin{array}{GC+23} a_1 \\ a_2 \end{array} \right)$

この式を見ると、左辺の最初の行列は VAR(A) そのものなので、 $VAR(A) \small %*% \Large a \ = \ \lambda \small %*% \Large a$ 　と書き換えることができる。

結局、

$VAR(A) \small %*% \Large a \ = \ \lambda \small %*% \Large a$ 　

という式を解くことで a を決めればよいことになる。 (%*% は行列、ベクトルのかけ算）

この形の式は様々な分野に出てきて、この文章でも最初の方に出てきた。同じことを写すと、

 行列 A があるとする。固有ベクトルを B、固有値を λ とすると、

$A \small %*% \Large B = \lambda \small %*% \Large B$ 　

   という関係が成り立つように、それらの値が決められる (%*% は行列、ベクトルのかけ
   算）。手で計算するのはとても大変だが、R 言語などを使うと簡単に計算できる。 左辺
   では行列 A だったところが、右辺では固有値 λ に置き換わって
   いる。 このことから、固有値は元の行列の性質を保持しながら縮約した、大切な部分を
   取り出したようなものと考えることもできる。

ここでも、

a は VAR(A) の固有ベクトル　（を並べた行列）

$\lambda$ 　は VAR(A) の固有値　（を並べたベクトル）

という形で解を求められる。

VAR(A) は、R 言語では cor(A) で計算できるが、上に書いたように SVD で求めた V, $\Sigma$ 　からも計算できる。

実際に、R 言語で上に書いた計算を試してみる。

 library(MASS)
 gyou <- 3
 retu <- 2
 mat1 <- matrix(0, gyou, retu)
 mat1[1, 1:2] <- c(1, 4)
 mat1[2, 1:2] <- c(2, 6)
 mat1[3, 1:2] <- c(9, 3)
 mat_scaled <- scale(mat1, scale = T)
 mat_svd <- svd(mat_scaled)

こうすると scale された行列は

$\left[ \begin{matrix} -0.6882 & -0.2182 \\ -0.4588 & 1.091 \\ 1.1471 & -0.8728 \end{matrix} \right ]$

それぞれの列の平均 = 0、不偏分散 = 1 にされている。

これを svd 関数で処理すると

 > mat_svd
 $d
 [1] 1.8307623 0.8051766
 
 $u
            [,1]       [,2]
 [1,] -0.1815424  0.7960584
 [2,] -0.5986357 -0.5552495
 [3,]  0.7801780 -0.2408089
 
 $v
            [,1]       [,2]
 [1,]  0.7071068 -0.7071068
 [2,] -0.7071068 -0.7071068

元の行列を再生するには　mat_svd$u %*% diag( mat_svd$d ) %*% t( mat_svd$V )　　で、

            [,1]       [,2]
 [1,] -0.6882472 -0.2182179
 [2,] -0.4588315  1.0910895
 [3,]  1.1470787 -0.8728716

となり、確かに元に戻っている。

$VAR(A) = \frac{1}{[N - 1]} A^T A$ 　を確かめると、

 > t(mat_scaled) %*% mat_scaled
           [,1]      [,2]
 [1,]  2.000000 -1.351691
 [2,] -1.351691  2.000000

で、これを　N - 1 = 2 で割るので、cor(mat_scaled) の値と一致する。

$VAR(A) = \frac{1}{[N - 1]} A^T A = \frac{1}{[N - 1]} V \Sigma^2 V^T$ 　を確かめると、

 > mat_svd$v %*% diag(mat_svd$d^2) %*% t(mat_svd$v)
           [,1]      [,2]
 [1,]  2.000000 -1.351691
 [2,] -1.351691  2.000000

で、これを　N - 1 = 2 で割るので合っている。

$VAR(A) \small %*% \Large a \ = \ \lambda \small %*% \Large a$ 　

で決められる、変換のためのベクトルを並べた行列 a は

 > eigen(cor(mat_scaled))$vectors
            [,1]       [,2]
 [1,] -0.7071068 -0.7071068
 [2,]  0.7071068 -0.7071068

$\lambda$ 　は

 > eigen(cor(mat_scaled))$values
 
 [1] 1.6758453 0.3241547

主成分分析では、主成分得点（スコア）、主成分負荷量、寄与率などの特有の値が出てくる。それぞれの意味と、上に書いたことの関係を考えてみる。

主成分得点（スコア）

この値は、「元の値を、変数をうまく混ぜ合わせることによって、データ同士ができるだけ分散するように変換した値」で、一番重要な値になる。上の例だと、元の値は scaleした値 mat_scaled で、それに変換のための（ベクトルを並べた)行列 a をかけ算することで計算できる。a は上に書いたように　eigen(cor(mat_scaled))$vectors　で計算できて、

 > mat_scaled %*% eigen(cor(mat_scaled))$vectors
            [,1]       [,2]
 [1,]  0.3323609  0.6409676
 [2,]  1.0959596 -0.4470739
 [3,] -1.4283205 -0.1938937

一列目は、最も分散が大きくなる変換によって計算された、データ同士を区別して特性を明らかにするために最も重要な値で、第一主成分得点（スコア）と呼ばれる。変換のための行列 a の一列目のベクトルで変換することで計算される。

今回の元の行列（データ）では、一行目、二行目は列2の方が値が大きい。3行目のデータは列1の方が値が大きい。そういう最も顕著な違いがあるので、第一主成分の値も3行目のデータがかけ離れた値になるように変換されている。

二列目は、変換のための行列 a の二列目で変換することで計算される。これは行1 > 行3 > 行2　　解釈すると、「列 1の値 + 列 2の値」の違いを表現しているのかもしれない。

群馬大学　青木先生が「R による統計処理」で公開されている PCA 関数での値と一致している。

 > pca(mat_scaled)$fs
           PC1        PC2
 #1  0.3323609  0.6409676
 #2  1.0959596 -0.4470739
 #3 -1.4283205 -0.1938937

主成分負荷量

これは、主成分分析で計算した主成分スコアと、元のデータの値との相関係数のことを言う。この値の絶対値が 1 に近ければ、元のデータの特徴を、その主成分がよく表現している・重要な主成分であるということになる。

 > cor(mat_scaled, mat_scaled %*% eigen(cor(mat_scaled))$vectors)
            [,1]       [,2]
 [1,] -0.9153812 -0.4025883
 [2,]  0.9153812 -0.4025883

群馬大学　青木先生が「R による統計処理」で公開されている PCA 関数での値と一致している。

 > pca(mat_scaled)$factor.loadings
           PC1        PC2
 X1 -0.9153812 -0.4025883
 X2  0.9153812 -0.4025883

寄与率

これは、変換のためのベクトル（を並べた行列）a を固有ベクトルとして計算した際に、セットとして出てくる固有値（のベクトル）の和を分母として、それぞれの固有値を分子として割合を計算し、100 をかけ算した値になる。

この値は主成分 1 から順番に出てきて、その値が大きいと、その主成分がデータの性質を強く表している（データの性質を表すのに寄与している）ことを示す。

 > eigen(cor(mat_scaled))$values / sum(eigen(cor(mat_scaled))$values) * 100

 [1] 83.79227 16.20773

群馬大学　青木先生が「R による統計処理」で公開されている PCA 関数での値と一致している。

 > pca(mat_scaled)$contribution
      PC1      PC2 
 83.79227 16.20773

 vim: set ts=8 sts=2 sw=2 et ft=a111_modified_flexwiki textwidth=0 lsp=12:

SVD -

目次

SVD 特異値分解

特異値分解と主成分分析 PCA

未定乗数法