Coverage for picos/expressions/exp_entropy.py : 49.79%

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
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:`Entropy` and :class:`NegativeEntropy`."""
22import math
23import operator
24from collections import namedtuple
26import cvxopt
27import numpy
29from .. import glyphs
30from ..apidoc import api_end, api_start
31from ..caching import cached_selfinverse_unary_operator, cached_unary_operator
32from ..constraints import KullbackLeiblerConstraint
33from .data import convert_and_refine_arguments, convert_operands, cvx2np
34from .exp_affine import AffineExpression
35from .expression import Expression, refine_operands, validate_prediction
37_API_START = api_start(globals())
38# -------------------------------
41class Entropy(Expression):
42 r"""Entropy or negative relative entropy of an affine expression.
44 Negative relative entropy is also known as the perspective of the logarithm.
46 :Definition:
48 Let :math:`x` be an :math:`n`-dimensional real affine expression.
50 1. If no additional expression :math:`y` is given, this is the entropy
52 .. math::
54 -\sum_{i = 1}^n \operatorname{vec}(x)_i
55 \log(\operatorname{vec}(x)_i).
57 2. If an additional affine expression :math:`y` of same shape as :math:`x`
58 is given, this is the negative relative entropy (or logarithmic
59 perspective)
61 .. math::
63 -\sum_{i = 1}^n \operatorname{vec}(x)_i
64 \log\left(
65 \frac{\operatorname{vec}(x)_i}{\operatorname{vec}(y)_i}
66 \right)
67 &= -\sum_{i = 1}^n \operatorname{vec}(x)_i
68 \left[
69 \log(\operatorname{vec}(x)_i) - \log(\operatorname{vec}(y)_i)
70 \right] \\
71 &= \sum_{i = 1}^n \operatorname{vec}(x)_i
72 \left[
73 \log(\operatorname{vec}(y)_i) - \log(\operatorname{vec}(x)_i
74 \right] \\
75 &= \sum_{i = 1}^n \operatorname{vec}(x)_i
76 \log\left(
77 \frac{\operatorname{vec}(y)_i}{\operatorname{vec}(x)_i}
78 \right).
80 .. warning::
82 When you pose a lower bound on this expression, then PICOS enforces
83 :math:`x \geq 0` through an auxiliary constraint during solution search.
84 When an additional expression :math:`y` is given, PICOS enforces
85 :math:`y \geq 0` as well.
86 """
88 # --------------------------------------------------------------------------
89 # Initialization and factory methods.
90 # --------------------------------------------------------------------------
92 @convert_and_refine_arguments("x", "y", allowNone=True)
93 def __init__(self, x, y=None):
94 """Construct an :class:`Entropy`.
96 :param x: The affine expression :math:`x`.
97 :type x: ~picos.expressions.AffineExpression
98 :param y: An additional affine expression :math:`y`. If necessary, PICOS
99 will attempt to reshape or broadcast it to the shape of :math:`x`.
100 :type y: ~picos.expressions.AffineExpression
101 """
102 if not isinstance(x, AffineExpression):
103 raise TypeError("Can only take the elementwise logarithm of a real "
104 "affine expression, not of {}.".format(type(x).__name__))
106 if y is not None:
107 if not isinstance(y, AffineExpression):
108 raise TypeError("The additional parameter y must be a real "
109 "affine expression, not {}.".format(type(x).__name__))
110 elif x.shape != y.shape:
111 y = y.reshaped_or_broadcasted(x.shape)
113 if y.is1:
114 y = None
116 self._x = x
117 self._y = y
119 if y is None:
120 typeStr = "Entropy"
121 if len(x) == 1:
122 symbStr = glyphs.neg(glyphs.mul(x.string, glyphs.log(x.string)))
123 else:
124 symbStr = glyphs.neg(
125 glyphs.dotp(x.string, glyphs.log(x.string)))
126 else:
127 typeStr = "Logarithmic Perspective"
128 if len(x) == 1:
129 symbStr = glyphs.mul(
130 x.string, glyphs.log(glyphs.div(y.string, x.string)))
131 else:
132 symbStr = glyphs.dotp(x.string,
133 glyphs.sub(glyphs.log(y.string), glyphs.log(x.string)))
135 Expression.__init__(self, typeStr, symbStr)
137 # --------------------------------------------------------------------------
138 # Abstract method implementations and method overridings, except _predict.
139 # --------------------------------------------------------------------------
141 def _get_refined(self):
142 if self._x.constant and (self._y is None or self._y.constant):
143 return AffineExpression.from_constant(self.value, 1, self._symbStr)
144 else:
145 return self
147 Subtype = namedtuple("Subtype", ("argdim", "y"))
149 def _get_subtype(self):
150 return self.Subtype(len(self._x), self._y is not None)
152 def _get_value(self):
153 x = numpy.ravel(cvx2np(self._x._get_value()))
155 if self._y is None:
156 s = -numpy.dot(x, numpy.log(x))
157 else:
158 y = numpy.ravel(cvx2np(self._y._get_value()))
159 s = numpy.dot(x, numpy.log(y / x))
161 return cvxopt.matrix(s)
163 @cached_unary_operator
164 def _get_mutables(self):
165 if self._y is None:
166 return self._x._get_mutables()
167 else:
168 return self._x._get_mutables().union(self._y.mutables)
170 def _is_convex(self):
171 return False
173 def _is_concave(self):
174 return True
176 def _replace_mutables(self, mapping):
177 return self.__class__(self._x._replace_mutables(mapping),
178 None if self._y is None else self._y._replace_mutables(mapping))
180 def _freeze_mutables(self, freeze):
181 return self.__class__(self._x._freeze_mutables(freeze),
182 None if self._y is None else self._y._freeze_mutables(freeze))
184 # --------------------------------------------------------------------------
185 # Python special method implementations, except constraint-creating ones.
186 # --------------------------------------------------------------------------
188 @cached_selfinverse_unary_operator
189 def __neg__(self):
190 return NegativeEntropy(self._x, self._y)
192 @convert_operands(scalarRHS=True)
193 @refine_operands()
194 def __add__(self, other):
195 if isinstance(other, AffineExpression):
196 other = Entropy(other, other * math.e)
198 if isinstance(other, Entropy):
199 if self._y is None and other._y is None:
200 entropy = Entropy(self._x.vec // other._x.vec)
201 elif self._y is not None and other._y is None:
202 one = AffineExpression.from_constant(1.0, (other.n, 1))
203 entropy = Entropy(
204 self._x.vec // other._x.vec, self._y.vec // one)
205 elif self._y is None and other._y is not None:
206 one = AffineExpression.from_constant(1.0, (self.n, 1))
207 entropy = Entropy(
208 self._x.vec // other._x.vec, one // other._y.vec)
209 else:
210 entropy = Entropy(
211 self._x.vec // other._x.vec, self._y.vec // other._y.vec)
213 entropy._symbStr = glyphs.clever_add(self.string, other.string)
215 return entropy
216 else:
217 return NotImplemented
219 @convert_operands(scalarRHS=True)
220 @refine_operands()
221 def __sub__(self, other):
222 if isinstance(other, AffineExpression):
223 other = NegativeEntropy(other, other / math.e)
225 if isinstance(other, NegativeEntropy):
226 return self + (-other)
227 else:
228 return NotImplemented
230 @convert_operands(scalarRHS=True)
231 @refine_operands()
232 def __radd__(self, other):
233 if isinstance(other, (AffineExpression, Entropy)):
234 entropy = self.__add__(other)
235 # NOTE: __add__ always creates a fresh expression.
236 entropy._symbStr = glyphs.clever_add(other.string, self.string)
237 return entropy
238 else:
239 return NotImplemented
241 # --------------------------------------------------------------------------
242 # Methods and properties that return expressions.
243 # --------------------------------------------------------------------------
245 @property
246 def x(self):
247 """The expression :math:`x`."""
248 return self._x
250 @property
251 def y(self):
252 """The additional expression :math:`y`, or :obj:`None`."""
253 return self._y
255 # --------------------------------------------------------------------------
256 # Methods and properties that describe the expression.
257 # --------------------------------------------------------------------------
259 @property
260 def n(self):
261 """Length of :attr:`x`."""
262 return len(self._x)
264 # --------------------------------------------------------------------------
265 # Constraint-creating operators, and _predict.
266 # --------------------------------------------------------------------------
268 @classmethod
269 def _predict(cls, subtype, relation, other):
270 assert isinstance(subtype, cls.Subtype)
272 if relation == operator.__ge__:
273 if issubclass(other.clstype, AffineExpression) \
274 and other.subtype.dim == 1:
275 return KullbackLeiblerConstraint.make_type(
276 argdim=subtype.argdim)
278 return NotImplemented
280 @convert_operands(scalarRHS=True)
281 @validate_prediction
282 @refine_operands()
283 def __ge__(self, other):
284 if isinstance(other, AffineExpression):
285 return KullbackLeiblerConstraint(-self, -other)
286 else:
287 return NotImplemented
290class NegativeEntropy(Expression):
291 r"""Negative or relative entropy of an affine expression.
293 Relative entropy is also known as the Kullback-Leibler divergence.
295 :Definition:
297 Let :math:`x` be an :math:`n`-dimensional real affine expression.
299 1. If no additional expression :math:`y` is given, this is the negative
300 entropy
302 .. math::
304 \sum_{i = 1}^n \operatorname{vec}(x)_i
305 \log(\operatorname{vec}(x)_i).
307 2. If an additional affine expression :math:`y` of same shape as :math:`x`
308 is given, this is the relative entropy (or Kullback-Leibler divergence)
310 .. math::
312 \sum_{i = 1}^n \operatorname{vec}(x)_i
313 \log\left(
314 \frac{\operatorname{vec}(x)_i}{\operatorname{vec}(y)_i}
315 \right)
316 = \sum_{i = 1}^n \operatorname{vec}(x)_i
317 \left[
318 \log(\operatorname{vec}(x)_i) - \log(\operatorname{vec}(y)_i)
319 \right].
321 .. warning::
323 When you pose an upper bound on this expression, then PICOS enforces
324 :math:`x \geq 0` through an auxiliary constraint during solution search.
325 When an additional expression :math:`y` is given, PICOS enforces
326 :math:`y \geq 0` as well.
327 """
329 # --------------------------------------------------------------------------
330 # Initialization and factory methods.
331 # --------------------------------------------------------------------------
333 @convert_and_refine_arguments("x", "y", allowNone=True)
334 def __init__(self, x, y=None):
335 """Construct a :class:`NegativeEntropy`.
337 :param x: The affine expression :math:`x`.
338 :type x: ~picos.expressions.AffineExpression
339 :param y: An additional affine expression :math:`y`. If necessary, PICOS
340 will attempt to reshape or broadcast it to the shape of :math:`x`.
341 :type y: ~picos.expressions.AffineExpression
342 """
343 if not isinstance(x, AffineExpression):
344 raise TypeError("Can only take the elementwise logarithm of a real "
345 "affine expression, not of {}.".format(type(x).__name__))
347 if y is not None:
348 if not isinstance(y, AffineExpression):
349 raise TypeError("The additional parameter y must be a real "
350 "affine expression, not {}.".format(type(x).__name__))
351 elif x.shape != y.shape:
352 y = y.reshaped_or_broadcasted(x.shape)
354 if y.is1:
355 y = None
357 self._x = x
358 self._y = y
360 if y is None:
361 typeStr = "Negative Entropy"
362 if len(x) == 1:
363 symbStr = glyphs.mul(x.string, glyphs.log(x.string))
364 else:
365 symbStr = glyphs.dotp(x.string, glyphs.log(x.string))
366 else:
367 typeStr = "Relative Entropy"
368 if len(x) == 1:
369 symbStr = glyphs.mul(
370 x.string, glyphs.log(glyphs.div(x.string, y.string)))
371 else:
372 symbStr = glyphs.dotp(x.string,
373 glyphs.sub(glyphs.log(x.string), glyphs.log(y.string)))
375 Expression.__init__(self, typeStr, symbStr)
377 # --------------------------------------------------------------------------
378 # Abstract method implementations and method overridings, except _predict.
379 # --------------------------------------------------------------------------
381 def _get_refined(self):
382 if self._x.constant and (self._y is None or self._y.constant):
383 return AffineExpression.from_constant(self.value, 1, self._symbStr)
384 else:
385 return self
387 Subtype = namedtuple("Subtype", ("argdim", "y"))
389 def _get_subtype(self):
390 return self.Subtype(len(self._x), self._y is not None)
392 def _get_value(self):
393 x = numpy.ravel(cvx2np(self._x._get_value()))
395 if self._y is None:
396 s = numpy.dot(x, numpy.log(x))
397 else:
398 y = numpy.ravel(cvx2np(self._y._get_value()))
399 s = numpy.dot(x, numpy.log(x / y))
401 return cvxopt.matrix(s)
403 @cached_unary_operator
404 def _get_mutables(self):
405 if self._y is None:
406 return self._x._get_mutables()
407 else:
408 return self._x._get_mutables().union(self._y.mutables)
410 def _is_convex(self):
411 return True
413 def _is_concave(self):
414 return False
416 def _replace_mutables(self, mapping):
417 return self.__class__(self._x._replace_mutables(mapping),
418 None if self._y is None else self._y._replace_mutables(mapping))
420 def _freeze_mutables(self, freeze):
421 return self.__class__(self._x._freeze_mutables(freeze),
422 None if self._y is None else self._y._freeze_mutables(freeze))
424 # --------------------------------------------------------------------------
425 # Python special method implementations, except constraint-creating ones.
426 # --------------------------------------------------------------------------
428 @cached_selfinverse_unary_operator
429 def __neg__(self):
430 return Entropy(self._x, self._y)
432 @convert_operands(scalarRHS=True)
433 @refine_operands()
434 def __add__(self, other):
435 if isinstance(other, AffineExpression):
436 other = NegativeEntropy(other, other / math.e)
438 if isinstance(other, NegativeEntropy):
439 if self._y is None and other._y is None:
440 negent = NegativeEntropy(self._x.vec // other._x.vec)
441 elif self._y is not None and other._y is None:
442 one = AffineExpression.from_constant(1.0, (other.n, 1))
443 negent = NegativeEntropy(
444 self._x.vec // other._x.vec, self._y.vec // one)
445 elif self._y is None and other._y is not None:
446 one = AffineExpression.from_constant(1.0, (self.n, 1))
447 negent = NegativeEntropy(
448 self._x.vec // other._x.vec, one // other._y.vec)
449 else:
450 negent = NegativeEntropy(
451 self._x.vec // other._x.vec, self._y.vec // other._y.vec)
453 negent._symbStr = glyphs.clever_add(self.string, other.string)
455 return negent
456 else:
457 return NotImplemented
459 @convert_operands(scalarRHS=True)
460 @refine_operands()
461 def __radd__(self, other):
462 if isinstance(other, (AffineExpression, NegativeEntropy)):
463 negent = self.__add__(other)
464 # NOTE: __add__ always creates a fresh expression.
465 negent._symbStr = glyphs.clever_add(other.string, self.string)
466 return negent
467 else:
468 return NotImplemented
470 @convert_operands(scalarRHS=True)
471 @refine_operands()
472 def __sub__(self, other):
473 if isinstance(other, AffineExpression):
474 other = Entropy(other, other * math.e)
476 if isinstance(other, Entropy):
477 return self + (-other)
478 else:
479 return NotImplemented
481 # --------------------------------------------------------------------------
482 # Methods and properties that return expressions.
483 # --------------------------------------------------------------------------
485 @property
486 def x(self):
487 """The expression :math:`x`."""
488 return self._x
490 @property
491 def y(self):
492 """The additional expression :math:`y`, or :obj:`None`."""
493 return self._y
495 # --------------------------------------------------------------------------
496 # Methods and properties that describe the expression.
497 # --------------------------------------------------------------------------
499 @property
500 def n(self):
501 """Length of :attr:`x`."""
502 return len(self._x)
504 # --------------------------------------------------------------------------
505 # Constraint-creating operators, and _predict.
506 # --------------------------------------------------------------------------
508 @classmethod
509 def _predict(cls, subtype, relation, other):
510 assert isinstance(subtype, cls.Subtype)
512 if relation == operator.__le__:
513 if issubclass(other.clstype, AffineExpression) \
514 and other.subtype.dim == 1:
515 return KullbackLeiblerConstraint.make_type(
516 argdim=subtype.argdim)
518 return NotImplemented
520 @convert_operands(scalarRHS=True)
521 @validate_prediction
522 @refine_operands()
523 def __le__(self, other):
524 if isinstance(other, AffineExpression):
525 return KullbackLeiblerConstraint(self, other)
526 else:
527 return NotImplemented
530# --------------------------------------
531__all__ = api_end(_API_START, globals())