MATLAB | 最小二乘法的两种解读

最小二乘法

大部分的最小二乘法公式推导，都是使用的 代价函数偏导 的方式来求得的，在这里首先展示如何通过代价函数求偏导的方式得到最小二乘公式，再展示李扬老师讲解的如何由向量到子空间的距离得来最小二乘法公式。

代价函数与最小二乘法

假设我们的拟合结果为：

h θ ( x ) = ∑ i = 0 n θ i x i = θ 0 + θ 1 x 1 + ⋯ + θ n x n h_{\theta}(x)=\sum_{i=0}^{n} \theta_{i} x_{i}=\theta_{0}+\theta_{1} x_{1}+\cdots+\theta_{n} x_{n} hθ(x)=i=0∑nθixi=θ0+θ1x1+⋯+θnxn

则平方损失函数为：

J ( θ ) = ∑ j = 1 m ( h θ ( x ( j ) ) − y ( j ) ) 2 = ∑ j = 1 m ( ∑ i = 0 n θ i x i ( j ) − y ( j ) ) 2 J(\theta)=\sum_{j=1}^{m}\left(h_{\theta}\left(x^{(j)}\right)-y^{(j)}\right)^{2}=\sum_{j=1}^{m}\left(\sum_{i=0}^{n} \theta_{i} x_{i}^{(j)}-y^{(j)}\right)^{2} J(θ)=j=1∑m(hθ(x(j))−y(j))2=j=1∑m(i=0∑nθixi(j)−y(j))2

其中 x ( j ) , y ( j ) x^{(j)},y^{(j)} x(j),y(j)为一组数据点, x ( j ) x^{(j)} x(j)是长度为 n + 1 n+1 n+1的向量， y ( j ) y^{(j)} y(j)是数值，而 x i → \overrightarrow{x_i} xi 我们定义为长为 m m m的向量，这在之后会被用到，注意这里为了追求整齐性，我们直接定义 x 0 x_0 x0衡为常数1：
x ( j ) = ( 1 , x 1 ( j ) , … , x n ( j ) ) x i → = ( x i ( 1 ) , x i ( 2 ) , … , x i ( m ) ) \begin{aligned} &x^{(j)}=\left(1, x_{1}^{(j)}, \ldots, x_{n}^{(j)} \right)\\ &\overrightarrow{x_{i}}=\left(x_{i}^{(1)}, x_{i}^{(2)}, \ldots, x_{i}^{(m)}\right) \end{aligned} x(j)=(1,x1(j),…,xn(j))xi =(xi(1),xi(2),…,xi(m))

平方损失函数的形式只有极小值，没有极大值，我们要使代价函数最小，我们要找到其极值点，即偏导均为0的点，代价函数对于各参数偏导如下：

∂ J ( θ ) θ i = 2 ∑ j = 1 m ( h θ ( x ( j ) ) − y ( j ) ) x i ( j ) , i = 0 , 1 , 2 , . . . n \frac{\partial J(\theta)}{\theta_{i}}=2 \sum_{j=1}^{m}\left(h_{\theta}\left(x^{(j)}\right)-y^{(j)}\right) x_i^{(j)},i=0,1,2,...n θi∂J(θ)=2j=1∑m(hθ(x(j))−y(j))xi(j),i=0,1,2,...n
令偏导为0得：
∑ j = 1 m ( h θ ( x ( j ) ) − y ( j ) ) x i ( j ) = 0 , i = 0 , 1 , 2 , . . . , n ⇒ ∑ j = 1 m ( h θ ( x ( j ) ) x i ( j ) ) = ∑ j = 1 m ( y ( j ) x i ( j ) ) , i = 0 , 1 , 2 , . . . , n \begin{aligned} &\sum_{j=1}^{m}\left(h_{\theta}\left(x^{(j)}\right)-y^{(j)}\right) x_i^{(j)}=0,i=0,1,2,...,n \\ \Rightarrow &\sum_{j=1}^{m}\left(h_{\theta}\left(x^{(j)}\right)x_i^{(j)}\right)=\sum_{j=1}^{m}\left(y^{(j)}x_i^{(j)}\right),i=0,1,2,...,n \end{aligned} ⇒j=1∑m(hθ(x(j))−y(j))xi(j)=0,i=0,1,2,...,nj=1∑m(hθ(x(j))xi(j))=j=1∑m(y(j)xi(j)),i=0,1,2,...,n

