深入机器学习系列之BFGS & L-BFGS

1val optimizer = new BreezeLBFGS[BDV[Double]]($(maxIter), 10, $(tol))
2val states =
3      optimizer.iterations(new CachedDiffFunction(costFun), initialWeights.toBreeze.toDenseVector)

 1def iterations(f: DF, init: T): Iterator[State] = {
 2    val adjustedFun = adjustFunction(f)
 3    infiniteIterations(f, initialState(adjustedFun, init)).takeUpToWhere(_.converged)
 4}
 5//调用infiniteIterations，其中State是一个样本类
 6def infiniteIterations(f: DF, state: State): Iterator[State] = {
 7    var failedOnce = false
 8    val adjustedFun = adjustFunction(f)
 9    //无限迭代
10    Iterator.iterate(state) { state => try {
11        //1 选择梯度下降方向
12        val dir = chooseDescentDirection(state, adjustedFun)
13        //2 计算步长
14        val stepSize = determineStepSize(state, adjustedFun, dir)
15        //3 更新权重
16        val x = takeStep(state,dir,stepSize)
17        //4 利用CostFun.calculate计算损失值和梯度
18        val (value,grad) = calculateObjective(adjustedFun, x, state.history)
19        val (adjValue,adjGrad) = adjust(x,grad,value)
20        val oneOffImprovement = (state.adjustedValue - adjValue)/(state.adjustedValue.abs max adjValue.abs max 1E-6 * state.initialAdjVal.abs)
21        //5 计算s和t
22        val history = updateHistory(x,grad,value, adjustedFun, state)
23        //6 只保存m个需要的s和t
24        val newAverage = updateFValWindow(state, adjValue)
25        failedOnce = false
26        var s = State(x,value,grad,adjValue,adjGrad,state.iter + 1, state.initialAdjVal, history, newAverage, 0)
27        val improvementFailure = (state.fVals.length >= minImprovementWindow && state.fVals.nonEmpty && state.fVals.last > state.fVals.head * (1-improvementTol))
28        if(improvementFailure)
29          s = s.copy(fVals = IndexedSeq.empty, numImprovementFailures = state.numImprovementFailures + 1)
30        s
31      } catch {
32        case x: FirstOrderException if !failedOnce =>
33          failedOnce = true
34          logger.error("Failure! Resetting history: " + x)
35          state.copy(history = initialHistory(adjustedFun, state.x))
36        case x: FirstOrderException =>
37          logger.error("Failure again! Giving up and returning. Maybe the objective is just poorly behaved?")
38          state.copy(searchFailed = true)
39      }
40    }
41  }

1protected def chooseDescentDirection(state: State, fn: DiffFunction[T]):T = {
2    state.history * state.grad
3}

 1def *(grad: T) = {
 2     val diag = if(historyLength > 0) {
 3       val prevStep = memStep.head
 4       val prevGradStep = memGradDelta.head
 5       val sy = prevStep dot prevGradStep
 6       val yy = prevGradStep dot prevGradStep
 7       if(sy < 0 || sy.isNaN) throw new NaNHistory
 8       sy/yy
 9     } else {
10       1.0
11     }
12     val dir = space.copy(grad)
13     val as = new Array[Double](m)
14     val rho = new Array[Double](m)
15     //第一次递归
16     for(i <- 0 until historyLength) {
17       rho(i) = (memStep(i) dot memGradDelta(i))
18       as(i) = (memStep(i) dot dir)/rho(i)
19       if(as(i).isNaN) {
20         throw new NaNHistory
21       }
22       axpy(-as(i), memGradDelta(i), dir)
23     }
24     dir *= diag
25     //第二次递归
26     for(i <- (historyLength - 1) to 0 by (-1)) {
27       val beta = (memGradDelta(i) dot dir)/rho(i)
28       axpy(as(i) - beta, memStep(i), dir)
29     }
30     dir *= -1.0
31     dir
32    }
33  }

 1protected def determineStepSize(state: State, f: DiffFunction[T], dir: T) = {
 2    val x = state.x
 3    val grad = state.grad
 4    val ff = LineSearch.functionFromSearchDirection(f, x, dir)
 5    val search = new StrongWolfeLineSearch(maxZoomIter = 10, maxLineSearchIter = 10) // TODO: Need good default values here.
 6    val alpha = search.minimize(ff, if(state.iter == 0.0) 1.0/norm(dir) else 1.0)
 7    if(alpha * norm(grad) < 1E-10)
 8      throw new StepSizeUnderflow
 9    alpha
10  }

