Coverage for picos/expressions/exp_norm.py: 88.46%
234 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-02-15 14:21 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2023-02-15 14:21 +0000
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# ------------------------------------------------------------------------------
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.
102 Special cases:
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}}`.
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).
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`.
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
124 .. math::
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}}
131 for :math:`p, q \in \mathbb{Q} \cup \{\infty\}` with :math:`p,q \geq 1`.
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).
139 Special cases:
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
146 .. math::
148 \lVert X \rVert_{1,\infty} =
149 \max_{j=1 \ldots n} \sum_{i=1}^m |X_{ij}|.
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`.
156 7. Complex matrix norms are not supported.
158 .. note::
160 You can write :math:`\infty` in Python as ``float("inf")``.
161 """
163 # --------------------------------------------------------------------------
164 # Initialization and factory methods.
165 # --------------------------------------------------------------------------
167 @convert_and_refine_arguments("x")
168 def __init__(self, x, p=2, q=None, denominator_limit=1000):
169 """Construct a :class:`Norm`.
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__))
187 complex = not isinstance(x, AffineExpression)
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
194 if q is None:
195 q = p
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)
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)
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")))
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)
258 if complex:
259 typeStr = "Complex " + typeStr
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
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
278 assert p == q
279 p, pNum, pDen = 2.0, 2, 1
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
286 assert isinstance(x, AffineExpression)
287 assert float(pNum) / float(pDen) == p and float(qNum) / float(qDen) == q
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
296 Expression.__init__(self, typeStr, symbStr)
298 # --------------------------------------------------------------------------
299 # Abstract method implementations and method overridings, except _predict.
300 # --------------------------------------------------------------------------
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
309 Subtype = namedtuple("Subtype", ("xShape", "pNum", "pDen", "qNum", "qDen"))
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)
317 def _get_value(self):
318 value = self._x._get_value()
320 if value.size == (1, 1):
321 return abs(value)
323 value = cvx2np(value)
325 p, q = self.p, self.q
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)
335 return cvxopt.matrix(value)
337 def _get_mutables(self):
338 return self._x._get_mutables()
340 def _is_convex(self):
341 return self.p >= 1 or len(self._x) == 1
343 def _is_concave(self):
344 return self.p <= 1 and len(self._x) != 1
346 def _replace_mutables(self, mapping):
347 return self.__class__(
348 self._x._replace_mutables(mapping), self.p, self.q, self._limit)
350 def _freeze_mutables(self, freeze):
351 return self.__class__(
352 self._x._freeze_mutables(freeze), self.p, self.q, self._limit)
354 # --------------------------------------------------------------------------
355 # Python special method implementations, except constraint-creating ones.
356 # --------------------------------------------------------------------------
358 @classmethod
359 def _mul(cls, self, other, forward):
360 if isinstance(other, AffineExpression) and other.constant:
361 factor = other.safe_value
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)
373 norm = cls(other*self._x, self.p, self.q, self._limit)
374 norm._typeStr = "Scaled " + norm._typeStr
375 norm._symbStr = string
377 return norm
379 if forward:
380 return Expression.__mul__(self, other)
381 else:
382 return Expression.__rmul__(self, other)
384 @convert_operands(scalarRHS=True)
385 @refine_operands()
386 def __mul__(self, other):
387 return Norm._mul(self, other, True)
389 @convert_operands(scalarRHS=True)
390 @refine_operands()
391 def __rmul__(self, other):
392 return Norm._mul(self, other, False)
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.")
402 p, q = self.p, self.q
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.")
412 return Expression.__pow__(self, other)
414 # --------------------------------------------------------------------------
415 # Properties and functions that describe the expression.
416 # --------------------------------------------------------------------------
418 @property
419 def p(self):
420 """The parameter :math:`p`.
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)
427 @property
428 def pnum(self):
429 """The limited precision fraction numerator of :math:`p`."""
430 return self._pNum
432 @property
433 def pden(self):
434 """The limited precision fraction denominator of :math:`p`."""
435 return self._pDen
437 @property
438 def q(self):
439 """The parameter :math:`q`.
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)
446 @property
447 def qnum(self):
448 """The limited precision fraction numerator of :math:`q`."""
449 return self._qNum
451 @property
452 def qden(self):
453 """The limited precision fraction denominator of :math:`q`."""
454 return self._qDen
456 # --------------------------------------------------------------------------
457 # Methods and properties that return modified copies.
458 # --------------------------------------------------------------------------
460 @property
461 def x(self):
462 """Real expression whose norm equals that of the original expression."""
463 return self._x
465 # --------------------------------------------------------------------------
466 # Constraint-creating operators, and _predict.
467 # --------------------------------------------------------------------------
469 @classmethod
470 def _predict(cls, subtype, relation, other):
471 assert isinstance(subtype, cls.Subtype)
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)
478 if relation == operator.__le__:
479 if not (xLen == 1 or p >= 1):
480 return NotImplemented # Not convex.
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.
498 assert p == q
500 if issubclass(other.clstype, AffineExpression) \
501 and other.subtype.dim == 1:
502 return VectorNormConstraint.make_type(
503 xLen, pNum, pDen, Constraint.GE)
505 return NotImplemented
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.")
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
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.")
537 assert self.p == self.q
539 if isinstance(other, AffineExpression):
540 return VectorNormConstraint(self, Constraint.GE, other)
541 else:
542 return NotImplemented
545# --------------------------------------
546__all__ = api_end(_API_START, globals())