Distributed Simulation of Statevectors and Density Matrices (Jones et al. 2023)

I. Introduction

A. Notation

Bit Extraction Notation

$i\left[t\right]$ represents the $t$ qubit of the binary basis $i$

$13={1101}_2=i$

bit:  3 2 1 0
      1 1 0 1

$i\left[2\right]=1$

Negated Bit Notation

If $i\left[t\right]=1$ represents the binary value of qubit $t$ in statevector $i$, then ! $i\left[t\right]$ represents flipping the value such that ! $i\left[t\right]=0$, and vice versa

Flipped Index Notation

If $i_t$ represents the statevector $i$ with $t$ unchanged, then $i_{\neg t}$ represents the same statevector with the $t$ bit flipped

Multiple-bit Flip Notation

Similar to flipped index notation, except $i_{\neg \{t_1,t_2\}}$ represents flipping both bit values in the binary string

Statevector Notation

For a $N$ qubit state

$|\psi ⟩={\sum }_i^{2^N}{\alpha }_i|i⟩$

Computational bases are represented using binary decimal notation. For a $3$ qubit state

${|0⟩}_3=|000⟩$, ${|1⟩}_3=|001⟩$, ${|2⟩}_3=|010⟩$, ${|3⟩}_3=|011⟩,\cdots ,{|7⟩}_3=|111⟩$

And it is stored as

${|\psi ⟩}_N=[{\alpha }_0,{\alpha }_1,\cdots ,{\alpha }_{2^N-1}]$

This is the statevector, where each index corresponds to the basis state

Prefix & Suffix Bits

Consider a node to be "one computer/machine participating in the simulation. Distributed simulation defines $W$ as the number of nodes being used, which assumes

$W=2^w$ (# of nodes available is a power of $2$)

Every amplitude index $i$ has $N$ bits total:

• Some of those bits identify the node

• The remaining bits identify the amplitudes position within that node

Suppose $N=6$. A statevector then would look like

q6 q5 q4 q3 q2 q1 q0

Now suppose $W=8$. That is, we have $8$ nodes to work with. Since $W=2^3$, that means we need $3$ bits to identify the node, i.e. the prefix bits. These prefix bits are defined as the rank bits

r1 r0 q4 q3 q2 q1 q0

The remaining suffix bits are local index bits, defined as $\lambda =N-w$. The key difference is

• Suffix/local bits can be manipulated entirely within one node

• Prefix/rank bits requires communication with other nodes, usually reassigning the statevector to a new node

global index $=$ (prefix bits)(suffix bits)

Moreover, global amplitude index $i$ can be decomposed into a combination of the node rank $r$ bit-shifted to the left by $\lambda$ and the local-index $j$

i = (r << lambda) + j

II. Local Statevector Algorithms

A. Local Costs

"Before we can discuss the challenges of distributed simulation, we must first understand the simpler problem of local (i.e. non-distributed), shared memory simulation." (Section II, first paragraph)

Recall that the simulation stores a full-state vector $|\psi ⟩={\sum }_i^{2^N}{\alpha }_i|i⟩$ as an array

$\psi =\{{\alpha }_0,{\alpha }_1,\cdots ,{\alpha }_{2^N-1}\}$

This notation gives a complete description of a quantum state, but the memory costs grow exponentially as $N$ increases

• Considering C/C$++$ stores complex doubles as $16$ bytes, for $N$ qubits we would need

$2^N$ states $\times$ $16$ bytes

• Dense mixed states are even worse

$2^N$ states $\times$ $2^N$ states $\times$ $16$ bytes

Full statevector simulators are useful for researchers because resource utilization and runtime costs are roughly the same across circuits of equal size, irregardless of state properties, making resource prediction trivial

But there are two main drawbacks of full state simulation, which is that it requires simultaneously

1. Storing $2^N$ complex doubles $\left(\alpha \right)$

2. Simulating $n$ qubit operators, which required $O\left(2^{N+n}\right)$ floating point operations

Thus, the authors use four local performance metrics to measure these costs

Basic Operations bops

Bitwise operations, integer arithmetic, logical checks, indexing, etc.

i ^ (1ULL << targ)
(i >> targ) & 1
if (ctrlBit == 1)

Floating Point Operations flops

Arithmetic operations on real or complex values

m00*beta + m01*gamma;
m10*beta + m11*gamma;

Memory memory

Memory overhead, i.e. the size of temporary data structures created during an algorithm

• This excludes persistent pre-allocated memory dedicated to the quantum state representation

  • one-target gate: $O\left(1\right)$ extra memory

  • many-target gate: $O\left(2^n\right)$ temporary array

Writes writes

The number of writes to heap/main memory, especially writes to statevectors

• Moving amplitudes to and from memory is more often the bottleneck as opposed to performing arithmetic operations

• Thus if an algorithm writes to fewer amplitudes, it may be much faster than other algorithms even if the math looks the same

Cost tags are summarized as:

[a bops][b flops][c memory][d writes]

B. One Target Gate

Let ${\hat{M}}_t$denote a single qubit gate operating on qubit $t$, described by a complex matrix

${\hat{M}}_t=\left[\begin{array}{cc}{m}_{00}& m_{01}\\ m_{10}& m_{11}\end{array}\right]$

${\hat{M}}_t$ gets applied to the target qubit $t$ as follows

${\hat{M}}_t|\psi ⟩={\sum }_i^{2^N}({\alpha }_i|i_{N-1}⟩\cdots {\hat{M}}_t|i_t⟩\cdots |i_0⟩)$

${\hat{M}}_t|i_t⟩=\{\begin{array}{ccc}{m}_{00}|0⟩+m_{01}|1⟩& & i_t=0\\ m_{10}|0⟩+m_{11}|1⟩& & i_t=1\end{array}$

${\hat{M}}_t|\psi ⟩={\sum }_i^{2^N}\left({\alpha }_i\right)|i_{N-1}⟩\cdots (m_{i_t{i}_t}|i_t⟩+m_{j_t{i}_t}|j_t⟩)\cdots |i_0⟩$

${\hat{M}}_t|\psi ⟩={\sum }_i^{2^N}({\alpha }_i{m}_{i_t{i}_t})|i⟩+{\alpha }_i{m}_{j_t{i}_t}|i_{\neg t}⟩$

This final equations reveals that applying a one target gate operator to a full-statevector linearly combines pairs of amplitudes. Thus, when applying the operator ${\hat{M}}_t$, the amplitudes of $|\psi ⟩$ are simultaneously modified to become

${\alpha }_i\to {\alpha }_i{m}_{i_t{i}_t}+{\alpha }_{i_{\neg t}}{m}_{i_t{j}_t}$

Since the pairs of amplitudes are $({\alpha }_i,{\alpha }_{i_{\neg t}})$, and since $i_{\neg t}=i\pm 2^t$, paired amplitudes are always a distance $2^t$ apart