Exemple și aplicații

În această lecție, vom explora câteva exemple de algoritmi variaționali și cum să îi aplicăm:

• Cum să scrii un algoritm variațional personalizat
• Cum să aplici un algoritm variațional pentru a găsi valorile proprii minime
• Cum să utilizezi algoritmi variaționali pentru a rezolva cazuri de utilizare din aplicații reale

Reține că framework-ul Qiskit patterns poate fi aplicat tuturor problemelor pe care le introducem aici. Cu toate acestea, pentru a evita repetiția, vom menționa explicit pașii framework-ului doar într-un singur exemplu, rulat pe hardware real.

Definiții ale problemelor

Imaginează-ți că vrem să folosim un algoritm variațional pentru a găsi valoarea proprie a următorului observabil:

$$\hat{O}_1 = 2 II - 2 XX + 3 YY - 3 ZZ,$$

Acest observabil are următoarele valori proprii:

$$\left\{ \begin{array}{c} \lambda_0 = -6 \\ \lambda_1 = 4 \\ \lambda_2 = 4 \\ \lambda_3 = 6 \end{array} \right\}$$

Și stările proprii:

$$\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)])

VQE personalizat

Vom explora mai întâi cum să construim o instanță VQE manual pentru a găsi cea mai mică valoare proprie a lui $\hat{O}_1$. Aceasta va încorpora o varietate de tehnici pe care le-am abordat pe parcursul acestui curs.

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

Vom începe depanarea pe simulatoare locale.

from qiskit.primitives import StatevectorEstimator as Estimator
from qiskit.primitives import StatevectorSampler as Sampler

estimator = Estimator()
sampler = Sampler()

Acum stabilim un set inițial de parametri:

x0 = np.ones(raw_ansatz.num_parameters)
print(x0)

[1. 1. 1. 1. 1. 1. 1. 1.]

Putem minimiza această funcție de cost pentru a calcula parametrii optimi

# 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

Deoarece această problemă de test folosește doar doi qubiți, putem verifica acest lucru utilizând rezolvitorul de valori proprii din algebra liniară al 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

După cum poți vedea, rezultatul este extrem de aproape de cel ideal.

Experimentând pentru a îmbunătăți viteza și acuratețea

Adaugă starea de referință

În exemplul anterior nu am folosit niciun operator de referință $U_R$. Acum să ne gândim cum poate fi obținută starea proprie ideală $\frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)$. Consideră următorul circuit.

from qiskit import QuantumCircuit

ideal_qc = QuantumCircuit(2)
ideal_qc.h(0)
ideal_qc.cx(0, 1)

ideal_qc.draw("mpl")

Putem verifica rapid că acest circuit ne oferă starea dorită.

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

Acum că am văzut cum arată un circuit care pregătește starea soluție, pare rezonabil să folosim o poartă Hadamard ca circuit de referință, astfel încât ansatz-ul complet devine:

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

Pentru acest nou circuit, soluția ideală ar putea fi atinsă cu toți parametrii setați la $0$, deci aceasta confirmă că alegerea circuitului de referință este rezonabilă.

Acum să comparăm numărul de evaluări ale funcției de cost, iterațiile optimizatorului și timpul necesar cu cele din tentativa anterioară.

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

Folosind parametrii optimi pentru a calcula valoarea proprie minimă:

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

În funcție de sistemul tău specific, acest lucru poate sau nu să ducă la o îmbunătățire a vitezei sau acurateței în acest exemplu la scară foarte mică. Ideea este că pornirea cu stări de referință motivate fizic devine din ce în ce mai importantă pentru îmbunătățirea vitezei și acurateței pe măsură ce problemele se scalează.

Schimbă punctul inițial

Acum că am văzut efectul adăugării stării de referință, vom analiza ce se întâmplă când alegem puncte inițiale diferite $\vec{\theta_0}$. În particular vom folosi $\vec{\theta_0}=(0,0,0,0,6,0,0,0)$ și $\vec{\theta_0}=(6,6,6,6,6,6,6,6,6)$.

