Diagonalisation quantique de Krylov

Dans cette leçon sur la diagonalisation quantique de Krylov (KQD), nous répondrons aux questions suivantes :

• Qu'est-ce que la méthode de Krylov, en général?
• Pourquoi la méthode de Krylov fonctionne-t-elle et dans quelles conditions?
• Quel est le rôle de l'informatique quantique?

La partie quantique des calculs est basée en grande partie sur les travaux de Ref [1].

La vidéo ci-dessous donne un aperçu des méthodes de Krylov dans l'informatique classique, motive leur utilisation et explique comment l'informatique quantique peut jouer un rôle dans ce domaine. Le texte suivant offre plus de détails et met en œuvre une méthode de Krylov à la fois de manière classique et en utilisant un ordinateur quantique.

1. Introduction aux méthodes de Krylov

Une méthode de sous-espace de Krylov peut se référer à n'importe laquelle des méthodes construites autour de ce qu'on appelle le sous-espace de Krylov. Un examen complet de ces derniers dépasse le cadre de cette leçon, mais les références [2-4] peuvent toutes fournir des informations plus détaillées. Nous nous concentrerons ici sur ce qu'est un sous-espace de Krylov, comment et pourquoi il est utile pour résoudre les problèmes de valeurs propres, et enfin comment il peut être mis en œuvre sur un ordinateur quantique.

Définition : Étant donné une matrice symétrique semi-définie positive $N\times N$ $A$, 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\leq N$, avec un vecteur de référence $\vert v \rangle$.

Bien que les vecteurs ci-dessus couvrent ce que nous appelons un sous-espace de Krylov, il n'y a aucune raison de penser qu'ils seront orthogonaux. On utilise souvent un processus d'orthonormalisation itératif similaire à l' orthogonalisation de Gram-Schmidt. Ici, le processus est légèrement différent puisque chaque nouveau vecteur est rendu orthogonal aux autres au fur et à mesure qu'il est généré. Dans ce contexte, on parle d' itération d'Arnoldi. A partir du vecteur initial $|v\rangle$, on génère le vecteur suivant $A|v\rangle$, puis on s'assure que ce second vecteur est orthogonal au premier en soustrayant sa projection sur $|v\rangle$. C'est-à-dire

On voit aisément que $|v_0\rangle \perp |v_1\rangle$, puisque

Nous procédons de la même manière pour le vecteur suivant, en nous assurant qu'il est orthogonal aux deux vecteurs précédents :

Si nous répétons ce processus pour tous les vecteurs $r$, nous obtenons une base orthonormée complète pour un espace de Krylov. Notez que le processus d'orthogonalisation ici produira zéro une fois $r>m$, puisque $m$ vecteur orthogonal couvre nécessairement tout l'espace. Le processus aboutira également à zéro si un vecteur est un vecteur propre de $A$ puisque tous les vecteurs suivants seront des multiples de ce vecteur.

1.1 Un exemple simple : Krylov à la main

Passons en revue la génération d'un sous-espace de Krylov sur une matrice trivialement petite, afin de voir le processus. Nous commençons par une matrice initiale $A$ qui nous intéresse :

Pour ce petit exemple, nous pouvons déterminer les vecteurs propres et les valeurs propres facilement, même à la main. Nous présentons ici la solution numérique.

# One might use linalg.eigh here, but later matrices may not be Hermitian. So we use
# linalg.eig in this lesson.

import numpy as np

A = np.array([[4, -1, 0], [-1, 4, -1], [0, -1, 4]])
eigenvalues, eigenvectors = np.linalg.eig(A)
print("The eigenvalues are ", eigenvalues)
print("The eigenvectors are ", eigenvectors)

Output:

The eigenvalues are  [2.58578644 4.         5.41421356]
The eigenvectors are  [[ 5.00000000e-01 -7.07106781e-01  5.00000000e-01]
 [ 7.07106781e-01  1.37464400e-16 -7.07106781e-01]
 [ 5.00000000e-01  7.07106781e-01  5.00000000e-01]]

Nous les enregistrons ici pour une comparaison ultérieure :

Nous aimerions étudier comment ce processus fonctionne (ou échoue) lorsque nous augmentons la dimension de notre sous-espace de Krylov, $r$. À cette fin, nous allons appliquer ce processus :

• Générer un sous-espace de l'espace vectoriel complet à partir d'un vecteur choisi aléatoirement $|v\rangle$ (l'appeler $|v_0\rangle$ s'il est déjà normalisé, comme ci-dessus).
• Projetez la matrice complète $A$ sur ce sous-espace et trouvez les valeurs propres de cette matrice projetée $\tilde{A}$.
• Augmenter la taille du sous-espace en générant davantage de vecteurs, en veillant à ce qu'ils soient orthonormés, à l'aide d'un processus similaire à l'orthogonalisation de Gram-Schmidt.
• Projetez $A$ sur le sous-espace le plus grand et trouvez les valeurs propres de la matrice résultante, $\tilde{A}$.
• Répétez cette opération jusqu'à ce que les valeurs propres convergent (ou, dans le cas présent, jusqu'à ce que vous ayez généré des vecteurs couvrant tout l'espace vectoriel de la matrice originale $A$ ).

Une mise en œuvre normale de la méthode de Krylov n'aurait pas besoin de résoudre le problème des valeurs propres de la matrice projetée sur chaque sous-espace de Krylov au fur et à mesure de sa construction. Vous pouvez construire le sous-espace de la dimension souhaitée, projeter la matrice sur ce sous-espace et diagonaliser la matrice projetée. La projection et la diagonalisation à chaque dimension du sous-espace ne sont effectuées que pour vérifier la convergence.

$r=1$ des dimensions :

Nous choisissons un vecteur aléatoire, disons

S'il n'est pas déjà normalisé, normalisez-le.

Nous projetons maintenant notre matrice $A$ sur le sous-espace de ce seul vecteur :

