Added method to HomQCQP for extracting the standard form of a problem for use with SCS, Clarabel, etc.

Added a function to run clarabel solver.
Added functions to convert between vectorized form of PSD matrices.

Author: holmesco

_scripts/rot_loop_admm.py

@@ -434,7 +434,7 @@ def get_junction_tree(self, plot=False):
             var2 = G.vs["name"][cliques[i][1]]
             clique_obj += [self.get_clique_obj(var1, var2, i)]
             for j in range(i + 1, len(cliques)):
-                # Get seperator set for list
+                # Get separator set for list
                 sepset = set(cliques[i]) & set(cliques[j])
                 if len(sepset) > 0:
                     junction.add_edge(i, j, weight=-len(sepset), sepset=sepset)

_test/test_homqcqp.py

@@ -6,8 +6,10 @@
 from poly_matrix import PolyMatrix
 
 from cert_tools import HomQCQP
+from cert_tools.linalg_tools import smat, svec
 from cert_tools.problems.rot_synch import RotSynchLoopProblem
-from cert_tools.sparse_solvers import solve_dsdp
+from cert_tools.sdp_solvers import solve_sdp_homqcqp
+from cert_tools.sparse_solvers import solve_clarabel, solve_dsdp
 
 
 def get_chain_rot_prob(N=10):
@@ -20,33 +22,32 @@ def get_loop_rot_prob():
 
 class TestHomQCQP(unittest.TestCase):
 
-    def test_solve():
+    def test_solve(self):
         # Create chain of rotations problem:
         problem = get_chain_rot_prob()
         assert isinstance(problem, HomQCQP), TypeError(
             "Problem should be homogenized qcqp object"
         )
         # Solve SDP via standard method
-        X, info, time = problem.solve_sdp(verbose=True)
+        X, info, time = solve_sdp_homqcqp(problem, verbose=True)
         # Convert solution
         R = problem.convert_sdp_to_rot(X)
 
-    def test_get_asg(plot=False):
-        """Test retreival of aggregate sparsity graph"""
+    def test_get_asg(self, plot=False):
+        """Test retrieval of aggregate sparsity graph"""
         # Test on chain graph
         problem = get_chain_rot_prob()
-        problem.get_asg(rm_homog=True)  # Get Graph
-        problem.triangulate_graph()  # Triangulate graph
+        problem.clique_decomposition()
         # No fill in expected
-        assert any(problem.asg.es["fill_edge"]) is False
+        assert problem.symb.fill[0] == 0, ValueError("Expected no fill in")
         if plot:
             problem.plot_asg()
 
         # Test on loop graph
         problem = get_loop_rot_prob()
-        problem.get_asg(rm_homog=True)
-        problem.triangulate_graph()
-        assert any(problem.asg.es["fill_edge"]) is True
+        problem.get_asg(rm_homog=False)
+        problem.clique_decomposition()
+        assert problem.symb.fill[0] > 0, ValueError("Expected fill in")
         if plot:
             problem.plot_asg()
 
@@ -77,8 +78,8 @@ def test_clique_decomp(self, rm_homog=False, plot=False):
             parent = problem.cliques[clique.parent]
 
             vertices = list(parent.var_sizes.keys()) + list(clique.var_sizes.keys())
