Coverage for qml_essentials/entanglement.py: 79%

1from typing import Optional, Any, List 

2import pennylane as qml 

3import pennylane.numpy as np 

4from copy import deepcopy 

5from qml_essentials.utils import logm_v 

6from qml_essentials.model import Model 

7import logging 

8 

9log = logging.getLogger(__name__) 

10 

11 

12class Entanglement: 

13 

14 @staticmethod 

15 def meyer_wallach( 

16 model: Model, 

17 n_samples: Optional[int | None], 

18 seed: Optional[int], 

19 scale: bool = False, 

20 **kwargs: Any, 

21 ) -> float: 

22 """ 

23 Calculates the entangling capacity of a given quantum circuit 

24 using Meyer-Wallach measure. 

25 

26 Args: 

27 model (Model): The quantum circuit model. 

28 n_samples (Optional[int]): Number of samples per qubit. 

29 If None or < 0, the current parameters of the model are used. 

30 seed (Optional[int]): Seed for the random number generator. 

31 scale (bool): Whether to scale the number of samples. 

32 kwargs (Any): Additional keyword arguments for the model function. 

33 

34 Returns: 

35 float: Entangling capacity of the given circuit, guaranteed 

36 to be between 0.0 and 1.0. 

37 """ 

38 if "noise_params" in kwargs: 

39 log.warning( 

40 "Meyer-Wallach measure not suitable for noisy circuits.\ 

41 Consider 'relative_entropy' instead." 

42 ) 

43 

44 if scale: 

45 n_samples = np.power(2, model.n_qubits) * n_samples 

46 

47 rng = np.random.default_rng(seed) 

48 if n_samples is not None and n_samples > 0: 

49 assert seed is not None, "Seed must be provided when samples > 0" 

50 # TODO: maybe switch to JAX rng 

51 model.initialize_params(rng=rng, repeat=n_samples) 

52 else: 

53 if seed is not None: 

54 log.warning("Seed is ignored when samples is 0") 

55 

56 # implicitly set input to none in case it's not needed 

57 kwargs.setdefault("inputs", None) 

58 # explicitly set execution type because everything else won't work 

59 rhos = model(execution_type="density", **kwargs).reshape( 

60 -1, 2**model.n_qubits, 2**model.n_qubits 

61 ) 

62 

63 mw_measure = np.zeros(len(rhos)) 

64 

65 for i, rho in enumerate(rhos): 

66 mw_measure[i] = Entanglement._compute_meyer_wallach_meas( 

67 rho, model.n_qubits 

68 ) 

69 

70 # Average all iterated states 

71 entangling_capability = min(max(mw_measure.mean(), 0.0), 1.0) 

72 log.debug(f"Variance of measure: {mw_measure.var()}") 

73 

74 # catch floating point errors 

75 return float(entangling_capability) 

76 

77 @staticmethod 

78 def _compute_meyer_wallach_meas(rho: np.ndarray, n_qubits: int): 

79 qb = list(range(n_qubits)) 

80 entropy = 0 

81 for j in range(n_qubits): 

82 # Formula 6 in https://doi.org/10.48550/arXiv.quant-ph/0305094 

83 density = qml.math.partial_trace(rho, qb[:j] + qb[j + 1 :]) 

84 # only real values, because imaginary part will be separate 

85 # in all following calculations anyway 

86 # entropy should be 1/2 <= entropy <= 1 

87 entropy += np.trace((density @ density).real) 

88 

89 # inverse averaged entropy and scale to [0, 1] 

90 return 2 * (1 - entropy / n_qubits) 

91 

92 @staticmethod 

93 def bell_measurements( 

94 model: Model, n_samples: int, seed: int, scale: bool = False, **kwargs: Any 

95 ) -> float: 

96 """ 

97 Compute the Bell measurement for a given model. 

98 

99 Args: 

100 model (Model): The quantum circuit model. 

101 n_samples (int): The number of samples to compute the measure for. 

102 seed (int): The seed for the random number generator. 

103 scale (bool): Whether to scale the number of samples 

104 according to the number of qubits. 

105 **kwargs (Any): Additional keyword arguments for the model function. 

106 

107 Returns: 

108 float: The Bell measurement value. 

109 """ 

