模型的优化与训练 --- 梯度下降法及其衍生

梯度下降（Grandient Descent）

梯度下降的核心原理：

• 函数的梯度方向表示了函数值增长速度最快的方向，那么和它相反的方向就可以看作是函数值减少速度最快的方向。对机器学习模型优化问题，当目标设定为求解目标函数最小值时，只要朝着梯度下降的方向前进，就能不断逼近最优值。

最简单的梯度下降算法 - 固定学习率的方法：

• 待优化的函数 $f(x)$
• 待优化函数的导数$g(x)$
• 变量$x$：保存当前优化过程中的参数值，优化开始时该变量将被初始化成某个数值，优化过程中这个变量会不断变化，直到它找到最小值
• 变量grad：保存变量$x$点处的梯度值
• 变量step：表示沿着梯度下降方向行进的步长，即学习率（Learining Rate），在优化过程中它将固定不变

def gd(x_start, step, g):		   #  Gradient Descent
	x = x_start
	for i in range(20):
		grad = g(x)
		x -= grad * step
		print '[ epoch {0} ]  grad = {1}, x = {2}'.format(i, grad, x)
		if abs(grad) < 1e-6:
			break;
	return x
	

由于优化的目标是寻找梯度为0的极值点，代码在每一轮迭代结束后衡量变量$x$所在的梯度值，因此当梯度值足够小的时，就认为$x$已经进入最优值附近一个极小的领域，$x$和最优值之间的差别不再明显，就可以停止优化了。$x$最后的值就是最优解的位置。

有兴趣的可以跑一下以下简单的二次函数的demo:

def f(x):
	return x * x - 2 * x + 1
def g(x):
	return 2 * x - 2
gd(5,0.1,g)

合适的步长非常重要

动量算法（Momentum）

在优化求解的过程中，动量代表了之前迭代优化量，它将在后面的优化过程中持续发挥作用，推动目标值前进。拥有了动量，一个已经结束的更新量不会立刻消失，只会以一定的形式衰减，剩下的能量将继续在优化过程中发挥作用。

def momentum(x_start, step, g ,discount = 0.7):
	x = np.array(x_start, dtype='float64')
	pre_grad = np.zeros_like(x)
	for i in range(50):
		grad = g(x)
		pre_grad = pre_grad * discount + grad * step
		x -= pre_grad
		print '[ epoch {0} ] grad = {1}, x = {2}'.format(i, grad ,x)
		if abs(sum(grad)) < 1e-6:
			break;
	return x

代码中多出了一个新变量 pre_grad，这个变量就是用于存储历史积累的动量，每一轮迭代动量都会乘以一个打折量（discount）做能量衰减，但它依然会被用于更新参数。
动量算法相比较梯度下降，能使得梯度在下降时减少左右震荡(震荡的方向是相反的，由于历史积累的动量，会相互抵消)，加快梯度下降速度。
但是动量优化存在一点问题，前面几轮的迭代过程中，梯度的震荡会比原来的还大，为了解决这个更强烈的抖动，就有了接下来的Nesterov算法

Nesterov算法

def nesterov(x_start, step, g, discount = 0.7)
	x = np.array(x_start,dtype='float64')
	pre_grad = np.zeros_like(x)
	for i in range(50):
		x_future = x - step * discount * pre_grad
		grad = g(x_future)
		pre_grad = pre_grad * 0.7 + grad
		x -= pre_grad * step
		print '[ Epoch {0} ] grad = {1} , x = {2} '.format(i, grad, x)
		if abs(sum(grad)) < 1e-6:
			break;
	return x

与动量算法相比，动量算法计算了当前目标点的梯度，而Nesterov算法计算了动量更新后优化点的梯度。当优化点已经积累了某个抖动方向的梯度后，这时对于动量算法来说，虽然当前点的梯度指向积累梯度的相反方向，但是量不够大，所以最终的优化方向还会在积累的方向上前进一段。对于Nesterov来说，如果按照积累方向再往前多走一段，这时梯度中指向积累梯度相反方向的量变大了许多，所以最终两个方向的梯度抵消，反而使得抖动方向的量迅速减少。Nesterov的衰减速度确实比动量方法要快不少。

很多科研人员已经给出了动量打折率的建议配置 – 0.9：
如果用$G$表示每一轮迭代的动量，$g$表示当前一轮迭代的更新量（方向 * 步长），$t$表示迭代轮数，$r$表示动量的打折率，那么对于时刻$t$的梯度更新量如下：
$G_t=r{G}_{t-1}+g_t$
$G_t=r(r{G}_{t-2}+g_{t-1})+g_t$
$G_t=r^2{G}_{t-2}+r{g}_{t-1}+g_t$
$G_t=r^2{g}_o+r^{t-1}{g}_1+...+g_t$
所以，对于第一轮迭代的更新$g_0$来说，从$G_0$到$G_T$，它的总贡献量为：
$$(r^t+r^{t-1}+...+r+1)g_0$$
它的贡献和为一个等比数列的和，比值为$r$。如果$r=0.9$，那么更新量在极限状态下贡献值：
$$\frac{g_0}{1-r}$$
当$r=0.9$时，它一共贡献了相当于自身10倍的能量。如果$r=0.99$，那就是100倍能量了。

