Retropropagação do operador (OBP) para estimativa de valores esperados

Estimativa de tempo de execução: 4 minutos em um processador Heron r3 (NOTA: Trata-se apenas de uma estimativa. (O tempo de execução pode variar.)

Resultados do aprendizado

Após concluir este tutorial, os usuários deverão compreender:

• Como usar qiskit-addon-obp para reduzir a profundidade do circuito quântico em troca de um maior número de execuções do circuito
• Como usar qiskit-addon-utils para construir hamiltonianos XYZ e seus circuitos de evolução temporal

Pré-requisitos

Recomendamos que os usuários estejam familiarizados com os seguintes tópicos antes de seguir com este tutorial:

• Utilização da primitiva Estimator para calcular os valores esperados de uma variável observável

Segundo plano

A retropropagação de operadores é uma técnica que consiste em incorporar operações do final de um circuito quântico no observável medido, reduzindo geralmente a profundidade do circuito à custa de termos adicionais no observável. O objetivo é realizar a retropropagação em o máximo possível do circuito, sem permitir que a observável se torne excessivamente grande. Uma implementação baseada no Qiskit está disponível no complemento OBP Qiskit. Leia a documentação correspondente para obter mais informações.

Considere um circuito de exemplo no qual se pretende medir uma observável $O = \sum_P c_P P$, sendo que $P$ são os estados de Pauli e $c_P$ são os coeficientes. Denotemos o circuito como um único " $U$ " unitário, que pode ser dividido logicamente em " $U = U_C U_Q$ ", conforme mostrado na figura abaixo.

A retropropagação do operador incorpora o operador unitário $U_C$ no observável, evoluindo-o como $O' = U_C^{\dagger}OU_C = \sum_P c_P U_C^{\dagger}PU_C$. Em outras palavras, parte do cálculo é realizada de forma clássica por meio da evolução do observável de $O$ para $O'$. O problema original pode agora ser reformulado como a medição do observável $O'$ para o novo circuito de menor profundidade, cujo operador unitário é $U_Q$.

O espaço unitário $U_C$ é representado por um conjunto de fatias $U_C = U_S U_{S-1}...U_2U_1$. Existem várias maneiras de definir uma fatia. Por exemplo, no circuito do exemplo acima, cada camada de portas do tipo “ $R_{zz}$ ” e cada camada de portas do tipo “ $R_x$ ” pode ser considerada como uma fatia individual. A retropropagação envolve, tradicionalmente, o cálculo de um $O' = \Pi_{s=1}^S \sum_P c_P U_s^{\dagger} P U_s$. Cada fatia $U_s$ pode ser representada como $U_s = exp(\frac{-i\theta_s P_s}{2})$, onde $P_s$ é um $n$ -qubit Pauli e $\theta_s$ é um escalar. É fácil verificar que

No exemplo acima, se $\{P,P_s\} = 0$, precisamos executar dois circuitos quânticos, em vez de um, para calcular o valor esperado. Portanto, a retropropagação pode aumentar o número de termos no observável, levando a um maior número de execuções do circuito. Uma maneira de permitir uma retropropagação mais profunda no circuito, evitando ao mesmo tempo que o operador se torne excessivamente grande, é truncar os termos com coeficientes pequenos, em vez de adicioná-los ao operador. Por exemplo, no exemplo acima, seria possível optar por truncar o termo que envolve $P_sP$, desde que $\theta_s$ seja suficientemente pequeno. O truncamento de termos pode reduzir o número de circuitos quânticos a serem executados, mas isso acarreta um erro no cálculo do valor esperado final, proporcional à magnitude dos coeficientes dos termos truncados.

Requisitos

Antes de iniciar este tutorial, verifique se você tem os seguintes itens instalados:

• Qiskit SDK v2.0 ou posterior, com suporte à visualização
• Qiskit Runtime v0.22 ou posterior (pip install qiskit-ibm-runtime)
• Complemento OBP Qiskit 0.3 ou versão posterior (pip install qiskit-addon-obp)
• Complemento Qiskit utils 0.3 ou posterior (pip install qiskit-addon-utils)

Instalação

import numpy as np
import matplotlib.pyplot as plt

from qiskit.primitives import StatevectorEstimator
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.quantum_info import SparsePauliOp
from qiskit.transpiler import CouplingMap
from qiskit.synthesis import LieTrotter

from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian
from qiskit_addon_utils.problem_generators import (
    generate_time_evolution_circuit,
)
from qiskit_addon_utils.slicing import slice_by_depth, combine_slices
from qiskit_addon_obp.utils.simplify import OperatorBudget
from qiskit_addon_obp import backpropagate
from qiskit_addon_obp.utils.truncating import setup_budget

from rustworkx.visualization import graphviz_draw

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import EstimatorV2, EstimatorOptions

Exemplo de simulador em pequena escala

Este tutorial implementa um padrão do Qiskit para simular a dinâmica quântica de uma cadeia de spins de Heisenberg utilizando o complemento OBP do Qiskit. Observe que, em um simulador sem ruído, o valor esperado obtido com e sem retropropagação será o mesmo.

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

Relacionar a evolução temporal de um modelo quântico de Heisenberg a um experimento quântico

Primeiro, utilizaremos a generate_xyz_hamiltonian função de qiskit-addon-utils para gerar um hamiltoniano do tipo de Heisenberg em um determinado gráfico de conectividade. Este gráfico pode ser um gráfico de barras ( rustworkx.PyGraph ) ou um gráfico de colunas ( CouplingMap ). A seguir, utilizaremos uma cadeia CouplingMap linear de 10 qubits.

num_qubits = 10
layout = [(i - 1, i) for i in range(1, num_qubits)]

# Instantiate a CouplingMap object
coupling_map = CouplingMap(layout)
graphviz_draw(coupling_map.graph, method="circo")

Output:

Em seguida, geramos um operador de Pauli que modela um hamiltoniano de Heisenberg XYZ:

onde $G(V,E)$ é a representação gráfica do mapa de acoplamento. Para este tutorial, definimos $J_x, J_y, J_z$ como $\frac{\pi}{8}, \frac{\pi}{4}, \frac{\pi}{2}$ e $h_x, h_y, h_z$ como $\frac{\pi}{3}, \frac{\pi}{6}, \frac{\pi}{9}$, respectivamente.

# Get a qubit operator describing the Heisenberg XYZ model
hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(np.pi / 8, np.pi / 4, np.pi / 2),
    ext_magnetic_field=(np.pi / 3, np.pi / 6, np.pi / 9),
)
print(hamiltonian)

