Now investigate and document positive initial null space with primal1.mps

Author: jajhall

check/TestQpSolver.cpp

@@ -6,7 +6,7 @@
 #include "catch.hpp"
 #include "io/FilereaderLp.h"
 
-const bool dev_run = true;//false;
+const bool dev_run = false;
 const double inf = kHighsInf;
 const double double_equal_tolerance = 1e-5;
 
@@ -954,7 +954,7 @@ TEST_CASE("test-qp-hot-start", "[qpsolver]") {
       const std::string filename =
           std::string(HIGHS_DIR) + "/check/instances/primal1.mps";
       REQUIRE(highs.readModel(filename) == HighsStatus::kOk);
-      highs.getBasis().print("test-qp-hot-start");
+      //      highs.getBasis().print("test-qp-hot-start");
       required_objective_function_value = -0.035012965733477348;
     } else if (k == 2) {
       // Not currently tested
@@ -1150,12 +1150,13 @@ TEST_CASE("rowless-qp", "[qpsolver]") {
 
 TEST_CASE("test-qp-atwood", "[qpsolver]") {
   Highs h;
-  //  highs.setOptionValue("output_flag", dev_run);
+  h.setOptionValue("output_flag", dev_run);
   std::string filename =
       std::string(HIGHS_DIR) + "/check/instances/atwood0.mps";
   REQUIRE(h.readModel(filename) == HighsStatus::kOk);
 
-  const double primal_feasibility_tolerance = h.getOptions().primal_feasibility_tolerance;
+  const double primal_feasibility_tolerance =
+      h.getOptions().primal_feasibility_tolerance;
 
   const double required_objective0 = 4.16347077e-02;
   const double required_objective1 = 2.91530651e-02;
@@ -1165,15 +1166,28 @@ TEST_CASE("test-qp-atwood", "[qpsolver]") {
   const HighsBasis& basis_ = h.getBasis();
   REQUIRE(model.lp_.row_lower_[1] == 0.26);
 
+  // After solving QP0, the lower bound on the second constraint is
+  // reduced to give QP1, so the optimal solution of QP0 is feasible
+  // for QP1. It should be possible to hot start QP1 using the optimal
+  // basis of QP0, but that basis has a null space dimension of 1, and
+  // the QP solver currently only hot starts from a vertex
+  // solution. This has a null space dimension of 0, although any free
+  // variables are considered as a special case (see primal1.mps),
+  // leading to a positive initial null space.
+  //
+  // Previously, this test case exposed a bug in hot start of the ASM
+  // solver, so it's a useful test case as/when hot start is improved.
+
   const bool hot_start = true;
-  for (HighsInt k= 0; k < 2; k++) {
+  for (HighsInt k = 0; k < 2; k++) {
     REQUIRE(h.run() == HighsStatus::kOk);
     HighsSolution solution = h.getSolution();
     HighsBasis basis = basis_;
 
     HighsInt num_basic = 0;
     HighsInt num_non_basic = 0;
     const HighsLp& lp = model.lp_;
+    /*
     printf("\nColumns\n");
     for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) {
       double lower = lp.col_lower_[iCol];
@@ -1182,16 +1196,13 @@ TEST_CASE("test-qp-atwood", "[qpsolver]") {
       double rsdu = std::min(value-lower, upper-value);
       HighsBasisStatus status = basis_.col_status[iCol];
       printf("%2d [%11.4g, %11.4g, %11.4g] %11.4g %s\n",
-	     int(iCol), lower, value, upper, rsdu, h.basisStatusToString(status).c_str());
-      if (status == HighsBasisStatus::kBasic) {
-	num_basic++;
-      } else if (status == HighsBasisStatus::kNonbasic) {
-	num_non_basic++;
-      } else {
-	REQUIRE(rsdu <= primal_feasibility_tolerance);
+             int(iCol), lower, value, upper, rsdu,
+    h.basisStatusToString(status).c_str()); if (status ==
+    HighsBasisStatus::kBasic) { num_basic++; } else if (status ==
+    HighsBasisStatus::kNonbasic) { num_non_basic++; } else { REQUIRE(rsdu <=
+    primal_feasibility_tolerance);
       }
     }
-  
     printf("Rows\n");
     for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) {
       double lower = lp.row_lower_[iRow];
@@ -1200,25 +1211,24 @@ TEST_CASE("test-qp-atwood", "[qpsolver]") {
       double rsdu = std::min(value-lower, upper-value);
       HighsBasisStatus status = basis_.row_status[iRow];
       printf("%2d [%11.4g, %11.4g, %11.4g] %11.4g %s\n",
-	   int(iRow), lower, value, upper, rsdu, h.basisStatusToString(status).c_str());
-      if (status == HighsBasisStatus::kBasic) {
-	num_basic++;
-      } else if (status == HighsBasisStatus::kNonbasic) {
-	num_non_basic++;
-      } else {
-	REQUIRE(rsdu <= primal_feasibility_tolerance);
+           int(iRow), lower, value, upper, rsdu,
+    h.basisStatusToString(status).c_str()); if (status ==
+    HighsBasisStatus::kBasic) { num_basic++; } else if (status ==
+    HighsBasisStatus::kNonbasic) { num_non_basic++; } else { REQUIRE(rsdu <=
+    primal_feasibility_tolerance);
       }
     }
-    printf("QP has %d basic and %d nonbasic variables\n", int(num_basic), int(num_non_basic));
-
+    printf("QP has %d basic and %d nonbasic variables\n", int(num_basic),
+    int(num_non_basic));
+    */
     if (k == 0) {
       const double objective0 = h.getInfo().objective_function_value;
       REQUIRE(std::fabs(objective0 - required_objective0) < 1e-4);
     } else {
       const double objective1 = h.getInfo().objective_function_value;
       REQUIRE(std::fabs(objective1 - required_objective1) < 1e-4);
-    }      
-  
+    }
+
     double lower = 0.25;
     double upper = kHighsInf;
     h.changeRowBounds(1, lower, upper);
@@ -1231,6 +1241,6 @@ TEST_CASE("test-qp-atwood", "[qpsolver]") {
       assert(basis_.valid);
     }
   }
- 
+
   h.resetGlobalScheduler(true);
 }

highs/lp_data/Highs.cpp

@@ -314,7 +314,6 @@ HighsStatus Highs::passModel(HighsModel model) {
   HighsStatus return_status = HighsStatus::kOk;
   // Clear the incumbent model and any associated data
   clearModel();
-  basis_.print("Highs::passModel(HighsModel model) after clearModel()");
   HighsLp& lp = model_.lp_;
   HighsHessian& hessian = model_.hessian_;
   // Move the model's LP and Hessian to the internal LP and Hessian
@@ -378,15 +377,12 @@ HighsStatus Highs::passModel(HighsModel model) {
   // Clear solver status, solution, basis and info associated with any
   // previous model; clear any HiGHS model object; create a HiGHS
   // model object for this LP
-  basis_.print("Highs::passModel(HighsModel model) before clearSolver()");
   return_status = interpretCallStatus(options_.log_options, clearSolver(),
                                       return_status, "clearSolver");
-  basis_.print("Highs::passModel(HighsModel model) after clearSolver()");
   // Apply any user scaling in call to optionChangeAction
   return_status =
       interpretCallStatus(options_.log_options, optionChangeAction(),
                           return_status, "optionChangeAction");
-  basis_.print("Highs::passModel(HighsModel model) exit");
   return returnFromHighs(return_status);
 }
 

highs/lp_data/HighsSolution.cpp

@@ -1815,25 +1815,7 @@ void HighsObjectiveSolution::clear() { this->col_value.clear(); }
 
 void HighsBasis::print(std::string message) const {
   this->printScalars(message);
-  if (!this->useful) {
-    if (this->col_status.size() > 0) {
-      HighsInt iCol = 0;
-      printf("Basis: col_status[%2d] = %d\n", int(iCol),
-	     int(this->col_status[iCol]));
-      iCol = this->col_status.size();
-      printf("Basis: col_status[%2d] = %d\n", int(iCol),
-	     int(this->col_status[iCol]));
-    }
-    if (this->row_status.size() > 0) {
-      HighsInt iRow = 0;
-      printf(" Basis: row_status[%2d] = %d\n", int(iRow),
-	     int(this->row_status[iRow]));
-      iRow = this->row_status.size();
-      printf(" Basis: row_status[%2d] = %d\n", int(iRow),
-	     int(this->row_status[iRow]));
-    }
-      return;
-  }
+  if (!this->useful) return;
   for (HighsInt iCol = 0; iCol < HighsInt(this->col_status.size()); iCol++)
     printf(" Basis: col_status[%2d] = %d\n", int(iCol),
            int(this->col_status[iCol]));

highs/qpsolver/a_quass.cpp

@@ -128,13 +128,15 @@ QpAsmStatus solveqp(Instance& instance, Settings& settings, Statistics& stats,
 
   QpSolution qp_solution(instance);
 
-  // presolve
-
-  // scale instance, store scaling factors
-
-  // perturb instance, store perturbance information
-
-  // regularize
+  // To be done!
+  //
+  // * Presolve
+  //
+  // * Scaling
+  //
+  // * Perturbation
+
+  // Regularize
   for (HighsInt i = 0; i < instance.num_var; i++) {
     for (HighsInt index = instance.Q.mat.start[i];
          index < instance.Q.mat.start[i + 1]; index++) {
@@ -144,7 +146,18 @@ QpAsmStatus solveqp(Instance& instance, Settings& settings, Statistics& stats,
     }
   }
 
-  // compute initial feasible point
+  // Determine an initial feasible point
+  //
+  // Cold start uses simplex to solve the LP feasibility problem to
+  // obtain a feasible vertex
+  //
+  // Hot start is currently limited to using the basis (and solution,
+  // if there is no basis) to hot start simplex when solving the LP
+  // feasibility problem
+  //
+  // It should be possible to hot start by using the QP basis to
+  // construct x = Yb, where Y and b correspond to the active
+  // constraints, and form the null space data structures
   QpHotstartInformation startinfo(instance.num_var, instance.num_con);
   if (instance.num_con == 0 && instance.num_var <= 15000) {
     computeStartingPointBounded(instance, settings, stats, qp_model_status,
@@ -162,10 +175,9 @@ QpAsmStatus solveqp(Instance& instance, Settings& settings, Statistics& stats,
     }
   } else {
     computeStartingPointByLp(instance, settings, stats, qp_model_status,
-			     startinfo,
-			     //			     highs_model_status,
-			     highs_basis,
-			     highs_solution, qp_timer);
+                             startinfo,
+                             //			     highs_model_status,
+                             highs_basis, highs_solution, qp_timer);
     if (qp_model_status != QpModelStatus::kNotset) {
       return quass2highs(instance, settings, stats, qp_model_status,
                          qp_solution, highs_model_status, highs_basis,
@@ -177,11 +189,13 @@ QpAsmStatus solveqp(Instance& instance, Settings& settings, Statistics& stats,
   solveqp_actual(instance, settings, startinfo, stats, qp_model_status,
                  qp_solution, qp_timer);
 
-  // undo perturbation and resolve
-
-  // undo scaling and resolve
-
-  // postsolve
+  // To be done!
+  //
+  // * Undo perturbation and resolve
+  //
+  // * Undo scaling and resolve
+  //
+  // * Postsolve
 
   // Transform QP status and qp_solution to HiGHS highs_basis and highs_solution
   return quass2highs(instance, settings, stats, qp_model_status, qp_solution,