Resource estimation for Hamiltonian simulation with GQSP

Hamiltonian simulation with GQSP

QSP is a method for performing polynomial transformations of block-encoded operators. It consists of a (i) sequence of interleaving signal operators that perform the block-encoding, and (ii) single-qubit signal-processing operators that define the polynomial. GQSP expands on the original approach by considering signal-processing operators that are general \(SU(2)\) transformations; this removes restrictions on available polynomial transformations and facilitates solving for the \(SU(2)\) phase factors, making it the modern method of choice [2].

Hamiltonian simulation is the task of implementing the time-evolution operator \(e^{-iHt}\). We do this for an input Hamiltonian \(H\) and a target evolution time \(t\), subject to a target error \(\varepsilon\). GQSP solves this challenge by expressing the complex exponential \(e^{-iHt}\) as a polynomial, truncated to a degree determined by \(t\) and \(\varepsilon\). The corresponding transformation is then implemented on a block-encoding of \(H\).

In practice, it is customary to instead encode the Hamiltonian as a walk operator \(\hat{W}\), and approximate the function \(e^{-i\cos (\hat{W})t}\) through the rapidly-converging Jacobi-Anger expansion [2]. This type of block-encoding is a qubitization of the Hamiltonian. It can be implemented by a sequence of Prepare and Select operators that are induced by a linear combination of unitaries (LCU) decomposition.

Let’s explore how to use PennyLane to estimate the cost of Hamiltonian simulation with GQSP.

Spin dynamics

We focus on the \(XX\) model Hamiltonian with no external field, defined as

where \(X,Y\) are Pauli matrices acting on a given spin site and \(N\) is the number of spins. The coefficients \(J_{ij}\) can be interpreted as the adjacency matrix of a graph that defines the coupling between spins. Generating samples from a time-evolved state under this Hamiltonian (as well as many other spin models) is widely believed to be an intractable classical problem [4]. But what is the cost of performing this simulation on a quantum computer? For example: how many qubits and gates are needed to perform Hamiltonian simulation with GQSP on a square grid of \(100\times 100\) spins?

To answer this, we first define the Hamiltonian. For the purpose of resource estimation, the specific coupling coefficients are unimportant since the algorithm works identically regardless of their concrete values. This allows us to define compact Hamiltonians that are easy to instantiate by specifying only the type and number of Pauli operators. With periodic boundary conditions, each spin site on the lattice is coupled to four nearest neighbours, so we have 10,000 qubits with 40,000 \(XX\) and \(YY\) couplings respectively. This information is defined as a dictionary and passed directly to the compact PauliHamiltonian class:

import pennylane.estimator as qre
import numpy as np

pauli_dictionary = {"XX": 40000, "YY": 40000}  # 4*100*100

xx_hamiltonian = qre.PauliHamiltonian(
    num_qubits=10000,  # 100*100
    pauli_terms=pauli_dictionary,
)

We now construct the walk operator, which consists of a sequence of Prepare and Select operators. For Prepare, we need extra qubits to load the coefficients, and will employ a standard state preparation algorithm called QROMStatePreparation, based on QROM, which is natively supported in PennyLane:

num_terms = xx_hamiltonian.num_terms  # number of terms in the Hamiltonian
num_state_prep_qubits = int(np.ceil(np.log2(num_terms)))

Prep = qre.QROMStatePreparation(
    num_state_qubits=num_state_prep_qubits,
)

print(f"Resources for Prepare")
print(qre.estimate(Prep))

Resources for Prepare
--- Resources: ---
 Total wires: 80
   algorithmic wires: 17
   allocated wires: 63
     zero state: 63
     any state: 0
 Total gates : 8.529E+6
   'Toffoli': 5.253E+5,
   'CNOT': 5.379E+6,
   'X': 1.049E+6,
   'Z': 32,
   'S': 64,
   'Hadamard': 1.577E+6

For Select, PennyLane provides the SelectPauli resource operator, tailored to Pauli Hamiltonians.

