QHack QML Challenge Walkthrough: Variational Quantum Eigensolver

QHack 2021, the quantum machine learning (QML) hackathon, ran earlier this year from 17-26 February. A big portion of the event was the QML Challenge Leaderboard, where hundreds of teams raced to solve QML programming problems in order to claim hardware credits to help with their Open Hackathon projects. Solutions to the challenges were not provided during the event, so in this series of blog posts, we’ll walk you through the steps to solving them.

In this first post, we’ll take a look at the problems in the Variational Quantum Eigensolver category. The category contained three problems, the difficulty of which was indicated by a point value (100, 200, and 500). The full set of problems is available at the QHack2021 Github repository.

Background: the VQE

The variational quantum eigensolver (VQE) is an algorithm that, in its most basic form, can be used to find the ground state energy of a quantum system. The VQE algorithm works by parametrizing the space of possible quantum states, and then solving the optimization problem

where α are the parameters, and H is the Hamiltonian of the system. This cost function is the expectation value, or energy of the system given a particular \alpha. The operation U(\alpha) is a special type of quantum circuit called an ansatz. It is specially created in that we expect there is an \tilde{\alpha} such that

where |\psi_g \rangle is the ground state.

The VQE algorithm consists, at a high level, of the following steps:

1. Choose a suitable ansatz circuit U(\alpha)
2. Choose a starting set of parameters \alpha
3. Apply U(\alpha) and measure the output state
4. Use measurement results to compute the numerical value of \langle 0 \vert U^\dagger(\alpha)HU(\alpha)\vert 0 \rangle (the energy)
5. Run an optimization routine to find a new \alpha that should produce a state nearer the ground state
6. Repeat steps 3-5 until the optimizer converges to a minimum value, or the number of iterations has exceeded a specified maximum

Optimization Orchestrator (100-point question)

Percentage solved: \approx 89% (top 100 teams) / \approx 44% (all teams with submissions)

The first problem served as a warmup: given a Hamiltonian and a variational ansatz, write PennyLane code for the classical optimization component (steps 3-6 above) that will find the ground state energy. The provided ansatz, variational_ansatz, consists of layers of arbitrary rotations followed by a ring pattern of CNOTs:

The coding task was to complete the missing parts of the following function:

def run_vqe(H):
    """Runs the variational quantum eigensolver on the problem Hamiltonian using the
    variational ansatz specified above.

    Fill in the missing parts between the # QHACK # markers below to run the VQE.
    Args:
        H (qml.Hamiltonian): The input Hamiltonian
    Returns:
        The ground state energy of the Hamiltonian.
    """
    # Initialize parameters
    num_qubits = len(H.wires)
    num_param_sets = (2 ** num_qubits) - 1
    params = np.random.uniform(low=-np.pi / 2, high=np.pi / 2, size=(num_param_sets, 3))

    energy = 0

    # QHACK #

    # QHACK #

    # Return the ground state energy
    return energy

There are many different ways to solve this problem; the solution involves initializing a device, building a cost function based on the expectation value of the Hamiltonian, and then running the optimization:

# QHACK #

# Initialize a quantum device
dev = qml.device('default.qubit', wires=num_qubits)

# Set up a cost function based on a provided ansatz
cost_fn = qml.ExpvalCost(variational_ansatz, H, dev)

# Set up an optimizer
opt = qml.GradientDescentOptimizer(stepsize=0.1)

# Run the VQE by iterating over many steps of the optimizer
max_iterations = 500

for n in range(max_iterations):
    params, energy = opt.step_and_cost(cost_fn, params)

# QHACK #

Ansatz Artistry (200-point question)

Percentage solved: \approx 56% (top 100 teams) / \approx 19% (all teams with submissions)

The 100 point question addressed steps 3-6 of the VQE steps, while the requirements of steps 1 and 2 were provided. In the 200 point problem, participants had to handle these steps themselves: they needed to design a specialized ansatz to create quantum states of a very particular form. Specifically, the n -qubit Hamiltonians used for this problem had eigenstates of the form

