机器学习之梯度下降算法(python实现)

一点关于学习梯度下降算法的感受

首先就个人体验来看，在本科期间在学习运筹与优化这门课程就接触过梯度下降算法，那个时候就是单一的求一个具体多变量函数的最小值问题。现在在机器学习里面的梯度下降算法更多的是根据训练数据集去找到一个合适的拟合函数，那么怎么样才能找到一个拟合度高的函数使之在测试集上有较好的泛化能力，此时就需要定义一个拟合函数模型和一个损失函数，当损失函数的取值最小时就可以得到这个在训练集拟合度好的函数了。得到这个函数模型之后就可以去做相关测试数据的预测了。其中求损失函数的最小值过程就是我们现在所讲的梯度下降算法迭代过程了。

一、梯度下降的相关概念

在详细了解梯度下降的算法之前，我们先看看相关的一些概念。

1. 步长（Learning rate）：步长决定了在梯度下降迭代的过程中，每一步沿梯度负方向前进的长度。用上面下山的例子，步长就是在当前这一步所在位置沿着最陡峭最易下山的位置走的那一步的长度。

2.特征（feature）：指的是样本中输入部分，比如2个单特征的样本 （ x ( 0 ) , y ( 0 ) ） , （ x ( 1 ) , y ( 1 ) ） （x^{(0)},y^{(0)}）,（x^{(1)},y^{(1)}） （x(0),y(0)）,（x(1),y(1)）,则第一个样本特征为 x ( 0 ) x^{(0)} x(0)，第一个样本输出为 y ( 0 ) y^{(0)} y(0)。

3. 假设函数（hypothesis function）即拟合函数：在监督学习中，为了拟合输入样本，而使用的假设函数，记为 h θ ( x ) h_{\theta}(x) hθ​(x)。比如对于单个特征的m个样本 （ x ( i ) , y ( i ) ） ( i = 1 , 2 , . . . m ) （x^{(i)},y^{(i)}）(i=1,2,...m) （x(i),y(i)）(i=1,2,...m),可以采用拟合函数如下： h θ ( x ) = θ 0 + θ 1 x h_{\theta}(x) = \theta_0+\theta_1x hθ​(x)=θ0​+θ1​x。

4. 损失函数（loss function）：为了评估模型拟合的好坏，通常用损失函数来度量拟合的程度。损失函数极小化，意味着拟合程度最好，对应的模型参数即为最优参数。在线性回归中，损失函数通常为样本输出和假设函数的差取平方。比如对于m个样本 （ x i , y i ） ( i = 1 , 2 , . . . m ) （x_i,y_i）(i=1,2,...m) （xi​,yi​）(i=1,2,...m),采用线性回归，损失函数为： J ( θ 0 , θ 1 ) = ∑ i = 1 m ( h θ ( x i ) − y i ) 2 J(\theta_0, \theta_1) = \sum\limits_{i=1}^{m}(h_\theta(x_i) - y_i)^2 J(θ0​,θ1​)=i=1∑m​(hθ​(xi​)−yi​)2
　　　　其中 x i x_i xi​表示第i个样本特征， y i y_i yi​表示第i个样本对应的输出， h θ ( x i ) h_\theta(x_i) hθ​(xi​)为假设函数。

二、梯度下降的详细算法(代数方式描述)

1. 先决条件： 确认优化模型的假设函数和损失函数。

比如对于线性回归，假设函数表示为 h θ ( x 1 , x 2 , . . . x n ) = θ 0 + θ 1 x 1 + . . . + θ n x n h_\theta(x_1, x_2, ...x_n) = \theta_0 + \theta_{1}x_1 + ... + \theta_{n}x_{n} hθ​(x1​,x2​,...xn​)=θ0​+θ1​x1​+...+θn​xn​, 其中 θ i \theta_i θi​(i = 0,1,2… n)为模型参数， x i x_i xi​ (i = 0,1,2… n)为每个样本的n个特征值。这个表示可以简化，我们增加一个特征 x 0 x_0 x0​=1 ，这样 h θ ( x 0 , x 1 , . . . x n ) = ∑ i = 0 n θ i x i h_\theta(x_0, x_1, ...x_n) = \sum\limits_{i=0}^{n}\theta_{i}x_{i} hθ​(x0​,x1​,...xn​)=i=0∑n​θi​xi​。

同样是线性回归，对应于上面的假设函数，损失函数为： J ( θ 0 , θ 1 . . . , θ n ) = 1 2 m ∑ j = 1 m ( h θ ( x 0 ( j ) , x 1 ( j ) , . . . x n ( j ) ) − y j ) 2 J(\theta_0, \theta_1..., \theta_n) = \frac{1}{2m}\sum\limits_{j=1}^{m}(h_\theta(x_0^{(j)}, x_1^{(j)}, ...x_n^{(j)}) - y_j)^2 J(θ0​,θ1​...,θn​)=2m1​j=1∑m​(hθ​(x0(j)​,x1(j)​,...xn(j)​)−yj​)2

