Diagonalizzazione quantistica di Krylov

In questa lezione sulla diagonalizzazione quantistica di Krylov (KQD) risponderemo a quanto segue:

• Che cos'è il metodo di Krylov, in generale?
• Perché il metodo di Krylov funziona e a quali condizioni?
• Che ruolo ha il calcolo quantistico?

La parte quantistica dei calcoli si basa in gran parte sul lavoro svolto in Ref [1].

Il video che segue offre una panoramica dei metodi di Krylov nell'informatica classica, ne motiva l'uso e spiega come l'informatica quantistica possa svolgere un ruolo in questo flusso di lavoro. Il testo successivo offre maggiori dettagli e implementa un metodo di Krylov sia in modo classico, sia utilizzando un computer quantistico.

1. Introduzione ai metodi di Krylov

Un metodo del sottospazio di Krylov può riferirsi a uno qualsiasi dei diversi metodi costruiti attorno al cosiddetto sottospazio di Krylov. Una rassegna completa di questi argomenti va al di là dello scopo di questa lezione, ma i Rif. [2-4] possono fornire un quadro sostanzialmente più completo. In questa sede ci concentreremo su cosa sia un sottospazio di Krylov, su come e perché sia utile per risolvere problemi di autovalori e infine su come possa essere implementato su un computer quantistico.

Definizione: Data una matrice $N\times N$ simmetrica, positiva e semidefinita $A$, 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\leq N$, con un vettore di riferimento $\vert v \rangle$.

Sebbene i vettori di cui sopra coprano quello che chiamiamo un sottospazio di Krylov, non c'è motivo di pensare che siano ortogonali. Spesso si utilizza un processo iterativo di ortonormalizzazione simile all' ortogonalizzazione di Gram-Schmidt. In questo caso il processo è leggermente diverso, poiché ogni nuovo vettore viene reso ortogonale agli altri man mano che viene generato. In questo contesto si parla di iterazione di Arnoldi. Partendo dal vettore iniziale $|v\rangle$, si genera il vettore successivo $A|v\rangle$, e si assicura che questo secondo vettore sia ortogonale al primo sottraendo la sua proiezione su $|v\rangle$. Ovvero

È ormai facile capire che $|v_0\rangle \perp |v_1\rangle$, poiché

Facciamo lo stesso per il vettore successivo, assicurandoci che sia ortogonale ai due precedenti:

Se si ripete questo procedimento per tutti i vettori $r$, si ottiene una base ortonormale completa per uno spazio di Krylov. Si noti che il processo di ortogonalizzazione in questo caso darà come risultato zero una volta $r>m$, poiché $m$ vettore ortogonale copre necessariamente l'intero spazio. Il processo darà anche zero se un vettore è un autovettore di $A$, poiché tutti i vettori successivi saranno multipli di quel vettore.

1.1 Un esempio semplice: Krylov a mano

Esaminiamo la generazione di un sottospazio di Krylov su una matrice banalmente piccola, in modo da vedere il processo. Partiamo da una matrice iniziale $A$ di nostro interesse:

Per questo piccolo esempio, possiamo determinare facilmente gli autovettori e gli autovalori anche a mano. Mostriamo qui la soluzione numerica.

# One might use linalg.eigh here, but later matrices may not be Hermitian. So we use
# linalg.eig in this lesson.

import numpy as np

A = np.array([[4, -1, 0], [-1, 4, -1], [0, -1, 4]])
eigenvalues, eigenvectors = np.linalg.eig(A)
print("The eigenvalues are ", eigenvalues)
print("The eigenvectors are ", eigenvectors)

Output:

The eigenvalues are  [2.58578644 4.         5.41421356]
The eigenvectors are  [[ 5.00000000e-01 -7.07106781e-01  5.00000000e-01]
 [ 7.07106781e-01  1.37464400e-16 -7.07106781e-01]
 [ 5.00000000e-01  7.07106781e-01  5.00000000e-01]]

Li registriamo qui per un successivo confronto:

Vorremmo studiare come questo processo funziona (o fallisce) quando aumentiamo la dimensione del nostro sottospazio di Krylov, $r$. A tal fine, applicheremo questo processo:

• Generare un sottospazio dello spazio vettoriale completo partendo da un vettore scelto a caso $|v\rangle$ (chiamarlo $|v_0\rangle$ se è già normalizzato, come sopra).
• Proiettare la matrice completa $A$ su quel sottospazio e trovare gli autovalori della matrice proiettata $\tilde{A}$.
• Aumentare la dimensione del sottospazio generando più vettori, assicurandosi che siano ortonormali, utilizzando un processo simile all'ortogonalizzazione di Gram-Schmidt.
• Proiettare $A$ sul sottospazio più grande e trovare gli autovalori della matrice risultante, $\tilde{A}$.
• Ripetere l'operazione finché gli autovalori non convergono (o, in questo caso giocattolo, finché non si sono generati vettori che coprono l'intero spazio vettoriale della matrice originale $A$ ).

Una normale implementazione del metodo di Krylov non dovrebbe risolvere il problema degli autovalori per la matrice proiettata su ogni sottospazio di Krylov man mano che viene costruito. Si può costruire il sottospazio della dimensione desiderata, proiettare la matrice su tale sottospazio e diagonalizzare la matrice proiettata. La proiezione e la diagonalizzazione in ogni dimensione del sottospazio vengono effettuate solo per verificare la convergenza.

$r=1$ i dimensionali:

Scegliamo un vettore casuale, ad esempio

Se non è già normalizzato, normalizzarlo.

Ora proiettiamo la nostra matrice $A$ sul sottospazio di questo unico vettore:

Questa è la proiezione della matrice sul nostro sottospazio di Krylov quando contiene un solo vettore, $|v_0\rangle$. L'autovalore di questa matrice è banalmente 4. Si tratta di una stima di ordine zerotico degli autovalori (in questo caso uno solo) di $A$. Sebbene si tratti di una stima insufficiente, è l'ordine di grandezza corretto.

