Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1# ------------------------------------------------------------------------------ 

2# Copyright (C) 2019 Maximilian Stahlberg 

3# Based on the original picos.expressions module by Guillaume Sagnol. 

4# 

5# This file is part of PICOS. 

6# 

7# PICOS is free software: you can redistribute it and/or modify it under the 

8# terms of the GNU General Public License as published by the Free Software 

9# Foundation, either version 3 of the License, or (at your option) any later 

10# version. 

11# 

12# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY 

13# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 

14# A PARTICULAR PURPOSE. See the GNU General Public License for more details. 

15# 

16# You should have received a copy of the GNU General Public License along with 

17# this program. If not, see <http://www.gnu.org/licenses/>. 

18# ------------------------------------------------------------------------------ 

19 

20"""Implements :class:`QuadraticExpression`.""" 

21 

22import operator 

23import sys 

24from collections import OrderedDict as odict 

25from collections import namedtuple 

26from types import MappingProxyType 

27 

28import cvxopt 

29import cvxopt.cholmod 

30import numpy 

31 

32from .. import glyphs 

33from ..apidoc import api_end, api_start 

34from ..caching import (borrow_cache, cached_property, 

35 cached_selfinverse_unary_operator, 

36 cached_unary_operator) 

37from ..constraints import (ConicQuadraticConstraint, Constraint, 

38 ConvexQuadraticConstraint, 

39 NonconvexQuadraticConstraint) 

40from ..legacy import deprecated 

41from .data import convert_operands, cvx2np, should_be_sparse 

42from .exp_affine import AffineExpression, ComplexAffineExpression, Constant 

43from .expression import Expression, refine_operands, validate_prediction 

44from .variables import BaseVariable 

45 

46_API_START = api_start(globals()) 

47# ------------------------------- 

48 

49 

50PSD_PERTURBATIONS = 3 

51r"""Maximum number of singular quadratic form perturbations. 

52 

53PICOS uses NumPy either or CHOLMOD to compute a Cholesky decomposition of a 

54positive semidefinite quadratic form :math:`Q`. Both libraries require 

55:math:`Q` to be nonsingular, which is not a requirement on PICOS' end. 

56If either library rejects :math:`Q`, then PICOS provides a sequence of 

57:data:`PSD_PERTURBATIONS` perturbed quadratic forms :math:`Q + \epsilon I` for 

58increasing :math:`\epsilon` until the perturbed matrix is found positive 

59definite or until the largest :math:`\epsilon` was tested unsuccessfully. 

60 

61If this is zero, PICOS will only decompose quadratic forms that are nonsingular. 

62If this is one, then only the largest epsilon is tested. 

63""" 

64 

65 

66class QuadraticExpression(Expression): 

67 """A scalar quadratic expression of the form :math:`x^TQx + a^Tx + b`.""" 

68 

69 # -------------------------------------------------------------------------- 

70 # Initialization and factory methods. 

71 # -------------------------------------------------------------------------- 

72 

73 def __init__( 

74 self, string, quadraticPart={}, affinePart=AffineExpression.zero(), 

75 scalarFactors=None, copyDecomposition=None): 

76 r"""Initialize a scalar quadratic expression. 

77 

78 This constructor is meant for internal use. As a user, you will most 

79 likely want to build expressions starting from 

80 :mod:`~picos.expressions.variables` or a :func:`~picos.Constant`. 

81 

82 :param str string: A symbolic string description. 

83 :param quadraticPart: A :class:`dict` mapping PICOS variable pairs to 

84 CVXOPT matrices. Each entry :math:`(x, y) \mapsto A_{xy}` represents 

85 a quadratic form 

86 :math:`\operatorname{vec_x}(x)^T A_{xy} \operatorname{vec_y}(y)` 

87 where :math:`\operatorname{vec_x}` and :math:`\operatorname{vec_y}` 

88 refer to the isometric real 

89 :mod:`~.picos.expressions.vectorizations` that are used by PICOS 

90 internally to store the variables :math:`x` and :math:`y`, 

91 respectively. The quadratic part :math:`Q` of the expression is then 

92 given as :math:`Q = \sum_{x, y} A_{xy}`. 

93 :param affinePart: The affine part :math:`a^Tx + b` of the expression. 

94 :type affinePart: ~picos.expressions.AffineExpression 

95 :param scalarFactors: A pair :math:`(u, v)` with both :math:`u` and 

96 :math:`v` scalar real affine expressions representing a known 

97 :math:`x^TQx + a^Tx + b = uv` factorization of the expression. 

98 :type factorization: 

99 tuple(~picos.expressions.AffineExpression) 

100 :param copyDecomposition: Another quadratic expression with equal 

101 quadratic part whose quadratic part decomposition shall be copied. 

102 :type copyDecomposition: 

103 ~picos.expressions.QuadraticExpression 

104 """ 

105 # Ensure correct usage within PICOS. 

