Diagonalisation quantique de Krylov des hamiltoniens de réseau

Estimation de l'utilisation : 20 minutes sur un Héron r2 (NOTE : Il s'agit uniquement d'une estimation. Votre durée d'exécution peut varier.)

Arrière-plan

Ce tutoriel montre comment mettre en œuvre l'algorithme de diagonalisation quantique de Krylov (KQD) dans le contexte des modèles Qiskit. Vous découvrirez d'abord la théorie qui sous-tend l'algorithme, puis vous assisterez à une démonstration de son exécution sur une QPU.

Toutes disciplines confondues, nous nous intéressons à l'apprentissage des propriétés de l'état fondamental des systèmes quantiques. Il s'agit par exemple de comprendre la nature fondamentale des particules et des forces, de prévoir et de comprendre le comportement de matériaux complexes et de comprendre les interactions et les réactions biochimiques. En raison de la croissance exponentielle de l'espace de Hilbert et des corrélations qui apparaissent dans les systèmes intriqués, les algorithmes classiques peinent à résoudre ce problème pour des systèmes quantiques de taille croissante. À une extrémité du spectre se trouve l'approche existante qui tire parti du matériel quantique en se concentrant sur les méthodes quantiques variationnelles (par exemple, l' eigensolver quantique variationnel ). Ces techniques se heurtent à des difficultés avec les dispositifs actuels en raison du nombre élevé d'appels de fonction requis dans le processus d'optimisation, ce qui ajoute un surcoût important en termes de ressources une fois que des techniques avancées d'atténuation des erreurs sont introduites, limitant ainsi leur efficacité aux systèmes de petite taille. À l'autre extrémité du spectre, il existe des méthodes quantiques tolérantes aux fautes avec des garanties de performance (par exemple, l' estimation quantique de la phase ), qui nécessitent des circuits profonds qui ne peuvent être exécutés que sur un dispositif tolérant aux fautes. Pour ces raisons, nous présentons ici un algorithme quantique basé sur des méthodes de sous-espace (telles que décrites dans cet article ), l'algorithme de diagonalisation quantique de Krylov (KQD). Cet algorithme fonctionne bien à grande échelle [1] sur le matériel quantique existant, offre des garanties de performance similaires à celles de l'estimation de phase, est compatible avec des techniques avancées d'atténuation des erreurs et pourrait fournir des résultats qui sont classiquement inaccessibles.

Exigences

Avant de commencer ce tutoriel, assurez-vous que les éléments suivants sont installés :

Qiskit SDK v2.0 ou plus tard, avec prise en charge de la visualisation
Qiskit Runtime v0.22 ou plus tard ( pip install qiskit-ibm-runtime )

Configuration

import numpy as np
import scipy as sp
import matplotlib.pylab as plt
from typing import Union, List
import itertools as it
import copy
from sympy import Matrix
import warnings

warnings.filterwarnings("ignore")

from qiskit.quantum_info import SparsePauliOp, Pauli, StabilizerState
from qiskit.circuit import Parameter, IfElseOp
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter
from qiskit.transpiler import Target, CouplingMap
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager


from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    EstimatorV2 as Estimator,
)


def solve_regularized_gen_eig(
    h: np.ndarray,
    s: np.ndarray,
    threshold: float,
    k: int = 1,
    return_dimn: bool = False,
) -> Union[float, List[float]]:
    """
    Method for solving the generalized eigenvalue problem with regularization

    Args:
        h (numpy.ndarray):
            The effective representation of the matrix in the Krylov subspace
        s (numpy.ndarray):
            The matrix of overlaps between vectors of the Krylov subspace
        threshold (float):
            Cut-off value for the eigenvalue of s
        k (int):
            Number of eigenvalues to return
        return_dimn (bool):
            Whether to return the size of the regularized subspace

    Returns:
        lowest k-eigenvalue(s) that are the solution of the
        regularized generalized eigenvalue problem


    """
    s_vals, s_vecs = sp.linalg.eigh(s)
    s_vecs = s_vecs.T
    good_vecs = np.array(
        [vec for val, vec in zip(s_vals, s_vecs) if val > threshold]
    )
    h_reg = good_vecs.conj() @ h @ good_vecs.T
    s_reg = good_vecs.conj() @ s @ good_vecs.T
    if k == 1:
        if return_dimn:
            return sp.linalg.eigh(h_reg, s_reg)[0][0], len(good_vecs)
        else:
            return sp.linalg.eigh(h_reg, s_reg)[0][0]
    else:
        if return_dimn:
            return sp.linalg.eigh(h_reg, s_reg)[0][:k], len(good_vecs)
        else:
            return sp.linalg.eigh(h_reg, s_reg)[0][:k]


def single_particle_gs(H_op, n_qubits):
    """
    Find the ground state of the single particle(excitation) sector
    """
    H_x = []
    for p, coeff in H_op.to_list():
        H_x.append(set([i for i, v in enumerate(Pauli(p).x) if v]))

    H_z = []
    for p, coeff in H_op.to_list():
        H_z.append(set([i for i, v in enumerate(Pauli(p).z) if v]))

    H_c = H_op.coeffs

    print("n_sys_qubits", n_qubits)

    n_exc = 1
    sub_dimn = int(sp.special.comb(n_qubits + 1, n_exc))
    print("n_exc", n_exc, ", subspace dimension", sub_dimn)

    few_particle_H = np.zeros((sub_dimn, sub_dimn), dtype=complex)

    # list all of the possible sets of n_exc indices of 1s in
    # n_exc-particle states
    sparse_vecs = [
        set(vec) for vec in it.combinations(range(n_qubits + 1), r=n_exc)
    ]

    m = 0
    for i, i_set in enumerate(sparse_vecs):
        for j, j_set in enumerate(sparse_vecs):
            m += 1

            if len(i_set.symmetric_difference(j_set)) <= 2:
                for p_x, p_z, coeff in zip(H_x, H_z, H_c):
                    if i_set.symmetric_difference(j_set) == p_x:
                        sgn = ((-1j) ** len(p_x.intersection(p_z))) * (
                            (-1) ** len(i_set.intersection(p_z))
                        )
                    else:
                        sgn = 0

                    few_particle_H[i, j] += sgn * coeff

    gs_en = min(np.linalg.eigvalsh(few_particle_H))
    print("single particle ground state energy: ", gs_en)
    return gs_en

Étape 1 : Mettre en correspondance les entrées classiques avec un problème quantique

L'espace de Krylov

L'espace de Krylov $\mathcal{K}^r$ d'ordre $r$ est l'espace couvert par les vecteurs obtenus en multipliant les puissances supérieures d'une matrice $A$ , jusqu'à $r-1$ , avec un vecteur de référence $\vert v \rangle$ .

\mathcal{K}^r = \left\{ \vert v \rangle, A \vert v \rangle, A^2 \vert v \rangle, ..., A^{r-1} \vert v \rangle \right\}

