Inaccurate solution for simple QP problem

[rienq]

Hi @goulart-paul, I was noticing some strange results with Clarabel (using v0.7.0) and I was able to reproduce it in the form of a simple QP with 2 variables, 1 equality constraint and 2 inequality constraints. The code in C++ looks like this:

  MatrixXd P_dense(2, 2);
  P_dense << 2., 0., 0., 4.;

  SparseMatrix<double> P = P_dense.sparseView();
  P.makeCompressed();

  Vector<double, 2> q = {-2., 0.};

  MatrixXd A_dense(3, 2);
  A_dense << 1., 1.,
    -1., 0.,
    0., -1.;

  SparseMatrix<double> A = A_dense.sparseView();
  A.makeCompressed();

  Vector<double, 3> b = {1., 0., 0.};

  vector<SupportedConeT<double>> cones{
    ZeroConeT<double>(1),
    NonnegativeConeT<double>(2),
  };

  DefaultSettings<double> settings = DefaultSettings<double>::default_settings();

  DefaultSolver<double> solver(P, q, A, b, cones, settings);

  solver.solve();
  DefaultSolution<double> solution = solver.solution();

The correct solution should be (1, 0), but the solution that is provided has an error of 1e-4 even though Clarabel reports primal and dual feasibility of 1e-12. Obtained solution from Clarabel = 0.999967 3.31364e-05

I attached the complete output as a txt file. Do you have any idea what is going wrong for this problem? simple_qp_output.txt

[rienq]

Update: I just tried the latest version (v0.7.1) in Python and I got a similar answer [0.9999668636034095, 3.313639612649994e-05]

P = sparse.csc_matrix([[2., 0.], [0., 4.]])
P = sparse.triu(P).tocsc()

q = np.array([-2., 0.])

A = sparse.csc_matrix(
    [[1., 1.],        # <-- LHS of equality constraint (lower bound)
     [-1.,  0.],       # <-- LHS of inequality constraint (lower bound)
     [0., -1.]])       # <-- LHS of inequality constraint (lower bound)

b = np.array([1., 0., 0.])

cones = [clarabel.ZeroConeT(1), clarabel.NonnegativeConeT(2)]
settings = clarabel.DefaultSettings()

solver = clarabel.DefaultSolver(P, q, A, b, cones, settings)
solution = solver.solve()
print(
    f"Solver terminated with solution"
    f"{dict(s=solution.s, x=solution.x, z=solution.z)}"
)

[goulart-paul]

I think that the issue is only in the interpretation of the error tolerances.

The primal and dual feasibility tolerances are checked against the residuals of the KKT conditions. In this case when I manually compute the dual residual $||Px + q + A^Tz||$ I get an error with very small inf-norm as expected. Likewise for the primal residual $||Ax + s - b||$. These residuals are what is reported in the solver output dres and pres columns.

The primal objective function at the solution produced by the solver is also really close, i.e. $-1 + (3\times 10^{-9})$, with similar accuracy for the dual objective value (true value is $-1$). The difference between them is reported in the gap column and shows high accuracy.

The error you are computing is the distance between the (approximate) optimizer produced by the solver and the true optimizer. This is not one of the stopping criteria, and indeed really can't be because there is no obvious way to compute this distance when the true optimizer is not known.

Normally you would not expect to get a substantial inaccuracy in the optimizer like that anyway. The reason it happens in this case is because the unconstrained minimizer for the problem happens to be at the boundary of the feasible set. The solution is therefore quite insensitive to small variations around that point, and the solver has produced a solution that is at a nearby point on the boundary with very similar objective value.

Part of the issue is perhaps also one of method. A simplex based solver would be more likely to produce a very high accuracy optimizer for this problem because the true solution is right on a vertex.

A possible fix for this would be for us to implement a "polishing" post-processing step as we did in OSQP. That would only be directly applicable to LP / QP problems, although I am aware that there is at least some work out there on similar polishing methods for problems with other constraint types.

[rienq]

Thank you so much for the quick response, @goulart-paul ! I agree with everything you said, and I had a suspicion it was related to the fact that the constraint for the second variable is only weakly active. The product of the slack variable and the multiplier is less than 1e-8, but that means that both the slack variable and the Lagrange multiplier are around 1e-4 even though they should both be zero in the optimal solution. The solution becomes more accurate if I reduce the tolerance value for the duality gap, for example, I get [0.999999968991507, 3.1e-08] for an absolute and relative gap tolerance of 1e-14. But I guess this could lead to numerical difficulties for other problems.

[rienq]

Feel free to close this issue. Or you can leave it open if you think that such a "polishing" post-processing step could be implemented also for Clarabel.