Bug fix of dual init point computation in barrier

cpp/include/cuopt/mathematical_optimization/constants.h

@@ -194,6 +194,11 @@
 #define CUOPT_METHOD_BARRIER      3
 #define CUOPT_METHOD_UNSET        4
 
+#define CUOPT_BARRIER_DUAL_INITIAL_POINT_AUTOMATIC             -1
+#define CUOPT_BARRIER_DUAL_INITIAL_POINT_LUSTIG_MARSTEN_SHANNO 0
+#define CUOPT_BARRIER_DUAL_INITIAL_POINT_DUAL_LEAST_SQUARES    1
+#define CUOPT_BARRIER_DUAL_INITIAL_POINT_SEDUMI_MU             2
+
 /* @brief PDLP precision mode constants */
 #define CUOPT_PDLP_DEFAULT_PRECISION -1
 #define CUOPT_PDLP_SINGLE_PRECISION  0

cpp/src/barrier/barrier.cu

@@ -94,6 +94,39 @@ bool validate_barrier_cone_layout(const lp_problem_t<i_t, f_t>& problem,
   return true;
 }
 
+// Push entries into interior of nonnegative orthant and SOC.
+template <typename i_t, typename f_t>
+static void ensure_initial_point_interior(dense_vector_t<i_t, f_t>& values,
+                                          f_t epsilon_adjust,
+                                          const std::vector<i_t>& linear_mask,
+                                          i_t linear_end,
+                                          const std::vector<i_t>& cone_dims)
+{
+  f_t min_linear = inf;
+  for (i_t j = 0; j < linear_end; ++j) {
+    if (linear_mask[j]) { min_linear = std::min(min_linear, values[j]); }
+  }
+  if (min_linear <= epsilon_adjust) {
+    const f_t delta = -min_linear + epsilon_adjust;
+    for (i_t j = 0; j < linear_end; ++j) {
+      if (linear_mask[j]) { values[j] += delta; }
+    }
+  }
+
+  i_t off = 0;
+  for (i_t q_k : cone_dims) {
+    const i_t base = linear_end + off;
+    f_t tail_sq    = 0.0;
+    for (i_t j = 1; j < q_k; ++j) {
+      const f_t t = values[base + j];
+      tail_sq += t * t;
+    }
+    const f_t tail_norm = std::sqrt(tail_sq);
+    if (values[base] <= tail_norm + epsilon_adjust) { values[base] = tail_norm + epsilon_adjust; }
+    off += q_k;
+  }
+}
+
 template <typename f_t>
 [[maybe_unused]] static void pairwise_multiply(
   f_t* a, f_t* b, f_t* out, int size, rmm::cuda_stream_view stream)
@@ -239,7 +272,6 @@ class iteration_data_t {
       complementarity_residual_norm_save(inf),
       diag(lp.num_cols),
       inv_diag(lp.num_cols),
-      inv_sqrt_diag(lp.num_cols),
       AD(lp.num_cols, lp.num_rows, 0),
       AT(lp.num_rows, lp.num_cols, 0),
       ADAT(lp.num_rows, lp.num_rows, 0),
@@ -540,12 +572,9 @@ class iteration_data_t {
       }
 
       inv_diag.set_scalar(1.0);
-      if (use_augmented) { diag.multiply_scalar(-1.0); }
       if (n_upper_bounds > 0 || (has_Q && !use_augmented)) { diag.inverse(inv_diag); }
       // TMP diag and inv_diag should directly created and filled on the GPU
       raft::copy(d_inv_diag.data(), inv_diag.data(), inv_diag.size(), stream_view_);
-      inv_sqrt_diag.set_scalar(1.0);
-      if (n_upper_bounds > 0 || (has_Q && !use_augmented)) { inv_diag.sqrt(inv_sqrt_diag); }
     }
 
     if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { return; }
