Trotterization tutorial¶
This tutorial notebook goes through the theoretical framework of Trotterization, as well as its implementation in Workbench algorithms. The first section provides a theoretical overview that will be useful to contextualize some of the challenges associated with implementing Trotterization for different applications — those already familiar with Trotterization may wish to skip to section 2 to get to the code implementation.
1. Theoretical overview¶
In quantum algorithms, we are often interested in implementing matrix exponentials of time-independent Hamiltonian operators,
$$ U(t) = \exp(-i H t) \ , $$
where we can express $H$ as a sum of Pauli matrices with real coefficients $H = \sum_{j=1}^m c_j \bigotimes_{k=1}^{n} P_k^{(j)} \ , \ P^{(j)}_k \in \{\mathbb{I}, X, Y, Z\} \ , \ c_j \in \mathbb{R}$.
In general, we do not have a way to directly implement this operation using the native basis gate set of a quantum computer. What we do have access to are Pauli product rotations (PPRs), which implement exponentials of single Pauli product terms as
$$ PPR(\theta) = \exp\left(-i \theta \bigotimes_{j=1}^n P_j\right) \ . $$
The starting point for the problem is the product of PPRs
$$ \tilde{U}(t) = \prod_{j=1}^m \exp\left(-i c_j t \bigotimes_{k=1}^{n} P_j^{(k)}\right) \ . $$
If all the terms in our Hamiltonian mutually commute, then $\tilde{U}(t) = U(t)$ and we are done. However, for many Hamiltonians (read: any Hamiltonians that are actually encoding anything interesting) this is not the case.
Trotterization¶
To begin, let us analyze the result of applying a product of PPRs for a simple Hamiltonian $H = A + B$ (with $A$ and $B$ being Pauli products). The result of applying the PPRs is given by the Baker-Campbell-Haussdorf (BCH) expansion:
$$ e^{-iAt}e^{-iBt} = \exp\left(-it\left(A + B + \frac{1}{2}[A, B] + \frac{1}{12}\left([A, [A, B]] + [B, [B, A]] \right) + \cdots \right)\right) $$
Here $[A, B] := AB - BA$ is the commutator of $A$ and $B$. All the terms apart from the individual $A$ and $B$ operators on the right-hand side are unwanted and are considered to be the Trotter error of the decomposition.
Trotterization is a mathematical trick that replaces the operators $A$ and $B$ with the scaled operators $A/N$ and $B/N$, for some integer $N$ (known as the number of Trotter steps). We can then apply the product of the scaled operator exponentials $N$ times to obtain:
$$ \left(e^{-iAt/N}e^{-iBt/N}\right)^N = \\ \exp\left(-it\left(A + B + \frac{1}{2N}[A, B] - \frac{1}{12 N^2}\left([A, [A, B]] + [B, [B, A]] \right) + \cdots\right)\right) $$
If it isn't clear how this follows, working through this by hand is a useful exercise (you may want to make use of the identity $[A, B] = -[B, A]$).
The main takeaway from doing this is that we are able to reduce the Trotter error by a factor of $N$, with the tradeoff being that the cost of the routine is increased by a factor of $N$. We typically denote this error as scaling with $O((t/N))$, truncating the expansion at the lowest order error term, although it should be noted that this often does not accurately predict the actual error rate of the decomposition, since the higher order commutator terms have non-trivial structure that can result in favorable cancellation, reducing the overall error. Obtaining accurate error predictions from product formulae such as that obtained through Trotterization is well known to be challenging in the literature.
For general Hamiltonians with, say $m$ terms acting on $n$ qubits, we can extend the Trotterization scheme to be
$$ \exp\Bigg(-i t \sum_{j=1}^m c_j \bigotimes_{k=1}^{n} P_j^{(k)}\Bigg) \approx \Bigg(\prod_{j=1}^m \exp\Bigg(-\frac{-ic_jt\bigotimes_{k=1}^n P^{(k)}_j}{N}\Bigg)\Bigg)^N $$
Higher order product formulae — the Suzuki-Trotter decomposition¶
We can exploit the structure and symmetries of the commutators in the BCH formula to improve the asymptotic scaling of the Trotterization protocol. The most common method to do this (and the only one we will discuss in this notebook) is the Suzuki-Trotter decomposition of order $k$. In this framework, the above decomposition would correspond to the first order Suzuki-Trotter decomposition (this framework is so ubiquitous that you will often hear this referred to simply as the "first order Trotterization"). The second order Trotterization is then defined to be the approximation
$$ \exp\Bigg(-i t \sum_{j=1}^m c_j \bigotimes_{k=1}^{n} P_j^{(k)}\Bigg) \approx S_2(t/N)^N := \\ \Bigg(\prod_{j=1}^m \exp\Bigg(-\frac{-ic_jt\bigotimes_{k=1}^n P^{(k)}_j}{2N}\Bigg) \prod_{j=m}^1 \exp\Bigg(-\frac{-ic_jt\bigotimes_{k=1}^n P^{(k)}_j}{2N}\Bigg)\Bigg)^N \ , $$
where the order of the operators matters — in the second product, the order of the operators is the reverse of the order of the first product. The Trotter error for second order Trotterization scales as $O((t/N)^2)$, which is acheieved do to the cancellation of the $\frac{1}{2N}[A, B]$ terms (again, if you're not feeling comfortable with the different product formulae, doing the expansion explicitly for a two-term Hamiltonian could be a useful exercise).
We can obtain higher order product formulae by using the first and second order expansions alongside the recursive formula
$$ S_{2k}(t) = S_{2k-2}(u_kt)^2S_{2k-2}((1-4u_k)t)S_{2k-2}(u_kt)^2 \ , $$
where
$$ u_k := \frac{1}{4-4^{1/(2k-1)}} \ . $$
Technically this denotes a Suzuki decomposition — to obtain the Suzuki-Trotter decomposition, we again divide $t$ by some integer number of Trotter steps $N$ and apply the resulting $S_{2k}(t/N)$ $N$ times. For each order $k$, the asymptotic scaling of the Trotter order is $O((t/N)^k)$, however the cost of applying higher order product formulae scales exponentially in $k$ (at least for this framework — more efficient methods have been proposed (see for example Greatly improved higher-order product formulae for quantum simulation ⧉), but we will not be discussing those here).
2. Workbench implementation¶
In Workbench, Trotterization is implemented using the TrotterQuery qubrick:
from psiqworkbench import QPU, Qubits
from workbench_algorithms import TrotterQuery
from workbench_algorithms.utils import PauliMask, PauliSum
import numpy as np
# set up and reset our QPU
qc = QPU()
qc.reset(2)
# define parameters for the routine
trotter_order = 2
trotter_steps = 3
evolution_time = 1
# initialize the Hamiltonian
hamiltonian = PauliSum(
[np.pi, PauliMask(3, 0)],
[np.pi, PauliMask(3, 3)],
[np.pi, PauliMask(0, 3)],
)
# set up a Qubits register to apply our answer onto
psi = Qubits(2, "psi", qc)
# initialize the Trotter Qubrick
trotter = TrotterQuery(hamiltonian, trotter_order)
# compute the Trotterization
trotter.compute(psi, trotter_steps, evolution_time)
# draw the circuit
qc.draw()
There's a lot going on here, so let's break it down. The first steps are just initializing the QPU and setting the parameters for the Trotterization. If either of these steps aren't familiar, then reading through some earlier tutorials and/or part 1 of this notebook could be useful.
Next we define a PauliSum representation of our Hamiltonian. A detailed overview of these objects is given in the PauliMasks and PauliSums tutorial notebook. For now, just note that for each PauliMask, the first integer corresponds to X terms on the qubits corresponding to the binary expression of that integer and the second corresponds to the Z mask (and when both an X and a Z term are present, the resulting operation is a Y).
Next, we initialize our Trotterization qubrick instance. We pass in the Hamiltonian and the order of the Trotterization, since these will generally be fixed for any instance of the problem. In compute, we pass the number of Trotter steps and the evolution time, since these could, in principle, vary at runtime.
We then draw the circuit. Note that we don't do a gate decomposition of the PPRs, but rather we leave them as basis gates. Ultimately, when we compile our circuits down to implementable gates, our final result will be in terms of PPRs and PPMs, and so it makes sense to go straight to this representation for Trotter. We can directly implement PPRs in Workbench:
from psiqworkbench import QPU
qc = QPU()
qc.reset(5)
qc.ppr(40, x_mask=17, z_mask=22)
print(f"{bin(17)=}")
print(f"{bin(22)=}")
qc.draw()
bin(17)='0b10001' bin(22)='0b10110'
Note that PPRs are defined using the same mask notation as PauliMasks. You should convince yourself that the PPR in the drawing matches what we asked for by referring to the binary values of the X and Z masks.
One of the more confusing aspects of PPRs is the question of how they are controlled. In Workbench, we don't have the ability to directly control PPRs, and instead we must use a decomposition into two PPRs. This is handled automatically by the TrotterQuery qubrick, you shouldn't have to write your own PPRs:
from psiqworkbench import QPU, Qubits
from workbench_algorithms import TrotterQuery
from workbench_algorithms.utils import PauliMask, PauliSum
import numpy as np
# set up and reset our QPU
qc = QPU()
qc.reset(5)
# define parameters for the routine
trotter_order = 2
trotter_steps = 3
evolution_time = 1
# initialize the Hamiltonian
hamiltonian = PauliSum(
[np.pi, PauliMask(3, 0)],
[np.pi, PauliMask(3, 3)],
[np.pi, PauliMask(0, 3)],
)
# set up a Qubits register to apply our answer onto and a control register
psi = Qubits(2, "psi", qc)
ctrl = Qubits(2, "ctrl", qc)
# initialize the Trotter Qubrick
trotter = TrotterQuery(hamiltonian, trotter_order)
# compute the Trotterization, controlled on ctrl
trotter.compute(psi, trotter_steps, evolution_time, ctrl)
# draw the circuit
qc.draw()
Note that there is some units conversion to be handled here. In Workbench, phase angles are given in degrees by default, however physicists typically are more comfortable working in radians. For the example above, the Hamiltonian coefficients (and therefore the rotation angles) are all set to $\pi$. This is then divided by $2*3=6$, since we are using second order Trotterization in this example with 3 Trotter steps to yield $\pi/6=30$ degrees (the conversion is handled under the hood). Note that in the controlled version, the angle is divided by 2 again, with one of the PPRs in each pair negated and extended such that it has an additional $Z$ on the control qubit — a good exercise to check your understanding is to explain why this works.
Trotterization in QPE¶
The result of Trotterization is a unitary operator that encodes the eigenvalues of a Hamiltonian (albeit approximately). It should therefore be of no surprise to anyone that we want to do QPE on this operator (you will get the most out of this tutorial if you're already familiar with QPE — check out that tutorial first if you haven't already).
One of the key advantages of qubricks is that we can directly take code for QPE and plug Trotter directly into it.
A note on examples¶
There is typically a notion that qubitization is a "difficult" subroutine that should be learned after Trotterization, which is a comparatively "easy" subroutine, but there are many conceptual difficulties surrounding Trotterization, in large part due to the fact that the unitary encoding of the Hamiltonian is only approximate, rather than exact as with qubitization. For this reason, tests for QPE can often become muddied by unnecessary details. Rather than get into the weeds of this problem here, we instead show a very simple Hamiltonian consisting of only mutually commuting terms. The mechanics of how you would set up the problem for an actual research task are identical, but the tasks of picking an appropriate order and number of Trotter steps, handling the error analysis and of preparing a valid eigenstate are left to individual research projects.
from workbench_algorithms import QPE, TrotterQuery
from workbench_algorithms.utils import PauliMask, PauliSum, pauli_sum_to_numpy
from psiqworkbench import QPU, Qubits, QUFixed
bits_of_precision = 3
t = 1
steps = 1
trotter_order = 1
val_idx = 2
hamiltonian_type = "Y"
qc = QPU()
qc.reset(bits_of_precision + 2)
phase_bits = [1, 0, 1]
phase = sum(bit / (1 << (i + 1)) for i, bit in enumerate(phase_bits))
print(f"{phase=}")
theta = 2 * np.pi * phase
if hamiltonian_type == "X":
ham = PauliSum(
[theta, PauliMask(1, 0)],
[theta, PauliMask(2, 0)],
[theta, PauliMask(3, 0)],
)
elif hamiltonian_type == "Y":
ham = PauliSum(
[theta, PauliMask(1, 1)],
[theta, PauliMask(2, 2)],
[theta, PauliMask(3, 3)],
)
elif hamiltonian_type == "Z":
ham = PauliSum(
[theta, PauliMask(0, 1)],
[theta, PauliMask(0, 2)],
[theta, PauliMask(0, 3)],
)
norm = ham.norm()
numpy_ham = pauli_sum_to_numpy(ham)
eig_values, eig_vectors = np.linalg.eigh(numpy_ham)
print(f"{eig_values=}")
print(f"{eig_vectors=}")
qreg = Qubits(2, "qreg", qc)
index = QUFixed(bits_of_precision, bits_of_precision, "index", qc)
trotter = TrotterQuery(hamiltonian=ham, trotter_order=trotter_order)
# get the last eigenstate
if hamiltonian_type == "X":
qreg.had()
qreg.z()
elif hamiltonian_type == "Y":
qreg.had()
qreg.s_inv()
else:
qreg.x()
qpe = QPE(bits_of_precision=bits_of_precision, unitary=trotter)
qpe.compute(psi=qreg, phase_reg=index, steps=steps, evo_time=t)
print(f"{index.read()=}")
energy = - index.read() * 2 * np.pi
print(f"Result: {energy}, target: {eig_values[0]}")
qc.draw()
phase=0.625
eig_values=array([-3.92699082, -3.92699082, -3.92699082, 11.78097245])
eig_vectors=array([[-4.00528394e-02+0.00000000e+00j, -8.65098705e-01+0.00000000e+00j,
0.00000000e+00+0.00000000e+00j, -5.00000000e-01+0.00000000e+00j],
[-5.66950282e-01-5.73000537e-01j, 2.62489915e-02+3.15513475e-01j,
-9.49976134e-17-9.18543318e-17j, 0.00000000e+00-5.00000000e-01j],
[ 2.83475141e-01+3.06526688e-01j, -1.31244958e-02+2.74792615e-01j,
-6.13969540e-01-3.50772581e-01j, 3.26633127e-18-5.00000000e-01j],
[-3.06526688e-01+2.83475141e-01j, -2.74792615e-01-1.31244958e-02j,
-3.50772581e-01+6.13969540e-01j, 5.00000000e-01+3.26633127e-18j]])
index.read()=0.625
Result: -3.9269908169872414, target: -3.926990816987242
Using an interface to implement more sophisticated Trotter products¶
The construction of the product of PPRs above is a convenient way of implementing the Trotterization subroutine for general Hamiltonians, but it is not the only way – in fact, the idea of iterating over products of exponential operators is completely general (as long as you have a way of constructing the individual exponentials accurately).
If you have a more complex example problem than the simple Hamiltonian examples outlined above, you can make use of a more general framework for implementing operations: you can implement a TrotterStrategy class.
What is a Strategy?¶
A strategy in this context refers to a software design pattern whereby one class defines an interface (which we refer to as the "strategy", as in "we'll accomplish x using strategy y") for a particular operation, and then accepts any class which implements that interface faithfully at runtime. In this case, rather than coupling the TrotterQuery Qubrick to the PPR implementation above, the actual Qubrick accepts a TrotterStrategy class that implements a apply_trotter_terms method that is called under the hood – this class is reproduced below:
from typing import Protocol, Union
class TrotterStrategy(Protocol):
"""Interface for applying unitary terms in Trotterization."""
def apply_trotter_terms(
self,
factor: Union[int, float, complex],
system: Qubits,
reverse: bool = False,
ctrl: Union[Qubits, int] = 0
):
"""Method describing how to apply the terms in Trotterization, to be implemented in concrete implementations."""
raise NotImplementedError
Don't be daunted by the Protocol inheritance, this is for typing and documentation. The key thing is that the only requirement for implementing a valid strategy is that the class passed to TrotterQuery has a method called apply_trotter_terms that matches the signature above – that's it! No importing or subclassing this class, just implement and go.
Note that in the previous cells we just passed Hamiltonians – if your use case is fairly simple such that everything you need can be conveniently captured by a PauliSum, you are still free to pass this to the TrotterQuery Qubrick. This way we have a balance of flexibility for folks who need it against simplicity for the easiest cases.
The Heisenberg Hamiltonian¶
One use case for this interface is to capture Hamiltonians that have well-defined structures such as spin systems. As a pedagogical example, here's an implementation of a strategy for implementing a Heisenberg model Hamiltonian:
from workbench_algorithms.utils import get_ppr_args_from_ham
class HeisenbergStrategy:
"""Implementation of a TrotterStrategy specifically tailored to Heisenberg Hamiltonians.
Args:
lattice (list): List of pairs of indices that define the lattice. In principle, would be handled robustly using
properly defined lattices, but for the following tests we just apply some random indices.
J (float): Coupling strength between spins in the lattice.
"""
def __init__(self, lattice, J=1):
self.coeffs = [J]
self.lattice = lattice
def hamiltonian_generator(self, reverse=False):
"""Loop through edges, yield appropriate Hamiltonian terms."""
if reverse:
for edge in self.lattice:
yield [self.coeffs[0], get_zz_term(edge)]
for edge in self.lattice:
yield [self.coeffs[0], get_yy_term(edge)]
for edge in self.lattice:
yield [self.coeffs[0], get_xx_term(edge)]
else:
for edge in self.lattice:
yield [self.coeffs[0], get_xx_term(edge)]
for edge in self.lattice:
yield [self.coeffs[0], get_yy_term(edge)]
for edge in self.lattice:
yield [self.coeffs[0], get_zz_term(edge)]
def apply_trotter_terms(self, factor, system, reverse=False, ctrl=0):
"""Apply the Trotterized unitary for the Hamiltonian."""
qc = system.qpu
if ctrl:
anc_mask = ctrl.mask()
for term in self.hamiltonian_generator(reverse=reverse):
base_angle, x_mask, z_mask = get_ppr_args_from_ham(term)
ham_x = system.mask(x_mask)
ham_z = system.mask(z_mask)
if ctrl:
qc.ppr(-factor * base_angle / 2, ham_x, ham_z | anc_mask)
qc.ppr(factor * base_angle / 2, ham_x, ham_z)
else:
qc.ppr(factor * base_angle, ham_x, ham_z)
# utility functions for this Hamiltonian
def get_xx_term(pair):
return PauliMask((1 << pair[0]) + (1 << pair[1]), 0)
def get_zz_term(pair):
return PauliMask(0, (1 << pair[0]) + (1 << pair[1]))
def get_yy_term(pair):
return PauliMask((1 << pair[0]) + (1 << pair[1]), (1 << pair[0]) + (1 << pair[1]))
Now here we're not doing much that's particularly exciting, it's still effectively a PauliSum under the hood, but by moving to a more abstract interface, we get a few advantages:
- We can pass around
HeisenbergStrategyobjects as a way of storing Hamiltonian instances in memory in a more intuitive way than we could using aPauliSum. - We're not storing all the terms of the Hamiltonian in memory (which in this case would be wasteful, since all the coefficients are the same and have a rigid pre-defined structure), we're using a generator to build them on-the-fly.
Although we're not doing it in this example, the flexibility of this approach allows users to do much more sophisticated things like call Qubricks etc.
To show that this works, and to demonstrate the interface, we can repeat our Trotter compute:
from psiqworkbench import QPU, Qubits
from workbench_algorithms import TrotterQuery
from workbench_algorithms.utils import PauliMask, PauliSum
import numpy as np
# set up and reset our QPU
qc = QPU()
qc.reset(5)
# define parameters for the routine
trotter_order = 2
trotter_steps = 3
evolution_time = 1
# initialize the strategy
# a lattice is just a list of pairs of edges
lattice = [(0, 1), (0, 2), (0, 3), (0, 4), (1, 2)]
# then we pass this to the strategy
strategy = HeisenbergStrategy(lattice, J=1)
# set up a Qubits register to apply our answer onto and a control register
psi = Qubits(2, "psi", qc)
ctrl = Qubits(2, "ctrl", qc)
# initialize the Trotter Qubrick
trotter = TrotterQuery(strategy, trotter_order)
# compute the Trotterization, controlled on ctrl
trotter.compute(psi, trotter_steps, evolution_time, ctrl)
# draw the circuit
qc.draw()
