在 Python 中构建并快速求解线性方程组的最佳方法是什么？

Question

假设我有一个大型线性方程组：

X_1*a_11 + X_2*a_21 + ... + X_n*a_n1 = b_1
.
.
.
X_1*a_1m + X_2*a_2m + ... + X_n*anm = b_m

尺寸n 和m 是固定的（不是符号）。 我可以找到并插入每个系数a_ij 的值，我想解决所有X_i。如果我可以将这些写成矩阵形式，例如，我可以使用numpy.linalg.solve。但是，我不想用系数手动构造这个矩阵。有没有一种有效的方法可以将此方程组转换为矩阵并使用numpy.linalg.solve 求解？

我知道我可以将所有 X 定义为 Sympy 符号，创建方程式列表并使用例如 Sympy 的 linear_eq_to_matrix 将符号方程式列表转换为符号矩阵。我也可以使用 Sympy 的 solve 或 linsolve 使用这个方程列表和一个未知数列表，但我认为与使用 Numpy 相比存在性能损失。

解决这个问题的最佳方法是什么？

编辑

一个非常简单的例子

from sympy import *
import numpy as np

dim = 4
# Getting coefficients
a = np.random.uniform(low = -3, high = 3, size= (dim,dim)) 
b = np.random.uniform(low = -3, high = 3, size= (dim,))
# Numerical array for solutions
n_X = np.empty(shape = (2,))
n_Y = np.empty(shape = (2,))

# Defining unknowns
s_X_1 = Symbol('X_1')
s_Y_1 = Symbol('Y_1')
s_X_2 = Symbol('X_2')
s_Y_2 = Symbol('Y_2')

#Defining equations
eq1 = s_X_1*a[0,0] + s_Y_1*a[1,0] + s_X_2*a[2,0] + s_Y_2*a[3,0]- b[0]
eq2 = s_X_1*a[0,1] + s_Y_1*a[1,1] + s_X_2*a[2,1] + s_Y_2*a[3,1]- b[1]
eq3 = s_X_1*a[0,2] + s_Y_1*a[1,2] + s_X_2*a[2,2] + s_Y_2*a[3,2]- b[2]
eq4 = s_X_1*a[0,3] + s_Y_1*a[1,3] + s_X_2*a[2,3] + s_Y_2*a[3,3]- b[3]

unknowns = [s_X_1, s_X_2] + [s_Y_1, s_Y_2]
eqs = [eq1, eq2, eq3, eq4]
# Transforming the system into matrix form and solving
A, b = linear_eq_to_matrix(eqs, unknowns)
result = list(linsolve((A,b), unknowns))[0]
# Filling the numerical array with the sympy result
n_X[0] = result[unknowns.index(s_X_1)]
n_X[1] = result[unknowns.index(s_X_2)]
n_Y[0] = result[unknowns.index(s_Y_1)]
n_Y[1] = result[unknowns.index(s_Y_2)]

但是，如果我已经拥有 A 和 B，我可以这样做：

A_num = np.array(A).astype(np.float64)
b_num = np.array(b).astype(np.float64)
result_num = np.linalg.solve(A_num, b_num)

衡量绩效：

%%timeit -r 100 -n 100
result_sym = list(linsolve((A,b), unknowns))[0]
1.95 ms ± 285 µs per loop (mean ± std. dev. of 100 runs, 100 loops each)

%timeit -r 100 -n 100
result_num = np.linalg.solve(A_num, b_num)
The slowest run took 9.04 times longer than the fastest. This could mean that an intermediate result is being cached.
15.5 µs ± 10.3 µs per loop (mean ± std. dev. of 100 runs, 100 loops each)

%%timeit -r 100 -n 100
A, b = linear_eq_to_matrix(eqs, unknowns)
A_num = np.array(A).astype(np.float64)
b_num = np.array(b).astype(np.float64)
441 µs ± 112 µs per loop (mean ± std. dev. of 100 runs, 100 loops each)

如您所见，sympy 解决方案比 numpy 解决方案慢几个数量级。同样，将我们的数值系统转换为矩阵形式，然后将其转换为 numpy 数组以使用numpy.linalg.solve 的操作比求解系统本身要慢得多！

Answer 1

首先，我假设您的输入是一个符号方程组，其中所有系数和独立项都是浮点数，就像您在 MWE 中显示的那样，但要大得多。

在这种情况下，也许这个方法对你有用。

from sympy import *
import numpy as np

dim = 4
# Getting coefficients
a = np.random.uniform(low = -3, high = 3, size= (dim,dim))
b = np.random.uniform(low = -3, high = 3, size= (dim,))
# Numerical array for solutions
n_X = np.empty(shape = (2,))
n_Y = np.empty(shape = (2,))

# Defining unknowns

s_X_1 = Symbol('X_1')
s_Y_1 = Symbol('Y_1')
s_X_2 = Symbol('X_2')
s_Y_2 = Symbol('Y_2')

