Variational Quantum Eigensolvers with SNES
Training a Variational Quantum Eigensolver with evotorch and pennylane¶
This example demonstrates how you can train variational quantum eigensolvers (VQEs) using EvoTorch and PennyLane. To execute this example, you will need to install PennyLane with:
This example is based on PennyLane's VQE example and we encourage the reader to read that tutorial for a deeper understanding of the problem setting.
Basics of variational quantum eigensolvers¶
VQE is a recently introduced algorithm which has the objective to train a quantum computer to prepare the ground state of a given molecule [1]. Finding this ground state is a central problem in computational physics and quantum chemistry, as knowledge of this state then enables parameterisation of further simulations.
In brief, a VQE consists of 3 components:
- A parameterisable quantum circuit
Qwithpparameters, which prepares the ground state of the molecule. - A cost function
Cthat computes the energy of a given ground state, which we want to minimise. - A classical optimisation algorithm, which searches for the optimal parameter vector which minimises the energy.
Here we will be using variation quantum circuits for Q, which means that the circuits are parameterised by a vector of angles.
Limitations of backpropagation¶
The most natural solution to this problem would be to use backpropagation to train the parameters of the quantum circuit to minimise the cost. However, this is only possible in simulation, as backpropagation through a physical quantum circuit would require observing and reusing the state values of the circuit, which is impossible. Because of this, alternative approaches to differentiation of quantum circuits have been proposed. A particularly prominent approach is the 'parameter-shift rule' [2] which, effectively allows us to compute analytic gradients on the parameters by evaluating precise shifts of the parameters, one-by-one. This, of course, means that for one update of a circuit with p parameters, we must perform 2p circuit evaluations.
A recent paper [3] introduced the idea that evolutionary algorithms such as SNES and XNES may serve as an alternative approach to quantum circuit optimisation, where the number of circuit evaluation in each iteration is simply the population size k. In practice, this can yield similar performance to gradient-descent based methods in less circuit evaluations. In this example, we'll be doing exactly that: training the parameters of a VQE using SNES.
Setting up the Cost Function¶
First, we will need to define the molecular structure. We will study the H\(_2\) molecule, which means that we have:
import torch
symbols = ["H", "H"] #H2 molecule
coordinates = torch.tensor([0.0, 0.0, -0.6614, 0.0, 0.0, 0.6614]) # nuclear coordinates in atomic units
Next, we need to calculate the electronic Hamiltonian:
import pennylane as qml
H, qubits = qml.qchem.molecular_hamiltonian(symbols, coordinates)
print("Number of qubits = ", qubits)
print("The Hamiltonian is ", H)
We'll be using the default PennyLane quantum simulator
Now we're ready to set up the circuit Q and cost function C. Note that due to the simplicity of the molecule, we have only a single parameter, p=1 to optimise!
electrons = 2 # H2 has 2 electrons
hf = qml.qchem.hf_state(electrons, qubits) # giving the Hartree-Fock state
def circuit(param, wires):
# Setting up the circuit to optimise, which simply consists of preparing the basis state as the Hartree-Fock state
# And then applying a double excitation parameterised by the parameter
qml.BasisState(hf, wires=wires)
qml.DoubleExcitation(param, wires=[0, 1, 2, 3])
@qml.qnode(dev, diff_method=None) # Disabling gradients -- we don't need them
def cost_fn(param):
# Defining the cost function: simply apply the parameterised circuit and take the exponent of the Hamiltonian
circuit(param, wires=range(qubits))
return qml.expval(H)
Creating a EvoTorch Problem instance¶
With the cost function well-defined, we're ready to create a Problem instance. Note that we will repeat the steps above, so that the model has its own internal definition of the cost function. This will allow us to exploit parallelization with Ray.
from evotorch import Problem, Solution
from evotorch.algorithms import SNES
from evotorch.logging import StdOutLogger, PandasLogger
from typing import Optional
class VGEMin(Problem):
def __init__(self, num_actors: Optional[int] = None):
super().__init__(
objective_sense='min', # Minimise the objective
solution_length = 1, # There is only a single parameter to optimise
initial_bounds = (-1e-6, 1e-6), # Start the search very close to zero
num_actors = num_actors, # Number of ray actors
)
symbols = ["H", "H"]
coordinates = torch.tensor([0.0, 0.0, -0.6614, 0.0, 0.0, 0.6614])
self._H, self._qubits = qml.qchem.molecular_hamiltonian(symbols, coordinates)
electrons = 2 # H2 has 2 electrons
self._hf = qml.qchem.hf_state(electrons, qubits)
def _prepare(self):
# Prepare function called by actors once instantaited
dev = qml.device("default.qubit", wires=self._qubits)
# Inline definition of cost function allows us to easily decorate it as a quantum node
@qml.qnode(dev, diff_method = None, interface = 'torch')
def actor_cost_fn(param):
with torch.no_grad():
wires = range(self._qubits)
qml.BasisState(self._hf, wires=wires)
qml.DoubleExcitation(param[0], wires=[0, 1, 2, 3])
return qml.expval(self._H)
self._cost_fn = actor_cost_fn
def _evaluate(self, individual: Solution):
x = individual.access_values() # Get the decision values -- in this case a vector of length 1
cost = self._cost_fn(x) # Evaluate the decision values
individual.set_evals(cost) # Update the fitness
problem = VGEMin(num_actors = 4) # Instantiate the problem class
population = problem.generate_batch(5) # Generate a population to test things out
problem.evaluate(population) # If we've set everything up correctly we should get no error
print(f'Initial fitness values {population.access_evals()}')
Training the VQE with EvoTorch¶
Now we're ready to train. Simply create the searcher and logger:
searcher = SNES(problem, stdev_init=0.1) # stdev_init=0.1 used in [3]
logger = PandasLogger(searcher)
And train for 100 generations
Taking a look at the mean of the searcher we have:
Where the learned mean is close to the known approximate global optima of 0.208.
print(f'Cost of learned mean is {cost_fn(searcher.status["center"][0].numpy())} vs. approx global optima {cost_fn(0.208)}')
A more challenging example¶
With the basics down, we can now think about a more challenging problem. Instead of the simple H\(_2\) molecule, we will instead consider the water molecule H2O. A closely related experiment was performed in [3], meaning that this is a demonstration of a relatively 'state of the art' result for QVE training using EvoTorch.
Now we have the molecule configuration:
symbols = ["H", "O", "H"] #H2O molecule
coordinates = torch.tensor([-0.0399, -0.0038, 0.0, 1.5780, 0.8540, 0.0, 2.7909, -0.5159, 0.0]) # nuclear coordinates in atomic units
And as before we need the electronic Hamiltonian
H, qubits = qml.qchem.molecular_hamiltonian(
symbols,
coordinates,
charge=0,
mult=1,
basis="sto-3g",
active_electrons=4,
active_orbitals=4,
)
print("Number of qubits = ", qubits)
print("The Hamiltonian is ", H)
Making a new cost function, this time using UCCSD ansatz
from functools import partial
electrons = 10
orbitals = 7
core, active = qml.qchem.active_space(electrons, orbitals, active_electrons=4, active_orbitals=4)
singles, doubles = qml.qchem.excitations(len(active), qubits)
hf = qml.qchem.hf_state(
len(active),
qubits,
) # giving the Hartree-Fock state
# Map excitations to the wires the UCCSD circuit will act on
s_wires, d_wires = qml.qchem.excitations_to_wires(singles, doubles)
# Define the device
dev = qml.device('default.qubit', wires=qubits)
def circuit2(param, wires):
# Setting up the circuit to optimise, which simply consists of preparing the basis state as the Hartree-Fock state
# And then applying a UCCSD ansatz
qml.UCCSD(param, wires=wires, s_wires = s_wires, d_wires = d_wires, init_state = hf)
@qml.qnode(dev, diff_method=None) # Disabling gradients -- we don't need them
def cost_fn(param):
# Defining the cost function: simply apply the parameterised circuit and take the exponent of the Hamiltonian
circuit2(param, wires=range(qubits))
return qml.expval(H)
Putting this together in a problem definition:
# Defining a new problem:
class VGEH2O(Problem):
def __init__(self, num_actors: Optional[int]):
super().__init__(
objective_sense='min', # Minimise the objective
solution_length = 26, # There are 26 parameters to optimise
initial_bounds = (-1e-6, 1e-6), # Start the search very close to zero
num_actors = num_actors,
)
def _prepare(self):
symbols = ["H", "O", "H"] #H2O molecule
coordinates = torch.tensor([-0.0399, -0.0038, 0.0, 1.5780, 0.8540, 0.0, 2.7909, -0.5159, 0.0])
self._H, self._qubits = qml.qchem.molecular_hamiltonian(
symbols,
coordinates,
charge=0,
mult=1,
basis="sto-3g",
active_electrons=4,
active_orbitals=4,
)
electrons = 10
orbitals = 7
core, active = qml.qchem.active_space(electrons, orbitals, active_electrons=4, active_orbitals=4)
singles, doubles = qml.qchem.excitations(len(active), qubits)
self._hf = qml.qchem.hf_state(
len(active),
qubits,
)
self._s_wires, self._d_wires = qml.qchem.excitations_to_wires(singles, doubles)
# Prepare function called by actors once instantaited
dev = qml.device("default.qubit", wires=self._qubits)
# Inline definition of cost function allows us to easily decorate it as a quantum node
@qml.qnode(dev, diff_method = None, interface = 'torch')
def actor_cost_fn(param):
with torch.no_grad():
wires = range(self._qubits)
qml.UCCSD(param, wires=wires, s_wires = self._s_wires, d_wires = self._d_wires, init_state = self._hf)
return qml.expval(self._H)
self._cost_fn = actor_cost_fn
def _evaluate(self, individual: Solution):
x = individual.access_values() # Get the decision values
cost = self._cost_fn(x) # Evaluate the parameter vector
individual.set_evals(cost) # Update the fitness
problem = VGEH2O(num_actors = 4)
population = problem.generate_batch(10)
problem.evaluate(population)
And from there, all we need to do is train!
searcher = SNES(problem, stdev_init=0.1) # stdev_init=0.1 used in [3]
pandas_logger = PandasLogger(searcher)
stdout_logger = StdOutLogger(searcher)
# Run for 200 generations
searcher.run(200)
And visualize the progress:
References¶
[1] Peruzzo, Alberto, et al. "A variational eigenvalue solver on a photonic quantum processor." Nature communications 5.1 (2014): 1-7.
[2] Schuld, Maria, et al. "Evaluating analytic gradients on quantum hardware." Physical Review A 99.3 (2019): 032331.
[3] Anand, Abhinav, Matthias Degroote, and Alán Aspuru-Guzik. "Natural evolutionary strategies for variational quantum computation.". Machine Learning: Science and Technology 2.4 (2021): 045012