Fix the bug detecting nonconvex quadratic constraints in the general … (#1438)

Cross-only indefinite quadratic constraints (e.g. `2·x₀·x₁ ≤ 0.5` with `H = [[0,2],[2,0]]`) previously passed convexity checks on the general SOC conversion path. Diagonal LDLT returned `rank = 0` without `INDEFINITE_MATRIX_RETURN`, and a degenerate `r = 0` SOC lift silently dropped the quadratic term, allowing wrong solutions to be reported as Optimal (e.g. obj=2 at (1,1) instead of rejecting the problem).

Fixes:
- `right_looking_ldlt`: return `INDEFINITE_MATRIX_RETURN` when factorization stalls at `rank = 0` on a nonzero matrix.
- `translate_soc.hpp` (general path): `cuopt_expects(rank >= 1, …)` before building the SOC lift.

Adds LDLT and SOC-conversion unit tests for the cross-only indefinite case.

## Issue

closes #1434

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

Approvers:
  - Rajesh Gandham (https://github.com/rg20)
  - Ramakrishna Prabhu (https://github.com/ramakrishnap-nv)
  - Miles Lubin (https://github.com/mlubin)

URL: #1438

Author: yuwenchen95

.gitignore

@@ -12,6 +12,7 @@ __pycache__
 *.dylib
 .cache
 .vscode
+.cursor/
 *.swp
 *.pytest_cache
 *.manifest

cpp/src/barrier/translate_soc.hpp

@@ -659,9 +659,15 @@ void convert_quadratic_constraints_to_second_order_cones(
                     "Quadratic constraint '%s' is non-convex (Q matrix is indefinite)",
                     qc.constraint_row_name.c_str());
 
-      // Since q_nnz >= 1 is enforced above, Q is nonzero and rank must be >= 1.
-      // (A nonzero Q entry produces a nonzero H diagonal or off-diagonal, guaranteeing rank > 0.)
-      assert(rank >= 1);
+      // q_nnz >= 1 implies H is nonzero, but diagonal LDLT may still return rank 0 (e.g. cross-only
+      // indefinite H with zero diagonals). Reject before building a degenerate r=0 SOC lift.
+      cuopt_expects(rank >= 1,
+                    error_type_t::ValidationError,
+                    "Quadratic constraint '%s' is non-convex or could not be converted to a "
+                    "second-order cone (LDLT rank %d; Q may be indefinite or have zero diagonal "
+                    "with cross terms)",
+                    qc.constraint_row_name.c_str(),
+                    static_cast<int>(rank));
 
       // Step 4: Build standard SOC of dimension rank + 2.
       // New variables: y_0,...,y_{r-1}, s_0 (head), s_{r+1} (tail)

cpp/src/dual_simplex/right_looking_lu.cpp

@@ -13,7 +13,6 @@
 
 #include <cassert>
 #include <cmath>
-#include <cstdio>
 
 using cuopt::ins_vector;
 
@@ -1684,7 +1683,8 @@ i_t right_looking_ldlt(const csc_matrix_t<i_t, f_t>& A,
                        f_t& work_estimate)
 {
   raft::common::nvtx::range scope("LU::right_looking_ldlt");
-  const i_t n = A.n;
+  const i_t n         = A.n;
+  const i_t input_nnz = A.nnz();
   assert(A.m == n);
 
   symmetric_trailing_matrix_t<i_t, f_t> trailing_matrix(A);
@@ -1720,7 +1720,12 @@ i_t right_looking_ldlt(const csc_matrix_t<i_t, f_t>& A,
     f_t pivot_val = 0;
     trailing_matrix.symmetric_markowitz_search(pivot_tol, pivot_p, pivot_val);
 
-    if (pivot_p == -1) { break; }  // No acceptable pivot found (remaining diagonals are zero/tiny)
+    if (pivot_p == -1) {
+      // No acceptable diagonal pivot. Zero matrix is rank 0; nonzero cross-only indefinite
+      // matrices (e.g. [[0, a]; [a, 0]]) also stall here and must not be treated as PSD.
+      if (pivots == 0 && input_nnz > 0) { return INDEFINITE_MATRIX_RETURN; }
+      break;
+    }
 
     // Check for indefiniteness: a negative pivot means the matrix is not PSD
     if (pivot_val < 0) { return INDEFINITE_MATRIX_RETURN; }

cpp/src/dual_simplex/right_looking_lu.hpp

@@ -43,7 +43,8 @@ i_t right_looking_lu_row_permutation_only(const csc_matrix_t<i_t, f_t>& A,
 //   D = diagonal factor (length = rank)
 // Returns:
 //   rank >= 0: number of successful pivots with D(k,k) >= pivot_tol (PSD case).
-//   INDEFINITE_MATRIX_RETURN (-4): a negative pivot was encountered (matrix is not PSD).
+//   INDEFINITE_MATRIX_RETURN (-4): a negative pivot was encountered, or the matrix is nonzero
+//     but no acceptable diagonal pivot was found (matrix is not PSD).
 //   CONCURRENT_HALT_RETURN (-2): concurrent halt requested.
 //   TIME_LIMIT_RETURN (-3): time limit exceeded.
 template <typename i_t, typename f_t>

cpp/src/pdlp/solve.cu

@@ -2193,6 +2193,13 @@ mps_data_model_to_optimization_problem(
                                               Q_offsets.size());
   }
 
+  // Preserve quadratic constraints.
+  if (data_model.has_quadratic_constraints()) {
+    static_cast<cuopt::mathematical_optimization::optimization_problem_interface_t<i_t, f_t>&>(
+      op_problem)
+      .set_quadratic_constraints(data_model.get_quadratic_constraints());
+  }
+
   return op_problem;
 }
 

cpp/src/pdlp/translate.hpp

@@ -44,8 +44,10 @@ static simplex::user_problem_t<i_t, f_t> cuopt_problem_to_user_problem(
   csr_A.x         = std::move(A_values);
   csr_A.j         = std::move(A_indices);
   csr_A.row_start = std::move(A_offsets);
-
-  csr_A.to_compressed_col(user_problem.A);
+  if (m == 0) {
+    csr_A.row_start.resize(1);
+    csr_A.row_start[0] = 0;
+  }
 
   user_problem.rhs.resize(m);
   user_problem.row_sense.resize(m);
@@ -96,6 +98,13 @@ static simplex::user_problem_t<i_t, f_t> cuopt_problem_to_user_problem(
   user_problem.Q_indices = problem.get_quadratic_objective_indices();
   user_problem.Q_values  = problem.get_quadratic_objective_values();
 
+  if (problem.has_quadratic_constraints()) {
+    barrier::convert_quadratic_constraints_to_second_order_cones<i_t, f_t>(
+      n, problem.get_quadratic_constraints(), csr_A, user_problem);
+  }
+
+  csr_A.to_compressed_col(user_problem.A);
+
   return user_problem;
 }
 

cpp/tests/dual_simplex/unit_tests/right_looking_ldlt.cpp

@@ -512,6 +512,26 @@ TEST(right_looking_ldlt, indefinite_2x2)
   EXPECT_EQ(result, INDEFINITE_MATRIX_RETURN);
 }
 
+// Test 14b: Cross-only indefinite matrix (zero diagonals, no 1x1 pivot).
+// A = [0, 2; 2, 0] has eigenvalues +2 and -2 (indefinite).
+TEST(right_looking_ldlt, indefinite_cross_only_2x2)
+{
+  const int n               = 2;
+  std::vector<double> dense = {0.0, 2.0, 2.0, 0.0};
+  auto A                    = dense_to_lower_csc(n, dense);
+
+  simplex_solver_settings_t<int, double> settings;
+  std::vector<int> perm;
+  csc_matrix_t<int, double> L(n, n, 1);
+  std::vector<double> D;
+  double work_estimate = 0;
+  double start_time    = tic();
+
+  int result = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate);
+
+  EXPECT_EQ(result, INDEFINITE_MATRIX_RETURN);
+}
+
 // Test 15: Larger indefinite matrix (4x4 with mixed eigenvalues).
 // A = [2, 1, 0, 1; 1, 0, 1, 0; 0, 1, 2, -1; 1, 0, -1, -1]
 TEST(right_looking_ldlt, indefinite_4x4)

cpp/tests/linear_programming/unit_tests/solution_interface_test.cu

@@ -400,6 +400,30 @@ TEST_F(SolutionInterfaceTest, mps_data_model_to_optimization_problem)
   }
 }
 
+TEST_F(SolutionInterfaceTest, mps_data_model_to_optimization_problem_quadratic_constraints)
+{
+  auto mps_data = cuopt::mathematical_optimization::io::read_lp_from_string<int, double>(R"LP(
+Minimize
+  obj: x + y
+Subject To
+  q0: [ 4 x * y ] <= 0.5
+Bounds
+  -1 <= x <= 1
+  -1 <= y <= 1
+End
+)LP");
+  ASSERT_TRUE(mps_data.has_quadratic_constraints());
+  ASSERT_EQ(mps_data.get_quadratic_constraints().size(), 1u);
+
+  raft::handle_t handle;
+  auto problem = mps_data_model_to_optimization_problem(&handle, mps_data);
+
+  ASSERT_TRUE(problem.has_quadratic_constraints());
+  ASSERT_EQ(problem.get_quadratic_constraints().size(), 1u);
+  EXPECT_EQ(problem.get_quadratic_constraints()[0].constraint_row_name, "q0");
+  EXPECT_NEAR(problem.get_quadratic_constraints()[0].rhs_value, 0.5, 1e-9);
+}
+
 // =============================================================================
 // Solution conversion tests (hand-constructed, known values)
 // =============================================================================

cpp/tests/socp/general_quadratic_test.cu

@@ -5,9 +5,12 @@
  */
 /* clang-format on */
 
+#include <gmock/gmock.h>
 #include <gtest/gtest.h>
 
 #include <barrier/translate_soc.hpp>
+#include <cuopt/error.hpp>
+#include <cuopt/mathematical_optimization/io/parser.hpp>
 #include <cuopt/mathematical_optimization/optimization_problem_interface.hpp>
 #include <dual_simplex/solve.hpp>
 #include <dual_simplex/sparse_matrix.hpp>
@@ -214,6 +217,90 @@ TEST(general_quadratic, rejects_non_convex)
     cuopt::logic_error);
 }
 
