Unexpected results from qiskit.aqua.algorithms.NumpyMinimumEigensolver

[gushuli]

Information

• Qiskit Aqua version: 0.7.3
• Python version: 3.7.7
• Operating system: MacOS Catalina, Version 10.15.5, Macbook Pro, 13 inch, 16GB

What is the current behavior?

The NumPyMinimumEigensolver in qiskit.aqua.algorithms is giving wrong results.

Steps to reproduce the problem

import numpy as np
import pylab
from qiskit.aqua import QuantumInstance
from qiskit.aqua.algorithms import NumPyMinimumEigensolver, VQE
from qiskit.aqua.components.optimizers import SLSQP
from qiskit.chemistry.drivers import PySCFDriver
from qiskit.chemistry.core import Hamiltonian, QubitMappingType

molecule = 'H .0 .0 -{0}; Na .0 .0 {0}'

pts  = [x * 0.1  for x in range(10, 25)]
pts += [x * 0.25 for x in range(10, 18)]
pts += [4.5]

energies = np.empty(len(pts))
distances = np.empty(len(pts))

for i, d in enumerate(pts):
    driver = PySCFDriver(atom=molecule.format(d/2), basis='sto3g')
    qmolecule = driver.run()
    operator =  Hamiltonian(qubit_mapping=QubitMappingType.JORDAN_WIGNER,
                                two_qubit_reduction=False,
                                freeze_core=True, orbital_reduction=[])
    qubit_op, aux_ops = operator.run(qmolecule)
    result = NumPyMinimumEigensolver(qubit_op, aux_operators=aux_ops).run()     
    result = operator.process_algorithm_result(result)
    energies[i] = result.energy
    distances[i] = d

pylab.plot(distances, energies, label="NumpyMinimumEigensolver")
pylab.xlabel('Interatomic distance')
pylab.ylabel('Energy')
pylab.title('NaH Ground State Energy')
pylab.legend(loc='upper right');

Run this code and the results is show in the following figure. The results after distance > 2.3 are not correct

What is the expected behavior?

This picture shows the expected results (the green line).

Suggested solutions

Look into the implementation of NumpyMinimumEigensolver. ``

[gushuli]

@ajavadia Hi Ali, I opened this issue for the bug I reported to you today.

[woodsp-ibm]

I confirm that if I run your code from the command line I get a plot like you show. Yet if I modify the nah_uccsd notebook to just run the NumpyMinimumEigensolver then the plot is as expected. I will need to investigate what is different when running that it behaves like this.

[gushuli]

Hi, thanks for your attention. Yes, I am using the code from nah_uccsd and I run it in Jupytor Notebook 6.0.3. The only change I made is that I changed the encoding from parity to jordan_wigner (and naturally the two-qubit reduction is set to False.) The expected results are from original program with parity encoding. But I believe encoding should not affect the final results here.

Please let me know if you need further information.

[woodsp-ibm]

Hi, some more information for you. If you set two_qubit_reduction to False when using PARITY it shows the same plot as JW mapping. If with JW mapping we set z2symmetry_reduction to 'auto' in the Hamiltonian (new capability in the latest 0.7 Aqua release) then the plot is again as expected.

Two qubit reduction is a type of reduction done by exploiting known symmetry. Z2Symmetry reduction reduces by finding symmetries too. I think what is happening is that, as the molecule is stretched, states are crossing over one another. If you print the result at each step i.e. the one returned process_algorithm_result (has some nice readable format) you will see the number of particles (electrons) measured changes from 2 to 1. When the operators are created, by 2nd quantization, they will end up with states ranging from all orbitals occupied, to all unoccupied. Of course NaH, with the core frozen, and only two particles, we are really wanting the state with two particles (2 orbitals occupied of the right spin too), but the NumpyMinimumEigensolver has no constraints over it and just picks off the lowest eigenvalue. When the symmetry reduction is done, which logically splits the operator up in multiple smaller ones, then the one we use (we know which one to select) seems to no longer contain the other state(s) that result in this crossover (presumably these states are in one of the other smaller operators we discard).

What I also did was to take the first 4 states from the matrix diagonalization that NumpyMinimumEigensolver uses. Here you can see the ground state in orange and another, that is labelled excited state 2 in red. At that distance of 2.3 you can see the states still exist but after the ordering has changed and what we consider the ground state is now what we started off with as Excited State 2.

When run with VQE+UCCSD the outcome of the above is correct since the wavefunction of UCCSD constrains the final state to have 2 particles in it. As NumpyMinimumEigensolver is a generic solver it is just picking the lowest eigenvalue no matter what and we see the behavior you observe.

I think this behavior may just be the 'expected' result with that molecule in this case. I will however try to check with one of our chemistry experts around this molecule's behavior to see if this can be explained too from a chemistry perspective.

[gushuli]

Thank you for the explaination. So I guess that we need to be careful when using the generic NumpyMinimumEigensolver and consider if the actual physical constraints are satisfied.

[woodsp-ibm]

Yes, it advisable to check the result. This is why we added the number of particles, spin etc as extra operators which we measure at for the eigenstate of the algorithm result. Checking is advised for VQE too, especially where a heuristic ansatz is used for the wavefunction where it is unconstrained by aspects such as number of particles to keep it in the proper solution space.

Given the above explanation etc I will close this off now.