Hadamard Test

Consider the observable \(O\) and two generic quantum states \(\ket{\psi}\) and \(\ket{\phi}\). We want to calculate the quantity

\[\bra{\psi} O \ket{\phi}.\]

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 + Re \bra{\psi} O \ket{\phi} \right]\]

\[P(1) = \frac{1}{2} \left[ I - Re \bra{\psi} O \ket{\phi} \right]\]

The difference between the probability of \(0\) and \(1\) gives

\[P(0)-P(1) = Re \bra{\psi} O \ket{\phi}\]

A- Numerical result as a reference:

import cudaq
import numpy as np
from functools import reduce

cudaq.set_target('nvidia')

qubit_num = 2

@cudaq.kernel
def psi(num:int):
    q = cudaq.qvector(num)
    h(q[1])

@cudaq.kernel
def phi(n:int):
    q = cudaq.qvector(n)
    x(q[0])

psi_state = cudaq.get_state(psi, qubit_num)
print('Psi state: ', psi_state)

phi_state=cudaq.get_state(phi, qubit_num)
print('Phi state: ', phi_state)

ham=cudaq.spin.x(0) * cudaq.spin.x(1)
ham_matrix = ham.to_matrix()
print('hamiltonian: ', np.array(ham_matrix), '\n')

exp_val=reduce(np.dot,(np.array(psi_state).conj().T, ham_matrix, phi_state))

print('Numerical expectation value: ', exp_val)

Psi state:  (0.707107,0)
(0,0)
(0.707107,0)
(0,0)

Phi state:  (0,0)
(1,0)
(0,0)
(0,0)

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)

B- Using sample algorithmic primitive to sample the ancilla qubit and compute the expectation value.

import cudaq

cudaq.set_target('nvidia')

@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 ham_cir(q:cudaq.qview):
    x(q[0])
    x(q[1])

@cudaq.kernel
def kernel(n:int):
    ancilla=cudaq.qubit()
    q = cudaq.qvector(n)
    h(ancilla)
    cudaq.control(U_phi,ancilla,q)
    cudaq.control(ham_cir,ancilla,q)
    cudaq.control(U_psi,ancilla,q)

    h(ancilla)

    mz(ancilla)

shots = 100000
qubit_num=2
count = cudaq.sample(kernel, qubit_num, shots_count = shots)
print(count)

mean_val = (count['0']-count['1']) / shots
error = np.sqrt(2* count['0'] * count['1'] / shots) / shots
print('Observable QC: ', mean_val,'+ -', error)
print('Numerical result', np.real(exp_val))

{ 0:85356 1:14644 }

Observable QC:  0.70712 + - 0.0015811092713661505
Numerical result 0.7071067690849304

C- Use multi-GPUs to compute the matrix elements

import cudaq

cudaq.set_target("nvidia-mqpu")

target = cudaq.get_target()
qpu_count = target.num_qpus()
print("Number of QPUs:", qpu_count)

@cudaq.kernel
def U_psi(q:cudaq.qview, theta:float):
    ry(theta, q[1])

@cudaq.kernel
def U_phi(q:cudaq.qview, theta: float):
    rx(theta, q[0])

@cudaq.kernel
def ham_cir(q:cudaq.qview):
    x(q[0])
    x(q[1])

@cudaq.kernel
def kernel(n:int, angle:float, theta:float):
    ancilla = cudaq.qubit()
    q = cudaq.qvector(n)
    h(ancilla)
    cudaq.control(U_phi, ancilla, q, theta)
    cudaq.control(ham_cir, ancilla, q)
    cudaq.control(U_psi, ancilla, q, angle)

    h(ancilla)

    mz(ancilla)

shots = 100000
angle = [0.0, 1.5,3.14,0.7]
theta = [0.6, 1.2 ,2.2 ,3.0]
qubit_num = 2

result = []
for i in range(4):
    count = cudaq.sample_async(kernel, qubit_num, angle[i], theta[i], shots_count = shots, qpu_id = i%qpu_count)
    result.append(count)

mean_val = np.zeros(len(angle))
i = 0
for count in result:
    print(i)
    i_result = count.get()
    print(i_result)
    mean_val[i] = (i_result['0'] - i_result['1']) / shots
    error = np.sqrt(2 * i_result['0'] * i_result['1'] / shots) / shots
    print('Observable QC: ',  mean_val[i],'+ -', error)
    i += 1

Number of QPUs: 5
0
{ 0:63807 1:36193 }

Observable QC:  0.27614 + - 0.0021491238917289066
1
{ 0:49929 1:50071 }

Observable QC:  -0.00142 + - 0.0022360657230949183
2
{ 0:50041 1:49959 }

Observable QC:  0.00082 + - 0.0022360672257336093
3
{ 0:50276 1:49724 }

Observable QC:  0.00552 + - 0.0022360339102974265

Diagonalize the matrix using for example Numpy or CuPy. In this example, since we are having 2x2 matrix, we use numpy.

import numpy as np

my_mat = np.zeros((2,2),dtype=float)
m = 0
for k in range(2):
    for j in range(2):
        my_mat[k,j] = mean_val[m]
        m += 1

print(my_mat)

E,V = np.linalg.eigh(my_mat)

print('Eigen values: ')
print(E)

print('Eigenvector: ')
print(V)

[[ 0.27614 -0.00142]
 [ 0.00082  0.00552]]
Eigen values:
[0.00551752 0.27614248]
Eigenvector:
[[ 0.00303004 -0.99999541]
 [-0.99999541 -0.00303004]]