计算大矩阵的均值和协方差（300000 x 70000）答案

【问题标题】：compute the mean and the covariance of a large matrix(300000 x 70000)计算大矩阵的均值和协方差（300000 x 70000）
【发布时间】：2018-07-16 23:46:32
【问题描述】：

我正在使用 Numpy 并尝试计算大型矩阵 (300000 x 70000) 的均值和协方差。我有 32GB 大小的可用内存。就计算效率和易于实施而言，这项任务的最佳实践是什么？

我目前的实现如下：

def compute_mean_variance(mat, chunk_size):
    row_count = mat.row_count
    col_count = mat.col_count
    # maintain the `x_sum`, `x2_sum` array
    # mean(x) = x_sum / row_count
    # var(x) = x2_sum / row_count - mean(x)**2
    x_sum = np.zeros([1, col_count])
    x2_sum = np.zeros([1, col_count])

    for i in range(0, row_count, chunk_size):
        sub_mat = mat[i:i+chunk_size, :]
        # in-memory sub_mat of size chunk_size x num_cols
        sub_mat = sub_mat.read().val
        x_sum += np.sum(sub_mat, 0)
        x2_sum += x2_sum + np.sum(sub_mat**2, 0)
    x_mean = x_sum / row_count
    x_var = x2_sum / row_count - x_mean ** 2
    return x_mean, x_var

有什么改进建议吗？

我发现下面的实现应该更容易理解。它还使用 numpy 计算列块的平均值和标准差。所以它应该更高效且数值稳定。

def compute_mean_std(mat, chunk_size):
    row_count = mat.row_count
    col_count = mat.col_count
    mean = np.zeros(col_count)
    std = np.zeros(col_count)

    for i in xrange(0, col_count, chunk_size):
        sub_mat = mat[:, i : i + chunk_size]
        # num_samples x chunk_size
        sub_mat = sub_mat.read().val
        mean[i : i + chunk_size] = np.mean(sub_mat, axis=0)
        std[i : i + chunk_size] = np.std(sub_mat, axis=0)

    return mean, std

【问题讨论】：

是python 2还是python 3？

标签： python numpy matrix linear-algebra

【解决方案1】：

我假设为了计算方差，您使用的是 Wiki 所称的 Naïve algorithm。但是，您可能会在那里找到：

因为x2_sum / row_count 和x_mean ** 2 可以是非常相似的数字，取消会导致结果的精度大大降低比用于浮点运算的固有精度执行计算。因此，该算法不应用于实践。如果标准差为相对于平均值来说很小。

作为替代方案，您可以使用two-pass algorithm，即首先计算均值，然后在计算方差时使用它。原则上，这似乎很浪费，因为必须对数据进行两次迭代。但是，用于计算方差的“均值”不必是真正的均值，一个合理的估计（可能仅从第一个块计算）就足够了。然后这将简化为assumed mean 的方法。

此外，一种可能性是将每个块的均值/方差的计算直接委托给 numpy，然后将它们组合起来，以便使用 parallel algorithm 获得整体均值/方差。

【讨论】：

【解决方案2】：

所以我在大学有一个项目，涉及不同矩阵乘法算法的时间复杂度测试。

我已经上传了源代码here。

我发现的优化之一是，您可以通过更改 for 循环的结构来优化数组访问，以一次只关注行而不是遍历列。这是由于缓存在空间局部性上的行为方式（即您的计算机尝试针对二维数组中并排而不是逐行的数组元素进行优化）

此外，如果这些是“稀疏”矩阵（大量归零元素），您可以更改数据结构以仅记录非归零元素。

显然，如果给定一个普通矩阵，将它们转换为稀疏矩阵的计算可能不值得，但我只是认为这些是值得分享的观察结果 :)

【讨论】：