Coverage for picos/expressions/exp_quantentr.py: 76.62%

1# ------------------------------------------------------------------------------ 

2# Copyright (C) 2024 Kerry He 

3# 

4# This file is part of PICOS. 

5# 

6# PICOS is free software: you can redistribute it and/or modify it under the 

7# terms of the GNU General Public License as published by the Free Software 

8# Foundation, either version 3 of the License, or (at your option) any later 

9# version. 

10# 

11# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY 

12# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 

13# A PARTICULAR PURPOSE. See the GNU General Public License for more details. 

14# 

15# You should have received a copy of the GNU General Public License along with 

16# this program. If not, see <http://www.gnu.org/licenses/>. 

17# ------------------------------------------------------------------------------ 

18 

19"""Implements :class:`QuantumEntropy`, :class:`NegativeQuantumEntropy`.""" 

20 

21# TODO: Common base class for QuantumEntropy and NegativeQuantumEntropy. 

22 

23import operator 

24from collections import namedtuple 

25 

26import cvxopt 

27import numpy 

28 

29from .. import glyphs 

30from ..apidoc import api_end, api_start 

31from ..caching import cached_selfinverse_unary_operator, cached_unary_operator 

32from ..constraints import ( 

33 QuantRelEntropyConstraint, 

34 ComplexQuantRelEntropyConstraint, 

35) 

36from .data import convert_and_refine_arguments, convert_operands, cvx2np 

37from .exp_affine import AffineExpression, ComplexAffineExpression 

38from .expression import Expression, refine_operands, validate_prediction 

39 

40_API_START = api_start(globals()) 

41# ------------------------------- 

42 

43 

44class QuantumEntropy(Expression): 

45 r"""Quantum or negative quantum relative entropy of an affine expression. 

46 

47 :Definition: 

48 

49 Let :math:`X` be an :math:`n \times n`-dimensional symmetric or hermitian 

50 matrix. 

51 

52 1. If no additional expression :math:`Y` is given, this is the quantum 

53 entropy 

54 

55 .. math:: 

56 

57 -\operatorname{Tr}(X \log(X)). 

58 

59 2. If an additional affine expression :math:`Y` of same shape as :math:`X` 

60 is given, this is the negative quantum relative entropy 

61 

62 .. math:: 

63 

64 \operatorname{Tr}(X \log(Y) - X\log(X)) 

65 

66 3. If an additional scalar valued real affine expression :math:`Y` 

67 is given, this is the homogenized quantum entropy 

68 

69 .. math:: 

70 

71 -\operatorname{Tr}(X \log(X/y)) 

72 

73 .. warning:: 

74 

75 When you pose a lower bound on this expression, then PICOS enforces 

76 :math:`X \succeq 0` through an auxiliary constraint during solution 

77 search. When an additional expression :math:`Y` is given, PICOS enforces 

78 :math:`Y \succeq 0` as well. 

79 """ 

80 

81 # -------------------------------------------------------------------------- 

82 # Initialization and factory methods. 

83 # -------------------------------------------------------------------------- 

84 

85 @convert_and_refine_arguments("X", "Y", allowNone=True) 

86 def __init__(self, X, Y=None): 

87 """Construct an :class:`QuantumEntropy`. 

88 

89 :param X: The affine expression :math:`X`. 

90 :type X: ~picos.expressions.AffineExpression 

91 :param Y: An additional affine expression :math:`Y`. If necessary, PICOS 

92 will attempt to reshape or broadcast it to the shape of :math:`X`. 

93 :type Y: ~picos.expressions.AffineExpression 

94 """ 

95 if not isinstance(X, ComplexAffineExpression): 

96 raise TypeError( 

97 "Can only take the matrix logarithm of a real " 

98 "or complex affine expression, not of {}.".format( 

99 type(X).__name__ 

100 ) 

101 ) 

102 if not X.hermitian: 

103 raise TypeError( 

104 "Can only take the matrix logarithm of a symmetric " 

105 "or Hermitian expression, not of {}.".format(type(X).__name__) 

106 ) 

107 

108 if Y is not None: 

109 if not isinstance(Y, ComplexAffineExpression): 

110 raise TypeError( 

111 "The additional parameter Y must be a real " 

112 "or complex affine expression, not {}.".format( 

113 type(Y).__name__ 

114 ) 

115 ) 

116 if not Y.hermitian: 

117 raise TypeError( 

118 "Can only take the matrix logarithm of a symmetric " 

119 "or Hermitian expression, not of {}.".format( 

120 type(Y).__name__ 

121 ) 

122 ) 

123 if X.shape != Y.shape and Y.shape != (1, 1): 

124 raise TypeError( 

125 "The additional parameter Y must either be the " 

126 "same shape as X, or be a real scalar expression, " 

127 "not {}.".format(type(Y).__name__) 

128 ) 