Si la matrice $A$ est l'hamiltonien $H$ , nous appellerons l'espace correspondant l'espace de Krylov des puissances $\mathcal{K}_P$ . Dans le cas où $A$ est l'opérateur d'évolution temporelle généré par l'hamiltonien $U=e^{-iHt}$ , nous appellerons l'espace l'espace de Krylov unitaire $\mathcal{K}_U$ . Le sous-espace de Krylov puissance que nous utilisons classiquement ne peut pas être généré directement sur un ordinateur quantique car $H$ n'est pas un opérateur unitaire. Au lieu de cela, nous pouvons utiliser l'opérateur d'évolution temporelle $U = e^{-iHt}$ , dont on peut montrer qu'il offre des garanties de convergence similaires à celles de la méthode de la puissance. Les puissances de $U$ deviennent alors des pas de temps différents $U^k = e^{-iH(kt)}$ .

\mathcal{K}_U^r = \left\{ \vert \psi \rangle, U \vert \psi \rangle, U^2 \vert \psi \rangle, ..., U^{r-1} \vert \psi \rangle \right\}

Voir l'annexe pour une dérivation détaillée de la manière dont l'espace de Krylov unitaire permet de représenter avec précision les états propres à faible énergie.

Algorithme de diagonalisation quantique de Krylov

Étant donné un hamiltonien $H$ que nous souhaitons diagonaliser, nous considérons d'abord l'espace de Krylov unitaire correspondant $\mathcal{K}_U$ . L'objectif est de trouver une représentation compacte de l'hamiltonien dans $\mathcal{K}_U$ , que nous appellerons $\tilde{H}$ . Les éléments de la matrice de $\tilde{H}$ , la projection de l'hamiltonien dans l'espace de Krylov, peuvent être calculés en calculant les valeurs d'espérance suivantes

\tilde{H}_{mn} = \langle \psi_m \vert H \vert \psi_n \rangle =

= \langle \psi \vert e^{i H t_m} H e^{-i H t_n} \vert \psi \rangle

= \langle \psi \vert e^{i H m dt} H e^{-i H n dt} \vert \psi \rangle

Où $\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle$ sont les vecteurs de l'espace de Krylov unitaire et $t_n = n dt$ sont les multiples du pas de temps $dt$ choisi. Sur un ordinateur quantique, le calcul de chaque élément de la matrice peut être effectué à l'aide de n'importe quel algorithme permettant d'obtenir un chevauchement entre les états quantiques. Ce tutoriel se concentre sur le test de Hadamard. Étant donné que le $\mathcal{K}_U$ a la dimension $r$ , l'hamiltonien projeté dans le sous-espace aura les dimensions $r \times r$ . Si $r$ est suffisamment petit (en général, $r<<100$ est suffisant pour obtenir une convergence des estimations des énergies propres), nous pouvons facilement diagonaliser l'hamiltonien projeté $\tilde{H}$ . Cependant, nous ne pouvons pas diagonaliser directement $\tilde{H}$ en raison de la non-orthogonalité des vecteurs de l'espace de Krylov. Nous devrons mesurer leurs chevauchements et construire une matrice $\tilde{S}$

\tilde{S}_{mn} = \langle \psi_m \vert \psi_n \rangle

Cela nous permet de résoudre le problème des valeurs propres dans un espace non orthogonal (également appelé problème généralisé des valeurs propres)

\tilde{H} \ \vec{c} = E \ \tilde{S} \ \vec{c}

On peut alors obtenir des estimations des valeurs propres et des états propres de $H$ en examinant ceux de $\tilde{H}$ . Par exemple, l'estimation de l'énergie de l'état fondamental est obtenue en prenant la plus petite valeur propre $c$ et l'état fondamental à partir du vecteur propre correspondant $\vec{c}$ . Les coefficients de $\vec{c}$ déterminent la contribution des différents vecteurs qui couvrent $\mathcal{K}_U$ .

La figure montre une représentation du circuit du test de Hadamard modifié, une méthode utilisée pour calculer le chevauchement entre différents états quantiques. Pour chaque élément de la matrice $\tilde{H}_{i,j}$ , un test de Hadamard entre les états $\vert \psi_i \rangle$ , $\vert \psi_j \rangle$ est effectué. Ceci est mis en évidence dans la figure par le schéma de couleurs pour les éléments de la matrice et les opérations $\text{Prep} \; \psi_i$ , $\text{Prep} \; \psi_j$ correspondantes. Ainsi, un ensemble de tests de Hadamard pour toutes les combinaisons possibles de vecteurs de l'espace de Krylov est nécessaire pour calculer tous les éléments de la matrice de l'hamiltonien projeté $\tilde{H}$ . Le fil supérieur du circuit de test de Hadamard est un qubit d'ancilla qui est mesuré soit dans la base X, soit dans la base Y. Sa valeur d'espérance détermine la valeur du chevauchement entre les états. Le fil du bas représente tous les qubits de l'hamiltonien du système. L'opération $\text{Prep} \; \psi_i$ prépare le qubit du système dans l'état $\vert \psi_i \rangle$ contrôlé par l'état du qubit ancillaire (de même pour $\text{Prep} \; \psi_j$ ) et l'opération $P$ représente la décomposition de Pauli de l'hamiltonien du système $H = \sum_i P_i$ . Une dérivation plus détaillée des opérations calculées par le test de Hadamard est donnée ci-dessous.

Définir Hamiltonien

Considérons l'hamiltonien de Heisenberg pour $N$ qubits sur une chaîne linéaire : $H= \sum_{i,j}^N X_i X_j + Y_i Y_j - J Z_i Z_j$

# Define problem Hamiltonian.
n_qubits = 30
J = 1  # coupling strength for ZZ interaction

# Define the Hamiltonian:
H_int = [["I"] * n_qubits for _ in range(3 * (n_qubits - 1))]
for i in range(n_qubits - 1):
    H_int[i][i] = "Z"
    H_int[i][i + 1] = "Z"
for i in range(n_qubits - 1):
    H_int[n_qubits - 1 + i][i] = "X"
    H_int[n_qubits - 1 + i][i + 1] = "X"
for i in range(n_qubits - 1):
    H_int[2 * (n_qubits - 1) + i][i] = "Y"
    H_int[2 * (n_qubits - 1) + i][i + 1] = "Y"
H_int = ["".join(term) for term in H_int]
H_tot = [(term, J) if term.count("Z") == 2 else (term, 1) for term in H_int]

# Get operator
H_op = SparsePauliOp.from_list(H_tot)
print(H_tot)

Output:

