Spark MLlib之線性迴歸源碼分析

1.理論基礎

    線性迴歸(Linear Regression)問題屬於監督學習（Supervised Learning）範疇，又稱分類（Classification）或歸納學習(Inductive Learning);這類分析中訓練資料集中給出的資料類標是確定的；機器學習的目標是，對於給定的一個訓練資料集，通過不斷的分析和學習產生一個聯絡屬性集合和類標集合的分類函數(Classification Function)或預測函數(Prediction Function),這個函數稱為分類模型(Classification Model)或預測模型(Prediction Model)；通過學習得到的模型可以是一個決策樹,規格集,貝葉斯模型或一個超平面。通過這個模型可以對輸入對象的特徵向量預測或對對象的類標進行分類。

   迴歸問題中通常使用最小二乘(Least Squares)法來迭代最優的特徵中每個屬性的比重，通過損失函數(Loss Function)或錯誤函數(Error Function)定義來設定收斂狀態，即作為剃度下降演算法的逼近參數因子。

2.矩陣向量運算庫jblas介紹

    由於spark MLlib中使用jlbas的線性代數運算庫，因此學習和掌握jlbas庫中基本的運算，對分析和學習spark中MLlib很多演算法很有協助；下面使用jlbas中DoubleMatrix矩陣對jlbas中基本運算進行簡單介紹：

    val matrix1 = DoubleMatrix.ones(10, 1) //建立所有值為1的10*1矩陣    val matrix2 = DoubleMatrix.zeros(10, 1) //建立所有值為0的10*1矩陣    matrix1.put(1, -10)    val absSum = matrix1.norm1() //絕對值之和    val euclideanNorm = matrix1.norm2() //歐幾裡德距離    val matrix3 = (matrix1.addi(matrix2))    val matrix4 = new DoubleMatrix(1, 10, (1 to 10).map(_.toDouble): _*) //建立Double向量對象    println("print init value:matrix3=" + matrix3)    println("print init value:matrix4=" + matrix4)    println("matrix sub matrix:" + matrix3.sub(matrix4) + "," + matrix4.sub(10)) //減法運算    println("matrix add matrix:" + matrix3.add(matrix4) + "," + matrix4.add(10)) //加法運算    println("matrix mul matrix:" + matrix3.mul(matrix4) + "," + matrix4.mul(10)) //乘法運算    println("matrix div matrix:" + matrix3.div(matrix4) + "," + matrix4.div(10)) //除法運算    println("matrix dot matrix:" + matrix3.dot(matrix4)) //向量積    val matrix5 = DoubleMatrix.ones(10, 10)    println("N*M Vector Matrix sub OP:\n" + matrix5.subRowVector(matrix4) + "\n" + matrix5.subColumnVector(matrix4)) //多個物件減法運算    println("N*M Vector Matrix add OP:\n" + matrix5.addRowVector(matrix4) + "\n" + matrix5.addColumnVector(matrix4)) //多個物件加法運算    println("N*M Vector Matrix mul OP:\n" + matrix5.mulRowVector(matrix4) + "\n" + matrix5.mulColumnVector(matrix4)) //多個物件乘法運算    println("N*M Vector Matrix div OP:\n" + matrix5.divRowVector(matrix4) + "\n" + matrix5.divColumnVector(matrix4)) //多個物件除法運算

