机器学习中导数最优化方法(基础篇)

作者：daniel-D

1. 前言

熟悉机器学习的童鞋都知道，优化方法是其中一个非常重要的话题，最常见的情形就是利用目标函数的导数通过多次迭代来求解无约束最优化问题。实现简单，coding 方便，是训练模型的必备利器之一。这篇博客主要总结一下使用导数的最优化方法的几个基本方法，梳理梳理相关的数学知识，本人也是一边写一边学，如有问题，欢迎指正，共同学习，一起进步。

2. 几个数学概念

1) 梯度（一阶导数）

考虑一座在 (x1, x2) 点高度是 f(x1, x2) 的山。那么，某一点的梯度方向是在该点坡度最陡的方向，而梯度的大小告诉我们坡度到底有多陡。注意，梯度也可以告诉我们不在最快变化方向的其他方向的变化速度（二维情况下，按照梯度方向倾斜的圆在平面上投影成一个椭圆）。对于一个含有 n 个变量的标量函数，即函数输入一个 n 维 的向量，输出一个数值，梯度可以定义为：

2) Hesse 矩阵（二阶导数）

Hesse 矩阵常被应用于 牛顿法解决的大规模优化问题 (后面会介绍)，主要形式如下：

3) Jacobi 矩阵

Jacobi 矩阵实际上是向量值函数的梯度矩阵，假设F:Rn→Rm 是一个从n维欧氏空间转换到m维欧氏空间的函数。这个函数由m个实函数组成: 。这些函数的偏导数(如果存在)可以组成一个m行n列的矩阵(m by n)，这就是所谓的雅可比矩阵：

a)如果 f(x) 是一个标量函数，那么雅克比矩阵是一个向量，等于 f(x) 的梯度， Hesse 矩阵是一个二维矩阵。如果 f(x) 是一个向量值函数，那么Jacobi 矩阵是一个二维矩阵，Hesse 矩阵是一个三维矩阵。

b)梯度是 Jacobian 矩阵的特例，梯度的 jacobian 矩阵就是 Hesse 矩阵（一阶偏导与二阶偏导的关系）。

3. 优化方法

1) Gradient Descent

Gradient descent 又叫 steepest descent，是利用一阶的梯度信息找到函数局部最优解的一种方法，也是机器学习里面最简单最常用的一种优化方法。Gradient descent 是 line search 方法中的一种，主要迭代公式如下：

# Code from Chapter 11 of Machine Learning: An Algorithmic Perspective # by Stephen Marsland (http://seat.massey.ac.nz/personal/s.r.marsland/MLBook.html) # Gradient Descent using steepest descent from numpy import * def Jacobian(x):  return array([x[0], 0.4*x[1], 1.2*x[2]]) def steepest(x0):  i = 0   iMax = 10  x = x0  Delta = 1  alpha = 1  while i<iMax and Delta>10**(-5):   p = -Jacobian(x)   xOld = x   x = x + alpha*p   Delta = sum((x-xOld)**2)   print 'epoch', i, ':'   print x, '/n'   i += 1 x0 = array([-2,2,-2]) steepest(x0) 

Steepest gradient 方法得到的是局部最优解，如果目标函数是一个凸优化问题，那么局部最优解就是全局最优解，理想的优化效果如下图，值得注意一点的是，每一次迭代的移动方向都与出发点的等高线垂直：

粗略来讲，在二次函数中，椭球面的形状受 hesse 矩阵的条件数影响，长轴与短轴对应矩阵的最小特征值和最大特征值的方向，其大小与特征值的平方根成反比，最大特征值与最小特征值相差越大，椭球面越扁，那么优化路径需要走很大的弯路，计算效率很低。

2) Newton’s method

在最速下降法中，我们看到，该方法主要利用的是目标函数的局部性质，具有一定的“盲目性”。牛顿法则是利用局部的一阶和二阶偏导信息，推测整个目标函数的形状，进而可以求得出近似函数的全局最小值，然后将当前的最小值设定近似函数的最小值。相比最速下降法，牛顿法带有一定对全局的预测性，收敛性质也更优良。牛顿法的主要推导过程如下：

第一步，利用 Taylor 级数求得原目标函数的二阶近似：

