Alink漫談(十) ：線性迴歸實現之資料預處理

羅西的思考發表於2020-07-11

Alink漫談(十) ：線性迴歸實現之資料預處理

0x00 摘要

Alink 是阿里巴巴基於實時計算引擎 Flink 研發的新一代機器學習演算法平臺，是業界首個同時支援批式演算法、流式演算法的機器學習平臺。本文和下文將介紹線性迴歸在Alink中是如何實現的，希望可以作為大家看線性迴歸程式碼的Roadmap。

因為Alink的公開資料太少，所以以下均為自行揣測，肯定會有疏漏錯誤，希望大家指出，我會隨時更新。

本系列目前已有十篇，歡迎大家指點

0x01 概念

1.1 線性迴歸

線性迴歸是利用數理統計中迴歸分析，來確定兩種或兩種以上變數間相互依賴的定量關係的一種統計分析方法，運用十分廣泛。其表達形式為y = w'x+e，e為誤差服從均值為0的正態分佈。

線上性迴歸中，目標值與特徵之間存在著線性相關的關係。即假設這個方程是一個線性方程，一個多元一次方程。

基本形式：給定由 d 個屬性描述的示例，線性模型試圖學得一個通過屬性的線性組合來進行預測的函式，即：

\[f(x)=w_1x_1 +w_2x_2 ... +w_dx_d+b \]

其中w為引數，也稱為權重，可以理解為x1，x2...和 xd 對f(x)的影響度。

一般形式為：

\[f(x)=w^Tx+b \]

假如我們依據這個公式來預測 f(x)，公式中的x是我們已知的，然而w，b的取值卻不知道，只要我們把w，b的取值求解出來，模型就得以確定。我們就可以依據這個公式來做預測了。

那麼如何依據訓練資料求解 w 和 b 的最優取值呢？關鍵是衡量 f 和 y 之間的差別。這就牽扯到另外一個概念：損失函式（Loss Function）。

1.2 優化模型

假如有一個模型 f(x)，如何判斷這個模型是否優秀？這種定性的判斷可以通過一個成為經驗誤差風險的數值來進行衡量，也就是模型 f 在所有訓練樣本上所犯錯誤的總和 E(x)。

我們通過在訓練集上最小化經驗損失來訓練模型。換言之，通過調節 f 的引數 w，使得經驗誤差風險 E(x) 不斷下降，最終達到最小值的時候，我們就獲得了一個 “最優” 的模型。

但是如果按照上面的定義，E(x) 是一組示性函式的和，因此是不連續不可導的函式，不易優化。為了解決這個問題，人們提出了“損失函式”的概念。損失函式就是和誤差函式有一定關係（比如是誤差函式的上界），但是具有更好的數學性質（比如連續，可導，凸性等），比較容易進行優化。所以我們就可以對損失函式來優化。

損失函式如果連續可導，所以我們可以用梯度下降法等一階演算法，也可以用牛頓法，擬牛頓法等二階演算法。當優化演算法收斂後，我們就得到一個不錯的模型。如果損失函式是一個凸函式，我們就可以得到最優模型。

典型的優化方法：

	一階演算法	二階演算法
確定性演算法	梯度下降法投影次梯度下降近端梯度下降 Frank-Wolfe演算法 Nesterov加速演算法座標下降法對偶座標上升法	牛頓法，擬牛頓法
隨機演算法	隨機梯度下降法隨機座標下降法隨機對偶座標上升法隨機方差減小梯度法	隨機擬牛頓法

所以我們可以知道，優化LinearRegression模型 f 的手段一定是：確定損失函式，用 x，y 作為輸入訓練以求得損失函式最小值，從而確定 f 的引數 w。過程大致如下：

處理輸入，把 x, y 轉換成演算法需要的格式。
找一個合適的預測函式，一般表示為 h 函式，該函式就是我們需要找的分類函式，它用來預測輸入資料的判斷結果。
構造一個Cost函式（損失函式），該函式表示預測的輸出（h）與訓練資料類別（y）之間的偏差，可以是二者之間的差（h-y）或者是其他的形式。綜合考慮所有訓練資料的 “損失”，將Cost求和或者求平均，記為J(θ)函式，表示所有訓練資料預測值與實際類別的偏差。
顯然，損失函式 J(θ) 函式的值越小表示預測函式越準確（即h函式越準確），所以這一步需要做的是找到 J(θ) 函式的最小值。注意，損失函式是關於 θ 的函式！也就是說，對於損失函式來講，θ不再是函式的引數，而是損失函式的自變數！
準備模型後設資料，建立模型。

1.3 損失函式&目標函式

先概括說明：

損失函式：計算的是一個樣本的誤差；
代價函式：是整個訓練集上所有樣本誤差的平均，經常和損失函式混用；
目標函式：代價函式 + 正則化項；

再詳細闡釋：

假設我們用 f(X) 來擬合真實值Y。這個輸出的f(X)與真實值Y可能是相同的，也可能是不同的，為了表示我們擬合的好壞，我們就用一個函式來度量擬合的程度。這個函式就稱為損失函式(loss function)，或者叫代價函式(cost function)。

損失函式用來衡量演算法的執行情況，估量模型的預測值與真實值的不一致程度，是一個非負實值函式，通常使用 L(Y,f(x)) 來表示。損失函式越小，模型的魯棒性就越好。損失函式是經驗風險函式的核心部分。

目標函式是一個相關但更廣的概念，對於目標函式來說在有約束條件下的最小化就是損失函式（loss function）。