Output:

SparsePauliOp(['IIIIIIIXXI', 'IIIIIIIYYI', 'IIIIIIIZZI', 'IIIIIXXIII', 'IIIIIYYIII', 'IIIIIZZIII', 'IIIXXIIIII', 'IIIYYIIIII', 'IIIZZIIIII', 'IXXIIIIIII', 'IYYIIIIIII', 'IZZIIIIIII', 'IIIIIIIIXX', 'IIIIIIIIYY', 'IIIIIIIIZZ', 'IIIIIIXXII', 'IIIIIIYYII', 'IIIIIIZZII', 'IIIIXXIIII', 'IIIIYYIIII', 'IIIIZZIIII', 'IIXXIIIIII', 'IIYYIIIIII', 'IIZZIIIIII', 'XXIIIIIIII', 'YYIIIIIIII', 'ZZIIIIIIII', 'IIIIIIIIIX', 'IIIIIIIIIY', 'IIIIIIIIIZ', 'IIIIIIIIXI', 'IIIIIIIIYI', 'IIIIIIIIZI', 'IIIIIIIXII', 'IIIIIIIYII', 'IIIIIIIZII', 'IIIIIIXIII', 'IIIIIIYIII', 'IIIIIIZIII', 'IIIIIXIIII', 'IIIIIYIIII', 'IIIIIZIIII', 'IIIIXIIIII', 'IIIIYIIIII', 'IIIIZIIIII', 'IIIXIIIIII', 'IIIYIIIIII', 'IIIZIIIIII', 'IIXIIIIIII', 'IIYIIIIIII', 'IIZIIIIIII', 'IXIIIIIIII', 'IYIIIIIIII', 'IZIIIIIIII', 'XIIIIIIIII', 'YIIIIIIIII', 'ZIIIIIIIII'],
              coeffs=[0.39269908+0.j, 0.78539816+0.j, 1.57079633+0.j, 0.39269908+0.j,
 0.78539816+0.j, 1.57079633+0.j, 0.39269908+0.j, 0.78539816+0.j,
 1.57079633+0.j, 0.39269908+0.j, 0.78539816+0.j, 1.57079633+0.j,
 0.39269908+0.j, 0.78539816+0.j, 1.57079633+0.j, 0.39269908+0.j,
 0.78539816+0.j, 1.57079633+0.j, 0.39269908+0.j, 0.78539816+0.j,
 1.57079633+0.j, 0.39269908+0.j, 0.78539816+0.j, 1.57079633+0.j,
 0.39269908+0.j, 0.78539816+0.j, 1.57079633+0.j, 1.04719755+0.j,
 0.52359878+0.j, 0.34906585+0.j, 1.04719755+0.j, 0.52359878+0.j,
 0.34906585+0.j, 1.04719755+0.j, 0.52359878+0.j, 0.34906585+0.j,
 1.04719755+0.j, 0.52359878+0.j, 0.34906585+0.j, 1.04719755+0.j,
 0.52359878+0.j, 0.34906585+0.j, 1.04719755+0.j, 0.52359878+0.j,
 0.34906585+0.j, 1.04719755+0.j, 0.52359878+0.j, 0.34906585+0.j,
 1.04719755+0.j, 0.52359878+0.j, 0.34906585+0.j, 1.04719755+0.j,
 0.52359878+0.j, 0.34906585+0.j, 1.04719755+0.j, 0.52359878+0.j,
 0.34906585+0.j])

