diff options
| -rw-r--r-- | sumofsquares/canon.py | 89 |
1 files changed, 28 insertions, 61 deletions
diff --git a/sumofsquares/canon.py b/sumofsquares/canon.py index 8186f59..bec9cf0 100644 --- a/sumofsquares/canon.py +++ b/sumofsquares/canon.py @@ -4,6 +4,7 @@ import polymatrix as poly from abc import abstractmethod from dataclassabc import dataclassabc from dataclasses import replace +from typing import Iterable from typing_extensions import override from polymatrix.expression.from_ import Expression @@ -22,40 +23,18 @@ class Canonicalization(Problem): """ Transform a problem into another problem. """ + def __post_init__(self): + """ TODO: explain """ + if self.cost is None: + object.__setattr__(self, "cost", self.problem.cost) + + if self.constraints is None: + object.__setattr__(self, "constraints", self.problem.constraints) def __str__(self): s = f"Apply canonicalization procedure {self.__class__.__qualname__} to\n" return s + str(self.problem) - def __call__(self, problem : SOSProblem) -> SOSProblem: - # TODO: is this a good idea? probably not but I don't want to introduce - # a new type of class for "horizontal" function in the commutative - # diagram. The idea is to make canon objects behave like functions - # so that given a problem P, you can do Canon2(Canon1(P)), but this - # needs a default implementation if the transformation is a stateful - # computation, which would need a statemonad, defeating the purpose of - # creating this class in the first place. - """ - Apply the canonicalization procedure to the given problem. - - Override this method to write transformation if it can be written - without having to apply the state. I.e. it is not stateful in the - problem. - """ - raise NotImplementedError - - @property - @override - def cost(self) -> Expression: - # TODO: error message - raise AttributeError - - @property - @override - def constraints(self) -> Expression: - # TODO: error message - raise AttributeError - @property @override def solver(self) -> Solver: @@ -90,12 +69,8 @@ class PutinarPSatz(Canonicalization): """ problem: SOSProblem - - @property - @override - def cost(self) -> Expression: - # PSatz does not change cost function - return self.problem.cost + cost: Expression | None = None + constraints: list[Expression] | None = None @override def apply(self, state: ExpressionState) -> tuple[ExpressionState, InternalSOSProblem]: @@ -155,7 +130,8 @@ class PutinarPSatz(Canonicalization): monomials = x.combinations(poly.arange(d + 1)).half_newton_polytope(x) coeffs = opt_variable_from_name(rf"\gamma_{i},c", shape=monomials.shape) - new_constr = constr.expression - coeffs.T @ monomials + multiplier = poly.give_name(coeffs.T @ monomials, rf"\gamma_{i},c") + new_constr = constr.expression - multiplier # Add other terms of the quadratic module for j, domain_poly in enumerate(constr.domain.polynomials): @@ -189,7 +165,7 @@ class PutinarPSatz(Canonicalization): monomials = x.combinations(poly.arange(d + 1)).half_newton_polytope(x) coeffs = opt_variable_from_name(rf"\gamma_{i},{j}", shape=monomials.shape) - multiplier = coeffs.T @ monomials + multiplier = poly.give_name(coeffs.T @ monomials, rf"\gamma_{i},{j}") # Subtract multiplier from expression and impose that it is also SOS new_constr -= multiplier * domain_poly @@ -204,12 +180,8 @@ class LogDet(Canonicalization): # TODO: reconsider extension of polymatrix Expression objects to introduce # logdet expression. The current solution works but is not ideal, as it may # be a bit counterintuitive. - - # Does this trick work with the det(A) instead of logdet? - # No, because then the inequality on the diagonal of Z is not a sum of log - # but a product of variables, i.e. non linear. """ - Create a new problem, whose cost function is the logdet of the given + Create a new problem, whose cost function is the -logdet of the given problem's cost function. The problem @@ -247,25 +219,19 @@ class LogDet(Canonicalization): (Z[i,i], 1, u[i]) in ExpCone for all i """ problem: SOSProblem + cost: Expression | None = None + constraints: list[Expression] | None = None - @property @override - def cost(self) -> Expression: - return self(self.problem).cost - - @property - @override - def constraints(self) -> Expression: - return self(self.problem).constraints - - @override - def __call__(self, problem : SOSProblem) -> SOSProblem: - A = problem.cost + def apply(self, state: ExpressionState) -> tuple[ExpressionState, InternalSOSProblem]: + A = self.problem.cost # FIXME: should check for name clashes in state object # i.e. what if there is already a variable named u_logdet / Z_logdet? # TODO: deterministically generate unique names + t = opt_variable_from_name('t_logdet') + n = A.shape[0] u = opt_variable_from_name('u_logdet', shape=(n, 1)) @@ -275,11 +241,12 @@ class LogDet(Canonicalization): Q = poly.concatenate(((A, Z), (Z.T, Z.diag()))) E = poly.h_stack((Z.diag(), poly.ones((n, 1)), u)) - new_cost = - u.T.sum() - new_constraints = tuple(problem.constraints) + (PositiveSemiDefinite(Q), ExponentialCone(E)) + new_cost = - t + new_constraints = tuple(self.problem.constraints) + ( + PositiveSemiDefinite(Q), + ExponentialCone(E), + NonNegative(u), + NonNegative(u.T.sum() - t) + ) - return replace(problem, cost=new_cost, constraints=new_constraints) - - @override - def apply(self, state: ExpressionState) -> tuple[ExpressionState, InternalSOSProblem]: - return self(self.problem).apply(state) + return replace(self.problem, cost=new_cost, constraints=new_constraints).apply(state) |