Exemples et applications

Dans cette leçon, nous allons explorer quelques exemples d'algorithmes variationnels et la façon de les appliquer :

• Comment écrire un algorithme variationnel personnalisé
• Comment appliquer un algorithme variationnel pour trouver les valeurs propres minimales
• Comment utiliser les algorithmes variationnels pour résoudre des cas d'usage applicatifs

Le cadre des patterns Qiskit peut être appliqué à tous les problèmes présentés ici. Cependant, afin d'éviter les répétitions, nous ne mentionnerons explicitement les étapes du cadre que dans un seul exemple, exécuté sur du matériel réel.

Définitions des problèmes

Imaginons que nous souhaitons utiliser un algorithme variationnel pour trouver la valeur propre de l'observable suivante :

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

Cet observable a les valeurs propres suivantes :

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

Et les états propres suivants :

$$\left\{ \begin{array}{c} |\phi_0\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)\\ |\phi_1\rangle = \frac{1}{\sqrt{2}}(|00\rangle - |11\rangle)\\ |\phi_2\rangle = \frac{1}{\sqrt{2}}(|01\rangle - |10\rangle)\\ |\phi_3\rangle = \frac{1}{\sqrt{2}}(|01\rangle + |10\rangle) \end{array} \right\}$$

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

from qiskit.quantum_info import SparsePauliOp

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

VQE personnalisé

Nous allons d'abord explorer comment construire manuellement une instance de VQE pour trouver la valeur propre la plus basse de $\hat{O}_1$. Cela incorporera une variété de techniques que nous avons abordées tout au long de ce cours.

