Shors Algorithmus

Jetzt wenden wir uns dem Problem der ganzzahligen Faktorisierung zu und sehen, wie es auf einem Quantencomputer mithilfe von Phasenschätzung effizient gelöst werden kann. Der Algorithmus, den wir erhalten werden, ist Shors Algorithmus zur ganzzahligen Faktorisierung. Shor hat seinen Algorithmus nicht explizit in Begriffen der Phasenschätzung beschrieben, aber es ist eine natürliche und anschauliche Möglichkeit, seine Funktionsweise zu erklären.

Wir beginnen mit einem Zwischenproblem, dem sogenannten Ordnungsbestimmungsproblem, und sehen, wie die Phasenschätzung eine Lösung dafür liefert. Anschließend sehen wir, wie eine effiziente Lösung des Ordnungsbestimmungsproblems uns eine effiziente Lösung des ganzzahligen Faktorisierungsproblems gibt. (Wenn die Lösung eines Problems die Lösung eines anderen Problems liefert, sagen wir, dass das zweite Problem auf das erste reduziert — in diesem Fall reduzieren wir also die ganzzahlige Faktorisierung auf die Ordnungsbestimmung.) Dieser zweite Teil von Shors Algorithmus macht überhaupt keinen Gebrauch von Quantencomputing; er ist vollständig klassisch. Quantencomputing wird nur für die Ordnungsbestimmung benötigt.

Das Ordnungsbestimmungsproblem

Ein wenig Zahlentheorie

Um das Ordnungsbestimmungsproblem und seine Lösung durch Phasenschätzung zu erklären, hilft es, mit ein paar grundlegenden zahlentheoretischen Konzepten zu beginnen und dabei hilfreiche Notation einzuführen.

Zunächst definieren wir für eine beliebige positive ganze Zahl $N$ die Menge $\mathbb{Z}_N$ wie folgt.

$$\mathbb{Z}_N = \{0,1,\ldots,N-1\}$$

Zum Beispiel: $\mathbb{Z}_1 = \{0\},\;$ $\mathbb{Z}_2 = \{0,1\},\;$ $\mathbb{Z}_3 = \{0,1,2\},\;$ und so weiter.

Das sind Mengen von Zahlen, aber wir können sie als mehr als nur Mengen betrachten. Insbesondere können wir über arithmetische Operationen auf $\mathbb{Z}_N$ nachdenken, wie Addition und Multiplikation — und wenn wir uns darauf einigen, die Ergebnisse stets modulo $N$ zu nehmen (d. h. durch $N$ zu teilen und den Rest als Ergebnis zu verwenden), bleiben wir bei diesen Operationen immer in dieser Menge. Die beiden konkreten Operationen Addition und Multiplikation, beide modulo $N$ genommen, machen $\mathbb{Z}_N$ zu einem Ring, einem grundlegend wichtigen Objekt in der Algebra.

Zum Beispiel sind $3$ und $5$ Elemente von $\mathbb{Z}_7$, und wenn wir sie multiplizieren, erhalten wir $3\cdot 5 = 15$, was bei Division durch $7$ den Rest $1$ lässt. Das schreibt man manchmal so:

$$3 \cdot 5 \equiv 1 \; (\textrm{mod } 7)$$

Wenn klar ist, dass wir in $\mathbb{Z}_7$ arbeiten, kann man aber auch einfach $3 \cdot 5 = 1$ schreiben, um die Notation möglichst einfach zu halten.

Hier sind die Additions- und Multiplikationstabellen für $\mathbb{Z}_6$ als Beispiel.

$$\begin{array}{c|cccccc} + & 0 & 1 & 2 & 3 & 4 & 5 \\\hline 0 & 0 & 1 & 2 & 3 & 4 & 5 \\ 1 & 1 & 2 & 3 & 4 & 5 & 0 \\ 2 & 2 & 3 & 4 & 5 & 0 & 1 \\ 3 & 3 & 4 & 5 & 0 & 1 & 2 \\ 4 & 4 & 5 & 0 & 1 & 2 & 3 \\ 5 & 5 & 0 & 1 & 2 & 3 & 4 \\ \end{array} \qquad \begin{array}{c|cccccc} \cdot & 0 & 1 & 2 & 3 & 4 & 5 \\\hline 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 2 & 3 & 4 & 5 \\ 2 & 0 & 2 & 4 & 0 & 2 & 4 \\ 3 & 0 & 3 & 0 & 3 & 0 & 3 \\ 4 & 0 & 4 & 2 & 0 & 4 & 2 \\ 5 & 0 & 5 & 4 & 3 & 2 & 1 \\ \end{array}$$

Unter den $N$ Elementen von $\mathbb{Z}_N$ sind die Elemente $a\in\mathbb{Z}_N$ mit $\gcd(a,N) = 1$ besonders. Die Menge dieser Elemente wird häufig mit einem Stern bezeichnet:

$$\mathbb{Z}_N^{\ast} = \{a\in \mathbb{Z}_N : \gcd(a,N) = 1\}$$

Betrachtet man nur die Multiplikation, bildet $\mathbb{Z}_N^{\ast}$ eine Gruppe — genauer eine abelsche Gruppe — ein weiteres wichtiges Objekt in der Algebra. Eine grundlegende Eigenschaft dieser Mengen (und endlicher Gruppen im Allgemeinen) ist, dass man für jedes $a\in\mathbb{Z}_N^{\ast}$, wenn man $a$ wiederholt mit sich selbst multipliziert, schließlich immer die Zahl $1$ erhält.

