diff options
| -rw-r--r-- | sumofsquares/solver/cvxopt.py | 2 | ||||
| -rw-r--r-- | sumofsquares/solver/scs.py | 35 |
2 files changed, 14 insertions, 23 deletions
diff --git a/sumofsquares/solver/cvxopt.py b/sumofsquares/solver/cvxopt.py index 11b7797..ef5aaaf 100644 --- a/sumofsquares/solver/cvxopt.py +++ b/sumofsquares/solver/cvxopt.py @@ -84,7 +84,7 @@ def solve_cone(prob: ConicProblem, verbose: bool = False, # pprint(f"{b_rows = }") if A_rows: - # b is multiplied by -1 because it is on the RHS + # multiplied by -1 because it is on the RHS b = cvxopt.matrix(np.vstack(b_rows)) A = cvxopt.matrix(np.vstack(A_rows)) * -1 diff --git a/sumofsquares/solver/scs.py b/sumofsquares/solver/scs.py index 1e474da..662f78e 100644 --- a/sumofsquares/solver/scs.py +++ b/sumofsquares/solver/scs.py @@ -29,38 +29,29 @@ class SCSInfo(SolverInfo, UserDict): raise KeyError(f"Key {key} was not found") return self[key] - -def vec(m: NDArray) -> NDArray: +def vec(S: NDArray) -> NDArray: r""" Vectorize a symmetric matrix for SCS by storing only a half of the entries and multiplying the off-diaognal elements by :math:`\sqrt{2}`. """ - n, _ = m.shape - # Since it is symmetric we only work with the lower triangular part - 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 + n = S.shape[0] + S = np.copy(S) + S *= np.sqrt(2) + S[range(n), range(n)] /= np.sqrt(2) + return S[np.triu_indices(n)].reshape((-1,1)) -def mat(v: NDArray) -> NDArray: +def mat(s: NDArray) -> NDArray: r""" Reconstruct a symmetric matrix from a vector that was created using :py:fn:`vec`. This scales the off-diagonal entries by :math:`1/\sqrt{2}`. """ - n_float = .5 * (-1 + np.sqrt(1 + 8 * v.shape[0])) - n = int(n_float) - if not math.isclose(n, n_float): - raise ValueError("To construct a n x n symmetric matrix, the vector " - "must be an n * (n + 1) / 2 dimensional column.") - - m = np.zeros((n, n)) - m[np.tril_indices(n)] = v - - 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 + n = int((np.sqrt(8 * len(s) + 1) - 1) / 2) + S = np.zeros((n, n)) + S[np.triu_indices(n)] = s / np.sqrt(2) + S = S + S.T + S[range(n), range(n)] /= np.sqrt(2) + return S def solve_cone(prob: ConicProblem, verbose: bool = False, |