Coverage for picos/expressions/exp_quadratic.py: 85.22%

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(x.scalar 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 @classmethod 

292 def _add_sub(cls, self, other, add, forward): 

293 def affine_part_and_string(other_aff): 

294 if add: 

295 if forward: 

296 affine = self._aff + other_aff 

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

298 else: 

299 affine = other_aff + self._aff 

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

301 else: 

302 if forward: 

303 affine = self._aff - other_aff 

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

305 else: 

306 affine = other_aff - self._aff 

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

308 

309 return affine, string 

310 

311 if isinstance(other, cls): 

312 affine, string = affine_part_and_string(other._aff) 

313 

314 quads = {} 

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

316 for vars in varPairs: 

317 if vars in self._quads: 

318 if vars in other._quads: 

319 if add: 

320 quads[vars] = self._quads[vars] + other._quads[vars] 

321 elif forward: 

322 quads[vars] = self._quads[vars] - other._quads[vars] 

323 else: 

324 quads[vars] = other._quads[vars] - self._quads[vars] 

325 elif not add and not forward: 

326 quads[vars] = -self._quads[vars] 

327 else: 

328 quads[vars] = self._quads[vars] 

329 elif not add and forward: 

330 quads[vars] = -other._quads[vars] 

331 else: 

332 quads[vars] = other._quads[vars] 

333 

334 return cls(string, quads, affine) 

335 elif isinstance(other, AffineExpression): 

336 affine, string = affine_part_and_string(other) 

337 

338 if add or forward: 

339 quads = self._quads 

340 copy_ = self 

341 else: 

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

343 copy_ = None 

344 

345 return cls(string, quads, affine, copyDecomposition=copy_) 

346 

347 if add: 

348 if forward: 

349 return Expression.__add__(self, other) 

350 else: 

351 return Expression.__radd__(self, other) 

352 else: 

353 if forward: 

354 return Expression.__sub__(self, other) 

355 else: 

356 return Expression.__rsub__(self, other) 

357 

358 @convert_operands(sameShape=True) 

359 def __add__(self, other): 

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

361 return QuadraticExpression._add_sub(self, other, True, True) 

362 

363 @convert_operands(sameShape=True) 

364 def __radd__(self, other): 

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

366 return QuadraticExpression._add_sub(self, other, True, False) 

367 

368 @convert_operands(sameShape=True) 

369 def __sub__(self, other): 

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

371 return QuadraticExpression._add_sub(self, other, False, True) 

372 

373 @convert_operands(sameShape=True) 

374 def __rsub__(self, other): 

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

376 return QuadraticExpression._add_sub(self, other, False, False) 

377 

378 @cached_selfinverse_unary_operator 

379 def __neg__(self): 

380 string = glyphs.clever_neg(self.string) 

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

382 affine = -self._aff 

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

384 

385 return QuadraticExpression(string, quads, affine, factors) 

386 

387 @classmethod 

388 def _mul_div(cls, self, other, div, forward): 

389 assert not div or forward 

390 

391 if isinstance(other, AffineExpression): 

392 if not other.constant: 

393 if div: 

394 # NIE as this would makes sense if dividing by a factor. 

395 raise NotImplementedError("You may only divide a quadratic " 

396 "expression by a constant term.") 

397 else: 

398 raise TypeError("You may only multiply a quadratic " 

399 "expression with a constant term.") 

400 

401 factor = other.safe_value 

402 

403 if not factor: 

404 if div: 

405 raise ZeroDivisionError( 

406 "Cannot divide {} by zero.".format(self.string)) 

407 else: 

408 return AffineExpression.zero() 

409 elif factor == 1: 

410 return self 

411 elif factor == -1: 

412 return -self 

413 

414 if div: 

415 factor = 1 / factor 

416 affine = self._aff / other 

417 string = glyphs.div(self.string, other.string) 

418 elif forward: 

419 affine = self._aff * other 

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

421 else: 

422 affine = other * self._aff 

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

424 

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

426 

427 if self._sf: 

428 r = factor**0.5 

429 if r.imag: 

430 factors = (r.imag*self._sf[0], -r.imag*self._sf[1]) 

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

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

433 else: 

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

435 else: 

436 factors = None 

437 

438 return cls(string, quads, affine, factors) 

439 

440 if div: 

441 return Expression.__truediv__(self, other) 

442 elif forward: 

443 return Expression.__mul__(self, other) 

444 else: 

445 return Expression.__rmul__(self, other) 

446 

447 @convert_operands(scalarRHS=True) 

448 def __mul__(self, other): 

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

450 return QuadraticExpression._mul_div(self, other, False, True) 

451 

452 @convert_operands(scalarRHS=True) 

453 def __rmul__(self, other): 

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

455 return QuadraticExpression._mul_div(self, other, False, False) 

456 

457 @convert_operands(scalarRHS=True) 

458 def __truediv__(self, other): 

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

460 return QuadraticExpression._mul_div(self, other, True, True) 

461 

462 # -------------------------------------------------------------------------- 

463 # Decomposition related methods. 

464 # -------------------------------------------------------------------------- 

465 

466 @cached_property 

467 def _sparse_quads(self): 

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

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

470 

471 def _Q_and_x_or_A_and_y(self, augmentedMatrix): 

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

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

474 # vectorization in the quadratic part. 

475 nonzero = {} 

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

477 nonzero.setdefault(x, set()) 

478 nonzero.setdefault(y, set()) 

479 

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

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

482 

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

484 if augmentedMatrix: 

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

486 nonzero.setdefault(x, set()) 

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

488 

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

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

491 

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

493 offset, offsets = 0, {} 

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

495 offsets[var] = offset 

496 offset += len(nz) 

497 

498 # Compute Q. 

499 V, I, J = [], [], [] 

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

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

502 

503 V.extend(B.V) 

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

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

506 

507 # Turn Q into A if requested. 

508 if augmentedMatrix: 

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

510 b = a[nonzero[x]] 

511 

512 V.extend(b.V) 

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

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

515 

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

517 I.append(offset) 

518 J.append(offset) 

519 

520 offset += 1 

521 

522 # Finalize Q/A. 

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

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

525 

526 if not Q: 

527 if augmentedMatrix: 

528 raise ValueError( 

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

530 else: 

531 raise ValueError( 

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

533 

534 # Compute x. 

535 x = None 

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

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

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

539 y = AffineExpression( 

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

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

542 

543 # Turn x into y if requested. 

544 if augmentedMatrix: 

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

546 

547 return Q, x 

548 

549 @cached_property 

550 def _Q_and_x(self): 

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

552 return self._Q_and_x_or_A_and_y(augmentedMatrix=False) 

553 

554 @cached_property 

555 def _A_and_y(self): 

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

557 return self._Q_and_x_or_A_and_y(augmentedMatrix=True) 

558 

559 @property 

560 def Q(self): 

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

562 

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

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

565 

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

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

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

569 

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

571 """ 

572 return self._Q_and_x[0] 

573 

574 @property 

575 def x(self): 

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

577 

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

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

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

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

582 

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

584 """ 

585 return self._Q_and_x[1] 

586 

587 @property 

588 def A(self): 

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

590 

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

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

593 but with zero rows and columns removed. 

594 

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

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

597 

598 :raises ValueError: When the expression is zero. 

599 """ 

600 return self._A_and_y[0] 

601 

602 @property 

603 def y(self): 

604 """See :attr:`A`. 

605 

606 :raises ValueError: When the expression is zero. 

607 """ 

608 return self._A_and_y[1] 

609 

610 def _L_or_M(self, augmentedMatrix): 

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

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

613 n = Q.size[0] 

614 

615 # Define a sequence of small numbers for perturbation. 

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

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

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

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

620 epsilons = set([0.0]) 

621 for i in range(PSD_PERTURBATIONS): 

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

623 epsilons = sorted(epsilons) 

624 

625 # Check if Q is really sparse. 

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

627 

628 # Attempt to find L. 

629 L = None 

630 if sparse: # Use CHOLMOD. 

631 F = None 

632 for eps in epsilons: 

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

634 

635 if F is None: 

636 F = cvxopt.cholmod.symbolic(P) 

637 

638 try: 

639 cvxopt.cholmod.numeric(P, F) 

640 except ArithmeticError: 

641 pass 

642 else: 

643 L = cvxopt.cholmod.getfactor(F) 

644 break 

645 else: # Use NumPy. 

646 Q = cvx2np(Q) 

647 for eps in epsilons: 

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

649 

650 try: 

651 L = numpy.linalg.cholesky(P) 

652 except numpy.linalg.LinAlgError: 

653 pass 

654 else: 

655 break 

656 

657 if L is None: 

658 if not PSD_PERTURBATIONS: 

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

660 "quadratic form that is either not positive semidefinite " 

661 "or singular. Try enabling PSD_PERTURBATIONS.") 

662 elif augmentedMatrix: 

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

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

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

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

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

668 )))) 

669 else: 

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

671 "quadratic part is not positive semidefinite." 

672 .format(self._symbStr)) 

673 

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

675 cutoff = 2 * n * maxEps**0.5 

676 if sparse: 

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

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

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

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

681 else: 

682 mask = abs(L) <= cutoff 

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

684 L[mask] = 0 

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

686 

687 return L 

688 

689 @cached_property 

690 def L(self): 

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

692 

693 :returns: A CVXOPT lower triangular sparse matrix. 

694 

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

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

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

698 """ 

699 return self._L_or_M(augmentedMatrix=False) 

700 

701 @cached_property 

702 def M(self): 

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

704 

705 :returns: A CVXOPT lower triangular sparse matrix. 

706 

707 :raises ValueError: When the expression is zero. 

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

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

710 numerically positive semidefinite. 

711 """ 

712 return self._L_or_M(augmentedMatrix=True) 

713 

714 @cached_property 

715 def quadroot(self): 

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

717 

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

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

720 

721 :Construction: 

722 

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

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

725 Then, 

726 

727 .. math:: 

728 x^TQx 

729 &= x^TLL^Tx \\ 

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

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

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

733 

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

735 

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

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

738 the quadratic part is not numerically positive semidefinite. 

739 """ 

740 L = self.L 

741 n = L.size[0] 

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

743 nnz = len(nz) 

744 

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

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

747 

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

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

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

751 return result 

752 

753 @cached_property 

754 def fullroot(self): 

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

756 

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

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

759 

760 :Construction: 

761 

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

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

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

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

766 

767 .. math:: 

768 x^TQx + a^Tx + b 

769 &= y^TAy \\ 

770 &= y^TMM^Ty \\ 

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

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

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

774 

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

776 

777 :raises ValueError: When the expression is zero. 

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

779 the quadratic part is not numerically positive semidefinite. 

780 """ 

781 M = self.M 

782 n = M.size[0] 

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

784 nnz = len(nz) 

785 

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

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

788 

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

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

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

792 return result 

793 

794 @property 

795 def rank(self): 

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

797 

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

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

800 

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

802 the quadratic part is not numerically positive semidefinite. 

803 """ 

804 try: 

805 return len(self.quadroot) 

806 except ValueError: 

807 return 0 # No quadratic part. 

808 

809 @cached_property 

810 def num_quad_terms(self): 

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

812 num_terms = 0 

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

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

815 if i >= j: 

816 num_terms += 1 

817 return num_terms 

818 

819 @property 

820 def is_squared_norm(self): 

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

822 

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

824 :math:`c` such that 

825 

826 .. math:: 

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

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

829 

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

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

832 removed. 

833 """ 

834 try: 

835 self.M