Jun 11, 2024 | 8 min read
Algoritmos de búsqueda#
\( \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} \)
Show code cell source
%run ../macro_tQ.py
import sys
sys.path.append('../')
import macro_tQ as tQ
Una aguja en un pajar#
Promesa: disponemos de un oráculo que implementa una función \(f: x\to f(x) = 0,1\) donde \(x \in S_n = 0,1,..., 2^{n-1}\), con la siguiente propiedad: existe un único valor \(x = w\) tal que
Reto: descubrir \(w\) con el mínimo número de consultas al oráculo.
En el peor de los casos, el proceso es análogo al de buscar una aguja en un pajar, o el nombre de un cliente de la guía telefónica asociado a un número concreto. La función \(f(x)\) es la función comparar con un patrón solución que da 1 sólo cuando se acierta. También puede tratarse de una función que da 1 cuando \(w\) es la solución de algún problema algebráico de fácil verificación (por ejemplo un sudoku).
Si \(N = 2^N\) es la cardinalidad de la base de datos, el algoritmos de Grover permite efectuar esta tarea en un número de pasos que crece como \({\cal O}(\sqrt{N})\).
La ganancia por tanto, no es exponencial, como en el caso de la transformada de Fourier cuántica. De hecho, no es fácil saber a priori si un algoritmo va a tener una ganancia exponencial o no.
El plano de Grover#
Nuestra base de datos \(x\in S_n\) se convierte en la base computacional de \(n\) cúbits, \(\{\ket{x}\}\). La aguja en el pajar que queremos encontrar será, ahora, un cierto vector \(\ket{w}\), tal que
Aquí \(U_f\) será el oráculo booleano asociado a la función de búsqueda
Nota
La implementación de esta función es trivial porque sólo contiene un min-term \(f(x=w)=1\). Por supuesto, este es desconocido para el usuario del oráculo. También podemos implementar esta función como un oráculo de fase, preparando la ancilla en el estado \(\ket{-}\)
En resumen cada reflexión \(U_f\) es equivalente a una consulta al oráculo.
Como de costumbre, la estrategia cuántica comienza introduciendo una superposición uniforme de elementos de la base
Definamos el vector unitario compuesto por la superposición uniforme de los elementos que no son solución
Claramente \(\braket{\omega}{s'} = 0\) son ortogonales. Ahora podemos reescribir \(\ket{s}\) en una forma más conveniente
Si ahora definimos el vector unitario $\( \ket{s'} = \frac{1}{\sqrt{N-1}}\sum_{x \neq w} \ket{x} \)$
podemos escribir el resultado en la forma siguiente
Eso quiere decir que:
podemos pensar en un subespacio de dimensión 2 (un plano) formado por los dos vectores \(\{\ket{w},\ket{s'}\}\)
en ese plano, nuestro estado \(\ket{s}\) es un vector
caracterizado por el ángulo \(\sin\theta = \frac{1}{\sqrt{N}}\), como se puede observar en la figura siguiente
Notar
el ángulo disminuye con el tamaño \(N\) de la base de datos. Este ángulo es la amplitud de probabilidad de encontrar la solución \(w\) después de una medida de \(\ket{s}\).
Para \(N\to \infty\) tenemos que \(\theta \sim 1/\sqrt{N}\)
Rotación de Grover#
La estrategia se basa en iterar una serie de operaciones que tengan como efecto una amplificación de la probabilidad de obtener el estado \(\ket{w}\).
Geométricamente, cada iteración de Grover es una composición de dos reflectores \(G = R_{s_\perp} R_w\), donde
el marcador, \(R_w\), produce una reflexión paralela a \(\ket{w}\).
Es decir, invierte la componente asociada a \(\ket{w}\) y mantiene intactas todas las demás. Un operador así tiene la siguiente forma
Efectivamente, operando se llega inmediatamente a
el difusor, \(R_{s_\perp}\) es también un reflector.
Esta vez, \(R_{s_\perp}\) produce una reflexión a lo largo de un vector \(\ket{s_\perp}\), perpendicular a \(\ket{s}\). Como hemos visto, el operador es
En las ilustraciones hemos representado el efecto de una primera iteración de Grover \(G = R_{s_\perp} R_w\) actuando sobre el vector \(\ket{s}\).
Observamos que:
el vector resultante sufre una rotación de ángulo \(2\theta\) en el plano \(\{\ket{w},\ket{s'}\}\).
este hecho no depende de cómo sea \(\ket{s}\), por tanto: es cierto para cada iteración \(G\) de Grover.
la amplitud del estado \(\ket{w}\) aumenta a expensas de las amplitudes de los demás estados.
El circuito de Grover#
La tarea ahora es crear sendos circuitos que implementen las reflexiones \(R_w\) y \(R_s\). Comenzaremos repasando la prescripción general para programar estos operadores.
Lema:
Sea \(\ket{\psi}\) un vector preparable en la forma \(\ket{\psi} = V\ket{0}\). El reflector paralelo a \(\ket{\psi}\) es
De modo que sólo necesitamos un circuito para el reflector \(~R_0 = I- 2\ket{0}\bra{0}~\) a lo largo del estado \(\ket{0}\).
Lema:
El siguiente circuito implementa la acción del reflector
donde la cadena \(XZX\) puede situarse en cualquier cúbit, no necesariamente el último.
Es evidente que cualquier estado que no sea \(\ket{00...0}\) pasará intacto por el circuito. El estado \(\ket{00...0}\) activa la puera \(Z\ket{1} = -\ket{1}\) en el último cúbit. Ese signo es el que estamos buscando.
from qiskit.circuit.library import ZGate
from qiskit import QuantumRegister, ClassicalRegister,QuantumCircuit
def R0(n):
qr = QuantumRegister(n)
qc = QuantumCircuit(qr)
qc.x(n-1)
qc.append(ZGate().control(n-1,ctrl_state= 0),qr)
qc.x(n-1)
qcR0 = qc.to_gate()
qcR0.name = 'R0'
return qcR0
R0(4).definition.draw('mpl')
/opt/anaconda3/envs/TalentQ/lib/python3.11/site-packages/qiskit/visualization/circuit/matplotlib.py:266: FutureWarning: The default matplotlib drawer scheme will be changed to "iqp" in a following release. To silence this warning, specify the current default explicitly as style="clifford", or the new default as style="iqp".
self._style, def_font_ratio = load_style(self._style)
Codificando el marcador \(O = R_w\)#
Comenzando por \(\ket{\psi} = \ket{w}\), el reflector \(R_w\) es un operador que depende del vector desconocido \(\ket{w}\). Su acción es la reflexión
Pero esta es, precisamente la acción del oráculo asociado a \(f(x)\), implementado como oráculo de fase. Por tanto no estámos usando información prohibida. La primera reflexión no es más que la invocación de oráculo.
def Rw(n,w):
qr = QuantumRegister(n)
qc = QuantumCircuit(qr)
wstring = format(int(w), 'b').zfill(n)
for i,wi in enumerate(reversed(wstring)):
if wi == '1':
qc.x(qr[i])
qc.append(R0(n),qr)
# qc.append(ZGate().control(n-1,ctrl_state= 0),qr) # el reflector R0
for i,wi in enumerate(reversed(wstring)):
if wi == '1':
qc.x(qr[i])
qcRw = qc.to_gate()
qcRw.name='Rw'
return qcRw
# Sobre un circuito de 4 cúbits, sea el estado w = 6 = 0110
Rw(4,6).definition.draw('mpl')
Codificando el difusor \(D = R_{s_\perp}\)#
Ahora el papel de \(\ket{\psi}\) lo jugará el vector perpendicular \(\ket{s_\perp}\). Pero ¡no sabemos preparar este estado!
Sin embargo es evidente que \( R_{s_\perp} R_{s} = -I\) porque el producto de reflexiones paralela y perpendicular a un vector acaba invirtiendo todas las componentes. Dicho de otra manera \(R_{s_\perp} = - R_s\). Y para \(\ket{s}\) si que tenemos una preparación \(\ket{s} = H^n\ket{0}\).
Por tanto, si elegimos codificar el reflector \(R_s\) en lugar de \(R_{s_\perp}\) sólo estaremos cometiendo un error en un signo. A priori los signos globales son inobservables. Por tanto esperamos que no habrá ninguna diferencia entre usar un reflector u otro.
def Rs(n):
qr = QuantumRegister(n)
qc = QuantumCircuit(qr)
qc.h(qr)
qc.append(R0(n),qr)
qc.h(qr)
qcRs = qc.to_gate()
qcRs.name='Rs'
return qcRs
Rs(4).definition.draw('mpl')
Con estos dos ingredientes ya podríamos fabricar el operador de Grover.
Sin embargo, como hemos dicho, no sabemos implementar \(R_s\). Con lo que sabemos hacer, sólo podemos implementar el opuesto del operador de Grover
def G(n,w):
qr = QuantumRegister(n)
qc = QuantumCircuit(qr)
qc.append(Rw(n,w),qr)
qc.append(Rs(n),qr)
qcG=qc.to_gate()
qcG.name='-G'
return qcG
G(4,6).definition.draw('mpl')
El circuito de Grover#
Finalmente podemos componer el circuito, que contiene la solución \(w\) como argumento de \(G\)
# el número de cúbits de la base de datos
n=4
# el número de iteraciones de Grover
niter = 5
# la solución oculta en decimal y en binario
w = 6
wstring = format(int(w), 'b').zfill(n)
print(wstring)
# el circuito de Grover
qr = QuantumRegister(n)
cr = ClassicalRegister(n)
qcGrover = QuantumCircuit(qr,cr)
qcGrover.h(qr)
qcGrover.barrier()
for i in range(niter):
qcGrover.append(G(n,w),qr)
qcGrover.barrier()
# medimos el resultado
qcGrover.measure(qr,cr)
qcGrover.draw('mpl')
0110
'ejecutamos el circuito'
from qiskit.primitives import Sampler
from qiskit.visualization import plot_histogram
nshots = 1000
dist = Sampler().run(qcGrover,shots=nshots).result().quasi_dists[0]
counts = tQ.dist2counts(dist,n,nshots)
plot_histogram(counts)
La rotación óptima#
Se sigue de la figura que ilustra la acción de \(G\) que cada iteración de produce una rotación de ángulo \(2\theta\) donde \(\theta = \arcsin 1/\sqrt{N}\) es el ángulo inicial
Esta rotación puede representarse en la forma de una matriz
La iteración del circuito de Grover conduce a una sucesión de vectores
que giran sin salirse del plano \((\ket{w},\ket{s'})\)
Esta rotación comienza aumentanto la amplitud del estado \(\ket{w}\) ya que \(\sin 3\theta > \sin \theta\). Éste es el efecto de amplificación de amplitud buscado
Habrá un valor de \(k\) para el cuál la iteración se aproxime mucho al valor \(\sin (2k+1)\theta \approx 1\).
Valores de \(k\) superiores conducen a una supresión de dicha amplitud al rotar el vector \(G^k\ket{s}\) alejándose del eje \(\ket{w}\).
Por para conocer el número óptimo de iteraciones, debemos solucionar la ecuación de máxima amplitud
Esta ecuación tendría fácil solución
si \(k\) fuese un número real.
Pero es un número entero \( \Rightarrow \) el siguiente teorema da la solución óptima
Teorema:
El valor entero de \(k\in{\mathbb Z}\) que maximiza la probabilidad \(P_k = \sin^2((2k+1)\theta)\) de medir la solución \(w\) es la parte entera
donde \(\theta = \arcsin \sqrt{\frac{1}{N}}\). \(~\)Se cumple en este caso que:
Ejercicio
Ejecuta el algoritmo de Grover variando el número niter de iteraciones de Grover y verifica que el valor de \(k\) óptimo es el enunciado en el teorema.
Cuando \(N\to \infty\), es decir, cuando nuestra base de datos es muy grande el vector inicial tiene un ángulo muy pequeño.
a orden más bajo en serie de potencias de \(1/\sqrt{N}\)
Teorema:
En el límite de \(N\to\infty\) el tiempo que necesita el algoritmo de Grover para aproximarse máximamente a la solución crece como
Esto supone un crecimiento \(\mathcal{O}(\sqrt{N})\) en lugar de \(\mathcal{O}(N)\) que es lo que esperaríamos clásicamente.
Varias agujas en un pajar#
Un usuario podría tener más de una línea telefónica contratada. En ese caso, varios números de teléfono servirán para dar con un sólo nombre. Supondremos ahora que la búsqueda en la base de datos admite \(M\) resultados satisfactorios distintos \(w_1,....,w_M\).
Ahora los \(M\) vectores de la base que son solución \(\{\ket{w_i} \}\) generan un subespacio \(\mathcal{B} \subset \mathcal{H}\) de dimension \(M\).
Podemos descomponer \(~\Hil = {\cal B} \oplus {\cal B}_\perp\) en suma directa de subespacios solución y su ortogonal. Un vector \(\ket{x}\notin {\cal B}~\) será \(~\ket{x} \in {\cal B}_\perp\). Definamos
El vector inicial se puede descomponer en suma de sus componentes a lo largo de \({\cal B}\) y de \({\cal B}_\perp\)
Con estas definiciones, las figuras en las que se representan las reflexiones son idénticas y todo ocurre en el plano de Grover generado por los vectores \(\{\ket{w},\ket{s'}\}\). La única diferencia está en el valor del ángulo inicial
Teorema:
El valor entero de \(k\in{\mathbb Z}\) que maximiza la probabilidad \(P_k = \sin^2 ((2k+1)\theta)\) de medir la solución \(w\) es la parte entera de \(\pi/4\theta\)
donde \(\theta = \arcsin \sqrt{\frac{M}{N}}.~\) En este caso, la probabilidad de medir el estado \(\ket{w}\) es
en el límite \(N\to \infty\) el coste computacional crece aun como la raíz cuadrada
Teorema:
En el límite de \(N\to\infty\) el tiempo que necesita el algoritmo de Grover para aproximarse máximamente a alguna de las soluciones \(w_i\), crece como
Vemos que el número de iteraciones que hay que realizar es menor. A cambio sólo obtendremos una solución \(w_i\) después de cada ejecución.
El circuito de Grover modificado#
El único operador que necesitamos modificar es \(R_w\)
def Rw(n,w_array): #w_array es una lista de enteros w_i en S_n
qr = QuantumRegister(n)
qc = QuantumCircuit(qr)
for w in w_array:
wstring = format(int(w), 'b').zfill(n)
for i,wi in enumerate(reversed(wstring)):
if wi == '1':
qc.x(qr[i])
qc.append(R0(n),qr)
for i,wi in enumerate(reversed(wstring)):
if wi == '1':
qc.x(qr[i])
qcRw = qc.to_gate()
qcRw.name='Rw'
return qcRw
Ahora podemos correr el mismo circuito con esta modificación. Primero debemos calcular el valor óptimo de iteraciones \(k_{op} = [\pi/(4 \arcsin \sqrt{M/N})]\)
# dimensión de la base de datos
n=4
# lista de soluciones ocultas
w_array = [1, 3, 6, 7, 15]
# valor de k óptimo
k_op = int(np.pi/(4*np.arcsin(np.sqrt(len(w_array)/2**n))))
print('k_op =',k_op)
k_op = 1
Ahora construimos el circuito de Grover con \(k_{op}\) iteraciones.
# pasamos las soluciones a cadenas binarias de n entradas
wstring = [format(int(w), 'b').zfill(n) for w in w_array]
print('w_i =',wstring)
qr = QuantumRegister(n)
cr = ClassicalRegister(n)
qcGrover = QuantumCircuit(qr,cr)
qcGrover.h(qr)
qcGrover.barrier()
for i in range(k_op):
qcGrover.append(G(n,w_array),qr)
qcGrover.barrier()
qcGrover.measure(qr,cr)
qcGrover.draw('mpl')
w_i = ['0001', '0011', '0110', '0111', '1111']
/opt/anaconda3/envs/TalentQ/lib/python3.11/site-packages/qiskit/visualization/circuit/matplotlib.py:266: FutureWarning: The default matplotlib drawer scheme will be changed to "iqp" in a following release. To silence this warning, specify the current default explicitly as style="clifford", or the new default as style="iqp".
self._style, def_font_ratio = load_style(self._style)
'ejecutamos el circuito'
from qiskit.primitives import Sampler
from qiskit.visualization import plot_histogram
nshots = 1000
dist = Sampler().run(qcGrover,shots=nshots).result().quasi_dists[0]
counts = tQ.dist2counts(dist,n,nshots)
plot_histogram(counts)
Extraemos los valores que han obtenido las máximas frecuencias de medida
majority_counts = {k: v for k, v in counts.items() if v>100 }
#print(projected_counts)
plot_histogram(majority_counts)
los convertimos a base decimal
wi_list = sorted([int(keys,2) for keys in majority_counts])
print('{Ns} ~ ', wi_list)
{Ns} ~ [1, 3, 6, 7, 15]
Contaje cuántico#
En el algoritmos de Grover es crucial saber cuántas iteracions \(k\) del operador \(G\) hay que realizar. En el caso de múltiples soluciones, esto depende de \(M\). Pero ¡\(M\) es desconocido de antemano! Por tanto, el algoritmo ve comprometida su eficiencia.
En esta sección vamos a explicar cómo deducir \(M\) como paso previo, invocando el algoritmo de Estimación de Fase Cuántica.
Ya hemos visto que, en la base \(\{|\omega\rangle, |s'\rangle\}\), podemos escribir la iteración de Grover como la matriz de rotación
donde
Conocer \(\theta\) nos permitiría tener acceso al valor de $\( M = N\sin\theta \)$
Los autovalores y autoestados de este operador son los siguientes
que están pidiendo a gritos el algoritmo QPE, de estimación cuántica de fase.
Ejercicio
prueba este resultado
Una sutileza crucial#
En su momento comentamos que el operador que codificamos en el circuito de Grover no es realmente \( R_{s_\perp}\), sino \(R_s = -R_{s_\perp}\). Por tanto, realmente estamos implementando la acción de \( R_s R_w = -R_{s_\perp}R_w = -G \).
Este signo global lo hemos despreciado en el algoritmo de Grover. Sin embargo, ahora, en la estimación de fase, ya no se trata de un signo global, al estar el operador controlado. Por tanto tenemos que tenerlo en cuenta. El operador \(-G = G e^{i\pi}\) tiene por autovalores y autovectores
Vamos a definir las cantidades \(\varphi_\pm\) relacionadas con los autovalores de \(-G\)
El algoritmo de estimación de fase nos debería dar una aproximación buena a \(\varphi_\pm\), y a partir de estos dos valores, a
Vamos a insertar \((-G)^{2^k}\) en un circuito de estimación de fase cuántica con \(t\) cúbits de contaje
A la salida de este circuito, el registro de phase, estará en una combinación de estados \(\ket{x_+}\) y \(\ket{x_-}\) donde \(x_\pm\) están distribuidas en torno a los estimadores de orden \(t\) para \( a_\pm = [\varphi_\pm]\)
El operador \(-G \in SU(2)\) puesto que su determinante es 1 Por tanto
El algoritmo de estimación de fase nos dará estimaciones para \(~~\varphi_+ = \varphi ~~\) , \(~~\varphi_- = 1 - \varphi~,~~\) ambos pertenecientes al intervalo \([0,1)\).
En ambos casos es inmediato verificar que se cumple
A partir de aquí obtenemos
Para generar el circuito de QPE, tenemos que programar la puerta \(-G^k\) controlada \(~\Rightarrow ~\cg{(-G)^{k}}\)
def c_G(n,w_array,power):
qrG = QuantumRegister(n)
qcG = QuantumCircuit(qrG)
for _ in range(power):
qcG.append(G(n,w_array),qrG)
c_mG=qcG.to_gate()
c_mG.name="(-G)^%i" % power
c_mG = c_mG.control()
return c_mG
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
t = 5 # la dimensión del registro de evaluación
n = 4
qrt=QuantumRegister(t,name="t")
qra=QuantumRegister(n,name="qr_U") # aquí 1 es la dimensión del espacio en el que opera U
crt=ClassicalRegister(t)
qpe = QuantumCircuit(qrt,qra,crt)
qpe.h(range(t+n))
qpe.barrier()
# aplicamos los operadores controlados c((-G)^(2^j))
for j in range(t):
qpe.append(c_G(n,w_array,2**j),[qrt[j]]+qra[:])
qpe.barrier()
# aplicamos la QFT^-1 que está definida en macro.py
qpe.append(TFC_adj(t),qrt)
qpe.measure(qrt,crt)
qpe.draw(output="mpl")
/opt/anaconda3/envs/TalentQ/lib/python3.11/site-packages/qiskit/visualization/circuit/matplotlib.py:266: FutureWarning: The default matplotlib drawer scheme will be changed to "iqp" in a following release. To silence this warning, specify the current default explicitly as style="clifford", or the new default as style="iqp".
self._style, def_font_ratio = load_style(self._style)
Ejecutamos el circuito de estimación de fase
'ejecutamos el circuito'
from qiskit.primitives import Sampler
from qiskit.visualization import plot_histogram
nshots = 1000
dist = Sampler().run(qpe,shots=nshots).result().quasi_dists[0]
counts = tQ.dist2counts(dist,n,nshots)
plot_histogram(counts)
Extraemos las salidas \(m_i\) que tienen máxima probabilidad
majority_counts = {k: v for k, v in counts.items() if v>300 }
#print(projected_counts)
plot_histogram(majority_counts)
Convertimos a base decimal y extraemos sucesivamente $\( x_i ~~\longrightarrow ~~ a_i = \left[\frac{x_i}{2^t}\right] ~\approx ~\varphi_i ~~\longrightarrow ~~ M_i = [N \cos^2(\pi \varphi_i)] \)$
# la lista de m_i
xi_list = sorted([int(keys,2) for keys in majority_counts])
print('{x_i} ~ ', xi_list)
# la lista de a_i
ai_list = [x/2**t for x in xi_list]
# la lista de M_i
Mi_list = [2**n*np.cos(np.pi*a)**2 for a in ai_list]
print('M_i=', Mi_list)
# redondeo de M_i al entero más próximo
rounded_Mi_list = [round(Mi) for Mi in Mi_list]
print(' M_i =', rounded_Mi_list)
{x_i} ~ [10, 22]
M_i= [4.938532541079283, 4.938532541079277]
M_i = [5, 5]
Así hemos obtenido el número de soluciones, \(M\), buscado,
print('Number of solutions =', round(np.mean(rounded_Mi_list)))
Number of solutions = 5