106 assert isinstance(quadraticPart, dict) and all( 

107 isinstance(xy, (tuple, list)) and len(xy) == 2 

108 and isinstance(xy[0], BaseVariable) 

109 and isinstance(xy[1], BaseVariable) 

110 and isinstance(A, (cvxopt.matrix, cvxopt.spmatrix)) 

111 and A.size == (xy[0].dim, xy[1].dim) 

112 for xy, A in quadraticPart.items()) 

113 assert isinstance(affinePart, ComplexAffineExpression) \ 

114 and len(affinePart) == 1 

115 assert not copyDecomposition \ 

116 or isinstance(copyDecomposition, QuadraticExpression) 

117 assert not scalarFactors or ( 

118 isinstance(scalarFactors, (tuple, list)) and len(scalarFactors) == 2 

119 and all(isinstance(x, AffineExpression) for x in scalarFactors) 

120 and all(len(x) == 1 for x in scalarFactors)) 

121 

122 Expression.__init__(self, "Quadratic Expression", string) 

123 

124 # Check affine part. 

125 affinePart = affinePart.refined 

126 if not isinstance(affinePart, AffineExpression): 

127 raise NotImplementedError( 

128 "PICOS does not support complex-valued quadratic expressions " 

129 "and the affine part of {} is not real.".format(self._symbStr)) 

130 elif len(affinePart) != 1: 

131 raise TypeError("Affine part of {} is not scalar." 

132 .format(self._symbStr)) 

133 

134 # Store quadratic part. 

135 self._quads = {} 

136 for xy, A in quadraticPart.items(): 

137 # Ignore coefficients that are already zero. 

138 if not A: 

139 continue 

140 

141 if xy[0] is xy[1]: 

142 # Store the unique symmetric form for two reasons: 

143 # - Hermitian forms concerning the original variables are 

144 # supplied as quadratic forms concerning the isometric real 

145 # vectorizations of the variables. These matrices are in 

146 # general still complex, but their symmetric form is not. 

147 # - Coefficients that cancel out are set to numeric zeros and 

148 # can be converted to structural zeros by _sparse_quads. 

149 A = 0.5*(A + A.T) 

150 if hash(xy[1]) < hash(xy[0]): 

151 # Store only one coefficient per unordered pair. 

152 xy, A = (xy[1], xy[0]), A.T 

153 

154 if xy in self._quads: 

155 self._quads[xy] = self._quads[xy] + A 

156 else: 

157 self._quads[xy] = A 

158 

159 # Remove coefficients of zero. 

160 self._quads = {xy: A for xy, A in self._quads.items() if A} 

161 

162 # Check if coefficients are real. 

163 for (x, y), A in self._quads.items(): 

164 if A.typecode == "z": 

165 if any(A.imag()): 

166 raise NotImplementedError( 

167 "PICOS does not support complex-valued quadratic " 

168 "expressions and the quadratic part of {} for variables" 

169 " {} and {} is not real or equivalent to a real form." 

170 .format(self._symbStr, x.string, y.string)) 

171 

172 self._quads[xy] = A.real() 

173 

174 # Store affine part. 

175 self._aff = affinePart 

176 

177 # Store a known factorization into real scalars. 

178 if scalarFactors: 

179 a, b = scalarFactors 

180 if a.equals(b): 

181 self._sf = (a, a) # Speedup future equality checks. 

182 else: 

183 self._sf = (a, b) 

184 else: 

185 self._sf = None 

186 

187 # Copy a decomposition from another quadratic expression. 

188 if copyDecomposition is not None: 

189 borrow_cache(self, copyDecomposition, ("_Q_and_x", "L", "quadroot")) 

190 

191 # -------------------------------------------------------------------------- 

192 # Abstract method implementations and method overridings, except _predict. 

193 # -------------------------------------------------------------------------- 

194 

195 @cached_unary_operator 

196 def _get_refined(self): 

197 if not self._quads: 

198 return self._aff.refined.renamed(self._symbStr) 

199 else: 

200 return self 

201 

202 Subtype = namedtuple("Subtype", ( 

203 "convex", "concave", "qterms", "rank", "haslin", "factors")) 

204 

205 @cached_unary_operator 

206 def _get_subtype(self): 

207 convex = self.convex 

208 concave = False if convex and self._quads else self.concave 

209 qterms = self.num_quad_terms 

210 haslin = not self._aff.constant 

211 

212 try: 

213 if convex: 

214 rank = self.rank 

215 elif concave: 

216 rank = (-self).rank 

217 else: 

218 rank = None 

219 except ValueError: # No quadratic part. 

220 rank = 0 

221 

222 factors = len(set(self._sf)) if self._sf else 0 

223 

224 return self.Subtype(convex, concave, qterms, rank, haslin, factors) 

225 

226 def _get_value(self): 

227 value = self._aff._get_value() 

228 

229 for xy, A in self._quads.items(): 

230 x, y = xy 

231 xVal = x._get_internal_value() 

232 yVal = y._get_internal_value() 

233 value = value + xVal.T * A * yVal 

234 

235 return value 

236 

237 @cached_unary_operator 

238 def _get_mutables(self): 

239 return self._aff._get_mutables().union( 

240 var for vars in self._quads for var in vars) 

241 

242 def _is_convex(self): 

243 try: 

