Fixed two errors in simplex solver when entries are close together or
when variables are not restricted to non-negative. Jira: MATH-434 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1090656 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
a3c552e324
commit
133cbc2dbf
3
pom.xml
3
pom.xml
|
|
@ -186,6 +186,9 @@
|
|||
<contributor>
|
||||
<name>J. Lewis Muir</name>
|
||||
</contributor>
|
||||
<contributor>
|
||||
<name>Thomas Neidhart</name>
|
||||
</contributor>
|
||||
<contributor>
|
||||
<name>Fredrik Norin</name>
|
||||
</contributor>
|
||||
|
|
|
|||
|
|
@ -22,6 +22,7 @@ import java.util.List;
|
|||
|
||||
import org.apache.commons.math.optimization.OptimizationException;
|
||||
import org.apache.commons.math.optimization.RealPointValuePair;
|
||||
import org.apache.commons.math.util.FastMath;
|
||||
import org.apache.commons.math.util.MathUtils;
|
||||
|
||||
|
||||
|
|
@ -31,26 +32,34 @@ import org.apache.commons.math.util.MathUtils;
|
|||
* @since 2.0
|
||||
*/
|
||||
public class SimplexSolver extends AbstractLinearOptimizer {
|
||||
|
||||
/** Default amount of error to accept in floating point comparisons. */
|
||||
|
||||
/** Default amount of error to accept for algorithm convergence. */
|
||||
private static final double DEFAULT_EPSILON = 1.0e-6;
|
||||
|
||||
/** Amount of error to accept in floating point comparisons. */
|
||||
|
||||
/** Amount of error to accept for algorithm convergence. */
|
||||
protected final double epsilon;
|
||||
|
||||
/** Default amount of error to accept in floating point comparisons (as ulps). */
|
||||
private static final int DEFAULT_ULPS = 10;
|
||||
|
||||
/** Amount of error to accept in floating point comparisons (as ulps). */
|
||||
protected final int maxUlps;
|
||||
|
||||
/**
|
||||
* Build a simplex solver with default settings.
|
||||
*/
|
||||
public SimplexSolver() {
|
||||
this(DEFAULT_EPSILON);
|
||||
this(DEFAULT_EPSILON, DEFAULT_ULPS);
|
||||
}
|
||||
|
||||
/**
|
||||
* Build a simplex solver with a specified accepted amount of error
|
||||
* @param epsilon the amount of error to accept in floating point comparisons
|
||||
* @param epsilon the amount of error to accept for algorithm convergence
|
||||
* @param maxUlps amount of error to accept in floating point comparisons
|
||||
*/
|
||||
public SimplexSolver(final double epsilon) {
|
||||
public SimplexSolver(final double epsilon, final int maxUlps) {
|
||||
this.epsilon = epsilon;
|
||||
this.maxUlps = maxUlps;
|
||||
}
|
||||
|
||||
/**
|
||||
|
|
@ -62,8 +71,9 @@ public class SimplexSolver extends AbstractLinearOptimizer {
|
|||
double minValue = 0;
|
||||
Integer minPos = null;
|
||||
for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) {
|
||||
if (MathUtils.compareTo(tableau.getEntry(0, i), minValue, epsilon) < 0) {
|
||||
minValue = tableau.getEntry(0, i);
|
||||
final double entry = tableau.getEntry(0, i);
|
||||
if (MathUtils.compareTo(entry, minValue, getEpsilon(entry)) < 0) {
|
||||
minValue = entry;
|
||||
minPos = i;
|
||||
}
|
||||
}
|
||||
|
|
@ -83,11 +93,13 @@ public class SimplexSolver extends AbstractLinearOptimizer {
|
|||
for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) {
|
||||
final double rhs = tableau.getEntry(i, tableau.getWidth() - 1);
|
||||
final double entry = tableau.getEntry(i, col);
|
||||
if (MathUtils.compareTo(entry, 0, epsilon) > 0) {
|
||||
|
||||
if (MathUtils.compareTo(entry, 0d, getEpsilon(entry)) > 0) {
|
||||
final double ratio = rhs / entry;
|
||||
if (MathUtils.equals(ratio, minRatio, epsilon)) {
|
||||
final int cmp = MathUtils.compareTo(ratio, minRatio, getEpsilon(ratio));
|
||||
if (cmp == 0) {
|
||||
minRatioPositions.add(i);
|
||||
} else if (ratio < minRatio) {
|
||||
} else if (cmp < 0) {
|
||||
minRatio = ratio;
|
||||
minRatioPositions = new ArrayList<Integer>();
|
||||
minRatioPositions.add(i);
|
||||
|
|
@ -103,7 +115,8 @@ public class SimplexSolver extends AbstractLinearOptimizer {
|
|||
for (Integer row : minRatioPositions) {
|
||||
for (int i = 0; i < tableau.getNumArtificialVariables(); i++) {
|
||||
int column = i + tableau.getArtificialVariableOffset();
|
||||
if (MathUtils.equals(tableau.getEntry(row, column), 1, epsilon) &&
|
||||
final double entry = tableau.getEntry(row, column);
|
||||
if (MathUtils.equals(entry, 1d, getEpsilon(entry)) &&
|
||||
row.equals(tableau.getBasicRow(column))) {
|
||||
return row;
|
||||
}
|
||||
|
|
@ -162,7 +175,7 @@ public class SimplexSolver extends AbstractLinearOptimizer {
|
|||
}
|
||||
|
||||
// if W is not zero then we have no feasible solution
|
||||
if (!MathUtils.equals(tableau.getEntry(0, tableau.getRhsOffset()), 0, epsilon)) {
|
||||
if (!MathUtils.equals(tableau.getEntry(0, tableau.getRhsOffset()), 0d, epsilon)) {
|
||||
throw new NoFeasibleSolutionException();
|
||||
}
|
||||
}
|
||||
|
|
@ -171,7 +184,8 @@ public class SimplexSolver extends AbstractLinearOptimizer {
|
|||
@Override
|
||||
public RealPointValuePair doOptimize() throws OptimizationException {
|
||||
final SimplexTableau tableau =
|
||||
new SimplexTableau(function, linearConstraints, goal, nonNegative, epsilon);
|
||||
new SimplexTableau(function, linearConstraints, goal, nonNegative,
|
||||
epsilon, maxUlps);
|
||||
|
||||
solvePhase1(tableau);
|
||||
tableau.dropPhase1Objective();
|
||||
|
|
@ -182,4 +196,12 @@ public class SimplexSolver extends AbstractLinearOptimizer {
|
|||
return tableau.getSolution();
|
||||
}
|
||||
|
||||
/**
|
||||
* Get an epsilon that is adjusted to the magnitude of the given value.
|
||||
* @param value the value for which to get the epsilon
|
||||
* @return magnitude-adjusted epsilon using {@link FastMath.ulp}
|
||||
*/
|
||||
private double getEpsilon(double value) {
|
||||
return FastMath.ulp(value) * (double) maxUlps;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -33,6 +33,7 @@ import org.apache.commons.math.linear.RealMatrix;
|
|||
import org.apache.commons.math.linear.RealVector;
|
||||
import org.apache.commons.math.optimization.GoalType;
|
||||
import org.apache.commons.math.optimization.RealPointValuePair;
|
||||
import org.apache.commons.math.util.FastMath;
|
||||
import org.apache.commons.math.util.MathUtils;
|
||||
|
||||
/**
|
||||
|
|
@ -65,6 +66,9 @@ class SimplexTableau implements Serializable {
|
|||
/** Column label for negative vars. */
|
||||
private static final String NEGATIVE_VAR_COLUMN_LABEL = "x-";
|
||||
|
||||
/** Default amount of error to accept in floating point comparisons (as ulps). */
|
||||
private static final int DEFAULT_ULPS = 10;
|
||||
|
||||
/** Serializable version identifier. */
|
||||
private static final long serialVersionUID = -1369660067587938365L;
|
||||
|
||||
|
|
@ -92,9 +96,12 @@ class SimplexTableau implements Serializable {
|
|||
/** Number of artificial variables. */
|
||||
private int numArtificialVariables;
|
||||
|
||||
/** Amount of error to accept in floating point comparisons. */
|
||||
/** Amount of error to accept when checking for optimality. */
|
||||
private final double epsilon;
|
||||
|
||||
/** Amount of error to accept in floating point comparisons. */
|
||||
private final int maxUlps;
|
||||
|
||||
/**
|
||||
* Build a tableau for a linear problem.
|
||||
* @param f linear objective function
|
||||
|
|
@ -102,16 +109,35 @@ class SimplexTableau implements Serializable {
|
|||
* @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 in floating point comparisons
|
||||
* @param epsilon amount of error to accept when checking for optimality
|
||||
*/
|
||||
SimplexTableau(final LinearObjectiveFunction f,
|
||||
final Collection<LinearConstraint> constraints,
|
||||
final GoalType goalType, final boolean restrictToNonNegative,
|
||||
final double epsilon) {
|
||||
this(f, constraints, goalType, restrictToNonNegative, epsilon, DEFAULT_ULPS);
|
||||
}
|
||||
|
||||
/**
|
||||
* 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
|
||||
*/
|
||||
SimplexTableau(final LinearObjectiveFunction f,
|
||||
final Collection<LinearConstraint> constraints,
|
||||
final GoalType goalType, final boolean restrictToNonNegative,
|
||||
final double epsilon,
|
||||
final int maxUlps) {
|
||||
this.f = f;
|
||||
this.constraints = normalizeConstraints(constraints);
|
||||
this.restrictToNonNegative = restrictToNonNegative;
|
||||
this.epsilon = epsilon;
|
||||
this.maxUlps = maxUlps;
|
||||
this.numDecisionVariables = f.getCoefficients().getDimension() +
|
||||
(restrictToNonNegative ? 0 : 1);
|
||||
this.numSlackVariables = getConstraintTypeCounts(Relationship.LEQ) +
|
||||
|
|
@ -172,7 +198,7 @@ class SimplexTableau implements Serializable {
|
|||
|
||||
if (!restrictToNonNegative) {
|
||||
matrix.setEntry(zIndex, getSlackVariableOffset() - 1,
|
||||
getInvertedCoeffiecientSum(objectiveCoefficients));
|
||||
getInvertedCoefficientSum(objectiveCoefficients));
|
||||
}
|
||||
|
||||
// initialize the constraint rows
|
||||
|
|
@ -188,7 +214,7 @@ class SimplexTableau implements Serializable {
|
|||
// x-
|
||||
if (!restrictToNonNegative) {
|
||||
matrix.setEntry(row, getSlackVariableOffset() - 1,
|
||||
getInvertedCoeffiecientSum(constraint.getCoefficients()));
|
||||
getInvertedCoefficientSum(constraint.getCoefficients()));
|
||||
}
|
||||
|
||||
// RHS
|
||||
|
|
@ -269,7 +295,7 @@ class SimplexTableau implements Serializable {
|
|||
* @param coefficients coefficients to sum
|
||||
* @return the -1 times the sum of all coefficients in the given array.
|
||||
*/
|
||||
protected static double getInvertedCoeffiecientSum(final RealVector coefficients) {
|
||||
protected static double getInvertedCoefficientSum(final RealVector coefficients) {
|
||||
double sum = 0;
|
||||
for (double coefficient : coefficients.getData()) {
|
||||
sum -= coefficient;
|
||||
|
|
@ -285,9 +311,10 @@ class SimplexTableau implements Serializable {
|
|||
protected Integer getBasicRow(final int col) {
|
||||
Integer row = null;
|
||||
for (int i = 0; i < getHeight(); i++) {
|
||||
if (MathUtils.equals(getEntry(i, col), 1.0, epsilon) && (row == null)) {
|
||||
final double entry = getEntry(i, col);
|
||||
if (MathUtils.equals(entry, 1d, getEpsilon(entry)) && (row == null)) {
|
||||
row = i;
|
||||
} else if (!MathUtils.equals(getEntry(i, col), 0.0, epsilon)) {
|
||||
} else if (!MathUtils.equals(entry, 0d, getEpsilon(entry))) {
|
||||
return null;
|
||||
}
|
||||
}
|
||||
|
|
@ -308,9 +335,10 @@ class SimplexTableau implements Serializable {
|
|||
|
||||
// positive cost non-artificial variables
|
||||
for (int i = getNumObjectiveFunctions(); i < getArtificialVariableOffset(); i++) {
|
||||
if (MathUtils.compareTo(tableau.getEntry(0, i), 0, epsilon) > 0) {
|
||||
columnsToDrop.add(i);
|
||||
}
|
||||
final double entry = tableau.getEntry(0, i);
|
||||
if (MathUtils.compareTo(entry, 0d, getEpsilon(entry)) > 0) {
|
||||
columnsToDrop.add(i);
|
||||
}
|
||||
}
|
||||
|
||||
// non-basic artificial variables
|
||||
|
|
@ -353,7 +381,8 @@ class SimplexTableau implements Serializable {
|
|||
*/
|
||||
boolean isOptimal() {
|
||||
for (int i = getNumObjectiveFunctions(); i < getWidth() - 1; i++) {
|
||||
if (MathUtils.compareTo(tableau.getEntry(0, i), 0, epsilon) < 0) {
|
||||
final double entry = tableau.getEntry(0, i);
|
||||
if (MathUtils.compareTo(entry, 0d, epsilon) < 0) {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
|
@ -382,7 +411,7 @@ class SimplexTableau implements Serializable {
|
|||
if (basicRows.contains(basicRow)) {
|
||||
// if multiple variables can take a given value
|
||||
// then we choose the first and set the rest equal to 0
|
||||
coefficients[i] = 0;
|
||||
coefficients[i] = 0 - (restrictToNonNegative ? 0 : mostNegative);
|
||||
} else {
|
||||
basicRows.add(basicRow);
|
||||
coefficients[i] =
|
||||
|
|
@ -545,6 +574,7 @@ class SimplexTableau implements Serializable {
|
|||
(numSlackVariables == rhs.numSlackVariables) &&
|
||||
(numArtificialVariables == rhs.numArtificialVariables) &&
|
||||
(epsilon == rhs.epsilon) &&
|
||||
(maxUlps == rhs.maxUlps) &&
|
||||
f.equals(rhs.f) &&
|
||||
constraints.equals(rhs.constraints) &&
|
||||
tableau.equals(rhs.tableau);
|
||||
|
|
@ -560,6 +590,7 @@ class SimplexTableau implements Serializable {
|
|||
numSlackVariables ^
|
||||
numArtificialVariables ^
|
||||
Double.valueOf(epsilon).hashCode() ^
|
||||
maxUlps ^
|
||||
f.hashCode() ^
|
||||
constraints.hashCode() ^
|
||||
tableau.hashCode();
|
||||
|
|
@ -585,5 +616,13 @@ class SimplexTableau implements Serializable {
|
|||
ois.defaultReadObject();
|
||||
MatrixUtils.deserializeRealMatrix(this, "tableau", ois);
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Get an epsilon that is adjusted to the magnitude of the given value.
|
||||
* @param value the value for which to get the epsilon
|
||||
* @return magnitude-adjusted epsilon using {@link FastMath.ulp}
|
||||
*/
|
||||
private double getEpsilon(double value) {
|
||||
return FastMath.ulp(value) * (double) maxUlps;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -52,6 +52,10 @@ The <action> type attribute can be add,update,fix,remove.
|
|||
If the output is not quite correct, check for invisible trailing spaces!
|
||||
-->
|
||||
<release version="3.0" date="TBD" description="TBD">
|
||||
<action dev="luc" type="fix" issue="MATH-434" due-to="Thomas Neidhart">
|
||||
Fixed two errors in simplex solver when entries are close together or
|
||||
when variables are not restricted to non-negative.
|
||||
</action>
|
||||
<action dev="luc" type="fix" issue="MATH-547" due-to="Thomas Neidhart">
|
||||
Improved robustness of k-means++ algorithm, by tracking changes in points assignments
|
||||
to clusters.
|
||||
|
|
|
|||
|
|
@ -25,10 +25,87 @@ import java.util.Collection;
|
|||
import org.apache.commons.math.optimization.GoalType;
|
||||
import org.apache.commons.math.optimization.OptimizationException;
|
||||
import org.apache.commons.math.optimization.RealPointValuePair;
|
||||
import org.apache.commons.math.util.MathUtils;
|
||||
import org.junit.Test;
|
||||
|
||||
public class SimplexSolverTest {
|
||||
|
||||
@Test
|
||||
public void test434NegativeVariable() throws OptimizationException
|
||||
{
|
||||
LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] {0.0, 0.0, 1.0}, 0.0d);
|
||||
ArrayList<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
|
||||
constraints.add(new LinearConstraint(new double[] {1, 1, 0}, Relationship.EQ, 5));
|
||||
constraints.add(new LinearConstraint(new double[] {0, 0, 1}, Relationship.GEQ, -10));
|
||||
|
||||
double epsilon = 1e-6;
|
||||
SimplexSolver solver = new SimplexSolver();
|
||||
RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MINIMIZE, false);
|
||||
|
||||
Assert.assertEquals(5.0, solution.getPoint()[0] + solution.getPoint()[1], epsilon);
|
||||
Assert.assertEquals(-10.0, solution.getPoint()[2], epsilon);
|
||||
Assert.assertEquals(-10.0, solution.getValue(), epsilon);
|
||||
|
||||
}
|
||||
|
||||
@Test(expected = NoFeasibleSolutionException.class)
|
||||
public void test434UnfeasibleSolution() throws OptimizationException
|
||||
{
|
||||
double epsilon = 1e-6;
|
||||
|
||||
LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] {1.0, 0.0}, 0.0);
|
||||
ArrayList<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
|
||||
constraints.add(new LinearConstraint(new double[] {epsilon/2, 0.5}, Relationship.EQ, 0));
|
||||
constraints.add(new LinearConstraint(new double[] {1e-3, 0.1}, Relationship.EQ, 10));
|
||||
|
||||
SimplexSolver solver = new SimplexSolver();
|
||||
// allowing only non-negative values, no feasible solution shall be found
|
||||
solver.optimize(f, constraints, GoalType.MINIMIZE, true);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void test434PivotRowSelection() throws OptimizationException
|
||||
{
|
||||
LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] {1.0}, 0.0);
|
||||
|
||||
double epsilon = 1e-6;
|
||||
ArrayList<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
|
||||
constraints.add(new LinearConstraint(new double[] {200}, Relationship.GEQ, 1));
|
||||
constraints.add(new LinearConstraint(new double[] {100}, Relationship.GEQ, 0.499900001));
|
||||
|
||||
SimplexSolver solver = new SimplexSolver();
|
||||
RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MINIMIZE, false);
|
||||
|
||||
Assert.assertTrue(MathUtils.compareTo(solution.getPoint()[0] * 200.d, 1.d, epsilon) >= 0);
|
||||
Assert.assertEquals(0.0050, solution.getValue(), epsilon);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void test434PivotRowSelection2() throws OptimizationException
|
||||
{
|
||||
LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] {0.0d, 1.0d, 1.0d, 0.0d, 0.0d, 0.0d, 0.0d}, 0.0d);
|
||||
|
||||
ArrayList<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
|
||||
constraints.add(new LinearConstraint(new double[] {1.0d, -0.1d, 0.0d, 0.0d, 0.0d, 0.0d, 0.0d}, Relationship.EQ, -0.1d));
|
||||
constraints.add(new LinearConstraint(new double[] {1.0d, 0.0d, 0.0d, 0.0d, 0.0d, 0.0d, 0.0d}, Relationship.GEQ, -1e-18d));
|
||||
constraints.add(new LinearConstraint(new double[] {0.0d, 1.0d, 0.0d, 0.0d, 0.0d, 0.0d, 0.0d}, Relationship.GEQ, 0.0d));
|
||||
constraints.add(new LinearConstraint(new double[] {0.0d, 0.0d, 0.0d, 1.0d, 0.0d, -0.0128588d, 1e-5d}, Relationship.EQ, 0.0d));
|
||||
constraints.add(new LinearConstraint(new double[] {0.0d, 0.0d, 0.0d, 0.0d, 1.0d, 1e-5d, -0.0128586d}, Relationship.EQ, 1e-10d));
|
||||
constraints.add(new LinearConstraint(new double[] {0.0d, 0.0d, 1.0d, -1.0d, 0.0d, 0.0d, 0.0d}, Relationship.GEQ, 0.0d));
|
||||
constraints.add(new LinearConstraint(new double[] {0.0d, 0.0d, 1.0d, 1.0d, 0.0d, 0.0d, 0.0d}, Relationship.GEQ, 0.0d));
|
||||
constraints.add(new LinearConstraint(new double[] {0.0d, 0.0d, 1.0d, 0.0d, -1.0d, 0.0d, 0.0d}, Relationship.GEQ, 0.0d));
|
||||
constraints.add(new LinearConstraint(new double[] {0.0d, 0.0d, 1.0d, 0.0d, 1.0d, 0.0d, 0.0d}, Relationship.GEQ, 0.0d));
|
||||
|
||||
double epsilon = 1e-7;
|
||||
SimplexSolver simplex = new SimplexSolver();
|
||||
RealPointValuePair solution = simplex.optimize(f, constraints, GoalType.MINIMIZE, false);
|
||||
|
||||
Assert.assertTrue(MathUtils.compareTo(solution.getPoint()[0], -1e-18d, epsilon) >= 0);
|
||||
Assert.assertEquals(1.0d, solution.getPoint()[1], epsilon);
|
||||
Assert.assertEquals(0.0d, solution.getPoint()[2], epsilon);
|
||||
Assert.assertEquals(1.0d, solution.getValue(), epsilon);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testMath272() throws OptimizationException {
|
||||
LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 2, 2, 1 }, 0);
|
||||
|
|
|
|||
Loading…
Reference in New Issue