-            assert set(clique.seperator).issubset(vertices), ValueError(
-                "seperator set should be in set of involved clique vertices"
+            assert set(clique.separator).issubset(vertices), ValueError(
+                "separator set should be in set of involved clique vertices"
             )
 
         # Check that mapping from variables to cliques is correct
@@ -98,7 +99,6 @@ def test_consistency_constraints(self):
         # Test chain topology
         nvars = 5
         problem = get_chain_rot_prob(N=nvars)
-        problem.get_asg(rm_homog=False)  # Get Graph
         problem.clique_decomposition()  # Run clique decomposition
         eq_list = problem.get_consistency_constraints()
 
@@ -115,11 +115,11 @@ def test_consistency_constraints(self):
         x_list, info = solve_dsdp(problem, verbose=True, tol=1e-8)
         # Verify that the clique variables are equal on overlaps
         for l, clique_l in enumerate(problem.cliques):
-            # seperator
-            sepset = clique_l.seperator
+            # separator
+            sepset = clique_l.separator
             if len(sepset) == 0:  # skip the root clique
                 continue
-            # fet parent clique and seperator set
+            # fet parent clique and separator set
             k = clique_l.parent
             clique_k = problem.cliques[k]
 
@@ -135,7 +135,6 @@ def test_decompose_matrix(self):
         # setup
         nvars = 5
         problem = get_chain_rot_prob(N=nvars)
-        problem.get_asg(rm_homog=False)  # Get Graph
         problem.clique_decomposition()  # get clique decomposition
         C = problem.C
 
@@ -173,10 +172,9 @@ def test_solve_dsdp(self):
         # Test chain topology
         nvars = 5
         problem = get_chain_rot_prob(N=nvars)
-        problem.get_asg(rm_homog=False)  # Get agg sparse graph
         problem.clique_decomposition()  # get cliques
         # Solve non-decomposed problem
-        X, info, time = problem.solve_sdp(verbose=True)
+        X, info, time = solve_sdp_homqcqp(problem, verbose=True)
         # get cliques from non-decomposed solution
         c_list_nd = problem.get_cliques_from_psd_mat(X)
         # Solve decomposed problem (Interior Point Version)
@@ -208,12 +206,52 @@ def test_solve_dsdp(self):
             err_msg="Completed and non-decomposed solutions differ",
         )
 
+    def test_standard_form(self):
+        """Test that the standard form problem definition is correct"""
+        nvars = 2
+        problem = get_chain_rot_prob(N=nvars)
+        problem.get_asg()
+        P, q, A, b = problem.get_standard_form()
+
+        # get solution from MOSEK
+        X, info, time = solve_sdp_homqcqp(problem, verbose=True)
+        x = svec(X)
+
+        # Check cost matrix
+        cost = np.dot(b, x)
+        np.testing.assert_allclose(
+            cost, info["cost"], atol=1e-12, err_msg="Cost incorrect"
+        )
+        # Check constraints
+        for i, vec in enumerate(A.T):
+            a = vec.toarray().squeeze(0)
+            value = np.dot(a, x)
+            np.testing.assert_allclose(
+                value, -q[i], atol=1e-10, err_msg=f"Constraint {i} has violation"
+            )
+
+    def test_clarabel(self):
+        nvars = 2
+        problem = get_chain_rot_prob(N=nvars)
+        problem.get_asg()
+        X_clarabel = solve_clarabel(problem)
+        X, info, time = solve_sdp_homqcqp(problem, verbose=True)
+
+        np.testing.assert_allclose(
+            X_clarabel,
+            X,
+            atol=1e-9,
+            err_msg="Clarabel and MOSEK solutions differ",
+        )
+
 
 if __name__ == "__main__":
     test = TestHomQCQP()
     # test.test_solve()
     # test.test_get_asg(plot=True)
     # test.test_clique_decomp(plot=False)
-    # test.test_consistency_constraints()
+    test.test_consistency_constraints()
     # test.test_decompose_matrix()
-    test.test_solve_dsdp()
+    # test.test_solve_dsdp()
+    # test.test_standard_form()
+    # test.test_clarabel()

_test/test_linalg_tools.py

