# Alink漫談(十一) ：線性迴歸 之 L-BFGS優化 [TOC] ## 0x00 摘要 Alink 是阿里巴巴基於實時計算引擎 Flink 研發的新一代機器學習演算法平臺，是業界首個同時支援批式演算法、流式演算法的機器學習平臺。本文介紹了線性迴歸的L-BFGS優化在Alink是如何實現的，希望可以作為大家看線性迴歸程式碼的Roadmap。 因為Alink的公開資料太少，所以以下均為自行揣測，肯定會有疏漏錯誤，希望大家指出，我會隨時更新。 本系列目前已有十一篇，歡迎大家指點 [Alink漫談(一) : 從KMeans演算法實現不同看Alink設計思想](https://www.cnblogs.com/rossiXYZ/p/12831175.html) [Alink漫談(二) : 從原始碼看機器學習平臺Alink設計和架構](https://www.cnblogs.com/rossiXYZ/p/12861848.html) [[Alink漫談之三\] AllReduce通訊模型](https://www.cnblogs.com/rossiXYZ/p/12898743.html) [Alink漫談(四) : 模型的來龍去脈](https://www.cnblogs.com/rossiXYZ/p/12940942.html) [Alink漫談(五) : 迭代計算和Superstep](https://www.cnblogs.com/rossiXYZ/p/12990632.html) [Alink漫談(六) : TF-IDF演算法的實現](https://www.cnblogs.com/rossiXYZ/p/13052449.html) [Alink漫談(七) : 如何劃分訓練資料集和測試資料集](https://www.cnblogs.com/rossiXYZ/p/13110960.html) [Alink漫談(八) : 二分類評估 AUC、K-S、PRC、Precision、Recall、LiftChart 如何實現](https://www.cnblogs.com/rossiXYZ/p/13194023.html) [Alink漫談(九) ：特徵工程 之 特徵雜湊/標準化縮放](https://www.cnblogs.com/rossiXYZ/p/13233067.html) [Alink漫談(十) ：線性迴歸實現 之 資料預處理](https://www.cnblogs.com/rossiXYZ/p/13282333.html) ## 0x01 回顧 到目前為止，已經處理完畢輸入，接下來就是優化。優化的主要目標是找到一個方向，引數朝這個方向移動之後使得損失函式的值能夠減小，這個方向往往由一階偏導或者二階偏導各種組合求得。 所以我們再次複習下基本思路。 ### 1.1 優化基本思路 對於線性迴歸模型 f(x) = w'x+e，我們構造一個Cost函式（損失函式）**J(θ)**，並且通過找到 **J(θ)** 函式的最小值，就可以確定線性模型的係數 w 了。 最終的優化函式是：min(L(Y, f(x)) + J(x)) ，即最優化經驗風險和結構風險，而這個函式就被稱為**目標函式**。 我們要做的是依據我們的訓練集，選取最優的θ，在我們的訓練集中讓f(x)儘可能接近真實的值。我們定義了一個函式來描述 “f(x)和真實的值y之間的差距“，這個函式稱為目標函式，表示式如下： $$ J(θ)≈\frac{1}{2} \sum_{i≈1}^m(f_θ(x^{(i)} — y^{(i)})^2\\ $$ 這裡的目標函式就是著名的**最小二乘函式**。 我們要選擇最優的θ，使得f(x)最近進真實值。這個問題就轉化為求解最優的θ，使目標函式 J(θ) 取最小值。 ### 1.2 各類優化方法 尋找合適的*W*令目標函式*f(W)* 最小，是一個無約束最優化問題，解決這個問題的通用做法是隨機給定一個初始的*W0*，通過迭代，在每次迭代中計算目標函式的下降方向並更新*W*，直到目標函式穩定在最小的點。 不同的優化演算法的區別就在於目標函式下降方向*Dt*的計算。下降方向是通過對目標函式在當前的*W*下求一階倒數（梯度，Gradient）和求二階導數（海森矩陣，Hessian Matrix）得到。常見的演算法有梯度下降法、牛頓法、擬牛頓法。 - 梯度下降法直接採用目標函式在當前*W*的梯度的反方向作為下降方向。 - 牛頓法是在當前*W*下，利用二次泰勒展開近似目標函式，然後利用該近似函式來求解目標函式的下降方向。其中*Bt*為目標函式*f(W)*在*Wt*處的海森矩陣。這個搜尋方向也稱作牛頓方向。 - 擬牛頓法只要求每一步迭代中計算目標函式的梯度，通過擬合的方式找到一個近似的海森矩陣用於計算牛頓方向。 - L-BFGS(Limited-memory BFGS)則是解決了BFGS中每次迭代後都需要儲存*N\*N*階海森逆矩陣的問題，只需要儲存每次迭代的兩組向量和一組標量即可。 **Alink中，UnaryLossObjFunc是目標函式，SquareLossFunc 是損失函式，使用L-BFGS演算法對模型進行優化**。 即 優化方向由**擬牛頓法L-BFGS**搞定(具體如何弄就是看f(x)的泰勒二階展開)，損失函式最小值是**平方損失函式**來計算。 ## 0x02 基本概念 因為L-BFGS演算法是擬牛頓法的一種，所以我們先從牛頓法的本質泰勒展開開始介紹。 ### 2.1 泰勒展開 泰勒展開是希望基於某區間一點x_0展開，用一組簡單的冪函式來近似一個複雜的函式f(x)在該區間的區域性。泰勒展開的應用場景例如：我們很難直接求得sin(1)的值，但我們可以通過sin的泰勒級數求得sin(1)的近似值，且展開項越多，精度越高。計算機一般都是把 sin(x)進行泰勒展開進行計算的。 泰勒當年為什麼要發明這條公式？ 因為當時數學界對簡單函式的研究和應用已經趨於成熟，而複雜函式，比如：f(x) = sin(x)ln(1+x) 這種一看就頭疼的函式，還有那種根本就找不到表示式的曲線。除了代入一個x可以得到它的y，就啥事都很難幹了。所以泰勒同學就迎難而上！決定讓這些式子統統現出原形，統統變簡單。 要讓一個複雜函式變簡單，能不能把它轉換成別的表示式？泰勒公式一句話描述：就是用多項式函式去逼近光滑函式。即，根據“以直代曲、化整為零”的數學思想，產生了泰勒公式。 泰勒公式**通過把【任意函式表示式】轉換（重寫）為【多項式】形式**，是一種極其強大的函式**近似工具**。為什麼說它強大呢？ - 多項式非常【友好】，三易，易計算，易求導，易積分 - 幾何感覺和計算感覺都很直觀，如拋物線和幾次方就是底數自己乘自己乘幾次 #### 如何通俗推理？ 泰勒公式乾的事情就是：**使用多項式表示式估計（近似） 在 附近的值**。 當我們想要仿造一個東西的時候，無形之中都會按照如下思路，即先保證大體上相似，再保證區域性相似，再保證細節相似，再保證更細微的地方相似……不斷地細化下去，無窮次細化以後，仿造的東西將無限接近真品。真假難辨。 物理學家得出結論：把生活中關於“仿造”的經驗運用到運動學問題中，如果想仿造一段曲線，那麼首先應該保證曲線的起始點一樣，其次保證起始點處位移隨時間的變化率一樣（速度相同），再次應該保證前兩者相等的同時關於時間的二階變化率一樣（加速度相同）……如果隨時間每一階變化率（每一階導數）都一樣，那這倆曲線肯定是完全等價的。 所以如果泰勒想一個辦法讓自己避免接觸 sin(x)這類函式，即**把這類函式替換掉。** 就可以根據這類函式的影象，仿造一個影象，與原來的影象相類似，這種行為在數學上叫近似。不扯這個名詞。講講如何仿造影象。 仿造的第一步，就是讓仿造的曲線也過這個點。 完成了仿造的第一步，很粗糙，甚至完全看不出來這倆有什麼相似的地方，那就繼續細節化。開始考慮曲線的變化趨勢，即導數，保證在此處的導數相等。 經歷了第二步，現在起始點相同了，整體變化趨勢相近了，可能看起來有那麼點意思了。想進一步精確化，應該考慮凹凸性。高中學過：表徵影象的凹凸性的引數為“導數的導數”。所以，下一步就讓二者的導數的導數相等。 起始點相同，增減性相同，凹凸性相同後，仿造的函式更像了。如果再繼續細化下去，應該會無限接近。所以泰勒認為“**仿造一段曲線，要先保證起點相同，再保證在此處導數相同，繼續保證在此處的導數的導數相同……**” 泰勒展開式就是把一個三角函式或者指數函式或者其他比較難纏的函式用多項式替換掉。 也就是說，有一個**原函式f(x)** ，我再造一個影象與原函式影象相似的**多項式函式 g(x)** ，為了保證相似，我只需要保證這倆函式在某一點的**初始值相等，1階導數相等，2階導數相等，……n階導數相等**。 ### 2.2 牛頓法 牛頓法的基本思路是，**在現有極小點估計值的附近對 f(x) 做二階泰勒展開，進而找到極小點的下一個估計值**。其核心思想是利用迭代點 x_k 處的一階導數(梯度)和二階導數(Hessen矩陣)對目標函式進行二次函式近似，然後把二次模型的極小點作為新的迭代點，並不斷重複這一過程，直至求得滿足精度的近似極小值。 梯度下降演算法是將函式在 xn 位置進行一次函式近似，也就是一條直線。計算梯度，從而決定下一步優化的方向是梯度的反方向。而牛頓法是將函式在 xn 位置進行二階函式近似，也就是二次曲線。計算梯度和二階導數，從而決定下一步的優化方向。 我們要優化的都是多元函式，x往往不是一個實數，而是一個向量。所以f(x) 的一階導數也是一個向量，再求導的二階導數是一個矩陣，就是Hessian矩陣。 #### 2.2.1 泰勒一階展開 牛頓法求根可以按照泰勒一階展開。例如對於方程 f(x) = 0，我們在任意一點 x0 處，進行一階泰勒展開： $$ f(x) = f(x_0) + (x - x_0)f^{ '}(x_0) $$ 令 f(x) = 0，帶入上式，即可得到： $$ x = x_0 - \frac{f(x_0)}{f^{'}(x_0)} $$ 注意，這裡使用了近似，得到的 x 並不是方程的根，只是近似解。但是，x 比 x0 更接近於方程的根。 然後，利用迭代方法求解，以 x 為 x0，求解下一個距離方程的根更近的位置。迭代公式可以寫成： $$ x_{n+1} = x_n - \frac{f(x_n)}{f^{'}(x_n)} $$ #### 2.2.2 泰勒二階展開 機器學習、深度學習中，損失函式的優化問題一般是基於一階導數梯度下降的。現在，從另一個角度來看，想要讓損失函式最小化，這其實是一個最值問題，對應函式的一階導數 f'(x) = 0。也就是說，如果我們找到了能讓 f'(x) = 0 的點 x，損失函式取得最小值，也就實現了模型優化目標。 現在的目標是計算 f’(x) = 0 對應的 x，即 f’(x) = 0 的根。轉化為求根問題，就可以利用上一節的牛頓法了。 與上一節有所不同，首先，對 f(x) 在 x0 處進行二階泰勒展開： $$ f(x) = f(x_0) + (x - x_0)f^{'}(x_0) + \frac{1}{2}(x-x_0)^2f^{''}(x_0) $$ 上式成立的條件是 f(x) 近似等於 f(x0)。令 f(x) = f(x0)，並對 (x - x0) 求導，可得： $$ x = x_0 - \frac{f^{'}(x_0)}{f^{''}(x_0)} $$ 同樣，雖然 x 並不是最優解點，但是 x 比 x0 更接近 f'(x) = 0 的根。這樣，就能得到最優化的迭代公式： $$ x_{n+1} = x_n - \frac{f^{'}(x_n)}{f^{''}(x_n)} $$ 通過迭代公式，就能不斷地找到 f'(x) = 0 的近似根，從而也就實現了損失函式最小化的優化目標。 #### 2.2.3 高維空間 在機器學習的優化問題中，我們要優化的都是多元函式，所以x往往不是一個實數，而是一個向量。所以將牛頓求根法利用到機器學習中時,x是一個向量，f(x) 的一階導數也是一個向量，再求導的二階導數是一個矩陣，就是Hessian矩陣。 在高維空間，我們用梯度替代一階導數，用Hessian矩陣替代二階導數，牛頓法的迭代公式不變： $$ x_{k+1} = x_k - [Hf(x_k)^{-1}].J_f(x_k) $$ 其中 J 定義為 **雅克比矩陣，對應一階偏導數。** 推導如下 ： 我們假設f(x)是二階可微實函式，把f(x)在xk處Taylor展開並取二階近似為 $$ \begin{aligned}&f(x)≈f(x^k)+∇f(x^k)T(x−x^k)+\frac{1}{2}(x−x^k)^T∇^2f(x^k)(x−x^k)\\&x^k為當前的極小值估計值\\&∇f(x^k)是f(x)在x^k處的一階導數\\&∇^2f(x) 是f(x)在x^k處的Hessen矩陣。\\\end{aligned} $$ 我們的目標是求f(x)的最小值，而導數為0的點極有可能為極值點，故在此對f(x)求導，並令其導數為0，即∇f(x)=0，可得 $$ ∇f(x)=∇f(x^k)+∇^2f(x^k)(x−x^k)=0 $$ 設∇2f(x)可逆，由(2)可以得到牛頓法的迭代公式 $$ \begin{aligned}&x^{k+1}=x^k−∇^2f(x^k)^{−1}∇f(x^k) \\&d= −∇^2f(x^k)^{−1}∇f(x^k)被稱為牛頓方向 \\\end{aligned} $$ 當原函式存在一階二階連續可導時，可以採用牛頓法進行一維搜尋，收斂速度快，具有區域性二階收斂速度。 #### 2.2.4 牛頓法基本流程 總結(模仿)一下使用牛頓法求根的步驟: ​ a.已知函式的情況下,隨機產生點. ​ b.由已知的點按照的公式進行k次迭代. ​ c.如果與上一次的迭代結果相同或者相差結果小於一個閾值時,本次結果就是函式的根. 虛擬碼如下 ```python def newton(feature, label, iterMax, sigma, delta): '''牛頓法 input: feature(mat):特徵 label(mat):標籤 iterMax(int):最大迭代次數 sigma(float), delta(float):牛頓法中的引數 output: w(mat):迴歸係數 ''' n = np.shape(feature)[1] w = np.mat(np.zeros((n, 1))) it = 0 while it <= iterMax: g = first_derivativ(feature, label, w) # 一階導數 G = second_derivative(feature) # 二階導數 d = -G.I * g m = get_min_m(feature, label, sigma, delta, d, w, g) # 得到步長中最小的值m w = w + pow(sigma, m) * d it += 1 return w ``` #### 2.2.5 問題點及解決 牛頓法不僅需要計算 Hessian 矩陣，而且需要計算 Hessian 矩陣的逆。當資料量比較少的時候，運算速度不會受到大的影響。但是，當資料量很大，特別在深度神經網路中，計算 Hessian 矩陣和它的逆矩陣是非常耗時的。從整體效果來看，牛頓法優化速度沒有梯度下降演算法那麼快。所以，目前神經網路損失函式的優化策略大多都是基於梯度下降。 另一個問題是，如果某個點的Hessian矩陣接近奇異（條件數過大），逆矩陣會導致數值不穩定，甚至迭代可能不會收斂。 當x的維度特別多的時候，我們想求得f(x) 的二階導數很困難，而牛頓法求駐點又是一個迭代演算法，所以這個困難我們還要面臨無限多次，導致了牛頓法求駐點在機器學習中無法使用。有沒有什麼解決辦法呢？ **實際應用中，我們通常不去求解逆矩陣，而是考慮求解Hessian矩陣的近似矩陣，最好是隻利用一階導數近似地得到二階導數的資訊，從而在較少的計算量下得到接近二階的收斂速率。這就是擬牛頓法。** 擬牛頓法的想法其實很簡單，就像是函式值的兩點之差可以逼近導數一樣，一階導數的兩點之差也可以逼近二階導數。幾何意義是求一階導數的“割線”，取極限時，割線會逼近切線，矩陣B就會逼近Hessian矩陣。 $$ 擬牛頓法，同樣可以顧名思義，就是模擬牛頓法，用一個近似於∇^2f(x)^{−1}的矩陣H_{k+1}來替代∇^2f(x)^{−1} $$ ### 2.3 擬牛頓法 擬牛頓法就是模擬出Hessen矩陣的構造過程，通過迭代來逼近。主要包括DFP擬牛頓法，BFGS擬牛頓法。擬牛頓法只要求每一步迭代時知道目標函式的梯度。 **各種擬牛頓法使用迭代法分別近似海森矩陣的逆和它自身；** 在各種擬牛頓法中，一般的構造Hk+1的策略是，H1通常選擇任意的一個n階對稱正定矩陣(一般為I)，然後通過不斷的修正Hk給出Hk+1，即 $$ H_{k+1}=H_k+ΔH_k \\ ΔH_k稱為校正矩陣 $$ 比如：BFGS法每次更新矩陣H(Hessian矩陣的逆矩陣)需要的是第k步的迭代點差s和梯度差y，第k+1步的H相當於需要從開始到第k步的所用s和y的值。 我們要通過牛頓求駐點法和BFGS演算法來求得一個函式的根，兩個演算法都需要迭代，所以我們乾脆讓他倆一起迭代就好了。兩個演算法都是慢慢逼近函式根，所以經過k次迭代以後，所得到的解就是機器學習中目標函式導函式的根。這種兩個演算法共同迭代的計算方式，我們稱之為On The Fly. 在BFGS演算法迭代的第一步，單位矩陣與梯度 g 相乘，就等於梯度 g，形式上同梯度下降的表示式是相同的。所以BFGS演算法可以理解為從梯度下降逐步轉換為牛頓法求函式解的一個演算法。 雖然我們使用了BFGS演算法來利用單位矩陣逐步逼近H矩陣，但是每次計算的時候都要儲存D矩陣，D矩陣有多大。呢。假設我們的資料集有十萬個維度(不算特別大)，那麼每次迭代所要儲存D矩陣的結果是74.5GB。我們無法儲存如此巨大的矩陣內容，如何解決呢？使用L-BFGS演算法. ### 2.4 L-BFGS演算法 L-BFGS演算法的基本思想是：演算法只儲存並利用最近m次迭代的曲率資訊來構造海森矩陣的近似矩陣。 我們要通過牛頓求駐點法和BFGS演算法來求得一個函式的根，兩個演算法都需要迭代，所以我們乾脆讓他倆一起迭代就好了。兩個演算法都是慢慢逼近函式根，所以經過k次迭代以後，所得到的解就是機器學習中目標函式導函式的根。這種兩個演算法共同迭代的計算方式，我們稱之為On The Fly. 在BFGS演算法迭代的第一步，單位矩陣與梯度g相乘，就等於梯度g，形式上同梯度下降的表示式是相同的。所以BFGS演算法可以理解為從梯度下降逐步轉換為牛頓法求函式解的一個演算法。 > 雖然我們使用了BFGS演算法來利用單位矩陣逐步逼近H矩陣，但是每次計算的時候都要儲存D矩陣，D矩陣有多大。呢。假設我們的資料集有十萬個維度(不算特別大)，那麼每次迭代所要儲存D矩陣的結果是74.5GB。我們無法儲存如此巨大的矩陣內容，如何解決呢？使用L-BFGS演算法. 我們每一次對D矩陣的迭代，都是通過迭代計算sk和yk得到的。我們的記憶體存不下時候只能丟掉一些存不下的資料。假設我們設定的儲存向量數為100，當s和y迭代超過100時，就會扔掉第一個s和y。每多一次迭代就對應的扔掉最前邊的s和y。這樣雖然損失了精度，但確可以保證使用有限的記憶體將函式的解通過BFGS演算法求得到。 所以L-BFGS演算法可以理解為對BFGS演算法的又一次近似或者逼近。 這裡不介紹數學論證，因為網上優秀文章有很多，這裡只是介紹工程實現。總結L-BFGS演算法的大致步驟如下： Step1: 選初始點x_0，儲存最近迭代次數m； Step2: k=0, H_0=I, r=f(x_0)； Step3: 根據更新的引數計算梯度和損失值，如果達到閾值，則返回最優解x_{k+1}，否則轉Step4； Step4: 計算本次迭代的可行梯度下降方向 p_k=-r_k； Step5: 計算步長alpha_k，進行一維搜尋; Step6: 更新權重x; Step7: 只保留最近m次的向量對； Step8: 計算並儲存 s_k, y_k Step9: 用two-loop recursion演算法求r_k； k=k+1，轉Step3。 ## 0x03 優化模型 -- L-BFGS演算法 回到程式碼，`BaseLinearModelTrainBatchOp.optimize`函式呼叫的是 ```java return new Lbfgs(objFunc, trainData, coefficientDim, params).optimize(); ``` 優化後將返回線性模型的係數。 ```java /** * optimizer api. * * @return the coefficient of linear problem. */ @Override public