Reference：ADADELTA: An Adaptive Learning Rate Method
超参数
超参数（Hyper-Parameter)是困扰神经网络训练的问题之一，因为这些参数不可通过常规方法学习获得。
神经网络经典五大超参数:
学习率(Leraning Rate)、权值初始化(Weight Initialization)、网络层数(Layers)
单层神经元数(Units)、正则惩罚项（Regularizer|Normalization)
这五大超参数使得神经网络更像是一门实践课，而不是理论课。
懂神经网络可能只要一小时，但是调神经网络可能要几天。
因此，后来Vapnik做SVM支持向量机的时候，通过巧妙的变换目标函数，避免传统神经网络的大部分超参数，
尤其是以自适应型的支持向量替代人工设置神经元，这使得SVM可以有效免于过拟合之灾。
传统对抗这些超参数的方法是经验规则（Rules of Thumb)。
这几年，随着深度学习的推进，全球神经网络研究者人数剧增，已经有大量研究组着手超参数优化问题：
★深度学习先锋的RBM就利用Pre-Traning自适应调出合适的权值初始化值。
★上个世纪末的LSTM长短期记忆网络，可视为“神经网络嵌套神经网络”，自适应动态优化层数。
★2010年Duchi et.al 则推出AdaGrad，自适应来调整学习率。
自适应调整学习率的方法，目前研究火热。一个经典之作，是 Matthew D. Zeiler 2012年在Google实习时，
提出的AdaDelta。
Matthew D. Zeiler亦是Hinton的亲传弟子之一，还是商业天才，大二时办了一个公司卖复习旧书。
Phd毕业之后，创办了Clarifai，估值五百万刀。参考[知乎专栏]
Clarifai的杰出成就是赢得了ImageNet 2013冠军，后来公布出CNN结构的时候，Caffe、Torch之类
的框架都仿真不出他在比赛时候跑的结果，应该是用了不少未公布的黑科技的。
再看他2012年提出的AdaDelta，肯定是用在的2013年的比赛当中，所以后来以普通方式才无法仿真的。
梯度更新
2.1 [一阶方法] 随机梯度
SGD(Stochastic Gradient Descent)是相对于BGD(Batch Gradient Descent)而生的。
BGD要求每次正反向传播，计算所有Examples的Error，这在大数据情况下是不现实的。
最初的使用的SGD，每次正反向传播，只计算一个Example，串行太明显，硬件利用率不高。
后续SGD衍生出Mini-Batch Gradient Descent，每次大概推进100个Example，介于BGD和SGD之间。
现在，SGD通常是指Mini-Batch方法，而不是早期单Example的方法。
一次梯度更新，可视为：
$x_{t+1}=x_{t}+Delta x_{t} quad where quad Delta x_{t}=-eta cdot g_{t}$
$x$为参数，$t$为时序，$Delta$为更新量，$eta$为学习率，$g$为梯度
2.2 [二阶方法] 牛顿法
二阶牛顿法替换梯度更新量：
$Delta x_{t}=H_{t}^{-1} cdot g_{t}$
$H$为参数的二阶导矩阵，称为Hessian矩阵。
牛顿法，用Hessian矩阵替代人工设置的学习率，在梯度下降的时候，可以完美的找出下降方向，
不会陷入局部最小值当中，是理想的方法。
但是，求逆矩阵的时间复杂度近似$O(n^{3})$，计算代价太高，不适合大数据。
常规优化方法
3.1 启发式模拟退火
早期最常见的手段之一就是模拟退火。当然这和模拟退火算法没有半毛钱关系。
引入一个超参数(常数)的退火公式：
$eta_{t}=frac{eta _{0}}{1+d imes t}$
$eta _{0}$为初始学习率，$d$为衰减常数，通常为$10^{-3}$
模拟退火基于一个梯度法优化的事实：
在优化过程中，Weight逐渐变大，因而需要逐渐减小学习率，保证更新平稳。
3.2 动量法
中期以及现在最普及的就是引入动量因子：
$Delta x_{t}= ho Delta x_{t-1}-eta cdot g_{t}$
$ ho$为动量因子，通常设为0.9
在更新中引入0.9这样的不平衡因子，使得：
★在下降初期，使用前一次的大比重下降方向，加速。
★在越过函数谷面时，异常的学习率，会使得两次更新方向基本相反，在原地”震荡“
此时，动量因子使得更新幅度减小，协助越过函数谷面。
★在下降中后期，函数面局部最小值所在的吸引盆数量较多，一旦陷进吸引盆当中，
$Gradient ightarrow 0$，但是前后两次更新方向基本相同。
此时，动量因子使得更新幅度增大，协助跃出吸引盆。
3.3 AdaGrad
AdaGrad思路基本是借鉴L2 Regularizer，不过此时调节的不是$W$，而是$Gradient$:
$Delta x_{t}=-frac{eta }{sqrt{sum_{ au=1}^{t}(g_{ au})^{2}}}cdot g_{t}$
AdaGrad过程，是一个递推过程，每次从$ au=1$，推到$ au=t$，把沿路的$Gradient$的平方根，作为Regularizer。
分母作为Regularizer项的工作机制如下：
★训练前期，梯度较小，使得Regularizer项很大，放大梯度。[激励阶段]
★训练后期，梯度较大，使得Regularizer项很小，缩小梯度。[惩罚阶段]
另外，由于Regularizer是专门针对Gradient的，所以有利于解决Gradient Vanish/Expoloding问题。
所以在深度神经网络中使用会非常不错。
当然，AdaGrad本身有不少缺陷：
★初始化W影响初始化梯度，初始化W过大，会导致初始梯度被惩罚得很小。
此时可以人工加大$eta$的值，但过大的$eta$会使得Regularizer过于敏感，调节幅度很大。
★训练到中后期，递推路径上累加的梯度平方和越打越多，迅速使得$Gradinet$被惩罚逼近0，提前结束训练。
AdaDelta
AdaDelta基本思想是用一阶的方法，近似模拟二阶牛顿法。
4.1 矩阵对角线近似逆矩阵
1988年，[Becker&LeCun]提出一种用矩阵对角线元素来近似逆矩阵的方法：
$Delta x_{t}=-frac{1}{left | diag(H_{t}) ight |+mu }cdot g_{t}$
$diag$指的是构造Hessian矩阵的对角矩阵，$mu$是常数项，防止分母为0。
2012年，[Schaul&S. Zhang&LeCun]借鉴了AdaGrad的做法，提出了更精确的近似：
$Delta x_{t}=-frac{1}{left | diag(H_{t}) ight |}frac{E[g_{t}-w:t]^{2}}{E[g_{t}^{2}-w:t]}cdot g_{t}$
$E[g_{t}-w:t]$指的是从当前t开始的前w个梯度状态的期望值。
$E[g_{t}^{2}-w:t]$指的是从当前t开始的前w个梯度状态的平方的期望值。
同样是基于Gradient的Regularizer，不过只取最近的w个状态，这样不会让梯度被惩罚至0。
4.2 窗口和近似概率期望
计算$E[g_{t}-w:t]$，需要存储前w个状态，比较麻烦。
AdaDelta使用了类似动量因子的平均方法：
$E[g^{2}]_{t}= ho E[g^{2}]_{t-1}+(1- ho )g_{t}^{2}$
当$ ho=0.5$时，这个式子就变成了求梯度平方和的平均数。
如果再求根的话，就变成了RMS(均方根)：
$RMS[g]_{t}=sqrt{E[g^{2}]_{t}+epsilon }$
再把这个RMS作为Gradient的Regularizer：
$Delta x_{t}=-frac{eta}{RMS[g]_{t}}cdot g_{t}$
其中，$epsilon$是防止分母爆0的常数。
这样，就有了一个改进版的AdaGrad。
该方法即Tieleman&Hinton的RMSProp，由于RMSProp和AdaDelta是同年出现的，
Matthew D. Zeiler并不知道这种改进的AdaGrad被祖师爷命名了。
RMSProp利用了二阶信息做了Gradient优化，在BatchNorm之后，对其需求不是很大。
但是没有根本实现自适应的学习率，依然需要线性搜索初始学习率，然后对其逐数量级下降。
另外，RMSProp的学习率数值与MomentumSGD差别甚大，需要重新线性搜索初始值。
注：$epsilon$的建议取值为1，出处是Inception V3，不要参考V3的初始学习率。
4.3 Hessian方法与正确的更新单元
Zeiler用了两个反复近似的式子来说明，一阶方法到底在哪里输给了二阶方法。
首先，考虑SGD和动量法：
$Delta x propto gpropto frac{partial f}{partial x} propto frac{1}{x}$
$Delta x$可以正比到梯度$g$问题，再正比到一阶导数。而$log$一阶导又可正比于$frac{1}{x}$。
再考虑二阶导Hessian矩阵法：
这里为了对比观察，使用了[Becker&LeCun 1988]的近似方法，让求逆矩阵近似于求对角阵的倒数：
$Delta x propto H^{-1}gpropto frac{frac{partial f}{partial x}}{frac{partial^{2}f}{partial x^{2}}}propto frac{frac{1}{x}}{frac{1}{x}*frac{1}{x}}propto x$
$Delta x$可以正比到Hessian逆矩阵$H^{-1}cdot g$问题，再正比到二阶导数。而$log$二阶导又可正比于$x$。
可以看到，一阶方法最终正比于$frac{1}{x}$，即与参数逆相关：参数逐渐变大的时候，梯度反而成倍缩小。
而二阶方法最终正比于$x$，即与参数正相关：参数逐渐变大的时候，梯度不受影响。
因此，Zeiler称Hessian方法得到了Correct Units(正确的更新单元)。
4.4 由Hessian方法推导出一阶近似Hessian方法
基于[Becker&LeCun 1988]的近似方法，有：
$Delta x approx frac{frac{partial f}{partial x}}{frac{partial^{2}f}{partial x^{2}}}$
进而又有：
$frac{frac{partial f}{partial x}}{frac{partial^{2}f}{partial x^{2}}}=frac{1}{frac{partial^{2}f}{partial x^{2}}}cdot frac{partial f}{partial x}=frac{1}{frac{partial^{2}f}{partial x^{2}}}cdot g_{t}$
简单收束变形一下, 然后用RMS来近似：
$frac{1}{frac{partial^{2}f}{partial x^{2}}}=frac{Delta x}{frac{partial f}{partial x}}approx -frac{RMS[Delta x]_{t-1}}{RMS[g]_{t}}$
最后，一阶完整近似式：
$Delta x= -frac{RMS[Delta x]_{t-1}}{RMS[g]_{t}}cdot g_t$
值得注意的是，使用了$RMS[Delta x]_{t-1}$而不是$RMS[Delta x]_{t}$，因为此时$Delta x_{t}$还没算出来。
4.5 算法流程
$quadquadquadqquadqquadqquad ALGORITHM:ADADELTA\\\\Require:DecayRate \, ho \, ,Constant \,\,epsilon \Require:InitialParam \,\,x_{1} \1: quad Initialize\,\,accumulation \,\,variables \,\,E[g^{2}]_{0}=E[Delta x^{2}]_{0=0} \2: quad For \,\,t=1:T \,\, do \,\, Loop \,\, all \,\,updates \3: quad quad Compute \,\,Gradients:g_{t} \4: quad quad Accumulate \,\, Gradient:E[g^{2}]_{t}= ho E[g^{2}]_{t-1}+(1- ho )g_{t}^{2} \5: quad quad Compute \,\,Update:Delta x= -frac{RMS[Delta x]_{t-1}}{RMS[g]_{t}}cdot g_t \6: quad quad Accumulate \,\, Updates:E[Delta x^{2}]_{t}= ho E[Delta x^{2}]_{t-1}+(1- ho )Delta x^{2} \7: quad quad Apply \,\,Update:x_{t+1}=x_{t}+Delta x_{t} \8: quad End \,\,For$
4.6 Theano实现
论文中，给出的两个超参数的合适实验值。
$ ho=0.95 quadquad epsilon=1e-6$
Theano的实现在LSTM的教学部分，个人精简了一下：
def AdaDelta(tparams,grads): p=0.95;e=1e-6 # init delta_x2=[theano.shared(p.get_value() * floatX(0.)) for k, p in tparams.iteritems()] g2 = [theano.shared(p.get_value() * floatX(0.)) for k, p in tparams.iteritems()] # first to update g2 update_g2=[(g2, p * g2 + (1-p) * (g ** 2)) for g2, g in zip(g2, grads)] fn_update_1=theano.function(inputs=[],updates=update_g2) #calc delta_x by RMS delta_x=[-T.sqrt(delta_x2_last + e) / T.sqrt(g2_now + e) * g for g, delta_x2_last, g2_now in zip(grads,delta_x2,g2)] # then to update delta_x2 and param update_delta_x2=[(delta_x2, p * delta_x2 + (1-p) * (delta_x ** 2)) for delta_x2, delta_x in zip(delta_x2, delta_x)] update_param=[(param, param + delta) for param, delta in zip(tparams.values(), delta_x)] fn_update_2=theano.function(inputs=[],updates=update_delta_x2+update_param) #return the update function of theano return fn_update_1, fn_update_2
4.7 Dragon(Caffe)实现
默认代码以我的Dragon框架为准，对Caffe代码进行了重写。
// hpp文件 template <typename Dtype> class AdaDeltaSolver :public SGDSolver < Dtype > { public: AdaDeltaSolver(const SolverParameter& param) :SGDSolver<Dtype>(param) { } AdaDeltaSolver(const string& param_file) :SGDSolver<Dtype>(param_file) { } protected: virtual void computeUpdateValue(int param_id, Dtype rate); virtual void applyUpdate(); }; // cpp文件 #include "gradient_solver.hpp" template <typename Dtype> void AdaDeltaSolver<Dtype>::computeUpdateValue(int param_id, Dtype rate){ Blob<Dtype>* net_param = net->getLearnableParams()[param_id]; const Dtype lr_mult = net->getLrMults()[param_id]; Dtype eps = param.delta(); Dtype momntum = param.momentum(); // adadelta will ignore base_lr Dtype lr = lr_mult; const int count = net_param->count(); switch (Dragon::get_mode()){ case Dragon::CPU: // history store for E[g^2] // update store for E[delta^2] // history=momentum*history + (1-momentum)*(diff^2) // 1. compute diff^2 in temp dragon_powx<Dtype>(count, net_param->cpu_diff(), Dtype(2), temp[param_id]->mutable_cpu_data()); // 2. compute history dragon_cpu_axpby<Dtype>(count, Dtype(1) - momntum, temp[param_id]->cpu_data(), momntum, history[param_id]->mutable_cpu_data()); // 3. compute RMS[history] as denominator in temp dragon_set<Dtype>(count, eps, temp[param_id]->mutable_cpu_data()); dragon_axpy<Dtype>(count, Dtype(1), history[param_id]->cpu_data(),temp[param_id]->mutable_cpu_data()); dragon_powx<Dtype>(count, temp[param_id]->cpu_data(), Dtype(0.5), temp[param_id]->mutable_cpu_data()); // 4. compute diff/RMS[history] in diff dragon_div<Dtype>(count, net_param->cpu_diff(), temp[param_id]->cpu_data(), net_param->mutable_cpu_diff()); // 5. compute RMS[update] as numerator in temp dragon_set<Dtype>(count, eps, temp[param_id]->mutable_cpu_data()); dragon_axpy<Dtype>(count, Dtype(1), update[param_id]->cpu_data(), temp[param_id]->mutable_cpu_data()); dragon_powx<Dtype>(count, temp[param_id]->cpu_data(), Dtype(0.5), temp[param_id]->mutable_cpu_data()); // 6. compute diff*RMS[update] in diff dragon_mul<Dtype>(count, net_param->cpu_diff(), temp[param_id]->cpu_data(), net_param->mutable_cpu_diff()); // 7. compute final diff^2 in temp dragon_powx<Dtype>(count, net_param->cpu_diff(), Dtype(2), temp[param_id]->mutable_cpu_data()); // 8. compute update dragon_cpu_axpby<Dtype>(count, (1 - momntum), temp[param_id]->cpu_data(), momntum, update[param_id]->mutable_cpu_data()); // 9. apply learning rate dragon_scal<Dtype>(count, lr, net_param->mutable_cpu_diff()); break; case Dragon::GPU: #ifndef CPU_ONLY dragon_gpu_powx<Dtype>(count, net_param->gpu_diff(), Dtype(2), temp[param_id]->mutable_gpu_data()); // 2. compute history dragon_gpu_axpby<Dtype>(count, Dtype(1) - momntum, temp[param_id]->gpu_data(), momntum, history[param_id]->mutable_gpu_data()); // 3. compute RMS[history] as denominator in temp dragon_gpu_set<Dtype>(count, eps, temp[param_id]->mutable_gpu_data()); dragon_gpu_axpy<Dtype>(count, Dtype(1), history[param_id]->gpu_data(), temp[param_id]->mutable_gpu_data()); dragon_gpu_powx<Dtype>(count, temp[param_id]->gpu_data(), Dtype(0.5), temp[param_id]->mutable_gpu_data()); // 4. compute diff/RMS[history] in diff dragon_gpu_div<Dtype>(count, net_param->gpu_diff(), temp[param_id]->gpu_data(), net_param->mutable_gpu_diff()); // 5. compute RMS[update] as numerator in temp dragon_gpu_set<Dtype>(count, eps, temp[param_id]->mutable_gpu_data()); dragon_gpu_axpy<Dtype>(count, Dtype(1), update[param_id]->gpu_data(), temp[param_id]->mutable_gpu_data()); dragon_gpu_powx<Dtype>(count, temp[param_id]->gpu_data(), Dtype(0.5), temp[param_id]->mutable_gpu_data()); // 6. compute diff*RMS[update] in diff dragon_gpu_mul<Dtype>(count, net_param->gpu_diff(), temp[param_id]->gpu_data(), net_param->mutable_gpu_diff()); // 7. compute final diff^2 in temp dragon_gpu_powx<Dtype>(count, net_param->gpu_diff(), Dtype(2), temp[param_id]->mutable_gpu_data()); // 8. compute update dragon_gpu_axpby<Dtype>(count, Dtype(1) - momntum, temp[param_id]->gpu_data(), momntum, update[param_id]->mutable_gpu_data()); // 9. apply learning rate dragon_gpu_scal<Dtype>(count, lr, net_param->mutable_gpu_diff()); #endif break; default:LOG(FATAL) << "Unknown mode: " << Dragon::get_mode(); } } template <typename Dtype> void AdaDeltaSolver<Dtype>::applyUpdate(){ CHECK(Dragon::get_root_solver()); Dtype rate = getLearningRate(); // AdaDelta do not need base lr if (param.display() && iter%param.display() == 0) LOG(INFO) << "Iteration " << iter << ", lr = AdaDelta"; clipGradients(); vector<Blob<Dtype>*> net_params = net->getLearnableParams(); for (int i = 0; i < net_params.size(); i++){ normalize(i); regularize(i); computeUpdateValue(i, rate); net_params[i]->update(); } } INSTANTIATE_CLASS(AdaDeltaSolver);
AdaDelta的缺陷
局部最小值
从多个数据集情况来看，AdaDelta在训练初期和中期，具有非常不错的加速效果。
但是到训练后期，进入局部最小值雷区之后，AdaDelta就会反复在局部最小值附近抖动。
主要体现在验证集错误率上，脱离不了局部最小值吸引盆。
这时候，切换成动量SGD，如果把学习率降低一个量级，就会发现验证集正确率有2%~5%的提升，
这与常规使用动量SGD，是一样的。
之后再切换成AdaDelta，发现正确率又退回去了。
再切换成动量SGD，发现正确率又回来了。
---------------------------------------------------------------------
注：使用Batch Norm之后，这样从AdaDelta切到SGD会导致数值体系崩溃，原因未知。
---------------------------------------------------------------------
个人猜测，人工学习率的量级降低，给训练造成一个巨大的抖动，从一个局部最小值，
抖动到了另一个局部最小值，而AdaDelta的二阶近似计算，或者说所有二阶方法，
则不会产生这么大的抖动，所以很难从局部最小值中抖出来。
这给追求state of art的结果带来灾难，因为只要你一直用AdaDelta，肯定是与state of art无缘的。
基本上state of art的结果，最后都是SGD垂死挣扎抖出来的。
这也是SGD为什么至今在state of art的论文中没有废除的原因，人家丑，但是实在。
精度
eps的数值不是固定的。
1e-6在Caffe Cifar10上就显得过小了，1e-8比较适合。
这意味着不同数值比例体系，精度需要人工注意。
paper里高精度反而没低精度好，说明精度也有比较大抖动。
so，究竟什么样的精度是最好的呢？
————————————————————————————————————
2016.5.19 更新：
在FCNN-AlexNet里，1e-8在epoch1之后就会产生数值问题。
原因是sqrt(1e-8)*grad很大，这时候1e-10是比较好的。
另外，DensePrediction一定要做normalize，否则也有可能让AdaDelta的迭代步长计算出现数值问题。
该问题在FCNN-AlexNet进行到epoch5左右时候开始明显化。
caffe默认给的1e-10实际上要比paper里的1e-6要相对robust。