Qiskit Analysis

Take a closer look at the behaviors of Qiskit

Qiskit offers functionality of ideal classical simulation using state vectors, shots-based noiseless simulation, and shots-based noisy simulation. Cross entropy loss serves as a good benchmark to gauge how close is an observed circuit output to the ideal classical simulation. The smaller the cross entropy loss, the better does the measured observation approximate the ground truth.

Obviously, a noisy simulation should produce a worse approximation than its noiseless shots-based counterpart. However, it is not so. I simulated a 9-qubit supremacy circuit in incremental batches of shots, and track its cross entropy loss changes. Intuitively, both the noisy and the noiseless simulations should eventually converge after enough shots, and noisy simulation should produce a worse cross entropy loss.

from qcg.generators import gen_supremacy, gen_hwea, gen_BV, gen_qft, gen_sycamore
import numpy as np
from utils.helper_fun import evaluate_circ, get_evaluator_info
import matplotlib.pyplot as plt
import math

def get_xticks(xvals):
    if len(xvals)<=10:
        return xvals
    else:
        x_ticks = []
        step = math.ceil(len(xvals)/10)
        for idx, x in enumerate(xvals):
            if idx%step==0 or idx==len(xvals)-1:
                x_ticks.append(x)
        return x_ticks

def make_plot(noiseless_ce_l,noisy_ce_l,full_circ_size,shots_increment):
    assert len(noiseless_ce_l) == len(noisy_ce_l)
    xvals = range(1,len(noiseless_ce_l)+1)
    plt.figure(figsize=(10,10))
    plt.plot(xvals,noiseless_ce_l,label='noiseless')
    plt.plot(xvals,noisy_ce_l,label='noisy')
    x_ticks = get_xticks(xvals=xvals)
    plt.xticks(ticks=x_ticks,labels=x_ticks)
    plt.ylabel('\u0394H, lower is better')
    plt.xlabel('shots [*%d]'%shots_increment)
    plt.title('%d qubit circuit'%(full_circ_size))
    plt.legend()
    plt.tight_layout()
    plt.show()

def cross_entropy(target,obs):
    assert len(target)==len(obs)
    epsilon = 1e-20
    obs = [abs(x) if x!=0 else epsilon for x in obs]
    sum_of_prob = sum(obs)
    obs = [x/sum_of_prob for x in obs]
    assert abs(sum(obs)-1) < 1e-10
    h = 0
    for p,q in zip(target,obs):
        if p==0:
            h += 0
        else:
            h += p*np.log(p/q)
    return h

def calculate_delta_H(circ,ground_truth,accumulated_prob,counter,shots_increment,evaluation_method):
    if evaluation_method == 'qasm_simulator':
        qasm_evaluator_info = {'num_shots':shots_increment}
        prob_batch = evaluate_circ(circ=circ,backend='noiseless_qasm_simulator',evaluator_info=qasm_evaluator_info)
    elif evaluation_method == 'noisy_qasm_simulator':
        qasm_noise_evaluator_info = get_evaluator_info(circ=circ,device_name='ibmq_boeblingen',
        fields=['device','basis_gates','coupling_map','properties','initial_layout','noise_model'])
        qasm_noise_evaluator_info['num_shots'] = shots_increment
        prob_batch = evaluate_circ(circ=circ,backend='noisy_qasm_simulator',evaluator_info=qasm_noise_evaluator_info)
    else:
        raise Exception('Illegal evaluation method:',evaluation_method)
    accumulated_prob = [(x*(counter-1)+y)/counter for x,y in zip(accumulated_prob,prob_batch)]
    assert abs(sum(accumulated_prob)-1)<1e-10
    accumulated_ce = cross_entropy(target=ground_truth,obs=accumulated_prob)
    return accumulated_ce, accumulated_prob

def circuit_decay(circuit,shots_increment):
    full_circ_size = len(circuit.qubits)
    ground_truth = evaluate_circ(circ=circuit,backend='statevector_simulator',evaluator_info=None)
    noiseless_accumulated_prob = [0 for i in range(np.power(2,full_circ_size))]
    noiseless_delta_H_l = []
    noisy_accumulated_prob = [0 for i in range(np.power(2,full_circ_size))]
    noisy_delta_H_l = []
    max_counter = max(20,int(20*np.power(2,full_circ_size)/shots_increment))
    for counter in range(1,max_counter+1):
        noiseless_accumulated_ce, noiseless_accumulated_prob = calculate_delta_H(circ=circuit,ground_truth=ground_truth,
        accumulated_prob=noiseless_accumulated_prob,counter=counter,shots_increment=shots_increment,evaluation_method='qasm_simulator')
        noiseless_delta_H_l.append(noiseless_accumulated_ce)

        noisy_accumulated_ce, noisy_accumulated_prob = calculate_delta_H(circ=circuit,ground_truth=ground_truth,
        accumulated_prob=noisy_accumulated_prob,counter=counter,shots_increment=shots_increment,evaluation_method='noisy_qasm_simulator')
        noisy_delta_H_l.append(noisy_accumulated_ce)
        
    return noiseless_delta_H_l, noisy_delta_H_l

circ = gen_supremacy(3,3,8)
noiseless_delta_H_l, noisy_delta_H_l = circuit_decay(circuit=circ,shots_increment=1024)
make_plot(noiseless_ce_l=noiseless_delta_H_l,noisy_ce_l=noisy_delta_H_l,full_circ_size=9,shots_increment=1024)
Generating 3x3, 1+8+1 supremacy circuit

png

The plot indicates that the noisy simulation is not always worse than noiseless.

Avatar
Wei Tang
Computer Science PhD Candidate