Exploring Trotterization

The commutation problem

For a completely isolated free particle experiencing no external potential, the system’s Hamiltonian can be defined in terms of solely kinetic energy and applied via a unitary gate representing the time evolution operator \(U_{free}(t)=e^{-iHt}\), where \(T\) is the kinetic energy term. As a Hamiltonian becomes more complex (e.g., additional terms are added), the chance that the time evolution operator \(e^{-iHt}\) can be implemented with available gates and computational resources diminishes.

If size is the issue, though, why not just split the Hamiltonian into smaller pieces that are executable and hardware compatible? If we take the representation

where \(H_n\) is an implementable Hamiltonian fragment, we could naïvely say that our time evolution operator becomes

In our naïveté, however, we would neglect to consider the commutation relations shared by the split Hamiltonian components. If certain fragments of the Hamiltonian do not commute (\([H_i,H_j] \neq 0\)), then this product would merely be an approximation that is likely not adequately representative of the system we are trying to simulate, resulting in systemic error and inaccurate outputs.

As we know, most physical systems require a combination of non-commuting operators to be fully described. If we took \(A\) and \(B\) to be non-commuting operators, our naïve approach would lead us to blindly exponentiate as

Applying these operators sequentially assumes that the system is accurately described by iterating \(A\) then iterating \(B\), one after the other. This completely ignores the influence the two operators have on each other, implied by their non-commutativity. Position and momentum, for example, do not evolve in time independently. Rats.

Luckily for us, the problem of exponentiating non-commuting equations is not new. The Lie-Trotter product formula is a well known means of navigating this issue. The Lie-Trotter theorem states that, for arbitrary \(m \times m\) matrices \(A\) and \(B\),

where \(r\) represents the number of time steps taken in the desired time evolution. Turning back to the Hamiltonian simulation picture, if we take a large, finite \(r\) we can approximate the Lie-Trotter formula to the first order Trotterization expression

So, by slicing the total time the simulation is trying to emulate, the dependency shared by non-commuting properties can be integrated via alternating applications of each operator within each time step. In other words, instead of taking one complete \(A\) step and one complete \(B\) step, we are now alternating small, partial steps in \(A\) and \(B\), approximating simultaneity to the best of our ability.

By splitting the Hamiltonian into fragments \(H_j\), we are generating a new effective Hamiltonian that becomes an approximation of the initial Hamiltonian (a true perspective is that \(e^{-i(A+B)t}\approx e^{-iAt}e^{-iBt}\)). Carrying out the above Trotterization allows us to execute an exact simulation of an approximate representation!

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.

The first step of Trotterization is splitting the Hamiltonian. As a rule of thumb, feasible algorithms should strive for a minimal gate count and reduced opportunity for error, meaning we should try to minimize the number of operators required. 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.

Once the Hamiltonian is appropriately split, we need to ensure that the fragments we are working with can actually be implemented using realistic quantum devices. This tends to require some kind of transformation between the native representation of the system and a computationally compatible representation. One relevant example is the mapping of Fermionic states to qubit states via the jordan_wigner() transformation. For the purposes of the following simple demonstration, we will assume our Hamiltonian was originally defined in the Pauli basis and forgo the transformation step.

Let our Hamiltonian be given by

where \(\alpha\) and \(\beta\) are arbitrary coefficients.

Recall that exponentiated Pauli operators are represented by rotation gates, where, letting \(\sigma_i\) be a Pauli operator,

So, taking \(H_1=\alpha X\) and \(H_2=\beta Z\)

Finally, the number of time steps \(r\) must be defined with the goal of 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 number of steps.

This is easily translated 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 [alpha, beta]
t = 10
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, coeffs, r):
    del_t = t / r

    # Apply the rotation r times
    for i in range(r):
        qp.RX(2 * coeffs[0] * del_t, wires=0)
        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))]


print(TrotterStepper(t, coeffs, R[4]))

[0.12542579222970107, -0.11296005772224645, 0.9856512547564343]

