无约束最优化问题求解

Posted on

无约束最优化问题求解方法的学习笔记

神经网络中的学习过程可以形式化为最小化损失函数问题, 该损失函数一般是由训练误差和正则项组成

L(w)=train_err(w)+norm(w)L (w)=train\_err(w) + norm(w)

损失函数的一阶偏导为

∇iL(w)=dfdwi,i=1,2,...,n\nabla_i L(w)=\frac {df}{dw_i},i=1,2,...,n

损失函数二阶偏导可以使用海塞矩阵 Hessian Matrix H\mathbf{H} 表示, 其中每个权重向量 ii 的元素 jj 的二阶偏导数为

Hi,jL(w)=d2fdwi⋅dwj\mathbf{H}_{i,j} L(w) = \frac {d^2 f} {dw_i \cdot dw_j}

一阶求解方法有 SGD Adam RMSProp 等，利用梯度（超平面）的信息求解，计算高效，收敛稍慢，需要超参数。

二阶求解方法有牛顿法，拟牛顿法，BFGS，L-BFGS 等，用二阶梯度（超曲面）的信息求解，计算复杂，收敛快，不需要超参数。

针对 LR + L1，求解方法有 OWL-QN （改进自 L-BFGS）/ grafting / GIS BLMVM / IRIS-LARS / PGD（Proximal Gradient Decent 近端梯度下降）/ ADMM

牛顿法

用损失函数的二阶偏导数寻找更好的训练方向. 记

Li=L(wi),gi=∇i,Hi=Hf(wi)Li = L(w_i), g_i=\nabla_i, \mathbf{H}_i=\mathbf{H}f(w_i)

在 w0w_0 点使用泰勒级数展开式二次逼近损失函数 LL

L=L0+g0⋅(w−w0)+12⋅(w−w0)2⋅H0L=L_0+g_0 \cdot (w-w_0) + \frac{1}{2} \cdot (w-w_0)^2 \cdot \mathbf{H}_0

另一阶导数 g=0g=0 , 有

g=g0+H0⋅(w−w0)=0g=g_0+\mathbf{H}_0 \cdot (w-w_0)=0

参数ww 迭代公式为

wi+1=wi−Hi−1⋅giw_{i+1}=w_i - {\mathbf{H}_i^{-1}} \cdot {g_i}

如果海塞矩阵非正定, 损失函数并不能保证在每一次迭代中都是减少的. 可以通过加上学习速率解决这个问题

wi+1=wi−Hi−1⋅gi⋅ηiw_{i+1}=w_i - {\mathbf{H}_i^{-1}} \cdot {g_i} \cdot \eta_i

优点: 比一阶导数更少迭代

缺点: 计算复杂度比一阶导数更高，约O(n3)O(n^3)，因为对海塞矩阵及其逆的精确求值在计算量复杂度是十分巨大的；海森矩阵非正定，容易失效。

拟牛顿法

Quasi-Newton method

拟牛顿法不直接计算海塞矩阵然后求其矩阵的逆, 而是在每次迭代的时候, 利用一阶偏导矩阵 Jacobian Matrix 或其他方法, 以逼近 Hessian Matrix 的逆, 复杂度约O(n2)O(n^2).

  1. 输入:

L,w0,H0−1,QuasiUpdateL,w_0,\mathbf{H}_0^{-1},\mathbf{QuasiUpdate}

  1. 对于n=0,1,...n = 0,1,... 更新参数, 直到收敛:

a) 计算搜索方向和步长

d=H−1gnα←minα≥0L(wn−αd)wn+1←wn−αd\begin{aligned} \\ d &= \mathbf{H}^{-1} g_n\\ \alpha &\leftarrow \min_{\alpha \geq 0} L(w_n - \alpha d) \\ w_{n+1} &\leftarrow w_n - \alpha d \\ \end{aligned}

b) 计算并保存

gn+1←∇L(wn+1)sn+1←wn+1−wnLn+1←gn+1−gn\begin{aligned}\\ g_{n+1} &\leftarrow \nabla L(w_{n+1}) \\ s_{n+1} &\leftarrow w_{n+1} - w_{n} \\ L_{n+1} &\leftarrow g_{n+1} - g_{n} \\ \end{aligned}

c) 更新

Hn+1−1←QuasiUpdate(Hn−1,sn+1,Ln+1)\mathbf{H}_{n+1}^{-1} \leftarrow \mathbf{QuasiUpdate}(\mathbf{H}_n^{-1}, s_{n+1}, L_{n+1})

