Krylov-Quantendiagonaliseren vun Gitterhamiltonianen

Bruukschattung: 20 Minuten op'n Heron r2 (HENWIES: Dit is blot'n Schattung. Dien Looptied kann anners sien.)

# 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

Achtergrund

Dit Tutorial wiest, wo man den Krylov Quantum Diagonalization Algorithm (KQD) binnen dat Qiskit-Muster ümsetten kann. Du lehrst toeerst över de Theorie achter den Algorithmus un denn seehst 'ne Vörföhrung vun sien Utföhren op'n QPU.

Över all Disziplinen sünd wi daran interesseert, Grundtostandseigenschappen vun Quantensystemen to lehren. Biespelen ümfaten dat Verstahn vun de fundamentale Natur vun Partikeln un Kräften, dat Vörütseggen un Verstahn vun dat Verholten vun komplexe Materialien un dat Verstahn vun biochemische Interaktionen un Reaktionen. Wegen dat exponentielle Wassen vun den Hilbert-Ruum un de Korrelation, de in verschränkte Systemen optreden, hebbt klassische Algorithmen Swierigkeden, dit Problem för Quantensystemen vun wassende Grött to lösen. An een Enn vun dat Spektrum is de bestahnde Angang, de vun de Quantenhardware profiteert un sik op variationelle Quantenmethoden konzentriert (to'n Bispill variational quantum eigensolver). Disse Techniken stahn vör Rutforderns mit aktuelle Geräten wegen de hoge Antall vun Funktioneopropen, de in'n Optimerungsprozess nödig sünd, wat 'n groten Overhead an Ressourcen bringt, wenn moderne Fehlerminnerungstechniken infohrt warrt, un dormit ehr Warksamkeit op lütte Systemen begrenzt. An dat anner Enn vun dat Spektrum gifft dat fehlertoleranteQuantenmethoden mit Leistungsgarantien (to'n Bispill quantum phase estimation), de deep Schaltkreisen bruken, de blot op'n fehlertoleranten Gerät utfohrt warrn köönt. Ut disse Grünn stellt wi hier'n Quantenalgorithmus vör, de op Deelruummethoden baseert (as beschreven in dit Översichtsppapier), den Krylov-Quantendiagonaliserungsalgorithmus (KQD). Disse Algorithmus warkt good bi groten Matstav [1] op bestahnde Quantenhardware, deelt ähnliche Leistungsgarantien as Phase Estimation, is kompatibel mit moderne Fehlerminnerungstechniken un künn Resultaten levern, de klassisch nich togänglich sünd.

Vöruttösettungen

Bevör du mit dit Tutorial anfängst, sörg dafür, dat du Folgendes installeert hest:

• Qiskit SDK v2.0 oder neger, mit Visualiseren
• Qiskit Runtime v0.22 oder neger ( pip install qiskit-ibm-runtime )

Opbau

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

Schritt 1: Kart klassische Ingänge op'n Quantenproblem

De Krylov-Ruum

De Krylov-Ruum $\mathcal{K}^r$ vun Ordnung $r$ is de Ruum, de vun Vektoren opspannt warrt, de wi kriegt, indem wi högere Potenzen vun'ne Matrix $A$, bet to $r-1$, mit'n Referenzvekter $\vert v \rangle$ multiplizeren.

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

Wenn de Matrix $A$ de Hamiltonian $H$ is, warrt wi den tohörigen Ruum as Potenz-Krylov-Ruum $\mathcal{K}_P$ bezeeken. In't Fall, wenn $A$ de Tiedentwikkelungsoperator is, de vun den Hamiltonian $U=e^{-iHt}$ erzeugt warrt, warrt wi den Ruum as unitären Krylov-Ruum $\mathcal{K}_U$ bezeeken. De Potenz-Krylov-Deelruum, den wi klassisch bruken, kann nich direkt op'n Quantencomputer erzeugt warrn, wieldat $H$ keen unitären Operator is. In steed darvun köön wi den Tiedentwikkelungsoperator $U = e^{-iHt}$ bruken, vun den man wiesen kann, dat he ähnliche Konvergenzgarantien as de Potenzmethode gifft. Potenzen vun $U$ warrt denn verscheden Tiedschritte $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\}$$

Seeh in'n Anhang för'ne detailleerte Herlieden, wo de unitäre Krylov-Ruum dat tolaten deit, neederenergetische Eigentostänn genau dartostelln.

