Mga Halimbawa at Aplikasyon
Sa araling ito, tatalakayin natin ang ilang halimbawa ng variational algorithm at kung paano ito ilalapat:
- Paano gumawa ng custom variational algorithm
- Paano gamitin ang variational algorithm para mahanap ang pinakamababang eigenvalue
- Paano gamitin ang mga variational algorithm para solusyunan ang mga praktikal na kaso ng paggamit
Tandaan na ang Qiskit patterns framework ay maaaring ilapat sa lahat ng problemang ipinapakilala natin dito. Gayunpaman, para maiwasan ang paulit-ulit, isa lamang na halimbawa ang ating tahasang itatampok sa framework na ito β ang isa na pinapatakbo sa tunay na hardware.
Mga Kahulugan ng Problemaβ
Isipin na gusto nating gumamit ng variational algorithm para mahanap ang eigenvalue ng sumusunod na observable:
Ang observable na ito ay may mga sumusunod na eigenvalue:
At mga eigenstate:
# Added by doQumentation β required packages for this notebook
!pip install -q numpy qiskit qiskit-ibm-runtime rustworkx scipy
from qiskit.quantum_info import SparsePauliOp
observable_1 = SparsePauliOp.from_list([("II", 2), ("XX", -2), ("YY", 3), ("ZZ", -3)])
Custom VQEβ
Una, tatalakayin natin kung paano manu-manong buuin ang isang VQE instance para mahanap ang pinakamababang eigenvalue ng . Gagamitin natin dito ang iba't ibang pamamaraan na tinalakay sa buong kursong ito.
def cost_func_vqe(params, ansatz, hamiltonian, estimator):
"""Return estimate of energy from estimator
Parameters:
params (ndarray): Array of ansatz parameters
ansatz (QuantumCircuit): Parameterized ansatz circuit
hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
estimator (Estimator): Estimator primitive instance
Returns:
float: Energy estimate
"""
pub = (ansatz, hamiltonian, params)
cost = estimator.run([pub]).result()[0].data.evs
return cost
from qiskit.circuit.library.n_local import n_local
from qiskit import QuantumCircuit
import numpy as np
reference_circuit = QuantumCircuit(2)
reference_circuit.x(0)
variational_form = n_local(
num_qubits=2,
rotation_blocks=["rz", "ry"],
entanglement_blocks="cx",
entanglement="linear",
reps=1,
)
raw_ansatz = reference_circuit.compose(variational_form)
raw_ansatz.decompose().draw("mpl")

Magsisimula tayo sa pag-debug gamit ang mga lokal na simulator.
from qiskit.primitives import StatevectorEstimator as Estimator
from qiskit.primitives import StatevectorSampler as Sampler
estimator = Estimator()
sampler = Sampler()
Itatatakda natin ngayon ang paunang hanay ng mga parameter:
x0 = np.ones(raw_ansatz.num_parameters)
print(x0)
[1. 1. 1. 1. 1. 1. 1. 1.]
Maaari nating i-minimize ang cost function na ito para makuha ang pinakamainam na mga parameter
# SciPy minimizer routine
from scipy.optimize import minimize
import time
start_time = time.time()
result = minimize(
cost_func_vqe,
x0,
args=(raw_ansatz, observable_1, estimator),
method="COBYLA",
options={"maxiter": 1000, "disp": True},
)
end_time = time.time()
execution_time = end_time - start_time
Return from COBYLA because the trust region radius reaches its lower bound.
Number of function values = 103 Least value of F = -5.999999998357189
The corresponding X is:
[2.27483579e+00 8.37593091e-01 1.57080508e+00 5.82932911e-06
2.49973063e+00 6.41884255e-01 6.33686904e-01 6.33688223e-01]
result
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: -5.999999998357189
x: [ 2.275e+00 8.376e-01 1.571e+00 5.829e-06 2.500e+00
6.419e-01 6.337e-01 6.337e-01]
nfev: 103
maxcv: 0.0
Dahil ang laruan na problemang ito ay gumagamit lamang ng dalawang qubit, maaari nating suriin ito gamit ang linear algebra eigensolver ng NumPy.
from numpy.linalg import eigvalsh
solution_eigenvalue = min(eigvalsh(observable_1.to_matrix()))
print(f"""Number of iterations: {result.nfev}""")
print(f"""Time (s): {execution_time}""")
print(
f"Percent error: {100*abs((result.fun - solution_eigenvalue)/solution_eigenvalue):.2e}"
)
Number of iterations: 103
Time (s): 0.4394676685333252
Percent error: 2.74e-08
Gaya ng makikita, ang resulta ay napakalapit sa ideal.
Pag-eeksperimento para Mapabuti ang Bilis at Katumpakanβ
Magdagdag ng Reference Stateβ
Sa nakaraang halimbawa, hindi tayo gumamit ng anumang reference operator na . Ngayon, pag-isipan natin kung paano makukuha ang ideal eigenstate na . Tingnan ang sumusunod na circuit.
from qiskit import QuantumCircuit
ideal_qc = QuantumCircuit(2)
ideal_qc.h(0)
ideal_qc.cx(0, 1)
ideal_qc.draw("mpl")
Mabilis nating masusuri na ang circuit na ito ay nagbibigay sa atin ng nais na estado.
from qiskit.quantum_info import Statevector
Statevector(ideal_qc)
Statevector([0.70710678+0.j, 0. +0.j, 0. +0.j,
0.70710678+0.j],
dims=(2, 2))
Ngayon na nakita natin kung ano ang hitsura ng isang circuit na naghahanda ng solution state, makatuwiran na gumamit ng Hadamard gate bilang reference circuit, upang ang buong ansatz ay maging:
reference = QuantumCircuit(2)
reference.h(0)
reference.cx(0, 1)
# Include barrier to separate reference from variational form
reference.barrier()
ref_ansatz = variational_form.decompose().compose(reference, front=True)
ref_ansatz.draw("mpl")

