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

Estimativa de uso: 37 minutos em um processador Heron (NOTA: Esta é apenas uma estimativa. Seu tempo de execução pode variar.)

Resultados de aprendizagem

Após concluir este tutorial, você pode esperar entender as seguintes informações:

• Como modelar uma cadeia de spins 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 hardware quântico IBM® usando primitivas e sessões do Qiskit Runtime

Pré-requisitos

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

• Fundamentos de informação quântica
• Introdução aos padrões Qiskit
• Design de algoritmos variacionais

Contexto

A cadeia de spins de Heisenberg é um dos modelos mais amplamente estudados na física da matéria condensada e no magnetismo quântico. Ela descreve uma rede unidimensional de spins quânticos em interação, onde spins de primeiros vizinhos são acoplados por meio de interações de troca. O Hamiltoniano para o 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 atuando no sítio $i$, a soma $\langle i,j \rangle$ percorre os pares de primeiros vizinhos, $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 do intervalo $[-1, 1]$. Note que na implementação abaixo, o conjunto de pares de "primeiros vizinhos" é determinado pelo acoplamento nativo do backend de hardware entre os primeiros $N$ qubits, o que pode não formar uma cadeia linear estrita dependendo da topologia do dispositivo.

Entender a energia do estado fundamental deste Hamiltoniano é de fundamental importância na física. O estado fundamental codifica informações sobre transições de fase quântica, estrutura de emaranhamento e ordenamento magnético. Classicamente, calcular a energia exata do estado fundamental torna-se intratável à medida que o número de spins cresce, pois a dimensão do espaço de Hilbert escala exponencialmente como $2^N$ para $N$ spins. Isso o torna um candidato natural para 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. Ele funciona preparando um estado quântico parametrizado $|\psi(\theta)\rangle$ (chamado ansatz) em um computador quântico e medindo o valor esperado $\langle \psi(\theta) | H | \psi(\theta) \rangle$. Um otimizador clássico então ajusta iterativamente os parâmetros $\theta$ para minimizar essa energia, aproveitando o princípio variacional, que garante que a energia medida é sempre um limite superior para a verdadeira energia do estado fundamental.

Neste tutorial, usamos o ansatz efficient_su2 da biblioteca de circuitos do Qiskit, que constrói camadas de rotações de qubit único e portas de emaranhamento. A otimização é realizada usando o algoritmo de Aproximação Estocástica por Perturbação Simultânea (SPSA), que é bem adequado para hardware quântico ruidoso porque estima gradientes usando apenas duas avaliações de função por iteração, independentemente do número de parâmetros.

Requisitos

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

• Qiskit SDK v2.0 ou posterior, com suporte para visualização
• Qiskit Runtime v0.44 ou posterior (pip install qiskit-ibm-runtime)

Configuração

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

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, percorremos cada etapa do padrão Qiskit em pequena escala, explicando os componentes principais conforme construímos o fluxo de trabalho.

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

• Entrada: Número de spins
• Saída: Ansatz e Hamiltoniano modelando a cadeia de Heisenberg

Construa um ansatz e Hamiltoniano que modelam uma cadeia de Heisenberg de 10 spins. Nesta etapa, vamos construir um Hamiltoniano de Heisenberg de 10 spins sobre o mapa de acoplamento do backend menos ocupado e preparar o ansatz efficient_su2.

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

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

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

Use a função generate_preset_pass_manager do Qiskit para gerar automaticamente uma rotina de otimização para nosso circuito em relação à QPU selecionada. Escolhemos optimization_level=3, que fornece o mais alto nível de otimização dos gerenciadores de passagem predefinidos. Também incluímos passes de agendamento ALAPScheduleAnalysis e PadDynamicalDecoupling 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)

Passo 3: Executar usando primitivas Qiskit

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

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

Como otimizamos o circuito para o backend no Passo 2, podemos evitar a transpilação no servidor Runtime definindo skip_transpilation=True e passando o circuito otimizado. Para esta demonstração, executaremos em uma QPU usando primitivas qiskit-ibm-runtime. Para executar com primitivas baseadas em statevector do qiskit, substitua o bloco de código que usa primitivas do Qiskit Runtime pelo bloco comentado.

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

Introduzindo o SPSA

