Allow slightly asymmetric square Hessians, and revert to using HighsHashTable (mainly) In assessMatrix

check/TestFilereader.cpp

@@ -566,3 +566,48 @@ TEST_CASE("lp-duplicate-variable", "[highs_filereader]") {
 
   std::remove(lp_file.c_str());
 }
+
+inline double getWallTime() {
+  using namespace std::chrono;
+  return duration_cast<duration<double> >(
+             std::chrono::high_resolution_clock::now().time_since_epoch())
+      .count();
+}
+
+TEST_CASE("efficient-add-row", "[highs_filereader]") {
+  std::string filename;
+  filename = std::string(HIGHS_DIR) + "/check/instances/adlittle.mps";
+  //    "/srv/mps_da/neos.mps.gz";
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  REQUIRE(h.readModel(filename) == HighsStatus::kOk);
+  HighsLp lp = h.getLp();
+  h.passModel(lp.num_col_, 0, 0, 0, 0, 0, 0, 0.0, lp.col_cost_.data(),
+              lp.col_lower_.data(), lp.col_upper_.data(), nullptr, nullptr,
+              nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr);
+  lp.a_matrix_.ensureRowwise();
+  double tt = 0;
+  if (dev_run) tt = -getWallTime();
+  for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) {
+    HighsInt iEl = lp.a_matrix_.start_[iRow];
+    h.addRow(lp.row_lower_[iRow], lp.row_upper_[iRow],
+             lp.a_matrix_.start_[iRow + 1] - iEl, &lp.a_matrix_.index_[iEl],
+             &lp.a_matrix_.value_[iEl]);
+  }
+  if (dev_run) {
+    tt += getWallTime();
+    printf("Added %d rows individually in %.2gs\n", int(lp.num_row_), tt);
+  }
+
+  REQUIRE(h.deleteRows(0, lp.num_row_ - 1) == HighsStatus::kOk);
+  REQUIRE(h.getLp().num_row_ == 0);
+
+  if (dev_run) tt = -getWallTime();
+  h.addRows(lp.num_row_, lp.row_lower_.data(), lp.row_upper_.data(),
+            lp.a_matrix_.start_[lp.num_row_], lp.a_matrix_.start_.data(),
+            lp.a_matrix_.index_.data(), lp.a_matrix_.value_.data());
+  if (dev_run) {
+    tt += getWallTime();
+    printf("Added %d rows together     in %.2gs\n", int(lp.num_row_), tt);
+  }
+}

check/TestHighsHessian.cpp

@@ -9,7 +9,7 @@ const bool dev_run = false;
 // No commas in test case name.
 TEST_CASE("HighsHessian", "[highs_hessian]") {
   HighsOptions options;
-  if (!dev_run) options.output_flag = false;
+  options.output_flag = dev_run;
 
   HighsHessian square_hessian;
   // .  0  1  2  3  4
@@ -327,3 +327,28 @@ TEST_CASE("HighsHessian", "[highs_hessian]") {
     triangular_hessian.print();
   }
 }
+
+TEST_CASE("square-near-symmetric-hessian", "[highs_hessian]") {
+  const double epsilon = 1e-15;
+  HighsOptions options;
+  options.output_flag = dev_run;
+
+  HighsHessian square_hessian;
+  // .  0  1  2
+  // 0  5  1   1+eps
+  // 1  1  4  0
+  // 2  1  0  3
+
+  square_hessian.dim_ = 3;
+  square_hessian.format_ = HessianFormat::kSquare;
+  square_hessian.start_ = {0, 4, 6, 8};
+  square_hessian.index_ = {0, 1, 2, 1, 0, 1, 0, 2};
+  square_hessian.value_ = {5, 1, 1, 0.1, 1.1, 4, 1 + epsilon, 3};
+  if (dev_run)
+    printf(
+        "square-near-symmetric-hessian: 0-1 and 1-0 values are %.24g and "
+        "%.24g, with difference %g\n",
+        square_hessian.value_[2], square_hessian.value_[6],
+        std::fabs(square_hessian.value_[2] - square_hessian.value_[6]));
+  REQUIRE(assessHessian(square_hessian, options) == HighsStatus::kOk);
+}

highs/model/HighsHessianUtils.cpp

