Quadratic Minimization solver using ADMM based Proximal Algorithm

LICENSE

@@ -174,3 +174,32 @@
       of your accepting any such warranty or additional liability.
 
    END OF TERMS AND CONDITIONS
+
+========================================================================
+BSD-style licenses
+========================================================================
+(BSD License) Proximal algorithms (https://github.com/cvxgrp/proximal)
+       Copyright (c) 2013, Neal Parikh and Stephen Boyd (Stanford University)
+       All rights reserved.
+
+       Redistribution and use in source and binary forms, with or without
+       modification, are permitted provided that the following conditions are met:
+       * Redistributions of source code must retain the above copyright
+         notice, this list of conditions and the following disclaimer.
+       * Redistributions in binary form must reproduce the above copyright
+         notice, this list of conditions and the following disclaimer in the
+         documentation and/or other materials provided with the distribution.
+       * Neither the name of Stanford University nor the
+         names of its contributors may be used to endorse or promote products
+         derived from this software without specific prior written permission.
+
+       THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+       ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+       WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+       DISCLAIMED. IN NO EVENT SHALL NEAL PARIKH OR STEPHEN BOYD BE LIABLE FOR ANY
+       DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+       (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+       LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+       ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+       (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+       SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

NOTICE

@@ -1 +1,13 @@
 Breeze is distributed under an Apache License V2.0 (See LICENSE)
+
+===============================================================================
+
+Proximal algorithms outlined in Proximal.scala (package breeze.optimize.quadratic)
+are based on https://github.com/cvxgrp/proximal (see LICENSE for details) and distributed with 
+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.
+

build.sbt

@@ -51,4 +51,3 @@ pomExtra := (
       <url>http://cs.berkeley.edu/~dlwh/</url>
     </developer>
   </developers>)
-

math/build.sbt

@@ -100,5 +100,3 @@ testOptions in Test += Tests.Argument("-oDF")
 
 //fork in Test := true
 javaOptions := Seq("-Xmx4g")
-
-jacoco.settings

math/src/main/scala/breeze/linalg/LinearAlgebraException.scala

@@ -31,3 +31,5 @@ class MatrixEmptyException extends IllegalArgumentException("Matrix is empty") w
  * @author dramage, dlwh
  */
 class MatrixSingularException(msg : String="") extends RuntimeException(msg)  with LinearAlgebraException
+
+class LapackException(msg: String="") extends RuntimeException(msg) with LinearAlgebraException

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

@@ -0,0 +1,189 @@
+package breeze.optimize.linear
+
+import breeze.linalg.{DenseMatrix, DenseVector, axpy}
+import breeze.optimize.proximal.QpGenerator
+import breeze.stats.distributions.Rand
+import breeze.util.Implicits._
+
+/**
+ * NNLS solves nonnegative least squares problems using a modified
+ * projected gradient method.
+ * @param maxIters user defined maximum iterations
+ * @author debasish83, coderxiang
+ */
+
+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)
+
+  // find the optimal unconstrained step
+  def steplen(ata: BDM, dir: BDV, res: BDV,
+              tmp: BDV): Double = {
+    val top = dir.dot(res)
+    tmp := ata * dir
+    // Push the denominator upward very slightly to avoid infinities and silliness
+    top / (tmp.dot(dir) + 1e-20)
+  }
+
+  // stopping condition
+  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
+      || (ndir < 1e-12 * nx) // gradient relatively too small
+      || (ndir < 1e-32) // gradient absolutely too small; numerical issues may lurk
+      )
+  }
+
+  /**
+   * Solve a least squares problem, possibly with nonnegativity constraints, by a modified
+   * 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
+   *
+   * 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 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 grad = 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)
+  }
+
+  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)
+
+      val iterMax = if (maxIters < 0) Math.max(400, 20 * n) else maxIters
+
+      // find the residual
+      res := ata * x
+      res -= atb
+      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
+        }
+        i = i + 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 (iter > lastWall + 1) {
+        val alpha = ngrad / lastNorm
+        axpy(alpha, lastDir, dir)
+        val dstep = steplen(ata, dir, res, tmp)
+        ndir = dir.dot(dir)
+        if (stop(dstep, ndir, nx)) {
+          // reject the CG step if it could lead to premature termination
+          dir := grad
+          ndir = dir.dot(dir)
+        } else {
+          step = dstep
+        }
+      } else {
+        ndir = dir.dot(dir)
+      }
+
+      // terminate?
+      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
+        }
+
+        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
+        }
+        lastDir := dir
+        val nextNorm = ngrad
+        State(x, grad, dir, lastDir, res, nextNorm, nextWall, iter + 1, false)
+      }
+    }.takeUpToWhere(_.converged)
+
+  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
+  }
+}
+
+object NNLS {
+  /** Compute the objective value */
+  def computeObjectiveValue(ata: DenseMatrix[Double], atb: DenseVector[Double], x: DenseVector[Double]): Double = {
+    val res = (x.t*ata*x)*0.5 - atb.dot(x)
+    res
+  }
+
+  def apply(iters: Int) = new NNLS(iters)
+
+  def main(args: Array[String]) {
+    if (args.length < 2) {
+      println("Usage: NNLS n s")
+      println("Test NNLS with quadratic function of dimension n for s consecutive solves")
+      sys.exit(1)
+    }
+
+    val problemSize = args(0).toInt
+    val numSolves = args(1).toInt
+    val nnls = new NNLS()
+
+    var i = 0
+    var nnlsTime = 0L
+    while(i < numSolves) {
+      val ata = QpGenerator.getGram(problemSize)
+      val atb = DenseVector.rand[Double](problemSize, Rand.gaussian(0,1))
+      val startTime = System.nanoTime()
+      nnls.minimize(ata, atb)
+      nnlsTime = nnlsTime + (System.nanoTime() - startTime)
+      i = i + 1
+    }
+    println(s"NNLS problemSize $problemSize solves $numSolves ${nnlsTime/1e6} ms")
+  }
+}

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