2. 算法相关参数初始化：主要是初始化 θ 0 , θ 1 . . . , θ n \theta_0, \theta_1..., \theta_n θ0​,θ1​...,θn​,算法终止距离ε以及步长α。在没有任何先验知识的时候，我喜欢将所有的θ初始化为0， 将步长初始化为0.001。在调优的时候再 优化。

3. 算法过程：

1）确定当前位置的损失函数的梯度，对于 θ i \theta_i θi​,其梯度表达式如下： ∂ ∂ θ i J ( θ 0 , θ 1 . . . , θ n ) \frac{\partial}{\partial\theta_i}J(\theta_0, \theta_1..., \theta_n) ∂θi​∂​J(θ0​,θ1​...,θn​)
2）用步长乘以损失函数的梯度，得到当前位置下降的距离，即 α ∂ ∂ θ i J ( θ 0 , θ 1 . . . , θ n ) \alpha\frac{\partial}{\partial\theta_i}J(\theta_0, \theta_1..., \theta_n) α∂θi​∂​J(θ0​,θ1​...,θn​)对应于前面登山例子中的某一步。

3）确定是否所有的 θ i \theta_i θi​,梯度下降的距离都小于ε，如果小于ε则算法终止，当前所有的θi(i=0,1,…n)即为最终结果。否则进入步骤4.

4）更新所有的 θ \theta θ，对于 θ i \theta_i θi​，其更新表达式如下。更新完毕后继续转入步骤1. θ i = θ i − α ∂ ∂ θ i J ( θ 0 , θ 1 . . . , θ n ) \theta_i = \theta_i - \alpha\frac{\partial}{\partial\theta_i}J(\theta_0, \theta_1..., \theta_n) θi​=θi​−α∂θi​∂​J(θ0​,θ1​...,θn​)
下面用线性回归的例子来具体描述梯度下降。假设我们的样本是 ( x 1 ( 0 ) , x 2 ( 0 ) , . . . x n ( 0 ) , y 0 ) , ( x 1 ( 1 ) , x 2 ( 1 ) , . . . x n ( 1 ) , y 1 ) , . . . ( x 1 ( m ) , x 2 ( m ) , . . . x n ( m ) , y m ) (x_1^{(0)}, x_2^{(0)}, ...x_n^{(0)}, y_0), (x_1^{(1)}, x_2^{(1)}, ...x_n^{(1)},y_1), ... (x_1^{(m)}, x_2^{(m)}, ...x_n^{(m)}, y_m) (x1(0)​,x2(0)​,...xn(0)​,y0​),(x1(1)​,x2(1)​,...xn(1)​,y1​),...(x1(m)​,x2(m)​,...xn(m)​,ym​),损失函数如前面先决条件所述： J ( θ 0 , θ 1 . . . , θ n ) = 1 2 m ∑ j = 0 m ( h θ ( x 0 ( j ) , x 1 ( j ) , . . . x n ( j ) ) − y j ) 2 J(\theta_0, \theta_1..., \theta_n) = \frac{1}{2m}\sum\limits_{j=0}^{m}(h_\theta(x_0^{(j)}, x_1^{(j)}, ...x_n^{(j)})- y_j)^2 J(θ0​,θ1​...,θn​)=2m1​j=0∑m​(hθ​(x0(j)​,x1(j)​,...xn(j)​)−yj​)2

则在算法过程步骤1中对于θi 的偏导数计算如下：

∂ ∂ θ i J ( θ 0 , θ 1 . . . , θ n ) = 1 m ∑ j = 0 m ( h θ ( x 0 ( j ) , x 1 ( j ) , . . . x n ( j ) ) − y j ) x i ( j ) \frac{\partial}{\partial\theta_i}J(\theta_0, \theta_1..., \theta_n)= \frac{1}{m}\sum\limits_{j=0}^{m}(h_\theta(x_0^{(j)}, x_1^{(j)}, ...x_n^{(j)}) - y_j)x_i^{(j)} ∂θi​∂​J(θ0​,θ1​...,θn​)=m1​j=0∑m​(hθ​(x0(j)​,x1(j)​,...xn(j)​)−yj​)xi(j)​
　　　　由于样本中没有 x 0 x_0 x0​上式中令所有的 x 0 j x_0^{j} x0j​为1.

步骤4中θi的更新表达式如下：
θ i = θ i − α 1 m ∑ j = 0 m ( h θ ( x 0 ( j ) , x 1 ( j ) , . . . x n j ) − y j ) x i ( j ) \theta_i = \theta_i - \alpha\frac{1}{m}\sum\limits_{j=0}^{m}(h_\theta(x_0^{(j)}, x_1^{(j)}, ...x_n^{j}) - y_j)x_i^{(j)} θi​=θi​−αm1​j=0∑m​(hθ​(x0(j)​,x1(j)​,...xnj​)−yj​)xi(j)​
　　　　从这个例子可以看出当前点的梯度方向是由所有的样本决定的，加 1 m \frac{1}{m} m1​是为了好理解。由于步长也为常数，他们的乘机也为常数，所以这里 α 1 m \alpha\frac{1}{m} αm1​可以用一个常数表示。

