1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
|
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.spark.ml.classification
import scala.collection.mutable
import breeze.linalg.{DenseVector => BDV, norm => brzNorm}
import breeze.optimize.{CachedDiffFunction, DiffFunction, LBFGS => BreezeLBFGS, OWLQN => BreezeOWLQN}
import org.apache.spark.{Logging, SparkException}
import org.apache.spark.annotation.Experimental
import org.apache.spark.ml.param._
import org.apache.spark.ml.param.shared._
import org.apache.spark.ml.util.Identifiable
import org.apache.spark.mllib.linalg._
import org.apache.spark.mllib.linalg.BLAS._
import org.apache.spark.mllib.regression.LabeledPoint
import org.apache.spark.mllib.stat.MultivariateOnlineSummarizer
import org.apache.spark.mllib.util.MLUtils
import org.apache.spark.rdd.RDD
import org.apache.spark.sql.DataFrame
import org.apache.spark.storage.StorageLevel
/**
* Params for logistic regression.
*/
private[classification] trait LogisticRegressionParams extends ProbabilisticClassifierParams
with HasRegParam with HasElasticNetParam with HasMaxIter with HasFitIntercept with HasTol
with HasThreshold
/**
* :: Experimental ::
* Logistic regression.
* Currently, this class only supports binary classification.
*/
@Experimental
class LogisticRegression(override val uid: String)
extends ProbabilisticClassifier[Vector, LogisticRegression, LogisticRegressionModel]
with LogisticRegressionParams with Logging {
def this() = this(Identifiable.randomUID("logreg"))
/**
* Set the regularization parameter.
* Default is 0.0.
* @group setParam
*/
def setRegParam(value: Double): this.type = set(regParam, value)
setDefault(regParam -> 0.0)
/**
* Set the ElasticNet mixing parameter.
* For alpha = 0, the penalty is an L2 penalty. For alpha = 1, it is an L1 penalty.
* For 0 < alpha < 1, the penalty is a combination of L1 and L2.
* Default is 0.0 which is an L2 penalty.
* @group setParam
*/
def setElasticNetParam(value: Double): this.type = set(elasticNetParam, value)
setDefault(elasticNetParam -> 0.0)
/**
* Set the maximal number of iterations.
* Default is 100.
* @group setParam
*/
def setMaxIter(value: Int): this.type = set(maxIter, value)
setDefault(maxIter -> 100)
/**
* Set the convergence tolerance of iterations.
* Smaller value will lead to higher accuracy with the cost of more iterations.
* Default is 1E-6.
* @group setParam
*/
def setTol(value: Double): this.type = set(tol, value)
setDefault(tol -> 1E-6)
/** @group setParam */
def setFitIntercept(value: Boolean): this.type = set(fitIntercept, value)
setDefault(fitIntercept -> true)
/** @group setParam */
def setThreshold(value: Double): this.type = set(threshold, value)
setDefault(threshold -> 0.5)
override protected def train(dataset: DataFrame): LogisticRegressionModel = {
// Extract columns from data. If dataset is persisted, do not persist oldDataset.
val instances = extractLabeledPoints(dataset).map {
case LabeledPoint(label: Double, features: Vector) => (label, features)
}
val handlePersistence = dataset.rdd.getStorageLevel == StorageLevel.NONE
if (handlePersistence) instances.persist(StorageLevel.MEMORY_AND_DISK)
val (summarizer, labelSummarizer) = instances.treeAggregate(
(new MultivariateOnlineSummarizer, new MultiClassSummarizer))(
seqOp = (c, v) => (c, v) match {
case ((summarizer: MultivariateOnlineSummarizer, labelSummarizer: MultiClassSummarizer),
(label: Double, features: Vector)) =>
(summarizer.add(features), labelSummarizer.add(label))
},
combOp = (c1, c2) => (c1, c2) match {
case ((summarizer1: MultivariateOnlineSummarizer,
classSummarizer1: MultiClassSummarizer), (summarizer2: MultivariateOnlineSummarizer,
classSummarizer2: MultiClassSummarizer)) =>
(summarizer1.merge(summarizer2), classSummarizer1.merge(classSummarizer2))
})
val histogram = labelSummarizer.histogram
val numInvalid = labelSummarizer.countInvalid
val numClasses = histogram.length
val numFeatures = summarizer.mean.size
if (numInvalid != 0) {
val msg = s"Classification labels should be in {0 to ${numClasses - 1} " +
s"Found $numInvalid invalid labels."
logError(msg)
throw new SparkException(msg)
}
if (numClasses > 2) {
val msg = s"Currently, LogisticRegression with ElasticNet in ML package only supports " +
s"binary classification. Found $numClasses in the input dataset."
logError(msg)
throw new SparkException(msg)
}
val featuresMean = summarizer.mean.toArray
val featuresStd = summarizer.variance.toArray.map(math.sqrt)
val regParamL1 = $(elasticNetParam) * $(regParam)
val regParamL2 = (1.0 - $(elasticNetParam)) * $(regParam)
val costFun = new LogisticCostFun(instances, numClasses, $(fitIntercept),
featuresStd, featuresMean, regParamL2)
val optimizer = if ($(elasticNetParam) == 0.0 || $(regParam) == 0.0) {
new BreezeLBFGS[BDV[Double]]($(maxIter), 10, $(tol))
} else {
// Remove the L1 penalization on the intercept
def regParamL1Fun = (index: Int) => {
if (index == numFeatures) 0.0 else regParamL1
}
new BreezeOWLQN[Int, BDV[Double]]($(maxIter), 10, regParamL1Fun, $(tol))
}
val initialWeightsWithIntercept =
Vectors.zeros(if ($(fitIntercept)) numFeatures + 1 else numFeatures)
if ($(fitIntercept)) {
/**
* For binary logistic regression, when we initialize the weights as zeros,
* it will converge faster if we initialize the intercept such that
* it follows the distribution of the labels.
*
* {{{
* P(0) = 1 / (1 + \exp(b)), and
* P(1) = \exp(b) / (1 + \exp(b))
* }}}, hence
* {{{
* b = \log{P(1) / P(0)} = \log{count_1 / count_0}
* }}}
*/
initialWeightsWithIntercept.toArray(numFeatures)
= math.log(histogram(1).toDouble / histogram(0).toDouble)
}
val states = optimizer.iterations(new CachedDiffFunction(costFun),
initialWeightsWithIntercept.toBreeze.toDenseVector)
var state = states.next()
val lossHistory = mutable.ArrayBuilder.make[Double]
while (states.hasNext) {
lossHistory += state.value
state = states.next()
}
lossHistory += state.value
// The weights are trained in the scaled space; we're converting them back to
// the original space.
val weightsWithIntercept = {
val rawWeights = state.x.toArray.clone()
var i = 0
// Note that the intercept in scaled space and original space is the same;
// as a result, no scaling is needed.
while (i < numFeatures) {
rawWeights(i) *= { if (featuresStd(i) != 0.0) 1.0 / featuresStd(i) else 0.0 }
i += 1
}
Vectors.dense(rawWeights)
}
if (handlePersistence) instances.unpersist()
val (weights, intercept) = if ($(fitIntercept)) {
(Vectors.dense(weightsWithIntercept.toArray.slice(0, weightsWithIntercept.size - 1)),
weightsWithIntercept(weightsWithIntercept.size - 1))
} else {
(weightsWithIntercept, 0.0)
}
new LogisticRegressionModel(uid, weights.compressed, intercept)
}
}
/**
* :: Experimental ::
* Model produced by [[LogisticRegression]].
*/
@Experimental
class LogisticRegressionModel private[ml] (
override val uid: String,
val weights: Vector,
val intercept: Double)
extends ProbabilisticClassificationModel[Vector, LogisticRegressionModel]
with LogisticRegressionParams {
/** @group setParam */
def setThreshold(value: Double): this.type = set(threshold, value)
/** Margin (rawPrediction) for class label 1. For binary classification only. */
private val margin: Vector => Double = (features) => {
BLAS.dot(features, weights) + intercept
}
/** Score (probability) for class label 1. For binary classification only. */
private val score: Vector => Double = (features) => {
val m = margin(features)
1.0 / (1.0 + math.exp(-m))
}
override val numClasses: Int = 2
/**
* Predict label for the given feature vector.
* The behavior of this can be adjusted using [[threshold]].
*/
override protected def predict(features: Vector): Double = {
if (score(features) > getThreshold) 1 else 0
}
override protected def raw2probabilityInPlace(rawPrediction: Vector): Vector = {
rawPrediction match {
case dv: DenseVector =>
var i = 0
val size = dv.size
while (i < size) {
dv.values(i) = 1.0 / (1.0 + math.exp(-dv.values(i)))
i += 1
}
dv
case sv: SparseVector =>
throw new RuntimeException("Unexpected error in LogisticRegressionModel:" +
" raw2probabilitiesInPlace encountered SparseVector")
}
}
override protected def predictRaw(features: Vector): Vector = {
val m = margin(features)
Vectors.dense(-m, m)
}
override def copy(extra: ParamMap): LogisticRegressionModel = {
copyValues(new LogisticRegressionModel(uid, weights, intercept), extra)
}
override protected def raw2prediction(rawPrediction: Vector): Double = {
val t = getThreshold
val rawThreshold = if (t == 0.0) {
Double.NegativeInfinity
} else if (t == 1.0) {
Double.PositiveInfinity
} else {
math.log(t / (1.0 - t))
}
if (rawPrediction(1) > rawThreshold) 1 else 0
}
override protected def probability2prediction(probability: Vector): Double = {
if (probability(1) > getThreshold) 1 else 0
}
}
/**
* MultiClassSummarizer computes the number of distinct labels and corresponding counts,
* and validates the data to see if the labels used for k class multi-label classification
* are in the range of {0, 1, ..., k - 1} in a online fashion.
*
* Two MultilabelSummarizer can be merged together to have a statistical summary of the
* corresponding joint dataset.
*/
private[classification] class MultiClassSummarizer extends Serializable {
private val distinctMap = new mutable.HashMap[Int, Long]
private var totalInvalidCnt: Long = 0L
/**
* Add a new label into this MultilabelSummarizer, and update the distinct map.
* @param label The label for this data point.
* @return This MultilabelSummarizer
*/
def add(label: Double): this.type = {
if (label - label.toInt != 0.0 || label < 0) {
totalInvalidCnt += 1
this
}
else {
val counts: Long = distinctMap.getOrElse(label.toInt, 0L)
distinctMap.put(label.toInt, counts + 1)
this
}
}
/**
* Merge another MultilabelSummarizer, and update the distinct map.
* (Note that it will merge the smaller distinct map into the larger one using in-place
* merging, so either `this` or `other` object will be modified and returned.)
*
* @param other The other MultilabelSummarizer to be merged.
* @return Merged MultilabelSummarizer object.
*/
def merge(other: MultiClassSummarizer): MultiClassSummarizer = {
val (largeMap, smallMap) = if (this.distinctMap.size > other.distinctMap.size) {
(this, other)
} else {
(other, this)
}
smallMap.distinctMap.foreach {
case (key, value) =>
val counts = largeMap.distinctMap.getOrElse(key, 0L)
largeMap.distinctMap.put(key, counts + value)
}
largeMap.totalInvalidCnt += smallMap.totalInvalidCnt
largeMap
}
/** @return The total invalid input counts. */
def countInvalid: Long = totalInvalidCnt
/** @return The number of distinct labels in the input dataset. */
def numClasses: Int = distinctMap.keySet.max + 1
/** @return The counts of each label in the input dataset. */
def histogram: Array[Long] = {
val result = Array.ofDim[Long](numClasses)
var i = 0
val len = result.length
while (i < len) {
result(i) = distinctMap.getOrElse(i, 0L)
i += 1
}
result
}
}
/**
* LogisticAggregator computes the gradient and loss for binary logistic loss function, as used
* in binary classification for samples in sparse or dense vector in a online fashion.
*
* Note that multinomial logistic loss is not supported yet!
*
* Two LogisticAggregator can be merged together to have a summary of loss and gradient of
* the corresponding joint dataset.
*
* @param weights The weights/coefficients corresponding to the features.
* @param numClasses the number of possible outcomes for k classes classification problem in
* Multinomial Logistic Regression.
* @param fitIntercept Whether to fit an intercept term.
* @param featuresStd The standard deviation values of the features.
* @param featuresMean The mean values of the features.
*/
private class LogisticAggregator(
weights: Vector,
numClasses: Int,
fitIntercept: Boolean,
featuresStd: Array[Double],
featuresMean: Array[Double]) extends Serializable {
private var totalCnt: Long = 0L
private var lossSum = 0.0
private val weightsArray = weights match {
case dv: DenseVector => dv.values
case _ =>
throw new IllegalArgumentException(
s"weights only supports dense vector but got type ${weights.getClass}.")
}
private val dim = if (fitIntercept) weightsArray.length - 1 else weightsArray.length
private val gradientSumArray = Array.ofDim[Double](weightsArray.length)
/**
* Add a new training data to this LogisticAggregator, and update the loss and gradient
* of the objective function.
*
* @param label The label for this data point.
* @param data The features for one data point in dense/sparse vector format to be added
* into this aggregator.
* @return This LogisticAggregator object.
*/
def add(label: Double, data: Vector): this.type = {
require(dim == data.size, s"Dimensions mismatch when adding new sample." +
s" Expecting $dim but got ${data.size}.")
val dataSize = data.size
val localWeightsArray = weightsArray
val localGradientSumArray = gradientSumArray
numClasses match {
case 2 =>
/**
* For Binary Logistic Regression.
*/
val margin = - {
var sum = 0.0
data.foreachActive { (index, value) =>
if (featuresStd(index) != 0.0 && value != 0.0) {
sum += localWeightsArray(index) * (value / featuresStd(index))
}
}
sum + { if (fitIntercept) localWeightsArray(dim) else 0.0 }
}
val multiplier = (1.0 / (1.0 + math.exp(margin))) - label
data.foreachActive { (index, value) =>
if (featuresStd(index) != 0.0 && value != 0.0) {
localGradientSumArray(index) += multiplier * (value / featuresStd(index))
}
}
if (fitIntercept) {
localGradientSumArray(dim) += multiplier
}
if (label > 0) {
// The following is equivalent to log(1 + exp(margin)) but more numerically stable.
lossSum += MLUtils.log1pExp(margin)
} else {
lossSum += MLUtils.log1pExp(margin) - margin
}
case _ =>
new NotImplementedError("LogisticRegression with ElasticNet in ML package only supports " +
"binary classification for now.")
}
totalCnt += 1
this
}
/**
* Merge another LogisticAggregator, and update the loss and gradient
* of the objective function.
* (Note that it's in place merging; as a result, `this` object will be modified.)
*
* @param other The other LogisticAggregator to be merged.
* @return This LogisticAggregator object.
*/
def merge(other: LogisticAggregator): this.type = {
require(dim == other.dim, s"Dimensions mismatch when merging with another " +
s"LeastSquaresAggregator. Expecting $dim but got ${other.dim}.")
if (other.totalCnt != 0) {
totalCnt += other.totalCnt
lossSum += other.lossSum
var i = 0
val localThisGradientSumArray = this.gradientSumArray
val localOtherGradientSumArray = other.gradientSumArray
val len = localThisGradientSumArray.length
while (i < len) {
localThisGradientSumArray(i) += localOtherGradientSumArray(i)
i += 1
}
}
this
}
def count: Long = totalCnt
def loss: Double = lossSum / totalCnt
def gradient: Vector = {
val result = Vectors.dense(gradientSumArray.clone())
scal(1.0 / totalCnt, result)
result
}
}
/**
* LogisticCostFun implements Breeze's DiffFunction[T] for a multinomial logistic loss function,
* as used in multi-class classification (it is also used in binary logistic regression).
* It returns the loss and gradient with L2 regularization at a particular point (weights).
* It's used in Breeze's convex optimization routines.
*/
private class LogisticCostFun(
data: RDD[(Double, Vector)],
numClasses: Int,
fitIntercept: Boolean,
featuresStd: Array[Double],
featuresMean: Array[Double],
regParamL2: Double) extends DiffFunction[BDV[Double]] {
override def calculate(weights: BDV[Double]): (Double, BDV[Double]) = {
val w = Vectors.fromBreeze(weights)
val logisticAggregator = data.treeAggregate(new LogisticAggregator(w, numClasses, fitIntercept,
featuresStd, featuresMean))(
seqOp = (c, v) => (c, v) match {
case (aggregator, (label, features)) => aggregator.add(label, features)
},
combOp = (c1, c2) => (c1, c2) match {
case (aggregator1, aggregator2) => aggregator1.merge(aggregator2)
})
// regVal is the sum of weight squares for L2 regularization
val norm = if (regParamL2 == 0.0) {
0.0
} else if (fitIntercept) {
brzNorm(Vectors.dense(weights.toArray.slice(0, weights.size -1)).toBreeze, 2.0)
} else {
brzNorm(weights, 2.0)
}
val regVal = 0.5 * regParamL2 * norm * norm
val loss = logisticAggregator.loss + regVal
val gradient = logisticAggregator.gradient
if (fitIntercept) {
val wArray = w.toArray.clone()
wArray(wArray.length - 1) = 0.0
axpy(regParamL2, Vectors.dense(wArray), gradient)
} else {
axpy(regParamL2, w, gradient)
}
(loss, gradient.toBreeze.asInstanceOf[BDV[Double]])
}
}
|