244 self.L 

245 except ValueError: 

246 return True # Expression is affine. 

247 except ArithmeticError: 

248 return False 

249 else: 

250 return True 

251 

252 def _is_concave(self): 

253 try: 

254 (-self).L 

255 except ValueError: 

256 return True # Expression is affine. 

257 except ArithmeticError: 

258 return False 

259 else: 

260 return True 

261 

262 def _replace_mutables(self, mapping): 

263 name_map = {old.name: new.name for old, new in mapping.items()} 

264 

265 string = self.string 

266 if isinstance(string, glyphs.GlStr): 

267 string = string.reglyphed(name_map) 

268 

269 quads = {(mapping[var1], mapping[var2]): coef 

270 for (var1, var2), coef in self._quads.items()} 

271 

272 if not self._sf: 

273 sf = None 

274 elif self._sf[0] is self._sf[1]: 

275 sf = (self._sf[0]._replace_mutables(mapping),)*2 

276 else: 

277 sf = tuple(f._replace_mutables(mapping) for f in self._sf) 

278 

279 return QuadraticExpression( 

280 string, quads, self.aff._replace_mutables(mapping), sf) 

281 

282 def _freeze_mutables(self, freeze): 

283 # TODO: Allow freezing of quadratic expressions. 

284 raise NotImplementedError( 

285 "Partially freezing quadratic expressions is not yet supported.") 

286 

287 # -------------------------------------------------------------------------- 

288 # Python special method implementations, except constraint-creating ones. 

289 # -------------------------------------------------------------------------- 

290 

291 def __len__(self): 

292 # Faster version that overrides Expression.__len__. 

293 return 1 

294 

295 @convert_operands(sameShape=True) 

296 def __add__(self, other): 

297 """Denote addition from the right hand side.""" 

298 if isinstance(other, QuadraticExpression): 

299 string = glyphs.clever_add(self.string, other.string) 

300 

301 quadPart = {} 

302 varPairs = set(self._quads.keys()).union(other._quads.keys()) 

303 for vars in varPairs: 

304 if vars in self._quads: 

305 if vars in other._quads: 

306 quadPart[vars] = self._quads[vars] + other._quads[vars] 

307 else: 

308 quadPart[vars] = self._quads[vars] 

309 else: 

310 quadPart[vars] = other._quads[vars] 

311 

312 affPart = self._aff + other._aff 

313 

314 return QuadraticExpression(string, quadPart, affPart) 

315 elif isinstance(other, AffineExpression): 

316 string = glyphs.clever_add(self.string, other.string) 

317 

318 return QuadraticExpression( 

319 string, self._quads, self._aff + other, copyDecomposition=self) 

320 else: 

321 return NotImplemented 

322 

323 @convert_operands(sameShape=True) 

324 def __radd__(self, other): 

325 """Denote addition from the left hand side.""" 

326 if isinstance(other, QuadraticExpression): 

327 assert False, "Impossible case." 

328 elif isinstance(other, AffineExpression): 

329 result = self.__add__(other) 

330 # NOTE: __add__ always creates a fresh expression. 

331 result._symbStr = glyphs.clever_add(other.string, self.string) 

332 return result 

333 else: 

334 return NotImplemented 

335 

336 @convert_operands(sameShape=True) 

337 def __sub__(self, other): 

338 """Denote substraction from the right hand side.""" 

339 if isinstance(other, QuadraticExpression): 

340 string = glyphs.clever_sub(self.string, other.string) 

341 

342 quadPart = {} 

343 varPairs = set(self._quads.keys()).union(other._quads.keys()) 

344 for vars in varPairs: 

345 if vars in self._quads: 

346 if vars in other._quads: 

347 quadPart[vars] = self._quads[vars] - other._quads[vars] 

348 else: 

349 quadPart[vars] = self._quads[vars] 

350 else: 

351 quadPart[vars] = -other._quads[vars] 

352 

353 affPart = self._aff - other._aff 

354 

355 return QuadraticExpression(string, quadPart, affPart) 

356 elif isinstance(other, AffineExpression): 

357 string = glyphs.clever_sub(self.string, other.string) 

358 

359 return QuadraticExpression( 

360 string, self._quads, self._aff - other, copyDecomposition=self) 

361 else: 

362 return NotImplemented 

363 

364 @convert_operands(sameShape=True) 

365 def __rsub__(self, other): 

366 """Denote substraction with self on the right hand side.""" 

367 if isinstance(other, QuadraticExpression): 

368 assert False, "Impossible case." 

369 elif isinstance(other, AffineExpression): 

370 result = (-self).__add__(other) 

371 # NOTE: __add__ always creates a fresh expression. 

372 result._symbStr = glyphs.clever_sub(other.string, self.string) 

373 return result 

374 else: 

375 return NotImplemented 

376 

377 @cached_selfinverse_unary_operator 

378 def __neg__(self): 

379 string = glyphs.clever_neg(self.string) 

380 quadPart = {xy: -A for xy, A in self._quads.items()} 

381 affPart = -self._aff 

382 factors = (self._sf[0], -self._sf[1]) if self._sf else None 