[('ZZIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IZZIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIZZIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIZZIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIZZIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIZZIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIZZIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIZZIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIZZIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIZZIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIZZIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIZZIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIZZIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIZZIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIZZIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIZZIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIZZIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIZZIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIZZIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIZZIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIZZIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIZZIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIZZIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIZZIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIZZIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIZZIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIZZII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIZZI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIZZ', 1), ('XXIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IXXIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIXXIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIXXIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIXXIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIXXIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIXXIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIXXIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIXXIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIXXIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIXXIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIXXIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIXXIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIXXIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIXXIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIXXIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIXXIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIXXIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIXXIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIXXIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIXXIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIXXIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIXXIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIXXIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIXXIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIXXIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIXXII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIXXI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIXX', 1), ('YYIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IYYIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIYYIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIYYIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIYYIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIYYIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIYYIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIYYIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIYYIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIYYIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIYYIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIYYIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIYYIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIYYIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIYYIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIYYIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIYYIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIYYIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIYYIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIYYIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIYYIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIYYIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIYYIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIYYIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIYYIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIYYIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIYYII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIYYI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIYY', 1)]

Définir les paramètres de l'algorithme

Nous choisissons de manière heuristique une valeur pour le pas de temps dt (sur la base des limites supérieures de la norme hamiltonienne). Ref [2] a montré qu'un pas de temps suffisamment petit est $\pi/\vert \vert H \vert \vert$ , et qu'il est préférable jusqu'à un certain point de sous-estimer cette valeur plutôt que de la surestimer, car une surestimation peut permettre aux contributions des états à haute énergie de corrompre même l'état optimal dans l'espace de Krylov. D'autre part, si l'on choisit $dt$ comme étant trop petit, le conditionnement du sous-espace de Krylov est moins bon, car les vecteurs de base de Krylov diffèrent moins d'un pas de temps à l'autre.

# Get Hamiltonian restricted to single-particle states
single_particle_H = np.zeros((n_qubits, n_qubits))
for i in range(n_qubits):
    for j in range(i + 1):
        for p, coeff in H_op.to_list():
            p_x = Pauli(p).x
            p_z = Pauli(p).z
            if all(
                p_x[k] == ((i == k) + (j == k)) % 2 for k in range(n_qubits)
            ):
                sgn = (
                    (-1j) ** sum(p_z[k] and p_x[k] for k in range(n_qubits))
                ) * ((-1) ** p_z[i])
            else:
                sgn = 0
            single_particle_H[i, j] += sgn * coeff
for i in range(n_qubits):
    for j in range(i + 1, n_qubits):
        single_particle_H[i, j] = np.conj(single_particle_H[j, i])

# Set dt according to spectral norm
dt = np.pi / np.linalg.norm(single_particle_H, ord=2)
dt

Output:

np.float64(0.10833078115826875)

Et définir d'autres paramètres de l'algorithme. Pour les besoins de ce tutoriel, nous nous contenterons d'utiliser un espace de Krylov à cinq dimensions seulement, ce qui est assez restrictif.

# Set parameters for quantum Krylov algorithm
krylov_dim = 5  # size of Krylov subspace
num_trotter_steps = 6
dt_circ = dt / num_trotter_steps

Préparation de l'État

Choisissez un état de référence $\vert \psi \rangle$ qui présente un certain chevauchement avec l'état fondamental. Pour cet hamiltonien, nous utilisons l'état a avec une excitation dans le qubit du milieu $\vert 00..010...00 \rangle$ comme état de référence.

qc_state_prep = QuantumCircuit(n_qubits)
qc_state_prep.x(int(n_qubits / 2) + 1)
qc_state_prep.draw("mpl", scale=0.5)

Output:

Évolution dans le temps

Nous pouvons réaliser l'opérateur d'évolution temporelle généré par un hamiltonien donné : $U=e^{-iHt}$ via l' approximation de Lie-Trotter.

t = Parameter("t")

## Create the time-evo op circuit
evol_gate = PauliEvolutionGate(
    H_op, time=t, synthesis=LieTrotter(reps=num_trotter_steps)
)

qr = QuantumRegister(n_qubits)
qc_evol = QuantumCircuit(qr)
qc_evol.append(evol_gate, qargs=qr)

Output:

<qiskit.circuit.instructionset.InstructionSet at 0x11eef9be0>

test de Hadamard

\begin{equation*} |0\rangle|0\rangle^N \quad\longrightarrow\quad \frac{1}{\sqrt{2}}\Big(|0\rangle + |1\rangle \Big)|0\rangle^N \quad\longrightarrow\quad \frac{1}{\sqrt{2}}\Big(|0\rangle|0\rangle^N+|1\rangle |\psi_i\rangle\Big) \quad\longrightarrow\quad \frac{1}{\sqrt{2}}\Big(|0\rangle |0\rangle^N+|1\rangle P |\psi_i\rangle\Big) \quad\longrightarrow\quad\frac{1}{\sqrt{2}}\Big(|0\rangle |\psi_j\rangle+|1\rangle P|\psi_i\rangle\Big) \end{equation*}

Où $P$ est l'un des termes de la décomposition de l'hamiltonien $H=\sum P$ et $\text{Prep} \; \psi_i$ , $\text{Prep} \; \psi_j$ sont des opérations contrôlées qui préparent $|\psi_i\rangle$ , $|\psi_j\rangle$ les vecteurs de l'espace de Krylov unitaire, avec $|\psi_k\rangle = e^{-i H k dt } \vert \psi \rangle = e^{-i H k dt } U_{\psi} \vert 0 \rangle^N$ . Pour mesurer $X$ , appliquez d'abord $H$ ...

\begin{equation*} \longrightarrow\quad\frac{1}{2}|0\rangle\Big( |\psi_j\rangle + P|\psi_i\rangle\Big) + \frac{1}{2}|1\rangle\Big(|\psi_j\rangle - P|\psi_i\rangle\Big) \end{equation*}

... puis mesurer :

\begin{equation*} \begin{split} \Rightarrow\quad\langle X\rangle &= \frac{1}{4}\Bigg(\Big\|| \psi_j\rangle + P|\psi_i\rangle \Big\|^2-\Big\||\psi_j\rangle - P|\psi_i\rangle\Big\|^2\Bigg) \\ &= \text{Re}\Big[\langle\psi_j| P|\psi_i\rangle\Big]. \end{split} \end{equation*}

D'après l'identité $|a + b\|^2 = \langle a + b | a + b \rangle = \|a\|^2 + \|b\|^2 + 2\text{Re}\langle a | b \rangle$ . De même, en mesurant $Y$ , on obtient

\begin{equation*} \langle Y\rangle = \text{Im}\Big[\langle\psi_j| P|\psi_i\rangle\Big]. \end{equation*}

## Create the time-evo op circuit
evol_gate = PauliEvolutionGate(
    H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)
)

## Create the time-evo op dagger circuit
evol_gate_d = PauliEvolutionGate(
    H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)
)
evol_gate_d = evol_gate_d.inverse()

