# import packages
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
# first variational form and Z measurement
def get_var_formZ(params):
qr = QuantumRegister(1, name="q")
cr = ClassicalRegister(1, name='c')
qc = QuantumCircuit(qr, cr)
qc.u(params[0], params[1], params[2], qr[0])
qc.measure(qr,cr)
return qc
# variational form and X measurement
def get_var_formX(params):
qr = QuantumRegister(1, name="q")
cr = ClassicalRegister(1, name='c')
qc = QuantumCircuit(qr, cr)
qc.u(params[0], params[1], params[2], qr[0])
qc.h(qr[0])
qc.measure(qr,cr)
return qc
from qiskit.circuit import Parameter
theta1 = Parameter('θ1')
theta2 = Parameter('θ2')
theta3 = Parameter('θ3')
ex_qc=get_var_formX([theta1,theta2,theta3])
ex_qc.draw('mpl')
ex_qc=get_var_formZ([theta1,theta2,theta3])
ex_qc.draw('mpl')
# proability distribution P(0) and P(1)
def get_probability_distribution(counts):
mysum=sum([v for v in counts.values()])
keys=['0','1']
output_distr =[]
for key in keys:
if key not in counts.keys():
output_distr.append(0)
else:
output_distr.append(counts[key]/mysum)
return output_distr
# choose backend to be the simulator
from qiskit import Aer, execute
backend = Aer.get_backend("qasm_simulator")
NUM_SHOTS = 10000
# define energy for Hamiltonian - B[0] sigma_z - B[1] sigma_x
# B is a list of two objects B=[Bz, Bx]
def energy(output_distrZ,output_distrX,B):
cost = -B[0]*(output_distrZ[0] -output_distrZ[1])-B[1]* (output_distrX[0] -output_distrX[1])
return cost
# the objective function is to evaluate the energy cost by running the quantum circuits
def objective_function(params,B):
# Obtain quantum circuit instances from the paramters
qc1 = get_var_formZ(params)
qc2 = get_var_formX(params)
# Execute the quantum circuit to obtain the probability distribution associated with the current parameters
runjob = execute([qc1, qc2],backend,shots=10000)
results=runjob.result()
# Obtain the counts for each measured state, and convert those counts into a probability vector
output_distr1 = get_probability_distribution(results.get_counts(qc1))
output_distr2 = get_probability_distribution(results.get_counts(qc2))
# Calculate the cost as the distance between the output distribution and the target distribution
cost = energy(output_distr1,output_distr2,B)
return cost
import numpy as np
from qiskit.algorithms.optimizers import COBYLA
nshots=10000
# Initialize the COBYLA optimizer via Qiskit
optimizer = COBYLA(maxiter=500, disp=True,tol=1e-6)
Bfield =np.random.rand(2)
print('(Bz,Bx)=',Bfield)
# Create the initial parameters (noting that our single qubit variational form has 3 parameters)
params = np.random.rand(3)
def myobject(params):
fvalue=objective_function(params,Bfield)
return fvalue
counts = []
values = []
myparams = []
def store_intermediate_result(eval_count, parameters, mean):
counts.append(eval_count)
values.append(mean)
myparams.append(parameters)
ret = optimizer.optimize(num_vars=3, objective_function=myobject, initial_point=params)
#
qc1 = get_var_formZ(ret[0])
qc2 = get_var_formX(ret[0])
#from qiskit.visualization import circuit_drawer
#circuit_drawer(qc1)
#circuit_drawer(qc2) cannot draw using the returned circuit!!
counts1 = execute(qc1, backend, shots=nshots).result().get_counts(qc1)
counts2 = execute(qc2, backend, shots=nshots).result().get_counts(qc2)
dst1=get_probability_distribution(counts1)
dst2=get_probability_distribution(counts2)
myEnergy=energy(dst1,dst2,Bfield)
exactEnergy=-np.sqrt(Bfield[0]**2+Bfield[1]**2)
print('Num of shots for each circuit=',nshots)
print('Estimated minimum energy=',myEnergy,' vs. exact=',exactEnergy,' error=',myEnergy-exactEnergy)
print("Sigma_z:", dst1)
print("Sigma_x:", dst2)
print("Parameters Found:", ret[0])
(Bz,Bx)= [0.03150905 0.86548563] Num of shots for each circuit= 10000 Estimated minimum energy= -0.8665165884764569 vs. exact= -0.8660590071530467 error= -0.0004575813234101167 Sigma_z: [0.5246, 0.4754] Sigma_x: [0.9997, 0.0003] Parameters Found: [1.52724859 0.01981231 1.26708084]
print(ret[0],ret[1],ret[2])
[1.52724859 0.01981231 1.26708084] -0.8650872934642916 51
from scipy.optimize import minimize
Bfield =np.random.rand(2)
print('(Bz,Bx)=',Bfield)
counts = []
values = []
myparams = []
def store_intermediate_result(parameters):
print(parameters)
myparams.append(parameters)
params = np.random.rand(3)
#Using directly SciPy optimizer
#ret2=minimize(myobject,params,method='COBYLA',options={'maxiter':500, 'disp':True,'tol':1e-6},callback=store_intermediate_result)
ret2=minimize(myobject,params,method='L-BFGS-B',options={'maxiter':500, 'disp':True,'maxcor': 10, 'ftol': 2.220446049250313e-09, 'gtol': 1e-05, 'eps': 1e-08, 'maxfun': 15000, 'maxiter': 15000, 'iprint': - 1, 'maxls': 20, 'finite_diff_rel_step': None},callback=store_intermediate_result)
(Bz,Bx)= [0.53013298 0.13439277] [0.36752161 0.52344559 0.3260579 ] [0.36752124 0.52344732 0.3260578 ] [0.36752124 0.52344732 0.3260578 ]
print(ret2)
exactEnergy=-np.sqrt(Bfield[0]**2+Bfield[1]**2)
print('Exact energy= ',exactEnergy)
fun: -0.5358098633484562 hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64> jac: array([-166944.56000269, 628522.01440866, -53310.82192104]) message: 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH' nfev: 120 nit: 3 njev: 30 status: 0 success: True x: array([0.36752124, 0.52344732, 0.3260578 ]) Exact energy= -0.5469025410884265
params = np.random.rand(3)
ret2=minimize(myobject,params,method='COBYLA',options={'maxiter':500, 'disp':True,'tol':1e-6},callback=store_intermediate_result)
#ret2=minimize(myobject,params,method='L-BFGS-B',options={'maxiter':500, 'disp':True,'maxcor': 10, 'ftol': 2.220446049250313e-09, 'gtol': 1e-05, 'eps': 1e-08, 'maxfun': 15000, 'maxiter': 15000, 'iprint': - 1, 'maxls': 20, 'finite_diff_rel_step': None},callback=store_intermediate_result)
/home/twei/.local/lib/python3.6/site-packages/scipy/optimize/_minimize.py:543: RuntimeWarning: Method COBYLA does not support callback. warn('Method %s does not support callback.' % method, RuntimeWarning)
print(ret2)
exactEnergy=-np.sqrt(Bfield[0]**2+Bfield[1]**2)
print(exactEnergy)
fun: -0.8058264460947901 maxcv: 0.0 message: 'Optimization terminated successfully.' nfev: 47 status: 1 success: True x: array([0.38331664, 0.20536547, 1.19112403]) -0.8106711884013778
ret2=minimize(myobject,ret2.x,method='L-BFGS-B',options={'maxiter':500, 'disp':True,'maxcor': 10, 'ftol': 2.220446049250313e-09, 'gtol': 1e-05, 'eps': 1e-08, 'maxfun': 15000, 'maxiter': 15000, 'iprint': - 1, 'maxls': 20},callback=store_intermediate_result)
[0.38331678 0.20536547 1.19112406] [0.38331676 0.20536549 1.19112406] [0.38331676 0.20536549 1.19112405]
print(ret2)
exactEnergy=-np.sqrt(Bfield[0]**2+Bfield[1]**2)
print(exactEnergy)
fun: -0.8117173256168172 hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64> jac: array([257313.89019313, -73373.54268975, -67806.98990773]) message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH' nfev: 96 nit: 3 njev: 24 status: 0 success: True x: array([0.38331676, 0.20536549, 1.19112405]) -0.8106711884013778
from qiskit.aqua.operators import WeightedPauliOperator # deprecated
from qiskit.opflow import X, Z, I
from qiskit.algorithms import NumPyEigensolver
#from qiskit.aqua.components.variational_forms import RY, RYRZ, SwapRZ #deprecated
from qiskit.algorithms import VQE
from qiskit.circuit.library import TwoLocal
from qiskit.algorithms.optimizers import COBYLA, L_BFGS_B, SLSQP, SPSA
/home/twei/virtual_qiskit29/lib/python3.8/site-packages/qiskit/aqua/__init__.py:86: DeprecationWarning: The package qiskit.aqua is deprecated. It was moved/refactored to qiskit-terra For more information see <https://github.com/Qiskit/qiskit-aqua/blob/main/README.md#migration-guide> warn_package('aqua', 'qiskit-terra')
def hamiltonian(a, b, c, d):
"""
Creates a*I + b*Z + c*X + d*Y pauli sum
that will be our Hamiltonian operator.
"""
pauli_dict = {
'paulis': [{"coeff": {"imag": 0.0, "real": a}, "label": "I"},
{"coeff": {"imag": 0.0, "real": b}, "label": "Z"},
{"coeff": {"imag": 0.0, "real": c}, "label": "X"},
{"coeff": {"imag": 0.0, "real": d}, "label": "Y"}
]
}
return WeightedPauliOperator.from_dict(pauli_dict)
def BzBxH(B):
return B[0]*Z + B[1]*X
H=BzBxH(Bfield)
def store_intermediate_result(eval_count, parameters, mean, std):
counts.append(eval_count)
values.append(mean)
params.append(parameters)
deviation.append(std)
depth=1
#var_name='RYRZ'
#var_form = RYRZ(H.num_qubits, depth=depth) # deprecated
var_form=TwoLocal(H.num_qubits,rotation_blocks=['ry','rz'],entanglement=None,reps=0)
print(var_form)
optimizer = COBYLA(maxiter=500, disp=True,tol=1e-6)
counts=[]
values=[]
params=[]
deviation=[]
vqe = VQE(var_form, optimizer, quantum_instance=backend, callback=store_intermediate_result)
┌──────────┐┌──────────┐ q_0: ┤ Ry(θ[0]) ├┤ Rz(θ[1]) ├ └──────────┘└──────────┘
/home/twei/virtual_qiskit29/lib/python3.8/site-packages/qiskit/aqua/components/optimizers/optimizer.py:49: DeprecationWarning: The package qiskit.aqua.components.optimizers is deprecated. It was moved/refactored to qiskit.algorithms.optimizers (pip install qiskit-terra). For more information see <https://github.com/Qiskit/qiskit-aqua/blob/main/README.md#migration-guide> warn_package('aqua.components.optimizers',
#backend = BasicAer.get_backend('statevector_simulator')
vqe_result = vqe.compute_minimum_eigenvalue(H)
print(vqe_result)
{ 'aux_operator_eigenvalues': None, 'cost_function_evals': 46, 'eigenstate': {'0': 0.0625, '1': 0.998044963916957}, 'eigenvalue': (-0.5233664536891796+0j), 'optimal_parameters': { ParameterVectorElement(θ[1]): -1.3906231144085042, ParameterVectorElement(θ[0]): 3.2795087358069264}, 'optimal_point': array([ 3.27950874, -1.39062311]), 'optimal_value': -0.5233664536891796, 'optimizer_evals': None, 'optimizer_time': 12.553154230117798}
exact_result = NumPyEigensolver().compute_eigenvalues(operator=H)
exact_energy = min(np.real(exact_result.eigenvalues))
print('VQE result: ',vqe_result.optimal_value,' vs Exact: ',exact_energy)
VQE result: -0.5233664536891796 vs Exact: -0.5469025410884265
# It allows to use callback to keep track of convergence
import matplotlib.pyplot as plt
plt.plot(counts, abs(exact_energy - values))
plt.xlabel('Eval count')
plt.ylabel('Energy difference from solution reference value')
plt.title('Energy convergence')
plt.yscale('log')
plt.legend(loc='upper right')
No handles with labels found to put in legend.
<matplotlib.legend.Legend at 0x7ff10d012c10>
import qiskit
qiskit.__qiskit_version__
{'qiskit-terra': '0.18.2', 'qiskit-aer': '0.9.0', 'qiskit-ignis': '0.6.0', 'qiskit-ibmq-provider': '0.16.0', 'qiskit-aqua': '0.9.5', 'qiskit': '0.30.0', 'qiskit-nature': None, 'qiskit-finance': None, 'qiskit-optimization': None, 'qiskit-machine-learning': None}