Result change betwen ojAlgo version 53 and 54 or later

[Programmer-Magnus]

When I update ojAlgo my tests stumble on a case that gives OPTIMAL in version 53 of ojAlgo but give INFEASIBLE in v54 and later.

I tried 55.1.2 with same result as 54.0.0

The solutions are very different. With 53.0.0 I get an optimal solution and 54.0.0 say it is infeasible.

version 53 gives OPTIMAL -1.4162052484896982E24 @ { 1.406005993853654E-7, -7.895855702762938E-6, 714.3828764355712, 1.682976677491224E12, -6017.883736346659, 4605.691058549535 }

version 54 & 55 gives INFEASIBLE 0.0 @ { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }

I have made two versions since ojAlgo API change slightly between them, quadraticTest53 usable in v53 and quadraticTest54 usable with v54. In the code I left some settings, which I tried, commented away.

I am curious if this is a bug or possible to resolve by changing the settings in a way I did not yet try.

Thank you, for a great library!

@Test //Quadratic model that fail in package org.ojalgo.matrix.store, ojAlgo, version 54.0 and later, but work with 53.0.0
   public void quadraticTest53() {

      double[][] q_ = new double[][]{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 1.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 1.0, 0.0, 0.0, 0.0}, {0.0
               , 0.0, 0.0, 1.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};

      MatrixStore<Double> q = Primitive64Store.FACTORY.makeWrapper(RawStore.wrap(q_));

      double[][] l_ = new double[][]{{-0.0}, {-0.017249993519915588}, {108459.67947957986}, {1.6829766774912231E12}, {114537.26244584078}
               , {139520.20603565968}};

      MatrixStore<Double> l = Primitive64Store.FACTORY.makeWrapper(RawStore.wrap(l_));

      double[][] ae_ = new double[][]{{1.0, 0.0, 0.0, -5.941854438982693E-13, 0.0, 0.0}, {0.0, 1.0, 0.0, 1.0244148343231408E-14, 0.0,
               0.0}, {-0.0, -0.0, 1.0, 3.0066278187986717E-9, 0.0, 0.0}, {0.0, 0.0, 0.0, -9.791261189309351E-9, 1.0, 0.0}, {0.0, 0.0, 0.0
               , -1.8266341604461006E-10, 0.0, 1.0}};

      MatrixStore<Double> ae = Primitive64Store.FACTORY.makeWrapper(RawStore.wrap(ae_));

      double[][] be_ = new double[][]{{-1.000000103584958}, {0.017232766886716062}, {5774.467373370047}, {-22496.347961179275},
               {4298.272789515592}};

      MatrixStore<Double> be = Primitive64Store.FACTORY.makeWrapper(RawStore.wrap(be_));

      double[][] ai_ = new double[][]{{-1.0, -0.0, -0.0, -0.0, -0.0, -0.0}, {-0.0, -1.0, -0.0, -0.0, -0.0, -0.0}, {-1.0, -1.0, -0.0, -0.0
               , -0.0, -0.0}};

      MatrixStore<Double> ai = Primitive64Store.FACTORY.makeWrapper(RawStore.wrap(ai_));

      double[][] bi_ = new double[][]{{0.0}, {0.017249993519915588}, {0.017248993519915587}};

      MatrixStore<Double> bi = Primitive64Store.FACTORY.makeWrapper(RawStore.wrap(bi_));

      ConvexSolver.Builder builder = ConvexSolver.newBuilder();
      builder.objective(q, l);
      if (ae != null && be != null) {
         builder.equalities(ae, be);
      }
      if (ai != null && bi != null) {
         builder.inequalities(ai, bi);
      }
      Optimisation.Options options = new Optimisation.Options();
      options.iterations_abort = 100000;
      options.iterations_suffice = 10000;
      options.time_abort = 10000;
      options.time_suffice = 1000;
      options.sparse = false;
      options.validate = true;
      ConvexSolver.Configuration convex = options.convex();
      NumberContext present = convex.iterative();
      convex.iterative(NumberContext.of(8));
//      convex.solverSPD(Cholesky.R064::make).solverGeneral(LU.R064::make).iterative(NumberContext.of(12));
      convex.extendedPrecision(false);
      final ConvexSolver convexModel = builder.build(options);
      Optimisation.Result startValue = Optimisation.Result.of(0, Optimisation.State.APPROXIMATE, new double[q.getColDim()]);

      Optimisation.Result firstResult = null;
      try {
         firstResult = convexModel.solve(startValue);
      } catch (Exception e) {
         firstResult = convexModel.solve(startValue);
      }
      assertTrue(firstResult.getState().isSuccess());

   }


   @Test //Quadratic model that fail in package org.ojalgo.matrix.store, ojAlgo, version 54.0
   public void quadraticTest54() {

      double[][] q_ = new double[][]{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 1.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 1.0, 0.0, 0.0, 0.0}, {0.0
               , 0.0, 0.0, 1.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};

      MatrixStore<Double> q = R064Store.FACTORY.makeWrapper(RawStore.wrap(q_));

      double[][] l_ = new double[][]{{-0.0}, {-0.017249993519915588}, {108459.67947957986}, {1.6829766774912231E12}, {114537.26244584078}
               , {139520.20603565968}};

      MatrixStore<Double> l = R064Store.FACTORY.makeWrapper(RawStore.wrap(l_));

      double[][] ae_ = new double[][]{{1.0, 0.0, 0.0, -5.941854438982693E-13, 0.0, 0.0}, {0.0, 1.0, 0.0, 1.0244148343231408E-14, 0.0,
               0.0}, {-0.0, -0.0, 1.0, 3.0066278187986717E-9, 0.0, 0.0}, {0.0, 0.0, 0.0, -9.791261189309351E-9, 1.0, 0.0}, {0.0, 0.0, 0.0
               , -1.8266341604461006E-10, 0.0, 1.0}};

      MatrixStore<Double> ae = R064Store.FACTORY.makeWrapper(RawStore.wrap(ae_));

      double[][] be_ = new double[][]{{-1.000000103584958}, {0.017232766886716062}, {5774.467373370047}, {-22496.347961179275},
               {4298.272789515592}};

      MatrixStore<Double> be = R064Store.FACTORY.makeWrapper(RawStore.wrap(be_));

      double[][] ai_ = new double[][]{{-1.0, -0.0, -0.0, -0.0, -0.0, -0.0}, {-0.0, -1.0, -0.0, -0.0, -0.0, -0.0}, {-1.0, -1.0, -0.0, -0.0
               , -0.0, -0.0}};

      MatrixStore<Double> ai = R064Store.FACTORY.makeWrapper(RawStore.wrap(ai_));

      double[][] bi_ = new double[][]{{0.0}, {0.017249993519915588}, {0.017248993519915587}};

      MatrixStore<Double> bi = R064Store.FACTORY.makeWrapper(RawStore.wrap(bi_));

      ConvexSolver.Builder builder = ConvexSolver.newBuilder();
      builder.objective(q, l);
      if (ae != null && be != null) {
         builder.equalities(ae, be);
      }
      if (ai != null && bi != null) {
         builder.inequalities(ai, bi);
      }
      Optimisation.Options options = new Optimisation.Options();
      options.iterations_abort = 100000;
      options.iterations_suffice = 10000;
      options.time_abort = 10000;
      options.time_suffice = 1000;
      options.sparse = false;
      options.validate = true;
      ConvexSolver.Configuration convex = options.convex();
      NumberContext present = convex.iterative();
      convex.iterative(NumberContext.of(8));
      convex.solverSPD(Cholesky.R064::make).solverGeneral(LU.R064::make).iterative(NumberContext.of(12));
      convex.extendedPrecision(false);
      final ConvexSolver convexModel = builder.build(options);
      Optimisation.Result startValue = Optimisation.Result.of(0, Optimisation.State.APPROXIMATE, new double[q.getColDim()]);

      Optimisation.Result firstResult = null;
      try {
         firstResult = convexModel.solve(startValue);
      } catch (Exception e) {
         firstResult = convexModel.solve(startValue);
      }
      assertTrue(firstResult.getState().isSuccess());

   }

