Alink漫談(十一) ：線性迴歸 之 L-BFGS優化
0x00 摘要
Alink 是阿里巴巴基於實時計算引擎 Flink 研發的新一代機器學習演算法平臺，是業界首個同時支援批式演算法、流式演算法的機器學習平臺。本文介紹了線性迴歸的L-BFGS優化在Alink是如何實現的，希望可以作為大家看線性迴歸程式碼的Roadmap。
因為Alink的公開資料太少，所以以下均為自行揣測，肯定會有疏漏錯誤，希望大家指出，我會隨時更新。
本系列目前已有十一篇，歡迎大家指點
Alink漫談(一) : 從KMeans演算法實現不同看Alink設計思想
Alink漫談(二) : 從原始碼看機器學習平臺Alink設計和架構
Alink漫談(八) : 二分類評估 AUC、K-S、PRC、Precision、Recall、LiftChart 如何實現
0x01 回顧
到目前為止，已經處理完畢輸入，接下來就是優化。優化的主要目標是找到一個方向，引數朝這個方向移動之後使得損失函式的值能夠減小，這個方向往往由一階偏導或者二階偏導各種組合求得。 所以我們再次複習下基本思路。
1.1 優化基本思路
對於線性迴歸模型 f(x) = w'x+e，我們構造一個Cost函式（損失函式）J(θ)，並且通過找到 J(θ) 函式的最小值，就可以確定線性模型的係數 w 了。
最終的優化函式是：min(L(Y, f(x)) + J(x)) ，即最優化經驗風險和結構風險，而這個函式就被稱為目標函式。
我們要做的是依據我們的訓練集，選取最優的θ，在我們的訓練集中讓f(x)儘可能接近真實的值。我們定義了一個函式來描述 “f(x)和真實的值y之間的差距“，這個函式稱為目標函式，表示式如下：
這裡的目標函式就是著名的最小二乘函式。
我們要選擇最優的θ，使得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) = 0，帶入上式，即可得到：
注意，這裡使用了近似，得到的 x 並不是方程的根，只是近似解。但是，x 比 x0 更接近於方程的根。
然後，利用迭代方法求解，以 x 為 x0，求解下一個距離方程的根更近的位置。迭代公式可以寫成：
2.2.2 泰勒二階展開
機器學習、深度學習中，損失函式的優化問題一般是基於一階導數梯度下降的。現在，從另一個角度來看，想要讓損失函式最小化，這其實是一個最值問題，對應函式的一階導數 f'(x) = 0。也就是說，如果我們找到了能讓 f'(x) = 0 的點 x，損失函式取得最小值，也就實現了模型優化目標。
現在的目標是計算 f’(x) = 0 對應的 x，即 f’(x) = 0 的根。轉化為求根問題，就可以利用上一節的牛頓法了。
與上一節有所不同，首先，對 f(x) 在 x0 處進行二階泰勒展開：
上式成立的條件是 f(x) 近似等於 f(x0)。令 f(x) = f(x0)，並對 (x - x0) 求導，可得：
同樣，雖然 x 並不是最優解點，但是 x 比 x0 更接近 f'(x) = 0 的根。這樣，就能得到最優化的迭代公式：
通過迭代公式，就能不斷地找到 f'(x) = 0 的近似根，從而也就實現了損失函式最小化的優化目標。
2.2.3 高維空間
在機器學習的優化問題中，我們要優化的都是多元函式，所以x往往不是一個實數，而是一個向量。所以將牛頓求根法利用到機器學習中時,x是一個向量，f(x) 的一階導數也是一個向量，再求導的二階導數是一個矩陣，就是Hessian矩陣。
在高維空間，我們用梯度替代一階導數，用Hessian矩陣替代二階導數，牛頓法的迭代公式不變：
其中 J 定義為 雅克比矩陣，對應一階偏導數。
推導如下 ：
我們假設f(x)是二階可微實函式，把f(x)在xk處Taylor展開並取二階近似為
我們的目標是求f(x)的最小值，而導數為0的點極有可能為極值點，故在此對f(x)求導，並令其導數為0，即∇f(x)=0，可得
設∇2f(x)可逆，由(2)可以得到牛頓法的迭代公式
當原函式存在一階二階連續可導時，可以採用牛頓法進行一維搜尋，收斂速度快，具有區域性二階收斂速度。
2.2.4 牛頓法基本流程
總結(模仿)一下使用牛頓法求根的步驟:
a.已知函式的情況下,隨機產生點.
b.由已知的點按照的公式進行k次迭代.
c.如果與上一次的迭代結果相同或者相差結果小於一個閾值時,本次結果就是函式的根.
虛擬碼如下
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矩陣。
2.3 擬牛頓法
擬牛頓法就是模擬出Hessen矩陣的構造過程，通過迭代來逼近。主要包括DFP擬牛頓法，BFGS擬牛頓法。擬牛頓法只要求每一步迭代時知道目標函式的梯度。
各種擬牛頓法使用迭代法分別近似海森矩陣的逆和它自身；
在各種擬牛頓法中，一般的構造Hk+1的策略是，H1通常選擇任意的一個n階對稱正定矩陣(一般為I)，然後通過不斷的修正Hk給出Hk+1，即
比如：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函式呼叫的是
return new Lbfgs(objFunc, trainData, coefficientDim, params).optimize();
優化後將返回線性模型的係數。
/**
* optimizer api.
*
* @return the coefficient of linear problem.
*/
@Override
public DataSet <Tuple2 <DenseVector, double[]>> optimize() {
/**
* solving problem using iteration.
* trainData is the distributed samples.
* initCoef is the initial model coefficient, which will be broadcast to every worker.
* objFuncSet is the object function in dataSet format
* .add(new PreallocateCoefficient(OptimName.currentCoef)) allocate memory for current coefficient
* .add(new PreallocateCoefficient(OptimName.minCoef)) allocate memory for min loss coefficient
* .add(new PreallocateLossCurve(OptimVariable.lossCurve)) allocate memory for loss values
* .add(new PreallocateVector(OptimName.dir ...)) allocate memory for dir
* .add(new PreallocateVector(OptimName.grad)) allocate memory for grad
* .add(new PreallocateSkyk()) allocate memory for sK yK
* .add(new CalcGradient(objFunc)) calculate local sub gradient
* .add(new AllReduce(OptimName.gradAllReduce)) sum all sub gradient with allReduce
* .add(new CalDirection()) get summed gradient and use it to calc descend dir
* .add(new CalcLosses(objFunc, OptimMethod.GD)) calculate local losses for line search
* .add(new AllReduce(OptimName.lossAllReduce)) sum all losses with allReduce
* .add(new UpdateModel(maxIter, epsilon ...)) update coefficient
* .setCompareCriterionOfNode0(new IterTermination()) judge stop of iteration
*/
DataSet <Row> model = new IterativeComQueue()
.initWithPartitionedData(OptimVariable.trainData, trainData)
.initWithBroadcastData(OptimVariable.model, coefVec)
.initWithBroadcastData(OptimVariable.objFunc, objFuncSet)
.add(new PreallocateCoefficient(OptimVariable.currentCoef))
.add(new PreallocateCoefficient(OptimVariable.minCoef))
.add(new PreallocateLossCurve(OptimVariable.lossCurve, maxIter))
.add(new PreallocateVector(OptimVariable.dir, new double[] {0.0, OptimVariable.learningRate}))
.add(new PreallocateVector(OptimVariable.grad))
.add(new PreallocateSkyk(OptimVariable.numCorrections))
.add(new CalcGradient())
.add(new AllReduce(OptimVariable.gradAllReduce))
.add(new CalDirection(OptimVariable.numCorrections))
.add(new CalcLosses(OptimMethod.LBFGS, numSearchStep))
.add(new AllReduce(OptimVariable.lossAllReduce))
.add(new UpdateModel(params, OptimVariable.grad, OptimMethod.LBFGS, numSearchStep))
.setCompareCriterionOfNode0(new IterTermination())
.closeWith(new OutputModel())
.setMaxIter(maxIter)
.exec();
return model.mapPartition(new ParseRowModel());
}
所以我們接下來的就是看Lbfgs，其重點就是引數更新的下降方向和搜尋步長。
3.1 如何分散式實施
如果由於所有輸入的資料都是相同維度的；演算法過程中不會對輸入修改，就可以將這些輸入資料進行切分。這樣的話，應該可以通過一次map-reduce來計算。
我們理一下L-BFGS中可以分散式平行計算的步驟：
- 計算梯度 可以並行化 ，比如在機器1計算梯度1，在機器2計算梯度2....，然後通過一個Reduce把這些合成一個完整的梯度向量。
- 計算方向 可以並行化，同樣可以通過把資料分割槽，然後Map算各自分割槽上的值，Reduce合起來得到方向。
- 計算損失 可以並行化，同樣可以通過把資料分割槽，然後Map算各自分割槽上的值，Reduce合起來得到損失。
Alink中，使用AllReduce功能而非Flink原生Map / Reduce來完成了以上三點的平行計算和通訊。
3.2 CalcGradient
線搜尋只有兩步，確定方向、確定步長。確定方向和模擬Hessian矩陣都需要計算梯度。目標函式的梯度向量計算中只需要進行向量間的點乘和相加，可以很容易將每個迭代過程拆分成相互獨立的計算步驟，由不同的節點進行獨立計算，然後歸併計算結果。
Alink將樣本特徵向量分佈到不同的計算節點，由各計算節點完成自己所負責樣本的點乘與求和計算，然後將計算結果進行歸併，則實現了按行並行的LR。
實際情況中也會存在針對高維特徵向量進行邏輯迴歸的場景，僅僅按行進行並行處理，無法滿足這類場景的需求，因此還需要按列將高維的特徵向量拆分成若干小的向量進行求解。這個也許是Alink以後需要優化的一個點吧 。
CalcGradient是Alink迭代運算元的派生類，函式總結如下：
- 獲取經過處理的輸入資料labledVectors，
- 從靜態記憶體中獲取 "迭代引數coef"，"優化函式objFunc" 和 "梯度"；
- 計算區域性梯度 objFunc.calcGradient(labledVectors, coef, grad.f0); 這裡呼叫到了目標函式的梯度相關API；
objFunc.calcGradient根據取樣點計算梯度。Alink這裡就是把x, y代入，求損失函式的梯度就是對 Coef求偏導數。具體我們在前文已經提到。- 對計算樣本中的每一個樣本，分別計算不同特徵的計算梯度。
- 將新計算出來的梯度乘以權重之後，存入靜態記憶體 gradAllReduce 中。
- 後續會通過聚合函式，對所有計算樣本的特徵的梯度進行累加，得到每一個特徵的累積梯度以及損失。
public class CalcGradient extends ComputeFunction {
/**
* object function class, it supply the functions to calc local gradient (or loss).
*/
private OptimObjFunc objFunc;
@Override
public void calc(ComContext context) {
Iterable<Tuple3<Double, Double, Vector>> labledVectors = context.getObj(OptimVariable.trainData);
// 經過處理的輸入資料
labledVectors = {ArrayList@9877} size = 4
0 = {Tuple3@9895} "(1.0,16.8,1.0 1.0 1.4657097546055162 1.4770978917519928)"
1 = {Tuple3@9896} "(1.0,6.7,1.0 1.0 -0.338240712601273 -0.7385489458759964)"
2 = {Tuple3@9897} "(1.0,6.9,1.0 1.0 -0.7892283294029703 -0.3692744729379982)"
3 = {Tuple3@9898} "(1.0,8.0,1.0 1.0 -0.338240712601273 -0.3692744729379982)"
// get iterative coefficient from static memory.
Tuple2<DenseVector, Double> state = context.getObj(OptimVariable.currentCoef);
int size = state.f0.size();
// 是Coef，1.7976931348623157E308是預設最大值
state = {Tuple2@9878} "(0.001 0.0 0.0 0.0,1.7976931348623157E308)"
f0 = {DenseVector@9879} "0.001 0.0 0.0 0.0"
f1 = {Double@9889} 1.7976931348623157E308
DenseVector coef = state.f0;
if (objFunc == null) {
objFunc = ((List<OptimObjFunc>)context.getObj(OptimVariable.objFunc)).get(0);
}
// 變數如下
objFunc = {UnaryLossObjFunc@9882}
unaryLossFunc = {SquareLossFunc@9891}
l1 = 0.0
l2 = 0.0
params = {Params@9892} "Params {featureCols=["f0","f1","f2"], labelCol="label", predictionCol="linpred"}"
Tuple2<DenseVector, double[]> grad = context.getObj(OptimVariable.dir);
// 變數如下
grad = {Tuple2@9952} "(0.0 0.0 0.0 0.0,[0.0, 0.1])"
f0 = {DenseVector@9953} "0.0 0.0 0.0 0.0"
f1 = {double[2]@9969}
coef = {DenseVector@9951} "0.001 0.0 0.0 0.0"
data = {double[4]@9982}
// calculate local gradient，使用目標函式
Double weightSum = objFunc.calcGradient(labledVectors, coef, grad.f0);
// prepare buffer vec for allReduce. the last element of vec is the weight Sum.
double[] buffer = context.getObj(OptimVariable.gradAllReduce);
if (buffer == null) {
buffer = new double[size + 1];
context.putObj(OptimVariable.gradAllReduce, buffer);
}
for (int i = 0; i < size; ++i) {
buffer[i] = grad.f0.get(i) * weightSum;
}
/* the last element is the weight value */
buffer[size] = weightSum;
}
}
// 最後結果是
buffer = {double[5]@9910}
0 = -38.396
1 = -38.396
2 = -14.206109929253465
3 = -14.364776997288134
4 = 0.0
3.3 AllReduce
這裡是前面提到的 "通過聚合函式，對所有計算樣本的特徵的梯度進行累加，得到每一個特徵的累積梯度以及損失"。
具體關於AllReduce如何運作，可以參見文章 [Alink漫談之三] AllReduce通訊模型
.add(new AllReduce(OptimVariable.gradAllReduce))
3.4 CalDirection
此時得到的梯度，已經是聚合之後的，所以可以開始計算方向。
3.4.1 預先分配
OptimVariable.grad 是預先分配的。
public class PreallocateSkyk extends ComputeFunction {
private int numCorrections;
/**
* prepare hessian matrix of lbfgs method. we allocate memory fo sK, yK at first iteration step.
*
* @param context context of iteration.
*/
@Override
public void calc(ComContext context) {
if (context.getStepNo() == 1) {
Tuple2<DenseVector, double[]> grad = context.getObj(OptimVariable.grad);
int size = grad.f0.size();
DenseVector[] sK = new DenseVector[numCorrections];
DenseVector[] yK = new DenseVector[numCorrections];
for (int i = 0; i < numCorrections; ++i) {
sK[i] = new DenseVector(size);
yK[i] = new DenseVector(size);
}
context.putObj(OptimVariable.sKyK, Tuple2.of(sK, yK));
}
}
}
3.4.2 計算方向
在計算的過程中，需要不斷的計算和儲存歷史的Hessian矩陣，在L-BFGS演算法，希望只保留最近的m次迭代資訊，便能夠擬合Hessian矩陣。在L-BFGS演算法中，不再儲存完整的Hk，而是儲存向量序列{sk}和{yk}，需要矩陣時Hk，使用向量序列{sk}和{yk}計算就可以得到，而向量序列{sk}和{yk}也不是所有都要儲存，只要儲存最新的m步向量即可。
具體原理和公式這裡不再贅述，網上很多文章講解非常好。
重點說明，dir的各個資料用途是：
dir = {Tuple2@9931} "(-9.599 -9.599 -3.5515274823133662 -3.5911942493220335,[4.0, 0.1])"
f0 = {DenseVector@9954} "-9.599 -9.599 -3.5515274823133662 -3.5911942493220335" //梯度
f1 = {double[2]@9938}
0 = 4.0 //權重
1 = 0.1 //學習率 learning rate，0.1是初始化數值，後續UpdateModel時候會更新
程式碼摘要如下：
@Override
public void calc(ComContext context) {
Tuple2 <DenseVector, double[]> grad = context.getObj(OptimVariable.grad);
Tuple2 <DenseVector, double[]> dir = context.getObj(OptimVariable.dir);
Tuple2 <DenseVector[], DenseVector[]> hessian = context.getObj(OptimVariable.sKyK);
int size = grad.f0.size();
// gradarr是上一階段CalcGradient的結果
double[] gradarr = context.getObj(OptimVariable.gradAllReduce);
// 變數為
gradarr = {double[5]@9962}
0 = -38.396
1 = -38.396
2 = -14.206109929253465
3 = -14.364776997288134
4 = 4.0
if (this.oldGradient == null) {
oldGradient = new DenseVector(size);
}
// hessian用來當作佇列，儲存sK，yK，只保留最近m個
DenseVector[] sK = hessian.f0;
DenseVector[] yK = hessian.f1;
for (int i = 0; i < size; ++i) {
//gradarr[size]是權重
grad.f0.set(i, gradarr[i] / gradarr[size]); //size = 4
}
// 賦值梯度，這裡都除以權重
grad = {Tuple2@9930} "(-9.599 -9.599 -3.5515274823133662 -3.5911942493220335,[0.0])"
f0 = {DenseVector@9937} "-9.599 -9.599 -3.5515274823133662 -3.5911942493220335"
data = {double[4]@9963}
0 = -9.599
1 = -9.599
2 = -3.5515274823133662
3 = -3.5911942493220335
f1 = {double[1]@9961}
0 = 0.0
dir.f1[0] = gradarr[size]; //權重
int k = context.getStepNo() - 1;
if (k == 0) { //首次迭代
dir.f0.setEqual(grad.f0); // 梯度賦予在這裡
oldGradient.setEqual(grad.f0);
} else {
yK[(k - 1) % m].setEqual(grad.f0);
yK[(k - 1) % m].minusEqual(oldGradient);
oldGradient.setEqual(grad.f0);
}
// copy g_k and store in qL
dir.f0.setEqual(grad.f0); //拷貝梯度到這裡
//
dir = {Tuple2@9931} "(-9.599 -9.599 -3.5515274823133662 -3.5911942493220335,[4.0, 0.1])"
f0 = {DenseVector@9954} "-9.599 -9.599 -3.5515274823133662 -3.5911942493220335" //梯度
f1 = {double[2]@9938}
0 = 4.0 //權重
1 = 0.1 //學習率 learning rate，0.1是初始化數值
// compute H^-1 * g_k
int delta = k > m ? k - m : 0;
int l = k <= m ? k : m; // m = 10
if (alpha == null) {
alpha = new double[m];
}
// two-loop的過程，通過擬牛頓法計算Hessian矩陣
for (int i = l - 1; i >= 0; i--) {
int j = (i + delta) % m;
double dot = sK[j].dot(yK[j]);
if (Math.abs(dot) > 0.0) {
double rhoJ = 1.0 / dot;
alpha[i] = rhoJ * (sK[j].dot(dir.f0)); // 計算alpha
dir.f0.plusScaleEqual(yK[j], -alpha[i]); // 重新修正d
}
}
for (int i = 0; i < l; i++) {
int j = (i + delta) % m;
double dot = sK[j].dot(yK[j]);
if (Math.abs(dot) > 0.0) {
double rhoJ = 1.0 / dot;
double betaI = rhoJ * (yK[j].dot(dir.f0)); // 乘以rho
dir.f0.plusScaleEqual(sK[j], (alpha[i] - betaI));// 重新修正d
}
}
}
//最後是儲存在 OptimVariable.dir
3.5 CalcLosses
根據更新的 dir 和 當前係數 計算損失函式誤差值，這個損失是為後續的線性搜尋準備的。目的是如果損失函式誤差值達到允許的範圍，那麼停止迭代，否則重複迭代。
CalcLosses基本邏輯如下：
- 1)得到本次步長 Double beta = dir.f1[1] / numSearchStep;
- 後續UpdateModel 中會對下一次計算的步長（learning rate）進行更新，比如 dir.f1[1] *= 1.0 / (numSearchStep * numSearchStep); 或者 dir.f1[1] = Math.min(dir.f1[1], numSearchStep);
- 2)呼叫目標函式的 calcSearchValues 來計算當前係數對應的損失；
- 3)calcSearchValues 遍歷輸入labelVectors，對於每個 labelVector 按照線性搜尋的步驟進行計算損失。vec[i] += weight * this.unaryLossFunc.loss(etaCoef - i * etaDelta, labelVector.f1); 迴圈內部如下：
- 3.1)用x-vec和coef計算出來的 Y ，etaCoef = getEta(labelVector, coefVector);
- 3.2)以x-vec和dirVec計算出來的 deltaY，etaDelta = getEta(labelVector, dirVec) * beta;
- 3.3)按照線性搜尋的步驟進行計算損失。vec[i] += weight * this.unaryLossFunc.loss(etaCoef - i * etaDelta, labelVector.f1); 聯絡損失函式可知，etaCoef - i * etaDelta, labelVector.f1 是 訓練資料預測值 與 實際類別 的偏差；
- 4)為後續聚合 lossAllReduce 準備資料；
程式碼如下：
public class CalcLosses extends ComputeFunction {
@Override
public void calc(ComContext context) {
Iterable<Tuple3<Double, Double, Vector>> labledVectors = context.getObj(OptimVariable.trainData);
Tuple2<DenseVector, double[]> dir = context.getObj(OptimVariable.dir);
Tuple2<DenseVector, Double> coef = context.getObj(OptimVariable.currentCoef);
if (objFunc == null) {
objFunc = ((List<OptimObjFunc>)context.getObj(OptimVariable.objFunc)).get(0);
}
/**
* calculate losses of current coefficient.
* if optimizer is owlqn, constraint search will used, else without constraint.
*/
Double beta = dir.f1[1] / numSearchStep;
double[] vec = method.equals(OptimMethod.OWLQN) ?
objFunc.constraintCalcSearchValues(labledVectors, coef.f0, dir.f0, beta, numSearchStep)
: objFunc.calcSearchValues(labledVectors, coef.f0, dir.f0, beta, numSearchStep);
// 變數為
dir = {Tuple2@9988} "(-9.599 -9.599 -3.5515274823133662 -3.5911942493220335,[4.0, 0.1])"
coef = {Tuple2@9989} "(0.001 0.0 0.0 0.0,1.7976931348623157E308)"
beta = {Double@10014} 0.025
vec = {double[5]@10015}
0 = 0.0
1 = 0.0
2 = 0.0
3 = 0.0
4 = 0.0
// prepare buffer vec for allReduce.
double[] buffer = context.getObj(OptimVariable.lossAllReduce);
if (buffer == null) {
buffer = vec.clone();
context.putObj(OptimVariable.lossAllReduce, buffer);
} else {
System.arraycopy(vec, 0, buffer, 0, vec.length);
}
}
}
其中搜尋是一元目標函式提供的，其又呼叫了損失函式。
public class UnaryLossObjFunc extends OptimObjFunc {
/**
* Calculate loss values for line search in optimization.
*
* @param labelVectors train data.
* @param coefVector coefficient of current time.
* @param dirVec descend direction of optimization problem.
* @param beta step length of line search.
* @param numStep num of line search step.
* @return double[] losses.
*/
@Override
public double[] calcSearchValues(Iterable<Tuple3<Double, Double, Vector>> labelVectors,
DenseVector coefVector,
DenseVector dirVec,
double beta,
int numStep) {
double[] vec = new double[numStep + 1];
// labelVector是三元組Tuple3<weight, label, feature vector>
labelVectors = {ArrayList@10007} size = 4
0 = {Tuple3@10027} "(1.0,16.8,1.0 1.0 1.4657097546055162 1.4770978917519928)"
1 = {Tuple3@10034} "(1.0,6.7,1.0 1.0 -0.338240712601273 -0.7385489458759964)"
2 = {Tuple3@10035} "(1.0,6.9,1.0 1.0 -0.7892283294029703 -0.3692744729379982)"
3 = {Tuple3@10036} "(1.0,8.0,1.0 1.0 -0.338240712601273 -0.3692744729379982)"
coefVector = {DenseVector@10008} "0.001 0.0 0.0 0.0"
dirVec = {DenseVector@10009} "-9.599 -9.599 -3.5515274823133662 -3.5911942493220335"
beta = 0.025
numStep = 4
vec = {double[5]@10026}
0 = 0.0
1 = 0.0
2 = 0.0
3 = 0.0
4 = 0.0
for (Tuple3<Double, Double, Vector> labelVector : labelVectors) {
double weight = labelVector.f0;
//用x-vec和coef計算出來的 Y
double etaCoef = getEta(labelVector, coefVector);
//以x-vec和dirVec計算出來的 deltaY
double etaDelta = getEta(labelVector, dirVec) * beta;
weight = 1.0
etaCoef = 0.001
etaDelta = -0.7427013482280431
for (int i = 0; i < numStep + 1; ++i) {
//labelVector.f1就是label y
//聯絡下面損失函式可知，etaCoef - i * etaDelta, labelVector.f1 是 訓練資料預測值 與 實際類別 的偏差
vec[i] += weight * this.unaryLossFunc.loss(etaCoef - i * etaDelta, labelVector.f1);
}
}
return vec;
}
private double getEta(Tuple3<Double, Double, Vector> labelVector, DenseVector coefVector) {
//labelVector.f2 = {DenseVector@9972} "1.0 1.0 1.4657097546055162 1.4770978917519928"
return MatVecOp.dot(labelVector.f2, coefVector);
}
}
vec = {double[5]@10026}
0 = 219.33160199999998
1 = 198.85962729259512
2 = 179.40202828917856
3 = 160.95880498975038
4 = 143.52995739431051
回顧損失函式如下
/**
* Squared loss function.
* https://en.wikipedia.org/wiki/Loss_functions_for_classification#Square_loss
*/
public class SquareLossFunc implements UnaryLossFunc {
public SquareLossFunc() { }
@Override
public double loss(double eta, double y) {
return 0.5 * (eta - y) * (eta - y);
}
@Override
public double derivative(double eta, double y) {
return eta - y;
}
@Override
public double secondDerivative(double eta, double y) {
return 1;
}
}
3.6 UpdateModel
本模組做兩件事
- 基於dir和step length來更新coefficient，即依據方向和步長計算。
- 如果簡化理解，引數更新公式為 ：下一刻引數 = 上一時刻引數 - 學習率 * (損失函式對這個引數的導數)。
- 判斷迴圈的收斂。
因為變數太多，所以有時候就忘記誰是誰了，所以再次標示。
- OptimVariable.dir是CalcGradient計算出來的梯度做修正之後的結果
- OptimVariable.lossAllReduce 這個會變化，此時是上一階段計算的損失
程式碼邏輯大致如下：
- 1)得出"最新損失"的最小值位置pos
- 2)得出學習率 beta = dir.f1[1] / numSearchStep;
- 3)根據"最新損失pos"和 上一個最小值 last loss value判斷來進行分別處理：
- 3.1)如果所有損失都比 last loss value 大，則
- 3.1.1)縮減學習率 multiply 1.0 / (numSearchStep*numSearchStep)
- 3.1.2)把 eta 設定為 0；這個就是步長了
- 3.1.3)把當前 loss 設定為 last loss value
- 3.2)如果losses[numSearchStep] 比 last loss value 小，則
- 3.2.1)增大學習率 multiply numSearchStep
- 3.2.2)設定eta是smallest value pos，eta = beta * pos; 這個eta就是步長了
- 3.2.3)把當前 loss 設定為當前loss最小值 losses[numSearchStep]
- 3.3)否則
- 3.3.1)學習率不更改
- 3.3.2)設定eta是smallest value pos，eta = beta * pos; 這個eta就是步長了
- 3.3.3)把當前 loss 設定為當前loss最小值 losses[pos]
- 3.1)如果所有損失都比 last loss value 大，則
- 4)修正Hessian矩陣
- 5)用方向向量和步長來來更新系數向量 curCoef.f0.plusScaleEqual(dir.f0, -eta);
- 6)如果當前loss比 min loss 小，則 用 current loss 更新 min loss
- 6.1)minCoef.f1 = curCoef.f1; 更新最小loss
- 6.2)minCoef.f0.setEqual(curCoef.f0); 更新最小loss所對應的Coef，即線性模型最後需要的係數
在這裡求步長，我沒有發現 Wolf-Powell 準則的使用，Alink做了某種優化。如果有朋友能看出來Wolf-Powell如何使用，還請留言指點 ，謝謝。
程式碼如下：
public class UpdateModel extends ComputeFunction {
@Override
public void calc(ComContext context) {
double[] losses = context.getObj(OptimVariable.lossAllReduce);
Tuple2<DenseVector, double[]> dir = context.getObj(OptimVariable.dir);
Tuple2<DenseVector, double[]> pseGrad = context.getObj(OptimVariable.pseGrad);
Tuple2<DenseVector, Double> curCoef = context.getObj(OptimVariable.currentCoef);
Tuple2<DenseVector, Double> minCoef = context.getObj(OptimVariable.minCoef);
double lossChangeRatio = 1.0;
double[] lossCurve = context.getObj(OptimVariable.lossCurve);
for (int j = 0; j < losses.length; ++j) {
losses[j] /= dir.f1[0]; //dir.f1[0]是權重
}
int pos = -1;
//get the min value of losses, and remember the position.
for (int j = 0; j < losses.length; ++j) {
if (losses[j] < losses[0]) {
losses[0] = losses[j];
pos = j;
}
}
// adaptive learningRate strategy
double beta = dir.f1[1] / numSearchStep;
double eta;
// 變數如下
losses = {double[5]@10001}
0 = 35.88248934857763
1 = 49.71490682314878
2 = 44.85050707229464
3 = 40.239701247437594
4 = 35.88248934857763
dir = {Tuple2@10002} "(-9.599 -9.599 -3.5515274823133662 -3.5911942493220335,[4.0, 0.1])"
pseGrad = null
curCoef = {Tuple2@10003} "(0.001 0.0 0.0 0.0,1.7976931348623157E308)"
minCoef = {Tuple2@10004} "(0.001 0.0 0.0 0.0,1.7976931348623157E308)"
lossChangeRatio = 1.0
lossCurve = {double[100]@10005}
pos = 4
beta = 0.025
curCoef.f1 = {Double@10006} 1.7976931348623157E308
if (pos == -1) {
/**
* if all losses larger than last loss value. we'll do the below things:
* 1. reduce learning rate by multiply 1.0 / (numSearchStep*numSearchStep).
* 2. set eta with zero.
* 3. set current loss equals last loss value.
*/
eta = 0;
dir.f1[1] *= 1.0 / (numSearchStep * numSearchStep); // 學習率
curCoef.f1 = losses[0]; // 最小loss
} else if (pos == numSearchStep) {
/**
* if losses[numSearchStep] smaller than last loss value. we'll do the below things:
* 1. enlarge learning rate by multiply numSearchStep.
* 2. set eta with the smallest value pos.
* 3. set current loss equals smallest loss value.
*/
eta = beta * pos;
dir.f1[1] *= numSearchStep;
dir.f1[1] = Math.min(dir.f1[1], numSearchStep); // 學習率
lossChangeRatio = Math.abs((curCoef.f1 - losses[pos]) / curCoef.f1);
curCoef.f1 = losses[numSearchStep]; // 最小loss
// 當前數值
numSearchStep = 4 是NUM_SEARCH_STEP預設值，是線性搜尋的引數，Line search parameter, which define the search value num of one step.
lossChangeRatio = 1.0
pos = 4
eta = 0.1
curCoef.f1 = {Double@10049} 35.88248934857763
dir.f1[1] = 0.4
} else {
/**
* else :
* 1. learning rate not changed.
* 2. set eta with the smallest value pos.
* 3. set current loss equals smallest loss value.
*/
eta = beta * pos;
lossChangeRatio = Math.abs((curCoef.f1 - losses[pos]) / curCoef.f1);
curCoef.f1 = losses[pos]; // 最小loss
}
/* update loss value in loss curve at this step */
lossCurve[context.getStepNo() - 1] = curCoef.f1;
lossCurve = {double[100]@9998}
0 = 35.88248934857763
1 = Infinity
if (method.equals(OptimMethod.OWLQN)) {
.....
} else if (method.equals(OptimMethod.LBFGS)) {
Tuple2<DenseVector[], DenseVector[]> sKyK = context.getObj(OptimVariable.sKyK);
int size = dir.f0.size();
int k = context.getStepNo() - 1;
DenseVector[] sK = sKyK.f0;
for (int s = 0; s < size; ++s) {
// 修正矩陣
sK[k % OptimVariable.numCorrections].set(s, dir.f0.get(s) * (-eta));
}
curCoef.f0.plusScaleEqual(dir.f0, -eta); // 這裡是用方向向量和步長來更新系數向量
sKyK = {Tuple2@10043} "([0.9599000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0],[0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0, 0.0 0.0 0.0 0.0])"
f0 = {DenseVector[10]@10044}
0 = {DenseVector@10074} "0.9599000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034"
}
/**
* if current loss is smaller than min loss, then update the min loss and min coefficient by current.
*/
if (curCoef.f1 < minCoef.f1) {
minCoef.f1 = curCoef.f1; // 最小loss
minCoef.f0.setEqual(curCoef.f0); // 最小loss所對應的Coef，即線性模型最後需要的係數
}
curCoef = {Tuple2@9996} "(0.9609000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034,35.88248934857763)"
f0 = {DenseVector@10059} "0.9609000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034"
f1 = {Double@10048} 35.88248934857763
minCoef = {Tuple2@9997} "(0.9609000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034,35.88248934857763)"
f0 = {DenseVector@10059} "0.9609000000000001 0.9599000000000001 0.35515274823133663 0.3591194249322034"
f1 = {Double@10048} 35.88248934857763
// judge the convergence of iteration.
filter(dir, curCoef, minCoef, context, lossChangeRatio);
}
}
filter 判斷是否收斂，裡面的列印log很清晰的說明了函式邏輯。
/**
* judge the convergence of iteration.
*/
public void filter(Tuple2<DenseVector, double[]> grad,
Tuple2<DenseVector, Double> c,
Tuple2<DenseVector, Double> m,
ComContext context,
double lossChangeRatio) {
double epsilon = params.get(HasEpsilonDv0000001.EPSILON);
int maxIter = params.get(HasMaxIterDefaultAs100.MAX_ITER);
double gradNorm = ((Tuple2<DenseVector, double[]>)context.getObj(gradName)).f0.normL2();
if (c.f1 < epsilon || gradNorm < epsilon) {
printLog(" method converged at step : ", c.f1, m.f1, grad.f1[1], gradNorm, context, lossChangeRatio);
grad.f1[0] = -1.0;
} else if (context.getStepNo() > maxIter - 1) {
printLog(" method stop at max step : ", c.f1, m.f1, grad.f1[1], gradNorm, context, lossChangeRatio);
grad.f1[0] = -1.0;
} else if (grad.f1[1] < EPS) {
printLog(" learning rate is too small, method stops at step : ", c.f1, m.f1, grad.f1[1], gradNorm,
context, lossChangeRatio);
grad.f1[0] = -1.0;
} else if (lossChangeRatio < epsilon && gradNorm < Math.sqrt(epsilon)) {
printLog(" loss change ratio is too small, method stops at step : ", c.f1, m.f1, grad.f1[1], gradNorm,
context, lossChangeRatio);
grad.f1[0] = -1.0;
} else {
printLog(" method continue at step : ", c.f1, m.f1, grad.f1[1], gradNorm, context, lossChangeRatio);
}
}
3.7 OutputModel
.closeWith(new OutputModel()) 這部分是每次迭代結束，臨時輸出模型，把資料轉換成Flink通用的Row型別，Transfer the state to model rows。
public class OutputModel extends CompleteResultFunction {
@Override
public List <Row> calc(ComContext context) {
// get the coefficient of min loss.
Tuple2 <DenseVector, double[]> minCoef = context.getObj(OptimVariable.minCoef);
double[] lossCurve = context.getObj(OptimVariable.lossCurve);
int effectiveSize = lossCurve.length;
for (int i = 0; i < lossCurve.length; ++i) {
if (Double.isInfinite(lossCurve[i])) {
effectiveSize = i;
break;
}
}
double[] effectiveCurve = new double[effectiveSize];
System.arraycopy(lossCurve, 0, effectiveCurve, 0, effectiveSize);
Params params = new Params();
params.set(ModelParamName.COEF, minCoef.f0);// 重點在這裡，minCoef是我們真正需要的
params.set(ModelParamName.LOSS_CURVE, effectiveCurve);
List <Row> model = new ArrayList <>(1);
model.add(Row.of(params.toJson()));
return model;
}
}
0x04 準備模型後設資料
這裡設定了並行度為1。
// Prepare the meta info of linear model.
DataSet<Params> meta = labelInfo.f0
.mapPartition(new CreateMeta(modelName, linearModelType, isRegProc, params))
.setParallelism(1);
具體程式碼
public static class CreateMeta implements MapPartitionFunction<Object, Params> {
@Override
public void mapPartition(Iterable<Object> rows, Collector<Params> metas) throws Exception {
Object[] labels = null;
if (!this.isRegProc) {
labels = orderLabels(rows);
}
Params meta = new Params();
meta.set(ModelParamName.MODEL_NAME, this.modelName);
meta.set(ModelParamName.LINEAR_MODEL_TYPE, this.modelType);
meta.set(ModelParamName.LABEL_VALUES, labels);
meta.set(ModelParamName.HAS_INTERCEPT_ITEM, this.hasInterceptItem);
meta.set(ModelParamName.VECTOR_COL_NAME, vectorColName);
meta.set(LinearTrainParams.LABEL_COL, labelName);
metas.collect(meta);
}
}
// 變數為
meta = {Params@9667} "Params {hasInterceptItem=true, vectorColName=null, modelName="Linear Regression", labelValues=null, labelCol="label", linearModelType="LinearReg"}"
0x05 建立模型
當迭代迴圈結束之後，Alink就根據Coef資料來建立模型。
/**
* build the linear model rows, the format to be output.
*/
public static class BuildModelFromCoefs extends AbstractRichFunction implements
@Override
public void mapPartition(Iterable<Tuple2<DenseVector, double[]>> iterable,
Collector<Row> collector) throws Exception {
for (Tuple2<DenseVector, double[]> coefVector : iterable) {
LinearModelData modelData = buildLinearModelData(meta,
featureNames,
labelType,
meanVar,
hasIntercept,
standardization,
coefVector);
new LinearModelDataConverter(this.labelType).save(modelData, collector);
}
}
}
得到模型資料為，裡面coef就是 f(x)=w^Tx+b 中的 w, b。是最終用來計算的。
modelData = {LinearModelData@10584}
featureNames = {String[3]@9787}
featureTypes = null
vectorColName = null
coefVector = {DenseVector@10485} "-3.938937407856857 4.799499941426075 0.8929571907809862 1.078169576770847"
coefVectors = null
vectorSize = 3
modelName = "Linear Regression"
labelName = null
labelValues = null
linearModelType = {LinearModelType@4674} "LinearReg"
hasInterceptItem = true
lossCurve = {double[12]@10593}
0 = 35.88248934857763
1 = 12.807366842002144
2 = 0.5228366663917704
3 = 0.031112070740366038
4 = 0.01098914933042993
5 = 0.009765757443537283
6 = 0.008750523231785415
7 = 0.004210085397869248
8 = 0.0039042232755530704
9 = 0.0038821509860327537
10 = 0.003882042680010676
11 = 0.0038820422536391033
labelType = {FractionalTypeInfo@9790} "Double"
0x06 使用模型預測
預測時候，使用的是LinearModelMapper，其內部部分變數列印出來如下，能夠看出來模型資料。
this = {LinearModelMapper@10704}
vectorColIndex = -1
model = {LinearModelData@10597}
featureNames = {String[3]@10632}
featureTypes = null
vectorColName = null
coefVector = {DenseVector@10612} "-3.938937407856857 4.799499941426075 0.8929571907809862 1.078169576770847"
coefVectors = null
vectorSize = 0
modelName = "Linear Regression"
labelName = null
labelValues = {Object[0]@10613}
linearModelType = {LinearModelType@10611} "LinearReg"
hasInterceptItem = true
lossCurve = null
labelType = null
具體預測是在 LinearModelMapper.predict 中完成，具體如下：
- 對應原始輸入
Row.of("$3$0:1.0 1:7.0 2:9.0", "1.0 7.0 9.0", 1.0, 7.0, 9.0, 16.8),， - 通過
FeatureLabelUtil.getFeatureVector處理之後， - 得到的四元組是
"1.0 1.0 7.0 9.0"，其中第一個 1.0 是通過aVector.set(0, 1.0);專門設定的固定值。比如模型是 f(x) = ax + b，這個固定值 1.0 就是 b 的初始化值，隨著優化過程會得到 b。所以這裡還是需要有一個 1.0 來進行預測。 - 模型係數是：
"-3.938937407856857 4.799499941426075 0.8929571907809862 1.078169576770847"。 - 四元組 和 模型係數 點積的結果就是
dotValue = 16.814789059973744。
這樣就能看出來模型係數如何使用的了。
public Object predict(Vector vector) throws Exception {
double dotValue = MatVecOp.dot(vector, model.coefVector);
switch (model.linearModelType) {
case LR:
case SVM:
case Perceptron:
return dotValue >= 0 ? model.labelValues[0] : model.labelValues[1];
case LinearReg:
case SVR:
return dotValue;
}
}
vector = {DenseVector@10610} "1.0 1.0 7.0 9.0"
data = {double[4]@10619}
0 = 1.0
1 = 1.0
2 = 7.0
3 = 9.0
model.coefVector = {DenseVector@10612} "-3.938937407856857 4.799499941426075 0.8929571907809862 1.078169576770847"
data = {double[4]@10620}
0 = -3.938937407856857
1 = 4.799499941426075
2 = 0.8929571907809862
3 = 1.078169576770847
0x07 本系列其他文章
Alink漫談(一) : 從KMeans演算法實現不同看Alink設計思想
Alink漫談(二) : 從原始碼看機器學習平臺Alink設計和架構
Alink漫談(八) : 二分類評估 AUC、K-S、PRC、Precision、Recall、LiftChart 如何實現
0xFF 參考
http://www.mamicode.com/info-detail-2508527.html)
機器學習演算法實現解析——liblbfgs之L-BFGS演算法
LogisticRegressionWithLBFGS--邏輯迴歸
[原創]用“人話”解釋不精確線搜尋中的Armijo-Goldstein準則及Wolfe-Powell準則
https://www.zhihu.com/question/36425542
https://zhuanlan.zhihu.com/p/32821110
CRF L-BFGS Line Search原理及程式碼分析