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 block-encode \(e^{-i\arccos (H)t}\) and approximate the function \(e^{-i\cos (H)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. We refer to this block encoding as a walk operator.

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

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? More specifically: 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 value. 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 PauliHamiltonian resource operator:

import pennylane.estimator as qre
import numpy as np

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

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

print(f"Compact spin Hamiltonian")
print(xx_hamiltonian.pauli_terms)

Compact spin Hamiltonian
{'XX': 40000, 'YY': 40000}

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 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,
    positive_and_real = True  # Can absorb negative coefficients into Pauli terms
)

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 : 4.269E+6
   'Toffoli': 2.631E+5,
   'CNOT': 2.691E+6,
   'X': 5.242E+5,
   'Z': 32,
   'S': 64,
   'Hadamard': 7.900E+5

For Select, PennyLane directly supports a 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 operation

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

Resources for Walk operator
--- Resources: ---
 Total wires: 1.008E+4
   algorithmic wires: 10017
   allocated wires: 63
     zero state: 63
     any state: 0
 Total gates : 9.578E+6
   'Toffoli': 6.063E+5,
   'CNOT': 5.703E+6,
   'X': 1.209E+6,
   'Z': 8.007E+4,
   'S': 1.601E+5,
   'Hadamard': 1.820E+6

Finally, these pieces are brought together to estimate the cost of performing Hamiltonian simulation with GQSP. This can be calculated with the built-in PennyLane function GQSPTimeEvolution. Under the hood, it constructs the GQSP sequence and determines the required polynomial degree in the GQSP transformation 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(W, 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 : 2.165E+9
   'Toffoli': 1.370E+8,
   'T': 2.996E+4,
   'CNOT': 1.289E+9,
   'X': 2.731E+8,
   'Z': 1.809E+7,
   'S': 3.619E+7,
   'Hadamard': 4.113E+8

This is a large system with non-trivial dynamics, yet we can analyze its requirements straightforardly 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

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 6 possible pairs of Pauli matrices when we account for anti-commutation. The last sum over \(k\) contains \(N\) terms of 3 different Paulis. We study a system of \(N=32\) spins, same as in the paper, which gives the following compact Hamiltonian:

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

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

print(f"Compact NMR Hamiltonian")
print(nmr_hamiltonian.pauli_terms)

Compact NMR Hamiltonian
{'XX': 496, 'YY': 496, 'ZZ': 496, 'XZ': 496, 'XY': 496, 'YZ': 496, 'X': 32, 'Y': 32, 'Z': 32}

We build Prepare and Select operators that are used to define the walk operator. For Prepare, PennyLane defaults to choices that minimize gate count at the expense of extra qubits, but this time we enforce a minimal use of ancilla qubits through 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)))

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

W_nmr = qre.Qubitization(Prep_nmr, Sel_nmr)


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

Resources for NMR Walk operator
--- Resources: ---
 Total wires: 107
   algorithmic wires: 44
   allocated wires: 63
     zero state: 63
     any state: 0
 Total gates : 3.200E+5
   'Toffoli': 2.086E+4,
   'CNOT': 1.865E+5,
   'X': 3.887E+4,
   'Z': 2.082E+3,
   'S': 4.160E+3,
   'Hadamard': 6.755E+4

Let’s explore how different choices of evolution time and error affect cost. From theoretical arguments, we expect linear growth in cost with time, and logarithmic increase in inverse error. We build a dedicated function that computes the total number of non-Clifford gates (T+Toffoli) depending on the choice of these parameters. The attribute gate_counts is a dictionary that can be used to extract specific gates.

def nmr_resources(time, one_norm, error):


    gqsp = qre.GQSPTimeEvolution(W_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 1/2.

import matplotlib
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)
plt.plot(time_array, gates2)
plt.xlabel("Time")
plt.ylabel("Non-Clifford gates")
plt.grid(True, which='both', linestyle='--')
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 elegantly aggregates them into easy-to-use operations. This library of operations allows us to rapidly and accurately quantify resources for specific use cases. Users also have the flexibility to customize the algorithm by leveraging the full breadth of capabilities offered as part of PennyLane’s estimator module, for example by constructing custom Prepare and Select operators and studying other systems of interest, such as electronic structure, vibrational, and vibronic Hamiltonians.

References