Krylov-Quantendiagonaliserungsalgorithmus

Geven'n Hamiltonian $H$, den wi diagonaliseren wüllt, betrachten wi toeerst den tohörigen unitären Krylov-Ruum $\mathcal{K}_U$. Dat Teel is dat, 'ne kompakte Darstellen vun den Hamiltonian in $\mathcal{K}_U$ to finnen, de wi as $\tilde{H}$ bezeeken warrt. De Matrixelemente vun $\tilde{H}$, de Projektion vun den Hamiltonian in'n Krylov-Ruum, köön berekent warrn, indem wi de folgenden Verwachtungsweerten berekenen

$$\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$$

Wo $\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle$ de Vektoren vun den unitären Krylov-Ruum sünd un $t_n = n dt$ de Veelfachen vun den Tiedschritt $dt$ sünd, de wählt wurrn. Op'n Quantencomputer kann de Bereknung vun jedes Matrixelement mit jeden Algorithmus makt warrn, de dat tolaten deit, Överlappungen twüschen Quantentostänn to kregen. Dit Tutorial konzentreert sik op den Hadamard-Test. Geven dat $\mathcal{K}_U$ de Dimension $r$ hett, hett de Hamiltonian, projizeert in den Deelruum, Dimensionen $r \times r$. Mit $r$ lütt noog (normalerwies langt $r<<100$ för de Konvergenz vun Schattungen vun Eigenenergien) köön wi denn eenfach den projizierten Hamiltonian $\tilde{H}$ diagonaliseren. Aber wi köön $\tilde{H}$ nich direkt diagonaliseren wegen de Nich-Orthogonalität vun de Krylov-Ruum-Vektoren. Wi mööt ehr Överlappungen meten un 'ne Matrix $\tilde{S}$ konstrueren

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

Dat erlööft uns, dat Eigenweertproblem in'n nich-orthogonalen Ruum to lösen (ook verallgemeinert Eigenweertproblem neemt)

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

Man kann denn Schattungen vun de Eigenweerten un Eigentostänn vun $H$ kregen, indem man op de vun $\tilde{H}$ kiekt. To'n Bispill kriggt man de Schatting vun de Grundtostandsenergie, indem man den lüttstesten Eigenweert $c$ nimmt un den Grundtostand ut den tohörigen Eigenvekter $\vec{c}$. De Koeffizienten in $\vec{c}$ bestimmen den Bidrag vun de verschedenen Vektoren, de $\mathcal{K}_U$ opspannen.

De Figur wiest 'ne Schaltkreisdarstellen vun den modifizierten Hadamard-Test, 'ne Methode, de bruukt warrt, üm de Överlappung twüschen verschedene Quantentostänn to berekenen. För jedes Matrixelement $\tilde{H}_{i,j}$ warrt'n Hadamard-Test twüschen den Tostand $\vert \psi_i \rangle$, $\vert \psi_j \rangle$ dörchmakt. Dat warrt in de Figur dörch dat Farvenschema för de Matrixelemente un de tohörigen $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$ Operationen hervörhaven. Also is'ne Menge vun Hadamard-Tests för all möglichen Kombinationen vun Krylov-Ruum-Vektoren nödig, üm all Matrixelemente vun den projizierten Hamiltonian $\tilde{H}$ to berekenen. De bovenste Lien in'n Hadamard-Test-Schaltkreis is'n Ancilla-Qubit, dat entweder in de X- oder Y-Basis meten warrt, sien Verwachtungsweert bestimmt den Weert vun de Överlappung twüschen de Tostänn. De unnerste Lien stellt all Qubits vun den System-Hamiltonian dar. De $\text{Prep} \; \psi_i$ Operation bereit dat System-Qubit in'n Tostand $\vert \psi_i \rangle$ vör, stüürt dörch den Tostand vun dat Ancilla-Qubit (ähnlich för $\text{Prep} \; \psi_j$) un de Operation $P$ stellt de Pauli-Zerleggen vun den System-Hamiltonian $H = \sum_i P_i$ dar. 'Ne detailleerte Herlieden vun de Operationen, de vun den Hadamard-Test berekent warrt, is ünnen geven.

Hamiltonian defineren

Lot uns den Heisenberg-Hamiltonian för $N$ Qubits op'ne lineere Kette betrachten: $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)]

