1. 使用

spark給出的example中涉及到LBFGS有兩個，分別是LBFGSExample.scala和LogisticRegressionWithLBFGSExample.scala，第一個是直接使用LBFGS直接訓練，需要指定一系列優化引數，優點是比較靈活，可以自己控制的引數較多。後者使用了LogisticRegressionWithLBFGS，只能設定class的個數，其他引數都是固定的，其實就是將第一個中自己能控制的引數，都指定了預設值，適合剛開始時學習。
下面會以第二個為例，因為其中封裝了第一個。

// Split data into training (60%) and test (40%)
 
.
val splits = data.randomSplit(Array(0.6, 0.4), seed = 11L)
val training = splits(0).cache()
val test = splits(1)
// Run training algorithm to build the model
val model = new LogisticRegressionWithLBFGS()
      .setNumClasses(10)
      .run(training)

設定包含10分類，調節訓練和驗證資料集的比例

2. 邏輯迴歸與L-BFGS演算法

這一節簡單介紹邏輯迴歸和L-BFGS演算法原理，以便於程式碼實現相互對照，但是不會做嚴格的數學推理。

2.1. softmax

這裡貼兩頁程式碼註釋中給出的PPT（https://www.slideshare.net/dbtsai/2014-0620-mlor-36132297）


簡單解釋下PPT中的內容，第二頁似然函式最後的形式中，α只有當樣本的label是0的時候取1，第一項為0，其他的label取0，也就是第一項是有值的。第一項雖然在累加號裡面，但是我們注意權重w的下標yk，這意味著當累加索引k從0到N迴圈時，只有其與樣本label y（label也是從1到N）相等時，x*w才會被計算，在程式碼中這項是marginY。
對於第二項，loss和gradient計算中含有指數運算，這部分累加在程式碼中是margins變數，如果資料中存在異常點，對應到指數如果超過了709.78，就會溢位，導致訓練失敗，這裡做了一點trick，在涉及到指數計算的地方，會先判斷計算出的指數部分是否大於0，如果大於0，直接在後面加指數部分，相當於在log裡面再除掉，這樣就不會溢位了。

log(1+∑i=1Nexp(xwi))=log(exp(−m)(1+∑i=1Nexp(xwi)))+m=log(1+exp(−m)+∑i=1Nexp(xwi−m)−1)+m
對比第一個式子和第三個式子，我們看到故意加1又在後面減1，這裡是為了保持形式上的一致性，程式碼實現的時候，當margin大於0，後面的累加是sum變數，根據margin是否大於0，分別計算，其他部分一樣計算。並且當margin大於0，會計算最大maxMargin，累加過程中，累加到max的時候，因為指數部分都會先減掉maxMargin，相當於結果是1，也就和後面的-1合併掉了，因此在程式碼計算的時候，當累加到maxMargin的時候，直接加上了前面那個-maxMargin了。
梯度的計算也是類似的。

2.2. L-BFGS

關於L-BFGS之前也貼過轉載的文章，詳細介紹了從牛頓法到L-BFGS是怎麼發展過來的，這裡僅僅貼出演算法部分，同時結合spark的實現簡單說下

這裡的演算法僅僅是計算搜尋方向，delta和L的作用主要是當k較小時（小於m），就迭代k輪，較大時（大於m）就迭代m次。
對比Wikipedia（https://en.wikipedia.org/wiki/Limited-memory_BFGS）， 這裡每輪的初值的設定應該是

D0k=sTk−1yk−1/yTk−1yk−1

3. 資料結構

3.1. LogisticRegressionWithLBFGS

上文有提到，這裡主要是封裝了LBFGS的一些引數，使得對外暴露的引數很少，值得注意的是這裡的正則化引數regParam預設是0，也就是不使用正則化。程式碼部分比較簡單，這裡簡單說下。程式碼中首先設定對特徵進行縮放，然後設定LBFGS的剃度使用LogistcGradient，正則化使用SquaredL2Updater，根據是2/多分類設定對label進行驗證的DataValidators，這裡值得注意的是，在run介面中，如果是二分類（numOfLinearPredictor == 1），這裡其實使用的是ml包中的ElasticNet，只有多分類時才是真正使用的LBFGS，如果想在二分類時使用，需要自己參考LBFGSExample.scala呼叫。
值得注意的是

def setNumClasses(numClasses: Int): this.type = { 
  require(numClasses > 1)
  numOfLinearPredictor = numClasses - 1 
  if (numClasses > 2) { 
    optimizer.setGradient(new LogisticGradient(numClasses))
  } 
  this
} 

