机器学习中导数最优化方法(基础篇)
1. 前言
熟悉機(jī)器學(xué)習(xí)的童鞋都知道，優(yōu)化方法是其中一個非常重要的話題，最常見的情形就是利用目標(biāo)函數(shù)的導(dǎo)數(shù)通過多次迭代來求解無約束最優(yōu)化問題。實(shí)現(xiàn)簡單，coding 方便，是訓(xùn)練模型的必備利器之一。這篇博客主要總結(jié)一下使用導(dǎo)數(shù)的最優(yōu)化方法的幾個基本方法，梳理梳理相關(guān)的數(shù)學(xué)知識，本人也是一邊寫一邊學(xué)，如有問題，歡迎指正，共同學(xué)習(xí)，一起進(jìn)步。
?
2. 幾個數(shù)學(xué)概念
1) 梯度（一階導(dǎo)數(shù)）
考慮一座在 (x1, x2) 點(diǎn)高度是 f(x1, x2) 的山。那么，某一點(diǎn)的梯度方向是在該點(diǎn)坡度最陡的方向，而梯度的大小告訴我們坡度到底有多陡。注意，梯度也可以告訴我們不在最快變化方向的其他方向的變化速度（二維情況下，按照梯度方向傾斜的圓在平面上投影成一個橢圓）。對于一個含有 n 個變量的標(biāo)量函數(shù)，即函數(shù)輸入一個 n 維 的向量，輸出一個數(shù)值，梯度可以定義為：
2) Hesse 矩陣（二階導(dǎo)數(shù)）
Hesse 矩陣常被應(yīng)用于牛頓法解決的大規(guī)模優(yōu)化問題(后面會介紹)，主要形式如下：
當(dāng) f(x) 為二次函數(shù)時，梯度以及 Hesse 矩陣很容易求得。二次函數(shù)可以寫成下列形式：
其中 A 是 n 階對稱矩陣，b 是 n 維列向量， c 是常數(shù)。f(x) 梯度是 Ax+b, Hesse 矩陣等于 A。
3) Jacobi 矩陣
Jacobi 矩陣實(shí)際上是向量值函數(shù)的梯度矩陣，假設(shè)F:Rn→Rm 是一個從n維歐氏空間轉(zhuǎn)換到m維歐氏空間的函數(shù)。這個函數(shù)由m個實(shí)函數(shù)組成:?。這些函數(shù)的偏導(dǎo)數(shù)(如果存在)可以組成一個m行n列的矩陣(m by n)，這就是所謂的雅可比矩陣：
總結(jié)一下,
a)?如果 f(x) 是一個標(biāo)量函數(shù)，那么雅克比矩陣是一個向量，等于 f(x) 的梯度， Hesse 矩陣是一個二維矩陣。如果 f(x) 是一個向量值函數(shù)，那么Jacobi 矩陣是一個二維矩陣，Hesse 矩陣是一個三維矩陣。
b)?梯度是 Jacobian 矩陣的特例，梯度的 jacobian 矩陣就是 Hesse 矩陣（一階偏導(dǎo)與二階偏導(dǎo)的關(guān)系）。
?
3. 優(yōu)化方法
1) Gradient Descent
Gradient descent 又叫 steepest descent，是利用一階的梯度信息找到函數(shù)局部最優(yōu)解的一種方法，也是機(jī)器學(xué)習(xí)里面最簡單最常用的一種優(yōu)化方法。Gradient descent 是 line search 方法中的一種，主要迭代公式如下：
其中，是第 k 次迭代我們選擇移動的方向，在 steepest descent 中，移動的方向設(shè)定為梯度的負(fù)方向，?是第 k 次迭代用 line search 方法選擇移動的距離，每次移動的距離系數(shù)可以相同，也可以不同，有時候我們也叫學(xué)習(xí)率（learning rate）。在數(shù)學(xué)上，移動的距離可以通過 line search 令導(dǎo)數(shù)為零找到該方向上的最小值，但是在實(shí)際編程的過程中，這樣計算的代價太大，我們一般可以將它設(shè)定位一個常量?？紤]一個包含三個變量的函數(shù)?，計算梯度得到?。設(shè)定 learning rate = 1，算法代碼如下：
# 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 descentfrom numpy import *def Jacobian(x):return array([x[0], 0.4*x[1], 1.2*x[2]])def steepest(x0):i = 0 iMax = 10x = x0Delta = 1alpha = 1while i<iMax and Delta>10**(-5):p = -Jacobian(x)xOld = xx = x + alpha*pDelta = sum((x-xOld)**2)print 'epoch', i, ':'print x, '\n'i += 1x0 = array([-2,2,-2]) steepest(x0)Steepest gradient 方法得到的是局部最優(yōu)解，如果目標(biāo)函數(shù)是一個凸優(yōu)化問題，那么局部最優(yōu)解就是全局最優(yōu)解，理想的優(yōu)化效果如下圖，值得注意一點(diǎn)的是，每一次迭代的移動方向都與出發(fā)點(diǎn)的等高線垂直：
需要指出的是，在某些情況下，最速下降法存在鋸齒現(xiàn)象（ zig-zagging）將會導(dǎo)致收斂速度變慢:
粗略來講，在二次函數(shù)中，橢球面的形狀受 hesse 矩陣的條件數(shù)影響，長軸與短軸對應(yīng)矩陣的最小特征值和最大特征值的方向，其大小與特征值的平方根成反比，最大特征值與最小特征值相差越大，橢球面越扁，那么優(yōu)化路徑需要走很大的彎路，計算效率很低。
2) Newton's method
在最速下降法中，我們看到，該方法主要利用的是目標(biāo)函數(shù)的局部性質(zhì)，具有一定的“盲目性”。牛頓法則是利用局部的一階和二階偏導(dǎo)信息，推測整個目標(biāo)函數(shù)的形狀，進(jìn)而可以求得出近似函數(shù)的全局最小值，然后將當(dāng)前的最小值設(shè)定近似函數(shù)的最小值。相比最速下降法，牛頓法帶有一定對全局的預(yù)測性，收斂性質(zhì)也更優(yōu)良。牛頓法的主要推導(dǎo)過程如下：
第一步，利用 Taylor 級數(shù)求得原目標(biāo)函數(shù)的二階近似：
第二步，把?x 看做自變量， 所有帶有 x^k 的項(xiàng)看做常量，令一階導(dǎo)數(shù)為 0 ，即可求近似函數(shù)的最小值：
即：
第三步，將當(dāng)前的最小值設(shè)定近似函數(shù)的最小值（或者乘以步長)。
與?1)?中優(yōu)化問題相同，牛頓法的代碼如下：
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 = 0iMax = 10x = x0Delta = 1alpha = 1while i<iMax and Delta>10**(-5):p = -dot(linalg.inv(Hessian(x)),Jacobian(x))xOld = xx = x + alpha*pDelta = sum((x-xOld)**2)i += 1print xx0 = array([-2,2,-2]) Newton(x0)上面例子中由于目標(biāo)函數(shù)是二次凸函數(shù)，Taylor 展開等于原函數(shù)，所以能一次就求出最優(yōu)解。
牛頓法主要存在的問題是：
3) Levenberg–Marquardt Algorithm
Levenberg–Marquardt algorithm 能結(jié)合以上兩種優(yōu)化方法的優(yōu)點(diǎn)，并對兩者的不足做出改進(jìn)。與 line search 的方法不同，LMA 屬于一種“信賴域法”(trust region)，牛頓法實(shí)際上也可以看做一種信賴域法，即利用局部信息對函數(shù)進(jìn)行建模近似，求取局部最小值。所謂的信賴域法，就是從初始點(diǎn)開始，先假設(shè)一個可以信賴的最大位移 s（牛頓法里面 s 為無窮大），然后在以當(dāng)前點(diǎn)為中心，以 s 為半徑的區(qū)域內(nèi)，通過尋找目標(biāo)函數(shù)的一個近似函數(shù)（二次的）的最優(yōu)點(diǎn)，來求解得到真正的位移。在得到了位移之后，再計算目標(biāo)函數(shù)值，如果其使目標(biāo)函數(shù)值的下降滿足了一定條件，那么就說明這個位移是可靠的，則繼續(xù)按此規(guī)則迭代計算下去；如果其不能使目標(biāo)函數(shù)值的下降滿足一定的條件，則應(yīng)減小信賴域的范圍，再重新求解。
LMA 最早提出是用來解決最小二乘法曲線擬合的優(yōu)化問題的，對于隨機(jī)初始化的已知參數(shù) beta， 求得的目標(biāo)值為：
對擬合曲線函數(shù)進(jìn)行一階 Jacobi 矩陣的近似：
進(jìn)而推測出 S 函數(shù)的周邊信息：
位移是多少時得到 S 函數(shù)的最小值呢？通過幾何的概念，當(dāng)殘差??垂直于 J 矩陣的 span 空間時， S 取得最小（至于為什么？請參考之前博客的最后一部分)
我們將這個公式略加修改，加入阻尼系數(shù)得到：
就是萊文貝格－馬夸特方法。這種方法只計算了一階偏導(dǎo)，而且不是目標(biāo)函數(shù)的 Jacobia 矩陣，而是擬合函數(shù)的 Jacobia 矩陣。當(dāng)??大的時候可信域小，這種算法會接近最速下降法，??小的時候可信域大，會接近高斯-牛頓方法。
算法過程如下：
- 算出移動向量
- 計算更新值：
- 計算目標(biāo)函數(shù)真實(shí)減少量與預(yù)測減少量的比率?
- if?，接受更新值
- else if?，說明近似效果很好，接受更新值，擴(kuò)大可信域（即減小阻尼系數(shù)）
- else: 目標(biāo)函數(shù)在變大，拒絕更新值，減小可信域（即增加阻尼系數(shù)）
維基百科在介紹 Gradient descent 時用包含了細(xì)長峽谷的 Rosenbrock function
展示了 zig-zagging 鋸齒現(xiàn)象：
用 LMA 優(yōu)化效率如何。套用到我們之前 LMA 公式中，有：
代碼如下：
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])**2J = (array([[-20*p[0],10],[-1,0]]))grad = dot(transpose(J),transpose(r))return fp,r,grad,Jdef lm(p0,tol=10**(-5),maxits=100):nvars=shape(p0)[0]nu=0.01p = p0fp,r,grad,J = function(p)e = sum(dot(transpose(r),r))nits = 0while nits<maxits and linalg.norm(grad)>tol:nits += 1fp,r,grad,J = function(p)H=dot(transpose(J),J) + nu*eye(nvars)pnew = zeros(shape(p))nits2 = 0while (p!=pnew).all() and nits2<maxits:nits2 += 1dp,resid,rank,s = linalg.lstsq(H,grad)pnew = p - dpfpnew,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 = 1p = pnewe = enewif rho>0.25:nu=nu/10else: nu=nu*10update = 0print fp, p, e, linalg.norm(grad), nup0 = array([-1.92,2]) lm(p0)大概 5 次迭代就可以得到最優(yōu)解 (1, 1).
Levenberg–Marquardt algorithm 對局部極小值很敏感，維基百科舉了一個二乘法曲線擬合的例子，當(dāng)使用不同的初始值時，得到的結(jié)果差距很大，我這里也有 python 代碼，就不細(xì)說了。
4) Conjugate Gradients
共軛梯度法也是優(yōu)化模型經(jīng)常經(jīng)常要用到的一個方法，背后的數(shù)學(xué)公式和原理稍微復(fù)雜一些，光這一個優(yōu)化方法就可以寫一篇很長的博文了，所以這里并不打算詳細(xì)講解每一步的推導(dǎo)過程，只簡單寫一下算法的實(shí)現(xiàn)過程。與最速梯度下降的不同，共軛梯度的優(yōu)點(diǎn)主要體現(xiàn)在選擇搜索方向上。在了解共軛梯度法之前，我們首先簡單了解一下共軛方向：
共軛方向和馬氏距離的定義有類似之處，他們都考慮了全局的數(shù)據(jù)分布。如上圖，d(1) 方向與二次函數(shù)的等值線相切，d(1) 的共軛方向 d(2) 則指向橢圓的中心。所以對于二維的二次函數(shù)，如果在兩個共軛方向上進(jìn)行一維搜索，經(jīng)過兩次迭代必然達(dá)到最小點(diǎn)。前面我們說過，等值線橢圓的形狀由 Hesse 矩陣決定，那么，上圖的兩個方向關(guān)于 Hessen 矩陣正交，共軛方向的定義如下：
如果橢圓是一個正圓， Hessen 矩陣是一個單位矩陣，上面等價于歐幾里得空間中的正交。
在優(yōu)化過程中，如果我們確定了移動方向（GD：垂直于等值線，CG：共軛方向），然后在該方向上搜索極小值點(diǎn)（恰好與該處的等值線相切），然后移動到最小值點(diǎn)，重復(fù)以上過程，那么 Gradient Descent 和 Conjugate gradient descent 的優(yōu)化過程可以用下圖的綠線與紅線表示：
講了這么多，共軛梯度算法究竟是如何算的呢？
- 用 Newton-Raphson 迭代計算移動距離，以便在該搜索方向移動到極小，公式就不寫了，具體思路就是利用一階梯度的信息向極小值點(diǎn)跳躍搜索
- 移動當(dāng)前的優(yōu)化解 x：?
- 用 Gram-Schmidt 方法構(gòu)造下一個共軛方向，即??, 按照??的確定公式又可以分為 FR 方法和 PR 和 HS 等。
在很多的資料中，介紹共軛梯度法都舉了一個求線性方程組 Ax = b 近似解的例子，實(shí)際上就相當(dāng)于這里所說的?
還是用最開始的目標(biāo)函數(shù) ?? ?來編寫共軛梯度法的優(yōu)化代碼：
?
# 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 algorithmfrom 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=0k=0r = -Jacobian(x0)p=rbetaTop = dot(r.transpose(),r)beta0 = betaTopiMax = 3epsilon = 10**(-2)jMax = 5# Restart every nDim iterationsnRestart = shape(x0)[0]x = x0while i < iMax and betaTop > epsilon**2*beta0:j=0dp = dot(p.transpose(),p)alpha = (epsilon+1)**2# Newton-Raphson iterationwhile j < jMax and alpha**2 * dp > epsilon**2:# Line searchalpha = -dot(Jacobian(x).transpose(),p) / (dot(p.transpose(),dot(Hessian(x),p)))print "N-R",x, alpha, px = x + alpha * pj += 1print x# Now construct betar = -Jacobian(x)print "r: ", rbetaBottom = betaTopbetaTop = dot(r.transpose(),r)beta = betaTop/betaBottomprint "Beta: ",beta# Update the estimatep = r + beta*pprint "p: ",pprint "----"k += 1if k==nRestart or dot(r.transpose(),p) <= 0:p = rk = 0print "Restarting"i +=1print xx0 = array([-2,2,-2]) CG(x0)?
??
?
?
參考資料：
[1] Machine Learning: An Algorithmic Perspective, chapter 11
[2] 最優(yōu)化理論與算法（第2版），陳寶林
[3] wikipedia
總結(jié)
以上是生活随笔為你收集整理的机器学习中导数最优化方法(基础篇)的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 漫谈：机器学习中距离和相似性度量方法
- 下一篇： Scala 中的函数式编程基础