Coverage for picos/expressions/exp_sqnorm.py: 89.61%
154 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
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 @convert_operands(sameShape=True)
190 @refine_operands()
191 def __add__(self, other):
192 """Denote addition from the right hand side."""
193 if isinstance(other, SquaredNorm): # No need to have __radd__ for this.
194 result = SquaredNorm(self._x // other._x)
195 result._symbStr = glyphs.clever_add(self.string, other.string)
197 return result
199 return QuadraticExpression.__add__(self, other)
201 @classmethod
202 def _mul_div(cls, self, other, div, forward):
203 assert not div or forward
205 if isinstance(other, AffineExpression) and other.constant:
206 factor = other.safe_value
208 if not factor:
209 if div:
210 raise ZeroDivisionError(
211 "Cannot divide {} by zero.".format(self.string))
212 else:
213 return AffineExpression.zero()
214 elif factor == 1:
215 return self
216 elif factor > 0:
217 if div:
218 factor = 1 / factor
219 string = glyphs.div(self.string, other.string)
220 elif forward:
221 string = glyphs.clever_mul(self.string, other.string)
222 else:
223 string = glyphs.clever_mul(other.string, self.string)
225 result = cls(self._x*factor**0.5)
226 result._typeStr = "Scaled " + result._typeStr
227 result._symbStr = string
229 return result
231 if div:
232 return QuadraticExpression.__truediv__(self, other)
233 elif forward:
234 return QuadraticExpression.__mul__(self, other)
235 else:
236 return QuadraticExpression.__rmul__(self, other)
238 @convert_operands(scalarRHS=True)
239 @refine_operands()
240 def __mul__(self, other):
241 """Denote scaling from the right hand side."""
242 return SquaredNorm._mul_div(self, other, div=False, forward=True)
244 @convert_operands(scalarRHS=True)
245 @refine_operands()
246 def __rmul__(self, other):
247 """Denote scaling from the left hand side."""
248 return SquaredNorm._mul_div(self, other, div=False, forward=False)
250 @convert_operands(scalarRHS=True)
251 @refine_operands()
252 def __truediv__(self, other):
253 """Denote division by a constant scalar."""
254 return SquaredNorm._mul_div(self, other, div=True, forward=True)
256 # --------------------------------------------------------------------------
257 # Method overridings for QuadraticExpression.
258 # --------------------------------------------------------------------------
260 @cached_property
261 def fullroot(self):
262 """Affine expression whose squared norm equals the expression.
264 Overrides :meth:`~.exp_quadratic.QuadraticExpression.fullroot` of
265 :class:`~.exp_quadratic.QuadraticExpression`.
266 """
267 return self._x.renamed(glyphs.Fn("fullroot({})")(self.string))
269 @property
270 def is_squared_norm(self):
271 """Always :obj:`True` for squared norm instances.
273 Overrides :meth:`~.exp_quadratic.QuadraticExpression.is_squared_norm`
274 of :class:`~.exp_quadratic.QuadraticExpression`.
275 """
276 return True
278 @property
279 def is0(self):
280 """Whether the expression is zero.
282 Overrides :meth:`~.exp_quadratic.QuadraticExpression.is0` of
283 :class:`~.exp_quadratic.QuadraticExpression`.
284 """
285 return self._x.is0
287 # --------------------------------------------------------------------------
288 # Constraint-creating operators, and _predict.
289 # --------------------------------------------------------------------------
291 @classmethod
292 def _predict(cls, subtype, relation, other):
293 assert isinstance(subtype, cls.Subtype)
295 if relation == operator.__le__:
296 if issubclass(other.clstype, AffineExpression) \
297 and other.subtype.dim == 1:
298 return SquaredNormConstraint.make_type(
299 subtype.argdim, other.subtype.constant)
301 return QuadraticExpression._predict(
302 subtype.quadratic_subtype, relation, other)
304 @convert_operands(sameShape=True)
305 @validate_prediction
306 @refine_operands()
307 def __le__(self, other):
308 if isinstance(other, AffineExpression):
309 return SquaredNormConstraint(self, other)
311 # NOTE: The following should handle the case where the upper bound has a
312 # scalar factorization efficiently by virtue of
313 # SquaredNorm.fullroot. See ConicQuadraticConstraint.
314 return QuadraticExpression.__le__(self, other)
317# --------------------------------------
318__all__ = api_end(_API_START, globals())