最小二乘拟合的三种算法推导及Python代码

1 最小二乘拟合（方法一）

本章介绍的是我在上研究生课程《数值分析》时所学的内容，也是最一般的解法，可以对任意函数进行拟合。

1.1 数学推导

对一组数据$(x_i,y_i),i=1,2\cdots ,m$，要在某个函数类 Φ = { ϕ 0 ( x ) , ϕ 1 ( x ) , ⋯   , ϕ n ( x ) } , n ≪ m \Phi=\{\phi_0(x),\phi_1(x),\cdots,\phi_n(x)\},n\ll m Φ={ϕ0​(x),ϕ1​(x),⋯,ϕn​(x)},n≪m中构造一个函数$\varphi (x)={\sum }_{k=0}^n{a}_k{\varphi }_k(x)$，使得$I={\sum }_{i=1}^m{\left[\varphi (x_i)-y_i\right]}^2$取得极小值。
此处的函数类$\Phi$就是指很多种类函数的集合，例如幂函数、三角函数、指数函数、有理函数、多项式函数等。例如，常用的多项式函数类有 Φ 1 = { 1 , x , x 2 , ⋯   , x n } \Phi_1=\{1,x,x^2,\cdots,x^n\} Φ1​={1,x,x2,⋯,xn}， Φ 2 = { 1 , x , x , ⋯   , x ⏟ n } \Phi_2=\{1,\underbrace{x,x,\cdots,x}_{n}\} Φ2​={1,n x,x,⋯,x​​}。其中${\Phi }_1$为n阶多项式，${\Phi }_2$为n元一阶多项式，分别对应着多项式拟合和多维线性回归。
残差函数$I$定义为
$$I(a_0,a_1,\cdots ,a_n)=\sum _{i=1}^m{\left[a_0{\varphi }_0(x_i)+a_1{\varphi }_1(x_i)+\cdots +a_n{\varphi }_n(x_i)-y_i\right]}^2$$
令$I$对系数$a_k$的偏导为0，即可求出使残差函数$I$最小的系数$a_k$的值
$$\begin{array}{rl}\frac{\partial I}{\partial a_k}& =2\sum _{i=1}^m\left[\varphi (x_i)-y_i\right]\frac{\partial \varphi (x_i)}{\partial a_k}\\ & =2\sum _{i=1}^m\left[2\sum _{j=0}^n{a}_j{\varphi }_j(x_i)-y_i\right]{\varphi }_k(x_i)\\ & =2\left\{\sum _{j=0}^n{a}_j\left[\sum _{i=1}^m{\varphi }_j(x_i){\varphi }_k(x_i)\right]-\sum _{i=1}^m{y}_i{\varphi }_k(x_i)\right\}\end{array}$$
简便起见，记
$$\begin{array}{rl}({\varphi }_j,{\varphi }_k)& =\sum _{i=1}^m{\varphi }_j(x_i){\varphi }_k(x_i)\\ (y,{\varphi }_k)& =\sum _{i=1}^m{y}_i{\varphi }_k(x_i)\end{array}$$
则$\frac{\partial I}{\partial a_k}$可写成如下形式
$$\frac{\partial I}{\partial a_k}=2\left[\sum _{j=0}^n{a}_j({\varphi }_j,{\varphi }_k)-(y,{\varphi }_k)\right]=0$$
化简为
$$\sum _{j=0}^n{a}_j({\varphi }_j,{\varphi }_k)=(y,{\varphi }_k)$$
注意到，公式中的系数$a_k$中$k=0,1,\cdots ,n$，则实际上的操作步骤为分别对每个系数求偏导，并令偏导为0，从而求出使残差函数$I$取到极小值的所有的系数。写成矩阵的形式
$$\left[\begin{array}{cccc}({\varphi }_0,{\varphi }_0)& ({\varphi }_0,{\varphi }_1)& \cdots & ({\varphi }_0,{\varphi }_n)\\ ({\varphi }_1,{\varphi }_0)& ({\varphi }_1,{\varphi }_1)& \cdots & ({\varphi }_1,{\varphi }_n)\\ ⋮& ⋮& & ⋮\\ ({\varphi }_n,{\varphi }_0)& ({\varphi }_n,{\varphi }_1)& \cdots & ({\varphi }_n,{\varphi }_n)\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n\end{array}\right]=\left[\begin{array}{c}(y,{\varphi }_0)\\ (y,{\varphi }_1)\\ ⋮\\ (y,{\varphi }_n)\end{array}\right]$$
上式方程组称为法方程组（Normal Equation），向量$\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n\end{array}\right]$称为回归系数（Regression Coefficients）。用更一般的形式将上式再次改写
$$Ax=b$$
则可通过矩阵求逆的方式求出回归系数
$$x=A^{-1}b$$