$r=2$ i dimensionali:

Ora generiamo il vettore successivo nel nostro sottospazio attraverso l'operazione con $A$ sul vettore precedente:

Ora sottraiamo la proiezione di questo vettore sul vettore precedente per garantire l'ortogonalità.

Se non è già normalizzato, normalizzarlo. In questo caso, il vettore era già normalizzato, quindi

Ora proiettiamo la nostra matrice A sul sottospazio di questi due vettori:

Rimane il problema di determinare gli autovalori di questa matrice. Ma questa matrice è leggermente più piccola della matrice completa. Nei problemi che coinvolgono matrici molto grandi, lavorare con questo sottospazio più piccolo può essere molto vantaggioso.

Sebbene non sia ancora una buona stima, è migliore di quella di ordine zero. Lo faremo per un'altra iterazione, per assicurarci che il processo sia chiaro. In questo modo, però, si vanifica l'utilità del metodo, poiché si finisce per diagonalizzare una matrice 3x3 nell'iterazione successiva, il che significa che non si è risparmiato né tempo né potenza di calcolo.

$r=3$ i dimensionali:

Ora generiamo il vettore successivo nel nostro sottospazio attraverso l'operazione con A sul vettore precedente:

Ora sottraiamo la proiezione di questo vettore su entrambi i vettori precedenti per garantire l'ortogonalità.

Se non è già normalizzato, normalizzarlo. In questo caso, il vettore era già normalizzato, quindi

Ora proiettiamo la nostra matrice $A$ sul sottospazio di questi vettori:

Determiniamo ora gli autovalori:

Questi autovalori sono esattamente gli autovalori della matrice originale $A$. Questo deve essere il caso, poiché abbiamo espanso il nostro sottospazio di Krylov per coprire l'intero spazio vettoriale della matrice originale $A$.

In questo esempio, il metodo di Krylov non sembra particolarmente più semplice della diagonalizzazione diretta. Infatti, come vedremo nelle sezioni successive, il metodo di Krylov è vantaggioso solo al di sopra di una certa dimensione della matrice; ciò ha lo scopo di aiutarci a risolvere problemi di autovalori/autovettori di matrici estremamente grandi.

Questo è l'unico esempio che mostreremo lavorato "a mano", ma la sezione 2 mostra esempi computazionali.

Chiarimento dei termini

Un'idea sbagliata comune è che esista un solo sottospazio di Krylov per un determinato problema. Ma naturalmente, poiché esistono molti vettori iniziali a cui applicare la nostra matrice, esistono molti possibili sottospazi di Krylov. Useremo l'espressione "il sottospazio di Krylov" solo per riferirci a uno specifico sottospazio di Krylov già definito per un esempio specifico. Per gli approcci generali alla risoluzione dei problemi ci riferiremo a "un sottospazio di Krylov". Un'ultima precisazione: è valido il riferimento a uno " spazio di Krylov". Spesso viene chiamato " sottospazio di Krylov" per il suo utilizzo nel contesto della proiezione di matrici da uno spazio iniziale in un sottospazio. In linea con questo contesto, in questa sede ci riferiremo a questo spazio come a un sottospazio.

Verifica la tua comprensione

Spiegare perché non è (a) utile e (b) possibile estendere la dimensione del sottospazio di Krylov $r$ oltre la dimensione $N$ della matrice di interesse.

• Risposta

  (a) Poiché ortonormalizziamo i vettori man mano che li generiamo, un insieme di $N$ i di tali vettori formerà una base completa, il che significa che una combinazione lineare di essi può essere utilizzata per ricreare qualsiasi vettore dello spazio.

  (b) Il processo di ortogonalizzazione consiste nel sottrarre la proiezione di un nuovo vettore su tutti i vettori precedenti. Se tutti i vettori precedenti generano l'intero spazio vettoriale, sottraendo le proiezioni sul sottospazio completo otterremo sempre un vettore nullo.

Supponiamo che un collega ricercatore stia illustrando il metodo di Krylov applicato a una piccola matrice di esempio. C'è qualcosa che non va nella loro scelta della matrice $A$ e del vettore iniziale $|\psi\rangle$?

e

• Risposta

  Il vostro collega ha accidentalmente scelto un autovettore per il suo vettore iniziale. Agendo con la matrice sul vettore iniziale, si restituirà semplicemente lo stesso vettore, scalato dell'autovalore. Questo non genererà un sottospazio di dimensione crescente. Consigliare al collega di scegliere un vettore iniziale diverso, assicurandosi che non sia un autovettore.

Applicare il metodo di Krylov alla matrice data, scegliendo un nuovo vettore iniziale adeguato. Indica le stime dell'autovalore minimo al 0° e al 1° ordine del tuo sottospazio di Krylov.