实际上若是令:

X = [ − ( x ( 1 ) ) − − ( x ( 2 ) ) − ⋮ − ( x ( m ) ) − ] Y = [ y ( 1 ) y ( 2 ) ⋮ y ( m ) ] θ = [ θ 0 θ 1 ⋮ θ n ] X=\left[\begin{array}{c} -\left(x^{(1)}\right)- \\ -\left(x^{(2)}\right)- \\ \vdots \\ -\left(x^{(m)}\right)- \end{array}\right] \quad Y=\left[\begin{array}{c} y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(m)} \end{array}\right] \quad \theta=\left[\begin{array}{c} \theta_0 \\ \theta_1\\ \vdots \\ \theta_n \end{array}\right] X=⎣⎢⎢⎢⎡−(x(1))−−(x(2))−⋮−(x(m))−⎦⎥⎥⎥⎤Y=⎣⎢⎢⎢⎡y(1)y(2)⋮y(m)⎦⎥⎥⎥⎤θ=⎣⎢⎢⎢⎡θ0θ1⋮θn⎦⎥⎥⎥⎤
为了尽量写的清晰明了，在这里对比一下 X X X与 X ′ X^{\prime} X′：
X = [ − ( x ( 1 ) ) − − ( x ( 2 ) ) − ⋮ − ( x ( m ) ) − ] X ′ = [ − ( x 0 → ) − − ( x 1 → ) − ⋮ − ( x n → ) − ] X=\left[\begin{array}{c} -\left(x^{(1)}\right)- \\ -\left(x^{(2)}\right)- \\ \vdots \\ -\left(x^{(m)}\right)- \end{array}\right] \quad X^{\prime}=\left[\begin{array}{c} -\left(\overrightarrow{x_0}\right)- \\ -\left(\overrightarrow{x_1}\right)- \\ \vdots \\ -\left(\overrightarrow{x_n}\right)- \end{array}\right] X=⎣⎢⎢⎢⎡−(x(1))−−(x(2))−⋮−(x(m))−⎦⎥⎥⎥⎤X′=⎣⎢⎢⎢⎢⎢⎢⎡−(x0 )−−(x1 )−⋮−(xn )−⎦⎥⎥⎥⎥⎥⎥⎤
则有：
∑ j = 1 m ( h θ ( x ( j ) ) x i ( j ) ) = ∑ j = 1 m ( y ( j ) x i ( j ) ) , i = 0 , 1 , 2 , . . . , n ⇒ x i → X θ = x i → Y , i = 0 , 1 , 2 , . . . , n ⇒ X ′ X θ = X ′ Y ⇒ θ = ( X ′ X ) − 1 X ′ Y \begin{aligned} \sum_{j=1}^{m}\left(h_{\theta}\left(x^{(j)}\right)x_i^{(j)}\right)&=\sum_{j=1}^{m}\left(y^{(j)}x_i^{(j)}\right),i=0,1,2,...,n \\ \Rightarrow \overrightarrow{x_i}X\theta&=\overrightarrow{x_i}Y,i=0,1,2,...,n \\ \Rightarrow X^{\prime}X\theta&=X^{\prime}Y \\ \Rightarrow \theta&=(X^{\prime}X)^{-1}X^{\prime}Y \end{aligned} j=1∑m(hθ(x(j))xi(j))⇒xi Xθ⇒X′Xθ⇒θ=j=1∑m(y(j)xi(j)),i=0,1,2,...,n=xi Y,i=0,1,2,...,n=X′Y=(X′X)−1X′Y