Ține minte că, după cum s-a discutat când a fost introdusă starea de referință, soluția ideală ar fi găsită când toți parametrii sunt $0$, deci primul punct inițial ar trebui să ducă la mai puține evaluări.

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

Ajustând punctul inițial la $\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

Experimentând cu puncte inițiale diferite, s-ar putea să poți atinge convergența mai rapid și cu mai puține evaluări ale funcției.

Experimentând cu optimizatori diferiți

Putem ajusta optimizatorul folosind argumentul method al funcției minimize din SciPy, cu mai multe opțiuni găsite aici. Am folosit inițial un minimizator cu constrângeri (COBYLA). În acest exemplu, vom explora folosirea unui minimizator fără constrângeri (BFGS)

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

Exemplu VQD

Aici implementăm cadrul Qiskit patterns în mod explicit.

Pasul 1: Maparea intrărilor clasice la o problemă cuantică

Acum, în loc să căutăm doar cea mai mică valoare proprie a observabilelor noastre, vom căuta toate cele $4$ (unde $k=4$).

Reamintește-ți că funcțiile de cost ale VQD sunt:

$$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\}$$

Acest lucru este deosebit de important deoarece un vector $\vec{\beta}=(\beta_0,\cdots,\beta_{k-1})$ (în acest caz $(\beta_0, \beta_1, \beta_2, \beta_3)$) trebuie transmis ca argument atunci când definim obiectul VQD.

De asemenea, în implementarea VQD din Qiskit, în loc să se ia în considerare observabilele efective descrise în lecția anterioară, fidelitățile $|\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2$ sunt calculate direct prin algoritmul ComputeUncompute, care utilizează o primitivă Sampler pentru a eșantiona probabilitatea de a obține $|0\rangle$ pentru circuitul $U_A^\dagger(\vec{\theta})U_A(\vec{\theta^j})$. Aceasta funcționează tocmai pentru că această probabilitate este

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

Să începem prin examinarea următorului observabil:

$$\hat{O}_2 := 2 II - 3 XX + 2 YY - 4 ZZ$$

Acest observabil are următoarele valori proprii:

$$\left\{ \begin{array}{c} \lambda_0 = -7 \\ \lambda_1 = 3\\ \lambda_2 = 5 \\ \lambda_3 = 7 \end{array} \right\}$$

Și stările proprii:

$$\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)])

Vom folosi următoarea funcție pentru a calcula penalitatea de suprapunere. Reține că aceasta face parte în continuare din maparea problemei la circuite cuantice. Cu toate acestea, după cum s-a discutat în lecția anterioară, această funcție calculează suprapunerea dintre un circuit variațional curent și circuitul optimizat dintr-o stare anterioară de energie/cost mai mică obținută. Noul circuit generat trebuie, de asemenea, transpilat pentru a rula pe hardware real. Am mai văzut această funcție, utilizată pe un simulator. Aici, trebuie să luăm deja în considerare transpilarea și optimizarea aferentă pentru atunci când folosim un backend real, de unde și liniile cu if realbackend == 1. Aceasta amestecă puțin pasul 2, dar vom preciza explicit pasul 2 mai târziu.

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)

Acum adăugăm funcția de cost a VQD. Reține că, față de lecția anterioară, avem acum două argumente suplimentare (realbackend și backend) pentru a ne ajuta cu transpilarea atunci când folosim backend-uri reale.

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

Vom folosi din nou simulatoare pentru depanare, apoi vom trece la hardware real.

from qiskit.primitives import StatevectorSampler
from qiskit.primitives import StatevectorEstimator

sampler = StatevectorSampler(default_shots=4092)
estimator = StatevectorEstimator()

Aici introducem numărul de stări pe care dorim să le calculăm, penalitățile și un set de parametri inițiali, x0.

from qiskit.quantum_info import SparsePauliOp

k = 4
betas = [50, 60, 40]
x0 = np.ones(8)

Vom testa acum algoritmul folosind simulatoare:

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

Aceste rezultate sunt destul de aproape de cele așteptate, cu excepția erorii de aproximare și a fazei globale. Am putea ajusta toleranța optimizatorului clasic și/sau penalitățile pentru suprapunerea vectorilor de stare pentru a obține valori mai precise.

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

