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

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:

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()

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)

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()

233 if isinstance(other, (AffineExpression, Entropy)):

235 # NOTE: __add__ always creates a fresh expression.

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:

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()

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)

455 return negent

456 else:

457 return NotImplemented

459 @convert_operands(scalarRHS=True)

460 @refine_operands()

462 if isinstance(other, (AffineExpression, NegativeEntropy)):

464 # NOTE: __add__ always creates a fresh expression.

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())