Diagonalizzazione quantistica di Krylov degli hamiltoniani reticolari

Stima di utilizzo: 20 minuti su un Heron r2 (NOTA: questa è solo una stima. Il tempo di esecuzione potrebbe variare)

Sfondo

Questa esercitazione mostra come implementare l'algoritmo di diagonalizzazione quantistica di Krylov (KQD) nel contesto dei modelli Qiskit. Verrà dapprima illustrata la teoria alla base dell'algoritmo e poi verrà mostrata una dimostrazione della sua esecuzione su una QPU.

In tutte le discipline, siamo interessati a conoscere le proprietà dello stato fondamentale dei sistemi quantistici. Tra gli esempi vi sono la comprensione della natura fondamentale delle particelle e delle forze, la previsione e la comprensione del comportamento di materiali complessi e la comprensione delle interazioni e delle reazioni biochimiche. A causa della crescita esponenziale dello spazio di Hilbert e delle correlazioni che sorgono nei sistemi entangled, gli algoritmi classici faticano a risolvere questo problema per sistemi quantistici di dimensioni crescenti. A un'estremità dello spettro c'è l'approccio esistente che sfrutta l'hardware quantistico concentrandosi sui metodi quantistici variazionali (ad esempio, l' eigensolver quantistico variazionale ). Queste tecniche si scontrano con i dispositivi attuali a causa dell'elevato numero di chiamate di funzione richieste nel processo di ottimizzazione, che aggiungono un grande overhead di risorse una volta introdotte le tecniche avanzate di mitigazione degli errori, limitando così la loro efficacia ai sistemi di piccole dimensioni. All'altra estremità dello spettro, ci sono i metodi quantistici a tolleranza di errore con garanzie di prestazioni (ad esempio, la stima quantistica della fase ), che richiedono circuiti profondi che possono essere eseguiti solo su un dispositivo a tolleranza di errore. Per questi motivi, introduciamo qui un algoritmo quantistico basato su metodi subspaziali (come descritto in questo articolo ), l'algoritmo di diagonalizzazione quantistica di Krylov (KQD). Questo algoritmo si comporta bene su larga scala [1] sull'hardware quantistico esistente, condivide garanzie di prestazione simili a quelle della stima di fase, è compatibile con tecniche avanzate di mitigazione degli errori e potrebbe fornire risultati classicamente inaccessibili.

Requisiti

Prima di iniziare questa esercitazione, assicuratevi di aver installato quanto segue:

• Qiskit SDK v2.0 o versioni successive, con supporto alla visualizzazione
• Qiskit Runtime v0.22 o successivamente ( pip install qiskit-ibm-runtime )

Configura

import numpy as np
import scipy as sp
import matplotlib.pylab as plt
from typing import Union, List
import itertools as it
import copy
from sympy import Matrix
import warnings

warnings.filterwarnings("ignore")

from qiskit.quantum_info import SparsePauliOp, Pauli, StabilizerState
from qiskit.circuit import Parameter, IfElseOp
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter
from qiskit.transpiler import Target, CouplingMap
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager


from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    EstimatorV2 as Estimator,
)


def solve_regularized_gen_eig(
    h: np.ndarray,
    s: np.ndarray,
    threshold: float,
    k: int = 1,
    return_dimn: bool = False,
) -> Union[float, List[float]]:
    """
    Method for solving the generalized eigenvalue problem with regularization

    Args:
        h (numpy.ndarray):
            The effective representation of the matrix in the Krylov subspace
        s (numpy.ndarray):
            The matrix of overlaps between vectors of the Krylov subspace
        threshold (float):
            Cut-off value for the eigenvalue of s
        k (int):
            Number of eigenvalues to return
        return_dimn (bool):
            Whether to return the size of the regularized subspace

    Returns:
        lowest k-eigenvalue(s) that are the solution of the
        regularized generalized eigenvalue problem


    """
    s_vals, s_vecs = sp.linalg.eigh(s)
    s_vecs = s_vecs.T
    good_vecs = np.array(
        [vec for val, vec in zip(s_vals, s_vecs) if val > threshold]
    )
    h_reg = good_vecs.conj() @ h @ good_vecs.T
    s_reg = good_vecs.conj() @ s @ good_vecs.T
    if k == 1:
        if return_dimn:
            return sp.linalg.eigh(h_reg, s_reg)[0][0], len(good_vecs)
        else:
            return sp.linalg.eigh(h_reg, s_reg)[0][0]
    else:
        if return_dimn:
            return sp.linalg.eigh(h_reg, s_reg)[0][:k], len(good_vecs)
        else:
            return sp.linalg.eigh(h_reg, s_reg)[0][:k]


