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, 2020 Guillaume Sagnol 

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:`Norm`.""" 

21 

22import operator 

23from collections import namedtuple 

24 

25import cvxopt 

26import numpy 

27 

28from .. import glyphs 

29from ..apidoc import api_end, api_start 

30from ..caching import cached_unary_operator 

31from ..constraints import (AbsoluteValueConstraint, Constraint, 

32 MatrixNormConstraint, NuclearNormConstraint, 

33 SOCConstraint, SpectralNormConstraint, 

34 VectorNormConstraint) 

35from ..legacy import throw_deprecation_warning 

36from .data import (convert_and_refine_arguments, convert_operands, cvx2np, 

37 make_fraction) 

38from .exp_affine import AffineExpression, ComplexAffineExpression 

39from .exp_sqnorm import SquaredNorm 

40from .expression import Expression, refine_operands, validate_prediction 

41 

42_API_START = api_start(globals()) 

43# ------------------------------- 

44 

45 

46class Norm(Expression): 

47 r"""Entrywise :math:`p`-norm or :math:`L_{p,q}`-norm of an expression. 

48 

49 This class can represent the absolute value, the modulus, a vector 

50 :math:`p`-norm and an entrywise matrix :math:`L_{p,q}`-norm of a (complex) 

51 affine expression. In addition to these convex norms, it can represent the 

52 concave *generalized vector* :math:`p`-*norm* with :math:`0 < p \leq 1`. 

53 

54 Not all of these norms are available on the complex field; see the 

55 definitions below to learn more. 

56 

57 :Definition: 

58 

59 If :math:`q` is not given (:obj:`None`), then it is set equal to :math:`p`. 

60 

61 1. If the normed expression is a real scalar :math:`x`, then this is the 

62 absolute value 

63 

64 .. math:: 

65 

66 |x| = \begin{cases} 

67 x, &\text{if }x \geq 0,\\ 

68 -x, &\text{otherwise.} 

69 \end{cases} 

70 

71 The parameters :math:`p` and :math:`q` are ignored in this case. 

72 

73 2. If the normed expression is a complex scalar :math:`z`, then this is the 

74 modulus 

75 

76 .. math:: 

77 

78 |z| = \sqrt{\operatorname{Re}(z)^2 + \operatorname{Im}(z)^2}. 

79 

80 The parameters :math:`p` and :math:`q` are ignored in this case. 

81 

82 3. If the normed expression is a real vector :math:`x` and :math:`p = q`, 

83 then this is the (generalized) vector :math:`p`-norm 

84 

85 .. math:: 

86 

87 \lVert x \rVert_p = \left(\sum_{i=1}^n |x_i|^p\right)^{\frac{1}{p}} 

88 

89 for :math:`p \in \mathbb{Q} \cup \{\infty\}` with :math:`p > 0` and 

90 :math:`|x_i|` the absolute value of the :math:`i`-th entry of :math:`x`. 

91 

92 Note that for :math:`p < 1` the expression is not convex and thus not a 

93 proper norm. However, it is concave over the nonnegative orthant and 

94 posing a lower bound on such a generalized norm yields a convex 

95 constraint for :math:`x \geq 0`. 

96 

97 .. warning:: 

98 

99 When you pose a lower bound on a concave generalized norm 

100 (:math:`p < 1`), then PICOS enforces :math:`x \geq 0` through an 

101 auxiliary constraint during solution search. 

102 

103 Special cases: 

104 

105 - For :math:`p = 1`, this is the *Manhattan* or *Taxicab* norm 

106 :math:`\lVert x \rVert_{\text{sum}}`. 

107 - For :math:`p = 2`, this is the Euclidean norm 

108 :math:`\lVert x \rVert = \lVert x \rVert_2`. 

109 - For :math:`p = \infty`, this is the *Maximum*, *Chebyshev*, or 

110 *Infinity* norm :math:`\lVert x \rVert_{\text{max}}`. 

111 

112 4. If the normed expression is a real vector :math:`x` and 

113 :math:`p \neq q`, then it is treated as a matrix with a single row or a 

114 single column, depending on the shape associated with :math:`x`. See 

115 case (5). 

116 

117 5. If the normed expression is a complex vector :math:`z` and 

118 :math:`p = q`, then the definition is the same as in case (3) but with 

119 :math:`x = z`, :math:`|x_i|` the modulus of the :math:`i`-th entry of 

120 :math:`x`, and :math:`p \geq 1`. 

121 

122 6. If the normed expression is a real :math:`m \times n` matrix :math:`X`, 

123 then this is the :math:`L_{p,q}`-norm 

124 

125 .. math:: 

126 

127 \lVert X \rVert_{p,q} = 

128 \left(\sum_{j = 1}^n \left(\sum_{i = 1}^m 

129 |X_{ij}|^p 

130 \right)^{\frac{q}{p}} \right)^{\frac{1}{q}} 

131 

132 for :math:`p, q \in \mathbb{Q} \cup \{\infty\}` with :math:`p,q \geq 1`. 

133 

134 If :math:`p = q`, then this is equal to the (generalized) vector 

135 :math:`p`-norm of the the vectorized matrix, that is 

136 :math:`\lVert X \rVert_{p,p} = \lVert \operatorname{vec}(X) \rVert_p`. 

137 In this case, the requirement :math:`p \geq 1` is relaxed to 

138 :math:`p > 0` and :math:`X` may be a complex matrix. See case (3). 

139 

140 Special cases: 

141 

142 - For :math:`p = q = 2`, this is the Frobenius norm 

143 :math:`\lVert X \rVert = \lVert X \rVert_F`. 

144 - For :math:`p = 1`, :math:`q = \infty`, this is the maximum absolute 

145 column sum 

146 

147 .. math:: 

148 

149 \lVert X \rVert_{1,\infty} = 

150 \max_{j=1 \ldots n} \sum_{i=1}^m |X_{ij}|. 

151 

152 This equals the operator norm induced by the vector :math:`1`-norm. 

153 You can obtain the maximum absolute row sum (the operator norm 

154 induced by the vector :math:`\infty`-norm) by first transposing 

155 :math:`X`. 

156 

157 7. Complex matrix norms are not supported. 

158 

159 .. note:: 

160 

161 You can write :math:`\infty` in Python as ``float("inf")``. 

162 """ 

163 

164 # -------------------------------------------------------------------------- 

165 # Initialization and factory methods. 

166 # -------------------------------------------------------------------------- 

167 

168 @convert_and_refine_arguments("x") 

169 def __init__(self, x, p=2, q=None, denominator_limit=1000): 

170 """Construct a :class:`Norm`. 

171 

172 :param x: The affine expression to take the norm of. 

173 :type x: ~picos.expressions.ComplexAffineExpression 

174 :param float p: The value for :math:`p`, which is cast to a limited 

175 precision fraction. 

176 :param float q: The value for :math:`q`, which is cast to a limited 

177 precision fraction. The default of :obj:`None` means *equal to* 

178 :math:`p`. 

179 :param int denominator_limit: The largest allowed denominator when 

180 casting :math:`p` and :math:`q` to a fraction. Higher values can 

181 yield a greater precision at reduced performance. 

182 """ 

183 # Validate x. 

184 if not isinstance(x, ComplexAffineExpression): 

185 raise TypeError("Can only form the norm of an affine expression, " 

186 "not of {}.".format(type(x).__name__)) 

187 

188 complex = not isinstance(x, AffineExpression) 

189 

190 if isinstance(p, tuple) and len(p) == 2: 

191 throw_deprecation_warning("Arguments 'p' and 'q' to Norm must be " 

192 "given separately in the future.", decoratorLevel=1) 

193 p, q = p 

194 

195 if q is None: 

196 q = p 

197 

198 # Load p as a limtied precision fraction. 

199 if p == float("inf"): 

200 pNum = p 

201 pDen = 1 

202 pStr = glyphs.infty() 

203 else: 

204 pNum, pDen, p, pStr = make_fraction(p, denominator_limit) 

205 

206 # Load q as a limtied precision fraction. 

207 if q == float("inf"): 

208 qNum = q 

209 qDen = 1 

210 qStr = glyphs.infty() 

211 else: 

212 qNum, qDen, q, qStr = make_fraction(q, denominator_limit) 

213 

214 # Validate that p and q are in the allowed range. 

215 if p == q: 

216 if complex and p < 1: 

217 raise NotImplementedError( 

218 "Complex p-norm requires {}.".format(glyphs.ge("p", "1"))) 

219 elif p <= 0: 

220 raise ValueError( 

221 "p-norm requires {}.".format(glyphs.gt("p", "0"))) 

222 elif p != q: 

223 if complex: 

224 raise NotImplementedError( 

225 "(p,q)-norm is not supported for complex expressions.") 

226 elif p < 1 or q < 1: 

227 raise ValueError("(p,q)-norm requires {} and {}." 

228 .format(glyphs.ge("p", "1"), glyphs.ge("q", "1"))) 

229 

230 # Build the string representations. 

231 if len(x) == 1: 

232 typeStr = "Modulus" if complex else "Absolute Value" 

233 symbStr = glyphs.abs(x.string) 

234 elif p == q: 

235 vec = 1 in x.shape 

236 if p == 1: 

237 typeStr = "Manhattan Norm" 

238 symbStr = glyphs.pnorm(x.string, "sum") 

239 elif p == 2: 

240 typeStr = "Euclidean Norm" if vec else "Frobenius Norm" 

241 symbStr = glyphs.norm(x.string) 

242 elif p == float("inf"): 

243 typeStr = "Maximum Norm" 

244 symbStr = glyphs.pnorm(x.string, "max") 

245 else: 

246 if p < 1: 

247 typeStr = "Generalized p-Norm" if vec else \ 

248 "Entrywise Generalized p-Norm" 

249 else: 

250 typeStr = "Vector p-Norm" if vec else "Entrywise p-Norm" 

251 symbStr = glyphs.pnorm(x.string, pStr) 

252 else: 

253 if p == 1 and q == float("inf"): 

254 typeStr = "Maximum Absolute Column Sum" 

255 else: 

256 typeStr = "Matrix (p,q)-Norm" 

257 symbStr = glyphs.pqnorm(x.string, pStr, qStr) 

258 

259 if complex: 

260 typeStr = "Complex " + typeStr 

261 

262 # Reduce the complex to the real case. 

263 if complex: 

264 if len(x) == 1: 

265 # From modulus to real Euclidean norm. 

266 x = x.real // x.imag 

267 

268 p, pNum, pDen = 2.0, 2, 1 

269 q, qNum, qDen = 2.0, 2, 1 

270 else: 

271 # From complex vector/entrywise p-norm to real matrix p,q-norm. 

272 if x.shape[0] == 1: 

273 x = x.real // x.imag 

274 elif x.shape[1] == 1: 

275 x = x.T.real // x.T.imag 

276 else: 

277 x = x.vec.T.real // x.vec.T.imag 

278 

279 assert p == q 

280 p, pNum, pDen = 2.0, 2, 1 

281 

282 # Reduce an entrywise p-norm to a vector p-norm and make all vector 

283 # norms refer to column vectors. 

284 if p == q and x.shape[1] != 1: 

285 x = x.vec 

286 

287 assert isinstance(x, AffineExpression) 

288 assert float(pNum) / float(pDen) == p and float(qNum) / float(qDen) == q 

289 

290 self._x = x 

291 self._pNum = pNum 

292 self._pDen = pDen 

293 self._qNum = qNum 

294 self._qDen = qDen 

295 self._limit = denominator_limit 

296 

297 Expression.__init__(self, typeStr, symbStr) 

298 

299 # -------------------------------------------------------------------------- 

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

301 # -------------------------------------------------------------------------- 

302 

303 @cached_unary_operator 

304 def _get_refined(self): 

305 if self._x.constant: 

306 return AffineExpression.from_constant(self.value, 1, self.string) 

307 else: 

308 return self 

309 

310 Subtype = namedtuple("Subtype", ("xShape", "pNum", "pDen", "qNum", "qDen")) 

311 

312 def _get_subtype(self): 

313 # NOTE: The xShape field refers to the internal (real and potentially 

314 # vectorized) representation of the normed expression. 

315 return self.Subtype( 

316 self._x.shape, self._pNum, self._pDen, self._qNum, self._qDen) 

317 

318 def _get_value(self): 

319 value = self._x._get_value() 

320 

321 if value.size == (1, 1): 

322 return abs(value) 

323 

324 value = cvx2np(value) 

325 

326 p, q = self.p, self.q 

327 

328 if p == q: 

329 value = numpy.linalg.norm(numpy.ravel(value), p) 

330 else: 

331 columns = value.shape[1] 

332 value = numpy.linalg.norm([ 

333 numpy.linalg.norm(numpy.ravel(value[:, j]), p) 

334 for j in range(columns)], q) 

335 

336 return cvxopt.matrix(value) 

337 

338 def _get_mutables(self): 

339 return self._x._get_mutables() 

340 

341 def _is_convex(self): 

342 return self.p >= 1 or len(self._x) == 1 

343 

344 def _is_concave(self): 

345 return self.p <= 1 and len(self._x) != 1 

346 

347 def _replace_mutables(self, mapping): 

348 return self.__class__( 

349 self._x._replace_mutables(mapping), self.p, self.q, self._limit) 

350 

351 def _freeze_mutables(self, freeze): 

352 return self.__class__( 

353 self._x._freeze_mutables(freeze), self.p, self.q, self._limit) 

354 

355 # -------------------------------------------------------------------------- 

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

357 # -------------------------------------------------------------------------- 

358 

359 @convert_operands(scalarRHS=True) 

360 @refine_operands() 

361 def __mul__(self, other): 

362 if isinstance(other, AffineExpression): 

363 if not other.constant: 

364 raise NotImplementedError("You may only multiply a nonconstant " 

365 "PICOS norm with a constant term.") 

366 

367 if other.value < 0: 

368 raise NotImplementedError("You may only multiply a nonconstant " 

369 "PICOS norm with a nonnegative term.") 

370 

371 norm = Norm(other*self._x, self.p, self.q, self._limit) 

372 norm._typeStr = "Scaled " + norm._typeStr 

373 norm._symbStr = glyphs.clever_mul(self.string, other.string) 

374 return norm 

375 else: 

376 return NotImplemented 

377 

378 @convert_operands(scalarRHS=True) 

379 @refine_operands() 

380 def __rmul__(self, other): 

381 if isinstance(other, AffineExpression): 

382 norm = self.__mul__(other) 

383 # NOTE: __mul__ always creates a fresh expression. 

384 norm._symbStr = glyphs.clever_mul(other.string, self.string) 

385 return norm 

386 else: 

387 return NotImplemented 

388 

389 @convert_operands(scalarRHS=True) 

390 @refine_operands() 

391 def __pow__(self, other): 

392 if isinstance(other, AffineExpression): 

393 if not other.constant or other.value != 2: 

394 raise NotImplementedError("You may only take a PICOS norm to " 

395 "the power of two.") 

396 

397 p, q = self.p, self.q 

398 

399 if p == q and p == 2: 

400 result = SquaredNorm(self._x) 

401 result._symbStr = glyphs.squared(self.string) 

402 return result 

403 else: 

404 raise NotImplementedError( 

405 "You may only square an Euclidean or Frobenius norm.") 

406 else: 

407 return NotImplemented 

408 

409 # -------------------------------------------------------------------------- 

410 # Properties and functions that describe the expression. 

411 # -------------------------------------------------------------------------- 

412 

413 @property 

414 def p(self): 

415 """The parameter :math:`p`. 

416 

417 This is a limited precision version of the parameter used when the norm 

418 was constructed. 

419 """ 

420 return float(self._pNum) / float(self._pDen) 

421 

422 @property 

423 def pnum(self): 

424 """The limited precision fraction numerator of :math:`p`.""" 

425 return self._pNum 

426 

427 @property 

428 def pden(self): 

429 """The limited precision fraction denominator of :math:`p`.""" 

430 return self._pDen 

431 

432 @property 

433 def q(self): 

434 """The parameter :math:`q`. 

435 

436 This is a limited precision version of the parameter used when the norm 

437 was constructed. 

438 """ 

439 return float(self._qNum) / float(self._qDen) 

440 

441 @property 

442 def qnum(self): 

443 """The limited precision fraction numerator of :math:`q`.""" 

444 return self._qNum 

445 

446 @property 

447 def qden(self): 

448 """The limited precision fraction denominator of :math:`q`.""" 

449 return self._qDen 

450 

451 # -------------------------------------------------------------------------- 

452 # Methods and properties that return modified copies. 

453 # -------------------------------------------------------------------------- 

454 

455 @property 

456 def x(self): 

457 """Real expression whose norm equals that of the original expression.""" 

458 return self._x 

459 

460 # -------------------------------------------------------------------------- 

461 # Constraint-creating operators, and _predict. 

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

463 

464 @classmethod 

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

466 assert isinstance(subtype, cls.Subtype) 

467 

468 xShape, pNum, pDen, qNum, qDen = subtype 

469 xLen = xShape[0] * xShape[1] 

470 p = float(pNum) / float(pDen) 

471 q = float(qNum) / float(qDen) 

472 

473 if relation == operator.__le__: 

474 if not (xLen == 1 or p >= 1): 

475 return NotImplemented # Not convex. 

476 

477 if issubclass(other.clstype, AffineExpression) \ 

478 and other.subtype.dim == 1: 

479 if xLen == 1: 

480 return AbsoluteValueConstraint.make_type() 

481 elif p == q == 2: 

482 return SOCConstraint.make_type(argdim=xLen) 

483 elif p == q: 

484 return VectorNormConstraint.make_type( 

485 xLen, pNum, pDen, Constraint.LE) 

486 else: 

487 return MatrixNormConstraint.make_type( 

488 xShape, pNum, pDen, qNum, qDen) 

489 elif relation == operator.__ge__: 

490 if not (xLen != 1 and p <= 1): 

491 return NotImplemented # Not concave. 

492 

493 assert p == q 

494 

495 if issubclass(other.clstype, AffineExpression) \ 

496 and other.subtype.dim == 1: 

497 return VectorNormConstraint.make_type( 

498 xLen, pNum, pDen, Constraint.GE) 

499 

500 return NotImplemented 

501 

502 @convert_operands(scalarRHS=True) 

503 @validate_prediction 

504 @refine_operands() 

505 def __le__(self, other): 

506 if not self.convex: 

507 raise TypeError("Cannot upper-bound a nonconvex generalized norm.") 

508 

509 if isinstance(other, AffineExpression): 

510 if len(self._x) == 1: 

511 return AbsoluteValueConstraint(self._x, other) 

512 elif self.p == self.q == 2: 

513 # NOTE: The custom string is necessary in case the original x 

514 # was complex; otherwise the string would be something 

515 # like "‖vec([Re(xᵀ); Im(xᵀ)])‖ ≤ b". 

516 return SOCConstraint(self._x, other, 

517 customString=glyphs.le(self.string, other.string)) 

518 elif self.p == self.q: 

519 return VectorNormConstraint(self, Constraint.LE, other) 

520 else: 

521 return MatrixNormConstraint(self, other) 

522 else: 

523 return NotImplemented 

524 

525 @convert_operands(scalarRHS=True) 

526 @validate_prediction 

527 @refine_operands() 

528 def __ge__(self, other): 

529 if not self.concave: 

530 raise TypeError("Cannot lower-bound a nonconcave norm.") 

531 

532 assert self.p == self.q 

533 

534 if isinstance(other, AffineExpression): 

535 return VectorNormConstraint(self, Constraint.GE, other) 

536 else: 

537 return NotImplemented 

538 

539 

540class SpectralNorm(Expression): 

541 r"""The spectral norm of a matrix. 

542 

543 This class can represent the spectral norm of a matrix-affine expression 

544 (real- or complex valued). The spectral norm is convex, so we can form 

545 expressions of the form ``SpectralNorm(X) <= t`` which are typically 

546 reformulated as LMIs that can be handled by SDP solvers. 

547 

548 :Definition: 

549 

550 If the normed expression is a matrix :math:`X`, then its spectral norm is 

551 

552 .. math:: 

553 

554 \|X\|_2 = \max \{ \|Xu\|_2 : \|u\| \leq 1\} 

555 = \sqrt{\lambda_{\max}(XX^*)}, 

556 

557 where :math:`\lambda_{\max}(\cdot)` denotes the largest eigenvalue of 

558 a matrix, and :math:`X^*` denotes the adjoint matrix of :math:`X` 

559 (i.e., the transposed matrix :math:`X^T` if :math:`X` is real-valued). 

560 

561 Special cases: 

562 

563 - If :math:`X` is scalar, then :math:`\|X\|_2` reduces to the the absolute 

564 value (or modulus) :math:`|X|`. 

565 - If :math:`X` is scalar, then :math:`\|X\|_2` coincides with the 

566 Euclidean norm of :math:`X`. 

567 

568 """ 

569 

570 @convert_and_refine_arguments("x") 

571 def __init__(self, x): 

572 """Construct a :class:`SpectralNorm`. 

573 

574 :param x: The affine expression to take the norm of. 

575 :type x: ~picos.expressions.ComplexAffineExpression 

576 """ 

577 # Validate x. 

578 if not isinstance(x, ComplexAffineExpression): 

579 raise TypeError("Can only form the spectral norm of an affine " 

580 "expression, not of {}.".format(type(x).__name__)) 

581 

582 complex = not isinstance(x, AffineExpression) 

583 

584 # Build the string representations. 

585 if len(x) == 1: 

586 typeStr = "Modulus" if complex else "Absolute Value" 

587 symbStr = glyphs.abs(x.string) 

588 elif 1 in x.shape: 

589 typeStr = "Euclidean Norm" 

590 symbStr = glyphs.norm(x.string) 

591 else: 

592 typeStr = "Spectral Norm" 

593 symbStr = glyphs.spnorm(x.string) 

594 

595 if complex: 

596 typeStr = "Complex " + typeStr 

597 

598 self._x = x 

599 self._complex = complex 

600 Expression.__init__(self, typeStr, symbStr) 

601 

602 # -------------------------------------------------------------------------- 

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

604 # -------------------------------------------------------------------------- 

605 

606 @cached_unary_operator 

607 def _get_refined(self): 

608 if self._x.constant: 

609 return AffineExpression.from_constant(self.value, 1, self.string) 

610 elif len(self._x) == 1 or (1 in self._x.shape): 

611 return Norm(self._x) 

612 else: 

613 return self 

614 

615 Subtype = namedtuple("Subtype", ("argshape", "complex", "hermitian")) 

616 

617 def _get_subtype(self): 

618 return self.Subtype(self._x.shape, self._complex, self._x.hermitian) 

619 

620 def _get_value(self): 

621 value = self._x._get_value() 

622 value = cvx2np(value) 

623 value = numpy.linalg.norm(value, 2) 

624 return cvxopt.matrix(value) 

625 

626 def _get_mutables(self): 

627 return self._x._get_mutables() 

628 

629 def _is_convex(self): 

630 return True 

631 

632 def _is_concave(self): 

633 return False 

634 

635 def _replace_mutables(self, mapping): 

636 return self.__class__(self._x._replace_mutables(mapping)) 

637 

638 def _freeze_mutables(self, freeze): 

639 return self.__class__(self._x._freeze_mutables(freeze)) 

640 

641 # -------------------------------------------------------------------------- 

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

643 # -------------------------------------------------------------------------- 

644 

645 @convert_operands(scalarRHS=True) 

646 @refine_operands() 

647 def __mul__(self, other): 

648 if isinstance(other, AffineExpression): 

649 if not other.constant: 

650 raise NotImplementedError("You may only multiply a nonconstant " 

651 "PICOS norm with a constant term.") 

652 

653 if other.value < 0: 

654 raise NotImplementedError("You may only multiply a nonconstant " 

655 "PICOS norm with a nonnegative term.") 

656 

657 norm = SpectralNorm(other*self._x) 

658 norm._typeStr = "Scaled " + norm._typeStr 

659 norm._symbStr = glyphs.clever_mul(self.string, other.string) 

660 return norm 

661 else: 

662 return NotImplemented 

663 

664 @convert_operands(scalarRHS=True) 

665 @refine_operands() 

666 def __rmul__(self, other): 

667 if isinstance(other, AffineExpression): 

668 norm = self.__mul__(other) 

669 # NOTE: __mul__ always creates a fresh expression. 

670 norm._symbStr = glyphs.clever_mul(other.string, self.string) 

671 return norm 

672 else: 

673 return NotImplemented 

674 

675 # -------------------------------------------------------------------------- 

676 # Methods and properties that return modified copies. 

677 # -------------------------------------------------------------------------- 

678 

679 @property 

680 def x(self): 

681 """Real expression whose norm equals that of the original expression.""" 

682 return self._x 

683 

684 # -------------------------------------------------------------------------- 

685 # Constraint-creating operators, and _predict. 

686 # -------------------------------------------------------------------------- 

687 

688 @classmethod 

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

690 assert isinstance(subtype, cls.Subtype) 

691 

692 arg_shape, arg_complex, arg_hermitian = subtype 

693 xLen = arg_shape[0] * arg_shape[1] 

694 

695 if relation == operator.__le__: 

696 if issubclass(other.clstype, AffineExpression) \ 

697 and other.subtype.dim == 1: 

698 if xLen == 1: 

699 return AbsoluteValueConstraint.make_type() 

700 elif 1 in arg_shape: 

701 assert False, "Unexpected case (should have been refined)" 

702 else: 

703 return SpectralNormConstraint.make_type( 

704 arg_shape, arg_complex, arg_hermitian) 

705 elif relation == operator.__ge__: 

706 return NotImplemented # Not concave. 

707 

708 return NotImplemented 

709 

710 @convert_operands(scalarRHS=True) 

711 @validate_prediction 

712 @refine_operands() 

713 def __le__(self, other): 

714 assert self.convex 

715 

716 if isinstance(other, AffineExpression): 

717 if len(self._x) == 1: 

718 return AbsoluteValueConstraint(self._x, other) 

719 elif 1 in self._x.shape: 

720 assert False, "Unexpected case (should have been refined)" 

721 else: 

722 return SpectralNormConstraint(self, other) 

723 else: 

724 return NotImplemented 

725 

726 

727class NuclearNorm(Expression): 

728 r"""The nuclear norm of a matrix. 

