Quantum Supremacy Is Both Closer and Farther than It Appears
Igor L. Markov Aneeqa Fatima Sergei V. Isakov Sergio Boixo




Quantum Supremacy Is Both Closer and Farther 

        than It Appears
    Igor L. Markov^1, Aneeqa Fatima^1, Sergei V. Isakov^2, and Sergio Boixo^3  

            ^1 University of Michigan,  2260 Hayward St, Ann Arbor, MI 48109 

            ^2 Google Inc., 8002 Zürich, Switzerland

            ^3 Google Inc., Venice, CA, 90291
            
    
================================================================================================================================================================================================================================================================================




As quantum computers improve in the number of qubits and fidelity, the question of when they surpass
state-of-the-art classical computation for a well-defined computational task is attracting much attention.
The leading candidate task for this 
milestone entails sampling from the output distribution defined by a random quantum circuit. We develop
a massively-parallel simulation tool Rollright that does not require inter-process communication (IPC) or proprietary hardware. We also develop two ways to trade circuit fidelity for computational speedups, so as to match the fidelity of a given quantum computer — a task previously thought impossible. We report massive speedups for the sampling task over prior software from Microsoft, IBM, Alibaba and Google, as well as supercomputer and GPU-based simulations. By using publicly available Google Cloud Computing, we price such simulations and enable comparisons by total cost across hardware platforms. We simulate approximate sampling from the output of a circuit with 7 × 8 qubits and depth 1+40+1 by producing one million bitstring probabilities with fidelity 0.5%, at an estimated cost of $35184. The simulation costs scale linearly with fidelity, and using this scaling 
we estimate that extending circuit depth to 1+48+1 increases costs to one million dollars.
Scaling the simulation to 10M bitstring probabilities needed for sampling 1M bitstrings helps comparing simulation to quantum computers. We describe refinements in benchmarks that slow down leading simulators, halving the circuit depth that can be simulated within the same time.
 



§ INTRODUCTION


 The promise of quantum computers inspired the pursuit of quantum computational supremacy, i.e.,
 using quantum computers to solve some task that is prohibitively hard for conventional computers <cit.>.
A key benchmark has emerged in terms of random, universal quantum circuits proposed by Google <cit.>,
and this inspired rapid progress in algorithms for simulating quantum computers <cit.>. Simulating sampling or performing cross-entropy benchmarking <cit.> requires calculating a large number of output probabilities requested at random. We focus on a method that essentially evolves the quantum wave function and outputs a large number of probabilities with small additional cost <cit.>. We show that two recent quantum simulations that required supercomputers with massive amounts of memory <cit.> can now be performed with modest amounts of space and time using publicly accessible cloud computing.
To match the accuracy of results in a given quantum computer with arbitrary quantum gates,
we developed an approximate simulation whose runtime scales linearly with circuit fidelity (a measure of 
how much quantum information is preserved).  In contrast to the recent use of the world's most powerful supercomputer (over 100 Pflops, 131K nodes) to simulate Google circuits <cit.>, we add a new dimension to quantum-versus-classical comparisons and report monetary cost of our simulations in terms of cloud-computing resources. This metric captures rapid progress in quantum-circuit simulation since 2016. On the other hand, we show that


  * Small changes to common quantum-supremacy circuits make them considerably harder to simulate.

  * The use of more sophisticated quantum gates substantially handicaps leading simulation methods.


In 2017 and 2018, IBM, Intel and Google announced new quantum computing chips that implement 50 and 72 qubits <cit.>, but no empirical results for these chips have yet been announced as the chips are being improved. These chips will have to compete with simulation software running on conventional computers. The most straightforward approach — Schrödinger-style simulation that maintains and modifies the wave-function — requires 2^n-sized memory for n qubits <cit.>. Therefore, it is possible to simulate 30-qubit circuits on a laptop, and 35-qubit circuits on a mid-range server. Further scaling would seem to require supercomputers, and a 2016 simulation <cit.> used 0.5 PB RAM on a Cori II supercomputer to simulate a depth-26 45-qubit circuit proposed by Google for quantum supremacy experiments <cit.>. These universal circuits were designed to run on planar √(n)×√(n) qubit-array architectures <cit.> and evade easy simulation.  They start with a Hadamard (H) gate on every qubit and arrange nearest-neighbor controlled-Z (CZ) gates in a repeated pattern, so that every possible CZ gate appears once every eight cycles. One-qubit , and T gates are randomly interspersed between CZ gates so as to prevent cancellations and induce chaotic quantum dynamics <cit.>, as illustrated in <cit.>. Each circuit ends with a measurement on every qubit. On a quantum computer this measurement produces an unpredictable bitstring and thus samples the distribution determined by the specific quantum circuit. Competing simulation techniques must calculate a large enough number of amplitudes (each corresponding to some bitstring) to simulate sampling from the output distribution <cit.>.

Common sense suggests that circuits with more than 49 qubits would require too much memory to simulate in practice. To the contrary, alternative methods related to
Feynman paths <cit.>, tensor network contractions <cit.>, and similar
approaches <cit.> can trade space for time complexity. Google researchers have shown that for circuits with low depth, one can quickly find any one amplitude on a single computer for over 100 qubits using a variable elimination algorithm <cit.>, and have illustrated the statistical distribution of the output of a circuit with 7 × 8 qubits and depth 30 by computing 2× 10^5 probabilities using multiple computers. 
In April 2018, a simulation on the world's largest supercomputer, Sunway TaihuLight <cit.>, computed 2^46 amplitudes for a circuit with 7 × 7 qubits and depth 1+39. Our notation for depth (introduced here) explicitly denotes layers of Hadamard gates with “1+” at the beginning of the circuit (clock cycle 0) and “+1” at the end. The remaining 39 layers (clock cycles) include CZ gates and other, one-qubit, gates. A related computation produced a single amplitude for a circuit of depth 55. The depth attained in these simulations includes 8 cycles gained by exploiting an unfortunate design choice — sequences of diagonal gates , — at the beginning of the circuit, specific to the circuits in Ref. <cit.> and avoided in revised benchmarks <cit.> as explained in Section
<ref>.

All simulations discussed so far seek exact output amplitudes, aside from negligible numerical errors. However, the computational task of interest is to approximately sample from the bitstring distribution defined by a random quantum circuit, and near-term quantum computers incur significant, unavoidable errors because gate errors accumulate exponentially <cit.>. In particular, for Google quantum supremacy circuits <cit.> with 7×7 qubits and depth 1+40+1, a reasonable goal is to achieve a 0.005 circuit fidelity, consistent with a two-qubit gate fidelity of 0.995, one-qubit gate fidelity of 0.999, initialization fidelity of 0.998 and measurement fidelity of 0.99 <cit.> (in this example, depth 1+40+1 entails initial and final Hadamard gates on each qubit). The same circuit fidelity at depth 1+48+1 would require a two-qubit gate fidelity of 0.996. In our simulations, we configure circuit fidelity values f=0.01 and f=0.005, and show how to revise resource estimates to any given f∈ (0, 1], using linear scaling that we prove. Significantly, rather than simulate individual gates to fixed fidelity, we control simulation fidelity for the entire circuit.
 
