Experimento em escala de utilidade I

nota

Tamiya Onodera (5 de julho de 2024)

Baixe o PDF da aula original. Alguns trechos de código podem ter ficado desatualizados, pois são imagens estáticas.

O tempo de QPU aproximado para executar este experimento é de 45 segundos.

1. Introdução ao artigo de utilidade

Nesta lição, executamos um circuito em escala de utilidade que aparece no que chamamos informalmente de "o artigo de utilidade", publicado na Nature Vol 618, 15 de junho de 2023. O artigo trata da evolução temporal do modelo de Ising 2D com campo transversal. Em particular, os autores estudam a dinâmica temporal do Hamiltoniano,

$$H = H_{ZZ} + H_X = - J \sum_{(i,j)} Z_i Z_j + h \sum_{i} X_i$$

em que $J > 0$ é o acoplamento entre spins de primeiros vizinhos com $i < j$ e $h$ é o campo transversal global. Eles simulam a dinâmica de spins a partir de um estado inicial por meio da decomposição de Trotter de primeira ordem do operador de evolução temporal,

$$\begin{aligned} \exp(-i H_{ZZ} \delta t) &= \prod_{(i,j)} \exp (i J \delta t Z_i Z_j) = \prod_{(i,j)} \mathrm{R}_{Z_i Z_j} ( - 2 J \delta t) \\ \exp(-i H_X \delta t) &= \prod_{i} \exp (-i h \delta t X_i ) = \prod_{i} \mathrm{R}_{X_i} ( 2 h \delta t) \end{aligned}$$

em que o tempo de evolução $T$ é discretizado em $T / \delta t$ passos de Trotter, e $\mathrm{R}_{Z_i Z_j}(\theta_J)$ e $\mathrm{R}_{X_i}(\theta_h)$ são as portas de rotação $ZZ$ e $X$, respectivamente.

Os experimentos foram realizados em um processador IBM Quantum® Eagle, um dispositivo de 127 qubits com conectividade heavy-hex, aplicando interações $X$ a todos os qubits e interações $ZZ$ a todas as arestas do mapa de acoplamento. Vale notar que todas as interações $ZZ$ não podem ser aplicadas simultaneamente devido a "dependências de dados". Por isso, os autores colorem o mapa de acoplamento para agrupá-las em camadas. As interações de uma mesma camada recebem a mesma cor e podem ser aplicadas em paralelo.

Além disso, por simplicidade experimental, eles se concentraram no caso $\theta_J=-\pi /2$.

A contribuição inovadora do artigo está no fato de que os autores construíram circuitos quânticos em uma escala além da simulação por vetor de estados, executaram-nos em computadores quânticos ruidosos e conseguiram extrair resultados confiáveis. Dessa forma, eles demonstraram a utilidade dos computadores quânticos ruidosos. Para isso, aplicaram extrapolação de zero ruído (ZNE) com amplificação probabilística de erros (PEA) para mitigar os erros dos dispositivos ruidosos.

A partir daí, passamos a chamar esses experimentos e circuitos de "escala de utilidade".

1.1 Seu objetivo

Seu objetivo nesta lição é construir um circuito em escala de utilidade e executá-lo em um processador Eagle. Está fora do escopo deste notebook extrair resultados confiáveis, em parte porque a PEA é um recurso experimental do Qiskit no momento em que este texto foi escrito, e em parte porque aplicar ZNE com PEA exige um tempo considerável.

Concretamente, você deve construir e executar o circuito correspondente à Figura 4b do artigo, e plotar os pontos "não mitigados" por conta própria. Como você verá, trata-se de um circuito de 127 qubits $\times$ 60 camadas (20 passos de Trotter) com $\langle Z_{62} \rangle$ como observável. Parece complicado?   Não se preocupe. As três últimas lições deste curso servem como degraus. Para começar, vamos demonstrar um experimento em escala menor: construir e executar em um dispositivo fake um circuito de 27 qubits $\times$ 6 camadas (2 passos de Trotter) com $\langle Z_{13} \rangle$ como observável.

