Simulare un modello di Ising forzato con la funzione TEM

Il metodo di mitigazione degli errori tramite rete tensoriale (TEM) di Algorithmiq è un algoritmo ibrido quantistico-classico progettato per eseguire la mitigazione del rumore interamente nella fase di post-elaborazione classica. Con TEM, puoi calcolare i valori di aspettativa degli osservabili, mitigando gli inevitabili errori indotti dal rumore sull'hardware quantistico con maggiore precisione ed efficienza dei costi, rendendolo un'opzione molto attraente per i ricercatori quantistici e i professionisti del settore.

Questo tutorial dimostra come TEM possa ottenere risultati significativi per la dinamica di un sistema quantistico, che sarebbe inaccessibile senza mitigazione degli errori e che richiederebbe sostanzialmente più risorse quantistiche se venissero utilizzati altri metodi di mitigazione degli errori come PEC e ZNE.

Stima di utilizzo: questo notebook utilizza circa 10 minuti QPU su dispositivi Heron r3. Il tempo di esecuzione può variare notevolmente in base al dispositivo scelto. Le stime di utilizzo per sezione sono riportate di seguito.

Eseguire esperimenti di fisica a molti corpi con mitigazione degli errori tramite la funzione TEM

Questo tutorial è basato sul seguente riferimento: L. E. Fischer et al., Nat. Phys. (2026). Questo riferimento discute una simulazione reale su hardware quantistico fino a 91 qubit. In questo tutorial, ricreamo una simulazione simile su un circuito di dimensioni più piccole.

Il modello di Ising forzato corrisponde al consueto modello di Ising:

$\hat{H}_{\text{I}} = J \sum_{n=0}^{N-2} \hat{Z}_n \hat{Z}_{n+1} + h \sum_{n=0}^{N-1} \hat{Z}_n$

al quale viene applicata una forzatura trasversale:

$\hat{H}_{K} = b \sum_{n=0}^{N-1} \hat{X}_n$

L'obiettivo è simulare la dinamica di uno stato sotto l'Hamiltoniano di Ising forzato trasversale, la cui evoluzione temporale può essere implementata da un unitario di Floquet $\hat{U}_{\text{KI}} = e^{-i \hat{H}_K} e^{-i \hat{H}_I}$. Lo stato iniziale da far evolvere è quello in cui il primo qubit è nello stato $|+\rangle$, mentre gli altri sono accoppiati e impostati nello stato di Bell $(|00\rangle + |11\rangle)/\sqrt{2}$.

La quantità che vogliamo osservare è la funzione di correlazione. L'articolo di riferimento discute come questa quantità possa essere riscritta come un operatore di Pauli $\hat{X}$ sull'$n^{esimo}$ qubit. Dopo un certo numero di passi temporali fisici $t$, calcoliamo il valore dell'operatore di Pauli $\hat{X}_{n=t}$. A seconda dei parametri del sistema, il valore di questo osservabile è uguale a un valore calcolabile esattamente, o solo simulabile mediante metodi approssimati. Specificamente, per $|J|=|b|=\pi/4$, è uguale a $[\cos(2h)]^t$, che è il valore che useremo per confrontare i risultati di questo tutorial. Inoltre, a un dato passo temporale $t$, $\langle\hat{X}_{n\neq t}\rangle$ è zero. Per i dettagli su come ottenere questi valori, e per un confronto con i risultati di simulazione classica approssimata al di fuori di questi parametri, vedere L. E. Fischer et al., Nat. Phys. (2026).

TEM funziona caratterizzando prima il rumore per ogni strato unico di gate a due qubit nel circuito, nonché caratterizzando l'errore di readout. Poi, il circuito viene eseguito sulla macchina quantistica. Infine, la mitigazione degli errori tramite rete tensoriale viene eseguita sulle risorse classiche in IBM Cloud® e viene restituito il valore mitigato. In questo esempio, il circuito ha due strati unici da caratterizzare.

Configurazione

Come prerequisito, assicurati che le dipendenze necessarie siano installate.

%pip install numpy matplotlib qiskit qiskit-ibm-catalog qiskit-ibm-runtime pylatexenc qiskit_qasm3_import

import os
from matplotlib import pyplot as plt
import numpy as np

from qiskit.quantum_info import SparsePauliOp
from qiskit.qasm3 import load

from qiskit_ibm_catalog import QiskitFunctionsCatalog

Mitigazione degli errori con TEM

