text:
\begin{document}

\title{Random circuit block-encoded matrix and a proposal of quantum LINPACK benchmark}

\author{Yulong Dong$^{1,2}$}
\author{Lin Lin$^{3,4}$}
\email[Electronic address: ]{linlin@math.berkeley.edu}

\affiliation{$^1$\BQIC}
\affiliation{$^2$\DeptChem}
\affiliation{$^3$\DeptMath}
\affiliation{$^4$\LBLMath}

\begin{abstract}

The LINPACK benchmark reports the performance of a computer for solving a system of linear equations with dense random matrices. Although this task was not designed with a real application directly in mind, the LINPACK benchmark has been used to define the list of TOP500 supercomputers since the debut of the list in 1993. We propose that a similar benchmark, called the quantum LINPACK benchmark, could be used to measure the whole machine performance of quantum computers. The success of the quantum LINPACK benchmark should be viewed as the minimal requirement for a quantum computer to perform a useful task of solving linear algebra problems, such as linear systems of equations. We propose an input model called the RAndom Circuit Block-Encoded Matrix (\bem), which is a proper generalization of a dense random matrix in the quantum setting. The RACBEM model is efficient to be implemented on a quantum computer, and can be designed to optimally adapt to any given quantum architecture, with relying on a black-box quantum compiler. Besides solving linear systems, the RACBEM model can be used to perform a variety of linear algebra tasks relevant to many physical applications, such as computing spectral measures, time series generated by a Hamiltonian simulation, and thermal averages of the energy. We implement these linear algebra operations on IBM Q quantum devices as well as quantum virtual machines, and demonstrate their performance in solving scientific computing problems.

\end{abstract}

\maketitle

