QuESTlink—Mathematica embiggened by a hardware-optimised quantum emulator^*

While the forthcoming sections give a thorough overview of QuESTlink's facilities, we here provide a short summary to acquaint the reader with the essentials. QuESTlink is, offers, features or enables:

• Multithreaded and GPU-accelerated emulation of quantum computers.
• State-vector and density matrix simulation.
• Off-loading of simulation to remote hardware, through a backend-agnostic interface.
• Seamless integration with Mathematica's powerful and comprehensive toolset, and interactive notebook programing style.
• Rapid development with Mathematica's concise, functional language.
• A concise but expressive circuit language, akin to the symbolic description of circuits in the literature.
• Stand-alone with no required external downloading, compiling, or installation of any kind whatsoever.
• Through Mathematica, is compatible with all major operating systems.
• Free and open-source, hosted publicly on Github [8].
• Rendering of circuit diagrams, which can be exported to any file format through Mathematica.
• A suite of functions for higher-level calculations, like computing Hamiltonian expectation values, and derivatives of unitary circuits.
• Analytic generation of matrices from circuit specifications, which can include symbols.
• A comprehensive suite of unitary gates and decoherence channels, including multi-controlled multi-qubit general unitaries, and multi-qubit Kraus maps.
• A rigorous user-input validation scheme to catch user errors, like supplying non-unitary matrices or physically invalid noise parameters.

QuESTlink consists of both a Mathematica package (QuESTlink.m) and an underlying C++ programme (quest_link.cpp), which interface using C/Link; an implementation of the Wolfram Symbolic Transfer Protocol (WSTP) [9] (formerly called MathLink [10]). QuESTlink emulates quantum computers using the Quantum Exact Simulation Toolkit (QuEST) [5], which is an open-source, high-performance emulator written in C. QuEST is multi-threaded, GPU-accelerated and distributed, and the first two of these facilities are made available to QuESTlink.

C/Link facilitates conversion of Mathematica types, utilised in the user's notebook, into C types, and maps Mathematica-callable functions to C functions. The executable quest_link (quest_link.c compiled with WSTP) sits below as a 'backend' process, performing intermediate processing and invoking QuEST's API. This software stack is visualised in figure 1.

Zoom In Zoom Out Reset image size

Through this stack, QuESTlink offers a high-level Mathematica interface to quantum emulation which is numerically performed in C/C++ (a significantly faster but lower level language than Mathematica), using memory persistent in the C process, and potentially using accelerated hardware. The expensive numerical representations of the quantum states reside only in the C backend (and through QuEST, may also be persistent also in GPU memory), and each are identified by a unique ID. Quantum circuits are represented symbolically in their entirety in Mathematica, and sent to the C backend (at a small runtime overhead) only at emulation-time, compactly encoded into arrays of real numbers. The circuit language, and details of the QuESTlink interface, are presented in section 4.2.

QuESTlink leverages QuEST's extensive user-input validation, so that user errors are detected early and reported. For example, this includes validation that input matrices are unitary, that a gate's control and target qubits are unique, and that parameters of a decoherence channel result in a completely positive map. As of QuESTlink v0.3 (integrating QuEST v3.1.0), there are 72 unique forms of input validation. This helps users avoid logical errors in their code, to mitigate the risk of drawing incorrect conclusions from the results of QuESTlink.

Core QuEST performs validation before any modification to a quantum register. If a problem is detected, a C++ exception is thrown and caught by QuESTlink, and a descriptive error message is propagated to the Mathematica kernel. Some QuESTlink validation is performed on the fly. For example, a problem in a gate may only be detected after applying the prior gates in a circuit. In these instances, QuESTlink always reverts the affected register to its original state before the circuit was applied. This means no register is left in an uncertain state.

QuESTlink can be obtained and deployed entirely within a Mathematica notebook, without any installation or configuration. This is done by hosting the package code and quest_link executable(s) online in a Github repository, which the qtechtheory.org server redirects to. Figure 2 presents the process of serving the package to a user's Mathematica kernel.

