Problem

In the dense backend with box constraints, update(..., update_preconditioner = false) (the default since #250) returns box multipliers scaled by a factor which compounds on every call. Re-solving an unchanged model gives a growing dual error while results.info.status stays PROXQP_SOLVED. The primal solution is unaffected, so the error is invisible unless the stationarity residual is checked explicitly.

// min 0.5 x' H x + g' x   s.t. -1 <= x <= 1,   H = diag(1, 100), g = (-3, -300)
// The data never changes, so every solve must return the same answer.
proxsuite::proxqp::dense::QP<double> qp(2, 0, 0, /* box_constraints = */ true);
qp.init(H, g, A, b, C, l, u, lb, ub);
for (int k = 0; k < 5; ++k) {
  if (k) qp.update(H, g, A, b, C, l, u, lb, ub, /* update_preconditioner = */ false);
  qp.solve();
  // print qp.work.i_scaled, qp.results.z, ||H x + g + z||
}

The exact multipliers are z = (2, 200):

solve 1  i_scaled = [1.000000 0.991046]  z = [2.00000 200.00000]  ||Hx+g+z|| = 1.023e-12
solve 2  i_scaled = [1.000000 0.982172]  z = [2.00000 201.80701]  ||Hx+g+z|| = 1.807e+00
solve 3  i_scaled = [1.000000 0.973377]  z = [2.00000 203.63034]  ||Hx+g+z|| = 3.630e+00
solve 4  i_scaled = [1.000000 0.964662]  z = [2.00000 205.47015]  ||Hx+g+z|| = 5.470e+00
solve 5  i_scaled = [1.000000 0.956024]  z = [2.00000 207.32659]  ||Hx+g+z|| = 7.327e+00

z(1) follows 200 / 0.991046^k, i.e. the returned multiplier is the true one divided by the
accumulated i_scaled drift.

Cause

dense/preconditioner/ruiz.hpp, scale_qp_in_place: the branch taken when execute_preconditioner == false scales i_scaled without resetting it first, while the true branch resets it via i_scaled.setOnes().

i_scaled is also the only scaled quantity that dense/helpers.hpp::setup does not refresh before scaling. The matrices/vectors H_scaled, g_scaled, A_scaled, b_scaled, C_scaled, l/u_scaled and l/u_box_scaled are all re-copied from the unscaled qpmodel on every update, so they start clean. i_scaled has no counterpart in the model to be copied from, and so carries the previous scaling into the next one.

After k updates it holds (delta.head * delta.tail)^k instead of delta.head * delta.tail.