Instances et extensions

Ce chapitre couvre plusieurs algorithmes quantiques variationnels, notamment

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

En utilisant ces algorithmes, tu apprendras plusieurs idées de conception pouvant être intégrées dans des algorithmes variationnels personnalisés, comme les poids, les pénalités, le sur-échantillonnage et le sous-échantillonnage. Nous t'encourageons à expérimenter ces concepts et à partager tes découvertes avec la communauté.

Le framework Qiskit patterns s'applique à tous ces algorithmes — mais nous ne mentionnerons explicitement les étapes que dans le premier exemple.

Variational Quantum Eigensolver (VQE)

VQE est l'un des algorithmes quantiques variationnels les plus utilisés, établissant un modèle sur lequel d'autres algorithmes peuvent s'appuyer.

Étape 1 : Mapper les entrées classiques vers un problème quantique

Organisation théorique

La structure de VQE est simple :

• Préparer les opérateurs de référence $U_R$
  • On part de l'état $|0\rangle$ et on atteint l'état de référence $|\rho\rangle$
• Appliquer la forme variationnelle $U_V(\vec\theta_{i,j})$ pour créer un ansatz $U_A(\vec\theta_{i,j})$
  • On passe de l'état $|\rho\rangle$ à $U_V(\vec\theta_{i,j})|\rho\rangle = |\psi(\vec\theta_{i,j})\rangle$
• Amorcer à $i=0$ si on dispose d'un problème similaire (généralement trouvé via simulation classique ou échantillonnage)
  • Chaque optimiseur sera amorcé différemment, produisant un ensemble initial de vecteurs de paramètres $\Theta_0 := \\{ {\vec\theta_{0,j} | j \in \mathcal{J}_\text{opt}^0} \\}$ (par exemple, à partir d'un point initial $\vec\theta_0$).
• Évaluer la fonction de coût $C(\vec\theta_{i,j}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ pour tous les états préparés sur un ordinateur quantique.
• Utiliser un optimiseur classique pour sélectionner le prochain ensemble de paramètres $\Theta_{i+1}$.
• Répéter le processus jusqu'à convergence.

Il s'agit d'une simple boucle d'optimisation classique dans laquelle on évalue la fonction de coût. Certains optimiseurs peuvent nécessiter plusieurs évaluations pour calculer un gradient, déterminer la prochaine itération ou évaluer la convergence.

Voici l'exemple pour l'observable suivante :

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

Implémentation

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

On peut utiliser cette fonction de coût pour calculer les paramètres optimaux

# 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

Étape 2 : Optimiser le problème pour l'exécution quantique

Nous allons sélectionner le backend le moins occupé et importer les composants nécessaires depuis 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')>

Nous allons transpiler le circuit en utilisant le gestionnaire de passes préconfiguré avec le niveau d'optimisation 3, et appliquer la disposition correspondante à l'observable.

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)

Étape 3 : Exécuter avec les primitives Qiskit Runtime

Nous sommes maintenant prêts à lancer notre calcul sur le matériel IBM Quantum®. Comme la minimisation de la fonction de coût est très itérative, nous allons démarrer une session Runtime. De cette façon, nous n'aurons à attendre dans la file d'attente qu'une seule fois. Une fois la tâche en cours d'exécution, chaque itération avec les mises à jour des paramètres s'exécutera immédiatement.

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)

Étape 4 : Post-traitement, retourner le résultat en format classique

On peut voir que la routine de minimisation s'est terminée avec succès, ce qui signifie que nous avons atteint la tolérance par défaut de l'optimiseur classique COBYLA. Si un résultat plus précis est requis, on peut spécifier une tolérance plus petite. Cela peut effectivement être nécessaire, car le résultat était décalé de plusieurs pourcents par rapport au résultat obtenu par le simulateur ci-dessus.

La valeur de x obtenue est la meilleure estimation actuelle des paramètres qui minimisent la fonction de coût. En cas d'itération pour obtenir une précision plus élevée, ces valeurs devraient être utilisées à la place du x0 initialement utilisé (un vecteur de uns).

