Pauli path integral

Main Reference:

• A polynomial-time classical algorithm for noisy quantum circuits

• quantum mean value

• A Polynomial-Time Classical Algorithm for Noisy Random Circuit Sampling

• non-unital …

Talk:

• STOC 2023 by Liu

• Simons by Liu

• QIP by Aharonov

…. todo

PauliPropagation.jl is a Julia package for Pauli propagation simulation of quantum circuits and quantum systems.

MSRudolph/PauliPropagation.jl

Pauli dynamics simulation in Python tbegusic/spd

Pauli algebra

def pauli_compose(a, b, a_num = 1, b_num = 1 ):
    # a_numb, b_num are numbers of bisection applied to a and b respectively
    n = a.num_qubits
    tot = SparsePauliOp(['I' * n], coeffs=[0. + 0.j])
    a_split = [a[i:i+ math.ceil(len(a) / (2**a_num))] for i in range(0, len(a), math.ceil(len(a) / (2**a_num)))]
    print(f'a_split: {len(a_split[0])}')
    print(f'a size: {len(a)}')
    print(f'a_split diff: {(a_split[0] - a).simplify()}')
    b_split = [b[i:i+ math.ceil(len(a) / (2**b_num))] for i in range(0, len(b), math.ceil(len(b) / (2**b_num)))]
    print(f'b_split diff: {(b_split[0] -b).simplify()}')
    for i in a_split:
        print('+', end='')
        for j in b_split:
            tot = tot + (i @ j).simplify()
            tot = tot.simplify()
    return tot

from qiskit.quantum_info import SparsePauliOp
import math, random
import numpy as np

def generate_rand_paulis(num_paulis, num_qubits):
    """
    Generates a list of random Pauli operators.
    """
    paulis, coeffs = [], []
    for _ in range(num_paulis):
        # Generate a random Pauli operator
        pauli_str = ''.join(['I', 'X', 'Y', 'Z'][random.randint(0, 3)] for _ in range(num_qubits))
        coeff = complex(random.uniform(-1, 1), random.uniform(-1, 1))
    #     paulis.append(SparsePauliOp(pauli_str, coeffs=[coeff]))
        paulis.append(pauli_str)
        coeffs.append(coeff)
    # Convert to SparsePauliOp
    pauli_op = SparsePauliOp(paulis, coeffs=coeffs).simplify()
    return pauli_op

num_paulis = 800
num_qubits = 11 
paulis1 = generate_rand_paulis(num_paulis, num_qubits)
paulis2 = generate_rand_paulis(num_paulis, num_qubits)
product = paulis1 @ paulis2
print(f'product size: {len(product)}')
product = product.simplify()
print(f'(#paulis)/(max #paulis): {len(product.coeffs)}/{4**num_qubits}')