Forniamo qui un circuito che implementa il modello di Ising forzato descritto sopra. Il circuito è preparato come segue. Prima c'è una fase di preparazione dello stato, in cui il primo qubit è nello stato $|+\rangle$, mentre gli altri sono in coppie di Bell $(|00\rangle + |11\rangle)/\sqrt{2}$. Questo è seguito dalla struttura a mattoni che implementa l'evoluzione unitaria $\hat{U}_{\text{KI}}$. Il numero di passi temporali fisici corrisponde a $t/2$ strati di circuito. Il seguente codice scarica i due file QASM necessari per questo tutorial.

# Download required QASM files
import urllib

urllib.request.urlretrieve(
    "https://ibm.box.com/shared/static/swy5jtq309b0xpzluzlmsmj908yphes8.qasm",
    "ki_30q.qasm",
)
urllib.request.urlretrieve(
    "https://ibm.box.com/shared/static/et3gkodonw6gsp2trs43lzaozrdtiu7s.qasm",
    "ki_12q.qasm",
)

Possiamo visualizzare una versione ridotta del circuito, con 12 qubit e sei passi temporali:

# Parameters of the kicked Ising model
h = 0.0
num_qubits = 12
t_steps = 6

# Load the circuit for the kicked Ising model
small_circuit = load("ki_12q.qasm")

# Draw the circuit
small_circuit.draw("mpl", scale=0.25, fold=-1)

Successivamente, costruisci l'osservabile $\hat{X}_{n=t}$. Viene costruito come una semplice stringa di Pauli con l'ordine che corrisponde a quello utilizzato da Qiskit:

def xt_observable(n_qubits, t_steps):
    pauli_str = "".join(["I" * t_steps, "X", "I" * (n_qubits - t_steps - 1)])
    pauli_str = pauli_str[::-1]  # Reverse the string to match qiskit order
    return SparsePauliOp(data=pauli_str, coeffs=1.0)

Nel nostro piccolo esempio a 12 qubit, l'osservabile appare così:

# Build the observable for the kicked Ising model
small_observable = xt_observable(n_qubits=12, t_steps=6)
print(small_observable)

SparsePauliOp(['IIIIIXIIIIII'],
              coeffs=[1.+0.j])

Le funzioni Qiskit utilizzano i PUB come mezzo per raccogliere gli input. Nel nostro caso, consideriamo un singolo circuito e osservabile come nostro PUB:

# Collect the input PUBs, in this case composed of a
# single circuit and observable
pubs = [(small_circuit, [small_observable])]

Successivamente, accediamo alla funzione TEM. Prima configuriamo l'autenticazione richiesta per IBM Cloud e selezioniamo un backend tra i dispositivi disponibili. Il token, i backend disponibili e i nomi delle risorse cloud corrispondenti (CRN) possono essere ottenuti accedendo al tuo account sul pannello di controllo di IBM Quantum Platform.

# Set IBM Quantum credentials and backend configuration
personal_token = os.environ.get(
    "QISKIT_IBM_TOKEN", "<API-KEY>"
)  # Replace with your personal token or set the environment variable
channel = "ibm_quantum_platform"
crn = "your_crn"  # Replace with the Cloud Resource Name (CRN)

# Select the QPU backend
backend_name = "ibm_qpu_name"  # Replace with your desired backend's name

Carica la funzione TEM dal Qiskit Functions Catalog:

# Load the TEM function from the Qiskit Functions Catalog
catalog = QiskitFunctionsCatalog(
    channel=channel,
    token=personal_token,
    instance=crn,
)
tem = catalog.load("algorithmiq/tem")

Possiamo ora eseguire un esperimento sul circuito di Ising forzato con la mitigazione degli errori fornita da TEM. Con le impostazioni predefinite, TEM può essere eseguito in modo semplice con un tempo di esecuzione QPU previsto di circa 2,5 minuti, a seconda del QPU:

tem_job = tem.run(pubs=pubs, backend_name=backend_name)

Con le opzioni predefinite, la funzione TEM esegue tre job sul computer quantistico: apprendimento del rumore, mitigazione del readout e campionamento del circuito. Il numero di shots utilizzati da ciascuno di questi può essere modificato nelle opzioni passate alla funzione. Per impostazione predefinita, questi parametri sono impostati per raggiungere una precisione di 0,05 nei valori di aspettativa mitigati. Puoi controllare lo stato del tuo job sul pannello di controllo di IBM Quantum Platform o con:

print(tem_job.status())

QUEUED

Quando lo stato è DONE, possiamo verificare i risultati grezzi e mitigati. I tem_evs definiti di seguito sono i valori di aspettativa degli osservabili richiesti, in questo caso un solo osservabile, $\langle \hat X_{n=t}\rangle$, e tem_std sono le corrispondenti deviazioni standard.

# Get the results of the TEM job
tem_results = tem_job.result()[
    0
]  # Get the first and only result from the job
tem_evs = tem_results.data.evs[0]
tem_std = tem_results.data.stds[0]