1protected def takeStep(state: State, dir: T, stepSize: Double) = state.x + dir * stepSize

1protected def calculateObjective(f: DF, x: T, history: History): (Double, T) = {
2     f.calculate(x)
3  }

 1//计算s和t
 2protected def updateHistory(newX: T, newGrad: T, newVal: Double,  f: DiffFunction[T], oldState: State): History = {
 3    oldState.history.updated(newX - oldState.x, newGrad :- oldState.grad)
 4}
 5//添加新的s和t，并删除过期的s和t
 6protected def updateFValWindow(oldState: State, newAdjVal: Double):IndexedSeq[Double] = {
 7    val interm = oldState.fVals :+ newAdjVal
 8    if(interm.length > minImprovementWindow) interm.drop(1)
 9    else interm
10  }

 1override protected def chooseDescentDirection(state: State, fn: DiffFunction[T]) = {
 2    val descentDir = super.chooseDescentDirection(state.copy(grad = state.adjustedGradient), fn)
 3
 4    // The original paper requires that the descent direction be corrected to be
 5    // in the same directional (within the same hypercube) as the adjusted gradient for proof.
 6    // Although this doesn't seem to affect the outcome that much in most of cases, there are some cases
 7    // where the algorithm won't converge (confirmed with the author, Galen Andrew).
 8    val correctedDir = space.zipMapValues.map(descentDir, state.adjustedGradient, { case (d, g) => if (d * g < 0) d else 0.0 })
 9
10    correctedDir
11  }

 1override protected def determineStepSize(state: State, f: DiffFunction[T], dir: T) = {
 2    val iter = state.iter
 3
 4    val normGradInDir = {
 5      val possibleNorm = dir dot state.grad
 6      possibleNorm
 7    }
 8    val ff = new DiffFunction[Double] {
 9       def calculate(alpha: Double) = {
10         val newX = takeStep(state, dir, alpha)
11         val (v, newG) =  f.calculate(newX)  // 计算梯度
12         val (adjv, adjgrad) = adjust(newX, newG, v) // 调整梯度
13         adjv -> (adjgrad dot dir)
14       }
15    }
16    val search = new BacktrackingLineSearch(state.value, shrinkStep= if(iter < 1) 0.1 else 0.5)
17    val alpha = search.minimize(ff, if(iter < 1) .5/norm(state.grad) else 1.0)
18
19    alpha
20  }

 1// projects x to be on the same orthant as y
 2  // this basically requires that x'_i = x_i if sign(x_i) == sign(y_i), and 0 otherwise.
 3
 4  override protected def takeStep(state: State, dir: T, stepSize: Double) = {
 5    val stepped = state.x + dir * stepSize
 6    val orthant = computeOrthant(state.x, state.adjustedGradient)
 7    space.zipMapValues.map(stepped, orthant, { case (v, ov) =>
 8      v * I(math.signum(v) == math.signum(ov))
 9    })
10  }

 1// Adds in the regularization stuff to the gradient
 2  override protected def adjust(newX: T, newGrad: T, newVal: Double): (Double, T) = {
 3    var adjValue = newVal
 4    val res = space.zipMapKeyValues.mapActive(newX, newGrad, {case (i, xv, v) =>
 5      val l1regValue = l1reg(i)
 6      require(l1regValue >= 0.0)
 7
 8      if(l1regValue == 0.0) {
 9        v
10      } else {
11        adjValue += Math.abs(l1regValue * xv)
12        xv match {
13          case 0.0 => {
14            val delta_+ = v + l1regValue   //计算左导数
15            val delta_- = v - l1regValue   //计算右导数
16            if (delta_- > 0) delta_- else if (delta_+ < 0) delta_+ else 0.0
17          }
18          case _ => v + math.signum(xv) * l1regValue
19        }
20      }
21    })
22    adjValue -> res
23  }
24