Isso encerra a introdução. Vamos embarcar nessa aventura em escala de utilidade!

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_aer": "qiskit-aer", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "rustworkx": "rustworkx"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import qiskit

qiskit.__version__

'2.0.2'

#!pip install qiskit_ibm_runtime
#!pip install qiskit_aer

import matplotlib.pyplot as plt
import numpy as np
import rustworkx as rx

from qiskit import QuantumCircuit, transpile
from qiskit.circuit import Parameter
from qiskit.circuit.library import YGate
from qiskit.quantum_info import SparsePauliOp
from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    fake_provider,
    EstimatorV2 as Estimator,
)
from qiskit_aer import AerSimulator

service = QiskitRuntimeService()

2. Preparação

2.1 Construindo RZZ(-$\pi$ / 2)

Primeiro, observe que a porta RZZ em geral requer duas portas $CX$.

from qiskit.circuit.library import RZZGate

θ_h = Parameter("$\\theta_h$")
qc1 = QuantumCircuit(2)
qc1.append(RZZGate(θ_h), [0, 1])
qc1.decompose(reps=1).draw("mpl")

Como mencionado acima, neste experimento nos concentramos na porta RZZ com um ângulo específico, -$\pi$ / 2. Conforme demonstrado no artigo, ela pode ser implementada com apenas um gate $CX$.

qc2 = QuantumCircuit(2)

qc2.sdg([0, 1])
qc2.append(YGate().power(1 / 2), [1])
qc2.cx(0, 1)
qc2.append(YGate().power(1 / 2).adjoint(), [1])

qc2.draw("mpl")

Definimos um gate a partir desse circuito para uso futuro.

rzz = qc2.to_gate(label="RZZ")

Vamos fazer um uso rápido do rzz recém-definido.

qc3 = QuantumCircuit(3)
qc3.append(rzz, [0, 1])
qc3.append(rzz, [0, 2])
display(qc3.draw("mpl"))
# display(qc.decompose(reps=1).draw("mpl"))

Antes de continuar, vamos verificar a equivalência lógica entre qc1 (a porta RZZ) para -pi/2 e o nosso rzz (ou qc2) recém-definido:

from qiskit.quantum_info import Operator

op1 = Operator(qc1.assign_parameters([-np.pi / 2]))
op2 = Operator(qc2)

op1.equiv(op2)

True

2.2 Colorindo o mapa de acoplamento

Vamos estudar como colorimos o mapa de acoplamento de um backend. Isso é necessário para agrupar as interações $ZZ$ em camadas.

Para começar, vamos visualizar o mapa de acoplamento de um backend. Note que os mapas de acoplamento de todos os dispositivos IBM Quantum atuais seguem a topologia heavy-hexagonal.

backend = service.least_busy(operational=True, simulator=False)

backend.coupling_map.draw()

Para colorir um mapa de acoplamento, usamos o rustworkx, um pacote Python para trabalhar com grafos e redes complexas. Ele oferece vários algoritmos de coloração, todos heurísticos e, portanto, sem garantia de encontrar a coloração mínima.

Dito isso, como os grafos heavy-hex são bipartidos, utilizamos graph_bipartite_edge_color, que deve encontrar a coloração mínima para esses grafos.

def color_coupling_map(backend):
    graph = backend.coupling_map.graph
    undirected_graph = graph.to_undirected(multigraph=False)
    edge_color_map = rx.graph_bipartite_edge_color(undirected_graph)
    if edge_color_map is None:
        edge_color_map = rx.graph_greedy_edge_color(undirected_graph)
    # build a map from color to a list of edges
    edge_index_map = undirected_graph.edge_index_map()
    color_edges_map = {color: [] for color in edge_color_map.values()}
    for edge_index, color in edge_color_map.items():
        color_edges_map[color].append(
            (edge_index_map[edge_index][0], edge_index_map[edge_index][1])
        )
    return edge_color_map, color_edges_map

Grafos heavy-hexagonais devem ser coloridos com três cores. Vamos verificar isso para o mapa de acoplamento acima.

