Rで見る主成分分析と特異値分解の関係

多次元主成分分析を行うにあたって、主成分分析と特異値分解の関係を復習してみた。

まず、データを作成する。

# x：1000個の正規分布N(0, 0.2)に従う乱数
# y : xを３倍し、そこにN(0, 0.3)に従うノイズを加えたもの

# x2, y2: x, yそれぞれについて平均を引いて標準偏差で割って標準化したもの

# x3, y3: x2, y2を大きさ１になるように規格化

# z2: cbind(x3, y3)で作成した1000 x 2の行列

x <- rnorm(1000, 0, 0.2)
y <- 3*x + rnorm(1000, 0, 0.3)
mx <- mean(x)
sdx <- sd(x)
my <- mean(y)
sdy <- sd(y)
x2 <- (x - mx)/sdx
y2 <- (y - my)/sdy

x2 %*% x2
[,1]
[1,] 999

y2 %*% y2
[,1]
[1,] 999

x3 <- x2/sqrt(999)
y3 <- y2/sqrt(999)

z2 <- cbind(x3, y3)

plot(x3, y3)

で確認できるように、第一主成分だけでほぼ説明できるデータである。

z2を特異値分解する。

# w2: z2のSVD
# w2$u 1000 x 2 主成分得点行列

w2 <- svd(z2)

w2の中身を見る。

str(w2)
List of 3
$ d: num [1:2] 1.376 0.327
$ u: num [1:1000, 1:2] 0.00732 0.01653 -0.00181 0.03669 -0.02696 ...
$ v: num [1:2, 1:2] 0.707 0.707 0.707 -0.707

# 主成分と元データの相関係数
cor(w2$u[,1], x3)
[1] 0.972893
cor(w2$u[,1], y3)
[1] 0.972893
cor(w2$u[,2], x3)
[1] 0.2312559
cor(w2$u[,2], y3)
[1] -0.2312559

# w2$uの直交性を利用して、データからw2$uを消して、因子負荷行列を作成
t(w2$u) %*% z2
x3 y3
[1,] 0.9728930 0.9728930
[2,] 0.2312559 -0.2312559

# 特異値の対角行列diad(w2$d)とt(w2$v)から、因子負荷行列を作成
diag(w2$d) %*% t(w2$v)
[,1] [,2]
[1,] 0.9728930 0.9728930
[2,] 0.2312559 -0.2312559

diag(w2$d) %*% t(w2$v) 2 x 2 因子負荷行列
因子負荷は主成分と元データの相関であることが確認できる。

主成分軸の表現
第一主成分のプロット
a <- seq(-0.2,0.2, 0.001)
b1 <- a * w2$v[1,1]
plot(x3, y3, xlim=c(-0.2, 0.2), ylim=c(-0.1, 0.15))
par(new=T)
plot(b1, b2, xlim=c(-0.2, 0.2), ylim=c(-0.1, 0.15), ty='l', col=rgb(1,0,0))b2 <- a * w2$v[2,1]

主成分分析では、w2＄vが相関係数行列の固有値問題を解いた場合の
固有ベクトルとして得られる。
vを主成分得点とした時の負荷行列を計算　
これが元データを回転させた主成分軸上の座標を与える
e <- w2$u %*% diag(w2$d)
max(e[,1])
[1] 0.174004
min(e[,1])
[1] -0.1284722
max(e[,2])
[1] 0.03881079
min(e[,2])
[1] -0.02979512
plot(e[,1],e[,2],xlim=c(-0.13,0.18), ylim=c(-0.13,0.18))

直交性を利用して、t(w2$v)を打ち消して、因子負荷行列を求める
上と同じプロットが得られる。
d <- z2 %*% w2$v
max(d[,1])
[1] 0.174004

max(d[,2])
[1] 0.03881079
min(d[,1])
[1] -0.1284722
min(d[,2])
[1] -0.02979512
plot(d[,1],d[,2],xlim=c(-0.13,0.18), ylim=c(-0.13,0.18))

因子負荷行列が主成分軸上での座標を表している。