CUDA Quantum in Python

Examples for running CUDA Quantum programs written in Python on the statevector simulator included in the cudaq package.

Introduction

Welcome to CUDA Quantum! We’re going to take a look at how to construct quantum programs through CUDA Quantum’s Kernel API.

When you create a Kernel and invoke its methods, a quantum program is constructed that can then be executed by calling, for example, cudaq::sample. Let’s take a closer look!

import cudaq

# We begin by defining the `Kernel` that we will construct our
# program with.
kernel = cudaq.make_kernel()

# Next, we can allocate qubits to the kernel via `qalloc(qubit_count)`.
# An empty call to `qalloc` will return a single qubit.
qubit = kernel.qalloc()

# Now we can begin adding instructions to apply to this qubit!
# Here we'll just add every non-parameterized
# single qubit gate that is supported by CUDA Quantum.
kernel.h(qubit)
kernel.x(qubit)
kernel.y(qubit)
kernel.z(qubit)
kernel.t(qubit)
kernel.s(qubit)

# Next, we add a measurement to the kernel so that we can sample
# the measurement results on our simulator!
kernel.mz(qubit)

# Finally, we can execute this kernel on the statevector simulator
# by calling `cudaq.sample`. This will execute the provided kernel
# `shots_count` number of times and return the sampled distribution
# as a `cudaq.SampleResult` dictionary.
result = cudaq.sample(kernel)

# Now let's take a look at the `SampleResult` we've gotten back!
print(result)  # or result.dump()

Bernstein-Vazirani

Bernstein Vazirani is an algorithm for finding the bitstring encoded in a given function.

For the original source of this algorithm, see: https://epubs.siam.org/doi/10.1137/S0097539796300921

In this example, we generate a random bitstring, encode it into an inner-product oracle, then we simulate the kernel and return the most probable bitstring from its execution.

If all goes well, the state measured with the highest probability should be our hidden bitstring!

import cudaq
import random


def random_bitstring(length: int):
    bitstring = ""
    for bit in range(length):
        bitstring += str(random.randint(0, 1))
    return bitstring


def oracle(kernel: cudaq.Kernel, register: cudaq.QuakeValue,
           auxillary_qubit: cudaq.QuakeValue, hidden_bitstring: str):
    """
    The inner-product oracle for Bernstein Vazirani.
    """
    for index, bit in enumerate(hidden_bitstring):
        if bit == "0":
            # Apply identity operation to the qubit if it's
            # to be in the 0-state.
            # In this case, we do nothing.
            pass
        else:
            # Otherwise, apply a `cx` gate with the current qubit as
            # the control and the auxillary qubit as the target.
            kernel.cx(control=register[index], target=auxillary_qubit)


def bernstein_vazirani(qubit_count: int):
    """
    Returns a kernel implementing the Bernstein Vazirani algorithm
    for a random, hidden bitstring.
    """
    kernel = cudaq.make_kernel()
    # Allocate the specified number of qubits - this
    # corresponds to the length of the hidden bitstring.
    qubits = kernel.qalloc(qubit_count)
    # Allocate an extra auxillary qubit.
    auxillary_qubit = kernel.qalloc()

    # Prepare the auxillary qubit.
    kernel.h(auxillary_qubit)
    kernel.z(auxillary_qubit)

    # Place the rest of the register in a superposition state.
    kernel.h(qubits)

    # Generate a random, hidden bitstring for the oracle
    # to encode. Note: we define the bitstring here so
    # as to be able to return it for verification.
    hidden_bitstring = random_bitstring(qubit_count)

    # Query the oracle.
    oracle(kernel, qubits, auxillary_qubit, hidden_bitstring)

    # Apply another set of Hadamards to the register.
    kernel.h(qubits)

    # Apply measurement gates to just the `qubits` (excludes the
    # auxillary qubit).
    kernel.mz(qubits)
    return kernel, hidden_bitstring


qubit_count = 5
kernel, hidden_bitstring = bernstein_vazirani(qubit_count)

result = cudaq.sample(kernel)
print(f"encoded bitstring = {hidden_bitstring}")
print(f"measured state = {result.most_probable()}")
print(f"Were we successful? {hidden_bitstring == result.most_probable()}")

Variational Quantum Eigensolver

Let’s take a look at how we can use CUDA Quantum’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 Quantum’s optimizers, cudaq.vqe will find and return the optimal set of parameters that minimize the energy, <Z>, of the system.

