Variational Quantum Eigensolver (VQE)

Note

Source code is embedded into the following sections.

Introduction

In 2013, Alberto Peruzzo et al. [20] proposed the Variational Quantum Eigensolver (VQE). VQE is a hybrid algorithm that runs partly on a quantum computer and partially on a classical computer. The algorithm efficiently calculates the expectation value of the Hermitian operator \(H\) in a quantum computer using a parametrization circuit of a polynomial depth. It works in tandem with a classical computer that calculates the eigenvalues and eigenvectors of \(H\). The algorithm proceeds by varying the parameters that control the quantum state of \(H\) using classical optimizers, systematically finding the minimum energy eigenstate. In this section, we shall examine how to implement a VQE using the Quantum Rings SDK.

The variational principle

Some of the problems that we want to solve on a quantum computer are not exactly solvable, and in many cases, an approximate solution is acceptable. Where we can’t find an analytic solution to the Schrödinger equation, two methods are employed in quantum mechanics – the perturbation theory and the variational principle.

The Hamiltonians corresponding to the hydrogen atom, the quantum harmonic oscillator etc., are the simplest form and hence can be solved exactly. The perturbation theory helps us to derive the Hamiltonians for most complex systems from such ideal forms by applying a small dimensionless term (the ‘perturbing’ Hamiltonian, which describes a weak disturbance). When this small term is zero, the system is exactly solvable. The system is studied by expanding the correctional term in a finite power series. The expansion is asymptotic, however, if the expansion term is small, the series can converge with a finite order. In most cases, after a certain order, the expansion is usually divergent.

With the variational theory, we can get convergent expansions. It is particularly useful in studying ground state of quantum systems. It is not required that the correction parameter is small or that the system is solvable in a certain limit.

The variational theorem states that:

Variational Theorem

Assume a system with a time-independent Hamiltonian \(H\). If \(𝜓\) is a normalized wavefunction (also known as the ansatz or the trial wavefunction) of the system that satisfies the boundary conditions of the problem, then,

\(\left\langle \psi\left| H \right| \psi \right\rangle \ge E_0\)

where \(E_0\) is the lowest eigenvalue (that is, the ground state energy) of \(H\).

According to the variational principle, by varying the trial wavefunction \(𝜓\) until the expectation value of \(𝐻\) is minimized, we can obtain an approximation to the wavefunction and energy of the ground-state. In the next few steps, we shall establish how this works. For now, we shall prove the variational theorem itself.

Let us assume that \(\{ |\psi_n \rangle \}\) is an orthonormal eigenbasis of \(H\), such that:

\[H |\psi_n \rangle = E_n |\psi_n \rangle -------(1)\]

Then, the trial wavefunction can be written as a linear combination of the orthonormal eigenbasis, using the quantum superposition principle.

\[\psi = \sum_{n}^{} c_n |\psi_n \rangle -------(2)\]

This leads us to the linear variational method, with \(𝑐_𝑛\) being the variational parameters. Note that, if the trial wavefunction is exact when \(𝑐_0=1\) and \(𝑐_𝑛=0\); \(∀ 𝑛 ≠0\). If the trial wavefunction is not exact, then it is a linear combination of different eigenstates. Note that normalization requires,

\[\psi = \sum_{n}^{} \left| c_n \right|^2 = 1 -------(3)\]

We can rewrite the above equation as:

\[|c_0|^2 = 1 -\sum_{n>0}^{} \left| c_n \right|^2 -------(4)\]

The expectation value of the energy of the trial wavefunction can be written as follows:

\[\left\langle H \right\rangle = \left\langle \psi\left| H \right| \psi \right\rangle = \ \left( \sum_{n}^{} c_n^* \langle \psi_n | \right) H \ \left( \sum_{m}^{} c_m |\psi_m \rangle \right)\]

with some algebra, this simplifies to:

\[\left\langle H \right\rangle = \sum_{n}^{} \left| c_n \right|^2 E_n -------(5)\]

Using (4),

\[ \begin{align}\begin{aligned}\left\langle \psi\left| H \right| \psi \right\rangle &=\ \left| c_0 \right|^2 E_0 + \sum_{n>0}^{} \left| c_n \right|^2 E_n\\& = E_0 + \sum_{n>0}^{} \left| c_n \right|^2 (E_n - E_0) -------(6)\end{aligned}\end{align} \]

