Simulação de Hamiltoniano Ising com chutes usando circuitos dinâmicos

Estimativa de uso: 7,5 minutos em um processador Heron r3. (NOTA: Esta é apenas uma estimativa. Seu tempo de execução pode variar.) Circuitos dinâmicos são circuitos com feedforward clássico - em outras palavras, são medições de meio de circuito seguidas por operações lógicas clássicas que determinam operações quânticas condicionadas à saída clássica. Neste tutorial, simulamos o modelo Ising com chutes em uma rede hexagonal de spins e usamos circuitos dinâmicos para realizar interações além da conectividade física do hardware.

O modelo Ising tem sido estudado extensivamente em várias áreas da física. Ele modela spins que sofrem interações Ising entre sítios da rede, bem como chutes do campo magnético local em cada sítio. A evolução temporal Trotterizada dos spins considerada neste tutorial, retirada de [1], é dada pelo seguinte unitário:

$$U(\theta)=\left(\prod_{\langle j, k\rangle} \exp \left(i \frac{\pi}{8} Z_j Z_k\right)\right)\left(\prod_j \exp \left(-i \frac{\theta}{2} X_j\right)\right)$$

Para investigar a dinâmica dos spins, estudamos a magnetização média dos spins em cada sítio como função dos passos de Trotter. Portanto, construímos o seguinte observável:

$$\langle O\rangle = \frac{1}{N} \sum_i \langle Z_i \rangle$$

Para realizar a interação ZZ entre sítios da rede, apresentamos uma solução usando o recurso de circuito dinâmico, levando a uma profundidade de dois qubits significativamente menor em comparação com o método de roteamento padrão com portas SWAP. Por outro lado, as operações de feedforward clássico em circuitos dinâmicos geralmente têm tempos de execução mais longos do que portas quânticas; portanto, circuitos dinâmicos têm limitações e compensações. Também apresentamos uma maneira de adicionar uma sequência de desacoplamento dinâmico em qubits ociosos durante a operação de feedforward clássico usando a duração stretch.

Requisitos

Antes de iniciar este tutorial, certifique-se de ter o seguinte instalado:

• Qiskit SDK v2.0 ou posterior com suporte de visualização
• Qiskit Runtime v0.37 ou posterior com suporte de visualização (pip install 'qiskit-ibm-runtime[visualization]')
• Biblioteca de grafos Rustworkx (pip install rustworkx)
• Qiskit Aer (pip install qiskit-aer)

Configuração

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-aer qiskit-ibm-runtime rustworkx

import numpy as np
from typing import List
import rustworkx as rx
import matplotlib.pyplot as plt
from rustworkx.visualization import mpl_draw
from qiskit.circuit import (
    Parameter,
    QuantumCircuit,
    QuantumRegister,
    ClassicalRegister,
)
from qiskit.transpiler import CouplingMap
from qiskit.quantum_info import SparsePauliOp
from qiskit.circuit.classical import expr
from qiskit.transpiler.preset_passmanagers import (
    generate_preset_pass_manager,
)
from qiskit.transpiler import PassManager
from qiskit.circuit.library import RZGate, XGate
from qiskit.transpiler.passes import (
    ALAPScheduleAnalysis,
    PadDynamicalDecoupling,
)

from qiskit.transpiler.basepasses import TransformationPass
from qiskit.circuit.measure import Measure
from qiskit.transpiler.passes.utils.remove_final_measurements import (
    calc_final_ops,
)
from qiskit.circuit import Instruction

from qiskit.visualization import plot_circuit_layout
from qiskit.circuit.tools import pi_check

from qiskit_aer import AerSimulator
from qiskit_aer.primitives import SamplerV2 as Aer_Sampler

from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    Batch,
    SamplerV2 as Sampler,
)
from qiskit_ibm_runtime.exceptions import QiskitBackendNotFoundError
from qiskit_ibm_runtime.visualization import (
    draw_circuit_schedule_timing,
)

Passo 1: Mapear entradas clássicas para um circuito quântico

Começamos definindo a rede a ser simulada. Escolhemos trabalhar com a rede em favo de mel (também chamada hexagonal), que é um grafo planar com nós de grau 3. Aqui, especificamos o tamanho da rede, os parâmetros de circuito relevantes de interesse na dinâmica Trotterizada. Simulamos a evolução temporal Trotterizada sob o modelo Ising sob três valores diferentes de $\theta$ do campo magnético local.

hex_rows = 3  # specify lattice size
hex_cols = 5
depths = range(9)  # specify Trotter steps
zz_angle = np.pi / 8  # parameter for ZZ interaction
max_angle = np.pi / 2  # max theta angle
points = 3  # number of theta parameters

θ = Parameter("θ")
params = np.linspace(0, max_angle, points)

def make_hex_lattice(hex_rows=1, hex_cols=1):
    """Define hexagon lattice."""
    hex_cmap = CouplingMap.from_hexagonal_lattice(
        hex_rows, hex_cols, bidirectional=False
    )
    data = list(hex_cmap.physical_qubits)
    graph = hex_cmap.graph.to_undirected(multigraph=False)
    edge_colors = rx.graph_misra_gries_edge_color(graph)
    layer_edges = {color: [] for color in edge_colors.values()}
    for edge_index, color in edge_colors.items():
        layer_edges[color].append(graph.edge_list()[edge_index])
    return data, layer_edges, hex_cmap, graph

Vamos começar com um pequeno exemplo de teste:

