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

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# ------------------------------------------------------------------------------ 

19 

20"""Implements :class:`Entropy` and :class:`NegativeEntropy`.""" 

21 

22# TODO: Common base class for Entropy and NegativeEntropy. 

23 

24import math 

25import operator 

26from collections import namedtuple 

27 

28import cvxopt 

29import numpy 

30 

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 

38 

39_API_START = api_start(globals()) 

40# ------------------------------- 

41 

42 

43class Entropy(Expression): 

44 r"""Entropy or negative relative entropy of an affine expression. 

45 

46 Negative relative entropy is also known as the perspective of the logarithm. 

47 

48 :Definition: 

49 

50 Let :math:`x` be an :math:`n`-dimensional real affine expression. 

51 

52 1. If no additional expression :math:`y` is given, this is the entropy 

53 

54 .. math:: 

55 

56 -\sum_{i = 1}^n \operatorname{vec}(x)_i 

57 \log(\operatorname{vec}(x)_i). 

58 

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) 

62 

63 .. math:: 

64 

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

81 

82 .. warning:: 

83 

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 """ 

89 

90 # -------------------------------------------------------------------------- 

91 # Initialization and factory methods. 

92 # -------------------------------------------------------------------------- 

93 

94 @convert_and_refine_arguments("x", "y", allowNone=True) 

95 def __init__(self, x, y=None): 

96 """Construct an :class:`Entropy`. 

97 

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

107 

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) 

114 

115 if y.is1: 

116 y = None 

117 

118 self._x = x 

119 self._y = y 

120 

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

136 

137 Expression.__init__(self, typeStr, symbStr) 

138 

139 # -------------------------------------------------------------------------- 

140 # Abstract method implementations and method overridings, except _predict. 

141 # -------------------------------------------------------------------------- 

142 

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 

148 

149 Subtype = namedtuple("Subtype", ("argdim", "y")) 

150 

151 def _get_subtype(self): 

152 return self.Subtype(len(self._x), self._y is not None) 

153 

154 def _get_value(self): 

155 x = numpy.ravel(cvx2np(self._x._get_value())) 

156 

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

162 

163 return cvxopt.matrix(s) 

164 

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) 

171 

172 def _is_convex(self): 

173 return False 

174 

175 def _is_concave(self): 

176 return True 

177 

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

181 

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

185 

186 # -------------------------------------------------------------------------- 

187 # Python special method implementations, except constraint-creating ones. 

188 # -------------------------------------------------------------------------- 

189 

190 @cached_selfinverse_unary_operator 

191 def __neg__(self): 

192 return NegativeEntropy(self._x, self._y) 

193 

194 @convert_operands(scalarRHS=True) 

195 @refine_operands() 

196 def __add__(self, other): 

197 other_str = other.string 

198 

199 if isinstance(other, AffineExpression): 

200 if other.is0: 

201 return self 

202 

203 other = Entropy(other, other * math.e) 

204 

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) 

219 

220 entropy._symbStr = glyphs.clever_add(self.string, other_str) 

221 

222 return entropy 

223 

224 return Expression.__add__(self, other) 

225 

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) 

231 

232 if entropy is not self: 

233 entropy._symbStr = glyphs.clever_add(other.string, self.string) 

234 

235 return entropy 

236 

237 return Expression.__radd__(self, other) 

238 

239 @convert_operands(scalarRHS=True) 

240 @refine_operands() 

241 def __sub__(self, other): 

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

243 return self + (-other) 

244 

245 return Expression.__sub__(self, other) 

246 

247 # -------------------------------------------------------------------------- 

248 # Methods and properties that return expressions. 

249 # -------------------------------------------------------------------------- 

250 

251 @property 

252 def x(self): 

253 """The expression :math:`x`.""" 

254 return self._x 

255 

256 @property 

257 def y(self): 

258 """The additional expression :math:`y`, or :obj:`None`.""" 

259 return self._y 

260 

261 # -------------------------------------------------------------------------- 

262 # Methods and properties that describe the expression. 

263 # -------------------------------------------------------------------------- 

264 

265 @property 

266 def n(self): 

267 """Length of :attr:`x`.""" 

268 return len(self._x) 

269 

270 # -------------------------------------------------------------------------- 

271 # Constraint-creating operators, and _predict. 

272 # -------------------------------------------------------------------------- 

273 

274 @classmethod 

275 def _predict(cls, subtype, relation, other): 

276 assert isinstance(subtype, cls.Subtype) 

277 

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) 

283 

284 return NotImplemented 

285 

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 

294 

295 

296class NegativeEntropy(Expression): 

297 r"""Negative or relative entropy of an affine expression. 