Als erstes Beispiel nehmen wir $N=6$. Es gilt $5\in\mathbb{Z}_6^{\ast}$, da $\gcd(5,6) = 1$, und wenn wir $5$ mit sich selbst multiplizieren, erhalten wir $1$, wie die obige Tabelle bestätigt.

$$5^2 = 1 \quad \text{(in $\mathbb{Z}_6$)}$$

Als zweites Beispiel nehmen wir $N = 21$. Die Zahlen von $0$ bis $20$, die mit $21$ den ggT $1$ haben, sind:

$$\mathbb{Z}_{21}^{\ast} = \{1,2,4,5,8,10,11,13,16,17,19,20\}$$

Für jedes dieser Elemente kann man die Zahl zu einer positiven ganzzahligen Potenz erheben und erhält $1$. Hier sind die kleinsten Potenzen, für die das funktioniert:

$$\begin{array}{ccc} 1^{1} = 1 \quad & 8^{2} = 1 \quad & 16^{3} = 1 \\[1mm] 2^{6} = 1 \quad & 10^{6} = 1 \quad & 17^{6} = 1 \\[1mm] 4^{3} = 1 \quad & 11^{6} = 1 \quad & 19^{6} = 1 \\[1mm] 5^{6} = 1 \quad & 13^{2} = 1 \quad & 20^{2} = 1 \end{array}$$

Wir arbeiten für alle diese Gleichungen in $\mathbb{Z}_{21}$, was wir nicht explizit hinschreiben — es ist implizit gemeint, um die Notation übersichtlich zu halten. Das werden wir im Rest der Lektion beibehalten.

Problemformulierung und Verbindung zur Phasenschätzung

Jetzt können wir das Ordnungsbestimmungsproblem formulieren.

Ordnungsbestimmung

Eingabe: positive ganze Zahlen $N$ und $a$ mit $\gcd(N,a) = 1$
Ausgabe: die kleinste positive ganze Zahl $r$ mit $a^r \equiv 1$ $(\textrm{mod } N)$

Anders ausgedrückt: In der Notation von oben ist $a \in \mathbb{Z}_N^{\ast}$ gegeben, und wir suchen die kleinste positive ganze Zahl $r$ mit $a^r = 1$. Diese Zahl $r$ heißt die Ordnung von $a$ modulo $N$.

Um das Ordnungsbestimmungsproblem mit der Phasenschätzung zu verbinden, betrachten wir die Operation auf einem System, dessen klassische Zustände $\mathbb{Z}_N$ entsprechen, bei der wir mit einem festen Element $a\in\mathbb{Z}_N^{\ast}$ multiplizieren.

$$M_a \vert x\rangle = \vert ax \rangle \qquad \text{(für jedes $x\in\mathbb{Z}_N$)}$$

Um das klarzustellen: Die Multiplikation findet in $\mathbb{Z}_N$ statt, also ist es implizit, dass wir das Produkt modulo $N$ im Ket auf der rechten Seite der Gleichung nehmen.

Nehmen wir zum Beispiel $N = 15$ und $a=2$. Die Wirkung von $M_2$ auf die Standardbasis $\{\vert 0\rangle,\ldots,\vert 14\rangle\}$ ist:

$$\begin{array}{ccc} M_{2} \vert 0 \rangle = \vert 0\rangle \quad & M_{2} \vert 5 \rangle = \vert 10\rangle \quad & M_{2} \vert 10 \rangle = \vert 5\rangle \\[1mm] M_{2} \vert 1 \rangle = \vert 2\rangle \quad & M_{2} \vert 6 \rangle = \vert 12\rangle \quad & M_{2} \vert 11 \rangle = \vert 7\rangle \\[1mm] M_{2} \vert 2 \rangle = \vert 4\rangle \quad & M_{2} \vert 7 \rangle = \vert 14\rangle \quad & M_{2} \vert 12 \rangle = \vert 9\rangle \\[1mm] M_{2} \vert 3 \rangle = \vert 6\rangle \quad & M_{2} \vert 8 \rangle = \vert 1\rangle \quad & M_{2} \vert 13 \rangle = \vert 11\rangle \\[1mm] M_{2} \vert 4 \rangle = \vert 8\rangle \quad & M_{2} \vert 9 \rangle = \vert 3\rangle \quad & M_{2} \vert 14 \rangle = \vert 13\rangle \end{array}$$

Diese Operation ist unitär, sofern $\gcd(a,N)=1$; sie permutiert die Elemente der Standardbasis $\{\vert 0\rangle,\ldots,\vert N-1\rangle\}$, ist also als Matrix eine Permutationsmatrix. Aus ihrer Definition ist offensichtlich, dass sie deterministisch ist, und eine einfache Möglichkeit, ihre Invertierbarkeit einzusehen, besteht darin, die Ordnung $r$ von $a$ modulo $N$ zu betrachten: Die Inverse von $M_a$ ist $M_a^{r-1}$.

$$M_a^{r-1} M_a = M_a^r = M_{a^r} = M_1 = \mathbb{I}$$

Es gibt noch eine andere Sichtweise auf die Inverse, die kein Wissen über $r$ voraussetzt (was schließlich das ist, was wir berechnen wollen). Für jedes Element $a\in\mathbb{Z}_N^{\ast}$ gibt es immer ein eindeutiges Element $b\in\mathbb{Z}_N^{\ast}$ mit $ab=1$. Dieses Element $b$ bezeichnen wir mit $a^{-1}$, und es kann effizient berechnet werden; eine Erweiterung von Euklids ggT-Algorithmus erledigt das mit quadratischem Aufwand in $\operatorname{lg}(N)$. Damit gilt:

$$M_{a^{-1}} M_a = M_{a^{-1}a} = M_1 = \mathbb{I}.$$

Die Operation $M_a$ ist also sowohl deterministisch als auch invertierbar. Das bedeutet, dass sie durch eine Permutationsmatrix beschrieben wird und daher unitär ist.

Betrachten wir nun die Eigenvektoren und Eigenwerte der Operation $M_a$, unter der Annahme $a\in\mathbb{Z}_N^{\ast}$. Wie gerade gezeigt, bedeutet diese Annahme, dass $M_a$ unitär ist.

Es gibt $N$ Eigenwerte von $M_a$, möglicherweise inklusive Wiederholungen, und es gibt gewisse Freiheit bei der Wahl zugehöriger Eigenvektoren — aber wir brauchen uns nicht um alle Möglichkeiten zu kümmern. Beginnen wir einfach und identifizieren wir zunächst nur einen Eigenvektor von $M_a$.

$$\vert \psi_0 \rangle = \frac{\vert 1 \rangle + \vert a \rangle + \cdots + \vert a^{r-1} \rangle}{\sqrt{r}}$$

Die Zahl $r$ ist die Ordnung von $a$ modulo $N$, hier und im weiteren Verlauf der Lektion. Der zugehörige Eigenwert ist $1$, da der Zustand bei Multiplikation mit $a$ unverändert bleibt.

$$M_a \vert \psi_0 \rangle = \frac{\vert a \rangle + \cdots + \vert a^{r-1} \rangle + \vert a^r \rangle}{\sqrt{r}} = \frac{\vert a \rangle + \cdots + \vert a^{r-1} \rangle + \vert 1 \rangle}{\sqrt{r}} = \vert \psi_0 \rangle$$

Das passiert, weil $a^r = 1$: jeder Standardbasiszustand $\vert a^k \rangle$ wird für $k\leq r-1$ auf $\vert a^{k+1} \rangle$ verschoben, und $\vert a^{r-1} \rangle$ wird zurück auf $\vert 1\rangle$ verschoben. Bildlich gesprochen rühren wir in $\vert \psi_0 \rangle$, aber da es bereits vollständig „durchgerührt" ist, ändert sich nichts.

Hier ist ein weiterer Eigenvektor von $M_a$. Dieser ist im Kontext von Ordnungsbestimmung und Phasenschätzung interessanter.

$$\vert \psi_1 \rangle = \frac{\vert 1 \rangle + \omega_r^{-1} \vert a \rangle + \cdots + \omega_r^{-(r-1)}\vert a^{r-1} \rangle}{\sqrt{r}}$$

Alternativ kann man diesen Vektor mit einer Summe schreiben:

$$\vert \psi_1 \rangle = \frac{1}{\sqrt{r}} \sum_{k = 0}^{r-1} \omega_r^{-k} \vert a^k \rangle$$

Hier taucht die komplexe Zahl $\omega_r = e^{2\pi i/r}$ auf natürliche Weise auf, bedingt durch die Weise, wie die Multiplikation mit $a$ modulo $N$ funktioniert. Der zugehörige Eigenwert ist dieses Mal $\omega_r$. Das sehen wir folgendermaßen:

$$M_a \vert \psi_1 \rangle = \frac{1}{\sqrt{r}}\sum_{k = 0}^{r-1} \omega_r^{-k} M_a\vert a^k \rangle = \frac{1}{\sqrt{r}}\sum_{k = 0}^{r-1} \omega_r^{-k} \vert a^{k+1} \rangle = \frac{1}{\sqrt{r}}\sum_{k = 1}^{r} \omega_r^{-(k - 1)} \vert a^{k} \rangle = \frac{1}{\sqrt{r}}\omega_r \sum_{k = 1}^{r} \omega_r^{-k} \vert a^{k} \rangle$$

Da $\omega_r^{-r} = 1 = \omega_r^0$ und $\vert a^r \rangle = \vert 1\rangle = \vert a^0\rangle$, ergibt sich:

$$\frac{1}{\sqrt{r}}\sum_{k = 1}^{r} \omega_r^{-k} \vert a^{k} \rangle = \frac{1}{\sqrt{r}}\sum_{k = 0}^{r-1} \omega_r^{-k} \vert a^k \rangle = \vert\psi_1\rangle,$$

also $M_a \vert\psi_1\rangle = \omega_r \vert\psi_1\rangle$.

Mit demselben Argument können wir weitere Eigenvektor/Eigenwert-Paare für $M_a$ identifizieren. Für jede Wahl von $j\in\{0,\ldots,r-1\}$ gilt:

$$\vert \psi_j \rangle = \frac{1}{\sqrt{r}} \sum_{k = 0}^{r-1} \omega_r^{-jk} \vert a^k \rangle$$

ist ein Eigenvektor von $M_a$ mit dem zugehörigen Eigenwert $\omega_r^j$.

$$M_a \vert \psi_j \rangle = \omega_r^j \vert \psi_j \rangle$$

Es gibt noch weitere Eigenvektoren von $M_a$, aber wir müssen uns nicht um sie kümmern — wir beschränken uns auf die Eigenvektoren $\vert\psi_0\rangle,\ldots,\vert\psi_{r-1}\rangle$, die wir gerade identifiziert haben.

Ordnungsbestimmung durch Phasenschätzung

Um das Ordnungsbestimmungsproblem für eine gegebene Wahl $a\in\mathbb{Z}_N^{\ast}$ zu lösen, können wir das Phasenschätzungsverfahren auf die Operation $M_a$ anwenden.