可以看到這裡同時設定了numOfLinearPredictor，也就是預測器的個數，比class個數少一個，這種多分類方法是有pivot class的，也就是優化時只需要優化k-1個係數，預設class 0的margin是0，對應演算法中P(y=0)的概率分子上是1，沒有係數向量。可以看到後面model的predict函式的做法，同時在後面初始化係數向量時，也是用了numOfLinearPredictor，可以注意下。

3.2. LogisticRegressionModel

都是設定其中的諸如weight，threshold等屬性的，比較簡單，這裡介紹predict函式，看是怎樣分類的，以便其他環境中自己實現分類函式

    override protected def predictPoint(
    //樣本特徵向量
      dataMatrix: Vector,  
      //weight
      weightMatrix: Vector,  
      //截距
      intercept: Double) = {  
    require(dataMatrix.size == numFeatures)

    // If dataMatrix and weightMatrix have the same dimension, it's binary logistic regression.
    if (numClasses == 2) {
      val margin = dot(weightMatrix, dataMatrix) + intercept
      val score = 1.0 / (1.0 + math.exp(-margin))
      threshold match {
        //分類
        case Some(t) => if (score > t) 1.0 else 0.0
        //迴歸
        case None => score
      }
    } else {
      /**
       * Compute and find the one with maximum margins. If the maxMargin is negative, then the
       * prediction result will be the first class.
       *
       * PS, if you want to compute the probabilities for each outcome instead of the outcome
       * with maximum probability, remember to subtract the maxMargin from margins if maxMargin
       * is positive to prevent overflow.
       */
      var bestClass = 0
      var maxMargin = 0.0
      val withBias = dataMatrix.size + 1 == dataWithBiasSize
      (0 until numClasses - 1).foreach { i =>
        var margin = 0.0
        dataMatrix.foreachActive { (index, value) =>
          if (value != 0.0) margin += value * weightsArray((i * dataWithBiasSize) + index)
        }
        // Intercept is required to be added into margin.
        if (withBias) {
          margin += weightsArray((i * dataWithBiasSize) + dataMatrix.size)
        }
        if (margin > maxMargin) {
          maxMargin = margin
          bestClass = i + 1
        }
      }
      bestClass.toDouble
    }
  }

從程式碼可以看到，多分類預測的過程就是分別將樣本的特徵向量於訓練得到的各個分類的weight相乘，取最大的值對應的label。需要注意的是weights向量裡連續存著1到numClasses-1分類的weight，從margin的計算就可以看出來。

3.3. LogisticGradient

邏輯迴歸的演算法前面有簡要介紹，這裡主要是要計算其剃度和loss，程式碼這裡定義
當有margin大於0，定義

sum=exp(−m)+∑i=1Nexp(xwi−m)−1
同時後面需要加m。所有margin都小於0，則
sum=∑i=1Nexp(xwi)
同時後面也不需要加m

/**
         * For Multinomial Logistic Regression.
         */
        val weightsArray = weights match {
          case dv: DenseVector => dv.values
          case _ =>
            throw new IllegalArgumentException(
              s"weights only supports dense vector but got type ${weights.getClass}.")
        }
        val cumGradientArray = cumGradient match {
          case dv: DenseVector => dv.values
          case _ =>
            throw new IllegalArgumentException(
              s"cumGradient only supports dense vector but got type ${cumGradient.getClass}.")
        }

        // marginY is margins(label - 1) in the formula.
        var marginY = 0.0
        var maxMargin = Double.NegativeInfinity
        var maxMarginIndex = 0

        val margins = Array.tabulate(numClasses - 1) { i =>
          var margin = 0.0
          data.foreachActive { (index, value) =>
            if (value != 0.0) margin += value * weightsArray((i * dataSize) + index)
          }
          //第一項
          if (i == label.toInt - 1) marginY = margin
          if (margin > maxMargin) {
            maxMargin = margin
            maxMarginIndex = i
          }
          margin
        }

        /**
         * When maxMargin > 0, the original formula will cause overflow as we discuss
         * in the previous comment.
         * We address this by subtracting maxMargin from all the margins, so it's guaranteed
         * that all of the new margins will be smaller than zero to prevent arithmetic overflow.
         */
        val sum = {
          var temp = 0.0
          if (maxMargin > 0) {
          //大於0時候的處理
            for (i <- 0 until numClasses - 1) {
              margins(i) -= maxMargin
              if (i == maxMarginIndex) {
                temp += math.exp(-maxMargin)
              } else {
                temp += math.exp(margins(i))
              }
            }
          } else {
          //不會溢位，不需要處理
            for (i <- 0 until numClasses - 1) {
              temp += math.exp(margins(i))
            }
          }
          temp
        }

        for (i <- 0 until numClasses - 1) {
          val multiplier = math.exp(margins(i)) / (sum + 1.0) - {
            if (label != 0.0 && label == i + 1) 1.0 else 0.0
          }
          data.foreachActive { (index, value) =>
            if (value != 0.0) cumGradientArray(i * dataSize + index) += multiplier * value
          }
        }
        //非0的時候，alpha=0，第一項有值
        val loss = if (label > 0.0) math.log1p(sum) - marginY else math.log1p(sum)

        if (maxMargin > 0) {
        //大於0，最後要加上m
          loss + maxMargin
        } else {
          loss
        }

