Hardware-efficient trial states for variational quantum eigensolvers

Cahyati Supriyati Sangaji (My Note)

h2_hamiltonian = (-1.0523732)  II + 
                 (0.39793742)  IZ + 
                 (-0.3979374)  ZI + 
                 (-0.0112801)  ZZ + 
                 (0.18093119)  XX

Note that these co-efficients are functions of the interatomic distance between the atoms. A term like IZ is shorthand notation for a tensor product of two-qubit operators -- the identity operator (𝐼) on one qubit and pauli-Z operator (𝑍) on the other qubit.

Installing necessary packages

Before we begin, you will need to install some prerequisites into your environment. Run the cell below to complete these installations. At the end, the cell outputs will be cleared.

!pip install -U qiskit==0.19
!pip install -U qiskit-ibmq-provider==0.7
from IPython.display import clear_output
clear_output()

Hardware efficient trial states

from numpy import pi
from qiskit import QuantumCircuit, Aer, execute
from qiskit.visualization import plot_histogram
def prepare_hets_circuit(depth, angle1, angle2):
hets_circ = QuantumCircuit(depth)
    hets_circ.ry(angle1, 0)
    hets_circ.rz(angle1, 0)
    hets_circ.ry(angle1, 1)
    hets_circ.rz(angle1, 1)
for ii in range(depth):
        hets_circ.cx(0,1)
        hets_circ.ry(angle2,0)
        hets_circ.rz(angle2,0)
        hets_circ.ry(angle2,1)
        hets_circ.rz(angle2,1)
    
    return hets_circ
hets_circuit = prepare_hets_circuit(2, pi/2, pi/2)
hets_circuit.draw()

Measuring expectation values

Next, we measure expectation values. We will begin by measuring the ZZ expectation value, or ⟨𝑍𝑍⟩. We will first create a copy of the hets_circ quantum circuit that we created above, and add measurements to it.

def measure_zz_circuit(given_circuit):
    zz_meas = given_circuit.copy()
    zz_meas.measure_all()
    return zz_meas
zz_meas = measure_zz_circuit(hets_circuit)
zz_meas.draw()

Next, let’s execute this quantum circuit and see the measurement outcomes.

simulator = Aer.get_backend('qasm_simulator')
result = execute(zz_meas, backend = simulator, shots=10000).result()
counts = result.get_counts(zz_meas)
plot_histogram(counts)

We can analyze the counts and calculate the ⟨𝑍𝑍⟩ as follows:

def measure_zz(given_circuit, num_shots = 10000):
zz_meas = measure_zz_circuit(given_circuit)
    
    result = execute(zz_meas, backend = simulator, shots = num_shots).result()
    counts = result.get_counts(zz_meas)
if '00' not in counts:
        counts['00'] = 0
    if '01' not in counts:
        counts['01'] = 0
    if '10' not in counts:
        counts['10'] = 0
    if '11' not in counts:
        counts['11'] = 0
total_counts = counts['00'] + counts['11'] + counts['01'] + counts['10']
    zz = counts['00'] + counts['11'] - counts['01'] - counts['10']
    zz = zz / total_counts
    
    return zz
zz = measure_zz(hets_circuit)
print("<ZZ> =", str(zz))

Out:

<ZZ> = -0.0036

What about ⟨𝑍𝐼⟩ and ⟨𝐼𝑍⟩? Do these need new circuits?

The answer is no, and they can be computed from the results above.

def measure_zi(given_circuit, num_shots = 10000):
    
    zz_meas = measure_zz_circuit(given_circuit)
    
    result = execute(zz_meas, backend = simulator, shots = num_shots).result()
    counts = result.get_counts(zz_meas)
if '00' not in counts:
        counts['00'] = 0
    if '01' not in counts:
        counts['01'] = 0
    if '10' not in counts:
        counts['10'] = 0
    if '11' not in counts:
        counts['11'] = 0
total_counts = counts['00'] + counts['11'] + counts['01'] + counts['10']
    
    zi = counts['00'] - counts['11'] + counts['01'] - counts['10']
    zi = zi / total_counts
    
    return zi
def measure_iz(given_circuit, num_shots = 10000):
    
    zz_meas = measure_zz_circuit(given_circuit)
    
    result = execute(zz_meas, backend = simulator, shots = num_shots).result()
    counts = result.get_counts(zz_meas)
if '00' not in counts:
        counts['00'] = 0
    if '01' not in counts:
        counts['01'] = 0
    if '10' not in counts:
        counts['10'] = 0
    if '11' not in counts:
        counts['11'] = 0
