Test circuit conversion time for large non-cirq circuits

[natestemen]

In an interview with a Mitiq user this week, we heard that Mitiq can be noticeably slower when using large non-Cirq circuits.

The scope of this issue is to diagnose the extent of this issue. This means answering the question how much more time does it take to use some of Mitiq's core functionality with non-Cirq circuits, than it does with Cirq circuits? As for "core functionality" I suggest we test both ZNE and PEC. Since we only want to test the conversion portion of this pipeline, executors and such can be mocked. In essence, fill out a table like this:

	ZNE	PEC
cirq	time (s)	time (s)
qiskit	time (s)	time (s)

If there are dramatic differences, we will create an issue(s) to help alleviate this problem to the extent that we can.

[purva-thakre]

In an interview with a Mitiq user this week, we heard that Mitiq can be noticeably slower when using large non-Cirq circuits.

Did the user provide an example circuit? What makes a large circuit from their perspective?

I am confused if their circuit had a large depth or if the circuit was built with a large number of qubits.

[natestemen]

Did the user provide an example circuit? What makes a large circuit from their perspective?

No example circuit, and I'm not sure if they meant depth or width, but I imagine both would show the slowdown. Both are worth testing, but depth is what I imagine they were providing feedback on.

[king-p3nguin]

@natestemen
I have created the benchmark script based on the tutorial to measure the circuit conversion time and the simulation time, and the results were as follows:
In the case of PEC, it appears that much of the simulation time is spent on circuit conversion.

Environment

• qiskit version: 1.1.0
• cirq version: 1.3.0
• mitiq version: 0.37.0dev (commit 03b3548)
• Python version: 3.11.7
• Operating system: Windows WSL2 (12th Gen Intel(R) Core(TM) i7-12700F)

Benchmark script:

from time import perf_counter
from mitiq.pec.representations.depolarizing import (
    represent_operations_in_circuit_with_local_depolarizing_noise,
)
from mitiq import benchmarks
import numpy as np
from cirq import DensityMatrixSimulator, depolarize
from mitiq.interface import convert_to_mitiq
from mitiq import zne
from mitiq import pec
from copy import copy
import matplotlib.pyplot as plt
from tqdm.auto import tqdm

SEED = 42
CONVERSION_TIME = 0.0
SIMULATION_TIME = 0.0


def generate_circuit(n_qubits, num_cliffords, return_type):
    return benchmarks.generate_rb_circuits(
        n_qubits=n_qubits,
        num_cliffords=num_cliffords,
        return_type=return_type,
        seed=SEED,
    )[0]


def execute(circuit, noise_level=0.01):
    """Returns Tr[ρ |0⟩⟨0|] where ρ is the state prepared by the circuit
    executed with depolarizing noise.
    """
    start = perf_counter()
    mitiq_circuit, _ = convert_to_mitiq(circuit)
    end = perf_counter()
    global CONVERSION_TIME
    CONVERSION_TIME += end - start

    start = perf_counter()
    noisy_circuit = mitiq_circuit.with_noise(depolarize(p=noise_level))
    rho = DensityMatrixSimulator().simulate(noisy_circuit).final_density_matrix
    end = perf_counter()
    global SIMULATION_TIME
    SIMULATION_TIME += end - start

    return rho[0, 0].real


def benchmark_zne(n_qubits, num_cliffords):
    global CONVERSION_TIME, SIMULATION_TIME
    CONVERSION_TIME, SIMULATION_TIME = 0.0, 0.0

    circuit = generate_circuit(n_qubits, num_cliffords, return_type="qiskit")
    circuit.measure_all()  # only for qiskit

    zne_value = zne.execute_with_zne(circuit, execute)
    ret_conv_time, ret_sim_time = copy(CONVERSION_TIME), copy(SIMULATION_TIME)
    CONVERSION_TIME, SIMULATION_TIME = 0.0, 0.0

    return ret_conv_time, ret_sim_time


def benchmark_pec(n_qubits, num_cliffords):
    global CONVERSION_TIME, SIMULATION_TIME
    CONVERSION_TIME, SIMULATION_TIME = 0.0, 0.0

    circuit = generate_circuit(n_qubits, num_cliffords, return_type="qiskit")

    noise_level = 0.01
    reps = represent_operations_in_circuit_with_local_depolarizing_noise(
        circuit, noise_level
    )
    pec_value = pec.execute_with_pec(circuit, execute, representations=reps)

    ret_conv_time, ret_sim_time = copy(CONVERSION_TIME), copy(SIMULATION_TIME)
    CONVERSION_TIME, SIMULATION_TIME = 0.0, 0.0

    return ret_conv_time, ret_sim_time