A partir do operador do qubit, podemos gerar um circuito quântico que modela sua evolução temporal. Utilizamos generate_time_evolution_circuit a decomposição de Lie-Trotter para construir o circuito de evolução temporal.

circuit = generate_time_evolution_circuit(
    hamiltonian,
    time=0.2,
    synthesis=LieTrotter(reps=2),
)
circuit.draw("mpl", style="iqp", fold=-1)

Output:

Etapa 2: Otimizar o problema para execução em hardware quântico

Criar fatias de circuito para retropropagação

A backpropagate função realiza a retropropagação de fatias inteiras do circuito de cada vez. Portanto, a escolha da divisão em fatias pode influenciar o desempenho da retropropagação em um determinado problema. Aqui, vamos agrupar portas do mesmo tipo em faixas usando a slice_by_depth função.

Para uma explicação mais detalhada sobre o corte de circuitos, consulte este guia prático do pacote qiskit-addon-utils .

slices = slice_by_depth(circuit, max_slice_depth=1)
print(f"Separated the circuit into {len(slices)} slices.")

Output:

Separated the circuit into 18 slices.

Limitar o tamanho que o operador pode atingir durante a retropropagação

Durante a retropropagação, o número de termos no operador geralmente se aproxima rapidamente de $2^L$, onde $L$ é o número de fatias. Quando dois termos do operador não comutam em nível de qubit, precisamos de circuitos separados para obter os valores esperados correspondentes a eles. Por exemplo, se tivermos um observável de dois qubits $O = 0.1 XX + 0.3 IZ - 0.5 IX$, então, como $[XX,IX] = 0$, a medição em uma única base é suficiente para calcular os valores esperados desses dois termos. No entanto, o operador $IZ$ é anticomutativo em relação aos outros dois termos; portanto, precisamos de uma medição de base separada para calcular o valor esperado de $IZ$. Em outras palavras, precisamos de dois circuitos em vez de um para calcular $\langle O \rangle$. À medida que o número de termos no operador aumenta, existe a possibilidade de que o número necessário de execuções do circuito também aumente.

O tamanho do operador pode ser limitado especificando o kwarg operator_budget da função backpropagate , que aceita uma OperatorBudget instância.

