逻辑回归梯度下降法详解

最新推荐文章于 2023-09-23 06:13:37 发布

cjneo

最新推荐文章于 2023-09-23 06:13:37 发布

阅读量779

点赞数

原文地址 http://blog.csdn.net/lookqlp/article/details/51161640

引言

逻辑回归常用于预测疾病发生的概率，例如因变量是是否恶性肿瘤，自变量是肿瘤的大小、位置、硬度、患者性别、年龄、职业等等（很多文章里举了这个例子，但现代医学发达，可以通过病理检查，即获取标本放到显微镜下观察是否恶变来判断）；广告界中也常用于预测点击率或者转化率(cvr/ctr)，例如因变量是是否点击，自变量是物料的长、宽、广告的位置、类型、用户的性别、爱好等等。
本章主要介绍逻辑回归算法推导、梯度下降法求最优值的推导及spark的源码实现。

常规方法

一般回归问题的步骤是：
1. 寻找预测函数（h函数，hypothesis）
2. 构造损失函数（J函数）
3. 使损失函数最小，获得回归系数θ

而第三步中常见的算法有：
1. 梯度下降
2. 牛顿迭代算法
3. 拟牛顿迭代算法（BFGS算法和L-BFGS算法）
其中随机梯度下降和L-BFGS在spark mllib中已经实现，梯度下降是最简单和容易理解的。

推导

二元逻辑回归

构造预测函数
$h θ (x) = g (θ T x) = 1 1 + e - θ T x$
其中：
$θ T x = \sum i = 1 n θ i x i = θ 0 + θ 1 x 1 + θ 2 x 2 + . . . + θ n x n θ = ⎡ ⎣ ⎢ ⎢ ⎢ θ 0 θ 1 . . . θ n ⎤ ⎦ ⎥ ⎥ ⎥, x = ⎡ ⎣ ⎢ ⎢ ⎢ x 0 x 1 . . . x n ⎤ ⎦ ⎥ ⎥ ⎥$
为何LR模型偏偏选择sigmoid 函数呢？逻辑回归不是回归问题，而是二分类问题，因变量不是0就是1，那么我们很自然的认为概率函数服从伯努利分布，而伯努利分布的指数形式就是个sigmoid 函数。
函数hθ(x)表示结果取1的概率，那么对于分类1和0的概率分别为：
$P (y = 1 | x; θ) = h θ (x) P (y = 0 | x; θ) = 1 - h θ (x)$
概率一般式为：
$P (y | x; θ) = (h θ (x)) y ((1 - h θ (x))) 1 - y$
最大似然估计的思想
当从模型总体随机抽取m组样本观测值后，我们的目标是寻求最合理的参数估计θ′使得从模型中抽取该m组样本观测值的概率最大。最大似然估计就是解决此类问题的方法。求最大似然函数的步骤是：
1. 写出似然函数
2. 对似然函数取对数
3. 对对数似然函数的参数求偏导并令其为0，得到方程组
4. 求方程组的参数
为什么第三步要取对数呢，因为取对数后，乘法就变成加法了，且单调性一致，不会改变极值的位置，后边就更好的求偏导。
构造损失函数
线性回归中的损失函数是：
$J (θ) = 1 2 m \sum i = 1 m (y i - h θ (x i)) 2$
线性回归损失函数有很明显的实际意义，就是平方损失。而逻辑回归却不是，它的预测函数hθ(x)明显是非线性的，如果类比的使用线性回归的损失函数于逻辑回归，那J(θ)很有可能就是非凸函数，即存在很多局部最优解，但不一定是全局最优解。我们希望构造一个凸函数，也就是一个碗型函数做为逻辑回归的损失函数。
按照求最大似然函数的方法，逻辑回归似然函数：
$L (θ) = \prod i = 1 m P (y i | x i; θ) = \prod i = 1 m (h θ (x i)) y i ((1 - h θ (x i))) 1 - y i$
其中m表示样本数量，取对数：
$l (θ) = l o g L (θ) = \sum i = 1 m (y i l o g h θ (x i) + (1 - y i) l o g (1 - h θ (x i)))$
我们的目标是求最大l(θ)时的θ，如上函数是一个上凸函数，可以使用梯度上升来求得最大似然函数值(最大值)。或者上式乘以-1，变成下凸函数，就可以使用梯度下降来求得最小负似然函数值（最小值）：
$J (θ) = - 1 m l (θ)$
同样是取极小值，思想与损失函数一致，即我们把如上的J(θ)作为逻辑回归的损失函数。Andrew Ng的课程中，上式乘了一个系数1/m，我怀疑就是为了和线性回归的损失函数保持一致吧。
求最小值时的参数
我们求最大似然函数参数的第三步时，令对参数θ偏导=0，然后求解方程组。考虑到参数数量的不确定，即参数数量很大，此时直接求解方程组的解变的很困难，或者根本就求不出精确的参数。于是，我们用随机梯度下降法，求解方程组的值。
当然也可以使用牛顿法、拟牛顿法。梯度下降法是最容易理解和推导的，如下是推导过程：
梯度下降θ的更新过程，走梯度方向的反方向：
$θ j : = θ j - α δ δ θ j J (θ)$
其中：
$δ δ θ j J (θ) = - 1 m \sum i = 1 m (y i 1 h θ ( x i ) δ δ θ j h θ (x i) - (1 - y i) 1 1 - h θ ( x i ) δ δ θ j h θ (x i)) = - 1 m \sum i = 1 m (y i 1 g ( θ T x i ) - (1 - y i) 1 1 - g ( θ T x i )) δ δ θ j g (θ T x i) = - 1 m \sum i = 1 m (y i 1 g ( θ T x i ) - (1 - y i) 1 1 - g ( θ T x i )) g (θ T x i) (1 - g (θ T x i)) δ δ θ j θ T x i = - 1 m \sum i = 1 m (y i (1 - g (θ T x i)) - (1 - y i) g (θ T x i)) x j i = - 1 m \sum i = 1 m (y i - g (θ T x i)) x j i = 1 m \sum i = 1 m (h θ (x i) - y i)) x j i$
第二步推导请注意：