Dazu müssen wir nicht nur $M_a$ effizient als Quantencircuit implementieren, sondern auch $M_a^2,$ $M_a^4,$ $M_a^8$ und so weiter, so weit wie nötig für eine ausreichend genaue Schätzung. Hier erklären wir, wie das gemacht werden kann, und wir werden später genau herausfinden, wie viel Präzision benötigt wird.

Fangen wir mit der Operation $M_a$ an sich an. Da wir natürlich im Quantencircuit-Modell arbeiten, verwenden wir Binärdarstellung, um die Zahlen zwischen $0$ und $N-1$ zu kodieren. Die größte Zahl, die wir kodieren müssen, ist $N-1$, die Anzahl der benötigten Bits ist also:

$$n = \operatorname{lg}(N-1) = \lfloor \log(N-1) \rfloor + 1.$$

Für $N = 21$ zum Beispiel gilt $n = \operatorname{lg}(N-1) = 5$. Hier ist die Kodierung der Elemente von $\mathbb{Z}_{21}$ als Binärstrings der Länge $5$:

$$\begin{gathered} 0 \mapsto 00000\\[1mm] 1 \mapsto 00001\\[1mm] \vdots\\[1mm] 20 \mapsto 10100 \end{gathered}$$

Und hier ist eine präzise Definition, wie $M_a$ als $n$-Qubit-Operation definiert wird:

$$M_a \vert x\rangle = \begin{cases} \vert ax \; (\textrm{mod}\;N)\rangle & 0\leq x < N\\[1mm] \vert x\rangle & N\leq x < 2^n \end{cases}$$

Es geht darum, dass wir zwar nur interessiert sind, wie $M_a$ auf $\vert 0\rangle,\ldots,\vert N-1\rangle$ wirkt, aber dennoch festlegen müssen, wie es auf die verbleibenden $2^n - N$ Standardbasiszustände wirkt — und das so, dass die Operation unitär bleibt. $M_a$ so zu definieren, dass es auf den verbleibenden Standardbasiszuständen nichts tut, erfüllt diese Anforderung.

Mithilfe der Algorithmen für ganzzahlige Multiplikation und Division, die in der vorherigen Lektion besprochen wurden, sowie der Methodik für reversible, müllfreie Implementierungen davon, können wir einen Quantencircuit für $M_a$ bauen — für jede Wahl von $a\in\mathbb{Z}_N^{\ast}$ — mit Aufwand $O(n^2)$. Eine Möglichkeit ist folgende:

1. Baue einen Circuit für die Operation

   $$\vert x \rangle \vert y \rangle \mapsto \vert x \rangle \vert y \oplus f_a(x)\rangle$$

   wobei

   $$f_a(x) = \begin{cases} ax \; (\textrm{mod}\;N) & 0\leq x < N\\[1mm] x & N\leq x < 2^n \end{cases}$$

   nach der in der vorherigen Lektion beschriebenen Methode. Das ergibt einen Circuit der Größe $O(n^2)$.

2. Tausche die beiden $n$-Qubit-Systeme mithilfe von $n$ Swap-Gates aus, indem du die Qubits einzeln tauschst.

3. Analog zum ersten Schritt, baue einen Circuit für die Operation

   $$\vert x \rangle \vert y \rangle \mapsto \vert x \rangle \bigl\vert y \oplus f_{a^{-1}}(x)\bigr\rangle$$

   wobei $a^{-1}$ die Inverse von $a$ in $\mathbb{Z}_N^{\ast}$ ist.

Durch Initialisierung der unteren $n$ Qubits und Komposition der drei Schritte erhält man folgende Transformation:

$$\vert x \rangle \vert 0^n \rangle \stackrel{\text{Schritt 1}}{\mapsto} \vert x \rangle \vert f_a(x)\rangle \stackrel{\text{Schritt 2}}{\mapsto} \vert f_a(x)\rangle \vert x \rangle \stackrel{\text{Schritt 3}}{\mapsto} \vert f_a(x)\rangle \bigl\vert x \oplus f_{a^{-1}}(f_a(x)) \bigr\rangle = \vert f_a(x)\rangle\vert 0^n \rangle$$

Die Methode benötigt Hilfs-Qubits, aber sie werden am Ende in ihren initialisierten Zustand zurückversetzt, was es uns erlaubt, diese Circuits für die Phasenschätzung zu verwenden. Der Gesamtaufwand des Circuits beträgt $O(n^2)$.

Um $M_a^2,$ $M_a^4,$ $M_a^8$ usw. zu implementieren, kann genau dieselbe Methode verwendet werden, außer dass wir $a$ durch $a^2,$ $a^4,$ $a^8$ usw. als Elemente von $\mathbb{Z}_N^{\ast}$ ersetzen. Das heißt: Für jede gewählte Potenz $k$ können wir einen Circuit für $M_a^k$ erstellen, nicht indem wir den Circuit für $M_a$ $k$-mal iterieren, sondern indem wir $b = a^k \in \mathbb{Z}_N^{\ast}$ berechnen und dann den Circuit für $M_b$ verwenden.

Die Berechnung von Potenzen $a^k \in \mathbb{Z}_N$ ist das modulare Exponentierungsproblem, das in der vorherigen Lektion erwähnt wurde. Diese Berechnung kann klassisch durchgeführt werden, mithilfe des dort erwähnten Algorithmus für modulare Exponentierung (in der rechnerischen Zahlentheorie oft Potenzalgorithmus genannt). Tatsächlich benötigen wir nur Zweierpotenzen von $a$, nämlich $a^2, a^4, \ldots a^{2^{m-1}} \in \mathbb{Z}_N^{\ast}$, und wir erhalten diese Potenzen durch iteriertes Quadrieren, $m-1$ Mal. Jedes Quadrieren kann durch einen booleschen Circuit der Größe $O(n^2)$ durchgeführt werden.