110 if "noise_params" in kwargs: 

111 log.warning( 

112 "Bell Measurements not suitable for noisy circuits.\ 

113 Consider 'relative_entropy' instead." 

114 ) 

115 

116 if scale: 

117 n_samples = np.power(2, model.n_qubits) * n_samples 

118 

119 def _circuit(params: np.ndarray, inputs: np.ndarray) -> List[np.ndarray]: 

120 """ 

121 Compute the Bell measurement circuit. 

122 

123 Args: 

124 params (np.ndarray): The model parameters. 

125 inputs (np.ndarray): The input to the model. 

126 

127 Returns: 

128 List[np.ndarray]: The probabilities of the Bell measurement. 

129 """ 

130 model._variational(params, inputs) 

131 

132 qml.map_wires( 

133 model._variational, 

134 {i: i + model.n_qubits for i in range(model.n_qubits)}, 

135 )(params, inputs) 

136 

137 for q in range(model.n_qubits): 

138 qml.CNOT(wires=[q, q + model.n_qubits]) 

139 qml.H(q) 

140 

141 obs_wires = [(q, q + model.n_qubits) for q in range(model.n_qubits)] 

142 return [qml.probs(wires=w) for w in obs_wires] 

143 

144 model.circuit = qml.QNode( 

145 _circuit, 

146 qml.device( 

147 "default.qubit", 

148 shots=model.shots, 

149 wires=model.n_qubits * 2, 

150 ), 

151 ) 

152 

153 rng = np.random.default_rng(seed) 

154 if n_samples is not None and n_samples > 0: 

155 assert seed is not None, "Seed must be provided when samples > 0" 

156 # TODO: maybe switch to JAX rng 

157 model.initialize_params(rng=rng, repeat=n_samples) 

158 params = model.params 

159 else: 

160 if seed is not None: 

161 log.warning("Seed is ignored when samples is 0") 

162 

163 if len(model.params.shape) <= 2: 

164 params = model.params.reshape(*model.params.shape, 1) 

165 else: 

166 log.info(f"Using sample size of model params: {model.params.shape[-1]}") 

167 params = model.params 

168 

169 n_samples = params.shape[-1] 

170 mw_measure = np.zeros(n_samples) 

171 

172 # implicitly set input to none in case it's not needed 

173 kwargs.setdefault("inputs", None) 

174 exp = model(params=params, **kwargs) 

175 exp = 1 - 2 * exp[:, :, -1] 

176 mw_measure = 2 * (1 - exp.mean(axis=0)) 

177 entangling_capability = min(max(mw_measure.mean(), 0.0), 1.0) 

178 log.debug(f"Variance of measure: {mw_measure.var()}") 

179 

180 return float(entangling_capability) 

181 

182 @staticmethod 

183 def relative_entropy( 

184 model: Model, 

185 n_samples: int, 

186 n_sigmas: int, 

187 seed: Optional[int], 

188 scale: bool = False, 

189 **kwargs: Any, 

190 ) -> float: 

191 """ 

192 Calculates the relative entropy of entanglement of a given quantum 

193 circuit. This measure is also applicable to mixed state, albeit it 

194 might me not fully accurate in this simplified case. 

195 

196 As the relative entropy is generally defined as the smallest relative 

197 entropy from the state in question to the set of separable states. 

198 However, as computing the nearest separable state is NP-hard, we select 

199 n_sigmas of random separable states to compute the distance to, which 

200 is not necessarily the nearest. Thus, this measure of entanglement 

201 presents an upper limit of entanglement. 

202 

203 As the relative entropy is not necessarily between zero and one, this 

204 function also normalises by the relative entroy to the GHZ state. 

205 

206 Args: 

207 model (Model): The quantum circuit model. 

208 n_samples (int): Number of samples per qubit. 

209 If <= 0, the current parameters of the model are used. 

210 n_sigmas (int): Number of random separable pure states to compare against. 

211 seed (Optional[int]): Seed for the random number generator. 

212 scale (bool): Whether to scale the number of samples. 

213 kwargs (Any): Additional keyword arguments for the model function. 

214 

215 Returns: 

216 float: Entangling capacity of the given circuit, guaranteed 

217 to be between 0.0 and 1.0. 

218 """ 

