Diagonalisasi kuantum Krylov untuk Hamiltonian kisi

Estimasi penggunaan: 20 menit di Heron r2 (CATATAN: Ini hanya estimasi. Waktu aktual bisa berbeda.)

# 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

Latar Belakang

Tutorial ini menunjukkan cara mengimplementasikan Algoritma Diagonalisasi Kuantum Krylov (KQD) dalam konteks pola Qiskit. Kamu akan mempelajari teori di balik algoritma ini terlebih dahulu, lalu melihat demonstrasi eksekusinya di QPU.

Di berbagai disiplin ilmu, kita tertarik untuk mempelajari properti keadaan dasar sistem kuantum. Contohnya meliputi pemahaman sifat fundamental partikel dan gaya, memprediksi dan memahami perilaku material kompleks, serta memahami interaksi dan reaksi biokimia. Karena pertumbuhan eksponensial ruang Hilbert dan korelasi yang muncul dalam sistem yang terjalin, algoritma klasik kesulitan memecahkan masalah ini untuk sistem kuantum yang semakin besar. Di satu sisi spektrum terdapat pendekatan yang memanfaatkan perangkat keras kuantum yang berfokus pada metode kuantum variasional (misalnya, variational quantum eigensolver). Teknik-teknik ini menghadapi tantangan pada perangkat saat ini karena jumlah pemanggilan fungsi yang tinggi dalam proses optimasi, yang menambah overhead besar dalam sumber daya begitu teknik mitigasi kesalahan tingkat lanjut diterapkan, sehingga membatasi efektivitasnya pada sistem kecil. Di sisi lain spektrum terdapat metode kuantum toleran terhadap kesalahan dengan jaminan kinerja (misalnya, quantum phase estimation), yang memerlukan sirkuit dalam yang hanya bisa dieksekusi di perangkat toleran kesalahan. Oleh karena itu, kami memperkenalkan algoritma kuantum berbasis metode subruang (seperti yang dijelaskan dalam makalah tinjauan ini), yaitu algoritma diagonalisasi kuantum Krylov (KQD). Algoritma ini bekerja baik pada skala besar [1] di perangkat keras kuantum yang ada, memiliki jaminan kinerja yang serupa dengan phase estimation, kompatibel dengan teknik mitigasi kesalahan tingkat lanjut, dan dapat memberikan hasil yang tidak bisa diakses secara klasik.

Persyaratan

Sebelum memulai tutorial ini, pastikan kamu sudah menginstal yang berikut:

• Qiskit SDK v2.0 atau lebih baru, dengan dukungan visualisasi
• Qiskit Runtime v0.22 atau lebih baru ( pip install qiskit-ibm-runtime )

Pengaturan

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

Langkah 1: Petakan input klasik ke masalah kuantum

Ruang Krylov

Ruang Krylov $\mathcal{K}^r$ berorde $r$ adalah ruang yang dibentangkan oleh vektor-vektor yang diperoleh dengan mengalikan pangkat-pangkat lebih tinggi dari matriks $A$, hingga $r-1$, dengan vektor referensi $\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\}$$

Jika matriks $A$ adalah Hamiltonian $H$, ruang yang bersangkutan disebut ruang Krylov daya $\mathcal{K}_P$. Dalam kasus di mana $A$ adalah operator evolusi waktu yang dihasilkan oleh Hamiltonian $U=e^{-iHt}$, ruangnya disebut ruang Krylov uniter $\mathcal{K}_U$. Subruang Krylov daya yang kita gunakan secara klasik tidak dapat dibangkitkan langsung di komputer kuantum karena $H$ bukan operator uniter. Sebagai gantinya, kita bisa menggunakan operator evolusi waktu $U = e^{-iHt}$ yang dapat ditunjukkan memberikan jaminan konvergensi serupa dengan metode daya. Pangkat-pangkat $U$ kemudian menjadi langkah waktu yang berbeda $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\}$$

Lihat Lampiran untuk derivasi rinci tentang bagaimana ruang Krylov uniter memungkinkan representasi keadaan eigen berenergi rendah secara akurat.

