CUDA-Q Dynamics

CUDA-Q enables the design, simulation and execution of quantum dynamics via the evolve API. Specifically, this API allows us to solve the time evolution of quantum systems or models. In the simulation mode, CUDA-Q provides the dynamics backend target, which is based on the cuQuantum library, optimized for performance and scale on NVIDIA GPU.

Quick Start

In the example below, we demonstrate a simple time evolution simulation workflow comprising of the following steps:

  1. Define a quantum system model

A quantum system model is defined by a Hamiltonian. For example, a superconducting transmon qubit can be modeled by the following Hamiltonian

\[H = \frac{\omega_z}{2} \sigma_z + \omega_x \cos(\omega_d t)\sigma_x,\]

where \(\sigma_z\) and \(\sigma_x\) are Pauli Z and X operators, respectively.

Using CUDA-Q operator, the above time-dependent Hamiltonian can be set up as follows.

import numpy as np
from cudaq.operator import *

# Qubit Hamiltonian
hamiltonian = 0.5 * omega_z * spin.z(0)
# Add modulated driving term to the Hamiltonian
hamiltonian += omega_x * ScalarOperator(lambda t: np.cos(omega_d * t)) * spin.x(
    0)

In particular, ScalarOperator provides an easy way to model arbitrary time-dependent control signals. Details about CUDA-Q operator, including builtin operators that it supports can be found here.

  1. Setup the evolution simulation

The below code snippet shows how to simulate the time-evolution of the above system with cudaq.evolve.

import cudaq
import cupy as cp

# Set the target to our dynamics simulator
cudaq.set_target("dynamics")

# Dimensions of sub-systems: a single two-level system.
dimensions = {0: 2}

# Initial state of the system (ground state).
rho0 = cudaq.State.from_data(
    cp.array([[1.0, 0.0], [0.0, 0.0]], dtype=cp.complex128))

# Schedule of time steps.
steps = np.linspace(0, t_final, n_steps)
schedule = Schedule(steps, ["t"])

# Run the simulation.
evolution_result = evolve(hamiltonian,
                          dimensions,
                          schedule,
                          rho0,
                          observables=[spin.x(0),
                                       spin.y(0),
                                       spin.z(0)],
                          collapse_operators=[],
                          store_intermediate_results=True)

Specifically, we need to set up the simulation by providing:

  • The system model in terms of a Hamiltonian as well as any decoherence terms, so-called collapse_operators.

  • The dimensionality of component systems in the model. CUDA-Q evolve allows users to model arbitrary multi-level systems, such as photonic Fock space.

  • The initial quantum state.

  • The time schedule, aka time steps, of the evolution.

  • Any ‘observable’ operator that we want to measure the expectation value with respect to the evolving state.

Note

By default, evolve will only return the final state and expectation values. To save intermediate results (at each time step specified in the schedule), the store_intermediate_results flag must be set to True.

  1. Retrieve and plot the results

Once the simulation is complete, we can retrieve the final state and the expectation values as well as intermediate values at each time step (with store_intermediate_results=True).

For example, we can plot the Pauli expectation value for the above simulation as follows.

get_result = lambda idx, res: [
    exp_vals[idx].expectation() for exp_vals in res.expectation_values()
]

import matplotlib.pyplot as plt

plt.plot(steps, get_result(0, evolution_result))
plt.plot(steps, get_result(1, evolution_result))
plt.plot(steps, get_result(2, evolution_result))
plt.ylabel("Expectation value")
plt.xlabel("Time")
plt.legend(("Sigma-X", "Sigma-Y", "Sigma-Z"))

In particular, for each time step, evolve captures an array of expectation values, one for each observable. Hence, we convert them into sequences for plotting purposes.

Operator

CUDA-Q provides builtin definitions for commonly-used operators, such as the ladder operators (\(a\) and \(a^\dagger\)) of a harmonic oscillator, the Pauli spin operators for a two-level system, etc.

Here is a list of those operators.

Builtin Operators

Name

Description

identity

Identity operator

zero