那么可以推导：
$δ δ θ j g (θ T x i) = - e - θ T x i ( 1 + e - θ T x i ) 2 δ δ θ j (- 1) θ T x i = g (θ T x i) (1 - g (θ T x i)) δ θ j θ T x i$
因此更新过程可以写成：
$θ j : = θ j - α 1 m \sum i = 1 m (h θ (x i) - y i)) x j i$
那迭代多少次停止呢，spark是指定迭代次数和比较两次梯度变化或者cost变化小于一定值时停止。
过拟合问题
过拟合问题，即我们求得的回归系数在实验集中效果很好，但之外的数据效果很差。机器学习中的特征基本上是靠人的经验选择的，有可能某一些特征或者特征组合与因变量没有任何关系，即某些θi≈0。所以我们需要把不必要的特征剔除，一般我们使用正则化来保留所有特征，并让它相应的系数≈0，L1范数正则化后θ的更新：
$θ j : = θ j - α 1 m \sum i = 1 m (h θ (x i) - y i)) x j i - λ m θ j$
λ越大，对模型的复杂度惩罚越大，有可能出现欠拟合现象。λ越小，惩罚越小，可能新出现过拟合现象。spark逻辑回归的随机梯度下降法中，使用的是L2范数正则化。

多元逻辑回归

推广到K元逻辑回归，即因变量为0、1、2、…、k-1。在二元逻辑回归中有这样的性质：

l o g P ( y = 1 | x , θ ) P ( y = 0 | x , θ ) = θ T x

推广至K元逻辑回归：

l o g P ( y = 1 | x , θ ) P ( y = 0 | x , θ ) = θ T 1 x l o g P ( y = 2 | x , θ ) P ( y = 0 | x , θ ) = θ T 2 x . . . l o g P ( y = K - 1 | x , θ ) P ( y = 0 | x , θ ) = θ T K - 1 x