219 dim = np.power(2, model.n_qubits) 

220 if scale: 

221 n_samples = dim * n_samples 

222 n_sigmas = dim * n_sigmas 

223 

224 rng = np.random.default_rng(seed) 

225 

226 # Random separable states 

227 log_sigmas = sample_random_separable_states( 

228 model.n_qubits, n_samples=n_sigmas, rng=rng, take_log=True 

229 ) 

230 

231 if n_samples > 0: 

232 assert seed is not None, "Seed must be provided when samples > 0" 

233 model.initialize_params(rng=rng, repeat=n_samples) 

234 else: 

235 if seed is not None: 

236 log.warning("Seed is ignored when samples is 0") 

237 

238 if len(model.params.shape) <= 2: 

239 model.params = model.params.reshape(*model.params.shape, 1) 

240 else: 

241 log.info(f"Using sample size of model params: {model.params.shape[-1]}") 

242 

243 ghz_model = Model(model.n_qubits, 1, "GHZ", data_reupload=False) 

244 

245 normalised_entropies = np.zeros((n_sigmas, n_samples)) 

246 for j, log_sigma in enumerate(log_sigmas): 

247 

248 # Entropy of GHZ states should be maximal 

249 ghz_entropy = Entanglement._compute_rel_entropies( 

250 ghz_model, 

251 log_sigma, 

252 ) 

253 

254 rel_entropy = Entanglement._compute_rel_entropies( 

255 model, log_sigma, **kwargs 

256 ) 

257 

258 normalised_entropies[j] = rel_entropy / ghz_entropy 

259 

260 # Average all iterated states 

261 entangling_capability = normalised_entropies.min(axis=0).mean() 

262 log.debug(f"Variance of measure: {normalised_entropies.var()}") 

263 

264 return entangling_capability 

265 

266 @staticmethod 

267 def _compute_rel_entropies( 

268 model: Model, 

269 log_sigma: np.ndarray, 

270 **kwargs, 

271 ) -> np.ndarray: 

272 """ 

273 Compute the relative entropy for a given model. 

274 

275 Args: 

276 model (Model): The model for which to compute entanglement 

277 log_sigma (np.ndarray): Density matrix of next separable state 

278 

279 Returns: 

280 np.ndarray: Relative Entropy for each sample 

281 """ 

282 # implicitly set input to none in case it's not needed 

283 kwargs.setdefault("inputs", None) 

284 # explicitly set execution type because everything else won't work 

285 rho = model(execution_type="density", **kwargs) 

286 rho = rho.reshape(-1, 2**model.n_qubits, 2**model.n_qubits) 

287 log_rho = logm_v(rho) / np.log(2) 

288 

289 rel_entropies = np.abs(np.trace(rho @ (log_rho - log_sigma), axis1=1, axis2=2)) 

290 

291 return rel_entropies 

292 

293 @staticmethod 

294 def entanglement_of_formation( 

295 model: Model, 

296 n_samples: int, 

297 seed: Optional[int], 

298 scale: bool = False, 

299 always_decompose: bool = False, 

300 **kwargs: Any, 

301 ) -> float: 

302 """ 

303 This function implements the entanglement of formation for mixed 

304 quantum systems. 

305 In that a mixed state gets decomposed into pure states with respective 

306 probabilities using the eigendecomposition of the density matrix. 

307 Then, the Meyer-Wallach measure is computed for each pure state, 

308 weighted by the eigenvalue. 

309 See e.g. https://doi.org/10.48550/arXiv.quant-ph/0504163 

310 

311 Note that the decomposition is *not unique*! Therefore, this measure 

312 presents the entanglement for *some* decomposition into pure states, 

313 not necessarily the one that is anticipated when applying the Kraus 

314 channels. 

315 If a pure state is provided, this results in the same value as the 

316 Entanglement.meyer_wallach function if `always_decompose` flag is not set. 

317 

318 Args: 

319 model (Model): The quantum circuit model. 

320 n_samples (int): Number of samples per qubit. 

321 seed (Optional[int]): Seed for the random number generator. 

322 scale (bool): Whether to scale the number of samples. 

323 always_decompose (bool): Whether to explicitly compute the 

324 entantlement of formation for the eigendecomposition of a pure 

325 state. 

326 kwargs (Any): Additional keyword arguments for the model function. 

327 

328 Returns: 

329 float: Entangling capacity of the given circuit, guaranteed 

330 to be between 0.0 and 1.0. 

331 """ 