# Put pieces together
qc_reg = QuantumRegister(n_qubits)
qc_temp = QuantumCircuit(qc_reg)
qc_temp.compose(qc_state_prep, inplace=True)
for _ in range(num_trotter_steps):
    qc_temp.append(evol_gate, qargs=qc_reg)
for _ in range(num_trotter_steps):
    qc_temp.append(evol_gate_d, qargs=qc_reg)
qc_temp.compose(qc_state_prep.inverse(), inplace=True)

# Create controlled version of the circuit
controlled_U = qc_temp.to_gate().control(1)

# Create hadamard test circuit for real part
qr = QuantumRegister(n_qubits + 1)
qc_real = QuantumCircuit(qr)
qc_real.h(0)
qc_real.append(controlled_U, list(range(n_qubits + 1)))
qc_real.h(0)

print(
    "Circuit for calculating the real part of the overlap in S via Hadamard test"
)
qc_real.draw("mpl", fold=-1, scale=0.5)

Output:

Circuit for calculating the real part of the overlap in S via Hadamard test

Le circuit de test de Hadamard peut être un circuit profond une fois que nous le décomposons en portes natives (ce qui augmentera encore si nous tenons compte de la topologie de l'appareil)

print(
    "Number of layers of 2Q operations",
    qc_real.decompose(reps=2).depth(lambda x: x[0].num_qubits == 2),
)

Output:

Number of layers of 2Q operations 112753

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

Test de Hadamard efficace

Nous pouvons optimiser les circuits profonds pour le test de Hadamard que nous avons obtenu en introduisant certaines approximations et en nous appuyant sur certaines hypothèses concernant l'hamiltonien du modèle. Par exemple, considérons le circuit suivant pour le test de Hadamard :

Supposons que nous puissions calculer classiquement $E_0$ , la valeur propre de $|0\rangle^N$ sous l'hamiltonien $H$ . Cette condition est remplie lorsque l'hamiltonien préserve la symétrie U(1). Bien que cette hypothèse puisse sembler forte, il existe de nombreux cas où l'on peut supposer qu'il existe un état de vide (dans ce cas, il s'agit de l'état $|0\rangle^N$ ) qui n'est pas affecté par l'action de l'hamiltonien. C'est le cas, par exemple, des hamiltoniens de chimie qui décrivent une molécule stable (où le nombre d'électrons est conservé). Étant donné que la porte $\text{Prep} \; \psi$ , prépare l'état de référence souhaité $\ket{psi} = \text{Prep} \; \psi \ket{0} = e^{-i H 0 dt} U_{\psi} \ket{0}$ , par exemple, préparer l'état HF pour la chimie $\text{Prep} \; \psi$ serait un produit de NOT à un seul qubit, de sorte que controlled- $\text{Prep} \; \psi$ est simplement un produit de CNOT. Le circuit ci-dessus met alors en œuvre l'état suivant avant la mesure :

\begin{equation} \begin{split} \ket{0} \ket{0}^N\xrightarrow{H}&\frac{1}{\sqrt{2}} \left( \ket{0}\ket{0}^N+ \ket{1} \ket{0}^N \right)\\ \xrightarrow{\text{1-ctrl-init}}&\frac{1}{\sqrt{2}}\left(|0\rangle|0\rangle^N+|1\rangle|\psi\rangle\right)\\ \xrightarrow{U}&\frac{1}{\sqrt{2}}\left(e^{i\phi}\ket{0}\ket{0}^N+\ket{1} U\ket{\psi}\right)\\ \xrightarrow{\text{0-ctrl-init}}&\frac{1}{\sqrt{2}} \left( e^{i\phi}\ket{0} \ket{\psi} +\ket{1} U\ket{\psi} \right)\\ =&\frac{1}{2} \left( \ket{+}\left(e^{i\phi}\ket{\psi}+U\ket{\psi}\right) +\ket{-}\left(e^{i\phi}\ket{\psi}-U\ket{\psi}\right) \right)\\ =&\frac{1}{2} \left( \ket{+i}\left(e^{i\phi}\ket{\psi}-iU\ket{\psi}\right) +\ket{-i}\left(e^{i\phi}\ket{\psi}+iU\ket{\psi}\right) \right) \end{split} \end{equation}

où nous avons utilisé le déphasage simulable classique $U\ket{0}^N = e^{i\phi}\ket{0}^N$ dans la troisième ligne. Par conséquent, les valeurs attendues sont obtenues comme suit

\begin{equation} \begin{split} \langle X\otimes P\rangle&=\frac{1}{4} \Big( \left(e^{-i\phi}\bra{\psi}+\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}+U\ket{\psi}\right) \\ &\qquad-\left(e^{-i\phi}\bra{\psi}-\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}-U\ket{\psi}\right) \Big)\\ &=\text{Re}\left[e^{-i\phi}\bra{\psi}PU\ket{\psi}\right], \end{split} \end{equation}

\begin{equation} \begin{split} \langle Y\otimes P\rangle&=\frac{1}{4} \Big( \left(e^{-i\phi}\bra{\psi}+i\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}-iU\ket{\psi}\right) \\ &\qquad-\left(e^{-i\phi}\bra{\psi}-i\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}+iU\ket{\psi}\right) \Big)\\ &=\text{Im}\left[e^{-i\phi}\bra{\psi}PU\ket{\psi}\right]. \end{split} \end{equation}

En utilisant ces hypothèses, nous avons pu écrire les valeurs attendues des opérateurs d'intérêt avec moins d'opérations contrôlées. En fait, nous ne devons mettre en œuvre que la préparation contrôlée de l'état $\text{Prep} \; \psi$ et non les évolutions temporelles contrôlées. En reformulant notre calcul comme indiqué ci-dessus, nous pourrons réduire considérablement la profondeur des circuits résultants.

Décomposer l'opérateur d'évolution temporelle avec la décomposition de Trotter

Au lieu d'implémenter exactement l'opérateur d'évolution temporelle, nous pouvons utiliser la décomposition de Trotter pour en implémenter une approximation. En répétant plusieurs fois une décomposition de Trotter d'un certain ordre, nous réduisons encore l'erreur introduite par l'approximation. Dans ce qui suit, nous construisons directement l'implémentation de Trotter de la manière la plus efficace pour le graphe d'interaction de l'hamiltonien que nous considérons (interactions avec les voisins les plus proches uniquement). En pratique, nous insérons les rotations de Pauli $R_{xx}$ , $R_{yy}$ , $R_{zz}$ avec un angle paramétré $t$ qui correspond à la mise en œuvre approximative de $e^{-i (XX + YY + ZZ) t}$ . Étant donné la différence de définition des rotations de Pauli et l'évolution temporelle que nous essayons de mettre en œuvre, nous devrons utiliser le paramètre $2*dt$ pour obtenir une évolution temporelle de $dt$ . En outre, nous inversons l'ordre des opérations pour un nombre impair de répétitions des étapes de Trotter, ce qui est fonctionnellement équivalent mais permet de synthétiser des opérations adjacentes dans une seule unité $SU(2)$ . Cela donne un circuit beaucoup moins profond que celui obtenu en utilisant la fonctionnalité générique PauliEvolutionGate() .