where all the \alpha_i are real numbers.

All that was provided for this problem was a function header, and a single hint:

“A good way to start is by determining the number of parameters you will need, their properties, and what quantum operations produce them. Work iteratively, starting with the n=2 case, then move to n=2, to see how the ansatz generalizes.”

Let’s address the hint step by step. First, we know there are n real-valued coefficients, and for the quantum state to be normalized, it must be the case that

Looking at some small cases of this, you should notice something familiar: the parameters satisfy the same properties as spherical coordinates! This means we can parametrize the state using generalized, n -dimensional spherical coordinates, of which there are n−1 . We will thus require n−1 single-parameter rotations in our quantum circuit; since the coefficients are all real, we’ll use y rotations.

The next step is to figure out how to combine y rotations and other circuit building blocks to produce states of the desired form. This is where the other aspect of the hint comes in: it is useful here to work iteratively. Consider the two-qubit case: using a single rotation, we need to produce a state of the form

where in the second expression, we’ve rewritten the coefficients in terms of a single parameter, \theta_1.

The ansatz acts on the |00\rangle state; to get started, we’ll have to apply an X gate in order to obtain a |1\rangle in the superposition. Let’s apply it to the first qubit:

Next, we add the rotation. In order to match the state above, we’ll apply it to the second qubit:

Note that the factor of 1/2 in the angle doesn’t really matter — the effective parametrization is still the same, just the variable is rescaled.

This looks closer to our desired state, however the second ket is wrong! It can have only a single qubit in state |1\rangle . Luckily, this is easily corrected by a CNOT gate from the second qubit to the first:

And there we have it! The circuit looks like the following:

Next, we move up to the 3-qubit case. You’ll notice that the 2-qubit case actually gets us most of the way there - if we apply the same operations as above to the first two of a set of three qubits, we’ll obtain

We need to do two things here:

• apply a second rotation with a different parameter, \theta_2,
• correct any spurious 1s that occur in the process.

To handle the first aspect, note that if we apply a rotation to any one qubit, we’ll get a linear combination of four kets in the process. To produc the combination of three kets, we thus want to apply a rotation in only one case. We can accomplish this by using a rotation controlled on the value of the second qubit so that only the second term in the superposition is affected:

We must then fix the middle |1\rangle in the third ket — we can do this in the same manner as the two-qubit case above: by applying a CNOT from the third qubit to the second to obtain

Note that these are exactly the generalized spherical coordinates for a 4-sphere, angles and all.

To generalize to higher numbers of qubits, we can take the same approach: apply a controlled y — rotation to the last qubit, and then apply a CNOT from the last qubit to the second last in order to cancel the extra |1\rangle. The circuits thus have the following form:

Implemented in PennyLane, they should look like this:

def variational_ansatz(params, wires):
    """ The variational ansatz circuit.

    Fill in the details of your ansatz between the # QHACK # comment markers. Your
    ansatz should produce an n-qubit state of the form

        a_0 |10...0> + a_1 |01..0> + ... + a_{n-2} |00...10> + a_{n-1} |00...01>

    where {a_i} are real-valued coefficients.

    Args:
         params (np.array): The variational parameters.
         wires (qml.Wires): The device wires that this circuit will run on.
    """
    # QHACK #
    qml.PauliX(wires=wires[0])

    qml.RY(params[0], wires=wires[1])
    qml.CNOT(wires=[wires[1], wires[0]])

    for wire_idx in range(1, len(params)):
        qml.CRY(params[wire_idx], wires=[wires[wire_idx], wires[wire_idx + 1]])
        qml.CNOT(wires=[wires[wire_idx + 1], wires[wire_idx]])
    # QHACK #

Moving on up (500-point question)

Percentage solved: \approx 30% (top 100 teams) / \approx 10% (all teams with submissions)

The final question posed quite the challenge for participants: given a Hamiltonian, use VQE methods to find not only the ground state energy, but also the energies of the first two excited states. Finding the excited state energies is an active area of research, and there are a couple different methods one could use. We solved this problem using the subspace-search VQE (SSVQE) algorithm. Surprisingly, it required adding only a few additional lines to the VQE cost function.

