Spark MLib 每周一算法 —— LogisticRegression

2018-02-19

LogisticRegression(逻辑回归)

LogisticRegression是一个Estimator，是常用的分类算法之一。输入带有label和features列的训练集得到分类器模型；LogisticRegressionModel是spark对训练集fit后的模型对象，该对象中存放了各个线性属性的系数值和LogisticRegressionSummary对象(Spark2.2.0暂时只实现了二分类的summary)，summary对象理存放了模型的一些性能指标信息，比如ROC等；最终transform新数据的输出包含了原有输入信息，线性回归的因变量值(raw prediction)，作用逻辑函数/Softmax函数后的几率值(probability)以及最终预测的分类(prediction)。

伪算法过程

模型
逻辑回归是广义线性模型的一种，当广义线性模型的链接函数(link function)为对数几率函数(logit)时，就被称之为逻辑回归。下式左侧既是logit函数，右侧是属性变量X和系数β的向量积。
$$\ln(\frac{p}{1-p})=X\beta$$
logit函数的反函数是logistic函数，公式如下：
$$\frac{e^{X\beta}}{1+e^{X\beta}}$$
logistic函数又被称作Sigmoid函数，它的特点就是函数值的变化在原点附近区域非常陡峭，迅速接近0或者1，相当于将回归后的结果映射到一个0/1的值域中，从而完成对数据的二分类。

迭代算法
Spark采用LBFGS/OWLQN迭代算法来进行回归，这2个算法都属于拟牛顿法在大规模数据时的优化，本质是利用牛顿法寻找目标函数的极小值，既凸优化问题的解。目标函数通常由损失函数(Loss function)和正则化惩罚函数(Regularization term)2部分组成，一般形式如下：
$$\frac{1}{n}\sum_{i=1}^{n}L(\beta x_{i}, y_{i})+\lambda R(\beta)$$
二分类逻辑回归时使用的损失函数为，其中y的取值为0或者1，这个函数又称最大似然估计函数(Maximum likelihood estimation)：
$$y_{i}log\frac{e^{X\beta}}{1+e^{X\beta}}+(1-y_{i})log(1-\frac{e^{X\beta}}{1+e^{X\beta}})$$
具体使用的迭代算法和正则化部分有关系，当不使用正则化(Regularization)或者使用L2作为正则化惩罚函数(L2 penalty，又称L2 regularization)的时候，采用LBFGS(Limited-memory Broyden–Fletcher–Goldfarb–Shanno)算法，其他情况比如L1正则化或者非0的elasticNet导致目标函数不可微，则使用OWLQN(Orthant-Wise Limited-memory Quasi-Newton)算法。

正则化
正则化(Regularization)是用来防止模型过拟合(overfitting)的。Spark提供了2个超参来针对正则化做模型调整——regParam和elasticNetParam。上面回归目标函数中正则化部分中的λ就是可以调整的参数regParam，越大表示对模型训练集的拟合系数的惩罚程度越大，系数越接近或者等于0，模型就越简单，但由于和原始的loss函数相差大了，精度就下降了，可能会导致欠拟合(underfitting)的结果。
同时Spark采用的elasticNet方法相当于同时集成了L1和L2正则化，并提供了动态调整的能力，elasticNet的惩罚函数如下：
$$R(\beta)=\alpha||\beta||+(1-\alpha)\frac{1}{2}||\beta||_{2}^{2}$$
上面的α就是参数elasticNetParam，当α是0的时候惩罚函数就是模型系数的平方和除以2n(系数向量的2范数平方的一半)，这时的回归也被称为岭回归(ridge regression)；当α是1的时候惩罚函数就是模型系数的绝对值之和除以n(系数向量的1范数)，这时的回归也被称为Lasso回归(least absolute shrinkage and selection operator)。我们可以动态在0～1之间调整α，越接近1的时候，模型稀疏化的能力越强，系数接近0的个数越多，相当于做了属性选择。