383 

384 return QuadraticExpression(string, quadPart, affPart, factors) 

385 

386 @convert_operands(scalarRHS=True) 

387 def __mul__(self, other): 

388 """Denote scaling from the right hand side.""" 

389 if isinstance(other, AffineExpression): 

390 if not other.constant: 

391 raise NotImplementedError("You may only multiply a PICOS " 

392 "quadratic expression with a constant term.") 

393 

394 factor = other.value 

395 string = glyphs.clever_mul(self.string, other.string) 

396 quadPart = {xy: factor*A for xy, A in self._quads.items()} 

397 affPart = self._aff*other 

398 

399 if self._sf: 

400 r = factor**0.5 

401 if self._sf[0] is self._sf[1]: 

402 factors = (r*self._sf[0],)*2 

403 else: 

404 factors = (r*self._sf[0], r*self._sf[1]) 

405 else: 

406 factors = None 

407 

408 return QuadraticExpression(string, quadPart, affPart, factors) 

409 else: 

410 return NotImplemented 

411 

412 @convert_operands(scalarRHS=True) 

413 def __rmul__(self, other): 

414 """Denote scaling from the left hand side.""" 

415 if isinstance(other, AffineExpression): 

416 result = self.__mul__(other) 

417 # NOTE: __mul__ always creates a fresh expression. 

418 result._symbStr = glyphs.clever_mul(other.string, self.string) 

419 return result 

420 else: 

421 return NotImplemented 

422 

423 @convert_operands(scalarRHS=True) 

424 def __truediv__(self, other): 

425 """Denote division by a constant scalar.""" 

426 if isinstance(other, AffineExpression): 

427 if not other.constant: 

428 raise NotImplementedError("You may only divide a PICOS " 

429 "quadratic expression by a constant term.") 

430 

431 result = self.__mul__(1.0 / other.value) 

432 # NOTE: __mul__ always creates a fresh expression. 

433 result._symbStr = glyphs.div(self.string, other.string) 

434 return result 

435 else: 

436 return NotImplemented 

437 

438 # -------------------------------------------------------------------------- 

439 # Decomposition related methods. 

440 # -------------------------------------------------------------------------- 

441 

442 @cached_property 

443 def _sparse_quads(self): 

444 """Quadratic coefficients (re-)cast to sparse matrices.""" 

445 return {xy: cvxopt.sparse(M) for xy, M in self._quads.items()} 

446 

447 def _Q_and_x_or_A_and_y(self, augmentedMatrix): 

448 """Either :attr:`_Q_and_x` or :attr:`_A_and_y`.""" 

449 # Find occurences of the scalar entries of the variable's isometric real 

450 # vectorization in the quadratic part. 

451 nonzero = {} 

452 for (x, y), A in self._sparse_quads.items(): 

453 nonzero.setdefault(x, set()) 

454 nonzero.setdefault(y, set()) 

455 

456 nonzero[x].update(A.I) 

457 nonzero[y].update(A.J) 

458 

459 # For the augmented matrix, find additional linear occurences. 

460 if augmentedMatrix: 

461 for x, a in self._aff._sparse_linear_coefs.items(): 

462 nonzero.setdefault(x, set()) 

463 nonzero[x].update(a.I) 

464 

465 nonzero = odict((var, sorted(nonzero[var])) 

466 for var in sorted(nonzero, key=lambda v: v.id)) 

467 

468 # Compute each (multidimensional) variable's offset in Q/A. 

469 offset, offsets = 0, {} 

470 for var, nz in nonzero.items(): 

471 offsets[var] = offset 

472 offset += len(nz) 

473 

474 # Compute Q. 

475 V, I, J = [], [], [] 

476 for (x, y), A in self._sparse_quads.items(): 

477 B = A[nonzero[x], nonzero[y]] 

478 

479 V.extend(B.V) 

480 I.extend(offsets[x] + i for i in B.I) 

481 J.extend(offsets[y] + j for j in B.J) 

482 

483 # Turn Q into A if requested. 

484 if augmentedMatrix: 

485 for x, a in self._aff._sparse_linear_coefs.items(): 

486 b = a[nonzero[x]] 

487 

488 V.extend(b.V) 

489 I.extend(offsets[x] + i for i in b.I) 

490 J.extend([offset]*len(b.J)) 

491 

492 V.append(self._aff._constant_coef[0]) 

493 I.append(offset) 

494 J.append(offset) 

495 

496 offset += 1 

497 

498 # Finalize Q/A. 

499 Q = cvxopt.spmatrix(V, I, J, (offset,)*2, tc="d") 

500 Q = 0.5*(Q + Q.T) 

501 

502 if not Q: 

503 if augmentedMatrix: 

504 raise ValueError( 

505 "The expression {} is zero.".format(self._symbStr)) 

506 else: 

507 raise ValueError( 

508 "The quadratic part of {} is zero.".format(self._symbStr)) 

509 

510 # Compute x. 

511 x = None 

512 for var, nz in nonzero.items(): 

513 n, r = len(nz), var.dim 

514 F = cvxopt.spmatrix([1]*n, range(n), nz, (n, r)) 

