Spark MLlib 源代码解析之逻辑回归LogisticRegression

最新推荐文章于 2020-05-05 00:59:11 发布

stevekangpei

最新推荐文章于 2020-05-05 00:59:11 发布

阅读量4.9k

点赞数 7

分类专栏： MLlib源代码解读文章标签： spark 源代码逻辑回归

本文链接：https://blog.csdn.net/stevekangpei/article/details/75331960

版权

MLlib源代码解读专栏收录该内容

6 篇文章 3 订阅

订阅专栏

Spark MLlib 逻辑回归 LogisticRegression

模型

逻辑回归的本质是线性回归，只是在特征到结果的过程上加上了一层映射。即首先需要把特征进行求和，然后将求和后的结果应用于一个g(z)函数。g(z)可以将值映射到0或者是1上面。

逻辑回归和多重线性回归有很多的相似之处。最大的区别是他们的因变量不同。这两个回归也可以统一归为广义线性模型。在spark mllib实现的过程中也是先定义好父类广义线性模型，然后让线性回归和逻辑回归去继承这个类，重新覆盖里面的一些参数，比如说GradientDescent，Updater等。

逻辑回归可以是二分类的，也可以是多分类的。本篇博客主要讨论的是二分类的模型。二分类也是最常用的逻辑回归模型。

逻辑回归主要是用于分类。

比如根据某些特点，来判断这个病人是不是得了某种疾病。
预测，根据模型，预测在不同的自变量情况下，发生某病或某种情况的概率有多大；

逻辑回归和线性回归的不同点在于，将线性回归的输出范围压缩到了0和1之间。
逻辑函数又称Sigmoid函数，函数形式如下所示：

g (z) = 1 1 + e - z

$g(z)=\frac{1}{1+e^{-z}}$
这个函数的图形大概如下所示：

Sigmoid

公式推导

假设样本有n个特征x=(x1,x2…xn), 设p(y=1|x)为观测样本y相对于事件x发生的概率，用Sigmoid函数表示为：

p (y = 1 | x) = 1 1 + e - g ( x )

$p(y=1|x)=\frac{1}{1+e^{-g(x)}}$

其中，g(x)=w0+w1x1+…..+wnxn.
在x条件下y不发生的概率为

p (y = 0 | x) = 1 - p (y = 1 | x) = 1 1 + e g ( x )

$p(y=0|x)=1-p(y=1|x)=\frac{1}{1+e^{g(x)}}$

如果目前有m个相互独立的事件，y=( $y^{(1)}, y^{(2)}, y^{(3)}, y^{(4)}, y^{(n)}$ ), 则一个事件y(i)发生的：

p (y (i)) = p y (i) * (1 - p) 1 - y (i)

$p(y^{(i)})=p^{y^{(i)}}*(1-p)^{1-y^{(i)}}$

当y=1的时候，后面的一项便消失了。当y=0的时候，前面这项便消失了。
对于整个样本集，即m个独立样本出现的似然函数(因为每个样本独立，所以m个样本出现的概率就是他们各自出现的概率相乘)。
其极大似然函数

L (Θ) = \prod i = 1 m f (x, Θ) = \prod i = 1 m (g (x)) y (i) * (1 - g (x)) 1 - y (i)

$L(\Theta )=\prod_{i=1}^{m}f(x,\Theta )=\prod_{i=1}^{m}(g(x))^{y^{(i)}}*(1-g(x))^{1-y^{(i)}}$

我们的目标是求解出这个似然函数的L(theta)取最大值的theta值。即求解出 $\Theta1,\Theta2, \Theta n$ ,使得 $L(\Theta)$ 取得极大值。

对 $L(\Theta)$ 取对数有：

L (Θ) = l o g (\prod p (y (i) = 1 | x (i)) (y (i)) * (1 - p (y (i) = 1 | x)) 1 - y (i))

$L(\Theta )=log(\prod p(y^{(i)}=1|x^{(i)})^{(y^{(i)})}*(1-p(y^{(i)}=1|x))^{1-y^{(i)}})$

= \sum i = 1 m y (i) l o g (p (y (i) = 1 | x (i))) + (1 - y i) l o g (1 - p (y (i) = 1 | x (i)))

$=\sum_{i=1}^{m}y^{(i)}log(p(y^{(i)}=1|x^{(i)}))+(1-y^{i})log(1-p(y^{(i)}=1|x^{(i)}))$

= \sum i = 1 m y (i) l o g p ( y ( i ) = 1 | x ( i ) ) 1 - p ( y ( i ) = 1 | x ( i ) ) + \sum i = 1 m l o g (1 - p (y (i) = 1 | x (i)))

