Python线性方程求解-矩阵左除“\“、右除“/“

目录

1 线性方程组求解方法

2 左除“\”→Ax=B

3 右除"/"→xA=B

4 其它说明

---

1 线性方程组求解方法

• 如果Ax=B，则x=A\B，称为左除；
• 如果xA=B，则x=B/A，称为右除。

        式中x为未知数。一般情况下，左除用的系比较多一些。在matlab里面实现左除或者右除会比较简单，直接有运算符号"\"和"/"。但是在Python里面就不能直接采用运算符号：

• Python里面"\"不是一个运算符号；
• Python直接采用B/A，表示的是矩阵B的每个元素除以矩阵A的每个元素，这并不是矩阵运算。

那在Python里面该如何实现矩阵的除法运算呢？

2 左除“\”→Ax=B

（1）当矩阵A是方阵，注：A的行和B的行相等

① 采用inv()函数，即：

② 采用solve()函数，即：

举例：

import numpy as np
A=np.array([[1,5,3],[4,8,6],[7,10,9]])
B=np.array([[6],[8],[10]])
inv_A=np.linalg.inv(A)
x_1=np.dot(inv_A,B)
x_2=np.linalg.solve(A,B)
print(f'A={A}')
print(f'B={B}')
print(f'x_1={x_1}')
print(f'x_2={x_2}')

运行结果：

A=[[ 1  5  3]
 [ 4  8  6]
 [ 7 10  9]]
B=[[ 6]
 [ 8]
 [10]]
x_1=[[-2.00000000e+00]
 [-1.77635684e-15]
 [ 2.66666667e+00]]
x_2=[[-2.00000000e+00]
 [ 1.49213975e-15]
 [ 2.66666667e+00]]

注：当A并不是方阵时，采用inv()和solve()会报错，无法求解。

（2）当A不是方阵，注：A的行和B的行相等

① 采用pinv()函数求解，即：

② 采用lstsq()函数求解：

语法：numpy.linalg.lstsq(A,B,rcond=“warn”)

• A是一个M行N列的系数矩阵；
• B是一个(M,)或者(M,K)，如果b是一个M行K列的二维矩阵，函数会逐个计算每一列的最小二乘法；
• rcond这个参数是可选的，是用于奇异矩阵的处理的，官方推荐我们一般用 rcond=None；
• 返回值：返回值的第一个元素即为我们想要的结果，所以一般的用法是lstsq()[0]；

举例说明：

import numpy as np
A=np.array([[1,8,5,6],[4,6,7,8],[9,10,12,14]])
B=np.array([[6],[8],[10]])
x_1=np.dot(np.linalg.pinv(A),B)
x_2=np.linalg.lstsq(A,B,rcond=None)[0]
print(f'A={A}')
print(f'B={B}')
print(f'x_1={x_1}')
print(f'x_2={x_2}')

运行结果：

A=[[ 1  8  5]
 [ 4  6  7]
 [ 9 10 12]
 [ 4 10 14]]
B=[[ 6]
 [ 8]
 [10]
 [ 9]]
x_1=[[0.32212839]
 [0.68526414]
 [0.08256129]]
x_2=[[0.32212839]
 [0.68526414]
 [0.08256129]]

3 右除"/"→xA=B

（1）B是方阵，注：A的列和B的列相等。

① 采用inv()函数实现，即：

② 采用solve()函数，即

举例：

import numpy as np
B=np.array([[1,5,3],[4,8,6],[7,10,9]])
A=np.array([6,8,10])
inv_B=np.linalg.inv(B)
x_1=np.dot(A,inv_B)
x_2=np.linalg.solve(B.T,A.T).T
print(f'A={A}')
print(f'B={B}')
print(f'x_1={x_1}')
print(f'x_2={x_2}')

运行结果：

A=[ 6  8 10]
B=[[ 1  5  3]
 [ 4  8  6]
 [ 7 10  9]]
x_1=[  6.66666667 -10.66666667   6.        ]
x_2=[  6.66666667 -10.66666667   6.        ]

（2）B不是方阵，注：A的列和B的列相等。

采用pinv()函数实现，即

$B/A=B*A^-^1=B*pinv(A)$

举例：

import numpy as np
B=np.array([[1,5],[4,8],[7,10]])
A=np.array([6,8])      #此时的A是(3,),是一维
A=A[:,np.newaxis].T    #此时的A是(3,1),是二维，必须是二维的才能用pinv()计算
print(f'A={A}')
print(f'B={B}')
x_1=np.dot(B,np.linalg.pinv(A))
print(f'A={A}')
print(f'B={B}')
print(f'x_1={x_1}')

运行结果（和matlab运行结果一样）：

A=[[6 8]]
B=[[ 1  5]
 [ 4  8]
 [ 7 10]]
x_1=[[0.46]
 [0.88]
 [1.22]]

4 其它说明

① pinv能够求解方阵，只是运算代价更大一点；

② 设矩阵的行数为m列数为n，对于左除(A\B)，当m≥n时，Python的求解结果和Matlab的求解结果一样；但是当m<n时，Python的求解结果不一样，原因是Matlab返回的解是尽可能多0值的解，Python返回得解是最小二乘解（范数最小）。不过这对于算法的应用区别不大，想返回何种解取决于实际问题，结果是否满足精度需求。