piday-code.ipynb

Kernel: Python 3

Estimación de pi ($\pi$) utilizando el Algoritmo de Estimación de Fase Cuántica

1. Resumen rápido del Algoritmo de Estimación de Fase Cuántica

La Estimación de Fase Cuántica (Quantum Phase Estimation, QPE) es un algoritmo cuántico que constituye la base de muchos algoritmos cuánticos más complejos. En esencia, QPE resuelve un problema bastante sencillo: dado un operador $U$ y un estado cuántico $\vert\psi\rangle$ que es un valor propio de $U$ con $U\vert\psi\rangle = \exp\left(2 \pi i \theta\right)\vert\psi\rangle$, ¿podemos obtener una estimación de $\theta$?

La respuesta es sí. El algoritmo QPE nos da $2^n\theta$, donde $n$ es el número de qubits que usamos para estimar la fase $\theta$.

2. Estimación de $\pi$

En esta demostración, elegimos $$U = p(\theta), \vert\psi\rangle = \vert1\rangle$$ donde $$p(\theta) = \begin{bmatrix} 1 & 0\ 0 & \exp(i\theta) \end{bmatrix}$$ es una de las compuertas cuánticas disponibles en Qiskit, y $$p(\theta)\vert1\rangle = \exp(i\theta)\vert1\rangle.$$

Al elegir la fase para que nuestra compuerta sea $\theta = 1$, podemos resolver para $\pi$ usando las siguientes dos relaciones:

1. A partir de la salida del algoritmo QPE, medimos una estimación para $2^n\theta$. Entonces, $\theta = \text{medido} / 2^n$

2. Apartir de la definición de la compuerta $p(\theta)$ anterior, sabemos que $2\pi\theta = 1 \Rightarrow \pi = 1 / 2\theta$

Combinando estas dos relaciones, $\pi = 1 / \left(2 \times (\text{(medido)}/2^n)\right)$.

Para obtener una comprensión detallada del algoritmo QPE, consulta el capítulo dedicado a este en el Libro de Texto de Qiskit que se encuentra en qiskit.org/textbook.

3. Hora de escribir código

Comenzamos importando las bibliotecas necesarias.

## importar las herramientas necesarias para nuestro trabajo
from IPython.display import clear_output
from  qiskit import *
from qiskit.visualization import plot_histogram
import numpy as np
import matplotlib.pyplot as plotter
from qiskit.tools.monitor import job_monitor
# Configuraciones de visualización
import seaborn as sns, operator
sns.set_style("dark")

pi = np.pi

La función qft_dagger calcula la Transformada Cuántica de Fourier inversa. Para obtener información detallada sobre este algoritmo, consulta el capítulo dedicado a este en el Libro de Texto de Qiskit.

## Código para la transformada cuántica de Fourier inversa
## adaptado del Libro de Texto de Qiskit en
## qiskit.org/textbook

def qft_dagger(circ_, n_qubits):
    """n-qubit QFTdagger los primeros n qubits en circ"""
    for qubit in range(int(n_qubits/2)):
        circ_.swap(qubit, n_qubits-qubit-1)
    for j in range(0,n_qubits):
        for m in range(j):
            circ_.cp(-np.pi/float(2**(j-m)), m, j)
        circ_.h(j)

La siguiente función, qpe_pre, prepara el estado inicial para la estimación. Ten en cuenta que el estado inicial se crea aplicando una compuerta Hadamard en todos los qubits excepto en el último y configurando el último qubit en $\vert 1\rangle$.

## Código para el estado inicial de la Quantum Phase Estimation
## adaptado del Libro de Texto de Qiskit en qiskit.org/textbook
## Ten en cuenta que el estado inicial se crea aplicando H
## en los primeros n_qubits y configurando el último qubit en |psi> = |1>

def qpe_pre(circ_, n_qubits):
    circ_.h(range(n_qubits))
    circ_.x(n_qubits)

    for x in reversed(range(n_qubits)):
        for _ in range(2**(n_qubits-1-x)):
            circ_.cp(1, n_qubits-1-x, n_qubits)

A continuación, escribimos una función rápida, run_job, para ejecutar un circuito cuántico y devolver los resultados.

## Ejecuta un job de Qiskit en hardware o simuladores

def run_job(circ, backend, shots=1000, optimization_level=0):
    t_circ = transpile(circ, backend, optimization_level=optimization_level)
    qobj = assemble(t_circ, shots=shots)
    job = backend.run(qobj)
    job_monitor(job)
    return job.result().get_counts()

