Added a proximal interface to ProjectQuasiNewton solver; Added projec…

…tion operators for ProbabilitySimplex and L1; Compared PQN with OWLQN and ADMM for L1 constraint

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

@@ -96,26 +96,6 @@ class LBFGS[T](maxIter: Int = -1, m: Int=10, tolerance: Double=1E-9)
     logger.info(s"LBFGS Min Eigen iterations ${eigenState.iter}")
     1.0/eigenState.eigenValue
   }
-
-  /**
-   * Get the maximum approximate eigen value of the quadratic approximation convex function specified
-   * through DiffFunction by doing a inverse power iteration using the CompactHessian representation
-   * generated from ApproximateInverseHessian, get the largest eigenvalue e
-   *
-   * @param init initial guess for the eigen vector
-   * @param state the current state of the optimization
-   * @return maximum eigen value
-   */
-  def maxEigen(state: State, init: T) : Double = {
-    val compactHessian = new CompactHessian(state.history.m)
-    val memStep = state.history.memStep.asInstanceOf[Seq[DenseVector[Double]]]
-    val memGradDelta = state.history.memGradDelta.asInstanceOf[Seq[DenseVector[Double]]]
-    (0 to state.history.historyLength - 1).map{i => compactHessian.updated(memStep(i), memGradDelta(i))}
-    val pm = new PowerMethod[T, CompactHessian]()(space, ProjectedQuasiNewton.multiplyCompactHessian)
-    val eigenState = pm.iterateAndReturnState(init, compactHessian)
-    logger.info(s"LBFGS Max Eigen iterations ${eigenState.iter}")
-    eigenState.eigenValue
-  }
 }
 
 object LBFGS {
@@ -180,6 +160,5 @@ object LBFGS {
       def apply(a: ApproximateInverseHessian[T], b: T): T = a * b
     }
   }
-
 }
 

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

@@ -2,8 +2,10 @@ package breeze.optimize
 
 import breeze.linalg._
 import breeze.collection.mutable.RingBuffer
-import breeze.linalg.operators.OpMulMatrix
 import breeze.math.{MutableInnerProductModule}
+import breeze.optimize.proximal.Constraint.Constraint
+import breeze.optimize.proximal.Constraint._
+import breeze.optimize.proximal._
 import breeze.util.SerializableLogging
 
 // Compact representation of an n x n Hessian, maintained via L-BFGS updates
@@ -84,7 +86,6 @@ class ProjectedQuasiNewton(tolerance: Double = 1e-6,
 
   type History = CompactHessian
 
-
   protected def initialHistory(f: DiffFunction[DenseVector[Double]], init: DenseVector[Double]): History = {
     new CompactHessian(m)
   }
@@ -109,7 +110,6 @@ class ProjectedQuasiNewton(tolerance: Double = 1e-6,
     }
   }
 
-
   protected def determineStepSize(state: State, fn: DiffFunction[DenseVector[Double]], dir: DenseVector[Double]): Double = {
     if (state.iter == 0)
       return scala.math.min(1.0, 1.0 / norm(state.grad,1.0))
@@ -153,14 +153,12 @@ class ProjectedQuasiNewton(tolerance: Double = 1e-6,
     projection(state.x + dir * stepSize)
   }
 
-
   protected def updateHistory(newX: DenseVector[Double], newGrad: DenseVector[Double], newVal: Double,  f: DiffFunction[DenseVector[Double]], oldState: State): History = {
     import oldState._
-    val s = newX - oldState.x
-    val y = newGrad - oldState.grad
+    val s = newX - x
+    val y = newGrad - grad
     oldState.history.updated(y, s)
   }
-
 }
 
 object ProjectedQuasiNewton {
@@ -186,12 +184,29 @@ object ProjectedQuasiNewton {
     }
   }
 
-  implicit def multiplyCompactHessian[T](implicit vspace: MutableInnerProductModule[T, Double]): OpMulMatrix.Impl2[CompactHessian, T, T] = {
-    new OpMulMatrix.Impl2[CompactHessian, T, T] {
-      def apply(a: CompactHessian, b: T): T = {
-        val result = a * b.asInstanceOf[DenseVector[Double]]
-        result.asInstanceOf[T]
+  case class Projection(proximal: Proximal) {
+    def project(x: DenseVector[Double]) : DenseVector[Double] = {
+      proximal.prox(x)
+      x
+    }
+  }
+
+  def proximal(proximal: Proximal) = new ProjectedQuasiNewton(projection = Projection(proximal).project)
+
+  def apply(ndim: Int, constraint: Constraint, lambda: Double=1.0): ProjectedQuasiNewton = {
+    constraint match {
+      case POSITIVE => proximal(ProjectPos())
+      case BOX => {
+        val lb = DenseVector.zeros[Double](ndim)
+        val ub = DenseVector.ones[Double](ndim)
+        proximal(ProjectBox(lb, ub))
+      }
+      case EQUALITY => {
+        val aeq = DenseVector.ones[Double](ndim)
+        proximal(ProjectHyperPlane(aeq, 1.0))
       }
+      case SPARSE => proximal(ProjectL1(lambda))
+      case PROBABILITYSIMPLEX => proximal(ProjectProbabilitySimplex(lambda))
     }
   }
 }

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

@@ -1,6 +1,6 @@
 package breeze.optimize
 
-import breeze.math.{Module, NormedVectorSpace}
+import breeze.math.{Module}
 
 trait Projecting[T] {
   def projection: T => T

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

@@ -6,5 +6,5 @@ package breeze.optimize.proximal
  */
 object Constraint extends Enumeration {
   type Constraint = Value
-  val SMOOTH, POSITIVE, BOX, SPARSE, EQUALITY = Value
+  val SMOOTH, POSITIVE, BOX, SPARSE, EQUALITY, PROBABILITYSIMPLEX = Value
 }