Add setting to change the dual initial point used by barrier

cpp/include/cuopt/linear_programming/constants.h

@@ -54,6 +54,7 @@
 #define CUOPT_AUGMENTED                   "augmented"
 #define CUOPT_DUALIZE                     "dualize"
 #define CUOPT_ORDERING                    "ordering"
+#define CUOPT_BARRIER_DUAL_INITIAL_POINT  "barrier_dual_initial_point"
 #define CUOPT_ELIMINATE_DENSE_COLUMNS     "eliminate_dense_columns"
 #define CUOPT_CUDSS_DETERMINISTIC         "cudss_deterministic"
 #define CUOPT_PRESOLVE                    "presolve"

cpp/include/cuopt/linear_programming/pdlp/solver_settings.hpp

@@ -214,6 +214,7 @@ class pdlp_solver_settings_t {
   i_t augmented{-1};
   i_t dualize{-1};
   i_t ordering{-1};
+  i_t barrier_dual_initial_point{-1};
   bool eliminate_dense_columns{true};
   bool save_best_primal_so_far{false};
   bool first_primal_feasible{false};

cpp/src/dual_simplex/barrier.cu

@@ -1619,7 +1619,8 @@ 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> dual_res(lp.num_cols);
-  if (1) {
+  float64_t epsilon_adjust = 10.0;
+  if (settings.barrier_dual_initial_point == -1 || settings.barrier_dual_initial_point == 0) {
     // 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
@@ -1686,6 +1687,8 @@ 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);
   } else {
     // First compute rhs = A*Dinv*c
     dense_vector_t<i_t, f_t> rhs(lp.num_rows);
@@ -1716,6 +1719,9 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
     // v = -E'*z
     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);
   }
 
   // Verify A'*y + z - E*v = c
@@ -1738,11 +1744,8 @@ 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 - c ||: %e\n", vector_norm2<i_t, f_t>(data.dual_residual));
 #endif
   // Make sure (w, x, v, z) > 0
-  float64_t epsilon_adjust = 10.0;
   data.w.ensure_positive(epsilon_adjust);
   data.x.ensure_positive(epsilon_adjust);
-  // data.v.ensure_positive(epsilon_adjust);
-  // data.z.ensure_positive(epsilon_adjust);
 #ifdef PRINT_INFO
   settings.log.printf("min v %e min z %e\n", data.v.minimum(), data.z.minimum());
 #endif

cpp/src/dual_simplex/simplex_solver_settings.hpp