product size: 640000
(#paulis)/(max #paulis): 593478/4194304

def pauli_to_dict(pauli_op, verbose=False):
    """
    Convert a SparsePauliOp to a dictionary.
    """
    pauli_op = pauli_op.simplify()
    pauli_labels = pauli_op.paulis.to_labels()
    pauli_coeffs = pauli_op.coeffs
    if verbose: print(pauli_labels, pauli_coeffs)
    pauli_dict = dict(zip(pauli_labels, pauli_coeffs))
    return pauli_dict

def dict_to_pauli(pauli_dict, verbose=True):
    """
    Convert a dictionary to a SparsePauliOp.
    """
    # paulis = list(pauli_dict.keys())
    # coeffs = list(pauli_dict.values())
    ## eliminiate all zeor values in the dictionary
    paulis = []
    coeffs = []
    if verbose: print(f'converting dict to pauli')
    n = len(list(pauli_dict.keys())[0])
    for key, value in pauli_dict.items():
        if value != 0+0.j:
            paulis.append(key)
            coeffs.append(value)

    if len(paulis) == 0:
        return SparsePauliOp(['I' * n], coeffs=[0. + 0.j])
    else:
        if verbose: print(f'pauli_dict (non-zero) size: {len(paulis)}')
        return SparsePauliOp(paulis, coeffs=coeffs)
    # return SparsePauliOp(paulis, coeffs=coeffs).simplify()

def update_pauli_dict(all_pauli_dict, new_pauli_dict):
    result_dict = all_pauli_dict.copy()
    for key, value in new_pauli_dict.items():
        result_dict[key] += value
    return result_dict

def add2pauli_dict(pauli_dict1, pauli_dict2):
    """
    Add two dictionaries of Pauli operators.
    """
    # Create a new dictionary to store the result
    result_dict = pauli_dict1.copy()
    for key, value in pauli_dict2.items():
        if key in result_dict:
            result_dict[key] += value
        else:
            result_dict[key] = value
    return result_dict

test =paulis1[:900]
# print(test)

test_dict = pauli_to_dict(test)
test_conv = dict_to_pauli(test_dict)
print('difference: ', (test - test_conv).simplify())
sum2dict = dict_to_pauli(add2pauli_dict(test_dict, pauli_to_dict(-test)))
print(sum2dict)
print('sum2dict: ', len(sum2dict.simplify()))
pauli_to_dict(SparsePauliOp(['X','X'], coeffs=[1, -2]), verbose=True)

converting dict to pauli
pauli_dict (non-zero) size: 900
difference:  SparsePauliOp(['IIIIIIIIIII'],
              coeffs=[0.+0.j])
converting dict to pauli
SparsePauliOp(['IIIIIIIIIII'],
              coeffs=[0.+0.j])
sum2dict:  1

import itertools
def pauli_dict_compose(a, b, verbose=False):
    n = a.num_qubits
    assert n == b.num_qubits, "Number of qubits must be the same for both operators"
    ## generate all pauli strings of length n each is from [I,X,Y,Z]
    all_paulis = [''.join(p) for p in itertools.product('IXYZ', repeat=n)]
    result_dict = dict(zip(all_paulis, [0. + 0.j] * len(all_paulis)))

    for index, pauli in enumerate(a):
        # print(result_dict)
        if (index+1) % 100 == 0: 
            # if verbose: print(f'#pauli:{len(temp_pauli_dict)}', end='')
            print(f' (i={index})')
        else:
            print('+', end='')
        temp_pauli_dict = pauli_to_dict(pauli @ b)
        result_dict = update_pauli_dict(result_dict, temp_pauli_dict)

    result_pauli = dict_to_pauli(result_dict)
    if verbose: print('#paulis:', len(result_pauli))

    return result_pauli

num_paulis = 8000
num_qubits = 11 
m = 1000
# num_paulis = 2
# num_qubits = 1 
# m = 2
paulis1 = generate_rand_paulis(num_paulis, num_qubits)
paulis2 = generate_rand_paulis(num_paulis, num_qubits)
# print(paulis1[:m])
# print(paulis2[:m])
dict_res = pauli_dict_compose(paulis1[:m], paulis2[:m], verbose=True)
allbyall = (paulis1[:m] @ paulis2[:m])
print('#direct compose:', len(allbyall))
allbyall = allbyall.simplify()
print(len(allbyall), len(dict_res), 4**num_qubits)
print(np.linalg.norm(dict_res), np.linalg.norm(allbyall))
# print('dict_res:', dict_res)
# print('allbyall:', allbyall)
# onebyone = pauli_compose(paulis1[:m], paulis2[:m], verbose=True)
# print(onebyone)
# print('diff ops (pauli_compose):', (onebyone-allbyall).simplify())
print('diff ops:', (dict_res-allbyall).simplify())
print('diff len:', len((dict_res-allbyall).simplify()))
print(np.linalg.norm((dict_res-allbyall).simplify()))

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ (i=99)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ (i=199)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ (i=299)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ (i=399)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ (i=499)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ (i=599)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ (i=699)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ (i=799)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ (i=899)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ (i=999)
converting dict to pauli
pauli_dict (non-zero) size: 889596
#paulis: 889596
#direct compose: 1000000
889596 889596 4194304
30185.407056470758 30185.407056470758
diff ops: SparsePauliOp(['IIIIIIIIIII'],
              coeffs=[0.+0.j])
diff len: 1
0.0

def pauli_compose(a, b, simplify_size=100, verbose=False):
    n = a.num_qubits
    tot = SparsePauliOp(['I' * n], coeffs=[0. + 0.j])

    for index, p in enumerate(a):
        if (index+1) % 100 == 0: 
            if verbose: print(f'#pauli:{len(tot)}', end='')
            print(f' ({index})')
        else:
            print('+', end='')
        tot = tot + (p @ b)
        if (index+1) % simplify_size == 0:
            tot = tot.simplify()
        tot = tot.simplify()
    return tot

m = 1000
onebyone = pauli_compose(paulis1[:m], paulis2[:m], verbose=True)
allbyall = (paulis1[:m] @ paulis2[:m])
print(len(allbyall))
allbyall = allbyall.simplify()
print(len(allbyall), len(onebyone))
print((onebyone-allbyall).simplify())
print(np.linalg.norm(onebyone-allbyall))

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#pauli:97845 (99)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#pauli:194259 (199)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#pauli:288537 (299)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#pauli:380516 (399)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#pauli:470299 (499)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#pauli:557943 (599)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#pauli:643496 (699)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#pauli:727169 (799)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#pauli:809056 (899)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#pauli:888561 (999)
1000000
889310 889310
SparsePauliOp(['IIIIIIIIIII'],
              coeffs=[0.+0.j])
0.0

paulis1[:8*m] @ paulis2[:8*m]

Commuting

from qiskit.quantum_info import SparsePauliOp
op = SparsePauliOp.from_list([("XX", 2), ("YY", 1), ("IZ",2j), ("ZZ",1j)])
print(op.group_commuting())
print(op.group_commuting(qubit_wise=True))
print(op.noncommutation_graph(qubit_wise=True))

Trotter circuit

using PauliPropagation

## number of qubits
nqubits = 32

## define the observable
# here I...IZI...I
observable = PauliString(nqubits, :Z, 16)

## define the circuit
# the number of layers
nlayers = 32

# bricklayertopology is also the default if you don't provide any
topology = bricklayertopology(nqubits; periodic=true)

# a circuit containing RX and RZZ Pauli gates on the topology
# derived from the Trotterization of a transverse field Ising Hamiltonian
circuit = tfitrottercircuit(nqubits, nlayers; topology=topology)

# time step
dt = 0.1
# count the number of parameters
nparams = countparameters(circuit)
# define the parameter vector
parameters = ones(nparams) * dt

## the truncations
# maximum Pauli weight
max_weight = 6
# minimal coefficient magnitude
min_abs_coeff = 1e-4

## propagate through the circuit with our best (and currently only propagation method)
pauli_sum = propagate(circuit, observable, parameters; max_weight=max_weight, min_abs_coeff=min_abs_coeff)

## overlap with the initial state
overlapwithzero(pauli_sum)
# yields 0.154596728241...