515 y = AffineExpression( 

516 glyphs.Fn("F")(var.string), n, coefficients={var: F}) 

517 x = y if x is None else x // y 

518 

519 # Turn x into y if requested. 

520 if augmentedMatrix: 

521 x = AffineExpression.from_constant(1) if x is None else x // 1 

522 

523 return Q, x 

524 

525 @cached_property 

526 def _Q_and_x(self): 

527 """Pair containing :attr:`Q` and :attr:`x`.""" 

528 return self._Q_and_x_or_A_and_y(augmentedMatrix=False) 

529 

530 @cached_property 

531 def _A_and_y(self): 

532 """Pair containing :attr:`R` and :attr:`y`.""" 

533 return self._Q_and_x_or_A_and_y(augmentedMatrix=True) 

534 

535 @property 

536 def Q(self): 

537 """The coefficient matrix :math:`Q` of the expression, condensed. 

538 

539 This equals the :math:`Q` of the quadratic expression 

540 :math:`x^TQx + a^Tx + b` but with zero rows and columns removed. 

541 

542 The vector :attr:`x` is a condensed version of :math:`x` that refers to 

543 this matrix, so that ``q.x.T*q.Q*q.x`` equals the quadratic part 

544 :math:`x^TQx` of the expression ``q``. 

545 

546 :raises ValueError: When the quadratic part is zero. 

547 """ 

548 return self._Q_and_x[0] 

549 

550 @property 

551 def x(self): 

552 """The stacked variable vector :math:`x` of the expression, condensed. 

553 

554 This equals the :math:`x` of the quadratic expression 

555 :math:`x^TQx + a^Tx + b` but entries corresponding to zero rows and 

556 columns in :attr:`Q` are removed, so that ``q.x.T*q.Q*q.x`` equals the 

557 :math:`x^TQx` part of the expression ``q``. 

558 

559 :raises ValueError: When the quadratic part is zero. 

560 """ 

561 return self._Q_and_x[1] 

562 

563 @property 

564 def A(self): 

565 r"""An affine-augmented quadratic coefficient matrix, condensed. 

566 

567 For a quadratic expression :math:`x^TQx + a^Tx + b`, this is 

568 :math:`A = \begin{bmatrix}Q&\frac{a}{2}\\\frac{a^T}{2}&b\end{bmatrix}` 

569 but with zero rows and columns removed. 

570 

571 The vector :attr:`y` is a condensed version of :math:`x` that refers to 

572 this matrix, so that ``q.y.T*q.A*q.y`` equals the expression ``q``. 

573 

574 :raises ValueError: When the expression is zero. 

575 """ 

576 return self._A_and_y[0] 

577 

578 @property 

579 def y(self): 

580 """See :attr:`A`. 

581 

582 :raises ValueError: When the expression is zero. 

583 """ 

584 return self._A_and_y[1] 

585 

586 def _L_or_M(self, augmentedMatrix): 

587 """Either :attr:`L` or :attr:`M`.""" 

588 Q = self.A if augmentedMatrix else self.Q 

589 n = Q.size[0] 

590 

591 # Define a sequence of small numbers for perturbation. 

592 # TODO: These numbers are empirical; find worst case values. 

593 spread = [abs(v) for v in Q.V if v] 

594 minEps = 1 * min(spread) * sys.float_info.epsilon 

595 maxEps = n * max(spread) * sys.float_info.epsilon 

596 epsilons = set([0.0]) 

597 for i in range(PSD_PERTURBATIONS): 

598 epsilons.add(minEps + i * (maxEps - minEps) / PSD_PERTURBATIONS) 

599 epsilons = sorted(epsilons) 

600 

601 # Check if Q is really sparse. 

602 sparse = should_be_sparse(Q.size, len(Q)) 

603 

604 # Attempt to find L. 

605 L = None 

606 if sparse: # Use CHOLMOD. 

607 F = None 

608 for eps in epsilons: 

609 P = Q + cvxopt.spdiag([eps]*n) if eps else Q # Perturbed Q. 

610 

611 if F is None: 

612 F = cvxopt.cholmod.symbolic(P) 

613 

614 try: 

615 cvxopt.cholmod.numeric(P, F) 

616 except ArithmeticError: 

617 pass 

618 else: 

619 L = cvxopt.cholmod.getfactor(F) 

620 break 

621 else: # Use NumPy. 

622 Q = cvx2np(Q) 

623 for eps in epsilons: 

624 P = Q + numpy.diag([eps]*n) if eps else Q # Perturbed Q. 

625 

626 try: 

627 L = numpy.linalg.cholesky(P) 

628 except numpy.linalg.LinAlgError: 

629 pass 

630 else: 

631 break 

632 

633 if L is None: 

634 if not PSD_PERTURBATIONS: 

635 raise ArithmeticError("The expression {} has an (augmented) " 

636 "quadratic form that is either not positive semidefinite " 

637 "or singular. Try enabling PSD_PERTURBATIONS.") 

638 elif augmentedMatrix: 

