Reduzindo o erro de Trotter da dinâmica Hamiltoniana com fórmulas multi-produto

Neste notebook, você aprenderá a usar uma Fórmula Multi-Produto (MPF) para alcançar um erro de Trotter menor em nosso observável em comparação com aquele incorrido pelo circuito de Trotter mais profundo que de fato executaremos. Você fará isso percorrendo as etapas de um Qiskit pattern:

• Etapa 1: Mapear para o problema quântico
  • Inicializar o Hamiltoniano do nosso problema
  • Usar uma MPF para gerar os circuitos de evolução temporal Trotterizados
• Etapa 2: Otimizar o problema
  • Aqui, transpilamos nossos circuitos para um GenericBackendV2
• Etapa 3: Executar experimentos
  • Usar um StatevectorEstimator por simplicidade neste notebook
• Etapa 4: Reconstruir resultados
  • Calcular o valor esperado da MPF

Etapa 1: Mapear para o problema quântico

1a: Configurando nosso Hamiltoniano

Usamos o modelo de Ising em uma linha de 10 sítios:

$$\hat{\mathcal{H}}_{\text{Ising}} = \sum_{i=1}^{9} J_{i,(i+1)} Z_i Z_{(i+1)} + \sum_{i=1}^{10} h_i X_i \, ,$$

onde $J$ é a intensidade do acoplamento entre dois sítios e $h$ é o campo magnético externo. O pacote qiskit_addon_utils fornece algumas funcionalidades reutilizáveis para diversos propósitos.

Seu módulo qiskit_addon_utils.problem_generators fornece funções para gerar Hamiltonianos do tipo Heisenberg em um determinado grafo de conectividade. Esse grafo pode ser tanto um rustworkx.PyGraph quanto um CouplingMap, facilitando o uso em fluxos de trabalho centrados em Qiskit.

A seguir, criamos uma linha simples de 10 qubits usando o método CouplingMap.from_line.

# Added by doQumentation — required packages for this notebook
!pip install -q numpy qiskit qiskit-addon-mpf qiskit-addon-utils rustworkx scipy

from qiskit.transpiler import CouplingMap

# Generate some coupling map to use for this example
coupling_map = CouplingMap.from_line(10, bidirectional=False)

from rustworkx.visualization import graphviz_draw

graphviz_draw(coupling_map.graph, method="circo")

Em seguida, geramos o SparsePauliOp na conectividade fornecida com as constantes desejadas.

from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian

# Get a qubit operator describing the Ising field model
hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(0.0, 0.0, 1.0),
    ext_magnetic_field=(0.4, 0.0, 0.0),
)
print(hamiltonian)

SparsePauliOp(['IIIIIIIZZI', 'IIIIIZZIII', 'IIIZZIIIII', 'IZZIIIIIII', 'IIIIIIIIZZ', 'IIIIIIZZII', 'IIIIZZIIII', 'IIZZIIIIII', 'ZZIIIIIIII', 'IIIIIIIIIX', 'IIIIIIIIXI', 'IIIIIIIXII', 'IIIIIIXIII', 'IIIIIXIIII', 'IIIIXIIIII', 'IIIXIIIIII', 'IIXIIIIIII', 'IXIIIIIIII', 'XIIIIIIIII'],
              coeffs=[1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j, 1. +0.j,
 1. +0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j,
 0.4+0.j, 0.4+0.j, 0.4+0.j])

O observável que mediremos é a magnetização total, que podemos construir simplesmente como mostrado abaixo:

from qiskit.quantum_info import SparsePauliOp

L = coupling_map.size()
observable = SparsePauliOp.from_sparse_list([("Z", [i], 1 / L / 2) for i in range(L)], num_qubits=L)
print(observable)

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

1b: Fórmulas Multi-Produto

As MPFs reduzem o erro de Trotter da dinâmica Hamiltoniana por meio de uma combinação ponderada de várias execuções de circuito.

Para tornar isso mais concreto, definimos uma MPF como:

$$\mu(t) = \sum_{j} x_j \rho^{k_j}_{j}\left(\frac{t}{k_j}\right) + \text{some remaining Trotter error} \, ,$$

onde $x_j$ são nossos coeficientes de ponderação, $\rho^{k_j}_j$ é a matriz densidade correspondente ao estado puro obtido pela evolução do estado inicial com a fórmula de produto, $S^{k_j}$, envolvendo $k_j$ etapas de Trotter, e $j$ indexa o número de fórmulas de produto que compõem a MPF.

A chave aqui é que o erro de Trotter remanescente é menor do que o erro de Trotter que se obteria simplesmente usando o maior valor de $k_j$!