def cost_func_vqe(params, ansatz, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """
    pub = (ansatz, hamiltonian, params)
    cost = estimator.run([pub]).result()[0].data.evs

    return cost

from qiskit.circuit.library.n_local import n_local
from qiskit import QuantumCircuit

import numpy as np

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

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

raw_ansatz = reference_circuit.compose(variational_form)
raw_ansatz.decompose().draw("mpl")

Nous allons commencer par déboguer sur des simulateurs locaux.

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

estimator = Estimator()
sampler = Sampler()

Nous définissons maintenant un ensemble de paramètres initiaux :

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

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

Nous pouvons minimiser cette fonction de coût pour calculer les paramètres optimaux

# SciPy minimizer routine
from scipy.optimize import minimize
import time

start_time = time.time()

result = minimize(
    cost_func_vqe,
    x0,
    args=(raw_ansatz, observable_1, estimator),
    method="COBYLA",
    options={"maxiter": 1000, "disp": True},
)

end_time = time.time()
execution_time = end_time - start_time

Return from COBYLA because the trust region radius reaches its lower bound.
Number of function values = 103   Least value of F = -5.999999998357189
The corresponding X is:
[2.27483579e+00 8.37593091e-01 1.57080508e+00 5.82932911e-06
 2.49973063e+00 6.41884255e-01 6.33686904e-01 6.33688223e-01]

result

message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: -5.999999998357189
       x: [ 2.275e+00  8.376e-01  1.571e+00  5.829e-06  2.500e+00
            6.419e-01  6.337e-01  6.337e-01]
    nfev: 103
   maxcv: 0.0

Ce problème jouet n'utilisant que deux qubits, nous pouvons le vérifier en utilisant le solveur d'algèbre linéaire de NumPy.

from numpy.linalg import eigvalsh

solution_eigenvalue = min(eigvalsh(observable_1.to_matrix()))

print(f"""Number of iterations: {result.nfev}""")
print(f"""Time (s): {execution_time}""")

print(
    f"Percent error: {100*abs((result.fun - solution_eigenvalue)/solution_eigenvalue):.2e}"
)

Number of iterations: 103
Time (s): 0.4394676685333252
Percent error: 2.74e-08

Comme vous pouvez le voir, le résultat est extrêmement proche de l'idéal.

Expérimenter pour améliorer la vitesse et la précision

Ajouter un état de référence

Dans l'exemple précédent, nous n'avons utilisé aucun opérateur de référence $U_R$. Réfléchissons maintenant à la façon dont l'état propre idéal $\frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)$ peut être obtenu. Considérons le circuit suivant.

from qiskit import QuantumCircuit

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

ideal_qc.draw("mpl")

Nous pouvons rapidement vérifier que ce circuit nous donne l'état désiré.

from qiskit.quantum_info import Statevector

Statevector(ideal_qc)

Statevector([0.70710678+0.j, 0.        +0.j, 0.        +0.j,
             0.70710678+0.j],
            dims=(2, 2))

Maintenant que nous avons vu à quoi ressemble un circuit préparant l'état solution, il semble raisonnable d'utiliser une porte de Hadamard comme circuit de référence, de sorte que l'ansatz complet devienne :

reference = QuantumCircuit(2)
reference.h(0)
reference.cx(0, 1)
# Include barrier to separate reference from variational form
reference.barrier()

ref_ansatz = variational_form.decompose().compose(reference, front=True)

ref_ansatz.draw("mpl")

Pour ce nouveau circuit, la solution idéale pourrait être atteinte avec tous les paramètres réglés à $0$, ce qui confirme que le choix du circuit de référence est raisonnable.

Comparons maintenant le nombre d'évaluations de la fonction de coût, les itérations de l'optimiseur et le temps passé avec ceux de la tentative précédente.

import time

start_time = time.time()

ref_result = minimize(
    cost_func_vqe, x0, args=(ref_ansatz, observable_1, estimator), method="COBYLA"
)

end_time = time.time()
execution_time = end_time - start_time

En utilisant nos paramètres optimaux pour calculer la valeur propre minimale :

experimental_min_eigenvalue_ref = cost_func_vqe(
    ref_result.x, ref_ansatz, observable_1, estimator
)
print(experimental_min_eigenvalue_ref)

-5.999999996759607

print("ADDED REFERENCE STATE:")
print(f"""Number of iterations: {ref_result.nfev}""")
print(f"""Time (s): {execution_time}""")
print(
    f"Percent error: {100*abs((experimental_min_eigenvalue_ref - solution_eigenvalue)/solution_eigenvalue):.2e}"
)

ADDED REFERENCE STATE:
Number of iterations: 127
Time (s): 0.5620882511138916
Percent error: 5.40e-08

Selon ton système spécifique, cela peut ou non entraîner une amélioration de la vitesse ou de la précision dans cet exemple à très petite échelle. L'essentiel est que le fait de partir d'états de référence physiquement motivés devient de plus en plus important pour améliorer la vitesse et la précision à mesure que les problèmes prennent de l'ampleur.

Changer le point initial

Maintenant que nous avons vu l'effet de l'ajout de l'état de référence, voyons ce qui se passe lorsque nous choisissons différents points initiaux $\vec{\theta_0}$. En particulier, nous utiliserons $\vec{\theta_0}=(0,0,0,0,6,0,0,0)$ et $\vec{\theta_0}=(6,6,6,6,6,6,6,6,6)$.

Rappelons que, comme nous l'avons évoqué lors de l'introduction de l'état de référence, la solution idéale serait trouvée lorsque tous les paramètres sont à $0$, donc le premier point initial devrait donner moins d'évaluations.

import time

start_time = time.time()

x0 = [0, 0, 0, 0, 6, 0, 0, 0]

x0_1_result = minimize(
    cost_func_vqe, x0, args=(raw_ansatz, observable_1, estimator), method="COBYLA"
)

end_time = time.time()
execution_time = end_time - start_time

print("INITIAL POINT 1:")
print(f"""Number of iterations: {x0_1_result.nfev}""")
print(f"""Time (s): {execution_time}""")

INITIAL POINT 1:
Number of iterations: 108
Time (s): 0.4492197036743164

Ajustement du point initial à $\vec{\theta_0}=(6,6,6,6,6,6,6,6,6)$ :

import time

start_time = time.time()

x0 = 6 * np.ones(raw_ansatz.num_parameters)

x0_2_result = minimize(
    cost_func_vqe, x0, args=(raw_ansatz, observable_1, estimator), method="COBYLA"
)

end_time = time.time()
execution_time = end_time - start_time

print("INITIAL POINT 2:")
print(f"""Number of iterations: {x0_2_result.nfev}""")
print(f"""Time (s): {execution_time}""")

INITIAL POINT 2:
Number of iterations: 107
Time (s): 0.40889453887939453

En expérimentant avec différents points initiaux, tu pourras peut-être obtenir une convergence plus rapide et avec moins d'évaluations de la fonction.

Expérimenter avec différents optimiseurs

Nous pouvons ajuster l'optimiseur en utilisant l'argument method de minimize de SciPy, avec plus d'options disponibles ici. Nous avons initialement utilisé un minimiseur sous contraintes (COBYLA). Dans cet exemple, nous allons explorer l'utilisation d'un minimiseur sans contraintes (BFGS) à la place.

import time

start_time = time.time()

result = minimize(
    cost_func_vqe, x0, args=(raw_ansatz, observable_1, estimator), method="BFGS"
)

end_time = time.time()
execution_time = end_time - start_time

print("CHANGED TO BFGS OPTIMIZER:")
print(f"""Number of iterations: {result.nfev}""")
print(f"""Time (s): {execution_time}""")

CHANGED TO BFGS OPTIMIZER:
Number of iterations: 117
Time (s): 0.31656408309936523

Exemple de VQD

Nous implémentons ici le cadre des patterns Qiskit de façon explicite.

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

Plutôt que de chercher uniquement la valeur propre la plus basse de nos observables, nous allons en chercher les $4$ (où $k=4$).

Rappelons que les fonctions de coût de VQD sont :

$$C_{l}(\vec{\theta}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle + \sum_{j=0}^{l-1}\beta_j |\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2 \quad \forall l\in\{1,\cdots,k\}=\{1,\cdots,4\}$$

Ceci est particulièrement important car un vecteur $\vec{\beta}=(\beta_0,\cdots,\beta_{k-1})$ (dans ce cas $(\beta_0, \beta_1, \beta_2, \beta_3)$) doit être passé en argument lors de la définition de l'objet VQD.

De plus, dans l'implémentation Qiskit de VQD, au lieu de considérer les observables effectifs décrits dans le notebook précédent, les fidélités $|\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2$ sont calculées directement via l'algorithme ComputeUncompute, qui exploite une primitive Sampler pour échantillonner la probabilité d'obtenir $|0\rangle$ pour le circuit $U_A^\dagger(\vec{\theta})U_A(\vec{\theta^j})$. Cela fonctionne précisément parce que cette probabilité est

$$\begin{aligned} p_0 & = |\langle 0|U_A^\dagger(\vec{\theta})U_A(\vec{\theta^j})|0\rangle|^2 \\[1mm] & = |\big(\langle 0|U_A^\dagger(\vec{\theta})\big)\big(U_A(\vec{\theta^j})|0\rangle\big)|^2 \\[1mm] & = |\langle \psi(\vec{\theta}) |\psi(\vec{\theta^j}) \rangle|^2 \\[1mm] \end{aligned}$$

ansatz = n_local(
    num_qubits=2,
    rotation_blocks=["ry", "rz"],
    entanglement_blocks="cz",
    #    entanglement="linear",
    reps=1,
)

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

Commençons par examiner l'observable suivante :

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

Cet observable a les valeurs propres suivantes :

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

Et les états propres :

$$\left\{ \begin{array}{c} |\phi_0\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)\\ |\phi_1\rangle = \frac{1}{\sqrt{2}}(|00\rangle - |11\rangle)\\ |\phi_2\rangle = \frac{1}{\sqrt{2}}(|01\rangle + |10\rangle)\\ |\phi_3\rangle = \frac{1}{\sqrt{2}}(|01\rangle - |10\rangle) \end{array} \right\}$$

from qiskit.quantum_info import SparsePauliOp

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

Nous allons utiliser la fonction suivante pour calculer la pénalité de chevauchement. Notez qu'il s'agit toujours d'une partie du mappage du problème vers des circuits quantiques. Cependant, comme discuté dans la leçon précédente, cette fonction calcule le chevauchement entre un circuit variationnel courant et le circuit optimisé d'un état précédent à énergie/coût plus faible obtenu. Le nouveau circuit généré doit également être transpilé pour fonctionner sur du matériel réel. Nous avons déjà vu cette fonction, utilisée sur un simulateur. Ici, nous devons déjà envisager la transpilation et l'optimisation associée pour quand nous utilisons un backend réel, d'où les lignes autour de if realbackend == 1. Cela mélange un peu l'étape 2, mais nous le signalerons explicitement plus tard.

import numpy as np
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

def calculate_overlaps(
    ansatz, prev_circuits, parameters, sampler, realbackend, backend
):
    def create_fidelity_circuit(circuit_1, circuit_2):
        if len(circuit_1.clbits) > 0:
            circuit_1.remove_final_measurements()
        if len(circuit_2.clbits) > 0:
            circuit_2.remove_final_measurements()

        circuit = circuit_1.compose(circuit_2.inverse())
        circuit.measure_all()
        return circuit

    overlaps = []

    for prev_circuit in prev_circuits:
        fidelity_circuit = create_fidelity_circuit(ansatz, prev_circuit)
        if realbackend == 1:
            pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
            fidelity_circuit = pm.run(fidelity_circuit)
        sampler_job = sampler.run([(fidelity_circuit, parameters)])
        meas_data = sampler_job.result()[0].data.meas

        counts_0 = meas_data.get_int_counts().get(0, 0)
        shots = meas_data.num_shots
        overlap = counts_0 / shots
        overlaps.append(overlap)

    return np.array(overlaps)

Nous ajoutons maintenant la fonction de coût de VQD. Notez que par rapport à la leçon précédente, nous avons maintenant deux arguments supplémentaires (realbackend et backend) pour nous aider avec la transpilation lorsque nous utilisons des backends réels.

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

    total_cost = 0

    if step > 1:
        overlaps = calculate_overlaps(
            ansatz, prev_states, parameters, sampler, realbackend, backend
        )
        total_cost = np.sum(
            [np.real(betas[state] * overlap) for state, overlap in enumerate(overlaps)]
        )

    estimator_result = estimator_job.result()[0]

    value = estimator_result.data.evs[0] + total_cost

    return value

Encore une fois, nous utiliserons des simulateurs pour le débogage, puis passerons au matériel réel.

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

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

Ici nous introduisons le nombre d'états que nous souhaitons calculer, les pénalités, et un ensemble de paramètres initiaux, x0.

from qiskit.quantum_info import SparsePauliOp

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

Nous allons maintenant tester l'algorithme en utilisant des simulateurs :

from scipy.optimize import minimize

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

realbackend = 0

for step in range(1, k + 1):
    if step > 1:
        prev_states.append(ansatz.assign_parameters(prev_opt_parameters))

    result = minimize(
        cost_func_vqd,
        x0,
        args=(
            ansatz,
            prev_states,
            step,
            betas,
            estimator,
            sampler,
            observable_2,
            realbackend,
            None,
        ),
        method="COBYLA",
        options={"maxiter": 200, "tol": 0.000001},
    )
    print(result)

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

message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: -6.9999999999996
       x: [ 1.571e+00  1.571e+00  2.519e+00  2.100e+00  1.242e+00
            6.935e-01  2.298e+00  1.991e+00]
    nfev: 151
   maxcv: 0.0
 message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: 3.698974255258432
       x: [ 1.269e+00  1.109e+00  1.080e+00  1.200e+00  1.094e+00
            1.163e+00  9.752e-01  9.519e-01]
    nfev: 103
   maxcv: 0.0
 message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: 4.731320121938101
       x: [ 1.533e+00  2.451e+00  2.526e+00  2.406e+00  1.968e+00
            2.105e+00  8.537e-01  8.442e-01]
    nfev: 110
   maxcv: 0.0
 message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: 7.008239313655201
       x: [ 4.150e+00  2.120e+00  3.495e+00  7.262e-01  1.953e+00
           -1.982e-01  3.263e-01  2.563e+00]
    nfev: 126
   maxcv: 0.0

eigenvalues

[np.float64(-6.9999999999996),
 np.float64(3.698974255258432),
 np.float64(4.731320121938101),
 np.float64(7.008239313655201)]

Ces résultats sont assez proches des valeurs attendues, à l'exception de l'erreur d'approximation et de la phase globale. Nous pourrions ajuster la tolérance de l'optimiseur classique et/ou les pénalités pour le chevauchement des vecteurs d'état afin d'obtenir des valeurs plus précises.

solution_eigenvalues = [-7, 3, 5, 7]

for index, experimental_eigenvalue in enumerate(eigenvalues):
    solution_eigenvalue = solution_eigenvalues[index]

    print(
        f"Percent error: {abs((experimental_eigenvalue - solution_eigenvalue)/solution_eigenvalue):.2e}"
    )

Percent error: 5.71e-14
Percent error: 2.33e-01
Percent error: 5.37e-02
Percent error: 1.18e-03

Changer les betas

Comme mentionné dans la leçon précédente, les valeurs de $\vec{\beta}$ doivent être supérieures à la différence entre les valeurs propres. Voyons ce qui se passe quand cette condition n'est pas respectée avec $\hat{O}_2$ :

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

avec les valeurs propres

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

from qiskit.quantum_info import SparsePauliOp

k = 4
betas = np.ones(3)
x0 = np.zeros(8)

from scipy.optimize import minimize

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

realbackend = 0

for step in range(1, k + 1):
    if step > 1:
        prev_states.append(ansatz.assign_parameters(prev_opt_parameters))

    result = minimize(
        cost_func_vqd,
        x0,
        args=(
            ansatz,
            prev_states,
            step,
            betas,
            estimator,
            sampler,
            observable_2,
            realbackend,
            None,
        ),
        method="COBYLA",
        options={"tol": 0.01, "maxiter": 200},
    )
    print(result)

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

message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: -6.999916534745094
       x: [ 1.568e+00 -1.569e+00  1.385e-01  1.398e-01 -7.972e-01
            7.835e-01 -2.375e-01  4.539e-02]
    nfev: 125
   maxcv: 0.0
 message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: -1.515139929812874
       x: [-5.317e-04 -2.514e-03  1.016e+00  9.998e-01  3.890e-04
            1.772e-04  1.568e-04  8.497e-04]
    nfev: 35
   maxcv: 0.0
 message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: -0.509948114293115
       x: [-3.796e-03  8.853e-03  3.015e-04  9.997e-01  6.271e-04
           -2.554e-03  1.017e-04  2.766e-04]
    nfev: 37
   maxcv: 0.0
 message: Return from COBYLA because the trust region radius reaches its lower bound.
 success: True
  status: 0
     fun: 0.4914672235935682
       x: [-7.178e-03 -8.652e-03  1.125e+00 -5.428e-02 -1.586e-03
            2.031e-03 -3.462e-03  5.734e-03]
    nfev: 35
   maxcv: 0.0

solution_eigenvalues = [-7, 3, 5, 7]

for index, experimental_eigenvalue in enumerate(eigenvalues):
    solution_eigenvalue = solution_eigenvalues[index]

    print(
        f"Percent error: {abs((experimental_eigenvalue - solution_eigenvalue)/solution_eigenvalue):.2e}"
    )

Percent error: 1.19e-05
Percent error: 1.51e+00
Percent error: 1.10e+00
Percent error: 9.30e-01

Cette fois, l'optimiseur renvoie le même état $|\phi_0\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)$ comme solution proposée pour tous les états propres : ce qui est clairement incorrect. Cela se produit parce que les betas étaient trop petits pour pénaliser l'état propre minimal dans les fonctions de coût successives. Par conséquent, il n'a pas été exclu de l'espace de recherche effectif dans les itérations ultérieures de l'algorithme, et a toujours été choisi comme la meilleure solution possible.

Nous recommandons d'expérimenter avec les valeurs de $\vec{\beta}$ et de s'assurer qu'elles sont supérieures à la différence entre les valeurs propres.

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

Pour exécuter ceci sur du matériel réel, nous devons optimiser les circuits quantiques pour l'ordinateur quantique de notre choix. Pour nos besoins ici, nous utiliserons simplement le backend le moins occupé.

from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit_ibm_runtime import Session, EstimatorOptions
from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
# Or use a specific backend
# backend = service.backend("ibm_brisbane")
print(backend)

<IBMBackend('ibm_brisbane')>

Nous allons transpiler notre circuit en utilisant un gestionnaire de passes préconfiguré et le niveau d'optimisation 3.

pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_ansatz = pm.run(ansatz)
isa_observable = observable_2.apply_layout(layout=isa_ansatz.layout)

Étape 3 : Exécuter en utilisant les primitives Qiskit

En veillant à remettre nos betas à des valeurs suffisamment élevées, nous pouvons maintenant exécuter notre calcul sur du vrai matériel quantique.

# Estimated compute resource usage: 25 minutes.
# Benchmarked at 24 min, 30 sec on an Eagle r3 processor on 5-30-24

k = 2
betas = [30, 50, 80]
x0 = np.zeros(8)

real_prev_states = []
real_prev_opt_parameters = []
real_eigenvalues = []

realbackend = 1

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

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

    for step in range(1, k + 1):
        if step > 1:
            real_prev_states.append(isa_ansatz.assign_parameters(prev_opt_parameters))

        result = minimize(
            cost_func_vqd,
            x0,
            args=(
                isa_ansatz,
                real_prev_states,
                step,
                betas,
                estimator,
                sampler,
                isa_observable,
                realbackend,
                backend,
            ),
            method="COBYLA",
            options={"maxiter": 200},
        )
        print(result)

        real_prev_opt_parameters = result.x
        real_eigenvalues.append(result.fun)

session.close()
print(real_eigenvalues)

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

Notre sortie est structurellement similaire à ce qui a été discuté dans les leçons et exemples précédents. Mais il y a quelque chose de problématique dans les résultats ci-dessus, duquel nous pouvons tirer un message d'avertissement pour le contexte des états excités. Pour limiter le temps de calcul utilisé sur cet exemple d'apprentissage, nous avons fixé un nombre maximum d'itérations pour l'optimiseur classique qui était potentiellement trop bas : 200 itérations. Un calcul précédent ci-dessus, sur un simulateur, n'a pas réussi à converger en 200 itérations. Ici, le nôtre a convergé... mais vers quelle tolérance ? Nous n'avons pas spécifié de tolérance pour que COBYLA se considère "convergé". Un coup d'œil à la valeur de la fonction et la comparaison avec les exécutions précédentes nous dit que COBYLA n'était pas proche de converger vers la précision que nous exigeons.

Il y a un autre problème : l'énergie du premier état excité semble être inférieure à l'énergie de l'état fondamental ! Essayez d'expliquer comment cela pourrait se produire. Indice : c'est lié au point de convergence que nous venons d'aborder. Ce comportement est expliqué en détail ci-dessous après l'application de VQD à la molécule H2.

Chimie quantique : solveur d'état fondamental et d'états excités

Notre objectif est de minimiser la valeur attendue de l'observable représentant l'énergie (Hamiltonien $\hat{\mathcal{H}}$) :

$$\min_{\vec\theta} \langle\psi(\vec\theta)|\hat{\mathcal{H}}|\psi(\vec\theta)\rangle$$

from qiskit.quantum_info import SparsePauliOp
from qiskit.circuit.library import efficient_su2

H2_op = SparsePauliOp.from_list(
    [
        ("II", -1.052373245772859),
        ("IZ", 0.39793742484318045),
        ("ZI", -0.39793742484318045),
        ("ZZ", -0.01128010425623538),
        ("XX", 0.18093119978423156),
    ]
)

chem_ansatz = efficient_su2(H2_op.num_qubits)

chem_ansatz.decompose().draw("mpl")

from qiskit import QuantumCircuit

def cost_func_vqe(params, ansatz, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """
    pub = (ansatz, hamiltonian, params)
    cost = estimator.run([pub]).result()[0].data.evs
    #    cost = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
    return cost

Nous définissons maintenant un ensemble de paramètres initiaux :

import numpy as np

x0 = np.ones(chem_ansatz.num_parameters)

Nous pouvons minimiser cette fonction de coût pour calculer les paramètres optimaux, et nous pouvons d'abord vérifier notre code en utilisant un simulateur local.

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

estimator = Estimator()
sampler = Sampler()

# SciPy minimizer routine
from scipy.optimize import minimize
import time

start_time = time.time()

result = minimize(
    cost_func_vqe, x0, args=(chem_ansatz, H2_op, estimator), method="COBYLA"
)

end_time = time.time()
execution_time = end_time - start_time

result

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -1.857275029048451
       x: [ 7.326e-01  1.354e+00 ...  1.040e+00  1.508e+00]
    nfev: 242
   maxcv: 0.0

La valeur minimale de la fonction de coût (-1.857...) est l'énergie de l'état fondamental de la molécule H2, en unités de hartrees.

États excités

Nous pouvons également exploiter VQD pour résoudre $k=2$ états au total (l'état fondamental et le premier état excité).

from qiskit.quantum_info import SparsePauliOp
import numpy as np

k = 2
betas = [33, 33]
# x0 = np.zeros(ansatz.num_parameters)
x0 = [
    1.164e00,
    -2.438e-01,
    9.358e-04,
    6.745e-02,
    1.990e00,
    9.810e-02,
    6.154e-01,
    5.454e-01,
]

Nous allons ajouter notre calcul de chevauchement :

from scipy.optimize import minimize

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

realbackend = 0

for step in range(1, k + 1):
    if step > 1:
        prev_states.append(ansatz.assign_parameters(prev_opt_parameters))

    result = minimize(
        cost_func_vqd,
        x0,
        args=(
            ansatz,
            prev_states,
            step,
            betas,
            estimator,
            sampler,
            H2_op,
            realbackend,
            None,
        ),
        method="COBYLA",
        options={"tol": 0.001, "maxiter": 2000},
    )
    print(result)

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

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -1.8572671093941977
       x: [ 1.164e+00 -2.437e-01  2.118e-03  6.448e-02  1.990e+00
            9.870e-02  6.167e-01  5.476e-01]
    nfev: 58
   maxcv: 0.0
 message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -1.0322873777662176
       x: [ 3.205e+00  1.502e+00  1.699e+00 -1.107e-02  3.086e+00
            1.530e+00  4.445e-02  7.013e-02]
    nfev: 99
   maxcv: 0.0

eigenvalues

[-1.8572671093941977, -1.0322873777662176]

Matériel réel et un dernier message d'avertissement

Pour exécuter ceci sur du matériel réel, nous devons optimiser les circuits quantiques pour l'ordinateur quantique de notre choix. Pour nos besoins ici, nous utiliserons simplement le backend le moins occupé.

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)

Nous allons utiliser un gestionnaire de passes préconfiguré pour la transpilation, et nous allons optimiser au maximum notre circuit en utilisant le niveau d'optimisation 3.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_ansatz = pm.run(ansatz)
isa_observable = H2_op.apply_layout(layout=isa_ansatz.layout)

Étant donné que VQD est très itératif, nous allons effectuer toutes les étapes dans une session Runtime, de sorte que nos tâches ne seront en file d'attente qu'au début, et non entre chaque mise à jour des paramètres. Rien d'autre ne change dans la syntaxe de la fonction de coût ou de l'estimateur.

x0 = [
    1.306e00,
    -2.284e-01,
    6.913e-02,
    -2.530e-02,
    1.849e00,
    7.433e-02,
    6.366e-01,
    5.600e-01,
]

# Estimated hardware usage: 20 min benchmarked on an Eagle r3 processor on 5-30-24

real_prev_states = []
real_prev_opt_parameters = []
real_eigenvalues = []

realbackend = 1

estimator_options = EstimatorOptions(resilience_level=1, default_shots=4096)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session)
    sampler = Sampler(mode=session)

    for step in range(1, k + 1):
        if step > 1:
            real_prev_states.append(
                isa_ansatz.assign_parameters(real_prev_opt_parameters)
            )

        result = minimize(
            cost_func_vqd,
            x0,
            args=(
                isa_ansatz,
                real_prev_states,
                step,
                betas,
                estimator,
                sampler,
                isa_observable,
                realbackend,
                backend,
            ),
            method="COBYLA",
            options={"tol": 0.001, "maxiter": 300},
        )
        print(result)

        real_prev_opt_parameters = result.x
        real_eigenvalues.append(result.fun)