edge_color_map, color_edges_map = color_coupling_map(backend)
print(
    f"{backend.name}, {backend.num_qubits}-qubit device, {len(color_edges_map.keys())} colors assigned."
)

ibm_strasbourg, 127-qubit device, 3 colors assigned.

Exatamente!

Por diversão, vamos pintar o mapa de acoplamento com as cores obtidas, usando o recurso de visualização do rustworkx.

color_str_map = {0: "green", 1: "red", 2: "blue"}

undirected_graph = backend.coupling_map.graph.to_undirected(multigraph=False)
for i in undirected_graph.edge_indices():
    undirected_graph.get_edge_data_by_index(i)["color"] = color_str_map[
        edge_color_map[i]
    ]

rx.visualization.graphviz_draw(
    undirected_graph, method="neato", edge_attr_fn=lambda edge: {"color": edge["color"]}
)

3. Resolvendo a evolução temporal trotterizada de um modelo de Ising 2D

Vamos definir uma rotina para construir o circuito do artigo de utilidade para a evolução temporal de um modelo de Ising 2D. A rotina recebe três parâmetros: um backend, um inteiro indicando o número de passos de Trotter e um booleano que controla a inserção de barreiras.

def get_utility_circuit(backend, num_steps: int, barrier: bool = False):
    num_qubits = backend.num_qubits
    _, color_edges_map = color_coupling_map(backend)
    θ_h = Parameter("$\\theta_h$")
    qc = QuantumCircuit(num_qubits)

    for i in range(num_steps):
        qc.rx(θ_h, range(num_qubits))

        for _, edge_list in color_edges_map.items():
            for edge in edge_list:
                qc.append(rzz, edge)

        if barrier:
            qc.barrier()
    return qc

Note que já realizamos manualmente o mapeamento e o roteamento de qubits para o circuito construído. Portanto, ao transpilar o circuito mais adiante, não pedimos (e não devemos pedir) ao transpilador que faça o mapeamento e o roteamento de qubits. Como você verá em breve, invocamos o transpilador com nível de otimização 1 e o método de layout "trivial".

Em seguida, definimos uma rotina simples para obter informações sobre o circuito construído, para uma verificação rápida.

def get_circuit_info(qc: QuantumCircuit, reps: int = 0):
    qc0 = qc.decompose(reps=reps)
    return (
        f"{qc0.num_qubits} qubits × {qc0.depth(lambda x: x.operation.num_qubits == 2)} layers ({qc0.depth()}-depth)"
        + ",  "
        + f"""Gate breakdown: {", ".join([f"{k.upper()} {v}" for k, v in qc0.count_ops().items()])}"""
    )

Vamos exercitar essas rotinas. Você deve ver um circuito de 27 qubits $\times$ 15 camadas (5 passos de Trotter). Como o dispositivo falso tem 28 arestas, deve haver 28*5 portas de emaranhamento.

backend = fake_provider.FakeTorontoV2()
num_steps = 5
qc = get_utility_circuit(backend, num_steps, True)

display(qc.draw(output="mpl", fold=-1))
print(get_circuit_info(qc, reps=0))
print(get_circuit_info(qc, reps=1))

27 qubits × 15 layers (20-depth),  Gate breakdown: CIRCUIT-165 140, RX 135, BARRIER 5
27 qubits × 15 layers (60-depth),  Gate breakdown: SDG 280, UNITARY 280, CX 140, R 135, BARRIER 5

4. Resolva a versão de 27 qubits do problema

Agora demonstramos uma versão em escala menor do experimento de utilidade. Construímos um circuito de 27 qubits $\times$ 6 camadas (2 passos de Trotter) com $\langle Z_{13} \rangle$ como observável, e o executamos tanto no AerSimulator quanto em um dispositivo falso.

Claro, seguimos nosso fluxo de trabalho de quatro etapas, o "padrão Qiskit", que consiste em Mapear, Otimizar, Executar e Pós-processar. De forma mais concreta:

• Mapear entradas clássicas para uma computação quântica.
• Otimizar circuitos para computação quântica.
• Executar circuitos usando primitivas.
• Pós-processar e retornar resultados em formato clássico.

