笔记1：A simple example of Gradient Descent using Python

引言

写这篇blog的目的主要是让自己熟悉一下在李宏毅老师的机器学习课程中学到的Gradient Descent的原理，以及记录下每一步是如何用编程实现的。

Linear Regression

对于Linear Regression，有如下形式的映射关系：
$$y=b+w\ast x\text{\hspace{1em}(1)}$$
事实上，整个Linear Regression的思想就是，在给定的一个function set中，找到与Training Data匹配最好的一个函数，并且这个函数用在Test Data上也非常准确。衡量一个函数的精确度，我们用Loss函数来表示。对于当前这个很简单的例子，我们的函数集就是形如(1)式的一群函数，我们的Loss函数用下式来表示：
$$L=\sum _n(y_{data}-(b+w{x}_{data}){)}^2\text{\hspace{1em}(2)}$$
我们要找出最精确的可以描述Training Data的函数，换言之，我们要找到合适的b和w使得Loss函数最小。具体实现方法就是接下来要讲到的Gradient Descent。

Gradient Descent原理

先从一维入手，假设Loss函数只是一个一元函数，即可表示为$L=L(w)$，那么我们要怎么找到其最小值呢？就是求导。我们随机取一个$w_0$，计算在这个点处的$\frac{dL}{dw}$，如果这个导数大于零，说明此时Loss函数随$w$呈上升趋势，那么我们需要向小的方向取新的$w$的值，反之，我们就需要向大的方向取值。那么，可以得到迭代公式：
$$w_n=w_{n-1}-\eta \frac{dL}{dw}{|}_{w=w_{n-1}}\text{\hspace{1em}(3)}$$
其中$\eta$为learning rate，可以人为设定。利用该公式进行迭代则可得到合适的值。
（PS：此处没有认真听讲，据李宏毅老师，经过一定的迭代，在Linear Regression中不会出现只是local optimal的答案）

将此迭代公式推广至二维，即$L=L(b,w)$，经过类似的思考，可以得到如下的迭代公式：
$$w_n=w_{n-1}-\eta \frac{\partial L}{\partial w}{|}_{w=w_{n-1},b=b_{n-1}}\text{\hspace{1em}(4)}$$
$$b_n=w_{n-1}-\eta \frac{\partial L}{\partial b}{|}_{w=w_{n-1},b=b_{n-1}}\text{\hspace{1em}(5)}$$
这就是Gradient Descent的原理了。接下来用一个简单的例子来熟悉一下这个思路转化成代码是怎么样的。

代码示例

导入相关函数库

import numpy as np
import matplotlib.pyplot as plt

我们需要numpy来进行矩阵的运算操作，以及matplotlib.pyplot进行结果的可视化。

导入Training Data

由于只是简单实现一下Gradient Descent，为了更方便理解原理，Training Data直接用定义的方式导入。

#Training data
#y_data = b + w * x_data
x_data = [338., 333., 328., 207., 226., 25., 179., 60., 208., 606.]
y_data = [640., 633., 619., 393., 428., 27., 193., 66., 226., 1591.]

建立model (function set)

x = np.arange(-200, -100, 1)  #[-200, -199, ... , -100]
y = np.arange(-5, 5, 0.1)  #[-5, -4.9, -4.8, ... , 5]
Z = np.zeros((len(x), len(y)))  #2-dimension matrix
X, Y = np.meshgrid(x, y)  #网格矩阵
for i in range(len(x)):
    for j in range (len(y)):
        b = x[i]
        w = y[j]
        Z[j][i] = 0  #事实上前面用zeros函数初始化过，这里不用再赋 0
        for n in range(len(x_data)):
            Z[j][i] = Z[j][i] + (y_data[n] - b - w*x_data[n])**2
        Z[j][i] = Z[j][i]/len(x_data)

用numpy的arange函数分别建立b和w可取值的数组，Z为一个二维矩阵，每个元素即为特定b和w对应的Loss值，meshgrid函数主要是形成两个二维矩阵，位置相同的两个元素组成一点的坐标。b和w遍历可取的所有值，并计算每一次取到的值对应的Loss值。

Iteration

#ydata = b + w * xdata
#初始化b, w, learning rate, iteration times
b = -120
w = -4
lr = 0.000001
iteration = 100000

b_history = [b]
w_history = [w]

#Iterations
for i in range(iteration):
    b_grad = 0.0
    w_grad = 0.0
    #计算偏导
    for n in range(len(x_data)):
        b_grad = b_grad - 2.0*(y_data[n] - b - w*x_data[n])*1.0
        w_grad = w_grad - 2.0*(y_data[n] - b - w*x_data[n])*x_data[n]
    #迭代主体
    b = b - lr*b_grad
    w = w - lr*w_grad
    
    #store parameters for plotting
    b_history.append(b)
    w_history.append(w)

这里我们就开始迭代了，首先要取定b, w, learning rate和迭代次数。b_history和w_history是记录b和w取到过的值。我们接下来就在迭代次数内进行迭代，用b_grad和w_grad记录Loss函数对b和w求偏导所得数值。由简单的数学知识，我们可以推导出来：
$$\frac{\partial L}{\partial b}=-2\sum _n(y_{data}-(b+w{x}_{data}))\text{\hspace{1em}(6)}$$
$$\frac{\partial L}{\partial w}=-2x_{data}\sum _n(y_{data}-(b+w{x}_{data}))\text{\hspace{1em}(7)}$$
经过迭代，b_history和w_history数组的最后一个数即为当前迭代次数可求出来的效果最好的b和w。

结论可视化

#plot the figure
plt.contourf(x, y, Z, 50, alpha=0.5, cmap=plt.get_cmap('jet'))
plt.plot([-188.4], [2.67], 'x', ms=12, markeredgewidth=3, color='orange')	#已知解
plt.plot(b_history, w_history, 'o-', ms=3, lw=1.5, color='black')			#每一次迭代的解
plt.xlim(-200, -100)	#横坐标
plt.ylim(-5, 5)			#纵坐标
plt.xlabel(r'$b$', fontsize=16)
plt.ylabel(r'$w$', fontsize=16)
plt.show()
print(str(b_history[len(b_history)-1]) + ', ' + str(w_history[len(w_history)-1]))

这一部分主要是将Linear Regression的效果进行一个可视化，颜色的冷暖程度体现了此处所取b和w对应Loss值的大小，颜色越暖，Loss值越大。至于contourf的函数，笔者并不是非常了解，只是照搬了李宏毅老师的代码而已，若有机会应该好好学习一下pyplot。