Krylov-Quanten-Diagonalisierung

In dieser Lektion über die Krylov-Quantendiagonalisierung (KQD) werden wir folgende Fragen beantworten:

• Was ist die Krylov-Methode im Allgemeinen?
• Warum funktioniert die Krylov-Methode und unter welchen Bedingungen?
• Welche Rolle spielt das Quantencomputing?

Der Quantenteil der Berechnungen stützt sich weitgehend auf die Arbeiten in [1].

Das folgende Video gibt einen Überblick über die Krylov-Methoden in der klassischen Datenverarbeitung, begründet ihre Verwendung und erklärt, wie die Quanteninformatik in diesem Arbeitsbereich eine Rolle spielen kann. Der folgende Text bietet mehr Details und implementiert eine Krylov-Methode sowohl klassisch als auch unter Verwendung eines Quantencomputers.

1. Einführung in die Krylov-Methoden

Eine Krylov-Unterraum-Methode kann sich auf eine von mehreren Methoden beziehen, die auf dem so genannten Krylov-Unterraum aufbauen. Eine vollständige Übersicht über diese Themen würde den Rahmen dieser Lektion sprengen, aber die Referenzen [2-4] können wesentlich mehr Hintergrundinformationen liefern. Hier werden wir uns darauf konzentrieren, was ein Krylov-Unterraum ist, wie und warum er bei der Lösung von Eigenwertproblemen nützlich ist, und schließlich, wie er auf einem Quantencomputer implementiert werden kann.

Definition: Bei einer symmetrischen, positiv semidefiniten $N\times N$ Matrix $A$ ist der Krylov-Raum $\mathcal{K}^r$ der Ordnung $r$ der Raum, der von Vektoren aufgespannt wird, die durch Multiplikation höherer Potenzen einer Matrix $A$ bis zu $r-1\leq N$ mit einem Referenzvektor $\vert v \rangle$ erhalten werden.

Obwohl die obigen Vektoren einen so genannten Krylov-Unterraum aufspannen, gibt es keinen Grund zu der Annahme, dass sie orthogonal sein werden. Häufig wird ein iteratives Orthonormalisierungsverfahren ähnlich der Gram-Schmidt-Orthogonalisierung verwendet. Hier ist der Prozess etwas anders, da jeder neue Vektor bei seiner Erzeugung orthogonal zu den anderen gemacht wird. In diesem Zusammenhang wird dies als Arnoldi-Iteration bezeichnet. Ausgehend von dem Anfangsvektor $|v\rangle$ erzeugt man den nächsten Vektor $A|v\rangle$ und stellt dann sicher, dass dieser zweite Vektor orthogonal zum ersten ist, indem man seine Projektion auf $|v\rangle$ abzieht. Das bedeutet

Es ist nun leicht zu erkennen, dass $|v_0\rangle \perp |v_1\rangle$ gilt, da

Das Gleiche gilt für den nächsten Vektor, wobei sichergestellt wird, dass er orthogonal zu den beiden vorherigen Vektoren ist:

Wenn wir diesen Vorgang für alle $r$ Vektoren wiederholen, haben wir eine vollständige orthonormale Basis für einen Krylov-Raum. Man beachte, dass der Orthogonalisierungsprozess hier einmal $r>m$ zu Null führt, da $m$ orthogonale Vektoren notwendigerweise den gesamten Raum umfassen. Der Prozess wird auch Null ergeben, wenn ein beliebiger Vektor ein Eigenvektor von $A$ ist, da alle nachfolgenden Vektoren Vielfache dieses Vektors sind.

1.1 Ein einfaches Beispiel: Krylov von Hand

Gehen wir schrittweise durch die Generierung eines Krylov-Unterraums auf einer trivial kleinen Matrix, damit wir den Prozess sehen können. Wir beginnen mit einer Ausgangsmatrix $A$, die für uns von Interesse ist:

Für dieses kleine Beispiel können wir die Eigenvektoren und Eigenwerte sogar leicht von Hand bestimmen. Wir zeigen hier die numerische Lösung.

# 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]]

Wir halten sie hier für einen späteren Vergleich fest:

Wir möchten untersuchen, wie dieser Prozess funktioniert (oder scheitert), wenn wir die Dimension unseres Krylov-Unterraums $r$ erhöhen. Zu diesem Zweck werden wir dieses Verfahren anwenden:

• Erzeugen Sie einen Unterraum des vollständigen Vektorraums, beginnend mit einem zufällig gewählten Vektor $|v\rangle$ (nennen Sie ihn $|v_0\rangle$, wenn er bereits normalisiert ist, wie oben).
• Projizieren Sie die vollständige Matrix $A$ auf diesen Unterraum, und finden Sie die Eigenwerte dieser projizierten Matrix $\tilde{A}$.
• Vergrößern Sie den Unterraum, indem Sie weitere Vektoren erzeugen und sicherstellen, dass sie orthonormal sind, indem Sie ein Verfahren ähnlich der Gram-Schmidt-Orthogonalisierung anwenden.
• Projizieren Sie $A$ auf den größeren Unterraum und finden Sie die Eigenwerte der resultierenden Matrix $\tilde{A}$.
• Wiederholen Sie dies so lange, bis die Eigenwerte konvergieren (oder in diesem Spielzeugfall, bis Sie Vektoren erzeugt haben, die den gesamten Vektorraum der ursprünglichen Matrix abdecken $A$ ).