@@ -1912,7 +1941,6 @@ class iteration_data_t {
 
   dense_vector_t<i_t, f_t> diag;
   pinned_dense_vector_t<i_t, f_t> inv_diag;
-  dense_vector_t<i_t, f_t> inv_sqrt_diag;
 
   rmm::device_uvector<f_t> d_original_A_values;
 
@@ -2155,41 +2183,55 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
   const bool use_augmented          = data.use_augmented;
   const bool has_direct_free_linear = data.n_direct_free_linear > 0;
 
-  // SOCP: data-dependent initial point following SeDuMi (Sturm, 1999).
-  //   mu = sqrt((1 + ||b||_inf) * (1 + ||c||_inf))
-  //   primal and dual: x = mu * e_K,  z = mu * e_K
+  const barrier_dual_initial_point_t input_strategy =
+    static_cast<barrier_dual_initial_point_t>(settings.barrier_dual_initial_point);
+
+  const barrier_dual_initial_point_t init_strategy =
+    (data.has_cones() && input_strategy == barrier_dual_initial_point_t::Automatic)
+      ? barrier_dual_initial_point_t::SedumiMu
+      : input_strategy;
+
+  // SedumiMu: Sturm/SeDuMi-style mu-based primal+dual initial point.
+  //   mu = sqrt((1 + ||b||_inf) * (1 + ||c||_inf)); x = z = mu * e_K.
   // where e_K is the identity of the symmetric cone:
   //   LP block: e = 1,  SOC block: e = (sqrt(2), 0, ..., 0)
-  if (data.has_cones()) {
-    const i_t cs     = data.cone_start();
-    const f_t norm_b = vector_norm_inf<i_t, f_t>(lp.rhs);
-    const f_t norm_c = vector_norm_inf<i_t, f_t>(lp.objective);
-    const f_t mu     = std::sqrt((1.0 + norm_b) * (1.0 + norm_c));
-    const f_t sqrt2  = std::sqrt(2.0);
-    const f_t x_soc  = mu * sqrt2;
-    const f_t z_soc  = mu * sqrt2;
-    // Linear orthant
-    for (i_t j = 0; j < cs; ++j) {
+  // Full primal+dual point; no factorization/solve (main loop factorizes later).
+  if (init_strategy == barrier_dual_initial_point_t::SedumiMu) {
+    const f_t norm_b     = vector_norm_inf<i_t, f_t>(lp.rhs);
+    const f_t norm_c     = vector_norm_inf<i_t, f_t>(lp.objective);
+    const f_t mu         = std::sqrt((1.0 + norm_b) * (1.0 + norm_c));
+    const f_t sqrt2      = std::sqrt(2.0);
+    const i_t linear_end = data.linear_xz_size(lp.num_cols);
+
+    // Linear orthant: x = z = mu * e, with e = 1
+    for (i_t j = 0; j < linear_end; ++j) {
       data.x[j] = mu;
       data.z[j] = mu;
     }
     if (has_direct_free_linear) {
       for (i_t j : presolve_info.direct_free_variables) {
-        if (j < cs) { data.z[j] = 0.0; }
+        if (j < linear_end) { data.z[j] = 0.0; }
       }
     }
-    // SOC blocks
-    i_t off = 0;
-    for (size_t k = 0; k < lp.second_order_cone_dims.size(); k++) {
-      i_t q_k          = lp.second_order_cone_dims[k];
-      data.x[cs + off] = x_soc;
-      data.z[cs + off] = z_soc;
-      for (i_t j = 1; j < q_k; ++j) {
-        data.x[cs + off + j] = 0.0;
-        data.z[cs + off + j] = 0.0;
+
+    // SOC blocks: x = z = mu * e, with e = (sqrt(2), 0, ..., 0)
+    if (data.has_cones()) {
+      const i_t cs    = data.cone_start();
+      const f_t x_soc = mu * sqrt2;
+      const f_t z_soc = mu * sqrt2;
+      i_t off         = 0;
+      for (size_t k = 0; k < lp.second_order_cone_dims.size(); k++) {
+        i_t q_k          = lp.second_order_cone_dims[k];
+        data.x[cs + off] = x_soc;
+        data.z[cs + off] = z_soc;
+        for (i_t j = 1; j < q_k; ++j) {
+          data.x[cs + off + j] = 0.0;
+          data.z[cs + off + j] = 0.0;
+        }
+        off += q_k;
       }
-      off += q_k;
     }
+
     data.y.set_scalar(0.0);
     if (data.n_upper_bounds > 0) {
       data.w.set_scalar(mu);
@@ -2229,15 +2271,18 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
   dense_vector_t<i_t, f_t> DinvFu(lp.num_cols);  // DinvFu = Dinv * Fu
   data.inv_diag.pairwise_product(Fu, DinvFu);
   dense_vector_t<i_t, f_t> q(lp.num_rows);
+  const i_t aug_base = lp.num_cols + lp.num_rows;
+  const i_t aug_size = use_augmented ? std::max<i_t>(aug_base, data.device_augmented.n) : aug_base;
   if (use_augmented) {
-    dense_vector_t<i_t, f_t> rhs(lp.num_cols + lp.num_rows);
+    dense_vector_t<i_t, f_t> rhs(aug_size);
+    rhs.set_scalar(0.0);
     for (i_t k = 0; k < lp.num_cols; k++) {
       rhs[k] = -Fu[k];
     }
     for (i_t k = 0; k < lp.num_rows; k++) {
       rhs[lp.num_cols + k] = rhs_x[k];
     }
-    dense_vector_t<i_t, f_t> soln(lp.num_cols + lp.num_rows);
+    dense_vector_t<i_t, f_t> soln(aug_size);
     i_t solve_status = data.chol->solve(rhs, soln);
     struct op_t {
       op_t(iteration_data_t<i_t, f_t>& data) : data_(data) {}
@@ -2255,7 +2300,11 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
       }
     } op(data);
 
-    if (settings.barrier_iterative_refinement) {
+    // Initial-point IR is LP/QP-only. Cone problems and direct-free linear vars are excluded.
+    if (settings.barrier_iterative_refinement && !data.has_cones() &&
+        data.n_direct_free_linear == 0) {
+      raft::copy(
+        data.d_diag_.data(), data.diag.data(), data.diag.size(), data.handle_ptr->get_stream());
       iterative_refinement<i_t, f_t, op_t>(op, rhs, soln);
     }
 
@@ -2330,8 +2379,21 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
   }
 
   float64_t epsilon_adjust = 10.0;
+  // Push entries into interior of nonnegative orthant and SOC.
+  const bool has_soc   = data.has_cones();
+  const i_t linear_end = has_soc ? data.cone_start() : lp.num_cols;
+  auto ensure_interior = [&](dense_vector_t<i_t, f_t>& values,
+                             const std::vector<i_t>& linear_mask) {
+    if (has_soc) {
+      ensure_initial_point_interior(
+        values, epsilon_adjust, linear_mask, linear_end, lp.second_order_cone_dims);
+    } else {
+      values.ensure_positive(epsilon_adjust, linear_mask);
+    }
+  };
 
-  if (settings.barrier_dual_initial_point == -1 || settings.barrier_dual_initial_point == 0) {
+  if (init_strategy == barrier_dual_initial_point_t::Automatic ||
+      init_strategy == barrier_dual_initial_point_t::LustigMarstenShanno) {
     // Use the dual starting point suggested by the paper
     // On Implementing Mehrotra’s Predictor–Corrector Interior-Point Method for Linear Programming
     // Irvin J. Lustig, Roy E. Marsten, and David F. Shanno
@@ -2380,12 +2442,12 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
       }
     }
   } else if (use_augmented) {
-    dense_vector_t<i_t, f_t> dual_rhs(lp.num_cols + lp.num_rows);
+    dense_vector_t<i_t, f_t> dual_rhs(aug_size);
     dual_rhs.set_scalar(0.0);
     for (i_t k = 0; k < lp.num_cols; k++) {
       dual_rhs[k] = data.c[k];
     }
-    dense_vector_t<i_t, f_t> py(lp.num_cols + lp.num_rows);
+    dense_vector_t<i_t, f_t> py(aug_size);
     data.chol->solve(dual_rhs, py);
     for (i_t k = 0; k < lp.num_cols; k++) {
       data.z[k] = py[k];
@@ -2399,7 +2461,6 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
     data.v.multiply_scalar(-1.0);
 
     data.v.ensure_positive(epsilon_adjust);
-    data.z.ensure_positive(epsilon_adjust, nonnegative_z);
   } else {
     // First compute rhs = A*Dinv*c
     dense_vector_t<i_t, f_t> rhs(lp.num_rows);
@@ -2423,7 +2484,6 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
     data.gather_upper_bounds(data.z, data.v);
     data.v.multiply_scalar(-1.0);
     data.v.ensure_positive(epsilon_adjust);
-    data.z.ensure_positive(epsilon_adjust, nonnegative_z);
   }
 
   // Verify A'*y + z - E*v  - Q*x = c
@@ -2441,6 +2501,7 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
   settings.log.printf("||A^T y + z - E*v - Q*x - c ||: %e\n",
                       vector_norm2<i_t, f_t>(init_dual_residual));
 #endif
+
   // Make sure (w, x, v, z) > 0. Skip free variables being handled directly.
   data.w.ensure_positive(epsilon_adjust);
   std::vector<i_t> nonnegative_variables(data.x.size(), 1);
@@ -2449,7 +2510,8 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
       nonnegative_variables[j] = 0;
     }
   }
-  data.x.ensure_positive(epsilon_adjust, nonnegative_variables);
+  ensure_interior(data.z, nonnegative_z);
+  ensure_interior(data.x, nonnegative_variables);
   // Direct free variables: reduced cost z = 0 (no complementarity condition).
   if (has_direct_free_linear) {
     for (i_t j : presolve_info.direct_free_variables) {

cpp/src/barrier/barrier.hpp

@@ -8,6 +8,7 @@
 
 #include <linear_algebra/dense_vector.hpp>
 
+#include <cuopt/mathematical_optimization/constants.h>
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
 #include <dual_simplex/solution.hpp>
@@ -18,6 +19,21 @@
 #include <rmm/device_uvector.hpp>
 namespace cuopt::mathematical_optimization::barrier {
 
+/**
+ * @brief Barrier primal-dual initial-point strategy.
+ *
+ * Automatic: use Lustig-Marsten-Shanno for LP/QP; Sturm/SeDuMi mu-based point for conic problems.
+ * LustigMarstenShanno: Mehrotra-style dual start (Lustig, Marsten, Shanno, SIAM J. Optim. 1992).
+ * DualLeastSquares: solve augmented or ADAT dual least-squares system.
+ * SedumiMu: Sturm/SeDuMi mu-based primal+dual point (no factorization).
+ */
+enum barrier_dual_initial_point_t : int {
+  Automatic           = CUOPT_BARRIER_DUAL_INITIAL_POINT_AUTOMATIC,
+  LustigMarstenShanno = CUOPT_BARRIER_DUAL_INITIAL_POINT_LUSTIG_MARSTEN_SHANNO,
+  DualLeastSquares    = CUOPT_BARRIER_DUAL_INITIAL_POINT_DUAL_LEAST_SQUARES,
+  SedumiMu            = CUOPT_BARRIER_DUAL_INITIAL_POINT_SEDUMI_MU
+};
+
 /** Validates SOC layout on an simplex::lp_problem_t before barrier presolve/solve. */
 template <typename i_t, typename f_t>
 bool validate_barrier_cone_layout(const simplex::lp_problem_t<i_t, f_t>& problem,

cpp/src/dual_simplex/simplex_solver_settings.hpp

@@ -167,8 +167,8 @@ struct simplex_solver_settings_t {
   i_t augmented;  // -1 automatic, 0 to solve with ADAT, 1 to solve with augmented system
   i_t dualize;    // -1 automatic, 0 to not dualize, 1 to dualize
   i_t ordering;   // -1 automatic, 0 to use nested dissection, 1 to use AMD
-  i_t barrier_dual_initial_point;  // -1 automatic, 0 to use Lustig, Marsten, and Shanno initial
-                                   // point, 1 to use initial point form dual least squares problem
+  i_t barrier_dual_initial_point;  // -1 automatic, 0 Lustig-Marsten-Shanno,
+                                   // 1 dual least squares, 2 SeDuMi mu-based
   bool check_Q;                    // true to check if Q is positive semidefinite
   bool crossover;                  // true to do crossover, false to not
   i_t refactor_frequency;          // number of basis updates before refactorization

cpp/src/math_optimization/solver_settings.cu

@@ -137,7 +137,7 @@ solver_settings_t<i_t, f_t>::solver_settings_t() : pdlp_settings(), mip_settings
     {CUOPT_FOLDING, &pdlp_settings.folding, -1, 1, -1},
     {CUOPT_DUALIZE, &pdlp_settings.dualize, -1, 1, -1},
     {CUOPT_ORDERING, &pdlp_settings.ordering, -1, 1, -1},
-    {CUOPT_BARRIER_DUAL_INITIAL_POINT, &pdlp_settings.barrier_dual_initial_point, -1, 1, -1},
+    {CUOPT_BARRIER_DUAL_INITIAL_POINT, &pdlp_settings.barrier_dual_initial_point, -1, 2, -1},
     {CUOPT_MIP_CUT_PASSES, &mip_settings.max_cut_passes, -1, std::numeric_limits<i_t>::max(), 10},
     {CUOPT_MIP_MIXED_INTEGER_ROUNDING_CUTS, &mip_settings.mir_cuts, -1, 1, -1},
     {CUOPT_MIP_MIXED_INTEGER_GOMORY_CUTS, &mip_settings.mixed_integer_gomory_cuts, -1, 1, -1},

python/cuopt_server/cuopt_server/tests/test_lp.py

@@ -180,7 +180,7 @@ def test_barrier_solver_options(
     - cudss_deterministic: True for deterministic, False for
       nondeterministic
     - barrier_dual_initial_point: (-1) automatic, (0) Lustig-Marsten-Shanno,
-      (1) dual least squares
+      (1) dual least squares, (2) Sturm/SeDuMi mu-based primal+dual
     """
     data = get_std_data_for_lp()
 

python/cuopt_server/cuopt_server/utils/linear_programming/data_definition.py

@@ -501,7 +501,8 @@ class SolverConfig(BaseModel):
         description="Set the type of dual initial point to use for the barrier"
         "solver. -1 for automatic, 0 to use Lustig, Marsten, and Shanno"
         "initial point, 1 to use initial point from a dual least squares"
-        "problem",
+        "problem, 2 to use Sturm/SeDuMi mu-based primal+dual"
+        "point",
     )
     eliminate_dense_columns: Optional[bool] = Field(
         default=True,