cuda OSQP examples

[Eqedos]

I was looking for resources to find examples where the CUDA back-end beats the CPU, unfortunately couldn't find any. I was wondering if anyone could provide examples/code for such problems.

[vineetbansal]

@Eqedos - the repo osqp-benchmarks has a lot of sample problems. In particular, the Huber class of problems benefit from the CUDA backend, for higher dimensions of problems. See a benchmark graphic below (blue is built-in backend, orange is cuda backend):

A sample standalone script representing the Huber problem would be:

import numpy as np
import scipy.sparse as spa
from osqp import OSQP


class HuberExample:

    def __init__(self, n, seed=1):
        np.random.seed(seed)

        self.n = int(n)  # Number of features
        self.m = int(self.n * 100)  # Number of data-points

        self.Ad = spa.random(self.m, self.n, density=0.15, data_rvs=np.random.randn)
        self.x_true = np.random.randn(n) / np.sqrt(n)
        ind95 = (np.random.rand(self.m) < 0.95).astype(float)
        self.bd = (
            self.Ad.dot(self.x_true)
            + np.multiply(0.5 * np.random.randn(self.m), ind95)
            + np.multiply(10.0 * np.random.rand(self.m), 1.0 - ind95)
        )

    def get_problem(self):
        # Construct the problem
        #       minimize    1/2 z.T * z + np.ones(m).T * (r + s)
        #       subject to  Ax - b - z = r - s
        #                   r >= 0
        #                   s >= 0
        # The problem reformulation follows from Eq. (24) of the following paper:
        # https://doi.org/10.1109/34.877518
        # x_solver = (x, z, r, s)
        Im = spa.eye(self.m)
        P = spa.block_diag(
            (
                spa.csc_matrix((self.n, self.n)),
                Im,
                spa.csc_matrix((2 * self.m, 2 * self.m)),
            ),
            format="csc",
        )
        q = np.hstack([np.zeros(self.n + self.m), np.ones(2 * self.m)])
        A = spa.bmat(
            [[self.Ad, -Im, -Im, Im], [None, None, Im, None], [None, None, None, Im]],
            format="csc",
        )
        l = np.hstack([self.bd, np.zeros(2 * self.m)])
        u = np.hstack([self.bd, np.inf * np.ones(2 * self.m)])

        # Constraints without bounds
        A_nobounds = spa.hstack([self.Ad, -Im, -Im, Im], format="csc")
        l_nobounds = self.bd
        u_nobounds = self.bd

        # Bounds
        lx = np.hstack([-np.inf * np.ones(self.n + self.m), np.zeros(2 * self.m)])
        ux = np.inf * np.ones(self.n + 3 * self.m)
        bounds_idx = np.arange(self.n + self.m, self.n + 3 * self.m)

        problem = {}
        problem["P"] = P
        problem["q"] = q
        problem["A"] = A
        problem["l"] = l
        problem["u"] = u
        problem["m"] = A.shape[0]
        problem["n"] = A.shape[1]
        problem["A_nobounds"] = A_nobounds
        problem["l_nobounds"] = l_nobounds
        problem["u_nobounds"] = u_nobounds
        problem["bounds_idx"] = bounds_idx
        problem["lx"] = lx
        problem["ux"] = ux

        return problem


if __name__ == "__main__":

    settings = {
        "eps_abs": 0.001,
        "eps_rel": 0.001,
        "polish": False,
        "max_iter": 1000000000,
        "eps_prim_inf": 1e-15,
        "eps_dual_inf": 1e-15,
        "check_dualgap": True,
        "restart_enable": True,
        "restart_sufficient": 0.3875,
        "restart_necessary": 0.825,
        "restart_artificial": 0.5,
        "verbose": False,
        "time_limit": 1000.0,
        "adaptive_rho_iterval": 5,
        "adaptive_rho_tolerance": 5,
        "verbose": True
    }

    problem = HuberExample(200).get_problem()

    solvers = [OSQP(algebra="builtin"), OSQP(algebra="cuda")]
    for solver in solvers:
        solver.setup(problem['P'], problem['q'], problem['A'], problem['l'],
                problem['u'],
                **settings)

        print(f"Solving using {solver}")
        solver.solve(problem)

Running on my local machine gives me:

-----------------------------------------------------------------
           OSQP v1.0.0.beta0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2021
-----------------------------------------------------------------
problem:  variables n = 60200, constraints m = 60000
          nnz(P) + nnz(A) = 720000
settings: algebra = Built-in,
          linear system solver = QDLDL v0.1.6,
          eps_abs = 1.0e-03, eps_rel = 1.0e-03,
          eps_prim_inf = 1.0e-15, eps_dual_inf = 1.0e-15,
          rho = 1.00e-01 (adaptive: 50 iterations),
          restart: on (sufficient: 0.39, necessary: 0.82, artificial: 0.36),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 1000000000
          check_termination: on (interval 25, duality gap: on),
          time_limit: 1.00e+03 sec,
          scaling: on, scaled_termination: off
          warm starting: on, polishing: off, 