@@ -0,0 +1,64 @@
+package breeze.optimize.linear
+
+import breeze.linalg.operators.OpMulMatrix
+import breeze.math.{MutableInnerProductModule}
+import breeze.numerics.abs
+import breeze.util.SerializableLogging
+import breeze.linalg.{DenseMatrix, DenseVector, norm}
+import breeze.util.Implicits._
+import breeze.optimize.proximal.QuadraticMinimizer
+/**
+ * Power Method/Inverse Power Method to compute maximum/minimum eigen value
+ * @param inverse if true run inverse power method for DenseMatrix[Double], DenseVector[Double]
+ * @author debasish83
+ */
+class PowerMethod[T, M](maxIterations: Int = 10,tolerance: Double = 1E-5, inverse: Boolean = false)
+                       (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)
+
+  def initialState(y: T, A: M): State = {
+    //Force y to be a vector of unit norm
+    val normInit = norm(y)
+    y *= 1.0 / normInit
+    //TO DO : How to fix the asInstanceOf, at least throw a exception/require
+    val ay = if (inverse)
+      QuadraticMinimizer.solveTriangular(A.asInstanceOf[DenseMatrix[Double]], y.asInstanceOf[DenseVector[Double]]).asInstanceOf[T]
+    else mult(A, y)
+
+    val lambda = y dot ay
+
+    y := ay
+    val norm1 = norm(ay)
+    y *= 1.0/norm1
+    if (lambda < 0.0) y *= -1.0
+    State(lambda, y, 0, false)
+  }
+
+  def iterations(y: T,
+                 A: M): Iterator[State] = Iterator.iterate(initialState(y, A)) { state =>
+    import state._
+    val ay = if (inverse)
+      QuadraticMinimizer.solveTriangular(A.asInstanceOf[DenseMatrix[Double]], eigenVector.asInstanceOf[DenseVector[Double]]).asInstanceOf[T]
+      else mult(A, eigenVector)
+    val lambda = eigenVector dot ay
+    val norm1 = norm(ay)
+    ay *= 1.0 / norm1
+    if (lambda < 0.0) ay *= -1.0
+
+    val val_dif = abs(lambda - eigenValue)
+    if (val_dif <= tolerance || iter > maxIterations) State(lambda, ay, iter + 1, true)
+    else State(lambda, ay, iter + 1, false)
+  }.takeUpToWhere(_.converged)
+
+  def iterateAndReturnState(y: T, A: M): State = {
+    iterations(y, A).last
+  }
+
+  def eigen(y: T, A: M): Double = {
+    iterateAndReturnState(y, A).eigenValue
+  }
+}

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

@@ -0,0 +1,10 @@
+package breeze.optimize.proximal
+
+/**
+ * Supported constraints by QuadraticMinimizer object
+ * @author debasish83
+ */
+object Constraint extends Enumeration {
+  type Constraint = Value
+  val SMOOTH, POSITIVE, BOX, SPARSE, EQUALITY = Value
+}