非凸優化之非限制優化

無限制優化

本篇文檔主要介紹無限制優化問題，同時初步介紹解該類問題目前常用的一種算法即

Quasi-Newton Method (拟牛頓法)。在介紹無限制優化問題之前，我們首先會從直覺上引入無限制優化的概念，并在此基礎上引入解這類問題的兩個重要概念：步長和方向。由步長的選擇引入重要概念 line search，由方向的選擇引入重要概念 Quasi-Newton Method。是以本篇介紹文檔主要分為以下幾個部分：無限制優化問題引入，Line Search，Quasi-Newton Method 和算法總結。

1．無限制最優化

對無限制優化不熟悉的讀者也許要問，什麼是無限制優化。這裡以一個例子來說明該問題。

上圖所示為一進制函數 f(x)的圖像，無限制優化問題，即不對定義域或值域做任何限制的情況下，求解函數 f(x)的小值，上面顯示兩個小值點：一個為全局小值點，另一個為局部小值點。受限于算法複雜度等問題，目前大部分無限制優化算法隻能保證求取局部小值點。這時讀者不免要問，既然隻能求取局部小值點，那為什麼這類算法還能應用呢。這是因為實際應用中，許多情形被抽象為函數形式後均為凸函數，對于凸函數來說局部小值點即為全局小值點，是以隻要能求得這類函數的一個小值點，該點一定為全局小值點。

了解了上面的無限制優化問題之後，我們就可以開始介紹無限制優化的求解過程

了，對于無限制優化的求解首先我們需要選擇一個初始點 x_0，如下所示：

初始點選擇好之後，就可以按照各種不同的無限制優化求解算法，求解小值點了。求解過程中主要涉及兩個概念，即從初始點開始沿“哪個方向”以及“走多遠”到達下一個點處。所謂“走多遠”即之前提的“步長”的概念，“哪個方向”即方向概念。

2．Line Search

Line search 主要用于解決之前提到的步長的概念，即方向确定好之後，需要确定從目前點 x_k 沿着該方向走多遠，以便走到下一個合适的點 x_k+1。若用 p_k 代表從第 k 個點走向第 k+1 點的方向，x_k 代表目前點，x_k+1 代表下一個點，a_k 代表步長，則存在如下的等式： x_k+1 = x_k + a_k * p_k (1)

這裡簡要介紹一下 p_k，大部分 line search 方法要求 p_k 為下降方向，即從目前點沿着 p_k 方向移動後導緻函數值減少。由于目标是求取一個函數的小值，是以優的情況是求取沿 p_k 方向滿足 f(x_k+1)為全局小的 a_k 值，可用下式表示為：

Ø (a_k) = f(x_k+1) = f(x_k + a_k * p_k) (2)

由于直接求取滿足Ø (a_k)為全局小值的 a_k 涉及到大量 f(x_k + a_k * p_k)的計算，若從求導角度計算小值，還會涉及到▽f_k+1的計算，計算量較大。是以從計算量角度考慮，可以采用如下較為折中的政策。 方向确定好之後，每一步的 line Search 主要涉及兩個問題：1）什麼樣的 a_k 是合理的 2）

如何選擇 a_k 的長度。下面将沿這兩個面展開讨論，首先讨論“什麼樣的 a_k 是合理的”，确定了該問題之後，我們就可以在此基礎上選擇 a_k 了。

2.1 a_k 合理性讨論 如下将要讨論關于a_k需要滿足的兩個條件，當a_k滿足這兩個條件後，就可以認為從x_k 點移動到x_k+1 點的步長已經确定下來了。第一個條件為sufficient decrease condition，從直

觀角度來看，該條件主要要用保證x_k+1 點的函數值要小于x_k點的函數值，滿足該條件後，才有全局收斂[1]的可能性。第二個條件為curvature condition，從直覺角度來看，該條件主要用于保證x_k點經過步長a_k的移動到達x_k+1 後，▽f_k+1小于▽f_k。

2.1.1 sufficient decrease condition

a_k 的選擇一定要使得函數值滿足 sufficient decrease condition，該條件可以用如下不等式描述：

f(x_k+1) ≤ f(x_k) + c_1 a_ k ▽f_k T * p_k (3)

将公式（1）代入上式，可得：