与 1) 中优化问题相同，牛顿法的代码如下：

Newton.py

# Code from Chapter 11 of Machine Learning: An Algorithmic Perspective # by Stephen Marsland (http://seat.massey.ac.nz/personal/s.r.marsland/MLBook.html) # Gradient Descent using Newton's method from numpy import * def Jacobian(x):  return array([x[0], 0.4*x[1], 1.2*x[2]]) def Hessian(x):  return array([[1,0,0],[0,0.4,0],[0,0,1.2]]) def Newton(x0):  i = 0  iMax = 10  x = x0  Delta = 1  alpha = 1  while i<iMax and Delta>10**(-5):   p = -dot(linalg.inv(Hessian(x)),Jacobian(x))   xOld = x   x = x + alpha*p   Delta = sum((x-xOld)**2)   i += 1  print x x0 = array([-2,2,-2]) Newton(x0) 

上面例子中由于目标函数是二次凸函数，Taylor 展开等于原函数，所以能一次就求出最优解。

牛顿法主要存在的问题是：

1. Hesse 矩阵不可逆时无法计算
2. 矩阵的逆计算复杂为 n 的立方，当问题规模比较大时，计算量很大，解决的办法是采用拟牛顿法如 BFGS, L-BFGS, DFP, Broyden’s Algorithm 进行近似。
3. 如果初始值离局部极小值太远，Taylor 展开并不能对原函数进行良好的近似

3) Levenberg–Marquardt Algorithm

Levenberg–Marquardt algorithm 能结合以上两种优化方法的优点，并对两者的不足做出改进。与 line search 的方法不同，LMA 属于一种“信赖域法”(trust region)，牛顿法实际上也可以看做一种信赖域法，即利用局部信息对函数进行建模近似，求取局部最小值。所谓的信赖域法，就是从初始点开始，先假设一个可以信赖的最大位移 s（牛顿法里面 s 为无穷大），然后在以当前点为中心，以 s 为半径的区域内，通过寻找目标函数的一个近似函数（二次的）的最优点，来求解得到真正的位移。在得到了位移之后，再计算目标函数值，如果其使目标函数值的下降满足了一定条件，那么就说明这个位移是可靠的，则继续按此规则迭代计算下去；如果其不能使目标函数值的下降满足一定的条件，则应减小信赖域的范围，再重新求解。

LMA 最早提出是用来解决最小二乘法曲线拟合的优化问题的，对于随机初始化的已知参数 beta， 求得的目标值为：

算法过程如下：

1. 给定一个初识值 x0
2. 当 并且没有到达最大迭代次数时
3. 重复执行:
   • 算出移动向量
   • 计算更新值：
   • 计算目标函数真实减少量与预测减少量的比率
   • if ，接受更新值
   • else if ，说明近似效果很好，接受更新值，扩大可信域（即减小阻尼系数）
   • else: 目标函数在变大，拒绝更新值，减小可信域（即增加阻尼系数）
4. 直到达到最大迭代次数

维基百科在介绍 Gradient descent 时用包含了细长峡谷的 Rosenbrock function

LevenbergMarquardt.py

# Code from Chapter 11 of Machine Learning: An Algorithmic Perspective # by Stephen Marsland (http://seat.massey.ac.nz/personal/s.r.marsland/MLBook.html) # The Levenberg Marquardt algorithm from numpy import * def function(p):  r = array([10*(p[1]-p[0]**2),(1-p[0])])  fp = dot(transpose(r),r) #= 100*(p[1]-p[0]**2)**2 + (1-p[0])**2  J = (array([[-20*p[0],10],[-1,0]]))  grad = dot(transpose(J),transpose(r))  return fp,r,grad,J def lm(p0,tol=10**(-5),maxits=100):  nvars=shape(p0)[0]  nu=0.01  p = p0  fp,r,grad,J = function(p)  e = sum(dot(transpose(r),r))  nits = 0  while nits<maxits and linalg.norm(grad)>tol:   nits += 1   fp,r,grad,J = function(p)   H=dot(transpose(J),J) + nu*eye(nvars)   pnew = zeros(shape(p))   nits2 = 0   while (p!=pnew).all() and nits2<maxits:    nits2 += 1    dp,resid,rank,s = linalg.lstsq(H,grad)    pnew = p - dp    fpnew,rnew,gradnew,Jnew = function(pnew)    enew = sum(dot(transpose(rnew),rnew))    rho = linalg.norm(dot(transpose(r),r)-dot(transpose(rnew),rnew))    rho /= linalg.norm(dot(transpose(grad),pnew-p))    if rho>0:     update = 1     p = pnew     e = enew     if rho>0.25:      nu=nu/10    else:      nu=nu*10     update = 0   print fp, p, e, linalg.norm(grad), nu p0 = array([-1.92,2]) lm(p0) 