Im Wesentlichen verlagern wir hier das Problem, $M_a$ bis zu $2^{m-1}$ Mal zu iterieren, auf eine effiziente klassische Berechnung. Und es ist ein glücklicher Umstand, dass das möglich ist! Bei einer beliebigen Wahl eines Quantencircuits im Phasenschätzungsproblem dürfte das nicht möglich sein — in diesem Fall wächst der Aufwand für die Phasenschätzung exponentiell in der Anzahl der Kontroll-Qubits $m$.

Lösung mit einem geeigneten Eigenvektor

Um zu verstehen, wie wir das Ordnungsbestimmungsproblem durch Phasenschätzung lösen können, nehmen wir zunächst an, dass wir das Phasenschätzungsverfahren auf die Operation $M_a$ mit dem Eigenvektor $\vert\psi_1\rangle$ anwenden. Diesen Eigenvektor zu beschaffen ist — wie sich herausstellt — nicht einfach, daher ist das nicht das Ende der Geschichte — aber es ist hilfreich, hier zu beginnen.

Der Eigenwert von $M_a$, der zum Eigenvektor $\vert \psi_1\rangle$ gehört, ist

$$\omega_r = e^{2\pi i \frac{1}{r}}.$$

Das heißt, $\omega_r = e^{2\pi i \theta}$ für $\theta = 1/r$. Wenn wir also das Phasenschätzungsverfahren auf $M_a$ mit dem Eigenvektor $\vert\psi_1\rangle$ anwenden, erhalten wir eine Approximation von $1/r$. Durch Berechnung des Kehrwerts können wir dann $r$ bestimmen — vorausgesetzt, unsere Approximation ist gut genug.

Genauer gesagt: Wenn wir das Phasenschätzungsverfahren mit $m$ Kontroll-Qubits durchführen, erhalten wir eine Zahl $y\in\{0,\ldots,2^m-1\}$. Wir nehmen $y/2^m$ als Schätzung für $\theta$, also $1/r$ in unserem Fall. Um $r$ aus dieser Approximation zu bestimmen, liegt es nahe, den Kehrwert der Schätzung zu berechnen und zur nächsten ganzen Zahl zu runden:

$$\left\lfloor \frac{2^m}{y} + \frac{1}{2} \right\rfloor$$

Nehmen wir zum Beispiel $r = 6$ und führen wir die Phasenschätzung auf $M_a$ mit dem Eigenvektor $\vert\psi_1\rangle$ und $m = 5$ Kontroll-Bits durch. Die beste 5-Bit-Approximation an $1/r = 1/6$ ist $5/32$, und wir erhalten mit ziemlich hoher Wahrscheinlichkeit (in diesem Fall etwa $68\%$) das Ergebnis $y=5$ aus der Phasenschätzung. Es gilt:

$$\frac{2^m}{y} = \frac{32}{5} = 6{,}4,$$

und Runden zur nächsten ganzen Zahl ergibt $6$, was die richtige Antwort ist.

Andererseits erhalten wir möglicherweise nicht die richtige Antwort, wenn wir nicht genug Präzision verwenden. Wenn wir zum Beispiel $m = 4$ Kontroll-Qubits bei der Phasenschätzung nehmen, erhalten wir möglicherweise die beste 4-Bit-Approximation an $1/r = 1/6$, nämlich $3/16$. Der Kehrwert ergibt:

$$\frac{2^m}{y} = \frac{16}{3} = 5{,}333 \cdots$$

und Runden zur nächsten ganzen Zahl ergibt die falsche Antwort $5$.

Wie viel Präzision benötigen wir also für die richtige Antwort? Wir wissen, dass die Ordnung $r$ eine ganze Zahl ist, und intuitiv brauchen wir genug Präzision, um $1/r$ von benachbarten Möglichkeiten zu unterscheiden, darunter $1/(r+1)$ und $1/(r-1)$. Die nächstgelegene Zahl zu $1/r$ ist $1/(r+1)$, und der Abstand zwischen diesen beiden Zahlen beträgt:

$$\frac{1}{r} - \frac{1}{r+1} = \frac{1}{r(r+1)}.$$

Um sicherzustellen, dass wir $1/r$ nicht mit $1/(r+1)$ verwechseln, reicht es daher aus, genug Präzision zu verwenden, damit eine beste Approximation $y/2^m$ an $1/r$ näher an $1/r$ liegt als an $1/(r+1)$. Wenn wir genug Genauigkeit verwenden, sodass

$$\left\vert \frac{y}{2^m} - \frac{1}{r} \right\vert < \frac{1}{2 r (r+1)},$$

der Fehler also kleiner als die Hälfte des Abstands zwischen $1/r$ und $1/(r+1)$ ist, dann liegt $y/2^m$ näher an $1/r$ als an jeder anderen Möglichkeit, einschließlich $1/(r+1)$ und $1/(r-1)$.

Zur Überprüfung: Angenommen,

$$\frac{y}{2^m} = \frac{1}{r} + \varepsilon$$

für $\varepsilon$ mit

$$\vert\varepsilon\vert < \frac{1}{2 r (r+1)}.$$

Beim Berechnen des Kehrwerts erhalten wir:

$$\frac{2^m}{y} = \frac{1}{\frac{1}{r} + \varepsilon} = \frac{r}{1+\varepsilon r} = r - \frac{\varepsilon r^2}{1+\varepsilon r}.$$