3.4. SquaredL2Updater

L2正則化，我們之前介紹過，這裡不再贅述，值得注意的是，正則化時需要更新優化物件weight，compute函式第一項正是更新過的weight。

4. 訓練

訓練主要在breeze包中的LBFGS.scala和其父類FirstOrderMinimizer的一些函式中

4.1. weight初始化

/**
   * Generate the initial weights when the user does not supply them
   */
  protected def generateInitialWeights(input: RDD[LabeledPoint]): Vector = {
    if (numFeatures < 0) {
      numFeatures = input.map(_.features.size).first()
    }

    /**
     * When `numOfLinearPredictor > 1`, the intercepts are encapsulated into weights,
     * so the `weights` will include the intercepts. When `numOfLinearPredictor == 1`,
     * the intercept will be stored as separated value in `GeneralizedLinearModel`.
     * This will result in different behaviors since when `numOfLinearPredictor == 1`,
     * users have no way to set the initial intercept, while in the other case, users
     * can set the intercepts as part of weights.
     *
     * TODO: See if we can deprecate `intercept` in `GeneralizedLinearModel`, and always
     * have the intercept as part of weights to have consistent design.
     */
    if (numOfLinearPredictor == 1) {
      Vectors.zeros(numFeatures)
    } else if (addIntercept) {
      Vectors.zeros((numFeatures + 1) * numOfLinearPredictor)
    } else {
      Vectors.zeros(numFeatures * numOfLinearPredictor)
    }
  }

根據二/多分類，是否有截距，初始化weight向量，多分類時，多個weight向量是連續存放在一個vector結構裡的，有過有截距，相當於每個weight多加了一維。

4.2. 樣本資料check

主要是檢查樣本label的範圍

/**
   * Function to check if labels used for k class multi-label classification are
   * in the range of {0, 1, ..., k - 1}.
   *
   * @return True if labels are all in the range of {0, 1, ..., k-1}, false otherwise.
   */
  @Since("1.3.0")
  def multiLabelValidator(k: Int): RDD[LabeledPoint] => Boolean = { data =>
    val numInvalid = data.filter(x =>
      x.label - x.label.toInt != 0.0 || x.label < 0 || x.label > k - 1).count()
    if (numInvalid != 0) {
      logError("Classification labels should be in {0 to " + (k - 1) + "}. " +
        "Found " + numInvalid + " invalid labels")
    }
    numInvalid == 0
  }

4.3. scaling

對樣本特徵值的範圍進行縮放，加快收斂速度，但是這裡僅僅對特徵值除以標準差，沒有做真正的歸一化。

val scaler = if (useFeatureScaling) {
      new StandardScaler(withStd = true, withMean = false).fit(input.map(_.features))
    } else {
      null
    }

// Prepend an extra variable consisting of all 1.0's for the intercept.
// TODO: Apply feature scaling to the weight vector instead of input data.
val data =
    if (addIntercept) {
        if (useFeatureScaling) {
          input.map(lp => (lp.label, appendBias(scaler.transform(lp.features)))).cache()
        } else {
          input.map(lp => (lp.label, appendBias(lp.features))).cache()
        }
    } else {
        if (useFeatureScaling) {
          input.map(lp => (lp.label, scaler.transform(lp.features))).cache()
        } else {
          input.map(lp => (lp.label, lp.features))
        }
    }

這段程式碼同時對有截距的情況進行了處理，因為我們之前把截距當成額外的一維放在了weight向量裡，這裡對樣本特徵也對應增加以為，且特徵值全部設定成1，就可以把截距也當成特徵統一處理了。