Note that, \(𝐸_0≤ 𝐸_1≤ 𝐸_2≤⋯ ≤𝐸_𝑛\), where \(𝐸_0\) is the energy corresponding to the ground state. Hence \((𝐸_𝑛−𝐸_0) ≥0\); \(∀ 𝑛>0\). Therefore, the second term of the above equation is positive definite. Thus, we can say that:

\[\left\langle \psi\left| H \right| \psi \right\rangle \ge E_0 -------(7)\]

which proves the theorem, establishing that the expectation value of the energy is always greater than or equal to the ground state energy. In other words, the lowest expectation values of the energy must have a good overlap with the ground state energy.

The quantum circuit

Let us now construct a quantum circuit. Let \(𝐻\) be the Hamiltonian representing the total energy of the quantum state of the quantum circuit. Let \(𝜃={𝜃_1,𝜃_2,..,𝜃_𝑛 }\) be a vector of a polynomial number of circuit parameters that can determine the quantum state of the quantum circuit. Assume that there exists a unitary \(𝑈(𝜃)\), such that a trial quantum state \(|\psi(\theta) \rangle = U(\theta) |\psi\rangle\) can be obtained by the action of the unitary \(𝑈(𝜃)\) on an initial state \(|\psi\rangle\). We can write the expectation value of energy over a state controlled by the parameter \(𝜃\) as follows:

\[E\left(\theta\right) = \left\langle H \right\rangle = \left\langle \psi(\theta)\left|H \right| \psi(\theta) \right \rangle -------(8)\]

The next step is to obtain a ground-state wavefunction by looking for the lowest possible expectation value by minimizing it. This process shall become clear in the following steps. Note that, equation (1) is an eigenvalue equation and it has several solutions \(|\psi_i\rangle\). For every solution \(|\psi_\rangle\), there is an eigenvalue \(𝐸_𝑖\), corresponding to the eigenstate \(|\psi_\rangle\). One of them corresponds to the ground-state. The spectral decomposition of \(𝐻\) can be written as:

\[H = \sum_{i=1}^{n}E_i \langle \psi_i \left| H \right| \psi_i \rangle -------(9)\]

Plugging (9) into (8), we get:

\[ \begin{align}\begin{aligned}\left\langle H \right\rangle &= \left\langle \psi(\theta)\left| \sum_{i}^{n} E_i \right| \psi_i \right\rangle \left\langle \psi_i| \psi(\theta) \right\rangle\\&= \sum_{i}^{n} E_i \left\langle \psi(\theta) | \psi_i \right\rangle \left\langle \psi_i| \psi(\theta) \right\rangle\\&= \sum_{i}^{n} E_i \left| \left\langle \psi_i| \psi(\theta) \right\rangle \right|^2 -------(10)\end{aligned}\end{align} \]

Hence the expectation value \(〈𝐻〉\) can be written as a linear combination of eigenvalues. Since \(\left| \left\langle \psi_i| \psi(\theta) \right\rangle \right|^2 \ge 0\), the minimum value of equation (10) must be higher than or equal to a certain \(𝐸_{𝑚𝑖𝑛}\), that is,

\[E_{min} \le \left\langle H \right\rangle = \left\langle \psi(\theta)\left| H \right| \psi(\theta) \right\rangle -------(11)\]

This means that the expectation value of the wavefunction is at least equal to the minimum eigenvalue. This also means that,

\[E_{min} = \left\langle \psi_{min}\left| H \right| \psi_{min} \right\rangle -------(12)\]

The variational approach starts with an arbitrary \(𝜃\) and computes \(\left\langle \psi(\theta)\left| H \right| \psi(\theta) \right\rangle\). By varying \(𝜃\), the variational approach attempts to identify the lowest energy eigenstate \(\left|\psi_{min}\right\rangle\) and performs the following minimization:

\[\underset{\theta}{min}= \left\langle \psi(\theta)\left| H \right| \psi(\theta) \right\rangle -------(13)\]

The unitary matrix \(𝑈(𝜃)\) is implemented as a Parameterized quantum circuit and it is often called the variational form or ansatz. Construction of the ansatz using a polynomial number of operations is essential so that the classical optimizing algorithms can define \(\left|\psi\right\rangle\) from a polynomial search space [22] [4].

The following graph illustrates the VQE algorithm.

VQE Algorithm.

VQE has several applications, such as molecular simulations in quantum chemistry and solving optimization problems. With this introduction, we shall now examine the inner workings of the VQE using an example.

Open a new Python 3 notebook in Jupyter and enter following code that imports the components we shall be using in this example.

Note

Be sure to use your API token and your account name. This module requires scipy optimizer. If you have not installed it before, use the following command from the terminal:

pip install scipy

import QuantumRingsLib
from QuantumRingsLib import QuantumRegister, AncillaRegister, ClassicalRegister, QuantumCircuit
from QuantumRingsLib import QuantumRingsProvider
from QuantumRingsLib import job_monitor
from QuantumRingsLib import JobStatus
from QuantumRingsLib import qasm2
from QuantumRingsLib import Parameter, ParameterVector

from random import random, uniform
from matplotlib import pyplot as plt
import numpy as np
import math
import time

import scipy
from scipy import optimize

provider = QuantumRingsProvider(token ="YOUR_TOKEN", name="YOUR_ACCOUNT_EMAIL")
backend = provider.get_backend("scarlet_quantum_rings")
shots = 2

print(provider.active_account())

Constructing the ansatz using RY gates

The first step is to construct a Parameterized quantum circuit (PQC), also called the ansatz. The PQC is a unitary operation \(𝑈(𝜃)\) on \(𝑛\) qubits, set in an initial state \(\left|\theta\right\rangle\). \(\theta\) refers to a tunable circuit parameter; in our example it is a vector of rotational angles of the RY gates used in the circuit. As shown in the workflow diagram, an objective function is evaluated based on the ansatz described by the corresponding Ising Hamiltonian. A classical minimizer optimizes \(\theta\) to get a minimal expectation value from the objective function. In literature [22], several Parameterized quantum circuits have been explored to encode the solution space with a reduced circuit depth.

In our example, we construct the ansatz using RY gates as shown in the following figure:

       ┌──────────────┐ ╎ ┌──────────────┐ ╎ ┌──────────────┐   ╎ ┌────────────────┐ ╎ ┌────────────────┐ ╎
q[0]: ■┤ RY(theta[0]) ├─╎─┤ RY(theta[4]) ├─╎─┤ RY(theta[8]) ├───╎─┤ RY(theta[12])  ├─╎─┤ RY(theta[16])  ├─╎─
       ├──────────────┤ ╎ ├──────────────┤ ╎ ├──────────────┤   ╎ ├────────────────┤ ╎ ├────────────────┤ ╎
q[1]: ■┤ RY(theta[1]) ├─╎─┤ RY(theta[5]) ├─╎─┤ RY(theta[9]) ├───╎─┤ RY(theta[13])  ├─╎─┤ RY(theta[17])  ├─╎─
       ├──────────────┤ ╎ ├──────────────┤ ╎ ├──────────────┴─┐ ╎ ├────────────────┤ ╎ ├────────────────┤ ╎
q[2]: ■┤ RY(theta[2]) ├─╎─┤ RY(theta[6]) ├─╎─┤ RY(theta[10])  ├─╎─┤ RY(theta[14])  ├─╎─┤ RY(theta[18])  ├─╎─
       ├──────────────┤ ╎ ├──────────────┤ ╎ ├────────────────┤ ╎ ├────────────────┤ ╎ ├────────────────┤ ╎
q[3]: ■┤ RY(theta[3]) ├─╎─┤ RY(theta[7]) ├─╎─┤ RY(theta[11])  ├─╎─┤ RY(theta[15])  ├─╎─┤ RY(theta[19])  ├─╎─
       └──────────────┘ ╎ └──────────────┘ ╎ └────────────────┘ ╎ └────────────────┘ ╎ └────────────────┘ ╎
c: 4/ ■═════════════════════════════════════════════════════════════════════════════════════════════════════