f(x_k + a_k p_k) ≤ f(x_k) + c_1 a_ k ▽f_k T p_k (4) 這裡有必要對上面的不等式做一些解釋：

f(x_k) 代表函數在第 x_k 點的值

▽f_k代表函數在第 x_k 點的梯度

p_k 代表從第 x_k 點的走到 x_k+1 點的方向

a_ k 代表從第 x_k 點沿着 p_k 方向走到 x_k+1 點步長

c_1[2]為常量，需滿足 0< c_1 < 1，一般取c_1 為 1E-4 當 p_k 為函數下降方向時，有：

▽f_k * p_k  <  0                                               (5) 是以不等式 4，即要求：            

f(x_k+1) ＜ f(x_k) (6) 從圖形角度看，函數位于第 k 點時，以上各參數中隻有 a_k 為變量，其他均為常量，因

此(4)可以用以下不等式重新描述為：

Ø(a_k) ≤ l(a_k) (7)

其中：

Ø (a_k) = f(x_k + a_k p_k) l(a_k) = f(x_k) + c_1 a_ k ▽f_k T p_k

以下為不等式（7）的圖形化表示：

是以隻要步長 a_k 的選擇使得函數 Ø(a_k)位于 acceptable 區間，就滿足 sufficient decrease condition。

2.1.2 curvature condition

a_k 的選擇一定要使得函數梯度值滿足 curvature condition，該條件可以用如下不等式描

述：

▽f(x_k + a_k * p_k) T *p_k ≥ c_2 * ▽f(x_k) T * p_k                         (8) 即為 
   ▽f _k+1 T *p_k ≥ c_2 * ▽f_k T * p_k                                           (9) 這裡有必要對上面的不等式做一些解釋：            

▽f_k+1代表函數在第 x_k+1 點的梯度

p_k 代表從第 k 點的走到 k+1 點的方向

c_2 為常量，需滿足 0< c_1 < c_2 < 1，一般取 c_2 為 0.9 當 p_k 為函數下降方向時，有：

▽f_k * p_k  <  0             

是以不等式 9，即要求：

▽f _k+1 ≥  c_2 * ▽f_k                                       (10) 從圖形角度看，不等式10即要求函數在第x_k+1點的變化速度要低于x_k點的變化速度，這一點可以從這兩點處的梯度處看出，如下圖所示。 
           

2.1.3 Wolfe conditions 所謂 Wolfe conditions 即 sufficient decrease condition 和 curvature condition 的綜合，即 a_k

需要同時滿足如下兩個條件：

f(x_k + a_k p_k) ≤ f(x_k) + c_1 a_ k ▽f_k T p_k (11)

▽f(x_k + a_k p_k) T p_k ≥ c_2 ▽f(x_k) T p_k (12)

圖形化表示後，如下圖所示：

實際應用中，常常會用到由 Wolfe conditions 引申出的 strong Wolfe conditions，即 a_k 需要

同時滿足如下兩個條件：

f(x_k + a_k p_k) ≤ f(x_k) + c_1 a_ k ▽f_k T p_k (13) |▽f(x_k + a_k p_k) T p_k| ≤ c_2 |▽f(x_k) T p_k| (14) 與 Wolfe condtions 唯一的差別是 strong Wolfe condtions 避免了▽f(x_k + a_k * p_k)取較大的正值情況。

2.2 a_k 步長的選擇

了解了 a_k 的合理性之後，就相當于獲得了标尺，在此基礎上我們可以選擇合适的政策來求取 a_k。所有的 line search 過程在計算每一步的 a_k 時，均需要提供一個初始點 a_0，然後再此基礎上生成一系列的{a_i}，直到 a_i 滿足 2.1 節所規定的條件為止，此時該 a_k 即被确定為 a_i，或者未找到一個合适的 a_k。這裡我們僅介紹目前常用的政策平方插值和立方插值法。是以本節内容分為兩部分，2.2.1 節介紹選擇 a_k 常用的平方插值和立方插值法，

2.2.2 節介紹由 x_k 點到 x_k+1 點，方向确定為 p_k 後，步長 a_k 具體計算過程。

2.2.1 平方立方插值法 當給定一個初始步長 a_0，若該初始步長滿足 Wolfe conditions（或 strong Wolfe