+// End-to-end: cross-only indefinite Q (issue #1434). H = [[0, 2]; [2, 0]] has zero diagonals so
+// LDLT returns rank 0 without a negative pivot; must still be rejected as non-convex via
+// solve_qcqp.
+TEST(general_quadratic, rejects_cross_only_indefinite)
+{
+  raft::handle_t handle{};
+  init_handler(&handle);
+
+  // H = [[0, 2]; [2, 0]] from [ 4 x * y ] in LP bracket notation.
+  auto lp = io::read_lp_from_string<i_t, f_t>(R"LP(
+Minimize
+  obj: x + y
+Subject To
+  q0: [ 4 x * y ] <= 0.5
+Bounds
+  -1 <= x <= 1
+  -1 <= y <= 1
+End
+)LP");
+
+  ASSERT_TRUE(lp.has_quadratic_constraints());
+  ASSERT_EQ(lp.get_quadratic_constraints().size(), 1u);
+
+  const i_t n = lp.get_n_variables();
+  const i_t m = lp.get_n_constraints();
+  EXPECT_EQ(n, 2);
+  EXPECT_EQ(m, 0);
+
+  user_problem_t<i_t, f_t> user_problem(&handle);
+  user_problem.num_rows  = m;
+  user_problem.num_cols  = n;
+  user_problem.objective = lp.get_objective_coefficients();
+
+  // Initialize the empty A matrix
+  user_problem.A.m      = m;
+  user_problem.A.n      = n;
+  user_problem.A.nz_max = 0;
+  user_problem.A.reallocate(0);
+  user_problem.A.col_start.assign(n + 1, 0);
+
+  user_problem.rhs.clear();
+  user_problem.row_sense.clear();
+  user_problem.lower          = lp.get_variable_lower_bounds();
+  user_problem.upper          = lp.get_variable_upper_bounds();
+  user_problem.num_range_rows = 0;
+  user_problem.var_types.assign(n, variable_type_t::CONTINUOUS);
+
+  const auto& src_qc = lp.get_quadratic_constraints()[0];
+  qc_t qc;
+  qc.constraint_row_index = src_qc.constraint_row_index;
+  qc.constraint_row_name  = src_qc.constraint_row_name;
+  qc.constraint_row_type  = src_qc.constraint_row_type;
+  qc.linear_values        = src_qc.linear_values;
+  qc.linear_indices       = src_qc.linear_indices;
+  qc.rhs_value            = src_qc.rhs_value;
+  qc.rows                 = src_qc.rows;
+  qc.cols                 = src_qc.cols;
+  qc.vals                 = src_qc.vals;
+
+  csr_matrix_t<i_t, f_t> csr_A(m, n, 0);
+  csr_A.m         = m;
+  csr_A.n         = n;
+  csr_A.row_start = {0};
+
+  std::vector<qc_t> qcs = {qc};
+
+  simplex_solver_settings_t<i_t, f_t> settings;
+  settings.barrier          = true;
+  settings.barrier_presolve = true;
+  settings.dualize          = 0;
+
+  try {
+    convert_quadratic_constraints_to_second_order_cones<i_t, f_t>(n, qcs, csr_A, user_problem);
+    csr_A.to_compressed_col(user_problem.A);
+    lp_solution_t<i_t, f_t> solution(user_problem.num_rows, user_problem.num_cols);
+    (void)solve_linear_program_with_barrier(user_problem, settings, solution);
+    FAIL() << "Expected ValidationError for cross-only indefinite Q";
+  } catch (const cuopt::logic_error& e) {
+    EXPECT_EQ(e.get_error_type(), cuopt::error_type_t::ValidationError);
+    EXPECT_THAT(e.what(), testing::HasSubstr("non-convex"));
+    EXPECT_THAT(e.what(), testing::HasSubstr("q0"));
+  }
+}
+
 // Test: rank-deficient PSD Q (e.g., Q = v*v^T with v = [1, 1])
 // minimize x0 + x1
 // subject to (x0 + x1)^2 <= 4   (i.e., |x0 + x1| <= 2)