Coverage for qml_essentials / algebra.py: 98%

1r"""Dynamical Lie algebra (DLA) closure helpers. 

2 

3The DLA :math:`\mathfrak{g} = \langle i G \rangle_{\mathrm{Lie}}` of a set of 

4Hermitian generators :math:`G` is the real Lie algebra spanned by repeated 

5brackets. Two representations are provided: 

6 

7- :func:`lie_closure_paulis` for generators given as Pauli words or strings. 

8 Every bracket :math:`[i P_a, i P_b]` is again proportional to a single Pauli 

9 string, so the closure is the set of Pauli strings reached by composing 

10 anticommuting pairs. 

11- :func:`lie_closure_matrices` for generators given as dense Hermitian matrices 

12 (for example sums of Paulis), returning a Hilbert-Schmidt-orthonormal 

13 Hermitian basis. 

14 

15In both cases the length of the returned basis equals 

16:math:`\dim \mathfrak{g}`. 

17 

18The matchgate algebra :math:`\mathfrak{so}(2n)` is also provided explicitly by 

19:func:`matchgate_generators`, :func:`matchgate_basis`, and :func:`dim_so2n`. 

20The g-purity of a state with respect to a DLA basis is given by 

21:func:`g_purity_from_basis` (Pauli-word basis) and :func:`g_purity_matrix` 

22(Hilbert-Schmidt-orthonormal matrix basis). 

23 

24Permutation-symmetric operators are built by :func:`symmetric_pauli_sum`, with 

25:func:`sn_equivariant_generators` and :func:`sn_equivariant_observable` giving the 

26generators and observable of the :math:`S_n`-equivariant ansatz; the generators 

27feed the matrix DLA via :func:`lie_closure_matrices`. 

28""" 

29 

30from itertools import combinations 

31from typing import List, Optional, Sequence, Union 

32 

33import numpy as np 

34 

35from jaqsi.paulis import PauliWord, state_expectation 

36 

37 

38def lie_closure_paulis( 

39 generators: Sequence[Union[str, PauliWord]], 

40 max_dim: Optional[int] = None, 

41) -> List[PauliWord]: 

42 r"""Return the Pauli-word spanning set of the DLA generated by *generators*. 

43 

44 Each generator is a Hermitian Pauli word :math:`P` (standing for 

45 :math:`i P` in :math:`\mathfrak{g}`). The closure is grown by repeatedly 

46 composing anticommuting pairs; every such product is again a single Pauli 

47 word. Deduplication is on the bare Pauli string (the global phase does not 

48 affect the span), and the returned words carry that bare string. The 

49 length of the result is :math:`\dim \mathfrak{g}`. 

50 

51 Args: 

52 generators: Hermitian generators, each a :class:`PauliWord` or a bare 

53 Pauli string over ``{'I', 'X', 'Y', 'Z'}`` (qubit 0 leftmost). 

54 max_dim: Optional cap (at least 1) on the number of words. Growth stops 

55 as soon as the basis holds ``max_dim`` words, which keeps an ansatz 

56 that saturates :math:`\mathfrak{su}(2^n)` from enumerating all 

57 :math:`4^n - 1` of them. A result of length ``max_dim`` therefore 

58 means :math:`\dim \mathfrak{g} \geq` ``max_dim`` and the basis may be 

59 partial: a truncated basis no longer spans an algebra, so use it for 

60 the dimension only, not for :func:`g_purity_from_basis`. 

61 

62 Returns: 

63 The list of Pauli words spanning :math:`i \mathfrak{g}`, truncated to 

64 ``max_dim`` entries if the cap was reached. 

65 """ 

66 

67 def _word(s: str) -> PauliWord: 

68 return PauliWord.from_pauli_string(s, list(range(len(s))), len(s)) 

69 

70 seen: set = set() 

71 basis: List[PauliWord] = [] 

72 limit = float("inf") if max_dim is None else max_dim 

73 

74 def _add(s: str) -> bool: 

75 if s in seen: 

76 return False 

77 seen.add(s) 

78 basis.append(_word(s)) 

79 return True 

80 

81 for g in generators: 

82 _add(g if isinstance(g, str) else g.to_pauli_string()) 

83 if len(basis) >= limit: 

84 return basis 

85 

86 frontier = list(basis) 

87 while frontier: 

88 new: List[PauliWord] = [] 

89 pool = list(basis) 

90 for a in frontier: 

91 for b in pool: 

92 if a.commutes_with(b): 

93 continue 

94 if _add(a.compose(b).to_pauli_string()): 

95 new.append(basis[-1]) 

96 if len(basis) >= limit: 

97 return basis 

98 frontier = new 

99 return basis 

100 

101 

102def lie_closure_matrices( 

103 gens: Sequence[np.ndarray], tol: float = 1e-9 

104) -> List[np.ndarray]: 