Modificarea valorilor beta

După cum s-a menționat în lecția anterioară, valorile lui $\vec{\beta}$ ar trebui să fie mai mari decât diferența dintre valorile proprii. Să vedem ce se întâmplă atunci când nu satisfac această condiție cu $\hat{O}_2$

$$\hat{O}_2 = 2 II - 3 XX + 2 YY - 4 ZZ$$

cu valorile proprii

$$\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

De data aceasta, optimizatorul returnează aceeași stare $|\phi_0\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)$ ca soluție propusă pentru toate stările proprii: ceea ce este în mod evident greșit. Acest lucru se întâmplă deoarece valorile beta erau prea mici pentru a penaliza starea proprie minimă în funcțiile de cost succesive. Prin urmare, aceasta nu a fost exclusă din spațiul de căutare efectiv în iterațiile ulterioare ale algoritmului și a fost aleasă întotdeauna ca cea mai bună soluție posibilă.

Recomandăm să experimentezi cu valorile lui $\vec{\beta}$ și să te asiguri că sunt mai mari decât diferența dintre valorile proprii.

Pasul 2: Optimizarea problemei pentru execuție cuantică

Pentru a rula aceasta pe hardware real, trebuie să optimizăm circuitele cuantice pentru calculatorul cuantic ales. În scopurile noastre, vom folosi pur și simplu backend-ul cel mai puțin ocupat.

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

Vom transpila circuitul nostru folosind un manager de pași preset și nivelul de optimizare 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)

Pasul 3: Execuție folosind primitivele Qiskit

Asigurând că resetăm valorile beta la valori suficient de mari, putem rula acum calculul pe hardware cuantic real.

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

Pasul 4: Post-procesare, returnarea rezultatului în format clasic

Rezultatul nostru este structural similar cu ceea ce s-a discutat în lecțiile și exemplele anterioare. Dar există ceva problematic în rezultatele de mai sus, din care putem deriva un mesaj de avertizare pentru contextul stărilor excitate. Pentru a limita timpul de calcul folosit în acest exemplu de învățare, am setat un număr maxim de iterații pentru optimizatorul clasic care era potențial prea mic: 200 de iterații. Un calcul anterior de mai sus, pe un simulator, nu a reușit să conveargă în 200 de iterații. Aici, al nostru a convergit... dar la ce toleranță? Nu am specificat o toleranță pentru ca COBYLA să se considere „convergit". O privire rapidă la valoarea funcției și compararea cu rulările anterioare ne spune că COBYLA nu era aproape de a converge la precizia de care avem nevoie.

Există o altă problemă: energia primei stări excitate pare a fi mai mică decât energia stării fundamentale! Încearcă să explici cum ar putea fi posibil acest lucru. Indiciu: este legat de punctul de convergență pe care tocmai l-am abordat. Acest comportament este explicat în detaliu mai jos după aplicarea VQD la molecula H2.

Chimie cuantică: rezolvator de stare fundamentală și stări excitate

Obiectivul nostru este să minimizăm valoarea de așteptare a observabilului care reprezintă energia (Hamiltonianul $\hat{\mathcal{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")

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

Acum stabilim un set inițial de parametri:

import numpy as np

x0 = np.ones(chem_ansatz.num_parameters)

Putem minimiza această funcție de cost pentru a calcula parametrii optimi și putem verifica mai întâi codul folosind un simulator local.

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

Valoarea minimă a funcției de cost (-1.857...) este energia stării fundamentale a moleculei H2, exprimată în unități de hartree.

Stări excitate

Putem folosi și VQD pentru a rezolva $k=2$ stări totale (starea fundamentală și prima stare excitată).

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,
]

Vom adăuga calculul de suprapunere:

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]

Hardware real și un avertisment final

Pentru a rula acest lucru pe hardware real, trebuie să optimizăm circuitele cuantice pentru calculatorul cuantic ales. În scopurile noastre, vom folosi pur și simplu backend-ul cel mai puțin ocupat.

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)