639 raise ArithmeticError("The expression {} is not a squared norm " 

640 "as its {} representation is not positive semidefinite." 

641 .format(self._symbStr, glyphs.matrix(glyphs.vertcat( 

642 glyphs.horicat("Q", glyphs.div("a", 2)), 

643 glyphs.horicat(glyphs.div(glyphs.transp("a"), 2), "b") 

644 )))) 

645 else: 

646 raise ArithmeticError("The expression {} is not convex as its " 

647 "quadratic part is not positive semidefinite." 

648 .format(self._symbStr)) 

649 

650 # Remove near-zeros in the order of n·sqrt(maxEps). 

651 cutoff = 2 * n * maxEps**0.5 

652 if sparse: 

653 VIJ = list(t for t in zip(L.V, L.I, L.J) if abs(t[0]) > cutoff) 

654 if VIJ: # Don't zero the matrix. 

655 V, I, J = zip(*VIJ) 

656 L = cvxopt.spmatrix(V, I, J, L.size) 

657 else: 

658 mask = abs(L) <= cutoff 

659 if not numpy.all(mask): # Don't zero the matrix. 

660 L[mask] = 0 

661 L = cvxopt.sparse(cvxopt.matrix(L)) 

662 

663 return L 

664 

665 @cached_property 

666 def L(self): 

667 """The :math:`L` of an :math:`LL^T` Cholesky decomposition of :attr:`Q`. 

668 

669 :returns: A CVXOPT lower triangular sparse matrix. 

670 

671 :raises ValueError: When the quadratic part is zero. 

672 :raises ArithmeticError: When the expression is not convex, that is when 

673 the matrix :attr:`Q` is not numerically positive semidefinite. 

674 """ 

675 return self._L_or_M(augmentedMatrix=False) 

676 

677 @cached_property 

678 def M(self): 

679 """The :math:`M` of an :math:`MM^T` Cholesky decomposition of :attr:`A`. 

680 

681 :returns: A CVXOPT lower triangular sparse matrix. 

682 

683 :raises ValueError: When the expression is zero. 

684 :raises ArithmeticError: When the expression can not be written as a 

685 squared norm, that is when the extended matrix :attr:`A` is not 

686 numerically positive semidefinite. 

687 """ 

688 return self._L_or_M(augmentedMatrix=True) 

689 

690 @cached_property 

691 def quadroot(self): 

692 r"""Affine expression whose squared norm equals :math:`x^TQx`. 

693 

694 For a convex quadratic expression ``q``, this is equal to the vector 

695 ``q.L.T*q.x`` with zero rows removed. 

696 

697 :Construction: 

698 

699 Let :math:`x^TQx` be the quadratic part of the expression with :math:`Q` 

700 positive semidefinite and :math:`Q = LL^T` a Cholesky decomposition. 

701 Then, 

702 

703 .. math:: 

704 x^TQx 

705 &= x^TLL^Tx \\ 

706 &= (L^Tx)^TL^Tx \\ 

707 &= \langle L^Tx, L^Tx \rangle \\ 

708 &= \lVert L^Tx \rVert^2. 

709 

710 Note that removing zero rows from :math:`L^Tx` does not affect the norm. 

711 

712 :raises ValueError: When the quadratic part is zero. 

713 :raises ArithmeticError: When the expression is not convex, that is when 

714 the quadratic part is not numerically positive semidefinite. 

715 """ 

716 L = self.L 

717 n = L.size[0] 

718 nz = sorted(set(L.J)) 

719 nnz = len(nz) 

720 

721 # F removes rows corresponding to zero columns in L. 

722 F = cvxopt.spmatrix([1.0]*nnz, range(nnz), nz, (nnz, n)) 

723 

724 result = (F*L.T)*self.x 

725 # NOTE: AffineExpression.__mul__ always creates a fresh expression. 

726 result._symbStr = glyphs.Fn("quadroot({})")(self.string) 

727 return result 

728 

729 @cached_property 

730 def fullroot(self): 

731 r"""Affine expression whose squared norm equals the expression. 

732 

733 For a convex quadratic expression ``q``, this is equal to the vector 

734 ``q.M.T*q.y`` with zero rows removed. 

735 

736 :Construction: 

737 

738 For a quadratic expression :math:`x^TQx + a^Tx + b` with 

739 :math:`A = \begin{bmatrix}Q&\frac{a}{2}\\\frac{a^T}{2}&b\end{bmatrix}` 

740 positive semidefinite, let :math:`A = MM^T` be a Cholesky decomposition. 

741 Let further :math:`y = \begin{bmatrix}x^T & 1\end{bmatrix}^T`. Then, 

742 

743 .. math:: 

744 x^TQx + a^Tx + b 

745 &= y^TAy \\ 

746 &= y^TMM^Ty \\ 

747 &= (M^Ty)^TM^Ty \\ 

748 &= \langle M^Ty, M^Ty \rangle \\ 

749 &= \lVert M^Ty \rVert^2. 

750 

751 Note that removing zero rows from :math:`M^Ty` does not affect the norm. 

752 

753 :raises ValueError: When the expression is zero. 

754 :raises ArithmeticError: When the expression is not convex, that is when 

755 the quadratic part is not numerically positive semidefinite. 

756 """ 

