# 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/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;
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.

# Since we'll be using a gradient-based optimizer, we can leverage the
# vector for us. The use of this class for gradient calculations is
# purely optional and can be replaced with your own custom gradient
# routine.

def objective_function(parameter_vector: List[float],
hamiltonian=hamiltonian,
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}")