Algoritma diagonalisasi kuantum Krylov

Diberikan Hamiltonian $H$ yang ingin kita diagonalisasikan, pertama kita pertimbangkan ruang Krylov uniter $\mathcal{K}_U$ yang bersesuaian. Tujuannya adalah menemukan representasi kompak dari Hamiltonian dalam $\mathcal{K}_U$, yang akan kita sebut $\tilde{H}$. Elemen matriks dari $\tilde{H}$, proyeksi Hamiltonian dalam ruang Krylov, dapat dihitung dengan menghitung nilai ekspektasi berikut

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

Di mana $\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle$ adalah vektor-vektor dari ruang Krylov uniter dan $t_n = n dt$ adalah kelipatan langkah waktu $dt$ yang dipilih. Di komputer kuantum, perhitungan setiap elemen matriks dapat dilakukan dengan algoritma apa pun yang memungkinkan untuk mendapatkan tumpang tindih antara keadaan kuantum. Tutorial ini berfokus pada uji Hadamard. Karena $\mathcal{K}_U$ memiliki dimensi $r$, Hamiltonian yang diproyeksikan ke subruang akan berdimensi $r \times r$. Dengan $r$ yang cukup kecil (umumnya $r<<100$ sudah cukup untuk mendapatkan konvergensi estimasi eigenenergy), kita kemudian dapat dengan mudah mendiagonalisasikan Hamiltonian yang diproyeksikan $\tilde{H}$. Namun, kita tidak bisa langsung mendiagonalisasikan $\tilde{H}$ karena ketidakortogonalan vektor-vektor ruang Krylov. Kita harus mengukur tumpang tindihnya dan membangun matriks $\tilde{S}$

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

Ini memungkinkan kita memecahkan masalah nilai eigen dalam ruang non-ortogonal (juga disebut masalah nilai eigen umum)

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

Selanjutnya kita bisa mendapatkan estimasi nilai eigen dan keadaan eigen dari $H$ dengan melihat yang ada di $\tilde{H}$. Misalnya, estimasi energi keadaan dasar diperoleh dengan mengambil nilai eigen terkecil $c$ dan keadaan dasar dari vektor eigen yang bersesuaian $\vec{c}$. Koefisien dalam $\vec{c}$ menentukan kontribusi vektor-vektor berbeda yang membentangkan $\mathcal{K}_U$.

Gambar tersebut menampilkan representasi Circuit dari uji Hadamard yang dimodifikasi, sebuah metode yang digunakan untuk menghitung tumpang tindih antara berbagai keadaan kuantum. Untuk setiap elemen matriks $\tilde{H}_{i,j}$, uji Hadamard antara keadaan $\vert \psi_i \rangle$, $\vert \psi_j \rangle$ dilakukan. Ini disorot dalam gambar dengan skema warna untuk elemen matriks dan operasi $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$ yang bersesuaian. Dengan demikian, serangkaian uji Hadamard untuk semua kemungkinan kombinasi vektor ruang Krylov diperlukan untuk menghitung semua elemen matriks dari Hamiltonian yang diproyeksikan $\tilde{H}$. Kawat atas dalam Circuit uji Hadamard adalah Qubit ancilla yang diukur baik dalam basis X maupun Y, nilai ekspektasinya menentukan nilai tumpang tindih antara keadaan-keadaan tersebut. Kawat bawah mewakili semua Qubit dari Hamiltonian sistem. Operasi $\text{Prep} \; \psi_i$ mempersiapkan Qubit sistem dalam keadaan $\vert \psi_i \rangle$ yang dikontrol oleh keadaan Qubit ancilla (demikian pula untuk $\text{Prep} \; \psi_j$) dan operasi $P$ mewakili dekomposisi Pauli dari Hamiltonian sistem $H = \sum_i P_i$. Derivasi lebih rinci tentang operasi yang dihitung oleh uji Hadamard diberikan di bawah.

Definisikan Hamiltonian

Mari kita pertimbangkan Hamiltonian Heisenberg untuk $N$ Qubit pada rantai linier: $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)]