Il s'agit de la projection de la matrice sur notre sous-espace de Krylov lorsqu'il ne contient qu'un seul vecteur, $|v_0\rangle$. La valeur propre de cette matrice est trivialement 4. Nous pouvons considérer cela comme notre estimation d'ordre zéro des valeurs propres (dans ce cas, une seule) de $A$. Bien qu'il s'agisse d'une estimation médiocre, c'est l'ordre de grandeur correct.

$r=2$ des dimensions :

Nous générons maintenant le vecteur suivant dans notre sous-espace en opérant avec $A$ sur le vecteur précédent :

Nous soustrayons maintenant la projection de ce vecteur sur notre vecteur précédent afin de garantir l'orthogonalité.

S'il n'est pas déjà normalisé, normalisez-le. Dans ce cas, le vecteur était déjà normalisé, donc

Nous projetons maintenant notre matrice A sur le sous-espace de ces deux vecteurs :

Il nous reste à déterminer les valeurs propres de cette matrice. Mais cette matrice est légèrement plus petite que la matrice complète. Dans les problèmes impliquant de très grandes matrices, il peut être très avantageux de travailler avec ce sous-espace plus petit.

Bien qu'il ne s'agisse pas encore d'une bonne estimation, elle est meilleure que l'estimation de second ordre. Nous procéderons à une nouvelle itération pour nous assurer que le processus est clair. Cependant, cela réduit l'intérêt de la méthode, puisque nous finirons par diagonaliser une matrice 3x3 lors de l'itération suivante, ce qui signifie que nous n'avons pas gagné de temps ni de puissance de calcul.

$r=3$ des dimensions :

Nous générons maintenant le vecteur suivant dans notre sous-espace en opérant avec A sur le vecteur précédent :

Nous soustrayons maintenant la projection de ce vecteur sur nos deux vecteurs précédents afin de garantir l'orthogonalité.

S'il n'est pas déjà normalisé, normalisez-le. Dans ce cas, le vecteur était déjà normalisé, donc

Nous projetons maintenant notre matrice $A$ sur le sous-espace de ces vecteurs :

Nous déterminons maintenant les valeurs propres :

Ces valeurs propres sont exactement les valeurs propres de la matrice originale $A$. Cela doit être le cas, puisque nous avons étendu notre sous-espace de Krylov à l'ensemble de l'espace vectoriel de la matrice originale $A$.

Dans cet exemple, la méthode de Krylov ne semble pas particulièrement plus facile que la diagonalisation directe. En effet, comme nous le verrons dans les sections suivantes, la méthode de Krylov n'est avantageuse qu'à partir d'une certaine dimension de la matrice; ceci a pour but de nous aider à résoudre les problèmes de valeurs propres/vecteurs propres de matrices extrêmement grandes.

C'est le seul exemple que nous montrerons "à la main", mais la section 2 ci-dessous présente des exemples de calcul.

Clarification des termes

Une idée fausse très répandue est qu'il n'existe qu'un seul sous-espace de Krylov pour un problème donné. Mais bien sûr, comme il existe de nombreux vecteurs initiaux auxquels notre matrice peut être appliquée, il existe de nombreux sous-espaces de Krylov possibles. Nous n'utiliserons l'expression "le sous-espace de Krylov" que pour faire référence à un sous-espace de Krylov spécifique déjà défini pour un exemple particulier. Pour les approches générales de résolution de problèmes, nous nous référerons à "a Krylov subspace " (sous-espace de Krylov). Une dernière précision : il est possible de se référer à un " espace de Krylov". Il est souvent appelé " sous-espace de Krylov" en raison de son utilisation dans le contexte de la projection de matrices d'un espace initial dans un sous-espace. Dans ce contexte, nous l'appellerons le plus souvent sous-espace.

Vérifiez votre compréhension

Expliquez pourquoi il n'est pas (a) utile et (b) possible d'étendre la dimension du sous-espace de Krylov $r$ au-delà de la dimension $N$ de la matrice concernée.

• Réponse

  (a) Étant donné que nous orthonormalisons les vecteurs au fur et à mesure que nous les générons, un ensemble d’ $N$ s de ces vecteurs formera une base complète, ce qui signifie qu’une combinaison linéaire de ceux-ci permet de créer n’importe quel vecteur de l’espace.

  (b) Le processus d'orthogonalisation consiste à soustraire la projection d'un nouveau vecteur sur tous les vecteurs précédents. Si tous les vecteurs précédents génèrent l'espace vectoriel entier, alors la soustraction des projections sur le sous-espace entier nous donnera toujours un vecteur nul.

Imaginons qu'un collègue chercheur fasse une démonstration de la méthode de Krylov appliquée à une petite matrice fictive. Y a-t-il un problème avec leur choix de matrice $A$ et de vecteur initial $|\psi\rangle$?

et

• Réponse

  Votre collègue a accidentellement choisi un vecteur propre pour son vecteur initial. Agir avec la matrice sur le vecteur initial renverra simplement le même vecteur, mis à l'échelle par la valeur propre. Cela ne génère pas de sous-espace de dimension croissante. Conseillez à votre collègue de choisir un vecteur initial différent, en vous assurant qu'il ne s'agit pas d'un vecteur propre.

Appliquez la méthode de Krylov à la matrice donnée, en choisissant un nouveau vecteur initial approprié. Notez les estimations de la valeur propre minimale aux ordres 0 et 1 de votre sous-espace de Krylov.