$=\sum_{i=1}^{m}y^{(i)}log\frac{p(y^{(i)}=1|x^{(i)})}{1-p(y^{(i)}=1|x^{(i)})}+\sum_{i=1}^{m}log(1-p(y^{(i)}=1|x^{(i)}))$

= \sum i = 1 m y (i) (Θ 0 + Θ 1 x 1 + . . . + Θ n x n) + \sum i = 1 m l o g (1 - p (y (i) = 1 | x (i)))

$=\sum_{i=1}^{m}y^{(i)}(\Theta _{0}+\Theta _{1} x_{1}+...+\Theta_{n} x_{n} )+\sum_{i=1}^{m}log(1-p(y^{(i)}=1|x^{(i)}))$

= \sum i = 1 m y (i) (Θ T x i) - \sum i = 1 m l o g (1 + e Θ T x i)

$=\sum_{i=1}^{m}y^{(i)}(\Theta ^{T}x^{i})-\sum_{i=1}^{m}log(1+e^{\Theta ^{T}x^{i}})$

求解最大似然估计其实就是求解这个的最大值时的 $\Theta$ 值，这里可以采用梯度上升法求解。也可以通过乘以一个-1/m，将问题转换为一个梯度下降法来求解。将 $L(\Theta)$ 转换为 $J(\Theta)$ .

J (Θ) = - 1 m * L (Θ)

$J(\Theta )=-\frac{1}{m}*L(\Theta )$

梯度下降算法可以继续沿用之前的线性回归中的广义梯度下降算法，只不过里面的参数有所改变。

梯度下降算法

$\Theta$ 更新过程:

Θ j : = Θ j - α * \partial J ( Θ ) \partial Θ j

$\Theta_{j} := \Theta_{j} -\alpha * \frac{\partial J(\Theta )}{\partial \Theta _{j}}$

\partial J ( Θ ) \partial Θ j = - 1 m \sum i = 1 m (y (i) x (i) j - g (Θ T x (i)) x (i) j)

$\frac{\partial J(\Theta )}{\partial \Theta _{j}} =-\frac{1}{m}\sum_{i=1}^{m}(y^{(i)}x_{j}^{(i)}-g(\Theta ^{T}x^{(i)})x_{j}^{(i)})$

= 1 m \sum i = 1 m (h Θ (x (i) - y (i)) x (i) j)

$=\frac{1}{m}\sum_{i=1}^{m}(h_{\Theta}(x^{(i)}-y^{(i)})x_{j}^{(i)})$

$\Theta$ 更新过程可以写成

Θ j : = Θ j - α 1 m \sum i = 1 m (h Θ (x (i)) - y (i)) x (i) j

$\Theta{j}:=\Theta{j}- \alpha \frac{1}{m}\sum_{i=1}^{m}(h_{\Theta}(x^{(i)})-y^{(i)})x_{j}^{(i)}$

向量化

上述过程是一个求和的过程，需要for循环，比较耗时。可以考虑将这个过程转换为一个向量或者矩阵的相乘，可以省去大量的时间。考虑训练数据的格式如下。x矩阵是一个m*n的矩阵。其每一行表示一条数据，每一列表示一行数据的特征。

x = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ x 1 . . . . . . x m ⎤ ⎦ ⎥ ⎥ ⎥ ⎥

$x=\begin{bmatrix} x1\\ ...\\ ...\\ xm\\ \end{bmatrix}$

x = ⎡ ⎣ ⎢ ⎢ x 11 x m 1 . . . . . . x 1 n x m n ⎤ ⎦ ⎥ ⎥

$x=\begin{bmatrix} x_{11}& ... &x_{1n} \\ & & \\ x_{m1}& ... &x_{mn} \end{bmatrix}$

y = ⎡ ⎣ ⎢ ⎢ y 1 . . . y m ⎤ ⎦ ⎥ ⎥

$y=\begin{bmatrix} y1\\ ...\\ ym\end{bmatrix}$ ,

Θ = ⎡ ⎣ ⎢ ⎢ Θ 1 . . . Θ n ⎤ ⎦ ⎥ ⎥

$\Theta =\begin{bmatrix} \Theta1 \\ ...\\ \Theta n \end{bmatrix}$

所以 A= x* $\Theta$ ,

E = h Θ (x) - Y

$E =h_{\Theta}(x)-Y$

= ⎡ ⎣ ⎢ ⎢ ⎢ g (A 1) - y 1 . . . g (A n) - y n ⎤ ⎦ ⎥ ⎥ ⎥

