Implementing the VQE

For a one-qubit system we can reach every point on the Bloch sphere (as discussed earlier) with a rotation about the $ x $-axis and the $ y $-axis.

We can express this mathematically through the following operations (see whiteboard for the drawing), giving us a new state $ \vert \psi\rangle $

$$ \vert\psi\rangle = R_y(\phi)R_x(\theta)\vert 0 \rangle. $$

We can produce multiple ansatzes for the new state in terms of the angles $ \theta $ and $ \phi $. With these ansatzes we can in turn calculate the expectation value of the above Hamiltonian, now rewritten in terms of various Pauli matrices (and thereby gates), that is compute

$$ \langle \psi \vert (c+\mathcal{E})\boldsymbol{I} + (\Omega+\omega_z)\boldsymbol{\sigma}_z + \omega_x\boldsymbol{\sigma}_x\vert \psi \rangle. $$

We can now set up a series of ansatzes for $ \vert \psi \rangle $ as function of the angles $ \theta $ and $ \phi $ and find thereafter the variational minimum using for example a gradient descent method.

To do so, we need to remind ourselves about the mathematical expressions for the rotational matrices/operators.

$$ R_x(\theta)=\cos{\frac{\theta}{2}}\boldsymbol{I}-\imath \sin{\frac{\theta}{2}}\boldsymbol{\sigma}_x, $$

and

$$ R_y(\phi)=\cos{\frac{\phi}{2}}\boldsymbol{I}-\imath \sin{\frac{\phi}{2}}\boldsymbol{\sigma}_y. $$

# define the rotation matrices
# Define angles theta and phi
theta = 0.5*np.pi; phi = 0.2*np.pi
Rx = np.cos(theta*0.5)*I-1j*np.sin(theta*0.5)*X
Ry = np.cos(phi*0.5)*I-1j*np.sin(phi*0.5)*Y
#define basis states
basis0 = np.array([1,0])
basis1 = np.array([0,1])

NewBasis = Ry @ Rx @ basis0
print(NewBasis)
# Compute the expectation value
#Note hermitian conjugation
Energy = NewBasis.conj().T @ Hamiltonian @ NewBasis
print(Energy)

Not an impressive results. We set up now a loop over many angles $ \theta $ and $ \phi $ and compute the energies

# define a number of angles
n = 20
angle = np.arange(0,180,10)
n = np.size(angle)
ExpectationValues = np.zeros((n,n))
for i in range (n):
    theta = np.pi*angle[i]/180.0
    Rx = np.cos(theta*0.5)*I-1j*np.sin(theta*0.5)*X
    for j in range (n):
        phi = np.pi*angle[j]/180.0
        Ry = np.cos(phi*0.5)*I-1j*np.sin(phi*0.5)*Y
        NewBasis = Ry @ Rx @ basis0
        Energy = NewBasis.conj().T @ Hamiltonian @ NewBasis
        Edifference=abs(np.real(EigValues[0]-Energy))
        ExpectationValues[i,j]=Edifference

print(np.min(ExpectationValues))

Clearly, this is not the very best way of proceeding. Rather, here we could try to find the optimal values for the parameters $ \theta $ and $ \phi $ through computation of their respective gradients and thereby find the minimum as function of the optimal angles $ \hat{\theta} $ and $ \hat{\phi} $.

Let us now implement a classical gradient descent algorithm to the computation of the energies. We will follow closely https://journals.aps.org/pra/abstract/10.1103/PhysRevA.99.032331 in order to calculate gradients of the Hamiltonian.