Exploring cuQuantum (1) — From GTC2022

Purpose

I briefly learned about cuQuantum at GTC2022, so I would like to summarize it in this article.

Specifically, I watched Scaling Quantum Circuit Simulations with cuQuantum for Quantum Algorithms [A41102].

I am writing this as a sort of personal note, so the intermediate calculations might be a bit messy and I haven't delved deeply into the underlying meanings, but I will consider those in the future. Also, as a bonus, I will take a brief look at the calculation of circuits with a large number of qubits.

What is cuQuantum?

From reference [CQ]:

NVIDIA cuQuantum is an SDK of optimized libraries and tools for accelerating quantum computing workflows. With NVIDIA Tensor Core GPUs, developers can use cuQuantum to speed up quantum circuit simulations based on state vector and tensor network methods by orders of magnitude.

Basically, it is an incredibly fast GPU simulator for quantum computers. For tensor networks, see for example reference [TN].

Let's install cuQuantum

Since my local environment is CUDA 11.0, I will use the following command:

pip install cuquantum cuquantum-python cupy-cuda110

That should do it. Incidentally, the content I will discuss below uses features from cuquantum-python 22.7.0 or later^[1], which require Python 3.8+. On the other hand, Google Colaboratory's Jupyter kernel seems to be 3.7.12, and only supports up to cuquantum-python 22.5.0. This means the API cannot be called, so I will give up on using Colab this time.

The experimental environment uses an NVIDIA T4, and I am not doing anything too elaborate. Since CuPy is used as the backend, the specific type of the tensors appearing below is cupy.ndarray.

What am I going to do?

I will convert a Qiskit quantum circuit into a tensor network, perform contraction calculations, and find the probability amplitudes of the circuit's output states.

Pros

• Since it doesn't maintain the state vector, there's no need to perform matrix calculations while holding a vector of length $2^n$ for $n$ qubits. In other words, it is memory-efficient.

Cons

• This is just a calculation in the classical supercomputer world and doesn't benefit from quantum mechanics, so, for example, superposition cannot be utilized, and all probability amplitudes cannot be obtained at once. Probably.

Let's look at a simple sample right away

First, let's import the necessary modules as usual.

from qiskit import QuantumCircuit, QuantumRegister
from cuquantum import CircuitToEinsum, contract
import numpy as np

Then, if we construct a quantum circuit and pass it to cuQuantum:

qc = QuantumCircuit(1)
qc.x(0)
converter = CircuitToEinsum(qc)
expr, operands = converter.amplitude(bitstring='1')
print(expr)
for op in operands:
    print(op)
    print()
print(contract(expr, *operands))

We get the following output.

a,ba,b->
[1.+0.j 0.+0.j]

[[0.+0.j 1.+0.j]
 [1.+0.j 0.+0.j]]

[0.+0.j 1.+0.j]

(1+0j)

The result of the contraction calculation contract(expr, *operands) is (1+0j), i.e., 1. What does this mean?

Let's do manual calculations

Writing operands[0], operands[1], and operands[2] as mathematical expressions gives the following:

$$\begin{align*} \psi = \begin{pmatrix} 1 \\ 0 \end{pmatrix},\ \xi = \begin{pmatrix} 0 \\ 1 \end{pmatrix},\ X = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \end{align*}$$

Since expr = 'a,ba,b->', writing it as a mathematical expression gives:

$$\begin{align*} \psi_{a} X_{ba} \xi_{b} &= \sum_{a,b} \psi_{a} X_{ba} \xi_{b} \\ &= (\xi_0\ \ \xi_1) \begin{pmatrix} x_{00} & x_{01} \\ x_{10} & x_{11} \end{pmatrix} \begin{pmatrix} \psi_0 \\ \psi_1 \end{pmatrix} \\ & = (0\ \ 1) \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} 1 \\ 0 \end{pmatrix} = 1 \end{align*}$$

In other words, it is calculating $\braket{\xi|X|\psi}$. This is exactly the calculation often seen in textbooks like, "given a quantum circuit with only an $X$ gate, calculating the probability amplitude of state $\ket{1}$ results in $p_{X\ket{0}}(\ket{1}) = 1$."

It turns out we aren't doing anything scary. Next, let's consider a slightly more complex circuit.

Let's consider an entangled state

A so-called Bell state is created with a circuit like this:

qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)

Let's put this into cuQuantum:

converter = CircuitToEinsum(qc, dtype='complex128', backend='cupy')
for bitstring in ['00', '01', '10', '11']:
    expr, operands = converter.amplitude(bitstring=bitstring)
    amplitude = contract(expr, *operands)
    print(amplitude)

(0.7071067811865475+0j)
0j
0j
(0.7071067811865475+0j)

What I want to say is that by calculating the probability amplitudes of the states $\ket{00}$, $\ket{01}$, $\ket{10}$, and $\ket{11}$, we can see that the circuit's output is $\frac{1}{\sqrt{2}}(\ket{00} + \ket{11})$.

Let's do manual calculations again

This time, let's find only the probability amplitude of state $\ket{11}$.

expr, operands = converter.amplitude(bitstring='11')
print(expr)
for op in operands:
    print(op)
    print()

a,b,ca,debc,e,d->
[1.+0.j 0.+0.j]

[1.+0.j 0.+0.j]

[[ 0.70710678+0.j  0.70710678+0.j]
 [ 0.70710678+0.j -0.70710678+0.j]]