$=\begin{bmatrix} g(A_{1})-y_{1}\\ ...\\ g(A_{n})-y_{n}\end{bmatrix}$

= ⎡ ⎣ ⎢ ⎢ e 1 . . . e m ⎤ ⎦ ⎥ ⎥

$=\begin{bmatrix} e1\\ ...\\ em\end{bmatrix}$

=g(A)-y.

g(A)的参数是一个列向量， $\Theta$ 更新过程可以改为：

Θ j : = Θ j - α 1 m \sum i = 1 m (h Θ (x (i)) - y (i)) x (i) j

$\Theta{j}:=\Theta{j}- \alpha \frac{1}{m}\sum_{i=1}^{m}(h_{\Theta}(x^{(i)})-y^{(i)})x_{j}^{(i)}$

= Θ j - α 1 m \sum i = 1 m e (i) x (i) j = Θ j - α * 1 m * x T E

$=\Theta{j}- \alpha \frac{1}{m}\sum_{i=1}^{m}e^{(i)}x_{j}^{(i)}=\Theta _{j}-\alpha *\frac{1}{m} * x^{T}E$

所以可以看到， $\Theta$ 更新过程

1) A=x* $\Theta$
2) E=g(A)-y
3) $Θ = Θ - α * x T E$ $\Theta =\Theta -\alpha *x^{T}E$

正则化

正则化主要是为了解决过拟合问题，防止高次项因子所带来的影响。比如说增大高次项因子的常数项的系数，进而来减小高次项所带来的影响。
关于正则化的具体的细节可以参考一下的链接 :

源码分析:

MLlib中的逻辑回归支持随机梯度下降和拟牛顿法下降算法来实现最优化。本篇博客先讨论随机梯度下降算法，下一篇博客会讨论拟牛顿法。

MLlib逻辑回归的方程:

h w (x) = 1 1 + e - w T x

$h_{w}(x)=\frac{1}{1+e^{-w^{T}x}}$

逻辑回归的损失函数是:

L (w, x, y) = 1 2 * (h w (x) - y) 2

$L(w,x,y)=\frac{1}{2}*(h_{w}(x)-y)^{2}$

逻辑回归使用L2正则化方法。

R (w) = 1 2 ∥ w ∥ 2

$R(w)=\frac{1}{2}\left \| w \right \|^{2}$

每一个样本的梯度的计算方法为：

margin= - w*x

multiplier= (1/ (1+ exp(margin))-y)

gradient = multiplier * x

每个样本的误差，

l o s s = 1 2 * (y - h w (x)) 2

$loss=\frac{1}{2}*(y-h_{w}(x))^{2}$

权重的更新方法：
weight = weight -alpha*(gradient + regParam * weight)

逻辑回归主要包含以下代码:

1) 首先是伴生对象类，LogisticRegressionWithSGD.(包含有静态train方法)
2) 然后是逻辑回归的主类，class LogisticRegressionWithSGD,这个类继承了GeneralizedLinearAlgorithm类。同时执行了父类的run方法。不过里面的部分参数，比如说梯度下降方法，权重更新方法在LogisticRegressionWithSGD有新的定义。父类包含有optimizer.optimize方法。用于执行梯度下降。权重的优化计算调用的是runMiniBatchWithSGD。梯度的计算调用的是Gradient.compute 方法。
3) 最后有一个逻辑回归模型，LogisticRegressionModel类。其里面也包含有predict方法来进行预测。

首先来看看伴生对象类。LogisticRegressionWithSGD(object)
这个是整个逻辑回归算法的入口。主要包含有train方法。train方法的参数包含

input— 训练样本，格式为RDD[LabeledPoint]，其中LabeledPoint格式为(label, features).
numIterations—迭代次数，默认为100次。
stepSize，每次的迭代步长，默认为1.
miniBatchFraction–每次迭代的时候，参与的样本的比例，默认为100%。
initialWeights–初始化权重。

object LogisticRegressionWithSGD {
  // NOTE(shivaram): We use multiple train methods instead of default arguments to support
  // Java programs.

   @Since("1.0.0")
  def train(         //静态的训练方法
      input: RDD[LabeledPoint], //训练样本，RDD格式为(label, features)，注意标签值仅限于0和1
      numIterations: Int,   //迭代次数
      stepSize: Double, //步长
      miniBatchFraction: Double, //表示每次参与迭代计算的样本的比例
      initialWeights: Vector): LogisticRegressionModel = { //初始化权重

    //静态的train方法里面初始化了一个LogisticRegressionWithSGD，通过run方法来进行计算
    new LogisticRegressionWithSGD(stepSize, numIterations, 0.0, miniBatchFraction)
      .run(input, initialWeights)
  }
}

