Quantum Circuit Simulation with _qsim_

• qsim is a high performance simulator of quantum circuits for multinode quantum trajectory simulations
• qsim can be used as a backend of Cirq, a Python software library for writing quantum circuits
• A delayed inner product algorithm for quantum trajectories can result in an order of magnitude speedup for low noise simulation
• This framework is usable in Google Cloud Platform, with high performance virtual machines in a single node or multinode setting
• Multinode configurations simulate noisy quantum circuits with quantum trajectories well
• An approximate noise model for Google's experimental quantum computing platform can be used and the results of noisy simulations can be compared with experiments for several quantum algorithms on Google's Quantum Computing Service

Introduction

• Classical software simulates quantum circuits with an approximate noise model of quantum hardware which enables the study of NISQ quantum algorithms and applications discovery
• Google Quantum AI recently launched qsim to allow users of its open source ecosystem to simulate quantum circuits more efficiently on classical processors
• Software tools include Cirq, a quantum programming framework, ReCirq, a repository of research examples, and application-specific libraries like OpenFermion and TensorFlow Quantum
• New features made quantum circuit simulations more performant and intuitive, and noise simulations more sophisticated
• The theory and software routines which underpin qsim's performance are described
• qsim implementations and workflows for various classical processor types and setups, including single and multinode CPU and GPU setups are outlined
• A generic noise model approximates Google's Quantum Computing Service (QCS)

Cirq: a programming framework for quantum circuits

• Cirq is a Python software library for writing, manipulating, optimizing and running quantum circuits on quantum computers and quantum simulators
• Cirq can be used with experimental quantum processors, such as Google's Quantum Computing Service, Alpine, Pasqal, Rigetti and IonQ
• It comes with built-in Python simulators for testing small circuits, and supports high performance simulators, such as Qulacs and quimb
• Cirq has also been integrated with other software libraries, such as QC Ware Forge, Xanadu Pennylane, Zapata Orquestra, Sandia National Lab pyGSTi, CQC t ket> and Quantum Benchmark True-Q
• Cirq is part of Google Quantum AI open source ecosystem, which includes ReCirq, OpenFermion and TensorFlow Quantum
• Focus is on the use of Cirq to simulate approximate experimental noise with qsim
• A Qubit in Cirq is an abstract object that has an identifier, its state is maintained in a quantum processor or a simulator
• A Gate in Cirq is an effect that can be applied to a collection of qubits and can be a unitary gate or a quantum channel
• Quantum channels can represent noise, such as amplitude or phase damping channels
• The primary representation of quantum programs in Cirq is the Circuit class
• A Circuit is a collection of Moments. Each Moment is a collection of Operations that act during the same time slice, but in different qubits.
• An Operation in Cirq is a Gate that has been applied to qubits
• Cirq includes tools to transform circuits, including adding quantum channels after unitary gates in a circuit to simulate noisy experimental quantum processors
• Noise can be added to a Cirq circuit before constructing a simulator object and simulating the circuit
• There are two procedures for adding noise to a Cirq circuit: adding individual noise events, or defining a global noise model
• Add an individual noise event as a cirq . Channel with corresponding noise parameters to the Cirq circuit (in the cirq. Circuit argument)
• Calling the cirq.kraus protocol on a channel returns the Kraus operators corresponding to that channel
• All channels are subclasses of cirq. Gate, they can act on qubits and be used in circuits
• Cirq has multiple common noise channel options built in, such as the depolarizing channel cirq.depolarize, the phase damping channel cirq.phase_damp and the bit flip channel cirq.bit_flip
• Channels can be controlled by appending .controlled. Custom channels can be defined using MixedUnitaryChannel or KrausChannel
• MixedUnitaryChannel takes a list of (probability, unitary) tuples and uses it to define the mixture_ method
• KrausChannel takes a list of Kraus operators and uses it to define the _kraus method
• A measurement key can be used as a parameter in a custom noise channel and store the index of the selected unitary or Kraus operator in the measurement results
• Noise that affects an entire circuit can be described with the cirq. NoiseModel type
• Objects of this type must define one of three methods to describe how to convert a "clean" circuit into a "noisy" circuit: (1) noisy_operation, which mutates each operation independently, (2) noisy_moment, which mutates each set of simultaneous operations, or “moment”, as a group, or (3) noisy_moments, which mutates the entire circuit at once
• A simple version is provided in cirq. Constant QubitNoiseModel, which applies a specified gate or channel to every qubit in the circuit at the start of each moment
• Implement own NoiseModel type for more complex behavior
• Once constructed, a NoiseModel can be applied to a circuit using the Circuit.with_noise method
• This generates a "noisy" version of the original circuit, which can then be simulated with qsim or one of the builtin Cirq simulators

