Callback function in VQE do not return the intermediate values of the parameters

Question

I tried to save all the intermediate values of the parameters when running the built-in VQE in qiskit. The callback function can save the mean values and counts, however it seems only save the very last set of parameters. Please see the codes below, which is really taken from the link here https://github.com/qiskit-community/qiskit-community-tutorials/blob/master/aqua/vqe_convergence.ipynb.

import numpy as np
import pylab

from qiskit import BasicAer
from qiskit.aqua.operators import WeightedPauliOperator
from qiskit.aqua import QuantumInstance, aqua_globals
from qiskit.aqua.algorithms import VQE, NumPyMinimumEigensolver
from qiskit.aqua.components.initial_states import Zero
from qiskit.aqua.components.optimizers import COBYLA, L_BFGS_B, SLSQP
from qiskit.circuit.library import TwoLocal

pauli_dict = {
    'paulis': [{"coeff": {"imag": 0.0, "real": -1.052373245772859}, "label": "II"},
              {"coeff": {"imag": 0.0, "real": 0.39793742484318045}, "label": "ZI"},
              {"coeff": {"imag": 0.0, "real": -0.39793742484318045}, "label": "IZ"},
              {"coeff": {"imag": 0.0, "real": -0.01128010425623538}, "label": "ZZ"},
              {"coeff": {"imag": 0.0, "real": 0.18093119978423156}, "label": "XX"}
              ]
}

qubit_op = WeightedPauliOperator.from_dict(pauli_dict)


num_qubits = qubit_op.num_qubits
init_state = Zero(num_qubits)
var_form = TwoLocal(num_qubits, 'ry', 'cz', initial_state=init_state)

counts = []
values = []
para = []
def store_intermediate_result(eval_count, parameters, mean, std):
    counts.append(eval_count)
    values.append(mean)
    para.append( parameters )

algo = VQE(qubit_op, var_form, COBYLA(), callback=store_intermediate_result)
backend = BasicAer.get_backend('statevector_simulator')
quantum_instance = QuantumInstance(backend=backend)  
algo_result = algo.run(quantum_instance)

counts = np.asarray(counts)
values = np.asarray(values)
para = np.asarray( para )

# print( 'counts' , counts )
# print( 'values' , values )
print( para) 

This is in fact exactly the same code as in this answer VQE restart from a previous computation (I have also browsed other answers with keyword "callback"), and I can get the intermediate values of the cost function. However, for the parameters, what I got is the following (just copy the first few lines), which seems to be the very last set of parameters. Any ideas which part went wrong here?