105 r"""Return an HS-orthonormal Hermitian basis of the real Lie algebra. 

106 

107 The algebra is generated by the Hermitian matrices *gens* under the bracket 

108 :math:`i [A, B]` (which preserves Hermiticity). Each candidate is 

109 orthonormalised against the current basis in the Hilbert-Schmidt inner 

110 product; candidates with norm below *tol* are discarded as dependent. The 

111 length of the result is :math:`\dim \mathfrak{g}`. 

112 

113 Args: 

114 gens: Hermitian generator matrices. 

115 tol: Linear-independence threshold on the Hilbert-Schmidt norm. 

116 

117 Returns: 

118 The HS-orthonormal Hermitian basis matrices of :math:`i \mathfrak{g}`. 

119 """ 

120 basis_vecs: List[np.ndarray] = [] 

121 basis_mats: List[np.ndarray] = [] 

122 

123 def add(M: np.ndarray): 

124 M = np.asarray(M, dtype=complex) 

125 v = M.reshape(-1) 

126 for b in basis_vecs: 

127 v = v - np.vdot(b, v) * b 

128 nrm = np.linalg.norm(v) 

129 if nrm < tol: 

130 return None 

131 v = v / nrm 

132 basis_vecs.append(v) 

133 basis_mats.append(v.reshape(M.shape)) 

134 return basis_mats[-1] 

135 

136 frontier = [m for g in gens if (m := add(g)) is not None] 

137 while frontier: 

138 current = list(basis_mats) 

139 new = [] 

140 for A in frontier: 

141 for B in current: 

142 m = add(1j * (A @ B - B @ A)) 

143 if m is not None: 

144 new.append(m) 

145 frontier = new 

146 return basis_mats 

147 

148 

149def matchgate_generators(n: int) -> List[str]: 

150 r"""Return the matchgate DLA generators :math:`\{Z_k\} \cup \{X_k X_{k+1}\}`. 

151 

152 These Hermitian Pauli strings generate the matchgate algebra 

153 :math:`\mathfrak{so}(2n)` under :func:`lie_closure_paulis`. 

154 

155 Args: 

156 n: Number of qubits. 

157 

158 Returns: 

159 The generator Pauli strings (qubit 0 leftmost). 

160 """ 

161 gens: List[str] = [] 

162 for k in range(n): 

163 s = ["I"] * n 

164 s[k] = "Z" 

165 gens.append("".join(s)) 

166 for k in range(n - 1): 

167 s = ["I"] * n 

168 s[k] = "X" 

169 s[k + 1] = "X" 

170 gens.append("".join(s)) 

171 return gens 

172 

173 

174def matchgate_basis(n: int) -> List[str]: 

175 r"""Return the explicit Pauli-string basis of :math:`\mathfrak{so}(2n)`. 

176 

177 The basis is the :math:`n` on-site :math:`Z_k`, and for each pair 

178 :math:`j < k` the four strings 

179 :math:`\sigma_j \left( \prod_{j < l < k} Z_l \right) \sigma'_k` with 

180 :math:`\sigma, \sigma' \in \{X, Y\}`. Its length is 

181 :math:`n + 4\binom{n}{2} = n(2n-1) = \dim \mathfrak{so}(2n)`. 

182 

183 Args: 

184 n: Number of qubits. 

185 

186 Returns: 

187 The basis Pauli strings (qubit 0 leftmost). 

188 """ 

189 basis: List[str] = [] 

190 for k in range(n): 

191 s = ["I"] * n 

192 s[k] = "Z" 

193 basis.append("".join(s)) 

194 for j, k in combinations(range(n), 2): 

195 for sj in ("X", "Y"): 

196 for sk in ("X", "Y"): 

197 s = ["I"] * n 

198 s[j] = sj 

199 s[k] = sk 

200 for ell in range(j + 1, k): 

201 s[ell] = "Z" 

202 basis.append("".join(s)) 

203 return basis 

204 

205 

206def dim_so2n(n: int) -> int: 

207 r"""Return :math:`\dim \mathfrak{so}(2n) = n(2n-1)`.""" 

208 return n * (2 * n - 1) 

209 

210 

211def g_purity_from_basis( 

212 state: np.ndarray, 

213 basis: Sequence[Union[str, PauliWord]], 

214) -> float: 

215 r"""Return the g-purity :math:`P_g = \sum_B \langle\psi|B|\psi\rangle^2`. 

216 

217 Each basis element :math:`B` is a Hermitian Pauli word, so its expectation 

218 is real and :math:`P_g` is the squared 2-norm of the expectation vector in 

219 the DLA basis. The basis is typically :func:`matchgate_basis` or the output 

220 of :func:`lie_closure_paulis`. 

221 

222 Args: 

223 state: Statevector of shape ``(2**n,)`` (qubit 0 leftmost). 

224 basis: Hermitian Pauli words, each a :class:`PauliWord` or a bare Pauli 

225 string over ``{'I', 'X', 'Y', 'Z'}``. 

226 

227 Returns: 

228 The g-purity :math:`P_g`. 

229 """ 

230 return float(sum(state_expectation(B, state).real ** 2 for B in basis)) 

231 

232 

