[BUG] get_lp_stats() reports garbage dual residual and objective for QPs and SOCPs

[mlubin]

Describe the bug

For quadratic problems (quadratic objective or quadratic constraints), which are solved via the barrier path, the convergence statistics surfaced through the solution — l2_dual_residual and gap (i.e. Python's get_lp_stats()["dual_residual"] / ["gap"]) — are wrong:

• l2_dual_residual omits the quadratic term. At a converged QP optimum it reports ‖Aᵀy + z − c‖ (the linear stationarity residual) instead of ‖Qx + c − Aᵀy − z‖. The Qx term is dropped, so a fully optimal solution reports a large dual residual.
• gap is a hardcoded 0.0, and primal_objective / dual_objective are both set to the same user_objective, regardless of the actual duality/complementarity gap the barrier computed.

l2_primal_residual and nb_iterations are correct. (For SOCP/quadratic constraints, dual variables are additionally filled with NaN, so any dual-derived stat is meaningless.)

Steps/Code to reproduce bug

// Expected behavior for a converged QP: the reported convergence stats should
// reflect the actual optimality of the solution. This test encodes that
// expectation and currently FAILS, demonstrating the bug.
//
// Minimize x1^2 + 4 x2^2 - 8 x1 - 16 x2 = (1/2) x^T H x + c^T x,
//   H = diag(2, 8), c = (-8, -16)
// s.t. x1 + x2 >= 5,  0 <= x1,x2 <= 10.
// Unconstrained optimum (4, 2): H x + c = (0,0), constraint slack (y=0),
// x interior to bounds (z=0). True stationarity residual Q x + c - A^T y - z ~ 0.
TEST(lp_parser_solve, qp_dual_residual_is_wrong)
{
  raft::handle_t handle;
  auto problem  = io::read_lp_from_string<int, double>(R"LP(
Minimize
  obj: -8 x1 - 16 x2 + [ 2 x1 ^ 2 + 8 x2 ^ 2 ] / 2
Subject To
  c1: x1 + x2 >= 5
Bounds
  0 <= x1 <= 10
  0 <= x2 <= 10
End
)LP");
  auto settings = pdlp_solver_settings_t<int, double>();
  auto solution = solve_lp(&handle, problem, settings);

  ASSERT_EQ(solution.get_termination_status(), pdlp_termination_status_t::Optimal);
  EXPECT_NEAR(solution.get_objective_value(), -32.0, 1e-4);

  auto info = solution.get_additional_termination_information();

  // Primal side is correct: the solution is genuinely feasible and converged.
  EXPECT_NEAR(info.l2_primal_residual, 0.0, 1e-4);

  // Expected: at the optimum the true stationarity residual is ~0.
  // BUG: reports ||A^T y + z - c|| = ||c|| = 16, dropping the Q x term.
  EXPECT_NEAR(info.l2_dual_residual, 0.0, 1e-4);

  // Expected: gap is the real primal-dual gap and the objectives are
  // independently computed bounds.
  // BUG: gap is hardcoded 0.0 and dual_objective is a copy of primal.
  EXPECT_NE(info.primal_objective, info.dual_objective)
    << "primal and dual objectives are set to the same value; gap is not computed";
}

[atiaomar1978-hub]

I reproduced the reported behavior analytically using the same QP as in the issue (no GPU run yet). The numbers match the bug description: correct x and objective, but wrong dual residual and gap in the surfaced stats.

Problem

Minimize

(1/2)(2·x1² + 8·x2²) − 8·x1 − 16·x2 = x1² + 4·x2² − 8·x1 − 16·x2

Same as (1/2) xᵀ H x + cᵀ x with H = diag(2, 8) and c = (−8, −16).

Subject to: x1 + x2 ≥ 5, 0 ≤ x1, x2 ≤ 10.

Optimum

Stationarity H x + c = 0 ⇒ x1 = 4, x2 = 2.

• Constraint slack: 4 + 2 − 5 = 1 > 0 (inactive).
• Both variables are strictly inside [0, 10].

So x* = (4, 2) and f(x*) = −32.

Dual residual (what it should be)

QP stationarity at a KKT point:

Q x + c − Aᵀ y − z = 0 (with Q = H)

At x* = (4, 2):

Q x* + c = (8, 16) + (−8, −16) = (0, 0)

With inactive constraint and interior bounds, y = 0 and z = 0 is consistent, so the true stationarity residual is 0. After an optimal solve, get_lp_stats()["dual_residual"] / l2_dual_residual should be ≈ 0, not large.

What the barrier export does (matches the report)

In barrier.cu (to_solution), the dual residual is built from z − c + Aᵀ y and does not include Q x. The value stored in l2_dual_residual uses an L∞ norm (despite the field name).

At the same optimum with y = 0, z = 0:

z − c + Aᵀ y = −c ⇒ ‖−c‖∞ = max(8, 16) = 16

So even when x and the objective are correct, the surfaced dual residual can be ≈ 16, matching the failing assertion EXPECT_NEAR(info.l2_dual_residual, 0.0, 1e-4).

In solve.cu, gap = 0.0 and dual_objective = primal_objective = user_objective on this path, which matches the gap / objective part of the report.

Summary

Quantity	At optimum (4, 2)	Issue expectation	Analytical “buggy stats”
Objective	−32	−32	−32
Primal feasibility	0	≈ 0	≈ 0
Stationarity / dual residual	0	≈ 0	≈ 16 (Qx term omitted)
Gap	meaningful	not identically 0	0 (hardcoded)

I can run the C++ test / Python get_lp_stats() on a CUDA machine next for an end-to-end check; this analytical case already matches the root cause described above.

[mlubin]

@atiaomar1978-hub your comment hasn't added any information to the issue that we didn't already know or couldn't easily ask an agent for.

[atiaomar1978-hub]

@mlubin Thanks for the feedback — understood.
My comment was only meant to clarify that without access to the specific GPU, I can’t reproduce the behavior directly, even though the underlying issue is mathematical. I’ll keep future comments focused strictly on new technical information that can help move the issue forward.
When I buy NVIDIA GPU , I will do full simulation/reproducing for the issue. Sorry for that.