NonlinearMinimizer using Projection and Proximal Operators

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

@@ -4,6 +4,7 @@ import breeze.linalg.{DenseMatrix, DenseVector, axpy}
 import breeze.optimize.proximal.QpGenerator
 import breeze.stats.distributions.Rand
 import breeze.util.Implicits._
+import spire.syntax.cfor._
 
 /**
  * NNLS solves nonnegative least squares problems using a modified
@@ -16,11 +17,10 @@ class NNLS(val maxIters: Int = -1) {
   type BDM = DenseMatrix[Double]
   type BDV = DenseVector[Double]
 
-  case class State(x: BDV, grad: BDV, dir: BDV, lastDir: BDV, res: BDV,
-                   lastNorm: Double, lastWall: Int, iter: Int, converged: Boolean)
+  case class State private[NNLS](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,
+  private def steplen(ata: BDM, dir: BDV, res: BDV,
               tmp: BDV): Double = {
     val top = dir.dot(res)
     tmp := ata * dir
@@ -29,7 +29,7 @@ class NNLS(val maxIters: Int = -1) {
   }
 
   // stopping condition
-  def stop(step: Double, ndir: Double, nx: Double): Boolean = {
+  private def stop(step: Double, ndir: Double, nx: Double): Boolean = {
     ((step.isNaN) // NaN
       || (step < 1e-7) // too small or negative
       || (step > 1e40) // too small; almost certainly numerical problems
@@ -52,7 +52,7 @@ class NNLS(val maxIters: Int = -1) {
    * 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 initialState(ata: DenseMatrix[Double], atb: DenseVector[Double], n: Int): State = {
+  private 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")
 
@@ -81,12 +81,10 @@ class NNLS(val maxIters: Int = -1) {
       grad := res
 
       // project the gradient
-      var i = 0
-      while (i < n) {
-        if (grad.data(i) > 0.0 && x.data(i) == 0.0) {
-          grad.data(i) = 0.0
+      cforRange(0 until n) { i =>
+        if (grad(i) > 0.0 && x(i) == 0.0) {
+          grad(i) = 0.0
         }
-        i = i + 1
       }
       val ngrad = grad.dot(grad)
       dir := grad
@@ -117,26 +115,22 @@ class NNLS(val maxIters: Int = -1) {
         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)
+        cforRange(0 until n){ i =>
+          if (step * dir(i) > x(i)) {
+            step = x(i) / dir(i)
           }
-          i = i + 1
         }
 
         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
+        cforRange(0 until n) { i =>
+          if (step * dir(i) > x(i) * (1 - 1e-14)) {
+            x(i) = 0
             nextWall = iter
           } else {
-            x.data(i) -= step * dir.data(i)
+            x(i) -= step * dir(i)
           }
-          i = i + 1
         }
         lastDir := dir
         val nextNorm = ngrad

math/src/main/scala/breeze/optimize/linear/PowerMethod.scala

@@ -8,16 +8,16 @@ import breeze.util.SerializableLogging
 import breeze.linalg.{DenseMatrix, DenseVector, norm}
 import breeze.util.Implicits._
 /**
- * Power Method to compute maximum eigen value
+ * Power Method to compute maximum eigen value and companion object to compute minimum eigen value through inverse power
+ * iterations
  * @author debasish83
  */
 class PowerMethod[T, M](maxIters: Int = 10,tolerance: Double = 1E-5)
