机器学习：线性回归

章节安排

1. 背景介绍
2. 均方根误差MSE
3. 最小二乘法
4. 梯度下降
5. 编程实现

背景

---

生活中大多数系统的输入输出关系为线性函数，或者在一定范围内可以近似为线性函数。在一些情形下，直接推断输入与输出的关系是较为困难的。因此，我们会从大量的采样数据中推导系统的输入输出关系。典型的单输入单输出线性系统可以用符号表示为：

\[y=f(x)=kx+b \]

其中，\(k\)为斜率，反应了当输入量\(x\)变化时，输出\(y\)的变化与输入\(x\)变化的比值；\(b\)反应了当系统没有输入（或输入为\(0\)）时，系统的输出值。

数据一般称观测数据或采样数据，这两种说法具有一定的侧重点，观测倾向于客观系统，例如每天的涨潮水深；采样倾向于主观系统，例如，对弹簧施加10N的压力，观察弹簧的形变量。

对于但输入单输出系统，数据可以表示为：

\[O=\{o_i\}_N=\{x_i,y_i\}_N \]

或

\[S=\{s_i\}_N=\{x_i,y_i\}_N \]

其中符号\(O\)对应observation(观测)、符号\(S\)对应sampling(采样),\(\{o_i\}_N\)中\(o_i\)表示采样序列中的每一个元素，\(N\)表示序列中元素的个数，\(x_i\)表示系统输入，\(y_i\)表示系统输出

在系统的推导过程中，一般称推导的结果为对实际系统的估计或近似，用符号记为\(\hat{y}=\hat{f}(x)\)。对于单个采样点，系统的误差定义为：对该采样输入，输出的真实值与输出的预测值的差为误差。用数据公式表示为：

\[\varepsilon_i = y_i-\hat{y_i}=y_i-\hat{f}(x_i) \]

对于整体采样序列，一种经典的误差是均方根误差（Mean Squared Error, MSE），其数学公式为：

\[\text{MSE}=\sum_{i=1}^{N}\varepsilon_i^2 \]

在推导系统输入输出关系，通常有两种方法，一种是基于数值推导的方法，一种是基于学习的方法。本文分别以最小二乘法和梯度下降为例讲解两种方法。

MSE

对于单个采样点的情形，MSE退化为方差的平方，即：

\[\text{MSE}=\varepsilon^2=(y-\hat y)^2 \]

假定参数\(b\)为常量，仅考虑MSE与参数的关系，有

\[\varepsilon^2=(kx+b-y)^2=x^2(k+\frac{b-y}{x})^2 \]

易得，MSE是关于\(k\)的二次函数，且该二次函数有唯一的零点：\(k_0=-(b-y)/x\)

对于多个点的情形，对每个点\(\{s_i\}=\{x_i,y_i\}\)，\(\varepsilon_i^2\)均可表示为关于\(k\)的二次函数，有：

\[\text{MSE}=\sum_{i=0}^{N}\varepsilon_i^2=\sum_{i=0}^{N}\big(x_i^2(k+\frac{b-y_i}{x_i})^2\big)=\sum_{i=0}^{N}\big(a_ik^2+b_ik+c_i\big)=Ak^2+Bk+C \]

即：序列的MSE也为关于参数\(k\)的二次函数，并且，\(MSE\geq0\)，当且仅当\((b-y_i)/x_i=M\)为常数时不等式取等。

可以很容易证明MSE也是关于参数\(b\)的二次函数

开口向上的二次函数有两个重要的性质：

1. 导数为\(0\)的点，为其最小值点。