Zero or null operator

annihilate

Bosonic annihilation operator (\(a\))

create

Bosonic creation operator (\(a^\dagger\))

number

Number operator of a bosonic mode (equivalent to \(a^\dagger a\))

parity

Parity operator of a bosonic mode (defined as \(e^{i\pi a^\dagger a}\))

displace

Displacement operator of complex amplitude \(\alpha\) (displacement). It is defined as \(e^{\alpha a^\dagger - \alpha^* a}\).

squeeze

Squeezing operator of complex squeezing amplitude \(z\) (squeezing). It is defined as \(\exp(\frac{1}{2}(z^*a^2 - z a^{\dagger 2}))\).

position

Position operator (equivalent to \((a^\dagger + a)/2\))

momentum

Momentum operator (equivalent to \(i(a^\dagger - a)/2\))

spin.x

Pauli \(\sigma_x\) operator

spin.y

Pauli \(\sigma_y\) operator

spin.z

Pauli \(\sigma_z\) operator

spin.plus

Pauli raising (\(\sigma_+\)) operator

spin.minus

Pauli lowering (\(\sigma_-\)) operator

As an example, let’s look at the Jaynes-Cummings model, which describes the interaction between a two-level atom and a light (Boson) field.

Mathematically, the Hamiltonian can be expressed as

\[H = \omega_c a^\dagger a + \omega_a \frac{\sigma_z}{2} + \frac{\Omega}{2}(a\sigma_+ + a^\dagger \sigma_-).\]

This Hamiltonian can be converted to CUDA-Q Operator representation with

hamiltonian = omega_c * operators.create(1) * operators.annihilate(1) \
                + (omega_a / 2) * spin.z(0) \
                + (Omega / 2) * (operators.annihilate(1) * spin.plus(0) + operators.create(1) * spin.minus(0))

In the above code snippet, we map the cavity light field to degree index 1 and the two-level atom to degree index 0. The description of composite quantum system dynamics is independent from the Hilbert space of the system components. The latter is specified by the dimension map that is provided to the cudaq.evolve call.

Time-Dependent Dynamics

In the previous examples of operator construction, we assumed that the systems under consideration were described by time-independent Hamiltonian. However, we may want to simulate systems whose Hamiltonian operators have explicit time dependence.

CUDA-Q provides multiple ways to construct time-dependent operators.

  1. Time-dependent coefficient

CUDA-Q ScalarOperator can be used to wrap a Python function that returns the coefficient value at a specific time.

As an example, we will look at a time-dependent Hamiltonian of the form \(H = H_0 + f(t)H_1\), where \(f(t)\) is the time-dependent driving strength given as \(cos(\omega t)\).

The following code sets up the problem

# Define the static (drift) and control terms
H0 = spin.z(0)
H1 = spin.x(0)
H = H0 + ScalarOperator(lambda t: np.cos(omega * t)) * H1

  1. Time-dependent operator

We can also construct a time-dependent operator from a function that returns a complex matrix representing the time dynamics of that operator.

As an example, let’s looks at the displacement operator. It can be defined as follows:

import numpy
import scipy
from cudaq.operator import *
from numpy.typing import NDArray


def displacement_matrix(
        dimension: int,
        displacement: NumericType) -> NDArray[numpy.complexfloating]:
    """
    Returns the displacement operator matrix.
    Args:
        displacement: Amplitude of the displacement operator.
            See also https://en.wikipedia.org/wiki/Displacement_operator.
    """
    displacement = complex(displacement)
    term1 = displacement * operators.create(0).to_matrix({0: dimension})
    term2 = numpy.conjugate(displacement) * operators.annihilate(0).to_matrix(
        {0: dimension})
    return scipy.linalg.expm(term1 - term2)


# The second argument here indicates the the defined operator
# acts on a single degree of freedom, which can have any dimension.
# An argument [2], for example, would indicate that it can only
# act on a single degree of freedom with dimension two.
ElementaryOperator.define("displace", [0], displacement_matrix)