t = Parameter("t")

# Create instruction for rotation about XX+YY-ZZ:
Rxyz_circ = QuantumCircuit(2)
Rxyz_circ.rxx(t, 0, 1)
Rxyz_circ.ryy(t, 0, 1)
Rxyz_circ.rzz(t, 0, 1)
Rxyz_instr = Rxyz_circ.to_instruction(label="RXX+YY+ZZ")

interaction_list = [
    [[i, i + 1] for i in range(0, n_qubits - 1, 2)],
    [[i, i + 1] for i in range(1, n_qubits - 1, 2)],
]  # linear chain

qr = QuantumRegister(n_qubits)
trotter_step_circ = QuantumCircuit(qr)
for i, color in enumerate(interaction_list):
    for interaction in color:
        trotter_step_circ.append(Rxyz_instr, interaction)
    if i < len(interaction_list) - 1:
        trotter_step_circ.barrier()
reverse_trotter_step_circ = trotter_step_circ.reverse_ops()

qc_evol = QuantumCircuit(qr)
for step in range(num_trotter_steps):
    if step % 2 == 0:
        qc_evol = qc_evol.compose(trotter_step_circ)
    else:
        qc_evol = qc_evol.compose(reverse_trotter_step_circ)

qc_evol.decompose().draw("mpl", fold=-1, scale=0.5)

Output:

Utiliser un circuit optimisé pour la préparation d'état

control = 0
excitation = int(n_qubits / 2) + 1
controlled_state_prep = QuantumCircuit(n_qubits + 1)
controlled_state_prep.cx(control, excitation)
controlled_state_prep.draw("mpl", fold=-1, scale=0.5)

Output:

Circuits modèles pour calculer les éléments matriciels de $\tilde{S}$ et $\tilde{H}$ via le test de Hadamard

La seule différence entre les circuits utilisés dans le test de Hadamard sera la phase de l'opérateur d'évolution temporelle et les observables mesurés. Nous pouvons donc préparer un circuit modèle qui représente le circuit générique pour le test de Hadamard, avec des espaces réservés pour les portes qui dépendent de l'opérateur d'évolution temporelle.

# Parameters for the template circuits
parameters = []
for idx in range(1, krylov_dim):
    parameters.append(2 * dt_circ * (idx))

# Create modified hadamard test circuit
qr = QuantumRegister(n_qubits + 1)
qc = QuantumCircuit(qr)
qc.h(0)
qc.compose(controlled_state_prep, list(range(n_qubits + 1)), inplace=True)
qc.barrier()
qc.compose(qc_evol, list(range(1, n_qubits + 1)), inplace=True)
qc.barrier()
qc.x(0)
qc.compose(
    controlled_state_prep.inverse(), list(range(n_qubits + 1)), inplace=True
)
qc.x(0)

qc.decompose().draw("mpl", fold=-1)

Output:

print(
    "The optimized circuit has 2Q gates depth: ",
    qc.decompose().decompose().depth(lambda x: x[0].num_qubits == 2),
)

Output:

The optimized circuit has 2Q gates depth:  74

Nous avons considérablement réduit la profondeur du test de Hadamard en combinant l'approximation de Trotter et les unitaires non contrôlés

Étape 3 : Exécutez à l'aide d' Qiskit primitives

Instanciation du backend et définition des paramètres d'exécution

service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
if (
    "if_else" not in backend.target.operation_names
):  # Needed as "op_name" could be "if_else"
    backend.target.add_instruction(IfElseOp, name="if_else")
print(backend.name)

Transpilation vers un QPU

Tout d'abord, sélectionnons des sous-ensembles de la carte de couplage avec des qubits « performants » (le terme « performant » étant ici assez arbitraire, nous voulons surtout éviter les qubits vraiment peu performants) et créons une nouvelle cible pour la transpilation

target = backend.target
cmap = target.build_coupling_map(filter_idle_qubits=True)
cmap_list = list(cmap.get_edges())

cust_cmap_list = copy.deepcopy(cmap_list)
for q in range(target.num_qubits):
    meas_err = target["measure"][(q,)].error
    t2 = target.qubit_properties[q].t2 * 1e6
    if meas_err > 0.02 or t2 < 100:
        for q_pair in cmap_list:
            if q in q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except:
                    continue

for q in cmap_list:
    op_name = list(target.operation_names_for_qargs(q))[0]
    twoq_gate_err = target[f"{op_name}"][q].error
    if twoq_gate_err > 0.005:
        for q_pair in cmap_list:
            if q == q_pair:
                try:
                    cust_cmap_list.remove(q)
                except:
                    continue


cust_cmap = CouplingMap(cust_cmap_list)
cust_target = Target.from_configuration(
    basis_gates=backend.configuration().basis_gates,
    coupling_map=cust_cmap,
)

Transpilez ensuite le circuit virtuel vers la meilleure disposition physique dans cette nouvelle cible

basis_gates = list(target.operation_names)
pm = generate_preset_pass_manager(
    optimization_level=3,
    target=cust_target,
    basis_gates=basis_gates,
)

qc_trans = pm.run(qc)

print("depth", qc_trans.depth(lambda x: x[0].num_qubits == 2))
print("num 2q ops", qc_trans.count_ops())
print(
    "physical qubits",
    sorted(
        [
            idx
            for idx, qb in qc_trans.layout.initial_layout.get_physical_bits().items()
            if qb._register.name != "ancilla"
        ]
    ),
)

Output:

depth 52
num 2q ops OrderedDict([('rz', 2058), ('sx', 1703), ('cz', 728), ('x', 84), ('barrier', 8)])
physical qubits [91, 92, 93, 94, 95, 98, 99, 108, 109, 110, 111, 113, 114, 115, 119, 127, 132, 133, 134, 135, 137, 139, 147, 148, 149, 150, 151, 152, 153, 154, 155]

Créer des PUB pour exécution avec Estimator

# Define observables to measure for S
observable_S_real = "I" * (n_qubits) + "X"
observable_S_imag = "I" * (n_qubits) + "Y"

observable_op_real = SparsePauliOp(
    observable_S_real
)  # define a sparse pauli operator for the observable
observable_op_imag = SparsePauliOp(observable_S_imag)

layout = qc_trans.layout  # get layout of transpiled circuit
observable_op_real = observable_op_real.apply_layout(
    layout
)  # apply physical layout to the observable
observable_op_imag = observable_op_imag.apply_layout(layout)
observable_S_real = (
    observable_op_real.paulis.to_labels()
)  # get the label of the physical observable
observable_S_imag = observable_op_imag.paulis.to_labels()

observables_S = [[observable_S_real], [observable_S_imag]]