print(f"TEM Result: {tem_evs:.3f} ± {tem_std:.3f}")

TEM Result: 1.031 ± 0.046

Possiamo anche verificare quanti tempi di esecuzione quantistica sono stati utilizzati per ogni chiamata su IBM Quantum Platform, o ispezionando i metadati dei risultati dal codice Python.

# Get the TEM job runtime
tem_runtime = tem_job.result().metadata["resource_usage"][
    "RUNNING: EXECUTING_QPU"
]["QPU_TIME"]

print(f"TEM Runtime: {tem_runtime} seconds")

TEM Runtime: 155.0 seconds

Personalizzare i parametri TEM e le opzioni avanzate

La funzione TEM fornisce diverse opzioni avanzate per personalizzare il tuo workflow di mitigazione degli errori. Queste opzioni ti consentono di controllare la precisione, il numero di shots, le strategie di apprendimento del rumore e altri parametri per adattarsi meglio ai requisiti del tuo esperimento e alle risorse quantistiche disponibili.

Le opzioni avanzate comuni sono:

• precision: Specifica la precisione target per i valori di aspettativa mitigati.
• default_shots: Invece di precision, puoi specificare il numero di shots utilizzati dal job di misurazione.
• tem_max_bond_dimension: La dimensione di legame massima utilizzata nella rete tensoriale.
• tem_compression_cutoff: Il valore di cutoff da utilizzare per la rete tensoriale.
• opzioni di apprendimento del rumore: Configura come viene caratterizzato il rumore, come il numero di ripetizioni o specifici circuiti di calibrazione.
• private: Garantisce che i circuiti e i risultati dell'esperimento siano privati e disabilita i download multipli dei risultati del job.

Consulta la documentazione TEM o il Qiskit Functions Catalog per un elenco completo delle opzioni supportate e le loro descrizioni. Puoi regolare questi parametri per bilanciare il tempo di esecuzione, l'utilizzo delle risorse e la precisione dei risultati. Puoi passare queste opzioni come dizionario all'argomento options quando esegui la funzione TEM:

options = {
    "default_shots": 10_000,
    "tem_max_bond_dimension": 512,
    "tem_compression_cutoff": 1e-16,
    # This option helps optimizing the measurement
    # stage since the observable is strongly biased
    # toward the X operator for all the qubits.
    "compute_shadows_bias_from_observable": True,
    # set to True to keep experiment results private,
    # recommended for confidential circuits
    "private": False,
}

Possono essere passate anche opzioni personalizzate per il noise learner. Seguono le definizioni utilizzate nel Qiskit Runtime NoiseLearnerOptions:

nl_options = {
    "num_randomizations": 32,
    "max_layers_to_learn": 2,
    "shots_per_randomization": 128,
    "layer_pair_depths": [0, 1, 2, 4, 16, 32],
}

# add noise learning options to the overall options
options |= nl_options

Riesegui l'esperimento con queste opzioni personalizzate ottimizzate per il nostro circuito. Il tempo di esecuzione previsto è di circa quattro minuti QPU.

tem_job_custom = tem.run(
    pubs=pubs, backend_name=backend_name, options=options
)

Se il job non è impostato come privato, possiamo recuperare il risultato in un secondo momento. Per farlo, salva l'ID del job stampato qui e usa tem_job_custom = catalog.get_job_by_id("your-job-id").

job_id = tem_job_custom.job_id
print(f"Job ID: {job_id}")

Job ID: 1ba10094-a541-457a-9287-dbd49306d12d

results_custom = tem_job_custom.result()
tem_evs = results_custom[0].data.evs[0]
tem_std = results_custom[0].data.stds[0]

print(f"TEM Result: {tem_evs:.3f} ± {tem_std:.3f}")

TEM Result: 0.956 ± 0.018

Possiamo ora ispezionare i risultati e i metadati per ottenere informazioni sull'esperimento:

metadata_custom = results_custom[0].metadata

unmitigated_evs = metadata_custom["evs_non_mitigated"][0]
unmitigated_stds = metadata_custom["stds_non_mitigated"][0]
print(f"Unmitigated Result: {unmitigated_evs:.3f} ± {unmitigated_stds:.3f}")

# Exact result for the kicked Ising model from the reference paper
exact_evs = np.cos(2 * h) ** t_steps
print("Exact Result:", exact_evs)

Unmitigated Result: 0.894 ± 0.015
Exact Result: 1.0

# Plot comparing the different expectation values
plt.bar(
    ["Unmitigated", "TEM"],
    [unmitigated_evs, tem_evs],
    yerr=[unmitigated_stds, tem_std],
    color=["grey", "c"],
)
plt.hlines(y=exact_evs, xmin=-0.5, xmax=1.5, colors="r", linestyles="dashed")
plt.ylabel("Expectation Value")
plt.ylim(0, 1.1)
plt.show()