Note that the circuit shown in the above figure can be optimized by summing up the rotational angles pertaining to the respective qubits.

Type the code in the following listing in a new cell to implement the ansatz circuit:

# Implements a simple ansatz circuit using RY rotations
# Inputs: qc – QuantumCircuit
# theta_list – parameter list must be (reps + 1) * n_qubits wide.
# q – the qubit register to use
# n_qubits – Number of qubits to use
# reps – repeats of the circuit, defaults to 5
# insert_barriers –  whether to create a barrier gate
# Returns: Inplace modified QuantumCircuit

def SimpleAnsatz ( qc, q, n_qubits, theta_list, reps = 5, insert_barriers=False):
    for i in range (reps+1):
        theta = 0
        for j in range (n_qubits):
            #theta += theta_list[(i * n_qubits)+j]
            qc.ry(theta_list[(i * n_qubits)+j], q[j])

        if ( True == insert_barriers):
                qc.barrier()
    return

Implementing a 2-local ansatz

Another way of implementing the ansatz is to use a 2-local circuit [1], as often illustrated in literature. Implementing the 2-local ansatz circuit is straightforward. The core circuit can be repeated by altering the input parameter reps, as shown in the listing below.

Type the code in the following listing in a new cell to implement the 2-local ansatz circuit:

# Implements the ansatz circuit using RY rotations and CZ gates in a liner arrangement
# Inputs: qc – QuantumCircuit
# theta_list – rotational angle parameter must be (reps + 1) * n_qubits wide.
# q – the qubit register to use
# n_qubits – Number of qubits to use
# reps – repeats of the circuit, defaults to 5
# insert_barriers – whether to create a barrier gate after each rep.
# Returns: Inplace modified QuantumCircuit

def TwoLocalAnsatz ( qc, q, n_qubits, theta_list, reps=5, insert_barriers=False):
    i_theta = 0

    for i in range(n_qubits):
        qc.ry(theta_list[i_theta], q[i])
        i_theta += 1

    for _ in range(reps):
        if ( True == insert_barriers):
            qc.barrier()

        for i in range(n_qubits-1):
            qc.cz(q[i], q[i+1])

        for i in range(n_qubits):
            qc.ry(theta_list[i_theta], q[i])
            i_theta += 1
    return

2-local ansatz circuit:

       ┌──────────────┐ ╎                ┌──────────────┐ ╎                ┌──────────────┐   ╎                »
q[0]: ■┤ RY(theta[0]) ├─╎───■────────────┤ RY(theta[4]) ├─╎───■────────────┤ RY(theta[8]) ├───╎───■────────────»
       ├──────────────┤ ╎ ┌─┴─┐          ├──────────────┤ ╎ ┌─┴─┐          ├──────────────┤   ╎ ┌─┴─┐          »
q[1]: ■┤ RY(theta[1]) ├─╎─┤ Z ├──■───────┤ RY(theta[5]) ├─╎─┤ Z ├──■───────┤ RY(theta[9]) ├───╎─┤ Z ├──■───────»
       ├──────────────┤ ╎ └───┘┌─┴─┐     ├──────────────┤ ╎ └───┘┌─┴─┐     ├──────────────┴─┐ ╎ └───┘┌─┴─┐     »
q[2]: ■┤ RY(theta[2]) ├─╎──────┤ Z ├──■──┤ RY(theta[6]) ├─╎──────┤ Z ├──■──┤ RY(theta[10])  ├─╎──────┤ Z ├──■──»
       ├──────────────┤ ╎      └───┘┌─┴─┐├──────────────┤ ╎      └───┘┌─┴─┐├────────────────┤ ╎      └───┘┌─┴─┐»
q[3]: ■┤ RY(theta[3]) ├─╎───────────┤ Z ├┤ RY(theta[7]) ├─╎───────────┤ Z ├┤ RY(theta[11])  ├─╎───────────┤ Z ├»
       └──────────────┘ ╎           └───┘└──────────────┘ ╎           └───┘└────────────────┘ ╎           └───┘»
c: 4/ ■════════════════════════════════════════════════════════════════════════════════════════════════════════»
                                                                                                               »
