直线拟合

目标

我们有很多数据,这些数据是有一系列的点构成。我们认为i这些点应该是有关系的,我们想要找到这种关系。我们假设,这些点是直线关系,现在我们要找到这条直线的参数。

如图,这些数据的分布如下,我们可以先剧透一下,这些数据是用 $y = (2x + 1) \pm 5$ 生成的随机数据。所以,我们应该拟合一条直线 $y = 2x + 1$。

机器学习系列-直线拟合-数据

定义直线与误差计算方法

首先,我们定义这条直线 $y = ax + b$。

我们定义 $y$ 代表数据集的真实值,使用 $\hat{y}$ 代表我们拟合出来的直线的计算值。

这里我们使用平方法计算误差,或者叫作二乘法。对于一个点而言,真实值与计算值的误差 $e = (y - \hat{y})^2$ 。那么,在整个数据集上,真实值与计算值的误差 $E = \sum_{i=0}^{n} (y_i - \hat{y_i})^2$ 。

直线的寻找方法

显然,我们寻找直线的过程就是最小化误差的过程。当我们的误差很小的时候,我们就可以认为,我们找到了这条直线。我们可以将误差计算方法看作一个函数,我们要最小化这个函数。关于如何最小化函数,我们可以考虑使用梯度下降。

梯度下降的计算公式为,$\eta$ 为学习率:

$$ x_{k+1} = x_k - \eta \nabla f(x_k) $$

对我们的误差函数使用梯度下降:

$$ \begin{bmatrix} a_{k+1} \\ b_{k+1} \end{bmatrix} = \begin{bmatrix} a_k \\ b_k \end{bmatrix} - \eta \nabla e(a_k, b_k) $$

误差函数的梯度向量:

$$ \nabla e(a, b) = \begin{bmatrix} \frac{\partial e(a, b)}{\partial a} \\ \frac{\partial e(a, b)}{\partial b} \end{bmatrix} $$

现在计算迭代公式,有了迭代公式,我们就可以使用计算机及进行拟合运算。我们需要分别考虑对于一个点和对于整个数据集的情况。

我们先来考虑,数据只有一个点的情况:

$$ \begin{aligned} e(a, b) &= (y - \hat{y})^2 \\ &= y^2 + \hat{y}^2 -2 y \hat{y} \\ &= y^2 + (ax + b)^2 - 2y(ax + b) \\ &= y^2 + (ax)^2 + b^2 + 2abx - 2y(ax + b) \\ &= y^2 + a^2x^2 + b^2 + 2abx - 2axy - 2by \\ &= x^2a^2 + b^2 + 2xab - 2xya - 2yb + y^2 \\ \nabla e(a, b) &= \begin{bmatrix} 2x^2a + 2xb - 2xy \\ 2xa + 2b - 2y \end{bmatrix} \\ \begin{bmatrix} a_{k+1} \\ b_{k+1} \end{bmatrix} &= \begin{bmatrix} a_k - \eta \frac{\partial e(a, b)}{\partial a} \\ b_k - \eta \frac{\partial e(a, b)}{\partial a} \end{bmatrix} \\ &= \begin{bmatrix} a_k - \eta (2x^2a_k + 2xb_k - 2xy) \\ b_k - \eta (2xa_k + 2b_k - 2y) \end{bmatrix} \end{aligned} $$

然后考虑对于整个数据集:

$$ \begin{aligned} E(a, b) &= \sum_{i=0}^{n} (y_i - \hat{y_i})^2 \\ &= \sum_{i=0}^{n} (x_i^2a^2 + b^2 + 2x_iab - 2x_iy_ia - 2y_ib + y_i^2) \\ &= a^2 \sum_{i=0}^{n} x_i^2 + b^2 \sum_{i=0}^{n} 1 + ab \sum_{i=0}^{n} 2x_i - 2a \sum_{i=0}^{n} x_iy_i - 2b \sum_{i=0}^{n} y_i + \sum_{i=0}^{n} y_i^2 \\ \nabla E(a, b) &= \begin{bmatrix} 2a \sum_{i=0}^{n} x_i^2 + b \sum_{i=0}^{n} 2x_i - 2 \sum_{i=0}^{n} x_iy_i \\ a \sum_{i=0}^{n} 2x_i + 2b \sum_{i=0}^{n} 1 - 2 \sum_{i=0}^{n} y_i \end{bmatrix} \\ &= \begin{bmatrix} \sum_{i=0}^{n} (2x_i^2a + 2x_ib - 2x_iy_i) \\ \sum_{i=0}^{n} (2x_ia + 2b - 2y_i) \end{bmatrix} \\ \begin{bmatrix} a_{k+1} \\ b_{k+1} \end{bmatrix} &= \begin{bmatrix} a_k - \eta \frac{\partial E(a, b)}{\partial a} \\ b_k - \eta \frac{\partial E(a, b)}{\partial a} \end{bmatrix} \\ &= \begin{bmatrix} a_k - \eta \sum_{i=0}^{n} (2x_i^2a + 2x_ib - 2x_iy_i) \\ b_k - \eta \sum_{i=0}^{n} (2x_ia + 2b - 2y_i) \end{bmatrix} \end{aligned} $$

