Fixed LBFGS.maxEigen using power law on CompactHessian

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

@@ -90,10 +90,11 @@ class LBFGS[T](maxIter: Int = -1, m: Int=10, tolerance: Double=1E-9)
    * @param state the current state of the optimization
    * @return minimium eigen eigen value
    */
-  def minEigen(init:T, state: State) : Double = {
+  def minEigen(state: State, init: T) : Double = {
     val pm = new PowerMethod[T, LBFGS.ApproximateInverseHessian[T]]()(space,LBFGS.multiplyInverseHessian)
-    val eigenMax = pm.eigen(init, state.history)
-    1.0/eigenMax
+    val eigenState = pm.iterateAndReturnState(init, state.history)
+    logger.info(s"LBFGS Min Eigen iterations ${eigenState.iter}")
+    1.0/eigenState.eigenValue
   }
 
   /**
@@ -105,14 +106,15 @@ class LBFGS[T](maxIter: Int = -1, m: Int=10, tolerance: Double=1E-9)
    * @param state the current state of the optimization
    * @return maximum eigen value
    */
-  def maxEigen(init:T, state: State) : Double = {
-    println(s"m ${state.history.m} historyLength ${state.history.historyLength}")
+  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)
-    pm.eigen(init, compactHessian)
+    val eigenState = pm.iterateAndReturnState(init, compactHessian)
+    logger.info(s"LBFGS Max Eigen iterations ${eigenState.iter}")
+    eigenState.eigenValue
   }
 }
 

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

@@ -188,7 +188,10 @@ 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 = a * b
+      def apply(a: CompactHessian, b: T): T = {
+        val result = a * b.asInstanceOf[DenseVector[Double]]
+        result.asInstanceOf[T]
+      }
     }
   }
 }

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

@@ -186,9 +186,6 @@ class NonlinearMinimizer(ndim: Int,
     }
 
     var initialState = primalSolve.minimizeAndReturnState(primal, init)
-
-    println(s"Quadratic Model minEigen ${primalSolve.minEigen(init, initialState)}")
-
     val proxPrimal = ProximalPrimal(primal, u, z, rho)
 
     val xHat = DenseVector.zeros[Double](ndim)

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

@@ -17,20 +17,14 @@
 
 package breeze.optimize.proximal
 
-import breeze.linalg.cholesky
-import breeze.linalg.LU
+import breeze.linalg._
 import breeze.util.SerializableLogging
-import scala.math.max
 import scala.math.sqrt
 import breeze.optimize.{LBFGS, OWLQN, DiffFunction}
 import org.netlib.util.intW
 import breeze.optimize.proximal.Constraint._
 import scala.math.abs
-import breeze.linalg.DenseVector
-import breeze.linalg.DenseMatrix
 import breeze.numerics._
-import breeze.linalg.LapackException
-import breeze.linalg.norm
 import com.github.fommil.netlib.LAPACK.{getInstance=>lapack}
 import breeze.optimize.linear.{PowerMethod, NNLS, ConjugateGradient}
 import breeze.stats.distributions.Rand
@@ -418,7 +412,16 @@ object QuadraticMinimizer {
                          H: DenseMatrix[Double],
                          q: DenseVector[Double]) = {
     val lbfgs = new LBFGS[DenseVector[Double]](-1, 7)
-    lbfgs.minimize(Cost(H, q), init)
+    val state = lbfgs.minimizeAndReturnState(Cost(H, q), init)
+    val approxMinEigen = lbfgs.minEigen(state, init)
+    val approxMaxEigen = lbfgs.maxEigen(state, init)
+    val eigs = eigSym(H).eigenvalues
+
+    val minEigen = min(eigs)
+    val maxEigen = max(eigs)
+    println(s"minEigen $minEigen approx $approxMinEigen maxEigen $maxEigen approx $approxMaxEigen")
+
+    state.x
   }
 
   def optimizeWithOWLQN(init: DenseVector[Double],