conditions），則 a_k 被确定為 a_0，目前點 x_k 步長計算過程結束，否則，我們可以在此基礎上，利用已知的三個資訊 Ø(a_0)、Ø(0)、Ø’(0)，建構一個二次（平方）插值多項式來拟合Ø(a_k)。該二次插值多項式如下所示：

Ø_q(a) =( ( Ø(a_0) - Ø(0) - a_0Ø’(0) ) / a_0 2 )a2 + Ø’(0)*a + Ø(0) (15) 觀察上面的二次插值多項式可知，其滿足如下插值條件：

Ø_q(0) = Ø(0) Ø_q’(0) = Ø’(0) Ø_q(a_0) = Ø(a_0)

對二次插值多項式(15)求導并令其為零，即可獲得使得該多項式取得小值的 a 值，如下所示：

a_1 = -1 Ø’(0) a_0 2 / ( 2 ( Ø(a_0) - Ø(0) – Ø’(0)a_0 ) ) (16) 若 a_1 滿足 Wolfe conditions（或 strong Wolfe conditions），則 a_k 被确定為 a_1。否則在此基礎上建構一個三次（立方）插值多項式，并求得使得該多項式取小值的 a 值，該 a 值的計算公式如下所示：

a_i+ 1= a_i – (a_i – a_i-1)(Ø’(a_i) + d_2 – d_1) / (Ø’(a_i) – Ø’(a_i-1) + 2d_2) (17) d_1 = Ø’(a_i-1) + Ø’(a_i) – 3( Ø(a_i-1) - Ø(a_i) ) / ( a_i-1 – a_i) d_2 = sign(a_i - a_i-1) ( d_12 - Ø’(a_i-1) * Ø’(a_i) ) 1/2

若 a_i+1 滿足 Wolfe conditions（或 strong Wolfe conditions），則 a_k 被确定為 a_i+1，否則一直使用三次插值多項式進行插值拟合。并且選擇使用 a_i+1 對應的Ø (a_i+1)和Ø’( a_i +1) 替換 a_i-1 或 a_i 相應的值，一旦确定好 a_i-1 或 a_i 中的一種後，每次都替換對應的值，如替換 a_i-1，則每次都使用新的 a_i+1 對應的值替換 a_i-1，直到找到滿足 Wolfe conditions （或 strong Wolfe conditions）的 a_i+1 為止，此時 a_k 被确定為 a_i+1，或者沒有找到合适的 a_i+1。

這裡有必要簡略解釋一下為何隻使用到三次插值多項式，而沒有使用更高階的插值多項式。原因是三次插值多項式對函數某個點處的具體值有較好的拟合效果，同時又有較好的抗過拟合作用。

後有必要解釋一下初始步長 a_0 的選擇，對于 Newton 或 quasi-Newton methods 來說，初始步長 a_0 總是确定為 1，該選擇確定當滿足 Wolfe conditions（或 strong Wolfe conditions），我們總是選擇機關 1 步長，因為該步長使得 Newton 或 quasi-Newton methods 達到較快的收斂速度。計算第零步長時，将初始步長 a_0 使用如下公式确定：

a_0 = 1.0 / (p_0 T* p_0) 1/2 (18)

第 1 步及其以後的初始步長 a_0，直接設定為 1。

2.2.2 步長 a_k 計算

本節主要做一個總結，即綜合前面步長需要滿足的條件及步長疊代計算公式給出步長計算的具體過程。下面假設我們處于 x_k 點，是以要從該點選擇一個滿足 Wolfe (或 strong

Wolfe) conditions 的步長 a_k，以便走到下一點 x_k+1。是以我們選擇初始點 a_0 = 1，Newton

或 quasi-Newton method 中初始步長總是選擇為 1。是以有：

1) 初始化 a_x  a_y  0 Ø(a_x) 、 Ø(a_y) 、Ø’(a_x) 、 Ø’(a_y)、a_min、a_max 2) 初始化 a_i  a_0、Ø’(a_i)、Ø(a_i)

若Ø(a_i) > Ø(a_x)