Atur parameter untuk algoritma

Kita secara heuristis memilih nilai untuk langkah waktu dt (berdasarkan batas atas dari norma Hamiltonian). Ref [2] menunjukkan bahwa langkah waktu yang cukup kecil adalah $\pi/\vert \vert H \vert \vert$, dan lebih baik sampai suatu titik untuk meremehkan nilai ini daripada melebih-lebihkannya, karena melebih-lebihkan dapat membiarkan kontribusi dari keadaan berenergi tinggi merusak bahkan keadaan optimal dalam ruang Krylov. Di sisi lain, memilih $dt$ yang terlalu kecil menghasilkan pengkondisian subruang Krylov yang lebih buruk, karena vektor basis Krylov berbeda lebih sedikit dari langkah waktu ke langkah waktu.

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

Dan atur parameter lain dari algoritma. Demi tutorial ini, kita akan membatasi diri pada penggunaan ruang Krylov dengan hanya lima dimensi, yang cukup terbatas.

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

Persiapan keadaan

Pilih keadaan referensi $\vert \psi \rangle$ yang memiliki suatu tumpang tindih dengan keadaan dasar. Untuk Hamiltonian ini, kita menggunakan keadaan dengan eksitasi di Qubit tengah $\vert 00..010...00 \rangle$ sebagai keadaan referensi kita.

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

Evolusi waktu

Kita dapat merealisasikan operator evolusi waktu yang dihasilkan oleh Hamiltonian tertentu: $U=e^{-iHt}$ melalui aproksimasi 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>

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

Di mana $P$ adalah salah satu suku dalam dekomposisi Hamiltonian $H=\sum P$ dan $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$ adalah operasi terkontrol yang mempersiapkan vektor $|\psi_i\rangle$, $|\psi_j\rangle$ dari ruang Krylov uniter, dengan $|\psi_k\rangle = e^{-i H k dt } \vert \psi \rangle = e^{-i H k dt } U_{\psi} \vert 0 \rangle^N$. Untuk mengukur $X$, pertama terapkan $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*}$$

... kemudian ukur:

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

Dari identitas $|a + b\|^2 = \langle a + b | a + b \rangle = \|a\|^2 + \|b\|^2 + 2\text{Re}\langle a | b \rangle$. Demikian pula, mengukur $Y$ menghasilkan

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

Circuit uji Hadamard bisa menjadi Circuit yang dalam setelah kita dekomposisi ke gate native (yang akan semakin dalam jika kita memperhitungkan topologi perangkat)

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

Langkah 2: Optimalkan masalah untuk eksekusi perangkat keras kuantum

Uji Hadamard yang efisien

Kita bisa mengoptimalkan Circuit yang dalam untuk uji Hadamard yang telah kita peroleh dengan memperkenalkan beberapa aproksimasi dan mengandalkan beberapa asumsi tentang Hamiltonian model. Misalnya, perhatikan Circuit berikut untuk uji Hadamard:

Asumsikan kita bisa menghitung $E_0$ secara klasik, nilai eigen dari $|0\rangle^N$ di bawah Hamiltonian $H$. Hal ini terpenuhi ketika Hamiltonian mempertahankan simetri U(1). Meskipun ini mungkin tampak seperti asumsi yang kuat, ada banyak kasus di mana aman untuk mengasumsikan bahwa ada keadaan vakum (dalam hal ini dipetakan ke keadaan $|0\rangle^N$) yang tidak terpengaruh oleh aksi Hamiltonian. Ini berlaku misalnya untuk Hamiltonian kimia yang mendeskripsikan molekul stabil (di mana jumlah elektron dipertahankan). Mengingat bahwa gate $\text{Prep} \; \psi$, mempersiapkan keadaan referensi yang diinginkan $\ket{psi} = \text{Prep} \; \psi \ket{0} = e^{-i H 0 dt} U_{\psi} \ket{0}$, misalnya, untuk mempersiapkan keadaan HF dalam kimia $\text{Prep} \; \psi$ adalah sebuah produk dari NOT Qubit tunggal, sehingga controlled-$\text{Prep} \; \psi$ hanyalah sebuah produk dari CNOT. Kemudian Circuit di atas mengimplementasikan keadaan berikut sebelum pengukuran:

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

