Evidence for the utility of quantum computing before fault tolerance or maybe not!

A useful application for 127-qubit quantum processors with error mitigation. A work by IBM and UC Berkeley shows the path toward useful quantum computing.

Monit Sharma
24 min readJun 29, 2023
The journal Nature featured the article “Evidence for the use of quantum computing before fault tolerance” on the cover of their June 15, 2023 issue.

This work gained a lot of hype, since it’s a realistic scenario using currently available IBM Quantum processors to explore meaningful computations and realistic applications before the era of fault tolerance. Although some papers published after this showed a better method for the classical simulation part.

In this article, I’ll go through the paper

See the article published in Nature: https://www.nature.com/articles/s41586-023-06096-3

If you’re looking for a good explanation of the work, I suggest reading the official IBM blog on it, since IBM's team is best at explaining what they do in quantum.

Quantum computing hopes to offer substantial speed-ups over their classical counterparts for certain problems, like the Grover search, which offers a quadratic speedup in search, the Shor’s algorithm that offers an exponential speedup, Amplitude estimation, which offers a quadratic speedup over classical Monte Carlo (which uses the same Quantum Phase Estimation as Shor’s algorithm ) and HHL algorithm. However, there’s a big obstacle in realizing or unlocking this potential of quantum computing, and that obstacle is NOISE. Current quantum devices are very noisy, and the widely accepted solution to this challenge is the use of fault-tolerant quantum circuits. But we are not there yet.

In the above paper, the team from UC Berkeley and IBM Quantum ran experiments on a noisy 127-qubit processor and demonstrated the measurement of accurate expectation values for the circuit volume, which is at scales beyond the brute-force method of classical computation. Hence, providing evidence for the utility of quantum computing in the pre-fault tolerance era or in the current noisy era.

It’s important to note that this isn’t a claim that today’s quantum computers exceed the abilities of classical computers — other classical methods and specialized computers already returned correct answers for the calculation they were testing. Like this one:

Link to the paper: https://arxiv.org/abs/2306.14887

We will be learning about both papers.

We know that advanced quantum algorithms like factoring or phase estimation (Shor’s factoring) will require quantum error correction, but even for the smaller circuits the current hardware is not reliable to provide an advantage over the classical counterparts. Current fidelity bounds are very far from the ideal one, like, for example, a quantum circuit with 100 qubits, and a 100 gate depth executed with 0.1% gate error yields a state fidelity of less than 0.0005.

Fidelity is a measure of how close the final quantum state of the real-life qubits is to the ideal case. If the fidelity of the logic gates is too low, calculations will fail because errors will accumulate faster than they can be corrected. Qualitatively, fidelity is the measure of the distance between two quantum states. Fidelity equal to 1, means that two states are equal. In the case of a density matrix, fidelity represents the overlap with a reference pure state. Fidelity equal to 1, means that the square of the density matrix is equal to the density matrix itself and equivalent to the pure reference state.

Using the techniques of quantum error correction can produce accurate results by using classical post-processing. As demonstrated in this article. Quantum Error Correction with Dynamic Circuits.

What is Quantum Advantage?

Quantum Advantage is a significant improvement in quantum algorithm runtime for practical cases over the best classical algorithm, or a stage where the quantum computer can surpass the performance of a traditional computer. It can be demonstrated in two steps:

  1. Demonstrating the ability of existing devices to perform accurate computations at a scale that lies beyond brute-force classical simulations.
  2. By finding problems with associated quantum circuits that derive an advantage from these devices.

This paper on ‘Utility of Quantum Computing in pre-fault tolerant era’ focuses on the first step. They used a quantum device with 127 qubits to run a quantum circuit with up to 60 layers of two-qubit gates, which is a total of 2,880 CNOT gates. These circuits are something not feasible with brute-force classical methods, since even with a supercomputer of 2,801,664 GB of RAM, you can best simulate 47 qubits. (want to know how I got this number, see my previous article on Quantum Volume).

The circuit used in this work is the Trotterized time evolution of a 2D transverse-field Ising model. This is how the circuit looks.

Noise characterization and scaling for 127-qubit Trotterized time-evolution circuit: a. Each trotter step of the Ising simulation includes single-qubit X and two-qubit ZZ rotations, Random Pauli gates are inserted to twirl and controllably scale the noise of each CNOT layer. The dagger indicates conjugation by the ideal layer. b. Three depth-1 layers of CNOT gates suffice to realize interactions between all neighbour pairs on ibm_kyiv. (Image taken from the paper mentioned above.)

