From c99f304114701ae7e9735bec775e4db18b77130a Mon Sep 17 00:00:00 2001 From: Nao Pross Date: Fri, 10 May 2024 15:51:08 +0200 Subject: Move cvxopt module into submodule for solvers --- sumofsquares/cvxopt.py | 176 ---------------------------------------- sumofsquares/problems.py | 2 +- sumofsquares/solver/__init__.py | 3 + sumofsquares/solver/cvxopt.py | 176 ++++++++++++++++++++++++++++++++++++++++ 4 files changed, 180 insertions(+), 177 deletions(-) delete mode 100644 sumofsquares/cvxopt.py create mode 100644 sumofsquares/solver/__init__.py create mode 100644 sumofsquares/solver/cvxopt.py diff --git a/sumofsquares/cvxopt.py b/sumofsquares/cvxopt.py deleted file mode 100644 index 6a1339f..0000000 --- a/sumofsquares/cvxopt.py +++ /dev/null @@ -1,176 +0,0 @@ -import cvxopt -import numpy as np - -from collections import UserDict -from numpy.typing import NDArray -from dataclasses import dataclass, fields -from pprint import pprint - -import polymatrix as poly -from polymatrix.expression.expression import Expression -from polymatrix.polymatrix.index import MonomialIndex, VariableIndex - -from .abc import Problem, SolverInfo -from .constraints import NonNegative, EqualToZero, PositiveSemiDefinite, ExponentialCone -from .variable import OptVariable - - -class CVXOPTInfo(SolverInfo, UserDict): - """ Dictionary returned by CVXOPT. """ - - def __getattr__(self, attr): - key = " ".join(attr.split("_")) - if key not in self.keys(): - raise KeyError(f"Key {key} was not found") - return self[key] - - -def solve_sos_cone(prob: Problem, verbose: bool = False, *args, **kwargs) -> tuple[dict[OptVariable, float], CVXOPTInfo]: - r""" - Solve a conic problem in the cone of SOS polynomials :math:`\mathbf{\Sigma}_d(x)`. - - Any `*args` and `**kwargs` other than `prob` are passed to the CVXOPT solver. - """ - # Check that the problem can be solved by CVXOpt, i.e. - # - Cost function is at most quadratic - # - Constraints are at most linear - - cost = poly.to_affine(prob.cost.read(prob.state)) - if cost.degree > 2: - raise ValueError("CVXOpt can solve at most quadratic conic programs, " - f"but the given problem has degree {cost.degree} > 2.") - - is_qp = (cost.degree == 2) - - # We need to convert problem to the primal canonical form (given in the - # CVXOPT doc) with x = vstack(prob.variables) - # - # minimize .5 * x.T @ P @ x + q.T @ x - # - # subject to G @ x + s = h - # A @ x = b - # s >= 0 - # - # If the problem is an LP, the CVXOPT doc call the coefficient in the cost - # functions `c`, here we will use `q` for both LP and QP. - - A_rows, b_rows = [], [] - G_rows, h_rows= [], [] - l_dims, q_dims, s_dims = [], [], [] - - # Collect variable indices in dictionary and sort them by value - # i.e. we want variable_indices.values() to be sorted - variable_indices: dict[OptVariable, VariableIndex] = dict(sorted({ - # FIXME there can be more than one indices associated to v - v: next(prob.state.get_indices_as_variable_index(v)) - for v in prob.variables - }.items(), key=lambda item: item[1])) - - # cost linear term - cost_coeffs = dict(cost.affine_coefficients_by_degrees(variable_indices.values())) - q = cost_coeffs[1].T - - # cost quadratic term - if is_qp: - n = len(prob.variables) - P = np.zeros((n, n)) - # This works because of how monomial indices are sorted - # see also polymatrix.index.MonomialIndex.__lt__ - P[np.triu_indices(n)] = cost_coeffs[2] - - # Make symmetric. There is no 0.5 factor because there it is already - # there in the canonical form of CVXOPT - P = (P + P.T) - - x = poly.v_stack(prob.polynomial_variables) - for c in prob.constraints: - if isinstance(c, PositiveSemiDefinite): - raise NotImplementedError - - elif isinstance(c, ExponentialCone): - raise NotImplementedError - - elif isinstance(c, EqualToZero | NonNegative): - # Convert constraints to affine expression in decision variables - constr = poly.to_affine(c.expression.quadratic_in(x).read(prob.state)) - if constr.degree > 1: - raise ValueError("CVXOpt can solve problem with constraints that are " - f"at most linear, but {str(c.expression)} has degree {constr.degree}.") - - # For each constraint we add PSD matrix - nrows, ncols = constr.shape - s_dims.append(nrows) - - # Constant term - c_coeff = constr.affine_coefficient(MonomialIndex.constant()) - c_stacked = c_coeff.T.reshape((nrows * ncols,)) - - # Linear term - l_stacked_rows = [] - for v in variable_indices.values(): - # Get the affine coefficient associated to this variable - linear = MonomialIndex((v,)) - v_coeff = constr.affine_coefficient(linear) - - # stack the columns of coeff matrix into a fat row vector (column major order). - v_stacked = v_coeff.T.reshape((nrows * ncols,)) - l_stacked_rows.append(v_stacked) - - l_stacked = np.vstack(l_stacked_rows) - - if isinstance(c, EqualToZero): - A_rows.append(l_stacked) - # Factor of -1 because it is on the right hand side - b_rows.append(-1 * c_stacked) - - elif isinstance(c, NonNegative): - if c.domain: - raise ValueError("Cannot convert non-negativity constraint to CVXOPT format. " - "Domain restricted non-negativity must be " - "preprocessed by using a Positivstellensatz!") - - # Factor of -1 because it is on the right hand side - G_rows.append(-1 * l_stacked) - h_rows.append(c_stacked) - - else: - raise NotImplementedError(f"Cannot convert constraint of type {type(c)} " - "into the canonical form of CVXOPT.") - - # Convert to CVXOPT matrices - q = cvxopt.matrix(q) - G = cvxopt.matrix(np.hstack(G_rows).T) if G_rows else None - h = cvxopt.matrix(np.hstack(h_rows)) if h_rows else None - A = cvxopt.matrix(np.hstack(A_rows).T) if A_rows else None - b = cvxopt.matrix(np.hstack(b_rows)) if b_rows else None - dims = {'l': l_dims, 'q': q_dims, 's': s_dims} - - # print(q) - # print(G) - # print(h) - # print(A) - # print(b) - # print(dims) - - cvxopt.solvers.options['show_progress'] = verbose - - if not any((G, h, A, b)): - raise ValueError("Optimization problem is unconstrained!") - - # Add matrices for debugging - # info = {} - info = {"q": q, "G": G, "h": h, "A": A, "b": b, "dims": dims} - - if is_qp: - P = cvxopt.matrix(P) - info |= cvxopt.solvers.coneqp(P=P, q=q, G=G, h=h, A=A, b=b, - dims=dims, *args, **kwargs) - else: - info |= cvxopt.solvers.conelp(c=q, G=G, h=h, A=A, b=b, - *args, **kwargs) - - results = {} - for i, v in enumerate(variable_indices.keys()): - results[v] = info['x'][i] - - return results, CVXOPTInfo(info) diff --git a/sumofsquares/problems.py b/sumofsquares/problems.py index 309b13c..14a9cd9 100644 --- a/sumofsquares/problems.py +++ b/sumofsquares/problems.py @@ -19,7 +19,7 @@ from polymatrix.variable.abc import Variable from .abc import Problem, Constraint, Solver, Result from .constraints import NonNegative -from .cvxopt import solve_sos_cone as cvxopt_solve_sos_cone +from .solver.cvxopt import solve_sos_cone as cvxopt_solve_sos_cone from .utils import partition from .variable import OptVariable, from_name as opt_variable_from_name diff --git a/sumofsquares/solver/__init__.py b/sumofsquares/solver/__init__.py new file mode 100644 index 0000000..68a2090 --- /dev/null +++ b/sumofsquares/solver/__init__.py @@ -0,0 +1,3 @@ +""" +This module contains the interfaces to various solvers. +""" diff --git a/sumofsquares/solver/cvxopt.py b/sumofsquares/solver/cvxopt.py new file mode 100644 index 0000000..0de8dea --- /dev/null +++ b/sumofsquares/solver/cvxopt.py @@ -0,0 +1,176 @@ +import cvxopt +import numpy as np + +from collections import UserDict +from numpy.typing import NDArray +from dataclasses import dataclass, fields +from pprint import pprint + +import polymatrix as poly +from polymatrix.expression.expression import Expression +from polymatrix.polymatrix.index import MonomialIndex, VariableIndex + +from ..abc import Problem, SolverInfo +from ..constraints import NonNegative, EqualToZero, PositiveSemiDefinite, ExponentialCone +from ..variable import OptVariable + + +class CVXOPTInfo(SolverInfo, UserDict): + """ Dictionary returned by CVXOPT. """ + + def __getattr__(self, attr): + key = " ".join(attr.split("_")) + if key not in self.keys(): + raise KeyError(f"Key {key} was not found") + return self[key] + + +def solve_sos_cone(prob: Problem, verbose: bool = False, *args, **kwargs) -> tuple[dict[OptVariable, float], CVXOPTInfo]: + r""" + Solve a conic problem in the cone of SOS polynomials :math:`\mathbf{\Sigma}_d(x)`. + + Any `*args` and `**kwargs` other than `prob` are passed to the CVXOPT solver. + """ + # Check that the problem can be solved by CVXOpt, i.e. + # - Cost function is at most quadratic + # - Constraints are at most linear + + cost = poly.to_affine(prob.cost.read(prob.state)) + if cost.degree > 2: + raise ValueError("CVXOpt can solve at most quadratic conic programs, " + f"but the given problem has degree {cost.degree} > 2.") + + is_qp = (cost.degree == 2) + + # We need to convert problem to the primal canonical form (given in the + # CVXOPT doc) with x = vstack(prob.variables) + # + # minimize .5 * x.T @ P @ x + q.T @ x + # + # subject to G @ x + s = h + # A @ x = b + # s >= 0 + # + # If the problem is an LP, the CVXOPT doc call the coefficient in the cost + # functions `c`, here we will use `q` for both LP and QP. + + A_rows, b_rows = [], [] + G_rows, h_rows= [], [] + l_dims, q_dims, s_dims = [], [], [] + + # Collect variable indices in dictionary and sort them by value + # i.e. we want variable_indices.values() to be sorted + variable_indices: dict[OptVariable, VariableIndex] = dict(sorted({ + # FIXME there can be more than one indices associated to v + v: next(prob.state.get_indices_as_variable_index(v)) + for v in prob.variables + }.items(), key=lambda item: item[1])) + + # cost linear term + cost_coeffs = dict(cost.affine_coefficients_by_degrees(variable_indices.values())) + q = cost_coeffs[1].T + + # cost quadratic term + if is_qp: + n = len(prob.variables) + P = np.zeros((n, n)) + # This works because of how monomial indices are sorted + # see also polymatrix.index.MonomialIndex.__lt__ + P[np.triu_indices(n)] = cost_coeffs[2] + + # Make symmetric. There is no 0.5 factor because there it is already + # there in the canonical form of CVXOPT + P = (P + P.T) + + x = poly.v_stack(prob.polynomial_variables) + for c in prob.constraints: + if isinstance(c, PositiveSemiDefinite): + raise NotImplementedError + + elif isinstance(c, ExponentialCone): + raise NotImplementedError + + elif isinstance(c, EqualToZero | NonNegative): + # Convert constraints to affine expression in decision variables + constr = poly.to_affine(c.expression.quadratic_in(x).read(prob.state)) + if constr.degree > 1: + raise ValueError("CVXOpt can solve problem with constraints that are " + f"at most linear, but {str(c.expression)} has degree {constr.degree}.") + + # For each constraint we add PSD matrix + nrows, ncols = constr.shape + s_dims.append(nrows) + + # Constant term + c_coeff = constr.affine_coefficient(MonomialIndex.constant()) + c_stacked = c_coeff.T.reshape((nrows * ncols,)) + + # Linear term + l_stacked_rows = [] + for v in variable_indices.values(): + # Get the affine coefficient associated to this variable + linear = MonomialIndex((v,)) + v_coeff = constr.affine_coefficient(linear) + + # stack the columns of coeff matrix into a fat row vector (column major order). + v_stacked = v_coeff.T.reshape((nrows * ncols,)) + l_stacked_rows.append(v_stacked) + + l_stacked = np.vstack(l_stacked_rows) + + if isinstance(c, EqualToZero): + A_rows.append(l_stacked) + # Factor of -1 because it is on the right hand side + b_rows.append(-1 * c_stacked) + + elif isinstance(c, NonNegative): + if c.domain: + raise ValueError("Cannot convert non-negativity constraint to CVXOPT format. " + "Domain restricted non-negativity must be " + "preprocessed by using a Positivstellensatz!") + + # Factor of -1 because it is on the right hand side + G_rows.append(-1 * l_stacked) + h_rows.append(c_stacked) + + else: + raise NotImplementedError(f"Cannot convert constraint of type {type(c)} " + "into the canonical form of CVXOPT.") + + # Convert to CVXOPT matrices + q = cvxopt.matrix(q) + G = cvxopt.matrix(np.hstack(G_rows).T) if G_rows else None + h = cvxopt.matrix(np.hstack(h_rows)) if h_rows else None + A = cvxopt.matrix(np.hstack(A_rows).T) if A_rows else None + b = cvxopt.matrix(np.hstack(b_rows)) if b_rows else None + dims = {'l': l_dims, 'q': q_dims, 's': s_dims} + + # print(q) + # print(G) + # print(h) + # print(A) + # print(b) + # print(dims) + + cvxopt.solvers.options['show_progress'] = verbose + + if not any((G, h, A, b)): + raise ValueError("Optimization problem is unconstrained!") + + # Add matrices for debugging + # info = {} + info = {"q": q, "G": G, "h": h, "A": A, "b": b, "dims": dims} + + if is_qp: + P = cvxopt.matrix(P) + info |= cvxopt.solvers.coneqp(P=P, q=q, G=G, h=h, A=A, b=b, + dims=dims, *args, **kwargs) + else: + info |= cvxopt.solvers.conelp(c=q, G=G, h=h, A=A, b=b, + *args, **kwargs) + + results = {} + for i, v in enumerate(variable_indices.keys()): + results[v] = info['x'][i] + + return results, CVXOPTInfo(info) -- cgit v1.2.1