Skip to content

Papilo primal/dual crush #1104

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
aliceb-nv wants to merge 7 commits into
base: main
Choose a base branch
from
Open

Papilo primal/dual crush #1104

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
871f60a
primal crush for papilo
aliceb-nv Apr 10, 2026
b27c1f0
tentative crash fix
aliceb-nv Apr 9, 2026
bb9c91b
crush dual vectors and reduced costs
aliceb-nv Apr 14, 2026
e709130
more tests, clenaup
aliceb-nv Apr 14, 2026
2b69c95
Merge branch 'main' into papilo-crush
aliceb-nv Apr 14, 2026
d69e8d0
cleanup the KKT check
aliceb-nv Apr 15, 2026
d47d709
cleanup
aliceb-nv Apr 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 26 additions & 5 deletions cpp/src/mip_heuristics/diversity/diversity_manager.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
#include "diversity_manager.cuh"

#include <mip_heuristics/mip_constants.hpp>
#include <mip_heuristics/presolve/third_party_presolve.hpp>

#include <mip_heuristics/presolve/conflict_graph/clique_table.cuh>
#include <mip_heuristics/presolve/probing_cache.cuh>
Expand Down Expand Up @@ -182,9 +183,26 @@ void diversity_manager_t<i_t, f_t>::add_user_given_solutions(
std::vector<solution_t<i_t, f_t>>& initial_sol_vector)
{
raft::common::nvtx::range fun_scope("add_user_given_solutions");
for (const auto& init_sol : context.settings.initial_solutions) {
const bool has_papilo = problem_ptr->has_papilo_presolve_data();
const i_t papilo_orig_n = problem_ptr->get_papilo_original_num_variables();
for (size_t sol_idx = 0; sol_idx < context.settings.initial_solutions.size(); ++sol_idx) {
const auto& init_sol = context.settings.initial_solutions[sol_idx];
solution_t<i_t, f_t> sol(*problem_ptr);
rmm::device_uvector<f_t> init_sol_assignment(*init_sol, sol.handle_ptr->get_stream());

if (has_papilo) {
cuopt_assert((i_t)init_sol_assignment.size() == papilo_orig_n,
"Initial solution size mismatch");
std::vector<f_t> h_original = host_copy(init_sol_assignment, sol.handle_ptr->get_stream());
std::vector<f_t> h_crushed;
problem_ptr->presolve_data.papilo_presolve_ptr->crush_primal_solution(h_original, h_crushed);
expand_device_copy(init_sol_assignment, h_crushed, sol.handle_ptr->get_stream());
CUOPT_LOG_DEBUG("Crushed initial solution %d through Papilo (%d -> %d vars)",
sol_idx,
papilo_orig_n,
h_crushed.size());
}

if (problem_ptr->pre_process_assignment(init_sol_assignment)) {
relaxed_lp_settings_t lp_settings;
lp_settings.time_limit = std::min(60., timer.remaining_time() / 2);
Expand All @@ -202,10 +220,10 @@ void diversity_manager_t<i_t, f_t>::add_user_given_solutions(
sol.handle_ptr->get_stream());
bool is_feasible = sol.compute_feasibility();
cuopt_func_call(sol.test_variable_bounds(true));
CUOPT_LOG_INFO("Adding initial solution success! feas %d objective %f excess %f",
is_feasible,
sol.get_user_objective(),
sol.get_total_excess());
CUOPT_LOG_DEBUG("Adding initial solution success! feas %d objective %f excess %f",
is_feasible,
sol.get_user_objective(),
sol.get_total_excess());
population.run_solution_callbacks(sol);
initial_sol_vector.emplace_back(std::move(sol));
} else {
Expand Down Expand Up @@ -418,6 +436,9 @@ solution_t<i_t, f_t> diversity_manager_t<i_t, f_t>::run_solver()
CUOPT_LOG_INFO("GPU heuristics disabled via CUOPT_DISABLE_GPU_HEURISTICS=1");
population.initialize_population();
population.allocate_solutions();
if (check_b_b_preemption()) { return population.best_feasible(); }
add_user_given_solutions(initial_sol_vector);
population.add_solutions_from_vec(std::move(initial_sol_vector));

while (!check_b_b_preemption()) {
std::this_thread::sleep_for(std::chrono::milliseconds(100));
Expand Down
249 changes: 245 additions & 4 deletions cpp/src/mip_heuristics/presolve/third_party_presolve.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -137,8 +137,10 @@ papilo::Problem<f_t> build_papilo_problem(const optimization_problem_t<i_t, f_t>
}
}

for (size_t i = 0; i < h_var_types.size(); ++i) {
builder.setColIntegral(i, h_var_types[i] == var_t::INTEGER);
if (category == problem_category_t::MIP) {
for (size_t i = 0; i < h_var_types.size(); ++i) {
builder.setColIntegral(i, h_var_types[i] == var_t::INTEGER);
}
}

if (!h_constr_lb.empty() && !h_constr_ub.empty()) {
Expand Down Expand Up @@ -562,8 +564,6 @@ void set_presolve_methods(papilo::Presolve<f_t>& presolver,
presolver.addPresolveMethod(uptr(new papilo::SimpleProbing<f_t>()));
presolver.addPresolveMethod(uptr(new papilo::ParallelRowDetection<f_t>()));
presolver.addPresolveMethod(uptr(new papilo::ParallelColDetection<f_t>()));

presolver.addPresolveMethod(uptr(new papilo::SingletonStuffing<f_t>()));
presolver.addPresolveMethod(uptr(new papilo::DualFix<f_t>()));
presolver.addPresolveMethod(uptr(new papilo::SimplifyInequalities<f_t>()));
presolver.addPresolveMethod(uptr(new papilo::CliqueMerging<f_t>()));
Expand All @@ -574,6 +574,11 @@ void set_presolve_methods(papilo::Presolve<f_t>& presolver,
presolver.addPresolveMethod(uptr(new papilo::Probing<f_t>()));

if (!dual_postsolve) {
// SingletonStuffing causes dual crushing failures on:
// tr12-30, ns1208400, gmu-35-50, dws008-01, neos-1445765,
// neos-5107597-kakapo, rocI-4-11, traininstance2, traininstance6,
// radiationm18-12-05, rococoB10-011000, b1c1s1
presolver.addPresolveMethod(uptr(new papilo::SingletonStuffing<f_t>()));
presolver.addPresolveMethod(uptr(new papilo::DualInfer<f_t>()));
presolver.addPresolveMethod(uptr(new papilo::SimpleSubstitution<f_t>()));
presolver.addPresolveMethod(uptr(new papilo::Sparsify<f_t>()));
Expand Down Expand Up @@ -876,6 +881,242 @@ void third_party_presolve_t<i_t, f_t>::uncrush_primal_solution(
full_primal = std::move(full_sol.primal);
}

template <typename i_t, typename f_t>
void third_party_presolve_t<i_t, f_t>::crush_primal_solution(
const std::vector<f_t>& original_primal, std::vector<f_t>& reduced_primal) const
{
cuopt_expects(presolver_ == cuopt::linear_programming::presolver_t::Papilo,
error_type_t::RuntimeError,
"Primal crushing is only supported for PaPILO presolve");
cuopt_assert(papilo_post_solve_storage_ != nullptr, "No postsolve storage available");
std::vector<f_t> unused_y, unused_z;
std::vector<f_t> empty_vals;
std::vector<i_t> empty_indices, empty_offsets;
crush_primal_dual_solution(original_primal,
{},
reduced_primal,
unused_y,
{},
unused_z,
empty_vals,
empty_indices,
empty_offsets);
}

/**
* Crush an original-space primal+dual solution into the presolved (reduced) space.
*
* This is the forward counterpart of Papilo's Postsolve::undo(). It replays
* each presolve reduction in forward order to transform variable/dual values,
* then projects onto the surviving columns/rows via origcol/origrow_mapping.
*
* Only two reductions actually transform survivor coordinates:
* kParallelCol — merges x[col1] into x[col2]
* kRowBoundChangeForcedByRow — conditionally transfers y[deleted_row] → y[kept_row]
*/
template <typename i_t, typename f_t>
void third_party_presolve_t<i_t, f_t>::crush_primal_dual_solution(
const std::vector<f_t>& x_original,
const std::vector<f_t>& y_original,
std::vector<f_t>& x_reduced,
std::vector<f_t>& y_reduced,
const std::vector<f_t>& z_original,
std::vector<f_t>& z_reduced,
const std::vector<f_t>& A_values,
const std::vector<i_t>& A_indices,
const std::vector<i_t>& A_offsets) const
{
cuopt_expects(presolver_ == cuopt::linear_programming::presolver_t::Papilo,
error_type_t::RuntimeError,
"Crushing is only supported for PaPILO presolve");
cuopt_assert(papilo_post_solve_storage_ != nullptr, "No postsolve storage available");

const auto& storage = *papilo_post_solve_storage_;
const auto& types = storage.types;
const auto& indices = storage.indices;
const auto& values = storage.values;
const auto& start = storage.start;
const auto& num = storage.num;

cuopt_assert((int)x_original.size() == (int)storage.nColsOriginal, "");

const bool crush_dual = !y_original.empty();
if (crush_dual) { cuopt_assert((int)y_original.size() == (int)storage.nRowsOriginal, ""); }

const bool crush_rc = !z_original.empty() && crush_dual;
if (crush_rc) { cuopt_assert((int)z_original.size() == (int)storage.nColsOriginal, ""); }

std::vector<f_t> x(x_original.begin(), x_original.end());
std::vector<f_t> y(y_original.begin(), y_original.end());
std::vector<f_t> z(z_original.begin(), z_original.end());

// Track current coefficient values for entries modified by kCoefficientChange,
// so repeated changes to the same (row, col) are handled correctly.
std::map<std::pair<int, int>, f_t> coeff_current;

// Look up the current value of a_{row,col}: use coeff_current if it was
// modified by a prior kCoefficientChange, otherwise fall back to the
// original CSR.
auto get_coeff = [&](int row, int col) -> f_t {
auto it = coeff_current.find({row, col});
if (it != coeff_current.end()) return it->second;
for (i_t p = A_offsets[row]; p < A_offsets[row + 1]; ++p) {
if (A_indices[p] == col) return A_values[p];
}
return 0;
};

for (int i = 0; i < (int)types.size(); ++i) {
int first = start[i];

switch (types[i]) {
case ReductionType::kParallelCol: {
// Storage layout: [orig_col1, flags1, orig_col2, flags2, -1]
// [col1lb, col1ub, col2lb, col2ub, col2scale]
int col1 = indices[first];
int col2 = indices[first + 2];
const f_t& scale = values[first + 4];
x[col2] += scale * x[col1];
break;
Comment on lines +973 to +980
Copy link
Copy Markdown

@coderabbitai coderabbitai bot Apr 15, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟠 Major

Fold reduced costs through kParallelCol as well.

This forward replay updates the survivor primal value, but it leaves z in the original basis. After the final projection, any eliminated parallel column contribution is dropped, so z_reduced is wrong whenever dual/reduced-cost crushing hits a parallel-column reduction.

🐛 Suggested fix 
       case ReductionType::kParallelCol: {
         // Storage layout: [orig_col1, flags1, orig_col2, flags2, -1]
         //                 [col1lb,    col1ub, col2lb,    col2ub, col2scale]
         int col1         = indices[first];
         int col2         = indices[first + 2];
         const f_t& scale = values[first + 4];
         x[col2] += scale * x[col1];
+        if (crush_rc) { z[col2] += scale * z[col1]; }
         break;
       }

As per coding guidelines, **/*.{cu,cuh,cpp,hpp,h}: Validate algorithm correctness in optimization logic and ensure variables and constraints are accessed from the correct problem context (original vs presolve vs folded vs postsolve); verify index mapping consistency across problem transformations.

🤖 Prompt for AI Agents 
Verify each finding against the current code and only fix it if needed.

In `@cpp/src/mip_heuristics/presolve/third_party_presolve.cpp` around lines 973 -
980, The kParallelCol case only updates the primal x but fails to fold reduced
costs; update the corresponding dual/reduced-cost vector (z) analogously so
eliminated parallel-column contributions aren't lost: inside the
ReductionType::kParallelCol branch (where col1 = indices[first], col2 =
indices[first+2], scale = values[first+4] and x[col2] += scale * x[col1]) also
perform z[col2] += scale * z[col1] (using the same index mapping and scale),
ensuring you access the z array from the same problem context and respect any
index-mapping helpers used elsewhere in this function.

}

case ReductionType::kRowBoundChangeForcedByRow: {
if (!crush_dual) break;
cuopt_assert(i >= 1 && types[i - 1] == ReductionType::kReasonForRowBoundChangeForcedByRow,
"kRowBoundChangeForcedByRow must be preceded by its reason record");

bool is_lhs = indices[first] == 1;
int row = (int)values[first];

int reason_first = start[i - 1];
int deleted_row = indices[reason_first + 1];
f_t factor = values[reason_first];
cuopt_assert(factor != 0, "parallel row factor must be nonzero");

// Forward rule: if the deleted row carried dual signal that the
// reverse would have attributed to the kept row, transfer it back.
f_t candidate = y[deleted_row] / factor;
bool sign_ok = is_lhs ? num.isGT(candidate, (f_t)0) : num.isLT(candidate, (f_t)0);

if (sign_ok) {
f_t y_old = y[row];
y[row] = candidate;
// Maintain z = c - A^T y: propagate the y change into reduced costs
if (crush_rc) {
f_t delta_y = candidate - y_old;
for (i_t p = A_offsets[row]; p < A_offsets[row + 1]; ++p) {
f_t a = get_coeff(row, A_indices[p]);
z[A_indices[p]] -= delta_y * a;
}
}
}
break;
}

case ReductionType::kCoefficientChange: {
if (!crush_rc) break;
int row = indices[first];
int col = indices[first + 1];
f_t a_new = values[first];
f_t a_old = get_coeff(row, col);
coeff_current[{row, col}] = a_new;
z[col] += (a_old - a_new) * y[row];
break;
}

case ReductionType::kSubstitutedColWithDual: {
// Singleton substitution: column j is expressed via equality row k as
// x_j = (rhs_k - Σ_{l≠j} a_kl·x_l) / a_kj
// This changes the objective for every column l in row k:
// c_red[l] = c_orig[l] - (c_j / a_kj) · a_kl
// Adjust z accordingly: Δz[l] = -(a_kl / a_kj)·z[j] - a_kl·y[k]
if (!crush_rc) break;
int row_k = indices[first]; // equality row (original space)
int row_length = (int)values[first];
// Row coefficients start at first+3
int row_coef_start = first + 3;
// Substituted column index is after the row coefficients
int col_j = indices[row_coef_start + row_length];

// Find a_kj (coefficient of col j in row k)
f_t a_kj = 0;
for (int p = 0; p < row_length; ++p) {
if (indices[row_coef_start + p] == col_j) {
a_kj = values[row_coef_start + p];
break;
}
}
if (a_kj == 0) break; // shouldn't happen

f_t z_j = z[col_j];
f_t y_k = y[row_k];

// Adjust z for each surviving column l in the equality row (l ≠ j)
for (int p = 0; p < row_length; ++p) {
int col_l = indices[row_coef_start + p];
if (col_l == col_j) continue;
f_t a_kl = values[row_coef_start + p];
z[col_l] -= (a_kl / a_kj) * z_j + a_kl * y_k;
}
break;
}

case ReductionType::kFixedCol: // Handled via projection
case ReductionType::kSubstitutedCol: // Col is dropped
case ReductionType::kFixedInfCol: // Col is dropped
case ReductionType::kVarBoundChange: // Noop
case ReductionType::kRedundantRow: // Noop
case ReductionType::kRowBoundChange: // Noop
case ReductionType::kReasonForRowBoundChangeForcedByRow: // Metadata for above
case ReductionType::kSaveRow: // Metadata
case ReductionType::kReducedBoundsCost: // Noop
case ReductionType::kColumnDualValue: // Column reduced-cost only
case ReductionType::kRowDualValue: // Handled via projection
break;
// no default: case to let the compiler yell at us if a new reduction is later introduced
}
}

const auto& col_map = storage.origcol_mapping;
const auto& row_map = storage.origrow_mapping;

// Cancel contributions from removed rows. The original-space z was
// computed as z = c - A^T y over ALL rows. The reduced-space stationarity
// only involves surviving rows, so we must add back the terms from removed
// rows: z[j] += y[i] * a_{i,j} for every removed row i with y[i] != 0.
if (crush_rc) {
std::vector<bool> row_survives((int)storage.nRowsOriginal, false);
for (size_t k = 0; k < row_map.size(); ++k) {
row_survives[row_map[k]] = true;
}
for (int i = 0; i < (int)storage.nRowsOriginal; ++i) {
if (row_survives[i] || y[i] == 0) continue;
for (i_t p = A_offsets[i]; p < A_offsets[i + 1]; ++p) {
z[A_indices[p]] += y[i] * get_coeff(i, A_indices[p]);
}
Comment on lines +1092 to +1096
Copy link
Copy Markdown

@coderabbitai coderabbitai bot Apr 15, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

Avoid exact-zero checks on removed-row duals.

y[i] == 0 is too strict for approximate PDLP/DualSimplex duals. Near-zero removed-row multipliers will still enter this correction path and inject noise into z_reduced. Use the same numeric tolerance machinery you already use elsewhere in this function.

As per coding guidelines, **/*.{cu,cuh,cpp,hpp,h}: Check numerical stability: prevent overflow/underflow, precision loss, division by zero/near-zero, and use epsilon comparisons for floating-point equality checks.

🤖 Prompt for AI Agents 
Verify each finding against the current code and only fix it if needed.

In `@cpp/src/mip_heuristics/presolve/third_party_presolve.cpp` around lines 1092 -
1096, The loop that updates z for removed rows uses an exact check y[i] == 0
which is too strict; replace that check with the function/epsilon-based test
used elsewhere in this function (i.e., treat y as zero when fabs(y[i]) <= the
existing numeric tolerance or by calling the existing isZero/is_near_zero
helper) so that only truly near-zero y[i] are skipped; locate the loop
referencing storage.nRowsOriginal, row_survives, y, A_offsets, A_indices, z and
get_coeff and change the condition to use the same tolerance variable or helper
used elsewhere in this file to avoid injecting noise into z (and therefore
z_reduced).

}
}

x_reduced.resize(col_map.size());
for (size_t k = 0; k < col_map.size(); ++k) {
x_reduced[k] = x[col_map[k]];
}

if (crush_dual) {
y_reduced.resize(row_map.size());
for (size_t k = 0; k < row_map.size(); ++k) {
y_reduced[k] = y[row_map[k]];
}
}

if (crush_rc) {
z_reduced.resize(col_map.size());
for (size_t k = 0; k < col_map.size(); ++k) {
z_reduced[k] = z[col_map[k]];
}
}
}

template <typename i_t, typename f_t>
third_party_presolve_t<i_t, f_t>::~third_party_presolve_t()
{
Expand Down
13 changes: 13 additions & 0 deletions cpp/src/mip_heuristics/presolve/third_party_presolve.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,19 @@ class third_party_presolve_t {

void uncrush_primal_solution(const std::vector<f_t>& reduced_primal,
std::vector<f_t>& full_primal) const;

void crush_primal_solution(const std::vector<f_t>& original_primal,
std::vector<f_t>& reduced_primal) const;

void crush_primal_dual_solution(const std::vector<f_t>& x_original,
const std::vector<f_t>& y_original,
std::vector<f_t>& x_reduced,
std::vector<f_t>& y_reduced,
const std::vector<f_t>& z_original,
std::vector<f_t>& z_reduced,
const std::vector<f_t>& A_values,
const std::vector<i_t>& A_indices,
const std::vector<i_t>& A_offsets) const;
const std::vector<i_t>& get_reduced_to_original_map() const { return reduced_to_original_map_; }
const std::vector<i_t>& get_original_to_reduced_map() const { return original_to_reduced_map_; }

Expand Down
Loading
Loading