Coverage for picos/constraints/con_quadratic.py: 92.09%

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# ------------------------------------------------------------------------------ 

18 

19"""Quadratic constraint types.""" 

20 

21from collections import namedtuple 

22 

23from .. import glyphs 

24from ..apidoc import api_end, api_start 

25from ..caching import cached_property 

26from .constraint import Constraint, ConstraintConversion 

27 

28_API_START = api_start(globals()) 

29# ------------------------------- 

30 

31 

32class NonconvexQuadraticConstraint(Constraint): 

33 """Bound on a nonconvex quadratic expression.""" 

34 

35 def __init__(self, lhs, relation, rhs): 

36 """Construct a :class:`NonconvexQuadraticConstraint`. 

37 

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 

47 

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 

53 

54 self.lhs = lhs 

55 self.rhs = rhs 

56 self.relation = relation 

57 

58 Constraint.__init__(self, self._get_type_term()) 

59 

60 def _get_type_term(self): 

61 return "Nonconvex Quadratic" 

62 

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 

71 

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 

76 

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 

81 

82 Subtype = namedtuple("Subtype", ("qterms")) 

83 

84 def _subtype(self): 

85 return self.Subtype(self.le0.num_quad_terms) 

86 

87 @classmethod 

88 def _cost(cls, subtype): 

89 return subtype.qterms 

90 

91 def _expression_names(self): 

92 yield "lhs" 

93 yield "rhs" 

94 

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) 

100 

101 def _get_size(self): 

102 return (1, 1) 

103 

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 

109 

110 

111class ConicQuadraticConstraint(NonconvexQuadraticConstraint): 

112 r"""Bound on a *nearly* convex quadratic expression. 

113 

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. 

121 

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`. 

127 

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. 

131 

132 :Example: 

133 

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> 

144 

145 .. note:: 

146 

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 """ 

151 

152 class Conversion(ConstraintConversion): 

153 """Nearly convex quadratic to (rotated) second order cone conversion.""" 

154 

155 @classmethod 

156 def predict(cls, subtype, options): 

157 """Implement :meth:`~.constraint.ConstraintConversion.predict`.""" 

158 from . import AbsoluteValueConstraint, SOCConstraint, RSOCConstraint 

159 

160 if options.assume_conic: 

161 n = subtype.conic_argdim 

162 

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) 

172 

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 

179 

180 P = Problem() 

181 

182 if options.assume_conic: 

183 q, r = con.smaller, con.greater 

184 root = AffineExpression.zero() if q.is0 else q.fullroot 

185 

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] 

190 

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)) 

198 

199 return P 

200 

201 def __init__(self, lhs, relation, rhs): 

202 """Construct a :class:`ConicQuadraticConstraint`. 

203 

204 See :meth:`NonconvexQuadraticConstraint.__init__` for more. 

205 """ 

206 from ..expressions import QuadraticExpression 

207 

208 NonconvexQuadraticConstraint.__init__(self, lhs, relation, rhs) 

209 

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 

215 

216 def _get_type_term(self): 

217 return "Conic Quadratic" 

218 

219 Subtype = namedtuple("Subtype", ("qterms", "conic_argdim", "rotated")) 

220 

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] 

226 

227 return self.Subtype(qterms, conic_argdim, rotated) 

228 

229 

230class ConvexQuadraticConstraint(NonconvexQuadraticConstraint): 

231 """Bound on a convex quadratic expression.""" 

232 

233 class ConicConversion(ConstraintConversion): 

234 """Convex quadratic to (rotated) second order cone conversion.""" 

235 

236 @classmethod 

237 def predict(cls, subtype, options): 

238 """Implement :meth:`~.constraint.ConstraintConversion.predict`.""" 

239 from . import AbsoluteValueConstraint, SOCConstraint, RSOCConstraint 

240 

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) 

247 

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 

254 

255 le0 = con.le0 

256 

257 P = Problem() 

258 

259 if con.subtype.haslin: 

260 one = AffineExpression.from_constant(1) 

261 

262 P.add_constraint(RSOCConstraint(le0.quadroot, -le0._aff, one)) 

263 else: 

264 assert le0._aff.constant 

265 value = -le0._aff.value 

266 

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)) 

271 

272 root = le0.quadroot 

273 bound = AffineExpression.from_constant(value**0.5, 

274 name=glyphs.sqrt(glyphs.scalar(value))) 

275 

276 if len(root) == 1: 

277 P.add_constraint(AbsoluteValueConstraint(root, bound)) 

278 else: 

279 P.add_constraint(SOCConstraint(root, bound)) 

280 

281 return P 

282 

283 def __init__(self, lhs, relation, rhs): 

284 """Construct a :class:`ConvexQuadraticConstraint`. 

285 

286 See :meth:`NonconvexQuadraticConstraint.__init__` for more. 

287 """ 

288 NonconvexQuadraticConstraint.__init__(self, lhs, relation, rhs) 

289 

290 def _get_type_term(self): 

291 return "Convex Quadratic" 

292 

293 Subtype = namedtuple("Subtype", ("qterms", "rank", "haslin")) 

294 

295 def _subtype(self): 

296 from ..expressions import QuadraticExpression 

297 

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 

305 

306 return self.Subtype(self.le0.num_quad_terms, self.le0.rank, haslin) 

307 

308 

309# -------------------------------------- 

310__all__ = api_end(_API_START, globals())