梯度下降法求解线性回归--Numpy实现

梯度下降法求解一元线性回归

依然是这个房价预测的任务，这是一个一元线性回归问题，这次我们采用梯度下降法来求解它可以分为5步

• 第1步加载样本数据x，y
• 第2步设置超参数，在这个例子中，超参数包括学习率和迭代次数
• 第3步设置模型参数的初值$w_0$，$b_0$，这个初值可以是任意的
• 第4步训练模型使用迭代公式更新模型，参数迭代完成之后，以可视化的形式输出结果

这是程序流程图，因为有迭代运算，所以需要通过循环来实现。

这部分是梯度下降法的实现，首先设置$w_0$，$b_0$，设置循环变量，然后利用迭代公式不断更新，w和b，并且计算每一次迭代的损失，直到循环结束为止。

在上节课中我们知道超参数的设置非常重要，对训练结果有很大的影响，在训练之前我们往往并不知道这个超参数应该设置成多少，一般需要根据经验反复尝试，同时观察算法是否收敛，并且达到了我们需要的精度，下面我们就编程实现以上步骤。

第1步导入需要的库加载数据

x和y分别是面积和房价，把它们放在Numpy数组中，它们都是长度为16的一维数组

import numpy as np
import matplotlib.pyplot as plt

x = np.array([137.97, 104.50, 100.00, 124.32, 79.20, 99.00, 124.00, 114.00,
              106.69, 138.05, 53.75, 46.91, 68.00, 63.02, 81.26, 86.21])

y = np.array([145.00, 110.00, 93.00, 116.00, 65.32, 104.00, 118.00, 91.00,
              62.00, 133.00, 51.00, 45.00, 78.50, 69.65, 75.69, 95.30])

第2步设置超参数

• learn_rate是学习率，通常是一个很小的常数。
• iter是迭代次数迭代100次
• 如果每次迭代都输出结果输出就会很长，也没有必要，因此我们可以每10次迭代输出一次结果这个display_step就是用来设置输出结果的。间隔它不属于超参数，因为它的取值完全不会影响模型的训练，只是会改变显示的效果，

learn_rate = 0.00001
iter = 100
display_step = 10

第3步给模型参数w和b设置初值

这是numpy的随机数生成函数。返回一个正态分布的浮点数组，当参数为空时，随机生成一个数字。

np.random.seed(612)
w = np.random.rand()
b = np.random.rand()

第4步训练模型

这个mse是一个Python列表，用来保存每次迭代后的损失值。
下面使用for循环实现迭代循环变量，从0开始，循环101次。为了描述方便，当i等于10时，我们就说第10次迭代。

mse = []
for i in range(0, iter + 1):
    dL_dw = np.mean(x*(w * x + b - y))
    dL_db = np.mean(w * x + b - y)
    
    w = w - learn_rate * dL_dw
    b = b - learn_rate * dL_db
    
    pred = w * x + b
    Loss = np.mean(np.square(y-pred))/2
    
    mse.append(Loss)
    
    if i % display_step == 0:
        print("i: %i, Loss: %f, w: %f, b: %f" % (i, mse[i], w, b))

在循环体中，首先计算损失函数队w和b的偏导数

然后使用迭代公式更新w和b

到这里就已经实现了梯度下降法的一次迭代，可以进入下一次循环了。但是我们希望能够观察到每次迭代的结果，判断是否收敛或者什么时候开始收敛，因此需要使用每次迭代后的w和b计算损失，并且把它显示出来。

这是使用当前这次循环得到的w和b计算所有样本的房价的估计值，这里的x是一个长度为16的一维数组，保存着所有样本的面积，因此这个pred也是一个长度为16的，一维数组，是所有的房价的预测值。

然后使用房价的实际值和预测值计算均方误差。y和pred都是长度为16的一维数组，这部分运算的结果仍然是一个一维数组，然后对这个数组中的所有元素，求评均值得到一个数字，最后再乘以1/2把得到的均方误差加入列表mse。mse一开始是一个空的列表，以后每执行一次，循环列表中就增加一个元素，整个循环执行完之后，其中就有101个元素，分别是每次迭代的损失值，如果当前的循环次数能够被10整除，就显示均方误差，w的值和b的值，这样做是为了能够动态的观察模型的训练过程。

这是i=0时的输出。

可以看到这个损失的值非常庞大。在前30次循环中，损失以很快的速度下降，这是因为在远离机制点的地方，损失函数曲线很陡峭，而且更新的步长比较大，因此损失的下降速度很快。

到第40次循环开始，损失下降的速度越来越慢，这是因为随着不断接近极值点，损失函数的曲线越来越平缓，步长也越来越小，因此损失的下降也越来越小。

在第100次迭代时，损失的值虽然还在下降，但是差值已经非常小了，这时候虽然还没有达到极值点，但是已经非常接近极值点了，大家可以把循环次数修改的更大一些，看看损失函数是不是还会下降。

梯度下降法得到的数值解是一个近似值，在收敛之后，只要达到了精度要求，就可以停止迭代，否则可以继续迭代，直到满足精度要求为止，下面把模型训练的结果以可视化的形式输出出来，

结果可视化–数据和模型

首先使用样本数据绘制销售记录散点图，然后绘制预测房价的散点图，这个pred是最后一次迭代之后计算出的房价的估计值，把它们连接在一起就是得到的线性模型。

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.figure()

plt.scatter(x, y, color="red", label="销售记录")
plt.scatter(x, pred, color="blue", label="梯度下降法")
plt.plot(x, pred, color="blue")

plt.xlabel("Area", fontsize=14)
plt.ylabel("Price", fontsize=14)

plt.legend(loc="upper left")
plt.show()

这是运行的结果。