Solving using OSQP with algebra=builtin (direct)
iter    objective    prim res   dual res   gap        rho         time
   1   -6.3999e+05   1.60e+01   4.36e+04  -3.36e+06   1.00e-01    2.56e+00s
  34*  -2.7620e+03   3.53e-01   2.33e+02  -2.37e+04   1.86e-02*   2.69e+00s
  38*   4.5432e+03   1.60e-01   4.99e+00  -2.17e+03   2.31e-03*   2.73e+00s
  50    5.5906e+03   1.35e-01   6.44e-01  -8.57e+02   4.49e-04*   2.81e+00s
  82*   5.7210e+03   1.28e-01   2.07e-05  -6.88e+02   1.10e-02*   2.96e+00s
 200    6.2265e+03   4.60e-02   6.91e-04  -2.54e+02   1.10e-02    3.45e+00s
 400    6.4441e+03   6.75e-03   5.85e-05  -4.52e+01   1.10e-02    4.22e+00s
 450    6.4604e+03   4.34e-03   3.92e-06  -2.91e+01   1.14e-01*   4.42e+00s
 475    6.4858e+03   6.65e-04   8.37e-05  -3.67e+00   1.14e-01    4.56e+00s

status:               solved
number of iterations: 475
optimal objective:    6485.7954
dual objective:       6489.4610
duality gap:          -3.6656e+00
primal-dual integral: 1.6106e+07
run time:             4.56e+00s
optimal rho estimate: 9.96e-02

-----------------------------------------------------------------
           OSQP v1.0.0.beta0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2021
-----------------------------------------------------------------
problem:  variables n = 60200, constraints m = 60000
          nnz(P) + nnz(A) = 720000
settings: algebra = CUDA 11.2,
          device = NVIDIA GeForce RTX 3050 Laptop GPU (Compute capability 8.6),
          linear system solver = CUDA Conjugate Gradient - Diagonal preconditioner,
          eps_abs = 1.0e-03, eps_rel = 1.0e-03,
          eps_prim_inf = 1.0e-15, eps_dual_inf = 1.0e-15,
          rho = 1.00e-01 (adaptive: 50 iterations),
          restart: on (sufficient: 0.39, necessary: 0.82, artificial: 0.36),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 1000000000
          check_termination: on (interval 5, duality gap: on),
          time_limit: 1.00e+03 sec,
          scaling: on, scaled_termination: off
          warm starting: on, polishing: off, 
Solving using OSQP with algebra=cuda (indirect)
iter    objective    prim res   dual res   gap        rho         time
   1   -6.3982e+05   1.63e+01   4.36e+01  -6.43e+05   1.00e-01    1.68e+00s
 170    6.4846e+03   1.06e-02   3.64e-04  -4.98e+00   1.00e-01    2.06e+00s

status:               solved
number of iterations: 170
optimal objective:    6484.6396
dual objective:       6489.6182
duality gap:          -4.9785e+00
primal-dual integral: 3.0667e+06
run time:             2.06e+00s
optimal rho estimate: 1.16e-01

[Eqedos]

Thank you! very helpful. Are there perhaps classes of problems the cuda OSQP backend can't solve at all or does not converge on? I've tried solving a reasonably small problem, and cuda osqp does not converge via CVXPY. I was wondering if CVXPYs canonicalization had anything to do with it, or perhaps it's the type of problem.

                                     CVXPY                                     
                                     v1.5.3                                    
===============================================================================
(CVXPY) Oct 06 08:19:20 PM: Your problem has 7326 variables, 2496 constraints, and 1656 parameters.
(CVXPY) Oct 06 08:19:21 PM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Oct 06 08:19:21 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
(CVXPY) Oct 06 08:19:21 PM: Your problem is compiled with the CPP canonicalization backend.
-------------------------------------------------------------------------------
                                  Compilation                                  
-------------------------------------------------------------------------------
(CVXPY) Oct 06 08:19:21 PM: Compiling problem (target solver=OSQP).
(CVXPY) Oct 06 08:19:21 PM: Reduction chain: EvalParams -> CvxAttr2Constr -> Qp2SymbolicQp -> QpMatrixStuffing -> OSQP
(CVXPY) Oct 06 08:19:21 PM: Applying reduction EvalParams
(CVXPY) Oct 06 08:19:21 PM: Applying reduction CvxAttr2Constr
(CVXPY) Oct 06 08:19:22 PM: Applying reduction Qp2SymbolicQp
(CVXPY) Oct 06 08:19:22 PM: Applying reduction QpMatrixStuffing
(CVXPY) Oct 06 08:19:24 PM: Applying reduction OSQP
(CVXPY) Oct 06 08:19:24 PM: Finished problem compilation (took 3.346e+00 seconds).
-------------------------------------------------------------------------------
                                Numerical solver                               
-------------------------------------------------------------------------------
(CVXPY) Oct 06 08:19:24 PM: Invoking solver OSQP  to obtain a solution.
-----------------------------------------------------------------
           OSQP v1.0.0.beta0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2021
-----------------------------------------------------------------
problem:  variables n = 7380, constraints m = 17178
          nnz(P) + nnz(A) = 38941
settings: algebra = CUDA 11.5,
          device = NVIDIA GeForce RTX 3060 Ti (Compute capability 8.6),
          linear system solver = CUDA Conjugate Gradient - Diagonal preconditioner,
          eps_abs = 1.0e-03, eps_rel = 1.0e-03,
          eps_prim_inf = 1.0e-15, eps_dual_inf = 1.0e-15,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 1000000000
          check_termination: on (interval 5),
          time_limit: 1.00e+03 sec,
          scaling: on, scaled_termination: off
          warm starting: on, polishing: off, ```

[imciner2]

Are you using a double or float version of the library (or just installing from pip - I think that is double)? We have noticed that the float version tends to perform worse on the CUDA backend currently, probably because of the conjugate gradient method we are using to solve the KKT system.

If you are able to share the problem, it would be a good test case for us to examine what is happening. We are actually working on some algorithmic improvements (restarting and different adaptive rho methods) that might also help make these converge.

[Eqedos]

Is there a way to check if I'm using float/double? I believe I just pip installed using the build instructions in the wiki, using wheels.