Aprimorando a estimativa de energia de um modelo de rede fermiônico com SQD
Neste tutorial, implementamos um padrão Qiskit que mostra como pós-processar amostras quânticas ruidosas para encontrar uma aproximação do estado fundamental de um Hamiltoniano de rede fermiônico conhecido como modelo de Anderson de impureza única. Seguiremos uma abordagem de diagonalização quântica baseada em amostras para processar amostras tomadas de um conjunto de estados de base de Krylov de 16 qubits sobre intervalos de tempo crescentes. Esses estados são realizados no dispositivo quântico usando Trotterização da evolução temporal. Para considerar o efeito do ruído quântico, é usada a técnica de recuperação de configuração. Assumindo um bom estado inicial e a esparsidade do estado fundamental, esta abordagem é provada convergir eficientemente.
O padrão pode ser descrito em quatro etapas:
- Etapa 1: Mapear para o problema quântico
- Gerar um conjunto de estados de base de Krylov (ou seja, circuitos de evolução temporal Trotterizados) sobre intervalos de tempo crescentes para estimar o estado fundamental
- Etapa 2: Otimizar o problema
- Transpilar os circuitos para o backend
- Etapa 3: Executar experimentos
- Tirar amostras dos circuitos usando a primitiva
Sampler
- Tirar amostras dos circuitos usando a primitiva
- Etapa 4: Pós-processar resultados
- Loop de recuperação de configuração autoconsistente
- Pós-processar o conjunto completo de amostras de bitstrings, usando conhecimento prévio do número de partículas e da ocupação orbital média calculada na iteração mais recente
- Criar probabilisticamente lotes de subamostras a partir das bitstrings recuperadas
- Projetar e diagonalizar o Hamiltoniano de rede fermiônico sobre cada subespaço amostrado
- Salvar a energia mínima do estado fundamental encontrada em todos os lotes e atualizar a ocupação orbital média
- Loop de recuperação de configuração autoconsistente
Etapa 1: Mapear o problema para um circuito quântico
Primeiro, criaremos os Hamiltonianos de um e dois corpos descrevendo o modelo de Anderson de impureza única (SIAM) unidimensional com 7 sítios de banho (8 elétrons em 8 orbitais). Este modelo é usado para descrever impurezas magnéticas embebidas em metais.
Em seguida, criaremos os circuitos Trotter de 16 qubits usados para gerar o subespaço quântico de Krylov.
# Added by doQumentation — required packages for this notebook
!pip install -q ffsim matplotlib numpy qiskit qiskit-addon-sqd qiskit-ibm-runtime scipy
import numpy as np
n_bath = 7 # number of bath sites
V = 1 # hybridization amplitude
t = 1 # bath hopping amplitude
U = 10 # Impurity onsite repulsion
eps = -U / 2 # Chemical potential for the impurity
# Place the impurity on the first qubit
impurity_index = 0
# One body matrix elements in the "position" basis
h1e = -t * np.diag(np.ones(n_bath), k=1) - t * np.diag(np.ones(n_bath), k=-1)
h1e[impurity_index, impurity_index + 1] = -V
h1e[impurity_index + 1, impurity_index] = -V
h1e[impurity_index, impurity_index] = eps
# Two body matrix elements in the "position" basis
h2e = np.zeros((n_bath + 1, n_bath + 1, n_bath + 1, n_bath + 1))
h2e[impurity_index, impurity_index, impurity_index, impurity_index] = U
Em seguida, geraremos o subespaço quântico de Krylov com um conjunto de circuitos quânticos Trotterizados. Aqui criamos auxiliares para gerar o estado inicial (de referência), bem como a evolução temporal para as partes de um e dois corpos do Hamiltoniano. Para uma descrição mais detalhada deste modelo e como os circuitos são projetados, consulte o artigo.
import ffsim
import scipy
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit.library import CPhaseGate, XGate, XXPlusYYGate
n_modes = n_bath + 1
nelec = (n_modes // 2, n_modes // 2)
dt = 0.2
Utar = scipy.linalg.expm(-1j * dt * h1e)
# The reference state
def initial_state(q_circuit, norb, nocc):
"""Prepare an initial state."""
for i in range(nocc):
q_circuit.append(XGate(), [i])
q_circuit.append(XGate(), [norb + i])
rot = XXPlusYYGate(np.pi / 2, -np.pi / 2)
for i in range(3):
for j in range(nocc - i - 1, nocc + i, 2):
q_circuit.append(rot, [j, j + 1])
q_circuit.append(rot, [norb + j, norb + j + 1])
q_circuit.append(rot, [j + 1, j + 2])
q_circuit.append(rot, [norb + j + 1, norb + j + 2])
# The one-body time evolution
free_fermion_evolution = ffsim.qiskit.OrbitalRotationJW(n_modes, Utar)
# The two-body time evolution
def append_diagonal_evolution(dt, U, impurity_qubit, num_orb, q_circuit):
"""Append two-body time evolution to a quantum circuit."""
if U != 0:
q_circuit.append(
CPhaseGate(-dt / 2 * U),
[impurity_qubit, impurity_qubit + num_orb],
)
Gere d estados evoluídos no tempo que especificam o subespaço quântico de Krylov. Aqui, escolhemos d=8. O erro proveniente da amostragem de estados de base de Krylov converge com o aumento de d. Note que as particularidades desta instância do problema permitem uma compilação eficiente da evolução de um corpo usando OrbitalRotationJW; no entanto, em geral, pode-se usar o PauliEvolutionGate do Qiskit para implementar a evolução temporal Trotterizada do Hamiltoniano completo.
# Generate the initial state
qubits = QuantumRegister(2 * n_modes, name="q")
init_state = QuantumCircuit(qubits)
initial_state(init_state, n_modes, n_modes // 2)
init_state.draw("mpl", scale=0.4, fold=-1)
d = 8 # Number of Krylov basis states
circuits = []
for i in range(d):
circ = init_state.copy()
circuits.append(circ)
for _ in range(i):
append_diagonal_evolution(dt, U, impurity_index, n_modes, circ)
circ.append(free_fermion_evolution, qubits)
append_diagonal_evolution(dt, U, impurity_index, n_modes, circ)
circ.measure_all()
circuits[0].draw("mpl", scale=0.4, fold=-1)
circuits[-1].draw("mpl", scale=0.4, fold=-1)
Etapa 2: Otimizar o problema
Após criar os circuitos Trotterizados, vamos otimizá-los para um hardware alvo. Precisamos escolher o dispositivo de hardware a ser usado antes da otimização. Usaremos um backend falso de 127 qubits do qiskit_ibm_runtime para emular um dispositivo real. Para executar em uma QPU real, basta substituir o backend falso por um backend real. Confira a documentação do Qiskit IBM Runtime para mais informações.
from qiskit_ibm_runtime.fake_provider.backends import FakeSherbrooke
backend = FakeSherbrooke()
Em seguida, transpilaremos os circuitos para o backend alvo usando o Qiskit.
from qiskit.transpiler import generate_preset_pass_manager
# The circuit needs to be transpiled to the AerSimulator target
pass_manager = generate_preset_pass_manager(optimization_level=3, backend=backend)
isa_circuits = pass_manager.run(circuits)
Etapa 3: Executar experimentos
Após otimizar os circuitos para a execução em hardware, estamos prontos para executá-los no hardware alvo e coletar amostras para a estimativa da energia do estado fundamental. Aqui usamos SamplerV2 do qiskit-ibm-runtime para simular amostras ruidosas tomadas do backend ibm_sherbrooke. Em seguida, combinamos as contagens de cada um dos estados de base de Krylov em um único dicionário de contagens e plotamos as 20 bitstrings amostradas com mais frequência.
Nota: O Qiskit Aer é necessário para simular amostras a partir de circuitos transpilados.
from qiskit.primitives import BitArray
from qiskit.visualization import plot_histogram
from qiskit_ibm_runtime import SamplerV2 as Sampler
# Sample from the circuits
noisy_sampler = Sampler(backend, options={"simulator": {"seed_simulator": 24}})
job = noisy_sampler.run(isa_circuits, shots=500)
# Combine the counts from the individual Trotter circuits
bit_array = BitArray.concatenate_shots([result.data.meas for result in job.result()])
plot_histogram(bit_array.get_counts(), number_to_keep=20)
Etapa 4: Pós-processar os resultados
Agora, executamos o algoritmo SQD usando a função diagonalize_fermionic_hamiltonian. Veja a documentação da API para explicações sobre os argumentos desta função.
from qiskit_addon_sqd.fermion import SCIResult, diagonalize_fermionic_hamiltonian
# List to capture intermediate results
result_history = []
def callback(results: list[SCIResult]):
result_history.append(results)
iteration = len(result_history)
print(f"Iteration {iteration}")
for i, result in enumerate(results):
print(f"\tSubsample {i}")
print(f"\t\tEnergy: {result.energy}")
print(f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}")
rng = np.random.default_rng(24)
result = diagonalize_fermionic_hamiltonian(
h1e,
h2e,
bit_array,
samples_per_batch=300,
norb=n_modes,
nelec=nelec,
num_batches=3,
max_iterations=10,
symmetrize_spin=True,
callback=callback,
seed=rng,
)
Iteration 1
Subsample 0
Energy: -13.257128325607997
Subspace dimension: 3969
Subsample 1
Energy: -13.257128325607997
Subspace dimension: 3969
Subsample 2
Energy: -13.257128325607997
Subspace dimension: 3969
Iteration 2
Subsample 0
Energy: -13.319666127542039
Subspace dimension: 4096
Subsample 1
Energy: -13.420534292304595
Subspace dimension: 4624
Subsample 2
Energy: -9.136171594591085
Subspace dimension: 4624
Iteration 3
Subsample 0
Energy: -13.422491814612833
Subspace dimension: 4900
Subsample 1
Energy: -13.422491814612833
Subspace dimension: 4900
Subsample 2
Energy: -13.422491814612833
Subspace dimension: 4900
Iteration 4
Subsample 0
Energy: -13.422491814612833
Subspace dimension: 4900
Subsample 1
Energy: -13.422491814612833
Subspace dimension: 4900
Subsample 2
Energy: -13.422491814612833
Subspace dimension: 4900
Agora, plotamos os resultados. O primeiro gráfico mostra que após algumas iterações obtemos a energia exata do estado fundamental.
Este exemplo é pequeno o suficiente para que sejamos capazes de explorar o espaço de Hilbert completo, como visto nas instruções de impressão acima. Lembre-se de que o espaço de Hilbert completo tem dimensão (num_orbitals choose nelec_a) * (num_orbitals choose nelec_b). Portanto, para este problema: (8 choose 4)**2 = 4900. As dimensões do subespaço aumentam devido à recuperação de configuração aprimorada e também ao fato de o algoritmo SQD carregar configurações importantes de uma iteração para a próxima. Na última iteração, estamos diagonalizando em todo o espaço de Hilbert.
O segundo gráfico mostra a ocupação média de cada orbital espacial em todas as soluções dos lotes. Vemos que com alta probabilidade cada orbital contém um elétron.
import matplotlib.pyplot as plt
exact_energy = -13.422491814605827
min_es = [min(result, key=lambda res: res.energy).energy for result in result_history]
min_id, min_e = min(enumerate(min_es), key=lambda x: x[1])
# Data for energies plot
x1 = range(len(result_history))
yt1 = list(np.arange(-13.5, -13.1, 0.1))
ytl = [f"{i:.1f}" for i in yt1]
# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))
fig, axs = plt.subplots(1, 2, figsize=(12, 6))
# Plot energies
axs[0].plot(x1, min_es, label="SQD energy", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yticks(yt1)
axs[0].set_yticklabels(ytl)
axs[0].axhline(y=exact_energy, color="#BF5700", linestyle="--", label="Exact energy")
axs[0].set_title("Approximated Ground State Energy vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy", fontdict={"fontsize": 12})
axs[0].legend()
# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})
print(f"Exact energy: {exact_energy:.5f}")
print(f"SQD energy: {min_e:.5f}")
print(f"Absolute error: {abs(min_e - exact_energy):.5f}")
plt.tight_layout()
plt.show()
Exact energy: -13.42249
SQD energy: -13.42249
Absolute error: 0.00000