def single_particle_gs(H_op, n_qubits):
    """
    Find the ground state of the single particle(excitation) sector
    """
    H_x = []
    for p, coeff in H_op.to_list():
        H_x.append(set([i for i, v in enumerate(Pauli(p).x) if v]))

    H_z = []
    for p, coeff in H_op.to_list():
        H_z.append(set([i for i, v in enumerate(Pauli(p).z) if v]))

    H_c = H_op.coeffs

    print("n_sys_qubits", n_qubits)

    n_exc = 1
    sub_dimn = int(sp.special.comb(n_qubits + 1, n_exc))
    print("n_exc", n_exc, ", subspace dimension", sub_dimn)

    few_particle_H = np.zeros((sub_dimn, sub_dimn), dtype=complex)

    # list all of the possible sets of n_exc indices of 1s in
    # n_exc-particle states
    sparse_vecs = [
        set(vec) for vec in it.combinations(range(n_qubits + 1), r=n_exc)
    ]

    m = 0
    for i, i_set in enumerate(sparse_vecs):
        for j, j_set in enumerate(sparse_vecs):
            m += 1

            if len(i_set.symmetric_difference(j_set)) <= 2:
                for p_x, p_z, coeff in zip(H_x, H_z, H_c):
                    if i_set.symmetric_difference(j_set) == p_x:
                        sgn = ((-1j) ** len(p_x.intersection(p_z))) * (
                            (-1) ** len(i_set.intersection(p_z))
                        )
                    else:
                        sgn = 0

                    few_particle_H[i, j] += sgn * coeff

    gs_en = min(np.linalg.eigvalsh(few_particle_H))
    print("single particle ground state energy: ", gs_en)
    return gs_en

Fase 1: mappare gli input classici su un problema quantistico

Lo spazio di Krylov

Lo spazio di Krylov $\mathcal{K}^r$ di ordine $r$ è lo spazio spaziato dai vettori ottenuti moltiplicando le potenze superiori di una matrice $A$, fino a $r-1$, con un vettore di riferimento $\vert v \rangle$.

Se la matrice $A$ è l'hamiltoniana $H$, ci riferiremo allo spazio corrispondente come spazio di Krylov di potenza $\mathcal{K}_P$. Nel caso in cui $A$ sia l'operatore di evoluzione temporale generato dall'hamiltoniana $U=e^{-iHt}$, ci riferiremo allo spazio come spazio di Krylov unitario $\mathcal{K}_U$. Il sottospazio di Krylov di potenza che utilizziamo classicamente non può essere generato direttamente su un computer quantistico, poiché $H$ non è un operatore unitario. Si può invece utilizzare l'operatore di evoluzione temporale $U = e^{-iHt}$, che può dare garanzie di convergenza simili a quelle del metodo delle potenze. Le potenze di $U$ diventano quindi passi temporali diversi $U^k = e^{-iH(kt)}$.

Si veda l'Appendice per una derivazione dettagliata di come lo spazio di Krylov unitario permetta di rappresentare accuratamente gli autostati a bassa energia.

Algoritmo di diagonalizzazione quantistica di Krylov

Data un'hamiltoniana $H$ che vogliamo diagonalizzare, consideriamo innanzitutto il corrispondente spazio di Krylov unitario $\mathcal{K}_U$. L'obiettivo è trovare una rappresentazione compatta dell'hamiltoniana in $\mathcal{K}_U$, che chiameremo $\tilde{H}$. Gli elementi della matrice di $\tilde{H}$, la proiezione dell'hamiltoniana nello spazio di Krylov, possono essere calcolati calcolando i seguenti valori di aspettativa

Dove $\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle$ sono i vettori dello spazio di Krylov unitario e $t_n = n dt$ sono i multipli del passo temporale $dt$ scelti. Su un computer quantistico, il calcolo degli elementi di ciascuna matrice può essere effettuato con qualsiasi algoritmo che consenta di ottenere la sovrapposizione tra stati quantistici. Questa esercitazione si concentra sul test di Hadamard. Dato che $\mathcal{K}_U$ ha dimensione $r$, l'hamiltoniana proiettata nel sottospazio avrà dimensione $r \times r$. Con $r$ sufficientemente piccola (in genere $r<<100$ è sufficiente per ottenere la convergenza delle stime delle autovalori) possiamo quindi facilmente diagonalizzare l'hamiltoniana proiettata $\tilde{H}$. Tuttavia, non possiamo diagonalizzare direttamente $\tilde{H}$ a causa della non ortogonalità dei vettori dello spazio di Krylov. Dovremo misurare le loro sovrapposizioni e costruire una matrice $\tilde{S}$

Questo ci permette di risolvere il problema degli autovalori in uno spazio non ortogonale (detto anche problema degli autovalori generalizzato)

Si possono quindi ottenere stime degli autovalori e degli autostati di $H$ osservando quelli di $\tilde{H}$. Ad esempio, la stima dell'energia dello stato fondamentale si ottiene prendendo l'autovalore più piccolo $c$ e lo stato fondamentale dal corrispondente autovettore $\vec{c}$. I coefficienti in $\vec{c}$ determinano il contributo dei diversi vettori che coprono $\mathcal{K}_U$.

La figura mostra una rappresentazione circuitale del test di Hadamard modificato, un metodo utilizzato per calcolare la sovrapposizione tra diversi stati quantistici. Per ogni elemento della matrice $\tilde{H}_{i,j}$, si esegue un test di Hadamard tra gli stati $\vert \psi_i \rangle$, $\vert \psi_j \rangle$. Ciò è evidenziato nella figura dallo schema di colori degli elementi della matrice e dalle corrispondenti operazioni $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$. Pertanto, per calcolare tutti gli elementi della matrice dell'hamiltoniana proiettata $\tilde{H}$ è necessario un insieme di test di Hadamard per tutte le possibili combinazioni di vettori dello spazio di Krylov. Il filo superiore del circuito del test di Hadamard è un qubit ancilla che viene misurato nella base X o Y; il suo valore di aspettativa determina il valore della sovrapposizione tra gli stati. Il filo inferiore rappresenta tutti i qubit dell'hamiltoniana del sistema. L'operazione $\text{Prep} \; \psi_i$ prepara il qubit del sistema nello stato $\vert \psi_i \rangle$ controllato dallo stato del qubit ancilla (analogamente per $\text{Prep} \; \psi_j$ ) e l'operazione $P$ rappresenta la decomposizione Pauli dell'hamiltoniana del sistema $H = \sum_i P_i$. Di seguito viene fornita una derivazione più dettagliata delle operazioni calcolate dal test di Hadamard.