Zoom In Zoom Out Reset image size

The quest_link executables for Windows, MacOS and Linux are automatically recompiled on the Github repository with every new release, using Github's integrated continuous-delivery service (Github actions [11]). This is for both security, ease of development, and its facility to push bug-fixes and updates seamlessly to the end user.

QuESTlink even enables quantum emulation using remote computing resources. Through WSTP, the quest_link process can run on a remote machine (e.g. a supercomputer) and communicate with the user's local Mathematica kernel via TCP/IP. The remote machine can employ more powerful hardware than available on the user's machine, and potentially emulate larger quantum states than can fit in the user's local memory. The protocol of off-loading a circuit emulation to remote hardware is outlined in figure 3. The tools and documentation for setting up a remote QuESTlink server are provided between the Github repo [8], and questlink.qtechtheory.org.

Zoom In Zoom Out Reset image size

In this remote configuration, Mathematica calls to QuESTlink functions will involve network communication with the remote environment, and hence incur overheads; this is worthwhile if the remote hardware sufficiently accelerates an expensive emulation which otherwise dominates runtime.

Despite QuESTlink's effort to minimise the communication cost, this network overhead will be prohibitively costly for some applications. For example, when large simulated quantum states undergo processing by the Mathematica kernel, and hence need to be copied back and forth between the localkernel and the remote environment. To mitigate this slowdown while still running QuESTlink remotely, one can launch their Mathematica kernel on the remote machine, using Mathematica's remote kernel facilities [12]. QuESTlink can also be used inside other Mathematica packages, and so be launched remotely and non-interactively, without a local notebook.

Before continuing, we offer a quick review of the Mathematica syntax used in this manuscript.

Evaluation of function f with input a is denoted by f[a], or equivalently f @ a. The expression g[f[a]] can be formed using prefix notation as g @ f @ a, or postfix notation as a //f //g. Matrix multiplication between (possibly complex) matrices a and b is denoted by a.b, and a^† denotes the conjugate transpose of a. Expressions with a trailing; suppress their otherwise displayed result. While f[] = a denotes immediate evaluation of a and assignment to f[], the syntax f[]:= a denotes delayed assignment, whereby later invoking f[] will evaluate a (which may have changed) each time. Elements of a list x = {a, b, c} are accessed as x[[i]] where index i ⩾ 1, and x[[m;;n]] returns the sublist spanning indices m to n. The shortcut ex /.a → b replaces sub-expressions of ex which match pattern a, with b. Comments appear in parentheses (* like this *).

Finally, for those copying code into Mathematica, subscripts can be quickly entered into a notebook with keyboard shortcut & .

In the code snippets featured in this manuscript, context should make clear what is input to the Mathematica notebook, and what is a rendered output.

The QuESTlink package requires no installation, and can be downloaded directly from within Mathematica:

The user then has a choice of several 'QuEST environments' in which to perform quantum emulation, all of which provide an identical user experience.

The first, CreateDownloadedQuESTEnv[], enables simulation on the user's machine by downloading a serial QuESTlink executable, pre-compiled for the user's operating system. This requires no apriori setup whatsoever.

CreateLocalQuESTEnv[fn] will attempt to launch an existing local QuESTlink executable, located at fn. This allows local simulation using serial, multi-core or GPU resources, depending on how the executable was compiled. Users can compile QuESTlink for their platform using the tools on the Github repo [8].

CreateRemoteQuESTEnv[ip, port1, port2] connects to a remote QuESTlink environment at the given ip address and ports. The remote machine can use serial, multithreading or GPU-acceleration to emulate quantum systems. The facilities to setup a remote QuESTlink server are provided in the Github repo [8].

Once connected to a QuEST environment, a full list of the supported QuESTlink facilities and gate symbols can be obtained by evaluating

and the documentation of a particular function or operator obtained similarly using ?.

Emulation begins by creating quantum registers (a 'Qureg'), each represented by a state-vector or density-matrix.