\section{Introduction}

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 \cite{AruteAryaBabbushEtAl2019}. 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 \cite{DongarraLuszczekPetitet2003}, 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 \cite{TOP500}.
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 \cite{HarrowHassidimLloyd2009,CaoPapageorgiouPetrasEtAl2013,ChildsKothariSomma2017,ChakrabortyGilyenJeffery2018,GilyenSuLowEtAl2019,SubasiSommaOrsucci2019,WossnigZhaoPrakash2018,AnLin2019,LinTong2019,XuSunEndoEtAl2019,Bravo-PrietoLaRoseCerezoEtAl2019,CasaresMartin-Delgado2019,TongAnWiebeEtAl2020}. Shortly speaking, given $A\in \CC^{2^n\times 2^n}$ and $\ket{b}\in\CC^{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 \cite{GilyenSuLowEtAl2019}, QLSP is the problem of finding an $(m+n)$-qubit unitary matrix $U$, such that
\begin{equation}
\ket{x}=(\bra{0^m}\otimes I_n)U(\ket{0^m}\otimes\ket{b})=\frac{A^{-1}\ket{b}}{\norm{A^{-1}\ket{b}}}.
\label{eqn:QLSP}
\end{equation}
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 \emph{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 \textit{whole machine performance} of quantum computers, rather than the performance of a few qubits as often measured by methods such as randomized benchmarking \cite{MagesanGambettaEmerson2011} and gateset tomography \cite{Blume-KohoutGambleNielsenEtAl2017}. There have been a number of whole machine quantum benchmarks proposed in the literature, such as quantum volume \cite{CrossBishopSheldonEtAl2019}, cycle benchmarking \cite{ErhardWallmanPostlerEtAl2019}, and the linear cross-entropy benchmarking as in the supremacy experiment \cite{BoixoIsakovSmelyanskiyEtAl2018,AruteAryaBabbushEtAl2019}. 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 \cite{ProctorRudingerYoungEtAl2020}.
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 \cite{GiovannettiLloydMaccone2008}. In fact, such a strong assumption on the input model often lead to dequantized classical algorithms \cite{Tang2019}. Instead we focus on matrices that are \emph{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 (\bem). We argue that the \bem model is a proper generalization of dense random matrices in the quantum setting, suitable for linear algebra tasks. The \bem model, and its Hermitian version called \hbem 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) \cite{GilyenSuLowEtAl2019}, we yield a practical algorithm for performing the quantum LINPACK benchmark on near-term devices with a shallow circuit depth.

With QSVT and the \hbem 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 \hbem 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 \hbem model can also significantly reduce the efforts needed to optimize and to compile the QSVT algorithm. The \bem 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 \emph{the same quantum circuit but with different parameters}, the combination of QSVT and the \hbem 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
\cite{White2009,StoudenmireWhite2010,MottaSunTanEtAl2020}, we also show how \hbem 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.

\section{Random circuit based block-encoding matrix}\label{sec:racbem}

\noindent\textit{Block-encoding:}
Inherently, quantum computers can only handle unitary operators. Hence any non-unitary operators must be encoded in terms of unitary operators. Let $A\in \CC^{N\times N}$ be an $n$-qubit Hermitian matrix ($N=2^n$). If we can find an $(n+1)$-qubit unitary matrix $U_A$ such that
\begin{equation}
U_A=\left(\begin{array}{cc}
{A} & {\cdot} \\
{\cdot} & {\cdot}
\end{array}\right)
\label{eqn:block_encode_exact_matrix}
\end{equation}
holds, i.e. $A$ is the upper-left matrix block of $U_A$, then we may get access to $A$ via the unitary matrix $U_A$ with $A=\left(\langle 0 | \otimes I_n\right) U_A \left( | 0 \rangle \otimes I_n \right)$. Clearly when the operator norm $\norm{A}_2$ is larger than $1$, $A$ cannot be encoded by any unitary $U_A$ as in \eqref{eqn:block_encode_exact_matrix}. Generally if we can find $\alpha, \epsilon \in \mathbb{R}_+$, and an $(m+n)$-qubit matrix $U_A$ such that
\begin{equation}
\Vert A - \alpha \left(\langle 0^m | \otimes I_n\right) U_A \left( | 0^m \rangle \otimes I_n \right) \Vert \leq \epsilon,
\label{eqn:block_encoding}
\end{equation}
then $U_A$ is called an $(\alpha, m, \epsilon)$-block-encoding of $A$. Here $m$ is called the number of ancilla qubits for block-encoding. The block-encoding is a powerful and versatile model, which can be used to efficiently encode density operators, Gram matrices, positive-operator valued measure (POVM), sparse-access matrices, as well as addition and multiplication of block-encoded matrices (we refer to \cite{GilyenSuLowEtAl2019} for a detailed illustration of such constructions).

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 \cite{GilyenSuLowEtAl2019}. 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 \cite{BarencoBennettCleveEtAl1995,SaeediPedram2013}. 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^{\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 \cite{BerryChildsKothari2015}, 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.

\vspace{1em}
\noindent\textit{\bem:}
To harness the power of the block-encoding model and to avoid its pitfalls, we propose the Random circuit based block-encoding matrix (\bem) model as follows. Instead of first identifying $A$ and then finding its block-encoding $U_A$, we reverse this thought process: we first identify a unitary $U_A$ that is easy to implement on a quantum computer, and then ask which matrix can be block-encoded by $U_A$.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.2\textwidth]{racbem_cartoon.png}
\caption{A cartoon illustration of the \bem model.}
\end{center}
\end{figure}
It turns out that any matrix $A$ with $\norm{A}_2\le 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^{\dag}$, where all singular values in $\Sigma$ belong to $[0,1]$. Then
we may construct an $(n+1)$-qubit unitary matrix
\[
\begin{split}
U_A:=&\left(\begin{array}{cc}
A &  W \sqrt{I_{n}-\Sigma^{2}} \\
\sqrt{I_{n}-\Sigma^{2}} V^{\dagger} & -\Sigma
\end{array}\right)\\
=&\left(\begin{array}{cc}
W & 0 \\
0 & I_{n}
\end{array}\right)\left(\begin{array}{cc}
\Sigma & \sqrt{I_{n}-\Sigma^{2}} \\
\sqrt{I_{n}-\Sigma^{2}} & -\Sigma
\end{array}\right)\left(\begin{array}{cc}
V^{\dagger} & 0 \\
0 & I_{n}
\end{array}\right),
\end{split}
\]
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 \cite{HarrowLow2009}, 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 \bem. Since the Haar measure is the uniform distribution of unitary matrices, we conclude that \bem 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 \textsf{ibmq\_burlington} device with its coupling map shown in \cref{fig:coupling_map}, we can choose qubit 1 as the \bem ancilla qubit, which results in a \bem 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.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.2\textwidth]{coupling_map.jpg}
    \caption{Coupling map of quantum computing backend \textsf{ibmq\_burlington} provided on IBM Q website \cite{IBMQ}. The encircled number stands for the label of qubit where the color represents the noise rate of one-qubit gates on it. If two qubits are related by an arrow, the CNOT gate is directly available between these two qubits.\label{fig:coupling_map}}
\end{figure}

\begin{figure*}[htbp]
    \centering
    \includegraphics[width=\textwidth]{racbem_qc_nsysq_3.pdf}
    \caption{A \bem circuit constructed using the basic gate set $\{ \mathrm{U}_1, \mathrm{U}_2, \mathrm{U}_3, \mathrm{CNOT} \}$. The CNOT gates are directly implementable according to the coupling map in (a). $q_0, q_1, q_2, q_3$ refer to qubits $1,2,3,4$ in (a), where qubit $0$ is excluded as a signal qubit. The circuit at the bottom is a continuation of the top circuit.\label{fig:racbem_3qubit}}
\end{figure*}

\vspace{1em}
\noindent\textit{\hbem:}
In many applications, such as Hamiltonian simulation, thermal state preparation etc, we are only interested in Hermitian matrices.  It is possible to find a general circuit $U_A$ that coincidentally block-encodes a Hermitian matrix, but this can become increasingly difficult as $n$ increases. A useful fact is that once a random circuit $U_A$ is given, its Hermitian conjugate $U_A^{\dag}$ is easily accessible by conjugating the each gate and reversing the gate sequence. We will show below that this allows us to get access to in principle any $n$-qubit Hermitian matrix.

\begin{figure*}[htbp]
\begin{center}
\[
\Qcircuit @C=0.8em @R=1.em {
 \lstick{\ket{0}}& \gate{\mathrm{H}} & \targ & \gate{e^{-\I \varphi_0 \mathrm{Z}}} & \targ & \qw & \targ & \gate{e^{-\I \varphi_{1} \mathrm{Z}}} & \targ & \qw & \targ & \gate{e^{-\I \varphi_0 \mathrm{Z}}} & \targ & \gate{\mathrm{H}}&\qw  \\
\lstick{\ket{0}}& \qw  &\ctrlo{-1} & \qw  & \ctrlo{-1} & \multigate{1}{U_A} & \ctrlo{-1} & \qw & \ctrlo{-1} & \multigate{1}{U^{\dag}_A}    &\ctrlo{-1} & \qw & \ctrlo{-1} & \qw &\qw\\
\lstick{\ket{\psi}}& \qw &\qw &\qw&\qw &\ghost{U_A} &\qw&\qw&\qw&\ghost{U^{\dag}_A}&\qw&\qw   & \qw&\qw& \qw
}
\]
\end{center}
\caption{Quantum circuit for generating a $(1,2,0)$-block-encoding of a \hbem from a $(1,1,0)$-block-encoding $U_A$ and its Hermitian conjugate. $\mathrm{H}$ is the Hadamard gate, and $\mathrm{Z}$ is the Pauli-Z gate.}
\label{fig:qsp_circuit_hermitian_racbem}
\end{figure*}

Consider the quantum circuit in \cref{fig:qsp_circuit_hermitian_racbem} denoted by $U_{\mathfrak{H}}$, where $\varphi_0,\varphi_1\in [-\pi,\pi)$. Direct calculation (see \cref{sec:hracbem}) shows that
\begin{equation}
\begin{split}
\mathfrak{H}=&(\bra{0^2}\otimes I_n)U_{\mathfrak{H}}(\ket{0^2}\otimes I_n)\\
=&\left[-2 \sin (2\varphi_0) \sin \varphi_1\right] A^{\dag}A+\cos(2\varphi_0-\varphi_1).
\end{split}
\label{eqn:hracbem}
\end{equation}
Here $\mathfrak{H}$ is a Hermitian matrix. We refer to it as a Hermitian \bem (\hbem), 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^{\dag} A$ is Hermitian positive semi-definite. This will be referred to as a \chbem.
In other words, a \chbem is constructed from its (non-unique) matrix square root $A$.

In \cref{fig:qsp_circuit_hermitian_racbem}, 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 \chbem, we only need 1 application of $U_A$, 1 application of $U_A^{\dag}$, 2 H gates, 2 standard CNOT gates, 1 $\mathrm{S}^{\dag}$ gate, and 2 T gates (see \cref{sec:hracbem}). Since any matrix with singular values bounded by $1$ can be represented as a \bem, we immediately have that any Hermitian positive semi-definite matrix with eigenvalues bounded by $1$ can be represented as a \chbem, with a sufficiently flexible $U_A$.

\vspace{1em}
\noindent\textit{Quantum singular value transformation:}
Given the SVD $A=W\Sigma V^{\dag}$, and a smooth function $f(x)$ of even parity, we define the quantum singular value transformation (QSVT) as
\begin{equation}
f^{\triangleright}(A):=V f(\Sigma)V^{\dag}.
\label{eqn:right_generalize}
\end{equation}
Here the right pointing triangle reflects that only the right singular vectors $V$ are kept. Clearly $f^{\triangleright}(A)=f(\sqrt{A^{\dag}A})$, where the right hand side is the standard matrix function.  Now let $f$ be a real even polynomial of degree $2d$ that satisfies $\abs{f(x)}\le 1$ for any $x\in[-1,1]$. Let $U_A$ be a $(1,m,0)$-block-encoding of $A$. Then following \cite[Corollary 11]{GilyenSuLowEtAl2019}, there exists a $(1, m+1,0)$-block-encoding of $f^{\triangleright}(A)$, denoted by $U_{f^{\triangleright}(A)}$.  The circuit to implement $U_{f^{\triangleright}(A)}$ is given in \cref{fig:qsp_circuit_real}, which can be constructed  using $d$ queries of $U_A$ and $d$ queries of $U_A^{\dag}$.
Here $\Phi:=(\varphi_0,\ldots,\varphi_{2d})$ are called the phase factors. One challenge in QSVT is to find the phase factors $\Phi$. Besides the methods for obtaining $\Phi$ by polynomial factorization \cite{GilyenSuLowEtAl2019,Haah2019}, recently an optimization based method is proposed to find $\Phi$ up to very high degrees \cite{DongMengWhaleyEtAl}. A brief description of the method is given in \cref{sec:optimization}.

\begin{figure*}[htbp]
\begin{center}
\[\scalebox{0.9}{
\Qcircuit @C=0.8em @R=1.em {
 \lstick{\ket{0}}& \gate{\mathrm{H}} & \targ & \gate{e^{-\I \varphi_{2d} \mathrm{Z}}} & \targ & \qw & \targ & \gate{e^{-\I \varphi_{2d-1} \mathrm{Z}}} & \targ & \qw & \qw &\raisebox{0em}{$\cdots$}&&\qw   &\targ & \gate{e^{-\I \varphi_0 \mathrm{Z}}} & \targ & \qw& \gate{\mathrm{H}}&\qw  \\
\lstick{\ket{0^m}}& \qw  &\ctrlo{-1} & \qw  & \ctrlo{-1} & \multigate{1}{U_A} & \ctrlo{-1} & \qw & \ctrlo{-1} & \multigate{1}{U^{\dag}_A} &\qw &\raisebox{0em}{$\cdots$} &&\qw    &\ctrlo{-1} & \qw & \ctrlo{-1} & \qw& \qw &\qw \\
\lstick{\ket{\psi}}& \qw &\qw &\qw&\qw &\ghost{U_A} &\qw&\qw&\qw&\ghost{U_A}&\qw &\raisebox{0em}{$\cdots$} &&\qw&\qw&\qw   & \qw&\qw& \qw &\qw
}}
\]
\end{center}
\caption{Quantum circuit for quantum singular value transformation of a real matrix polynomial of degree $2d$. For QSVT of a \hbem, the number of ancilla qubits $m$ can be set to $1$.}
\label{fig:qsp_circuit_real}
\end{figure*}

Therefore \cref{fig:qsp_circuit_hermitian_racbem} 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
\begin{equation}
h(x)=\left[-2 \sin (2\varphi_0) \sin \varphi_1\right] x^2+\cos(2\varphi_0-\varphi_1).
\label{eqn:h_polynomial}
\end{equation}
A \chbem is given by $h(x)=x^2$.

Consider any real polynomial $g(x)$ of degree $d$ without parity constraint, satisfying $\abs{g(x)}\le 1$ for any $x\in[-1,1]$.  Then using the identity
\begin{equation}
g(\mathfrak{H})=(g\circ h)^{\triangleright}(A):=f^{\triangleright}(A),
\label{eqn:poly_hracbem}
\end{equation}
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$ \cite{GilyenSuLowEtAl2019}. Normally such controlled operations have significant overhead. For instance, if we would like to implement a controlled-\bem, 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 \hbem $\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 \hbem). Furthermore, each controlled gate in \cref{fig:qsp_circuit_real} is a standard CNOT gate rather than a Toffoli gate.  Therefore the \hbem 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 \hbem is also helpful: Note that if we remove all Z-rotations in the first line of \cref{fig:qsp_circuit_real}, the circuit implements an identity operator since $U_A,U_A^{\dag}$ 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 \cite{GilyenSuLowEtAl2019}.

