# Is it possible to measure "aux_operators" in built-in VQE with backend "qasm_simulator"

I am trying to measure the number operator, together with the ground state energy with the built-in VQE on Qiskit. However, if I pick the backend to be Aer.get_backend('qasm_simulator'), it seems to give an error 'numpy.float64' object is not iterable. Please see the minimal working code and the error message below. This is rather confusing, because the backend BasicAer.get_backend('statevector_simulator') works perfectly fine. Thanks for the help!

from qiskit.aqua.algorithms import VQE, NumPyEigensolver

from qiskit.chemistry.components.variational_forms import UCCSD
from qiskit.chemistry.components.initial_states import HartreeFock
from qiskit.chemistry.drivers import PySCFDriver, UnitsType
from qiskit import Aer, BasicAer
from qiskit.chemistry import FermionicOperator
from qiskit.aqua.operators import Z2Symmetries
from qiskit.aqua.components.optimizers import L_BFGS_B

optimizer = L_BFGS_B()
backend = Aer.get_backend('qasm_simulator')
# backend = BasicAer.get_backend('statevector_simulator')

atom='H .0 .0 .0; H .0 .0 0.74'
map_type = 'parity'

driver = PySCFDriver(atom=atom, unit=UnitsType.ANGSTROM, basis='sto3g')
molecule = driver.run()

num_alpha = molecule.num_alpha
num_beta = molecule.num_beta
#     num_particles = molecule.num_alpha + molecule.num_beta
num_particles = [ num_alpha , num_beta ]
num_spin_orbitals = molecule.num_orbitals * 2
ferOp = FermionicOperator(h1=molecule.one_body_integrals, h2=molecule.two_body_integrals)

numOp = ferOp.total_particle_number()

qubitOp = ferOp.mapping(map_type=map_type)

qubitNumOp = numOp.mapping(map_type=map_type )

qubitOp = Z2Symmetries.two_qubit_reduction(qubitOp, num_particles)
qubitNumOp = Z2Symmetries.two_qubit_reduction(qubitNumOp, num_particles)

print('Ground state energy without shift is ' , NumPyEigensolver( qubitOp , k=2 ).run().eigenvalues.real )

init_state = HartreeFock( num_spin_orbitals , num_particles , map_type )
print( 'HF =  ' , init_state.bitstr )

# setup the variational form for VQE

var_form_vqe = UCCSD(
num_orbitals=num_spin_orbitals,
num_particles=num_particles,
initial_state=init_state,
qubit_mapping=map_type ,
two_qubit_reduction = True ,
)

algorithm_vqe = VQE(qubitOp, var_form_vqe, optimizer , aux_operators = [qubitNumOp] )

result_vqe = algorithm_vqe.run(backend)

print( 'para_vqe = ' , result_vqe['optimal_point' ] )

print( 'eigenvalue_vqe = ' , result_vqe['eigenvalue' ].real )

print( result_vqe  )


error message is attached below

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-1-5fb4d9ac0b3c> in <module>
59 algorithm_vqe = VQE(qubitOp, var_form_vqe, optimizer , aux_operators = [qubitNumOp] )
60
---> 61 result_vqe = algorithm_vqe.run(backend)
62
63

~/<redacted>/qiskit/aqua/algorithms/quantum_algorithm.py in run(self, quantum_instance, **kwargs)
68                 self.quantum_instance = quantum_instance
69
---> 70         return self._run()
71
72     @abstractmethod

~/<redacted>/qiskit/aqua/algorithms/minimum_eigen_solvers/vqe.py in _run(self)
425
426         if self.aux_operators:
--> 427             self._eval_aux_ops()
428             # TODO remove when ._ret is deprecated
429             result.aux_operator_eigenvalues = self._ret['aux_ops'][0]

~/<redacted>/qiskit/aqua/algorithms/minimum_eigen_solvers/vqe.py in _eval_aux_ops(self, threshold)
445         # Deal with the aux_op behavior where there can be Nones or Zero qubit Paulis in the list
446         self._ret['aux_ops'] = [None if is_none else [result]
--> 447                                 for (is_none, result) in zip(self._aux_op_nones, aux_op_results)]
448         self._ret['aux_ops'] = np.array([self._ret['aux_ops']])
449

TypeError: 'numpy.float64' object is not iterable
$$$$

• Yes. I don't see why not. If you can evaluate $\langle \psi | H | \psi \rangle$ on qasm_simulator then why not $\langle \psi | O | \psi \rangle$ for some other operators... you just need to make sure the operator get mapped to qubit operator correctly. Dec 15 '20 at 23:57
• @KAJ226 I think the number operator is mapped correctly, because 'statevector_simulator' works fine. Also, if I get rid of the number operator, then "qasm_simulator" also works fine...
– fagd
Dec 16 '20 at 0:02
• interesting. I need to look into it a little more careful but in the mean time I posted a quick way to side-step this issue... probably not the best but just a quick way to side step the issue for the moment :) Dec 16 '20 at 0:57

