VQE parameter bounds

[MariaSapova]

Information

• Qiskit Aqua version: 0.19.0
• Python version: 3.6
• Operating system: Ubuntu

What is the current behavior?

In VQE implementation circuit parameters are bounded automatically (-pi, pi). During optimization process parameter might stuck in the bound even if it is not an optimal point. Below I demonstrate it on the simple example of H2 molecule. If I simply remove the bounds then SLSQP optimizer does not work anymore (parameters go to infinity). I also checked the latest Qiskit version 0.23.0. However, seems like optimization does not work at all: initial [-pi, 1.0] -> optimal [-pi, 1.0], initial [-pi, -1.0] -> optimal [-pi, -0.99] (see the gist)

Steps to reproduce the problem

Run the gist

[woodsp-ibm]

Since the parameters end up in gate rotations having a larger range (bound) will just replicate the core bound - e.g. rotation x is same as 2pi + x. An optimizer might get stuck in a local minima anywhere of course.

Your gist uses the qasm_simulator which will have shot noise which will give the optimizer difficulty - especially if its trying to compute a finite difference gradient. The SPSA optimizer is more suitable in this case. Now I will note that when we first released 0.7.0, by default VQE would use a so-called snapshot instruction method of computing the expectation values which is done by Aer itself. This is much faster, and behaves like the statevector in giving an ideal outcome, but with no need to convert the Paulis to a matrix operator as is normally done when using statevector. But we got feedback that the qasm_simulator was normally chosen by people when they wanted behavior that was more like a real device and expected to see shot noise. So it was changed so the the default behavior is no longer like statevector. Maybe you ran your sample when 0.7.0 was released. To run it the same way use include_custom=True on VQE which will allow this custom expectation using the snapshot to be included in the set VQE will chose from and when using Aer qasm_simulator it will indeed chose this and the result will be like the ideal statevector. More info can be seen in the online docs for VQE

[MariaSapova]

@woodsp-ibm thank you for the answer. I still think that the current implementation has a bug. In my simple H2 example (see the plot), if we fix the first parameter, there are two initial points for the second parameter: -1.0 (green) and 1.0 (red). The optimal values for the second parameter are roughly -2.9 and 3.4. As you can see on the plot, from the -1.0 initial point, gradient descent optimizer can simply decrease the parameter value and reach -2.9 which happens to be optimal. When starting from 1.0, the optimizer wants to increase the second parameter and, ideally, reach 3.4. However, it cannot reach 3.4 because it cannot get past the pi bound. When the parameter equals pi, gradient points toward the increase of it, but the bound restricts further optimization. If, for example, the parameters were cyclic, we would move from pi to -pi and successfully reach -2.9.

[woodsp-ibm]

Currently bounds are informed from the variational forms. I think the general ones we have supply -pi,pi for bound of each parameter so that is passed onto the optimizer. In your example above, with the minimum being reproducible on a cyclic basis if the optimizer does not search the other way from the bound then it will get stuck - as it seems you are seeing.

I will note that when no initial point is given then VQE will generate a random initial point within the bounds (or bounded by -2pi, 2pi if the lower and/or upper bound respectively is None).

So if in general if these minima are cyclic over the parameter space then arguably we would want -2pi to 2pi as bounds and have an initial point bounded by -pi, pi so there is room for it to move. Of course if the initial point values are given still as say 2pi and 2pi would be the point in your picture then things will still get stuck. Removing the bounds would give it freedom to go anywhere but you say the optimization fails then. I guess this is something we need to investigate and think about further.

@Cryoris I tagged you since you have been working on the var forms that are now in the circuit library.

[MariaSapova]

@woodsp-ibm, @Cryoris thank you! I should note that the optimization fails with SLSQP optimizer if I remove bounds completely, but COBYLA works fine in that case. Also for larger systems this issue complicates calculations with heuristic ansatzes, because parameters frequently stuck in the bounds and it is not obvious whether the points -pi or pi are optimal or it is just a failure.

[woodsp-ibm]

With VQE and chemistry some evaluation of the state that is arrived at is possible. The aux operators are used for this. The chemistry module generates operators to measure total spin, num particles and total magnetization so knowing what these should be allows us evaluate the state based on the expectation values produced from these operators. It of course does not prevent things getting stuck but at least gives a way to look at the result and check it.

[Cryoris]

I'm going ahead and close this since I don't think it exists anymore. @MariaSapova feel free to re-open if you still encounter it!

Below is a minimal example to check there are no automatic parameter bounds.

from qiskit.algorithms import VQE
from qiskit.algorithms.optimizers import *
from qiskit.providers.aer import AerSimulator
from qiskit.opflow import X, Z
from qiskit.circuit import QuantumCircuit, Parameter
from qiskit.circuit.library import RealAmplitudes

x = Parameter("x")
# ansatz = QuantumCircuit(1)
# ansatz.ry(0.1 * x, 0)

ansatz = RealAmplitudes(1, reps=0)
ansatz.assign_parameters([0.1 * x], inplace=True)

op = Z # min. eigval: -sqrt(2)

sim = AerSimulator(method="statevector")

opt = SPSA(maxiter=1000)
# opt = L_BFGS_B()
vqe = VQE(ansatz, optimizer=opt, quantum_instance=sim)
res = vqe.compute_minimum_eigenvalue(op)
print(res.optimal_point)  # [31.41...]
print(res.eigenvalue) # -1