Para controlar a quantidade de recursos adicionais (número de execuções do circuito e, consequentemente, o tempo de QPU necessário) alocados, restringimos o número máximo de grupos de Pauli comutativos por qubit que o observável retropropagado pode ter. Especificamos aqui que a retropropagação deve ser interrompida quando o número de grupos de Pauli comutativos por qubit no operador ultrapassar oito.

op_budget = OperatorBudget(max_qwc_groups=8)

Retropropagação de fatias do circuito

Primeiro, especificamos o observável a ser $M_Z = \frac{1}{N} \sum_{i=1}^N \langle Z_i \rangle$, sendo $N$ o número de qubits. Faremos a retropropagação das fatias do circuito de evolução temporal até que os termos no observável não possam mais ser combinados em oito ou menos grupos Pauli comutáveis de qubit-wise.

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

Output:

SparsePauliOp(['IIIIIIIIIZ', 'IIIIIIIIZI', 'IIIIIIIZII', 'IIIIIIZIII', 'IIIIIZIIII', 'IIIIZIIIII', 'IIIZIIIIII', 'IIZIIIIIII', 'IZIIIIIIII', 'ZIIIIIIIII'],
              coeffs=[0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j,
 0.1+0.j, 0.1+0.j])

Abaixo, você verá que realizamos a retropropagação em seis fatias, e os termos foram agrupados em seis, e não em oito grupos. Isso significa que a retropropagação de mais uma fatia faria com que o número de grupos de Pauli ultrapassasse oito. Podemos verificar se é esse o caso analisando os metadados retornados. Observe também que, nesta parte, a transformação do circuito é exata. Ou seja, nenhum termo do novo observável $O’$ foi truncado. O circuito com retropropagação e o operador com retropropagação fornecem o mesmo resultado exato que o circuito e o operador originais.

# Backpropagate slices onto the observable
bp_obs, remaining_slices, metadata = backpropagate(
    observable, slices, operator_budget=op_budget
)
# Recombine the slices remaining after backpropagation
bp_circuit = combine_slices(remaining_slices)

print(f"Backpropagated {metadata.num_backpropagated_slices} slices.")
print(
    f"New observable has {len(bp_obs.paulis)} terms, which can be combined into "
    f"{len(bp_obs.group_commuting(qubit_wise=True))} groups."
)
print(
    f"Note that backpropagating one more slice would result in "
    f"{metadata.backpropagation_history[-1].num_paulis[0]} terms "
    f"across {metadata.backpropagation_history[-1].num_qwc_groups} groups."
)
print("The remaining circuit after backpropagation looks as follows:")
bp_circuit.draw("mpl", fold=-1, scale=0.6)

Output:

Backpropagated 6 slices.
New observable has 60 terms, which can be combined into 6 groups.
Note that backpropagating one more slice would result in 114 terms across 12 groups.
The remaining circuit after backpropagation looks as follows:

Para o exemplo em pequena escala no simulador, não utilizaremos o truncamento. Isso ocorre porque, na ausência de ruído, o circuito com e sem retropropagação leva ao mesmo resultado, e o truncamento piora o resultado devido à aproximação adicional.

Transpilar os circuitos para o conjunto de portas de base

Agora, compilamos tanto o circuito original quanto o circuito resultante da retropropagação para a porta de base do backend. Não precisamos fazer a transpilagem no backend propriamente dito, já que vamos rodar em um simulador para a instância pequena.

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=133
)
print(backend)

Output:

<IBMBackend('ibm_kingston')>

pm_basis = generate_preset_pass_manager(
    optimization_level=3, basis_gates=backend.configuration().basis_gates
)
isa_circuit = pm_basis.run(circuit)
isa_bp_circuit = pm_basis.run(bp_circuit)

Passo 3: Execute usando Qiskit primitives

Primeiro, criamos dois Blocos Primitivos Unificados (PUBs) correspondentes ao circuito original e ao circuito resultante da retropropagação. Em seguida, aplicamos os intervalos de confiança a um estimador ideal para obter os valores esperados.

pubs = [(isa_circuit, observable), (isa_bp_circuit, bp_obs)]