因為f(x)可能會過度學習歷史資料，導致它在真正預測時效果會很不好，這種情況稱為過擬合(over-fitting)。這樣得到的函式會過於複雜。所以我們不僅要讓經驗風險最小化，還要讓結構風險最小化。這個時候就定義了一個函式 J(x)，這個函式專門用來度量模型的複雜度，在機器學習中也叫正則化(regularization)。常用的有 L1, L2範數。

L1 正則的本質是為模型增加了“模型引數服從零均值拉普拉斯分佈”這一先驗知識。

L2 正則的本質是為模型增加了“模型引數服從零均值正態分佈”這一先驗知識。

L1 正則化增加了所有權重 w 引數的絕對值之和逼迫更多 w 為零，也就是變稀疏（ L2 因為其導數也趨 0, 奔向零的速度不如 L1 給力了）。L1 正則化的引入就是為了完成特徵自動選擇的光榮使命，它會學習地去掉無用的特徵，也就是把這些特徵對應的權重置為 0。

L2 正則化中增加所有權重 w 引數的平方之和，逼迫所有 w 儘可能趨向零但不為零（L2 的導數趨於零）。因為在未加入 L2 正則化發生過擬合時，擬合函式需要顧忌每一個點，最終形成的擬合函式波動很大，在某些很小的區間裡，函式值的變化很劇烈，也就是某些 w 值非常大。為此，L2 正則化的加入就懲罰了權重變大的趨勢。

到這一步我們就可以說我們最終的優化函式是：min(L(Y, f(x) + J(x)) ，即最優化經驗風險和結構風險，而這個函式就被稱為目標函式。

在迴歸問題中，通過目標函式來求解最優解，常用的是平方誤差（最小二乘線性迴歸）代價函式。損失函式則是平方損失函式。

1.4 最小二乘法

均方誤差是迴歸任務中最常用的效能度量，因此可以使均方誤差最小。基於均方誤差最小化來進行模型求解的方法稱為“最小二乘法”。線上性迴歸中，最小二乘法就是找到一條直線，使所有樣本到直線的 "歐式距離和" 最小。於是線性迴歸中損失函式就是平方損失函式。

有了這些基礎概念，下面我們就開始動手分析Alink的程式碼。

0x02 示例程式碼

首先，我們給出線性迴歸的示例。

public class LinearRegressionExample {
    static Row[] vecrows = new Row[] {
            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),
            Row.of("$3$0:1.0 1:3.0 2:3.0", "1.0 3.0 3.0", 1.0, 3.0, 3.0, 6.7),
            Row.of("$3$0:1.0 1:2.0 2:4.0", "1.0 2.0 4.0", 1.0, 2.0, 4.0, 6.9),
            Row.of("$3$0:1.0 1:3.0 2:4.0", "1.0 3.0 4.0", 1.0, 3.0, 4.0, 8.0)
    };
    static String[] veccolNames = new String[] {"svec", "vec", "f0", "f1", "f2", "label"};
    static BatchOperator vecdata = new MemSourceBatchOp(Arrays.asList(vecrows), veccolNames);
    static StreamOperator svecdata = new MemSourceStreamOp(Arrays.asList(vecrows), veccolNames);

    public static void main(String[] args) throws Exception {
        String[] xVars = new String[] {"f0", "f1", "f2"};
        String yVar = "label";
        String vec = "vec";
        String svec = "svec";
        LinearRegression linear = new LinearRegression()
                .setLabelCol(yVar)  // 這裡把變數都設定好了，後續會用到
                .setFeatureCols(xVars)
                .setPredictionCol("linpred");

        Pipeline pl = new Pipeline().add(linear);
        PipelineModel model = pl.fit(vecdata);

        BatchOperator result = model.transform(vecdata).select(
                new String[] {"label", "linpred"});

        List<Row> data = result.collect();
    }
}

輸出是

svec|vec|f0|f1|f2|label|linpred
----|---|--|--|--|-----|-------
$3$0:1.0 1:7.0 2:9.0|1.0 7.0 9.0|1.0000|7.0000|9.0000|16.8000|16.8148
$3$0:1.0 1:3.0 2:4.0|1.0 3.0 4.0|1.0000|3.0000|4.0000|8.0000|7.8521
$3$0:1.0 1:3.0 2:3.0|1.0 3.0 3.0|1.0000|3.0000|3.0000|6.7000|6.7739
$3$0:1.0 1:2.0 2:4.0|1.0 2.0 4.0|1.0000|2.0000|4.0000|6.9000|6.959

根據前文我們可以知道，在迴歸問題中，通過優化目標函式來求解最優解，常用的是平方誤差（最小二乘線性迴歸）代價函式。損失函式則是平方損失函式。

對應到Alink，優化函式或者優化器是擬牛頓法的L-BFGS演算法，目標函式是UnaryLossObjFunc，損失函式是SquareLossFunc。線性迴歸訓練總體邏輯是LinearRegTrainBatchOp。所以我們下面一一論述。

0x03 整體概述

LinearRegression 訓練用到LinearRegTrainBatchOp，而LinearRegTrainBatchOp的基類是BaseLinearModelTrainBatchOp。所以我們來看BaseLinearModelTrainBatchOp。

public class LinearRegression extends Trainer <LinearRegression, LinearRegressionModel> implements LinearRegTrainParams <LinearRegression>, LinearRegPredictParams <LinearRegression> {
   @Override
   protected BatchOperator train(BatchOperator in) {
      return new LinearRegTrainBatchOp(this.getParams()).linkFrom(in);
   }
}