757 M = self.M 

758 n = M.size[0] 

759 nz = sorted(set(M.J)) 

760 nnz = len(nz) 

761 

762 # F removes rows corresponding to zero columns in M. 

763 F = cvxopt.spmatrix([1.0]*nnz, range(nnz), nz, (nnz, n)) 

764 

765 result = (F*M.T)*self.y 

766 # NOTE: AffineExpression.__mul__ always creates a fresh expression. 

767 result._symbStr = glyphs.Fn("fullroot({})")(self.string) 

768 return result 

769 

770 @property 

771 def rank(self): 

772 """The length of the vector :attr:`quadroot`. 

773 

774 Up to numerical considerations, this is the rank of the (convex) 

775 quadratic coefficient matrix :math:`Q` of the expression. 

776 

777 :raises ArithmeticError: When the expression is not convex, that is when 

778 the quadratic part is not numerically positive semidefinite. 

779 """ 

780 try: 

781 return len(self.quadroot) 

782 except ValueError: 

783 return 0 # No quadratic part. 

784 

785 @cached_property 

786 def num_quad_terms(self): 

787 """The number of terms in the simplified quadratic form.""" 

788 num_terms = 0 

789 for A in self._sparse_quads.values(): 

790 for i, j in zip(A.I, A.J): 

791 if i >= j: 

792 num_terms += 1 

793 return num_terms 

794 

795 @property 

796 def is_squared_norm(self): 

797 r"""Whether the expression can be written as a squared norm. 

798 

799 If this is :obj:`True`, then the there is a coefficient vector 

800 :math:`c` such that 

801 

802 .. math:: 

803 \lVert c^T \begin{bmatrix} x \\ 1 \end{bmatrix} \rVert^2 

804 = x^TQx + a^Tx + b. 

805 

806 If the expression is also nonzero, then :attr:`fullroot` is 

807 :math:`c^T \begin{bmatrix} x^T & 1 \end{bmatrix}^T` with zero entries 

808 removed. 

809 """ 

810 try: 

811 self.M 

812 except ValueError: 

813 assert self.is0 

814 return True 

815 except ArithmeticError: 

816 return False 

817 else: 

818 return True 

819 

820 # -------------------------------------------------------------------------- 

821 # Properties and functions that describe the expression. 

822 # -------------------------------------------------------------------------- 

823 

824 @property 

825 @deprecated("2.2", "This property will be removed in a future release.") 

826 def quadratic_forms(self): 

827 """The quadratic forms as a map from variable pairs to sparse matrices. 

828 

829 .. warning:: 

830 

831 Do not modify the returned matrices. 

832 """ 

833 return MappingProxyType(self._sparse_quads) 

834 

835 @property 

836 def is0(self): 

837 """Whether the quadratic expression is zero.""" 

838 return not self._quads and self._aff.is0 

839 

840 @property 

841 def scalar_factors(self): 

842 r"""Decomposition into scalar real affine expressions. 

843 

844 If the expression is known to be equal to :math:`ab` for scalar real 

845 affine expressions :math:`a` and :math:`b`, this is the pair ``(a, b)``. 

846 Otherwise, this is :obj:`None`. 

847 

848 Note that if :math:`a = b`, then also ``a is b`` in the returned tuple. 

849 """ 

850 return self._sf 

851 

852 # -------------------------------------------------------------------------- 

853 # Methods and properties that return modified copies. 

854 # -------------------------------------------------------------------------- 

855 

856 @cached_property 

857 def quad(self): 

858 """Quadratic part :math:`x^TQx` of the quadratic expression.""" 

859 factors = self._sf if self._sf and self._aff.is0 else None 

860 

861 return QuadraticExpression(glyphs.quadpart(self._symbStr), self._quads, 

862 scalarFactors=factors, copyDecomposition=self) 

863 

864 @cached_property 

865 def aff(self): 

866 """Affine part :math:`a^Tx + b` of the quadratic expression.""" 

867 return self._aff.renamed(glyphs.affpart(self._symbStr)) 

868 

869 @cached_property 

870 def lin(self): 

871 """Linear part :math:`a^Tx` of the quadratic expression.""" 

872 return self._aff.lin.renamed(glyphs.linpart(self._symbStr)) 

873 

874 @cached_property 

875 def cst(self): 

876 """Constant part :math:`b` of the quadratic expression.""" 

877 return self._aff.cst.renamed(glyphs.cstpart(self._symbStr)) 

878 

879 # -------------------------------------------------------------------------- 

880 # Constraint-creating operators, and _predict. 

881 # -------------------------------------------------------------------------- 

882 

883 @classmethod 

884 def _predict(cls, subtype, relation, other): 

885 assert isinstance(subtype, cls.Subtype) 

886 

887 if issubclass(other.clstype, QuadraticExpression): 

888 raise NotImplementedError("Constraint outcome prediction not " 

889 "supported when comparing two quadratic expressions.") 

890 

891 if relation == operator.__le__: 

892 if issubclass(other.clstype, AffineExpression) \ 