多分类
Spark同时支持多分类的逻辑回归(Multinomial logistic regression)，将二分类逻辑回归做扩展，把一组线性系数和自变量计算后的值通过Softmax函数映射成在0，1范围内，并且所有分类和为1的数值，从而决定这组自变量属于哪个分类。所以又称为Softmax回归，Softmax函数形式如下，其中K为待分类的类别数，这时的系数β不再是个向量，而是一个包含K个向量矩阵，最终的计算结果是一个K维的向量，每个分量代表一个分类label的可能值：
$$\sigma(X\beta)=\frac{e^{X\beta_{j}}}{\sum_{k=1}^{K}e^{X\beta_{k}}}\qquad j=1, …, K$$

类似的损失函数如下形式，其中1{y=j}为示性函数(y等于j时为1，否则为0)：
$$\sum_{j=1}^{K}1\{y=j\}log\frac{e^{X\beta_{j}}}{\sum_{k=1}^{K}e^{X\beta_{k}}}$$

截距
Spark默认使用带截距(intercept term)的线性回归，就是说模型系数里会带有一个常量值，当然也可以在确保实际数据必过原点的情况下设置采用不带截距的线性回归。

权重
Spark支持带权重的回归模型，权重值相当于对训练集数据做了采样(Sampling)，比如对于比较常见的不平衡(imbalanced)数据的二分类逻辑回归(比如疾病诊断)，训练集中label=1的数据占了95%，那就可以设置label=1的数据的权重为0.05，而label=0的数据权重为0.95，这样可以避免模型对应该为label=0的数据大量误判。

标准化
Spark支持对训练集的属性数据先做标准化(Standardization)处理，线性回归的解本身对数据的scale程度没有依赖，做标准化处理后的解也与原始数据等价。单通常会做标准化处理，这样会使得迭代算法达到更好的收敛效果。Spark给出的最终线性系数是基于原始输入数据的。

对系数的限制
Spark可以在回归时对系数的取值范围做出限制，这时采用LBFGS的优化算法LBFGS-B(bound-constrained optimization)。

参数

输入输出相关:

labelCol: 输入训练集中标记字段的名称 (默认值: label)
featuresCol: 输入训练集和预测数据集中的属性向量的字段名称(默认值: features)
weightCol: 输入训练集中权重字段的名称，如果输入的训练集包含权重字段，则采用带权重的回归模型，如果不包含，相当于权重都是1.0 (没有默认值)
rawPredictionCol: 输出结果数据中回归因变量值的字段名称 (默认值: rawPrediction)
probabilityCol: 输出结果数据中逻辑函数/Softmax函数几率值的字段名称 (默认值: probability)
predictionCol: 输出结果数据中最终判别类别的字段名称 (默认值: prediction)

算法模型相关:

family: 采用二分类还是多分类(默认值：auto)
可选项1. binomial-二分类逻辑回归；
可选项2. multinomial-多分类逻辑回归；
可选项3. auto根据输入训练集的label字段中的数据自动选择，如果label字段最大数据max>=2，则采用0～max的多分类逻辑回归，否则就是0或者1的二分类逻辑回归。
fitIntercept: 是否使用带截距的回归。(boolean类型，默认值: true)
maxIter: 迭代算法最大迭代次数，和下面的tol值一起来控制拟牛顿算法的结束条件。(>0的整数，默认值: 100)
tol: 收敛容忍系数(convergence tolerance)，算法每次迭代后的比较阈值来确定算法是否结束，值越小，可能执行的迭代次数就越多。(>0的实数，默认值: 1.0E-6)
regParam: 正则化惩罚程度参数，具体使用方式见上面正则化相关介绍。(>=0的实数，默认值：0)
elasticNetParam: 正则化时弹性网络调整参数，具体使用方式见上面正则化相关介绍。([0,1]之间，默认值：0)
threshold: 二分类算法中根据该值将计算得到逻辑函数几率值映射到分类类别上。([0,1]之间，默认值：0.5)
thresholds: 多分类算法中根据该组数组值将计算得到Softmax函数几率值映射到分类类别上，该数组长度必须和类别数量吻合，最终选择的分类是几率值除以该分类的threshold值最大的那个，如果没有设置，相当于直接取几率最大的。(无默认值)
standardization: 是否在回归前对训练集中的属性数据做标准化处理。(boolean类型，默认值: true)
lowerBoundsOnCoefficients: 系数的下限约束，能够支持设置多分类(矩阵中每个向量代表一个分类系数)的系数。(矩阵类型，无默认值)
upperBoundsOnCoefficients: 系数的上限约束，能够支持设置多分类(矩阵中每个向量代表一个分类系数)的系数。(矩阵类型，无默认值)
lowerBoundsOnIntercepts: 常量系数(截距)的下限约束，能够支持设置多分类(向量中每个分量代表一个分类系数)的截距。(向量类型，无默认值)
upperBoundsOnIntercepts: 常量系数(截距)的上限约束，能够支持设置多分类(向量中每个分量代表一个分类系数)的截距。(向量类型，无默认值)

