Estimativa da energia do estado fundamental da cadeia de Heisenberg com VQE

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

Resultados do aprendizado

Ao concluir este tutorial, você deverá compreender as seguintes informações:

• Como modelar uma cadeia de espín de Heisenberg como um hamiltoniano quântico usando o Qiskit
• Como usar o otimizador SPSA para estimar a energia do estado fundamental de um sistema quântico
• Como executar fluxos de trabalho variacionais em um hardware quântico d IBM®, utilizando primitivas e sessões d Qiskit Runtime

Pré-requisitos

Recomenda-se que você se familiarize com estes tópicos:

• Noções básicas de informação quântica
• Introdução aos padrões do Qiskit
• Design de algoritmos variacionais

Segundo plano

A cadeia de espín de Heisenberg é um dos modelos mais estudados na física da matéria condensada e no magnetismo quântico. Descreve uma rede unidimensional de spins quânticos em interação, na qual os spins vizinhos mais próximos estão acoplados por meio de interações de troca. O hamiltoniano do modelo de Heisenberg isotrópico com um campo magnético externo é dado por:

$H = \sum_{\langle i,j \rangle} \left( J_x X_i X_j + J_y Y_i Y_j + J_z Z_i Z_j \right) + \sum_{i} h_i Z_i,$

onde $X_i$, $Y_i$ e $Z_i$ são os operadores de Pauli que atuam no sítio $i$, a soma $\langle i,j \rangle$ abrange os pares de vizinhos mais próximos, $J_x = J_y = J_z = 0.5$ são as constantes de acoplamento de troca (isotrópicas neste tutorial) e $h_i$ representa um campo magnético externo dependente do sítio. Neste tutorial, os valores do campo magnético são amostrados aleatoriamente no intervalo $[-1, 1]$. Observe que, na implementação abaixo, o conjunto de pares de “vizinhos mais próximos” é determinado pelo acoplamento nativo do backend de hardware entre os primeiros qubits $N$, o que pode não formar uma cadeia linear estrita, dependendo da topologia do dispositivo.

Compreender a energia do estado fundamental desse hamiltoniano é de importância fundamental na física. O estado fundamental contém informações sobre transições de fase quânticas, estrutura de entrelaçamento e ordenação magnética. Tradicionalmente, o cálculo da energia exata do estado fundamental torna-se impraticável à medida que o número de spins aumenta, uma vez que a dimensão do espaço de Hilbert varia exponencialmente como $2^N$ para $N$ spins. Isso o torna um candidato natural para a simulação quântica.

O Variational Quantum Eigensolver (VQE) é um algoritmo híbrido quântico-clássico projetado para estimar a energia do estado fundamental de um hamiltoniano. O método funciona preparando-se um estado quântico parametrizado $|\psi(\theta)\rangle$ (denominado ansatz) em um computador quântico e medindo-se o valor esperado $\langle \psi(\theta) | H | \psi(\theta) \rangle$. Em seguida, um otimizador clássico ajusta iterativamente os parâmetros $\theta$ para minimizar essa energia, aproveitando o princípio variacional, que garante que a energia medida seja sempre um limite superior da verdadeira energia do estado fundamental.

Neste tutorial, utilizamos o efficient_su2 ansatz da biblioteca de circuitos do Qiskit, que constrói camadas de rotações de qubits únicos e portas de entrelaçamento. A otimização é realizada utilizando o algoritmo de Aproximação Estocástica por Perturbação Simultânea (SPSA), que é adequado para hardware quântico sujeito a ruídos, pois estima gradientes utilizando apenas duas avaliações da função por iteração, independentemente do número de parâmetros.

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.44 ou posterior (pip install qiskit-ibm-runtime)

Instalação

import numpy as np
import matplotlib.pyplot as plt
from typing import Sequence

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.primitives import BaseEstimatorV2
from qiskit.circuit.library import XGate
from qiskit.circuit.library import efficient_su2
from qiskit.transpiler import PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.transpiler.passes.scheduling import (
    ALAPScheduleAnalysis,
    PadDynamicalDecoupling,
)
from qiskit_ibm_runtime import QiskitRuntimeService, Session, EstimatorV2


def visualize_results(results):
    plt.plot(results["cost_history"], lw=2)
    plt.xlabel("Number of function evaluations")
    plt.ylabel("Energy")
    plt.show()

Exemplo em pequena escala

Nesta seção, vamos percorrer cada etapa do padrão Qiskit em pequena escala, explicando os principais componentes à medida que construímos o fluxo de trabalho.

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

• Entrada: Número de giros
• Resultado: Ansatz e Hamiltoniano modelando a cadeia de Heisenberg

Construa um ansatz e um hamiltoniano que modelem uma cadeia de Heisenberg de 10 spins. Nesta etapa, construiremos um hamiltoniano de Heisenberg de 10 spins sobre o mapa de acoplamento do backend menos ocupado e prepararemos o efficient_su2 ansatz.

num_spins = 10
ansatz = efficient_su2(num_qubits=num_spins, reps=2)

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

coupling = backend.target.build_coupling_map()
reduced_coupling = coupling.reduce(list(range(num_spins)))

