深入機器學習系列17-BFGS & L-BFGS

05-04

轉載請註明出處，該文章的官方來源：

最優化演算法 | Teaching ML

1 牛頓法

設f(x)是二次可微實函數，又設 $x^{(k)}$ 是f(x)一個極小點的估計，我們把f(x)在 $x^{(k)}$ 處展開成Taylor級數，並取二階近似。

上式中最後一項的中間部分表示f(x)在 $x^{(k)}$ 處的Hesse矩陣。對上式求導並令其等於0，可以的到下式：

設Hesse矩陣可逆，由上式可以得到牛頓法的迭代公式如下**(1.1)**

值得注意，當初始點遠離極小點時，牛頓法可能不收斂。原因之一是牛頓方向不一定是下降方向，經迭代，目標函數可能上升。此外，即使目標函數下降，得到的點也不一定是沿牛頓方向最好的點或極小點。因此，我們在牛頓方向上增加一維搜索，提出阻尼牛頓法。其迭代公式是**(1.2)**：

其中，lambda是由一維搜索（參考文獻【1】了解一維搜索）得到的步長，即滿足

2 擬牛頓法

2.1 擬牛頓條件

前面介紹了牛頓法，它的突出優點是收斂很快，但是運用牛頓法需要計算二階偏導數，而且目標函數的Hesse矩陣可能非正定。為了克服牛頓法的缺點，人們提出了擬牛頓法，它的基本思想是用不包含二階導數的矩陣近似牛頓法中的Hesse矩陣的逆矩陣。由於構造近似矩陣的方法不同，因而出現不同的擬牛頓法。

下面分析怎樣構造近似矩陣並用它取代牛頓法中的Hesse矩陣的逆。上文**(1.2)**已經給出了牛頓法的迭代公式，為了構造Hesse矩陣逆矩陣的近似矩陣 $H_{(k)}$ ，需要先分析該逆矩陣與一階導數的關係。

設在第k次迭代之後，得到 $x^{(k+1)}$ ，我們將目標函數f(x)在點 $x^{(k+1)}$ 展開成Taylor級數，並取二階近似，得到

由此可知，在 $x^{(k+1)}$ 附近有，

記

則有

又設Hesse矩陣可逆，那麼上式可以寫為如下形式。

這樣，計算出p和q之後，就可以通過上面的式子估計Hesse矩陣的逆矩陣。因此，為了用不包含二階導數的矩陣 $H_{(k+1)}$ 取代牛頓法中Hesse矩陣的逆矩陣，有理由令 $H_{(k+1)}$ 滿足公式**(2.1)**：

公式**(2.1)**稱為擬牛頓條件。

2.2 秩1校正

當Hesse矩陣的逆矩陣是對稱正定矩陣時，滿足擬牛頓條件的矩陣 $H_{(k)}$ 也應該是對稱正定矩陣。構造這樣近似矩陣的一般策略是， $H_{(1)}$ 取為任意一個n階對稱正定矩陣，通常選擇n階單位矩陣I，然後通過修正 $H_{(k)}$ 給定 $H_{(k+1)}$ 。令，

秩1校正公式寫為如下公式**(2.2)**形式。

2.3 DFP演算法

著名的DFP方法是Davidon首先提出，後來又被Feltcher和Powell改進的演算法，又稱為變尺度法。在這種方法中，定義校正矩陣為公式**(2.3)**

那麼得到的滿足擬牛頓條件的DFP公式如下**(2.4)**

查看文獻【1】，了解DFP演算法的計算步驟。

2.4 BFGS演算法

前面利用擬牛頓條件**(2.1)推導出了DFP公式(2.4)。下面我們用不含二階導數的矩陣 $B_{(k+1)}$ 近似Hesse矩陣，從而給出另一種形式的擬牛頓條件(2.5)**:

將公式**(2.1)的H換為B，p和q互換正好可以得到公式(2.5)。所以我們可以得到B的修正公式(2.6)**:

這個公式稱關於矩陣B的BFGS修正公式，也稱為DFP公式的對偶公式。設 $B_{(k+1)}$ 可逆，由公式**(2.1)以及(2.5)**可以推出：

這樣可以得到關於H的BFGS公式為下面的公式**(2.7)**:

這個重要公式是由Broyden,Fletcher,Goldfard和Shanno於1970年提出的，所以簡稱為BFGS。數值計算經驗表明，它比DFP公式還好，因此目前得到廣泛應用。

2.5 L-BFGS（限制內存BFGS）演算法

在BFGS演算法中，仍然有缺陷，比如當優化問題規模很大時，矩陣的存儲和計算將變得不可行。為了解決這個問題，就有了L-BFGS演算法。L-BFGS即Limited-memory BFGS。 L-BFGS的基本思想是只保存最近的m次迭代信息，從而大大減少數據的存儲空間。對照BFGS，重新整理一下公式：

之前的BFGS演算法有如下公式**(2.8)**

那麼同樣有

將該式子帶入到公式**(2.8)**中，可以推導出如下公式

假設當前迭代為k，只保存最近的m次迭代信息，按照上面的方式迭代m次，可以得到如下的公式**(2.9)**

上面迭代的最終目的就是找到k次迭代的可行方向，即

為了求可行方向r，可以使用two-loop recursion演算法來求。該演算法的計算過程如下，演算法中出現的y即上文中提到的t：

演算法L-BFGS的步驟如下所示。

2.6 OWL-QN演算法

2.6.1 L1 正則化

在機器學習演算法中，使用損失函數作為最小化誤差，而最小化誤差是為了讓我們的模型擬合我們的訓練數據，此時，若參數過分擬合我們的訓練數據就會有過擬合的問題。正則化參數的目的就是為了防止我們的模型過分擬合訓練數據。此時，我們會在損失項之後加上正則化項以約束模型中的參數：

$J(x) = l(x) + r(x)$

公式右邊的第一項是損失函數，用來衡量當訓練出現偏差時的損失，可以是任意可微凸函數（如果是非凸函數該演算法只保證找到局部最優解）。第二項是正則化項。用來對模型空間進行限制，從而得到一個更「簡單」的模型。

根據對模型參數所服從的概率分布的假設的不同，常用的正則化一般有L2正則化（模型參數服從Gaussian分布）、L1正則化（模型參數服從Laplace分布）以及它們的組合形式。

L1正則化的形式如下

$J(x) = l(x) + C ||x||_{1}$

L2正則化的形式如下

$J(x) = l(x) + C ||x||_{2}$

L1正則化和L2正則化之間的一個最大區別在於前者可以產生稀疏解，這使它同時具有了特徵選擇的能力，此外，稀疏的特徵權重更具有解釋意義。如下圖：

圖左側是L2正則，右側為L1正則。當模型中只有兩個參數，即 $w_1$ 和 $w_2$ 時，L2正則的約束空間是一個圓，而L1正則的約束空間為一個正方形，這樣，基於L1正則的約束會產生稀疏解，即圖中某一維( $w_2$ )為0。而L2正則只是將參數約束在接近0的很小的區間里，而不會正好為0(不排除有0的情況)。對於L1正則產生的稀疏解有很多的好處，如可以起到特徵選擇的作用，因為有些維的係數為0，說明這些維對於模型的作用很小。

這裡有一個問題是，L1正則化項不可微，所以無法像求L-BFGS那樣去求。微軟提出了OWL-QN(Orthant-Wise Limited-Memory Quasi-Newton)演算法，該演算法是基於L-BFGS演算法的可用於求解L1正則的演算法。簡單來講，OWL-QN演算法是指假定變數的象限確定的條件下使用L-BFGS演算法來更新，同時，使得更新前後變數在同一個象限中(使用映射來滿足條件)。

2.6.2 OWL-QN演算法的具體過程

1 次微分