Spark实现相关:

aggregationDepth: Spark迭代计算过程中才用treeAggregate的方式来reduce各个分区数据，这里的参数可以用来调整treeAggregate的层次深度，层次越深reduce时用到的执行节点更多，能并行的数据量就越多。(>=2的整数，默认值2)

model对象的成员

numClasses: 返回模型有多少个分类。
numFeatures: 返回模型有多少个属性。
coefficients: 二分类逻辑回归模型的系数向量。
coefficientMatrix: 多分类逻辑回归模型的系数矩阵。
intercept: 二分类逻辑回归模型的常量系数(截距)。
interceptVector: 多分类逻辑回归模型的常量系数向量。
evaluate(testData): 针对某个带label字段的测试数据集应用回归模型，生成并返回LogisticRegressionSummary对象，相当于针对测试集查看模型性能。(多分类回归模型在当前Spark2.2.0版本里返回的summary对象依旧是二分类的对象，因为整个threshold的情况不同，所以这个summary对象里面除了看模型输出外其他性能属性基本没有意义)
summary: 返回针对训练集的LogisticRegressionTrainingSummary对象，是LogisticRegressionSummary的子类，相当于针对训练集查看模型性能。(Spark2.2.0暂时只有二分类的summary对象，对多分类的模型还没有实现该对象)
hasSummary(): 如果该模型有summary对象，那么该方法返回true(二分类时)，反之返回false(多分类时，后面应该会实现)。

summary对象的成员

不管是evaluate方法返回的summary对象还是model自带的summary属性对象，在现阶段2.2.0版本都是针对二分类BinaryLogisticRegressionSummary或者子类 BinaryLogisticRegressionTrainingSummary的实例，所以具体使用的时候都可以转换(asInstanceOf)成BinaryLogisticRegressionSummary类来使用。

