GF2 presolve rank-deficiency fix

cpp/src/mip_heuristics/presolve/gf2_presolve.cpp

@@ -8,8 +8,10 @@
 #include "gf2_presolve.hpp"
 
 #include <mip_heuristics/mip_constants.hpp>
+#include <utilities/macros.cuh>
 
 #include <cmath>
+#include <cstdint>
 #include <unordered_map>
 
 #if GF2_PRESOLVE_DEBUG
@@ -33,55 +35,103 @@ static inline i_t positive_modulo(i_t i, i_t n)
   return (i % n + n) % n;
 }
 
+static constexpr int GF2_WORD_BITS = 64;
+
+static inline int gf2_nwords(int N) { return (N + GF2_WORD_BITS - 1) / GF2_WORD_BITS; }
+
+static inline bool gf2_test_bit(const std::vector<uint64_t>& row, int col)
+{
+  return (row[col / GF2_WORD_BITS] >> (col % GF2_WORD_BITS)) & uint64_t{1};
+}
+
+static inline void gf2_set_bit(std::vector<uint64_t>& row, int col)
+{
+  row[col / GF2_WORD_BITS] |= (uint64_t{1} << (col % GF2_WORD_BITS));
+}
+
 // this is kind-of a stopgap implementation (as in practice MIPLIB2017 only contains a couple of GF2
 // problems and they're small) but cuDSS could be used for this since A is likely to be sparse and
 // low-bandwidth (i think?) unlikely to occur in real-world problems however. doubt it'd be worth