total_counts = counts['00'] + counts['11'] + counts['01'] + counts['10']
    
    iz = counts['00'] - counts['11'] - counts['01'] + counts['10']
    iz = iz / total_counts
    
    return iz
zi = measure_zi(hets_circuit)
print("<ZI> =", str(zi))
iz = measure_iz(hets_circuit)
print("<IZ> =", str(iz))

Out:

<ZI> = -0.0072
<IZ> = -0.0046

Next, we measure ⟨𝑋𝑋⟩

def measure_xx_circuit(given_circuit):
    xx_meas = given_circuit.copy()
    
    ### WRITE YOUR CODE BETWEEN THESE LINES - START
    given_circuit.h(0)
    given_circuit.h(1)
    xx_meas.measure_all()
    ### WRITE YOUR CODE BETWEEN THESE LINES - END
    return xx_meas
xx_meas = measure_xx_circuit(hets_circuit)
xx_meas.draw()

def measure_xx(given_circuit, num_shots = 10000):
    
    xx_meas = measure_xx_circuit(given_circuit)
    
    result = execute(xx_meas, backend = simulator, shots = num_shots).result()
    counts = result.get_counts(xx_meas)
if '00' not in counts:
        counts['00'] = 0
    if '01' not in counts:
        counts['01'] = 0
    if '10' not in counts:
        counts['10'] = 0
    if '11' not in counts:
        counts['11'] = 0
total_counts = counts['00'] + counts['11'] + counts['01'] + counts['10']
    xx = counts['00'] + counts['11'] - counts['01'] - counts['10']
    xx = xx / total_counts
    
    return xx
xx = measure_xx(hets_circuit)
print("<XX> =", str(xx))

Out:

<XX> = -1.0

Now we evaluate the energy of the trial state

def get_energy(given_circuit, num_shots = 10000):
    
    zz = measure_zz(given_circuit, num_shots = num_shots)
    iz = measure_iz(given_circuit, num_shots = num_shots)
    zi = measure_zi(given_circuit, num_shots = num_shots)
    xx = measure_xx(given_circuit, num_shots = num_shots)
    
    energy = (-1.0523732)*1 + (0.39793742)*iz + (-0.3979374)*zi + (-0.0112801)*zz + (0.18093119)*xx
    
    return energy
energy = get_energy(hets_circuit)
print("The energy of the trial state is", str(energy))

Out:

The energy of the trial state is -1.0541510196219999

Computing gradients

hets_circuit_plus = None
hets_circuit_minus = None
### WRITE YOUR CODE BETWEEN THESE LINES - START
def hets_circuit_plus(depth, angle1, angle2):
hets_circ_plus = QuantumCircuit(depth)
    hets_circ_plus.ry(angle1, 0)
    hets_circ_plus.rz(angle1, 0)
    hets_circ_plus.ry(angle1, 1)
    hets_circ_plus.rz(angle1, 1)
for ii in range(depth):
        hets_circ_plus.cx(0,1)
        hets_circ_plus.ry(angle2,0)
        hets_circ_plus.rz(angle2,0)
        hets_circ_plus.ry(angle2,1)
        hets_circ_plus.rz(angle2,1)
    
    return hets_circ_plus
def hets_circuit_minus(depth, angle1, angle2):
hets_circ_minus = QuantumCircuit(depth)
    hets_circ_minus.ry(angle1, 0)
    hets_circ_minus.rz(angle1, 0)
    hets_circ_minus.ry(angle1, 1)
    hets_circ_minus.rz(angle1, 1)
for ii in range(depth):
        hets_circ_minus.cx(0,1)
        hets_circ_minus.ry(angle2,0)
        hets_circ_minus.rz(angle2,0)
        hets_circ_minus.ry(angle2,1)
        hets_circ_minus.rz(angle2,1)
    
    return hets_circ_minus
hets_circuit_plus = hets_circuit_plus(depth=2,angle1 = pi/2 + 0.1*pi/2, angle2 = pi/2)
hets_circuit_minus = hets_circuit_minus(depth=2,angle1 = pi/2 - 0.1*pi/2, angle2 = pi/2)
print(hets_circuit_plus)
print(hets_circuit_minus)
### WRITE YOUR CODE BETWEEN THESE LINES - END

energy_plus = get_energy(hets_circuit_plus, num_shots=100000)
energy_minus = get_energy(hets_circuit_minus, num_shots=100000)
print(energy_plus, energy_minus)

Out:

-0.9600031637752001 -1.1304885633485997

As you can see, one of these is certainly lower energy than the other, and is also lower energy than the case when angle1 = pi/2. This is a suitable next point for our iteration of a variational eigensolver.

Bonus 1

While this is not graded, explore whether the decision above would be easy if your execution ran different numbers of shots. In particular, measure energy_plus and energy_minus again with 100, 1000 and 10000 shots to explore how easy or difficult this decision gets with each one.

energy_plus_100, energy_plus_1000, energy_plus_10000 = 0, 0, 0
energy_minus_100, energy_minus_1000, energy_minus_10000 = 0, 0, 0
### WRITE YOUR CODE BETWEEN THESE LINES - START
def hets_circuit_plus_100(depth, angle1, angle2):
hets_circ_plus_100 = QuantumCircuit(depth)
    hets_circ_plus_100.ry(angle1, 0)
    hets_circ_plus_100.rz(angle1, 0)
    hets_circ_plus_100.ry(angle1, 1)
    hets_circ_plus_100.rz(angle1, 1)
for ii in range(depth):
        hets_circ_plus_100.cx(0,1)
        hets_circ_plus_100.ry(angle2,0)
        hets_circ_plus_100.rz(angle2,0)
        hets_circ_plus_100.ry(angle2,1)
        hets_circ_plus_100.rz(angle2,1)
    
    return hets_circ_plus_100
def hets_circuit_plus_1000(depth, angle1, angle2):
hets_circ_plus_1000 = QuantumCircuit(depth)
    hets_circ_plus_1000.ry(angle1, 0)
    hets_circ_plus_1000.rz(angle1, 0)
    hets_circ_plus_1000.ry(angle1, 1)
    hets_circ_plus_1000.rz(angle1, 1)
for ii in range(depth):
        hets_circ_plus_1000.cx(0,1)
        hets_circ_plus_1000.ry(angle2,0)
        hets_circ_plus_1000.rz(angle2,0)
        hets_circ_plus_1000.ry(angle2,1)
        hets_circ_plus_1000.rz(angle2,1)
    
    return hets_circ_plus_1000
def hets_circuit_plus_10000(depth, angle1, angle2):
hets_circ_plus_10000 = QuantumCircuit(depth)
    hets_circ_plus_10000.ry(angle1, 0)
    hets_circ_plus_10000.rz(angle1, 0)
    hets_circ_plus_10000.ry(angle1, 1)
    hets_circ_plus_10000.rz(angle1, 1)
for ii in range(depth):
        hets_circ_plus_10000.cx(0,1)
        hets_circ_plus_10000.ry(angle2,0)
        hets_circ_plus_10000.rz(angle2,0)
        hets_circ_plus_10000.ry(angle2,1)
        hets_circ_plus_10000.rz(angle2,1)
    
    return hets_circ_plus_10000
def hets_circuit_minus_100(depth, angle1, angle2):
hets_circ_minus_100 = QuantumCircuit(depth)
    hets_circ_minus_100.ry(angle1, 0)
    hets_circ_minus_100.rz(angle1, 0)
    hets_circ_minus_100.ry(angle1, 1)
    hets_circ_minus_100.rz(angle1, 1)
for ii in range(depth):
        hets_circ_minus_100.cx(0,1)
        hets_circ_minus_100.ry(angle2,0)
        hets_circ_minus_100.rz(angle2,0)
        hets_circ_minus_100.ry(angle2,1)
        hets_circ_minus_100.rz(angle2,1)
    
    return hets_circ_minus_100
def hets_circuit_minus_1000(depth, angle1, angle2):
hets_circ_minus_1000 = QuantumCircuit(depth)
    hets_circ_minus_1000.ry(angle1, 0)
    hets_circ_minus_1000.rz(angle1, 0)
    hets_circ_minus_1000.ry(angle1, 1)
    hets_circ_minus_1000.rz(angle1, 1)
for ii in range(depth):
        hets_circ_minus_1000.cx(0,1)
        hets_circ_minus_1000.ry(angle2,0)
        hets_circ_minus_1000.rz(angle2,0)
        hets_circ_minus_1000.ry(angle2,1)
        hets_circ_minus_1000.rz(angle2,1)
    
    return hets_circ_minus_1000
def hets_circuit_minus_10000(depth, angle1, angle2):
hets_circ_minus_10000 = QuantumCircuit(depth)
    hets_circ_minus_10000.ry(angle1, 0)
    hets_circ_minus_10000.rz(angle1, 0)
    hets_circ_minus_10000.ry(angle1, 1)
    hets_circ_minus_10000.rz(angle1, 1)
for ii in range(depth):
        hets_circ_minus_10000.cx(0,1)
        hets_circ_minus_10000.ry(angle2,0)
        hets_circ_minus_10000.rz(angle2,0)
        hets_circ_minus_10000.ry(angle2,1)
        hets_circ_minus_10000.rz(angle2,1)
    
    return hets_circ_minus_10000
hets_circuit_plus_100 = hets_circuit_plus_100(depth=2,angle1 = pi/2 + 0.1*pi/2, angle2 = pi/2)
hets_circuit_plus_1000 = hets_circuit_plus_1000(depth=2,angle1 = pi/2 + 0.1*pi/2, angle2 = pi/2)
hets_circuit_plus_10000 = hets_circuit_plus_10000(depth=2,angle1 = pi/2 + 0.1*pi/2, angle2 = pi/2)
hets_circuit_minus_100 = hets_circuit_minus_100(depth=2,angle1 = pi/2 - 0.1*pi/2, angle2 = pi/2)
hets_circuit_minus_1000 = hets_circuit_minus_1000(depth=2,angle1 = pi/2 - 0.1*pi/2, angle2 = pi/2)
hets_circuit_minus_10000 = hets_circuit_minus_10000(depth=2,angle1 = pi/2 - 0.1*pi/2, angle2 = pi/2)
print(hets_circuit_plus_100)
print(hets_circuit_plus_1000)
print(hets_circuit_plus_10000)
print(hets_circuit_minus_100)
print(hets_circuit_minus_1000)
print(hets_circuit_minus_10000)
### WRITE YOUR CODE BETWEEN THESE LINES - END

energy_plus_100 = get_energy(hets_circuit_plus_100, num_shots=100)
energy_plus_1000 = get_energy(hets_circuit_plus_100, num_shots=1000)
energy_plus_10000 = get_energy(hets_circuit_plus_100, num_shots=10000)
energy_minus_100 = get_energy(hets_circuit_minus_100, num_shots=100)
energy_minus_1000 = get_energy(hets_circuit_minus_100, num_shots=1000)
energy_minus_10000 = get_energy(hets_circuit_minus_100, num_shots=10000)
print(energy_plus_100, energy_minus_100, "difference = ", energy_minus_100 - energy_plus_100)
print(energy_plus_1000, energy_minus_1000, "difference = ", energy_minus_1000 - energy_plus_1000)
print(energy_plus_10000, energy_minus_10000, "difference = ", energy_minus_10000 - energy_plus_10000)

Out:

-0.9122630381999999 -1.1198517056 difference =  -0.20758866740000004
-0.4401202818199999 -0.45621825961999984 difference =  -0.01609797779999994
-0.96466067423 -1.1187116705179998 difference =  -0.1540509962879998

Bonus 2

While this is not graded, diagonalize the Hamiltonian by writing down the matrices for the Pauli operators I, X and Z, and find the exact ground state energy.

### WRITE YOUR CODE BETWEEN THESE LINES - START
import numpy as np
from qiskit.quantum_info.operators import Operator, Pauli
h = [[-1.0523732, Pauli(label='II')], 
     [0.39793742, Pauli(label='IZ')], 
     [-0.3979374, Pauli(label='ZI')], 
     [-0.0112801, Pauli(label='ZZ')],
     [0.18093119, Pauli(label='XX')]
    ]
H_m = []
for hi in h:
    H_m.append(hi[0]*Operator(hi[1]))
np.linalg.eigvalsh((H_m[0]+H_m[1]+H_m[2]+H_m[3]+H_m[4]).data)
from qiskit.aqua.algorithms import NumPyEigensolver
from qiskit.aqua.operators import WeightedPauliOperator
H = WeightedPauliOperator(h)
print(H.print_details())
ee = NumPyEigensolver(H)
result = ee.run()
ref = result['eigenvalues']
print(ref)
### WRITE YOUR CODE BETWEEN THESE LINES - END

Out:

II	(-1.0523732+0j)
IZ	(0.39793742+0j)
ZI	(-0.3979374+0j)
ZZ	(-0.0112801+0j)
XX	(0.18093119+0j)

[-1.85727496-1.35566063e-16j]

References:

Qiskit (Global Summer School), Introduction to Quantum Computing and Quantum Hardware — Lab 9.