predictions: LogisticRegressionSummary的成员属性，返回模型对数据集做transform的输出结果，是个DataFrame对象，包含原有数据的label，features和结果数据的raw prediction，probability和prediction。
objectiveHistory: LogisticRegressionTrainingSummary的成员属性，返回每次迭代的目标函数取值组成而成的数组。
totalIterations: LogisticRegressionTrainingSummary的类属性，返回了实际的迭代次数。
precisionByThreshold: BinaryLogisticRegressionSummary的成员属性，模型准确率数据，按threshold排列，在二分类问题中特指混淆矩阵(Confusion Matrix)中的真阳性(TP)在所有模型判断为阳性(TP+FP)的比例。
recallByThreshold: BinaryLogisticRegressionSummary的成员属性，模型识别率数据，按threshold排列，在二分类问题中特指混淆矩阵(Confusion Matrix)中的真阳性(TP)在所有实际阳性数据(TP+FN)中的比例。
fMeasureByThreshold: BinaryLogisticRegressionSummary的成员属性，模型准确率(precision)和识别率(recall)综合考量数据，公式如下，其中P为准确率，R为识别率：
$$F_{\beta}=\frac{(\beta^2+1)PR}{\beta^2P+R}$$
Spark的fMearsure为F1，既上式中β=1的情况：$\frac{2PR}{P+R}$ 返回的值同样是按对应的threshold排列。
pr: BinaryLogisticRegressionSummary的成员属性，同一个threshold下模型识别率(recall)和准确率(precision)的对照数据，可以根据这个数据以reacll为横坐标，precision为纵坐标画出每一个不同threashold的点(默认通过点(0,1))，然后近似的画出pr曲线(P/R Curve)，越接近点(1,1)的曲线通常被认为模型性能越好，可以用曲线下的面积(AUPR)来衡量曲线接近点(1,1)的程度。
roc: BinaryLogisticRegressionSummary的成员属性，又称为接收者操作特征曲线(receiver operating characteristic curve)，和pr曲线非常类似，也是按threshold排列的横纵坐标2个值的对照数据，区别就是它的横坐标是FPR——混淆矩阵(Confusion Matrix)中的假阳性(FP)在所有实际阴性数据(FP+TN)中的比例，纵坐标是TPR(和recall定义完全一致)——真阳性(TP)在所有实际阳性数据(TP+FN)中的比例，相当于阴性数据判错率和阳性数据盘对率之间的对照。该曲线数据点(默认通过点(0,0)和点(1,1))都应该在三角形<(0,0)~(0,1)~(1,1)>中，因为最差的模型，既所谓的不看属性随机分类的模型对于不同的threshold得到的分类结果，最终数据会按照阳性和阴性数据分布情况落在直线<(0,0)~(1,1)>上。同样的我们也可以用曲线下的面积来衡量模型的性能，越接近点(0,1)，面积越大，模型性能越好。Spark直接提供了返回该值的参数，见下面的成员属性。
areaUnderROC: BinaryLogisticRegressionSummary的成员属性，roc曲线下面积的值。

例子

val training = spark.createDataFrame(Seq(
  (1.0, Vectors.dense(0.0, 1.1, 0.1)),
  (0.0, Vectors.dense(2.0, 1.0, -1.0)),
  (0.0, Vectors.dense(2.0, 1.3, 1.0)),
  (1.0, Vectors.dense(0.0, 1.2, -0.5))
)).toDF("label", "features")
training.show(false)
//    +-----+--------------+
//    |label|features      |
//    +-----+--------------+
//    |1.0  |[0.0,1.1,0.1] |
//    |0.0  |[2.0,1.0,-1.0]|
//    |0.0  |[2.0,1.3,1.0] |
//    |1.0  |[0.0,1.2,-0.5]|
//    +-----+--------------+

val lr = new LogisticRegression().setMaxIter(10).setRegParam(0.01)

val binomialModel = lr.fit(training)
println("Binomial model's numClasses: " + binomialModel.numClasses)
//2
println("Binomial model's numFeatures: " + binomialModel.numFeatures)
//3
println("Binomial model's coefficients: " + binomialModel.coefficients)
//[-3.1009356010205327,2.60821473832145,-0.3801791225430309]
println("Binomial model's intercept: " + binomialModel.intercept)
//0.06817659473873576
println("Binomial model's hasSummary: " + binomialModel.hasSummary)
//true

val summary = binomialModel.summary
summary.predictions.show(false)
//    +-----+--------------+----------------------------------------+----------------------------------------+----------+
//    |label|features      |rawPrediction                           |probability                             |prediction|
//    +-----+--------------+----------------------------------------+----------------------------------------+----------+
//    |1.0  |[0.0,1.1,0.1] |[-2.8991948946380277,2.8991948946380277]|[0.05219337666300758,0.9478066233369924]|1.0       |
//    |0.0  |[2.0,1.0,-1.0]|[3.1453007464378486,-3.1453007464378486]|[0.9587231582899965,0.04127684171000351]|0.0       |
//    |0.0  |[2.0,1.3,1.0] |[3.1231945700274752,-3.1231945700274752]|[0.9578394235295768,0.04216057647042307]|0.0       |
//    |1.0  |[0.0,1.2,-0.5]|[-3.3881238419959914,3.3881238419959914]|[0.03266869266264558,0.9673313073373544]|1.0       |
//    +-----+--------------+----------------------------------------+----------------------------------------+----------+