hex_rows_test = 1
hex_cols_test = 2

data_test, layer_edges_test, hex_cmap_test, graph_test = make_hex_lattice(
    hex_rows=hex_rows_test, hex_cols=hex_cols_test
)

# display a small example for illustration
node_colors_test = ["lightblue"] * len(graph_test.node_indices())
pos = rx.graph_spring_layout(
    graph_test,
    k=5 / np.sqrt(len(graph_test.nodes())),
    repulsive_exponent=1,
    num_iter=150,
)
mpl_draw(graph_test, node_color=node_colors_test, pos=pos)

Usaremos o exemplo pequeno para ilustração e simulação. Abaixo também construímos um exemplo grande para mostrar que o fluxo de trabalho pode ser estendido para tamanhos grandes.

data, layer_edges, hex_cmap, graph = make_hex_lattice(
    hex_rows=hex_rows, hex_cols=hex_cols
)
num_qubits = len(data)
print(f"num_qubits = {num_qubits}")

# display the honeycomb lattice to simulate
node_colors = ["lightblue"] * len(graph.node_indices())
pos = rx.graph_spring_layout(
    graph,
    k=5 / np.sqrt(num_qubits),
    repulsive_exponent=1,
    num_iter=150,
)
mpl_draw(graph, node_color=node_colors, pos=pos)
plt.show()

num_qubits = 46

Construir circuitos unitários

Com o tamanho do problema e os parâmetros especificados, estamos prontos para construir o circuito parametrizado que simula a evolução temporal Trotterizada de $U(\theta)$ com diferentes passos de Trotter, especificados pelo argumento depth. O circuito que construímos tem camadas alternadas de portas Rx($\theta$) e portas Rzz. As portas Rzz realizam as interações ZZ entre spins acoplados, que serão colocadas entre cada sítio da rede especificado pelo argumento layer_edges.

def gen_hex_unitary(
    num_qubits=6,
    zz_angle=np.pi / 8,
    layer_edges=[
        [(0, 1), (2, 3), (4, 5)],
        [(1, 2), (3, 4), (5, 0)],
    ],
    θ=Parameter("θ"),
    depth=1,
    measure=False,
    final_rot=True,
):
    """Build unitary circuit."""
    circuit = QuantumCircuit(num_qubits)
    # Build trotter layers
    for _ in range(depth):
        for i in range(num_qubits):
            circuit.rx(θ, i)
        circuit.barrier()
        for coloring in layer_edges.keys():
            for e in layer_edges[coloring]:
                circuit.rzz(zz_angle, e[0], e[1])
        circuit.barrier()
    # Optional final rotation, set True to be consistent with Ref. [1]
    if final_rot:
        for i in range(num_qubits):
            circuit.rx(θ, i)
    if measure:
        circuit.measure_all()

    return circuit

Visualizar o pequeno circuito de teste:

circ_unitary_test = gen_hex_unitary(
    num_qubits=len(data_test),
    layer_edges=layer_edges_test,
    θ=Parameter("θ"),
    depth=1,
    measure=True,
)
circ_unitary_test.draw(output="mpl", fold=-1)

Da mesma forma, construir os circuitos unitários do exemplo grande em diferentes passos de Trotter e o observável para estimar o valor esperado.

circuits_unitary = []
for depth in depths:
    circ = gen_hex_unitary(
        num_qubits=num_qubits,
        layer_edges=layer_edges,
        θ=Parameter("θ"),
        depth=depth,
        measure=True,
    )
    circuits_unitary.append(circ)

observables_unitary = SparsePauliOp.from_sparse_list(
    [("Z", [i], 1 / num_qubits) for i in range(num_qubits)],
    num_qubits=num_qubits,
)

Construir implementação de circuito dinâmico

Esta seção demonstra a principal implementação de circuito dinâmico para simular a mesma evolução temporal Trotterizada. Observe que a rede em favo de mel que queremos simular não corresponde à rede pesada dos qubits de hardware. Uma maneira direta de mapear o circuito para o hardware é introduzir uma série de operações SWAP para trazer qubits interagentes próximos uns dos outros, para realizar a interação ZZ. Aqui destacamos uma abordagem alternativa usando circuitos dinâmicos como solução, que ilustra que podemos usar a combinação de computação quântica e clássica em tempo real dentro de um circuito no Qiskit para realizar interações além do vizinho mais próximo.

Na implementação do circuito dinâmico, a interação ZZ é efetivamente implementada usando qubits ancilla, medição de meio de circuito e feedforward. Para entender isso, observe que as rotações ZZ aplicam um fator de fase $e^{i\theta}$ ao estado com base em sua paridade. Para dois qubits, os estados da base computacional são $|00\rangle$, $|01\rangle$, $|10\rangle$ e $|11\rangle$. A porta de rotação ZZ aplica um fator de fase aos estados $|01\rangle$ e $|10\rangle$ cuja paridade (o número de uns no estado) é ímpar e deixa os estados de paridade par inalterados. O seguinte descreve como podemos implementar efetivamente interações ZZ em dois qubits usando circuitos dinâmicos.

1. Calcular paridade em um qubit ancilla: em vez de aplicar ZZ diretamente a dois qubits, introduzimos um terceiro qubit, o qubit ancilla, para armazenar as informações de paridade dos dois qubits de dados. Emaranhamos o ancilla com cada qubit de dados usando portas CX do qubit de dados para o qubit ancilla.

