chapter2 最小二乘法完整实践理论

Posted on | In | | Visitors: ℃
Words count in article: 3k | Reading time ≈ 13

chapter2- 最小二乘法完整实践理论

1.概要

​ 此章，通过最小二乘法，在统计学习，机器学习，深度学习的轨道上进行讨论，它们叁之间的共性与异性。

2.什么是最小二乘法

​ 最小二乘法（least squares method），又称最小平方法，是一种数学优化方法。它通过最小化误差的平方和寻找数据的最佳函数匹配。

​ 利用最小二乘法可以简便地求得未知的数据，并使得求得的数据与实际数据之间误差的平方和为最小。

​ 高斯于1823年在误差$e_1,…,e_n$独立同分布的假定下,证明了最小二乘方法的一个最优性质: 在所有无偏的线性估计类中,最小二乘方法是其中方差最小的。 对于数据$(x_i,y_i)(i = 1,2,3,…,m)$,能够拟合出下列公式（1），即是一个多项式函数
$$
h_\theta(x_1,x_2,…x_n)=\theta_0+\theta_1x_1+…+\theta_{n−1}x_{n−1}
$$

​ 利用数据$(x_i,y_i)$代入公式1,利用数学知识是能够解出最佳拟合参数$W_i$。即能够得到$H(x)$使得与真实值相差最小—残差$e_i$。因此，得到一个代价函数(Cost Function)，亦为损失函数(公式2)。

