Coverage for picos/expressions/exp_quantcondentr.py: 88.80%

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:`QuantumConditionalEntropy`.""" 

20 

21import math 

22import operator 

23from collections import namedtuple 

24from functools import reduce 

25 

26import cvxopt 

27import numpy 

28 

29from .. import glyphs 

30from ..apidoc import api_end, api_start 

31from ..caching import cached_unary_operator 

32from ..constraints import ( 

33 QuantCondEntropyConstraint, 

34 ComplexQuantCondEntropyConstraint, 

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 QuantumConditionalEntropy(Expression): 

45 r"""Quantum conditional 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. Then this is defined as 

51 

52 .. math:: 

53 

54 -S(X) + S(\operatorname{Tr}_i(X)), 

55 

56 where :math:`S(X)=-\operatorname{Tr}(X\log(X))` is the quantum entropy, and 

57 :math:`\operatorname{Tr}_i` denotes the partial trace with respect to the 

58 :math:`i`-th subsystem. 

59 

60 .. warning:: 

61 

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

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

64 search. 

65 """ 

66 

67 # -------------------------------------------------------------------------- 

68 # Initialization and factory methods. 

69 # -------------------------------------------------------------------------- 

70 

71 @convert_and_refine_arguments("X") 

72 def __init__(self, X, subsystems, dimensions=2): 

73 r"""Construct an :class:`QuantumConditionalEntropy`. 

74 

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

76 :type X: ~picos.expressions.AffineExpression 

77 

78 :param subsystems: A collection of or a single subystem number, indexed 

79 from zero, corresponding to subsystems that shall be traced over. 

80 The value :math:`-1` refers to the last subsystem. 

81 :type subsystems: int or tuple or list 

82 

83 :param dimensions: Either an integer :math:`d` so that the subsystems 

84 are assumed to be all of shape :math:`d \times d`, or a sequence of 

85 subsystem shapes where an integer :math:`d` within the sequence is 

86 read as :math:`d \times d`. In any case, the elementwise product 

87 over all subsystem shapes must equal the expression's shape. 

88 :type dimensions: int or tuple or list 