說明步長 a_i 選擇的過大，是以使得Ø(a_i)滿足 Wolfe(或 strong Wolfe) conditions 的步長 a_k 應位于 a_x 和 a_i 之間，使用平方和立方插值法對 a_x 和 a_i 分别進行插值，求取兩個新的步長 a_quadratic 和 a_cubic，并取兩者之中和 a_x 較接近的步長，這裡假設為

a_quadratic，将新的步長設為 a_qudratic，同時進行以下操作： a_y  a_i

Ø(a_y)  Ø(a_i) Ø’(a_y)  Ø’(a_i) a_i+1  a_quadrati

同時在此情況下，說明 a_k 位于 a_x 和 a_y 之間

若Ø’(a_i)* Ø’(a_x) < 0

說明步長 a_i 選擇的過大，是以使得Ø(a_i)滿足 Wolfe(或 strong Wolfe) conditions 的步長 a_k 應位于 a_x 和 a_i 之間，使用平方和立方插值法對 a_x 和 a_i 分别進行插值，求取兩個

新的步長a_quadratic和a_cubic，并取兩者之中和a_i較接近的步長，這裡假設為a_quadratic，

将新的步長設為 a_qudratic，同時進行以下操作：

a_y  a_x

Ø(a_y)  Ø(a_x)

Ø’(a_y)  Ø’(a_x)

a_x  a_i

Ø(a_x)  Ø(a_i)

Ø’(a_x)  Ø’(a_i)

a_i+1  a_quadrati

同時在此情況下，說明 a_k 位于 a_x 和 a_y 之間

5) 若|Ø’(a_i)| <|Ø’(a_x) |

說明步長 a_i 選擇的過小，是以使得Ø(a_i)滿足 Wolfe(或 strong Wolfe) conditions 的步長 a_k 應位于 a_i 和 a_y 之間，使用平方和立方插值法對 a_x 和 a_i 分别進行插值，求取兩個

将新的步長設為 a_qudratic，同時進行以下操作： a_x  a_i

5) 若|Ø’(a_i)| ≥|Ø’(a_x) |

說明步長a_i選擇的過小，若之前已經确定出了a_k所屬的範圍a_x和a_y，則使得Ø( a_i ) 滿足 Wolfe(或 strong Wolfe) conditions 的步長 a_k 應位于 a_i 和 a_y 之間，使用立方插值法對 a_i 和 a_y 進行插值，求取新的步長 a_cubic，新的步長設為 a_cubic；否則若 a_x 小于 a_i，說明的新的步長不夠大，是以将新的步長設定為 a_max，若 a_x 大于等于 a_i，則說明新的步長不夠小，将其設為 a_min，同時進行以下操作： a_x  a_i

Ø’(a_x)  Ø’(a_i)

if 之前已經确定出了 a_k 所屬的範圍 a_x 和 a_y

a_i+1  a_cubic

else if a_x ＜ a_i a_i+1  a_max

else if a_x ≥ a_i a_i+1  a_min

6) 計算判斷若 a_i+1 使得以下兩式(strong Wolfe condition)均成立： f(x_k + a_i+1 p_k) ≤ f(x_k) + c_1 a_i+1 ▽f_k T p_k

|▽f(x_k + a_i+1 p_k) T p_k| ≤ c_2 |▽f(x_k) T p_k|

則找到合理的步長 a_k，将其設定為 a_i+1 a_k  a_i+1

則 x_k 點步長 a_k 計算結束。 否則轉到 2）繼續計算合理步長。

3．Quasi-Newton Method

在第 2 節中我們了解了步長的概念，以及從 x_k 走到 x_k+1 點使用 line search 方法計算步長的方法。不過我們在那裡忽略了一個重要的概念，即“方向”。從第 2 節，我們了解到從每一點 x_k 走到下一點 x_k+1 時，需要給出要走的“方向”，隻有“方向”确定好之後，才能在此基礎上應用 line search 方法找到對應的“步長”，是以在解決了“步長”計算問題之後，這裡我們将和大家一起了解一下每一步的“方向”如何确定。本節分為 2 大部分，首先我們通過 newton method 引入方向的概念，在此基礎上引入 quasi-newton method。然後引入 quasi-newton method 中的一種重要方法 BFGS method，并在 BFGS method 的基礎上介紹用于大規模計算的 LBFGS method 算法，同時以此結束本節的所有内容。            

3.1 Newton Method