BaseLinearModelTrainBatchOp.linkFrom 程式碼如下，註釋中給出了清晰的邏輯：

大體是：

獲取演算法引數，label資訊；
準備，轉換資料到 Tuple3 format <weight, label, feature vector>；
獲得統計資訊，比如向量大小，均值和方差；
對訓練資料做標準化和插值；
使用L-BFGS演算法，通過對損失函式求最小值從而對模型優化；
準備模型後設資料；
建立模型；

public T linkFrom(BatchOperator<?>... inputs) {
    BatchOperator<?> in = checkAndGetFirst(inputs);
    // Get parameters of this algorithm.
    Params params = getParams();
    // Get type of processing: regression or not
    boolean isRegProc = getIsRegProc(params, linearModelType, modelName);
    // Get label info : including label values and label type.
    Tuple2<DataSet<Object>, TypeInformation> labelInfo = getLabelInfo(in, params, isRegProc);
    // Transform data to Tuple3 format.//weight, label, feature vector.
    DataSet<Tuple3<Double, Double, Vector>> initData = transform(in, params, labelInfo.f0, isRegProc);
    // Get statistics variables : including vector size, mean and variance of train data.
    Tuple2<DataSet<Integer>, DataSet<DenseVector[]>>
        statInfo = getStatInfo(initData, params.get(LinearTrainParams.STANDARDIZATION));
    // Do standardization and interception to train data.
    DataSet<Tuple3<Double, Double, Vector>> trainData = preProcess(initData, params, statInfo.f1);
    // Solve the optimization problem.
    DataSet<Tuple2<DenseVector, double[]>> coefVectorSet = optimize(params, statInfo.f0,
        trainData, linearModelType, MLEnvironmentFactory.get(getMLEnvironmentId()));
    // Prepare the meta info of linear model.
    DataSet<Params> meta = labelInfo.f0
        .mapPartition(new CreateMeta(modelName, linearModelType, isRegProc, params))
        .setParallelism(1);
    // Build linear model rows, the format to be output.
    DataSet<Row> modelRows;
    String[] featureColTypes = getFeatureTypes(in, params.get(LinearTrainParams.FEATURE_COLS));
    modelRows = coefVectorSet
        .mapPartition(new BuildModelFromCoefs(labelInfo.f1,
            params.get(LinearTrainParams.FEATURE_COLS),
            params.get(LinearTrainParams.STANDARDIZATION),
            params.get(LinearTrainParams.WITH_INTERCEPT), featureColTypes))
        .withBroadcastSet(meta, META)
        .withBroadcastSet(statInfo.f1, MEAN_VAR)
        .setParallelism(1);
    // Convert the model rows to table.
    this.setOutput(modelRows, new LinearModelDataConverter(labelInfo.f1).getModelSchema());
    return (T)this;
}

我們後續還會對此邏輯進行細化。

0x04 基礎功能

我們首先介紹下相關基礎功能和相關概念，比如損失函式，目標函式，梯度等。

4.1 損失函式

損失函式涉及到若干概念。

4.1.1 導數和偏導數

導數也是函式，是函式的變化率與位置的關係。導數代表了在自變數變化趨於無窮小的時候，函式值的變化與自變數的變化的比值。幾何意義是這個點的切線。物理意義是該時刻的（瞬時）變化率。

導數反映的是函式y=f(x)在某一點處沿x軸正方向的變化率。直觀地看，也就是在x軸上某一點處，如果f’(x)>0，說明f(x)的函式值在x點沿x軸正方向是趨於增加的；如果f’(x)<0，說明f(x)的函式值在x點沿x軸正方向是趨於減少的。

一元導數表徵的是：一元函式 f(x)與自變數 x 在某點附近變化的比率（變化率，斜率）。

如果是多元函式呢？則為偏導數。偏導數是多元函式“退化”成一元函式時的導數，這裡“退化”的意思是固定其他變數的值，只保留一個變數，依次保留每個變數，則N元函式有N個偏導數。偏導數為函式在每個位置處沿著自變數座標軸方向上的導數（切線斜率）。二元函式的偏導數表徵的是：函式 F(x,y) 與自變數 x（或y）在某點附近變化的比率（變化率）。

4.1.2 方向導數

導數和偏導數的定義中，均是沿座標軸正方向討論函式的變化率。那麼當我們討論函式沿任意方向的變化率時，也就引出了方向導數的定義，即：某一點在某一趨近方向上的導數值。

方向導數就是偏導數合成向量與方向向量的內積。方向導數的本質是一個數值，簡單來說其定義為：一個函式沿指定方向的變化率。

4.1.3 Hessian矩陣

在一元函式求解的問題中，我們可以很愉快的使用牛頓法求駐點。但在機器學習的優化問題中，我們要優化的都是多元函式，x往往不是一個實數，而是一個向量，所以將牛頓求根法利用到機器學習中時，x 是一個向量, y 也是一個向量，對 x 求導以後得到的是一個矩陣，就是Hessian矩陣。

在數學中，海森矩陣（Hessian matrix 或 Hessian）是一個自變數為向量的實值函式的二階偏導陣列成的方塊矩陣。多元函式的二階導數就是一個海森矩陣。

4.1.4 平方損失函式 in Alink

前面提到，線性迴歸中損失函式就是平方損失函式。我們來看看實現。後續實現將呼叫此類的 loss 和 derivative，具體遇到時候再講。

UnaryLossFunc是介面，代表一元損失函式。它定義的每個函式都有兩個輸入 (eta and y)，Alink把這兩個輸入的差作為損失函式的一元變數。基本API是求損失，求導數，求二階導數。

