[Python知识库]线性代数Python计算：线性方程组的最小二乘解

给定?上无解线性方程组$\mathbf{A}\mathbf{x}=\mathbf{b}$，构造${\mathbf{A}}^{\text{T}}\mathbf{A}$及${\mathbf{A}}^{\text{T}}\mathbf{b}$，然后调用博文《线性方程组的通解》定义的mySolve函数，解方程组${\mathbf{A}}^{\text{T}}\mathbf{A}\mathbf{x}={\mathbf{A}}^{\text{T}}\mathbf{b}$。取任一特解${\mathbf{x}}_0$即为解线性方程组$\mathbf{A}\mathbf{x}=\mathbf{b}$的一个最小二乘解（即$\mathbf{A}$的列向量生成空间中距$\mathbf{b}$最近的向量）。
例1 用Python计算无解方程组$?\begin{array}{l}4x_1+2x_2?x_3=2\\ 3x_1?x_2+2x_3=10\\ 11x_1+3x_2=8\end{array}$的最小二乘解。

import numpy as np                              #导入numpy
from fractions import Fraction as F             #导入Fraction
np.set_printoptions(formatter=                  #设置输出数据格式
                    {'all':lambda x:str(F(x).limit_denominator())})
A=np.array([[4,2,-1],                           #设置系数矩阵A
            [3,-1,2],
            [11,3,0]],dtype='float')
b=np.array([2,10,8])                            #常数项量b
B=np.matmul(A.T,A)                              #A的转置与A的积
c=np.matmul(A.T,b.reshape(3,1))                 #A的转置与b的积
X=mySolve(B,c)                                  #解最小二乘方程组
print(X[:,0])

程序的第5~7行设置原方程组的系数矩阵A，第8行设置原方程组的常数向量b。第9行调用numpy的matmul函数计算${\mathbf{A}}^{\text{T}}\mathbf{A}$，存于B。第10行计算${\mathbf{A}}^{\text{T}}\mathbf{A}\mathbf{b}$，存于c。第11行调用函数mySolve（见博文《线性方程组的通解》）解方程组${\mathbf{A}}^{\text{T}}\mathbf{A}\mathbf{x}={\mathbf{A}}^{\text{T}}\mathbf{b}$，解集记为X。注意，X中第1列（X[:,0]）存储的是方程组的特解。运行程序，输出

[9/5 -18/5 0]

即原方程组$\mathbf{A}\mathbf{x}=\mathbf{b}$的最小二乘解为${\mathbf{x}}_0=?\begin{array}{c}\frac{9}{5}\\ ?\frac{18}{5}\\ 0\end{array}?$。它是${\mathbf{\alpha }}_1=?\begin{array}{c}4\\ 3\\ 1\end{array}?,{\mathbf{\alpha }}_2=?\begin{array}{c}2\\ ?1\\ 3\end{array}?,{\mathbf{\alpha }}_3=?\begin{array}{c}1\\ 2\\ 0\end{array}?$的生成子空间$W$中距$\mathbf{b}=?\begin{array}{c}2\\ 10\\ 8\end{array}?$最近的向量。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！