Introducing the Variational Quantum Eigensolver

In the previous blogs, we have learned how to map the max-cut and the Traveling Salesman Problem to the Ising Model. We also solved the problems by computing the eigenvalues of the Hamiltonian matrix and finding the eigenstate of the corresponding minimum eigenvalue. We are able to do that because our problem size is small. For example, for the Traveling Salesman Problem, a problem size of n=3 cities needs 3^2 = 9 qubits. This corresponds to a matrix size of (2^9)\cdot (2^9)=512\times 512 = 262144. If n=4, the matrix size becomes 4294967296 elements! There is another way to tackle this complexity. It’s a heuristic called the Variational Quantum Eigensolver.

VQE Algorithm

The Variational Quantum Eigensolver algorithm looks like this:

1. Using a set of parameters \theta, generate a trial state |\psi\rangle using the quantum computer.
2. Compute the expectation value


3. If the result above is less than the current minimum value, update the minimum value.
4. Use a classical optimizer to drive the generation of the \theta parameters. For example, the LBFGS as the classic optimizer can be used.
5. Repeat the process until you converge to the minimum value.

Step 4 is standard, so we will just need to understand how qubits are generated and how the expectation value is computed. Finally, we will use the VQE to solve an instance of the TSP.

Generating Trial States

Given the qubit |0\rangle, we can generate a superposition of |0\rangle and |1\rangle by performing a rotation along the Y axis. This Y rotation is defined by

R_y(\alpha)=\exp(-i\alpha Y/2) = \left( \begin{array}{cc}  \cos(\frac{\alpha}{2}) & -\sin(\frac{\alpha}{2})\\  \sin(\frac{\alpha}{2}) & \cos(\frac{\alpha}{2})  \end{array}\right)

where Y is the Pauli Y operator given by the matrix

\sigma_y = \left(  \begin{array}{cc}  0 & -i\\  i & 0  \end{array}  \right)

We can also entangle qubits using pairwise Control-Z gates. The resulting circuit looks like the below, with the one in the RED box applied d times (d is called the depth.) In QISKit, this is implemented via the Ry variational form.

Computing the Expectation Value

Let |\psi\rangle a state vector and H be a Hermitian Operator, then the expectation value of H on the state |\psi\rangle is given by


Let \{\lambda_0,\ldots,\lambda_{n-1}\} be the eigenvalues of H and \{|\phi_0\rangle,\ldots|\phi_{n-1}\rangle\} be the corresponding eigenvectors. Then we can write the above as

\begin{array}{rl}  \displaystyle \langle\psi|H\psi\rangle &= \displaystyle \langle\psi|\sum_{i=0}^{n-1} \lambda_i |\phi_i\rangle\langle\phi_i|\psi\rangle\\  &= \displaystyle \sum_{i=0}^{n-1} \lambda_i \big|\langle\psi|\phi\rangle\big|^2  \end{array}

The above expression is similar to the arithmetic mean where \big|\langle\psi|\phi\rangle\big|^2 is the probability of the eigenvalue \lambda_i occurring.

Let’s see how we can use this in computing the expectation value of H in the Ising Model.

Suppose we have the following circuit,

Using QISKit, we can simulate this circuit:


backend_sim = Aer.get_backend('qasm_simulator')
shots = 1000
job_sim = execute(qc, backend=backend_sim, shots=shots)
result_sim = job_sim.result()
answer = result_sim.get_counts()

for i in answer:

The simulation above gives us the probability of a given eigenstate to occur. To get the expectation value, we just need to multiply the eigenvalue with the probability of it to occur and sum it for all eigenvalues. In the Ising Model, we can easily read-off the eigenvalue of a given eigenstate. Since the eigenvalues/eigenvector pairs of Z are given:

Z |0\rangle = |0\rangle


Z|1\rangle = -|1\rangle

Then eigenvalue of a state in the Ising Model is equal to

(-1)^\text{Number of occurrences of 1}

For example, for the state 1100, we have \lambda = -1^2 = 1.

Using this formula, we can compute for the expectation value of the above circuit:

\begin{array}{ccrr}  \mathrm{State} & \mathrm{Probability} & \mathrm{Eigenvalue} & \mathrm{Probability}\times\mathrm{Eigenvalue} \\\hline  0010 & 0.135 & -1 & -0.135 \\  0000 & 0.720 & 1 & 0.720 \\  0001 & 0.123 & -1 & -0.123 \\  0011 & 0.022 & 1 & 0.022 \\\hline  \mathrm{Expectation Value} & & & 0.484  \end{array}

Applying to the Traveling Salesman Problem

Now that we understand how the VQE works, we can just apply this to the Traveling Salesman Problem. We assume you have read this blog to make sense of the code below (lifted directly from the QISKit tutorial):

seed = 10598

spsa = SPSA(max_trials=300)
ry = RY(qubitOp.num_qubits, depth=5, entanglement='linear')
vqe = VQE(qubitOp, ry, spsa, 'matrix')

backend = Aer.get_backend('statevector_simulator')
quantum_instance = QuantumInstance(backend=backend, seed=seed)

result =

print('energy:', result['energy'])
print('time:', result['eval_time'])
x = tsp.sample_most_likely(result['eigvecs'][0])
print('feasible:', tsp.tsp_feasible(x))
z = tsp.get_tsp_solution(x)
print('solution:', z)
print('solution objective:', tsp.tsp_value(z, tspdata.w))
draw_tsp_solution(G, z, colors, pos)

This gives the following result:

energy: -597002.6480513655
time: 19.69760489463806
feasible: True
solution: [1, 2, 0]
solution objective: 223.0


This completes the series on how Quantum Computing can be used to tackle NP-Complete problems. More information can be found from here and here.


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s