• Réponse

  De nombreuses réponses sont possibles en fonction du choix du vecteur initial. Nous choisirons :

  $$|v_0\rangle=\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix}.$$

  Pour obtenir $|v_1\rangle$, on applique une fois $A$ à $|v_0\rangle$, puis on rend $|v_1\rangle$ orthogonal à $|v_0\rangle$.

  $$A|v_0\rangle=\begin{pmatrix}1&1&0\\1&1&1\\0&1&1\end{pmatrix}\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix} = \frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix}$$$$A|v_0\rangle - \langle v_0|A|v_0\rangle |v_0\rangle=\frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix} - \frac{1}{\sqrt{3}}\begin{pmatrix}1&1&1\end{pmatrix}\frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix}\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix} = \frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix}-\frac{7}{3}\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix}=\sqrt{\frac{3}{2}}\begin{pmatrix}-1/3\\2/3\\-1/3\end{pmatrix}$$

  Au 0ème ordre, la projection sur notre sous-espace de Krylov est la suivante

  $$\langle v_0|A|v_0\rangle=\frac{1}{\sqrt{3}}\begin{pmatrix}1&1&1\end{pmatrix} \begin{pmatrix}1&1&0\\1&1&1\\0&1&1\end{pmatrix} \frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix} = \frac{7}{3}$$

  Au premier ordre, la projection sur ce sous-espace de Krylov est la suivante

  $$\langle V^1|A|V^1\rangle=\begin{pmatrix}\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}\\-\sqrt{\frac{1}{6}}&\sqrt{\frac{2}{3}}&-\sqrt{\frac{1}{6}}\end{pmatrix} \begin{pmatrix}1&1&0\\1&1&1\\0&1&1\end{pmatrix} \begin{pmatrix}\frac{1}{\sqrt{3}}&-\sqrt{\frac{1}{6}}\\\frac{1}{\sqrt{3}}& \sqrt{\frac{2}{3}} \\ \frac{1}{\sqrt{3}}&-\sqrt{\frac{1}{6}}\end{pmatrix}$$

  Cette opération peut être effectuée à la main, mais elle est plus facile à réaliser avec numpy :

  import numpy as np
  vstar = np.array([[1/np.sqrt(3),1/np.sqrt(3),1/np.sqrt(3)],[-1/np.sqrt(6),np.sqrt(2/3),-1/np.sqrt(6)]]
  )
  A = np.array([[1, 1, 0],
                [1, 1, 1],
                [0, 1, 1]])
  v = np.array([[1/np.sqrt(3),-1/np.sqrt(6)],[1/np.sqrt(3),np.sqrt(2/3)],[1/np.sqrt(3),-1/np.sqrt(6)]])
  proj = vstar@A@v
  print(proj)
  eigenvalues, eigenvectors = np.linalg.eig(proj)
  print("The eigenvalues are ", eigenvalues)
  print("The eigenvectors are ", eigenvectors)

  outputs:

  [[ 2.33333333  0.47140452]
   [ 0.47140452 -0.33333333]]
  The eigenvalues are  [ 2.41421356 -0.41421356]
  The eigenvectors are  [[ 0.98559856 -0.16910198]
   [ 0.16910198  0.98559856]]

  L'estimation de la valeur propre minimale est -0.414.

1.2 Types de méthodes de Krylov

les "méthodes du sous-espace de Krylov" peuvent faire référence à plusieurs techniques itératives utilisées pour résoudre les grands systèmes linéaires et les problèmes de valeurs propres. Toutes ces méthodes ont en commun de construire une solution approximative à partir d'un sous-espace de Krylov

$\mathcal{K}^r(A,|v\rangle ) = \text{span}\{|v\rangle, A|v\rangle, A^2|v\rangle, ..., A^{r-1}|v\rangle\},$

où $|v\rangle$ est l'estimation initiale (voir Ref [5]). Elles diffèrent par la manière dont elles choisissent la meilleure approximation à partir de ce sous-espace, en équilibrant des facteurs tels que le taux de convergence, l'utilisation de la mémoire et le coût de calcul global. L'objectif de cette leçon est d'exploiter l'informatique quantique dans le contexte des méthodes du sous-espace de Krylov; une discussion exhaustive de ces méthodes dépasse le cadre de cette leçon. Les brèves définitions ci-dessous ne sont données qu'à titre indicatif et comprennent quelques références permettant d'approfondir ces méthodes.

La méthode du gradient conjugué (CG) : Cette méthode est utilisée pour résoudre les systèmes linéaires symétriques et définis positifs [6]. Il minimise la norme A de l'erreur à chaque itération, ce qui le rend particulièrement efficace pour les systèmes issus d'EDP elliptiques discrétisées [7]. Nous utiliserons cette approche dans la section suivante pour expliquer pourquoi un sous-espace de Krylov serait un sous-espace efficace pour rechercher des solutions améliorées aux systèmes linéaires.

La méthode du résidu minimal généralisé (GMRES) : Cette méthode est conçue pour résoudre des systèmes linéaires généraux non symétriques. Il minimise la norme résiduelle sur un espace de Krylov à chaque itération, ce qui le rend robuste mais potentiellement gourmand en mémoire pour les grands systèmes [7].

La méthode du résidu minimal (MINRES) : Cette méthode est utilisée pour résoudre les systèmes linéaires indéfinis symétriques. Il est similaire au GMRES mais tire parti de la symétrie de la matrice pour réduire les coûts de calcul [8].

Parmi les autres approches intéressantes, on peut citer la méthode d'orthogonalisation complète (FOM), qui est étroitement liée à la méthode d'Arnoldi pour les problèmes de valeurs propres, la méthode du gradient bi-conjugué ( BiCG ) et la méthode de réduction de la dimension induite (IDR).

1.3 Pourquoi la méthode du sous-espace de Krylov fonctionne-t-elle?

Nous expliquerons ici que la méthode du sous-espace de Krylov devrait être un moyen efficace d'approximer les valeurs propres des matrices par un raffinement itératif des approximations des vecteurs propres, à travers la lentille de la descente la plus abrupte. Nous démontrerons qu'étant donné une supposition initiale d'un état fondamental, l'espace des corrections successives de cette supposition initiale qui produit la convergence la plus rapide est un sous-espace de Krylov. Nous n'apportons pas de preuve rigoureuse du comportement de convergence.

Supposons que notre matrice d'intérêt $A$ soit symétrique et définie positive. C'est ce qui rend notre argument le plus pertinent pour la méthode CG ci-dessus. Nous ne faisons ici aucune hypothèse sur l'éparpillement; nous ne prétendons pas non plus que $A$ doit être un hermitien (ce qu'il doit être s'il s'agit d'un hamiltonien).

Nous souhaitons généralement résoudre un problème de la forme suivante

On peut imaginer que $|b\rangle=c|x\rangle$ où $c$ est une constante, comme dans un problème de valeurs propres. Mais l'énoncé de notre problème reste plus général pour l'instant.