session.close()
print(real_eigenvalues)

L'énergie de l'état fondamental obtenue (-1,83 hartrees) ne s'éloigne pas trop de la valeur correcte (-1,85 hartrees). Cependant, l'énergie de l'état excité est assez éloignée. Ceci est similaire au comportement erroné que nous avons vu plus tôt dans cette leçon. L'énergie rapportée pour l'état excité est presque la même que celle de l'état fondamental. Dans le cas précédent, nous avons même vu une énergie de l'état excité qui était inférieure à l'énergie de l'état fondamental rapportée.

Il n'est pas possible qu'un calcul variationnel produise une énergie inférieure à la vraie énergie de l'état fondamental. Dans le cas précédent, l'énergie de l'état fondamental que nous avons obtenue n'était pas très proche du vrai état fondamental. Comme nous n'avons pas obtenu la vraie énergie de l'état fondamental dans ce cas, il n'y a pas de contradiction. Dans le cas présent, l'énergie de l'état fondamental était assez proche de la valeur correcte, et pourtant l'énergie de l'état excité semble étrangement proche de cette même valeur.

Pour mieux comprendre comment cela s'est produit, rappelons que la façon dont nous trouvons un état excité est d'exiger que l'état variationnel soit orthogonal à l'état fondamental (en utilisant les circuits de chevauchement et les termes de pénalité). Si nous ne parvenons pas à obtenir une énergie précise de l'état fondamental (ou si nous nous trompons de quelques pour cent), alors nous ne parvenons pas non plus à obtenir un vecteur précis de l'état fondamental ! Ainsi, lorsque nous exigeons que l'état excité soit orthogonal au premier état que nous avons trouvé, nous n'imposions pas l'orthogonalité avec le vrai état fondamental, mais plutôt avec une approximation de celui-ci (parfois une approximation médiocre). Ainsi, l'état excité n'a pas été contraint d'être orthogonal au vrai état fondamental, et nos estimations d'énergie pour les états excités étaient en fait assez proches de l'énergie de l'état fondamental.

