Random circuit block-encoded matrix and a proposal of quantum LINPACK benchmark

Quantum computers hold the promise of dramatically accelerating calculations in a wide range of fields, and quantum supremacy was achieved in 2019 via sampling random quantum circuits Arute et al. (2019). Assume that there are ten thousand quantum computers (or many more) available now, how should we select the top 500 best performing computers for scientific computing applications? The answer in the context of classical supercomputers is given by the LINPACK benchmark Dongarra et al. (2003), which measures the floating point computing power of a classical computer via its performance for solving linear systems of equations $Ax=b$. The input matrix $A$ is a dense pseudo-random matrix, and there is no immediate application associated with such a matrix (similar to a quantum supremacy experiment in this sense). There has been much controversy over its effectiveness in measuring the capability of classical computers in scientific computing applications since the very beginning. However, LINPACK is widely used and performance numbers are available for almost all relevant systems. The LINPACK benchmark has been used as the defining criterion of TOP500 supercomputers since the debut of the list in 1993 TOP . One important reason is that dense matrices and dense matrix operations are relatively easy to implement and to optimize on classical computers. These operations have been tuned to be highly scalable, which enabled the performance benchmark of systems that cover a performance range of 12 orders of magnitude in the past 20 years.

In order to mimic the success of the LINPACK benchmark on quantum computers, we consider the problem of solving the quantum linear system problem (QLSP). Many challenging high-dimensional problems in physics, such as computing Green’s functions for a quantum many-body system, can be formulated in terms of QLSP. This field has witnessed significant progresses in the past few years Harrow et al. (2009); Cao et al. (2013); Childs et al. (2017); Chakraborty et al. (2018); Gilyén et al. (2019); Subaşı et al. (2019); Wossnig et al. (2018); An and Lin (2019); Lin and Tong (2019); Xu et al. (2019); Bravo-Prieto et al. (2019); Casares and Martin-Delgado (2019); Tong et al. (2020). Shortly speaking, given $A\in\mathbb{C}^{2^{n}\times 2^{n}}$ and $\ket{b}\in\mathbb{C}^{2^{n}}$, QLSP aims at obtaining an $n$-qubit solution vector $\ket{x}\propto A^{-1}\ket{b}$. More precisely and using the language of block-encoding Gilyén et al. (2019), QLSP is the problem of finding an $(m+n)$-qubit unitary matrix $U$, such that

In other words, the solution is obtained upon measuring $0$ for all $m$ ancilla qubits (with a success probability $p)$.

In this paper, we propose the quantum LINPACK benchmark, which can be concisely stated as the problem of using the quantum computer to evaluate the success probability $p$ for a certain random matrix $A$. While the rationale of such a task will be discussed in detail in the main text, we first emphasize that the quantum LINPACK benchmark examines the whole machine performance of quantum computers, rather than the performance of a few qubits as often measured by methods such as randomized benchmarking Magesan et al. (2011) and gateset tomography Blume-Kohout et al. (2017). There have been a number of whole machine quantum benchmarks proposed in the literature, such as quantum volume Cross et al. (2019), cycle benchmarking Erhard et al. (2019), and the linear cross-entropy benchmarking as in the supremacy experiment Boixo et al. (2018); Arute et al. (2019). However, those benchmark methods are proposed for generic settings and therefore they are less representative for the performance on the user-specified applications. In particular, the error of a structured circuit can be significantly larger than that of a fully randomized one Proctor et al. (2020). The quantum LINPACK benchmark targets directly at the performance of quantum computers for scientific computing applications, as in the case of the LINPACK benchmark for classical supercomputers. We emphasize that the success of the quantum LINPACK benchmark does not guarantee that the quantum computer has solved $\ket{x}$ accurately, but should be viewed as the minimal requirement for a quantum computer to solve QLSP. Given the wide range of potential applications of QLSP from quantum many body problems to quantum machine learning, it is important for future quantum computers to first meet the criterion of the quantum LINPACK benchmark, in order to achieve quantum advantage via the path of solving linear algebra problems.

In order to perform the quantum LINPACK benchmark, note that it would be highly inefficient if we first generate a dense pseudo-random matrix $A$ classically and then feed it into the quantum computer using e.g. QRAM Giovannetti et al. (2008). In fact, such a strong assumption on the input model often lead to dequantized classical algorithms Tang (2019). Instead we focus on matrices that are inherently easy to generate on quantum computers. In particular, the supremacy experiment inspires us to generate a random matrix directly using a random quantum circuit.

