How to use quadprog - 10 common examples

# Set options of the OSQP solver
options = {'eps_abs':       1e-4,
           'eps_rel':       1e-4,
           'alpha':         1.6,
           'scale_problem': True,
           'scale_steps':   4,
           'polish':        False}

# Create an svm object
basis_pursuit_obj = basis_pursuit(m, n, dens_lvl=0.3, osqp_opts=options)

# Solve with different solvers
resultsCPLEX = basis_pursuit_obj.solve(solver=CPLEX)
resultsGUROBI = basis_pursuit_obj.solve(solver=GUROBI)
resultsOSQP = basis_pursuit_obj.solve(solver=OSQP)


# Print objective values
print "CPLEX  Objective Value: %.3f" % resultsCPLEX.objval
print "GUROBI Objective Value: %.3f" % resultsGUROBI.objval
print "OSQP   Objective Value: %.3f" % resultsOSQP.objval
print "\n"

# Print timings
print "CPLEX  CPU time: %.3f" % resultsCPLEX.cputime
print "GUROBI CPU time: %.3f" % resultsGUROBI.cputime
print "OSQP   CPU time: %.3f" % resultsOSQP.cputime

# Set options of the OSQP solver
options = {'eps_abs':       1e-4,
           'eps_rel':       1e-4,
           'alpha':         1.6,
           'scale_problem': True,
           'scale_steps':   4,
           'polish':        False}

# Create a lasso object
huber_fit_obj = huber_fit(m, n, dens_lvl=0.50, osqp_opts=options)

# Solve with different solvers
resultsCPLEX = huber_fit_obj.solve(solver=CPLEX)
resultsGUROBI = huber_fit_obj.solve(solver=GUROBI)
resultsOSQP = huber_fit_obj.solve(solver=OSQP)

# Print objective values
print "CPLEX  Objective Value: %.3f" % resultsCPLEX.objval
print "GUROBI Objective Value: %.3f" % resultsGUROBI.objval
print "OSQP   Objective Value: %.3f" % resultsOSQP.objval
print "\n"

# Print timings
print "CPLEX  CPU time: %.3f" % resultsCPLEX.cputime
print "GUROBI CPU time: %.3f" % resultsGUROBI.cputime
print "OSQP   CPU time: %.3f" % resultsOSQP.cputime

# ipdb.set_trace()

# Recover A, x and b from the problem
A = huber_fit_obj._osqp.problem.A[:m, :n]

m = 200

# Set options of the OSQP solver
options = {'eps_abs':       1e-4,
           'eps_rel':       1e-4,
           'alpha':         1.6,
           'scale_problem': True,
           'scale_steps':   4,
           'polish':        False}

# Create an svm object
basis_pursuit_obj = basis_pursuit(m, n, dens_lvl=0.3, osqp_opts=options)

# Solve with different solvers
resultsCPLEX = basis_pursuit_obj.solve(solver=CPLEX)
resultsGUROBI = basis_pursuit_obj.solve(solver=GUROBI)
resultsOSQP = basis_pursuit_obj.solve(solver=OSQP)


# Print objective values
print "CPLEX  Objective Value: %.3f" % resultsCPLEX.objval
print "GUROBI Objective Value: %.3f" % resultsGUROBI.objval
print "OSQP   Objective Value: %.3f" % resultsOSQP.objval
print "\n"

# Print timings
print "CPLEX  CPU time: %.3f" % resultsCPLEX.cputime
print "GUROBI CPU time: %.3f" % resultsGUROBI.cputime
print "OSQP   CPU time: %.3f" % resultsOSQP.cputime

m = 10*n

# Set options of the OSQP solver
options = {'eps_abs':       1e-4,
           'eps_rel':       1e-4,
           'alpha':         1.6,
           'scale_problem': True,
           'scale_steps':   4,
           'polish':        False}

# Create a lasso object
huber_fit_obj = huber_fit(m, n, dens_lvl=0.50, osqp_opts=options)

# Solve with different solvers
resultsCPLEX = huber_fit_obj.solve(solver=CPLEX)
resultsGUROBI = huber_fit_obj.solve(solver=GUROBI)
resultsOSQP = huber_fit_obj.solve(solver=OSQP)

# Print objective values
print "CPLEX  Objective Value: %.3f" % resultsCPLEX.objval
print "GUROBI Objective Value: %.3f" % resultsGUROBI.objval
print "OSQP   Objective Value: %.3f" % resultsOSQP.objval
print "\n"

# Print timings
print "CPLEX  CPU time: %.3f" % resultsCPLEX.cputime
print "GUROBI CPU time: %.3f" % resultsGUROBI.cputime
print "OSQP   CPU time: %.3f" % resultsOSQP.cputime

# ipdb.set_trace()

# Recover A, x and b from the problem

P = spspa.block_diag((2*D, 2*spspa.eye(k)), format='csc')
            q = np.append(-mu / gamma, np.zeros(k))
            A = spspa.vstack([
                    spspa.hstack([spspa.csc_matrix(np.ones((1, n))),
                                  spspa.csc_matrix((1, k))]),
                    spspa.hstack([F.T, -spspa.eye(k)]),
                    spspa.hstack([spspa.eye(n), spspa.csc_matrix((n, k))])
                ]).tocsc()
            lA = np.hstack([1., np.zeros(k), np.zeros(n)])
            uA = np.hstack([1., np.zeros(k), np.ones(n)])

        else:
            assert False, "Unhandled version"

        # Create a quadprogProblem and store it in a private variable
        self._prob = qp.quadprogProblem(P, q, A, lA, uA)
        # Create an OSQP object and store it in a private variable
        self._osqp = osqp.OSQP(**osqp_opts)
        self._osqp.problem(P, q, A, lA, uA)

