[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

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;
  * <ul>
  *   <li>Algorithm convergence: 1e-6</li>
  *   <li>Floating-point comparisons: 10 ulp</li>
- *   <li>Cut-Off value: 1e-12</li>
+ *   <li>Cut-Off value: 1e-10</li>
  * </ul>
  * <p>
  * 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,
+                               maxUlps);
-                               cutOff);
 
         solvePhase1(tableau);
         tableau.dropPhase1Objective();

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,
+        this(f, constraints, goalType, restrictToNonNegative, epsilon, SimplexSolver.DEFAULT_ULPS);
-             SimplexSolver.DEFAULT_ULPS, SimplexSolver.DEFAULT_CUT_OFF);
     }
 
     /**
@@ -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<LinearConstraint> 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;
+            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;
         }
     }
 

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<LinearConstraint> 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