We propose an input model called the RAndom Circuit Block-Encoded Matrix (RACBEM). We argue that the RACBEM model is a proper generalization of dense random matrices in the quantum setting, suitable for linear algebra tasks. The RACBEM model, and its Hermitian version called H-RACBEM model, are simple to construct and allow us to get access to in principle any $n$-qubit matrix and $n$-qubit Hermitian matrix, respectively (up to a scaling factor) by adding only one ancilla qubit. Together with the recently developed technique of quantum singular value transformation (QSVT) Gilyén et al. (2019), we yield a practical algorithm for performing the quantum LINPACK benchmark on near-term devices with a shallow circuit depth.

With QSVT and the H-RACBEM model, the circuit used in the quantum LINPACK benchmark can be designed to adapt to the coupling map of almost any given gate-based quantum architecture. All operations can be carried out with straightforward usage of basic one-qubit gates and CNOT gates, and there is no complex controlled unitaries involved. Due to the use of the basic gate set and the adaptivity to the quantum architecture, the quantum LINPACK benchmark does not require the explicit use of the compiler, while the randomized benchmarking requires gate compilation. Furthermore, the number of ancilla qubits needed is minimal (usually 2). By using the H-RACBEM model, the condition number of the random matrices is fully controllable which is crucial for reducing the circuit depth. While the development of quantum hardware has been very rapid in the past few years, quantum resources are expected to remain costly and limited for some time, with or without the fault-tolerant capability. The H-RACBEM model can also significantly reduce the efforts needed to optimize and to compile the QSVT algorithm. The RACBEM and its Hermitian version provides a simple and reliable way to generate random matrices on a quantum computer with the minimal use of the ancilla qubits and the controllability to the circuit depth. Therefore, we expect that the quantum LINPACK benchmark uses the minimal circuit to gauge the performance of a quantum computer for solving linear algebra problems.

Furthermore, we demonstrate that using the same quantum circuit but with different parameters, the combination of QSVT and the H-RACBEM model can be used to perform many other linear algebra tasks, such as computing spectral measures and performing time series analysis (without Trotter splitting). Using the minimally entangled typical thermal state (METTS) algorithm White (2009); Stoudenmire and White (2010); Motta et al. (2020), we also show how H-RACBEM simplifies the computation of the thermal average of the energy. These linear algebra tasks can also be used to construct benchmarks similar to the quantum LINPACK benchmark, whose performance reflects the minimal requirement for a quantum computer in solving corresponding scientific computing problems. We implement all algorithms on the IBM Q quantum architecture. Due to the high noise level of the currently available quantum architecture, we also demonstrate the numerical performance using quantum virtual machines (QVM) with tunable, approximate error models derived from quantum devices.

The block-encoding model also has its limitation. Take an $n$-qubit, $d$-sparse matrix $A$ (i.e. the number of nonzero entries in each row/column does not exceed $d$) for example, assuming access to its row/column entries via certain sparse-access oracles, then we can construct a block-encoding using $n+3$ ancilla qubits Gilyén et al. (2019). Any further manipulation of $U_{A}$, such as quantum signal processing, would require using $(n+4)$-qubit Toffoli gates, which are relatively expensive to implement Barenco et al. (1995); Saeedi and Pedram (2013). Another example is that $A$ is given by a linear combination of $K$ terms, each term being a tensor product of Pauli matrices. In the setting of the Hamiltonian simulation, $e^{\mathrm{i}At}$ can be simply approximated by exponentiating each term individually following a certain order via the Trotter-Suzuki formula. However, the block-encoding model would essentially require using linear combination of unitaries Berry et al. (2015), which not only requires $\lceil\log K\rceil$ ancilla qubits, but also usage of $(\lceil\log K\rceil+1)$-qubit Toffoli gates to implement the prepare and select oracles needed to obtain the linear combination. Such operations are essentially forbidden for near-term applications due to high error rates, and are still challenging when the number of qubits and the gate depth remain a limitation in reaching the desired accuracy on fault-tolerant devices.

It turns out that any matrix $A$ with $\left\lVert A\right\rVert_{2}\leq 1$ can be given by a $(1,1,0)$-block-encoding. Consider any $n$-qubit matrix with its singular value decomposition (SVD) $A=W\Sigma V^{{\dagger}}$, where all singular values in $\Sigma$ belong to $[0,1]$. Then we may construct an $(n+1)$-qubit unitary matrix

