使用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万行和列]。它也不是稀疏矩阵。大多数条目很小，但非零。当然，仅尝试构建该矩阵会引发内存错误。

我如何修改上面的代码以处理巨大的矩阵？

我像PyTables一样heard of solutions，但是我不确定如何应用它们，我也不知道它们是否适用于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

这应该为临时文件使用较少的内存，以为您可以通过使用out的np.matmul参数来做得更好。我应该添加一些内容来处理每个循环中的最后一个片段，以防dim无法被block_*整除，但希望您能理解。

对于不适合内存使用的数组，您可以通过上面注释中的链接来应用工具。