Bei einer normalen Implementierung der Krylov-Methode müsste das Eigenwertproblem für die auf jeden Krylov-Unterraum projizierte Matrix nicht gelöst werden, während sie aufgebaut wird. Man könnte den Unterraum der gewünschten Dimension konstruieren, die Matrix auf diesen Unterraum projizieren und die projizierte Matrix diagonalisieren. Die Projektion und Diagonalisierung in jeder Unterraumdimension wird nur zur Überprüfung der Konvergenz durchgeführt.

$r=1$ der Dimension:

Wir wählen einen Zufallsvektor, sagen wir

Wenn er nicht bereits normalisiert ist, normalisieren Sie ihn.

Wir projizieren nun unsere Matrix $A$ auf den Unterraum dieses einen Vektors:

Dies ist unsere Projektion der Matrix auf unseren Krylov-Unterraum, wenn dieser nur einen einzigen Vektor enthält, $|v_0\rangle$. Der Eigenwert dieser Matrix ist trivialerweise 4. Wir können uns dies als unsere Schätzung der Eigenwerte (in diesem Fall nur einer) von $A$ vorstellen. Obwohl es eine schlechte Schätzung ist, ist es die richtige Größenordnung.

$r=2$ der Dimension:

Wir erzeugen nun den nächsten Vektor in unserem Unterraum durch eine Operation mit $A$ auf dem vorherigen Vektor:

Nun subtrahieren wir die Projektion dieses Vektors auf unseren vorherigen Vektor, um die Orthogonalität zu gewährleisten.

Wenn er nicht bereits normalisiert ist, normalisieren Sie ihn. In diesem Fall wurde der Vektor bereits normalisiert, so dass

Wir projizieren nun unsere Matrix A auf den Unterraum dieser beiden Vektoren:

Bleibt noch das Problem, die Eigenwerte dieser Matrix zu bestimmen. Diese Matrix ist jedoch etwas kleiner als die Vollmatrix. Bei Problemen mit sehr großen Matrizen kann die Arbeit mit diesem kleineren Unterraum von großem Vorteil sein.

Dies ist zwar immer noch keine gute Schätzung, aber besser als die Schätzung nullter Ordnung. Wir werden dies in einer weiteren Iteration durchführen, um sicherzustellen, dass der Prozess klar ist. Dies untergräbt jedoch den Sinn der Methode, da wir in der nächsten Iteration eine 3x3 -Matrix diagonalisieren werden, was bedeutet, dass wir weder Zeit noch Rechenleistung gespart haben.

$r=3$ der Dimension:

Wir erzeugen nun den nächsten Vektor in unserem Unterraum durch eine Operation mit A auf dem vorherigen Vektor:

Nun subtrahieren wir die Projektion dieses Vektors auf die beiden vorherigen Vektoren, um die Orthogonalität sicherzustellen.

Wenn er nicht bereits normalisiert ist, normalisieren Sie ihn. In diesem Fall wurde der Vektor bereits normalisiert, so dass

Wir projizieren nun unsere Matrix $A$ auf den Unterraum dieser Vektoren:

Wir bestimmen nun die Eigenwerte:

Diese Eigenwerte sind genau die Eigenwerte der Originalmatrix $A$. Dies muss der Fall sein, da wir unseren Krylov-Unterraum so erweitert haben, dass er den gesamten Vektorraum der Originalmatrix $A$ umfasst.

In diesem Beispiel mag die Krylov-Methode nicht besonders einfacher erscheinen als die direkte Diagonalisierung. Wie wir in späteren Abschnitten sehen werden, ist die Krylov-Methode erst ab einer bestimmten Matrixdimension vorteilhaft; dies soll uns helfen, Eigenwert-/Eigenvektor-Probleme von extrem großen Matrizen zu lösen.

Dies ist das einzige Beispiel, das wir zeigen werden, das "von Hand" bearbeitet wurde, aber Abschnitt 2 unten zeigt rechnerische Beispiele.

Begriffsklärung

Ein häufiges Missverständnis ist, dass es nur einen einzigen Krylov-Unterraum für ein bestimmtes Problem gibt. Da es aber viele Ausgangsvektoren gibt, auf die unsere Matrix angewendet werden könnte, gibt es natürlich auch viele mögliche Krylov-Unterräume. Wir werden den Ausdruck "die Krylov-Unterraum" nur verwenden, um auf einen bestimmten Krylov-Unterraum zu verweisen, der bereits für ein bestimmtes Beispiel definiert wurde. Für allgemeine Problemlösungsansätze sprechen wir von einem " Krylov-Unterraum". Eine letzte Klarstellung ist, dass die Bezeichnung " Krylov-Raum " zulässig ist. Er wird oft als " Krylov-Unterraum " bezeichnet, weil er im Zusammenhang mit der Projektion von Matrizen aus einem Ausgangsraum in einen Unterraum verwendet wird. In diesem Zusammenhang werden wir ihn hier meist als Unterraum bezeichnen.

Überprüfen Sie Ihr Verständnis

Erklären Sie, warum es (a) nicht sinnvoll und (b) nicht möglich ist, die Dimension des Krylov-Unterraums $r$ über die Dimension $N$ der interessierenden Matrix hinaus zu erweitern.

• Antwort

  (a) Da wir die Vektoren bei ihrer Erzeugung orthonormalisieren, bildet eine Menge von $N$ solchen Vektoren eine vollständige Basis, was bedeutet, dass sich jeder Vektor im Raum als Linearkombination dieser Vektoren darstellen lässt.

  (b) Der Orthogonalisierungsprozess besteht darin, die Projektion eines neuen Vektors auf alle vorherigen Vektoren zu subtrahieren. Wenn alle vorherigen Vektoren den gesamten Vektorraum aufspannen, erhalten wir durch Subtraktion der Projektionen auf den vollständigen Unterraum immer einen Nullvektor.

