diff --git a/LICENSE b/LICENSE index bb5e4b237..2d2f1a9dc 100644 --- a/LICENSE +++ b/LICENSE @@ -174,3 +174,32 @@ of your accepting any such warranty or additional liability. END OF TERMS AND CONDITIONS + +======================================================================== +BSD-style licenses +======================================================================== +(BSD License) Proximal algorithms (https://github.com/cvxgrp/proximal) + Copyright (c) 2013, Neal Parikh and Stephen Boyd (Stanford University) + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of Stanford University nor the + names of its contributors may be used to endorse or promote products + derived from this software without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL NEAL PARIKH OR STEPHEN BOYD BE LIABLE FOR ANY + DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. \ No newline at end of file diff --git a/NOTICE b/NOTICE index b7127a937..bd6ae4083 100644 --- a/NOTICE +++ b/NOTICE @@ -1 +1,13 @@ Breeze is distributed under an Apache License V2.0 (See LICENSE) + +=============================================================================== + +Proximal algorithms outlined in Proximal.scala (package breeze.optimize.quadratic) +are based on https://github.com/cvxgrp/proximal (see LICENSE for details) and distributed with +Copyright (c) 2014 by Debasish Das (Verizon), all rights reserved. + +=============================================================================== + +QuadraticMinimizer class in package breeze.optimize.proximal is distributed with Copyright (c) +2014, Debasish Das (Verizon), all rights reserved. + diff --git a/build.sbt b/build.sbt index eb9ee761c..7f10ab090 100644 --- a/build.sbt +++ b/build.sbt @@ -51,4 +51,3 @@ pomExtra := ( http://cs.berkeley.edu/~dlwh/ ) - diff --git a/math/build.sbt b/math/build.sbt index 7d951fd1f..e83d447e5 100644 --- a/math/build.sbt +++ b/math/build.sbt @@ -100,5 +100,3 @@ testOptions in Test += Tests.Argument("-oDF") //fork in Test := true javaOptions := Seq("-Xmx4g") - -jacoco.settings diff --git a/math/src/main/scala/breeze/linalg/LinearAlgebraException.scala b/math/src/main/scala/breeze/linalg/LinearAlgebraException.scala index d05909f61..905cb09f1 100644 --- a/math/src/main/scala/breeze/linalg/LinearAlgebraException.scala +++ b/math/src/main/scala/breeze/linalg/LinearAlgebraException.scala @@ -31,3 +31,5 @@ class MatrixEmptyException extends IllegalArgumentException("Matrix is empty") w * @author dramage, dlwh */ class MatrixSingularException(msg : String="") extends RuntimeException(msg) with LinearAlgebraException + +class LapackException(msg: String="") extends RuntimeException(msg) with LinearAlgebraException \ No newline at end of file diff --git a/math/src/main/scala/breeze/optimize/linear/NNLS.scala b/math/src/main/scala/breeze/optimize/linear/NNLS.scala new file mode 100644 index 000000000..f9ae7eccb --- /dev/null +++ b/math/src/main/scala/breeze/optimize/linear/NNLS.scala @@ -0,0 +1,183 @@ +package breeze.optimize.linear + +import breeze.linalg.{DenseMatrix, DenseVector, axpy} +import breeze.optimize.proximal.QpGenerator +import breeze.stats.distributions.Rand +import breeze.util.Implicits._ +import spire.syntax.cfor._ + +/** + * NNLS solves nonnegative least squares problems using a modified + * projected gradient method. + * @param maxIters user defined maximum iterations + * @author debasish83, coderxiang + */ + +class NNLS(val maxIters: Int = -1) { + type BDM = DenseMatrix[Double] + type BDV = DenseVector[Double] + + case class State private[NNLS](x: BDV, grad: BDV, dir: BDV, lastDir: BDV, res: BDV, lastNorm: Double, lastWall: Int, iter: Int, converged: Boolean) + + // find the optimal unconstrained step + private def steplen(ata: BDM, dir: BDV, res: BDV, + tmp: BDV): Double = { + val top = dir.dot(res) + tmp := ata * dir + // Push the denominator upward very slightly to avoid infinities and silliness + top / (tmp.dot(dir) + 1e-20) + } + + // stopping condition + private def stop(step: Double, ndir: Double, nx: Double): Boolean = { + ((step.isNaN) // NaN + || (step < 1e-7) // too small or negative + || (step > 1e40) // too small; almost certainly numerical problems + || (ndir < 1e-12 * nx) // gradient relatively too small + || (ndir < 1e-32) // gradient absolutely too small; numerical issues may lurk + ) + } + + /** + * Solve a least squares problem, possibly with nonnegativity constraints, by a modified + * projected gradient method. That is, find x minimising ||Ax - b||_2 given A^T A and A^T b. + * + * We solve the problem + * min_x 1/2 x' ata x' - x'atb + * subject to x >= 0 + * + * The method used is similar to one described by Polyak (B. T. Polyak, The conjugate gradient + * method in extremal problems, Zh. Vychisl. Mat. Mat. Fiz. 9(4)(1969), pp. 94-112) for bound- + * constrained nonlinear programming. Polyak unconditionally uses a conjugate gradient + * direction, however, while this method only uses a conjugate gradient direction if the last + * iteration did not cause a previously-inactive constraint to become active. + */ + private def initialState(ata: DenseMatrix[Double], atb: DenseVector[Double], n: Int): State = { + require(ata.cols == ata.rows, s"NNLS:iterations gram matrix must be symmetric") + require(ata.rows == atb.length, s"NNLS:iterations gram matrix rows must be same as length of linear term") + + val grad = DenseVector.zeros[Double](n) + val x = DenseVector.zeros[Double](n) + val dir = DenseVector.zeros[Double](n) + val lastDir = DenseVector.zeros[Double](n) + val res = DenseVector.zeros[Double](n) + val lastNorm = 0.0 + val lastWall = 0 + State(x, grad, dir, lastDir, res, lastNorm, lastWall, 0, false) + } + + def iterations(ata: DenseMatrix[Double], + atb: DenseVector[Double]): Iterator[State] = + Iterator.iterate(initialState(ata, atb, atb.length)) { state => + import state._ + val n = atb.length + val tmp = DenseVector.zeros[Double](atb.length) + + val iterMax = if (maxIters < 0) Math.max(400, 20 * n) else maxIters + + // find the residual + res := ata * x + res -= atb + grad := res + + // project the gradient + cforRange(0 until n) { i => + if (grad(i) > 0.0 && x(i) == 0.0) { + grad(i) = 0.0 + } + } + val ngrad = grad.dot(grad) + dir := grad + + // use a CG direction under certain conditions + var step = steplen(ata, grad, res, tmp) + var ndir = 0.0 + val nx = x.dot(x) + + if (iter > lastWall + 1) { + val alpha = ngrad / lastNorm + axpy(alpha, lastDir, dir) + val dstep = steplen(ata, dir, res, tmp) + ndir = dir.dot(dir) + if (stop(dstep, ndir, nx)) { + // reject the CG step if it could lead to premature termination + dir := grad + ndir = dir.dot(dir) + } else { + step = dstep + } + } else { + ndir = dir.dot(dir) + } + + // terminate? + if (stop(step, ndir, nx) || iter > iterMax) + State(x, grad, dir, lastDir, res, lastNorm, lastWall, iter + 1, true) + else { + // don't run through the walls + cforRange(0 until n){ i => + if (step * dir(i) > x(i)) { + step = x(i) / dir(i) + } + } + + var nextWall = lastWall + + // take the step + cforRange(0 until n) { i => + if (step * dir(i) > x(i) * (1 - 1e-14)) { + x(i) = 0 + nextWall = iter + } else { + x(i) -= step * dir(i) + } + } + lastDir := dir + val nextNorm = ngrad + State(x, grad, dir, lastDir, res, nextNorm, nextWall, iter + 1, false) + } + }.takeUpToWhere(_.converged) + + def minimizeAndReturnState(ata: DenseMatrix[Double], + atb: DenseVector[Double]): State = { + iterations(ata, atb).last + } + + def minimize(ata: DenseMatrix[Double], atb: DenseVector[Double]): DenseVector[Double] = { + minimizeAndReturnState(ata, atb).x + } +} + +object NNLS { + /** Compute the objective value */ + def computeObjectiveValue(ata: DenseMatrix[Double], atb: DenseVector[Double], x: DenseVector[Double]): Double = { + val res = (x.t*ata*x)*0.5 - atb.dot(x) + res + } + + def apply(iters: Int) = new NNLS(iters) + + def main(args: Array[String]) { + if (args.length < 2) { + println("Usage: NNLS n s") + println("Test NNLS with quadratic function of dimension n for s consecutive solves") + sys.exit(1) + } + + val problemSize = args(0).toInt + val numSolves = args(1).toInt + val nnls = new NNLS() + + var i = 0 + var nnlsTime = 0L + while(i < numSolves) { + val ata = QpGenerator.getGram(problemSize) + val atb = DenseVector.rand[Double](problemSize, Rand.gaussian(0,1)) + val startTime = System.nanoTime() + nnls.minimize(ata, atb) + nnlsTime = nnlsTime + (System.nanoTime() - startTime) + i = i + 1 + } + println(s"NNLS problemSize $problemSize solves $numSolves ${nnlsTime/1e6} ms") + } +} diff --git a/math/src/main/scala/breeze/optimize/linear/PowerMethod.scala b/math/src/main/scala/breeze/optimize/linear/PowerMethod.scala new file mode 100644 index 000000000..1d8fd0640 --- /dev/null +++ b/math/src/main/scala/breeze/optimize/linear/PowerMethod.scala @@ -0,0 +1,83 @@ +package breeze.optimize.linear + +import breeze.linalg.operators.OpMulMatrix +import breeze.math.MutableInnerProductModule +import breeze.numerics.abs +import breeze.optimize.proximal.QuadraticMinimizer +import breeze.util.SerializableLogging +import breeze.linalg.{DenseMatrix, DenseVector, norm} +import breeze.util.Implicits._ +/** + * Power Method to compute maximum eigen value and companion object to compute minimum eigen value through inverse power + * iterations + * @author debasish83 + */ +class PowerMethod[T, M](maxIters: Int = 10,tolerance: Double = 1E-5) + (implicit space: MutableInnerProductModule[T, Double], mult: OpMulMatrix.Impl2[M, T, T]) extends SerializableLogging { + + import space._ + + case class State private[PowerMethod] (eigenValue: Double, eigenVector: T, iter: Int, converged: Boolean) + + //memory allocation for the eigen vector result + def normalize(y: T) : T = { + val normInit = norm(y) + val init = copy(y) + init *= 1.0/normInit + } + + def initialState(y: T, A: M): State = { + val ynorm = normalize(y) + val ay = mult(A, ynorm) + val lambda = nextEigen(ynorm, ay) + State(lambda, ynorm, 0, false) + } + + //in-place modification of eigen vector + def nextEigen(eigenVector: T, ay: T) = { + val lambda = eigenVector dot ay + eigenVector := ay + val norm1 = norm(ay) + eigenVector *= 1.0/norm1 + if (lambda < 0.0) eigenVector *= -1.0 + lambda + } + + def iterations(y: T, + A: M): Iterator[State] = Iterator.iterate(initialState(y, A)) { state => + import state._ + val ay = mult(A, eigenVector) + val lambda = nextEigen(eigenVector, ay) + val val_dif = abs(lambda - eigenValue) + if (val_dif <= tolerance || iter > maxIters) State(lambda, eigenVector, iter + 1, true) + else State(lambda, eigenVector, iter + 1, false) + }.takeUpToWhere(_.converged) + + def iterateAndReturnState(y: T, A: M): State = { + iterations(y, A).last + } + + def eigen(y: T, A: M): Double = { + iterateAndReturnState(y, A).eigenValue + } +} + +object PowerMethod { + def inverse(maxIters: Int = 10, tolerance: Double = 1E-5) : PowerMethod[DenseVector[Double], DenseMatrix[Double]] = new PowerMethod[DenseVector[Double], DenseMatrix[Double]](maxIters, tolerance) { + override def initialState(y: DenseVector[Double], A: DenseMatrix[Double]): State = { + val ynorm = normalize(y) + val ay = QuadraticMinimizer.solveTriangular(A, ynorm) + val lambda = nextEigen(ynorm, ay) + State(lambda, ynorm, 0, false) + } + override def iterations(y: DenseVector[Double], + A: DenseMatrix[Double]): Iterator[State] = Iterator.iterate(initialState(y, A)) { state => + import state._ + val ay = QuadraticMinimizer.solveTriangular(A, eigenVector) + val lambda = nextEigen(eigenVector, ay) + val val_dif = abs(lambda - eigenValue) + if (val_dif <= tolerance || iter > maxIters) State(lambda, eigenVector, iter + 1, true) + else State(lambda, eigenVector, iter + 1, false) + }.takeUpToWhere(_.converged) + } +} \ No newline at end of file diff --git a/math/src/main/scala/breeze/optimize/proximal/Constraint.scala b/math/src/main/scala/breeze/optimize/proximal/Constraint.scala new file mode 100644 index 000000000..c3f647551 --- /dev/null +++ b/math/src/main/scala/breeze/optimize/proximal/Constraint.scala @@ -0,0 +1,10 @@ +package breeze.optimize.proximal + +/** + * Supported constraints by QuadraticMinimizer object + * @author debasish83 + */ +object Constraint extends Enumeration { + type Constraint = Value + val SMOOTH, POSITIVE, BOX, SPARSE, EQUALITY = Value +} \ No newline at end of file diff --git a/math/src/main/scala/breeze/optimize/proximal/Proximal.scala b/math/src/main/scala/breeze/optimize/proximal/Proximal.scala new file mode 100644 index 000000000..ecb7c4dd5 --- /dev/null +++ b/math/src/main/scala/breeze/optimize/proximal/Proximal.scala @@ -0,0 +1,197 @@ +/* + * Library of Proximal Algorithms adapted from https://github.com/cvxgrp/proximal + * In-place modifications which later should be BLAS-ed when applicable for more efficiency + * @author debasish83 + */ + +package breeze.optimize.proximal + +import breeze.numerics.pow + +import scala.math.max +import scala.math.min +import scala.math.sqrt +import scala.math.abs +import scala.Double.NegativeInfinity +import scala.Double.PositiveInfinity +import breeze.linalg._ +import spire.syntax.cfor._ +import breeze.linalg.norm + +trait Proximal { + def prox(x: DenseVector[Double], rho: Double) +} + +case class ProjectBox(l: DenseVector[Double], u: DenseVector[Double]) extends Proximal { + def prox(x: DenseVector[Double], rho: Double = 0.0) = { + cforRange(0 until x.length) { i => x.update(i, max(l(i), min(x(i), u(i))))} + } +} + +case class ProjectPos() extends Proximal { + def prox(x: DenseVector[Double], rho: Double = 0.0) = { + cforRange(0 until x.length) { i => x.update(i, max(0, x(i)))} + } +} + +case class ProjectSoc() extends Proximal { + def prox(x: DenseVector[Double], rho: Double = 0.0) = { + var nx: Double = 0.0 + val n = x.length + cforRange(1 until n) { i => + nx += x(i) * x(i) + } + nx = sqrt(nx) + + if (nx > x(0)) { + if (nx <= -x(0)) { + cforRange(0 until n ) { + i => x(i) = 0 + } + } else { + val alpha = 0.5 * (1 + x(0) / nx) + x.update(0, alpha * nx) + cforRange(1 until n) { + i => x.update(i, alpha * x(i)) + } + } + } + } +} + +//Projection onto Affine set +//Let C = { x \in R^{n} | Ax = b } where A \in R^{m x n} +//If A is full rank matrix then the projection is given by v - A'(Av - b) where A' is the cached Moore-Penrose pseudo-inverse of A +case class ProjectEquality(Aeq: DenseMatrix[Double], beq: DenseVector[Double]) extends Proximal { + val invAeq = pinv(Aeq) + def prox(x: DenseVector[Double], rho: Double = 0.0) = { + val Av = Aeq*x + Av -= beq + x += invAeq*Av + } +} + +//Projection onto hyper-plane is a special case of projection onto affine set and is given by +//x + ((b - a'x)/||a||_2^2)a +case class ProjectHyperPlane(a: DenseVector[Double], b: Double) extends Proximal { + val at = a.t + + def prox(x: DenseVector[Double], rho: Double = 0.0) = { + val atx = at * x + val anorm = norm(a, 2) + val scale = (b - atx) / (anorm * anorm) + val ascaled = a * scale + x += ascaled + } +} + +case class ProximalL1() extends Proximal { + var lambda = 1.0 + + def setLambda(lambda: Double) = { + this.lambda = lambda + this + } + + def prox(x: DenseVector[Double], rho: Double) = { + cforRange(0 until x.length) { i => + x.update(i, max(0, x(i) - lambda / rho) - max(0, -x(i) - lambda / rho)) + } + } +} + +case class ProximalL2() extends Proximal { + def prox(x: DenseVector[Double], rho: Double) = { + val xnorm = norm(x) + cforRange(0 until x.length) { i => + if (xnorm >= 1 / rho) x.update(i, x(i) * (1 - 1 / (rho * xnorm))) + else x.update(i, 0) + } + } +} + +// f = (1/2)||.||_2^2 +case class ProximalSumSquare() extends Proximal { + def prox(x: DenseVector[Double], rho: Double) = { + cforRange(0 until x.length) { i => + x.update(i, x(i) * (rho / (1 + rho))) + } + } +} + +// f = -sum(log(x)) +case class ProximalLogBarrier() extends Proximal { + def prox(x: DenseVector[Double], rho: Double) = { + cforRange(0 until x.length) { i => + x.update(i, 0.5 * (x(i) + sqrt(x(i) * x(i) + 4 / rho))) + } + } +} + +// f = huber = x^2 if |x|<=1, 2|x| - 1 otherwise +case class ProximalHuber() extends Proximal { + def proxScalar(v: Double, rho: Double, oracle: Double => Double, l: Double, u: Double, x0: Double): Double = { + val MAX_ITER = 1000 + val tol = 1e-8 + + var g: Double = 0.0 + var x = max(l, min(x0, u)) + + var lIter = l + var uIter = u + var iter = 0 + + while (iter < MAX_ITER && u - l > tol) { + g = -1 / x + rho * (x - v) + + if (g > 0) { + lIter = max(lIter, x - g / rho) + uIter = x + } else if (g < 0) { + lIter = x + uIter = min(uIter, x - g / rho) + } + x = (lIter + uIter) / 2 + iter = iter + 1 + } + x + } + + def proxSeparable(x: DenseVector[Double], rho: Double, oracle: Double => Double, l: Double, u: Double) = { + x.map(proxScalar(_, rho, oracle, l, u, 0)) + cforRange(0 until x.length) { i => + x.update(i, proxScalar(x(i), rho, oracle, l, u, 0)) + } + } + + def subgradHuber(x: Double): Double = { + if (abs(x) <= 1) { + 2 * x + } else { + val projx = if (x > 0) x else -x + 2 * projx + } + } + + def prox(x: DenseVector[Double], rho: Double) = { + proxSeparable(x, rho, subgradHuber, NegativeInfinity, PositiveInfinity) + } +} + +// f = c'*x +case class ProximalLinear(c: DenseVector[Double]) extends Proximal { + def prox(x: DenseVector[Double], rho: Double) = { + cforRange(0 until x.length) { i => + x.update(i, x(i) - c(i) / rho) + } + } +} + +// f = c'*x + I(x >= 0) +case class ProximalLp(c: DenseVector[Double]) extends Proximal { + def prox(x: DenseVector[Double], rho: Double) = { + cforRange(0 until x.length) { i => + x.update(i, max(0, x(i) - c(i) / rho)) + } + } +} diff --git a/math/src/main/scala/breeze/optimize/proximal/QuadraticMinimizer.scala b/math/src/main/scala/breeze/optimize/proximal/QuadraticMinimizer.scala new file mode 100644 index 000000000..52ec56146 --- /dev/null +++ b/math/src/main/scala/breeze/optimize/proximal/QuadraticMinimizer.scala @@ -0,0 +1,538 @@ +package breeze.optimize.proximal + +import breeze.linalg.cholesky +import breeze.linalg.LU +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 +import breeze.util.Implicits._ +import spire.syntax.cfor._ + +/** + * Proximal operators and ADMM based Primal-Dual QP Solver + * + * Reference: http://www.stanford.edu/~boyd/papers/admm/quadprog/quadprog.html + * + * + * It solves problem that has the following structure + * + * 1/2 x'Hx + f'x + g(x) + * s.t Aeqx = b + * + * g(x) represents the following constraints which covers ALS based matrix factorization use-cases + * + * 1. x >= 0 + * 2. lb <= x <= ub + * 3. L1(x) + * 4. L2(x) + * 5. Generic regularization on x + * + * @param nGram rank of dense gram matrix + * @param proximal proximal operator to be used + * @param Aeq rhs matrix for equality constraints + * @param beq lhs constants for equality constraints + * @author debasish83 + */ +class QuadraticMinimizer(nGram: Int, + proximal: Proximal = null, + Aeq: DenseMatrix[Double] = null, + beq: DenseVector[Double] = null, + maxIters: Int = -1, abstol: Double = 1e-6, reltol: Double = 1e-4) + extends SerializableLogging { + + type BDM = DenseMatrix[Double] + type BDV = DenseVector[Double] + + case class State private[QuadraticMinimizer](x: BDV, u: BDV, z: BDV, R: BDM, pivot: Array[Int], xHat: BDV, zOld: BDV, residual: BDV, s: BDV, iter: Int, converged: Boolean) + + val linearEquality = if (Aeq != null) Aeq.rows else 0 + + if(linearEquality > 0) + require(beq.length == linearEquality, s"QuadraticMinimizer linear equalities should match beq vector") + + val n = nGram + linearEquality + + //TO DO: Tune alpha based on Nesterov's acceleration + val alpha: Double = 1.0 + + /** + * wsH is the workspace for gram matrix / quasi definite system based on the problem definition + * Quasi definite system is formed if we have a set of affine constraints in quadratic minimization + * minimize 0.5x'Hx + c'x + g(x) + * s.t Aeq x = beq + * + * Proximal formulations can handle different g(x) as specified above but if Aeq x = beq is not + * a probability simplex (1'x = s, x >=0) and arbitrary equality constraint, then we form the + * quasi definite system by adding affine constraint Aeq x = beq in x-update and handling bound + * constraints like lb <= z <= ub in g(z) update. + * + * Capability of adding affine constraints let us support all the features that an interior point + * method like Matlab QuadProg and Mosek can handle + * + * The affine system looks as follows: + * + * Aeq is l x rank + * H is rank x rank + * + * wsH is a quasi-definite matrix + * [ P + rho*I, Aeq' ] + * [ Aeq , 0 ] + * + * if wsH is positive definite then we use cholesky factorization and cache the factorization + * if wsH is quasi definite then we use lu factorization and cache the factorization + * + * Preparing the wsH workspace is expensive and so memory is allocated at the construct time and + * an API updateGram is exposed to users so that the workspace can be used by many solves that show + * up in distributed frameworks like Spark mllib ALS + */ + + private val wsH = + if (linearEquality > 0) { + val ws = DenseMatrix.zeros[Double](n, n) + val transAeq = Aeq.t + ws(nGram until (nGram + Aeq.rows), 0 until Aeq.cols) := Aeq + ws(0 until nGram, nGram until (nGram + Aeq.rows)) := transAeq + ws + } else { + DenseMatrix.zeros[Double](n, n) + } + + val admmIters = if (maxIters < 0) Math.max(400, 20 * n) else maxIters + + def getProximal = proximal + + def updateGram(row: Int, col: Int, value: Double) { + if (row < 0 || row >= n) + throw new IllegalArgumentException("QuadraticMinimizer row out of bounds for gram matrix update") + if (col < 0 || col >= n) + throw new IllegalArgumentException("QuadraticMinimizer column out of bounds for gram matrix update") + wsH.update(row, col, value) + } + + //u is the langrange multiplier + //z is for the proximal operator application + + private def initialState(nGram: Int) = { + var R: DenseMatrix[Double] = null + var pivot: Array[Int] = null + + //Dense cholesky factorization if the gram matrix is well defined + if (linearEquality > 0) { + val lu= LU(wsH) + R = lu._1 + pivot = lu._2 + } else { + R = cholesky(wsH).t + } + + val x = DenseVector.zeros[Double](nGram) + + val z = DenseVector.zeros[Double](nGram) + val u = DenseVector.zeros[Double](nGram) + + val xHat = DenseVector.zeros[Double](nGram) + val zOld = DenseVector.zeros[Double](nGram) + + val residual = DenseVector.zeros[Double](nGram) + val s = DenseVector.zeros[Double](nGram) + + State(x, u, z, R, pivot, xHat, zOld, residual, s, 0, false) + } + + def iterations(q: DenseVector[Double], + rho: Double): Iterator[State] = Iterator.iterate(initialState(nGram)) { state => + import state._ + + //scale will hold q + linearEqualities + val convergenceScale = sqrt(n) + + //scale = rho*(z - u) - q + val scale = DenseVector.zeros[Double](n) + cforRange(0 until z.length) { i => + val entryScale = rho * (z(i) - u(i)) - q(i) + scale.update(i, entryScale) + } + if (linearEquality > 0) + cforRange(0 until beq.length) { i => scale.update(nGram + i, beq(i)) } + + //TO DO : Use LDL' decomposition for efficiency if the Gram matrix is sparse + // If the Gram matrix is positive definite then use Cholesky else use LU Decomposition + val xlambda = if (linearEquality > 0) { + // x = U \ (L \ q) + QuadraticMinimizer.solveTriangularLU(R, pivot, scale) + } else { + // x = R \ (R' \ scale) + //Step 1 : R' * y = scale + //Step 2 : R * x = y + QuadraticMinimizer.solveTriangular(R, scale) + } + cforRange(0 until x.length) {i => x.update(i, xlambda(i))} + + //Unconstrained Quadratic Minimization does need any proximal step + if (proximal == null) { + State(x, u, z, R, pivot, xHat, zOld, residual, s, iter + 1, true) + } + else { + //z-update with relaxation + + //zold = (1-alpha)*z + //x_hat = alpha*x + zold + zOld := z + zOld *= 1 - alpha + + xHat := 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 := x + residual -= z + val residualNorm = norm(residual, 2) + + //history.s_norm(k) = norm(-rho*(z - zold)) + s := z + s -= zOld + s *= -rho + val sNorm = norm(s, 2) + + //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 = convergenceScale * abstol + reltol * max(norm(x, 2), norm(residual, 2)) + val epsDual = convergenceScale * abstol + reltol * norm(s, 2) + + val converged = residualNorm < epsPrimal && sNorm < epsDual || iter > admmIters + + State(x, u, z, R, pivot, xHat, zOld, residual, s, iter + 1, converged) + } + }.takeUpToWhere(_.converged) + + private def computeRho(H: DenseMatrix[Double]): Double = { + proximal match { + case null => 0.0 + case ProximalL1() => { + val eigenMax = QuadraticMinimizer.normColumn(H) + val eigenMin = QuadraticMinimizer.approximateMinEigen(H) + sqrt(eigenMin * eigenMax) + } + case _ => sqrt(QuadraticMinimizer.normColumn(H)) + } + } + + def minimize(H: DenseMatrix[Double], q: DenseVector[Double]): DenseVector[Double] = { + minimizeAndReturnState(H, q).x + } + + def minimizeAndReturnState(H: DenseMatrix[Double], q: DenseVector[Double]): State = { + iterations(H, q).last + } + + def minimizeAndReturnState(q: DenseVector[Double]) : State = { + iterations(q).last + } + + def minimize(q: DenseVector[Double]): DenseVector[Double] = { + minimizeAndReturnState(q).x + } + + def iterations(q: DenseVector[Double]) : Iterator[State] = { + val rho = computeRho(wsH) + cforRange(0 until q.length) {i => wsH.update(i, i, wsH(i, i) + rho)} + iterations(q, rho) + } + + def iterations(H: DenseMatrix[Double], q: DenseVector[Double]): Iterator[State] = { + wsH(0 until H.rows, 0 until H.cols) := H + iterations(q) + } +} + +/** + * 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 = { + var absColSum = 0.0 + var maxColSum = 0.0 + for (c <- 0 until H.cols) { + for (r <- 0 until H.rows) { + absColSum += abs(H(r, c)) + } + if (absColSum > maxColSum) maxColSum = absColSum + absColSum = 0.0 + } + maxColSum + } + + //approximate max eigen using inverse power method + def approximateMaxEigen(H: DenseMatrix[Double]) : Double = { + val pm = new PowerMethod[DenseVector[Double], DenseMatrix[Double]]() + val init = DenseVector.rand[Double](H.rows, Rand.gaussian(0, 1)) + pm.eigen(init, H) + } + + //approximate min eigen using inverse power method + def approximateMinEigen(H: DenseMatrix[Double]) : Double = { + val R = cholesky(H).t + val pmInv = PowerMethod.inverse() + val init = DenseVector.rand[Double](H.rows, Rand.gaussian(0, 1)) + 1.0/pmInv.eigen(init, R) + } + + /* + * Triangular LU solve for A*X = B + * TO DO : Add appropriate exception from LAPACK + */ + def solveTriangularLU(A: DenseMatrix[Double], pivot: Array[Int], B: DenseVector[Double]) : DenseVector[Double] = { + require(A.rows == A.cols) + + val X = new DenseMatrix(B.length, 1, B.data.clone) + + val n = A.rows + val nrhs = X.cols + var info: intW = new intW(0) + + lapack.dgetrs("No transpose", n, nrhs, A.data, 0, A.rows, pivot, 0, X.data, 0, X.rows, info) + + if (info.`val` > 0) throw new LapackException("DGETRS: LU solve unsuccessful") + + DenseVector(X.data) + } + + /*Triangular cholesky solve for A*X = B */ + def solveTriangular(A: DenseMatrix[Double], B: DenseVector[Double]) : DenseVector[Double] = { + require(A.rows == A.cols) + + val X = new DenseMatrix(B.length, 1, B.data.clone) + + val n = A.rows + val nrhs = X.cols + var info: intW = new intW(0) + + lapack.dpotrs("L", n, nrhs, A.data, 0, A.rows, X.data, 0, X.rows, info) + + if (info.`val` > 0) throw new LapackException("DPOTRS : Leading minor of order i of A is not positive definite.") + + DenseVector(X.data) + } + + def apply(rank: Int, + constraint: Constraint, + lambda: Double): QuadraticMinimizer = { + constraint match { + case POSITIVE => new QuadraticMinimizer(rank, ProjectPos()) + case BOX => { + //Direct QP with bounds + val lb = DenseVector.zeros[Double](rank) + val ub = DenseVector.ones[Double](rank) + new QuadraticMinimizer(rank, ProjectBox(lb, ub)) + } + case EQUALITY => { + //Direct QP with equality and positivity constraint + val Aeq = DenseMatrix.ones[Double](1, rank) + val beq = DenseVector.ones[Double](1) + new QuadraticMinimizer(rank, ProjectPos(), Aeq, beq) + } + case SPARSE => new QuadraticMinimizer(rank, ProximalL1().setLambda(lambda)) + } + } + + def computeObjective(h: DenseMatrix[Double], q: DenseVector[Double], x: DenseVector[Double]): Double = { + val res = (x.t*h*x)*0.5 + q.dot(x) + res + } + + 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) + } + + def main(args: Array[String]) { + if (args.length < 4) { + println("Usage: QpSolver n m lambda beta") + println("Test QpSolver with a simple quadratic function of dimension n and m equalities lambda beta for elasticNet") + sys.exit(1) + } + + val problemSize = args(0).toInt + val nequalities = args(1).toInt + + val lambda = args(2).toDouble + val beta = args(3).toDouble + + println(s"Generating randomized QPs with rank ${problemSize} equalities ${nequalities}") + val (aeq, b, bl, bu, q, h) = QpGenerator(problemSize, nequalities) + + println(s"Test QuadraticMinimizer, CG , BFGS and OWLQN with $problemSize variables and $nequalities equality constraints") + + val luStart = System.nanoTime() + val luResult = h \ q:*(-1.0) + val luTime = System.nanoTime() - luStart + + val cg = new ConjugateGradient[DenseVector[Double], DenseMatrix[Double]]() + + val startCg = System.nanoTime() + val cgResult = cg.minimize(q:*(-1.0), h) + val cgTime = System.nanoTime() - startCg + + val qpSolver = new QuadraticMinimizer(problemSize) + val qpStart = System.nanoTime() + val result = qpSolver.minimize(h, q) + val qpTime = System.nanoTime() - qpStart + + val startBFGS = System.nanoTime() + val bfgsResult = optimizeWithLBFGS(DenseVector.rand[Double](problemSize), h, q) + val bfgsTime = System.nanoTime() - startBFGS + + println(s"||qp - lu|| norm ${norm(result - luResult, 2)} max-norm ${norm(result - luResult, inf)}") + println(s"||cg - lu|| norm ${norm(cgResult - luResult,2)} max-norm ${norm(cgResult - luResult, inf)}") + println(s"||bfgs - lu|| norm ${norm(bfgsResult - luResult, 2)} max-norm ${norm(bfgsResult - luResult, inf)}") + + val luObj = computeObjective(h, q, luResult) + val bfgsObj = computeObjective(h, q, bfgsResult) + val qpObj = computeObjective(h, q, result) + + println(s"Objective lu $luObj bfgs $bfgsObj qp $qpObj") + + println(s"dim ${problemSize} lu ${luTime/1e6} ms qp ${qpTime/1e6} ms cg ${cgTime/1e6} ms bfgs ${bfgsTime/1e6} ms") + + val lambdaL1 = lambda * beta + val lambdaL2 = lambda * (1 - beta) + + 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 owlqnTime = System.nanoTime() - startOWLQN + + println(s"||owlqn - sparseqp|| norm ${norm(owlqnResult.x - sparseQpResult.x, 2)} inf-norm ${norm(owlqnResult.x - sparseQpResult.x, inf)}") + println(s"sparseQp ${sparseQpTime/1e6} ms iters ${sparseQpResult.iter} owlqn ${owlqnTime/1e6} ms iters ${owlqnResult.iter}") + + val posQp = QuadraticMinimizer(h.rows, POSITIVE, 0.0) + val posQpStart = System.nanoTime() + val posQpResult = posQp.minimizeAndReturnState(h, q) + val posQpTime = System.nanoTime() - posQpStart + + val nnls = new NNLS() + val nnlsStart = System.nanoTime() + val nnlsResult = nnls.minimizeAndReturnState(h, q) + val nnlsTime = System.nanoTime() - nnlsStart + + println(s"posQp ${posQpTime/1e6} ms iters ${posQpResult.iter} nnls ${nnlsTime/1e6} ms iters ${nnlsResult.iter}") + + val boundsQp = new QuadraticMinimizer(h.rows, ProjectBox(bl, bu)) + val boundsQpStart = System.nanoTime() + val boundsQpResult = boundsQp.minimizeAndReturnState(h, q) + val boundsQpTime = System.nanoTime() - boundsQpStart + + println(s"boundsQp ${boundsQpTime/1e6} ms iters ${boundsQpResult.iter} converged ${boundsQpResult.converged}") + + val qpEquality = new QuadraticMinimizer(h.rows, ProjectPos(), aeq, b) + val qpEqualityStart = System.nanoTime() + val qpEqualityResult = qpEquality.minimizeAndReturnState(h, q) + val qpEqualityTime = System.nanoTime() - qpEqualityStart + + println(s"Qp Equality ${qpEqualityTime/1e6} ms iters ${qpEqualityResult.iter} converged ${qpEqualityResult.converged}") + } +} \ No newline at end of file diff --git a/math/src/test/scala/breeze/optimize/linear/NNLSTest.scala b/math/src/test/scala/breeze/optimize/linear/NNLSTest.scala new file mode 100644 index 000000000..76ad22923 --- /dev/null +++ b/math/src/test/scala/breeze/optimize/linear/NNLSTest.scala @@ -0,0 +1,81 @@ +package breeze.optimize.linear + +import breeze.linalg.{DenseMatrix, DenseVector, norm} +import breeze.numerics._ +import breeze.optimize.OptimizeTestBase + +class NNLSTest extends OptimizeTestBase { + /** Generate an NNLS problem whose optimal solution is the all-ones vector. */ + def genOnesData(n: Int): (DenseMatrix[Double], DenseVector[Double]) = { + val A = DenseMatrix.rand[Double](n, n) + val b = A*DenseVector.ones[Double](n) + + val ata = A.t*A + val atb = A.t*b + + (ata, atb) + } + + test("NNLS: exact solution cases") { + val n = 20 + val nnls = new NNLS() + var numSolved = 0 + + // About 15% of random 20x20 [-1,1]-matrices have a singular value less than 1e-3. NNLS + // can legitimately fail to solve these anywhere close to exactly. So we grab a considerable + // sample of these matrices and make sure that we solved a substantial fraction of them. + for (k <- 0 until 100) { + val (ata, atb) = genOnesData(n) + val x = nnls.minimize(ata, atb) + atb *= -1.0 + val golden = DenseVector.ones[Double](n) + x -= golden + if ((norm(x, 2) < 1e-2) && (norm(x, inf) < 1e-3)) numSolved = numSolved + 1 + } + assert(numSolved > 50) + } + + test("NNLS: nonnegativity constraint active") { + val n = 5 + 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[Double](-1.632, 2.115, 1.094, -1.025, -0.636) + + val goodx = Array(0.13025, 0.54506, 0.2874, 0.0, 0.028628) + + val nnls = new NNLS() + val x = nnls.minimize(ata, atb) + for (i <- 0 until n) { + assert(abs(x(i) - goodx(i)) < 1E-3) + assert(x(i) >= 0) + } + } + + test("NNLS: objective value test") { + val n = 5 + val ata = new DenseMatrix[Double](5, 5, + Array(517399.13534, 242529.67289, -153644.98976, 130802.84503, -798452.29283 + , 242529.67289, 126017.69765, -75944.21743, 81785.36128, -405290.60884 + , -153644.98976, -75944.21743, 46986.44577, -45401.12659, 247059.51049 + , 130802.84503, 81785.36128, -45401.12659, 67457.31310, -253747.03819 + , -798452.29283, -405290.60884, 247059.51049, -253747.03819, 1310939.40814) + ) + + val atb = DenseVector(-31755.05710, 13047.14813, -20191.24443, 25993.77580, 11963.55017) + + /** reference solution obtained from matlab function quadprog */ + val refx = DenseVector(34.90751, 103.96254, 0.00000, 27.82094, 58.79627) + val refObj = NNLS.computeObjectiveValue(ata, atb, refx) + + val nnls = new NNLS() + val x = nnls.minimize(ata, atb) + val obj = NNLS.computeObjectiveValue(ata, atb, x) + + assert(obj < refObj + 1E-5) + } +} diff --git a/math/src/test/scala/breeze/optimize/linear/PowerMethodTest.scala b/math/src/test/scala/breeze/optimize/linear/PowerMethodTest.scala new file mode 100644 index 000000000..cf54f5045 --- /dev/null +++ b/math/src/test/scala/breeze/optimize/linear/PowerMethodTest.scala @@ -0,0 +1,42 @@ +package breeze.optimize.linear + +import breeze.numerics.abs +import org.junit.runner.RunWith +import org.scalatest._ +import org.scalatest.junit._ +import breeze.linalg._ + +@RunWith(classOf[JUnitRunner]) +class PowerMethodTest extends FunSuite { + val n = 5 + val gram = new DenseMatrix[Double](n, n, Array(5.6880,-4.5286,6.7923,0.3049,-6.0388, + -4.5286,8.0638,-6.9012,-2.6776,6.1795, + 6.7923,-6.9012,12.5510,-1.1917,-8.3500, + 0.3049,-2.6776,-1.1917,4.0684,-1.7535, + -6.0388,6.1795,-8.3500,-1.7535,8.2831)) + val init = DenseVector(0.1770,0.2505,1.5957,0.7204,0.9246) + val eigs = eigSym(gram) + + test("max eigen value from power method approximately equal to eigSym max") { + val eigenGold = max(eigs.eigenvalues) + val pm = new PowerMethod[DenseVector[Double], DenseMatrix[Double]]() + val eigenApprox = pm.eigen(init, gram) + assert(abs(eigenGold - eigenApprox) < 1e-3) + } + + test("min eigen value from power method approximately equal to eigSym min") { + val eigenGold = min(eigs.eigenvalues) + val pm = new PowerMethod[DenseVector[Double], DenseMatrix[Double]]() + val inverseGram = gram \ DenseMatrix.eye[Double](gram.rows) + val eigenApprox = 1.0/pm.eigen(init, inverseGram) + assert(abs(eigenGold - eigenApprox) < 1e-3) + } + + test("min eigen value from inverse power method approximately equal to eigSym min") { + val eigenGold = min(eigs.eigenvalues) + val pmInv = PowerMethod.inverse(10, 1e-5) + val R = cholesky(gram).t + val eigenApprox = 1.0/pmInv.eigen(init, R) + assert(abs(eigenGold - eigenApprox) < 1e-3) + } +} diff --git a/math/src/test/scala/breeze/optimize/proximal/QuadraticMinimizerTest.scala b/math/src/test/scala/breeze/optimize/proximal/QuadraticMinimizerTest.scala new file mode 100644 index 000000000..ca345f8c1 --- /dev/null +++ b/math/src/test/scala/breeze/optimize/proximal/QuadraticMinimizerTest.scala @@ -0,0 +1,312 @@ +package breeze.optimize.proximal + +import breeze.numerics.abs +import breeze.optimize.linear.NNLS +import org.scalatest._ +import org.scalatest.junit._ +import org.junit.runner.RunWith +import breeze.linalg._ +import breeze.optimize._ +import breeze.optimize.proximal.Constraint._ +import breeze.numerics._ + +@RunWith(classOf[JUnitRunner]) +class QuadraticMinimizerTest extends OptimizeTestBase with Matchers { + def matricesNearlyEqual(A: DenseMatrix[Double], B: DenseMatrix[Double], threshold: Double = 1E-6) { + for (i <- 0 until A.rows; j <- 0 until A.cols) + A(i, j) should be(B(i, j) +- threshold) + } + + test("Randomly generated matrix vector solve") { + val H = DenseMatrix((2.4779, 1.9284, 1.4078), (1.9284, 3.0218, 2.2968), (1.4078, 2.2968, 2.2033)) + val q = DenseVector(0.14678, 0.32895, 0.71402) + val x = H \ q + + assert(norm(x - DenseVector(-0.0035577, -0.6592646, 1.0135739), 2) < 1e-4) + } + + test("cholesky factorization based forward-backward solve") { + val H = QpGenerator.getGram(3) + val R = cholesky(H) + val q = DenseVector.rand[Double](3) + val goldenBefore = H \ q + val x = QuadraticMinimizer.solveTriangular(R, q) + val goldenAfter = H \ q + + assert(norm(goldenBefore.toDenseVector - x, inf) < 1E-5) + assert(norm(goldenAfter.toDenseVector - x, inf) < 1E-5) + } + + test("lu factorization based forward-backward solve") { + val H = QpGenerator.getGram(5) + val lu = LU(H) + val q = DenseVector.rand[Double](5) + val x = QuadraticMinimizer.solveTriangularLU(lu._1, lu._2, q) + val golden = H \ q + assert(norm(golden - x) < 1E-8) + } + + val n = 5 + val gram = new DenseMatrix[Double](n, n, Array(5.6880,-4.5286,6.7923,0.3049,-6.0388, + -4.5286,8.0638,-6.9012,-2.6776,6.1795, + 6.7923,-6.9012,12.5510,-1.1917,-8.3500, + 0.3049,-2.6776,-1.1917,4.0684,-1.7535, + -6.0388,6.1795,-8.3500,-1.7535,8.2831)) + val init = DenseVector(0.1770,0.2505,1.5957,0.7204,0.9246) + val eigs = eigSym(gram) + + test("min eigen computed using inverse power law approximately same as min eigen") { + val eigs = eigSym(gram) + val eigenMin = min(eigs.eigenvalues) + val approxEigenMin = QuadraticMinimizer.approximateMinEigen(gram) + assert(abs(eigenMin - approxEigenMin) < 1e-3) + } + + test("max eigen computed using power law approximately same as max eigen") { + val eigs = eigSym(gram) + val eigenMax = max(eigs.eigenvalues) + val approxEigenMax = QuadraticMinimizer.approximateMaxEigen(gram) + assert(abs(eigenMax - approxEigenMax) < 1e-3) + } + + test("Unconstrained Quadratic Minimization compared to LU Solve") { + val problemSize = 10 + + val H = DenseMatrix.eye[Double](problemSize) :* 2.0 + val f = DenseVector.ones[Double](problemSize) :* 6.0 + + val dposvResult = H \ f + + val qpSolver = new QuadraticMinimizer(problemSize) + val result = qpSolver.minimize(H, f :* (-1.0)) + + assert(norm(result - dposvResult, 2) < 1E-4) + } + + test("Unconstrained Quadratic Minimization compared to BFGS") { + val problemSize = 10 + val lbfgs = new LBFGS[DenseVector[Double]](-1, 7) + + def optimizeWithLBFGS(init: DenseVector[Double]) = { + val f = new DiffFunction[DenseVector[Double]] { + def calculate(x: DenseVector[Double]) = { + (norm((x - 3.0) :^ 2.0, 1), (x * 2.0) - 6.0) + } + } + lbfgs.minimize(f, init) + } + + val init = DenseVector.zeros[Double](problemSize) + init(0 until init.length by 2) := -1.2 + init(1 until init.length by 2) := 1.0 + + val bfgsResult = optimizeWithLBFGS(init) + + val H = DenseMatrix.eye[Double](problemSize) :* 2.0 + val f = DenseVector.ones[Double](problemSize) :* 6.0 + + val qpSolver = new QuadraticMinimizer(problemSize) + val result = qpSolver.minimize(H, f :* (-1.0)) + + assert(norm(result - bfgsResult, 2) < 1E-4) + } + + test("Quadratic Minimization with L1 compared to OWLQN") { + val problemSize = 10 + + val owlqn = new OWLQN[Int, DenseVector[Double]](100, 4) + + def optimizeWithOWLQN(init: DenseVector[Double]) = { + val f = new DiffFunction[DenseVector[Double]] { + def calculate(x: DenseVector[Double]) = { + (norm((x - 3.0) :^ 2.0, 1), (x * 2.0) - 6.0) + } + } + val result = owlqn.minimize(f, init) + norm(result - 2.5, 2) < 1E-10 + } + + val init = DenseVector.zeros[Double](problemSize) + init(0 until init.length by 2) := -1.2 + init(1 until init.length by 2) := 1.0 + + assert(optimizeWithOWLQN(init)) + + val H = DenseMatrix.eye[Double](problemSize) :* 2.0 + val f = DenseVector.ones[Double](problemSize) :* (-6.0) + + val qpSolverL1 = QuadraticMinimizer(problemSize, SPARSE, 1.0) + + val l1Result = qpSolverL1.minimize(H, f) + + val normL1 = norm(l1Result - 2.5, 2) + assert(normL1 < 1E-3) + } + + test("Quadratic Minimization with positivity compared to NNLS") { + val n = 5 + 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 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 + val posResult = qpSolverPos.minimizeAndReturnState(ata, atb) + assert(posResult.converged) + assert(norm(posResult.x - goodx, 2) < 1E-3) + } + + test("Quadratic Minimization with bounds compared to Octave with NNLS Example") { + val n = 5 + 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 goodBounds = DenseVector(0.0, 0.25000000000236045, 0.2499999999945758, 0.0, 0.0) + + val lb = DenseVector.zeros[Double](5) + val ub = DenseVector.ones[Double](5) :* 0.25 + + val qpSolverBounds = new QuadraticMinimizer(5, ProjectBox(lb, ub)) + val boundsResult = qpSolverBounds.minimizeAndReturnState(ata, atb :* (-1.0)) + + assert(boundsResult.converged) + assert(norm(boundsResult.x - goodBounds) < 1E-4) + } + + test("Quadratic Minimization with Equality compared to Octave MovieLens Example") { + 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 ml = 25 + + val directQpMl = QuadraticMinimizer(ml, EQUALITY, 0.0) + val directQpResult = directQpMl.minimizeAndReturnState(Hml, fml) + + 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) + + assert(directQpResult.converged) + assert(norm(directQpResult.x - golden) < 1e-3) + assert(abs(sum(directQpResult.x) - 1.0) < 1e-4) + } + + test("Quadratic Minimization with L1 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) + + val qpSolverMlL1 = QuadraticMinimizer(25, SPARSE, 2.0) + + val qpMlL1Result = qpSolverMlL1.minimizeAndReturnState(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) + + assert(qpMlL1Result.converged) + assert(norm(qpMlL1Result.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)) + val q1 = DenseVector(-1880.240029, -1920.949941, -2030.476172, -1956.642164, -2021.502985, -2090.503157, -1965.934820, -2072.855628, -2098.075034, -2035.059185, -1999.005923, -2121.515181, -1944.759586, -2035.397706, -1939.872057, -1888.635008, -1986.031605, -1973.738457, -2024.468051, -2003.765736) + + val qpIters = QuadraticMinimizer(20, POSITIVE, 0.0) + val qpItersResult = qpIters.minimizeAndReturnState(P1, q1) + + val nnls = new NNLS() + val nnlsResult = nnls.minimize(P1, q1 :* (-1.0)) + + //NNLS test + val P2 = new DenseMatrix[Double](20, 20, Array(333907.312770, -60814.043975, 207935.829941, -162881.367739, -43730.396770, 17511.428983, -243340.496449, -225245.957922, 104700.445881, 32430.845099, 336378.693135, -373497.970207, -41147.159621, 53928.060360, -293517.883778, 53105.278068, 0.000000, -85257.781696, 84913.970469, -10584.080103, -60814.043975, 13826.806664, -38032.612640, 33475.833875, 10791.916809, -1040.950810, 48106.552472, 45390.073380, -16310.282190, -2861.455903, -60790.833191, 73109.516544, 9826.614644, -8283.992464, 56991.742991, -6171.366034, 0.000000, 19152.382499, -13218.721710, 2793.734234, 207935.829941, -38032.612640, 129661.677608, -101682.098412, -27401.299347, 10787.713362, -151803.006149, -140563.601672, 65067.935324, 20031.263383, 209521.268600, -232958.054688, -25764.179034, 33507.951918, -183046.845592, 32884.782835, 0.000000, -53315.811196, 52770.762546, -6642.187643, -162881.367739, 33475.833875, -101682.098412, 85094.407608, 25422.850782, -5437.646141, 124197.166330, 116206.265909, -47093.484134, -11420.168521, -163429.436848, 189574.783900, 23447.172314, -24087.375367, 148311.355507, -20848.385466, 0.000000, 46835.814559, -38180.352878, 6415.873901, -43730.396770, 10791.916809, -27401.299347, 25422.850782, 8882.869799, 15.638084, 35933.473986, 34186.371325, -10745.330690, -974.314375, -43537.709621, 54371.010558, 7894.453004, -5408.929644, 42231.381747, -3192.010574, 0.000000, 15058.753110, -8704.757256, 2316.581535, 17511.428983, -1040.950810, 10787.713362, -5437.646141, 15.638084, 2794.949847, -9681.950987, -8258.171646, 7754.358930, 4193.359412, 18052.143842, -15456.096769, -253.356253, 4089.672804, -12524.380088, 5651.579348, 0.000000, -1513.302547, 6296.461898, 152.427321, -243340.496449, 48106.552472, -151803.006149, 124197.166330, 35933.473986, -9681.950987, 182931.600236, 170454.352953, -72361.174145, -19270.461728, -244518.179729, 279551.060579, 33340.452802, -37103.267653, 219025.288975, -33687.141423, 0.000000, 67347.950443, -58673.009647, 8957.800259, -225245.957922, 45390.073380, -140563.601672, 116206.265909, 34186.371325, -8258.171646, 170454.352953, 159322.942894, -66074.960534, -16839.743193, -226173.967766, 260421.044094, 31624.194003, -33839.612565, 203889.695169, -30034.828909, 0.000000, 63525.040745, -53572.741748, 8575.071847, 104700.445881, -16310.282190, 65067.935324, -47093.484134, -10745.330690, 7754.358930, -72361.174145, -66074.960534, 35869.598076, 13378.653317, 106033.647837, -111831.682883, -10455.465743, 18537.392481, -88370.612394, 20344.288488, 0.000000, -22935.482766, 29004.543704, -2409.461759, 32430.845099, -2861.455903, 20031.263383, -11420.168521, -974.314375, 4193.359412, -19270.461728, -16839.743193, 13378.653317, 6802.081898, 33256.395091, -30421.985199, -1296.785870, 7026.518692, -24443.378205, 9221.982599, 0.000000, -4088.076871, 10861.014242, -25.092938, 336378.693135, -60790.833191, 209521.268600, -163429.436848, -43537.709621, 18052.143842, -244518.179729, -226173.967766, 106033.647837, 33256.395091, 339200.268106, -375442.716811, -41027.594509, 54636.778527, -295133.248586, 54177.278365, 0.000000, -85237.666701, 85996.957056, -10503.209968, -373497.970207, 73109.516544, -232958.054688, 189574.783900, 54371.010558, -15456.096769, 279551.060579, 260421.044094, -111831.682883, -30421.985199, -375442.716811, 427793.208465, 50528.074431, -57375.986301, 335203.382015, -52676.385869, 0.000000, 102368.307670, -90679.792485, 13509.390393, -41147.159621, 9826.614644, -25764.179034, 23447.172314, 7894.453004, -253.356253, 33340.452802, 31624.194003, -10455.465743, -1296.785870, -41027.594509, 50528.074431, 7255.977434, -5281.636812, 39298.355527, -3440.450858, 0.000000, 13717.870243, -8471.405582, 2071.812204, 53928.060360, -8283.992464, 33507.951918, -24087.375367, -5408.929644, 4089.672804, -37103.267653, -33839.612565, 18537.392481, 7026.518692, 54636.778527, -57375.986301, -5281.636812, 9735.061160, -45360.674033, 10634.633559, 0.000000, -11652.364691, 15039.566630, -1202.539106, -293517.883778, 56991.742991, -183046.845592, 148311.355507, 42231.381747, -12524.380088, 219025.288975, 203889.695169, -88370.612394, -24443.378205, -295133.248586, 335203.382015, 39298.355527, -45360.674033, 262923.925938, -42012.606885, 0.000000, 79810.919951, -71657.856143, 10464.327491, 53105.278068, -6171.366034, 32884.782835, -20848.385466, -3192.010574, 5651.579348, -33687.141423, -30034.828909, 20344.288488, 9221.982599, 54177.278365, -52676.385869, -3440.450858, 10634.633559, -42012.606885, 13238.686902, 0.000000, -8739.845698, 16511.872845, -530.252003, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 118.430000, 0.000000, 0.000000, 0.000000, -85257.781696, 19152.382499, -53315.811196, 46835.814559, 15058.753110, -1513.302547, 67347.950443, 63525.040745, -22935.482766, -4088.076871, -85237.666701, 102368.307670, 13717.870243, -11652.364691, 79810.919951, -8739.845698, 0.000000, 26878.133950, -18588.407734, 3894.934299, 84913.970469, -13218.721710, 52770.762546, -38180.352878, -8704.757256, 6296.461898, -58673.009647, -53572.741748, 29004.543704, 10861.014242, 85996.957056, -90679.792485, -8471.405582, 15039.566630, -71657.856143, 16511.872845, 0.000000, -18588.407734, 23649.538538, -1951.083671, -10584.080103, 2793.734234, -6642.187643, 6415.873901, 2316.581535, 152.427321, 8957.800259, 8575.071847, -2409.461759, -25.092938, -10503.209968, 13509.390393, 2071.812204, -1202.539106, 10464.327491, -530.252003, 0.000000, 3894.934299, -1951.083671, 738.955915)) + val q2 = DenseVector(31755.057100, -13047.148129, 20191.244430, -25993.775800, -11963.550172, -4272.425977, -33569.856044, -33451.387021, 2320.764250, -5333.136834, 30633.764272, -49513.939049, -10351.230305, 872.276714, -37634.078430, -4628.338543, -0.000000, -18109.093788, 1856.725521, -3397.693211) + + val P3 = new DenseMatrix[Double](25, 25, Array(0.986619, 0.639909, 0.748906, 0.900377, 0.688079, 0.734711, 0.835164, 0.723973, 0.822436, 0.852591, 0.699979, 0.609533, 0.559504, 0.708015, 0.544744, 0.658359, 0.632510, 0.751316, 0.653993, 0.642734, 0.799106, 0.898689, 0.712825, 0.878405, 0.849565, 0.639909, 1.055175, 0.884940, 0.975502, 0.815121, 0.845699, 0.899780, 0.709264, 0.960949, 1.021108, 0.896508, 0.692635, 0.659746, 0.809355, 0.539466, 0.730501, 0.639971, 0.881502, 0.840159, 0.628515, 0.917052, 0.950677, 0.823301, 1.022355, 0.994935, 0.748906, 0.884940, 1.421868, 1.175203, 0.986093, 1.028669, 1.091437, 0.849192, 1.204776, 1.249037, 1.115160, 0.815680, 0.772715, 0.971263, 0.621197, 0.875925, 0.757031, 1.034483, 1.022001, 0.705077, 1.115237, 1.164796, 0.983688, 1.226847, 1.180170, 0.900377, 0.975502, 1.175203, 1.676666, 1.065369, 1.109756, 1.181625, 0.944176, 1.220418, 1.328803, 1.156144, 0.984967, 0.916500, 1.046903, 0.728221, 0.991042, 0.855095, 1.181719, 1.095485, 0.901193, 1.214851, 1.277434, 1.077374, 1.372354, 1.356724, 0.688079, 0.815121, 0.986093, 1.065369, 1.249908, 0.953824, 1.004765, 0.771115, 1.082255, 1.161322, 1.021400, 0.757474, 0.736266, 0.923406, 0.598005, 0.812629, 0.706870, 1.011984, 0.968135, 0.682813, 1.034818, 1.039625, 0.937088, 1.152792, 1.121475, 0.734711, 0.845699, 1.028669, 1.109756, 0.953824, 1.334869, 1.091564, 0.824709, 1.157992, 1.226413, 1.045522, 0.731169, 0.709382, 0.980030, 0.634635, 0.853988, 0.758256, 1.070744, 0.997542, 0.692832, 1.118828, 1.119519, 0.977618, 1.202464, 1.152186, 0.835164, 0.899780, 1.091437, 1.181625, 1.004765, 1.091564, 1.566492, 0.994066, 1.307932, 1.301575, 1.075206, 0.697477, 0.670553, 1.068449, 0.725819, 0.908391, 0.877151, 1.083915, 0.993647, 0.735423, 1.202928, 1.275592, 1.040266, 1.242863, 1.150213, 0.723973, 0.709264, 0.849192, 0.944176, 0.771115, 0.824709, 0.994066, 1.291185, 1.084907, 0.934737, 0.802436, 0.559117, 0.510659, 0.801938, 0.624240, 0.682461, 0.761152, 0.633965, 0.633248, 0.622034, 0.826798, 1.084260, 0.777897, 0.868127, 0.764337, 0.822436, 0.960949, 1.204776, 1.220418, 1.082255, 1.157992, 1.307932, 1.084907, 1.821667, 1.386725, 1.220398, 0.710929, 0.683088, 1.116057, 0.719380, 0.925482, 0.907590, 1.010205, 1.041999, 0.673723, 1.224967, 1.373886, 1.082813, 1.261720, 1.130003, 0.852591, 1.021108, 1.249037, 1.328803, 1.161322, 1.226413, 1.301575, 0.934737, 1.386725, 1.831419, 1.289097, 0.885418, 0.868685, 1.190024, 0.740180, 1.034238, 0.884778, 1.342867, 1.252762, 0.811797, 1.374303, 1.319088, 1.190385, 1.478896, 1.427420, 0.699979, 0.896508, 1.115160, 1.156144, 1.021400, 1.045522, 1.075206, 0.802436, 1.220398, 1.289097, 1.501713, 0.831799, 0.807718, 0.994071, 0.599941, 0.876646, 0.730630, 1.074102, 1.088886, 0.677923, 1.126406, 1.127440, 1.013104, 1.257194, 1.217591, 0.609533, 0.692635, 0.815680, 0.984967, 0.757474, 0.731169, 0.697477, 0.559117, 0.710929, 0.885418, 0.831799, 1.206659, 0.820404, 0.673313, 0.477307, 0.698907, 0.529141, 0.853290, 0.812957, 0.721698, 0.775190, 0.790003, 0.746242, 0.994816, 1.053857, 0.559504, 0.659746, 0.772715, 0.916500, 0.736266, 0.709382, 0.670553, 0.510659, 0.683088, 0.868685, 0.807718, 0.820404, 1.103586, 0.666515, 0.455429, 0.667077, 0.500250, 0.855895, 0.807716, 0.674035, 0.759131, 0.732206, 0.732658, 0.965557, 1.021876, 0.708015, 0.809355, 0.971263, 1.046903, 0.923406, 0.980030, 1.068449, 0.801938, 1.116057, 1.190024, 0.994071, 0.673313, 0.666515, 1.293629, 0.631919, 0.822011, 0.745899, 1.060809, 0.965440, 0.672388, 1.090059, 1.069719, 0.959136, 1.162012, 1.109763, 0.544744, 0.539466, 0.621197, 0.728221, 0.598005, 0.634635, 0.725819, 0.624240, 0.719380, 0.740180, 0.599941, 0.477307, 0.455429, 0.631919, 0.802572, 0.551340, 0.549393, 0.654652, 0.567237, 0.531461, 0.684655, 0.749979, 0.627718, 0.744096, 0.710310, 0.658359, 0.730501, 0.875925, 0.991042, 0.812629, 0.853988, 0.908391, 0.682461, 0.925482, 1.034238, 0.876646, 0.698907, 0.667077, 0.822011, 0.551340, 1.075722, 0.642962, 0.955470, 0.862026, 0.658120, 0.955639, 0.944118, 0.835790, 1.057347, 1.042521, 0.632510, 0.639971, 0.757031, 0.855095, 0.706870, 0.758256, 0.877151, 0.761152, 0.907590, 0.884778, 0.730630, 0.529141, 0.500250, 0.745899, 0.549393, 0.642962, 0.976258, 0.726041, 0.658693, 0.579842, 0.811285, 0.917944, 0.731996, 0.859613, 0.798469, 0.751316, 0.881502, 1.034483, 1.181719, 1.011984, 1.070744, 1.083915, 0.633965, 1.010205, 1.342867, 1.074102, 0.853290, 0.855895, 1.060809, 0.654652, 0.955470, 0.726041, 1.779551, 1.210339, 0.817870, 1.289494, 1.030990, 1.082950, 1.419900, 1.451896, 0.653993, 0.840159, 1.022001, 1.095485, 0.968135, 0.997542, 0.993647, 0.633248, 1.041999, 1.252762, 1.088886, 0.812957, 0.807716, 0.965440, 0.567237, 0.862026, 0.658693, 1.210339, 1.442263, 0.688333, 1.139661, 0.988895, 0.991869, 1.276105, 1.281761, 0.642734, 0.628515, 0.705077, 0.901193, 0.682813, 0.692832, 0.735423, 0.622034, 0.673723, 0.811797, 0.677923, 0.721698, 0.674035, 0.672388, 0.531461, 0.658120, 0.579842, 0.817870, 0.688333, 1.044625, 0.753991, 0.788593, 0.710741, 0.906241, 0.935628, 0.799106, 0.917052, 1.115237, 1.214851, 1.034818, 1.118828, 1.202928, 0.826798, 1.224967, 1.374303, 1.126406, 0.775190, 0.759131, 1.090059, 0.684655, 0.955639, 0.811285, 1.289494, 1.139661, 0.753991, 1.621034, 1.201765, 1.081841, 1.365902, 1.323655, 0.898689, 0.950677, 1.164796, 1.277434, 1.039625, 1.119519, 1.275592, 1.084260, 1.373886, 1.319088, 1.127440, 0.790003, 0.732206, 1.069719, 0.749979, 0.944118, 0.917944, 1.030990, 0.988895, 0.788593, 1.201765, 1.692948, 1.052744, 1.265542, 1.169388, 0.712825, 0.823301, 0.983688, 1.077374, 0.937088, 0.977618, 1.040266, 0.777897, 1.082813, 1.190385, 1.013104, 0.746242, 0.732658, 0.959136, 0.627718, 0.835790, 0.731996, 1.082950, 0.991869, 0.710741, 1.081841, 1.052744, 1.291439, 1.189060, 1.159221, 0.878405, 1.022355, 1.226847, 1.372354, 1.152792, 1.202464, 1.242863, 0.868127, 1.261720, 1.478896, 1.257194, 0.994816, 0.965557, 1.162012, 0.744096, 1.057347, 0.859613, 1.419900, 1.276105, 0.906241, 1.365902, 1.265542, 1.189060, 1.847657, 1.521522, 0.849565, 0.994935, 1.180170, 1.356724, 1.121475, 1.152186, 1.150213, 0.764337, 1.130003, 1.427420, 1.217591, 1.053857, 1.021876, 1.109763, 0.710310, 1.042521, 0.798469, 1.451896, 1.281761, 0.935628, 1.323655, 1.169388, 1.159221, 1.521522, 1.887527)) + val q3 = DenseVector(-3.640583, -5.563638, -7.040787, -7.387618, -6.410455, -6.452543, -5.698284, -2.604581, -6.184568, -8.450254, -7.746985, -6.173699, -6.072464, -5.954310, -2.867163, -5.751653, -3.362224, -8.631145, -8.238493, -4.198840, -7.653769, -5.764195, -6.387377, -8.917469, -9.348541) + + val nnlsResult2 = nnls.minimize(P2, q2 :* (-1.0)) + val posResult2 = qpIters.minimizeAndReturnState(P2, q2) + + assert(posResult2.converged) + assert(norm(posResult2.x - nnlsResult2, 2) < 1e-4) + + val nnls2 = new NNLS() + val qpIters2 = QuadraticMinimizer(25, POSITIVE, 0.0) + val nnlsResult3 = nnls2.minimize(P3, q3 :* (-1.0)) + + val posResult3 = qpIters2.minimizeAndReturnState(P3, q3) + + assert(posResult3.converged) + assert(norm(posResult3.x - nnlsResult3, 2) <= 1e-2) + + val qpObj = QuadraticMinimizer.computeObjective(P1, q1, qpItersResult.x) + val nnlsObj = NNLS.computeObjectiveValue(P1, q1 :* (-1.0), nnlsResult) + + assert(abs(qpObj - nnlsObj) < 1e-4) + assert(qpItersResult.converged) + assert(norm(qpItersResult.x - nnlsResult, 2) < 1e-4) + } + + test("Quadratic Minimization consistency with multiple updates to proximal operator") { + //Ill conditioned problems + //ALS: Diagnosing userOrProduct 34 lambdaL2 6.500000000000006E-4 lambdaL1 0.06435 + val data = Array(0.172933, 0.121672, 0.083271, 0.081677, 0.053990, 0.060374, 0.122440, 0.057874, 0.202366, 0.058154, 0.221297, 0.098373, 0.027677, 0.208701, 0.174751, 0.133861, 0.143158, 0.221107, 0.180647, 0.059440, 0.121672, 0.107064, 0.044376, 0.055041, 0.020705, 0.040884, 0.080482, 0.049619, 0.105660, 0.054398, 0.195126, 0.069424, 0.024926, 0.116618, 0.075244, 0.115016, 0.125965, 0.148266, 0.150126, 0.043844, 0.083271, 0.044376, 0.057854, 0.042197, 0.040014, 0.023171, 0.074746, 0.018888, 0.114484, 0.012764, 0.071563, 0.057906, 0.007638, 0.137569, 0.114753, 0.043443, 0.046860, 0.106226, 0.069941, 0.033199, 0.081677, 0.055041, 0.042197, 0.041480, 0.028744, 0.029312, 0.059437, 0.026284, 0.103280, 0.025698, 0.099548, 0.046935, 0.012378, 0.104832, 0.091798, 0.060753, 0.064436, 0.107269, 0.082718, 0.028049, 0.053990, 0.020705, 0.040014, 0.028744, 0.036006, 0.022218, 0.043569, 0.010584, 0.102747, 0.006677, 0.034772, 0.030206, 0.003981, 0.094862, 0.103568, 0.024249, 0.022693, 0.079784, 0.036521, 0.016478, 0.060374, 0.040884, 0.023171, 0.029312, 0.022218, 0.041842, 0.021271, 0.025349, 0.110247, 0.030934, 0.087805, 0.012667, 0.012995, 0.053306, 0.089922, 0.058710, 0.056007, 0.099385, 0.063017, 0.008061, 0.122440, 0.080482, 0.074746, 0.059437, 0.043569, 0.021271, 0.120124, 0.031163, 0.115960, 0.022294, 0.127774, 0.098052, 0.013227, 0.190554, 0.114633, 0.071859, 0.083771, 0.136461, 0.119578, 0.057799, 0.057874, 0.049619, 0.018888, 0.026284, 0.010584, 0.025349, 0.031163, 0.027251, 0.062063, 0.029397, 0.096479, 0.025963, 0.012966, 0.048467, 0.044052, 0.058668, 0.062030, 0.077193, 0.071494, 0.016748, 0.202366, 0.105660, 0.114484, 0.103280, 0.102747, 0.110247, 0.115960, 0.062063, 0.385655, 0.064201, 0.210133, 0.076002, 0.029096, 0.267930, 0.351138, 0.143235, 0.134931, 0.318630, 0.172262, 0.043164, 0.058154, 0.054398, 0.012764, 0.025698, 0.006677, 0.030934, 0.022294, 0.029397, 0.064201, 0.038538, 0.110689, 0.018742, 0.015550, 0.033089, 0.039631, 0.068239, 0.070888, 0.081952, 0.077663, 0.013057, 0.221297, 0.195126, 0.071563, 0.099548, 0.034772, 0.087805, 0.127774, 0.096479, 0.210133, 0.110689, 0.375374, 0.109388, 0.049333, 0.187174, 0.143602, 0.223740, 0.240412, 0.283103, 0.278892, 0.070417, 0.098373, 0.069424, 0.057906, 0.046935, 0.030206, 0.012667, 0.098052, 0.025963, 0.076002, 0.018742, 0.109388, 0.085273, 0.011164, 0.149983, 0.074667, 0.059915, 0.071766, 0.103003, 0.101268, 0.049539, 0.027677, 0.024926, 0.007638, 0.012378, 0.003981, 0.012995, 0.013227, 0.012966, 0.029096, 0.015550, 0.049333, 0.011164, 0.008698, 0.019803, 0.019141, 0.030083, 0.031669, 0.037444, 0.035675, 0.007410, 0.208701, 0.116618, 0.137569, 0.104832, 0.094862, 0.053306, 0.190554, 0.048467, 0.267930, 0.033089, 0.187174, 0.149983, 0.019803, 0.342885, 0.268144, 0.111535, 0.122612, 0.258983, 0.181204, 0.086551, 0.174751, 0.075244, 0.114753, 0.091798, 0.103568, 0.089922, 0.114633, 0.044052, 0.351138, 0.039631, 0.143602, 0.074667, 0.019141, 0.268144, 0.338802, 0.101363, 0.092574, 0.275155, 0.129493, 0.040859, 0.133861, 0.115016, 0.043443, 0.060753, 0.024249, 0.058710, 0.071859, 0.058668, 0.143235, 0.068239, 0.223740, 0.059915, 0.030083, 0.111535, 0.101363, 0.137992, 0.143845, 0.178542, 0.165663, 0.038686, 0.143158, 0.125965, 0.046860, 0.064436, 0.022693, 0.056007, 0.083771, 0.062030, 0.134931, 0.070888, 0.240412, 0.071766, 0.031669, 0.122612, 0.092574, 0.143845, 0.156765, 0.182358, 0.180036, 0.046104, 0.221107, 0.148266, 0.106226, 0.107269, 0.079784, 0.099385, 0.136461, 0.077193, 0.318630, 0.081952, 0.283103, 0.103003, 0.037444, 0.258983, 0.275155, 0.178542, 0.182358, 0.314775, 0.224864, 0.061995, 0.180647, 0.150126, 0.069941, 0.082718, 0.036521, 0.063017, 0.119578, 0.071494, 0.172262, 0.077663, 0.278892, 0.101268, 0.035675, 0.181204, 0.129493, 0.165663, 0.180036, 0.224864, 0.218079, 0.063453, 0.059440, 0.043844, 0.033199, 0.028049, 0.016478, 0.008061, 0.057799, 0.016748, 0.043164, 0.013057, 0.070417, 0.049539, 0.007410, 0.086551, 0.040859, 0.038686, 0.046104, 0.061995, 0.063453, 0.031610) + val Psparse = new DenseMatrix(20, 20, data) + val qsparse = DenseVector(-2.631990, -1.894166, -1.225159, -1.254382, -0.816634, -1.018898, -1.763319, -0.929908, -3.243997, -0.969543, -3.519254, -1.402682, -0.451124, -3.055922, -2.747532, -2.151673, -2.272190, -3.494880, -2.813748, -0.852818) + + val proxL1 = ProximalL1() + val qpSparse = new QuadraticMinimizer(20, proxL1.setLambda(0.06435)) + val qpSparseGold = qpSparse.minimize(Psparse, qsparse) + + + proxL1.setLambda(0.0) + + val P1 = new DenseMatrix(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)) + val q1 = DenseVector(-1880.240029, -1920.949941, -2030.476172, -1956.642164, -2021.502985, -2090.503157, -1965.934820, -2072.855628, -2098.075034, -2035.059185, -1999.005923, -2121.515181, -1944.759586, -2035.397706, -1939.872057, -1888.635008, -1986.031605, -1973.738457, -2024.468051, -2003.765736) + + qpSparse.minimize(P1, q1) + proxL1.setLambda(0.06435) + + val qpSparseResult = qpSparse.minimize(Psparse, qsparse) + assert(norm(qpSparseResult - qpSparseGold, 2) < 1E-4) + } + + test("minimize API should generate identical answers") { + val n = 5 + 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) + atb *= -1.0 + + val qpSolverPos = QuadraticMinimizer(n, POSITIVE, 0.0) + val posResult = qpSolverPos.minimize(ata, atb) + + val qpSolverPosTest = QuadraticMinimizer(n, POSITIVE, 0.0) + ata.iterator.foreach{ + case ((row, col), entry) => qpSolverPosTest.updateGram(row, col, entry) + } + val posResultTest = qpSolverPosTest.minimize(atb) + assert(norm(posResult - posResultTest, inf) < 1E-6) + } +} \ No newline at end of file diff --git a/project/plugins.sbt b/project/plugins.sbt index d512d6f23..cd79a74a1 100644 --- a/project/plugins.sbt +++ b/project/plugins.sbt @@ -7,3 +7,6 @@ addSbtPlugin("com.eed3si9n" % "sbt-assembly" % "0.10.0") addSbtPlugin("com.github.mpeltonen" % "sbt-idea" % "1.6.0") addSbtPlugin("de.johoop" % "jacoco4sbt" % "2.1.6") + +addSbtPlugin("com.typesafe.sbteclipse" % "sbteclipse-plugin" % "2.2.0") +