@@ -0,0 +1,68 @@
+import unittest
+
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy.sparse as sp
+
+from cert_tools.linalg_tools import smat, svec
+
+
+class TestLinAlg(unittest.TestCase):
+
+    def get_psd_mat(self, n=3):
+        # get random symmetric matrix
+        A = np.random.random((n, n))
+        U, S, V = np.linalg.svd(A)
+        S = np.abs(S)
+        return (U * S) @ U.T
+
+    def test_svec(self):
+        # fix seed
+        np.random.seed(0)
+        # generate matrices
+        S1 = self.get_psd_mat()
+        S2 = self.get_psd_mat()
+        # Vectorize
+        s1 = svec(S1)
+        s2 = svec(S2)
+        # test mapping
+        np.testing.assert_almost_equal(smat(s1), S1)
+        np.testing.assert_almost_equal(smat(s2), S2)
+        # products should be equal
+        prod_mat = np.trace(S1 @ S2)
+        prod_vec = np.dot(s1, s2)
+        assert prod_mat == prod_vec, "PSD Inner product not equal"
+
+    def test_svec_sparse(self):
+        # fix seed
+        np.random.seed(0)
+        # generate matrices
+        S1_dense = self.get_psd_mat()
+        S2_dense = self.get_psd_mat()
+        # Remove element to make sure still works when not dense
+        # S1_dense[4, 5] = 0.0
+        # S1_dense[5, 4] = 0.0
+        S1 = sp.csc_matrix(S1_dense)
+        S2 = sp.csc_matrix(S2_dense)
+        S1.eliminate_zeros()
+        S2.eliminate_zeros()
+        # Vectorize
+        s1 = svec(S1)
+        s2 = svec(S2)
+        s1_dense = svec(S1_dense)
+        s2_dense = svec(S2_dense)
+        np.testing.assert_almost_equal(s1_dense, s1.toarray().squeeze(0))
+        np.testing.assert_almost_equal(s2_dense, s2.toarray().squeeze(0))
+        # test mapping
+        np.testing.assert_almost_equal(smat(s1), S1.toarray())
+        np.testing.assert_almost_equal(smat(s2), S2.toarray())
+        # products should be equal
+        prod_mat = np.trace(S1.toarray() @ S2.toarray())
+        prod_vec = (s1 @ s2.T).toarray()
+        assert prod_mat == prod_vec, "PSD Inner product not equal"
+
+
+if __name__ == "__main__":
+    test = TestLinAlg()
+    test.test_svec()
+    test.test_svec_sparse()

cert_tools/base_clique.py

@@ -27,25 +27,25 @@ def __init__(
         index,
         var_sizes,
         parent,
-        seperator,
+        separator,
         Q: PolyMatrix = PolyMatrix(),
         A_list=[],
         b_list=[],
     ):
         self.index = index
-        if var_sizes is not None:
-            assert "h" in var_sizes, f"Each clique must have a homogenizing variable"
+        # if var_sizes is not None:
+        #     assert "h" in var_sizes, f"Each clique must have a homogenizing variable"
         self.var_sizes = var_sizes
         self.var_list = list(self.var_sizes.keys())
         self.var_inds, self.size = self._get_start_indices()
         # Store clique tree information
-        for key in seperator:
+        for key in separator:
             assert (
                 key in self.var_sizes.keys()
-            ), f"seperator element {key} not contained in clique"
-        self.seperator = seperator  # seperator set between this clique and its parent
+            ), f"separator element {key} not contained in clique"
+        self.separator = separator  # separator set between this clique and its parent
         self.residual = self.var_list.copy()
-        for varname in seperator:
+        for varname in separator:
             self.residual.remove(varname)
         self.parent = parent  # index of the parent clique
         self.children = set()  # set of children of this clique in the clique tree
@@ -81,7 +81,7 @@ def _get_indices(self, var_list):
             _type_: _description_
         """
         if type(var_list) is not list:
-            var_list = list(var_list)
+            var_list = [var_list]
         # Get index slices for the rows
         slices = []
         for varname in var_list:
@@ -96,10 +96,10 @@ def get_slices(self, mat, var_list_row, var_list_col=[]):
         If one list provided then slices are assumed to be symmetric. If two lists are provided, they are interpreted as the row and column lists, respectively.
         """
         # Get index slices for the rows
-        inds1 = self._get_indices(var_list_col)
+        inds1 = self._get_indices(var_list_row)
         # Get index slices for the columns
         if len(var_list_col) > 0:
-            inds2 = self._get_indices(var_list_row)
+            inds2 = self._get_indices(var_list_col)
         else:
             # If not defined use the same list as rows
             inds2 = inds1