rng = np.random.default_rng()
estimator = StatevectorEstimator(seed=rng)
job = estimator.run(pubs)

Etapa 4: Realizar o pós-processamento e apresentar o resultado no formato clássico desejado

Agora, calculamos os valores esperados dos circuitos original e retropropagado.

primitive_result = job.result()
circuit_expval = primitive_result[0].data.evs.item()
bp_circuit_expval = primitive_result[1].data.evs.item()

methods = [
    "No backpropagation",
    "Backpropagation",
]
values = [circuit_expval, bp_circuit_expval]

ax = plt.gca()
plt.bar(methods, values, color="#a56eff", width=0.4, edgecolor="#8a3ffc")
ax.set_ylim([0.6, 0.92])
ax.set_ylabel(r"$M_Z$", fontsize=12)

Output:

Text(0, 0.5, '$M_Z$')

Como era de se esperar, os dois valores esperados coincidem. Como estamos utilizando um simulador de vetor de estado sem ruído, a retropropagação é uma transformação exata do par circuito-observável; portanto, os fluxos de trabalho original e retropropagado devem produzir o mesmo valor de $M_Z$. A vantagem da retropropagação só se torna evidente em hardware com ruído, onde o circuito retropropagado mais curto acumula menos erro, conforme ilustrado no exemplo de hardware em grande escala abaixo.

Exemplo de hardware em grande escala

Ao desenvolver um experimento, é útil começar com um circuito simples para facilitar as visualizações e simulações. Agora, analisamos a retropropagação de operadores para um hamiltoniano de Heisenberg de 50 qubits com o mesmo conjunto de valores para os parâmetros $J$ e $h$ e o mesmo observável $M_Z$, mas para quatro passos de Trotter. O valor esperado ideal nessa escala não pode ser calculado por um método de força bruta; por isso, utilizamos uma rede tensorial e obtemos o valor esperado ideal como sendo $\simeq 0.89$.

Além da retropropagação, neste exemplo em grande escala, também apresentamos a retropropagação com truncamento. O ideal é realizar a retropropagação o máximo possível para reduzir a profundidade do circuito efetivo. No entanto, isso muitas vezes resulta em um grande número de termos não comutativos no observável atualizado, aumentando a sobrecarga quântica. Portanto, podemos eliminar os termos observáveis com coeficientes pequenos por meio de uma técnica chamada truncamento. Embora o truncamento permita uma maior propagação ao reduzir o número de termos no observável atualizado, ele também introduz alguma aproximação. Portanto, é necessário restringir o truncamento dentro de certos limites, para que o erro de aproximação não supere a redução de ruído obtida com uma retropropagação mais profunda.

Para limitar o grau de truncamento, atribuímos um limite de erro para cada fatia, bem como um limite de erro total para todo o circuito submetido à retropropagação, utilizando a setup_budget função. Isso garante que o truncamento seja controlado tanto para cada fatia quanto para todo o circuito. Consulte também este guia para conhecer outras formas de distribuir o orçamento.

num_qubits = 50
layout = [(i - 1, i) for i in range(1, num_qubits)]

# Instantiate a CouplingMap object
coupling_map = CouplingMap(layout)

hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(np.pi / 8, np.pi / 4, np.pi / 2),
    ext_magnetic_field=(np.pi / 3, np.pi / 6, np.pi / 9),
)

# Generate a time evolution circuit for the Hamiltonian
circuit = generate_time_evolution_circuit(
    hamiltonian,
    time=0.2,
    synthesis=LieTrotter(reps=4),
)

# Define the observable to measure
observable = SparsePauliOp.from_sparse_list(
    [("Z", [i], 1 / num_qubits) for i in range(num_qubits)],
    num_qubits,
)

slices = slice_by_depth(circuit, max_slice_depth=1)

# Define the maximum number of qwc groups allowed in the
# backpropagated observable,
# and the truncation error budget
op_budget = OperatorBudget(max_qwc_groups=15)
truncation_error_budget = setup_budget(
    max_error_total=0.03, max_error_per_slice=0.005
)

# First backpropagation without truncation
bp_obs, remaining_slices, metadata = backpropagate(
    observable, slices, operator_budget=op_budget
)
bp_circuit = combine_slices(remaining_slices)