println("Binomial summary's objectiveHistory: " + summary.objectiveHistory.mkString(","))
//0.6931471805599453,0.4274232924955569,0.12582029480507145,0.1144102682046959,
//0.10884647747685207,0.10840715948469526,0.10831890205457553,0.10826850363446151,
//0.10823620528587521,0.10822178855356387,0.10821879266352699
println("Binomial summary's totalIterations: " + summary.totalIterations)
//11

val binomialSummary = summary.asInstanceOf[BinaryLogisticRegressionSummary]
binomialSummary.precisionByThreshold.show
//    +-------------------+------------------+
//    |          threshold|         precision|
//    +-------------------+------------------+
//    | 0.9673313073373544|               1.0|
//    | 0.9478066233369924|               1.0|
//    |0.04216057647042307|0.6666666666666666|
//    |0.04127684171000351|               0.5|
//    +-------------------+------------------+

binomialSummary.recallByThreshold.show
//    +-------------------+------+
//    |          threshold|recall|
//    +-------------------+------+
//    | 0.9673313073373544|   0.5|
//    | 0.9478066233369924|   1.0|
//    |0.04216057647042307|   1.0|
//    |0.04127684171000351|   1.0|
//    +-------------------+------+

binomialSummary.fMeasureByThreshold.show
//    +-------------------+------------------+
//    |          threshold|         F-Measure|
//    +-------------------+------------------+
//    | 0.9673313073373544|0.6666666666666666|
//    | 0.9478066233369924|               1.0|
//    |0.04216057647042307|               0.8|
//    |0.04127684171000351|0.6666666666666666|
//    +-------------------+------------------+

binomialSummary.pr.show
//    +------+------------------+
//    |recall|         precision|
//    +------+------------------+
//    |   0.0|               1.0|
//    |   0.5|               1.0|
//    |   1.0|               1.0|
//    |   1.0|0.6666666666666666|
//    |   1.0|               0.5|
//    +------+------------------+

binomialSummary.roc.show
//    +---+---+
//    |FPR|TPR|
//    +---+---+
//    |0.0|0.0|
//    |0.0|0.5|
//    |0.0|1.0|
//    |0.5|1.0|
//    |1.0|1.0|
//    |1.0|1.0|
//    +---+---+

println("Binomial summary's areaUnderROC: " + binomialSummary.areaUnderROC)
//1.0

val test = spark.createDataFrame(Seq(
  (1.0, Vectors.dense(-1.0, 1.5, 1.3)),
  (0.0, Vectors.dense(3.0, 2.0, -0.1)),
  (1.0, Vectors.dense(0.0, 2.2, -1.5))
)).toDF("label", "features")
test.show
//    +-----+--------------+
//    |label|      features|
//    +-----+--------------+
//    |  1.0|[-1.0,1.5,1.3]|
//    |  0.0|[3.0,2.0,-0.1]|
//    |  1.0|[0.0,2.2,-1.5]|
//    +-----+--------------+

binomialModel.evaluate(test).predictions.show(false)
//    +-----+--------------+--------------------------------------+------------------------------------------+----------+
//    |label|features      |rawPrediction                         |probability                               |prediction|
//    +-----+--------------+--------------------------------------+------------------------------------------+----------+
//    |1.0  |[-1.0,1.5,1.3]|[-6.587201443935503,6.587201443935503]|[0.0013759947069214356,0.9986240052930786]|1.0       |
//    |0.0  |[3.0,2.0,-0.1]|[3.980182819425659,-3.980182819425659]|[0.9816604009374171,0.01833959906258293]  |0.0       |
//    |1.0  |[0.0,2.2,-1.5]|[-6.376517702860472,6.376517702860472]|[0.0016981475578358176,0.9983018524421641]|1.0       |
//    +-----+--------------+--------------------------------------+------------------------------------------+----------+

