diff --git a/NOTICE b/NOTICE index bd6ae4083..e40679a31 100644 --- a/NOTICE +++ b/NOTICE @@ -2,7 +2,7 @@ Breeze is distributed under an Apache License V2.0 (See LICENSE) =============================================================================== -Proximal algorithms outlined in Proximal.scala (package breeze.optimize.quadratic) +Proximal algorithms outlined in Proximal.scala (package breeze.optimize.proximal) 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. @@ -11,3 +11,7 @@ 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.proximal is distributed with Copyright (c) +2015, Debasish Das (Verizon), all rights reserved. diff --git a/README.md b/README.md index 5ae23e6b9..2cfb02ccd 100644 --- a/README.md +++ b/README.md @@ -100,7 +100,7 @@ Contributions from: * Chris Stucchio (@stucchio) * Xiangrui Meng (@mengxr) * Gabriel Schubiner (@gabeos) - +* Debasish Das (@debasish83) And others (contact David Hall if you've contributed code and aren't listed). diff --git a/math/src/main/scala/breeze/linalg/operators/DenseMatrixOps.scala b/math/src/main/scala/breeze/linalg/operators/DenseMatrixOps.scala index 9a0fb6213..2dc5137de 100644 --- a/math/src/main/scala/breeze/linalg/operators/DenseMatrixOps.scala +++ b/math/src/main/scala/breeze/linalg/operators/DenseMatrixOps.scala @@ -124,7 +124,7 @@ trait DenseMatrixMultiplyStuff extends DenseMatrixOps // square: LUSolve val X = DenseMatrix.zeros[Double](V.rows, V.cols) X := V - LUSolve(X,A) + LUSolve(X, A) X } else { // non-square: QRSolve @@ -137,15 +137,13 @@ trait DenseMatrixMultiplyStuff extends DenseMatrixOps /** X := A \ X, for square A */ def LUSolve(X: DenseMatrix[Double], A: DenseMatrix[Double]): DenseMatrix[Double] = { - require(X.offset == 0) - require(A.offset == 0) val piv = new Array[Int](A.rows) val newA = A.copy assert(!newA.isTranspose) val info: Int = { val info = new intW(0) - lapack.dgesv(A.rows, X.cols, newA.data, newA.majorStride, piv, X.data, X.majorStride, info) + lapack.dgesv(A.rows, X.cols, newA.data, newA.offset, newA.majorStride, piv, 0, X.data, X.offset, X.majorStride, info) info.`val` } diff --git a/math/src/main/scala/breeze/linalg/support/CanMapValues.scala b/math/src/main/scala/breeze/linalg/support/CanMapValues.scala index 9031d1173..979a9939b 100644 --- a/math/src/main/scala/breeze/linalg/support/CanMapValues.scala +++ b/math/src/main/scala/breeze/linalg/support/CanMapValues.scala @@ -51,14 +51,14 @@ trait CanMapValuesLowPrio { object CanMapValues extends CanMapValuesLowPrio { class HandHold[From, ValueType] - /* - implicit def canMapSelf[V, V2]: CanMapValues[V, V, V2, V2] = { - new CanMapValues[V, V, V2, V2] { - def map(from: V, fn: (V) => V2) = fn(from) - def mapActive(from: V, fn: (V) => V2) = fn(from) - } - } - */ + implicit def canMapSelfDouble[V2]: CanMapValues[Double, Double, V2, V2] = canMapSelf[Double, V2] + implicit def canMapSelfInt[V2]: CanMapValues[Int, Int, V2, V2] = canMapSelf[Int, V2] + implicit def canMapSelfFloat[V2]: CanMapValues[Float, Float, V2, V2] = canMapSelf[Float, V2] + implicit def canMapSelfLong[V2]: CanMapValues[Long, Long, V2, V2] = canMapSelf[Long, V2] + implicit def canMapSelfShort[V2]: CanMapValues[Short, Short, V2, V2] = canMapSelf[Short, V2] + implicit def canMapSelfByte[V2]: CanMapValues[Byte, Byte, V2, V2] = canMapSelf[Byte, V2] + implicit def canMapSelfChar[V2]: CanMapValues[Char, Char, V2, V2] = canMapSelf[Char, V2] + type Op[From, A, B, To] = CanMapValues[From, A, B, To] diff --git a/math/src/main/scala/breeze/optimize/BacktrackingLineSearch.scala b/math/src/main/scala/breeze/optimize/BacktrackingLineSearch.scala index 9863d5d33..a8e463086 100644 --- a/math/src/main/scala/breeze/optimize/BacktrackingLineSearch.scala +++ b/math/src/main/scala/breeze/optimize/BacktrackingLineSearch.scala @@ -8,7 +8,8 @@ package breeze.optimize * * @author dlwh */ -class BacktrackingLineSearch(maxIterations: Int = 20, +class BacktrackingLineSearch(initfval: Double, + maxIterations: Int = 20, shrinkStep: Double = 0.5, growStep: Double = 2.1, cArmijo: Double = 1E-4, @@ -24,7 +25,8 @@ class BacktrackingLineSearch(maxIterations: Int = 20, require(cWolfe < 1.0) def iterations(f: DiffFunction[Double], init: Double = 1.0): Iterator[State] = { val (f0, df0) = f.calculate(0.0) - val (initfval, initfderiv) = f.calculate(init) + val initfderiv = f.calculate(init)._2 + //val (initfval, initfderiv) = f.calculate(init) Iterator.iterate( (State(init, initfval, initfderiv), false, 0)) { case (state@State(alpha, fval, fderiv), _, iter) => val multiplier = if(fval > f0 + alpha * df0 * cArmijo) { shrinkStep diff --git a/math/src/main/scala/breeze/optimize/FirstOrderMinimizer.scala b/math/src/main/scala/breeze/optimize/FirstOrderMinimizer.scala index c34ce0fa3..07479c4c5 100644 --- a/math/src/main/scala/breeze/optimize/FirstOrderMinimizer.scala +++ b/math/src/main/scala/breeze/optimize/FirstOrderMinimizer.scala @@ -1,7 +1,7 @@ package breeze.optimize import breeze.linalg.norm -import breeze.math.{MutableEnumeratedCoordinateField, MutableCoordinateField, MutableFiniteCoordinateField, NormedModule} +import breeze.math.{MutableEnumeratedCoordinateField, MutableFiniteCoordinateField, NormedModule} import breeze.optimize.FirstOrderMinimizer.ConvergenceReason import breeze.stats.distributions.{RandBasis, ThreadLocalRandomGenerator} import breeze.util.Implicits._ @@ -108,10 +108,11 @@ abstract class FirstOrderMinimizer[T, DF<:StochasticDiffFunction[T]](maxIter: In f.calculate(x) } - def infiniteIterations(f: DF, init: T): Iterator[State] = { + def infiniteIterations(f: DF, state: State): Iterator[State] = { var failedOnce = false val adjustedFun = adjustFunction(f) - Iterator.iterate(initialState(adjustedFun,init)) { state => try { + + Iterator.iterate(state) { state => try { val dir = chooseDescentDirection(state, adjustedFun) val stepSize = determineStepSize(state, adjustedFun, dir) logger.info(f"Step Size: $stepSize%.4g") @@ -141,7 +142,8 @@ abstract class FirstOrderMinimizer[T, DF<:StochasticDiffFunction[T]](maxIter: In } def iterations(f: DF, init: T): Iterator[State] = { - infiniteIterations(f, init).takeUpToWhere(_.converged) + val adjustedFun = adjustFunction(f) + infiniteIterations(f, initialState(adjustedFun, init)).takeUpToWhere(_.converged) } def minimize(f: DF, init: T): T = { diff --git a/math/src/main/scala/breeze/optimize/LBFGS.scala b/math/src/main/scala/breeze/optimize/LBFGS.scala index 148e5311a..09e7bc432 100644 --- a/math/src/main/scala/breeze/optimize/LBFGS.scala +++ b/math/src/main/scala/breeze/optimize/LBFGS.scala @@ -19,12 +19,12 @@ package breeze.optimize import breeze.linalg._ import breeze.linalg.operators.OpMulMatrix import breeze.math.MutableInnerProductModule +import breeze.optimize.linear.PowerMethod import breeze.util.SerializableLogging - /** * Port of LBFGS to Scala. - * + * * Special note for LBFGS: * If you use it in published work, you must cite one of: * * J. Nocedal. Updating Quasi-Newton Matrices with Limited Storage @@ -80,7 +80,6 @@ class LBFGS[T](maxIter: Int = -1, m: Int=10, tolerance: Double=1E-9) throw new StepSizeUnderflow alpha } - } object LBFGS { @@ -140,13 +139,10 @@ object LBFGS { } } - implicit def multiplyInverseHessian[T](implicit vspace: MutableInnerProductModule[T, Double]):OpMulMatrix.Impl2[ApproximateInverseHessian[T], T, T] = { new OpMulMatrix.Impl2[ApproximateInverseHessian[T], T, T] { def apply(a: ApproximateInverseHessian[T], b: T): T = a * b } - } - } diff --git a/math/src/main/scala/breeze/optimize/OWLQN.scala b/math/src/main/scala/breeze/optimize/OWLQN.scala index 9b140a387..4aa090c58 100644 --- a/math/src/main/scala/breeze/optimize/OWLQN.scala +++ b/math/src/main/scala/breeze/optimize/OWLQN.scala @@ -70,7 +70,7 @@ class OWLQN[K, T](maxIter: Int, m: Int, l1reg: K => Double, tolerance: Double)(i adjv -> (adjgrad dot dir) } } - val search = new BacktrackingLineSearch(shrinkStep= if(iter < 1) 0.1 else 0.5) + val search = new BacktrackingLineSearch(state.value, shrinkStep= if(iter < 1) 0.1 else 0.5) val alpha = search.minimize(ff, if(iter < 1) .5/norm(state.grad) else 1.0) alpha diff --git a/math/src/main/scala/breeze/optimize/ProjectedQuasiNewton.scala b/math/src/main/scala/breeze/optimize/ProjectedQuasiNewton.scala index f27ba49c0..65d2e2a13 100644 --- a/math/src/main/scala/breeze/optimize/ProjectedQuasiNewton.scala +++ b/math/src/main/scala/breeze/optimize/ProjectedQuasiNewton.scala @@ -1,8 +1,24 @@ package breeze.optimize +/* + Copyright 2015 David Hall, Debasish Das + + Licensed under the Apache License, Version 2.0 (the "License") + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. +*/ + import breeze.linalg._ import breeze.collection.mutable.RingBuffer -import breeze.math.{MutableInnerProductModule, MutableVectorField} +import breeze.math.MutableInnerProductModule import breeze.util.SerializableLogging // Compact representation of an n x n Hessian, maintained via L-BFGS updates @@ -66,16 +82,18 @@ class ProjectedQuasiNewton(tolerance: Double = 1e-6, val m: Int = 10, val initFeas: Boolean = false, val testOpt: Boolean = true, - val maxNumIt: Int = 500, + maxIter: Int = -1, val maxSrchIt: Int = 50, val gamma: Double = 1e-4, val projection: DenseVector[Double] => DenseVector[Double] = identity) (implicit space: MutableInnerProductModule[DenseVector[Double],Double]) - extends FirstOrderMinimizer[DenseVector[Double], DiffFunction[DenseVector[Double]]](maxIter = maxNumIt, tolerance = tolerance) with Projecting[DenseVector[Double]] with SerializableLogging { - val innerOptimizer = new SpectralProjectedGradient[DenseVector[Double], DiffFunction[DenseVector[Double]]]( - testOpt = true, + extends FirstOrderMinimizer[DenseVector[Double], DiffFunction[DenseVector[Double]]](maxIter = maxIter, tolerance = tolerance) with Projecting[DenseVector[Double]] with SerializableLogging { + type BDV = DenseVector[Double] + + val innerOptimizer = new SpectralProjectedGradient[BDV]( tolerance = tolerance, - maxIter = 500, + maxIter = 50, + bbMemory = 5, initFeas = true, minImprovementWindow = 10, projection = projection @@ -83,15 +101,13 @@ class ProjectedQuasiNewton(tolerance: Double = 1e-6, type History = CompactHessian - protected def initialHistory(f: DiffFunction[DenseVector[Double]], init: DenseVector[Double]): History = { new CompactHessian(m) } - override protected def adjust(newX: DenseVector[Double], newGrad: DenseVector[Double], newVal: Double):(Double,DenseVector[Double]) = (newVal,-projectedVector(newX, -newGrad)) + override protected def adjust(newX: DenseVector[Double], newGrad: DenseVector[Double], newVal: Double):(Double,DenseVector[Double]) = (newVal,projectedVector(newX, -newGrad)) private def computeGradient(x: DenseVector[Double], g: DenseVector[Double]): DenseVector[Double] = projectedVector(x, -g) - private def computeGradientNorm(x: DenseVector[Double], g: DenseVector[Double]): Double = norm(computeGradient(x, g),Double.PositiveInfinity) protected def chooseDescentDirection(state: State, fn: DiffFunction[DenseVector[Double]]): DenseVector[Double] = { import state._ @@ -101,72 +117,54 @@ class ProjectedQuasiNewton(tolerance: Double = 1e-6, // Update the limited-memory BFGS approximation to the Hessian //B.update(y, s) // Solve subproblem; we use the current iterate x as a guess - val subprob = new ProjectedQuasiNewton.QuadraticSubproblem(fn, state.adjustedValue, x, grad, history) - val p = innerOptimizer.minimize(new CachedDiffFunction(subprob), x) - p - x + val subprob = new ProjectedQuasiNewton.QuadraticSubproblem(state.adjustedValue, x, grad, history) + val spgResult = innerOptimizer.minimizeAndReturnState(new CachedDiffFunction(subprob), x) + logger.info(f"ProjectedQuasiNewton: outerIter ${state.iter} innerIters ${spgResult.iter}") + spgResult.x - x // time += subprob.time } } - - 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)) - val dirnorm = norm(dir, Double.PositiveInfinity) - if(dirnorm < 1E-10) return 0.0 - import state._ - // Backtracking line-search - var accepted = false - var lambda = 1.0 - val gTd = grad dot dir - var srchit = 0 - - do { - val candx = x + dir * lambda - val candf = fn.valueAt(candx) - val suffdec = gamma * lambda * gTd - - if (testOpt && srchit > 0) { - logger.debug(f"PQN: SrchIt $srchit%4d: f $candf%-10.4f t $lambda%-10.4f\n") - } - - if (candf < state.adjustedValue + suffdec) { - accepted = true - } else if (srchit >= maxSrchIt) { - accepted = true - } else { - lambda *= 0.5 - srchit = srchit + 1 - } - } while (!accepted) - - if (srchit >= maxSrchIt) { - logger.info("PQN: Line search cannot make further progress") - throw new LineSearchFailed(norm(state.grad,Double.PositiveInfinity), norm(dir, Double.PositiveInfinity)) - } - lambda + /** + * Given a direction, perform a Strong Wolfe Line Search + * + * TO DO: Compare performance with Cubic Interpolation based line search from Mark's PQN paper + * + * @param state the current state + * @param f The objective + * @param dir The step direction + * @return stepSize + */ + protected def determineStepSize(state: State, f: DiffFunction[DenseVector[Double]], dir: DenseVector[Double]) = { + val x = state.x + val grad = state.grad + + val ff = LineSearch.functionFromSearchDirection(f, x, dir) + val search = new BacktrackingLineSearch(state.value, maxIterations = maxSrchIt, shrinkStep= if(state.iter < 1) 0.1 else 0.5) + var alpha = if(state.iter == 0.0) min(1.0, 1.0/norm(dir)) else 1.0 + alpha = search.minimize(ff, alpha) + + if(alpha * norm(grad) < 1E-10) throw new StepSizeUnderflow + + alpha } - protected def takeStep(state: State, dir: DenseVector[Double], stepSize: Double): DenseVector[Double] = { 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 { +object ProjectedQuasiNewton extends SerializableLogging { // Forms a quadratic model around fun, the argmin of which is then a feasible // quasi-Newton descent direction - class QuadraticSubproblem(fun: DiffFunction[DenseVector[Double]], - fk: Double, + class QuadraticSubproblem(fk: Double, xk: DenseVector[Double], gk: DenseVector[Double], B: CompactHessian) extends DiffFunction[DenseVector[Double]] { diff --git a/math/src/main/scala/breeze/optimize/Projecting.scala b/math/src/main/scala/breeze/optimize/Projecting.scala index 5effdb154..b54c9cc5d 100644 --- a/math/src/main/scala/breeze/optimize/Projecting.scala +++ b/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 diff --git a/math/src/main/scala/breeze/optimize/SpectralProjectedGradient.scala b/math/src/main/scala/breeze/optimize/SpectralProjectedGradient.scala index dbc079bee..4c5afae16 100644 --- a/math/src/main/scala/breeze/optimize/SpectralProjectedGradient.scala +++ b/math/src/main/scala/breeze/optimize/SpectralProjectedGradient.scala @@ -1,84 +1,131 @@ package breeze.optimize -import breeze.math.{MutableInnerProductModule, MutableVectorField} +/* + Copyright 2015 David Hall, Daniel Ramage, Debasish Das + + Licensed under the Apache License, Version 2.0 (the "License") + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. +*/ + +import breeze.math.{InnerProductModule, MutableVectorField} import breeze.linalg.norm import breeze.util.SerializableLogging - /** * SPG is a Spectral Projected Gradient minimizer; it minimizes a differentiable - * function subject to the optimum being in some set, given by the projection operator projection + * function subject to the optimum being in some set, given by the projection operator projection * @tparam T vector type - * @param optTol termination criterion: tolerance for norm of projected gradient - * @param gamma sufficient decrease parameter - * @param M number of history entries for linesearch + * @param tolerance termination criterion: tolerance for norm of projected gradient + * @param suffDec sufficient decrease parameter + * @param bbMemory number of history entries for linesearch * @param alphaMax longest step * @param alphaMin shortest step - * @param maxNumIt maximum number of iterations - * @param testOpt perform optimality check based on projected gradient at each iteration + * @param maxIter maximum number of iterations + * @param maxSrcht maximum number of iterations inside line search * @param initFeas is the initial guess feasible, or should it be projected? - * @param maxSrchIt maximum number of line search attempts * @param projection projection operations + * @param curvilinear if curvilinear true, do the projection inside line search in place of doing it in chooseDescentDirection */ -class SpectralProjectedGradient[T, DF <: DiffFunction[T]]( +class SpectralProjectedGradient[T]( val projection: T => T = { (t: T) => t }, tolerance: Double = 1e-6, - val suffDec: Double = 1e-4, - minImprovementWindow: Int = 10, - val alphaMax: Double = 1e10, - val alphaMin: Double = 1e-10, - maxIter: Int = 500, - val testOpt: Boolean = true, + suffDec: Double = 1e-4, + minImprovementWindow: Int = 30, + alphaMax: Double = 1e10, + alphaMin: Double = 1e-10, + bbMemory: Int = 10, + maxIter: Int = -1, val initFeas: Boolean = false, - val maxSrchIt: Int = 30)(implicit space: MutableVectorField[T, Double]) extends FirstOrderMinimizer[T, DF](minImprovementWindow = minImprovementWindow, maxIter = maxIter, tolerance = tolerance) with Projecting[T] with SerializableLogging { + val curvilinear: Boolean = false, + val bbType: Int = 1, + val maxSrcht: Int = 30)(implicit space: MutableVectorField[T, Double]) extends FirstOrderMinimizer[T, DiffFunction[T]](minImprovementWindow = minImprovementWindow, maxIter = maxIter, tolerance = tolerance, improvementTol = tolerance) with Projecting[T] with SerializableLogging { import space._ - type History = Double - protected def initialHistory(f: DF, init: T): History = 1.0 - protected def chooseDescentDirection(state: State, f: DF): T = projectedVector(state.x, state.grad :* -state.history) - override protected def adjust(newX: T, newGrad: T, newVal: Double):(Double,T) = (newVal,-projectedVector(newX, - newGrad)) - protected def takeStep(state: State, dir: T, stepSize: Double): T = projection(state.x + dir :* stepSize) - protected def updateHistory(newX: T, newGrad: T, newVal: Double, f: DF, oldState: State): History = { - val y = newGrad - oldState.grad + case class History(alphaBB: Double, fvals: IndexedSeq[Double]) + + override protected def initialHistory(f: DiffFunction[T], init: T): History = { + History(0.1, IndexedSeq.empty) + } + + /** + * From Mark Schmidt's Matlab code + * if bbType == 1 + * alpha = (s'*s)/(s'*y); + * else + * alpha = (s'*y)/(y'*y); + */ + protected def bbAlpha(s: T, y: T) : Double = { + var alpha = + if (bbType == 1) (s dot s) / (s dot y) + else (s dot y) / (y dot y) + if (alpha <= alphaMin || alpha > alphaMax) alpha = 1.0 + if (alpha.isNaN) alpha = 1.0 + alpha + } + + override protected def updateHistory(newX: T, newGrad: T, newVal: Double, f: DiffFunction[T], oldState: State): History = { val s = newX - oldState.x - val alpha = s.dot(s) / s.dot(y) - if (alpha.isNaN()) - 0.0 - else if (alpha < alphaMin || alpha > alphaMax) - 1 - else - alpha + val y = newGrad - oldState.grad + History(bbAlpha(s, y), (newVal +: oldState.history.fvals).take(bbMemory)) } - protected def determineStepSize(state: State, f: DF, direction: T): Double = { - import state._ - val funRef = if (fVals.isEmpty) Double.PositiveInfinity else fVals.max - val t = if (iter == 0) { - scala.math.min(1.0, 1.0 / norm(grad)) - } else { - 1.0 + override protected def takeStep(state: State, dir: T, stepSize: Double): T = { + val qq = projection(state.x + dir * stepSize) + assert(projection(qq) == qq) + qq + } + + override protected def chooseDescentDirection(state: State, f: DiffFunction[T]): T = { + if (curvilinear) state.x - state.grad * state.history.alphaBB + else projection(state.x - state.grad * state.history.alphaBB) - state.x + } + + override protected def determineStepSize(state: State, f: DiffFunction[T], direction: T): Double = { + val fb = if (state.history.fvals.isEmpty) state.value else state.value max state.history.fvals.max + val normGradInDir = state.grad dot direction + + var gamma = + if (state.iter == 0) scala.math.min(1.0, 1.0 / norm(state.grad)) + else 1.0 + + val searchFun = + if (curvilinear) functionFromSearchDirection(f, state.x, direction, projection) + else LineSearch.functionFromSearchDirection(f, state.x, direction) + + //TO DO : + // 1. Add cubic interpolation and see it's performance. Bisection did not work for L1 projection + val search = new BacktrackingLineSearch(fb, maxIterations=maxSrcht) + gamma = search.minimize(searchFun, gamma) + + if (gamma < 1e-10) { + throw new LineSearchFailed(normGradInDir, norm(direction)) } - val searchStep = direction * t - val sufficientDecrease = grad.dot(searchStep) * suffDec - val requiredValue = funRef + sufficientDecrease - val lineSearchFunction = LineSearch.functionFromSearchDirection(f, x, direction) - val ls = new SimpleLineSearch(requiredValue, maxSrchIt) - ls.minimize(lineSearchFunction, t) + gamma } - class SimpleLineSearch(requiredValue: Double, maxIterations: Int) extends ApproximateLineSearch { - def iterations(f: DiffFunction[Double], init: Double = 1.0): Iterator[State] = { - val (initfval, initfderiv) = f.calculate(init) - Iterator.iterate((State(init, initfval, initfderiv), 0)) { - case (State(alpha, fval, fderiv), iter) => - val newAlpha = alpha / 2.0 - val (fvalnew, fderivnew) = f.calculate(newAlpha) - (State(newAlpha, fvalnew, fderivnew), iter + 1) - }.takeWhile { - case (state, iterations) => - (iterations == 0) || - (iterations < maxIterations && - state.value > requiredValue) - }.map(_._1) + + // because of the projection, we have to do our own verstion + private def functionFromSearchDirection[T, I](f: DiffFunction[T], x: T, direction: T, project: T=>T)(implicit prod: InnerProductModule[T, Double]):DiffFunction[Double] = new DiffFunction[Double] { + import prod._ + + /** calculates the value at a point */ + override def valueAt(alpha: Double): Double = f.valueAt(project(x + direction * alpha)) + + /** calculates the gradient at a point */ + override def gradientAt(alpha: Double): Double = f.gradientAt(project(x + direction * alpha)) dot direction + + /** Calculates both the value and the gradient at a point */ + def calculate(alpha: Double): (Double, Double) = { + val (ff, grad) = f.calculate(x + direction * alpha) + ff -> (grad dot direction) } } } diff --git a/math/src/main/scala/breeze/optimize/proximal/Constraint.scala b/math/src/main/scala/breeze/optimize/proximal/Constraint.scala index c3f647551..d7729565f 100644 --- a/math/src/main/scala/breeze/optimize/proximal/Constraint.scala +++ b/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 IDENTITY, SMOOTH, POSITIVE, BOX, SPARSE, EQUALITY, PROBABILITYSIMPLEX = Value } \ No newline at end of file diff --git a/math/src/main/scala/breeze/optimize/proximal/LinearGenerator.scala b/math/src/main/scala/breeze/optimize/proximal/LinearGenerator.scala new file mode 100644 index 000000000..a89abb3c2 --- /dev/null +++ b/math/src/main/scala/breeze/optimize/proximal/LinearGenerator.scala @@ -0,0 +1,37 @@ +package breeze.optimize.proximal + +import breeze.linalg.{DenseVector, DenseMatrix} +import breeze.optimize.DiffFunction +import breeze.stats.distributions.Rand + +/** + * @author debasish83 + */ +object LinearGenerator { + case class Cost(data: DenseMatrix[Double], labels: DenseVector[Double]) extends DiffFunction[DenseVector[Double]] { + def calculate(x: DenseVector[Double]) = { + val cumGradient = DenseVector.zeros[Double](x.length) + var cumLoss = 0.0 + var i = 0 + while (i < data.rows) { + val brzData = data(i, ::).t + val diff = x.dot(brzData) - labels(i) + cumGradient += brzData * (2.0 * diff) + cumLoss += diff * diff + i = i + 1 + } + (cumLoss, cumGradient) + } + } + + def apply(ndim: Int) : (DiffFunction[DenseVector[Double]], DenseMatrix[Double], DenseVector[Double]) = { + val rand = Rand.gaussian(0, 1) + val data = DenseMatrix.rand[Double](ndim, ndim, rand) + val labels = DenseVector.rand[Double](ndim, rand).map { x => if (x > 0.5) 1.0 else 0.0} + //||ax - b||_2^{2} = x'a'ax - 2*x'a'*b + c + val h = (data.t*data)*2.0 + val q = (data.t*labels) + q *= -2.0 + (Cost(data, labels), h, q) + } +} \ No newline at end of file diff --git a/math/src/main/scala/breeze/optimize/proximal/LogisticGenerator.scala b/math/src/main/scala/breeze/optimize/proximal/LogisticGenerator.scala new file mode 100644 index 000000000..f07f46fbb --- /dev/null +++ b/math/src/main/scala/breeze/optimize/proximal/LogisticGenerator.scala @@ -0,0 +1,44 @@ +package breeze.optimize.proximal + +import breeze.linalg.{DenseVector, DenseMatrix} +import breeze.optimize.DiffFunction +import breeze.stats.distributions.Rand + +/** + * @author debasish83 + */ +object LogisticGenerator { + + case class Cost(data: DenseMatrix[Double], + labels: DenseVector[Double]) extends DiffFunction[DenseVector[Double]] { + def calculate(x: DenseVector[Double]) = { + val cumGradient = DenseVector.zeros[Double](x.length) + var cumLoss = 0.0 + + var i = 0 + while (i < data.rows) { + val brzData = data(i, ::).t + val margin: Double = -1.0 * x.dot(brzData) + val gradientMultiplier = (1.0 / (1.0 + math.exp(margin))) - labels(i) + val gradient = brzData * gradientMultiplier + val loss = + if (labels(i) > 0) { + math.log1p(math.exp(margin)) // log1p is log(1+p) but more accurate for small p + } else { + math.log1p(math.exp(margin)) - margin + } + cumGradient += gradient + cumLoss += loss + i = i + 1 + } + (cumLoss, cumGradient) + } + } + + def apply(ndim: Int): DiffFunction[DenseVector[Double]] = { + val rand = Rand.gaussian(0, 1) + val data = DenseMatrix.rand[Double](ndim, ndim, rand) + val labels = DenseVector.rand[Double](ndim, rand).map { x => if (x > 0.5) 1.0 else 0.0} + Cost(data, labels) + } +} \ No newline at end of file diff --git a/math/src/main/scala/breeze/optimize/proximal/NonlinearMinimizer.scala b/math/src/main/scala/breeze/optimize/proximal/NonlinearMinimizer.scala new file mode 100644 index 000000000..9a3cdf74e --- /dev/null +++ b/math/src/main/scala/breeze/optimize/proximal/NonlinearMinimizer.scala @@ -0,0 +1,380 @@ +package breeze.optimize.proximal + +import breeze.linalg.{norm, DenseVector} +import breeze.math.MutableInnerProductModule +import breeze.optimize._ +import breeze.util.SerializableLogging +import breeze.optimize.proximal.Constraint._ +import scala.math._ +import breeze.linalg.DenseMatrix +import breeze.linalg.sum +import breeze.linalg.max +import breeze.util.Implicits._ + +/** + * + * NonlinearMinimizer solves the problem that has the following structure + * minimize f(x) + g(x) + * + * g(x) represents the following constraints + * + * 1. x >= 0 + * 2. lb <= x <= ub + * 3. L1(x) + * 4. Aeq*x = beq + * 5. aeq'x = beq + * 6. 1'x = s, x >= 0 ProbabilitySimplex from the reference Proximal Algorithms by Boyd et al, Duchi et al + * + * f(x) can be a smooth convex function defined by DiffFunction or a proximal operator. For now the exposed API + * takes DiffFunction + * + * g(x) is defined by a proximal operator + * + * For proximal algorithms like L1 through soft-thresholding and Huber Loss (look into library of proximal + * algorithms for further details) we provide ADMM based Proximal algorithm based on the following + * reference: https://web.stanford.edu/~boyd/papers/admm/logreg-l1/distr_l1_logreg.html + * + * A subset of proximal operators are projection operators. For projection, NonlinearMinimizer companion + * object provides project API which generates a Spectral Projected Gradient (SPG) or Projected Quasi Newton (PQN) + * solver. For projection operators like positivity, bounds, probability simplex etc, these algorithms converges + * faster as compared to ADMM based proximal algorithm. + * + * TO DO + * + * 1. Implement FISTA / Nesterov's accelerated method and compare with ADMM + * 2. For non-convex function experiment with TRON-like Primal solver + * + * @author debasish83 + */ + +class NonlinearMinimizer(proximal: Proximal, maxIters: Int = -1, innerIters: Int = 3, bfgsMemory: Int = 7, + rho: Double = 1.0, alpha: Double = 1.0, + abstol: Double = 1e-6, reltol: Double = 1e-4) extends SerializableLogging { + import NonlinearMinimizer.BDV + + val lbfgs = new LBFGS[BDV](m=bfgsMemory,tolerance=abstol,maxIter=innerIters) + + case class State private[NonlinearMinimizer](bfgsState: lbfgs.State, u: BDV, z: BDV, xHat: BDV, zOld: BDV, residual: BDV, s: BDV, admmIters: Int, iter: Int, converged: Boolean) + + private def initialState(primal: DiffFunction[BDV], init: BDV) = { + val z = init.copy + val u = init.copy + + val xHat = init.copy + val zOld = init.copy + + val residual = init.copy + val s = init.copy + + val resultState = lbfgs.minimizeAndReturnState(primal, xHat) + val admmIters = if (maxIters < 0) max(400, 20*z.length) else maxIters + State(resultState, u, z, xHat, zOld, residual, s, admmIters, 0, false) + } + + def iterations(primal: DiffFunction[BDV], init: BDV): Iterator[State] = Iterator.iterate(initialState(primal, init)) { state => + import state._ + + val scale = sqrt(init.size) * abstol + val proxPrimal = NonlinearMinimizer.ProximalPrimal(primal, u, z, rho) + + val resultState = lbfgs.minimizeAndReturnState(proxPrimal, bfgsState.x) + //z-update with relaxation + + //zold = (1-alpha)*z + //x_hat = alpha*x + zold + + zOld := z + zOld *= 1 - alpha + + xHat := resultState.x + xHat *= alpha + xHat += zOld + + //zold = z + zOld := z + + //z = xHat + u + z := xHat + z += u + + //Apply proximal operator + proximal.prox(z, rho) + + //z has proximal(x_hat) + + //Dual (u) update + xHat -= z + u += xHat + + //Convergence checks + //history.r_norm(k) = norm(x - z) + residual := resultState.x + residual -= z + val residualNorm = norm(residual) + + //history.s_norm(k) = norm(-rho*(z - zold)) + s := z + s -= zOld + s *= -rho + val sNorm = norm(s) + + //TO DO : Make sure z.muli(-1) is actually needed in norm calculation + residual := z + residual *= -1.0 + + //s = rho*u + s := u + s *= rho + + val epsPrimal = scale + reltol * max(norm(resultState.x), norm(residual)) + val epsDual = scale + reltol * norm(s) + + if (residualNorm < epsPrimal && sNorm < epsDual || iter > admmIters) { + State(resultState, u, z, xHat, zOld, residual, s, admmIters, iter + 1, true) + } else { + State(resultState, u, z, xHat, zOld, residual, s, admmIters, iter + 1, false) + } + }.takeUpToWhere{_.converged} + + def minimize(primal: DiffFunction[BDV], init: BDV): BDV = { + minimizeAndReturnState(primal, init).z + } + + def minimizeAndReturnState(primal: DiffFunction[BDV], init: BDV): State = { + iterations(primal, init).last + } +} + +object NonlinearMinimizer { + type BDV = DenseVector[Double] + + /** + * Proximal modifications to Primal algorithm for scaled ADMM formulation + * AdmmObj(x, u, z) = f(x) + rho/2*||x - z + u||2 + * dAdmmObj/dx = df/dx + rho*(x - z + u) + */ + case class ProximalPrimal[T](primal: DiffFunction[T], + u: T, + z: T, + rho: Double) + (implicit space: MutableInnerProductModule[T, Double]) extends DiffFunction[T] { + + import space._ + + override def calculate(x: T) = { + val (f, g) = primal.calculate(x) + val scale = x - z + u + val proxObj = f + 0.5 * rho * pow(norm(scale), 2) + val proxGrad = g + scale :* rho + (proxObj, proxGrad) + } + } + + case class Projection(proximal: Proximal) { + def project(x: BDV): BDV = { + proximal.prox(x) + x + } + } + + /** + * A subset of proximal operators can be represented as Projection operators and for those + * operators, we give an option to the user to choose a projection based algorithm. The options + * available for users are SPG (Spectral Projected Gradient) and PQN (Projected Quasi Newton) + * + * @param proximal operator that defines proximal algorithm + * @return FirstOrderMinimizer to optimize on f(x) and proximal operator + */ + def project(proximal: Proximal, maxIter: Int = -1, m: Int = 10, tolerance: Double = 1e-6, usePQN: Boolean = false): FirstOrderMinimizer[BDV, DiffFunction[BDV]] = { + val projectionOp = Projection(proximal) + if (usePQN) new ProjectedQuasiNewton(projection = projectionOp.project, tolerance = 1e-6, maxIter = maxIter, m = m) + else new SpectralProjectedGradient(projection = projectionOp.project, tolerance = tolerance, maxIter = maxIter, bbMemory = m) + } + + /** + * A compansion object to generate projection based minimizer that can use SPG/PQN as the solver + * @param ndim the problem dimension + * @param constraint one of the available constraint, possibilities are x>=0; lb<=x<=ub;aeq*x = beq; + * 1'x = s, x >= 0; ||x||1 <= s + * @param lambda the regularization parameter for most of the constraints + * @return FirstOrderMinimizer to optimize on f(x) and proximal operator + */ + def apply(ndim: Int, constraint: Constraint, lambda: Double, usePQN: Boolean = false): FirstOrderMinimizer[BDV, DiffFunction[BDV]] = { + constraint match { + case IDENTITY => project(ProjectIdentity()) + case POSITIVE => project(ProjectPos()) + case BOX => { + val lb = DenseVector.zeros[Double](ndim) + val ub = DenseVector.ones[Double](ndim) + project(ProjectBox(lb, ub)) + } + case EQUALITY => { + val aeq = DenseVector.ones[Double](ndim) + project(ProjectHyperPlane(aeq, 1.0)) + } + case PROBABILITYSIMPLEX => project(ProjectProbabilitySimplex(lambda)) + case SPARSE => project(ProjectL1(lambda)) + case _ => throw new IllegalArgumentException("NonlinearMinimizer does not support the Projection Operator") + } + } + + def main(args: Array[String]) { + if (args.length < 3) { + println("Usage: ProjectedQuasiNewton n lambda beta") + println("Test NonlinearMinimizer with a quadratic function of dimenion n and m equalities with lambda beta for elasticNet") + sys.exit(1) + } + + val problemSize = args(0).toInt + val lambda = args(1).toDouble + val beta = args(2).toDouble + + println(s"Generating Linear and Logistic Loss with rank ${problemSize}") + + val (quadraticCost, h, q) = LinearGenerator(problemSize) + + val lambdaL1 = lambda * beta + val lambdaL2 = lambda * (1 - beta) + + val owlqn = new OWLQN[Int, DenseVector[Double]](-1, 10, lambdaL1, 1e-6) + + val regularizedGram = h + (DenseMatrix.eye[Double](h.rows) :* lambdaL2) + + val sparseQp = QuadraticMinimizer(h.rows, SPARSE, lambdaL1) + val sparseQpStart = System.nanoTime() + val sparseQpResult = sparseQp.minimizeAndReturnState(regularizedGram, q) + val sparseQpTime = System.nanoTime() - sparseQpStart + + val init = DenseVector.zeros[Double](problemSize) + + val owlqnStart = System.nanoTime() + val owlqnResult = QuadraticMinimizer.optimizeWithOWLQN(init, regularizedGram, q, lambdaL1) + val owlqnTime = System.nanoTime() - owlqnStart + + println("ElasticNet Formulation") + + println("Linear Regression") + + val owlqnObj = QuadraticMinimizer.computeObjective(regularizedGram, q, owlqnResult.x) + lambdaL1 * owlqnResult.x.foldLeft(0.0) { (agg, entry) => agg + abs(entry)} + val sparseQpL1Obj = sparseQpResult.x.foldLeft(0.0) { (agg, entry) => agg + abs(entry)} + val sparseQpObj = QuadraticMinimizer.computeObjective(regularizedGram, q, sparseQpResult.x) + lambdaL1 * sparseQpL1Obj + val quadraticCostWithL2 = QuadraticMinimizer.Cost(regularizedGram, q) + + init := 0.0 + val projectL1Linear = ProjectL1(sparseQpL1Obj) + val nlSparseStart = System.nanoTime() + val nlSparseResult = NonlinearMinimizer.project(projectL1Linear).minimizeAndReturnState(quadraticCostWithL2, init) + val nlSparseTime = System.nanoTime() - nlSparseStart + val nlSparseL1Obj = nlSparseResult.x.foldLeft(0.0) { (agg, entry) => agg + abs(entry)} + val nlSparseObj = QuadraticMinimizer.computeObjective(regularizedGram, q, nlSparseResult.x) + lambdaL1 * nlSparseL1Obj + + init := 0.0 + val nlProx = new NonlinearMinimizer(proximal = ProximalL1(lambdaL1)) + val nlProxStart = System.nanoTime() + val nlProxResult = nlProx.minimizeAndReturnState(quadraticCostWithL2, init) + val nlProxTime = System.nanoTime() - nlProxStart + val nlProxObj = QuadraticMinimizer.computeObjective(regularizedGram, q, nlProxResult.z) + lambdaL1 * nlProxResult.z.foldLeft(0.0) { (agg, entry) => agg + abs(entry)} + + println(s"owlqn ${owlqnTime / 1e6} ms iters ${owlqnResult.iter} sparseQp ${sparseQpTime / 1e6} ms iters ${sparseQpResult.iter}") + println(s"nlSparseTime ${nlSparseTime / 1e6} ms iters ${nlSparseResult.iter}") + println(s"nlProxTime ${nlProxTime / 1e6} ms iters ${nlProxResult.iter}") + println(s"owlqnObj $owlqnObj sparseQpObj $sparseQpObj nlSparseObj $nlSparseObj nlProxObj $nlProxObj") + + val logisticLoss = LogisticGenerator(problemSize) + val elasticNetLoss = DiffFunction.withL2Regularization(logisticLoss, lambdaL2) + + println("Linear Regression with Bounds") + + init := 0.0 + val nlBox = NonlinearMinimizer(problemSize, BOX, 0.0) + + val nlBoxStart = System.nanoTime() + val nlBoxResult = nlBox.minimizeAndReturnState(quadraticCostWithL2, init) + val nlBoxTime = System.nanoTime() - nlBoxStart + val nlBoxObj = quadraticCostWithL2.calculate(nlBoxResult.x)._1 + + val qpBox = QuadraticMinimizer(problemSize, BOX, 0.0) + val qpBoxStart = System.nanoTime() + val qpBoxResult = qpBox.minimizeAndReturnState(regularizedGram, q) + val qpBoxTime = System.nanoTime() - qpBoxStart + val qpBoxObj = QuadraticMinimizer.computeObjective(regularizedGram, q, qpBoxResult.x) + + println(s"qpBox ${qpBoxTime / 1e6} ms iters ${qpBoxResult.iter}") + println(s"nlBox ${nlBoxTime / 1e6} ms iters ${nlBoxResult.iter}") + println(s"qpBoxObj $qpBoxObj nlBoxObj $nlBoxObj") + + println("Logistic Regression with Bounds") + + init := 0.0 + val nlBoxLogisticStart = System.nanoTime() + val nlBoxLogisticResult = nlBox.minimizeAndReturnState(elasticNetLoss, init) + val nlBoxLogisticTime = System.nanoTime() - nlBoxLogisticStart + val nlBoxLogisticObj = elasticNetLoss.calculate(nlBoxLogisticResult.x)._1 + println(s"Objective nl ${nlBoxLogisticObj} time ${nlBoxLogisticTime / 1e6} ms") + + println("Linear Regression with ProbabilitySimplex") + + val nlSimplex = NonlinearMinimizer(problemSize, PROBABILITYSIMPLEX, 1.0) + + init := 0.0 + val nlSimplexStart = System.nanoTime() + val nlSimplexResult = nlSimplex.minimizeAndReturnState(quadraticCost, init) + val nlSimplexTime = System.nanoTime() - nlSimplexStart + val nlSimplexObj = quadraticCost.calculate(nlSimplexResult.x)._1 + + val qpSimplexStart = System.nanoTime() + val qpSimplexResult = QuadraticMinimizer(problemSize, EQUALITY).minimizeAndReturnState(h, q) + val qpSimplexTime = System.nanoTime() - qpSimplexStart + val qpSimplexObj = quadraticCost.calculate(qpSimplexResult.x)._1 + + println(s"Objective nl $nlSimplexObj qp $qpSimplexObj") + println(s"Constraint nl ${sum(nlSimplexResult.x)} qp ${sum(qpSimplexResult.x)}") + println(s"time nl ${nlSimplexTime / 1e6} ms qp ${qpSimplexTime / 1e6} ms") + + println("Logistic Regression with ProbabilitySimplex") + + init := 0.0 + val nlLogisticSimplexStart = System.nanoTime() + val nlLogisticSimplexResult = nlSimplex.minimizeAndReturnState(elasticNetLoss, init) + val nlLogisticSimplexTime = System.nanoTime() - nlLogisticSimplexStart + val nlLogisticSimplexObj = elasticNetLoss.calculate(nlLogisticSimplexResult.x)._1 + + init := 0.0 + + val nlProxSimplex = new NonlinearMinimizer(proximal = ProjectProbabilitySimplex(1.0)) + val nlProxLogisticSimplexStart = System.nanoTime() + val nlProxLogisticSimplexResult = nlProxSimplex.minimizeAndReturnState(elasticNetLoss, init) + val nlProxLogisticSimplexTime = System.nanoTime() - nlProxLogisticSimplexStart + val nlProxLogisticSimplexObj = elasticNetLoss.calculate(nlProxLogisticSimplexResult.z)._1 + + println(s"Objective nl ${nlLogisticSimplexObj} admm ${nlProxLogisticSimplexObj}") + println(s"Constraint nl ${sum(nlLogisticSimplexResult.x)} admm ${sum(nlProxLogisticSimplexResult.z)}") + println(s"time nlProjection ${nlLogisticSimplexTime / 1e6} ms nlProx ${nlProxLogisticSimplexTime / 1e6} ms") + + println("Logistic Regression with ProximalL1 and ProjectL1") + + init := 0.0 + val owlqnLogisticStart = System.nanoTime() + val owlqnLogisticResult = owlqn.minimizeAndReturnState(elasticNetLoss,init) + val owlqnLogisticTime = System.nanoTime() - owlqnLogisticStart + val owlqnLogisticObj = elasticNetLoss.calculate(owlqnLogisticResult.x)._1 + val s = owlqnLogisticResult.x.foldLeft(0.0) { case (agg, entry) => agg + abs(entry)} + + init := 0.0 + val proximalL1 = ProximalL1().setLambda(lambdaL1) + val nlLogisticProximalL1Start = System.nanoTime() + val nlLogisticProximalL1Result = new NonlinearMinimizer(proximalL1).minimizeAndReturnState(elasticNetLoss, init) + val nlLogisticProximalL1Time = System.nanoTime() - nlLogisticProximalL1Start + + init := 0.0 + val nlLogisticProjectL1Start = System.nanoTime() + val nlLogisticProjectL1Result = NonlinearMinimizer(problemSize, SPARSE, s).minimizeAndReturnState(elasticNetLoss, init) + val nlLogisticProjectL1Time = System.nanoTime() - nlLogisticProjectL1Start + + val proximalL1Obj = elasticNetLoss.calculate(nlLogisticProximalL1Result.z)._1 + val projectL1Obj = elasticNetLoss.calculate(nlLogisticProjectL1Result.x)._1 + + println(s"Objective proximalL1 $proximalL1Obj projectL1 $projectL1Obj owlqn $owlqnLogisticObj") + println(s"time proximalL1 ${nlLogisticProximalL1Time / 1e6} ms projectL1 ${nlLogisticProjectL1Time / 1e6} ms owlqn ${owlqnLogisticTime / 1e6} ms") + } +} diff --git a/math/src/main/scala/breeze/optimize/proximal/Proximal.scala b/math/src/main/scala/breeze/optimize/proximal/Proximal.scala index ecb7c4dd5..863b34487 100644 --- a/math/src/main/scala/breeze/optimize/proximal/Proximal.scala +++ b/math/src/main/scala/breeze/optimize/proximal/Proximal.scala @@ -6,8 +6,7 @@ package breeze.optimize.proximal -import breeze.numerics.pow - +import breeze.numerics.signum import scala.math.max import scala.math.min import scala.math.sqrt @@ -19,7 +18,45 @@ import spire.syntax.cfor._ import breeze.linalg.norm trait Proximal { - def prox(x: DenseVector[Double], rho: Double) + def prox(x: DenseVector[Double], rho: Double = 1.0) + def valueAt(x: DenseVector[Double]) = 0.0 +} + +case class ProjectIdentity() extends Proximal { + def prox(x: DenseVector[Double], rho: Double = 1.0) {} +} + +//TO DO: +//1. Implement the binary search algorithm from http://see.stanford.edu/materials/lsocoee364b/hw4sol.pdf and compare performance +//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 + cforRange(0 until x.length) { i => + x.update(i, max(x(i) - cs(ndx), 0.0)) + } + } +} + +/** + * Projection formula from Duchi et al's paper Efficient Projections onto the l1-Ball for Learning in High Dimensions + * */ +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 + } + projectSimplex.prox(u, rho) + cforRange(0 until x.length) { i => + x.update(i, signum(x(i)) * u(i)) + } + } } case class ProjectBox(l: DenseVector[Double], u: DenseVector[Double]) extends Proximal { @@ -85,8 +122,7 @@ case class ProjectHyperPlane(a: DenseVector[Double], b: Double) extends Proximal } } -case class ProximalL1() extends Proximal { - var lambda = 1.0 +case class ProximalL1(var lambda: Double = 1.0) extends Proximal { def setLambda(lambda: Double) = { this.lambda = lambda @@ -98,6 +134,10 @@ case class ProximalL1() extends Proximal { x.update(i, max(0, x(i) - lambda / rho) - max(0, -x(i) - lambda / rho)) } } + + override def valueAt(x: DenseVector[Double]) = { + lambda * x.foldLeft(0.0) { (agg, entry) => agg + abs(entry)} + } } case class ProximalL2() extends Proximal { diff --git a/math/src/main/scala/breeze/optimize/proximal/QpGenerator.scala b/math/src/main/scala/breeze/optimize/proximal/QpGenerator.scala new file mode 100644 index 000000000..8aef3c552 --- /dev/null +++ b/math/src/main/scala/breeze/optimize/proximal/QpGenerator.scala @@ -0,0 +1,50 @@ +package breeze.optimize.proximal + +import breeze.linalg.{DenseVector, DenseMatrix} +import breeze.stats.distributions.Rand +/** + * PDCO dense quadratic program generator + * + * Reference + * + * Generates random instances of Quadratic Programming Problems + * 0.5x'Px + q'x + * s.t Ax = b + * lb <= x <= ub + * + * nGram rank of quadratic problems to be generated + * @author debasish83 + * @return A is the equality constraint + * @return b is the equality parameters + * @return lb is vector of lower bounds (default at 0.0) + * @return ub is vector of upper bounds (default at 10.0) + * @return q is linear representation of the function + * @return H is the quadratic representation of the function + */ +object QpGenerator { + def getGram(nGram: Int) = { + val hrand = DenseMatrix.rand[Double](nGram, nGram, Rand.gaussian(0, 1)) + val hrandt = hrand.t + val hposdef = hrandt * hrand + val H = hposdef.t + hposdef + H + } + + def apply(nGram: Int, nEqualities: Int) = { + val en = DenseVector.ones[Double](nGram) + val zn = DenseVector.zeros[Double](nGram) + + val A = DenseMatrix.rand[Double](nEqualities, nGram) + val x = en + + val b = A * x + val q = DenseVector.rand[Double](nGram) + + val lb = zn.copy + val ub = en :* 10.0 + + val H = getGram(nGram) + + (A, b, lb, ub, q, H) + } +} \ No newline at end of file diff --git a/math/src/main/scala/breeze/optimize/proximal/QuadraticMinimizer.scala b/math/src/main/scala/breeze/optimize/proximal/QuadraticMinimizer.scala index 52ec56146..576846ec9 100644 --- a/math/src/main/scala/breeze/optimize/proximal/QuadraticMinimizer.scala +++ b/math/src/main/scala/breeze/optimize/proximal/QuadraticMinimizer.scala @@ -1,19 +1,13 @@ 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 @@ -244,7 +238,7 @@ class QuadraticMinimizer(nGram: Int, private def computeRho(H: DenseMatrix[Double]): Double = { proximal match { case null => 0.0 - case ProximalL1() => { + case ProximalL1(lambda:Double) => { val eigenMax = QuadraticMinimizer.normColumn(H) val eigenMin = QuadraticMinimizer.approximateMinEigen(H) sqrt(eigenMin * eigenMax) @@ -281,52 +275,6 @@ class QuadraticMinimizer(nGram: Int, } } -/** - * PDCO dense quadratic program generator - * - * Reference - * - * Generates random instances of Quadratic Programming Problems - * 0.5x'Px + q'x - * s.t Ax = b - * lb <= x <= ub - * - * nGram rank of quadratic problems to be generated - * @return A is the equality constraint - * @return b is the equality parameters - * @return lb is vector of lower bounds (default at 0.0) - * @return ub is vector of upper bounds (default at 10.0) - * @return q is linear representation of the function - * @return H is the quadratic representation of the function - */ -object QpGenerator { - def getGram(nGram: Int) = { - val hrand = DenseMatrix.rand[Double](nGram, nGram, Rand.gaussian(0, 1)) - val hrandt = hrand.t - val hposdef = hrandt * hrand - val H = hposdef.t + hposdef - H - } - - def apply(nGram: Int, nEqualities: Int) = { - val en = DenseVector.ones[Double](nGram) - val zn = DenseVector.zeros[Double](nGram) - - val A = DenseMatrix.rand[Double](nEqualities, nGram) - val x = en - - val b = A * x - val q = DenseVector.rand[Double](nGram) - - val lb = zn.copy - val ub = en :* 10.0 - - val H = getGram(nGram) - - (A, b, lb, ub, q, H) - } -} - object QuadraticMinimizer { //upper bound on max eigen value def normColumn(H: DenseMatrix[Double]): Double = { @@ -396,8 +344,9 @@ object QuadraticMinimizer { def apply(rank: Int, constraint: Constraint, - lambda: Double): QuadraticMinimizer = { + lambda: Double=1.0): QuadraticMinimizer = { constraint match { + case SMOOTH => new QuadraticMinimizer(rank) case POSITIVE => new QuadraticMinimizer(rank, ProjectPos()) case BOX => { //Direct QP with bounds @@ -420,16 +369,27 @@ object QuadraticMinimizer { res } + case class Cost(H: DenseMatrix[Double], + q: DenseVector[Double]) extends DiffFunction[DenseVector[Double]] { + def calculate(x: DenseVector[Double]) = { + (computeObjective(H, q, x), H * x + q) + } + } + def optimizeWithLBFGS(init: DenseVector[Double], H: DenseMatrix[Double], q: DenseVector[Double]) = { val lbfgs = new LBFGS[DenseVector[Double]](-1, 7) - val f = new DiffFunction[DenseVector[Double]] { - def calculate(x: DenseVector[Double]) = { - (computeObjective(H, q, x), H * x + q) - } - } - lbfgs.minimize(f, init) + val state = lbfgs.minimizeAndReturnState(Cost(H, q), init) + state.x + } + + def optimizeWithOWLQN(init: DenseVector[Double], + regularizedGram: DenseMatrix[Double], + q: DenseVector[Double], + lambdaL1: Double) = { + val owlqn = new OWLQN[Int, DenseVector[Double]](-1, 7, lambdaL1, 1e-6) + owlqn.minimizeAndReturnState(Cost(regularizedGram, q), init) } def main(args: Array[String]) { @@ -486,24 +446,13 @@ object QuadraticMinimizer { val regularizedGram = h + (DenseMatrix.eye[Double](h.rows) :* lambdaL2) - val owlqn = new OWLQN[Int, DenseVector[Double]](-1, 7, lambdaL1) - - def optimizeWithOWLQN(init: DenseVector[Double]) = { - val f = new DiffFunction[DenseVector[Double]] { - def calculate(x: DenseVector[Double]) = { - (computeObjective(regularizedGram, q, x), regularizedGram * x + q) - } - } - owlqn.minimizeAndReturnState(f, init) - } - val sparseQp = QuadraticMinimizer(h.rows, SPARSE, lambdaL1) val sparseQpStart = System.nanoTime() val sparseQpResult = sparseQp.minimizeAndReturnState(regularizedGram, q) val sparseQpTime = System.nanoTime() - sparseQpStart val startOWLQN = System.nanoTime() - val owlqnResult = optimizeWithOWLQN(DenseVector.rand[Double](problemSize)) + val owlqnResult = optimizeWithOWLQN(DenseVector.rand[Double](problemSize), regularizedGram, q, lambdaL1) val owlqnTime = System.nanoTime() - startOWLQN println(s"||owlqn - sparseqp|| norm ${norm(owlqnResult.x - sparseQpResult.x, 2)} inf-norm ${norm(owlqnResult.x - sparseQpResult.x, inf)}") diff --git a/math/src/test/scala/breeze/optimize/ProjectedQuasiNewtonTest.scala b/math/src/test/scala/breeze/optimize/ProjectedQuasiNewtonTest.scala index 2dfdaec31..c95a7d83f 100644 --- a/math/src/test/scala/breeze/optimize/ProjectedQuasiNewtonTest.scala +++ b/math/src/test/scala/breeze/optimize/ProjectedQuasiNewtonTest.scala @@ -1,7 +1,7 @@ package breeze.optimize /* - Copyright 2009 David Hall, Daniel Ramage - + Copyright 2015 David Hall, Daniel Ramage, Debasish Das + Licensed under the Apache License, Version 2.0 (the "License") you may not use this file except in compliance with the License. You may obtain a copy of the License at @@ -15,15 +15,14 @@ package breeze.optimize limitations under the License. */ +import breeze.optimize.proximal.{ProjectL1, QuadraticMinimizer} import org.scalatest._ import org.scalatest.junit._ import org.scalatest.prop._ -import org.scalacheck._ import org.junit.runner.RunWith -import org.scalatest.matchers.ShouldMatchers - import breeze.linalg._ import breeze.numerics._ +import breeze.optimize.proximal.NonlinearMinimizer.Projection @RunWith(classOf[JUnitRunner]) class ProjectedQuasiNewtonTest extends PropSpec with PropertyChecks with OptimizeTestBaseTrait with VectorMatchers with Matchers { @@ -46,6 +45,7 @@ class ProjectedQuasiNewtonTest extends PropSpec with PropertyChecks with Optimiz val optimizer = new ProjectedQuasiNewton(tolerance = 1.0E-5, projection = _.map(scala.math.min(_, 2.0))) forAll { init: DenseVector[Double] => + init := clip(init, Double.NegativeInfinity, 2.0) val f = new DiffFunction[DenseVector[Double]] { def calculate(x: DenseVector[Double]) = { (sum((x - 3.0) :^ 4.0), (x - 3.0) :^ 3.0 :* 4.0) @@ -91,4 +91,43 @@ class ProjectedQuasiNewtonTest extends PropSpec with PropertyChecks with Optimiz } } } + + property("simple linear solve without projection") { + val n = 5 + val H = new DenseMatrix(n, n, Array(1.8984250861699135,0.5955576666769438,-1.484430453342902,-1.0434994471390804,-3.675310432634351, + 0.5955576666769438,0.9090751938470876,-2.146380947361661,-0.13037609428980368,-0.40639564652095117, + -1.484430453342902,-2.146380947361661,10.262733520770384,-6.097698907163584,2.29625304115155, + -1.0434994471390804,-0.13037609428980368,-6.097698907163584,27.775920405610677,-5.574220233644466, + -3.675310432634351,-0.40639564652095117,2.29625304115155,-5.574220233644466,12.21329172136971)) + val f = DenseVector(-1.2320199653150048, -0.14220655875869606, 0.38477404739124765, -0.3480575854151014, -0.4729810900829228) + + val cost = QuadraticMinimizer.Cost(H, f:*(-1.0)) + val init = DenseVector.zeros[Double](n) + + init := 0.0 + val x = H \ f + + init := 0.0 + val nlResult = new ProjectedQuasiNewton().minimize(cost, init) + assert(norm(x - nlResult, inf) < 1e-4) + } + + property("simple linear solve with L1 projection compared to Octave") { + val Hl1 = new DenseMatrix(25, 25, Array(253.535098, 236.477785, 234.421906, 223.374867, 241.007512, 204.695511, 226.465507, 223.351032, 249.179386, 221.411909, 238.679352, 203.579010, 217.564498, 243.852681, 266.607649, 213.994496, 241.620759, 223.602907, 220.038678, 264.704959, 240.341716, 223.672378, 244.323303, 223.216217, 226.074990, 236.477785, 278.862035, 245.756639, 237.489890, 252.783139, 214.077652, 241.816953, 238.790633, 260.536460, 228.121417, 255.103936, 216.608405, 237.150426, 258.933231, 281.958112, 228.971242, 252.508513, 234.242638, 240.308477, 285.416390, 254.792243, 240.176223, 259.048267, 235.566855, 236.277617, 234.421906, 245.756639, 269.162882, 231.416867, 251.321527, 208.134322, 236.567647, 236.558029, 255.805108, 226.535825, 251.514713, 212.770208, 228.565362, 261.748652, 273.946966, 227.411615, 252.767900, 232.823977, 233.084574, 278.315614, 250.872786, 235.227909, 255.104263, 238.931093, 235.402356, 223.374867, 237.489890, 231.416867, 254.771963, 241.703229, 209.028084, 231.517998, 228.768510, 250.805315, 216.548935, 245.473869, 207.687875, 222.036114, 250.906955, 263.018181, 216.128966, 244.445283, 227.436840, 231.369510, 270.721492, 242.475130, 226.471530, 248.130112, 225.826557, 228.266719, 241.007512, 252.783139, 251.321527, 241.703229, 285.702320, 219.051868, 249.442308, 240.400187, 264.970407, 232.503138, 258.819837, 220.160683, 235.621356, 267.743972, 285.795029, 229.667231, 260.870105, 240.751687, 247.183922, 289.044453, 260.715749, 244.210258, 267.159502, 242.992822, 244.070245, 204.695511, 214.077652, 208.134322, 209.028084, 219.051868, 210.164224, 208.151908, 201.539036, 226.373834, 192.056565, 219.950686, 191.459568, 195.982460, 226.739575, 240.677519, 196.116652, 217.352348, 203.533069, 204.581690, 243.603643, 217.785986, 204.205559, 223.747953, 203.586842, 200.165867, 226.465507, 241.816953, 236.567647, 231.517998, 249.442308, 208.151908, 264.007925, 227.080718, 253.174653, 220.322823, 248.619983, 210.100242, 223.279198, 254.807401, 269.896959, 222.927882, 247.017507, 230.484479, 233.358639, 274.935489, 249.237737, 235.229584, 253.029955, 228.601700, 230.512885, 223.351032, 238.790633, 236.558029, 228.768510, 240.400187, 201.539036, 227.080718, 258.773479, 249.471480, 215.664539, 243.078577, 202.337063, 221.020998, 249.979759, 263.356244, 213.470569, 246.182278, 225.727773, 229.873732, 266.295057, 242.954024, 225.510760, 249.370268, 227.813265, 232.141964, 249.179386, 260.536460, 255.805108, 250.805315, 264.970407, 226.373834, 253.174653, 249.471480, 302.360150, 237.902729, 265.769812, 224.947876, 243.088105, 273.690377, 291.076027, 241.089661, 267.772651, 248.459822, 249.662698, 295.935799, 267.460908, 255.668926, 275.902272, 248.495606, 246.827505, 221.411909, 228.121417, 226.535825, 216.548935, 232.503138, 192.056565, 220.322823, 215.664539, 237.902729, 245.154567, 234.956316, 199.557862, 214.774631, 240.339217, 255.161923, 209.328714, 232.277540, 216.298768, 220.296241, 253.817633, 237.638235, 220.785141, 239.098500, 220.583355, 218.962732, 238.679352, 255.103936, 251.514713, 245.473869, 258.819837, 219.950686, 248.619983, 243.078577, 265.769812, 234.956316, 288.133073, 225.087852, 239.810430, 268.406605, 283.289840, 233.858455, 258.306589, 240.263617, 246.844456, 290.492875, 267.212598, 243.218596, 265.681905, 244.615890, 242.543363, 203.579010, 216.608405, 212.770208, 207.687875, 220.160683, 191.459568, 210.100242, 202.337063, 224.947876, 199.557862, 225.087852, 217.501685, 197.897572, 229.825316, 242.175607, 201.123644, 219.820165, 202.894307, 211.468055, 246.048907, 225.135194, 210.076305, 226.806762, 212.014431, 205.123267, 217.564498, 237.150426, 228.565362, 222.036114, 235.621356, 195.982460, 223.279198, 221.020998, 243.088105, 214.774631, 239.810430, 197.897572, 244.439113, 241.621129, 260.400953, 216.482178, 236.805076, 216.680343, 223.816297, 263.188711, 236.311810, 222.950152, 244.636356, 219.121372, 219.911078, 243.852681, 258.933231, 261.748652, 250.906955, 267.743972, 226.739575, 254.807401, 249.979759, 273.690377, 240.339217, 268.406605, 229.825316, 241.621129, 302.928261, 288.344398, 238.549018, 267.239982, 248.073140, 254.230916, 296.789984, 267.158551, 252.226496, 271.170860, 248.325354, 253.694013, 266.607649, 281.958112, 273.946966, 263.018181, 285.795029, 240.677519, 269.896959, 263.356244, 291.076027, 255.161923, 283.289840, 242.175607, 260.400953, 288.344398, 343.457361, 257.368309, 284.795470, 263.122266, 271.239770, 320.209823, 283.933299, 264.416752, 292.035194, 268.764031, 265.345807, 213.994496, 228.971242, 227.411615, 216.128966, 229.667231, 196.116652, 222.927882, 213.470569, 241.089661, 209.328714, 233.858455, 201.123644, 216.482178, 238.549018, 257.368309, 239.295031, 234.913508, 218.066855, 219.648997, 257.969951, 231.243624, 224.657569, 238.158714, 217.174368, 215.933866, 241.620759, 252.508513, 252.767900, 244.445283, 260.870105, 217.352348, 247.017507, 246.182278, 267.772651, 232.277540, 258.306589, 219.820165, 236.805076, 267.239982, 284.795470, 234.913508, 289.709239, 241.312315, 247.249491, 286.702147, 264.252852, 245.151647, 264.582984, 240.842689, 245.837476, 223.602907, 234.242638, 232.823977, 227.436840, 240.751687, 203.533069, 230.484479, 225.727773, 248.459822, 216.298768, 240.263617, 202.894307, 216.680343, 248.073140, 263.122266, 218.066855, 241.312315, 255.363057, 230.209787, 271.091482, 239.220241, 225.387834, 247.486715, 226.052431, 224.119935, 220.038678, 240.308477, 233.084574, 231.369510, 247.183922, 204.581690, 233.358639, 229.873732, 249.662698, 220.296241, 246.844456, 211.468055, 223.816297, 254.230916, 271.239770, 219.648997, 247.249491, 230.209787, 264.014907, 271.938970, 246.664305, 227.889045, 249.908085, 232.035369, 229.010298, 264.704959, 285.416390, 278.315614, 270.721492, 289.044453, 243.603643, 274.935489, 266.295057, 295.935799, 253.817633, 290.492875, 246.048907, 263.188711, 296.789984, 320.209823, 257.969951, 286.702147, 271.091482, 271.938970, 352.825726, 286.200221, 267.716897, 297.182554, 269.776351, 266.721561, 240.341716, 254.792243, 250.872786, 242.475130, 260.715749, 217.785986, 249.237737, 242.954024, 267.460908, 237.638235, 267.212598, 225.135194, 236.311810, 267.158551, 283.933299, 231.243624, 264.252852, 239.220241, 246.664305, 286.200221, 294.042749, 246.504021, 269.570596, 243.980697, 242.690997, 223.672378, 240.176223, 235.227909, 226.471530, 244.210258, 204.205559, 235.229584, 225.510760, 255.668926, 220.785141, 243.218596, 210.076305, 222.950152, 252.226496, 264.416752, 224.657569, 245.151647, 225.387834, 227.889045, 267.716897, 246.504021, 259.897656, 251.730847, 229.335712, 229.759185, 244.323303, 259.048267, 255.104263, 248.130112, 267.159502, 223.747953, 253.029955, 249.370268, 275.902272, 239.098500, 265.681905, 226.806762, 244.636356, 271.170860, 292.035194, 238.158714, 264.582984, 247.486715, 249.908085, 297.182554, 269.570596, 251.730847, 303.872223, 251.585636, 247.878402, 223.216217, 235.566855, 238.931093, 225.826557, 242.992822, 203.586842, 228.601700, 227.813265, 248.495606, 220.583355, 244.615890, 212.014431, 219.121372, 248.325354, 268.764031, 217.174368, 240.842689, 226.052431, 232.035369, 269.776351, 243.980697, 229.335712, 251.585636, 257.544914, 228.810942, 226.074990, 236.277617, 235.402356, 228.266719, 244.070245, 200.165867, 230.512885, 232.141964, 246.827505, 218.962732, 242.543363, 205.123267, 219.911078, 253.694013, 265.345807, 215.933866, 245.837476, 224.119935, 229.010298, 266.721561, 242.690997, 229.759185, 247.878402, 228.810942, 253.353769)) + val fl1 = DenseVector(-892.842851, -934.071560, -932.936015, -888.124343, -961.050207, -791.191087, -923.711397, -904.289301, -988.384984, -883.909133, -959.465030, -798.551172, -871.622303, -997.463289, -1043.912620, -863.013719, -976.975712, -897.033693, -898.694786, -1069.245497, -963.491924, -901.263474, -983.768031, -899.865392, -902.283567) + val cost = QuadraticMinimizer.Cost(Hl1, fl1) + + val octaveL1 = DenseVector(0.18611, 0.00000, 0.06317, -0.10417, 0.11262, + -0.20495, 0.52668, 0.32790, 0.19421, 0.72180, + 0.06309, -0.41326, -0.00000, 0.52078, -0.00000, + 0.18040, 0.62915, 0.16329, -0.06424, 0.37539, + 0.01659, 0.00000, 0.11215, 0.24778, 0.04082) + + val s = octaveL1.foldLeft(0.0) { case (agg, entry) => agg + abs(entry)} + val projection = Projection(ProjectL1(s)) + val pqnResult = new ProjectedQuasiNewton(projection=projection.project).minimizeAndReturnState(cost, DenseVector.zeros[Double](25)) + + println(s"L1 projection iter ${pqnResult.iter}") + assert(norm(pqnResult.x - octaveL1, 2) < 1e-4) + } } diff --git a/math/src/test/scala/breeze/optimize/SpectralProjectedGradientTest.scala b/math/src/test/scala/breeze/optimize/SpectralProjectedGradientTest.scala new file mode 100644 index 000000000..72eb326ae --- /dev/null +++ b/math/src/test/scala/breeze/optimize/SpectralProjectedGradientTest.scala @@ -0,0 +1,113 @@ +/* + * + * Copyright 2015 David Hall + * + * Licensed under the Apache License, Version 2.0 (the "License") + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * / + */ + +package breeze.optimize + +/* + Copyright 2009 David Hall, Daniel Ramage + + Licensed under the Apache License, Version 2.0 (the "License") + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. +*/ + +import breeze.linalg._ +import breeze.numerics._ +import breeze.optimize.proximal.NonlinearMinimizer.Projection +import breeze.optimize.proximal.{ProjectL1, QuadraticMinimizer} +import org.junit.runner.RunWith +import org.scalatest._ +import org.scalatest.junit._ +import org.scalatest.prop._ + +@RunWith(classOf[JUnitRunner]) +class SpectralProjectedGradientTest extends PropSpec with PropertyChecks with OptimizeTestBaseTrait with VectorMatchers with Matchers { + + + property("optimize a simple multivariate gaussian") { + val optimizer = new SpectralProjectedGradient[DenseVector[Double]](tolerance = 1.0E-9) + forAll { init: DenseVector[Double] => + val f = new DiffFunction[DenseVector[Double]] { + def calculate(x: DenseVector[Double]) = { + (sum((x - 3.0) :^ 2.0), (x * 2.0) - 6.0) + } + } + + val result = optimizer.minimize(f, init) + result should beSimilarTo(DenseVector.fill(result.size)(3.0), allowedDeviation = 1E-5) + } + } + + property("optimize a simple multivariate gaussian with projection") { + val optimizer = new SpectralProjectedGradient[DenseVector[Double]](tolerance = 1.0E-5, projection = _.map(scala.math.min(_, 2.0))) + + forAll { init: DenseVector[Double] => + init := clip(init, Double.NegativeInfinity, 2.0) + val f = new DiffFunction[DenseVector[Double]] { + def calculate(x: DenseVector[Double]) = { + (sum((x - 3.0) :^ 4.0), (x - 3.0) :^ 3.0 :* 4.0) + } + } + + val result = optimizer.minimize(f, init) + result should beSimilarTo(DenseVector.fill(result.size)(2.0), allowedDeviation = 1E-10) + } + } + + property("simple linear solve without projection") { + val n = 5 + val H = new DenseMatrix(n, n, Array(1.8984250861699135,0.5955576666769438,-1.484430453342902,-1.0434994471390804,-3.675310432634351,0.5955576666769438,0.9090751938470876,-2.146380947361661,-0.13037609428980368,-0.40639564652095117,-1.484430453342902,-2.146380947361661,10.262733520770384,-6.097698907163584,2.29625304115155,-1.0434994471390804,-0.13037609428980368,-6.097698907163584,27.775920405610677,-5.574220233644466,-3.675310432634351,-0.40639564652095117,2.29625304115155,-5.574220233644466,12.21329172136971)) + val f = DenseVector(-1.2320199653150048, -0.14220655875869606, 0.38477404739124765, -0.3480575854151014, -0.4729810900829228) + + val cost = QuadraticMinimizer.Cost(H, f:*(-1.0)) + val init = DenseVector.zeros[Double](n) + + init := 0.0 + val x = H \ f + + init := 0.0 + val spgResult = new SpectralProjectedGradient[DenseVector[Double]]().minimize(cost, init) + assert(norm(x - spgResult, inf) < 1e-4, s"$x $spgResult") + } + + property("simple linear solve with L1 projection") { + val Hl1 = new DenseMatrix(25, 25, Array(253.535098, 236.477785, 234.421906, 223.374867, 241.007512, 204.695511, 226.465507, 223.351032, 249.179386, 221.411909, 238.679352, 203.579010, 217.564498, 243.852681, 266.607649, 213.994496, 241.620759, 223.602907, 220.038678, 264.704959, 240.341716, 223.672378, 244.323303, 223.216217, 226.074990, 236.477785, 278.862035, 245.756639, 237.489890, 252.783139, 214.077652, 241.816953, 238.790633, 260.536460, 228.121417, 255.103936, 216.608405, 237.150426, 258.933231, 281.958112, 228.971242, 252.508513, 234.242638, 240.308477, 285.416390, 254.792243, 240.176223, 259.048267, 235.566855, 236.277617, 234.421906, 245.756639, 269.162882, 231.416867, 251.321527, 208.134322, 236.567647, 236.558029, 255.805108, 226.535825, 251.514713, 212.770208, 228.565362, 261.748652, 273.946966, 227.411615, 252.767900, 232.823977, 233.084574, 278.315614, 250.872786, 235.227909, 255.104263, 238.931093, 235.402356, 223.374867, 237.489890, 231.416867, 254.771963, 241.703229, 209.028084, 231.517998, 228.768510, 250.805315, 216.548935, 245.473869, 207.687875, 222.036114, 250.906955, 263.018181, 216.128966, 244.445283, 227.436840, 231.369510, 270.721492, 242.475130, 226.471530, 248.130112, 225.826557, 228.266719, 241.007512, 252.783139, 251.321527, 241.703229, 285.702320, 219.051868, 249.442308, 240.400187, 264.970407, 232.503138, 258.819837, 220.160683, 235.621356, 267.743972, 285.795029, 229.667231, 260.870105, 240.751687, 247.183922, 289.044453, 260.715749, 244.210258, 267.159502, 242.992822, 244.070245, 204.695511, 214.077652, 208.134322, 209.028084, 219.051868, 210.164224, 208.151908, 201.539036, 226.373834, 192.056565, 219.950686, 191.459568, 195.982460, 226.739575, 240.677519, 196.116652, 217.352348, 203.533069, 204.581690, 243.603643, 217.785986, 204.205559, 223.747953, 203.586842, 200.165867, 226.465507, 241.816953, 236.567647, 231.517998, 249.442308, 208.151908, 264.007925, 227.080718, 253.174653, 220.322823, 248.619983, 210.100242, 223.279198, 254.807401, 269.896959, 222.927882, 247.017507, 230.484479, 233.358639, 274.935489, 249.237737, 235.229584, 253.029955, 228.601700, 230.512885, 223.351032, 238.790633, 236.558029, 228.768510, 240.400187, 201.539036, 227.080718, 258.773479, 249.471480, 215.664539, 243.078577, 202.337063, 221.020998, 249.979759, 263.356244, 213.470569, 246.182278, 225.727773, 229.873732, 266.295057, 242.954024, 225.510760, 249.370268, 227.813265, 232.141964, 249.179386, 260.536460, 255.805108, 250.805315, 264.970407, 226.373834, 253.174653, 249.471480, 302.360150, 237.902729, 265.769812, 224.947876, 243.088105, 273.690377, 291.076027, 241.089661, 267.772651, 248.459822, 249.662698, 295.935799, 267.460908, 255.668926, 275.902272, 248.495606, 246.827505, 221.411909, 228.121417, 226.535825, 216.548935, 232.503138, 192.056565, 220.322823, 215.664539, 237.902729, 245.154567, 234.956316, 199.557862, 214.774631, 240.339217, 255.161923, 209.328714, 232.277540, 216.298768, 220.296241, 253.817633, 237.638235, 220.785141, 239.098500, 220.583355, 218.962732, 238.679352, 255.103936, 251.514713, 245.473869, 258.819837, 219.950686, 248.619983, 243.078577, 265.769812, 234.956316, 288.133073, 225.087852, 239.810430, 268.406605, 283.289840, 233.858455, 258.306589, 240.263617, 246.844456, 290.492875, 267.212598, 243.218596, 265.681905, 244.615890, 242.543363, 203.579010, 216.608405, 212.770208, 207.687875, 220.160683, 191.459568, 210.100242, 202.337063, 224.947876, 199.557862, 225.087852, 217.501685, 197.897572, 229.825316, 242.175607, 201.123644, 219.820165, 202.894307, 211.468055, 246.048907, 225.135194, 210.076305, 226.806762, 212.014431, 205.123267, 217.564498, 237.150426, 228.565362, 222.036114, 235.621356, 195.982460, 223.279198, 221.020998, 243.088105, 214.774631, 239.810430, 197.897572, 244.439113, 241.621129, 260.400953, 216.482178, 236.805076, 216.680343, 223.816297, 263.188711, 236.311810, 222.950152, 244.636356, 219.121372, 219.911078, 243.852681, 258.933231, 261.748652, 250.906955, 267.743972, 226.739575, 254.807401, 249.979759, 273.690377, 240.339217, 268.406605, 229.825316, 241.621129, 302.928261, 288.344398, 238.549018, 267.239982, 248.073140, 254.230916, 296.789984, 267.158551, 252.226496, 271.170860, 248.325354, 253.694013, 266.607649, 281.958112, 273.946966, 263.018181, 285.795029, 240.677519, 269.896959, 263.356244, 291.076027, 255.161923, 283.289840, 242.175607, 260.400953, 288.344398, 343.457361, 257.368309, 284.795470, 263.122266, 271.239770, 320.209823, 283.933299, 264.416752, 292.035194, 268.764031, 265.345807, 213.994496, 228.971242, 227.411615, 216.128966, 229.667231, 196.116652, 222.927882, 213.470569, 241.089661, 209.328714, 233.858455, 201.123644, 216.482178, 238.549018, 257.368309, 239.295031, 234.913508, 218.066855, 219.648997, 257.969951, 231.243624, 224.657569, 238.158714, 217.174368, 215.933866, 241.620759, 252.508513, 252.767900, 244.445283, 260.870105, 217.352348, 247.017507, 246.182278, 267.772651, 232.277540, 258.306589, 219.820165, 236.805076, 267.239982, 284.795470, 234.913508, 289.709239, 241.312315, 247.249491, 286.702147, 264.252852, 245.151647, 264.582984, 240.842689, 245.837476, 223.602907, 234.242638, 232.823977, 227.436840, 240.751687, 203.533069, 230.484479, 225.727773, 248.459822, 216.298768, 240.263617, 202.894307, 216.680343, 248.073140, 263.122266, 218.066855, 241.312315, 255.363057, 230.209787, 271.091482, 239.220241, 225.387834, 247.486715, 226.052431, 224.119935, 220.038678, 240.308477, 233.084574, 231.369510, 247.183922, 204.581690, 233.358639, 229.873732, 249.662698, 220.296241, 246.844456, 211.468055, 223.816297, 254.230916, 271.239770, 219.648997, 247.249491, 230.209787, 264.014907, 271.938970, 246.664305, 227.889045, 249.908085, 232.035369, 229.010298, 264.704959, 285.416390, 278.315614, 270.721492, 289.044453, 243.603643, 274.935489, 266.295057, 295.935799, 253.817633, 290.492875, 246.048907, 263.188711, 296.789984, 320.209823, 257.969951, 286.702147, 271.091482, 271.938970, 352.825726, 286.200221, 267.716897, 297.182554, 269.776351, 266.721561, 240.341716, 254.792243, 250.872786, 242.475130, 260.715749, 217.785986, 249.237737, 242.954024, 267.460908, 237.638235, 267.212598, 225.135194, 236.311810, 267.158551, 283.933299, 231.243624, 264.252852, 239.220241, 246.664305, 286.200221, 294.042749, 246.504021, 269.570596, 243.980697, 242.690997, 223.672378, 240.176223, 235.227909, 226.471530, 244.210258, 204.205559, 235.229584, 225.510760, 255.668926, 220.785141, 243.218596, 210.076305, 222.950152, 252.226496, 264.416752, 224.657569, 245.151647, 225.387834, 227.889045, 267.716897, 246.504021, 259.897656, 251.730847, 229.335712, 229.759185, 244.323303, 259.048267, 255.104263, 248.130112, 267.159502, 223.747953, 253.029955, 249.370268, 275.902272, 239.098500, 265.681905, 226.806762, 244.636356, 271.170860, 292.035194, 238.158714, 264.582984, 247.486715, 249.908085, 297.182554, 269.570596, 251.730847, 303.872223, 251.585636, 247.878402, 223.216217, 235.566855, 238.931093, 225.826557, 242.992822, 203.586842, 228.601700, 227.813265, 248.495606, 220.583355, 244.615890, 212.014431, 219.121372, 248.325354, 268.764031, 217.174368, 240.842689, 226.052431, 232.035369, 269.776351, 243.980697, 229.335712, 251.585636, 257.544914, 228.810942, 226.074990, 236.277617, 235.402356, 228.266719, 244.070245, 200.165867, 230.512885, 232.141964, 246.827505, 218.962732, 242.543363, 205.123267, 219.911078, 253.694013, 265.345807, 215.933866, 245.837476, 224.119935, 229.010298, 266.721561, 242.690997, 229.759185, 247.878402, 228.810942, 253.353769)) + val fl1 = DenseVector(-892.842851, -934.071560, -932.936015, -888.124343, -961.050207, -791.191087, -923.711397, -904.289301, -988.384984, -883.909133, -959.465030, -798.551172, -871.622303, -997.463289, -1043.912620, -863.013719, -976.975712, -897.033693, -898.694786, -1069.245497, -963.491924, -901.263474, -983.768031, -899.865392, -902.283567) + val cost = QuadraticMinimizer.Cost(Hl1, fl1) + + val octaveL1 = DenseVector(0.18611, 0.00000, 0.06317, -0.10417, 0.11262, + -0.20495, 0.52668, 0.32790, 0.19421, 0.72180, + 0.06309, -0.41326, -0.00000, 0.52078, -0.00000, + 0.18040, 0.62915, 0.16329, -0.06424, 0.37539, + 0.01659, 0.00000, 0.11215, 0.24778, 0.04082) + + val s = octaveL1.foldLeft(0.0) { case (agg, entry) => agg + abs(entry)} + val projectL1 = ProjectL1(s) + val spgResult = new SpectralProjectedGradient[DenseVector[Double]](Projection(projectL1).project).minimizeAndReturnState(cost, DenseVector.zeros[Double](25)) + println(s"SpectralProjectedGradient L1 projection iters ${spgResult.iter}") + assert(norm(spgResult.x - octaveL1, 2) < 1e-4) + } +} diff --git a/math/src/test/scala/breeze/optimize/proximal/NonlinearMinimizerTest.scala b/math/src/test/scala/breeze/optimize/proximal/NonlinearMinimizerTest.scala new file mode 100644 index 000000000..55d861cc1 --- /dev/null +++ b/math/src/test/scala/breeze/optimize/proximal/NonlinearMinimizerTest.scala @@ -0,0 +1,128 @@ +package breeze.optimize.proximal + +import breeze.linalg.{sum, DenseMatrix, norm, DenseVector} +import breeze.optimize.proximal.Constraint._ +import breeze.optimize.OptimizeTestBase +import org.scalatest.junit.JUnitRunner +import org.junit.runner.RunWith +import breeze.numerics._ + +import org.scalatest.Matchers + +@RunWith(classOf[JUnitRunner]) +class NonlinearMinimizerTest extends OptimizeTestBase with Matchers { + val n = 5 + val H = new DenseMatrix(n, n, Array(1.8984250861699135,0.5955576666769438,-1.484430453342902,-1.0434994471390804,-3.675310432634351,0.5955576666769438,0.9090751938470876,-2.146380947361661,-0.13037609428980368,-0.40639564652095117,-1.484430453342902,-2.146380947361661,10.262733520770384,-6.097698907163584,2.29625304115155,-1.0434994471390804,-0.13037609428980368,-6.097698907163584,27.775920405610677,-5.574220233644466,-3.675310432634351,-0.40639564652095117,2.29625304115155,-5.574220233644466,12.21329172136971)) + val f = DenseVector(-1.2320199653150048, -0.14220655875869606, 0.38477404739124765, -0.3480575854151014, -0.4729810900829228) + + val cost = QuadraticMinimizer.Cost(H, f:*(-1.0)) + val init = DenseVector.zeros[Double](n) + + test("Nonlinear Minimization with Identity constraint compared to BFGS") { + init := 0.0 + val x = H \ f + + init := 0.0 + val nlResult = NonlinearMinimizer(n,IDENTITY,0.0).minimize(cost, init) + assert(norm(x - nlResult, inf) < 1e-4) + } + + test("Nonlinear Minimization with positivity constraint compared to Octave") { + val ata = new DenseMatrix[Double](5, 5, Array( + 4.377, -3.531, -1.306, -0.139, 3.418, + -3.531, 4.344, 0.934, 0.305, -2.140, + -1.306, 0.934, 2.644, -0.203, -0.170, + -0.139, 0.305, -0.203, 5.883, 1.428, + 3.418, -2.140, -0.170, 1.428, 4.684)) + + val atb = DenseVector(-1.632, 2.115, 1.094, -1.025, -0.636) + val goodx = DenseVector(0.13025, 0.54506, 0.2874, 0.0, 0.028628) + val nlResult = NonlinearMinimizer(n,POSITIVE,0.0).minimizeAndReturnState(QuadraticMinimizer.Cost(ata, atb:*(-1.0)), DenseVector.zeros[Double](n)) + println(s"Positivity projection iter ${nlResult.iter}") + assert(norm(nlResult.x - goodx, 2) < 1E-3) + } + + test("Nonlinear Minimization with positivity proximal compared to Octave") { + val ata = new DenseMatrix[Double](5, 5, Array( + 4.377, -3.531, -1.306, -0.139, 3.418, + -3.531, 4.344, 0.934, 0.305, -2.140, + -1.306, 0.934, 2.644, -0.203, -0.170, + -0.139, 0.305, -0.203, 5.883, 1.428, + 3.418, -2.140, -0.170, 1.428, 4.684)) + + val atb = DenseVector(-1.632, 2.115, 1.094, -1.025, -0.636) + val goodx = DenseVector(0.13025, 0.54506, 0.2874, 0.0, 0.028628) + val nlResult = new NonlinearMinimizer(ProjectPos()).minimizeAndReturnState(QuadraticMinimizer.Cost(ata, atb:*(-1.0)), DenseVector.zeros[Double](n)) + println(s"Positivity proximal iter ${nlResult.iter}") + assert(norm(nlResult.z - goodx, 2) < 1E-3) + } + + test("Nonlinear Minimization with bounds constraint compared to QuadraticMinimizer") { + init := 0.0 + val gold = QuadraticMinimizer(n, BOX).minimize(H, f:*(-1.0)) + val nlResult = NonlinearMinimizer(n, BOX, 0.0).minimizeAndReturnState(cost, init) + println(s"Bounds projection iter ${nlResult.iter}") + assert(norm(nlResult.x - gold) < 1E-4) + } + + test("Nonlinear Minimization with bounds proximal compared to QuadraticMinimizer") { + init := 0.0 + val gold = QuadraticMinimizer(n, BOX).minimize(H, f:*(-1.0)) + val lb = DenseVector.zeros[Double](n) + val ub = DenseVector.ones[Double](n) + val nlResult = new NonlinearMinimizer(ProjectBox(lb, ub)).minimizeAndReturnState(cost, init) + println(s"Bounds proximal iter ${nlResult.iter}") + assert(norm(nlResult.z - gold) < 1E-4) + } + + test("Nonlinear Minimization with probability simplex compared to Octave") { + val Hml = new DenseMatrix(25, 25, Array(112.647378, 44.962984, 49.127829, 43.708389, 43.333008, 46.345827, 49.581542, 42.991226, 43.999341, 41.336724, 42.879992, 46.896465, 41.778920, 46.333559, 51.168782, 44.800998, 43.735417, 42.672057, 40.024492, 48.793499, 48.696170, 45.870016, 46.398093, 44.305791, 41.863013, 44.962984, 107.202825, 44.218178, 38.585858, 36.606830, 41.783275, 44.631314, 40.883821, 37.948817, 34.908843, 38.356328, 43.642467, 36.213124, 38.760314, 43.317775, 36.803445, 41.905953, 40.238334, 42.325769, 45.853665, 46.601722, 40.668861, 49.084078, 39.292553, 35.781804, 49.127829, 44.218178, 118.264304, 40.139032, 43.741591, 49.387932, 45.558785, 40.793703, 46.136010, 41.839393, 39.544248, 43.161644, 43.361811, 43.832852, 50.572459, 42.036961, 47.251940, 45.273068, 42.842437, 49.323737, 52.125739, 45.831747, 49.466716, 44.762183, 41.930313, 43.708389, 38.585858, 40.139032, 94.937989, 36.562570, 41.628404, 38.604965, 39.080500, 37.267530, 34.291272, 34.891704, 39.216238, 35.970491, 40.733288, 41.872521, 35.825264, 38.144457, 41.368293, 40.751910, 41.123673, 41.930358, 41.002915, 43.099168, 36.018699, 33.646602, 43.333008, 36.606830, 43.741591, 36.562570, 105.764912, 42.799031, 38.215171, 42.193565, 38.708056, 39.448031, 37.882184, 40.172339, 40.625192, 39.015338, 36.433413, 40.848178, 36.480813, 41.439981, 40.797598, 40.325652, 38.599119, 42.727171, 39.382845, 41.535989, 41.518779, 46.345827, 41.783275, 49.387932, 41.628404, 42.799031, 114.691992, 43.015599, 42.688570, 42.722905, 38.675192, 38.377970, 44.656183, 39.087805, 45.443516, 50.585268, 40.949970, 41.920556, 43.711898, 41.463472, 51.248836, 46.869144, 45.178199, 45.709593, 42.402465, 44.097412, 49.581542, 44.631314, 45.558785, 38.604965, 38.215171, 43.015599, 114.667896, 40.966284, 37.748084, 39.496813, 40.534741, 42.770125, 40.628678, 41.194251, 47.837969, 44.596875, 43.448257, 43.291878, 39.005953, 50.493111, 46.296591, 43.449036, 48.798961, 42.877859, 37.055014, 42.991226, 40.883821, 40.793703, 39.080500, 42.193565, 42.688570, 40.966284, 106.632656, 37.640927, 37.181799, 40.065085, 38.978761, 36.014753, 38.071494, 41.081064, 37.981693, 41.821252, 42.773603, 39.293957, 38.600491, 43.761301, 42.294750, 42.410289, 40.266469, 39.909538, 43.999341, 37.948817, 46.136010, 37.267530, 38.708056, 42.722905, 37.748084, 37.640927, 102.747189, 34.574727, 36.525241, 39.839891, 36.297838, 42.756496, 44.673874, 38.350523, 40.330611, 42.288511, 39.472844, 45.617102, 44.692618, 41.194977, 41.284030, 39.580938, 42.382268, 41.336724, 34.908843, 41.839393, 34.291272, 39.448031, 38.675192, 39.496813, 37.181799, 34.574727, 94.205550, 37.583319, 38.504211, 36.376976, 34.239351, 39.060978, 37.515228, 37.079566, 37.317791, 38.046576, 36.112222, 39.406838, 39.258432, 36.347136, 38.927619, 41.604838, 42.879992, 38.356328, 39.544248, 34.891704, 37.882184, 38.377970, 40.534741, 40.065085, 36.525241, 37.583319, 98.109622, 39.428284, 37.518381, 39.659011, 38.477483, 40.547021, 42.678061, 42.279104, 41.515782, 43.478416, 45.003800, 42.433639, 42.757336, 35.814356, 39.017848, 46.896465, 43.642467, 43.161644, 39.216238, 40.172339, 44.656183, 42.770125, 38.978761, 39.839891, 38.504211, 39.428284, 103.478446, 39.984358, 40.587958, 44.490750, 40.600474, 40.698368, 42.296794, 41.567854, 47.002908, 43.922434, 43.479144, 44.291425, 43.352951, 42.613649, 41.778920, 36.213124, 43.361811, 35.970491, 40.625192, 39.087805, 40.628678, 36.014753, 36.297838, 36.376976, 37.518381, 39.984358, 99.799628, 38.027891, 44.268308, 36.202204, 39.921811, 38.668774, 36.832286, 45.833218, 43.228963, 36.833273, 44.787401, 38.176476, 39.062471, 46.333559, 38.760314, 43.832852, 40.733288, 39.015338, 45.443516, 41.194251, 38.071494, 42.756496, 34.239351, 39.659011, 40.587958, 38.027891, 114.304283, 46.958354, 39.636801, 40.927870, 49.118094, 43.093642, 50.196436, 45.535041, 43.087415, 48.540036, 35.942528, 37.962886, 51.168782, 43.317775, 50.572459, 41.872521, 36.433413, 50.585268, 47.837969, 41.081064, 44.673874, 39.060978, 38.477483, 44.490750, 44.268308, 46.958354, 122.935323, 39.948695, 46.801841, 44.455283, 40.160668, 54.193098, 49.678271, 41.834745, 47.227606, 42.214571, 42.598524, 44.800998, 36.803445, 42.036961, 35.825264, 40.848178, 40.949970, 44.596875, 37.981693, 38.350523, 37.515228, 40.547021, 40.600474, 36.202204, 39.636801, 39.948695, 97.365126, 40.163209, 39.177628, 38.935283, 41.465246, 40.962743, 40.533287, 43.367907, 38.723316, 36.312733, 43.735417, 41.905953, 47.251940, 38.144457, 36.480813, 41.920556, 43.448257, 41.821252, 40.330611, 37.079566, 42.678061, 40.698368, 39.921811, 40.927870, 46.801841, 40.163209, 110.416786, 46.843429, 41.834126, 46.788801, 46.983780, 43.511429, 47.291825, 40.023523, 40.581819, 42.672057, 40.238334, 45.273068, 41.368293, 41.439981, 43.711898, 43.291878, 42.773603, 42.288511, 37.317791, 42.279104, 42.296794, 38.668774, 49.118094, 44.455283, 39.177628, 46.843429, 107.474576, 44.590023, 48.333476, 44.059916, 42.653703, 44.171623, 39.363181, 41.716539, 40.024492, 42.325769, 42.842437, 40.751910, 40.797598, 41.463472, 39.005953, 39.293957, 39.472844, 38.046576, 41.515782, 41.567854, 36.832286, 43.093642, 40.160668, 38.935283, 41.834126, 44.590023, 105.140579, 43.149105, 41.516560, 43.494333, 45.664210, 36.466241, 37.477898, 48.793499, 45.853665, 49.323737, 41.123673, 40.325652, 51.248836, 50.493111, 38.600491, 45.617102, 36.112222, 43.478416, 47.002908, 45.833218, 50.196436, 54.193098, 41.465246, 46.788801, 48.333476, 43.149105, 123.746816, 53.234332, 44.633908, 53.537592, 43.196327, 42.747181, 48.696170, 46.601722, 52.125739, 41.930358, 38.599119, 46.869144, 46.296591, 43.761301, 44.692618, 39.406838, 45.003800, 43.922434, 43.228963, 45.535041, 49.678271, 40.962743, 46.983780, 44.059916, 41.516560, 53.234332, 125.202062, 43.967875, 52.416619, 39.937196, 39.775405, 45.870016, 40.668861, 45.831747, 41.002915, 42.727171, 45.178199, 43.449036, 42.294750, 41.194977, 39.258432, 42.433639, 43.479144, 36.833273, 43.087415, 41.834745, 40.533287, 43.511429, 42.653703, 43.494333, 44.633908, 43.967875, 107.336922, 44.396001, 39.819884, 38.676633, 46.398093, 49.084078, 49.466716, 43.099168, 39.382845, 45.709593, 48.798961, 42.410289, 41.284030, 36.347136, 42.757336, 44.291425, 44.787401, 48.540036, 47.227606, 43.367907, 47.291825, 44.171623, 45.664210, 53.537592, 52.416619, 44.396001, 114.651847, 40.826050, 37.634130, 44.305791, 39.292553, 44.762183, 36.018699, 41.535989, 42.402465, 42.877859, 40.266469, 39.580938, 38.927619, 35.814356, 43.352951, 38.176476, 35.942528, 42.214571, 38.723316, 40.023523, 39.363181, 36.466241, 43.196327, 39.937196, 39.819884, 40.826050, 96.128345, 40.788606, 41.863013, 35.781804, 41.930313, 33.646602, 41.518779, 44.097412, 37.055014, 39.909538, 42.382268, 41.604838, 39.017848, 42.613649, 39.062471, 37.962886, 42.598524, 36.312733, 40.581819, 41.716539, 37.477898, 42.747181, 39.775405, 38.676633, 37.634130, 40.788606, 97.605849)) + val fml = DenseVector(-1219.296604, -1126.029219, -1202.257728, -1022.064083, -1047.414836, -1155.507387, -1169.502847, -1091.655366, -1063.832607, -1015.829142, -1075.864072, -1101.427162, -1058.907539, -1115.171116, -1205.015211, -1090.627084, -1143.206126, -1140.107801, -1100.285642, -1198.992795, -1246.276120, -1159.678276, -1194.177391, -1056.458015, -1058.791892) + val cost = QuadraticMinimizer.Cost(Hml, fml) + + val nlResult = NonlinearMinimizer(25, PROBABILITYSIMPLEX, 1.0).minimizeAndReturnState(cost, DenseVector.zeros[Double](25)) + + val golden = DenseVector(0.3131862265452959, 0.0, 0.01129486116330884, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.060642310566736704, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6151756449091074, 0.0, 0.0, 0.0, 0.0) + + println(s"Probability Simplex projection iter ${nlResult.iter}") + + assert(norm(nlResult.x - golden) < 1e-3) + assert(abs(sum(nlResult.x) - 1.0) < 1e-4) + } + + test("Nonlinear Minimization with L1 projection compared to Octave") { + val Hl1 = new DenseMatrix(25, 25, Array(253.535098, 236.477785, 234.421906, 223.374867, 241.007512, 204.695511, 226.465507, 223.351032, 249.179386, 221.411909, 238.679352, 203.579010, 217.564498, 243.852681, 266.607649, 213.994496, 241.620759, 223.602907, 220.038678, 264.704959, 240.341716, 223.672378, 244.323303, 223.216217, 226.074990, 236.477785, 278.862035, 245.756639, 237.489890, 252.783139, 214.077652, 241.816953, 238.790633, 260.536460, 228.121417, 255.103936, 216.608405, 237.150426, 258.933231, 281.958112, 228.971242, 252.508513, 234.242638, 240.308477, 285.416390, 254.792243, 240.176223, 259.048267, 235.566855, 236.277617, 234.421906, 245.756639, 269.162882, 231.416867, 251.321527, 208.134322, 236.567647, 236.558029, 255.805108, 226.535825, 251.514713, 212.770208, 228.565362, 261.748652, 273.946966, 227.411615, 252.767900, 232.823977, 233.084574, 278.315614, 250.872786, 235.227909, 255.104263, 238.931093, 235.402356, 223.374867, 237.489890, 231.416867, 254.771963, 241.703229, 209.028084, 231.517998, 228.768510, 250.805315, 216.548935, 245.473869, 207.687875, 222.036114, 250.906955, 263.018181, 216.128966, 244.445283, 227.436840, 231.369510, 270.721492, 242.475130, 226.471530, 248.130112, 225.826557, 228.266719, 241.007512, 252.783139, 251.321527, 241.703229, 285.702320, 219.051868, 249.442308, 240.400187, 264.970407, 232.503138, 258.819837, 220.160683, 235.621356, 267.743972, 285.795029, 229.667231, 260.870105, 240.751687, 247.183922, 289.044453, 260.715749, 244.210258, 267.159502, 242.992822, 244.070245, 204.695511, 214.077652, 208.134322, 209.028084, 219.051868, 210.164224, 208.151908, 201.539036, 226.373834, 192.056565, 219.950686, 191.459568, 195.982460, 226.739575, 240.677519, 196.116652, 217.352348, 203.533069, 204.581690, 243.603643, 217.785986, 204.205559, 223.747953, 203.586842, 200.165867, 226.465507, 241.816953, 236.567647, 231.517998, 249.442308, 208.151908, 264.007925, 227.080718, 253.174653, 220.322823, 248.619983, 210.100242, 223.279198, 254.807401, 269.896959, 222.927882, 247.017507, 230.484479, 233.358639, 274.935489, 249.237737, 235.229584, 253.029955, 228.601700, 230.512885, 223.351032, 238.790633, 236.558029, 228.768510, 240.400187, 201.539036, 227.080718, 258.773479, 249.471480, 215.664539, 243.078577, 202.337063, 221.020998, 249.979759, 263.356244, 213.470569, 246.182278, 225.727773, 229.873732, 266.295057, 242.954024, 225.510760, 249.370268, 227.813265, 232.141964, 249.179386, 260.536460, 255.805108, 250.805315, 264.970407, 226.373834, 253.174653, 249.471480, 302.360150, 237.902729, 265.769812, 224.947876, 243.088105, 273.690377, 291.076027, 241.089661, 267.772651, 248.459822, 249.662698, 295.935799, 267.460908, 255.668926, 275.902272, 248.495606, 246.827505, 221.411909, 228.121417, 226.535825, 216.548935, 232.503138, 192.056565, 220.322823, 215.664539, 237.902729, 245.154567, 234.956316, 199.557862, 214.774631, 240.339217, 255.161923, 209.328714, 232.277540, 216.298768, 220.296241, 253.817633, 237.638235, 220.785141, 239.098500, 220.583355, 218.962732, 238.679352, 255.103936, 251.514713, 245.473869, 258.819837, 219.950686, 248.619983, 243.078577, 265.769812, 234.956316, 288.133073, 225.087852, 239.810430, 268.406605, 283.289840, 233.858455, 258.306589, 240.263617, 246.844456, 290.492875, 267.212598, 243.218596, 265.681905, 244.615890, 242.543363, 203.579010, 216.608405, 212.770208, 207.687875, 220.160683, 191.459568, 210.100242, 202.337063, 224.947876, 199.557862, 225.087852, 217.501685, 197.897572, 229.825316, 242.175607, 201.123644, 219.820165, 202.894307, 211.468055, 246.048907, 225.135194, 210.076305, 226.806762, 212.014431, 205.123267, 217.564498, 237.150426, 228.565362, 222.036114, 235.621356, 195.982460, 223.279198, 221.020998, 243.088105, 214.774631, 239.810430, 197.897572, 244.439113, 241.621129, 260.400953, 216.482178, 236.805076, 216.680343, 223.816297, 263.188711, 236.311810, 222.950152, 244.636356, 219.121372, 219.911078, 243.852681, 258.933231, 261.748652, 250.906955, 267.743972, 226.739575, 254.807401, 249.979759, 273.690377, 240.339217, 268.406605, 229.825316, 241.621129, 302.928261, 288.344398, 238.549018, 267.239982, 248.073140, 254.230916, 296.789984, 267.158551, 252.226496, 271.170860, 248.325354, 253.694013, 266.607649, 281.958112, 273.946966, 263.018181, 285.795029, 240.677519, 269.896959, 263.356244, 291.076027, 255.161923, 283.289840, 242.175607, 260.400953, 288.344398, 343.457361, 257.368309, 284.795470, 263.122266, 271.239770, 320.209823, 283.933299, 264.416752, 292.035194, 268.764031, 265.345807, 213.994496, 228.971242, 227.411615, 216.128966, 229.667231, 196.116652, 222.927882, 213.470569, 241.089661, 209.328714, 233.858455, 201.123644, 216.482178, 238.549018, 257.368309, 239.295031, 234.913508, 218.066855, 219.648997, 257.969951, 231.243624, 224.657569, 238.158714, 217.174368, 215.933866, 241.620759, 252.508513, 252.767900, 244.445283, 260.870105, 217.352348, 247.017507, 246.182278, 267.772651, 232.277540, 258.306589, 219.820165, 236.805076, 267.239982, 284.795470, 234.913508, 289.709239, 241.312315, 247.249491, 286.702147, 264.252852, 245.151647, 264.582984, 240.842689, 245.837476, 223.602907, 234.242638, 232.823977, 227.436840, 240.751687, 203.533069, 230.484479, 225.727773, 248.459822, 216.298768, 240.263617, 202.894307, 216.680343, 248.073140, 263.122266, 218.066855, 241.312315, 255.363057, 230.209787, 271.091482, 239.220241, 225.387834, 247.486715, 226.052431, 224.119935, 220.038678, 240.308477, 233.084574, 231.369510, 247.183922, 204.581690, 233.358639, 229.873732, 249.662698, 220.296241, 246.844456, 211.468055, 223.816297, 254.230916, 271.239770, 219.648997, 247.249491, 230.209787, 264.014907, 271.938970, 246.664305, 227.889045, 249.908085, 232.035369, 229.010298, 264.704959, 285.416390, 278.315614, 270.721492, 289.044453, 243.603643, 274.935489, 266.295057, 295.935799, 253.817633, 290.492875, 246.048907, 263.188711, 296.789984, 320.209823, 257.969951, 286.702147, 271.091482, 271.938970, 352.825726, 286.200221, 267.716897, 297.182554, 269.776351, 266.721561, 240.341716, 254.792243, 250.872786, 242.475130, 260.715749, 217.785986, 249.237737, 242.954024, 267.460908, 237.638235, 267.212598, 225.135194, 236.311810, 267.158551, 283.933299, 231.243624, 264.252852, 239.220241, 246.664305, 286.200221, 294.042749, 246.504021, 269.570596, 243.980697, 242.690997, 223.672378, 240.176223, 235.227909, 226.471530, 244.210258, 204.205559, 235.229584, 225.510760, 255.668926, 220.785141, 243.218596, 210.076305, 222.950152, 252.226496, 264.416752, 224.657569, 245.151647, 225.387834, 227.889045, 267.716897, 246.504021, 259.897656, 251.730847, 229.335712, 229.759185, 244.323303, 259.048267, 255.104263, 248.130112, 267.159502, 223.747953, 253.029955, 249.370268, 275.902272, 239.098500, 265.681905, 226.806762, 244.636356, 271.170860, 292.035194, 238.158714, 264.582984, 247.486715, 249.908085, 297.182554, 269.570596, 251.730847, 303.872223, 251.585636, 247.878402, 223.216217, 235.566855, 238.931093, 225.826557, 242.992822, 203.586842, 228.601700, 227.813265, 248.495606, 220.583355, 244.615890, 212.014431, 219.121372, 248.325354, 268.764031, 217.174368, 240.842689, 226.052431, 232.035369, 269.776351, 243.980697, 229.335712, 251.585636, 257.544914, 228.810942, 226.074990, 236.277617, 235.402356, 228.266719, 244.070245, 200.165867, 230.512885, 232.141964, 246.827505, 218.962732, 242.543363, 205.123267, 219.911078, 253.694013, 265.345807, 215.933866, 245.837476, 224.119935, 229.010298, 266.721561, 242.690997, 229.759185, 247.878402, 228.810942, 253.353769)) + val fl1 = DenseVector(-892.842851, -934.071560, -932.936015, -888.124343, -961.050207, -791.191087, -923.711397, -904.289301, -988.384984, -883.909133, -959.465030, -798.551172, -871.622303, -997.463289, -1043.912620, -863.013719, -976.975712, -897.033693, -898.694786, -1069.245497, -963.491924, -901.263474, -983.768031, -899.865392, -902.283567) + val cost = QuadraticMinimizer.Cost(Hl1, fl1) + + val octaveL1 = DenseVector(0.18611, 0.00000, 0.06317, -0.10417, 0.11262, + -0.20495, 0.52668, 0.32790, 0.19421, 0.72180, + 0.06309, -0.41326, -0.00000, 0.52078, -0.00000, + 0.18040, 0.62915, 0.16329, -0.06424, 0.37539, + 0.01659, 0.00000, 0.11215, 0.24778, 0.04082) + + val s = octaveL1.foldLeft(0.0) { case (agg, entry) => agg + abs(entry)} + val nlResult = NonlinearMinimizer(25, SPARSE, s).minimizeAndReturnState(cost, DenseVector.zeros[Double](25)) + + println(s"L1 projection iter ${nlResult.iter}") + assert(norm(nlResult.x - octaveL1, 2) < 1e-4) + } + + test("Nonlinear Minimization with L1 proximal compared to Octave") { + val Hl1 = new DenseMatrix(25, 25, Array(253.535098, 236.477785, 234.421906, 223.374867, 241.007512, 204.695511, 226.465507, 223.351032, 249.179386, 221.411909, 238.679352, 203.579010, 217.564498, 243.852681, 266.607649, 213.994496, 241.620759, 223.602907, 220.038678, 264.704959, 240.341716, 223.672378, 244.323303, 223.216217, 226.074990, 236.477785, 278.862035, 245.756639, 237.489890, 252.783139, 214.077652, 241.816953, 238.790633, 260.536460, 228.121417, 255.103936, 216.608405, 237.150426, 258.933231, 281.958112, 228.971242, 252.508513, 234.242638, 240.308477, 285.416390, 254.792243, 240.176223, 259.048267, 235.566855, 236.277617, 234.421906, 245.756639, 269.162882, 231.416867, 251.321527, 208.134322, 236.567647, 236.558029, 255.805108, 226.535825, 251.514713, 212.770208, 228.565362, 261.748652, 273.946966, 227.411615, 252.767900, 232.823977, 233.084574, 278.315614, 250.872786, 235.227909, 255.104263, 238.931093, 235.402356, 223.374867, 237.489890, 231.416867, 254.771963, 241.703229, 209.028084, 231.517998, 228.768510, 250.805315, 216.548935, 245.473869, 207.687875, 222.036114, 250.906955, 263.018181, 216.128966, 244.445283, 227.436840, 231.369510, 270.721492, 242.475130, 226.471530, 248.130112, 225.826557, 228.266719, 241.007512, 252.783139, 251.321527, 241.703229, 285.702320, 219.051868, 249.442308, 240.400187, 264.970407, 232.503138, 258.819837, 220.160683, 235.621356, 267.743972, 285.795029, 229.667231, 260.870105, 240.751687, 247.183922, 289.044453, 260.715749, 244.210258, 267.159502, 242.992822, 244.070245, 204.695511, 214.077652, 208.134322, 209.028084, 219.051868, 210.164224, 208.151908, 201.539036, 226.373834, 192.056565, 219.950686, 191.459568, 195.982460, 226.739575, 240.677519, 196.116652, 217.352348, 203.533069, 204.581690, 243.603643, 217.785986, 204.205559, 223.747953, 203.586842, 200.165867, 226.465507, 241.816953, 236.567647, 231.517998, 249.442308, 208.151908, 264.007925, 227.080718, 253.174653, 220.322823, 248.619983, 210.100242, 223.279198, 254.807401, 269.896959, 222.927882, 247.017507, 230.484479, 233.358639, 274.935489, 249.237737, 235.229584, 253.029955, 228.601700, 230.512885, 223.351032, 238.790633, 236.558029, 228.768510, 240.400187, 201.539036, 227.080718, 258.773479, 249.471480, 215.664539, 243.078577, 202.337063, 221.020998, 249.979759, 263.356244, 213.470569, 246.182278, 225.727773, 229.873732, 266.295057, 242.954024, 225.510760, 249.370268, 227.813265, 232.141964, 249.179386, 260.536460, 255.805108, 250.805315, 264.970407, 226.373834, 253.174653, 249.471480, 302.360150, 237.902729, 265.769812, 224.947876, 243.088105, 273.690377, 291.076027, 241.089661, 267.772651, 248.459822, 249.662698, 295.935799, 267.460908, 255.668926, 275.902272, 248.495606, 246.827505, 221.411909, 228.121417, 226.535825, 216.548935, 232.503138, 192.056565, 220.322823, 215.664539, 237.902729, 245.154567, 234.956316, 199.557862, 214.774631, 240.339217, 255.161923, 209.328714, 232.277540, 216.298768, 220.296241, 253.817633, 237.638235, 220.785141, 239.098500, 220.583355, 218.962732, 238.679352, 255.103936, 251.514713, 245.473869, 258.819837, 219.950686, 248.619983, 243.078577, 265.769812, 234.956316, 288.133073, 225.087852, 239.810430, 268.406605, 283.289840, 233.858455, 258.306589, 240.263617, 246.844456, 290.492875, 267.212598, 243.218596, 265.681905, 244.615890, 242.543363, 203.579010, 216.608405, 212.770208, 207.687875, 220.160683, 191.459568, 210.100242, 202.337063, 224.947876, 199.557862, 225.087852, 217.501685, 197.897572, 229.825316, 242.175607, 201.123644, 219.820165, 202.894307, 211.468055, 246.048907, 225.135194, 210.076305, 226.806762, 212.014431, 205.123267, 217.564498, 237.150426, 228.565362, 222.036114, 235.621356, 195.982460, 223.279198, 221.020998, 243.088105, 214.774631, 239.810430, 197.897572, 244.439113, 241.621129, 260.400953, 216.482178, 236.805076, 216.680343, 223.816297, 263.188711, 236.311810, 222.950152, 244.636356, 219.121372, 219.911078, 243.852681, 258.933231, 261.748652, 250.906955, 267.743972, 226.739575, 254.807401, 249.979759, 273.690377, 240.339217, 268.406605, 229.825316, 241.621129, 302.928261, 288.344398, 238.549018, 267.239982, 248.073140, 254.230916, 296.789984, 267.158551, 252.226496, 271.170860, 248.325354, 253.694013, 266.607649, 281.958112, 273.946966, 263.018181, 285.795029, 240.677519, 269.896959, 263.356244, 291.076027, 255.161923, 283.289840, 242.175607, 260.400953, 288.344398, 343.457361, 257.368309, 284.795470, 263.122266, 271.239770, 320.209823, 283.933299, 264.416752, 292.035194, 268.764031, 265.345807, 213.994496, 228.971242, 227.411615, 216.128966, 229.667231, 196.116652, 222.927882, 213.470569, 241.089661, 209.328714, 233.858455, 201.123644, 216.482178, 238.549018, 257.368309, 239.295031, 234.913508, 218.066855, 219.648997, 257.969951, 231.243624, 224.657569, 238.158714, 217.174368, 215.933866, 241.620759, 252.508513, 252.767900, 244.445283, 260.870105, 217.352348, 247.017507, 246.182278, 267.772651, 232.277540, 258.306589, 219.820165, 236.805076, 267.239982, 284.795470, 234.913508, 289.709239, 241.312315, 247.249491, 286.702147, 264.252852, 245.151647, 264.582984, 240.842689, 245.837476, 223.602907, 234.242638, 232.823977, 227.436840, 240.751687, 203.533069, 230.484479, 225.727773, 248.459822, 216.298768, 240.263617, 202.894307, 216.680343, 248.073140, 263.122266, 218.066855, 241.312315, 255.363057, 230.209787, 271.091482, 239.220241, 225.387834, 247.486715, 226.052431, 224.119935, 220.038678, 240.308477, 233.084574, 231.369510, 247.183922, 204.581690, 233.358639, 229.873732, 249.662698, 220.296241, 246.844456, 211.468055, 223.816297, 254.230916, 271.239770, 219.648997, 247.249491, 230.209787, 264.014907, 271.938970, 246.664305, 227.889045, 249.908085, 232.035369, 229.010298, 264.704959, 285.416390, 278.315614, 270.721492, 289.044453, 243.603643, 274.935489, 266.295057, 295.935799, 253.817633, 290.492875, 246.048907, 263.188711, 296.789984, 320.209823, 257.969951, 286.702147, 271.091482, 271.938970, 352.825726, 286.200221, 267.716897, 297.182554, 269.776351, 266.721561, 240.341716, 254.792243, 250.872786, 242.475130, 260.715749, 217.785986, 249.237737, 242.954024, 267.460908, 237.638235, 267.212598, 225.135194, 236.311810, 267.158551, 283.933299, 231.243624, 264.252852, 239.220241, 246.664305, 286.200221, 294.042749, 246.504021, 269.570596, 243.980697, 242.690997, 223.672378, 240.176223, 235.227909, 226.471530, 244.210258, 204.205559, 235.229584, 225.510760, 255.668926, 220.785141, 243.218596, 210.076305, 222.950152, 252.226496, 264.416752, 224.657569, 245.151647, 225.387834, 227.889045, 267.716897, 246.504021, 259.897656, 251.730847, 229.335712, 229.759185, 244.323303, 259.048267, 255.104263, 248.130112, 267.159502, 223.747953, 253.029955, 249.370268, 275.902272, 239.098500, 265.681905, 226.806762, 244.636356, 271.170860, 292.035194, 238.158714, 264.582984, 247.486715, 249.908085, 297.182554, 269.570596, 251.730847, 303.872223, 251.585636, 247.878402, 223.216217, 235.566855, 238.931093, 225.826557, 242.992822, 203.586842, 228.601700, 227.813265, 248.495606, 220.583355, 244.615890, 212.014431, 219.121372, 248.325354, 268.764031, 217.174368, 240.842689, 226.052431, 232.035369, 269.776351, 243.980697, 229.335712, 251.585636, 257.544914, 228.810942, 226.074990, 236.277617, 235.402356, 228.266719, 244.070245, 200.165867, 230.512885, 232.141964, 246.827505, 218.962732, 242.543363, 205.123267, 219.911078, 253.694013, 265.345807, 215.933866, 245.837476, 224.119935, 229.010298, 266.721561, 242.690997, 229.759185, 247.878402, 228.810942, 253.353769)) + val fl1 = DenseVector(-892.842851, -934.071560, -932.936015, -888.124343, -961.050207, -791.191087, -923.711397, -904.289301, -988.384984, -883.909133, -959.465030, -798.551172, -871.622303, -997.463289, -1043.912620, -863.013719, -976.975712, -897.033693, -898.694786, -1069.245497, -963.491924, -901.263474, -983.768031, -899.865392, -902.283567) + val cost = QuadraticMinimizer.Cost(Hl1, fl1) + + val octaveL1 = DenseVector(0.18611, 0.00000, 0.06317, -0.10417, 0.11262, + -0.20495, 0.52668, 0.32790, 0.19421, 0.72180, + 0.06309, -0.41326, -0.00000, 0.52078, -0.00000, + 0.18040, 0.62915, 0.16329, -0.06424, 0.37539, + 0.01659, 0.00000, 0.11215, 0.24778, 0.04082) + val proximalL1 = ProximalL1().setLambda(2.0) + val nl = new NonlinearMinimizer(proximalL1) + val nlResult = nl.minimizeAndReturnState(cost, DenseVector.zeros[Double](25)) + + println(s"L1 Proximal iter ${nlResult.iter}") + assert(norm(nlResult.z - octaveL1, Inf) < 1e-4) + } +} diff --git a/math/src/test/scala/breeze/optimize/proximal/QuadraticMinimizerTest.scala b/math/src/test/scala/breeze/optimize/proximal/QuadraticMinimizerTest.scala index ca345f8c1..4fb86791e 100644 --- a/math/src/test/scala/breeze/optimize/proximal/QuadraticMinimizerTest.scala +++ b/math/src/test/scala/breeze/optimize/proximal/QuadraticMinimizerTest.scala @@ -154,10 +154,7 @@ class QuadraticMinimizerTest extends OptimizeTestBase with Matchers { val atb = DenseVector(-1.632, 2.115, 1.094, -1.025, -0.636) - val nnls = new NNLS() - val nnlsBounds = nnls.minimize(ata, atb) val goodx = DenseVector(0.13025, 0.54506, 0.2874, 0.0, 0.028628) - val qpSolverPos = QuadraticMinimizer(n, POSITIVE, 0.0) atb *= -1.0 @@ -204,7 +201,7 @@ class QuadraticMinimizerTest extends OptimizeTestBase with Matchers { assert(abs(sum(directQpResult.x) - 1.0) < 1e-4) } - test("Quadratic Minimization with L1 compared to Octave MovieLens Example") { + test("Quadratic Minimization with L1 Proximal compared to Octave MovieLens Example") { val Hl1 = new DenseMatrix(25, 25, Array(253.535098, 236.477785, 234.421906, 223.374867, 241.007512, 204.695511, 226.465507, 223.351032, 249.179386, 221.411909, 238.679352, 203.579010, 217.564498, 243.852681, 266.607649, 213.994496, 241.620759, 223.602907, 220.038678, 264.704959, 240.341716, 223.672378, 244.323303, 223.216217, 226.074990, 236.477785, 278.862035, 245.756639, 237.489890, 252.783139, 214.077652, 241.816953, 238.790633, 260.536460, 228.121417, 255.103936, 216.608405, 237.150426, 258.933231, 281.958112, 228.971242, 252.508513, 234.242638, 240.308477, 285.416390, 254.792243, 240.176223, 259.048267, 235.566855, 236.277617, 234.421906, 245.756639, 269.162882, 231.416867, 251.321527, 208.134322, 236.567647, 236.558029, 255.805108, 226.535825, 251.514713, 212.770208, 228.565362, 261.748652, 273.946966, 227.411615, 252.767900, 232.823977, 233.084574, 278.315614, 250.872786, 235.227909, 255.104263, 238.931093, 235.402356, 223.374867, 237.489890, 231.416867, 254.771963, 241.703229, 209.028084, 231.517998, 228.768510, 250.805315, 216.548935, 245.473869, 207.687875, 222.036114, 250.906955, 263.018181, 216.128966, 244.445283, 227.436840, 231.369510, 270.721492, 242.475130, 226.471530, 248.130112, 225.826557, 228.266719, 241.007512, 252.783139, 251.321527, 241.703229, 285.702320, 219.051868, 249.442308, 240.400187, 264.970407, 232.503138, 258.819837, 220.160683, 235.621356, 267.743972, 285.795029, 229.667231, 260.870105, 240.751687, 247.183922, 289.044453, 260.715749, 244.210258, 267.159502, 242.992822, 244.070245, 204.695511, 214.077652, 208.134322, 209.028084, 219.051868, 210.164224, 208.151908, 201.539036, 226.373834, 192.056565, 219.950686, 191.459568, 195.982460, 226.739575, 240.677519, 196.116652, 217.352348, 203.533069, 204.581690, 243.603643, 217.785986, 204.205559, 223.747953, 203.586842, 200.165867, 226.465507, 241.816953, 236.567647, 231.517998, 249.442308, 208.151908, 264.007925, 227.080718, 253.174653, 220.322823, 248.619983, 210.100242, 223.279198, 254.807401, 269.896959, 222.927882, 247.017507, 230.484479, 233.358639, 274.935489, 249.237737, 235.229584, 253.029955, 228.601700, 230.512885, 223.351032, 238.790633, 236.558029, 228.768510, 240.400187, 201.539036, 227.080718, 258.773479, 249.471480, 215.664539, 243.078577, 202.337063, 221.020998, 249.979759, 263.356244, 213.470569, 246.182278, 225.727773, 229.873732, 266.295057, 242.954024, 225.510760, 249.370268, 227.813265, 232.141964, 249.179386, 260.536460, 255.805108, 250.805315, 264.970407, 226.373834, 253.174653, 249.471480, 302.360150, 237.902729, 265.769812, 224.947876, 243.088105, 273.690377, 291.076027, 241.089661, 267.772651, 248.459822, 249.662698, 295.935799, 267.460908, 255.668926, 275.902272, 248.495606, 246.827505, 221.411909, 228.121417, 226.535825, 216.548935, 232.503138, 192.056565, 220.322823, 215.664539, 237.902729, 245.154567, 234.956316, 199.557862, 214.774631, 240.339217, 255.161923, 209.328714, 232.277540, 216.298768, 220.296241, 253.817633, 237.638235, 220.785141, 239.098500, 220.583355, 218.962732, 238.679352, 255.103936, 251.514713, 245.473869, 258.819837, 219.950686, 248.619983, 243.078577, 265.769812, 234.956316, 288.133073, 225.087852, 239.810430, 268.406605, 283.289840, 233.858455, 258.306589, 240.263617, 246.844456, 290.492875, 267.212598, 243.218596, 265.681905, 244.615890, 242.543363, 203.579010, 216.608405, 212.770208, 207.687875, 220.160683, 191.459568, 210.100242, 202.337063, 224.947876, 199.557862, 225.087852, 217.501685, 197.897572, 229.825316, 242.175607, 201.123644, 219.820165, 202.894307, 211.468055, 246.048907, 225.135194, 210.076305, 226.806762, 212.014431, 205.123267, 217.564498, 237.150426, 228.565362, 222.036114, 235.621356, 195.982460, 223.279198, 221.020998, 243.088105, 214.774631, 239.810430, 197.897572, 244.439113, 241.621129, 260.400953, 216.482178, 236.805076, 216.680343, 223.816297, 263.188711, 236.311810, 222.950152, 244.636356, 219.121372, 219.911078, 243.852681, 258.933231, 261.748652, 250.906955, 267.743972, 226.739575, 254.807401, 249.979759, 273.690377, 240.339217, 268.406605, 229.825316, 241.621129, 302.928261, 288.344398, 238.549018, 267.239982, 248.073140, 254.230916, 296.789984, 267.158551, 252.226496, 271.170860, 248.325354, 253.694013, 266.607649, 281.958112, 273.946966, 263.018181, 285.795029, 240.677519, 269.896959, 263.356244, 291.076027, 255.161923, 283.289840, 242.175607, 260.400953, 288.344398, 343.457361, 257.368309, 284.795470, 263.122266, 271.239770, 320.209823, 283.933299, 264.416752, 292.035194, 268.764031, 265.345807, 213.994496, 228.971242, 227.411615, 216.128966, 229.667231, 196.116652, 222.927882, 213.470569, 241.089661, 209.328714, 233.858455, 201.123644, 216.482178, 238.549018, 257.368309, 239.295031, 234.913508, 218.066855, 219.648997, 257.969951, 231.243624, 224.657569, 238.158714, 217.174368, 215.933866, 241.620759, 252.508513, 252.767900, 244.445283, 260.870105, 217.352348, 247.017507, 246.182278, 267.772651, 232.277540, 258.306589, 219.820165, 236.805076, 267.239982, 284.795470, 234.913508, 289.709239, 241.312315, 247.249491, 286.702147, 264.252852, 245.151647, 264.582984, 240.842689, 245.837476, 223.602907, 234.242638, 232.823977, 227.436840, 240.751687, 203.533069, 230.484479, 225.727773, 248.459822, 216.298768, 240.263617, 202.894307, 216.680343, 248.073140, 263.122266, 218.066855, 241.312315, 255.363057, 230.209787, 271.091482, 239.220241, 225.387834, 247.486715, 226.052431, 224.119935, 220.038678, 240.308477, 233.084574, 231.369510, 247.183922, 204.581690, 233.358639, 229.873732, 249.662698, 220.296241, 246.844456, 211.468055, 223.816297, 254.230916, 271.239770, 219.648997, 247.249491, 230.209787, 264.014907, 271.938970, 246.664305, 227.889045, 249.908085, 232.035369, 229.010298, 264.704959, 285.416390, 278.315614, 270.721492, 289.044453, 243.603643, 274.935489, 266.295057, 295.935799, 253.817633, 290.492875, 246.048907, 263.188711, 296.789984, 320.209823, 257.969951, 286.702147, 271.091482, 271.938970, 352.825726, 286.200221, 267.716897, 297.182554, 269.776351, 266.721561, 240.341716, 254.792243, 250.872786, 242.475130, 260.715749, 217.785986, 249.237737, 242.954024, 267.460908, 237.638235, 267.212598, 225.135194, 236.311810, 267.158551, 283.933299, 231.243624, 264.252852, 239.220241, 246.664305, 286.200221, 294.042749, 246.504021, 269.570596, 243.980697, 242.690997, 223.672378, 240.176223, 235.227909, 226.471530, 244.210258, 204.205559, 235.229584, 225.510760, 255.668926, 220.785141, 243.218596, 210.076305, 222.950152, 252.226496, 264.416752, 224.657569, 245.151647, 225.387834, 227.889045, 267.716897, 246.504021, 259.897656, 251.730847, 229.335712, 229.759185, 244.323303, 259.048267, 255.104263, 248.130112, 267.159502, 223.747953, 253.029955, 249.370268, 275.902272, 239.098500, 265.681905, 226.806762, 244.636356, 271.170860, 292.035194, 238.158714, 264.582984, 247.486715, 249.908085, 297.182554, 269.570596, 251.730847, 303.872223, 251.585636, 247.878402, 223.216217, 235.566855, 238.931093, 225.826557, 242.992822, 203.586842, 228.601700, 227.813265, 248.495606, 220.583355, 244.615890, 212.014431, 219.121372, 248.325354, 268.764031, 217.174368, 240.842689, 226.052431, 232.035369, 269.776351, 243.980697, 229.335712, 251.585636, 257.544914, 228.810942, 226.074990, 236.277617, 235.402356, 228.266719, 244.070245, 200.165867, 230.512885, 232.141964, 246.827505, 218.962732, 242.543363, 205.123267, 219.911078, 253.694013, 265.345807, 215.933866, 245.837476, 224.119935, 229.010298, 266.721561, 242.690997, 229.759185, 247.878402, 228.810942, 253.353769)) val fl1 = DenseVector(-892.842851, -934.071560, -932.936015, -888.124343, -961.050207, -791.191087, -923.711397, -904.289301, -988.384984, -883.909133, -959.465030, -798.551172, -871.622303, -997.463289, -1043.912620, -863.013719, -976.975712, -897.033693, -898.694786, -1069.245497, -963.491924, -901.263474, -983.768031, -899.865392, -902.283567) @@ -222,6 +219,26 @@ class QuadraticMinimizerTest extends OptimizeTestBase with Matchers { assert(norm(qpMlL1Result.x - octaveL1, 2) < 1e-4) } + test("Quadratic Minimization with L1 projection compared to Octave") { + val Hl1 = new DenseMatrix(25, 25, Array(253.535098, 236.477785, 234.421906, 223.374867, 241.007512, 204.695511, 226.465507, 223.351032, 249.179386, 221.411909, 238.679352, 203.579010, 217.564498, 243.852681, 266.607649, 213.994496, 241.620759, 223.602907, 220.038678, 264.704959, 240.341716, 223.672378, 244.323303, 223.216217, 226.074990, 236.477785, 278.862035, 245.756639, 237.489890, 252.783139, 214.077652, 241.816953, 238.790633, 260.536460, 228.121417, 255.103936, 216.608405, 237.150426, 258.933231, 281.958112, 228.971242, 252.508513, 234.242638, 240.308477, 285.416390, 254.792243, 240.176223, 259.048267, 235.566855, 236.277617, 234.421906, 245.756639, 269.162882, 231.416867, 251.321527, 208.134322, 236.567647, 236.558029, 255.805108, 226.535825, 251.514713, 212.770208, 228.565362, 261.748652, 273.946966, 227.411615, 252.767900, 232.823977, 233.084574, 278.315614, 250.872786, 235.227909, 255.104263, 238.931093, 235.402356, 223.374867, 237.489890, 231.416867, 254.771963, 241.703229, 209.028084, 231.517998, 228.768510, 250.805315, 216.548935, 245.473869, 207.687875, 222.036114, 250.906955, 263.018181, 216.128966, 244.445283, 227.436840, 231.369510, 270.721492, 242.475130, 226.471530, 248.130112, 225.826557, 228.266719, 241.007512, 252.783139, 251.321527, 241.703229, 285.702320, 219.051868, 249.442308, 240.400187, 264.970407, 232.503138, 258.819837, 220.160683, 235.621356, 267.743972, 285.795029, 229.667231, 260.870105, 240.751687, 247.183922, 289.044453, 260.715749, 244.210258, 267.159502, 242.992822, 244.070245, 204.695511, 214.077652, 208.134322, 209.028084, 219.051868, 210.164224, 208.151908, 201.539036, 226.373834, 192.056565, 219.950686, 191.459568, 195.982460, 226.739575, 240.677519, 196.116652, 217.352348, 203.533069, 204.581690, 243.603643, 217.785986, 204.205559, 223.747953, 203.586842, 200.165867, 226.465507, 241.816953, 236.567647, 231.517998, 249.442308, 208.151908, 264.007925, 227.080718, 253.174653, 220.322823, 248.619983, 210.100242, 223.279198, 254.807401, 269.896959, 222.927882, 247.017507, 230.484479, 233.358639, 274.935489, 249.237737, 235.229584, 253.029955, 228.601700, 230.512885, 223.351032, 238.790633, 236.558029, 228.768510, 240.400187, 201.539036, 227.080718, 258.773479, 249.471480, 215.664539, 243.078577, 202.337063, 221.020998, 249.979759, 263.356244, 213.470569, 246.182278, 225.727773, 229.873732, 266.295057, 242.954024, 225.510760, 249.370268, 227.813265, 232.141964, 249.179386, 260.536460, 255.805108, 250.805315, 264.970407, 226.373834, 253.174653, 249.471480, 302.360150, 237.902729, 265.769812, 224.947876, 243.088105, 273.690377, 291.076027, 241.089661, 267.772651, 248.459822, 249.662698, 295.935799, 267.460908, 255.668926, 275.902272, 248.495606, 246.827505, 221.411909, 228.121417, 226.535825, 216.548935, 232.503138, 192.056565, 220.322823, 215.664539, 237.902729, 245.154567, 234.956316, 199.557862, 214.774631, 240.339217, 255.161923, 209.328714, 232.277540, 216.298768, 220.296241, 253.817633, 237.638235, 220.785141, 239.098500, 220.583355, 218.962732, 238.679352, 255.103936, 251.514713, 245.473869, 258.819837, 219.950686, 248.619983, 243.078577, 265.769812, 234.956316, 288.133073, 225.087852, 239.810430, 268.406605, 283.289840, 233.858455, 258.306589, 240.263617, 246.844456, 290.492875, 267.212598, 243.218596, 265.681905, 244.615890, 242.543363, 203.579010, 216.608405, 212.770208, 207.687875, 220.160683, 191.459568, 210.100242, 202.337063, 224.947876, 199.557862, 225.087852, 217.501685, 197.897572, 229.825316, 242.175607, 201.123644, 219.820165, 202.894307, 211.468055, 246.048907, 225.135194, 210.076305, 226.806762, 212.014431, 205.123267, 217.564498, 237.150426, 228.565362, 222.036114, 235.621356, 195.982460, 223.279198, 221.020998, 243.088105, 214.774631, 239.810430, 197.897572, 244.439113, 241.621129, 260.400953, 216.482178, 236.805076, 216.680343, 223.816297, 263.188711, 236.311810, 222.950152, 244.636356, 219.121372, 219.911078, 243.852681, 258.933231, 261.748652, 250.906955, 267.743972, 226.739575, 254.807401, 249.979759, 273.690377, 240.339217, 268.406605, 229.825316, 241.621129, 302.928261, 288.344398, 238.549018, 267.239982, 248.073140, 254.230916, 296.789984, 267.158551, 252.226496, 271.170860, 248.325354, 253.694013, 266.607649, 281.958112, 273.946966, 263.018181, 285.795029, 240.677519, 269.896959, 263.356244, 291.076027, 255.161923, 283.289840, 242.175607, 260.400953, 288.344398, 343.457361, 257.368309, 284.795470, 263.122266, 271.239770, 320.209823, 283.933299, 264.416752, 292.035194, 268.764031, 265.345807, 213.994496, 228.971242, 227.411615, 216.128966, 229.667231, 196.116652, 222.927882, 213.470569, 241.089661, 209.328714, 233.858455, 201.123644, 216.482178, 238.549018, 257.368309, 239.295031, 234.913508, 218.066855, 219.648997, 257.969951, 231.243624, 224.657569, 238.158714, 217.174368, 215.933866, 241.620759, 252.508513, 252.767900, 244.445283, 260.870105, 217.352348, 247.017507, 246.182278, 267.772651, 232.277540, 258.306589, 219.820165, 236.805076, 267.239982, 284.795470, 234.913508, 289.709239, 241.312315, 247.249491, 286.702147, 264.252852, 245.151647, 264.582984, 240.842689, 245.837476, 223.602907, 234.242638, 232.823977, 227.436840, 240.751687, 203.533069, 230.484479, 225.727773, 248.459822, 216.298768, 240.263617, 202.894307, 216.680343, 248.073140, 263.122266, 218.066855, 241.312315, 255.363057, 230.209787, 271.091482, 239.220241, 225.387834, 247.486715, 226.052431, 224.119935, 220.038678, 240.308477, 233.084574, 231.369510, 247.183922, 204.581690, 233.358639, 229.873732, 249.662698, 220.296241, 246.844456, 211.468055, 223.816297, 254.230916, 271.239770, 219.648997, 247.249491, 230.209787, 264.014907, 271.938970, 246.664305, 227.889045, 249.908085, 232.035369, 229.010298, 264.704959, 285.416390, 278.315614, 270.721492, 289.044453, 243.603643, 274.935489, 266.295057, 295.935799, 253.817633, 290.492875, 246.048907, 263.188711, 296.789984, 320.209823, 257.969951, 286.702147, 271.091482, 271.938970, 352.825726, 286.200221, 267.716897, 297.182554, 269.776351, 266.721561, 240.341716, 254.792243, 250.872786, 242.475130, 260.715749, 217.785986, 249.237737, 242.954024, 267.460908, 237.638235, 267.212598, 225.135194, 236.311810, 267.158551, 283.933299, 231.243624, 264.252852, 239.220241, 246.664305, 286.200221, 294.042749, 246.504021, 269.570596, 243.980697, 242.690997, 223.672378, 240.176223, 235.227909, 226.471530, 244.210258, 204.205559, 235.229584, 225.510760, 255.668926, 220.785141, 243.218596, 210.076305, 222.950152, 252.226496, 264.416752, 224.657569, 245.151647, 225.387834, 227.889045, 267.716897, 246.504021, 259.897656, 251.730847, 229.335712, 229.759185, 244.323303, 259.048267, 255.104263, 248.130112, 267.159502, 223.747953, 253.029955, 249.370268, 275.902272, 239.098500, 265.681905, 226.806762, 244.636356, 271.170860, 292.035194, 238.158714, 264.582984, 247.486715, 249.908085, 297.182554, 269.570596, 251.730847, 303.872223, 251.585636, 247.878402, 223.216217, 235.566855, 238.931093, 225.826557, 242.992822, 203.586842, 228.601700, 227.813265, 248.495606, 220.583355, 244.615890, 212.014431, 219.121372, 248.325354, 268.764031, 217.174368, 240.842689, 226.052431, 232.035369, 269.776351, 243.980697, 229.335712, 251.585636, 257.544914, 228.810942, 226.074990, 236.277617, 235.402356, 228.266719, 244.070245, 200.165867, 230.512885, 232.141964, 246.827505, 218.962732, 242.543363, 205.123267, 219.911078, 253.694013, 265.345807, 215.933866, 245.837476, 224.119935, 229.010298, 266.721561, 242.690997, 229.759185, 247.878402, 228.810942, 253.353769)) + val fl1 = DenseVector(-892.842851, -934.071560, -932.936015, -888.124343, -961.050207, -791.191087, -923.711397, -904.289301, -988.384984, -883.909133, -959.465030, -798.551172, -871.622303, -997.463289, -1043.912620, -863.013719, -976.975712, -897.033693, -898.694786, -1069.245497, -963.491924, -901.263474, -983.768031, -899.865392, -902.283567) + + val octaveL1 = DenseVector(0.18611, 0.00000, 0.06317, -0.10417, 0.11262, + -0.20495, 0.52668, 0.32790, 0.19421, 0.72180, + 0.06309, -0.41326, -0.00000, 0.52078, -0.00000, + 0.18040, 0.62915, 0.16329, -0.06424, 0.37539, + 0.01659, 0.00000, 0.11215, 0.24778, 0.04082) + + val s = octaveL1.foldLeft(0.0) { case (agg, entry) => agg + abs(entry)} + val projectL1 = ProjectL1(s) + + val qpResultProject = new QuadraticMinimizer(25, proximal = projectL1).minimizeAndReturnState(Hl1, fl1) + val qpResultProximal = QuadraticMinimizer(25, SPARSE, 2.0).minimizeAndReturnState(Hl1, fl1) + + assert(norm(qpResultProximal.x - octaveL1, 2) < 1e-4) + assert(norm(qpResultProject.x - octaveL1, 2) < 1e-4) + } + test("Quadratic Minimization with Positivity High Iterations compared to NNLS") { //Movielens Positive/Bounds high iterations val P1 = new DenseMatrix[Double](20, 20, Array(539.101887, 494.598042, 505.700671, 505.846716, 504.700928, 516.629473, 507.958246, 514.096818, 514.801371, 505.735357, 507.322795, 522.547578, 498.320793, 502.829895, 505.847128, 488.934012, 516.116942, 501.906569, 505.627629, 496.409513, 494.598042, 565.723334, 513.239749, 519.155649, 514.070934, 524.154020, 521.694985, 523.512877, 521.122745, 513.862711, 518.653059, 530.426712, 511.054588, 510.096410, 521.373582, 503.132142, 531.736861, 514.161101, 515.005997, 500.799198, 505.700671, 513.239749, 602.920633, 524.780488, 547.978722, 558.807137, 526.999189, 553.273432, 552.657103, 547.690555, 537.912646, 563.616990, 527.634170, 541.947698, 524.060188, 507.650395, 534.403391, 534.406246, 546.625588, 535.221534, 505.846716, 519.155649, 524.780488, 585.686194, 522.548624, 537.124362, 534.911663, 530.505003, 533.364761, 525.544862, 530.149606, 543.063850, 518.884670, 517.707324, 531.252004, 511.635097, 541.217141, 522.706817, 526.063019, 513.574796, 504.700928, 514.070934, 547.978722, 522.548624, 602.711620, 557.824564, 523.454584, 556.453086, 551.975932, 548.308951, 537.908414, 562.811202, 530.685949, 544.687533, 523.961961, 507.966023, 534.868428, 535.823631, 546.491009, 534.209994, 516.629473, 524.154020, 558.807137, 537.124362, 557.824564, 624.863357, 539.305550, 566.795282, 566.357133, 560.044835, 550.021868, 578.832952, 538.509515, 554.979692, 535.459614, 518.812719, 546.530562, 546.168894, 558.181118, 548.792405, 507.958246, 521.694985, 526.999189, 534.911663, 523.454584, 539.305550, 591.413305, 532.696713, 536.852122, 527.232598, 531.751544, 545.578374, 520.832790, 521.212336, 533.128446, 513.908219, 544.084087, 525.037440, 527.089172, 515.549361, 514.096818, 523.512877, 553.273432, 530.505003, 556.453086, 566.795282, 532.696713, 621.615462, 562.996181, 554.493454, 545.658444, 574.476466, 538.757428, 555.556328, 533.365769, 517.189079, 543.997925, 544.795736, 554.632016, 544.157448, 514.801371, 521.122745, 552.657103, 533.364761, 551.975932, 566.357133, 536.852122, 562.996181, 621.891951, 552.016433, 544.239294, 573.926455, 532.639189, 553.095240, 530.294369, 516.269094, 544.643533, 537.922422, 549.757806, 545.057129, 505.735357, 513.862711, 547.690555, 525.544862, 548.308951, 560.044835, 527.232598, 554.493454, 552.016433, 603.860301, 539.638699, 564.591708, 529.636837, 543.558799, 524.759472, 506.968931, 534.506284, 537.156481, 547.911779, 536.271289, 507.322795, 518.653059, 537.912646, 530.149606, 537.908414, 550.021868, 531.751544, 545.658444, 544.239294, 539.638699, 589.023204, 554.623574, 528.460278, 532.243589, 530.559000, 512.720156, 539.864087, 531.591296, 538.972277, 527.533747, 522.547578, 530.426712, 563.616990, 543.063850, 562.811202, 578.832952, 545.578374, 574.476466, 573.926455, 564.591708, 554.623574, 638.367746, 544.269457, 561.616144, 541.245807, 524.820018, 552.498719, 552.102404, 563.397777, 554.806495, 498.320793, 511.054588, 527.634170, 518.884670, 530.685949, 538.509515, 520.832790, 538.757428, 532.639189, 529.636837, 528.460278, 544.269457, 576.086269, 524.503702, 521.427126, 503.797286, 531.685079, 525.230922, 529.783815, 520.222911, 502.829895, 510.096410, 541.947698, 517.707324, 544.687533, 554.979692, 521.212336, 555.556328, 553.095240, 543.558799, 532.243589, 561.616144, 524.503702, 597.742352, 520.329969, 504.306186, 530.688065, 531.679956, 541.233986, 531.982099, 505.847128, 521.373582, 524.060188, 531.252004, 523.961961, 535.459614, 533.128446, 533.365769, 530.294369, 524.759472, 530.559000, 541.245807, 521.427126, 520.329969, 586.613767, 513.362642, 542.438988, 524.690396, 526.354706, 511.243067, 488.934012, 503.132142, 507.650395, 511.635097, 507.966023, 518.812719, 513.908219, 517.189079, 516.269094, 506.968931, 512.720156, 524.820018, 503.797286, 504.306186, 513.362642, 550.229150, 523.966242, 506.244433, 507.769334, 495.459277, 516.116942, 531.736861, 534.403391, 541.217141, 534.868428, 546.530562, 544.084087, 543.997925, 544.643533, 534.506284, 539.864087, 552.498719, 531.685079, 530.688065, 542.438988, 523.966242, 606.957114, 533.719037, 534.295092, 522.409881, 501.906569, 514.161101, 534.406246, 522.706817, 535.823631, 546.168894, 525.037440, 544.795736, 537.922422, 537.156481, 531.591296, 552.102404, 525.230922, 531.679956, 524.690396, 506.244433, 533.719037, 583.853405, 536.033903, 522.009144, 505.627629, 515.005997, 546.625588, 526.063019, 546.491009, 558.181118, 527.089172, 554.632016, 549.757806, 547.911779, 538.972277, 563.397777, 529.783815, 541.233986, 526.354706, 507.769334, 534.295092, 536.033903, 599.942302, 535.001472, 496.409513, 500.799198, 535.221534, 513.574796, 534.209994, 548.792405, 515.549361, 544.157448, 545.057129, 536.271289, 527.533747, 554.806495, 520.222911, 531.982099, 511.243067, 495.459277, 522.409881, 522.009144, 535.001472, 593.140288))