Angenommen, ein Kollege demonstriert die Krylov-Methode anhand einer kleinen Beispielmatrix. Stimmt etwas nicht mit ihrer Wahl der Matrix $A$ und des Startvektors $|\psi\rangle$?

und

• Antwort

  Ihr Kollege hat versehentlich einen Eigenvektor für seinen Ausgangsvektor gewählt. Wenn man mit der Matrix auf den Ausgangsvektor einwirkt, erhält man einfach denselben Vektor zurück, skaliert um den Eigenwert. Dadurch wird kein Unterraum mit zunehmender Dimension erzeugt. Raten Sie Ihrem Kollegen, einen anderen Ausgangsvektor zu wählen und darauf zu achten, dass es sich nicht um einen Eigenvektor handelt.

Wende die Krylov-Methode auf die angegebene Matrix an und wähle dabei einen geeigneten neuen Startvektor aus. Schreiben Sie die Schätzungen des minimalen Eigenwerts für die 0. und 1. Ordnung Ihres Krylov-Unterraums auf.

• Antwort

  Je nach Wahl des Ausgangsvektors gibt es viele mögliche Antworten. Wir werden wählen:

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

  Um „ $|v_1\rangle$ “ zu erhalten, wenden wir „ $A$ “ einmal auf „ $|v_0\rangle$ “ an und machen dann „ $|v_1\rangle$ “ orthogonal zu „ $|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}$$

  Bei 0-ter Ordnung ist die Projektion auf unseren Krylov-Unterraum

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

  Bei 1. Ordnung ist die Projektion auf diesen Krylov-Unterraum

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

  Dies kann von Hand geschehen, ist aber am einfachsten mit Numpy zu bewerkstelligen:

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

  Die minimale Eigenwertschätzung lautet -0.414.

1.2 Arten von Krylov-Methoden

"Krylov-Unterraum-Methoden" können sich auf verschiedene iterative Techniken beziehen, die zur Lösung großer linearer Systeme und Eigenwertprobleme verwendet werden. Sie alle haben gemeinsam, dass sie eine Näherungslösung aus einem Krylov-Unterraum konstruieren

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

wobei $|v\rangle$ die anfängliche Schätzung ist (siehe Ref [5] ). Sie unterscheiden sich darin, wie sie die beste Annäherung aus diesem Unterraum auswählen und dabei Faktoren wie Konvergenzrate, Speicherverbrauch und Gesamtberechnungskosten ausgleichen. Der Schwerpunkt dieser Lektion liegt auf der Nutzung von Quantencomputern im Zusammenhang mit Krylov-Unterraummethoden; eine ausführliche Diskussion dieser Methoden würde den Rahmen sprengen. Die nachstehenden kurzen Definitionen dienen nur der Veranschaulichung und enthalten einige Verweise für die weitere Untersuchung dieser Methoden.

Die Methode des konjugierten Gradienten (CG) : Diese Methode wird zur Lösung symmetrischer, positiv definitiver linearer Systeme verwendet [6]. Es minimiert die A-Norm des Fehlers bei jeder Iteration, was es besonders effektiv für Systeme macht, die aus diskretisierten elliptischen PDEs entstehen [7]. Wir werden diesen Ansatz im nächsten Abschnitt verwenden, um zu begründen, warum ein Krylov-Unterraum ein effektiver Unterraum ist, in dem man nach besseren Lösungen für lineare Systeme suchen kann.

Die verallgemeinerte Methode der minimalen Residuen (GMRES) : Diese Methode ist für die Lösung allgemeiner, nicht symmetrischer linearer Systeme konzipiert. Es minimiert die Restnorm über einen Krylov-Raum bei jeder Iteration, was es robust, aber potenziell speicherintensiv für große Systeme macht [7].

Die Methode der minimalen Residuen (MINRES) : Diese Methode wird für die Lösung symmetrischer unbestimmter linearer Systeme verwendet. Es ist ähnlich wie GMRES, nutzt aber die Matrixsymmetrie, um die Rechenkosten zu senken [8].

Weitere erwähnenswerte Ansätze sind die Methode der vollständigen Orthogonalisierung (FOM ), die eng mit der Arnoldi-Methode für Eigenwertprobleme verwandt ist, die Methode des bi-konjugierten Gradienten ( BiCG ) und die Methode der induzierten Dimensionsreduktion (IDR).

1.3 Warum die Krylov-Unterraummethode funktioniert

Hier werden wir motivieren, dass die Krylov-Unterraum-Methode ein effizienter Weg sein sollte, Matrixeigenwerte durch iterative Verfeinerung von Eigenvektor-Näherungen durch die Linse des steilsten Abstiegs zu approximieren. Wir werden argumentieren, dass bei einer anfänglichen Schätzung eines Grundzustands der Raum der aufeinanderfolgenden Korrekturen dieser anfänglichen Schätzung, der die schnellste Konvergenz ergibt, ein Krylov-Unterraum ist. Wir verzichten auf einen rigorosen Nachweis des Konvergenzverhaltens.

Nehmen wir an, unsere interessierende Matrix $A$ ist symmetrisch und positiv definit. Dies macht unsere Argumentation für die obige CG-Methode besonders relevant. Wir machen hier keine Annahmen über die Sparsamkeit; wir behaupten auch nicht, dass $A$ hermitesch sein muss (was es sein muss, wenn es ein Hamiltonian ist).

Normalerweise wollen wir ein Problem der folgenden Form lösen

Man könnte sich vorstellen, dass $|b\rangle=c|x\rangle$, wobei $c$ eine Konstante ist, wie bei einem Eigenwertproblem. Aber unsere Problemstellung bleibt vorerst allgemeiner.