2. Aplicar uma rotação Z de um único qubit ao qubit ancilla: isso ocorre porque o ancilla tem as informações de paridade dos dois qubits de dados, o que efetivamente implementa a rotação ZZ nos qubits de dados.

3. Medir o qubit ancilla na base X: este é o passo chave que colapsa o estado do qubit ancilla, e o resultado da medição nos diz o que aconteceu:

   • Medir 0: quando um resultado 0 é observado, na verdade aplicamos corretamente uma rotação $ZZ(\theta)$ aos nossos qubits de dados.

   • Medir 1: quando um resultado 1 é observado, aplicamos $ZZ(\theta + \pi)$ em vez disso.

4. Aplicar porta de correção ao medir 1: Se medimos 1, aplicamos portas Z aos qubits de dados para "corrigir" a fase extra $\pi$.

O circuito resultante é o seguinte:

Quando adotamos essa abordagem para simular uma rede em favo de mel, o circuito resultante se encaixa perfeitamente no hardware com uma rede heavy-hex: todos os qubits de dados residem nos sítios de grau 3 da rede, que forma uma rede hexagonal. Cada par de qubits de dados compartilha um qubit ancilla residindo em um sítio de grau 2. Abaixo, construímos a rede de qubits para a implementação do circuito dinâmico, introduzindo qubits ancilla (mostrados nos círculos roxos mais escuros).

def make_lattice(hex_rows=1, hex_cols=1):
    """Define heavy-hex lattice and corresponding lists of data and ancilla nodes."""
    hex_cmap = CouplingMap.from_hexagonal_lattice(
        hex_rows, hex_cols, bidirectional=False
    )
    data = list(hex_cmap.physical_qubits)

    heavyhex_cmap = CouplingMap()
    for d in data:
        heavyhex_cmap.add_physical_qubit(d)

    # make coupling map
    a = len(data)
    for edge in hex_cmap.get_edges():
        heavyhex_cmap.add_physical_qubit(a)
        heavyhex_cmap.add_edge(edge[0], a)
        heavyhex_cmap.add_edge(edge[1], a)
        a += 1
    ancilla = list(range(len(data), a))
    qubits = data + ancilla

    # color edges
    graph = heavyhex_cmap.graph.to_undirected(multigraph=False)
    edge_colors = rx.graph_misra_gries_edge_color(graph)
    layer_edges = {color: [] for color in edge_colors.values()}
    for edge_index, color in edge_colors.items():
        layer_edges[color].append(graph.edge_list()[edge_index])

    # construct observable
    obs_hex = SparsePauliOp.from_sparse_list(
        [("Z", [i], 1 / len(data)) for i in data],
        num_qubits=len(qubits),
    )

    return (data, qubits, ancilla, layer_edges, heavyhex_cmap, graph, obs_hex)

Visualizar a rede heavy-hex para qubits de dados e qubits ancilla em pequena escala:

(data, qubits, ancilla, layer_edges, heavyhex_cmap, graph, obs_hex) = (
    make_lattice(hex_rows=hex_rows, hex_cols=hex_cols)
)

print(f"number of data qubits = {len(data)}")
print(f"number of ancilla qubits = {len(ancilla)}")

node_colors = []
for node in graph.node_indices():
    if node in ancilla:
        node_colors.append("purple")
    else:
        node_colors.append("lightblue")

pos = rx.graph_spring_layout(
    graph,
    k=1 / np.sqrt(len(qubits)),
    repulsive_exponent=2,
    num_iter=200,
)

# Visualize the graph, blue circles are data qubits and purple circles are ancillas
mpl_draw(graph, node_color=node_colors, pos=pos)
plt.show()

number of data qubits = 46
number of ancilla qubits = 60

Abaixo, construímos o circuito dinâmico para a evolução temporal Trotterizada. As portas RZZ são substituídas pela implementação de circuito dinâmico usando os passos descritos acima.

def gen_hex_dynamic(
    depth=1,
    zz_angle=np.pi / 8,
    θ=Parameter("θ"),
    hex_rows=1,
    hex_cols=1,
    measure=False,
    add_dd=True,
):
    """Build dynamic circuits."""
    (data, qubits, ancilla, layer_edges, heavyhex_cmap, graph, obs_hex) = (
        make_lattice(hex_rows=hex_rows, hex_cols=hex_cols)
    )
    # Initialize circuit
    qr = QuantumRegister(len(qubits), "qr")
    cr = ClassicalRegister(len(ancilla), "cr")
    circuit = QuantumCircuit(qr, cr)

    for k in range(depth):
        # Single-qubit Rx layer
        for d in data:
            circuit.rx(θ, d)
        circuit.barrier()

        # CX gates from data qubits to ancilla qubits
        for same_color_edges in layer_edges.values():
            for e in same_color_edges:
                circuit.cx(e[0], e[1])
        circuit.barrier()

        # Apply Rz rotation on ancilla qubits and rotate into X basis
        for a in ancilla:
            circuit.rz(zz_angle, a)
            circuit.h(a)
        # Add barrier to align terminal measurement
        circuit.barrier()

        # Measure ancilla qubits
        for i, a in enumerate(ancilla):
            circuit.measure(a, i)
        d2ros = {}
        a2ro = {}
        # Retrieve ancilla measurement outcomes
        for a in ancilla:
            a2ro[a] = cr[ancilla.index(a)]

        # For each data qubit, retrieve measurement outcomes of neighboring ancilla qubits
        for d in data:
            ros = [a2ro[a] for a in heavyhex_cmap.neighbors(d)]
            d2ros[d] = ros

        # Build classical feedforward operations (optionally add DD on idling data qubits)
        for d in data:
            if add_dd:
                circuit = add_stretch_dd(circuit, d, f"data_{d}_depth_{k}")

            # # XOR the neighboring readouts of the data qubit; if True, apply Z to it
            ros = d2ros[d]
            parity = ros[0]
            for ro in ros[1:]:
                parity = expr.bit_xor(parity, ro)
            with circuit.if_test(expr.equal(parity, True)):
                circuit.z(d)

        # Reset the ancilla if its readout is 1
        for a in ancilla:
            with circuit.if_test(expr.equal(a2ro[a], True)):
                circuit.x(a)
        circuit.barrier()

    # Final single-qubit Rx layer to match the unitary circuits
    for d in data:
        circuit.rx(θ, d)

    if measure:
        circuit.measure_all()
    return circuit, obs_hex