三、代码实现(python)和相关效果图

说不多说，到了大家最关心的代码实现的部分了，其实上面所说的算法过程既是批量梯度下降算法的过程。还有就是说实话这是我认真写过的有详细注释的代码了，欧克，肝就完了！
下面展示全部代码 。

import numpy as np
import matplotlib.pyplot as plt

# 定义数据集的大小为30
m = 30
# x0为1
x = np.ones((m, 1))
# 设置三个特征向量取值数组的拼接
x = np.concatenate((x, np.arange(1, m + 1).reshape(m, 1), np.arange(2, 2 * m + 1, 2).reshape(m, 1)), axis=1)
y_data = np.loadtxt('data.txt', delimiter=',')
# 从数组提取y的取值
y = np.array(y_data[:, 1][0:30]).reshape(m, 1)
# 设置学习效率
alpha = 0.001
# 设置模型参数theta的初始迭代点(一般都是从0开始)
theta = np.zeros((3, 1))
# 分别定义一个存储损失函数值和迭代次数的列表
costlist, num = [], []


# 定义一个代价函数
def costfunction(theta, x, y):
    error = np.dot(x, theta) - y
    return (1 / 2 * m) * (np.dot(error.T, error))


# 定义一个梯度迭代的函数，即偏导数部分
def gradientdescent(theta, x, y):
    error = np.dot(x, theta) - y
    return (1 / m) * (np.dot(x.T, error))


# 下降迭代过程函数
def GDworkfunction(theta, alpha, x, y):
    gd = gradientdescent(theta, x, y)
    # 迭代10万次终止
    for i in range(100000):
        theta = theta - alpha * gd
        cost1 = costfunction(theta, x, y)
        costlist.append(cost1)
        num.append(i)
        gd = gradientdescent(theta, x, y)
    return theta


consulttheta = GDworkfunction(theta, alpha, x, y)
print('cost值取得最小时的函数为：y={}+{}*x1+{}*x2'.format(consulttheta[0], consulttheta[1], consulttheta[2]))
cost = costfunction(consulttheta, x, y)[0][0]
print('此时的最小损失值为{}'.format(cost))
costlist1 = []
# 对cost值进行特征缩小，缩小的范围自己主观决定，怎么方便怎么来！
for j in costlist:
    costlist1.append(j / 4500)

plt.figure(num=3, figsize=(8, 5))
# 在绘图时，需要对cost值进行降维成一维，因为保存的是二维的数组
plt.plot(num, np.squeeze(np.array(costlist1)))
plt.xlim((1, 80000))
plt.ylim(1, 10)
plt.xlabel('iterate_nums')
plt.ylabel('costvalue')
# 设置x轴刻度表示范围
my_x_ticks = np.arange(1, 80000, 6000)
plt.xticks(my_x_ticks)
plt.savefig('F:\\DataBuilding\\BGD\\BGD.png')  # 保存图片
plt.show()

import sys; print('Python %s on %s' % (sys.version, sys.platform))
sys.path.extend(['F:\\pycharmspace\\pycharm-project\\日常dome', 'F:/pycharmspace/pycharm-project/日常dome'])
PyDev console: starting.
Python 3.8.3 (tags/v3.8.3:6f8c832, May 13 2020, 22:37:02) [MSC v.1924 64 bit (AMD64)] on win32
runfile('F:/pycharmspace/pycharm-project/日常dome/Algorithm/MygradientDecentdome.py', wdir='F:/pycharmspace/pycharm-project/日常dome/Algorithm')
cost值取得最小时的函数为：y=[7.51769165]+[0.01716882]*x1+[0.03433763]*x2
此时的最小损失值为18629.353835283197

emmm，做完之后就觉得这个损失值的最小值大的惊人，可能是数据选取的数据跨度大，不过不影响实际算法目的。接下来上迭代次数和损失函数的效果图：

就这张图来说，当下降迭代次数到达6000次的时候就开始逐渐趋向最小损失值了。下面附上 y i y_i yi​的数据集

17.592,9.1302,11.854,6.8233,11.886,4.3483,12,6.5987,3.8166,3.2522,15.505,3.1551,7.2258,0.71618,3.51295.3048,0.56077,3.6518,5.3893,3.1386,21.767,4.263,5.1875,3.0825,22.638,13.501,7.0467,14.692,24.147

四、总结

梯度下降算法机器学习的求最优值这方面无论怎么改动都有不可代替的优势，因为其稳定，效率较高，且结果一定会得出来且基本都接近或者就是最优解，当然我们也可以用数据分析神器SPSS来做一个的分析结果来验证代码跑出来的结果是否一致。这个大家可以做做看，感兴趣的可以把我这个代码跑跑看。
这是我第一次自写的博客，写的比较不怎么成熟，以后我会继续更新一些有关机器学习的算法，比如说求全局最优值的随机游走算法，和遗传算法等…