From 881a4ee8db1e1826b749579d2284ace982c953ca Mon Sep 17 00:00:00 2001 From: Thomas Neidhart Date: Wed, 18 Dec 2013 17:41:26 +0000 Subject: [PATCH] [MATH-1079] Further improvements to SimplexSolver. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1552046 13f79535-47bb-0310-9956-ffa450edef68 --- .../math3/optim/linear/SimplexSolver.java | 14 ++++--- .../math3/optim/linear/SimplexTableau.java | 39 +++---------------- .../math3/optim/linear/SimplexSolverTest.java | 7 ++-- 3 files changed, 18 insertions(+), 42 deletions(-) diff --git a/src/main/java/org/apache/commons/math3/optim/linear/SimplexSolver.java b/src/main/java/org/apache/commons/math3/optim/linear/SimplexSolver.java index d0569f735..4af68ee08 100644 --- a/src/main/java/org/apache/commons/math3/optim/linear/SimplexSolver.java +++ b/src/main/java/org/apache/commons/math3/optim/linear/SimplexSolver.java @@ -22,6 +22,7 @@ import java.util.List; import org.apache.commons.math3.exception.TooManyIterationsException; import org.apache.commons.math3.optim.OptimizationData; import org.apache.commons.math3.optim.PointValuePair; +import org.apache.commons.math3.util.FastMath; import org.apache.commons.math3.util.Precision; /** @@ -49,7 +50,7 @@ import org.apache.commons.math3.util.Precision; * *

* The cut-off value has been introduced to zero out very small numbers in the Simplex tableau, @@ -66,7 +67,7 @@ public class SimplexSolver extends LinearOptimizer { static final int DEFAULT_ULPS = 10; /** Default cut-off value. */ - static final double DEFAULT_CUT_OFF = 1e-12; + static final double DEFAULT_CUT_OFF = 1e-10; /** Default amount of error to accept for algorithm convergence. */ private static final double DEFAULT_EPSILON = 1.0e-6; @@ -249,7 +250,11 @@ public class SimplexSolver extends LinearOptimizer { final double rhs = tableau.getEntry(i, tableau.getWidth() - 1); final double entry = tableau.getEntry(i, col); - if (Precision.compareTo(entry, 0d, maxUlps) > 0) { + // zero-out tableau entries that are too close to zero to avoid + // numerical instabilities as these entries might be used as divisor + if (FastMath.abs(entry) < cutOff) { + tableau.setEntry(i, col, 0); + } else if (Precision.compareTo(entry, 0d, maxUlps) > 0) { final double ratio = rhs / entry; // check if the entry is strictly equal to the current min ratio // do not use a ulp/epsilon check @@ -371,8 +376,7 @@ public class SimplexSolver extends LinearOptimizer { getGoalType(), isRestrictedToNonNegative(), epsilon, - maxUlps, - cutOff); + maxUlps); solvePhase1(tableau); tableau.dropPhase1Objective(); diff --git a/src/main/java/org/apache/commons/math3/optim/linear/SimplexTableau.java b/src/main/java/org/apache/commons/math3/optim/linear/SimplexTableau.java index 372dd35d1..04b19a24a 100644 --- a/src/main/java/org/apache/commons/math3/optim/linear/SimplexTableau.java +++ b/src/main/java/org/apache/commons/math3/optim/linear/SimplexTableau.java @@ -33,7 +33,6 @@ import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealVector; import org.apache.commons.math3.optim.nonlinear.scalar.GoalType; import org.apache.commons.math3.optim.PointValuePair; -import org.apache.commons.math3.util.FastMath; import org.apache.commons.math3.util.Precision; /** @@ -99,9 +98,6 @@ class SimplexTableau implements Serializable { /** Amount of error to accept in floating point comparisons. */ private final int maxUlps; - /** Cut-off value for entries in the tableau. */ - private final double cutOff; - /** Maps basic variables to row they are basic in. */ private int[] basicVariables; @@ -123,8 +119,7 @@ class SimplexTableau implements Serializable { final GoalType goalType, final boolean restrictToNonNegative, final double epsilon) { - this(f, constraints, goalType, restrictToNonNegative, epsilon, - SimplexSolver.DEFAULT_ULPS, SimplexSolver.DEFAULT_CUT_OFF); + this(f, constraints, goalType, restrictToNonNegative, epsilon, SimplexSolver.DEFAULT_ULPS); } /** @@ -142,32 +137,11 @@ class SimplexTableau implements Serializable { final boolean restrictToNonNegative, final double epsilon, final int maxUlps) { - this(f, constraints, goalType, restrictToNonNegative, epsilon, maxUlps, SimplexSolver.DEFAULT_CUT_OFF); - } - - /** - * Build a tableau for a linear problem. - * @param f linear objective function - * @param constraints linear constraints - * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE} or {@link GoalType#MINIMIZE} - * @param restrictToNonNegative whether to restrict the variables to non-negative values - * @param epsilon amount of error to accept when checking for optimality - * @param maxUlps amount of error to accept in floating point comparisons - * @param cutOff the cut-off value for tableau entries - */ - SimplexTableau(final LinearObjectiveFunction f, - final Collection constraints, - final GoalType goalType, - final boolean restrictToNonNegative, - final double epsilon, - final int maxUlps, - final double cutOff) { this.f = f; this.constraints = normalizeConstraints(constraints); this.restrictToNonNegative = restrictToNonNegative; this.epsilon = epsilon; this.maxUlps = maxUlps; - this.cutOff = cutOff; this.numDecisionVariables = f.getCoefficients().getDimension() + (restrictToNonNegative ? 0 : 1); this.numSlackVariables = getConstraintTypeCounts(Relationship.LEQ) + getConstraintTypeCounts(Relationship.GEQ); @@ -516,7 +490,9 @@ class SimplexTableau implements Serializable { for (int i = 0; i < getHeight(); i++) { if (i != pivotRow) { final double multiplier = getEntry(i, pivotCol); - subtractRow(i, pivotRow, multiplier); + if (multiplier != 0.0) { + subtractRow(i, pivotRow, multiplier); + } } } @@ -557,12 +533,7 @@ class SimplexTableau implements Serializable { final double[] minuendRow = getRow(minuendRowIndex); final double[] subtrahendRow = getRow(subtrahendRowIndex); for (int i = 0; i < getWidth(); i++) { - double result = minuendRow[i] - subtrahendRow[i] * multiplier; - // cut-off values smaller than the cut-off threshold, otherwise may lead to numerical instabilities - if (result != 0.0 && FastMath.abs(result) < cutOff) { - result = 0.0; - } - minuendRow[i] = result; + minuendRow[i] -= subtrahendRow[i] * multiplier; } } diff --git a/src/test/java/org/apache/commons/math3/optim/linear/SimplexSolverTest.java b/src/test/java/org/apache/commons/math3/optim/linear/SimplexSolverTest.java index 1d1090948..2c4aa8982 100644 --- a/src/test/java/org/apache/commons/math3/optim/linear/SimplexSolverTest.java +++ b/src/test/java/org/apache/commons/math3/optim/linear/SimplexSolverTest.java @@ -750,7 +750,7 @@ public class SimplexSolverTest { @Test public void testSolutionCallback() { // re-use the problem from testcase for MATH-930 - // it normally requires 113 iterations + // it normally requires 144 iterations final List constraints = createMath930Constraints(); double[] objFunctionCoeff = new double[33]; @@ -777,9 +777,10 @@ public class SimplexSolverTest { // 2. iteration limit allows to reach phase 2, but too low to find an optimal solution try { // we need to use a DeterministicLinearConstraintSet to always get the same behavior - solver.optimize(new MaxIter(112), f, new LinearConstraintSet(constraints), + solver.optimize(new MaxIter(140), f, new LinearConstraintSet(constraints), GoalType.MINIMIZE, new NonNegativeConstraint(true), callback, - PivotSelectionRule.BLAND); + PivotSelectionRule.DANTZIG); + System.out.println(solver.getIterations()); Assert.fail("expected TooManyIterationsException"); } catch (TooManyIterationsException ex) { // expected