大概 5 次迭代就可以得到最优解 (1, 1).

Levenberg–Marquardt algorithm 对局部极小值很敏感，维基百科举了一个二乘法曲线拟合的例子，当使用不同的初始值时，得到的结果差距很大，我这里也有 python 代码，就不细说了。

4) Conjugate Gradients

共轭梯度法也是优化模型经常经常要用到的一个方法，背后的数学公式和原理稍微复杂一些，光这一个优化方法就可以写一篇很长的博文了，所以这里并不打算详细讲解每一步的推导过程，只简单写一下算法的实现过程。与最速梯度下降的不同，共轭梯度的优点主要体现在选择搜索方向上。在了解共轭梯度法之前，我们首先简单了解一下共轭方向：

在优化过程中，如果我们确定了移动方向（GD：垂直于等值线，CG：共轭方向），然后在该方向上搜索极小值点（恰好与该处的等值线相切），然后移动到最小值点，重复以上过程，那么 Gradient Descent 和 Conjugate gradient descent 的优化过程可以用下图的 绿线 与 红线 表示：

1. 给定一个出发点 x0 和一个停止参数 e, 第一次移动方向为 最速下降方向 :
2. 
3. while :
   • 用 Newton-Raphson 迭代计算移动距离，以便在该搜索方向移动到极小，公式就不写了，具体思路就是利用一阶梯度的信息向极小值点跳跃搜索
   • 移动当前的优化解 x：
   • 用 Gram-Schmidt 方法构造下一个共轭方向，即 , 按照 的确定公式又可以分为 FR 方法和 PR 和 HS 等。

在很多的资料中，介绍共轭梯度法都举了一个求线性方程组 Ax = b 近似解的例子，实际上就相当于这里所说的

还是用最开始的目标函数 来编写共轭梯度法的优化代码：

# Code from Chapter 11 of Machine Learning: An Algorithmic Perspective # by Stephen Marsland (http://seat.massey.ac.nz/personal/s.r.marsland/MLBook.html) # The conjugate gradients algorithm from numpy import * def Jacobian(x):  #return array([.4*x[0],2*x[1]])  return array([x[0], 0.4*x[1], 1.2*x[2]]) def Hessian(x):  #return array([[.2,0],[0,1]])  return array([[1,0,0],[0,0.4,0],[0,0,1.2]]) def CG(x0):  i=0  k=0  r = -Jacobian(x0)  p=r  betaTop = dot(r.transpose(),r)  beta0 = betaTop  iMax = 3  epsilon = 10**(-2)  jMax = 5  # Restart every nDim iterations  nRestart = shape(x0)[0]  x = x0  while i < iMax and betaTop > epsilon**2*beta0:   j=0   dp = dot(p.transpose(),p)   alpha = (epsilon+1)**2   # Newton-Raphson iteration   while j < jMax and alpha**2 * dp > epsilon**2:    # Line search    alpha = -dot(Jacobian(x).transpose(),p) / (dot(p.transpose(),dot(Hessian(x),p)))    print "N-R",x, alpha, p    x = x + alpha * p    j += 1   print x   # Now construct beta   r = -Jacobian(x)   print "r: ", r   betaBottom = betaTop   betaTop = dot(r.transpose(),r)   beta = betaTop/betaBottom   print "Beta: ",beta   # Update the estimate   p = r + beta*p   print "p: ",p   print "----"   k += 1   if k==nRestart or dot(r.transpose(),p) <= 0:    p = r    k = 0    print "Restarting"   i +=1  print x x0 = array([-2,2,-2]) CG(x0) 

参考资料：