Variational Quantum Eigensolver¶
The Variational Quantum Eigensolver (VQE) algorithm, originally proposed in this publication, is a hybrid algorithm that can make use of both quantum and classical resources.
Let’s take a look at how we can use CUDA-Q’s built-in vqe module to run our own custom VQE routines!
Given a parameterized quantum kernel, a system spin Hamiltonian, and one of CUDA-Q’s optimizers,
cudaq.vqe will find and return the optimal set of parameters that minimize the energy, <Z>, of the system.
The code block below represents the contents of a file titled simple_vqe.py.
import cudaq
from cudaq import spin
from typing import List
# We begin by defining the spin Hamiltonian for the system that we are working
# with. This is achieved through the use of `cudaq.SpinOperator`'s, which allow
# for the convenient creation of complex Hamiltonians out of Pauli spin operators.
hamiltonian = 5.907 - 2.1433 * spin.x(0) * spin.x(1) - 2.1433 * spin.y(
0) * spin.y(1) + .21829 * spin.z(0) - 6.125 * spin.z(1)
# Next, using the `cudaq.Kernel`, we define the variational quantum circuit
# that we'd like to use as an ansatz.
# Create a kernel that takes a list of floats as a function argument.
@cudaq.kernel
def kernel(angles: List[float]):
# Allocate 2 qubits.
qubits = cudaq.qvector(2)
x(qubits[0])
# Apply an `ry` gate that is parameterized by the first value
# of our `angles`.
ry(angles[0], qubits[1])
x.ctrl(qubits[1], qubits[0])
# Note: the kernel must not contain measurement instructions.
# The last thing we need is to pick an optimizer from the suite of `cudaq.optimizers`.
# We can optionally tune this optimizer through its initial parameters, iterations,
# optimization bounds, etc. before passing it to `cudaq.vqe`.
optimizer = cudaq.optimizers.COBYLA()
# optimizer.max_iterations = ...
# optimizer...
# Finally, we pass all of that into `cudaq.vqe`, and it will automatically run our
# optimization loop, returning a tuple of the minimized eigenvalue of our `spin_operator`
# and the list of optimal variational parameters.
energy, parameter = cudaq.vqe(
kernel=kernel,
spin_operator=hamiltonian,
optimizer=optimizer,
# list of parameters has length of 1:
parameter_count=1)
print(f"\nminimized <H> = {round(energy,16)}")
print(f"optimal theta = {round(parameter[0],16)}")
// Compile and run with:
// ```
// nvq++ vqe_h2.cpp -o vqe.x && ./vqe.x
// ```
#include <cudaq.h>
#include <cudaq/algorithm.h>
#include <cudaq/builder.h>
#include <cudaq/gradients.h>
#include <cudaq/optimizers.h>
// Here we build up a CUDA-Q kernel with N layers and each
// layer containing an arrangement of random SO(4) rotations. The algorithm
// leverages the CUDA-Q VQE support to compute the ground state of the
// Hydrogen atom.
// The SO4 random entangler written as a CUDA-Q kernel free function
// since this is a pure-device quantum kernel
__qpu__ void so4(cudaq::qubit &q, cudaq::qubit &r,
const std::vector<double> &thetas) {
ry(thetas[0], q);
ry(thetas[1], r);
h(r);
x<cudaq::ctrl>(q, r);
h(r);
ry(thetas[2], q);
ry(thetas[3], r);
h(r);
x<cudaq::ctrl>(q, r);
h(r);
ry(thetas[4], q);
ry(thetas[5], r);
h(r);
x<cudaq::ctrl>(q, r);
h(r);
}
// The SO4 fabric CUDA-Q kernel. Keeps track of simple
// arithmetic class members controlling the number of qubits and
// entangling layers.
struct so4_fabric {
void operator()(std::vector<double> params, int n_qubits,
int n_layers) __qpu__ {
cudaq::qvector q(n_qubits);
x(q[0]);
x(q[2]);
const int block_size = 2;
int counter = 0;
for (int i = 0; i < n_layers; i++) {
// first layer of so4 blocks (even)
for (int k = 0; k < n_qubits; k += 2) {
auto subq = q.slice(k, block_size);
auto so4_params = cudaq::slice_vector(params, counter, 6);
so4(subq[0], subq[1], so4_params);
counter += 6;
}
// second layer of so4 blocks (odd)
for (int k = 1; k + block_size < n_qubits; k += 2) {
auto subq = q.slice(k, block_size);
auto so4_params = cudaq::slice_vector(params, counter, 6);
so4(subq[0], subq[1], so4_params);
counter += 6;
}
}
}
};
int main() {
// Read in the spin op from file
std::vector<double> h2_data{0, 0, 0, 0, -0.10647701149499994, 0.0,
1, 1, 1, 1, 0.0454063328691, 0.0,
1, 1, 3, 3, 0.0454063328691, 0.0,
3, 3, 1, 1, 0.0454063328691, 0.0,
3, 3, 3, 3, 0.0454063328691, 0.0,
2, 0, 0, 0, 0.170280101353, 0.0,
2, 2, 0, 0, 0.120200490713, 0.0,
2, 0, 2, 0, 0.168335986252, 0.0,
2, 0, 0, 2, 0.165606823582, 0.0,
0, 2, 0, 0, -0.22004130022499996, 0.0,
0, 2, 2, 0, 0.165606823582, 0.0,
0, 2, 0, 2, 0.174072892497, 0.0,
0, 0, 2, 0, 0.17028010135300004, 0.0,
0, 0, 2, 2, 0.120200490713, 0.0,
0, 0, 0, 2, -0.22004130022499999, 0.0,
15};
cudaq::spin_op H(h2_data, /*nQubits*/ 4);
// For 8 qubits, 36 parameters per layer
int n_layers = 2, n_qubits = H.num_qubits(), block_size = 2, p_counter = 0;
int n_blocks_per_layer = 2 * (n_qubits / block_size) - 1;
int n_params = n_layers * 6 * n_blocks_per_layer;
printf("%d qubit hamiltonian -> %d parameters\n", n_qubits, n_params);
// Define the initial parameters and ansatz.
auto init_params =
cudaq::random_vector(-1, 1, n_params, std::mt19937::default_seed);
so4_fabric ansatz;
auto argMapper = [&](std::vector<double> x) {
return std::make_tuple(x, n_qubits, n_layers);
};
// Run VQE.
cudaq::optimizers::lbfgs optimizer;
optimizer.initial_parameters = init_params;
optimizer.max_eval = 20;
optimizer.max_line_search_trials = 10;
cudaq::gradients::central_difference gradient(ansatz, argMapper);
auto [opt_val, opt_params] =
cudaq::vqe(ansatz, gradient, H, optimizer, n_params, argMapper);
printf("Optimal value = %.16lf\n", opt_val);
}
Let’s look at a more advanced variation of the previous example.
As an alternative to cudaq.vqe, we can also use the cudaq.optimizers suite on its own to write custom variational algorithm routines. Much of this can be slightly modified for use with third-party optimizers, such as scipy.
import cudaq
from cudaq import spin
from typing import List, Tuple
# We will be optimizing over a custom objective function that takes a vector
# of parameters as input and returns either the cost as a single float,
# or a tuple of (cost, gradient_vector) depending on the optimizer used.
# In this example, we will use the spin Hamiltonian and ansatz from `simple_vqe.py`
# and find the `angles` that minimize the expectation value of the system.
hamiltonian = 5.907 - 2.1433 * spin.x(0) * spin.x(1) - 2.1433 * spin.y(
0) * spin.y(1) + .21829 * spin.z(0) - 6.125 * spin.z(1)
@cudaq.kernel
def kernel(angles: List[float]):
qvector = cudaq.qvector(2)
x(qvector[0])
ry(angles[0], qvector[1])
x.ctrl(qvector[1], qvector[0])
# Define the optimizer that we'd like to use.
optimizer = cudaq.optimizers.Adam()
# Since we'll be using a gradient-based optimizer, we can leverage the
# CUDA-Q gradient helper class to automatically compute the gradient
# vector for us. The use of this class for gradient calculations is
# purely optional and can be replaced with your own custom gradient
# routine.
gradient = cudaq.gradients.CentralDifference()
def objective_function(parameter_vector: List[float],
hamiltonian=hamiltonian,
gradient_strategy=gradient,
kernel=kernel) -> Tuple[float, List[float]]:
"""
Note: the objective function may also take extra arguments, provided they
are passed into the function as default arguments in python.
"""
# Call `cudaq.observe` on the spin operator and ansatz at the
# optimizer provided parameters. This will allow us to easily
# extract the expectation value of the entire system in the
# z-basis.
# We define the call to `cudaq.observe` here as a lambda to
# allow it to be passed into the gradient strategy as a
# function. If you were using a gradient-free optimizer,
# you could purely define `cost = cudaq.observe().expectation()`.
get_result = lambda parameter_vector: cudaq.observe(
kernel, hamiltonian, parameter_vector).expectation()
# `cudaq.observe` returns a `cudaq.ObserveResult` that holds the
# counts dictionary and the `expectation`.
cost = get_result(parameter_vector)
print(f"<H> = {cost}")
# Compute the gradient vector using `cudaq.gradients.STRATEGY.compute()`.
gradient_vector = gradient_strategy.compute(parameter_vector, get_result,
cost)
# Return the (cost, gradient_vector) tuple.
return cost, gradient_vector
cudaq.set_random_seed(13) # make repeatable
energy, parameter = optimizer.optimize(dimensions=1,
function=objective_function)
print(f"\nminimized <H> = {round(energy,16)}")
print(f"optimal theta = {round(parameter[0],16)}")