凸函数性质：
- 它通常有唯一的极值点。
- Hessian矩阵正定。
求解非凸函数： \(f(x,y)=x^2-2x+100xy+10y^2+20y\)
Gradient Descent
\[\theta = \theta - \eta \cdot \nabla_\theta J( \theta)\]def batch_gd(x_start, step, g):
x = np.array(x_start, dtype='float64')
for i in range(50):
grad = g(x)
x -= grad * step
if abs(sum(grad)) < 1e-6:
break;
return x
GD算法的挑战：
- 难以挑选合适的学习率 $\eta$，太小收敛慢，太大容易振荡，甚至不收敛
- 期望学习率可以随着迭代次数变化
- 学习率对所有对参数是固定的，如果数据稀疏或特征具备不同的频率，会不合理
- 对非凸的损失函数（如NN），GD可能被困在鞍点（saddle point）
Mementum
动量算法把上一次迭代的梯度以一定权重（即 $\gamma$ 通常取0.9）叠加到本轮的梯度上，这样可以加速相应方向的SGD，避免振荡。
\[v_t = \gamma v_{t-1} + \eta \nabla_\theta J( \theta) \\ \theta = \theta - v_t\]def momentum(x_start, step, g, discount = 0.9):
x = np.array(x_start, dtype='float64')
passing_dot = [x.copy()]
for i in range(50):
grad = g(x)
pre_grad = pre_grad * discount + grad
x -= pre_grad * step
if abs(sum(grad)) < 1e-6:
break;
return x
Nesterov accelerated gradient
动量方法首先计算本轮的梯度，然后沿着更新后的累积梯度方向前进一大步。 NAG首先在上轮梯度方向前进一大步，然后计算梯度再做微调，这样避免跑得太快，反应过激。
\[v_t = \gamma v_{t-1} + \eta \nabla_\theta J( \theta - \gamma v_{t-1} ) \\ \theta = \theta - v_t\]可以看出Mementum其实并没有改变本轮梯度，而Nesterov利用之前的动量校正了本轮的梯度。
def nesterov(x_start, step, g, discount = 0.7):
x = np.array(x_start, dtype='float64')
passing_dot = [x.copy()]
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
if abs(sum(grad)) < 1e-6:
break;
return x
Adagrad
动量和NAG方法都是在梯度上做文章，还没有涉及到学习率，即对所有参数学习率都一样。
Adagrad可以让学习率自适应参数，对低频的参数加大更新，对高频的参数减小更新，因此它非常适合处理稀疏数据。
定义 $g_{t,i}$ 为参数 $\theta_i$ 第 $t$ 步的梯度：
\[g_{t, i} = \nabla_\theta J( \theta_i )\]则，SGD对每个参数 $\theta_i$ 在第 $t$ 步的更新公式（update rule）为：
\[\theta_{t+1, i} = \theta_{t, i} - \eta \cdot g_{t, i}\]而Adagrad调整相应的学习率 $\eta$，变成这样：
\[\theta_{t+1, i} = \theta_{t, i} - \dfrac{\eta}{\sqrt{G_{t, i} + \epsilon}} \cdot g_{t, i}\]整体用矩阵表示，$G_{t} \in \mathbb{R}^{d \times d}$ 是一个对角矩阵，对角线上的元素是参数 $\theta_i$ 从开始到第 $t$ 步的所有梯度的平方和，即 $G_{t,i}=\sum_t g_{t,i}^2$；$\epsilon$ 是个平滑因子，避免分母为0。
\[\theta_{t+1} = \theta_{t} - \dfrac{\eta}{\sqrt{G_{t} + \epsilon}} \odot g_{t}\]因为，$G_{t}$积累了参数过去所有梯度的平方和，所以如果初期某参数的梯度变化很小，则会加大该参数的学习率；但也由于 $G_{t}$ 是始终递增的，中后期学习率会越来越小，这有可能导致迭代提前结束（因为每轮之间的参数基本没有变化了）。
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
Adadelta
为了解决adagrad的学习率单调递减问题，adadelta并不累加过去所有梯度平方 $g^2$，而是约束在一个窗口 $w$ 里。
但直接保存过去 $w$ 步的所有 $g^2$ 并不是个高效的方法（尤其对于DNN这种大模型），adadelta的办法是计算过去 $g^2$ 的期望，并和当前梯度做一个折中（类似动量方法的思路，这里 $\gamma$ 一般也取 0.9）：
\[E[g^2]_t = \gamma E[g^2]_{t-1} + (1 - \gamma) g^2_t\]这里直接用 $E[g^2]_t$ 替代adagrad的 $G_t$ 就得到了 RMSprop 算法：
\[\theta_{t+1} = \theta_{t} - \dfrac{\eta}{\sqrt{E[g^2]_t + \epsilon}} \odot g_{t}\]进一步推导
把SGD参数更新部分用 $\Delta \theta$ 表示。
adadelta作者发现 $\Delta \theta$ 的量纲不对（包括之前的SGD、动量、adagrad），至少应该包含和参数 $\theta$ 同样的单位。类似 $E[g^2]_t$ 的方法 ，作者定义了参数增量的平方：
\[E[\Delta \theta^2]_t = \gamma E[\Delta \theta^2]_{t-1} + (1 - \gamma) \Delta \theta^2_t\]然后，然后把它开方再添加到 $\Delta \theta$ 的分子上，这样量纲就对了。而平方和开方其实就是计算RMS：
\[RMS[\Delta \theta]_{t} = \sqrt{E[\Delta \theta^2]_t + \epsilon}\]这样就得到了adadelta算法：
\[\Delta \theta_t = - \dfrac{RMS[\Delta \theta]_{t-1}}{RMS[g]_{t}} \odot g_{t} \\ \theta_{t+1} = \theta_t + \Delta \theta_t\]def adadelta(x_start, step, g, momentum = 0.9, delta=1e-1):
x = np.array(x_start, dtype='float64')
sum_grad = np.zeros_like(x)
sum_diff = np.zeros_like(x)
for i in range(50):
grad = g(x)
sum_grad = momentum * sum_grad + (1 - momentum) * grad * grad
diff = np.sqrt((sum_diff + delta) / (sum_grad + delta)) * grad
x -= step * diff
sum_diff = momentum * sum_diff + (1 - momentum) * (diff * diff)
if abs(sum(grad)) < 1e-6:
break;
return x
RMSprop
RMSprop是adadelta的特例（$\gamma=0.5$），不赘述。
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)
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)
if abs(sum(grad)) < 1e-6:
break;
return x
Adam
Adam（Adaptive Moment Estimation）是另一种自适应学习率算法。它综合adadelta和动量算法，既使用历史梯度平方和，也使用历史梯度和：
\[m_t = \beta_1 m_{t-1} + (1 - \beta_1) g_t \\ v_t = \beta_2 v_{t-1} + (1 - \beta_2) g_t^2\]这里 $m_t$ 和 $v_t$ 对应均值和方差，其实是梯度的估计值，这也是adam算法命名的原因。由于 $m_t$ 和 $v_t$ 均初始化为 0 向量，这里做了一下平滑：
\[\hat{m}_t = \dfrac{m_t}{1 - \beta^t_1} \\ \hat{v}_t = \dfrac{v_t}{1 - \beta^t_2}\]这样得到最终的adam算法：
\[\theta_{t+1} = \theta_{t} - \dfrac{\eta}{\sqrt{\hat{v}_t} + \epsilon} \hat{m}_t\]其中，一般 $\beta_1$ 取值0.9， $\beta_2$ 取值 0.999， $\epsilon$ 取值 1e-8。
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)
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 * sum_m / (np.sqrt(sum_v + delta))
if abs(sum(grad)) < 1e-6:
break;
return x
参考
- http://sebastianruder.com/optimizing-gradient-descent/
- http://blog.csdn.net/luo123n/article/details/48239963
- https://github.com/hsmyy/zhihuzhuanlan/blob/master/momentum.ipynb