What is the Trotterized Time Evolution of a 2D Transverse-Field Ising Model?

In the field of quantum computing, one of the important tasks is to simulate the behavior of quantum systems. The Trotterization technique provides a way to approximate the time evolution of a quantum system by breaking it down into smaller steps.

Transverse-Field Ising Model

Let’s start with a simpler term, a Quantum Spin Model. If a single body quantum problem is described by a Hilbert space H of dimension dimH=d then N distinguishable quantum particles are described by the tensor product of N Hilbert spaces.

So, a quantum spin -½ particle has a Hilbert space H of dimension 2, but N spin -½ particles have a Hilbert space of H^N of dimension 2^N. This exponential scaling of the Hilbert space dimension is a big challenge. For N=30 spins is already a size of 2³⁰ = 10⁹. A single complex vector needs 16GB of memory and may just barely fit into the memory of your personal computer.

The simplest quantum spin model is of the quantum transverse field Ising model (TFIM), which adds a magnetic field in the x direction to a lattice of spin-½ particles coupled by an Ising interaction:

Hamiltonian for a transverse field Ising model

Here the symbol <i,j> means the sum of all bonds in the lattice. In the absence of the second term, the first term is nothing but a classical Ising model. The second term does not exist in the classical Ising model, where the spins point only in the z-direction.

The 2D transverse-field Ising model describes a lattice of spins with interactions between neighboring spins and an external magnetic field applied perpendicular to the spin orientation. J represents the strength of the spin-spin interaction, σz_i is the Pauli-Z operator acting on spin i, σx_i is the Pauli-X operator acting on spin i, Γ is the strength of the transverse magnetic field, and the sum ∑<i,j> runs over all pairs of neighboring spins.

Time Evolution

The time evolution of a quantum system is governed by the time-dependent Schrödinger equation:

Time-Dependent Schrodinger Equation

where H is the Hamiltonian of the system, ψ represents the quantum state, t is time, and ħ is the reduced Planck constant. To simulate the time evolution of the system, we need to calculate the time evolution operator, U(t), which is given by:

How time evolution operator acts on a state.

Trotteriztion

Calculating the exact time evolution operator for a many-body quantum system is generally intractable due to its exponential complexity. Trotterization provides a technique to approximate this operator.

The Trotter-Suzuki decomposition is a popular technique to approximate the time evolution operator by breaking it down into a product of simpler operators. For the 2D transverse-field Ising model, we can split the Hamiltonian into two parts:

where

and

Using the Trotter-Suzuki decomposition, we can approximate the time evolution operator as:

The trotterization of time evolution.

where N is the number of time steps and

is the time step size.

Spin dynamics from an initial state can be simulated by means of the first-order Trotter decomposition of the time-evolution operator, which yields:

The RZZ and RX are the ZZ and X rotation gates respectively. We take this trotterization to be ideal, hence no error from here. For the sake of the experiment, we limit ourselves to the case of

This one assumption makes our life easier since it reduces the ZZ rotation gate to needing only one CNOT (up to a global phase).

Click here to learn more about RZZ gate.

So each trotter step adds a layer of single-qubit rotation RX, followed by commuting layers of parallelized two-qubit rotations RZZ.

This study was done on ibm_kyiv , the 127 fixed frequency transmon qubits, with median T_1 and T_2 times 288 μs and 127 μs, respectively.

What are these T1 and T2 times?

In the field of quantum computing, T1 and T2 times refer to two important parameters related to the coherence and decoherence of quantum states within a quantum system.

T1 time characterizes the decay of the quantum state’s energy or population due to relaxation processes in a quantum system. It measures the timescale over which the system loses information and returns to its equilibrium state. T1 time is associated with the decay of the excited states in the system to lower energy states. In practical terms, T1 time represents the timescale over which the quantum state undergoes spontaneous transitions and loses coherence. It determines the lifetime of quantum information stored in a qubit or quantum register. Longer T1 times are desirable in quantum systems as they allow for longer computation times and reduce the impact of errors.

T2 time refers to the timescale over which the quantum state loses its phase coherence or information due to interactions with the environment. It characterizes the decay of quantum superposition and is influenced by various factors such as noise, magnetic field fluctuations, and interactions with surrounding particles. T2 time primarily represents the decay of off-diagonal elements of the density matrix, which correspond to the coherence between quantum states. During T2 time, quantum states become susceptible to de-phasing, causing errors in quantum computations. Longer T2 times are desirable as they allow for more precise manipulation and preservation of quantum superposition.