Wir beginnen mit einem Vektor $|x_0\rangle$, der eine Näherungslösung darstellt. Obwohl es Parallelen zwischen dieser Vermutung $|x_0\rangle$ und $|v_0\rangle$ im Abschnitt 1.1 gibt, ziehen wir hier keinen Nutzen daraus. Unser Vermutungs $|x_0\rangle$ enthält einen Fehler, den wir als „ $|e_0\rangle$ “ bezeichnen:

Außerdem definieren wir die Rest $R_0$ wie folgt:

Hier verwenden wir die Großbuchstaben $R$, um das Residuum von der Dimension unseres Krylov-Unterraums $r$ zu unterscheiden.

Wir wollen nun einen Korrekturschritt der Form

was hoffentlich unsere Annäherung verbessert. Hier ist $|p_0\rangle$ ein noch zu bestimmender Vektor. $|e_1\rangle$ sei der Fehler nach der Korrektur. Dann

Wir sind daran interessiert, wie sich unser Fehler verhält, wenn er durch unsere Matrix transformiert wird. Berechnen wir also die $A$ -norm des Fehlers. Das heißt

Dabei haben wir die Symmetrie von $A$ sowie die Tatsache, dass $A |e_0\rangle = |R_0\rangle$, genutzt. Dabei ist $d$ eine von $|p_0\rangle$ unabhängige Konstante. Wie in Abschnitt 1.2 erwähnt, ist die $A$ -Norm des Fehlers nicht die einzige Größe, die wir zur Minimierung wählen könnten, aber sie ist eine gute Wahl. Wir möchten untersuchen, wie sich diese Größe in Abhängigkeit von der Wahl unserer Korrekturvektoren $|p_0\rangle$ verändert. Dazu definieren wir die Funktion $f$ durch Setzen von

$f$ ist einfach der Fehler $|e_1\rangle$ in Abhängigkeit von der Korrektur $|p_0\rangle$ gemessen in der $A$ -Norm. Daher wollen wir $|p_0\rangle$ so wählen, dass $f(|p_0\rangle)$ so klein wie möglich ist. Zu diesem Zweck berechnen wir den Gradienten von $f$. Unter Verwendung der Symmetrie von $A$ ergibt sich

Die Steigung zeigt in die Richtung des steilsten Anstiegs, was bedeutet, dass ihre entgegengesetzte Richtung die Richtung angibt, in der die Funktion am stärksten abfällt: die Richtung des steilsten Abfalls. Bei unserer anfänglichen Vermutung $|x_0\rangle$, wobei $|p_0\rangle=0$ gilt, ergibt sich $\nabla f(0) = -2|R_0\rangle$. Somit nimmt die Funktion $f$ in Richtung des Residuals $|R_0\rangle$ am stärksten ab. Daher würde unsere anfängliche Wahl am meisten von der Addition des Vektors $|p_0\rangle=\alpha_0 |R_0\rangle$ mit einem bestimmten Skalar $\alpha_0$ profitieren.

Im nächsten Schritt wählen wir wieder einen Vektor $|p_1\rangle$ und addieren seinen Wert zur aktuellen Approximation. Mit demselben Argument wie zuvor wählen wir $|p_1\rangle = \alpha_1 |R_1\rangle$ für einen Skalar $\alpha_1$. Wir fahren auf diese Weise fort, so dass die $k^\text{th}$ Iteration unseres Vektors ist

Entsprechend wollen wir den Raum, aus dem wir unsere verbesserten Schätzungen auswählen, durch Hinzufügen von $|R_0\rangle$, $|R_1\rangle$, usw. in der Reihenfolge aufbauen. Der geschätzte Vektor $k^\text{th}$ liegt in

Unter Verwendung der Beziehung, dass

wir sehen, dass

Das heißt, der Raum, den wir aufbauen und der die richtige Lösung $|x\rangle$ am effizientesten approximiert, ist genau der Raum, der durch die aufeinanderfolgende Anwendung der Matrix $A$ auf $|R_0\rangle$ entsteht. Ein Krylov-Unterraum ist der Raum, der durch die Vektoren der aufeinanderfolgenden Richtungen des steilsten Abstiegs aufgespannt wird.

Abschließend möchten wir noch einmal darauf hinweisen, dass wir keine numerischen Aussagen über die Skalierung dieses Ansatzes gemacht haben, noch haben wir den vergleichbaren Nutzen für dünn besetzte Matrizen diskutiert. Dies soll nur dazu dienen, die Verwendung von Krylov-Unterraum-Methoden zu motivieren und ihnen einen gewissen intuitiven Sinn zu verleihen. Wir werden nun das Verhalten dieser Methoden numerisch untersuchen.

Überprüfen Sie Ihr Verständnis

In dem obigen Arbeitsablauf haben wir vorgeschlagen, die $A$ -Norm des Fehlers zu minimieren. Welche anderen Größen könnte man bei der Suche nach dem Grundzustand und seinem Eigenwert minimieren?

• Antwort

  Man könnte sich vorstellen, den Residualvektor anstelle der $A$ -Norm des Fehlers zu verwenden. Es kann Fälle geben, in denen die Betrachtung des Fehlervektors selbst sinnvoll ist.

2. Krylov-Methoden in der klassischen Berechnung

In diesem Abschnitt implementieren wir Arnoldi-Iterationen rechnerisch, so dass wir einen Krylov-Unterraum bei der Lösung von Eigenwertproblemen nutzen können. Wir werden dies zunächst auf ein kleines Beispiel anwenden und dann untersuchen, wie die Berechnungszeit mit zunehmender Größe der interessierenden Matrix steigt. Ein Schlüsselgedanke dabei ist, dass die Erzeugung der Vektoren, die den Krylov-Raum aufspannen, einen großen Anteil an der benötigten Gesamtrechenzeit hat. Der benötigte Speicherplatz variiert zwischen den einzelnen Krylov-Methoden. Allerdings können Speicherbeschränkungen die Verwendung traditioneller Krylov-Methoden einschränken.