Extending recent progress in quantum circuit simulation algorithms, coauthors from the University of 
Michigan implemented Rollright — a massively-parallel simulation tool that does not require inter-process communication (IPC) or proprietary hardware. Rollright can run in commercial cloud-based preemptible virtual machines (VMs).[A preemptible VM is cheaper than an on-demand VM because it does not guarantee real-time execution. A process may be terminated (preempted) by a higher-priority task at any time, but can then be restarted later. Pricing for preemptible VMs in a commercial cloud varies for different hardware configurations and may fluctuate hourly with demand for computational resources.] Table <ref> summarizes advantages of this approach over supercomputer simulations <cit.>.
 We report the following results:


 
  * Using a MacBook Pro laptop with 16GiB RAM, Rollright can simulate circuits with 5 × 5 qubits
         to any depth, with 16× and 10× speedups over the latest versions of simulators from
         Microsoft QDK (Appendix <ref>) and IBM QISKit / QASM (Appendix <ref>),
         respectively.  Rollright uses 3.27× and 1.7× less memory respectively and can also
         simulate 6×5-qubit circuits of any depth on the MacBook. On a single 72-thread server,
         Rollright in its depth-limited mode simulates a depth-28 circuit with 6×5 qubits 31×
         and 16× faster than Microsoft and IBM simulators, respectively. Rollright uses 886×
         and 590× less memory for this simulation.[These results were presented at Microsoft and IBM.]
         Running QISKit on GPU-based servers through the IBM Q service is 12× slower than Rollright,
         which uses no GPUs. The speedup increases to 877× for a 32-qubit circuit.

  * Using a single n1-highcpu-96 server
       from Google Cloud, we calculate one million probabilities (requested at random) per simulation, enough to simulate sampling <cit.> for circuits used in two recent supercomputer simulations:
       
       
  * For a circuit with 6 × 7 qubits and depth 1+25 <cit.>,
       our simulation took 4.7 hours, with estimated cost $3.34 in the cloud.
       This is significantly cheaper and more accessible than simulations in Ref. <cit.>,
       but we recall
       that those simulations study the convergence to the Porter-Thomas distribution as a function of depth
       with exponential precision.
       
  * For a circuit with 9 × 5 qubits and depth 1+25 <cit.> our simulation took 20 minutes,
       with estimated cost $0.24 in the cloud.
       We use 17.4GiB peak memory where Ref. <cit.> used 0.5PB — a 28600 times
       improvement.[Our peak memory usage can be improved or traded for runtime at the same
       cost point.]  Notably, for methods that compute one amplitude at a time,
       9 × 5 circuits are not substantially harder than 
       5 × 5 circuits <cit.>.
   
   These results are reported only to compare our simulation to prior work. We also calculate one million probabilities for the circuit  <cit.> with 7 × 6 qubits and depth 1+25+1. This refined benchmark is substantially harder to simulate for some methods (see Section <ref>). The runtime was 6 hours (as apposed to 4.7 hours for the previous circuit of similar size), corresponding to a cost of $4.32.


 
  * We propose two complementary new approaches to approximate simulation of quantum circuits
         with arbitrary gate libraries (Sections <ref> and <ref>):
         
         
  * When simulating a quantum circuit, we obtain arbitrary requested sets of output
         amplitudes with prescribed fidelity.  Our ability to produce
         millions of amplitudes in one run implies massive speed-ups over prior works that
         produce one amplitude at a time <cit.>.
         
  * Given a set of output probabilities, drawn at random, we produce a bitstring
         with a tunable statistical error. Ten randomly chosen probabilities suffice
         to produce one bitstring sample with negligible error. 
         
         Numerically, we compute one million probabilities with fidelity 0.5%  for the
          circuit <cit.> with
         7 × 7 qubits, depth 1+40+1 and controlled-Z (CZ) two-qubit gates. For this,
         we use 625 n1-highcpu-32 servers from Google Cloud, with estimated cost $8734.
         We study three types of scaling for resource requirements
         
         
  * Linear scaling with fidelity f.
         
  * Nonlinear scaling with the number of amplitudes (Table <ref>).
         
  * Exponential scaling with circuit depth. 
         Extending circuit depth to 1+48+1 increases cost estimates to one million US dollars. 
         

  * We produce as many bitstring samples as required
         for cross-entropy benchmarking with the fidelities chosen, enabling cost and runtime comparisons between simulators
         and quantum computers.

  * We propose an interactive protocol to validate the results of quantum-supremacy
      simulations without knowing the correct answer (see Appendix <ref>).
      The protocol uses (i) approximate simulation with prescribed fidelity, and
      (ii) inner-product estimation.


Our simulation produces a large number of requested amplitudes at once (cf. <cit.>) and is not specialized to individual circuits (cf. <cit.>). In Section <ref>, we summarize weaknesses found in quantum supremacy circuits <cit.> simulated in prior work and introduce new benchmarks that are harder to simulate. Section <ref> shows that replacing CZ gates with iSWAP gates makes simulation much more difficult yet.







§ OUR QUANTUM-CIRCUIT SIMULATION FRAMEWORK


Our algorithms can be characterized as Schrödinger-Feynman hybrids <cit.>. In the context of nearest-neighbor quantum architectures, we partition a given qubit layout into blocks. We then decompose quantum gates acting across the partition into sums of separable terms, such that, for each term, each block can be simulated independently and the results can be added up. For example, using a pair of half-sized qubit blocks instead of representing a full wave function reduces memory requirements for k qubits from 2^k to 2^k/2 + 1, but introduces a dependency on the number of decomposed gates. For CZ gates the decomposition has two terms[A two-qubit gate can be decomposed into at most four tensor products, e.g., with one-side operator basis _0, _1, |0⟩⟨ 1|, |1⟩⟨ 0|.]

    = diag(1,1,1,-1) = (
       [ 1 0; 0 0;   ]) ⊗(
       [ 1 0; 0 1;   ])
       +
       (
       [ 0 0; 0 1;   ])⊗(
       [  1  0;  0 -1;    ])
       = _0⊗ + _1⊗ .

Thus, applying each cross-block CZ gate (xCZ) to a tensor term produces two tensor product terms, doubling runtime.[Moreover, new pairs of tensor terms are orthogonal because their first components (produced by P_0 and P_1) are. While unitary operators preserve their orthogonality, projections can increase their cosine similarity.] In comparison, traditional Feynman-style path summation <cit.> uses very small amounts of memory, but doubles its runtime on every (branching) gate, resulting in much longer
runtime and not being able to use available memory fully. Our simulator combines highly-optimized Schrödinger-style simulation within each qubit block and simulates xCZ gates with Feynman-style path summation, to limit memory use. Unlike in Feynman-style simulation, runtime scales with the number of
xCZ gates, which is very limited in planar qubit-array architectures with nearest-neighbor gates <cit.>. Unlike traditional Schrödinger-style simulation, the resulting algorithms are depth-limited, and supercomputer simulations may hold some advantage for very deep circuits. However, near-term quantum computers rely on noisy gates <cit.> that also limit circuit depth.

 

Our Schrödinger-style simulation includes
optimizations that help with arbitrary gate libraries and some
that help with small gate libraries, while others target the Clifford+T and related gate
 libraries[Clifford+T and related gate libraries are common because they promise compatibility with
quantum error-correction, they are used in Google and IBM quantum computers.] These optimizations
equally apply when using supercomputers, laptops and other hardware.


  * Clustering gates of a kind (with reordering), rather than gates acting on the same qubits as in prior simulations. For example, in Google quantum supremacy circuits <cit.>, we collect separate clusters of CZ, T, and gates.[This technique might not help in circuits that use a large number of different gate types.]

  * Compact encoding of gate positions in each cluster into bitmasks, so that efficient bitwise operations (parity, population count, as well as counts of trailing and leading zero bits) and mod-8 arithmetics (for T gates) can apply all diagonal gates in a single pass over the wave-function or its slices.[Some of these techniques exploit specific gate types, and our simulator may run 2-3× slower on a different gate library.]

  * Cache-efficient algorithms to simulate large gate clusters. Instead of accessing
an entire wave-function, we apply gate clusters on amplitude slices sized
to fit in the L2 cache.[This technique does not depend on specific gate types.]

Compared to the state of the art described in <cit.>, we can eliminate floating-point multiplications when simulating Google circuits (and only use them for convenience when simulating T gates). The number of floating-point additions is reduced by a factor of three, while fully benefiting from vectorized instruction extensions (AVX2) in modern CPUs. This puts emphasis on cache performance and
limits thread scalability. By streamlining memory accesses, we improve single-thread performance as well
as scaling to multiple threads and multiple processes using the same memory bank. However, the impact of subsequent improvements is even greater.

 

For Feynman-style path summation, our simulator optimizes performance in several ways.
When simulating a circuit with x>0 xCZ gates, we use the technique suggested in <cit.> and demonstrated in <cit.> where each process receives a bitstring that encodes a unique path, then save the simulation results in a file for subsequent summation.[Our simulator was implemented before <cit.> was described publicly.] This is convenient because no interprocess communication is required, but many processes simulate the same path prefixes. We therefore subdivide x=x_p + x_b and use only x_p-bit prefixes to spawn separate processes. After the first x_p bits, each process checkpoints its simulation state (doubling its memory usage once) and simulates each of 2^x_b remaining paths starting from this checkpoint. A second optimization for Feynman-style path summation is supported by the Schrödinger-style simulator algorithm: we noticed that the use of projection operators in gate decompositions such as Equation (<ref>) leads to a large number of zeros in block-local wave-functions. Our simulator skips blocks of zero amplitudes.

 