-// the effort trashes A and b, return true if solved
-static bool gf2_solve(std::vector<std::vector<int>>& A, std::vector<int>& b, std::vector<int>& x)
+// the effort
+gf2_status_t gf2_solve(std::vector<std::vector<uint64_t>>& A,
+                       int n_cols,
+                       std::vector<int>& b,
+                       std::vector<int>& x,
+                       std::vector<uint8_t>& determined)
 {
-  int i, j, k;
-  const int N = A.size();
-  for (i = 0; i < N; i++) {
-    // Find pivot
+  const int m      = (int)A.size();
+  const int n      = n_cols;
+  const int nwords = gf2_nwords(n);
+  cuopt_assert(m > 0, "");
+  cuopt_assert(n >= 0, "");
+  cuopt_assert((int)b.size() == m, "");
+  cuopt_assert((int)A[0].size() == nwords, "");
+
+  // pivot_row_of_col[c] = row holding the pivot for column c, or -1 if free
+  std::vector<int> pivot_row_of_col(n, -1);
+  int next_pivot_row = 0;
+
+  for (int col = 0; col < n; col++) {
     int pivot = -1;
-    for (j = i; j < N; j++) {
-      if (A[j][i]) {
-        pivot = j;
+    for (int r = next_pivot_row; r < m; r++) {
+      if (gf2_test_bit(A[r], col)) {
+        pivot = r;
         break;
       }
     }
-    if (pivot == -1) return false;  // No solution
-
-    // Swap current row with pivot row if needed
-    if (pivot != i) {
-      for (k = 0; k < N; k++) {
-        int temp    = A[i][k];
-        A[i][k]     = A[pivot][k];
-        A[pivot][k] = temp;
-      }
-      int temp = b[i];
-      b[i]     = b[pivot];
-      b[pivot] = temp;
+    if (pivot == -1) continue;  // free column
+
+    if (pivot != next_pivot_row) {
+      std::swap(A[next_pivot_row], A[pivot]);
+      std::swap(b[next_pivot_row], b[pivot]);
     }
 
-    // Eliminate downwards
-    for (j = i + 1; j < N; j++) {
-      if (A[j][i]) {
-        for (k = i; k < N; k++)
-          A[j][k] ^= A[i][k];
-        b[j] ^= b[i];
+    // Eliminate column from all other rows (RREF)
+    for (int r = 0; r < m; r++) {
+      if (r != next_pivot_row && gf2_test_bit(A[r], col)) {
+        for (int w = 0; w < nwords; w++)
+          A[r][w] ^= A[next_pivot_row][w];
+        b[r] ^= b[next_pivot_row];
       }
     }
+
+    pivot_row_of_col[col] = next_pivot_row;
+    next_pivot_row++;
   }
 
-  // Back-substitution
-  for (i = N - 1; i >= 0; i--) {
-    x[i] = b[i];
-    for (j = i + 1; j < N; j++)
-      x[i] ^= (A[i][j] & x[j]);
-    if (!A[i][i] && x[i]) return false;  // No solution
+  const int rank = next_pivot_row;
+  for (int r = rank; r < m; r++) {
+    for (int w = 0; w < nwords; w++) {
+      cuopt_assert(A[r][w] == 0, "RREF unused row must be zero");
+    }
+    if (b[r]) return gf2_status_t::Infeasible;
   }
-  return true;  // Success
+
+  std::vector<uint64_t> free_mask(nwords, 0);
+  for (int c = 0; c < n; c++) {
+    if (pivot_row_of_col[c] == -1) gf2_set_bit(free_mask, c);
+  }
+
+  determined.assign(n, 0);
+  x.assign(n, 0);
+
+  for (int col = 0; col < n; col++) {
+    int row = pivot_row_of_col[col];
+    if (row == -1) continue;  // free: x=0, determined=false
+
+    bool has_free_support = false;
+    for (int w = 0; w < nwords; w++) {
+      if (A[row][w] & free_mask[w]) {
+        has_free_support = true;
+        break;
+      }
+    }
+    // Particular solution with free vars = 0: x[pivot] = b[row]
+    x[col]          = b[row];
+    determined[col] = !has_free_support;
+  }
+
+  return gf2_status_t::Feasible;
 }
 
 template <typename f_t>
@@ -161,16 +211,20 @@ papilo::PresolveStatus GF2Presolve<f_t>::execute(const papilo::Problem<f_t>& pro
         if (key_var_idx != -1) { NOT_GF2("multiple key variables", var_idx); }
         key_var_idx   = var_idx;
         key_var_coeff = coeff;
-        gf2_key_vars.insert({var_idx, gf2_key_vars.size()});
       } else {
         // Binary variable
         constraint_bin_vars.push_back({var_idx, coeff});
-        gf2_bin_vars.insert({var_idx, gf2_bin_vars.size()});
       }
     }
 
     if (key_var_idx == -1) NOT_GF2("missing key variable");
 
+    // Commit to global maps only after the row is fully accepted
+    gf2_key_vars.insert({(size_t)key_var_idx, gf2_key_vars.size()});
+    for (auto [bin_var, _] : constraint_bin_vars) {
+      gf2_bin_vars.insert({bin_var, gf2_bin_vars.size()});
+    }
+
     gf2_constraints.emplace_back((size_t)cstr_idx,
                                  std::move(constraint_bin_vars),
                                  std::pair<size_t, f_t>{key_var_idx, key_var_coeff},
@@ -183,12 +237,11 @@ papilo::PresolveStatus GF2Presolve<f_t>::execute(const papilo::Problem<f_t>& pro
   // If no GF2 constraints found, return unchanged
   if (gf2_constraints.empty()) { return papilo::PresolveStatus::kUnchanged; }
 
-  // Skip if that would cause computational explosion (O(n^3) with simple gaussian elimination)
-  if (gf2_constraints.size() > 1000) { return papilo::PresolveStatus::kUnchanged; }
+  // one unique key per GF2 row. #bins may differ from #rows.
+  if (gf2_key_vars.size() != gf2_constraints.size()) { return papilo::PresolveStatus::kUnchanged; }
 
-  // Validate structure
-  if (gf2_key_vars.size() != gf2_constraints.size() ||
-      gf2_bin_vars.size() != gf2_constraints.size()) {
+  // Skip if that would cause computational explosion (dense GE ~ O(m * n * min(m,n)))
+  if (gf2_constraints.size() > 1000 || gf2_bin_vars.size() > 1000) {
     return papilo::PresolveStatus::kUnchanged;
   }
 
@@ -198,40 +251,76 @@ papilo::PresolveStatus GF2Presolve<f_t>::execute(const papilo::Problem<f_t>& pro
     gf2_bin_vars_invmap.insert({gf2_idx, var_idx});
   }
 
-  // Build binary matrix
-  // Could be a flat vector but. oh well. in practice N is small
-  std::vector<std::vector<int>> A(gf2_constraints.size(),
-                                  std::vector<int>(gf2_constraints.size(), 0));
-  std::vector<int> b(gf2_constraints.size());
-  for (size_t gf2_cstr_idx = 0; gf2_cstr_idx < gf2_constraints.size(); ++gf2_cstr_idx) {
+  // Build binary matrix as packed uint64_t words
+  const int m      = (int)gf2_constraints.size();
+  const int n      = (int)gf2_bin_vars.size();
+  const int nwords = gf2_nwords(n);
+  std::vector<std::vector<uint64_t>> A(m, std::vector<uint64_t>(nwords, 0));
+  std::vector<int> b(m);
+  for (int gf2_cstr_idx = 0; gf2_cstr_idx < m; ++gf2_cstr_idx) {
     const auto& cons = gf2_constraints[gf2_cstr_idx];
     for (auto [bin_var, _] : cons.bin_vars) {
-      A[gf2_cstr_idx][gf2_bin_vars[bin_var]] = 1;
+      gf2_set_bit(A[gf2_cstr_idx], (int)gf2_bin_vars[bin_var]);
     }
     b[gf2_cstr_idx] = cons.rhs;
   }
 
-  std::vector<int> solution(gf2_constraints.size());
-  bool feasible = gf2_solve(A, b, solution);
-  if (!feasible) { return papilo::PresolveStatus::kInfeasible; }
+  std::vector<int> solution(n);
+  std::vector<uint8_t> determined(n);
+  gf2_status_t gf2_status = gf2_solve(A, n, b, solution, determined);
+  if (gf2_status == gf2_status_t::Infeasible) { return papilo::PresolveStatus::kInfeasible; }
 
   std::unordered_map<size_t, f_t> fixings;
-  // Fix binary variables
-  for (size_t sol_idx = 0; sol_idx < gf2_constraints.size(); ++sol_idx) {
-    fixings[gf2_bin_vars_invmap[sol_idx]] = solution[sol_idx];
+
+  // Fix only uniquely determined binaries
+  for (int sol_idx = 0; sol_idx < n; ++sol_idx) {
+    if (determined[sol_idx]) { fixings[gf2_bin_vars_invmap[sol_idx]] = solution[sol_idx]; }
   }
 
-  // Compute fixings for key variables by solving for the constraint
+  // Fix key only when every binary in that constraint is uniquely determined
   for (const auto& cons : gf2_constraints) {
+    bool all_bins_determined = true;
+    for (auto [bin_var, _] : cons.bin_vars) {
+      cuopt_assert(gf2_bin_vars.count(bin_var), "");
+      if (!determined[gf2_bin_vars[bin_var]]) {
+        all_bins_determined = false;
+        break;
+      }
+    }
+    if (!all_bins_determined) continue;
+
     auto [key_var_idx, key_var_coeff] = cons.key_var;
-    f_t constraint_rhs                = lhs_values[cons.cstr_idx];  // equality constraint
+    const f_t constraint_rhs          = std::round(lhs_values[cons.cstr_idx]);
     f_t lhs                           = -constraint_rhs;
     for (auto [bin_var, coeff] : cons.bin_vars) {
+      cuopt_assert(fixings.count(bin_var), "");
       lhs += fixings[bin_var] * coeff;
     }
-    fixings[key_var_idx] = std::round(-lhs / key_var_coeff);
+    const f_t key_val = std::round(-lhs / key_var_coeff);
+
+    // Residual must be exactly 0 after rounding (rejects half-integer / inconsistent carry)
+    if (!num.isEq(lhs + key_val * key_var_coeff, f_t{0})) {
+      return papilo::PresolveStatus::kInfeasible;
+    }
+    // Dual-role: same var already fixed as a GF(2) binary
+    if (fixings.count(key_var_idx) && !num.isEq(fixings[key_var_idx], key_val)) {
+      return papilo::PresolveStatus::kInfeasible;
+    }
+    if (!col_flags[key_var_idx].test(papilo::ColFlag::kLbInf) &&
+        key_val < lower_bounds[key_var_idx] - integrality_tolerance) {
+      return papilo::PresolveStatus::kInfeasible;
+    }
+    if (!col_flags[key_var_idx].test(papilo::ColFlag::kUbInf) &&
+        key_val > upper_bounds[key_var_idx] + integrality_tolerance) {
+      return papilo::PresolveStatus::kInfeasible;
+    }
+
+    fixings[key_var_idx] = key_val;
   }
 
+  // necessary because Papilo asserts on empty TransactionGuard
+  if (fixings.empty()) { return papilo::PresolveStatus::kUnchanged; }
+
   papilo::PresolveStatus status = papilo::PresolveStatus::kUnchanged;
   papilo::TransactionGuard rg{reductions};
   for (const auto& [var_idx, fixing] : fixings) {

cpp/src/mip_heuristics/presolve/gf2_presolve.hpp

@@ -20,8 +20,22 @@
 #pragma GCC diagnostic pop
 #endif
 
+#include <cstdint>
+#include <vector>
+
 namespace cuopt::mathematical_optimization::mip {
 
+enum class gf2_status_t { Feasible, Infeasible };
+
+// Solves A x = b over GF(2). A is m x n, each row packed into ceil(n/64) words (column c lives at
+// word c/64, bit c%64). Trashes A and b. On Feasible, x is the solution obtained by setting the
+// free variables to 0, and determined[c] is set iff x[c] is the same in every solution.
+gf2_status_t gf2_solve(std::vector<std::vector<uint64_t>>& A,
+                       int n_cols,
+                       std::vector<int>& b,
+                       std::vector<int>& x,
+                       std::vector<uint8_t>& determined);
+
 template <typename f_t>
 class GF2Presolve : public papilo::PresolveMethod<f_t> {
  public: