Coverage for picos/expressions/exp

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

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

5# This file is part of PICOS.

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

20"""Implements :class:`Norm`."""

22import operator

23from collections import namedtuple

25import cvxopt

26import numpy

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

41_API_START = api_start(globals())

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

45class Norm(Expression):

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

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

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

54 definitions below to learn more.

56 :Definition:

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

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

61 absolute value

63 .. math::

65 |x| = \begin{cases}

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

67 -x, &\text{otherwise.}

68 \end{cases}

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

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

73 modulus

75 .. math::

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

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

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

84 .. math::

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

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

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

96 .. warning::

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

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

234 statements