import cudaq
from cudaq import spin

# 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.
kernel, thetas = cudaq.make_kernel(list)
# Allocate 2 qubits.
qubits = kernel.qalloc(2)
kernel.x(qubits[0])
# Apply an `ry` gate that's parameterized by the first
# `QuakeValue` entry of our list, `thetas`.
kernel.ry(thetas[0], qubits[1])
kernel.cx(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 = XXXX
# optimizer...

# Finally, we can pass all of that into `cudaq.vqe` and it will automatically run our
# optimization loop and return 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,3)}")
print(f"optimal theta = {round(parameter[0],3)}")

Let’s look at a more advanced examples.

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 in a tuple of (cost, gradient_vector) depending on the optimizer used.

# In this case, we will use the spin Hamiltonian and ansatz from `simple_vqe.py`
# and find the `thetas` 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)

kernel, thetas = cudaq.make_kernel(list)
qubits = kernel.qalloc(2)
kernel.x(qubits[0])
kernel.ry(thetas[0], qubits[1])
kernel.cx(qubits[1], qubits[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
# CUDA Quantum's 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_z()`.
    get_result = lambda parameter_vector: cudaq.observe(
        kernel, hamiltonian, parameter_vector, shots_count=100).expectation_z()
    # `cudaq.observe` returns a cudaq.ObserveResult` that holds the
    # counts dictionary and the `expectation_z`.
    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=parameter_vector, function=get_result)

    # Return the (cost, gradient_vector) tuple.
    return cost, gradient_vector


energy, parameter = optimizer.optimize(dimensions=1,
                                       function=objective_function)

print(f"\nminimized <H> = {round(energy,3)}")
print(f"optimal theta = {round(parameter[0],3)}")

Quantum Approximate Optimization Algorithm

Let’s now see how we can leverage the VQE algorithm to compute the maxcut of a rectangular graph.

import cudaq
from cudaq import spin

import numpy as np

# Here we build up a kernel for QAOA with p layers, with each layer
# containing the alternating set of unitaries corresponding to the problem
# and the mixer Hamiltonians. The algorithm leverages the VQE algorithm
# to compute the maxcut of a rectangular graph illustrated below.

#       v0  0---------------------0 v1
#           |                     |
#           |                     |
#           |                     |
#           |                     |
#       v3  0---------------------0 v2
# The maxcut for this problem is 0101 or 1010.

# The problem Hamiltonian
hamiltonian = 0.5 * spin.z(0) * spin.z(1) + 0.5 * spin.z(1) * spin.z(2) \
       + 0.5 * spin.z(0) * spin.z(3) + 0.5 * spin.z(2) * spin.z(3)

# Problem parameters.
qubit_count: int = 4
layer_count: int = 2
parameter_count: int = 2 * layer_count


def kernel_qaoa() -> cudaq.Kernel:
    """QAOA ansatz for maxcut"""
    kernel, thetas = cudaq.make_kernel(list)
    qreg = kernel.qalloc(qubit_count)

    # Create superposition
    kernel.h(qreg)

    # Loop over the layers
    for i in range(layer_count):
        # Loop over the qubits
        # Problem unitary
        for j in range(qubit_count):
            kernel.cx(qreg[j], qreg[(j + 1) % qubit_count])
            kernel.rz(2.0 * thetas[i], qreg[(j + 1) % qubit_count])
            kernel.cx(qreg[j], qreg[(j + 1) % qubit_count])

        # Mixer unitary
        for j in range(qubit_count):
            kernel.rx(2.0 * thetas[i + layer_count], qreg[j])

    return kernel


# Specify the optimizer and its initial parameters.
optimizer = cudaq.optimizers.COBYLA()
optimizer.initial_parameters = np.random.uniform(-np.pi / 8.0, np.pi / 8.0,
                                                 parameter_count)
print("Initial parameters = ", optimizer.initial_parameters)

# Pass the kernel, spin operator, and optimizer to `cudaq.vqe`.
optimal_expectation, optimal_parameters = cudaq.vqe(
    kernel=kernel_qaoa(),
    spin_operator=hamiltonian,
    optimizer=optimizer,
    parameter_count=parameter_count)

# Print the optimized value and its parameters
print("Optimal value = ", optimal_expectation)
print("Optimal parameters = ", optimal_parameters)

# Sample the circuit using the optimized parameters
counts = cudaq.sample(kernel_qaoa(), optimal_parameters)
counts.dump()