It’s important to note that T2 time is generally shorter than T1 time, as various sources of noise and interactions impact coherence more rapidly than energy relaxation.

As each qubit has at most three neighbors, all ZZ interactions can be performed in three layers of parallelized CNOT gates.

b) Figure shows the architecture of 127 qubits (numbered) of ibm_kyiv. Since the ZZ gate is a two-qubit gate, so three layers of depth-1 CNOT gates are sufficient to realize the interactions between all neighbour pairs. As we can see from the figure, qubit 0 to qubit 1 has the third layer of ZZ, qubit 1 to qubit 2 has the first layer of ZZ and qubit 2 to qubit 3 has the second layer of ZZ gates, and so on. The red shows the first ZZ layer, the blue shows the second ZZ layer and the green shows the third ZZ layer.

These kinds of hardware improvements enable us to implement problems on an even larger scale and get a successful execution with techniques of error mitigation.

Techniques of Error Mitigation

Noise in existing quantum processors only enables an approximation to ideal quantum computation. However, these approximations can be
vastly improved by error mitigation, for the computation of expectation value. However, the practical scaling of these methods to larger system sizes remains unknown. The efficacy of the error mitigation is greatly enhanced by additional error suppression techniques and native gate decomposition that reduce the circuit time.

Link to the paper:

Zero-Noise Extrapolation

In this, expectation values are measured for circuits run at varying noise levels to extrapolate to their zero-noise limit. The noise levels in the circuit are typically varied either by stretching the control pulses in time or by insertion of noisy identity-equivalent operations. It is widely used for gate-error mitigation on near-term quantum computers because it can be implemented in software and does not require knowledge of the quantum computer noise parameters.

Read more about Zero Noise extrapolation and how to implement it in code in my previous blog from the same series.

Zero-noise extrapolation for intermediate-scale short-depth quantum circuits: An example of a 26-qubit short-depth quantum circuit composed of layers of single-qubit RX and two-qubit RZZ gates that we employ to study the performance of error mitigation. The noise is amplified for zero-noise extrapolation by appropriately scaling the duration and amplitude of the microwave pulses that compose the gates. A number of additional techniques are employed to enhance the performance of zero-noise extrapolation, that are not incorporated in the bare circuit. The insertion of the dynamic coupling is to suppress dephasing during qubit idle times. Pauli Twirling averages out the off-diagonal coherent error in the Pauli basis and is implemented by sandwiching the two-qubit gates with an additional single-qubit gate drawn from a set of Pauli Operations. Taken from the paper : [2108.09197] Scalable error mitigation for noisy quantum circuits produces competitive expectation values (arxiv.org)

Probabilistic Error Cancellation

Probabilistic error cancellation (PEC) is an error mitigation technique in which ideal operations are represented as linear combinations of noisy operations. In PEC, unbiased estimates of expectation values are obtained by Monte Carlo averaging over different noisy circuits.

Probabilistic error cancellation (PEC) is a noise-aware error mitigation technique which is based on two main ideas:

  • The first idea is to express ideal gates as linear combinations of implementable noisy gates. These linear combinations are called quasi-probability representations
  • The second idea is to probabilistically sample from the previous quasi-probability representations to approximate quantum expectation values via a Monte Carlo average.

A representative noise model is learned and effectively inverted by sampling from a distribution of noisy circuits related to the learned model. Yet, for the current error rate, the sampling overhead is an obstacle, hence the method of Zero Noise extrapolation is used.

You can read more about Zero Noise Extrapolation in my previous blog here.

Link to the paper : [2201.09866] Probabilistic error cancellation with sparse Pauli-Lindblad models on noisy quantum processors (arxiv.org)
Here Λ is the noise affecting our circuits, so we learn the noise and do something equivalent to applying an inverse noise gate, which effectively cancels the noise. The only issue is that this inverse is not a Trace Preserving map and has negative eigenvalues, so basically nonphysical. The above paper suggests a method to achieve this.

Why is the inverse Noise operation non-physical?

Let’s think of an example:

Suppose a quantum circuit, with an Identity gate and the error to be a bit flip error, that means there’s a certain probability p, by which the state of the qubit will be flipped, and a probability 1-p by which it will stay the same. As shown in the figure:

Image Source : Zlatkov Minev’s Lecture on PEC in QIskit Seminar Series.

Visualizing this on the Bloch sphere:

Image Source : Zlatkov Minev’s Lecture on PEC in QIskit Seminar Series.
Bloch Sphere visualization for the above circuit.