729 

730 This class can represent the nuclear norm of a matrix-affine expression 

731 (real- or complex valued). The nuclear norm is convex, so we can form 

732 expressions of the form ``NuclearNorm(X) <= t`` which are typically 

733 reformulated as LMIs that can be handled by SDP solvers. 

734 

735 :Definition: 

736 

737 If the normed expression is a matrix :math:`X`, then its nuclear norm is 

738 

739 .. math:: 

740 

741 \|X\|_* = \operatorname{trace}\ (X^*X)^{1/2} 

742 = \sum_{i=1}^{\min(n,m)} \sigma_i(X) 

743 

744 where the :math:`\sigma_i(X)` denote the singular values of 

745 a :math:`X`, and :math:`X^*` denotes the adjoint matrix of :math:`X` 

746 (i.e., the transposed matrix :math:`X^T` if :math:`X` is real-valued). 

747 

748 Special cases: 

749 

750 - If :math:`X` is scalar, then :math:`\|X\|_*` reduces to the the absolute 

751 value (or modulus) :math:`|X|`. 

752 - If :math:`X` is scalar, then :math:`\|X\|_*` coincides with the 

753 Euclidean norm of :math:`X`. 

754 

755 """ 

756 

757 @convert_and_refine_arguments("x") 

758 def __init__(self, x): 

759 """Construct a :class:`NuclearNorm`. 

760 

761 :param x: The affine expression to take the norm of. 

762 :type x: ~picos.expressions.ComplexAffineExpression 

