具有多个输入的Scipy线性操作符

Question

我需要反转一个大而密集的矩阵，我希望使用Scipy的gmres来完成。幸运的是，密集矩阵A遵循一种模式，我不需要在内存中存储矩阵。 LinearOperator类允许我们构建一个对象，它可以作为GMRES的矩阵，并直接计算矩阵向量积A*v。也就是说，我们编写一个函数mv(v)，它以向量v作为输入，并返回mv(v) = A*v。然后，我们可以使用LinearOperator类创建A_LinOp = LinearOperator(shape = shape, matvec = mv)。我们可以将线性算子放入Scipy的gmres命令中，以评估矩阵向量乘积，而无需完全加载A到内存中。 LinearOperator的文档在这里：LinearOperator Documentation。

这是我的问题：为了编写计算矩阵向量积mv(v) = A*v的例程，我需要另一个输入向量C。 A中的条目形式为A[i,j] = f(C[i] - C[j])。因此，我真正想要的是，mv有两个输入，一个是固定向量输入C，另一个是可变输入v，我们希望计算A*v。

MATLAB有类似的设置，其中我们编写x = gmres(@(v) mv(v,C),b)，其中b是问题Ax = b的右手边，mv是以可变输入v作为变量的函数，我们希望计算A*v，而C是已知的固定向量，我们需要用它来组装A。

我的问题是，我无法弄清楚如何让LinearOperator类接受两个输入，一个可变和一个“固定”，就像在MATLAB中一样。是否有一种方法在SciPy中执行类似的操作？或者，如果有人知道更好的方法来反转一个大而密集的矩阵(50000, 50000)，其中条目遵循一种模式，我将非常感激任何建议。

谢谢！

编辑：我应该明确这个信息。该矩阵实际上是（分块形式）[A C; C^T 0]，其中A是N x N（N很大），C是N x 3，0是3 x 3，而C^T是C的转置。这个数组C与上面提到的那个数组是相同的。A的条目遵循模式A[i,j] = f(C[i]-C[j])。

我写了mv(v,C)来逐行构造A*v[i]，其中i=0,N，通过计算求和f(C[i]-C[j)*v[j]（实际上我做的是numpy.dot(FC,v)，其中FC[j] = f(C[i]-C[j])，这个方法效果不错）。然后，在结束时完成对C^T行的计算。我希望最终用multiprocessing调用来并行化for循环，但这是考虑未来的事情。我还将研究使用Cython加速计算。

- user35959

---

你不一定需要让mv接受两个参数，你可以使用一个闭包（甚至是functools.partial）替代它：make_mv = lambda c: lambda v: sum(f(c[i] - c[j])...(但为了可读性最好使用def)，而且这个函数可能应该用Cython编写，否则对于这么大的矩阵来说速度会太慢。 - jorgeca

谢谢！我认为lambda结构正是我想要的。我可以在我的脚本中定义mv2 = lambda v: mv(v,C)，然后将mv2(v)输入gmres。由于我要处理的大矩阵，我的方法速度相当慢，所以我会研究Cython。 - user35959

要小心，因为那不是一个闭包，而且在某些时候会让你吃亏。你的 mv2 函数将使用当前作用域中的任何 C，而不是创建函数时所在作用域中的 C。如果你想让 mv2 引用特定的 C，你绝对需要上述两种方法之一，然后使用 mv = make_mv(C)。 - jorgeca

Answer 1

虽然有些晚了，但如果您还感兴趣...

由于A矩阵是一个非线性变换的秩-2矩阵，因此它的秩非常低。而且它是对称的。这意味着它很容易求逆：使用截断特征值分解，例如5个特征值：A = U*S*U'，然后求逆：A^-1 = U*S^-1*U'。S是对角线矩阵，因此这很便宜。您可以使用eigh获得截断特征值分解。

这样就处理了A。然后对于其余部分：使用block matrix inversion formula。看起来很复杂，但我敢打赌，它比您之前使用的直接方法快50倍，可以赢得100,000,000普鲁士法郎的赌注。

- Patrick Mineault

---

谢谢您的回复。我一定会研究这个，因为任何比直接方法更好的改进都将受到极大的赞赏。 - user35959

Answer 2

我曾经遇到过与你相同的情况（比你晚几年），尝试使用多个参数传递给LinearOperator，但是为了解决另一个问题。我找到的解决方案是使用全局变量，以避免将变量作为参数传递给函数。

- Fernando