Warm start after changing row range

[TendTo]

In my use case, I need to solve the same LP problem many times, enabling and disabling the same constraints as required. I noticed that more often than not the results after a previous solve were incorrect, as if they were stuck in the previous basis.

In fact, forcing a refresh with clearBasis between each solve fixes the issue, but that means having a cold start every time. I was wondering if this behaviour is intended or something is going wrong within the internals of the solver.

The following is minimal snippet to reproduce the problem. I have just started investigating it, so I probably missed something.

#include "soplex.h"

#include <iostream>

using namespace soplex;

void check_and_print_solution(SoPlex &soplex, int expected) {
  SPxSolver::Status stat;
  DVectorRational prim(2);
  DVectorRational dual(1);

  stat = soplex.optimize();
  if (stat != SPxSolver::OPTIMAL || soplex.objValueRational() != expected) {
    throw std::runtime_error((std::stringstream{}
                              << "SoPlex returned with status "
                              << std::to_string(stat) << " and objective value "
                              << soplex.objValueRational() << " instead of "
                              << expected)
                                 .str());
  }

  soplex.getPrimalRational(prim);
  soplex.getDualRational(dual);
  std::cout << "LP solved to optimality.\n";
  std::cout << "Objective value is " << soplex.objValueRational() << ".\n";
  std::cout << "Primal solution is [" << prim[0] << ", " << prim[1] << "].\n";
  std::cout << "Dual solution is [" << dual[0] << "].\n";
}

void test_rational() {
  SoPlex mysoplex;

  /* set parameters for exact solving */
  mysoplex.setIntParam(SoPlex::READMODE, SoPlex::READMODE_RATIONAL);
  mysoplex.setIntParam(SoPlex::SOLVEMODE, SoPlex::SOLVEMODE_RATIONAL);
  mysoplex.setIntParam(SoPlex::CHECKMODE, SoPlex::CHECKMODE_RATIONAL);
  mysoplex.setIntParam(SoPlex::SYNCMODE, SoPlex::SYNCMODE_AUTO);
  mysoplex.setIntParam(SoPlex::SIMPLIFIER, SoPlex::SIMPLIFIER_OFF);
  mysoplex.setRealParam(SoPlex::FEASTOL, 0.0);
  mysoplex.setRealParam(SoPlex::OPTTOL, 0.0);

  /* set the objective sense (MAXIMIZE) */
  mysoplex.setIntParam(SoPlex::OBJSENSE, SoPlex::OBJSENSE_MAXIMIZE);

  /* we first add variables (the integer data is converted to type Rational) */
  mysoplex.addColRational(LPColRational(0, DSVectorRational(), infinity, 0)); // 0 <= x0        | 0 coeff
  mysoplex.addColRational(LPColRational(1, DSVectorRational(), 1, 0));        // 0 <= x1 <= 1   | 1 coeff

  /* then constraints one by one (here we show how Rationals can be used
   * directly) */
  DSVectorRational row1(2);
  row1.add(0, 1);
  row1.add(1, 1);
  mysoplex.addRowRational(LPRowRational(-infinity, row1, infinity)); // C0: x0 + x1 <= infinity

  // Maximize
  //   obj: 1 x1
  // Subject To
  //  C0 : 1 x0 + 1 x1 <= infinity
  // Bounds
  //   x1 <= 1
  // End

  check_and_print_solution(mysoplex, 1);

  // Maximize
  //   obj: 1 x1
  // Subject To
  //  C0 : 1 x0 + 1 x1 <= 0
  // Bounds
  //   x1 <= 1
  // End

  // Enable C0
  mysoplex.changeRangeRational(0, -infinity, 0); // C0: x0 + x1 <= 0
  check_and_print_solution(mysoplex, 0);

  // Maximize
  //   obj: 1 x1
  // Subject To
  //  C0 : 1 x0 + 1 x1 <= infinity
  // Bounds
  //   x1 <= 1
  // End

  // Disable C0 again
  mysoplex.changeRangeRational(0, -infinity, infinity); // C0: x0 + x1 <= infinity

  // mysoplex.clearBasis(); // Adding this fixes the issue, but it's very expensive, especially for repeated solves
  //
  // mysoplex.syncLPRational(); // Gets ignored because of SYNCMODE_AUTO
  //
  // mysoplex._syncLPRational(); // Does nothing
  //
  // mysoplex._solver.reLoad();  // Does nothing
  //
  // mysoplex._hasBasis = false; // Solve the issue, but requires private access and still expensive

  // ERROR: we don't get 1 anymore (like in the original problem), but 0 (like the last one)
  check_and_print_solution(mysoplex, 1); }

int main() {
  test_rational();
  return 0;
}

[leoneifler]