我想我們還依稀記得，微積分中用于求一進制函數根的 Newton 法。該方法可用如下所示

的圖形來描述：

首先我們選擇一個初始點 x_0 并計算獲得其對應的 f(x_0)，然後該方法用曲線 y = f(x)在 (x_0,f(x_0))的切線近似該曲線，把切線和 x 軸的交點記作 x_1，點 x_1 通常 x_0 更接近所要求得根，同樣方法用點 x_1 處的切線近似該曲線，并求取切線和 x 軸的交點 x_2，一直疊代下去，直到找到滿足我們所需的充分接近真實跟的解為止。是以 Newton 法從第 k 次近似值 x_k 求得第 k+1 次近似值 x_k+1 即為求 x_k 點處的切線和 x 軸的交點。 切線方程為：

y – f(x_k) = ▽f (x_k)(x – x_k)

=> 0 – f(x_k) =▽ f(x_k)(x – x_k)

=> ▽f(x_k)x = ▽f(x_k)x_k – f(x_k) => x = x_k – f(x_k)/▽f (x_k) 是以 x_k+1 = x_k – f(x_k)/▽f(x_k) (19) 從上面使用 Newton Method 求函數根的過程可以發現，首先需要選擇一個初始點，并

在點處建構一個模型來近似該函數。

切線模型： M_k(x_k+p_k) = f(x_k) + ▽f (x_k) T *p_k

上面使用了相應點處的切線模型來近似函數，然後求取該近似模型的根以便求得更接近函數根的下一個點，該過程一直疊代下去，直到找到根為止。從上面建構每個點處的近似模型可以發現，該模型相對原函數來說簡化了很多，是以求解要容易一些。

現在來考慮求取函數的小值問題，方法類似。首先開始我們選擇一個初始點 x_0，并建構函數在該點處的一個近似模型，上面求函數根時，我們建構的近似模型為切線模型。這裡我們建構一個抛物線模型：

抛物線模型：M_k(x_k + p_k) = f(x_k) + ▽f(x_k) T p_k + 1/2 p_k T ▽2f(x_k) p_k

并求解該模型的梯度，同時令其為零，即：▽M_k(x+) = 0，在此基礎上求得 x+，該值

即使得近似模型取得小值的點。

對抛物線模型 M_k 求導并令其為零後，可得以下公式：

p_k = –▽2f(x_k) -1▽f(x_k) (20) 是以 x_k+1 = x_k –▽2f(x_k) -1▽f(x_k) (21)

從上式可見，使用Newton法時，每一步的方向為p_k，而步長為 1。從Newton法方向公式我們不難看出，若使用Newton法，則從每個小值近似點x_k走到下一個近似點x_k+1 的過程中将涉及函數Hessian矩陣▽▽f(x_k) (二階導)計算，而該Hessian矩陣無法保證在每個點處都是正定[3]的，對正定矩陣來說存在如下不等式：

p_k T ▽2f(x_k) p_k > 0 (22) 由于矩陣無法保證為正定矩陣，是以下式中

–▽f(x_k) T p_k = p_k T ▽2f(x_k) *p_k (23) p_k 無法保證總是為下降方向，即負梯度方向。 從上面的分析可見，Newton Method 面臨兩個問題：

Hessian 矩陣計算量較大

Hessian 矩陣無法保證在疊代的過程中始終處于正定狀态

為了解決這兩個問題，我們引入 Quasi-Newton Method。

3.2 Quasi-Newton Method

Quasi-Newton Method 每一步計算過程中僅涉及到函數值和函數梯度值計算，這樣有效避免了 Newton Method 中涉及到的 Hessian 矩陣計算問題。于 Newton Method 不同的是 Quasi-Newton Method 在每點處建構一個如下的近似模型：

M_k(x_k + p_k) = f(x_k) + ▽f(x_k) T p_k + 1/2 p_k T B_kp_k

從上面的近似模型我們可以看出，該模型用 B_k 代替了 Newton Method 中近似模型中涉及到的 Hessian 矩陣。是以 Quasi-Newton Method 中方向計算公式如下所示：

p_k = –B_k -1*▽f(x_k) (24)

這裡有必要解釋一下用于近似 Hessian 矩陣的 B_k 可行性，及一個指導性方案。根據