298 

299 Relative entropy is also known as the Kullback-Leibler divergence. 

300 

301 :Definition: 

302 

303 Let :math:`x` be an :math:`n`-dimensional real affine expression. 

304 

305 1. If no additional expression :math:`y` is given, this is the negative 

306 entropy 

307 

308 .. math:: 

309 

310 \sum_{i = 1}^n \operatorname{vec}(x)_i 

311 \log(\operatorname{vec}(x)_i). 

312 

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) 

315 

316 .. math:: 

317 

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]. 

326 

327 .. warning:: 

328 

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 """ 

334 

335 # -------------------------------------------------------------------------- 

336 # Initialization and factory methods. 

337 # -------------------------------------------------------------------------- 

338 

339 @convert_and_refine_arguments("x", "y", allowNone=True) 

340 def __init__(self, x, y=None): 

341 """Construct a :class:`NegativeEntropy`. 

342 

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

352 

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) 

359 

360 if y.is1: 

361 y = None 

362 

363 self._x = x 

364 self._y = y 

365 

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

380 

381 Expression.__init__(self, typeStr, symbStr) 

382 

383 # -------------------------------------------------------------------------- 

384 # Abstract method implementations and method overridings, except _predict. 

385 # -------------------------------------------------------------------------- 

386 

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 

392 

393 Subtype = namedtuple("Subtype", ("argdim", "y")) 

394 

395 def _get_subtype(self): 

396 return self.Subtype(len(self._x), self._y is not None) 

397 

398 def _get_value(self): 

399 x = numpy.ravel(cvx2np(self._x._get_value())) 

400 

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

406 

407 return cvxopt.matrix(s) 

408 

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) 

415 

416 def _is_convex(self): 

417 return True 

418 

419 def _is_concave(self): 

420 return False 

421 

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

425 

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

429 

430 # -------------------------------------------------------------------------- 

431 # Python special method implementations, except constraint-creating ones. 

432 # -------------------------------------------------------------------------- 

433 

434 @cached_selfinverse_unary_operator 

435 def __neg__(self): 

436 return Entropy(self._x, self._y) 

437 

438 @convert_operands(scalarRHS=True) 

439 @refine_operands() 

440 def __add__(self, other): 

441 other_str = other.string 

442 

443 if isinstance(other, AffineExpression): 

444 if other.is0: 

445 return self 

446 

447 other = NegativeEntropy(other, other / math.e) 

448 

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) 

463 

464 negent._symbStr = glyphs.clever_add(self.string, other_str) 

465 

466 return negent 

467 

468 return Expression.__add__(self, other) 

469 

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) 

475 

476 if negent is not self: 

477 negent._symbStr = glyphs.clever_add(other.string, self.string) 

478 

479 return negent 

480 

481 return Expression.__radd__(self, other) 

482 

483 @convert_operands(scalarRHS=True) 

484 @refine_operands() 

485 def __sub__(self, other): 

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

487 return self + (-other) 

488 

489 return Expression.__sub__(self, other) 

490 

491 # -------------------------------------------------------------------------- 

492 # Methods and properties that return expressions. 

493 # -------------------------------------------------------------------------- 

494 

495 @property 

496 def x(self): 

497 """The expression :math:`x`.""" 

498 return self._x 

499 

500 @property 

501 def y(self): 

502 """The additional expression :math:`y`, or :obj:`None`.""" 

503 return self._y 

504 

505 # -------------------------------------------------------------------------- 

506 # Methods and properties that describe the expression. 

507 # -------------------------------------------------------------------------- 

508 

509 @property 

510 def n(self): 

511 """Length of :attr:`x`.""" 

512 return len(self._x) 

513 

514 # -------------------------------------------------------------------------- 

515 # Constraint-creating operators, and _predict. 

516 # -------------------------------------------------------------------------- 

517 

518 @classmethod 

519 def _predict(cls, subtype, relation, other): 

520 assert isinstance(subtype, cls.Subtype) 

521 

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) 

527 

528 return NotImplemented 

529 

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 

538 

539 

540# -------------------------------------- 

541__all__ = api_end(_API_START, globals())