Istanze ed estensioni

Questo capitolo tratta diversi algoritmi variazionali quantistici, tra cui:

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

Utilizzando questi algoritmi, impareremo diverse idee progettuali che possono essere incorporate in algoritmi variazionali personalizzati, come pesi, penalità, sovra-campionamento e sotto-campionamento. Ti incoraggiamo a sperimentare con questi concetti e a condividere i tuoi risultati con la comunità.

Il framework dei pattern Qiskit si applica a tutti questi algoritmi — ma indicheremo esplicitamente i passaggi solo nel primo esempio.

Variational Quantum Eigensolver (VQE)

Il VQE è uno degli algoritmi quantistici variazionali più diffusi e stabilisce un modello su cui altri algoritmi possono basarsi.

Passo 1: Mappare gli input classici in un problema quantistico

Schema teorico

La struttura del VQE è semplice:

• Preparare gli operatori di riferimento $U_R$
  • Si parte dallo stato $|0\rangle$ e si arriva allo stato di riferimento $|\rho\rangle$
• Applicare la forma variazionale $U_V(\vec\theta_{i,j})$ per creare un ansatz $U_A(\vec\theta_{i,j})$
  • Si passa dallo stato $|\rho\rangle$ a $U_V(\vec\theta_{i,j})|\rho\rangle = |\psi(\vec\theta_{i,j})\rangle$
• Effettuare il bootstrap a $i=0$ se si dispone di un problema simile (tipicamente trovato tramite simulazione classica o campionamento)
  • Ogni ottimizzatore viene inizializzato in modo diverso, producendo un insieme iniziale di vettori di parametri $\Theta_0 := \\{ {\vec\theta_{0,j} | j \in \mathcal{J}_\text{opt}^0} \\}$ (ad esempio, a partire da un punto iniziale $\vec\theta_0$).
• Valutare la funzione di costo $C(\vec\theta_{i,j}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ per tutti gli stati preparati su un computer quantistico.
• Utilizzare un ottimizzatore classico per selezionare il prossimo insieme di parametri $\Theta_{i+1}$.
• Ripetere il processo fino al raggiungimento della convergenza.

Questo è un semplice ciclo di ottimizzazione classica in cui si valuta la funzione di costo. Alcuni ottimizzatori possono richiedere più valutazioni per calcolare un gradiente, determinare l'iterazione successiva o verificare la convergenza.

Ecco l'esempio per il seguente osservabile:

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

Implementazione

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

Possiamo usare questa funzione di costo per calcolare i parametri ottimali

# 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

Passo 2: Ottimizzare il problema per l'esecuzione quantistica

Selezioneremo il backend meno occupato e importeremo i componenti necessari da 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')>

Traspileremo il circuito usando il pass manager preimpostato con livello di ottimizzazione 3, e applicheremo il layout corrispondente all'osservabile.

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)

Passo 3: Eseguire usando le primitive di Qiskit Runtime

Siamo ora pronti per eseguire il calcolo sull'hardware IBM Quantum®. Poiché la minimizzazione della funzione di costo è altamente iterativa, avvieremo una sessione Runtime. In questo modo, dovremo aspettare in coda una sola volta. Una volta che il job inizia l'esecuzione, ogni iterazione con aggiornamenti dei parametri verrà eseguita immediatamente.

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)

Passo 4: Post-elaborazione, restituzione del risultato in formato classico

Possiamo vedere che la routine di minimizzazione è terminata con successo, il che significa che abbiamo raggiunto la tolleranza predefinita dell'ottimizzatore classico COBYLA. Se richiedessimo un risultato più preciso, potremmo specificare una tolleranza più piccola. Questo potrebbe effettivamente essere il caso, dato che il risultato si discostava di diversi punti percentuali rispetto al risultato ottenuto dal simulatore sopra.

Il valore di x ottenuto è la migliore stima corrente per i parametri che minimizzano la funzione di costo. Se si itera per ottenere una precisione maggiore, quei valori dovrebbero essere usati al posto dell'x0 inizialmente utilizzato (un vettore di uni).