public interface UnaryLossFunc extends Serializable {
	// Loss function.
	double loss(double eta, double y);
	// The derivative of loss function.
	double derivative(double eta, double y);
	// The second derivative of the loss function.
	double secondDerivative(double eta, double y);
}

平方損失函式具體實現如下：

public class SquareLossFunc implements UnaryLossFunc {

   @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;
   }
}

4.2 目標函式

這裡涉及的概念是梯度，梯度下降法。

4.2.1 梯度

對於模型優化，我們要選擇最優的 θ，使得 f(x) 最接近真實值。這個問題就轉化為求解最優的 θ，使損失函式 J(θ) 取最小值。那麼如何解決這個轉化後的問題呢？這又牽扯到一個概念：梯度下降（Radient Descent）。

所以我們首先要溫習下梯度。

向量的定義是有方向（direction）有大小（magnitude）的量。
梯度其實是一個向量，即有方向有大小；其定義為：一個多元函式對於其自變數分別求偏導數，這些偏導數所組成的向量就是函式的梯度。
梯度即函式在某一點最大的方向導數，函式沿梯度方向函式有最大的變化率。
梯度的第一層含義就是“方向導數的最大值”。
當前位置的梯度方向，為函式在該位置處方向導數最大的方向，也是函式值上升最快的方向，反方向為下降最快的方向；
梯度的幾何含義就是：沿向量所在直線的方向變化率最大。

4.2.2 梯度下降法

梯度下降法是一個一階最優化演算法，它的核心思想是：要想最快找到一個函式的區域性極小值，必須沿函式當前點對應“梯度”（或者近似梯度）的反方向（下降）進行規定步長“迭代”搜尋。沿梯度（斜率）的反方向移動，這就是“梯度下降法”。

既然在變數空間的某一點處，函式沿梯度方向具有最大的變化率，那麼在優化目標函式的時候，自然是沿著負梯度方向去減小函式值，以此達到我們的優化目標。

梯度下降中的下降，意思是讓函式的未知數隨著梯度的方向運動。什麼是梯度的方向呢？把這一點帶入到梯度函式中，結果為正，那我們就把這一點的值變小一些，同時就是讓梯度變小些；當這一點帶入梯度函式中的結果為負的時候，就給這一點的值增大一些。

如何沿著負梯度方向減小函式值呢？既然梯度是偏導數的集合，同時梯度和偏導數都是向量，那麼參考向量運演算法則，我們在每個變數軸上減小對應變數值即可。

梯度下降就是讓梯度中所有偏導函式都下降到最低點的過程.(劃重點:下降)。都下降到最低點了,那每個未知數(或者叫維度)的最優解就得到了，所以他是解決函式最優化問題的演算法。

“最小二乘法”和“梯度下降法”，前者用於“搜尋最小誤差”，後者用於“用最快的速度搜尋”，二者常常配合使用。對最小二乘法的引數調優就轉變為了求這個二元函式的極值問題，也就是說可以應用“梯度下降法”了。

在最小二乘函式中，已擁有的條件是一些樣本點和樣本點的結果，就是矩陣X和每一條X樣本的lable值y。X是矩陣，y是向量。所以我們要知道，梯度下降中求偏導數的未知數不是x和y，而是x的引數w。

4.2.3 目標函式 in Alink

目標函式的基類是OptimObjFunc，其提供API 比如計算梯度，損失，hessian矩陣，以及依據取樣點更新梯度和hessian矩陣。其幾個派生類如下，從註釋中可以看到使用範圍。

我們可以看到正則化(regularization) L1, L2範數，這是相比損失函式增加的模組。

public abstract class OptimObjFunc implements Serializable {
    protected final double l1;
    protected final double l2; // 正則化(regularization) L1, L2範數。
    protected Params params;   
  .....
}

// Unary loss object function.
public class UnaryLossObjFunc extends OptimObjFunc

// The OptimObjFunc for multilayer perceptron.
public class AnnObjFunc extends OptimObjFunc
  
// Accelerated failure time Regression object function.
public class AftRegObjFunc extends OptimObjFunc

// Softmax object function.
public class SoftmaxObjFunc extends OptimObjFunc

對於線性模型，BaseLinearModelTrainBatchOp 中會根據模型型別來生成目標函式，可以看到在生成目標函式同時，也相應設定了不同的損失函式，其中 SquareLossFunc 就是我們之前提到的。

public static OptimObjFunc getObjFunction(LinearModelType modelType, Params params) {
    OptimObjFunc objFunc;
    // For different model type, we must set corresponding loss object function.
    switch (modelType) {
        case LinearReg:
            // 我們這裡！
            objFunc = new UnaryLossObjFunc(new SquareLossFunc(), params);
            break;
        case SVR:
            double svrTau = params.get(LinearSvrTrainParams.TAU);
            objFunc = new UnaryLossObjFunc(new SvrLossFunc(svrTau), params);
            break;
        case LR:
            objFunc = new UnaryLossObjFunc(new LogLossFunc(), params);
            break;
        case SVM:
            objFunc = new UnaryLossObjFunc(new SmoothHingeLossFunc(), params);
            break;
        case Perceptron:
            objFunc = new UnaryLossObjFunc(new PerceptronLossFunc(), params);
            break;
        case AFT:
            objFunc = new AftRegObjFunc(params);
            break;
        default:
            throw new RuntimeException("Not implemented yet!");
    }
    return objFunc;
}

