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

osqp-eigen/include/OsqpEigen/Solver.tpp

Lines 163 to 296 in a752c5e

	if (!m_isSolverInitialized)
	{
	debugStream() << "[OsqpEigen::Solver::updateLinearConstraintsMatrix] The solver has not "
	"been initialized."
	<< std::endl;
	return false;
	}
	if (((c_int)linearConstraintsMatrix.rows() != getData()->m)
	|| ((c_int)linearConstraintsMatrix.cols() != getData()->n))
	{
	debugStream() << "[OsqpEigen::Solver::updateLinearConstraintsMatrix] The constraints "
	"matrix has to be a mxn matrix"
	<< std::endl;
	return false;
	}
	// evaluate the triplets from old and new hessian sparse matrices
	if (!OsqpEigen::SparseMatrixHelper::osqpSparseMatrixToTriplets(getData()->A,
	m_oldLinearConstraintsTriplet))
	{
	debugStream() << "[OsqpEigen::Solver::updateLinearConstraintsMatrix] Unable to evaluate "
	"triplets from the old hessian matrix."
	<< std::endl;
	return false;
	}
	if (!OsqpEigen::SparseMatrixHelper::eigenSparseMatrixToTriplets(linearConstraintsMatrix,
	m_newLinearConstraintsTriplet))
	{
	debugStream() << "[OsqpEigen::Solver::updateLinearConstraintsMatrix] Unable to evaluate "
	"triplets from the old hessian matrix."
	<< std::endl;
	return false;
	}
	// try to update the linear constraints matrix without reinitialize the solver
	// according to the osqp library it can be done only if the sparsity pattern of the
	// matrix does not change.
	if (evaluateNewValues(m_oldLinearConstraintsTriplet,
	m_newLinearConstraintsTriplet,
	m_constraintsNewIndices,
	m_constraintsNewValues))
	{
	if (m_constraintsNewValues.size() > 0)
	{
	#ifdef OSQP_EIGEN_OSQP_IS_V1
	if (osqp_update_data_mat(m_solver.get(),
	nullptr,
	nullptr,
	0,
	m_constraintsNewValues.data(),
	m_constraintsNewIndices.data(),
	m_constraintsNewIndices.size())
	!= 0)
	{
	#else
	if (osqp_update_A(m_workspace.get(),
	m_constraintsNewValues.data(),
	m_constraintsNewIndices.data(),
	m_constraintsNewIndices.size())
	!= 0)
	{
	#endif
	debugStream() << "[OsqpEigen::Solver::updateLinearConstraintsMatrix] Unable to "
	"update linear constraints matrix."
	<< std::endl;
	return false;
	}
	}
	} else
	{
	// the sparsity pattern has changed
	// the solver has to be setup again
	// get the primal and the dual variables
	if (!getPrimalVariable(m_primalVariables))
	{
	debugStream() << "[OsqpEigen::Solver::updateLinearConstraintsMatrix] Unable to get the "
	"primal variable."
	<< std::endl;
	return false;
	}
	if (!getDualVariable(m_dualVariables))
	{
	debugStream() << "[OsqpEigen::Solver::updateLinearConstraintsMatrix] Unable to get the "
	"dual variable."
	<< std::endl;
	return false;
	}
	// clear old linear constraints matrix
	m_data->clearLinearConstraintsMatrix();
	// set new linear constraints matrix
	if (!m_data->setLinearConstraintsMatrix(linearConstraintsMatrix))
	{
	debugStream() << "[OsqpEigen::Solver::updateLinearConstraintsMatrix] Unable to update "
	"the hessian matrix in "
	<< "Data object." << std::endl;
	return false;
	}
	// clear the old solver
	clearSolver();
	if (!initSolver())
	{
	debugStream() << "[OsqpEigen::Solver::updateLinearConstraintsMatrix] Unable to "
	"Initialize the solver."
	<< std::endl;
	return false;
	}
	// set the old primal and dual variables
	if (!setPrimalVariable(m_primalVariables))
	{
	debugStream() << "[OsqpEigen::Solver::updateLinearConstraintsMatrix] Unable to set the "
	"primal variable."
	<< std::endl;
	return false;
	}
	if (!setDualVariable(m_dualVariables))
	{
	debugStream() << "[OsqpEigen::Solver::updateLinearConstraintsMatrix] Unable to set the "
	"dual variable."
	<< std::endl;
	return false;
	}
	}

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.