• Risposta

  Ci sono molte risposte possibili a seconda della scelta del vettore iniziale. Sceglieremo:

  $$|v_0\rangle=\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix}.$$

  Per ottenere $|v_1\rangle$ applichiamo una volta $A$ a $|v_0\rangle$, quindi rendiamo $|v_1\rangle$ ortogonale a $|v_0\rangle$.

  $$A|v_0\rangle=\begin{pmatrix}1&1&0\\1&1&1\\0&1&1\end{pmatrix}\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix} = \frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix}$$$$A|v_0\rangle - \langle v_0|A|v_0\rangle |v_0\rangle=\frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix} - \frac{1}{\sqrt{3}}\begin{pmatrix}1&1&1\end{pmatrix}\frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix}\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix} = \frac{1}{\sqrt{3}}\begin{pmatrix}2\\3\\2\end{pmatrix}-\frac{7}{3}\frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix}=\sqrt{\frac{3}{2}}\begin{pmatrix}-1/3\\2/3\\-1/3\end{pmatrix}$$

  All'ordine 0, la proiezione sul nostro sottospazio di Krylov è

  $$\langle v_0|A|v_0\rangle=\frac{1}{\sqrt{3}}\begin{pmatrix}1&1&1\end{pmatrix} \begin{pmatrix}1&1&0\\1&1&1\\0&1&1\end{pmatrix} \frac{1}{\sqrt{3}}\begin{pmatrix}1\\1\\1\end{pmatrix} = \frac{7}{3}$$

  Al 1° ordine, la proiezione su questo sottospazio di Krylov è

  $$\langle V^1|A|V^1\rangle=\begin{pmatrix}\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}\\-\sqrt{\frac{1}{6}}&\sqrt{\frac{2}{3}}&-\sqrt{\frac{1}{6}}\end{pmatrix} \begin{pmatrix}1&1&0\\1&1&1\\0&1&1\end{pmatrix} \begin{pmatrix}\frac{1}{\sqrt{3}}&-\sqrt{\frac{1}{6}}\\\frac{1}{\sqrt{3}}& \sqrt{\frac{2}{3}} \\ \frac{1}{\sqrt{3}}&-\sqrt{\frac{1}{6}}\end{pmatrix}$$

  Questa operazione può essere eseguita a mano, ma è più facile da realizzare utilizzando numpy:

  import numpy as np
  vstar = np.array([[1/np.sqrt(3),1/np.sqrt(3),1/np.sqrt(3)],[-1/np.sqrt(6),np.sqrt(2/3),-1/np.sqrt(6)]]
  )
  A = np.array([[1, 1, 0],
                [1, 1, 1],
                [0, 1, 1]])
  v = np.array([[1/np.sqrt(3),-1/np.sqrt(6)],[1/np.sqrt(3),np.sqrt(2/3)],[1/np.sqrt(3),-1/np.sqrt(6)]])
  proj = vstar@A@v
  print(proj)
  eigenvalues, eigenvectors = np.linalg.eig(proj)
  print("The eigenvalues are ", eigenvalues)
  print("The eigenvectors are ", eigenvectors)

  output:

  [[ 2.33333333  0.47140452]
   [ 0.47140452 -0.33333333]]
  The eigenvalues are  [ 2.41421356 -0.41421356]
  The eigenvectors are  [[ 0.98559856 -0.16910198]
   [ 0.16910198  0.98559856]]

  La stima dell'autovalore minimo è -0.414.

1.2 Tipi di metodi di Krylov

i "metodi dei sottospazi di Krylov" possono riferirsi a una qualsiasi delle numerose tecniche iterative utilizzate per risolvere grandi sistemi lineari e problemi di autovalori. Ciò che hanno in comune è che costruiscono una soluzione approssimativa da un sottospazio di Krylov

$\mathcal{K}^r(A,|v\rangle ) = \text{span}\{|v\rangle, A|v\rangle, A^2|v\rangle, ..., A^{r-1}|v\rangle\},$

dove $|v\rangle$ è l'ipotesi iniziale (cfr. Rif. [5] ). Si differenziano per il modo in cui scelgono la migliore approssimazione da questo sottospazio, bilanciando fattori come la velocità di convergenza, l'uso della memoria e il costo computazionale complessivo. L'obiettivo di questa lezione è quello di sfruttare il calcolo quantistico nel contesto dei metodi dei sottospazi di Krylov; una discussione esaustiva di questi metodi esula dal suo scopo. Le brevi definizioni riportate di seguito servono solo come contesto e includono alcuni riferimenti per approfondire questi metodi.

Il metodo del gradiente coniugato (CG) : Questo metodo viene utilizzato per risolvere sistemi lineari simmetrici e definiti positivamente [6]. Minimizza la norma A dell'errore a ogni iterazione, il che lo rende particolarmente efficace per i sistemi derivanti da PDE ellittiche discretizzate [7]. Nella prossima sezione utilizzeremo questo approccio per motivare il motivo per cui un sottospazio di Krylov sarebbe un sottospazio efficace in cui sondare le soluzioni migliori dei sistemi lineari.

Il metodo del residuo minimo generalizzato (GMRES) : È progettato per risolvere sistemi lineari generali non simmetrici. Minimizza la norma residua su uno spazio di Krylov a ogni iterazione, il che lo rende robusto ma potenzialmente impegnativo in termini di memoria per sistemi di grandi dimensioni [7].

Il metodo del residuo minimo (MINRES) : Questo metodo viene utilizzato per risolvere sistemi lineari indefiniti simmetrici. È simile a GMRES, ma sfrutta la simmetria della matrice per ridurre il costo computazionale [8].

Altri approcci degni di nota includono il metodo di ortogonalizzazione completa (FOM), che è strettamente correlato al metodo di Arnoldi per i problemi agli autovalori, il metodo del gradiente biconiugato ( BiCG ) e il metodo di riduzione dimensionale indotta (IDR).

1.3 Perché funziona il metodo del sottospazio di Krylov

In questa sede motiveremo che il metodo del sottospazio di Krylov dovrebbe essere un modo efficiente per approssimare gli autovalori delle matrici attraverso il perfezionamento iterativo delle approssimazioni degli autovettori, attraverso la lente della discesa più ripida. Si sosterrà che, data un'ipotesi iniziale di stato fondamentale, lo spazio delle correzioni successive a tale ipotesi iniziale che produce la convergenza più rapida è un sottospazio di Krylov. Ci fermiamo a una prova rigorosa del comportamento di convergenza.

Supponiamo che la matrice di interesse $A$ sia simmetrica e definita positivamente. Ciò rende la nostra argomentazione più pertinente al metodo CG di cui sopra. Non facciamo alcuna ipotesi sulla sparsità, né sosteniamo che $A$ debba essere un'Hermitiana (che deve essere tale se è un'Hamiltoniana).

In genere si vuole risolvere un problema della forma

Si potrebbe immaginare che $|b\rangle=c|x\rangle$ dove $c$ è una costante, come in un problema di autovalori. Ma il nostro problema rimane più generale per ora.