I agree that the code you provided only work with statevector_simulator. I ran it, get the output of particle number to be 2 on the statevector_simulator.

Interestingly, you can get around this problem by just running another VQE execution with the found optimal parameters and not doing any optimization (by set the max_iteration =0 in the optimizer ), to get the particle number... so essentially you are just calculating $$\langle \psi^{optimal} | O_{particle \ number} |\psi^{optimal} \rangle$$.

The code below does what I just mentioned above. I kept it as you wrote it, just added the VQE energy to the end of what you wrote, and the part to find the particle number from what I described.

from qiskit.aqua.algorithms import VQE, NumPyEigensolver

from qiskit.chemistry.components.variational_forms import UCCSD
from qiskit.chemistry.components.initial_states import HartreeFock
from qiskit.chemistry.drivers import PySCFDriver, UnitsType
from qiskit import Aer, BasicAer
from qiskit.chemistry import FermionicOperator
from qiskit.aqua.operators import Z2Symmetries
from qiskit.aqua.components.optimizers import L_BFGS_B

optimizer = L_BFGS_B()
backend = Aer.get_backend('qasm_simulator')
# backend = BasicAer.get_backend('statevector_simulator')

atom='H .0 .0 .0; H .0 .0 0.74'
map_type = 'parity'

driver = PySCFDriver(atom=atom, unit=UnitsType.ANGSTROM, basis='sto3g')
molecule = driver.run()

num_alpha = molecule.num_alpha
num_beta = molecule.num_beta
#     num_particles = molecule.num_alpha + molecule.num_beta
num_particles = [ num_alpha , num_beta ]
num_spin_orbitals = molecule.num_orbitals * 2
ferOp = FermionicOperator(h1=molecule.one_body_integrals, h2=molecule.two_body_integrals)

numOp = ferOp.total_particle_number()

qubitOp = ferOp.mapping(map_type=map_type)

qubitNumOp = numOp.mapping(map_type=map_type )

qubitOp = Z2Symmetries.two_qubit_reduction(qubitOp, num_particles)
qubitNumOp = Z2Symmetries.two_qubit_reduction(qubitNumOp, num_particles)

print('Ground state energy without shift is ' , NumPyEigensolver( qubitOp , k=2 ).run().eigenvalues.real )

init_state = HartreeFock( num_spin_orbitals , num_particles , map_type )
print( 'HF =  ' , init_state.bitstr )

# setup the variational form for VQE

var_form_vqe = UCCSD(
num_orbitals=num_spin_orbitals,
num_particles=num_particles,
initial_state=init_state,
qubit_mapping=map_type ,
two_qubit_reduction = True ,
)

algorithm_vqe = VQE(qubitOp, var_form_vqe, optimizer ,  include_custom = True )
from qiskit.aqua import QuantumInstance
quantum_instance = QuantumInstance(backend = backend, shots= 10000, optimization_level= 3)
result_vqe = algorithm_vqe.run(quantum_instance)

print( 'para_vqe = ' , result_vqe['optimal_point' ] )
print('VQE energy =', result_vqe['optimal_value'] )
print('eigenstate:', result_vqe['eigenstate'])

print('**************** Finding Particle Number *************************')
from qiskit.aqua.components.optimizers import COBYLA
optimizer = COBYLA(maxiter= 0,tol=0.0001)
initial_point = result_vqe['optimal_point' ]
particle_number_vqe = VQE(qubitNumOp, var_form_vqe, optimizer , include_custom = True, initial_point = initial_point )
particle_number = particle_number_vqe.run(quantum_instance)
print( 'Particle_number = ' , particle_number['optimal_value'] )


The output would be something like:

Ground state energy without shift is  [-1.85238817 -1.2458777 ]
HF =   [False  True]
para_vqe =  [-7.22409784e-09  1.35083082e-08 -1.12782817e-01]
VQE energy = -1.852388173569581
eigenstate: {'01': 9888, '10': 112}
**************** Finding Particle Number *************************
Particle_number =  1.999999999999998

• LOL... That is a really smart work-around. I don't know we can run VQE zero times. Thanks a lot, you really help me twice!
– fagd
Dec 16 '20 at 0:57
• @fagd no problem! Glad I was able to help. :) Dec 16 '20 at 1:01
• One further related question. I notice "qasm_simulator" is significantly slower than "statevector_simulator". Is it something normal?
– fagd
Dec 16 '20 at 1:10
• @fagd I never compare them directly but one thing for sure, if you simulate something really large, like more than 15 qubits (and also depending on what local machine you have but I am speaking in term of average machine) you will find statevector to be a problem as it tries to save all the data of the density matrix, which grows exponentially. In qasm_simulator you can use the snap-shot mode. By setting include_custom = True in your VQE execution and shots = 1 in your quantum instance. This will makes it run much faster!! Dec 16 '20 at 1:18
• quantum_instance = QuantumInstance(backend = backend, shots= 1, optimization_level= 3) and algorithm_vqe = VQE(qubitOp, var_form_vqe, optimizer , include_custom = True )` Dec 16 '20 at 1:20