Skip to main content

Statistical Phase Estimation

In this tutorial, we will go over how to execute statistical phase estimation (SPE) using QURI Algo and QURI VM.

Overview​

SPE is an algorithm that is used to estimate the eigenenergies of a problem Hamiltonian. It does so by performing a Hadamard test on a trial wave-function, resulting in a measurable signal that can be used to infer the wave-functions spectral density. If this wave-function has a significant overlap with the true ground state of the system, the ground state energy should make a significant contribution to the spectral density. The following propoerties makes SPE an attractive algorithm in EFTQC

  • Good circuit depth scaling compared to FTQC algorithms
  • The run-time scales with the energy accuracy ϵ\epsilon as O(ϵ−1)\mathcal{O}(\epsilon^{-1})
  • Noise resilience

In this notebook we will introduce all of the pieces needed to run SPE. In order, this notebook will

  1. Explain the LT22 and Gaussian filter variants of SPE and showcase their implementations in QURI Algo
  2. Instruct in the use of QURI Parts' noise models and showcase the robustness of SPE to noise based on the STAR architecture
  3. Estimate the Fidelity of the Hadamard test circuits used in SPE after transpilation to the STAR architecture and introduce circuit_cost_estimator

Prerequisites​

We do not go into great detail about the following topics, but they will be part of this notebook.

Set up the system​

Before running SPE we briefly introduce some of the building blocks we need. We will here show how to

  • Define a H2\textup{H}_2 molecule
  • Set up programming abstractions needed for SPE using QURI Algo
  • Verify that they behave as expected by conducting the Hadamard test

We first set up the problem using PySCF and quri-parts. PySCF is used to define the molecular geometry and then give us the spin-restricted Hartree-Fock ground state solution to the molecular Hamiltonian. Then the second quantized Hamiltonian is generated by calculating the electron integrals, which is then converted to a qubit operator. For this we use the quri_parts.openfermion and quri_parts.pyscf modules.

from pyscf import gto, scf
from quri_parts.pyscf.mol import get_spin_mo_integrals_from_mole
from quri_parts.openfermion.mol import get_qubit_mapped_hamiltonian
from quri_algo.problem import QubitHamiltonianInput


# Prepare Hamiltonian
mole = gto.M(atom="H 0 0 0; H 0 0 1")
mf = scf.RHF(mole).run(verbose=0)
hamiltonian, mapping = get_qubit_mapped_hamiltonian(
    *get_spin_mo_integrals_from_mole(mole, mf.mo_coeff)
)

Before moving on, we compute some exact system characteristics, e.g. the exact ground state energy and the exact energy gap between the ground state and the first excited state. To this end, we first convert the Hamiltonian to a sparse matrix that we can diagonalize using scipy's sparse linear algebra library.

import numpy as np
from quri_parts.core.operator import get_sparse_matrix

vals, vecs = np.linalg.eigh(get_sparse_matrix(hamiltonian).toarray())

EXACT_GS_ENERGY = np.min(vals)
EXACT_GAP = np.abs(vals[np.argsort(vals)][:2] @ [1, -1])

print("E_0:", EXACT_GS_ENERGY)
print("Delta:", EXACT_GAP)
# output:
    E_0: -1.1011503302326187
    Delta: 0.3552785371970515

To perform further computations, we need to encode our problem into some circuit that represents some operator function of our problem. For example in SPE, to estimate the ground state of a Hamiltonian, we need to encode our Hamiltonian into a controlled time evolution circuit. Here, we choose to do the time evolution with Trotterization. To do this, we wrap the hamiltonian into QubitHamiltonianInput for encoding into a circuit later.

QubitHamiltonianInput can be seen as a generalization of the Operator class from QURI Parts, containing both the terms of a qubit operator and the number of qubits. This makes it convenient when contructing a class hierarchy with quantum circuit factories intended to represent time-evolution and other problem defined quantum circuit constructs.

hamiltonian_input = QubitHamiltonianInput(mapping.n_qubits, hamiltonian)

Circuit factories​

As a direct example of the aforementioned circuit factories we consider the TrotterControlledTimeEvolutionCircuitFactory. As the name suggests this is a circuit factory that instance time-evolution circuits based on the Trotter-Suzuki decomposition. It needs a QubitHamiltonianInput which serves as the generator of time evolution, as well as the number of Trotter steps. Details of this factory are explained in circuit factories.

