numpy.linalg.solve 具有多个相关的右侧答案

【问题标题】：numpy.linalg.solve with multiple dependent right-hand sidenumpy.linalg.solve 具有多个相关的右侧
【发布时间】：2019-12-12 05:54:49
【问题描述】：

我正在尝试求解像 Ax=b 这样的方程组，它具有满秩矩阵 A 和右侧 b，包括两个具有相同左侧的不同系统的两列。（Ax=b1，Ax=b2，b=连接（（b1，b2），轴=1））。关键是我们对两个系统有相同的 A。因此，我们应该能够将第一个系统的信息用于第二个系统，即 A 的倒数。

numpy.linalg.solve 如果 b 的列是独立的，那么很容易解决这个系统，这不是我的情况。在我的例子中，b 的第二列取决于第一个系统的解。

我尝试分解矩阵 A 并使用此分解来求解两个系统。但是，它根本没有效率。考虑到 A 不是对称的，我使用了 RQ 和 LU 分解。

from scipy.linalg import lu
import numpy as np
from scipy.linalg import solve_triangular
def qr_solver(a, b):
   q, r = np.linalg.qr(a)
   z = solve_triangular(r, b, lower=False)
   ans = np.dot(q.T, z)
   return ans
def plu_solver(a, b):
   per, l, u = lu(kkt_matrix)
   z = np.dot(per.T , b)
   x = solve_triangular(l, z, lower=True)
   ans = solve_triangular(u, x)
   return ans

【问题讨论】：

这个问题有什么问题吗？我们有可以在这里使用的 LAPACK 模块吗？
你有什么问题？你是知道如何修改b2的人。你不能只做一个solve，创建新的b 并重复吗？
考虑到解决两个系统的时间，总时间将增加两倍！当你有相同的左侧时，对第二个方程也使用第一个分解更合理。
这里的重点是，对于两个系统，我们有相同的 A 和不同的 b。如果左侧是独立的，Numpy 可以有效地解决这些系统。但就我而言，第二个方程的“b”取决于第一个系统的答案。现在清楚了吗？
为了更清楚地查看这个等式：stackoverflow.com/questions/48387261/…

标签： python-3.x numpy scipy linear-algebra

【解决方案1】：

scipy 针对此类问题公开了lu_factor 和lu_solve：

from scipy.linalg import lu_factor, lu_solve

# Solving Ax = b1, Ay = f(x) with same A

lu, pivot = lu_factor(A)
x = lu_solve((lu, pivot), b1)
b2 = f(x)
y = lu_solve((lu, pivot), b2)

因此，如果 RHS 向量不是线性独立的（隐式 Runge-Kutta 方案就是一个很好的例子），您可以对 LHS 进行一次因式分解，然后根据需要重复使用它来求解。

【讨论】：

我只是尝试实现和比较你的漂亮代码。它似乎工作得很好。

【解决方案2】：

我只是尝试编写“Talonmies”所说的内容：

    from scipy.linalg import lu_factor, lu_solve
    from scipy.linalg import cho_factor
    from scipy.linalg import solve
    import numpy as np
    import time
    # Solving Ax = b1, Ay = f(x) with same A
    A = np.random.normal(1, 10, (2000, 2000))
    b1 = np.random.normal(1, 10, (2000, 1))
    b2 = np.random.normal(1, 10, (2000, 1))
    start = time.time()
    lu, pivot = lu_factor(A)
    x = lu_solve((lu, pivot), b1)
    y = lu_solve((lu, pivot), b2)
    end = time.time()
    time_factorization = start-end
    start = time.time()
    x = np.linalg.solve(A, b1)
    y = np.linalg.solve(A, b1)
   end = time.time()
   time_solve = start-end
   print('time_factorization:', time_factorization, 'time_solve:', time_solve)

结果如下： time_factorization: 0.6500372886657715 time_solve: 0.9420537948608398

看来这种分解（得到 lu 和 pivots）确实很有效。

【讨论】：