\[f(x_i)= \min{f(x)}\iff f'(x_i)=0 \]

2. 任意点距离最小值点的距离与其导数值成正比，方向为导数方向的反方向

\[x_i-x_{\min}\propto -f'(x_i) \]

性质1、2分别是最小二乘法、梯度下降法的理论基础/依据。

最小二乘法

---

最小二乘法基于MSE进行设计，其思想为，找到一组参数，使得MSE关于每个参数的偏导为0，对于一元输入的情形，即：

\[\begin{align} \frac{\partial\text{MSE}}{\partial k}&=0 \tag{3.1}\\ \frac{\partial\text{MSE}}{\partial b}&=0 \tag{3.2} \end{align} \]

首先化简公式\((3.2)\)

\[\begin{align*} \frac{\partial\text{MSE}}{\partial b}&=\frac{1}{N}\sum_{i=1}^{N}\frac{\partial (\varepsilon_i^2)}{\partial b}\\ &=\frac{1}{N}\sum_{i=1}^{N}2\epsilon_i\cdot\frac{\partial}{\partial b} (\varepsilon_i)\\ &=\frac{2}{N}\sum_{i=1}^{N}\epsilon_i\cdot\frac{\partial}{\partial b} (kx_i+b-y_i)\\ &=\frac{2}{N}\sum_{i=1}^{N}(kx_i+b-y_i)\\ &=\frac{2}{N}\Big(k\sum_{i=1}^{N}x_i+Nb-\sum_{i=1}^{N}y_i\Big) \end{align*} \]

由公式\((3.2)\)有：

\[\begin{align*} \frac{2}{N}\Big(k\sum_{i=1}^{N}x_i+Nb-\sum_{i=1}^{N}y_i\Big)&=0\\ b&=\frac{1}{N}\Big(\sum_{i=1}^{N}y_i-k\sum_{i=1}^{N}x_i\Big) \tag{3.3} \end{align*} \]

其次化简公式\(3.1\)

\[\begin{align*} \frac{\partial\text{MSE}}{\partial k}&=\frac{1}{N}\sum_{i=1}^{N}\frac{\partial (\varepsilon_i^2)}{\partial k}\\ &=\frac{1}{N}\sum_{i=1}^{N}2\epsilon_i\cdot\frac{\partial}{\partial k} (\varepsilon_i)\\ &=\frac{2}{N}\sum_{i=1}^{N}\epsilon_i\cdot\frac{\partial}{\partial k} (kx_i+b-y_i)\\ &=\frac{2}{N}\sum_{i=1}^{N}x_i(kx_i+b-y_i)\\ &=\frac{2}{N}\Big(k\sum_{i=1}^{N}x_i^2+b\sum_{i=1}^{N}x_i-\sum_{i=1}^{N}x_iy_i\Big) \end{align*} \]

代入公式\((3.1),(3.3)\)有：

\[\begin{align*} \frac{2}{N}\Big(k\sum_{i=1}^{N}x_i^2+b\sum_{i=1}^{N}x_i-\sum_{i=1}^{N}x_iy_i\Big)&=0\\ k\sum_{i=1}^{N}x_i^2+\frac{1}{N}\sum_{i=1}^{N}x_i\Big(\sum_{i=1}^{N}y_i-k\sum_{i=1}^{N}x_i\Big)-\sum_{i=1}^{N}x_iy_i&=0\\ k\Big(\sum_{i=1}^{N}x_i^2-\frac{1}{N}\big(\sum_{i=1}^{N}x_i\big)^2\Big)&=\sum_{i=1}^{N}x_iy_i-\frac{1}{N}\sum_{i=1}^{N}x_i\sum_{i=1}^{N}y_i\\ k&=\frac{N\sum x_i^2-\big(\sum x_i\big)^2} {N\sum x_iy_i-\sum x_i\sum y_i} \tag{3.4} \end{align*} \]

公式\((3.3),(3.4)\)即为最小二乘法的参数公式

梯度下降

---

对于学习机器学习的初学者，我们首先讨论最简单的情形：基于单个采样点的学习。

二次函数具有重要性质：任意点距离最小值点的距离与其导数值成正比

\[x_i-x_{\min}\propto -f'(x_i) \]

基于该性质，我们可以可以设计参数更新公式如下

\[\begin{align*} \Delta k_t&=-\lambda\frac{\partial\varepsilon_i^2}{\partial k}\\ &=-\lambda(2\varepsilon_i\frac{\partial\varepsilon_i}{\partial k})\\ &=-\lambda(2\varepsilon_i x_i) \end{align*} \]

\[\begin{align*} \Delta b_t&=-\lambda\frac{\partial\varepsilon_i^2}{\partial b}\\ &=-\lambda(2\varepsilon_i\frac{\partial\varepsilon_i}{\partial b})\\ &=-\lambda(2\varepsilon_i) \end{align*} \]

故有参数更新公式：

\[\begin{align*} \varepsilon_i&=y-(kx_i+b_i)\tag{4.1}\\ k&:=k-\lambda(2\varepsilon_i x_i) \tag{4.2}\\ b&:=v--\lambda(2\varepsilon_i)\tag{4.3} \end{align*} \]

其中\(\lambda\)为学习率，一般取\(0.1\sim10^{-6}\)

常数\(2\)是可以缺省的，可以视为学习率放大了两倍。

编程实现

建议读者按照如下方法创建头文件、定义函数
typedef.h ：定义变量类型
random_point.h：生成随机点
least_square.h：最小二乘法的实现
gradient_descent.h：梯度下降方法的实现

类型定义

---

首先我们需要定义采样点，以及采样点序列类型。
采样点是包含\(x\)、\(y\)两个值的数据类型。同时，为方便使用，定义别名Point
采样点序列，或者称数据，可以存储为类型为Point的vector

struct SamplePoint{
  float x;
  float y;
}
using Point = SamplePoint;

using Data = std::vector<Point>;

对于直线，其包含\(k\)，\(b\)两个参数，同时，为了方便调用，定义括号运算符()重载

struct LinearFunc{
  float k;
  float b;
  float operator()(float x){
    return k*x+b;
  }
}
using Line = LinearFunc;
using Func = LinearFunc;

数据生成

---

采用random库中的normal_distribution随机数引擎

#include <random>
#include <cmath>
#include "typedef.h"

Data generatePoints(const Func& func, float sigma, float a, float b, int numPoints) {
    Data points;
    std::random_device rd;
    std::mt19937 gen(rd());
    // std::uniform_real_distribution<> distX(a, b); // 均匀分布
    std::normal_distribution<> distX((a + b) / 2, (b - a) / 2.8); // 正态分布
    std::normal_distribution<> distY(0, sigma);

    for (int i = 0; i < numPoints; ++i) {
        float x = distX(gen);
        float y = func(x) + distY(gen);
        points.push_back({ x, y });
    }

    return points;
}

该方法接受五个输入，分别是：

1. func：函数，自变量\(x\)与自变量\(y\)的关系
2. sigma：\(y\)的观测值与真实值的误差的方差
3. a、b：生成的数据范围的参考上下界，决定了生成数据的宽度，同时，绝大多数数据将位于此区间
4. numPoints：点的个数

最小二乘法

---

最小二乘法仅需接受一个输入：数据Data，同时返回数据。

\[\begin{align*} k&=\frac{N\sum x_i^2-\big(\sum x_i\big)^2} {N\sum x_iy_i-\sum x_i\sum y_i} \tag{3.4}\\ b&=\frac{1}{N}\Big(\sum_{i=1}^{N}y_i-k\sum_{i=1}^{N}x_i\Big) \tag{3.3} \end{align*} \]

在实现中，需要遍历采样数据，并分别进行累加计算\(\sum x_i\)、\(\sum y_i\)、\(\sum x_i^2\)和\(\sum x_iy_i\)

Line Least_Square(const Data& data) {
  Line line;

  float s_x = 0.0f;
  float s_y = 0.0f;
  float s_xx = 0.0f;
  float s_xy = 0.0f;

  float n = static_cast<float>(data.size());

  for (const auto& p : data) {
    s_x += p.x;
    s_y += p.y;
    s_xx += p.x * p.x;
    s_xy += p.x * p.y;
  }
	
  line.k = (n * s_xy - s_x * s_y) / (n * s_xx - s_x * s_x);
  line.b = (s_y - line.k * s_x) / n;

  return line;
}

梯度下降

---

梯度下降法是一种学习方法。对参数的估计逐渐向最优估计靠近。在本例中表现为，MSE逐渐降低。
首先实现单步的迭代，在该过程中，遍历所有的采样数据，依据参数更新公式对参数进行修正。

\[\begin{align*} \varepsilon_i&=y-(kx_i+b_i)\tag{4.1}\\ k&:=k-\lambda(2\varepsilon_i x_i) \tag{4.2}\\ b&:=v--\lambda(2\varepsilon_i)\tag{4.3} \end{align*} \]

梯度下降法需要一个给定的初值，对于线性函数，除了人工生成、随机初值外，一种方式是，假定为正比例函数，以估计\(k\)，假定为常函数，以估计\(b\)，公式如下：

\[\begin{align*} k_0&=\sum y_i/\sum x_i \tag{5.1}\\ b_0&=\sum y_i/ N \tag{5.2} \end{align*} \]

在本例中，设定为对初值进行100次迭代后得到最终估计，读者可根据实际情况调整，在学习度设计的合适的情况下，一般迭代次数在\(50\sim200\)次

#include "typedef.h"

constexpr float eps = 1e-1;
constexpr float lambda = 1e-5;

void GD_step(Func& func, const Data& data) {
  for (const auto& p : data) {
    float error = func(p.x) - p.y;
    func.k -= lambda * error * p.x;
    func.b -= lambda * error;
  }
}

Func Gradient_Descent(Func& func, const Data& data) {
  float s_x = 0, s_y = 0;
  for (const auto& p : data) {
    s_x += p.x;
    s_y += p.y;
  }

  Line line;
  line.k = s_y / s_x;
  line.b = s_y / data.size();

  float lambda = 1e-5f;

  for (size_t _ = 0; _ < 100; _++) {
    GD_step(line, data);
  }

  return line;
}

附录

nan问题

该问题有两种产生的原因，参数更新符号错误及学习率过高。

参数更新符号错误
在更新公式中，如果错误的使用+号，或者采用\(\hat y-y\)计算\(\varepsilon_i\)，都将会导致参数向误差更大的方向更新，经过了数次迭代后，与真实值的距离越来越远，最终产生nan。

\[k:=k-\lambda(2\varepsilon_i x_i) \]

学习率过高
如下图，当学习率设置的过高时，新的参数组\(\{k_{t+1},b_{t+1}\}\)将比旧参数\(\{k_{t},b_{t}\}\)带来更大的估计误差（红色箭头），而良好的学习率是使得估计误差逐渐下降的