@@ -70,6 +70,7 @@ struct simplex_solver_settings_t {
       augmented(0),
       dualize(-1),
       ordering(-1),
+      barrier_dual_initial_point(-1),
       crossover(false),
       refactor_frequency(100),
       iteration_log_frequency(1000),
@@ -128,17 +129,19 @@ struct simplex_solver_settings_t {
   bool barrier;               // true to use barrier method, false to use dual simplex method
   bool eliminate_dense_columns;  // true to eliminate dense columns from A*D*A^T
   i_t folding;                   // -1 automatic, 0 don't fold, 1 fold
-  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
-  bool crossover;          // true to do crossover, false to not
-  i_t refactor_frequency;  // number of basis updates before refactorization
-  i_t iteration_log_frequency;  // number of iterations between log updates
-  i_t first_iteration_log;      // number of iterations to log at beginning of solve
-  i_t num_threads;              // number of threads to use
-  i_t random_seed;              // random seed
-  i_t num_bfs_threads;          // number of threads dedicated to the best-first search
-  i_t num_diving_threads;       // number of threads dedicated to diving
+  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
+  bool crossover;                  // true to do crossover, false to not
+  i_t refactor_frequency;          // number of basis updates before refactorization
+  i_t iteration_log_frequency;     // number of iterations between log updates
+  i_t first_iteration_log;         // number of iterations to log at beginning of solve
+  i_t num_threads;                 // number of threads to use
+  i_t random_seed;                 // random seed
+  i_t num_bfs_threads;             // number of threads dedicated to the best-first search
+  i_t num_diving_threads;          // number of threads dedicated to diving
   i_t inside_mip;  // 0 if outside MIP, 1 if inside MIP at root node, 2 if inside MIP at leaf node
   std::function<void(std::vector<f_t>&, f_t)> solution_callback;
   std::function<void()> heuristic_preemption_callback;

cpp/src/linear_programming/solve.cu

@@ -407,6 +407,7 @@ run_barrier(dual_simplex::user_problem_t<i_t, f_t>& user_problem,
   barrier_settings.augmented                       = settings.augmented;
   barrier_settings.dualize                         = settings.dualize;
   barrier_settings.ordering                        = settings.ordering;
+  barrier_settings.barrier_dual_initial_point      = settings.barrier_dual_initial_point;
   barrier_settings.barrier                         = true;
   barrier_settings.crossover                       = settings.crossover;
   barrier_settings.eliminate_dense_columns         = settings.eliminate_dense_columns;

cpp/src/linear_programming/solver_settings.cu

@@ -49,6 +49,7 @@ pdlp_solver_settings_t<i_t, f_t>::pdlp_solver_settings_t(const pdlp_solver_setti
     augmented(other.augmented),
     dualize(other.dualize),
     ordering(other.ordering),
+    barrier_dual_initial_point(other.barrier_dual_initial_point),
     cudss_deterministic(other.cudss_deterministic),
     eliminate_dense_columns(other.eliminate_dense_columns),
     save_best_primal_so_far(other.save_best_primal_so_far),

cpp/src/math_optimization/solver_settings.cu

@@ -95,7 +95,8 @@ solver_settings_t<i_t, f_t>::solver_settings_t() : pdlp_settings(), mip_settings
     {CUOPT_AUGMENTED, &pdlp_settings.augmented, -1, 1, -1},
     {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_ORDERING, &pdlp_settings.ordering, -1, 1, -1},
+    {CUOPT_BARRIER_DUAL_INITIAL_POINT, &pdlp_settings.barrier_dual_initial_point, -1, 1, -1}
   };
 
     // Bool parameters

python/cuopt/cuopt/linear_programming/solver/solver_parameters.pyx

@@ -77,6 +77,7 @@ cdef extern from "cuopt/linear_programming/constants.h": # noqa
     cdef const char* c_CUOPT_ELIMINATE_DENSE_COLUMNS "CUOPT_ELIMINATE_DENSE_COLUMNS" # noqa
     cdef const char* c_CUOPT_CUDSS_DETERMINISTIC "CUOPT_CUDSS_DETERMINISTIC" # noqa
     cdef const char* c_CUOPT_ORDERING "CUOPT_ORDERING" # noqa
+    cdef const char* c_CUOPT_BARRIER_DUAL_INITIAL_POINT "CUOPT_BARRIER_DUAL_INITIAL_POINT" # noqa
 
 # Create Python string constants from C string literals
 CUOPT_ABSOLUTE_DUAL_TOLERANCE = c_CUOPT_ABSOLUTE_DUAL_TOLERANCE.decode('utf-8') # noqa
@@ -117,3 +118,4 @@ CUOPT_DUALIZE = c_CUOPT_DUALIZE.decode('utf-8') # noqa
 CUOPT_ELIMINATE_DENSE_COLUMNS = c_CUOPT_ELIMINATE_DENSE_COLUMNS.decode('utf-8') # noqa
 CUOPT_CUDSS_DETERMINISTIC = c_CUOPT_CUDSS_DETERMINISTIC.decode('utf-8') # noqa
 CUOPT_ORDERING = c_CUOPT_ORDERING.decode('utf-8') # noqa
+CUOPT_BARRIER_DUAL_INITIAL_POINT = c_CUOPT_BARRIER_DUAL_INITIAL_POINT.decode('utf-8') # noqa

python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.py

@@ -20,6 +20,7 @@
     CUOPT_ABSOLUTE_GAP_TOLERANCE,
     CUOPT_ABSOLUTE_PRIMAL_TOLERANCE,
     CUOPT_AUGMENTED,
+    CUOPT_BARRIER_DUAL_INITIAL_POINT,
     CUOPT_CROSSOVER,
     CUOPT_CUDSS_DETERMINISTIC,
     CUOPT_DUAL_INFEASIBLE_TOLERANCE,
@@ -390,6 +391,9 @@ def toDict(self):
             "folding": self.get_parameter(CUOPT_FOLDING),
             "dualize": self.get_parameter(CUOPT_DUALIZE),
             "ordering": self.get_parameter(CUOPT_ORDERING),
+            "barrier_dual_initial_point": self.get_parameter(
+                CUOPT_BARRIER_DUAL_INITIAL_POINT
+            ),
             "eliminate_dense_columns": self.get_parameter(
                 CUOPT_ELIMINATE_DENSE_COLUMNS
             ),

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

@@ -555,6 +555,13 @@ class SolverConfig(StrictModel):
         description="Set the type of ordering to use for the barrier solver."
         "-1 for automatic, 0 to use cuDSS default ordering, 1 to use AMD",
     )
+    barrier_dual_initial_point: Optional[int] = Field(
+        default=-1,
+        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",
+    )
     eliminate_dense_columns: Optional[bool] = Field(
         default=True,
         description="Set if dense columns should be eliminated from the "

python/cuopt_server/cuopt_server/utils/linear_programming/solver.py

@@ -26,6 +26,7 @@
     CUOPT_ABSOLUTE_GAP_TOLERANCE,
     CUOPT_ABSOLUTE_PRIMAL_TOLERANCE,
     CUOPT_AUGMENTED,
+    CUOPT_BARRIER_DUAL_INITIAL_POINT,
     CUOPT_CROSSOVER,
     CUOPT_CUDSS_DETERMINISTIC,
     CUOPT_DUAL_INFEASIBLE_TOLERANCE,
@@ -451,6 +452,11 @@ def is_mip(var_types):
             solver_settings.set_parameter(
                 CUOPT_ORDERING, solver_config.ordering
             )
+        if solver_config.barrier_dual_initial_point is not None:
+            solver_settings.set_parameter(
+                CUOPT_BARRIER_DUAL_INITIAL_POINT,
+                solver_config.barrier_dual_initial_point,
+            )
         if solver_config.eliminate_dense_columns is not None:
             solver_settings.set_parameter(
                 CUOPT_ELIMINATE_DENSE_COLUMNS,