Luego, carga tu cuenta para usar el simulador en la nube o dispositivos reales.

## Carga tu cuenta de IBMQ si deseas utilizar 
## el simulador en la nube o dispositivos cuánticos reales
my_provider = IBMQ.load_account()
simulator_cloud = my_provider.get_backend('ibmq_qasm_simulator')
device = my_provider.get_backend('ibmq_16_melbourne')

ibmqfactory.load_account:WARNING:2021-02-18 11:32:59,136: Credentials are already in use. The existing account in the session will be replaced.

simulator = Aer.get_backend('qasm_simulator')

Finalmente, juntamos todo en una función llamada get_pi_estimate que usa n_qubits para obtener una estimación para $\pi$.

## Función para estimar pi
## Resumen: utilizando la notación del Libro de Texto de Qiskit (qiskit.org/textbook),
## haz una estimación de fase cuántica con el operador de 'fase' U = p(theta) y |psi> = |1>
## tal que p(theta)|1> = exp(2 x pi x i x theta)|1>
## Al establecer theta = 1 radianes, podemos resolver para pi
## usando 2^n x 1 radianes = conteo medido con mayor frecuencia = 2 x pi

def get_pi_estimate(n_qubits):

    # crear el circuito
    circ = QuantumCircuit(n_qubits + 1, n_qubits)
    # crear el estado de entrada
    qpe_pre(circ, n_qubits)
    # aplicar una barrera
    circ.barrier()
    # aplicar la transformada inversa de fourier
    qft_dagger(circ, n_qubits)
    # aplicar una barrera
    circ.barrier()
    # medir todos menos los últimos qubits
    circ.measure(range(n_qubits), range(n_qubits))

    # ejecutar el job y obtener los resultados
    counts = run_job(circ, backend=simulator, shots=10000, optimization_level=0)
    # print(counts) 

    # obtener el conteo que ocurrió con mayor frecuencia
    max_counts_result = max(counts, key=counts.get)
    max_counts_result = int(max_counts_result, 2)
    
    # resolver para pi a partir de los recuentos medidos
    theta = max_counts_result/2**n_qubits
    return (1./(2*theta))

Ahora, ejecuta la función get_pi_estimate con diferentes números de qubits e imprime las estimaciones.

# estimar pi usando diferentes números de qubits
nqs = list(range(2,12+1))
pi_estimates = []
for nq in nqs:
    thisnq_pi_estimate = get_pi_estimate(nq)
    pi_estimates.append(thisnq_pi_estimate)
    print(f"{nq} qubits, pi ≈ {thisnq_pi_estimate}")

Job Status: job has successfully run
2 qubits, pi ≈ 2.0
Job Status: job has successfully run
3 qubits, pi ≈ 4.0
Job Status: job has successfully run
4 qubits, pi ≈ 2.6666666666666665
Job Status: job has successfully run
5 qubits, pi ≈ 3.2
Job Status: job has successfully run
6 qubits, pi ≈ 3.2
Job Status: job has successfully run
7 qubits, pi ≈ 3.2
Job Status: job has successfully run
8 qubits, pi ≈ 3.1219512195121952
Job Status: job has successfully run
9 qubits, pi ≈ 3.1604938271604937
Job Status: job has successfully run
10 qubits, pi ≈ 3.1411042944785277
Job Status: job has successfully run
11 qubits, pi ≈ 3.1411042944785277
Job Status: job has successfully run
12 qubits, pi ≈ 3.1411042944785277

Y grafica todos los resultados.

plotter.plot(nqs, [pi]*len(nqs), '--r')
plotter.plot(nqs, pi_estimates, '.-', markersize=12)
plotter.xlim([1.5, 12.5])
plotter.ylim([1.5, 4.5])
plotter.legend(['$\pi$', 'estimado de $\pi$'])
plotter.xlabel('Número de qubits', fontdict={'size':20})
plotter.ylabel('$\pi$ y estimado de $\pi$', fontdict={'size':20})
plotter.tick_params(axis='x', labelsize=12)
plotter.tick_params(axis='y', labelsize=12)
plotter.show()

import qiskit
qiskit.__qiskit_version__

{'qiskit-terra': '0.16.1',
 'qiskit-aer': '0.7.2',
 'qiskit-ignis': '0.5.1',
 'qiskit-ibmq-provider': '0.11.1',
 'qiskit-aqua': '0.8.1',
 'qiskit': '0.23.2'}