@@ -403,10 +403,17 @@ HighsStatus normaliseHessian(const HighsOptions& options,
   upper_off_diagonal.assign(dim, 0);
   // Should be no duplicates
   HighsInt debug_num_duplicate = 0;
-  HighsInt num_non_symmetric = 0;
   HighsInt num_summation = 0;
   HighsInt num_upper_triangle = 0;
   HighsInt num_hessian_el = 0;
+  const double kSquareHessianAsymmetryTolerance = 1e-10;
+  HighsInt num_illegal_asymmetry = 0;
+  double min_illegal_asymmetry = kHighsInf;
+  double max_illegal_asymmetry = 0;
+  HighsInt num_ok_asymmetry = 0;
+  double min_ok_asymmetry = kHighsInf;
+  double max_ok_asymmetry = 0;
+
   //  const bool expensive_2821_check = true;
   for (HighsInt iCol = 0; iCol < dim; iCol++) {
     for (HighsInt iEl = from_hessian.start_[iCol];
@@ -427,9 +434,23 @@ HighsStatus normaliseHessian(const HighsOptions& options,
       upper_off_diagonal[iRow] = upper.value_[iEl];
       if (square) {
         // When square, ensure that the upper triangular value matches
-        // the corresponding lower triangular value
-        if (upper_off_diagonal[iRow] != lower_on_below_diagonal[iRow])
-          num_non_symmetric++;
+        // the corresponding lower triangular value to within the
+        // tolaernace used by JuMP
+        const double asymmetry =
+            std::fabs(upper_off_diagonal[iRow] - lower_on_below_diagonal[iRow]);
+        if (asymmetry > kSquareHessianAsymmetryTolerance) {
+          num_illegal_asymmetry++;
+          min_illegal_asymmetry = std::min(asymmetry, min_illegal_asymmetry);
+          max_illegal_asymmetry = std::max(asymmetry, max_illegal_asymmetry);
+        } else if (asymmetry) {
+          num_ok_asymmetry++;
+          min_ok_asymmetry = std::min(asymmetry, min_ok_asymmetry);
+          max_ok_asymmetry = std::max(asymmetry, max_ok_asymmetry);
+          double average =
+              (upper_off_diagonal[iRow] + lower_on_below_diagonal[iRow]) * 0.5;
+          upper_off_diagonal[iRow] = average;
+          lower_on_below_diagonal[iRow] = average;
+        }
         // Don't zero the upper off diagonal entry, so that it's
         // possible to check that nonzeros in the lower off diagonal
         // match the upper off diagonal entry
@@ -505,11 +526,24 @@ HighsStatus normaliseHessian(const HighsOptions& options,
 
   bool warning_found = false;
   bool error_found = false;
-  if (num_non_symmetric) {
+  if (num_ok_asymmetry) {
+    assert(square);
+    highsLogUser(options.log_options, HighsLogType::kInfo,
+                 "Square Hessian contains %d non-symmetr%s in [%.2g, %.2g] "
+                 "within tolerance of %.1g\n",
+                 int(num_ok_asymmetry), num_ok_asymmetry == 1 ? "y" : "ies",
+                 min_ok_asymmetry, max_ok_asymmetry,
+                 kSquareHessianAsymmetryTolerance);
+  }
+  if (num_illegal_asymmetry) {
     assert(square);
     highsLogUser(options.log_options, HighsLogType::kError,
-                 "Square Hessian contains %d non-symmetr%s\n",
-                 int(num_non_symmetric), num_non_symmetric == 1 ? "y" : "ies");
+                 "Square Hessian contains %d non-symmetr%s in [%.2g, %.2g] "
+                 "exceeding tolerance of %.1g\n",
+                 int(num_illegal_asymmetry),
+                 num_illegal_asymmetry == 1 ? "y" : "ies",
+                 min_illegal_asymmetry, max_illegal_asymmetry,
+                 kSquareHessianAsymmetryTolerance);
     error_found = true;
   }
   if (num_upper_triangle) {

highs/util/HighsMatrixUtils.cpp

@@ -143,11 +143,14 @@ HighsStatus assessMatrix(
   double max_large_value = 0;
   double min_large_value = kHighsInf;
   HighsInt num_duplicate = 0;
-  // Use index_map to identify duplicates.
-  //  HighsHashTable<HighsInt> index_set;
-  std::vector<HighsInt> el_in_vec;
   const HighsInt illegal_el = -1;
-  el_in_vec.assign(vec_dim, illegal_el);
+  // When duplicates only have to be identified (and not summed) use
+  // highs_hash
+  const bool use_highs_hash = !sum_duplicates;
+  HighsHashTable<HighsInt> highs_hash;
+  // Otherwise, use index-key map to identify entry number of
+  // duplicates so they can be summed
+  std::unordered_map<HighsInt, HighsInt> index_el_map;
 
   for (HighsInt ix = 0; ix < num_vec; ix++) {
     HighsInt from_el = matrix_start[ix];
@@ -180,26 +183,22 @@ HighsStatus assessMatrix(
                      matrix_name.c_str(), ix, el, component, vec_dim);
         return HighsStatus::kError;
       }
-      // Check that the index has not already occurred.
-      HighsInt previous_el = el_in_vec[component];
-      bool is_duplicate = previous_el > illegal_el;
-      // 2821 eliminates need for index_set
-      /*
-      legal_component = index_set.find(component) == nullptr;
-      if (legal_component != !is_duplicate) {
-        printf(
-            "assessMatrix: ix = %d; el_in_vec[%d/%d] = %d; legal_component = "
-            "%d\n",
-            int(ix), int(component), int(vec_dim), int(previous_el),
-            legal_component);
+      // Check whether the index has already occurred.
+      HighsInt previous_el = illegal_el;
+      bool is_duplicate = false;
+      if (use_highs_hash) {
+        is_duplicate = highs_hash.find(component) != nullptr;
+      } else {
+        auto found_component = index_el_map.find(component);
+        is_duplicate = found_component != index_el_map.end();
+        if (is_duplicate) previous_el = found_component->second;
       }
-      assert(legal_component == !is_duplicate);
-      */
       if (is_duplicate) {
         if (sum_duplicates) {
           num_duplicate++;
-          // Sum the duplicate entry
-          assert(matrix_index[previous_el] == component);
+          // Sum the duplicate entry, making sure that it's been
+          // assigned
+          assert(previous_el != illegal_el);
           matrix_value[previous_el] += matrix_value[el];
           continue;
         }
@@ -216,35 +215,21 @@ HighsStatus assessMatrix(
       // the new number of nonzeros
       matrix_index[num_new_nz] = matrix_index[el];
       matrix_value[num_new_nz] = matrix_value[el];
-      // Record where the index has occurred
-      //      index_set.insert(component);
-      el_in_vec[component] = num_new_nz;
+      if (use_highs_hash) {
+        // Record that the index has occurred
+        highs_hash.insert(component);
+      } else {
+        // Record where the index has occurred
+        index_el_map.insert({component, num_new_nz});
+      }
       num_new_nz++;
     }
     from_el = matrix_start[ix];
     to_el = num_new_nz;
     // Reset num_new_nz
     num_new_nz = matrix_start[ix];
-    // Reset el_in_vec
-    for (HighsInt el = from_el; el < to_el; el++)
-      el_in_vec[matrix_index[el]] = illegal_el;
-    //    const bool expensive_2821_check = true;
-    /*
-    if (expensive_2821_check) {
-      // Check el_in_vec !! Remove later!
-      for (HighsInt lc_ix = 0; lc_ix < vec_dim; lc_ix++)
-        assert(el_in_vec[lc_ix] == illegal_el);
-    }
-    */
     for (HighsInt el = from_el; el < to_el; el++) {
       HighsInt component = matrix_index[el];
-      /*
-      if (expensive_2821_check) {
-        // Ensure that duplicates have been eliminated
-        HighsInt previous_el = el_in_vec[component];
-        assert(previous_el == illegal_el);
-      }
-      */
       // Check the value
       double abs_value = fabs(matrix_value[el]);
       // Check that the value is not too large
@@ -266,26 +251,14 @@ HighsStatus assessMatrix(
         // the new number of nonzeros
         matrix_index[num_new_nz] = matrix_index[el];
         matrix_value[num_new_nz] = matrix_value[el];
-        /*
-        if (expensive_2821_check) {
-          // Record where the index has occurred
-          el_in_vec[component] = num_new_nz;
-        }
-        */
         num_new_nz++;
       }
     }  // Loop from_el; to_el
-    //    index_set.clear();
-    /*
-    if (expensive_2821_check) {
-      // Reset el_in_vec
-      for (HighsInt el = from_el; el < to_el; el++)
-        el_in_vec[matrix_index[el]] = illegal_el;
-      // Check el_in_vec !! Remove later!
-      for (HighsInt lc_ix = 0; lc_ix < vec_dim; lc_ix++)
-        assert(el_in_vec[lc_ix] == illegal_el);
+    if (use_highs_hash) {
+      highs_hash.clear();
+    } else {
+      index_el_map.clear();
     }
-    */
   }  // Loop 0; num_vec
   if (num_duplicate) {
     highsLogUser(log_options, HighsLogType::kInfo,