English
Languages
English
Bengali
French
German
Japanese
Korean
Portuguese
Spanish
Tamil

Source code for qiskit.quantum_info.synthesis.clifford_decompose

# This code is part of Qiskit.
#
# (C) Copyright IBM 2017, 2021
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.
"""
Circuit synthesis for the Clifford class.
"""
# pylint: disable=invalid-name

from itertools import product

import numpy as np

from qiskit.circuit import QuantumCircuit
from qiskit.exceptions import QiskitError
from qiskit.quantum_info.operators.symplectic.clifford_circuits import (
    _append_cx,
    _append_h,
    _append_s,
    _append_swap,
    _append_v,
    _append_w,
    _append_x,
    _append_z,
)
from qiskit.quantum_info.operators.symplectic.pauli import Pauli


[docs]def decompose_clifford(clifford, method=None): """Decompose a Clifford operator into a QuantumCircuit. For N <= 3 qubits this is based on optimal CX cost decomposition from reference [1]. For N > 3 qubits this is done using the general non-optimal greedy compilation routine from reference [3], which typically yields better CX cost compared to the AG method in [2]. Args: clifford (Clifford): a clifford operator. method (str): Optional, a synthesis method ('AG' or 'greedy'). If set this overrides optimal decomposition for N <=3 qubits. Return: QuantumCircuit: a circuit implementation of the Clifford. References: 1. S. Bravyi, D. Maslov, *Hadamard-free circuits expose the structure of the Clifford group*, `arXiv:2003.09412 [quant-ph] <https://arxiv.org/abs/2003.09412>`_ 2. S. Aaronson, D. Gottesman, *Improved Simulation of Stabilizer Circuits*, Phys. Rev. A 70, 052328 (2004). `arXiv:quant-ph/0406196 <https://arxiv.org/abs/quant-ph/0406196>`_ 3. Sergey Bravyi, Shaohan Hu, Dmitri Maslov, Ruslan Shaydulin, *Clifford Circuit Optimization with Templates and Symbolic Pauli Gates*, `arXiv:2105.02291 [quant-ph] <https://arxiv.org/abs/2105.02291>`_ """ num_qubits = clifford.num_qubits if method == "AG": return decompose_clifford_ag(clifford) if method == "greedy": return decompose_clifford_greedy(clifford) if num_qubits <= 3: return decompose_clifford_bm(clifford) return decompose_clifford_greedy(clifford)
# --------------------------------------------------------------------- # Synthesis based on Bravyi & Maslov decomposition # --------------------------------------------------------------------- def decompose_clifford_bm(clifford): """Decompose a clifford""" num_qubits = clifford.num_qubits if num_qubits == 1: return _decompose_clifford_1q(clifford.tableau) clifford_name = str(clifford) # Inverse of final decomposed circuit inv_circuit = QuantumCircuit(num_qubits, name="inv_circ") # CNOT cost of clifford cost = _cx_cost(clifford) # Find composition of circuits with CX and (H.S)^a gates to reduce CNOT count while cost > 0: clifford, inv_circuit, cost = _reduce_cost(clifford, inv_circuit, cost) # Decompose the remaining product of 1-qubit cliffords ret_circ = QuantumCircuit(num_qubits, name=clifford_name) for qubit in range(num_qubits): pos = [qubit, qubit + num_qubits] circ = _decompose_clifford_1q(clifford.tableau[pos][:, pos + [-1]]) if len(circ) > 0: ret_circ.append(circ, [qubit]) # Add the inverse of the 2-qubit reductions circuit if len(inv_circuit) > 0: ret_circ.append(inv_circuit.inverse(), range(num_qubits)) return ret_circ.decompose() # --------------------------------------------------------------------- # Synthesis based on Aaronson & Gottesman decomposition # --------------------------------------------------------------------- def decompose_clifford_ag(clifford): """Decompose a Clifford operator into a QuantumCircuit. Args: clifford (Clifford): a clifford operator. Return: QuantumCircuit: a circuit implementation of the Clifford. """ # Use 1-qubit decomposition method if clifford.num_qubits == 1: return _decompose_clifford_1q(clifford.tableau) # Compose a circuit which we will convert to an instruction circuit = QuantumCircuit(clifford.num_qubits, name=str(clifford)) # Make a copy of Clifford as we are going to do row reduction to # reduce it to an identity clifford_cpy = clifford.copy() for i in range(clifford.num_qubits): # put a 1 one into position by permuting and using Hadamards(i,i) _set_qubit_x_true(clifford_cpy, circuit, i) # make all entries in row i except ith equal to 0 # by using phase gate and CNOTS _set_row_x_zero(clifford_cpy, circuit, i) # treat Zs _set_row_z_zero(clifford_cpy, circuit, i) for i in range(clifford.num_qubits): if clifford_cpy.destab_phase[i]: _append_z(clifford_cpy, i) circuit.z(i) if clifford_cpy.stab_phase[i]: _append_x(clifford_cpy, i) circuit.x(i) # Next we invert the circuit to undo the row reduction and return the # result as a gate instruction return circuit.inverse() # --------------------------------------------------------------------- # 1-qubit Clifford decomposition # --------------------------------------------------------------------- def _decompose_clifford_1q(tableau): """Decompose a single-qubit clifford""" circuit = QuantumCircuit(1, name="temp") # Add phase correction destab_phase, stab_phase = tableau[:, 2] if destab_phase and not stab_phase: circuit.z(0) elif not destab_phase and stab_phase: circuit.x(0) elif destab_phase and stab_phase: circuit.y(0) destab_phase_label = "-" if destab_phase else "+" stab_phase_label = "-" if stab_phase else "+" destab_x, destab_z = tableau[0, 0], tableau[0, 1] stab_x, stab_z = tableau[1, 0], tableau[1, 1] # Z-stabilizer if stab_z and not stab_x: stab_label = "Z" if destab_z: destab_label = "Y" circuit.s(0) else: destab_label = "X" # X-stabilizer elif not stab_z and stab_x: stab_label = "X" if destab_x: destab_label = "Y" circuit.sdg(0) else: destab_label = "Z" circuit.h(0) # Y-stabilizer else: stab_label = "Y" if destab_z: destab_label = "Z" else: destab_label = "X" circuit.s(0) circuit.h(0) circuit.s(0) # Add circuit name name_destab = f"Destabilizer = ['{destab_phase_label}{destab_label}']" name_stab = f"Stabilizer = ['{stab_phase_label}{stab_label}']" circuit.name = f"Clifford: {name_stab}, {name_destab}" return circuit # --------------------------------------------------------------------- # Helper functions for Bravyi & Maslov decomposition # --------------------------------------------------------------------- def _reduce_cost(clifford, inv_circuit, cost): """Two-qubit cost reduction step""" num_qubits = clifford.num_qubits for qubit0 in range(num_qubits): for qubit1 in range(qubit0 + 1, num_qubits): for n0, n1 in product(range(3), repeat=2): # Apply a 2-qubit block reduced = clifford.copy() for qubit, n in [(qubit0, n0), (qubit1, n1)]: if n == 1: _append_v(reduced, qubit) elif n == 2: _append_w(reduced, qubit) _append_cx(reduced, qubit0, qubit1) # Compute new cost new_cost = _cx_cost(reduced) if new_cost == cost - 1: # Add decomposition to inverse circuit for qubit, n in [(qubit0, n0), (qubit1, n1)]: if n == 1: inv_circuit.sdg(qubit) inv_circuit.h(qubit) elif n == 2: inv_circuit.h(qubit) inv_circuit.s(qubit) inv_circuit.cx(qubit0, qubit1) return reduced, inv_circuit, new_cost # If we didn't reduce cost raise QiskitError("Failed to reduce Clifford CX cost.") def _cx_cost(clifford): """Return the number of CX gates required for Clifford decomposition.""" if clifford.num_qubits == 2: return _cx_cost2(clifford) if clifford.num_qubits == 3: return _cx_cost3(clifford) raise Exception("No Clifford CX cost function for num_qubits > 3.") def _rank2(a, b, c, d): """Return rank of 2x2 boolean matrix.""" if (a & d) ^ (b & c): return 2 if a or b or c or d: return 1 return 0 def _cx_cost2(clifford): """Return CX cost of a 2-qubit clifford.""" U = clifford.tableau[:, :-1] r00 = _rank2(U[0, 0], U[0, 2], U[2, 0], U[2, 2]) r01 = _rank2(U[0, 1], U[0, 3], U[2, 1], U[2, 3]) if r00 == 2: return r01 return r01 + 1 - r00 def _cx_cost3(clifford): """Return CX cost of a 3-qubit clifford.""" # pylint: disable=too-many-return-statements,too-many-boolean-expressions U = clifford.tableau[:, :-1] n = 3 # create information transfer matrices R1, R2 R1 = np.zeros((n, n), dtype=int) R2 = np.zeros((n, n), dtype=int) for q1 in range(n): for q2 in range(n): R2[q1, q2] = _rank2(U[q1, q2], U[q1, q2 + n], U[q1 + n, q2], U[q1 + n, q2 + n]) mask = np.zeros(2 * n, dtype=int) mask[[q2, q2 + n]] = 1 isLocX = np.array_equal(U[q1, :] & mask, U[q1, :]) isLocZ = np.array_equal(U[q1 + n, :] & mask, U[q1 + n, :]) isLocY = np.array_equal((U[q1, :] ^ U[q1 + n, :]) & mask, (U[q1, :] ^ U[q1 + n, :])) R1[q1, q2] = 1 * (isLocX or isLocZ or isLocY) + 1 * (isLocX and isLocZ and isLocY) diag1 = np.sort(np.diag(R1)).tolist() diag2 = np.sort(np.diag(R2)).tolist() nz1 = np.count_nonzero(R1) nz2 = np.count_nonzero(R2) if diag1 == [2, 2, 2]: return 0 if diag1 == [1, 1, 2]: return 1 if ( diag1 == [0, 1, 1] or (diag1 == [1, 1, 1] and nz2 < 9) or (diag1 == [0, 0, 2] and diag2 == [1, 1, 2]) ): return 2 if ( (diag1 == [1, 1, 1] and nz2 == 9) or ( diag1 == [0, 0, 1] and (nz1 == 1 or diag2 == [2, 2, 2] or (diag2 == [1, 1, 2] and nz2 < 9)) ) or (diag1 == [0, 0, 2] and diag2 == [0, 0, 2]) or (diag2 == [1, 2, 2] and nz1 == 0) ): return 3 if diag2 == [0, 0, 1] or ( diag1 == [0, 0, 0] and ( (diag2 == [1, 1, 1] and nz2 == 9 and nz1 == 3) or (diag2 == [0, 1, 1] and nz2 == 8 and nz1 == 2) ) ): return 5 if nz1 == 3 and nz2 == 3: return 6 return 4 # --------------------------------------------------------------------- # Helper functions for Aaronson & Gottesman decomposition # --------------------------------------------------------------------- def _set_qubit_x_true(clifford, circuit, qubit): """Set destabilizer.X[qubit, qubit] to be True. This is done by permuting columns l > qubit or if necessary applying a Hadamard """ x = clifford.destab_x[qubit] z = clifford.destab_z[qubit] if x[qubit]: return # Try to find non-zero element for i in range(qubit + 1, clifford.num_qubits): if x[i]: _append_swap(clifford, i, qubit) circuit.swap(i, qubit) return # no non-zero element found: need to apply Hadamard somewhere for i in range(qubit, clifford.num_qubits): if z[i]: _append_h(clifford, i) circuit.h(i) if i != qubit: _append_swap(clifford, i, qubit) circuit.swap(i, qubit) return def _set_row_x_zero(clifford, circuit, qubit): """Set destabilizer.X[qubit, i] to False for all i > qubit. This is done by applying CNOTS assumes k<=N and A[k][k]=1 """ x = clifford.destab_x[qubit] z = clifford.destab_z[qubit] # Check X first for i in range(qubit + 1, clifford.num_qubits): if x[i]: _append_cx(clifford, qubit, i) circuit.cx(qubit, i) # Check whether Zs need to be set to zero: if np.any(z[qubit:]): if not z[qubit]: # to treat Zs: make sure row.Z[k] to True _append_s(clifford, qubit) circuit.s(qubit) # reverse CNOTS for i in range(qubit + 1, clifford.num_qubits): if z[i]: _append_cx(clifford, i, qubit) circuit.cx(i, qubit) # set row.Z[qubit] to False _append_s(clifford, qubit) circuit.s(qubit) def _set_row_z_zero(clifford, circuit, qubit): """Set stabilizer.Z[qubit, i] to False for all i > qubit. Implemented by applying (reverse) CNOTS assumes qubit < num_qubits and _set_row_x_zero has been called first """ x = clifford.stab_x[qubit] z = clifford.stab_z[qubit] # check whether Zs need to be set to zero: if np.any(z[qubit + 1 :]): for i in range(qubit + 1, clifford.num_qubits): if z[i]: _append_cx(clifford, i, qubit) circuit.cx(i, qubit) # check whether Xs need to be set to zero: if np.any(x[qubit:]): _append_h(clifford, qubit) circuit.h(qubit) for i in range(qubit + 1, clifford.num_qubits): if x[i]: _append_cx(clifford, qubit, i) circuit.cx(qubit, i) if z[qubit]: _append_s(clifford, qubit) circuit.s(qubit) _append_h(clifford, qubit) circuit.h(qubit) # --------------------------------------------------------------------- # Synthesis based on Bravyi et. al. greedy clifford compiler # --------------------------------------------------------------------- def decompose_clifford_greedy(clifford): """Decompose a Clifford operator into a QuantumCircuit. Args: clifford (Clifford): a clifford operator. Return: QuantumCircuit: a circuit implementation of the Clifford. Raises: QiskitError: if symplectic Gaussian elimination fails. """ num_qubits = clifford.num_qubits circ = QuantumCircuit(num_qubits, name=str(clifford)) qubit_list = list(range(num_qubits)) clifford_cpy = clifford.copy() # Reducing the original Clifford to identity # via symplectic Gaussian elimination while len(qubit_list) > 0: # Calculate the adjoint of clifford_cpy without the phase clifford_adj = clifford_cpy.copy() tmp = clifford_adj.destab_x.copy() clifford_adj.destab_x = clifford_adj.stab_z.T clifford_adj.destab_z = clifford_adj.destab_z.T clifford_adj.stab_x = clifford_adj.stab_x.T clifford_adj.stab_z = tmp.T list_greedy_cost = [] for qubit in qubit_list: pauli_x = Pauli("I" * (num_qubits - qubit - 1) + "X" + "I" * qubit) pauli_x = pauli_x.evolve(clifford_adj, frame="s") pauli_z = Pauli("I" * (num_qubits - qubit - 1) + "Z" + "I" * qubit) pauli_z = pauli_z.evolve(clifford_adj, frame="s") list_pairs = [] pauli_count = 0 # Compute the CNOT cost in order to find the qubit with the minimal cost for i in qubit_list: typeq = _from_pair_paulis_to_type(pauli_x, pauli_z, i) list_pairs.append(typeq) pauli_count += 1 cost = _compute_greedy_cost(list_pairs) list_greedy_cost.append([cost, qubit]) _, min_qubit = (sorted(list_greedy_cost))[0] # Gaussian elimination step for the qubit with minimal CNOT cost pauli_x = Pauli("I" * (num_qubits - min_qubit - 1) + "X" + "I" * min_qubit) pauli_x = pauli_x.evolve(clifford_adj, frame="s") pauli_z = Pauli("I" * (num_qubits - min_qubit - 1) + "Z" + "I" * min_qubit) pauli_z = pauli_z.evolve(clifford_adj, frame="s") # Compute the decoupling operator of cliff_ox and cliff_oz decouple_circ, decouple_cliff = _calc_decoupling( pauli_x, pauli_z, qubit_list, min_qubit, num_qubits, clifford_cpy ) circ = circ.compose(decouple_circ) # Now the clifford acts trivially on min_qubit clifford_cpy = decouple_cliff.adjoint().compose(clifford_cpy) qubit_list.remove(min_qubit) # Add the phases (Pauli gates) to the Clifford circuit for qubit in range(num_qubits): stab = clifford_cpy.stab_phase[qubit] destab = clifford_cpy.destab_phase[qubit] if destab and stab: circ.y(qubit) elif not destab and stab: circ.x(qubit) elif destab and not stab: circ.z(qubit) return circ # --------------------------------------------------------------------- # Helper functions for Bravyi et. al. greedy clifford compiler # --------------------------------------------------------------------- # Global arrays of the 16 pairs of Pauli operators # divided into 5 equivalence classes under the action of single-qubit Cliffords # Class A - canonical representative is 'XZ' A_class = [ [[False, True], [True, True]], # 'XY' [[False, True], [True, False]], # 'XZ' [[True, True], [False, True]], # 'YX' [[True, True], [True, False]], # 'YZ' [[True, False], [False, True]], # 'ZX' [[True, False], [True, True]], ] # 'ZY' # Class B - canonical representative is 'XX' B_class = [ [[True, False], [True, False]], # 'ZZ' [[False, True], [False, True]], # 'XX' [[True, True], [True, True]], ] # 'YY' # Class C - canonical representative is 'XI' C_class = [ [[True, False], [False, False]], # 'ZI' [[False, True], [False, False]], # 'XI' [[True, True], [False, False]], ] # 'YI' # Class D - canonical representative is 'IZ' D_class = [ [[False, False], [False, True]], # 'IX' [[False, False], [True, False]], # 'IZ' [[False, False], [True, True]], ] # 'IY' # Class E - only 'II' E_class = [[[False, False], [False, False]]] # 'II' def _from_pair_paulis_to_type(pauli_x, pauli_z, qubit): """Converts a pair of Paulis pauli_x and pauli_z into a type""" type_x = [pauli_x.z[qubit], pauli_x.x[qubit]] type_z = [pauli_z.z[qubit], pauli_z.x[qubit]] return [type_x, type_z] def _compute_greedy_cost(list_pairs): """Compute the CNOT cost of one step of the algorithm""" A_num = 0 B_num = 0 C_num = 0 D_num = 0 for pair in list_pairs: if pair in A_class: A_num += 1 elif pair in B_class: B_num += 1 elif pair in C_class: C_num += 1 elif pair in D_class: D_num += 1 if (A_num % 2) == 0: raise QiskitError("Symplectic Gaussian elimination fails.") # Calculate the CNOT cost cost = 3 * (A_num - 1) / 2 + (B_num + 1) * (B_num > 0) + C_num + D_num if list_pairs[0] not in A_class: # additional SWAP cost += 3 return cost def _calc_decoupling(pauli_x, pauli_z, qubit_list, min_qubit, num_qubits, cliff): """Calculate a decoupling operator D: D^{-1} * Ox * D = x1 D^{-1} * Oz * D = z1 and reduce the clifford such that it will act trivially on min_qubit """ circ = QuantumCircuit(num_qubits) # decouple_cliff is initialized to an identity clifford decouple_cliff = cliff.copy() num_qubits = decouple_cliff.num_qubits decouple_cliff.phase = np.zeros(2 * num_qubits) decouple_cliff.symplectic_matrix = np.eye(2 * num_qubits) qubit0 = min_qubit # The qubit for the symplectic Gaussian elimination # Reduce the pair of Paulis to a representative in the equivalence class # ['XZ', 'XX', 'XI', 'IZ', 'II'] by adding single-qubit gates for qubit in qubit_list: typeq = _from_pair_paulis_to_type(pauli_x, pauli_z, qubit) if typeq in [ [[True, True], [False, False]], # 'YI' [[True, True], [True, True]], # 'YY' [[True, True], [True, False]], ]: # 'YZ': circ.s(qubit) _append_s(decouple_cliff, qubit) elif typeq in [ [[True, False], [False, False]], # 'ZI' [[True, False], [True, False]], # 'ZZ' [[True, False], [False, True]], # 'ZX' [[False, False], [False, True]], ]: # 'IX' circ.h(qubit) _append_h(decouple_cliff, qubit) elif typeq in [ [[False, False], [True, True]], # 'IY' [[True, False], [True, True]], ]: # 'ZY' circ.s(qubit) circ.h(qubit) _append_s(decouple_cliff, qubit) _append_h(decouple_cliff, qubit) elif typeq == [[True, True], [False, True]]: # 'YX' circ.h(qubit) circ.s(qubit) _append_h(decouple_cliff, qubit) _append_s(decouple_cliff, qubit) elif typeq == [[False, True], [True, True]]: # 'XY' circ.s(qubit) circ.h(qubit) circ.s(qubit) _append_s(decouple_cliff, qubit) _append_h(decouple_cliff, qubit) _append_s(decouple_cliff, qubit) # Reducing each pair of Paulis (except of qubit0) to 'II' # by adding two-qubit gates and single-qubit gates A_qubits = [] B_qubits = [] C_qubits = [] D_qubits = [] for qubit in qubit_list: typeq = _from_pair_paulis_to_type(pauli_x, pauli_z, qubit) if typeq in A_class: A_qubits.append(qubit) elif typeq in B_class: B_qubits.append(qubit) elif typeq in C_class: C_qubits.append(qubit) elif typeq in D_class: D_qubits.append(qubit) if len(A_qubits) % 2 != 1: raise QiskitError("Symplectic Gaussian elimination fails.") if qubit0 not in A_qubits: # SWAP qubit0 and qubitA qubitA = A_qubits[0] circ.swap(qubit0, qubitA) _append_swap(decouple_cliff, qubit0, qubitA) if qubit0 in B_qubits: B_qubits.remove(qubit0) B_qubits.append(qubitA) A_qubits.remove(qubitA) A_qubits.append(qubit0) elif qubit0 in C_qubits: C_qubits.remove(qubit0) C_qubits.append(qubitA) A_qubits.remove(qubitA) A_qubits.append(qubit0) elif qubit0 in D_qubits: D_qubits.remove(qubit0) D_qubits.append(qubitA) A_qubits.remove(qubitA) A_qubits.append(qubit0) else: A_qubits.remove(qubitA) A_qubits.append(qubit0) # Reduce pairs in Class C to 'II' for qubit in C_qubits: circ.cx(qubit0, qubit) _append_cx(decouple_cliff, qubit0, qubit) # Reduce pairs in Class D to 'II' for qubit in D_qubits: circ.cx(qubit, qubit0) _append_cx(decouple_cliff, qubit, qubit0) # Reduce pairs in Class B to 'II' if len(B_qubits) > 1: for qubit in B_qubits[1:]: qubitB = B_qubits[0] circ.cx(qubitB, qubit) _append_cx(decouple_cliff, qubitB, qubit) if len(B_qubits) > 0: qubitB = B_qubits[0] circ.cx(qubit0, qubitB) circ.h(qubitB) circ.cx(qubitB, qubit0) _append_cx(decouple_cliff, qubit0, qubitB) _append_h(decouple_cliff, qubitB) _append_cx(decouple_cliff, qubitB, qubit0) # Reduce pairs in Class A (except of qubit0) to 'II' Alen = int((len(A_qubits) - 1) / 2) if Alen > 0: A_qubits.remove(qubit0) for qubit in range(Alen): circ.cx(A_qubits[2 * qubit + 1], A_qubits[2 * qubit]) circ.cx(A_qubits[2 * qubit], qubit0) circ.cx(qubit0, A_qubits[2 * qubit + 1]) _append_cx(decouple_cliff, A_qubits[2 * qubit + 1], A_qubits[2 * qubit]) _append_cx(decouple_cliff, A_qubits[2 * qubit], qubit0) _append_cx(decouple_cliff, qubit0, A_qubits[2 * qubit + 1]) return circ, decouple_cliff