其中，

θ=(θ1,θ2,...,θK−1)T ，是个（k-1）*（n+1）的矩阵，n为特征的个数，加1是增加截距项。去除对数则得到概率分布：

P (y = 0 | x, θ) = 1 1 + \sum K - 1 i = 1 e θ T i x P (y = 1 | x, θ) = e θ T 1 x 1 + \sum K - 1 i = 1 e θ T i x . . . P (y = K - 1 | x, θ) = e θ T K - 1 x 1 + \sum K - 1 i = 1 e θ T i x

K元逻辑回归似然函数：

L (θ) = \prod i = 1 m P (y | x, θ)

定义：

α (y i) = 1 i f y i = 0 α (y i) = 0 i f y i \neq 0

取对数：

l (θ, x) = \sum i = 1 m l o g P (y i | x i, θ) = \sum i = 1 m α (y i) l o g P (y = 0 | x i, θ) + (1 - α (y i)) l o g P (y i | x i, θ) = \sum i = 1 m α (y i) l o g 1 1 + \sum K - 1 k = 1 e θ T k x + (1 - α (y i)) l o g e θ T y i x 1 + \sum K - 1 k = 1 e θ T k x = \sum i = 1 m (1 - α (y i)) θ T y i x - l o g (1 + \sum k = 1 K - 1 e θ T k x)

同样的我们得到损失函数：

J (θ, x) = - 1 m l (θ, x)

θ 更新过程：

θ j : = θ j - α δ δ θ j J (θ, x)

对

θ 求偏导得到梯度：

G k j (θ, x) = - 1 m δ l ( θ , x ) δ θ k j = - 1 m (\sum i = 1 m (1 - α (y i)) x i j δ k, y i - e θ T x i 1 + e θ T x i x i j)

其中k表示因变量，j表示特征数量，i表示实验数。
spark源码注释中，稍稍不一样，

l(w,x) 乘以了-1，其实与我们上边推导的

−1m 意义一样。我们来看看spark的推导过程。