Ce sera toujours une préoccupation dans VQD. Mais en principe, cela peut être corrigé en augmentant le nombre maximum d'itérations pour l'optimiseur classique, en imposant une tolérance plus faible pour l'optimiseur classique, et éventuellement en essayant un ansatz différent si nous manquons habituellement le vrai état fondamental. Comme nous l'avons vu, il peut également être nécessaire de modifier les pénalités de chevauchement (betas). Mais c'est vraiment un problème séparé. Aucune pénalité pour le chevauchement ne t'éloignera du vrai état fondamental si tu n'as pas trouvé une très bonne estimation du vrai état fondamental pour le circuit de chevauchement.

Optimisation : Max-Cut

Le problème de coupe maximale (max-cut) est un problème d'optimisation combinatoire qui consiste à diviser les sommets d'un graphe en deux ensembles disjoints de façon à maximiser le nombre d'arêtes entre les deux ensembles. Plus formellement, étant donné un graphe non orienté $G=(V,E)$, où $V$ est l'ensemble des sommets et $E$ est l'ensemble des arêtes, le problème max-cut demande de partitionner les sommets en deux sous-ensembles disjoints, $S$ et $T$, de façon à maximiser le nombre d'arêtes ayant une extrémité dans $S$ et l'autre dans $T$.

