Hamiltonian Simulation

View on GitHub

Open this notebook in GitHub to run it yourself

Quantum Hamiltonian simulation is one of the most important problems in quantum computing. It consists of modeling the time evolution of a complex quantum system using more controllable devices; the quantum computers. Simulations of this type are essential to comprehend quantum systems that are described by complex dynamics, such as molecules in chemistry or materials in condensed matter physics, where classical computers fail due to the exponential scale up required by larger quantum systems. In this guide, you will be using Classiq to work with simple problems of Hamiltonian simulation using two different methods:

Suzuki-Trotter decomposition
qDRIFT

Apart from these methods, known as product formulas, this guide also briefly discusses block-encoding methods and links to a specific guide on this topic.

Intoduction
Suzuki-Trotter decomposition
- Simple Exponentiation of Suzuki Trotter formula of order 1 with 1 repetition
qDRIFT
Hamiltonian simulation with block encoding
Measuring the expected magnetization as a function of time
Summary and Exercises

Introduction

Quantum Hamiltonian simulation has several significant applications. Some of them are directly related to physics and chemistry, while others extend to fields such as optimization and machine learning:

Quantum Chemistry: It helps predict the properties of molecules using quantum methods such as the Quantum Variational Eigensolver (VQE), which can be more efficient than classical methods in certain conditions.
Materials Science: It aids in designing new materials with desired properties by simulating their quantum mechanical behavior.
Optimization and Machine Learning: The use of quantum Hamiltonian simulation algorithms in optimization and machine learning is an active area of research.

These algorithms are often incorporated as subroutines within broader quantum algorithms, combined with other quantum techniques tailored to optimization and machine learning tasks. Within the framework of product formulas, we can mathematically approach the Hamiltonian Simulation problem as follows: Suppose we have a quantum computer that can run the evolution operators of the set of Hamiltonians

\{H_j\}

natively. Now, suppose you want to understand the dynamics of the following, more complex, Hamiltonian:

What are the operations that this quantum computer can effectively execute?

If this quantum computer can only execute the evolution operators of the set of Hamiltonians

\{H_j\}

, then you could only have operations in the form:

O = \prod_{j} e^{-i \,\tau_j H_j}.

\begin{equation} H_{full} = \sum_{j} h_j H_j, \end{equation}

For some real values

h_j

. Therefore, a decomposition of the operator into the set of operations is needed, formed by

\{e^{-itH_j}\}

, where

e^{-itH_j}

is the evolution operator of the Hamiltonian

H

Suzuki-Trotter Decomposition

Now that the problem is stated, a possible solution is the Suzuki-Trotter decomposition, also known as Trotterization. Simply put, this method involves ‘breaking’ the evolution operator into the evolution operators of its components. In the first order, it looks like this:

\begin{equation*} e^{-itH}=\exp\left\{-it\sum_{j=1}^N h_j H_j\right\} \approx \left(\prod_{j}^N e^{-it h_j H_j/r}\right)^r + \mathcal{O}(t^2/r). \end{equation*}

Details about this decomposition are in [1]. This formula is important and has applications in several scenarios. Additionally, if you need an error that scales better than

t^2

, it is possible to achieve higher-order Suzuki-Trotter formulas. For example, the second-order Suzuki-Trotter would look like this:

\begin{equation*} e^{-itH}=\exp\left\{-i\left(\sum_{j=1}^N \frac{t_j}{2}\,H_j + \sum_{j=N}^1 \frac{t_j}{2}\,H_j\right)\right\} \approx \left(\prod_{j=1}^N e^{-it_j\, H_j/(2r)}\prod_{j=N}^1 e^{-it_j\, H_j/(2r)}\right)^r + \mathcal{O}(t^3/r). \end{equation*}

Higher orders for the Suzuki-Trotter formula are in [1]. However, as the number of terms in the Hamiltonian and the approximation order increases, it gets more complicated to construct the Suzuki-Trotter formulas. To this end, use the Classiq suzuki_trotter() function to execute an operator’s Trotter decomposition, specifying the order and number of repetitions. The function’s inputs are:

pauli_operator: SparsePauliOp: The Pauli operator to be exponentiated, representing the term $H_j$ and its respective coefficients $h_j$ ;
evolution_coefficient: CArray: A global evolution coefficient multiplying the Pauli operator, representing the value $t$ for the coefficient in the time evolution operator;
order: CInt: The order of the Suzuki-Trotter decomposition;
repetitions: CInt: The number of repetitions of the Suzuki-Trotter decomposition;
qbv: QArray: The target quantum variable of the exponentiation, representing the qubits on which the operation will be applied.

The following example demonstrates how to use it: Now that the inputs for the Suzuki-Trotter method are defined, apply it to the example of approximating the dynamics of the following two-qubits Hamiltonian:

H = 0.3 \, Z\otimes Z + 0.7 \,X\otimes I + 0.2\, I\otimes X,

which represents two spins interacting in a transversal field. Identify the operators

\{H_j\}

and their respective coefficients

\{h_j\}

\begin{split} H_1 = Z \otimes Z,\text{ coefficient: }h_1 = 0.3; \\ H_2 = X \otimes I,\text{ coefficient: }h_2 = 0.7; \\ H_3 = I \otimes X,\text{ coefficient: }h_3 = 0.2. \\ \end{split}

Define the number of

r

repetitions into which to “break” the evolution, in this case,

r=10

. The order of the Suzuki-Trotter decomposition can be altered in the parameter order. In the current case, apply it to the first order. Define the evolution coefficient

t

, stating how much time the quantum system will be evolving under this Hamiltonian, according to the specific application you want. Set the evolution coefficient in the evolution_coefficient variable, considering it unitary in this example. In this part of the guide, you can follow an example where the evolution coefficient varies.

from classiq import *

# Defining the Hamiltonian for product formulas:
HAMILTONIAN = 0.3 * Pauli.Z(0) * Pauli.Z(1) + 0.7 * Pauli.X(0) + 0.2 * Pauli.X(1)


@qfunc
def main(qba: Output[QArray]):
    allocate(2, qba)
    suzuki_trotter(
        HAMILTONIAN,
        evolution_coefficient=1.0,
        order=1,
        repetitions=10,
        qbv=qba,
    )


qprog_trotter = synthesize(main)
show(qprog_trotter)

Output:

Quantum program link: https://platform.classiq.io/circuit/3DtdllWxDeJxISWzBmUSbC06a9G
  

Exponentiate

It is also possible to avoid choosing the Suzuki Trotter parameters, and simply use the default ones of order=1 and repetitions=1. For this we can call Classiq exponentiate function.

@qfunc
def main(qba: Output[QArray]):
    allocate(2, qba)
    exponentiate(
        HAMILTONIAN,
        evolution_coefficient=1,
        qbv=qba,
    )


qprog_exponentiation = synthesize(main)
show(qprog_exponentiation)

Output:

Quantum program link: https://platform.classiq.io/circuit/3Dtdm975OAGJk01tNb5xQL3MbKm
  

The number of repetitions can also be understood from the IDE, it is the number of repetitive blocks you see. In the current case, it applies the first-order decomposition, denoted in the IDE by the name single_trotter_suzuki_layer_xxxx. In the case of higher orders, the output would identify which order is being used. Trotter_gif

qDRIFT

The quantum stochastic drift protocol (qDRIFT) [2] is similar to the first-order Suzuki-Trotter decomposition; however, it relies on a stochastic distribution of the evolution operators. The algorithm operates by sampling unitaries from the set

\{e^{-itH_j}\}

according to a probability distribution defined by the weights

\{h_j\}

, normalized by a factor

\lambda = \sum_j h_j

. In a more structured manner, the qDRIFT protocol functions as follows: Input: A list of Hamiltonian terms

\{H_j\}

, a classical oracle function SAMPLE() that returns a value

j

according to the probability distribution

p_j = h_j/\lambda

, and a target precision

\epsilon

. Output: An ordered list of evolution operators from the set

\{e^{-itH_j}\}

that approximates the unitary

e^{-it\,H}

with a bound error

\epsilon

Define a normalization to transform $\{h_j\}$ into a probability distribution.

For this, define

\sum_j h_j = \lambda

The program depth is defined according to the required precision.

For this, set

N = \lceil 2\lambda^2t^2/\epsilon \rceil

Now generate the ordered list of $N$ evolution operators according to the probability distribution.

The qDRIFT has proved to be a good alternative for the Suzuki-Trotter decomposition when the number of terms in the Hamiltonian is not Pauli sparse, i.e., the number of terms on the expansion

H = \sum_j h_j H_j

is not small when compared to all possible terms in it. This interesting method can be applied using the Classiq qdrift() function. Its inputs are:

pauli_operator: CArray[PauliTerm]: A list of Pauli operators that represent the term $H_j$ and its respective coefficients $h_j$ ;
evolution_coefficient : CArray The value $t$ for the coefficient in the time evolution operator;
num_qdrift:CInt: The number $N$ of unitary operators in the list of gates given by the qDRIFT;
qbv:QArray: The target qubits.

In the current example, fix

N=288

, which corresponds to error

\epsilon = 0.01

in the approximation.

@qfunc
def main(qba: Output[QArray]):
    allocate(2, qba)
    qdrift(
        HAMILTONIAN,
        evolution_coefficient=1,
        num_qdrift=288,
        qbv=qba,
    )


