Coverage for picos/expressions/exp_norm.py : 77.27%
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) 2019 Maximilian Stahlberg, 2020 Guillaume Sagnol
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, NuclearNormConstraint,
33 SOCConstraint, SpectralNormConstraint,
34 VectorNormConstraint)
35from ..legacy import throw_deprecation_warning
36from .data import (convert_and_refine_arguments, convert_operands, cvx2np,
37 make_fraction)
38from .exp_affine import AffineExpression, ComplexAffineExpression
39from .exp_sqnorm import SquaredNorm
40from .expression import Expression, refine_operands, validate_prediction
42_API_START = api_start(globals())
43# -------------------------------
46class Norm(Expression):
47 r"""Entrywise :math:`p`-norm or :math:`L_{p,q}`-norm of an expression.
49 This class can represent the absolute value, the modulus, a vector
50 :math:`p`-norm and an entrywise matrix :math:`L_{p,q}`-norm of a (complex)
51 affine expression. In addition to these convex norms, it can represent the
52 concave *generalized vector* :math:`p`-*norm* with :math:`0 < p \leq 1`.
54 Not all of these norms are available on the complex field; see the
55 definitions below to learn more.
57 :Definition:
59 If :math:`q` is not given (:obj:`None`), then it is set equal to :math:`p`.
61 1. If the normed expression is a real scalar :math:`x`, then this is the
62 absolute value
64 .. math::
66 |x| = \begin{cases}
67 x, &\text{if }x \geq 0,\\
68 -x, &\text{otherwise.}
69 \end{cases}
71 The parameters :math:`p` and :math:`q` are ignored in this case.
73 2. If the normed expression is a complex scalar :math:`z`, then this is the
74 modulus
76 .. math::
78 |z| = \sqrt{\operatorname{Re}(z)^2 + \operatorname{Im}(z)^2}.
80 The parameters :math:`p` and :math:`q` are ignored in this case.
82 3. If the normed expression is a real vector :math:`x` and :math:`p = q`,
83 then this is the (generalized) vector :math:`p`-norm
85 .. math::
87 \lVert x \rVert_p = \left(\sum_{i=1}^n |x_i|^p\right)^{\frac{1}{p}}
89 for :math:`p \in \mathbb{Q} \cup \{\infty\}` with :math:`p > 0` and
90 :math:`|x_i|` the absolute value of the :math:`i`-th entry of :math:`x`.
92 Note that for :math:`p < 1` the expression is not convex and thus not a
93 proper norm. However, it is concave over the nonnegative orthant and
94 posing a lower bound on such a generalized norm yields a convex
95 constraint for :math:`x \geq 0`.
97 .. warning::
99 When you pose a lower bound on a concave generalized norm
100 (:math:`p < 1`), then PICOS enforces :math:`x \geq 0` through an
101 auxiliary constraint during solution search.
103 Special cases:
105 - For :math:`p = 1`, this is the *Manhattan* or *Taxicab* norm
106 :math:`\lVert x \rVert_{\text{sum}}`.
107 - For :math:`p = 2`, this is the Euclidean norm
108 :math:`\lVert x \rVert = \lVert x \rVert_2`.
109 - For :math:`p = \infty`, this is the *Maximum*, *Chebyshev*, or
110 *Infinity* norm :math:`\lVert x \rVert_{\text{max}}`.
112 4. If the normed expression is a real vector :math:`x` and
113 :math:`p \neq q`, then it is treated as a matrix with a single row or a
114 single column, depending on the shape associated with :math:`x`. See
115 case (5).
117 5. If the normed expression is a complex vector :math:`z` and
118 :math:`p = q`, then the definition is the same as in case (3) but with
119 :math:`x = z`, :math:`|x_i|` the modulus of the :math:`i`-th entry of
120 :math:`x`, and :math:`p \geq 1`.
122 6. If the normed expression is a real :math:`m \times n` matrix :math:`X`,
123 then this is the :math:`L_{p,q}`-norm
125 .. math::
127 \lVert X \rVert_{p,q} =
128 \left(\sum_{j = 1}^n \left(\sum_{i = 1}^m
129 |X_{ij}|^p
130 \right)^{\frac{q}{p}} \right)^{\frac{1}{q}}
132 for :math:`p, q \in \mathbb{Q} \cup \{\infty\}` with :math:`p,q \geq 1`.
134 If :math:`p = q`, then this is equal to the (generalized) vector
135 :math:`p`-norm of the the vectorized matrix, that is
136 :math:`\lVert X \rVert_{p,p} = \lVert \operatorname{vec}(X) \rVert_p`.
137 In this case, the requirement :math:`p \geq 1` is relaxed to
138 :math:`p > 0` and :math:`X` may be a complex matrix. See case (3).
140 Special cases:
142 - For :math:`p = q = 2`, this is the Frobenius norm
143 :math:`\lVert X \rVert = \lVert X \rVert_F`.
144 - For :math:`p = 1`, :math:`q = \infty`, this is the maximum absolute
145 column sum
147 .. math::
149 \lVert X \rVert_{1,\infty} =
150 \max_{j=1 \ldots n} \sum_{i=1}^m |X_{ij}|.
152 This equals the operator norm induced by the vector :math:`1`-norm.
153 You can obtain the maximum absolute row sum (the operator norm
154 induced by the vector :math:`\infty`-norm) by first transposing
155 :math:`X`.
157 7. Complex matrix norms are not supported.
159 .. note::
161 You can write :math:`\infty` in Python as ``float("inf")``.
162 """
164 # --------------------------------------------------------------------------
165 # Initialization and factory methods.
166 # --------------------------------------------------------------------------
168 @convert_and_refine_arguments("x")
169 def __init__(self, x, p=2, q=None, denominator_limit=1000):
170 """Construct a :class:`Norm`.
172 :param x: The affine expression to take the norm of.
173 :type x: ~picos.expressions.ComplexAffineExpression
174 :param float p: The value for :math:`p`, which is cast to a limited
175 precision fraction.
176 :param float q: The value for :math:`q`, which is cast to a limited
177 precision fraction. The default of :obj:`None` means *equal to*
178 :math:`p`.
179 :param int denominator_limit: The largest allowed denominator when
180 casting :math:`p` and :math:`q` to a fraction. Higher values can
181 yield a greater precision at reduced performance.
182 """
183 # Validate x.
184 if not isinstance(x, ComplexAffineExpression):
185 raise TypeError("Can only form the norm of an affine expression, "
186 "not of {}.".format(type(x).__name__))
188 complex = not isinstance(x, AffineExpression)
190 if isinstance(p, tuple) and len(p) == 2:
191 throw_deprecation_warning("Arguments 'p' and 'q' to Norm must be "
192 "given separately in the future.", decoratorLevel=1)
193 p, q = p
195 if q is None:
196 q = p
198 # Load p as a limtied precision fraction.
199 if p == float("inf"):
200 pNum = p
201 pDen = 1
202 pStr = glyphs.infty()
203 else:
204 pNum, pDen, p, pStr = make_fraction(p, denominator_limit)
206 # Load q as a limtied precision fraction.
207 if q == float("inf"):
208 qNum = q
209 qDen = 1
210 qStr = glyphs.infty()
211 else:
212 qNum, qDen, q, qStr = make_fraction(q, denominator_limit)
214 # Validate that p and q are in the allowed range.
215 if p == q:
216 if complex and p < 1:
217 raise NotImplementedError(
218 "Complex p-norm requires {}.".format(glyphs.ge("p", "1")))
219 elif p <= 0:
220 raise ValueError(
221 "p-norm requires {}.".format(glyphs.gt("p", "0")))
222 elif p != q:
223 if complex:
224 raise NotImplementedError(
225 "(p,q)-norm is not supported for complex expressions.")
226 elif p < 1 or q < 1:
227 raise ValueError("(p,q)-norm requires {} and {}."
228 .format(glyphs.ge("p", "1"), glyphs.ge("q", "1")))
230 # Build the string representations.
231 if len(x) == 1:
232 typeStr = "Modulus" if complex else "Absolute Value"
233 symbStr = glyphs.abs(x.string)
234 elif p == q:
235 vec = 1 in x.shape
236 if p == 1:
237 typeStr = "Manhattan Norm"
238 symbStr = glyphs.pnorm(x.string, "sum")
239 elif p == 2:
240 typeStr = "Euclidean Norm" if vec else "Frobenius Norm"
241 symbStr = glyphs.norm(x.string)
242 elif p == float("inf"):
243 typeStr = "Maximum Norm"
244 symbStr = glyphs.pnorm(x.string, "max")
245 else:
246 if p < 1:
247 typeStr = "Generalized p-Norm" if vec else \
248 "Entrywise Generalized p-Norm"
249 else:
250 typeStr = "Vector p-Norm" if vec else "Entrywise p-Norm"
251 symbStr = glyphs.pnorm(x.string, pStr)
252 else:
253 if p == 1 and q == float("inf"):
254 typeStr = "Maximum Absolute Column Sum"
255 else:
256 typeStr = "Matrix (p,q)-Norm"
257 symbStr = glyphs.pqnorm(x.string, pStr, qStr)
259 if complex:
260 typeStr = "Complex " + typeStr
262 # Reduce the complex to the real case.
263 if complex:
264 if len(x) == 1:
265 # From modulus to real Euclidean norm.
266 x = x.real // x.imag
268 p, pNum, pDen = 2.0, 2, 1
269 q, qNum, qDen = 2.0, 2, 1
270 else:
271 # From complex vector/entrywise p-norm to real matrix p,q-norm.
272 if x.shape[0] == 1:
273 x = x.real // x.imag
274 elif x.shape[1] == 1:
275 x = x.T.real // x.T.imag
276 else:
277 x = x.vec.T.real // x.vec.T.imag
279 assert p == q
280 p, pNum, pDen = 2.0, 2, 1
282 # Reduce an entrywise p-norm to a vector p-norm and make all vector
283 # norms refer to column vectors.
284 if p == q and x.shape[1] != 1:
285 x = x.vec
287 assert isinstance(x, AffineExpression)
288 assert float(pNum) / float(pDen) == p and float(qNum) / float(qDen) == q
290 self._x = x
291 self._pNum = pNum
292 self._pDen = pDen
293 self._qNum = qNum
294 self._qDen = qDen
295 self._limit = denominator_limit
297 Expression.__init__(self, typeStr, symbStr)
299 # --------------------------------------------------------------------------
300 # Abstract method implementations and method overridings, except _predict.
301 # --------------------------------------------------------------------------
303 @cached_unary_operator
304 def _get_refined(self):
305 if self._x.constant:
306 return AffineExpression.from_constant(self.value, 1, self.string)
307 else:
308 return self
310 Subtype = namedtuple("Subtype", ("xShape", "pNum", "pDen", "qNum", "qDen"))
312 def _get_subtype(self):
313 # NOTE: The xShape field refers to the internal (real and potentially
314 # vectorized) representation of the normed expression.
315 return self.Subtype(
316 self._x.shape, self._pNum, self._pDen, self._qNum, self._qDen)
318 def _get_value(self):
319 value = self._x._get_value()
321 if value.size == (1, 1):
322 return abs(value)
324 value = cvx2np(value)
326 p, q = self.p, self.q
328 if p == q:
329 value = numpy.linalg.norm(numpy.ravel(value), p)
330 else:
331 columns = value.shape[1]
332 value = numpy.linalg.norm([
333 numpy.linalg.norm(numpy.ravel(value[:, j]), p)
334 for j in range(columns)], q)
336 return cvxopt.matrix(value)
338 def _get_mutables(self):
339 return self._x._get_mutables()
341 def _is_convex(self):
342 return self.p >= 1 or len(self._x) == 1
344 def _is_concave(self):
345 return self.p <= 1 and len(self._x) != 1
347 def _replace_mutables(self, mapping):
348 return self.__class__(
349 self._x._replace_mutables(mapping), self.p, self.q, self._limit)
351 def _freeze_mutables(self, freeze):
352 return self.__class__(
353 self._x._freeze_mutables(freeze), self.p, self.q, self._limit)
355 # --------------------------------------------------------------------------
356 # Python special method implementations, except constraint-creating ones.
357 # --------------------------------------------------------------------------
359 @convert_operands(scalarRHS=True)
360 @refine_operands()
361 def __mul__(self, other):
362 if isinstance(other, AffineExpression):
363 if not other.constant:
364 raise NotImplementedError("You may only multiply a nonconstant "
365 "PICOS norm with a constant term.")
367 if other.value < 0:
368 raise NotImplementedError("You may only multiply a nonconstant "
369 "PICOS norm with a nonnegative term.")
371 norm = Norm(other*self._x, self.p, self.q, self._limit)
372 norm._typeStr = "Scaled " + norm._typeStr
373 norm._symbStr = glyphs.clever_mul(self.string, other.string)
374 return norm
375 else:
376 return NotImplemented
378 @convert_operands(scalarRHS=True)
379 @refine_operands()
380 def __rmul__(self, other):
381 if isinstance(other, AffineExpression):
382 norm = self.__mul__(other)
383 # NOTE: __mul__ always creates a fresh expression.
384 norm._symbStr = glyphs.clever_mul(other.string, self.string)
385 return norm
386 else:
387 return NotImplemented
389 @convert_operands(scalarRHS=True)
390 @refine_operands()
391 def __pow__(self, other):
392 if isinstance(other, AffineExpression):
393 if not other.constant or other.value != 2:
394 raise NotImplementedError("You may only take a PICOS norm to "
395 "the power of two.")
397 p, q = self.p, self.q
399 if p == q and p == 2:
400 result = SquaredNorm(self._x)
401 result._symbStr = glyphs.squared(self.string)
402 return result
403 else:
404 raise NotImplementedError(
405 "You may only square an Euclidean or Frobenius norm.")
406 else:
407 return NotImplemented
409 # --------------------------------------------------------------------------
410 # Properties and functions that describe the expression.
411 # --------------------------------------------------------------------------
413 @property
414 def p(self):
415 """The parameter :math:`p`.
417 This is a limited precision version of the parameter used when the norm
418 was constructed.
419 """
420 return float(self._pNum) / float(self._pDen)
422 @property
423 def pnum(self):
424 """The limited precision fraction numerator of :math:`p`."""
425 return self._pNum
427 @property
428 def pden(self):
429 """The limited precision fraction denominator of :math:`p`."""
430 return self._pDen
432 @property
433 def q(self):
434 """The parameter :math:`q`.
436 This is a limited precision version of the parameter used when the norm
437 was constructed.
438 """
439 return float(self._qNum) / float(self._qDen)
441 @property
442 def qnum(self):
443 """The limited precision fraction numerator of :math:`q`."""
444 return self._qNum
446 @property
447 def qden(self):
448 """The limited precision fraction denominator of :math:`q`."""
449 return self._qDen
451 # --------------------------------------------------------------------------
452 # Methods and properties that return modified copies.
453 # --------------------------------------------------------------------------
455 @property
456 def x(self):
457 """Real expression whose norm equals that of the original expression."""
458 return self._x
460 # --------------------------------------------------------------------------
461 # Constraint-creating operators, and _predict.
462 # --------------------------------------------------------------------------
464 @classmethod
465 def _predict(cls, subtype, relation, other):
466 assert isinstance(subtype, cls.Subtype)
468 xShape, pNum, pDen, qNum, qDen = subtype
469 xLen = xShape[0] * xShape[1]
470 p = float(pNum) / float(pDen)
471 q = float(qNum) / float(qDen)
473 if relation == operator.__le__:
474 if not (xLen == 1 or p >= 1):
475 return NotImplemented # Not convex.
477 if issubclass(other.clstype, AffineExpression) \
478 and other.subtype.dim == 1:
479 if xLen == 1:
480 return AbsoluteValueConstraint.make_type()
481 elif p == q == 2:
482 return SOCConstraint.make_type(argdim=xLen)
483 elif p == q:
484 return VectorNormConstraint.make_type(
485 xLen, pNum, pDen, Constraint.LE)
486 else:
487 return MatrixNormConstraint.make_type(
488 xShape, pNum, pDen, qNum, qDen)
489 elif relation == operator.__ge__:
490 if not (xLen != 1 and p <= 1):
491 return NotImplemented # Not concave.
493 assert p == q
495 if issubclass(other.clstype, AffineExpression) \
496 and other.subtype.dim == 1:
497 return VectorNormConstraint.make_type(
498 xLen, pNum, pDen, Constraint.GE)
500 return NotImplemented
502 @convert_operands(scalarRHS=True)
503 @validate_prediction
504 @refine_operands()
505 def __le__(self, other):
506 if not self.convex:
507 raise TypeError("Cannot upper-bound a nonconvex generalized norm.")
509 if isinstance(other, AffineExpression):
510 if len(self._x) == 1:
511 return AbsoluteValueConstraint(self._x, other)
512 elif self.p == self.q == 2:
513 # NOTE: The custom string is necessary in case the original x
514 # was complex; otherwise the string would be something
515 # like "‖vec([Re(xᵀ); Im(xᵀ)])‖ ≤ b".
516 return SOCConstraint(self._x, other,
517 customString=glyphs.le(self.string, other.string))
518 elif self.p == self.q:
519 return VectorNormConstraint(self, Constraint.LE, other)
520 else:
521 return MatrixNormConstraint(self, other)
522 else:
523 return NotImplemented
525 @convert_operands(scalarRHS=True)
526 @validate_prediction
527 @refine_operands()
528 def __ge__(self, other):
529 if not self.concave:
530 raise TypeError("Cannot lower-bound a nonconcave norm.")
532 assert self.p == self.q
534 if isinstance(other, AffineExpression):
535 return VectorNormConstraint(self, Constraint.GE, other)
536 else:
537 return NotImplemented
540class SpectralNorm(Expression):
541 r"""The spectral norm of a matrix.
543 This class can represent the spectral norm of a matrix-affine expression
544 (real- or complex valued). The spectral norm is convex, so we can form
545 expressions of the form ``SpectralNorm(X) <= t`` which are typically
546 reformulated as LMIs that can be handled by SDP solvers.
548 :Definition:
550 If the normed expression is a matrix :math:`X`, then its spectral norm is
552 .. math::
554 \|X\|_2 = \max \{ \|Xu\|_2 : \|u\| \leq 1\}
555 = \sqrt{\lambda_{\max}(XX^*)},
557 where :math:`\lambda_{\max}(\cdot)` denotes the largest eigenvalue of
558 a matrix, and :math:`X^*` denotes the adjoint matrix of :math:`X`
559 (i.e., the transposed matrix :math:`X^T` if :math:`X` is real-valued).
561 Special cases:
563 - If :math:`X` is scalar, then :math:`\|X\|_2` reduces to the the absolute
564 value (or modulus) :math:`|X|`.
565 - If :math:`X` is scalar, then :math:`\|X\|_2` coincides with the
566 Euclidean norm of :math:`X`.
568 """
570 @convert_and_refine_arguments("x")
571 def __init__(self, x):
572 """Construct a :class:`SpectralNorm`.
574 :param x: The affine expression to take the norm of.
575 :type x: ~picos.expressions.ComplexAffineExpression
576 """
577 # Validate x.
578 if not isinstance(x, ComplexAffineExpression):
579 raise TypeError("Can only form the spectral norm of an affine "
580 "expression, not of {}.".format(type(x).__name__))
582 complex = not isinstance(x, AffineExpression)
584 # Build the string representations.
585 if len(x) == 1:
586 typeStr = "Modulus" if complex else "Absolute Value"
587 symbStr = glyphs.abs(x.string)
588 elif 1 in x.shape:
589 typeStr = "Euclidean Norm"
590 symbStr = glyphs.norm(x.string)
591 else:
592 typeStr = "Spectral Norm"
593 symbStr = glyphs.spnorm(x.string)
595 if complex:
596 typeStr = "Complex " + typeStr
598 self._x = x
599 self._complex = complex
600 Expression.__init__(self, typeStr, symbStr)
602 # --------------------------------------------------------------------------
603 # Abstract method implementations and method overridings, except _predict.
604 # --------------------------------------------------------------------------
606 @cached_unary_operator
607 def _get_refined(self):
608 if self._x.constant:
609 return AffineExpression.from_constant(self.value, 1, self.string)
610 elif len(self._x) == 1 or (1 in self._x.shape):
611 return Norm(self._x)
612 else:
613 return self
615 Subtype = namedtuple("Subtype", ("argshape", "complex", "hermitian"))
617 def _get_subtype(self):
618 return self.Subtype(self._x.shape, self._complex, self._x.hermitian)
620 def _get_value(self):
621 value = self._x._get_value()
622 value = cvx2np(value)
623 value = numpy.linalg.norm(value, 2)
624 return cvxopt.matrix(value)
626 def _get_mutables(self):
627 return self._x._get_mutables()
629 def _is_convex(self):
630 return True
632 def _is_concave(self):
633 return False
635 def _replace_mutables(self, mapping):
636 return self.__class__(self._x._replace_mutables(mapping))
638 def _freeze_mutables(self, freeze):
639 return self.__class__(self._x._freeze_mutables(freeze))
641 # --------------------------------------------------------------------------
642 # Python special method implementations, except constraint-creating ones.
643 # --------------------------------------------------------------------------
645 @convert_operands(scalarRHS=True)
646 @refine_operands()
647 def __mul__(self, other):
648 if isinstance(other, AffineExpression):
649 if not other.constant:
650 raise NotImplementedError("You may only multiply a nonconstant "
651 "PICOS norm with a constant term.")
653 if other.value < 0:
654 raise NotImplementedError("You may only multiply a nonconstant "
655 "PICOS norm with a nonnegative term.")
657 norm = SpectralNorm(other*self._x)
658 norm._typeStr = "Scaled " + norm._typeStr
659 norm._symbStr = glyphs.clever_mul(self.string, other.string)
660 return norm
661 else:
662 return NotImplemented
664 @convert_operands(scalarRHS=True)
665 @refine_operands()
666 def __rmul__(self, other):
667 if isinstance(other, AffineExpression):
668 norm = self.__mul__(other)
669 # NOTE: __mul__ always creates a fresh expression.
670 norm._symbStr = glyphs.clever_mul(other.string, self.string)
671 return norm
672 else:
673 return NotImplemented
675 # --------------------------------------------------------------------------
676 # Methods and properties that return modified copies.
677 # --------------------------------------------------------------------------
679 @property
680 def x(self):
681 """Real expression whose norm equals that of the original expression."""
682 return self._x
684 # --------------------------------------------------------------------------
685 # Constraint-creating operators, and _predict.
686 # --------------------------------------------------------------------------
688 @classmethod
689 def _predict(cls, subtype, relation, other):
690 assert isinstance(subtype, cls.Subtype)
692 arg_shape, arg_complex, arg_hermitian = subtype
693 xLen = arg_shape[0] * arg_shape[1]
695 if relation == operator.__le__:
696 if issubclass(other.clstype, AffineExpression) \
697 and other.subtype.dim == 1:
698 if xLen == 1:
699 return AbsoluteValueConstraint.make_type()
700 elif 1 in arg_shape:
701 assert False, "Unexpected case (should have been refined)"
702 else:
703 return SpectralNormConstraint.make_type(
704 arg_shape, arg_complex, arg_hermitian)
705 elif relation == operator.__ge__:
706 return NotImplemented # Not concave.
708 return NotImplemented
710 @convert_operands(scalarRHS=True)
711 @validate_prediction
712 @refine_operands()
713 def __le__(self, other):
714 assert self.convex
716 if isinstance(other, AffineExpression):
717 if len(self._x) == 1:
718 return AbsoluteValueConstraint(self._x, other)
719 elif 1 in self._x.shape:
720 assert False, "Unexpected case (should have been refined)"
721 else:
722 return SpectralNormConstraint(self, other)
723 else:
724 return NotImplemented
727class NuclearNorm(Expression):
728 r"""The nuclear norm of a matrix.
730 This class can represent the nuclear norm of a matrix-affine expression
731 (real- or complex valued). The nuclear norm is convex, so we can form
732 expressions of the form ``NuclearNorm(X) <= t`` which are typically
733 reformulated as LMIs that can be handled by SDP solvers.
735 :Definition:
737 If the normed expression is a matrix :math:`X`, then its nuclear norm is
739 .. math::
741 \|X\|_* = \operatorname{trace}\ (X^*X)^{1/2}
742 = \sum_{i=1}^{\min(n,m)} \sigma_i(X)
744 where the :math:`\sigma_i(X)` denote the singular values of
745 a :math:`X`, and :math:`X^*` denotes the adjoint matrix of :math:`X`
746 (i.e., the transposed matrix :math:`X^T` if :math:`X` is real-valued).
748 Special cases:
750 - If :math:`X` is scalar, then :math:`\|X\|_*` reduces to the the absolute
751 value (or modulus) :math:`|X|`.
752 - If :math:`X` is scalar, then :math:`\|X\|_*` coincides with the
753 Euclidean norm of :math:`X`.
755 """
757 @convert_and_refine_arguments("x")
758 def __init__(self, x):
759 """Construct a :class:`NuclearNorm`.
761 :param x: The affine expression to take the norm of.
762 :type x: ~picos.expressions.ComplexAffineExpression
763 """
764 # Validate x.
765 if not isinstance(x, ComplexAffineExpression):
766 raise TypeError("Can only form the nuclear norm of an affine "
767 "expression, not of {}.".format(type(x).__name__))
769 complex = not isinstance(x, AffineExpression)
771 # Build the string representations.
772 if len(x) == 1:
773 typeStr = "Modulus" if complex else "Absolute Value"
774 symbStr = glyphs.abs(x.string)
775 elif 1 in x.shape:
776 typeStr = "Euclidean Norm"
777 symbStr = glyphs.norm(x.string)
778 else:
779 typeStr = "Nuclear Norm"
780 symbStr = glyphs.ncnorm(x.string)
782 if complex:
783 typeStr = "Complex " + typeStr
785 self._x = x
786 self._complex = complex
787 Expression.__init__(self, typeStr, symbStr)
789 # --------------------------------------------------------------------------
790 # Abstract method implementations and method overridings, except _predict.
791 # --------------------------------------------------------------------------
793 @cached_unary_operator
794 def _get_refined(self):
795 if self._x.constant:
796 return AffineExpression.from_constant(self.value, 1, self.string)
797 elif len(self._x) == 1 or (1 in self._x.shape):
798 return Norm(self._x)
799 else:
800 return self
802 Subtype = namedtuple("Subtype", ("argshape", "complex", "hermitian"))
804 def _get_subtype(self):
805 return self.Subtype(self._x.shape, self._complex, self._x.hermitian)
807 def _get_value(self):
808 value = self._x._get_value()
809 value = cvx2np(value)
810 value = numpy.linalg.norm(value, 'nuc')
811 return cvxopt.matrix(value)
813 def _get_mutables(self):
814 return self._x._get_mutables()
816 def _is_convex(self):
817 return True
819 def _is_concave(self):
820 return False
822 def _replace_mutables(self, mapping):
823 return self.__class__(self._x._replace_mutables(mapping))
825 def _freeze_mutables(self, freeze):
826 return self.__class__(self._x._freeze_mutables(freeze))
828 # --------------------------------------------------------------------------
829 # Python special method implementations, except constraint-creating ones.
830 # --------------------------------------------------------------------------
832 @convert_operands(scalarRHS=True)
833 @refine_operands()
834 def __mul__(self, other):
835 if isinstance(other, AffineExpression):
836 if not other.constant:
837 raise NotImplementedError("You may only multiply a nonconstant "
838 "PICOS norm with a constant term.")
840 if other.value < 0:
841 raise NotImplementedError("You may only multiply a nonconstant "
842 "PICOS norm with a nonnegative term.")
844 norm = NuclearNorm(other*self._x)
845 norm._typeStr = "Scaled " + norm._typeStr
846 norm._symbStr = glyphs.clever_mul(self.string, other.string)
847 return norm
848 else:
849 return NotImplemented
851 @convert_operands(scalarRHS=True)
852 @refine_operands()
853 def __rmul__(self, other):
854 if isinstance(other, AffineExpression):
855 norm = self.__mul__(other)
856 # NOTE: __mul__ always creates a fresh expression.
857 norm._symbStr = glyphs.clever_mul(other.string, self.string)
858 return norm
859 else:
860 return NotImplemented
862 # --------------------------------------------------------------------------
863 # Methods and properties that return modified copies.
864 # --------------------------------------------------------------------------
866 @property
867 def x(self):
868 """Real expression whose norm equals that of the original expression."""
869 return self._x
871 # --------------------------------------------------------------------------
872 # Constraint-creating operators, and _predict.
873 # --------------------------------------------------------------------------
875 @classmethod
876 def _predict(cls, subtype, relation, other):
877 assert isinstance(subtype, cls.Subtype)
879 arg_shape, arg_complex, arg_hermitian = subtype
880 xLen = arg_shape[0] * arg_shape[1]
882 if relation == operator.__le__:
883 if issubclass(other.clstype, AffineExpression) \
884 and other.subtype.dim == 1:
885 if xLen == 1:
886 return AbsoluteValueConstraint.make_type()
887 elif 1 in arg_shape:
888 assert False, "Unexpected case (should have been refined)"
889 else:
890 return NuclearNormConstraint.make_type(
891 arg_shape, arg_complex, arg_hermitian)
892 elif relation == operator.__ge__:
893 return NotImplemented # Not concave.
895 return NotImplemented
897 @convert_operands(scalarRHS=True)
898 @validate_prediction
899 @refine_operands()
900 def __le__(self, other):
901 assert self.convex
903 if isinstance(other, AffineExpression):
904 if len(self._x) == 1:
905 return AbsoluteValueConstraint(self._x, other)
906 elif 1 in self._x.shape:
907 assert False, "Unexpected case (should have been refined)"
908 else:
909 return NuclearNormConstraint(self, other)
910 else:
911 return NotImplemented
914# --------------------------------------
915__all__ = api_end(_API_START, globals())