Nous partons d'un vecteur $|x_0\rangle$ qui constitue une solution approximative. Bien qu'il existe des similitudes entre cette hypothèse $|x_0\rangle$ et $|v_0\rangle$ dans la section 1.1, nous ne nous en servons pas ici. Nous pensons que l' $|x_0\rangle$ contient une erreur, que nous appelons « $|e_0\rangle$ » :

Nous définissons également l' $R_0$ résiduelle :

Nous utilisons ici la majuscule $R$ pour distinguer le résidu de la dimension de notre sous-espace de Krylov $r$.

Nous voulons maintenant procéder à une correction de la forme suivante

ce qui, nous l'espérons, améliorera notre approximation. Ici, $|p_0\rangle$ est un vecteur qui reste à déterminer. Soit $|e_1\rangle$ l'erreur après la correction. Ensuite

Nous nous intéressons à la manière dont notre erreur se comporte lorsqu'elle est transformée par notre matrice. Calculons donc la norme $A$ de l'erreur. C'est-à-dire

où nous avons utilisé la symétrie de $A$ ainsi que le fait que $A |e_0\rangle = |R_0\rangle$. Ici, $d$ est une constante indépendante de $|p_0\rangle$. Comme mentionné dans la section 1.2, la norme $A$ de l'erreur n'est pas la seule grandeur que l'on pourrait choisir de minimiser, mais c'est un bon choix. Nous souhaitons observer comment cette grandeur varie en fonction du choix des vecteurs de correction $|p_0\rangle$. Nous définissons donc la fonction $f$ en posant

$f$ est simplement l'erreur $|e_1\rangle$ en fonction de la correction $|p_0\rangle$ mesurée dans la norme $A$. Nous voulons donc choisir $|p_0\rangle$ de telle sorte que $f(|p_0\rangle)$ soit aussi petit que possible. Pour ce faire, nous calculons le gradient de $f$. En utilisant la symétrie de $A$, nous avons

Le gradient indique la direction de la pente la plus raide, ce qui signifie que sa direction opposée correspond à la direction dans laquelle la fonction diminue le plus : la direction de la pente la plus raide en sens inverse. À partir de notre hypothèse initiale $|x_0\rangle$, où $|p_0\rangle=0$, on a que $\nabla f(0) = -2|R_0\rangle$. Ainsi, la fonction $f$ diminue le plus fortement dans la direction du résidu $|R_0\rangle$. Notre choix initial serait donc le plus favorisé par l'ajout du vecteur $|p_0\rangle=\alpha_0 |R_0\rangle$ pour un certain scalaire $\alpha_0$.

Dans l'étape suivante, nous choisissons à nouveau un vecteur $|p_1\rangle$ et ajoutons sa valeur à l'approximation actuelle. En utilisant le même argument que précédemment, nous choisissons $|p_1\rangle = \alpha_1 |R_1\rangle$ pour un certain scalaire $\alpha_1$. Nous continuons ainsi, de sorte que l'itération $k^\text{th}$ de notre vecteur est

De manière équivalente, nous voulons construire l'espace dans lequel nous choisissons nos estimations améliorées en ajoutant $|R_0\rangle$, $|R_1\rangle$, et ainsi de suite, dans l'ordre. Le vecteur estimé $k^\text{th}$ se situe dans

Maintenant, en utilisant la relation que

nous constatons que

Autrement dit, l'espace que nous construisons et qui se rapproche le plus efficacement de la solution correcte $|x\rangle$ est précisément l'espace obtenu par l'application successive de la matrice $A$ sur $|R_0\rangle$. Un sous-espace de Krylov est l'espace engendré par les vecteurs représentant les directions successives de la descente la plus raide.

Enfin, nous rappelons que nous n'avons fait aucune déclaration numérique sur la mise à l'échelle de cette approche et que nous n'avons pas non plus discuté de l'avantage comparatif pour les matrices éparses. Il s'agit uniquement de motiver l'utilisation des méthodes du sous-espace de Krylov et de leur donner un sens intuitif. Nous allons maintenant explorer le comportement de ces méthodes numériquement.

Vérifiez votre compréhension

Dans le processus ci-dessus, nous avons proposé de minimiser la norme $A$ de l'erreur. Quelles autres quantités pourrait-on envisager de minimiser dans la recherche de l'état fondamental et de sa valeur propre?

• Réponse

  On pourrait imaginer d'utiliser le vecteur résiduel au lieu de la norme $A$ de l'erreur. Dans certains cas, il peut être utile de considérer le vecteur d'erreur lui-même.

2. Méthodes de Krylov dans le calcul classique

Dans cette section, nous mettons en œuvre les itérations d'Arnoldi de manière informatique afin de pouvoir exploiter un sous-espace de Krylov pour résoudre les problèmes de valeurs propres. Nous l'appliquerons d'abord à un exemple à petite échelle, puis nous examinerons comment le temps de calcul s'adapte à l'augmentation de la taille de la matrice concernée. Une idée clé ici est que la génération des vecteurs couvrant l'espace de Krylov contribuera largement au temps de calcul total requis. La mémoire nécessaire varie selon les méthodes de Krylov. Mais les contraintes de mémoire peuvent limiter l'utilisation des méthodes traditionnelles de Krylov.

2.1 Exemple simple à petite échelle

Lors de la création d'un sous-espace de Krylov, nous devons orthonormer les vecteurs de notre sous-espace. Définissons une fonction qui prend un vecteur établi dans notre sous-espace vknown (non supposé être normalisé) et un vecteur candidat à ajouter à notre sous-espace vnext et qui rend vnext orthogonal à vknown et normalisé. Définissons également une fonction qui suit ce processus pour tous les vecteurs établis dans notre sous-espace de Krylov afin de garantir un ensemble totalement orthonormé.

# vknown is some established vector in our subspace. vnext is one we wish to add,
# which must be orthogonal to vknown.


def orthog_pair(vknown, vnext):
    vknown = vknown / np.sqrt(vknown.T @ vknown)
    diffvec = vknown.T @ vnext * vknown
    vnext = vnext - diffvec
    return vnext


