diff options
Diffstat (limited to '')
| -rw-r--r-- | sumofsquares/__init__.py | 71 | ||||
| -rw-r--r-- | sumofsquares/abc.py | 80 | ||||
| -rw-r--r-- | sumofsquares/canon.py | 130 | ||||
| -rw-r--r-- | sumofsquares/constraints.py | 105 | ||||
| -rw-r--r-- | sumofsquares/expression.py | 90 | ||||
| -rw-r--r-- | sumofsquares/problems.py | 339 | ||||
| -rw-r--r-- | sumofsquares/solver/cvxopt.py | 179 |
7 files changed, 592 insertions, 402 deletions
diff --git a/sumofsquares/__init__.py b/sumofsquares/__init__.py index 6ca2b39..6c02f24 100644 --- a/sumofsquares/__init__.py +++ b/sumofsquares/__init__.py @@ -9,6 +9,7 @@ Module Overview Class diagram: Arrows denote inheritance, diamonds indicate composition. :: + ┏━━sumofsquares package━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ ┃ ┃ ┌─.abc──────────────┐ ┌─.problems────────────────────────────────────────────────────┐ ┃ @@ -18,12 +19,12 @@ Class diagram: Arrows denote inheritance, diamonds indicate composition. ┃ │ │ │ │ ┌───────────────────┐ ┌────────────┐ ┌──────────────────┐ │ ┃ ┃ │ ◇ │ │ │ ConicProblem │ │ SOSProblem │ │InternalSOSProblem│ │ ┃ ┃ │ ┌──────────────┐ │ │ └───────────────────┘ └────────────┘ └──────────────────┘ │ ┃ - ┃ │ │ Solver(Enum) │ │ │ ▲ │ ┃ - ┃ │ └──────────────┘ │ │ ┌────┴────────────────┐ │ ┃ - ┃ │ │ │ │ │ │ ┃ - ┃ │ ┌──────────────┐ │ │ ┌─────────────┐ ┌─────────────────┐ ┌────────────────────┐ │ ┃ - ┃ │ │ Result │◀─┼───┼─│ ConicResult │ │PutinarSOSProblem│ │SchmuedgenSOSProblem│ │ ┃ - ┃ │ └──────────────┘ │ │ └─────────────┘ └─────────────────┘ └────────────────────┘ │ ┃ + ┃ │ │ Solver(Enum) │ │ │ │ ┃ + ┃ │ └──────────────┘ │ │ │ ┃ + ┃ │ │ │ │ ┃ + ┃ │ ┌──────────────┐ │ │ ┌─────────────┐ │ ┃ + ┃ │ │ Result │◀─┼───┼─│ ConicResult │ │ ┃ + ┃ │ └──────────────┘ │ │ └─────────────┘ │ ┃ ┃ │ │ │ └──────────────────────────────────────────────────────────────┘ ┃ ┃ │ │ │ ┃ ┃ │ ◇ │ ┌─.solvers.cvxopt──┐ ┌─.solvers.scs─────┐ ┌─.solvers.mosek───┐ ┃ @@ -50,18 +51,18 @@ Class diagram: Arrows denote inheritance, diamonds indicate composition. ┃ │ │ │ │ ┃ ┃ └───────────────────┘ └───────────────────────────────────────────────────────────────────┘ ┃ ┃ ┃ - ┃ ┌─.variable─────────────────┐ ┃ ┏━━━polymatrix package━━━━━━━━━┓ - ┃ │ │ ┃ ┃ ┃ - ┃ │ ┌───────────┐ │ ┃ ┃ ┌──────────┐ ┌──────────┐ ┃ - ┃ │ │OptVariable│──────────┼──╋───╋─▶│ Variable │ │Expression│ ┃ - ┃ │ └───────────┘ │ ┃ ┃ └──────────┘ └──────────┘ ┃ - ┃ │ ▲ │ ┃ ┃ │ ┃ - ┃ │ │ │ ┃ ┃ ┌─────────┘ ┃ - ┃ │ │ │ ┃ ┃ ◇ ┃ - ┃ │ ┌────────────────┐ │ ┃ ┃ ┌─────────────────────┐ ┃ - ┃ │ │OptVariableMixin│────────┼──╋───╋─▶│ ExpressionBaseMixin │ ┃ - ┃ │ └────────────────┘ │ ┃ ┃ └─────────────────────┘ ┃ - ┃ └───────────────────────────┘ ┃ ┃ ┃ + ┃ ┌─.canon───────────────────────────────────────────────────┐ ┌─.variable─────────────────┐ ┃ ┏━━━polymatrix package━━━━━━━━━┓ + ┃ │ │ │ │ ┃ ┃ ┃ + ┃ │ ┌────────────────────┐ │ │ ┌───────────┐ │ ┃ ┃ ┌──────────┐ ┌──────────┐ ┃ + ┃ │ │ Canonicalization │ │ │ │OptVariable│──────────┼──╋───╋─▶│ Variable │ │Expression│ ┃ + ┃ │ └────────────────────┘ │ │ └───────────┘ │ ┃ ┃ └──────────┘ └──────────┘ ┃ + ┃ │ ▲ │ │ ▲ │ ┃ ┃ │ ┃ + ┃ │ ├────────────────────┐ │ │ │ │ ┃ ┃ ┌─────────┘ ┃ + ┃ │ │ │ │ │ │ │ ┃ ┃ ◇ ┃ + ┃ │ ┌───────────────┐ ┌───────────┐ │ │ ┌────────────────┐ │ ┃ ┃ ┌─────────────────────┐ ┃ + ┃ │ │ PutinarPSatz │ │ LogDet │ │ │ │OptVariableMixin│────────┼──╋───╋─▶│ ExpressionBaseMixin │ ┃ + ┃ │ └───────────────┘ └───────────┘ │ │ └────────────────┘ │ ┃ ┃ └─────────────────────┘ ┃ + ┃ └──────────────────────────────────────────────────────────┘ └───────────────────────────┘ ┃ ┃ ┃ ┗━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┛ ┗━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┛ """ @@ -72,8 +73,9 @@ from polymatrix.expression.expression import Expression from polymatrix.expressionstate import ExpressionState from .abc import Problem, Result, Set, Constraint, Solver +from .canon import Canonicalization, PutinarPSatz from .constraints import BasicSemialgebraicSet, NonNegative -from .problems import PutinarSOSProblem +from .problems import SOSProblem from .utils import partition from .variable import ( OptVariable, @@ -97,41 +99,26 @@ def make_problem( cost: Expression, constraints: Iterable[Constraint] = (), solver: Solver = Solver.CVXOPT, - state: ExpressionState | None = None - ) -> Problem: + psatz: Canonicalization = PutinarPSatz + ) -> SOSProblem: """ Create a sum-of-squares optimization problem. """ - if not state: - state = poly.make_state() - - # Convert to polymatrix object to extract variables - state, _ = cost.apply(state) - for constr in constraints: - state, _ = constr.expression.apply(state) - - # Collect decision variables from state object - def is_optvariable(v): - return isinstance(v, OptVariable) - - poly_variables, opt_variables = partition(is_optvariable, state.indices.keys()) - return PutinarSOSProblem( - cost=cost, - constraints=constraints, - variables=tuple(opt_variables), - polynomial_variables=tuple(poly_variables), - solver=solver, - state=state) + return psatz(SOSProblem(cost, constraints, solver)) def solve_problem( cost: Expression, constraints: Iterable[Constraint] = (), solver: Solver = Solver.CVXOPT, - verbose: bool = False + verbose: bool = False, + state: ExpressionState | None = None ) -> tuple[Problem, Result]: """ Solve a sum-of-squares optimization problem. """ + if state is None: + state = poly.make_state() + prob = make_problem(cost, constraints, solver) return prob, prob.solve(verbose) diff --git a/sumofsquares/abc.py b/sumofsquares/abc.py index f0e95bb..20f4bd6 100644 --- a/sumofsquares/abc.py +++ b/sumofsquares/abc.py @@ -1,101 +1,85 @@ """ Abstract base classes for Sum of Squares package """ from __future__ import annotations -from typing import Sequence, Any from abc import ABC, abstractmethod from enum import Enum, auto - -from polymatrix.expression.expression import Expression -from polymatrix.expressionstate import ExpressionState -from polymatrix.variable import Variable +from typing import Any, Generic, TypeVar from sumofsquares.variable import OptVariable +# ┏━┓┏━┓╻ ╻ ╻┏━╸┏━┓ +# ┗━┓┃ ┃┃ ┃┏┛┣╸ ┣┳┛ +# ┗━┛┗━┛┗━╸┗┛ ┗━╸╹┗╸ + class Solver(Enum): """ Enumeration for the supported solvers. """ CVXOPT = auto() - # MOSEK = auto() + MOSEK = auto() + SCS = auto() class SolverInfo(ABC): """ Type that information returned by a specific solver. """ -class Result(ABC): - """ Result of an optimization problem. """ - @abstractmethod - def value_of(self, var: OptVariable) -> float: - """ Retrieve value of variable. """ +# ┏━╸┏━┓┏┓╻┏━┓╺┳╸┏━┓┏━┓╻┏┓╻╺┳╸┏━┓ +# ┃ ┃ ┃┃┗┫┗━┓ ┃ ┣┳┛┣━┫┃┃┗┫ ┃ ┗━┓ +# ┗━╸┗━┛╹ ╹┗━┛ ╹ ╹┗╸╹ ╹╹╹ ╹ ╹ ┗━┛ class Set(ABC): """ A set. In this context always a subset of :math:`\mathbf{R}^q` """ -class Constraint(ABC): - r""" - Optimization constraint. - """ +E = TypeVar("E") +class Constraint(ABC, Generic[E]): + """ Optimization constraint. """ @property @abstractmethod - def expression(self) -> Expression: + def expression(self) -> E: """ Expression under the constraint. """ def __str__(self): return f"{self.__class__.__qualname__}({str(self.expression)})" -class Problem(ABC): - """ Optimization Problem. """ - - # --- dunder methods --- +# ┏━┓┏━┓┏━┓┏┓ ╻ ┏━╸┏┳┓┏━┓ +# ┣━┛┣┳┛┃ ┃┣┻┓┃ ┣╸ ┃┃┃┗━┓ +# ╹ ╹┗╸┗━┛┗━┛┗━╸┗━╸╹ ╹┗━┛ - def __str__(self): - vs = ", ".join(map(str, self.variables)) - n = len(self.constraints) - s = f"SumOfSquares Problem ({self.__class__.__qualname__}):\n\n " \ - f"find {vs}\n minimize {self.cost}\n " \ - f"subject to ({n} constraint{'s' if n > 1 else ''})\n" - for c in self.constraints: - s += f"\t* {str(c)}\n" +class Result(ABC): + """ Result of an optimization problem. """ + @abstractmethod + def value_of(self, var: OptVariable) -> float: + """ Retrieve value of variable. """ + - return s +class Problem(ABC): + """ Optimization Problem. """ # --- Properties --- @property @abstractmethod - def cost(self) -> Expression: - """ Expression of cost function. """ + def cost(self) -> Any: + """ Cost function of the problem. """ + # The type here is unspecified on purpose, see + # programs module for why is it Any @property @abstractmethod - def constraints(self) -> Sequence[Constraint]: + def constraints(self) -> Any: """ Optimization constraints. """ - - @property - @abstractmethod - def variables(self) -> Sequence[OptVariable]: - """ Optimization variables of the problem. """ - - @property - @abstractmethod - def polynomial_variables(self) -> Sequence[Variable]: - """ Polynomial variables, or varaibles that are not optimization varaibles. """ + # The type here is unspecified on purpose, see + # programs module for why is it Any @property @abstractmethod def solver(self) -> Solver: """ Which solver was / will be used to solve the problem. """ - - @property - @abstractmethod - def state(self) -> ExpressionState: - """ Expression state used to evaluate constraints and costs. """ - # --- Methods --- @abstractmethod diff --git a/sumofsquares/canon.py b/sumofsquares/canon.py index 4918758..e5960e7 100644 --- a/sumofsquares/canon.py +++ b/sumofsquares/canon.py @@ -1,11 +1,131 @@ """ Canonicalization of Problems """ -from .abc import Problem +import math +import polymatrix as poly +from abc import abstractmethod +from dataclassabc import dataclassabc +from dataclasses import replace +from typing import Sequence +from typing_extensions import override -def log_canon(prob: Problem) -> Problem: - raise NotImplementedError +from polymatrix.expression.from_ import Expression, from_statemonad +from polymatrix.expressionstate import ExpressionState +from polymatrix.statemonad import init_state_monad +from polymatrix.polymatrix.mixins import PolyMatrixMixin +from polymatrix.variable import Variable +from polymatrix.expression.init import ( + init_concatenate_expr, + init_diag_expr, + init_lower_triangular_expr, + init_transpose_expr) + +from .abc import Constraint, Result +from .constraints import NonNegative, PositiveSemiDefinite, ExponentialCone +from .error import AlgebraicError +from .problems import SOSProblem, InternalSOSProblem +from .variable import OptVariable, init_opt_variable_expr, from_name as opt_variable_from_name +from .utils import partition + + +class Canonicalization: + """ Transform a problem into another problem. """ + + def __str__(self): + s = f"Apply canonicalization procedure {self.__class__.__qualname__} to\n" + return s + str(self.problem) + + @property + @abstractmethod + def problem(self) -> SOSProblem: + """ The problem that will be canonicalized """ + + def solve(self, verbose: bool = False, state: ExpressionState | None = None) -> Result: + if state is None: + state = poly.make_state() + + state, internal_prob = self.apply(state) + return internal_prob.solve(verbose) + + @abstractmethod + def apply(self, state: ExpressionState) -> tuple[ExpressionState, InternalSOSProblem]: + """ Apply the canonicalization procedure """ + + +@dataclassabc(frozen=True) +class PutinarPSatz(Canonicalization): + problem: SOSProblem + + @staticmethod + def make_multiplier(variables: Sequence[OptVariable], + polynomial_variables: Sequence[Variable], + constr: Expression, + domain_poly: Expression, + multiplier_name: str) -> Expression: + r""" + Make a multiplier polynomial for a domain polynomial of a constraint. + + Since Putinar's PSatz does not explicity tell us the order of the + multiplier this function will make a multiplier, let's call it + :math:`\gamma`, such that :math:`\deg(\gamma * p) = \deg(c)`, wherein + :math:`p` is the domain polynomial and :math:`c` is the contraint + polynomial. + """ + # Number of variables in the original problem + q = sum(math.prod(v.shape) for v in variables) + x = poly.v_stack((1,) + tuple(polynomial_variables)) + + def make_multiplier_later(state: ExpressionState) -> tuple[ExpressionState, PolyMatrixMixin]: + state, constr_pm = constr.expression.degree().apply(state) + state, domain_poly_pm = domain_poly.degree().apply(state) + + constr_deg = constr_pm.at(0, 0).constant() + domain_poly_deg = domain_poly_pm.at(0, 0).constant() + + # degree of multiplier + d = constr_deg - domain_poly_deg + if d < 0: + raise AlgebraicError("Cannot create Positivstellensatz multiplier because " + f"constraint polynomial has degree {constr_deg} and domain " + f"polynomial {domain_poly} has degree {domain_poly_deg}") + + n = sum(math.comb(k + q - 1, k) for k in range(0, d+1)) + + # TODO: check that there is not a variable with this name already + # Coefficients of the multiplier polynomial + c = opt_variable_from_name(multiplier_name, shape=(1,n)) + multiplier = c @ x.combinations(tuple(range(d +1))) + + return multiplier.apply(state) + return poly.from_statemonad(init_state_monad(make_multiplier_later)) + + # TODO: rewrite this function, this is adapted from how it was done + # previously, and it works but it's not efficient. + @abstractmethod + def apply(self, state: ExpressionState) -> tuple[ExpressionState, InternalSOSProblem]: + + def need_positivstellensatz(c: Constraint) -> bool: + return isinstance(c, NonNegative) and (c.domain is not None) + + constraints, to_process = partition(need_positivstellensatz, self.problem.constraints) + constraints = tuple(constraints) # because it is a generator + + state, prob = self.problem.apply(state) + + new_constraints: list[Constraint] = list(constraints) + for i, constr in enumerate(to_process): + new_constr = constr.expression + for domain_poly in constr.domain.polynomials: + m = PutinarPSatz.make_multiplier(prob.variables, + prob.polynomial_variables, + constr, domain_poly, rf"\gamma_{i}") + + # Subtract multiplier from expression and impose that it is also SOS + new_constr -= m * domain_poly + new_constraints.append(NonNegative(m)) + + new_constraints.append(NonNegative(new_constr)) + + return replace(self.problem, constraints=new_constraints).apply(state) -def log_det_canon(prob: Problem) -> Problem: - raise NotImplementedError diff --git a/sumofsquares/constraints.py b/sumofsquares/constraints.py index 71dbb06..3d681c5 100644 --- a/sumofsquares/constraints.py +++ b/sumofsquares/constraints.py @@ -1,30 +1,24 @@ -import math - from dataclasses import dataclass from dataclassabc import dataclassabc from typing_extensions import override -from typing import Iterable - -import polymatrix as poly +from typing import Generic, TypeVar from polymatrix.expression.expression import Expression -from polymatrix.expressionstate import ExpressionState -from polymatrix.polymatrix.abc import PolyMatrix -from polymatrix.statemonad import init_state_monad -from polymatrix.variable import Variable +from polymatrix.polymatrix.mixins import PolyMatrixMixin from .abc import Set, Constraint -from .error import AlgebraicError -from .variable import OptVariable, from_name as opt_variable_from_name + + +# Type variable for expression or polymatrix +E = TypeVar("E", bound=Expression | PolyMatrixMixin) # ┏━┓┏━╸╺┳╸┏━┓ # ┗━┓┣╸ ┃ ┗━┓ # ┗━┛┗━╸ ╹ ┗━┛ - @dataclass(frozen=True) -class BasicSemialgebraicSet(Set): +class BasicSemialgebraicSet(Set, Generic[E]): r""" A set that is described by the intersection of the positive loci of a finite number of polynomials (under some conditions). Roughly speaking this @@ -41,7 +35,7 @@ class BasicSemialgebraicSet(Set): This tuple stores the :math:`g_i(x)`. """ - polynomials: tuple[Expression, ...] + polynomials: tuple[E, ...] def __str__(self): polynomials = ', '.join(map(str, self.polynomials)) @@ -54,9 +48,9 @@ class BasicSemialgebraicSet(Set): @dataclassabc(frozen=True) -class EqualToZero(Constraint): +class EqualToZero(Constraint[E]): """ Constrain an expression to be equal to zero. """ - expression: Expression + expression: E @override def __str__(self): @@ -64,10 +58,10 @@ class EqualToZero(Constraint): @dataclassabc(frozen=True) -class NonNegative(Constraint): +class NonNegative(Constraint[E]): """ - Constrain an expression to be non-negative in a certain domain. If the - domain in set to `None` then it is interpreted as being non-negative + Constrain a scalar expression to be non-negative in a certain domain. If + the domain in set to `None` then it is interpreted as being non-negative everywhere. **Note:** In this package, non-negativitiy constraint will always @@ -75,8 +69,8 @@ class NonNegative(Constraint): be written as a sum-of-squares, hence this constraint can also be understood as `expression` being SOS. """ - expression: Expression - domain: BasicSemialgebraicSet | None = None + expression: E + domain: BasicSemialgebraicSet[E] | None = None @override def __str__(self): @@ -87,72 +81,11 @@ class NonNegative(Constraint): @dataclassabc(frozen=True) -class PositiveSemiDefinite(Constraint): +class PositiveSemiDefinite(Constraint[E]): """ Constrain a matrix to be positive semidefinite. """ - expression: Expression + expression: E @dataclassabc(frozen=True) -class ExponentialCone(Constraint): - expression: Expression - - -# ┏━┓┏━┓┏━┓╻╺┳╸╻╻ ╻┏━┓╺┳╸┏━╸╻ ╻ ┏━╸┏┓╻┏━┓┏━┓┏━╸╺┳╸╺━┓┏━╸ -# ┣━┛┃ ┃┗━┓┃ ┃ ┃┃┏┛┗━┓ ┃ ┣╸ ┃ ┃ ┣╸ ┃┗┫┗━┓┣━┫┣╸ ┃ ┏━┛┣╸ -# ╹ ┗━┛┗━┛╹ ╹ ╹┗┛ ┗━┛ ╹ ┗━╸┗━╸┗━╸┗━╸╹ ╹┗━┛╹ ╹┗━╸ ╹ ┗━╸┗━╸ - -def putinar_psatz( - constr: NonNegative, - variables: Iterable[OptVariable], - polynomial_variables: Iterable[Variable], - multiplier_name: str - ) -> Iterable[Constraint]: - """ - Apply Putinar's positivstellensatz to a non-negativity constraint over an - archimedean domain. - """ - if not constr.domain: - # TODO: maybe raise an exception to warn the user? - return (constr,) - - # Number of variables in the original problem - q = sum(math.prod(v.shape) for v in variables) - x = poly.v_stack((1,) + tuple(polynomial_variables)) - - new_constraints = [] - new_constr = constr.expression - - for domain_poly in constr.domain.polynomials: - # To know the degree of the domain polynomial we need to evaluate the - # expression, hence we use a monad here - def make_multiplier(state: ExpressionState) -> tuple[ExpressionState, PolyMatrix]: - state, constr_pm = constr.expression.degree().apply(state) - state, domain_poly_pm = domain_poly.degree().apply(state) - - constr_deg = constr_pm.at(0, 0).constant() - domain_poly_deg = domain_poly_pm.at(0, 0).constant() - - # degree of multiplier - d = constr_deg - domain_poly_deg - if d < 0: - raise AlgebraicError("Cannot create Positivstellensatz multiplier because " - f"constraint polynomial has degree {constr_deg} and domain " - f"polynomial {domain_poly} has degree {domain_poly_deg}") - - n = sum(math.comb(k + q - 1, k) for k in range(0, d+1)) - - # TODO: check that there is not a variable with this name already - # Coefficients of the multiplier polynomial - c = opt_variable_from_name(multiplier_name, shape=(1,n)) - multiplier = c @ x.combinations(tuple(range(d +1))) - - return multiplier.apply(state) - - m = poly.from_statemonad(init_state_monad(make_multiplier)) - - # Subtract multiplier from expression and impose that it is also SOS - new_constr -= m * domain_poly - new_constraints.append(NonNegative(m)) - - new_constraints.append(NonNegative(new_constr)) - return new_constraints +class ExponentialCone(Constraint[E]): + expression: E diff --git a/sumofsquares/expression.py b/sumofsquares/expression.py index 08b4fbd..28f1794 100644 --- a/sumofsquares/expression.py +++ b/sumofsquares/expression.py @@ -6,18 +6,30 @@ from typing_extensions import override from dataclassabc import dataclassabc from polymatrix.expression.expression import ExpressionBaseMixin -from polymatrix.expressionstate.abc import ExpressionState +from polymatrix.expression.from_ import from_statemonad +from polymatrix.expressionstate import ExpressionState +from polymatrix.statemonad import init_state_monad +from polymatrix.polymatrix.mixins import PolyMatrixMixin + +from polymatrix.expression.init import ( + init_concatenate_expr, + init_diag_expr, + init_lower_triangular_expr, + init_transpose_expr, + init_slice_expr) from .abc import Constraint -from .constraints import NonNegative +from .constraints import PositiveSemidefinite, ExponentialCone +from .error import SolverError +from .variable import init_opt_variable_expr class SOSExpressionBaseMixin(ExpressionBaseMixin): @override def apply(self, state: ExpressionState) -> tuple[ExpressionState, ExpressionBaseMixin]: - # FIXME: think of name for this exception - raise RuntimeError - + raise SolverError(f"Expression containing {self.__class__.__qualname__} " + "cannot be used directly, they need to be rewritten " + "into an equivalent form using a canonicalization function.") @abstractmethod def recast(self) -> tuple[ExpressionBaseMixin, Iterable[Constraint]]: @@ -42,17 +54,79 @@ class LogDetMixin(SOSExpressionBaseMixin): def recast(self) -> tuple[ExpressionBaseMixin, Iterable[Constraint]]: # The problem # - # minimize logdet(A) + # maximize logdet(A) # # is equivalent to solving # - # minimize -t + # maximize t + # + # subject to [ A Z ] + # [ Z.T diag(Z) ] >= 0 + # + # Z lower triangular + # t <= sum_i log(Z[i,i]) + # + # and the last constraint of the above is equivalent to + # + # t <= sum_i u[i] + # u_i <= log(Z[i, i]) for all i + # + # And finally to get rid of the log the latter constraint one is + # equivalent to + # + # (Z[i,i], 1, u[i]) in Exponential Cone for all i + # + # Hence we can replace the original problem with + # + # minimize - sum_i u[i] # # subject to [ A Z ] # [ Z.T diag(Z) ] >= 0 # # Z lower triangular - # t <= sum_i log(Z[i, i]) + # (Z[i,i], 1, u[i]) in ExpCone for all i + + A = self.underlying + + # FIXME: get rid of these functions, create ShapeExprMixin in polymatrix? + def make_u(state: ExpressionState) -> tuple[ExpressionState, PolyMatrixMixin]: + state, pm = A.apply(state) + n, m = pm.shape + if n != m: + raise ValueError(f"Matrix A of logdet(A) must be square, " + f"but it has shape {pm.shape}") + + # FIXME: should check to avoid name clashes + u = init_opt_variable_expr("u_logdet", shape=(n, 1)) + return u.apply(state) + + def make_z(state: ExpressionState) -> tuple[ExpressionState, PolyMatrixMixin]: + state, pm = A.apply(state) + n, m = pm.shape + if n != m: + raise ValueError(f"Matrix A of logdet(A) must be square, " + f"but it has shape {pm.shape}") + + # FIXME: should check to avoid name clashes + Z = init_lower_triangular_expr(init_opt_variable_expr("Z_logdet", shape=(n * (n + 1) // 2, 1))) + return Z.apply(state) + + Z = from_statemonad(init_state_monad(make_z)) + Z_T = init_transpose_expr(Z) + Z_diag = init_diag_expr(Z) + + u = from_statemonad(init_state_monad(make_u)) + + # we call the new big matrix Q + Q = init_concatenate_expr(((A, Z), (Z_T, Z_diag))) + + def make_expcones(state: ExpressionState) -> tuple[ExpressionState, tuple[Constraint]]: + pass + + + constraints = [PositiveSemidefinite(Q)] + + raise NotImplementedError @dataclassabc(froze=True) diff --git a/sumofsquares/problems.py b/sumofsquares/problems.py index 500d06a..6161641 100644 --- a/sumofsquares/problems.py +++ b/sumofsquares/problems.py @@ -1,31 +1,91 @@ """ Specific optimization problems that are related to Sum of Squares. + +TODO: explanation for InternalSOSProblem """ from __future__ import annotations -import math import polymatrix as poly +import numpy as np from dataclassabc import dataclassabc from dataclasses import replace -from itertools import chain +from numpy.typing import NDArray from typing import Any, Sequence from typing_extensions import override from polymatrix.expression.expression import Expression, VariableExpression, init_expression from polymatrix.expressionstate import ExpressionState +from polymatrix.polymatrix.mixins import PolyMatrixMixin +from polymatrix.polymatrix.index import MonomialIndex, VariableIndex from polymatrix.variable import Variable from .abc import Problem, Constraint, Solver, Result -from .constraints import NonNegative, putinar_psatz -from .error import AlgebraicError +from .constraints import NonNegative, EqualToZero, PositiveSemiDefinite, ExponentialCone 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 +from .variable import OptVariable + + +# ┏━╸┏━┓┏┓╻╻┏━╸ ┏━┓┏━┓┏━┓┏┓ ╻ ┏━╸┏┳┓ +# ┃ ┃ ┃┃┗┫┃┃ ┣━┛┣┳┛┃ ┃┣┻┓┃ ┣╸ ┃┃┃ +# ┗━╸┗━┛╹ ╹╹┗━╸ ╹ ╹┗╸┗━┛┗━┛┗━╸┗━╸╹ ╹ + + +@dataclassabc(frozen=True) +class ConicProblem(Problem): + """ + Conic program with linear or quadratic cost functions. + All values stored in this class are numerical (numpy arrays). + + The pproblem has the following form + :: + + minimize .5 * x.T @ P @ x + q.T @ x + + subject to G @ x <= h + A @ x = b + + """ + P: NDArray + q: NDArray + A: NDArray + b: NDArray + G: NDArray + h: NDArray + + dims: dict + is_qp: bool + solver: Solver + variables: Sequence[OptVariable] + + @property + @override + def cost(self) -> tuple[NDArray, NDArray]: + return self.P, self.q + + @property + @override + def constraints(self) -> dict[str, tuple[NDArray, NDArray]]: + return { + "eq": (self.A, self.b), + "ineq": (self.G, self.h) + } + + @override + def solve(self, verbose: bool = False) -> ConicResult: + match self.solver: + case Solver.CVXOPT: + result, info = cvxopt_solve_sos_cone(self, verbose) + case _: + raise NotImplementedError + + return ConicResult(values=result, solver_info=info) + @dataclassabc(frozen=True) -class SOSResult(Result): +class ConicResult(Result): values: dict[OptVariable, float] solver_info: Any @@ -46,75 +106,234 @@ class SOSResult(Result): return self.values[var] +# ┏━┓╻ ╻┏┳┓ ┏━┓┏━╸ ┏━┓┏━┓╻ ╻┏━┓┏━┓┏━╸┏━┓ ┏━┓┏━┓┏━┓┏━╸┏━┓┏━┓┏┳┓ +# ┗━┓┃ ┃┃┃┃ ┃ ┃┣╸ ┗━┓┃┓┃┃ ┃┣━┫┣┳┛┣╸ ┗━┓ ┣━┛┣┳┛┃ ┃┃╺┓┣┳┛┣━┫┃┃┃ +# ┗━┛┗━┛╹ ╹ ┗━┛╹ ┗━┛┗┻┛┗━┛╹ ╹╹┗╸┗━╸┗━┛ ╹ ╹┗╸┗━┛┗━┛╹┗╸╹ ╹╹ ╹ + + @dataclassabc(frozen=True) -class SchmuedgenSOSProblem(Problem): - r""" - SOS problem with quadratic cost, linear constraints and optinally with - positivity constraints over a basic closed semialgebraic set - :math:`\mathbf{K}`. To constrain the positivity over :math:`\mathbf{K}` - this problem uses Schmüdgen's Positivstellensatz. +class SOSProblem(Problem): + """ + Generic sum of squares problem. + This problem contains expression objects. """ cost: Expression - constraints: Sequence[Constraint] - variables: Sequence[OptVariable] - polynomial_variables: Sequence[Variable] + constraints: Sequence[Constraint[Expression]] solver: Solver - state: ExpressionState + # --- dunder methods --- + + def __str__(self): + n = len(self.constraints) + s = f"SumOfSquares Problem ({self.__class__.__qualname__}):\n\n " \ + f"minimize {self.cost}\n " \ + f"subject to ({n} constraint{'s' if n > 1 else ''})\n" + + for c in self.constraints: + s += f"\t* {str(c)}\n" + + return s + + @override + def solve(self, verbose: bool = False, state: ExpressionState | None = None) -> Result: + if state is None: + state = poly.make_state() + + state, internal_prob = self.apply(state) + return internal_prob.solve(verbose) + + def apply(self, state: ExpressionState) -> tuple[ExpressionState, InternalSOSProblem]: + """ Convert to internal SOS problem by applying state to the expressions. """ + + constraints: list[Constraint[PolyMatrixMixin]] = [] + state, cost = self.cost.apply(state) + + # Compute the polymatrix of each constraint expression. Even though the + # result may not be used, this is necessary to account for the case of + # variables that are only present in the constraints, which the state + # object may not contain yet. + # TODO: Is there a more efficient way? There has to be a better way. + for c in self.constraints: + state, _ = c.expression.apply(state) + + # Collect variables + def is_optvariable(v): + return isinstance(v, OptVariable) + + polynomial_variables, variables = partition(is_optvariable, state.indices.keys()) + polynomial_variables = tuple(polynomial_variables) # because it is a generator + x = poly.v_stack(polynomial_variables) + for c in self.constraints: + # Scalar polynomial equality and non-negativity contraints need to + # be converted to (linear) coefficient matching constraints + if isinstance(c, NonNegative | EqualToZero): + state, pm = c.expression.quadratic_in(x).apply(state) + + elif isinstance(c, PositiveSemiDefinite): + state, pm = c.expression.apply(state) + + elif isinstance(c, ExponentialCone): + raise NotImplementedError + + else: + raise NotImplementedError + + constraints.append(replace(c, expression=pm)) + + return state, InternalSOSProblem(cost, tuple(constraints), + tuple(variables), polynomial_variables, + self.solver, state) + + +# FIXME: I hate this name, and the fact that this class needs to exist @dataclassabc(frozen=True) -class PutinarSOSProblem(Problem): - r""" - SOS problem with quadratic cost, linear constraints and optinally with - positivity constraints over an archimedean set :math:`\mathbf{K}`. To - constrain the positivity over :math:`\mathbf{K}` this problem uses - Putinar's positivstellensatz. +class InternalSOSProblem(Problem): """ - cost: Expression - constraints: Sequence[Constraint] + This is a class for internal use only. It is a SOS problem that contains + polymatrix objects. + + It exists only because expressions are lazily evaluated, so if we want to + perform transformations on problems (see canon module), we may need to + apply the state to the expressions to convert the expressions to polymatrix + objects. + """ + cost: PolyMatrixMixin + constraints: Sequence[Constraint[PolyMatrixMixin]] variables: Sequence[OptVariable] polynomial_variables: Sequence[Variable] solver: Solver - state: ExpressionState - - def after_positivstellensatz(self) -> Problem: - """ - Apply Putinar's positivstellensatz to convert the SOS constraints - into coefficient matching constraints - """ - def need_positivstellensatz(c: Constraint) -> bool: - return isinstance(c, NonNegative) and (c.domain is not None) - - # Add langrange-multiplier-like polynomials to SOS constraints, while - # leaving the rest untouched. - constraints, to_process = partition(need_positivstellensatz, self.constraints) - new_constraints: list[Constraint] = [] - for i, c in enumerate(to_process): - new_constraints.extend(putinar_psatz(c, self.variables, self.polynomial_variables, r"\gamma_{i}")) - - # Convert to polymatrix object to extract variables created by p-satz - # FIXME: this is a dirty trick to make it work for now, it is not - # efficient, is there a better way? - new_state = self.state - for c in new_constraints: - new_state, _ = c.expression.apply(new_state) + # TODO: remove state field from this class, it is redundant + state: ExpressionState - def is_optvariable(v): - return isinstance(v, OptVariable) + def to_conic_problem(self) -> ConicProblem: + cost = poly.to_affine(self.cost) + if cost.degree > 2: + raise ValueError("This package can solve at most quadratic conic programs, " + f"but the given problem has degree {cost.degree} > 2.") + + is_qp = (cost.degree == 2) + + A_rows, b_rows = [], [] + G_rows, h_rows= [], [] + l_dims, q_dims, s_dims = [], [], [] + + # indices of variables in the optimization problem + variable_indices = sum(sorted(tuple(self.state.get_indices_as_variable_index(v)) + for v in self.variables), ()) + + # cost linear term + cost_coeffs = dict(cost.affine_coefficients_by_degrees(variable_indices)) + q = cost_coeffs.get(1) + if q is None and is_qp: + q = np.zeros((len(variable_indices), 1)) + + # cost quadratic term + if is_qp: + n = len(variable_indices) + 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) + + else: + P = None + + for c in self.constraints: + if isinstance(c, PositiveSemiDefinite): + constr = poly.to_affine(c.expression) - new_variables = filter(is_optvariable, new_state.indices.keys()) + nrows, ncols = constr.shape + if nrows != ncols: + raise ValueError(f"PSD constraint cannot contain non-square matrix of shape ({nrows, ncols})!") - return replace(self, - state=new_state, - variables=tuple(new_variables), - constraints=tuple(constraints) + tuple(new_constraints)) + s_dims.append(nrows) + + # constant term + c_coeff = constr.affine_coefficient(MonomialIndex.constant()) + c_stacked = c_coeff.T.reshape((nrows * ncols,)) + + # linear terms + l_stacked_rows = [] + for v in variable_indices: + v_coeff = constr.affine_coefficient(v) + v_stacked = v_coeff.T.reshape((nrows * ncols,)) + l_stacked_rows.append(v_stacked) + + l_stacked = np.vstack(l_stacked_rows) + + h_rows.append(c_stacked) + G_rows.append(-1 * l_stacked) + + elif isinstance(c, ExponentialCone): + # TODO: implement + raise NotImplementedError + + elif isinstance(c, EqualToZero | NonNegative): + constr = poly.to_affine(c.expression) + if constr.degree > 1: + raise ValueError("To convert to conic constraints must be linear or affine " + f"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: + # Get the affine coefficient associated to this variable + v_coeff = constr.affine_coefficient(v) + + # 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 conic constraint. " + "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 a conic constriant.") + + G = np.hstack(G_rows).T if G_rows else None + h = np.hstack(h_rows) if h_rows else None + A = np.hstack(A_rows).T if A_rows else None + b = np.hstack(b_rows) if b_rows else None + + dims = {'l': l_dims, 'q': q_dims, 's': s_dims} + + if all((G is None, h is None, A is None, b is None)): + raise ValueError("Optimization problem is unconstrained!") + + return ConicProblem(P=P, q=q, A=A, b=b, G=G, h=h, + dims=dims, is_qp=is_qp, + solver=self.solver, + variables=self.variables) @override def solve(self, verbose: bool = False) -> Result: - if not self.solver == Solver.CVXOPT: - raise NotImplementedError(f"Cannot solve problem {self} using solver {str(self.solver)}.") - - prob = self.after_positivstellensatz() - result, info = cvxopt_solve_sos_cone(prob, verbose) - return SOSResult(result, info) + return self.to_conic_problem().solve(verbose) diff --git a/sumofsquares/solver/cvxopt.py b/sumofsquares/solver/cvxopt.py index 24ec413..6418a2d 100644 --- a/sumofsquares/solver/cvxopt.py +++ b/sumofsquares/solver/cvxopt.py @@ -1,11 +1,14 @@ """ Solve sumofsquares problems using CVXOPT """ +from __future__ import annotations + import cvxopt import numpy as np import math from collections import UserDict +from typing import TYPE_CHECKING import polymatrix as poly from polymatrix.expression.expression import Expression @@ -16,39 +19,28 @@ from ..constraints import NonNegative, EqualToZero, PositiveSemiDefinite, Expone from ..error import NotSupportedBySolver, SolverError from ..variable import OptVariable +if TYPE_CHECKING: + from ..problems import ConicProblem + class CVXOPTInfo(SolverInfo, UserDict): """ Dictionary returned by CVXOPT. """ def __getattr__(self, attr): + # trick to make it behave like a dataclass 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, +def solve_sos_cone(prob: ConicProblem, 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)` using CVXOPT. - 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) + # CVXOPT can solve problems that have the form # # minimize .5 * x.T @ P @ x + q.T @ x # @@ -56,147 +48,28 @@ def solve_sos_cone(prob: Problem, verbose: bool = False, # 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 = [], [], [] - - # indices of variables in the optimization problem - variable_indices = sum(sorted(tuple(prob.state.get_indices_as_variable_index(v)) - for v in prob.variables), ()) - - # cost linear term - cost_coeffs = dict(cost.affine_coefficients_by_degrees(variable_indices)) - q = cost_coeffs.get(1) - if q is None and is_qp: - q = np.zeros((len(variable_indices), 1)) - - # cost quadratic term - if is_qp: - n = len(variable_indices) - 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): - constr = poly.to_affine(c.expression.read(prob.state)) - - nrows, ncols = constr.shape - if nrows != ncols: - raise ValueError(f"PSD constraint cannot contain non-square matrix of shape ({nrows, ncols})!") - - s_dims.append(nrows) - - # constant term - c_coeff = constr.affine_coefficient(MonomialIndex.constant()) - c_stacked = c_coeff.T.reshape((nrows * ncols,)) - - # linear terms - l_stacked_rows = [] - for v in variable_indices: - v_coeff = constr.affine_coefficient(v) - v_stacked = v_coeff.T.reshape((nrows * ncols,)) - l_stacked_rows.append(v_stacked) - - l_stacked = np.vstack(l_stacked_rows) - - h_rows.append(c_stacked) - G_rows.append(-1 * l_stacked) - - elif isinstance(c, ExponentialCone): - raise NotSupportedBySolver("CVXOpt cannot solve problems with the exponential cone.") - - elif isinstance(c, EqualToZero | NonNegative): - # Convert constraints to affine expression in decision variables - # TODO: Use half Newton polytope optimization to reduce size of marix problem - 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: - # Get the affine coefficient associated to this variable - v_coeff = constr.affine_coefficient(v) - - # 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.T) if q is not None else None - 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} - + # This matches the structure of ConicProblem, so the values can be used + # directly. try: - 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) + # Some are transposed because CVXOPT uses column major arrays + q = cvxopt.matrix(prob.q.T) + A = cvxopt.matrix(prob.A) if prob.A is not None else None + b = cvxopt.matrix(prob.b.T) if prob.b is not None else None + G = cvxopt.matrix(prob.G) if prob.G is not None else None + h = cvxopt.matrix(prob.h.T) if prob.h is not None else None + + if prob.is_qp: + P = cvxopt.matrix(prob.P) + + info = cvxopt.solvers.coneqp(P=P, q=q, G=G, h=h, A=A, b=b, + dims=prob.dims, *args, **kwargs) else: - info |= cvxopt.solvers.conelp(c=q, G=G, h=h, A=A, b=b, + info = cvxopt.solvers.conelp(c=q, G=G, h=h, A=A, b=b, *args, **kwargs) except Exception as e: - raise SolverError("CVXOpt can't solve this problem.") from e + raise SolverError("CVXOpt can't solve this problem, " + "see previous exception why.") from e results = {} |