其中红色的点是实际的销售房价，蓝色的点是预测出的房价。蓝色的直线是训练得到的模型，那么这个模型是否准确呢？这是在上一讲中我们计算出的解。解析解是一个精确的结果。现在我们可以把解析解对应的线性模型也绘制出来，进行比较。

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.figure()

plt.scatter(x, y, color="red", label="销售记录")
plt.scatter(x, pred, color="blue")
plt.plot(x, pred, color="blue", label="梯度下降法")
plt.plot(x, 0.89*x+5.41, color="green", label="解析法")

plt.xlabel("Area", fontsize=14)
plt.ylabel("Price", fontsize=14)

plt.legend(loc="upper left")
plt.show()

这里的w和b保留两位小数，这条绿色的直线就是解析法得到的模型。可以看到采用梯度下降法得到的模型和它有一定的偏差，但是在可以接受的范围之内，大家也可以尝试增加迭代次数，继续更新权值，让w和b更接近极值点。

这张图展示了整个迭代过程中模型直线的变化过程，这是使用w和b的初始值得到的直线，这是迭代100次之后得到的直线。后面的这些直线非常接近，很多都重叠在一起了，要做出这个图，只要在每次迭代之后，都增加一条plot函数绘制直线就可以了。

mse = []
for i in range(0, iter + 1):
    dL_dw = np.mean(x*(w * x + b - y))
    dL_db = np.mean(w * x + b - y)

    w = w - learn_rate * dL_dw
    b = b - learn_rate * dL_db

    pred = w * x + b
    Loss = np.mean(np.square(y-pred))/2

    mse.append(Loss)
    
    plt.plot(x, pred)

    if i % display_step == 0:
        print("i: %i, Loss: %f, w: %f, b: %f" % (i, mse[i], w, b))

结果可视化–损失变化

通过这张图，我们可以更加清楚的观察损失值的变化，图中的横坐标是迭代次数，纵坐标是损失值，然后逐渐减缓，会是这个图也非常简单，因为每次迭代的损失值都已经被存放在列表mse中了，现在只要把它们取出来，连在一起，就可以得到损失值变化的曲线图。

plt.figure()
plt.plot(mse)

plt.xlabel("Iteration", fontsize=14)
plt.ylabel("Loss", fontsize=14)

plt.show()

在这张图中，因为开始时的损失值非常大，所以纵坐标的刻度也很大，导致从第20次迭代以后损失的下降很难直接在这张图中看出来。如果想要观察到第20次迭代之后损失，可以把plot函数的参数修改一下第20次迭代之后的。

plt.figure()
plt.plot(range(20, 100), mse[20:100])

plt.xlabel("Iteration", fontsize=14)
plt.ylabel("Loss", fontsize=14)

plt.show()

可以看到现在从第20次迭代时的损失值148，开始显示损失变化曲线，在第20~40次迭代时损失下降的很快，40次之后就逐渐平坦了。

采用同样的方法也可以绘制出第40次迭代之后损失变化的曲线，大家可以继续尝试一下观察和总结损失下降的规律。

plt.figure()
plt.plot(range(40, 100), mse[40:100])

plt.xlabel("Iteration", fontsize=14)
plt.ylabel("Loss", fontsize=14)

plt.show()

另外要注意的是，这个图是在训练模型的过程中，损失函数的值变化的曲线，不是损失函数的图。损失函数本身应该和右边的图近似。

结果可视化–估计值&标签值

为了更加直观的展示预测值和实际值之间的差距，可以使用这样的图。这个图中的横坐标是样本序号，一共有16个点，每个点对应一套商品房，纵坐标是房价。

plt.plot(y, color="red", marker="o", label="销售记录")
plt.plot(pred, color="blue", marker="o", label="梯度下降法")

plt.legend()
plt.xlabel("Sample", fontsize=14)
plt.ylabel("PRICE", fontsize=14)
plt.show()

红色的数据点是样本数据，它们的纵坐标是每套商品房的实际销售价格。蓝色的点是我们通过模型预测出来的房价。
可以看到有些房价的估计很准确，有些存在着一定的偏差，这个图形也是使用plot函数绘制的，把参数mark的值设置为o，数据点就以这种大圆点的形式显示出来。

plt.plot(y, color="red", marker="o", label="销售记录")
plt.plot(pred, color="blue", marker="o", label="梯度下降法")
plt.plot(0.89*x+5.41, color="green", marker="o", label="解析法")

plt.legend()
plt.xlabel("Sample", fontsize=14)
plt.ylabel("PRICE", fontsize=14)
plt.show()

也可以把解析法得到的结果绘制出来进行对比，这条绿色的线就是解析法画出来的线。

为了使代码更加的简洁紧凑，结果更加便于观察，我们把这些图放在同一个画布中，显示首先设置画布尺寸，然后划分子图，在每个子图中分别绘制不同的图形。

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.figure(figsize=(20, 4))

plt.subplot(1, 3, 1)
plt.scatter(x, y, color="red", label="销售记录")
plt.plot(x, pred, color="blue", label="预测记录")
plt.xlabel("Area", fontsize=14)
plt.ylabel("Price", fontsize=14)
plt.legend(loc="upper left")

plt.subplot(1, 3, 2)
plt.plot(mse)
plt.xlabel("Iteration", fontsize=14)
plt.ylabel("Loss", fontsize=14)

plt.subplot(1, 3, 3)
plt.plot(y, color="red", marker="o", label="销售记录")
plt.plot(pred, color="blue", marker="o", label="预测记录")
plt.legend()
plt.xlabel("Sample", fontsize=14)
plt.ylabel("PRICE", fontsize=14)
plt.legend(loc="upper left")

plt.show()

这是运行的结果。

现在我们已经使用梯度下降法完成了对一元线性回归模型的训练，下面就可以使用这个模型对未知的数据进行预测了。