def add_stretch_dd(qc, q, name):
    """Add XpXm DD sequence."""
    s = qc.add_stretch(name)
    qc.delay(s, q)
    qc.x(q)
    qc.delay(s, q)
    qc.delay(s, q)
    qc.rz(np.pi, q)
    qc.x(q)
    qc.rz(-np.pi, q)
    qc.delay(s, q)
    return qc

Desacoplamento dinâmico (DD) e suporte para duração stretch

Uma ressalva do uso da implementação de circuito dinâmico para realizar a interação ZZ é que a medição de meio de circuito e as operações de feedforward clássico normalmente levam mais tempo para executar do que portas quânticas. Para suprimir a decoerência de qubit durante o tempo ocioso para as operações clássicas acontecerem, adicionamos uma sequência de desacoplamento dinâmico (DD) após a operação de medição nos qubits ancilla, e antes da operação Z condicional no qubit de dados, antes da instrução if_test.

A sequência DD é adicionada pela função add_stretch_dd(), que usa as durações stretch para determinar os intervalos de tempo entre as portas DD. Uma duração stretch é uma maneira de especificar uma duração de tempo extensível para a operação delay de modo que a duração do atraso possa crescer para preencher o tempo ocioso do qubit. As variáveis de duração especificadas por stretch são resolvidas em tempo de compilação em durações desejadas que satisfazem uma certa restrição. Isso é muito útil quando o timing das sequências DD é essencial para alcançar um bom desempenho de supressão de erros. Para mais detalhes sobre o tipo stretch, consulte a documentação OpenQASM. Atualmente, o suporte para o tipo stretch no Qiskit Runtime é experimental. Para detalhes sobre suas restrições de uso, consulte a seção de limitações da documentação stretch.

Usando as funções definidas acima, construímos os circuitos de evolução temporal Trotterizada, com e sem DD, e os observáveis correspondentes. Começamos visualizando o circuito dinâmico de um pequeno exemplo:

hex_rows_test = 1
hex_cols_test = 1

(
    data_test,
    qubits_test,
    ancilla_test,
    layer_edges_test,
    heavyhex_cmap_test,
    graph_test,
    obs_hex_test,
) = make_lattice(hex_rows=hex_rows_test, hex_cols=hex_cols_test)

node_colors = []
for node in graph_test.node_indices():
    if node in ancilla_test:
        node_colors.append("purple")
    else:
        node_colors.append("lightblue")
pos = rx.graph_spring_layout(
    graph_test,
    k=5 / np.sqrt(len(qubits_test)),
    repulsive_exponent=2,
    num_iter=150,
)

# display a small example for illustration
node_colors_test = ["lightblue"] * len(graph_test.node_indices())
mpl_draw(graph_test, node_color=node_colors, pos=pos)

circuit_dynamic_test, obs_dynamic_test = gen_hex_dynamic(
    depth=1,
    θ=Parameter("θ"),
    hex_rows=hex_rows_test,
    hex_cols=hex_cols_test,
    measure=False,
    add_dd=False,
)
circuit_dynamic_test.draw("mpl", fold=-1)

circuit_dynamic_dd_test, _ = gen_hex_dynamic(
    depth=1,
    θ=Parameter("θ"),
    hex_rows=hex_rows_test,
    hex_cols=hex_cols_test,
    measure=False,
    add_dd=True,
)
circuit_dynamic_dd_test.draw("mpl", fold=-1)

Da mesma forma, construir os circuitos dinâmicos para o exemplo grande:

circuits_dynamic = []
circuits_dynamic_dd = []
observables_dynamic = []
for depth in depths:
    circuit, obs = gen_hex_dynamic(
        depth=depth,
        θ=Parameter("θ"),
        hex_rows=hex_rows,
        hex_cols=hex_cols,
        measure=True,
        add_dd=False,
    )
    circuits_dynamic.append(circuit)

    circuit_dd, _ = gen_hex_dynamic(
        depth=depth,
        θ=Parameter("θ"),
        hex_rows=hex_rows,
        hex_cols=hex_cols,
        measure=True,
        add_dd=True,
    )
    circuits_dynamic_dd.append(circuit_dd)
    observables_dynamic.append(obs)

Passo 2: Otimizar problema para execução em hardware

Agora estamos prontos para transpilar o circuito para o hardware. Iremos transpilar tanto a implementação padrão unitária quanto a implementação de circuito dinâmico para o hardware.

Para transpilar para hardware, primeiro instanciamos o backend. Se disponível, escolheremos um backend onde a instrução MidCircuitMeasure (measure_2) é suportada.