edge_list = reduced_coupling.graph.edge_list()
ham_list = []

for edge in edge_list:
    ham_list.append(("ZZ", edge, 0.5))
    ham_list.append(("YY", edge, 0.5))
    ham_list.append(("XX", edge, 0.5))

for qubit in reduced_coupling.physical_qubits:
    ham_list.append(("Z", [qubit], np.random.random() * 2 - 1))

hamiltonian = SparsePauliOp.from_sparse_list(ham_list, num_qubits=num_spins)

ansatz.draw("mpl", style="iqp")

Output:

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

• Entrada: Circuito abstrato, observável
• Saída: Circuito-alvo e observável, otimizado para a QPU selecionada

Use a função generate_preset_pass_manager do Qiskit para gerar automaticamente uma rotina de otimização para o nosso circuito com relação à QPU selecionada. Escolhemos optimization_level=3, que oferece o mais alto nível de otimização dos gerenciadores de passagem predefinidos. Incluímos também ALAPScheduleAnalysis e PadDynamicalDecoupling passes de agendamento para suprimir erros de decoerência.

target = backend.target
pm = generate_preset_pass_manager(optimization_level=3, target=target)
pm.scheduling = PassManager(
    [
        ALAPScheduleAnalysis(durations=target.durations()),
        PadDynamicalDecoupling(
            durations=target.durations(),
            dd_sequence=[XGate(), XGate()],
            pulse_alignment=target.pulse_alignment,
        ),
    ]
)
isa_ansatz = pm.run(ansatz)
isa_observable = hamiltonian.apply_layout(isa_ansatz.layout)
isa_ansatz.draw("mpl", scale=0.6, style="iqp", fold=-1, idle_wires=False)

Output:

Passo 3: Execute usando Qiskit primitives

• Entrada: Circuito alvo e observável
• Saída: Resultados da otimização

Minimize a energia estimada do estado fundamental do sistema otimizando os parâmetros do circuito. Use a Estimator primitiva disponível em Qiskit Runtime para calcular a função de custo durante a otimização.

Como otimizamos o circuito para o backend na Etapa 2, podemos evitar a transpilagem no servidor de tempo de execução definindo skip_transpilation=True e passando o circuito otimizado. Para esta demonstração, vamos executar em uma QPU usando qiskit-ibm-runtime primitivas. Para executar com primitivas qiskit baseadas em vetores de estado, substitua o bloco de código que utiliza primitivas do tipo Qiskit Runtime pelo bloco comentado.

Neste tutorial, utilizamos a Aproximação Estocástica por Perturbação Simultânea (SPSA), que é um otimizador baseado no gradiente. A seguir, apresentamos uma breve introdução ao assunto e fornecemos o código para implementar o SPSA usando o Qiskit v2.0.

Apresentando a SPSA

A Aproximação Estocástica por Perturbação Simultânea (SPSA) [1] é um algoritmo de otimização que aproxima todo o vetor gradiente utilizando apenas duas chamadas de função em cada iteração. Seja $f:\mathbb{R}^p\rightarrow \mathbb{R}$ a função de custo com $p$ parâmetros a serem otimizados, e $x_i\in \mathbb{R}^p$ o vetor de parâmetros no $i^{th}$ passo da iteração. Para calcular o gradiente, é criado um vetor aleatório $\Delta_i$ de dimensão $p$, em que cada elemento $\Delta_{ij}$, $\forall$ $j\in \{1,2,...,p\}$ é amostrado uniformemente a partir de $\{-1, 1\}$. Em seguida, cada elemento do vetor aleatório $\Delta_i$ é multiplicado por um pequeno valor $c_i$ para criar uma perturbação aleatória. O gradiente é então estimado como

$[\nabla f(x_i)]_j \approx \frac{f(x_i + c_i \Delta_i) - f(x_i - c_i \Delta_i)}{2c_i\Delta_{ij}}.$

Intuitivamente, uma vez que uma perturbação aleatória é aplicada durante a estimativa do gradiente, espera-se que pequenos desvios nos valores exatos de um $f$, decorrentes do ruído, possam ser tolerados e levados em conta. Na verdade, o SPSA é especialmente conhecido por sua robustez contra ruídos e requer apenas duas chamadas de hardware por iteração. É, portanto, um dos otimizadores mais utilizados para a implementação de algoritmos variacionais.

Neste tutorial, os hiperparâmetros para a iteração do algoritmo de otimização por gradiente descendente ( $i^{th}$ ), $a_i$ e $c_i$, são calculados da seguinte forma:

$a_i = \frac{a}{(A + i + 1)^\alpha} \quad \text{and} \quad c_i = \frac{c}{(i+1)^\gamma},$

onde os valores constantes são tomados como $A = 30$, $\alpha = 0.9$, $a = 0.3$, $c = 0.1$ e $\gamma = 0.4$. Esses valores foram selecionados a partir de [2]. É necessário ajustar adequadamente os hiperparâmetros para obter um bom desempenho do SPSA.