其他的方法基本都是对train方法的重载。

接下来看看逻辑回归类。这个类继承了GeneralizedLinearAlgorithm广义回归类。该类主要初始化梯度下降的方法，梯度更新方法和优化计算方法，然后调用了父类的run方法来执行。
代码如下:

class LogisticRegressionWithSGD private[mllib] (
    private var stepSize: Double,  //步长
    private var numIterations: Int,  //迭代次数
    private var regParam: Double,  //正则化参数
    private var miniBatchFraction: Double)  //每次参与迭代计算的比例

//同样可以看到这个类也继承自广义线性回归，同时调用的也是广义线性回归里面的run方法。
  extends GeneralizedLinearAlgorithm[LogisticRegressionModel] with Serializable {

/// 这里定义了逻辑回归的梯度下降算法，为LogisticGradient
  private val gradient = new LogisticGradient()
  //这里定义了逻辑回归的更新方法，为SquaredL2Updater，为L2正则化
  private val updater = new SquaredL2Updater()

  @Since("0.8.0")
  //根据梯度下降方法，梯度更新方法，新建梯度优化算法。
  override val optimizer = new GradientDescent(gradient, updater)
    .setStepSize(stepSize)
    .setNumIterations(numIterations)
    .setRegParam(regParam)
    .setMiniBatchFraction(miniBatchFraction)
  override protected val validators = List(DataValidators.binaryLabelValidator)

  /**
   * Construct a LogisticRegression object with default parameters: {stepSize: 1.0,
   * numIterations: 100, regParm: 0.01, miniBatchFraction: 1.0}.
   */
  @Since("0.8.0")
  def this() = this(1.0, 100, 0.01, 1.0)

  override protected[mllib] def createModel(weights: Vector, intercept: Double) = {
    new LogisticRegressionModel(weights, intercept)
  }
}

模型调用的run方法是在广义回归算法类里面。run方法首先做一些初始化处理，比如增加偏置项，初始化权重。然后调用optimizer.optimize方法进行计算。这个代码块与之前的线性回归是类似的。

可以看到，在这个方法里面进行了特征维度的检测，然后是数据是否缓存，是否降维

是否需要增加偏置项，最后初始化权重。然后调用了optimizer的optimize方法来进行计算。
最后调用createModel方法，返回结果。

其中，在optimizer是一个GradientDescent类的对象。所以我们之后可以进一步看这个方法，这个方法也是整个线性回归最核心的方法

/**
   * Run the algorithm with the configured parameters on an input RDD
   * of LabeledPoint entries starting from the initial weights provided.
   *执行run方法
   */
  @Since("1.0.0")
  def run(input: RDD[LabeledPoint], initialWeights: Vector): M = { //样本训练的run方法。

    if (numFeatures < 0) { //特征的维度，如果特征的维度被设置为小于0，则取出第一个特征的特征的维度
      numFeatures = input.map(_.features.size).first()
    }
      //看看输入样本有没有缓存。
    if (input.getStorageLevel == StorageLevel.NONE) {
      logWarning("The input data is not directly cached, which may hurt performance if its"
        + " parent RDDs are also uncached.")
    }

    // Check the data properties before running the optimizer
    //检查数据的属性。
    if (validateData && !validators.forall(func => func(input))) {
      throw new SparkException("Input validation failed.")
    }

    /**
     * Scaling columns to unit variance as a heuristic to reduce the condition number:
     * 数据的降维。
     *在优化过程中，收敛率取决于训练数据的维度。
     *通过降维，改变了收敛速度。
     */
    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.
    //是否需要增加偏置项。即theta0的常数项。
    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))
        }
      }

    /**
     * TODO: For better convergence, in logistic regression, the intercepts should be computed
     * from the prior probability distribution of the outcomes; for linear regression,
     * the intercept should be set as the average of response.
     */
     //初始的权重和偏置项。
    val initialWeightsWithIntercept = if (addIntercept && numOfLinearPredictor == 1) {
      appendBias(initialWeights)
    } else {
      /** If `numOfLinearPredictor > 1`, initialWeights already contains intercepts. */
      initialWeights
    }

     //利用了optimizer的optimize方法进行梯度下降。返回最优权重，调用的是GradientDescent的optimize方法。
     //这一行很重要，其最核心的计算在这个optimize方法里面
    val weightsWithIntercept = optimizer.optimize(data, initialWeightsWithIntercept)

    val intercept = if (addIntercept && numOfLinearPredictor == 1) {
      weightsWithIntercept(weightsWithIntercept.size - 1)
    } else {
      0.0
    }

    var weights = if (addIntercept && numOfLinearPredictor == 1) {
      Vectors.dense(weightsWithIntercept.toArray.slice(0, weightsWithIntercept.size - 1))
    } else {
      weightsWithIntercept
    }

   if (useFeatureScaling) {
      if (numOfLinearPredictor == 1) {
        weights = scaler.transform(weights)
      } else {

        var i = 0
        val n = weights.size / numOfLinearPredictor
        val weightsArray = weights.toArray
        while (i < numOfLinearPredictor) {
          val start = i * n
          val end = (i + 1) * n - { if (addIntercept) 1 else 0 }

          val partialWeightsArray = scaler.transform(
            Vectors.dense(weightsArray.slice(start, end))).toArray

          System.arraycopy(partialWeightsArray, 0, weightsArray, start, partialWeightsArray.size)
          i += 1
        }
        weights = Vectors.dense(weightsArray)
      }
    }

    // Warn at the end of the run as well, for increased visibility.
    if (input.getStorageLevel == StorageLevel.NONE) {
      logWarning("The input data was not directly cached, which may hurt performance if its"
        + " parent RDDs are also uncached.")
    }

    // Unpersist cached data
    if (data.getStorageLevel != StorageLevel.NONE) {
      data.unpersist(false)
    }

    createModel(weights, intercept)
  }