[[ 1.82181708 -4.4450779   6.76132624  2.43831885  1.15624184  7.51034512
   6.81040957  3.23410408]
 [ 1.82181708 -4.4450779   6.76132624  2.43831885  1.15624184  7.51034512
   6.81040957  3.23410408]
 [ 1.82181708 -4.4450779   6.76132624  2.43831885  1.15624184  7.51034512
   6.81040957  3.23410408]
 [ 1.82181708 -4.4450779   6.76132624  2.43831885  1.15624184  7.51034512
   6.81040957  3.23410408]
 [ 1.82181708 -4.4450779   6.76132624  2.43831885  1.15624184  7.51034512
   6.81040957  3.23410408] ...

Answer 1

Not an answer to your issue but it is too long for comment so I put it here.

---

Very interesting. One thing for sure, if you print out the parameters at each optimization step, then you can see it.

You can do this by adding a an extra line to the store_intermediate_result function definition.

def store_intermediate_result(eval_count, parameters, mean, std):
    counts.append(eval_count)
    values.append(mean)
    para.append( parameters )
    print(parameters) 

If I do this then I can see the changes in the parameters. For example, with two iterations:

[-4.59856165  5.39733825  0.30296211  5.05318614  5.98066982 -4.51917001
  4.19536837  4.18612614]
[-3.59856165  5.39733825  0.30296211  5.05318614  5.98066982 -4.51917001
  4.19536837  4.18612614]

But if I print-out the para variable then I get the same thing you got:

[[-3.59856165  5.39733825  0.30296211  5.05318614  5.98066982 -4.51917001
   4.19536837  4.18612614]
 [-3.59856165  5.39733825  0.30296211  5.05318614  5.98066982 -4.51917001
   4.19536837  4.18612614]]

Answer 2

In case someone sees this, and would like to have a workaround, the following code is inspired by the answer due to @KAJ226. Essentially, instead of saving the parameters into a list, we print them into a file so that we can extract them if needed. I compare the following three energies, the values saved by the call-back, the values recalculated with the parameters printed out to the file, and the values recalculated with the parameters saved by the call-back. It is shown that the first and the second values agree with each other. I attached the figure below.

import numpy as np
import pylab

from qiskit import BasicAer
from qiskit.aqua.operators import WeightedPauliOperator
from qiskit.aqua import QuantumInstance, aqua_globals
from qiskit.aqua.algorithms import VQE, NumPyMinimumEigensolver
from qiskit.aqua.components.initial_states import Zero
from qiskit.aqua.components.optimizers import COBYLA, L_BFGS_B, SLSQP
from qiskit.circuit.library import TwoLocal

pauli_dict = {
    'paulis': [{"coeff": {"imag": 0.0, "real": -1.052373245772859}, "label": "II"},
              {"coeff": {"imag": 0.0, "real": 0.39793742484318045}, "label": "ZI"},
              {"coeff": {"imag": 0.0, "real": -0.39793742484318045}, "label": "IZ"},
              {"coeff": {"imag": 0.0, "real": -0.01128010425623538}, "label": "ZZ"},
              {"coeff": {"imag": 0.0, "real": 0.18093119978423156}, "label": "XX"}
              ]
}

qubit_op = WeightedPauliOperator.from_dict(pauli_dict)


num_qubits = qubit_op.num_qubits
init_state = Zero(num_qubits)
var_form = TwoLocal(num_qubits, 'ry', 'cz', initial_state=init_state)

counts = []
values = []
para3 = []
with open('out.txt', 'w') as f: 
    pass
def store_intermediate_result(eval_count, parameters, mean, std):
    counts.append(eval_count)
    values.append(mean)
    para3.append( parameters ) 
#     print( parameters )
    with open('out.txt', 'a') as f:
        for item in parameters:
            f.write( "%f " % item )
        f.write( "\n" )
    

algo = VQE(qubit_op, var_form, COBYLA(), include_custom=True , callback=store_intermediate_result)
backend = Aer.get_backend('qasm_simulator')
algo_result = algo.run(backend)

with open('out.txt' ) as f:
    para2 = [line.split() for line in f]        

para = []
for x in para2:
    para.append( [float(y) for y in x ] )


counts = np.asarray(counts)
values = np.asarray(values)

values2 = []
for x in counts:
    algo = VQE(qubit_op, var_form, COBYLA(maxiter= 0), include_custom=True , initial_point = para[x-1])
    backend = Aer.get_backend('qasm_simulator')
    algo_result = algo.run(backend)

    values2.append( algo_result['eigenvalue'] )

values3 = []    
for x in para3:
    algo = VQE(qubit_op, var_form, COBYLA(maxiter= 0), include_custom=True , initial_point = x )
    backend = Aer.get_backend('qasm_simulator')
    algo_result = algo.run(backend)

    values3.append( algo_result['eigenvalue'] )
    
fig,ax = pylab.subplots()
Numval = 60
ax.plot( counts[:Numval]  , values[:Numval] , 'ro' , label ='energy output by call-back')
ax.plot( counts[:Numval]  , values2[:Numval] , 'bx' , label = 'energy recalculated with para print-out (correct)')
ax.plot( counts[:Numval]  , values3[:Numval] , 'k*' , label = 'energy recalculated with para by the calle-back (wrong)')
pylab.legend(loc='upper right')