def plot_benchmark(n_qubits, num_cliffords, mitigation_type):
    conv_times = []
    sim_times = []
    for i in tqdm(range(len(num_cliffords))):
        if mitigation_type == "ZNE":
            conv_time, sim_time = benchmark_zne(n_qubits, num_cliffords[i])
        elif mitigation_type == "PEC":
            conv_time, sim_time = benchmark_pec(n_qubits, num_cliffords[i])
        conv_times.append(conv_time)
        sim_times.append(sim_time)
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))
    fig.suptitle(f"{mitigation_type} Benchmark for {n_qubits} qubits")
    ax1.plot(num_cliffords, conv_times, label="Qiskit to Cirq Conversion Time")
    ax1.plot(num_cliffords, sim_times, label="Simulation Time")
    ax1.set_xlabel("Number of Cliffords")
    ax1.set_ylabel("Time (s)")
    ax1.legend()
    ax2.bar(
        list(range(len(num_cliffords))),
        np.array(conv_times) / (np.array(conv_times) + np.array(sim_times)),
        tick_label=num_cliffords,
        align="center",
        label="Conversion Time / (Conversion Time + Simulation Time)",
    )
    ax2.set_xlabel("Number of Cliffords")
    ax2.set_ylabel("Conversion Time / (Conversion Time + Simulation Time)")
    ax2.legend()
    plt.show()


n_qubits = 2
num_cliffords = list(range(1, 100, 10))
plot_benchmark(n_qubits, num_cliffords, "ZNE")


n_qubits = 1
num_cliffords = list(range(1, 10, 1))
plot_benchmark(n_qubits, num_cliffords, "PEC")

(I am contributing to this project as a unitaryHACK participant)

[natestemen]

Nice work @king-p3nguin! This is looking really good. Thanks for sharing the results, and the script.

My biggest question is what happens to the CONVERSION_TIME if we start the time after the circuit has already been generated? Since users are mostly using Mitiq with their circuit pre-defined I don't think we should count the amount of time it takes to generate the circuit as part of conversion.

[king-p3nguin]

@natestemen
Thank you for the reply! I just realized that I don't need to convert the circuit in the execute function.

I wrote a benchmark script according to the tutorial because I am using mitiq for the first time, and in the tutorial, the execute function includes circuit conversion process as follows:

https://mitiq.readthedocs.io/en/stable/guide/zne-1-intro.html#:~:text=import%20numpy%20as,%5D.real

import numpy as np
from cirq import DensityMatrixSimulator, depolarize
from mitiq.interface import convert_to_mitiq

def execute(circuit, noise_level=0.01):
    """Returns Tr[ρ |0⟩⟨0|] where ρ is the state prepared by the circuit
    executed with depolarizing noise.
    """
    # Replace with code based on your frontend and backend.
    mitiq_circuit, _ = convert_to_mitiq(circuit)
    noisy_circuit = mitiq_circuit.with_noise(depolarize(p=noise_level))
    rho = DensityMatrixSimulator().simulate(noisy_circuit).final_density_matrix
    return rho[0, 0].real

Maybe that mitiq user was not aware that the code the user had written was not efficient.

Here is the improved benchmark (I just added two lines in the code)

from time import perf_counter
from mitiq.pec.representations.depolarizing import (
    represent_operations_in_circuit_with_local_depolarizing_noise,
)
from mitiq import benchmarks
import numpy as np
from cirq import DensityMatrixSimulator, depolarize
from mitiq.interface import convert_to_mitiq
from mitiq import zne
from mitiq import pec
from copy import copy
import matplotlib.pyplot as plt
from tqdm.auto import tqdm

SEED = 42
CONVERSION_TIME = 0.0
SIMULATION_TIME = 0.0


def generate_circuit(n_qubits, num_cliffords, return_type):
    return benchmarks.generate_rb_circuits(
        n_qubits=n_qubits,
        num_cliffords=num_cliffords,
        return_type=return_type,
        seed=SEED,
    )[0]


def execute(circuit, noise_level=0.01):
    """Returns Tr[ρ |0⟩⟨0|] where ρ is the state prepared by the circuit
    executed with depolarizing noise.
    """
    start = perf_counter()
    circuit, _ = convert_to_mitiq(circuit)
    end = perf_counter()
    global CONVERSION_TIME
    CONVERSION_TIME += end - start

    start = perf_counter()
    noisy_circuit = circuit.with_noise(depolarize(p=noise_level))
    rho = DensityMatrixSimulator().simulate(noisy_circuit).final_density_matrix
    end = perf_counter()
    global SIMULATION_TIME
    SIMULATION_TIME += end - start

    return rho[0, 0].real


