CUDA Quantum in Python

Welcome to CUDA Quantum! This is a introduction by example for using CUDA Quantum in Python.

Introduction

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 state vector 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 this publication.

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 is 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 = ...
# 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, get_result,
                                                cost)

    # 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 Max-Cut 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 Max-Cut of a rectangular graph illustrated below.

#       v0  0---------------------0 v1
#           |                     |
#           |                     |
#           |                     |
#           |                     |
#       v3  0---------------------0 v2
# The Max-Cut 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 Max-Cut"""
    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()

Using Quantum Hardware Providers

CUDA Quantum contains support for using a set of hardware providers. For more information about executing quantum kernels on different hardware backends, please take a look at CUDA Quantum Hardware Backends.

The following code illustrates how run kernels on Quantinuum’s backends.

import cudaq

# You only have to set the target once! No need to redefine it
# for every execution call on your kernel.
# By default, we will submit to the Quantinuum syntax checker.
cudaq.set_target("quantinuum")

# Create the kernel we'd like to execute on Quantinuum.
kernel = cudaq.make_kernel()
qubits = kernel.qalloc(2)
kernel.h(qubits[0])
kernel.cx(qubits[0], qubits[1])
kernel.mz(qubits[0])
kernel.mz(qubits[1])

# Submit to Quantinuum's endpoint and confirm the program is valid.

# Option A:
# By using the synchronous `cudaq.sample`, the execution of
# any remaining classical code in the file will occur only
# after the job has been executed by the Quantinuum service.
# We will use the synchronous call to submit to the syntax
# checker to confirm the validity of the program.
syntax_check = cudaq.sample(kernel)
if (syntax_check):
    print("Syntax check passed! Kernel is ready for submission.")

# Now we can update the target to the Quantinuum emulator and
# execute our program.
cudaq.set_target("quantinuum", machine="H1-2E")

# Option B:
# By using the asynchronous `cudaq.sample_async`, the remaining
# classical code will be executed while the job is being handled
# by Quantinuum. This is ideal when submitting via a queue over
# the cloud.
async_results = cudaq.sample_async(kernel)
# ... more classical code to run ...

# We can either retrieve the results later in the program with
# ```
# async_counts = async_results.get()
# ```
# or wee can also write the job reference (`async_results`) to
# a file and load it later or from a different process.
file = open("future.txt", "w")
file.write(str(async_results))
file.close()

# We can later read the file content and retrieve the job
# information and results.
same_file = open("future.txt", "r")
retrieved_async_results = cudaq.AsyncSampleResult(str(same_file.read()))

counts = retrieved_async_results.get()
print(counts)

The following code illustrates how run kernels on IonQ’s backends.

import cudaq

# You only have to set the target once! No need to redefine it
# for every execution call on your kernel.
# To use different targets in the same file, you must update
# it via another call to `cudaq.set_target()`
cudaq.set_target("ionq")

# Create the kernel we'd like to execute on IonQ.
kernel = cudaq.make_kernel()
qubits = kernel.qalloc(2)
kernel.h(qubits[0])
kernel.cx(qubits[0], qubits[1])

# Note: All qubits will be measured at the end upon performing
# the sampling. You may encounter a pre-flight error on IonQ
# backends if you include explicit measurements.

# Execute on IonQ and print out the results.

# Option A:
# By using the asynchronous `cudaq.sample_async`, the remaining
# classical code will be executed while the job is being handled
# by IonQ. This is ideal when submitting via a queue over
# the cloud.
async_results = cudaq.sample_async(kernel)
# ... more classical code to run ...

# We can either retrieve the results later in the program with
# ```
# async_counts = async_results.get()
# ```
# or wee can also write the job reference (`async_results`) to
# a file and load it later or from a different process.
file = open("future.txt", "w")
file.write(str(async_results))
file.close()

# We can later read the file content and retrieve the job
# information and results.
same_file = open("future.txt", "r")
retrieved_async_results = cudaq.AsyncSampleResult(str(same_file.read()))

counts = retrieved_async_results.get()
print(counts)

# Option B:
# By using the synchronous `cudaq.sample`, the execution of
# any remaining classical code in the file will occur only
# after the job has been returned from IonQ.
counts = cudaq.sample(kernel)
print(counts)