在求解机器学习算法的模型参数,即无约束优化问题时,梯度下降(Gradient Descent)是最常采用的方法之一,另一种常用的方法是最小二乘法。
1. 梯度
在微积分里面,对多元函数的参数求∂偏导数,把求得的各个参数的偏导数以向量的形式写出来,就是梯度。比如函数f(x,y), 分别对x,y求偏导数,求得的梯度向量就是 $(∂f/∂x, ∂f/∂y)^T$ ,简称 $grad f(x,y)$ 或者$▽f(x,y)$ 。对于在点 $(x_0,y_0)$ 的具体梯度向量就是 $(∂f/∂x_0, ∂f/∂y_0)^T$ .或者 $▽f(x_0,y_0)$ ,如果是3个参数的向量梯度,就是 $(∂f/∂x, ∂f/∂y,∂f/∂z)^T$ ,以此类推 。
那么这个梯度向量求出来有什么意义呢?它的意义从几何意义上讲,就是函数变化增加最快的地方。具体来说,对于函数 $f(x,y)$ ,在点 $(x_0,y_0)$ ,沿着梯度向量的方向就是 $(∂f/∂x_0, ∂f/∂y_0)^T$ 的方向是$f(x,y)$ 增加最快的地方。或者说,沿着梯度向量的方向,更加容易找到函数的最大值。反过来说,沿着梯度向量相反的方向,也就是$ -(∂f/∂x_0, ∂f/∂y_0)^T$ 的方向,梯度减少最快,也就是更加容易找到函数的最小值。
2. 梯度下降与梯度上升
在机器学习算法中,在最小化损失函数时,可以通过梯度下降法来一步步的迭代求解,得到最小化的损失函数,和模型参数值。反过来,如果我们需要求解损失函数的最大值,这时就需要用梯度上升法来迭代了。
梯度下降法和梯度上升法是可以互相转化的。比如我们需要求解损失函数 $f(θ)$ 的最小值,这时我们需要用梯度下降法来迭代求解。但是实际上,我们可以反过来求解损失函数 $-f(θ)$ 的最大值,这时梯度上升法就派上用场了。
3. 梯度下降法算法详解
3.1 梯度下降的直观解释
首先来看看梯度下降的一个直观的解释。比如我们在一座大山上的某处位置,由于我们不知道怎么下山,于是决定走一步算一步,也就是在每走到一个位置的时候,求解当前位置的梯度,沿着梯度的负方向,也就是当前最陡峭的位置向下走一步,然后继续求解当前位置梯度,向这一步所在位置沿着最陡峭最易下山的位置走一步。这样一步步的走下去,一直走到觉得我们已经到了山脚。当然这样走下去,有可能我们不能走到山脚,而是到了某一个局部的山峰低处。
从上面的解释可以看出,梯度下降不一定能够找到全局的最优解,有可能是一个局部最优解。当然,如果损失函数是凸函数,梯度下降法得到的解就一定是全局最优解。
3.2 梯度下降的相关概念
在详细了解梯度下降的算法之前,我们先看看相关的一些概念。
-
步长(Learning rate)
步长决定了在梯度下降迭代的过程中,每一步沿梯度负方向前进的长度。用上面下山的例子,步长就是在当前这一步所在位置沿着最陡峭最易下山的位置走的那一步的长度。
-
特征(feature)
指的是样本中输入部分,比如2个单特征的样本 $(x^{(0)},y^{(0)}),(x^{(1)},y^{(1)})$ ,则第一个样本特征为 $x^{(0)}$ , 第一个样本输出为 $y^{(0)}$ 。
-
假设函数(hypothesis function)
在监督学习中,为了拟合输入样本,而使用的假设函数,记为 $ h_{\theta}(x)$ 。比如对于单个特征的m个样本 $(x^{(i)},y^{(i)})(i=1,2,…m)$ , 可以采用拟合函数如下:$h_{\theta}(x) = \theta_0+\theta_1x$ .
-
损失函数(loss function)
为了评估模型拟合的好坏,通常用损失函数来度量拟合的程度。损失函数极小化,意味着拟合程度最好,对应的模型参数即为最优参数。在线性回归中,损失函数通常为样本输出和假设函数的差取平方。比如对于 m 个样本 $(x_i,y_i)(i=1,2,…m)$ , 采用线性回归,损失函数为:
\[J(\theta_0, \theta_1) = \sum\limits_{i=1}^{m}(h_\theta(x_i) - y_i)^2\]其中 $x_i$ 表示第i个样本特征,$y_i$ 表示第 $i$ 个样本对应的输出,$h_θ(x_i)$ 为假设函数。
3.3 梯度下降的详细算法
梯度下降法的算法可以有代数法和矩阵法(也称向量法)两种表示,如果对矩阵分析不熟悉,则代数法更加容易理解。不过矩阵法更加的简洁,且由于使用了矩阵,实现逻辑更加的一目了然。这里先介绍代数法,后介绍矩阵法。
3.3.1 梯度下降法的代数方式描述
1.先决条件: 确认优化模型的假设函数和损失函数。
比如对于线性回归,假设函数表示为 $h_\theta(x_1, x_2, …x_n) = \theta_0 + \theta_{1}x_1 + … + \theta_{n}x_{n}$ , 其中 $θ_i (i = 0,1,2… n)$ 为模型参数,$x_i (i = 0,1,2… n)$ 为每个样本的n个特征值。这个表示可以简化,我们增加一个特征$x_0=1$ ,这样 $h_\theta(x_0, x_1, …x_n) = \sum\limits_{i=0}^{n}\theta_{i}x_{i}$ 。
同样是线性回归,对应于上面的假设函数,损失函数为:
\[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\]2.算法相关参数初始化
主要是初始化 $\theta_0, \theta_1…, \theta_n$ ,算法终止距离 $ε$ 以及步长 $α$ 。在没有任何先验知识的时候,我喜欢将所有的 $θ$ 初始化为0, 将步长初始化为 1。在调优的时候再优化。
3.算法过程
1)确定当前位置的损失函数的梯度,对于 $θ_i$ ,其梯度表达式如下:
\[\frac{\partial}{\partial\theta_i}J(\theta_0, \theta_1..., \theta_n)\]2)用步长乘以损失函数的梯度,得到当前位置下降的距离,即 $\alpha\frac{\partial}{\partial\theta_i}J(\theta_0, \theta_1…, \theta_n)$ 对应于前面登山例子中的某一步。
3)确定是否所有的 $\theta_i$ , 梯度下降的距离都小于 $ε$ , 如果小于 $ε$ 则算法终止,当前所有的$θ_i(i=0,1,…n)$ 即为最终结果。否则进入步骤4.
4)更新所有的 $θ$,对于$θi$ ,其更新表达式如下。更新完毕后继续转入步骤1.
\[\theta_i = \theta_i - \alpha\frac{\partial}{\partial\theta_i}J(\theta_0, \theta_1..., \theta_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)$ ,
损失函数如前面先决条件所述:
\[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\]则在算法过程步骤1中对于$θ_i$ 的偏导数计算如下:
\[\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)}\]由于样本中没有 $x_0$上式中令所有的 $x_0^{j}$ 为1.
步骤4中$θ_i$ 的更新表达式如下:
\[\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)}\]从这个例子可以看出当前点的梯度方向是由所有的样本决定的,加 $\frac{1}{m}$ 是为了好理解。由于步长也为常数,他们的乘机也为常数,所以这里 $\alpha\frac{1}{m}$ 可以用一个常数表示。
3.3.2 梯度下降法的矩阵方式描述
这一部分主要讲解梯度下降法的矩阵方式表述,相对于3.3.1的代数法,要求有一定的矩阵分析的基础知识,尤其是矩阵求导的知识。
1.先决条件:
和3.3.1类似, 需要确认优化模型的假设函数和损失函数。对于线性回归,假设函数
\[h_\theta(x_1, x_2, ...x_n) = \theta_0 + \theta_{1}x_1 + ... + \theta_{n}x_{n}\]的矩阵表达方式为:$h_\mathbf{\theta}(\mathbf{X}) = \mathbf{X\theta}$ , 其中, 假设函数 $h_\mathbf{\theta}(\mathbf{X})$ 为$m\times1$的向量 , $θ$ 为$(n+1) \times 1$ 的向量 , 里面有n+1个代数法的模型参数。$\mathbf{X}$ 为 $m\times (n+1)$ 维的矩阵。m 代表样本的个数,n+1代表样本的特征数。
2.算法相关参数初始化
$θ$ 向量可以初始化为默认值,或者调优后的值。算法终止距离 $ε$ ,步长 $α$ 和3.3.1比没有变化。
3.算法过程:
1)确定当前位置的损失函数的梯度,对于 $θ$ 向量,其梯度表达式如下:
\[\frac{\partial}{\partial\mathbf\theta}J(\mathbf\theta)\]2)用步长乘以损失函数的梯度,得到当前位置下降的距离,即 $\alpha\frac{\partial}{\partial\theta}J(\theta)$ 对应于前面登山例子中的某一步。
3)确定 $θ$ 向量里面的每个值,梯度下降的距离都小于 $ε$ ,如果小于$ε$ 则算法终止,当前 $θ$ 向量即为最终结果。否则进入步骤4.
4)更新 $θ$ 向量,其更新表达式如下。更新完毕后继续转入步骤1.
\[\mathbf\theta= \mathbf\theta - \alpha\frac{\partial}{\partial\theta}J(\mathbf\theta)\]线性回归的例子来描述具体的算法过程。
损失函数对于 $θ$ 向量的偏导数计算如下:
\[\frac{\partial}{\partial\mathbf\theta}J(\mathbf\theta) = \mathbf{X}^T(\mathbf{X\theta} - \mathbf{Y})\]步骤4中 $ θ$ 向量的更新表达式如下:$\mathbf\theta= \mathbf\theta - \alpha\mathbf{X}^T(\mathbf{X\theta} - \mathbf{Y})$
对于3.3.1的代数法,可以看到矩阵法要简洁很多。这里面用到了矩阵求导链式法则,和两个矩阵求导的公式。
这里面用到了矩阵求导链式法则,和两个个矩阵求导的公式。
公式1:$\frac{\partial}{\partial\mathbf{x}}(\mathbf{x^Tx}) =2\mathbf{x}\;\;x为向量 $
公式2:$ \nabla_Xf(AX+B) = A^T\nabla_Yf,\;\; Y=AX+B,\;\;f(Y)为标量$
3.4 梯度下降的算法调优
在使用梯度下降时,需要进行调优。哪些地方需要调优呢?
-
算法的步长选择 : 在前面的算法描述中,我提到取步长为1,但是实际上取值取决于数据样本,可以多取一些值,从大到小,分别运行算法,看看迭代效果,如果损失函数在变小,说明取值有效,否则要增大步长。前面说了。步长太大,会导致迭代过快,甚至有可能错过最优解。步长太小,迭代速度太慢,很长时间算法都不能结束。所以算法的步长需要多次运行后才能得到一个较为优的值。
-
算法参数的初始值选择。 初始值不同,获得的最小值也有可能不同,因此梯度下降求得的只是局部最小值;当然如果损失函数是凸函数则一定是最优解。由于有局部最优解的风险,需要多次用不同初始值运行算法,关键损失函数的最小值,选择损失函数最小化的初值。
-
归一化。由于样本不同特征的取值范围不一样,可能导致迭代很慢,为了减少特征取值的影响,可以对特征数据归一化,也就是对于每个特征x,求出它的期望 $\overline{x}$ 和标准差std(x),然后转化为:
\[\frac{x - \overline{x}}{std(x)}\]这样特征的新期望为0,新方差为1,迭代速度可以大大加快。
4. 梯度下降法大家族(BGD,SGD,MBGD)
4.1 批量梯度下降法(Batch Gradient Descent)
批量梯度下降法,是梯度下降法最常用的形式,具体做法也就是在更新参数时使用所有的样本来进行更新,这个方法对应于前面3.3.1的线性回归的梯度下降算法,也就是说3.3.1的梯度下降算法就是批量梯度下降法。
\[\theta_i = \theta_i - \alpha\sum\limits_{j=1}^{m}(h_\theta(x_0^{(j)}, x_1^{(j)}, ...x_n^{(j)}) - y_j)x_i^{(j)}\]由于我们有m个样本,这里求梯度的时候就用了所有m个样本的梯度数据。
4.2 随机梯度下降法(Stochastic Gradient Descent)
随机梯度下降法,其实和批量梯度下降法原理类似,区别在与求梯度时没有用所有的m个样本的数据,而是仅仅选取一个样本j来求梯度。对应的更新公式是:
\[\theta_i = \theta_i - \alpha (h_\theta(x_0^{(j)}, x_1^{(j)}, ...x_n^{(j)}) - y_j)x_i^{(j)}\]机梯度下降法,和4.1的批量梯度下降法是两个极端,一个采用所有数据来梯度下降,一个用一个样本来梯度下降。自然各自的优缺点都非常突出。对于训练速度来说,随机梯度下降法由于每次仅仅采用一个样本来迭代,训练速度很快,而批量梯度下降法在样本量很大的时候,训练速度不能让人满意。对于准确度来说,随机梯度下降法用于仅仅用一个样本决定梯度方向,导致解很有可能不是最优。对于收敛速度来说,由于随机梯度下降法一次迭代一个样本,导致迭代方向变化很大,不能很快的收敛到局部最优解。
那么,有没有一个中庸的办法能够结合两种方法的优点呢?有!这就是4.3的小批量梯度下降法。
4.3 小批量梯度下降法(Mini-batch Gradient Descent)
小批量梯度下降法是批量梯度下降法和随机梯度下降法的折衷,也就是对于m个样本,我们采用x个样子来迭代,1<x<m。一般可以取x=10,当然根据样本的数据,可以调整这个x的值。对应的更新公式是:
\[\theta_i = \theta_i - \alpha \sum\limits_{j=t}^{t+x-1}(h_\theta(x_0^{(j)}, x_1^{(j)}, ...x_n^{(j)}) - y_j)x_i^{(j)}\]5. 梯度下降法和其他无约束优化算法的比较
在机器学习中的无约束优化算法,除了梯度下降以外,还有前面提到的最小二乘法,此外还有牛顿法和拟牛顿法。
梯度下降法和最小二乘法相比,梯度下降法需要选择步长,而最小二乘法不需要。梯度下降法是迭代求解,最小二乘法是计算解析解。如果样本量不算很大,且存在解析解,最小二乘法比起梯度下降法要有优势,计算速度很快。但是如果样本量很大,用最小二乘法由于需要求一个超级大的逆矩阵,这时就很难或者很慢才能求解解析解了,使用迭代的梯度下降法比较有优势。
梯度下降法和牛顿法/拟牛顿法相比,两者都是迭代求解,不过梯度下降法是梯度求解,而牛顿法/拟牛顿法是用二阶的海森矩阵的逆矩阵或伪逆矩阵求解。相对而言,使用牛顿法/拟牛顿法收敛更快。但是每次迭代的时间比梯度下降法长。
6. python 实现梯度下降
import numpy as np
import matplotlib.pyplot as plt
# Size of the points dataset.
m = 20
# Points x-coordinate and dummy value (x0, x1).
X0 = np.ones((m, 1))
X1 = np.arange(1, m+1).reshape(m, 1)
X = np.hstack((X0, X1))
# Points y-coordinate
y = np.array([
3, 4, 5, 5, 2, 4, 7, 8, 11, 8, 12,
11, 13, 13, 16, 17, 18, 17, 19, 21
]).reshape(m, 1)
# The Learning Rate alpha.
alpha = 0.01
def error_function(theta, X, y):
'''Error function J definition.'''
diff = np.dot(X, theta) - y
return (1./2*m) * np.dot(diff.T, diff)
def gradient_function(theta, X, y):
'''损失函数求偏导'''
diff = np.dot(X, theta) - y
return (1./m) * np.dot(X.T, diff)
def gradient_descent(X, y, alpha):
'''Perform gradient descent.'''
# theta 初始化的值
theta = np.array([1, 1]).reshape(2, 1)
# gradient 损失函数求偏导后的结果
gradient = gradient_function(theta, X, y)
while not np.all(np.absolute(gradient) <= 1e-5):
# 更新theta
theta = theta - alpha * gradient
gradient = gradient_function(theta, X, y)
return theta
optimal = gradient_descent(X, y, alpha)
print('optimal:', optimal)
print('error function:', error_function(optimal, X, y)[0,0])
x_v = np.linspace(1, 21, 1000)
# 拟合直线
y_v = optimal[0][0] + optimal[1][0]*x_v
plt.scatter(X1.reshape(1, -1),y.reshape(1,-1),)
plt.plot(x_v, y_v, color = "r")
plt.savefig("GD.png", dpi = 600)
plt.show()
效果:
源码:传送门