2.1 Einfaches Beispiel im kleinen Maßstab

Bei der Erstellung eines Krylov-Unterraums müssen wir die Vektoren in unserem Unterraum orthonormalisieren. Definieren wir eine Funktion, die einen etablierten Vektor aus unserem Unterraum vknown (der nicht normalisiert sein soll) und einen Kandidatenvektor nimmt, den wir zu unserem Unterraum vnext hinzufügen und vnext orthogonal zu vknown und normalisiert machen. Definieren wir ferner eine Funktion, die diesen Prozess für alle etablierten Vektoren in unserem Krylov-Unterraum durchläuft, um eine vollständig orthonormale Menge zu gewährleisten.

# 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

Definieren wir nun eine Funktion, die iterativ einen immer größeren Krylov-Unterraum aufbaut, bis der Raum der Krylov-Vektoren den gesamten Raum der ursprünglichen Matrix umfasst. So können wir sehen, wie gut die mit unserer Krylov-Unterraum-Methode erhaltenen Eigenwerte mit den exakten Werten übereinstimmen, und zwar in Abhängigkeit von der Krylov-Unterraum-Dimension. Wichtig ist, dass unsere Funktion krylov_full_build die Krylov-Vektoren, die projizierten Hamiltonianer, die Eigenwerte und die benötigte Zeit liefert.

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

Wir werden dies an einer Matrix testen, die zwar immer noch recht klein ist, aber größer, als wir sie vielleicht von Hand erstellen wollen.

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

Wir können unsere Funktionen überprüfen, indem wir sicherstellen, dass im letzten Schritt (wenn der Krylov-Raum der vollständige Vektorraum der ursprünglichen Matrix ist) die Eigenwerte der Krylov-Methode genau mit denen der exakten numerischen Diagonalisierung übereinstimmen:

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]

Das war erfolgreich. Entscheidend ist natürlich, wie gut unsere Annäherung in Abhängigkeit von der Dimension unseres Krylov-Unterraums ist. Da wir uns häufig mit der Suche nach Grundzuständen und anderen minimalen Eigenwerten beschäftigen (und aus anderen, eher algebraischen Gründen, die weiter unten erläutert werden), betrachten wir unsere Schätzung des niedrigsten Eigenwerts als Funktion der Dimension des Krylov-Unterraums. Das heißt

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:

Wir sehen, dass der minimale Eigenwert ziemlich genau erreicht wird, sobald der Krylov-Unterraum auf $\mathcal{K}^2$ angewachsen ist, und gemäß $\mathcal{K}^3$ perfekt ist.

2.2 Zeitskalierung mit Matrixdimension

Wir wollen uns davon überzeugen, dass die Krylov-Methode gegenüber exakten numerischen Eigenlösungen auf folgende Weise vorteilhaft sein kann:

• Konstruktion von Zufallsmatrizen (nicht spärlich, nicht die ideale Anwendung für KQD)
• Bestimmen Sie Eigenwerte mit zwei Methoden: direkt mit NumPy und mit einem Krylov-Unterraum.
• Wir wählen einen Grenzwert dafür, wie genau unsere Eigenwerte sein müssen, bevor wir die Krylov-Schätzungen akzeptieren.
• Vergleichen Sie die Wandzeit, die für diese beiden Lösungswege benötigt wird.

Vorsichtsmaßnahmen: Wie wir weiter unten ausführlich erörtern werden, lässt sich die Krylov-Quantendiagonalisierung am besten auf Operatoren anwenden, deren Matrixdarstellungen spärlich sind und/oder mit einer kleinen Anzahl von Gruppen von kommutierenden Pauli-Operatoren geschrieben werden können. Die Zufallsmatrizen, die wir hier verwenden, entsprechen nicht dieser Beschreibung. Diese sind nur nützlich, um die Skala auszuloten, bei der klassische Krylov-Methoden nützlich sein könnten. Zweitens werden bei der Anwendung der Krylov-Methode die Eigenwerte unter Verwendung vieler unterschiedlich großer Krylov-Unterräume berechnet. Wir werden die Zeit angeben, die für den minimal dimensionierten Krylov-Unterraum benötigt wird, der die von uns geforderte Genauigkeit für den Grundzustandseigenwert erreicht. Auch dies ist etwas anderes als die Lösung eines Problems, das für exakte Eigenlöser unlösbar ist, da wir die exakte Lösung verwenden, um die benötigte Dimension zu bestimmen.

Wir beginnen mit der Erstellung einer Reihe von Zufallsmatrizen.

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

Wir diagonalisieren nun jede Matrix direkt mit Numpy. Zum späteren Vergleich berechnen wir die für die Diagonalisierung benötigte Zeit.

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:

Man beachte, dass in der obigen Abbildung die anomal hohe Zeit um eine Dimension von 125 herum auf die zufällige Natur der Matrizen oder auf die Implementierung auf dem verwendeten klassischen Prozessor zurückzuführen sein könnte, aber nicht reproduzierbar ist. Wenn Sie den Code erneut ausführen, erhalten Sie ein anderes Profil mit anderen anomalen Spitzenwerten.