def spsa(
    fun, x0, args=(), A=30, alpha=0.9, a=0.3, c=0.1, gamma=0.4, maxiter=100
):
    nparams = len(x0)
    x = np.copy(x0)

    for i in range(maxiter):
        a_i = a / (A + i + 1) ** alpha
        c_i = c / (i + 1) ** gamma
        delta_i = np.random.choice([-1, 1], nparams)

        # two hardware calls
        eval_1 = fun(x + c_i * delta_i, *args)
        eval_2 = fun(x - c_i * delta_i, *args)

        # compute the gradient and update the parameters
        grad = (eval_1 - eval_2) / (2 * c_i) * np.reciprocal(delta_i)
        x = x - a_i * grad

    return x

def cost_func(
    params: Sequence,
    ansatz: QuantumCircuit,
    hamiltonian: SparsePauliOp,
    estimator: BaseEstimatorV2,
    cost_history_dict: dict,
) -> float:
    """Ground state energy evaluation."""
    energy = (
        estimator.run([(ansatz, hamiltonian, [params])]).result()[0].data.evs
    )

    cost_history_dict["iters"] += 1
    cost_history_dict["prev_vector"] = list(params)
    cost_history_dict["cost_history"].append(float(energy[0]))

    print(
        f"Fx Iters. done: {cost_history_dict['iters']} [Current cost: {round(energy[0], 5)}]",
        end="\r",
    )

    return energy


def solve(x0, isa_ansatz, isa_observable, maxiter=150):
    cost_history_dict = {
        "prev_vector": None,
        "iters": 0,
        "cost_history": [],
        "y_min": None,
    }

    # Evaluate the problem using a QPU via Qiskit IBM Runtime
    with Session(backend=backend) as session:
        estimator = EstimatorV2(mode=session)
        estimator.skip_transpilation = True
        estimator.options.environment.job_tags = ["TUT_HSVQE"]
        x_opt = spsa(
            cost_func,
            x0=x0,
            args=(isa_ansatz, isa_observable, estimator, cost_history_dict),
            maxiter=maxiter,
        )

        y_min = cost_func(
            x_opt, isa_ansatz, isa_observable, estimator, cost_history_dict
        )

    return y_min, cost_history_dict

np.random.seed(42)
num_params = ansatz.num_parameters
params = 2 * np.pi * np.random.random(num_params)

Aqui definimos o maxiter = 50. Observe que, como cada iteração requer duas chamadas à função para calcular o gradiente, o número total de chamadas à função será de $2 \times \text{maxiter}$. O valor de maxiter pode ser aumentado para qualquer valor maior, a fim de obter uma estimativa de energia mais precisa.

maxiter = 50
spsa_min, spsa_history = solve(
    params, isa_ansatz, isa_observable, maxiter=maxiter
)

Output:

Fx Iters. done: 101 [Current cost: -3.03843]

Etapa 4: Pós-processamento e retorno do resultado no formato clássico desejado

• Entrada: Estimativas da energia do estado fundamental durante a otimização
• Resultado: Energia estimada do estado fundamental

print(f"Estimated ground state energy: {spsa_min}")

Output:

Estimated ground state energy: [-3.03842968]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

Output:

Exemplo de hardware em grande escala

Este tutorial não inclui um exemplo de hardware em grande escala. À medida que o número de qubits aumenta, a VQE enfrenta desafios significativos devido ao fenômeno do platô estéril : o gradiente da função de custo desaparece exponencialmente com o aumento do tamanho do sistema, tornando a otimização praticamente inviável para circuitos de grande porte. Quando combinado com o ruído do hardware, isso significa que a aplicação do VQE a cadeias de spin maiores não produz resultados confiáveis e reproduzíveis. Para conhecer abordagens que superam essas limitações, consulte a seção “Próximos passos” abaixo.

Desafio

Agora que você tem uma implementação funcional do VQE para a cadeia de Heisenberg, tente o seguinte:

1. Experimente alterar a profundidade do ansatz: modifique o reps parâmetro em efficient_su2 (por exemplo, tente reps=1 e reps=3). Como a profundidade do ansatz afeta a estimativa da energia do estado fundamental e a velocidade de convergência? Em que momento você percebe uma diminuição do rendimento ou instabilidade?
2. Ajuste os hiperparâmetros do SPSA: ajuste os parâmetros da programação da taxa de aprendizagem (a, c, alpha, gamma, A) e observe como eles afetam a convergência. Você consegue encontrar uma configuração que converja mais rapidamente do que as predefinições usadas aqui?
3. Compare topologias de acoplamento: em vez de usar o mapa de acoplamento nativo do backend, tente construir uma cadeia linear simples de vizinhos mais próximos e compare os resultados. De que forma a conectividade do hardware físico afeta a profundidade do circuito transpilado e a estimativa final de energia?

Referências

[1] Spall, J. C. (2002). Implementação do algoritmo de perturbação simultânea para otimização estocástica. IEEE Transactions on Aerospace and Electronic Systems, 34(3), 817-823.

[2] Sahin, M. Emre, et al. (2025). Machine Learning do Qiskit: uma biblioteca de código aberto para tarefas de aprendizado de máquina quântico em grande escala em hardware quântico e simuladores clássicos. arXiv:2505.17756.

Próximas etapas