在不反转矩阵的情况下求解 Ax=By答案

【问题标题】：Solving Ax=By without inverting matrices在不反转矩阵的情况下求解 Ax=By
【发布时间】：2019-12-15 07:36:16
【问题描述】：

我需要以 Ax=By 的形式求解 x 的方程。我知道我不应该通过反转 B 来解决它，但我无法用 scipy.gmres 或 linalg.solve 解决 B^-1Ax=y，因为当我尝试用 linalg.inv 反转 B 时它失败了。它返回错误消息“Singular matrix”。

还有其他方法可以反转矩阵吗？效率并不重要，因为我只需要做一次。我不想像 T=Ax 和 x 一样先求解方程两次。

【问题讨论】：

np.linalg.solve(A, B @ y)?
A,B= 2^16x2^16 x,y=2^16x1。给定矩阵 A、B 和 y。然后我必须将 B 和 y 相乘 2^12 次。这是一个循环。
np.linalg.solve 当 A 为 NxN 时为 O(N^3)，因此在您的情况下，这转化为至少 2^54 次操作，在 3 GHz 下至少需要 70 天，可能很多更长。如果 A 是稀疏的、对称的等，您可以寻找专用的求解器。
它很稀疏，这就是为什么我想将 scipy.gmre 用于 (B^-1*A)x=y。
或scipy.sparse.linalg.gmres(A, B.multiply(y))。

标签： python numpy matrix scipy linear-algebra

【解决方案1】：

我会加入用户hilberts-drinking-problem 的建议，而不是反转B，而是遵循两步方法：首先将B 和y 相乘以产生向量By，然后求解系统@987654326 @。下面用小数组作为测试数据来说明这种方法：

from numpy import array;
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import gmres

A = coo_matrix((3, 3), dtype=float)
A.setdiag( [ 2, 4, -1 ] )
A.setdiag( [ 2, -0.5 ], 1 )
print( "A", A )

B = coo_matrix((3, 2), dtype=float)
B.setdiag( [ 1, -2 ] )
B.setdiag( [ -0.5, 4 ], -1 )
print( "B", B )

y = array( [ 13, 5 ] )
print( "y", y )

By = B * y
print( "By", By )

x = gmres( A, By )
print( "x", x )

【讨论】：