Nun werden wir für jede Matrix einen Krylov-Unterraum aufbauen und schrittweise Eigenwerte berechnen. Bei jedem Schritt wird geprüft, ob der niedrigste Eigenwert innerhalb des von uns angegebenen absoluten Fehlers erzielt wurde. Der Unterraum, der als erster Eigenwerte innerhalb des angegebenen Fehlers liefert, ist der Unterraum, für den wir die Berechnungszeiten aufzeichnen werden. Die Ausführung dieser Zelle kann je nach Prozessorgeschwindigkeit mehrere Minuten dauern. Es steht Ihnen frei, die Auswertung zu überspringen oder die maximale Dimension der diagonalisierten Matrizen zu verringern. Es genügt, sich die vorberechneten Ergebnisse anzusehen.

# 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

Wir wollen die Zeiten, die wir für diese beiden Methoden erhalten haben, zum Vergleich aufzeichnen:

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:

Dies sind die tatsächlich benötigten Zeiten, aber für die Zwecke der Diskussion glätten wir diese Kurven, indem wir den Durchschnitt über einige benachbarte Punkte / Matrixdimensionen bilden. Dies geschieht im Folgenden:

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:

Beachten Sie, dass die Zeit, die für die Erstellung eines Krylov-Unterraums benötigt wird, zunächst die Zeit übersteigt, die für die vollständige Diagonalisierung in Numpy benötigt wird. Mit zunehmender Größe der Matrix wird die Krylov-Methode jedoch vorteilhaft. Dies gilt auch dann, wenn wir unseren akzeptablen Fehler verringern, aber der Vorteil setzt bei einer größeren Matrixgröße ein. Es lohnt sich, dies auseinander zu nehmen.

Die Zeitkomplexität der numerischen Diagonalisierung ist $O(n^3)$ (mit einigen Abweichungen zwischen den Algorithmen). Die Zeitkomplexität der Erzeugung einer orthonormalen Basis von $n$ Vektoren ist ebenfalls $O(n^3)$. Der Vorteil der Krylov-Methode liegt also nicht in der Verwendung der orthonormalen Basis $\it{some}$, sondern in der Verwendung einer bestimmten orthonormalen Basis, die effektiv die Eigenwerte von Interesse herausgreift. Wir haben dies bereits in der Skizze eines Beweises im ersten Abschnitt dieser Lektion gesehen, und dies ist entscheidend für die Konvergenzgarantien in Krylov-Verfahren.

Lassen Sie uns einen Blick auf unsere bisherigen Fortschritte werfen:

• Bei sehr großen Matrizen kann die Krylov-Unterraum-Methode schneller als herkömmliche Diagonalisierungsalgorithmen approximative Eigenwerte innerhalb der erforderlichen Toleranzen liefern.
• Bei solch großen Matrizen ist die Erzeugung eines Krylov-Unterraums der zeitaufwändigste Teil der Krylov-Unterraum-Methode.
• Daher wäre ein effizientes Verfahren zur Erzeugung eines Krylov-Unterraums von großem Nutzen. Hier kommt endlich der Quantencomputer ins Spiel.

Überprüfen Sie Ihr Verständnis

Siehe die obenstehende geglättete Darstellung der Diagonalisierungszeiten in Abhängigkeit von der Matrixdimension.

(a) Ab welcher Matrixdimension wurde die Krylov-Methode laut dieser Grafik schneller?

(b) Welche Aspekte der Berechnung könnten den Bereich beeinflussen, ab dem die Krylov-Methode schneller wird?

• Antwort

  (a) Die Ergebnisse können variieren, wenn man die Berechnung erneut durchführt, doch die Krylov-Methode wird ab einer Dimension von etwa 80–85 schneller.

  (b) Es gibt viele mögliche Antworten. Wichtige Faktoren sind die von uns geforderte Genauigkeit und die Sparsität der zu diagonalizierenden Matrizen.

3. Krylov über Zeitentwicklung

Alles, was wir bisher beschrieben haben, kann klassisch durchgeführt werden. Wie und wann würden wir also einen Quantencomputer verwenden? Bei sehr großen Matrizen kann die Krylov-Methode lange Rechenzeiten und große Mengen an Speicherplatz erfordern. Die für die Matrixoperation von $H$ auf $|v\rangle$ benötigte Zeit skaliert im schlimmsten Fall wie $O(N^2)$. Selbst die Multiplikation dünn besetzter Matrizen auf einem Vektor (der typische Fall für klassische Krylov-Löser) hat eine Zeitkomplexität, die wie $O(N)$ skaliert. Dies wird für jeden gewünschten Vektor in unserem Unterraum durchgeführt. Die Unterraumdimension $r$ ist normalerweise kein signifikanter Bruchteil von $N$ und skaliert oft wie $\log(N)$. Die Erzeugung aller Vektoren skaliert also im schlimmsten Fall wie $O(N^2 \log(N))$. Obwohl es noch andere Schritte gibt, wie z. B. die Orthogonalisierung, ist dies die wichtigste Skalierung, die man im Auge behalten sollte.

Das Quantencomputing ermöglicht es uns, die Eigenschaften des Problems zu ändern, die die Skalierung der benötigten Zeit und Ressourcen bestimmen. Anstelle der generellen Abhängigkeit von der Matrixgröße $N$ werden wir Dinge wie die Anzahl der Schüsse und die Anzahl der nicht-kommutierenden Pauli-Terme sehen, aus denen sich der Hamiltonian zusammensetzt. Sehen wir uns an, wie das funktioniert.

3.1 Zeitliche Entwicklung

