Hamiltonian simulation with quantum computation
One of the most sought-after applications of quantum computing is the simulation of physical processes, e.g., simulating molecules for drug development and many-body problems. These simulations generally involve taking a Hamiltonian operator H and a time t which is used to evolve an arbitrary state from |ψ⟩ to e^[iHt]|ψ⟩. Creating a circuit that applies the operator e^[iHt] is not immediately obvious, but we can start off with a two-qubit example.
Simulating σ_i exponentials
Suppose we want to simulate e^[i σ_z⊗σ_z t]. The circuit to do this is known to be [1]:
where the argument to the rotation gate R_z is −2t.
Showing that this circuit is equivalent to e^[i σ_z⊗σ_z t] is not too hard. Remember that we have the Euler-like decomposition of the operator as e^[i σ_z⊗σ_z t]=cos(t)I+i sin(t) σ_z⊗σ_z (which works since we have that σ_z^2=I). Thus, we have to show that the circuit above simplifies to this same expression. Doing this is a bit long but not hard, as shown below:
This can easily be extended to any number of qubits, i.e., e^[i σ_z⊗σ_z⊗⋯⊗σ_z t] [3]. Now, what happens if we have a Pauli matrix that is not σ_z in the exponent? Well, just apply a self-inverse gate U before and after the chain of CXs such that σ_k=Uσ_zU, where k∈{x,y}. This gate would be the Hadamard gate for σ_x and
for σ_y. And if its the identity, you can just leave that qubit unaffected!
Simulating general exponentials
But this just tells us how to simulate Pauli strings, and we may not be given the Hamiltonian H on this form. However, it is a fact that for any matrix A with its dimension being a power of 2 (which would be the case here since we are working with qubits), there exists some decomposition into a Pauli string of the following form [2]:
where
Thus, knowing how to simulate Pauli strings, we can basically simulate any Hamiltonian H. Using the Pauli string of a Hamiltonian, we can construct the Trotter decomposition for which we can create a circuit based of the simple circuit we saw above. The Trotter decomposition looks as follows:
where P_k is a tensor product of Pauli matrices, c_k the coefficient corresponding to each P_k, t the time introduced at the beginning, and N a number called the Trotter number. The bigger N is, the more accurate our simulation is. The specifics of Trotter decomposition are out of the scope of this article, but a good start is chapter 4.1 of [3].
Qiskit implementation
With this background in place, we can start implementing the simulation in Qiskit. First we are going to define a simple function that creates the circuit scheme introduced at the beginsning for operators of the form e^[i σ_z⊗σ_z⊗⋯⊗σ_z t].
from qiskit import QuantumCircuit
import numpy as np
def sim_z(t, qc, qubits):
"""
Add gates to simulate a Pauli Z string exponential
Parameters:
t: float
Time parameter to simulate for
qc: QuantumCircuit
Circuit to append gates to
qubits: array
Array indicating qubits indeces (in order) to
append the gates to
"""
for i in range(len(qubits) - 1):
qc.cx(qubits[i], qubits[i + 1])
qc.rz(-2 * t, qubits[-1])
for i in range(len(qubits) - 1, 0, -1):
qc.cx(qubits[i - 1], qubits[i])qc = QuantumCircuit(3)
sim_z(np.pi, qc, [0, 1, 2])
qc.draw()
Which outputs:
q_0: ──■─────────────────────────■──
┌─┴─┐ ┌─┴─┐
q_1: ┤ X ├──■───────────────■──┤ X ├
└───┘┌─┴─┐┌─────────┐┌─┴─┐└───┘
q_2: ─────┤ X ├┤ Rz(-2π) ├┤ X ├─────
└───┘└─────────┘└───┘Now, we can expand this so we can simulate any Pauli string.
def sim_pauli(arr, t, qc, qubits):
"""
Append gates to simulate any Pauli string
Parameters:
arr: array
Array encoding the Pauli string
t: float
Time parameter to simulate for
qc: QuantumCircuit
Circuit to append gates to
qubits: array
Array indicating qubits indeces (in order) to
append the gates to
"""
new_arr = []
new_qub = []
for idx in range(len(arr)):
if arr[idx] != 'I':
new_arr.append(arr[idx])
new_qub.append(qubits[idx])
h_y = 1 / np.sqrt(2) * np.array([[1, -1j], [1j, -1]])
for i in range(len(new_arr)):
if new_arr[i] == 'X':
qc.h(new_qub[i])
elif new_arr[i] == 'Y':
qc.unitary(h_y, [new_qub[i]], r'$H_y$')
sim_z(t, qc, new_qub)
for i in range(len(new_arr)):
if new_arr[i] == 'X':
qc.h(new_qub[i])
elif new_arr[i] == 'Y':
qc.unitary(h_y, [new_qub[i]], r'$H_y$')And finally, we write a small function that takes in a Hamiltonian H as its Trotter decomposition and simulates it in a given quantum circuit.
def sim_ham(hamiltonian, t, qc, qubits, trotter=1):
"""
Simulates Hamiltonian given as Pauli string
Parameters:
hamiltonian: dict
Dictionary encoding the hamiltonian with each
Pauli product as a key with the coefficient as value
t: float
Time parameter to simulate for
qc: QuantumCircuit
Circuit to append gates to
qubits: array
Array indicating qubits indeces (in order) to
append the gates to
"""
temp = QuantumCircuit(len(qubits))
delta_t = t / trotter
for pauli in hamiltonian:
sim_pauli(pauli, hamiltonian[pauli] * delta_t, temp, range(len(qubits)))
for i in range(trotter):
qc.compose(temp, qubits, inplace=True)Analyzing this function, we can see that it basically resembles this formula that we introduced in the background:
First, we see that delta_t corresponds to t/N. Then, for each Pauli term in the Pauli string representing the Hamiltonian (the first loop corresponds to the product ∏_k, since consecutive gates correspond to matrix multiplication), we call our function sim_pauli that adds the corresponding gates, using time = hamiltonian[pauli] * delta_t, which corresponds to c_k t/N in the formula. Finally, we simply append the resulting circuit N times in the second loop, which resembles the exponentiation.
Comparing to exact calculation
To wrap everything up, we can see how our simulation circuit compares to directly computing e^[iHt]|ψ⟩. Here, we test this with the Hamiltonian H=2XZY+5ZXX+2YXZ, time 1/(2∗π) and N=50, which gives us high accuracy.
from qiskit import Aer, execute
from sympy import Matrix
qc = QuantumCircuit(3)
sim_ham({"XZY": 2, "ZXX": 5, "YXZ": 2}, 1 / (2 * np.pi), qc, [0, 1, 2], trotter=50)
qc = qc.reverse_bits()
backend = Aer.get_backend('statevector_simulator')
result = execute(qc, backend).result()
vec = result.get_statevector()
vec = vec / vec[0] # global phase
qvec = vec / np.linalg.norm(vec) # normalize
Matrix(np.round(qvec, 5))
Now that we have our result with a quantum circuit, we can directly compute the expected result and compare.
import scipy
# Start with |0> state
start = np.zeros(2 ** 3)
start[0] = 1
# Get the matrix corresponding to some Pauli product
# This function can be optimized, but it is left as it
# is for clarity purposes
def get_pauli(string):
init = string[0]
string = string[1:]
if init == 'X':
out = np.array([[0, 1], [1, 0]])
elif init == 'Z':
out = np.array([[1, 0], [0, -1]])
elif init == 'Y':
out = np.array([[0, -1j], [1j, 0]])
else:
out = np.eye(2)
for p in string:
if p == 'X':
out = np.kron(out, np.array([[0, 1], [1, 0]]))
elif p == 'Z':
out = np.kron(out, np.array([[1, 0], [0, -1]]))
elif p == 'Y':
out = np.kron(out, np.array([[0, -1j], [1j, 0]]))
else:
out = np.kron(out, np.eye(2))
return out# Hamiltonian is calculated from the Pauli decomposition
decomp = {"XZY": 2, "ZXX": 5, "YXZ": 2}
H = np.zeros((2 ** 3, 2 ** 3))
H = H.astype('complex128')
for pauli in decomp:
H += decomp[pauli] * get_pauli(pauli)
# Hamiltonian is exponentiated and we multiply the starting
# vector by it to get our result
simul = scipy.linalg.expm(1j * H * (1 / (2 * np.pi)))
vec = simul @ start
vec = vec / vec[0] # global phase
cvec = vec / np.linalg.norm(vec) # normalize
Matrix(np.round(cvec, 5))
By simply comparing the two vectors, we can say that our simulation was pretty accurate! As a more mathematically rigorous comparison, we can use state fidelity, which also shows us that our simulation is pretty good!
from qiskit.quantum_info import state_fidelity
state_fidelity(qvec, cvec)Which gives:
0.9999834345555809There is a Jupyter notebook version of this article available here. If you enjoyed this post, remember to follow Quantum Untangled for more similar content!
References
- [1] An approach for Hamiltonian simulation by Davit Khachatryan on Quantum Computing StackExchange.
- [2] Can arbitrary matrices be decomposed using the Pauli basis? by Danylo Y on Quantum Computing StackExchange.
- [3] Simulation of electronic structure Hamiltonians using quantum computers by James D. Whitfield ,Jacob Biamonte and Alán Aspuru-Guzik.