Thanks for posting this well-defined issue. I will look into it. I'm guessing this is on the current master?

Just to note: Warm-starting (while still beneficial) is not as impactful in rational solving as in normal double-precision mode.

[leoneifler]

Yeah this is definitely a bug in SoPlex (not even the rational part, the same happens if you solve with tolerances). It sets the row to be nonbasic-free after you deactivate it, but the implementation does not consider these types for entering the basis apparently. I will try to fix this, but it might take some time. As a workaround in the meantime, you could set the rhs to some big M (effectively also deactivating it), and it will work correctly.

[TendTo]

Thanks, I'll quickly update you:

• Yes, the version I'm using in this example is the public master branch here on GitHub
• I noticed the error still happens even in real mode. I have added the slightly modified functions. It may be simpler to start looking at this example
• While the performance gain may not be very noticeable, even a tiny improvement can make a huge difference on repeated solves. In some instances, re-initialising the base is quite expensive. Reducing the cost by a small percentage would help, without even taking into account the ideally faster successive solves

If I find any other insights, I'll share them here.

void check_and_print_solution_real(SoPlex &soplex, double expected) {
  SPxSolver::Status stat;
  VectorReal prim(2);
  VectorReal dual(1);

  stat = soplex.optimize();
  if (stat != SPxSolver::OPTIMAL || soplex.objValueRational() != expected) {
    throw std::runtime_error((std::stringstream{} << "SoPlex returned with status " << std::to_string(stat)
                                                  << " and objective value " << soplex.objValueRational()
                                                  << " instead of " << expected)
                                 .str());
  }

  soplex.getPrimal(prim);
  soplex.getDual(dual);
  std::cout << "LP solved to optimality.\n";
  std::cout << "Objective value is " << soplex.objValueRational() << ".\n";
  std::cout << "Primal solution is [" << prim[0] << ", " << prim[1] << "].\n";
  std::cout << "Dual solution is [" << dual[0] << "].\n";
}

void test_real() {
  SoPlex mysoplex;
  /* set the objective sense (MAXIMIZE) */
  mysoplex.setIntParam(SoPlex::OBJSENSE, SoPlex::OBJSENSE_MAXIMIZE);

  /* we first add variables (the integer data is converted to type Rational) */
  mysoplex.addColReal(LPColReal(0, DSVectorReal(), infinity, 0)); // 0 <= x0        | 0 coeff
  mysoplex.addColReal(LPColReal(1, DSVectorReal(), 1, 0));        // 0 <= x1 <= 1   | 1 coeff

  /* then constraints one by one (here we show how Rationals can be used
   * directly) */
  DSVectorReal row1(2);
  row1.add(0, 1);
  row1.add(1, 1);
  mysoplex.addRowReal(LPRowReal(-infinity, row1, infinity)); // C0: 1 + 1 <= infinity

  // Maximize
  //   obj: 1 x1
  // Subject To
  //  C0 : 1 x0 + 1 x1 <= infinity
  // Bounds
  //   x1 <= 1
  // End

  check_and_print_solution_real(mysoplex, 1);

  // Maximize
  //   obj: 1 x1
  // Subject To
  //  C0 : 1 x0 + 1 x1 <= 0
  // Bounds
  //   x1 <= 1
  // End

  // Enable C0
  mysoplex.changeRangeReal(0, -infinity, 0); // C0: 1 + 1 <= 0
  check_and_print_solution_real(mysoplex, 0);

  // Maximize
  //   obj: 1 x1
  // Subject To
  //  C0 : 1 x0 + 1 x1 <= infinity
  // Bounds
  //   x1 <= 1
  // End

  // Disable C0 again
  mysoplex.changeRangeReal(0, -infinity, infinity); // C1: 1 + 0 <= infinity

  // ATTEMPTS

  // Adding this fixes the issue, but it's very expensive, especially for repeated solves
  // mysoplex.clearBasis();

  // Gets ignored because of SYNCMODE_AUTO
  // mysoplex.syncLPReal();

  // Does nothing
  // mysoplex._syncLPReal();

  // Does nothing
  // mysoplex._solver.reLoad();

  // Solve the issue, but requires private access and still expensive
  // mysoplex._hasBasis = false;

  // ERROR: we don't get 1 anymore (like in the original problem), but 0 (like the last one)
  check_and_print_solution_real(mysoplex, 1);
}

[ambros-gleixner]

You could try the following workaround (not even a bad one, maybe): Set the rhs to the max. activity of the constraint + small positive value (0.1, e.g.). By max. activity I mean any upper bound on a^T x over the feasible region (assuming the constraint reads a^T x <= b). This effectively removes the constraint, should avoid SoPlex's problems with infinity, and give you a hotstart. It is trivially possible if all your variables have finite bounds.