Partiamo da un vettore $|x_0\rangle$ che rappresenta una soluzione approssimativa. Sebbene vi siano dei parallelismi tra questa ipotesi $|x_0\rangle$ e $|v_0\rangle$ nella sezione 1.1, in questa sede non li stiamo prendendo in considerazione. A nostro avviso, l' $|x_0\rangle$ e presenta un errore, che chiamiamo “ $|e_0\rangle$ ”:

Definiamo inoltre l' $R_0$ residuale:

Qui usiamo la maiuscola $R$ per distinguere il residuo dalla dimensione del nostro sottospazio di Krylov $r$.

Si vuole ora effettuare un passo di correzione della forma

che speriamo migliori la nostra approssimazione. Qui $|p_0\rangle$ è un vettore ancora da determinare. Sia $|e_1\rangle$ l'errore dopo la correzione. Allora

Ci interessa sapere come si comporta il nostro errore quando viene trasformato dalla nostra matrice. Calcoliamo quindi la $A$ -norm dell'errore. Cioè

dove abbiamo sfruttato la simmetria di $A$ e anche il fatto che $A |e_0\rangle = |R_0\rangle$. Qui $d$ è una costante indipendente da $|p_0\rangle$. Come menzionato nella Sezione 1.2, la norma $A$ dell'errore non è l'unica grandezza che potremmo scegliere di minimizzare, ma è una buona scelta. Vogliamo vedere come questa grandezza vari in funzione della scelta dei vettori di correzione $|p_0\rangle$. Definiamo quindi la funzione $f$ ponendo

$f$ è solo l'errore $|e_1\rangle$ in funzione della correzione $|p_0\rangle$ misurata nella norma $A$. Pertanto, vogliamo scegliere $|p_0\rangle$ in modo che $f(|p_0\rangle)$ sia il più piccolo possibile. A questo scopo, calcoliamo il gradiente di $f$. Utilizzando la simmetria di $A$ abbiamo

Il gradiente indica la direzione di massima crescita; ciò significa che la sua opposta ci indica la direzione in cui la funzione diminuisce maggiormente: la direzione di massima diminuzione. Partendo dalla nostra ipotesi iniziale $|x_0\rangle$, dove $|p_0\rangle=0$, si ha che $\nabla f(0) = -2|R_0\rangle$. Pertanto, la funzione $f$ diminuisce maggiormente nella direzione del residuo $|R_0\rangle$. Quindi la nostra scelta iniziale trarrebbe il massimo beneficio dall’aggiunta del vettore $|p_0\rangle=\alpha_0 |R_0\rangle$ per un certo scalare $\alpha_0$.

Nella fase successiva, scegliamo nuovamente un vettore $|p_1\rangle$ e aggiungiamo il suo valore all'approssimazione corrente. Utilizzando lo stesso ragionamento di prima, scegliamo $|p_1\rangle = \alpha_1 |R_1\rangle$ per un qualche scalare $\alpha_1$. Continuiamo in questo modo, in modo tale che l'iterazione $k^\text{th}$ del nostro vettore sia

In modo equivalente, vogliamo costruire lo spazio da cui scegliere le nostre stime migliorate aggiungendo $|R_0\rangle$, $|R_1\rangle$, e così via, in ordine. Il vettore stimato da $k^\text{th}$ si trova in

Ora, utilizzando la relazione che

vediamo che

In altre parole, lo spazio che costruiamo e che si avvicina nel modo più efficiente alla soluzione corretta $|x\rangle$ è esattamente lo spazio costruito mediante l’applicazione successiva della matrice $A$ su $|R_0\rangle$. Un sottospazio di Krylov è lo spazio generato dai vettori delle direzioni successive di discesa più ripida.

Infine, ribadiamo che non abbiamo fatto alcuna affermazione numerica sulla scalabilità di questo approccio, né abbiamo discusso i vantaggi comparativi per le matrici rade. Questo ha solo lo scopo di motivare l'uso dei metodi dei sottospazi di Krylov e di aggiungere un senso intuitivo a questi metodi. Ora esploreremo il comportamento di questi metodi dal punto di vista numerico.

Verifica la tua comprensione

Nel flusso di lavoro precedente, abbiamo proposto di minimizzare la $A$ -norm dell'errore. Quali altre grandezze si potrebbero minimizzare nella ricerca dello stato fondamentale e del suo autovalore?

• Risposta

  Si potrebbe immaginare di utilizzare il vettore dei residui invece della $A$ -norm dell'errore. Ci possono essere casi in cui è utile considerare il vettore di errore stesso.

2. Metodi di Krylov nel calcolo classico

In questa sezione implementiamo le iterazioni di Arnoldi dal punto di vista computazionale, in modo da poter sfruttare un sottospazio di Krylov nella risoluzione di problemi di autovalori. Applicheremo questo metodo dapprima a un esempio su piccola scala, quindi esamineremo come il tempo di calcolo scala all'aumentare delle dimensioni della matrice di interesse. L'idea chiave è che la generazione dei vettori che coprono lo spazio di Krylov contribuirà notevolmente al tempo di calcolo totale richiesto. La memoria richiesta varia a seconda dei metodi di Krylov. Ma i vincoli di memoria possono limitare l'uso dei metodi di Krylov tradizionali.

2.1 Semplice esempio su piccola scala

Nel processo di creazione di un sottospazio di Krylov, dovremo ortonormalizzare i vettori nel nostro sottospazio. Definiamo una funzione che prenda un vettore stabilito dal nostro sottospazio vknown (che non si presume sia normalizzato) e un vettore candidato da aggiungere al nostro sottospazio vnext e che renda vnext ortogonale a vknown e normalizzato. Definiamo inoltre una funzione che esegue questo processo per tutti i vettori stabiliti nel nostro sottospazio di Krylov per garantire un insieme completamente ortonormale.

# vknown is some established vector in our subspace. vnext is one we wish to add,
# which must be orthogonal to vknown.


def orthog_pair(vknown, vnext):
    vknown = vknown / np.sqrt(vknown.T @ vknown)
    diffvec = vknown.T @ vnext * vknown
    vnext = vnext - diffvec
    return vnext


