VQE简单实现

最新推荐文章于 2025-11-12 15:33:40 发布

原创最新推荐文章于 2025-11-12 15:33:40 发布 · 1.5k 阅读

23 ·

CC 4.0 BY-SA版权

文章标签：

#量子力学

实现简单的VQE。
导入包。

from random import random

import numpy as np
import qiskit
from numpy import pi
# importing Qiskit
from qiskit import QuantumCircuit, execute, Aer, IBMQ, QuantumRegister, ClassicalRegister
from qiskit.aqua.algorithms import NumPyEigensolver
from qiskit.aqua.operators import CircuitOp, WeightedPauliOperator
from scipy.optimize import minimize
from qiskit.circuit.library.standard_gates import U2Gate
from qiskit.providers.ibmq import least_busy
from qiskit.quantum_info import Operator
from qiskit.tools.monitor import job_monitor
from qiskit.visualization import plot_histogram, plot_bloch_multivector
import matplotlib.pyplot as plt

初始化一些数据以及输出用数值方法算出的标准下限。

scale = 10
a = random()*scale
b = random()*scale
c = random()*scale
d = random()*scale
# yy = Operator([
#     [1/np.sqrt(2),1/np.sqrt(2)],
#     [complex(0,1)/np.sqrt(2),complex(0,-1)/np.sqrt(2)]
# ])

yy = Operator([
    [1/np.sqrt(2),complex(0,-1)/np.sqrt(2)],
    [1/np.sqrt(2),complex(0,1)/np.sqrt(2)]
])

# qc.hamiltonian(0)

#the standard result(classical) for comparison
def hamiltonian_operator(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)
scale = 10
a, b, c, d = (scale*random(), scale*random(),
              scale*random(), scale*random())
H = hamiltonian_operator(a, b, c, d)
H = hamiltonian_operator(a, b, c, d)
exact_result = NumPyEigensolver(H).run()
reference_energy = min(np.real(exact_result.eigenvalues))
print('The exact ground state energy is: {}'.format(reference_energy))

利用变分电路计算 $<\psi |H|\psi>$
其中 $H$ 是一系列的 $\sigma$ Pauli算符的和。
演示代码的哈密顿量是 $b\sigma _x + c\sigma _y + d\sigma _z$

def U_theta(circuit, parameters):
    circuit.rx(parameters[0],0)
    circuit.ry(parameters[1],0)

def measure(parameters, measure):

    q = QuantumRegister(1)
    classical = ClassicalRegister(1)
    qc = QuantumCircuit(q,classical)
    U_theta(qc,parameters)
    action = a

    if measure == 'X':
        qc.h(0)
        action = b
    elif measure == 'Y':
        # qc.U2Gate(0, np.pi/2)
        qc.unitary(yy, 0, label='yy')
        action = c

    elif measure == 'Z':
        action = d
    elif measure == 'I':
        return a

    qc.measure(0,0)
    # qc.draw(output='mpl')
    # plt.show()
    backend = Aer.get_backend("qasm_simulator")
    shots = 8192
    results = execute(qc, backend=backend, shots=shots).result()
    answer = results.get_counts()

    if('0' in answer):
        res = action*(answer['0'] - (shots - answer['0']))/shots
    else:
        res = -1
    return res

def vqe(parameters):
    adder_res = measure(parameters,'I') + measure(parameters,'X') + measure(parameters, 'Y') + measure(parameters, 'Z')
    return adder_res





parameters = [np.pi, np.pi]
tol = 1e-3 # tolerance for optimization precision.
vqe_result = minimize(vqe, parameters, method="Powell", tol=tol)
print(vqe_result.fun)



# backend = Aer.get_backend("qasm_simulator")
# shots = 1000
# results = execute(qc, backend=backend, shots=shots).result()
# answer = results.get_counts()
# print(answer)
# plot_histogram(answer)
# plt.show()

其中 $a, b, c, d$ 是一系列是初始值，代码中随机生成，电路如下图所示:

最终总代码:

from random import random

import numpy as np
import qiskit
from numpy import pi
# importing Qiskit
from qiskit import QuantumCircuit, execute, Aer, IBMQ, QuantumRegister, ClassicalRegister
from qiskit.aqua.algorithms import NumPyEigensolver
from qiskit.aqua.operators import CircuitOp, WeightedPauliOperator
from scipy.optimize import minimize
from qiskit.circuit.library.standard_gates import U2Gate
from qiskit.providers.ibmq import least_busy
from qiskit.quantum_info import Operator
from qiskit.tools.monitor import job_monitor
from qiskit.visualization import plot_histogram, plot_bloch_multivector
import matplotlib.pyplot as plt

scale = 10
a = random()*scale
b = random()*scale
c = random()*scale
d = random()*scale
# yy = Operator([
#     [1/np.sqrt(2),1/np.sqrt(2)],
#     [complex(0,1)/np.sqrt(2),complex(0,-1)/np.sqrt(2)]
# ])

yy = Operator([
    [1/np.sqrt(2),complex(0,-1)/np.sqrt(2)],
    [1/np.sqrt(2),complex(0,1)/np.sqrt(2)]
])

# qc.hamiltonian(0)

#the standard result(classical) for comparison
def hamiltonian_operator(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)
scale = 10
a, b, c, d = (scale*random(), scale*random(),
              scale*random(), scale*random())
H = hamiltonian_operator(a, b, c, d)
H = hamiltonian_operator(a, b, c, d)
exact_result = NumPyEigensolver(H).run()
reference_energy = min(np.real(exact_result.eigenvalues))
print('The exact ground state energy is: {}'.format(reference_energy))

def U_theta(circuit, parameters):
    circuit.rx(parameters[0],0)
    circuit.ry(parameters[1],0)

def measure(parameters, measure):

    q = QuantumRegister(1)
    classical = ClassicalRegister(1)
    qc = QuantumCircuit(q,classical)
    U_theta(qc,parameters)
    action = a

    if measure == 'X':
        qc.h(0)
        action = b
    elif measure == 'Y':
        # qc.U2Gate(0, np.pi/2)
        qc.unitary(yy, 0, label='yy')
        action = c

    elif measure == 'Z':
        action = d
    elif measure == 'I':
        return a

    qc.measure(0,0)
    # qc.draw(output='mpl')
    # plt.show()
    backend = Aer.get_backend("qasm_simulator")
    shots = 8192
    results = execute(qc, backend=backend, shots=shots).result()
    answer = results.get_counts()

    if('0' in answer):
        res = action*(answer['0'] - (shots - answer['0']))/shots
    else:
        res = -1
    return res

def vqe(parameters):
    adder_res = measure(parameters,'I') + measure(parameters,'X') + measure(parameters, 'Y') + measure(parameters, 'Z')
    return adder_res





parameters = [np.pi, np.pi]
tol = 1e-3 # tolerance for optimization precision.
vqe_result = minimize(vqe, parameters, method="Powell", tol=tol)
print(vqe_result.fun)



# backend = Aer.get_backend("qasm_simulator")
# shots = 1000
# results = execute(qc, backend=backend, shots=shots).result()
# answer = results.get_counts()
# print(answer)
# plot_histogram(answer)
# plt.show()