qsim: a quantum circuit simulator

• qsim is a full state vector quantum circuit simulator that computes all the 2n amplitudes of the state-vector, where n is the number of qubits
• Gates or operators are applied to the state vector, the simulator performs matrix-vector multiplications repeatedly
• Single precision arithmetic and gate fusion are used to speed up the simulation
• SIMD instructions (single instruction/multiple data) are used for vectorization and OpenMP for multi-threading on CPUs
• SIMD versions are available: SSE, AVX2/FMA, and AVX512. CUDA is used in the GPU implementation

Matrix-vector multiplication

• For a q-qubit gate and an n-qubit state vector, the full matrix-vector multiplication can be block diagonalized into 2n-q (gate matrix)-subvector multiplications, where the gate matrix is of size 2q × 2q and each subvector is of size 2q
• The total number of flops is = 29+2 · 2n+1 and the total number of bytes to read and write is 2n+4 (single precision)
• The arithmetic intensity (the ratio of the number of flops to the number of bytes to read and write) is 29-1
• The arithmetic intensity is small for small values of q, performance is usually limited by the memory bandwidth
• Fusion of small gates into larger gates increases the arithmetic intensity and better utilizes compute power of modern CPUs and GPUs

Gate fusion

• A quantum circuit can be considered as a lattice structure that has spatial and time directions.
• The time direction corresponds to the order in which gates are applied.
• Gate fusion combines gates that are close in space and time into larger gates
• Gate fusion increases the arithmetic intensity, decreases the number of gates, and leads to a significant speedup
• There are two steps in the fusion algorithm: First, large gates and small gates that are neighbors in time and act on the same qubits are combined. Second, the algorithm greedily combines gates that are close in space and time
• Essentially, the resulting gates from the first step are unmarked. The unmarked gates are picked for processing in increasing time order.
• The first unmarked gate is picked and marked.
• The nearest (unmarked) neighbors (the gates that share qubits with the picked gate) forth in time and the next nearest unmarked neighbors back in time (if they do not have unmarked neighbors further back in time) are added while the resulting fused gate is no greater than the specified maximum fuse size f
• All added gates are marked/ This procedure is repeated until all the unmarked gates are exhausted
• The optimal value of the maximum fuse size f is 4 for large numbers of threads (or on GPUs) and large circuits
• Smaller fuse size, f = 2 or f = 3, can be optimal for small numbers of threads and/or small circuits

SIMD implementation

• SIMD instructions (single instruction/multiple data) are used to make the most of the compute power of CPUs
• Avoid the usage of horizontal SIMD instructions by keeping the real and imaginary parts of k state-vector amplitudes in separate SIMD registers, where k is the SIMD register size in floats
• In single precision, k = 4 for SSE, k = 8 for AVX, and k = 16 for AVX512
• Separate blocks of szie k can be optionally stored in memory
• Parallel matrix-vector multiplications of up to k subvectors
• “high” qubit indices are qubit indices that are larger than or equal to log2(k) and the qubit indices that are smaller than log2(k) as “low” qubit indices
• Make use of SIMD instructions if all the gate qubits are high with first loop runs from 0 to 2n-q-log2(k) 1 and k2q state-vector amplitudes are loaded into 2·2q SIMD registers, and SIMD arithmetic instructions calculate k matrix-vector products simultaneously
• For some low gate qubits then the matrix elements to calculate the sum in line 6 in Algorithm 1 are loaded into SIMD registers accordingly