Vom folosi un manager de pași presetat pentru transpilare și vom optimiza maximal circuitul nostru folosind nivelul de optimizare 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)

Deoarece VQD este extrem de iterativ, vom realiza toți pașii în interiorul unei sesiuni Runtime, astfel încât joburile noastre să fie puse în coadă doar la început, nu între fiecare actualizare de parametri. Nimic altceva nu se schimbă în sintaxa funcției de cost sau a estimatorului.

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)

Energia stării fundamentale obținută (-1.83 hartree) nu este prea departe de valoarea corectă (-1.85 hartree). Totuși, energia stării excitate este destul de diferită. Aceasta seamănă cu comportamentul eronat pe care l-am văzut mai devreme în această lecție. Energia raportată pentru starea excitată este aproape aceeași cu cea a stării fundamentale. În cazul anterior, am văzut chiar o energie a stării excitate mai mică decât energia raportată a stării fundamentale.

Este imposibil ca un calcul variațional să producă o energie mai mică decât energia adevărată a stării fundamentale. În instanța anterioară, energia stării fundamentale obținută nu era foarte apropiată de cea adevărată. Deoarece nu am obținut energia adevărată a stării fundamentale în acel caz, nu există nicio contradicție. În cazul de față, energia stării fundamentale era destul de apropiată de valoarea corectă, și totuși energia stării excitate pare ciudat de apropiată de aceeași valoare.

Pentru a înțelege mai bine cum s-a întâmplat asta, reamintește-ți că modul în care găsim o stare excitată este prin a impune că starea variațională fie ortogonală cu starea fundamentală (folosind circuite de suprapunere și termeni de penalizare). Dacă nu reușim să obținem o energie a stării fundamentale precisă (sau suntem cu câteva procente departe), atunci nu reușim nici să obținem un vector al stării fundamentale precis! Astfel, atunci când impunem că starea excitată fie ortogonală cu prima stare găsită, nu impuneam ortogonalitate cu adevărata stare fundamentală, ci cu o aproximare a ei (uneori o aproximare slabă). Prin urmare, starea excitată nu a fost forțată să fie ortogonală cu adevărata stare fundamentală, iar estimările noastre de energie pentru stările excitate au fost de fapt foarte apropiate de energia stării fundamentale.

Aceasta va fi întotdeauna o preocupare în VQD. Dar în principiu, acest lucru poate fi corectat prin creșterea numărului maxim de iterații pentru optimizatorul clasic, impunând o toleranță mai mică pentru optimizatorul clasic și posibil încercând un ansatz diferit dacă ratăm în mod constant adevărata stare fundamentală. Cum am văzut, poate fi necesar și să modifici penalizările de suprapunere (betas). Dar aceasta este cu adevărat o problemă separată. Nicio penalizare pentru suprapunere nu te va ține departe de adevărata stare fundamentală, dacă nu ai găsit o estimare foarte bună a adevăratei stări fundamentale pentru circuitul de suprapunere.

Optimizare: Max-Cut

Problema tăieturii maxime (Max-Cut) este o problemă de optimizare combinatorie care implică împărțirea vârfurilor unui graf în două mulțimi disjuncte astfel încât numărul de muchii dintre cele două mulțimi să fie maximizat. Mai formal, dat un graf neorientat $G=(V,E)$, unde $V$ este mulțimea vârfurilor și $E$ este mulțimea muchiilor, problema Max-Cut cere partiționarea vârfurilor în două submulțimi disjuncte, $S$ și $T$, astfel încât numărul de muchii cu un capăt în $S$ și celălalt în $T$ să fie maximizat.

Putem aplica Max-Cut pentru a rezolva diverse probleme, cum ar fi clusterizarea, proiectarea rețelelor și tranzițiile de fază. Vom începe prin crearea unui graf de problemă:

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