893 and other.subtype.dim == 1: 

894 if subtype.convex: 

895 return ConvexQuadraticConstraint.make_type( 

896 qterms=subtype.qterms, rank=subtype.rank, 

897 haslin=(subtype.haslin or not other.subtype.constant)) 

898 else: 

899 return NonconvexQuadraticConstraint.make_type( 

900 qterms=subtype.qterms) 

901 elif relation == operator.__ge__: 

902 if issubclass(other.clstype, AffineExpression) \ 

903 and other.subtype.dim == 1: 

904 if subtype.concave: 

905 return ConvexQuadraticConstraint.make_type( 

906 qterms=subtype.qterms, rank=subtype.rank, 

907 haslin=(subtype.haslin or not other.subtype.constant)) 

908 elif subtype.factors and other.subtype.constant \ 

909 and other.subtype.nonneg: 

910 # See __ge__ for the variants below. 

911 

912 # Variant 1: 

913 return ConicQuadraticConstraint.make_type( 

914 qterms=subtype.qterms, conic_argdim=1, 

915 rotated=(subtype.factors == 2)) 

916 

917 # Variant 2: 

918 # return RSOCConstraint.make_type(argdim=1) 

919 

920 # Variant 3: 

921 # if subtype.factors == 1: 

922 # return AffineConstraint.make_type(dim=1, eq=False) 

923 # else: 

924 # return RSOCConstraint.make_type(argdim=1) 

925 else: 

926 return NonconvexQuadraticConstraint.make_type( 

927 qterms=subtype.qterms) 

928 

929 return NotImplemented 

930 

931 @convert_operands(sameShape=True) 

932 @validate_prediction 

933 @refine_operands() 

934 def __le__(self, other): 

935 LE = Constraint.LE 

936 

937 if isinstance(other, QuadraticExpression): 

938 le0 = self - other # Unpredictable outcome. 

939 if le0.convex: 

940 return ConvexQuadraticConstraint(self, LE, other) 

941 elif self.is_squared_norm and other.scalar_factors: 

942 return ConicQuadraticConstraint(self, LE, other) 

943 else: 

944 return NonconvexQuadraticConstraint(self, LE, other) 

945 elif isinstance(other, AffineExpression): 

946 if self.convex: 

947 return ConvexQuadraticConstraint(self, LE, other) 

948 else: 

949 return NonconvexQuadraticConstraint(self, LE, other) 

950 else: 

951 return NotImplemented 

952 

953 @convert_operands(sameShape=True) 

954 @validate_prediction 

955 @refine_operands() 

956 def __ge__(self, other): 

957 GE = Constraint.GE 

958 

959 if isinstance(other, QuadraticExpression): 

960 le0 = other - self # Unpredictable outcome. 

961 if le0.convex: 

962 return ConvexQuadraticConstraint(self, GE, other) 

963 elif other.is_squared_norm and self.scalar_factors: 

964 return ConicQuadraticConstraint(self, GE, other) 

965 else: 

966 return NonconvexQuadraticConstraint(self, GE, other) 

967 elif isinstance(other, AffineExpression): 

968 if self.concave: 

969 return ConvexQuadraticConstraint(self, GE, other) 

970 elif self._sf and other.constant and other.value >= 0: 

971 # Handle the case of c <= x*x and c <= x*y with c pos. constant. 

972 root = Constant(glyphs.sqrt(other.string), other.value**0.5) 

973 

974 # Variant 1: Return a ConicQuadraticConstraint that can be 

975 # reformulated to an (R)SOCConstraint. This is consistent 

976 # with all other quadratic expression comparison outcomes. 

977 other = QuadraticExpression(other.string, affinePart=other, 

978 scalarFactors=(root, root)) # Unrefined QE. 

979 return ConicQuadraticConstraint(self, GE, other) 

980 

981 # Variant 2: Return an RSOCConstraint. This is somewhat 

982 # consistent with Norm returning an SOCConstraint and it is 

983 # also legacy behavior. However this may not be what the 

984 # user wants, because it silently forces x, y >= 0! 

985 # return RSOCConstraint(root, self._sf[0], self._sf[1], 

986 # customString = glyphs.le(other.string, self.string)) 

987 

988 # Variant 3: Return either an RSOCConstraint or an affine 

989 # constraint that equals an SOCConstraint. This is faster 

990 # than Variant 2 but even more confusing, because the 

991 # constraint would transform from e.g. 1 <= x*x to 

992 # sqrt(1) <= x without any mention of cones. 

993 # if self._sf[0] is self._sf[1]: 

994 # # NOTE: The following is not allowed by SOCConstraint 

995 # # since this is effectively an affine constraint. 

996 # # return SOCConstraint(root, self._sf[0]) 

997 # 

998 # return root <= self._sf[0] 

999 # else: 

1000 # return RSOCConstraint(root, *self._sf) 

1001 else: 

1002 return NonconvexQuadraticConstraint(self, GE, other) 

1003 else: 

1004 return NotImplemented 

1005 

1006 

1007# -------------------------------------- 

1008__all__ = api_end(_API_START, globals())