P (y = 0 | x, w) = 1 / (1 + \sum i K - 1 exp (x w i)) P (y = 1 | x, w) = e x p (x w 1) / (1 + \sum i K - 1 exp (x w i)) . . . P (y = K - 1 | x, w) = e x p (x w K - 1) / (1 + \sum i K - 1 exp (x w i)

取对数：

l (w, x) = - l o g P (y | x, w) = - α (y) l o g P (y = 0 | x, w) - (1 - α (y)) l o g P (y | x, w) = l o g (1 + \sum i K - 1 exp (x w i)) - (1 - α (y)) x w y - 1 = l o g (1 + \sum i K - 1 exp (m a r g i n s i)) - (1 - α (y)) m a r g i n s y - 1

其中：

α (i) = 1 i f i! = 0, α (i) = 0 i f i = = 0, m a r g i n s i = x w i .

求偏导：

\partial l ( w , x ) \partial w i j = (exp (x w i) / (1 + \sum k K - 1 exp (x w k)) - (1 - α (y) δ y, i + 1)) * x j = m u l t i p l i e r i * x j

其中：

δ i, j = 1 i f i = = j, δ i, j = 0 i f i! = j, m u l t i p l i e r = exp (m a r g i n s i) / (1 + \sum k K - 1 exp (m a r g i n s i)) - (1 - α (y) δ y, i + 1)

为了不让数值溢出,xw项减了maxMargin，

l(w,x) 改写为：

l (w, x) = l o g (1 + \sum i K - 1 exp (m a r g i n s i)) - (1 - α (y)) m a r g i n s y - 1 = l o g (exp (- m a x M a r g i n) + \sum i K - 1 exp (m a r g i n s i - m a x M a r g i n)) + m a x M a r g i n - (1 - α (y)) m a r g i n s y - 1 = l o g (1 + s u m) + m a x M a r g i n - (1 - α (y)) m a r g i n s y - 1

其中：

s u m = exp (- m a x M a r g i n) + \sum i K - 1 exp (m a r g i n s i - m a x M a r g i n) - 1

而multiplier可以表示为：

m u l t i p l i e r = exp (m a r g i n s i) / (1 + \sum k K - 1 exp (m a r g i n s i)) - (1 - α (y) δ y, i + 1) = exp (m a r g i n s i - m a x M a r g i n) / (1 + s u m) - (1 - α (y) δ y, i + 1)

spark源码

先看实例代码:

import org.apache.spark.SparkContext
import org.apache.spark.mllib.classification.{LogisticRegressionWithLBFGS, LogisticRegressionModel}
import org.apache.spark.mllib.evaluation.MulticlassMetrics
import org.apache.spark.mllib.regression.LabeledPoint
import org.apache.spark.mllib.linalg.Vectors
import org.apache.spark.mllib.util.MLUtils

// Load training data in LIBSVM format.
//样例数据格式:
//1 特征id1:值id1 特征id2:值id2 ...
//0 特征id1:值id3 特征id4:值id4 ...
//特征和特征对应的值都使用数值一一标示了
val data = MLUtils.loadLibSVMFile(sc, "data/mllib/sample_libsvm_data.txt")

// 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
// 官方样例分类数设置为10,但样例数据因变量是0和1,所以这里应该时设置错了.
// 梯度下降法每次迭代都会变量整个样本集,推荐使用拟牛顿法LBFGS,后续文章中继续介绍
val model = new LogisticRegressionWithLBFGS()
  .setNumClasses(10)
  .run(training)

// Compute raw scores on the test set.
val predictionAndLabels = test.map { case LabeledPoint(label, features) =>
  val prediction = model.predict(features)
  (prediction, label)
}

// Get evaluation metrics.
val metrics = new MulticlassMetrics(predictionAndLabels)
val precision = metrics.precision
println("Precision = " + precision)

// Save and load model
// 输出是个模型,就是一个向量$\theta$,带入概率分布函数求得类型的概率
model.save(sc, "myModelPath")
val sameModel = LogisticRegressionModel.load(sc, "myModelPath")

随机梯度下降调用:

 /**
   * Train a logistic regression model given an RDD of (label, features) pairs. We run a fixed
   * number of iterations of gradient descent using the specified step size. Each iteration uses
   * `miniBatchFraction` fraction of the data to calculate the gradient. The weights used in
   * gradient descent are initialized using the initial weights provided.
   * NOTE: Labels used in Logistic Regression should be {0, 1}
   *
   * @param input RDD of (label, array of features) pairs.
   * @param numIterations Number of iterations of gradient descent to run.迭代次数
   * @param stepSize Step size to be used for each iteration of gradient descent.步长
   * @param miniBatchFraction Fraction of data to be used per iteration.用于模型预估数据的比例
   * @param initialWeights Initial set of weights to be used. Array should be equal in size to the number of features in the data.初始化权重
   */
  @Since("1.0.0")
  def train(
      input: RDD[LabeledPoint],
      numIterations: Int,
      stepSize: Double,
      miniBatchFraction: Double,
      initialWeights: Vector): LogisticRegressionModel = {
    new LogisticRegressionWithSGD(stepSize, numIterat2 Aions, 0.0, miniBatchFraction)
      .run(input, initialWeights)
  }

LogisticRegressionWithLBFGS和LogisticRegressionWithSGD都继承于GeneralizedLinearModel,它的run方法:

 def run(input: RDD[LabeledPoint], initialWeights: Vector): M = {

    if (numFeatures < 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:
     *
     * During the optimization process, the convergence (rate) depends on the condition number of
     * the training dataset. Scaling the variables often reduces this condition number
     * heuristically, thus improving the convergence rate. Without reducing the condition number,
     * some training datasets mixing the columns with different scales may not be able to converge.
     *
     * GLMNET and LIBSVM packages perform the scaling to reduce the condition number, and return
     * the weights in the original scale.
     * See page 9 in http://cran.r-project.org/web/packages/glmnet/glmnet.pdf
     *
     * Here, if useFeatureScaling is enabled, we will standardize the training features by dividing
     * the variance of each column (without subtracting the mean), and train the model in the
     * scaled space. Then we transform the coefficients from the scaled space to the original scale
     * as GLMNET and LIBSVM do.
     *通过每一列除以这一列的标准差,将数据标准化.LBFGS算法中可以启用.
     * Currently, it's only enabled in LogisticRegressionWithLBFGS
     */
    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))
        }
      }

    /**
     * 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
    }
    //SGD 或者 LBFGS算法
    val weightsWithIntercept = optimizer.optimize(data, initialWeightsWithIntercept)
    ...
    createModel(weights, intercept)
  }

梯度下降SGD实现:

def runMiniBatchSGD(
      data: RDD[(Double, Vector)],
      gradient: Gradient,
      updater: Updater,
      stepSize: Double,
      numIterations: Int,
      regParam: Double,
      miniBatchFraction: Double,
      initialWeights: Vector,
      convergenceTol: Double): (Vector, Array[Double]) = {
    ...
      //不知道此数组干啥用的
    val stochasticLossHistory = new ArrayBuffer[Double](numIterations)
    ...
    // Initialize weights as a column vector
    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
    while (!converged && i <= numIterations) {
      val bcWeights = data.context.broadcast(weights)
      // Sample a subset (fraction miniBatchFraction) of the total data
      // compute and sum up the subgradients on this subset (this is one map-reduce)
      val (gradientSum, lossSum, miniBatchSize) = data.sample(false, miniBatchFraction, 42 + i)
        .treeAggregate((BDV.zeros[Double](n), 0.0, 0L))(
          seqOp = (c, v) => {
            // c: (grad, loss, count), v: (label, features)
            // 返回损失loss,没看明白为何要算loss,及loss为何这么算log(1 + exp(margin))
            // 主要目的时计算c._1梯度向量
            val l = gradient.compute(v._2, v._1, bcWeights.value, Vectors.fromBreeze(c._1))
            (c._1, c._2 + l, c._3 + 1)
          },
          combOp = (c1, c2) => {
            // c: (grad, loss, count)
            (c1._1 += c2._1, c1._2 + c2._2, c1._3 + c2._3)
          })

      if (miniBatchSize > 0) {
        /**
         * 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)
        // 正则化
        val update = updater.compute(
          weights, Vectors.fromBreeze(gradientSum / miniBatchSize.toDouble),
          stepSize, i, regParam)
        weights = update._1
        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)

  }

combOp是θ的更新过程中的∑过程.在二元逻辑回归情况下:

case 2 =>
        /**
         * For Binary Logistic Regression.
         *
         * Although the loss and gradient calculation for multinomial one is more generalized,
         * and multinomial one can also be used in binary case, we still implement a specialized
         * binary version for performance reason.
         */
        val margin = -1.0 * dot(data, weights)
        val multiplier = (1.0 / (1.0 + math.exp(margin))) - label
        axpy(multiplier, data, cumGradient)
        if (label > 0) {
          // The following is equivalent to log(1 + exp(margin)) but more numerically stable.
          MLUtils.log1pExp(margin)
        } else {
          MLUtils.log1pExp(margin) - margin
        }

margin 就是(−θTx),而multiplier就是hθ(xi)−yi.axpy方法就是(hθ(xi)−yi))xi.

cjneo

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
逻辑回归梯度下降法详解

原文地址 http://blog.csdn.net/lookqlp/article/details/51161640引言逻辑回归常用于预测疾病发生的概率，例如因变量是是否恶性肿瘤，自变量是肿瘤的大小、位置、硬度、患者性别、年龄、职业等等（很多文章里举了这个例子，但现代医学发达，可以通过病理检查，即获取标本放到显微镜下观察是否恶变来判断）；广告界中也常用于预测点击率或者转化率(cvr/ctr)，例如...
复制链接

扫一扫