Der Operator, der einen Quantenzustand zeitlich entwickelt, ist $e^{-iHt/\hbar}$ (und es ist sehr üblich, insbesondere in der Quanteninformatik, $\hbar$ aus der Notation zu streichen). Eine Möglichkeit, eine solche Exponentialfunktion eines Operators zu verstehen und sogar zu realisieren, ist die Betrachtung ihrer Taylorreihenentwicklung. Man beachte, dass diese Operation, die auf einen Anfangsvektor $|v\rangle$ wirkt, eine Summe von Termen mit zunehmenden Potenzen von $H$ ergibt, die auf den Anfangszustand angewendet werden. Es sieht so aus, als könnten wir unseren Krylov-Unterraum einfach durch Zeitentwicklung unseres anfänglichen Schätzzustandes erstellen!

Die Schwierigkeit besteht darin, die Zeitentwicklung auf einem echten Quantencomputer zu realisieren. Viele der Terme in der Hamiltonformel sind nicht miteinander vertauschbar. Während also einige einfache exponentielle Operatoren wie $e^{-iZ}$ einfachen Schaltkreisen entsprechen, tun dies allgemeine Hamiltonianer nicht. Und da sie nicht-kommutierende Terme enthalten, können wir den Exponentialwert nicht einfach in ein Produkt einfacher Terme zerlegen, wie wir es bei Zahlen können.

Dies ist also nicht trivial, aber ein gut untersuchter Prozess in der Quanteninformatik. Wir führen die Zeitevolution auf Quantencomputern mit Hilfe eines Prozesses durch, der Trotterisierung genannt wird, was an sich schon ein umfangreiches Thema ist [10]. Aber auf einer sehr hohen Ebene begrenzen wir die Auswirkungen der Nicht-Kommutativität von Termen, indem wir die zeitliche Entwicklung in sehr kleine Schritte unterteilen, z. B. $m$ mit der Größe $dt$.

wobei $dt = t/m$.

Wir nennen einen Krylov-Unterraum der Ordnung r, den wir im klassischen Kontext mit Potenzen von H direkt erzeugt haben, einen "Potenz-Krylov-Unterraum".

Jetzt erzeugen wir einen ähnlichen Raum mit dem unitären Zeit-Evolutions-Operator $U \equiv e^{-iHt}$; wir nennen ihn den "unitären Krylov-Raum" $\mathcal{K}_U^r$. Der Power-Krylov-Unterraum $\mathcal{K}_P^r$, den wir klassischerweise verwenden, kann nicht direkt auf einem Quantencomputer erzeugt werden, da $H$ kein unitärer Operator ist. Es kann gezeigt werden, dass die Verwendung des unitären Krylov-Unterraums ähnliche Konvergenzgarantien bietet wie der Potenz-Krylov-Unterraum, nämlich dass der Grundzustandsfehler effizient konvergiert, solange der Anfangszustand $|v\rangle$ eine Überlappung mit dem wahren Grundzustand hat, die nicht exponentiell verschwindet, und solange eine ausreichende Lücke zwischen den Eigenwerten besteht. Für eine genauere Diskussion der Konvergenz siehe [1].

Hier werden die Potenzen von $U$ zu verschiedenen Zeitschritten (die Potenz $k^\text{th}$ von $U$ schreitet um eine Zeit $k \times dt$ voran). Wir können das Element des Unterraums, das sich für die Gesamtzeit $k dt$ zeitlich entwickelt, als $|\psi_k\rangle$ bezeichnen.

Wir können unseren Hamiltonian H auf den unitären Krylov-Unterraum $\mathcal{K}_U^r$ projizieren. Mit anderen Worten: Wir berechnen jedes Matrixelement von $H$ in der Basis $\mathcal{K}_U^r$. Wir bezeichnen diese projizierte Matrix als $\tilde{H}$.

3.2 Wie man es auf einem Quantencomputer implementiert

Die Matrixelemente von $\tilde{H}$ sind durch die Erwartungswerte $\langle \psi_m |H| \psi_n\rangle$ gegeben, die mit dem Quantencomputer geschätzt werden können. Beachten Sie, dass $H$ als eine Summe von Pauli-Operatoren auf verschiedenen Qubits geschrieben werden kann und dass nicht alle Pauli-Operatoren gleichzeitig gemessen werden können. Wir können die Pauli-Terme in Gruppen von kommutierenden Termen sortieren und alle diese Terme auf einmal messen. Aber wir brauchen vielleicht viele solcher Gruppen, um alle Begriffe abzudecken. Daher wird die Anzahl der verschiedenen Pendelgruppen, in die die Begriffe unterteilt werden können, $N_\text{GCP}$, wichtig.

Dabei ist $P_\alpha$ ein Pauli-String der Form $P_\alpha \sim IZIXII...YZXIX$ oder eine Menge solcher Pauli-Strings, die miteinander kommutieren. Da wir „ $H$ “ als eine solche Summe messbarer Operatoren darstellen können, lassen sich die folgenden Ausdrücke für Matrixelemente von „ $\tilde{H}$ “ mithilfe des primitiven Schätzers „ IBM Quantum “ realisieren.

Dabei sind $\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle$ die Vektoren des unitären Krylov-Raums und $t_n = n dt$ die Vielfachen des gewählten Zeitschritts $dt$. Auf einem Quantencomputer kann die Berechnung der einzelnen Matrixelemente mit jedem beliebigen Algorithmus erfolgen, der es ermöglicht, eine Überlappung zwischen Quantenzuständen zu erzielen. In dieser Lektion werden wir uns mit dem Hadamard-Test beschäftigen. Da der Hamilton-Operator „ $\mathcal{K}_U$ “ die Dimension „ $r$ “ hat, hat der auf den Unterraum projizierte Hamilton-Operator die Dimension „ $r \times r$ “. Ist „ $r$ “ klein genug (im Allgemeinen reicht „ $r<<100$ “ aus, um die Konvergenz der Schätzungen der Eigenwerte zu gewährleisten), können wir den projizierten Hamilton-Operator „ $\tilde{H}$ “ klassisch problemlos diagonalisieren. Allerdings können wir $\tilde{H}$ nicht direkt diagonalisieren, da die Vektoren des Krylov-Raums nicht orthogonal sind. Wir müssen ihre Überschneidungen messen und eine Matrix erstellen: $\tilde{S}$

