fix(socp): canonicalize QC Q COO for rotated SOC detection (#1435)

cpp/include/cuopt/linear_programming/io/data_model_view.hpp

@@ -417,10 +417,14 @@ class data_model_view_t {
   bool is_Q_symmetrized() const noexcept;
 
   /**
-   * @brief Quadratic constraints (MPS QCMATRIX); owned copy for writers when not using spans.
+   * @brief Set quadratic constraints (linear part + Q in COO) on the view.
+   *
+   * Stores an owned copy retrievable via get_quadratic_constraints().
    */
   void set_quadratic_constraints(
     std::vector<typename mps_data_model_t<i_t, f_t>::quadratic_constraint_t> constraints);
+  /** @copydoc set_quadratic_constraints(std::vector<typename mps_data_model_t<i_t,
+   * f_t>::quadratic_constraint_t>) */
   template <typename qc_t>
   void set_quadratic_constraints(const std::vector<qc_t>& constraints)
   {

cpp/include/cuopt/linear_programming/io/mps_data_model.hpp

@@ -239,7 +239,9 @@ class mps_data_model_t {
    * - row identity and type (from ROWS),
    * - sparse linear coefficients (from COLUMNS),
    * - RHS value (from RHS),
-   * - quadratic matrix Q in COO (SoA: row, col, value) from QCMATRIX — one triplet per nonzero.
+   * - quadratic matrix Q in canonical COO (SoA: row, col, value): one triplet per variable pair,
+   *   sorted by (row, col). Off-diagonal cross terms store the full coefficient on x_i*x_j
+   *   (e.g. -2*d for rotated SOC), not as symmetric halves.
    */
   struct quadratic_constraint_t {
     /** ROWS declaration index (among all constraint rows), not an index into the linear CSR. */
@@ -251,7 +253,7 @@ class mps_data_model_t {
     std::vector<f_t> linear_values{};
     std::vector<i_t> linear_indices{};
     f_t rhs_value{f_t(0)};
-    /** Q nonzeros: parallel arrays, same length (COO / SoA). Sorted by (row, col) in append. */
+
     std::vector<i_t> rows{};
     std::vector<i_t> cols{};
     std::vector<f_t> vals{};
@@ -262,7 +264,11 @@ class mps_data_model_t {
    * @note All span inputs are host memory; the model copies this data.
    * @param linear_values, linear_indices Same nnz; can be empty for a purely quadratic row (rare).
    * @param vals, rows, cols COO triplets for Q; same length; may all be empty if Q is empty.
-   *        Stored sorted by (row, col).
+   *        Canonicalized on ingest to one triplet per variable pair (full x^T Q x coefficient).
+   * @param require_symmetric_q_offdiagonal Input validation only (default false). When false,
+   *        a single off-diagonal (i,j,v) per cross term is accepted (API/LP/C style). When true
+   *        (MPS QCMATRIX), each off-diagonal pair must appear as both (i,j,v) and (j,i,v) with
+   *        equal v before the halves merge to one stored entry; does not change canonical output.
    * @param constraint_row_type MPS ROWS type: 'L' (<=) or 'G' (>=). Stored as given; 'G' rows are
    *        converted to '<=' form when building the SOCP for the barrier solver. Equality ('E') is
    *        not supported.
@@ -275,7 +281,8 @@ class mps_data_model_t {
                                    f_t rhs_value,
                                    std::span<const f_t> vals,
                                    std::span<const i_t> rows,
-                                   std::span<const i_t> cols);
+                                   std::span<const i_t> cols,
+                                   bool require_symmetric_q_offdiagonal = false);
 
   const std::vector<quadratic_constraint_t>& get_quadratic_constraints() const;
 
@@ -385,4 +392,14 @@ class mps_data_model_t {
 
 };  // class mps_data_model_t
 
+/**
+ * @brief Canonicalize Q COO on each quadratic constraint in place.
+ *
+ * Used for raw API / data_model_view input (single cross term per pair). MPS QCMATRIX
+ * symmetric-half validation is done in append_quadratic_constraint instead.
+ */
+template <typename i_t, typename f_t>
+void canonicalize_quadratic_constraints(
+  std::vector<typename mps_data_model_t<i_t, f_t>::quadratic_constraint_t>& constraints);
+
 }  // namespace cuopt::linear_programming::io

cpp/include/cuopt/linear_programming/optimization_problem_interface.hpp

@@ -84,6 +84,7 @@ class optimization_problem_interface_t {
    * Quadratic matrix Q is COO (row_index, col_index, coeff spans). Linear term d uses parallel
    * linear_values and linear_indices (empty allowed). constraint_row_index is assigned
    * automatically as n_linear_constraints + n_existing_quadratic_constraints.
+   * Q is canonicalized on ingest (merge duplicates, collapse matching symmetric pairs).
    */
   virtual void add_quadratic_constraint(char constraint_row_type,
                                         f_t rhs_value,

cpp/include/cuopt/linear_programming/optimization_problem_utils.hpp

@@ -134,7 +134,7 @@ void populate_from_mps_data_model(optimization_problem_interface_t<i_t, f_t>* pr
                                             q_offsets.data(),
                                             n_vars + 1);
   }
-  // Handle quadratic constraints if present
+  // Quadratic constraints from mps_data_model are already canonical (append_quadratic_constraint).
   if (data_model.has_quadratic_constraints()) {
     problem->set_quadratic_constraints(data_model.get_quadratic_constraints());
   }
@@ -295,8 +295,11 @@ void populate_from_data_model_view(
     problem->set_row_names(data_model->get_row_names());
   }
 
+  // Raw Q COO from data_model_view is canonicalized here before solver storage.
   if (data_model->has_quadratic_constraints()) {
-    problem->set_quadratic_constraints(data_model->get_quadratic_constraints());
+    auto qcs = data_model->get_quadratic_constraints();
+    io::canonicalize_quadratic_constraints<i_t, f_t>(qcs);
+    problem->set_quadratic_constraints(std::move(qcs));
   }
 }
 

cpp/src/barrier/translate_soc.hpp

@@ -92,8 +92,7 @@ void convert_quadratic_constraints_to_second_order_cones(
   //        -s*x_head^2 + sum_i s*x_tail_i^2 <= 0   (any common s > 0; divide by s to normalize)
   //   2) rotated SOC rows:
   //        -2*d*x_head0*x_head1 + sum_i s*x_tail_i^2 <= 0   (d>0, s>0; canonical d=s)
-  //      symmetric Q off-diagonals (-d,-d) give x^T Q x cross term -2*d*x0*x1, i.e. a*x0*x1
-  //      in the inequality 2*d*x0*x1 >= s*||tail||^2 with a = 2*d. Lift uses sqrt(d/s) on heads.
+  //      stored as one Q cross term (head0, head1, -2*d). Lift uses sqrt(d/s) on heads.
   //   3) quadratic rows with linear part:
   //        sum_i s*x_tail_i^2 + a^T x <= 0
   //      represented as diagonal +s QCMATRIX entries plus linear terms in COLUMNS.
@@ -108,7 +107,8 @@ void convert_quadratic_constraints_to_second_order_cones(
     i_t head1{};
     std::vector<i_t> tails{};
     bool head1_is_constant_half{false};
-    /// For two-head rotated SOC: sqrt(d/s) where Q_off = -d and tail diagonals +s (canonical 1).
+    /// For two-head rotated SOC: sqrt(d/s) where Q cross = -2*d and tail diagonals +s (canonical
+    /// 1).
     f_t head_lift_sqrt_ratio{1};
   };
   // This is the index of the auxiliary variable for the linear part of the quadratic constraint.
@@ -188,7 +188,7 @@ void convert_quadratic_constraints_to_second_order_cones(
 
     // Verify Q as either:
     // - standard SOC: one diagonal -s (head), tail diagonals +s for a common s > 0,
-    // - rotated SOC: symmetric (-s,-s) off-diagonal pair on the two heads, tails +s,
+    // - rotated SOC: one off-diagonal cross term (-2*d) on the two heads, tails +s,
     // - affine SOC: tail diagonals +s and linear terms (no Q off-diagonals).
     // Feasibility is unchanged after dividing the quadratic row by s; affine rows also scale
     // linear coefficients when forming the auxiliary t = -(1/s) a^T x.
@@ -274,11 +274,21 @@ void convert_quadratic_constraints_to_second_order_cones(
         }
       }
     }
-    const bool use_general_path = has_duplicate_rows || has_near_zero_diag || has_nonzero_rhs ||
-                                  has_nonuniform_diag || offdiag_entries.size() > 2 ||
-                                  (offdiag_entries.size() == 1) || (neg_diag_rows.size() > 1) ||
-                                  (!neg_diag_rows.empty() && has_linear_part) ||
-                                  (!neg_diag_rows.empty() && !offdiag_entries.empty());
+    const bool rotated_soc_cross_eligible = [&]() {
+      if (offdiag_entries.size() != 1 || has_linear_part || !neg_diag_rows.empty()) {
+        return false;
+      }
+      const i_t a = std::get<0>(offdiag_entries[0]);
+      const i_t b = std::get<1>(offdiag_entries[0]);
+      const f_t v = std::get<2>(offdiag_entries[0]);
+      // Match cross_d = -v/2 > tol validated on the rotated SOC fast path below.
+      return a != b && (-v / f_t(2)) > tol;
+    }();
+    const bool use_general_path =
+      has_duplicate_rows || has_near_zero_diag || has_nonzero_rhs || has_nonuniform_diag ||
+      offdiag_entries.size() > 1 || (offdiag_entries.size() == 1 && !rotated_soc_cross_eligible) ||
+      (neg_diag_rows.size() > 1) || (!neg_diag_rows.empty() && has_linear_part) ||
+      (!neg_diag_rows.empty() && !offdiag_entries.empty());
 
     f_t uniform_s        = 0;
     bool have_uniform_s  = false;
@@ -456,52 +466,28 @@ void convert_quadratic_constraints_to_second_order_cones(
                       "Quadratic constraint '%s' rotated SOC Q: could not infer uniform scale s",
                       qc.constraint_row_name.c_str());
         cuopt_expects(
-          offdiag_entries.size() == 2,
+          offdiag_entries.size() == 1,
           error_type_t::ValidationError,
-          "Quadratic constraint '%s' rotated SOC Q must contain exactly one symmetric off-diagonal "
-          "pair (-d,-d); found %zu off-diagonal entries",
+          "Quadratic constraint '%s' rotated SOC Q must contain exactly one cross term with "
+          "coefficient -2*d on the head variable pair; found %zu off-diagonal entries",
           qc.constraint_row_name.c_str(),
           offdiag_entries.size());
 
         const i_t a  = std::get<0>(offdiag_entries[0]);
         const i_t b  = std::get<1>(offdiag_entries[0]);
         const f_t v0 = std::get<2>(offdiag_entries[0]);
-        cuopt_expects(
-          v0 < -tol,
-          error_type_t::ValidationError,
-          "Quadratic constraint '%s' rotated SOC Q off-diagonal must be negative; got %.17g",
-          qc.constraint_row_name.c_str(),
-          static_cast<double>(v0));
+
         cuopt_expects(a != b,
                       error_type_t::ValidationError,
-                      "Quadratic constraint '%s' rotated SOC Q off-diagonal pair must use distinct "
+                      "Quadratic constraint '%s' rotated SOC Q cross term must use distinct head "
                       "variables",
                       qc.constraint_row_name.c_str());
-        cuopt_expects(std::get<0>(offdiag_entries[1]) == b && std::get<1>(offdiag_entries[1]) == a,
-                      error_type_t::ValidationError,
-                      "Quadratic constraint '%s' rotated SOC Q must have symmetric entries (a,b) "
-                      "and (b,a) with the same value",
-                      qc.constraint_row_name.c_str());
-        const f_t v1 = std::get<2>(offdiag_entries[1]);
-        cuopt_expects(
-          v1 < -tol,
-          error_type_t::ValidationError,
-          "Quadratic constraint '%s' rotated SOC Q off-diagonal must be negative; got %.17g",
-          qc.constraint_row_name.c_str(),
-          static_cast<double>(v1));
-        cuopt_expects(
-          approx_eq_scaled(v0, v1),
-          error_type_t::ValidationError,
-          "Quadratic constraint '%s' rotated SOC Q symmetric off-diagonals must match; got %.17g "
-          "and %.17g",
-          qc.constraint_row_name.c_str(),
-          static_cast<double>(v0),
-          static_cast<double>(v1));
-        const f_t cross_d = -v0;
+        const f_t cross_d = -v0 / f_t(2);
         cuopt_expects(
           cross_d > tol,
           error_type_t::ValidationError,
-          "Quadratic constraint '%s' rotated SOC Q cross coefficient d = -Q_off must be positive",
+          "Quadratic constraint '%s' rotated SOC Q cross coefficient d = -Q_cross/2 must be "
+          "positive",
           qc.constraint_row_name.c_str());
         const f_t head_lift_sqrt_ratio = std::sqrt(cross_d / uniform_s);
         cuopt_expects(std::isfinite(static_cast<double>(head_lift_sqrt_ratio)),
@@ -511,12 +497,12 @@ void convert_quadratic_constraints_to_second_order_cones(
                       qc.constraint_row_name.c_str(),
                       static_cast<double>(cross_d),
                       static_cast<double>(uniform_s));
-        cuopt_expects(static_cast<i_t>(tail_vars.size()) == q_nnz - 2,
+        cuopt_expects(static_cast<i_t>(tail_vars.size()) == q_nnz - 1,
                       error_type_t::ValidationError,
                       "Quadratic constraint '%s' rotated SOC Q: expected %d diagonal +s entries "
                       "(tails), found %zu",
                       qc.constraint_row_name.c_str(),
-                      static_cast<int>(q_nnz - 2),
+                      static_cast<int>(q_nnz - 1),
                       tail_vars.size());
         cuopt_expects(q_nnz >= 3,
                       error_type_t::ValidationError,

cpp/src/grpc/codegen/generate_conversions.py

@@ -2215,6 +2215,10 @@ def _gen_proto_to_problem(registry, indent="  "):
                 f'{b["name"]} size mismatch");'
             )
             lines.append(f"{ind}    }}")
+        if name == "quadratic_constraints":
+            lines.append(
+                f"{ind}    io::canonicalize_coo_matrix(_entry.rows, _entry.cols, _entry.vals);"
+            )
         lines.append(f"{ind}    _entries.push_back(std::move(_entry));")
         lines.append(f"{ind}  }}")
         lines.append(f"{ind}  cpu_problem.{setter}(std::move(_entries));")
@@ -2811,6 +2815,10 @@ def _gen_chunked_arrays_to_problem(registry, indent="  "):
                     f'{b["name"]} size mismatch");'
                 )
                 lines.append(f"{ind}    }}")
+            if name == "quadratic_constraints":
+                lines.append(
+                    f"{ind}    io::canonicalize_coo_matrix(_entry.rows, _entry.cols, _entry.vals);"
+                )
             lines.append(f"{ind}    _entries.push_back(std::move(_entry));")
             lines.append(f"{ind}  }}")
             lines.append(f"{ind}  cpu_problem.{setter}(std::move(_entries));")

cpp/src/grpc/codegen/generated/generated_chunked_arrays_to_problem.inc

@@ -194,6 +194,7 @@
       if (_entry.cols.size() != _entry.vals.size()) {
         throw std::invalid_argument("set_quadratic_constraints: cols/vals size mismatch");
       }
+      io::canonicalize_coo_matrix(_entry.rows, _entry.cols, _entry.vals);
       _entries.push_back(std::move(_entry));
     }
     cpu_problem.set_quadratic_constraints(std::move(_entries));

cpp/src/grpc/codegen/generated/generated_proto_to_problem.inc

@@ -96,6 +96,7 @@
       if (_entry.cols.size() != _entry.vals.size()) {
         throw std::invalid_argument("set_quadratic_constraints: cols/vals size mismatch");
       }
+      io::canonicalize_coo_matrix(_entry.rows, _entry.cols, _entry.vals);
       _entries.push_back(std::move(_entry));
     }
     cpu_problem.set_quadratic_constraints(std::move(_entries));

cpp/src/grpc/grpc_problem_mapper.cpp

@@ -11,6 +11,7 @@
 #include <cuopt/linear_programming/cpu_optimization_problem.hpp>
 #include <cuopt/linear_programming/mip/solver_settings.hpp>
 #include <cuopt/linear_programming/pdlp/solver_settings.hpp>
+#include <quadratic_constraint_coo.hpp>
 #include "grpc_settings_mapper.hpp"
 
 #include <algorithm>

cpp/src/io/CMakeLists.txt

@@ -10,6 +10,7 @@ set(PARSERS_SRC_FILES
   ${CMAKE_CURRENT_SOURCE_DIR}/mps_data_model.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/mps_parser.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/mps_writer.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/quadratic_constraint_coo.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/parser.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/writer.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/utilities/cython_parser.cpp

cpp/src/io/data_model_view.cpp

@@ -6,6 +6,7 @@
 /* clang-format on */
 
 #include <cuopt/linear_programming/io/data_model_view.hpp>
+
 #include <utilities/error.hpp>
 
 #include <span>

cpp/src/io/lp_parser.cpp

@@ -118,43 +118,6 @@ bool is_free_keyword(std::string_view lower) { return lower == "free"; }
 
 bool is_infinity_text(std::string_view lower) { return lower == "inf" || lower == "infinity"; }
 
-// Builds the symmetric Q in COO from LP-format raw upper-triangular triples.
-// Each input triple (i, j, c) with i <= j represents `c * x_i * x_j` in the
-// LP source. The output Q satisfies x^T Q x = sum of those terms.
-//   Diagonal (i == j): Q[i,i] = c (one entry).
-//   Off-diagonal (i != j): Q[i,j] = Q[j,i] = c/2 (two entries; symmetric split).
-template <typename i_t, typename f_t>
-void build_symmetric_q_coo(const coo_entries_t<i_t, f_t>& raw_triples,
-                           std::vector<i_t>& out_row_indices,
-                           std::vector<i_t>& out_col_indices,
-                           std::vector<f_t>& out_values)
-{
-  out_row_indices.clear();
-  out_col_indices.clear();
-  out_values.clear();
-  out_row_indices.reserve(raw_triples.size() * 2);
-  out_col_indices.reserve(raw_triples.size() * 2);
-  out_values.reserve(raw_triples.size() * 2);
-
-  for (size_t p = 0; p < raw_triples.size(); p++) {
-    const i_t i = raw_triples.rows[p];
-    const i_t j = raw_triples.cols[p];
-    const f_t c = raw_triples.vals[p];
-    if (i == j) {
-      out_row_indices.push_back(i);
-      out_col_indices.push_back(i);
-      out_values.push_back(c);
-    } else {
-      out_row_indices.push_back(i);
-      out_col_indices.push_back(j);
-      out_values.push_back(c / 2);
-      out_row_indices.push_back(j);
-      out_col_indices.push_back(i);
-      out_values.push_back(c / 2);
-    }
-  }
-}
-
 // ===========================================================================
 // Token stream
 // ===========================================================================
@@ -964,9 +927,8 @@ void LpParseEngine<i_t, f_t>::parse_quadratic_bracket(int outer_sign,
                        "constraint must contain at least one quadratic term",
                        peek().line);
 
-    // Coefficients are at face value — the post-pass that flushes the
-    // quadratic_constraint_block_t to the data model handles the symmetric
-    // expansion and the /2 split for off-diagonals.
+    // Coefficients are face value for x^T Q x; canonicalization runs in
+    // append_quadratic_constraint() when the block is flushed to the data model.
     out_quad_entries.rows.insert(
       out_quad_entries.rows.end(), raw_quad.rows.begin(), raw_quad.rows.end());
     out_quad_entries.cols.insert(
@@ -1514,19 +1476,15 @@ void flush_quadratic_constraints(mps_data_model_t<i_t, f_t>& problem,
   const i_t linear_row_count = static_cast<i_t>(parser.row_names.size());
   for (i_t k = 0; k < static_cast<i_t>(parser.quadratic_constraint_blocks.size()); k++) {
     const auto& block = parser.quadratic_constraint_blocks[k];
-    std::vector<i_t> q_row_indices;
-    std::vector<i_t> q_col_indices;
-    std::vector<f_t> q_values;
-    build_symmetric_q_coo(block.quad_triples, q_row_indices, q_col_indices, q_values);
     problem.append_quadratic_constraint(linear_row_count + k,
                                         block.row_name,
                                         static_cast<char>(block.row_type),
                                         block.linear_values,
                                         block.linear_indices,
                                         block.rhs_value,
-                                        q_values,
-                                        q_row_indices,
-                                        q_col_indices);
+                                        block.quad_triples.vals,
+                                        block.quad_triples.rows,
+                                        block.quad_triples.cols);
   }
 }