Sel = qre.SelectPauli(xx_hamiltonian)
print(f"Resources for Select")
print(qre.estimate(Sel))

Resources for Select
--- Resources: ---
 Total wires: 1.003E+4
   algorithmic wires: 10017
   allocated wires: 16
     zero state: 16
     any state: 0
 Total gates : 1.040E+6
   'Toffoli': 8.000E+4,
   'CNOT': 3.200E+5,
   'X': 1.600E+5,
   'Z': 8.000E+4,
   'S': 1.600E+5,
   'Hadamard': 2.400E+5

We use Prepare and Select to construct the walk operator, which can be built directly using the dedicated Qubitization operator:

Walk = qre.Qubitization(Prep, Sel)
print(f"Resources for Walk operator")
print(qre.estimate(Walk))

Resources for Walk operator
--- Resources: ---
 Total wires: 1.008E+4
   algorithmic wires: 10017
   allocated wires: 63
     zero state: 63
     any state: 0
 Total gates : 1.810E+7
   'Toffoli': 1.131E+6,
   'CNOT': 1.108E+7,
   'X': 2.257E+6,
   'Z': 8.007E+4,
   'S': 1.601E+5,
   'Hadamard': 3.393E+6

Finally, these pieces are brought together to estimate the cost of performing Hamiltonian simulation with GQSP. This can be calculated directly with the GQSPTimeEvolution operator. Under the hood, it constructs the GQSP sequence and determines the required polynomial degree to simulate the desired dynamics.

As an example, we assume the Hamiltonian is normalized and calculate the degree needed to evolve for

\(t=100\) and a target error of \(\epsilon=0.1\%\).

HamSim = qre.GQSPTimeEvolution(Walk, time=100, one_norm=1, poly_approx_precision=0.001)

print(f"Resources for Hamiltonian simulation with GQSP")
print(qre.estimate(HamSim))

Resources for Hamiltonian simulation with GQSP
--- Resources: ---
 Total wires: 1.008E+4
   algorithmic wires: 10018
   allocated wires: 63
     zero state: 63
     any state: 0
 Total gates : 4.090E+9
   'Toffoli': 2.555E+8,
   'T': 2.996E+4,
   'CNOT': 2.504E+9,
   'X': 5.101E+8,
   'Z': 1.809E+7,
   'S': 3.619E+7,
   'Hadamard': 7.669E+8

This is a large system with non-trivial dynamics, yet we can analyze its requirements in a straightforward manner using PennyLane. We will now study a more practical example applying spin dynamics to NMR spectroscopy.

Heisenberg model for NMR spectral prediction

We follow work presented in Ref. [3] describing quantum simulation of NMR in the zero-to-ultralow field regime. The main computational task is to simulate time evolution under a Heisenberg Hamiltonian of the form

where the sums over \(k,l\) run over spin sites, and the sum over \(\alpha, \beta\) run through spatial directions \(x,y\), and \(z\). We’re using the notation \(\vec{\sigma}\cdot \vec{\sigma} = \sigma_{x}\sigma_{x}+\sigma_{y}\sigma_{y}+\sigma_{z}\sigma_{z}\).

As before, we build a compact Hamiltonian representation by counting the number of Pauli operators of each kind, which is straightforward from expanding the Hamiltonian. For \(N\) spins, sums over \(k\neq l\) run over \(N(N-1)/2\) terms, and there are 9 possible pairs of Pauli matrices. The last sum over \(k\) contains \(N\) terms of 3 different Paulis. We study a system of \(N=32\) spins as in the paper, which gives the following description of the Hamiltonian:

pauli_dictionary = {
    "XX": 496,  # 32*31/2 = 496
    "YY": 496,
    "ZZ": 496,
    "XZ": 496 * 2,  # accounting for "ZX"
    "XY": 496 * 2,  # accounting for "YX"
    "YZ": 496 * 2,  # accounting for "ZY"
    "X": 32,
    "Y": 32,
    "Z": 32,
}