233def g_purity_matrix( 

234 state: np.ndarray, 

235 basis_mats: Sequence[np.ndarray], 

236) -> float: 

237 r"""Return the g-purity against an HS-orthonormal Hermitian matrix basis. 

238 

239 Computes :math:`P_g^{\mathrm{HS}} = \sum_\mu |\langle\psi|B_\mu|\psi\rangle|^2` 

240 for a Hilbert-Schmidt-orthonormal Hermitian basis :math:`\{B_\mu\}` of 

241 :math:`i\mathfrak{g}`, for example the output of :func:`lie_closure_matrices`. 

242 

243 Args: 

244 state: Statevector of shape ``(2**n,)``. 

245 basis_mats: HS-orthonormal Hermitian basis matrices. 

246 

247 Returns: 

248 The g-purity :math:`P_g^{\mathrm{HS}}`. 

249 """ 

250 psi = np.asarray(state, dtype=complex) 

251 return float(sum(abs(np.vdot(psi, B @ psi)) ** 2 for B in basis_mats)) 

252 

253 

254def symmetric_pauli_sum(pauli: str, n: int, locality: int = 1) -> np.ndarray: 

255 r"""Return the symmetric sum of a Pauli over all ``locality``-qubit subsets. 

256 

257 The permutation-symmetric operator 

258 :math:`\sum_{S} \prod_{q \in S} P_q` where :math:`P` is the single-qubit 

259 Pauli *pauli* and :math:`S` ranges over all size-``locality`` subsets of the 

260 :math:`n` qubits. For example ``locality=1`` gives :math:`\sum_k P_k` and 

261 ``locality=2`` gives :math:`\sum_{j<k} P_j P_k`. The result is Hermitian. 

262 

263 Args: 

264 pauli: Single-qubit Pauli, one of ``{'I', 'X', 'Y', 'Z'}``. 

265 n: Number of qubits. 

266 locality: Subset size, with :math:`1 \leq \text{locality} \leq n`. 

267 

268 Returns: 

269 The dense Hermitian operator of shape :math:`(2^n, 2^n)`. 

270 

271 Raises: 

272 ValueError: If ``pauli`` is not a Pauli label or ``locality`` is outside 

273 ``[1, n]``. 

274 """ 

275 if pauli not in {"I", "X", "Y", "Z"}: 

276 raise ValueError(f"pauli must be one of I, X, Y, Z, got {pauli!r}") 

277 if not 1 <= locality <= n: 

278 raise ValueError( 

279 f"locality must satisfy 1 <= locality <= n, " 

280 f"got locality={locality} for n={n}" 

281 ) 

282 dim = 2**n 

283 acc = np.zeros((dim, dim), dtype=complex) 

284 for subset in combinations(range(n), locality): 

285 s = ["I"] * n 

286 for q in subset: 

287 s[q] = pauli 

288 acc += np.asarray( 

289 PauliWord.from_pauli_string("".join(s), list(range(n)), n).to_matrix() 

290 ) 

291 return acc 

292 

293 

294def sn_equivariant_generators(n: int) -> List[np.ndarray]: 

295 r"""Return the :math:`S_n`-equivariant generators 

296 :math:`\{\sum_k X_k, \sum_k Y_k, \sum_{j<k} Z_j Z_k\}`. 

297 

298 The three Hermitian generators of the permutation-equivariant ansatz of 

299 Schatzki et al. (arXiv:2210.09974). They are unnormalised (the scaling is 

300 irrelevant for the Lie closure) and feed :func:`lie_closure_matrices`. 

301 

302 Args: 

303 n: Number of qubits, with :math:`n \geq 2` (the pair term needs two). 

304 

305 Returns: 

306 The list of three generator matrices, each of shape :math:`(2^n, 2^n)`. 

307 

308 Raises: 

309 ValueError: If ``n < 2``. 

310 """ 

311 if n < 2: 

312 raise ValueError(f"n must be at least 2, got {n}") 

313 return [ 

314 symmetric_pauli_sum("X", n, 1), 

315 symmetric_pauli_sum("Y", n, 1), 

316 symmetric_pauli_sum("Z", n, 2), 

317 ] 

318 

319 

320def sn_equivariant_observable(n: int) -> np.ndarray: 

321 r"""Return the :math:`S_n`-equivariant observable 

322 :math:`O = \frac{2}{n(n-1)} \sum_{j<k} X_j X_k`. 

323 

324 The permutation-symmetric observable of Schatzki et al. (arXiv:2210.09974), 

325 normalised by the number of pairs. 

326 

327 Args: 

328 n: Number of qubits, with :math:`n \geq 2`. 

329 

330 Returns: 

331 The dense Hermitian observable of shape :math:`(2^n, 2^n)`. 

332 

333 Raises: 

334 ValueError: If ``n < 2``. 

335 """ 

336 if n < 2: 

337 raise ValueError(f"n must be at least 2, got {n}") 

338 return (2.0 / (n * (n - 1))) * symmetric_pauli_sum("X", n, 2)