# Define observables to measure for H
# Hamiltonian terms to measure
observable_list = []
for pauli, coeff in zip(H_op.paulis, H_op.coeffs):
    # print(pauli)
    observable_H_real = pauli[::-1].to_label() + "X"
    observable_H_imag = pauli[::-1].to_label() + "Y"
    observable_list.append([observable_H_real])
    observable_list.append([observable_H_imag])

layout = qc_trans.layout

observable_trans_list = []
for observable in observable_list:
    observable_op = SparsePauliOp(observable)
    observable_op = observable_op.apply_layout(layout)
    observable_trans_list.append([observable_op.paulis.to_labels()])

observables_H = observable_trans_list


# Define a sweep over parameter values
params = np.vstack(parameters).T


# Estimate the expectation value for all combinations of
# observables and parameter values, where the pub result will have
# shape (# observables, # parameter values).
pub = (qc_trans, observables_S + observables_H, params)

Circuits de course

Les circuits pour $t=0$ sont calculables de manière classique

qc_cliff = qc.assign_parameters({t: 0})


# Get expectation values from experiment
S_expval_real = StabilizerState(qc_cliff).expectation_value(
    Pauli("I" * (n_qubits) + "X")
)
S_expval_imag = StabilizerState(qc_cliff).expectation_value(
    Pauli("I" * (n_qubits) + "Y")
)

# Get expectation values
S_expval = S_expval_real + 1j * S_expval_imag

H_expval = 0
for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):
    # Get expectation values from experiment
    expval_real = StabilizerState(qc_cliff).expectation_value(
        Pauli(pauli[::-1].to_label() + "X")
    )
    expval_imag = StabilizerState(qc_cliff).expectation_value(
        Pauli(pauli[::-1].to_label() + "Y")
    )
    expval = expval_real + 1j * expval_imag

    # Fill-in matrix elements
    H_expval += coeff * expval


print(H_expval)

Output:

(25+0j)

Exécuter des circuits pour $S$ et $\tilde{H}$ à l'aide d'Estimator

# Experiment options
num_randomizations = 300
num_randomizations_learning = 30
shots_per_randomization = 100
noise_factors = [1, 1.2, 1.4]
learning_pair_depths = [0, 4, 24, 48]


experimental_opts = {}
experimental_opts["resilience"] = {
    "measure_mitigation": True,
    "measure_noise_learning": {
        "num_randomizations": num_randomizations_learning,
        "shots_per_randomization": shots_per_randomization,
    },
    "zne_mitigation": True,
    "zne": {"noise_factors": noise_factors},
    "layer_noise_learning": {
        "max_layers_to_learn": 10,
        "layer_pair_depths": learning_pair_depths,
        "shots_per_randomization": shots_per_randomization,
        "num_randomizations": num_randomizations_learning,
    },
    "zne": {
        "amplifier": "pea",
        "extrapolated_noise_factors": [0] + noise_factors,
    },
}
experimental_opts["twirling"] = {
    "num_randomizations": num_randomizations,
    "shots_per_randomization": shots_per_randomization,
    "strategy": "all",
}

estimator = Estimator(mode=backend, options=experimental_opts)


job = estimator.run([pub])

Étape 4 : Post-traitement et restitution du résultat dans le format classique souhaité

results = job.result()[0]

Calculer les matrices hamiltoniennes et de chevauchement effectives

Calculer d'abord la phase accumulée par l'état $\vert 0 \rangle$ au cours de l'évolution temporelle incontrôlée

prefactors = [
    np.exp(-1j * sum([c for p, c in H_op.to_list() if "Z" in p]) * i * dt)
    for i in range(1, krylov_dim)
]

Une fois que nous avons les résultats de l'exécution des circuits, nous pouvons post-traiter les données pour calculer les éléments de la matrice de $S$

# Assemble S, the overlap matrix of dimension D:
S_first_row = np.zeros(krylov_dim, dtype=complex)
S_first_row[0] = 1 + 0j

# Add in ancilla-only measurements:
for i in range(krylov_dim - 1):
    # Get expectation values from experiment
    expval_real = results.data.evs[0][0][
        i
    ]  # automatic extrapolated evs if ZNE is used
    expval_imag = results.data.evs[1][0][
        i
    ]  # automatic extrapolated evs if ZNE is used

    # Get expectation values
    expval = expval_real + 1j * expval_imag
    S_first_row[i + 1] += prefactors[i] * expval

S_first_row_list = S_first_row.tolist()  # for saving purposes


S_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)

# Distribute entries from first row across matrix:
for i, j in it.product(range(krylov_dim), repeat=2):
    if i >= j:
        S_circ[j, i] = S_first_row[i - j]
    else:
        S_circ[j, i] = np.conj(S_first_row[j - i])

Matrix(S_circ)

Output:

\displaystyle \left[\begin{matrix}1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i & -0.180546747798251 - 0.492624093654174 i & 0.0012070853532697 + 0.312052218182462 i\\-0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i & -0.180546747798251 - 0.492624093654174 i\\0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i\\-0.180546747798251 + 0.492624093654174 i & 0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i\\0.0012070853532697 - 0.312052218182462 i & -0.180546747798251 + 0.492624093654174 i & 0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0\end{matrix}\right]

Et les éléments de la matrice de $\tilde{H}$

# Assemble S, the overlap matrix of dimension D:
H_first_row = np.zeros(krylov_dim, dtype=complex)
H_first_row[0] = H_expval

for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):
    # Add in ancilla-only measurements:
    for i in range(krylov_dim - 1):
        # Get expectation values from experiment
        expval_real = results.data.evs[2 + 2 * obs_idx][0][
            i
        ]  # automatic extrapolated evs if ZNE is used
        expval_imag = results.data.evs[2 + 2 * obs_idx + 1][0][
            i
        ]  # automatic extrapolated evs if ZNE is used

        # Get expectation values
        expval = expval_real + 1j * expval_imag
        H_first_row[i + 1] += prefactors[i] * coeff * expval

H_first_row_list = H_first_row.tolist()

H_eff_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)

# Distribute entries from first row across matrix:
for i, j in it.product(range(krylov_dim), repeat=2):
    if i >= j:
        H_eff_circ[j, i] = H_first_row[i - j]
    else:
        H_eff_circ[j, i] = np.conj(H_first_row[j - i])

Matrix(H_eff_circ)

Output:

\displaystyle \left[\begin{matrix}25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i & -5.15587257589417 - 8.88280836036843 i & 1.98818301405581 + 5.8897614762563 i\\-14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i & -5.15587257589417 - 8.88280836036843 i\\10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i\\-5.15587257589417 + 8.88280836036843 i & 10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i\\1.98818301405581 - 5.8897614762563 i & -5.15587257589417 + 8.88280836036843 i & 10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0\end{matrix}\right]

Enfin, nous pouvons résoudre le problème des valeurs propres généralisées pour $\tilde{H}$ :

