【数学基础】最小二乘法

【fishing-pan：https://blog.csdn.net/u013921430 转载请注明出处】

前言

当我们有NNN组数据，希望能够用一个函数来拟合这组数据的分布情况时，首先想到的就是最小二乘法。最小二乘法是线性回归中最基础也是最常用模型，它也可以用来拟合高次多项式。

什么是最小二乘法

那么，什么是最小二乘法呢？将我们手中已有的数据表示成一个集合，（为了简化模型，暂且只考虑输入属性只有xxx一个的情况）D={(x1;y1),(x2;y2),⋯ ,(xn;yn)}D=\left \{ (x_{1};y_{1}),(x_{2};y_{2}),\cdots ,(x_{n};y_{n}) \right \}D={(x1;y1),(x2;y2),⋯,(xn;yn)}
首先假设最终得到的线性模型为；
y=ax+by=ax+by=ax+b

如上图所示，绿色的线条为所求的回归模型，黑色点为实际样本输入数据，红色线段为实际值与模型预测值之间的欧氏距离，用以下公式表示；
ei=∥axi+b−yi∥e_{i}=\left \| ax_{i}+b-y_{i} \right \|ei=∥axi+b−yi∥
为了让最终所求的模型尽量准确地表现数据分布规律，我们需要让eie_{i}ei的和尽可能的小，所以(a∗,b∗)(a^{*},b^{*})(a∗,b∗)在eie_{i}ei的和取最小值时取得；
(a∗,b∗)=argmin(a,b)∑i=1n∥axi+b−yi∥(a^{*},b^{*})=\underset{(a,b)}{argmin}\sum_{i=1}^{n}\left \| ax_{i}+b-y_{i} \right \|(a∗,b∗)=(a,b)argmini=1∑n∥axi+b−yi∥
但是上述式子中的绝对值是不可导的，求解困难。因此采用均方误差来代替绝对值。
(a∗,b∗)=argmin(a,b)∑i=1n(axi+b−yi)2(a^{*},b^{*})=\underset{(a,b)}{argmin}\sum_{i=1}^{n} (ax_{i}+b-y_{i} )^{2}(a∗,b∗)=(a,b)argmini=1∑n(axi+b−yi)2
这种基于均方误差来进行模型求解的方法就叫最小二乘法，也叫最小平方法。

参数推导

