Get Unstable Answer When Solve a simple QP using Update Matrix Functions

[little-bird-go]

I don't know the problem is caused by the sparse matrix or the algorithm, hoping get an explaination. I run the code in ubuntu18.04 with the latest version of osqp and osqp-eigen. The whole problem and code is in following file. issure.md

[traversaro]

Can you post the issue directly as text, instead of posting it as an image and then a markdown attached to the issue? Thanks!

Furthermore, does this issue also happens if you specifying the problem for example in Python? In that case, the problem may not be in the osqp-eigen package, but rather in the osqp library itself: https://github.com/osqp .

[little-bird-go]

The codes that are not show in above issures are directly below:

• The code with wrong answer

#include <Eigen/Dense>
#include <OsqpEigen/OsqpEigen.h>
#include <iostream>

using namespace std;

int main(void)
{
    // construct matrix
    Eigen::SparseMatrix<double> H(2,2);
    Eigen::VectorXd f(2);
    Eigen::SparseMatrix<double> A(3,2);
    Eigen::VectorXd lb(3);
    Eigen::VectorXd ub(3);

    // initialize matrix
    H.insert(0,0) = 4.0;
    H.insert(0,1) = 1.0; 
    H.insert(1,0) = 1.0;
    H.insert(1,1) = 2.0;

    f << 1.0, 1.0;

    A.insert(0,0) = 1.0;
    A.insert(0,1) = 1.0;
    A.insert(1,0) = 1.0;
    A.insert(1,1) = 0.0;
    A.insert(2,0) = 0.0;
    A.insert(2,1) = 1.0;

    lb << 1.0, 0.0, 0.0;
    ub << 1.0, 0.7, 0.7;

    // create solver
    OsqpEigen::Solver solver;

    //settings
    solver.settings()->setVerbosity(true);
    solver.settings()->setWarmStart(true);

    // set the dimension
    solver.data()->setNumberOfVariables(2);
    solver.data()->setNumberOfConstraints(3);

    // check before solve
    if(!solver.data()->setHessianMatrix(H)){
        cout << "H matrix error!" << endl;
        return 1;
    } 
    if(!solver.data()->setGradient(f)){
        cout << "f matrix error!" << endl;
        return 1;
    }
    if(!solver.data()->setLinearConstraintsMatrix(A)){
        cout << "A matrix error!" << endl;
        return 1;
    }
    if(!solver.data()->setLowerBound(lb)){
        cout << "lb matrix error!" << endl;
        return 1;
    }
    if(!solver.data()->setUpperBound(ub)){
        cout << "ub matrix error!" << endl;
        return 1;
    }

    // initial the solver
    if(!solver.initSolver()){
        cout << "solver initial error!" << endl;
        return 1;
    }

    // solve the problem
    auto err_code = solver.solveProblem();
    if(err_code != OsqpEigen::ErrorExitFlag::NoError){
        cout << "solve err: " << static_cast<int>(err_code) << endl;
        return 1;
    }

    Eigen::Vector2d solution = solver.getSolution();
    double obj = solver.getObjValue();

    cout << "the solution is: " << endl;
    cout << solution << endl;
    cout << "the obj value is: " << obj << endl;

    ///////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////
    // update H and A
    Eigen::SparseMatrix<double> H_update(2,2);
    Eigen::SparseMatrix<double> A_update(3,2);

    H_update.insert(0,0) = 5.0;
    H_update.insert(0,1) = 1.5; 
    H_update.insert(1,0) = 1.5;
    H_update.insert(1,1) = 1.0;

    A_update.insert(0,0) = 1.2;
    A_update.insert(0,1) = 1.1;
    A_update.insert(1,0) = 1.5;
    // A_update.insert(1,1) = 0.0;
    // A_update.insert(2,0) = 0.0;
    A_update.insert(2,1) = 0.8;

    cout << A_update << endl;

    solver.updateHessianMatrix(H_update);
    solver.updateLinearConstraintsMatrix(A_update);

    err_code = solver.solveProblem();
    if(err_code != OsqpEigen::ErrorExitFlag::NoError){
        cout << "solve err: " << static_cast<int>(err_code) << endl;
        return 1;
    }

    solution = solver.getSolution();
    obj = solver.getObjValue();

    cout << "the solution after updating is: " << endl;
    cout << solution << endl;
    cout << "the obj value after updating is: " << obj << endl;
    return 0;
}

• The code with right answer

#include <Eigen/Dense>
#include <OsqpEigen/OsqpEigen.h>
#include <iostream>

using namespace std;

