奇异值分解（SVD）：使用R中的例子

原文

如果你有深入研究过多元变量分析的统计计算，那么，你肯定也遇到过奇异值分解（SVD）。 SVD是分解数据的一种主要技术，类似于对因分析和主成分分析。在这一帖子中，作者以一种直观的视角来解释SVD到底是怎样的？通过使用R中的例子进行验证。如果你之前还未曾学习过SVD，那就略过这帖子吧！这只是写给听说过SVD，想要更好地理解SVD，但又讨厌数学（公式、符号）的少数人。

下表是一张关于标准化残差的列联表，它列出一份报纸的读者受教育水平与阅读本报的深入程度两者之间的关系。

生成这张表的R代码如下：

education.by.readership = matrix(c(5, 18, 19, 12, 3, 7, 46, 29, 40, 7, 2, 20, 39, 49, 16),     nrow = 5,    dimnames = list(        "Level of education" = c("Some primary", "Primary completed", "Some secondary", "Secondary completed", "Some tertiary"),        "Category of readership" = c("Glance", "Fairly thorough", "Very thorough"))) O = education.by.readership / sum(education.by.readership)E = rowSums(O) %o% colSums(O)Z = (O - E) / sqrt(E)

---

如何计算SVD

上表是一个数值矩阵。作者将其称为Z。通过使用svd函数来实现奇异值分解。下面的代码计算了矩阵Z的奇异值分解，并将其赋值给了称为SVD的一个新对象，它包含有向量d，以及两个矩阵u和v。向量d包含奇异值。第一个矩阵u包含左奇异向量，第二个矩阵v包含右奇异向量。左奇异向量对应矩阵Z的“行”，右奇异向量对应矩阵Z的“列”。

SVD = svd(Z)

$d
[1] 2.652708e-01 1.135421e-01 2.718254e-17

$u
[,1] [,2] [,3]
[1,] -0.4386666 -0.42375592 -0.3714480
[2,] -0.6516462 0.35501142 -0.4906752
[3,] 0.1603076 -0.67246939 -0.2423522
[4,] 0.3711005 0.48847409 -0.3785281
[5,] 0.4685240 -0.05979793 -0.6474922

$v
[,1] [,2] [,3]
[1,] -0.4097795 -0.80584644 -0.4274252
[2,] -0.4887795 0.58960413 -0.6430097
[3,] 0.7701788 -0.05457549 -0.6354889

---

复原矩阵Z

奇异值分解（SVD）有4个有用的特性。第一个特性，上文的两个矩阵和向量可以进行乘法运算从而重新创建原始矩阵Z。矩阵Z中有一个值，-0.064751，它位于第五行第二列。从SVD的计算结果中，我们可以将其计算得到，这可以通过将d的每一元素与矩阵u的第五行和矩阵v的第二行的每一元素进行相乘，进而累加得到。也就是说：

-0.064751 = 0.2652708*0.468524*(-0.4887795) + 0.1135421*(-0.0597979)0.5896041 + 0(-0.6474922)*(-0.6430097)

在R中可以这样实现：

sum(SVD$d * SVD$u[5, ] * SVD$v[2, ])

更佳的是，如果我们想要一次性重新计算矩阵Z中的所有数值，可以借助矩阵代数的一点知识，在R中的实现如下：

SVD$u %*% diag(SVD$d) %*% t(SVD$v)

现在乍一看，这一特性似乎看起来没什么用。实际上，这看起来似乎甚至不是一种聪明的做法。因为，一开始，我们的一个矩阵中只有15个数值。现在，我们有一个向量和两个矩阵，共包含有27个数值。这似乎是在“开倒车”。

---

数据降维

SVD的第二个有用特性与向量d中的数值有关，它们以降序排列（可能与ties参数有关，类似地，参考rank函数）。为什么这很重要？让我们来看看向量d中的最后一个值，2.71825390754254E-17。实际上，这可以看作是0（因为计算机“很难”精确地计算出0）。当我们考虑复原矩阵Z时，我们可以忽略向量d中的最后一值，进而也可以忽略矩阵u和v的最后一列。矩阵u和v最后一列的值被0相乘，因而是不相关的。现在，我们可以看到剩下18个数值了（原本27个）。这仍然比最初的15个数值还要多。

其实，在描述原始矩阵Z方面，向量d中的值告诉了我们矩阵u和v的列之间的相对重要性。我们可以计算矩阵Z中的方差，首先对向量d中的值求平方，然后将其表达成比例的形式。可以运行下面R代码来实现，它表明第一维可以解释矩阵Z中85%的变异（方差）。

variance.explained = prop.table(svd(Z)$d^2)

所以，如果可以接受忽略掉原始矩阵Z中15%的变异信息，那我们只需要查看矩阵u和v中的第一列就可以了。现在，可以看到，只剩下不多于原始数值个数的一半了。

看到数值个数减少一半似乎不太具有充分的说服力。然而，如果数据集越大，那么利用SVD将节省（存储）的越多。例如，如果有一个20行20列的矩阵，利用SVD，或许我们只需要查看矩阵u和v的第一列，从而只需要考虑原始矩阵中数值个数的10%就可以了。这是与主成分分析和对因分析相类似分解技术的基本逻辑。除了减少需要查看的数值个数之外，这也能减少用图来可视化这些数值的所花费的时间。基于一个矩阵的20个列来画图真的不是一个好的方法，但是，只考虑2列来画图就直观得多了。

第三和第四个特性
SVD的第三个特性是矩阵u的行对应原始矩阵Z的行变量的各类别（这里是”水平”），矩阵v的行对应原始矩阵Z的列变量的各类别（这里是”水平”）。第四个特性是矩阵u的列之间相互正交，矩阵v的列之间相互正交。

---