# v is the candidate vector to be added to our subspace. s is the existing subspace.


def orthoset(v, s):
    v = v / np.sqrt(v.T @ v)
    temp = v
    for i in range(len(s)):
        temp = orthog_pair(s[i], temp)
    v = temp / np.sqrt(temp.T @ temp)
    return v

Definiamo ora una funzione che costruisce un sottospazio di Krylov iterativamente sempre più grande, fino a quando lo spazio dei vettori di Krylov si estende all'intero spazio della matrice originale. Questo ci permetterà di vedere quanto gli autovalori ottenuti con il nostro metodo del sottospazio di Krylov corrispondano ai valori esatti, in funzione della dimensione del sottospazio di Krylov. È importante notare che la nostra funzione krylov_full_build restituisce i vettori di Krylov, gli hamiltoniani proiettati, gli autovalori e il tempo richiesto.

# Necessary imports and definitions to track time in microseconds
import time


def time_mus():
    return int(time.time() * 1000000)


# This function constructs a Krylov subspace that spans the whole space of the original matrix.
#     Input:
#       v0          : initial vector
#       matrix      : original matrix to be diagonalized
#     Output:
#       ks          : Krylov vectors
#       Hs          : projected Hamiltonians
#       eigs        : eigenvalues
#       k_tot_times : time required for the operation


def krylov_full_build(v0, matrix):
    t0 = time_mus()
    b = v0 / np.sqrt(v0 @ v0.T)
    A = matrix
    ks = []
    ks.append(b)
    Hs = []
    eigs = []
    Hs.append(b.T @ A @ b)
    eigs.append(np.array([b.T @ A @ b]))
    k_tot_times = []

    for j in range(len(A) - 1):
        vec = A @ ks[j].T
        ortho = orthoset(vec, ks)
        ks.append(ortho)
        ksarray = np.array(ks)
        Hs.append(ksarray @ A @ ksarray.T)
        eigs.append(np.linalg.eig(Hs[j + 1]).eigenvalues)
        k_tot_times.append(time_mus() - t0)

    # Return the Krylov vectors, the projected Hamiltonians, the eigenvalues,
    # and the total time required.
    return (ks, Hs, eigs, k_tot_times)

Lo proveremo su una matrice ancora piuttosto piccola, ma più grande di quella che si potrebbe fare a mano.

# Define our small test matrix
test_matrix = np.array(
    [
        [4, -1, 0, 1, 0],
        [-1, 4, -1, 2, 1],
        [0, -1, 4, 3, 3],
        [1, 2, 3, 4, 0],
        [0, 1, 3, 0, 4],
    ]
)

# Give the test matrix and an initial guess as arguments in the function defined above.
# Calculate outputs.
test_ks, test_Hs, test_eigs, text_k_tot_times = krylov_full_build(
    np.array([0.5, 0.5, 0, 0.5, 0.5]), test_matrix
)

Possiamo verificare le nostre funzioni assicurandoci che nell'ultimo passo (quando lo spazio di Krylov è lo spazio vettoriale completo della matrice originale) gli autovalori del metodo di Krylov corrispondano esattamente a quelli della diagonalizzazione numerica esatta:

print(np.linalg.eig(test_matrix).eigenvalues)
print(test_eigs[len(test_matrix) - 1])

Output:

[-1.36956923  8.43756009  2.9040308   5.34436028  4.68361806]
[-1.36956923  8.43756009  2.9040308   4.68361806  5.34436028]

È stato un successo. Naturalmente, ciò che conta davvero è la bontà della nostra approssimazione in funzione della dimensione del nostro sottospazio di Krylov. Poiché spesso ci preoccupiamo di trovare stati fondamentali e altri autovalori minimi (e per altre ragioni più algebriche spiegate più avanti), esaminiamo la nostra stima dell'autovalore più basso in funzione della dimensione del sottospazio di Krylov. Cioè

def errors(matrix, krylov_eigs):
    targ_min = min(np.linalg.eig(matrix).eigenvalues)
    err = []
    for i in range(len(matrix)):
        err.append(min(krylov_eigs[i]) - targ_min)
    return err

import matplotlib.pyplot as plt

krylov_error = errors(test_matrix, test_eigs)

plt.plot(krylov_error)
plt.axhline(y=0, color="red", linestyle="--")  # Add dashed red line at y=0
plt.xlabel("Order of Krylov subspace")  # Add x-axis label
plt.ylabel("Error in minimum eigenvalue")  # Add y-axis label
plt.show()

Output:

Si osserva che l’autovalore minimo viene raggiunto con discreta precisione una volta che il sottospazio di Krylov ha raggiunto una dimensione pari a $\mathcal{K}^2$ ed è perfetto quando $\mathcal{K}^3$.

2.2 Scalatura temporale con dimensione della matrice

Per convincerci che il metodo di Krylov può essere vantaggioso rispetto agli eigensolver numerici esatti, procediamo nel modo seguente:

• Costruire matrici casuali (non rade, non l'applicazione ideale per KQD)
• Determinare gli autovalori con due metodi: direttamente con NumPy e con un sottospazio di Krylov.
• Scegliamo un limite per la precisione degli autovalori prima di accettare le stime di Krylov.
• Confrontate il tempo di parete richiesto per risolvere questi due modi.

Avvertenze: Come vedremo in dettaglio più avanti, la diagonalizzazione quantistica di Krylov si applica meglio agli operatori le cui rappresentazioni matriciali sono rade e/o possono essere scritte utilizzando un piccolo numero di gruppi di operatori di Pauli commutativi. Le matrici casuali che stiamo utilizzando non rientrano in questa descrizione. Questi sono utili solo per sondare la scala alla quale i metodi di Krylov classici potrebbero essere utili. In secondo luogo, utilizzando il metodo di Krylov, calcoleremo gli autovalori utilizzando molti sottospazi di Krylov di dimensioni diverse. Riporteremo il tempo necessario per ottenere il sottospazio di Krylov di dimensione minima che raggiunge l'accuratezza richiesta per l'autovalore dello stato fondamentale. Anche in questo caso, si tratta di una soluzione un po' diversa da quella di un problema intrattabile per gli eigensolver esatti, poiché stiamo usando la soluzione esatta per valutare la dimensione necessaria.

Iniziamo generando il nostro insieme di matrici casuali.

import numpy as np

# Set the random seed
np.random.seed(42)

# how many random matrices will we make
num_matrix = 200

matrices = []
for m in range(1, num_matrix):
    matrices.append(np.random.rand(m, m))

Ora diagonalizziamo ogni matrice direttamente, utilizzando numpy. Per un successivo confronto, calcoliamo il tempo necessario per la diagonalizzazione.

matrix_numpy_times = []
matrix_numpy_eigs = []
for mm in range(num_matrix - 1):
    t0 = time_mus()
    matrix_numpy_eigs.append(min(np.linalg.eig(matrices[mm]).eigenvalues))
    matrix_numpy_times.append(time_mus() - t0)

plt.plot(matrix_numpy_times)
plt.xlabel("Dimension of matrix")  # Add x-axis label
plt.ylabel("Time to diagonalize (microsec)")  # Add y-axis label
plt.show()

Output:

Si noti che nell'immagine precedente, il tempo anomalo intorno a una dimensione di 125 può essere dovuto alla natura casuale delle matrici o all'implementazione sul processore classico utilizzato, ma non è riproducibile. Eseguendo nuovamente il codice si otterrà un profilo diverso con diversi picchi anomali.

Ora per ogni matrice costruiremo un sottospazio di Krylov e calcoleremo gli autovalori per gradi. A ogni passo, controlleremo se l'autovalore più basso è stato ottenuto entro l'errore assoluto specificato. Il sottospazio che per primo ci fornisce autovalori entro l'errore specificato è il sottospazio per il quale registreremo i tempi di calcolo. L'esecuzione di questa cella può richiedere alcuni minuti, a seconda della velocità del processore. È possibile saltare la valutazione o ridurre la dimensione massima delle matrici diagonalizzate. È sufficiente guardare i risultati precalcolati.

# Choose the absolute error you can tolerate, and make a list for tracking the Krylov subspace size
# at which that error is achieved.
abserr = 0.05
accept_subspace_size = []

# Lists to store total time spent on the Krylov method, and the subset of that time spent on
# diagonalizing the projected matrix.
matrix_krylov_tot_times = []
matrix_krylov_dim = []

# Step through all our random matrices
for mm in range(0, num_matrix - 1):
    test_ks, test_Hs, test_eigs, test_k_tot_times = krylov_full_build(
        np.ones(len(matrices[mm])), matrices[mm]
    )
    # We have not yet found a Krylov subspace that produces our minimum eigenvalue to
    # within the required error.
    found = 0
    for j in range(0, len(matrices[mm]) - 1):
        # If we still haven't found the desired subspace...
        if found == 0:
            # ...but if this one satisfies the requirement, then record everything
            if (
                abs((min(test_eigs[j]) - matrix_numpy_eigs[mm]) / matrix_numpy_eigs[mm])
                < abserr
            ):
                accept_subspace_size.append(j)
                matrix_krylov_tot_times.append(test_k_tot_times[j])
                matrix_krylov_dim.append(mm)
                found = 1

Tracciamo per confronto i tempi ottenuti con questi due metodi:

plt.plot(matrix_numpy_times, color="blue")
plt.plot(matrix_krylov_dim, matrix_krylov_tot_times, color="green")
plt.xlabel("Dimension of matrix")  # Add x-axis label
plt.ylabel("Time to diagonalize (microsec)")  # Add y-axis label
plt.show()

Output:

Questi sono i tempi effettivi richiesti, ma ai fini della discussione, smussiamo queste curve calcolando la media su alcuni punti adiacenti / dimensioni della matrice. Questo viene fatto di seguito:

smooth_numpy_times = []
smooth_krylov_times = []

# Choose the number of adjacent points over which to average forward;
# the same will be used backward.
smooth_steps = 10

# We will do this smoothing for all points/matrix dimensions
for i in range(len(matrix_krylov_tot_times)):
    # Ensure we don't exceed the boundaries of our lists
    start = max(0, i - smooth_steps)
    end = min(len(matrix_krylov_tot_times) - 1, i + smooth_steps)

    # Dummy variables for accumulating an average over adjacent points. This is done for both Krylov
    # and the NumPy calculations.
    smooth_count = 0
    smooth_numpy_sum = 0
    smooth_krylov_sum = 0

    for j in range(start, end):
        smooth_numpy_sum = smooth_numpy_sum + matrix_numpy_times[j]
        smooth_krylov_sum = smooth_krylov_sum + matrix_krylov_tot_times[j]
        smooth_count = smooth_count + 1

    # Appending the averaged adjacent values to our new smooth lists
    smooth_numpy_times.append(smooth_numpy_sum / smooth_count)
    smooth_krylov_times.append(smooth_krylov_sum / smooth_count)

plt.plot(smooth_numpy_times, color="blue")
plt.plot(smooth_krylov_times, color="green")
plt.xlabel("Dimension of matrix")  # Add x-axis label
plt.ylabel("Time to diagonalize (smoothed, microsec)")  # Add y-axis label
plt.show()

Output:

Si noti che il tempo richiesto per la costruzione di un sottospazio di Krylov inizialmente supera il tempo richiesto per la diagonalizzazione completa di numpy. Ma quando le dimensioni della matrice crescono, il metodo di Krylov diventa vantaggioso. Ciò è vero anche se riduciamo l'errore accettabile, ma il vantaggio si manifesta con una matrice di dimensioni maggiori. Vale la pena di analizzare questo aspetto.

La complessità temporale della diagonalizzazione numerica è $O(n^3)$ (con qualche variazione tra gli algoritmi). Anche la complessità temporale della generazione di una base ortonormale di $n$ vettori è $O(n^3)$. Quindi il vantaggio del metodo di Krylov non è legato all'uso di $\it{some}$ basi ortonormali, ma all'uso di una particolare base ortonormale che seleziona efficacemente gli autovalori di interesse. Lo abbiamo già visto nello schizzo di una prova nella prima sezione di questa lezione, e questo è fondamentale per le garanzie di convergenza dei metodi di Krylov.

Ripercorriamo i progressi compiuti finora:

• Per matrici molto grandi, il metodo del sottospazio di Krylov può produrre autovalori approssimati entro le tolleranze richieste, più velocemente degli algoritmi di diagonalizzazione tradizionali.
• Per matrici così grandi, la generazione di un sottospazio di Krylov è la parte più lunga del metodo dei sottospazi di Krylov.
• Pertanto, un modo efficiente di generare un sottospazio di Krylov sarebbe molto utile. È qui che entra in gioco il computer quantistico.

Verifica la tua comprensione

Si veda il grafico livellato dei tempi di diagonalizzazione in funzione della dimensione della matrice riportato sopra.

(a) In base a questo grafico, a partire da quale dimensione approssimativa della matrice il metodo di Krylov ha iniziato a essere più veloce?

(b) Quali aspetti del calcolo potrebbero modificare la dimensione a partire dalla quale il metodo di Krylov diventa più veloce?

• Risposta

  (a) I risultati possono variare se si esegue nuovamente il calcolo, ma il metodo di Krylov diventa più veloce a partire da una dimensione di circa 80-85.

  (b) Le risposte possibili sono molte. Alcuni fattori importanti sono il livello di precisione richiesto e la sparsità delle matrici da diagonalizzare.

3. Krylov tramite evoluzione temporale

Tutto ciò che abbiamo descritto finora può essere fatto in modo classico. Come e quando potremmo usare un computer quantistico? Per matrici molto grandi, il metodo di Krylov può richiedere lunghi tempi di calcolo e grandi quantità di memoria. Il tempo richiesto per l'operazione matriciale di $H$ su $|v\rangle$ scala come $O(N^2)$ nel caso peggiore. Anche la moltiplicazione di matrici rade su un vettore (il caso tipico dei solutori classici di tipo Krylov) ha una complessità temporale che scala come $O(N)$. Questo viene fatto per ogni vettore che vogliamo nel nostro sottospazio. La dimensione del sottospazio $r$ di solito non è una frazione significativa di $N$, e spesso scala come $\log(N)$. Quindi, nel caso peggiore, la generazione di tutti i vettori scala come $O(N^2 \log(N))$. Anche se ci sono altri passaggi, come l'ortogonalizzazione, questa è la scalatura dominante da tenere a mente.

L'informatica quantistica ci permette di cambiare gli attributi del problema che determinano la scalabilità del tempo e delle risorse necessarie. Invece della dipendenza dalla dimensione della matrice $N$ in generale, vedremo cose come il numero di scatti e il numero di termini di Pauli non commutativi che compongono l'Hamiltoniana. Vediamo come funziona.

3.1 Evoluzione temporale

Ricordiamo che l'operatore che fa evolvere nel tempo uno stato quantistico è $e^{-iHt/\hbar}$ (ed è molto comune, soprattutto nell'informatica quantistica, eliminare $\hbar$ dalla notazione). Un modo per comprendere e persino realizzare una funzione esponenziale di un operatore è osservare la sua espansione in serie di Taylor. Si noti che questa operazione che agisce su un vettore iniziale $|v\rangle$ produce una somma di termini con potenze crescenti di $H$ applicate allo stato iniziale. Sembra che possiamo creare il nostro sottospazio di Krylov facendo evolvere nel tempo il nostro stato di ipotesi iniziale!