$\tilde{H} \vec{c} = c S \vec{c}$

et obtenir une estimation de l'énergie de l'état fondamental $c_{min}$

gnd_en_circ_est_list = []
for d in range(1, krylov_dim + 1):
    # Solve generalized eigenvalue problem for different size of the Krylov space
    gnd_en_circ_est = solve_regularized_gen_eig(
        H_eff_circ[:d, :d], S_circ[:d, :d], threshold=9e-1
    )
    gnd_en_circ_est_list.append(gnd_en_circ_est)
    print("The estimated ground state energy is: ", gnd_en_circ_est)

Output:

The estimated ground state energy is:  25.0
The estimated ground state energy is:  22.572154819954875
The estimated ground state energy is:  21.691509219286587
The estimated ground state energy is:  21.23882298756386
The estimated ground state energy is:  20.965499325470294

Pour un secteur à une seule particule, nous pouvons calculer efficacement l'état fondamental de ce secteur de l'hamiltonien de manière classique

gs_en = single_particle_gs(H_op, n_qubits)

Output:

n_sys_qubits 30
n_exc 1 , subspace dimension 31
single particle ground state energy:  21.021912418526906

plt.plot(
    range(1, krylov_dim + 1),
    gnd_en_circ_est_list,
    color="blue",
    linestyle="-.",
    label="KQD estimate",
)
plt.plot(
    range(1, krylov_dim + 1),
    [gs_en] * krylov_dim,
    color="red",
    linestyle="-",
    label="exact",
)
plt.xticks(range(1, krylov_dim + 1), range(1, krylov_dim + 1))
plt.legend()
plt.xlabel("Krylov space dimension")
plt.ylabel("Energy")
plt.title(
    "Estimating Ground state energy with Krylov Quantum Diagonalization"
)
plt.show()

Output:

Annexe : Sous-espace de Krylov à partir d'évolutions en temps réel

L'espace de Krylov unitaire est défini comme suit

\mathcal{K}_U(H, |\psi\rangle) = \text{span}\left\{ |\psi\rangle, e^{-iH\,dt} |\psi\rangle, \dots, e^{-irH\,dt} |\psi\rangle \right\}

pour un certain pas de temps $dt$ que nous déterminerons plus tard. Supposons temporairement que $r$ est pair : définissons alors $d=r/2$ . Remarquez que lorsque nous projetons l'hamiltonien dans l'espace de Krylov ci-dessus, il est indiscernable de l'espace de Krylov

\mathcal{K}_U(H, |\psi\rangle) = \text{span}\left\{ e^{i\,d\,H\,dt}|\psi\rangle, e^{i(d-1)H\,dt} |\psi\rangle, \dots, e^{-i(d-1)H\,dt} |\psi\rangle, e^{-i\,d\,H\,dt} |\psi\rangle \right\},

c'est-à-dire lorsque toutes les évolutions temporelles sont décalées vers l'arrière de $d$ pas de temps. La raison pour laquelle il n'est pas possible de les distinguer est que les éléments de la matrice

\tilde{H}_{j,k} = \langle\psi|e^{i\,j\,H\,dt}He^{-i\,k\,H\,dt}|\psi\rangle=\langle\psi|He^{i(j-k)H\,dt}|\psi\rangle

sont invariants en cas de décalage global du temps d'évolution, puisque les évolutions temporelles commutent avec l'hamiltonien. Pour les $r$ impairs, nous pouvons utiliser l'analyse pour les $r-1$ .

Nous voulons montrer que quelque part dans cet espace de Krylov, il est garanti qu'il existe un état de basse énergie. Nous le faisons au moyen du résultat suivant, qui est dérivé du théorème 3.1 dans [3] :

Affirmation 1 : il existe une fonction $f$ telle que pour les énergies $E$ dans le domaine spectral de l'hamiltonien (c'est-à-dire entre l'énergie de l'état fondamental et l'énergie maximale),...

$f(E_0)=1$
$|f(E)|\le2\left(1 + \delta\right)^{-d}$ pour toutes les valeurs de $E$ qui se situent à $\ge\delta$ de $E_0$ , c'est-à-dire qu'elle est exponentiellement supprimée
$f(E)$ est une combinaison linéaire de $e^{ijE\,dt}$ pour $j=-d,-d+1,...,d-1,d$

Nous donnons une preuve ci-dessous, mais elle peut être ignorée sans risque, à moins que l'on ne veuille comprendre l'argument complet et rigoureux. Pour l'instant, nous nous concentrons sur les implications de l'affirmation ci-dessus. En vertu de la propriété 3 ci-dessus, nous pouvons voir que l'espace de Krylov décalé ci-dessus contient l'état $f(H)|\psi\rangle$ . Il s'agit de notre état de basse énergie. Pour comprendre pourquoi, il faut écrire $|\psi\rangle$ dans la base propre de l'énergie :

|\psi\rangle = \sum_{k=0}^{N}\gamma_k|E_k\rangle,

où $|E_k\rangle$ est le kème état propre énergétique et $\gamma_k$ est son amplitude dans l'état initial $|\psi\rangle$ . Exprimé en termes de ceci, $f(H)|\psi\rangle$ est donné par

f(H)|\psi\rangle = \sum_{k=0}^{N}\gamma_kf(E_k)|E_k\rangle,

en utilisant le fait que l'on peut remplacer $H$ par $E_k$ lorsqu'il agit sur l'état propre $|E_k\rangle$ . L'erreur énergétique de cet état est donc

\text{energy error} = \frac{\langle\psi|f(H)(H-E_0)f(H)|\psi\rangle}{\langle\psi|f(H)^2|\psi\rangle}

= \frac{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2}.

Pour transformer ceci en une borne supérieure plus facile à comprendre, nous séparons d'abord la somme du numérateur en termes avec $E_k-E_0\le\delta$ et en termes avec $E_k-E_0>\delta$ :