Enfin, notons que la fonction a été évaluée 96 fois au cours du processus d'optimisation. Ce nombre peut différer du nombre d'étapes d'optimisation, car certains optimiseurs nécessitent plusieurs évaluations de la fonction en une seule étape, par exemple lors de l'estimation d'un gradient.

Subspace Search VQE (SSVQE)

SSVQE est une variante de VQE qui permet d'obtenir les $k$ premières valeurs propres d'une observable $\hat{H}$ avec les valeurs propres $\{\lambda_0, \lambda_1,...,\lambda_{N-1}\}$, où $N\geq k$. Sans perte de généralité, on suppose que $\lambda_0<\lambda_1<...<\lambda_{N-1}$. SSVQE introduit une nouvelle idée en ajoutant des poids pour aider à prioriser l'optimisation du terme ayant le poids le plus élevé.

Pour implémenter cet algorithme, nous avons besoin de $k$ états de référence mutuellement orthogonaux $\{ |\rho_j\rangle \}_{j=0}^{k-1}$, c'est-à-dire $\langle \rho_j | \rho_l \rangle = \delta_{jl}$ pour $j,l<k$. Ces états peuvent être construits à partir d'opérateurs de Pauli. La fonction de coût de cet algorithme est alors :

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

où $w_j$ est un nombre positif arbitraire tel que si $j<l<k$ alors $w_j>w_l$, et $U_V(\vec{\theta})$ est la forme variationnelle définie par l'utilisateur.

L'algorithme SSVQE repose sur le fait que les états propres correspondant à des valeurs propres différentes sont mutuellement orthogonaux. Plus précisément, le produit scalaire de $U_V(\vec{\theta})|\rho_j\rangle$ et $U_V(\vec{\theta})|\rho_l\rangle$ peut s'exprimer comme suit :

$$\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 première égalité tient parce que $U_{V}(\vec{\theta})$ est un opérateur quantique et est donc unitaire. La dernière égalité tient en raison de l'orthogonalité des états de référence $|\rho_j\rangle$. Le fait que l'orthogonalité soit préservée par les transformations unitaires est profondément lié au principe de conservation de l'information, tel qu'il est exprimé en science de l'information quantique. Sous cet angle, les transformations non unitaires représentent des processus où l'information est soit perdue, soit injectée.