which is a $(1,1,0)$-block-encoding of $A$. Since a random circuit with $\text{poly}(n)$ depth can approximate an $n$-qubit Haar measure at least according to the criterion of the $2$-design Harrow and Low (2009), a sufficiently general $(n+1)$-qubit unitary $U_{A}$ can give access to in principle any $n$-qubit non-unitary matrices $A$ (up to a scaling factor). Furthermore, such a random circuit $U_{A}$ can be constructed using only basic one-qubit unitaries and CNOT gates. The matrix $A$ obtained by measuring the first qubit (or in fact, any qubit used as the ancilla) is called a RACBEM. Since the Haar measure is the uniform distribution of unitary matrices, we conclude that RACBEM is a proper generalization of dense matrices on quantum computers suitable for performing linear algebra tasks. The layout of the two-qubit operations can be designed to be compatible with the coupling map of the hardware.

For instance, for the ibmq_burlington device with its coupling map shown in Fig. 2, we can choose qubit 1 as the RACBEM ancilla qubit, which results in a RACBEM model with 3 system qubits $2,3,4$. The qubit 0 is not used here, and is reserved as the signal qubit for quantum singular value transformation to be discussed later.

Consider the quantum circuit in Fig. 4 denoted by $U_{\mathfrak{H}}$, where $\varphi_{0},\varphi_{1}\in[-\pi,\pi)$. Direct calculation (see Appendix A) shows that

Here $\mathfrak{H}$ is a Hermitian matrix. We refer to it as a Hermitian RACBEM (H-RACBEM), and $U_{\mathfrak{H}}$ is its $(1,2,0)$-block-encoding. In particular, choosing $\varphi_{0}=\pi/8$, $\varphi_{1}=-\pi/4$, then $\mathfrak{H}=A^{{\dagger}}A$ is Hermitian positive semi-definite. This will be referred to as a canonical H-RACBEM. In other words, a canonical H-RACBEM is constructed from its (non-unique) matrix square root $A$.

In Fig. 4, the CNOT gate controlling on $0$ instead of $1$ is mainly for notational convenience, and in fact not all CNOT gates are necessary here. For example, in order to implement a canonical H-RACBEM, we only need 1 application of $U_{A}$, 1 application of $U_{A}^{{\dagger}}$, 2 H gates, 2 standard CNOT gates, 1 $\mathrm{S}^{{\dagger}}$ gate, and 2 T gates (see Appendix A). Since any matrix with singular values bounded by $1$ can be represented as a RACBEM, we immediately have that any Hermitian positive semi-definite matrix with eigenvalues bounded by $1$ can be represented as a canonical H-RACBEM, with a sufficiently flexible $U_{A}$.

Therefore Fig. 4 implements a QSVT for a second order polynomial $\mathfrak{H}=h^{\triangleright}(A)$ with a symmetric choice of phase factors $\Phi=(\varphi_{0},\varphi_{1},\varphi_{0})$, where

A canonical H-RACBEM is given by $h(x)=x^{2}$.

Consider any real polynomial $g(x)$ of degree $d$ without parity constraint, satisfying $\left\lvert g(x)\right\rvert\leq 1$ for any $x\in[-1,1]$. Then using the identity

any matrix function $g(\mathfrak{H})$ can be expressed as a QSVT with respect to an even polynomial $f=g\circ h$ of degree $2d$. We remark that when $g$ does not have a definite parity, the associated QSVT of $A$ is much more involved. It generally requires using a linear combination of block-encoding of the even and odd parts, which in turn requires implementing controlled $U_{A}$ Gilyén et al. (2019). Normally such controlled operations have significant overhead. For instance, if we would like to implement a controlled-RACBEM, generally we need to convert all quantum gates in the circuit of $U_{A}$ to the controlled version.

By expressing $g(\mathfrak{H})$ as a QSVT associated with an even polynomial, we have not only eliminated the need of controlled unitary operations, but also saved one additional qubit. This is because if we first construct a H-RACBEM $\mathfrak{H}$ and then construct $g(\mathfrak{H})$, we need $3$ ancilla qubits in total, and possibly a LCU circuit when $g$ does not have a definite parity. On the other hand, by considering the composite function $f$, the parity constraint on $g$ is completely removed, and we only need $2$ ancilla qubits (the same as that needed for a H-RACBEM). Furthermore, each controlled gate in Fig. 5 is a standard CNOT gate rather than a Toffoli gate. Therefore the H-RACBEM model has a salient advantage and significantly simplifies the implementation.