Definire l'Hamiltoniano

Consideriamo l'hamiltoniana di Heisenberg per $N$ qubit su una catena lineare: $H= \sum_{i,j}^N X_i X_j + Y_i Y_j - J Z_i Z_j$

# Define problem Hamiltonian.
n_qubits = 30
J = 1  # coupling strength for ZZ interaction

# Define the Hamiltonian:
H_int = [["I"] * n_qubits for _ in range(3 * (n_qubits - 1))]
for i in range(n_qubits - 1):
    H_int[i][i] = "Z"
    H_int[i][i + 1] = "Z"
for i in range(n_qubits - 1):
    H_int[n_qubits - 1 + i][i] = "X"
    H_int[n_qubits - 1 + i][i + 1] = "X"
for i in range(n_qubits - 1):
    H_int[2 * (n_qubits - 1) + i][i] = "Y"
    H_int[2 * (n_qubits - 1) + i][i + 1] = "Y"
H_int = ["".join(term) for term in H_int]
H_tot = [(term, J) if term.count("Z") == 2 else (term, 1) for term in H_int]

# Get operator
H_op = SparsePauliOp.from_list(H_tot)
print(H_tot)

Output:

[('ZZIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IZZIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIZZIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIZZIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIZZIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIZZIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIZZIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIZZIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIZZIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIZZIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIZZIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIZZIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIZZIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIZZIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIZZIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIZZIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIZZIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIZZIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIZZIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIZZIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIZZIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIZZIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIZZIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIZZIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIZZIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIZZIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIZZII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIZZI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIZZ', 1), ('XXIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IXXIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIXXIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIXXIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIXXIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIXXIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIXXIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIXXIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIXXIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIXXIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIXXIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIXXIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIXXIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIXXIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIXXIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIXXIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIXXIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIXXIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIXXIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIXXIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIXXIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIXXIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIXXIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIXXIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIXXIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIXXIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIXXII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIXXI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIXX', 1), ('YYIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IYYIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIYYIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIYYIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIYYIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIYYIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIYYIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIYYIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIYYIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIYYIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIYYIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIYYIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIYYIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIYYIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIYYIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIYYIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIYYIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIYYIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIYYIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIYYIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIYYIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIYYIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIYYIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIYYIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIYYIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIYYIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIYYII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIYYI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIYY', 1)]

Imposta i parametri per l'algoritmo

Scegliamo euristicamente un valore per il passo temporale dt (basato sui limiti superiori della norma hamiltoniana). Ref [2] ha dimostrato che un passo temporale sufficientemente piccolo è $\pi/\vert \vert H \vert \vert$, e che è preferibile fino a un certo punto sottostimare questo valore piuttosto che sovrastimarlo, poiché una sovrastima può permettere ai contributi degli stati ad alta energia di corrompere anche lo stato ottimale nello spazio di Krylov. D'altra parte, la scelta di $dt$ troppo piccola porta a un condizionamento peggiore del sottospazio di Krylov, poiché i vettori base di Krylov differiscono meno da un passo all'altro.

# Get Hamiltonian restricted to single-particle states
single_particle_H = np.zeros((n_qubits, n_qubits))
for i in range(n_qubits):
    for j in range(i + 1):
        for p, coeff in H_op.to_list():
            p_x = Pauli(p).x
            p_z = Pauli(p).z
            if all(
                p_x[k] == ((i == k) + (j == k)) % 2 for k in range(n_qubits)
            ):
                sgn = (
                    (-1j) ** sum(p_z[k] and p_x[k] for k in range(n_qubits))
                ) * ((-1) ** p_z[i])
            else:
                sgn = 0
            single_particle_H[i, j] += sgn * coeff
for i in range(n_qubits):
    for j in range(i + 1, n_qubits):
        single_particle_H[i, j] = np.conj(single_particle_H[j, i])

# Set dt according to spectral norm
dt = np.pi / np.linalg.norm(single_particle_H, ord=2)
dt

Output:

np.float64(0.10833078115826875)

E impostare altri parametri dell'algoritmo. Ai fini di questa esercitazione, ci limiteremo a utilizzare uno spazio di Krylov con solo cinque dimensioni, che è abbastanza limitante.

# Set parameters for quantum Krylov algorithm
krylov_dim = 5  # size of Krylov subspace
num_trotter_steps = 6
dt_circ = dt / num_trotter_steps

Preparazione dello Stato

Scegliere uno stato di riferimento $\vert \psi \rangle$ che abbia una certa sovrapposizione con lo stato fondamentale. Per questa hamiltoniana, utilizziamo come stato di riferimento lo stato a con un'eccitazione nel qubit centrale $\vert 00..010...00 \rangle$.

qc_state_prep = QuantumCircuit(n_qubits)
qc_state_prep.x(int(n_qubits / 2) + 1)
qc_state_prep.draw("mpl", scale=0.5)

Output:

Evoluzione temporale

Possiamo realizzare l'operatore di evoluzione temporale generato da una data hamiltoniana: $U=e^{-iHt}$ tramite l' approssimazione di Lie-Trotter.

t = Parameter("t")