Para sa bagong circuit na ito, ang ideal na solusyon ay maaabot kapag lahat ng parameter ay nakatakda sa , kaya nakukumpirma nito na ang pagpili ng reference circuit ay makatuwiran.
Ngayon, ikumpara natin ang bilang ng mga pagsusuri ng cost function, ang mga iterasyon ng optimizer, at ang oras na kinuha kumpara sa nakaraang pagtatangka.
import time
start_time = time.time()
ref_result = minimize(
cost_func_vqe, x0, args=(ref_ansatz, observable_1, estimator), method="COBYLA"
)
end_time = time.time()
execution_time = end_time - start_time
Gamit ang ating pinakamainam na mga parameter para kalkulahin ang pinakamababang eigenvalue:
experimental_min_eigenvalue_ref = cost_func_vqe(
ref_result.x, ref_ansatz, observable_1, estimator
)
print(experimental_min_eigenvalue_ref)
-5.999999996759607
print("ADDED REFERENCE STATE:")
print(f"""Number of iterations: {ref_result.nfev}""")
print(f"""Time (s): {execution_time}""")
print(
f"Percent error: {100*abs((experimental_min_eigenvalue_ref - solution_eigenvalue)/solution_eigenvalue):.2e}"
)
ADDED REFERENCE STATE:
Number of iterations: 127
Time (s): 0.5620882511138916
Percent error: 5.40e-08
Depende sa iyong partikular na sistema, maaaring magresulta o hindi ito sa pagpapabuti ng bilis o katumpakan sa napakaliit na halimbawang ito. Ang punto ay ang pagsisimula sa mga reference state na may pisikang batayan ay nagiging lalong mahalaga sa pagpapabuti ng bilis at katumpakan habang lumalaki ang mga problema.
Baguhin ang Paunang Puntoβ
Ngayon na nakita na natin ang epekto ng pagdaragdag ng reference state, tatalakayin natin kung ano ang nangyayari kapag pumili tayo ng iba't ibang paunang punto na . Sa partikular, gagamitin natin ang at .
Tandaan na, gaya ng tinalakay nang ipinakilala ang reference state, ang ideal na solusyon ay mahahanap kapag lahat ng parameter ay , kaya ang unang paunang punto ay dapat magbigay ng mas kaunting pagsusuri.
import time
start_time = time.time()
x0 = [0, 0, 0, 0, 6, 0, 0, 0]
x0_1_result = minimize(
cost_func_vqe, x0, args=(raw_ansatz, observable_1, estimator), method="COBYLA"
)
end_time = time.time()
execution_time = end_time - start_time
print("INITIAL POINT 1:")
print(f"""Number of iterations: {x0_1_result.nfev}""")
print(f"""Time (s): {execution_time}""")
INITIAL POINT 1:
Number of iterations: 108
Time (s): 0.4492197036743164
Inaayos ang paunang punto sa :
import time
start_time = time.time()
x0 = 6 * np.ones(raw_ansatz.num_parameters)
x0_2_result = minimize(
cost_func_vqe, x0, args=(raw_ansatz, observable_1, estimator), method="COBYLA"
)
end_time = time.time()
execution_time = end_time - start_time
print("INITIAL POINT 2:")
print(f"""Number of iterations: {x0_2_result.nfev}""")
print(f"""Time (s): {execution_time}""")
INITIAL POINT 2:
Number of iterations: 107
Time (s): 0.40889453887939453
Sa pamamagitan ng pag-eeksperimento sa iba't ibang paunang punto, maaari kang makamit ng mas mabilis na convergence at mas kaunting pagsusuri ng function.
Pag-eeksperimento sa Iba't Ibang Optimizerβ
Maaari nating ayusin ang optimizer gamit ang method argument ng SciPy minimize, na may higit pang mga opsyon na mahahanap dito. Ginamit natin noong una ang isang constrained minimizer (COBYLA). Sa halimbawang ito, tatalakayin natin ang paggamit ng unconstrained minimizer (BFGS) sa halip.
import time
start_time = time.time()
result = minimize(
cost_func_vqe, x0, args=(raw_ansatz, observable_1, estimator), method="BFGS"
)
end_time = time.time()
execution_time = end_time - start_time
print("CHANGED TO BFGS OPTIMIZER:")
print(f"""Number of iterations: {result.nfev}""")
print(f"""Time (s): {execution_time}""")
CHANGED TO BFGS OPTIMIZER:
Number of iterations: 117
Time (s): 0.31656408309936523
Halimbawa ng VQDβ
Dito isinasagawa natin ang Qiskit patterns framework nang hayagan.
Hakbang 1: I-map ang mga klasikal na input sa isang quantum na problemaβ
Sa halip na hanapin lamang ang pinakamababang eigenvalue ng ating mga observable, hahanapin natin ang lahat ng , (kung saan ang ).
Alalahanin na ang mga cost function ng VQD ay:
Mahalaga ito lalo na dahil kailangang ipasa ang isang vector na (sa kasong ito, ) bilang argumento kapag tinukoy natin ang VQD na object.
Gayundin, sa implementasyon ng VQD sa Qiskit, sa halip na isaalang-alang ang mga epektibong observable na inilarawan sa nakaraang notebook, ang mga fidelity na ay kinakalkula nang direkta sa pamamagitan ng ComputeUncompute na algorithm, na gumagamit ng Sampler primitive para i-sample ang probabilidad ng pagkuha ng para sa circuit na
. Gumagana ito nang tama dahil ang probabilidad na ito ay
ansatz = n_local(
num_qubits=2,
rotation_blocks=["ry", "rz"],
entanglement_blocks="cz",
# entanglement="linear",
reps=1,
)
ansatz.decompose().draw("mpl")

