Instanțe și extensii

Acest capitol va acoperi mai mulți algoritmi cuantici variațional, printre care:

• Variational Quantum Eigensolver (VQE)
• Subspace Search VQE (SSVQE)
• Variational Quantum Deflation (VQD)
• Quantum Sampling Regression (QSR)

Folosind acești algoritmi, vom învăța mai multe idei de design care pot fi integrate în algoritmi variationali personalizați, cum ar fi ponderi, penalizări, supra-eșantionare și sub-eșantionare. Te încurajăm să experimentezi cu aceste concepte și să îți împărtășești descoperirile cu comunitatea.

Cadrul Qiskit patterns se aplică tuturor acestor algoritmi — dar vom evidenția explicit pașii doar în primul exemplu.

Variational Quantum Eigensolver (VQE)

VQE este unul dintre cei mai utilizați algoritmi cuantici variationali, oferind un șablon pe care ceilalți algoritmi îl pot extinde.

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

Prezentare teoretică

Structura VQE este simplă:

• Pregătește operatorii de referință $U_R$
  • Pornim din starea $|0\rangle$ și ajungem la starea de referință $|\rho\rangle$
• Aplică forma variațională $U_V(\vec\theta_{i,j})$ pentru a crea un ansatz $U_A(\vec\theta_{i,j})$
  • Trecem din starea $|\rho\rangle$ în $U_V(\vec\theta_{i,j})|\rho\rangle = |\psi(\vec\theta_{i,j})\rangle$
• Bootstrap la $i=0$ dacă avem o problemă similară (găsită de obicei prin simulare clasică sau eșantionare)
  • Fiecare optimizer va fi bootstrapped diferit, rezultând un set inițial de vectori de parametri $\Theta_0 := \\{ {\vec\theta_{0,j} | j \in \mathcal{J}_\text{opt}^0} \\}$ (de exemplu, dintr-un punct inițial $\vec\theta_0$).