Final-state amplitudes are calculated for any requested subset of indices, unlike simulators in
<cit.>. This is necessary to simulate sampling <cit.> and to perform cross-entropy benchmarking <cit.>. Subsets of millions of amplitudes are calculated almost as quickly as single amplitudes (c.f. <cit.>). However, if a very large subset is requested, such calculations may dominate the cost of “easy” simulations. 
Contributions to requested amplitudes are accumulated over simulation paths. When used in its exact mode, the simulator calculates each amplitude with a negligibly small numerical error.[As a check, we calculate the norm of the state vector using higher-precision accumulator variables and underflow mitigation. Comparing the result to 1.0, we observe a very close match in all cases when such a match is expected.]



§ TRADING OFF SIMULATION FIDELITY FOR COMPUTATIONAL RESOURCES


We are interested in performing the following computational task proposed for a near-term quantum supremacy demonstration: approximately sampling from the output distribution defined by a random quantum circuit <cit.>. Quantum computers experience
errors in qubit initialization, gates and measurements. This suggests speeding up the simulation by reducing the accuracy of results to the level attained by quantum computers. We are not aware of such attempts for simulating quantum supremacy circuits, but the work in <cit.> developed approximate simulation of circuits dominated by Clifford gates.

 

Our approach to approximate simulation
can handle arbitrary quantum gates, and our fidelity scaling is very different from error scaling in Refs. <cit.>, as shown below. Equation (<ref>) shows that each xCZ results in two tensor-product terms (up to four terms for other gates). Simulating x such xCZ gates (among other gates) produces 2^x simulation paths, corresponding to all possible choices of either _0 ⊗ I or _1 ⊗ for each xCZ, see Equation (<ref>) . Therefore we can write the output of the simulation as

    |ψ⟩ = ∑_j=0^2^x -1 |φ_j⟩ ,

where j enumerates the different paths.
To approximate this sum, one can drop some of the terms. In the context of multiprocess simulation,
one can simply skip some of the processes.
When simulating quantum-supremacy circuits <cit.>, the 2^x terms tend
to have nearly-identical norms because odd amplitudes and even amplitudes selected by projections in Equation (<ref>) are equally distributed. This is the worst case (otherwise, we would have dropped
terms with lower norms), so we assume |φ_j ⟩ = 2^-x/2  ∀ j. Furthermore, each path corresponds to a different quantum random circuit, and therefore their output states are almost orthogonal, |⟨φ_j | φ_k ⟩|^2 ≃ 2^- x · n for j  k.
 For requested simulation fidelity 0<f≤1, we include f 2^x terms out of 2^x
 
    |ψ_a⟩ = ∑_j=0^f 2^x -1  |φ_j⟩,
       so  ⟨ψ_a |ψ_a⟩= f 2^-x(see Figure <ref>) .

Fidelity can then be estimated as

    ⟨ψ |ψ_a⟩⟨ψ_a |ψ⟩/⟨ψ_a |ψ_a⟩ = f^2 2^-x/f 2^-x = f ,

see Figure <ref>.
Attaining fidelity f with only an f fraction of work is remarkable, e.g.,
quantum simulation with fidelity 0.1 is 10× faster than exact simulation.
This scaling holds up well in multiprocess simulation, as we confirmed empirically.
The specific selection as to which
 |φ_j⟩ terms to leave in Equation (<ref>) does not matter
for Google quantum supremacy circuits because the norms |φ_j ⟩ are nearly identical (we checked this). In general, dropping terms with smaller norms can yield a better-than-linear
tradeoff between fidelity and runtime.

Approximate simulation methods in <cit.> do not  appear particularly promising on quantum-supremacy circuits <cit.>, which have numerous T gates (see Tables <ref> and <ref>)
. The results in <cit.> differ from ours in another essential way. Their runtime scales with 1/ϵ^2 for circuit error-rate ϵ>0, which is beneficial for small-error simulation. Our scaling, linear in f, is beneficial for small-fidelity simulation and
 near-term quantum-supremacy experiments. 
 Cloud-based approximate simulation is discussed in Section <ref>.

 

Any cross-block multiqubit gate can be handled with our techniques by replacing the decomposition in Equation (<ref>) by an operator Schmidt decomposition of the gate, ensuring the minimal number of terms. The linear scaling of runtime with fidelity does not depend on this Schmidt rank.





§ HOW MANY AMPLITUDES MUST BE CALCULATED?



The computational task for quantum supremacy explored in this paper entails sampling from a distribution
of bistrings that is defined by the output of a given random quantum circuit <cit.>. 
To perform this task, we now leverage our simulation techniques that provide the output probabilities {p(x_j)} for any chosen set of bistrings {x_j} (see Table <ref>). Our frugal rejection sampling yields ℓ bistrings from the appropriate distribution (on average) using only 10 ℓ output probabilities. For the fidelity considered here, a quantum computer must sample ∼ 10^6 bitstrings to enable cross-entropy benchmarking <cit.>, so a simulation must find ∼ 10^7 amplitudes/probabilities. 

Basic rejection sampling.
A quantum computer running an n-qubit circuit from <cit.> may yield any one of N=2^n bitstrings on the output. However, outcome probabilities {p(x_j)} for sufficiently deep circuits <cit.> follow an exponential (Porter-Thomas) form 1/N∑_x δ (p(x) - p ) = N e^-N p <cit.>. Figures <ref> and <ref> show that our simulations — both exact and approximate — reproduce this trend. As a first approximation of quantum sampling, a competing classical computer might select a bitstring x_j uniformly at random, calculate its probability p(x_j) by simulation, and accept x_j with probability p(x_j)N/M≤ 1, repeating the entire process until acceptance. This rejection sampling assumes a bound M ≥max_x p(x) N and requires an average of M probabilities p(x) per accepted bitstring. It can be implemented by Algorithm 1 in Figure <ref>.
We can estimate the expected number ϵ of probabilities p(x) > M/N using the Porter-Thomas distribution as follows



    Exp[ ∑_p(x) ≥ M/N  1 ] =
        ϵ =  N ∫_M^∞ e^-t d t = N e^- M .

Hence for given ϵ we can choose M = ln(N/ϵ). For ϵ = 10^-3 and n = 49, we need M = ⌈ln(10^3 2^49) ⌉ =  41 probabilities per bitstring. 