A Aproximação Estocástica por Perturbação Simultânea (SPSA) [1] é um algoritmo de otimização que aproxima o vetor gradiente completo usando 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, cria-se um vetor aleatório $\Delta_i$ de tamanho $p$, onde cada elemento $\Delta_{ij}$, $\forall$ $j\in \{1,2,...,p\}$, é amostrado uniformemente de $\{-1, 1\}$. Em seguida, cada elemento do vetor aleatório $\Delta_i$ é multiplicado por um valor pequeno $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, como uma perturbação aleatória é aplicada durante a estimativa do gradiente, espera-se que pequenos desvios nos valores exatos de $f$ provenientes de ruído possam ser tolerados e contabilizados. De fato, o SPSA é particularmente conhecido por ser robusto contra ruído e requer apenas duas chamadas de hardware por iteração. Portanto, é um dos otimizadores altamente preferidos para implementar algoritmos variacionais.

Neste tutorial, os hiperparâmetros para a $i^{th}$ iteração, $a_i$ e $c_i$, são calculados como

$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 $A = 30$, $\alpha = 0.9$, $a = 0.3$, $c = 0.1$ e $\gamma = 0.4$. Esses valores são selecionados de [2]. O ajuste apropriado dos hiperparâmetros é necessário para extrair 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 maxiter = 50. Note que como cada iteração requer duas chamadas à função para calcular o gradiente, o número total de chamadas de função será $2 \times \text{maxiter}$. O maxiter pode ser aumentado para qualquer valor mais alto para uma melhor estimativa de energia.

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

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

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

• Entrada: Estimativas de energia de estado fundamental durante a otimização
• Saída: Energia de estado fundamental estimada

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

Estimated ground state energy: [-3.03842968]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

Exemplo em hardware de grande escala

Um exemplo em hardware de grande escala não está incluído neste tutorial. À medida que o número de qubits aumenta, o VQE encontra desafios significativos devido ao fenômeno do platô estéril: o gradiente da função de custo desaparece exponencialmente com o tamanho do sistema, tornando a otimização praticamente inviável para circuitos grandes. Combinado com o ruído do hardware, isso significa que escalar o VQE para cadeias de spins maiores não produz resultados reprodutíveis de forma confiável. Para abordagens que superam essas limitações, consulte a seção de Próximos passos abaixo.

Desafio

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

1. Experimente com a profundidade do ansatz: Modifique o parâmetro reps no efficient_su2 (por exemplo, tente reps=1 e reps=3). Como a profundidade do ansatz afeta a energia estimada do estado fundamental e a velocidade de convergência? Em que ponto você observa retornos decrescentes ou instabilidade?
2. Ajuste os hiperparâmetros do SPSA: Ajuste os parâmetros da programação de taxa de aprendizado (a, c, alpha, gamma, A) e observe como eles impactam a convergência. Você consegue encontrar uma configuração que converge mais rapidamente do que os padrões usados aqui?
3. Compare topologias de acoplamento: Em vez de usar o mapa de acoplamento nativo do backend, tente construir uma cadeia linear simples de primeiros vizinhos e compare os resultados. Como 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). Implementation of the simultaneous perturbation algorithm for stochastic optimization. IEEE Transactions on Aerospace and Electronic Systems, 34(3), 817-823.

[2] Sahin, M. Emre, et al. (2025). Qiskit Machine Learning: an open-source library for quantum machine learning tasks at scale on quantum hardware and classical simulators. arXiv:2505.17756.

Próximos passos

Recomendações

Se você achou este trabalho interessante, pode se interessar pelo seguinte material:

• Experimente a Diagonalização Quântica Baseada em Amostras (SQD): Como demonstrado neste tutorial, o VQE enfrenta desafios em escala devido a platôs estéreis e alta sobrecarga de medição. A IBM desenvolveu a Diagonalização Quântica Baseada em Amostras (SQD) como uma alternativa mais escalável. Ao contrário do VQE, o SQD evita totalmente a otimização variacional; em vez disso, um computador quântico gera amostras e um computador clássico projeta o Hamiltoniano em um subespaço gerado por essas amostras e o diagonaliza. Isso fornece um limite superior para a energia do estado fundamental com significativamente menos medições e sem suscetibilidade a platôs estéreis. Siga o tutorial de SQD para ver essa abordagem em ação.
• Explore o curso de Algoritmos de Diagonalização Quântica: Aprofunde seu entendimento sobre VQE e SQD, incluindo suas vantagens e desvantagens, no curso Algoritmos de diagonalização quântica no IBM Quantum Learning.