SGD的变种算法

Adagrad

Adagrad是一种自适应的梯度下降方法。何为自适应呢？在梯度下降法中，参数的更新量等于梯度乘以学习率，也就是说，更新量和梯度是正相关的；而在实际应用中，每个参数的梯度各有不同，有的梯度大，有的梯度比较小，那么就有可能遇到参数优化不均衡的情况。

参数优化不均衡对模型训练来说不是件好事，这意味着不同的参数更新适用于不同的学习率。而Adagrad的自适应算法也正是要解决这个问题。算法希望不同参数的更新量能够比较均衡。对于已经更新比较多的参数，它的更新量要适当衰减，而更新比较少的参数，它的更新量要尽量多一些，它的参数更新公式如下：
$$x-=lr\cdot \frac{g}{\sqrt{\sum g^2}+\epsilon }$$

其中$\epsilon$的取值一般比较小，它只是为了防止分母为0.

def adagrad(x_start, step, g , delta=1e-8):
	x.np.array(x_start,dtype='float64')
	sum_grad = np.zeros_like(x)
	for i in range(50):
		grad = g(x)
		sum_grad += grad * grad
		x -= step* grad /(np.sqrt(sum_grad)+delta)
		if abs(sum(grad)) < 1e-6:
			break;
	return x

从公式和代码中可以发现，算法积累了历史的梯度值的和，并用这个加和来调整每个参数的更新量——对于之前更新量大的参数，分母也会比较大，于是未来它的更新量会比较小；对于之前更新量小的参数，分母也相对小一些，于是未来它的更新量会相对大一些。

Rmsprop

Adagrad算法有一个很大的问题，那就是随着优化的迭代次数不断增加， 更新公式的分母项会变得越来越大。所以理论上更新量也会越来越小，这对优化十分不利。Rmsprop就试图解决这个问题，在它的算法中，分母的梯度平方和不再随优化而增加，而是做加权平均。更新公式如下：
$$G_{t+1}=\beta G_t+(1-\beta )g^2$$
$$x_{t+1}=x_t-lr\frac{g}{{\sqrt{G}}_{t+1}+\epsilon }$$

def rmsprop(x_start, step, g , rms_decay = 0.9, delta=1e-8):
	x = np.array(x_start, dtype = 'float64')
	sum_grad = np.zeros_like(x)
	passing_dot = [x.copy()]
	for i in range(50):
		grad = g(x)
		sum_grad = rms_decay * sum_grad + (1 - rms_decay) * grad * grad
		x -= step * grad /(np.sqrt(sum_grad) + delta)
		passing_dot.append(x.copy())
		if abs(sum(grad)) < 1e-6:
			break;
	return x , passing_dot

Adam

Adam既包含了动量算法的思想，也包含了RmsProp的自适应梯度的思想。在计算过程汇总，Adam既要像动量算法那样计算累计的动量：
$$m_{t+1}={\beta }_1{m}_t+(1-{\beta }_1)g_t$$
又要像RmsProp那样计算梯度的滑动平方和：
$$v_{t+1}={\beta }_2{v}_t+(1-{\beta }_2)g_t^2$$
作者没有直接把这两个计算值加入最终计算的公式中，作者推导了两个计算量与期望的差距，于是给这两个变量加上了修正量：
$${\hat{m}}_t=\frac{m_t}{1-{\beta }_1^t}$$
$${\hat{v}}_t=\frac{v_t}{1-{\beta }_2^t}$$
最后，两个计算量将融合到最终的公式中：
$$x_{t+1}=x_t-lr\frac{{\hat{m}}_t}{\sqrt{{\hat{v}}_t}+\epsilon }$$

def adam(x_start, step, g, beta1 = 0.9, beta2 = 0.999, delta=1e-8):
	x = np.array(x_start, dtype='float64')
	sum_m = np.zeros_like(x)
	sum_v = np.zeros_like(x)
	passing_dot = [x.copy()]
	for i in range (50):
		grad = g(x)
		sum_m = beta1 * sum_m + (1 - beta1) * grad
		sum_v = beta2 * sum_v + (1 - beta2) * grad * grad
		correction = np.sqrt(1 - beta2 ** i) / (1 - beta1 ** i)
		x -= step * correction * sum_m / (np.sqrt(sum_v) + delta)
		passing_dot.append(x.copy())
		if abs(sum(grad)) < 1e-6:
			break;
	return x, passing_dot

Adam算法结合了动量的“惯性”和自适应的“起步快”这两个特点。综合来看，RmsProp和Adam的表现更平稳，现在大部分科研人员都在使用这两种优化方法。