service = QiskitRuntimeService()
try:
    backend = service.least_busy(
        operational=True,
        simulator=False,
        use_fractional_gates=True,
        filters=lambda b: "measure_2" in b.supported_instructions,
    )
except QiskitBackendNotFoundError:
    backend = service.least_busy(
        operational=True,
        simulator=False,
        use_fractional_gates=True,
    )

Transpilação para circuitos dinâmicos

Primeiro, transpilamos os circuitos dinâmicos, com e sem adicionar a sequência DD. Para garantir que usamos o mesmo conjunto de qubits físicos em todos os circuitos para resultados mais consistentes, primeiro transpilamos o circuito uma vez, e então usamos seu layout para todos os circuitos subsequentes, especificado por initial_layout no gerenciador de passes. Em seguida, construímos os blocos unificados primitivos (PUBs) como entrada primitiva do Sampler.

pm_temp = generate_preset_pass_manager(
    optimization_level=3,
    backend=backend,
)
isa_temp = pm_temp.run(circuits_dynamic[-1])
dynamic_layout = isa_temp.layout.initial_index_layout(filter_ancillas=True)

pm = generate_preset_pass_manager(
    optimization_level=3, backend=backend, initial_layout=dynamic_layout
)

dynamic_isa_circuits = [pm.run(circ) for circ in circuits_dynamic]
dynamic_pubs = [(circ, params) for circ in dynamic_isa_circuits]

dynamic_isa_circuits_dd = [pm.run(circ) for circ in circuits_dynamic_dd]
dynamic_pubs_dd = [(circ, params) for circ in dynamic_isa_circuits_dd]

Podemos visualizar o layout de qubit do circuito transpilado abaixo. Os círculos pretos mostram os qubits de dados e os qubits ancilla usados na implementação de circuito dinâmico.

def _heron_coords_r2():
    cord_map = np.array(
        [
            [
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
                3,
                7,
                11,
                15,
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
                1,
                5,
                9,
                13,
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
                3,
                7,
                11,
                15,
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
                1,
                5,
                9,
                13,
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
                3,
                7,
                11,
                15,
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
                1,
                5,
                9,
                13,
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
                3,
                7,
                11,
                15,
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
            ],
            -1
            * np.array([j for i in range(15) for j in [i] * [16, 4][i % 2]]),
        ],
        dtype=int,
    )

    hcords = []
    ycords = cord_map[0]
    xcords = cord_map[1]
    for i in range(156):
        hcords.append([xcords[i] + 1, np.abs(ycords[i]) + 1])

    return hcords

plot_circuit_layout(
    dynamic_isa_circuits_dd[8],
    backend,
    qubit_coordinates=_heron_coords_r2(),
    view="virtual",
)

nota

Se você receber erros sobre neato não encontrado de plot_circuit_layout(), certifique-se de ter o pacote graphviz instalado e disponível em seu PATH. Se ele instalar em um local não padrão (por exemplo, usando homebrew no MacOS), pode ser necessário atualizar sua variável de ambiente PATH. Isso pode ser feito dentro deste notebook usando o seguinte:

import os
os.environ['PATH'] = f"path/to/neato{os.pathsep}{os.environ['PATH']}"

dynamic_isa_circuits[1].draw(fold=-1, output="mpl", idle_wires=False)

dynamic_isa_circuits_dd[1].draw(fold=-1, output="mpl", idle_wires=False)

Transpilar usando MidCircuitMeasure

MidCircuitMeasure é uma adição às operações de medição disponíveis, calibrada especificamente para realizar medições de meio de circuito. A instrução MidCircuitMeasure mapeia para a instrução measure_2 suportada pelos backends. Note que measure_2 não é suportado em todos os backends. Você pode usar service.backends(filters=lambda b: "measure_2" in b.supported_instructions) para encontrar backends que o suportem. Aqui, mostramos como transpilar o circuito para que as medições de meio de circuito definidas no circuito sejam executadas usando a operação MidCircuitMeasure, se o backend a suportar.

Abaixo, imprimimos a duração para a instrução measure_2 e a instrução measure padrão.

print(
    f'Mid-circuit measurement `measure_2` duration: {backend.instruction_durations.get('measure_2',0) * backend.dt * 1e9/1e3} μs'
)
print(
    f'Terminal measurement `measure` duration: {backend.instruction_durations.get('measure',0) * backend.dt *1e9/1e3} μs'
)

Mid-circuit measurement `measure_2` duration:  1.624 μs
Terminal measurement `measure` duration: 2.2 μs

"""Pass that replaces terminal measures in the middle of the circuit with
MidCircuitMeasure instructions."""

class ConvertToMidCircuitMeasure(TransformationPass):
    """This pass replaces terminal measures in the middle of the circuit with
    MidCircuitMeasure instructions.
    """

    def __init__(self, target):
        super().__init__()
        self.target = target

    def run(self, dag):
        """Run the pass on a dag."""
        mid_circ_measure = None
        for inst in self.target.instructions:
            if isinstance(inst[0], Instruction) and inst[0].name.startswith(
                "measure_"
            ):
                mid_circ_measure = inst[0]
                break
        if not mid_circ_measure:
            return dag

        final_measure_nodes = calc_final_ops(dag, {"measure"})
        for node in dag.op_nodes(Measure):
            if node not in final_measure_nodes:
                dag.substitute_node(node, mid_circ_measure, inplace=True)

        return dag

