SVD初探（1）

时间：2015-03-22 22:12:00 收藏：0 阅读：479

谈到对高维度数据降维，SVD和PCA是最基本最常用的降维工具，在探索性数据分析中也有诸多妙用。SVD常常和PCA一同提起，至于他们之间的区别，哈佛在edX上的线上课PH525.3x Advanced Statistics for the Life Sciences有提到仅仅是去均值化的区别（PCA本身包含normalization过程），不过我也查到SVD比PCA更为稳定一说（比如处理Läuchli matrix，http://math.stackexchange.com/questions/3869/what-is-the-intuitive-relationship-between-svd-and-pca）。我本身对线性代数不是太熟悉，就经验来说，我曾遇见过在column远远大于row的时候PCA无法进行而SVD可行的情况。

13年上Coursera的data analysis的时候（JHU课程），论坛上对SVD究竟对数据是个怎样的过程进行过热烈地讨论，有一个跟帖的解析非常棒（非数学解析更直观），我当时简单地做了笔记。时至今日已经查不到原来的帖子了，故结合该贴精要部分，重新对SVD的性质作进一步探索。

假设对一个矩阵M（n*p）进行SVD，我们得到U，D，V三个矩阵，有M = UDV‘：

其中U的维度是n*p，D的维度是p*p，V的维度是p*p。U和V是正交矩阵，D是对角矩阵。D代表了每一个主成分包含的方差。

这里大部分数据分析书本会告诉你V是eigenvector，D是eigenvalue，DV‘代表变化后的新坐标轴，而U是M在新坐标轴的坐标，我会在后面解析。

我们用R进行SVD操作，为了便于说明使用二位数据，并假设变量（特征）间高度相关：

rho <- 0.95
nn <- 1000
set.seed(1)
xx <- rnorm(nn)
yy <- rho*xx + sqrt(1-rho^2)*rnorm(nn)

为了让数据看上去更有特征性，我们对前半部分数据和后半部分数据添加颜色。

ord <- order(xx)
xx <- xx[ord]
yy <- yy[ord]
plot(xx, yy, pch=19, xlab="", ylab="", col = rep(3:4, each = nn/2))
abline(v=0, col="gray")
abline(h=0, col="gray")

技术分享

我们可以目测由于x和y两个变量高度相关，从x可以得到很多y的信息，反之亦然。

我们运行SVD计算：

SVD <- svd(cbind(xx, yy))
SVD$d^2/sum(SVD$d^2)
#[1] 0.97487248 0.02512752

得到的D是一个长度p的向量，可以用diag命令转换进行定义下的计算：cbind(xx, yy) = SVD$u %*% diag(SVD$d) %*% t(SVD%v)

如同预计，第一个主成分（PC）解释了约97.5%的方差。

那么V代表什么意思呢，我们可以画出来：

arrows(x0=0, y0=0, x1=SVD$v[1,1], y1=SVD$v[2,1], lwd=3, col="red")
arrows(x0=0, y0=0, x1=SVD$v[1,2], y1=SVD$v[2,2], lwd=3, col="red")

技术分享

连接原点和V的第一个column的坐标代表了右上方向，长度为1（模）的箭头，沿着箭头指向数据方差最大的方向。同理，V的第二个column指向左上，指向剩余方差里方差第二大的方向，并与第一个箭头垂直（正交）。注意这里的箭头是只有方向的单位向量，至于坐标轴怎么伸缩的信息包含在D。

我们接着看看U的第一个column和第二个column代表了什么。

plot(SVD$u[, 1], SVD$u[, 2], pch=19, xlab="", ylab="", col = rep(3:4, each = nn/2))

技术分享

我们可以看到整个数据集沿着上图的箭头发生了旋转，但好像数据尺度不对？

zz <- diag(SVD$d) %*% t(SVD$u)
plot(zz[1,], zz[2,], pch=19, xlab="", ylab="", col = rep(3:4, each = nn/2))

技术分享

在基因分析中，往往row代表了基因（变量、特征），而column是不同样本，那么这个时候SVD又怎么处理呢？

同样用刚才的数据构建矩阵，只不过x和y变量变成了row，进行SVD后沿箭头方向重建坐标轴：

SVD_r <- svd(rbind(xx, yy))
zz_r <- diag(SVD_r$d) %*% t(SVD_r$v)
plot(-zz_r[1, ], zz_r[2, ], pch=19,xlab="", ylab="", col = rep(3:4, each = nn/2))

技术分享

注意由于SVD对正负号不敏感（即可有多个解），我对第一个PC添加了符号使其看起来更像一个完整的旋转过程。

原文：http://www.cnblogs.com/nutastray/p/4357927.html