Help us understand the problem. What is going on with this article?

特異値分解と固有値分解が同じ結果を出していることの確認(Rのsvd、eigen関数)

More than 3 years have passed since last update.

ある行列A(サンプルN×変数 D)があったとして、
$A = U \Sigma V^T$
のように行列分解するのが、特異値分解(SVD)。

Aの共分散行列
$S = \frac{1}{N-1}A^T A$
に対して、
$SV = V\Lambda$
となるような、固有ベクトル$V$を探す、すなわち、
$S = V\Lambda V^T$
のように行列分解するのが、固有値分解(EVD、Vは特異値分解のものと同じ)。

$S = \frac{1}{N-1}A^T A = \frac{1}{N-1}(U \Sigma V^T)^T U \Sigma V^T = \frac{1}{N-1}V \Sigma U^T U \Sigma V^T = \frac{1}{N-1}V \Sigma^2 V^T$
となることから、

SVDの特異値とEVDの固有値には以下のような関係があるはず。
$\frac{1}{N-1}\Sigma^2 = \Lambda$

というのは、数式上ではわかるが、実際にRの関数が同じ値を返すか不安だったので、確認してみた。

準備

まずデータを用意する。

# 10サンプル × 15次元のデータとする
set.seed(123)
A <- matrix(runif(150), nrow=10)
A <- scale(A, scale=FALSE) # あらかじめmean-centeringはしておく

続いて、SVD、EVDを実行

# SVD
res.svd <- svd(A)
# EVD
S <- t(A) %*% A / (10 - 1)
res.evd <- eigen(S)

固有値は同じか

まずは、SVDが返す特異値の自乗$\Sigma^2$と、EVDが返す固有値$\Lambda$が同じかを調べた

RではSVDはsvd関数、EVDはeigen関数があるので、これらを使った

# SVDの特異値^2 = EVDの固有値
plot(1:10, res.svd$d^2/(10 - 1), col="black", type="l", ylim=c(0,4), xlim=c(1,10), main="Eigen Value")
par(new=TRUE)
plot(1:10, res.evd$values[1:10], col="blue", type="l", ylim=c(0,4), xlim=c(1,10), ann=FALSE)
legend("topright", legend=c("SVD", "EVD"), col=c("black", "blue"), pch=16)

Rplot.png

かなり一致しているのが見える

固有ベクトルは同じか

次に、固有ベクトル$V$がsvdとeigen関数で同じかを調べた

固有ベクトルは長さ1のベクトルなので、同じ方向を向いているほど、内積が1に近づく

180度反対の向きのものも正解なので、内積の絶対値を確認した

evec.svd <- res.svd$v
evec.evd <- res.evd$vectors[,1:10]
image(abs(t(evec.svd) %*% evec.evd), xlab="SVD", ylab="EVD", main="Eigen Vector")

Rplot.png

最後が怪しいが、それ以外は一致した

ローディングは同じか

最後に、長さの違うもう一つの固有ベクトル$U$が同じになるかを調べた

SVDでは、
$A = U \Sigma V^T$
の$U$ですでに計算してある。

EVDでは、データ行列$A$と固有ベクトル$V$と、固有値$\Lambda$しかないが、$A$に$V$をかければ、$U$がもとまる

$A V = U \Sigma$

ただし、$\Sigma$分の重さがかかってるので、長さを1にして$\Sigma$を取り除く

loading.evd <- A %*% evec.evd
loading.svd <- res.svd$u
apply(loading.evd, 2, function(x){norm(x, "2")}) # この固有ベクトルにはまだ、固有値分の重みがかかっている
apply(loading.svd, 2, function(x){norm(x, "2")}) # SVDのVはΣと分離されているからOk
for(i in 1:ncol(loading.evd)){
    loading.evd[, i] <- loading.evd[, i] / norm(loading.evd[, i], "2")
}
image(abs(t(loading.evd) %*% loading.svd), xlab="EVD", ylab="SVD", main="Loading")

Rplot.png

やはり、最後(10番目の固有ベクトル)が何かおかしい

固有ベクトルなら、互いに直行しているため、内積が0になる

そのため、SVDの$U$と、EVDの$U$で各々内積をみた

image(crossprod(loading.svd), main="SVD")
image(crossprod(loading.evd), main="EVD")

Rplot*-1.png
Rplot*-0.png

これを見る限り、EVDの10番目の軸が他の軸と直行していない

このように、原因はよくわからないが、eigen関数の高次の部分が不安定になる事例が見受けられた

antiplastics
書くぜーどんどん書くぜー(書いてない
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away