These functions return a unique ID for each Qureg, the memory for which is stored in the QuEST environment. Once created, QuESTlink provides a few functions to initialise Quregs:

or directly modify them (here, creating a random density matrix):

A quantum circuit u can be applied to a Qureg, agnostic of whether it is a state vector (effecting u |ψ⟩) or a density matrix (effecting u ρ u^†), using ApplyCircuit[u, qureg]. QuESTlink features a concise and expressive language for specifying gates and decoherence processes, where target qubits are denoted with subscript integers to gate symbols, and control gates 'wrap' the base gate. The below example makes use of the Circuit function, which disables commutation and allows a concise product representation of the circuit. In combination, this enables a syntax akin to how circuits are denoted in the quantum computing literature.

The result of Circuit[] is just a list of operators, allowing easy circuit manipulation. Mathematica features many tools for such lists allowing the user to easily extend, alter, join, compare, etc, their quantum circuit. For example:

These circuits can be swiftly visualised with DrawCircuit[u], which supports saving rendered circuit diagrams to raster and vector images.

After applying a unitary circuit, one can of course opt to obtain the entire state vector or density matrix, as for example:

However, it is typically more efficient to calculate desired quantities, like the expectation value of a Hamiltonian in the Pauli-basis, in the QuEST environment itself.

At any time, Quregs can be individually (or, all-together) destroyed to free up memory in the QuEST environment.

By leveraging Mathematica's powerful suite of symbolic calculations, QuESTlink can even build analytic matrices from circuit specifications.

With these facilities, QuESTlink offers a seamless integration with Mathematica's comprehensive range of computational and graphical tools, as illustrated in the following demonstrations.

To begin, we very compactly demonstrate the effect of two-qubit depolarising noise on the expected measurement of a simple Hamiltonian h. Starting in a random pure state ψ, a depolarising channel with probability 0.1 of any Pauli error occurring is repeatedly applied, in total 100 times.

Note that undisclosed variable opts contained additional code for customising the plot.

In this demonstration, we emulate the quantum variational imaginary-time simulation routine [14] to approximate the ground-state of a molecular Hamiltonian.

We first download a reduced 6-qubit representation of the electronic structure Hamiltonian of Lithium Hydride (LiH).

We create a simple 6-qubit ansatz circuit featuring 39 parameters, denoted with variables $\to {\theta }=\left\{{\theta }_{1},\dots ,\;{\theta }_{39}\right\}$.

This ansatz circuit consists of one and two qubit rotation gates; $\mathrm{exp}\left(-\mathrm{i}{\theta }_{j}\hat{\sigma }/2\right)$, for $\hat{\sigma }\in \left\{X,Y,Z,X\otimes X,Y\otimes Y,Z\otimes Z\right\}$, which in QuESTlink are denoted by Rx_q[θ] (etc) and R[θ, X_q1X_q2]. In the diagram below, each such paired two-qubit rotation is marked by a vertical link.

We now create several Quregs with which to emulate the quantum algorithm.

ψ will be maintained as the output state of the ansatz circuit, and hψ will store the result of applying the Hamiltonian h to ψ. ϕ will merely provide intermediate work-space for calculations, and dψ will store a Qureg for each parameter in the ansatz (a total of nθ).

Next, we choose a random initial assignment of the ansatz parameters, and measure the energy of the resulting quantum state.

The variational imaginary-time algorithm [14] involves repeatedly measuring a matrix and vector of quantities,

$$\begin{equation}{A}_{ij}=\text{Re}\left\{\left\langle \frac{\partial \psi \left(\to {\theta }\right)}{\partial {\theta }_{i}}\vert \frac{\partial \psi \left(\to {\theta }\right)}{\partial {\theta }_{j}}\right\rangle \right\},\end{equation} \tag{ 1 }$$

$$\begin{equation}{C}_{i}=-\text{Re}\left\{\left\langle \psi \left(\to {\theta }\right)\vert h\vert \frac{\partial \psi \left(\to {\theta }\right)}{\partial {\theta }_{i}}\right\rangle \right\}\end{equation} \tag{ 2 }$$