## Create the time-evo op circuit
evol_gate = PauliEvolutionGate(
    H_op, time=t, synthesis=LieTrotter(reps=num_trotter_steps)
)

qr = QuantumRegister(n_qubits)
qc_evol = QuantumCircuit(qr)
qc_evol.append(evol_gate, qargs=qr)

Output:

<qiskit.circuit.instructionset.InstructionSet at 0x11eef9be0>

test di Hadamard

Dove $P$ è uno dei termini della decomposizione dell'hamiltoniana $H=\sum P$ e $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$ sono operazioni controllate che preparano $|\psi_i\rangle$, $|\psi_j\rangle$ vettori dello spazio di Krylov unitario, con $|\psi_k\rangle = e^{-i H k dt } \vert \psi \rangle = e^{-i H k dt } U_{\psi} \vert 0 \rangle^N$. Per misurare $X$, applicare prima $H$...

... poi misurare:

Dall'identità $|a + b\|^2 = \langle a + b | a + b \rangle = \|a\|^2 + \|b\|^2 + 2\text{Re}\langle a | b \rangle$. Analogamente, misurando $Y$ si ottiene

## Create the time-evo op circuit
evol_gate = PauliEvolutionGate(
    H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)
)

## Create the time-evo op dagger circuit
evol_gate_d = PauliEvolutionGate(
    H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)
)
evol_gate_d = evol_gate_d.inverse()

# Put pieces together
qc_reg = QuantumRegister(n_qubits)
qc_temp = QuantumCircuit(qc_reg)
qc_temp.compose(qc_state_prep, inplace=True)
for _ in range(num_trotter_steps):
    qc_temp.append(evol_gate, qargs=qc_reg)
for _ in range(num_trotter_steps):
    qc_temp.append(evol_gate_d, qargs=qc_reg)
qc_temp.compose(qc_state_prep.inverse(), inplace=True)

# Create controlled version of the circuit
controlled_U = qc_temp.to_gate().control(1)

# Create hadamard test circuit for real part
qr = QuantumRegister(n_qubits + 1)
qc_real = QuantumCircuit(qr)
qc_real.h(0)
qc_real.append(controlled_U, list(range(n_qubits + 1)))
qc_real.h(0)

print(
    "Circuit for calculating the real part of the overlap in S via Hadamard test"
)
qc_real.draw("mpl", fold=-1, scale=0.5)

Output:

Circuit for calculating the real part of the overlap in S via Hadamard test

Il circuito di test Hadamard può essere un circuito profondo una volta decomposto in porte native (il che aumenterà ancora di più se teniamo conto della topologia del dispositivo)

print(
    "Number of layers of 2Q operations",
    qc_real.decompose(reps=2).depth(lambda x: x[0].num_qubits == 2),
)

Output:

Number of layers of 2Q operations 112753

Fase 2: Ottimizzazione del problema per l'esecuzione su hardware quantistico

Test di Hadamard efficiente

Possiamo ottimizzare i circuiti profondi per il test di Hadamard che abbiamo ottenuto introducendo alcune approssimazioni e facendo affidamento su alcune ipotesi sull'hamiltoniana del modello. Ad esempio, si consideri il seguente circuito per il test di Hadamard:

Supponiamo di poter calcolare classicamente $E_0$, l'autovalore di $|0\rangle^N$ sotto l'hamiltoniana $H$. Questo è soddisfatto quando l'hamiltoniana preserva la simmetria U(1). Anche se questa può sembrare un'ipotesi forte, in molti casi è lecito supporre che esista uno stato di vuoto (che in questo caso corrisponde allo stato $|0\rangle^N$ ) che non è influenzato dall'azione dell'hamiltoniana. Questo è vero, ad esempio, per gli hamiltoniani della chimica che descrivono molecole stabili (dove il numero di elettroni è conservato). Dato che il gate $\text{Prep} \; \psi$, prepara lo stato di riferimento desiderato $\ket{psi} = \text{Prep} \; \psi \ket{0} = e^{-i H 0 dt} U_{\psi} \ket{0}$, ad esempio, preparare lo stato HF per la chimica $\text{Prep} \; \psi$ sarebbe un prodotto di NOT a singolo-qubit, quindi controllato- $\text{Prep} \; \psi$ è solo un prodotto di CNOT. Il circuito di cui sopra implementa il seguente stato prima della misurazione:

dove abbiamo utilizzato il classico spostamento di fase simulabile $U\ket{0}^N = e^{i\phi}\ket{0}^N$ nella terza riga. Pertanto i valori di aspettativa si ottengono come

Utilizzando queste ipotesi, siamo stati in grado di scrivere i valori di aspettativa degli operatori di interesse con un minor numero di operazioni controllate. In effetti, dobbiamo implementare solo la preparazione controllata degli stati $\text{Prep} \; \psi$ e non le evoluzioni temporali controllate. Riformulare il calcolo come sopra ci permetterà di ridurre notevolmente la profondità dei circuiti risultanti.

Scomporre l'operatore di evoluzione temporale con la decomposizione di Trotter

