<jupyter_start><jupyter_text>Magic states and T-gadgets with Qiskit  TheoryAssume that we want to prepare (or simulate) $$\rho=U |0\dots0 \rangle_{\rm sys} \langle 0\dots 0|_{\rm sys} U^\dagger,$$where $U$ is an arbitrary circuit compiled into the $H$, $S$, $CNOT$, $T$ gateset.Assume that the number of $T$ gates, also known as T-count, is $t$.Using magic states and T-gadgets, it can be shown [<NAME> (2016)](https://arxiv.org/abs/1601.07601) that:$$\rho = 2^t \, {\rm tr_{aux}}[\rho' (I \otimes |0\dots0 \rangle_{\rm aux} \langle 0\dots 0|_{\rm aux}) ],$$where $$\rho'= U' (|0\dots0 \rangle_{\rm sys}\langle 0\dots 0|_{\rm sys} \otimes |T\dots T \rangle_{\rm aux}  \langle T\dots T|_{\rm aux}) U'^\dagger,$$$|T\dots T \rangle_{\rm aux}  \langle T\dots T|_{\rm aux}$ is the tensor product of $t$ magic states.$$|T\rangle = T\, H\, |0\rangle, \text{ (magic state)}$$  and, **importantly**, $U'$ is a Clifford circuit.More precisely, $U'$ is obtained from $U$ by replacing all $T$ gates with $CNOT$ gates connecting qubits of the system and magic states.**NOTE:** This is just a particular way of using T-gadgets and it corresponds to just a particular post-selected string of results on the auxiliary qubits (the string with $0 \dots 0$).  Setting<jupyter_code>from copy import deepcopy

import qiskit
from qiskit import QuantumCircuit
from qiskit.quantum_info import Pauli

import numpy as np
from matplotlib import pyplot as plt

seed = 42
n_qubits = 8
n_bits = n_qubits
num_t = 9
num_paulis = 4 ** 5
def chop_and_real(x, eps = 1e-10):
    """Chop values of Numpy arrays."""
    y = np.copy(x)
    y[np.abs(x) < eps] = 0
    return np.real(y)

def random_pauli_string(n_qubits, seed=None):
    """Sample a random pauli string. """
    if seed is None:
        rnd_state = np.random
    else:
        rnd_state = np.random.RandomState(seed)  
    rnd_paulis = rnd_state.choice(["I", "X", "Y", "Z"], n_qubits, replace=True)
    return "".join(rnd_paulis)

random_pauli_string(5)
def sample_pauli_strings(num_paulis, n_qubits, with_replacement=False, seed=0):
    """Returns a list of random Pauli strings sorted from many Is to few Is."""
    max_iterations = 4 ** (n_qubits + 2)
    sampled_paulis = []

    for j in range(max_iterations):
        if len(sampled_paulis) == num_paulis:
            sampled_paulis = sorted(
                sampled_paulis, key=lambda x: -sum([int(c=="I") for c in x])
            )
            return sampled_paulis
        
        sampled_paulis.append(random_pauli_string(n_qubits, seed=seed+j))
        
        if not with_replacement:
            sampled_paulis = list(set(sampled_paulis))
            
    raise RuntimeError("Max number of iterations reached.")



sample_pauli_strings(4 ** 1, n_qubits=n_qubits)<jupyter_output><empty_output><jupyter_text>Generate a random Clifford+T circuit.<jupyter_code>def rnd_clifford_t_circuit(n_qubits, n_bits, num_t, t_gadget=False, randomize=False, seed=None, magic_prep=True, seed_s=None):
    """Makes a random circuit with Clifford gates and T gates."""
    
    if n_qubits < 5:
        raise ValueError("The number of qubits must be at least 5.")

    # init rnd number generator
    if seed is None:
        rnd_state = np.random
    else:
        rnd_state = np.random.RandomState(seed)

    if seed is None:
        rnd_state_s = np.random
    else:
        rnd_state_s = np.random.RandomState(seed_s)

    # Init circuit
    if not t_gadget:
        circuit = QuantumCircuit(n_qubits)
    else:
        circuit = QuantumCircuit(n_qubits + num_t)
        # Prepare magic states
        if magic_prep:
            for idx in range(num_t):
                circuit.h(n_qubits + idx)
                circuit.t(n_qubits + idx)

    # Random circuit of CNOT, H, S and T gates
    qubit_indices = list(range(n_qubits))

    # Deterministic layer of Hadamards
    for idx in range(n_qubits):
        circuit.h(idx)

    # Random part of the circuit
    for idx in range(num_t):
        cnot_a, cnot_b, h, s, t = rnd_state.choice(qubit_indices, 5, replace=False)
        circuit.cx(cnot_a, cnot_b)
        circuit.h(h)
        circuit.s(s)
        if not t_gadget:
            circuit.t(t)
        elif not randomize:
            circuit.cnot(t, n_qubits + idx)
        else:
            circuit.cnot(t, n_qubits + idx)
            if rnd_state_s.randint(2): # apply t-gadget adaptive correction
                circuit.s(t)
                circuit.x(n_qubits + idx)

    return circuit<jupyter_output><empty_output><jupyter_text>Let us print a random Clifford+T circuit:<jupyter_code>circuit = rnd_clifford_t_circuit(n_qubits=n_qubits, n_bits=n_bits, num_t=num_t, t_gadget=False, seed=seed)
print(circuit)<jupyter_output>┌───┐┌───┐┌───┐               ┌───┐                                   »
q_0: ┤ H ├┤ H ├┤ H ├────────────■──┤ T ├─────────────────■─────────────────»
     ├───┤└───┘├───┤            │  └───┘                 │                 »
q_1: ┤ H ├──■──┤ S ├────────────┼──────────────■─────────┼────■────────────»
     ├───┤  │  ├───┤            │  ┌───┐┌───┐  │  ┌───┐  │    │  ┌───┐┌───┐»
q_2: ┤ H ├──┼──┤ T ├────────────┼──┤ X ├┤ S ├──┼──┤ T ├──┼────┼──┤ X ├┤ H ├»
     ├───┤  │  └───┘     ┌───┐  │  └─┬─┘├───┤  │  ├───┤  │    │  └─┬─┘└───┘»
q_3: ┤ H ├──┼─────────■──┤ H ├──┼────┼──┤ H ├──┼──┤ T ├──┼────┼────┼───────»
     ├───┤  │  ┌───┐  │  ├───┤  │    │  ├───┤  │  ├───┤  │    │    │  ┌───┐»
q_4: ┤ H ├──┼──┤ S ├──┼──┤ T ├──┼────┼──┤ H ├──┼──┤ T ├──┼────┼────■──┤ T ├»
     ├───┤┌─┴─┐├───┤  │  └───┘  │    │  ├───┤  │  ├───┤  │    │  ┌───┐├───┤»
q_5: ┤ H ├┤ X ├┤ T ├──┼─────────┼────■──┤ T ├──┼──┤ S ├──┼────┼──┤ S ├┤ H ├»
     ├───┤└───┘└───┘  │       ┌─┴─┐┌───┐├───┤  │  └───┘  │  ┌─┴─┐└───┘└───┘»[...]<jupyter_text>Let us print the corresponding T-gadget version (assuming auxiliary qubits are projected on the zero state).<jupyter_code>circuit_with_gadgets = rnd_clifford_t_circuit(
    n_qubits=n_qubits, n_bits=n_bits, num_t=num_t, t_gadget=True, seed=seed
)
# print(circuit_with_gadgets)<jupyter_output><empty_output><jupyter_text>Execute the circuit<jupyter_code>def compute_pauli_expval_with_state(circuit, pauli, backend = None, rho=None):
    """Compute the expectation value of a Pauli observable.
    Returns also the density matrix.
    """
    obs = qiskit.quantum_info.Pauli(pauli).to_matrix()
    if rho is not None:
        return chop_and_real(np.trace(rho @ obs)), rho
  
    circuit_to_measure = circuit.copy()

    if backend is None:
        backend = qiskit.Aer.get_backend('aer_simulator')
        circuit_to_measure.save_density_matrix()

    job = backend.run(circuit_to_measure)
    result = job.result()
      
    if not result.success:
        print("Simulation failed")
    
    if backend.name() == 'aer_simulator_extended_stabilizer':
        # Fake rho just to test speed
        rho = np.eye(2 ** (n_qubits))
    else:
        rho = result.data()["density_matrix"]
    
    return chop_and_real(np.trace(rho @ obs)), rho

pauli_to_measure = "Z"+ "".join((n_qubits - 1) * ["I"])
circuit = rnd_clifford_t_circuit(n_qubits=n_qubits, n_bits=n_bits, num_t=num_t, t_gadget=False, seed=seed)
compute_pauli_expval_with_state(circuit, pauli=pauli_to_measure)[0]<jupyter_output><empty_output><jupyter_text>T-gadget simulation We remove all T gates to get a clifford circuit.<jupyter_code>def compute_pauli_expval_with_gadget(circuit, pauli):
    """Compute the expectation value of a Pauli observable.
    """

    circuit_to_measure = circuit.copy()
    circuit_to_measure.save_density_matrix()

    #backend = AerSimulator(method='extended_stabilizer')
    backend = qiskit.Aer.get_backend('aer_simulator_statevector')
    job = backend.run(circuit_to_measure)
    result = job.result()

    if not result.success:
        print("Simulation failed")

    rho = result.data()["density_matrix"]
    
    extended_pauli = num_t*"I" + pauli
    obs = qiskit.quantum_info.Pauli(extended_pauli).to_matrix()
    
    # Define projector on gadget qubits
    id_on_sys = np.eye(2 ** n_qubits)
    proj_on_zero_aux = np.zeros((2 ** num_t, 2 ** num_t))
    proj_on_zero_aux[0, 0] = 1
    
    gadjet_projector = np.kron(proj_on_zero_aux, id_on_sys)
    return (2 ** num_t) * chop_and_real(np.trace(rho @ obs @ gadjet_projector))

circuit_without_t = rnd_clifford_t_circuit(
    n_qubits=n_qubits, n_bits=n_bits, num_t=num_t, t_gadget=True, randomize=True, seed=seed, seed_s = 1
)
# compute_pauli_expval_with_gadget(circuit_without_t, pauli="Z" + (n_qubits - 1) * "I")<jupyter_output><empty_output><jupyter_text>Linear combination of paulisWe want to compute $E= {\rm tr}(\rho P)$, where $\rho=U |0\dots0 \rangle_{\rm sys} \langle 0\dots 0|_{\rm sys} U^\dagger$ and $P$ is a Pauli operator.The previous code shows that, through T-gadgets, we can evaluate the same expectation value as:$$E= 2^t \, {\rm tr}[\rho' (P \otimes |0\dots0 \rangle_{\rm aux} \langle 0\dots 0|_{\rm aux}) ],$$where $$\rho'= U' (|0\dots0 \rangle_{\rm sys}\langle 0\dots 0|_{\rm sys} \otimes |T\dots T \rangle_{\rm aux}  \langle T\dots T|_{\rm aux}) U'^\dagger,$$$|T\rangle = T\, H\, |0\rangle$ is a magic state and, **importantly**, $U'$ is a Clifford circuit..If the observable is a Pauli operator, can we try to use an Heisenberg-picture-like approach to evaluate its expectation value?This requires expanding final zero-state-projector on the auxiliary qubits as a linear combination of Pauli strings. We use:$$ |0\dots 0 \rangle \langle 0\dots 0| = \prod_{j=1}^{t} \frac{(I_j + Z_j)}{2}=\frac{1}{2^t}\sum_{b \in \{0, 1\}^{\otimes t}} Z^{b_1}Z^{b_2}\dots Z^{b_t}. $$So it is possible to expand $P \otimes |0\dots0 \rangle_{\rm aux} \langle 0\dots 0|_{\rm aux}$ as a linear combination of Pauli strings. Unfortunately the number of terms is $2^t$, so exponential in $t$.<jupyter_code>def paulis_of_zero_projector(num_t, max_z_weight):
    """Returns a list of pauli strings corresponding (up to normalization)
    to the projector on the zero state of n_qubits"""
    paulis = ["I", "Z"]
    for _ in range(num_t - 1):
        new_paulis = []
        for pauli in paulis:
            new_paulis.append(pauli + "I")
            new_paulis.append(pauli + "Z")
        paulis = new_paulis

    # Discard pauli operators with high too many Zs
    truncated_paulis = [p for p in paulis if sum([int(c=="Z") for c in p]) <= max_z_weight]
    
    # Sort from few Zs to many Zs
    return sorted(truncated_paulis, key=lambda x: sum([int(c=="Z") for c in x]))

paulis_of_zero_projector(3, max_z_weight=2)
def paulis_of_zero_projector_evolved(circuit, n_qubits, num_t, max_z_weight):
    """Returns the pauli operators of the magic state projector (on the zero state) 
    evolved by the input circuit (in the Heisenberg picture)."""
    
    # Evolve generators
    generators = paulis_of_zero_projector(num_t, max_z_weight=1)[1:]
    evolved_generators = [Pauli(p + n_qubits * "I").evolve(circuit) for p in generators]
    
    # Multiply evolved generators
    projector_paulis = paulis_of_zero_projector(num_t, max_z_weight)
    evolved = []
    for p in projector_paulis:
        evolved_pauli = Pauli((n_qubits + num_t) * "I")
        for local_p, gen in zip(p, evolved_generators):
            if local_p == "Z":
                evolved_pauli = evolved_pauli.dot(gen)
        evolved.append(evolved_pauli)    
    return evolved

# Example
circuit = rnd_clifford_t_circuit(
    n_qubits=n_qubits, n_bits=n_bits, num_t=num_t, t_gadget=True, randomize=False, magic_prep=False, seed=seed
)
proj_evolved = paulis_of_zero_projector_evolved(circuit, n_qubits, num_t, num_t)
proj_evolved<jupyter_output><empty_output><jupyter_text>Magic state |T> in the Pauli basisThe magic sate $|T\rangle=TH|0\rangle$ has known expectation values with respect to the Pauli basis: $$\langle T | I |T\rangle=1,$$ $$\langle T | Y |T\rangle=\langle T | X |T\rangle=1/\sqrt{2},$$$$\langle T | Z |T\rangle=0.$$We numerically verify this fact.<jupyter_code>t_circuit = QuantumCircuit(1)
t_circuit.h(0)
t_circuit.t(0)

assert np.isclose(compute_pauli_expval_with_state(t_circuit, "I")[0], 1)
assert np.isclose(compute_pauli_expval_with_state(t_circuit, "X")[0], 1 / np.sqrt(2))
assert np.isclose(compute_pauli_expval_with_state(t_circuit, "Y")[0], 1 / np.sqrt(2))
assert np.isclose(compute_pauli_expval_with_state(t_circuit, "Z")[0], 0)<jupyter_output><empty_output><jupyter_text>Linear combination of paulis using Clifford rules<jupyter_code>def expval_with_init_state(pauli, n_qubits, num_t):
    """Evaluates the expectation value of the input Pauli with respect to the initial state of the
    system i.e., |psi> = |0_{1} ... O_{num_qubits}> |T_{1} ... T_{num_t}>.
    """
    
    string = pauli.to_label()
    expval = 1.0

    # Check extract sign from string
    if string[0] == "-":
        expval = -1.0
        string = string[1:]

    # The right part of the string is associated to the system
    sys_string = string[-n_qubits:]
    # The left part of the string is associated to the magic states
    magic_string = string[:num_t]
    assert len(sys_string) + len(magic_string) == len(string)

    # First consider all the cases in which the result is zero
    if "X" in sys_string or "Y" in sys_string or "Z" in magic_string:
        return 0.0
    
    # Evaluate expectation as a product of single-qubit non-zero terms
    for p in magic_string:
        if p in ["X", "Y"]:
            expval *= 1 / np.sqrt(2)
    
    return expval

test_string = "".join(num_t * ["X"] + n_qubits * ["Z"])
assert np.isclose(
    expval_with_init_state(Pauli(test_string), n_qubits=n_qubits, num_t=num_t),
    np.sqrt(2) ** -num_t
)
def compute_pauli_expval_with_gadget_and_linear_comb_heisenberg(circuit, pauli, n_qubits, num_t, max_z_weight, proj_paulis=None):
    """Computes the expectation value of a Pauli observable expanding the projector on the magic states
    as a linear combination of Pauli operators. Each Pauli term is evolved by the circuit
    and evaluated using the Heisenberg picture.
    
    Returns: The expectation value and also individual terms before averaging.
    """
    if proj_paulis is None:
        proj_paulis = paulis_of_zero_projector_evolved(circuit, n_qubits, num_t, max_z_weight)
    
    #before = time()
    p_evolved = Pauli(num_t * "I" + p).evolve(circuit)
    #print("p_evolved time", time() - before)
    #before = time()
    evolved_paulis = [p_evolved.dot(p) for p in proj_paulis]
    #print("evolved_paulis time", time() - before)

    #before = time()
    all_expvals = np.array([expval_with_init_state(p, n_qubits, num_t) for p in evolved_paulis])
    #print("all_expvals time", time() - before)

    return np.sum(all_expvals), all_expvals<jupyter_output><empty_output><jupyter_text>Direct simulation<jupyter_code>sampled_paulis = sample_pauli_strings(num_paulis, n_qubits=n_qubits, with_replacement=False, seed=seed)
print(f"{len(sampled_paulis)} random Paulis sampled.")
%%time

circuit = rnd_clifford_t_circuit(
    n_qubits=n_qubits, n_bits=n_bits, num_t=num_t, t_gadget=False, randomize=False, magic_prep=True, seed=seed
)

heavy_paulis = []
direct_expvals = []
rho = None

for p in sampled_paulis:
    estimate, rho = compute_pauli_expval_with_state(circuit, pauli=p, rho=rho)
    
    # Skip zero expectation values
    #if np.isclose(estimate, 0.0):
    #    continue

    heavy_paulis.append(p)
    direct_expvals.append(estimate)
    # print(f"Pauli: {p} \texpval_heisenberg: {estimate:.4}")<jupyter_output>CPU times: user 16.6 s, sys: 11.8 s, total: 28.4 s
Wall time: 2.39 s<jupyter_text>Qiskit extended stabilizer backend<jupyter_code>%%time

circuit = rnd_clifford_t_circuit(
    n_qubits=n_qubits, n_bits=n_bits, num_t=num_t, t_gadget=False, randomize=False, magic_prep=False, seed=seed
)

heavy_paulis_extended = []
direct_expvals_extended = []

backend = qiskit.Aer.aer.get_backend("aer_simulator_extended_stabilizer")
for p in sampled_paulis:
    estimate = compute_pauli_expval_with_state(circuit, pauli=p, backend=backend)
    
    # Skip zero expectation values
    #if np.isclose(estimate, 0.0):
    #    continue

    heavy_paulis_extended.append(p)
    direct_expvals_extended.append(estimate)
    # print(f"Pauli: {p} \texpval_heisenberg: {estimate:.4}")<jupyter_output>CPU times: user 35.3 s, sys: 38.5 s, total: 1min 13s
Wall time: 6.28 s<jupyter_text>Truncated Heisenberg-only simulation<jupyter_code>from time import time
%%time

truncation_factor = 1.0
max_z_weight = int(num_t * truncation_factor)

circuit = rnd_clifford_t_circuit(
    n_qubits=n_qubits, n_bits=n_bits, num_t=num_t, t_gadget=True, randomize=False, magic_prep=False, seed=seed
)

heavy_paulis_approx = []
heis_expvals_approx = []
all_terms_list = []
# clifford_circuit = qiskit.quantum_info.Clifford(circuit)

# Pre execute Pauli operators of evolved projector
before = time()
proj_paulis = paulis_of_zero_projector_evolved(circuit, n_qubits, num_t, max_z_weight)
print("proj_paulis time:\t", time() - before)


before = time()
for p in heavy_paulis:
    estimate, all_terms = compute_pauli_expval_with_gadget_and_linear_comb_heisenberg(
        circuit, pauli=p, n_qubits=n_qubits, num_t=num_t, max_z_weight=max_z_weight, proj_paulis=proj_paulis,
    )
    
    heavy_paulis_approx.append(p)
    heis_expvals_approx.append(estimate)
    all_terms_list.append(all_terms)
    # print(f"Pauli: {p} \texpval_heisenberg: {estimate:.4}")
print("Expvals for all Paulis time:\t", time() - before)
mean_abs_terms = np.sum(np.abs(np.array(all_terms_list)), axis=0)
plt.plot(mean_abs_terms)<jupyter_output><empty_output><jupyter_text>Comparison between exact and approximated results<jupyter_code># Sort both ideal and approximate results with respect to ideal results
zipped_expvals = np.array(sorted(zip(direct_expvals, heis_expvals_approx), key=lambda x: -abs(x[0])))
direct_expvals_sorted = zipped_expvals[::, 0]
heis_expvals_approx_sorted = zipped_expvals[::, 1]
plt.plot(direct_expvals_sorted, marker = ".", linestyle="None", label="exact")
plt.plot(np.array(heis_expvals_approx_sorted) + 0.02, marker = ".", linestyle="None", label="approx")
plt.title("Comparison between ideal and approximate expectation values")
plt.legend()
plt.show()
# Mean squared error
np.square(np.subtract(direct_expvals_sorted, heis_expvals_approx_sorted)).mean()
# Mean squared error w.r.t. heavy paulis only
non_zero_mask = np.invert(np.isclose(np.abs(direct_expvals_sorted) + np.abs(heis_expvals_approx_sorted), 0.0))
np.sum(non_zero_mask * np.square(np.subtract(direct_expvals_sorted, heis_expvals_approx_sorted)))/np.sum(non_zero_mask)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14445103217808442, 'fraction_numerical': 0.009540236080418261, 'mean_word_length': 3.8800631246712256, 'pattern_counts': {'":': 2, '<': 54, '<?xml version=': 0, '>': 55, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '91653', 'license': 'MIT', 'n_tokens_mistral': 7411, 'n_tokens_neox': 7022, 'n_words': 1679, 'path': 'clifford/magic_state_similation3.ipynb', 'repo_name': 'andreamari/dropbox', 'size': '18561'}