最小二乘法

最小二乘法原理与编程实现

背景

数据与数据（变量与变量）之间，很多时候是存在一些关系的，如线性关系和非线性关系。我们常常会希望找到数据之间的关系，用一个函数（或者一条曲线）去描述两个变量之间的关系。
然而因为各种原因，测量得到的数据是会存在误差的，于是我们要一种方法去减少误差带来的干扰，尽可能的描绘出数据之间的关系，这个方法就是最小二乘法，也称最小方差法。

原理

比如说，我们得到了四对数据，分别是（1，0.9），（2, 2.3），（3.5, 2.6），（4, 3.8）把它们画在平面直角坐标系上

从图像上我们很容易发现，这些点大致分布在一条直线上面，于是我们大胆的猜测y与x呈线性关系，于是我们很自然地想要用一条直线去拟合他们。
也就是

[y = kx + b ]

然而，事实是它们只是看起来像一条直线，但实际上并不是一条直线。把方程组列出来，

[egin{cases} 1k + b = 0.9 \ 2k + b = 2.3 \ 3.5k + b = 2.6 \ 4k + b = 3.8 \ end{cases} ]

求解一下就会发现，这个方程组无解！从图像上来看就是我们找不到一条直线通过所有的点。
那该怎么办？那我们只好退一步，不要求找一条完全经过所有点的直线，只要求找一条能够大致刻画它们关系的直线，并且使得误差最小。那误差怎么衡量呢？
我们随便画一条直线，比如说

[y = kx + b\(图例k=1,b=0) ]

因为有些点没有落到直线上，于是我们把误差定义为每一个观测点的y值和我们预测的真实值之间的距离的平方,也就是

[e_i = (f(x_i)-y_i)^2=(kx_i + b - y_i)^2 ]

我们的目标是使得总体的误差最小，也就是

[min(sum_{i=1}^n{e_i}) ]

这是一条关于k和b的二元函数，我们求偏导数并找到导数为0的点就可以使其最小，也即是令

[egin{cases}frac{partial{e}}{partial{k}}=2sum{(kx_i + b-y_i)x_i} = 0 \frac{partial{e}}{partial{b}}=2sum{(kx_i + b-y_i)} = 0 \end{cases} ]

我们把原来的方程组写成矩阵的形式(这里把k,b当成待求参数，写成(x_1,x_2))

[egin{bmatrix} 1 & 1 \ 2 & 1 \ 3.5 & 1 \ 4 & 1 end{bmatrix} egin{bmatrix} x_1 \ x_2 end{bmatrix} =egin{bmatrix} 0.6\3\2.6\3.8 end{bmatrix}]

然后我们把误差写成列向量，也就是

[egin{bmatrix} e_1\e_2\e_3\e_4 end{bmatrix} = egin{bmatrix} 1 & 1 \ 2 & 1 \ 3.5 & 1 \ 4 & 1 end{bmatrix} egin{bmatrix} x_1 \ x_2 end{bmatrix} -egin{bmatrix} 0.6\3\2.6\3.8 end{bmatrix} ]

记

[E =egin{bmatrix}e_1\e_2\e_3\e_4end{bmatrix}, A =egin{bmatrix}1 & 1 \2 & 1 \3.5 & 1 \4 & 1 end{bmatrix}, X = egin{bmatrix}x_1 \ x_2end{bmatrix}, b = egin{bmatrix}0.6\3\2.6\3.8end{bmatrix} ]

则得到

[E = AX - b ]

要让(E^2)最小，即

[E^TE = (AX-b)^T(AX-b) ]

我们对X求导

[frac{partial{E^2}}{partial{X}}= frac{partial{(AX-b)^T(AX-b)}}{partial{X}}\= frac{partial{(AX-b)^T}}{partial{X}} *(AX-b)+ frac{partial{(AX-b)^T}}{partial{X}}*(AX-b)\= 2frac{partial{(AX-b)^T}}{partial{X}}*(AX-b)\ =2frac{partial{(AX)}^T}{partial{X}}*(AX-b)\ =2A^T(AX-b) ]

我们令其等于0，就可以解出X

[X = (A^TA)^{-1}A^Tb ]

这是经过计算分析得到的，那有没有一些更加直觉化的解析呢？

有！我们可以通过向量的角度去理解它的原理！

刚刚那个矩阵方程可以写成向量方程的形式，

[egin{bmatrix}1 & 1 \2 & 1 \3.5 & 1 \4 & 1 end{bmatrix}egin{bmatrix}x_1 \ x_2end{bmatrix}=x_1egin{bmatrix}1 \2 \3.5 \4 \end{bmatrix}+x_2egin{bmatrix}1 \1 \1 \1 end{bmatrix}=egin{bmatrix}0.6\3\2.6\3.8end{bmatrix} ]

我们记

[v_1 = egin{bmatrix} 1 \ 2 \ 3.5 \ 4 \ end{bmatrix}, v_2 = egin{bmatrix} 1 \ 1 \ 1 \ 1 end{bmatrix} b = egin{bmatrix} 0.6\3\2.6\3.8 end{bmatrix} ]

容易看到，(v_1)和(v_2)是线性无关的，也就是它们的线性组合可以张成四维空间里的一个平面。

我们想要找的是它的某种组合使得其等于b，但这种组合找不到，于是我们就去找那个平面里的一条线，使得它到b的距离最短。容易想象那条线就是b在span((v_1,v_2))下的投影，我们记该投影向量为c ，则b和c之间的距离就是(||b-c||),我们记为e，这时e最小，且满足e垂直与span((v_1,v_2))。所以e肯定也垂直于(v_1,v_2)，即是

[ev_1=0,ev_2=0\ev_1+ev_2=0\egin{bmatrix}1\2\3.5\4end{bmatrix}*e+egin{bmatrix}1\1\1\1end{bmatrix}*e=0\ ]

这里是点乘，写成矩阵相乘就是

[egin{bmatrix}1&2&3.5&4\1&1&1&1end{bmatrix}e=0\ ]

注意到这里是A的转置(A^T),再把e=(AX-b)代入得

[A^T(AX-b)=0 ]

结果是和上面推导的一样的

以上的方法对于二次，三次拟合都是成立的，具体可以参考其他资料。

拟合的结果

一次拟合

###二次拟合 ###三次拟合 ##代码(Python)

# -*- coding: utf-8 -*-
"""
Created on Mon Feb 17 13:51:43 2020
最小二乘法
@author: urahyou
"""
import numpy as np
import matplotlib.pyplot as plt

x = np.array([3.0, 5, 6, 8, 10])
y = np.array([5.0, 2, 1, 2, 4])

p1 = plt.scatter(x,y,c='red')

def LSD(x, y, n):
    N = x.size  #获取方程组个数
    A = np.ones(N)
    
    for i in range(1, n+1):
        A = np.vstack((A,x**i))  #垂直拼接
    A = A.T  #转置回来
    #求解
    B = np.linalg.inv(A.T@A)@A.T
    #求出解系数
    sol = np.dot(B, y)
    return sol

def poly(x,sol):
    y =  np.zeros_like(x)  #每一个x,对应一个y
    n = sol.size
    for i in range(n):
        y += sol[i]*x**i
    return y


sol = LSD(x,y,2)
X = np.arange(0, 14, 0.1)
Y = poly(X,sol)
p2 = plt.plot(X,Y)