# v is the candidate vector to be added to our subspace. s is the existing subspace.


def orthoset(v, s):
    v = v / np.sqrt(v.T @ v)
    temp = v
    for i in range(len(s)):
        temp = orthog_pair(s[i], temp)
    v = temp / np.sqrt(temp.T @ temp)
    return v

Définissons maintenant une fonction qui construit un sous-espace de Krylov de plus en plus grand, jusqu'à ce que l'espace des vecteurs de Krylov couvre l'espace complet de la matrice originale. Cela nous permettra de voir dans quelle mesure les valeurs propres obtenues à l'aide de notre méthode de sous-espace de Krylov correspondent aux valeurs exactes, en fonction de la dimension du sous-espace de Krylov. Il est important de noter que notre fonction krylov_full_build renvoie les vecteurs de Krylov, les hamiltoniens projetés, les valeurs propres et le temps nécessaire.

# Necessary imports and definitions to track time in microseconds
import time


def time_mus():
    return int(time.time() * 1000000)


# This function constructs a Krylov subspace that spans the whole space of the original matrix.
#     Input:
#       v0          : initial vector
#       matrix      : original matrix to be diagonalized
#     Output:
#       ks          : Krylov vectors
#       Hs          : projected Hamiltonians
#       eigs        : eigenvalues
#       k_tot_times : time required for the operation


def krylov_full_build(v0, matrix):
    t0 = time_mus()
    b = v0 / np.sqrt(v0 @ v0.T)
    A = matrix
    ks = []
    ks.append(b)
    Hs = []
    eigs = []
    Hs.append(b.T @ A @ b)
    eigs.append(np.array([b.T @ A @ b]))
    k_tot_times = []

    for j in range(len(A) - 1):
        vec = A @ ks[j].T
        ortho = orthoset(vec, ks)
        ks.append(ortho)
        ksarray = np.array(ks)
        Hs.append(ksarray @ A @ ksarray.T)
        eigs.append(np.linalg.eig(Hs[j + 1]).eigenvalues)
        k_tot_times.append(time_mus() - t0)

    # Return the Krylov vectors, the projected Hamiltonians, the eigenvalues,
    # and the total time required.
    return (ks, Hs, eigs, k_tot_times)

Nous allons tester cette méthode sur une matrice qui est encore assez petite, mais plus grande que ce que nous pourrions vouloir faire à la main.

# Define our small test matrix
test_matrix = np.array(
    [
        [4, -1, 0, 1, 0],
        [-1, 4, -1, 2, 1],
        [0, -1, 4, 3, 3],
        [1, 2, 3, 4, 0],
        [0, 1, 3, 0, 4],
    ]
)

# Give the test matrix and an initial guess as arguments in the function defined above.
# Calculate outputs.
test_ks, test_Hs, test_eigs, text_k_tot_times = krylov_full_build(
    np.array([0.5, 0.5, 0, 0.5, 0.5]), test_matrix
)

Nous pouvons vérifier nos fonctions en nous assurant qu'à la dernière étape (lorsque l'espace de Krylov est l'espace vectoriel complet de la matrice originale), les valeurs propres de la méthode de Krylov correspondent exactement à celles de la diagonalisation numérique exacte :

print(np.linalg.eig(test_matrix).eigenvalues)
print(test_eigs[len(test_matrix) - 1])

Output:

[-1.36956923  8.43756009  2.9040308   5.34436028  4.68361806]
[-1.36956923  8.43756009  2.9040308   4.68361806  5.34436028]