3.梯度下降(Gradient Descent)演算法介紹

    梯度下降演算法用於在迭代過程中逐漸降階，不斷更新特徵權重向量，從而得到無限接近或擬合的最優特徵權重向量 ；梯度下降演算法主要有兩種,第一種是批量梯度下降(Batch Gradient Descent)演算法，此種方式實現過程是對權重向量進行累加，然後批次更新的一種方式，一般不實用於大規模資料集處理；另外一種是隨機梯度下降(Stochastic Gradient Descent)演算法，這種方式對給定訓練資料集中每個對象都進行權重計算和更新，在某些情況下容易收斂到局部最優解上。Spark MLlib庫中主要使用隨機梯度下降演算法。為了更好的理解MLlib庫中隨機梯度演算法(MLlib庫中類尾碼名以SGD結尾的所有演算法)實現，下面是使用隨機梯度演算法進行直線擬合的Demo：

  def sgdDemo {    val featuresMatrix: List[List[Double]] = List(List(1, 4), List(2, 5), List(5, 1), List(4, 2))//特徵矩陣    val labelMatrix: List[Double] = List(19, 26, 19, 20)//真實值向量    var theta: List[Double] = List(0, 0)    var loss: Double = 10.0    for {      i <- 0 until 1000 //迭代次數      if (loss > 0.01) //收斂條件loss<=0.01    } {      var error_sum = 0.0 //總誤差      var j = i % 4      var h = 0.0      for (k <- 0 until 2) {        h += featuresMatrix(j)(k) * theta(k)      } //計算給出的測試資料集中第j個對象的計算類標籤      error_sum = labelMatrix(j) - h //計算給出的測試資料集中類標籤與計算的類標籤的誤差值      var cacheTheta: List[Double] = List()      for (k <- 0 until 2) {        val updaterTheta = theta(k) + 0.001 * (error_sum) * featuresMatrix(j)(k)        cacheTheta = updaterTheta +: cacheTheta      } //更新權重向量      cacheTheta.foreach(t => print(t + ","))      print(error_sum + "\n")      theta = cacheTheta      //更新誤差率      var currentLoss: Double = 0      for (j <- 0 until 4) {        var sum = 0.0        for (k <- 0 until 2) {          sum += featuresMatrix(j)(k) * theta(k)        }        currentLoss += (sum - labelMatrix(j)) * (sum - labelMatrix(j))      }      loss = currentLoss      println("loss->>>>" + loss / 4 + ",i>>>>>" + i)    }  }

4.MLlib線性迴歸源碼分析

     MLlib中可用的線性迴歸演算法有：LinearRegressionWithSGD，RidgeRegressionWithSGD，LassoWithSGD；MLlib迴歸分析中涉及到的主要類有，GeneralizedLinearAlgorithm，GradientDescent。下面以對LinearRegressionWithSGD實現進行主要分析。

    第一步：在使用LinearRegressionWithSGD演算法之前首先將輸入資料解析成包含類標和特徵向量的LabeledPoint對象的RDD彈性分布式資料集合。

    第二步：調用LinearRegressionWithSGD伴生對象中的train方法傳輸第一步建立的RDD集合和最大迭代次數，在train中主要實現建立一個新的LinearRegressionWithSGD對象，初始化梯度下降算使用使用最小二乘梯度下降演算法SquaredGradient以及更新權重向量使用SimpleUpdater，執行父類GeneralizedLinearAlgorithm中的run方法進行權重向量和攔截參數計算，並返回訓練得到的模型屬性權重向量

    LinearRegressionWithSGD伴生對象中train方法實現

  def train(      input: RDD[LabeledPoint],      numIterations: Int,      stepSize: Double,//預設步長為1      miniBatchFraction: Double)//每次跌打使用的batch因子，預設為1    : LinearRegressionModel =  {    new LinearRegressionWithSGD(stepSize, numIterations, miniBatchFraction).run(input)  }

   LinearRegressionWithSGD中run方法實現

  def run(input: RDD[LabeledPoint], initialWeights: Array[Double]) : M = {    // Check the data properties before running the optimizer    if (validateData && !validators.forall(func => func(input))) {//預驗證輸入資料的合法性，validators中儲存驗證的所有方法      throw new SparkException("Input validation failed.")    }    // Prepend an extra variable consisting of all 1.0's for the intercept.    val data = if (addIntercept) {//判讀是否需要添加intercept      input.map(labeledPoint => (labeledPoint.label, 1.0 +: labeledPoint.features))    } else {      input.map(labeledPoint => (labeledPoint.label, labeledPoint.features))    }//對象轉化為元組(類標籤,特徵)    val initialWeightsWithIntercept = if (addIntercept) {      0.0 +: initialWeights    } else {      initialWeights    }//初始化權重特徵    val weightsWithIntercept = optimizer.optimize(data, initialWeightsWithIntercept)//返回最佳化後的權重    val (intercept, weights) = if (addIntercept) {      (weightsWithIntercept(0), weightsWithIntercept.tail)    } else {      (0.0, weightsWithIntercept)    }    logInfo("Final weights " + weights.mkString(","))    logInfo("Final intercept " + intercept)    createModel(weights, intercept)//使用計算後的權重向量和截距建立模型  }

其中optimizer.optimize(data, initialWeightsWithIntercept)是LinearRegressionWithSGD實現的核心，oprimizer的類型為GradientDescent，optimize方法中主要調用GradientDescent伴生對象的runMiniBatchSGD方法，返回當前迭代產生的最優特徵權重向量。

GradientDescentd對象中optimize方法實現

  def optimize(data: RDD[(Double, Array[Double])], initialWeights: Array[Double])    : Array[Double] = {//返回最佳化後的權重向量，和迭代過程中誤差    val (weights, stochasticLossHistory) = GradientDescent.runMiniBatchSGD(        data,        gradient,        updater,        stepSize,        numIterations,        regParam,        miniBatchFraction,        initialWeights)    weights  }

GradientDescent伴生對象中runMiniBatchSGD方法實現

  def runMiniBatchSGD(    data: RDD[(Double, Array[Double])],    gradient: Gradient,//SquaredGradient—平方剃度下降演算法    updater: Updater,//SimpleUpdater    stepSize: Double,//1.0    numIterations: Int,//100    regParam: Double,//0.0    miniBatchFraction: Double,//1.0    initialWeights: Array[Double]) : (Array[Double], Array[Double]) = {    val stochasticLossHistory = new ArrayBuffer[Double](numIterations)    val nexamples: Long = data.count()    val miniBatchSize = nexamples * miniBatchFraction    // Initialize weights as a column vector／／建立一維向量，第一個參數為行數，第二個參數為列數，第三個參數開始為值    var weights = new DoubleMatrix(initialWeights.length, 1, initialWeights:_*)    var regVal = 0.0        for (i <- 1 to numIterations) {      /**       * 使用平方梯度下降演算法       * gradientSum:本次選擇迭代樣本的梯度總和，       * lossSum:本次選擇迭代樣本的誤差總和       */      val (gradientSum, lossSum) = data.sample(false, miniBatchFraction, 42+i).map {        case (y, features) =>//(標籤，特徵)          val featuresCol = new DoubleMatrix(features.length, 1, features:_*)          val (grad, loss) = gradient.compute(featuresCol, y, weights)//(特徵，標籤，特徵屬性權重向量)          /**           * class SquaredGradient extends Gradient {  override def compute(data: DoubleMatrix, label: Double, weights: DoubleMatrix):       (DoubleMatrix, Double) = {    val diff: Double = data.dot(weights) - label//計算當前計算對象的類標籤與實際類標籤值之差    val loss = 0.5 * diff * diff//當前平方梯度下降值    val gradient = data.mul(diff)    (gradient, loss)  }}           */          (grad, loss)//當前訓練對象的特徵權重向量和誤差      }.reduce((a, b) => (a._1.addi(b._1), a._2 + b._2)) //計算本次迭代所選訓練資料新的特徵權重向量之和和誤差總和           /**       * NOTE(Xinghao): lossSum is computed using the weights from the previous iteration       * and regVal is the regularization value computed in the previous iteration as well.       */      stochasticLossHistory.append(lossSum / miniBatchSize + regVal)//miniBatchSize=樣本中對象數量＊batch因子，regVal為迴歸因子      val update = updater.compute(        weights, gradientSum.div(miniBatchSize), stepSize, i, regParam)//weights：屬性向量中設定的權重因子，regParam：為迴歸參數，stepSize：計算步長，i：當前迭代次數      /**       * class SimpleUpdater extends Updater {       * /**          * weihtsOld:上一次迭代計算後的特徵權重向量          * gradient:本次迭代計算的特徵權重向量          * stepSize:迭代步長          * iter:當前迭代次數          * regParam:迴歸參數           */  override def compute(weightsOld: DoubleMatrix, gradient: DoubleMatrix,      stepSize: Double, iter: Int, regParam: Double): (DoubleMatrix, Double) = {    val thisIterStepSize = stepSize / math.sqrt(iter)//以當前迭代次數的平方根的倒數作為本次迭代趨近(下降)的因子    val normGradient = gradient.mul(thisIterStepSize)    (weightsOld.sub(normGradient), 0)//返回本次剃度下降後更新的特徵權重向量  }}   *       */      weights = update._1      regVal = update._2//使用SimpleUpdater值為0    }    logInfo("GradientDescent finished. Last 10 stochastic losses %s".format(      stochasticLossHistory.takeRight(10).mkString(", ")))    (weights.toArray, stochasticLossHistory.toArray)  }

  在MiniBatchSGD中主要實現對輸入資料集進行迭代抽樣，通過使用SquaredGradient作為梯度下降演算法，使用SimpleUpdater作為更新演算法，不斷對抽樣資料集進行迭代計算從而找出最優的特徵權重向量解。

  使用官方測試代碼如下：

  def linearRegressionAPITest(sc: SparkContext) {    val url = "/Users/yangguo/hadoop/spark/mllib/data/ridge-data/lpsa.data"    val data = sc.textFile(url)    val parseData = data.map { line =>      val parts = line.split(',')      LabeledPoint(parts(0).toDouble, parts(1).split(' ').map(x => x.toDouble).toArray)    }    val numIterations = 20    val model = LinearRegressionWithSGD.train(parseData, numIterations)    val valuesAndPreds = parseData.map { point =>      val prediction = model.predict(point.features)      (point, prediction)    }    valuesAndPreds.foreach { case (v, p) => print("[" + v.label + "," + p + "]"); v.features.foreach(base => print(base + "--")); println("\n") }    val isSuccessed = valuesAndPreds.map { case (v, p) => math.pow((p - v.label), 2) }.reduce(_ + _) / valuesAndPreds.count    println(isSuccessed)  }

參考:

[1].http://rdc.taobao.org/?p=2163

[2].http://cs229.stanford.edu/notes/cs229-notes1.pdf

[3].http://blog.sina.com.cn/s/blog_62339a2401015jyq.html

[4].http://blog.csdn.net/pennyliang/article/details/6998517

[5].http://en.wikipedia.org/wiki/Lasso_(statistics)#Lasso_method