$$
min\sum_{i=1}^{m}e_i = \min\sum_{i=1}^n(y_i-h_θ(x^{(i)})^2
$$

3.统计学中如何去解最小二乘法

​ 我们以一次线性函数模型来进行用数学公式求解，让大家初始一下统计学如何去拟合数据的。对于高阶函数，统计学知识求解，感兴趣的同学，可以去Google一下。

3.1 代数法求解最小二乘法

​ 最简单的线性式$h_\theta= \theta_1x + \theta_0$有代价函数为：
$$
J(\theta_1,\theta_0) =\sum_{i=1}^n(y_i-h_θ(x^{(i)})^2=\sum_{i=1}^n(y_i-\theta_1x - \theta_0)^2
$$
​ 有了公式(3)代价函数，我们如何去使它最小了，最常用的解法就是对$\theta_1$和$\theta_0$分别求编导，令编导为 零，就可以解出$\theta_1$和$\theta_0$的值，即过程为：

​ 　$J(\theta_1,\theta_0)$对$\theta_0 $求导，得到如下方程：
$$
\sum_{i=1}^m(y^{i}-\theta_0-\theta_1x^{(i)}) = 0
$$
​ $J(\theta_1,\theta_0)$对$\theta_1 $求导，得到如下方程：
$$
\sum_{i=1}^m(y^{i}-\theta_0-\theta_1x^{(i)})x^{(i)} = 0
$$
​ 公式(4)、(5)组成一个二元方程组，容易解出$\theta_1$和$\theta_0$的值

​
$$
\theta_0=\frac{\sum_{i=1}^m(x^{(i)})^2\sum_{i=1}^my^{x^{(i)}}-\sum_{i=1}^mx^{(i)}\sum_{i=1}^mx^{(i)}y^{(i)}}
{m\sum_{i=1}^m(x^{(i)})^2 - (\sum_{i=1}^mx^{(i)})^2}
$$

$$
\theta_1=\frac{m\sum_{i=1}^mx^{(i)}y^{(i)} - \sum_{i=1}^mx^{(i)}\sum_{i=1}^my^{(i)}}
{m\sum_{i=1}^m(x^{(i)})^2 - (\sum_{i=1}^mx^{(i)})^2}
$$

​ 利用求导的方式,很容推广到多个样本特征的多项式拟合当中。即在拟合函数公式(1)当中，求解参数 $\theta_i(i=0,1,2,3,\cdots,m)$,可以简化拟合函数为：
$$
h_θ(x_1,x_2,…x_n)=\sum_{i=0}^{n}\theta_ix_i
$$
​ 损失函数就可以表示为：
$$
J(\theta_0,\theta_1,\cdots,\theta_j) =\sum_{i=1}^m(h_θ(x_0^{(i)},x_1^{(i)},\cdots,x_j^{(i)}) - y^{(j)})^2=\sum_{i=1}^n(\sum_{j=1}^m\theta_jx_j^{(i)}-y^{(i)})^2
$$
​ 利用损失函数分别对$\theta_i(i=0,1,2,3,\cdots,m)$求导，并令导数为0可得：
$$
\sum_{i=1}^n(\sum_{j=1}^m\theta_jx_j^{(i)}-y^{(i)})x_j^=0 \ (j=0,\cdots,m)
$$
​ 就可以得到一个M+1个一元方程组，求解这个方程组，就可以得到所有的$\theta_i(i=0,1,2,3,\cdots,m)$。代 入参数到$h_\theta$就可以得到拟合函数了

3.2 矩阵法求解最小二乘法

​ 对于上述$h_\theta$,矩阵表达式形式为：

$$
h_\theta(X) = \theta X
$$
​ 即$h_{\theta}=y= \theta_1 x + \theta_0$,写成代价矩阵式为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
$$
\min\limits_{x_0,x_1}
\begin{Vmatrix} 
    \begin{pmatrix} 
          \theta_0\\
          \theta_1
    \end{pmatrix}
    \begin{pmatrix} 
           1&x_1 \\
           .&. \\
           .&. \\
           .&. \\
           1&x_m \\
    \end{pmatrix}
    -
    \begin{pmatrix} 
           y_1\\
           .\\
           .\\
           .\\
           y_m \\
	\end{pmatrix}
\end{Vmatrix}_2 =\  
\min\limits_{x}||\theta X-Y||_2
$$

​ 化简为矩阵损失函数 $J(θ)=\frac{1}{2}(θX−Y)^T(θX−Y)$，其中， 假设函数$h_{\theta}(X)$为mx1的向量,$\theta$为nx1的 向量，里面有n个代数法的模型参数。$X$为$m\times n$维的矩阵。m代表样本的个数，n代表样本的特征数.

​ 对于损失函数求导，并且令为0，可能

​
$$
\frac{\partial}{\partial \theta}J{(\theta)} = X^T(X\theta-Y) = 0
$$
​ 这里用到了矩阵求导链式法则，和两个矩阵求导的公式：

​ 公式1: $\frac{\partial}{\partial x}(X^TX) = 2X \ x为向量$

​ 公式2: $ \Delta_X f(AX+B)=A^T\Delta_Y f,Y=AX+B,\ f(Y)为标量$

​ 对上述求导等式整理后可得：
$$
X^TX\theta = X^TY
$$
​ 两边同时左乘$(X^TX)$可得：
$$
\theta = (X^TX)^{-1}X^TY
$$
​ 这样就可以求出$\theta$向量表达式的公式，免去了代数法一个个去求导的麻烦.只要给了数据，就可以用公式(15)算出$\theta$,得到数据拟合函数$h_{\theta}(X)$了。

4.机器学习实现最小二乘法

现在我们使用机器学习方法去实现最小二乘法。这里选择Python3.7进行实现。

逻辑步骤如下：

a.导入相关的工具包以及定义相关函数
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
import numpy as np

import scipy as sp
from scipy.optimize import leastsq

import matplotlib.pyplot as plt
%matplotlib inline

# 目标函数
def real_func(x):
    return np.sin(2*np.pi*x)

# 多项式
def fit_func(p, x):
    f = np.poly1d(p)
    return f(x)

# 残差
def residuals_func(p, x, y):
    ret = fit_func(p, x) - y
    return ret
b.生成线性函数数据
1
2
3
4
5
6
7
8
9
# 十个点
nums = 10
#等差
x = np.linspace(0, 1, nums)
# print(x)
x_points = np.linspace(0, 1, 1000)
# 加上正态分布噪音的目标函数的值
y_ = real_func(x)
y = [np.random.normal(0, 0.1) + y1 for y1 in y_]
c.创建拟合函数并且可视化数据
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
def fitting(M=0):
    """
    M    为多项式的次数
    """
    # 随机初始化多项式参数
    p_init = np.random.rand(M + 1)
    # 最小二乘法 求导得到参数theta   
    p_lsq = leastsq(residuals_func, p_init, args=(x, y))
    print('Fitting Parameters:', p_lsq[0])

    # 可视化
    plt.plot(x_points, real_func(x_points), label='real')
    plt.plot(x_points, fit_func(p_lsq[0], x_points), label='fitted curve')
    plt.plot(x, y, 'bo', label='noise')
    plt.legend()
    return p_lsq
d.训练数据
1
2
# M=0
p_lsq_0 = fitting(M=0)

1
2
# M=1
p_lsq_1 = fitting(M=1)

1
2
# M=3
p_lsq_3 = fitting(M=3)

1
2
# M=9
p_lsq_9 = fitting(M=9)

​ 随着M阶数的增加，函数对点的拟合程度也随之提高，当M=3时，函数拟合程度的泛化能力算是最好的。当M=9的时候，函数把每一个点都进行了拟合，并且都在线上，也就是说在训练集当中，准确率为100%，但是如果设置测试集，M=9拟合的函数一定没有M=3的函数的泛化能力强。下面，介绍正则化。当我们追求准确率时，又要保证函数的泛化能力，正则化是一个很好的技巧，可以有效的防止过拟合问题。当然不同的正则化，适应着不同的数据分布，随着我们篇章的书写，会继续讨论相关问题。

1
2
3
4
5
6
regularization = 0.0001

def residuals_func_regularization(p, x, y):
    ret = fit_func(p, x) - y
    ret = np.append(ret,np.sqrt(0.5 * regularization * np.square(p)))  # L2范数作为正则化项
    return ret
1
2
3
# 最小二乘法,加正则化项
p_init = np.random.rand(9 + 1)
p_lsq_regularization = leastsq(residuals_func_regularization, p_init, args=(x, y))
1
2
3
4
5
6
7
8
9
#可视化正则化后与不适用正则化之间的效果的效果
plt.plot(x_points, real_func(x_points), label='real')
plt.plot(x_points, fit_func(p_lsq_9[0], x_points), label='fitted curve')
plt.plot(
    x_points,
    fit_func(p_lsq_regularization[0], x_points),
    label='regularization')
plt.plot(x, y, 'bo', label='noise')
plt.legend()

从上图，可以看出，正则化对于模型的泛化能力还是很有效果的，当然增加拟合数据x，y的数量，拟合程度降维更好，有兴趣的同学可以试一试。

5. 深度学习实现最小二乘法

现在我们使用深度学习方法去实现最小二乘法。这里选择TensorFlow 2.0进行实现。可以使用在线编译器google-colab-notebook.

逻辑步骤如下：

a.导入工具包
1
2
%tensorflow_version 2.x
import tensorflow as tf
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#目标函数
def real_func(x):
  y = tf.sin(2*np.pi*x)
  return y

class Model(object):
    def __init__(self):
        self.W_1 = tf.Variable(tf.random.uniform([1]))  # 随机初始化参数
        self.W_2 = tf.Variable(tf.random.uniform([1]))
        self.W_3 = tf.Variable(tf.random.uniform([1]))
        self.b = tf.Variable(tf.random.uniform([1]))

    def __call__(self,x):
      y_pred = self.W_1*x + self.W_2*x**2 + self.W_3*x**3 + self.b
      return y_pred

def loss(predicted_y, desired_y):
  return tf.reduce_mean(tf.square(predicted_y - desired_y))
1
2
3
4
5
6
7
8
9
10
#拟合数据
nums = 100000 #产生数据的点数
t_x = tf.linspace(-1.0,1.0,nums)
#加入正太分布噪声
noise = tf.random.truncated_normal([nums,], mean = .0, stddev = .1, dtype=tf.float32, seed = 38)
t_y = real_func(t_x) + noise
# print(t_y)
#原图像
t_x_points = tf.linspace(-1.0,1.0,10000)
t_y_points = real_func(t_x_points)
1
2
3
4
5
6
7
8
9
10
11
12
#这里建立3阶拟合函数，可以自行增加拟合函数的阶数，注意需要更新的参数W.
def train_step(model, inputs, outputs, learning_rate):
  with tf.GradientTape() as tape:
    current_loss = loss(model(inputs), outputs)

  dW_1, dW_2, dW_3, db = tape.gradient(current_loss, [model.W_1, model.W_2, model.W_3, model.b])
  #加入优化器 能够提升 0.550 -- 0.506
  optimizer.apply_gradients(zip((dW_1, dW_2, dW_3, db), [model.W_1, model.W_2, model.W_3, model.b]))
  model.W_1.assign_sub(learning_rate * dW_1)
  model.W_2.assign_sub(learning_rate * dW_2)
  model.W_3.assign_sub(learning_rate * dW_3)
  model.b.assign_sub(learning_rate * db)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#加载模型变量
model = Model()
#加载优化算法 
optimizer = tf.keras.optimizers.Adam()

# 收集 W 和 b 的历史数值，用于显示
current_loss = 0
Ws1, Ws2, Ws3, bs = [], [], [], []
epochs = range(100)
for epoch in epochs:
  Ws1.append(model.W_1.numpy())
  Ws2.append(model.W_2.numpy())
  Ws3.append(model.W_3.numpy())
  bs.append(model.b.numpy())
  current_loss = loss(model(t_x), t_y)

  train_step(model, t_x, t_y, learning_rate=0.01)
  print('Epoch %2d: W1=%1.2f W2=%1.2f W3=%1.2f b=%1.2f, loss=%2.5f' %
         (epoch, Ws1[-1], Ws2[-1], Ws3[-1], bs[-1], current_loss))
print("current_loss %.3f" % current_loss)
1
2
3
4
5
6
7
8
9
#通过获得参数，实现拟合函数
def pre_f(x):
  f = Ws1[-1]*x + Ws2[-1]*x**2 + Ws3[-1]*x**3 + bs[-1]
  return f
#可视化拟合函数
plt.plot(t_x_points.numpy(), t_y_points.numpy(),label = 'real')
plt.plot(t_x.numpy(), t_y.numpy(),label = 'noise')
plt.plot(t_x.numpy(), pre_f(t_x).numpy(), label = 'pre')
plt.legend()

​ 可见，在小量的数据当中，深度学习的效果，不尽人意，和ML实现的效果相比，相差还是较大的。你可以把数据量以及迭代次数增加，可以增加拟合效果，特别是数据量的增加，可以有效的提升模型的效果。如果感兴趣，可以试一试。当然深度学习的正则化亦有效果的。

6. 总结三者之间的联系以及区别。

​ 通过统计学，机器学习，深度学习上理论和实践的讲述。我们很容易从中总结如下几点它们之间的联系与区别。

a.联系

​ 1.每一种学习都离不开坚实数学知识基础。

​ 2.机器学习和深度学习代码实践上都由统计学学习理论支撑，即三者在实现最小二乘法上的逻辑一致。

b.区别

​ 1.可以发现统计学学习在参数，数据量较大和函数阶数较大后，实现最小二乘法基本是不可能的事情，进行对数值较大的计算，并且在生活中，数据量是较大的，人工完成，工具量极大，并且准确率不能够保证，所以产生机器学习。

​ 2.机器学习没有把统计学理论上的结果，完完全全的实现，只是拟合出来的函数，能够逼近正确结果，也就是近似值。所以在机器学习当中，模型的鲁棒性(robust)显得格外重要。当然传统的机器学习在数据量上的运行亦是有限制的。所以在一些方面产生了深度学习去训练。比如nlp，cv…

3.由于数据量已经超级巨大(部署与分布式)，并且在某些方面深度学习有较大的优势，因为它在数据量极大的时候，效果更好，更快并且自动化提取特征。

​ 总之，我们从统计学学习->机器学习->深度学习,依次完成了最小二乘法。希望大家能够在流程与实践当中，能够体会到数学基础的重要性，以及三者之间的实质联系与区别。

参考文献

@(*.*)@ @(^.^)@ @(*.*)@ @(^.^)@ @(*.*)@ @(^.^)@ @(*.*)@ @(^.^)@ @(*.*)@ @(^.^)@

统计学习方法概论

最小二乘法原理

TensorFlow 2.0基础