pm = PassManager(ConvertToMidCircuitMeasure(backend.target))

dynamic_isa_circuits_meas2 = [pm.run(circ) for circ in dynamic_isa_circuits]
dynamic_pubs_meas2 = [(circ, params) for circ in dynamic_isa_circuits_meas2]

dynamic_isa_circuits_dd_meas2 = [
    pm.run(circ) for circ in dynamic_isa_circuits_dd
]
dynamic_pubs_dd_meas2 = [
    (circ, params) for circ in dynamic_isa_circuits_dd_meas2
]

Transpilação para circuitos unitários

Para estabelecer uma comparação justa entre os circuitos dinâmicos e sua contraparte unitária, usamos o mesmo conjunto de qubits físicos usados nos circuitos dinâmicos para os qubits de dados como o layout para transpilar os circuitos unitários.

init_layout = [
    dynamic_layout[ind] for ind in range(circuits_unitary[0].num_qubits)
]

pm = generate_preset_pass_manager(
    target=backend.target,
    initial_layout=init_layout,
    optimization_level=3,
)

def transpile_minimize(circ: QuantumCircuit, pm: PassManager, iterations=10):
    """Transpile circuits for specified number of iterations and return the one with smallest two-qubit gate depth"""
    circs = [pm.run(circ) for i in range(iterations)]
    circs_sorted = sorted(
        circs,
        key=lambda x: x.depth(lambda x: x.operation.num_qubits == 2),
    )
    return circs_sorted[0]

unitary_isa_circuits = []
for circ in circuits_unitary:
    circ_t = transpile_minimize(circ, pm, iterations=100)
    unitary_isa_circuits.append(circ_t)

unitary_pubs = [(circ, params) for circ in unitary_isa_circuits]

Visualizamos o layout de qubit dos circuitos unitários transpilados. Os círculos pretos indicam os qubits físicos usados para transpilar os circuitos unitários e seus índices correspondem aos índices de qubit virtuais. Ao comparar isso com o layout plotado para os circuitos dinâmicos, podemos confirmar que os circuitos unitários usam o mesmo conjunto de qubits físicos que os qubits de dados nos circuitos dinâmicos.

plot_circuit_layout(
    unitary_isa_circuits[-1],
    backend,
    qubit_coordinates=_heron_coords_r2(),
    view="virtual",
)

Agora adicionamos a sequência DD aos circuitos transpilados e construímos os PUBs correspondentes para submissão de trabalho.

pm_dd = PassManager(
    [
        ALAPScheduleAnalysis(target=backend.target),
        PadDynamicalDecoupling(
            dd_sequence=[
                XGate(),
                RZGate(np.pi),
                XGate(),
                RZGate(-np.pi),
            ],
            spacing=[1 / 4, 1 / 2, 0, 0, 1 / 4],
            target=backend.target,
        ),
    ]
)

unitary_isa_circuits_dd = pm_dd.run(unitary_isa_circuits)
unitary_pubs_dd = [(circ, params) for circ in unitary_isa_circuits_dd]

Comparar profundidade de porta de dois qubits de circuitos unitários e dinâmicos

# compare circuit depth of unitary and dynamic circuit implementations
unitary_depth = [
    unitary_isa_circuits[i].depth(lambda x: x.operation.num_qubits == 2)
    for i in range(len(unitary_isa_circuits))
]

dynamic_depth = [
    dynamic_isa_circuits[i].depth(lambda x: x.operation.num_qubits == 2)
    for i in range(len(dynamic_isa_circuits))
]

plt.plot(
    list(range(len(unitary_depth))),
    unitary_depth,
    label="unitary circuits",
    color="#be95ff",
)
plt.plot(
    list(range(len(dynamic_depth))),
    dynamic_depth,
    label="dynamic circuits",
    color="#ff7eb6",
)
plt.xlabel("Trotter steps")
plt.ylabel("Two-qubit depth")
plt.legend()

<matplotlib.legend.Legend at 0x374225760>

O principal benefício do circuito baseado em medição é que, ao implementar múltiplas interações ZZ, as camadas CX podem ser paralelizadas, e as medições podem ocorrer simultaneamente. Isso ocorre porque todas as interações ZZ comutam, então o cálculo pode ser realizado com profundidade de medição 1. Após transpilar os circuitos, observamos que a abordagem de circuito dinâmico produz uma profundidade de dois qubits significativamente menor do que a abordagem unitária padrão, com a ressalva de que a medição de meio de circuito adicional e o feedforward clássico em si levam tempo e introduzem suas próprias fontes de erros.

Passo 3: Executar usando primitivos Qiskit

Modo de teste local

Antes de submeter os trabalhos ao hardware, podemos executar uma pequena simulação de teste do circuito dinâmico usando o modo de teste local.

aer_sim = AerSimulator()
pm = generate_preset_pass_manager(backend=aer_sim, optimization_level=1)
circuit_dynamic_test.measure_all()
isa_qc = pm.run(circuit_dynamic_test)
with Batch(backend=aer_sim) as batch:
    sampler = Sampler(mode=batch)
    result = sampler.run([(isa_qc, params)]).result()

print(
    "Simulated average magnetization at trotter step = 1 at three theta values"
)
result[0].data["meas"].expectation_values(obs_dynamic_test[0])

Simulated average magnetization at trotter step = 1 at three theta values

array([ 0.16666667,  0.01855469, -0.13476562])

Simulação MPS