DFP

BFGS

是一种拟牛顿法.

记

gn+1=∇L(wn+1)yn+1=gn+1−gnsn+1=wn+1−wnB≈HD≈H−1\begin{aligned} g_{n+1} &= \nabla L(w_{n+1}) \\ y_{n+1} &= g_{n+1} - g_{n} \\ s_{n+1} &= w_{n+1} - w_{n} \\ \mathbf{B} &\approx \mathbf{H} \\ \mathbf{D} &\approx \mathbf{H^{-1}} \end{aligned}

由中值定理我们知道

gn−gn−1=Bn(wn−wn−1)g_n - g_{n-1} = \mathbf{B}_n(w_n - w_{n-1})

所以有

Bn−1yn=sn\mathbf{B}_n^{-1} y_n = s_n

同时，Hessian矩阵是对称矩阵。在这两个条件的基础上，我们希望Bn\mathbf{B}_n 相对于\mathbf{B}_{n−1} 的变化并不大, 即

minB−1 ∣∣Bn−1−Bn−1−1∣∣s.t.B−1yn=snB−1is symmetric\begin{aligned} \min_{\mathbf{B}^{-1}}\ &||\mathbf{B}^{-1}_n - \mathbf{B}^{-1}_{n-1}||\\ \text{s.t.} & \mathbf{B}^{-1}y_n=s_n\\ & \mathbf{B}^{-1} \text{is symmetric}\\ \end{aligned}

其中范式为 Weighted Frobenius Norm, 这个式子的解, 即 BFGS 的QuasiUpdate\mathbf{QuasiUpdate} 定义为

Bn+1−1=(I−pnynsnT)Bn−1(I−pnsnynT)+pnsnsnT\mathbf{B}^{-1}_{n+1}=(I-p_n y_n s_n^T) \mathbf{B}^{-1}_n (I - p_n s_n y_n^T)+p_n s_n s_n^T

其中

pn=(ynTsn)−1p_n = (y_n^T s_n)^{-1}

其中，只要 Hn−1\mathbf{H}^{-1}_{n} 是正定的，Hn+1−1\mathbf{H}^{-1}_{n+1} 就一定是正定的，所以需选择一个正定的 H0−1\mathbf{H}^{-1}_{0} 或单位矩阵。

完整算法 wiki：

  1. 给定初始值 w_0 和精度阈值 ϵ\epsilon，另 D0=I,n=0\mathbf{D_0}=I, n=0
  2. 确定搜索方向 dn=−Dngnd_n=-\mathbf{D_n} g_n
  3. 通过 line search 方法得到步长 λn\lambda_n，令$ s_n = \lambda_n d_n$
  4. 参数更新 wn+1=wn+snw_{n+1} = w_n + s_n
  5. 若 ∣∣gn∣∣≤ϵ||g_n|| \le \epsilon 则算法结束
  6. 计算 yn=gn+1−gny_n = g_{n+1} - g_n
  7. 计算 Dn+1=(I−pnsnynT)Dn(I−pnynsnT)+pnsnsnT\mathbf{D_{n+1}}=\big(I - p_n s_n y_n^T\big) \mathbf{D_n} \big(I - p_n y_n s_n^T\big) + p_n s_n s_n^T
  8. 令 n=n+1n=n+1，执行 2

L-BFGS

理解L-BFGS算法

BFGS 每次迭代的复杂度 O(n2)O(n^2)。而L-BFGS, 即 limited-mem BFGS, 是 O(nm),m=5→40O(nm), m=5 \to 40。

BFGS 需要保存每次迭代的 sn,yns_n, y_n 值。而 L-BFGS 只使用最近的 mm 个 sn,yns_n,y_n 用于计算，减轻内存负担

在 full batch 非常有效，但是在 mini batch 收敛较慢。在大型数据集常用 SGD Adam 或 FTLR 等在线算法求解。

img

img

其中 line search 方法有 backtracking line search wiki 等。

OWL-QN

为解决 LR + L1 的正则项在 0 点不可导的问题，OWL-QN 增加以下几点，改进了 L-BFGS：

  1. 次梯度 ⋄f\diamond f 替代梯度，解决如 L1 在零点不连续、不可微的问题

  2. 限定搜索方向象限：pk⇐π(dk;vk)p^k \Leftarrow \pi(d^k;v^k) ，即不跨越象限（Orthant-Wise）

  3. YY 向量（loss 梯度 delta，不包括正则项）