Il problema sta nel realizzare l'evoluzione temporale su un vero computer quantistico. Molti dei termini dell'Hamiltoniana non si scambiano tra loro. Quindi, mentre alcuni semplici operatori esponenziali come $e^{-iZ}$ corrispondono a circuiti semplici, gli hamiltoniani generali non lo sono. E poiché contengono termini non commutativi, non possiamo semplicemente decomporre l'esponenziale in un prodotto di termini semplici, come possiamo fare con i numeri.

Non si tratta quindi di un processo banale, ma è un processo ben studiato nell'informatica quantistica. Realizziamo l'evoluzione del tempo sui computer quantistici utilizzando un processo chiamato trotterizzazione, che di per sé è un argomento molto ricco [10]. Ma a un livello molto alto, spezzando l'evoluzione temporale in passi molto piccoli, ad esempio $m$ passi di dimensione $dt$, limitiamo gli effetti della non commutatività dei termini.

dove $dt = t/m$.

Chiamiamo un sottospazio di Krylov di ordine r che abbiamo generato nel contesto classico usando direttamente potenze di H un "sottospazio di Krylov di potenza".

Ora generiamo uno spazio simile usando l'operatore unitario di evoluzione temporale $U \equiv e^{-iHt}$; ci riferiremo a questo come allo "spazio unitario di Krylov" $\mathcal{K}_U^r$. Il sottospazio di Krylov di potenza $\mathcal{K}_P^r$ che usiamo classicamente non può essere generato direttamente su un computer quantistico, poiché $H$ non è un operatore unitario. È possibile dimostrare che l'utilizzo del sottospazio unitario di Krylov fornisce garanzie di convergenza simili a quelle del sottospazio di Krylov di potenza, ossia che l'errore di stato a terra converge in modo efficiente finché lo stato iniziale $|v\rangle$ ha una sovrapposizione con lo stato a terra vero che non è esponenzialmente svanente e finché c'è un intervallo sufficiente tra gli autovalori. Per una discussione più precisa sulla convergenza, si veda Ref [1].

Qui le potenze di $U$ diventano passi temporali diversi (la potenza $k^\text{th}$ di $U$ avanza di un tempo $k \times dt$ ). Possiamo etichettare l'elemento del sottospazio che si evolve nel tempo per il tempo totale $k dt$ come $|\psi_k\rangle$.