#Defining equations
eq1 = s_X_1*a[0,0] + s_Y_1*a[1,0] + s_X_2*a[2,0] + s_Y_2*a[3,0] - b[0]
eq2 = s_X_1*a[0,1] + s_Y_1*0 + s_X_2*a[2,1] + s_Y_2*a[3,1]- b[1]
eq3 = s_X_1*a[0,2] + s_Y_1*a[1,2] + s_X_2*a[2,2] + s_Y_2*a[3,2]- b[2]
eq4 = s_X_1*a[0,3] + s_Y_1*a[1,3] + s_X_2*0 + s_Y_2*a[3,3]- b[3]

# List of equations and unknowns

eqs = [eq1,eq2,eq3,eq4]
unknowns = ["X_1","Y_1","X_2","Y_2"]

b = []                                      # Independent terms
A = []                                      # Coefficients matrix

for eq in eqs:

    a = [0.]*len(unknowns)                  # List of coefficients for equation

    # Equation to list of strings

    l = str(eq).strip().replace(" + ","*").replace(" - ","*-").split("*")

    # Dictionary where keys are variables and values coefficients ( X_1: a[0,0] )

    dic = {l[i+1]:l[i] for i in range(0,len(l)-1,2)}
    for key in dic:
        a[unknowns.index(key)] = float(dic[key])

    # Add elements to A and b

    b.append(float(l[-1]))
    A.append(a)

# Solving the system

x = np.linalg.solve(A,b)

它基本上将每个方程转换为一个字符串，然后创建一个字典，其中键是未知数，值是系数，因此可以轻松识别它们。

为了检查该方法是否比仅使用 sympy 更快，我编写了一个程序，该程序允许我创建相同的示例，但系统很大，而且除了需要大量时间阅读所有内容之外符号方程（如果您已经在使用它们，这应该不是问题），它会创建 A 和 b 并在大约 3 秒内计算出 100x100 案例的解（sympy 大约花费 13 秒）。

我把代码贴在这里，以防你想自己检查。

with open("out.py","w") as file:

    dim = 100

    file.write("from sympy import * \n")
    file.write("import time\n")
    file.write("import numpy as np \n\n")


    for i in range(0,dim):

        file.write('s_X_%d = symbols("X_%d")\n'%(i+1,i+1))

    file.write("\n")
    file.write("dim = %d \n\n"%dim)

    file.write("a = np.random.uniform(low = -3, high = 3, size= (dim,dim))\n")
    file.write("b = np.random.uniform(low = -3, high = 3, size= (dim,))\n\n")

    for i in range(dim):
        s = "eq%d = "%(i+1)
        for j in range(dim):

            s += " + s_X_%d*a[%d,%d]"%(j+1,j,i)

        s += " - b[%d]\n"%i
        file.write(s)

    st = "eqs = [eq1"
    #uk = 'unknowns = [s_X_1'                           # Uncomment for sympy
    uk = 'unknowns = ["X_1"'                            # Uncomment for my method

    for i in range(2,dim+1):
        st += ",eq%d"%i
        uk += ',"X_%d"'%i                               # Uncomment for my method
        #uk += ',s_X_%d'%i                              # Uncomment for sympy

    file.write("\n")
    file.write(st + "]\n")
    file.write(uk + "]\n")

    # This uses sympy
    """
    file.write("print('Starts running')\n")
    file.write("t0 = time.time()\n\n")
    file.write('A, b = linear_eq_to_matrix(eqs, unknowns)\n')
    file.write('result = list(linsolve((A,b), unknowns))[0]\n')
    file.write('n_X = np.empty(shape = (%d,))\n\n'%dim)

    for i in range(dim):
        file.write('n_X[%d] = result[unknowns.index(s_X_%d)]\n'%(i,i+1))

    file.write('print("Elapsed time:",time.time() - t0)')
    """
    # This uses my method
    file.write("\nb = []\nA = []\n")
    file.write('print("Starts running")\nt0 = time.time()\n')
    file.write("for eq in eqs:\n\n")
    file.write("\ta = [0]*len(unknowns)\n")
    file.write('\tl = str(eq).strip().replace(" + ","*").replace(" - ","*-").split("*")\n')
    file.write('\tdic = {l[i+1]:l[i] for i in range(0,len(l)-1,2)}\n\tfor key in dic:\n')
    file.write('\t\ta[unknowns.index(key)] = float(dic[key])\n')
    file.write('\tb.append(float(l[-1]))\n\tA.append(a)\n')
    file.write('x = np.linalg.solve(A,b)\n')
    file.write('\nprint("Elapsed time:",time.time() - t0)')

这是我在这里的第一个答案，所以希望它有用！

编辑

我意识到您实际上并不需要使用字典来存储系数值。可以这样做。

# Old code
dic = {l[i+1]:l[i] for i in range(0,len(l)-1,2)}
for key in dic:
    a[unknowns.index(key)] = float(dic[key])

# New code
for i in range(0,len(l)-1,2):
    a[unknowns.index(l[i+1])] = float(l[i])

也许节省大量时间无济于事，但更简单。