Skip to content

Fix a bug in dual variables and reduced costs when we add implied bounds on free variables #1133

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
chris-maes wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Fix a bug in dual variables and reduced costs when we add implied bounds on free variables #1133

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f540331
Fix a bug in dual variables and reduced costs when we add implied bou…
chris-maes Apr 22, 2026
c29cc5e
Merge branch 'main' into dual_post_solve_bounded_free_vars
chris-maes Apr 28, 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
76 changes: 68 additions & 8 deletions cpp/src/dual_simplex/presolve.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -851,6 +851,12 @@ i_t presolve(const lp_problem_t<i_t, f_t>& original,

i_t removed_free_variables = 0;

// Track which constraint provided each implied bound for dual correction
std::vector<i_t> lower_bound_constraint(problem.num_cols, -1);
std::vector<f_t> lower_bound_coefficient(problem.num_cols, 0.0);
std::vector<i_t> upper_bound_constraint(problem.num_cols, -1);
std::vector<f_t> upper_bound_coefficient(problem.num_cols, 0.0);

if (constraints_to_check.size() > 0) {
// Check if the constraints are feasible
csr_matrix_t<i_t, f_t> Arow(0, 0, 0);
Expand Down Expand Up @@ -928,30 +934,38 @@ i_t presolve(const lp_problem_t<i_t, f_t>& original,
if (lower_inf_i == 1) {
const f_t new_upper = 1.0 / a_ij * (rhs - lower_activity_i);
if (new_upper < max_bound) {
problem.upper[j] = new_upper;
bounded = true;
problem.upper[j] = new_upper;
upper_bound_constraint[j] = i;
upper_bound_coefficient[j] = a_ij;
bounded = true;
}
}
if (upper_inf_i == 1) {
const f_t new_lower = 1.0 / a_ij * (rhs - upper_activity_i);
if (new_lower > -max_bound) {
problem.lower[j] = new_lower;
bounded = true;
problem.lower[j] = new_lower;
lower_bound_constraint[j] = i;
lower_bound_coefficient[j] = a_ij;
bounded = true;
}
}
} else if (a_ij < 0) {
if (lower_inf_i == 1) {
const f_t new_lower = 1.0 / a_ij * (rhs - lower_activity_i);
if (new_lower > -max_bound) {
problem.lower[j] = new_lower;
bounded = true;
problem.lower[j] = new_lower;
lower_bound_constraint[j] = i;
lower_bound_coefficient[j] = a_ij;
bounded = true;
}
}
if (upper_inf_i == 1) {
const f_t new_upper = 1.0 / a_ij * (rhs - upper_activity_i);
if (new_upper < max_bound) {
problem.upper[j] = new_upper;
bounded = true;
problem.upper[j] = new_upper;
upper_bound_constraint[j] = i;
upper_bound_coefficient[j] = a_ij;
bounded = true;
}
}
}
Expand All @@ -973,6 +987,24 @@ i_t presolve(const lp_problem_t<i_t, f_t>& original,
}
}

// Record bounded free variables for dual correction in uncrush.
// After the keep-one-bound logic, each bounded variable has exactly one finite bound.
for (i_t j : current_free_variables) {
i_t bounding_constraint = -1;
f_t bounding_coefficient = 0.0;
if (problem.lower[j] > -inf && lower_bound_constraint[j] != -1) {
bounding_constraint = lower_bound_constraint[j];
bounding_coefficient = lower_bound_coefficient[j];
} else if (problem.upper[j] < inf && upper_bound_constraint[j] != -1) {
bounding_constraint = upper_bound_constraint[j];
bounding_coefficient = upper_bound_coefficient[j];
}
if (bounding_constraint != -1) {
presolve_info.bounded_free_variables.push_back(
{j, bounding_constraint, bounding_coefficient});
}
}

i_t new_free_variables = 0;
for (i_t j = 0; j < problem.num_cols; j++) {
if (problem.lower[j] == -inf && problem.upper[j] == inf) { new_free_variables++; }
Expand Down Expand Up @@ -1562,6 +1594,7 @@ void uncrush_dual_solution(const user_problem_t<i_t, f_t>& user_problem,
template <typename i_t, typename f_t>
void uncrush_solution(const presolve_info_t<i_t, f_t>& presolve_info,
const simplex_solver_settings_t<i_t, f_t>& settings,
const lp_problem_t<i_t, f_t>& original_problem,
const std::vector<f_t>& crushed_x,
const std::vector<f_t>& crushed_y,
const std::vector<f_t>& crushed_z,
Expand Down Expand Up @@ -1711,6 +1744,32 @@ void uncrush_solution(const presolve_info_t<i_t, f_t>& presolve_info,
}
}

// Dual correction for originally-free variables that received implied bounds.
// Barrier produced (u, w) with w_j != 0 satisfying A^T u + w = c + Qx.
// We need corrected (y, z) with z_j = 0: set du = w_j / a_{i*,j}, then
// y_{i*} += du and z_k -= a_{i*,k} * du for all k in constraint i*.
if (!presolve_info.bounded_free_variables.empty()) {
settings.log.printf("Post-solve: Correcting duals for %d bounded free variables\n",
static_cast<i_t>(presolve_info.bounded_free_variables.size()));
const csc_matrix_t<i_t, f_t>& A = original_problem.A;
for (const auto& bfv : presolve_info.bounded_free_variables) {
const f_t w_j = input_z[bfv.variable];
if (w_j == 0.0) { continue; }
const f_t du = w_j / bfv.coefficient;
input_y[bfv.constraint] += du;
Comment on lines +1751 to +1759
Copy link
Copy Markdown

@coderabbitai coderabbitai Bot Apr 22, 2026 •
edited
Loading

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟠 Major

Guard the dual update against tiny a_{i*,j} and tiny w_j.

This divides by bfv.coefficient directly and only skips exact-zero w_j. On ill-conditioned rows, a near-zero coefficient will blow up du, and exact floating-point equality here is too strict for postsolve residuals.

As per coding guidelines, "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/dual_simplex/presolve.cpp` around lines 1751 - 1759, The loop
updating input_y for presolve_info.bounded_free_variables divides by
bfv.coefficient and only skips when w_j == 0.0, which is numerically unsafe;
change the guard to check for near-zero values using an epsilon (e.g.,
fabs(bfv.coefficient) < eps or fabs(w_j) < eps) and skip or clamp the update
when either the coefficient or w_j is too small to avoid huge du; update the
conditional around computing du (referencing
presolve_info.bounded_free_variables, bfv.coefficient, input_z, input_y) to use
an absolute epsilon comparison and optionally log or count skipped
tiny-coefficient cases instead of performing the division.

Copy link
Copy Markdown
Collaborator

@rgsl888prabhu rgsl888prabhu Apr 23, 2026

Choose a reason for hiding this comment

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

@chris-maes is this a good suggestion ?

Copy link
Copy Markdown

@coderabbitai coderabbitai Bot Apr 23, 2026

Choose a reason for hiding this comment

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

Seems like the humans are having a chat. I'll hop back into my burrow for now. If you need me again, just tag @coderabbitai in a new comment, and I'll come hopping out!

Copy link
Copy Markdown
Contributor Author

@chris-maes chris-maes Apr 23, 2026

Choose a reason for hiding this comment

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

We remove small coefficients in A in presolve. I'm not sure if this is necessary.

for (i_t j = 0; j < A.n; j++) {
const i_t col_start = A.col_start[j];
const i_t col_end = A.col_start[j + 1];
for (i_t p = col_start; p < col_end; p++) {
if (A.i[p] == bfv.constraint) {
input_z[j] -= A.x[p] * du;
break;
}
}
}
}
}

assert(uncrushed_x.size() == input_x.size());
assert(uncrushed_y.size() == input_y.size());
assert(uncrushed_z.size() == input_z.size());
Expand Down Expand Up @@ -1769,6 +1828,7 @@ template void uncrush_dual_solution<int, double>(const user_problem_t<int, doubl

template void uncrush_solution<int, double>(const presolve_info_t<int, double>& presolve_info,
const simplex_solver_settings_t<int, double>& settings,
const lp_problem_t<int, double>& original_problem,
const std::vector<double>& crushed_x,
const std::vector<double>& crushed_y,
const std::vector<double>& crushed_z,
Expand Down
13 changes: 13 additions & 0 deletions cpp/src/dual_simplex/presolve.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -129,6 +129,15 @@ struct folding_info_t {
bool is_folded;
};

// Free variable that received an implied bound during presolve.
// Stores the bounding constraint and coefficient for dual correction in uncrush.
template <typename i_t, typename f_t>
struct bounded_free_var_t {
i_t variable; // j: the originally-free variable
i_t constraint; // i*: the constraint that implied the bound
f_t coefficient; // a_{i*,j}: the coefficient of x_j in constraint i*
};

template <typename i_t, typename f_t>
struct presolve_info_t {
// indices of variables in the original problem that remain in the presolved problem
Expand All @@ -153,6 +162,9 @@ struct presolve_info_t {

// Variables that were negated to handle -inf < x_j <= u_j
std::vector<i_t> negated_variables;

// Originally-free variables that received implied bounds, with the constraint used
std::vector<bounded_free_var_t<i_t, f_t>> bounded_free_variables;
};

template <typename i_t, typename f_t>
Expand Down Expand Up @@ -241,6 +253,7 @@ void uncrush_dual_solution(const user_problem_t<i_t, f_t>& user_problem,
template <typename i_t, typename f_t>
void uncrush_solution(const presolve_info_t<i_t, f_t>& presolve_info,
const simplex_solver_settings_t<i_t, f_t>& settings,
const lp_problem_t<i_t, f_t>& original_problem,
const std::vector<f_t>& crushed_x,
const std::vector<f_t>& crushed_y,
const std::vector<f_t>& crushed_z,
Expand Down
2 changes: 2 additions & 0 deletions cpp/src/dual_simplex/solve.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -297,6 +297,7 @@ lp_status_t solve_linear_program_with_advanced_basis(
unscale_solution<i_t, f_t>(column_scales, solution.x, solution.z, unscaled_x, unscaled_z);
uncrush_solution(presolve_info,
settings,
original_lp,
unscaled_x,
solution.y,
unscaled_z,
Expand Down Expand Up @@ -439,6 +440,7 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t<i_t, f_t>& us
// Undo presolve
uncrush_solution(presolve_info,
barrier_settings,
original_lp,
unscaled_x,
barrier_solution.y,
unscaled_z,
Expand Down
5 changes: 5 additions & 0 deletions cpp/tests/dual_simplex/unit_tests/solve.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,11 +17,13 @@
#include <dual_simplex/user_problem.hpp>

#include <mps_parser/parser.hpp>
#include <utilities/logger.hpp>

namespace cuopt::linear_programming::dual_simplex::test {

TEST(dual_simplex, chess_set)
{
cuopt::init_logger_t log("", true);
namespace dual_simplex = cuopt::linear_programming::dual_simplex;
raft::handle_t handle{};
dual_simplex::user_problem_t<int, double> user_problem(&handle);
Expand Down Expand Up @@ -95,6 +97,7 @@ TEST(dual_simplex, chess_set)

TEST(dual_simplex, burglar)
{
cuopt::init_logger_t log("", true);
constexpr int num_items = 8;
constexpr double max_weight = 102;

Expand Down Expand Up @@ -169,6 +172,7 @@ TEST(dual_simplex, burglar)

TEST(dual_simplex, empty_columns)
{
cuopt::init_logger_t log("", true);
// Same as burglar problem above but with an empty column inserted
constexpr int num_items = 9;
constexpr double max_weight = 102;
Expand Down Expand Up @@ -257,6 +261,7 @@ TEST(dual_simplex, empty_columns)

TEST(dual_simplex, dual_variable_greater_than)
{
cuopt::init_logger_t log("", true);
// minimize 3*x0 + 2 * x1
// subject to x0 + x1 >= 1
// x0 + 2x1 >= 3
Expand Down
42 changes: 42 additions & 0 deletions cpp/tests/dual_simplex/unit_tests/solve_barrier.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,9 +8,13 @@
#include <cstdio>

#include <utilities/common_utils.hpp>
#include <utilities/copy_helpers.hpp>

#include <gtest/gtest.h>

#include <cuopt/linear_programming/constants.h>
#include <cuopt/linear_programming/pdlp/solver_settings.hpp>
#include <cuopt/linear_programming/solve.hpp>
#include <dual_simplex/presolve.hpp>
#include <dual_simplex/solve.hpp>
#include <dual_simplex/tic_toc.hpp>
Expand All @@ -20,6 +24,7 @@
#include <raft/core/cusparse_macros.hpp>

#include <mps_parser/parser.hpp>
#include <utilities/logger.hpp>

namespace cuopt::linear_programming::dual_simplex::test {

Expand All @@ -35,6 +40,7 @@ static void init_handler(const raft::handle_t* handle_ptr)

TEST(barrier, chess_set)
{
cuopt::init_logger_t log("", true);
namespace dual_simplex = cuopt::linear_programming::dual_simplex;
raft::handle_t handle{};
init_handler(&handle);
Expand Down Expand Up @@ -104,6 +110,7 @@ TEST(barrier, chess_set)

TEST(barrier, dual_variable_greater_than)
{
cuopt::init_logger_t log("", true);
// minimize 3*x0 + 2 * x1
// subject to x0 + x1 >= 1
// x0 + 2x1 >= 3
Expand Down Expand Up @@ -174,4 +181,39 @@ TEST(barrier, dual_variable_greater_than)
EXPECT_NEAR(solution.z[1], 0.0, 1e-5);
}

TEST(barrier, min_x_squared_free_variable_dual_correction)
{
// minimize x^2 (Q = [2.0], so 0.5 * x^T Q x = x^2)
// subject to x >= 1
// x is free
//
// Optimal: x = 1, obj = 1, y[0] = 2, z[0] = 0
// This tests the dual correction for originally-free variables that
// received implied bounds during presolve.

const raft::handle_t handle{};
init_handler(&handle);

auto path = cuopt::test::get_rapids_dataset_root_dir() +
"/quadratic_programming/min_x_squared.mps";
auto mps_data = cuopt::mps_parser::parse_mps<int, double>(path);

auto settings = cuopt::linear_programming::pdlp_solver_settings_t<int, double>{};

auto solution = cuopt::linear_programming::solve_lp(&handle, mps_data, settings);

EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERMINATION_STATUS_OPTIMAL);

auto h_x = cuopt::host_copy(solution.get_primal_solution(), handle.get_stream());
auto h_y = cuopt::host_copy(solution.get_dual_solution(), handle.get_stream());
auto h_z = cuopt::host_copy(solution.get_reduced_cost(), handle.get_stream());

printf("x %e y %e z %e\n", h_x[0], h_y[0], h_z[0]);

const double tol = 1e-5;
EXPECT_NEAR(h_x[0], 1.0, tol);
EXPECT_NEAR(h_y[0], 2.0, tol);
EXPECT_NEAR(h_z[0], 0.0, tol);
}

} // namespace cuopt::linear_programming::dual_simplex::test
13 changes: 13 additions & 0 deletions datasets/quadratic_programming/min_x_squared.mps
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
NAME min_x_squared
ROWS
N obj
G c1
COLUMNS
x c1 1
RHS
RHS_V c1 1
BOUNDS
FR BOUND x
QUADOBJ
x x 2
ENDATA
Loading