diff options
| -rw-r--r-- | sumofsquares/error.py | 3 | ||||
| -rw-r--r-- | sumofsquares/problems.py | 105 |
2 files changed, 56 insertions, 52 deletions
diff --git a/sumofsquares/error.py b/sumofsquares/error.py index 6fe04fd..75bdfc4 100644 --- a/sumofsquares/error.py +++ b/sumofsquares/error.py @@ -13,3 +13,6 @@ class SolverError(Exception): class AlgebraicError(Exception): """ TODO """ + +class ConstraintWarning(UserWarning): + """ A constraint was removed from the problem. """ diff --git a/sumofsquares/problems.py b/sumofsquares/problems.py index 6737341..808e51c 100644 --- a/sumofsquares/problems.py +++ b/sumofsquares/problems.py @@ -10,10 +10,10 @@ import numpy as np from dataclassabc import dataclassabc from dataclasses import replace -from itertools import groupby from numpy.typing import NDArray from typing import Any, Sequence from typing_extensions import override +from warnings import warn from polymatrix.expression.expression import Expression, VariableExpression from polymatrix.expression.mixins.variableexprmixin import VariableExprMixin @@ -25,6 +25,7 @@ from polymatrix.symbol import Symbol from .abc import Problem, Constraint, Solver, Result from .constraints import NonNegative, EqualToZero, PositiveSemiDefinite, ExponentialCone +from .error import ConstraintWarning from .solver.cvxopt import solve_cone as cvxopt_solve_cone from .solver.scs import solve_cone as scs_solve_cone from .utils import partition @@ -207,7 +208,7 @@ class SOSProblem(Problem): for i, c in enumerate(self.constraints): if isinstance(c, EqualToZero): - state, deg = c.expression.degree().apply(state) + state, deg = c.expression.linear_in(x).degree().max().apply(state) # Polynomial equality must be converted into coefficient # matching condition @@ -226,28 +227,43 @@ class SOSProblem(Problem): "to conic constraint. Domain restricted non-negativity " "must be preprocessed by using a Positivstellensatz!") - state, deg = c.expression.degree().apply(state) + if isinstance(c.expression.shape, tuple): + nrows, ncols = c.expression.shape - if isinstance(c.expression.shape[0], int): - nrows = c.expression.shape[0] + # it is an expression else: - state, shape = c.expression.shape[0].apply(state) + state, shape = c.expression.shape.apply(state) nrows = shape.at(0, 0).constant() - - # It is a vector column goes in the non-negative orthant - if nrows > 1: - constraints.append(c) - - # Polynomial non-negativity contraints is converted to PSD - # constraint of SOS quadratic form - elif deg.scalar().constant() > 1: - # TODO: it seems to work fine even without .symmetric(). Why? - cnew = c.expression.quadratic_in(x).symmetric() - constraints.append(PositiveSemiDefinite(cnew)) - - # A normal (affine) constraint - else: - constraints.append(c) + ncols = shape.at(1, 0).constant() + + # Row wise interpretation, we need to know the shape + if ncols > 1: + raise ValueError(f"Cannot convert non-negativity constraint {str(c)} (nr {i}) " + "because it has more than one column. If you want a " + "matrix to be PSD use PositiveSemiDefinite, if you want " + "to have a vector in the non-negative orthant of R^n the " + "interpretation is row-wise, not column wise (transpose " + "your vector)") + + for row in range(nrows): + state, deg = c.expression.linear_monomials(x).degree().max().apply(state) + d = deg.scalar().constant() + + # Polynomial non-negativity contraints is converted to PSD + # constraint of SOS quadratic form + if d > 1: + if d % 2 != 0: + raise ValueError("Cannot convert expression inside non-negativity " + f"constraint {str(c)} (nr. {i}) into quadratic form for PSD condition " + f"because it has degree {deg.scalar().constant()}, which " + "is not even.") + + cnew = c.expression.quadratic_in(x).symmetric() + constraints.append(PositiveSemiDefinite(cnew)) + + # A normal (affine) constraint + else: + constraints.append(c) elif isinstance(c, PositiveSemiDefinite): state, pm = c.expression.cache().apply(state) @@ -273,28 +289,24 @@ class SOSProblem(Problem): # Convert Expressions into PolyMatrix objects pm_constraints: list[Constraint[PolyMatrixMixin]] = [] - for c in constraints: + for i, c in enumerate(constraints): state, pm = c.expression.apply(state) - pm_constraints.append(replace(c, expression=pm)) - - # TODO: can we get rid of for loop inside InternalSOSProblem.to_conic_problem? - # Concatenate constraints so that there is only a big constraint per cone. - # for (ctype, group) in groupby(constraints, key=type): - # if ctype in (EqualToZero, NonNegative): - # state, pm = poly.v_stack((c.expression for c in group)).apply(state) - # pm_constraints.append(ctype(pm)) - - # elif ctype is PositiveSemiDefinite: - # expressions = (c.expression for c in group) - # state, pm = poly.block_diag(expressions).apply(state) - # pm_constraints.append(ctype(pm)) + constr = poly.to_affine(pm) - # elif ctype is ExponentialCone: - # state, pm = poly.v_stack((c.expression for c in group)).apply(state) - # pm_constraints.append(ctype(pm)) + # Remove constraints that have disappeared + # There is only one "variable" in the affine expression + if len(constr.slices.keys()) == 1: + # That variable is actually not a variable + if MonomialIndex.constant() in constr.slices.keys(): + # constraint has disappeared + warn(f"Constraint on {str(c.expression)} was removed.", ConstraintWarning) + continue + # FIXME: need to check for each type of constraint + # that it is ok to remove constant value. Example of + # a problematic edge case is NonNegative(a - a - 5) + # where a is an optimization variable. - # else: - # raise NotImplementedError(f"Cannot process constraint of type {ctype} (yet).") + pm_constraints.append(replace(c, expression=constr)) return state, InternalSOSProblem(cost, tuple(pm_constraints), symbols, polynomial_symbols, @@ -383,20 +395,9 @@ class InternalSOSProblem(Problem): P = None for c in self.constraints: - constr = poly.to_affine(c.expression) + constr = c.expression nrows, ncols = constr.shape - # There is only one "variable" in the affine expression - if len(constr.slices.keys()) == 1: - # That variable is actually not a variable - if MonomialIndex.constant() in constr.slices.keys(): - # constraint has disappeared - continue - # FIXME: need to check for each type of constraint - # that it is ok to remove constant value. Example of - # a problematic edge case is NonNegative(a - a - 5) - # where a is an optimization variable. - if constr.degree > 1: # If this error occurs an it is not the user's fault, there is # a bug in SOSProblem.apply |