Coverage for picos/expressions/exp_norm.py: 88.46%

Shortcuts on this page

r m x   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

234 statements  

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:`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, SOCConstraint, 

33 VectorNormConstraint) 

34from ..legacy import throw_deprecation_warning 

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

36 make_fraction) 

37from .exp_affine import AffineExpression, ComplexAffineExpression 

38from .exp_sqnorm import SquaredNorm 

39from .expression import Expression, refine_operands, validate_prediction 

40 

41_API_START = api_start(globals()) 

42# ------------------------------- 

43 

44 

45class Norm(Expression): 

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

47 

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

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

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

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

52 

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

54 definitions below to learn more. 

55 

56 :Definition: 

57 

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

59 

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

61 absolute value 

62 

63 .. math:: 

64 

65 |x| = \begin{cases} 

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

67 -x, &\text{otherwise.} 

68 \end{cases} 

69 

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

71 

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

73 modulus 

74 

75 .. math:: 

76 

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

78 

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

80 

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

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

83 

84 .. math:: 

85 

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

87 

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

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

90 

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

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

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

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

95 

96 .. warning:: 

97 

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

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

100 auxiliary constraint during solution search. 

101 

102 Special cases: 

103 

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

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

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

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

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

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

110 

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

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

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

114 case (5). 

115 

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

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

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

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

120 

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

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

123 

124 .. math:: 

125 

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

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

128 |X_{ij}|^p 

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

130 

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

132 

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

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

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

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

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

138 

139 Special cases: 

140 

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

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

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

144 column sum 

145 

146 .. math:: 

147 

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

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

150 

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

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

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

154 :math:`X`. 

155 

156 7. Complex matrix norms are not supported. 

157 

158 .. note:: 

159 

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

161 """ 

162 

163 # -------------------------------------------------------------------------- 

164 # Initialization and factory methods. 

165 # -------------------------------------------------------------------------- 

166 

167 @convert_and_refine_arguments("x") 

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

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

170 

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

172 :type x: ~picos.expressions.ComplexAffineExpression 

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

174 precision fraction. 

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

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

177 :math:`p`. 

178 :param int denominator_limit: The largest allowed denominator when 

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

180 yield a greater precision at reduced performance. 

181 """ 

182 # Validate x. 

183 if not isinstance(x, ComplexAffineExpression): 

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

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

186 

187 complex = not isinstance(x, AffineExpression) 

188 

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

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

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

192 p, q = p 

193 

194 if q is None: 

195 q = p 

196 

197 # Load p as a limtied precision fraction. 

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

199 pNum = p 

200 pDen = 1 

201 pStr = glyphs.infty() 

202 else: 

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

204 

205 # Load q as a limtied precision fraction. 

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

207 qNum = q 

208 qDen = 1 

209 qStr = glyphs.infty() 

210 else: 

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

212 

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

214 if p == q: 

215 if complex and p < 1: 

216 raise NotImplementedError( 

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

218 elif p <= 0: 

219 raise ValueError( 

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

221 elif p != q: 

222 if complex: 

223 raise NotImplementedError( 

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

225 elif p < 1 or q < 1: 

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

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

228 

229 # Build the string representations. 

230 if len(x) == 1: 

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

232 symbStr = glyphs.abs(x.string) 

233 elif p == q: 

234 vec = 1 in x.shape 

235 if p == 1: 

236 typeStr = "Manhattan Norm" 

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

238 elif p == 2: 

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

240 symbStr = glyphs.norm(x.string) 

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

242 typeStr = "Maximum Norm" 

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

244 else: 

245 if p < 1: 

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

247 "Entrywise Generalized p-Norm" 

248 else: 

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

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

251 else: 

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

253 typeStr = "Maximum Absolute Column Sum" 

254 else: 

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

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

257 

258 if complex: 

259 typeStr = "Complex " + typeStr 

260 

261 # Reduce the complex to the real case. 

262 if complex: 

263 if len(x) == 1: 

264 # From modulus to real Euclidean norm. 

265 x = x.real // x.imag 

266 

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

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

269 else: 

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

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

272 x = x.real // x.imag 

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

274 x = x.T.real // x.T.imag 

275 else: 

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

277 

278 assert p == q 

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

280 

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

282 # norms refer to column vectors. 

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

284 x = x.vec 

285 

286 assert isinstance(x, AffineExpression) 

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

288 

289 self._x = x 

290 self._pNum = pNum 

291 self._pDen = pDen 

292 self._qNum = qNum 

293 self._qDen = qDen 

294 self._limit = denominator_limit 

295 

296 Expression.__init__(self, typeStr, symbStr) 

297 

298 # -------------------------------------------------------------------------- 

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

300 # -------------------------------------------------------------------------- 

301 

302 @cached_unary_operator 

303 def _get_refined(self): 

304 if self._x.constant: 

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

306 else: 

307 return self 

308 

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

310 

311 def _get_subtype(self): 

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

313 # vectorized) representation of the normed expression. 

314 return self.Subtype( 

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

316 

317 def _get_value(self): 

318 value = self._x._get_value() 

319 

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

321 return abs(value) 

322 

323 value = cvx2np(value) 

324 

325 p, q = self.p, self.q 

326 

327 if p == q: 

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

329 else: 

330 columns = value.shape[1] 

331 value = numpy.linalg.norm([ 

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

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

334 

335 return cvxopt.matrix(value) 

336 

337 def _get_mutables(self): 

338 return self._x._get_mutables() 

339 

340 def _is_convex(self): 

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

342 

343 def _is_concave(self): 

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

345 

346 def _replace_mutables(self, mapping): 

347 return self.__class__( 

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

349 

350 def _freeze_mutables(self, freeze): 

351 return self.__class__( 

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

353 

354 # -------------------------------------------------------------------------- 

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

356 # -------------------------------------------------------------------------- 

357 

358 @classmethod 

359 def _mul(cls, self, other, forward): 

360 if isinstance(other, AffineExpression) and other.constant: 

361 factor = other.safe_value 

362 

363 if not factor: 

364 return AffineExpression.zero() 

365 elif factor == 1: 

366 return self 

367 elif factor > 0: 

368 if forward: 

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

370 else: 

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

372 

373 norm = cls(other*self._x, self.p, self.q, self._limit) 

374 norm._typeStr = "Scaled " + norm._typeStr 

375 norm._symbStr = string 

376 

377 return norm 

378 

379 if forward: 

380 return Expression.__mul__(self, other) 

381 else: 

382 return Expression.__rmul__(self, other) 

383 

384 @convert_operands(scalarRHS=True) 

385 @refine_operands() 

386 def __mul__(self, other): 

387 return Norm._mul(self, other, True) 

388 

389 @convert_operands(scalarRHS=True) 

390 @refine_operands() 

391 def __rmul__(self, other): 

392 return Norm._mul(self, other, False) 

393 

394 @convert_operands(scalarRHS=True) 

395 @refine_operands() 

396 def __pow__(self, other): 

397 if isinstance(other, AffineExpression): 

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

399 raise NotImplementedError( 

400 "You may only take a norm to the power of two.") 

401 

402 p, q = self.p, self.q 

403 

404 if p == q and p == 2: 

405 result = SquaredNorm(self._x) 

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

407 return result 

408 else: 

409 raise NotImplementedError( 

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

411 

412 return Expression.__pow__(self, other) 

413 

414 # -------------------------------------------------------------------------- 

415 # Properties and functions that describe the expression. 

416 # -------------------------------------------------------------------------- 

417 

418 @property 

419 def p(self): 

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

421 

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

423 was constructed. 

424 """ 

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

426 

427 @property 

428 def pnum(self): 

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

430 return self._pNum 

431 

432 @property 

433 def pden(self): 

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

435 return self._pDen 

436 

437 @property 

438 def q(self): 

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

440 

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

442 was constructed. 

443 """ 

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

445 

446 @property 

447 def qnum(self): 

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

449 return self._qNum 

450 

451 @property 

452 def qden(self): 

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

454 return self._qDen 

455 

456 # -------------------------------------------------------------------------- 

457 # Methods and properties that return modified copies. 

458 # -------------------------------------------------------------------------- 

459 

460 @property 

461 def x(self): 

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

463 return self._x 

464 

465 # -------------------------------------------------------------------------- 

466 # Constraint-creating operators, and _predict. 

467 # -------------------------------------------------------------------------- 

468 

469 @classmethod 

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

471 assert isinstance(subtype, cls.Subtype) 

472 

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

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

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

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

477 

478 if relation == operator.__le__: 

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

480 return NotImplemented # Not convex. 

481 

482 if issubclass(other.clstype, AffineExpression) \ 

483 and other.subtype.dim == 1: 

484 if xLen == 1: 

485 return AbsoluteValueConstraint.make_type() 

486 elif p == q == 2: 

487 return SOCConstraint.make_type(argdim=xLen) 

488 elif p == q: 

489 return VectorNormConstraint.make_type( 

490 xLen, pNum, pDen, Constraint.LE) 

491 else: 

492 return MatrixNormConstraint.make_type( 

493 xShape, pNum, pDen, qNum, qDen) 

494 elif relation == operator.__ge__: 

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

496 return NotImplemented # Not concave. 

497 

498 assert p == q 

499 

500 if issubclass(other.clstype, AffineExpression) \ 

501 and other.subtype.dim == 1: 

502 return VectorNormConstraint.make_type( 

503 xLen, pNum, pDen, Constraint.GE) 

504 

505 return NotImplemented 

506 

507 @convert_operands(scalarRHS=True) 

508 @validate_prediction 

509 @refine_operands() 

510 def __le__(self, other): 

511 if not self.convex: 

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

513 

514 if isinstance(other, AffineExpression): 

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

516 return AbsoluteValueConstraint(self._x, other) 

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

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

519 # was complex; otherwise the string would be something 

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

521 return SOCConstraint(self._x, other, 

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

523 elif self.p == self.q: 

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

525 else: 

526 return MatrixNormConstraint(self, other) 

527 else: 

528 return NotImplemented 

529 

530 @convert_operands(scalarRHS=True) 

531 @validate_prediction 

532 @refine_operands() 

533 def __ge__(self, other): 

534 if not self.concave: 

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

536 

537 assert self.p == self.q 

538 

539 if isinstance(other, AffineExpression): 

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

541 else: 

542 return NotImplemented 

543 

544 

545# -------------------------------------- 

546__all__ = api_end(_API_START, globals())