? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?幾種常用的優(yōu)化方法

1. 前言

熟悉機(jī)器學(xué)習(xí)的童鞋都知道，優(yōu)化方法是其中一個(gè)非常重要的話題，最常見的情形就是利用目標(biāo)函數(shù)的導(dǎo)數(shù)通過多次迭代來求解無約束最優(yōu)化問題。實(shí)現(xiàn)簡單，coding 方便，是訓(xùn)練模型的必備利器之一。

?

2. 幾個(gè)數(shù)學(xué)概念

1) 梯度（一階導(dǎo)數(shù)）

考慮一座在 (x1, x2) 點(diǎn)高度是 f(x1, x2) 的山。那么，某一點(diǎn)的梯度方向是在該點(diǎn)坡度最陡的方向，而梯度的大小告訴我們坡度到底有多陡。注意，梯度也可以告訴我們不在最快變化方向的其他方向的變化速度（二維情況下，按照梯度方向傾斜的圓在平面上投影成一個(gè)橢圓）。對于一個(gè)含有 n 個(gè)變量的標(biāo)量函數(shù)，即函數(shù)輸入一個(gè) n 維 的向量，輸出一個(gè)數(shù)值，梯度可以定義為：

2) Hesse 矩陣（二階導(dǎo)數(shù)）

Hesse 矩陣常被應(yīng)用于牛頓法解決的大規(guī)模優(yōu)化問題(后面會(huì)介紹)，主要形式如下：

當(dāng) f(x) 為二次函數(shù)時(shí)，梯度以及 Hesse 矩陣很容易求得。二次函數(shù)可以寫成下列形式：

其中 x為列向量，A 是 n 階對稱矩陣，b 是 n 維列向量， c 是常數(shù)。f(x) 梯度是 Ax+b, Hesse 矩陣等于 A。

3) Jacobi 矩陣

Jacobi 矩陣實(shí)際上是向量值函數(shù)的梯度矩陣，假設(shè)F:Rn→Rm 是一個(gè)從n維歐氏空間轉(zhuǎn)換到m維歐氏空間的函數(shù)。這個(gè)函數(shù)由m個(gè)實(shí)函數(shù)組成:。這些函數(shù)的偏導(dǎo)數(shù)(如果存在)可以組成一個(gè)m行n列的矩陣(m by n)，這就是所謂的雅可比矩陣：

總結(jié)一下,

a)?如果 f(x) 是一個(gè)標(biāo)量函數(shù)，那么雅克比矩陣是一個(gè)向量，等于 f(x) 的梯度， Hesse 矩陣是一個(gè)二維矩陣。如果 f(x) 是一個(gè)向量值函數(shù)，那么Jacobi 矩陣是一個(gè)二維矩陣，Hesse 矩陣是一個(gè)三維矩陣。

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 次迭代我們選擇移動(dòng)的方向，在 steepest descent 中，移動(dòng)的方向設(shè)定為梯度的負(fù)方向，是第 k 次迭代用 line search 方法選擇移動(dòng)的距離，每次移動(dòng)的距離系數(shù)可以相同，也可以不同，有時(shí)候我們也叫學(xué)習(xí)率（learning rate）。 在數(shù)學(xué)上，移動(dòng)的距離可以通過 line search 令導(dǎo)數(shù)為零找到該方向上的最小值，但是在實(shí)際編程的過程中，這樣計(jì)算的代價(jià)太大，我們一般可以將它設(shè)定位一個(gè)常量。考慮一個(gè)包含三個(gè)變量的函數(shù)，計(jì)算梯度得到。設(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ù)是一個(gè)凸優(yōu)化問題，那么局部最優(yōu)解就是全局最優(yōu)解，理想的優(yōu)化效果如下圖，值得注意一點(diǎn)的是，每一次迭代的移動(dòng)方向都與出發(fā)點(diǎn)的等高線垂直：

需要指出的是，在某些情況下，最速下降法存在鋸齒現(xiàn)象（ zig-zagging）將會(huì)導(dǎo)致收斂速度變慢:

粗略來講，在二次函數(shù)中，橢球面的形狀受 hesse 矩陣的條件數(shù)影響，長軸與短軸對應(yīng)矩陣的最小特征值和最大特征值的方向，其大小與特征值的平方根成反比，最大特征值與最小特征值相差越大，橢球面越扁，那么優(yōu)化路徑需要走很大的彎路，計(jì)算效率很低。

2) Newton’s method

在最速下降法中，我們看到，該方法主要利用的是目標(biāo)函數(shù)的局部性質(zhì)，具有一定的“盲目性”。牛頓法則是利用局部的一階和二階偏導(dǎo)信息，推測整個(gè)目標(biāo)函數(shù)的形狀，進(jìn)而可以求得出近似函數(shù)的全局最小值，然后將當(dāng)前的最小值設(shè)定近似函數(shù)的最小值。相比最速下降法，牛頓法帶有一定對全局的預(yù)測性，收斂性質(zhì)也更優(yōu)良。牛頓法的主要推導(dǎo)過程如下：

第一步，利用 Taylor 級(jí)數(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)解。

牛頓法主要存在的問題是：

Hesse 矩陣不可逆時(shí)無法計(jì)算

矩陣的逆計(jì)算復(fù)雜為 n 的立方，當(dāng)問題規(guī)模比較大時(shí)，計(jì)算量很大，解決的辦法是采用擬牛頓法如 BFGS, L-BFGS, DFP, Broyden’s Algorithm 進(jìn)行近似。