«       ┌────────────────┐ ╎                ┌────────────────┐ ╎                ┌────────────────┐
«q[0]: ■┤ RY(theta[12])  ├─╎───■────────────┤ RY(theta[16])  ├─╎───■────────────┤ RY(theta[20])  ├
«       ├────────────────┤ ╎ ┌─┴─┐          ├────────────────┤ ╎ ┌─┴─┐          ├────────────────┤
«q[1]: ■┤ RY(theta[13])  ├─╎─┤ Z ├──■───────┤ RY(theta[17])  ├─╎─┤ Z ├──■───────┤ RY(theta[21])  ├
«       ├────────────────┤ ╎ └───┘┌─┴─┐     ├────────────────┤ ╎ └───┘┌─┴─┐     ├────────────────┤
«q[2]: ■┤ RY(theta[14])  ├─╎──────┤ Z ├──■──┤ RY(theta[18])  ├─╎──────┤ Z ├──■──┤ RY(theta[22])  ├
«       ├────────────────┤ ╎      └───┘┌─┴─┐├────────────────┤ ╎      └───┘┌─┴─┐├────────────────┤
«q[3]: ■┤ RY(theta[15])  ├─╎───────────┤ Z ├┤ RY(theta[19])  ├─╎───────────┤ Z ├┤ RY(theta[23])  ├
«       └────────────────┘ ╎           └───┘└────────────────┘ ╎           └───┘└────────────────┘
«c: 4/ ■══════════════════════════════════════════════════════════════════════════════════════════
«

The 2-local circuit described above uses a CZ gate in a linear arrangement to create the entanglement between the alternating layers of RY gates (see figure). Other gate combinations and entanglement arrangements together with randomized order of rotational angles can be tried for specific problems.

Both the ansatzes can be equivalently tried in this example.

Setting up the Ising Hamiltonian

To illustate the VQE example, we shall be directly creating the Pauli operators corresponding to a toy Ising Hamiltonian to illustrate the workflow.

Consider the following example problem Hamiltonian which is a sum of tensor products of Pauli operators.

\[𝐻 = 2.0*𝐼𝐼𝐼𝑍 + 4.0*𝐼𝐼𝑍𝑍 + 8.0∗𝐼𝑍𝐼𝐼 + 16.0∗𝑍𝐼𝐼𝑍 -------(14)\]

Recall equations (11), (12), and (13). Suppose we are able to construct an ansatz \(\left|\psi(\theta)\right\rangle\) which drives the state corresponding to the minimum eigenvalue. In that case, the expectation value of the minimum Hamiltonian can be obtained by measuring the quantum circuit in the measurement basis several times and averaging out the results.

Having constructed the ansatz, the second step of the algorithm is to construct the Pauli operators from the Ising Hamiltonian. Type the following code in the Jupyter notebook:

# setup the Pauli operators corresponding to the problem Hamiltonian
# Inputs: nothing
# Returns: a list of Pauli operators and the corresponding weighing coefficients corresponding to the Hamiltonian


def get_Paulioperator():
    # Assume the following Hamiltonian
    #<H> = 2*<ψ(θ)|H1|ψ(θ)> + 4*<ψ(θ)|H2|ψ(θ)> + 8*<ψ(θ)|H3|ψ(θ)> + 16*<ψ(θ)|H4|ψ(θ)>
    #H1 = IIIZ
    #H2 = IIZZ
    #H3 = IZII
    #H4 = ZIIZ

    pauli_list = []
    w=[2+0j, 4+0j, 8+0j, 16+0j]

    pauli_list.append([w[0], "IIIZ"])
    pauli_list.append([w[1], "IIZZ"])
    pauli_list.append([w[2], "IZII"])
    pauli_list.append([w[3], "ZIIZ"])

    return pauli_list

Note that, there are some tools to create the Ising Hamiltonian and derive the Pauli operators. This is a toy example and we shall code them directly as illustrated above.

Recall the Pauli operator group \(𝐺=𝐼,𝑋,𝑌,𝑍\). For a quantum system of \(𝑛\) qubits, the Pauli operators \(𝑃_𝑛\) are a collection of \(4^{𝑛+1}\) elements of the form:

\[P_n = i^l \sigma_1 \otimes \sigma_2 \otimes ... \otimes \sigma_n -------(15)\]