Without the noise, the Bloch sphere is intact, but after the noise, the Bloch sphere contracts, and all of the states shrink to the middle of the sphere. So what should be the inverse?

The inverse should be a dilation, an information gain which has therefore negative eigenvalues. It’s a quantum operation that is not allowed.

Inverse of noise map is not physical.

The gist of this thing is that:

If we know the noise, we can cancel the noise with more noise with the tailored noise, using random quantum circuits.

But this thing is difficult to implement. For starters we have to completely know our noise, what it means by that is for a two qubit noise model the matrix for noise is 4^n * 4^n size, hence requires 255 parameters to be learned in general.

As the number increases, the parameters increase, for a 10 qubits circuit, there are 10^12 parameters and for a 50 qubit circuit, we have 10^60 parameters. Storing 10^60 parameters on a classical computer requires 10^50 GB of memory, and that’s not going to happen. Then the question arises:

Is it even possible to learn the noise?

Yes!

Step 1: Simplify the noise (Twirl)

Pauli twirling reduces the 4^n * 4^n noise matrix to a diagonal one with 4^n entries in Pauli basis.

Pauli Twirling (read the paper for more details on it)

Step 2: Amplify and Learn the noise per layer

As the name suggests, do repetitive implementations of the noise which will amplify the noise. But a simple repetition would not be possible since our circuit is made of many entangling layers.

Amplification and Layer wise learning of Noise

Step 3: Twirl readout-error mitigation

We would like to minimize the effect of readout errors, hence we use the twirling-based readout-error mitigation approach.

Sparse Pauli-Lindblad Model

As we discussed above, our noise model was a matrix of 4^n * 4^n size, which by Pauli Twirling we made to a diagonal with 4^n parameters. Our noise model is then represented as:

Noise Map

There are still 4^n Pauli terms here, but there are very few parameters that we actually need to learn. For example, if you have a noise in qubit 1 and noise in qubit 2, that is just the product of noise in q1 and noise in q2. So, there are really fewer parameters.

So, instead of this cumbersome representation, we look at the generator for this, which is termed Liouvillian:

For a typical physical system, the Liouvillian has a Lindblad form:

This model really helps us know what the noise is. We can implement this on a simple quantum circuit, get the Lindblad tomogram, and can see the cause of our error.

Lindblad tomogram before the Dynamical Decoupling

Since the Z error is dominating, we can use the Dynamic Decoupling technique to cancel out these Z errors. After the implementation of DD on the idle qubits, that are showing the most Z errors.

Lindblad tomogram after the Dynamical Decoupling.

The dominant Z errors are now completely suppressed. So we can also use this technique to diagnose and optimize our experiment, which is essentially what they did in this paper.

Can we scale this to bigger circuits?

Yes! We can. As shown in the figure below a 20 qubit layer tomogram (taken from the slides of Zlatkov Minev)

20 qubit layer tomogram

The circled pairs represent a CNOT gate, like 1–4 has a CNOT gate, 7–10 has one and so on. There are single-qubit errors on each qubit which are characterized by the circles, for instance, we can see that qubit 1 has a lot of Z error, same with qubit 18 and qubit 5 and so on. Qubit 23 is dominated by Y error. The same goes for the two body terms, they are shown in boxes, like the XX, XZ, ZX and all others. They all are colour coded with respect to the amplitude of the error they have.

Ising Model Trotterized Time Evolution

Now let’s look at the Ising model, which was used in the paper, which generally simulates a spin system.

Source: Quanta Magazine

We can simulate the spin chain by trotterize the Hamiltonian, which will result in circuits like these:

Now for the error mitigation part, the first step is to take the circuit and amplify the noise ( as discussed above)

Then for each of the layers create a Lindblad sparse tomogram:

Then we use the data from this Lindblad sparse tomogram to randomly sample the inverse noise map and then to average out the circuit to recover the data.

Zero-noise extrapolation requires a controlled amplification of the intrinsic hardware noise by a known gain factor G to extrapolate back to the ideal G=0 result. (The same parameter is called λ in the ZNE paper and my blog article). More precise noise amplification can, however, enable substantial reductions in the bias of the extrapolated estimator.

Noise Cancelling Headphones: Effectively we are learning the noise, then tailoring the noise as per our requirements, and then cancelling noise with noise. (image source: Taken from Zlatkov Minev’s slides from Qiskit Lecture Series).

The sparse Pauli–Lindblad noise model proposed above turns out to be especially well suited for noise shaping in ZNE. It will be similar to that of PEC. The Pauli errors inserted at proportional rates can be used to either cancel (PEC) or amplify (ZNE) the intrinsic noise.