Infine, possiamo verificare l'impatto delle opzioni personalizzate sul tempo di esecuzione QPU e classico:

# Get the metadata of the TEM job
job_metadata = results_custom.metadata

# Get the runtime of the TEM job
qpu_runtime = job_metadata["resource_usage"]["RUNNING: EXECUTING_QPU"][
    "QPU_TIME"
]
classical_runtime = (
    job_metadata["resource_usage"]["RUNNING: OPTIMIZING_FOR_HARDWARE"][
        "CPU_TIME"
    ]
    + job_metadata["resource_usage"]["RUNNING: POST_PROCESSING"]["CPU_TIME"]
)

print(f"QPU Runtime: {qpu_runtime} seconds")
print(f"Classical Runtime: {classical_runtime} seconds")

QPU Runtime: 342.0 seconds
Classical Runtime: 107.632604 seconds

Scalare TEM a circuiti grandi

I circuiti grandi possono, in linea di principio, essere eseguiti con la funzione TEM. Tuttavia, è importante essere consapevoli delle limitazioni delle risorse classiche, poiché TEM viene eseguito su runner IBM Cloud con tempi di esecuzione potenzialmente molto lunghi. Per circuiti estremamente grandi, contatta il team di supporto TEM all'indirizzo qiskit_ibm@algorithmiq.fi.

Qui eseguiamo un esempio con un circuito più grande a scala utilitaristica di 30 qubit, ottimizzando i parametri TEM per la velocità piuttosto che per la precisione.

# Kicked Ising model parameters
n_qubits = 30
t_steps = 15
h = 0.0

# Load the circuit for the kicked Ising model
circuit = load("ki_30q.qasm")

# Build the observable for the kicked Ising model
observable = xt_observable(n_qubits=n_qubits, t_steps=t_steps)

# Collect the input PUBs, in this case composed of a
# single circuit and observable
pubs = [(circuit, [observable])]

Definiamo alcune opzioni orientate alle prestazioni:

options = {
    "num_randomizations": 32,
    "max_layers_to_learn": 2,
    "shots_per_randomization": 128,
    "layer_pair_depths": [0, 1, 2, 4, 16, 32, 64],
    "default_shots": 5_000,
    "tem_max_bond_dimension": 128,
    "tem_compression_cutoff": 1e-10,
    "compute_shadows_bias_from_observable": True,
    "private": False,
}

Infine, esegui l'esperimento, ottieni il risultato e visualizzalo. Questo richiederà circa 3,5 minuti QPU.

tem_job_large = tem.run(pubs=pubs, backend_name=backend_name, options=options)

job_id = tem_job_large.job_id
print(f"Job ID: {job_id}")

Job ID: 9f3f190f-f4b0-4dcb-bb83-5f71f37d0d77

results_large = tem_job_large.result()
tem_evs = results_large[0].data.evs[0]
tem_std = results_large[0].data.stds[0]

print(f"TEM Result: {tem_evs:.3f} ± {tem_std:.3f}")

# Get the metadata of the TEM job
job_metadata = tem_job_large.result().metadata

# Get the runtime of the TEM job
qpu_runtime = job_metadata["resource_usage"]["RUNNING: EXECUTING_QPU"][
    "QPU_TIME"
]
classical_runtime = (
    job_metadata["resource_usage"]["RUNNING: OPTIMIZING_FOR_HARDWARE"][
        "CPU_TIME"
    ]
    + job_metadata["resource_usage"]["RUNNING: POST_PROCESSING"]["CPU_TIME"]
)

print(f"QPU Runtime: {qpu_runtime} seconds")
print(f"Classical Runtime: {classical_runtime} seconds")

TEM Result: 0.794 ± 0.026
QPU Runtime: 203.0 seconds
Classical Runtime: 251.71805499999996 seconds

# Plot comparing the different expectation values
metadata_large = results_large[0].metadata
unmitigated_evs = metadata_large["evs_non_mitigated"][0]
unmitigated_stds = metadata_large["stds_non_mitigated"][0]

exact_evs = np.cos(2 * h) ** t_steps

plt.bar(
    ["Unmitigated", "TEM"],
    [unmitigated_evs, tem_evs],
    yerr=[unmitigated_stds, tem_std],
    color=["grey", "c"],
)
plt.hlines(y=exact_evs, xmin=-0.5, xmax=1.5, colors="r", linestyles="dashed")
plt.ylabel("Expectation Value")
plt.ylim(0, 1.1)
plt.show()