Possiamo proiettare la nostra hamiltoniana H sul sottospazio unitario di Krylov, $\mathcal{K}_U^r$. In altre parole, calcoliamo ogni elemento della matrice di $H$ nella base di $\mathcal{K}_U^r$. Ci riferiremo a questa matrice proiettata come $\tilde{H}$.

3.2 Come implementare su un computer quantistico

Gli elementi della matrice $\tilde{H}$ sono dati dai valori di aspettativa $\langle \psi_m |H| \psi_n\rangle$, che possono essere stimati utilizzando il computer quantistico. Si tenga presente che $H$ può essere scritto come una somma di operatori Pauli su diversi qubit e che non tutti gli operatori Pauli possono essere misurati contemporaneamente. Possiamo ordinare i termini di Pauli in gruppi di termini commutativi e misurarli tutti insieme. Ma potrebbero essere necessari molti gruppi di questo tipo per coprire tutti i termini. Diventa quindi importante il numero di gruppi pendolari distinti in cui i termini possono essere suddivisi, $N_\text{GCP}$.

Qui $P_\alpha$ indica una stringa di Pauli della forma $P_\alpha \sim IZIXII...YZXIX$ oppure un insieme di tali stringhe di Pauli che commutano tra loro. Dato che possiamo scrivere $H$ come somma di operatori misurabili, le seguenti espressioni per gli elementi di matrice di $\tilde{H}$ possono essere realizzate utilizzando l’Estimator primitivo IBM Quantum.

Dove $\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle$ sono i vettori dello spazio unitario di Krylov e $t_n = n dt$ sono i multipli del passo temporale $dt$ scelto. Su un computer quantistico, il calcolo di ciascun elemento della matrice può essere effettuato con qualsiasi algoritmo che consenta di ottenere una sovrapposizione tra stati quantistici. In questa lezione ci concentreremo sul test di Hadamard. Dato che l’ $\mathcal{K}_U$ e ha dimensione $r$, l’Hamiltoniano proiettato nel sottospazio avrà dimensione $r \times r$. Se $r$ è sufficientemente piccolo (in genere $r<<100$ è sufficiente per ottenere la convergenza delle stime degli autovalori), possiamo quindi facilmente diagonalizzare l’Hamiltoniano proiettato $\tilde{H}$, in termini classici. Tuttavia, non è possibile diagonalizzare direttamente l' $\tilde{H}$ e 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 le soluzioni di questo problema generalizzato degli autovalori. Ad esempio, la stima dell'energia dello stato fondamentale si ottiene prendendo l'autovalore più piccolo $E$ 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$.

Problema generalizzato degli autovalori

Perché non possiamo semplicemente diagonalizzare $\tilde{H}$? Poiché $\tilde{S}$ contiene le informazioni sulla geometria della base di Krylov (che è non ortogonale in tutti i casi tranne quelli molto speciali), $\tilde{H}$ da solo non descrive una proiezione dell'hamiltoniana completa, quindi i suoi autovalori non hanno alcuna relazione particolare con quelli dell'hamiltoniana completa: potrebbero essere qualsiasi valore casuale. La soluzione del problema degli autovalori generalizzati è necessaria per ottenere gli autovalori e gli autovettori approssimati corrispondenti alla proiezione dell'hamiltoniana completa nello spazio di Krylov.

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$. L'implementazione di questa operazione su un computer quantistico è illustrata in dettaglio qui di seguito.

4. Diagonalizzazione quantistica di Krylov su un computer quantistico

Ora implementeremo la diagonalizzazione quantistica di Krylov su un computer quantistico reale. Iniziamo importando alcuni pacchetti utili.

import numpy as np
import scipy as sp
import matplotlib.pylab as plt
from typing import Union, List
import warnings

from qiskit.quantum_info import SparsePauliOp, Pauli
from qiskit.circuit import Parameter
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter

# from qiskit.providers.fake_provider import Fake20QV1
from qiskit_ibm_runtime import QiskitRuntimeService, EstimatorV2 as Estimator, Batch

import itertools as it

warnings.filterwarnings("ignore")

Definiamo la funzione che segue per risolvere il problema generalizzato degli autovalori appena spiegato sopra.

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 our Krylov subspace
        s (numpy.ndarray):
            The matrix of overlaps between vectors of our 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]

Almeno nel benchmarking iniziale, è utile conoscere una soluzione classica esatta per verificare il comportamento della convergenza. La funzione seguente calcola l'energia di stato fondamentale di un'hamiltoniana, utilizzando come argomenti l'hamiltoniana e il numero di qubit.

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)

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

    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

4.1 Fase 1: Mappare il problema su circuiti quantistici e operatori

Ora definiremo un'hamiltoniana. Questa funzione si distingue dalla funzione precedente in quanto quest'ultima prende come argomento un'hamiltoniana e restituisce solo lo stato fondamentale, e lo fa in modo classico. L'hamiltoniana che definiamo qui determina i livelli energetici di tutti gli autostati energetici e può essere costruita utilizzando gli operatori di Pauli e implementata su un computer quantistico.

Scegliamo un'hamiltoniana corrispondente a una catena di spin che può avere qualsiasi orientamento nello spazio, chiamata "catena di Heisenberg". Assumiamo che lo spin $i^\text{th}$ possa essere influenzato dai suoi vicini più prossimi (gli spin $(i-1)^\text{th}$ e $(i+1)^\text{th}$ ) ma non da quelli più lontani. Si ammette anche la possibilità che l'interazione tra gli spin sia diversa quando gli spin puntano lungo assi diversi. Questa asimmetria si verifica talvolta, ad esempio, a causa della struttura del reticolo cristallino in cui sono inseriti gli spin.

# Define problem Hamiltonian.
n_qubits = 10
# coupling strength for XX, YY, and ZZ interactions
JX = 1
JY = 3
JZ = 2

# 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, JZ)
    if term.count("Z") == 2
    else (term, JY)
    if term.count("Y") == 2
    else (term, JX)
    for term in H_int
]

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