1.2 算例

利用一阶多项式 Φ 1 = { 1 , x } \Phi_1=\{1,x\} Φ1​={1,x}通过1.1节的方法对如下数据进行最小二乘拟合。

Independent variable x	Dependent variable y
0	394.33
4	329.50
8	291.00
12	255.17
16	229.33
20	204.83
24	179.00
28	163.83
32	150.33

则法方程组为
$$\left[\begin{array}{cc}9& \sum x_i\\ \sum x_i& \sum x_i^2\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\end{array}\right]=\left[\begin{array}{c}\sum y_i\\ \sum x_i{y}_i\end{array}\right]$$
带入数据
$$\left[\begin{array}{cc}9& 144\\ 144& 3264\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\end{array}\right]=\left[\begin{array}{c}2197.32\\ 28167.72\end{array}\right]$$
求解得到
$$\left[\begin{array}{c}{a}_0\\ a_1\end{array}\right]=\left[\begin{array}{c}360.636667\\ -7.280625\end{array}\right]$$
则可得到线性回归方程
$$y=-7.280625x+360.636667$$
拟合结果如下

1.3 Python 代码

#%% import neccesary packages

import numpy as np
import matplotlib.pyplot as plt
#%% generate the dataset

x = np.arange(0,33,4)
y = np.array([394.33, 329.50, 291.00, 255.17, 229.33, 204.83, 179.00, 163.83, 150.33])
print(x)
#%% calculate the normal equation

A11 = 9
A12 = np.sum(x)
A21 = A12
A22 = np.sum(np.power(x,2))
A = np.array([[A11,A12],[A21,A22]])
print(A)

b1 = np.sum(y)
b2 = np.sum(x*y)
b=np.array([b1, b2])
print(b)
#%% solve the polynamial coefficient a0 & a1

a = np.linalg.inv(A).dot(b)
a0 = a[0]
a1 = a[1]
print(a)
#%% plot the picture
fig, ax = plt.subplots(1, 1, figsize=(6, 3))
ax.scatter(x, y, color='black')
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
ax.plot(x, a0 + a1 * x)
plt.show()

---

2.最小二乘拟合（方法二）

本章介绍我自己的推导思路，本文以多项式拟合为例说明思路，其他函数的拟合思路相同。

2.1 数学推导

对一组数据$(x_i,y_i),i=1,2\cdots ,m$进行多项式拟合，假设多项式为
$$f(x)=a_n{x}^n+a_{n-1}{x}^{n-1}+\cdots +a_{1}x+a_0$$
对任意数据点$(x_k,y_k)$，它与多项式$f(x)$在$y$轴方向上的距离称为残差${\delta }_k$，定义为
$${\delta }_k=f(x_k)-y_k$$
拟合的目标是寻找一组多项式系数 { a 0 , a 1 , ⋯   , a n } \{a_0,a_1,\cdots,a_n\} {a0​,a1​,⋯,an​}使得残差函数$I$取最小值，本章将直接利用最小二乘法求解。将每个数据点的残差写成矩阵形式
$$\left[\begin{array}{ccccc}1& x_1& x_1^2& \cdots & x_1^n\\ 1& x_2& x_2^2& \cdots & x_2^n\\ ⋮& ⋮& ⋮& & ⋮\\ 1& x_m& x_m^2& \cdots & x_m^n\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n\end{array}\right]-\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_m\end{array}\right]=\left[\begin{array}{c}{\delta }_1\\ {\delta }_2\\ ⋮\\ {\delta }_m\end{array}\right]$$
定义残差函数$I(a_0,a_1,\cdots ,a_n)$为
$$I=\parallel {\delta }_i{\parallel }_2^2=\parallel f(x_i)-y_i{\parallel }_2^2$$
要使残差函数$I$最小，本质上是求如下方程的最小二乘解
$$\left[\begin{array}{ccccc}1& x_1& x_1^2& \cdots & x_1^n\\ 1& x_2& x_2^2& \cdots & x_2^n\\ ⋮& ⋮& ⋮& & ⋮\\ 1& x_m& x_m^2& \cdots & x_m^n\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n\end{array}\right]=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_m\end{array}\right]$$
上式可写成一般形式
$$Ax=b$$
（1）当$A$为方阵时，其最小二乘解为
$$x=-(A^{T}A{)}^{-1}{A}^{T}b$$
（2）在大多数情况下，$A$不为方阵，可用Moore-Penrose逆$A^{+}$求解
定理1： 设$A\in {\mathbb{C}}^{m\times n}$是秩为$r(r>0)$的矩阵，且$A$的满秩分解为
$$A=FG(F\in {\mathbb{C}}^{m\times r},G\in {\mathbb{C}}^{r\times n}秩都为r)$$
则
$$A^{+}=G^H(G{G}^H{)}^{-1}(F^{H}F{)}^{-1}{F}^H$$
定理2： 设$A\in {\mathbb{C}}^{m\times n},b\in {\mathbb{C}}^m$有解的充要条件是
$$A{A}^{+}b=b$$
其通解为
$$x=A^{+}b+(I-A^{+}A)y$$
其唯一极小范数最小二乘解为
$$x=A^{+}b$$

