Merge pull request #61 from cvxgrp/add_clarabel

add clarabel

.gitignore

@@ -515,4 +515,6 @@ MigrationBackup/
 # Ionide (cross platform F# VS Code tools) working folder
 .ionide/
 
-.vscode
+.vscode
+.direnv/
+.envrc

README.md

@@ -22,6 +22,7 @@ You will need a C++11-capable compiler to build `diffcp`.
 * [pybind11](https://github.com/pybind/pybind11/tree/stable) >= 2.4
 * [threadpoolctl](https://github.com/joblib/threadpoolctl) >= 1.1
 * [ECOS](https://github.com/embotech/ecos-python) >= 2.0.10
+* [Clarabel](https://github.com/oxfordcontrol/Clarabel.rs) >= 0.5.1
 * Python >= 3.7
 
 `diffcp` uses Eigen; Eigen operations can be automatically vectorized by compilers. To enable vectorization, install with
@@ -73,7 +74,7 @@ functions for evaluating the derivative and its adjoint (transpose).
 These functions respectively compute right and left multiplication of the derivative
 of the solution map at `A`, `b`, and `c` by a vector.
 The `solver` argument determines which solver to use; the available solvers
-are `solver="SCS"` and `solver="ECOS"`.
+are `solver="SCS"`, `solver="ECOS"`, and `solver="Clarabel"`.
 If no solver is specified, `diffcp` will choose the solver itself.
 In the case that the problem is not solved, i.e. the solver fails for some reason, we will raise
 a `SolverError` Exception.

diffcp/cone_program.py

@@ -3,6 +3,7 @@
 from multiprocessing.pool import ThreadPool
 
 import ecos
+import clarabel
 import numpy as np
 import scipy.sparse as sparse
 import scs
@@ -192,7 +193,7 @@ def solve_and_derivative(A, b, c, cone_dict, warm_start=None, mode='lsqr',
       warm_start: (optional) A tuple (x, y, s) at which to warm-start SCS.
       mode: (optional) Which mode to compute derivative with, options are
           ["dense", "lsqr", "lsmr"].
-      solve_method: (optional) Name of solver to use; SCS or ECOS.
+      solve_method: (optional) Name of solver to use; SCS, ECOS, or Clarabel.
       kwargs: (optional) Keyword arguments to send to the solver.
 
     Returns:
@@ -248,6 +249,8 @@ def solve_and_derivative_internal(A, b, c, cone_dict, solve_method=None,
     if solve_method is None:
         psd_cone = ('s' in cone_dict) and (cone_dict['s'] != [])
         ed_cone = ('ed' in cone_dict) and (cone_dict['ed'] != 0)
+
+        # TODO(sbarratt): consider setting default to clarabel
         if psd_cone or ed_cone:
             solve_method = "SCS"
         else:
@@ -386,7 +389,65 @@ def solve_and_derivative_internal(A, b, c, cone_dict, solve_method=None,
                           'setupTime': solution['info']['timing']['tsetup'],
                           'iter': solution['info']['iter'],
                           'pobj': solution['info']['pcost']}
+    elif solve_method == "Clarabel":
+        # for now set P to 0
+        P = sparse.csc_matrix((c.size, c.size))
+
+        cones = []
+        if "z" in cone_dict:
+            v = cone_dict["z"]
+            if v > 0:
+                cones.append(clarabel.ZeroConeT(v))
+        if "f" in cone_dict:
+            v = cone_dict["f"]
+            if v > 0:
+                cones.append(clarabel.ZeroConeT(v))
+        if "l" in cone_dict:
+            v = cone_dict["l"]
+            if v > 0:
+                cones.append(clarabel.NonnegativeConeT(v))
+        if "q" in cone_dict:
+            for v in cone_dict["q"]:
+                cones.append(clarabel.SecondOrderConeT(v))
+        if "s" in cone_dict:
+            for v in cone_dict["s"]:
+                cones.append(clarabel.PSDTriangleConeT(v))
+        if "ep" in cone_dict:
+            v = cone_dict["ep"]
+            cones += [clarabel.ExponentialConeT()] * v
+
+        kwargs.setdefault("verbose", False)
+        settings = clarabel.DefaultSettings()
+
+        for key, value in kwargs.items():
+            setattr(settings, key, value)
+
+        solver = clarabel.DefaultSolver(P,c,A,b,cones,settings)
+        solution = solver.solve()
 
+        result = {}
+        result["x"] = np.array(solution.x)
+        result["y"] = np.array(solution.z)
+        result["s"] = np.array(solution.s)
+
+        x, y, s = result["x"], result["y"], result["s"]
+
+        CLARABEL2SCS_STATUS_MAP = {
+            "Solved": "Solved",
+            "PrimalInfeasible": "Infeasible",
+            "DualInfeasible": "Unbounded",
+            "AlmostSolved": "Optimal Inaccurate",
+            "AlmostPrimalInfeasible": "Infeasible Inaccurate",
+            "AlmostDualInfeasible": "Unbounded Inaccurate",
+        }
+
+        result["info"] = {
+            "status": CLARABEL2SCS_STATUS_MAP.get(str(solution.status), "Failure"),
+            "solveTime": solution.solve_time,
+            "setupTime": -1,
+            "iter": solution.iterations,
+            "pobj": solution.obj_val,
+        }
     else:
         raise ValueError("Solver %s not supported." % solve_method)
 

examples/hello_world.py

@@ -6,7 +6,7 @@
 
 
 cone_dict = {
-    'f': 3,
+    'z': 3,
     'l': 3,
     'q': [5]
 }
@@ -19,7 +19,7 @@
 A, b, c = utils.random_cone_prog(m, n, cone_dict)
 
 m, n = A.shape
-x, y, s, D, DT = diffcp.solve_and_derivative(A, b, c, cone_dict)
+x, y, s, D, DT = diffcp.solve_and_derivative(A, b, c, cone_dict, solve_method="Clarabel")
 
 # evaluate the derivative
 nonzeros = A.nonzero()

setup.py

@@ -102,7 +102,9 @@ def is_platform_mac():
         "scipy >= 1.1.0",
         "pybind11 >= 2.4",
         "threadpoolctl >= 1.1",
-        "ecos >= 2.0.10"],
+        "ecos >= 2.0.10",
+        "clarabel >= 0.5.1"
+        ],
     url="http://github.com/cvxgrp/diffcp/",
     ext_modules=ext_modules,
     license="Apache License, Version 2.0",

tests/test_clarabel.py

@@ -0,0 +1,122 @@
+import cvxpy as cp
+import numpy as np
+
+import diffcp.cone_program as cone_prog
+import diffcp.utils as utils
+
+
+def test_solve_and_derivative():
+    np.random.seed(0)
+    m = 20
+    n = 10
+
+    A, b, c, cone_dims = utils.least_squares_eq_scs_data(m, n)
+    for mode in ["lsqr", "dense"]:
+        x, y, s, derivative, adjoint_derivative = cone_prog.solve_and_derivative(
+            A, b, c, cone_dims, mode=mode, solve_method="Clarabel")
+
+        dA = utils.get_random_like(
+            A, lambda n: np.random.normal(0, 1e-6, size=n))
+        db = np.random.normal(0, 1e-6, size=b.size)
+        dc = np.random.normal(0, 1e-6, size=c.size)
+
+        dx, dy, ds = derivative(dA, db, dc)
+
+        x_pert, y_pert, s_pert, _, _ = cone_prog.solve_and_derivative(
+            A + dA, b + db, c + dc, cone_dims, solve_method="Clarabel")
+
+        np.testing.assert_allclose(x_pert - x, dx, atol=1e-8)
+        np.testing.assert_allclose(y_pert - y, dy, atol=1e-8)
+        np.testing.assert_allclose(s_pert - s, ds, atol=1e-8)
+
+        objective = c.T @ x
+        dA, db, dc = adjoint_derivative(
+            c, np.zeros(y.size), np.zeros(s.size))
+
+        x_pert, _, _, _, _ = cone_prog.solve_and_derivative(
+            A + 1e-6 * dA, b + 1e-6 * db, c + 1e-6 * dc, cone_dims, solve_method="Clarabel")
+        objective_pert = c.T @ x_pert
+
+        np.testing.assert_allclose(
+            objective_pert - objective,
+            1e-6 * dA.multiply(dA).sum() + 1e-6 * db @ db + 1e-6 * dc @ dc, atol=1e-8)
+
+
+def test_threading():
+    np.random.seed(0)
+    test_rtol = 1e-3
+    test_atol = 1e-8
+    m = 20
+    n = 10
+    As, bs, cs, cone_dicts = [], [], [], []
+    results = []
+
+    for _ in range(50):
+        A, b, c, cone_dims = utils.least_squares_eq_scs_data(m, n)
+        As += [A]
+        bs += [b]
+        cs += [c]
+        cone_dicts += [cone_dims]
+        results.append(cone_prog.solve_and_derivative(A, b, c, cone_dims))
+
+    for n_jobs in [1, -1]:
+        xs, ys, ss, _, DT_batch = cone_prog.solve_and_derivative_batch(
+            As, bs, cs, cone_dicts, n_jobs_forward=n_jobs, n_jobs_backward=n_jobs)
+
+        for i in range(50):
+            np.testing.assert_allclose(results[i][0], xs[i], rtol=test_rtol, atol=test_atol)
+            np.testing.assert_allclose(results[i][1], ys[i], rtol=test_rtol, atol=test_atol)
+            np.testing.assert_allclose(results[i][2], ss[i], rtol=test_rtol, atol=test_atol)
+
+        dAs, dbs, dcs = DT_batch(xs, ys, ss)
+        for i in range(50):
+            dA, db, dc = results[
+                i][-1](results[i][0], results[i][1], results[i][2])
+            np.testing.assert_allclose(dA.todense(), dAs[i].todense(), rtol=test_rtol, atol=test_atol)
+            np.testing.assert_allclose(dbs[i], db, rtol=test_rtol, atol=test_atol)
+            np.testing.assert_allclose(dcs[i], dc, rtol=test_rtol, atol=test_atol)
+
+
+def test_expcone():
+    np.random.seed(0)
+    n = 10
+    y = cp.Variable(n)
+    obj = cp.Minimize(- cp.sum(cp.entr(y)))
+    const = [cp.sum(y) == 1]
+    prob = cp.Problem(obj, const)
+    A, b, c, cone_dims = utils.scs_data_from_cvxpy_problem(prob)
+    for mode in ["lsqr", "lsmr", "dense"]:
+        x, y, s, D, DT = cone_prog.solve_and_derivative(
+            A,
+            b,
+            c,
+            cone_dims,
+            solve_method="Clarabel",
+            mode=mode,
+            tol_gap_abs=1e-13,
+            tol_gap_rel=1e-13,
+            tol_feas=1e-13,
+            tol_ktratio=1e-13,
+        )
+        dA = utils.get_random_like(A, lambda n: np.random.normal(0, 1e-6, size=n))
+        db = np.random.normal(0, 1e-6, size=b.size)
+        dc = np.random.normal(0, 1e-6, size=c.size)
+        dx, dy, ds = D(dA, db, dc)
+        x_pert, y_pert, s_pert, _, _ = cone_prog.solve_and_derivative(
+            A + dA,
+            b + db,
+            c + dc,
+            cone_dims,
+            solve_method="Clarabel",
+            mode=mode,
+            tol_gap_abs=1e-13,
+            tol_gap_rel=1e-13,
+            tol_feas=1e-13,
+            tol_infeas_abs=1e-13,
+            tol_infeas_rel=1e-13,
+            tol_ktratio=1e-13,
+        )
+
+        np.testing.assert_allclose(x_pert - x, dx, atol=1e-8)
+        np.testing.assert_allclose(y_pert - y, dy, atol=1e-8)
+        np.testing.assert_allclose(s_pert - s, ds, atol=1e-8)

tests/test_ecos.py

@@ -87,6 +87,7 @@ def test_expcone():
                                                                       feastol=1e-10,
                                                                       abstol=1e-10,
                                                                       reltol=1e-10)
+
         np.testing.assert_allclose(x_pert - x, dx, atol=1e-8)
         np.testing.assert_allclose(y_pert - y, dy, atol=1e-8)
         np.testing.assert_allclose(s_pert - s, ds, atol=1e-8)

tests/test_scs.py

@@ -119,6 +119,7 @@ def test_expcone():
                                                                       solve_method="SCS",
                                                                       mode=mode,
                                                                       eps=1e-10)
+
         np.testing.assert_allclose(x_pert - x, dx, atol=1e-8)
         np.testing.assert_allclose(y_pert - y, dy, atol=1e-8)
-        np.testing.assert_allclose(s_pert - s, ds, atol=1e-8)
+        np.testing.assert_allclose(s_pert - s, ds, atol=1e-8)