GPU implementation

• GPU implementation is very similar to the SIMD implementation and uses CUDA
• 32-thread warps are used instead of SIMD registers and SIMD instructions, k = 32
• A single thread in a warp performs the same role as a single data element in a SIMD register
• k/2l matrix-vector products are calculated in parallel by a single warp, where l is the number of low qubits
• The implementation efficiently utilizes the GPU compute resources and memory bandwidth

Quantum trajectories in qsim

• qsim supports noisy circuit simulations using quantum trajectories, implemented by choosing one Kraus operator for each quantum channel with Kraus operators {K}
• The probability to sample the Kraus operator Ki is pi = (Ψ|KiKi|Ψ)
• Probabilities pi sum to unity for each channel, ∑pi = 1
• A Kraus operator is sampled and applied to the state vector for each channel sequentially
• This procedure is typically repeated many times, once per quantum trajectory
• The Monte Carlo statistical error for an observable estimated with quantum trajectories goes like 1/√r, where r is the number of trajectories
• The number of trajectories is typically in the thousands or higher
• A delayed inner product algorithm for quantum trajectories can result in an order of magnitude speedup for low noise simulation
• In the conventional quantum trajectory algorithm, at least one Kraus operator is applied immediately for each channel.
• Improved delayed inner product algorithm uses a lower bound pi for each sampling probability pi, given by the smallest singular value (squared) of the operator matrix Ki
• Note that in general, the bounds pi sum to a value s that is smaller than unity
• Sample the Kraus operator by drawing a random number r from the range [0,1). If r < s then there is no need to apply any Kraus operator immediately and we avoid computing any inner product ⟨Ψ|KiKi|Ψ⟩
• The operator can be sampled just by using the lower bounds, and the application of the picked operator can be deferred
• Can make use of gate fusion by deferring the operator application
• First, fuse and apply all the operators that were deferred in the previous steps
• Second, use the conventional sampling procedure
• Expectation value can be calculated in place without copying the state vector to a temporary vector
• This reduces the memory usage
• In the case of weak noise, the sum s of lower bounds is typically close to one and the operators get deferred with a high probability, giving rise to a significant speedup
• The runtime is linear in noise strength for typical noise values
• The runtime of the conventional algorithm is weakly dependent on the noise probability, so the runtime at a noise probability of 0.1 gives a lower bound for the runtime of the conventional algorithm
• If all Kraus operators in a channel are proportional to unitary matrices then s = 1 and we always defer the application of such a channel

Main qsim runtime factors

• The main factor in the runtime of a circuit simulation is the number of qubits
• The size of the state vector for n qubits is 2n, and therefore the runtime is also exponential in the number of qubits
• For best performance, threads should equal the number of cores in the machine
• If the maximum number of threads is not used on multi-socket machines, then it is advisable to distribute threads evenly to all sockets or to run all threads within a single socket
• The runtime is linear in the circuit depth, as the number of matrix-vector multiplications is linear in the depth
• The runtime increases linearly with the noise strength for a quantum trajectory
• The performance does not depend on the noise strength if all Kraus operators in the quantum channel are proportional to unitary matrices

Simulating quantum circuits on the Google Cloud Platform

• First consideration when choosing hardware is it's the memory required to simulate the circuit
• The memory required to simulate an n qubit circuit is 8 · 2n bytes
• The maximum number of qubits that can be simulated on a given machine is limited by its RAM memory
• The maximum is currently 32 qubits in a Google cloud GPU (on an NVIDIA A100 GPU with 40GB of memory), and 40 qubits on a virtual machine (on an m2-ultramem-416)