from quri_parts.circuit.utils.circuit_drawer import draw_circuit
from quri_algo.circuit.time_evolution.trotter_time_evo import (
    TrotterControlledTimeEvolutionCircuitFactory,
)

trotter_concotrolled_time_evo_circuit_factory = (
    TrotterControlledTimeEvolutionCircuitFactory(hamiltonian_input, n_trotter=50)
)

# Time-evolution with t=10
c_time_evo = trotter_concotrolled_time_evo_circuit_factory(evolution_time=1)

# draw_circuit(c_time_evo, line_length=10000)

As seen above, the time-evolution is composed of a series of Pauli rotations. Since we chose a single Trotter step and since it is a controlled time-evolution, each term in the Hamiltonian roughly results in two Pauli rotations. We represent these Pauli rotations as single quantum gates, however to run them on a real device, in most cases, they need to first be transpiled to a series of two-qubit gates and single qubit rotations. Later in this notebook, we will show how this is done to conform to the STAR architecture.

Estimator​

The SPE algorithm requires an esimtation of ⟨e−iHt⟩\langle e^{-iHt} \rangle on different evolution time tt. Follow the estimator tutorial, we build the exact estimator and another one that based on Hadamard test with Trotterized time evolution operator.

from quri_algo.core.estimator.time_evolution_estimator.exact_spectrum import (
    ExactTimeEvolutionExpectationValueEstimator
)
from quri_algo.core.estimator.time_evolution_estimator.trotter import (
    TrotterTimeEvolutionHadamardTest
)
from quri_parts.qulacs.sampler import create_qulacs_vector_sampler

sampler = create_qulacs_vector_sampler()
exact_time_evolution_estimator = ExactTimeEvolutionExpectationValueEstimator(
    hamiltonian_input, vals, vecs
)
trotter_time_evolution_estimator = TrotterTimeEvolutionHadamardTest(
    hamiltonian_input, sampler, 50
)

We then build a state as an example. We choose a computational basis superposition state.

from quri_parts.core.state import quantum_state, comp_basis_superposition

state = comp_basis_superposition(
    quantum_state(mapping.n_qubits, bits=0b0011),
    quantum_state(mapping.n_qubits, bits=0b1100),
    7 / 20,
    0,
)

Let's compute ⟨e−iHt⟩\langle e^{-iHt} \rangle with the state above. First we use the exact estimator, which is based on matrix multiplication

print(
    f"The time-evolution expectation value is {exact_time_evolution_estimator(state, evolution_time=5.).value}"
)
print(
    f"The trotter time-evolution expectation value is {trotter_time_evolution_estimator(state, evolution_time=5., n_shots=int(1e6)).value}"
)
# output:
    The time-evolution expectation value is (0.7803970841187445-0.5734363633936669j)
    The trotter time-evolution expectation value is (0.780776-0.573478j)

As we shall see below, this Hadamard test is essential to SPE. The tools provided above can be used to dictate whether the Hadamard test should be done exactly or with a Trotterized quantum circuit implementation based on sampling.

Executing statistical phase estimation​

SPE has two main variants, which we will explore in this notebook. The first one is by Lin and Tong (LT22) which introduced SPE in the context of early fault tolerant devices for the first time. We also explore a variant of it by Wang et. al. which we refer to as Gaussian filter. The chief difference between the two is that the signal that is generated uses a convolution function that is either a Heaviside step functino (LT22) or a Gaussian (Gaussian filter).

Introducing the signal and convolution functions​