Luckily for us, PennyLane has the tools to make this much simpler. Using trotterize(), the exact same procedure can be carried out on the target Hamiltonian.

def FirstOrderExpansion(time, theta, phi, wires=None):
    qp.RX(2 * time * theta, wires=0)
    qp.RZ(2 * time * phi, wires=0)


@qp.qnode(dev)
def BuiltInTrotter(time, theta, phi, num_trotter_steps):
    qp.trotterize(FirstOrderExpansion, n=num_trotter_steps, order=1)(time, theta, phi, wires=[0])
    return [qp.expval(qp.PauliX(0)), qp.expval(qp.PauliY(0)), qp.expval(qp.PauliZ(0))]


# Compare results for 50 Trotter steps
print(BuiltInTrotter(t, coeffs[0], coeffs[1], R[4]))

[0.12542579222970107, -0.11296005772224645, 0.9856512547564343]

That’s a match! Another option would be to use TrotterProduct, which is an operator class that handles the fragmentation and exponentiation itself. These tools are advantageous for dealing with complex Hamiltonians or working Trotterization into a larger PennyLane workflow. For the purposes of this demonstration, we will continue using our DIY solution, but the choice is yours.

For this simple Hamiltonian, fragmentation is handled by definition. In more complex scenarios, fragmentation must be carried out as its own preliminary step. 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 yield \(H_1=Z_1Z_3+Z_2Z_4\) and \(H_2=X_1X_2\). This is often not obvious in complex physical systems, so thorough analysis should be carried out to ensure optimal operator pairing. Tools such as group_observables() are very helpful in this case.

Trotter error

Understanding how Hamiltonian simulations deviate from theoretically expected system behaviour is a developing field of research. Intuitively, we expect the size of the time step to be a major contributor to the error, since time steps that are too large tend to obscure system behaviour. A nuanced exploration of exact Trotter error typically begins with consideration of the Baker-Campbell-Hausdorff formula [1], which describes the true expansion of an exponentiated non-commuting operator pair as

As is familiar when handling series expansions, the degree to which the BCH formula is truncated dictates the amount of error to expect. 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 second-order term \(-\frac{t^2}{2}[B,A]\). Taking \(t=\Delta t=\frac{t}{r}\) for each time step, we can reason that, in this case, the error is proportional to \(\frac{t^2}{r}\) after the operator is applied \(r\) times.

In the simple example implemented in this demo, an observed Trotter error can be achieved by calculating the time evolution of the system analytically and carrying out a comparison. This would, of course, be unfeasible for large, complicated (in other words, useful) models or long time scales, but it is sufficient for this example.

# Approximate evolution for Trotter error calculation
# 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])

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, coeffs, r)
    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.

Higher-order Trotterizations

Not too long ago we mentioned that we were working with the BCH formula truncated to the first order term. In other words, we have been taking a first-order Trotterization approach. While this is sufficient for simple, short time scale problems (e.g., systems with weak interaction), ignoring higher order terms and, therefore, reducing precision can result in ignorance of important system characteristics and unnecessarily high error that requires increased resources to account for [5]. Thus, the use of higher-order Trotterizations is often necessary to achieve realistic simulation.

A simple, baseline approach to achieving high-order Trotterizations is introducing symmetries into the system that move us closer to reality. Revisiting the analogy used above, the first-order Trotterization alternates between two non-commuting operators and incrementally steps each term in alternating time slices. In a second-order approach, one of the operator terms is further divided into two half-steps to be applied before and after the other operator. In the case of our Hamiltonian, the second-order Trotter product, discovered simultaneously by Strang [6] and Verlet [7] in the 1960s, would be

which can be implemented as

A benefit of imposing a symmetric formula such as this is that it cancels out the dominant error term discussed previously. The symmetry, essentially, does not allow for even powers to survive the operator application process, eliminating the dominant \(t^2\) term.

Altering the TrotterStepper() function to a second-order Trotterization shows the impact this symmetry has on the Trotter error.

