R对稀疏矩阵的内部处理

Question

我一直在比较Python和R中几个PCA实现的性能，并注意到一个有趣的行为：
虽然似乎无法在Python中计算稀疏矩阵的PCA（唯一的方法是scikit-learn的TruncatedSVD，但它不支持所需的均值居中以等效于PCA的协方差解决方案。 他们的论点是，这将破坏矩阵的稀疏属性。 Facebook的PCA算法或scikit learn中的PCA/randomPCA方法等其他实现由于类似的原因不支持稀疏矩阵。

虽然我都理解，但是几个R包（如irlba，rsvd等）能够处理稀疏矩阵（例如使用rsparsematrix生成），并且甚至允许特定的center=True参数。

我的问题是，R是如何处理内部的，因为它似乎比可比的Python实现更有效率。 R是否仍通过进行绝对缩放来保持稀疏性（这理论上会使结果失真，但至少保持稀疏性）？ 或者是否有任何方式可以显式地存储零值的平均值，并仅存储一次（而不是为每个值单独存储）？

为了解除保持状态： R如何在不爆炸式使用RAM的情况下内部存储具有均值居中的矩阵。 希望这足够简洁....

- dennlinger

---

这是一个有趣的问题，但我不确定SO是否是询问它的最佳场所。您可以考虑在交叉验证上提问，那里更有可能得到答案。 - piman314

谢谢你的提示。我正在考虑使用SO，因为它可能会被标记为离题在Cross Validated中。如果问题仍未得到解答，也许我会在那里发问。 - dennlinger

我认为答案可以在 ?irlba 中找到：“使用可选的‘center’参数，隐式从‘A’的每一列中减去‘center’向量的值，计算'sweep(A, 2, center, FUN=-)'的截断SVD，而不需要显式形成居中矩阵”（强调添加；换句话说，这是一种算法技巧而不是存储技巧）。然后你需要查看代码：https://github.com/bwlewis/irlba/blob/master/R/irlba.R，以了解`center`参数如何在算法中实际使用。 - Ben Bolker

也许你可以看一下这个链接。 - rpanai

谢谢提供链接，但我不太确定这个链接如何帮助我？文章中甚至没有提到稀疏矩阵，而且代码完全基于Python... 我已经知道Python不支持稀疏处理（至少不是scikit-learn的“高效”包）。 - dennlinger

Answer 1

重点在于部分奇异值分解的底层实现（重新启动的Lanczos双对角化C代码）不会存储矩阵。相反，您记录从前一次迭代获得的一小组向量应用于矩阵的线性操作的结果。

与其解释C代码中使用的具体方法，这是相当高级的（请参见论文进行描述），我将用一个更简单的算法来解释它，该算法捕捉了如何保持稀疏性效率的关键思想：幂法（或子空间迭代法，以便将其推广到多个特征值）。该算法通过迭代应用线性算子，然后归一化（或正交化一小组向量，在子空间迭代的情况下）来返回矩阵A的最大特征值。

每次迭代时所做的是：

v=A*v
v=v/norm(v)

矩阵乘法步骤是关键的一步，因此让我们看看在使用居中A时会发生什么。 居中A的矩阵公式（其中 center 为具有平均列值的向量， ones 为全1向量）为：

A_center=A-ones*transpose(center)

因此，如果我们将迭代算法应用于这个新矩阵，我们将得到

v=A*v-dotproduct(center,v)*ones

由于A是稀疏矩阵，我们可以在(A,v)上使用稀疏矩阵-向量乘积，而-dotproduct(center,v)*ones仅涉及从结果向量中减去中心点和v的点积，该向量在A的维度上是线性的。

- Juan Carlos Ramirez