Diagonalizarea cuantică Krylov a Hamiltonienilor de rețea

Estimare de utilizare: 20 de minute pe un Heron r2 (NOTĂ: Aceasta este doar o estimare. Durata de execuție poate varia.)

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

# This cell is hidden from users – it disables some lint rules
# ruff: noqa: E402 E722 F601

Context

Acest tutorial demonstrează cum să implementezi algoritmul de diagonalizare cuantică Krylov (KQD) în contextul tiparelor Qiskit. Vei învăța mai întâi teoria din spatele algoritmului, apoi vei vedea o demonstrație a execuției sale pe un QPU.

În multe discipline, suntem interesați să cunoaștem proprietățile stării fundamentale ale sistemelor cuantice. Exemple includ înțelegerea naturii fundamentale a particulelor și forțelor, predicția și înțelegerea comportamentului materialelor complexe, precum și înțelegerea interacțiunilor și reacțiilor biochimice. Din cauza creșterii exponențiale a spațiului Hilbert și a corelațiilor care apar în sistemele entanglate, algoritmii clasici se confruntă cu dificultăți în rezolvarea acestei probleme pentru sisteme cuantice de dimensiuni tot mai mari. La un capăt al spectrului se află abordarea existentă care valorifică hardware-ul cuantic, concentrată pe metodele cuantice variaționale (de exemplu, eigensolver-ul cuantic variațional). Aceste tehnici se confruntă cu provocări pe dispozitivele actuale din cauza numărului mare de apeluri de funcție necesare în procesul de optimizare, ceea ce adaugă un supracost mare de resurse odată ce sunt introduse tehnici avansate de mitigare a erorilor, limitând astfel eficacitatea lor la sisteme mici. La celălalt capăt al spectrului se află metodele cuantice tolerante la defecte cu garanții de performanță (de exemplu, estimarea fazei cuantice), care necesită circuite adânci ce pot fi executate doar pe un dispozitiv tolerant la defecte. Din aceste motive, introducem aici un algoritm cuantic bazat pe metode de subspațiu (descris în această lucrare de sinteză), algoritmul de diagonalizare cuantică Krylov (KQD). Acest algoritm funcționează bine la scară largă [1] pe hardware-ul cuantic existent, are garanții de performanță similare cu estimarea fazei, este compatibil cu tehnici avansate de mitigare a erorilor și ar putea furniza rezultate inaccesibile pe cale clasică.

Cerințe

Înainte de a începe acest tutorial, asigură-te că ai instalate următoarele:

• Qiskit SDK v2.0 sau versiune ulterioară, cu suport pentru vizualizare
• Qiskit Runtime v0.22 sau versiune ulterioară ( pip install qiskit-ibm-runtime )

Configurare

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)

    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

Pasul 1: Maparea intrărilor clasice la o problemă cuantică

Spațiul Krylov

Spațiul Krylov $\mathcal{K}^r$ de ordin $r$ este spațiul generat de vectorii obținuți prin înmulțirea puterilor mai mari ale unei matrice $A$, până la $r-1$, cu un vector de referință $\vert v \rangle$.

$$\mathcal{K}^r = \left\{ \vert v \rangle, A \vert v \rangle, A^2 \vert v \rangle, ..., A^{r-1} \vert v \rangle \right\}$$

Dacă matricea $A$ este Hamiltonianul $H$, vom numi spațiul corespunzător spațiul Krylov al puterilor $\mathcal{K}_P$. În cazul în care $A$ este operatorul de evoluție temporală generat de Hamiltonian $U=e^{-iHt}$, vom numi spațiul spațiu Krylov unitar $\mathcal{K}_U$. Subspațiul Krylov al puterilor pe care îl folosim clasic nu poate fi generat direct pe un calculator cuantic, deoarece $H$ nu este un operator unitar. În schimb, putem folosi operatorul de evoluție temporală $U = e^{-iHt}$, care poate fi demonstrat că oferă garanții de convergență similare cu metoda puterilor. Puterile lui $U$ devin atunci diferiți pași de timp $U^k = e^{-iH(kt)}$.

$$\mathcal{K}_U^r = \left\{ \vert \psi \rangle, U \vert \psi \rangle, U^2 \vert \psi \rangle, ..., U^{r-1} \vert \psi \rangle \right\}$$

Consultă Anexa pentru o derivare detaliată a modului în care spațiul Krylov unitar permite reprezentarea cu precizie a stărilor proprii de energie joasă.

Algoritmul de diagonalizare cuantică Krylov

Dat un Hamiltonian $H$ pe care dorim să-l diagonalizăm, considerăm mai întâi spațiul Krylov unitar corespunzător $\mathcal{K}_U$. Scopul este de a găsi o reprezentare compactă a Hamiltonianului în $\mathcal{K}_U$, pe care o vom numi $\tilde{H}$. Elementele de matrice ale $\tilde{H}$, proiecția Hamiltonianului în spațiul Krylov, pot fi calculate prin calcularea valorilor de așteptare următoare

$$\tilde{H}_{mn} = \langle \psi_m \vert H \vert \psi_n \rangle =$$ $$= \langle \psi \vert e^{i H t_m} H e^{-i H t_n} \vert \psi \rangle$$ $$= \langle \psi \vert e^{i H m dt} H e^{-i H n dt} \vert \psi \rangle$$

Unde $\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle$ sunt vectorii spațiului Krylov unitar și $t_n = n dt$ sunt multiplii pasului de timp $dt$ ales. Pe un calculator cuantic, calculul fiecărui element de matrice poate fi realizat cu orice algoritm care permite obținerea suprapunerii dintre stările cuantice. Acest tutorial se concentrează pe testul Hadamard. Dat că $\mathcal{K}_U$ are dimensiunea $r$, Hamiltonianul proiectat în subspațiu va avea dimensiunile $r \times r$. Cu $r$ suficient de mic (în general $r<<100$ este suficient pentru a obține convergența estimărilor valorilor proprii) putem apoi să diagonalizăm ușor Hamiltonianul proiectat $\tilde{H}$. Cu toate acestea, nu putem diagonaliza direct $\tilde{H}$ din cauza ne-ortogonalității vectorilor spațiului Krylov. Va trebui să măsurăm suprapunerile acestora și să construim o matrice $\tilde{S}$

$$\tilde{S}_{mn} = \langle \psi_m \vert \psi_n \rangle$$

Aceasta ne permite să rezolvăm problema valorilor proprii într-un spațiu ne-ortogonal (numită și problemă generalizată a valorilor proprii)

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

Se pot obține astfel estimări ale valorilor proprii și stărilor proprii ale $H$ examinând valorile proprii ale $\tilde{H}$. De exemplu, estimarea energiei stării fundamentale se obține luând cea mai mică valoare proprie $c$, iar starea fundamentală din vectorul propriu corespunzător $\vec{c}$. Coeficienții din $\vec{c}$ determină contribuția diferiților vectori care generează $\mathcal{K}_U$.

Figura prezintă o reprezentare circuitară a testului Hadamard modificat, o metodă folosită pentru a calcula suprapunerea dintre diferite stări cuantice. Pentru fiecare element de matrice $\tilde{H}_{i,j}$, se efectuează un test Hadamard între stările $\vert \psi_i \rangle$, $\vert \psi_j \rangle$. Aceasta este evidențiată în figură prin schema de culori pentru elementele de matrice și operațiile corespunzătoare $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$. Astfel, este necesar un set de teste Hadamard pentru toate combinațiile posibile de vectori ai spațiului Krylov pentru a calcula toate elementele de matrice ale Hamiltonianului proiectat $\tilde{H}$. Firul de sus din circuitul testului Hadamard este un Qubit ansilă care este măsurat fie în baza X, fie în baza Y; valoarea sa de așteptare determină valoarea suprapunerii dintre stări. Firul de jos reprezintă toți Qubiții Hamiltonianului sistemului. Operația $\text{Prep} \; \psi_i$ pregătește Qubitul sistemului în starea $\vert \psi_i \rangle$ controlată de starea Qubitului ansilă (similar pentru $\text{Prep} \; \psi_j$), iar operația $P$ reprezintă descompunerea Pauli a Hamiltonianului sistemului $H = \sum_i P_i$. O derivare mai detaliată a operațiilor calculate de testul Hadamard este prezentată mai jos.

Definirea Hamiltonianului

Să considerăm Hamiltonianul Heisenberg pentru $N$ Qubiți pe un lanț liniar: $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)

[('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)]

Stabilirea parametrilor algoritmului

Alegem euristic o valoare pentru pasul de timp dt (bazată pe limite superioare ale normei Hamiltonianului). Ref [2] a arătat că un pas de timp suficient de mic este $\pi/\vert \vert H \vert \vert$, și că este preferabil, până la un punct, să subestimăm această valoare mai degrabă decât să o supraestimăm, deoarece supraestimarea poate permite contribuțiilor din stările de energie înaltă să corupă chiar și starea optimă din spațiul Krylov. Pe de altă parte, alegerea unui $dt$ prea mic duce la un condiționament mai slab al subspațiului Krylov, deoarece vectorii bazei Krylov diferă mai puțin de la un pas de timp la altul.

# 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

np.float64(0.10833078115826875)

Și setăm ceilalți parametri ai algoritmului. De dragul acestui tutorial, ne vom limita la utilizarea unui spațiu Krylov cu doar cinci dimensiuni, ceea ce este destul de restrictiv.

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

Pregătirea stării

Alege o stare de referință $\vert \psi \rangle$ care are o suprapunere cu starea fundamentală. Pentru acest Hamiltonian, folosim o stare cu o excitație în Qubitul din mijloc $\vert 00..010...00 \rangle$ ca stare de referință.

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

Evoluția temporală

Putem realiza operatorul de evoluție temporală generat de un Hamiltonian dat: $U=e^{-iHt}$ prin aproximarea 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)

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

Testul Hadamard

$$\begin{equation*} |0\rangle|0\rangle^N \quad\longrightarrow\quad \frac{1}{\sqrt{2}}\Big(|0\rangle + |1\rangle \Big)|0\rangle^N \quad\longrightarrow\quad \frac{1}{\sqrt{2}}\Big(|0\rangle|0\rangle^N+|1\rangle |\psi_i\rangle\Big) \quad\longrightarrow\quad \frac{1}{\sqrt{2}}\Big(|0\rangle |0\rangle^N+|1\rangle P |\psi_i\rangle\Big) \quad\longrightarrow\quad\frac{1}{\sqrt{2}}\Big(|0\rangle |\psi_j\rangle+|1\rangle P|\psi_i\rangle\Big) \end{equation*}$$

Unde $P$ este unul dintre termenii din descompunerea Hamiltonianului $H=\sum P$, iar $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$ sunt operații controlate care pregătesc vectorii $|\psi_i\rangle$, $|\psi_j\rangle$ ai spațiului Krylov unitar, cu $|\psi_k\rangle = e^{-i H k dt } \vert \psi \rangle = e^{-i H k dt } U_{\psi} \vert 0 \rangle^N$. Pentru a măsura $X$, se aplică mai întâi $H$...

$$\begin{equation*} \longrightarrow\quad\frac{1}{2}|0\rangle\Big( |\psi_j\rangle + P|\psi_i\rangle\Big) + \frac{1}{2}|1\rangle\Big(|\psi_j\rangle - P|\psi_i\rangle\Big) \end{equation*}$$

... apoi se măsoară:

$$\begin{equation*} \begin{split} \Rightarrow\quad\langle X\rangle &= \frac{1}{4}\Bigg(\Big\|| \psi_j\rangle + P|\psi_i\rangle \Big\|^2-\Big\||\psi_j\rangle - P|\psi_i\rangle\Big\|^2\Bigg) \\ &= \text{Re}\Big[\langle\psi_j| P|\psi_i\rangle\Big]. \end{split} \end{equation*}$$

Din identitatea $|a + b\|^2 = \langle a + b | a + b \rangle = \|a\|^2 + \|b\|^2 + 2\text{Re}\langle a | b \rangle$. Similar, măsurarea $Y$ produce

$$\begin{equation*} \langle Y\rangle = \text{Im}\Big[\langle\psi_j| P|\psi_i\rangle\Big]. \end{equation*}$$

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

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

Circuitul testului Hadamard poate fi un circuit adânc odată ce îl descompunem în porți native (care va crește și mai mult dacă luăm în considerare topologia dispozitivului)

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

Number of layers of 2Q operations 112753

Pasul 2: Optimizează problema pentru execuția pe hardware cuantic

Test Hadamard eficient

Putem optimiza circuitele adânci pentru testul Hadamard pe care le-am obținut prin introducerea unor aproximații și prin utilizarea unor ipoteze despre Hamiltonianul modelului. De exemplu, consideră următorul circuit pentru testul Hadamard:

Să presupunem că putem calcula clasic $E_0$, valoarea proprie a lui $|0\rangle^N$ sub Hamiltonianul $H$. Această condiție este satisfăcută atunci când Hamiltonianul păstrează simetria U(1). Deși poate părea o ipoteză puternică, există multe cazuri în care este sigur să presupunem că există o stare de vid (în acest caz aceasta corespunde stării $|0\rangle^N$) care nu este afectată de acțiunea Hamiltonianului. Acest lucru este adevărat, de exemplu, pentru Hamiltonieni de chimie care descriu molecule stabile (unde numărul de electroni este conservat). Dat fiind că poarta $\text{Prep} \; \psi$ pregătește starea de referință dorită $\ket{psi} = \text{Prep} \; \psi \ket{0} = e^{-i H 0 dt} U_{\psi} \ket{0}$, de exemplu, pentru a pregăti starea HF în chimie, $\text{Prep} \; \psi$ ar fi un produs de operații NOT pe qubiți individuali, deci $\text{Prep} \; \psi$ controlată este doar un produs de porți CNOT. Astfel, circuitul de mai sus implementează următoarea stare înainte de măsurare:

$$\begin{equation} \begin{split} \ket{0} \ket{0}^N\xrightarrow{H}&\frac{1}{\sqrt{2}} \left( \ket{0}\ket{0}^N+ \ket{1} \ket{0}^N \right)\\ \xrightarrow{\text{1-ctrl-init}}&\frac{1}{\sqrt{2}}\left(|0\rangle|0\rangle^N+|1\rangle|\psi\rangle\right)\\ \xrightarrow{U}&\frac{1}{\sqrt{2}}\left(e^{i\phi}\ket{0}\ket{0}^N+\ket{1} U\ket{\psi}\right)\\ \xrightarrow{\text{0-ctrl-init}}&\frac{1}{\sqrt{2}} \left( e^{i\phi}\ket{0} \ket{\psi} +\ket{1} U\ket{\psi} \right)\\ =&\frac{1}{2} \left( \ket{+}\left(e^{i\phi}\ket{\psi}+U\ket{\psi}\right) +\ket{-}\left(e^{i\phi}\ket{\psi}-U\ket{\psi}\right) \right)\\ =&\frac{1}{2} \left( \ket{+i}\left(e^{i\phi}\ket{\psi}-iU\ket{\psi}\right) +\ket{-i}\left(e^{i\phi}\ket{\psi}+iU\ket{\psi}\right) \right) \end{split} \end{equation}$$

unde am folosit deplasarea de fază calculabilă clasic $U\ket{0}^N = e^{i\phi}\ket{0}^N$ pe a treia linie. Prin urmare, valorile de așteptare se obțin ca

$$\begin{equation} \begin{split} \langle X\otimes P\rangle&=\frac{1}{4} \Big( \left(e^{-i\phi}\bra{\psi}+\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}+U\ket{\psi}\right) \\ &\qquad-\left(e^{-i\phi}\bra{\psi}-\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}-U\ket{\psi}\right) \Big)\\ &=\text{Re}\left[e^{-i\phi}\bra{\psi}PU\ket{\psi}\right], \end{split} \end{equation}$$ $$\begin{equation} \begin{split} \langle Y\otimes P\rangle&=\frac{1}{4} \Big( \left(e^{-i\phi}\bra{\psi}+i\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}-iU\ket{\psi}\right) \\ &\qquad-\left(e^{-i\phi}\bra{\psi}-i\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}+iU\ket{\psi}\right) \Big)\\ &=\text{Im}\left[e^{-i\phi}\bra{\psi}PU\ket{\psi}\right]. \end{split} \end{equation}$$

Folosind aceste ipoteze am reușit să scriem valorile de așteptare ale operatorilor de interes cu mai puține operații controlate. De fapt, trebuie să implementăm doar pregătirea controlată a stării $\text{Prep} \; \psi$, nu și evoluțiile temporale controlate. Reformulând calculul nostru ca mai sus, vom putea reduce considerabil adâncimea circuitelor rezultante.

Descompune operatorul de evoluție temporală cu descompunerea Trotter

În loc să implementăm operatorul de evoluție temporală exact, putem folosi descompunerea Trotter pentru a implementa o aproximație a acestuia. Repetând de mai multe ori o descompunere Trotter de un anumit ordin obținem o reducere suplimentară a erorii introduse de aproximație. În continuare, construim direct implementarea Trotter în cel mai eficient mod pentru graful de interacțiune al Hamiltonianului pe care îl considerăm (numai interacțiuni între vecini direcți). În practică, inserăm rotații Pauli $R_{xx}$, $R_{yy}$, $R_{zz}$ cu un unghi parametrizat $t$, care corespund implementării aproximative a lui $e^{-i (XX + YY + ZZ) t}$. Dată diferența dintre definiția rotațiilor Pauli și evoluția temporală pe care încercăm să o implementăm, va trebui să folosim parametrul $2*dt$ pentru a obține o evoluție temporală de $dt$. În plus, inversăm ordinea operațiilor pentru un număr impar de repetări ale pașilor Trotter, ceea ce este funcțional echivalent, dar permite sinteza operațiilor adiacente într-un singur unitar $SU(2)$. Aceasta produce un circuit mult mai puțin adânc decât cel obținut cu funcționalitatea generică 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)

Folosește un circuit optimizat pentru pregătirea stării

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)

Circuite șablon pentru calculul elementelor de matrice ale $\tilde{S}$ și $\tilde{H}$ prin testul Hadamard

Singura diferență dintre circuitele utilizate în testul Hadamard va fi faza din operatorul de evoluție temporală și observabilele măsurate. Prin urmare, putem pregăti un circuit șablon care să reprezinte circuitul generic pentru testul Hadamard, cu marcatori pentru porțile care depind de operatorul de evoluție temporală.

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

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

The optimized circuit has 2Q gates depth:  74

Am redus considerabil adâncimea testului Hadamard printr-o combinație de aproximație Trotter și unitare necontrolate

Pasul 3: Execuție folosind primitivele Qiskit

Instanțiază Backend-ul și setează parametrii de execuție

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)

Transpilare către un QPU

Mai întâi, să selectăm subseturi ale hărții de cuplare cu qubiți cu performanțe „bune" (unde „bun" este destul de arbitrar — dorim în principal să evităm qubiți cu performanțe foarte slabe) și să creăm o țintă nouă pentru transpilare

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

Apoi transpilăm circuitul virtual la cel mai bun layout fizic din această nouă țintă

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

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]

Creează PUB-uri pentru execuție cu 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)

Rularea circuitelor

Circuitele pentru $t=0$ sunt calculabile clasic

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)

(25+0j)

Execută circuitele pentru $S$ și $\tilde{H}$ cu 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])

Pasul 4: Post-procesare și returnarea rezultatului în formatul clasic dorit

results = job.result()[0]

Calcularea matricelor Hamiltonianului efectiv și de suprapunere

Mai întâi calculăm faza acumulată de starea $\vert 0 \rangle$ în timpul evoluției temporale necontrolate

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

Odată ce avem rezultatele execuțiilor circuitelor, putem post-procesa datele pentru a calcula elementele matricei $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)

$$\displaystyle \left[\begin{matrix}1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i & -0.180546747798251 - 0.492624093654174 i & 0.0012070853532697 + 0.312052218182462 i\\-0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i & -0.180546747798251 - 0.492624093654174 i\\0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i\\-0.180546747798251 + 0.492624093654174 i & 0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i\\0.0012070853532697 - 0.312052218182462 i & -0.180546747798251 + 0.492624093654174 i & 0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0\end{matrix}\right]$$

Și elementele matricei $\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)

$$\displaystyle \left[\begin{matrix}25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i & -5.15587257589417 - 8.88280836036843 i & 1.98818301405581 + 5.8897614762563 i\\-14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i & -5.15587257589417 - 8.88280836036843 i\\10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i\\-5.15587257589417 + 8.88280836036843 i & 10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i\\1.98818301405581 - 5.8897614762563 i & -5.15587257589417 + 8.88280836036843 i & 10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0\end{matrix}\right]$$

În final, putem rezolva problema valorilor proprii generalizate pentru $\tilde{H}$:

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

și obținem o estimare a energiei stării fundamentale $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)

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

Pentru un sector cu o singură particulă, putem calcula eficient clasic starea fundamentală a acestui sector al Hamiltonianului

gs_en = single_particle_gs(H_op, n_qubits)

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

Anexă: Subspațiul Krylov din evoluții în timp real

Spațiul Krylov unitar este definit ca

$$\mathcal{K}_U(H, |\psi\rangle) = \text{span}\left\{ |\psi\rangle, e^{-iH\,dt} |\psi\rangle, \dots, e^{-irH\,dt} |\psi\rangle \right\}$$

pentru un pas de timp $dt$ pe care îl vom determina ulterior. Presupunem temporar că $r$ este par: definim $d=r/2$. Observăm că atunci când proiectăm Hamiltonianul în spațiul Krylov de mai sus, acesta este indistinguibil de spațiul Krylov

$$\mathcal{K}_U(H, |\psi\rangle) = \text{span}\left\{ e^{i\,d\,H\,dt}|\psi\rangle, e^{i(d-1)H\,dt} |\psi\rangle, \dots, e^{-i(d-1)H\,dt} |\psi\rangle, e^{-i\,d\,H\,dt} |\psi\rangle \right\},$$

adică acolo unde toate evoluțiile temporale sunt deplasate înapoi cu $d$ pași de timp. Motivul pentru care este indistinguibil constă în faptul că elementele de matrice

$$\tilde{H}_{j,k} = \langle\psi|e^{i\,j\,H\,dt}He^{-i\,k\,H\,dt}|\psi\rangle=\langle\psi|He^{i(j-k)H\,dt}|\psi\rangle$$

sunt invariante la deplasări globale ale timpului de evoluție, deoarece evoluțiile temporale comutează cu Hamiltonianul. Pentru $r$ impar, putem folosi analiza pentru $r-1$.

Dorim să arătăm că undeva în acest spațiu Krylov există garantat o stare de energie scăzută. Facem asta prin intermediul următorului rezultat, derivat din Teorema 3.1 din [3]:

Afirmația 1: există o funcție $f$ astfel încât pentru energii $E$ în intervalul spectral al Hamiltonianului (adică între energia stării fundamentale și energia maximă)...

1. $f(E_0)=1$
2. $|f(E)|\le2\left(1 + \delta\right)^{-d}$ pentru toate valorile $E$ aflate la distanță $\ge\delta$ față de $E_0$, adică este suprimată exponențial
3. $f(E)$ este o combinație liniară de $e^{ijE\,dt}$ pentru $j=-d,-d+1,...,d-1,d$

Oferim mai jos o demonstrație, care poate fi sărită în siguranță dacă nu dorim să înțelegem argumentul complet și riguros. Acum ne concentrăm pe implicațiile afirmației de mai sus. Prin proprietatea 3 de mai sus, putem vedea că spațiul Krylov deplasat de mai sus conține starea $f(H)|\psi\rangle$. Aceasta este starea noastră de energie scăzută. Pentru a înțelege de ce, scriem $|\psi\rangle$ în baza proprie de energie:

$$|\psi\rangle = \sum_{k=0}^{N}\gamma_k|E_k\rangle,$$

unde $|E_k\rangle$ este cea de-a k-a stare proprie de energie și $\gamma_k$ este amplitudinea sa în starea inițială $|\psi\rangle$. Exprimat în acești termeni, $f(H)|\psi\rangle$ este dat de

$$f(H)|\psi\rangle = \sum_{k=0}^{N}\gamma_kf(E_k)|E_k\rangle,$$

folosind faptul că putem înlocui $H$ cu $E_k$ atunci când acționează asupra stării proprii $|E_k\rangle$. Eroarea de energie a acestei stări este, prin urmare,

$$\text{energy error} = \frac{\langle\psi|f(H)(H-E_0)f(H)|\psi\rangle}{\langle\psi|f(H)^2|\psi\rangle}$$ $$= \frac{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2}.$$

Pentru a transforma aceasta într-o margine superioară mai ușor de înțeles, separăm mai întâi suma din numărător în termeni cu $E_k-E_0\le\delta$ și termeni cu $E_k-E_0>\delta$:

$$\text{energy error} = \frac{\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} + \frac{\sum_{E_k> E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2}.$$

Putem margini superior primul termen prin $\delta$,

$$\frac{\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} < \frac{\delta\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} \le \delta,$$

unde primul pas urmează deoarece $E_k-E_0\le\delta$ pentru fiecare $E_k$ din sumă, iar al doilea pas urmează deoarece suma din numărător este un subset al sumei din numitor. Pentru al doilea termen, mărginim inferior numitorul prin $|\gamma_0|^2$, deoarece $f(E_0)^2=1$: adunând totul înapoi, obținem

$$\text{energy error} \le \delta + \frac{1}{|\gamma_0|^2}\sum_{E_k>E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0).$$

Pentru a simplifica restul, observăm că pentru toți acești $E_k$, prin definiția lui $f$ știm că $f(E_k)^2 \le 4\left(1 + \delta\right)^{-2d}$. Mărginind suplimentar $E_k-E_0<2\|H\|$ și mărginind $\sum_{E_k>E_0+\delta}|\gamma_k|^2<1$ obținem

$$\text{energy error} \le \delta + \frac{8}{|\gamma_0|^2}\|H\|\left(1 + \delta\right)^{-2d}.$$

Aceasta este valabilă pentru orice $\delta>0$, deci dacă setăm $\delta$ egal cu eroarea noastră țintă, atunci marginea de eroare de mai sus converge exponențial spre aceasta cu dimensiunea Krylov $2d=r$. De asemenea, rețineți că dacă $\delta<E_1-E_0$ atunci termenul $\delta$ dispare complet în marginea de mai sus.

Pentru a completa argumentul, menționăm că cele de mai sus reprezintă doar eroarea de energie a stării particulare $f(H)|\psi\rangle$, nu eroarea de energie a stării cu energia cea mai mică din spațiul Krylov. Cu toate acestea, prin principiul variațional (Rayleigh-Ritz), eroarea de energie a stării cu energia cea mai mică din spațiul Krylov este mărginită superior de eroarea de energie a oricărei stări din spațiul Krylov, deci cele de mai sus reprezintă și o margine superioară pentru eroarea de energie a stării cu energia cea mai mică, adică rezultatul algoritmului de diagonalizare cuantică Krylov.

O analiză similară cu cea de mai sus poate fi efectuată ținând cont suplimentar de zgomot și de procedura de prag discutată în notebook. Consultați [2] și [4] pentru această analiză.

Anexă: demonstrația Afirmației 1

Următoarele sunt derivate în mare parte din [3], Teorema 3.1: fie $0 < a < b$ și fie $\Pi^*_d$ spațiul polinoamelor reziduale (polinoame a căror valoare în 0 este 1) de grad cel mult $d$. Soluția la

$$\beta(a, b, d) = \min_{p \in \Pi^*_d} \max_{x \in [a, b]} |p(x)| \quad$$

este

$$p^*(x) = \frac{T_d\left(\frac{b + a - 2x}{b - a}\right)}{T_d\left(\frac{b + a}{b - a}\right)}, \quad$$

și valoarea minimă corespunzătoare este

$$\beta(a, b, d) = T_d^{-1}\left(\frac{b + a}{b - a}\right).$$

Dorim să convertim aceasta într-o funcție care poate fi exprimată natural în termeni de exponențiale complexe, deoarece acestea sunt evoluțiile în timp real care generează spațiul cuantic Krylov. Pentru a face aceasta, este convenabil să introducem următoarea transformare a energiilor din intervalul spectral al Hamiltonianului în numere din intervalul $[0,1]$: definim

$$g(E) = \frac{1-\cos\big((E-E_0)dt\big)}{2},$$

unde $dt$ este un pas de timp astfel încât $-\pi < E_0dt < E_\text{max}dt < \pi$. Observăm că $g(E_0)=0$ și $g(E)$ crește pe măsură ce $E$ se îndepărtează de $E_0$.

Acum folosind polinomul $p^*(x)$ cu parametrii a, b, d setați la $a = g(E_0 + \delta)$, $b = 1$, și d = int(r/2), definim funcția:

$$f(E) = p^* \left( g(E) \right) = \frac{T_d\left(1 + 2\frac{\cos\big((E-E_0)dt\big) - \cos\big(\delta\,dt\big)}{1 +\cos\big(\delta\,dt\big)}\right)}{T_d\left(1 + 2\frac{1-\cos\big(\delta\,dt\big)}{1 + \cos\big(\delta\,dt\big)}\right)}$$

unde $E_0$ este energia stării fundamentale. Putem vedea prin inserarea $\cos(x)=\frac{e^{ix}+e^{-ix}}{2}$ că $f(E)$ este un polinom trigonometric de grad $d$, adică o combinație liniară de $e^{ijE\,dt}$ pentru $j=-d,-d+1,...,d-1,d$. Mai mult, din definiția lui $p^*(x)$ de mai sus avem că $f(E_0)=p(0)=1$ și pentru orice $E$ din intervalul spectral astfel încât $\vert E-E_0 \vert > \delta$ avem

$$|f(E)| \le \beta(a, b, d) = T_d^{-1}\left(1 + 2\frac{1-\cos\big(\delta\,dt\big)}{1 + \cos\big(\delta\,dt\big)}\right)$$ $$\leq 2\left(1 + \delta\right)^{-d} = 2\left(1 + \delta\right)^{-\lfloor k/2\rfloor}.$$

Referințe

[1] N. Yoshioka, M. Amico, W. Kirby et al. "Diagonalization of large many-body Hamiltonians on a quantum processor". arXiv:2407.14431

[2] Ethan N. Epperly, Lin Lin, and Yuji Nakatsukasa. "A theory of quantum subspace diagonalization". SIAM Journal on Matrix Analysis and Applications 43, 1263–1290 (2022).

[3] Å. Björck. "Numerical methods in matrix computations". Texts in Applied Mathematics. Springer International Publishing. (2014).

[4] William Kirby. "Analysis of quantum Krylov algorithms with errors". Quantum 8, 1457 (2024).

Sondaj tutorial

Te rog să completezi acest scurt sondaj pentru a oferi feedback despre acest tutorial. Părerile tale ne vor ajuta să îmbunătățim conținutul și experiența utilizatorilor.

Link to survey

Note: This survey is provided by IBM Quantum and relates to the original English content. To give feedback on doQumentation's website, translations, or code execution, please open a GitHub issue.