Nous pouvons appliquer le max-cut pour résoudre divers problèmes, tels que le regroupement, la conception de réseaux et les transitions de phase. Nous commencerons par créer un graphe du problème :

import rustworkx as rx
from rustworkx.visualization import mpl_draw

n = 4
G = rx.PyGraph()
G.add_nodes_from(range(n))
# The edge syntax is (start, end, weight)
edges = [(0, 1, 1.0), (0, 2, 1.0), (0, 3, 1.0), (1, 2, 1.0), (2, 3, 1.0)]
G.add_edges_from(edges)

mpl_draw(
    G, pos=rx.shell_layout(G), with_labels=True, edge_labels=str, node_color="#1192E8"
)

Ce problème peut être exprimé comme un problème d'optimisation binaire. Pour chaque nœud $0 \leq i < n$, où $n$ est le nombre de nœuds du graphe (dans ce cas $n=4$), nous allons considérer la variable binaire $x_i$. Cette variable aura la valeur $1$ si le nœud $i$ appartient à l'un des groupes que nous allons appeler $1$ et $0$ s'il est dans l'autre groupe, que nous allons appeler $0$. Nous noterons également $w_{ij}$ (élément $(i,j)$ de la matrice d'adjacence $w$) le poids de l'arête qui va du nœud $i$ au nœud $j$. Parce que le graphe est non orienté, $w_{ij}=w_{ji}$. Nous pouvons alors formuler notre problème comme la maximisation de la fonction de coût suivante :

$$\begin{aligned} C(\vec{x}) & =\sum_{i,j=0}^n w_{ij} x_i(1-x_j)\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i,j=0}^n w_{ij} x_ix_j\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i=0}^n \sum_{j=0}^i 2w_{ij} x_ix_j \end{aligned}$$

