Reintroduced spurious unboundedness test in feasibility_bounded.cpp

Author: Julian Hall

check/TestQpSolver.cpp

@@ -57,24 +57,41 @@ void testPrimalDualObjective(Highs& h,
   REQUIRE(fabs(info.primal_dual_objective_error) < optimality_tolerance);
 }
 
-/*
 TEST_CASE("qp-unbounded", "[qpsolver]") {
   std::string filename;
   filename = std::string(HIGHS_DIR) + "/check/instances/qpunbounded.lp";
-
+  // Maximize
+  //  obj: x2 + [ - x1^2 ]/2
+  // Bounds
+  //  0 <= x1 <= 40
+  // End
   Highs highs;
   highs.setOptionValue("output_flag", dev_run);
   REQUIRE(highs.readModel(filename) == HighsStatus::kOk);
 
-  for (auto& solver : solvers) {
-    highs.setOptionValue("solver", solver);
-    REQUIRE(highs.run() == HighsStatus::kOk);
-    REQUIRE(highs.getModelStatus() == unboundedStatus(solver));
+  for (HighsInt k = 0; k < 2; k++) {
+    for (auto& solver : solvers) {
+      highs.setOptionValue("solver", solver);
+      REQUIRE(highs.run() == HighsStatus::kOk);
+      REQUIRE(highs.getModelStatus() == unboundedStatus(solver));
+      // Test the oracle Hessian
+      highs.setOptionValue("test_qp_oracle", true);
+      REQUIRE(highs.run() == HighsStatus::kOk);
+      REQUIRE(highs.getModelStatus() == unboundedStatus(solver));
+      highs.setOptionValue("test_qp_oracle", false);
+    }
+    // Now flip the cost and bounds of x2 (which has index 0 since it
+    // appears first in LP file objective) for code coverage
+    REQUIRE(highs.changeColCost(0, -highs.getLp().col_cost_[0]) ==
+            HighsStatus::kOk);
+    REQUIRE(highs.changeColBounds(0, -highs.getLp().col_upper_[0],
+                                  -highs.getLp().col_lower_[0]) ==
+            HighsStatus::kOk);
+    highs.clearSolver();
   }
 
   highs.resetGlobalScheduler(true);
 }
-*/
 
 TEST_CASE("qp-infeasible", "[qpsolver]") {
   std::string filename;
@@ -216,7 +233,6 @@ TEST_CASE("qpsolver", "[qpsolver]") {
   highs.resetGlobalScheduler(true);
 }
 
-/*
 TEST_CASE("test-qod", "[qpsolver]") {
   HighsStatus return_status;
   HighsModelStatus model_status;
@@ -225,6 +241,7 @@ TEST_CASE("test-qod", "[qpsolver]") {
   double required_x1;
 
   Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
 
   for (auto& solver : solvers) {
     HighsModel local_model;
@@ -250,7 +267,6 @@ TEST_CASE("test-qod", "[qpsolver]") {
     hessian.index_ = {0};
     hessian.value_ = {2.0};
 
-    highs.setOptionValue("output_flag", dev_run);
     const HighsInfo& info = highs.getInfo();
     const HighsSolution& solution = highs.getSolution();
     const double& objective_function_value = info.objective_function_value;
@@ -369,7 +385,6 @@ TEST_CASE("test-qod", "[qpsolver]") {
 
   highs.resetGlobalScheduler(true);
 }
-*/
 
 TEST_CASE("test-qjh", "[qpsolver]") {
   // Test passing/reading and solving the problem qjh
@@ -382,6 +397,7 @@ TEST_CASE("test-qjh", "[qpsolver]") {
   double required_objective_function_value;
 
   Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
 
   for (auto& solver : solvers) {
     HighsModel local_model;
@@ -408,7 +424,6 @@ TEST_CASE("test-qjh", "[qpsolver]") {
     hessian.index_ = {0, 2, 1, 2};
     hessian.value_ = {2.0, -1.0, 0.2, 2.0};
 
-    highs.setOptionValue("output_flag", dev_run);
     highs.setOptionValue("solver", solver);
     const HighsInfo& info = highs.getInfo();
     const double& objective_function_value = info.objective_function_value;

highs/lp_data/HStruct.h

@@ -257,8 +257,8 @@ struct HessianOracle {
   // form products with Q and then apply any non-trivial multiplier_
   // or shift_
   HighsInt dim_ = 0;
-  HighsInt multiplier_ = 1.0;  // Minimize
-  HighsInt shift_ = 0.0;       // No regularization
+  double multiplier_ = 1.0;  // Minimize
+  double shift_ = 0.0;       // No regularization
   HighsHessianFunctionType call_ = nullptr;
   void* data_ = nullptr;
   void clear();

highs/lp_data/Highs.cpp

@@ -4323,6 +4323,8 @@ HighsStatus Highs::callSolveQp() {
     instance.con_up = lp.row_upper_;
     instance.var_lo = lp.col_lower_;
     instance.var_up = lp.col_upper_;
+    // Clear the instance Hessian data
+    instance.Q.mat.clear();
     HighsHessian oracle_hessian;
     if (hessian.dim_ > 0) {
       if (options_.test_qp_oracle) {
@@ -4331,9 +4333,6 @@ HighsStatus Highs::callSolveQp() {
         instance.Q.mat.oracle_.dim_ = hessian.dim_;
         instance.Q.mat.oracle_.call_ = oracleCallSquareHessian;
         instance.Q.mat.oracle_.data_ = &oracle_hessian;
-        instance.Q.mat.start.clear();
-        instance.Q.mat.index.clear();
-        instance.Q.mat.value.clear();
       } else {
         assert(instance.Q.mat.oracle_.call_ == nullptr);
         instance.Q.mat.num_col = lp.num_col_;
@@ -4344,8 +4343,6 @@ HighsStatus Highs::callSolveQp() {
     } else {
       assert(hessian.isOracle());
       instance.Q.mat.oracle_ = hessian.oracle_;
-      instance.Q.mat.num_col = 0;
-      instance.Q.mat.num_row = 0;
     }
 
     for (HighsInt i = 0; i < (HighsInt)instance.c.value.size(); i++) {

highs/qpsolver/a_quass.cpp

@@ -52,6 +52,9 @@ static QpAsmStatus quass2highs(Instance& instance, Settings& settings,
     case QpModelStatus::kLargeNullspace:
       highs_model_status = HighsModelStatus::kSolveError;
       return QpAsmStatus::kError;
+    case QpModelStatus::kNonConvex:
+      highs_model_status = HighsModelStatus::kSolveError;
+      return QpAsmStatus::kError;
     case QpModelStatus::kError:
       highs_model_status = HighsModelStatus::kSolveError;
       return QpAsmStatus::kError;
@@ -60,6 +63,7 @@ static QpAsmStatus quass2highs(Instance& instance, Settings& settings,
       return QpAsmStatus::kError;
     default:
       highs_model_status = HighsModelStatus::kNotset;
+      assert(highs_model_status != HighsModelStatus::kNotset);
       return QpAsmStatus::kError;
   }
 

highs/qpsolver/feasibility_bounded.hpp

@@ -18,7 +18,7 @@ static void computeStartingPointBounded(Instance& instance, Settings& settings,
                                         HighsTimer& timer) {
   // compute initial feasible point for problems with bounds only (no general
   // linear constraints)
-
+  const bool debug_printing = false;
   // Solve  Qx + c = 0 --> x = -Q^-1c
   HighsInt dim = instance.num_var;
   QpVector res = -instance.c;
@@ -40,6 +40,7 @@ static void computeStartingPointBounded(Instance& instance, Settings& settings,
     L.resize(dim * dim);
 
     auto printL = [&]() {
+      if (!debug_printing) return;
       for (HighsInt iRow = 0; iRow < dim; iRow++) {
         for (HighsInt iCol = 0; iCol <= iRow; iCol++)
           printf(" %11.4g", L[iRow * dim + iCol]);
@@ -71,17 +72,20 @@ static void computeStartingPointBounded(Instance& instance, Settings& settings,
           sum += L[iRow * dim + k] * L[iCol * dim + k];
         if (iCol < iRow) {
           double value = (L[iRow * dim + iCol] - sum) / L[iCol * dim + iCol];
-          printf("Computed L[%1d, %1d] = (%11.4g - %11.4g) / %11.4g = %11.4g\n",
-                 int(iRow), int(iCol), L[iRow * dim + iCol], sum,
-                 L[iCol * dim + iCol], value);
+          if (debug_printing)
+            printf(
+                "Computed L[%1d, %1d] = (%11.4g - %11.4g) / %11.4g = %11.4g\n",
+                int(iRow), int(iCol), L[iRow * dim + iCol], sum,
+                L[iCol * dim + iCol], value);
           L[iRow * dim + iCol] = value;
         } else {
           double value = L[iCol * dim + iCol] - sum;
-          printf(
-              "Computed L[%1d, %1d] = sqrt(%11.4g - %11.4g) = sqrt(%11.4g) = "
-              "%11.4g\n",
-              int(iCol), int(iCol), L[iCol * dim + iCol], sum, value,
-              std::sqrt(value));
+          if (debug_printing)
+            printf(
+                "Computed L[%1d, %1d] = sqrt(%11.4g - %11.4g) = sqrt(%11.4g) = "
+                "%11.4g\n",
+                int(iCol), int(iCol), L[iCol * dim + iCol], sum, value,
+                std::sqrt(value));
           if (value <= 0) {
             modelstatus = QpModelStatus::kNonConvex;
             return;
@@ -115,12 +119,21 @@ static void computeStartingPointBounded(Instance& instance, Settings& settings,
   std::vector<BasisStatus> atlower;
 
   for (int i = 0; i < instance.num_var; i++) {
-    // Deleted spurious identification of unboundedness based on
-    // res.value[i] exceeding
-    // 0.5/settings.hessian_regularization_value, where the
-    // denominator can be zero
-    //
-    // Clearly need something for qp-unbounded and test-qod
+    // Check for unboundedness so that qp-unbounded and test-qod pass:
+    // assumes diagonal Hessian and
+    // settings.hessian_regularization_value > 0
+    if (settings.hessian_regularization_value > 0) {
+      if (res.value[i] > 0.5 / settings.hessian_regularization_value &&
+          instance.var_up[i] == kHighsInf && instance.c.value[i] < 0.0) {
+        modelstatus = QpModelStatus::kUnbounded;
+        return;
+      } else if (res.value[i] < -0.5 / settings.hessian_regularization_value &&
+                 instance.var_lo[i] == -kHighsInf &&
+                 instance.c.value[i] > 0.0) {
+        modelstatus = QpModelStatus::kUnbounded;
+        return;
+      }
+    }
     if (res.value[i] <
         instance.var_lo[i] - settings.primal_feasibility_tolerance) {
       res.value[i] = instance.var_lo[i];

highs/qpsolver/matrix.hpp

@@ -22,6 +22,15 @@ struct MatrixBase {
   std::vector<double> value;
   HessianOracle oracle_;
 
+  void clear() {
+    num_row = 0;
+    num_col = 0;
+    start.clear();
+    index.clear();
+    value.clear();
+    oracle_.clear();
+  }
+
   bool isOracle() const { return oracle_.call_ != nullptr; }
 
   double diagonal(const HighsInt iCol) const {