\section{Proposal of quantum LINPACK benchmark}\label{sec:linpack}

The \bem, as well as the \hbem 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 \hbem model for solving QLSP. We assume $A=\mathfrak{H}$ is a \hbem, 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 $\Or(\kappa^{1-\delta})$ for any $\delta>0$ \cite{HarrowHassidimLloyd2009}, and the cost of using QSVT to solve general linear systems scales as $\wt{\Or}(\kappa^2\log(1/\epsilon))$ \cite{GilyenSuLowEtAl2019}.
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 \chbem $\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 \hbem model offers a new and natural way to solve this problem. From \cref{eqn:hracbem}, assuming $\cos(2\varphi_0-\varphi_1)>0$ and  $-2\sin(2\varphi_0)\sin\varphi_1>0$, and use that $0\preceq A^{\dag}A\preceq 1$, the condition number of $\mathfrak{H}$ can be bounded from above:
\[
\kappa(\mathfrak{H})\le \frac{\cos(2\varphi_0+\varphi_1)}{\cos(2\varphi_0-\varphi_1)}.
\]
Therefore the condition number of a \hbem is fully tunable by changing the phase factors $\varphi_0,\varphi_1$. According to \cref{eqn:h_polynomial}, 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
\[
|g(x)-x^{-1}|\le \epsilon, \quad x\in[\kappa^{-1}, 1].
\]
Following \cite[Corollary 69]{,GilyenSuLowEtAl2018}, there can be satisfied by an odd polynomial $g(x)$ with degree $d\sim\Or(\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 \cite{DongMengWhaleyEtAl}, and the polynomial can be chosen to be either even or odd.  \cref{fig:xinv} 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.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=5cm]{xinv-approximation.pdf}
    \caption{Approximate polynomials obtained by the Remez algorithm. The even and odd polynomials are generated such that the maximal approximation error to $1/(\kappa x)$ on $[\kappa^{-1}, 1]$ is $10^{-3}$. We find that the polynomial degree of the even approximation is lower than that of the odd approximation.}
    \label{fig:xinv}
\end{figure}