val data = Array(
  Vectors.dense(-1.0, 1.5, 1.3),
  Vectors.dense(3.0, 2.0, -0.1),
  Vectors.dense(0.0, 2.2, -1.5)
)
val df = spark.createDataFrame(data.map(Tuple1.apply)).toDF("features")
binomialModel.transform(df).show(false)
//    +--------------+--------------------------------------+------------------------------------------+----------+
//    |features      |rawPrediction                         |probability                               |prediction|
//    +--------------+--------------------------------------+------------------------------------------+----------+
//    |[-1.0,1.5,1.3]|[-6.587201443935503,6.587201443935503]|[0.0013759947069214356,0.9986240052930786]|1.0       |
//    |[3.0,2.0,-0.1]|[3.980182819425659,-3.980182819425659]|[0.9816604009374171,0.01833959906258293]  |0.0       |
//    |[0.0,2.2,-1.5]|[-6.376517702860472,6.376517702860472]|[0.0016981475578358176,0.9983018524421641]|1.0       |
//    +--------------+--------------------------------------+------------------------------------------+----------+

val multinomialTraining = spark.createDataFrame(Seq(
  (3.0, Vectors.dense(0.0, 1.1, 0.1)),
  (0.0, Vectors.dense(2.0, 1.0, -1.0)),
  (0.0, Vectors.dense(2.0, 1.3, 1.0)),
  (1.0, Vectors.dense(0.0, 1.2, -0.5))
)).toDF("label", "features")
multinomialTraining.show(false)
//    +-----+--------------+
//    |label|features      |
//    +-----+--------------+
//    |3.0  |[0.0,1.1,0.1] |
//    |0.0  |[2.0,1.0,-1.0]|
//    |0.0  |[2.0,1.3,1.0] |
//    |1.0  |[0.0,1.2,-0.5]|
//    +-----+--------------+

val multinomialModel = lr.fit(multinomialTraining)
println("Multinomial model's numClasses: " + multinomialModel.numClasses)
//4
println("Multinomial model's numFeatures: " + multinomialModel.numFeatures)
//3
println("Multinomial model's coefficients: " + multinomialModel.coefficientMatrix)
//    2.3755319833863706    4.288895693579067   0.3392571745662531
//    -1.116821154666315    6.986561348601089   -1.1694540068639707
//    -0.19800783089303808  -18.26543596057775  0.04667029423568621
//    -1.0607029978270186   6.989978918397579   0.7835265380620305
println("Multinomial model's intercept: " + multinomialModel.interceptVector)
//[0.5485507718070327,0.05028251887521283,-0.9026627343287896,0.30382944364654396]
println("Multinomial model's hasSummary: " + multinomialModel.hasSummary)
//false

multinomialModel.transform(df).show(false)
//    +--------------+-----------------------------------------------------------------------------+-------------------------------------------------------------------------------------+----------+
//    |features      |rawPrediction                                                                |probability                                                                          |prediction|
//    +--------------+-----------------------------------------------------------------------------+-------------------------------------------------------------------------------------+----------+
//    |[-1.0,1.5,1.3]|[5.047396655725392,10.126655487519999,-28.04213746179599,12.868085318550571] |[3.768926166805804E-4,0.06054966051507028,1.6055193854749213E-18,0.9390734468682492] |3.0       |
//    |[3.0,2.0,-0.1]|[16.219012391667654,10.789887152764843,-38.03222517758697,11.023325633154444]|[0.9901702416440773,0.004343809163875811,2.7208008666051956E-24,0.005485949192046843]|0.0       |
//    |[0.0,2.2,-1.5]|[9.475235535831603,17.17489849609357,-41.15662728895338,14.506493257028174]  |[4.2341855234635453E-4,0.9347404120410062,4.341414472675871E-26,0.06483616940664733] |1.0       |
//    +--------------+-----------------------------------------------------------------------------+-------------------------------------------------------------------------------------+----------+