# Now backpropagate with truncation, using the same operator budget and
# the defined truncation error budget
bp_obs_trunc, remaining_slices_trunc, metadata = backpropagate(
    observable,
    slices,
    operator_budget=op_budget,
    truncation_error_budget=truncation_error_budget,
)
bp_circuit_trunc = combine_slices(
    remaining_slices_trunc, include_barriers=False
)

# Now we transpile the original circuit and the two backpropagated circuits,
# and apply the layout to the corresponding observables
pm = generate_preset_pass_manager(optimization_level=3, backend=backend)

isa_circuit = pm.run(circuit)
isa_bp_circuit = pm.run(bp_circuit)
isa_bp_circuit_trunc = pm.run(bp_circuit_trunc)

isa_observable = observable.apply_layout(isa_circuit.layout)
isa_bp_observable = bp_obs.apply_layout(isa_bp_circuit.layout)
isa_bp_observable_trunc = bp_obs_trunc.apply_layout(
    isa_bp_circuit_trunc.layout
)

# Compare the 2-qubit depth of each transpiled circuit to see how much
# depth backpropagation saved
print(
    f"2-qubit depth without backpropagation: "
    f"{isa_circuit.depth(lambda x: x.operation.num_qubits == 2)}"
)
print(
    f"2-qubit depth with backpropagation: "
    f"{isa_bp_circuit.depth(lambda x: x.operation.num_qubits == 2)}"
)
print(
    f"2-qubit depth with backpropagation and truncation: "
    f"{isa_bp_circuit_trunc.depth(lambda x: x.operation.num_qubits == 2)}"
)

pubs = [
    (isa_circuit, isa_observable),
    (isa_bp_circuit, isa_bp_observable),
    (isa_bp_circuit_trunc, isa_bp_observable_trunc),
]

# Now we instantiate the Estimator primitive for the hardware with
# ZNE and measurement error
# mitigation and compute the three circuits and observables
options = EstimatorOptions()
options.default_precision = 0.01
options.resilience_level = 2
options.resilience.zne.noise_factors = [1, 1.2, 1.4]
options.resilience.zne.extrapolator = ["linear"]
estimator = EstimatorV2(mode=backend, options=options)

estimator.options.environment.job_tags = ["TUT_OBP"]
job = estimator.run(pubs)

# Retrieve the results and the standard deviations
result_no_bp = job.result()[0].data.evs.item()
result_bp = job.result()[1].data.evs.item()
result_bp_trunc = job.result()[2].data.evs.item()

std_no_bp = job.result()[0].data.stds.item()
std_bp = job.result()[1].data.stds.item()
std_bp_trunc = job.result()[2].data.stds.item()

Output:

2-qubit depth without backpropagation: 24
2-qubit depth with backpropagation: 20
2-qubit depth with backpropagation and truncation: 18

print(f"Expectation value without backpropagation: {result_no_bp}")
print(f"Backpropagated expectation value: {result_bp}")
print(f"Backpropagated expectation value with truncation: {result_bp_trunc}")

Output:

Expectation value without backpropagation: 0.9543907942381811
Backpropagated expectation value: 0.9445337385406468
Backpropagated expectation value with truncation: 0.934050286970965

# Plot the results
methods = [
    "No backpropagation",
    "Backpropagation",
    "Backpropagation w/ truncation",
]
values = [result_no_bp, result_bp, result_bp_trunc]
error_bars = [std_no_bp, std_bp, std_bp_trunc]

ax = plt.gca()
plt.bar(methods, values, color="#a56eff", width=0.4, edgecolor="#8a3ffc")
plt.errorbar(methods, values, yerr=error_bars, fmt="o", color="r", capsize=5)
plt.axhline(0.89)
ax.set_ylim([0.8, 0.98])
plt.text(0.25, 0.895, "Exact result")
ax.set_ylabel(r"$M_Z$", fontsize=12)

Output:

Text(0, 0.5, '$M_Z$')

Próximas etapas

Se você achou este trabalho interessante, talvez se interesse pelo seguinte material: