You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by tn...@apache.org on 2012/07/28 18:37:40 UTC
svn commit: r1366707 - in /commons/proper/math/trunk/src:
main/java/org/apache/commons/math3/optimization/linear/
test/java/org/apache/commons/math3/optimization/linear/
Author: tn
Date: Sat Jul 28 16:37:40 2012
New Revision: 1366707
URL: http://svn.apache.org/viewvc?rev=1366707&view=rev
Log:
[MATH-828] Fixed numerical instabilities in SimplexSolver leading to unbounded solutions.
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/linear/SimplexSolver.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/linear/SimplexTableau.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/linear/SimplexSolverTest.java
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/linear/SimplexSolver.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/linear/SimplexSolver.java?rev=1366707&r1=1366706&r2=1366707&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/linear/SimplexSolver.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/linear/SimplexSolver.java Sat Jul 28 16:37:40 2012
@@ -71,7 +71,9 @@ public class SimplexSolver extends Abstr
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 Abstr
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 Abstr
}
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 Abstr
// 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);
}
}
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/linear/SimplexTableau.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/linear/SimplexTableau.java?rev=1366707&r1=1366706&r2=1366707&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/linear/SimplexTableau.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/linear/SimplexTableau.java Sat Jul 28 16:37:40 2012
@@ -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 Serializ
/** 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 Serializ
*/
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);
+ }
}
/**
Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/linear/SimplexSolverTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/linear/SimplexSolverTest.java?rev=1366707&r1=1366706&r2=1366707&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/linear/SimplexSolverTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/linear/SimplexSolverTest.java Sat Jul 28 16:37:40 2012
@@ -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;
@@ -30,6 +31,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 <LinearConstraint>constraints = new ArrayList<LinearConstraint>();
+
+ 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<LinearConstraint> 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;
+ }
+
}