@qp.qnode(dev)
def TrotterStepperSO(t, coeffs, r):
    del_t = t / r

    # Apply the rotation r times
    for _ in range(r):
        qp.RX(coeffs[0] * del_t, wires=0)
        qp.RZ(2 * coeffs[1] * del_t, wires=0)
        qp.RX(coeffs[0] * del_t, wires=0)

    return [qp.expval(qp.PauliX(0)), qp.expval(qp.PauliY(0)), qp.expval(qp.PauliZ(0))]


print(f"{'r':>5} | {'X error':>8} | {'Y error':>8} | {'Z error':>8} | {'Total error':>11}")
print("-" * 51)

for r in R:
    resultSO = TrotterStepperSO(t, coeffs, r)
    X_errorSO = abs(resultSO[0] - exact_exp_X)
    Y_errorSO = abs(resultSO[1] - exact_exp_Y)
    Z_errorSO = abs(resultSO[2] - exact_exp_Z)
    total_errorSO = np.sqrt(X_errorSO**2 + Y_errorSO**2 + Z_errorSO**2)
    print(
        f"{r:>5} | {X_errorSO:>8.5f} | {Y_errorSO:>8.5f} | {Z_errorSO:>8.5f} | {total_errorSO:>11.5f}"
    )

r |  X error |  Y error |  Z error | Total error
---------------------------------------------------
   10 |  0.06783 |  0.09452 |  0.01283 |     0.11705
   20 |  0.01819 |  0.02249 |  0.00446 |     0.02927
   30 |  0.00815 |  0.00988 |  0.00208 |     0.01298
   40 |  0.00460 |  0.00554 |  0.00119 |     0.00730
   50 |  0.00295 |  0.00354 |  0.00077 |     0.00467
  100 |  0.00074 |  0.00088 |  0.00019 |     0.00117
  200 |  0.00018 |  0.00022 |  0.00005 |     0.00029
  400 |  0.00005 |  0.00006 |  0.00001 |     0.00007

Comparing to the first-order results where, for example, the error when \(r=10\) is approximately 20%, updating to the second-order Trotterization reduces this error to approximately 12%. Not too shabby!

Beyond second-order, higher-order Trotterizations can be achieved via the nested application of the second-order Trotterization sequence. A well-known example is the Suzuki five-step ladder [8], which achieves a fourth-order implementation represented as

Where \(s_1=(4-4^{1/3})^{-1}\) is the first-order Suzuki constant. It should be noted that this is not the only method of achieving higher-order representations and that various methods exist with the goal of reducing Trotter error, minimizing gate count, and achieving high accuracy. To further optimize, some approaches allow the selective application of high-order Trotter products to dominant terms in a simulation, allowing for resources to be allocated to only the instances they are justified.

Gate synthesis considerations

Before we can officially add Trotterization to our list of capabilities, we must consider how our simulations will interface with eventual fault-tolerant quantum hardware. As upsetting as it can be, the resources provided by quantum hardware will be, in some way, limited. With this in mind, it is important that the resource requirements of the systems we work with are quantified and evaluated. One metric of interest is the number of T-gates required to synthesize our operations.

The PennyLane estimate() tool can be implemented to approximate the number of gates used to carry out a given process, taking a naïve approach. Alternatively, the multiplexed phase gradient method, which takes advantage of a static register that holds spatially dependent phase values which can be added to a target state as needed, can be used for synthesis.

When we select a gate synthesis method, our initial instinct is to minimize the T count in favour of implementability and efficiency above all. From this perspective, it seems obvious to always select a low-cost method, such as the aforementioned 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. Since this method only transforms \(R_z\) operations, we can perform our desired \(R_x\) rotations using an \(R_z\) rotation sandwiched between two Hadamard gates. For this example, a step count of \(r=200\) will be taken in the resource estimation step.

# Convert to phase gradient approach
prec = 0.05
b = int(np.ceil(np.log2(1 / prec)))