Para circuitos grandes, podemos usar o simulador matrix_product_state (MPS), que fornece um resultado aproximado para o valor esperado de acordo com a dimensão de ligação escolhida. Posteriormente, usamos os resultados da simulação MPS como linha de base para comparar os resultados do hardware.

# The MPS simulation below took approximately 7 minutes to run on a laptop with Apple M1 chip

mps_backend = AerSimulator(
    method="matrix_product_state",
    matrix_product_state_truncation_threshold=1e-5,
    matrix_product_state_max_bond_dimension=100,
)
mps_sampler = Aer_Sampler.from_backend(mps_backend)

shots = 4096

data_sim = []
for j in range(points):
    circ_list = [
        circ.assign_parameters([params[j]]) for circ in circuits_unitary
    ]

    mps_job = mps_sampler.run(circ_list, shots=shots)
    result = mps_job.result()

    point_data = [
        result[d].data["meas"].expectation_values(observables_unitary)
        for d in depths
    ]

    data_sim.append(point_data)  # data at one theta value

data_sim = np.array(data_sim)

Com os circuitos e observáveis preparados, agora os executamos em hardware usando o primitivo Sampler.

Aqui submetemos três trabalhos para unitary_pubs, dynamic_pubs e dynamic_pubs_dd. Cada um é uma lista de circuitos parametrizados correspondendo a nove passos de Trotter diferentes com três parâmetros $\theta$ diferentes.

shots = 10000

with Batch(backend=backend) as batch:
    sampler = Sampler(mode=batch)

    sampler.options.experimental = {
        "execution": {
            "scheduler_timing": True
        },  # set to True to retrieve circuit timing info
    }

    job_unitary = sampler.run(unitary_pubs, shots=shots)
    print(f"unitary: {job_unitary.job_id()}")

    job_unitary_dd = sampler.run(unitary_pubs_dd, shots=shots)
    print(f"unitary_dd: {job_unitary_dd.job_id()}")

    job_dynamic = sampler.run(dynamic_pubs, shots=shots)
    print(f"dynamic: {job_dynamic.job_id()}")

    job_dynamic_dd = sampler.run(dynamic_pubs_dd, shots=shots)
    print(f"dynamic_dd: {job_dynamic_dd.job_id()}")

    job_dynamic_meas2 = sampler.run(dynamic_pubs_meas2, shots=shots)
    print(f"dynamic_meas2: {job_dynamic_meas2.job_id()}")

    job_dynamic_dd_meas2 = sampler.run(dynamic_pubs_dd_meas2, shots=shots)
    print(f"dynamic_dd_meas2: {job_dynamic_dd_meas2.job_id()}")

unitary: d5dtt0ldq8ts73fvbhj0
unitary: d5dtt11smlfc739onuag
dynamic: d5dtt1hsmlfc739onuc0
dynamic_dd: d5dtt25jngic73avdne0
dynamic_meas2: d5dtt2ldq8ts73fvbhm0
dynamic_dd_meas2: d5dtt2tjngic73avdnf0

Passo 4: Pós-processar e retornar resultados no formato clássico desejado

Após a conclusão dos trabalhos, podemos recuperar a duração do circuito dos metadados dos resultados do trabalho e visualizar as informações de agendamento do circuito. Para ler mais sobre visualização de informações de agendamento de um circuito, consulte esta página.

# Circuit durations is reported in the unit of `dt` which can be retrieved from `Backend` object
unitary_durations = [
    job_unitary.result()[i].metadata["compilation"]["scheduler_timing"][
        "circuit_duration"
    ]
    for i in depths
]

dynamic_durations = [
    job_dynamic.result()[i].metadata["compilation"]["scheduler_timing"][
        "circuit_duration"
    ]
    for i in depths
]

dynamic_durations_meas2 = [
    job_dynamic_meas2.result()[i].metadata["compilation"]["scheduler_timing"][
        "circuit_duration"
    ]
    for i in depths
]

result_dd = job_dynamic_dd.result()[1]
circuit_schedule_dd = result_dd.metadata["compilation"]["scheduler_timing"][
    "timing"
]

# to visualize the circuit schedule, one can show the figure below
fig_dd = draw_circuit_schedule_timing(
    circuit_schedule=circuit_schedule_dd,
    included_channels=None,
    filter_readout_channels=False,
    filter_barriers=False,
    width=1000,
)

# Save to a file since the figure is large
fig_dd.write_html("scheduler_timing_dd.html")

Plotamos as durações de circuito para circuitos unitários e os circuitos dinâmicos. Do gráfico abaixo, podemos ver que, apesar do tempo necessário para medições de meio de circuito e operações clássicas, a implementação de circuito dinâmico com measure_2 resulta em durações de circuito comparáveis à implementação unitária.

# visualize circuit durations

def convert_dt_to_microseconds(circ_duration: List, backend_dt: float):
    dt = backend_dt * 1e6  # dt in microseconds
    return list(map(lambda x: x * dt, circ_duration))

dt = backend.target.dt
plt.plot(
    depths,
    convert_dt_to_microseconds(unitary_durations, dt),
    color="#be95ff",
    linestyle=":",
    label="unitary",
)
plt.plot(
    depths,
    convert_dt_to_microseconds(dynamic_durations, dt),
    color="#ff7eb6",
    linestyle="-.",
    label="dynamic",
)
plt.plot(
    depths,
    convert_dt_to_microseconds(dynamic_durations_meas2, dt),
    color="#ff7eb6",
    linestyle="-.",
    marker="s",
    mfc="none",
    label="dynamic w/ meas2",
)

