2. Implementación simple

\( \newcommand{\bra}[1]{\langle #1|} \) \( \newcommand{\ket}[1]{|#1\rangle} \) \( \newcommand{\braket}[2]{\langle #1|#2\rangle} \) \( \newcommand{\i}{{\color{blue} i}} \) \( \newcommand{\Hil}{{\cal H}} \) \( \newcommand{\cg}[1]{{\rm C}#1} \) \( \newcommand{\lp}{\left(} \) \( \newcommand{\rp}{\right)} \) \( \newcommand{\lc}{\left[} \) \( \newcommand{\rc}{\right]} \) \( \newcommand{\lch}{\left\{} \) \( \newcommand{\rch}{\right\}} \) \( \newcommand{\Lp}{\Bigl(} \) \( \newcommand{\Rp}{\Bigr)} \) \( \newcommand{\Lc}{\Bigl[} \) \( \newcommand{\Rc}{\Bigr]} \) \( \newcommand{\Lch}{\Bigl\{} \) \( \newcommand{\Rch}{\Bigr\}} \) \( \newcommand{\rqa}{\quad \Rightarrow \quad} \)

2. Implementación simple#

import matplotlib.pyplot as plt
import numpy as np
from qiskit import QuantumCircuit, Aer, transpile, assemble
from qiskit.visualization import plot_histogram
# import qiskit.tools.jupyter
# import math
from math import gcd
import pandas as pd
from numpy.random import randint
from fractions import Fraction
%matplotlib inline

# Import qibo libraries
from qibo.models import Circuit
from qibo import  gates

2.1 - Hallar del periodo de una función (Period Finding)
2.2 - Implementación (ad hoc) en Qiskit para \(N=15\)

El algoritmo de Shor se basa en el algoritmo de Estimación de Fase Cuántica (Quantum Phase Estimation) ([6] y [13]), que a su vez usa la Transformada de Fourier Cuántica (Quantum Fourier Transform) ([7] y [12]). El algoritmo de Shor lo que hace es convertir el problema de la factorización de un número en un problema de encontrar el periodo de una función, el cual puede ser implementado en un tiempo polinómico.

Para factorizar un número \(N\) básicamente el algoritmo de Shor lo que hace calcular el periodo de una función (periódica) de la forma

(3)#\[\begin{equation} f(x) = a^x \text{mod} N \end{equation}\]

donde \(a\) y \(N\) son enteros positivos mayores que 1, siendo además \(a\) < \(N\) y no teniendo factores comunes. La operación (\(z\) mod\(N\)) a lo que se refiere es a quedarnos con el resto de dividir el número que \(z\) por \(N\). El periodo \(r\) de esta función se calcula mediante el algoritmo de estimación de fase cuántica. Una vez se tiene el periodo \(r\), si este es par (sino hay que probar con otro valor de \(a\)) se pueden calcular los factores de \(N\) ya que existe una alta probabilidad de que el máximo común divisor de \(N\) y \(a^{\,r/2}-1\) o \(a^{\,r/2}+1\) sea un factor propio de \(N\).

En esta guía vamos a seguir la explicación del Notebook de Qiskit [14].

2.1 - Hallar del periodo de una función (Period Finding)#

2.1.1 - La función#

Como comentamos en la introducción, lo que queremos es hallar el periodo de la función

\[ f(x) = a^x \text{mod } N \]

donde \(a\) y \(N\) son enteros positivos mayores que 1, siendo además \(a\) < \(N\) y no teniendo factores comunes. La operación (\(z\) mod\(N\)) a lo que se refiere es a quedarnos con el resto de dividir el número que \(z\) por \(N\). A este tipo de funciones se las denomina exponenciales moduladas.

Denominaremos \(r\) al valor del periodo de la función \(f(x)\), es decir, \(r\) es el mínimo valor entero para que se cumple:

\[ f(x+r) = f(x) \]

En la siguiente figura vemos un ejemplo de este tipo de funciones con \(a = 3\) y \(N=35\). Vemos que para este caso el periodo es \(r = 12\). (Véase que las lineas puntuadas que unen las cruces son solo por estética).

N = 35
a = 3

# Calculate the plotting data
xvals = np.arange(35, dtype = float)
yvals = [np.mod(a**x, N) for x in xvals]

# Use matplotlib to display it nicely
fig, ax = plt.subplots()
ax.plot(xvals, yvals, linewidth=1, linestyle='dotted', marker='x')
ax.set(xlabel='$x$', ylabel='$%i^x$ mod $%i$' % (a, N),
       title="Example of Periodic Function in Shor's Algorithm")
try: # plot r on the graph
    r = yvals[1:].index(1) +1 
    plt.annotate('', xy=(0,1), xytext=(r,1), arrowprops=dict(arrowstyle='<->'))
    plt.annotate('$r=%i$' % r, xy=(r/3,1.5))
    print('Value from f(0) to f(r):')
    print(yvals[:r+1])
except ValueError:
    print('Could not find period, check a < N and have no common factors.')

Value from f(0) to f(r):
[1.0, 3.0, 9.0, 27.0, 11.0, 33.0, 29.0, 17.0, 16.0, 13.0, 4.0, 12.0, 1.0]

../../../_images/109ebf69c93c17ac2e3569e95e49beef807585eff151c3303d2f049c1b28a98e.png

2.1.2 - Solución: Estimación de fase de un operador U#

El algoritmo de Shor se basa en implementar el algoritmo de estimación de fases al operador unitario

\[ U |y \rangle \equiv | ay \, \text{mod} N \rangle \]

Al aplicar sucesivas veces el operador \(U\) sobre el estado \(|1 \rangle\) vamos obteniendo los valores de \(f(x)\) con \(x \in \mathbb{N}\), esto es,

\[ U^x |1 \rangle = | f(x) \rangle \]

Por ejemplo, para el caso que vimos en la gráfica anterior (\(a=3\) y \(N=35\)) tenemos

(4)#\[\begin{align} U^0 |1\rangle & = |1\rangle \\ U |1\rangle & = |3\rangle \\ U^2 |1\rangle & = |9\rangle \\ \vdots \\ U^{r-1} |1\rangle & = |12\rangle \\ U^r |1\rangle & = |1\rangle \end{align}\]

(Recordemos que dado un estado de \(n\) qubits \(\left| x \right\rangle\) tenemos que \(\left| x \right\rangle = \left| x_1 x_2 \dots x_n \right\rangle\) donde \(x_1\) es el bit más significativo.)

Como podemos ver, aplicar una vez más el operador \(U\) significa pasar de un número al siguiente de la lista periódica. Veámoslo explícitamente:

(5)#\[\begin{align} U(U^0 |1\rangle) & = U(|1\rangle) = |3\rangle \\ U(U |1\rangle) & = U(|3\rangle) = |9\rangle \\ U(U^2 |1\rangle) & = U(|9\rangle) = |27\rangle\\ \vdots \\ U(U^{r-1} |1\rangle) & = U(|12\rangle) =|1\rangle \\ U(U^r |1\rangle) & = U(|1\rangle) = |3\rangle \end{align}\]

Con esto se entiende fácilmente que la superposición equiprobable de todos los estados es un autoestado del operador \(U\) con autovalor 1:

\[ | u_0 \rangle = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} |a^k \, \text{mod} N \rangle, \qquad \text{donde} \quad U|u_0 \rangle = |u_0\rangle \]

Un autoestado de autovalor 1 no nos es muy interesante a la hora de aplicar el algoritmo de estimación fase. Otro conjunto de autoestados mucho más interesantes son aquellos de la forma:

\[ | u_s \rangle = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} e^{- \boldsymbol{2 \pi i}k \frac{s}{r}} |a^k \, \text{mod} N \rangle, \qquad \text{donde} \quad U|u_s \rangle = e^{\boldsymbol{2 \pi i} \frac{s}{r}}|u_s\rangle \]

