Coverage for qml_essentials/entanglement.py: 90%

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 measure = np.zeros(len(rhos)) 

64 

65 for i, rho in enumerate(rhos): 

66 measure[i] = Entanglement._compute_meyer_wallach_meas(rho, model.n_qubits) 

67 

68 # Average all iterated states 

69 entangling_capability = min(max(measure.mean(), 0.0), 1.0) 

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

71 

72 # catch floating point errors 

73 return float(entangling_capability) 

74 

75 @staticmethod 

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

77 qb = list(range(n_qubits)) 

78 entropy = 0 

79 for j in range(n_qubits): 

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

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

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

83 # in all following calculations anyway 

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

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

86 

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

88 return 2 * (1 - entropy / n_qubits) 

89 

90 @staticmethod 

91 def bell_measurements( 

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

93 ) -> float: 

94 """ 

95 Compute the Bell measurement for a given model. 

96 

97 Args: 

98 model (Model): The quantum circuit model. 

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

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

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

102 according to the number of qubits. 

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

104 

105 Returns: 

106 float: The Bell measurement value. 

107 """ 

108 if "noise_params" in kwargs: 

109 log.warning( 

110 "Bell Measurements not suitable for noisy circuits.\ 

111 Consider 'relative_entropy' instead." 

112 ) 

113 

114 if scale: 

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

116 

117 def _circuit( 

118 params: np.ndarray, 

119 inputs: np.ndarray, 

120 enc_params: Optional[np.ndarray] = None, 

121 ) -> List[np.ndarray]: 

122 """ 

123 Compute the Bell measurement circuit. 

124 

125 Args: 

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

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

128 enc_params (Optional[np.ndarray]): The frequency encoding parameters. 

129 

130 Returns: 

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

132 """ 

133 model._variational(params, inputs, enc_params) 

134 

135 qml.map_wires( 

136 model._variational, 

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

138 )(params, inputs) 

139 

140 for q in range(model.n_qubits): 

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

142 qml.H(q) 

143 

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

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

146 

147 model.circuit = qml.QNode( 

148 _circuit, 

149 qml.device( 

150 "default.qubit", 

151 shots=model.shots, 

152 wires=model.n_qubits * 2, 

153 ), 

154 ) 

155 

156 rng = np.random.default_rng(seed) 

157 if n_samples is not None and n_samples > 0: 

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

159 # TODO: maybe switch to JAX rng 

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

161 params = model.params 

162 else: 

163 if seed is not None: 

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

165 

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

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

168 else: 

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

170 params = model.params 

171 

172 n_samples = params.shape[-1] 

173 measure = np.zeros(n_samples) 

174 

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

176 kwargs.setdefault("inputs", None) 

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

178 exp = 1 - 2 * exp[..., -1] 

179 measure = 2 * (1 - exp.mean(axis=0)) 

180 entangling_capability = min(max(measure.mean(), 0.0), 1.0) 

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

182 

183 return float(entangling_capability) 

184 

185 @staticmethod 

186 def relative_entropy( 

187 model: Model, 

188 n_samples: int, 

189 n_sigmas: int, 

190 seed: Optional[int], 

191 scale: bool = False, 

192 **kwargs: Any, 

193 ) -> float: 

194 """ 

195 Calculates the relative entropy of entanglement of a given quantum 

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

197 might me not fully accurate in this simplified case. 

198 

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

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

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

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

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

204 presents an upper limit of entanglement. 

205 

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

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

208 

209 Args: 

210 model (Model): The quantum circuit model. 

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

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

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

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

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

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

217 

218 Returns: 

219 float: Entangling capacity of the given circuit, guaranteed 

220 to be between 0.0 and 1.0. 

221 """ 

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

223 if scale: 

224 n_samples = dim * n_samples 

225 n_sigmas = dim * n_sigmas 

226 

227 rng = np.random.default_rng(seed) 

228 

229 # Random separable states 

230 log_sigmas = sample_random_separable_states( 

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

232 ) 

233 

234 if n_samples is not None and n_samples > 0: 

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

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

237 else: 

238 if seed is not None: 

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

240 

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

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

243 else: 

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

245 

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

247 

248 normalised_entropies = np.zeros((n_sigmas, model.params.shape[-1])) 

249 for j, log_sigma in enumerate(log_sigmas): 

250 

251 # Entropy of GHZ states should be maximal 

252 ghz_entropy = Entanglement._compute_rel_entropies( 

253 ghz_model, 

254 log_sigma, 

255 ) 

256 

257 rel_entropy = Entanglement._compute_rel_entropies( 

258 model, log_sigma, **kwargs 

259 ) 

260 

261 normalised_entropies[j] = rel_entropy / ghz_entropy 

262 

263 # Average all iterated states 

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

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

266 

267 return entangling_capability 

268 

269 @staticmethod 

270 def _compute_rel_entropies( 

271 model: Model, 

272 log_sigma: np.ndarray, 

273 **kwargs, 

274 ) -> np.ndarray: 

275 """ 

276 Compute the relative entropy for a given model. 

277 

278 Args: 

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

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

281 

282 Returns: 

283 np.ndarray: Relative Entropy for each sample 

284 """ 

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

286 kwargs.setdefault("inputs", None) 

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

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

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

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

291 

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

293 

294 return rel_entropies 

295 

296 @staticmethod 

297 def entanglement_of_formation( 

298 model: Model, 

299 n_samples: int, 

300 seed: Optional[int], 

301 scale: bool = False, 

302 always_decompose: bool = False, 

303 **kwargs: Any, 

304 ) -> float: 

