When performing a VQE algorithm, the electronic problem Hamiltonian of the physical system under study needs to be mapped to a qubit Hamiltonian written as a sum of tensor products of Pauli operators (X, Y, Z) and the Identity operator forming the so-called Pauli strings. When measuring the expectation value of a Hamiltonian during running VQE calculation, it is desirable to group the list of Pauli strings into commuting subsets in order to minimize the number of measurements. I am trying to understand how this can be achieved in Qiskit. I have looked over the Qiskit code and the conclusions, to me at least, are a bit surprising. Here is the code I have used to run a VQE calculation:
from qiskit.algorithms.optimizers import SLSQP
from qiskit.primitives import Estimator
from qiskit_nature.second_q.algorithms.ground_state_solvers import GroundStateEigensolver, VQEUCCFactory
from qiskit_nature.second_q.circuit.library.ansatzes import UCC
from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.formats.molecule_info import MoleculeInfo
from qiskit_nature.second_q.mappers import JordanWignerMapper, ParityMapper, QubitConverter
from qiskit_nature.second_q.algorithms.initial_points import HFInitialPoint
qubit_converter = QubitConverter(ParityMapper(), two_qubit_reduction=True)
molecule = MoleculeInfo(["H", "H"], [(0.0, 0.0, 0.0), (0.0, 0.0, 0.74)])
driver = PySCFDriver.from_molecule(molecule, basis="sto3g")
electronic_structure_problem = driver.run()
second_quantized_hamiltonian = electronic_structure_problem.second_q_ops()
pauli_sum_operator = qubit_converter.convert(second_quantized_hamiltonian[0], num_particles=electronic_structure_problem.num_particles)
num_particles = electronic_structure_problem.num_particles
num_spatial_orbitals = electronic_structure_problem.num_spatial_orbitals
ucc = UCC(num_spatial_orbitals, num_particles, excitations='sd', qubit_converter=qubit_converter)
vqe_factory = VQEUCCFactory(Estimator(), ucc, SLSQP(), initial_point = HFInitialPoint())
gse = GroundStateEigensolver(qubit_converter, vqe_factory)
# comment out line where solve method is invoked:
# result = gse.solve(electronic_structure_problem)
# the commented line above is actually calling:
main_operator, aux_ops = gse.get_qubit_operators(electronic_structure_problem)
raw_mes_result = gse.solver.compute_minimum_eigenvalue(main_operator, aux_ops)
print("\n grouping type: ", main_operator.grouping_type)
This is pretty much generic Qiskit code. I might be wrong here but, it looks to me that when running this piece of code, no grouping algorithm is applied by default to the Pauli strings. The last print statement in the code above outputs: 'None'. Also, I cannot find an obvious way to select a grouping method. How can one do that? It would be great if somebody could provide an answer along with some pointers to the Qiskit code in order to understand what is actually going on.
After a grouping method is applied when measuring the expectation value of Pauli terms, a rotation must be made to a basis where X and Y operators are changed into Zs in order to perform measurements on a real quantum computer, which typically makes measurements only in the Z basis. From looking at Qiskit code, I see two options for grouping: qubit-wise and a second graph based approach. After using the qubit-wise method, finding a rotation basis is rather easy. The downside is that the number of Pauli groups is relatively large. The other method produces a significantly smaller number of Pauli groups but finding a rotation basis is not so easy. I am wondering what is the computational cost for the second approach? Is it feasible for large molecules (100 orbitals) or is it one of those problems that grow exponentially with the Pauli weight of the list of Pauli strings (the maximum number of non-identity elements in a given spin operator)?
Qiskit code also makes reference to a so-called TPB (Tensor Product Basis) grouping method. What is this and how it relates to the discussion above?
Here is the software version I am using:
