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

Quadratic-constraint `Q` is now stored in **canonical COO** internally: one coefficient per variable pair (e.g. a single `-2` for `2·x₀·x₁`), with symmetric MPS halves merged at ingest. Canonicalization runs at ingest boundaries (MPS/LP parser, C API, Python→solver, gRPC, PDLP CPU/GPU problem setup). MPS `QCMATRIX` still accepts symmetric halves; the MPS writer expands canonical cross terms back to symmetric form on export.

The RSOC fast path now accepts a single eligible cross term (e.g. `-2·t·u` for `||tail||² ≤ 2·t·u`) instead of requiring two symmetric off-diagonal entries. Previously, natural Python/API forms were misrouted to the general path or failed pattern matching, producing wrong optima (e.g. ~0 instead of √2).

## Issue

closes #1435, #1354

Authors:
  - Yuwen Chen (https://github.com/yuwenchen95)

Approvers:
  - Ramakrishna Prabhu (https://github.com/ramakrishnap-nv)
  - Miles Lubin (https://github.com/mlubin)
  - Chris Maes (https://github.com/chris-maes)

URL: #1439

Author: yuwenchen95

cpp/include/cuopt/mathematical_optimization/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/mathematical_optimization/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 upper-triangular COO format.
+   *   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,8 @@ 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),
+   *        stored in upper-triangular (min row, max col) form.
    * @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.
@@ -385,4 +388,17 @@ class mps_data_model_t {
 
 };  // class mps_data_model_t
 
+/**
+ * @brief Canonicalize each quadratic constraint's Q matrix.
+ *
+ * Merge duplicate COO entries and off-diagonal variable pairs, then store Q in
+ * upper-triangular form with one coefficient per variable pair (the full x^T Q x
+ * contribution for that pair).
+ *
+ * @param constraints Quadratic constraints whose rows, cols, and vals are updated.
+ */
+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::mathematical_optimization::io

cpp/include/cuopt/mathematical_optimization/io/parser.hpp

@@ -75,8 +75,8 @@ mps_data_model_t<i_t, f_t> read_mps_from_string(std::string_view mps_contents,
  * between objective and constraints:
  *   - Objective bracket: MUST be followed by `/ 2` (the LP file states
  *     coefficients in the `0.5 x^T Q x` convention).
- *   - Constraint bracket: MUST NOT be followed by `/ 2`; coefficients are
- *     taken at face value (`x^T Q x`).
+ *   - Constraint bracket: MUST NOT be followed by `/ 2`; Q coefficients are
+ *     stored as upper triangular form.
  *
  * SOS constraints, PWL objectives, general constraints, and user cuts cause
  * a ValidationError when encountered.

cpp/include/cuopt/mathematical_optimization/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/mathematical_optimization/optimization_problem_utils.hpp

@@ -130,7 +130,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());
   }
@@ -291,8 +291,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,14 +497,14 @@ 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,
+        cuopt_expects(q_nnz >= 2,
                       error_type_t::ValidationError,
                       "Quadratic constraint '%s' rotated SOC Q must have at least 1 tail entry",
                       qc.constraint_row_name.c_str());

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/mathematical_optimization/cpu_optimization_problem.hpp>
 #include <cuopt/mathematical_optimization/mip/solver_settings.hpp>
 #include <cuopt/mathematical_optimization/pdlp/solver_settings.hpp>
+#include <mps_parser_internal.hpp>
 #include "grpc_settings_mapper.hpp"
 
 #include <algorithm>