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