Add more strategies of init point computation in barrier

cpp/include/cuopt/mathematical_optimization/constants.h

@@ -156,6 +156,9 @@
 /* @brief QCQP (barrier) scaling hyper-parameters */
 #define CUOPT_QCQP_HYPER_RUIZ_EQUILIBRATION "qcqp_hyper_ruiz_equilibration"
 
+/* @brief Barrier initial point safeguard */
+#define CUOPT_BARRIER_INITIAL_POINT_SAFEGUARD "barrier_initial_point_safeguard"
+
 /* @brief MIP determinism mode constants */
 #define CUOPT_MODE_OPPORTUNISTIC 0
 #define CUOPT_MODE_DETERMINISTIC 1
@@ -209,6 +212,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_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/include/cuopt/mathematical_optimization/pdlp/solver_settings.hpp

@@ -294,13 +294,16 @@ class pdlp_solver_settings_t {
   i_t augmented{-1};
   i_t dualize{-1};
   i_t ordering{-1};
-  i_t barrier_dual_initial_point{-1};
+  barrier_dual_initial_point_t barrier_dual_initial_point{barrier_dual_initial_point_t::Automatic};
   i_t postsolve_info{-1};
   i_t barrier_presolve_bound_free_variables{-1};  // -1 automatic, 0 disabled, 1 enabled
   // Ruiz equilibration for QCQP (barrier) scaling: -1 automatic (row/column
   // imbalance heuristic), 0 disabled, 1 enabled. Distinct from PDLP's own Ruiz
   // scaling in pdlp_hyper_params_t.
   i_t qcqp_ruiz_equilibration{-1};
+  // Margin used to push the barrier method's initial iterate into the interior of the
+  // nonnegative orthant / SOC (values are shifted to be at least this far from the boundary).
+  f_t barrier_initial_point_safeguard{10.0};
   bool eliminate_dense_columns{true};
   pdlp_precision_t pdlp_precision{pdlp_precision_t::DefaultPrecision};
   bool barrier_iterative_refinement{true};

cpp/include/cuopt/mathematical_optimization/utilities/internals.hpp

@@ -142,5 +142,20 @@ enum presolver_t : int {
   PSLP    = CUOPT_PRESOLVE_PSLP
 };
 
+/**
+ * @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_LEAST_SQUARES,
+  SedumiMu            = CUOPT_BARRIER_DUAL_INITIAL_POINT_SEDUMI_MU
+};
+
 }  // namespace mathematical_optimization
 }  // namespace cuopt

cpp/src/barrier/barrier.cu

@@ -97,6 +97,34 @@ 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)
+{
+  // Linear shift
+  std::vector<i_t> linear_only_mask(values.size(), 0);
+  std::copy(linear_mask.begin(), linear_mask.begin() + linear_end, linear_only_mask.begin());
+  values.ensure_positive(epsilon_adjust, linear_only_mask);
+
+  // Cone shift
+  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;
+  }
+}
+
 // -1 automatic: enable for cones, disable otherwise; 0 off; 1 on
 template <typename i_t, typename f_t>
 bool should_use_adaptive_regularization(const simplex_solver_settings_t<i_t, f_t>& settings,
@@ -2182,41 +2210,54 @@ 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 = 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);
@@ -2359,9 +2400,18 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
 #endif
   }
 
-  float64_t epsilon_adjust = 10.0;
+  const f_t epsilon_adjust = settings.barrier_initial_point_safeguard;
+  // 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) {
+    ensure_initial_point_interior(
+      values, epsilon_adjust, linear_mask, linear_end, lp.second_order_cone_dims);
+  };
 
-  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
@@ -2430,7 +2480,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);
@@ -2454,7 +2503,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
@@ -2472,6 +2520,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);
@@ -2480,7 +2529,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,8 @@
 
 #include <linear_algebra/dense_vector.hpp>
 
+#include <cuopt/mathematical_optimization/constants.h>
+#include <cuopt/mathematical_optimization/utilities/internals.hpp>
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
 #include <dual_simplex/solution.hpp>

cpp/src/dual_simplex/simplex_solver_settings.hpp

@@ -9,6 +9,7 @@
 
 #include <cuopt/mathematical_optimization/mip/diving_hyper_params.hpp>
 #include <cuopt/mathematical_optimization/mip/submip_hyper_params.hpp>
+#include <cuopt/mathematical_optimization/utilities/internals.hpp>
 
 #include <dual_simplex/logger.hpp>
 #include <math_optimization/types.hpp>
@@ -78,10 +79,11 @@ struct simplex_solver_settings_t {
       augmented(0),
       dualize(-1),
       ordering(-1),
-      barrier_dual_initial_point(-1),
+      barrier_dual_initial_point(barrier_dual_initial_point_t::Automatic),
       postsolve_info(-1),
       barrier_presolve_bound_free_variables(-1),
       qcqp_ruiz_equilibration(-1),
+      barrier_initial_point_safeguard(10.0),
       check_Q(false),
       crossover(false),
       refactor_frequency(100),
@@ -174,11 +176,14 @@ 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 postsolve_info;              // -1 automatic (disabled), 0 disabled, 1 enabled
+  barrier_dual_initial_point_t
+    barrier_dual_initial_point;               // -1 automatic, 0 Lustig-Marsten-Shanno,
+                                              // 1 dual least squares, 2 SeDuMi mu-based
+  i_t postsolve_info;                         // -1 automatic (disabled), 0 disabled, 1 enabled
   i_t barrier_presolve_bound_free_variables;  // -1 automatic, 0 disabled, 1 enabled
-  i_t qcqp_ruiz_equilibration;     // -1 automatic (imbalance heuristic), 0 disabled, 1 enabled
+  i_t qcqp_ruiz_equilibration;          // -1 automatic (imbalance heuristic), 0 disabled, 1 enabled
+  f_t barrier_initial_point_safeguard;  // margin pushing the barrier initial iterate into
+  // the interior of the nonnegative orthant / SOC
   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/grpc/codegen/field_registry.yaml

@@ -558,6 +558,7 @@ pdlp_settings:
   - barrier_dual_initial_point:
       field_num: 26
       type: int32
+      from_proto_cast: "barrier_dual_initial_point_t"
       optional: true
   - eliminate_dense_columns:
       field_num: 27

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

@@ -31,7 +31,7 @@
   pb_settings->set_augmented(settings.augmented);
   pb_settings->set_dualize(settings.dualize);
   pb_settings->set_ordering(settings.ordering);
-  pb_settings->set_barrier_dual_initial_point(settings.barrier_dual_initial_point);
+  pb_settings->set_barrier_dual_initial_point(static_cast<int32_t>(settings.barrier_dual_initial_point));
   pb_settings->set_eliminate_dense_columns(settings.eliminate_dense_columns);
   pb_settings->set_barrier_iterative_refinement(settings.barrier_iterative_refinement);
   pb_settings->set_barrier_adaptive_regularization(settings.barrier_adaptive_regularization);

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

@@ -68,7 +68,7 @@
     settings.ordering = pb_settings.ordering();
   }
   if (pb_settings.has_barrier_dual_initial_point()) {
-    settings.barrier_dual_initial_point = pb_settings.barrier_dual_initial_point();
+    settings.barrier_dual_initial_point = static_cast<barrier_dual_initial_point_t>(pb_settings.barrier_dual_initial_point());
   }
   if (pb_settings.has_eliminate_dense_columns()) {
     settings.eliminate_dense_columns = pb_settings.eliminate_dense_columns();

cpp/src/math_optimization/solver_settings.cu

@@ -104,6 +104,7 @@ solver_settings_t<i_t, f_t>::solver_settings_t() : pdlp_settings(), mip_settings
     {CUOPT_MIP_CUT_CHANGE_THRESHOLD, &mip_settings.cut_change_threshold, f_t(-1.0), std::numeric_limits<f_t>::infinity(), f_t(-1.0)},
     {CUOPT_MIP_CUT_MIN_ORTHOGONALITY, &mip_settings.cut_min_orthogonality, f_t(0.0), f_t(1.0), f_t(0.5)},
     {CUOPT_BARRIER_STEP_SCALE, &pdlp_settings.barrier_step_scale, f_t(0.5), f_t(0.9999), f_t(0.9)},
+    {CUOPT_BARRIER_INITIAL_POINT_SAFEGUARD, &pdlp_settings.barrier_initial_point_safeguard, f_t(0.0), std::numeric_limits<f_t>::infinity(), f_t(10.0), "margin pushing the barrier initial iterate into the interior of the nonnegative orthant / SOC"},
     // MIP heuristic hyper-parameters (hidden from default --help: name contains "hyper_")
     {CUOPT_MIP_HYPER_HEURISTIC_ROOT_LP_TIME_RATIO, &mip_settings.heuristic_params.root_lp_time_ratio, f_t(0.0), f_t(1.0), f_t(0.1), "fraction of total time for root LP"},
     {CUOPT_MIP_HYPER_HEURISTIC_ROOT_LP_MAX_TIME, &mip_settings.heuristic_params.root_lp_max_time, f_t(0.0), std::numeric_limits<f_t>::infinity(), f_t(15.0), "hard cap on root LP seconds"},
@@ -137,7 +138,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, reinterpret_cast<int*>(&pdlp_settings.barrier_dual_initial_point), -1, 2, -1},
     {CUOPT_POSTSOLVE_INFO, &pdlp_settings.postsolve_info, -1, 1, -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},

cpp/src/pdlp/solve.cu

@@ -510,6 +510,7 @@ std::tuple<simplex::lp_solution_t<i_t, f_t>, simplex::lp_status_t, f_t, f_t, f_t
   barrier_settings.postsolve_info             = settings.postsolve_info;
   barrier_settings.barrier_presolve_bound_free_variables =
     settings.barrier_presolve_bound_free_variables;
+  barrier_settings.barrier_initial_point_safeguard = settings.barrier_initial_point_safeguard;
   barrier_settings.barrier                         = true;
   barrier_settings.barrier_presolve                = true;
   barrier_settings.crossover                       = settings.crossover;

cpp/tests/linear_programming/grpc/grpc_client_test.cpp

@@ -2224,24 +2224,25 @@ TEST(MapperRoundtrip, PDLPSettingsAllFields)
   orig.tolerances.absolute_primal_tolerance   = 5e-7;
   orig.tolerances.relative_primal_tolerance   = 6e-7;
 
-  orig.time_limit                   = 99.5;
-  orig.iteration_limit              = 10000;
-  orig.log_to_console               = false;
-  orig.detect_infeasibility         = true;
-  orig.strict_infeasibility         = true;
-  orig.pdlp_solver_mode             = pdlp_solver_mode_t::Fast1;
-  orig.method                       = method_t::Barrier;
-  orig.presolver                    = presolver_t::Default;
-  orig.dual_postsolve               = true;
-  orig.crossover                    = true;
-  orig.num_gpus                     = 4;
-  orig.per_constraint_residual      = true;
-  orig.cudss_deterministic          = true;
-  orig.folding                      = 1;
-  orig.augmented                    = 1;
-  orig.dualize                      = 1;
-  orig.ordering                     = 2;
-  orig.barrier_dual_initial_point   = 1;
+  orig.time_limit              = 99.5;
+  orig.iteration_limit         = 10000;
+  orig.log_to_console          = false;
+  orig.detect_infeasibility    = true;
+  orig.strict_infeasibility    = true;
+  orig.pdlp_solver_mode        = pdlp_solver_mode_t::Fast1;
+  orig.method                  = method_t::Barrier;
+  orig.presolver               = presolver_t::Default;
+  orig.dual_postsolve          = true;
+  orig.crossover               = true;
+  orig.num_gpus                = 4;
+  orig.per_constraint_residual = true;
+  orig.cudss_deterministic     = true;
+  orig.folding                 = 1;
+  orig.augmented               = 1;
+  orig.dualize                 = 1;
+  orig.ordering                = 2;
+  orig.barrier_dual_initial_point =
+    cuopt::mathematical_optimization::barrier_dual_initial_point_t::LustigMarstenShanno;
   orig.eliminate_dense_columns      = true;
   orig.barrier_iterative_refinement = false;  // not the default true, to detect overwrite-on-decode
   orig.barrier_step_scale           = 0.75;   // not the default 0.9
@@ -2282,7 +2283,8 @@ TEST(MapperRoundtrip, PDLPSettingsAllFields)
   EXPECT_EQ(restored.augmented, 1);
   EXPECT_EQ(restored.dualize, 1);
   EXPECT_EQ(restored.ordering, 2);
-  EXPECT_EQ(restored.barrier_dual_initial_point, 1);
+  EXPECT_EQ(restored.barrier_dual_initial_point,
+            cuopt::mathematical_optimization::barrier_dual_initial_point_t::LustigMarstenShanno);
   EXPECT_EQ(restored.eliminate_dense_columns, true);
   EXPECT_EQ(restored.barrier_iterative_refinement, false);
   EXPECT_DOUBLE_EQ(restored.barrier_step_scale, 0.75);

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()