Hadamard Test and Application¶
Consider the observable \(O\) and two generic quantum states \(\ket{\psi}\) and \(\ket{\phi}\). We want to calculate the quantity
\[\braket{\psi | O | \psi}.\]
where \(O\) is a Pauli operator.
First of all we shall prepare the states \(\ket{\psi}\) and \(\ket{\phi}\) using a quantum circuit for each of them. So we have
\[\ket{\psi} = U_{\psi}\ket{0} \qquad \ket{\phi} = U_{\phi}\ket{0}\]
Let’s define an observable we want to use:
\[O = X_1X_2\]
Now we can evaluate the matrix element using the following fact:
\[\bra{\psi}O\ket{\phi} = \bra{0}U_\psi^\dagger O U_\phi\ket{0}\]
This is just an expectation value which can be solved with a simple Hadamard test. The probability to measure \(0\) or \(1\) in the ancilla qubit is
\[P(0) = \frac{1}{2} \left[ I + \operatorname{Re} \bra{\psi} O \ket{\phi} \right]\]
\[P(1) = \frac{1}{2} \left[ I - \operatorname{Re} \bra{\psi} O \ket{\phi} \right]\]
The difference between the probability of \(0\) and \(1\) gives
\[\braket{X} = P(0)-P(1) = \operatorname{Re} \braket{\psi | O |
\phi}.\]
Similarly, the imaginary part can be obtained from Y measurement
\[\braket{Y} = \operatorname{Im} \braket{\psi | O | \phi}.\]
Combining these results, the quantity \(\braket{\psi | O | \psi}\) is obtained.
Numerical result as a reference:¶
[1]:
from functools import reduce
import numpy as np
import cudaq
cudaq.set_target("nvidia")
num_qubits = 2
@cudaq.kernel
def psi(num_qubits: int):
q = cudaq.qvector(num_qubits)
h(q[1])
@cudaq.kernel
def phi(num_qubits: int):
q = cudaq.qvector(num_qubits)
x(q[0])
psi_state = cudaq.get_state(psi, num_qubits)
print("Psi state: ", np.array(psi_state))
print(cudaq.draw(psi, 2))
phi_state = cudaq.get_state(phi, num_qubits)
print("Phi state: ", np.array(phi_state))
print(cudaq.draw(phi, 2))
ham = cudaq.spin.x(0) * cudaq.spin.x(1)
ham_matrix = ham.to_matrix()
print("hamiltonian: ", np.array(ham_matrix), "\n")
ev_numerical = np.array(psi_state).conj() @ ham_matrix @ np.array(phi_state).T
print("Numerical expectation value: ", ev_numerical)
Psi state: [0.70710677+0.j 0. +0.j 0.70710677+0.j 0. +0.j]
q0 : ─────
╭───╮
q1 : ┤ h ├
╰───╯
Phi state: [0.+0.j 1.+0.j 0.+0.j 0.+0.j]
╭───╮
q0 : ┤ x ├
╰───╯
hamiltonian: [[0.+0.j 0.+0.j 0.+0.j 1.+0.j]
[0.+0.j 0.+0.j 1.+0.j 0.+0.j]
[0.+0.j 1.+0.j 0.+0.j 0.+0.j]
[1.+0.j 0.+0.j 0.+0.j 0.+0.j]]
Numerical expectation value: (0.7071067690849304+0j)
Using observe
algorithmic primitive to compute the expectation value for ancilla qubits.¶
[2]:
@cudaq.kernel
def u_psi(q: cudaq.qview):
h(q[1])
@cudaq.kernel
def u_phi(q: cudaq.qview):
x(q[0])
@cudaq.kernel
def apply_pauli(q: cudaq.qview):
x(q[0])
x(q[1])
@cudaq.kernel
def kernel(num_qubits: int):
ancilla = cudaq.qubit()
q = cudaq.qvector(num_qubits)
h(ancilla)
cudaq.control(u_phi, ancilla, q)
cudaq.control(apply_pauli, ancilla, q)
cudaq.control(u_psi, ancilla, q)
num_qubits = 2
shots = 100000
x_0 = cudaq.spin.x(0)
y_0 = cudaq.spin.y(0)
results = cudaq.observe(kernel, [x_0, y_0], num_qubits, shots_count=shots)
evs = np.array([result.expectation() for result in results])
std_errs = np.sqrt((1 - evs**2) / shots)
print(f"QC result: {evs[0]:.3f}+{evs[1]:.3f}i ± {std_errs[0]:.3f}+{std_errs[1]:.3f}i")
print("Numerical result", ev_numerical)
QC result: 0.705+-0.007i ± 0.002+0.003i
Numerical result (0.7071067690849304+0j)
Use multi-GPUs to compute multiple Hadamard test in parallel¶
[3]:
# Use multi-QPUs
cudaq.set_target("nvidia-mqpu")
target = cudaq.get_target()
num_qpus = target.num_qpus()
print("Number of QPUs:", num_qpus)
@cudaq.kernel
def u_psi(q: cudaq.qview, theta: float):
ry(theta, q[1])
@cudaq.kernel
def u_phi(q: cudaq.qview, theta: float):
s(q[0])
rx(theta, q[0])
s(q[0])
@cudaq.kernel
def ham_circuit(q: cudaq.qview):
x(q[0])
x(q[1])
@cudaq.kernel
def kernel(angle0: float, angle1: float):
ancilla = cudaq.qubit()
q = cudaq.qvector(2)
h(ancilla)
cudaq.control(u_phi, ancilla, q, angle0)
cudaq.control(ham_circuit, ancilla, q)
cudaq.control(u_psi, ancilla, q, angle1)
angles = 2 * np.pi * np.random.rand(10, 2)
print(f"{angles=}")
async_results = [
cudaq.observe_async(
kernel,
x_0,
float(angle[0]),
float(angle[1]),
qpu_id=i % num_qpus,
)
for i, angle in enumerate(angles)
]
for i, async_result in enumerate(async_results):
ev = async_result.get().expectation()
print(f"{i}-th {ev=}")
Number of QPUs: 5
angles=array([[1.56322878, 3.09176639],
[0.04025496, 5.59986135],
[1.87024074, 0.93078226],
[4.44015281, 5.05675948],
[1.92402471, 2.12981374],
[0.49704605, 3.6020906 ],
[4.50280746, 2.78988978],
[4.006956 , 3.7581442 ],
[3.00524035, 3.10937881],
[3.13405202, 1.33235091]])
0-th ev=-0.7042075991630554
1-th ev=-0.006743329111486673
2-th ev=-0.36111390590667725
3-th ev=-0.45839524269104004
4-th ev=-0.7175908088684082
5-th ev=-0.23948131501674652
6-th ev=-0.765204668045044
7-th ev=-0.865047037601471
8-th ev=-0.9975475072860718
9-th ev=-0.6179792881011963