We may then construct a degree $2d$ polynomial $f=g\circ h$ as in \cref{sec:racbem}. Once we find the associated phase factors, the circuit in \cref{fig:qsp_circuit_real} implements $g(\mathfrak{H})\ket{b}$, which satisfies the error bound $\norm{g(\mathfrak{H})\ket{b}-\kappa^{-1} \mathfrak{H}^{-1}\ket{b}}_2\le \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
\[
p=\norm{g(\mathfrak{H})\ket{b}}_2^2\ge (\kappa^{-1}-\epsilon)^2,
\]
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 \emph{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 \hbem are given, the quantum LINPACK benchmark reports the relative error $\abs{p-p_{\mathrm{exact}}}/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 \emph{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 \emph{minimal requirement} in order to achieve quantum advantages by solving QLSP or similar linear algebra tasks.

\section{Other quantum linear algebra applications}

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 \hbem defined via a unitary matrix $U_A$ and $h(x)$ is a second order polynomial.

\vspace{1em}
\noindent\textit{Time series analysis:} Given an $n$-qubit state $\ket{\psi}$, consider the computation of the following time series
\begin{equation}
s(t)=\braket{\psi|e^{\I \mathfrak{H}t}|\psi}.
\label{eqn:signal_proc}
\end{equation}
When $\mathfrak{H}$ is a \hbem, we can evaluate $s(t)$ by measuring the its real and imaginary components separately:
\begin{equation}
\begin{aligned}
s(t)=&\braket{\psi|\cos (\mathfrak{H}t)|\psi}+\I\braket{\psi|\sin (\mathfrak{H}t)|\psi}\\
=&2\Braket{\psi|\frac12(\cos(\mathfrak{H}t)+I_n)|\psi}\\
&+2\I\Braket{\psi|\frac12(\sin(\mathfrak{H}t)+I_{n})|\psi}-(1+\I)\\
=&:2\norm{g_{c,t}(\mathfrak{H})\ket{\psi}}_2^2+2\I\norm{g_{s,t}(\mathfrak{H})\ket{\psi}}_2^2-(1+\I).
\end{aligned}
\label{eqn:signal_proc_identity}
\end{equation}
Here we introduce the functions for the cosine and sine part, respectively
\[
g_{c,t}(x)=\sqrt{\frac{\cos(xt)+1}{2}}, \quad
g_{s,t}(x)=\sqrt{\frac{\sin(xt)+1}{2}}.
\]
Now the quantities in \cref{eqn:signal_proc_identity} can be directly obtained via the success probability of the QSVT circuit with $f_{c,t}:=g_{c,t}\circ h$ and $f_{s,t}:=g_{s,t}\circ h$, respectively (with a suitable polynomial approximation of $g_{c,t},g_{s,t}$). Note that the access to the matrix square root of $\mathfrak{H}$ here is crucial: though $g_{c,t}$ is an even function itself, $g_{s,t}$ does not have a definite parity, and therefore the direct implementation of $g_{s,t}(\mathfrak{H})$ would require a LCU circuit and hence controlled unitary operations. In the setting of \hbem, the treatment of $g_{c,t}$ and $g_{s,t}$ can be put on the same footing due to the composition with the even polynomial $h$.

In order to accelerate the convergence of the polynomial approximation, in practice, we may introduce another tunable parameter $\eta \ge 1$ in the formulation so that
\[
g_{c,t,\eta}(x) = \sqrt{\frac{\cos(xt)+\eta}{2}}, \quad
g_{s,t,\eta}(x) = \sqrt{\frac{\sin(xt)+\eta}{2}}
\]
and then
\[
s(t) = 2\norm{g_{c,t,\eta}(\mathfrak{H})\ket{\psi}}_2^2 + 2\I\norm{g_{s,t,\eta}(\mathfrak{H})\ket{\psi}}_2^2 - \eta(1+\I).
\]
For instance, $\eta$ can be set to $1.0$ or $1.5$ and the details are given in \cref{tab:param-ts}.

In \cite{Somma2019}, 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.

\vspace{1em}
\noindent\textit{Spectral measure:} Given an $n$-qubit state $\ket{\psi}$, in order to approximately evaluate the spectral measure, we may use the Plemelj formula
\begin{equation}
\begin{split}
s(E)=&\braket{\psi|\delta(\mathfrak{H}-E)|\psi}\\
=&\lim_{\eta\to 0^+} \frac{1}{\pi} \Im \braket{\psi|(\mathfrak{H}-E-\I \eta)^{-1}|\psi}, \quad E\in[-1, 1].
\end{split}
\label{eqn:spec_measure}
\end{equation}
By choosing a finite $\eta$ as the broadening parameter, we need to evaluate \[
s_{\eta}(E)=\frac{\eta}{\pi}\Braket{\psi|((\mathfrak{H}-E)^2+\eta^2)^{-1}|\psi}.
\]
Now define
\[
g_{\eta,E}(x)=\left(\frac{\eta^2}{(x-E)^2+\eta^2}\right)^{\frac12},
\]
which satisfies $|g_{\eta,E}|\le 1$ for $x\in[-1,1]$, then $g_{\eta,E}\circ h$ can be approximated by an even order polynomial. Therefore
\[
s_{\eta}(E)=\frac{1}{\eta\pi}\norm{g_{\eta,E}(\mathfrak{H})\ket{\psi}}^2_2,
\]
and it can be obtained via the success probability of the QSVT circuit with $f=g_{\eta,E}\circ h$ (with a suitable polynomial approximation of $g_{\eta,E}$).

\vspace{1em}
\noindent\textit{Thermal averages:} In order to evaluate the thermal average of the energy ($\beta$ is the inverse temperature)
\[
E(\beta)=\frac{\Tr[\mathfrak{H} e^{-\beta \mathfrak{H}}]}{\Tr[e^{-\beta \mathfrak{H}}]},
\]
we may use the  minimally entangled typical thermal state (METTS) algorithm \cite{White2009,StoudenmireWhite2010,MottaSunTanEtAl2020}. Let $Z=\Tr[e^{-\beta \mathfrak{H}}]$, and
\[
E(\beta)=\frac{1}{Z}\sum_{i\in[N]}\braket{i|e^{-\beta \mathfrak{H}/2} \mathfrak{H} e^{-\beta \mathfrak{H}/2}|i}=\sum_{i\in[N]} p_{i}\braket{\phi_i|\mathfrak{H}|\phi_i}.
\]
Here $i$ loops over all states of the computational basis, each represented by an $n$-bit string. From the unnormalized states $\ket{\wt{\phi}_i}= e^{-\beta \mathfrak{H}/2}\ket{i}$, we define a probability distribution  $p_i=\braket{\wt{\phi_i}|\wt{\phi_i}}/Z$, and corresponding normalized states
$\ket{\phi_i}=\ket{\wt{\phi}_i}/\sqrt{\braket{\wt{\phi_i}|\wt{\phi_i}}}$. For each $i$, the contribution to the energy can be expressed as
\begin{equation}
\braket{\phi_i|\mathfrak{H}|\phi_i}=\frac{\braket{i|e^{-\beta \mathfrak{H}} \mathfrak{H} |i}}{\braket{i| e^{-\beta \mathfrak{H}}|i}}
=\frac{\norm{g_{n,\beta}(\mathfrak{H})\ket{i}}_2^2}{\norm{g_{d,\beta}(\mathfrak{H})\ket{i}}_2^2}.
\label{eqn:thermal_eachterm}
\end{equation}
Here we define
\[
g_{n,\beta}(x)=e^{-\beta x/2} \sqrt{x}, \quad g_{d,\beta}(x)=e^{-\beta x/2},
\]
for the numerator and the denominator, respectively, and without loss of generality we assume $\mathfrak{H}\succeq 0$. The numerator and the denominator can be obtained via the success probability of the QSVT circuit with $f_{n,\beta}:=g_{n,\beta}\circ h$ and $f_{d,\beta}:=g_{d,\beta}\circ h$, respectively (with a suitable polynomial approximation of $g_{n,\beta},g_{d,\beta}$).

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
\begin{enumerate}
\item Evaluate the contribution to the energy from the state $\ket{\phi_i}$ via \cref{eqn:thermal_eachterm}.

\item Collapse $\ket{\phi_i}$ to a new state in the computational basis $\ket{i'}$ with probability $|\braket{i'|\phi_i}|^2$, and repeat step 1.
\end{enumerate}

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{\wt{\phi}_i}$, which is readily available when computing the denominator in \cref{eqn:thermal_eachterm}. Then collapsing to $\ket{i'}$ 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 \chbem with $h(x) = x^2$, an alternative method is to approximate an odd function $\wt{f}_{n,\beta}(x) = x e^{-\beta x^2 / 2}$ instead of $f_{n,\beta}$. Then the complexity for thermal average calculation can be only $\wt{\Or}(\sqrt{\beta}\log(1/\epsilon))$ \cite[Corollary 64]{GilyenSuLowEtAl2018}.

\section{Numerical results}
The source code of \bem is available in the Github repository\footnote{\url{https://github.com/qsppack/RACBEM}}. We use the optimization-based algorithm which is described in \cref{sec:optimization} to generate phase factors and the source code is available in the Github repository\footnote{\url{https://github.com/qsppack/QSPPACK}}. We demonstrate the numerical performance of the \bem  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 \textsf{python3.7} and \textsf{Qiskit} \cite{Qiskit}. 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 \cref{app:model} 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 \cref{alg:custom-rqc} (in \cref{app:model}) 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 \ge 3$. We report the statistics of singular-value distributions of the block-encoded matrix in \cref{fig:svd-dist} (\cref{app:details}) 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 \cref{app:model} 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 \bem whose depth is $\ell$ and with $n$ system qubits, to implement the QSVT of a real polynomial of degree $\qspdeg$, the total gate count for the circuit is upper bounded by $2 + 7(\qspdeg+1) + \qspdeg  \ell n$.  The details about QSVT phase factors used in numerical experiments can be found in \cref{app:details}. Unless otherwise noted, the input quantum state of the quantum circuit is set to $\ket{0^{n+2}}$.

\vspace{1em}
\noindent\textit{\bem:}
Before presenting results of various numerical linear algebra problems, we measure the effect of noise on \bem directly on quantum computing backends provided by IBM Q. The numerical results are displayed in \cref{fig:be-ibmq-benchmark}. Each data point is obtained by generating a \bem using \cref{alg:custom-rqc}, and we measure the success probability of obtaining the block-encoded matrix applied to $\ket{0^n}$, which is equal to $\norm{A\ket{0^n}}_2^2$. The number of repeated measurements (shots) is 8192. \cref{fig:be-ibmq-benchmark} shows that the relative error of the quantum device can be considerable, ranging  from $10\%$ to $30\%$ as $n$ increases. Since the \bem is the building block of all subsequent quantum linear algebra tasks, we expect that the relative error of such tasks on quantum computing backends provided by IBM Q should be at least around the same level.

\begin{figure}[htbp]
    \centering
    \includegraphics[width = 8cm]{be_ibmq_benchmark.pdf}
    \caption{Relative error in success probability of obtaining $A\ket{0^n}$ for different number of system qubits.  The box ranges from the first to the third quartiles (25\% and 75\% percentiles respectively) of the dataset. The horizontal line inside the box stands for the median of the dataset. The whisker extending from the box indicates the rest of the dataset. The \bem's are generated using our custom random circuit generator in which the coupling map, basic gates and logical-to-physical layout are specified. When the number of system qubits is less than or equal to $3$, the results are obtained on the $5$-qubit backend \textsf{ibmq\_burlington}. Other results are obtained using $15$-qubit backend \textsf{ibmq\_16\_melbourne}. The dataset consists of results from $200$ \bem's when $n=1,2,3$, $300$ \bem's when $n=4,5,6$, $400$ \bem's when $n=7$ and $500$ \bem's when $n=8$.}
\label{fig:be-ibmq-benchmark}
\end{figure}

\vspace{1em}
\noindent\textit{QLSP:} According to the discussion in \cref{sec:linpack}, it is possible to measure the accuracy of the solution to QLSP by sampling the output distribution of $g(\mathfrak{H})\ket{0^n}$ using e.g. a cross-entropy test similar to that in \cite{AruteAryaBabbushEtAl2019}. Here for simplicity we only measure $\norm{g(\mathfrak{H})\ket{0^n}}_2^2$ as a success probability (of measuring all ancilla qubits and obtain 0's), and evaluate the relative error compared to the exact value $\norm{\scalefac^{-1} \mathfrak{H}^{-1}\ket{0^n}}_2^2$ evaluated on a classical computer where $\scalefac$ is a scaling factor (see \cref{tab:param-qlsp} for details). The small relative error in success probability is a necessary condition to ensure the solution is correct which yields the benchmark in a weak sense. The condition number of \hbem is controlled by a second order polynomial $h_\kappa(x) := (1-\kappa^{-1})x^2 + \kappa^{-1}$ as in \cref{sec:linpack}. Another polynomial $g_\kappa(x)$ which approximates $x^{-1}$ is chosen to perform matrix inversion. The composite polynomial $g_\kappa \circ h_\kappa$ can be implemented by an even order QSVT circuit to carry out the overall procedure with only two ancilla qubits.

We report the performance of solving linear systems on the IBM Q device using the \hbem model in \cref{fig:ls-ibmq-benchmark}. The architecture of the five backends are identical, and therefore we may draw a \hbem 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 \cref{fig:ls-ibmq-benchmark} indicate that for the shallow circuit,  the relative error of the success probability is similar to that observed in \cref{fig:be-ibmq-benchmark}, 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.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=8cm]{ls_ibmq_benchmark.pdf}
    \caption{Relative error of solving QLSP using the \hbem model on IBM Q 5-qubit backends with identical quantum architectures. The condition number of each \hbem is upper bounded by $2$ and the number of system qubits is $3$. The size of the box indicates the first and the third quartiles and the extending whisker indicates the rest of the dataset. The phase factor sequences of length $3$ and $11$ are used to carry out QSVT to solve QLSP. Each dataset is obtained using $100$ samples at random. The details about phase factors are given in \cref{tab:param-qlsp}.}
    \label{fig:ls-ibmq-benchmark}
\end{figure}

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 \cref{fig:ls-qasm-kappa}. 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 $\Or(\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.

\begin{figure}[htbp]
    \centering
    \includegraphics[width = 8cm]{ls_qasm_kappa.pdf}
    \caption{Relative error of solving QLSP for different number of system qubits, condition number and noise magnitude. The numerical results are obtained by running on QVM with architecture and noise model retrieved from IBM Q backend \textsf{ibmq\_16\_melbourne}. Each kind of parameters is computed by solving QLSP of $300$ \bem's at random. The marker is the average of each dataset and the shaded area is the 95\% confidence interval inferred from the dataset. The details about phase factors are given in \cref{tab:param-qlsp}.}
    \label{fig:ls-qasm-kappa}
\end{figure}

\vspace{1em}
\noindent\textit{Time series analysis:}
The numerical results regarding the time series analysis are shown in \cref{fig:ts}. The results in \cref{fig:ts-ibmq} are obtained on the IBM Q backends. When the number of system qubit is $1$ (the circuit uses 3 qubits in total), the features of the time series obtained from the quantum device can agree qualitatively with the exact solution. However, as the number of system qubits increases to $5$, the result is dominated by the quantum noise. In order to investigate the performance of larger systems and the effects of noise, we use the tunable QVM instead in \cref{fig:ts-qasm}.  The length of phase factors is chosen adaptively as $t$ increase in order to reduce the error of the polynomial approximation (details in \cref{tab:param-ts}). When the noise level $\sigma$ is tuned to $0$, the results from the QSVT circuit is uniformly accurate for all $t$ in $[1,10]$. However, since the circuit depth needs to increase proportionally to $t$, when the noise level is not $0$, the error is significantly larger as $t$ increases. We also remark that our custom noise model with $\sigma=1$ only contains partial kinds of noise derived from the real IBM Q device (details in \cref{app:model}). Therefore, in \cref{fig:ts}, the time series computed on the real device behaves much worse than that computed from the simulation using the noise model with $\sigma=1$ on the QVM.

\begin{figure*}[htbp]
    \centering
    \vspace{-10pt}
    \subfigure[\label{fig:ts-ibmq}]{
        \includegraphics[width=11cm]{sp_ibmq.pdf}
    }

    \subfigure[\label{fig:ts-qasm}]{
        \includegraphics[width=14cm]{sp_qasm.pdf}
    }
    \vspace{-10pt}
    \caption{Time series computed by using the \chbem model and QSVT. The details about phase factors are given in \cref{tab:param-ts} and the details about block-encoded matrices are given in \cref{fig:svd-matrix}. (a) Results obtained on IBM Q backends. The probability of CNOT is 0.1 when generating random circuits. (b) Results obtained on QVM with architecture and noise model retrieved from the IBM Q backend \textsf{ibmq\_16\_melbourne}. The probability of CNOT is set to 0.5.}
    \label{fig:ts}
\end{figure*}

\vspace{1em}
\noindent\textit{Spectral measure:} In the limit when $\eta\to 0$, the spectral measure is given by the summation of Dirac-$\delta$ functions. Hence when $n$ is small, the spectral measure has sharp features even when $\eta$ takes a finite value, which in turn requires a polynomial of higher degree to resolve. So we only consider  spectral measures for $n=7,8,9,10$ on QVM, and the size of $\mathfrak{H}$ ranges from $128$ to $1024$. The numerical results of spectral measures are presented in \cref{fig:sm-qasm}. The spectral measures exhibit rather different features. This is not so much related to the number of system qubits, and is mostly due to the specific instance of the \hbem. In all cases, the quantum algorithm can capture at least the qualitative features of the spectral measures, though the noise plays an important role particularly for the instance $n=9$. We also remark that the functionality of the QSVT circuit depends only on the set of phase factors as $E$ sweeps across the spectrum.  The length of corresponding QSVT phase factors is set to $11$ for each point. In \cref{tab:param-sm} and \cref{tab:param-sm-large}, for points closer to the middle of the spectrum ($E=0.5$), we observe that a polynomial of larger degree is needed to achieve the same accuracy in approximation. Therefore, it can be seen from \cref{fig:sm-qasm} that the deviation between the value given by noiseless QSVT and the exact one increases as the parameter moves towards the middle of the spectrum. The choice of the circuit depth reflects the tradeoff between approximation error and the effects of noise. To illustrate this, we also compute the spectral measures in \cref{app:add-num} by using a deeper QSVT circuits. Although the deeper circuit can produce better results in the noiseless setting, the quality of the spectral measure can  significantly deteriorate in the presence of the quantum noise.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=8cm]{sd_qasm.pdf}
    \caption{Spectral measure computed by using the \chbem model and QSVT. The numerical results are obtained by running on QVM with architecture and noise model retrieved from IBM Q backend \textsf{ibmq\_16\_melbourne}. For each number of system qubits, we draw a \chbem randomly according to the backend architecture. The probability of CNOT is 0.5 when generating random circuits. The details about phase factors are given in \cref{tab:param-sm} and the details about block-encoded matrices are given in \cref{fig:svd-matrix}.}
    \label{fig:sm-qasm}
\end{figure}

\vspace{1em}
\noindent\textit{Thermal average of the energy:}
In \cref{fig:ta-ibmq}, we compute the thermal average of the energy of \hbem's on IBM Q backends for $n$ from $1$ to $5$. We reduce the length of QSVT phase factors as much as possible while balancing the approximation error and the quantum error, and the details can be found in \cref{tab:param-ta}. Since METTS is a Monte Carlo algorithm, we perform a sufficiently large number of steps to reduce the error due to the METTS algorithm, by monitoring the cumulative moving average of the thermal average of the energy in \cref{appfig:ta-cma} in \cref{app:add-num}.
Compared to the results above obtained on IBM Q, the results for the thermal average of the energy are somewhat surprisingly accurate in all cases. We also compute the thermal average of the energy on QVM to further investigate the effect of the noise for $n=3, 5$ and $7$. In \cref{fig:ta-qasm}(i) and (ii), the identical quantum computing task is emulated on different backends, and the only difference is due to the noise model. It is evident that the behaviour of the solution depends sensitive on the noise. On the other hand, we find that the solution remains remarkably accurate as $n$ increases, even when $\beta$ is relatively large. For most data points, the thermal average of the energy decreases monotonically with respect to $\beta$, and we observe that as $\sigma$ increases, the energy curve shifts downwards monotonically.

 \begin{figure}[htbp]
    \centering
    \vspace{5pt}
    \subfigure[\label{fig:ta-ibmq}]{
        \includegraphics[width=.45\textwidth]{ta_ibmq.pdf}
    }
    \subfigure[\label{fig:ta-qasm}]{
        \includegraphics[width=.45\textwidth]{ta_qasm.pdf}
    }
    \vspace{-10pt}
     \caption{Thermal average of the energy computed by using the \chbem model and QSVT. The details about phase factors are given in \cref{tab:param-ta} and the details about block-encoded matrices are given in \cref{fig:svd-matrix}. (a) Results obtained on IBM Q backends. The probability of CNOT is 0.1 when generating random circuits. (b) Results obtained on QVM. The probability of CNOT is set to 0.5.}
     \label{fig:ta}
 \end{figure}

\section{Conclusion}
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 (\hbem) input model and the quantum singular value transformation (QSVT). The flexibility provided by the \hbem 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 \cite{EmersonWeinsteinSaracenoEtAl2003,HarrowLow2009,NahumRuhmanVijayEtAl2017,AaronsonArkhipov2011,BremnerMontanaroShepherd2016,Fujii2016,AaronsonChen2016,BoulandFeffermanNirkheEtAl2019,AaronsonGunn2019}. 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 \emph{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 \emph{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.

\section*{Acknowledgements}
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.

\bibliographystyle{abbrvnat}

config: arxiv


TITLE: Solving second order differential equation involving sine around function
QUESTION [0 upvotes]: I didn't really know how to put this in the title as I am not too familiar with solving differential equations. I am trying to solve:
$x''(t) = -Asin(x(t))$
I'm not sure how to even begin solving this equation and have been searching online and can't find anything (although I may be searching the wrong thing).
Many thanks.

REPLY [1 votes]: It's easy enough to integrate once: multiply both sides by $2x'(t)$ and integrating gives
$$ x'(t)^2 = 2A(\cos{x(t)}-\cos{x_0}). $$
Taking the square root and dividing gives
$$ 1 = \frac{x'(t)}{\sqrt{2A(\cos{x}-\cos{x_0})}}. $$
This doesn't have an elementary integral: one can express $t$ in terms of the elliptic integral of the first kind,
$$ \pm \sqrt{A}(t-t_0) = (\csc{(x_0/2)})F( x/2 , \csc{(x_0/2)} ) $$
This has an inverse, the Jacobi amplitude $\operatorname{am}$. Specifically,
$$ x = \pm 2\operatorname{am}(\sqrt{A}(t-t_0)\sin{(x_0/2)},\csc{(x_0/2)}). $$


TITLE: Lagrangian Oscilattor
QUESTION [1 upvotes]: I want to know how to calculate the normal modes from a Lagrangian.
I make the T (kinetic energy matrix) and U (potential energy matrix), and then I calculate the determinant of $|T-\omega ^2 U|$ (with omega as the angular speed). This way I find the $\omega$'s but I dont know how to find the normal modes.
Are they just the eigenvectors of the $T-\omega ^2 U$?

REPLY [2 votes]: $\newcommand{\oh}[0]{\frac{1}{2}}$You're basically right. It is the eigenvector times $\cos(\omega t)$. I will put in more details,plaigarizing from a pervious answer of mine.
I will adapt the discussion in landau mechanics 3rd edition starting at the bottom of page 65. In general the hamiltonian is $L = \oh \dot{x}_i m_{ij} \dot{x}_j - \oh x_i V_{ij} x_j$. The equation of motion (EOM) is $m_{ij} \ddot{x}_j = - V_{ij} x_j$. We seek a solution where each coordinate oscillates with the same frequency: $x_j(t) = A_j \sin(\omega t + \phi_j)$. Plugging this into the EOM we get $\omega^2 m_{ij} A_j - V_{ij} A_j = 0$. i.e. $A_j$ is a null vector of $\omega^2 m_{ij} - V_{ij}$. This is only possible if $\omega$ takes on such a value where the determinant of $\omega^2 m_{ij}  - V_{ij}$ is zero. The number of $\omega^2$'s which give zero determinant will be equal to the dimensionality of our system. We will index this set of $\omega$'s by the greek letter $\alpha$.
Now say we have found a $\omega_\alpha$ where the determinant of $\omega_\alpha^2 m_{ij}  - V_{ij}$ is zero. Then we must find the $A_{k\alpha}$ which is a null vector of $\omega_\alpha^2 m_{ij}  - V_{ij}$. Having done this for each $\alpha$, we find that the most general solution is $x_k(t) = A_{k \alpha} q_\alpha \sin(\omega_\alpha t - \phi_\alpha)$. The $q_\alpha$ are seen to be the normal mode coordinates.


TITLE: Shorter notation for $\{1,...,N\}$ and $\{-N,...,N\}$
QUESTION [0 upvotes]: Often I need to write something like $n \in \{1,...,N\}$ or $n \in \{-N,...,N\}$ in subscripts and sums. I'd like to use something a bit shorter. Looking around at all the notation questions on Math.SE I came up with the following:
$$ n \in \{1,...,N\} \rightarrow n \in \mathbb{Z}^+_{\le N}$$
$$ n \in \{-N,...,N\} \rightarrow |n| \in \mathbb{Z}_{\le N}$$
Is there any better notation I could use that is also short? Or should I just define my own notation like $\{1,...,N\} = \mathbb{Z}^N, \{-N,...,N\} = \mathbb{Z}^{|N|}$

REPLY [2 votes]: I've seen $[n]$ to mean $\{1,\dots,n\}$. This is used, for example, in the wonderful book generatingfunctionology. Don't assume your readers will be familiar with it, though; you should explain your notation in the text or a footnote.


TITLE: How can a physical quantity be neither vector nor scalar?
QUESTION [22 upvotes]: I was going through the definition of stress. I came to know that stress is a tensor quantity. When I checked the net about - 'What is a tensor quantity', I came to know that it is neither a scalar quantity nor a vector quantity. So how can a quantity be neither vector nor scalar?
Physicists on this site, please explain your answers in simple ways too, that would be just enough for others (who are school students, etc.) to understand. You can do so by keeping the heading "SIMPLER WAY" on top of your answers.

REPLY [0 votes]: I'm wondering just how you came across a stress tensor without meeting a linear algebra course first. Okay, never mind, I'll try to give you some hint of an answer, in accessible terms only.
Basically, a tensor is an even more intricate thingie than a vector.
(That's not entirely accurate, actually scalars and vectors are a special case of tensors, but you won't really hear people refer to them that way.)
Where a scalar has one component and a vector in your usual euclidean 3D space has three (x, y, z), most tensors you come across will have nine (even if probably not all of them independent). When writing these down, such a tensor will be represented by a 3x3 matrix of components:
xx   xy   xz
yx   yy   yz
zx   zy   zz
(Sorry, I don't know how to format this properly here.)
There are even more monstrous kinds of tensors than this, but you don't find a rank 3 or higher tensor all that often.


TITLE: Tossing a biased coin
QUESTION [1 upvotes]: We toss a biased coin $k$ times with the probability of tossing a head being $1/m$. We know that in the first two attempts there was at least one tail. What the probability that we tossed exactly $n$ heads?
Well, the tosses are independent so getting at least one tail in the first two attempts is the sum of $2\cdot (\tfrac{m-1}{m})^1 \cdot (\tfrac{1}{m})^2 = \tfrac{2m-2}{m^3}$ and $1\cdot (\tfrac{m-1}{m})^2 \cdot (\tfrac{1}{m})^0 = \tfrac{(m-1)^2}{m^2}$, which is $\tfrac{(m-1)(m^2-m+2)}{m^3}$. Thus getting at least on head in the first two attempts is $1-\tfrac{(m-1)(m^2-m+2)}{m^3}$.
Now we are left with $k-2$ tosses so we want to hit $n-1$ heads. The chance is $${k-2 \choose n-1}\tfrac{1}{m^{n-1}}\cdot (\tfrac{m-1}{m})^{k-2 - n +1}.$$
So the answer should be
$$1-\tfrac{(m-1)(m^2-m+2)}{m^3} + {k-2 \choose n-1}\tfrac{1}{m^{n-1}}\cdot (\tfrac{m-1}{m})^{k-1 - n}$$
Is this right?

REPLY [1 votes]: The probability that we tossed exactly n heads, given there was exactly 1 tail in the first two tries, is $Q_{1}={{k-2}\choose{n-1}}(\frac{m-1}{m})^{n-1}(\frac{1}{m})^{k-n-1}$, and the probability that we tossed exactly n heads, given there was exactly 2 tails in the first two tries, is $Q_{2}={{k-2}\choose{n}}(\frac{m-1}{m})^{n}(\frac{1}{m})^{k-n-2}$.
The probability that we had exactly $1$ tail in the first two tries, given that we know there was at least $1$ tail, is the conditional probability given by $P_{1}=\frac{P((exactly \ one \ tail)\cap(at \ least \ one \ tail))}{P(at \ least \ one \ tail)}=\frac{P(exactly \ one \ tail) }{P(at \ least \ one \ tail)}=\frac{2(\frac{1}{m})(\frac{m-1}{m})}{2(\frac{1}{m})(\frac{m-1}{m})+(\frac{m-1}{m})^{2}}$.
Similarly, the probability that we had exactly $2$ tails in the first two tries, given that we know there was at least $1$ tail, is given by $P_{2}=\frac{(\frac{m-1}{m})^{2}}{2(\frac{1}{m})(\frac{m-1}{m})+(\frac{m-1}{m})^{2}}$.
The total probability is $P_{1}Q_{1}+P_{2}Q_{2}$.


TITLE: When are $\frac{1}{|x|^s}$ and $\log|x|$ integrable near the origin?
QUESTION [6 upvotes]: When are $\frac{1}{|x|^s}$ for $s>0$ and $\log|x|$ integrable near the origin? I'm reading Evans PDE and in the construction of the fundamental solution of Poisson's equation, he defines
$$
\Phi(x) = C \log|x|
$$
for $\mathbb{R}^2$ and a suitable $C$, and
$$
\Phi(x) = C \frac{1}{|x|^{n-2}}
$$
for $\mathbb{R}^n$ with $n\geq 3$ and $C=C(n)$ a suitable constant.
The construction then goes about defining $\Phi * f$ for an $f \in C^2_c$. Clearly if $\Phi \in L^1_{\mathrm{loc}}$, then the convolution will be finite a.e., but I worry about  whether this holds because of the blow up near $0$. In general in what dimensions are
$$
\log |x|,~~~~~\frac{1}{|x|^s}
$$
integrable in a ball around $0$?

REPLY [4 votes]: For the logarithm, the answer is already implicitly given in the comments above. On $\mathbb{R}^{n}$ with $n \geq 2,$ we have, using polar co-ordinates, $$\int_{B_{1}(0)} \log |x| dx = \lim_{\varepsilon \rightarrow 0} \int_{B_{1}(0)\setminus B_{\varepsilon}(0)} \log |x| dx = \lim_{\varepsilon \rightarrow 0} \int_{B_{1}(0)\setminus B_{\varepsilon}(0)}  r^{n-1} \log r \ dr d\Omega_{S^{n-1}}, $$
where $d\Omega_{S^{n-1}}$ is the volume form on $S^{n-1}.$ Now, integrate by parts to explicitly integrate this and check this will always be finite as long as $n \geq 2.$
For $\frac{1}{|x|^{s}}$ would be integrable near the origin of $\mathbb{R}^{n}$ when $s < n. $ This can be easily checked using polar co-ordinates again.
$$ \int_{B_{1}(0)} \frac{1}{|x|^{s}} dx = \lim_{\varepsilon \rightarrow 0} \int_{B_{1}(0)\setminus B_{\varepsilon}(0)} \frac{1}{|x|^{s}} dx = \lim_{\varepsilon \rightarrow 0} \int_{B_{1}(0)\setminus B_{\varepsilon}(0)} \frac{1}{r^{s}} r^{n-1} dr d\Omega_{S^{n-1}}.$$
 If $s <n,$ integrating, we obtain,
$$ \int_{B_{1}(0)} \frac{1}{|x|^{s}} dx = \lim_{\varepsilon \rightarrow 0} \left[  \frac{r^{n-s}}{n-s}\right]_{\varepsilon}^{1} .\operatorname{Vol}(S^{n-1}).$$ This limit would be finite since $n-s >0$ and thus $\varepsilon^{n-s} \rightarrow 0$ as $\varepsilon \rightarrow 0.$ The integral will blow up for $s=n,$ since in that case after integrating we would obtain,
$$ \int_{B_{1}(0)} \frac{1}{|x|^{n}} dx = \lim_{\varepsilon \rightarrow 0} \left[  \log r \right]_{\varepsilon}^{1} .\operatorname{Vol}(S^{n-1}).$$ Now since $-\log \varepsilon \rightarrow + \infty$ as $\varepsilon \rightarrow 0,$ the inetgral would not be finite for $s=n$ ( and obviously also for $s > n$ as well ).


TITLE: Comparison Test for series.
QUESTION [0 upvotes]: I need to use the comparison test to determine whether the series:
$\sum\frac{3^n+7^n}{3^n+8^n}$
converges or diverges.
So far, I've tried splitting it into $\frac{3^n}{3^n + 8^n} + \frac{7^n}{3^n + 8^n} $ but have gotten nowhere

REPLY [0 votes]: It will be simpler to use  the limit comparison test via *asymptotice equivalence of functions:
$$3^n+7^n\sim_\infty 7^n,\quad3^n+8^n\sim_\infty 8^n,\quad\text{ hence }\quad \frac{3^n+7^n}{3^n+8^n}\sim_\infty \frac{7^n}{8^n}=\Bigl(\frac78\Bigr)^n,$$
which converges.


TITLE: Non finitely based varieties of groups defined by finitely many variables
QUESTION [3 upvotes]: A set $\Sigma$ of group identities is called bounded if there is $n\geq 1$ such that for any $(w\approx 1)\in \Sigma$, we have $w\in F(x_1, \ldots, x_n)$. A variety $\mathbf{V}$ is called bounded defined if $\mathbf{V}=Mod(\Sigma)$ for some bounded set $\Sigma$.
Question: Is there an example of a bounded defined non-finitely based variety of groups?
P.S. $F(x_1, \ldots, x_n)$ is the free group of rank $n$.

REPLY [5 votes]: Yes. It is usually called a non-finitely based variety of finite axiomatic rank.
Actually, one of the first examples of non-finitely based variety of groups has axiomatic rank two; Adian proved that the following set of identities is independent:
$$
\{[x^{pn},y^{pn}]^n=1\;|\; p \hbox{ is prime}\},
$$
where $n\ge 1003$ is any given odd number. Independence means that none of the laws is  a consequence of the other ones.


TITLE: Verify if a span of vectors is $\mathbb{R}^2$
QUESTION [1 upvotes]: if $H = \mathrm{Span}\{v_1, v_2\} $
and
$v_1 = (2,3)$
$v_2 = (-1,4)$
Is knowing $v_1$ and $v_2$ are linearly independent since $\det \begin{bmatrix} 2 & -1 \\ 3 & 4\end{bmatrix} = 11$ enough to say that $ H = \mathbb{R}^2$

REPLY [1 votes]: Also if don't want to use a theorem and can prove this directly like this:
Let $x \in \mathbb{R}^2$ such that $x=(x_1,x_2)$ and $x=av_1+bv_2$
$(x_1,x_2)=a(2,3)+b(-1,4) \Longleftrightarrow (x_1,x_2)=(2a-b,3a+4b)$
Thus  $$2a-b=x_1$$ $$3a+4b=x_2$$
So we have by solving this system: $$b=\frac{2x_2-3x_1}{11}$$ $$a=\frac{8x_1+2x_2}{22}$$
Thus for any $x=(x_1,x_2)$ you can find $a,b \in \mathbb{R}$ such that $x=av_1+bv_2$


TITLE: solve equation in positive integers
QUESTION [3 upvotes]: Can anybody help me with this equation?
Solve in $\mathbb{N}$:
$$3x^2 - 7y^2 +1=0$$
One solution is the pair $(3,2)$, and i think this is the only pair of positive integers that can be a solution. Any idea?

REPLY [1 votes]: For some Pell equations. For special cases, there is a formula describing their solutions.
If the equation: $aX^2-qY^2=f$ rational root $\sqrt{\frac{f}{a-q}}$
Solutions can be written using the following equation Pell:  $p^2-aqs^2=1$
Solutions have the form:
$Y=(2aps\pm(p^2+aqs^2))\sqrt{\frac{f}{a-q}}$
$X=(2qps\pm(p^2+aqs^2))\sqrt{\frac{f}{a-q}}$
According to these decisions can be found double this solution:
$Y_2=Y+2as(Yqs-Xp)$
$X_2=X+2p(Yqs-Xp)$
If the other root is rational: $\sqrt{\frac{f}{a}}$
Solution has the form:
$Y=2ps\sqrt{fa}$
$X=(p^2+aqs^2)\sqrt{\frac{f}{a}}$
Must take into account that the number: $p,s$ can be any character.


TITLE: How to plot and shade on Wolfram Alpha?
QUESTION [1 upvotes]: Many Mathematica parameters don't work on Wolfram Alpha. I tried filling -> f which it does not recognize.
I am trying to plot a simple region bounded by $y = x^2$ and $y = x$ and shade that region in red.
I typed in
 plot[x^2,x,{x,0,1}] red filling

and I got something really dumb. I don't want to make those inequality plots either.
Thanks.

REPLY [2 votes]: How about region x^2<x:

Arbitrary characters to exceed arbitrary limit.


TITLE: whether a polynomial is irreducible or not
QUESTION [1 upvotes]: Let $$I=\langle x^4+3x^2+2\rangle$$ and $F=\mathbb{Q}[x]/I$
The question is whether $F$ is a field or not.
I know how to proceed this problem. I'm confusing with irreducibility of the given polynomial
Here, $x^4+3x^2+2=0$ implies $x^4+3x^2=-2$. This happens only $x \in \mathbb{C}$. So the polynomial has no root in $\mathbb{Q}$ and hence irreducible.
However, $x^4+3x^2+2=(x^2+2)(x^2+1)$, so the given polynomial is reducible.
Where I'm doing wrong ?

REPLY [0 votes]: You are assuming that $p(x)$ has no roots in $\mathbb Q$ implies that $p(x)$ is irreducible in $\mathbb Q$. This is true if the degree of $p(x)$ is $2$ or $3$, but not in general, as your own example shows.

REPLY [0 votes]: Irreducibility over $\Bbb Q$ does not mean "has no rational root". As this
example shows, a polynomial may be written as a non-trivial product
but with no roots over $\Bbb Q$. The smallest degree where his happens is $4$ where a polynomial, as here, can be a product of two irreducible
quadratics.


TITLE: Why are there no ellipsoidal drums?
QUESTION [7 upvotes]: It occurred to me today that all drums I could think of have circular heads. It then occurred to me that perhaps an elliptical drumhead would produce different overtones. Do the overtones produced by an elliptical drumhead not sound pleasant? And as a more general question, how does a general drumhead reverberate? What are the equations for a drumhead described by a certain polar curve's, its standing waves?

REPLY [6 votes]: I suspect the actual answer is something boring like ease of manufacture and tuning.
However, one can work out the solutions for the wave equation in elliptic coordinates, perform a separation of variables and end up with a system of differential equations admitting Mathieu functions as solutions. The boundary value problem can be solved relatively straightforwardly [1-4].
(Following [3]) In elliptic coordinates $\xi,\eta$, $x=f \cosh(\xi)\cos(\eta), y=f\sinh(\xi)\sin(\eta)$ where $f$ is the distance from the origin to the foci $(\pm f,0)$. $0\leq \xi < \infty$ is the "radial" coordinate constant on ellipses and $0\leq \eta < 2\pi$ is the "polar" coordinate constant on hyperbolas. The solutions of the wave equation $$\psi(\xi,\eta,t)=T(t)R(\xi)\Theta(\eta)$$ split into the time part $$T''(t)+k^2\nu^2T=0,$$ and the two spatial parts $$R''(\xi)-(\alpha-2q\cosh(2\xi))R(\xi)=0$$ (modified Mathieu equation) and $$\Theta''(\eta)+(\alpha-2q\cos(2\eta))\Theta(\eta)=0$$ (ordinary Mathieu equation) where $-k^2, \alpha$ are the constants of separation and $q=k^2f^2/4$.
One interesting difference from circular membranes is that for each mode there is an even and odd mode, and they oscillate with different frequencies [3]. So if you stimulate one of the modes it is likely to produce a mix of two frequencies that likely do not have a nice rational ratio, and hence does not sound very harmonious. By changing the eccentricity they can likely be made to fit [4], but I suspect this will not make the whole spectrum harmonious.
[1] E. Mathieu, M´e moire sur le mouvement vibratoire d’une membrane de forme elliptique, J. Math. Pures Appl., vol. 13, pp. 137-203, 1868. http://sites.mathdoc.fr/JMPA/PDF/JMPA_1868_2_13_A8_0.pdf
[2] http://booksite.elsevier.com/9780123846549/Chap_Mathieu.pdf
[3] http://optica.mty.itesm.mx/pmog/Papers/P001.pdf
[4] http://www.altenberg.com/peter/pdfs/mathJournalSamp.pdf


TITLE: Power series of analytic function
QUESTION [0 upvotes]: Suppose $f(z)=\sum\limits_{n=0}^{\infty}a_n z^n$ for $|z|<r$, where $r>0$.
Also $f(z)$ is continuous on $|z|\leq r$.
My question is whether the function $f(z)$ is analytic in some "big" disk $|z|<R(R>r)$?
I think the answer is "yes",  because there is no singularity on $|z|=r$, so it can be analytic in a big disk! But I can not give strict proof.
Or some theorem which I do not know can imply my question!
Any help or hint  will welcome!

REPLY [2 votes]: Take $\displaystyle f(z)=\sum_{n=1}^\infty\frac{z^n}{n^2}$. Then $f$ is continuous in the region $\{z\in\mathbb{C}\,|\,|z|\leqslant1\}$. And analytic on $D(0,1)$. However, you cannot extend it to an analytic function on a larger circle.


TITLE: Combination - Inverse Way of solving this problem
QUESTION [1 upvotes]: Regarding my previous post  , I'll repeat the question

A five member committee is to be selected from among four Math teachers and five English teachers. In how many different ways can the committee be formed if the committee must contain at least three Math teachers.

I know it could be solved like this
$\binom{4}{3}\binom{5}{2} + \binom{4}{4} \binom{5}{1} = 45 $ Ans
I wanted to know How I would solve this the other way round. Like for example if it was for at least 1 math teacher it would be
$\binom{9}{5} - \binom{5}{5} $ how would I use the same method but instead calculating for at least 1 I would be calculating for at least 3 ?

REPLY [1 votes]: If I understand you correctly, you want to use that sort of subtraction technique, where you overcount and then take out the complement of the desired probability, to solve the "at least 3" math teacher case.
Then you might do ${9 \choose 5} - \left( {\color{#10C}{5 \choose 5}} +  \color{#070}{{5 \choose 4}{4 \choose 1}} + \color{#C01}{{5 \choose 3}{4 \choose 2}}\right)$, where $\color{#10C}{blue}$ counts the number of ways of choosing 5 English people, $\color{#070}{green}$ counts the number of ways of choosing 4 English and 1 mathie, and $\color{#C01}{red}$ counts the number of ways of choosing 3 English and only 2 mathies.
Of course, this is nothing more than saying that ${9 \choose 5} = {5 \choose 5} + {5 \choose 4}{4 \choose 1} + { 5 \choose 3}{4 \choose 2} + {5 \choose 2}{4 \choose 3} + { 5 \choose 1}{4 \choose 4}$.


text:
\begin{document}

\title{Variation diminishing linear time-invariant systems}
\author{Christian Grussler \thanks{The author is with the University of California Berkeley, Berkeley, CA 94720, USA, 	{\tt\small christian.grussler@berkeley.edu}.} and Rodolphe Sepulchre
	\thanks{The author is with the Control Group at the Department of Engineering, University of Cambridge, Cambridge CB2 1PZ, United Kingdom,
		{\tt\small r.sepulchre @eng.cam.ac.uk}.}}

\maketitle

\input{sec_abstract_RS.tex}
\input{sec_intro.tex}

\input{sec_prelim_merge_revRS.tex}
\input{sec_mot_revised.tex}
\input{sec_vardim_hankel.tex}

\input{sec_vardim_toep}

\input{sec_examples_alt.tex}
\input{sec_discuss.tex}

\section*{Acknowledgment}
The research leading to these results was completed while the first author was a postdoctoral research associate at the University of Cambridge. The research has received funding from the European Research Council under the Advanced ERC Grant Agreement Switchlet n.670645.

\input{sec_appendix.tex}

\bibliographystyle{plain}

\bibliography{refkpos,refopt,refpos,science}

\end{document}