The key idea in the VQE is to find an ansatz circuit that, given just the right set of parameters, transforms the |0{\cdots}0\rangle state into the ground state. Now, this transformation is unitary, and unitary transformations preserve the angles between vectors. So everything that is orthogonal to |0{\cdots}0\rangle before the transformation is still orthogonal to it afterwards. This leads to an interesting question: can we find a unitary that simultaneously transforms |0{\cdots}0\rangle to the ground state, and another basis state to the first excited state? Can we find a unitary that does this for the first two excited states? What about more, is there a unitary transformation that sends a subset of the computational basis states to the complete set of eigenvectors of the Hamiltonian?

This is the task accomplished using SSVQE. In particular, in our solution we use a variant called the weighted SSVQE that enables you to simultaneously find multiple eigenstates. Rather than considering a cost function that involves only one term, the cost function involves computing the weighted expectation values for multiple candidate eigenstates. That is, our optimization problem becomes:

where the \{\vert\varphi_j \rangle\} are three orthogonal starting states (we simply use three computational basis states), and the w_j are different weights that ensure we find the right eigenstates. The code below implements exactly this algorithm - it uses a very generic ansatz, and the weight vector is simply a list of integers.

def find_excited_states(H):
    """
    Fill in the missing parts between the # QHACK # markers below. Implement
    a variational method that can find the three lowest energies of the provided
    Hamiltonian.

    Args:
        H (qml.Hamiltonian): The input Hamiltonian

    Returns:
        The lowest three eigenenergies of the Hamiltonian as a comma-separated string,
        sorted from smallest to largest.
    """

    energies = np.zeros(3)

    # QHACK #

    # Get the number of qubits in the Hamiltonian
    num_qubits = len(H.wires)

    # Initialize the device
    dev = qml.device('default.qubit', wires=num_qubits)

    # Implements the weighted subspace search VQE [1810.09434]

    # This ansatz works well enough, it's not the same as the one in the paper though
    def ansatz(params, wires, state_idx=0):
        # Need to prepare a different orthogonal state each time
        qml.PauliX(wires=wires[state_idx])
        qml.templates.StronglyEntanglingLayers(params, wires=wires)

    # The total cost is a sum of expectation values for different orthogonal
    # starting states. Essentially the ansatz turns into something that can
    # simultaneously map the original starting states to the eigenstates.
    single_cost = qml.ExpvalCost(ansatz, H, dev)

    # Weight vector - weights the different eigenstates in the cost function so
    # that it's the lowest ones that are found
    w = np.arange(num_qubits, 0, -1)

    # The full cost - computes single_cost for each starting state
    def total_cost(params):
        cost = 0
        for state_idx in range(num_qubits):
            cost += w[state_idx] * single_cost(params, state_idx=state_idx)
        return cost

   # Set up a cost function and optimizer, and run the SSVQE
    opt = qml.GradientDescentOptimizer(stepsize=0.05)
    max_iterations = 200
    costs = []

    # Initial parameters for the ansatz
    num_layers = 8
    params = np.random.uniform(low=0, high=2*np.pi, size=(num_layers, num_qubits, 3))

    # Optimize!
    for _ in range(max_iterations):
        params = opt.step(total_cost, params)

    # After optimization, get the energies for the original starting states
    for state_idx in range(3):
        energies[state_idx] = single_cost(params, state_idx=state_idx)

    # QHACK #

    return ",".join([str(E) for E in energies])

Some sample Hamiltonians are provided in the problem template. Try playing around with the weight vector, ansatz structure, and the optimizer. Can you tweak this in order to find all the eigenstates of these Hamiltonians?

Concluding thoughts

We encourage you to try out these solutions yourselves, and explore the many options available. These solutions are far from unique, and there are plenty of different methods to explore, especially for the 500-point problem. If you find a particularly interesting or optimized solution, we’d love to see it!