設 $f:I ightarrow R$ 是一個實變數凸函數，定義在實數軸上的開區間內。這種函數不一定是處處可導的，例如絕對值函數 $f(x)=|x|$ 。但是，從下面的圖中可以看出（也可以嚴格地證明），對於定義域中的任何 $x_0$ ，我們總可以作出一條直線，它通過點( $x_0$ , $f(x_0)$ )，並且要麼接觸f的圖像，要麼在它的下方。這條直線的斜率稱為函數的次導數。推廣到多元函數就叫做次梯度。

凸函數 $f:I ightarrow R$ 在點 $x_0$ 的次導數，是實數c使得：

對於所有I內的x。我們可以證明，在點 $x_0$ 的次導數的集合是一個非空閉區間 $[a, b]$ ，其中a和b是單側極限。

它們一定存在，且滿足 $a leqslant b$ 。所有次導數的集合 $[a, b]$ 稱為函數f在 $x_0$ 的次微分。

2 偽梯度

利用次梯度的概念推廣了梯度，定義了一個符合上述原則的偽梯度，求一維搜索的可行方向時用偽梯度來代替L-BFGS中的梯度。

其中

我們要如何理解這個偽梯度呢？對於不是處處可導的凸函數，可以分為下圖所示的三種情況。

左側極限小於0：

右側極限大於0：

其它情況：

結合上面的三幅圖表示的三種情況以及偽梯度函數公式，我們可以知道，偽梯度函數保證了在 $x_0$ 處取得的方嚮導數是最小的。

3 映射

有了函數的下降的方向，接下來必須對變數的所屬象限進行限制，目的是使得更新前後變數在同一個象限中，定義函數： $pi: mathbb{R}^{n} ightarrow mathbb{R}^{n}$

上述函數 $pi$ 直觀的解釋是若x和y在同一象限則取x，若兩者不在同一象限中，則取0。

4 線搜索

上述的映射是防止更新後的變數的坐標超出象限，而對坐標進行的一個約束，具體的約束的形式如下：

其中 $x^{k} + alpha p _{k}$ 是更新公式， $zeta$ 表示 $x^k$ 所在的象限， $p^k$ 表示偽梯度下降的方向，它們具體的形式如下：

上面的公式中， $v^k$ 為負偽梯度方向， $d^k = H_{k}v^{k}$ 。

選擇 $alpha$ 的方式有很多種，在OWL-QN中，使用了backtracking line search的一種變種。選擇常數 $eta, gamma subset (0,1)$ ，對於 $n=0,1,2,...$ ，使得 $alpha = eta^{n}$ 滿足：

5 演算法流程

與L-BFGS相比，第一步用偽梯度代替梯度，第二、三步要求一維搜索不跨象限，也就是迭代前的點與迭代後的點處於同一象限，第四步要求估計Hessian矩陣時依然使用損失函數的梯度。

3 源碼解析

3.1 BreezeLBFGS

spark Ml調用breeze中實現的BreezeLBFGS來解最優化問題。

val optimizer = new BreezeLBFGS[BDV[Double]]($(maxIter), 10, $(tol))val states = optimizer.iterations(new CachedDiffFunction(costFun), initialWeights.toBreeze.toDenseVector)

下面重點分析lbfgs.iterations的實現。