Parameter för den Algorithmus fastleggen

Wi wählen heuristisch'n Weert för den Tiedschritt dt (baseert op bovenste Grenzen vun de Hamiltonian-Norm). Ref [2] hett wiest, dat'n genoog lütten Tiedschritt $\pi/\vert \vert H \vert \vert$ is, un dat dat bet to'n Punkt beter is, dissen Weert to ünnerschatten in steed vun överschatten, wieldat Överschatten Bidräg vun hoogenenergetische Tostänn tolatenlaten kann, de sölfst den optimalen Tostand in'n Krylov-Ruum verdarven köön. Annerersied föhrt de Wahl vun $dt$ to lütt darto, dat de Konditioneren vun den Krylov-Deelruum slechter warrt, wieldat de Krylov-Basisvektoren sik weniger vun Tiedschritt to Tiedschritt ünnerscheden.

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

Un sett anner Parameter vun den Algorithmus. För dat Tutorial beschränken wi uns op de Nutzung vun'n Krylov-Ruum mit blot fief Dimensionen, wat temlich beschränkend is.

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

Tostandsvörbereeden

Wähl'n Referenztostand $\vert \psi \rangle$, de irgendwie Överlappung mit den Grundtostand hett. För dissen Hamiltonian bruken wi den Tostand mit'ne Anregung in't middelste Qubit $\vert 00..010...00 \rangle$ as unsen Referenztostand.

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

Tiedentwikkelung

Wi köön den Tiedentwikkelungsoperator, de vun'n geven Hamiltonian erzeugt warrt, realiseren: $U=e^{-iHt}$ över de Lie-Trotter-Approximation.

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>

Hadamard-Test

$$\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*}$$

Wo $P$ een vun de Termen in de Zerleggen vun den Hamiltonian $H=\sum P$ is un $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$ stüürte Operationen sünd, de $|\psi_i\rangle$, $|\psi_j\rangle$ Vektoren vun den unitären Krylov-Ruum vörbereeden, mit $|\psi_k\rangle = e^{-i H k dt } \vert \psi \rangle = e^{-i H k dt } U_{\psi} \vert 0 \rangle^N$. Üm $X$ to meten, wendt toeerst $H$ an...

$$\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*}$$

... denn meet:

$$\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*}$$

Ut de Identität $|a + b\|^2 = \langle a + b | a + b \rangle = \|a\|^2 + \|b\|^2 + 2\text{Re}\langle a | b \rangle$. Ähnlich levvt de Metung vun $Y$

$$\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

De Hadamard-Test-Schaltkreis kann'n depen Schaltkreis sien, sobold wi in native Gates zerleggen (wat noch mehr tonimmt, wenn wi de Topologie vun dat Gerät berücksichtigen)

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

Schritt 2: Optimeer dat Problem för de Utföhren op Quantenhardware

Effizienten Hadamard-Test

Wi köön de depen Schaltkreisen för den Hadamard-Test optimeren, de wi kregen hebbt, indem wi'n poor Approximationen infohren un uns op'n poor Annahmen över den Modell-Hamiltonian verlaten. To'n Bispill betrach den folgenden Schaltkreis för den Hadamard-Test:

Nimm an, wi köön klassisch $E_0$ berekenen, den Eigenweert vun $|0\rangle^N$ ünner den Hamiltonian $H$. Dat is erfüllt, wenn de Hamiltonian de U(1)-Symmetrie erhollt. Obschoonst dat 'ne starke Annahm schinkt, gifft dat vele Fäll, wo dat seker is antonahmen, dat dat'n Vakuumtostand gifft (in dissen Fall kümmt dat to den $|0\rangle^N$ Tostand), de vun de Warken vun den Hamiltonian unberührt blifft. Dat stimmt to'n Bispill för Chemie-Hamiltoniane, de stabiele Moleküle beschrieven (wo de Antall vun Elektronen erholden blifft). Geven dat de Poort $\text{Prep} \; \psi$ den wünschten Referenztostand vörbereidt $\ket{psi} = \text{Prep} \; \psi \ket{0} = e^{-i H 0 dt} U_{\psi} \ket{0}$, to'n Bispill üm den HF-Tostand för Chemie vörtobereeden, wör $\text{Prep} \; \psi$'n Produkt vun Enkel-Qubit-NOTs sien, also is stüürt-$\text{Prep} \; \psi$ blot'n Produkt vun CNOTs. Denn implementeert de Schaltkreis baven den folgenden Tostand vör de Metung:

$$\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}$$

wo wi de klassisch simuleerbare Phasenverschuven $U\ket{0}^N = e^{i\phi}\ket{0}^N$ in de drütte Reeg bruukt hebbt. Dorher warrt de Verwachtungsweerten kregen as

$$\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}$$

Mit disse Annahmen kunnen wi de Verwachtungsweerten vun Operatoren vun Interesse mit wenijer stüürte Operationen schrieven. Tatsächlich mööt wi blot de stüürte Tostandsvörbereeden $\text{Prep} \; \psi$ implementeren un nich de stüürten Tiedentwikkelungen. Unse Bereknung as baven neü to formuleren erlööft uns, de Deep vun de resultierenden Schaltkreisen stark to reduzeren.

Zerlegg den Tiedentwikkelungsoperator mit de Trotter-Zerleggen

Anstatt den Tiedentwikkelungsoperator genau to implementeren, köön wi de Trotter-Zerleggen bruken, üm 'ne Approximation darvun to implementeren. Wenn wi 'ne bestimmte Ordnung vun de Trotter-Zerleggen mehrmals weerholen, kriegen wi 'ne wiederere Reduktion vun den Fehler, de dörch de Approximation infohrt warrt. In't Folgende boen wi de Trotter-Implementierung direkt op de effizienteste Wies för den Interaktionsgraphen vun den Hamiltonian, den wi betrachten (blot Nächst-Nachbar-Interaktionen). In de Praxis setten wi Pauli-Rotationen $R_{xx}$, $R_{yy}$, $R_{zz}$ mit'n parametrisierten Winkel $t$ in, de de approximierte Implementierung vun $e^{-i (XX + YY + ZZ) t}$ entspreken. Geven den Ünnerscheed in de Definition vun de Pauli-Rotationen un de Tiedentwikkelung, de wi implementeren wüllt, mööt wi den Parameter $2*dt$ bruken, üm 'ne Tiedentwikkelung vun $dt$ to erreken. Buten dat kehren wi de Regenfolg vun de Operationen för ungerade Antallen vun Trotter-Schritten üm, wat funktionell äquivalent is, aver dat Synthetiseren vun benachbarte Operationen in 'ne enkel $SU(2)$ Unitäre erlööft. Dat gifft 'n veel flaacheren Schaltkreis as dat, wat man mit de generische PauliEvolutionGate()-Funktionalität kriegt.

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)

Bruuk 'n optimierten Schaltkreis för de Tostandsvörbereeden

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)

Schablonenschaltkreisen för de Bereknung vun Matrixelemente vun $\tilde{S}$ un $\tilde{H}$ över den Hadamard-Test

De eenzige Ünnerscheed twüschen de Schaltkreisen, de in'n Hadamard-Test bruukt warrt, is de Phase in'n Tiedentwikkelungsoperator un de Observablen, de meten warrt. Dorher köön wi 'n Schablonenschaltkreis vörbereeden, de den generischen Schaltkreis för den Hadamard-Test darstellt, mit Platzholders för de Gates, de vun den Tiedentwikkelungsoperator afhängt.

# 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

Wi hebbt de Deep vun den Hadamard-Test bedüüdend reduziert mit 'ne Kombination vun Trotter-Approximation un unstüürte Unitäre.

Schritt 3: Utföhren mit Qiskit Primitives

Instanzieer dat Backend un sett Runtime-Parameter

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)

Transpilierung to 'ne QPU

Eerst wählen wi Deelmengen vun de Koppelungskart mit "gode" Qubits ut (wobei "good" hier temlich willkürlich is, wi wüllt vör allem echt slecht performende Qubits vermieden) un maakt 'n nee't Target för de Transpilierung.

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

Denn transpileer den virtuellen Schaltkreis to dat beste physische Layout in dit nee't Target.

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

qc_trans = pm.run(qc)

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

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]

PUBs för de Utföhrung mit den Estimator erstellen

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

Schaltkreisen utföhren

Schaltkreisen för $t=0$ sünd klassisch berekenbor.

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)

Föhr de Schaltkreisen för $S$ un $\tilde{H}$ mit den Estimator ut.

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

Schritt 4: Nabearbeiten un Törüchgeven vun't Resultat in't wünschte klassische Format