• Evaluează funcția de cost $C(\vec\theta_{i,j}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ pentru toate stările pregătite pe un calculator cuantic.
• Folosește un optimizer clasic pentru a selecta următorul set de parametri $\Theta_{i+1}$.
• Repetă procesul până la convergență.

Acesta este un simplu ciclu de optimizare clasică în care evaluăm funcția de cost. Unii optimizatori pot necesita evaluări multiple pentru a calcula un gradient, a determina următoarea iterație sau a evalua convergența.

Iată exemplul pentru următorul operator observabil:

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

Implementare

# Added by doQumentation — required packages for this notebook
!pip install -q numpy qiskit qiskit-ibm-runtime scipy

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.circuit.library import TwoLocal
import numpy as np

theta_list = (2 * np.pi * np.random.rand(1, 8)).tolist()
observable = SparsePauliOp.from_list([("II", 2), ("XX", -2), ("YY", 3), ("ZZ", -3)])

reference_circuit = QuantumCircuit(2)
reference_circuit.x(0)

variational_form = TwoLocal(
    2,
    rotation_blocks=["rz", "ry"],
    entanglement_blocks="cx",
    entanglement="linear",
    reps=1,
)

ansatz = reference_circuit.compose(variational_form)

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

def cost_func_vqe(parameters, 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
    """

    estimator_job = estimator.run([(ansatz, hamiltonian, [parameters])])
    estimator_result = estimator_job.result()[0]

    cost = estimator_result.data.evs[0]
    return cost

from qiskit.primitives import StatevectorEstimator

estimator = StatevectorEstimator()

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

# SciPy minimizer routine
from scipy.optimize import minimize

x0 = np.ones(8)

result = minimize(
    cost_func_vqe, x0, args=(ansatz, observable, estimator), method="COBYLA"
)

result

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -5.999999982445723
       x: [ 1.741e+00  9.606e-01  1.571e+00  2.115e-05  1.899e+00
            1.243e+00  6.063e-01  6.063e-01]
    nfev: 136
   maxcv: 0.0

Pasul 2: Optimizarea problemei pentru execuție cuantică

Vom selecta Backend-ul cel mai puțin ocupat și vom importa componentele necesare din qiskit_ibm_runtime.

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

<IBMBackend('ibm_brisbane')>

Vom transpila Circuit-ul folosind managerul de pași preset cu nivelul de optimizare 3 și vom aplica layout-ul corespunzător operatorului observabil.

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 = observable.apply_layout(layout=isa_ansatz.layout)

Pasul 3: Executare folosind primitivele Qiskit Runtime

Suntem acum pregătiți să rulăm calculul pe hardware IBM Quantum®. Deoarece minimizarea funcției de cost este puternic iterativă, vom porni o sesiune Runtime. Astfel, va trebui să așteptăm în coadă o singură dată. Odată ce job-ul începe să ruleze, fiecare iterație cu actualizări ale parametrilor va rula imediat.

x0 = np.ones(8)

estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    result = minimize(
        cost_func_vqe,
        x0,
        args=(isa_ansatz, isa_observable, estimator),
        method="COBYLA",
        options={"maxiter": 200, "disp": True},
    )
session.close()
print(result)

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

Putem observa că rutina de minimizare s-a terminat cu succes, ceea ce înseamnă că am atins toleranța implicită a optimizatorului clasic COBYLA. Dacă avem nevoie de un rezultat mai precis, putem specifica o toleranță mai mică. Aceasta poate fi într-adevăr situația, deoarece rezultatul a diferit cu câteva procente față de cel obținut de simulator mai sus.

Valoarea lui x obținută este cea mai bună estimare curentă a parametrilor care minimizează funcția de cost. Dacă iterăm pentru a obține o precizie mai mare, acele valori ar trebui folosite în locul lui x0 utilizat inițial (un vector de unu).

În final, remarcăm că funcția a fost evaluată de 96 de ori în procesul de optimizare. Aceasta poate diferi de numărul de pași de optimizare, deoarece unii optimizatori necesită evaluări multiple ale funcției într-un singur pas, cum ar fi atunci când estimează un gradient.

Subspace Search VQE (SSVQE)

SSVQE este o variantă a VQE care permite obținerea primelor $k$ valori proprii ale unui operator observabil $\hat{H}$ cu valorile proprii $\{\lambda_0, \lambda_1,...,\lambda_{N-1}\}$, unde $N\geq k$. Fără pierdere de generalitate, presupunem că $\lambda_0<\lambda_1<...<\lambda_{N-1}$. SSVQE introduce o idee nouă prin adăugarea de ponderi pentru a prioritiza optimizarea termenului cu ponderea cea mai mare.

Pentru a implementa acest algoritm, avem nevoie de $k$ stări de referință mutual ortogonale $\{ |\rho_j\rangle \}_{j=0}^{k-1}$, adică $\langle \rho_j | \rho_l \rangle = \delta_{jl}$ pentru $j,l<k$. Aceste stări pot fi construite folosind operatori Pauli. Funcția de cost a acestui algoritm este:

$$\begin{aligned} C(\vec{\theta}) & := \sum_{j=0}^{k-1} w_j \langle \rho_j | U_{V}^{\dagger}(\vec{\theta})\hat{H} U_{V}(\vec{\theta})|\rho_j \rangle \\[1mm] & := \sum_{j=0}^{k-1} w_j \langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle \\[1mm] \end{aligned}$$

unde $w_j$ este un număr pozitiv arbitrar astfel încât dacă $j<l<k$ atunci $w_j>w_l$, iar $U_V(\vec{\theta})$ este forma variațională definită de utilizator.

Algoritmul SSVQE se bazează pe faptul că stările proprii corespunzătoare unor valori proprii diferite sunt mutual ortogonale. Mai precis, produsul intern dintre $U_V(\vec{\theta})|\rho_j\rangle$ și $U_V(\vec{\theta})|\rho_l\rangle$ poate fi exprimat ca:

$$\begin{aligned} \langle \rho_j | U_{V}^{\dagger}(\vec{\theta})U_{V}(\vec{\theta})|\rho_l \rangle & = \langle \rho_j | I |\rho_l \rangle \\[1mm] & = \langle \rho_j | \rho_l \rangle \\[1mm] & = \delta_{jl} \end{aligned}$$

Prima egalitate este valabilă deoarece $U_{V}(\vec{\theta})$ este un operator cuantic și este, prin urmare, unitar. Ultima egalitate este valabilă datorită ortogonalității stărilor de referință $|\rho_j\rangle$. Faptul că ortogonalitatea este păstrată prin transformări unitare este strâns legat de principiul conservării informației, exprimat în știința informației cuantice. Din această perspectivă, transformările neunitare reprezintă procese în care informația fie se pierde, fie este injectată.

Ponderile $w_j$ ajută la asigurarea că toate stările sunt stări proprii. Dacă ponderile sunt suficient de diferite, termenul cu ponderea cea mai mare (adică $w_0$) va fi prioritizat în timpul optimizării față de ceilalți. Ca rezultat, starea $U_{V}(\vec{\theta})|\rho_0 \rangle$ va deveni starea proprie corespunzătoare lui $\lambda_0$. Deoarece $\{ U_{V}(\vec{\theta})|\rho_j\rangle \}_{j=0}^{k-1}$ sunt mutual ortogonale, stările rămase vor fi ortogonale față de ea și, prin urmare, conținute în subspațiul corespunzător valorilor proprii $\{\lambda_1,...,\lambda_{N-1}\}$.

Aplicând același argument pentru restul termenilor, următoarea prioritate va fi termenul cu ponderea $w_1$, astfel că $U_{V}(\vec{\theta})|\rho_1 \rangle$ va fi starea proprie corespunzătoare lui $\lambda_1$, iar ceilalți termeni vor fi conținuți în spațiul propriu al $\{\lambda_2,...,\lambda_{N-1}\}$.

Raționând inductiv, deducem că $U_{V}(\vec{\theta})|\rho_j \rangle$ va fi o stare proprie aproximativă a lui $\lambda_j$ pentru $0\leq j < k.$

Prezentare teoretică

SSVQE poate fi rezumat astfel:

• Pregătește mai multe stări de referință aplicând un unitar U_R la k stări diferite din baza computațională
  • Acest algoritm necesită utilizarea a $k$ stări de referință mutual ortogonale $\{ |\rho_j\rangle \}_{j=0}^{k-1}$, astfel încât $\langle \rho_j | \rho_l \rangle = \delta_{jl}$ pentru $j,l<k$.
• Aplică forma variațională $U_V(\vec\theta_{i,j})$ fiecărei stări de referință, rezultând următorul ansatz $U_A(\vec\theta_{i,j})$.
• Bootstrap la $i=0$ dacă o problemă similară este disponibilă (de obicei găsită prin simulare clasică sau eșantionare).
• Evaluează funcția de cost $C(\vec\theta_{i,j}) := \sum_{j=0}^{k-1} w_j \langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ pentru toate stările pregătite pe un calculator cuantic.
  • Aceasta poate fi separată în calcularea valorii de așteptare pentru un operator observabil $\langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ și înmulțirea acelui rezultat cu $w_j$.
  • Ulterior, funcția de cost returnează suma tuturor valorilor de așteptare ponderate.
• Folosește un optimizer clasic pentru a determina următorul set de parametri $\Theta_{i+1}$.
• Repetă pașii de mai sus până la convergență.

Vei reconstrui funcția de cost SSVQE în evaluare, dar avem următorul fragment de cod pentru a-ți motiva soluția:

import numpy as np

def cost_func_ssvqe(
    params, initialized_anastz_list, weights, ansatz, hamiltonian, estimator
):
    # """Return estimate of energy from estimator

    # Parameters:
    #     params (ndarray): Array of ansatz parameters
    #     initialized_anastz_list (list QuantumCircuit): Array of initialised ansatz with reference
    #     weights (list): List of weights
    #     ansatz (QuantumCircuit): Parameterized ansatz circuit
    #     hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
    #     estimator (Estimator): Estimator primitive instance

    # Returns:
    #     float: Weighted energy estimate
    # """

    energies = []

    # Define SSVQE

    weighted_energy_sum = np.dot(energies, weights)
    return weighted_energy_sum

Variational Quantum Deflation (VQD)

VQD este o metodă iterativă care extinde VQE pentru a obține primele $k$ valori proprii ale unui operator observabil $\hat{H}$ cu valorile proprii $\{\lambda_0, \lambda_1,...,\lambda_{N-1}\}$, unde $N\geq k$, în loc de doar prima. Pentru restul acestei secțiuni, vom presupune, fără pierdere de generalitate, că $\lambda_0\leq\lambda_1\leq...\leq\lambda_{N-1}$. VQD introduce noțiunea de cost de penalizare pentru a ghida procesul de optimizare.

VQD introduce un termen de penalizare, notat cu $\beta$, pentru a echilibra contribuția fiecărui termen de suprapunere la cost. Acest termen de penalizare servește la penalizarea procesului de optimizare dacă ortogonalitatea nu este atinsă. Impunem această constrângere deoarece stările proprii ale unui operator observabil sau ale unui operator hermitic, corespunzătoare unor valori proprii diferite, sunt întotdeauna mutual ortogonale, sau pot fi făcute să fie astfel în cazul degenerării sau al valorilor proprii repetate. Astfel, prin impunerea ortogonalității față de starea proprie corespunzătoare lui $\lambda_0$, optimizăm efectiv în subspațiul corespunzător restului valorilor proprii $\{\lambda_1, \lambda_2,..., \lambda_{N-1}\}$. Aici, $\lambda_1$ este cea mai mică valoare proprie din restul valorilor proprii și, prin urmare, soluția optimă a noii probleme poate fi obținută folosind teorema variațională.

Ideea generală din spatele VQD este de a folosi VQE ca de obicei pentru a obține cea mai mică valoare proprie $\lambda_0 := C_0(\vec\theta^0) \equiv C_\text{VQE}(\vec\theta^0)$ împreună cu starea proprie (aproximativă) corespunzătoare $|\psi(\vec{\theta^0})\rangle$ pentru un vector de parametri optim $\vec{\theta^0}$. Apoi, pentru a obține valoarea proprie următoare $\lambda_1 > \lambda_0$, în loc să minimizăm funcția de cost $C_0(\vec{\theta}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$, optimizăm:

$$C_1(\vec{\theta}) := C_0(\vec{\theta})+ \beta_0 |\langle \psi(\vec{\theta})| \psi(\vec{\theta^0})\rangle |^2$$

Valoarea pozitivă $\beta_0$ ar trebui să fie ideal mai mare decât $\lambda_1-\lambda_0$.

Aceasta introduce o nouă funcție de cost care poate fi privită ca o problemă cu constrângeri, unde minimizăm $C_\text{VQE}(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ sub constrângerea că starea trebuie să fie ortogonală față de $|\psi(\vec{\theta^0})\rangle$ obținut anterior, cu $\beta_0$ acționând ca termen de penalizare dacă constrângerea nu este satisfăcută.

Alternativ, această nouă problemă poate fi interpretată ca rularea VQE pe noul operator observabil:

$$\hat{H_1} := \hat{H} + \beta_0 |\psi(\vec{\theta^0})\rangle \langle \psi(\vec{\theta^0})| \quad \Rightarrow \quad C_1(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H_1} | \psi(\vec{\theta})\rangle,$$

Presupunând că soluția la noua problemă este $|\psi(\vec{\theta^1})\rangle$, valoarea așteptată a lui $\hat{H}$ (nu $\hat{H_1}$) ar trebui să fie $\langle \psi(\vec{\theta^1}) | \hat{H} | \psi(\vec{\theta^1})\rangle = \lambda_1$. Pentru a obține a treia valoare proprie $\lambda_2$, funcția de cost de optimizat este:

$$C_2(\vec{\theta}) := C_1(\vec{\theta}) + \beta_1 |\langle \psi(\vec{\theta})| \psi(\vec{\theta^1})\rangle |^2$$

unde $\beta_1$ este o constantă pozitivă suficient de mare pentru a impune ortogonalitatea stării soluție față de atât $|\psi(\vec{\theta^0})\rangle$, cât și $|\psi(\vec{\theta^1})\rangle$. Aceasta penalizează stările din spațiul de căutare care nu îndeplinesc această cerință, restricționând efectiv spațiul de căutare. Astfel, soluția optimă a noii probleme ar trebui să fie starea proprie corespunzătoare lui $\lambda_2$.

Ca și în cazul anterior, această nouă problemă poate fi interpretată și ca VQE cu operatorul observabil:

$$\hat{H_2} := \hat{H_1} + \beta_1 |\psi(\vec{\theta^1})\rangle \langle \psi(\vec{\theta^1})| \quad \Rightarrow \quad C_2(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H_2} | \psi(\vec{\theta})\rangle.$$

Dacă soluția la această nouă problemă este $|\psi(\vec{\theta^2})\rangle$, valoarea așteptată a lui $\hat{H}$ (nu $\hat{H_2}$) ar trebui să fie $\langle \psi(\vec{\theta^2}) | \hat{H} | \psi(\vec{\theta^2})\rangle = \lambda_2$. Prin analogie, pentru a obține a $k$-a valoare proprie $\lambda_{k-1}$, ai minimiza funcția de cost:

$$C_{k-1}(\vec{\theta}) := C_{k-2}(\vec{\theta}) + \beta_{k-2} |\langle \psi(\vec{\theta})| \psi(\vec{\theta^{k-2}})\rangle |^2,$$

Reamintim că am definit $\vec{\theta^j}$ astfel încât $\langle \psi(\vec{\theta^j}) | \hat{H} | \psi(\vec{\theta^j})\rangle = \lambda_j, \forall j<k$. Această problemă este echivalentă cu minimizarea lui $C(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ dar cu constrângerea că starea trebuie să fie ortogonală față de $|\psi(\vec{\theta^j})\rangle ; \forall j \in {0, \cdots, k-1}$, restricționând astfel spațiul de căutare la subspațiul corespunzător valorilor proprii $\{\lambda_{k-1},\cdots,\lambda_{N-1}\}$.

Această problemă este echivalentă cu un VQE cu operatorul observabil:

$$\hat{H}_{k-1} := \hat{H}_{k-2} + \beta_{k-2} |\psi(\vec{\theta^{k-2}})\rangle \langle \psi(\vec{\theta^{k-2}})| \quad \Rightarrow \quad C_{k-1}(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H}_{k-1} | \psi(\vec{\theta})\rangle,$$

După cum se poate observa din proces, pentru a obține a $k$-a valoare proprie, ai nevoie de stările proprii (aproximative) ale celor $k-1$ valori proprii anterioare, deci ar trebui să rulezi VQE de $k$ ori în total. Prin urmare, funcția de cost a VQD este:

$$C_k(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle + \sum_{j=0}^{k-1}\beta_j |\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2$$

Prezentare teoretică

Structura VQD poate fi rezumată astfel:

• Pregătește un operator de referință $U_R$
• Aplică forma variațională $U_V(\vec\theta_{i,j})$ la starea de referință, creând următoarele ansaetze $U_A(\vec\theta_{i,j})$
• Bootstrap la $i=0$ dacă avem o problemă similară (de obicei găsită prin simulare clasică sau eșantionare).
• Evaluează funcția de cost $C_k(\vec{\theta})$, care implică calcularea a $k$ stări excitate și un șir de $\beta$-uri care definesc penalizarea de suprapunere pentru fiecare termen de suprapunere.
  • Calculează valoarea de așteptare pentru un operator observabil $\langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ pentru fiecare $k$
  • Calculează penalizarea $\sum_{j=0}^{k-1}\beta_j |\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2$.
  • Funcția de cost ar trebui să returneze suma acestor doi termeni
• Folosește un optimizer clasic pentru a alege următorul set de parametri $\Theta_{i+1}$.
• Repetă acest proces până la convergență.

Implementare

Pentru această implementare, vom crea o funcție pentru o penalizare de suprapunere. Această penalizare va fi folosită în funcția de cost la fiecare iterație. Acest proces va fi repetat pentru fiecare stare excitată.

from qiskit.circuit.library import TwoLocal

ansatz = TwoLocal(2, rotation_blocks=["ry", "rz"], entanglement_blocks="cz", reps=1)

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

Mai întâi, vom configura o funcție care calculează fidelitatea stării — un procent de suprapunere între două stări, pe care îl vom folosi ca penalizare pentru VQD:

import numpy as np

def calculate_overlaps(ansatz, prev_circuits, parameters, sampler):
    def create_fidelity_circuit(circuit_1, circuit_2):
        """
        Constructs the list of fidelity circuits to be evaluated.
        These circuits represent the state overlap between pairs of input circuits,
        and their construction depends on the fidelity method implementations.
        """

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

Este momentul să scriem funcția de cost a VQD. Ca și înainte, când am calculat doar starea fundamentală, vom determina cea mai mică stare de energie folosind primitiva Estimator. Cu toate acestea, după cum s-a descris mai sus, vom adăuga acum un termen de penalizare pentru a asigura ortogonalitatea stărilor de energie mai înaltă. Adică, pentru fiecare nouă stare excitată, se adaugă o penalizare pentru orice suprapunere între starea variațională curentă și stările proprii de energie mai mică deja găsite.

def cost_func_vqd(
    parameters, ansatz, prev_states, step, betas, estimator, sampler, hamiltonian
):
    estimator_job = estimator.run([(ansatz, hamiltonian, [parameters])])

    total_cost = 0

    if step > 1:
        overlaps = calculate_overlaps(ansatz, prev_states, parameters, sampler)
        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

Remarcă în special că funcția de cost de mai sus face referire la funcția calculate_overlaps, care creează de fapt un nou Circuit cuantic. Dacă vrem să rulăm pe hardware real, acel nou Circuit trebuie să fie transpilat, de preferință într-un mod optim, pentru a rula pe Backend-ul pe care îl selectăm. Remarcă că transpilarea a fost inclusă în funcțiile calculate_overlaps sau cost_func_vqd. Simte-te liber să încerci să modifici codul pentru a include această transpilare suplimentară (și condiționată) — dar aceasta va fi realizată și în lecția următoare.

În această lecție, vom rula algoritmul VQD folosind Statevector Sampler și Statevector Estimator:

from qiskit.primitives import StatevectorEstimator as Estimator

sampler = Sampler()
estimator = Estimator()

Vom introduce un operator observabil de estimat. În lecția următoare, vom adăuga un context fizic, cum ar fi starea excitată a unei molecule. Poate fi util să te gândești la acest operator observabil ca la Hamiltonianul unui sistem care poate avea stări excitate, chiar dacă acesta nu a fost ales să corespundă niciunei molecule sau atom particular.

from qiskit.quantum_info import SparsePauliOp

observable = SparsePauliOp.from_list([("II", 2), ("XX", -2), ("YY", 3), ("ZZ", -3)])

Aici stabilim numărul total de stări pe care dorim să le calculăm (starea fundamentală și stările excitate, k), precum și penalizările (betas) pentru suprapunerea dintre vectorii de stare care ar trebui să fie ortogonali. Consecințele alegerii unor betas prea mari sau prea mici vor fi explorate în lecția următoare. Deocamdată, vom folosi pur și simplu pe cele furnizate mai jos. Vom începe prin a folosi toți zero ca parametri inițiali. În propriile calcule, poate vei dori să folosești parametri inițiali mai inteligenți bazați pe cunoașterea sistemului sau pe calcule anterioare.

k = 3
betas = [33, 33, 33]
x0 = np.zeros(8)

Putem acum rula calculul:

from scipy.optimize import minimize

prev_states = []
prev_opt_parameters = []
eigenvalues = []

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),
        method="COBYLA",
        options={
            "maxiter": 200,
        },
    )
    print(result)

    prev_opt_parameters = result.x
    eigenvalues.append(result.fun)

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -5.999999979545955
       x: [-5.150e-01 -5.452e-02 -1.571e+00 -2.853e-05  2.671e-01
           -2.672e-01 -8.509e-01 -8.510e-01]
    nfev: 131
   maxcv: 0.0
 message: Optimization terminated successfully.
 success: True
  status: 1
     fun: 4.024550284767612
       x: [-3.745e-01  1.041e+00  8.637e-01  1.202e+00 -8.847e-02
            1.181e-02  7.611e-01 -3.006e-01]
    nfev: 110
   maxcv: 0.0
 message: Optimization terminated successfully.
 success: True
  status: 1
     fun: 5.608925562838559
       x: [-2.670e-01  1.280e+00  1.070e+00 -8.031e-01 -1.524e-01
           -6.956e-02  7.018e-01  1.514e+00]
    nfev: 90
   maxcv: 0.0

Valorile obținute din funcția de cost sunt aproximativ -6,00, 4,02 și 5,61. Lucrul important în legătură cu aceste rezultate este că valorile funcției sunt crescătoare. Dacă am fi obținut o primă stare excitată cu o energie mai mică decât calculul inițial, neconstrâns al stării fundamentale, asta ar fi indicat o eroare undeva în codul nostru.

Valorile lui x sunt parametrii care au produs un vector de stare corespunzător fiecăruia dintre aceste costuri (energii).

În final, remarcăm că toate trei minimizările au conversat în toleranța implicită a optimizatorului clasic (în acest caz COBYLA). Acestea au necesitat 131, 110 și, respectiv, 90 de evaluări ale funcției.

Quantum Sampling Regression (QSR)

Una dintre principalele probleme ale VQE o reprezintă apelurile multiple către un calculator cuantic necesare pentru a obține parametrii la fiecare pas, de exemplu $k$, $k-1$ și așa mai departe. Aceasta este deosebit de problematică atunci când accesul la dispozitivele cuantice este în coadă de așteptare. Deși o Session poate fi folosită pentru a grupa mai multe apeluri iterative, o abordare alternativă este utilizarea eșantionării. Prin utilizarea mai multor resurse clasice, putem finaliza întregul proces de optimizare într-un singur apel. Aceasta este situația în care intervine Quantum Sampling Regression. Deoarece accesul la calculatoarele cuantice este încă o marfă cu ofertă redusă/cerere ridicată, considerăm că acest compromis este atât posibil, cât și convenabil pentru multe studii actuale. Această abordare valorifică toate capacitățile clasice disponibile, capturând în același timp multe din mecanismele interne și proprietățile intrinseci ale calculelor cuantice care nu apar în simulare.

Ideea din spatele QSR este că funcția de cost $C(\theta) := \langle \psi(\theta) | \hat{H} | \psi(\theta)\rangle$ poate fi exprimată ca o serie Fourier în felul următor:

$$\begin{aligned} C(\vec{\theta}) & := \langle \psi(\theta) | \hat{H} | \psi(\theta)\rangle \\[1mm] & := a_0 + \sum_{k=1}^S[a_k\cos(k\theta)+ b_k\sin(k\theta)] \\[1mm] \end{aligned}$$

În funcție de periodicitatea și lățimea de bandă a funcției originale, mulțimea $S$ poate fi finită sau infinită. În scopul acestei discuții, vom presupune că este infinită. Următorul pas este să eșantionăm funcția de cost $C(\theta)$ de mai multe ori pentru a obține coeficienții Fourier $\{a_0, a_k, b_k\}_{k=1}^S$. Mai precis, deoarece avem $2S+1$ necunoscute, va trebui să eșantionăm funcția de cost de $2S+1$ ori.

Dacă eșantionăm funcția de cost pentru $2S+1$ valori de parametri $\{\theta_1,...,\theta_{2S+1}\}$, putem obține următorul sistem:

$$\begin{pmatrix} 1 & \cos(\theta_1) & \sin(\theta_1) & \cos(2\theta_1) & ... & \sin(S\theta_1) \\ 1 & \cos(\theta_2) & \sin(\theta_2) & \cos(2\theta_2) & \cdots & \sin(S\theta_2)\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & \cos(\theta_{2S+1}) & \sin(\theta_{2S+1}) & \cos(2\theta_{2S+1}) & \cdots & \sin(S\theta_{2S+1}) \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ b_1 \\ a_2 \\ \vdots \\ b_S \end{pmatrix} = \begin{pmatrix} C(\theta_1) \\ C(\theta_2) \\ \vdots \\ C(\theta_{2S+1}) \end{pmatrix},$$

pe care îl vom rescrie ca

$$Fa=c.$$

În practică, acest sistem nu este în general consistent deoarece valorile funcției de cost $c$ nu sunt exacte. Prin urmare, este de obicei o idee bună să le normalizăm înmulțindu-le la stânga cu $F^\dagger$, ceea ce rezultă în:

$$F^\dagger Fa = F^\dagger c.$$

Acest nou sistem este întotdeauna consistent, iar soluția lui este o soluție în sensul celor mai mici pătrate pentru problema originală. Dacă avem $k$ parametri în loc de unul singur, și fiecare parametru $\theta^i$ are propriul $S_i$ pentru $i \in {1,...,k}$, atunci numărul total de eșantioane necesar este:

$$T=\prod_{i=1}^k(2S_i+1)\leq \prod_{i=1}^k(2S_{max}+1) = (2S_{max}+1)^n,$$

unde $S_{\max} = \max_i(S_i)$. În plus, ajustarea lui $S_{\max}$ ca parametru reglabil (în loc să fie inferit) deschide noi posibilități, cum ar fi:

• Supra-eșantionarea pentru a îmbunătăți precizia.
• Sub-eșantionarea pentru a crește performanța prin reducerea supraîncărcărilor de rulare sau eliminarea minimelor locale.

Prezentare teoretică

Structura QSR poate fi rezumată astfel:

• Pregătește operatorii de referință $U_R$.
  • Vom trece din starea $|0\rangle$ la starea de referință $|\rho\rangle$
• Aplică forma variațională $U_V(\vec\theta_{i,j})$ pentru a crea un ansatz $U_A(\vec\theta_{i,j})$.
  • Determină lățimea de bandă asociată fiecărui parametru din ansatz. O limită superioară este suficientă.
• Bootstrap la $i=0$ dacă avem o problemă similară (de obicei găsită prin simulare clasică sau eșantionare).
• Eșantionează funcția de cost $C(\vec\theta) := a_0 + \sum_{k=1}^S[a_k\cos(k\theta)+ b_k\sin(k\theta)]$ de cel puțin $T$ ori.
  • $T=\prod_{i=1}^k(2S_i+1)\leq \prod_{i=1}^k(2S_{max}+1) = (2S_{max}+1)^n$
  • Decide dacă să supra-eșantionezi/sub-eșantionezi pentru a echilibra viteza față de precizie prin ajustarea lui $T$.
• Calculează coeficienții Fourier din eșantioane (adică, rezolvă sistemul normalizat de ecuații liniare).
• Rezolvă pentru minimul global al funcției de regresie rezultante pe un calculator clasic.

Rezumat

Cu această lecție, ai aflat despre mai multe instanțe variaționale disponibile:

• Structura generală
• Introducerea ponderilor și penalizărilor pentru a ajusta o funcție de cost
• Explorarea sub-eșantionării față de supra-eșantionarea pentru a face compromisul între viteză și precizie