Coverage for picos/expressions/exp_entropy.py: 61.00%
241 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
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`."""
22# TODO: Common base class for Entropy and NegativeEntropy.
24import math
25import operator
26from collections import namedtuple
28import cvxopt
29import numpy
31from .. import glyphs
32from ..apidoc import api_end, api_start
33from ..caching import cached_selfinverse_unary_operator, cached_unary_operator
34from ..constraints import KullbackLeiblerConstraint
35from .data import convert_and_refine_arguments, convert_operands, cvx2np
36from .exp_affine import AffineExpression
37from .expression import Expression, refine_operands, validate_prediction
39_API_START = api_start(globals())
40# -------------------------------
43class Entropy(Expression):
44 r"""Entropy or negative relative entropy of an affine expression.
46 Negative relative entropy is also known as the perspective of the logarithm.
48 :Definition:
50 Let :math:`x` be an :math:`n`-dimensional real affine expression.
52 1. If no additional expression :math:`y` is given, this is the entropy
54 .. math::
56 -\sum_{i = 1}^n \operatorname{vec}(x)_i
57 \log(\operatorname{vec}(x)_i).
59 2. If an additional affine expression :math:`y` of same shape as :math:`x`
60 is given, this is the negative relative entropy (or logarithmic
61 perspective)
63 .. math::
65 -\sum_{i = 1}^n \operatorname{vec}(x)_i
66 \log\left(
67 \frac{\operatorname{vec}(x)_i}{\operatorname{vec}(y)_i}
68 \right)
69 &= -\sum_{i = 1}^n \operatorname{vec}(x)_i
70 \left[
71 \log(\operatorname{vec}(x)_i) - \log(\operatorname{vec}(y)_i)
72 \right] \\
73 &= \sum_{i = 1}^n \operatorname{vec}(x)_i
74 \left[
75 \log(\operatorname{vec}(y)_i) - \log(\operatorname{vec}(x)_i
76 \right] \\
77 &= \sum_{i = 1}^n \operatorname{vec}(x)_i
78 \log\left(
79 \frac{\operatorname{vec}(y)_i}{\operatorname{vec}(x)_i}
80 \right).
82 .. warning::
84 When you pose a lower bound on this expression, then PICOS enforces
85 :math:`x \geq 0` through an auxiliary constraint during solution search.
86 When an additional expression :math:`y` is given, PICOS enforces
87 :math:`y \geq 0` as well.
88 """
90 # --------------------------------------------------------------------------
91 # Initialization and factory methods.
92 # --------------------------------------------------------------------------
94 @convert_and_refine_arguments("x", "y", allowNone=True)
95 def __init__(self, x, y=None):
96 """Construct an :class:`Entropy`.
98 :param x: The affine expression :math:`x`.
99 :type x: ~picos.expressions.AffineExpression
100 :param y: An additional affine expression :math:`y`. If necessary, PICOS
101 will attempt to reshape or broadcast it to the shape of :math:`x`.
102 :type y: ~picos.expressions.AffineExpression
103 """
104 if not isinstance(x, AffineExpression):
105 raise TypeError("Can only take the elementwise logarithm of a real "
106 "affine expression, not of {}.".format(type(x).__name__))
108 if y is not None:
109 if not isinstance(y, AffineExpression):
110 raise TypeError("The additional parameter y must be a real "
111 "affine expression, not {}.".format(type(x).__name__))
112 elif x.shape != y.shape:
113 y = y.reshaped_or_broadcasted(x.shape)
115 if y.is1:
116 y = None
118 self._x = x
119 self._y = y
121 if y is None:
122 typeStr = "Entropy"
123 if len(x) == 1:
124 symbStr = glyphs.neg(glyphs.mul(x.string, glyphs.log(x.string)))
125 else:
126 symbStr = glyphs.neg(
127 glyphs.dotp(x.string, glyphs.log(x.string)))
128 else:
129 typeStr = "Logarithmic Perspective"
130 if len(x) == 1:
131 symbStr = glyphs.mul(
132 x.string, glyphs.log(glyphs.div(y.string, x.string)))
133 else:
134 symbStr = glyphs.dotp(x.string,
135 glyphs.sub(glyphs.log(y.string), glyphs.log(x.string)))
137 Expression.__init__(self, typeStr, symbStr)
139 # --------------------------------------------------------------------------
140 # Abstract method implementations and method overridings, except _predict.
141 # --------------------------------------------------------------------------
143 def _get_refined(self):
144 if self._x.constant and (self._y is None or self._y.constant):
145 return AffineExpression.from_constant(self.value, 1, self._symbStr)
146 else:
147 return self
149 Subtype = namedtuple("Subtype", ("argdim", "y"))
151 def _get_subtype(self):
152 return self.Subtype(len(self._x), self._y is not None)
154 def _get_value(self):
155 x = numpy.ravel(cvx2np(self._x._get_value()))
157 if self._y is None:
158 s = -numpy.dot(x, numpy.log(x))
159 else:
160 y = numpy.ravel(cvx2np(self._y._get_value()))
161 s = numpy.dot(x, numpy.log(y / x))
163 return cvxopt.matrix(s)
165 @cached_unary_operator
166 def _get_mutables(self):
167 if self._y is None:
168 return self._x._get_mutables()
169 else:
170 return self._x._get_mutables().union(self._y.mutables)
172 def _is_convex(self):
173 return False
175 def _is_concave(self):
176 return True
178 def _replace_mutables(self, mapping):
179 return self.__class__(self._x._replace_mutables(mapping),
180 None if self._y is None else self._y._replace_mutables(mapping))
182 def _freeze_mutables(self, freeze):
183 return self.__class__(self._x._freeze_mutables(freeze),
184 None if self._y is None else self._y._freeze_mutables(freeze))
186 # --------------------------------------------------------------------------
187 # Python special method implementations, except constraint-creating ones.
188 # --------------------------------------------------------------------------
190 @cached_selfinverse_unary_operator
191 def __neg__(self):
192 return NegativeEntropy(self._x, self._y)
194 @convert_operands(scalarRHS=True)
195 @refine_operands()
196 def __add__(self, other):
197 other_str = other.string
199 if isinstance(other, AffineExpression):
200 if other.is0:
201 return self
203 other = Entropy(other, other * math.e)
205 if isinstance(other, Entropy):
206 if self._y is None and other._y is None:
207 entropy = Entropy(self._x.vec // other._x.vec)
208 elif self._y is not None and other._y is None:
209 one = AffineExpression.from_constant(1.0, (other.n, 1))
210 entropy = Entropy(
211 self._x.vec // other._x.vec, self._y.vec // one)
212 elif self._y is None and other._y is not None:
213 one = AffineExpression.from_constant(1.0, (self.n, 1))
214 entropy = Entropy(
215 self._x.vec // other._x.vec, one // other._y.vec)
216 else:
217 entropy = Entropy(
218 self._x.vec // other._x.vec, self._y.vec // other._y.vec)
220 entropy._symbStr = glyphs.clever_add(self.string, other_str)
222 return entropy
224 return Expression.__add__(self, other)
226 @convert_operands(scalarRHS=True)
227 @refine_operands()
228 def __radd__(self, other):
229 if isinstance(other, (AffineExpression, Entropy)):
230 entropy = Entropy.__add__(self, other)
232 if entropy is not self:
233 entropy._symbStr = glyphs.clever_add(other.string, self.string)
235 return entropy
237 return Expression.__radd__(self, other)
239 @convert_operands(scalarRHS=True)
240 @refine_operands()
241 def __sub__(self, other):
242 if isinstance(other, (AffineExpression, NegativeEntropy)):
243 return self + (-other)
245 return Expression.__sub__(self, other)
247 # --------------------------------------------------------------------------
248 # Methods and properties that return expressions.
249 # --------------------------------------------------------------------------
251 @property
252 def x(self):
253 """The expression :math:`x`."""
254 return self._x
256 @property
257 def y(self):
258 """The additional expression :math:`y`, or :obj:`None`."""
259 return self._y
261 # --------------------------------------------------------------------------
262 # Methods and properties that describe the expression.
263 # --------------------------------------------------------------------------
265 @property
266 def n(self):
267 """Length of :attr:`x`."""
268 return len(self._x)
270 # --------------------------------------------------------------------------
271 # Constraint-creating operators, and _predict.
272 # --------------------------------------------------------------------------
274 @classmethod
275 def _predict(cls, subtype, relation, other):
276 assert isinstance(subtype, cls.Subtype)
278 if relation == operator.__ge__:
279 if issubclass(other.clstype, AffineExpression) \
280 and other.subtype.dim == 1:
281 return KullbackLeiblerConstraint.make_type(
282 argdim=subtype.argdim)
284 return NotImplemented
286 @convert_operands(scalarRHS=True)
287 @validate_prediction
288 @refine_operands()
289 def __ge__(self, other):
290 if isinstance(other, AffineExpression):
291 return KullbackLeiblerConstraint(-self, -other)
292 else:
293 return NotImplemented
296class NegativeEntropy(Expression):
297 r"""Negative or relative entropy of an affine expression.
299 Relative entropy is also known as the Kullback-Leibler divergence.
301 :Definition:
303 Let :math:`x` be an :math:`n`-dimensional real affine expression.
305 1. If no additional expression :math:`y` is given, this is the negative
306 entropy
308 .. math::
310 \sum_{i = 1}^n \operatorname{vec}(x)_i
311 \log(\operatorname{vec}(x)_i).
313 2. If an additional affine expression :math:`y` of same shape as :math:`x`
314 is given, this is the relative entropy (or Kullback-Leibler divergence)
316 .. math::
318 \sum_{i = 1}^n \operatorname{vec}(x)_i
319 \log\left(
320 \frac{\operatorname{vec}(x)_i}{\operatorname{vec}(y)_i}
321 \right)
322 = \sum_{i = 1}^n \operatorname{vec}(x)_i
323 \left[
324 \log(\operatorname{vec}(x)_i) - \log(\operatorname{vec}(y)_i)
325 \right].
327 .. warning::
329 When you pose an upper bound on this expression, then PICOS enforces
330 :math:`x \geq 0` through an auxiliary constraint during solution search.
331 When an additional expression :math:`y` is given, PICOS enforces
332 :math:`y \geq 0` as well.
333 """
335 # --------------------------------------------------------------------------
336 # Initialization and factory methods.
337 # --------------------------------------------------------------------------
339 @convert_and_refine_arguments("x", "y", allowNone=True)
340 def __init__(self, x, y=None):
341 """Construct a :class:`NegativeEntropy`.
343 :param x: The affine expression :math:`x`.
344 :type x: ~picos.expressions.AffineExpression
345 :param y: An additional affine expression :math:`y`. If necessary, PICOS
346 will attempt to reshape or broadcast it to the shape of :math:`x`.
347 :type y: ~picos.expressions.AffineExpression
348 """
349 if not isinstance(x, AffineExpression):
350 raise TypeError("Can only take the elementwise logarithm of a real "
351 "affine expression, not of {}.".format(type(x).__name__))
353 if y is not None:
354 if not isinstance(y, AffineExpression):
355 raise TypeError("The additional parameter y must be a real "
356 "affine expression, not {}.".format(type(x).__name__))
357 elif x.shape != y.shape:
358 y = y.reshaped_or_broadcasted(x.shape)
360 if y.is1:
361 y = None
363 self._x = x
364 self._y = y
366 if y is None:
367 typeStr = "Negative Entropy"
368 if len(x) == 1:
369 symbStr = glyphs.mul(x.string, glyphs.log(x.string))
370 else:
371 symbStr = glyphs.dotp(x.string, glyphs.log(x.string))
372 else:
373 typeStr = "Relative Entropy"
374 if len(x) == 1:
375 symbStr = glyphs.mul(
376 x.string, glyphs.log(glyphs.div(x.string, y.string)))
377 else:
378 symbStr = glyphs.dotp(x.string,
379 glyphs.sub(glyphs.log(x.string), glyphs.log(y.string)))
381 Expression.__init__(self, typeStr, symbStr)
383 # --------------------------------------------------------------------------
384 # Abstract method implementations and method overridings, except _predict.
385 # --------------------------------------------------------------------------
387 def _get_refined(self):
388 if self._x.constant and (self._y is None or self._y.constant):
389 return AffineExpression.from_constant(self.value, 1, self._symbStr)
390 else:
391 return self
393 Subtype = namedtuple("Subtype", ("argdim", "y"))
395 def _get_subtype(self):
396 return self.Subtype(len(self._x), self._y is not None)
398 def _get_value(self):
399 x = numpy.ravel(cvx2np(self._x._get_value()))
401 if self._y is None:
402 s = numpy.dot(x, numpy.log(x))
403 else:
404 y = numpy.ravel(cvx2np(self._y._get_value()))
405 s = numpy.dot(x, numpy.log(x / y))
407 return cvxopt.matrix(s)
409 @cached_unary_operator
410 def _get_mutables(self):
411 if self._y is None:
412 return self._x._get_mutables()
413 else:
414 return self._x._get_mutables().union(self._y.mutables)
416 def _is_convex(self):
417 return True
419 def _is_concave(self):
420 return False
422 def _replace_mutables(self, mapping):
423 return self.__class__(self._x._replace_mutables(mapping),
424 None if self._y is None else self._y._replace_mutables(mapping))
426 def _freeze_mutables(self, freeze):
427 return self.__class__(self._x._freeze_mutables(freeze),
428 None if self._y is None else self._y._freeze_mutables(freeze))
430 # --------------------------------------------------------------------------
431 # Python special method implementations, except constraint-creating ones.
432 # --------------------------------------------------------------------------
434 @cached_selfinverse_unary_operator
435 def __neg__(self):
436 return Entropy(self._x, self._y)
438 @convert_operands(scalarRHS=True)
439 @refine_operands()
440 def __add__(self, other):
441 other_str = other.string
443 if isinstance(other, AffineExpression):
444 if other.is0:
445 return self
447 other = NegativeEntropy(other, other / math.e)
449 if isinstance(other, NegativeEntropy):
450 if self._y is None and other._y is None:
451 negent = NegativeEntropy(self._x.vec // other._x.vec)
452 elif self._y is not None and other._y is None:
453 one = AffineExpression.from_constant(1.0, (other.n, 1))
454 negent = NegativeEntropy(
455 self._x.vec // other._x.vec, self._y.vec // one)
456 elif self._y is None and other._y is not None:
457 one = AffineExpression.from_constant(1.0, (self.n, 1))
458 negent = NegativeEntropy(
459 self._x.vec // other._x.vec, one // other._y.vec)
460 else:
461 negent = NegativeEntropy(
462 self._x.vec // other._x.vec, self._y.vec // other._y.vec)
464 negent._symbStr = glyphs.clever_add(self.string, other_str)
466 return negent
468 return Expression.__add__(self, other)
470 @convert_operands(scalarRHS=True)
471 @refine_operands()
472 def __radd__(self, other):
473 if isinstance(other, (AffineExpression, NegativeEntropy)):
474 negent = NegativeEntropy.__add__(self, other)
476 if negent is not self:
477 negent._symbStr = glyphs.clever_add(other.string, self.string)
479 return negent
481 return Expression.__radd__(self, other)
483 @convert_operands(scalarRHS=True)
484 @refine_operands()
485 def __sub__(self, other):
486 if isinstance(other, (AffineExpression, Entropy)):
487 return self + (-other)
489 return Expression.__sub__(self, other)
491 # --------------------------------------------------------------------------
492 # Methods and properties that return expressions.
493 # --------------------------------------------------------------------------
495 @property
496 def x(self):
497 """The expression :math:`x`."""
498 return self._x
500 @property
501 def y(self):
502 """The additional expression :math:`y`, or :obj:`None`."""
503 return self._y
505 # --------------------------------------------------------------------------
506 # Methods and properties that describe the expression.
507 # --------------------------------------------------------------------------
509 @property
510 def n(self):
511 """Length of :attr:`x`."""
512 return len(self._x)
514 # --------------------------------------------------------------------------
515 # Constraint-creating operators, and _predict.
516 # --------------------------------------------------------------------------
518 @classmethod
519 def _predict(cls, subtype, relation, other):
520 assert isinstance(subtype, cls.Subtype)
522 if relation == operator.__le__:
523 if issubclass(other.clstype, AffineExpression) \
524 and other.subtype.dim == 1:
525 return KullbackLeiblerConstraint.make_type(
526 argdim=subtype.argdim)
528 return NotImplemented
530 @convert_operands(scalarRHS=True)
531 @validate_prediction
532 @refine_operands()
533 def __le__(self, other):
534 if isinstance(other, AffineExpression):
535 return KullbackLeiblerConstraint(self, other)
536 else:
537 return NotImplemented
540# --------------------------------------
541__all__ = api_end(_API_START, globals())