Pour résoudre ce problème avec un ordinateur quantique, nous allons exprimer la fonction de coût comme la valeur attendue d'un observable. Cependant, les observables que Qiskit admet nativement consistent en opérateurs de Pauli, qui ont des valeurs propres $1$ et $-1$ au lieu de $0$ et $1$. C'est pourquoi nous allons effectuer le changement de variable suivant :

Où $\vec{x}=(x_0,x_1,\cdots ,x_{n-1})$. Nous pouvons utiliser la matrice d'adjacence $w$ pour accéder facilement aux poids de toutes les arêtes. Cela sera utilisé pour obtenir notre fonction de coût :

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

Ce qui implique que :

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

Donc la nouvelle fonction de coût que nous voulons maximiser est :

$$\begin{aligned} C(\vec{z}) & = \sum_{i,j=0}^n w_{ij} \bigg(\frac{1-z_i}{2}\bigg)\bigg(1-\frac{1-z_j}{2}\bigg)\\[1mm] & = \sum_{i,j=0}^n \frac{w_{ij}}{4} - \sum_{i,j=0}^n \frac{w_{ij}}{4} z_iz_j\\[1mm] & = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j \end{aligned}$$

De plus, la tendance naturelle d'un ordinateur quantique est de trouver des minima (généralement l'énergie la plus basse) au lieu de maxima, donc au lieu de maximiser $C(\vec{z})$, nous allons minimiser :