4.2.4 一元目標函式 in Alink

一元目標函式就是我們線性迴歸用到的目標函式，其只有一個新增變數：unaryLossFunc。就是一元損失函式。

/**
 * Unary loss object function.
 */
public class UnaryLossObjFunc extends OptimObjFunc {
    private UnaryLossFunc unaryLossFunc;
}

一元目標函式提供了很多功能，我們這裡用到主要是：

calcGradient ：根據一組取樣點計算梯度，這是從基類OptimObjFunc整合的。
updateGradient ：根據一個取樣點更新梯度；
calcSearchValues ：為線性搜尋計算損失；

4.2.4.1 依據一組取樣點計算梯度

對於本文，這裡更新的是損失函式的梯度。

再次囉嗦下，損失函式用來度量擬合的程度，從而評估模型擬合的好壞，記為 J(θ)。注意，損失函式是關於 θ 的函式！也就是說，對於損失函式來講，θ不再是函式的引數，而是損失函式的自變數！

當我們計算損失時，是將每個樣本中的特徵 xi 和對應的目標變數真實值 yi 帶入損失函式，此時，損失函式中就只剩下 θ 是未知的。

損失函式的梯度即對 θi 求偏導，由於損失函式是關於 θ 的函式，因此，θ 的取值不同，得出來的的梯度向量也是不同的。借用“下山”的比喻來解釋，θ 的不同取值，相當於處於山上的不同位置，每一個位置都會計算出一個梯度向量▽J(θ)。

這裡的 l1, l2 就是之前提到的正則化(regularization) L1, L2範數。

/**
 * Calculate gradient by a set of samples.
 *
 * @param labelVectors train data.
 * @param coefVector   coefficient of current time.
 * @param grad         gradient.
 * @return weight sum
 */
public double calcGradient(Iterable<Tuple3<Double, Double, Vector>> labelVectors,
                           DenseVector coefVector, DenseVector grad) {
    double weightSum = 0.0;
    for (int i = 0; i < grad.size(); i++) {
        grad.set(i, 0.0);
    }
  
// 對輸入的樣本集合labelVectors逐個計算梯度  
    for (Tuple3<Double, Double, Vector> labelVector : labelVectors) {
        if (labelVector.f2 instanceof SparseVector) {
           ((SparseVector)(labelVector.f2)).setSize(coefVector.size());
        }
      
// 以這個樣本為例 
labelVector = {Tuple3@9895} "(1.0,16.8,1.0 1.0 1.4657097546055162 1.4770978917519928)"
 f0 = {Double@9903} 1.0
 f1 = {Double@9904} 16.8
 f2 = {DenseVector@9905} "1.0 1.0 1.4657097546055162 1.4770978917519928"
  
        weightSum += labelVector.f0; // labelVector.f0是權重
        updateGradient(labelVector, coefVector, grad);
    }
    if (weightSum > 0.0) {
        grad.scaleEqual(1.0 / weightSum);
    }
// l2正則化    
    if (0.0 != this.l2) {
        grad.plusScaleEqual(coefVector, this.l2 * 2);
    }
// l1正則化   
    if (0.0 != this.l1) {
        double[] coefArray = coefVector.getData();
        for (int i = 0; i < coefVector.size(); i++) {
            grad.add(i, Math.signum(coefArray[i]) * this.l1);
        }
    }
    return weightSum;
}

4.2.4.2 根據一個取樣點更新梯度

這裡 labelVector.f0是權重，labelVector.f1是 y，labelVector.f2是 x-vec 四維向量，coefVector是w係數向量。

getEta是點積，即 x向量與當前w係數的點積，就是當前計算的 y。
labelVector.f0 * unaryLossFunc.derivative(eta, labelVector.f1); 就是呼叫SquareLossFunc.derivative 函式來計算一階導數。
updateGrad.plusScaleEqual(labelVector.f2, div); 就是在原有梯度基礎上更新梯度

public class UnaryLossObjFunc extends OptimObjFunc {
    /**
     * Update gradient by one sample.
     *
     * @param labelVector a sample of train data.
     * @param coefVector  coefficient of current time.
     * @param updateGrad  gradient need to update.
     */
    @Override
    protected void updateGradient(Tuple3<Double, Double, Vector> labelVector, DenseVector coefVector, DenseVector updateGrad) {
        // 點積，就是當前計算出來的y
        double eta = getEta(labelVector, coefVector); 
        // 一階導數。labelVector.f0是權重
        double div = labelVector.f0 * unaryLossFunc.derivative(eta, labelVector.f1); 
        // 點乘之後還需要相加。labelVector.f2 就是x—vec，比如 1.0 1.0 1.4657097546055162 1.4770978917519928
        updateGrad.plusScaleEqual(labelVector.f2, div);
    }
  
    private double getEta(Tuple3<Double, Double, Vector> labelVector, DenseVector coefVector) {
        // 點積，表示第 i 次迭代中節點上的第 k 個特徵向量與特徵權重分量的點乘。coefVector中第 c 項表示為第 i 次迭代中特徵權重向量在第 c 列節點上的分量
        return MatVecOp.dot(labelVector.f2, coefVector);
    }
}

/**
* Plus with another vector scaled by "alpha".
*/
public void plusScaleEqual(Vector other, double alpha) {
	if (other instanceof DenseVector) {
		BLAS.axpy(alpha, (DenseVector) other, this);
	} else {
		BLAS.axpy(alpha, (SparseVector) other, this);
	}
}

4.3 優化函式