763 """ 

764 # Validate x. 

765 if not isinstance(x, ComplexAffineExpression): 

766 raise TypeError("Can only form the nuclear norm of an affine " 

767 "expression, not of {}.".format(type(x).__name__)) 

768 

769 complex = not isinstance(x, AffineExpression) 

770 

771 # Build the string representations. 

772 if len(x) == 1: 

773 typeStr = "Modulus" if complex else "Absolute Value" 

774 symbStr = glyphs.abs(x.string) 

775 elif 1 in x.shape: 

776 typeStr = "Euclidean Norm" 

777 symbStr = glyphs.norm(x.string) 

778 else: 

779 typeStr = "Nuclear Norm" 

780 symbStr = glyphs.ncnorm(x.string) 

781 

782 if complex: 

783 typeStr = "Complex " + typeStr 

784 

785 self._x = x 

786 self._complex = complex 

787 Expression.__init__(self, typeStr, symbStr) 

788 

789 # -------------------------------------------------------------------------- 

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

791 # -------------------------------------------------------------------------- 

792 

793 @cached_unary_operator 

794 def _get_refined(self): 

795 if self._x.constant: 

796 return AffineExpression.from_constant(self.value, 1, self.string) 

797 elif len(self._x) == 1 or (1 in self._x.shape): 

798 return Norm(self._x) 

799 else: 

800 return self 

801 

802 Subtype = namedtuple("Subtype", ("argshape", "complex", "hermitian")) 

803 

804 def _get_subtype(self): 

805 return self.Subtype(self._x.shape, self._complex, self._x.hermitian) 

806 

807 def _get_value(self): 

808 value = self._x._get_value() 

809 value = cvx2np(value) 

810 value = numpy.linalg.norm(value, 'nuc') 

811 return cvxopt.matrix(value) 

812 

813 def _get_mutables(self): 

814 return self._x._get_mutables() 

815 

816 def _is_convex(self): 

817 return True 

818 

819 def _is_concave(self): 

820 return False 

821 

822 def _replace_mutables(self, mapping): 

823 return self.__class__(self._x._replace_mutables(mapping)) 

824 

825 def _freeze_mutables(self, freeze): 

826 return self.__class__(self._x._freeze_mutables(freeze)) 

827 

828 # -------------------------------------------------------------------------- 

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

830 # -------------------------------------------------------------------------- 

831 

832 @convert_operands(scalarRHS=True) 

833 @refine_operands() 

834 def __mul__(self, other): 

835 if isinstance(other, AffineExpression): 

836 if not other.constant: 

837 raise NotImplementedError("You may only multiply a nonconstant " 

838 "PICOS norm with a constant term.") 

839 

840 if other.value < 0: 

841 raise NotImplementedError("You may only multiply a nonconstant " 

842 "PICOS norm with a nonnegative term.") 

843 

844 norm = NuclearNorm(other*self._x) 

845 norm._typeStr = "Scaled " + norm._typeStr 

846 norm._symbStr = glyphs.clever_mul(self.string, other.string) 

847 return norm 

848 else: 

849 return NotImplemented 

850 

851 @convert_operands(scalarRHS=True) 

852 @refine_operands() 

853 def __rmul__(self, other): 

854 if isinstance(other, AffineExpression): 

855 norm = self.__mul__(other) 

856 # NOTE: __mul__ always creates a fresh expression. 

857 norm._symbStr = glyphs.clever_mul(other.string, self.string) 

858 return norm 

859 else: 

860 return NotImplemented 

861 

862 # -------------------------------------------------------------------------- 

863 # Methods and properties that return modified copies. 

864 # -------------------------------------------------------------------------- 

865 

866 @property 

867 def x(self): 

868 """Real expression whose norm equals that of the original expression.""" 

869 return self._x 

870 

871 # -------------------------------------------------------------------------- 

872 # Constraint-creating operators, and _predict. 

873 # -------------------------------------------------------------------------- 

874 

875 @classmethod 

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

877 assert isinstance(subtype, cls.Subtype) 

878 

879 arg_shape, arg_complex, arg_hermitian = subtype 

880 xLen = arg_shape[0] * arg_shape[1] 

881 

882 if relation == operator.__le__: 

883 if issubclass(other.clstype, AffineExpression) \ 

884 and other.subtype.dim == 1: 

885 if xLen == 1: 

886 return AbsoluteValueConstraint.make_type() 

887 elif 1 in arg_shape: 

888 assert False, "Unexpected case (should have been refined)" 

889 else: 

890 return NuclearNormConstraint.make_type( 

891 arg_shape, arg_complex, arg_hermitian) 

892 elif relation == operator.__ge__: 

893 return NotImplemented # Not concave. 

894 

895 return NotImplemented 

896 

897 @convert_operands(scalarRHS=True) 

898 @validate_prediction 

899 @refine_operands() 

900 def __le__(self, other): 

901 assert self.convex 

902 

903 if isinstance(other, AffineExpression): 

904 if len(self._x) == 1: 

905 return AbsoluteValueConstraint(self._x, other) 

906 elif 1 in self._x.shape: 

907 assert False, "Unexpected case (should have been refined)" 

908 else: 

909 return NuclearNormConstraint(self, other) 

910 else: 

911 return NotImplemented 

912 

913 

914# -------------------------------------- 

915__all__ = api_end(_API_START, globals())