def benchmark_zne(n_qubits, num_cliffords):
    global CONVERSION_TIME, SIMULATION_TIME
    CONVERSION_TIME, SIMULATION_TIME = 0.0, 0.0

    circuit = generate_circuit(n_qubits, num_cliffords, return_type="qiskit")
    circuit.measure_all()  # only for qiskit
    circuit, _ = convert_to_mitiq(circuit) ### ADDED ###

    zne_value = zne.execute_with_zne(circuit, execute)
    ret_conv_time, ret_sim_time = copy(CONVERSION_TIME), copy(SIMULATION_TIME)
    CONVERSION_TIME, SIMULATION_TIME = 0.0, 0.0

    return ret_conv_time, ret_sim_time


def benchmark_pec(n_qubits, num_cliffords):
    global CONVERSION_TIME, SIMULATION_TIME
    CONVERSION_TIME, SIMULATION_TIME = 0.0, 0.0

    circuit = generate_circuit(n_qubits, num_cliffords, return_type="qiskit")
    circuit, _ = convert_to_mitiq(circuit) ### ADDED ###

    noise_level = 0.01
    reps = represent_operations_in_circuit_with_local_depolarizing_noise(
        circuit, noise_level
    )
    pec_value = pec.execute_with_pec(circuit, execute, representations=reps)

    ret_conv_time, ret_sim_time = copy(CONVERSION_TIME), copy(SIMULATION_TIME)
    CONVERSION_TIME, SIMULATION_TIME = 0.0, 0.0

    return ret_conv_time, ret_sim_time


def plot_benchmark(n_qubits, num_cliffords, mitigation_type):
    conv_times = []
    sim_times = []
    for i in tqdm(range(len(num_cliffords))):
        if mitigation_type == "ZNE":
            conv_time, sim_time = benchmark_zne(n_qubits, num_cliffords[i])
        elif mitigation_type == "PEC":
            conv_time, sim_time = benchmark_pec(n_qubits, num_cliffords[i])
        conv_times.append(conv_time)
        sim_times.append(sim_time)
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))
    fig.suptitle(f"{mitigation_type} Benchmark for {n_qubits} qubits")
    ax1.plot(num_cliffords, conv_times, label="Qiskit to Cirq Conversion Time")
    ax1.plot(num_cliffords, sim_times, label="Simulation Time")
    ax1.set_xlabel("Number of Cliffords")
    ax1.set_ylabel("Time (s)")
    ax1.legend()
    ax2.bar(
        list(range(len(num_cliffords))),
        np.array(conv_times) / (np.array(conv_times) + np.array(sim_times)),
        tick_label=num_cliffords,
        align="center",
        label="Conversion Time / (Conversion Time + Simulation Time)",
    )
    ax2.set_xlabel("Number of Cliffords")
    ax2.set_ylabel("Conversion Time / (Conversion Time + Simulation Time)")
    ax2.legend()
    plt.show()


n_qubits = 2
num_cliffords = list(range(1, 100, 10))
plot_benchmark(n_qubits, num_cliffords, "ZNE")


n_qubits = 1
num_cliffords = list(range(1, 10, 1))
plot_benchmark(n_qubits, num_cliffords, "PEC")

[natestemen]

Sorry for not catching this earlier, but I think we need to make another modification. When using Mitiq with a non-Cirq frontend the internal workflow will look something like this.

flowchart LR
  subgraph Mitiq
    direction LR
    subgraph Conversions
        direction TB
        i1[convert to cirq] --> f1["manipulate cirq circuit (potentially generating many cirq circuits)"]
        f1 --> f3[convert back to qiskit]
    end
    subgraph QiskitExecutor
    end
  end
  A[Qiskit Circuit] --> Conversions --> QiskitExecutor --> B[Result]

Loading

In particular there are two stages of conversion. I think to get the most accurate numbers, we should have an executor function written in Qiskit. We can keep the circuit that you've been using the same, however.

[king-p3nguin]

@natestemen
Thank you for the clarification.
I have modified the code on my fork to record the execution time of convert_to_mitiq and convert_from_mitiq and ran a benchmark script. The ratio of time spent on circuit conversion to overall simulation time appears to be constant.

Link to diff of my fork: king-p3nguin/mitiq@03b3548...c868b5b

Here is the result:

[natestemen]

Great work here @king-p3nguin! This indeed shows that for non-cirq circuits the workflow is slower because of the conversions. It's not, IMO, horrendous, but as we push into larger regimes (depth >1000) this is potentially more worrying.

Thanks for going back and forth on this a few times. These are great plots to help us figure out where to go from here. Congratulations on claiming the first Mitiq bounty for unitaryHACK 2024!!