In many applications, such as a LCU circuit and the Hadamard test, we do need to get access to controlled $U_{\mathfrak{H}}$ or $U_{g(\mathfrak{H})}$. In such a case, the fact that $\mathfrak{H}$ is a H-RACBEM is also helpful: Note that if we remove all Z-rotations in the first line of Fig. 5, the circuit implements an identity operator since $U_{A},U_{A}^{{\dagger}}$ always appear in pairs thanks to that $f=g\circ h$ is an even polynomial. Therefore a controlled version $U_{g(\mathfrak{H})}$ can be simply implemented by controlling all the Z-rotation gates Gilyén et al. (2019).

The RACBEM, as well as the H-RACBEM model provides a solution to the read-in problem using only basic quantum gates, and we can design them to be optimally adapted to the hardware architecture without resorting to complex quantum compilers. Hence they can be regarded as the proper generalization of “test dense matrices” in the quantum setting.

In this section, we demonstrate that the usage of the H-RACBEM model for solving QLSP. We assume $A=\mathfrak{H}$ is a H-RACBEM, and without loss of generality we may take $\ket{b}=\ket{0^{n}}$. Let the condition number of $\mathfrak{H}$ be denoted by $\kappa$, which is the ratio between the maximum and the minimum of the singular values of $\mathfrak{H}$. It is believed that the computational complexity for solving QLSP cannot generally be better than $\mathcal{O}(\kappa^{1-\delta})$ for any $\delta>0$ Harrow et al. (2009), and the cost of using QSVT to solve general linear systems scales as $\widetilde{\mathcal{O}}(\kappa^{2}\log(1/\epsilon))$ Gilyén et al. (2019). So the treatment of ill-conditioned matrices is very difficult especially on near-term devices. To reduce the circuit depth, in the near future it may be more productive to focus on solving well conditioned linear systems.

Note that if $A$ has at least one singular value that is zero (or near zero), a canonical H-RACBEM $\mathfrak{H}$ is not invertible (or very ill-conditioned). Such events can occur more frequently as $n$ becomes large. It can be difficult to diagnose such a problem without first obtaining some spectral information of $A$, which is perhaps a more difficult task than solving the linear system problem itself.

The H-RACBEM model offers a new and natural way to solve this problem. From Eq. 4, assuming $\cos(2\varphi_{0}-\varphi_{1})>0$ and $-2\sin(2\varphi_{0})\sin\varphi_{1}>0$, and use that $0\preceq A^{{\dagger}}A\preceq 1$, the condition number of $\mathfrak{H}$ can be bounded from above:

Therefore the condition number of a H-RACBEM is fully tunable by changing the phase factors $\varphi_{0},\varphi_{1}$. According to Eq. 6, this changes the second order polynomial function $h(x)$, so that $\mathfrak{H}:=h^{\triangleright}(A)\succ 0$ has a tunable, bounded condition number.

In order to solve QLSP, we are interested in finding a polynomial $g(x)$ of degree $d$ so that

Following (Gilyén et al., 2018, Corollary 69), there can be satisfied by an odd polynomial $g(x)$ with degree $d\sim\mathcal{O}(\kappa\log(1/\epsilon))$, which gives an upper bound on $d$. Numerical results shows that a better approximation to $g(x)$ can be obtained by solving a minimax problem using the Remez algorithm Dong et al. (2020), and the polynomial can be chosen to be either even or odd. Fig. 6 shows the shape of the optimal polynomial even/odd approximation to $x^{-1}$ on the interval $[\kappa^{-1},1]$ found by the Remez algorithm to reach a target accuracy of $10^{-3}$ with $\kappa=10$. In particular, the usage of an even polynomial can further reduce the polynomial degree, which will be used in the numerical tests below.

We may then construct a degree $2d$ polynomial $f=g\circ h$ as in Section II. Once we find the associated phase factors, the circuit in Fig. 5 implements $g(\mathfrak{H})\ket{b}$, which satisfies the error bound $\left\lVert g(\mathfrak{H})\ket{b}-\kappa^{-1}\mathfrak{H}^{-1}\ket{b}\right\rVert_{2}\leq\epsilon$ for any normalized vector $\ket{b}$. The success probability of measuring the ancilla qubits and obtaining an all $0$ output (will be referred to as the success probability for short) is