results = job.result()[0]

Effektiven Hamiltonian un Överlappungsmatrizen bereken

Eerst bereken de Phase, de vun den $\vert 0 \rangle$ Tostand wäährend de unstüürte Tiedentwikkelung opsammelt warrt.

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

Sobald wi de Resultaten vun de Schaltkreisutföhrungen hebbt, köön wi de Daten nabearbeiten, üm de Matrixelemente vun $S$ to bereken.

# 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]$$

Un de Matrixelemente vun $\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]$$

Toletzt köön wi dat verallgemeineert Eigenweertproblem för $\tilde{H}$ lösen:

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

un 'ne Schatting vun de Grundtostandsenergie $c_{min}$ kregen.

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

För 'n Eendeelchen-Sektor köön wi den Grundtostand vun dissen Sektor vun den Hamiltonian effizient klassisch bereken.

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

Anhang: Krylov-Deelruum ut reelle Tiedentwikkelungen

De unitäre Krylov-Ruum is defineert as

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

för 'n Tiedschritt $dt$, den wi later bestimmen warrt. Nimm vöreerst an, dat $r$ gerade is: denn defineer $d=r/2$. Acht dorop, dat wenn wi den Hamiltonian in den Krylov-Ruum baven projiziern, dat nich to ünnerscheden is vun den Krylov-Ruum

$$\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\},$$

dat hett, wo all de Tiedentwikkelungen üm $d$ Tiedschritte na achtern verschaven sünd. De Grund, worüm dat nich to ünnerscheden is, is, dat de Matrixelemente

$$\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$$

invariant ünner Gesamtverschiebungen vun de Entwikkelungstied sünd, wieldat de Tiedentwikkelungen mit den Hamiltonian kommuteren. För ungerade $r$ köön wi de Analyse för $r-1$ bruken.

Wi wüllt wiesen, dat irgendwo in dissen Krylov-Ruum garantiert 'n neederenergetischer Tostand is. Dat doon wi mit Hülp vun dat folgende Resultat, wat ut Theorem 3.1 in [3] afleidt is:

Behauptung 1: Dat gifft 'ne Funktion $f$, so dat för Energien $E$ in'n Spektralbereich vun den Hamiltonian (dat hett, twüschen de Grundtostandsenergie un de maximale Energie)...

1. $f(E_0)=1$
2. $|f(E)|\le2\left(1 + \delta\right)^{-d}$ för all Weerten vun $E$, de $\ge\delta$ weg vun $E_0$ liggt, dat hett, dat is exponentiell ünnerdröckt
3. $f(E)$ is 'ne Linearkombination vun $e^{ijE\,dt}$ för $j=-d,-d+1,...,d-1,d$

Wi geven ünnen 'n Bewies, aver den kann man gerust överspringen, wenn man nich dat vullständige, rigorose Argument verstahn will. För nu konzentreern wi uns op de Implikationen vun de baven nömmte Behauptung. Dörch Eigenschap 3 baven köön wi sehn, dat de verschavene Krylov-Ruum baven den Tostand $f(H)|\psi\rangle$ enthölt. Dat is unsen neederenergetische Tostand. Üm to sehn, worüm, schrieven wi $|\psi\rangle$ in de Energie-Eigenbasis:

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

wo $|E_k\rangle$ de k-te Energie-Eigentostand is un $\gamma_k$ sien Amplitude in'n Anfangstostand $|\psi\rangle$. Utdrückt in düsse Termen is $f(H)|\psi\rangle$ geven as

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

mit de Tatsach, dat wi $H$ dörch $E_k$ ersetten köön, wenn dat op den Eigentostand $|E_k\rangle$ warkt. De Energiefehler vun dissen Tostand is dorher

$$\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}.$$

Üm dat in 'ne bovenste Grenz ümtowandelln, de leechter to verstahn is, trennden wi eerst de Summ in'n Teller in Termen mit $E_k-E_0\le\delta$ un Termen mit $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}.$$

Wi köön den eersten Term na baven dörch $\delta$ begrenzen,

$$\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,$$

wo de eerste Schritt folgt, wieldat $E_k-E_0\le\delta$ för jeden $E_k$ in de Summ gellt, un de tweete Schritt folgt, wieldat de Summ in'n Teller 'ne Deelmenge vun de Summ in'n Nenner is. För den tweeten Term begrenzen wi eerst den Nenner na ünnen dörch $|\gamma_0|^2$, wieldat $f(E_0)^2=1$: wenn wi allens weerher tosamenföhrt, gifft dat

$$\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).$$

Üm dat, wat övrig blifft, to vereenfachen, acht dorop, dat för all disse $E_k$, na de Definition vun $f$, wi weten, dat $f(E_k)^2 \le 4\left(1 + \delta\right)^{-2d}$. Wenn wi buten dat $E_k-E_0<2\|H\|$ na baven begrenzen un $\sum_{E_k>E_0+\delta}|\gamma_k|^2<1$ na baven begrenzen, gifft dat

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

Dit gellt för jeden $\delta>0$, also wenn wi $\delta$ gliek unsen Teelfehler setten, denn konvergeert de Fehlergrenz baven exponentiell mit de Krylov-Dimension $2d=r$ dorhen. Acht ok dorop, dat wenn $\delta<E_1-E_0$ is, denn föllt de $\delta$-Term tatsächlich ganz weg in de baven nömmte Grenz.

Üm dat Argument to vervullständigen, merken wi eerst an, dat dat baven blot de Energiefehler vun den speziellen Tostand $f(H)|\psi\rangle$ is, un nich de Energiefehler vun den neederigsten Energietostand in'n Krylov-Ruum. Aver na dat (Rayleigh-Ritz) variationelle Prinzip is de Energiefehler vun den neederigsten Energietostand in'n Krylov-Ruum na baven begrenzt dörch de Energiefehler vun jeden Tostand in'n Krylov-Ruum, also is dat baven ook 'ne bovenste Grenz för de Energiefehler vun den neederigsten Energietostand, dat hett, de Utgav vun den Krylov-Quantendiagonaliserungsalgorithmus.

'Ne ähnliche Analyse as baven kann dörchmakt warrn, de ook Rusch un de Drempelweert-Prozedur, de in dit Notebook bespraken is, berücksichtigt. Seeh [2] un [4] för disse Analyse.

Anhang: Bewies vun Behauptung 1

Dat Folgende is grötendeels ut [3], Theorem 3.1 afleidt: laat $0 < a < b$ un laat $\Pi^*_d$ de Ruum vun Residualpolynomen sien (Polynomen, de bi 0 den Weert 1 hebbt) vun Graad höchstens $d$. De Lösung vun

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

is

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

un de tohörige minimale Weert is

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

Wi wüllt dat in 'ne Funktion ümwandeln, de natüürlich in Termen vun komplexe Exponentialfunktionen utdrückt warrn kann, wieldat dat de reelle Tiedentwikkelungen sünd, de den Quanten-Krylov-Ruum erzeugt. Üm dat to doon, is dat praktisch, de folgende Transformation vun Energien binnen den Spektralbereich vun den Hamiltonian to Tallen in'n Bereich $[0,1]$ intofohren: defineer

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

wo $dt$ 'n Tiedschritt is, so dat $-\pi < E_0dt < E_\text{max}dt < \pi$. Acht dorop, dat $g(E_0)=0$ un $g(E)$ wasst, wenn $E$ sik vun $E_0$ wegbewegt.

Nu bruken wi dat Polynom $p^*(x)$ mit de Parameter a, b, d sett op $a = g(E_0 + \delta)$, $b = 1$, un d = int(r/2), un definiern de Funktion:

$$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)}$$

wo $E_0$ de Grundtostandsenergie is. Wi köön sehn, indem wi $\cos(x)=\frac{e^{ix}+e^{-ix}}{2}$ insetten, dat $f(E)$ 'n trigonometrisches Polynom vun Graad $d$ is, dat hett, 'ne Linearkombination vun $e^{ijE\,dt}$ för $j=-d,-d+1,...,d-1,d$. Buten dat, ut de Definition vun $p^*(x)$ baven, hebbt wi, dat $f(E_0)=p(0)=1$ un för jeden $E$ in'n Spektralbereich, so dat $\vert E-E_0 \vert > \delta$, hebbt wi

$$|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}.$$

Referenzen

[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).

Tutorial-Ümmfraag

Nimm geern an disse korte Ümmfraag deel, üm Feedback to dit Tutorial to geven. Dien Insichten helpt uns, uns Inholtsanbod un Brukerfarung to verbeten.

Link to survey