Les poids $w_j$ aident à garantir que tous les états sont des états propres. Si les poids sont suffisamment différents, le terme ayant le plus grand poids (c'est-à-dire $w_0$) sera prioritaire lors de l'optimisation par rapport aux autres. Par conséquent, l'état résultant $U_{V}(\vec{\theta})|\rho_0 \rangle$ deviendra l'état propre correspondant à $\lambda_0$. Comme $\{ U_{V}(\vec{\theta})|\rho_j\rangle \}_{j=0}^{k-1}$ sont mutuellement orthogonaux, les états restants lui seront orthogonaux et, par conséquent, contenus dans le sous-espace correspondant aux valeurs propres $\{\lambda_1,...,\lambda_{N-1}\}$.

En appliquant le même raisonnement au reste des termes, la prochaine priorité serait alors le terme de poids $w_1$, de sorte que $U_{V}(\vec{\theta})|\rho_1 \rangle$ serait l'état propre correspondant à $\lambda_1$, et les autres termes seraient contenus dans l'espace propre de $\{\lambda_2,...,\lambda_{N-1}\}$.

En raisonnant par induction, on déduit que $U_{V}(\vec{\theta})|\rho_j \rangle$ sera un état propre approché de $\lambda_j$ pour $0\leq j < k.$

Organisation théorique

SSVQE peut se résumer comme suit :

• Préparer plusieurs états de référence en appliquant un unitaire U_R à k états de base computationnels différents
  • Cet algorithme nécessite l'utilisation de $k$ états de référence mutuellement orthogonaux $\{ |\rho_j\rangle \}_{j=0}^{k-1}$, tels que $\langle \rho_j | \rho_l \rangle = \delta_{jl}$ pour $j,l<k$.
• Appliquer la forme variationnelle $U_V(\vec\theta_{i,j})$ à chaque état de référence, donnant l'ansatz $U_A(\vec\theta_{i,j})$.
• Amorcer à $i=0$ si un problème similaire est disponible (généralement trouvé via simulation classique ou échantillonnage).
• Évaluer la fonction de coût $C(\vec\theta_{i,j}) := \sum_{j=0}^{k-1} w_j \langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ pour tous les états préparés sur un ordinateur quantique.
  • Cela peut se décomposer en calculant la valeur d'espérance d'une observable $\langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ et en multipliant ce résultat par $w_j$.
  • Ensuite, la fonction de coût retourne la somme de toutes les valeurs d'espérance pondérées.
• Utiliser un optimiseur classique pour déterminer le prochain ensemble de paramètres $\Theta_{i+1}$.
• Répéter les étapes ci-dessus jusqu'à convergence.

Tu vas reconstruire la fonction de coût de SSVQE dans l'évaluation, mais nous avons l'extrait suivant pour t'aider à élaborer ta solution :

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

Déflexion Quantique Variationnelle (VQD)

VQD est une méthode itérative qui étend VQE pour obtenir les $k$ premières valeurs propres d'un observable $\hat{H}$ aux valeurs propres $\{\lambda_0, \lambda_1,...,\lambda_{N-1}\}$, où $N\geq k$, et non plus seulement la première. Pour la suite de cette section, on supposera, sans perte de généralité, que $\lambda_0\leq\lambda_1\leq...\leq\lambda_{N-1}$. VQD introduit la notion de coût de pénalité pour guider le processus d'optimisation.

VQD introduit un terme de pénalité, noté $\beta$, pour équilibrer la contribution de chaque terme de chevauchement au coût. Ce terme de pénalité sert à pénaliser le processus d'optimisation si l'orthogonalité n'est pas atteinte. On impose cette contrainte car les états propres d'un observable, ou d'un opérateur hermitien, correspondant à des valeurs propres différentes sont toujours mutuellement orthogonaux, ou peuvent l'être en cas de dégénérescence ou de valeurs propres répétées. Ainsi, en imposant l'orthogonalité avec l'état propre correspondant à $\lambda_0$, on optimise effectivement sur le sous-espace correspondant aux valeurs propres restantes $\{\lambda_1, \lambda_2,..., \lambda_{N-1}\}$. Ici, $\lambda_1$ est la valeur propre la plus basse parmi les valeurs propres restantes et, par conséquent, la solution optimale du nouveau problème peut être obtenue à l'aide du théorème variationnel.

L'idée générale derrière VQD est d'utiliser VQE comme d'habitude pour obtenir la valeur propre la plus basse $\lambda_0 := C_0(\vec\theta^0) \equiv C_\text{VQE}(\vec\theta^0)$ ainsi que l'état propre (approximatif) correspondant $|\psi(\vec{\theta^0})\rangle$ pour un vecteur de paramètres optimal $\vec{\theta^0}$. Puis, pour obtenir la valeur propre suivante $\lambda_1 > \lambda_0$, au lieu de minimiser la fonction de coût $C_0(\vec{\theta}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$, on optimise :

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

La valeur positive $\beta_0$ devrait idéalement être supérieure à $\lambda_1-\lambda_0$.

Ceci introduit une nouvelle fonction de coût qui peut être vue comme un problème contraint, où on minimise $C_\text{VQE}(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ sous la contrainte que l'état doit être orthogonal au $|\psi(\vec{\theta^0})\rangle$ précédemment obtenu, avec $\beta_0$ agissant comme terme de pénalité si la contrainte n'est pas satisfaite.

Autrement dit, ce nouveau problème peut être interprété comme l'exécution de VQE sur le nouvel observable :

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

En supposant que la solution du nouveau problème est $|\psi(\vec{\theta^1})\rangle$, la valeur d'espérance de $\hat{H}$ (pas $\hat{H_1}$) devrait être $\langle \psi(\vec{\theta^1}) | \hat{H} | \psi(\vec{\theta^1})\rangle = \lambda_1$. Pour obtenir la troisième valeur propre $\lambda_2$, la fonction de coût à optimiser est :

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

où $\beta_1$ est une constante positive suffisamment grande pour imposer l'orthogonalité de l'état solution à la fois à $|\psi(\vec{\theta^0})\rangle$ et à $|\psi(\vec{\theta^1})\rangle$. Cela pénalise les états de l'espace de recherche qui ne satisfont pas cette exigence, ce qui restreint efficacement l'espace de recherche. Ainsi, la solution optimale du nouveau problème devrait être l'état propre correspondant à $\lambda_2$.

Comme dans le cas précédent, ce nouveau problème peut également être interprété comme VQE avec l'observable :

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

Si la solution de ce nouveau problème est $|\psi(\vec{\theta^2})\rangle$, la valeur d'espérance de $\hat{H}$ (pas $\hat{H_2}$) devrait être $\langle \psi(\vec{\theta^2}) | \hat{H} | \psi(\vec{\theta^2})\rangle = \lambda_2$. De façon analogue, pour obtenir la $k$-ième valeur propre $\lambda_{k-1}$, tu minimiserais la fonction de coût :

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

Rappelons que nous avons défini $\vec{\theta^j}$ tel que $\langle \psi(\vec{\theta^j}) | \hat{H} | \psi(\vec{\theta^j})\rangle = \lambda_j, \forall j<k$. Ce problème est équivalent à minimiser $C(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ mais avec la contrainte que l'état doit être orthogonal à $|\psi(\vec{\theta^j})\rangle ; \forall j \in {0, \cdots, k-1}$, restreignant ainsi l'espace de recherche au sous-espace correspondant aux valeurs propres $\{\lambda_{k-1},\cdots,\lambda_{N-1}\}$.

Ce problème est équivalent à un VQE avec l'observable :

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

Comme tu peux le voir dans ce processus, pour obtenir la $k$-ième valeur propre, tu as besoin des états propres (approximatifs) des $k-1$ valeurs propres précédentes, tu devras donc exécuter VQE un total de $k$ fois. Par conséquent, la fonction de coût de VQD est la suivante :

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

Structure théorique

La structure de VQD peut être résumée comme suit :

• Préparer un opérateur de référence $U_R$
• Appliquer la forme variationnelle $U_V(\vec\theta_{i,j})$ à l'état de référence, en créant les ansatze suivants $U_A(\vec\theta_{i,j})$
• Initialiser à $i=0$ si on a un problème similaire (généralement trouvé par simulation classique ou échantillonnage).
• Évaluer la fonction de coût $C_k(\vec{\theta})$, ce qui implique de calculer $k$ états excités et un tableau de $\beta$ définissant la pénalité de chevauchement pour chaque terme de chevauchement.
  • Calculer la valeur d'espérance pour un observable $\langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ pour chaque $k$
  • Calculer la pénalité $\sum_{j=0}^{k-1}\beta_j |\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2$.
  • La fonction de coût doit alors retourner la somme de ces deux termes
• Utiliser un optimiseur classique pour choisir le prochain ensemble de paramètres $\Theta_{i+1}$.
• Répéter ce processus jusqu'à convergence.

Implémentation

Pour cette implémentation, nous allons créer une fonction pour une pénalité de chevauchement. Cette pénalité sera utilisée dans la fonction de coût à chaque itération. Ce processus sera répété pour chaque état excité.

from qiskit.circuit.library import TwoLocal

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

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

Tout d'abord, nous allons définir une fonction qui calcule la fidélité d'état — un pourcentage de chevauchement entre deux états que nous utiliserons comme pénalité pour 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 est temps d'écrire la fonction de coût de VQD. Comme précédemment, lorsque nous calculions uniquement l'état fondamental, nous allons déterminer l'état d'énergie la plus basse en utilisant la primitive Estimator. Cependant, comme décrit ci-dessus, nous allons maintenant ajouter un terme de pénalité pour assurer l'orthogonalité des états de plus haute énergie. C'est-à-dire que pour chaque nouvel état excité, une pénalité est ajoutée pour tout chevauchement entre l'état variationnel courant et les états propres de plus basse énergie déjà trouvés.

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

Note bien que la fonction de coût ci-dessus fait référence à la fonction calculate_overlaps, qui crée en réalité un nouveau circuit quantique. Si tu veux exécuter sur du vrai matériel, ce nouveau circuit doit également être transpilé, idéalement de façon optimale, pour fonctionner sur le backend choisi. Note que la transpilation n'a pas été intégrée aux fonctions calculate_overlaps ou cost_func_vqd. N'hésite pas à modifier toi-même le code pour y intégrer cette transpilation supplémentaire (et conditionnelle) — mais cela sera aussi fait pour toi dans la prochaine leçon.

Dans cette leçon, nous allons exécuter l'algorithme VQD en utilisant le Statevector Sampler et le Statevector Estimator :

from qiskit.primitives import StatevectorEstimator as Estimator

sampler = Sampler()
estimator = Estimator()

Nous allons introduire un observable à estimer. Dans la prochaine leçon, nous ajouterons un contexte physique, comme l'état excité d'une molécule. Il peut être utile de penser à cet observable comme au Hamiltonien d'un système pouvant avoir des états excités, même si cet observable n'a pas été choisi pour correspondre à une molécule ou un atome particulier.

from qiskit.quantum_info import SparsePauliOp

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

Ici, nous fixons le nombre total d'états que nous souhaitons calculer (état fondamental et états excités, k), ainsi que les pénalités (betas) pour le chevauchement entre vecteurs d'état qui devraient être orthogonaux. Les conséquences de choisir des betas trop élevés ou trop bas seront un peu explorées dans la prochaine leçon. Pour l'instant, nous utiliserons simplement ceux fournis ci-dessous. Nous commencerons en utilisant des zéros comme paramètres. Dans tes propres calculs, tu voudras peut-être utiliser des paramètres de départ plus judicieux basés sur ta connaissance du système ou sur des calculs antérieurs.

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

Nous pouvons maintenant lancer le calcul :

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

Les valeurs que nous avons obtenues de la fonction de coût sont approximativement -6,00, 4,02 et 5,61. Ce qui importe dans ces résultats, c'est que les valeurs de la fonction sont croissantes. Si nous avions obtenu un premier état excité d'énergie inférieure à notre calcul initial non contraint de l'état fondamental, cela aurait indiqué une erreur quelque part dans notre code.

Les valeurs de x sont les paramètres qui ont produit un vecteur d'état correspondant à chacun de ces coûts (énergies).

Enfin, nous notons que les trois minimisations ont convergé dans la tolérance par défaut de l'optimiseur classique (ici COBYLA). Elles ont nécessité respectivement 131, 110 et 90 évaluations de la fonction.

Régression par Échantillonnage Quantique (QSR)

L'un des principaux problèmes de VQE est le grand nombre d'appels à un ordinateur quantique requis pour obtenir les paramètres à chaque étape, par exemple, $k$, $k-1$, etc. C'est particulièrement problématique lorsque l'accès aux dispositifs quantiques est soumis à une file d'attente. Bien qu'une Session puisse être utilisée pour regrouper plusieurs appels itératifs, une approche alternative consiste à utiliser l'échantillonnage. En exploitant davantage de ressources classiques, on peut compléter l'intégralité du processus d'optimisation en un seul appel. C'est là qu'intervient la Régression par Échantillonnage Quantique. Puisque l'accès aux ordinateurs quantiques reste une ressource rare très demandée, ce compromis nous semble à la fois possible et avantageux pour de nombreuses études actuelles. Cette approche exploite toutes les capacités classiques disponibles tout en capturant une grande partie du fonctionnement interne et des propriétés intrinsèques des calculs quantiques qui n'apparaissent pas dans les simulations.

L'idée derrière QSR est que la fonction de coût $C(\theta) := \langle \psi(\theta) | \hat{H} | \psi(\theta)\rangle$ peut être exprimée comme une série de Fourier de la manière suivante :

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

Selon la périodicité et la largeur de bande de la fonction originale, l'ensemble $S$ peut être fini ou infini. Pour les besoins de cette discussion, nous supposerons qu'il est infini. L'étape suivante consiste à échantillonner la fonction de coût $C(\theta)$ plusieurs fois afin d'obtenir les coefficients de Fourier $\{a_0, a_k, b_k\}_{k=1}^S$. Plus précisément, puisque nous avons $2S+1$ inconnues, nous devrons échantillonner la fonction de coût $2S+1$ fois.

Si nous échantillonnons ensuite la fonction de coût pour $2S+1$ valeurs de paramètres $\{\theta_1,...,\theta_{2S+1}\}$, nous pouvons obtenir le système suivant :

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

que nous réécrirons comme

$$Fa=c.$$

En pratique, ce système n'est généralement pas consistant car les valeurs de la fonction de coût $c$ ne sont pas exactes. Il est donc généralement judicieux de les normaliser en les multipliant par $F^\dagger$ à gauche, ce qui donne :

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

Ce nouveau système est toujours consistant, et sa solution est une solution aux moindres carrés du problème original. Si nous avons $k$ paramètres au lieu d'un seul, et que chaque paramètre $\theta^i$ a son propre $S_i$ pour $i \in {1,...,k}$, alors le nombre total d'échantillons requis est :

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

où $S_{\max} = \max_i(S_i)$. De plus, ajuster $S_{\max}$ en tant que paramètre réglable (plutôt que de l'inférer) ouvre de nouvelles possibilités, telles que :

• Le sur-échantillonnage pour améliorer la précision.
• Le sous-échantillonnage pour améliorer les performances en réduisant la surcharge d'exécution ou en éliminant les minima locaux.

Structure théorique

La structure de QSR peut être résumée comme suit :

• Préparer les opérateurs de référence $U_R$.
  • On passera de l'état $|0\rangle$ à l'état de référence $|\rho\rangle$
• Appliquer la forme variationnelle $U_V(\vec\theta_{i,j})$ pour créer un ansatz $U_A(\vec\theta_{i,j})$.
  • Déterminer la largeur de bande associée à chaque paramètre dans l'ansatz. Une borne supérieure suffit.
• Initialiser à $i=0$ si on a un problème similaire (généralement trouvé par simulation classique ou échantillonnage).
• Échantillonner la fonction de coût $C(\vec\theta) := a_0 + \sum_{k=1}^S[a_k\cos(k\theta)+ b_k\sin(k\theta)]$ au moins $T$ fois.
  • $T=\prod_{i=1}^k(2S_i+1)\leq \prod_{i=1}^k(2S_{max}+1) = (2S_{max}+1)^n$
  • Décider de sur-échantillonner ou sous-échantillonner pour équilibrer vitesse et précision en ajustant $T$.
• Calculer les coefficients de Fourier à partir des échantillons (c'est-à-dire, résoudre le système d'équations linéaires normalisé).
• Trouver le minimum global de la fonction de régression résultante sur une machine classique.

Résumé

Avec cette leçon, tu as découvert plusieurs instances variationnelles disponibles :

• Structure générale
• Introduction de poids et de pénalités pour ajuster une fonction de coût
• Explorer le sous-échantillonnage et le sur-échantillonnage pour trouver un compromis entre vitesse et précision

Ces idées peuvent être adaptées pour former un algorithme variationnel personnalisé correspondant à ton problème. Nous t'encourageons à partager tes résultats avec la communauté. La prochaine leçon explorera comment utiliser un algorithme variationnel pour résoudre une application.