权重优化计算。

梯度下降法求解权重.run方法中调用的是optimizer.optimize方法来进行计算。optimizer的类型为GradientDescent类。所以说optimize方法其实是GradientDescent的optimize方法。optimize方法又调用了runMiniBatchWithSGD。

这个是GradientDescent类的optimize方法，其内部又调用了一个runMiniBatchSGD方法
runMiniBatchSGD返回的结果是权重

 //data为RDD格式，其类型为RDD[(Double,Vector)] 训练的数据，initialWeights, 初始化的权重。
  //其返回值为更新的权重，类型为Vector类型。
  //这个方法里面又再一次调用了GradientDescent.runMiniBatchSGD方法。
  def optimize(data: RDD[(Double, Vector)], initialWeights: Vector): Vector = {
    val (weights, _) = GradientDescent.runMiniBatchSGD(
      data,
      gradient,
      updater,
      stepSize,
      numIterations,
      regParam,
      miniBatchFraction,
      initialWeights,
      convergenceTol)
    weights
  }
}

接下来看这个runMiniBatchSGD方法,这个方法是整个线性回归最核心的方法

它大体的思路是，首先初始化好初始权重参数和历史迭代误差的可变数组，然后在每次迭代的时候，广播这个更新的权重到每个rdd。调用treeAggregate算子，每次对数据进行随机采样(无放回采样)，然后先对每个分区的数据进行计算梯度值和误差值，然后接下来对每个分区的计算好的梯度值和误差值进行累加。最后更新权重值。