Invece di implementare esattamente l'operatore di evoluzione temporale, possiamo utilizzare la decomposizione di Trotter per implementarne un'approssimazione. Ripetendo più volte una decomposizione di Trotter di un certo ordine si ottiene un'ulteriore riduzione dell'errore introdotto dall'approssimazione. Di seguito, costruiamo direttamente l'implementazione di Trotter nel modo più efficiente per il grafo delle interazioni dell'hamiltoniana che stiamo considerando (solo interazioni di prossimità). In pratica inseriamo le rotazioni Pauli $R_{xx}$, $R_{yy}$, $R_{zz}$ con un angolo parametrizzato $t$ che corrispondono all'implementazione approssimata di $e^{-i (XX + YY + ZZ) t}$. Data la differenza di definizione delle rotazioni Pauli e dell'evoluzione temporale che stiamo cercando di implementare, dovremo utilizzare il parametro $2*dt$ per ottenere un'evoluzione temporale di $dt$. Inoltre, invertiamo l'ordine delle operazioni per un numero dispari di ripetizioni dei passi di Trotter, il che è funzionalmente equivalente ma permette di sintetizzare operazioni adiacenti in un unico $SU(2)$ unitario. In questo modo si ottiene un circuito molto meno profondo di quello che si ottiene utilizzando la funzionalità generica PauliEvolutionGate() .

t = Parameter("t")

# Create instruction for rotation about XX+YY-ZZ:
Rxyz_circ = QuantumCircuit(2)
Rxyz_circ.rxx(t, 0, 1)
Rxyz_circ.ryy(t, 0, 1)
Rxyz_circ.rzz(t, 0, 1)
Rxyz_instr = Rxyz_circ.to_instruction(label="RXX+YY+ZZ")

interaction_list = [
    [[i, i + 1] for i in range(0, n_qubits - 1, 2)],
    [[i, i + 1] for i in range(1, n_qubits - 1, 2)],
]  # linear chain

qr = QuantumRegister(n_qubits)
trotter_step_circ = QuantumCircuit(qr)
for i, color in enumerate(interaction_list):
    for interaction in color:
        trotter_step_circ.append(Rxyz_instr, interaction)
    if i < len(interaction_list) - 1:
        trotter_step_circ.barrier()
reverse_trotter_step_circ = trotter_step_circ.reverse_ops()

qc_evol = QuantumCircuit(qr)
for step in range(num_trotter_steps):
    if step % 2 == 0:
        qc_evol = qc_evol.compose(trotter_step_circ)
    else:
        qc_evol = qc_evol.compose(reverse_trotter_step_circ)

qc_evol.decompose().draw("mpl", fold=-1, scale=0.5)

Output:

Utilizzare un circuito ottimizzato per la preparazione dello stato

control = 0
excitation = int(n_qubits / 2) + 1
controlled_state_prep = QuantumCircuit(n_qubits + 1)
controlled_state_prep.cx(control, excitation)
controlled_state_prep.draw("mpl", fold=-1, scale=0.5)

Output:

Circuiti modello per il calcolo degli elementi di matrice di $\tilde{S}$ e $\tilde{H}$ tramite il test di Hadamard

L'unica differenza tra i circuiti utilizzati nel test di Hadamard sarà la fase dell'operatore di evoluzione temporale e le osservabili misurate. Possiamo quindi preparare un circuito modello che rappresenti il circuito generico per il test di Hadamard, con segnaposti per le porte che dipendono dall'operatore di evoluzione temporale.

# Parameters for the template circuits
parameters = []
for idx in range(1, krylov_dim):
    parameters.append(2 * dt_circ * (idx))

# Create modified hadamard test circuit
qr = QuantumRegister(n_qubits + 1)
qc = QuantumCircuit(qr)
qc.h(0)
qc.compose(controlled_state_prep, list(range(n_qubits + 1)), inplace=True)
qc.barrier()
qc.compose(qc_evol, list(range(1, n_qubits + 1)), inplace=True)
qc.barrier()
qc.x(0)
qc.compose(
    controlled_state_prep.inverse(), list(range(n_qubits + 1)), inplace=True
)
qc.x(0)

qc.decompose().draw("mpl", fold=-1)

Output:

print(
    "The optimized circuit has 2Q gates depth: ",
    qc.decompose().decompose().depth(lambda x: x[0].num_qubits == 2),
)

Output:

The optimized circuit has 2Q gates depth:  74

Abbiamo ridotto considerevolmente la profondità del test di Hadamard con una combinazione di approssimazione di Trotter e unitarietà non controllata

Passaggio 3: eseguire utilizzando Qiskit primitives

Istanziare il backend e impostare i parametri di runtime

service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
if (
    "if_else" not in backend.target.operation_names
):  # Needed as "op_name" could be "if_else"
    backend.target.add_instruction(IfElseOp, name="if_else")
print(backend.name)

Transpiling su una QPU

Per prima cosa, selezioniamo sottoinsiemi della mappa di accoppiamento con qubit dalle prestazioni "buone" (dove "buone" è piuttosto arbitrario in questo caso, vogliamo principalmente evitare qubit dalle prestazioni davvero scadenti) e creiamo un nuovo obiettivo per la transpilazione

target = backend.target
cmap = target.build_coupling_map(filter_idle_qubits=True)
cmap_list = list(cmap.get_edges())

cust_cmap_list = copy.deepcopy(cmap_list)
for q in range(target.num_qubits):
    meas_err = target["measure"][(q,)].error
    t2 = target.qubit_properties[q].t2 * 1e6
    if meas_err > 0.02 or t2 < 100:
        for q_pair in cmap_list:
            if q in q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except:
                    continue

for q in cmap_list:
    op_name = list(target.operation_names_for_qargs(q))[0]
    twoq_gate_err = target[f"{op_name}"][q].error
    if twoq_gate_err > 0.005:
        for q_pair in cmap_list:
            if q == q_pair:
                try:
                    cust_cmap_list.remove(q)
                except:
                    continue


