merge with qp branch; NOTICE file updated

NOTICE

@@ -8,6 +8,10 @@ Copyright (c) 2014 by Debasish Das (Verizon), all rights reserved.
 
 ===============================================================================
 
+QuadraticMinimizer class in package breeze.optimize.proximal is distributed with Copyright (c)
+2014, Debasish Das (Verizon), all rights reserved.
+
+===============================================================================
+
 NonlinearMinimizer class in package breeze.optimize.quadratic is distributed with Copyright (c)
 2015, Debasish Das (Verizon), all rights reserved.
-

math/src/main/scala/breeze/optimize/quadratic/NNLS.scala → math/src/main/scala/breeze/optimize/linear/NNLS.scala

@@ -15,12 +15,12 @@
  * limitations under the License.
  */
 
-package breeze.optimize.quadratic
+package breeze.optimize.linear
 
-import breeze.linalg.DenseVector
-import breeze.linalg.DenseMatrix
-import breeze.linalg.axpy
+import breeze.linalg.{DenseMatrix, DenseVector, axpy}
+import breeze.optimize.proximal.QpGenerator
 import breeze.stats.distributions.Rand
+import breeze.util.Implicits._
 /**
  * Object used to solve nonnegative least squares problems using a modified
  * projected gradient method.
@@ -30,13 +30,14 @@ class NNLS(val maxIters: Int = -1) {
   type BDM = DenseMatrix[Double]
   type BDV = DenseVector[Double]
 
-  case class State(x: BDV, grad: BDV, res: BDV, iterations: Int, converged: Boolean)
+  case class State(x: BDV, grad: BDV, dir: BDV, lastDir: BDV, res: BDV,
+                   lastNorm: Double, lastWall: Int, iter: Int, converged: Boolean)
 
   // find the optimal unconstrained step
   def steplen(ata: BDM, dir: BDV, res: BDV,
               tmp: BDV): Double = {
     val top = dir.dot(res)
-    tmp := ata*dir
+    tmp := ata * dir
     // Push the denominator upward very slightly to avoid infinities and silliness
     top / (tmp.dot(dir) + 1e-20)
   }
@@ -56,47 +57,45 @@ class NNLS(val maxIters: Int = -1) {
    * projected gradient method.  That is, find x minimising ||Ax - b||_2 given A^T A and A^T b.
    *
    * We solve the problem
-   *   min_x      1/2 x' ata x' - x'atb
-   *   subject to x >= 0
+   * min_x      1/2 x' ata x' - x'atb
+   * subject to x >= 0
    *
    * The method used is similar to one described by Polyak (B. T. Polyak, The conjugate gradient
    * method in extremal problems, Zh. Vychisl. Mat. Mat. Fiz. 9(4)(1969), pp. 94-112) for bound-
    * constrained nonlinear programming.  Polyak unconditionally uses a conjugate gradient
    * direction, however, while this method only uses a conjugate gradient direction if the last
    * iteration did not cause a previously-inactive constraint to become active.
    */