[[[[1.+0.j 0.+0.j]
   [0.+0.j 0.+0.j]]

  [[0.+0.j 0.+0.j]
   [0.+0.j 1.+0.j]]]]

The results are quite daunting. Let's write them as mathematical expressions again:

$$\begin{align*} \psi &= \begin{pmatrix} 1 \\ 0 \end{pmatrix},\ \varphi = \begin{pmatrix} 1 \\ 0 \end{pmatrix},\ H = (h_{ij}) = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} \\ S &= (s_{ij}^{k\ell}) = \Bigg[ \left[ \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}, \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix} \right], \left[ \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix}, \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} \right] \Bigg], \\ \xi &= \begin{pmatrix} 0 \\ 1 \end{pmatrix},\ \eta = \begin{pmatrix} 0 \\ 1 \end{pmatrix} \end{align*}$$

These might not be the best symbols, but I've chosen them this way. Since expr = 'a,b,ca,debc,e,d->', writing this calculation as a mathematical expression gives:

$$\begin{align*} \psi_a \varphi_b h_{ca} s_{bc}^{de} \xi_e \eta_d &= \sum_{a,b,c,d,e} \psi_a \varphi_b h_{ca} s_{bc}^{de} \xi_e \eta_d \\ &= \sum_{a,b,c} \psi_a \varphi_b h_{ca} \sum_{de} s_{bc}^{de} \xi_e \eta_d \\ &= \! \sum_{a,b,c} \psi_a \varphi_b h_{ca} \! \left[ \! \begin{pmatrix} 1 \!&\! 0 \\ 0 \!&\! 0 \end{pmatrix} \! \eta_0 \xi_0 \!+\! \begin{pmatrix} 0 \!&\! 0 \\ 0 \!&\! 1 \end{pmatrix} \! \eta_0 \xi_1 \!+\! \begin{pmatrix} 0 \!&\! 0 \\ 1 \!&\! 0 \end{pmatrix} \! \eta_1 \xi_0 \!+\! \begin{pmatrix} 0 \!&\! 1 \\ 0 \!&\! 0 \end{pmatrix} \! \eta_1 \xi_1 \! \right] \\ &= \sum_{a,b,c} \psi_a \varphi_b h_{ca} \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} = \sum_{b,c} \varphi_b \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} \sum_{a} h_{ca} \psi_a = \frac{1}{\sqrt{2}} \end{align*}$$

And so on. As for why the CX gate ended up like this, I can't write it cleanly yet and just followed the calculation, but the result seems to be like calculating $\braket{\varphi| 0}\braket{1|H|0}$ ^[2]. I haven't followed it properly, but I think the tensor decomposition of CNOT in reference [BQ] appears as the 4th-order tensor $S$ written above, but I'll leave that as homework for now...

What I want to say here is that even complex quantum circuit calculations like entanglement have become just a calculation of tensor products.

Circuits with a Large Number of Qubits

In a previous article Playing with Qiskit (3), I concluded by saying:

This clarifies why simulating state vectors only works for a few dozen qubits. In this experiment, a VM with 26GB of memory was used, and it's suspected that it exceeded the limit at 31 qubits, which happened to be just barely manageable. This seems to be the limit in a typically available experimental environment.

This time, I want to challenge 200 qubits^[3]. With this many qubits, attempting to find probability amplitudes using the settings from the previous article would have completely failed.

n_qubits = 200

qr = QuantumRegister(n_qubits, 'q')
qc = QuantumCircuit(qr)
qc.h(qr[:]);

The circuit diagram is quite plain, with an H gate placed after each input qubit, so I'll omit it. Or rather, it's so long that I don't feel like taking a screenshot. And then,

converter = CircuitToEinsum(qc, dtype='complex128', backend='cupy')
bitstring = '0' * n_qubits
expr, operands = converter.amplitude(bitstring=bitstring)
amplitude = contract(expr, *operands)

...with this, the probability amplitude for the state $\ket{000\cdots 000}$ should have been calculated... probably.

amplitude = contract(expr, *operands)
print(np.isclose(amplitude, np.sqrt(1/2**n_qubits)))
print(amplitude)

True
(7.888609052209983e-31+0j)

And so, the probability amplitude $\sqrt{\frac{1}{2^{200}}}$ was successfully found.

Summary

I tried using cuQuantum, which allows for GPU-accelerated tensor calculations via tensor networks for quantum circuit simulations. While this doesn't explain why tensor networks are used or their specific benefits, since the scope of this article is limited to "using cuQuantum," I'll address those topics separately in the future... For now, I'll wrap up with: I'm glad I could calculate probability amplitudes for a 200-qubit circuit!

On a related note, one could perform measurements using Qiskit's MPS simulator based on the matrix product concept discussed in Thinking about Matrix Products and following the Matrix product state simulation method as shown below. However, since this would require something on the scale of $2^{200}$ shots, it seems impractical to derive probability amplitudes from there... but I wonder.

qc.measure_all()
simulator = AerSimulator(method='matrix_product_state')

tcirc = transpile(qc, simulator)
result = simulator.run(tcirc).result()
counts = result.get_counts(0)

References

[CQ] cuQuantum
[TN] Tensor Network
[BQ] Yuichiro Minato, Decomposition of CNOT Gate, 2022

脚注

1. Specifically, cuquantum.CircuitToEinsum is used.

2. When calculating the probability amplitudes of other states, the inner $\ket{0}\bra{1}$ part should change.

3. There's no deep meaning; it's just a number of qubits sufficiently larger than in the previous article.