[apete]

Hej Magnus,

I briefly tried the v54 code and can confirm that the solver returns INFEASIBLE.

It is the LP solver, used to find the initial feasible solution, that concludes the problem to be infeasible.

I believe there's only minor tweaking to the LP solver between v53 and v54. The problem parameters are of greatly varying magnitude – 1.6829766774912231E12 and 1.0244148343231408E-14. Most likely some minor numerical difference cause the algorithm to interpret the numbers differently. I'll have a look at this in more detail later.

[Programmer-Magnus]

Thanks for your answer.

I modified the test above and try several modified model versions, by limiting small values in two different ways, using either the setSmallValuesToEpsilon or setSmallValuesToZero function, included here at the end. Setting small values to zero do not make this example work. But setting them to Epsilon give a result if epsilon is 1E-8.

Do you know a better way to get an approximate solution when the initial model fail?
If I use an Integration (e.g. CPLEX) will it be used in this case, or are they only used by ExpressionsBasedModel's?

Best Regards,
Magnus

` @test //Quadratic model that fail in package org.ojalgo.matrix.store, ojAlgo, version 54.0
public void quadraticTest6() {

  double[][] q_ = new double[][]{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 1.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 1.0, 0.0, 0.0, 0.0}, {0.0
           , 0.0, 0.0, 1.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};

  MatrixStore<Double> q = R064Store.FACTORY.makeWrapper(RawStore.wrap(q_));

  double[][] l_ = new double[][]{{-0.0}, {-0.017249993519915588}, {108459.67947957986}, {1.6829766774912231E12}, {114537.26244584078}
           , {139520.20603565968}};

  MatrixStore<Double> l = R064Store.FACTORY.makeWrapper(RawStore.wrap(l_));

  double[][] ae_ = new double[][]{{1.0, 0.0, 0.0, -5.941854438982693E-13, 0.0, 0.0}, {0.0, 1.0, 0.0, 1.0244148343231408E-14, 0.0,
           0.0}, {-0.0, -0.0, 1.0, 3.0066278187986717E-9, 0.0, 0.0}, {0.0, 0.0, 0.0, -9.791261189309351E-9, 1.0, 0.0}, {0.0, 0.0, 0.0
           , -1.8266341604461006E-10, 0.0, 1.0}};

  MatrixStore<Double> ae = R064Store.FACTORY.makeWrapper(RawStore.wrap(ae_));

  double[][] be_ = new double[][]{{-1.000000103584958}, {0.017232766886716062}, {5774.467373370047}, {-22496.347961179275},
           {4298.272789515592}};

  MatrixStore<Double> be = R064Store.FACTORY.makeWrapper(RawStore.wrap(be_));

  double[][] ai_ = new double[][]{{-1.0, -0.0, -0.0, -0.0, -0.0, -0.0}, {-0.0, -1.0, -0.0, -0.0, -0.0, -0.0}, {-1.0, -1.0, -0.0, -0.0
           , -0.0, -0.0}};

  MatrixStore<Double> ai = R064Store.FACTORY.makeWrapper(RawStore.wrap(ai_));

  double[][] bi_ = new double[][]{{0.0}, {0.017249993519915588}, {0.017248993519915587}};

  MatrixStore<Double> bi = R064Store.FACTORY.makeWrapper(RawStore.wrap(bi_));


  Optimisation.Options options = new Optimisation.Options();
  options.iterations_abort = 100000;
  options.iterations_suffice = 10000;
  options.time_abort = 10000;
  options.time_suffice = 1000;
  options.sparse = false;
  options.validate = true;
  ConvexSolver.Configuration convex = options.convex();
  NumberContext present = convex.iterative();
  convex.iterative(NumberContext.of(8));
  convex.solverSPD(Cholesky.R064::make).solverGeneral(LU.R064::make).iterative(NumberContext.of(12));
  convex.extendedPrecision(false);

  Optimisation.Result firstResult;
  double epsilon = 1.0E-20;
  do {
     epsilon *= 10;
     ConvexSolver.Builder builder = ConvexSolver.newBuilder();
     builder.objective(q, l);
     if (ae != null && be != null) {
        MatrixStore<Double> ae_nice = setSmallValuesToEpsilon(ae, epsilon);
        builder.equalities(ae_nice, be);
     }
     if (ai != null && bi != null) {
        MatrixStore<Double> ai_nice = setSmallValuesToEpsilon(ai, epsilon);
        builder.inequalities(ai_nice, bi);
     }
     java.util.logging.Logger.getLogger("org.ojalgo.optimisation").log(Level.INFO," ePSILON= " + epsilon);
     final ConvexSolver convexModel = builder.build(options);
     Optimisation.Result startValue = Optimisation.Result.of(0, Optimisation.State.APPROXIMATE, new double[q.getColDim()]);
     firstResult = convexModel.solve(startValue);
  } while (firstResult.getState().isFailure());

  assertTrue(firstResult.getState().isSuccess());

}

private MatrixStore setSmallValuesToZero(MatrixStore ae, double epsilon) {
long r = ae.countRows();
long c = ae.countColumns();
R064Store modified = FACTORY.newBuilder(r, c);
for (int row = 0; row < ae.countRows(); row++) {
for (int col = 0; col < ae.countColumns(); col++) {
Double value = ae.get(row, col);
if(value == 0) continue;
if (abs(value) < epsilon) {
modified.set(row, col, 0.0);
}else {
modified.set(row, col, value);
}
}
}
return modified;
}

private MatrixStore setSmallValuesToEpsilon(MatrixStore ae, double epsilon) {
long r = ae.countRows();
long c = ae.countColumns();
R064Store modified = FACTORY.newBuilder(r, c);
for (int row = 0; row < ae.countRows(); row++) {
for (int col = 0; col < ae.countColumns(); col++) {
Double value = ae.get(row, col);
if(value == 0) continue;
if (abs(value) < epsilon) {
modified.set(row, col, epsilon);
}else {
modified.set(row, col, value);
}
}
}
return modified;
}
`

[apete]

The 3:d party solver integrations (like CPLEX) are only available when using ExpressionsBasedModel.

I did try to build a model and solve the problem that way. Then the result is a solution very similar to what you got in your test case. BUT, if I solve that same model using CPLEX, CPLEX states the problem to be infeasible.