Durch Maximierung im Zähler und Minimierung im Nenner können wir die Abweichung von $r$ begrenzen:

$$\left\vert \frac{\varepsilon r^2}{1+\varepsilon r} \right\vert \leq \frac{ \frac{r^2}{2 r(r+1)}}{1 - \frac{r}{2r(r+1)}} = \frac{r}{2 r + 1} < \frac{1}{2}$$

Wir weichen weniger als $1/2$ von $r$ ab, sodass wir beim Runden $r$ erhalten, wie erwartet.

Leider können wir, da wir $r$ noch nicht kennen, es nicht verwenden, um zu bestimmen, wie viel Genauigkeit wir benötigen. Stattdessen können wir die Tatsache nutzen, dass $r$ kleiner als $N$ sein muss, um sicherzustellen, dass wir genug Präzision verwenden. Genauer: Wenn wir genug Genauigkeit verwenden, um sicherzustellen, dass die beste Approximation $y/2^m$ an $1/r$ erfüllt:

$$\left\vert \frac{y}{2^m} - \frac{1}{r} \right\vert \leq \frac{1}{2N^2},$$

dann haben wir genug Präzision, um $r$ bei Berechnung des Kehrwerts korrekt zu bestimmen. Die Wahl $m = 2\operatorname{lg}(N)+1$ stellt sicher, dass wir mit der oben beschriebenen Methode eine hohe Chance haben, eine Schätzung mit dieser Präzision zu erhalten. (Die Wahl $m = 2\operatorname{lg}(N)$ reicht aus, wenn man mit einer unteren Schranke von 40% für die Erfolgswahrscheinlichkeit zufrieden ist.)

Allgemeine Lösung

Wie wir gerade gesehen haben, können wir $r$ durch Phasenschätzung bestimmen, wenn wir den Eigenvektor $\vert \psi_1 \rangle$ von $M_a$ besitzen — vorausgesetzt, wir verwenden genug Kontroll-Qubits für ausreichende Präzision. Leider ist es nicht einfach, diesen Eigenvektor zu beschaffen, daher müssen wir einen Weg finden, ohne ihn auszukommen.

Angenommen, wir gehen genauso vor wie oben, aber mit dem Eigenvektor $\vert\psi_k\rangle$ anstelle von $\vert\psi_1\rangle$, für eine beliebige Wahl von $k\in\{0,\ldots,r-1\}$. Das Ergebnis des Phasenschätzungsverfahrens ist eine Approximation:

$$\frac{y}{2^m} \approx \frac{k}{r}.$$

Unter der Annahme, dass wir weder $k$ noch $r$ kennen, kann das uns $r$ verraten oder auch nicht. Wenn zum Beispiel $k = 0$, erhalten wir eine Approximation $y/2^m$ an $0$, was uns leider nichts sagt. Das ist jedoch ein Sonderfall; für andere Werte von $k$ können wir zumindest etwas über $r$ lernen.

Wir können einen Algorithmus namens Kettenbruchalgorithmus verwenden, um aus unserer Approximation $y/2^m$ benachbarte Brüche zu berechnen — einschließlich $k/r$, wenn die Approximation gut genug ist. Den Kettenbruchalgorithmus erklären wir hier nicht. Stattdessen geben wir eine bekannte Aussage über diesen Algorithmus an.

Fakt

Für eine ganze Zahl $N\geq 2$ und eine reelle Zahl $\alpha\in(0,1)$ gibt es höchstens eine Wahl von ganzen Zahlen $u,v\in\{0,\ldots,N-1\}$ mit $v\neq 0$ und $\gcd(u,v)=1$, die $\vert \alpha - u/v\vert < \frac{1}{2N^2}$ erfüllen. Gegeben $\alpha$ und $N$, findet der Kettenbruchalgorithmus $u$ und $v$, oder meldet, dass sie nicht existieren. Dieser Algorithmus kann als boolescher Circuit der Größe $O((\operatorname{lg}(N))^3)$ implementiert werden.

Haben wir eine sehr genaue Approximation $y/2^m$ an $k/r$ und führen den Kettenbruchalgorithmus für $N$ und $\alpha = y/2^m$ aus, erhalten wir $u$ und $v$ wie im Fakt beschrieben. Eine Analyse des Fakts erlaubt uns zu schließen, dass:

$$\frac{u}{v} = \frac{k}{r}.$$

Beachte insbesondere, dass wir nicht notwendigerweise $k$ und $r$ bestimmen, sondern nur $k/r$ in gekürzter Form.

Wie bereits festgestellt, werden wir aus $k=0$ nichts lernen. Aber das ist der einzige Wert von $k$, bei dem das passiert. Wenn $k$ von null verschieden ist, kann es gemeinsame Teiler mit $r$ haben, aber die Zahl $v$ aus dem Kettenbruchalgorithmus muss zumindest ein Teiler von $r$ sein.

Es ist nicht offensichtlich, aber es stimmt: Wenn wir $u$ und $v$ für $u/v = k/r$ mit gleichmäßig zufällig gewähltem $k\in\{0,\ldots,r-1\}$ bestimmen können, können wir $r$ mit hoher Wahrscheinlichkeit nach nur wenigen Stichproben ermitteln. Insbesondere: Wenn unsere Schätzung für $r$ das kleinste gemeinsame Vielfache aller beobachteten Nennerwerte $v$ ist, sind wir mit hoher Wahrscheinlichkeit richtig. Intuitiv sind einige Werte von $k$ problematisch, da sie gemeinsame Teiler mit $r$ haben, und diese gemeinsamen Teiler werden versteckt, wenn wir $u$ und $v$ erhalten. Aber zufällige Wahlen von $k$ verbergen Faktoren von $r$ nicht lange, und die Wahrscheinlichkeit, dass wir $r$ nicht korrekt bestimmen, indem wir das kleinste gemeinsame Vielfache der beobachteten Nenner nehmen, fällt exponentiell in der Anzahl der Stichproben.