def iterations(f: DF, init: T): Iterator[State] = { val adjustedFun = adjustFunction(f) infiniteIterations(f, initialState(adjustedFun, init)).takeUpToWhere(_.converged)}//調用infiniteIterations，其中State是一個樣本類def infiniteIterations(f: DF, state: State): Iterator[State] = { var failedOnce = false val adjustedFun = adjustFunction(f) //無限迭代 Iterator.iterate(state) { state => try { //1 選擇梯度下降方向 val dir = chooseDescentDirection(state, adjustedFun) //2 計算步長 val stepSize = determineStepSize(state, adjustedFun, dir) //3 更新權重 val x = takeStep(state,dir,stepSize) //4 利用CostFun.calculate計算損失值和梯度 val (value,grad) = calculateObjective(adjustedFun, x, state.history) val (adjValue,adjGrad) = adjust(x,grad,value) val oneOffImprovement = (state.adjustedValue - adjValue)/(state.adjustedValue.abs max adjValue.abs max 1E-6 * state.initialAdjVal.abs) //5 計算s和t val history = updateHistory(x,grad,value, adjustedFun, state) //6 只保存m個需要的s和t val newAverage = updateFValWindow(state, adjValue) failedOnce = false var s = State(x,value,grad,adjValue,adjGrad,state.iter + 1, state.initialAdjVal, history, newAverage, 0) val improvementFailure = (state.fVals.length >= minImprovementWindow && state.fVals.nonEmpty && state.fVals.last > state.fVals.head * (1-improvementTol)) if(improvementFailure) s = s.copy(fVals = IndexedSeq.empty, numImprovementFailures = state.numImprovementFailures + 1) s } catch { case x: FirstOrderException if !failedOnce => failedOnce = true logger.error("Failure! Resetting history: " + x) state.copy(history = initialHistory(adjustedFun, state.x)) case x: FirstOrderException => logger.error("Failure again! Giving up and returning. Maybe the objective is just poorly behaved?") state.copy(searchFailed = true) } } }

看上面的代碼注釋，它的流程可以分五步來分析。

3.1.1 選擇梯度下降方向

protected def chooseDescentDirection(state: State, fn: DiffFunction[T]):T = { state.history * state.grad}

這裡的*是重寫的方法，它的實現如下：

