Coverage for picos/constraints/con_quadratic.py: 92.09%
139 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
1# ------------------------------------------------------------------------------
2# Copyright (C) 2019 Maximilian Stahlberg
3#
4# This file is part of PICOS.
5#
6# PICOS is free software: you can redistribute it and/or modify it under the
7# terms of the GNU General Public License as published by the Free Software
8# Foundation, either version 3 of the License, or (at your option) any later
9# version.
10#
11# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY
12# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
14#
15# You should have received a copy of the GNU General Public License along with
16# this program. If not, see <http://www.gnu.org/licenses/>.
17# ------------------------------------------------------------------------------
19"""Quadratic constraint types."""
21from collections import namedtuple
23from .. import glyphs
24from ..apidoc import api_end, api_start
25from ..caching import cached_property
26from .constraint import Constraint, ConstraintConversion
28_API_START = api_start(globals())
29# -------------------------------
32class NonconvexQuadraticConstraint(Constraint):
33 """Bound on a nonconvex quadratic expression."""
35 def __init__(self, lhs, relation, rhs):
36 """Construct a :class:`NonconvexQuadraticConstraint`.
38 :param lhs: Left hand side quadratic or affine expression.
39 :type lhs: ~picos.expressions.QuadraticExpression or
40 ~picos.expressions.AffineExpression
41 :param str relation: Constraint relation symbol.
42 :param rhs: Right hand side quadratic or affine expression.
43 :type rhs: ~picos.expressions.QuadraticExpression or
44 ~picos.expressions.AffineExpression
45 """
46 from ..expressions import AffineExpression, QuadraticExpression
48 assert isinstance(lhs, (AffineExpression, QuadraticExpression))
49 assert isinstance(rhs, (AffineExpression, QuadraticExpression))
50 assert any(isinstance(exp, QuadraticExpression) for exp in (lhs, rhs))
51 assert relation in self.LE + self.GE
52 assert lhs.size == rhs.size
54 self.lhs = lhs
55 self.rhs = rhs
56 self.relation = relation
58 Constraint.__init__(self, self._get_type_term())
60 def _get_type_term(self):
61 return "Nonconvex Quadratic"
63 # TODO: Add common interface for all the le0, ge0, smaller, greater methods?
64 @cached_property
65 def le0(self):
66 """Quadratic expression constrained to be at most zero."""
67 if self.relation == self.GE:
68 return self.rhs - self.lhs
69 else:
70 return self.lhs - self.rhs
72 @property
73 def smaller(self):
74 """Smaller-or-equal side of the constraint."""
75 return self.rhs if self.relation == self.GE else self.lhs
77 @property
78 def greater(self):
79 """Greater-or-equal side of the constraint."""
80 return self.lhs if self.relation == self.GE else self.rhs
82 Subtype = namedtuple("Subtype", ("qterms"))
84 def _subtype(self):
85 return self.Subtype(self.le0.num_quad_terms)
87 @classmethod
88 def _cost(cls, subtype):
89 return subtype.qterms
91 def _expression_names(self):
92 yield "lhs"
93 yield "rhs"
95 def _str(self):
96 if self.relation == self.LE:
97 return glyphs.le(self.lhs.string, self.rhs.string)
98 else:
99 return glyphs.ge(self.lhs.string, self.rhs.string)
101 def _get_size(self):
102 return (1, 1)
104 def _get_slack(self):
105 if self.relation == self.LE:
106 return self.rhs.safe_value - self.lhs.safe_value
107 else:
108 return self.lhs.safe_value - self.rhs.safe_value
111class ConicQuadraticConstraint(NonconvexQuadraticConstraint):
112 r"""Bound on a *nearly* convex quadratic expression.
114 Nearly convex means that the bilinear form representing the quadratic part
115 of the expression that the constraint poses to be at most zero has exactly
116 one negative eigenvalue. More precisely, if the constraint is of the form
117 :math:`x^TQx + a^Tx + b \leq x^TRx + b^Tx + c`, then `Q - R` needs to have
118 exactly one negative eigenvalue for the constraint to be nearly convex. Such
119 constraints can be posed as conic constraints under an additional assumption
120 that some affine term is nonnegative.
122 At this point, this class may only be used for nonconvex constraints of the
123 form :math:`x^TQx + p^Tx + q \leq (a^Tx + b)(c^Tx + d)` with
124 :math:`x^TQx + p^Tx + q` representable as a squared norm. In this case, the
125 additional assumptions required for a conic reformulation are
126 :math:`a^Tx + b \geq 0` and :math:`b^Tx + c \geq 0`.
128 Whether a constraint of this type is strengthened to a conic constraint or
129 relabeled as a :class:`NonconvexQuadraticConstraint` depends on the
130 :ref:`assume_conic <option_assume_conic>` option.
132 :Example:
134 >>> from picos import Options, RealVariable
135 >>> x, y, z = RealVariable("x"), RealVariable("y"), RealVariable("z")
136 >>> C = x**2 + 1 <= y*z; C
137 <Conic Quadratic Constraint: x² + 1 ≤ y·z>
138 >>> P = C.__class__.Conversion.convert(C, Options(assume_conic=True))
139 >>> list(P.constraints.values())[0]
140 <4×1 RSOC Constraint: ‖fullroot(x² + 1)‖² ≤ y·z ∧ y, z ≥ 0>
141 >>> Q = C.__class__.Conversion.convert(C, Options(assume_conic=False))
142 >>> list(Q.constraints.values())[0]
143 <Nonconvex Quadratic Constraint: x² + 1 ≤ y·z>
145 .. note::
147 Solver implementations must not support this constraint type so that
148 the user's choice for the :ref:`assume_conic <option_assume_conic>`
149 option is respected.
150 """
152 class Conversion(ConstraintConversion):
153 """Nearly convex quadratic to (rotated) second order cone conversion."""
155 @classmethod
156 def predict(cls, subtype, options):
157 """Implement :meth:`~.constraint.ConstraintConversion.predict`."""
158 from . import AbsoluteValueConstraint, SOCConstraint, RSOCConstraint
160 if options.assume_conic:
161 n = subtype.conic_argdim
163 if subtype.rotated:
164 yield ("con", RSOCConstraint.make_type(argdim=n), 1)
165 elif n == 1:
166 yield ("con", AbsoluteValueConstraint.make_type(), 1)
167 else:
168 yield ("con", SOCConstraint.make_type(argdim=n), 1)
169 else:
170 yield ("con", NonconvexQuadraticConstraint.make_type(
171 qterms=subtype.qterms), 1)
173 @classmethod
174 def convert(cls, con, options):
175 """Implement :meth:`~.constraint.ConstraintConversion.convert`."""
176 from ..expressions import AffineExpression
177 from ..modeling import Problem
178 from . import AbsoluteValueConstraint, SOCConstraint, RSOCConstraint
180 P = Problem()
182 if options.assume_conic:
183 q, r = con.smaller, con.greater
184 root = AffineExpression.zero() if q.is0 else q.fullroot
186 if r.scalar_factors[0] is not r.scalar_factors[1]:
187 P.add_constraint(RSOCConstraint(root, *r.scalar_factors))
188 else:
189 bound = r.scalar_factors[0]
191 if len(root) == 1:
192 P.add_constraint(AbsoluteValueConstraint(root, bound))
193 else:
194 P.add_constraint(SOCConstraint(root, bound))
195 else:
196 P.add_constraint(NonconvexQuadraticConstraint(
197 con.lhs, con.relation, con.rhs))
199 return P
201 def __init__(self, lhs, relation, rhs):
202 """Construct a :class:`ConicQuadraticConstraint`.
204 See :meth:`NonconvexQuadraticConstraint.__init__` for more.
205 """
206 from ..expressions import QuadraticExpression
208 NonconvexQuadraticConstraint.__init__(self, lhs, relation, rhs)
210 # Additional usage restrictions for the near-convex case.
211 assert isinstance(self.lhs, QuadraticExpression)
212 assert isinstance(self.rhs, QuadraticExpression)
213 assert self.smaller.is_squared_norm
214 assert self.greater.scalar_factors
216 def _get_type_term(self):
217 return "Conic Quadratic"
219 Subtype = namedtuple("Subtype", ("qterms", "conic_argdim", "rotated"))
221 def _subtype(self):
222 qterms = self.le0.num_quad_terms
223 conic_argdim = 1 if self.smaller.is0 else len(self.smaller.fullroot)
224 sf = self.greater.scalar_factors
225 rotated = sf[0] is not sf[1]
227 return self.Subtype(qterms, conic_argdim, rotated)
230class ConvexQuadraticConstraint(NonconvexQuadraticConstraint):
231 """Bound on a convex quadratic expression."""
233 class ConicConversion(ConstraintConversion):
234 """Convex quadratic to (rotated) second order cone conversion."""
236 @classmethod
237 def predict(cls, subtype, options):
238 """Implement :meth:`~.constraint.ConstraintConversion.predict`."""
239 from . import AbsoluteValueConstraint, SOCConstraint, RSOCConstraint
241 if subtype.haslin:
242 yield ("con", RSOCConstraint.make_type(argdim=subtype.rank), 1)
243 elif subtype.rank == 1:
244 yield ("con", AbsoluteValueConstraint.make_type(), 1)
245 else:
246 yield ("con", SOCConstraint.make_type(argdim=subtype.rank), 1)
248 @classmethod
249 def convert(cls, con, options):
250 """Implement :meth:`~.constraint.ConstraintConversion.convert`."""
251 from ..expressions import AffineExpression
252 from ..modeling import Problem
253 from . import AbsoluteValueConstraint, SOCConstraint, RSOCConstraint
255 le0 = con.le0
257 P = Problem()
259 if con.subtype.haslin:
260 one = AffineExpression.from_constant(1)
262 P.add_constraint(RSOCConstraint(le0.quadroot, -le0._aff, one))
263 else:
264 assert le0._aff.constant
265 value = -le0._aff.value
267 if value < 0:
268 raise ValueError("The constraint {} is infeasible as it "
269 "upper-bounds a positive semidefinite quadratic form by"
270 " a negative constant.".format(con))
272 root = le0.quadroot
273 bound = AffineExpression.from_constant(value**0.5,
274 name=glyphs.sqrt(glyphs.scalar(value)))
276 if len(root) == 1:
277 P.add_constraint(AbsoluteValueConstraint(root, bound))
278 else:
279 P.add_constraint(SOCConstraint(root, bound))
281 return P
283 def __init__(self, lhs, relation, rhs):
284 """Construct a :class:`ConvexQuadraticConstraint`.
286 See :meth:`NonconvexQuadraticConstraint.__init__` for more.
287 """
288 NonconvexQuadraticConstraint.__init__(self, lhs, relation, rhs)
290 def _get_type_term(self):
291 return "Convex Quadratic"
293 Subtype = namedtuple("Subtype", ("qterms", "rank", "haslin"))
295 def _subtype(self):
296 from ..expressions import QuadraticExpression
298 # HACK: Be consistent with the prediction, which cannot foresee whether
299 # the linear terms of both sides cancel out.
300 lhs_aff_is_cst = self.lhs._aff.constant \
301 if isinstance(self.lhs, QuadraticExpression) else self.lhs.constant
302 rhs_aff_is_cst = self.rhs._aff.constant \
303 if isinstance(self.rhs, QuadraticExpression) else self.rhs.constant
304 haslin = not lhs_aff_is_cst or not rhs_aff_is_cst
306 return self.Subtype(self.le0.num_quad_terms, self.le0.rank, haslin)
309# --------------------------------------
310__all__ = api_end(_API_START, globals())