Question: Does there exist functionality in Qiskit to compute expectation values of ListOps as efficiently as PauliSumOps?

[JoelHBierman]

This could possibly be classified as a feature request, but that depends entirely on whether or not what I am looking for 1. Already exists and I just don't know it. and 2. Whether or not it is even possible. (Because I am unsure of this, for now I have simply classified it as a question.)

My observations:

If you want to evaluate the expectation value of a sum of a large number of Pauli operators, Qiskit can handle this very efficiently. If you want to calculate the expectation values of those Pauli operators in that same summation individually (i.e. perform an analogous set of measurements on the same collection of operators, but not sum them up), then as far as I can tell, this cannot be done efficiently in Qiskit.

Motivation:

The 1 and 2 electron reduced density matrices are important quantities in computational chemistry that require taking expectation values of large numbers of operators separately without summing them up. For example, the number of operator evaluations needed in calculating the 2-RDM scales as O(N^4) where N is the number of qubits if we are assuming no qubit reductions via symmetry considerations. The number of terms in the electronic structure Hamiltonian also scales as O(N^4). Qiskit seems to be able to efficiently calculate electronic structure Hamiltonian expectation values efficiently, but not RDM calculations.

Here is a simple demonstration of what I mean. This is a simple script which constructs the electronic structure Hamiltonian of the LiH molecule in a minimum basis. At the end we define a simple function which takes its expectation value with respect to a UCCSD ansatz circuit with random parameter values.

from time import perf_counter
import numpy as np
from qiskit import transpile
from qiskit_nature.converters.second_quantization import QubitConverter
from qiskit_nature.mappers.second_quantization import JordanWignerMapper
from qiskit_nature.circuit.library import UCCSD, HartreeFock
from qiskit_nature.drivers import PySCFDriver, UnitsType, Molecule
from qiskit_nature.problems.second_quantization import ElectronicStructureProblem
from qiskit.opflow import CircuitStateFn, CircuitSampler, AerPauliExpectation, StateFn, ListOp
from qiskit.providers.aer import AerSimulator
from qiskit.utils import QuantumInstance
from qiskit.utils.backend_utils import is_aer_provider

interatomic_distance = 0.735
basis = 'sto-3g'

backend = AerSimulator(method='statevector', device='CPU')
quantum_instance = QuantumInstance(backend=backend, seed_transpiler=1)

from qiskit_nature.drivers import PySCFDriver, UnitsType, Molecule
molecule = Molecule(geometry=[['Li', [0., 0., 0.]],
                            ['H', [0., 0., interatomic_distance]]],
                    charge=0, multiplicity=1)
driver = PySCFDriver(molecule = molecule, unit=UnitsType.ANGSTROM, basis=basis)
q_molecule = driver.run()
qubit_converter = QubitConverter(JordanWignerMapper())
es_problem = ElectronicStructureProblem(driver=driver)
second_q_op = es_problem.second_q_ops()
qubit_op = qubit_converter.convert(second_q_op[0])

Hf_state = HartreeFock(num_spin_orbitals=2*q_molecule.num_molecular_orbitals,
                        num_particles=es_problem.num_particles,
                        qubit_converter=qubit_converter)

ansatz = UCCSD(num_spin_orbitals=2*q_molecule.num_molecular_orbitals,
                num_particles=es_problem.num_particles,
                qubit_converter=qubit_converter,
                initial_state=Hf_state)
params = np.random.uniform(low=-2*np.pi, high=2*np.pi, size=ansatz.num_parameters)
param_bindings = dict(zip(ansatz.parameters, params.transpose().tolist()))

ansatz = transpile(ansatz, backend=backend)
circuit_op = CircuitStateFn(ansatz)
circuit_sampler = CircuitSampler(quantum_instance, param_qobj=is_aer_provider(quantum_instance.backend))
expectation = AerPauliExpectation()

def measure_expectation(op):

    start_time = perf_counter()
    expect_op = expectation.convert(StateFn(op, is_measurement=True)).compose(circuit_op).reduce()
    sampled_expect_op = circuit_sampler.convert(expect_op, params=param_bindings)
    mean = sampled_expect_op.eval()
    stop_time = perf_counter()
    return stop_time - start_time

We can run the following function calls:

measure_expectation(qubit_op)
measure_expectation(qubit_op[0])

The first function call measures how long it takes to measure the entire Pauli summation. The second measures how long it takes to measure just the first term in this summation. On my local machine these two evaluate to almost exactly the same. (About one second). This makes intuitive sense. What this indicates to me is that the primary bottleneck is not the number of terms in a group of operators being evaluated on the same circuit instance. The bottleneck is primarily other things like sampling the circuit, allocating memory for the circuit object, ect...

This has the implication that if we wanted to calculate the expectation value of a large summed operator like a Hamiltonian, we could do so efficiently. If we wanted to do something like calculate the RDMs for a chemical system, we would have to loop over a huge number of operators and take their expectation values individually, which could could take several orders of magnitude more time than the Hamiltonian evaluation, despite the fact that both tasks involve O(N^4) operators. My guess is that Qiskit sidesteps this vast inefficiency by somehow grouping all of the terms in the summation and associating them with one circuit instance and simulates the circuit once, not O(N^4) times.

A naive guess for how one might replicate this efficiency in the case where we want the expectation values of the individual operators, but not sum them up, would be to group these operators in a similar fashion and take their expectation values:

measure_expectation(ListOp(oplist=[qubit_op[n] for n in range(len(qubit_op))]))

but this does not seem to be the case. This seems to take so long that I just stopped the program when it allocated all of the memory on my machine to this task and does not seem to offer an advantage over just looping over the terms individually. It would seem that Qiskit associates each element in the ListOp with its own separate circuit object and evaluates their expectation values one by one.

The main question: Using existing Qiskit functionality, is it possible to somehow replicate the efficiency of how Qiskit measures expectation values of large PauliSumOps, but with large lists of individual operators? Is it possible for functionality of this kind to be added to Qiskit or is there something which prohibits this? It seems that these two procedures are very close to being the same, except with one you are taking the sum of the expectation values and the other you want to return a list of the individual expectation values.

[Cryoris]

This is an interesting observation @JoelHBierman!

So the reason the expectation value calculation is so fast in your code, and why it takes the same time for the full Hamiltonian as a single Pauli, I'm pretty sure is that you're using AerPauliExpectation and a statevector backend. In this case, Qiskit doesn't need to do any basis transformations but can just compute the Hamiltonian matrix and multiply it with the state. Since the expensive part is computing the state this takes roughly the same time for a sum of Paulis as for a single Pauli.

To perform the basis transformations which you would have to do on a quantum computer, you can use the PauliExpectation and the QasmSimulator (or just the plain AerSimulator). Here you will see that computing a single operator is much much faster than the whole Hamiltonian.

That being said -- there's still an optimization that Qiskit doesn't do yet. If you compute the expectation of the full Hamiltonian as sum of Paulis, Qiskit will group commuting terms, and in your case reduce from 631 to 136 Pauli terms. This is not the case if you evaluate the Pauli terms individually, so the latter will be ~5 times slower. It would be nice if Qiskit would still be able to group the terms if you evaluate the Paulis individually (as ListOp instead of SummedOp/PauliSumOp) and measure 136 circuits instead of 631.

[JoelHBierman]

I see, thank you! This sufficiently answers my question. I'll close this issue for now.