Taylor(泰勒)級數可知如下公式：

▽f(x_k + p_k) = ▽f(x_k) +▽2f(x_k) T p_k + ∫[▽2f(x_k+tp_k)-▽2f(x_k)]p_kdt

由于函數▽f(.)連續，是以上式可以表示為：

▽f(x_k+1) = ▽f(x_k) + ▽2f(x_k)*(x_k+1 – x_k) + o(||x_k+1 – x_k||)

▽2f(x_k)*(x_k+1 – x_k) ≈ ▽f(x_k+1) - ▽f(x_k) (25)

是以每一選擇Hessian矩陣的近似B_ k+1時，可以像式(24)那樣模仿真實的Hessian

矩陣的性質。得到下式：

B_k+1*s_k = y_k (26)

其中：

s_k = x_k+1 – x_k     y_k = ▽f(x_k+1) - ▽f(x_k)        (27）            

同時要求 B_k+1 為對稱正定矩陣。

3.2.1 BFGS Method

從 Quasi-Newton Method 方向公式 (24) 中，可以看到每一步計算方向的過程中均涉及到 B_k+1 矩陣求逆的問題，為了避免該計算，通過分析公式（26）可知，我們可以建構一個近似 H_k+1，該近視滿足如下方程：

H_k+1*y_k = s_k (28)

同時要求 H_k+1 為對稱正定矩陣。是以 BFGS Method 中，每個點處的方向由如下公式計算： p_k = –H_k*▽f(x_k) (29)

在此基礎上，BFGS 方向疊代公式如下所示：

H_k+1 = (1 – ρ_k s_ky_k T) H_k (1-ρ_k y_k s_k T) +ρ_k s_k s_k T (30)

其中ρ_k 為一個标量：

ρ_k  = 1 / (y_k T * s_k)            

有了上面(30)的 H_k 疊代公式後，還有一個問題就是初始的 H_0 如何計算，目前常用的方法是初始的 H_0 直接設為機關矩陣 I。是以 BFGS Method 用于解無限制優化的過程可以表示為如下過程：

選擇一個初始點 x_0，并選擇收斂判斷條件 ε> 0，及初始 H_0 (機關矩陣 I)

k  0

while ||▽f(x_k)|| > ε

計算從目前點 x_k 走到下一個點 x_k+1 的方向

p_k = –H_k*▽f(x_k)

采用 line search 政策計算步長 a_k

x_k+1 = x_k + a_k * p_k

s_k = x_k+1 – x_k y_k = ▽f(x_k+1) - ▽f(x_k)

H_k+1 = (1 – ρ_k s_ky_k T) H_k (1-ρ_k y_k s_k T) +ρ_k s_k s_k T

k  k+1

3.2.2 LBFGS Method

上一節所介紹的 BFGS Method 比較适合解決中小規模無限制優化問題，但是 BFGS 算法産生的 Hessian 近似矩陣 H_k 為 n * n 的，同時該矩陣非稀疏，是以當 n 的規模較大時将面臨兩個問題：

1） 存儲問題：n 規模較大時，nn 矩陣對記憶體的消耗将較大； 2） 計算問題：n 規模較大，同時 nn 矩陣非稀疏時，計算複雜度将較高；

為了解決以上問題，引申出了 Limited-Memory Quasi-Newton Method，目前使用較多的

LBFGS 算法即屬于該類算法。為了減少 H_k 矩陣的存儲，LBFGS 算法利用近幾代的

curvature 資訊來建構 Hessian 矩陣的近似。由 BFGS Method 我們知道：

x_k+1 = x_k + a_k H_k▽f(x_k)

其中 a_k 為步長，H_k 為 Hessian 矩陣的近似，可以通過如下疊代公式計算：

H_k+1 = V_k H_kV_k+ρ_k s_k s_k (31)

ρ_k  = 1 / (y_k T * s_k)   V_k = 1-ρ_k * y_k *s_k T            s_k = x_k+1 – x_k     y_k = ▽f(x_k+1) - ▽f(x_k)            

從上面的 H_k 的疊代計算公式可知，H_k 會慢慢由稀疏矩陣轉變為稠密矩陣，是以存儲該矩陣以及進行該矩陣和向量的相乘運算的消耗将較大。為了避免該問題，LBFGS 算法在 BFGS 算法的基礎上從兩點進行了改進：

1）估算每一步對應的 Hessian 近似矩陣時，給出一個目前步的初始 Hessian 矩陣估計 H_k0

2) 利用過去目前代及過去 m-1 代的 curvature 資訊修正初始 Hessian 矩陣估計 H_k0，得到終的 Hessian 矩陣近似估計 H_k。 計算式如下所示：