di mana kita telah menggunakan pergeseran fase yang bisa disimulasikan secara klasik $U\ket{0}^N = e^{i\phi}\ket{0}^N$ pada baris ketiga. Oleh karena itu nilai-nilai ekspektasi diperoleh sebagai

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

Dengan menggunakan asumsi-asumsi ini kita berhasil menuliskan nilai-nilai ekspektasi dari operator yang diminati dengan lebih sedikit operasi terkontrol. Bahkan, kita hanya perlu mengimplementasikan persiapan keadaan terkontrol $\text{Prep} \; \psi$ dan bukan evolusi waktu terkontrol. Merumuskan ulang perhitungan kita seperti di atas akan memungkinkan kita mengurangi kedalaman Circuit yang dihasilkan secara signifikan.

Dekomposisi operator evolusi waktu dengan dekomposisi Trotter

Alih-alih mengimplementasikan operator evolusi waktu secara tepat, kita bisa menggunakan dekomposisi Trotter untuk mengimplementasikan aproksimasinya. Mengulangi beberapa kali dekomposisi Trotter dengan orde tertentu memberikan kita pengurangan lebih lanjut dari kesalahan yang diperkenalkan oleh aproksimasi. Berikut ini, kita langsung membangun implementasi Trotter dengan cara paling efisien untuk graf interaksi Hamiltonian yang sedang kita pertimbangkan (hanya interaksi tetangga terdekat). Dalam praktiknya kita menyisipkan rotasi Pauli $R_{xx}$, $R_{yy}$, $R_{zz}$ dengan sudut berparameter $t$ yang sesuai dengan implementasi aproksimasi dari $e^{-i (XX + YY + ZZ) t}$. Mengingat perbedaan dalam definisi rotasi Pauli dan evolusi waktu yang kita coba implementasikan, kita harus menggunakan parameter $2*dt$ untuk mencapai evolusi waktu sebesar $dt$. Selanjutnya, kita membalik urutan operasi untuk jumlah pengulangan langkah Trotter yang ganjil, yang secara fungsional setara namun memungkinkan sintesis operasi yang bersebelahan dalam satu uniter $SU(2)$. Ini menghasilkan Circuit yang jauh lebih dangkal daripada yang diperoleh menggunakan fungsionalitas generik 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)

Gunakan Circuit yang dioptimalkan untuk persiapan keadaan

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)

Template Circuit untuk menghitung elemen matriks dari $\tilde{S}$ dan $\tilde{H}$ melalui uji Hadamard

Satu-satunya perbedaan antara Circuit yang digunakan dalam uji Hadamard adalah fase dalam operator evolusi waktu dan observabel yang diukur. Oleh karena itu kita bisa mempersiapkan Circuit template yang mewakili Circuit generik untuk uji Hadamard, dengan placeholder untuk gate yang bergantung pada operator evolusi waktu.

# 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

Kita telah mengurangi kedalaman uji Hadamard secara signifikan dengan kombinasi aproksimasi Trotter dan uniter yang tidak terkontrol

Langkah 3: Eksekusi menggunakan Qiskit primitives

Buat instance backend dan atur parameter runtime

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

Transpilasi ke QPU

Pertama, mari kita pilih subset dari coupling map dengan qubit yang berkinerja "baik" (di sini "baik" cukup arbitrer, yang paling penting kita ingin menghindari qubit dengan performa sangat buruk) dan membuat target baru untuk transpilasi

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

Kemudian transpilasi Circuit virtual ke layout fisik terbaik pada target baru ini

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]

Buat PUB untuk eksekusi dengan 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)

Jalankan Circuit

Circuit untuk $t=0$ dapat dihitung secara klasik

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)

Eksekusi Circuit untuk $S$ dan $\tilde{H}$ dengan 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])

Langkah 4: Pasca-proses dan kembalikan hasil dalam format klasik yang diinginkan