def displacement(degree: int) -> ElementaryOperator:
    """
    Instantiates a displacement operator acting on the given degree of freedom.
    """
    return ElementaryOperator("displace", [degree])

The defined operator is parameterized by the displacement amplitude. To create simulate the evolution of an operator under a time dependent displacement amplitude, we can define how the amplitude changes in time:

import cudaq

# Define a system consisting of a single degree of freedom (0) with dimension 3.
system_dimensions = {0: 3}
system_operator = displacement(0)

# Define the time dependency of the system operator as a schedule that linearly
# increases the displacement parameter from 0 to 1.
time_dependence = Schedule(numpy.linspace(0, 1, 100), ['displacement'])
initial_state = cudaq.State.from_data(
    numpy.ones(3, dtype=numpy.complex128) / numpy.sqrt(3))

# Simulate the evolution of the system under this time dependent operator.
cudaq.evolve(system_operator, system_dimensions, time_dependence, initial_state)

Let’s say we want to add a squeezing term to the system operator. We can independently vary the squeezing amplitude and the displacement amplitude by instantiating a schedule with a custom function that returns the desired value for each parameter:

system_operator = displacement(0) + operators.squeeze(0)


# Define a schedule such that displacement amplitude increases linearly in time
# but the squeezing amplitude decreases, that is follows the inverse schedule.
def parameter_values(time_steps):

    def compute_value(param_name, step_idx):
        match param_name:
            case 'displacement':
                return time_steps[int(step_idx)]
            case 'squeezing':
                return time_steps[-int(step_idx + 1)]
            case _:
                raise ValueError(f"value for parameter {param_name} undefined")

    return Schedule(range(len(time_steps)), system_operator.parameters.keys(),
                    compute_value)


time_dependence = parameter_values(numpy.linspace(0, 1, 100))
cudaq.evolve(system_operator, system_dimensions, time_dependence, initial_state)

Numerical Integrators

CUDA-Q provides a set of numerical integrators, to be used with the dynamics backend target.

Numerical Integrators

Name

Description

RungeKuttaIntegrator

Explicit 4th-order Runge-Kutta method (default integrator)

ScipyZvodeIntegrator

Complex-valued variable-coefficient ordinary differential equation solver (provided by SciPy)

CUDATorchDiffEqDopri5Integrator

Runge-Kutta of order 5 of Dormand-Prince-Shampine (provided by torchdiffeq)

CUDATorchDiffEqAdaptiveHeunIntegrator

Runge-Kutta of order 2 (provided by torchdiffeq)

CUDATorchDiffEqBosh3Integrator

Runge-Kutta of order 3 of Bogacki-Shampine (provided by torchdiffeq)

CUDATorchDiffEqDopri8Integrator

Runge-Kutta of order 8 of Dormand-Prince-Shampine (provided by torchdiffeq)

CUDATorchDiffEqEulerIntegrator

Euler method (provided by torchdiffeq)

CUDATorchDiffEqExplicitAdamsIntegrator

Explicit Adams-Bashforth method (provided by torchdiffeq)

CUDATorchDiffEqImplicitAdamsIntegrator

Implicit Adams-Bashforth-Moulton method (provided by torchdiffeq)

CUDATorchDiffEqMidpointIntegrator

Midpoint method (provided by torchdiffeq)

CUDATorchDiffEqRK4Integrator

Fourth-order Runge-Kutta with 3/8 rule (provided by torchdiffeq)

Note

To use Torch-based integrators, users need to install torchdiffeq (e.g., with pip install torchdiffeq). This is an optional dependency of CUDA-Q, thus will not be installed by default.

Warning

Torch-based integrators require a CUDA-enabled Torch installation. Depending on your platform (e.g., aarch64), the default Torch pip package may not have CUDA support.

The below command can be used to verify your installation:

python3 -c "import torch; print(torch.version.cuda)"

If the output is a ‘None’ string, it indicates that your Torch installation does not support CUDA. In this case, you need to install a CUDA-enabled Torch package via other mechanisms, e.g., building Torch from source or using their Docker images.