qprog_qdrift = synthesize(main)
show(qprog_qdrift)

Output:

Quantum program link: https://platform.classiq.io/circuit/3DtdmurE56VJL70WHs0ycPSPMMu
  

Hamiltonian Simulation with Block Encoding

Other methods of Hamiltonian simulation outside of the product formulas leverage the idea of Block Encoding and eigenvalue transformations. These approaches allow for the efficient simulation of quantum systems by encoding operators into unitary matrices and applying transformations to the eigenvalues of these matrices.

Qubitization

The qubitization [3] method enables the efficient simulation of quantum Hamiltonians. The key idea is to embed the Hamiltonian into a larger unitary operator, a method known as block encoding, where the block-encoded Hamiltonian is simulated using a quantum walk operator. This allows precise control over the evolution of quantum states and can significantly reduce the resources required for Hamiltonian simulation. For details on implementing the Qubitization method, visit the Hamiltonian simulation with qubitization.

QSVT

The general Quantum Singular Value Transformation (QSVT) [4] method transforms the singular values of matrices encoded in quantum states. In Hamiltonian simulation, we need to combine two QSVT blocks, one for a sine operation and another for the cosine operation, for constructing the exponential of the matrix. For details on implementing the QSVT method, visit the Hamiltonian simulation with QSVT.

GQSP

Generalized Quantum Signal Processing (GQSP) [5] achieves Hamiltonian simulation by applying a Laurent polynomial to the Szegedy walk operator of the block-encoded Hamiltonian. The polynomial coefficients are derived from the Jacobi–Anger expansion and approximate the time-evolution phase

e^{-i\lambda t}

on each eigenvalue of the Hamiltonian, requiring only one auxiliary qubit and no amplitude amplification. For details on implementing the GQSP method, visit the Hamiltonian simulation with GQSP.

Comparison of Block-Encoding Methods

Method	Extra block qubits	Controlled $U_H$ ?	Amplitude amplification?	Classical preprocessing
GQSP	1	Yes	No	Angle computation
QSVT	2	No	Yes (for a factor of 2)	Angle computation
Qubitization	$O(\log d)$	Yes	Yes (for the sum of Cheb. coefficients)	None

All three methods share the same asymptotic query complexity.

Measuring the Expected Magnetization as a Function of Time

Learn how to apply these different methods and evaluate their performance in comparison with the exact evolution. You will see the evolution of the expected value of a specific operator and compare it to the exact solution. This section uses the same Hamiltonian from the examples:

H = 0.3 \cdot \, Z\otimes Z + 0.7 \cdot \,X\otimes I + 0.2 \cdot \, I\otimes X

Now, analyze a quantity that varies with different values of time. For this, consider the behavior of the magnetization of the system,

M = (\langle I\otimes Z \rangle + \langle Z \otimes I\rangle)/2

[6], as a function of time. Precisely evaluate the result by exponentiating the Hamiltonian

H

directly, and then comparing the results:

import numpy as np
from scipy.linalg import expm

time_list = np.linspace(0, 2, 100).tolist()

magnetization_hamiltonian = [
    PauliTerm(pauli=[Pauli.Z, Pauli.I], coefficient=0.5),
    PauliTerm(pauli=[Pauli.I, Pauli.Z], coefficient=0.5),
]

# Magnetization observable for observe:
magnetization_HAMILTONIAN = 0.5 * Pauli.Z(0) + 0.5 * Pauli.Z(1)

# Hamiltonians to matrices:
magnetization_matrix = hamiltonian_to_matrix(magnetization_hamiltonian)
Hamiltonian_matrix = hamiltonian_to_matrix(HAMILTONIAN)

initial_state = np.zeros(4)
initial_state[0] = 1.0


def expected_value(state, operator):
    state_H = np.conj(state.T)
    return state_H @ operator @ state


ideal_magnetization = []
for t in time_list:
    state = expm(-1j * t * Hamiltonian_matrix) @ initial_state
    ideal_magnetization.append(expected_value(state, magnetization_matrix))

To evaluate the expectation values of magnetization for all three methods, use the observe function:

@qfunc
def main(t: CReal, qba: Output[QArray]):
    allocate(2, qba)
    suzuki_trotter(
        HAMILTONIAN,
        evolution_coefficient=t,
        order=1,
        repetitions=30,
        qbv=qba,
    )


qprog_magnetization_trotter = synthesize(main)
t_values = [{"t": times} for times in time_list]

magnetization_ST_results = observe(
    qprog_magnetization_trotter, magnetization_HAMILTONIAN, parameters=t_values
)

Output:

Submitting job to simulator
  Job: https://platform.classiq.io/jobs/04251f67-6b2c-4e3f-b4ae-38ff7e7e0c1d
  