Această problemă poate fi exprimată ca o problemă de optimizare binară. Pentru fiecare nod $0 \leq i < n$, unde $n$ este numărul de noduri ale grafului (în acest caz $n=4$), vom considera variabila binară $x_i$. Această variabilă va avea valoarea $1$ dacă nodul $i$ se află în unul dintre grupuri, pe care îl vom eticheta $1$, și $0$ dacă se află în celălalt grup, pe care îl vom eticheta $0$. Vom nota de asemenea cu $w_{ij}$ (elementul $(i,j)$ al matricei de adiacență $w$) ponderea muchiei care merge de la nodul $i$ la nodul $j$. Deoarece graful este neorientat, $w_{ij}=w_{ji}$. Atunci putem formula problema noastră ca maximizarea următoarei funcții de cost:

$$\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}$$

Pentru a rezolva această problemă cu un calculator cuantic, vom exprima funcția de cost ca valoarea de așteptare a unui observabil. Totuși, observabilele pe care Qiskit le acceptă nativ constau din operatori Pauli, care au valorile proprii $1$ și $-1$ în loc de $0$ și $1$. De aceea vom face următoarea schimbare de variabilă:

Unde $\vec{x}=(x_0,x_1,\cdots ,x_{n-1})$. Putem folosi matricea de adiacență $w$ pentru a accesa comod ponderile tuturor muchiilor. Aceasta va fi utilizată pentru a obține funcția noastră de cost:

$$z_i = 1-2x_i \rightarrow x_i = \frac{1-z_i}{2}$$

Aceasta implică că:

$$\begin{array}{lcl} x_i=0 & \rightarrow & z_i=1 \\ x_i=1 & \rightarrow & z_i=-1.\end{array}$$

Deci noua funcție de cost pe care vrem să o maximizăm este:

$$\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}$$

Mai mult, tendința naturală a unui calculator cuantic este de a găsi minime (de obicei, energia cea mai mică) în loc de maxime, deci în loc să maximizăm $C(\vec{z})$ vom minimiza:

$$-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}$$

Acum că avem o funcție de cost de minimizat ale cărei variabile pot lua valorile $-1$ și $1$, putem face următoarea analogie cu Pauli $Z$:

$$z_i \equiv Z_i = \overbrace{I}^{n-1}\otimes ... \otimes \overbrace{Z}^{i} \otimes ... \otimes \overbrace{I}^{0}$$

Cu alte cuvinte, variabila $z_i$ va fi echivalentă cu o poartă $Z$ care acționează pe qubit-ul $i$. Mai mult:

$$Z_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$$

Atunci observabilul pe care îl vom considera este:

$$\hat{H} = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} Z_iZ_j$$

la care va trebui să adăugăm termenul independent ulterior:

$$\texttt{offset} = - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}$$

Operatorul este o combinație liniară de termeni cu operatori Z pe noduri conectate printr-o muchie (reamintește-ți că qubit-ul 0 se află cel mai în dreapta): $IIZZ + IZIZ + IZZI + ZIIZ + ZZII$. Odată ce operatorul este construit, ansatz-ul pentru algoritmul QAOA poate fi ușor construit folosind circuitul QAOAAnsatz din biblioteca de circuite Qiskit.

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

Acum stabilim un set inițial de parametri aleatori:

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]

Orice optimizator clasic poate fi folosit pentru a minimiza funcția de cost. Pe un sistem cuantic real, un optimizator proiectat pentru peisaje de funcții de cost non-netede funcționează de obicei mai bine. Aici folosim rutina COBYLA din SciPy prin funcția minimize.

Deoarece executăm iterativ multe apeluri la Runtime, folosim o sesiune pentru a executa toate apelurile în cadrul unui singur bloc. Mai mult, pentru QAOA, soluția este codificată în distribuția de ieșire a circuitului ansatz legat cu parametrii optimi din minimizare. Prin urmare, vom avea nevoie de primitiva Sampler și o vom instanția cu aceeași sesiune. Și rulăm rutina noastră de minimizare:

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

Vectorul soluție al unghiurilor de parametri (x), când este introdus în circuitul ansatz, produce partiționarea grafului pe care o căutam.

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
)

Rezumat

Cu această lecție, ai învățat:

• Cum să scrii un algoritm variațional personalizat
• Cum să aplici un algoritm variațional pentru a găsi valorile proprii minime
• Cum să utilizezi algoritmi variațional pentru a rezolva cazuri de utilizare a aplicațiilor