Alink中提供給了一系列並行優化函式，比如GD, SGD, LBFGS, OWLQN, NEWTON method。

其基類是Optimizer。

public abstract class Optimizer {
    protected final DataSet<?> objFuncSet; // 具體目標函式，計算梯度和損失
    protected final DataSet<Tuple3<Double, Double, Vector>> trainData; //訓練資料
    protected final Params params; //引數
    protected DataSet<Integer> coefDim; //dimension of features.
    protected DataSet<DenseVector> coefVec = null; //最終係數w
    .......
}

線性迴歸主要用到了LBFGS演算法。

public class Lbfgs extends Optimizer

具體呼叫如下

public static DataSet<Tuple2<DenseVector, double[]>> optimize(.....) {
    // Loss object function
    DataSet<OptimObjFunc> objFunc = session.getExecutionEnvironment()
        .fromElements(getObjFunction(modelType, params));

    if (params.contains(LinearTrainParams.OPTIM_METHOD)) {
        LinearTrainParams.OptimMethod method = params.get(LinearTrainParams.OPTIM_METHOD);
        return OptimizerFactory.create(objFunc, trainData, coefficientDim, params, method).optimize();
    } else if (params.get(HasL1.L_1) > 0) {
        return new Owlqn(objFunc, trainData, coefficientDim, params).optimize();
    } else {
        // 我們的程式將執行到這裡
        return new Lbfgs(objFunc, trainData, coefficientDim, params).optimize();
    }
}

機器學習基本優化套路是：

準備資料 ----> 優化函式 ----> 目標函式 ----> 損失函式

對應我們這裡是

BaseLinearModelTrainBatchOp.linkFrom(整體邏輯) -----> Lbfgs(繼承Optimizer) ----> UnaryLossObjFunc(繼承OptimObjFunc) ----> SquareLossFunc(繼承UnaryLossFunc)

0x05 資料準備

看完完底層功能，我們再次回到線性迴歸總體流程。

總結 BaseLinearModelTrainBatchOp.linkFrom 的基本流程如下：(發現某些媒體對於列表排版支援不好，所以加上序號)。

首先再給出輸入一個例子：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),這裡後面 4 項對應列名是 "f0", "f1", "f2", "label"。

1)獲取到label的資訊，包括label數值和種類。 labelInfo = getLabelInfo() 這裡有一個 distinct 操作，所以會去重。最後得到label的可能取值範圍：0，1，型別是 Double。
2)用transform函式把輸入轉換成三元組Tuple3<weight, label, feature vector>。具體說，會把輸入中的三個特徵"f0", "f1", "f2" 轉換為一個向量 vec, 我們以後稱之為x-vec。重點就在於特徵變成了一個向量。所以這個三元組可以認為是 <權重, y-value, x-vec>。
3)用statInfo = getStatInfo() 獲取統計變數，包括vector size, mean和variance。這裡流程比較複雜。
- 3.1)用trainData.map{return value.f2;}來獲取訓練資料中的 x-vec。
- 3.2)呼叫StatisticsHelper.summary來對 x-vec 做處理
  - 3.2.1)呼叫 summarizer
    - 3.2.1.1)呼叫 mapPartition(new VectorSummarizerPartition(bCov))
      - 3.2.1.1.1)呼叫VectorSummarizerPartition.mapPartition，其遍歷列表，列表中的每一個變數 sv 是 x-vec。srt = srt.visit(sv)，會根據每一個新輸入重新計算count，sum，squareSum，normL1..，這樣就得到了本partiton中輸入每列的這些統計數值。
    - 3.2.1.2)呼叫 reduce(VectorSummarizerUtil.merge(value1, value2)) 來歸併每一個partition的結果。
  - 3.2.2)呼叫map(BaseVectorSummarizer summarizer)，其實呼叫到DenseVectorSummarizer，就是生成一個DenseVectorSummary向量，裡面是count，sum，squareSum，normL1，min，max，numNonZero。
- 3.3)呼叫 coefficientDim = summary.map
- 3.4)呼叫 meanVar = coefficientDim.map，最後得到 Tuple2.of(coefficientDim, meanVar)
4)preProcess(initData, params, statInfo.f1) 用3) 計算的結果對輸入資料做標準化和插值 standardization and interception。上面得到的 meanVar 將會作為引數傳入。這裡是對 x-vec 做標準化。比如原始輸入Row是"(1.0,16.8,1.0 7.0 9.0)"，其中 x-vec 是"1.0 7.0 9.0"，進行標準化之後，x-vec 變成了 4 項：{ 第1項是固定值 "1.0 "，所以4 項是 "1.0 1.0 1.4657097546055162 1.4770978917519928" }，所以轉換後的Row是"(1.0,16.8,1.0 1.0 1.4657097546055162 1.4770978917519928)"。即weight 是1.0，y-value是16.8，後續4個是x-vec。
以上完成了對資料的處理。
5)呼叫 optimize(params, statInfo.f0, trainData, linearModelType) 通過對損失函式求最小值從而對模型優化。(使用L-BFGS演算法，會單獨拿出來講解)
6)呼叫 mapPartition(new CreateMeta()) 來準備模型後設資料。
7)呼叫 mapPartition(new BuildModelFromCoefs) 來建立模型。

可以看到，資料準備佔據了很大部分，下面我們看看資料準備的幾個步驟。

5.1 獲取label資訊

此處程式碼對應上面基本流程的 1)

因為之前有一個distinct操作，所以會去重。最後得到label的可能取值範圍：0，1，型別是 Double。

