Skip to content

Extend the barrier solver to SOCP problems#1051

Closed
yan-zaretskiy wants to merge 94 commits into
NVIDIA:release/26.06from
yan-zaretskiy:socp-barrier-solver
Closed

Extend the barrier solver to SOCP problems#1051
yan-zaretskiy wants to merge 94 commits into
NVIDIA:release/26.06from
yan-zaretskiy:socp-barrier-solver

Conversation

@yan-zaretskiy
Copy link
Copy Markdown

This is a draft PR that extends the 2x2 augmented system formulation of the barrier solver to solve SOCP problems. It supports both explicit cone variables as well as the cone rows in the Ax=b constraint, both of which can be specified in the CBF format. It uses Nesterov-Todd scaling to make the conic diagonal blocks symmetric.

Build the cone-local Jordan-product, scaling, and corrector kernels needed
for the SOCP barrier path before augmented-system integration.
Wire the SOCP cone Hessian updates into the augmented solve path and fold in
the follow-up cleanup that removes now-redundant kernel plumbing.
Tighten the augmented-system SOCP updates so the barrier path converges on the
mixed LP/QP cone cases covered by the new regression tests.
Flatten the remaining SOCP barrier plumbing, remove superseded cone kernels,
and harden presolve so linear columns stay ahead of the trailing cone block.
…hecks

Add row-cone and presolve regression cases around the SOCP barrier path, clean
up the PR-facing test wording, and align the cone layout validation with the
shared infinity convention.

Signed-off-by: Yan Zaretskiy <yzaretskiy@nvidia.com>
@copy-pr-bot
Copy link
Copy Markdown

copy-pr-bot Bot commented Apr 7, 2026

This pull request requires additional validation before any workflows can run on NVIDIA's runners.

Pull request vetters can view their responsibilities here.

Contributors can view more details about this message here.

@rg20 rg20 added feature request New feature or request non-breaking Introduces a non-breaking change labels Apr 23, 2026
@rg20 rg20 added this to the 26.06 milestone Apr 23, 2026
@rg20 rg20 requested review from chris-maes, rg20 and yuwenchen95 April 23, 2026 22:08
// x_j <= 1/a_ij * (rhs - upper_activity_i)
const i_t j = last_free_i;
const f_t a_ij = last_free_coeff_i;
if (a_ij == 0) { continue; }
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Unnecessary as we have removed all 0 coefficients from the matrix.

}
}
if (!(problem.lower[j] == -inf && problem.upper[j] == inf)) {
presolve_info.phase1_bounded_linear_indices.push_back(j);
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

We should not have a variable named phase_1_bounded_linear_indices

}
}

// Barrier presolve phase 2: negate one-sided bounds (-inf < x <= u -> -u <= x < inf).
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Remove "Barrier presolve phase 2" from comment.


// The original problem may have nonzero lower bounds
// 0 != l_j <= x_j <= u_j
// Barrier presolve phase 3: shift nonzero lower bounds to zero (cone/stack columns not shifted).
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Remove "Barrier phase 3". Put back the original comment that was deleted.
Fine to add something about "Cone variables should not be shifted"

}

// Check for empty rows
// Barrier presolve phase 4: remove empty rows and empty linear columns.
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Remove "Barrier presolve phase 4"