@qfunc
def main(t: CReal, qba: Output[QArray]):
    allocate(2, qba)
    qdrift(
        HAMILTONIAN,
        evolution_coefficient=t,
        num_qdrift=288,
        qbv=qba,
    )


qprog_magnetization_qdrift = synthesize(main)

magnetization_qdrift_results = observe(
    qprog_magnetization_qdrift, magnetization_HAMILTONIAN, parameters=t_values
)

Output:

Submitting job to simulator
  Job: https://platform.classiq.io/jobs/33a215be-ff75-4d3f-a3b3-2487dff359ff
  

@qfunc
def main(t: CReal, qba: Output[QArray]):
    allocate(2, qba)
    exponentiate(
        HAMILTONIAN,
        evolution_coefficient=t,
        qbv=qba,
    )


qprog_magnetization_exponentiate = synthesize(main)

magnetization_exponentiate_results = observe(
    qprog_magnetization_exponentiate, magnetization_HAMILTONIAN, parameters=t_values
)

Output:

Submitting job to simulator
  Job: https://platform.classiq.io/jobs/1b73b90f-9638-4013-af04-a768535dbd47
  

You can plot the magnetization values of the three different methods with the exact values:

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.set_xlabel("Evolution coefficient")
ax.set_ylabel("Magnetization")

magnetization_ST = magnetization_ST_results
magnetization_qdrift = magnetization_qdrift_results
magnetization_exponentiate = magnetization_exponentiate_results

plt.plot(
    time_list,
    magnetization_ST,
    label="Suzuki-Trotter 1st order",
)
plt.plot(time_list, magnetization_exponentiate, label="Exponentiate")
plt.plot(time_list, magnetization_qdrift, label="qDRIFT")
plt.plot(time_list, np.real(ideal_magnetization), label="Ideal Simulation")
plt.legend()
plt.show()

It’s possible to observe that the larger t is, the worse the approximation gets for the same set of arguments. To get a better approximation, there is a need for a deeper quantum program with more layers and or with a higher order of approximation.

Summary and Exercises

This guide introduced the Suzuki-Trotter and qDRIFT methods for Hamiltonian simulation, highlighting the simplicity of Classiq’s high level functional design in tackling complex problems. To become more familiar with the methods and see how the synthesis engine is able to optimize quantum models for Hamiltonian simulation, try to apply the qDrift and Suzuki-Trotter methods to this tutorial as well as different values of the evolution_coefficient:

Hamiltonian Evolution for a Water Molecule

Another good exercise is solving molecules other than the water molecule. With Classiq, you can generate the Hamiltonian of any valid molecule structure, following the same process for the water molecule.

Besides the methods presented in this guide, there are also mixed approaches. One of them, proposed by Matthew Hagan and Nathan Wiebe [7], offers an interesting combination of the Suzuki-Trotter and qDRIFT methods.

References

[1]: Finding exponential product formulas of higher orders (Naomichi Hatano and Masuo Suzuki) [2]: A random compiler for fast Hamiltonian simulation (Earl Campbell) [3]: Hamiltonian simulation by qubitization (Guang Hao Low and Isaac L. Chuang) [4]: Quantum singular value transformation and beyond: Exponential improvements for quantum matrix arithmetics (András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe) [5]: Motlagh, D., and Wiebe, N. Generalized quantum signal processing. PRX Quantum 5, 020368 (2024). [6]: The microscopic magnetization: Concept and application (L. L. Hirst) [7]: Composite quantum simulations (Matthew Hagan and Nathan Wiebe)

Algorithms

Applications

Functions

Tutorials

View on GitHub

Table of Contents

Introduction

Suzuki-Trotter Decomposition

Exponentiate

qDRIFT

Hamiltonian Simulation with Block Encoding

Qubitization

QSVT

GQSP

Comparison of Block-Encoding Methods

Measuring the Expected Magnetization as a Function of Time

Summary and Exercises

Read More

References

Algorithms

Applications

Functions

Tutorials

Documentation Index

View on GitHub

​Table of Contents

​Introduction

​Suzuki-Trotter Decomposition

​

​Exponentiate

​qDRIFT

​Hamiltonian Simulation with Block Encoding

​Qubitization

​QSVT

​GQSP

​

​Comparison of Block-Encoding Methods

​Measuring the Expected Magnetization as a Function of Time

​Summary and Exercises

​Read More

​References

Table of Contents

Introduction

Suzuki-Trotter Decomposition

Exponentiate

qDRIFT

Hamiltonian Simulation with Block Encoding

Qubitization

QSVT

GQSP

Comparison of Block-Encoding Methods

Measuring the Expected Magnetization as a Function of Time

Summary and Exercises

Read More

References