Lumaktaw sa pangunahing nilalaman

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:

O^1=2IIβˆ’2XX+3YYβˆ’3ZZ,\hat{O}_1 = 2 II - 2 XX + 3 YY - 3 ZZ,

Ang observable na ito ay may mga sumusunod na eigenvalue:

{Ξ»0=βˆ’6Ξ»1=4Ξ»2=4Ξ»3=6}\left\{ \begin{array}{c} \lambda_0 = -6 \\ \lambda_1 = 4 \\ \lambda_2 = 4 \\ \lambda_3 = 6 \end{array} \right\}

At mga eigenstate:

{βˆ£Ο•0⟩=12(∣00⟩+∣11⟩)βˆ£Ο•1⟩=12(∣00βŸ©βˆ’βˆ£11⟩)βˆ£Ο•2⟩=12(∣01βŸ©βˆ’βˆ£10⟩)βˆ£Ο•3⟩=12(∣01⟩+∣10⟩)}\left\{ \begin{array}{c} |\phi_0\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)\\ |\phi_1\rangle = \frac{1}{\sqrt{2}}(|00\rangle - |11\rangle)\\ |\phi_2\rangle = \frac{1}{\sqrt{2}}(|01\rangle - |10\rangle)\\ |\phi_3\rangle = \frac{1}{\sqrt{2}}(|01\rangle + |10\rangle) \end{array} \right\}
# 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 O^1\hat{O}_1. 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")

Output of the previous code cell

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 URU_R. Ngayon, pag-isipan natin kung paano makukuha ang ideal eigenstate na 12(∣00⟩+∣11⟩)\frac{1}{\sqrt{2}}(|00\rangle + |11\rangle). 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")

Output of the previous code cell

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")

Output of the previous code cell

Para sa bagong circuit na ito, ang ideal na solusyon ay maaabot kapag lahat ng parameter ay nakatakda sa 00, 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 ΞΈ0βƒ—\vec{\theta_0}. Sa partikular, gagamitin natin ang ΞΈ0βƒ—=(0,0,0,0,6,0,0,0)\vec{\theta_0}=(0,0,0,0,6,0,0,0) at ΞΈ0βƒ—=(6,6,6,6,6,6,6,6,6)\vec{\theta_0}=(6,6,6,6,6,6,6,6,6).

Tandaan na, gaya ng tinalakay nang ipinakilala ang reference state, ang ideal na solusyon ay mahahanap kapag lahat ng parameter ay 00, 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 ΞΈ0βƒ—=(6,6,6,6,6,6,6,6,6)\vec{\theta_0}=(6,6,6,6,6,6,6,6,6):

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 44, (kung saan ang k=4k=4).

Alalahanin na ang mga cost function ng VQD ay:

Cl(ΞΈβƒ—):=⟨ψ(ΞΈβƒ—)∣H^∣ψ(ΞΈβƒ—)⟩+βˆ‘j=0lβˆ’1Ξ²j∣⟨ψ(ΞΈβƒ—)∣ψ(ΞΈjβƒ—)⟩∣2βˆ€l∈{1,⋯ ,k}={1,⋯ ,4}C_{l}(\vec{\theta}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle + \sum_{j=0}^{l-1}\beta_j |\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2 \quad \forall l\in\{1,\cdots,k\}=\{1,\cdots,4\}

Mahalaga ito lalo na dahil kailangang ipasa ang isang vector na Ξ²βƒ—=(Ξ²0,⋯ ,Ξ²kβˆ’1)\vec{\beta}=(\beta_0,\cdots,\beta_{k-1}) (sa kasong ito, (Ξ²0,Ξ²1,Ξ²2,Ξ²3)(\beta_0, \beta_1, \beta_2, \beta_3)) 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 ∣⟨ψ(ΞΈβƒ—)∣ψ(ΞΈjβƒ—)⟩∣2|\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2 ay kinakalkula nang direkta sa pamamagitan ng ComputeUncompute na algorithm, na gumagamit ng Sampler primitive para i-sample ang probabilidad ng pagkuha ng ∣0⟩|0\rangle para sa circuit na UA†(ΞΈβƒ—)UA(ΞΈjβƒ—)U_A^\dagger(\vec{\theta})U_A(\vec{\theta^j}). Gumagana ito nang tama dahil ang probabilidad na ito ay

p0=∣⟨0∣UA†(ΞΈβƒ—)UA(ΞΈjβƒ—)∣0⟩∣2=∣(⟨0∣UA†(ΞΈβƒ—))(UA(ΞΈjβƒ—)∣0⟩)∣2=∣⟨ψ(ΞΈβƒ—)∣ψ(ΞΈjβƒ—)⟩∣2\begin{aligned} p_0 & = |\langle 0|U_A^\dagger(\vec{\theta})U_A(\vec{\theta^j})|0\rangle|^2 \\[1mm] & = |\big(\langle 0|U_A^\dagger(\vec{\theta})\big)\big(U_A(\vec{\theta^j})|0\rangle\big)|^2 \\[1mm] & = |\langle \psi(\vec{\theta}) |\psi(\vec{\theta^j}) \rangle|^2 \\[1mm] \end{aligned}
ansatz = n_local(
num_qubits=2,
rotation_blocks=["ry", "rz"],
entanglement_blocks="cz",
# entanglement="linear",
reps=1,
)

ansatz.decompose().draw("mpl")

Output of the previous code cell

Magsimula tayo sa pagsusuri ng sumusunod na observable:

O^2:=2IIβˆ’3XX+2YYβˆ’4ZZ\hat{O}_2 := 2 II - 3 XX + 2 YY - 4 ZZ

Ang observable na ito ay may mga sumusunod na eigenvalue:

{Ξ»0=βˆ’7Ξ»1=3Ξ»2=5Ξ»3=7}\left\{ \begin{array}{c} \lambda_0 = -7 \\ \lambda_1 = 3\\ \lambda_2 = 5 \\ \lambda_3 = 7 \end{array} \right\}

At mga eigenstate:

{βˆ£Ο•0⟩=12(∣00⟩+∣11⟩)βˆ£Ο•1⟩=12(∣00βŸ©βˆ’βˆ£11⟩)βˆ£Ο•2⟩=12(∣01⟩+∣10⟩)βˆ£Ο•3⟩=12(∣01βŸ©βˆ’βˆ£10⟩)}\left\{ \begin{array}{c} |\phi_0\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)\\ |\phi_1\rangle = \frac{1}{\sqrt{2}}(|00\rangle - |11\rangle)\\ |\phi_2\rangle = \frac{1}{\sqrt{2}}(|01\rangle + |10\rangle)\\ |\phi_3\rangle = \frac{1}{\sqrt{2}}(|01\rangle - |10\rangle) \end{array} \right\}
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 Ξ²βƒ—\vec{\beta} 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 O^2\hat{O}_2

O^2=2IIβˆ’3XX+2YYβˆ’4ZZ\hat{O}_2 = 2 II - 3 XX + 2 YY - 4 ZZ

na may mga eigenvalue na

{Ξ»0=βˆ’7Ξ»1=3Ξ»2=5Ξ»3=7}\left\{ \begin{array}{c} \lambda_0 = -7 \\ \lambda_1 = 3\\ \lambda_2 = 5 \\ \lambda_3 = 7 \end{array} \right\}
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 βˆ£Ο•0⟩=12(∣00⟩+∣11⟩)|\phi_0\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle) 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 Ξ²βƒ—\vec{\beta}, 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 H^\hat{\mathcal{H}}):

minβ‘ΞΈβƒ—βŸ¨Οˆ(ΞΈβƒ—)∣H^∣ψ(ΞΈβƒ—)⟩\min_{\vec\theta} \langle\psi(\vec\theta)|\hat{\mathcal{H}}|\psi(\vec\theta)\rangle
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")

Output of the previous code cell

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 k=2k=2 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 G=(V,E)G=(V,E), kung saan VV ang set ng mga vertice at EE ang set ng mga gilid, tinatanong ng Max-Cut na problema kung paano hahatiin ang mga vertice sa dalawang magkahiwalay na subset, SS at TT, upang ma-maximize ang bilang ng mga gilid na may isang dulo sa SS at ang isa pa sa TT.

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"
)

Output of the previous code cell

Maaaring ipahayag ang problemang ito bilang isang binary optimization na problema. Para sa bawat node 0≀i<n0 \leq i < n, kung saan nn ang bilang ng mga node ng graph (sa kasong ito n=4n=4), isasaalang-alang natin ang binary variable xix_i. Ang variable na ito ay magkakaroon ng halaga 11 kung ang node ii ay kabilang sa grupong tatawaging 11 at 00 kung ito ay nasa kabilang grupo, na tatawaging 00. Itatanda rin natin bilang wijw_{ij} (elemento (i,j)(i,j) ng adjacency matrix ww) ang bigat ng gilid na pumupunta mula sa node ii patungong node jj. Dahil ang graph ay undirected, wij=wjiw_{ij}=w_{ji}. Maaari nating buuhin ang ating problema bilang pag-maximize ng sumusunod na cost function:

C(xβƒ—)=βˆ‘i,j=0nwijxi(1βˆ’xj)=βˆ‘i,j=0nwijxiβˆ’βˆ‘i,j=0nwijxixj=βˆ‘i,j=0nwijxiβˆ’βˆ‘i=0nβˆ‘j=0i2wijxixj\begin{aligned} C(\vec{x}) & =\sum_{i,j=0}^n w_{ij} x_i(1-x_j)\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i,j=0}^n w_{ij} x_ix_j\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i=0}^n \sum_{j=0}^i 2w_{ij} x_ix_j \end{aligned}

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 11 at βˆ’1-1 sa halip na 00 at 11. Kaya naman gagawa tayo ng sumusunod na pagbabago ng variable:

Kung saan xβƒ—=(x0,x1,⋯ ,xnβˆ’1)\vec{x}=(x_0,x_1,\cdots ,x_{n-1}). Maaari nating gamitin ang adjacency matrix ww para madaling ma-access ang mga bigat ng lahat ng gilid. Gagamitin ito para makuha ang ating cost function:

zi=1βˆ’2xiβ†’xi=1βˆ’zi2z_i = 1-2x_i \rightarrow x_i = \frac{1-z_i}{2}

Nangangahulugan ito na:

xi=0β†’zi=1xi=1β†’zi=βˆ’1.\begin{array}{lcl} x_i=0 & \rightarrow & z_i=1 \\ x_i=1 & \rightarrow & z_i=-1.\end{array}

Kaya ang bagong cost function na gusto nating i-maximize ay:

C(zβƒ—)=βˆ‘i,j=0nwij(1βˆ’zi2)(1βˆ’1βˆ’zj2)=βˆ‘i,j=0nwij4βˆ’βˆ‘i,j=0nwij4zizj=βˆ‘i=0nβˆ‘j=0iwij2βˆ’βˆ‘i=0nβˆ‘j=0iwij2zizj\begin{aligned} C(\vec{z}) & = \sum_{i,j=0}^n w_{ij} \bigg(\frac{1-z_i}{2}\bigg)\bigg(1-\frac{1-z_j}{2}\bigg)\\[1mm] & = \sum_{i,j=0}^n \frac{w_{ij}}{4} - \sum_{i,j=0}^n \frac{w_{ij}}{4} z_iz_j\\[1mm] & = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j \end{aligned}

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 C(z⃗)C(\vec{z}), i-minimize na lang natin ang:

βˆ’C(zβƒ—)=βˆ‘i=0nβˆ‘j=0iwij2zizjβˆ’βˆ‘i=0nβˆ‘j=0iwij2-C(\vec{z}) = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}

Ngayong mayroon na tayong cost function na i-minimize na ang mga variable ay maaaring magkaroon ng halaga na βˆ’1-1 at 11, maaari tayong gumawa ng sumusunod na analohiya sa Pauli ZZ:

zi≑Zi=I⏞nβˆ’1βŠ—...βŠ—Z⏞iβŠ—...βŠ—I⏞0z_i \equiv Z_i = \overbrace{I}^{n-1}\otimes ... \otimes \overbrace{Z}^{i} \otimes ... \otimes \overbrace{I}^{0}

Sa ibang salita, ang variable ziz_i ay katumbas ng isang ZZ gate na kumikilos sa qubit ii. Bukod pa rito:

Zi∣xnβˆ’1β‹―x0⟩=zi∣xnβˆ’1β‹―x0βŸ©β†’βŸ¨xnβˆ’1β‹―x0∣Zi∣xnβˆ’1β‹―x0⟩=ziZ_i|x_{n-1}\cdots x_0\rangle = z_i|x_{n-1}\cdots x_0\rangle \rightarrow \langle x_{n-1}\cdots x_0 |Z_i|x_{n-1}\cdots x_0\rangle = z_i

Kaya ang observable na ating isasaalang-alang ay:

H^=βˆ‘i=0nβˆ‘j=0iwij2ZiZj\hat{H} = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} Z_iZ_j

na kung saan kailangan pa nating idagdag ang independyenteng termino sa bandang huli:

offset=βˆ’βˆ‘i=0nβˆ‘j=0iwij2\texttt{offset} = - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}

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): IIZZ+IZIZ+IZZI+ZIIZ+ZZIIIIZZ + IZIZ + IZZI + ZIIZ + ZZII. 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")

Output of the previous code cell

# 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")

Output of the previous code cell

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)

Output of the previous code cell

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
)

Output of the previous code cell

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