H_k = (V_k-1 T … V_k-m T) H_k0(V_k-m … V_k-1)

+ρ_k-m * (V_k-1 T … V_k-m+1 T)* s_k-m*s_k-m T *(V_k-m+1 … V_k-1) 
 +ρ_k-m+1*(V_k-1 T … V_k-m+2 T)*s_k-m+1*s_k-m+1 T *(V_k-m+2 … V_k-1) 
 + … 
 +ρ_k-1*s_k-1*s_k-1 T  (32)            

上述計算式(32），可以通過公式(31)遞歸計算擷取。公式(32)可以用以下算法表示：

q  ▽f(x_k)

for i = k-1, k-2, …, k-m

α_i ρ_i s_i T q;

q  q -α_i *y_i

end

r  H_k0 * q;

for i = k-m, k-m+1, …, k-1

β ρ_i y_i T r

r  r + s_i*(α_i –β)

stop with result H_k*▽f(x_k) = r

從上面計算 H_k 的公式(32)可知，要估算每個點 x_k 處的 Hessian 矩陣近似，需要給出

初始估計 H_k0，H_k0 一般通過以下公式計算：

H_k0 = r_kI r_k =( s_k-1Ty_k-1)(/y_k-1 T *y_k-1)

有了上面的方向計算算法後，LBFGS 算法用于解無限制優化問題，可以表示為如下算法：

選擇一個初始點 x_0，并選擇收斂判斷條件 ε> 0，以及常量 m(代表過去代數)一般為 6

k  0 H_0  I，是以 r = H_0 *▽f(x_0) =▽f(x_0)

p_k = –r

if k > m

删除 LBFGS 計算 H_k 時用不上的向量對(s_k-m, y_k-m)            

計算并儲存 s_k = x_k+1 – x_k y_k = ▽f(x_k+1) - ▽f(x_k)

采用 LBFGS Hessian 矩陣近似算法計算 r

4．算法總結 用于解無限制優化算法的 Quasi-Newton Method 中的 LBFGS 算法到這裡總算初步介紹完了，不過這裡筆者要承認的是這篇文檔省略了許多内容，包括算法收斂性的證明以及收斂速度證明等許多内容。是以讀者若希望對這一塊有一個更深入的認識可以參考以下兩本書：

Numerical Methods for Unconstrained Optimization and Nonlinear Equations

J.E. Dennis Jr. Robert B. Schnabel

Numerical Optimization

Jorge Nocedal Stephen J. Wright

同時我想引用一下侯捷老師的話：

大哉問。學習需要明師。但明師可遇不可求，是以退而求其次你需要好書，并盡早建立自修的基礎。迷時師渡，悟了自渡，尋好書看好書，就是你的自渡法門。切記，徒學不足以自行，計算機是實作性很強的一門科技，你一定要動手做，忌諱眼高手低。學而不思則罔，思而不學則殆，一定要思考、沉澱、整理。 是以這裡附上一個我閱讀後感覺代碼還比較美觀的 LBFGS 實作（C++）：

http://www.chokkan.org/software/liblbfgs/

Naoaki Okazaki

後我不得不承認這篇文檔，對許多讀者來說也許是幾張廢紙，甚至使讀者昏昏欲睡，

或煩躁。因為該文檔從頭至尾并未涉及到一個例子。-_- ||

老天啊，饒恕我吧，畢竟願望是美好的，但願它能對你了解 LBFGS 算法有一點幫助，

哪怕是一點點…

jianzhu                                                                            [email protected] 
                                                        2009-11-1-2009-12-1            

注：這裡以及後面所說的全局收斂，均指的是可以收斂到一個局部小值處 ↑

注：Quasi-Newton Method 中要求 0< c_1< 0.5 ↑

關于矩陣正定可以查閱 Linear Algebra and Its Applications 一書的後一章 ↑