$$-C(\vec{z}) = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}$$

Maintenant que nous avons une fonction de coût à minimiser dont les variables peuvent prendre les valeurs $-1$ et $1$, nous pouvons faire l'analogie suivante avec le Pauli $Z$ :

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

En d'autres termes, la variable $z_i$ sera équivalente à une porte $Z$ agissant sur le qubit $i$. De plus :

$$Z_i|x_{n-1}\cdots x_0\rangle = z_i|x_{n-1}\cdots x_0\rangle \rightarrow \langle x_{n-1}\cdots x_0 |Z_i|x_{n-1}\cdots x_0\rangle = z_i$$

Alors l'observable que nous allons considérer est :

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

auquel nous devrons ajouter le terme indépendant par la suite :

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

L'opérateur est une combinaison linéaire de termes avec des opérateurs Z sur les nœuds connectés par une arête (rappelons que le qubit 0 est le plus à droite) : $IIZZ + IZIZ + IZZI + ZIIZ + ZZII$. Une fois l'opérateur construit, l'ansatz pour l'algorithme QAOA peut facilement être construit en utilisant le circuit QAOAAnsatz de la bibliothèque de circuits Qiskit.

from qiskit.circuit.library import QAOAAnsatz
from qiskit.quantum_info import SparsePauliOp