305 """ 

306 This function implements the entanglement of formation for mixed 

307 quantum systems. 

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

309 probabilities using the eigendecomposition of the density matrix. 

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

311 weighted by the eigenvalue. 

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

313 

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

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

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

317 channels. 

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

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

320 

321 Args: 

322 model (Model): The quantum circuit model. 

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

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

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

326 always_decompose (bool): Whether to explicitly compute the 

327 entantlement of formation for the eigendecomposition of a pure 

328 state. 

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

330 

331 Returns: 

332 float: Entangling capacity of the given circuit, guaranteed 

333 to be between 0.0 and 1.0. 

334 """ 

335 

336 if scale: 

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

338 

339 rng = np.random.default_rng(seed) 

340 if n_samples is not None and n_samples > 0: 

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

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

343 else: 

344 if seed is not None: 

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

346 

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

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

349 else: 

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

351 

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

353 kwargs.setdefault("inputs", None) 

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

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

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

357 for i, rho in enumerate(rhos): 

358 entanglement[i] = Entanglement._compute_entanglement_of_formation( 

359 rho, model.n_qubits, always_decompose 

360 ) 

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

362 return float(entangling_capability) 

363 

364 @staticmethod 

365 def _compute_entanglement_of_formation( 

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

367 ) -> float: 

368 """ 

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

370 

371 Args: 

372 rho (np.ndarray): The density matrix 

373 n_qubits (int): Number of qubits 

374 always_decompose (bool): Whether to explicitly compute the 

375 entantlement of formation for the eigendecomposition of a pure 

376 state. 

377 

378 Returns: 

379 float: Entanglement for the provided state. 

380 """ 

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

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

383 return Entanglement._compute_meyer_wallach_meas(rho, n_qubits) 

384 ent = 0 

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

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

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

388 measure = Entanglement._compute_meyer_wallach_meas(rho, n_qubits) 

389 ent += prob * measure 

390 return ent 

391 

392 @staticmethod 

393 def concentratable_entanglement( 

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

395 ) -> float: 

396 """ 

397 Computes the concentratable entanglement of a given model. 

398 

399 This method utilizes the Concentratable Entanglement measure from 

400 https://arxiv.org/abs/2104.06923. 

401 

402 Args: 

403 model (Model): The quantum circuit model. 

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

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

406 scale (bool): Whether to scale the number of samples according to 

407 the number of qubits. 

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

409 

410 Returns: 

411 float: Entangling capability of the given circuit, guaranteed 

412 to be between 0.0 and 1.0. 

413 """ 

414 if "noise_params" in kwargs: 

415 log.warning( 

416 "Concentratable entanglement is not suitable for noisy circuits.\ 

417 Consider 'relative_entropy' instead." 

418 ) 

419 

420 n = model.n_qubits 

421 

422 if scale: 

423 n_samples = np.power(2, model.n) * n_samples 

424 

425 def _circuit( 

426 params: np.ndarray, 

427 inputs: np.ndarray, 

428 enc_params: Optional[np.ndarray] = None, 

429 ) -> List[np.ndarray]: 

430 """ 

431 Constructs a circuit to compute the concentratable entanglement using the 

432 swap test by creating two copies of the models circuit and map the output 

433 wires accordingly 

434 

435 Args: 

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

437 inputs (np.ndarray): The input data for the model. 

438 enc_params (Optional[np.ndarray]): Optional encoding parameters. 

439 

440 Returns: 

441 List[np.ndarray]: Probabilities obtained from the swap test circuit. 

442 """ 

443 

444 qml.map_wires(model._variational, {i: i + n for i in range(n)})( 

445 params, inputs, enc_params 

446 ) 

447 qml.map_wires(model._variational, {i: i + 2 * n for i in range(n)})( 

448 params, inputs, enc_params 

449 ) 

450 

451 # Perform swap test 

452 for i in range(n): 

453 qml.H(i) 

454 

455 for i in range(n): 

456 qml.CSWAP([i, i + n, i + 2 * n]) 

457 

458 for i in range(n): 

459 qml.H(i) 

460 

461 return qml.probs(wires=[i for i in range(n)]) 

462 

463 model.circuit = qml.QNode( 

464 _circuit, 

465 qml.device( 

466 "default.qubit", 

467 shots=model.shots, 

468 wires=n * 3, 

469 ), 

470 ) 

471 

472 rng = np.random.default_rng(seed) 

473 if n_samples is not None and n_samples > 0: 

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

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

476 params = model.params 

477 else: 

478 if seed is not None: 

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

480 

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

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

483 else: 

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

485 params = model.params 

486 

487 n_samples = params.shape[-1] 

488 

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

490 kwargs.setdefault("inputs", None) 

491 

492 samples_probs = model(params=params, execution_type="probs", **kwargs) 

493 if n_samples == 1: 

494 samples_probs = [samples_probs] 

495 

496 ce_measure = np.zeros(len(samples_probs)) 

497 

498 for i, probs in enumerate(samples_probs): 

499 ce_measure[i] = 1 - probs[0] 

500 

501 # Average all iterated states 

502 entangling_capability = min(max(ce_measure.mean(), 0.0), 1.0) 

503 log.debug(f"Variance of measure: {ce_measure.var()}") 

504 

505 # catch floating point errors 

506 return float(entangling_capability) 

507 

508 

509def sample_random_separable_states( 

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

511) -> np.ndarray: 

512 """ 

513 Sample random separable states (density matrix). 

514 

515 Args: 

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

517 n_samples (int): number of states 

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

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

520 

521 Returns: 

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

523 """ 

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

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

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

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

528 if take_log: 

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

530 

531 return sigmas