which can be of modest size given $\kappa$ is not too large.

Since measuring the accuracy of all entries of the solution is not a practical goal, our proposal of the quantum LINPACK benchmark is to measure the success probability $p$, i.e. the probability of measuring the $2$ ancilla qubits and obtaining $\ket{00}$, and to compare it with the numerically exact probability (denoted by $p_{\mathrm{exact}}$) computed using a classical computer. When $\kappa,\epsilon$ and H-RACBEM are given, the quantum LINPACK benchmark reports the relative error $\left\lvert p-p_{\mathrm{exact}}\right\rvert/p_{\mathrm{exact}}$ to describe the accuracy of a quantum computer for solving QLSP. In the future we may also take into account the wall clock time, or a quantity analogous to the floating point operations per second (FLOPS) for classical computers to measure the efficiency of a quantum computer. The quantum LINPACK benchmark uses only basic single-qubit gates and CNOT gates, and hence is easy to implement even on near-term devices.

Clearly, the success of the quantum LINPACK benchmark is only a necessary condition for the accurate solution of QLSP. However, as will be shown in numerical experiments, this task can already be challenging for near-term devices due to the presence of noise, and therefore the benchmark provides a meaningful and easily implementable criterion for measuring the accuracy of quantum computers. Furthermore, due to the very flexible construction of $U_{A}$ and the near-minimal usage ancilla qubits and other primitive gates, we also expect that the success of quantum LINPACK benchmark is the minimal requirement in order to achieve quantum advantages by solving QLSP or similar linear algebra tasks.

Since QSVT implements a general matrix function, its application is therefore not limited to solving QLSP. Here we demonstrate a few more applications. Throughout the section $\mathfrak{H}$ is a H-RACBEM defined via a unitary matrix $U_{A}$ and $h(x)$ is a second order polynomial.

In order to accelerate the convergence of the polynomial approximation, in practice, we may introduce another tunable parameter $\eta\geq 1$ in the formulation so that

and then

For instance, $\eta$ can be set to $1.0$ or $1.5$ and the details are given in Table A2.

In Somma (2019), the time series analysis is used for eigenvalue estimation, which is implemented via the Hadamard test and requires a controlled Hamiltonian evolution procedure, followed by a Fourier transform. The procedure above removes the need of performing controlled Hamiltonian evolution. We also remark that in terms of quantum estimation of eigenvalues, the method using the spectral measures to be discussed below can be more advantageous, which does not require a subsequent Fourier transform, and the resulting spectral measure is guaranteed to be non-negative.

Unlike Metropolis type algorithms which follows an acceptance / rejection procedure, the METTS algorithm samples the states $\{\ket{\phi_{i}}\}$ as follows. We start from a computational basis state $\ket{i}$ (e.g. $\ket{i}=\ket{0^{n}}$). Then

1. 1.

   Evaluate the contribution to the energy from the state $\ket{\phi_{i}}$ via Eq. 11.

2. 2.

   Collapse $\ket{\phi_{i}}$ to a new state in the computational basis $\ket{i^{\prime}}$ with probability $|\braket{i^{\prime}}{\phi_{i}}|^{2}$, and repeat step 1.

Note that step $2$ of the METTS algorithm is very simple to implement: we only need to construct a QSVT circuit for preparing the unnormalized state $\ket{\widetilde{\phi}_{i}}$, which is readily available when computing the denominator in Eq. 11. Then collapsing to $\ket{i^{\prime}}$ can be implemented by measuring all $n$ system qubits in the computational basis.

The evaluation of $f_{n,\beta}:=g_{n,\beta}\circ h$ requires approximating the square root function. When $\mathfrak{H}$ is a canonical H-RACBEM with $h(x)=x^{2}$, an alternative method is to approximate an odd function $\widetilde{f}_{n,\beta}(x)=xe^{-\beta x^{2}/2}$ instead of $f_{n,\beta}$. Then the complexity for thermal average calculation can be only $\widetilde{\mathcal{O}}(\sqrt{\beta}\log(1/\epsilon))$ (Gilyén et al., 2018, Corollary 64).