max_hamiltonian = SparsePauliOp.from_list(
    [("IIZZ", 1), ("IZIZ", 1), ("IZZI", 1), ("ZIIZ", 1), ("ZZII", 1)]
)

max_ansatz = QAOAAnsatz(max_hamiltonian, reps=2)
# Draw
max_ansatz.decompose(reps=3).draw("mpl")

# Sum the weights, and divide by 2

offset = -sum(edge[2] for edge in edges) / 2
print(f"""Offset: {offset}""")

Offset: -2.5

def cost_func(params, ansatz, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """
    pub = (ansatz, hamiltonian, params)
    cost = estimator.run([pub]).result()[0].data.evs
    #    cost = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
    return cost

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

estimator = Estimator()
sampler = Sampler()

Nous définissons maintenant un ensemble de paramètres aléatoires initiaux :

import numpy as np

x0 = 2 * np.pi * np.random.rand(max_ansatz.num_parameters)
print(x0)

[6.0252949  0.58448176 2.15785731 1.13646074]

N'importe quel optimiseur classique peut être utilisé pour minimiser la fonction de coût. Sur un système quantique réel, un optimiseur conçu pour des paysages de fonctions de coût non lisses fait généralement mieux. Ici, nous utilisons la routine COBYLA de SciPy via la fonction minimize.

Étant donné que nous exécutons de nombreux appels à Runtime de façon itérative, nous utilisons une session pour exécuter tous les appels dans un seul bloc. De plus, pour QAOA, la solution est encodée dans la distribution de sortie du circuit ansatz lié aux paramètres optimaux de la minimisation. Par conséquent, nous aurons besoin d'une primitive Sampler, et nous l'instancierons avec la même session Et nous lançons notre routine de minimisation :

result = minimize(
    cost_func, x0, args=(max_ansatz, max_hamiltonian, estimator), method="COBYLA"
)

print(result)

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -2.585287311689236
       x: [ 7.332e+00  3.904e-01  2.045e+00  1.028e+00]
    nfev: 80
   maxcv: 0.0

Le vecteur solution des angles de paramètres (x), lorsqu'il est introduit dans le circuit ansatz, produit la partition du graphe que nous recherchions.

eigenvalue = cost_func(result.x, max_ansatz, max_hamiltonian, estimator)
print(f"""Eigenvalue: {eigenvalue}""")
print(f"""Max-Cut Objective: {eigenvalue + offset}""")

Eigenvalue: -2.585287311689236
Max-Cut Objective: -5.085287311689235

from qiskit.result import QuasiDistribution
from qiskit.primitives import StatevectorSampler

sampler = StatevectorSampler()

# Assign solution parameters to ansatz
qc = max_ansatz.assign_parameters(result.x)

# Add measurements to our circuit
qc.measure_all()

# Sample ansatz at optimal parameters
# samp_dist = sampler.run(qc).result().quasi_dists[0]

shots = 1024
job = sampler.run([qc], shots=shots)

qc.decompose().draw("mpl")

data_pub = job.result()[0].data
bitstrings = data_pub.meas.get_bitstrings()
counts = data_pub.meas.get_counts()
quasi_dist = QuasiDistribution(
    {outcome: freq / shots for outcome, freq in counts.items()}
)
probabilities = quasi_dist

# Close the session since we are now done with it
# session.close()

from qiskit.visualization import plot_distribution

plot_distribution(counts)

binary_string = max(counts.items(), key=lambda kv: kv[1])[0]
x = np.asarray([int(y) for y in reversed(list(binary_string))])

colors = ["r" if x[i] == 0 else "c" for i in range(n)]
mpl_draw(
    G, pos=rx.shell_layout(G), with_labels=True, edge_labels=str, node_color=colors
)

Résumé

Avec cette leçon, tu as appris :

• Comment écrire un algorithme variationnel personnalisé
• Comment appliquer un algorithme variationnel pour trouver les valeurs propres minimales
• Comment utiliser les algorithmes variationnels pour résoudre des cas d'usage applicatifs

Passe à la leçon finale pour passer ton évaluation et obtenir ton badge !