Es bleibt die Frage, wie wir einen Eigenvektor $\vert\psi_k\rangle$ von $M_a$ erhalten, auf dem wir das Phasenschätzungsverfahren ausführen. Wie sich herausstellt, müssen wir ihn gar nicht explizit erzeugen!

Stattdessen führen wir das Phasenschätzungsverfahren auf dem Zustand $\vert 1\rangle$ aus — gemeint ist die $n$-Bit-Binärkodierung der Zahl $1$ — anstelle eines Eigenvektors $\vert\psi\rangle$ von $M_a$. Bislang haben wir das Phasenschätzungsverfahren nur für einen konkreten Eigenvektor beschrieben, aber nichts hindert uns daran, es auf einem Eingabezustand auszuführen, der kein Eigenvektor von $M_a$ ist — und genau das tun wir hier mit dem Zustand $\vert 1\rangle$. (Das ist kein Eigenvektor von $M_a$, es sei denn $a=1$, was uns nicht interessiert.)

Der Grund für die Wahl von $\vert 1\rangle$ anstelle eines Eigenvektors von $M_a$ ist, dass die folgende Gleichung gilt:

$$\vert 1\rangle = \frac{1}{\sqrt{r}} \sum_{k = 0}^{r-1} \vert \psi_k\rangle$$

Eine Möglichkeit, diese Gleichung zu überprüfen, besteht darin, die Skalarprodukte beider Seiten mit jedem Standardbasiszustand zu vergleichen und dabei frühere Formeln der Lektion zu verwenden. Dadurch erhalten wir genau dieselben Messergebnisse, als hätten wir $k\in\{0,\ldots,r-1\}$ gleichmäßig zufällig gewählt und $\vert\psi_k\rangle$ als Eigenvektor verwendet.

Im Detail: Stellen wir uns vor, wir führen das Phasenschätzungsverfahren mit dem Zustand $\vert 1\rangle$ anstelle eines der Eigenvektoren $\vert\psi_k\rangle$ aus. Nach der inversen Quanten-Fourier-Transformation verbleiben wir im Zustand:

$$\frac{1}{\sqrt{r}} \sum_{k = 0}^{r-1} \vert \psi_k\rangle \vert \gamma_k\rangle,$$

wobei

$$\vert\gamma_k\rangle = \frac{1}{2^m} \sum_{y=0}^{2^m - 1} \sum_{x=0}^{2^m-1} e^{2\pi i x (k/r - y/2^m)} \vert y\rangle.$$

Der Vektor $\vert\gamma_k\rangle$ beschreibt den Zustand der oberen $m$ Qubits nach der inversen Quanten-Fourier-Transformation.

Da $\{\vert\psi_0\rangle,\ldots,\vert\psi_{r-1}\rangle\}$ eine Orthonormalbasis ist, liefert eine Messung der oberen $m$ Qubits eine Approximation $y/2^m$ an den Wert $k/r$, wobei $k\in\{0,\ldots,r-1\}$ gleichmäßig zufällig gewählt wird. Wie bereits diskutiert, erlaubt uns das, $r$ mit hoher Zuverlässigkeit nach mehreren unabhängigen Läufen zu bestimmen — was unser Ziel war.

Gesamtkosten

Der Aufwand zur Implementierung jedes gesteuerten Unitären $M_a^k$ beträgt $O(n^2)$. Es gibt $m$ gesteuerte unitäre Operationen, und da $m = O(n)$, beträgt der Gesamtaufwand für die gesteuerten unitären Operationen $O(n^3)$. Hinzu kommen $m$ Hadamard-Gates (die $O(n)$ zum Aufwand beitragen) und die inverse Quanten-Fourier-Transformation, die $O(n^2)$ beiträgt. Damit dominieren die gesteuerten unitären Operationen den Gesamtaufwand des Verfahrens — der also $O(n^3)$ beträgt.

Neben dem Quantencircuit selbst gibt es einige klassische Berechnungen, die unterwegs durchgeführt werden müssen. Das umfasst die Berechnung der Potenzen $a^k$ in $\mathbb{Z}_N$ für $k = 2, 4, 8, \ldots, 2^{m-1}$, die zur Erstellung der gesteuerten unitären Gates benötigt werden, sowie den Kettenbruchalgorithmus, der Approximationen von $\theta$ in Brüche umwandelt. Diese Berechnungen können durch boolesche Circuits mit einem Gesamtaufwand von $O(n^3)$ durchgeführt werden.

Wie üblich können all diese Schranken durch asymptotisch schnellere Algorithmen verbessert werden; diese Schranken setzen Standard-Algorithmen für grundlegende arithmetische Operationen voraus.

Faktorisierung durch Ordnungsbestimmung

Das Letzte, was wir besprechen müssen, ist, wie die Lösung des Ordnungsbestimmungsproblems bei der Faktorisierung hilft. Dieser Teil ist vollständig klassisch — er hat nichts mit Quantencomputing zu tun.

