different result in c++ and matlab

[baharehhj]

Hello,

I have implemented the "quadprog" function in MATLAB using the library you provided. I expect that my C++ code will yield the same answer as the MATLAB function. I anticipate that the solution in C++ will match the 'x' variable in MATLAB, but it is different. However, the optimal objective is close to the 'fval' value in MATLAB. It is important to have same answer as matlab. Am I doing something wrong? I appreciate any help.

matlab code

 `[x,fval,eflag,output] = quadprog(H,f,[],[],[],[],lb,ub,[],options);`

part of expected answer: fval:-65600.5031027019

 x:
5.8216955543117365e+02
5.7973844230561019e+02
5.8027260369019268e+02
5.8568208435242616e+02
5.8418146146394179e+02
5.8550733274921367e+02
5.7973844230561019e+02
5.7973844230561019e+02
6.3045238512548747e+02
6.5479734584311518e+02
6.5513552885269223e+02
5.9150783180361020e+02
5.9061530879065367e+02
5.8687045068447662e+02
5.8099606779562760e+02
5.7973844230561019e+02
5.7973844230561019e+02
1.8728350720004594e+03
2.6054743225167676e+03
2.4078617399097943e+03
2.2801859321432567e+03
1.0711290870444805e+03
5.9734933802758121e+02
5.9550936224737495e+02
5.8435715656327272e+02
5.7973844230561019e+02

c++ code

OsqpEigen::Solver solver;
//Eigen::VectorXd lb;
Eigen::VectorXd lb(721);
lb.setOnes();
//lb = lb * 500;
Eigen::MatrixXd denseMatrix = Eigen::MatrixXd::Identity(721, 721);
Eigen::SparseMatrix<double> diagonalMatrix(721, 721);
diagonalMatrix.reserve(Eigen::VectorXi::Constant(721, 1));

for (int i = 0; i < 721; ++i) {
    diagonalMatrix.insert(i, i) = 1;
}
diagonalMatrix.makeCompressed();

Eigen::VectorXd ub = Eigen::VectorXd::Constant(721, std::numeric_limits<double>::infinity());

solver.data()->setNumberOfVariables(H.rows());
solver.data()->setNumberOfConstraints(721);
solver.data()->setLinearConstraintsMatrix(diagonalMatrix);
solver.data()->setHessianMatrix(H);
solver.data()->setGradient(fval);
solver.data()->setLowerBound(lb);
solver.data()->setUpperBound(ub);

std::cout << 1;
if (!solver.initSolver()) {
    std::cerr << "Error initializing solver." << std::endl;
    return 1;
}

solver.settings()->setVerbosity(true);
//solver.solveProblem() == OsqpEigen::ErrorExitFlag::NoError;
if (solver.solveProblem() != OsqpEigen::ErrorExitFlag::NoError)
    return 1;

Eigen::Matrix<c_float, 721, 1> solution = solver.getSolution();

std::cout << "Solution:" << std::endl;
std::cout << solution << std::endl;

answer in terminal

-----------------------------------------------------------------
           OSQP v0.6.3  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2021
-----------------------------------------------------------------
problem:  variables n = 721, constraints m = 721
          nnz(P) + nnz(A) = 149524
settings: linear system solver = qdldl,
          eps_abs = 1.0e-03, eps_rel = 1.0e-03,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 4000
          check_termination: on (interval 25),
          scaling: on, scaled_termination: off
          warm start: on, polish: off, time_limit: off

iter  objective    pri res    dua res    rho        time
   1  -5.7952e+02   1.00e+00   4.87e-01   1.00e-01   5.83e-02s
  75  -6.5601e+04   1.82e-12   1.63e-05   1.00e-06   1.20e-01s

status:               solved
number of iterations: 75
optimal objective:    -65600.9526
run time:             1.21e-01s
optimal rho estimate: 1.00e-06

 Solution:

4.78641 1.06484 27.2492 20.226 13.5775 10.0481 1.06484 1.06484 106.544 204.121 151.687 45.4514 23.29 14.5737 2.85424 1.06484 1.06484 1651.04 2681.07 2667.71 2488.53 746.985 47.8401 27.2042 8.38749 1.06484 1.06484 1025.8 2629.15 2638.96 2218.94 2187.01 2092.31 71.3086 23.6676 11.7751 1.06484 1.06484 525.9 2690.49 2671.57 2533.4 2345.48 2101.57

[S-Dafarra]

Hi @baharehhj, thanks for the issue.

First of all, let me remark that this repo is only a C++ interface to the ospq solver. So for a more dedicated support, I would suggest you to open a ticket in the upstream repo: https://github.com/osqp/osqp

Nonetheless, I would suggest you to check the rank of H. If the Hessian is singular, there could be infinite minima. This might explain why the optimal cost value is very close in the two cases but with very different values for the variables.

[baharehhj]

Hi, @S-Dafarra, thank you for your reply.

I apologize for the delay. I will consider asking it in the repository you recommended.

I will check if the Hessian matrix is singular. I corrected the lower bound constraint but I still did not get the same answer as matlab and it is important for me to get the same answer.

Thanks again for your assistance.

[baharehhj]

disp(options); quadprog options:

Options used by current Algorithm ('interior-point-convex'): (Other available algorithms: 'active-set', 'trust-region-reflective')

Set properties: No options set.

Default properties: Algorithm: 'interior-point-convex' ConstraintTolerance: 1.0000e-08 Display: 'final' LinearSolver: 'auto' MaxIterations: 200 OptimalityTolerance: 1.0000e-08 StepTolerance: 1.0000e-12

these are setting for guadprog function in matlab