\text{energy error} = \frac{\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} + \frac{\sum_{E_k> E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2}.

Nous pouvons limiter le premier terme par $\delta$ ,

\frac{\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} < \frac{\delta\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} \le \delta,

où la première étape suit parce que $E_k-E_0\le\delta$ pour chaque $E_k$ dans la somme, et la deuxième étape suit parce que la somme dans le numérateur est un sous-ensemble de la somme dans le dénominateur. Pour le second terme, on commence par abaisser le dénominateur par $|\gamma_0|^2$ , puisque $f(E_0)^2=1$ : en additionnant le tout, on obtient

\text{energy error} \le \delta + \frac{1}{|\gamma_0|^2}\sum_{E_k>E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0).

Pour simplifier ce qui reste, remarquez que pour tous ces $E_k$ , par la définition de $f$ nous savons que $f(E_k)^2 \le 4\left(1 + \delta\right)^{-2d}$ . De plus, la borne supérieure de $E_k-E_0<2\|H\|$ et la borne supérieure de $\sum_{E_k>E_0+\delta}|\gamma_k|^2<1$ donnent

\text{energy error} \le \delta + \frac{8}{|\gamma_0|^2}\|H\|\left(1 + \delta\right)^{-2d}.

Ceci est valable pour n'importe quel $\delta>0$ , donc si nous fixons $\delta$ égal à notre erreur cible, alors la limite d'erreur ci-dessus converge vers cela de manière exponentielle avec la dimension de Krylov $2d=r$ . Notez également que si $\delta<E_1-E_0$ , le terme $\delta$ disparaît entièrement dans la limite ci-dessus.

Pour compléter l'argument, nous notons tout d'abord que ce qui précède n'est que l'erreur énergétique de l'état particulier $f(H)|\psi\rangle$ , plutôt que l'erreur énergétique de l'état de plus faible énergie dans l'espace de Krylov. Toutefois, en vertu du principe variationnel (Rayleigh-Ritz), l'erreur énergétique de l'état de plus faible énergie dans l'espace de Krylov est limitée par l'erreur énergétique de tout état dans l'espace de Krylov, de sorte que ce qui précède est également une limite supérieure de l'erreur énergétique de l'état de plus faible énergie, c'est-à-dire la sortie de l'algorithme de diagonalisation quantique de Krylov.

Il est possible d'effectuer une analyse similaire à la précédente en tenant compte du bruit et de la procédure de seuillage décrite dans le manuel. Voir [2] et [4] pour cette analyse.

Annexe : preuve de la revendication 1

Ce qui suit est en grande partie dérivé de [3], Théorème 3.1: Soit $0 < a < b$ et $\Pi^*_d$ l'espace des polynômes résiduels (polynômes dont la valeur en 0 est 1) de degré au plus $d$ . La solution de

\beta(a, b, d) = \min_{p \in \Pi^*_d} \max_{x \in [a, b]} |p(x)| \quad

est

p^*(x) = \frac{T_d\left(\frac{b + a - 2x}{b - a}\right)}{T_d\left(\frac{b + a}{b - a}\right)}, \quad

et la valeur minimale correspondante est

\beta(a, b, d) = T_d^{-1}\left(\frac{b + a}{b - a}\right).

Nous voulons convertir cette fonction en une fonction qui peut être exprimée naturellement en termes d'exponentielles complexes, car ce sont les évolutions temporelles réelles qui génèrent l'espace de Krylov quantique. Pour ce faire, il est commode d'introduire la transformation suivante des énergies dans le domaine spectral de l'hamiltonien en nombres dans l'intervalle $[0,1]$ : définir

g(E) = \frac{1-\cos\big((E-E_0)dt\big)}{2},

où $dt$ est un pas de temps tel que $-\pi < E_0dt < E_\text{max}dt < \pi$ . Remarquez que $g(E_0)=0$ et $g(E)$ croissent au fur et à mesure que $E$ s'éloigne de $E_0$ .

En utilisant maintenant le polynôme $p^*(x)$ avec les paramètres a, b, d fixés à $a = g(E_0 + \delta)$ , $b = 1$ , et d = int( r/2 ), nous définissons la fonction :

f(E) = p^* \left( g(E) \right) = \frac{T_d\left(1 + 2\frac{\cos\big((E-E_0)dt\big) - \cos\big(\delta\,dt\big)}{1 +\cos\big(\delta\,dt\big)}\right)}{T_d\left(1 + 2\frac{1-\cos\big(\delta\,dt\big)}{1 + \cos\big(\delta\,dt\big)}\right)}

où $E_0$ est l'énergie de l'état fondamental. Nous pouvons voir en insérant $\cos(x)=\frac{e^{ix}+e^{-ix}}{2}$ que $f(E)$ est un polynôme trigonométrique de degré $d$ , c'est-à-dire une combinaison linéaire de $e^{ijE\,dt}$ pour $j=-d,-d+1,...,d-1,d$ . De plus, d'après la définition de $p^*(x)$ ci-dessus, nous avons que $f(E_0)=p(0)=1$ et pour tout $E$ dans le domaine spectral tel que $\vert E-E_0 \vert > \delta$ nous avons

|f(E)| \le \beta(a, b, d) = T_d^{-1}\left(1 + 2\frac{1-\cos\big(\delta\,dt\big)}{1 + \cos\big(\delta\,dt\big)}\right)

\leq 2\left(1 + \delta\right)^{-d} = 2\left(1 + \delta\right)^{-\lfloor k/2\rfloor}.

Références

[1] N. Yoshioka, M. Amico, W. Kirby et al. "Diagonalization of large many-body Hamiltonians on a quantum processor". arXiv:2407.14431

[2] Ethan N. Epperly, Lin Lin et Yuji Nakatsukasa. "Une théorie de la diagonalisation du sous-espace quantique". SIAM Journal on Matrix Analysis and Applications 43, 1263-1290 (2022).

[3] Å. Björck. "Méthodes numériques dans les calculs matriciels". Textes en mathématiques appliquées. Springer International Publishing. (2014).

[4] William Kirby. "Analyse des algorithmes de Krylov quantiques avec erreurs". Quantum 8, 1457 (2024).

Enquête tutorielle

Veuillez répondre à cette courte enquête pour nous faire part de vos commentaires sur ce didacticiel. Vos commentaires nous aideront à améliorer nos offres de contenu et l'expérience des utilisateurs.

Lien vers l'enquête

Cette page a-t-elle été utile ?

Signaler un bogue, une coquille ou proposer du contenu sur GitHub.

Diagonalisation quantique de Krylov des hamiltoniens de réseau

Arrière-plan

Exigences

Configuration

Étape 1 : Mettre en correspondance les entrées classiques avec un problème quantique

L'espace de Krylov

Algorithme de diagonalisation quantique de Krylov

Définir Hamiltonien

Définir les paramètres de l'algorithme

Préparation de l'État

Évolution dans le temps

test de Hadamard

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

Test de Hadamard efficace

Décomposer l'opérateur d'évolution temporelle avec la décomposition de Trotter

Utiliser un circuit optimisé pour la préparation d'état

Circuits modèles pour calculer les éléments matriciels de S~\tilde{S}S~ et H~\tilde{H}H~ via le test de Hadamard

Étape 3 : Exécutez à l'aide d' Qiskit primitives

Transpilation vers un QPU

Créer des PUB pour exécution avec Estimator

Circuits de course

Étape 4 : Post-traitement et restitution du résultat dans le format classique souhaité

Calculer les matrices hamiltoniennes et de chevauchement effectives

Annexe : Sous-espace de Krylov à partir d'évolutions en temps réel

Annexe : preuve de la revendication 1

Références

Enquête tutorielle

Circuits modèles pour calculer les éléments matriciels de $\tilde{S}$ et $\tilde{H}$ via le test de Hadamard