在 numpy 中处理非常大的矩阵

Question

我有一个转换矩阵，我想为其计算一个稳态向量。我使用的代码改编自this question，它适用于正常大小的矩阵：

def steady_state(matrix):
    dim = matrix.shape[0]
    q = (matrix - np.eye(dim))
    ones = np.ones(dim)
    q = np.c_[q, ones]
    qtq = np.dot(q, q.T)
    bqt = np.ones(dim)
    return np.linalg.solve(qtq, bqt)

但是，我正在使用的矩阵大约有150 万行和列。它也不是稀疏矩阵；大多数条目很小但非零。当然，只是尝试构建该矩阵会引发内存错误。

我怎样才能修改上面的代码来处理巨大的矩阵？我有 heard of solutions 喜欢 PyTables，但我不知道如何应用它们，我不知道是否他们会为np.linalg.solve 之类的任务工作。

对于 numpy 非常陌生并且对线性代数非常缺乏经验，我非常感谢我的案例中的一个示例。我愿意使用 numpy 以外的东西，如果需要的话，甚至可以使用 Python 以外的东西。

Answer 1

以下是一些想法：

我们可以利用这样一个事实，即任何初始概率向量在时间演化下都会收敛于稳态（假设它是遍历的、非周期性的、规则的等）。

对于我们可以使用的小矩阵

def steady_state(matrix):
    dim = matrix.shape[0]
    prob = np.ones(dim) / dim
    other = np.zeros(dim)
    while np.linalg.norm(prob - other) > 1e-3:
        other = prob.copy()
        prob = other @ matrix
    return prob

（我认为问题中函数假定的约定是分布按行排列）。

现在我们可以利用矩阵乘法和norm 可以逐块完成这一事实：

def steady_state_chunk(matrix, block_in=100, block_out=10):
    dim = matrix.shape[0]
    prob = np.ones(dim) / dim
    error = 1.
    while error > 1e-3:
        error = 0.
        other = prob.copy()
        for i in range(0, dim, block_out):
            outs = np.s_[i:i+block_out]
            vec_out = np.zeros(block_out)
            for j in range(0, dim, block_in):
                ins = np.s_[j:j+block_in]
                vec_out += other[ins] @ matrix[ins, outs]
            error += np.linalg.norm(vec_out - prob[outs])**2
            prob[outs] = vec_out
        error = np.sqrt(error)
    return prob

这应该使用更少的临时内存，认为使用np.matmul 的out 参数可以做得更好。 我应该添加一些东西来处理每个循环中的最后一个切片，以防dim 不能被block_* 整除，但我希望你明白这一点。

对于不适合内存的数组，您可以应用上面 cmets 中链接中的工具。