In PEC, we choose to obtain an overall zero-gain noise level. In ZNE, we instead amplify the noise to different gain levels and estimate the zero-noise limit using extrapolation.

Clifford circuits serve as useful benchmarks of estimates produced by error mitigation, as they can be efficiently simulated classically. They are a specific class of quantum circuits that play a fundamental role in quantum computing. These circuits consist entirely of Clifford gates, which are a set of quantum gates that can be efficiently simulated classically.

The Clifford gate set includes gates such as the Hadamard gate (H), Pauli gates (X, Y, Z), phase gates (S and S†), and the controlled-NOT gate (CNOT). These gates form a universal set in combination with the single-qubit Pauli gates (X, Y, Z) and the CNOT gate.

The entire Ising Trotter circuit becomes Clifford when θh is chosen as a multiple of π/2. As a first example, they set the transverse field to zero (RX(0) = I) and evolve the initial state |0⟩⊗127. As shown here:

After running the circuit for four trotter steps, that makes it 12 CNOT layers, and the noise gain level G ∈ {1, 1.2, 1.6}. For each G, they generated 2,000 circuit instances in which, before each layer l, they inserted products of one-qubit and two-qubit Pauli errors.

Zero Noise Extrapolation with Probabilistic error amplification. Mitigated expectation values from Trotter circuits at the Clifford condition.

Convergence of the unmitigated (G=1), noise amplified d (G > 1) and noise mitigated (ZNE) estimates of ⟨Z_106⟩ after four Trotter steps. In all panels, error bars indicate 68% confidence intervals obtained by means of percentile bootstrap. Exponential extrapolation (exp, dark blue) tends to outperform linear extrapolation (linear, light blue) when differences between the converged estimates of ⟨Z_106⟩_G≠0 are well resolved.

The procedure outlined in the above figure was applied to the measurement results from each qubit q to estimate all N = 127 Pauli expectations ⟨Z_q⟩_0

The variation of unmitigated and mitigated observables (as shown in the figure below) is an indication of the non-uniformity in the error rates of the entire processor. Although the unmitigated results show gradual decay from 1 with an increasing deviation for deeper circuits.

Magnetization is computed as the mean of the individual estimates of <zq> for all qubits. As circuit depth is increased, unmitigated estimates of the average magnetization decay monotonically from the ideal value of 1. ZNE greatly improves the estimates even after 20 Trotter steps or 60 CNOT depth.

Now the method was tested on non-Clifford circuits, with non-trivial entangling dynamics. They are very important to test as the validity of the exponential extrapolation is no longer guaranteed. In this case, we restricted the circuit to 5 Trotter steps, which is equivalent to 15 CNOT layers and wisely chose the observables that are exactly verifiable.

The figure below shows the result as the θh is swept between 0 and π/2 for three such observables of increasing weight. The first plot shows Mz as before, an average of weight-1<Z> observable.

Expectation value estimates for θh sweeps at a fixed depth of five Trotter steps for the circuit. The considered circuits are non-Clifford except at θh = 0, π/2. Light-cone and depth reductions of respective circuits enable exact classical simulation of the observables for all θh. For all three plotted quantities (panel titles), mitigated experimental results (blue) closely track the exact behaviour (grey). In all panels, error bars indicate 68% confidence intervals obtained by means of percentile bootstrap.

whereas the next shows the weight-10 and weight-17 observable.

The weight-10 and weight-17 observables are stabilizers of the circuit at θh = π/2 with respective eigenvalues +1 and −1; all values in c have been negated for visual simplicity. The lower inset a depicts a variation of ⟨Z q⟩ at θh = 0.2 across the device before and after mitigation and compares with exact results. Upper insets in all panels illustrate causal light cones, indicating in blue the final qubits measured (top) and the nominal set of initial qubits that can influence the state of the final qubits (bottom). M z also depends on 126 other cones besides the example shown. Although in all panels exact results are obtained from simulations of only causal qubits, we include tensor network simulations of all 127 qubits (MPS, isoTNS) to help gauge the domain of validity for those techniques, as discussed in the main text. isoTNS results for the weight-17 operator in c are not accessible with current methods

Although the entire 127-qubit circuit is executed experimentally, light-cone and depth-reduced (LCDR) circuits enable brute-force classical simulation of the magnetization and weight-10 operator at this depth.