C'est un succès. Bien entendu, ce qui compte vraiment, c'est la qualité de notre approximation en fonction de la dimension de notre sous-espace de Krylov. Étant donné que nous sommes souvent préoccupés par la recherche d'états fondamentaux et d'autres valeurs propres minimales (et pour d'autres raisons plus algébriques expliquées ci-dessous), examinons notre estimation de la valeur propre la plus faible en fonction de la dimension du sous-espace de Krylov. C'est-à-dire

def errors(matrix, krylov_eigs):
    targ_min = min(np.linalg.eig(matrix).eigenvalues)
    err = []
    for i in range(len(matrix)):
        err.append(min(krylov_eigs[i]) - targ_min)
    return err

import matplotlib.pyplot as plt

krylov_error = errors(test_matrix, test_eigs)

plt.plot(krylov_error)
plt.axhline(y=0, color="red", linestyle="--")  # Add dashed red line at y=0
plt.xlabel("Order of Krylov subspace")  # Add x-axis label
plt.ylabel("Error in minimum eigenvalue")  # Add y-axis label
plt.show()

Output:

On constate que la valeur propre minimale est atteinte avec une assez grande précision dès que le sous-espace de Krylov a atteint la taille $\mathcal{K}^2$, et qu’elle est parfaite selon $\mathcal{K}^3$.

2.2 Échelle de temps avec dimension matricielle

Convainquons-nous que la méthode de Krylov peut être avantageuse par rapport aux résolveurs numériques exacts de la manière suivante :

• Construire des matrices aléatoires (pas très clairsemées, pas l'application idéale pour la KQD)
• Déterminer les valeurs propres à l'aide de deux méthodes : directement en utilisant NumPy et en utilisant un sous-espace de Krylov.
• Nous choisissons un seuil de précision pour nos valeurs propres, avant d'accepter les estimations de Krylov.
• Comparez le temps nécessaire pour résoudre le problème de ces deux manières.

Mises en garde : Comme nous le verrons en détail ci-dessous, la diagonalisation quantique de Krylov s'applique mieux aux opérateurs dont les représentations matricielles sont peu nombreuses et/ou peuvent être écrites à l'aide d'un petit nombre de groupes d'opérateurs de Pauli commutatifs. Les matrices aléatoires que nous utilisons ici ne correspondent pas à cette description. Celles-ci ne sont utiles que pour sonder l'échelle à laquelle les méthodes de Krylov classiques pourraient être utiles. Deuxièmement, en utilisant la méthode de Krylov, nous calculerons les valeurs propres en utilisant de nombreux sous-espaces de Krylov de tailles différentes. Nous indiquerons le temps nécessaire pour obtenir le sous-espace de Krylov de dimension minimale qui permet d'obtenir la précision requise pour la valeur propre de l'état fondamental. Là encore, il s'agit d'une situation un peu différente de la résolution d'un problème insoluble pour les résolveurs exacts, puisque nous utilisons la solution exacte pour évaluer la dimension nécessaire.

Nous commençons par générer notre ensemble de matrices aléatoires.

import numpy as np

# Set the random seed
np.random.seed(42)

# how many random matrices will we make
num_matrix = 200

matrices = []
for m in range(1, num_matrix):
    matrices.append(np.random.rand(m, m))

Nous allons maintenant diagonaliser chaque matrice directement, en utilisant numpy. Nous calculons le temps nécessaire à la diagonalisation pour une comparaison ultérieure.

matrix_numpy_times = []
matrix_numpy_eigs = []
for mm in range(num_matrix - 1):
    t0 = time_mus()
    matrix_numpy_eigs.append(min(np.linalg.eig(matrices[mm]).eigenvalues))
    matrix_numpy_times.append(time_mus() - t0)

plt.plot(matrix_numpy_times)
plt.xlabel("Dimension of matrix")  # Add x-axis label
plt.ylabel("Time to diagonalize (microsec)")  # Add y-axis label
plt.show()

Output:

Notez que dans l'image ci-dessus, le temps anormalement élevé autour d'une dimension de 125 peut être dû à la nature aléatoire des matrices ou à l'implémentation sur le processeur classique utilisé, mais il n'est pas reproductible. Une nouvelle exécution du code produira un profil différent avec des pics anormaux différents.

Pour chaque matrice, nous allons maintenant construire un sous-espace de Krylov et calculer les valeurs propres par étapes. À chaque étape, nous vérifions si la valeur propre la plus basse a été obtenue avec l'erreur absolue spécifiée. Le sous-espace qui nous donne en premier des valeurs propres dans les limites de l'erreur spécifiée est le sous-espace pour lequel nous enregistrerons les temps de calcul. L'exécution de cette cellule peut prendre plusieurs minutes, en fonction de la vitesse du processeur. N'hésitez pas à sauter l'évaluation ou à réduire la dimension maximale des matrices diagonalisées. Il suffit de regarder les résultats précalculés.

# Choose the absolute error you can tolerate, and make a list for tracking the Krylov subspace size
# at which that error is achieved.
abserr = 0.05
accept_subspace_size = []

# Lists to store total time spent on the Krylov method, and the subset of that time spent on
# diagonalizing the projected matrix.
matrix_krylov_tot_times = []
matrix_krylov_dim = []

# Step through all our random matrices
for mm in range(0, num_matrix - 1):
    test_ks, test_Hs, test_eigs, test_k_tot_times = krylov_full_build(
        np.ones(len(matrices[mm])), matrices[mm]
    )
    # We have not yet found a Krylov subspace that produces our minimum eigenvalue to
    # within the required error.
    found = 0
    for j in range(0, len(matrices[mm]) - 1):
        # If we still haven't found the desired subspace...
        if found == 0:
            # ...but if this one satisfies the requirement, then record everything
            if (
                abs((min(test_eigs[j]) - matrix_numpy_eigs[mm]) / matrix_numpy_eigs[mm])
                < abserr
            ):
                accept_subspace_size.append(j)
                matrix_krylov_tot_times.append(test_k_tot_times[j])
                matrix_krylov_dim.append(mm)
                found = 1

Traçons les temps que nous avons obtenus pour ces deux méthodes à titre de comparaison :

plt.plot(matrix_numpy_times, color="blue")
plt.plot(matrix_krylov_dim, matrix_krylov_tot_times, color="green")
plt.xlabel("Dimension of matrix")  # Add x-axis label
plt.ylabel("Time to diagonalize (microsec)")  # Add y-axis label
plt.show()

Output:

Il s'agit des temps réels nécessaires, mais pour les besoins de la discussion, nous allons lisser ces courbes en faisant la moyenne sur quelques points adjacents / dimensions de la matrice. C'est ce qui est fait ci-dessous :

smooth_numpy_times = []
smooth_krylov_times = []

# Choose the number of adjacent points over which to average forward;
# the same will be used backward.
smooth_steps = 10

# We will do this smoothing for all points/matrix dimensions
for i in range(len(matrix_krylov_tot_times)):
    # Ensure we don't exceed the boundaries of our lists
    start = max(0, i - smooth_steps)
    end = min(len(matrix_krylov_tot_times) - 1, i + smooth_steps)

    # Dummy variables for accumulating an average over adjacent points. This is done for both Krylov
    # and the NumPy calculations.
    smooth_count = 0
    smooth_numpy_sum = 0
    smooth_krylov_sum = 0

    for j in range(start, end):
        smooth_numpy_sum = smooth_numpy_sum + matrix_numpy_times[j]
        smooth_krylov_sum = smooth_krylov_sum + matrix_krylov_tot_times[j]
        smooth_count = smooth_count + 1

    # Appending the averaged adjacent values to our new smooth lists
    smooth_numpy_times.append(smooth_numpy_sum / smooth_count)
    smooth_krylov_times.append(smooth_krylov_sum / smooth_count)

plt.plot(smooth_numpy_times, color="blue")
plt.plot(smooth_krylov_times, color="green")
plt.xlabel("Dimension of matrix")  # Add x-axis label
plt.ylabel("Time to diagonalize (smoothed, microsec)")  # Add y-axis label
plt.show()

Output:

Notez que le temps nécessaire à la construction d'un sous-espace de Krylov dépasse initialement le temps nécessaire à la diagonalisation complète de numpy. Mais lorsque la taille de la matrice augmente, la méthode de Krylov devient avantageuse. Cela est vrai même si nous abaissons notre erreur acceptable, mais l'avantage apparaît lorsque la taille de la matrice est plus grande. Cela vaut la peine d'être décortiqué.

La complexité temporelle de la diagonalisation numérique est de $O(n^3)$ (avec quelques variations d'un algorithme à l'autre). La complexité temporelle de la génération d'une base orthonormée de vecteurs $n$ est également de $O(n^3)$. L'avantage de la méthode de Krylov ne réside donc pas dans l'utilisation d'une base orthonormée $\it{some}$, mais dans l'utilisation d'une base orthonormée particulière qui sélectionne efficacement les valeurs propres d'intérêt. Nous l'avons déjà vu dans l'esquisse d'une preuve dans la première section de cette leçon, et cela est essentiel pour les garanties de convergence dans les méthodes de Krylov.

Passons en revue les progrès réalisés jusqu'à présent :

• Pour les très grandes matrices, la méthode du sous-espace de Krylov peut donner des valeurs propres approximatives dans les tolérances requises plus rapidement que les algorithmes de diagonalisation traditionnels.
• Pour des matrices aussi grandes, la génération d'un sous-espace de Krylov est la partie la plus longue de la méthode du sous-espace de Krylov.
• Il serait donc très utile de disposer d'un moyen efficace de générer un sous-espace de Krylov. C'est enfin là que l'ordinateur quantique entre en scène.

Vérifiez votre compréhension

Reportez-vous au graphique lissé ci-dessus représentant les temps de diagonalisation en fonction de la dimension de la matrice.

(a) D'après ce graphique, à partir de quelle dimension approximative de la matrice la méthode de Krylov est-elle devenue plus rapide?

(b) Quels aspects du calcul pourraient modifier la dimension à partir de laquelle la méthode de Krylov devient plus rapide?

• Réponse

  (a) Les résultats peuvent varier si vous relancez le calcul, mais la méthode de Krylov devient plus rapide à partir d'une dimension d'environ 80 à 85.

  (b) Il existe de nombreuses réponses possibles. Parmi les facteurs importants, on peut citer la précision requise et la稀疏性 des matrices à diagonaliser.

3. Krylov via l'évolution temporelle

Tout ce que nous avons décrit jusqu'à présent peut être réalisé de manière classique. Comment et quand utiliserions-nous un ordinateur quantique? Pour les très grandes matrices, la méthode de Krylov peut nécessiter de longs temps de calcul et de grandes quantités de mémoire. Le temps nécessaire à l'opération matricielle de $H$ sur $|v\rangle$ est équivalent à celui de $O(N^2)$ dans le pire des cas. Même la multiplication de matrices éparses sur un vecteur (le cas typique des solveurs classiques de type Krylov) a une échelle de complexité temporelle de l'ordre de $O(N)$. Ceci est fait pour chaque vecteur que nous voulons dans notre sous-espace. La dimension du sous-espace $r$ n'est généralement pas une fraction significative de $N$, et s'échelonne souvent comme $\log(N)$. Par conséquent, la génération de tous les vecteurs est comparable à $O(N^2 \log(N))$ dans le pire des cas. Bien qu'il y ait d'autres étapes, comme l'orthogonalisation, c'est l'échelle dominante qu'il faut garder à l'esprit.

L'informatique quantique nous permet de modifier les attributs du problème qui déterminent l'échelle du temps et des ressources nécessaires. Au lieu de dépendre de la taille de la matrice $N$, nous verrons des éléments tels que le nombre de tirs et le nombre de termes de Pauli non commutatifs qui composent l'hamiltonien. Voyons comment cela fonctionne.

3.1 Évolution dans le temps

Rappelons que l'opérateur qui fait évoluer dans le temps un état quantique est $e^{-iHt/\hbar}$ (et il est très courant, en particulier en informatique quantique, de supprimer le $\hbar$ de la notation). Une façon de comprendre et même de réaliser une telle fonction exponentielle d'un opérateur est d'examiner son expansion en série de Taylor. Notez que cette opération agissant sur un vecteur initial $|v\rangle$ produit une somme de termes avec des puissances croissantes de $H$ appliquées à l'état initial. Il semble que nous puissions créer notre sous-espace de Krylov en faisant évoluer dans le temps l'état de notre supposition initiale!

La mise en garde concerne la réalisation de l'évolution temporelle sur un véritable ordinateur quantique. De nombreux termes de l'hamiltonien ne commuteront pas entre eux. Ainsi, si certains opérateurs exponentiels simples comme $e^{-iZ}$ correspondent à des circuits simples, ce n'est pas le cas des hamiltoniens généraux. Et comme ils contiennent des termes non commutatifs, nous ne pouvons pas simplement décomposer l'exponentielle en un produit de termes simples, comme nous pouvons le faire avec les nombres.

Ce n'est donc pas trivial, mais c'est un processus bien étudié dans l'informatique quantique. Nous réalisons l'évolution temporelle sur les ordinateurs quantiques à l'aide d'un processus appelé trotterisation, qui est en soi un sujet très riche [10]. Mais à un niveau très élevé, en divisant l'évolution temporelle en très petites étapes, par exemple $m$ étapes de taille $dt$, nous limitons les effets de la non-commutativité des termes.

où $dt = t/m$.

Appelons "sous-espace de Krylov d'ordre r" un sous-espace de Krylov que nous avons généré dans le contexte classique en utilisant directement des puissances de H.

Nous générons maintenant un espace similaire en utilisant l'opérateur d'évolution temporelle unitaire $U \equiv e^{-iHt}$; nous l'appellerons "espace de Krylov unitaire" $\mathcal{K}_U^r$. Le sous-espace de Krylov puissance $\mathcal{K}_P^r$ 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. On peut montrer que l'utilisation du sous-espace de Krylov unitaire donne des garanties de convergence similaires à celles du sous-espace de Krylov puissance, à savoir que l'erreur de l'état fondamental converge efficacement tant que l'état initial $|v\rangle$ a un chevauchement avec l'état fondamental réel qui ne s'évanouit pas exponentiellement, et tant qu'il y a un écart suffisant entre les valeurs propres. Voir la référence [1] pour une discussion plus précise sur la convergence.

Ici, les puissances de $U$ deviennent des pas de temps différents (la puissance $k^\text{th}$ de $U$ avance d'un temps $k \times dt$ ). Nous pouvons appeler $|\psi_k\rangle$ l'élément du sous-espace qui évolue dans le temps pour la durée totale $k dt$.

Nous pouvons projeter notre hamiltonien H sur le sous-espace de Krylov unitaire, $\mathcal{K}_U^r$. En d'autres termes, nous calculons chaque élément de la matrice de $H$ dans la base $\mathcal{K}_U^r$. Nous appellerons cette matrice projetée $\tilde{H}$.

3.2 Comment implémenter sur un ordinateur quantique

Les éléments de la matrice de $\tilde{H}$ sont donnés par les valeurs d'espérance $\langle \psi_m |H| \psi_n\rangle$, qui peuvent être estimées à l'aide de l'ordinateur quantique. N'oubliez pas que $H$ peut être écrit comme une somme d'opérateurs de Pauli sur différents qubits, et que tous les opérateurs de Pauli ne peuvent pas être mesurés simultanément. Nous pouvons trier les termes de Pauli en groupes de termes commutatifs et les mesurer tous en même temps. Mais nous pourrions avoir besoin de plusieurs groupes de ce type pour couvrir tous les termes. Le nombre de groupes distincts de navetteurs dans lesquels les termes peuvent être répartis ( $N_\text{GCP}$ ) devient donc important.

Ici, $P_\alpha$ désigne une chaîne de Pauli de la forme $P_\alpha \sim IZIXII...YZXIX$ ou un ensemble de telles chaînes de Pauli qui commutent entre elles. Étant donné que l'on peut écrire $H$ sous la forme d'une somme d'opérateurs mesurables, les expressions suivantes pour les éléments de matrice de $\tilde{H}$ peuvent être mises en œuvre à l'aide de l'estimateur primitif IBM Quantum.

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 matrice peut être effectué à l'aide de n'importe quel algorithme permettant d'obtenir un chevauchement entre les états quantiques. Dans cette leçon, nous allons nous intéresser au test de Hadamard. Étant donné que l’ $\mathcal{K}_U$ e a pour dimension $r$, l’hamiltonien projeté dans le sous-espace aura pour dimension $r \times r$. Si $r$ est suffisamment petit (en général, $r<<100$ suffit pour garantir la convergence des estimations des valeurs propres), on peut alors facilement diagonaliser l’hamiltonien projeté $\tilde{H}$, de manière classique. Cependant, nous ne pouvons pas diagonaliser directement $\tilde{H}$ en raison de la non-orthogonalité des vecteurs de l'espace de Krylov. Il va falloir mesurer leurs chevauchements et construire une matrice $\tilde{S}$

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)

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

Problème général des valeurs propres

Pourquoi ne pouvons-nous pas simplement diagonaliser $\tilde{H}$? Étant donné que $\tilde{S}$ contient des informations sur la géométrie de la base de Krylov (qui n'est pas orthogonale sauf dans des cas très particuliers), $\tilde{H}$ ne décrit pas à lui seul une projection de l'hamiltonien complet, de sorte que ses valeurs propres n'ont aucune relation particulière avec celles de l'hamiltonien complet - elles pourraient être n'importe quelles valeurs aléatoires. La résolution du problème des valeurs propres généralisées est nécessaire pour obtenir les valeurs propres et les vecteurs propres approximatifs correspondant à la projection de l'hamiltonien complet dans l'espace de Krylov...

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$. La mise en œuvre de cette opération sur un ordinateur quantique est présentée plus en détail ci-dessous.

4. Diagonalisation quantique de Krylov sur un ordinateur quantique

Nous allons maintenant mettre en œuvre la diagonalisation quantique de Krylov sur un véritable ordinateur quantique. Commençons par importer quelques paquets utiles.

import numpy as np
import scipy as sp
import matplotlib.pylab as plt
from typing import Union, List
import warnings

from qiskit.quantum_info import SparsePauliOp, Pauli
from qiskit.circuit import Parameter
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter

# from qiskit.providers.fake_provider import Fake20QV1
from qiskit_ibm_runtime import QiskitRuntimeService, EstimatorV2 as Estimator, Batch

import itertools as it

warnings.filterwarnings("ignore")

Nous définissons la fonction ci-dessous pour résoudre le problème des valeurs propres généralisées que nous venons d'expliquer.

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 our Krylov subspace
        s (numpy.ndarray):
            The matrix of overlaps between vectors of our 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]

Au moins dans le cadre d'une analyse comparative initiale, il est utile de connaître une solution classique exacte pour vérifier le comportement de la convergence. La fonction ci-dessous calcule l'énergie de l'état fondamental d'un hamiltonien, en utilisant le hamiltonien et le nombre de qubits comme arguments.

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)

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

    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

4.1 Étape 1 : Mettre en correspondance le problème avec les circuits quantiques et les opérateurs

Nous allons maintenant définir un hamiltonien. Cette fonction se distingue de la fonction ci-dessus en ce sens qu'elle prend un hamiltonien comme argument et ne renvoie que l'état fondamental, et ce de manière classique. L'hamiltonien que nous définissons ici détermine les niveaux d'énergie de tous les états propres énergétiques, et cet hamiltonien peut être construit à l'aide d'opérateurs de Pauli et mis en œuvre sur un ordinateur quantique.

Nous choisissons un hamiltonien correspondant à une chaîne de spins qui peut avoir n'importe quelle orientation dans l'espace, appelée "chaîne de Heisenberg". Nous supposons que le spin $i^\text{th}$ peut être influencé par ses voisins les plus proches (les spins $(i-1)^\text{th}$ et $(i+1)^\text{th}$ ) mais pas par des voisins plus éloignés. Nous tenons également compte de la possibilité que l'interaction entre les spins soit différente lorsque les spins pointent le long d'axes différents. Cette asymétrie est parfois due, par exemple, à la structure du réseau cristallin dans lequel les spins sont intégrés.

# Define problem Hamiltonian.
n_qubits = 10
# coupling strength for XX, YY, and ZZ interactions
JX = 1
JY = 3
JZ = 2

# 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, JZ)
    if term.count("Z") == 2
    else (term, JY)
    if term.count("Y") == 2
    else (term, JX)
    for term in H_int
]

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