Coverage for picos/expressions/exp_sqnorm.py : 77.70%
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) 2020 Maximilian Stahlberg
3#
4# This file is part of PICOS.
5#
6# PICOS is free software: you can redistribute it and/or modify it under the
7# terms of the GNU General Public License as published by the Free Software
8# Foundation, either version 3 of the License, or (at your option) any later
9# version.
10#
11# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY
12# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
14#
15# You should have received a copy of the GNU General Public License along with
16# this program. If not, see <http://www.gnu.org/licenses/>.
17# ------------------------------------------------------------------------------
19"""Implements :class:`SquaredNorm`."""
21import operator
22from collections import namedtuple
24import cvxopt
26from .. import glyphs
27from ..apidoc import api_end, api_start
28from ..caching import cached_property, cached_unary_operator
29from ..constraints import SquaredNormConstraint
30from .data import convert_and_refine_arguments, convert_operands, cvxopt_hcat
31from .exp_affine import AffineExpression, ComplexAffineExpression
32from .exp_quadratic import QuadraticExpression
33from .expression import Expression, refine_operands, validate_prediction
35_API_START = api_start(globals())
36# -------------------------------
39class SquaredNorm(QuadraticExpression):
40 """A squared Euclidean or Frobenius norm.
42 This is a lightweight wrapper around
43 :class:`~picos.expressions.QuadraticExpression` that can handle common
44 constraint formulations more efficiently.
45 """
47 # --------------------------------------------------------------------------
48 # Initialization and factory methods.
49 # --------------------------------------------------------------------------
51 @convert_and_refine_arguments("x")
52 def __init__(self, x):
53 """Create a squared Euclidean or Frobenius norm.
55 :param x:
56 The (complex) affine expression under the squared norm.
57 :type affinePart:
58 ~picos.expressions.ComplexAffineExpression
59 """
60 # Validate x.
61 if not isinstance(x, ComplexAffineExpression):
62 raise TypeError("Can only form the squared norm of an affine "
63 "expression, not of {}.".format(type(x).__name__))
65 if len(x) == 1:
66 typeStr = "Squared Scalar"
67 symbStr = glyphs.squared(x.string)
68 else:
69 typeStr = "Squared Norm"
70 symbStr = glyphs.squared(glyphs.norm(x.string))
72 Expression.__init__(self, typeStr, symbStr)
74 # TODO: Add a nonzero-vectorization to BiaffineExpression that returns
75 # a scalar zero for an all-zero expression and the vectorization
76 # of the expression with zero rows removed otherwise.
77 if x.is0:
78 self._x = AffineExpression.zero()
79 else:
80 # Vectorize and stack real and imaginary parts.
81 vec = x.vec if x.isreal else x.vec.real // x.vec.imag
83 # Remove zero rows from the vectorization.
84 A = abs(cvxopt_hcat(vec._coefs.values()))
85 a = cvxopt.sparse(sum(A[:, j] for j in range(A.size[1])))
86 nonzero = a.I
87 nnz = len(nonzero)
88 B = cvxopt.spmatrix([1.0]*nnz, range(nnz), nonzero, (nnz, len(x)))
90 self._x = B*vec
92 # --------------------------------------------------------------------------
93 # Allow inheriting from QuadraticExpression.
94 # --------------------------------------------------------------------------
96 @cached_property
97 def _quadratic_form(self):
98 """The squared norm as a pure quadratic expression.
100 If the expression under the norm is constant, then this is :obj:`None`.
101 """
102 # HACK: Make a shallow copy of self._x so that the product does not
103 # recognize that the operation below represents a squared norm.
104 # This only works as long as the product checks for operand
105 # equality with the "is" keyword.
106 result = (self._x.renamed("HACK") | self._x)
108 if isinstance(result, AffineExpression):
109 return None
110 else:
111 assert isinstance(result, QuadraticExpression)
112 result._symbStr = self.string
113 return result
115 @cached_property
116 def _quads(self):
117 if self._quadratic_form:
118 return self._quadratic_form._quads
119 else:
120 return {}
122 @cached_property
123 def _aff(self):
124 if self._quadratic_form:
125 return self._quadratic_form._aff
126 else:
127 refined = self.refined
128 assert isinstance(refined, ComplexAffineExpression)
129 return refined
131 @cached_property
132 def _sf(self):
133 if len(self._x) == 1:
134 return (self._x, self._x)
135 else:
136 return None
138 # --------------------------------------------------------------------------
139 # Squared norm specific properties.
140 # --------------------------------------------------------------------------
142 @property
143 def argdim(self):
144 """Number of nonzero elements of the expression under the norm."""
145 return len(self._x)
147 # --------------------------------------------------------------------------
148 # Abstract method implementations and method overridings, except _predict.
149 # --------------------------------------------------------------------------
151 @cached_unary_operator
152 def _get_refined(self):
153 if self._x.constant:
154 value = self._x.value_as_matrix
155 return AffineExpression.from_constant(
156 value.T*value, (1, 1), self._symbStr)
157 else:
158 return self
160 Subtype = namedtuple("Subtype", ("argdim", "quadratic_subtype"))
162 @cached_unary_operator
163 def _get_subtype(self):
164 return self.Subtype(
165 len(self._x),
166 self._quadratic_form.subtype if self._quadratic_form else None)
168 def _get_value(self):
169 value = self._x._get_value()
170 return value.T*value
172 @cached_unary_operator
173 def _get_variables(self):
174 return self._x.variables
176 def _is_convex(self):
177 return True
179 def _is_concave(self):
180 return self._x.constant
182 def _replace_variables(self, var_map):
183 return self.__class__(self._x._replace_variables(var_map))
185 # --------------------------------------------------------------------------
186 # Python special method implementations, except constraint-creating ones.
187 # --------------------------------------------------------------------------
189 def __len__(self):
190 # Faster version that overrides Expression.__len__.
191 return 1
193 @convert_operands(sameShape=True)
194 @refine_operands()
195 def __add__(self, other):
196 """Denote addition from the right hand side."""
197 if isinstance(other, SquaredNorm):
198 result = SquaredNorm(self._x // other._x)
199 result._symbStr = glyphs.clever_add(self.string, other.string)
200 return result
201 else:
202 return QuadraticExpression.__add__(self, other)
204 @convert_operands(scalarRHS=True)
205 @refine_operands()
206 def __mul__(self, other):
207 """Denote scaling from the right hand side."""
208 if isinstance(other, AffineExpression):
209 if not other.constant:
210 raise TypeError("You may only multiply a squared norm with a "
211 "constant term.")
213 value = other.value
215 if value < 0:
216 return QuadraticExpression.__mul__(self, other)
218 result = SquaredNorm(self._x*value**0.5)
219 result._symbStr = glyphs.clever_mul(self.string, other.string)
220 return result
221 else:
222 return QuadraticExpression.__mul__(self, other)
224 @convert_operands(scalarRHS=True)
225 @refine_operands()
226 def __rmul__(self, other):
227 """Denote scaling from the left hand side."""
228 if isinstance(other, AffineExpression):
229 result = self.__mul__(other)
230 # NOTE: __mul__ always creates a fresh expression.
231 result._symbStr = glyphs.clever_mul(other.string, self.string)
232 return result
233 else:
234 return QuadraticExpression.__rmul__(self, other)
236 @convert_operands(scalarRHS=True)
237 @refine_operands()
238 def __truediv__(self, other):
239 """Denote division by a constant scalar."""
240 if isinstance(other, AffineExpression):
241 if not other.constant:
242 raise TypeError("You may only divide a squared norm by a "
243 "constant term.")
245 result = self.__mul__(1.0 / other.value)
246 # NOTE: __mul__ always creates a fresh expression.
247 result._symbStr = glyphs.div(self.string, other.string)
248 return result
249 else:
250 return QuadraticExpression.__truediv__(self, other)
252 # --------------------------------------------------------------------------
253 # Method overridings for QuadraticExpression.
254 # --------------------------------------------------------------------------
256 @cached_property
257 def fullroot(self):
258 """Affine expression whose squared norm equals the expression.
260 Overrides :meth:`~.exp_quadratic.QuadraticExpression.fullroot` of
261 :class:`~.exp_quadratic.QuadraticExpression`.
262 """
263 return self._x.renamed(glyphs.Fn("fullroot({})")(self.string))
265 @property
266 def is_squared_norm(self):
267 """Always :obj:`True` for squared norm instances.
269 Overrides :meth:`~.exp_quadratic.QuadraticExpression.is_squared_norm`
270 of :class:`~.exp_quadratic.QuadraticExpression`.
271 """
272 return True
274 @property
275 def is0(self):
276 """Whether the expression is zero.
278 Overrides :meth:`~.exp_quadratic.QuadraticExpression.is0` of
279 :class:`~.exp_quadratic.QuadraticExpression`.
280 """
281 return self._x.is0
283 # --------------------------------------------------------------------------
284 # Constraint-creating operators, and _predict.
285 # --------------------------------------------------------------------------
287 @classmethod
288 def _predict(cls, subtype, relation, other):
289 assert isinstance(subtype, cls.Subtype)
291 if relation == operator.__le__:
292 if issubclass(other.clstype, AffineExpression) \
293 and other.subtype.dim == 1:
294 return SquaredNormConstraint.make_type(
295 subtype.argdim, other.subtype.constant)
297 return QuadraticExpression._predict(
298 subtype.quadratic_subtype, relation, other)
300 @convert_operands(sameShape=True)
301 @validate_prediction
302 @refine_operands()
303 def __le__(self, other):
304 if isinstance(other, AffineExpression):
305 return SquaredNormConstraint(self, other)
307 # NOTE: The following should handle the case where the upper bound has a
308 # scalar factorization efficiently by virtue of
309 # SquaredNorm.fullroot. See ConicQuadraticConstraint.
310 return QuadraticExpression.__le__(self, other)
313# --------------------------------------
314__all__ = api_end(_API_START, globals())