Você pode encarar a utilidade da MPF de duas perspectivas:

1. Para um orçamento fixo de etapas de Trotter que você é capaz de executar, você pode obter resultados com um erro de Trotter total menor.
2. Para um número de etapas de Trotter que resultam em circuitos profundos, você pode usar a MPF para encontrar vários circuitos de menor profundidade para executar que resultam em um erro de Trotter semelhante.

Uma introdução às MPFs estáticas

MPFs estáticas são aquelas em que os valores de $x_j$ NÃO dependem do tempo de evolução, $t$.

Determinar os coeficientes da MPF estática para um dado conjunto de valores de $k_j$ equivale a resolver um sistema linear de equações: $Ax=b$, onde $x$ são nossos coeficientes de interesse, $A$ é uma matriz que depende de $k_j$ e do tipo de PF que usamos ($S$), e $b$ é um vetor de restrições. Por brevidade, não vamos entrar em mais detalhes aqui e, em vez disso, encaminhamos você para a documentação de LSE.

Podemos encontrar uma solução para $x$ analiticamente como $x = A^{-1}b$, ver, por exemplo, Carrera Vazquez et al., 2023 ou Zhuk et al., 2023. No entanto, essa solução exata pode ser "mal condicionada", resultando em normas L1 muito grandes de nossos coeficientes, $x$, o que pode levar a um mau desempenho da MPF. Em vez disso, também é possível obter uma solução aproximada que minimize a norma L1 de $x$ a fim de tentar otimizar o comportamento da MPF.

A seguir, você aprenderá a fazer tudo isso.

Escolhendo $k_j$

A escolha de $k_j$ cabe ao usuário final. Em princípio, quaisquer valores podem ser escolhidos, mas alguns $k_j$ levarão a uma maior amplificação de ruído em dispositivos reais do que outras escolhas. Assim, é importante tentar encontrar valores "bons" de $k_j$.

Aqui, simplesmente escolheremos alguns valores fixos para $k_j$. O menor valor é motivado pelo tempo de evolução-alvo de $t=8.0$, que normalmente nos diz para satisfazer $t/k_{\text{min}} \lt 1$, mas empiricamente sabemos que defini-lo igual a $1$ geralmente também funciona. Se quiser saber mais sobre isso e como escolher seus outros valores de $k_j$, consulte o respectivo guia: How to choose the Trotter steps for an MPF.

time = 8.0
trotter_steps = (8, 12, 19)

Configurando o LSE

Agora que escolhemos nossos $k_j$s, devemos primeiro construir o LSE, $Ax=b$, conforme explicado acima. A matriz $A$ depende não apenas de $k_j$, mas também da nossa escolha de fórmula de produto (PF) — em particular, de sua ordem. Adicionalmente, pode-se levar em conta se a PF é simétrica ou não (ver Carrera Vazquez et al., 2023), definindo symmetric=True. No entanto, isso não é obrigatório, conforme demonstrado por Zhuk et al., 2023.

Aqui, vamos usar uma fórmula Suzuki-Trotter de segunda ordem, resultando em order=2, e definiremos symmetric=True.

from qiskit_addon_mpf.static import setup_static_lse

lse = setup_static_lse(trotter_steps, order=2, symmetric=True)
print(lse)

LSE(A=array([[1.00000000e+00, 1.00000000e+00, 1.00000000e+00],
       [1.56250000e-02, 6.94444444e-03, 2.77008310e-03],
       [2.44140625e-04, 4.82253086e-05, 7.67336039e-06]]), b=array([1., 0., 0.]))

Resolvendo $x$ analiticamente

Como mencionado anteriormente, podemos encontrar $x$ analiticamente:

import numpy as np

coeffs_analytical = lse.solve()
print(coeffs_analytical)

[ 0.17239057 -1.19447005  2.02207947]

Otimizando $x$ usando um modelo exato

Como alternativa a calcular $x=A^{-1}b$, você também pode usar setup_exact_problem para construir uma instância de cvxpy.Problem que utiliza o LSE como restrições e cuja solução ótima produzirá $x$.

Na próxima seção, ficará claro por que essa interface existe.

from qiskit_addon_mpf.costs import setup_exact_problem

model_exact, coeffs_exact = setup_exact_problem(lse)
model_exact.solve()
print(coeffs_exact.value)

[ 0.17239057 -1.19447005  2.02207947]

Como indicador de se uma MPF construída com esses coeficientes produzirá bons resultados, podemos usar a norma L1 (ver também Carrera Vazquez et al., 2023).

print(np.linalg.norm(coeffs_exact.value, ord=1))

3.3889400921655914