results = job.result()[0]

Hitung matriks Hamiltonian Efektif dan Overlap

Pertama hitung fase yang terakumulasi oleh keadaan $\vert 0 \rangle$ selama evolusi waktu yang tidak terkontrol

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

Setelah kita memiliki hasil eksekusi Circuit, kita dapat memproses data untuk menghitung elemen matriks $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]$$

Dan elemen matriks dari $\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]$$

Terakhir, kita dapat menyelesaikan masalah nilai eigen umum untuk $\tilde{H}$:

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

dan mendapatkan estimasi energi ground state $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

Untuk sektor partikel tunggal, kita dapat menghitung ground state dari sektor Hamiltonian ini secara efisien secara klasik

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

Lampiran: Krylov subspace dari evolusi waktu nyata

Krylov space uniter didefinisikan sebagai

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

untuk suatu langkah waktu $dt$ yang akan kita tentukan nanti. Sementara itu asumsikan $r$ genap: kemudian definisikan $d=r/2$. Perhatikan bahwa ketika kita memproyeksikan Hamiltonian ke Krylov space di atas, hasilnya tidak dapat dibedakan dari Krylov space

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

yakni ketika semua evolusi waktu digeser mundur sebesar $d$ langkah waktu. Alasan mengapa tidak dapat dibedakan adalah karena elemen matriks

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

tidak berubah di bawah pergeseran keseluruhan waktu evolusi, karena evolusi waktu berkomutasi dengan Hamiltonian. Untuk $r$ ganjil, kita dapat menggunakan analisis untuk $r-1$.

Kita ingin menunjukkan bahwa di suatu tempat dalam Krylov space ini, dijamin ada keadaan berenergi rendah. Kita lakukan dengan cara berikut, yang diturunkan dari Teorema 3.1 dalam [3]:

Klaim 1: terdapat fungsi $f$ sedemikian rupa sehingga untuk energi $E$ dalam rentang spektrum Hamiltonian (yaitu, antara energi ground state dan energi maksimum)...

1. $f(E_0)=1$
2. $|f(E)|\le2\left(1 + \delta\right)^{-d}$ untuk semua nilai $E$ yang berada $\ge\delta$ jauh dari $E_0$, artinya nilainya disupresi secara eksponensial
3. $f(E)$ adalah kombinasi linier dari $e^{ijE\,dt}$ untuk $j=-d,-d+1,...,d-1,d$

Kita berikan bukti di bawah ini, tetapi bagian itu dapat dilewati jika tidak ingin memahami argumen yang lengkap dan ketat. Untuk saat ini kita fokus pada implikasi dari klaim di atas. Berdasarkan properti 3, kita dapat melihat bahwa Krylov space yang telah digeser di atas mengandung keadaan $f(H)|\psi\rangle$. Inilah keadaan berenergi rendah kita. Untuk memahami alasannya, tuliskan $|\psi\rangle$ dalam basis eigenenergy:

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

di mana $|E_k\rangle$ adalah eigenstate energi ke-k dan $\gamma_k$ adalah amplitudonya dalam keadaan awal $|\psi\rangle$. Dinyatakan dalam hal ini, $f(H)|\psi\rangle$ diberikan oleh

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

menggunakan fakta bahwa kita dapat mengganti $H$ dengan $E_k$ ketika ia bekerja pada eigenstate $|E_k\rangle$. Kesalahan energi dari keadaan ini oleh karena itu adalah

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

Untuk mengubah ini menjadi batas atas yang lebih mudah dipahami, pertama-tama kita pisahkan jumlah dalam pembilang menjadi suku dengan $E_k-E_0\le\delta$ dan suku dengan $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}.$$

Kita dapat membatasi suku pertama dari atas dengan $\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,$$

di mana langkah pertama mengikuti karena $E_k-E_0\le\delta$ untuk setiap $E_k$ dalam jumlah, dan langkah kedua mengikuti karena jumlah dalam pembilang adalah bagian dari jumlah dalam penyebut. Untuk suku kedua, pertama kita batasi penyebut dari bawah dengan $|\gamma_0|^2$, karena $f(E_0)^2=1$: menggabungkan semuanya kembali bersama, ini memberikan

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