Magsimula tayo sa pagsusuri ng sumusunod na observable:
Ang observable na ito ay may mga sumusunod na eigenvalue:
At mga eigenstate:
from qiskit.quantum_info import SparsePauliOp
observable_2 = SparsePauliOp.from_list([("II", 2), ("XX", -3), ("YY", 2), ("ZZ", -4)])
Gagamitin natin ang sumusunod na function para kalkulahin ang overlap penalty. Tandaan na ito ay bahagi pa rin ng pag-map ng problema sa mga quantum circuit. Gayunpaman, tulad ng tinalakay sa nakaraang aralin, kinakalkula ng function na ito ang overlap sa pagitan ng kasalukuyang variational circuit at ng na-optimize na circuit mula sa isang nakaraang estado na may mas mababang energy/cost. Ang bagong circuit na nilikha ay kailangang i-transpile din para mapatakbo sa tunay na hardware. Nakita na natin ang function na ito noon, ginamit sa isang simulator. Dito, kailangan na nating isaalang-alang ang transpiling at kaugnay na optimization para sa paggamit ng tunay na backend, kaya naman may mga linya tungkol sa if realbackend == 1. Bahagyang hinahalo nito ang hakbang 2, ngunit tahasang banggitin natin ang hakbang 2 sa ibang pagkakataon.
import numpy as np
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
def calculate_overlaps(
ansatz, prev_circuits, parameters, sampler, realbackend, backend
):
def create_fidelity_circuit(circuit_1, circuit_2):
if len(circuit_1.clbits) > 0:
circuit_1.remove_final_measurements()
if len(circuit_2.clbits) > 0:
circuit_2.remove_final_measurements()
circuit = circuit_1.compose(circuit_2.inverse())
circuit.measure_all()
return circuit
overlaps = []
for prev_circuit in prev_circuits:
fidelity_circuit = create_fidelity_circuit(ansatz, prev_circuit)
if realbackend == 1:
pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
fidelity_circuit = pm.run(fidelity_circuit)
sampler_job = sampler.run([(fidelity_circuit, parameters)])
meas_data = sampler_job.result()[0].data.meas
counts_0 = meas_data.get_int_counts().get(0, 0)
shots = meas_data.num_shots
overlap = counts_0 / shots
overlaps.append(overlap)
return np.array(overlaps)
Idadagdag natin ngayon ang cost function ng VQD. Pansinin na kumpara sa nakaraang aralin, mayroon na tayong dalawang karagdagang argumento (realbackend at backend) para tulungan tayo sa transpilation kapag gumagamit ng tunay na mga backend.
def cost_func_vqd(
parameters,
ansatz,
prev_states,
step,
betas,
estimator,
sampler,
hamiltonian,
realbackend,
backend,
):
estimator_job = estimator.run([(ansatz, hamiltonian, [parameters])])
total_cost = 0
if step > 1:
overlaps = calculate_overlaps(
ansatz, prev_states, parameters, sampler, realbackend, backend
)
total_cost = np.sum(
[np.real(betas[state] * overlap) for state, overlap in enumerate(overlaps)]
)
estimator_result = estimator_job.result()[0]
value = estimator_result.data.evs[0] + total_cost
return value
Muli, gagamitin natin ang mga simulator para sa pag-debug, at pagkatapos ay lilipat sa tunay na hardware.
from qiskit.primitives import StatevectorSampler
from qiskit.primitives import StatevectorEstimator
sampler = StatevectorSampler(default_shots=4092)
estimator = StatevectorEstimator()
Dito ipinakilala natin ang bilang ng mga estado na nais nating kalkulahin, ang mga penalty, at isang hanay ng mga paunang parameter, x0.
from qiskit.quantum_info import SparsePauliOp
k = 4
betas = [50, 60, 40]
x0 = np.ones(8)
Susubukan na natin ang algorithm gamit ang mga simulator:
from scipy.optimize import minimize
prev_states = []
prev_opt_parameters = []
eigenvalues = []
realbackend = 0
for step in range(1, k + 1):
if step > 1:
prev_states.append(ansatz.assign_parameters(prev_opt_parameters))
result = minimize(
cost_func_vqd,
x0,
args=(
ansatz,
prev_states,
step,
betas,
estimator,
sampler,
observable_2,
realbackend,
None,
),
method="COBYLA",
options={"maxiter": 200, "tol": 0.000001},
)
print(result)
prev_opt_parameters = result.x
eigenvalues.append(result.fun)
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: -6.9999999999996
x: [ 1.571e+00 1.571e+00 2.519e+00 2.100e+00 1.242e+00
6.935e-01 2.298e+00 1.991e+00]
nfev: 151
maxcv: 0.0
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: 3.698974255258432
x: [ 1.269e+00 1.109e+00 1.080e+00 1.200e+00 1.094e+00
1.163e+00 9.752e-01 9.519e-01]
nfev: 103
maxcv: 0.0
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: 4.731320121938101
x: [ 1.533e+00 2.451e+00 2.526e+00 2.406e+00 1.968e+00
2.105e+00 8.537e-01 8.442e-01]
nfev: 110
maxcv: 0.0
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: 7.008239313655201
x: [ 4.150e+00 2.120e+00 3.495e+00 7.262e-01 1.953e+00
-1.982e-01 3.263e-01 2.563e+00]
nfev: 126
maxcv: 0.0
eigenvalues
[np.float64(-6.9999999999996),
np.float64(3.698974255258432),
np.float64(4.731320121938101),
np.float64(7.008239313655201)]
Ang mga resultang ito ay medyo malapit sa mga inaasahang halaga, maliban sa approximation error at global phase. Maaari nating ayusin ang tolerance ng klasikal na optimizer at/o ang mga penalty para sa statevector overlap para makakuha ng mas tumpak na mga halaga.
solution_eigenvalues = [-7, 3, 5, 7]
for index, experimental_eigenvalue in enumerate(eigenvalues):
solution_eigenvalue = solution_eigenvalues[index]
print(
f"Percent error: {abs((experimental_eigenvalue - solution_eigenvalue)/solution_eigenvalue):.2e}"
)
Percent error: 5.71e-14
Percent error: 2.33e-01
Percent error: 5.37e-02
Percent error: 1.18e-03
Baguhin ang mga betaβ
Tulad ng nabanggit sa nakaraang aralin, ang mga halaga ng ay dapat na mas malaki kaysa sa pagkakaiba ng mga eigenvalue. Tingnan natin kung ano ang mangyayari kapag hindi nila natutugunan ang kundisyong iyon gamit ang
na may mga eigenvalue na
from qiskit.quantum_info import SparsePauliOp
k = 4
betas = np.ones(3)
x0 = np.zeros(8)
from scipy.optimize import minimize
prev_states = []
prev_opt_parameters = []
eigenvalues = []
realbackend = 0
for step in range(1, k + 1):
if step > 1:
prev_states.append(ansatz.assign_parameters(prev_opt_parameters))
result = minimize(
cost_func_vqd,
x0,
args=(
ansatz,
prev_states,
step,
betas,
estimator,
sampler,
observable_2,
realbackend,
None,
),
method="COBYLA",
options={"tol": 0.01, "maxiter": 200},
)
print(result)
prev_opt_parameters = result.x
eigenvalues.append(result.fun)
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: -6.999916534745094
x: [ 1.568e+00 -1.569e+00 1.385e-01 1.398e-01 -7.972e-01
7.835e-01 -2.375e-01 4.539e-02]
nfev: 125
maxcv: 0.0
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: -1.515139929812874
x: [-5.317e-04 -2.514e-03 1.016e+00 9.998e-01 3.890e-04
1.772e-04 1.568e-04 8.497e-04]
nfev: 35
maxcv: 0.0
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: -0.509948114293115
x: [-3.796e-03 8.853e-03 3.015e-04 9.997e-01 6.271e-04
-2.554e-03 1.017e-04 2.766e-04]
nfev: 37
maxcv: 0.0
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: 0.4914672235935682
x: [-7.178e-03 -8.652e-03 1.125e+00 -5.428e-02 -1.586e-03
2.031e-03 -3.462e-03 5.734e-03]
nfev: 35
maxcv: 0.0
solution_eigenvalues = [-7, 3, 5, 7]
for index, experimental_eigenvalue in enumerate(eigenvalues):
solution_eigenvalue = solution_eigenvalues[index]
print(
f"Percent error: {abs((experimental_eigenvalue - solution_eigenvalue)/solution_eigenvalue):.2e}"
)
Percent error: 1.19e-05
Percent error: 1.51e+00
Percent error: 1.10e+00
Percent error: 9.30e-01
Sa pagkakataong ito, ibinabalik ng optimizer ang parehong estado na bilang iminungkahing solusyon sa lahat ng eigenstate β na malinaw na mali. Nangyayari ito dahil masyadong maliit ang mga beta para parusahan ang minimum eigenstate sa mga sunod na cost function. Kaya naman, hindi ito nai-exclude mula sa epektibong search space sa mga susunod na iteration ng algorithm, at palaging pinipili bilang pinakamahusay na posibleng solusyon.
Inirerekomenda namin na mag-eksperimento sa mga halaga ng , at tiyaking mas malaki ang mga ito kaysa sa pagkakaiba ng mga eigenvalue.
Hakbang 2: I-optimize ang problema para sa quantum executionβ
Para mapatakbo ito sa tunay na hardware, kailangan nating i-optimize ang mga quantum circuit para sa quantum computer na ating pipiliin. Para sa ating layunin dito, gagamitin natin ang pinaka-hindi abala na backend.
from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit_ibm_runtime import Session, EstimatorOptions
from qiskit_ibm_runtime import QiskitRuntimeService
service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
# Or use a specific backend
# backend = service.backend("ibm_brisbane")
print(backend)
<IBMBackend('ibm_brisbane')>
Ita-transpile natin ang ating circuit gamit ang isang preset pass manager at optimization level 3.
pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_ansatz = pm.run(ansatz)
isa_observable = observable_2.apply_layout(layout=isa_ansatz.layout)
Hakbang 3: I-execute gamit ang mga Qiskit primitiveβ
Siguraduhing naibalik ang mga beta sa sapat na mataas na halaga, maaari na nating patakbuhin ang ating kalkulasyon sa tunay na quantum hardware.
# Estimated compute resource usage: 25 minutes. Benchmarked at 24 min, 30 sec on an Eagle r3 processor on 5-30-24
k = 2
betas = [30, 50, 80]
x0 = np.zeros(8)
real_prev_states = []
real_prev_opt_parameters = []
real_eigenvalues = []
realbackend = 1
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)
with Session(backend=backend) as session:
estimator = Estimator(mode=session, options=estimator_options)
sampler = Sampler(mode=session)
for step in range(1, k + 1):
if step > 1:
real_prev_states.append(isa_ansatz.assign_parameters(prev_opt_parameters))
result = minimize(
cost_func_vqd,
x0,
args=(
isa_ansatz,
real_prev_states,
step,
betas,
estimator,
sampler,
isa_observable,
realbackend,
backend,
),
method="COBYLA",
options={"maxiter": 200},
)
print(result)
real_prev_opt_parameters = result.x
real_eigenvalues.append(result.fun)
session.close()
print(real_eigenvalues)
Hakbang 4: Post-process, ibalik ang resulta sa klasikal na formatβ
Ang ating output ay structurally katulad ng tinalakay sa mga nakaraang aralin at halimbawa. Ngunit may isang problemadong bagay sa mga resulta sa itaas, mula sa kung saan maaari tayong makakuha ng isang babala para sa konteksto ng mga excited state. Para malimitahan ang oras ng computing na ginagamit sa learning example na ito, nagtakda tayo ng maximum na bilang ng mga iteration para sa klasikal na optimizer na posibleng masyadong mababa: 200 iterations. Ang isang nakaraang kalkulasyon sa itaas, sa isang simulator, ay nabigong mag-converge sa loob ng 200 iterations. Dito, nag-converge ang atin... ngunit sa anong tolerance? Hindi natin tinukoy ang isang tolerance para ikonsidera ng COBYLA ang sarili nitong "converged". Isang tingin sa halaga ng function at paghahambing sa mga nakaraang run ay nagsasabi sa atin na hindi pa malapit ang COBYLA sa pag-converge sa katumpakan na kailangan natin.
May isa pang isyu: ang energy ng unang excited state ay tila mas mababa kaysa sa energy ng ground state! Subukan mong ipaliwanag kung paano ito maaaring mangyari. Pahiwatig: ito ay may kaugnayan sa convergence point na ating tinalakay lamang. Ang ganitong gawi ay inipaliwanag nang detalyado sa ibaba pagkatapos mailapat ang VQD sa H2 molecule.
Quantum chemistry: ground state at excited energy solverβ
Ang layunin natin ay i-minimize ang expectation value ng observable na kumakatawan sa enerhiya (Hamiltonian ):
from qiskit.quantum_info import SparsePauliOp
from qiskit.circuit.library import efficient_su2
H2_op = SparsePauliOp.from_list(
[
("II", -1.052373245772859),
("IZ", 0.39793742484318045),
("ZI", -0.39793742484318045),
("ZZ", -0.01128010425623538),
("XX", 0.18093119978423156),
]
)
chem_ansatz = efficient_su2(H2_op.num_qubits)
chem_ansatz.decompose().draw("mpl")