The source code of RACBEM is available in the Github repository¹¹1https://github.com/qsppack/RACBEM. We use the optimization-based algorithm which is described in Appendix B to generate phase factors and the source code is available in the Github repository²²2https://github.com/qsppack/QSPPACK. We demonstrate the numerical performance of the RACBEM model for solving various numerical quantum linear algebra tasks. The algorithms are performed on the IBM Q quantum device, as well as QVM with approximate noisy running environment retrieved from quantum devices. All numerical tests are implemented in python3.7 and Qiskit Abraham et al. (2019). In order to construct an adjustable noise model, we retrieve the noise model from real quantum devices provided by IBM Q backends. The magnitude of the noise level is then made to be fully tunable through a single parameter $\sigma$. When $\sigma=0$, the only noise contribution comes from the Monte Carlo error in measurements, which can be systematically reduced by increasing the number of samples. When $\sigma=1$, the noise model contains all the readout errors and quantum errors associated with a probability distribution given by the retrieved noise model (see Appendix C for details). Unless otherwise noted, the number of measurements (shots) is fixed to be $8192$ throughout this section. As the quantum linear algebra tasks are solved via QSVT, the overall error consists of contributions from the polynomial approximation, the noise from the device, and the Monte Carlo sampling error. We remark that in numerical results we distinguish the setup of “QSVT without error” (only taking into account the polynomial approximation error) from that of “$\sigma=0$” (including both the Monte Carlo sampling error and the polynomial approximation error).

We use Algorithm 1 (in Appendix C) to generate custom random quantum circuits with respect to a given coupling map. In each layer of the quantum circuit, we apply a one-qubit (or two-qubit) gate to each qubit (or a pair of qubits) randomly selected from the basic gate set. Though CNOT gates can increase the entanglement among the qubits, their error rate is much higher than that of one qubit gates on IBM Q backends. So we control the number of CNOT gates via a parameter to determine the probability that a CNOT gate is drawn. Unless otherwise noted, this probability is set to be $p=0.5$. The circuit depth is the same as the number of layers in the generating circuit. We would like to emphasize that the circuit depth must not be too small. Otherwise the block-encoded matrix can sometimes become degenerate (such as a scaled identity matrix). For an $n$-qubit system, we empirically set the depth $\ell$ to be $\ell=3$ when $n=1$, $\ell=7$ when $n=2$, and $\ell=15+2(n-3)$ when $n\geq 3$. We report the statistics of singular-value distributions of the block-encoded matrix in Fig. A5 (Appendix E) to justify this adaptive choice of the circuit depth.

According to the random quantum circuit generation algorithm, we measure the total gate count of a quantum circuit by its logical gate count with respect to the basic gate set, which is upper bounded by its depth $\ell$ times the number of system qubits $n$. The available basic gates provided by IBM Q backends are the CNOT gate and U1, U2, U3 gates, which are parameterized families of generic one-qubit gates (see Appendix C for details). In order to reduce the noise due to U3 gates, we restrict the basic gate set in the custom random circuit generator to be {CNOT, U1, U2}, which is still a universal gate set. Each controlled rotation in QSVT circuit costs $7$ logical gates (4 X gates implemented by U2 gates, 2 CNOT gates, and 1 U1 gate). Therefore, given a RACBEM whose depth is $\ell$ and with $n$ system qubits, to implement the QSVT of a real polynomial of degree $d$, the total gate count for the circuit is upper bounded by $2+7(d+1)+d\ell n$. The details about QSVT phase factors used in numerical experiments can be found in Appendix E. Unless otherwise noted, the input quantum state of the quantum circuit is set to $\ket{0^{n+2}}$.

We report the performance of solving linear systems on the IBM Q device using the H-RACBEM model in Fig. 8. The architecture of the five backends are identical, and therefore we may draw a H-RACBEM at random and test it on all five backends. By tuning $h(x)$, the condition number of $\mathfrak{H}$ is upper bounded by $2$. Therefore we may even use a very short QSVT circuit with $d=2$, and the corresponding number of phase factors is $3$. This is essentially a linear approximation to the inverse, and the accuracy measured by the $L^{\infty}$ error is less than $0.03$. In such a case, the total logical gate count is upper bounded by $113$. We can refine the polynomial approximation by using a more accurate, and deeper QSVT circuit with $d=10$ (the number of phase factors is $11$, and the $L^{\infty}$ error is less than $3\times 10^{-5}$). In a case when $d=10$, the total logical gate count is upper bounded by $529$. The results in Fig. 8 indicate that for the shallow circuit, the relative error of the success probability is similar to that observed in Fig. 7, which only measures the success probability of the block-encoding. However, the relative error is significantly larger using the deeper circuit, despite that the QSVT circuit implements a more accurate polynomial approximation to the matrix inverse. Thus, we conclude that when designing the quantum circuit, a proper choice of QSVT phase factors is needed which reflects the tradeoff between the polynomial approximation error and the error caused by the noisy running environment.