cust_cmap = CouplingMap(cust_cmap_list)
cust_target = Target.from_configuration(
    basis_gates=backend.configuration().basis_gates,
    coupling_map=cust_cmap,
)

Quindi trasporre il circuito virtuale nel miglior layout fisico in questo nuovo target

basis_gates = list(target.operation_names)
pm = generate_preset_pass_manager(
    optimization_level=3,
    target=cust_target,
    basis_gates=basis_gates,
)

qc_trans = pm.run(qc)

print("depth", qc_trans.depth(lambda x: x[0].num_qubits == 2))
print("num 2q ops", qc_trans.count_ops())
print(
    "physical qubits",
    sorted(
        [
            idx
            for idx, qb in qc_trans.layout.initial_layout.get_physical_bits().items()
            if qb._register.name != "ancilla"
        ]
    ),
)

Output:

depth 52
num 2q ops OrderedDict([('rz', 2058), ('sx', 1703), ('cz', 728), ('x', 84), ('barrier', 8)])
physical qubits [91, 92, 93, 94, 95, 98, 99, 108, 109, 110, 111, 113, 114, 115, 119, 127, 132, 133, 134, 135, 137, 139, 147, 148, 149, 150, 151, 152, 153, 154, 155]

Creare PUB per l'esecuzione con Estimator

# Define observables to measure for S
observable_S_real = "I" * (n_qubits) + "X"
observable_S_imag = "I" * (n_qubits) + "Y"

observable_op_real = SparsePauliOp(
    observable_S_real
)  # define a sparse pauli operator for the observable
observable_op_imag = SparsePauliOp(observable_S_imag)

layout = qc_trans.layout  # get layout of transpiled circuit
observable_op_real = observable_op_real.apply_layout(
    layout
)  # apply physical layout to the observable
observable_op_imag = observable_op_imag.apply_layout(layout)
observable_S_real = (
    observable_op_real.paulis.to_labels()
)  # get the label of the physical observable
observable_S_imag = observable_op_imag.paulis.to_labels()

observables_S = [[observable_S_real], [observable_S_imag]]


# Define observables to measure for H
# Hamiltonian terms to measure
observable_list = []
for pauli, coeff in zip(H_op.paulis, H_op.coeffs):
    # print(pauli)
    observable_H_real = pauli[::-1].to_label() + "X"
    observable_H_imag = pauli[::-1].to_label() + "Y"
    observable_list.append([observable_H_real])
    observable_list.append([observable_H_imag])

layout = qc_trans.layout

observable_trans_list = []
for observable in observable_list:
    observable_op = SparsePauliOp(observable)
    observable_op = observable_op.apply_layout(layout)
    observable_trans_list.append([observable_op.paulis.to_labels()])

observables_H = observable_trans_list


# Define a sweep over parameter values
params = np.vstack(parameters).T


# Estimate the expectation value for all combinations of
# observables and parameter values, where the pub result will have
# shape (# observables, # parameter values).
pub = (qc_trans, observables_S + observables_H, params)

Esegui circuiti

I circuiti per $t=0$ sono calcolabili in modo classico

qc_cliff = qc.assign_parameters({t: 0})


# Get expectation values from experiment
S_expval_real = StabilizerState(qc_cliff).expectation_value(
    Pauli("I" * (n_qubits) + "X")
)
S_expval_imag = StabilizerState(qc_cliff).expectation_value(
    Pauli("I" * (n_qubits) + "Y")
)

# Get expectation values
S_expval = S_expval_real + 1j * S_expval_imag

H_expval = 0
for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):
    # Get expectation values from experiment
    expval_real = StabilizerState(qc_cliff).expectation_value(
        Pauli(pauli[::-1].to_label() + "X")
    )
    expval_imag = StabilizerState(qc_cliff).expectation_value(
        Pauli(pauli[::-1].to_label() + "Y")
    )
    expval = expval_real + 1j * expval_imag

    # Fill-in matrix elements
    H_expval += coeff * expval


print(H_expval)

Output:

(25+0j)

Esegui i circuiti per $S$ e $\tilde{H}$ con Estimator

# Experiment options
num_randomizations = 300
num_randomizations_learning = 30
shots_per_randomization = 100
noise_factors = [1, 1.2, 1.4]
learning_pair_depths = [0, 4, 24, 48]


experimental_opts = {}
experimental_opts["resilience"] = {
    "measure_mitigation": True,
    "measure_noise_learning": {
        "num_randomizations": num_randomizations_learning,
        "shots_per_randomization": shots_per_randomization,
    },
    "zne_mitigation": True,
    "zne": {"noise_factors": noise_factors},
    "layer_noise_learning": {
        "max_layers_to_learn": 10,
        "layer_pair_depths": learning_pair_depths,
        "shots_per_randomization": shots_per_randomization,
        "num_randomizations": num_randomizations_learning,
    },
    "zne": {
        "amplifier": "pea",
        "extrapolated_noise_factors": [0] + noise_factors,
    },
}
experimental_opts["twirling"] = {
    "num_randomizations": num_randomizations,
    "shots_per_randomization": shots_per_randomization,
    "strategy": "all",
}

estimator = Estimator(mode=backend, options=experimental_opts)


job = estimator.run([pub])

Fase 4: Post-elaborazione e restituzione del risultato nel formato classico desiderato

results = job.result()[0]

Calcolare le matrici Hamiltoniane efficaci e di sovrapposizione