img

PGD

Projected Gradient Descent

目标函数

minwf(w)+∣∣w∣∣22+∑i=1dλi∣∣wi∣∣1\min_w f(w) + ||\mathbf{w}||_2^2 + \sum_{i=1}^d\lambda_i||w_i||_1

每次迭代更新

wik+1=Proxλi,η(wik−μwik−η∇fi(wk))w^{k+1}_i = Prox_{\lambda_i,\eta}(w^k_i - \mu w_i^k - \eta \nabla f_i(\mathbf{w}^k))

其中

Proxλi,η(z)={z−λiη, if z−λiη>0z+λiη, if z−λiη<00, otherwise\begin{aligned} Prox_{\lambda_i,\eta}(z) = \begin{cases} z - \lambda_i \eta, \text{ if } z - \lambda_i \eta > 0 \\ z + \lambda_i \eta, \text{ if } z - \lambda_i \eta < 0 \\ 0, \text{ otherwise} \end{cases} \end{aligned}

共轭梯度法

Conjugate gradient, 可认为是梯度下降法和牛顿法的中间物, 希望能加速梯度下降的收敛速度, 同时避免使用海塞矩阵进行求值、储存和求逆获得必要的优化信息. 每次迭代, 沿着共轭方向 (conjugate directions) 执行搜索的, 所以通常该算法要比沿着梯度下降方向优化收敛得更迅速. 共轭梯度法的训练方向是与海塞矩阵共轭的.

TODO

梯度下降

wi+1:=wi−di⋅ηiw_{i+1} := w_i - d_i \cdot \eta_i

优点: 使用一阶导数计算, 复杂度小于二阶导数

缺点: 变量没有归一化, 锯齿下降现象, 因为非线性函数局部的梯度方向并不一定就是朝着最优点

SGD

Stochastic Gradient Descent

每次迭代, 选取部分样本进行计算

相对于梯度下降，loss 函数更加波动，能帮助函数跳入另一个局部最优解。

Momentum

An overview of gradient descent optimization algorithms

为解决 SGD 在沟壑（有一维梯度值特别大）的 Z 字形游走问题，引入动量，减少 Z 字震荡

vt=γvt−1+ηgtθ=θ−vtγ≈0.9\begin{aligned} v_t &= \gamma v_{t-1} + \eta g_t \\ \theta &= \theta - v_t \\ \gamma &\approx 0.9 \end{aligned}

Nesterov Momentum

NAG ref

NAG 是一种能给予梯度项「预测」功能的方法。我们通过计算θ−γvt−1\theta - \gamma v_{t-1} 给我们一个对下个参数θ\theta 的估计，加速收敛。

NAG 的更新规则是

vt=γvt−1+ηg(θ−γvt−1)θ=θ−vt\begin{aligned} v_t &= \gamma v_{t-1} + \eta g(\theta - \gamma v_{t-1}) \\ \theta &= \theta - v_t \\ \end{aligned}

NAG 比 Momentum 收敛得更快。

AdaGrad

使用了自适应技术来更新学习率：对变化[小|大]的参数进行更[大|小]的更新（通过 g）。epsilon 是一个用于抑制学习率产生变动率的常量。ref

比如 Dense / Sparse 特征共存的情况会有所优化。

Dean [3] 发现它改进 SGD 的鲁棒性，将其应用在大规模神经网络训练。

gt,i=ΔθJ(θi)θt+1,i=θt,i−η⋅gt,i(SGD)θt+1,i=θt,i−ηGt,ii+ϵ⋅gt,i(Adagrad)θt+1=θt−ηGt+ϵ⊙gt\begin{aligned} g_{t,i} &= \Delta_\theta J(\theta_i) \\ \theta_{t+1, i} &= \theta_{t,i} - \eta \cdot g_{t,i} (SGD) \\ \theta_{t+1,i} &= \theta_{t,i} - \frac{\eta}{\sqrt{G_{t,ii} + \epsilon}}\cdot g_{t,i} (Adagrad) \\ \theta_{t+1} &= \theta_{t} - \frac{\eta}{\sqrt{G_t + \epsilon}} \odot g_t \\ \end{aligned}

其中 Gt∈Rd×dG_t \in \mathbb{R}^{d \times d} 是对角矩阵，元素i,ii,i 是 θi\theta_i 从t0t^{0} 到 tit^{i} 的平方和

代码表示

1
2
cache += dx ** 2
x += - learning_rate * dx / (np.sqrt(cache) + 1e-7)

Adagrad 有个缺点：随着迭代次数的增多，学习率项ηGt,ii+ϵ\frac{\eta}{\sqrt{G_{t,ii} + \epsilon}} 会急剧递减，可能会掉入局部最优点。Adadelta 和 RMSprop 尝试解决这个问题。

Adadelta

是 Adagrad 的扩展，减少 Adagrad 快速下降的学习率。把 Adagrad 的梯度平方和GtG_t 限制在时间窗口内

Δθt=−ηE[g2]t+ϵ∗gtE[g2]t=γE[g2]t−1+(1−γ)gt2E[Δθ2]t=γE[Δθ2]t−1+(1−γ)Δθt2RMS[Δθ]t=E[Δθ2]t+ϵΔθt=−RMS[Δθ]t−1RMS[g]tgtθt+1=θt+Δθt\begin{aligned} \Delta\theta_t &= - \frac{\eta}{\sqrt{ E[g^2]_t + \epsilon}} * g_t \\ E[g^2]_{t} &= \gamma E[g^2]_{t-1} + (1-\gamma) g^2_t \\ E[\Delta\theta^2]_t &= \gamma E[\Delta\theta^2]_{t-1} + (1-\gamma) \Delta\theta_t^2 \\ RMS[\Delta\theta]_t &= \sqrt{E[\Delta\theta^2]_t + \epsilon} \\ \Delta\theta_t &= - \frac{RMS[\Delta\theta]_{t-1}}{RMS[g]_t} g_t \\ \theta_{t+1} &=\theta_{t} + \Delta\theta_t \\ \end{aligned}

RMSprop

由 Hinton 在 coursera 提出 [2]。类似 Adadelta，解决 Adagrad 快速降低的学习率

E[g2]t=0.9E[g2]t−1+0.1gt2θt+1=θt−ηE[g2]t+ϵgt\begin{aligned} E[g^2]_t &= 0.9 E[g^2]_{t-1} + 0.1 g_t^2 \\ \theta_{t+1} &= \theta_t - \frac{\eta}{\sqrt{E[g^2]_t} + \epsilon} g_t \end{aligned}

Hinton 建议 γ=0.9,η=0.001\gamma = 0.9, \eta=0.001

代码

1
2
cache = decay_rate * cache + (1 - decay_rate) * dx**2
x += - learning_rate * dx / (np.sqrt(cache) + 1e-7)

Adam

Adaptive Moment Estimation, 除了像 Adadelta RMSprop 存储之前迭代的快速下降的梯度平方信息vtv_t ，它还存储之前迭代的快速下降的梯度信息mtm_t ，类似 momentum

mt=β1mt−1+(1−β1)gtvt=β2vt−1+(1−β2)gt2m^t=mt1−β1tv^t=vt1−β2tθt+1=θt−ηv^t+ϵ\begin{aligned} 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 \\ \hat m_t &= \frac{m_t}{1-\beta_1^t} \\ \hat v_t &= \frac{v_t}{1-\beta_2^t} \\ \theta_{t+1} &= \theta_t - \frac{\eta}{\sqrt{\hat v_t + \epsilon}} \\ \end{aligned}

其中，[ mtm_t | vtv_t ] 是第 [1|2] 个时刻的梯度估计 [the mean | uncenter variance]，也是 Adam 名字的由来

1
2
3
4
5
6
7
8
m = beta1 * m + (1 - beta1) * dx
v = beta2 * v + (1 - beta2) * dx**2
a
# bias correlation, only relavant in first few iter when m,v close to 0
m /= 1 - beta1**t
v /= 1 - beta2**t

x += - learning_rate * m / (np.sqrt(v) + 1e-7)

SGD 技巧

  • shuffling training data
  • curriculumn learning Bengio09 循序渐进学习
  • batch normalization
  • early stopping
  • gradient noise Neelakantan

gt,i=gt,i+N(0,σt2)σt2=η(1+t)γ\begin{aligned} g_{t,i} &= g_{t,i} + N(0,\sigma_t^2) \\ \sigma_t^2 &= \frac{\eta}{(1+t)^\gamma} \end{aligned}

总结

opt-update-rule

gif demo

SGD, Momentum, NAG 花更多时间跳出对称区域

如何选择？

  • 如果数据稀疏，或者训练深度神经网络，用 adative learning-rate methods，不需要手动调整学习率，超参数默认值 (0.9) 很容易达到较好的效果
  • Adam 可能是较好的选择，Kingma 指出 Adam 的 bias-correction 的特性帮助 Adam 稍微优于 RMSProp Adadelta
  • 有趣的是，很多最新的论文，都直接使用了（不带动量项的）Vanilla SGD 法，配合一个简单的学习率（退火）列表。这些 SGD 最终都能帮助他们找到一个最小值，但会花费远多于上述方法的时间。

Learning Rate Decay

step decay

αt=βαt−1\alpha_t = \beta \alpha_{t-1}

exp decay

α=α0∗e−kt\alpha = \alpha_0 * e^{-kt}

1/t decay

a=a0/(1+kt)a = a_0/(1+kt)

Dropout

由 Hinton 等人在 [4] 提出。做法是在 forward pass / backward prop 随机地把部分神经元置 0。

dropout 能使网络

  • redundant representation
  • ensemble of models
  • 提升成绩

Weight Init

参考 [0]

在小规模网络，可以这样初始化 W，但可能导致 non-homogeneous distribution of activations across the layers of a network

1
W = 0.01 * np.random.randn(D, H)

Pre-train

Bengio 在 2009 提出的 [9] ，用 unsupervised 方法 pre-training 网络，如 stacked-auto encoder，达到初始化网络的目的。

Xavier Initialization

Glorot 等人在 2010 发表了 [5] 探讨了随机初始化参数来训练 SGD 效果不好的原因，并提出 Xavier Initialization 的参数初始化方法

1
W = np.random.randn(fan_in, fan_out) / np.sqrt(fan_in)

可用于 tanh 网络，但不适用于 ReLU

He et al., 2015

[7] 改进的 Xavier Initialization 方法，可用于 ReLU

1
W = np.random.randn(fan_in, fan_out) / np.sqrt(fan_in / 2)

Batch Normalization

Ioffe et al. 2015 [6] 等人提出，对每个 mini batch，归一化每个样本，作为网络的输入

\begin{align} \hat x^{(k)} = \frac{x^{(k)} - E[x^{(k)}]}{\sqrt{Var[x^{(k)}]}} \\ y^{(k)} = \gamma^{(k)} \hat x^{(k)} + \beta^{(k)} \end{align}

优点

  • 改善梯度流
  • 允许更高的学习率
  • 降低对 weight init 的依赖
  • 有 regulization 的作用，可以减少 dropout

Data-dependent Initialization

[8] 将 pre-train 时间减少 3 个数量级。

Appendix

https://nlperic.github.io/line-search/

Reference

[0] CS231n Winter 2016: Lecture 5: Neural Networks Part 2 VIDEO PDF

[1] CS231n Winter 2016 Lecture 6 Neural Networks Part 3 VIDEO PDF

[2] T. Tieleman and G. E. Hinton. Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude. Coursera. 2012. PDF

[3] Dean, Jeffrey, et al. “Large scale distributed deep networks.” Advances in neural information processing systems. 2012. PDF

[4] Srivastava N, Hinton G E, Krizhevsky A, et al. Dropout: a simple way to prevent neural networks from overfitting[J]. Journal of machine learning research, 2014, 15(1): 1929-1958. PDF

[5] Glorot X, Bengio Y. Understanding the difficulty of training deep feedforward neural networks[C]//Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics. 2010: 249-256. PDF

[6] Ioffe S, Szegedy C. Batch normalization: Accelerating deep network training by reducing internal covariate shift[C]//International Conference on Machine Learning. 2015: 448-456. PDF

[7] He K, Zhang X, Ren S, et al. Delving deep into rectifiers: Surpassing human-level performance on imagenet classification[C]//Proceedings of the IEEE international conference on computer vision. 2015: 1026-1034. PDF

[8] Krähenbühl P, Doersch C, Donahue J, et al. Data-dependent initializations of convolutional neural networks[J]. arXiv preprint arXiv:1511.06856, 2015. PDF

[9] Erhan D, Manzagol P A, Bengio Y, et al. The difficulty of training deep architectures and the effect of unsupervised pre-training[C]//Artificial Intelligence and Statistics. 2009: 153-160. PDF

本文有帮助？