Over the full extent of the θh sweep, the error-mitigated observables show good agreement with the exact evolution(see the above figures). However, for the weight-17 operator, the light cone expands to 68 qubits, a scale beyond brute-force classical simulation, so we turn to tensor network methods.

Tensor Networks

They have been used widely to approximate and compress quantum state vectors that arise in the study of low-energy eigenstates of and time evolution by local Hamiltonians, and also have been shown successfully simulate low-depth noisy quantum circuits. Simulation accuracy can be improved by increasing the bond dimension χ, which constrains the amount of entanglement of the represented quantum state, at a computational cost scaling polynomially with χ. As the entanglement (bond dimension) of a generic state grows linearly (exponentially) with time evolution until it saturates the volume law, deep quantum circuits are inherently difficult for tensor networks.

In this work both the 1D matrix product states (MPS) and the 2D isometric tensor network states (isoTNS) that have O (χ³ )and O (χ⁷) scaling of time-evolution complexity, respectively are used. For the case of weight-17 observable, an MPS simulation of the LCDR circuit at χ = 2,048 is sufficient to obtain the exact evolution. This comparison suggests that the domain of experimental accuracy could extend beyond the scale of the exact classical simulation.

Finally, the experiment is stretched to regimes where the exact solution is not available with the classical methods that we discussed above. As shown in the figure below, we do something similar to the one we did before, but we add a further final layer of single-qubit Pauli rotations that interrupt the circuit-depth reduction that previously enabled exact verification for any θh.

At the verifiable Clifford point θh = π/2, the mitigated results agree again with the ideal value, whereas the χ = 3,072 MPS simulation of the 68-qubit LCDR circuit markedly fails in the strongly entangling regime of interest.

Estimating expectation values beyond exact verification. Plot markers, confidence intervals and causal light cones appear as defined in the previous plot, Estimates of a weight-17 observable (panel title) after five Trotter steps for several values of θh.

As a final example, they extended the circuit depth to 20 Trotter steps
(60 CNOT layers) and estimate the θh dependence of a weight-1 observable, ⟨Z_62⟩, as shown in the figure below.

Despite the greater depth, the MPS simulations of the LCDR circuit agree well with the experiment in the weakly entangling regime of small θh. Although deviations from the experimental trace emerge with increasing θh, we note that the MPS simulations slowly move in the direction of the experimental data with increasing χ.

This tensor network was run on a 64-core,2.45 GHz processor with 128 GB of memory, in which the run time to access an individual data point at fixed θh was 8 h for the first plot and 30h for the second plot. The corresponding quantum wall-clock run time was approximately 4 h for the first and 9.5 h for the second.

The observation that a noisy quantum processor, even before the advent of fault-tolerant quantum computing, produces reliable expectation values at a scale beyond 100 qubits and non-trivial circuit depth leads to the conclusion that there is indeed merit to pursuing research towards deriving a practical computational advantage from noise-limited quantum circuits.

But wait!

Just after 2 weeks after this paper was published in Nature, a paper on Efficient tensor network simulation of IBM’s kicked Ising experiment was made available on arxiv.

Link to the paper: https://arxiv.org/abs/2306.14887

They showed an accurate, memory and time-efficient classical simulation of a 127-qubit kicked Ising quantum system on the heavy-hexagonal lattice. They showed that, by adopting a tensor network approach that reflects the qubit connectivity of the device, we can perform a classical simulation that is significantly more accurate than the results obtained from the quantum device in the verifiable regime and comparable to the quantum simulation results for larger depths. The tensor network approach used will likely have
broader applications for simulating the dynamics of quantum systems with tree-like correlations.

Qubit layout for the Eagle quantum processor and the tensor network structure used for this simulation.
The structure directly reflects the processor layout. On-site tensors Γv are coloured in blue and possess physical, uncontracted, indices of dimension 2 (represented by their dangling legs) and virtual indices of dimension χ (represented by the edges of the network) which are shared with neighboring tensors. Positive, diagonal bond tensors Λe live on the edges e between the site tensors and are coloured in grey.

In the figure below, they have compared the results from the IBM paper. The availability of exact results in this case, based on brute force light-cone simulation techniques, allows us to directly assess the errors. The tensor network state (TNS) and optimization methods used to simulate the full 127-qubit system and result in highly accurate expectation values — the expected value of the average single-site magnetization can be computed to an accuracy of ∼ 10-¹⁴ with a simulation that runs in less than 2 minutes, In this classically verifiable regime we systematically outperform the quantum processor results reported in the IBM paper to orders of magnitude and with very modest resources.

The result from the IBM paper:

The result from this paper:

Comparison for classically verifiable systems of our tensor network state using the belief propagation approximation (BP-TNS) approach to simulating the dynamics of the 2D transverse field Ising model versus the Eagle quantum processor and alternative tensor network methods. Expectation values with respect to the state |ψ(θh, 5)⟩ (i.e. following 5 Trotter steps of the dynamics of the model), alongside exact results
determined from light-cone simulations. a) Average magnetization. b) Weight-10 observable. c) Weight-17 observable. The bottom plots show errors defined as the absolute difference between the simulation result and the exact result. For some data points the error from our BP-TNS simulation is too small to fit on the scale of the plot and so these points are not marked. Circled, annotated points denote, for a given θh, the memory required to store the final state and the wall time associated with performing the simulation and calculating the relevant observable on a single Intel Skylake CPU.

Specifically, for every expectation value obtained by the processor and plotted in IBM paper, they obtained a corresponding value to orders of
magnitude better accuracy with a simulation that takes less than 7 minutes to run on a single Intel Skylake CPU and a state that takes up, at most, 0.3GB of memory.

In the next figure, they compared their method to IBM's result, showing properties after larger numbers of Trotter evolution steps. In this regime, the increased depth of the circuit means exact data is no longer available for benchmarking. For both the n = 6 and n = 20 Trotter step simulations they captured the known exact values at the Clifford points θh = 0 and θh = π/2 — with the weight-17 stabilizer demonstrating the anticipated saturation to unity that is not captured by the matrix product state (MPS) approach.

Comparison for non-classically-verifiable systems of our tensor network state using the belief propagation
approximation (BP-TNS) approach to simulating the dynamics of the 2D transverse field Ising model versus the Eagle quantum processor and alternative tensor network methods. Expectation values calculated following a number of Trotter steps of the dynamics of the model are plotted. a) Weight-17 stabilizer after 6 steps
of evolution. b) Weight-1 observable after 20 steps. c) Top and bottom correspond to scaling, with the inverse bond dimension of the heavy-hex TNS ansatz, of the observables in a) and b) respectively at select values of θh.

The accuracy of their results can be understood from the locally tree-like nature of the heavy-hex lattice, with the smallest loop having a path length of 12. Their chosen methods for performing the time evolution and taking expectation values are optimal if the network is a tree. The heavy-hex lattice is close enough to a tree for such methods to still be highly effective.

Their work opens up new directions in which highly flexible and computationally inexpensive tensor network approaches can be used to benchmark new quantum processor designs and can better delineate which many-body quantum systems could become difficult for classical computing techniques.

For the code and data for the IBM paper click here.

The plots from the code are:

Expectation value estimates for θh sweeps at a fixed depth of five Trotter steps for the circuit. The weight-10 and weight-17 observables are stabilizers of the circuit at θh = π/2 with respective eigenvalues +1 and −1; all values in c have been negated for visual simplicity
Estimates of a weight-17 observable (panel title) after 20 trotter steps and five Trotter steps for several values of θh.

Although the simulation presented by the paper is very complex and requires very heavy hardware to run it. We’ll try a simple small-scale simulation of the same model using Pennylane.

Instead of 127 qubits, we will only work with 9 qubits, which means the spins will be placed in a 3x3 grid with nearest neighbour interaction. The circuits for the time evolution will be set in DepolarizingChannel which will be the cause of the error.

# import basic pennylane and jax
import pennylane as qml
import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt
jax.config.update("jax_enable_x64", True)
jax.config.update('jax_platform_name', 'cpu')

Decide on the number of qubits and add the noise channel

n_wires = 9
# Describe noise
noise_gate = qml.DepolarizingChannel
p = 0.005

Load the device, one error-free, one error-prone

# Load devices
dev_ideal = qml.device("default.mixed", wires=n_wires)
dev_noisy = qml.transforms.insert(noise_gate, p, position="all")(dev_ideal)

How the connections of qubit will be?

# 3x3 grid with nearest neighbors
connections = [(0, 1), (1, 2),
(3, 4), (4, 5),
(6, 7), (7, 8),
(0, 3), (3, 6),
(1, 4), (4, 7),
(2, 5), (5, 8)]

Define the time evolution operator, which returns the expectation value

def time_evolution(theta_h, n_layers = 10, obs = qml.PauliZ(4)):
for _ in range(n_layers):
for i, j in connections:
qml.IsingZZ(-jnp.pi/2, wires=(i, j))
for i in range(n_wires):
qml.RX(theta_h, wires=i)
return qml.expval(obs)
qnode_ideal = qml.QNode(time_evolution, dev_ideal, interface="jax")
qnode_noisy = qml.QNode(time_evolution, dev_noisy, interface="jax")