-  def minimize(ata: DenseMatrix[Double], atb: DenseVector[Double]) : DenseVector[Double] = {
-    iterations(ata, atb).x
-  }
-
-  def iterations(ata: DenseMatrix[Double], atb: DenseVector[Double]) : State = {
+  def initialState(ata: DenseMatrix[Double], atb: DenseVector[Double], n: Int): State = {
     require(ata.cols == ata.rows, s"NNLS:iterations gram matrix must be symmetric")
     require(ata.rows == atb.length, s"NNLS:iterations gram matrix rows must be same as length of linear term")
 
-    val n = atb.length
-
-    val tmp = DenseVector.zeros[Double](n)
     val grad = DenseVector.zeros[Double](n)
-    val x =  DenseVector.zeros[Double](n)
+    val x = DenseVector.zeros[Double](n)
     val dir = DenseVector.zeros[Double](n)
     val lastDir = DenseVector.zeros[Double](n)
     val res = DenseVector.zeros[Double](n)
+    val lastNorm = 0.0
+    val lastWall = 0
+    State(x, grad, dir, lastDir, res, lastNorm, lastWall, 0, false)
+  }
 
-    val iterMax = if (maxIters < 0) Math.max(400, 20 * n) else maxIters
+  def iterations(ata: DenseMatrix[Double],
+                 atb: DenseVector[Double]): Iterator[State] =
+    Iterator.iterate(initialState(ata, atb, atb.length)) { state =>
+      import state._
+      val n = atb.length
+      val tmp = DenseVector.zeros[Double](atb.length)
 
-    var lastNorm = 0.0
-    var iterno = 0
-    var lastWall = 0 // Last iteration when we hit a bound constraint.
-    var i = 0
+      val iterMax = if (maxIters < 0) Math.max(400, 20 * n) else maxIters
 
-    while (iterno < iterMax) {
       // find the residual
-      res := ata*x
+      res := ata * x
       res -= atb
       grad := res
-      
+
       // project the gradient
-      i = 0
+      var i = 0
       while (i < n) {
         if (grad.data(i) > 0.0 && x.data(i) == 0.0) {
           grad.data(i) = 0.0
@@ -105,13 +104,13 @@ class NNLS(val maxIters: Int = -1) {
       }
       val ngrad = grad.dot(grad)
       dir := grad
-          
+
       // use a CG direction under certain conditions
       var step = steplen(ata, grad, res, tmp)
       var ndir = 0.0
       val nx = x.dot(x)
-      
-      if (iterno > lastWall + 1) {
+
+      if (iter > lastWall + 1) {
         val alpha = ngrad / lastNorm
         axpy(alpha, lastDir, dir)
         val dstep = steplen(ata, dir, res, tmp)
@@ -128,36 +127,44 @@ class NNLS(val maxIters: Int = -1) {
       }
 
       // terminate?
-      if (stop(step, ndir, nx)) {
-        return State(x, grad, res, iterno, true)
-      }
-
-      // don't run through the walls
-      i = 0
-      while (i < n) {
-        if (step * dir.data(i) > x.data(i)) {
-          step = x.data(i) / dir.data(i)
+      if (stop(step, ndir, nx) || iter > iterMax)
+        State(x, grad, dir, lastDir, res, lastNorm, lastWall, iter + 1, true)
+      else {
+        // don't run through the walls
+        i = 0
+        while (i < n) {
+          if (step * dir.data(i) > x.data(i)) {
+            step = x.data(i) / dir.data(i)
+          }
+          i = i + 1
         }
-        i = i + 1
-      }
 
-      // take the step
-      i = 0
-      while (i < n) {
-        if (step * dir.data(i) > x.data(i) * (1 - 1e-14)) {
-          x.data(i) = 0
-          lastWall = iterno
-        } else {
-          x.data(i) -= step * dir.data(i)
+        var nextWall = lastWall
+
+        // take the step
+        i = 0
+        while (i < n) {
+          if (step * dir.data(i) > x.data(i) * (1 - 1e-14)) {
+            x.data(i) = 0
+            nextWall = iter
+          } else {
+            x.data(i) -= step * dir.data(i)
+          }
+          i = i + 1
         }
-        i = i + 1
+        lastDir := dir
+        val nextNorm = ngrad
+        State(x, grad, dir, lastDir, res, nextNorm, nextWall, iter + 1, false)
       }
+    }.takeUpToWhere(_.converged)
 
-      iterno = iterno + 1
-      lastDir := dir
-      lastNorm = ngrad
-    }
-    State(x, grad, res, iterno, false)
+  def minimizeAndReturnState(ata: DenseMatrix[Double],
+                             atb: DenseVector[Double]): State = {
+    iterations(ata, atb).last
+  }
+
+  def minimize(ata: DenseMatrix[Double], atb: DenseVector[Double]): DenseVector[Double] = {
+    minimizeAndReturnState(ata, atb).x
   }
 }
 

math/src/main/scala/breeze/optimize/quadratic/Constraint.scala → math/src/main/scala/breeze/optimize/proximal/Constraint.scala

@@ -14,9 +14,9 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package breeze.optimize.quadratic
+package breeze.optimize.proximal
 
 object Constraint extends Enumeration {
   type Constraint = Value
-  val SMOOTH, POSITIVE, BOUNDS, SPARSE, EQUALITY = Value
+  val SMOOTH, POSITIVE, BOX, SPARSE, EQUALITY = Value
 }

math/src/main/scala/breeze/optimize/quadratic/NonlinearMinimizer.scala → math/src/main/scala/breeze/optimize/proximal/NonlinearMinimizer.scala

@@ -15,12 +15,12 @@
  * limitations under the License.
  */
 
-package breeze.optimize.quadratic
+package breeze.optimize.proximal
 
-import breeze.linalg.{norm, DenseVector}
-import breeze.optimize.quadratic.Constraint._
-import breeze.optimize.{LBFGS, DiffFunction}
+import breeze.linalg.{DenseVector, norm}
+import breeze.optimize.{DiffFunction, LBFGS}
 import breeze.stats.distributions.Rand
+
 import scala.math._
 
 /**

math/src/main/scala/breeze/optimize/quadratic/Proximal.scala → math/src/main/scala/breeze/optimize/proximal/Proximal.scala

@@ -19,7 +19,7 @@
  * Library of Proximal Algorithms adapted from https://github.com/cvxgrp/proximal
  */
 
-package breeze.optimize.quadratic
+package breeze.optimize.proximal
 
 import scala.math.max
 import scala.math.min
@@ -122,7 +122,7 @@ case class ProximalL1() extends Proximal {
   def prox(x: DenseVector[Double], rho: Double) = {
     var i = 0
     while (i < x.length) {
-      x.update(i, max(0, x(i) - rho) - max(0, -x(i) - rho))
+      x.update(i, max(0, x(i) - lambda/rho) - max(0, -x(i) - lambda/rho))
       i = i + 1
     }
   }
@@ -243,4 +243,4 @@ case class ProximalLp(c: DenseVector[Double]) extends Proximal {
       i = i + 1
     }
   }
-}
+}