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

Question

我正在尝试计算一个表达式，形式是 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⁶，S 大约 100x100，C 大约 100x10⁶。P 是对角矩阵，S 是对称矩阵，而 C 每行只有一个元素。这个方法将用于实现一个使用稀疏矩阵的卡尔曼滤波器，详情见

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

Answer 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 矩阵不大，你也可以使用密集矩阵的逆来解决。