Infine, notiamo che la funzione è stata valutata 96 volte nel processo di ottimizzazione. Questo potrebbe differire dal numero di passi di ottimizzazione, poiché alcuni ottimizzatori richiedono più valutazioni della funzione in un singolo passo, ad esempio quando si stima un gradiente.

Subspace Search VQE (SSVQE)

Lo SSVQE è una variante del VQE che permette di ottenere i primi $k$ autovalori di un osservabile $\hat{H}$ con autovalori $\{\lambda_0, \lambda_1,...,\lambda_{N-1}\}$, dove $N\geq k$. Senza perdita di generalità, assumiamo che $\lambda_0<\lambda_1<...<\lambda_{N-1}$. SSVQE introduce una nuova idea aggiungendo pesi per aiutare a dare priorità all'ottimizzazione del termine con il peso maggiore.

Per implementare questo algoritmo, abbiamo bisogno di $k$ stati di riferimento mutuamente ortogonali $\{ |\rho_j\rangle \}_{j=0}^{k-1}$, ovvero tali che $\langle \rho_j | \rho_l \rangle = \delta_{jl}$ per $j,l<k$. Questi stati possono essere costruiti usando operatori di Pauli. La funzione di costo di questo algoritmo è quindi:

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

dove $w_j$ è un numero positivo arbitrario tale che se $j<l<k$ allora $w_j>w_l$, e $U_V(\vec{\theta})$ è la forma variazionale definita dall'utente.

L'algoritmo SSVQE si basa sul fatto che gli autostati corrispondenti a diversi autovalori sono mutuamente ortogonali. In particolare, il prodotto interno di $U_V(\vec{\theta})|\rho_j\rangle$ e $U_V(\vec{\theta})|\rho_l\rangle$ può essere espresso come:

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

La prima uguaglianza vale perché $U_{V}(\vec{\theta})$ è un operatore quantistico ed è quindi unitario. L'ultima uguaglianza vale per l'ortogonalità degli stati di riferimento $|\rho_j\rangle$. Il fatto che l'ortogonalità venga preservata dalle trasformazioni unitarie è profondamente legato al principio di conservazione dell'informazione, come espresso nella scienza dell'informazione quantistica. In questa prospettiva, le trasformazioni non unitarie rappresentano processi in cui l'informazione viene persa o iniettata.

I pesi $w_j$ aiutano a garantire che tutti gli stati siano autostati. Se i pesi sono sufficientemente diversi, il termine con il peso maggiore (ovvero $w_0$) riceverà priorità durante l'ottimizzazione rispetto agli altri. Di conseguenza, lo stato risultante $U_{V}(\vec{\theta})|\rho_0 \rangle$ diventerà l'autostato corrispondente a $\lambda_0$. Poiché $\{ U_{V}(\vec{\theta})|\rho_j\rangle \}_{j=0}^{k-1}$ sono mutuamente ortogonali, gli stati rimanenti saranno ortogonali ad esso e, quindi, contenuti nel sottospazio corrispondente agli autovalori $\{\lambda_1,...,\lambda_{N-1}\}$.

Applicando lo stesso ragionamento al resto dei termini, la priorità successiva sarebbe il termine con peso $w_1$, quindi $U_{V}(\vec{\theta})|\rho_1 \rangle$ sarebbe l'autostato corrispondente a $\lambda_1$, e gli altri termini sarebbero contenuti nell'autospazio di $\{\lambda_2,...,\lambda_{N-1}\}$.

Ragionando per induzione, deduciamo che $U_{V}(\vec{\theta})|\rho_j \rangle$ sarà un autostato approssimato di $\lambda_j$ per $0\leq j < k.$

Schema teorico

Lo SSVQE può essere riassunto come segue:

• Preparare diversi stati di riferimento applicando una unitaria U_R a $k$ diversi stati della base computazionale
  • Questo algoritmo richiede l'uso di $k$ stati di riferimento mutuamente ortogonali $\{ |\rho_j\rangle \}_{j=0}^{k-1}$, tali che $\langle \rho_j | \rho_l \rangle = \delta_{jl}$ per $j,l<k$.
• Applicare la forma variazionale $U_V(\vec\theta_{i,j})$ a ciascuno stato di riferimento, ottenendo il seguente ansatz $U_A(\vec\theta_{i,j})$.
• Effettuare il bootstrap a $i=0$ se è disponibile un problema simile (di solito trovato tramite simulazione classica o campionamento).
• Valutare la funzione di costo $C(\vec\theta_{i,j}) := \sum_{j=0}^{k-1} w_j \langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ per tutti gli stati preparati su un computer quantistico.
  • Questo può essere suddiviso nel calcolo del valore di aspettazione per un osservabile $\langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ e nella moltiplicazione di quel risultato per $w_j$.
  • Successivamente, la funzione di costo restituisce la somma di tutti i valori di aspettazione pesati.
• Utilizzare un ottimizzatore classico per determinare il prossimo insieme di parametri $\Theta_{i+1}$.
• Ripetere i passi precedenti fino al raggiungimento della convergenza.

Ricostruirai la funzione di costo di SSVQE nella valutazione, ma ecco il seguente frammento di codice per motivare la tua soluzione:

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)

Il VQD è un metodo iterativo che estende il VQE per ottenere i primi $k$ autovalori di un osservabile $\hat{H}$ con autovalori $\{\lambda_0, \lambda_1,...,\lambda_{N-1}\}$, dove $N\geq k$, invece del solo primo. Per il resto di questa sezione, assumeremo, senza perdita di generalità, che $\lambda_0\leq\lambda_1\leq...\leq\lambda_{N-1}$. Il VQD introduce la nozione di costo di penalità per guidare il processo di ottimizzazione.

Il VQD introduce un termine di penalità, denotato come $\beta$, per bilanciare il contributo di ciascun termine di sovrapposizione al costo. Questo termine di penalità serve a penalizzare il processo di ottimizzazione se l'ortogonalità non viene raggiunta. Imponiamo questo vincolo perché gli autostati di un osservabile, o di un operatore hermitiano, corrispondenti a diversi autovalori sono sempre mutuamente ortogonali, o possono essere resi tali nel caso di degenerazione o autovalori ripetuti. Pertanto, imponendo l'ortogonalità con l'autostato corrispondente a $\lambda_0$, stiamo effettivamente ottimizzando sul sottospazio che corrisponde al resto degli autovalori $\{\lambda_1, \lambda_2,..., \lambda_{N-1}\}$. Qui, $\lambda_1$ è il più basso autovalore tra gli autovalori rimanenti e, quindi, la soluzione ottima del nuovo problema può essere ottenuta usando il teorema variazionale.

L'idea generale alla base del VQD è quella di usare il VQE come di consueto per ottenere il più basso autovalore $\lambda_0 := C_0(\vec\theta^0) \equiv C_\text{VQE}(\vec\theta^0)$ insieme al corrispondente autostato (approssimato) $|\psi(\vec{\theta^0})\rangle$ per qualche vettore di parametri ottimale $\vec{\theta^0}$. Poi, per ottenere il prossimo autovalore $\lambda_1 > \lambda_0$, invece di minimizzare la funzione di costo $C_0(\vec{\theta}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$, si ottimizza:

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

Il valore positivo $\beta_0$ dovrebbe idealmente essere maggiore di $\lambda_1-\lambda_0$.