nmr_hamiltonian = qre.PauliHamiltonian(
    num_qubits=32,
    pauli_terms=pauli_dictionary,
)

We build Prepare and Select operators that are used to define the walk operator. For QROMStatePreparation, we enforce a minimal use of auxiliary qubits via the select_swap_depths argument.

num_terms = nmr_hamiltonian.num_terms  # number of terms in the Hamiltonian
num_state_prep_qubits = int(np.ceil(np.log2(num_terms)))

Sel_nmr = qre.SelectPauli(nmr_hamiltonian)
Prep_nmr = qre.QROMStatePreparation(
    num_state_qubits=num_state_prep_qubits,
    positive_and_real=True,
    select_swap_depths=1,  # this minimizes auxiliary qubits
)

Walk_nmr = qre.Qubitization(Prep_nmr, Sel_nmr)

print(f"Resources for NMR Walk operator")
print(qre.estimate(Walk_nmr))

Resources for NMR Walk operator
--- Resources: ---
 Total wires: 108
   algorithmic wires: 45
   allocated wires: 63
     zero state: 63
     any state: 0
 Total gates : 6.071E+5
   'Toffoli': 3.885E+4,
   'CNOT': 3.609E+5,
   'X': 7.461E+4,
   'Z': 3.074E+3,
   'S': 6.144E+3,
   'Hadamard': 1.236E+5

Let’s explore how different choices of evolution time and error affect cost. From theoretical arguments, one would expect linear growth in cost with time, and logarithmic growth with inverse error. We analyze the total number of non-Clifford gates (T+Toffoli) as a function of these parameters:

def nmr_resources(time, one_norm, error):

    gqsp = qre.GQSPTimeEvolution(Walk_nmr, time, one_norm, error)
    resources = qre.estimate(gqsp)
    T_gates = resources.gate_counts["T"]
    Toffoli_gates = resources.gate_counts["Toffoli"]

    return int(T_gates + Toffoli_gates)

We plot the non-Clifford gate cost of the algorithm for different values of total evolution time. This includes two cases where the one-norm differs by a factor of 2, to illustrate the linear increase in cost as a function of one-norm. This is equivalent to rescaling the units of time by a factor of \(\frac{1}{2}\).

import matplotlib.pyplot as plt

one_norm = 1.0
error = 0.001
time_array = np.linspace(1, 1000, 100)  # time from 1 to 1,000
gates = [nmr_resources(t, one_norm, 0.001) for t in time_array]
gates2 = [nmr_resources(t, 2 * one_norm, 0.001) for t in time_array]  # double the one-norm

plt.plot(time_array, gates, label=f"one_norm: {one_norm}")
plt.plot(time_array, gates2, label=f"one_norm: {2*one_norm}")
plt.xlabel("Time")
plt.ylabel("Non-Clifford gates")
plt.grid(True, which="both", linestyle="--")
plt.legend()
plt.show()

Finally, we analyze resources as a function of error to highlight how Hamiltonian simulation with GQSP can rapidly converge to very small errors with minimal overhead.

error_array = np.logspace(2, 9, 100)
gates = [nmr_resources(10, one_norm, 1 / error) for error in error_array]
plt.plot(error_array, gates)
plt.xlabel("Inverse Error")
plt.xscale("log")
plt.ylabel("Non-Clifford gates")
plt.grid(True, which="both", linestyle="--")
plt.show()

Conclusion

Hamiltonian simulation with GQSP is a well-established quantum algorithm with many useful applications. It consists of multiple subroutines that may require time to master, but PennyLane’s quantum resource estimator elegantly aggregates them into easy-to-use operations. This library of operations allows us to both rapidly and accurately quantify resources for specific use cases, and assess their hardware-readiness. Try using estimator yourself, and see how much it would cost for you to study other systems using a quantum computer, such as electronic structure, vibrational, and vibronic Hamiltonians.

References