where \(\sigma_i∈𝐼,𝑋,𝑌,𝑍\) and \(𝑙 ∈{1,2,3,4}\). For a 1-qubit system, the Pauli operators are defined by two-dimensional vectors:

\[ \begin{align}\begin{aligned}& I \to \sigma_{00} \equiv (0,0)\\& X \to \sigma_{01} \equiv (0,1)\\& Y \to \sigma_{10} \equiv (1,0)\\& Z \to \sigma_{11} \equiv (1,1) -------(16)\end{aligned}\end{align} \]

For a \(𝑛\)-qubit system, it is a \(2𝑛\) dimensional vector:

\[\sigma_{z_1x_1} \otimes \sigma_{z_2x_2} \otimes ...\otimes \sigma_{z_nx_n} \equiv \left\{ z_1, z_2, ..., z_n, x_1, x_2, ..., x_n \right\} -------(17)\]

The action of Pauli operators on tensor products can be derived easily using Dirac’s bra-ket notation. The following equation provides an example.

\[IXZX (\left| 0101 \right\rangle - \left| 1010 \right\rangle) = \left| 0000 \right\rangle - \left| 1111 \right\rangle -------(18)\]

To perform Pauli measurements, we have to transform the Pauli operators into the computation basis. Assume that we are measuring a certain qubit in state \(\left|\psi\right\rangle\) in the computation basis. The probability of projecting the qubit to a \(\left|0\right\rangle\) or \(\left|1\right\rangle\) state is:

\[ \begin{align}\begin{aligned}& P(0) = \left| \left\langle 0 | \psi\right\rangle \right|^2 = \left\langle \psi|0 \right\rangle \left\langle 0 | \psi \right\rangle\\& P(1) = \left| \left\langle 1 | \psi\right\rangle \right|^2 = \left\langle \psi|1 \right\rangle \left\langle 1 | \psi \right\rangle -------(19)\end{aligned}\end{align} \]

Recall that the Pauli matrices \(\sigma_I\) and \(\sigma_Z\) can be written in terms of outer products:

\[ \begin{align}\begin{aligned}& \sigma_I = I = \left|0\right\rangle \left\langle 0 \right| + \left|1\right\rangle \left\langle 1 \right|\\& \sigma_Z = Z = \left|0\right\rangle \left\langle 0 \right| - \left|1\right\rangle \left\langle 1 \right| -------(20)\end{aligned}\end{align} \]

Their expectation values are:

\[ \begin{align}\begin{aligned}\left\langle \psi \left| I \right| \psi\right\rangle &=\ \left\langle \psi \right|\left( \left|0\right\rangle \left\langle 0 \right| + \left|1\right\rangle \left\langle 1 \right| \right) \left| \psi \right\rangle\\& = \left\langle \psi|0 \right\rangle \left\langle 0 | \psi \right\rangle + \left\langle \psi|1 \right\rangle \left\langle 1 | \psi \right\rangle\\& = P\left( 0 \right) + P\left( 1 \right) -------(21)\end{aligned}\end{align} \]

Simillarly,

\[ \begin{align}\begin{aligned}\left\langle \psi \left| Z \right| \psi\right\rangle &=\ \left\langle \psi \right|\left( \left|0\right\rangle \left\langle 0 \right| - \left|1\right\rangle \left\langle 1 \right| \right) \left| \psi \right\rangle\\& = \left\langle \psi|0 \right\rangle \left\langle 0 | \psi \right\rangle - \left\langle \psi|1 \right\rangle \left\langle 1 | \psi \right\rangle\\& = P\left( 0 \right) - P\left( 1 \right) -------(22)\end{aligned}\end{align} \]

Ignoring the phase factor, which cannot be measured, the expectation values of the operators \(𝐼\) and \(𝑍\) are the probabilities of measuring the state \(\left|\psi\right\rangle\) in the computation basis. Hence, to measure the operators \(𝐼\) and \(𝑍\), we prepare the system using the ansatz and measure them directly in the computational basis.