and iteratively updating the parameters under

$$\begin{equation}\to {\theta }\to \to {\theta }+{\Delta}t\;{A}^{-1}\;\to {C},\end{equation} \tag{ 3 }$$

which we now concisely emulate for nt iterations.

The energy of the output quantum state, as produced by the reached assignment of the parameters, is close to the true groundstate found through matrix diagonalisation.

Interestingly, a significantly slower but direct minimisation in Mathematica reveals the energy reached by imaginary-time was not quite the best possible of our chosen ansatz.

We invite the interested reader to compare the Mathematica code above to a native C implementation of imaginary time evolution, as hosted on Github [15].

In this demonstration, we emulate Trotterisation of a spin-ring Hamiltonian, with and without noise, and compare it to direct numerical solving of the Schrödinger equation. Our Hamiltonian is the one-dimensional nearest-neighbour (periodic boundary conditions) nQb-spin Heisenberg model with a random magnetic field in the z direction,

$$\begin{equation}\sum _{j}^{\mathtt{\text{nQb}}}{\to {\sigma }}_{j}{\to {\sigma }}_{j+1}+{r}_{j}{\sigma }_{j}^{z},\;\;\;\;{r}_{j}\sim \mathcal{U}\left[-1,1\right].\end{equation} \tag{ 4 }$$

Real-time simulation of this model to time t = nQb is nominated by Childs et al [16] as an early practical application of a quantum computer. In this example, we will study a five-spin chain Hamiltonian h.

To simulate evolution on a quantum computer, we will emulate circuits formed by the Suzuki-Trotter decompositions [17] of the unitary evolution operator of varying order n. Below, r is the number of repetitions of the order-n circuit to perform, to ultimately reach time t.

Note the original ordering of the Hamiltonian terms is arbitrary, and should ideally be optimised into commuting groups, or at least randomly shuffled. In particular, Childs et al's uniform randomisation scheme [18], whereby each repetition sees a random Hamiltonian ordering, is trivially implemented as:

Even Campbell's qDRIFT routine [19] can be given a compact implementation.

However, we opt instead to utilise only deterministic Trotterization for clarity. For example, the first-order (n = 1) single-repetition (r = 1) Trotter circuit (arbitrarily to time t = 1) has the form:

We compare the states produced from Trotter circuits with the 'true' (to numerical precision) time evolution, found by numerically solving the Schrödinger equation using Mathematica's in-built NDSolveValue routine.

In the proceeding solutions, the Qureg ψ0 will store the initial state (a random pure state), ψt will store the true state, and ψ will store outputs of the Trotter circuits.

We can now commence the computation by generating circuits of a varying Trotter order and number of repetitions, emulating their execution, and computing the fidelity of their output state with the true state.

We then plot the result below, using some undisclosed parameters opts to tweak the plot style. The horizontal axis is the total number of gates in the Trotter circuit.

We can even imitate an noisy quantum device by inserting decoherence operators into our circuit. In this demonstration, we will follow each Trotter-prescribed gate with one or two qubit depolarising operators, which effect the channel

$$\begin{equation}\mathtt{\text{Depol}}\left[\rho \right]=\left(1-p\right)\rho +\frac{p}{3}\sum _{\sigma }\sigma \rho \sigma ,\end{equation} \tag{ 5 }$$

(noting σ = σ^† for Paulis) and similarly for the two-qubit analogue. We choose an error rate of p = 10⁻⁴, below the rates recently demonstrated in Google's Sycamore machine [2]. Below is a visualisation of the first five unitary gates of the resulting noisy circuit.

Our simulation of the noisy circuit is near-identical to our simulation of pure states, with the exception that we now operate upon density matrices.

As one might expect, the accuracy of the Trotter circuits have waned, and eventually decrease with increased circuit depth, due to the opportunity for additional errors to accrue.

Complete notebooks demonstrating some extended computations in QuESTlink are available at questlink.qtechtheory.org.