SAMPLED and EXACT coefficients in `generate_qpd_weights`: How they are both handled ?

[Cmurilochem]

Dear colleagues.

Wanted to kindly request some clarification regarding the num_samples parameter and how the joint QPD is calculated/estimated within the code.

I understand that if num_samples >= 6^{Ms}, where Ms corresponds to, say, the number of space-like cuts, the resulting joint QPD and all its 6^{Ms} terms/subexperiments will be calculated. So, in principle, if the proper number of shots are considered per subexperiment, the final expectation value reconstructed from the QPD expansion is nearly exact.

However, in the case where num_samples < 6^{Ms}, the terms with probabilities < 1/ num_samples, according to the docstrings, are sampled by Monte Carlo, while the remaining terms are calculated with the "exact" coefficient.

I am interested in understanding this more in details, so I have created a simple example, see below:

import numpy as np

from qiskit import QuantumCircuit
from qiskit.providers import BackendV2 as Backend
from qiskit.quantum_info import SparsePauliOp
from qiskit.providers import BackendV2 as Backend
from qiskit.quantum_info import SparsePauliOp
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

from qiskit_aer import AerSimulator

from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit_ibm_runtime import Batch
from qiskit_ibm_runtime import Session

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_addon_cutting import (
    generate_cutting_experiments,
    partition_problem,
    reconstruct_expectation_values,
)
from qiskit_addon_cutting.qpd import generate_qpd_weights


def simple_quantum_circuit():
    """Fixture to create a simple Qiskit Circuit for testing."""
    qc = QuantumCircuit(2)
    qc.h(0)
    qc.cx(0, 1)
    qc.cx(0, 1)
    qc.h(0)
    return qc


def get_expectation_value_without_circuit_cutting(
    qc: QuantumCircuit,
    observable: SparsePauliOp,
    shots: int,
    backend: Backend,
) -> float:
    """Get the expectation value of a simple quantum circuit without cutting."""
    pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
    isa_circuit = pm.run(qc)
    isa_observable = observable.apply_layout(isa_circuit.layout)

    with Session(backend=backend) as session:
        estimator = Estimator(mode=session, options={"default_shots": shots})
        job = estimator.run([(isa_circuit, isa_observable)])

    results = job.result()
    return results[0].data.evs


if __name__ == "__main__":

    backend = AerSimulator(method="statevector")
    circuit = simple_quantum_circuit()
    observable = SparsePauliOp(["II", "IZ", "ZI", "ZZ"])
    shots = 2**14
    partition_labels = "AB"
    num_samples = 20


    # Partition the circuit into subcircuits based on the partition labels
    partitioned_problem = partition_problem(
        circuit=circuit,
        partition_labels=partition_labels,
        observables=observable.paulis,
    )

    # Extract subcircuits, observables, and QPD bases from the partitioned problem
    subcircuits = partitioned_problem.subcircuits
    subobservables = partitioned_problem.subobservables
    bases = partitioned_problem.bases


    # Prepare subexperiments and calculate joint QPD coefficients
    subexperiments, coefficients = generate_cutting_experiments(
        circuits=subcircuits, observables=subobservables, num_samples=num_samples
    )

    print("Joint QPD coefficients:")
    for i, coeff in enumerate(coefficients):
        print(f"  Coefficient {i}: {coeff}")

    #random_samples = generate_qpd_weights(bases, num_samples=num_samples)
    #print(random_samples)
    #print(sum([value[0] for value in random_samples.values()]))

    # Transpile the subexperiments to ISA circuits
    pass_manager = generate_preset_pass_manager(optimization_level=1, backend=backend)
    isa_subexperiments = {
        label: pass_manager.run(partition_subexpts)
        for label, partition_subexpts in subexperiments.items()
    }

    # Submit each partition's subexperiments to the Qiskit Runtime Sampler primitive, in a single batch so that the jobs will run back-to-back.
    with Batch(backend=backend) as batch:
        sampler = Sampler(
            mode=batch,
        )
        jobs = {
            label: sampler.run(subsystem_subexpts, shots=shots)
            for label, subsystem_subexpts in isa_subexperiments.items()
        }

    # Retrieve results
    results = {label: job.result() for label, job in jobs.items()}

    # Get expectation values for each observable term
    reconstructed_expvals = reconstruct_expectation_values(
        results,
        coefficients,
        subobservables,
    )

    # Reconstruct the total expectation value
    circuit_cut_expval = np.dot(reconstructed_expvals, observable.coeffs).real

    
    # Calculate expectation value without circuit cutting
    expected_expval = get_expectation_value_without_circuit_cutting(
        circuit, observable, shots, backend
    )
    
    print(f"Expectation value without circuit cutting: {expected_expval }")
    print(f"Expectation value with circuit cutting: {circuit_cut_expval }")

Note that, in my environment, I have installed, the following packages versions:

qiskit                                     2.0.1
qiskit-addon-cutting             0.10.0
qiskit-aer                                0.17.0
qiskit-ibm-runtime                 0.39.0

The output of the above run is:

Joint QPD coefficients:
  Coefficient 0: (np.float64(0.9), <WeightType.SAMPLED: 2>)
  Coefficient 1: (np.float64(0.9), <WeightType.SAMPLED: 2>)
  Coefficient 2: (np.float64(0.9), <WeightType.SAMPLED: 2>)
  Coefficient 3: (np.float64(-0.9), <WeightType.SAMPLED: 2>)
  Coefficient 4: (np.float64(0.9), <WeightType.SAMPLED: 2>)
  Coefficient 5: (np.float64(0.45), <WeightType.SAMPLED: 2>)
  Coefficient 6: (np.float64(-0.45), <WeightType.SAMPLED: 2>)
  Coefficient 7: (np.float64(-0.45), <WeightType.SAMPLED: 2>)
  Coefficient 8: (np.float64(0.45), <WeightType.SAMPLED: 2>)
  Coefficient 9: (np.float64(-0.45), <WeightType.SAMPLED: 2>)
  Coefficient 10: (np.float64(0.45), <WeightType.SAMPLED: 2>)
  Coefficient 11: (np.float64(-0.45), <WeightType.SAMPLED: 2>)
  Coefficient 12: (np.float64(0.45), <WeightType.SAMPLED: 2>)
  Coefficient 13: (np.float64(0.45), <WeightType.SAMPLED: 2>)
  Coefficient 14: (np.float64(-0.45), <WeightType.SAMPLED: 2>)

Expectation value without circuit cutting: 4.0
Expectation value with circuit cutting: 2.6943627610802654

As you can see, It looks like that 15-20 terms will be sampled from the joint QPD and calculated using the Sampler primitive. But:

• How are the other <WeightType.EXACT: 1> coefficients are calculated ?
• Where they are defined then ?
• Are we just performing here a MC sampling out of the ~num_samples terms in the current version?

I would appreciate very much your guidance in this matter as this will certainly help other users as well.
Kind regards,
Carlos

[LucaCimino]

Hi,

I am going through your same questions and this is what I understood so far:
Considering the circuit in your example, which has 2 CNOT gates cutted into their quasiprobability distributions, this should lead to the generation of 36 ($= 6^{Ms}$) different 2-qubit circuits. But if the num_samples < 6^{Ms}, instead of generating all the 36 circuits, only n <= num_samples are actually generated.
Which circuits are generated is based on a random sampling performed on the 36 circuits, based on their QPD probabilities.
Considering the quasiprobability decomposition of a channel $\Phi = \sum_i c_i \Phi_i$, a CNOT gate is decomposed into a sum of 6 channels where all the $c_i = \pm 0.5$. So, the probability of each channel is $p_i = |c_i| / \sum_i |c_i| = 0.166667$.
If we have 2 CNOT gates, each of the 36 subcircuits has a probability of $0.02777$. And these are the probabilities on which the Monte Carlo sampling is based on.
The sampling procedure consists of randomly picking one of the 36 circuits for num_samples times. It may happen that a circuit is picked more than one time, and this caused an increased of its output weight. For example, in your code you have num_samples = 20 but only 15 circuits in the output. This is why the first 5 circuits are picked 2 times each during the sampling, indeed their output weight is $\pm 0.9$ compared to $\pm 0.45$ of the others. This random sampling causes the fact that every time the program is executed, the number of circuits generated changes, and also their weights.

Regarding the WeightType.EXACT coefficients: in your example, none are generated because the probability of each individual subcircuit ($1/36 = 0.02777$) is lower than the threshold 1 / num_samples $=1/20 = 0.05$. Since all terms are below this threshold, the entire distribution is handled via Monte Carlo sampling.

Now, my further questions for this discussion are about the performance of this method. The described procedure has the pro of not generating all the possible circuits, but it seems performing bad in the estimation of the expectation value. I still need to run some test, but could be better to always generate all the circuits and then try to optimize the overhead decreasing the number of shot for each circuits? In the literature, the Monte Carlo sampling is usually performed on all the circuits to decide which one to execute, but the coefficients are never changed.

Link to other issues #411 and #284 opened by @garrison about clarification between weights, coefficients and samples.

Best,
Luca

[garrison]

Does this section of the doc help clarify: https://qiskit.github.io/qiskit-addon-cutting/explanation/index.html#sample-weights-in-the-qiskit-addon-for-circuit-cutting ?

The most obvious method is to sample all coefficients (with replacement), but a lower variance can be achieved by first identifying the coefficients with absolute weight greater than 1/num_samples and evaluating those coefficients exactly. The function returns EXACT vs SAMPLED so you can identify which is which. If you call the function multiple times with the same arguments, the EXACT coefficients are expected to always be the same, but not the SAMPLED coefficients.

[LucaCimino]

Yes, thank you, now the doc is much clearer than when I read it the first time.

However, reducing the num_samples only decreases the number of generated subexperiments. The sampling overhead required to achieve a desired accuracy remains the same, because the sum of the abs value of all the coefficients $\sum_i |c_i| = \gamma$ is independent of the number of subexperiments.
Less subexperiments results in higher coefficients which means a higher number of shots for each subexperiment.
In general the number of shots for each subexperiments can be calculated as

samples_shots = [math.ceil(abs(coeff[0]) * SHOTS * gamma) for coeff in coefficients]

where SHOTS is the number of shots required to obtain a given accuracy if running the large circuit.

For balanced quasiprobability distributions, such as CNOT gates, is it effective to use a num_samples value different from np.inf? Or is this useful only for a large number of cuts?

Best,
Luca