Pular para o conteúdo principal

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=HZZ+HX=−J∑(i,j)ZiZj+h∑iXiH = H_{ZZ} + H_X = - J \sum_{(i,j)} Z_i Z_j + h \sum_{i} X_i

em que J>0J > 0 é o acoplamento entre spins de primeiros vizinhos com i<ji < j e hh é 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,

exp⁡(−iHZZδt)=∏(i,j)exp⁡(iJδtZiZj)=∏(i,j)RZiZj(−2Jδt)exp⁡(−iHXδt)=∏iexp⁡(−ihδtXi)=∏iRXi(2hδt)\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 TT é discretizado em T/δtT / \delta t passos de Trotter, e RZiZj(θJ)\mathrm{R}_{Z_i Z_j}(\theta_J) e RXi(θh)\mathrm{R}_{X_i}(\theta_h) são as portas de rotação ZZZZ e XX, respectivamente.

Os experimentos foram realizados em um processador IBM Quantum® Eagle, um dispositivo de 127 qubits com conectividade heavy-hex, aplicando interações XX a todos os qubits e interações ZZZZ a todas as arestas do mapa de acoplamento. Vale notar que todas as interações ZZZZ 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 θJ=−π/2\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 ⟨Z62⟩\langle Z_{62} \rangle como observável. image.png 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 ⟨Z13⟩\langle Z_{13} \rangle como observável.

É só isso para a introdução. Vamos embarcar nessa aventura em escala de utilidade!

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-aer qiskit-ibm-runtime rustworkx
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 dois gates CXCX.

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")

Output of the previous code cell

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 CXCX.

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")

Output of the previous code cell

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"))

Output of the previous code cell

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 ZZZZ 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()

Output of the previous code cell

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"]}
)

Output of the previous code cell

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 devemos (e não devemos) solicitar 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 fake tem 28 arestas, deve haver 28*5 gates 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))

Output of the previous code cell

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 ⟨Z13⟩\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 de parede. 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()

Saída da célula de código anterior

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 ⟨Z62⟩\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()

Saída da célula de código anterior

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​