plt.xlabel("Trotter steps")
plt.ylabel(r"Circuit durations in $\mu$s")
plt.legend()

<matplotlib.legend.Legend at 0x17f73c6e0>

Após a conclusão dos trabalhos, recuperamos os dados abaixo e calculamos a magnetização média estimada pelos observáveis observables_unitary ou observables_dynamic que construímos anteriormente.

runs = {
    "unitary": (
        job_unitary,
        [observables_unitary] * len(circuits_unitary),
    ),
    "unitary_dd": (
        job_unitary_dd,
        [observables_unitary] * len(circuits_unitary),
    ),
    # Omitting Dyn w/o DD and Dynamic w/ DD plots for better readability
    # "dynamic": (job_dynamic, observables_dynamic),
    # "dynamic_dd": (job_dynamic_dd, observables_dynamic),
    "dynamic_meas2": (job_dynamic_meas2, observables_dynamic),
    "dynamic_dd_meas2": (
        job_dynamic_dd_meas2,
        observables_dynamic,
    ),
}

data_dict = {}
for key, (job, obs) in runs.items():
    data = []
    for i in range(points):
        data.append(
            [
                job.result()[ind].data["meas"].expectation_values(obs[ind])[i]
                for ind in depths
            ]
        )
    data_dict[key] = data

Abaixo plotamos a magnetização de spin como função dos passos de Trotter em diferentes valores de $\theta$, correspondendo a diferentes intensidades do campo magnético local. Plotamos tanto os resultados de simulação MPS pré-computados para os circuitos ideais unitários, juntamente com os resultados experimentais do seguinte:

1. executando os circuitos unitários com DD
2. executando os circuitos dinâmicos com DD e MidCircuitMeasure

plt.figure(figsize=(10, 6))

colors = ["#0f62fe", "#be95ff", "#ff7eb6"]
for i in range(points):
    plt.plot(
        depths,
        data_sim[i],
        color=colors[i],
        linestyle="solid",
        label=f"θ={pi_check(i*max_angle/(points-1))} (MPS)",
    )
    # plt.plot(
    #     depths,
    #     data_dict["unitary"][i],
    #     color=colors[i],
    #     linestyle=":",
    #     label=f"θ={pi_check(i*max_angle/(points-1))} (Unitary)",
    # )

    plt.plot(
        depths,
        data_dict["unitary_dd"][i],
        color=colors[i],
        marker="o",
        mfc="none",
        linestyle=":",
        label=f"θ={pi_check(i*max_angle/(points-1))} (Unitary w/DD)",
    )

    # Omitting Dyn w/o DD and Dynamic w/ DD plots for better readability
    # plt.plot(
    #     depths,
    #     data_dict["dynamic"][i],
    #     color=colors[i],
    #     linestyle="-.",
    #     label=f"θ={pi_check(i*max_angle/(points-1))} (Dyn w/o DD)",
    # )
    # plt.plot(
    #     depths,
    #     data_dict["dynamic_dd"][i],
    #     marker="D",
    #     mfc="none",
    #     color=colors[i],
    #     linestyle="-.",
    #     label=f"θ={pi_check(i*max_angle/(points-1))} (Dynamic w/ DD)",
    # )

    # plt.plot(
    #     depths,
    #     data_dict["dynamic_meas2"][i],
    #     color=colors[i],
    #     marker="s",
    #     mfc="none",
    #     linestyle=':',
    #     label=f"θ={pi_check(i*max_angle/(points-1))} (Dynamic w/ MidCircuitMeas)",
    # )

    plt.plot(
        depths,
        data_dict["dynamic_dd_meas2"][i],
        color=colors[i],
        marker="*",
        markersize=8,
        linestyle=":",
        label=f"θ={pi_check(i*max_angle/(points-1))} (Dynamic w/ DD & MidCircuitMeas)",
    )

plt.xlabel("Trotter steps", fontsize=16)
plt.ylabel("Average magnetization", fontsize=16)
plt.xticks(rotation=45)
handles, labels = plt.gca().get_legend_handles_labels()
plt.legend(
    handles,
    labels,
    loc="upper right",
    bbox_to_anchor=(1.46, 1.0),
    shadow=True,
    ncol=1,
)
plt.title(
    f"{hex_rows}x{hex_cols} hex ring, {num_qubits} data qubits, {len(ancilla)} ancilla qubits \n{backend.name}: Sampler"
)
plt.show()

Quando comparamos os resultados experimentais com a simulação, vemos que a implementação de circuito dinâmico (linha pontilhada com estrelas) no geral tem melhor desempenho do que a implementação unitária padrão (linha pontilhada com círculos). Em resumo, apresentamos circuitos dinâmicos como uma solução para simular modelos de spin Ising em uma rede em favo de mel, uma topologia que não é nativa do hardware. A solução de circuito dinâmico permite interações ZZ entre qubits que não são vizinhos mais próximos, com uma profundidade de porta de dois qubits menor do que usar portas SWAP, ao custo de introduzir qubits ancilla extras e operações de feedforward clássico.

Referências

[1] Quantum computing with Qiskit, by Javadi-Abhari, A., Treinish, M., Krsulich, K., Wood, C.J., Lishman, J., Gacon, J., Martiel, S., Nation, P.D., Bishop, L.S., Cross, A.W. and Johnson, B.R., 2024. arXiv preprint arXiv:2405.08810 (2024)