332 

333 if scale: 

334 n_samples = np.power(2, model.n_qubits) * n_samples 

335 

336 rng = np.random.default_rng(seed) 

337 if n_samples > 0: 

338 assert seed is not None, "Seed must be provided when samples > 0" 

339 model.initialize_params(rng=rng, repeat=n_samples) 

340 else: 

341 if seed is not None: 

342 log.warning("Seed is ignored when samples is 0") 

343 

344 if len(model.params.shape) <= 2: 

345 model.params = model.params.reshape(*model.params.shape, 1) 

346 else: 

347 log.info(f"Using sample size of model params: {model.params.shape[-1]}") 

348 

349 # implicitly set input to none in case it's not needed 

350 kwargs.setdefault("inputs", None) 

351 rhos = model(execution_type="density", **kwargs) 

352 rhos = rhos.reshape(-1, 2**model.n_qubits, 2**model.n_qubits) 

353 entanglement = np.zeros(len(rhos)) 

354 for i, rho in enumerate(rhos): 

355 entanglement[i] = Entanglement._compute_entanglement_of_formation( 

356 rho, model.n_qubits, always_decompose 

357 ) 

358 entangling_capability = min(max(entanglement.mean(), 0.0), 1.0) 

359 return float(entangling_capability) 

360 

361 @staticmethod 

362 def _compute_entanglement_of_formation( 

363 rho: np.ndarray, n_qubits: int, always_decompose: bool 

364 ) -> float: 

365 """ 

366 Computes the entanglement of formation for a given density matrix rho. 

367 

368 Args: 

369 rho (np.ndarray): The density matrix 

370 n_qubits (int): Number of qubits 

371 always_decompose (bool): Whether to explicitly compute the 

372 entantlement of formation for the eigendecomposition of a pure 

373 state. 

374 

375 Returns: 

376 float: Entanglement for the provided state. 

377 """ 

378 eigenvalues, eigenvectors = np.linalg.eigh(rho) 

379 if any(np.isclose(eigenvalues, 1.0)) and not always_decompose: # Pure state 

380 return Entanglement._compute_meyer_wallach_meas(rho, n_qubits) 

381 ent = 0 

382 for prob, ev in zip(eigenvalues, eigenvectors): 

383 ev = ev.reshape(-1, 1) 

384 rho = ev @ np.conjugate(ev).T 

385 mw_measure = Entanglement._compute_meyer_wallach_meas(rho, n_qubits) 

386 ent += prob * mw_measure 

387 return ent 

388 

389 

390def sample_random_separable_states( 

391 n_qubits: int, n_samples: int, rng: np.random.Generator, take_log: bool = False 

392) -> np.ndarray: 

393 """ 

394 Sample random separable states (density matrix). 

395 

396 Args: 

397 n_qubits (int): number of qubits in the state 

398 n_samples (int): number of states 

399 rng (np.random.Generator): random number generator 

400 take_log (bool): if the matrix logarithm of the density matrix should be taken. 

401 

402 Returns: 

403 np.ndarray: Density matrices of shape (n_samples, 2**n_qubits, 2**n_qubits) 

404 """ 

405 model = Model(n_qubits, 1, "No_Entangling", data_reupload=False) 

406 model.initialize_params(rng=rng, repeat=n_samples) 

407 # explicitly set execution type because everything else won't work 

408 sigmas = model(execution_type="density", inputs=None) 

409 if take_log: 

410 sigmas = logm_v(sigmas) / np.log(2) 

411 

412 return sigmas