if (settings.barrier_presolve && free_variables > 0 && problem.Q.n > 0) {
presolve_info.free_variable_indices.clear();
for (i_t j = 0; j < problem.num_cols; j++) {
// Barrier presolve phase 5: free linear variables — 5a native (QP/SOCP) or 5b v-w split (LP).
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Remove "Barrier presolve phase 5" comment. And "5a, 5b"

for (i_t j = 0; j < linear_cols; j++) {
if (problem.lower[j] == -inf && problem.upper[j] == inf) {
presolve_info.free_variable_indices.push_back(j);
presolve_info.native_free_linear_indices.push_back(j);
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Why are we changing from free_variable_indices to native_free_linear_indices? Let's keep the old name

}
} else if (settings.barrier_presolve && free_variables > 0) {
settings.log.printf(
"Keeping %d native free linear variables for augmented-system barrier (QP/SOCP)\n",
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

"Keeping %d free variables"

"Keeping %d native free linear variables for augmented-system barrier (QP/SOCP)\n",
native_free_count);
} else if (settings.barrier_presolve && !has_cones && free_variables > 0) {
// Phase 5b: x_j = v - w with v, w >= 0 (LP without cones; SOCP/QP use phase 5a).
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Remove "Phase 5" and "phase 5a" from comment. You could delete this comment.

// becomes
// sum_{k != j} c_k x_k + c_j v - c_j w

std::vector<i_t> pair_index(problem.num_cols, -1);
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Why does any of this code need to change. This is for the LP case. Just make sure we aren't in an SOCP or QP and then use this code. I'd suggest reverting these changes.

Q_j[row_starts[partner_row]] = partner_col;
Q_x[row_starts[partner_row]] = qij;
row_starts[partner_row]++;
}
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Please revert all changes from line 1196/1135 to here.

problem.A.i[q] = i;
problem.A.x[q] = -aij;
q++;
csc_matrix_t<i_t, f_t> expanded_A(problem.A.m, num_cols, nnz);
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Can you do this expansion only in the SOCP case?

}

if (!presolve_info.phase1_bounded_linear_indices.empty()) {
settings.log.printf("Post-solve: %d linear column(s) had phase-1 bound tightening (x unchanged)\n",
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

The user should not be expected to know the internal phases of presolve!

Remove this print and this block of code entirely.

const i_t u = free_variable_pairs[k];
const i_t v = free_variable_pairs[k + 1];
input_x[u] -= input_x[v];
remove_partner[v] = true;
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Why is this code needed?

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;
csr_matrix_t<i_t, f_t> A_row(A.m, A.n, A.nnz());
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Revert these changes. Converting to A_row is not needed

if (w_j == 0.0) { continue; }
const f_t du = w_j / bfv.coefficient;
input_y[bfv.constraint] += du;
for (i_t j = 0; j < A.n; j++) {
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Revert these changes.

Comment thread cpp/src/dual_simplex/presolve.hpp
Comment thread cpp/src/dual_simplex/simplex_solver_settings.hpp
Comment thread cpp/src/dual_simplex/simplex_solver_settings.hpp
Comment thread cpp/src/math_optimization/solver_settings.cu
Comment thread cpp/src/pdlp/solve.cu
Comment thread cpp/src/pdlp/solve.cu
Comment thread cpp/src/pdlp/solve.cu
Comment thread cpp/src/pdlp/solve.cu
}
settings.method = method_t::Barrier;
settings.presolver = presolver_t::None;
// Quadratic objective support is minimization-only.
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Is this correct? I thought we supported maximizing with quadratic objectives.

user_problem.Q_values = model.Q_values;

if (model.original_problem_ptr->has_quadratic_constraints()) {
detail::apply_soc_qcmatrix_conversion_for_simplex<i_t, f_t>(
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

for_simplex seems like a bad name here. As we aren't using the simplex solver. Maybe this was picked up because of the simplex namespace? Anyway, better to rename the function to avoid confusion.

Comment thread .claude-plugin/marketplace.json
@@ -61,6 +65,29 @@ namespace cuopt::linear_programming::dual_simplex {
cuda::std::make_tuple(a, b), out, size, cuda::std::multiplies<>{}, stream.value());
}

// out[i] = is_native_free_linear[i] ? 0 : a[i] * b[i]
[[maybe_unused]] static void pairwise_multiply_skip_native_free_linear(
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Nit: any reason not to template this on f_t to avoid having two identical versions?

}
} else if (settings.check_Q) {
// TODO: Check to ensure that Q is positive semi-definite
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Consider adding a warning that check_Q is ignored.

RAFT_CUDA_TRY(cub::DeviceSelect::Flagged(
nullptr,
flag_buffer_size,
d_inv_diag_prime.data(), // Not the actual input but just to allcoate the memory
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Typo: allocate

Comment thread cpp/src/io/mps_parser.cpp
} else if (str == "UI") {
return UpperBoundIntegerVariable;
} else if (str == "SC") {
} else if (str == "SC" || str == "LC") {
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

This diff shouldn't be here, please revert this line.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Comment thread cpp/src/io/mps_parser.cpp

namespace {

/** @brief ceil(log2(x)) for x >= 1; used for sparse-vs-dense CSR heuristic only. */
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

This code is pretty confusing with three code paths(?) for processing CSR. Let's just have one code path until there's a clear reason to have more? LP also needs COO to CSR for the objective, right? let's share the code.

Comment thread cpp/src/io/mps_parser.cpp
csr_result.offsets);
++quadratic_row_id;
const i_t nnz = static_cast<i_t>(block.entries.size());
std::vector<i_t> qc_rows(static_cast<size_t>(nnz));
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

I don't think we need a static_cast<size_t> everywhere, please remove.

Comment thread cpp/src/io/mps_parser.cpp
row_id,
row_names[row_id],
static_cast<char>(row_types[row_id]),
std::span<const f_t>(A_values[row_id].data(), A_values[row_id].size()),
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

I think you can pass A_values (etc) directly, explicit span construction shouldn't be needed.

Comment thread cpp/tests/CMakeLists.txt
"${dejavu_SOURCE_DIR}"
)
# Third-party headers: SYSTEM matches libcuopt and suppresses NVCC sign-conversion noise.
target_include_directories(cuopttestutils SYSTEM PRIVATE "${dejavu_SOURCE_DIR}")
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Unrelated diff? please revert

Apache 2.0 <br>
## Use Case: <br>
Developers and engineers use this skill to capture generalizable learnings from agent interactions and propose updates to existing skills, ensuring continuous improvement of agent guidance without blocking the user's primary task. <br>
Developers and engineers use this skill to continuously improve agent skill files by capturing generalizable learnings from problem-solving interactions and proposing structured updates. <br>
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Please revert all of these

ConfigureTest(DUAL_SIMPLEX_TEST
${CMAKE_CURRENT_SOURCE_DIR}/unit_tests/solve.cpp
${CMAKE_CURRENT_SOURCE_DIR}/unit_tests/solve_barrier.cu
LABELS numopt)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Do not delete the test label

Copy link
Copy Markdown
Contributor

@rg20 rg20 left a comment

Choose a reason for hiding this comment

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

I reviewed only a subset of files. I will continue to review.

std::vector<i_t> linear_indices{};
f_t rhs_value{f_t(0)};
/** Q in COO: parallel arrays, same length. */
std::vector<i_t> quadratic_row_indices{};
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Suggested change
std::vector<i_t> quadratic_row_indices{};
std::vector<i_t> rows{};

lp_problem_t<i_t, f_t>& scaled,
std::vector<f_t>& column_scaling)
std::vector<f_t>& col_scale,
std::vector<f_t>& row_scale)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Looks like row_scaling is not used to scale matrix A?

row_scale.assign(m, 1.0);

// =========================================================================
// Ruiz equilibration for SOCP (and QP) problems
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Can we remove this scaling?

We should do this cleanly on the KKT system.

}

template <typename i_t, typename f_t>
bool validate_barrier_cone_layout(const lp_problem_t<i_t, f_t>& problem,
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

This function does not belong here. Move this to src/barrier

@chris-maes
Copy link
Copy Markdown
Contributor

This PR is moved to #1290 . Please make all comments and changes on that PR.

@chris-maes chris-maes closed this May 23, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

feature request New feature or request non-breaking Introduces a non-breaking change P0

Projects

None yet

Development

Successfully merging this pull request may close these issues.

10 participants