Calcolare innanzitutto la fase accumulata dallo stato $\vert 0 \rangle$ durante l'evoluzione temporale non controllata

prefactors = [
    np.exp(-1j * sum([c for p, c in H_op.to_list() if "Z" in p]) * i * dt)
    for i in range(1, krylov_dim)
]

Una volta ottenuti i risultati delle esecuzioni dei circuiti, possiamo post-elaborare i dati per calcolare gli elementi della matrice di $S$

# Assemble S, the overlap matrix of dimension D:
S_first_row = np.zeros(krylov_dim, dtype=complex)
S_first_row[0] = 1 + 0j

# Add in ancilla-only measurements:
for i in range(krylov_dim - 1):
    # Get expectation values from experiment
    expval_real = results.data.evs[0][0][
        i
    ]  # automatic extrapolated evs if ZNE is used
    expval_imag = results.data.evs[1][0][
        i
    ]  # automatic extrapolated evs if ZNE is used

    # Get expectation values
    expval = expval_real + 1j * expval_imag
    S_first_row[i + 1] += prefactors[i] * expval

S_first_row_list = S_first_row.tolist()  # for saving purposes


S_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)

# Distribute entries from first row across matrix:
for i, j in it.product(range(krylov_dim), repeat=2):
    if i >= j:
        S_circ[j, i] = S_first_row[i - j]
    else:
        S_circ[j, i] = np.conj(S_first_row[j - i])

Matrix(S_circ)

Output:

E gli elementi della matrice di $\tilde{H}$

# Assemble S, the overlap matrix of dimension D:
H_first_row = np.zeros(krylov_dim, dtype=complex)
H_first_row[0] = H_expval

for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):
    # Add in ancilla-only measurements:
    for i in range(krylov_dim - 1):
        # Get expectation values from experiment
        expval_real = results.data.evs[2 + 2 * obs_idx][0][
            i
        ]  # automatic extrapolated evs if ZNE is used
        expval_imag = results.data.evs[2 + 2 * obs_idx + 1][0][
            i
        ]  # automatic extrapolated evs if ZNE is used

        # Get expectation values
        expval = expval_real + 1j * expval_imag
        H_first_row[i + 1] += prefactors[i] * coeff * expval

H_first_row_list = H_first_row.tolist()

H_eff_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)

# Distribute entries from first row across matrix:
for i, j in it.product(range(krylov_dim), repeat=2):
    if i >= j:
        H_eff_circ[j, i] = H_first_row[i - j]
    else:
        H_eff_circ[j, i] = np.conj(H_first_row[j - i])

Matrix(H_eff_circ)

Output:

Infine, possiamo risolvere il problema generalizzato degli autovalori per $\tilde{H}$ :

$\tilde{H} \vec{c} = c S \vec{c}$

e ottenere una stima dell'energia dello stato fondamentale $c_{min}$

gnd_en_circ_est_list = []
for d in range(1, krylov_dim + 1):
    # Solve generalized eigenvalue problem for different size of the Krylov space
    gnd_en_circ_est = solve_regularized_gen_eig(
        H_eff_circ[:d, :d], S_circ[:d, :d], threshold=9e-1
    )
    gnd_en_circ_est_list.append(gnd_en_circ_est)
    print("The estimated ground state energy is: ", gnd_en_circ_est)

Output:

The estimated ground state energy is:  25.0
The estimated ground state energy is:  22.572154819954875
The estimated ground state energy is:  21.691509219286587
The estimated ground state energy is:  21.23882298756386
The estimated ground state energy is:  20.965499325470294

Per un settore a una sola particella, possiamo calcolare in modo efficiente lo stato fondamentale di questo settore dell'hamiltoniana in modo classico

gs_en = single_particle_gs(H_op, n_qubits)

Output:

n_sys_qubits 30
n_exc 1 , subspace dimension 31
single particle ground state energy:  21.021912418526906

plt.plot(
    range(1, krylov_dim + 1),
    gnd_en_circ_est_list,
    color="blue",
    linestyle="-.",
    label="KQD estimate",
)
plt.plot(
    range(1, krylov_dim + 1),
    [gs_en] * krylov_dim,
    color="red",
    linestyle="-",
    label="exact",
)
plt.xticks(range(1, krylov_dim + 1), range(1, krylov_dim + 1))
plt.legend()
plt.xlabel("Krylov space dimension")
plt.ylabel("Energy")
plt.title(
    "Estimating Ground state energy with Krylov Quantum Diagonalization"
)
plt.show()

Output:

Appendice: sottospazio di Krylov dalle evoluzioni in tempo reale

Lo spazio di Krylov unitario è definito come

per qualche passo temporale $dt$ che determineremo in seguito. Assumiamo temporaneamente che $r$ sia pari: allora definiamo $d=r/2$. Si noti che quando proiettiamo l'Hamiltoniana nello spazio di Krylov di cui sopra, essa è indistinguibile dallo spazio di Krylov

cioè dove tutte le evoluzioni temporali sono spostate all'indietro di $d$ passi temporali. Il motivo per cui è indistinguibile è che gli elementi della matrice

sono invarianti sotto spostamenti complessivi del tempo di evoluzione, poiché le evoluzioni temporali commutano con l'hamiltoniana. Per $r$ dispari, possiamo utilizzare l'analisi per $r-1$.

Vogliamo dimostrare che da qualche parte in questo spazio di Krylov è garantita l'esistenza di uno stato a bassa energia. Lo facciamo attraverso il seguente risultato, derivato dal teorema 3.1 di [3] :