如果初始值離局部極小值太遠(yuǎn)，Taylor 展開并不能對原函數(shù)進(jìn)行良好的近似

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è)一個(gè)可以信賴的最大位移 s（牛頓法里面 s 為無窮大），然后在以當(dāng)前點(diǎn)為中心，以 s 為半徑的區(qū)域內(nèi)，通過尋找目標(biāo)函數(shù)的一個(gè)近似函數(shù)（二次的）的最優(yōu)點(diǎn)，來求解得到真正的位移。在得到了位移之后，再計(jì)算目標(biāo)函數(shù)值，如果其使目標(biāo)函數(shù)值的下降滿足了一定條件，那么就說明這個(gè)位移是可靠的，則繼續(xù)按此規(guī)則迭代計(jì)算下去；如果其不能使目標(biāo)函數(shù)值的下降滿足一定的條件，則應(yīng)減小信賴域的范圍，再重新求解。

LMA 最早提出是用來解決最小二乘法曲線擬合的優(yōu)化問題的，對于隨機(jī)初始化的已知參數(shù) beta， 求得的目標(biāo)值為：

對擬合曲線函數(shù)進(jìn)行一階 Jacobi 矩陣的近似：

進(jìn)而推測出 S 函數(shù)的周邊信息：

位移是多少時(shí)得到 S 函數(shù)的最小值呢？通過幾何的概念，當(dāng)殘差垂直于 J 矩陣的 span 空間時(shí)， S 取得最小（至于為什么？請參考之前博客的最后一部分)

我們將這個(gè)公式略加修改，加入阻尼系數(shù)得到：

就是萊文貝格－馬夸特方法。這種方法只計(jì)算了一階偏導(dǎo)，而且不是目標(biāo)函數(shù)的 Jacobia 矩陣，而是擬合函數(shù)的 Jacobia 矩陣。當(dāng)大的時(shí)候可信域小，這種算法會(huì)接近最速下降法，小的時(shí)候可信域大，會(huì)接近高斯-牛頓方法。

算法過程如下：

給定一個(gè)初識(shí)值 x0

當(dāng)并且沒有到達(dá)最大迭代次數(shù)時(shí)

重復(fù)執(zhí)行:

• 算出移動(dòng)向量
• 計(jì)算更新值：
• 計(jì)算目標(biāo)函數(shù)真實(shí)減少量與預(yù)測減少量的比率
• if，接受更新值
• else if，說明近似效果很好，接受更新值，擴(kuò)大可信域（即減小阻尼系數(shù)）
• else: 目標(biāo)函數(shù)在變大，拒絕更新值，減小可信域（即增加阻尼系數(shù)）

直到達(dá)到最大迭代次數(shù)

維基百科在介紹 Gradient descent 時(shí)用包含了細(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 對局部極小值很敏感，維基百科舉了一個(gè)二乘法曲線擬合的例子，當(dāng)使用不同的初始值時(shí)，得到的結(jié)果差距很大，我這里也有 python 代碼，就不細(xì)說了。

4) Conjugate Gradients

共軛梯度法也是優(yōu)化模型經(jīng)常經(jīng)常要用到的一個(gè)方法，背后的數(shù)學(xué)公式和原理稍微復(fù)雜一些，光這一個(gè)優(yōu)化方法就可以寫一篇很長的博文了，所以這里并不打算詳細(xì)講解每一步的推導(dǎo)過程，只簡單寫一下算法的實(shí)現(xiàn)過程。與最速梯度下降的不同，共軛梯度的優(yōu)點(diǎn)主要體現(xiàn)在選擇搜索方向上。在了解共軛梯度法之前，我們首先簡單了解一下共軛方向：

共軛方向和馬氏距離的定義有類似之處，他們都考慮了全局的數(shù)據(jù)分布。如上圖，d(1) 方向與二次函數(shù)的等值線相切，d(1) 的共軛方向 d(2) 則指向橢圓的中心。所以對于二維的二次函數(shù)，如果在兩個(gè)共軛方向上進(jìn)行一維搜索，經(jīng)過兩次迭代必然達(dá)到最小點(diǎn)。前面我們說過，等值線橢圓的形狀由 Hesse 矩陣決定，那么，上圖的兩個(gè)方向關(guān)于 Hessen 矩陣正交，共軛方向的定義如下：

如果橢圓是一個(gè)正圓， Hessen 矩陣是一個(gè)單位矩陣，上面等價(jià)于歐幾里得空間中的正交。

在優(yōu)化過程中，如果我們確定了移動(dòng)方向（GD：垂直于等值線，CG：共軛方向），然后在該方向上搜索極小值點(diǎn)（恰好與該處的等值線相切），然后移動(dòng)到最小值點(diǎn)，重復(fù)以上過程，那么 Gradient Descent 和 Conjugate gradient descent 的優(yōu)化過程可以用下圖的綠線與紅線表示：

講了這么多，共軛梯度算法究竟是如何算的呢？

給定一個(gè)出發(fā)點(diǎn) x0 和一個(gè)停止參數(shù) e, 第一次移動(dòng)方向?yàn)?strong>最速下降方向:

while:

• 用 Newton-Raphson 迭代計(jì)算移動(dòng)距離，以便在該搜索方向移動(dòng)到極小，公式就不寫了，具體思路就是利用一階梯度的信息向極小值點(diǎn)跳躍搜索
• 移動(dòng)當(dāng)前的優(yōu)化解 x：
• 用 Gram-Schmidt 方法構(gòu)造下一個(gè)共軛方向，即, 按照的確定公式又可以分為 FR 方法和 PR 和 HS 等。

在很多的資料中，介紹共軛梯度法都舉了一個(gè)求線性方程組 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)

總結(jié)

以上是生活随笔為你收集整理的几种常用的优化方法梯度下降法、牛顿法、）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。