首先还是来看上面提到的最简单的情况。令；
E(a,b)=∑i=1n(axi+b−yi)2E(a,b)=\sum_{i=1}^{n} (ax_{i}+b-y_{i} )^{2}E(a,b)=i=1∑n(axi+b−yi)2
EEE是关于aaa与bbb的二次函数，只会有一个极值点。而很明显，让EEE取极大值的(a,b)(a,b)(a,b)出现在空间中无限远处的，是不存在的。所以当(a,b)(a,b)(a,b)让EEE取得极值时，EEE取得极小值，分别对aaa、bbb求偏导。
∂E∂a=2∑i=1n(axi+b−yi)xi=2(∑i=1naxi2+∑i=1nbxi−∑i=1nyixi)=2(∑i=1naxi2+b∗nxˉ−∑i=1nyixi)\begin{aligned} \frac{\partial E}{\partial a}&=2\sum_{i=1}^{n} (ax_{i}+b-y_{i} )x_{i}\\ &=2(\sum_{i=1}^{n} ax_{i}^{2}+\sum_{i=1}^{n}bx_{i}-\sum_{i=1}^{n}y_{i}x_{i})\\&=2(\sum_{i=1}^{n} ax_{i}^{2}+b*n\bar{x}-\sum_{i=1}^{n}y_{i}x_{i}) \end{aligned} ∂a∂E=2i=1∑n(axi+b−yi)xi=2(i=1∑naxi2+i=1∑nbxi−i=1∑nyixi)=2(i=1∑naxi2+b∗nxˉ−i=1∑nyixi)
∂E∂b=2∑i=1n(axi+b−yi)=2(a∗nxˉ+nb−nyˉ)\frac{\partial E}{\partial b}=2\sum_{i=1}^{n} (ax_{i}+b-y_{i} )=2(a*n\bar{x}+nb-n\bar{y}) ∂b∂E=2i=1∑n(axi+b−yi)=2(a∗nxˉ+nb−nyˉ)
令；
∂E∂a=0，∂E∂b=0\frac{\partial E}{\partial a}=0，\frac{\partial E}{\partial b}=0∂a∂E=0，∂b∂E=0
可得；
a=∑i=1nxiyi−nxˉyˉ∑i=1nxi2−nxˉ2a=\frac{\sum_{i=1}^{n}x_{i}y_{i}-n\bar{x}\bar{y}}{\sum_{i=1}^{n}x_{i}^{2}-n\bar{x}^{2}}a=∑i=1nxi2−nxˉ2∑i=1nxiyi−nxˉyˉ
b=yˉ−axˉ=yˉ∑i=1nxi2−xˉ∑i=1nxiyi∑i=1nxi2−nxˉ2b=\bar{y}-a\bar{x}=\frac{\bar{y}\sum_{i=1}^{n}x_{i}^{2}-\bar{x}\sum_{i=1}^{n}x_{i}y_{i}}{\sum_{i=1}^{n}x_{i}^{2}-n\bar{x}^{2}} b=yˉ−axˉ=∑i=1nxi2−nxˉ2yˉ∑i=1nxi2−xˉ∑i=1nxiyi
由于，
∑i=1n(xi−xˉ)(yi−yˉ)=∑i=1n(xiyi−xiyˉ−yixˉ+xˉyˉ)=∑i=1nxiyi−yˉ∑i=1nxi−xˉ∑i=1nyi+nxˉyˉ=∑i=1nxiyi−nxˉyˉ\begin{aligned} \sum_{i=1}^{n}(x_{i}-\bar{x})(y_{i}-\bar{y})&=\sum_{i=1}^{n}(x_{i}y_{i}-x_{i}\bar{y}-y_{i}\bar{x}+\bar{x}\bar{y})\\ &=\sum_{i=1}^{n}x_{i}y_{i}-\bar{y}\sum_{i=1}^{n}x_{i}-\bar{x}\sum_{i=1}^{n}y_{i}+n\bar{x}\bar{y}\\ &=\sum_{i=1}^{n}x_{i}y_{i}-n\bar{x}\bar{y} \end{aligned} i=1∑n(xi−xˉ)(yi−yˉ)=i=1∑n(xiyi−xiyˉ−yixˉ+xˉyˉ)=i=1∑nxiyi−yˉi=1∑nxi−xˉi=1∑nyi+nxˉyˉ=i=1∑nxiyi−nxˉyˉ
同理，
∑i=1n(xi−xˉ)2=∑i=1nxi2−nxˉ2\sum_{i=1}^{n}(x_{i}-\bar{x})^{2}=\sum_{i=1}^{n}x_{i}^{2}-n\bar{x}^{2} i=1∑n(xi−xˉ)2=i=1∑nxi2−nxˉ2
所以aaa还可以表示成，
a=∑i=1n(xi−xˉ)(yi−yˉ)∑i=1n(xi−xˉ)2a=\frac{\sum_{i=1}^{n}(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sum_{i=1}^{n}(x_{i}-\bar{x})^{2}} a=∑i=1n(xi−xˉ)2∑i=1n(xi−xˉ)(yi−yˉ)

多元线性回归

上述的推导中输入属性数目为一，在一般的情况下，样本的属性会有多个，所以通过矩阵的形式将上述推导过程进行推广。假设数据集DDD中样本属性有ddd个；将第iii个样本可以表示成以下形式；
(xi1;xi2;⋯ ;xid;y)(x_{i}^{1};x_{i}^{2};\cdots ;x_{i}^{d};y)(xi1;xi2;⋯;xid;y)
用向量 xi{\textbf{x}}_{i}xi表示 ddd 个属性，w{\textbf{w}}w表示参数的向量;
xi=(xi1;xi2;⋯ ;xid){\textbf{x}}_{i}=(x_{i}^{1};x_{i}^{2};\cdots ;x_{i}^{d})xi=(xi1;xi2;⋯;xid)
w=(w1;w2;⋯ ;wd){\textbf{w}}=(w_{1};w_{2};\cdots ;w_{d})w=(w1;w2;⋯;wd)
则回归模型可以表示成；
f(xi)=wTxi+bf({\textbf{x}}_{i})={\textbf{w}}^{T}{\textbf{x}}_{i}+bf(xi)=wTxi+b
为了便于讨论，令
w^=(w;b)\hat{\textbf{w}}=({\textbf{w}};b)w^=(w;b)
对应的将x{\textbf{x}}x也增加一项；
xi=(xi1;xi2;⋯ ;xid;1){\textbf{x}}_{i}=(x_{i}^{1};x_{i}^{2};\cdots ;x_{i}^{d};1)xi=(xi1;xi2;⋯;xid;1)
f(xi)=w^Txif({\textbf{x}}_{i})=\hat{\textbf{w}}^{T}{\textbf{x}}_{i}f(xi)=w^Txi
X=(x11x12⋯x1d1x21x22⋯x2d1⋮⋱⋱⋮⋮xn1xn2⋯xnd1)=(x1T1x2T1⋮⋮xnT1){\textbf{X}}=\begin{pmatrix} x_{1}^{1} & x_{1}^{2} & \cdots & x_{1}^{d} &1 \\ x_{2}^{1} & x_{2}^{2} & \cdots & x_{2}^{d} &1\\ \vdots & \ddots & \ddots &\vdots & \vdots \\ x_{n}^{1} & x_{n}^{2} & \cdots & x_{n}^{d} &1 \end{pmatrix}=\begin{pmatrix} {\textbf{x}}_{1}^{T} &1 \\ {\textbf{x}}_{2}^{T} &1 \\ \vdots &\vdots \\ {\textbf{x}}_{n}^{T} & 1 \end{pmatrix}X=⎝⎜⎜⎜⎛x11x21⋮xn1x12x22⋱xn2⋯⋯⋱⋯x1dx2d⋮xnd11⋮1⎠⎟⎟⎟⎞=⎝⎜⎜⎜⎛x1Tx2T⋮xnT11⋮1⎠⎟⎟⎟⎞
同样地； y=(y1;y2;⋯ ;yn){\textbf{y}}=({y}_{1};{y}_{2};\cdots;{y}_{n})y=(y1;y2;⋯;yn)，那么与上面的类似，有;
w^∗=arg minw^(y−Xw^)T(y−Xw^)\hat{\textbf{w}}^{*}=\underset{\hat{\textbf{w}}}{arg\, min}({\textbf{y}}-{\textbf{X}}\hat{\textbf{w}})^{T}({\textbf{y}}-{\textbf{X}}\hat{\textbf{w}})w^∗=w^argmin(y−Xw^)T(y−Xw^)
令E=(y−Xw^)T(y−Xw^)E=({\textbf{y}}-{\textbf{X}}\hat{\textbf{w}})^{T}({\textbf{y}}-{\textbf{X}}\hat{\textbf{w}})E=(y−Xw^)T(y−Xw^)；对w^\hat{\textbf{w}}w^求偏导数；
∂E∂w^=2XT(Xw^−y)\frac{\partial E}{\partial \hat{\textbf{w}}}=2{\textbf{X}}^{T}({\textbf{X}}\hat{\textbf{w}}-\textbf{y})∂w^∂E=2XT(Xw^−y)
最终所求
w^∗=(XTX)−1XTy\hat{\textbf{w}}^{*}=({\textbf{X}}^{T}{\textbf{X}})^{-1}{\textbf{X}}^{T}\textbf{y}w^∗=(XTX)−1XTy

继续推广延伸

更一般的，可以将这种线性回归用于非线性函数的拟合；例如对于对数函数；
f(xi)=y=ln(wTxi+b)f({\textbf{x}}_{i})=y=ln({\textbf{w}}^{T}{\textbf{x}}_{i}+b)f(xi)=y=ln(wTxi+b)
可以利用反函数，将其转换成h=ey=wTxi+bh=e^{y}={\textbf{w}}^{T}{\textbf{x}}_{i}+bh=ey=wTxi+b；成为一个线性回归的问题，最终再转换成对数形式；

甚至可以用最小二乘法拟合多项式曲线。
y=a0+a1x+a2x2+⋯+amxmy=a_{0}+a_{1}x+a_{2}x^{2}+\cdots +a_{m}x^{m}y=a0+a1x+a2x2+⋯+amxm
X=(x1mx1m−1⋯x111x2mx2m−1⋯x211⋮⋱⋱⋮⋮xnmxnm−1⋯xn11){\textbf{X}}=\begin{pmatrix} x_{1}^{m} & x_{1}^{m-1} & \cdots & x_{1}^{1} &1 \\ x_{2}^{m} & x_{2}^{m-1} & \cdots & x_{2}^{1} &1\\ \vdots & \ddots & \ddots &\vdots & \vdots \\ x_{n}^{m} & x_{n}^{m-1} & \cdots & x_{n}^{1} &1 \end{pmatrix}X=⎝⎜⎜⎜⎛x1mx2m⋮xnmx1m−1x2m−1⋱xnm−1⋯⋯⋱⋯x11x21⋮xn111⋮1⎠⎟⎟⎟⎞
此时上标mmm不再表示第mmm个属性，而是表示mmm次幂；在修改过X{\textbf{X}}X后，最终w^∗\hat{\textbf{w}}^{*}w^∗依然表示为；
w^∗=(XTX)−1XTy\hat{\textbf{w}}^{*}=({\textbf{X}}^{T}{\textbf{X}})^{-1}{\textbf{X}}^{T}\textbf{y}w^∗=(XTX)−1XTy

Matlab 实现高次多项式拟合

用Matlab实现了最小二乘法一元高次多项式拟合；代码如下；

**%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%不用先生，2018.11.22，最小二乘法一元高次多项式拟合
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
clear all;
xi=[-1,0,2,3,4,5,8,9,10,15,20];
yi=[-2,-0.5,1.5,2.8,2.8,3.9,7.6,7.5,8,13,16.5];[~,n]=size(xi);
m=8;     %%最高8次多项式；
W=zeros(m+1,m);for i=1:mX=zeros(n,i+1);    %%创建一个n*+1的矩阵；for j=1:nfor s=1:i          X(j,s)=xi(j)^(i+1-s);   %%计算矩阵中每一个元素的值endX(j,i+1)=1;endXX=inv(X'*X)*X'*yi';  %%得到系数并赋值for j=1:i+1W(j,i)=XX(j);end
end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%构造方程，画图
x=-5:0.01:25;
for i=1:my=W(1,i)*x.^i;for j=2:i+1y=y+W(j,i)*x.^(i-j+1);end   subplot(2,4,i);plot(xi,yi,'*');hold on;plot(x,y);hold off
end
suptitle('用1到8次多项式拟合')
**

很明显，在高次多项式拟合时，曲线边缘出现了龙格现象。

参考

《机器学习》，周志华著，清华大学出版社；

已完。。