We can now simulate the final expectation value, both with noise and without noise

thetas = jnp.linspace(0, jnp.pi/2, 50)

res_ideal = jax.vmap(qnode_ideal)(thetas)
res_noisy = jax.vmap(qnode_noisy)(thetas)

plt.plot(thetas, res_ideal, label="exact")
plt.plot(thetas, res_noisy, label="noisy")
plt.xticks([0, jnp.pi/8, jnp.pi/4, 3*jnp.pi/8, jnp.pi/2], ["0", "$\\pi$/8", "$\\pi/4$", "$3\\pi/4$", "$\\pi/2$"])
plt.xlabel("$\\theta_h$")
plt.ylabel("$\\langle Z_4 \\rangle$")
plt.legend()
plt.show()
We see that the fidelity of the result is decreased by the noise.

Next, we show how this noise can be effectively mitigated.

Error Mitigation via ZNE

First learn the parameters of the noise model, which is the Pauli Lindblad model. The noise model of our simulation is relatively simple and we have full control over it. This means that we can simply attenuate the noise of our model by an appropriate gain factor. Here, G=(1,1.2,1.6) in accordance with the IBM paper

dev_noisy1 = qml.transforms.insert(noise_gate, p*1.2, position="all")(dev_ideal)
dev_noisy2 = qml.transforms.insert(noise_gate, p*1.6, position="all")(dev_ideal)

qnode_noisy1 = qml.QNode(time_evolution, dev_noisy1, interface="jax")
qnode_noisy2 = qml.QNode(time_evolution, dev_noisy2, interface="jax")

res_noisy1 = jax.vmap(qnode_noisy1)(thetas)
res_noisy2 = jax.vmap(qnode_noisy2)(thetas)

We can take these results and simply extrapolate them back to G=0 with a polynomial fit. We can visualize this by plotting the noisy, exact and extrapolated results.

Gs = jnp.array([1., 1.2, 1.6])
y = jnp.array([res_noisy[0], res_noisy1[0], res_noisy2[0]])
coeff = jnp.polyfit(Gs, y, 2)
x = jnp.linspace(0, 1.6, 100)

plt.plot(x, jnp.polyval(coeff, x), label="fit")
plt.plot(Gs, y, "x", label="noisy results")
plt.plot([0], res_ideal[0], "X", label="exact result")
plt.xlabel("noise gain G")
plt.ylabel("$\\langle Z_4 \\rangle$")
plt.legend()
plt.show()

This was for a single value of θh, we will repeat it for all values of it and use the richardson_extrapolate() to see our performance:

res_mitigated = [qml.transforms.richardson_extrapolate(Gs, [res_noisy[i], res_noisy1[i], res_noisy2[i]]) for i in range(len(res_ideal))]

plt.plot(thetas, res_ideal, label="exact")
plt.plot(thetas, res_mitigated, label="mitigated")
plt.plot(thetas, res_noisy, label="noisy")
plt.xticks([0, jnp.pi/8, jnp.pi/4, 3*jnp.pi/8, jnp.pi/2], ["0", "$\\pi$/8", "$\\pi/4$", "$3\\pi/4$", "$\\pi/2$"])
plt.xlabel("$\\theta_h$")
plt.ylabel("$\\langle Z_4 \\rangle$")
plt.legend()
plt.show()

The observation that a noisy quantum processor, even before the advent of fault-tolerant quantum computing, produces reliable expectation values at a scale beyond 100 qubits and non-trivial circuit depth leads to the conclusion that there is indeed merit in pursuing research towards deriving a practical computational advantage from noise-limited quantum circuits.

Although we were still able to find a better classical simulation that was able to match and even outperform this quantum solution.

Ref:

[1] : Evidence for the utility of quantum computing before fault tolerance

[2] :van den Berg, E., Minev, Z.K., Kandala, A. et al. Probabilistic error cancellation with sparse Pauli–Lindblad models on noisy quantum processors. Nat. Phys. https://doi.org/10.1038/s41567-023-02042-2 (2023).

[3]: Temme, K., Bravyi, S. & Gambetta, J. M. Error mitigation for short-depth quantum circuits. Phys. Rev. Lett. 119, 180509 (2017).

[4]: Efficient tensor network simulation of IBM’s kicked Ising experiment

[5] Is quantum computing useful before fault tolerance?- Pennylane Blog

--

--