向量到子空间的距离与最小二乘法

注意，为了迎合(线代/高代)的习惯，以下内容推导过程依旧使用常用的 A X = b AX=b AX=b符号，但在最后使用时，为了对多项式函数或多元函数各变量进行对应，将其修改为 X θ = Y X\theta=Y Xθ=Y的符号体系，请大家注意对应。

实数域上的线性方程组 A s × n X n × 1 = b s × 1 A_{s \times n}X_{n\times1}=b_{s\times1} As×nXn×1=bs×1大部分情况下 s > n s>n s>n无解，现在想找一个 X 0 X_0 X0使得 ∣ A X 0 − b ∣ \left|A X_{0}-b\right| ∣AX0−b∣最小，这就是所谓的最小二乘法。

∣ A X 0 − b ∣ 最小 ⟺ 对于任意的 X 都有 ∣ A X 0 − b ∣ ≤ ∣ A X − b ∣ ⟺ A X 0 + b ⊥ U 其中 U = { A X ∣ X ∈ R n } = L ( α 1 , ⋯ , α n ) ⟺ A X 0 − b ⊥ α i ( i = 1 , 2 , ⋯ , n ) ⟺ ( α 1 ′ ⋮ α n ′ ) ( A X 0 − b ) = 0 ( i = 1 , 2 , ⋯ , n ) ⟺ A ′ ( A X 0 − b ) = 0 ⟺ A ′ A X 0 = A ′ b . \begin{aligned} \left|A X_{0}-b\right| \text { 最小 } & \Longleftrightarrow \text { 对于任意的 } X \text { 都有 }\left|A X_{0}-b\right| \leq|A X-b| \\ & \Longleftrightarrow A X_{0}+b \perp U \text { 其中 } U=\left\{A X \mid X \in \mathbb{R}^{n}\right\}=L\left(\alpha_{1}, \cdots, \alpha_{n}\right) \\ & \Longleftrightarrow A X_{0}-b \perp \alpha_{i}(i=1,2, \cdots, n) \\ & \Longleftrightarrow\left(\begin{array}{c} \alpha_{1}^{\prime} \\ \vdots \\ \alpha_{n}^{\prime} \end{array}\right)\left(A X_{0}-b\right)=0(i=1,2, \cdots, n) \\ & \Longleftrightarrow A^{\prime}\left(A X_{0}-b\right)=0 \\ & \Longleftrightarrow A^{\prime} A X_{0}=A^{\prime} b . \end{aligned} ∣AX0−b∣ 最小 ⟺ 对于任意的 X 都有 ∣AX0−b∣≤∣AX−b∣⟺AX0+b⊥U 其中 U={AX∣X∈Rn}=L(α1,⋯,αn)⟺AX0−b⊥αi(i=1,2,⋯,n)⟺⎝⎜⎛α1′⋮αn′⎠⎟⎞(AX0−b)=0(i=1,2,⋯,n)⟺A′(AX0−b)=0⟺A′AX0=A′b.

若最小的 X 0 X_0 X0存在，则一定满足 A ′ A X 0 = A ′ b A^{\prime} A X_{0}=A^{\prime} b A′AX0=A′b形式，现证明该方程组一定有解，原因是：

r ( A ′ A , A ′ b ) = r ( A ′ ( A , b ) ) ≤ r ( A ′ ) = r ( A ) . r\left(A^{\prime} A, A^{\prime} b\right)=r\left(A^{\prime}(A, b)\right) \leq r\left(A^{\prime}\right)=r(A) . r(A′A,A′b)=r(A′(A,b))≤r(A′)=r(A).

同时 r ( A ′ A , A ′ b ) ≥ r ( A ′ A ) = r ( A ) r\left(A^{\prime} A, A^{\prime} b\right) \geq r\left(A^{\prime} A\right)=r(A) r(A′A,A′b)≥r(A′A)=r(A),这就说明：