In order to further demonstrate how various parameters can affect the accuracy of the QLSP solver, we vary the number of system qubits, the condition number and the noise magnitude, and compute the relative error of success probability under these different settings on QVM. The numerical results are presented in Fig. 9. In all cases, we find that the QLSP solver can perform very well when $\sigma$ is small. The error is only due to the polynomial approximation and the Monte Carlo sampling error. This corresponds to the setting of fault-tolerant quantum devices. However, the accuracy rapidly deteriorates as $\sigma$ increases. In the noisy setup, the error also increases nearly proportionally to the condition number $\kappa$. This is because the polynomial degree $d$ should increase as $\mathcal{O}(\kappa)$ in order to achieve constant accuracy, and so is the circuit depth. The error also increases with respect to the number of system qubits, but the effect is less significant compared to that due to $\kappa$ which increases the circuit depth.

Analogous to the LINPACK benchmark for measuring the performance of classical supercomputers, we have proposed a quantum LINPACK benchmark to measure the performance of current and future quantum computers for scientific computing applications. The quantum LINPACK benchmark solves a well conditioned quantum linear system problem, which is enabled by the Hermitian RAndom Circuit Block-Encoded Matrix (H-RACBEM) input model and the quantum singular value transformation (QSVT). The flexibility provided by the H-RACBEM model also allows us to perform other linear algebra tasks already on near-term devices, such as computing spectral measures, and time series generated by a Hamiltonian simulation. We can also compute the thermal average of the energy, using a quantum version of the minimally entangled typical thermal state (METTS) algorithm.

We perform these linear algebra tasks on IBM Q quantum devices, and quantum virtual machines with a tunable noise model. Although present quantum devices still suffer from high noise levels, it is nonetheless encouraging to observe that the solutions can already be qualitatively obtained in the noisy environment. Among all numerical tests, the thermal average of the energy computed via the METTS algorithm appears to be particularly robust with respect to noise. When designing the quantum circuit, we need to carefully choose the length of QSVT phase factors, by balancing the polynomial approximation error and the additional quantum error caused by the increase of the circuit depth. Our numerical tests are currently limited to matrices up to $10$ qubits (the corresponding matrix size is $1024$) in order to compare with numerically exact results obtained on classical computers. However, the number of qubits can be directly increased to be much larger, especially on future devices with reduced noise level.

Note that Google’s quantum supremacy experiment is a benchmark problem (the linear cross-entropy test) and is also hard for any classical computer. The classical hardness is justified by that sampling a certain random quantum circuit (e.g. one that implements a chaotic quantum evolution) and generating heavy outputs is a difficult task on any classical computer Emerson et al. (2003); Harrow and Low (2009); Nahum et al. (2017); Aaronson and Arkhipov (2011); Bremner et al. (2016); Fujii (2016); Aaronson and Chen (2016); Bouland et al. (2019); Aaronson and Gunn (2019). Although the quantum LINPACK benchmark also involves a supremacy-type random circuits, the success of the benchmark is defined by measuring only two ancilla qubits, and the task could therefore possibly be “classically spoofed” when the block-encoding matrix $U_{A}$ is drawn from a known distribution, such as the Haar measure. We would like to therefore clarify that a quantum benchmark problem does not need to be classically hard in order to provide useful information on the performance of quantum computers. Indeed, we expect that the success of quantum LINPACK benchmark is the minimal requirement in order to achieve quantum advantages via linear algebra tasks such as QLSP. On the other hand, it is possible to go beyond the quantum LINPACK benchmark, and to formulate a more elaborate cross-entropy test related to QLSP that is also classically hard. This requires detailed studied of the statistical properties of truncated unitaries in the block-encoding model and QSVT problem, which will be our future work.

This work was partially supported by a Google Quantum Research Award (Y.D.,L.L.), by the Department of Energy under Grant No. DE-AC02-05CH11231, No. DE-SC0017867, and the Quantum System Accelerator project (L.L.). This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC05-00OR22725. We thank Yunchao Liu, Yu Tong and K. Birgitta Whaley for useful discussions.