-                       (implicit space: MutableInnerProductModule[T, Double],
-                        mult: OpMulMatrix.Impl2[M, T, T]) extends SerializableLogging {
+                       (implicit space: MutableInnerProductModule[T, Double], mult: OpMulMatrix.Impl2[M, T, T]) extends SerializableLogging {
 
   import space._
 
-  case class State(eigenValue: Double, eigenVector: T, iter: Int, converged: Boolean)
+  case class State private[PowerMethod] (eigenValue: Double, eigenVector: T, iter: Int, converged: Boolean)
 
   //memory allocation for the eigen vector result
   def normalize(y: T) : T = {
@@ -26,6 +26,13 @@ class PowerMethod[T, M](maxIters: Int = 10,tolerance: Double = 1E-5)
     init *= 1.0/normInit
   }
 
+  def initialState(y: T, A: M): State = {
+    val ynorm = normalize(y)
+    val ay = mult(A, ynorm)
+    val lambda = nextEigen(ynorm, ay)
+    State(lambda, ynorm, 0, false)
+  }
+
   //in-place modification of eigen vector
   def nextEigen(eigenVector: T, ay: T) = {
     val lambda = eigenVector dot ay
@@ -36,13 +43,6 @@ class PowerMethod[T, M](maxIters: Int = 10,tolerance: Double = 1E-5)
     lambda
   }
 
-  def initialState(y: T, A: M): State = {
-    val ynorm = normalize(y)
-    val ay = mult(A, ynorm)
-    val lambda = nextEigen(ynorm, ay)
-    State(lambda, ynorm, 0, false)
-  }
-
   def iterations(y: T,
                  A: M): Iterator[State] = Iterator.iterate(initialState(y, A)) { state =>
     import state._
@@ -63,7 +63,7 @@ class PowerMethod[T, M](maxIters: Int = 10,tolerance: Double = 1E-5)
 }
 
 object PowerMethod {
-  def inverse(maxIters: Int = 10, tolerance: Double = 1E-5) = new PowerMethod[DenseVector[Double], DenseMatrix[Double]](maxIters, tolerance) {
+  def inverse(maxIters: Int = 10, tolerance: Double = 1E-5) : PowerMethod[DenseVector[Double], DenseMatrix[Double]] = new PowerMethod[DenseVector[Double], DenseMatrix[Double]](maxIters, tolerance) {
     override def initialState(y: DenseVector[Double], A: DenseMatrix[Double]): State = {
       val ynorm = normalize(y)
       val ay = QuadraticMinimizer.solveTriangular(A, ynorm)

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

@@ -7,14 +7,15 @@
 package breeze.optimize.proximal
 
 import breeze.numerics.signum
-
 import scala.math.max
 import scala.math.min
 import scala.math.sqrt
 import scala.math.abs
 import scala.Double.NegativeInfinity
 import scala.Double.PositiveInfinity
 import breeze.linalg._
+import spire.syntax.cfor._
+import breeze.linalg.norm
 
 trait Proximal {
   def prox(x: DenseVector[Double], rho: Double = 1.0)
@@ -27,18 +28,16 @@ case class ProjectIdentity() extends Proximal {
 
 //TO DO:
 //1. Implement the binary search algorithm from http://see.stanford.edu/materials/lsocoee364b/hw4sol.pdf and compare performance
-//2. Implement randomized algorithm from Duchi et al's paper Efficient Projections onto the l1-Ball for Learning in High Dimensions
+//2. Implement randomized O(n) algorithm from Duchi et al's paper Efficient Projections onto the l1-Ball for Learning in High Dimensions
 case class ProjectProbabilitySimplex(s: Double) extends Proximal {
   require(s > 0, s"Proximal:ProjectProbabilitySimplex Radius s must be strictly positive")
   def prox(x: DenseVector[Double], rho: Double = 1.0)  = {
     val sorted = x.data.sorted(Ordering[Double].reverse)
     val cum = sorted.scanLeft(0.0)(_ + _).slice(1, x.length + 1)
-    val cs = DenseVector(cum.zipWithIndex.map{elem => (elem._1 - s)/(elem._2 + 1)})
-    val ndx = (DenseVector(sorted) - cs).data.filter{elem => elem >= 0.0}.length - 1
-    var i = 0
-    while(i < x.length) {
+    val cs = DenseVector(cum.zipWithIndex.map { elem => (elem._1 - s) / (elem._2 + 1)})
+    val ndx = (DenseVector(sorted) - cs).data.filter { elem => elem >= 0.0}.length - 1
+    cforRange(0 until x.length) { i =>
       x.update(i, max(x(i) - cs(ndx), 0.0))
-      i = i + 1
     }
   }
 }
@@ -49,63 +48,48 @@ case class ProjectProbabilitySimplex(s: Double) extends Proximal {
 case class ProjectL1(s: Double) extends Proximal {
   val projectSimplex = ProjectProbabilitySimplex(s)
 
-  def prox(x: DenseVector[Double], rho: Double = 1.0) : Unit = {
-    val u = x.mapValues {_.abs}
+  def prox(x: DenseVector[Double], rho: Double = 1.0): Unit = {
+    val u = x.mapValues {
+      _.abs
+    }
     projectSimplex.prox(u, rho)
-    var i = 0
-    while (i < x.length) {
+    cforRange(0 until x.length) { i =>
       x.update(i, signum(x(i)) * u(i))
-      i = i + 1
     }
   }
 }
 
 case class ProjectBox(l: DenseVector[Double], u: DenseVector[Double]) extends Proximal {
   def prox(x: DenseVector[Double], rho: Double = 0.0) = {
-    var i = 0
-    while (i < x.length) {
-      x.update(i, max(l(i), min(x(i), u(i))))
-      i = i + 1
-    }
+    cforRange(0 until x.length) { i => x.update(i, max(l(i), min(x(i), u(i))))}
   }
 }
 
 case class ProjectPos() extends Proximal {
   def prox(x: DenseVector[Double], rho: Double = 0.0) = {
-    var i = 0
-    while (i < x.length) {
-      x.update(i, max(0, x(i)))
-      i = i + 1
-    }
+    cforRange(0 until x.length) { i => x.update(i, max(0, x(i)))}
   }
 }
 
 case class ProjectSoc() extends Proximal {
   def prox(x: DenseVector[Double], rho: Double = 0.0) = {
     var nx: Double = 0.0
-    var i: Int = 1
     val n = x.length
-
-    while (i < n) {
+    cforRange(1 until n) { i =>
       nx += x(i) * x(i)
-      i = i + 1
     }
     nx = sqrt(nx)
 
     if (nx > x(0)) {
       if (nx <= -x(0)) {
-        i = 0
-        while (i < n) {
-          x(i) = 0
-          i = i + 1
+        cforRange(0 until n ) {
+          i => x(i) = 0
         }
       } else {
         val alpha = 0.5 * (1 + x(0) / nx)
         x.update(0, alpha * nx)
-        i = 1
-        while (i < n) {
-          x.update(i, alpha * x(i))
-          i = i + 1
+        cforRange(1 until n) {
+          i => x.update(i, alpha * x(i))
         }
       }
     }
@@ -147,10 +131,8 @@ 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) - lambda/rho) - max(0, -x(i) - lambda/rho))
-      i = i + 1
+    cforRange(0 until x.length) { i =>
+      x.update(i, max(0, x(i) - lambda / rho) - max(0, -x(i) - lambda / rho))
     }
   }
 
@@ -161,42 +143,28 @@ case class ProximalL1() extends Proximal {
 
 case class ProximalL2() extends Proximal {
   def prox(x: DenseVector[Double], rho: Double) = {
-    var normSquare: Double = 0.0
-    var i = 0
-
-    while (i < x.length) {
-      normSquare = normSquare + x(i) * x(i)
-      i = i + 1
-    }
-
-    val norm = sqrt(normSquare)
-    i = 0
-    while (i < x.length) {
-      if (norm >= 1 / rho) x.update(i, x(i) * (1 - 1 / (rho * norm)))
+    val xnorm = norm(x)
+    cforRange(0 until x.length) { i =>
+      if (xnorm >= 1 / rho) x.update(i, x(i) * (1 - 1 / (rho * xnorm)))
       else x.update(i, 0)
-      i = i + 1
     }
   }
 }
 
 // f = (1/2)||.||_2^2
 case class ProximalSumSquare() extends Proximal {
   def prox(x: DenseVector[Double], rho: Double) = {
-    var i = 0
-    while (i < x.length) {
+    cforRange(0 until x.length) { i =>
       x.update(i, x(i) * (rho / (1 + rho)))
-      i = i + 1
     }
   }
 }
 
 // f = -sum(log(x))
 case class ProximalLogBarrier() extends Proximal {
   def prox(x: DenseVector[Double], rho: Double) = {
-    var i = 0
-    while (i < x.length) {
+    cforRange(0 until x.length) { i =>
       x.update(i, 0.5 * (x(i) + sqrt(x(i) * x(i) + 4 / rho)))
-      i = i + 1
     }
   }
 }
@@ -232,11 +200,8 @@ case class ProximalHuber() extends Proximal {
 
   def proxSeparable(x: DenseVector[Double], rho: Double, oracle: Double => Double, l: Double, u: Double) = {
     x.map(proxScalar(_, rho, oracle, l, u, 0))
-
-    var i = 0
-    while (i < x.length) {
+    cforRange(0 until x.length) { i =>
       x.update(i, proxScalar(x(i), rho, oracle, l, u, 0))
-      i = i + 1
     }
   }
 
@@ -257,21 +222,17 @@ case class ProximalHuber() extends Proximal {
 // f = c'*x
 case class ProximalLinear(c: DenseVector[Double]) extends Proximal {
   def prox(x: DenseVector[Double], rho: Double) = {
-    var i = 0
-    while (i < x.length) {
+    cforRange(0 until x.length) { i =>
       x.update(i, x(i) - c(i) / rho)
-      i = i + 1
     }
   }
 }
 
 // f = c'*x + I(x >= 0)
 case class ProximalLp(c: DenseVector[Double]) extends Proximal {
   def prox(x: DenseVector[Double], rho: Double) = {
-    var i = 0
-    while (i < x.length) {
+    cforRange(0 until x.length) { i =>
       x.update(i, max(0, x(i) - c(i) / rho))
-      i = i + 1
     }
   }
 }