Frugal rejection sampling. To reduce the number M of probabilities per output bistring to M' without significantly increasing the error, we set M'(ε) so that 
∑_x : p(x) > M' /N^N p(x) ≤ε. The distribution p̃(x) from which Algorithm 1  samples has statistical variational distance from the true distribution <cit.> 



    1/2∑_x | p(x) - p̃  (x)| = 
          Exp[∑_x : p(x) > M' /N p(x)]  ≤ ε .

To get a tight bound on the number of probabilities needed per bitstring,
we use Porter-Thomas statistics [In comparison, boson sampling requires
∼ 100 probabilities per bistring <cit.>, despite the
use of Metropolised independence sampling.]

    Exp[
    ∑_x : p(x) > M' /N p(x)  ] = ∫_M'^∞ t e^-t d t = e^- M' (1+M') .

The statistical distance decreases exponentially with the number of probabilities M' per bitstring sample. M' = 10 gives statistical distance equal to ε = 5 · 10^-4, independent of the number 
of qubits. 
To further decrease ε, we skip Step 4 in Algorithm 1 in Figure <ref> 
and always accept x_j when p(x_j) ≥ M' / N in Step 4 in Algorithm 2. In other words,
instead of ignoring high-probability bitstrings, we sample them with a lower probability. 
Extending earlier calculations gives



    1/2∑_x | p(x) - p̃  (x)| = 2 exp( -  M'  1 - e^- M') .

For M' = 10, this reduces the statistical distance (between p(x) and p̃(x))
to 9 · 10^-5. 

Our frugal rejection-sampling can be applied without assuming Porter-Thomas statistics.
For such uses, one cannot rely on analytical error bounds, but can calculate the sampling error
(for a given M') for each batch of ℓ M' probability values. 
We numerically obtain
ε = ∑_x : p(x) > 10 /N p(x) = (4.4 ± 0.7) · 10^-4 in the case of Figure <ref> (exact simulation) [The errors were estimated by bootstrapping.], and ε = (5.3 ± 0.7) · 10^-4 in the case of
Figure <ref> (approximate simulation).
Since each figure shows a good fit to Porther-Thomas, we can compare these error estimates
to the analytical prediction 5 · 10^-4 based on the Porter-Thomas distribution.
These errors and the runtime of rejection sampling are negligible in the context of
approximate circuit simulation, where ℓ M' probabilities are produced per run for Algorithm 2 
from Figure 1. Then ℓ±√(ℓ M) bitstrings are generated, as seen from McDiarmid's 
inequality.






§ SIMULATION RESULTS



The infrastructure used by our simulation is described in 
Sec. <ref>. 
To put our results in perspective, Appendices <ref> and <ref> compare Rollright (without approximation) to the most recent simulators from Microsoft and IBM on circuits up to 30 qubits. Rollright shows speed advantage of 31× and 16× over Microsoft QDK and IBM QISKit-Terra/QASM, respectively. Running QISKit-Terra on Power9 servers with NVIDIA GPUs through the IBM Q service decreases Rollright's speed advantage to 12.8×, but for a 32-qubit circuit Rollright is 877× faster. Without using GPUs, Rollright simulates 30-qubit circuits on a MacBook Pro at least twice 
as fast as QISKit on GPUs. Rollright's memory advantage over Microsoft and IBM is 886× and 590×. Appendix <ref> shows that Rollright considerably outperforms
the simulator from <cit.> on 64-qubit circuits of depth 22.

For simulations of “easy” quantum circuits (42 and 45 qubits), we use 4 threads per process with 22 parallel processes on one n1-highcpu-96 server with 96 vCPUs. [On Google Cloud, a virtual CPU is implemented as a single hardware hyper-thread.] 
Simulations for up to 45 qubits were performed in exact mode,
 in configurations that match those previously reported, so as to facilitate comparisons (the gate counts
 do not always match exactly, but these differences have only minor impact on results). As implied by
 Table <ref>, we use much more modest resources than supercomputer-based simulations from 2016 and 2017, and achieve lower runtimes than cluster-based simulations. In particular, a 45-qubit simulation of depth 1+25, previously performed on a supercomputer with 0.5PB RAM <cit.>, took $0.24 and 17.4 GiB RAM with our algorithms. This illustrates rapid progress in the understanding and availability of quantum circuit simulation. 



 

Runtime variation is important to track for two reasons.


  * Simulators in some prior work <cit.> appear exponentially sensitive to quantum-supremacy circuit instances, possibly because they exploit easy-to-simulate gate sequences that can be removed (see Section <ref>). This is not the case in our work.
For different circuits with the same depth and numbers of qubits, same placement of xCZs and similar number of gates, the runtime of our simulator is also similar. For a given circuit,
the amount of work performed by individual branches of simulation and their runtimes on a given CPU do not exceed mean values for that CPU by more than 10-20%.

  * Using
 servers, especially in a shared Cloud Computing environment, exposes simulations to hardware-performance variations possibly due to CPU diversity and transient loads at compute nodes (jobs by other users can schedule on the same compute nodes), see Figure <ref>.


Approximate simulation of “hard” quantum circuits.
Simulations for 7×7 and 7×8 qubit arrays are more difficult than simulations we used
 for benchmarking purposes above, and are performed to a greater circuit depth in our work,  1+39+1 and 1+40+1, including
 the initial and the final layers of Hadamard gates <cit.>. Depth of at least 40 (plus Hadamards) was suggested <cit.> as sufficient to be hard for simulations optimized for low depth circuits. It is easier to simulate circuits on oblong arrays than on square arrays, although specific algorithms may be affected in different ways <cit.>.
 For our simulation,
 

  * we choose a cut with 7 xCZ gates every 8 cycles,

  * the 7×7 array does not admit a balanced cut, and
 the blocks formed by the most balanced, smallest cut contain 28 and 21 qubits, whereas the 7×8 array admits an even 28 + 28 partition.
 
 This level of simulation difficulty offers a good opportunity to evaluate our approximate simulation technique.
 Therefore, we use f=10^-2 for the 7×7 array depth  1+39+1 , and f=1/196 for depth 1+40+1 .
Distributed simulations were performed on shared physical hardware in Google Cloud,
using up to 625 n1-highcpu-32 virtual servers for 7 × 7 qubits and 625 n1-highmem-32 servers for 7×8 qubits, with 32 vCPUs each. Key results are reported in Table <ref>,
and details are available in appendices. These results can be used to estimate resources for different fidelity values f because runtime scales linearly with f. The linear scaling does not depend on the number of qubits or circuit depth.

We also performed a simulation of a circuit with 7× 7 qubits and depth 1+48+1 with fidelity 2^-22 using 512 n1-highcpu-32 virtual servers with an estimated cost $52. Using the linear scaling of computational cost with fidelity, this implies a cost of one million dollars for 0.5% fidelity.





 

Scaling simulation to produce multiple bitstring samples. As discussed in Section <ref>, we strive to emulate a quantum computer that samples bitstrings from a circuit-defined
distribution. The fidelity f of the output of a noisy quantum computer can be estimated with cross-entropy benchmarking <cit.>, which requires sampling ∼ 1/f^2 bitstrings. For the fidelities considered in this paper, consistent with the noise estimates for near-term quantum computers, this requires to produce at least one million bitstrings or repeated runs.

To produce bitstring samples, we first select a number of bitstrings uniformly at random, then simulate the circuit to save the bitstring probabilities, and then use methods from Section <ref>
to produce a bitstring per ten probabilities on average. Therefore, to match  a million runs of a quantum computer, a simulator would produce ten million amplitudes with randomly selected indices. For algorithms that produce one amplitude at a time <cit.>, this implies a 10^7× slow-down to produce 10^6 bitstrings. 
In contrast, our simulator produces multiple amplitudes in a batch. Table <ref> shows
the scaling of runtimes when the number of amplitudes saved (with random indices) increases
from 1 to 10^8. Saving 10^6 amplitudes leads to a 17% slowdown because the bottleneck
is in the simulation itself. Saving 10^7 amplitudes leads to a 2.76× slowdown. Producing amplitudes in a batch gives us a 10^7/2.76 = 3.6×10^6 speedup over producing one amplitude at a time.



§ REFINED QUANTUM-SUPREMACY BENCHMARKS



We now review the difficulty of Quantum Supremacy benchmarks for sampling the output distribution of random quantum circuits. Analysis in <cit.> suggests circuits with at least 7 × 7 qubits, as the quantum state of such a system is too large for typical Schrödinger-style simulations. With the same 7 × 7 array,
circuit depth must be at least 1+40+1, given the exponential growth of simulation difficulty with depth for  simulation by variable-elimination <cit.>. Note that the simulation of sampling requires calculating a large number of amplitudes requested at random <cit.>, whereas variable elimination tends to produce one amplitude at a time. We have shown in previous sections how to simulate 7 × 7 circuits with CZ gates of depth 1+40+1 trading computational cost for fidelity. This suggests revising the initially proposed benchmarks <cit.>.

 Quantum Supremacy entails an inherently adversarial protocol that asymmetrically favors quantum computers — a computational problem is being selected that can be solved by quantum yet not classical computers <cit.>. Reliably defeating Quantum Supremacy requires more than a handful of opportunistic simulations as one has to anticipate modifications of problem instances that complicate simulation more than they complicate quantum evolution. It is sometimes easy to confuse Quantum Supremacy with almost the opposite, i.e., showing that classical computers can solve tasks that present-day quantum computers cannot. The later is trivial and can be demonstrated on many common tasks for which classical software and hardware excel, while quantum computers have no algorithmic advantage and are currently much smaller and error-prone. For example, sorting 2^n numbers requires Ω(n2^n) time on both classical and quantum computers <cit.> and can be accomplished just as quickly. Of course, such statements do not preclude Quantum Supremacy. On the other hand, being able to quickly simulate only some circuits of a kind <cit.> would not negate Quantum Supremacy, e.g., the work in <cit.> attempts to simulate a number of similar circuits but then abandons the ones that time out. The authors of <cit.> change the layout of CZ gates to make the simulation approximately 16 times faster at depth 1+22. 
Simulations in <cit.> attain additional depth by omitting the final layer of Hadamard gates, which makes the CZs in preceding cycles unnecessary. Moreover, the work in <cit.> calculates a single amplitude at a time, which is insufficient for sampling.
Finding unexpected simulation shortcuts in existing quantum-supremacy benchmarks <cit.> is also of limited value if these shortcuts can be invalidated.  To this end, the work in <cit.> describes how to simulate eight additional cycles in some cases by exploiting sequences of diagonal gates .  Avoiding such “easy” sequences during benchmark generation does not affect the execution on a quantum computer. Thus, Google benchmarks have been revised to
 

  * replace every T gate appearing after a CZ gate with a non-diagonal gate,

  * explicitly include a cycle of Hadamard gates before measurement,

  * reorder cycles so that “horizontal” two-qubit gates alternate with “vertical” two-qubit gates.

While keeping the same gate family, these minimal changes complicate many approaches to simulation. 
As seen in Tables <ref> and <ref>, the benchmarks retain numerous non-Clifford gates to hamper stabilizer-based simulation <cit.>.
Revised circuit benchmarks are publicly available <cit.> to facilitate fair runtime comparisons among different simulation methods. To ensure the replicability of results, we recommend mentioning specific circuit files, as we do in our work. The performance of our simulator on revised benchmarks
is compared to its performance on older benchmarks (see Table <ref>). The slowdown is moderate and primarily due to the slight increase in the number of gates, especially the Hadamards at
the end of the circuit.



§ HARDER BENCHMARKS THROUGH TWO-QUBIT GATES



In considering more drastic changes to circuit benchmarks, we estimate the impact of upgrading the gate library on simulation difficulty and on specific algorithms used in our work.
 To this end, replacing existing quantum gates within tensor blocks might require additional provisions for simulators that implement speedups for particular gate types, but would result only in a constant-times runtime difference. However, replacing cross-block CZ (xCZ) gates with more complicated gates can increase the branching factor in our multiprocessor simulation.  The 4x4 matrix of a generic two-qubit gate can be decomposed into a sum of four tensor products of one-qubit gates

    (
       [ A_00 A_01; A_10 A_11;      ])
      = _0⊗ A_00 +  _1⊗ A_11 +  |0⟩⟨ 1| ⊗ A_01 +|1⟩⟨ 0|⊗ A_10

where the A terms represent 2x2 blocks of A.[A similar decomposition can be written out with the projections appearing in the second tensor component.]
 For example,

    iSWAP = _0⊗_0 + _1⊗_1+  i |0⟩⟨ 1|⊗ |1⟩⟨ 0| + i |1⟩⟨ 0|⊗ |0⟩⟨ 1| .

Note that non-unitary gates pose no difficulties to common simulation algorithms, and projections
are particularly convenient because they create numerous intermediate zero amplitudes that can be skipped
by simulators that track patterns of nonzeros (cf. Equation (<ref>)). Lemma 1 in <cit.> shows that such decompositions cannot be shorter than the operator Schmidt decomposition (which must have rank 1, 2 or 4 for a unitary operator). The operator Schmidt decomposition for iSWAP

    2 · iSWAP = ⊗ + i ⊗ +
     + i ⊗ + ⊗

proves that the decomposition in Equation (<ref>) is minimal.
If all CZ gates are replaced with iSWAPs, the branching factor in multiprocess simulation increases from two to four, halving the circuit depth that can be handled within a given amount of time. Similar, if not greater, complications hamper algorithms which map the circuit to undirected graphical models and use variable elimination <cit.>. 

Using an arbitrary two-qubit gate can slow down simulation of individual qubit blocks and also decrease the number of intermediate zero amplitudes, but the iSWAP gate has been well-studied and implemented in leading quantum chips <cit.>. Moreover, an arbitrary two-qubit gate may be easier to approximate than an iSWAP by truncating the least-norm terms in its operator Schmidt decomposition. To create a new benchmark suite, while preserving the structure of existing benchmarks, we replace every CZ gate with an iSWAP gate. This benchmark suite is also available at <cit.>.

In addition to considering more general one-qubit and two-qubit gates, the benchmarks can be modified
to include two-qubit gates that couple arbitrary (non-adjacent) qubits. When such gates are applied within tensor partitions, our simulator will treat them just as efficiently as it treats gates in existing benchmarks. However, when such gates cross partitions, they will contribute to the branching factor of simulation. Nevertheless, our work focuses on simulating realistic quantum-circuit architectures, such as those pursued by IBM, Google, Rigetti and others. To this end, existing benchmarks are sufficient.



§ CONCLUSIONS

Universal quantum computers based on several different technologies started scaling beyond a handful of qubits around 2011. Today, detailed reports are published of programmable quantum-circuit computers with 22 qubits, while 50-qubit systems and larger are under testing <cit.>. Two-qubit gate fidelities have improved within the 0.9 - 0.995 range <cit.>, but continue to limit implementable circuit depth. Numerical simulations traditionally focused on straightforward full-wave-function kernels, where progress from 42 qubits in 2010 to 45 qubits in 2018 <cit.> had been bounded by exponential memory scaling. The more recent focus on limited-depth quantum circuits extended supercomputer-based simulations well past 50 qubits <cit.>, facilitated the use of GPUs and distributed cluster without interprocess communication <cit.>, and enabled meaningful simulations on single servers <cit.> for up to 100 qubits.

We show that full-circuit fidelity in quantum-circuit simulation can be traded off for computational resources, and a simulation can match requested circuit fidelity to compete with a given quantum computer. We also show that (i) competitive simulations of quantum-supremacy circuits do not require specialized systems, and (ii) monetary cost can be estimated for cloud-based simulations. Within a given physical-time span, supercomputers might be able to simulate larger quantum circuits than those we simulated, but perhaps at a much greater cost. 

The Schrödinger-Feynman simulation that we developed compares favorably to competitors:


  * Stabilizer-based methods <cit.>
         are exponentially sensitive to the number of non-Clifford gates, whereas diagonal gates are easy
         for our approach. To this end, simulations in <cit.> involve circuits with up to 60 such
         gates, whereas the circuits we simulate include hundreds of them. 

  * Bucket-elimination methods  <cit.> handle diagonal gates well,
        but so far produce one amplitude at a time, whereas we can produce millions. Additionally,
        bucket-elimination methods have not yet been evaluated on revised quantum-supremacy
        benchmarks from Section <ref>, which seem harder for this approach.

Stabilizer-based techniques have the advantage of better handling non-nearest-neighbor circuits, whereas
this generalization would impact bucket elimination methods more than it would impact our approach. However, leading quantum computers are consistent with the nearest-neighbor architecture.

Our work emphasizes rapid progress in simulation algorithms
and its implications for quantum supremacy experiments, as follows.


  * Limited-depth quantum circuits now admit surprisingly efficient simulation on cloud computers.

  * Deep quantum circuits without error correction appear out of reach for both simulation
(due to exponentially-growing resource requirements) and quantum chips (due to the exponential accumulation of gate errors). Quantum error correction is currently out of reach for near-term quantum computers <cit.>.

  * Schrödinger-Feynman hybrid simulation used in our work favors conventional computers over GPUs
and does not require fast interconnect found in supercomputers.
It can be cast into an Internet-scale -style effort, where transient participants all over the world contribute CPU time on their computers to a joint effort.

  * Approximate simulation of quantum circuits is surprisingly effective in practice [Complexity-theoretic arguments against the existence of classical polynomial-time approximation algorithms for this sampling problem <cit.> do not conflict with the algorithm in our work, as our algorithm requires exponential time.], and this raises the bar for quantum supremacy experiments because low quantum circuit fidelity (due to gate errors) facilitates efficiency gains in approximate simulation. To this end, quantum computers incur stochastic errors in qubit initialization, measurement, one- and two-qubit gates, and gate errors tend to increase with the number of qubits. In contrast, our simulation algorithms incur deterministic inaccuracies in only a very small fraction of two-qubit gates, and the complexity of simulation scales linearly with circuit fidelity f.

  * We show how to interactively validate the results 
        simulations with lower-fidelity simulations.

  * Second-generation benchmarks, with iSWAP gates instead of CZ gates, appear
  significantly harder to simulate, yet implementable on a quantum computer. We estimate that changing CZ gates to iSWAP gates reduces circuit depth attainable by leading simulation methods approximately two-fold. In this context, two-qubit gate fidelity remains a key parameter to demonstrate during many-qubit operation.

Pricing quantum-supremacy simulations on cloud-computing services in terms of cost per amplitude should facilitate fair comparisons between different hardware types and can reflect the lower cost of idling general-purpose resources. Comparisons with quantum computers by cost are also possible and can help (i) making a long-term case for commercial quantum computers, (ii) pricing quantum computing resources. Such comparisons must account not only for the resources used by each simulation and each quantum computer run, but also the number of runs required to produce trustworthy results.
When sampling outputs on a quantum computer based on superconducting qubits, reported times are around 45 ns per cycle with CZ gates, 25 ns per cycle with only single-qubit gates <cit.>. Readout can be as fast as 140 ns <cit.>, and qubits can be initialized to a known state in a few hundred ns <cit.>. Thus, comparisons to simulation can be performed once we know circuit fidelity values for 49 or 56 qubits, as well as the cost of quantum-computing resources. As both quantum and classical computing technologies continue to improve and become less costly, quantum supremacy remains a moving target and will require careful comparisons.



§ INFRASTRUCTURE AND PERFORMANCE ESTIMATION


Performance estimation for our simulation is based on the following parameters.
The qubit array is partitioned in two blocks q = q_1 + q_2, attempting to minimize
the larger block and reduce the number of xCZ gates (x). We use a single straight line cut in the qubit array for each simulation.
As explained before, we split x=x_p + x_b. The corresponding split for circuit depth
is d = d_p + d_b. The simulation is configured to save n_a > 0 amplitudes with specified indices.
The total runtime of all simulation processes (however many threads each process may use)
can then be estimated as

    T_tot (f, q_1,q_2,d_p,d_b, x_p, x_b, n_a) = C_1 f 2^x_p (q_1 2^q_1 + q_2 2^q_2) ( d_p + C_2 2^x_b d_b) + C_3 2^x_p + x_b n_a

  where C_1, C_2, C_3 > 0 are implementation-specific constants.

 In each of the d = d_p + d_b cycles, we simulate on the order of q gates, of which q_1
 gates act on 2^q_1 amplitudes and q_2 gates act on 2^q_2 amplitudes. After the first
 d_p cycles, the simulation state is saved and reused in 2^x_b branches over d_b cycles.
 At the end of each branch, n_a requested amplitudes are collected (unlike the work in <cit.>,
 our simulator can save any subset of amplitudes specified by indices). For n_a= 10000, the last term
 is significant for q≤ 36, but can be neglected when n_a ≪ 2^q_1 + 2^q_2. Given a simulation
 cluster that can run p simulation processes per node on N nodes, we distinguish billable time and wallclock time


    T_ bill = ω(p) T_ tot/p         
      T_ clock = T_ bill / N

where  ω(p)≥ 1 reflects the slowdown of an individual process sharing a node with other p-1 such processes,
e.g., due to memory contention. While ω(1)=1, for the more common case 16 ≤ p ≤ 32, ω(p) can be in the 1.5-2 range, depending on
how much memory bandwidth a particular system offers and how much bandwidth is required by each process.
In terms of RAM usage, a single process requires

    M_ proc =  (C_4 (2^q_1 + 2^q_2) + n_a) ·

bytes, where C_4>0 is a small integer. Each compute node requires p M_ proc bytes,
and the entire cluster peaks at p N M_ proc bytes. Circuit depth does not affect RAM requirements.

For the simulations reported below, the billable time and wallclock time are measured directly.

 

Infrastructure used in our simulations has been as a follows. Key algorithms are implemented  in C++11 in a package called Rollright (University of Michigan) and compiled with g++ 7.2. Our kernels use AVX-2 instructions, but not AVX-512 instructions used in <cit.>.  We also use python scripts for several tasks. Empirical performance is evaluated on quantum supremacy circuits <cit.>, and results are shown in Tables <ref> and <ref>, with additional details given in the appendices. Unlike prior published efforts, we use no GPUs, no supercomputers, and no high-performance node-to-node interconnect.  Software development was performed on a MacBook Pro 2017 with 16 GiB RAM, where our baseline Schrödinger simulation
completes a 30-qubit circuit with depth 1+26+1 in 71 s using a little over 8 GiB of RAM. For circuits with more than 35 qubits, we use the Google Cloud, a common cloud-computing service where users rent networked servers of various types. We used several virtual server types (specified for each simulation below) with Ubuntu Linux, available to the public for fair comparisons. Most of the RAM available on those servers was not used in our simulations. For each server type, we report hourly rental cost in preemptible mode as of June 2018. Cost comparisons can be performed with other server types, adjusted for other price points, and extended to include quantum computers, so as to establish or refute quantum supremacy.


 

Practical validation of results is a critically-important step in quantum-supremacy simulations because
subtle bugs or misguided optimizations can significantly improve runtime, while producing incorrect results.
For circuits of up to 32 qubits we validated our main Schrödinger simulator against several independently-developed simulators based on different algorithms, to ensure that the output amplitudes match with good accuracy (here we used the publicly available QuIDDPro <cit.> simulator and the simulators from
IBM QISKit-Terra, as well as several in-house simulators). We then simulated circuits using our Schödinger-Feynman simulator and validated results against a Schrödinger-only simulation, and  against a variable elimination simulation <cit.> on low-depth circuits with up to 64 qubits. Computing state norms offers a good sanity-check for both exact and approximate simulation. We also validated approximate simulation against exact simulation for low-depth circuits by producing amplitudes for the same subset of indices and estimating state-overlap (fidelity) based on these amplitudes.[We also calculated cross-entropy difference <cit.>. This gave similar fidelity estimates, but with a greater variance of
estimation errors.] Another sanity check is to plot the distribution of probability values for a given set of
amplitudes. For exact simulation, the results should closely match the Porter-Thomas distribution, see Figure <ref>. While this is not required for approximate simulation, our method also produces a Porter-Thomas distribution in this case, see Figure <ref>. 






To validate the results of  approximate simulation in this work, we performed an exact simulation of  circuit inst_6x6_41_0.txt (from <cit.>) with 6 × 6 qubits  and depth 1+40+1, and saved 10^5 amplitudes. We calculated the same amplitudes using 40 different random fractions for each f = 1/16, 1/32, 1/64, and 1/128 of paths. Figure <ref> shows that empirical fidelity matches the fraction of paths used. The runtimes scale linearly with f. We also performed an exact 7×7 qubit depth 17 simulation and saved 1M random amplitudes. Then we performed approximate simulation with requested fidelities f=10^-4, f=10^-2 and f=1/64 and saved the same amplitudes. We estimated the fidelity of approximate state by computing the inner product, and the fidelities matched the requested values to within 2%. 
Additionally, Appendix <ref> describes our novel technique for validating results of higher-fidelity simulation using randomized lower-fidelity simulation.

naturemag
10
url<#>1urlprefixURL 


aaronson2011computational
authorAaronson, S. & authorArkhipov, A.
titleThe computational complexity of linear optics.
In booktitleSTOC, pages333–342
  (organizationACM, year2011).

bremner2015average
authorBremner, M. J., authorMontanaro, A. &
  authorShepherd, D. J.
titleAverage-case complexity versus approximate simulation
  of commuting quantum computations.
journalPhys. Rev. Lett.
  volume117, pages080501
  (year2016).

Preskill
authorPreskill, J.
titleQuantum computing and the entanglement frontier
  (year2012).
note25th Solvay Conf.

bremner16
authorBremner, M. J., authorMontanaro, A. &
  authorShepherd, D. J.
titleAchieving quantum supremacy with sparse and noisy
  commuting quantum computations.
journalQuantum volume1,
  pages8 (year2017).

GoogleSuprem
authorBoixo, S. et al.
titleCharacterizing quantum supremacy in near-term
  devices.
journalNat. Phys. volume14,
  pages595 (year2018).
notearXiv:1608.00263.

NatureSupremacy
authorHarrow, A. W. & authorMontanaro, A.
titleQuantum computational supremacy.
journalNature volume549,
  pages203 (year2017).

aaronson2016complexity
authorAaronson, S. & authorChen, L.
titleComplexity-theoretic foundations of quantum supremacy
  experiments.
journalarXiv:1612.05903  (year2016).

neill_blueprint_2017
authorNeill, C. et al.
titleA blueprint for demonstrating quantum supremacy with
  superconducting qubits.
journalarXiv:1709.06678  (year2017).

bouland_quantum_2018
authorBouland, A., authorFefferman, B.,
  authorNirkhe, C. & authorVazirani, U.
titleQuantum Supremacy and the Complexity of Random
  Circuit Sampling.
journalarXiv:1803.04402  (year2018).
noteArXiv:1803.04402.

de_raedt_massively_2007
authorDe Raedt, K. et al.
titleMassively parallel quantum computer simulator.
journalComputer Physics Communications
  volume176, pages121–136
  (year2007).

PB
authorHäner, T. & authorSteiger, D. S.
title0.5 Petabyte Simulation of a 45-Qubit Quantum
  Circuit.
journalarXiv:1704.01127  (year2017).

IBM
authorPednault, E. et al.
titleBreaking the 49-qubit barrier in the simulation of
  quantum circuits.
journalarXiv:1710.05867  (year2017).

GoogleBucket
authorBoixo, S., authorIsakov, S. V.,
  authorSmelyanskiy, V. N. & authorNeven, H.
titleSimulation of low-depth quantum circuits as complex
  undirected graphical models.
journalarXiv:1712.05384  (year2017).

64q
authorChen, Z.-Y. et al.
title64-qubit quantum circuit simulation.
journalarXiv:1802.06952  (year2018).

Taihu
authorLi, R., authorWu, B., authorYing, M.,
  authorSun, X. & authorYang, G.
titleQuantum supremacy circuit simulation on Sunway
  TaihuLight.
journalarXiv:804.04797  (year2018).

chen_classical_2018
authorChen, J., authorZhang, F.,
  authorHuang, C., authorNewman, M. &
  authorShi, Y.
titleClassical Simulation of Intermediate-Size
  Quantum Circuits.
journalarXiv:1805.01450  (year2018).

de_raedt_massively_2018
authorDe Raedt, H. et al.
titleMassively parallel quantum computer simulator, eleven
  years later.
journalarXiv:1805.04708  (year2018).

ibm_2017
titleIBM Announces Advances to IBM Q Systems &
  Ecosystem, IBM Press Release (year2017).

intel_18
titleQuantum Computing.
<https://newsroom.intel.com/press-kits/quantum-computing/>.

google_2018
titleA Preview of Bristlecone, Google’s New Quantum
  Processor, Google AI Blog (year2018).

neville_classical_2017
authorNeville, A. et al.
titleClassical boson sampling algorithms with superior
  performance to near-term experiments.
journalNat. Phys. volume13,
  pages1153–1157 (year2017).

bernstein1997quantum
authorBernstein, E. & authorVazirani, U.
titleQuantum complexity theory.
journalSIAM Journal on Computing
  volume26, pages1411–1473
  (year1997).

markov_simulating_2008
authorMarkov, I. L. & authorShi, Y.
titleSimulating quantum computation by contracting tensor
  networks.
journalSICOMP volume38,
  pages963–981 (year2008).

new_benchmarks
authorBoixo, S. & authorNeill, C.
titleThe question of quantum supremacy.
journalGoogle AI Blog  (year2018).
noteAvailable on GitHub at
  <https://github.com/sboixo/GRCS>.

barends_superconducting_2014
authorBarends, R. et al.
titleSuperconducting quantum circuits at the surface code
  threshold for fault tolerance.
journalNature volume508,
  pages500–503 (year2014).

barends_digitized_2015
authorBarends, R. et al.
titleDigitized adiabatic quantum computing with a
  superconducting circuit.
journalNature volume534,
  pages222–226 (year2016).

kelly_state_2015
authorKelly, J. et al.
titleState preservation by repetitive error detection in a
  superconducting quantum circuit.
journalNature volume519,
  pages66–69 (year2015).

bravyi_improved_2016
authorBravyi, S. & authorGosset, D.
titleImproved classical simulation of quantum circuits
  dominated by Clifford gates.
journalPhys. Rev. Lett.
  volume116, pages250501
  (year2016).
  
bravyi_2018
authorBravyi, S. et al.
titleSimulation of quantum circuits by low-rank stabilizer decompositions.
journalarXiv:1808.00128 (year2018).


Viamontes
authorViamontes, G. F., authorMarkov, I. L. &
  authorHayes, J. P.
titleQuantum circuit simulation
  (publisherSpringer Science & Business Media,
  year2009).

hoyer2002quantum
authorHøyer, P., authorNeerbek, J. &
  authorShi, Y.
titleQuantum complexities of ordered searching, sorting,
  and element distinctness.
journalAlgorithmica volume34,
  pages429–448 (year2002).

gottesman1998heisenberg
authorGottesman, D.
titleThe heisenberg representation of quantum computers.
journalarXiv: quant-ph/9807006
  (year1998).

nielsen2003quantum
authorNielsen, M. A. et al.
titleQuantum dynamics as a physical resource.
journalPhys. Rev. A volume67,
  pages052301 (year2003).

neeley2010generation
authorNeeley, M. et al.
titleGeneration of three-qubit entangled states using
  superconducting phase qubits.
journalNature volume467,
  pages570 (year2010).

mckay2016universal
authorMcKay, D. C. et al.
titleUniversal gate for fixed-frequency qubits via a
  tunable bus.
journalPhys. Rev. Applied
  volume6, pages064007 (year2016).

jeffrey2014fast
authorJeffrey, E. et al.
titleFast accurate state measurement with superconducting
  qubits.
journalPhys. Rev. Lett.
  volume112, pages190504
  (year2014).

reed2010fast
authorReed, M. et al.
titleFast reset and suppressing spontaneous emission of a
  superconducting qubit.
journalAppl. Phys. Lett.
  volume96, pages203110
  (year2010).

magnard2018fast
authorMagnard, P. et al.
titleFast and unconditional all-microwave reset of a
  superconducting qubit.
journalarXiv:1801.07689  (year2018).





 §.§ Acknowledgments

We acknowledge support from Kevin Kissell on running computations in Google Cloud and reviewing the manuscript,
as well as helpful comments from Edward Farhi, Héctor Garcia, Riling Li, Dmitri Maslov, Hartmut Neven,
David Wecker, Robert Wille and Alwin Zulehner.




§ COMPARISONS TO MICROSOFT QDK


To put the capabilities of our simulator in perspective, we compared it to the Microsoft Quantum Development Kit (QDK) v0.2.1806.3001 (June 30 2018).[Downloaded from https://www.microsoft.com/en-us/quantum/development-kit in July 2018.]
This software includes a back-end circuit simulator and front-end support for the Q# language, integrated with Microsoft Visual Studio and available on Windows, MacOS and Linux. Originally released to the public in the end of 2017, the package underwent significant improvements in February and June 2018, particularly in simulation speed. 

To compare the QDK simulator to our simulator Rollright, we used revised (v2) circuit benchmarks described in Section <ref> and converted circuit descriptions into the Q# language. In particular, circuit depth 1+26+1 indicates layers of Hadamard gates at the beginning and at the end. Comparisons were performed on a MacBook Pro 2017 with 16 GiB RAM and Intel Core i7-7700HQ CPU (2.80GHz). Results are reported in Table <ref>. To exclude code segments from memory comparisons, we first measured max resident memory for each simulator on the 16-qubit benchmark and then used those measurements as baselines. For 24- and 25-qubit benchmarks we report the difference between max resident memory and the respective baseline. Given that the baseline includes memory used by 16-qubit data structures, these estimates are a little lower than the actual data structures, but the inaccuracy is close to 1 MiB and removed by rounding, as can be seen from Rollright data — 128 MiB and 256 MiB are the actual sizes of amplitude vectors in Rollright. We note that Q# uses double-precision floating-point numbers to represent amplitudes, whereas Rollright uses single-precision numbers (while keeping track of accuracy through higher-accuracy norm computations).
Hence, the memory-usage ratio should be at least 2.0 on any platform. The runtime ratios are more remarkable and are explained by two factors: (i) a certain slowdown due to the use of the Microsoft .NET framework, (ii) algorithmic differences. In particular, algorithmic differences explain why the runtime ratio grows with the number of qubits.

We note that for circuits with 30 qubits and fewer, the best runtimes on a laptop are attained with the basic Schrödinger capability of our simulator. However, the Microsoft QDK simulator required more than the 16 GiB memory available (Rollright used a little over 8 GiB). We therefore installed Microsoft QDK on an 18-core Linux server with sufficient memory and observed that it used all 72 available threads. 
We first configured Rollright in the Schrödinger mode.  As shown in Table <ref>, Rollright is 19× faster than the QDK simulator with a 3× memory advantage. The ratios remain similar for circuits of larger depth, and the memory ratio does not depend on how many threads Rollright uses. We then evaluated Rollright in the Schrödinger-Feynman mode using the even (30=15+15) cut, with 18 parallel processes and four threads per process.
Rollright was 31× faster and used 886× less memory. The memory ratio depends linearly on how many simultaneous processes Rollright uses. The runtime ratio grows exponentially with circuit depth (specifically, with the number of xCZ gates). Detailed logs of Rollright simulations
are provided in Appendix <ref>.

While Rollright significantly outperforms the Microsoft SDK simulator in our comparisons, we consider only quantum-supremacy benchmarks, whereas Microsoft developers may have had different priorities and continue improving their software. Our algorithmic contributions can be useful when simulating many families of quantum circuits, but the technical effort necessary to demonstrate such applications is beyond the scope of this paper. Moreover, Microsoft SDK offers value in areas unrelated to circuit simulation.





§ COMPARISONS TO IBM QISKIT AND IBM Q


IBM QISKit (https://qiskit.org/) is an open-source quantum computing framework, which includes several simulators in the software toolchain. We used QISKit-Terra v0.5.7 (July 19 2018) and followed the comparison protocol in Appendix <ref>. We found the QASM simulator to be faster than the state-vector simulator on quantum-supremacy circuits. Comparisons of our simulator Rollright to QASM 
are reported in Table <ref> and follow the trends in comparisons to Microsoft QDK 
in Table <ref>. On a depth 1+26+1 30-qubit circuit, Rollright runs 16× faster and uses
590× less memory than QASM. On a depth 1+26+1 32-qubit circuit, Rollright runs 430.6×
faster and uses 1331× less memory than QASM, even though QASM exhibits a high degree 
of parallelism with 6010% CPU utilization rate. The 32-qubit circuit uses the oblong 4x8 
qubit array, and the Schrödinger-Feynman mode of our simulator is able to exploit this shape.
Therefore, we consider results for both the 5x6 and the 4x8 qubit arrays important.
Detailed logs of Rollright simulations are provided in Appendix <ref>.

QISKit-Terra offers the option to run simulations on remote servers
with Power9 CPUs and NVIDIA GPUs through the IBM Q service. The last four rows of Table <ref> compare such simulations with single-server CPU-based Rollright simulations in two modes, with considerable wins for Rollright. Enabling GPU support in the IBM simulator does not impact memory usage (we do not measure memory usage inside the GPU), but improves runtime by almost 30% on 30-qubit circuits.
However, the GPU-based version is twice as slow on 32-qubit circuits. This may be due to the limitations
of GPUs in terms of on-board memory and memory bandwidth, in which case future generations of GPUs
may be helpful. In any case, Rollright exhibits a 430.6× speedup and 1331× memory advantage for 32-qubit circuits, likely beyond upcoming GPU improvements.

We also compared the performance of QISKit on IBM Q to that of Rollright on a MacBook Pro. Even in such an unfair comparison, Rollright is twice as fast, while using 30-70% less memory.
Implementing our algorithms on GPUs is likely to give a significant additional speedup. 
On the other hand, modern GPU resources can be costly. Whether they improve
the total cost of simulation remains an open question.

 




§ VALIDATION OF SIMULATION RESULTS
 


For Google quantum-supremacy circuits, simulation results can be checked using approximate simulation. Consider the correct final state |ψ⟩ of simulation and a norm-1 approximate state |ψ_a⟩ with |⟨ψ|ψ_a⟩|^2 =f. As shown above,
multiple |ψ_a⟩ states can be produced f times faster (each) than the |ψ⟩ state.
Such approximate states help validate claims of producing |ψ⟩. In the following two-party validation protocol, we assume a quantum circuit that the Verifier can simulate to fidelity f<1
in the sense that the Verifier can approximate the final simulation state with fidelity f and do so with
a sufficient variety of states at random. For each state, the Verifier and the Claimant save k amplitudes
so as to estimate fidelities between their states to accuracy δ≤ f.

 
  * The Verifier gives the Claimant k pseudo-randomly generated amplitude indices
          (we used k=1M and specified the indices by a PRNG seed).
 
  * The Claimant simulates the circuit and saves k amplitudes with given indices.
 
  * The Verifier picks a fidelity value 0 ≪ f_1 < 1 and simulates the circuit with fidelity f_1
         producing one of many possible | ψ_a ⟩ states and saving k amplitudes with the specified indices.
 
  * The Claimant shares their amplitudes with the Verifier.
 
  * The Verifier estimates fidelity f_e between the Verifier's and the Claimant's states.
 
  * If |f_e - f_1| > δ, the test fails. Else, the test may be repeated to increase
         the certainty of passing.

A handful of sufficiently unpredictable approximate states
suffice to validate |ψ⟩ because for a generic approximate state with fidelity f_1 ≪ 1,
the expected state fidelity w.r.t |ψ_a⟩ is f_1 f ≪ f. Moreover, the Verifier can estimate the fidelity
of the Claimant's state. The Verifier can select circuits from some family and repeat the test for different circuits.

 


§ 30- AND 32-QUBIT SIMULATIONS



   Below we include detailed simulation logs in support of comparisons to Microsoft and IBM simulators
   in Appendices <ref> and <ref>.



 §.§ 30-qubit depth 1+26+1 simulation on a laptop

This Schrödinger-mode simulation uses the circuit  publicly available from <cit.>.




 §.§ 30-qubit depth 1+26+1 simulation on a server

This Schrödinger-mode simulation uses the circuit  publicly available from <cit.>.




 §.§ 30-qubit depth 1+26+1 cloud simulation

This Schrödinger-Feynman simulation uses the circuit  publicly available from <cit.>.




 §.§ 32-qubit depth 1+26+1 simulation on a server

This Schrödinger-mode simulation uses the circuit  publicly available from <cit.>.





 §.§ 32-qubit depth 1+26+1 cloud simulation

This Schrödinger-Feynman simulation uses the circuit  publicly available from <cit.>.




§  “EASY” 42-QUBIT DEPTH 1+25+1 SIMULATION

This Schrödinger-Feynman simulation uses the circuit  publicly available from <cit.>, Such revised benchmarks  
are harder to simulate than older quantum-supremacy benchmarks
(see Section <ref>).





§  64-QUBIT DEPTH 1+22 SIMULATION


Recent work related to ours <cit.> used GPUs and a distributed cluster with 128 nodes to simulate quantum circuits. The layout of CZ two-qubit gates used in that work was changed relative to the prior work to lower the cost of classical simulation by approximately 16 times, which prevents a direct comparison of runtimes in general: we focus instead on the circuits from Ref. <cit.>, which are significantly harder to simulate. Nevertheless, we performed a simulation of a circuit with 8 × 8 qubits and depth 1+22 with the choice for the layout of CZs from Ref. <cit.>.  Our simulation took 18.6 hours on 64 n1-highmem-32 Google Cloud virtual machines, a 1.7× runtime improvement using half as many nodes and running in a commercial cloud environment.[The nodes in the simulation of Ref. <cit.> have each a Xeon E5-2690V4 processor with 14 cores. We use n1-highmem-32 virtual machines in Google Cloud, with 32 vCPUs or hyper-threads, equivalent to 16 cores 2.0 GHz Intel (Skylake) platform.] 






§ 49-QUBIT LARGE-DEPTH SIMULATIONS

Here we report simulations of circuits on 7× 7 qubits of depth 1+39+1,
1+40+1 and 1+48+1.



 §.§  49-qubit depth 1+39+1 simulation with fidelity 1%

This simulation uses the circuit  publicly available from <cit.>.






 §.§  49-qubit depth 1+40+1 simulation with fidelity 0.51%

This simulation uses the circuit  publicly available from <cit.>.





 §.§ 49-qubit depth 1+48+1 simulation with fidelity 2^-22

This simulation uses the circuit  publicly available from <cit.>.




§  56-QUBIT DEPTH 1+40+1 SIMULATION WITH FIDELITY 0.51%

This simulation uses the circuit  publicly available from <cit.>.