from qiskit import QuantumCircuit
def cost_func_vqe(params, ansatz, hamiltonian, estimator):
"""Return estimate of energy from estimator
Parameters:
params (ndarray): Array of ansatz parameters
ansatz (QuantumCircuit): Parameterized ansatz circuit
hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
estimator (Estimator): Estimator primitive instance
Returns:
float: Energy estimate
"""
pub = (ansatz, hamiltonian, params)
cost = estimator.run([pub]).result()[0].data.evs
# cost = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
return cost
Magtatakda tayo ngayon ng paunang set ng mga parameter:
import numpy as np
x0 = np.ones(chem_ansatz.num_parameters)
Maaari nating i-minimize ang cost function na ito para makuha ang pinakamainam na mga parameter, at maaari nating subukan muna ang ating code gamit ang isang lokal na simulator.
from qiskit.primitives import StatevectorEstimator as Estimator
from qiskit.primitives import StatevectorSampler as Sampler
estimator = Estimator()
sampler = Sampler()
# SciPy minimizer routine
from scipy.optimize import minimize
import time
start_time = time.time()
result = minimize(
cost_func_vqe, x0, args=(chem_ansatz, H2_op, estimator), method="COBYLA"
)
end_time = time.time()
execution_time = end_time - start_time
result
message: Optimization terminated successfully.
success: True
status: 1
fun: -1.857275029048451
x: [ 7.326e-01 1.354e+00 ... 1.040e+00 1.508e+00]
nfev: 242
maxcv: 0.0
Ang pinakamababang halaga ng cost function (-1.857...) ay ang ground state energy ng molekula ng H2, sa yunit ng hartrees.
Mga Excited Stateβ
Maaari rin nating gamitin ang VQD para malutas ang na kabuuang estado (ang ground state at ang unang excited state).
from qiskit.quantum_info import SparsePauliOp
import numpy as np
k = 2
betas = [33, 33]
# x0 = np.zeros(ansatz.num_parameters)
x0 = [
1.164e00,
-2.438e-01,
9.358e-04,
6.745e-02,
1.990e00,
9.810e-02,
6.154e-01,
5.454e-01,
]
Idadagdag natin ang ating overlap na kalkulasyon:
from scipy.optimize import minimize
prev_states = []
prev_opt_parameters = []
eigenvalues = []
realbackend = 0
for step in range(1, k + 1):
if step > 1:
prev_states.append(ansatz.assign_parameters(prev_opt_parameters))
result = minimize(
cost_func_vqd,
x0,
args=(
ansatz,
prev_states,
step,
betas,
estimator,
sampler,
H2_op,
realbackend,
None,
),
method="COBYLA",
options={"tol": 0.001, "maxiter": 2000},
)
print(result)
prev_opt_parameters = result.x
eigenvalues.append(result.fun)
message: Optimization terminated successfully.
success: True
status: 1
fun: -1.8572671093941977
x: [ 1.164e+00 -2.437e-01 2.118e-03 6.448e-02 1.990e+00
9.870e-02 6.167e-01 5.476e-01]
nfev: 58
maxcv: 0.0
message: Optimization terminated successfully.
success: True
status: 1
fun: -1.0322873777662176
x: [ 3.205e+00 1.502e+00 1.699e+00 -1.107e-02 3.086e+00
1.530e+00 4.445e-02 7.013e-02]
nfev: 99
maxcv: 0.0
eigenvalues
[-1.8572671093941977, -1.0322873777662176]
Tunay na hardware at isang pangwakas na babalaβ
Para patakbuhin ito sa tunay na hardware, kailangan nating i-optimize ang mga quantum circuit para sa ating piling quantum computer. Para sa ating layunin dito, gagamitin natin ang pinaka-hindi abala na backend.
from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit_ibm_runtime import Session, EstimatorOptions
from qiskit_ibm_runtime import QiskitRuntimeService
service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
Gagamit tayo ng preset pass manager para sa transpilation, at ino-optimize natin nang husto ang ating circuit gamit ang optimization level 3.
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_ansatz = pm.run(ansatz)
isa_observable = H2_op.apply_layout(layout=isa_ansatz.layout)
Dahil ang VQD ay lubhang iteratibo, isasagawa natin ang lahat ng hakbang sa loob ng isang Runtime session, upang ang ating mga trabaho ay ma-queue lamang sa simula, at hindi sa pagitan ng bawat pag-update ng parameter. Wala nang ibang nagbabago sa syntax ng cost function o estimator.
x0 = [
1.306e00,
-2.284e-01,
6.913e-02,
-2.530e-02,
1.849e00,
7.433e-02,
6.366e-01,
5.600e-01,
]
# Estimated hardware usage: 20 min benchmarked on an Eagle r3 processor on 5-30-24
real_prev_states = []
real_prev_opt_parameters = []
real_eigenvalues = []
realbackend = 1
estimator_options = EstimatorOptions(resilience_level=1, default_shots=4096)
with Session(backend=backend) as session:
estimator = Estimator(mode=session)
sampler = Sampler(mode=session)
for step in range(1, k + 1):
if step > 1:
real_prev_states.append(
isa_ansatz.assign_parameters(real_prev_opt_parameters)
)
result = minimize(
cost_func_vqd,
x0,
args=(
isa_ansatz,
real_prev_states,
step,
betas,
estimator,
sampler,
isa_observable,
realbackend,
backend,
),
method="COBYLA",
options={"tol": 0.001, "maxiter": 300},
)
print(result)
real_prev_opt_parameters = result.x
real_eigenvalues.append(result.fun)
session.close()
print(real_eigenvalues)
Ang nakuhang ground state energy (-1.83 hartrees) ay hindi masyadong malayo sa tamang halaga (-1.85 hartrees). Gayunpaman, ang excited state energy ay medyo malayo sa tamang sagot. Katulad ito ng maling gawi na nakita natin kanina sa araling ito. Ang energy na iniulat para sa excited state ay halos kapareho ng sa ground state. Sa nakaraang kaso, nakita pa natin ang isang excited state energy na mas mababa kaysa sa iniulat na ground state energy.
Hindi posible para sa isang variational na kalkulasyon na makakuha ng energy na mas mababa kaysa sa tunay na ground state energy. Sa nakaraang pagkakataon, ang ground state energy na ating nakuha ay hindi masyadong malapit sa tunay na ground state. Dahil hindi natin nakuha ang tunay na ground state energy sa kasong iyon, walang kontradiksyon. Sa kasalukuyang kaso, ang ground state energy ay medyo malapit sa tamang halaga, at gayunpaman ang excited state energy ay tila kakatuwang malapit sa parehong halaga.
Para mas maunawaan kung paano ito nangyari, alalahanin na ang paraan ng paghahanap natin ng excited state ay sa pamamagitan ng pag-require na ang variational state ay ortogonal sa ground state (gamit ang mga overlap circuit at penalty term). Kung hindi tayo nakakakuha ng tumpak na ground state energy (o nangangahulugang may ilang porsyentong pagkakamali), hindi rin tayo nakakakuha ng tumpak na ground state vector! Kaya nang i-require natin na ang excited state ay ortogonal sa unang estado na ating nahanap, hindi tayo nagpapataw ng ortogonalidad sa tunay na ground state, kundi sa isang approximation nito (minsan isang mahinang approximation). Kaya naman, ang excited state ay hindi napilitang maging ortogonal sa tunay na ground state, at ang ating mga energy estimate para sa mga excited state ay talagang malapit sa ground state energy.
Ito ay palaging magiging alalahanin sa VQD. Ngunit sa prinsipyo, maaari itong itama sa pamamagitan ng pagdaragdag ng maximum na bilang ng mga iteration para sa klasikal na optimizer, pagpapataw ng mas mababang tolerance para sa klasikal na optimizer, at posibleng subukan din ang ibang ansatz kung lagi tayong hindi nakakatama sa tunay na ground state. Tulad ng ating nakita, maaaring kailangan ding baguhin ang mga overlap penalty (betas). Ngunit iyon ay isang hiwalay na isyu talaga. Walang penalty para sa overlap ang makakapagsiguro na mananatili kang malayo sa tunay na ground state, kung hindi ka pa nakakakita ng napakagandang estimate ng tunay na ground state para sa overlap circuit.
Optimization: Max-Cutβ
Ang maximum cut (Max-Cut) na problema ay isang combinatorial optimization na problema na kinabibilangan ng paghahati ng mga vertice ng isang graph sa dalawang magkahiwalay na set upang ma-maximize ang bilang ng mga gilid sa pagitan ng dalawang set. Sa mas pormal na paraan, sa isang undirected graph , kung saan ang set ng mga vertice at ang set ng mga gilid, tinatanong ng Max-Cut na problema kung paano hahatiin ang mga vertice sa dalawang magkahiwalay na subset, at , upang ma-maximize ang bilang ng mga gilid na may isang dulo sa at ang isa pa sa .
Maaari nating ilapat ang Max-Cut para malutas ang iba't ibang problema, tulad ng clustering, network design, at phase transition. Magsisimula tayo sa pamamagitan ng paggawa ng problem graph:
import rustworkx as rx
from rustworkx.visualization import mpl_draw
n = 4
G = rx.PyGraph()
G.add_nodes_from(range(n))
# The edge syntax is (start, end, weight)
edges = [(0, 1, 1.0), (0, 2, 1.0), (0, 3, 1.0), (1, 2, 1.0), (2, 3, 1.0)]
G.add_edges_from(edges)
mpl_draw(
G, pos=rx.shell_layout(G), with_labels=True, edge_labels=str, node_color="#1192E8"
)
Maaaring ipahayag ang problemang ito bilang isang binary optimization na problema. Para sa bawat node , kung saan ang bilang ng mga node ng graph (sa kasong ito ), isasaalang-alang natin ang binary variable . Ang variable na ito ay magkakaroon ng halaga kung ang node ay kabilang sa grupong tatawaging at kung ito ay nasa kabilang grupo, na tatawaging . Itatanda rin natin bilang (elemento ng adjacency matrix ) ang bigat ng gilid na pumupunta mula sa node patungong node . Dahil ang graph ay undirected, . Maaari nating buuhin ang ating problema bilang pag-maximize ng sumusunod na cost function:
Para malutas ang problemang ito gamit ang isang quantum computer, ipapahayag natin ang cost function bilang expected value ng isang observable. Gayunpaman, ang mga observable na tinatanggap ng Qiskit nang likas ay binubuo ng mga Pauli operator, na may eigenvalue na at sa halip na at . Kaya naman gagawa tayo ng sumusunod na pagbabago ng variable:
Kung saan . Maaari nating gamitin ang adjacency matrix para madaling ma-access ang mga bigat ng lahat ng gilid. Gagamitin ito para makuha ang ating cost function:
Nangangahulugan ito na:
Kaya ang bagong cost function na gusto nating i-maximize ay:
Bukod pa rito, ang likas na ugali ng isang quantum computer ay maghanap ng minima (kadalasan ang pinakamababang enerhiya) sa halip na maxima, kaya sa halip na i-maximize ang , i-minimize na lang natin ang:
Ngayong mayroon na tayong cost function na i-minimize na ang mga variable ay maaaring magkaroon ng halaga na at , maaari tayong gumawa ng sumusunod na analohiya sa Pauli :
Sa ibang salita, ang variable ay katumbas ng isang gate na kumikilos sa qubit . Bukod pa rito:
Kaya ang observable na ating isasaalang-alang ay:
na kung saan kailangan pa nating idagdag ang independyenteng termino sa bandang huli:
Ang operator ay isang linear na kombinasyon ng mga termino na may Z operator sa mga node na konektado ng isang gilid (alalahanin na ang 0th qubit ay nasa pinakakanang bahagi): . Kapag naitayo na ang operator, ang ansatz para sa QAOA algorithm ay madaling mabubuo sa pamamagitan ng paggamit ng QAOAAnsatz circuit mula sa Qiskit circuit library.
from qiskit.circuit.library import QAOAAnsatz
from qiskit.quantum_info import SparsePauliOp
max_hamiltonian = SparsePauliOp.from_list(
[("IIZZ", 1), ("IZIZ", 1), ("IZZI", 1), ("ZIIZ", 1), ("ZZII", 1)]
)
max_ansatz = QAOAAnsatz(max_hamiltonian, reps=2)
# Draw
max_ansatz.decompose(reps=3).draw("mpl")
# Sum the weights, and divide by 2
offset = -sum(edge[2] for edge in edges) / 2
print(f"""Offset: {offset}""")
Offset: -2.5
def cost_func(params, ansatz, hamiltonian, estimator):
"""Return estimate of energy from estimator
Parameters:
params (ndarray): Array of ansatz parameters
ansatz (QuantumCircuit): Parameterized ansatz circuit
hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
estimator (Estimator): Estimator primitive instance
Returns:
float: Energy estimate
"""
pub = (ansatz, hamiltonian, params)
cost = estimator.run([pub]).result()[0].data.evs
# cost = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
return cost
from qiskit.primitives import StatevectorEstimator as Estimator
from qiskit.primitives import StatevectorSampler as Sampler
estimator = Estimator()
sampler = Sampler()
Magtatakda tayo ngayon ng paunang set ng mga random na parameter:
import numpy as np
x0 = 2 * np.pi * np.random.rand(max_ansatz.num_parameters)
print(x0)
[6.0252949 0.58448176 2.15785731 1.13646074]
Maaaring gamitin ang anumang klasikal na optimizer para i-minimize ang cost function. Sa isang tunay na quantum system, ang isang optimizer na dinisenyo para sa mga hindi maayos na cost function landscape ay karaniwang mas epektibo. Dito ginagamit natin ang COBYLA routine mula sa SciPy sa pamamagitan ng minimize function.
Dahil paulit-ulit tayong nagsasagawa ng maraming tawag sa Runtime, gumagamit tayo ng session para isagawa ang lahat ng tawag sa loob ng isang bloke. Bukod pa rito, para sa QAOA, ang solusyon ay nakalagak sa output distribution ng ansatz circuit na nakatali sa pinakamainam na mga parameter mula sa minimization. Kaya naman kailangan natin ng Sampler primitive, at ibibigay natin ito sa parehong session. At patakbuhin ang ating minimization routine:
result = minimize(
cost_func, x0, args=(max_ansatz, max_hamiltonian, estimator), method="COBYLA"
)
print(result)
message: Optimization terminated successfully.
success: True
status: 1
fun: -2.585287311689236
x: [ 7.332e+00 3.904e-01 2.045e+00 1.028e+00]
nfev: 80
maxcv: 0.0
Ang solution vector ng mga parameter angle (x), kapag ipinasok sa ansatz circuit, ay nagbubunga ng graph partitioning na ating hinahanap.
eigenvalue = cost_func(result.x, max_ansatz, max_hamiltonian, estimator)
print(f"""Eigenvalue: {eigenvalue}""")
print(f"""Max-Cut Objective: {eigenvalue + offset}""")
Eigenvalue: -2.585287311689236
Max-Cut Objective: -5.085287311689235
from qiskit.result import QuasiDistribution
from qiskit.primitives import StatevectorSampler
sampler = StatevectorSampler()
# Assign solution parameters to ansatz
qc = max_ansatz.assign_parameters(result.x)
# Add measurements to our circuit
qc.measure_all()
# Sample ansatz at optimal parameters
# samp_dist = sampler.run(qc).result().quasi_dists[0]
shots = 1024
job = sampler.run([qc], shots=shots)
qc.decompose().draw("mpl")
data_pub = job.result()[0].data
bitstrings = data_pub.meas.get_bitstrings()
counts = data_pub.meas.get_counts()
quasi_dist = QuasiDistribution(
{outcome: freq / shots for outcome, freq in counts.items()}
)
probabilities = quasi_dist
# Close the session since we are now done with it
# session.close()
from qiskit.visualization import plot_distribution
plot_distribution(counts)
binary_string = max(counts.items(), key=lambda kv: kv[1])[0]
x = np.asarray([int(y) for y in reversed(list(binary_string))])
colors = ["r" if x[i] == 0 else "c" for i in range(n)]
mpl_draw(
G, pos=rx.shell_layout(G), with_labels=True, edge_labels=str, node_color=colors
)
Buodβ
Sa araling ito, natutunan mo ang:
- Kung paano magsulat ng custom na variational algorithm
- Kung paano ilapat ang isang variational algorithm para mahanap ang mga pinakamababang eigenvalue
- Kung paano gamitin ang mga variational algorithm para malutas ang mga praktikal na kaso ng aplikasyon