89 """ 

90 # Check that X is an affine Hermitian expression 

91 if not isinstance(X, ComplexAffineExpression): 

92 raise TypeError( 

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

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

95 type(X).__name__ 

96 ) 

97 ) 

98 if not X.hermitian: 

99 raise TypeError( 

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

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

102 ) 

103 

104 self._X = X 

105 

106 self._iscomplex = not isinstance(X, AffineExpression) 

107 

108 # Check that subsystems and dimension are compatible with X. 

109 if isinstance(dimensions, int): 

110 dimensions = self._square_equal_subsystem_dims(dimensions) 

111 else: 

112 dimensions = [ 

113 (d, d) if isinstance(d, int) else d for d in dimensions 

114 ] 

115 

116 total_dim = reduce(lambda x, y: (x[0] * y[0], x[1] * y[1]), dimensions) 

117 if total_dim != X.shape: 

118 raise TypeError("Subsystem dimensions do not match expression.") 

119 

120 if isinstance(subsystems, int): 

121 subsystems = (subsystems,) 

122 

123 numSys = len(dimensions) 

124 subsystems = set(numSys - 1 if sys == -1 else sys for sys in subsystems) 

125 

126 for sys in subsystems: 

127 if not isinstance(sys, int): 

128 raise IndexError( 

129 "Subsystem indices must be integer, not {}.".format( 

130 type(sys).__name__ 

131 ) 

132 ) 

133 elif sys < 0: 

134 raise IndexError("Subsystem indices must be nonnegative.") 

135 elif sys >= numSys: 

136 raise IndexError( 

137 "Subsystem index {} out of range for {} " 

138 "systems total.".format(sys, numSys) 

139 ) 

140 elif dimensions[sys][0] != dimensions[sys][1]: 

141 raise TypeError( 

142 "Subsystem index {} refers to a non-square subsystem that " 

143 "cannot be traced over.".format(sys) 

144 ) 

145 

146 self._subsystems = subsystems 

147 self._dimensions = dimensions 

148 

149 sysStrings = None 

150 for sys in range(numSys): 

151 # Shape of current system. 

152 p, q = dimensions[sys] 

153 sysString = glyphs.matrix(glyphs.shape((p, q))) 

154 

155 # Only trace over selected subsystems. 

156 if sys not in subsystems: 

157 sysStrings = glyphs.kron(sysStrings, sysString) \ 

158 if sysStrings else sysString 

159 continue 

160 else: 

161 sysStrings = glyphs.kron(sysStrings, glyphs.trace(sysString)) \ 

162 if sysStrings else glyphs.trace(sysString) 

163 

164 typeStr = "Quantum Conditional Entropy" 

165 pxStr = glyphs.ptrace_(X.string, sysStrings) 

166 symbStr = glyphs.sub(glyphs.qe(X.string), glyphs.qe(pxStr)) 

167 

168 Expression.__init__(self, typeStr, symbStr) 

169 

170 def _square_equal_subsystem_dims(self, diagLen): 

171 m, n = self._X.shape 

172 k = math.log(m, diagLen) 

173 

174 if m != n or int(k) != k: 

175 raise TypeError( 

176 "The expression has shape {} so it cannot be " 

177 "decomposed into subsystems of shape {}.".format( 

178 glyphs.shape(self._X.shape), glyphs.shape((diagLen,) * 2) 

179 ) 

180 ) 

181 

182 return ((diagLen,) * 2,) * int(k) 

183 

184 # -------------------------------------------------------------------------- 

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

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

187 

188 def _get_refined(self): 

189 if self._X.constant: 

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

191 else: 

192 return self 

193 

194 Subtype = namedtuple("Subtype", ("argdim", "iscomplex")) 

195 

196 def _get_subtype(self): 

197 return self.Subtype(len(self._X), self._iscomplex) 

198 

199 def _get_value(self): 

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

201 eigX = numpy.linalg.eigvalsh(X) 

202 eigX = eigX[eigX > 1e-12] 

203 

204 pX = self._X.partial_trace(self.subsystems, self.dimensions) 

205 pX = cvx2np(pX._get_value()) 

206 eigpX = numpy.linalg.eigvalsh(pX) 

207 eigpX = eigpX[eigpX > 1e-12] 

208 

209 s = numpy.dot(eigpX, numpy.log(eigpX)) 

210 s -= numpy.dot(eigX, numpy.log(eigX)) 

211 

212 return cvxopt.matrix(s) 

213 

214 @cached_unary_operator 

215 def _get_mutables(self): 

216 return self._X._get_mutables() 

217 

218 def _is_convex(self): 

219 return False 

220 

221 def _is_concave(self): 

222 return True 

223 

224 def _replace_mutables(self, mapping): 

225 return self.__class__(self._X._replace_mutables(mapping)) 

226 

227 def _freeze_mutables(self, freeze): 

228 return self.__class__(self._X._freeze_mutables(freeze)) 

229 

230 # -------------------------------------------------------------------------- 

231 # Methods and properties that return expressions. 

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

233 

234 @property 

235 def X(self): 

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

237 return self._X 

238 

239 # -------------------------------------------------------------------------- 

240 # Methods and properties that describe the expression. 

241 # -------------------------------------------------------------------------- 

242 

243 @property 

244 def subsystems(self): 

245 """The subsystems being traced out of :math:`X`.""" 

246 return self._subsystems 

247 

248 @property 

249 def dimensions(self): 

250 """The dimensions of the subsystems of :math:`X`.""" 

251 return self._dimensions 

252 

253 @property 

254 def n(self): 

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

256 return len(self._X) 

257 

258 @property 

259 def iscomplex(self): 

260 """Whether :attr:`X` is a complex expression or not.""" 

261 return self._iscomplex 

262 

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

264 # Constraint-creating operators, and _predict. 

265 # -------------------------------------------------------------------------- 

266 

267 @classmethod 

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

269 assert isinstance(subtype, cls.Subtype) 

270 

271 if relation == operator.__ge__: 

272 if ( 

273 issubclass(other.clstype, AffineExpression) 

274 and other.subtype.dim == 1 

275 ): 

276 if subtype.iscomplex: 

277 return ComplexQuantCondEntropyConstraint.make_type( 

278 argdim=subtype.argdim 

279 ) 

280 else: 

281 return QuantCondEntropyConstraint.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 if self.iscomplex: 

292 return ComplexQuantCondEntropyConstraint(self, other) 

293 else: 

294 return QuantCondEntropyConstraint(self, other) 

295 else: 

296 return NotImplemented 

297 

298 

299# -------------------------------------- 

300__all__ = api_end(_API_START, globals())