129 if Y.is1: 

130 Y = None 

131 

132 self._X = X 

133 self._Y = Y 

134 

135 self._iscomplex = not isinstance(X, AffineExpression) or \ 

136 not isinstance(Y, AffineExpression) 

137 

138 if Y is None: 

139 typeStr = "Quantum Entropy" 

140 symbStr = glyphs.qe(X.string) 

141 elif Y.shape == (1, 1): 

142 typeStr = "Homogenized Quantum Entropy" 

143 symbStr = glyphs.neg( 

144 glyphs.qre(X.string, glyphs.mul(Y.string, glyphs.idmatrix))) 

145 else: 

146 typeStr = "Negative Quantum Relative Entropy" 

147 symbStr = glyphs.neg(glyphs.qre(X.string, Y.string)) 

148 

149 Expression.__init__(self, typeStr, symbStr) 

150 

151 # -------------------------------------------------------------------------- 

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

153 # -------------------------------------------------------------------------- 

154 

155 def _get_refined(self): 

156 if self._X.constant and (self._Y is None or self._Y.constant): 

157 return AffineExpression.from_constant(self.value, 1, self._symbStr) 

158 else: 

159 return self 

160 

161 Subtype = namedtuple("Subtype", ("argdim", "Y", "iscomplex")) 

162 

163 def _get_subtype(self): 

164 return self.Subtype(len(self._X), self._Y is not None, self._iscomplex) 

165 

166 def _get_value(self): 

167 X = cvx2np(self._X._get_value()) 

168 eigvalsX, eigvecsX = numpy.linalg.eigh(X) 

169 

170 if self._Y is None: 

171 eigvalsX = eigvalsX[eigvalsX > 1e-12] 

172 s = -numpy.sum(eigvalsX * numpy.log(eigvalsX)) 

173 else: 

174 Y = eigvecsX.conj().T @ cvx2np(self._Y._get_value()) @ eigvecsX 

175 

176 Dy, Uy = numpy.linalg.eigh(Y) 

177 logY = Uy @ numpy.diag(numpy.log(Dy)) @ Uy.conj().T 

178 

179 s = numpy.sum(numpy.diag(eigvalsX) * logY.conj()).real 

180 eigvalsX = eigvalsX[eigvalsX > 1e-12] 

181 s -= numpy.sum(eigvalsX * numpy.log(eigvalsX)) 

182 

183 return cvxopt.matrix(s) 

184 

185 @cached_unary_operator 

186 def _get_mutables(self): 

187 if self._Y is None: 

188 return self._X._get_mutables() 

189 else: 

190 return self._X._get_mutables().union(self._Y.mutables) 

191 

192 def _is_convex(self): 

193 return False 

194 

195 def _is_concave(self): 

196 return True 

197 

198 def _replace_mutables(self, mapping): 

199 return self.__class__( 

200 self._X._replace_mutables(mapping), 

201 None if self._Y is None else self._Y._replace_mutables(mapping), 

202 ) 

203 

204 def _freeze_mutables(self, freeze): 

205 return self.__class__( 

206 self._X._freeze_mutables(freeze), 

207 None if self._Y is None else self._Y._freeze_mutables(freeze), 

208 ) 

209 

210 # -------------------------------------------------------------------------- 

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

212 # -------------------------------------------------------------------------- 

213 

214 @cached_selfinverse_unary_operator 

215 def __neg__(self): 

216 return NegativeQuantumEntropy(self._X, self._Y) 

217 

218 # -------------------------------------------------------------------------- 

219 # Methods and properties that return expressions. 

220 # -------------------------------------------------------------------------- 

221 

222 @property 

223 def X(self): 

224 """The expression :math:`X`.""" 

225 return self._X 

226 

227 @property 

228 def Y(self): 

229 """The additional expression :math:`Y`, or :obj:`None`.""" 

230 return self._Y 

231 

232 # -------------------------------------------------------------------------- 

233 # Methods and properties that describe the expression. 

234 # -------------------------------------------------------------------------- 

235 

236 @property 

237 def n(self): 

238 """Length of :attr:`X`.""" 

239 return len(self._X) 

240 

241 @property 

242 def iscomplex(self): 

243 """Whether :attr:`X` and :attr:`Y` are complex expressions or not.""" 

244 return self._iscomplex 

245 

246 # -------------------------------------------------------------------------- 

247 # Constraint-creating operators, and _predict. 

248 # -------------------------------------------------------------------------- 

249 

250 @classmethod 

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

252 assert isinstance(subtype, cls.Subtype) 

253 

254 if relation == operator.__ge__: 

255 if ( 

256 issubclass(other.clstype, AffineExpression) 

257 and other.subtype.dim == 1 

258 ): 