int main(void)
{
    // construct matrix
    Eigen::SparseMatrix<double> H(2,2);
    Eigen::VectorXd f(2);
    Eigen::SparseMatrix<double> A(3,2);
    Eigen::VectorXd lb(3);
    Eigen::VectorXd ub(3);

    // initialize matrix
    H.insert(0,0) = 4.0;
    H.insert(0,1) = 1.0; 
    H.insert(1,0) = 1.0;
    H.insert(1,1) = 2.0;

    f << 1.0, 1.0;

    A.insert(0,0) = 1.0;
    A.insert(0,1) = 1.0;
    A.insert(1,0) = 1.0;
    A.insert(1,1) = 0.0;
    A.insert(2,0) = 0.0;
    A.insert(2,1) = 1.0;

    lb << 1.0, 0.0, 0.0;
    ub << 1.0, 0.7, 0.7;

    // create solver
    OsqpEigen::Solver solver;

    //settings
    solver.settings()->setVerbosity(true);
    solver.settings()->setWarmStart(true);

    // set the dimension
    solver.data()->setNumberOfVariables(2);
    solver.data()->setNumberOfConstraints(3);

    // check before solve
    if(!solver.data()->setHessianMatrix(H)){
        cout << "H matrix error!" << endl;
        return 1;
    } 
    if(!solver.data()->setGradient(f)){
        cout << "f matrix error!" << endl;
        return 1;
    }
    if(!solver.data()->setLinearConstraintsMatrix(A)){
        cout << "A matrix error!" << endl;
        return 1;
    }
    if(!solver.data()->setLowerBound(lb)){
        cout << "lb matrix error!" << endl;
        return 1;
    }
    if(!solver.data()->setUpperBound(ub)){
        cout << "ub matrix error!" << endl;
        return 1;
    }

    // initial the solver
    if(!solver.initSolver()){
        cout << "solver initial error!" << endl;
        return 1;
    }

    // solve the problem
    auto err_code = solver.solveProblem();
    if(err_code != OsqpEigen::ErrorExitFlag::NoError){
        cout << "solve err: " << static_cast<int>(err_code) << endl;
        return 1;
    }

    Eigen::Vector2d solution = solver.getSolution();
    double obj = solver.getObjValue();

    cout << "the solution is: " << endl;
    cout << solution << endl;
    cout << "the obj value is: " << obj << endl;

    ///////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////
    // update H and A
    Eigen::SparseMatrix<double> H_update(2,2);
    Eigen::SparseMatrix<double> A_update(3,2);

    H_update.insert(0,0) = 5.0;
    H_update.insert(0,1) = 1.5; 
    H_update.insert(1,0) = 1.5;
    H_update.insert(1,1) = 1.0;

    A_update.insert(0,0) = 1.2;
    A_update.insert(0,1) = 1.1;
    A_update.insert(1,0) = 1.5;
    // A_update.insert(1,1) = 0.0;
    // A_update.insert(2,0) = 0.0;
    A_update.insert(2,1) = 0.8;

    cout << A_update << endl;

    solver.updateHessianMatrix(H_update);
    solver.updateLinearConstraintsMatrix(A_update);

    err_code = solver.solveProblem();
    if(err_code != OsqpEigen::ErrorExitFlag::NoError){
        cout << "solve err: " << static_cast<int>(err_code) << endl;
        return 1;
    }

    solution = solver.getSolution();
    obj = solver.getObjValue();

    cout << "the solution after updating is: " << endl;
    cout << solution << endl;
    cout << "the obj value after updating is: " << obj << endl;
    return 0;
}

The different between these two version of code is the Constrains Matrix construction.

In the wrong answer version, I delete two lines of code, which has no affect to the matrix.

A_update.insert(1,1) = 0.0;
A_update.insert(2,0) = 0.0;

[little-bird-go]

Can you post the issue directly as text, instead of posting it as an image and then a markdown attached to the issue? Thanks!

Furthermore, does this issue also happens if you specifying the problem for example in Python? In that case, the problem may not be in the osqp-eigen package, but rather in the osqp library itself: https://github.com/osqp .

I have tried to use osqp C interface to test the problem, and got the right solution. I wonder if there is a bug in updateLinearConstraintsMatrix function.

[traversaro]

Thanks, this is useful to isolate the problem. Indeed this seems something is going wrong in the updateLinearConstraintsMatrix method. Perhaps it was never tested with sparse matrix that have zero values in their elements? Can you provide the ucode that you used to test the osqp C interface? In particular, which triplets did you passed?

[traversaro]

It would be also interesting to debug a bit in https://github.com/robotology/osqp-eigen/blob/a752c5e6353cae6176578488a1fa641d3b9cf0b9/include/OsqpEigen/Solver.tpp#L163-L296 to understand if the code thinks that the sparsity pattern changed or not (it is a bit of a corner case as we are inserting 0.0 as explicit element in a sparse matrix).

[traversaro]