Die Grundidee ist folgende. Wir möchten die Zahl $N$ faktorisieren und können das rekursiv tun. Konkret können wir uns auf die Aufgabe konzentrieren, $N$ zu teilen, d. h. zwei ganze Zahlen $b,c\geq 2$ zu finden, für die $N = bc$ gilt. Das ist nicht möglich, wenn $N$ eine Primzahl ist, aber wir können mit einem Primzahltest zunächst effizient prüfen, ob $N$ prim ist, und wenn nicht, versuchen wir $N$ zu teilen. Sobald wir $N$ geteilt haben, können wir einfach rekursiv auf $b$ und $c$ vorgehen, bis alle Faktoren prim sind und wir die Primfaktorzerlegung von $N$ erhalten.

Gerade Zahlen zu teilen ist einfach: wir geben einfach $2$ und $N/2$ aus.

Auch vollständige Potenzen lassen sich leicht teilen — also Zahlen der Form $N = s^j$ für ganze Zahlen $s,j\geq 2$ — indem man die Wurzeln $N^{1/2},$ $N^{1/3},$ $N^{1/4}$ usw. approximiert und nahe gelegene ganze Zahlen als Kandidaten für $s$ prüft. Man muss dabei nicht weiter als $\log(N)$ Schritte in dieser Folge gehen, denn ab dann fällt die Wurzel unter $2$ und liefert keine weiteren Kandidaten.

Es ist gut, dass wir beides können, denn die Ordnungsbestimmung hilft weder bei der Faktorisierung gerader Zahlen noch bei Primzahlpotenzen, wo $s$ eine Primzahl ist. Ist $N$ jedoch ungerade und keine Primzahlpotenz, erlaubt uns die Ordnungsbestimmung, $N$ zu teilen.

Probabilistischer Algorithmus zum Teilen einer ungeraden zusammengesetzten ganzen Zahl N, die keine Primzahlpotenz ist

1. Wähle zufällig $a\in\{2,\ldots,N-1\}$.

2. Berechne $d=\gcd(a,N)$.

3. Falls $d > 1$, gib $b = d$ und $c = N/d$ aus und stoppe. Andernfalls fahre weiter, in dem Wissen, dass $a\in\mathbb{Z}_N^{\ast}$.

4. Sei $r$ die Ordnung von $a$ modulo $N$. (Hier wird die Ordnungsbestimmung benötigt.)

5. Falls $r$ gerade ist:

   5.1 Berechne $x = a^{r/2} - 1$ modulo $N$
   5.2 Berechne $d = \gcd(x,N).$
   5.3 Falls $d>1$, gib $b=d$ und $c = N/d$ aus und stoppe.

6. Wurde dieser Punkt erreicht, konnte der Algorithmus keinen Faktor von $N$ finden.

Ein Durchlauf dieses Algorithmus kann scheitern, einen Faktor von $N$ zu finden. Das passiert genau in zwei Situationen:

• Die Ordnung von $a$ modulo $N$ ist ungerade.
• Die Ordnung von $a$ modulo $N$ ist gerade und $\gcd\bigl(a^{r/2} - 1, N\bigr) = 1$.

Mit elementarer Zahlentheorie lässt sich beweisen, dass für eine zufällige Wahl von $a$ mit Wahrscheinlichkeit mindestens $1/2$ keines dieser Ereignisse eintritt. Genauer gesagt ist die Wahrscheinlichkeit, dass eines der Ereignisse eintritt, höchstens $2^{-(m-1)}$, wobei $m$ die Anzahl der verschiedenen Primfaktoren von $N$ ist — deshalb ist die Annahme, dass $N$ keine Primzahlpotenz ist, erforderlich. (Die Annahme, dass $N$ ungerade ist, wird für die Gültigkeit dieser Aussage ebenfalls benötigt.)

Das bedeutet, dass jeder Durchlauf mit mindestens 50% Wahrscheinlichkeit $N$ teilt. Wenn wir den Algorithmus $t$ Mal ausführen, wobei wir jedes Mal $a$ zufällig wählen, werden wir $N$ mit Wahrscheinlichkeit mindestens $1 - 2^{-t}$ erfolgreich teilen.

Die Grundidee des Algorithmus ist folgende. Wenn wir ein $a$ haben, für das die Ordnung $r$ von $a$ modulo $N$ gerade ist, dann ist $r/2$ eine ganze Zahl und wir können die Zahlen

$$a^{r/2} - 1\; (\textrm{mod}\; N) \quad \text{und} \quad a^{r/2} + 1\; (\textrm{mod}\; N)$$

betrachten. Mit der Formel $Z^2 - 1 = (Z+1)(Z-1)$ folgt:

$$\bigl(a^{r/2} - 1\bigr) \bigl(a^{r/2} + 1\bigr) = a^r - 1.$$

Wir wissen, dass $a^r \; (\textrm{mod}\; N) = 1$ per Definition der Ordnung — was gleichbedeutend damit ist, dass $N$ das Produkt $a^r - 1$ teilt. Das bedeutet, dass $N$ das Produkt

$$\bigl(a^{r/2} - 1\bigr) \bigl(a^{r/2} + 1\bigr)$$

teilt.

Damit das stimmt, müssen alle Primfaktoren von $N$ auch Primfaktoren von $a^{r/2} - 1$ oder $a^{r/2} + 1$ (oder beider) sein — und bei einer zufälligen Wahl von $a$ ist es unwahrscheinlich, dass alle Primfaktoren von $N$ nur einen der Terme teilen und keinen den anderen. Solange also einige der Primfaktoren von $N$ den ersten Term teilen und einige den zweiten, können wir durch Berechnung des ggT mit dem ersten Term einen nicht-trivialen Faktor von $N$ finden.