Python左乘矩阵与稀疏矩阵的逆答案

【问题标题】：Python left multiplication of of a matrix with inverse of a sparse matrixPython左乘矩阵与稀疏矩阵的逆
【发布时间】：2011-10-14 12:59:27
【问题描述】：

我正在尝试计算 K = P*C.T*S^-1 形式的表达式（卡尔曼滤波器的实现）

所有涉及的矩阵都是稀疏的，我当然想避免计算实际的逆矩阵。

我尝试过使用

import scipy.sparse.linalg as spln

self.K = self.P.dot(spln.spsolve(S.T, C).T)

问题在于 spsolve 期望它的第二个参数是向量而不是矩阵。

编辑：澄清一下，这个问题可以在 Matlab 中通过 K = P * (C / S) 来解决，所以我正在寻找一种类似于 spsolve 但可以接受矩阵作为其第二个参数的方法。这当然可以通过将 C 拆分为多个列向量 c1..cn 并解决每个列向量的问题，然后将它们重新组合成一个矩阵来完成，但我怀疑这样做既麻烦又低效。

编辑2&3：矩阵的尺寸通常约为 P~10⁶x10^6，S~100x100，C=100x10⁶。 P对角线和S对称，C每行只有一个元素。它将用于使用稀疏矩阵实现卡尔曼滤波器，请参阅

http://en.wikipedia.org/wiki/Kalman_filter#The_Kalman_filter

【问题讨论】：

如果不计算逆，就无法计算K。您可以在不计算逆的情况下为某个向量x 计算Kx，这将涉及求解线性系统。
我不同意，在 matlab 中，我的问题的解决方案很简单： K = P*(S' \ C)' 或等价的 K = P*(C / S)矩阵而不是向量不会改变推理，您可以按照您所说的对 C 中的每一列求解一次。我的问题是 spsolve 将我限制为 C 作为向量，而在 Matlab 中它可以也可以用于矩阵。根据矩阵的维度，这仍然比计算实际的逆矩阵要有效得多。
对不起，我隐含地假设所有矩阵都是正方形的。那么为什么不简单地遍历 C 的列并求解每一列呢？由于您每次都必须求解线性系统，因此循环开销可以忽略不计。
刚刚看到您的编辑。不，我认为迭代 C 的 colmuns 不会是低效的。不过，根据线性系统的大小，请考虑使用迭代求解器而不是 spsolve()。
如果没有更好的方法，我肯定会尝试通过遍历列来解决它，但我有一种预感，它不会那么有效。矩阵的尺寸通常约为 P~10⁶x10^6，S~100x100，C=100x10⁶。 P 和 S 将是对角线，而 C 每行只有一个元素。我也会用这些信息更新我的问题。

标签： python numpy scipy sparse-matrix linear-equation

【解决方案1】：

作为一种解决方法可以做

import numpy as np
from scipy.sparse.linalg import splu

def spsolve2(a, b):
    a_lu = splu(a)
    out = np.empty((A.shape[1], b.shape[1]))
    for j in xrange(b.shape[1]):
        out[:,j] = a_lu.solve(b[:,j])
    return out

self.K = self.P.dot(spsolve2(S.T, C).T)

但是是的，spsolve 不接受矩阵是一个错误。

但是，由于您的 S 不是很大，您也可以使用密集逆。

【讨论】：