The operators \(𝑋\) and \(𝑌\) cannot be directly written as outer products in the computation basis. However, we know that \(\sigma_x=𝐻𝑍𝐻\) and \(\sigma_y=𝑆𝑋𝑆^{†}= 𝑆𝐻𝑍𝐻𝑆^{†}\). Therefore, the expectation values of the operators \(𝑋\) and \(𝑌\) can be written in terms of the operators \(𝐼\) and \(𝑍\). This is explained in the following steps:

\[ \begin{align}\begin{aligned}& \left\langle \psi \left| X \right| \psi\right\rangle = \left\langle \psi \left| HZH \right| \psi\right\rangle\\& \left\langle \psi \left| Y \right| \psi\right\rangle = \left\langle \psi \left| SHZHS^{\dagger} \right| \psi\right\rangle -------(23)\end{aligned}\end{align} \]

Therefore, to measure \(𝑋\) in the computation basis, we must apply an H-gate to the corresponding qubit. Similarly, to measure \(𝑌\) in the computation basis, we must apply an \(𝑆^{†}\) gate followed by an \(H\)-gate to the corresponding qubit. The following table summarizes these transformations, which we shall use in the variational circuit to measure the observable \(\left\langle H \right\rangle\).

Table 1. Pauli measurements in computational basis

Operator

Unitary Transformations

\(I\)

\(I\)

\(𝑋\)

\(H\)

\(Y\)

\(HS^{\dagger}\)

\(Z\)

\(I\)

Performing the Pauli measurements

The next step is to measure \(\left\langle H \right\rangle\). This step is done by applying each of the weighted Pauli operators to the ansatz and measuring the circuit in the computational basis. If any of the Pauli operators are in \(𝑌\) or \(𝑋\) basis, we need to perform a basis conversion as illustrated in Table 1 before the respective qubits can be measured in the computational basis.

The iterative method of Pauli measurements adds up a complexity of \(𝑂(𝑝)\), \(𝑝\) being the number of Pauli operator vectors. While this is not so efficient, a direct optimization can be done by performing this step in parallel (using multitasking), and we leave it as an exercise. However, this optimization is an area of research [3], and some performance can be obtained by grouping the operators into commuting families [7]. For example, the operators \(𝐼\) and \(𝑍\) share the same eigenbasis (since they commute), and some reductions can be applied.

For now, we shall follow the original research paper. Type the following code fragments in the Jupyter notebook in a new cell.

# Performs pauli measurements
# qubitop: The Pauli operator
# param_dict: list of theta's for the parametrization circuit
# SHOTS: number of times, the circuit is to be repeated
def perform_pauli_measurements( qubitOp, param_dict, SHOTS=1024):

    avg = 0.0
    n_qubits = len(qubitOp[0][1])
    pauli_list = qubitOp


    # for each Pauli operator
    for p in pauli_list:
        weight = p[0].real
        pauli  = p[1]

        # assign parameters to the pqc and clone it
        qc = vqe_pqc.assign_parameters(param_dict)

        # Apply the Pauli operators.
        # no actions for "I" or "Z"
        for i in range(n_qubits):
            if (pauli[i] == "Y"):
                qc.sdg(q[i])
                qc.h(q[i])
            elif (pauli[i] == "X"):
                qc.h(q[i])

        # We should measure this circuit in the computation basis now
        qc.measure_all()

        job = backend.run(qc, shots= SHOTS, mode="sync", performance="HighestEfficiency", quiet=True)
        job_monitor(job, quiet=True)

        results = job.result()
        result_dict = results.get_counts()

        # perform the pauli measurement
        # convert the operator into binary
        measurement = 0.0
        pauli_int = int (p[1].replace("I","0").replace("Z","1").replace("X","1").replace("Y","1"),2)

        for key, value in result_dict.items():
            sign = -1.0 if ( bin(int(key,2) & pauli_int).count("1") & 1 ) else 1.0
            measurement += sign * value
        measurement /= SHOTS

        measurement = measurement * weight
        avg = avg + measurement

    return avg


