Python中大型数组对Matlab的svds(A, k)的等价实现？

Question

我正在尝试把一些代码从Matlab转到Python，但遇到了一个问题。我找不到svds的对应功能。

我试过用numpy.corrcoef和numpy.linalg.eig，但numpy.corrcoef在处理大数组时（比如500 x 20000的数组）不太好用。

这是Matlab里的代码，可能会有帮助：

s = size(data, 2)
mean = sum(data, 2)/s
m_data = ( data - repmat(mean, 1, s) ) / sqrt(s - 1)
[res_u,res_s] = svds(m_data, s)
eigenvals = diag(res_s).^2
eigenvecs = res_u

Answer 1

如果你的矩阵是稀疏的，可以使用最近版本的 scipy.sparse.linalg.svds。这个工具是用来处理稀疏矩阵的奇异值分解（SVD），它是基于ARPACK的。不过根据我的经验，这个工具有点问题（不确定是scipy的问题还是ARPACK的问题），所以我建议你在使用时先做一些测试，看看得到的奇异值是否符合预期。

如果你的矩阵是密集的，500行乘以20000列虽然比较大，但在现在的普通电脑上处理起来还是可以的。你有使用 numpy.linalg.svd 并且设置了 full_matrices=False 吗？

你用 numpy.corrcoef 是为了什么呢？在你发的MATLAB代码片段中并没有类似的内容。

Answer 2

你想要的东西在 python 和 numpy 中大概是这样的（我稍微调整了一下，不是直接把 matlab 的代码翻译成 python 和 numpy，而是让它看起来更像 python 的风格，当然 matlab 的代码也可以做类似的调整）：

import numpy as np

def _cas(D):
    """Center at mean and standardize."""
    return (D- D.mean(1)[:, None])/ (D.shape[1]- 1)** .5

def example(D):
    """Eigenvalues and -vectors, based on SVD."""
    u, s, v= np.linalg.svd(D, full_matrices= False);
    return np.diag(s)** 2, u

if __name__ == '__main__':
    data= np.random.rand(5, 20)
    data= _cas(data) # preprocess data according to your requirements
    eigenvals, eigenvecs= example(data)
    print eigenvals
    print eigenvecs

但是你在这方面遇到了性能问题吗？

能不能具体说说你现在的性能情况，以及你希望提升多少？在我这台普通的电脑上，执行 example(.) 对一个随机生成的 (500, 20000) 矩阵大约需要 20 秒。

我可以通过基本的线性代数把执行时间减少到大约 2.5 秒（几乎提高了 10 倍）！如果你希望的性能比这个更好，那请详细说明一下你的 data 的“性质”！

你的数据来自哪里？你具体是用来计算特征值和特征向量的什么情况？也就是说，你的主要目标是什么？