259 if subtype.iscomplex: 

260 return ComplexQuantRelEntropyConstraint.make_type( 

261 argdim=subtype.argdim 

262 ) 

263 else: 

264 return QuantRelEntropyConstraint.make_type( 

265 argdim=subtype.argdim 

266 ) 

267 return NotImplemented 

268 

269 @convert_operands(scalarRHS=True) 

270 @validate_prediction 

271 @refine_operands() 

272 def __ge__(self, other): 

273 if isinstance(other, AffineExpression): 

274 if self.iscomplex: 

275 return ComplexQuantRelEntropyConstraint(-self, -other) 

276 else: 

277 return QuantRelEntropyConstraint(-self, -other) 

278 else: 

279 return NotImplemented 

280 

281 

282class NegativeQuantumEntropy(Expression): 

283 r"""Negative or quantum relative entropy of an affine expression. 

284 

285 :Definition: 

286 

287 Let :math:`X` be an :math:`n \times n`-dimensional symmetric or hermitian 

288 matrix. 

289 

290 1. If no additional expression :math:`Y` is given, this is the negative 

291 quantum entropy 

292 

293 .. math:: 

294 

295 \operatorname{Tr}(X \log(X)). 

296 

297 2. If an additional affine expression :math:`Y` of same shape as :math:`X` 

298 is given, this is the quantum relative entropy 

299 

300 .. math:: 

301 

302 \operatorname{Tr}(X\log(X) - X\log(Y)). 

303 

304 3. If an additional scalar valued real affine expression :math:`Y` 

305 is given, this is the homogenized negative quantum entropy 

306 

307 .. math:: 

308 

309 \operatorname{Tr}(X \log(X/y)). 

310 

311 .. warning:: 

312 

313 When you pose an upper bound on this expression, then PICOS enforces 

314 :math:`X \succeq 0` through an auxiliary constraint during solution 

315 search. When an additional expression :math:`Y` is given, PICOS enforces 

316 :math:`Y \succeq 0` as well. 

317 """ 

318 

319 # -------------------------------------------------------------------------- 

320 # Initialization and factory methods. 

321 # -------------------------------------------------------------------------- 

322 

323 @convert_and_refine_arguments("X", "Y", allowNone=True) 

324 def __init__(self, X, Y=None): 

325 """Construct a :class:`NegativeQuantumEntropy`. 

326 

327 :param X: The affine expression :math:`X`. 

328 :type X: ~picos.expressions.AffineExpression 

329 :param Y: An additional affine expression :math:`Y`. If necessary, PICOS 

330 will attempt to reshape or broadcast it to the shape of :math:`X`. 

331 :type Y: ~picos.expressions.AffineExpression 

332 """ 

333 if not isinstance(X, ComplexAffineExpression): 

334 raise TypeError( 

335 "Can only take the matrix logarithm of a real " 

336 "or complex affine expression, not of {}.".format( 

337 type(X).__name__ 

338 ) 

339 ) 

340 if not X.hermitian: 

341 raise TypeError( 

342 "Can only take the matrix logarithm of a symmetric " 

343 "or Hermitian expression, not of {}.".format(type(X).__name__) 

344 ) 

345 

346 if Y is not None: 

347 if not isinstance(Y, ComplexAffineExpression): 

348 raise TypeError( 

349 "The additional parameter Y must be a real " 

350 "or complex affine expression, not {}.".format( 

351 type(Y).__name__ 

352 ) 

353 ) 

354 if not Y.hermitian: 

355 raise TypeError( 

356 "Can only take the matrix logarithm of a symmetric " 

357 "or Hermitian expression, not of {}.".format( 

358 type(Y).__name__ 

359 ) 

360 ) 

361 if X.shape != Y.shape and Y.shape != (1, 1): 

362 raise TypeError( 

363 "The additional parameter Y must either be the " 

364 "same shape as X, or be a real scalar expression, " 

365 "not {}.".format(type(Y).__name__) 

366 ) 

367 if Y.is1: 

368 Y = None 

369 

370 self._X = X 

371 self._Y = Y 

372 

373 self._iscomplex = not isinstance(X, AffineExpression) or ( 

374 not isinstance(Y, AffineExpression) 

375 ) 

376 

377 if Y is None: 

378 typeStr = "Negative Quantum Entropy" 

379 symbStr = glyphs.neg(glyphs.qe(X.string)) 

380 elif Y.shape == (1, 1): 

381 typeStr = "Negative Homogenized Quantum Entropy" 

382 symbStr = glyphs.qre( 

383 X.string, glyphs.mul(Y.string, glyphs.idmatrix)) 

384 else: 

385 typeStr = "Quantum Relative Entropy" 

386 symbStr = glyphs.qre(X.string, Y.string) 