# Given a theta_list, calculate the eigenstate
# Performs pauli measurements
# param_dict: list of theta's for the parametrization circuit
# SHOTS: number of times, the circuit is to be repeated
def find_eigenstate (param_dict, SHOTS=1024 ):

    # assign parameters to the pqc and clone it
    qc = vqe_pqc.assign_parameters(param_dict)

    # We should measure the circuit now
    qc.measure_all()

    job = backend.run(qc, shots= SHOTS, mode="sync", performance="HighestEfficiency", quiet=True)
    job_monitor(job, quiet=True)
    results = job.result()
    result_dict = results.get_counts()

    sorted_list = sorted(result_dict.items(), key=lambda value: value[1])

    return sorted_list[-1][0]

Putting everything together

Let us now complete the rest of the code by creating the Parameterized VQE circuit and the main routine that runs the classical optimizer. Type in the code contained in the following fragments in the Jupyter notebook.

# Create the Parameterized Quantum Circuit
# Input:
# qubitop: The Pauli Operator
# theta_list: Parameter vector
# layers: number of circuit layers to use in the ansatz
# Returns: the Parameterized quantum circuit
# try both the ansatzes by uncommenting one of them

def create_ParameterizedVQECircuit(qubitOp, theta_list, layers):
    n_qubits = len(qubitOp[0][1])

    # construct the ansatz
    q = QuantumRegister(n_qubits, "q")
    c = ClassicalRegister(n_qubits, "c")
    circ = QuantumCircuit(q,c)

    #SimpleAnsatz (circ, q, n_qubits, theta_list, reps = layers)
    TwoLocalAnsatz (circ, q, n_qubits, theta_list, reps = layers)


    return circ

The main routine is listed below:

# Define the Pauli operator for the toy Hamiltonian described above
qubitOp = get_Paulioperator() # Setup the Pauli operators

# declare some constants
n_layers = 5 # Number of layers in the ansatz
n_calls = 300 # Number of optimizer iterations
n_qubits = len(qubitOp[0][1]) # Number of qubits required
n_params = (n_layers + 1) * n_qubits # Size of Parameter array
theta_list = ParameterVector("theta", n_params)

# Initial theta_list (parameters)
initial_parameters_list = []
for i in range (n_params):
    initial_parameters_list.append(uniform(-2*math.pi, 2*math.pi))

# create the Parameterized circuit
vqe_pqc = create_ParameterizedVQECircuit(qubitOp, theta_list, n_layers)

avg_list = []

# The VQE routine called by the optimizer
# theta_list: the parameter list to be tried in the current cycle

def vqe(param_list) -> float:
    # create the parameter dictionary
    param_dict = {}
    for index in range (0, n_params):
        param_dict[theta_list[index].name()] = param_list[index]

    avg = perform_pauli_measurements (qubitOp, param_dict)
    #print(theta_list)
    #print(avg)
    avg_list.append(avg)
    return avg

# Method to find the eigenstate, for a given theta_list
# theta_list: The parameter array
def vqe_eigenstate(param_list):
    # create the parameter dictionary
    param_dict = {}
    for index in range (0, n_params):
        param_dict[theta_list[index].name()] = param_list[index]

    return find_eigenstate(param_dict)

print("Using scipy minimize, method: COBYLA:")

vqe_result_scipy = scipy.optimize.minimize(vqe, initial_parameters_list, method="COBYLA")

print ("The estimated ground state energy: ", vqe_result_scipy.fun)

eigenstate_scipy = vqe_eigenstate(vqe_result_scipy.x)
print ("Eigenstate: ", eigenstate_scipy)

This code produced the following output upon execution on a Google Colab notebook and it took about a minute to converge:

Using scipy minimize, method: COBYLA:
The estimated ground state energy:  -25.55078125
Eigenstate:  1110

On repeated execution, the Google Colab notebook also produced the following output:

Using scipy minimize, method: COBYLA:
The estimated ground state energy:  -29.4453125
Eigenstate:  0101

Note that the ground state energy may be slightly varying in your implementation. You may try other optimizers such as the SPSA for better results.

It is a good idea to see how the classical optimizer minimized the ground state energy, which corresponds to the general flow of the VQE algoritm. Execute the following code now:

# Plot the energies
plt.plot(avg_list)
plt.xlabel('Iteration')
plt.ylabel('<H>')
plt.grid(alpha=0.4,linestyle='--')
plt.show()

A graph simillar to the following figure can be expected:

VQE Optimizer.

Footnotes