我们仔细观察,对于整个数据集的迭代公式,这个式子应该是可以化简的。$2x_i^2a + 2x_ib - 2x_iy_i$ 就是 $\frac {\partial e(a, b)} {\partial a}$ 在点 $i$ 的表达,$2x_ia + 2b - 2y_i$ 就是 $\frac {\partial e(a, b)} {\partial b}$ 在点 $i$ 的表达。

换言之,$\sum_{i=0}^{n} (2x_i^2a + 2x_ib - 2x_iy_i)$,就是将整个数据集上的 $\frac {\partial e(a, b)} {\partial a}$ 进行求和,$\sum_{i=0}^{n} (2x_ia + 2b - 2y_i)$,就是将整个数据集上的 $\frac {\partial e(a, b)} {\partial b}$ 进行求和。

现在,对于整个数据集上的迭代公式,我们可以有以下形式:

$$ \begin{aligned} \begin{bmatrix} a_{k+1} \\ b_{k+1} \end{bmatrix} &= \begin{bmatrix} a_k - \eta \sum_{i=0}^{n} \frac {\partial e_i(a_k, b_k)} {\partial a_k} \\ b_k - \eta \sum_{i=0}^{n} \frac {\partial e_i(a_k, b_k)} {\partial b_k} \end{bmatrix} \end{aligned} $$

编写数据生成代码

import random


class Generator():
    def gen_data(self):
        data = []
        x = 0
        y = 0
        while x <= 50:
            y = (2 * x + 1) + random.uniform(-5, 5)
            data.append((x, y))
            x += 0.1
        return data

    def __gen_expectLine(self):
        data = []
        x = 0
        y = 0
        while x <= 50:
            y = (2 * x + 1)
            data.append((x, y))
            x += 0.1
        return data

    def __draw_write(self, points, filename):
        f = open(filename, "w")
        for point in points:
            f.write("%.4f  %.4f\n" % point)
        f.close()

    def draw(self, points, filename):
        self.__draw_write(points, filename)

    def draw_expectLine(self, filename):
        a = self.__gen_expectLine()
        self.__draw_write(a, filename)

编写误差函数代码

class E():
    def e(self, dataset, a, b):
        e_sum = 0
        for data in dataset:
            x = data[0]
            y = data[1]
            e_sum += (y - (a * x + b)) ** 2
        return e_sum

    def partial_a(self, dataset, a, b):
        partial_a_sum = 0
        for data in dataset:
            x = data[0]
            y = data[1]
            partial_a_sum += 2 * x ** 2 * a + 2 * x * b - 2 * x * y
        return partial_a_sum

    def partial_b(self, dataset, a, b):
        partial_b_sum = 0
        for data in dataset:
            x = data[0]
            y = data[1]
            partial_b_sum += 2 * x * a + 2 * b - 2 * y
        return partial_b_sum

    def __draw_gen(self, dataset):
        points = []
        a = -50
        b = -50
        z = 0
        while a <= 50:
            while b <= 50:
                z = self.e(dataset, a, b)
                points.append((a, b, z))
                b += 0.1
            a += 0.1
            points.append("\n")
            b = -50
        return points

    def __draw_write(self, points, filename):
        f = open(filename, "w")
        for point in points:
            if point == "\n":
                f.write("\n")
            else:
                f.write("%.4f  %.4f  %.4f\n" % point)
        f.close()

    def draw(self, dataset, filename):
        a = self.__draw_gen(dataset)
        self.__draw_write(a, filename)

编写直线代码

class MyLine():
    a = 0
    b = 0

    def func(self, x):
        return self.a * x + self.b

    def train(self, dataset, e, lr):
        change_a = lr * e.partial_a(dataset, self.a, self.b)
        change_b = lr * e.partial_b(dataset, self.a, self.b)
        self.a -= change_a
        self.b -= change_b
        return e.e(dataset, self.a, self.b)

    def __draw_gen(self):
        points = []
        x = 0
        while x <= 50:
            y = self.func(x)
            points.append((x, y))
            x = x + 0.1
        return points

    def __draw_write(self, points, filename):
        f = open(filename, "w")
        for point in points:
            f.write("%.4f  %.4f\n" % point)
        f.close()

    def draw(self, filename):
        a = self.__draw_gen()
        self.__draw_write(a, filename)

数据可视化代码

请注意,我们使用gnuplot进行图片的绘制,请确保您的系统已经安装了该软件。

set terminal png small size 900,900
set output "机器学习系列-直线拟合-dataset.png"
set title  "dataset"            #图片标题
set xlabel "x"                  #X轴标题
set ylabel "y"                  #Y轴标题
set grid                        #显示网格