Dies ermöglicht es uns, das Eigenwertproblem in einem nicht-orthogonalen Raum zu lösen (auch verallgemeinertes Eigenwertproblem genannt)

Anhand der Lösungen dieses verallgemeinerten Eigenwertproblems können dann die Eigenwerte und Eigenzustände von $H$ geschätzt werden. Die Schätzung der Energie des Grundzustands wird beispielsweise durch die Ermittlung des kleinsten Eigenwerts $E$ und des Grundzustands aus dem entsprechenden Eigenvektor $\vec{c}$ gewonnen. Die Koeffizienten in $\vec{c}$ bestimmen den Beitrag der verschiedenen Vektoren, die $\mathcal{K}_U$ aufspannen.

Verallgemeinertes Eigenwertproblem

Warum können wir nicht einfach $\tilde{H}$ diagonalisieren? Da $\tilde{S}$ die Information über die Geometrie der Krylov-Basis enthält (die in allen außer in sehr speziellen Fällen nicht orthogonal ist), beschreibt $\tilde{H}$ für sich genommen keine Projektion des vollständigen Hamiltonianers, so dass seine Eigenwerte keine besondere Beziehung zu denen des vollständigen Hamiltonianers haben - sie könnten beliebige Werte sein. Die Lösung des verallgemeinerten Eigenwertproblems ist erforderlich, um die approximativen Eigenwerte und Eigenvektoren zu erhalten, die der Projektion des vollständigen Hamiltonianers in den Krylov-Raum entsprechen.

Die Abbildung zeigt eine Schaltungsdarstellung des modifizierten Hadamard-Tests, einer Methode, die zur Berechnung der Überlappung zwischen verschiedenen Quantenzuständen verwendet wird. Für jedes Matrixelement $\tilde{H}_{i,j}$ wird ein Hadamard-Test zwischen den Zuständen $\vert \psi_i \rangle$, $\vert \psi_j \rangle$ durchgeführt. Dies wird in der Abbildung durch das Farbschema für die Matrixelemente und die entsprechenden Operationen $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$ hervorgehoben. Daher ist eine Reihe von Hadamard-Tests für alle möglichen Kombinationen von Krylov-Raum-Vektoren erforderlich, um alle Matrixelemente des projizierten Hamiltonian $\tilde{H}$ zu berechnen. Der oberste Draht in der Hadamard-Testschaltung ist ein Ancilla-Qubit, das entweder in der X- oder Y-Basis gemessen wird; sein Erwartungswert bestimmt den Wert der Überlappung zwischen den Zuständen. Der untere Draht stellt alle Qubits des System-Hamiltonianers dar. Die Operation $\text{Prep} \; \psi_i$ bereitet das System-Qubit auf den Zustand $\vert \psi_i \rangle$ vor, der durch den Zustand des Ancilla-Qubits gesteuert wird (ähnlich für $\text{Prep} \; \psi_j$ ), und die Operation $P$ stellt die Pauli-Zerlegung des System-Hamiltonschen $H = \sum_i P_i$ dar. Die Umsetzung dieser Operation auf einem Quantencomputer wird im Folgenden näher erläutert.

4. Krylov-Quanten-Diagonalisierung auf einem Quantencomputer

Wir werden nun die Krylov-Quanten-Diagonalisierung auf einem echten Quantencomputer implementieren. Beginnen wir damit, einige nützliche Pakete zu importieren.

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")

Wir definieren die folgende Funktion, um das oben erläuterte verallgemeinerte Eigenwertproblem zu lösen.

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]

Zumindest für ein erstes Benchmarking ist es nützlich, eine exakte klassische Lösung zu kennen, um das Konvergenzverhalten zu überprüfen. Die folgende Funktion berechnet die Grundzustandsenergie eines Hamiltonianers, indem sie den Hamiltonianer und die Anzahl der Qubits als Argumente verwendet.

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 Schritt 1: Problem auf Quantenschaltungen und Operatoren abbilden

Jetzt werden wir einen Hamiltonianer definieren. Sie unterscheidet sich von der obigen Funktion dadurch, dass die obige Funktion einen Hamiltonian als Argument nimmt und nur den Grundzustand zurückgibt, und zwar auf klassische Weise. Dieser Hamiltonian, den wir hier definieren, bestimmt die Energieniveaus aller Energieeigenzustände, und dieser Hamiltonian kann mit Pauli-Operatoren konstruiert und auf einem Quantencomputer implementiert werden.

Wir wählen einen Hamiltonian, der einer Kette von Spins entspricht, die eine beliebige Orientierung im Raum haben können, eine sogenannte "Heisenberg-Kette". Wir gehen davon aus, dass der Spin $i^\text{th}$ von seinen nächsten Nachbarn (den Spins $(i-1)^\text{th}$ und $(i+1)^\text{th}$ ) beeinflusst werden kann, nicht aber von weiter entfernten Nachbarspins. Wir berücksichtigen auch die Möglichkeit, dass die Wechselwirkung zwischen Spins unterschiedlich ist, wenn die Spins entlang verschiedener Achsen zeigen. Diese Asymmetrie kann beispielsweise durch die Struktur des Kristallgitters entstehen, in das die Spins eingebettet sind.

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