private Tuple2<DataSet<Object>, TypeInformation> getLabelInfo(BatchOperator in,
                                                                  Params params,
                                                                  boolean isRegProc) {
        String labelName = params.get(LinearTrainParams.LABEL_COL);
        // Prepare label values
        DataSet<Object> labelValues;
        TypeInformation<?> labelType = null;
        if (isRegProc) {
            // 因為是迴歸，所以是這裡
            labelType = Types.DOUBLE;
            labelValues = MLEnvironmentFactory.get(in.getMLEnvironmentId())
                .getExecutionEnvironment().fromElements(new Object());
        } else {
          .....
        }
        return Tuple2.of(labelValues, labelType);
}

5.2 把輸入轉換成三元組

此處程式碼對應上面基本流程的 2) 。

用transform函式把輸入轉換成三元組Tuple3<weight, label, feature vector>。具體說，會把輸入中的三個特徵"f0", "f1", "f2" 轉換為一個向量 vec, 我們以後稱之為x-vec。重點就在於特徵變成了一個向量。所以這個三元組可以認為是 <權重, y-value, x-vec>。

private DataSet<Tuple3<Double, Double, Vector>> transform(BatchOperator in,
                                                          Params params,
                                                          DataSet<Object> labelValues,
                                                          boolean isRegProc) {
    ......
    // 獲取Schema
    TableSchema dataSchema = in.getSchema();
    // 獲取各種index 
    int labelIdx = TableUtil.findColIndexWithAssertAndHint(dataSchema.getFieldNames(), labelName);
    ......
    int weightIdx = weightColName != null ? TableUtil.findColIndexWithAssertAndHint(in.getColNames(), weightColName) : -1;
    int vecIdx = vectorColName != null ? TableUtil.findColIndexWithAssertAndHint(in.getColNames(), vectorColName) : -1;
    // 用transform函式把輸入轉換成三元組Tuple3<weight, label, feature vector>
    return in.getDataSet().map(new Transform(isRegProc, weightIdx, vecIdx, featureIndices, labelIdx)).withBroadcastSet(labelValues, LABEL_VALUES);
}

這裡對應的變數列印出來為

params = {Params@2745} "Params {featureCols=["f0","f1","f2"], labelCol="label", predictionCol="linpred"}"
labelValues = {DataSource@2845} 
isRegProc = true
featureColNames = {String[3]@2864} 
 0 = "f0"
 1 = "f1"
 2 = "f2"
labelName = "label"
weightColName = null
vectorColName = null
dataSchema = {TableSchema@2866} "root\n |-- svec: STRING\n |-- vec: STRING\n |-- f0: DOUBLE\n |-- f1: DOUBLE\n |-- f2: DOUBLE\n |-- label: DOUBLE\n"
featureIndices = {int[3]@2878} 
 0 = 2
 1 = 3
 2 = 4
labelIdx = 5
weightIdx = -1
vecIdx = -1

具體在runtime時候，會進入到Transform.map函式。我們可以看到，會把輸入中的三個特徵"f0", "f1", "f2",轉換為一個向量 vec, 我們以後稱之為x-vec。

private static class Transform extends RichMapFunction<Row, Tuple3<Double, Double, Vector>> {
    @Override
    public Tuple3<Double, Double, Vector> map(Row row) throws Exception {
        // 獲取權重
        Double weight = weightIdx != -1 ? ((Number)row.getField(weightIdx)).doubleValue() : 1.0;
        // 獲取label
        Double val = FeatureLabelUtil.getLabelValue(row, this.isRegProc,
            labelIdx, this.positiveLableValueString);
        if (featureIndices != null) {
            // 獲取x-vec
            DenseVector vec = new DenseVector(featureIndices.length);
            for (int i = 0; i < featureIndices.length; ++i) {
                vec.set(i, ((Number)row.getField(featureIndices[i])).doubleValue());
            }
            // 構建三元組
            return Tuple3.of(weight, val, vec);
        } else {
            Vector vec = VectorUtil.getVector(row.getField(vecIdx));
            return Tuple3.of(weight, val, vec);
        }
    }
}

如果對應原始輸入 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), ，則程式中各種變數為：

row = {Row@9723} "$3$0:1.0 1:7.0 2:9.0,1.0 7.0 9.0,1.0,7.0,9.0,16.8"
weight = {Double@9724} 1.0
val = {Double@9725} 16.8
vec = {DenseVector@9729} "1.0 7.0 9.0"
vecIdx = -1
featureIndices = {int[3]@9726} 
 0 = 2
 1 = 3
 2 = 4

5.3 獲取統計變數

用getStatInfo() 對輸入資料做標準化和插值 standardization and interception。

此處程式碼對應上面基本流程的 3)

用statInfo = getStatInfo() 獲取統計變數，包括vector size, mean和variance。這裡流程比較複雜。

3.1)用trainData.map{return value.f2;}來獲取訓練資料中的 x-vec。
3.2)呼叫StatisticsHelper.summary來對 x-vec 做處理
- 3.2.1)呼叫 summarizer
  - 3.2.1.1)呼叫 mapPartition(new VectorSummarizerPartition(bCov))
    - 3.2.1.1.1)呼叫VectorSummarizerPartition.mapPartition，其遍歷列表，列表中的每一個變數 sv 是 x-vec。srt = srt.visit(sv)，會根據每一個新輸入重新計算count，sum，squareSum，normL1..，這樣就得到了本partiton中輸入每列的這些統計數值。
  - 3.2.1.2)呼叫 reduce(VectorSummarizerUtil.merge(value1, value2)) 來歸併每一個partition的結果。
