Fix lbfgsb linesearch out of bound and findAlpha method

math/src/main/scala/breeze/optimize/LBFGSB.scala

@@ -77,14 +77,47 @@ class LBFGSB(lowerBounds: DenseVector[Double],
 
 
   override protected def determineStepSize(state: State, f: DiffFunction[DenseVector[Double]], direction: DenseVector[Double]): Double =  {
-    val x = state.x
-    val ff = LineSearch.functionFromSearchDirection(f, x, direction)
+    val ff = new DiffFunction[Double] {
+      def calculate(alpha: Double) = {
+        val newX = takeStep(state, direction, alpha)
+        val (ff, grad) = f.calculate(newX)
+        ff -> (grad dot direction)
+      }
+    }
     val wolfeRuleSearch = new StrongWolfeLineSearch(maxZoomIter, maxLineSearchIter) // TODO: Need good default values here.
-    wolfeRuleSearch.minimize(ff, 1.0)
+
+    var minStepBound = Double.PositiveInfinity
+    var i = 0
+    while (i < lowerBounds.length) {
+      val dir = direction(i)
+      if (dir != 0.0) {
+        val bound = if (dir < 0.0) lowerBounds(i) else upperBounds(i)
+        val stepBound = (bound - state.x(i)) / dir
+        assert(stepBound >= 0.0)
+        if (stepBound < minStepBound) {
+          minStepBound = stepBound
+        }
+      }
+      i += 1
+    }
+
+    val initStep = if (minStepBound < 1.0) minStepBound else 1.0
+    wolfeRuleSearch.minimizeWithBound(ff, initStep, minStepBound)
   }
 
   override protected def takeStep(state: State, dir: DenseVector[Double], stepSize: Double) = {
-    state.x + (dir *:* stepSize)
+    val newX = state.x + (dir :* stepSize)
+    var i = 0
+    while (i < newX.length) {
+      if (newX(i) > upperBounds(i)) {
+        newX(i) = upperBounds(i)
+      }
+      if (newX(i) < lowerBounds(i)) {
+        newX(i) = lowerBounds(i)
+      }
+      i += 1
+    }
+    newX
   }
 
   private def initialize(f: DiffFunction[DenseVector[Double]], x0: DenseVector[Double]) = {

math/src/main/scala/breeze/optimize/StrongWolfe.scala

@@ -58,12 +58,19 @@ class StrongWolfeLineSearch(maxZoomIter: Int, maxLineSearchIter: Int) extends Cu
   val c1 = 1e-4
   val c2 = 0.9
 
+  def minimize(f: DiffFunction[Double], init: Double = 1.0): Double = {
+    minimizeWithBound(f, init = 1.0, bound = Double.PositiveInfinity)
+  }
+
   /**
-   * Performs a line search on the function f, returning a point satisfying
-   * the Strong Wolfe conditions. Based on the line search detailed in
-   * Nocedal & Wright Numerical Optimization p58.
+   * Performs a line search on the function f with bound, returning a point satisfying
+   * the Strong Wolfe conditions OR satisfying sufficient decrease condition and hit bound.
+   * Based on the line search detailed in Nocedal & Wright Numerical Optimization p58.
+   * BUT add some modification for bound checking.
    */
-  def minimize(f: DiffFunction[Double], init: Double = 1.0):Double = {
+  def minimizeWithBound(f: DiffFunction[Double], init: Double = 1.0, bound: Double = 1.0): Double = {
+
+    require(init <= bound, "init value should <= bound")
 
     def phi(t: Double): Bracket = {
       val (pval, pdd) = f.calculate(t)
@@ -171,8 +178,17 @@ class StrongWolfeLineSearch(maxZoomIter: Int, maxLineSearchIter: Int) extends Cu
         }
 
         low = c
-        t *= 1.5
-        logger.debug("Sufficent Decrease condition but not curvature condition satisfied. Increased t to: " + t)
+        if (t == bound) {
+          logger.debug("Reach bound, satisfy sufficent decrease condition," +
+            " but not curvature condition satisfied.")
+          return bound
+        } else {
+          t *= 1.5
+          if (t > bound) {
+            t = bound
+          }
+          logger.debug("Sufficent Decrease condition but not curvature condition satisfied. Increased t to: " + t)
+        }
       }
     }
 

math/src/test/scala/breeze/optimize/LBFGSBTest.scala

@@ -97,4 +97,23 @@ class LBFGSBTest extends OptimizeTestBase{
     assert(ures.value < res.value)
   }
 
+  test("issue 572") {
+    val solver = new LBFGSB(DenseVector[Double](1E-12), DenseVector[Double](Double.MaxValue))
+
+    val f = new DiffFunction[DenseVector[Double]] {
+      override def calculate(x: DenseVector[Double]): (Double, DenseVector[Double]) = {
+        val cost = x(0) + 1.0/x(0)
+        val grad = DenseVector(1.0 - 1.0/(x(0)*x(0)))
+        (cost, grad)
+      }
+    }
+
+    val nearX0 = DenseVector[Double](1.5)
+    val nearRes = solver.minimizeAndReturnState(f, nearX0)
+    assert(abs(nearRes.x(0) - 1.0) < EPS)
+
+    val farX0 = DenseVector[Double](1500)
+    val farRes = solver.minimizeAndReturnState(f, farX0)
+    assert(abs(farRes.x(0) - 1.0) < EPS)
+  }
 }