【最小二乘法】矩阵法解法及Python代码_使用矩阵分解求最小二乘法python

当线性方程组无解时，如何寻求一个最精确解，以贴合在坐标系中离散的坐标点。
给出三个笛卡尔坐标系下的坐标点：
$$（1,0）（2,0）（1,2）$$
找到一条直线$y=kx+b$，能够贴合这三点（显然不可能）。为表述方便，以 $x_1$代 $k$，$x_2$ 代 $b$，
得到以下方程组：
{ 1 × k + b = 0 2 × k + b = 0 1 × k + b = 2 ⇔ { 1 × x 1 + x 2 = 0 2 × x 1 + x 2 = 0 1 × x 1 + x 2 = 2 \left\{

$$\begin{array}{c} 1×k+b=0 \\ 2×k+b=0 \\ 1×k+b=2 \end{array}$$ \right. ⇔ \left\{ $$\begin{array}{c} 1×x_1+x_2=0 \\ 2×x_1+x_2=0 \\ 1×x_1+x_2=2 \end{array}$$ \right. ⎩ ⎨ ⎧​1×k+b=02×k+b=01×k+b=2​⇔⎩ ⎨ ⎧​1×x1​+x2​=02×x1​+x2​=01×x1​+x2​=2​
将线性方程组写成矩阵形式：
$$\begin{array}{c}A=\left[\begin{array}{cc}1& 1\\ 2& 1\\ 1& 1\end{array}\right]\end{array}\begin{array}{c}b=\left[\begin{array}{c}0\\ 0\\ 2\end{array}\right]\end{array}$$
得出矩阵符合关系式：$Ax=b$，而我们的目的则是求得一个原始列向量:
$$\begin{array}{c}x=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]\end{array}$$

分析线性方程组解的情况

（1）

根据矩阵 A = [ 1 1 2 1 1 1 ]

$$\begin{gathered} A= \begin{bmatrix} 1&1 \\ 2&1 \\ 1&1 \end{bmatrix} \end{gathered}$$ A=⎣ ⎡​121​111​⎦ ⎤​​列向量线性无关,可得$m=3,n=2,r=2$。
由$r<m,且n=r$可知：线性方程组无解，或存在唯一解。

（2）

由于向量b= [ 0 0 2 ]

$$\begin{gathered} \begin{bmatrix} 0 \\ 0 \\ 2 \end{bmatrix} \end{gathered}$$ ⎣ ⎡​002​⎦ ⎤​​，无法由$\begin{array}{c}A=\left[\begin{array}{cc}1& 1\\ 2& 1\\ 1& 1\end{array}\right]\end{array}$线性组合得到,即向量b不在列空间C(A)上。

综上所述可知，线性方程组无解，准确来说，无精确解。

寻找近似解

由线性方程组$Ax=b，即可以得到矩阵$ A = [ 1 1 2 1 1 1 ]

$$\begin{gathered} A= \begin{bmatrix} 1&1 \\ 2&1 \\ 1&1 \end{bmatrix} \end{gathered}$$ A=⎣ ⎡​121​111​⎦ ⎤​​的两个线性无关向量：
$$\begin{array}{c}{a}_1=\left[\begin{array}{c}1\\ 2\\ 1\end{array}\right]\end{array}\begin{array}{c}{a}_2=\left[\begin{array}{c}1\\ 1\\ 1\end{array}\right]\end{array}$$
这两个向量张成$R^3$空间中的一个二维空间，即矩阵A的列空间$C（A）$。由于向量$\begin{array}{c}b=\left[\begin{array}{c}0\\ 0\\ 2\end{array}\right]\end{array}$在列空间$C（A）$外，因此线性方程组$Ax=b$是无解的。

虽然线性方程组$Ax=b$不存在精确解，但我们可以寻求距离目标最近的近似解。根据上图，我们可以直观的思考到一个合理的方法——从向量$b$向二维空间即列空间$C(A)$上引垂线，得到$b$的投影向量$p$。由$p$和$b$可得误差向量$e=b-p$。
向量$b$经“某一矩阵”变换后，可得位于列空间$C(A)$上的投影向量$p$。而“某一矩阵”则称为“投影矩阵”。我们可以得到关系式：
$$p=Pb$$
引出投影向量$p$后，问题就清晰了，我们可以求出$p$向量来近似求解方程$A\hat{x}=p$对应 x = [ x 1 x 2 ]

$$\begin{gathered} x= \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \end{gathered}$$ x=[x1​x2​​]​。

求解的思路与步骤

由上图可知$p$向量的两个性质：
（1）
$p$向量位于列空间$C（A）$上，即$A\hat{x}=p$。
A x ^ = [ a 1 a 2 ] [ x 1 ^ x 2 ^ ] = p A\hat{x}=

$$\begin{gathered} \begin{bmatrix} a_1&a_2\end{bmatrix} \end{gathered}$$ $$\begin{gathered} \begin{bmatrix} \hat{x_1} \\ \hat{x_2} \end{bmatrix} \end{gathered}$$ =p Ax^=[a1​​a2​​]​[x1​^​x2​^​​]​=p
（2）
误差向量垂直于二维空间$C(A)$的任一组向量。
$$a_1・e=0和a_2・e=0$$
将$\{\begin{array}{c}e=b-p\\ p=A\hat{x}\end{array}$代入上述等式：
$$\{\begin{array}{c}{a}_{1}e=a_1(b-p)=a_1(b-A\hat{x})=0\\ a_{2}e=a_2(b-p)=a_2(b-A\hat{x})=0\end{array}$$

将向量内积转化为矩阵乘法

{ a 1 T ( b − A x ^ ) = 0 a 2 T ( b − A x ^ ) = 0 \left \{

$$\begin{array}{c} a^T_1(b-A\hat{x})=0 \\ a^T_2(b-A\hat{x})=0 \end{array}$$ \right. {a1T​(b−Ax^)=0a2T​(b−Ax^)=0​
整理得：
$$\begin{array}{c}\left[\begin{array}{c}{a}_1^T\\ a_2^T\end{array}\right](b-A\hat{x})=0\end{array}$$
$$\begin{array}{c}A=\left[\begin{array}{cc}{a}_1& a_2\end{array}\right]\end{array}则\begin{array}{c}\left[\begin{array}{c}{a}_1^T\\ a_2^T\end{array}\right]=A^T\end{array}\Rightarrow$$
$$\hat{x}=(A^{T}A{)}^{-1}{A}^{T}b$$
将上式代入$p=A\hat{x}$,得：
$$p=A・(A^{T}A{)}^{-1}{A}^{T}b$$
代入$p=Pb$，得：
$$P=A・(A^{T}A{)}^{-1}A$$
公式汇总：
$$\{\begin{array}{c}\hat{x}=(A^{T}A{)}^{-1}{A}^{T}b\\ p=A・(A^{T}A{)}^{-1}{A}^{T}b\\ P=A・(A^{T}A{)}^{-1}A\end{array}$$

结论

最小二乘法的几何意义是高维空间中的一个向量在低维空间上的投影。

Python代码求解

代码

import numpy as np
from scipy import linalg
A=np.array([[1,1],
            [2,1],
            [1,1]])
b=np.array([0,0,2])
A_T_A=np.dot(A.T,A)
A_T_A_in=linalg.inv(A_T_A)
x_hat=np.dot(np.dot(A_T_A_in,A.T),b)
p=np.dot(A,x_hat)
p_matrix=np.dot(A,np.dot(A_T_A_in,A.T))
print("近似解x为：")
print(x_hat)
print("\n投影向量p为：")
print(p)
print("\n投影矩阵P为：")
print(p_matrix)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

运行结果

近似解x为：
[-1.  2.]    

投影向量p为：
[1. 0. 1.]   

投影矩阵P为：
[[0.5 0.  0.5]
 [0.  1.  0. ]
 [0.5 0.  0.5]]
1
2
3
4
5
6
7
8
9
10