I tried to make a cirq program calculating an eigenvalue of the observable by VQE. Inspired by the qulacs VQE tutorial, I defined a cost function from the expectation value of the observable with hardware efficient ansatz.
import cirq
import sympy
import numpy as np
def hea_multilayer(qubits, depth):
n_qubits = len(qubits)
n_params = 2*n_qubits*(depth + 1)
theta = sympy.symbols(f'a:{n_params}')
for j in range(depth):
for i in range(n_qubits):
yield cirq.ry(theta[2*j*n_qubits + i])(qubits[i])
yield cirq.rz(theta[2*j*n_qubits + i + n_qubits])(qubits[i])
for i in range(n_qubits - 1):
yield cirq.CNOT(qubits[i], qubits[i+1])
for i in range(n_qubits):
yield cirq.ry(theta[2*depth*n_qubits + i])(qubits[i])
yield cirq.rz(theta[2*depth*n_qubits + n_qubits + i])(qubits[i])
n_qubits = 2
depth = 1
qubits = cirq.LineQubit.range(n_qubits)
qc = cirq.Circuit()
qc.append(hea_multilayer(qubits, depth))
observable = cirq.PauliSum.from_pauli_strings([
cirq.PauliString(1.0, cirq.I(qubits[0]), cirq.I(qubits[1])),
cirq.PauliString(1.0, cirq.X(qubits[0]), cirq.X(qubits[1])),
cirq.PauliString(-1.0, cirq.Y(qubits[0]), cirq.Y(qubits[1])),
cirq.PauliString(1.0, cirq.Z(qubits[0]), cirq.Z(qubits[1]))
])
print(observable)
print(qc)
sim = cirq.Simulator()
def cost(parameter):
subs = {}
for i in range(2*n_qubits*(depth+1)):
subs[f'a{i}'] = parameter[i]
result = sim.simulate_expectation_values(qc, observable, subs)
return result[0].real
After that, I minimized this cost function using scipy minimize function.
np.random.seed(2023)
parameter = np.random.random(8)
from scipy.optimize import minimize
cost_history = []
cost_history.append(cost(parameter))
min_result = minimize(cost, parameter, method="BFGS", callback=lambda x: cost_history.append(cost(x)))
print(min_result)
print(cost_history)
However, this minimization did not work with the following message:
message: Desired error not necessarily achieved due to precision loss.
success: False
status: 2
fun: 0.2531989514827728
x: [ 3.220e-01 8.904e-01 5.881e-01 -1.969e-01 -6.673e-01
6.296e-01 2.209e-02 2.421e-01]
nit: 7
jac: [ 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
2.000e+00 0.000e+00 0.000e+00]
hess_inv: [[ 1.000e+00 0.000e+00 ... 0.000e+00 0.000e+00]
[ 0.000e+00 1.000e+00 ... 0.000e+00 0.000e+00]
...
[ 0.000e+00 0.000e+00 ... 2.805e-06 -5.530e-06]
[ 0.000e+00 0.000e+00 ... -5.530e-06 1.104e-05]]
nfev: 993
njev: 109
[0.36777588725090027, 0.25321047008037567, 0.2532104179263115, 0.25321032106876373, 0.2532094195485115, 0.25320935994386673, 0.2532089278101921, 0.2531989514827728]
I implemented this program using qiskit, and compared the result to the qiskit VQE solver:
from qiskit.primitives import Estimator
from qiskit.quantum_info import SparsePauliOp
from qiskit.circuit.library import EfficientSU2
import numpy as np
np.random.seed(2023)
qc = EfficientSU2(2, reps = 1, insert_barriers=True, flatten=True)
observable = SparsePauliOp(["II", "XX", "YY", "ZZ"], coeffs = [1, 1, -1, 1])
estimator = Estimator()
parameter = np.random.random(8)
job = estimator.run(qc, observable, parameter_values = parameter)
print(f'Parameter: {parameter}')
print(f'Circuit:\n{qc}')
print(f'Job result: {job.result()}')
def cost(parameter):
job = estimator.run(qc, observable, parameter_values=parameter)
return job.result().values[0]
from scipy.optimize import minimize
cost_history = []
cost_history.append(cost(parameter))
min_result = minimize(cost, parameter, method="L-BFGS-B", callback=lambda x: cost_history.append(cost(x)))
print(min_result)
print(cost_history)
from qiskit.algorithms.minimum_eigensolvers import VQE
from qiskit.algorithms.optimizers import L_BFGS_B
cost_history = []
cost_history.append(cost(parameter))
vqe = VQE(estimator=estimator, ansatz=qc, optimizer=L_BFGS_B(), initial_point=parameter)
result = vqe.compute_minimum_eigenvalue(observable)
print(result)
print(cost_history)
And my result is similar to the qiskit VQE solver result.
Job result: EstimatorResult(values=array([1.43214656]), metadata=[{}])
message: CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
success: True
status: 0
fun: 5.907418998418734e-11
x: [ 2.681e-01 1.710e+00 8.852e-01 1.201e-02 -1.934e-01
1.231e+00 3.499e-01 1.055e+00]
nit: 6
jac: [-2.187e-06 -9.281e-06 -2.776e-06 -4.619e-06 6.217e-06
-9.692e-06 -1.099e-06 -1.099e-06]
nfev: 63
njev: 7
hess_inv: <8x8 LbfgsInvHessProduct with dtype=float64>
[1.432146558079639, 0.2513704685300491, 0.15156282475910865, 0.0019124923978873776, 3.6796306554398583e-05, 1.1992205672939349e-08, 5.907418998418734e-11]
/home/wleelinux/sources/qiskit_expect/hea_estimator.py:31: DeprecationWarning: ``qiskit.algorithms`` has been migrated to an independent package: https://github.com/qiskit-community/qiskit-algorithms. The ``qiskit.algorithms`` import path is deprecated as of qiskit-terra 0.25.0 and will be removed no earlier than 3 months after the release date. Please run ``pip install qiskit_algorithms`` and use ``import qiskit_algorithms`` instead.
from qiskit.algorithms.minimum_eigensolvers import VQE
{ 'aux_operators_evaluated': None,
'cost_function_evals': 63,
'eigenvalue': 2.642741581126984e-11,
'optimal_circuit': <qiskit.circuit.library.n_local.efficient_su2.EfficientSU2 object at 0x7ff1ac903590>,
'optimal_parameters': { ParameterVectorElement(θ[0]): -0.08608117349727416,
ParameterVectorElement(θ[1]): 1.9438040339495675,
ParameterVectorElement(θ[2]): 1.1497492983826083,
ParameterVectorElement(θ[3]): -0.1815167532084752,
ParameterVectorElement(θ[4]): -0.23455036012933897,
ParameterVectorElement(θ[5]): 1.413331393301791,
ParameterVectorElement(θ[6]): 0.7537967436651564,
ParameterVectorElement(θ[7]): 1.4589817775219438},
'optimal_point': array([-0.08608117, 1.94380403, 1.1497493 , -0.18151675, -0.23455036,
1.41333139, 0.75379674, 1.45898178]),
'optimal_value': 2.642741581126984e-11,
'optimizer_evals': None,
'optimizer_result': <qiskit.algorithms.optimizers.optimizer.OptimizerResult object at 0x7ff1a92233d0>,
'optimizer_time': 0.2594926357269287}
Why my cirq code failed? Is it wrong to use simulate_expectation_values method and should I implement the expectation value from the measure of the quantum circuit in cirq?