2.2 算例

本章算例与1.2节相同，以作对照。
利用一阶多项式$y=a_0+a_{1}x$对如下数据进行最小二乘拟合。

Independent variable x	Dependent variable y
0	394.33
4	329.50
8	291.00
12	255.17
16	229.33
20	204.83
24	179.00
28	163.83
32	150.33

列写$Ax=b$的一般形式
$$\left[\begin{array}{cc}1& 0\\ 1& 4\\ 1& 8\\ 1& 12\\ 1& 16\\ 1& 20\\ 1& 24\\ 1& 28\\ 1& 32\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\end{array}\right]=\left[\begin{array}{c}394.33\\ 329.50\\ 291.00\\ 255.17\\ 229.33\\ 204.83\\ 179.00\\ 163.83\\ 150.33\end{array}\right]$$
利用$Python$命令$numpy.linalg.pinv()$求解矩阵$A$的$A^{+}$逆，并求其唯一极小范数最小二乘解
$$\left[\begin{array}{c}{a}_0\\ a_1\end{array}\right]=\left[\begin{array}{c}360.636667\\ -7.280625\end{array}\right]$$
求解结果与第一章相同。

2.3 Python 代码

#%% import neccesary packages

import numpy as np
import matplotlib.pyplot as plt
#%% generate the dataset

x = np.arange(0,33,4)
y = np.array([394.33, 329.50, 291.00, 255.17, 229.33, 204.83, 179.00, 163.83, 150.33])
print(x)
#%% solve the polynamial coefficient a0 & a1

A = np.vstack((np.array([1]*9), x))
print(A)
pinv_A = np.linalg.pinv(A).T
print(pinv_A)
a = pinv_A.dot(y)
print(a)
a0 = a[0]
a1 = a[1]
#%% plot the picture

fig, ax = plt.subplots(1, 1, figsize=(6, 3))
ax.scatter(x, y, color='black')
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
ax.plot(x, a0 + a1 * x)
plt.show()

---

3 最小二乘法拟合（方法三）

本章重点介绍线性回归算法。其实该算法本质上和第一章方法一样，只不过将其限定为一阶多项式，推导出更直观的表达方法。

3.1 数学推导

假设线性回归模型如下
$$y_i={\beta }_0+{\beta }_1{x}_i+{ϵ}_i$$
式中${ϵ}_i$为模型噪声服从高斯分布：$ϵ\sim N(0,{\sigma }^2)$。
假设需要对一组数据$(x_i,y_i),i=1,2\cdots ,m$进行线性回归，定义残差函数
$$Q=\sum _{i=1}^n{\left[y_i-({\beta }_0+{\beta }_1{x}_i)\right]}^2$$
优化目标为最小化残差函数
$$\underset{{\beta }_0,{\beta }_1}{\min }:Q$$
和第一章方法一样，令残差函数$Q$对系数${\beta }_0、{\beta }_1$的偏导为0
$$\begin{array}{rl}\frac{\partial Q}{\partial {\beta }_0}=-2\sum _{i=1}^n\left[y_i-({\beta }_0+{\beta }_1{x}_i)\right]& =0\\ \Rightarrow n{\beta }_0+{\beta }_1\sum _{i=1}^n{x}_i& =\sum _{i=1}^n{y}_i(1)\end{array}$$
$$\begin{array}{rl}\frac{\partial Q}{\partial {\beta }_1}=-2\sum _{i=1}^n{x}_i\left[y_i-({\beta }_0+{\beta }_1{x}_i)\right]& =0\\ \Rightarrow {\beta }_0\sum _{i=1}^n{x}_i+{\beta }_1\sum _{i=1}^n{x}_i^2& ={\beta }_1\sum _{i=1}^n{x}_i{y}_i(2)\end{array}$$
联立式（1）和式（2）解方程得到${\beta }_0、{\beta }_1$的代数表达式
$${\hat{\beta }}_0=\frac{({\sum }_{i=1}^n{x}_i^2)({\sum }_{i=1}^n{y}_i)-({\sum }_{i=1}^n{x}_i)({\sum }_{i=1}^n{x}_i{y}_i)}{n{\sum }_{i=1}^n{x}_i^2-({\sum }_{i=1}^n{x}_i{)}^2}$$
$${\hat{\beta }}_1=\frac{{\sum }_{i=1}^n{x}_i{y}_i-({\sum }_{i=1}^n{y}_i)({\sum }_{i=1}^n{x}_i)}{n{\sum }_{i=1}^n{x}_i^2-({\sum }_{i=1}^n{x}_i{)}^2}$$
为表达简便，记
$$S_{xy}=\sum _{i=1}^n{x}_i{y}_i-\frac{1}{n}(\sum _{i=1}^n{x}_i)(\sum _{i=1}^n{y}_i)=\sum _{i=1}^n(x_i-\overline{x})(y_i-\overline{y})$$
$$S_{xx}=\sum _{i=1}^n{x}_i^2-\frac{1}{n}(\sum _{i=1}^n{x}_i{)}^2=\sum _{i=1}^n(x_i-\overline{x}{)}^2$$
$$S_{yy}=\sum _{i=1}^n{y}_i^2-\frac{1}{n}(\sum _{i=1}^n{y}_i{)}^2=\sum _{i=1}^n(y_i-\overline{y}{)}^2$$
其中，$\overline{x}、\overline{y}$分别为拟合数据$x_i、y_i$的平均值，$i=1,2,\cdots ,m$
$$\begin{array}{rl}\overline{x}& =\frac{1}{n}\sum _{i=1}^n{x}_i\\ \overline{y}& =\frac{1}{n}\sum _{i=1}^n{y}_i\end{array}$$
则${\beta }_0、{\beta }_1$的表达式可改写为
$$\begin{array}{rl}\widehat{{\beta }_0}& =\overline{y}-{\beta }_1\overline{x}\\ \widehat{{\beta }_1}& =\frac{S_{xy}}{S_{xx}}\end{array}$$

3.2 算例

本章算例与1.2节相同，以作对照。
对如下数据进行线性回归。

Independent variable x	Dependent variable y
0	394.33
4	329.50
8	291.00
12	255.17
16	229.33
20	204.83
24	179.00
28	163.83
32	150.33

计算系数${\hat{\beta }}_0、{\hat{\beta }}_1$
$$\begin{array}{rl}\overline{x}& =16\\ \overline{y}& =244.15\\ S_{xy}& =-6989.40\\ S_{xx}& =960\\ \widehat{{\beta }_1}& =\frac{S_{xy}}{S_{xx}}=-7.280625\\ \widehat{{\beta }_0}& =\overline{y}-{\beta }_1\overline{x}=360.636667\end{array}$$
计算结果和第一章、第二章结果一致。
$$y=360.636667-7.280625x$$

3.3 Python 代码

#%% import neccesary packages

import numpy as np
import matplotlib.pyplot as plt
#%% generate the dataset

x = np.arange(0,33,4)
y = np.array([394.33, 329.50, 291.00, 255.17, 229.33, 204.83, 179.00, 163.83, 150.33])
print(x)
#%%

x_bar = np.sum(x)/9
y_bar = np.sum(y)/9
Sxy = np.sum((x-x_bar)*(y-y_bar))
Sxx = np.sum(np.power(x-x_bar,2))
beta_1 = Sxy/Sxx
beta_0 = y_bar - x_bar * beta_1
print(beta_0)
print(beta_1)
#%%

fig, ax = plt.subplots(1, 1, figsize=(6, 3))
ax.scatter(x, y, color='black')
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
ax.plot(x, beta_0 + beta_1 * x)
plt.show()

---

4 利用sklearn.linear_model()

4.1 参考资料

官网链接：Scikit-learn官网链接
教学1：python-sklearn学习笔记 第一节 linear_model
教学2：numpy中newaxis函数的基本用法及其理解（傻瓜版）

4.2 Python 代码

#%% import neccesary packages

import numpy as np
import matplotlib.pyplot as plt
import sklearn.linear_model as lm

#%% generate the dataset

x = np.arange(0,33,4)
y = np.array([394.33, 329.50, 291.00, 255.17, 229.33, 204.83, 179.00, 163.83, 150.33])
print(x)
#%%

x_lr = np.linspace(0,40,500)
lr = lm.LinearRegression()
lr.fit(x[:, np.newaxis],y)
y_lr = lr.predict(x_lr[:, np.newaxis])
#%%

fig, ax = plt.subplots(1, 1, figsize=(6, 3))
ax.scatter(x, y, color='black')
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
ax.plot(x_lr, y_lr)