机器学习 --超参数调优

 机器学习中，绝大部分模型没有解析解，需要采用梯度下降法求解最有参数，各种各样的梯度下降法都会遇到一个问题，就是如何设置学习率，是一个技术活，更是一个运气活。

一 学习率参数调优的原理

超参数调优，经常会遇到两个问题：
1、 模型发散，参数随着迭代数值绝对值越来越大，甚至发散到无穷，从损失函数来看，误差也会越来越大。
2、 震荡，从损失函数来看，误差出现震荡，模型在局部最优解附近徘徊。

1.1 模型参数发散的原因

假定损失函数为$Q(\beta)\gt 0$，通常梯度下降法可表示为

$$\beta^{new}=\beta^{old}-\lambda \nabla Q(\beta)|_{\beta=\beta^{old}}$$
当 $Q(\beta)$满足条件：
$$\Vert \nabla Q(\beta)\Vert\ge s\Vert \beta\Vert;\; \forall \Vert \beta\Vert>M >0\;and \;s>0$$

几乎所有损失函数都满足这一条件，比如最简单的$Q(\theta)=a\theta^2,a>0,\forall \theta \ne0$,有$c=2a$。

有：

$$\begin{array}{left} \Vert\beta-\lambda \nabla Q(\beta)\Vert^2 &= \Vert\beta\Vert^2-2cos(\theta)\lambda \Vert\nabla Q(\beta)\Vert\Vert\beta\Vert+\lambda^2\Vert\nabla Q(\beta)\Vert\\ &=\Vert\beta\Vert^2+\lambda\Vert\nabla Q(\beta)\Vert(\Vert\nabla Q(\beta)\Vert-2cos(\theta)\Vert\beta\Vert)\\ &\ge \Vert\beta\Vert^2+\lambda\Vert\nabla Q(\beta)\Vert(\lambda c-2cos(\theta)\Vert\beta\Vert\\ \end{array}$$
当 $\lambda c>2$时必有：
$$\begin{array}{left} \Vert\beta-\lambda \nabla Q(\beta)\Vert^2 &\ge \Vert\beta\Vert^2+\lambda\Vert\nabla Q(\beta)\Vert(\lambda c-2cos(\theta)\Vert\beta\Vert\\ &\gt \lambda c\epsilon\Vert\beta\Vert^2+ \Vert\beta\Vert^2\\ &= (1+\lambda t\epsilon)\Vert\beta\Vert^2 \end{array}$$
即 $\Vert\beta^{k}\Vert\gt (1+\lambda c\epsilon)^{k/2}\Vert\beta^{0}\Vert$,此时迭代将使得参数指数增长，参数必发散且必有损失函数趋于无穷大。

 改变损失函数本身性质，必须改变损失函数，即便这样做也很难保证不出现参数爆炸的情况。剩下的就是改变超参数$\lambda$.
 任何时候，当你设定的学习率导致参数发散的时候，简单的办法就是降低学习率$\lambda$，破坏条件$\lambda c>2$从而避免参数发散。

1.2 震荡的原因

 不同于参数发散，损失函数值震荡，主要是因为在局部最优解附近学习率过大，导致迭代无法收敛到局部最优值。
假定$\beta^*$是局部最优值，在局部最优值附近有$\beta=t\alpha+\beta^*$,简单起见，不防设$\phi\alpha=\nabla^2 Q(\beta^*)\alpha$,即恰好$\alpha$是损失函数在局部最优解附近的海塞矩阵的特征方向。
则得到过程可以表示为

$$\begin{array}{left} \lambda \nabla Q(\beta) &\approx t\lambda \nabla^2 Q(\beta^*)\alpha =t\lambda \phi\alpha \\ &\Downarrow\\ \beta^{(0)}&\approx t\alpha+\beta^*\\ \beta^{(1)}&\approx \beta^{(0)}-t\lambda \phi\alpha=\beta^*+t(1-\lambda \phi)\alpha\\ \beta^{(2)}&\approx \beta^{(1)}-t(1-\lambda \phi)\lambda \phi\alpha=\beta^*+t(1-\lambda \phi)^2\alpha\\ \beta^{(k)}&\approx \beta^*+t(1-\lambda \phi)^k\alpha \end{array}$$
可以看到：
1. 若 $0\lt(1-\lambda \phi)<1$则迭代逐步趋向于 $\beta^*$;
2. 若 $0\gt(1-\lambda \phi)\gt -1$则迭代也会逐步趋向于 $\beta^*$，但每一步会跳过局部最优值;
3. 若 $(1-\lambda \phi)\gt 1$则迭代会逐步离开 $\beta^*$，有可能离开该局部最优值，也可能既无法离开也无法收敛，与损失函数更大范围的性质有关系，损失曲线一般表现出不规则震动;
4. 若 $(1-\lambda \phi)= -1$则迭代表现出在两个值之间周期性跳转，损失曲线上可能看不到变化，也可能表现为震动。
因此若要保证收敛，还是需要缩小学习率 $\lambda$.

1.3 其他情况

 复杂模型，还会遇到梯度爆炸和梯度消失的问题，梯度消失则表明模型已经训练到达一个局部最优解。而梯度爆炸的问题是训练过程中更为麻烦的问题，它会迫使搜索脱离当前区域。
 样本数量改变经常导致损失函数变化，因此在小样本上训练收敛的参数，在大样本下可能仍然不收敛，这导致工程应用中有麻烦，没法设定固定的学习率。

二 实战研究

 下面以简单回归模型为例，说明学习率参数对机器学习收敛的影响。一般而言学习率由大到小，先导致发散，缩小学习率后误差震荡，继续缩小则训练过程收敛。 跟混沌理论差不多，导致震荡学习率有一个临界区间。同时学习率并不是越小越好，缩小学习率需要更多的迭代次数才能。
 给定模型$y=5x_1+6x_2+abx_3$,产生随机样本：

Num=1000
set.seed(1000)
x1 = runif(n = Num,min = -10,max = 3)
x2 = runif(n = Num,min = -5,max = 5)
x3 = runif(n = Num,min = -5,max = 10)
X=cbind(x1,x2,x3)
y=5*x1+6*x2+5*6*x3
sd(y)
# 添加随机误差
y=y+rnorm(Num,mean = 0,sd = 1)

模型的损失函数为：$Q(a,b)=\frac{1}{2}\sum_{i=1}^n(y_i-ax_{i1}-bx_{i2}-abx_{i3})^2$，

梯度为:

$$\begin{array}{left} \frac{\partial Q(a)}{ \partial a}&=-\sum_{i=1}^n(y_i-ax_{i1}-bx_{i2}-abx_{i3})(x_{i1}+bx_{i3})\\ &=-\sum_{i=1}^n(y_ix_{i1}-ax_{i1}^2-bx_{i1}x_{i2}-abx_{i1}x_{i3}+b(y_ix_{i3}-ax_{i1}x_{i3}-bx_{i2}x_{i3}-abx_{i3}^2))\\ \frac{\partial Q(a)}{ \partial b}&=-\sum_{i=1}^n(y_i-ax_{i1}-bx_2-abx_{i3})(x_{i2}+ax_{i3})\\ &=-\sum_{i=1}^n(y_ix_{i2}-ax_{i1}x_{i1}x_{i2}-bx_{i2}^2-abx_{i2}x_{i3}+a(y_ix_{i3}-ax_{i1}x_{i3}-bx_{i2}x_{i3}-abx_{i3}^2))\\ \end{array}$$
梯度下降迭代公式为：
$$a^{(k+1)}=a^{(k)}+\lambda\sum_{i=1}^n(y_i-ax_{i1}-bx_{i2}-abx_{i3})(x_{i1}+bx_{i3})\\ b^{(k+1)}=b^{(k)}+\lambda\sum_{i=1}^n(y_i-ax_{i1}-bx_2-abx_{i3})(x_{i2}+ax_{i3})$$
R语言实现迭代过程：

olm<-function(y,x,lambda=0.02,iter=100,beta=runif(n = 2,min = 0.01,max = 0.2)){
  stopifnot(length(y)==dim(x)[1])
   stopifnot(3==dim(x)[2])
  xy=c(t(y)%*%x)
  xx=t(x)%*%x
  #为了观察内部过程，需要将每一步迭代后的参数保留下来
  betas=data.frame(a=rep(0,iter+1),b=rep(0,iter+1))
  betas[1,]=beta
  for(i in 1:iter){
    b0=c(beta,beta[1]*beta[2])
    b1=matrix(data=c(1,0,0,1,beta[2:1]),nrow = 2)
    #browser()
    beta=beta +lambda*b1%*%(xy-xx %*%b0)
    betas[i+1,]=beta
  }
  betas
}
#betas=olm(y,X,lambda = 0.00056,iter = 1000)

计算每一步迭代后个损失函数值：

betas=olm(y,X,lambda =0.0000028,iter = 500)
beta2=betas
beta2[,3]=beta2[,1]*beta2[,2]
err=0.5*colSums((y-X%*%t(beta2))^2)/Num
plot(err,type = "l")

plot(betas,type="b")

 从上图可以看出，此参数使得训练在若干步后发散，损失值突然发散到无穷大，训练模型参数由于最后太大，导致之前的模型参数被压缩在一个极小的区域，从而在图上显示为一个点。

 重新调整初始化参数，缩小学习率参数,同时加大迭代次数,可以从下图中看到，这一次明显有误差震荡。

betas1=olm(y,X,lambda = 0.000002,iter = 1000)
beta2=betas1
beta2[,3]=beta2[,1]*beta2[,2]
err1=0.5*colSums((y-X%*%t(beta2))^2)/Num
plot(err1,type = "l")

plot(betas1,type="b")

 上面的参数，使得参数在最优值附近跳转，在损失函数值上看来，表现为震荡。迭代参数基本在一条直线上，这并不是偶然的，这条直线会逐步逼近损失函数在改局部最小值的Hassan矩阵的特征方向，一般而言是有最大特征值的特征方向。来看最后50次迭代情况。

plot(betas1[950:1000,],type="b")

继续缩小学习率。

betas1=olm(y,X,lambda =0.0000014,iter = 1000)
beta2=betas1
beta2[,3]=beta2[,1]*beta2[,2]
err1=0.5*colSums((y-X%*%t(beta2))^2)/Num
plot(err1[-c(1:2)],type = "l")

plot(betas1,type="b")

plot(betas1[950:1000,],type="o")

 这时，训练误差迅速逼近最小值，从模型参数迭代过程来看，参数迅速靠近（5,6）附近从最后一图来看，模型参数值移动范围越来越小。

三、梯度下降法的改进思路

 不管是参数发散还是震荡，归根揭底都是由于给定学习率后，跟梯度模长有关系，实际上在梯度越大的地方，损失函数变化越大，这时我们恰好需要更小的学习率以免模型参数下一步跑过最小值。事实上梯度下降法的关键是梯度方向，而不是梯度的模长，而模长依赖于样本数据，尤其当样本数量变化时，对梯度影响更加明显，更多的样本量意味着需要更小的参数，而机器学习需要的样本越多越好。
 随机梯度下降法可以控制每次参与，但每次参与计算样本数量仍然对学习率参数有显著的影响。比较更好的办法是控制梯度模长，因为对于模型训练结果而言，尽需要保证每次迭代的梯度方向大致正确。
比随机梯度更彻底的办法是，是控制梯度模长。保障$V= t\nabla Q(\beta)|_{\beta=\beta^{old}} \;;\quad\beta^{new}=\beta^{old}-\lambda V$迭代参数仍然可以做到收敛。若能任何时候保障$|V|<M$，有

$$\begin{array}{left} \Vert\beta-\lambda V\Vert^2 &= \Vert\beta\Vert^2-2cos(\theta)\lambda \Vert V\Vert\Vert\beta\Vert+\lambda^2\Vert V\Vert^2\\ &=\Vert\beta\Vert^2+\lambda\Vert V\Vert(\lambda\Vert V\Vert-2cos(\theta)\Vert\beta\Vert)\\ \end{array}$$
 只要给定 $\lambda\Vert V\Vert<\lambda M<2$就能破坏参数发散的条件，也能避免当局部梯度太大时，模型参数 $\beta$脱离当前搜索区域，使得梯度下降法搜索趋于稳定。
这个改动可以避免参数发散的问题，也能避免梯度爆炸的问题，代价是训练过程需要更多的迭代次数，增加计算复杂。
一般我们可以取:
$$\begin{array}{left} V=\frac{\nabla Q(\beta)}{max(\Vert \nabla Q(\beta)\Vert ),1)} \end{array}$$
 此时，从而有 $||V||\le1$取范数取 $p^{-\infty}$,计算过程相对简单。
 令一种办法是取 $V=sign(\nabla Q(\beta))$，能大致保持梯度方向，并且有 $||V||\le m$， m m <script type="math/tex" id="MathJax-Element-4404">m</script>是参数个数，这个计算也很简单，尤其是在分布式集群上运算的时候。