A seguir, temos a etapa de Mapeamento para criar um circuito para um experimento em escala menor. Depois, temos um conjunto de Otimização e Execução para o AerSimulator e outro para um dispositivo falso. Por fim, temos a etapa de Pós-processamento para plotar os resultados.

4.1 Etapa 1: Mapeamento

backend = fake_provider.FakeTorontoV2()  # a 27 qubit fake device.
num_steps = 2
qc = get_utility_circuit(backend, num_steps)
obs = SparsePauliOp.from_sparse_list(
    [("Z", [13], 1)], num_qubits=backend.num_qubits
)  # Falcon
angles = [
    0,
    0.1,
    0.2,
    0.3,
    0.4,
    0.5,
    0.6,
    0.7,
    0.8,
    1.0,
    np.pi / 2,
]  # We try 11 angles for theta_h.

4.2 Etapas 2 e 3: Otimização e Execução (Simulador)

backend_sim = AerSimulator()
transpiled_qc_sim = transpile(
    qc, backend_sim, optimization_level=1, layout_method="trivial"
)
transpiled_obs_sim = obs.apply_layout(layout=transpiled_qc_sim.layout)

print(get_circuit_info(qc, reps=1))
print(get_circuit_info(transpiled_qc_sim, reps=1))

27 qubits × 6 layers (23-depth),  Gate breakdown: SDG 112, UNITARY 112, CX 56, R 54
27 qubits × 6 layers (16-depth),  Gate breakdown: U3 80, CX 56, R 54, U1 32, U 28

Um usuário executou a próxima célula em um MacBook Pro com processador Intel Core i7 quad-core de 2,3 GHz equipado com 32 GB de RAM 3LPDDR4X, rodando macOS 14.5. Levou 161 ms de tempo real (wall time). Cada laptop pode apresentar resultados ligeiramente diferentes.

%%time
params = [[p] for p in angles]
estimator = Estimator(mode=backend_sim)
pub = (transpiled_qc_sim, transpiled_obs_sim, params)
result_sim = estimator.run([pub]).result()

CPU times: user 231 ms, sys: 186 ms, total: 417 ms
Wall time: 111 ms

4.3 Etapas 2 e 3: Otimização e Execução (dispositivo falso)

backend_fake = fake_provider.FakeTorontoV2()
transpiled_qc_fake = transpile(
    qc, backend_fake, optimization_level=1, layout_method="trivial"
)
transpiled_obs_fake = obs.apply_layout(layout=transpiled_qc_fake.layout)

print(get_circuit_info(qc, reps=1))
print(get_circuit_info(transpiled_qc_fake, reps=1))

27 qubits × 6 layers (23-depth),  Gate breakdown: SDG 112, UNITARY 112, CX 56, R 54
27 qubits × 6 layers (49-depth),  Gate breakdown: SDG 324, U1 274, H 162, CX 56, U3 14

Quando o mesmo usuário executou a próxima célula no mesmo ambiente descrito acima, levou 2 min 19 s de tempo de parede. Executar um circuito em um dispositivo falso invoca uma simulação com ruído, o que demora muito mais do que uma simulação exata. Recomendamos que você não execute um circuito maior (como 27 qubits $\times$ 9 camadas com 3 passos de Trotter) em um dispositivo falso.

%%time
params = [[p] for p in angles]
estimator = Estimator(mode=backend_fake)
pub = (transpiled_qc_fake, transpiled_obs_fake, params)
result_fake = estimator.run([pub]).result()

CPU times: user 4min 42s, sys: 9.35 s, total: 4min 51s
Wall time: 38.3 s

4.4 Etapa 4: Pós-processamento

Plotamos os resultados das simulações exata e com ruído. Você verá os efeitos severos do ruído no FakeToronto.

plt.plot(angles, result_fake[0].data.evs, "o", label="Fake Device")
plt.plot(angles, result_sim[0].data.evs, "o", label="AerSimulator")
plt.xlabel("$\\mathrm{R_x}$ angle $\\theta_h$")
plt.title("$\\langle Z_{13} \\rangle$")
plt.legend()
plt.show()