r ( A ′ A , A ′ b ) = r ( A ) = r ( A ′ A ) . r\left(A^{\prime} A, A^{\prime} b\right)=r(A)=r\left(A^{\prime} A\right) . r(A′A,A′b)=r(A)=r(A′A).

因此要求解让 ∣ X θ − Y ∣ \left|X \theta-Y\right| ∣Xθ−Y∣取最小值的 θ \theta θ 只需求解 ( X ′ X ) − 1 X ′ Y (X^{\prime}X)^{-1}X^{\prime}Y (X′X)−1X′Y 。

最小二乘法与多项式拟合

以下展示自己编写最小二乘法拟合多项式与MATLAB自带函数 polyfit 拟合多项式的参数对比，注意，为了和MATLAB自带函数保持一致， θ \theta θ 向量变为第一个参数为 θ n \theta_n θn ,最后一个参数为 θ 0 \theta_0 θ0 , X X X 矩阵也做了相应的调整：

% 最小二乘法多项式拟合% 原三次函数+随机噪声
f=@(x)x.^3+6.*x.^2-2.*x+4+(rand(size(x))-.5).*20;% 构造原始数据
x=-5:.1:5;
y=f(x);% 自己写一个最小二乘
n=3;% 最高次数为三次
X=(x').^(n:-1:0);
theta1=((X'*X)\X'*y')';% MATLAB自带多项式拟合
theta2=polyfit(x,y,n);% 输出对比
disp(theta1)
disp(theta2)% 一个小技巧，下面的写法能够快速将
% 参数向量变成有关x的多项式匿名函数
func=matlabFunction(poly2sym(theta1));

theta1=
0.9686 6.0178 -1.8845 4.3362
theta2=
0.9686 6.0178 -1.8845 4.3362

多项式拟合结果绘图：

% 绘图部分% 保持坐标区域不刷新并添加网格
ax=gca;hold(ax,'on');grid(ax,'on');% 绘制原数据点和拟合结果
plot(x,y,'o','MarkerFaceColor',[94,142,179]./255);
plot(x,func(x),'Color',[0,64,115]./255,'LineWidth',2);% 修饰一下
ax.FontName='cambria';
ax.LineWidth=1.5;
ax.GridLineStyle='--';
ax.XColor=[1,1,1].*.3;
ax.YColor=[1,1,1].*.3;
ax.ZColor=[1,1,1].*.3;

最小二乘法与多元线性回归

以下展示自己编写最小二乘法进行多元线性回归与MATLAB自带函数 regress 进行多元线性回归的参数对比：

% 最小二乘法多元线性回归% 原二元函数+随机噪声
f=@(x1,x2) 3.*x1+4.*x2+5+(rand(size(x1))-.5).*10;% 构造原始数据
[x1,x2]=meshgrid(-5:.5:5,-5:.5:5);
y=f(x1,x2);% 自己写一个最小二乘
X=[x1(:),x2(:),ones(size(x1(:)))];
theta1=((X'*X)\X'*y(:));% MATLAB多元线性回归
theta2=regress(y(:),X);% 输出对比
disp(theta1)
disp(theta2)% 构造拟合结果的二元匿名函数
func=matlabFunction([sym('x1'),sym('x2'),1]*theta1);

theta1=
2.9285 4.0688 4.7520
theta2=
2.9285 4.0688 4.7520

多元线性回归结果绘图：

% 绘图部分% 保持坐标区域不刷新并添加网格
ax=gca;hold(ax,'on');grid(ax,'on');% 绘制原数据点和拟合结果
mesh(x1,x2,func(x1,x2),'FaceColor','flat','FaceAlpha',.8)
scatter3(x1(:),x2(:),y(:),20,'filled')% 修饰一下
ax.FontName='cambria';
ax.LineWidth=1.5;
ax.GridLineStyle='--';
ax.XColor=[1,1,1].*.3;
ax.YColor=[1,1,1].*.3;
ax.ZColor=[1,1,1].*.3;
view(30,20)