IBM Quantum Challenge 2021 — Part 3
Last but not least, we have the most recent major advances in quantum computing, with the Variational Quantum Eigensolver. This algorithm is a realization of physicist Richard Feynman’s dream of simulating large and complex molecules. As with our journey in the IBM Quantum Challenge, we are on our fourth day of the challenge, with two full days to get a place on the leaderboard.
Now with Quantum Error Correction and Transmon Qubits completed, we have one more exercise more to go. VQE was one of those things, along QML, that I’ve been meaning to dive into for so long. And when I said meaning, I mean procrastinating, partly because the Qiskit textbook looks so complicated, but also because I never really liked application-specific algorithms. Nevertheless, I kept my fingers crossed and proceeded to watch the 6 hours long lecture of last year’s Qiskit Summer School. And boy, did I change my mind quickly. So here we go, VQE 101, explained in less than 6 hours.
Motivating the Molecular Hamiltonian
The world around us is made up of atoms, and atoms made up of molecules. If we want to manipulate and use these molecules in interesting ways, we have to study and understand them. Now one of the important properties of molecular systems is its Hamiltonian, which is the system’s total energy. Below is the formula of the molecular Hamiltonian.
I think it’s important not to be overwhelmed by all the novel symbols and focus on the intuitions behind all the variables and how they relate to one another. I don’t know about you but when I look at the kinetic energy terms, they just remind me so much of the ½ mv² formula from Physics class. From previous exposure, I learned that nabla (the inverted triangle) is a vector of the first derivative. Combined with the fact that the “i” and “A” represent electrons and nuclear coordinates, I have very strong reason to believe the nabla is the velocity vector. Similarly, with the potential energy term, I think it’s quite nice that it is Coulomb’s law in action where E = kqQ/r.
In summary, the Hamiltonian of a molecule takes into account the kinetic energy that comes from the vibration of the nucleus and electrons, and the potential energy of the nucleus and electrons interact with themselves and with one another.
The molecular Hamiltonian is a very rich quantity because it contains so much information about the interactions of the constituents of the molecules. But exactly why is this important? The answer is the same with everything in life. Not all Hamiltonians are created equal. And in the world of quantum noise and uncertainty fluctuations, the ground-state Hamiltonian comes out on top, or rather the bottom.
The lowest energy is the state where the molecules are most stable in. A chemical reaction happens because it results in its constituent having lower energy. So, to predict if a chemical reaction will occur, we check if the final energy is lower than the initial energy. In technical terms, we want to find the lowest eigenvalues of the Hamiltonian, which itself is a square matrix. For further reading on this eigen-jargon, check out this article by Quantum Untangled.
However, finding this energy is hard, so we will simplify the problem by considering the nuclear term and the electronic term separately.
Yet the process of determining the ground state is still very much a treacherous one. But at least for now, we can stay calm knowing that quantum computers and VQE have it under control.
From Molecular Hamiltonian to Qubits
The first step to any quantum algorithm is encoding the inputs as qubits. For example, with Simon, or Shor, we have the input from 0 to 2^n — 1 represented as binary superpositions by taking the Hadamard of all the inputs. Similarly, for VQE, recall the Born-Oppenheimer approximation, we can now only consider the electrons self-interaction, which in itself is a reduced form of the general molecular Hamiltonian.
More explicitly, we set the nuclear KE term to 0, and treat the Electron-Nuclear Coulomb interaction as an energy shift in the zero point fluctuation.
The next step is rewriting the variables in terms of quantum mechanical vectors and operators. We end up with the below equation.
The parts that we are interested in here, are the a’s, and a-daggers, which stand for creation and annihilation operators respectively. However, the full technicality of this equation is beyond the scope of learning to use VQE. Understanding the steps of the derivation is not too difficult per se, so here’s the link to the lecture that I got these equations.
If you have time on your hand, I would recommend checking it out. Again, it’s not too difficult, it’s just long.
There are a few methods to map the Hamiltonian inputs into the qubits input. The goal of this is to rewrite the creation and annihilation operator to Pauli spin matrices, the familiar I, X, Y, Z we use in quantum circuits!
In the challenge, we are provided with 3 functions:
- The Jordan-Wigner map uses O(N) Pauli-Z matrix for N creation or annihilation operators.
- The Parity map also scales O(N), and it uses fewer qubits than the Jordan-Wigner map.
- The Bravyi-Kitaev map is O(log(N)), however, the advantage is only seen when N is sufficiently large.
The Ansatz
After the mapping to qubits is done, the molecular Hamiltonian is now a Hermitian matrix that has below form.
The expectation value on an arbitrary state |ψ ⟩ is
By definition of eigenvalues,
We can derive
Or using an arbitrary state, by the variational principle,
The physical implementation of ansatz is called the variational form, which is a parameterized circuit that initializes a starting point.
Ideally, we want the variational form to be able to generate any state |ψ ⟩ when |ψ ⟩ is an element of the complex vector space of n qubits but using as few parameters as possible. For n = 1 and n=2, we have some simple examples of variational forms.
So to reiterate, our objective is to find the ground state eigenvalue λ, and to find this, we want to minimize λ(θ) as much as possible, and this is where optimization algorithms come in.
Parameter Optimization
The goal of the optimization algorithms is to find the set of parameters that minimizes the eigenvalues of our Hamiltonian Matrix. For VQE, I learned that general gradient descent is sub-optimal because it is expensive and can run into problems with local minima. In the challenge, we are provided with a few options, each with its pros and cons.
SPSA approximates the gradient and updates the gradient randomly. SLSQP & COBYLA can be used when noise is minimal. COBYLA is especially fast because it operates independently of the dataset size. Finally, L_BFGS_B uses the second derivative in limited memory space.
To recap on everything, here’s a flowchart explanation of VQE. After the molecular Hamiltonian is mapped into qubits and turned into a Hermitian matrix.
- The ansatz initialized qubits to a general state.
- Run and measure the circuit multiple times, to find the expectation value, which we can use to determine the eigenvalues.
- Use a classical optimization algorithm to update the parameters to lower the minimize the eigenvalues.
- Repeat step 1–3
The cost of the solution is defined as the number of CNOT gates used in the entire circuit. This includes the initialization circuit and the ansatz. Now, with this information in mind, we can start going through the challenge notebook.
Cost 172 Solution
So after looking through the example of H2, I understand that the task is to build an ansatz for which we are provided with 4 options: TwoLocal, UCCSD, PUCCD, and SUCCD. I’m just going to explain the one I use for my solutions, which in this case is UCCSD. The motivation for UCCSD comes from coupled-cluster theory which is among the most accurate classical methods for many-body simulation. In our case, the many-body we are working with are the electrons orbiting around our Li and H atoms. This is the paper I’ve read on UCCSD.
The author goes into great depth explaining all the technicalities behind UCCSD. For this challenge, we only need to know that UCCSD is a very robust ansatz, which means that it can represent a wide range of chemical states. However, we trade off this representation for its relatively big size. On my first implementation of UCCSD with only the “two_qubits” options set as True, I discovered that my circuit is too big that the grader doesn’t even accept it. After some decomposition of the ansatz, I found out that my intuition was correct.
This is a technical problem I was not ready for. So I once again went back to the problem statement and tried to implement the hints that Qiskit provided.
Before getting into implementing these tips, we need to initialize our molecule.
molecule = ‘Li 0.0 0.0 0.0; H 0.0 0.0 1.5474’driver = PySCFDriver(atom=molecule)qmolecule = driver.run()Matrix(np.round(qmolecule.one_body_integrals, 10))
Now, we can look at its one-body integrals that will help us with the first tip about freezing core electrons.
The molecular orbitals that are going to be frozen are determined by the transformer when the “freeze core” parameter is set to True, which is the default. Therefore, we only need to worry about which orbitals to further remove. Following the method the lecturer uses in this video, we know that we can remove orbitals 3 and 4. Now that we know which orbitals to remove, we do so while initializing the electronic structure problem object. We do so with the following lines of code.
problem = ElectronicStructureProblem(driver, [FreezeCoreTransformer(freeze_core = True, remove_orbitals=[3,4])]) # Freeze orbitals 3 and 4second_q_ops = problem.second_q_ops() # Second-quantized operatorsmain_op = second_q_ops[0] # Hamiltonian
To reduce the number of qubits further, we use a special mapper that allows us to get rid of two more qubits. On top of that, we use Z2 symmetries to reduce the number of qubits even more.
mapper = ParityMapper()converter = QubitConverter(mapper=mapper, two_qubit_reduction=True, z2symmetry_reduction=[1,1])num_particles = (problem.molecule_data_transformed.num_alpha,problem.molecule_data_transformed.num_beta)qubit_op = converter.convert(main_op, num_particles=num_particles)
We see that the number of qubits in our circuit is just 4, which allows us to connect all qubits using fewer CNOT gates.
Now, we need to declare the backend and optimizer. I chose to use the “COBYLA” optimizer because I found it to be the best for my specific ansatz. It was also the one that got me the lowest error along with the lowest cost in the tutorial, so I thought it would also work in this case. Moreover, I set “maxiter=5000” because I found that the optimizer halted around that number even if I had the maximum set higher.
backend = Aer.get_backend(‘statevector_simulator’)optimizer = COBYLA(maxiter=5000)
We can now go ahead and run the VQE algorithm. We can use the piece of code that was given to us in the tutorial since nothing important is hardcoded.
try:
initial_point = [0.01] * len(ansatz.ordered_parameters)
except:
initial_point = [0.01] * ansatz.num_parametersalgorithm = VQE(ansatz, optimizer=optimizer, quantum_instance=backend, callback=callback,initial_point=initial_point)result = algorithm.compute_minimum_eigenvalue(qubit_op)print(“Calculated eigenvalue:”, result.eigenvalue)
print(“Optimizer evaluations:”, result.optimizer_evals)
Using the following Numpy code we can compare our result with the expected value.
solver = NumPyMinimumEigensolverFactory()calc = GroundStateEigensolver(converter, solver)result_exact = calc.solve(problem)exact_energy = np.real(result_exact.eigenenergies[0])print("Exact ground energy:", exact_energy)
Here, we can see that 𝛿𝐸 ≈ 2.339 < 4 mHa, meaning that we got the desired chemical accuracy.
Now that I officially have a working solution under my belt, I went on to look for ways to optimize my circuit even more.
Cost 3 Solution
For starters, we get the “HartreeFock” state to use as our initialization circuit. Notice that it is different from the initial state used in the tutorial part for the H2 molecule.
num_spin_orbitals = 2 * problem.molecule_data_transformed.num_molecular_orbitalsinit_state = HartreeFock(num_spin_orbitals, num_particles, converter)init_state.draw(‘mpl’)
Next, we need to build our ansatz. We are going to exploit the fact that we managed to bring down the number of qubits to 4 by using the “TwoLocal” circuit and setting “entanglement” to “linear”. This makes the total number of CNOT gates three, which makes a ladder from qubit 0 through qubit 3. At first, I was using the “rotation_blocks” array given in the tutorial part, which just consisted of two layers: one Rz followed by one Ry. However, simulating with this ansatz caused the optimizer to halt after very few iterations, resulting in a really inaccurate result. Thus, I decided to add another layer of Rz and Ry. Also, notice that “skip_final_rotation_layer” is set to “False”; I did this to have rotation blocks on both sides of the CNOT stair and thus allow for better optimization.
rotation_blocks = [‘rz’, ‘rz’, ‘ry’, ‘ry’]entanglement_blocks = ‘cx’entanglement = ‘linear’repetitions = 1ansatz = TwoLocal(qubit_op.num_qubits, rotation_blocks, entanglement_blocks, reps=repetitions,entanglement=entanglement, skip_final_rotation_layer=False)ansatz.compose(init_state, front=True, inplace=True)getScore(ansatz)ansatz.draw(‘mpl’)
We follow the same optimization step as the above solution, and we got the lowest cost possible of 3 CNOT gates. We’ve solved exercise 5 with the lowest cost possible!
Cost 1 Solution
So something unexpected happens during the course of the challenge. The best cost of exercise 5 keeps going down and down. At the beginning of the first day, it was 225, and at the end of the day saw a significant drop to 12. Some time on the third day, something ridiculous happened: the cost was 1. Even though there are some significant red flags for solutions with such a solution, I, amongst many people, begin to pursue the “Holy One”. However, on the fourth day of the challenge, there was a lot of controversy within the community about the plausibility of cost 1 and 2 solutions. After reviewing all of such solutions, IBM decided to make the minimum cost 3, which we have demonstrated above.
Using all the same circuit settings from above, i.e freeze-core, Z2 Symmetry, and two-qubits reduction let us begin our journey towards the forbidden cost 1 solution. This is in fact quite a simple solution, yet somehow it’s not as elegant as we have seen with cost 3 from above.
To do this, we need to use the “Custom” ansatz option. This is because the smallest “TwoLocal” circuit is size 3. We will need to initialize 18 circuit parameters, which we will use on the four U gates at the beginning, and 2 U gates after our single CNOT.
# Create parameters
a = Parameter('a')
b = Parameter('b')
c = Parameter('c')
d = Parameter('d')
e = Parameter('e')
f = Parameter('f')
g = Parameter('g')
h = Parameter('h')
i = Parameter('i')
j = Parameter('j')
k = Parameter('k')
l = Parameter('l')
m = Parameter('m')
n = Parameter('n')
o = Parameter('o')
p = Parameter('p')
q = Parameter('q')
r = Parameter('r')# Pre-CX U gates
qc.u(a, b, c, 0)
qc.u(d, e, f, 1)
qc.u(g, h, i, 2)
qc.u(j, k, l, 3)# CX
qc.cx(0, 1)# Post-CX U gates
qc.u(m, n, o, 0)
qc.u(p, q, r, 1)
Because we applied the U gates before the CNOT on qubits 0 and 1, we have basically nullified the effect of the X gate. This is because any states that can be represented with XU can be represented with just U. The control and target bit of the CNOT is also quite trivial. This is because the circuit is invariant with respect to CNOT transformation.
After we created our ansatz, we can train it using any of the optimizers provided. I myself use 5000 iterations of “L_BFGS_B” and here is the result.
Recall from above, the optimal value with freeze core is supposed to be -1.0887. We are only 0.0184 away, yet this is still outside the threshold of the grader of 4 mHa. I tried using different optimizers and letting each optimizer run for longer, but none of these fixes seems to help with the convergence. I looked on Slack to see what other cost-1 solvers have to say. It didn’t help at all. I even tried out all the combinations of circuit settings, yet none of the settings seems to work. The more I experimented and failed, the more I realized something fishy is going on. But in the end, I found a way to get my value to converge. I started to think about the molecular Hamiltonian. I reasoned that if the atoms were a little closer to one another, the potential energy of the electron’s repulsion would go up and the ground state energy would go down. The way to “move the atoms closer”, was to modify the driver’s distance of the electrons. With trial and error, I found the intermolecular distance from 1.468, and I got a good eigenvalue within the accepted range of the solution.
molecule = ‘H 0.0 0.0 0.0; Li 0.0 0.0 1.468’This modification yields the ground state eigenvalue of -1.08871051
Although as I found out later on, what I did was illegal, and subsequently I had to update my solution to a cost 3 ansatz (as previously shown). Oh well…
And that’s it, ladies and gentlemen, that is IBM Quantum Challenge 2021 in a nutshell for you. In the last week, we’ve taken a tour through the history of quantum computing landmarks. We started in Part 1 with the Toffoli gate and Shor’s algorithm. We then visited Quantum Error Correction and Transmon Qubits in Part 2. And today, we concluded our journey with the Variational Quantum Eigensolver.
The editors of the IQC-21 series sincerely hope that you have enjoyed the discussion of the ideas behind each quantum breakthroughs and the walkthrough of the solutions (legal or otherwise). We look forward to bringing you more content in the future. If you have any comments on the series or any ideas of what you would like us to cover in the future, feel free to leave a comment or join us at our Quantum Untangled discord server.
We wish you all the best in your future quantum computing endeavors.
Extra special thanks to Emilio Peláez for his collaboration on IQC-21, and for his help co-authoring the IQC-21 series.
