1# Licensed under the Apache License, Version 2.0 (the "License"); 2# you may not use this file except in compliance with the License. 3# You may obtain a copy of the License at 4# 5# http://www.apache.org/licenses/LICENSE-2.0 6# 7# Unless required by applicable law or agreed to in writing, software 8# distributed under the License is distributed on an "AS IS" BASIS, 9# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 10# See the License for the specific language governing permissions and 11# limitations under the License. 12"""bravyi_kitaev_fast transform on fermionic operators.""" 13 14from typing import Optional 15 16import numpy 17import networkx 18 19from openfermion.ops.operators import QubitOperator 20from openfermion.ops.representations import InteractionOperator 21from openfermion.utils.operator_utils import count_qubits 22 23 24def bravyi_kitaev_fast(operator: InteractionOperator) -> QubitOperator: 25 """ 26 Find the Pauli-representation of InteractionOperator for Bravyi-Kitaev 27 Super fast (BKSF) algorithm. Pauli-representation of general 28 FermionOperator is not possible in BKSF. Also, the InteractionOperator 29 given as input must be Hermitian. In future we might provide a 30 transformation for a restricted set of fermion operator. 31 32 Args: 33 operator: InteractionOperator. 34 35 Returns: 36 transformed_operator: An instance of the QubitOperator class. 37 38 Raises: 39 TypeError: If operator is not an InteractionOperator 40 41 """ 42 if isinstance(operator, InteractionOperator): 43 return bravyi_kitaev_fast_interaction_op(operator) 44 else: 45 raise TypeError("operator must be an InteractionOperator.") 46 47 48def bravyi_kitaev_fast_interaction_op(iop: InteractionOperator 49 ) -> QubitOperator: 50 r""" 51 Transform from InteractionOperator to QubitOperator for Bravyi-Kitaev fast 52 algorithm. 53 54 The electronic Hamiltonian is represented in terms of creation and 55 annihilation operators. These creation and annihilation operators could be 56 used to define Majorana modes as follows: 57 c_{2i} = a_i + a^{\dagger}_i, 58 c_{2i+1} = (a_i - a^{\dagger}_{i})/(1j) 59 These Majorana modes can be used to define edge operators B_i and A_{ij}: 60 B_i=c_{2i}c_{2i+1}, 61 A_{ij}=c_{2i}c_{2j} 62 using these edge operators the fermionic algebra can be generated and 63 hence all the terms in the electronic Hamiltonian can be expressed in 64 terms of edge operators. The terms in electronic Hamiltonian can be 65 divided into five types (arXiv 1208.5986). We can find the edge operator 66 expression for each of those five types. For example, the excitation 67 operator term in Hamiltonian when represented in terms of edge operators 68 becomes: 69 a_i^{\dagger}a_j+a_j^{\dagger}a_i = (-1j/2)*(A_ij*B_i+B_j*A_ij) 70 For the sake of brevity the reader is encouraged to look up the 71 expressions of other terms from the code below. The variables for edge 72 operators are chosen according to the nomenclature defined above 73 (B_i and A_ij). A detailed description of these operators and the terms 74 of the electronic Hamiltonian are provided in (arXiv 1712.00446). 75 76 Args: 77 iop (InteractionOperator): 78 n_qubit (int): Number of qubits 79 80 Returns: 81 qubit_operator: An instance of the QubitOperator class. 82 """ 83 n_qubits = count_qubits(iop) 84 85 # Initialize qubit operator as constant. 86 qubit_operator = QubitOperator((), iop.constant) 87 edge_matrix = bravyi_kitaev_fast_edge_matrix(iop) 88 edge_matrix_indices = numpy.array( 89 numpy.nonzero( 90 numpy.triu(edge_matrix) - numpy.diag(numpy.diag(edge_matrix)))) 91 # Loop through all indices. 92 for p in range(n_qubits): 93 for q in range(n_qubits): 94 # Handle one-body terms. 95 coefficient = complex(iop[(p, 1), (q, 0)]) 96 if coefficient and p >= q: 97 qubit_operator += (coefficient * 98 _one_body(edge_matrix_indices, p, q)) 99 100 # Keep looping for the two-body terms. 101 for r in range(n_qubits): 102 for s in range(n_qubits): 103 coefficient = complex(iop[(p, 1), (q, 1), (r, 0), (s, 0)]) 104 105 # Skip zero terms. 106 if (not coefficient) or (p == q) or (r == s): 107 # coverage: ignore 108 continue 109 110 # Identify and skip one of the complex conjugates. 111 if [p, q, r, s] != [s, r, q, p]: 112 if len(set([p, q, r, s])) == 4: 113 if min(r, s) < min(p, q): 114 continue 115 # Handle case of 3 unique indices 116 elif len(set([p, q, r, s])) == 3: 117 transformed_term = _two_body( 118 edge_matrix_indices, p, q, r, s) 119 transformed_term *= .5 * coefficient 120 qubit_operator += transformed_term 121 continue 122 elif p != r and q < p: 123 # TODO: remove pragma if reachable continue 124 continue # pragma: no cover 125 126 # Handle the two-body terms. 127 transformed_term = _two_body(edge_matrix_indices, p, q, r, 128 s) 129 transformed_term *= coefficient 130 qubit_operator += transformed_term 131 return qubit_operator 132 133 134def bravyi_kitaev_fast_edge_matrix(iop: InteractionOperator, 135 n_qubits: Optional[int] = None 136 ) -> numpy.ndarray: 137 """ 138 Use InteractionOperator to construct edge matrix required for the algorithm 139 140 Edge matrix contains the information about the edges between vertices. 141 Edge matrix is required to build the operators in bravyi_kitaev_fast model. 142 143 Args: 144 iop (InteractionOperator): 145 n_qubits (int): number of qubits 146 147 Returns: 148 edge_matrix (Numpy array): A square numpy array containing information 149 about the edges present in the model. 150 """ 151 n_qubits = count_qubits(iop) 152 edge_matrix = 1j * numpy.zeros((n_qubits, n_qubits)) 153 # Loop through all indices. 154 for p in range(n_qubits): 155 for q in range(n_qubits): 156 157 # Handle one-body terms. 158 coefficient = complex(iop[(p, 1), (q, 0)]) 159 if coefficient and p >= q: 160 edge_matrix[p, q] = bool(complex(iop[(p, 1), (q, 0)])) 161 162 # Keep looping for the two-body terms. 163 for r in range(n_qubits): 164 for s in range(n_qubits): 165 coefficient2 = complex(iop[(p, 1), (q, 1), (r, 0), (s, 0)]) 166 167 # Skip zero terms. 168 if (not coefficient2) or (p == q) or (r == s): 169 # TODO: remove pragma if this is a reachable continue 170 continue # pragma: no cover 171 172 # Identify and skip one of the complex conjugates. 173 if [p, q, r, s] != [s, r, q, p]: 174 if len(set([p, q, r, s])) == 4: 175 if min(r, s) < min(p, q): 176 continue 177 elif p != r and q < p: 178 continue 179 180 # Handle case of four unique indices. 181 if len(set([p, q, r, s])) == 4: 182 183 if coefficient2 and p >= q: 184 edge_matrix[p, q] = bool( 185 complex(iop[(p, 1), (q, 1), (r, 0), (s, 0)])) 186 a, b = sorted([r, s]) 187 edge_matrix[b, a] = bool( 188 complex(iop[(p, 1), (q, 1), (r, 0), (s, 0)])) 189 190 # Handle case of three unique indices. 191 elif len(set([p, q, r, s])) == 3: 192 193 # Identify equal tensor factors. 194 if p == r: 195 a, b = sorted([q, s]) 196 edge_matrix[b, a] = bool( 197 complex(iop[(p, 1), (q, 1), (r, 0), (s, 0)])) 198 elif p == s: 199 a, b = sorted([q, r]) 200 edge_matrix[b, a] = bool( 201 complex(iop[(p, 1), (q, 1), (r, 0), (s, 0)])) 202 elif q == r: 203 a, b = sorted([p, s]) 204 edge_matrix[b, a] = bool( 205 complex(iop[(p, 1), (q, 1), (r, 0), (s, 0)])) 206 elif q == s: 207 a, b = sorted([p, r]) 208 edge_matrix[b, a] = bool( 209 complex(iop[(p, 1), (q, 1), (r, 0), (s, 0)])) 210 211 return edge_matrix.transpose() 212 213 214def _one_body(edge_matrix_indices: numpy.ndarray, p: int, 215 q: int) -> QubitOperator: 216 r""" 217 Map the term a^\dagger_p a_q + a^\dagger_q a_p to QubitOperator. 218 219 The definitions for various operators will be presented in a paper soon. 220 221 Args: 222 edge_matrix_indices(numpy array): Specifying the edges 223 p and q (int): specifying the one body term. 224 225 Return: 226 An instance of QubitOperator() 227 """ 228 # Handle off-diagonal terms. 229 qubit_operator = QubitOperator() 230 if p != q: 231 a, b = sorted([p, q]) 232 B_a = edge_operator_b(edge_matrix_indices, a) 233 B_b = edge_operator_b(edge_matrix_indices, b) 234 A_ab = edge_operator_aij(edge_matrix_indices, a, b) 235 qubit_operator += -1j / 2. * (A_ab * B_b + B_a * A_ab) 236 237 # Handle diagonal terms. 238 else: 239 B_p = edge_operator_b(edge_matrix_indices, p) 240 qubit_operator += (QubitOperator((), 1) - B_p) / 2. 241 242 return qubit_operator 243 244 245def _two_body(edge_matrix_indices: numpy.ndarray, p: int, q: int, r: int, 246 s: int) -> QubitOperator: 247 r""" 248 Map the term a^\dagger_p a^\dagger_q a_r a_s + h.c. to QubitOperator. 249 250 The definitions for various operators will be covered in a paper soon. 251 252 Args: 253 edge_matrix_indices (numpy array): specifying the edges 254 p, q, r and s (int): specifying the two body term. 255 256 Return: 257 An instance of QubitOperator() 258 """ 259 # Initialize qubit operator. 260 qubit_operator = QubitOperator() 261 262 # Handle case of four unique indices. 263 if len(set([p, q, r, s])) == 4: 264 B_p = edge_operator_b(edge_matrix_indices, p) 265 B_q = edge_operator_b(edge_matrix_indices, q) 266 B_r = edge_operator_b(edge_matrix_indices, r) 267 B_s = edge_operator_b(edge_matrix_indices, s) 268 A_pq = edge_operator_aij(edge_matrix_indices, p, q) 269 A_rs = edge_operator_aij(edge_matrix_indices, r, s) 270 qubit_operator += 1 / 8. * A_pq * A_rs * (-QubitOperator( 271 ()) - B_p * B_q + B_p * B_r + B_p * B_s + B_q * B_r + B_q * B_s - 272 B_r * B_s + 273 B_p * B_q * B_r * B_s) 274 return qubit_operator 275 276 # Handle case of three unique indices. 277 elif len(set([p, q, r, s])) == 3: 278 # Identify equal tensor factors. 279 if p == r: 280 B_p = edge_operator_b(edge_matrix_indices, p) 281 B_q = edge_operator_b(edge_matrix_indices, q) 282 B_s = edge_operator_b(edge_matrix_indices, s) 283 A_qs = edge_operator_aij(edge_matrix_indices, q, s) 284 qubit_operator += 1j * (A_qs * B_s + B_q * A_qs) * (QubitOperator( 285 ()) - B_p) / 4. 286 287 elif p == s: 288 B_p = edge_operator_b(edge_matrix_indices, p) 289 B_q = edge_operator_b(edge_matrix_indices, q) 290 B_r = edge_operator_b(edge_matrix_indices, r) 291 A_qr = edge_operator_aij(edge_matrix_indices, q, r) 292 qubit_operator += -1j * (A_qr * B_r + B_q * A_qr) * (QubitOperator( 293 ()) - B_p) / 4. 294 295 elif q == r: 296 B_p = edge_operator_b(edge_matrix_indices, p) 297 B_q = edge_operator_b(edge_matrix_indices, q) 298 B_s = edge_operator_b(edge_matrix_indices, s) 299 A_ps = edge_operator_aij(edge_matrix_indices, p, s) 300 qubit_operator += -1j * (A_ps * B_s + B_p * A_ps) * (QubitOperator( 301 ()) - B_q) / 4. 302 303 elif q == s: 304 B_p = edge_operator_b(edge_matrix_indices, p) 305 B_q = edge_operator_b(edge_matrix_indices, q) 306 B_r = edge_operator_b(edge_matrix_indices, r) 307 A_pr = edge_operator_aij(edge_matrix_indices, p, r) 308 qubit_operator += 1j * (A_pr * B_r + B_p * A_pr) * (QubitOperator( 309 ()) - B_q) / 4. 310 311 # Handle case of two unique indices. 312 elif len(set([p, q, r, s])) == 2: 313 # Get coefficient. 314 if p == s: 315 B_p = edge_operator_b(edge_matrix_indices, p) 316 B_q = edge_operator_b(edge_matrix_indices, q) 317 qubit_operator += (QubitOperator(()) - B_p) * (QubitOperator( 318 ()) - B_q) / 4. 319 320 else: 321 B_p = edge_operator_b(edge_matrix_indices, p) 322 B_q = edge_operator_b(edge_matrix_indices, q) 323 qubit_operator += -(QubitOperator(()) - B_p) * (QubitOperator( 324 ()) - B_q) / 4. 325 326 return qubit_operator 327 328 329def edge_operator_b(edge_matrix_indices: numpy.ndarray, 330 i: int) -> QubitOperator: 331 """Calculate the edge operator B_i. The definitions used here are 332 consistent with arXiv:quant-ph/0003137 333 334 Args: 335 edge_matrix_indices(numpy array(square and symmetric): 336 i (int): index for specifying the edge operator B. 337 338 Returns: 339 An instance of QubitOperator 340 """ 341 B_i = QubitOperator() 342 qubit_position_matrix = numpy.array(numpy.where(edge_matrix_indices == i)) 343 qubit_position = qubit_position_matrix[1][:] 344 qubit_position = numpy.sort(qubit_position) 345 operator = tuple() 346 for d1 in qubit_position: 347 operator += ((int(d1), 'Z'),) 348 B_i += QubitOperator(operator) 349 return B_i 350 351 352def edge_operator_aij(edge_matrix_indices: numpy.ndarray, i: int, 353 j: int) -> QubitOperator: 354 """Calculate the edge operator A_ij. The definitions used here are 355 consistent with arXiv:quant-ph/0003137 356 357 Args: 358 edge_matrix_indices(numpy array): specifying the edges 359 i,j (int): specifying the edge operator A 360 361 Returns: 362 An instance of QubitOperator 363 """ 364 a_ij = QubitOperator() 365 operator = tuple() 366 position_ij = -1 367 qubit_position_i = numpy.array(numpy.where(edge_matrix_indices == i)) 368 for edge_index in range(numpy.size(edge_matrix_indices[0, :])): 369 if set((i, j)) == set(edge_matrix_indices[:, edge_index]): 370 position_ij = edge_index 371 operator += ((int(position_ij), 'X'),) 372 373 for edge_index in range(numpy.size(qubit_position_i[0, :])): 374 if edge_matrix_indices[int(not (qubit_position_i[0, edge_index]))][ 375 qubit_position_i[1, edge_index]] < j: 376 operator += ((int(qubit_position_i[1, edge_index]), 'Z'),) 377 qubit_position_j = numpy.array(numpy.where(edge_matrix_indices == j)) 378 for edge_index in range(numpy.size(qubit_position_j[0, :])): 379 if edge_matrix_indices[int(not (qubit_position_j[0, edge_index]))][ 380 qubit_position_j[1, edge_index]] < i: 381 operator += ((int(qubit_position_j[1, edge_index]), 'Z'),) 382 a_ij += QubitOperator(operator, 1) 383 if j < i: 384 a_ij = -1 * a_ij 385 return a_ij 386 387 388def vacuum_operator(edge_matrix_indices: numpy.ndarray) -> QubitOperator: 389 """Use the stabilizers to find the vacuum state in bravyi_kitaev_fast. 390 391 Args: 392 edge_matrix_indices(numpy array): specifying the edges 393 394 Return: 395 An instance of QubitOperator 396 397 """ 398 # Initialize qubit operator. 399 g = networkx.Graph() 400 g.add_edges_from(tuple(edge_matrix_indices.transpose())) 401 stabs = numpy.array(networkx.cycle_basis(g)) 402 vac_operator = QubitOperator(()) 403 for stab in stabs: 404 405 A = QubitOperator(()) 406 stab = numpy.array(stab) 407 for i in range(numpy.size(stab)): 408 if i == (numpy.size(stab) - 1): 409 A = (1j) * A * edge_operator_aij(edge_matrix_indices, stab[i], 410 stab[0]) 411 else: 412 A = (1j) * A * edge_operator_aij(edge_matrix_indices, stab[i], 413 stab[i + 1]) 414 vac_operator = vac_operator * (QubitOperator(()) + A) / numpy.sqrt(2) 415 416 return vac_operator 417 418 419def number_operator(iop, mode_number=None): 420 """Find the qubit operator for the number operator in bravyi_kitaev_fast 421 representation 422 423 Args: 424 iop (InteractionOperator): 425 mode_number: index mode_number corresponding to the mode 426 for which number operator is required. 427 428 Return: 429 A QubitOperator 430 431 """ 432 n_qubit = iop.n_qubits 433 num_operator = QubitOperator() 434 edge_matrix = bravyi_kitaev_fast_edge_matrix(iop) 435 edge_matrix_indices = numpy.array( 436 numpy.nonzero( 437 numpy.triu(edge_matrix) - numpy.diag(numpy.diag(edge_matrix)))) 438 if mode_number is None: 439 for i in range(n_qubit): 440 num_operator += (QubitOperator( 441 ()) - edge_operator_b(edge_matrix_indices, i)) / 2. 442 443 else: 444 num_operator += (QubitOperator( 445 ()) - edge_operator_b(edge_matrix_indices, mode_number)) / 2. 446 447 return num_operator 448 449 450def generate_fermions(edge_matrix_indices: numpy.ndarray, i: int, 451 j: int) -> QubitOperator: 452 """The QubitOperator for generating fermions in bravyi_kitaev_fast 453 representation 454 455 Args: 456 edge_matrix_indices(numpy array): specifying the edges 457 i and j (int) 458 459 Return: 460 A QubitOperator 461 """ 462 gen_fer_operator = (-1j / 463 2.) * (edge_operator_aij(edge_matrix_indices, i, j) * 464 edge_operator_b(edge_matrix_indices, j) - 465 edge_operator_b(edge_matrix_indices, i) * 466 edge_operator_aij(edge_matrix_indices, i, j)) 467 return gen_fer_operator 468