plot "data_dataset.txt" using 1:2 title "dataset" with points pointtype 7 linetype 0
set terminal png small size 900,900
set output "机器学习系列-直线拟合-result.png"
set title  "result"             #图片标题
set xlabel "x"                  #X轴标题
set ylabel "y"                  #Y轴标题
set grid                        #显示网格

plot \
    "data_dataset.txt" using 1:2 title "dataset" with points pointtype 7 linetype 0 ,\
    "data_expectLine.txt" using 1:2 title "y = 2x + 1" with lines linetype 1 ,\
    "data_myLine.txt" using 1:2 title "myline" with lines linetype 2
set terminal png small size 900,900
set output "机器学习系列-直线拟合-loss.png"
set title  "loss"                   #图片标题
set xlabel "count"                  #X轴标题
set ylabel "loss"                   #Y轴标题
set grid                            #显示网格

plot "data_loss.txt" using 1:2 title "loss" with lines linetype 3
set terminal png small size 900,900
set output "机器学习系列-直线拟合-e.png"
set title  "e"                  # 图片标题
set xlabel "a"                  # X轴标题
set ylabel "b"                  # Y轴标题
set zlabel "e"                  # z轴标题

set pm3d at b
set ticslevel 1.5
set isosample 100,100

splot \
    "data_e.txt" using 1:2:3 title "e",\
    "data_parameter.txt" using 1:2:3 title "parameter" with points pointtype 7 linetype 2

运行代码

if __name__ == "__main__":
    generator = Generator()
    myLine = MyLine()
    e = E()

    dataset = generator.gen_data()
    loss_data = []
    parameter_data = []

    for i in range(200):
        loss = myLine.train(dataset, e, 0.000001)
        loss_data.append((i, loss))
        parameter_data.append((myLine.a, myLine.b, loss))
        print("%.4f  %.4f  %.4f" % (myLine.a, myLine.b, loss))

    generator.draw(dataset, "data_dataset.txt")
    generator.draw_expectLine("data_expectLine.txt")
    myLine.draw("data_myLine.txt")
    e.draw(dataset, "data_e.txt")

    f = open("data_loss.txt", "w")
    for point in loss_data:
        f.write("%.4f  %.4f\n" % point)
    f.close()

    f = open("data_parameter.txt", "w")
    for point in parameter_data:
        f.write("%.4f  %.4f  %.4f\n" % point)
    f.close()
#!/bin/bash

python main.py

cat gnuplot_dataset.conf | gnuplot
cat gnuplot_result.conf  | gnuplot
cat gnuplot_loss.conf    | gnuplot
cat gnuplot_e.conf       | gnuplot

获取我们的结果

请注意,由于数据使用随机数进行生成,每次运行的结果都不尽相同,但应该是类似的。

1.6842  0.0508  52644.7306
1.9678  0.0595  5600.2232
2.0156  0.0612  4265.8164
2.0236  0.0617  4227.9213
2.0250  0.0619  4226.8002
2.0252  0.0621  4226.7221
2.0252  0.0624  4226.6737

...

2.0240  0.1026  4218.0977
2.0240  0.1028  4218.0545
2.0240  0.1030  4218.0112
2.0240  0.1032  4217.9680
2.0240  0.1034  4217.9248

机器学习系列-直线拟合-结果

机器学习系列-直线拟合-loss

机器学习系列-直线拟合-e

结果分析

我们的直线最终收敛到 $y = 2.0240x + 0.1034$,此时的误差为 $4217.9248$。为什么不是 $y = 2.000x + 1.000$,为什么没有把误差降为 $0$,为什么我们选择了 $0.000001$ 作为学习率。

首先我们考虑,为什么学习率设置为了 $0.000001$ 。我们先试一试学习率为 $0.0001$ 的情况。Numerical result out of range 直接数字溢出了。我们继续尝试 $0.000003$ ,这次没有溢出,但是,我们获取的结果显然不太对。

5.0715  0.1531  3846876.1614
-2.5092  -0.0739  8589715.5986
8.8220  0.2672  19186584.1247
-8.1153  -0.2409  42863038.4463
17.2017  0.5203  95763051.7882
-20.6411  -0.6157  213956906.6297
35.9245  1.0842  478035995.1913

...

111544799665044376218191084161859584.0000  3350704033848358181964575412322304.0000  5178055970267548973924251206062346011916491500681666419082537824927317229568.0000
-166731971421014279783506711335141376.0000  -5008476333181802238281115072724992.0000  11569267301314213277921265094749927721535367024547556014368944683295801606144.0000

让我们看一下图:

机器学习系列-直线拟合-e-2

很明显,误差的值太大了,进行梯度下降的时候,算出的 $\Delta a$ 和 $\Delta b$ 太大,以至于更新的时候出现了摇摆。

到底哪里出了问题,是我们的学习率选取出了问题还是误差函数的问题。应该是误差函数。根据我们的公式,我们会将整个数据集上各个点的误差进行累加。随着数据集的增加,最终的误差也会越来越大,直到超出计算机的表示的范围,超出梯度下降合理的更新范围。

也许,我们将误差求个算术平均会好一点。