By the way, by looking in the documentation of insert's Eigen method in https://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html#ae2d8f72ff86a300b76f9edd67df8d8fd, I am not sure if the usecase insert([..]) = 0.0 is actually supported by Eigen itself.

[little-bird-go]

#include <iostream>
#include <osqp/osqp.h>
#include <stdlib.h>
#include <stdio.h>

using namespace std;

int main(void)
{
    // Hessian Matrix
    OSQPFloat P_x[3] = {4.0, 1.0, 2.0};
    OSQPInt P_nnz = 3;
    OSQPInt P_i[3] = {0, 0, 1};
    OSQPInt P_p[3] = {0, 1, 3};

    // Gradient Vector
    OSQPFloat q[2] = {1.0, 1.0};

    // Linear Constraints Matrix
    OSQPFloat A_x[4] = {1.0, 1.0, 1.0, 1.0};
    OSQPInt A_nnz = 4;
    OSQPInt A_i[4] = {0, 1, 0, 2};
    OSQPInt A_p[3] = {0, 2, 4};

    // Lower Bounds
    OSQPFloat lb[3] = {1.0, 0.0, 0.0};
    OSQPFloat ub[3] = {1.0, 0.7, 0.7};

    // Problem dimention
    OSQPInt n = 2;
    OSQPInt m = 3;

        /* Exitflag */
    OSQPInt exitflag;

    /* Solver, settings, matrices */
    OSQPSolver*   solver   = NULL;
    OSQPSettings* settings = NULL;
    OSQPCscMatrix* P = static_cast<OSQPCscMatrix*> (malloc(sizeof(OSQPCscMatrix)));
    OSQPCscMatrix* A = static_cast<OSQPCscMatrix*> (malloc(sizeof(OSQPCscMatrix)));

    /* Populate matrices */
    csc_set_data(A, m, n, A_nnz, A_x, A_i, A_p);
    csc_set_data(P, n, n, P_nnz, P_x, P_i, P_p);

    /* Set default settings */
    settings = (OSQPSettings *)malloc(sizeof(OSQPSettings));
    if (settings) {
        osqp_set_default_settings(settings);
        //settings->polishing = 1;
        //settings->linsys_solver = OSQP_DIRECT_SOLVER;
        //settings->linsys_solver = OSQP_INDIRECT_SOLVER;
    }

    OSQPInt cap = osqp_capabilities();

    printf("This OSQP library supports:\n");
    if(cap & OSQP_CAPABILITY_DIRECT_SOLVER) {
        printf("    A direct linear algebra solver\n");
    }
    if(cap & OSQP_CAPABILITY_INDIRECT_SOLVER) {
        printf("    An indirect linear algebra solver\n");
    }
    if(cap & OSQP_CAPABILITY_CODEGEN) {
        printf("    Code generation\n");
    }
    if(cap & OSQP_CAPABILITY_DERIVATIVES) {
        printf("    Derivatives calculation\n");
    }
    printf("\n");

    /* Setup solver */
    exitflag = osqp_setup(&solver, P, q, A, lb, ub, m, n, settings);

    /* Solve problem */
    if (!exitflag) exitflag = osqp_solve(solver);

    cout << solver->solution->x[0] << " " << solver->solution->x[1] << endl;

    OSQPFloat P_xnew[3] = {5.0, 1.5, 1.0};
    OSQPInt P_new_n = 3;
    OSQPInt P_new_i[3] = {0, 1, 2};

    OSQPFloat A_xnew[4] = {1.2, 1.5, 1.1, 0.8};
    OSQPInt A_new_n = 4;
    OSQPInt A_new_i[4] = {0, 1, 2, 3};

    exitflag = osqp_update_data_mat(solver, P_xnew, P_new_i, P_new_n, A_xnew, A_new_i, A_new_n);
    if (!exitflag) exitflag = osqp_solve(solver);

    cout << solver->solution->x[0] << " " << solver->solution->x[1] << endl;

    return 0;
}

Above is the code I test with osqp. In this case or let A_new_idx and P_new_idx be OSQP_NULL are both OK.

[traversaro]

In this case or let A_new_idx and P_new_idx be OSQP_NULL are both OK.

Sorry, I am not sure what you mean with this sentence.

[little-bird-go]

I also think the problem is caused by inserting zeros in the sparematrix with eigen, which may create different matrix in the memory. But I don't know how eigen works with sparematrix. When I insert zeros to the matrix, I get When I do not insert zeros, I get

[little-bird-go]

I mean that using osqp_update_data_mat(solver, P_xnew, P_new_i, P_new_n, A_xnew, A_new_i, A_new_n); or osqp_update_data_mat(solver, P_xnew, OSQP_NULL, P_new_n, A_xnew,OSQP_NULL, A_new_n); are both OK.