From 574c9abbf2b9311e8678ac4c857241c6328d89b6 Mon Sep 17 00:00:00 2001 From: Thomas Neidhart Date: Sat, 28 Jul 2012 16:37:40 +0000 Subject: [PATCH] [MATH-828] Fixed numerical instabilities in SimplexSolver leading to unbounded solutions. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1366707 13f79535-47bb-0310-9956-ffa450edef68 --- .../optimization/linear/SimplexSolver.java | 54 +++++++++++++----- .../optimization/linear/SimplexTableau.java | 15 ++++- .../linear/SimplexSolverTest.java | 55 +++++++++++++++++++ 3 files changed, 107 insertions(+), 17 deletions(-) diff --git a/src/main/java/org/apache/commons/math3/optimization/linear/SimplexSolver.java b/src/main/java/org/apache/commons/math3/optimization/linear/SimplexSolver.java index c3d90890f..c2fa14dad 100644 --- a/src/main/java/org/apache/commons/math3/optimization/linear/SimplexSolver.java +++ b/src/main/java/org/apache/commons/math3/optimization/linear/SimplexSolver.java @@ -71,7 +71,9 @@ public class SimplexSolver extends AbstractLinearOptimizer { Integer minPos = null; for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) { final double entry = tableau.getEntry(0, i); - if (Precision.compareTo(entry, minValue, maxUlps) < 0) { + // check if the entry is strictly smaller than the current minimum + // do not use a ulp/epsilon check + if (entry < minValue) { minValue = entry; minPos = i; } @@ -95,7 +97,9 @@ public class SimplexSolver extends AbstractLinearOptimizer { if (Precision.compareTo(entry, 0d, maxUlps) > 0) { final double ratio = rhs / entry; - final int cmp = Precision.compareTo(ratio, minRatio, maxUlps); + // check if the entry is strictly equal to the current min ratio + // do not use a ulp/epsilon check + final int cmp = Double.compare(ratio, minRatio); if (cmp == 0) { minRatioPositions.add(i); } else if (cmp < 0) { @@ -107,20 +111,40 @@ public class SimplexSolver extends AbstractLinearOptimizer { } if (minRatioPositions.size() == 0) { - return null; + return null; } else if (minRatioPositions.size() > 1) { - // there's a degeneracy as indicated by a tie in the minimum ratio test - // check if there's an artificial variable that can be forced out of the basis - for (Integer row : minRatioPositions) { - for (int i = 0; i < tableau.getNumArtificialVariables(); i++) { - int column = i + tableau.getArtificialVariableOffset(); - final double entry = tableau.getEntry(row, column); - if (Precision.equals(entry, 1d, maxUlps) && - row.equals(tableau.getBasicRow(column))) { - return row; - } + // there's a degeneracy as indicated by a tie in the minimum ratio test + + // 1. check if there's an artificial variable that can be forced out of the basis + for (Integer row : minRatioPositions) { + for (int i = 0; i < tableau.getNumArtificialVariables(); i++) { + int column = i + tableau.getArtificialVariableOffset(); + final double entry = tableau.getEntry(row, column); + if (Precision.equals(entry, 1d, maxUlps) && row.equals(tableau.getBasicRow(column))) { + return row; + } + } } - } + + // 2. apply Bland's rule to prevent cycling: + // take the row for which the corresponding basic variable has the smallest index + // + // see http://www.stanford.edu/class/msande310/blandrule.pdf + // see http://en.wikipedia.org/wiki/Bland%27s_rule (not equivalent to the above paper) + Integer minRow = null; + int minIndex = tableau.getWidth(); + for (Integer row : minRatioPositions) { + for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1 && minRow != row; i++) { + if (row == tableau.getBasicRow(i)) { + if (i < minIndex) { + minIndex = i; + minRow = row; + } + } + } + } + + return minRow; } return minRatioPositions.get(0); } @@ -149,7 +173,7 @@ public class SimplexSolver extends AbstractLinearOptimizer { // set the rest of the pivot column to 0 for (int i = 0; i < tableau.getHeight(); i++) { if (i != pivotRow) { - double multiplier = tableau.getEntry(i, pivotCol); + final double multiplier = tableau.getEntry(i, pivotCol); tableau.subtractRow(i, pivotRow, multiplier); } } diff --git a/src/main/java/org/apache/commons/math3/optimization/linear/SimplexTableau.java b/src/main/java/org/apache/commons/math3/optimization/linear/SimplexTableau.java index 327b2ae65..e07cfc5e5 100644 --- a/src/main/java/org/apache/commons/math3/optimization/linear/SimplexTableau.java +++ b/src/main/java/org/apache/commons/math3/optimization/linear/SimplexTableau.java @@ -27,6 +27,8 @@ import java.util.HashSet; import java.util.List; import java.util.Set; +import odk.lang.FastMath; + import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealMatrix; @@ -68,6 +70,9 @@ class SimplexTableau implements Serializable { /** Default amount of error to accept in floating point comparisons (as ulps). */ private static final int DEFAULT_ULPS = 10; + /** The cut-off threshold to zero-out entries. */ + private static final double CUTOFF_THRESHOLD = 1e-12; + /** Serializable version identifier. */ private static final long serialVersionUID = -1369660067587938365L; @@ -453,8 +458,14 @@ class SimplexTableau implements Serializable { */ protected void subtractRow(final int minuendRow, final int subtrahendRow, final double multiple) { - tableau.setRowVector(minuendRow, tableau.getRowVector(minuendRow) - .subtract(tableau.getRowVector(subtrahendRow).mapMultiply(multiple))); + for (int i = 0; i < getWidth(); i++) { + double result = tableau.getEntry(minuendRow, i) - tableau.getEntry(subtrahendRow, i) * multiple; + // cut-off values smaller than the CUTOFF_THRESHOLD, otherwise may lead to numerical instabilities + if (FastMath.abs(result) < CUTOFF_THRESHOLD) { + result = 0.0; + } + tableau.setEntry(minuendRow, i, result); + } } /** diff --git a/src/test/java/org/apache/commons/math3/optimization/linear/SimplexSolverTest.java b/src/test/java/org/apache/commons/math3/optimization/linear/SimplexSolverTest.java index ae4d706b3..b4bf643ad 100644 --- a/src/test/java/org/apache/commons/math3/optimization/linear/SimplexSolverTest.java +++ b/src/test/java/org/apache/commons/math3/optimization/linear/SimplexSolverTest.java @@ -21,6 +21,7 @@ import org.junit.Assert; import java.util.ArrayList; import java.util.Collection; +import java.util.List; import org.apache.commons.math3.optimization.GoalType; import org.apache.commons.math3.optimization.PointValuePair; @@ -29,6 +30,27 @@ import org.junit.Test; public class SimplexSolverTest { + @Test + public void testMath828() { + LinearObjectiveFunction f = new LinearObjectiveFunction( + new double[] { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, 0.0); + + ArrayList constraints = new ArrayList(); + + constraints.add(new LinearConstraint(new double[] {0.0, 39.0, 23.0, 96.0, 15.0, 48.0, 9.0, 21.0, 48.0, 36.0, 76.0, 19.0, 88.0, 17.0, 16.0, 36.0,}, Relationship.GEQ, 15.0)); + constraints.add(new LinearConstraint(new double[] {0.0, 59.0, 93.0, 12.0, 29.0, 78.0, 73.0, 87.0, 32.0, 70.0, 68.0, 24.0, 11.0, 26.0, 65.0, 25.0,}, Relationship.GEQ, 29.0)); + constraints.add(new LinearConstraint(new double[] {0.0, 74.0, 5.0, 82.0, 6.0, 97.0, 55.0, 44.0, 52.0, 54.0, 5.0, 93.0, 91.0, 8.0, 20.0, 97.0,}, Relationship.GEQ, 6.0)); + constraints.add(new LinearConstraint(new double[] {8.0, -3.0, -28.0, -72.0, -8.0, -31.0, -31.0, -74.0, -47.0, -59.0, -24.0, -57.0, -56.0, -16.0, -92.0, -59.0,}, Relationship.GEQ, 0.0)); + constraints.add(new LinearConstraint(new double[] {25.0, -7.0, -99.0, -78.0, -25.0, -14.0, -16.0, -89.0, -39.0, -56.0, -53.0, -9.0, -18.0, -26.0, -11.0, -61.0,}, Relationship.GEQ, 0.0)); + constraints.add(new LinearConstraint(new double[] {33.0, -95.0, -15.0, -4.0, -33.0, -3.0, -20.0, -96.0, -27.0, -13.0, -80.0, -24.0, -3.0, -13.0, -57.0, -76.0,}, Relationship.GEQ, 0.0)); + constraints.add(new LinearConstraint(new double[] {7.0, -95.0, -39.0, -93.0, -7.0, -94.0, -94.0, -62.0, -76.0, -26.0, -53.0, -57.0, -31.0, -76.0, -53.0, -52.0,}, Relationship.GEQ, 0.0)); + + double epsilon = 1e-6; + PointValuePair solution = new SimplexSolver().optimize(f, constraints, GoalType.MINIMIZE, true); + Assert.assertEquals(1.0d, solution.getValue(), epsilon); + Assert.assertTrue(validSolution(solution, constraints, epsilon)); + } + @Test public void testMath781() { LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 2, 6, 7 }, 0); @@ -560,4 +582,37 @@ public class SimplexSolverTest { return new LinearConstraint(lhs, relationship, rhs); } + private static boolean validSolution(PointValuePair solution, List constraints, double epsilon) { + double[] vals = solution.getPoint(); + for (LinearConstraint c : constraints) { + double[] coeffs = c.getCoefficients().toArray(); + double result = 0.0d; + for (int i = 0; i < vals.length; i++) { + result += vals[i] * coeffs[i]; + } + + switch (c.getRelationship()) { + case EQ: + if (!Precision.equals(result, c.getValue(), epsilon)) { + return false; + } + break; + + case GEQ: + if (Precision.compareTo(result, c.getValue(), epsilon) < 0) { + return false; + } + break; + + case LEQ: + if (Precision.compareTo(result, c.getValue(), epsilon) > 0) { + return false; + } + break; + } + } + + return true; + } + }