Otimizando $x$ usando um modelo aproximado

Pode acontecer de a norma L1 para o conjunto escolhido de valores de $k_j$ ser considerada muito alta. Se for esse o caso e você não puder escolher um conjunto diferente de valores de $k_j$, pode usar uma solução aproximada para o LSE em vez de uma exata.

Para isso, basta usar setup_sum_of_squares_problem para construir uma instância diferente de cvxpy.Problem que restringe a norma L1 a um limite escolhido enquanto minimiza a diferença entre $Ax$ e $b$.

from qiskit_addon_mpf.costs import setup_sum_of_squares_problem

model_approx, coeffs_approx = setup_sum_of_squares_problem(lse, max_l1_norm=3.0)
model_approx.solve()
print(coeffs_approx.value)
print(np.linalg.norm(coeffs_approx.value, ord=1))

[-0.40454257  0.57553173  0.8290123 ]
1.8090865903790838

Observe que você tem total liberdade sobre como resolver esse problema de otimização, ou seja, pode alterar o solver de otimização, seus limiares de convergência, e assim por diante. Confira o respectivo guia em How to use the approximate model.

1c: Configurando os circuitos de Trotter

Neste ponto, encontramos nossos coeficientes de expansão, $x$, e tudo o que resta a fazer é gerar os circuitos quânticos Trotterizados. Mais uma vez, o módulo qiskit_addon_utils.problem_generators vem ao resgate para fazer exatamente isso:

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import generate_time_evolution_circuit

circuits = []
for k in trotter_steps:
    circ = generate_time_evolution_circuit(
        hamiltonian,
        synthesis=SuzukiTrotter(order=2, reps=k),
        time=time,
    )
    circuits.append(circ)

circuits[0].draw("mpl", fold=-1)

circuits[1].draw("mpl", fold=-1)

circuits[2].draw("mpl", fold=-1)

Etapa 2: Otimizar o problema

Normalmente, esta é a etapa do pattern durante a qual você otimiza seus circuitos para execução em hardware. Aqui, como usamos apenas um simulador sem ruído, simplesmente transpilamos nosso circuito para um GenericBackendV2.

from qiskit.providers.fake_provider import GenericBackendV2
from qiskit.transpiler import generate_preset_pass_manager

backend = GenericBackendV2(num_qubits=10)
transpiler = generate_preset_pass_manager(optimization_level=2, backend=backend)

transpiled_circuits = [transpiler.run(circ) for circ in circuits]

Etapa 3: Executar experimentos quânticos

Como explicado no início, vamos pular a etapa 2 de otimização porque iremos simplesmente calcular os valores esperados do nosso observável-alvo usando um simulador livre de ruído, a saber, o StatevectorEstimator.

from qiskit.primitives import StatevectorEstimator

estimator = StatevectorEstimator()
job = estimator.run([(circ, observable) for circ in transpiled_circuits])
result = job.result()

Etapa 4: Reconstruir resultados

Primeiro, extraímos os valores esperados individuais obtidos para cada um dos circuitos de Trotter:

evs = [res.data.evs for res in result]
print(evs)

[array(0.23799162), array(0.35754312), array(0.38649906)]

Em seguida, simplesmente os recombinamos com nossos coeficientes da MPF para produzir os valores esperados totais da MPF. Abaixo, fazemos isso para cada uma das diferentes formas pelas quais calculamos $x$.

print("Analytical    solution:", evs @ coeffs_analytical)
print("Exact model   solution:", evs @ coeffs_exact.value)
print("Approx. model solution:", evs @ coeffs_approx.value)

Analytical    solution: 0.3954847855980006
Exact model   solution: 0.39548478559800204
Approx. model solution: 0.42991214253489807

Por fim, para este pequeno problema, podemos calcular o valor de referência exato usando scipy.linalg.expm da seguinte forma:

from scipy.linalg import expm

exp_H = expm(-1j * time * hamiltonian.to_matrix())

initial_state = np.zeros(exp_H.shape[0])
initial_state[0] = 1.0

time_evolved_state = exp_H @ initial_state

exact_obs = time_evolved_state.conj() @ observable.to_matrix() @ time_evolved_state
print(exact_obs.real)

0.40060242487899755

Podemos ver claramente que a MPF reduziu o erro de Trotter em comparação com aquele obtido com a PF individual mais profunda, com $k_j=19$. No entanto, também vemos que o modelo aproximado não é perfeito, pois, na verdade, resultou em um valor esperado pior do que a solução exata. Isso mostra a importância de usar critérios de convergência rigorosos no modelo aproximado, conforme você aprenderá no guia How to use the approximate model.