SPE tries to obtain the spectrum of a Hamiltonian H=∑iλi∣λi⟩⟨λi∣H = \sum_{i} \lambda_i \ket{\lambda_i} \bra{\lambda_i} by using a trial state ∣ψ⟩|\psi\rangle. ∣ψ⟩|\psi\rangle`s overlap with the eigenstates ∣λi⟩|\lambda_i\rangle of the Hamiltonian HH is pi=∣⟨λi∣ψ⟩∣2p_{i} = |\langle \lambda_i| \psi\rangle|^2. In terms of this overlap the spectral density function is

ϱ(x)=∑ipiδ(x−λiτ),\varrho(x) = \sum_{i} p_i\delta(x-\lambda_i\tau),

where δ(x)\delta(x) is the Dirac delta function and τ\tau is a dimensionless scaling parameter applied to the Hamiltonian, that we introduce through the algorithm itself.

Depending on the concrete implementation of SPE, we may never directly resolve the spectral density function exactly, but we construct a signal that will have features in common with it. This is done by first picking a convolution function F(x)F(x) such that when the spectral density ϱ(x)\varrho(x) is convolved by F(x)F(x), the eigenvalues λi\lambda_i manifest as special points in the resulting signal

(F∗ϱ)(x)=∫−TTF(y)ϱ(y−x)dy=∑ipiF(x−λiτ).(F * \varrho)(x) = \int_{-T}^T F(y)\varrho(y-x)dy =\sum_{i}p_{i}F(x-\lambda_i \tau).

Here TT defines the support of F(x)F(x), while τ\tau determines where the eigenvalues are laid out in the resulting signal. We here skip the full derivation, but the interested reader is directed to our paper summarizing these ideas while providing benchmarks of SPE. Here it suffices to say that we can express this convolution in terms of the Fourier components F~\tilde{F} as

(F∗ϱ)(x)≈∑n=0N−1eiknxF~(kn)⟨e−iknτH⟩.(F * \varrho)(x)\approx \sum_{n=0}^{N-1}e^{ik_n x} \widetilde{F}(k_n) \langle e^{-ik_n\tau H}\rangle.

Here NN is a parameter that gives the upper bound on the number of Fourier modes included in the algorithm. From this expression and the preceeding discussion we see that this signal can be recovered using a quantum processor. The expectation value of the time-evolution operator can be evaluated using the Hadamard test. The only other component we need is a Fourier transform of the convolution function. The convolution function itself can be chosen depending on desired algorithm scaling.

The convolution function is only approximately represented in the convolution integral because of the cut-off NN. This defines a trade-off between the accuracy of the convolution function and the run-time needed to accurately represent time-evolution by kNτk_N \tau.

Finally, the scaling of the integrand with F~(kn)\tilde{F}(k_n) is not carried out by direct multiplication. Instead F~(kn)\tilde{F}(k_n) is sampled to determine the number of samples used to evaluate ⟨e−iknτH⟩\langle e^{-ik_n\tau H}\rangle.

In the following we will show how to proceed with QURI Algo's SPE implementation with the LT22 variant is chosen in which FF is a Heaviside step-function.

The LT22 algorithm​

We start with the LT22 variant. First we have to specify the maximum index of the Fourier modes considered, δ\delta, a parameter that controls the interval within which the error of the Fourier transformed Heaviside step function is bounded. In addition, we specify τ\tau. As mentioned, τ\tau effectively scales first axis of the resulting signal, thus acting as a parameter for fine-graining the signal. The consequence of choosing a small τ\tau is that the signal itself may not include the ground state energy or other eigenenergies of interest, however, choosing it too large will render the resolution too low to effectively distinguish different peaks in the signal.

max_fourier = 1000
delta = 1e-4
n_sample = 10000
tau = 1 / 20
from quri_algo.algo.phase_estimation.spe import StepFunctionParam
from quri_algo.algo.phase_estimation.spe.lt22 import SingleSignalLT22GSEE
from pprint import pprint

eta = 0.4

signal_param = StepFunctionParam(d=1000, delta=1e-4, n_sample=10000)

lt22_algorithm = SingleSignalLT22GSEE(trotter_time_evolution_estimator, tau=tau)

spe_result = lt22_algorithm(state, signal_param, eta)
lt22_gs_energy = spe_result.phase / tau
print(f"The obtained ground state energy is {lt22_gs_energy}")
print(f"The error is = {abs(lt22_gs_energy - np.min(vals))}")
# output:
    The obtained ground state energy is -1.1067624476952391
    The error is = 0.0056121174626204695

TmaxT_{\text{max}} and TtotalT_{\text{total}} can be computed by passing in the list of signal functions to the SPEGSEERecorder.

from quri_algo.algo.phase_estimation.spe.utils.recorder import SPEGSEERecorder

resource = SPEGSEERecorder(spe_result.signal_functions, tau)
pprint(resource)
# output:
    SPEGSEERecorder(max_evolution_time=np.float64(49.95),
                    total_evolution_time=np.float64(102955.90000000005),
                    n_shots=np.int64(20000))

When the algorithm is executed, the signal is automatically saved in the algorithm object, where we can plot out the signal at any derivative order

import matplotlib.pyplot as plt

xs = np.linspace(-np.pi / 10, np.pi / 10, 10000)
signal_func = spe_result.signal_functions[0]
ys = np.abs(signal_func(xs, derivative_order=0))
yds = np.abs(signal_func(xs, derivative_order=1))

fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot(
    xs,
    ys,
    lw=3,
    label=r"$Z_F(x; \tau)$",
)
ax1.axvline(np.min(vals) * tau, color="red", linestyle="--", label="Exact")
ax1.legend(fontsize=20)

ax2.plot(xs, yds, lw=3, label=r"$Z'_F(x; \tau)$")
ax2.axvline(np.min(vals) * tau, color="red", linestyle="--", label="Exact")
ax2.legend()
ax2.legend(fontsize=20)

for ax in [ax1, ax2]:
    ax.tick_params(axis="both", which="major", labelsize=16)

fig.set_size_inches(20, 8)
plt.show()

png

As seen above, SPE with the step function convolution function results in fairly distinguishable signatures of the spectral density function peaks. These appear as steps in the raw signal and as peaks in the first order derivative. If the calculation is carried out exactly, the first order derivative is seen to result in the actual spectral density. However, due to sampling noise and the limit on Fourier modes the resulting signal carries significant noise.

That being said, for the above signal these errors do not appear to significantly alter the centroids of the resulting peaks. As long as the overlap is such that each peak retains a height well above the noise floor it appears that LT22 is robust to noise.

Gaussian​

The Gaussian filter variant of SPE is overall similar to LT22 and the flow is the same apart from some classical post-processing steps. However, Gaussian filter is sensitive to the parameter choices made since the eigenenergies manifest as Gaussians in the returned signal. The standard deviation of each Gaussian for instance influences the energy resolution possible. Below we define a function that returns optimal metaparameters based on the paper by Wang et. al.

from typing import Optional
from quri_algo.algo.phase_estimation.spe import GaussianParam
from quri_algo.algo.phase_estimation.spe.gaussian import GaussianFittingGSEE
from quri_algo.algo.phase_estimation.spe.utils.fourier.gaussian_coefficient import (
    GaussianSampler,
)


def get_recommended_gaussian_parameter(
    gap: float,
    eps: float,
    overlap: float,
    delta: float,
    n_discretize: int,
    max_shot_limit: Optional[int] = None,
) -> tuple[GaussianParam, int]:
    Delta = gap
    sigma = np.min(
        [
            0.9 * Delta / np.sqrt(2 * np.log(9 * Delta * eps**-1 * overlap**-1)),
            0.2 * Delta,
        ]
    )
    eps_tilde = 0.1 * eps * overlap / (np.sqrt(2 * np.pi) * sigma**3)
    T = (
        np.pi**-1
        * sigma**-1
        * np.sqrt(2 * np.log(8 * np.pi**-1 * eps_tilde**-1 * sigma**-2))
    )
    M = np.ceil(sigma / eps) + 1
    fk = GaussianSampler(T, n_discretize, sigma).fourier_coefficients
    shot_const = np.abs(fk).sum()
    S = int(
        shot_const**2 * np.ceil(np.log(4 * M / (delta / 2)) / (eps_tilde / 2) ** 2)
    )
    if max_shot_limit is not None:
        S = np.min([S, max_shot_limit])
    return (
        GaussianParam(T, n_discretize, sigma, S),
        M,
    )

In order to use this, we must provide an upper bound on the excitation gap of the system. Although we can obtain it fairly easily for the above H2\textup{H}_2 Hamiltonian, we do not in general expect to be able to solve the problem exactly. So we will instead use the coupled cluster excitation energy as a rough estimate. We can calculate it by using the equation of motion CCSD using PySCF

from pyscf import cc

ccsd = cc.CCSD(mf).run()
e_ee, _ = ccsd.eeccsd(nroots=2)
print("EOMCCSD excitation energy is:", e_ee[0])
# output:
    E(CCSD) = -1.101150330244479  E_corr = -0.03504168092654223
    EOMCCSD excitation energy is: 0.3552785372089112


    <class 'pyscf.cc.ccsd.CCSD'> does not have attributes  converged

We now obtain the parameters

gaussian_param, n_scanned_energy = get_recommended_gaussian_parameter(
    e_ee[0],
    eps=1e-6,
    overlap=0.5,
    delta=1e-4,
    n_discretize=1000,
    max_shot_limit=1e7,
)
pprint(gaussian_param)
# output:
    GaussianParam(T=np.float64(31.331137453515634),
                  N=1000,
                  sigma=np.float64(0.057114676949872885),
                  n_sample=np.float64(10000000.0))

Before launching the Gaussian filter variation of SPE let's see if we can get a rough idea of what the eigenenergies are with LT22 and then based on that refine the search criteria of the Gaussian filter method.

rough_signal_param = StepFunctionParam(d=100, delta=1e-4, n_sample=10000)
rough_energy = lt22_algorithm(state, rough_signal_param, eta=0.5).phase / tau
print(f"Rough energy estimation by LT22 = {rough_energy}")
print(f"Energy error = {rough_energy - np.min(vals)}")
# output:
    Rough energy estimation by LT22 = -1.1673438044412443
    Energy error = -0.06619347420862565

This should help narrow down the interval to search for the true ground state energy.

gaussian_spe_algorithm = GaussianFittingGSEE(trotter_time_evolution_estimator, tau)
search_range = gaussian_spe_algorithm.get_recommended_search_range(
    rough_energy * tau, n_scanned_energy, gaussian_param.sigma, 10
)
gaussian_result = gaussian_spe_algorithm(state, gaussian_param, search_range, eta=0.5)
gaussian_energy = gaussian_result.phase / tau
print("Gaussian energy:", gaussian_energy)
print("Energy error:", np.abs(gaussian_energy - EXACT_GS_ENERGY))
# output:
    Gaussian energy: -0.8452554910894953
    Energy error: 0.2558948391431234
gaussian_resource = SPEGSEERecorder(gaussian_result.signal_functions, tau)
pprint(gaussian_resource)
# output:
    SPEGSEERecorder(max_evolution_time=np.float64(0.7550804126297271),
                    total_evolution_time=np.float64(2223559.6273274785),
                    n_shots=np.int64(20000000))

As opposed to LT22, the Gaussian filter method ideally results in a signal which contains a set of Gaussians that are centered on each eigenvalue. Here we have done enough pre-processing to narrow down the centroid quite a lot and will only focus on the ground state energy. The Gaussian will be flat at its peak, and because of that we can obtain the ground state energy by finding the zeros of the first derivative of the signal.

xs = search_range
signal_func = gaussian_result.signal_functions[0]
ys = np.abs(signal_func(xs, derivative_order=0))
yds = np.real(signal_func(xs, derivative_order=1))

fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot(
    xs,
    ys,
    lw=3,
    label=r"$Z_F(x; \tau)$",
)
ax1.axvline(np.min(vals) * tau, color="red", linestyle="--", label="Exact")
ax1.legend(fontsize=20)

ax2.plot(xs, yds, lw=3, label=r"$Z'_F(x; \tau)$")
ax2.axvline(np.min(vals) * tau, color="red", linestyle="--", label="Exact")
ax2.legend()
ax2.legend(fontsize=20)
# plt.hlines()

for ax in [ax1, ax2]:
    ax.tick_params(axis="both", which="major", labelsize=16)

fig.set_size_inches(20, 8)
fig.suptitle(
    "$T_{max}$="
    f"{gaussian_resource.max_evolution_time: .2f}, "
    "$T_{total}$="
    f"{gaussian_resource.total_evolution_time: .2e}, "
    "$N_{shots}$="
    f"{gaussian_resource.n_shots: .2e}",
    size=20,
)
plt.show()

png

Executing statistical phase estimation on noisy devices​

Finally, we show how to obtain estimates of the performance in a noisy environment. Here we will use transpilers from QURI Parts to optimize the circuit somewhat. Most of these implement some minor circuit optimizations, such as removing gates or combinations of gates that amount to identity operations, but some of them are essential for device simulation, such as the decomposition of Pauli rotations to actual gates implemntable by a quantum processor. This is done through PauliRotationDecomposeTranspiler. Finally we combine them into a single SequentialTranspiler. The transpiler combination we choose is tailored to the STAR architecture, which is an EFTQC architecture with partial error correction. Later we evaluate the cost of running the time-evolution circuit needed for SPE on the STAR architecture.

from quri_vm import VM
from quri_parts.backend.devices import star_device
from quri_parts.backend.units import TimeUnit, TimeValue

logical_qubit_count = 5
p_phys = 1e-4
code_distance = 9
qec_cycle = TimeValue(1, TimeUnit.MICROSECOND)

star_property = star_device.generate_device_property(
    logical_qubit_count, code_distance, qec_cycle, p_phys
)

star_vm = VM.from_device_prop(prop=star_property)

The circuit factory for the time-evolution accepts a transpiler argument which is applied to the circuit. With this in mind we instantiate the circuit factory below and obtain the resulting set of gates with

star_evolution_circuit_factory = TrotterControlledTimeEvolutionCircuitFactory(
    hamiltonian_input, n_trotter=50, transpiler=star_vm.transpile
)

print("The set of gates used by SPE after transpilation to the STAR architecture is")
set([g.name for g in star_evolution_circuit_factory(10).gates])
# output:
    The set of gates used by SPE after transpilation to the STAR architecture is
    {'CNOT', 'H', 'RZ', 'S'}

Noise model​

For the noise model we choose a simplified model with the assumption that the Clifford gates, which undergo surface code error correction under the STAR architecture, can be assumed to be ideal, while the non-Clifford Rz gate is subject only to post-selection, whereby the resulting phase flip error probability based on device error probability pphysp_\textup{phys} becomes:

PL=(1−2pphys15)2P_L = \left(1 - \frac{2 p_{\textup{phys}}}{15} \right)^2

Noisy LT22 simulation​

Finally, we perform the noisy simulation using QURI Parts

rough_signal_param = StepFunctionParam(d=100, delta=1e-4, n_sample=10000)

noisy_trotter_time_evolution_estimator = TrotterTimeEvolutionHadamardTest(
    encoded_problem=hamiltonian_input,
    sampler=star_vm.sample,
    n_trotter=50,
    transpiler=star_vm.transpile,
)

Comparison with the exact time evolution and the noise-free Trotterized time-evolution shows that noise has a significant effect on the estimation of the time-evolution unitary expectation value.

evolution_time = 5.0
print(
    "Exact time evolution <e^{-iHt}>:",
    exact_time_evolution_estimator(state, evolution_time).value,
)
print(
    "Trotter time evolution <e^{-iHt}>:",
    trotter_time_evolution_estimator(state, evolution_time, 10000).value,
)
print(
    "Noisy Trotter time evolution <e^{-iHt}>:",
    noisy_trotter_time_evolution_estimator(state, evolution_time, 10000).value,
)
# output:
    Exact time evolution <e^{-iHt}>: (0.7803970841187445-0.5734363633936669j)
    Trotter time evolution <e^{-iHt}>: (0.7732-0.56j)
    Noisy Trotter time evolution <e^{-iHt}>: (0.7792-0.5636j)

Based on this inaccuracy it is interesting to see what the resulting signal will look like when used in SPE. All we have to do is pass the noisy trotter time evolution estimator to the LTE algorithm.

noisy_lt22 = SingleSignalLT22GSEE(noisy_trotter_time_evolution_estimator, tau=tau)
noisy_lt22_result = noisy_lt22(state, rough_signal_param, eta=0.5)
noisy_lt22_energy = noisy_lt22_result.phase / tau
print("Noisy LT22 energy:", noisy_lt22_energy)
print("Energy error:", noisy_lt22_energy - EXACT_GS_ENERGY)
# output:
    Noisy LT22 energy: -1.1031857807105505
    Energy error: -0.002035450477931855
noisy_lt22_resource = SPEGSEERecorder(noisy_lt22_result.signal_functions, tau)
pprint(noisy_lt22_resource)
# output:
    SPEGSEERecorder(max_evolution_time=np.float64(4.95),
                    total_evolution_time=np.float64(13016.700000000003),
                    n_shots=np.int64(20000))
xs = np.linspace(-np.pi / 10, np.pi / 10, 1000)
signal_func = noisy_lt22_result.signal_functions[0]
ys = np.abs(signal_func(xs, derivative_order=0))
yds = np.real(signal_func(xs, derivative_order=1))

fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot(
    xs,
    ys,
    lw=3,
    label=r"$Z_F(x; \tau)$",
)
ax1.axvline(np.min(vals) * tau, color="red", linestyle="--", label="Exact")
ax1.legend(fontsize=20)

ax2.plot(xs, np.abs(yds), lw=3, label=r"$Z'_F(x; \tau)$")
ax2.axvline(np.min(vals) * tau, color="red", linestyle="--", label="Exact")
ax2.legend()
ax2.legend(fontsize=20)

for ax in [ax1, ax2]:
    ax.tick_params(axis="both", which="major", labelsize=16)

fig.set_size_inches(20, 8)
fig.suptitle(
    "$T_{max}$="
    f"{noisy_lt22_resource.max_evolution_time: .2f}, "
    "$T_{total}$="
    f"{noisy_lt22_resource.total_evolution_time: .2e}, "
    "$N_{shots}$="
    f"{noisy_lt22_resource.n_shots: .2e}",
    size=20,
)
plt.show()

png

This result is remarkable as it shows that even in an environment with significant noise, it is still possible to accurately determine the ground state energy.

Resource estimation and circuit evaluation​

The noise simulations carried out so far have been based on a simplified model of the partial error correction of the STAR architecture. While it is enough for demonstration purposes we can utilize a much more accurate device characteristic by using the cost_estimator. First we import the device objects we need, in this case star_device, but also the estimators. We will be estimating the circuit fidelity and latency, so we are going to use estimate_circuit_fidelity and estimate_circuit_latency respectively. For convenience we define the function get_star_cost_fidelity below which returns both of these estimates.

In the SPE implementation we have used, the time-evolution circuits have all been implemented with the same number of Trotter steps, the only difference between them resulting from different angles chosen for the Pauli rotations that were since transpiled into a combination of CNOT, H, S and Rz gates. Because of this we need only consider a single one of these circuits. We note that the Hadamard test is executed twice to get both the real and imaginary components of ⟨e−ikjτH⟩\braket{e^{-ik_j\tau H}}, and because of this the execution time is doubled.

spe_samples = noisy_lt22_result.signal_functions[0].spe_samples

evo_time = spe_samples[0].classical_sample.k * tau
circuit = star_evolution_circuit_factory(evo_time)
star_analysis_result = star_vm.analyze(circuit)
star_fidelity = star_analysis_result.fidelity
star_latency = star_analysis_result.latency.in_ns()
star_latency *= 1e-9

print("Hadamard test circuit fidelity:", star_fidelity)
print("circuit execution time per circuit:", star_latency)
total_execution_time = (
    star_latency * len(spe_samples) * 2
)  # factor of 2 for sampling real and imaginary parts
print(f"total execution time: {total_execution_time} seconds")
print(f"energy error: {(noisy_lt22_energy - EXACT_GS_ENERGY) * 1000: .2f} mHa")
# output:
    Hadamard test circuit fidelity: 0.9565507633471438
    circuit execution time per circuit: 0.19710000000000003
    total execution time: 39.81420000000001 seconds
    energy error: -2.04 mHa

Based on the above we obtain the circuit execution time, as well as the energy error and circuit fidelity.

Summary​

In this notebook we showed how to

  • Execute SPE with QURI Algo and QURI Parts
  • Transpile SPE circuits for the STAR architecture
  • Build a simple noise model for the STAR architecture in QURI Parts and simulate SPE using it
  • Build a faithful noise model and use it to estimate run-time as well as fidelity

SPE is an algorithm that appears promising in the context of EFTQC. Using it, it appears that we can make use of its noise robustness in the EFTQC era on devices such as the STAR architecture.