- 3.2.2)呼叫map(BaseVectorSummarizer summarizer)，其實呼叫到DenseVectorSummarizer，就是生成一個DenseVectorSummary向量，裡面是count，sum，squareSum，normL1，min，max，numNonZero。
3.3)呼叫 coefficientDim = summary.map
3.4)呼叫 meanVar = coefficientDim.map，最後得到 Tuple2.of(coefficientDim, meanVar)

private Tuple2<DataSet<Integer>, DataSet<DenseVector[]>> getStatInfo(
    DataSet<Tuple3<Double, Double, Vector>> trainData, final boolean standardization) {
    if (standardization) {
        DataSet<BaseVectorSummary> summary = StatisticsHelper.summary(trainData.map(
            new MapFunction<Tuple3<Double, Double, Vector>, Vector>() {
                @Override
                public Vector map(Tuple3<Double, Double, Vector> value) throws Exception {
                    return value.f2; //獲取訓練資料中的 x-vec
                }
            }).withForwardedFields());

        DataSet<Integer> coefficientDim = summary.map(new MapFunction<BaseVectorSummary, Integer>() {
            public Integer map(BaseVectorSummary value) throws Exception {
                return value.vectorSize(); // 獲取dimension
            }
        });

        DataSet<DenseVector[]> meanVar = summary.map(new MapFunction<BaseVectorSummary, DenseVector[]>() {
            public DenseVector[] map(BaseVectorSummary value) {
                if (value instanceof SparseVectorSummary) {
                    // 計算min, max
                    DenseVector max = ((SparseVector)value.max()).toDenseVector();
                    DenseVector min = ((SparseVector)value.min()).toDenseVector();
                    for (int i = 0; i < max.size(); ++i) {
                        max.set(i, Math.max(Math.abs(max.get(i)), Math.abs(min.get(i))));
                        min.set(i, 0.0);
                    }
                    return new DenseVector[] {min, max};
                } else {
                    // 計算standardDeviation
                    return new DenseVector[] {(DenseVector)value.mean(),
                        (DenseVector)value.standardDeviation()};
                }
            }
        });
        return Tuple2.of(coefficientDim, meanVar);
    } 
}

5.4 對輸入資料做標準化和插值

這裡對應基本流程的 4) 。

對輸入資料做標準化和插值 standardization and interception。上面得到的 meanVar 作為引數傳入。這裡是對 x-vec 做標準化。

比如原始輸入Row是"(1.0,16.8,1.0 7.0 9.0)"，其中 x-vec 是"1.0 7.0 9.0"，進行標準化之後，x-vec 變成了 4 項，第一項是固定值 "1.0 "， 4 項是 "1.0 1.0 1.4657097546055162 1.4770978917519928"，所以轉換後的Row是"(1.0,16.8,1.0 1.0 1.4657097546055162 1.4770978917519928)"。

為什麼第一項是固定值 "1.0 " ？因為按照線性模型 f(x)=w^Tx+b，我們應該得出一個常數 b，這裡設定 "1.0 "，就是 b 的初始值。

private DataSet<Tuple3<Double, Double, Vector>> preProcess(
    return initData.map(
        new RichMapFunction<Tuple3<Double, Double, Vector>, Tuple3<Double, Double, Vector>>() {
            private DenseVector[] meanVar;

            @Override
            public Tuple3<Double, Double, Vector> map(Tuple3<Double, Double, Vector> value){
// value = {Tuple3@9791} "(1.0,16.8,1.0 7.0 9.0)"
                Vector aVector = value.f2;
// aVector = {DenseVector@9792} "1.0 7.0 9.0"
                if (aVector instanceof DenseVector) {
                    DenseVector bVector;
                    if (standardization) {
                        if (hasInterceptItem) {
                            bVector = new DenseVector(aVector.size() + 1);
                            bVector.set(0, 1.0); // 設定了固定值
                            for (int i = 0; i < aVector.size(); ++i) {
                                // 對輸入資料做標準化和插值
                                bVector.set(i + 1, (aVector.get(i) - meanVar[0].get(i)) / meanVar[1].get(i));
                            }
                        } 
                    } 
// bVector = {DenseVector@9814} "1.0 1.0 1.4657097546055162 1.4770978917519928"
                    return Tuple3.of(value.f0, value.f1, bVector);
                } 
            }
        }).withBroadcastSet(meanVar, MEAN_VAR);
}

// 這裡是對 x-vec 做標準化。比如原始輸入Row是"(1.0,16.8,1.0 7.0 9.0)"，其中 x-vec 是"1.0 7.0 9.0"，進行標準化之後，x-vec 變成了 4 項，第一項是 "1.0 "，是 "1.0 1.0 1.4657097546055162 1.4770978917519928"，所以轉換後的Row是"(1.0,16.8,1.0 1.0 1.4657097546055162 1.4770978917519928)"

至此，輸入處理完畢。

比如原始輸入Row是"(1.0,16.8,1.0 7.0 9.0)"，其中 x-vec 是"1.0 7.0 9.0"。

進行標準化之後，x-vec 變成了 4 項：{ 第1項是固定值 "1.0 "，所以4 項是 "1.0 1.0 1.4657097546055162 1.4770978917519928" }，

轉換後的Row是"(1.0,16.8,1.0 1.0 1.4657097546055162 1.4770978917519928)"。即weight 是1.0，y-value是16.8，後續4個是x-vec。

下面我們可以開始進行優化模型了，敬請期待下文。