commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From t.@apache.org
Subject 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/
Date Sat, 28 Jul 2012 16:37:40 GMT
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;
+    }
+
 }



Mime
View raw message