def runMiniBatchSGD(
      data: RDD[(Double, Vector)],  //输入样本
      gradient: Gradient,      //梯度函数对象，(用于计算损失函数的梯度的一个单一的例子。)
      updater: Updater,    //梯度更新的函数的对象。
      stepSize: Double,      //步长
      numIterations: Int,    //迭代次数
      regParam: Double,       //正则化参数
      miniBatchFraction: Double, //每次迭代参与计算的样本的比例，默认这个比例是1.0
      initialWeights: Vector,  //初始化权重
      convergenceTol: Double): (Vector, Array[Double]) = { //返回为两个元素
    //第一个元素是一个列矩阵，表示的是每一个特征的权重。第二个元素表示的是迭代的损失值。

    if (miniBatchFraction < 1.0 && convergenceTol > 0.0) {
      logWarning("Testing against a convergenceTol when using miniBatchFraction " +
        "< 1.0 can be unstable because of the stochasticity in sampling.")
    }

       //历史迭代的误差数组。存储的是每次迭代的误差值。
    val stochasticLossHistory = new ArrayBuffer[Double](numIterations)
    // Record previous weight and current one to calculate solution vector difference

    var previousWeights: Option[Vector] = None //之前的权重
    var currentWeights: Option[Vector] = None   //当前的权重

    //训练的样本数量。
    val numExamples = data.count()

    // if no data, return initial weights to avoid NaNs
    //如果数据为空，则返回初始的输入参数，即初始的权重和一个误差数组。因为没有找到数据
    if (numExamples == 0) {
      logWarning("GradientDescent.runMiniBatchSGD returning initial weights, no data found")
      return (initialWeights, stochasticLossHistory.toArray)
    }


    //如果数据量乘以采样比例小于1的话，说明miniBatchFraction设置的太小了。弹出警告需要设置的大一点。
    if (numExamples * miniBatchFraction < 1) {
      logWarning("The miniBatchFraction is too small")
    }


    var weights = Vectors.dense(initialWeights.toArray) //将权重初始化，转换为密集向量。
    val n = weights.size  //表示参数的个数

    /**
     * For the first iteration, the regVal will be initialized as sum of weight squares
     * if it's L2 updater; for L1 updater, the same logic is followed.
     *第一次迭代，正则化值初始化为权重的加权平方和。
     */
    var regVal = updater.compute(
      weights, Vectors.zeros(weights.size), 0, 1, regParam)._2

    //这个参数用于表明是否收敛
    var converged = false // indicates whether converged based on convergenceTol
    var i = 1 //i等于1表明第一次迭代
    //


    //接下来就是真个梯度下降法的核心代码。
    //weights权重的迭代计算
    while (!converged && i <= numIterations) {
      //首先广播权重, 注意在每次迭代的开始的时候都需要广播更新的权重值
      val bcWeights = data.context.broadcast(weights)      
      //聚合的时候利用的是treeAggregate方法进行聚合。聚合后返回值的类型为
     //(gradientSum(表示的是梯度的和),lossSum(表示的是损失和),miniBatchSize(表示的是采样比例)

     //treeAggregate算子的执行逻辑如下：
     //treeAggregate的逻辑和aggregate相似，不过它是采用一种多层树结构的模式进行聚合。
     //和aggregate不一样的另一个区别是它的初始值不会被应用到第二个reduce函数上面去。
     //默认的这个tree的深度是2.
     //举个简单的例子。
     //val z = sc.parallelize(List(1,2,3,4,5,6), 2)
     //z.treeAggregate(0)(math.max(_, _), _ + _)
     //res40: Int = 9
     //注意，这个初始值不会作用到第二个reduce函数。s
     //z.treeAggregate(5)(math.max(_, _), _ + _)
     //res42: Int = 11
     // reduce of partition 0 will be max(5, 1, 2, 3) = 5
     // reduce of partition 1 will be max(4, 5, 6) = 6
     // final reduce across partitions will be 5 + 6 = 11

     //梯度计算采用的是随机梯度下降方法。false表示的是不放回抽样
    //随机抽取样本自己，采样时采用不放回采样。每次采样比例为miniBatchFraction。最后一个参数表示为随机种子，每次的值都不一样。
    //保证每次抽样是随机的
      val (gradientSum, lossSum, miniBatchSize) = data.sample(false, miniBatchFraction, 42 + i)
        .treeAggregate((BDV.zeros[Double](n), 0.0, 0L))( //调用BDV.zeros方法初始化一个长度为n的0向量。
            //初始值为一个长度为n的0向量，初始的误差值设为0，

           //计算每一个样本的梯度，然后对所有的样本进行累加。  
          seqOp = (c, v) => {
            // c: (grad, loss, count), v: (label, features)
            //第一个seqOp函数输入为(c,v)类型，返回的是一个c类型。
            //通过调用gradient.compute方法来计算误差值。这个方法输入参数为features，label，权重值，以及得到的梯度值
            //返回的类型为(梯度值，误差值，计数值，样本数+1)
            //默认调用的是LeastSquaresGradient的compute方法。
            val l = gradient.compute(v._2, v._1, bcWeights.value, Vectors.fromBreeze(c._1))
            (c._1, c._2 + l, c._3 + 1)
          },



          //这个表示对于所有的处理好的样本（均为c类型）进行聚合。
          combOp = (c1, c2) => {
            // c: (grad, loss, count)
            //即对应的梯度向量值相加，对应的损失和相加，对应的计数值相加。最后一个参数表示的是样本数量
            (c1._1 += c2._1, c1._2 + c2._2, c1._3 + c2._3)
          })


      if (miniBatchSize > 0) { //当样本数量大于0的时候。

        /**
         *保存误差，迭代误差=平均损失+正则误差。
         */
        stochasticLossHistory.append(lossSum / miniBatchSize + regVal)  //这个表示迭代完成后将误差加入误差数组。
         //其中的损失为平均损失，即总的损失除以总数量的计数和。

        //调用updater的compute方法来更新梯度值。
        val update = updater.compute(
          weights, Vectors.fromBreeze(gradientSum / miniBatchSize.toDouble),
          //Vectors.fromBreeze(gradientSum / miniBatchSize.toDouble)表示总的梯度和除以数据量表示平均梯度。
          stepSize, i, regParam)  //stepSize表示步长，i表示第i次迭代，regParam表示正则化参数。

        weights = update._1  //将权重更新为update的第一个值,表示的是权重因子
        regVal = update._2  //表示的是正则值

        previousWeights = currentWeights
        currentWeights = Some(weights)
        if (previousWeights != None && currentWeights != None) {
          converged = isConverged(previousWeights.get,
            currentWeights.get, convergenceTol)
        }
      } else {
        logWarning(s"Iteration ($i/$numIterations). The size of sampled batch is zero")
      }
      i += 1
    }

    logInfo("GradientDescent.runMiniBatchSGD finished. Last 10 stochastic losses %s".format(
      stochasticLossHistory.takeRight(10).mkString(", ")))

    (weights, stochasticLossHistory.toArray)  //迭代完成之后，返回的是一个迭代的初始权重和每次迭代的损失数组。

  }