Untuk menyederhanakan sisanya, perhatikan bahwa untuk semua $E_k$ ini, berdasarkan definisi $f$ kita tahu bahwa $f(E_k)^2 \le 4\left(1 + \delta\right)^{-2d}$. Selain itu dengan membatasi $E_k-E_0<2\|H\|$ dari atas dan $\sum_{E_k>E_0+\delta}|\gamma_k|^2<1$ dari atas memberikan

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

Ini berlaku untuk $\delta>0$ apa pun, sehingga jika kita menetapkan $\delta$ sama dengan kesalahan target kita, maka batas kesalahan di atas konvergen ke sana secara eksponensial terhadap dimensi Krylov $2d=r$. Perhatikan juga bahwa jika $\delta<E_1-E_0$ maka suku $\delta$ sebenarnya hilang sepenuhnya dari batas di atas.

Untuk melengkapi argumen, pertama-tama kita catat bahwa yang di atas hanyalah kesalahan energi dari keadaan khusus $f(H)|\psi\rangle$, bukan kesalahan energi dari keadaan berenergi terendah dalam Krylov space. Namun, berdasarkan prinsip variasional (Rayleigh-Ritz), kesalahan energi dari keadaan berenergi terendah dalam Krylov space dibatasi dari atas oleh kesalahan energi dari keadaan mana pun dalam Krylov space, sehingga yang di atas juga merupakan batas atas pada kesalahan energi dari keadaan berenergi terendah, yaitu, output dari algoritma Krylov quantum diagonalization.

Analisis serupa seperti di atas dapat dilakukan yang juga memperhitungkan noise dan prosedur thresholding yang dibahas dalam notebook. Lihat [2] dan [4] untuk analisis ini.

Lampiran: Bukti Klaim 1

Berikut ini sebagian besar diturunkan dari [3], Teorema 3.1: misalkan $0 < a < b$ dan misalkan $\Pi^*_d$ adalah ruang polinomial residual (polinomial yang nilainya di 0 adalah 1) dengan derajat paling banyak $d$. Solusi untuk

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

adalah

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

dan nilai minimal yang bersesuaian adalah

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

Kita ingin mengubah ini menjadi fungsi yang dapat diekspresikan secara alami dalam bentuk eksponen kompleks, karena itulah evolusi waktu nyata yang menghasilkan Krylov space kuantum. Untuk melakukannya, akan mudah jika kita memperkenalkan transformasi energi berikut dalam rentang spektrum Hamiltonian ke angka dalam rentang $[0,1]$: definisikan

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

di mana $dt$ adalah langkah waktu sedemikian rupa sehingga $-\pi < E_0dt < E_\text{max}dt < \pi$. Perhatikan bahwa $g(E_0)=0$ dan $g(E)$ bertumbuh seiring $E$ menjauh dari $E_0$.

Sekarang menggunakan polinomial $p^*(x)$ dengan parameter a, b, d diset ke $a = g(E_0 + \delta)$, $b = 1$, dan d = int(r/2), kita definisikan fungsi:

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

di mana $E_0$ adalah energi ground state. Kita dapat melihat dengan memasukkan $\cos(x)=\frac{e^{ix}+e^{-ix}}{2}$ bahwa $f(E)$ adalah polinomial trigonometri berderajat $d$, yaitu kombinasi linier dari $e^{ijE\,dt}$ untuk $j=-d,-d+1,...,d-1,d$. Selanjutnya, dari definisi $p^*(x)$ di atas kita memiliki bahwa $f(E_0)=p(0)=1$ dan untuk $E$ apa pun dalam rentang spektrum sedemikian rupa sehingga $\vert E-E_0 \vert > \delta$ kita miliki

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

Referensi

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

Survei tutorial

Silakan isi survei singkat ini untuk memberikan masukan tentang tutorial ini. Wawasan Anda akan membantu kami meningkatkan konten dan pengalaman pengguna.

Tautan ke survei

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.