def *(grad: T) = { val diag = if(historyLength > 0) { val prevStep = memStep.head val prevGradStep = memGradDelta.head val sy = prevStep dot prevGradStep val yy = prevGradStep dot prevGradStep if(sy < 0 || sy.isNaN) throw new NaNHistory sy/yy } else { 1.0 } val dir = space.copy(grad) val as = new Array[Double](m) val rho = new Array[Double](m) //第一次遞歸 for(i <- 0 until historyLength) { rho(i) = (memStep(i) dot memGradDelta(i)) as(i) = (memStep(i) dot dir)/rho(i) if(as(i).isNaN) { throw new NaNHistory } axpy(-as(i), memGradDelta(i), dir) } dir *= diag //第二次遞歸 for(i <- (historyLength - 1) to 0 by (-1)) { val beta = (memGradDelta(i) dot dir)/rho(i) axpy(as(i) - beta, memStep(i), dir) } dir *= -1.0 dir } }

非常明顯，該方法就是實現了上文提到的two-loop recursion演算法。

3.1.2 計算步長

protected def determineStepSize(state: State, f: DiffFunction[T], dir: T) = { val x = state.x val grad = state.grad val ff = LineSearch.functionFromSearchDirection(f, x, dir) val search = new StrongWolfeLineSearch(maxZoomIter = 10, maxLineSearchIter = 10) // TODO: Need good default values here. val alpha = search.minimize(ff, if(state.iter == 0.0) 1.0/norm(dir) else 1.0) if(alpha * norm(grad) < 1E-10) throw new StepSizeUnderflow alpha }

這一步對應L-BFGS的步驟的Step 5，通過一維搜索計算步長。

3.1.3 更新權重

protected def takeStep(state: State, dir: T, stepSize: Double) = state.x + dir * stepSize

這一步對應L-BFGS的步驟的Step 5，更新權重。

3.1.4 計算損失值和梯度

protected def calculateObjective(f: DF, x: T, history: History): (Double, T) = { f.calculate(x) }

這一步對應L-BFGS的步驟的Step 7，使用傳人的CostFun.calculate方法計算梯度和損失值。並計算出s和t。

3.1.5 計算s和t，並更新history

//計算s和tprotected def updateHistory(newX: T, newGrad: T, newVal: Double, f: DiffFunction[T], oldState: State): History = { oldState.history.updated(newX - oldState.x, newGrad :- oldState.grad)}//添加新的s和t，並刪除過期的s和tprotected def updateFValWindow(oldState: State, newAdjVal: Double):IndexedSeq[Double] = { val interm = oldState.fVals :+ newAdjVal if(interm.length > minImprovementWindow) interm.drop(1) else interm }

3.2 BreezeOWLQN

BreezeOWLQN的實現與BreezeLBFGS的實現主要有下面一些不同點。

3.2.1 選擇梯度下降方向

override protected def chooseDescentDirection(state: State, fn: DiffFunction[T]) = { val descentDir = super.chooseDescentDirection(state.copy(grad = state.adjustedGradient), fn) // The original paper requires that the descent direction be corrected to be // in the same directional (within the same hypercube) as the adjusted gradient for proof. // Although this doesnt seem to affect the outcome that much in most of cases, there are some cases // where the algorithm wont converge (confirmed with the author, Galen Andrew). val correctedDir = space.zipMapValues.map(descentDir, state.adjustedGradient, { case (d, g) => if (d * g < 0) d else 0.0 }) correctedDir }

此處調用了BreezeLBFGS的chooseDescentDirection方法選擇梯度下降的方向，然後調整該下降方向為正確的方向（方向必須一致）。

3.2.2 計算步長 $alpha$

override protected def determineStepSize(state: State, f: DiffFunction[T], dir: T) = { val iter = state.iter val normGradInDir = { val possibleNorm = dir dot state.grad possibleNorm } val ff = new DiffFunction[Double] { def calculate(alpha: Double) = { val newX = takeStep(state, dir, alpha) val (v, newG) = f.calculate(newX) // 計算梯度 val (adjv, adjgrad) = adjust(newX, newG, v) // 調整梯度 adjv -> (adjgrad dot dir) } } val search = new BacktrackingLineSearch(state.value, shrinkStep= if(iter < 1) 0.1 else 0.5) val alpha = search.minimize(ff, if(iter < 1) .5/norm(state.grad) else 1.0) alpha }

takeStep方法用於更新參數。

// projects x to be on the same orthant as y // this basically requires that x_i = x_i if sign(x_i) == sign(y_i), and 0 otherwise. override protected def takeStep(state: State, dir: T, stepSize: Double) = { val stepped = state.x + dir * stepSize val orthant = computeOrthant(state.x, state.adjustedGradient) space.zipMapValues.map(stepped, orthant, { case (v, ov) => v * I(math.signum(v) == math.signum(ov)) }) }

calculate方法用於計算梯度，adjust方法用於調整梯度。

// Adds in the regularization stuff to the gradient override protected def adjust(newX: T, newGrad: T, newVal: Double): (Double, T) = { var adjValue = newVal val res = space.zipMapKeyValues.mapActive(newX, newGrad, {case (i, xv, v) => val l1regValue = l1reg(i) require(l1regValue >= 0.0) if(l1regValue == 0.0) { v } else { adjValue += Math.abs(l1regValue * xv) xv match { case 0.0 => { val delta_+ = v + l1regValue //計算左導數 val delta_- = v - l1regValue //計算右導數 if (delta_- > 0) delta_- else if (delta_+ < 0) delta_+ else 0.0 } case _ => v + math.signum(xv) * l1regValue } } }) adjValue -> res }

參考文獻

【1】陳寶林，最優化理論和演算法

【2】[Updating Quasi-Newton Matrices with Limited Storage](docs/Updating Quasi-Newton Matrices with Limited Storage.pdf)

【3】[On the Limited Memory BFGS Method for Large Scale Optimization](docs/On the Limited Memory BFGS Method for Large Scale Optimization.pdf)

【4】L-BFGS演算法

【5】BFGS演算法

【6】邏輯回歸模型及LBFGS的Sherman Morrison(SM) 公式推導

【7】Scalable Training of L1-Regularized Log-Linear Models