Asserzione 1: esiste una funzione $f$ tale che per le energie $E$ nell'intervallo spettrale dell'hamiltoniana (cioè tra l'energia dello stato fondamentale e l'energia massima)...

1. $f(E_0)=1$
2. $|f(E)|\le2\left(1 + \delta\right)^{-d}$ per tutti i valori di $E$ che si trovano a $\ge\delta$ distanza da $E_0$, cioè è soppresso esponenzialmente
3. $f(E)$ è una combinazione lineare di $e^{ijE\,dt}$ per $j=-d,-d+1,...,d-1,d$

Di seguito viene fornita una dimostrazione, che può essere tranquillamente saltata a meno che non si voglia comprendere l'argomentazione completa e rigorosa. Per ora ci concentriamo sulle implicazioni di questa affermazione. In base alla proprietà 3, possiamo vedere che lo spazio di Krylov spostato sopra contiene lo stato $f(H)|\psi\rangle$. Questo è il nostro stato a bassa energia. Per capire perché, scriviamo $|\psi\rangle$ nell'autobase dell'energia:

dove $|E_k\rangle$ è l'autostato di energia k e $\gamma_k$ è la sua ampiezza nello stato iniziale $|\psi\rangle$. Espresso in termini di questo, $f(H)|\psi\rangle$ è dato da

utilizzando il fatto che possiamo sostituire $H$ con $E_k$ quando agisce sull'autostato $|E_k\rangle$. L'errore energetico di questo stato è quindi

Per trasformare questo risultato in un limite superiore di più facile comprensione, separiamo prima la somma nel numeratore in termini con $E_k-E_0\le\delta$ e termini con $E_k-E_0>\delta$ :

Il primo termine può essere limitato in alto da $\delta$,

dove il primo passo segue perché $E_k-E_0\le\delta$ per ogni $E_k$ nella somma, e il secondo passo segue perché la somma nel numeratore è un sottoinsieme della somma nel denominatore. Per il secondo termine, prima di tutto abbassiamo il denominatore di $|\gamma_0|^2$, dato che $f(E_0)^2=1$ : sommando il tutto, si ottiene

Per semplificare ciò che rimane, si noti che per tutti questi $E_k$, dalla definizione di $f$ sappiamo che $f(E_k)^2 \le 4\left(1 + \delta\right)^{-2d}$. Inoltre, l'upper bounding di $E_k-E_0<2\|H\|$ e l'upper bounding di $\sum_{E_k>E_0+\delta}|\gamma_k|^2<1$ dà

Questo vale per qualsiasi $\delta>0$, quindi se impostiamo $\delta$ uguale all'errore che ci prefiggiamo, il limite dell'errore di cui sopra converge verso tale valore in modo esponenziale con la dimensione di Krylov $2d=r$. Si noti inoltre che se $\delta<E_1-E_0$ allora il termine $\delta$ scompare del tutto nel limite di cui sopra.

Per completare l'argomentazione, notiamo innanzitutto che quanto sopra è solo l'errore energetico del particolare stato $f(H)|\psi\rangle$, piuttosto che l'errore energetico dello stato a più bassa energia nello spazio di Krylov. Tuttavia, in base al principio variazionale (di Rayleigh-Ritz), l'errore di energia dello stato a più bassa energia nello spazio di Krylov è limitato dall'errore di energia di qualsiasi stato nello spazio di Krylov, quindi quanto sopra è anche un limite superiore sull'errore di energia dello stato a più bassa energia, cioè l'uscita dell'algoritmo di diagonalizzazione quantistica di Krylov.

È possibile eseguire un'analisi simile a quella precedente, che tenga conto anche del rumore e della procedura di sogliatura descritta nel quaderno. Per questa analisi si veda [2] e [4].

Appendice: prova della rivendicazione 1

Quanto segue deriva in gran parte da [3], Teorema 3.1: Sia $0 < a < b$ e sia $\Pi^*_d$ lo spazio dei polinomi residui (polinomi il cui valore a 0 è 1) di grado massimo $d$. La soluzione di

è

e il corrispondente valore minimo è

Vogliamo convertire questa funzione in una funzione che possa essere espressa naturalmente in termini di esponenziali complessi, perché queste sono le reali evoluzioni temporali che generano lo spazio di Krylov quantistico. Per fare ciò, è conveniente introdurre la seguente trasformazione delle energie all'interno dell'intervallo spettrale dell'hamiltoniano in numeri nell'intervallo $[0,1]$ : definire

dove $dt$ è un passo temporale tale che $-\pi < E_0dt < E_\text{max}dt < \pi$. Si noti che $g(E_0)=0$ e $g(E)$ crescono quando $E$ si allontana da $E_0$.

Ora, utilizzando il polinomio $p^*(x)$ con i parametri a, b, d impostati a $a = g(E_0 + \delta)$, $b = 1$, e d = int( r/2 ), definiamo la funzione:

dove $E_0$ è l'energia dello stato fondamentale. Inserendo $\cos(x)=\frac{e^{ix}+e^{-ix}}{2}$ si vede che $f(E)$ è un polinomio trigonometrico di grado $d$, cioè una combinazione lineare di $e^{ijE\,dt}$ per $j=-d,-d+1,...,d-1,d$. Inoltre, dalla definizione di $p^*(x)$ abbiamo che $f(E_0)=p(0)=1$ e per qualsiasi $E$ nell'intervallo spettrale tale che $\vert E-E_0 \vert > \delta$ abbiamo