uA = np.array([np.inf, np.inf])
        p = qp.quadprogProblem(P, q, A, lA, uA)
    elif example == 'random':
        # Random Example
        n = 30
        m = 50
        # Generate random Matrices
        Pt = sp.randn(n, n)
        P = spspa.csc_matrix(np.dot(Pt.T, Pt))
        q = sp.randn(n)
        A = spspa.csc_matrix(sp.randn(m, n))
        uA = 3 + sp.randn(m)
        # lA = uA
        lA = -3 + sp.randn(m)

        p = qp.quadprogProblem(P, q, A, lA, uA)
    elif example == 'lp':
        # Random Example
        n = 100
        m = 50
        # Generate random Matrices
        P = spspa.csc_matrix(np.zeros((n, n)))
        q = sp.randn(n)
        A = spspa.vstack([spspa.csc_matrix(sp.randn(m, n)), spspa.eye(n)])
        lA = np.append(- 3 + sp.randn(m), - 3 + sp.randn(n))
        uA = np.append(3 + sp.randn(m), 3 + sp.randn(n))
        p = qp.quadprogProblem(P, q, A, lA, uA)
    else:
        assert False, "Unknown example"

    # Solve with CPLEX
    resultsCPLEX = p.solve(solver=CPLEX, verbose=1)

if example == 'random':
        # Random Example
        nx = 100
        neq = 10
        nineq = 20
        # Generate random Matrices
        Qt = sp.randn(nx, nx)
        Q = spspa.csc_matrix(np.dot(Qt.T, Qt))
        c = sp.randn(nx)
        Aeq = spspa.csc_matrix(sp.randn(neq, nx))
        beq = sp.randn(neq)
        Aineq = spspa.csc_matrix(sp.randn(nineq, nx))
        bineq = 100 * sp.rand(nineq)
        lb = 0. * np.ones(nx)
        ub = 5. * np.ones(nx)
        p = qp.quadprogProblem(Q, c, Aeq, beq, Aineq, bineq, lb, ub)
    else:
        assert False, "Unknown example"

    for i in range(2):
        if example == 'random' and i >= 1:
            # Reuse factorizations
            c = sp.randn(nx)
            beq = sp.randn(neq)
            bineq = 100 * sp.rand(nineq)
            p = qp.quadprogProblem(Q, c, Aeq, beq, Aineq, bineq, lb, ub)

        # Solve with GUROBI
        resultsGUROBI = p.solve(solver=GUROBI, OutputFlag=0)

        # Solve with OSQP
        if i == 0:

#                   v >= 0
        Im = spspa.eye(m)
        P = spspa.block_diag((spspa.csc_matrix((n, n)), Im,
                              spspa.csc_matrix((m, m))), format='csc')
        q = np.append(np.zeros(m + n), np.ones(m))
        A = spspa.vstack([
                spspa.hstack([A, Im, Im]),
                spspa.hstack([A, -Im, -Im]),
                spspa.hstack([spspa.csc_matrix((m, n)), Im,
                              spspa.csc_matrix((m, m))]),
                spspa.hstack([spspa.csc_matrix((m, n + m)), Im])]).tocsc()
        lA = np.hstack([b, -np.inf*np.ones(m), np.zeros(2*m)])
        uA = np.hstack([np.inf*np.ones(m), b, np.ones(m), np.inf*np.ones(m)])

        # Create a quadprogProblem and store it in a private variable
        self._prob = qp.quadprogProblem(P, q, A, lA, uA)
        # Create an OSQP object and store it in a private variable
        self._osqp = osqp.OSQP(**osqp_opts)
        self._osqp.problem(P, q, A, lA, uA)

# Infeasible example
        # P = spspa.eye(2)
        P = spspa.csc_matrix((2, 2))
        q = np.ones(2)
        A = spspa.csc_matrix(np.array([[1, 0], [0, 1], [1, 1]]))
        lA = np.array([0., 0., -1.])
        uA = np.array([np.inf, np.inf, -1.])
        p = qp.quadprogProblem(P, q, A, lA, uA)
    elif example == 'unbounded':
        # Unbounded example
        P = spspa.csc_matrix((2, 2))
        q = np.array([2, -1])
        A = spspa.eye(2)
        lA = np.array([0., 0.])
        uA = np.array([np.inf, np.inf])
        p = qp.quadprogProblem(P, q, A, lA, uA)
    elif example == 'random':
        # Random Example
        n = 30
        m = 50
        # Generate random Matrices
        Pt = sp.randn(n, n)
        P = spspa.csc_matrix(np.dot(Pt.T, Pt))
        q = sp.randn(n)
        A = spspa.csc_matrix(sp.randn(m, n))
        uA = 3 + sp.randn(m)
        # lA = uA
        lA = -3 + sp.randn(m)

        p = qp.quadprogProblem(P, q, A, lA, uA)
    elif example == 'lp':
        # Random Example

# Load file
    m = spio.loadmat(f)

    # Convert matrices
    Q = m['Q'].astype(float)
    c = m['c'].T.flatten().astype(float)
    Aeq = m['A'].astype(float)
    beq = m['ru'].T.flatten().astype(float)
    lb = m['lb'].T.flatten().astype(float)
    ub = m['ub'].T.flatten().astype(float)
    nx = Q.shape[0]
    Aineq = spspa.csc_matrix(np.zeros((1, nx)))
    bineq = np.array([0.0])

    # Define problem
    p = qp.quadprogProblem(Q, c, Aeq, beq, Aineq, bineq, lb, ub)

    return p