Source code for picos.expressions.exp_norm

# ------------------------------------------------------------------------------
# Copyright (C) 2019 Maximilian Stahlberg
# Based on the original picos.expressions module by Guillaume Sagnol.
#
# This file is part of PICOS.
#
# PICOS is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along with
# this program.  If not, see <http://www.gnu.org/licenses/>.
# ------------------------------------------------------------------------------

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

import operator
from collections import namedtuple

import cvxopt
import numpy

from .. import glyphs
from ..apidoc import api_end, api_start
from ..caching import cached_unary_operator
from ..constraints import (AbsoluteValueConstraint, Constraint,
                           MatrixNormConstraint, SOCConstraint,
                           VectorNormConstraint)
from ..legacy import throw_deprecation_warning
from .data import (convert_and_refine_arguments, convert_operands, cvx2np,
                   make_fraction)
from .exp_affine import AffineExpression, ComplexAffineExpression
from .exp_sqnorm import SquaredNorm
from .expression import Expression, refine_operands, validate_prediction

_API_START = api_start(globals())
# -------------------------------


[docs]class Norm(Expression): r"""Entrywise :math:`p`-norm or :math:`L_{p,q}`-norm of an expression. This class can represent the absolute value, the modulus, a vector :math:`p`-norm and an entrywise matrix :math:`L_{p,q}`-norm of a (complex) affine expression. In addition to these convex norms, it can represent the concave *generalized vector* :math:`p`-*norm* with :math:`0 < p \leq 1`. Not all of these norms are available on the complex field; see the definitions below to learn more. :Definition: If :math:`q` is not given (:obj:`None`), then it is set equal to :math:`p`. 1. If the normed expression is a real scalar :math:`x`, then this is the absolute value .. math:: |x| = \begin{cases} x, &\text{if }x \geq 0,\\ -x, &\text{otherwise.} \end{cases} The parameters :math:`p` and :math:`q` are ignored in this case. 2. If the normed expression is a complex scalar :math:`z`, then this is the modulus .. math:: |z| = \sqrt{\operatorname{Re}(z)^2 + \operatorname{Im}(z)^2}. The parameters :math:`p` and :math:`q` are ignored in this case. 3. If the normed expression is a real vector :math:`x` and :math:`p = q`, then this is the (generalized) vector :math:`p`-norm .. math:: \lVert x \rVert_p = \left(\sum_{i=1}^n |x_i|^p\right)^{\frac{1}{p}} for :math:`p \in \mathbb{Q} \cup \{\infty\}` with :math:`p > 0` and :math:`|x_i|` the absolute value of the :math:`i`-th entry of :math:`x`. Note that for :math:`p < 1` the expression is not convex and thus not a proper norm. However, it is concave over the nonnegative orthant and posing a lower bound on such a generalized norm yields a convex constraint for :math:`x \geq 0`. .. warning:: When you pose a lower bound on a concave generalized norm (:math:`p < 1`), then PICOS enforces :math:`x \geq 0` through an auxiliary constraint during solution search. Special cases: - For :math:`p = 1`, this is the *Manhattan* or *Taxicab* norm :math:`\lVert x \rVert_{\text{sum}}`. - For :math:`p = 2`, this is the Euclidean norm :math:`\lVert x \rVert = \lVert x \rVert_2`. - For :math:`p = \infty`, this is the *Maximum*, *Chebyshev*, or *Infinity* norm :math:`\lVert x \rVert_{\text{max}}`. 4. If the normed expression is a real vector :math:`x` and :math:`p \neq q`, then it is treated as a matrix with a single row or a single column, depending on the shape associated with :math:`x`. See case (5). 5. If the normed expression is a complex vector :math:`z` and :math:`p = q`, then the definition is the same as in case (3) but with :math:`x = z`, :math:`|x_i|` the modulus of the :math:`i`-th entry of :math:`x`, and :math:`p \geq 1`. 6. If the normed expression is a real :math:`m \times n` matrix :math:`X`, then this is the :math:`L_{p,q}`-norm .. math:: \lVert X \rVert_{p,q} = \left(\sum_{j = 1}^n \left(\sum_{i = 1}^m |X_{ij}|^p \right)^{\frac{q}{p}} \right)^{\frac{1}{q}} for :math:`p, q \in \mathbb{Q} \cup \{\infty\}` with :math:`p,q \geq 1`. If :math:`p = q`, then this is equal to the (generalized) vector :math:`p`-norm of the the vectorized matrix, that is :math:`\lVert X \rVert_{p,p} = \lVert \operatorname{vec}(X) \rVert_p`. In this case, the requirement :math:`p \geq 1` is relaxed to :math:`p > 0` and :math:`X` may be a complex matrix. See case (3). Special cases: - For :math:`p = q = 2`, this is the Frobenius norm :math:`\lVert X \rVert = \lVert X \rVert_F`. - For :math:`p = 1`, :math:`q = \infty`, this is the maximum absolute column sum .. math:: \lVert X \rVert_{1,\infty} = \max_{j=1 \ldots n} \sum_{i=1}^m |X_{ij}|. This equals the operator norm induced by the vector :math:`1`-norm. You can obtain the maximum absolute row sum (the operator norm induced by the vector :math:`\infty`-norm) by first transposing :math:`X`. 7. Complex matrix norms are not supported. .. note:: You can write :math:`\infty` in Python as ``float("inf")``. """ # -------------------------------------------------------------------------- # Initialization and factory methods. # --------------------------------------------------------------------------
[docs] @convert_and_refine_arguments("x") def __init__(self, x, p=2, q=None, denominator_limit=1000): """Construct a :class:`Norm`. :param x: The affine expression to take the norm of. :type x: ~picos.expressions.ComplexAffineExpression :param float p: The value for :math:`p`, which is cast to a limited precision fraction. :param float q: The value for :math:`q`, which is cast to a limited precision fraction. The default of :obj:`None` means *equal to* :math:`p`. :param int denominator_limit: The largest allowed denominator when casting :math:`p` and :math:`q` to a fraction. Higher values can yield a greater precision at reduced performance. """ # Validate x. if not isinstance(x, ComplexAffineExpression): raise TypeError("Can only form the norm of an affine expression, " "not of {}.".format(type(x).__name__)) complex = not isinstance(x, AffineExpression) if isinstance(p, tuple) and len(p) == 2: throw_deprecation_warning("Arguments 'p' and 'q' to Norm must be " "given separately in the future.", decoratorLevel=1) p, q = p if q is None: q = p # Load p as a limtied precision fraction. if p == float("inf"): pNum = p pDen = 1 pStr = glyphs.infty() else: pNum, pDen, p, pStr = make_fraction(p, denominator_limit) # Load q as a limtied precision fraction. if q == float("inf"): qNum = q qDen = 1 qStr = glyphs.infty() else: qNum, qDen, q, qStr = make_fraction(q, denominator_limit) # Validate that p and q are in the allowed range. if p == q: if complex and p < 1: raise NotImplementedError( "Complex p-norm requires {}.".format(glyphs.ge("p", "1"))) elif p <= 0: raise ValueError( "p-norm requires {}.".format(glyphs.gt("p", "0"))) elif p != q: if complex: raise NotImplementedError( "(p,q)-norm is not supported for complex expressions.") elif p < 1 or q < 1: raise ValueError("(p,q)-norm requires {} and {}." .format(glyphs.ge("p", "1"), glyphs.ge("q", "1"))) # Build the string representations. if len(x) == 1: typeStr = "Modulus" if complex else "Absolute Value" symbStr = glyphs.abs(x.string) elif p == q: vec = 1 in x.shape if p == 1: typeStr = "Manhattan Norm" symbStr = glyphs.pnorm(x.string, "sum") elif p == 2: typeStr = "Euclidean Norm" if vec else "Frobenius Norm" symbStr = glyphs.norm(x.string) elif p == float("inf"): typeStr = "Maximum Norm" symbStr = glyphs.pnorm(x.string, "max") else: if p < 1: typeStr = "Generalized p-Norm" if vec else \ "Entrywise Generalized p-Norm" else: typeStr = "Vector p-Norm" if vec else "Entrywise p-Norm" symbStr = glyphs.pnorm(x.string, pStr) else: if p == 1 and q == float("inf"): typeStr = "Maximum Absolute Column Sum" else: typeStr = "Matrix (p,q)-Norm" symbStr = glyphs.pqnorm(x.string, pStr, qStr) if complex: typeStr = "Complex " + typeStr # Reduce the complex to the real case. if complex: if len(x) == 1: # From modulus to real Euclidean norm. x = x.real // x.imag p, pNum, pDen = 2.0, 2, 1 q, qNum, qDen = 2.0, 2, 1 else: # From complex vector/entrywise p-norm to real matrix p,q-norm. if x.shape[0] == 1: x = x.real // x.imag elif x.shape[1] == 1: x = x.T.real // x.T.imag else: x = x.vec.T.real // x.vec.T.imag assert p == q p, pNum, pDen = 2.0, 2, 1 # Reduce an entrywise p-norm to a vector p-norm and make all vector # norms refer to column vectors. if p == q and x.shape[1] != 1: x = x.vec assert isinstance(x, AffineExpression) assert float(pNum) / float(pDen) == p and float(qNum) / float(qDen) == q self._x = x self._pNum = pNum self._pDen = pDen self._qNum = qNum self._qDen = qDen self._limit = denominator_limit Expression.__init__(self, typeStr, symbStr)
# -------------------------------------------------------------------------- # Abstract method implementations and method overridings, except _predict. # -------------------------------------------------------------------------- @cached_unary_operator def _get_refined(self): if self._x.constant: return AffineExpression.from_constant(self.value, 1, self.string) else: return self Subtype = namedtuple("Subtype", ("xShape", "pNum", "pDen", "qNum", "qDen")) def _get_subtype(self): # NOTE: The xShape field refers to the internal (real and potentially # vectorized) representation of the normed expression. return self.Subtype( self._x.shape, self._pNum, self._pDen, self._qNum, self._qDen) def _get_value(self): value = self._x._get_value() if value.size == (1, 1): return abs(value) value = cvx2np(value) p, q = self.p, self.q if p == q: value = numpy.linalg.norm(numpy.ravel(value), p) else: columns = value.shape[1] value = numpy.linalg.norm([ numpy.linalg.norm(numpy.ravel(value[:, j]), p) for j in range(columns)], q) return cvxopt.matrix(value) def _get_mutables(self): return self._x._get_mutables() def _is_convex(self): return self.p >= 1 or len(self._x) == 1 def _is_concave(self): return self.p <= 1 and len(self._x) != 1 def _replace_mutables(self, mapping): return self.__class__( self._x._replace_mutables(mapping), self.p, self.q, self._limit) def _freeze_mutables(self, freeze): return self.__class__( self._x._freeze_mutables(freeze), self.p, self.q, self._limit) # -------------------------------------------------------------------------- # Python special method implementations, except constraint-creating ones. # -------------------------------------------------------------------------- @classmethod def _mul(cls, self, other, forward): if isinstance(other, AffineExpression) and other.constant: factor = other.safe_value if not factor: return AffineExpression.zero() elif factor == 1: return self elif factor > 0: if forward: string = glyphs.clever_mul(self.string, other.string) else: string = glyphs.clever_mul(other.string, self.string) norm = cls(other*self._x, self.p, self.q, self._limit) norm._typeStr = "Scaled " + norm._typeStr norm._symbStr = string return norm if forward: return Expression.__mul__(self, other) else: return Expression.__rmul__(self, other)
[docs] @convert_operands(scalarRHS=True) @refine_operands() def __mul__(self, other): return Norm._mul(self, other, True)
[docs] @convert_operands(scalarRHS=True) @refine_operands() def __rmul__(self, other): return Norm._mul(self, other, False)
[docs] @convert_operands(scalarRHS=True) @refine_operands() def __pow__(self, other): if isinstance(other, AffineExpression): if not other.constant or other.value != 2: raise NotImplementedError( "You may only take a norm to the power of two.") p, q = self.p, self.q if p == q and p == 2: result = SquaredNorm(self._x) result._symbStr = glyphs.squared(self.string) return result else: raise NotImplementedError( "You may only square an Euclidean or Frobenius norm.") return Expression.__pow__(self, other)
# -------------------------------------------------------------------------- # Properties and functions that describe the expression. # -------------------------------------------------------------------------- @property def p(self): """The parameter :math:`p`. This is a limited precision version of the parameter used when the norm was constructed. """ return float(self._pNum) / float(self._pDen) @property def pnum(self): """The limited precision fraction numerator of :math:`p`.""" return self._pNum @property def pden(self): """The limited precision fraction denominator of :math:`p`.""" return self._pDen @property def q(self): """The parameter :math:`q`. This is a limited precision version of the parameter used when the norm was constructed. """ return float(self._qNum) / float(self._qDen) @property def qnum(self): """The limited precision fraction numerator of :math:`q`.""" return self._qNum @property def qden(self): """The limited precision fraction denominator of :math:`q`.""" return self._qDen # -------------------------------------------------------------------------- # Methods and properties that return modified copies. # -------------------------------------------------------------------------- @property def x(self): """Real expression whose norm equals that of the original expression.""" return self._x # -------------------------------------------------------------------------- # Constraint-creating operators, and _predict. # -------------------------------------------------------------------------- @classmethod def _predict(cls, subtype, relation, other): assert isinstance(subtype, cls.Subtype) xShape, pNum, pDen, qNum, qDen = subtype xLen = xShape[0] * xShape[1] p = float(pNum) / float(pDen) q = float(qNum) / float(qDen) if relation == operator.__le__: if not (xLen == 1 or p >= 1): return NotImplemented # Not convex. if issubclass(other.clstype, AffineExpression) \ and other.subtype.dim == 1: if xLen == 1: return AbsoluteValueConstraint.make_type() elif p == q == 2: return SOCConstraint.make_type(argdim=xLen) elif p == q: return VectorNormConstraint.make_type( xLen, pNum, pDen, Constraint.LE) else: return MatrixNormConstraint.make_type( xShape, pNum, pDen, qNum, qDen) elif relation == operator.__ge__: if not (xLen != 1 and p <= 1): return NotImplemented # Not concave. assert p == q if issubclass(other.clstype, AffineExpression) \ and other.subtype.dim == 1: return VectorNormConstraint.make_type( xLen, pNum, pDen, Constraint.GE) return NotImplemented
[docs] @convert_operands(scalarRHS=True) @validate_prediction @refine_operands() def __le__(self, other): if not self.convex: raise TypeError("Cannot upper-bound a nonconvex generalized norm.") if isinstance(other, AffineExpression): if len(self._x) == 1: return AbsoluteValueConstraint(self._x, other) elif self.p == self.q == 2: # NOTE: The custom string is necessary in case the original x # was complex; otherwise the string would be something # like "‖vec([Re(xᵀ); Im(xᵀ)])‖ ≤ b". return SOCConstraint(self._x, other, customString=glyphs.le(self.string, other.string)) elif self.p == self.q: return VectorNormConstraint(self, Constraint.LE, other) else: return MatrixNormConstraint(self, other) else: return NotImplemented
[docs] @convert_operands(scalarRHS=True) @validate_prediction @refine_operands() def __ge__(self, other): if not self.concave: raise TypeError("Cannot lower-bound a nonconcave norm.") assert self.p == self.q if isinstance(other, AffineExpression): return VectorNormConstraint(self, Constraint.GE, other) else: return NotImplemented
# -------------------------------------- __all__ = api_end(_API_START, globals())