Implement semidefinite constraints in SCS

1 files changed, 22 insertions, 12 deletions

diff --git a/sumofsquares/solver/scs.py b/sumofsquares/solver/scs.py
index a476700..71c2248 100644
--- a/sumofsquares/solver/scs.py
+++ b/sumofsquares/solver/scs.py
@@ -37,7 +37,7 @@ def vec(m: NDArray) -> NDArray:
     """
     n, _ = m.shape
     # Since it is symmetric we only work with the lower triangular part
-    scaling = (np.tril(n, k=-1) * np.sqrt(2) + np.eye(n))
+    scaling = np.tril(np.ones(n), k=-1) * np.sqrt(2) + np.eye(n)
     v = (m * scaling)[np.tril_indices(n)].reshape((-1, 1))
     assert v.shape[0] == n * (n + 1) // 2, "SCS Matrix vectorization is incorrect!"
     return v
@@ -57,7 +57,7 @@ def mat(v: NDArray) -> NDArray:
     m = np.zeros((n, n))
     m[np.tril_indices(n)] = v
 
-    scaling = (np.tril(n, k=-1) / np.sqrt(2) + np.eye(n))
+    scaling = (np.tril(np.ones(n), k=-1) / np.sqrt(2) + np.eye(n))
     m = m * scaling
     m = m + m.T - np.diag(np.diag(m))
     return m
@@ -87,11 +87,17 @@ def solve_cone(prob: ConicProblem, verbose: bool = False,
     #   ep*  dual exponential cone
     #   p    power cone
     #   p*   dual power cone
+    #
+    # In their documentation they call the linear coefficient c. Here we use q.
 
     P, q = None, None
     A_rows, b_rows = [], []
 
-    # for (linear, constant) in prob.constrains["z"]:
+    if prob.q is not None:
+        q = prob.q
+
+    if prob.constraints["z"]:
+        raise NotImplementedError
 
     if prob.constraints["l"]:
         raise NotImplementedError
@@ -102,13 +108,15 @@ def solve_cone(prob: ConicProblem, verbose: bool = False,
     if prob.constraints["q"]:
         raise NotImplementedError
 
-    if prob.constraints["s"]:
-        raise NotImplementedError
+    for (linear, constant) in prob.constraints["s"]:
+        # -1 because RHS
+        b_rows.append(vec(constant) * -1)
+        A_rows.append(np.hstack(tuple(vec(m) for m in linear)))
 
-    if prob.constraints["ep"]:
+    if prob.constraints["e"]:
         raise NotImplementedError
 
-    if prob.constraints["ep*"]:
+    if prob.constraints["e*"]:
         raise NotImplementedError
 
     if prob.constraints["p"]:
@@ -117,15 +125,16 @@ def solve_cone(prob: ConicProblem, verbose: bool = False,
     if prob.constraints["p*"]:
         raise NotImplementedError
 
-    assert len(A_rows) > 0, "Problem is unconstrained! Something is very wrong"
+    assert q is None or P is None or len(A_rows) > 0, "Problem is unconstrained! Something is very wrong"
 
-    A = np.hstack(A_rows)
-    b = np.hstack(b_rows)
+    A = np.vstack(A_rows)
+    b = np.vstack(b_rows)
 
     data = {
         "P": csc_matrix(P) if P is not None else None,
-        "q": q,
-        "A": csc_matrix(A), "b": b,
+        "c": q.reshape((-1,)),
+        "A": csc_matrix(A),
+        "b": b.reshape((-1,)),
     }
 
     cone = {
@@ -155,6 +164,7 @@ def solve_cone(prob: ConicProblem, verbose: bool = False,
     results, i = {}, 0
     for variable, shape in prob.variables.items():
         num_indices = math.prod(shape)
+        # FIXME: need to restore shape using mat()
         values = np.array(sol["x"][i:i+num_indices]).reshape(shape)
         if values.shape == (1, 1):
             values = values[0, 0]