Coverage for qml_essentials/entanglement.py: 79%

145 statements  

« prev     ^ index     » next       coverage.py v7.8.0, created at 2025-06-10 13:05 +0000

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( 

120 params: np.ndarray, inputs: np.ndarray, 

121 enc_params: Optional[np.ndarray] = None 

122 ) -> List[np.ndarray]: 

123 """ 

124 Compute the Bell measurement circuit. 

125 

126 Args: 

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

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

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

130 

131 Returns: 

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

133 """ 

134 model._variational(params, inputs, enc_params) 

135 

136 qml.map_wires( 

137 model._variational, 

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

139 )(params, inputs) 

140 

141 for q in range(model.n_qubits): 

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

143 qml.H(q) 

144 

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

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

147 

148 model.circuit = qml.QNode( 

149 _circuit, 

150 qml.device( 

151 "default.qubit", 

152 shots=model.shots, 

153 wires=model.n_qubits * 2, 

154 ), 

155 ) 

156 

157 rng = np.random.default_rng(seed) 

158 if n_samples is not None and n_samples > 0: 

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

160 # TODO: maybe switch to JAX rng 

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

162 params = model.params 

163 else: 

164 if seed is not None: 

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

166 

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

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

169 else: 

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

171 params = model.params 

172 

173 n_samples = params.shape[-1] 

174 mw_measure = np.zeros(n_samples) 

175 

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

177 kwargs.setdefault("inputs", None) 

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

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

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

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

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

183 

184 return float(entangling_capability) 

185 

186 @staticmethod 

187 def relative_entropy( 

188 model: Model, 

189 n_samples: int, 

190 n_sigmas: int, 

191 seed: Optional[int], 

192 scale: bool = False, 

193 **kwargs: Any, 

194 ) -> float: 

195 """ 

196 Calculates the relative entropy of entanglement of a given quantum 

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

198 might me not fully accurate in this simplified case. 

199 

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

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

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

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

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

205 presents an upper limit of entanglement. 

206 

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

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

209 

210 Args: 

211 model (Model): The quantum circuit model. 

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

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

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

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

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

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

218 

219 Returns: 

220 float: Entangling capacity of the given circuit, guaranteed 

221 to be between 0.0 and 1.0. 

222 """ 

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

224 if scale: 

225 n_samples = dim * n_samples 

226 n_sigmas = dim * n_sigmas 

227 

228 rng = np.random.default_rng(seed) 

229 

230 # Random separable states 

231 log_sigmas = sample_random_separable_states( 

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

233 ) 

234 

235 if n_samples > 0: 

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

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

238 else: 

239 if seed is not None: 

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

241 

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

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

244 else: 

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

246 

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

248 

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

250 for j, log_sigma in enumerate(log_sigmas): 

251 

252 # Entropy of GHZ states should be maximal 

253 ghz_entropy = Entanglement._compute_rel_entropies( 

254 ghz_model, 

255 log_sigma, 

256 ) 

257 

258 rel_entropy = Entanglement._compute_rel_entropies( 

259 model, log_sigma, **kwargs 

260 ) 

261 

262 normalised_entropies[j] = rel_entropy / ghz_entropy 

263 

264 # Average all iterated states 

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

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

267 

268 return entangling_capability 

269 

270 @staticmethod 

271 def _compute_rel_entropies( 

272 model: Model, 

273 log_sigma: np.ndarray, 

274 **kwargs, 

275 ) -> np.ndarray: 

276 """ 

277 Compute the relative entropy for a given model. 

278 

279 Args: 

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

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

282 

283 Returns: 

284 np.ndarray: Relative Entropy for each sample 

285 """ 

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

287 kwargs.setdefault("inputs", None) 

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

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

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

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

292 

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

294 

295 return rel_entropies 

296 

297 @staticmethod 

298 def entanglement_of_formation( 

299 model: Model, 

300 n_samples: int, 

301 seed: Optional[int], 

302 scale: bool = False, 

303 always_decompose: bool = False, 

304 **kwargs: Any, 

305 ) -> float: 

306 """ 

307 This function implements the entanglement of formation for mixed 

308 quantum systems. 

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

310 probabilities using the eigendecomposition of the density matrix. 

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

312 weighted by the eigenvalue. 

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

314 

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

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

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

318 channels. 

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

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

321 

322 Args: 

323 model (Model): The quantum circuit model. 

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

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

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

327 always_decompose (bool): Whether to explicitly compute the 

328 entantlement of formation for the eigendecomposition of a pure 

329 state. 

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

331 

332 Returns: 

333 float: Entangling capacity of the given circuit, guaranteed 

334 to be between 0.0 and 1.0. 

335 """ 

336 

337 if scale: 

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

339 

340 rng = np.random.default_rng(seed) 

341 if n_samples > 0: 

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

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

344 else: 

345 if seed is not None: 

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

347 

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

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

350 else: 

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

352 

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

354 kwargs.setdefault("inputs", None) 

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

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

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

358 for i, rho in enumerate(rhos): 

359 entanglement[i] = Entanglement._compute_entanglement_of_formation( 

360 rho, model.n_qubits, always_decompose 

361 ) 

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

363 return float(entangling_capability) 

364 

365 @staticmethod 

366 def _compute_entanglement_of_formation( 

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

368 ) -> float: 

369 """ 

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

371 

372 Args: 

373 rho (np.ndarray): The density matrix 

374 n_qubits (int): Number of qubits 

375 always_decompose (bool): Whether to explicitly compute the 

376 entantlement of formation for the eigendecomposition of a pure 

377 state. 

378 

379 Returns: 

380 float: Entanglement for the provided state. 

381 """ 

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

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

384 return Entanglement._compute_meyer_wallach_meas(rho, n_qubits) 

385 ent = 0 

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

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

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

389 mw_measure = Entanglement._compute_meyer_wallach_meas(rho, n_qubits) 

390 ent += prob * mw_measure 

391 return ent 

392 

393 

394def sample_random_separable_states( 

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

396) -> np.ndarray: 

397 """ 

398 Sample random separable states (density matrix). 

399 

400 Args: 

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

402 n_samples (int): number of states 

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

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

405 

406 Returns: 

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

408 """ 

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

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

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

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

413 if take_log: 

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

415 

416 return sigmas