线性最小二乘法

UTF8gbsn

本文介绍线性最小二乘法.
线性最小二乘法比较简单并且在很多场景里面都有应用.
下面我们就来举几个例子.

例子

已知散点集合 X = { x 1 , x 2 , ⋯   , x n } , Y = { y 1 , y 2 , ⋯   , y n } X=\{x_1,x_2,\cdots, x_n\}, Y=\{y_1,y_2,\cdots, y_n\} X={x1​,x2​,⋯,xn​},Y={y1​,y2​,⋯,yn​},
求参数方程$y=ax+b$.

已知是$X,Y$未知是$a,b$. 这个问题是一个标准的可以使用最小二乘方法的问题.
我们希望残差的平方和最小. 也就是$y=ax+b$更好的贴近数据集$X,Y$.

$$r_i=y_i-(a{x}_i+b)$$

我们的代价函数定义为 $$F(a,b)=\sum _{i=1}^n[y_i-(a{x}_i+b){]}^2$$

于是最小二乘问题变化为

$$(a^{\ast },b^{\ast })=\arg \underset{a,b}{\min }F(a,b)=\arg \underset{a,b}{\min }\sum _{i=1}^n[y_i-(a{x}_i+b){]}^2$$

我们把上式改写成向量的形式.

$$v_i=(x_i,1{)}^T,g=(a,b{)}^T,Y=(y_1,y_2,\cdots ,y_n{)}^T,V=(v_1,v_2,\cdots ,v_n)$$

那么,向量的形式来描述上线的方程为
$$g^{\ast }=\arg \underset{g}{\min }F(g)=\arg \underset{g}{\min }(Y-V^{T}g{)}^T(Y-V^{T}g)$$

那么对上式求导数并等于0可得

$$V{V}^{T}g=V^{T}Y\Rightarrow g=(V{V}^T{)}^{+}{V}^{T}Y$$

其中符号$(V{V}^T{)}^{+}$,表示$V{V}^T$的伪逆.

python实现代码如下.

import numpy as np
from matplotlib import pyplot as plt

N = 100
a = 0
b = 10
s = (b-a)/N
noise = 2

def linear_function(x):
    return 3.1*x+2.5


def plt_scatter(points):
    x = [cell[0] for cell in points]
    y = [cell[1] for cell in points]
    plt.scatter(x,y, s=5)

def linear_least_square(X, Y):
    M=X.dot(X.T)
    invM = np.linalg.inv(M)
    FY = X.dot(Y)
    return invM.dot(FY)

points = []
for k in range(N):
    r = 1
    if np.random.randint(2) == 0:
        r = -1
    x = a + k*s + np.random.rand()*noise * r
    r = 1
    if np.random.randint(2) == 0:
        r = -1
    y = linear_function(x) + np.random.rand()*noise *r
    points.append((x,y))

plt_scatter(points)

X = [(cell[0], 1) for cell in points]
X = np.array(X)
X = X.T
Y = [cell[1] for cell in points]
Y = np.array(Y)
Y = Y.T

g=linear_least_square(X,Y)
x1, x2 = g[0], g[1]


x_arr = []
y_arr = []

for k in range(N):
    x = a + k*s
    y = x*x1+x2
    x_arr.append(x)
    y_arr.append(y)

plt.plot(x_arr,y_arr, '-r')
plt.show()