- Demos/
- Algorithms/
Exploring Trotterization
Whether we like it or not, time is always moving forward. Even more frustrating, things tend to change with time, meaning we cannot neglect time evolution as we set out to model realistic systems. To execute this successfully, though, tends to be incredibly computationally intense. In the case of particle systems, which constitute much of the interesting simulation work being done in cutting edge science, the exponential nature of the system’s configuration scaling (ie. for :math`n` particles, there are \(2^n\) possible configurations, meaning the Hamiltonian energy matrix of the system used in the time evolution operator \(U(t)=e^{-iHt}\) would be \(2^n \times 2^n\)) very quickly renders classical, sequential simulation impossible, even for incredibly short time scales.
Luckily, quantum computers have shown promise in addressing this issue since, instead of having to represent each possible state contained in the problem’s Hilbert space, qubits can themselves act as analogues for the particles that make up the system, enabling polynomial scaling with particle count rather than exponential scaling. Maybe Feynman was onto something! Carrying out time evolution on these more-efficient systems, however, is not easy, meaning Hamiltonian simulation techniques such as Trotterization need to be implemented to achieve legible and realistic results. This demo will explore how we can do exactly that!
The Commutation Problem
For a completely isolated free particle experiencing no potential energy, the Hamiltonian of the system can be simply defined in terms of kinetic energy and applied to the system via a unitary gate representing \(U_{free}(t)=e^{-iHt}=e^{-iTt}. Sticking to this idealized case would, of course, be useless. With each additional term added to the Hamiltonian, though, the feasibility of the time evolution operator :math:`e^{-iHt}\) being executable on hardware diminishes. This brings to mind an easy fix, why not just split the Hamiltonian into pieces? If we taken the representation \(H=H_1+H_2+...+H_n\), we could naïvely say that our time evolution operator could become \(e^{-i(H_1+H_2+...+H_n)t}=e^{-iH_1t}e^{-iH_2t}...e^{-iH_nt}\). In our naïveté, however, we would neglect to consider the commutation relations that the components of the split Hamiltonian share.
Recall
For commuting matrices, where \(AB=BA\) and it is implied that A and B are completely independent,
\((A+B)^2 = A^2+2AB+B^2\).
For non-commuting matrices, where \(AB \neq BA\) and it is implied that A and B have some level of dependency,
\((A+B)^2 = A^2+AB+BA+B^2\)
So, the Taylor expansion \(e^{A+B}=I+(A+B)+\frac{1}{2}(A+B)^2+...\) differs in the two cases!
As we know, most physical systems require a combination of non-commuting operators to be fully described. Thus, splitting the Hamiltonian this way in a non-commuting scenario would not produce an accurate picture of reality in the time evolution operator. To imagine this more physically, consider what it would mean to exponentiate two dependent properties in this way. Taking A and B to be position and momentum, for example, and letting :math:`e^{iHt}=e^{iAt}e^{iBt}, applying these operators seperately assumes that the system is accurately described by iterating the position then iterating the momentum sequentially. This, of course, ignores the nuance of how the two properties influence each other and renders the simulation inaccurate. Rats!
Luckily for us, this is an issue that some brilliant people gave thought to in order to solve. The Lie-Trotter product formula is the result of such effort. The theorem states that for arbitrary mxm matrices A and B
Turning back to the Hamiltonian picture, if we take a large, finite \(r\), letting \(r\) represent the time step of the evolution, we can approximate the Lie-Trotter formula to
So, by slicing the total time the simulation is trying to emulate, the dependency between non-commuting properties can be integrated via alternating applications of each operator. To turn back to our earlier example, instead of taking one complete position step and one complete momentum step, we are now alternating small, partial steps in position and momentum, simulating simultaneity to the best of our ability.
Implementing the Trotter Method
The form of the Lie-Trotter approximation implies a clear order of operations that must be executed to simulate the time evolution of a non-commuting system.
As implied, the first step of carrying out Hamiltonian simulation using Trotterization is splitting the Hamiltonian. To jump a few steps ahead, implementing the Trotterized time evolution operators will require the construction of equivalent unitary operators using universal gate sets. As a rule of thumb, feasible algorithms should strive for a minimum gate count, meaning we should strive for a minimum number of operators and, therefore, a minumum number of Hamiltonian fragments. To achieve this, the Hamiltonian should be split only where mathematically necessary, meaning commuting terms should always be grouped together to the greatest extent possible. Take, for example, a Hamiltonian containing the terms \(Z_1Z_3\), \(Z_2Z_4\), and \(X_1X_2\), where \(Z_i\) and \(X_j\) are Pauli operators. The first and second terms will always commute since they consist completely of Z operators, but neither of the first two operators commute with the third operator. Therefore, the most optimal splitting would yeild \(H_1=Z_1Z_3+Z_2Z_4\) and \(H_2=X_1X_2\). This is not always so obvious in complex physical systems, so thorough analysis should be carried out to ensure optimal operator pairing.
Once the Hamiltonian is appropriately split, we need to ensure that the fragmented Hamiltonian we are working with can be implemented on realistic hardware. This tends to require some kind of transformation to take place. Recalling the relevance of this method in simulating particle systems, one relevant example is the mapping of Fermionic states to qubit states via the Jordan-Wigner transformation, which translates fermionic creation and annihilation operators to Pauli operators. For the purposes of this demonstration, we will assume our Hamiltonian was originally defined in the Pauli basis and forgo the transformation step.
To carry out the split-time method, the number of time steps \(r\) must be defined, ideally achieving a high level of accuracy while remaining within reasonable computational resource bounds. From there, the circuit should simply alternate between the time evolution operator unitary gates for the desired length of time. Once again considering the Pauli basis, recall that exponentiated Pauli operators are represented by rotation gates, where, letting \(\sigma_i\) be a Pauli operator,
So, taking the Hamiltonian \(H=\alpha X+\beta Z\), where \(\alpha\) and \(\beta\) are probability weights, and recognizing \(H_1=\alpha Z\) and \(H_2=\beta X\)
This is easily translatable to code.
import pennylane as qp
import numpy as np
from scipy.linalg import expm
import pennylane.estimator as qre
from pennylane.transforms.rz_phase_gradient import rz_phase_gradient
#Define System Parameters
coeffs = [0.2, 1.3] #Define Hamiltonian coefficients
observables = [qp.PauliX(0), qp.PauliZ(0)] #Define observables
t = 10 #Define time in seconds
R = [10, 20, 30, 40, 50, 100, 200, 400] #Trotter steps
dev = qp.device("lightning.qubit", wires=1)
#Define Trotterization Function
@qp.qnode(dev)
def TrotterStepper(t,r,coeffs):
del_t = t/r
#Apply the rotation r times
for i in range(r):
U_A = qp.RX(2*coeffs[0]*del_t, wires=0)
U_B = qp.RZ(2*coeffs[1]*del_t, wires=0)
return [qp.expval(qp.PauliX(0)), qp.expval(qp.PauliY(0)), qp.expval(qp.PauliZ(0))]
Trotter Error
Understanding how Hamiltonian simulations deviate from the theory of the system it is trying to recreate is a field of study in itself. Intuitively, we expect the size of the time step being taken to dictate the amount of error in the simulation. A nuanced exploration tends to begin expansion via the Baker-Campbell-Hausdorff formula [#Su2020]_.
As is familiar when handling series expansions, the degree to which the BCH formula is truncated dictates the amount of error to expect in the Hamiltonian. Comparing the previous definition of the approximated Trotter formula to this expansion expression, we can see that we are only concerned with the first term and, therefore, our error is dominated by the term \(-\frac{t^2}{2}[B,A]\) (also known as the First-Order Trotter error). Taking \(t=\Delta t\) for each time step, we can reason out that the error is proportional to \(\frac{t^2}{r}\). This result aligns with our earlier intuition, implying that error will increase with length of time and decrease with number of steps.
In the simple example implemented in this demo, an observed Trotter error can be achieved by calculating time evolution of the system mathematically. This would, of course, be unfeasible for large, complicated models or long time scales, but it is sufficient for this example.
#Approximate evolution for Trotter error calculation
H = qp.Hamiltonian(coeffs, observables)
#Pauli matrix representations
X = np.array([[0, 1], [1, 0]])
Z = np.array([[1, 0], [0, -1]])
Y = np.array([[0, -1j], [1j, 0]])
H = coeffs[0]*X + coeffs[1]*Z
U = expm(-1j*H*t)
state_0 = np.array([1,0])
state_1 = np.array([0,1])
evolved_H = U @ state_0
exact_exp_X = np.real(evolved_H.conj() @ X @ evolved_H)
exact_exp_Y = np.real(evolved_H.conj() @ Y @ evolved_H)
exact_exp_Z = np.real(evolved_H.conj() @ Z @ evolved_H)
By keeping the simulation duration fixed and varying the number of steps taken, the error at each resolution can be obtained.
print(f"{'r':>5} | {'X error':>8} | {'Y error':>8} | {'Z error':>8} | {'Total error':>11}")
print("-" * 51)
for r in R:
result = TrotterStepper(t,r,coeffs)
X_error = abs(result[0]-exact_exp_X)
Y_error = abs(result[1]-exact_exp_Y)
Z_error = abs(result[2]-exact_exp_Z)
total_error = np.sqrt(X_error**2+Y_error**2+Z_error**2)
print(f"{r:>5} | {X_error:>8.5f} | {Y_error:>8.5f} | {Z_error:>8.5f} | {total_error:>11.5f}")
r | X error | Y error | Z error | Total error
---------------------------------------------------
10 | 0.09619 | 0.18017 | 0.00453 | 0.20429
20 | 0.07092 | 0.07839 | 0.00115 | 0.10572
30 | 0.05195 | 0.04863 | 0.00051 | 0.07116
40 | 0.04066 | 0.03498 | 0.00028 | 0.05364
50 | 0.03334 | 0.02723 | 0.00018 | 0.04305
100 | 0.01746 | 0.01282 | 0.00005 | 0.02166
200 | 0.00892 | 0.00620 | 0.00001 | 0.01087
400 | 0.00451 | 0.00305 | 0.00000 | 0.00544
As expected, the simulation error decreases with increasing time steps.
Gate Synthesis Considerations
In fault tolerant architectures, arbitrary rotations must be decomposed into a series of Clifford+T gates via some method of gate synthesis <https://pennylane.ai/compilation/two-qubit-synthesis>_. The PennyLane estimate() tool can be used to determine how many gates are used to implement the synthesized rotation and takes a Repeat Until Success (RUS) approach, which tends to cost \(\mathcal{O}(2.4\log_2(1/\epsilon))\) T-gates per rotation [#Paetznick2014]_. An alternative approach is the phase gradient method (check out the linked demo if you need to brush up on T-counts and quantifying efficiency before going further), which takes advantage of a static register that holds spatially dependent phase values that can be added to a target state as needed. This method has a total cost of \(\mathcal{O}(4log_2(1/\epsilon)+4N)\), where \(N\) is the number of rotations.
When we discuss selecting a gate synthesis method, our goal tends to be to minimize the T-count in favour of implementability and efficiency. From this perspective, it seems obvious to always select the method with the lowest cost, here being the phase gradient approach. Using PennyLane’s rz_phase_gradient() transformation, the expensive \(R_z{\phi}\) rotations implemented in the above Trotterization can be translated into phase gradient additions and compared to the naïve approach. For example, a step count of \(r=200\) will be taken.
#Convert to phase gradient approach
prec = 0.1
b = int(np.ceil(np.log2(1/prec)))
angle_wires = list(range(1,1+b)) #Since we are only dealing with one logical wire, we can just start indexing the auxillary wires at 1
gradient_wires = list(range(1+b,1+2*b))
work_wires = list(range(1+2*b,1+3*b))
dev2 = qp.device("lightning.qubit", wires=(1+3*b))
@qp.qnode(dev2)
@rz_phase_gradient(angle_wires,gradient_wires,work_wires)
def TrotterStepperPG(t,r,coeffs):
del_t = t/r
#Apply the rotation r times
for i in range(r):
U_A = qp.RX(2*coeffs[0]*del_t, wires=0)
U_B = qp.RZ(2*coeffs[1]*del_t, wires=0)
return [qp.expval(qp.PauliX(0)), qp.expval(qp.PauliY(0)), qp.expval(qp.PauliZ(0))]
#Resource estimation
test_r = 200
Trotter_resources = qre.estimate(TrotterStepper)(t,test_r,coeffs) #Repeat Until Success!
Trotter_resources_PG = qre.estimate(TrotterStepperPG)(t,test_r,coeffs)
print("Repeat Until Success")
print(Trotter_resources)
print("Phase Gradient")
print(Trotter_resources_PG)
Repeat Until Success
--- Resources: ---
Total wires: 1
algorithmic wires: 1
allocated wires: 0
zero state: 0
any state: 0
Total gates : 1.760E+4
'T': 1.760E+4
Phase Gradient
--- Resources: ---
Total wires: 16
algorithmic wires: 13
allocated wires: 3
zero state: 3
any state: 0
Total gates : 1.560E+4
'Toffoli': 600,
'T': 8.800E+3,
'CNOT': 4.400E+3,
'Hadamard': 1.800E+3
As expected, the ‘T’-count is much lower in the phase gradient method. This, however, is not a definitive affirmation that the phase gradient method is the ideal gate synthesis method here. Let’s compare the trotter error in the two methods.
print("Repeat Until Success")
print(f"{'r':>5} | {'X error':>8} | {'Y error':>8} | {'Z error':>8} | {'Total error':>11}")
print("-" * 51)
for r in R:
result = TrotterStepper(t,r,coeffs)
X_error = abs(result[0]-exact_exp_X)
Y_error = abs(result[1]-exact_exp_Y)
Z_error = abs(result[2]-exact_exp_Z)
total_error = np.sqrt(X_error**2+Y_error**2+Z_error**2)
print(f"{r:>5} | {X_error:>8.5f} | {Y_error:>8.5f} | {Z_error:>8.5f} | {total_error:>11.5f}")
print("Phase Gradient")
print(f"{'r':>5} | {'X error':>8} | {'Y error':>8} | {'Z error':>8} | {'Total error':>11}")
print("-" * 51)
for r in R:
resultPG = TrotterStepperPG(t,r,coeffs)
X_error_PG = abs(resultPG[0]-exact_exp_X)
Y_error_PG = abs(resultPG[1]-exact_exp_Y)
Z_error_PG = abs(resultPG[2]-exact_exp_Z)
total_error_PG = np.sqrt(X_error_PG**2+Y_error_PG**2+Z_error_PG**2)
print(f"{r:>5} | {X_error_PG:>8.5f} | {Y_error_PG:>8.5f} | {Z_error_PG:>8.5f} | {total_error_PG:>11.5f}")
Repeat Until Success
r | X error | Y error | Z error | Total error
---------------------------------------------------
10 | 0.09619 | 0.18017 | 0.00453 | 0.20429
20 | 0.07092 | 0.07839 | 0.00115 | 0.10572
30 | 0.05195 | 0.04863 | 0.00051 | 0.07116
40 | 0.04066 | 0.03498 | 0.00028 | 0.05364
50 | 0.03334 | 0.02723 | 0.00018 | 0.04305
100 | 0.01746 | 0.01282 | 0.00005 | 0.02166
200 | 0.00892 | 0.00620 | 0.00001 | 0.01087
400 | 0.00451 | 0.00305 | 0.00000 | 0.00544
Phase Gradient
r | X error | Y error | Z error | Total error
---------------------------------------------------
10 | 0.09209 | 0.37899 | 0.18924 | 0.43350
20 | 0.09209 | 0.27772 | 0.16079 | 0.33385
30 | 0.09209 | 0.32937 | 0.20337 | 0.39790
40 | 0.09209 | 0.22663 | 0.19267 | 0.31139
50 | 0.09209 | 0.19622 | 0.09665 | 0.23732
100 | 0.09209 | 0.19863 | 0.09182 | 0.23742
200 | 0.09209 | 0.89700 | 1.63948 | 1.87109
400 | 0.09209 | 0.89700 | 1.63948 | 1.87109
# Ah ha! A tradeoff has made itself clear! In the phase gradient implementation, the Trotter error is univerally higher than in the RUS case. What is also (maybe even more) interesting is that, after a certain :math:`r` threshold is reached, the phase gradient system no longer evolves. This is an example of `underflow <https://en.wikipedia.org/wiki/Arithmetic_underflow>`_ in `quantum arithmetic <https://pennylane.ai/demos/tutorial_how_to_use_quantum_arithmetic_operators>`_, where, put simply, the simulation has reached a computational limit and is stuck rounding to the same value each pass. So, when we choose which techniques to use for gate synthesis, we must consider the needs of our system in tandem with the cost of our system. Sometimes an investment is necessary!
#
# Fast-Fowarding
# --------------
# To make a quantum simulation truly justifiable, one would hope that the required number of time steps would be minimal. If the depth of a circuit maintains proportionality to the length of the time interval being simulated, for example, it runs the risk of exceeding the coherence time of the system. In a completely ideal case, running a simulation for time :math:`t` would require a complexity less than :math:`\mathcal{O}(t)` or, in other words, with a complexity that is sublinear in :math:`t`. This, in theory, can be achieved by employing a "fast-forwarding" technique, which is an umbrella term that refers to techniques used to reduce the depth of a simulation circuit below the time step threshold.
#
# The technique contained under this umbrella that should be used to reduce the depth of a specific simulation must be decided based on the structure of the Hamiltonian itself. If your system has an analytical solution, your life is easy. Fast-forwarding in the analytical case can be carried out by computing the unitary transformation that diagonalizes the system's Hamiltonian and executing it in a quantum circuit. Since the entire evolution of the system is known, the phase angles used in the time-evolution operator (as shown above) can be altered to obtain different time lengths, essentially carrying out simulation with :math:`\mathcal{O}(1)` complexity. This is an incredibly ideal result that applies to a select few, generally not very interesting scenarios, such as a fixed pendulum.
#
# When more complexity is added and the analytical solution of the system is no longer obtainable, we encounter the *no-fast-forwarding theorem* [#Childs2010]_.
#
# .. admonition:: The No-Fast-Forwarding Theorem, Put Simply
# :class: tip
# The time evolution of a generic physical system, often described by a general sparse Hamiltonian, cannot be carried out in sublinear time since there does not exist a generic method for fast-forwarding Hamiltonian simulations.
#
# In the general case, then, one can either accept that the resource requirements are at least :math:`\mathcal{O}(t)` or approximate methods can be employed to carry out fast-forwarding at the cost of some error. One such approximation method is called **Variational Fast Forwarding (VFF)**, in which feedback is used to approximately diagonalize the general Hamiltonian [#Cirstoiu2020]_. In the case of Trotterization, this implies that the unitary describing the Trotter step :math:`e^{-iH\Delta t}` is diagonalized via a gradient descent optimization that targets the diagonalized form. When this is adequately approximated, the diagonal unitary can be applied to the system and fast-forwarding can be carried out, again, by altering the phase angles of the operator [#Cirstoiu2020]_.
#
# .. _references:
#
# References
# ----------
# .. [#Su2020] Y. Su, "A Theory of Trotter Error," presented at the Simons Institute for the Theory of Computing, UC Berkeley, Berkeley, CA, USA, Apr. 2020. [Online]. Available: https://simons.berkeley.edu/sites/default/files/docs/15639/trottererrortheorysimons.pdf
#
# .. [#Paetznick2014] A. Paetznick and K. M. Svore, "Repeat-Until-Success: Non-deterministic decomposition of single-qubit unitaries," *Quantum Inf. Comput.*, vol. 14, no. 15-16, pp. 1277–1301, Nov. 2014, arXiv: 1311.1074 [quant-ph].
#
# .. [#Childs2010] A. M. Childs and R. Kothari, "Limitations on the simulation of non-sparse Hamiltonians," *Quantum Inf. Comput.*, vol. 10, no. 7, pp. 669–684, Jul. 2010, arXiv: `0908.4398 <https://arxiv.org/abs/0908.4398>`_ [quant-ph].
#
# .. [#Cirstoiu2020] C. Cîrstoiu, Z. Holmes, J. Iosue, L. Cincio, P. J. Coles, and A. Sornborger, "Variational fast forwarding for quantum simulation beyond the coherence time," *npj Quantum Inf.*, vol. 6, no. 1, p. 82, Sep. 2020, arXiv: `1910.04292 <https://arxiv.org/abs/1910.04292>`_ [quant-ph].
About the author
Emily Nobes
Emily is a Master of Science in Physics student at McGill University. She works on integrated photonic hardware and tends to dabble in quantum encryption. Emily is joining Xanadu as a resident for summer 2026.
Total running time of the script: (0 minutes 6.020 seconds)