387 

388 Expression.__init__(self, typeStr, symbStr) 

389 

390 # -------------------------------------------------------------------------- 

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

392 # -------------------------------------------------------------------------- 

393 

394 def _get_refined(self): 

395 if self._X.constant and (self._Y is None or self._Y.constant): 

396 return AffineExpression.from_constant(self.value, 1, self._symbStr) 

397 else: 

398 return self 

399 

400 Subtype = namedtuple("Subtype", ("argdim", "Y", "iscomplex")) 

401 

402 def _get_subtype(self): 

403 return self.Subtype(len(self._X), self._Y is not None, self._iscomplex) 

404 

405 def _get_value(self): 

406 X = cvx2np(self._X._get_value()) 

407 eigvalsX, eigvecsX = numpy.linalg.eigh(X) 

408 

409 if self._Y is None: 

410 eigvalsX = eigvalsX[eigvalsX > 1e-12] 

411 s = numpy.sum(eigvalsX * numpy.log(eigvalsX)) 

412 else: 

413 Y = eigvecsX.conj().T @ cvx2np(self._Y._get_value()) @ eigvecsX 

414 

415 Dy, Uy = numpy.linalg.eigh(Y) 

416 logY = Uy @ numpy.diag(numpy.log(Dy)) @ Uy.conj().T 

417 

418 s = -numpy.sum(numpy.diag(eigvalsX) * logY.conj()).real 

419 eigvalsX = eigvalsX[eigvalsX > 1e-12] 

420 s += numpy.sum(eigvalsX * numpy.log(eigvalsX)) 

421 

422 return cvxopt.matrix(s) 

423 

424 @cached_unary_operator 

425 def _get_mutables(self): 

426 if self._Y is None: 

427 return self._X._get_mutables() 

428 else: 

429 return self._X._get_mutables().union(self._Y.mutables) 

430 

431 def _is_convex(self): 

432 return True 

433 

434 def _is_concave(self): 

435 return False 

436 

437 def _replace_mutables(self, mapping): 

438 return self.__class__( 

439 self._X._replace_mutables(mapping), 

440 None if self._Y is None else self._Y._replace_mutables(mapping), 

441 ) 

442 

443 def _freeze_mutables(self, freeze): 

444 return self.__class__( 

445 self._X._freeze_mutables(freeze), 

446 None if self._Y is None else self._Y._freeze_mutables(freeze), 

447 ) 

448 

449 # -------------------------------------------------------------------------- 

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

451 # -------------------------------------------------------------------------- 

452 

453 @cached_selfinverse_unary_operator 

454 def __neg__(self): 

455 return QuantumEntropy(self._X, self._Y) 

456 

457 # -------------------------------------------------------------------------- 

458 # Methods and properties that return expressions. 

459 # -------------------------------------------------------------------------- 

460 

461 @property 

462 def X(self): 

463 """The expression :math:`X`.""" 

464 return self._X 

465 

466 @property 

467 def Y(self): 

468 """The additional expression :math:`Y`, or :obj:`None`.""" 

469 return self._Y 

470 

471 # -------------------------------------------------------------------------- 

472 # Methods and properties that describe the expression. 

473 # -------------------------------------------------------------------------- 

474 

475 @property 

476 def n(self): 

477 """Length of :attr:`X`.""" 

478 return len(self._X) 

479 

480 @property 

481 def iscomplex(self): 

482 """Whether :attr:`X` and :attr:`Y` are complex expressions or not.""" 

483 return self._iscomplex 

484 

485 # -------------------------------------------------------------------------- 

486 # Constraint-creating operators, and _predict. 

487 # -------------------------------------------------------------------------- 

488 

489 @classmethod 

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

491 assert isinstance(subtype, cls.Subtype) 

492 

493 if relation == operator.__le__: 

494 if ( 

495 issubclass(other.clstype, AffineExpression) 

496 and other.subtype.dim == 1 

497 ): 

498 if subtype.iscomplex: 

499 return ComplexQuantRelEntropyConstraint.make_type( 

500 argdim=subtype.argdim 

501 ) 

502 else: 

503 return QuantRelEntropyConstraint.make_type( 

504 argdim=subtype.argdim 

505 ) 

506 

507 return NotImplemented 

508 

509 @convert_operands(scalarRHS=True) 

510 @validate_prediction 

511 @refine_operands() 

512 def __le__(self, other): 

513 if isinstance(other, AffineExpression): 

514 if self.iscomplex: 

515 return ComplexQuantRelEntropyConstraint(self, other) 

516 else: 

517 return QuantRelEntropyConstraint(self, other) 

518 else: 

519 return NotImplemented 

520 

521 

522# -------------------------------------- 

523__all__ = api_end(_API_START, globals())