Questo introduce una nuova funzione di costo che può essere vista come un problema vincolato, in cui si minimizza $C_\text{VQE}(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ soggetto al vincolo che lo stato deve essere ortogonale al $|\psi(\vec{\theta^0})\rangle$ precedentemente ottenuto, con $\beta_0$ che agisce come termine di penalità se il vincolo non è soddisfatto.

In alternativa, questo nuovo problema può essere interpretato come l'esecuzione del VQE sul nuovo osservabile:

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

Assumendo che la soluzione al nuovo problema sia $|\psi(\vec{\theta^1})\rangle$, il valore atteso di $\hat{H}$ (non $\hat{H_1}$) dovrebbe essere $\langle \psi(\vec{\theta^1}) | \hat{H} | \psi(\vec{\theta^1})\rangle = \lambda_1$. Per ottenere il terzo autovalore $\lambda_2$, la funzione di costo da ottimizzare è:

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

dove $\beta_1$ è una costante positiva sufficientemente grande da imporre l'ortogonalità dello stato soluzione rispetto sia a $|\psi(\vec{\theta^0})\rangle$ che a $|\psi(\vec{\theta^1})\rangle$. Questo penalizza gli stati nello spazio di ricerca che non soddisfano questo requisito, restringendo efficacemente lo spazio di ricerca. Pertanto, la soluzione ottima del nuovo problema dovrebbe essere l'autostato corrispondente a $\lambda_2$.

Come nel caso precedente, questo nuovo problema può essere interpretato anche come un VQE con l'osservabile:

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

Se la soluzione a questo nuovo problema è $|\psi(\vec{\theta^2})\rangle$, il valore atteso di $\hat{H}$ (non $\hat{H_2}$) dovrebbe essere $\langle \psi(\vec{\theta^2}) | \hat{H} | \psi(\vec{\theta^2})\rangle = \lambda_2$. In modo analogo, per ottenere il $k$-esimo autovalore $\lambda_{k-1}$, si minimizza la funzione di costo:

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

Ricorda che abbiamo definito $\vec{\theta^j}$ tale che $\langle \psi(\vec{\theta^j}) | \hat{H} | \psi(\vec{\theta^j})\rangle = \lambda_j, \forall j<k$. Questo problema è equivalente a minimizzare $C(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ ma con il vincolo che lo stato deve essere ortogonale a $|\psi(\vec{\theta^j})\rangle ; \forall j \in {0, \cdots, k-1}$, restringendo così lo spazio di ricerca al sottospazio corrispondente agli autovalori $\{\lambda_{k-1},\cdots,\lambda_{N-1}\}$.

Questo problema è equivalente a un VQE con l'osservabile:

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

Come si può vedere dal processo, per ottenere il $k$-esimo autovalore, si hanno bisogno degli autostati (approssimati) dei precedenti $k-1$ autovalori, quindi sarebbe necessario eseguire il VQE un totale di $k$ volte. Pertanto, la funzione di costo del VQD è la seguente:

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

Schema teorico

La struttura del VQD può essere riassunta come segue:

• Preparare un operatore di riferimento $U_R$
• Applicare la forma variazionale $U_V(\vec\theta_{i,j})$ allo stato di riferimento, creando i seguenti ansatz $U_A(\vec\theta_{i,j})$
• Effettuare il bootstrap a $i=0$ se si dispone di un problema simile (tipicamente trovato tramite simulazione classica o campionamento).
• Valutare la funzione di costo $C_k(\vec{\theta})$, che comporta il calcolo di $k$ stati eccitati e un array di $\beta$ che definisce la penalità di sovrapposizione per ciascun termine.
  • Calcolare il valore di aspettazione per un osservabile $\langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ per ciascun $k$
  • Calcolare la penalità $\sum_{j=0}^{k-1}\beta_j |\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2$.
  • La funzione di costo deve quindi restituire la somma di questi due termini
• Utilizzare un ottimizzatore classico per scegliere il prossimo insieme di parametri $\Theta_{i+1}$.
• Ripetere questo processo fino al raggiungimento della convergenza.

Implementazione

Per questa implementazione, creeremo una funzione per una penalità di sovrapposizione. Questa penalità verrà usata nella funzione di costo ad ogni iterazione. Questo processo verrà ripetuto per ciascuno stato eccitato.

from qiskit.circuit.library import TwoLocal

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

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

Per prima cosa, definiremo una funzione che calcola la fedeltà degli stati — una percentuale di sovrapposizione tra due stati che useremo come penalità per il 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)

È il momento di scrivere la funzione di costo del VQD. Come prima quando calcolavamo solo lo stato fondamentale, determineremo lo stato a energia più bassa usando la primitiva Estimator. Tuttavia, come descritto sopra, aggiungeremo ora un termine di penalità per garantire l'ortogonalità degli stati a energia più alta. Ovvero, per ogni nuovo stato eccitato, viene aggiunta una penalità per qualsiasi sovrapposizione tra lo stato variazionale corrente e gli autostati a energia più bassa già trovati.

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

Nota in particolare che la funzione di costo sopra fa riferimento alla funzione calculate_overlaps, che in realtà crea un nuovo circuito quantistico. Se vogliamo eseguire su hardware reale, quel nuovo circuito deve anche essere traspilato, possibilmente in modo ottimale, per l'esecuzione sul backend scelto. Nota che la traspilazione non è stata integrata nelle funzioni calculate_overlaps o cost_func_vqd. Sentiti libero di provare a modificare il codice per integrare questa traspilazione aggiuntiva (e condizionale) — ma questo verrà fatto anche per te nella lezione successiva.

In questa lezione, eseguiremo l'algoritmo VQD usando lo Statevector Sampler e lo Statevector Estimator:

from qiskit.primitives import StatevectorEstimator as Estimator

sampler = Sampler()
estimator = Estimator()

Introdurremo un osservabile da stimare. Nella prossima lezione aggiungeremo un contesto fisico, come lo stato eccitato di una molecola. Può essere utile pensare a questo osservabile come all'Hamiltoniano di un sistema che può avere stati eccitati, anche se questo osservabile non è stato scelto per corrispondere a nessuna molecola o atomo in particolare.

from qiskit.quantum_info import SparsePauliOp

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

Qui, impostiamo il numero totale di stati che vogliamo calcolare (stato fondamentale e stati eccitati, k) e le penalità (betas) per la sovrapposizione tra statevector che dovrebbero essere ortogonali. Le conseguenze della scelta di betas troppo alti o troppo bassi saranno esplorate nella prossima lezione. Per ora, useremo semplicemente quelli forniti di seguito. Inizieremo usando tutti zeri come parametri. Nei tuoi calcoli potresti voler usare parametri iniziali più intelligenti basati sulla tua conoscenza del sistema o su calcoli precedenti.

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

Possiamo ora eseguire il calcolo:

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

I valori ottenuti dalla funzione di costo sono approssimativamente -6,00, 4,02 e 5,61. La cosa importante di questi risultati è che i valori della funzione sono crescenti. Se avessimo ottenuto un primo stato eccitato con energia inferiore al nostro calcolo iniziale e non vincolato dello stato fondamentale, ciò avrebbe indicato un errore da qualche parte nel codice.

I valori di x sono i parametri che hanno prodotto uno statevector corrispondente a ciascuno di questi costi (energie).

Infine, notiamo che tutte e tre le minimizzazioni sono convergite entro la tolleranza predefinita dell'ottimizzatore classico (qui COBYLA). Hanno richiesto rispettivamente 131, 110 e 90 valutazioni della funzione.

Quantum Sampling Regression (QSR)

Uno dei principali problemi del VQE è il numero elevato di chiamate a un computer quantistico necessarie per ottenere i parametri a ogni passo, ad esempio $k$, $k-1$, e così via. Ciò è particolarmente problematico quando l'accesso ai dispositivi quantistici è in coda. Sebbene una Session possa essere usata per raggruppare più chiamate iterative, un approccio alternativo è quello di usare il campionamento. Sfruttando più risorse classiche, è possibile completare l'intero processo di ottimizzazione con una singola chiamata. È qui che entra in gioco la Quantum Sampling Regression. Poiché l'accesso ai computer quantistici è ancora una risorsa con bassa offerta e alta domanda, questo compromesso risulta sia praticabile che conveniente per molti studi attuali. Questo approccio sfrutta tutte le capacità classiche disponibili pur catturando molte delle funzionalità interne e delle proprietà intrinseche dei calcoli quantistici che non emergono nella simulazione.

L'idea alla base del QSR è che la funzione di costo $C(\theta) := \langle \psi(\theta) | \hat{H} | \psi(\theta)\rangle$ può essere espressa come una serie di Fourier nel seguente modo:

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

A seconda della periodicità e della larghezza di banda della funzione originale, l'insieme $S$ può essere finito o infinito. Ai fini di questa discussione, assumeremo che sia infinito. Il passo successivo è campionare la funzione di costo $C(\theta)$ più volte per ottenere i coefficienti di Fourier $\{a_0, a_k, b_k\}_{k=1}^S$. In particolare, poiché abbiamo $2S+1$ incognite, sarà necessario campionare la funzione di costo $2S+1$ volte.

Se si campiona la funzione di costo per $2S+1$ valori di parametri $\{\theta_1,...,\theta_{2S+1}\}$, si ottiene il seguente sistema:

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

che riscriviamo come

$$Fa=c.$$

In pratica, questo sistema generalmente non è consistente perché i valori della funzione di costo $c$ non sono esatti. Pertanto, di solito è una buona idea normalizzarli moltiplicando per $F^\dagger$ a sinistra, il che porta a:

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

Questo nuovo sistema è sempre consistente, e la sua soluzione è una soluzione ai minimi quadrati del problema originale. Se abbiamo $k$ parametri invece di uno solo, e ciascun parametro $\theta^i$ ha il proprio $S_i$ per $i \in {1,...,k}$, allora il numero totale di campioni richiesti è:

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

dove $S_{\max} = \max_i(S_i)$. Inoltre, regolare $S_{\max}$ come parametro sintonizzabile (invece di inferirlo) apre nuove possibilità, come:

• Sovra-campionamento per migliorare la precisione.
• Sotto-campionamento per migliorare le prestazioni riducendo il sovraccarico di runtime o eliminando i minimi locali.

Schema teorico

La struttura del QSR può essere riassunta come segue:

• Preparare gli operatori di riferimento $U_R$.
  • Si passa dallo stato $|0\rangle$ allo stato di riferimento $|\rho\rangle$
• Applicare la forma variazionale $U_V(\vec\theta_{i,j})$ per creare un ansatz $U_A(\vec\theta_{i,j})$.
  • Determinare la larghezza di banda associata a ciascun parametro nell'ansatz. È sufficiente un limite superiore.
• Effettuare il bootstrap a $i=0$ se si dispone di un problema simile (tipicamente trovato tramite simulazione classica o campionamento).
• Campionare la funzione di costo $C(\vec\theta) := a_0 + \sum_{k=1}^S[a_k\cos(k\theta)+ b_k\sin(k\theta)]$ almeno $T$ volte.
  • $T=\prod_{i=1}^k(2S_i+1)\leq \prod_{i=1}^k(2S_{max}+1) = (2S_{max}+1)^n$
  • Decidere se effettuare sovra-campionamento/sotto-campionamento per bilanciare velocità e precisione regolando $T$.
• Calcolare i coefficienti di Fourier dai campioni (ovvero, risolvere il sistema lineare normalizzato di equazioni).
• Trovare il minimo globale della funzione di regressione risultante su una macchina classica.

Riepilogo

Con questa lezione, hai appreso diversi algoritmi variazionali disponibili:

• Schema generale
• Introduzione di pesi e penalità per modificare una funzione di costo
• Esplorazione del sotto-campionamento vs sovra-campionamento per bilanciare velocità e precisione

Queste idee possono essere adattate per formare un algoritmo variazionale personalizzato adatto al tuo problema. Ti incoraggiamo a condividere i tuoi risultati con la comunità. La prossima lezione esplorerà come usare un algoritmo variazionale per risolvere un'applicazione.