angle_wires = list(range(1, 1 + b))
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, coeffs, r):
    del_t = t / r

    # Apply the rotation r times
    for i in range(r):
        # Perform X rotation in terms of RZ gates using basis change
        qp.Hadamard(wires=0)
        qp.RZ(2 * coeffs[0] * del_t, wires=0)
        qp.Hadamard(wires=0)

        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, coeffs, test_r)  # Naïve approach
Trotter_resources_PG = qre.estimate(TrotterStepperPG)(t, coeffs, test_r)  # Phase gradient approach

So, in the case of the naïve approach, the estimated requirements are:

print(Trotter_resources)

--- 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

In the case of the phase gradient method, the estimated requirements are:

print(Trotter_resources_PG)

--- Resources: ---
 Total wires: 20
   algorithmic wires: 16
   allocated wires: 4
     zero state: 4
     any state: 0
 Total gates : 1.840E+4
   'Toffoli': 1.600E+3,
   'CNOT': 1.160E+4,
   'Hadamard': 5.200E+3

As expected, the T count is much lower in the phase gradient case, even when adding in the Toffoli gate decomposition taking 1 Toffoli = 4 T gates [4]. Things are looking good, but let us compare the Trotter error in the two methods before we jump to any conclusions.

Recalling our naïve approach accuracy,

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, coeffs, r)
    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

Using the same comparison to the analytical solution, we can estimate the error in the phase gradient approach.

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, coeffs, r)
    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}"
    )

r |  X error |  Y error |  Z error | Total error
---------------------------------------------------
   10 |  0.25203 |  0.14019 |  0.48879 |     0.56753
   20 |  0.18798 |  0.14019 |  0.64247 |     0.68393
   30 |  0.03503 |  0.14019 |  0.43239 |     0.45590
   40 |  0.24055 |  0.14019 |  0.69749 |     0.75101
   50 |  0.09209 |  0.14019 |  0.01417 |     0.16833
  100 |  0.09209 |  0.14019 |  0.01417 |     0.16833
  200 |  0.09209 |  0.14019 |  0.01417 |     0.16833
  400 |  0.09209 |  0.14019 |  0.01417 |     0.16833

Ah ha! A trade-off has made itself clear! In the phase gradient implementation, the Trotter error is universally higher. What is also (maybe even more) interesting is that, after a certain \(r\) threshold is reached, the phase gradient system no longer evolves. This is an example of underflow in quantum arithmetic, where the simulation has reached a limit and is stuck rounding to the same value each pass. Since the phase gradient method relies on quantum addition, it may not be the right tool here if our main goal is to reduce error. 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.

Another thing to consider is the potential fast-forwardability of the system [3]. As shown, carrying out quantum simulation on hardware requires a series of gates to be implemented for each time step, meaning the depth of the circuit grows notably with increasing time. If the depth of a simulation circuit maintains proportionality to the length of the time interval being simulated, for example, it runs the (very real) risk of exceeding the coherence time of the system. Ideally, running a simulation for time \(t\) would require a complexity less than \(\mathcal{O}(t)\) or, in other words, a complexity that is sublinear in \(t\). This, in theory, can be achieved by employing a fast-forwarding technique, in which the phase angles used in the time evolution operator can be strategically altered to cover more time in a single step to achieve sublinear complexity. This technique requires the Hamiltonian have a known analytical solution (or, in other words, that it is “exactly integrable”), which is often not the case.

Conclusion

There is a phenomenon that is renewed over and over again in which the field of mathematics is continually years ahead of physics (especially applied physics). The use of product formulas as tools for time evolution in Hamiltonian simulation is a beautiful example of this, in which a purely mathematical description of how to exponentiate non-commuting operators has become a defining method for time-evolution simulation in today’s quantum pursuits. Understanding the basics of how Trotter products can be used and optimized for various applications of quantum simulation opens the door to not only increased utility but a heightened awareness of how the input of many fields is required to achieve viable outcomes. If you are, wisely, looking to add Trotterization to your quantum skill set, check out our Codebook chapter on product formulas. Keep calm and Trotter on!

References