逻辑回归的梯度计算方法

在上述方法中的每次迭代都会调用gradient.compute方法。这个方法用来计算每个样本的梯度和误差。这个compute方法是LogisticGradient类的compute方法。该方法基于最小二乘计算梯度值和损失值。计算每个样本损失值。
就像上面所讲的
每一个样本的梯度的计算方法为：

margin= - w*x

multiplier= (1/ (1+ exp(margin))-y)

gradient = multiplier * x

每个样本的误差，

l o s s = 1 2 * (y - h w (x)) 2

$loss=\frac{1}{2}*(y-h_{w}(x))^{2}$

权重的更新方法：
weight = weight -alpha*(gradient + regParam * weight)

下面的代码只包含二元逻辑回归

class LogisticGradient(numClasses: Int) extends Gradient {

  def this() = this(2)  //默认的逻辑回归的类别标签数是2

  override def compute(data: Vector, label: Double, weights: Vector): (Vector, Double) = {
    val gradient = Vectors.zeros(weights.size)  //初始化一个梯度向量
    val loss = compute(data, label, weights, gradient) //调用下面的compute方法，最后返回一个元祖
    (gradient, loss)
  }

  //梯度，损失计算，计算好的梯度值，会被添加到这个cumGradient向量中。，然后返回损失值
  override def compute(
      data: Vector,
      label: Double,
      weights: Vector,
      cumGradient: Vector): Double = {
    val dataSize = data.size

    // (weights.size / dataSize + 1) is number of classes
    require(weights.size % dataSize == 0 && numClasses == weights.size / dataSize + 1)
    numClasses match {

      case 2 =>  //二元逻辑回归

               val margin = -1.0 * dot(data, weights)  //首先是计算margin
        //margin=-w*x。即权重和数据相乘。

        val multiplier = (1.0 / (1.0 + math.exp(margin))) - label
        //multiplier=(1/(1+exp(margin))-y)，表示估计值减去实际值

        axpy(multiplier, data, cumGradient) //然后进行梯度的计算cumGradient=multiplier*x+cumGradient

        if (label > 0) {
          // The following is equivalent to log(1 + exp(margin)) but more numerically stable.
        //用来计算log(1+exp(margin))
          MLUtils.log1pExp(margin) //表示用来计算损失值
        } else {
          //log(1+exp(margin))-margin
          MLUtils.log1pExp(margin) - margin
        }
       }

最后是权重的更新，Updater。LogisticRegression实现的是L2的正则化更新
SquaredL2Updater。这个类在optimization的updater类里面。
权重的更新方法：
weight = weight -alpha*(gradient + regParam * weight)

/*
* L2正则化更新梯度，L2正则化：
*   
*    R(W)=1/2*||w||^2
* 利用step-size/sqrt(iterations)老作为更新的系数
*L2正则化的更新公式为：weight=weight-(stepsize/sqrt(iters))*(gradient+regParam*weight)
*/
class SquaredL2Updater extends Updater {
  override def compute(
      weightsOld: Vector,
      gradient: Vector,
      stepSize: Double,
      iter: Int,
      regParam: Double): (Vector, Double) = {

    // add up both updates from the gradient of the loss (= step) as well as
    // the gradient of the regularizer (= regParam * weightsOld)
    // w' = w - thisIterStepSize * (gradient + regParam * w)
    // w' = (1 - thisIterStepSize * regParam) * w - thisIterStepSize * gradient

    val thisIterStepSize = stepSize / math.sqrt(iter) //首先是去计算alpha值
    //thisIterStepSize表示的是alpha值，它是迭代次数的(-1/2)次方，因此会随着迭代次数的增加而逐渐减小，这样也保证了刚开始迭代的时候学习速率比较快，后期比较慢。

    val brzWeights: BV[Double] = weightsOld.toBreeze.toDenseVector //转变权重为密集向量

    //注意下面这个公式，mllib在这里首先做了拆分，首先计算公式(2)的第一部分，然后是第二部分
      // w' = w - thisIterStepSize * (gradient + regParam * w) (1)
    // w' = (1 - thisIterStepSize * regParam) * w - thisIterStepSize * gradient (2)

    brzWeights :*= (1.0 - thisIterStepSize * regParam) //首先计算的是w=(1 - thisIterStepSize * regParam) * w
    //冒号相乘表示的是追乘操作
    brzAxpy(-thisIterStepSize, gradient.toBreeze, brzWeights) //接下来计算的是brzWeights-thisIterStepSize * gradient

    val norm = brzNorm(brzWeights, 2.0) //接下来计算L2范数

   //最终返回的是权重向量和L2范数乘以0.5倍的正则因子
    (Vectors.fromBreeze(brzWeights), 0.5 * regParam * norm * norm)

  }

最后是逻辑回归模型类LogisticRegressionModel.这个类包含有predictPoint方法。用来进行节点的预测。然后还包含有最基本的模型加载方法和模型保存方法。这个博客只包含预测方法。

class LogisticRegressionModel @Since("1.3.0") (
    @Since("1.0.0") override val weights: Vector, //权重每个特征的权重
    @Since("1.0.0") override val intercept: Double,  //偏置项
    @Since("1.3.0") val numFeatures: Int, ///特征的维度
    @Since("1.3.0") val numClasses: Int)  //标签的类别数，也即分类数，默认是二元分类
  extends GeneralizedLinearModel(weights, intercept) with ClassificationModel with Serializable
  with Saveable with PMMLExportable

这个model类继承了广义线性模型GeneralizedLinearModel的predict方法。

 def predict(testData: RDD[Vector]): RDD[Double] = { //testData为测试数据集，
    // A small optimization to avoid serializing the entire model. Only the weightsMatrix
    // and intercept is needed.
    val localWeights = weights 
    val bcWeights = testData.context.broadcast(localWeights) //获取权重，广播权重
    val localIntercept = intercept
    testData.mapPartitions { iter =>   //在每个分区得带这个权重值，调用predictPoint方法进行预测
      val w = bcWeights.value
      iter.map(v => predictPoint(v, w, localIntercept))
    } 
  }

注意这个predictpoint方法，这个predictpoint方法只是GeneralizedLinearModel广义线性模型这个抽象类的一个抽象方法，而LogisticRegressionModel实现了这个predictPoint方法。
当大于我们设置的阈值的时候，我们将其视作1，否则视作0.

//这个模型类继承自广义线性回归的抽象模型类，实现了里面的predicPoint方法
  override protected def predictPoint(
      dataMatrix: Vector,
      weightMatrix: Vector,
      intercept: Double) = {
    //保证dataMatrix和features的特征的维度相同。
    require(dataMatrix.size == numFeatures)

    // If dataMatrix and weightMatrix have the same dimension, it's binary logistic regression.
    //如果是二元回归的话，
    if (numClasses == 2) {
      //这个是W*X+intercept(加上偏置)做向量的点乘
      val margin = dot(weightMatrix, dataMatrix) + intercept
      //这个是计算s函数 1/(1+exp(-margin))
      val score = 1.0 / (1.0 + math.exp(-margin))

      threshold match {
         //如果得到的这个score大于阈值，则设置为1，否则设置为0
        case Some(t) => if (score > t) 1.0 else 0.0
        case None => score //如果没有设置阈值，则返回计算的score，默认的阈值为0.5.
      }
    }

参考链接:
* 逻辑回归
* coursera公开笔记-逻辑回归
* 逻辑回归理论简介

stevekangpei

关注

7
点赞
踩
12

收藏

觉得还不错? 一键收藏
1
评论
Spark MLlib 源代码解析之逻辑回归LogisticRegression

Spark MLlib 逻辑回归 LogisticRegression模型逻辑回归的本质是线性回归，只是在特征到结果的过程上加上了一层映射。即首先需要把特征进行求和，然后将求和后的结果应用于一个g(z)函数。g(z)可以将值映射到0或者是1上面。逻辑回归和多重线性回归有很多的相似之处。最大的区别是他们的因变量不同。这两个回归也可以统一归为广义线性模型。在spark mllib实现的过程中也是先定义好
复制链接

扫一扫