donde \(0 \leq s \leq r-1\). Si ahora aplicamos el algoritmo de estimación de fase cuántica a uno de estos autoestados \(|u_s \rangle\), lo que obtendremos en el registro de conteo es \(|2^n s/r \rangle\). De aquí podemos extraer el valor de \(r\). Sin embargo, para preparar el estado \(|u_s \rangle\) tenemos que conocer \(r\), es decir, lo que queremos calcular.

Un solución elegante y fácil de implementar es darnos cuenta de que la suma de todos estos estados \(|u_s \rangle\) nos da el estado \(|1\rangle\), esto es,

\[ \frac{1}{r} \sum_{s=0}^{r-1} |u_s \rangle = \frac{1}{r} \left(|u_0 \rangle + |u_1 \rangle + \dots +|u_{r-1} \rangle \right) = | 1 \rangle \]

Si ahora aplicamos el algoritmo de estimación de fase cuántico (QPS) al estado |1> (un estado fácilmente implementable) obtenemos una superposición equiprobable de estados de la forma \(|2^n s/r \rangle\), es decir:

\[ |0 \rangle|1 \rangle \xrightarrow{QPS} \frac{1}{\sqrt{r}} \left( \left|2^t \frac{1}{r} \right\rangle + \left|2^t \frac{2}{r} \right\rangle + \dots + \left|2^t \frac{r-1}{r} \right\rangle \right) |1 \rangle \]

donde \(t\) es el número de qubits del registro de conteo. Usando el las fracciones continuas [17] podemos calcular \(r\) a partir de los cocientes \(s/r\).

En la siguiente imagen podemos ver el circuito (en el orden de Qiskit para los qubits) que implementa la estimación de fase cuántica (figura tomada de [14])

Nota

Denominamos \(n\) al número de qubits que necesitamos para codificar el número \(N\) (que queremos factorizar) en un registro cuántico. Para aplicar el algoritmo de Shor se suelen usar \(t=2n\) qubits en el registro de conteo.

Nota

Véase que este circuito sigue el convenio más extendido en la literatura, es decir, aquel en el que el bit más significativo es el de más arriba y el bit menos significativo el de más abajo. Qiskit sigue el convenio contrario. Esto es, dado un estado

\[ \left| x \right\rangle = \left| a_{n-1} a_{n-2} \dots a_1 a_0 \right\rangle \]

que podemos escribirlo en la base decimal como

\[ x = a_{n-1} 2^{n-1} + a_{n-2} 2^{n-2} + \dots + a_0 2^{0} \]

tenemos que al escribir los circuitos

En este caso, para escribir el circuito en el convenio de qiskit lo que habría que hacer es invertir los qubits que controlan las puertas \(U\), es decir, aplicar \(U^{2^{t-1}}\) controlada por el último qubit, \(U^{2^{t-2}}\) controlada el penúltimo qubit, … y \(U^{2^{0}}\) controlada por el primer qubit.

2. Implementación simple

Contents

2. Implementación simple#

2.1 - Hallar del periodo de una función (Period Finding)#

2.1.1 - La función#

2.1.2 - Solución: Estimación de fase de un operador U#

2.2 - Implementación (ad hoc) en Qiskit para \(N=15\)#

2.2.1 - Caso con muchos shots (ineficiente pero didáctico)#

2.2.2 - Caso shot a shot (óptimo)#