5. Resolva a versão de 127 qubits do problema

Seu objetivo é executar o experimento em escala de utilidade mencionado no início. Você criará e executará um circuito de 127 qubits e 60 camadas (20 passos de Trotter) com $\langle Z_{62} \rangle$ como observável. Recomendamos que você tente fazer isso por conta própria, aproveitando o código da versão de 27 qubits quando apropriado. Mas a solução também é fornecida aqui.

Solução:

5.1 Etapa 1: Mapeamento

# backend_map = service.backend("ibm_brisbane")
backend_map = service.least_busy(operational=True, simulator=False)

num_steps = 20
qc = get_utility_circuit(backend_map, num_steps)
obs = SparsePauliOp.from_sparse_list(
    [("Z", [62], 1)], num_qubits=backend_map.num_qubits
)  # Eagle
angles = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0, np.pi / 2]

5.2 Etapas 2 e 3: Otimização e Execução

Observamos que o mapa de acoplamento do processador Eagle possui 144 arestas.

# backend = service.backend("ibm_brisbane")
backend = backend_map

transpiled_qc = transpile(qc, backend, optimization_level=1, layout_method="trivial")
transpiled_obs = obs.apply_layout(layout=transpiled_qc.layout)

print(get_circuit_info(qc, reps=1))
print(get_circuit_info(transpiled_qc))

156 qubits × 60 layers (221-depth),  Gate breakdown: SDG 7040, UNITARY 7040, CX 3520, R 3120
156 qubits × 60 layers (201-depth),  Gate breakdown: RZ 11933, SX 6240, CZ 3520

params = [[p] for p in angles]
estimator = Estimator(mode=backend)
pub = (transpiled_qc, transpiled_obs, params)
job = estimator.run([pub])

job_id = job.job_id()
print(f"job id={job_id}")

job id=d1479n6qf56g0081sxa0

5.3 Pós-processamento

Fornecemos os valores para os pontos "mitigados" da Figura 4b do artigo de utilidade. Plote-os junto com seus resultados.

result_paper = [
    1.0171,
    1.0044,
    0.9563,
    0.9602,
    0.8394,
    0.8120,
    0.5466,
    0.4556,
    0.1953,
    0.0141,
    0.0117,
]

# REPLACE WITH YOUR OWN JOB ID
job = service.job(job_id)

plt.plot(angles, job.result()[0].data.evs, "o", label=f"{job.backend().name}")
plt.plot(angles, result_paper, "o", label="Utility Paper")
plt.xlabel("$\\mathrm{R_x}$ angle $\\theta_h$")
plt.title("$\\langle Z_{62} \\rangle$")
plt.legend()
plt.show()

Seus resultados são semelhantes ao "não mitigado" da Figura 4b?   Eles podem ser bem diferentes, dependendo do dispositivo e de sua condição no momento do experimento. Não se preocupe com os resultados em si. O que vamos verificar é se você fez a codificação corretamente. Se sim, parabéns — você chegou à linha de largada da era de utilidade.

Assim como no artigo de Utilidade, cientistas de todo o mundo demonstraram uma enorme criatividade para extrair resultados significativos mesmo na presença de ruído. O objetivo final desse esforço coletivo é a vantagem quântica: um estado em que computadores quânticos conseguem resolver alguns problemas relevantes para a indústria de forma mais rápida, com maior fidelidade ou a um custo menor do que os computadores clássicos. É provável que isso não seja um evento único, mas sim uma era em que a reprodução clássica de resultados quânticos se torna progressivamente mais lenta, até que, em algum momento, essa vantagem temporal se torne criticamente importante. Uma coisa é certa sobre a vantagem quântica: só chegamos lá por meio de experimentos em escala de utilidade. Se este curso fizer com que você se junte a essa busca — cheia de desafios e diversão — ficaremos muito felizes.

Referência

• Kim, Y., Eddins, A., Anand, S. et al. Evidence for the utility of quantum computing before fault tolerance. Nature 618, 500–505 (2023). https://doi.org/10.1038/s41586-023-06096-3