Spectral Gaps with Quantum Counting Queries and Oblivious State Preparation

TL;DR

This paper introduces a quantum algorithm for efficiently approximating spectral gaps and midpoints of Hermitian matrices, offering potential speed-ups over classical methods, especially for large gaps.

Contribution

The work presents a novel quantum algorithm that approximates spectral gaps with logarithmic qubits and analyzes its complexity, including in the black-box model.

Findings

Quantum algorithm approximates spectral gaps up to additive error with logarithmic qubits.

Complexity bounds show potential speed-up over classical algorithms for large spectral gaps.

Established an  lower bound for spectral gap decision in the black-box model.

Abstract

Approximating the $k$-th spectral gap $\Delta_k=|\lambda_k-\lambda_{k+1}|$ and the corresponding midpoint $\mu_k=\frac{\lambda_k+\lambda_{k+1}}{2}$ of an $N\times N$ Hermitian matrix with eigenvalues $\lambda_1\geq\lambda_2\geq\ldots\geq\lambda_N$, is an important special case of the eigenproblem with numerous applications in science and engineering. In this work, we present a quantum algorithm which approximates these values up to additive error $\epsilon\Delta_k$ using a logarithmic number of qubits. Notably, in the QRAM model, its total complexity (queries and gates) is bounded by $O\left( \frac{N^2}{\epsilon^{2}\Delta_k^2}\mathrm{polylog}\left( N,\frac{1}{\Delta_k},\frac{1}{\epsilon},\frac{1}{\delta}\right)\right)$, where $\epsilon,\delta\in(0,1)$ are the accuracy and the failure probability, respectively. For large gaps $\Delta_k$, this provides a speed-up against the best-known complexities of classical algorithms, namely, $O \left( N^{\omega}\mathrm{polylog} \left( N,\frac{1}{\Delta_k},\frac{1}{\epsilon}\right)\right)$, where $\omega\lesssim 2.371$ is the matrix multiplication exponent. A key technical step in the analysis is the preparation of a suitable random initial state, which ultimately allows us to efficiently count the number of eigenvalues that are smaller than a threshold, while maintaining a quadratic complexity in $N$. In the black-box access model, we also report an $\Omega(N^2)$ query lower bound for deciding the existence of a spectral gap in a binary (albeit non-symmetric) matrix.

Equations

$$\mathrm{A}{\mathrm{v}}_i={\lambda }_i{\mathrm{v}}_i.$$

$${\tilde{\Delta }}_k\in (1\pm ϵ){\Delta }_k,\text{and }{\tilde{\mu }}_k\in {\mu }_k\pm ϵ\frac{{\Delta }_k}{2}\subseteq [{\lambda }_{k+1},{\lambda }_k].$$

$$\Vert \mathrm{A}-a\tilde{\mathrm{A}}\Vert \le ϵ,$$

$${\tilde{\lambda }}_{k+1}$$

$$O\left(\left[N^2+N\cdot \frac{a^{2}b\mathcal{C}({\mathrm{U}}_{\mathrm{H}};{ϵ}_{\mathrm{enc}})}{{ϵ}^2{\Delta }_k^2}\right]\mathrm{polylog}\left(N,a,{ϵ}^{-1},{\Delta }_k^{-1},{\delta }^{-1}\right)\right),$$

$$\text{absolutevalue}⟨x|y⟩\ge \frac{1}{N^{O(1)}}.$$

$$\mathrm{Pr}\left[\text{absolutevalue}⟨x|y⟩<\gamma \right]\le \delta .$$

$$\mathrm{Pr}\left[\left|⟨x|y⟩\right|>\frac{1}{2N}\right]\ge \frac{3}{5},$$

$$\mathrm{Pr}[|(\mathrm{H}\mathrm{D}\mathrm{y}{)}_i|\ge t]\le 2\exp (-N{t}^2/2).$$

$$\mathrm{Pr}[\Vert \mathrm{H}\mathrm{D}\mathrm{y}{\Vert }_{\infty }\ge t]\le \delta .$$

$$|{\mathrm{v}}_{m^{\prime }}{|}^2$$

$$|{\mathrm{z}}_A^{\top }\mathrm{v}-{\mathrm{z}}_B^{\top }\mathrm{v}|$$

$$|{\mathrm{z}}_A^{\top }\mathrm{v}-{\mathrm{z}}_C^{\top }\mathrm{v}|$$

$$|{\mathrm{z}}_A^{\top }\mathrm{v}-{\mathrm{z}}_D^{\top }\mathrm{v}|$$

$$|{\mathrm{z}}_B^{\top }\mathrm{v}-{\mathrm{z}}_C^{\top }\mathrm{v}|$$

$$|{\mathrm{z}}_B^{\top }\mathrm{v}-{\mathrm{z}}_D^{\top }\mathrm{v}|$$

$$|{\mathrm{z}}_C^{\top }\mathrm{v}-{\mathrm{z}}_D^{\top }\mathrm{v}|$$

$$A\wedge B\wedge C\wedge D$$

$$\neg$$

$$1-{\delta }_{\mathrm{osp}}-\delta =3/4-{\delta }_{\mathrm{osp}}/2\in (1/2,3/4].$$

$$X_i=\{\begin{array}{l}1,\text{ if }|{\tilde{\sigma }}^2-{\sigma }_N^2|\le {ϵ}_{\sigma },\\ 0,\text{else},\end{array}$$

$$\mathrm{Pr}[X\le r/2]<\exp (-r(p-1/2))\le \exp (-r(1-2{\delta }_{\mathrm{osp}})/4).$$

$$\mathcal{C}({\mathrm{U}}_{\mathrm{sgn}};\Delta ,{ϵ}_{\mathrm{enc}})$$

$$O\left(\mathcal{C}({\mathrm{U}}_{\mathrm{sgn}};\Delta ,{ϵ}_{\mathrm{enc}})\frac{aN}{{ϵ}_{\mathrm{tr}}}\log (\frac{1}{\delta })\right)$$

$$=O\left(\left(\mathcal{C}({\mathrm{U}}_{\mathrm{A}};{ϵ}_{\mathrm{enc}})+b\right)\frac{a^{2}N}{\Delta {ϵ}_{\mathrm{tr}}}\log \left(\frac{1}{{ϵ}_{\mathrm{sgn}}}\right)\log (\frac{1}{\delta })\right).$$

$$O\left(\left(\mathcal{C}({\mathrm{U}}_{\mathrm{A}};{ϵ}_{\mathrm{enc}})+b\right)\cdot \frac{a^{2}N}{\Delta }\log \left(N\right)\log (\frac{1}{\delta })\right).$$

$$1/2\ge {\lambda }_1\ge {\lambda }_2\ge \dots \ge {\lambda }_N\ge -1/2.$$

$$\left|{\tilde{\mu }}_k-{\mu }_k\right|\le \frac{7}{16}{\Delta }_k,{\tilde{\Delta }}_k\in [{\Delta }_k/2,4{\Delta }_k].$$

$$Q_1:=O\left(\frac{a^2}{\gamma {\Delta }_k^2(1-2{\delta }_{\mathrm{osp}})}\log (\frac{1}{\gamma })\log (\frac{a}{{\Delta }_k})\log (\frac{\log (a/{\Delta }_k)}{1-2{\delta }_{\mathrm{osp}}})\log (\frac{1}{\delta })\log (\frac{1}{{\Delta }_k})\right),$$

$$O\left(\frac{1}{\gamma (1-2{\delta }_{\mathrm{osp}})}\log (\frac{a}{{\Delta }_k})\log (\frac{\log (a/{\Delta }_k)}{1-2{\delta }_{\mathrm{osp}}})\log (\frac{1}{\delta })\log (\frac{1}{{\Delta }_k})\right),$$

$$Q_2:=O\left(\frac{a^{2}N}{{\Delta }_k}\log \left(N\right)\log (\frac{1}{{\Delta }_k\delta })\log (\frac{1}{{\Delta }_k})\right),$$

$$\mathrm{Pr}\left[\left|{\tilde{\sigma }}_i^2-{\sigma }_{\min }^2(\mathrm{H}-h_i)\right|>\frac{{\Delta }^{2j}}{64}\right]\le \frac{{\delta }_j}{2}.$$

Full text

1 Introduction

Eigenvalue problems lie in the heart of computational science, engineering, and machine learning. Given an $N\times N$ input matrix $\boldsymbol{\mathrm{A}}$, the goal is to compute (a subset of) the eigenvalues $\lambda_{i}$ and/or the eigenvectors $\boldsymbol{\mathrm{v}}_{i}$, for $i=1,\ldots,N$, which satisfy:

[TABLE]

In this work we focus on algorithms for an important special case of the eigenproblem, namely, the approximation of the so-called spectral gap (or eigengap) and the corresponding midpoint between pairs of consecutive eigenvalues of Hermitian matrices, which we formally define as follows.

Problem 1.1** (Spectral gap and midpoint approximation).**

Let $\boldsymbol{\mathrm{H}}$ be an $N\times N$ Hermitian matrix, with eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{N}$, and $\operatorname{\Delta}_{k}:=|\lambda_{k}-\lambda_{k+1}|$ be the $k$-th spectral gap and $\mu_{k}:=\frac{\lambda_{k}+\lambda_{k+1}}{2}$ is the corresponding midpoint between $\lambda_{k}$ and $\lambda_{k+1}$. Given $\boldsymbol{\mathrm{H}}$, an integer $k\in[N-1]$, and an accuracy $\epsilon\in(0,1)$, compute two values $\widetilde{\operatorname{\Delta}}_{k},\widetilde{\mu}_{k}$ such that:

[TABLE]

No prior information regarding the eigenvalues or the eigenvectors of $\boldsymbol{\mathrm{H}}$ is assumed.

Problem 1.1 has numerous applications in different areas. In numerical analysis, estimating spectral gaps is a key subroutine of so-called spectral divide-and-conquer algorithms for diagonalizing matrices [1, 2, 3, 4]. These methods rely on finding a “splitting” point $\mu$ that does not coincide with an eigenvalue; our midpoint estimator $\mu_{k}=(\lambda_{k}+\lambda_{k+1})/2$ provides a natural and well-conditioned choice for such a split, as it maximizes the distance to the nearest eigenvalues. In data analytics, a large gap in the spectrum of a data matrix indicates a separation between the actual information and the “noise” in the data. This phenomenon is often used in Principal Component Analysis [5, 6] to identify good splitting points, or “elbows”, for dimensionality reduction, data visualization, and noise filtering. By directly targeting the $k$-th gap, our algorithm enables the automated detection of these elbows without requiring the explicit computation of (potentially large) subsets of singular values or a full singular value decomposition (SVD). In spectral graph theory, gaps between the eigenvalues of the graph Laplacian reveal connectivity properties, and they can be used, for example, to identify communities, [7, 8, 9, 10].

In the study of quantum phase transitions, the closing of the $k$-th gap, where $\lambda_{k+1}-\lambda_{k}\to 0$ as the system size increases, serves as a key indicator of topological or magnetic phase transitions [11]. In the context of quantum computing, the size of the spectral gap also determines the runtime of adiabatic state preparation, since the adiabatic theorem requires the evolution time to scale inversely with (a power of) the minimum spectral gap encountered [12]. Last but not least, spectral gaps are also important in electronic structure calculations, for example, to estimate the energy gap (Fermi level) that separates the occupied and unoccupied states of the underlying system [13, 14, 15, 16]. In that direction, our gap estimator could be combined with the techniques of [17] to calculate the electron density of a system at zero temperature.

The aim of this work is to provide a detailed analysis of quantum algorithms for solving Problem 1.1, aiming for provable approximation guarantees and end-to-end worst-case complexity bounds.

1.1 Existing algorithms and complexity bounds

Here we give a brief overview of the related algorithms literature. We first recall that, from the Abel-Ruffini theorem, computing eigenvalues of general matrices is intractable in the algebraic model; see e.g. [18]. Besides trivial cases, such as triangular matrices, algorithms for eigenvalue problems in general return approximate solutions. The computational complexity depends on several different factors, namely, the underlying model of computation, the properties of the input matrix (for example, Hermitian, sparse, structured, or low-rank), the quantities that need to be computed (subset of eigenvalues/eigenvectors), the accuracy, and, for randomized algorithms, the success probability.

In the classical (non-quantum) regime, recent works have established that both the arithmetic and boolean worst-case complexity of diagonalization (approximation of all the eigenvalues and eigenvectors), is upper bounded by $\widetilde{O}(N^{\omega})$ operations [3, 19, 20, 21, 22, 4]. Here $\omega$ is the matrix multiplication exponent, currently upper bounded by $\omega<2.371339$ [23], and $\widetilde{O}$ hides polylogarithmic factors in $N$ and potentially other parameters, such as the accuracy. In terms of solving Problem 1.1 per-se, [20] reports a complexity of $O\left(N^{\omega}\operatorname{polylog}\left(N,\frac{1}{\operatorname{\Delta}_{k}},\frac{1}{\epsilon}\right)\right)$ operations in the Boolean RAM model.

Beyond the matrix multiplication-type worst-case complexity, fast algorithms have been studied for special classes of matrices. For instance, the eigenvalues of banded symmetric matrices with $k$ bands can be computed in $\widetilde{O}(N^{2}k^{\omega-2})$ time via tridiagonal reduction [20, 24]. For symmetric diagonally dominant matrices (SDD) and graph Laplacians, the smallest $k$ eigenvalues of an $N\times N$ SDD matrix with $M$ non-zeros can be approximated in $\widetilde{O}(M+kN/\epsilon^{2})$ time [25, 26]. This is faster than the aforementioned $\widetilde{O}(N^{\omega})$ bounds even for $k\in\Theta(N)$ and $M\in\Theta(N^{2})$, albeit with increased dependence on the accuracy $\epsilon$ (polynomial instead of logarithmic). The recent work of [27] described an elegant approach specifically targeting spectral gaps of general (non-SDD) sparse symmetric matrices based on stochastic trace estimation [28, 29, 30, 31], albeit, the current worst-case complexity analysis can grow larger than $\widetilde{O}(N^{\omega})$.

The literature on quantum algorithms for eigenvalue problems is vast, and, to date, it keeps growing rapidly. The majority of existing works has focused on the following problems: (i) estimation of extremal eigenvalues, i.e., the ground state energy, (ii) approximating an eigenvalue, given access to the corresponding eigenvector, or (iii) preparation of quantum states that encode eigenvalues as amplitudes. Such methods include the so-called Variational Quantum Eigensolvers [32, 33, 34], and methods based on Quantum Phase Estimation (QPE) [35, 36, 37], such as [38, 39, 40, 41, 42, 43, 44].

Here we highlight recent works that can be considered for solving Problem 1.1, with explicit worst-case complexity bounds. Low and Su [39] (see e.g. Theorem 3) recently proved sharp quantum query bounds for estimating any eigenvalue of general matrices. The query bound is, $O((a\epsilon)^{-1}\log(1/\delta))$ for Hermitian matrices, where $a$ is the block-encoding scaling factor. The algorithm assumes as input an initial state that is $\epsilon$-close (in the Euclidean distance) to the corresponding eigenvector. Such a state is hard to obtain in the general case, as it entails to approximate an (unknown) eigenvector.

State preparation assumptions can be relaxed for the problem of estimating the smallest eigenvalues. Recent works of [38, 45, 44, 40] suggest that $1/\operatorname{poly}(N)$ overlap between the initial state and the target bottom-$k$ states is sufficient. If an initial state that satisfies this requirement is readily available, then a straightforward approach to solve Problem 1.1 is to directly [45, 44] to approximate all $N-k$ smallest eigenvalues, in a “brute-force” manner. As it is clearly noted in these works, however, this approach is targeted for approximating few eigenvalues close to the end of the spectrum, in this case, when $N-k$ is small. Moreover, beyond the initial state preparation assumptions, it is required that the $N-k$ smallest eigenvalues are strictly separated by a well-defined minimum gap. For large $N-k\in\Theta(N)$, and for general Hermitian matrices, it becomes challenging to ensure or to override these assumptions and to guarantee competitive end-to-end complexity bounds in the worst-case.

End-to-end worst-case complexity bounds were proved in recent works by Apers and de Wolf [46], and by Chen, Gilyén and de Wolf [47]. The former showed how to approximate the top-$k$ eigenvalues and eigenvectors of SDD matrices with $M$ non-zeros in $\widetilde{O}(\sqrt{MN}+kN/\epsilon^{2})$ time in the QRAM model. It is faster than the $O(M+kN/\epsilon^{2})$ classical SDD algorithms for small $k$ and large $M$. For general square matrices, [47] showed how to approximate the top-$k$ eigenvectors in $\widetilde{O}(kN^{1.5+o(1)})$, using quantum versions of the power method. As also noted by the authors, however, the algorithm is less competitive for large $k\in\Theta(N)$.

As a last remark, we note that obtaining initial states which sufficiently approximate the target eigenvectors is an interesting research topic on its own. For certain applications, an initial state can be obtained by solving an “easier” or “smaller” version of the problem as a pre-processing step. As [38] notes, this is particularly interesting for quantum chemistry applications [48]. Random states are often recommended for the general case, but they can be also expensive to prepare, as we also discuss here.

1.2 Contributions and methodology

Our main contribution is an algorithm (Algorithm 2) to solve Problem 1.1. In a high-level, it consists of two phases:

• •

In the first phase (Alg. 1), the spectral gap $\operatorname{\Delta}_{k}$ and the midpoint $\mu_{k}$ are computed to constant accuracy (i.e., a crude initial approximation).

• •

The second phase (Alg. 2) uses the values returned in the first phase in order to perform a refined search in the restricted domain and increase the solution accuracy.

Both phases use the following steps to locate gaps and/or eigenvalues:

Encode the given Hermitian matrix in a quantum circuit. 2. 2.

Split the spectrum in two parts, by shifting the diagonal elements to a carefully chosen point. 3. 3.

Compute the spectral projector on each part of the spectrum with a polynomial approximation of the sign function. 4. 4.

Compute the trace of the projector. 5. 5.

Repeat 2-4 as a search mechanism (“counting queries”) over a grid to locate eigenvalues and/or gaps.

For Step 1, we use the block-encoding framework [49, 50], which has recently received significant attention for matrix problems. We use the following definition.

Definition 1.1** (Block-Encoding).**

Let $\boldsymbol{\mathrm{A}}\in\mathbb{C}^{N\times N}$ be an $s$-qubit matrix, $a\geq\|\boldsymbol{\mathrm{A}}\|,\epsilon>0$, and $b\in\mathbb{N}$. An $(s+b)$-Qubit Unitary $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}}$ is an $(a,b,\epsilon)$-Block-Encoding of $\boldsymbol{\mathrm{A}}$, which is also abbreviated as $(a,b,\epsilon)$-QUBE for short, if

[TABLE]

where $\widetilde{\boldsymbol{\mathrm{A}}}:=(\bra{0}^{\otimes b}\otimes\boldsymbol{\mathrm{I}})\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}}(\ket{0}^{\otimes b}\otimes\boldsymbol{\mathrm{I}})$.

Our main result is summarized in the following theorem.

Theorem 1.1** (Informal, see Theorem 3.1).**

Let $\boldsymbol{\mathrm{H}}$ be an $N\times N$ Hermitian matrix with eigenvalues $-1/2\leq\lambda_{N}\leq\lambda_{N-1}\leq\ldots\leq\lambda_{1}\leq 1/2$, $\operatorname{\Delta}_{k}:=|\lambda_{k}-\lambda_{k+1}|$, and $\mu_{k}:=\frac{\lambda_{k}+\lambda_{k+1}}{2}$, and let $\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon_{\rm enc})$ denote the cost (gate complexity) of the unitary $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}$ which forms an $(a,b,\epsilon_{\rm enc})$-block-encoding of $\boldsymbol{\mathrm{H}}$. The proposed Algorithm 2 takes as input $\boldsymbol{\mathrm{H}}$, an integer $k\in[N-1]$, an accuracy $\epsilon\in(0,1)$, and a failure probability parameter $\delta\in(0,1/2)$, and it returns values $\widetilde{\lambda}_{k},\widetilde{\lambda}_{k+1},\widetilde{\mu}_{k},$ and $\widetilde{\operatorname{\Delta}}_{k}$, which satisfy:

[TABLE]

The total complexity, including all queries to unitaries and auxiliary gates, is bounded by:

[TABLE]

where $\epsilon_{\rm enc}\in\Theta\left(\frac{\epsilon^{2}\operatorname{\Delta}_{k}^{2}}{a^{2}N}\right)$.

Notably, in the QRAM model, the complexity becomes $\widetilde{O}(N^{2}\operatorname{\Delta}_{k}^{-2}\epsilon^{-2})$. This is faster than the $\widetilde{O}(N^{\omega})$ arithmetic complexity of the best-known classical algorithms, for all matrices with sufficiently large gaps. This provides evidence that quantum algorithms can beat “matrix multiplication-time” for these problems, at least for the current bound of $\omega$.

While the general idea is simple, there are several non-trivial parts that need to be addressed to obtain competitive end-to-end complexities. Two major challenges, which can be of independent interest, were the choice of the splitting points, and the preparation of appropriate initial states, which we discuss further below.

Informed choice of splitting points via condition number estimation

As mentioned, the splitting points to split the spectrum and count eigenvalues. If the chosen point is close to an actual eigenvalue, the condition number of the shifted matrix can grow arbitrarily, and it becomes prohibitively expensive to approximate the corresponding spectral projector. Indeed, the convergence of iterative methods, such as the Newton-Schultz iteration or other rational approximations [51], scales as $\Omega(1/\operatorname{\Delta})$, where $\operatorname{\Delta}$ is the distance of the splitting point from the closest eigenvalue [49]. For this reason, many existing algorithms that follow this approach require prior knowledge of $\operatorname{\Delta}$ (or at least a suitable lower bound), which is generally not trivial to obtain.

To avoid this problem, we use a dedicated subroutine which tests whether the gap $\operatorname{\Delta}$ is larger than a well-defined threshold, before executing any counting query. This is achieved by conveniently applying the ground-energy estimation algorithm of [38] in order to estimate the smallest singular value of the matrix, after shifting its diagonal elements to the chosen splitting point. If this singular value is “large enough”, it means that we can safely count eigenvalues. If not, we keep searching for smaller gaps. This way we have full control of the value of $\operatorname{\Delta}$, and, consequently, of the cost of eigenvalue counting.

Initial state preparation

In order to execute our singular value estimation subroutine efficiently, we need to address a second major challenge, which involves the preparation of an appropriate initial state. The algorithm of [38] requires as input an initial state $\ket{x}$ which will overlap appropriately with the corresponding singular vector. In particular, this target singular vector can be any fixed state $\ket{y}$, and we need to ensure that the initial guess $\ket{x}$ will satisfy

[TABLE]

As the authors of [38, 45] note, such an initial state can be efficiently prepared for certain classes of problems. For instance, if $\ket{y}$ is a binary string with exactly one nonzero entry, a Hadamard state is sufficient, or, in quantum chemistry applications, the Hartree-Fock state can be pre-computed as a “cheap” initial approximation to the ground state [48].

In the general case, where $\ket{y}$ can be any arbitrary state without any explicit prior information, the analysis becomes more involved. In the context of QPE, a common recommendation is to use a random state. However, if the random distribution is not chosen carefully, it can be harmful rather than beneficial. Contributing to this direction, in Section 2.1 we analyze a randomized state preparation method, which achieves all the desired properties. It is oblivious to the target state, in a sense that it guarantees the desired overlap $1/N$, for any fixed target state, with probability at least $3/5$. Our construction is an ensemble of a random Rademacher state, together with a randomized Hadamard transform [52, 53, 54]. The latter is a well-established tool in Randomized Numerical Linear Algebra, and it is crucial for our analysis do to its “flattening” properties, commonly attributed as a consequence of the Heisenberg principle.

The size of the final quantum circuit scales as $O(N)$, and it requires a total of $O(N)$ random classical bits to prepare. While this cost is prohibitive for applications that target exponential quantum speed-ups, it is “sufficiently cheap” for our analysis, and it can be useful for other works that target polynomial quantum speed-ups.

Lower bounds

The aforementioned upper bounds are also supported with a $\Omega(N^{2})$ query lower bound for the problem of deciding whether there exists a gap between two consecutive eigenvalues of a (non-symmetric) binary matrix in the black-box access model. The techniques used here are based on “majority” arguments, inspired by [55, 56, 57]. Evidently, these lower bounds are not tight, as they do not match the reported upper bound. In order to close this complexity gap, we think that one should use entirely different techniques, both for the upper and the lower bounds. This is left as a major open question.

1.3 Notation

$\mathcal{H}_{S}$ denotes an explicit Hilbert space and $\mathcal{H}_{A}$ an ancillary space for post-selection. When we want to emphasize its dimension $d$, we write $\mathcal{H}_{A}\cong\mathbb{C}^{d}$. The total Hilbert space is given by $\mathcal{H}=\mathcal{H}_{S}\otimes\mathcal{H}_{A}$. Matrices and vectors are denoted with bold capital and small letters, respectively. For a matrix $\boldsymbol{\mathrm{A}}$, $\boldsymbol{\mathrm{A}}^{\dagger}$ is the conjugate transpose, $\lVert\boldsymbol{\mathrm{A}}\rVert$ denotes its spectral norm, $\lVert\boldsymbol{\mathrm{A}}\rVert_{F}$ its Frobenius norm, and $\Lambda(\boldsymbol{\mathrm{A}})$ denotes its spectrum. $\mathcal{C}(\boldsymbol{\mathrm{U}})$ denotes the gate complexity of the unitary $\boldsymbol{\mathrm{U}}$, and we will write $\mathcal{C}(\boldsymbol{\mathrm{U}};\epsilon,\delta)$ when the complexity of $\boldsymbol{\mathrm{U}}$ is parametrized by $\epsilon,\delta$.

1.4 Model of computation

The complexity of a quantum algorithm is two-fold. We define the query complexity to be the total number of calls to the block-encoding of the original matrix. The gate complexity is the total number of elementary one- and two-qubit quantum gates, i.e., the size of the quantum circuit. For all algorithms, we explicitly mention both the query and the gate complexities, taking into account all the involved parameters, such as the matrix size, the accuracy, and the success probabilities. For our state-preparation analysis, we also assume access to a (classical) random number generator that produces independent, uniformly random bits.

In Section 4, we specialize the analysis in the so-called Quantum-RAM (or QRAM) model [58]. In this model the quantum processor has access to a memory that stores classical bits. These bits can be read in superposition, but we can only write zeros and ones, just as in the classical case. This model is not necessary for our main analysis, but we dedicate Section 4 to it since it is used in theoretical works, e.g., for graph and matrix computations [46, 47, 41]. For a discussion on QRAM architectures we refer to [59, 60].

1.5 Outline

Section˜2 contains the analysis of the main subroutines that are used by the main algorithms, namely, for state preparation, computing the minimum singular value, and counting eigenvalues that are smaller than a threshold. Section˜3 contains the main algorithm for computing spectral gaps and eigenvalue pairs. In Section˜4 we report the analysis in the QRAM model. In Section˜5 we prove a query lower bound for deciding the existence of a spectral gap in binary matrices in the black-box access model. We finally provide some concluding remarks in Section˜6.

Acknowledgements

We are grateful to Ivano Tavernelli and Simon Apers for valuable discussions.

2 Main Subroutines

Here we describe the main subroutines that are used as building blocks for our main algorithms. Some of the techniques reported here rely on existing results on block-encodings and quantum eigenvalue transformation from [49, 61, 38, 50], which are summarized in Appendix A.

2.1 Oblivious State Preparation

In many applications of QPE, it is required to prepare an initial state $\ket{x}$ that is not orthogonal to the target state $\ket{y}$ (e.g. the ground state). A sufficient condition is typically $\absolutevalue{\braket{x|y}}^{2}\geq\Omega(1/N)$. The state $\ket{x}=\ket{+}:=\boldsymbol{\mathrm{H}}\ket{0}$ is sufficient for many applications, e.g., when the desired state is a binary string. For arbitrary $\ket{x}$, however, there is no theoretical guarantee that $\ket{+}$ is not orthogonal to $\ket{x}$. For example, it is known that $\ket{+}$ spans the kernel of a connected graph Laplacian, and it is therefore orthogonal to every eigenvector that corresponds to a nonzero eigenvalue.

A reasonable approach to circumvent this limitation is to use random states. However, as already mentioned previously, this needs to be handled with care. Here we describe a random state preparation, inspired by well-known techniques from randomized matrix algorithms [52]. Before stating the main result, we define a state-preparation unitary.

Definition 2.1**.**

A distribution $\boldsymbol{\mathrm{\Pi}}$ of random quantum circuits forms a $(\gamma,\delta)$-Oblivious State Preparation, or $(\gamma,\delta)$-OSP, if the unitary $\boldsymbol{\mathrm{U}}_{S}\sim\boldsymbol{\mathrm{\Pi}}$ operates on $O(\log(N))$ qubits and it prepares a state $\ket{x}=\boldsymbol{\mathrm{U}}_{S}\ket{0}$ such that, for any fixed state $\ket{y}\in\mathbb{C}^{N}$:

[TABLE]

Lemma 2.1**.**

Let $\ket{y}\in\mathbb{C}^{N}$ be any fixed state (e.g. the target ground-state). Let $\ket{x}=\boldsymbol{\mathrm{D}}\boldsymbol{\mathrm{H}}\ket{z}$. Here $\boldsymbol{\mathrm{D}}$ is diagonal with independent $\{\pm 1\}$ diagonal elements, $\boldsymbol{\mathrm{H}}:=\begin{pmatrix}\tfrac{1}{\sqrt{2}}&\tfrac{1}{\sqrt{2}}\\ \tfrac{1}{\sqrt{2}}&-\tfrac{1}{\sqrt{2}}\end{pmatrix}^{\otimes\log(N)}$ is a Hadamard unitary, and $\ket{z}\in\left\{\pm\frac{1}{\sqrt{N}}\right\}^{N}$ has independent Rademacher elements. Let $\boldsymbol{\mathrm{Z}}$ be a unitary such that $\boldsymbol{\mathrm{Z}}\ket{0}=\ket{z}$. Then, for all $N\geq C$, where $C\approx 22$ is a fixed constant, the distribution of unitaries $\boldsymbol{\mathrm{U}}_{S}=\boldsymbol{\mathrm{D}}\boldsymbol{\mathrm{H}}\boldsymbol{\mathrm{Z}}$ is an $(\frac{1}{2N},\frac{2}{5})$-OSP, that is,

[TABLE]

The state $\ket{x}$ can be prepared with $O(N)$ elementary gates, i.e., $\mathcal{C}(\boldsymbol{\mathrm{U}}_{S})\in O(N)$.

Proof.

Let us denote $\boldsymbol{\mathrm{x}}=\ket{x},\boldsymbol{\mathrm{z}}=\ket{z}$, and $\boldsymbol{\mathrm{y}}=\ket{y}$ to simplify the notation. The operator $\boldsymbol{\mathrm{H}}\boldsymbol{\mathrm{D}}$ is a well-studied tool in randomized linear algebra (see e.g. [52]), which has the property that it “flattens” the elements of any fixed vector with high probability. Indeed, using Hoeffding’s inequality, we can see that for any index $i\in[N]$ it holds that

[TABLE]

Taking a union bound for all $i\in[N]$ and setting $t=\sqrt{\frac{2\ln(2N/\delta)}{N}}$, we can conclude that

[TABLE]

Denote $\boldsymbol{\mathrm{v}}:=\boldsymbol{\mathrm{H}}\boldsymbol{\mathrm{D}}\boldsymbol{\mathrm{y}}$. It remains to get a lower bound on the overlap $|\boldsymbol{\mathrm{x}}^{\top}\boldsymbol{\mathrm{y}}|=|\boldsymbol{\mathrm{z}}^{\top}\boldsymbol{\mathrm{v}}|$. Since $\boldsymbol{\mathrm{H}}\boldsymbol{\mathrm{D}}$ is unitary, $\boldsymbol{\mathrm{v}}$ is a unit vector with flattened entries, which is the key observation to obtain the desired bound. Let $\boldsymbol{\mathrm{v}}_{m}$ and $\boldsymbol{\mathrm{v}}_{m^{\prime}}$ be the largest and second largest elements of $\boldsymbol{\mathrm{v}}$ in modulus, where $m,m^{\prime}\in[N]$. It holds that $|\boldsymbol{\mathrm{v}}_{m}|^{2}\in[1/N,t^{2}]$. For $\boldsymbol{\mathrm{v}}_{m^{\prime}}$, observe that

[TABLE]

For $\delta=1/5$ and $N\gtrsim 22$, we have that $|\boldsymbol{\mathrm{v}}_{m^{\prime}}|^{2}\geq 1/(2N-2)>1/(2N)$. We now condition on the success of the flattening, and proceed to get a bound for the overlap. Let $\boldsymbol{\mathrm{z}}_{A},\boldsymbol{\mathrm{z}}_{B},\boldsymbol{\mathrm{z}}_{C},\boldsymbol{\mathrm{z}}_{D}$ be four vectors that have exactly the same signs at all the positions except $\{m,m^{\prime}\}$, where in these two positions they have the signs $\{+1,+1\},\{+1,-1\},\{-1,+1\},\{-1,-1\}$, respectively. The following relations hold:

[TABLE]

Now, denote $L:=1/(2N)$, and let $A$ be the event that $|\boldsymbol{\mathrm{z}}_{A}^{\top}\boldsymbol{\mathrm{v}}|>L$, $B:|\boldsymbol{\mathrm{z}}_{B}^{\top}\boldsymbol{\mathrm{v}}|>L$, $C:|\boldsymbol{\mathrm{z}}_{C}^{\top}\boldsymbol{\mathrm{v}}|>L$, and $D:|\boldsymbol{\mathrm{z}}_{D}^{\top}\boldsymbol{\mathrm{v}}|>L$. From the relations above, we can see that (deterministically), we are in one of the following five cases:

[TABLE]

In all cases, at least three out of four vectors $\{\boldsymbol{\mathrm{z}}_{A},\boldsymbol{\mathrm{z}}_{B},\boldsymbol{\mathrm{z}}_{C},\boldsymbol{\mathrm{z}}_{D}\}$ satisfy the desired lower bound (deterministically). Then, since we pick any of the four vectors uniformly at random, the success probability is at least $3/4$.

Since the flattening succeeds with probability $4/5$, both events succeed with probability at least $\frac{3}{4}\cdot\frac{4}{5}=\frac{3}{5}$, which concludes the proof. ∎

2.2 Computing the smallest singular value

Here we describe the subroutine which computes the smallest singular value of a block-encoded matrix. [38] gave an algorithm with near optimal dependence on the precision $\epsilon$ to find the smallest eigenvalue of a Hermitian matrix $\boldsymbol{\mathrm{H}}$.

Theorem 2.1** (Ground energy, Imported Thm. 8 from [38]).**

Let $\boldsymbol{\mathrm{H}}=\sum_{k}\lambda_{k}\ket{\psi_{k}}\bra{\psi_{k}}\in\mathbb{C}^{N\times N}$, where $\lambda_{k}\geq\lambda_{k+1}$, be a Hermitian matrix given through its $(\alpha,m,0)$-block-encoding $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}$. Let $\ket{\phi_{0}}$ be an initial state prepared by circuit $\boldsymbol{\mathrm{U}}_{S}$, with the promise $\absolutevalue{\braket{\phi_{0}|\psi_{0}}}\geq\gamma$. Then the ground energy $\lambda_{N}$ can be estimated to additive precision $\epsilon$ with probability $1-\delta$ with the following costs:

Queries to $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}$: $O\left(\frac{\alpha}{\gamma\epsilon}\log(\frac{\alpha}{\epsilon})\log(\frac{1}{\gamma})\log(\frac{\log(\alpha/\epsilon)}{\delta})\right)$, 2. 2.

Queries to $\boldsymbol{\mathrm{U}}_{S}$: $O\left(\frac{1}{\gamma}\log(\frac{\alpha}{\epsilon})\log(\frac{\log(\alpha/\epsilon)}{\delta})\right)$, 3. 3.

Number of qubits: $O\left(\log(N)+m+\log(\frac{1}{\gamma})\right)$, 4. 4.

Other elementary gates: $O\left(\frac{m\alpha}{\gamma\epsilon}\log(\frac{\alpha}{\epsilon})\log(\frac{1}{\gamma})\log(\frac{\log(\alpha/\epsilon)}{\delta})\right)$.

This can be applied appropriately on the Gramian $\boldsymbol{\mathrm{H}}^{\dagger}\boldsymbol{\mathrm{H}}$ to compute the minimum singular value.

Lemma 2.2** ($\mathsf{QSMIN}$).**

Let $\boldsymbol{\mathrm{H}}\in\mathbb{C}^{N\times N}$ be a Hermitian matrix, with singular values $0\leq\sigma_{N}\leq\sigma_{N-1}\leq\ldots\leq\sigma_{1}\leq 1$. Let $\epsilon_{\sigma}\in(0,1/2]$ be an accuracy parameter, and assume that we have access to:

• •

$\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}$, an $(a,b,\epsilon_{\rm enc})$-QUBE for $\boldsymbol{\mathrm{H}}$, $\epsilon_{\rm enc}\in[0,\epsilon_{\sigma}/4a]$, with cost $\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon_{\rm enc})$

• •

$\boldsymbol{\mathrm{\Pi}}$, a $(\gamma,\delta_{\rm osp})$-OSP, where $\delta_{\rm osp}\in[0,1/2)$, and $\boldsymbol{\mathrm{U}}_{S}\sim\boldsymbol{\mathrm{\Pi}}$ has cost $\mathcal{C}(\boldsymbol{\mathrm{U}}_{S})$.

We can compute a boolean string $X\leftarrow\mathsf{QSMIN}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}},\epsilon_{\sigma},\delta_{\sigma})$, such that $\widetilde{\sigma}^{2}:=0.X\in[0,1]$ satisfies $|\widetilde{\sigma}^{2}-\sigma_{N}^{2}|\leq\epsilon_{\sigma}$ with probability at least $1-\delta_{\sigma}$, with the following costs:

Queries to $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}$: $O\left(\frac{a^{2}}{\gamma\epsilon_{\sigma}}\log(\frac{1}{\gamma})\log(\frac{a}{\epsilon_{\sigma}})\log(\frac{\log(a/\epsilon_{\sigma})}{1-2\delta_{\rm osp}})\frac{\log(1/\delta_{\sigma})}{1-2\delta_{\rm osp}}\right)$, 2. 2.

Queries to random circuits $\boldsymbol{\mathrm{U}}_{S}\sim\boldsymbol{\mathrm{\Pi}}$: $O\left(\frac{1}{\gamma}\log(\frac{a}{\epsilon_{\sigma}})\log(\frac{\log(a/\epsilon_{\sigma})}{1-2\delta_{\rm osp}})\frac{\log(1/\delta_{\sigma})}{1-2\delta_{\rm osp}}\right)$, 3. 3.

Number of qubits: $O\left(\log(N)+b+\log(\frac{1}{\gamma})\right)$, 4. 4.

Other elementary gates: $O\left(\frac{a^{2}b}{\gamma\epsilon_{\sigma}}\log(\frac{1}{\gamma})\log(\frac{a}{\epsilon_{\sigma}})\log(\frac{\log(a/\epsilon_{\sigma})}{1-2\delta_{\rm osp}})\frac{\log(1/\delta_{\sigma})}{1-2\delta_{\rm osp}}\right)$.

Proof.

Given $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}$, an $(a,b,\epsilon_{\rm enc})$-QUBE for $\boldsymbol{\mathrm{H}}$ with cost $\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon_{\rm enc})$, it holds that $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}^{\dagger}$ is an $(a,b,\epsilon_{\rm enc})$-QUBE for $\boldsymbol{\mathrm{H}}^{\dagger}$, also with cost $\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon_{\rm enc})$. We can combine these two block-encodings as in [50, Lemma 4] to create an $(a^{2},2b,2a\epsilon_{\rm enc})$-QUBE for $\boldsymbol{\mathrm{H}^{\dagger}\mathrm{H}}$, denoted by $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}^{\dagger}\mathrm{H}}}$,with cost equal to $\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}^{\dagger}\mathrm{H}}};\epsilon_{\rm enc})=2\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon_{\rm enc})$. The minimum eigenvalue of $\boldsymbol{\mathrm{H}^{\dagger}\mathrm{H}}$ is the square of the minimum singular value of $\boldsymbol{\mathrm{H}}$.

The following steps are performed. We first prepare a state $\ket{\psi}=\boldsymbol{\mathrm{U}}_{S}\ket{0}$ where $\boldsymbol{\mathrm{U}}_{S}\sim\boldsymbol{\mathrm{\Pi}}$. Then, we use Theorem 2.1 with inputs $\ket{\psi}$, $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}^{\dagger}\mathrm{H}}}$, $\epsilon=\epsilon_{\sigma}/2$, and $\delta=1/4-\delta_{\rm osp}/2\in(0,1/4]$. The algorithm returns a bit-string $\widetilde{\sigma}^{2}\in[0,1]$ which satisfies $|\widetilde{\sigma}^{2}-\sigma_{N}(\widetilde{\boldsymbol{\mathrm{H}^{\dagger}\mathrm{H}}})|\leq\epsilon_{\sigma}/2$ with probability at least:

[TABLE]

Here $\widetilde{\boldsymbol{\mathrm{H}^{\dagger}\mathrm{H}}}$ is the block-encoded version of $\boldsymbol{\mathrm{H}}^{\dagger}\boldsymbol{\mathrm{H}}$, which satisfies $\|\widetilde{\boldsymbol{\mathrm{H}^{\dagger}\mathrm{H}}}-a^{2}\boldsymbol{\mathrm{H}^{\dagger}\mathrm{H}}\|\leq 2a\epsilon_{H}$.

From the stability of singular values, which is a consequence of Weyl’s inequality, we conclude that $|\widetilde{\sigma}^{2}-\sigma_{N}^{2}(\boldsymbol{\mathrm{H}})|\leq 2a\epsilon_{H}+\epsilon_{\sigma}/2\leq\epsilon_{\sigma}$, due to the assumption that $\epsilon_{\rm enc}\leq\epsilon_{\sigma}/4a$. From Theorem 2.1, the cost of such a run is bounded by

• •

$O\left(\frac{a^{2}}{\gamma\epsilon_{\sigma}}\log(\frac{1}{\gamma})\log(\frac{a}{\epsilon_{\sigma}})\log(\frac{\log(a/\epsilon_{\sigma})}{1-2\delta_{\rm osp}})\right)$ queries to $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}^{\dagger}\boldsymbol{\mathrm{H}}}$,

• •

$O\left(\frac{1}{\gamma}\log(\frac{a}{\epsilon_{\sigma}})\log(\frac{\log(a/\epsilon_{\sigma})}{1-2\delta_{\rm osp}})\right)$ queries to $\boldsymbol{\mathrm{U}}_{S}$,

• •

$O\left(\log(N)+b+\log(\frac{1}{\gamma})\right)$ qubits,

• •

$O\left(\frac{a^{2}b}{\gamma\epsilon_{\sigma}}\log(\frac{1}{\gamma})\log(\frac{a}{\epsilon_{\sigma}})\log(\frac{\log(a/\epsilon_{\sigma})}{1-2\delta_{\rm osp}})\right)$ other elementary gates.

To boost the success probability we can perform multiple runs and keep the median (often referred to as the “median trick”). Note that the returned value $\widetilde{\sigma}^{2}$ is always in $[0,1]$, regardless whether the algorithm succeeds or not. Now, assume that we execute $r$ runs, and consider the following random variables for each run:

[TABLE]

where $i=1,\ldots,r$. $X_{i}$ are independent Bernoulli variables which succeed with probability $p$, where $p\geq 3/4-\delta_{\rm osp}/2$. Let $X=\sum_{i=1}^{r}X_{i}$. We need to bound the probability that $X>r/2$. It holds that $\mathbb{E}[X]=rp$. Since, by assumption, $\delta_{\rm osp}<1/2,$ we can use a multiplicative Chernoff bound (see e.g. [62]), which gives $\Pr[X<(1-\eta)rp]<\exp(-\eta rp).$ For $\eta=1-\frac{1}{2p}$, we have that

[TABLE]

Setting $r=O\left(\frac{\log(1/\delta)}{1-2\delta_{\rm osp}}\right)$ we obtain the desired success probability $1-\delta$. ∎

2.3 Counting eigenvalues

In this section we describe the second subroutine of our algorithm, which is used to count the eigenvalues of a matrix that are smaller than a given threshold. This is a quantum variant of a standard subroutine in spectral divide-and-conquer eigensolvers [63, 64, 2, 3, 21, 19], often referred to as “NEGCOUNT.” Formally, we have the following lemma.

Lemma 2.3** ($\mathsf{QCOUNT}$).**

Let $\boldsymbol{\mathrm{A}}$ be a Hermitian matrix with $\|\boldsymbol{\mathrm{A}}\|\leq 1$, and a promise that $\sigma_{\min}(\boldsymbol{\mathrm{A}})\geq\Delta$, for some $\Delta>0$. Given $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}}$, an $(a,b,\epsilon_{\rm enc})$-QUBE for $\boldsymbol{\mathrm{A}}$, if $\epsilon_{\rm enc}\leq\tfrac{\Delta^{2}}{32a^{2}C_{\rm sgn}}N$ for some absolute constant $C_{\rm sgn}$, then there exists a quantum algorithm which returns $z\leftarrow\mathsf{QCOUNT}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}},\Delta,\delta)$, the true number of eigenvalues of $\boldsymbol{\mathrm{A}}$ that are smaller than zero, with probability at least $1-\delta$, with a cost:

Queries to $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}}$: $O\left(\tfrac{a^{2}N}{\Delta}\log\left(N\right)\log(\tfrac{1}{\delta})\right)$, 2. 2.

Number of qubits: $O\left(\log(N)+b\right)$, 3. 3.

Other elementary gates: $O\left(b\tfrac{a^{2}N}{\Delta}\log\left(N\right)\log(\tfrac{1}{\delta})\right)$, 4. 4.

Total cost: $O\left(\left(\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}};\epsilon_{\rm enc})+b\right)\cdot\tfrac{a^{2}N}{\Delta}\log\left(N\right)\log(\tfrac{1}{\delta})\right)$.

Proof.

The following steps achieve the result.

We use Corollary A.1 to prepare a block-encoding for $\operatorname{sgn}(\boldsymbol{\mathrm{A}})$. As long as $\epsilon_{\rm enc}\leq\epsilon_{\rm sgn}\tfrac{\Delta^{2}}{C_{\rm sgn}}a^{2}$ for some absolute constant $C_{\rm sgn}$, the corollary returns an $(1,b+2,\epsilon_{\rm sgn})$-QUBE $\boldsymbol{\mathrm{U}}_{\rm sgn}$ for $\operatorname{sgn}(\boldsymbol{\mathrm{A}})$ with cost

[TABLE] 2. 2.

Thereafter, Corollary A.2 is applied to compute the trace of $\operatorname{sgn}(\boldsymbol{\mathrm{A}})$ with accuracy $\epsilon_{\rm tr}$, which takes time

[TABLE] 3. 3.

We now set $\epsilon_{\rm sgn}=1/(8N)$, $\epsilon_{\rm tr}=1/4$, and $\epsilon_{\rm enc}=\epsilon_{\rm sgn}\tfrac{\Delta^{2}}{C_{\rm sgn}a^{2}}=\tfrac{\Delta^{2}}{8a^{2}CN}$. Note that this value of $\epsilon_{\rm enc}$ satisfies all the requirements of the previous steps. Plugging-in the parameters, the total complexity becomes

[TABLE]

It remains to verify the correctness. Since $\epsilon_{\rm enc}<\Delta/2$, the number of eigenvalues of $\boldsymbol{\mathrm{A}}$ (and also $\widetilde{\boldsymbol{\mathrm{A}}}$) that are smaller than zero, is equal to $z=\frac{1}{2}\tr(\boldsymbol{\mathrm{I}}-\operatorname{sgn}({\widetilde{\boldsymbol{\mathrm{A}}}}))=\frac{N}{2}-\frac{1}{2}\tr(\operatorname{sgn}(\widetilde{\boldsymbol{\mathrm{A}}}))=\frac{N}{2}-\frac{1}{2}\tr(\operatorname{sgn}(\boldsymbol{\mathrm{A}}))$.

From Corollary A.2, the above steps return a value $\widetilde{t}$ which satisfies $\left|\widetilde{t}-\tr(\operatorname{sgn}(\boldsymbol{\mathrm{A}}))\right|\leq 1/4$. This means that the value $\widetilde{z}=N/2-\widetilde{t}/2$ satisfies $|z-\widetilde{z}|\leq 1/8$. Therefore, we can round $\widetilde{z}$ to the nearest integer to recover the true $z$. ∎

3 Main algorithms

We now have all the required tools to describe our main algorithms for the $k$-th spectral gap and midpoint, namely, $\operatorname{\Delta}_{k}=\lambda_{k}-\lambda_{k+1}$ and $\mu_{k}=\frac{\lambda_{k}+\lambda_{k+1}}{2}$. Two algorithms are described. The first one, Algorithm 1 takes as input a Hermitian matrix $\boldsymbol{\mathrm{H}}$ (scaled to have $\|\boldsymbol{\mathrm{H}}\|\leq 1/2$), an integer $k\in[N-1]$, and a failure probability parameter $\delta$. It returns constant factor approximations to $\operatorname{\Delta}_{k}$ and $\mu_{k}$.

The second Algorithm 2, takes as input $\boldsymbol{\mathrm{H}}$, $k$, $\delta$, and an accuracy parameter $\epsilon$. It first uses Algorithm 1 as a subroutine to compute a coarse (constant factor) approximation to $\operatorname{\Delta}_{k}$ and $\mu_{k}$, and then uses these approximations to refine the accuracy of the solution to the given $\epsilon$. As a side-result, it also returns $\epsilon$-approximations to the actual eigenvalues $\lambda_{k}$ and $\lambda_{k+1}$ around the gap.

The complexity of these algorithms depends on the value of $\operatorname{\Delta}_{k}$. However, we highlight that no prior knowledge or lower bound on $\operatorname{\Delta}_{k}$ is required. The only required parameters are the ones given as input to the algorithms. All other parameters that are involved in the complexity are computed/approximated on-the-fly, if required.

Lemma 3.1** (Alg. 1, $\mathsf{QGAPconst}$).**

Let $\boldsymbol{\mathrm{H}}$ be an $N\times N$ Hermitian matrix with eigenvalues

[TABLE]

Let:

• •

$\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon_{\rm enc})$* be the cost of the unitary $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}$ which forms an $(a,b,\epsilon_{\rm enc})$-QUBE of $\boldsymbol{\mathrm{H}}$.*

• •

$\boldsymbol{\mathrm{\Pi}}$* be a $(\gamma,\delta_{\rm osp})$-OSP distribution, where $\boldsymbol{\mathrm{U}}_{S}\sim\boldsymbol{\mathrm{\Pi}}$ has cost $\mathcal{C}(\boldsymbol{\mathrm{U}}_{S})$.*

Then, given $\boldsymbol{\mathrm{H}}$ and an integer $k\in[N-1]$, Algorithm 1 ($\mathsf{QGAPconst}(\boldsymbol{\mathrm{H}},k,\delta)$) computes an approximate midpoint $\widetilde{\mu}_{k}$ and gap $\widetilde{\operatorname{\Delta}}_{k}$, such that:

[TABLE]

The algorithm succeeds with probability at least $1-\delta$. It makes at most:

Queries to unitaries with gate complexity $O\left(\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon_{\rm enc})\right)$, where $\epsilon_{\rm enc}\in\Theta\left(\frac{\operatorname{\Delta}_{k}^{2}}{a}\right)$:

[TABLE] 2. 2.

Queries to $\boldsymbol{\mathrm{U}}_{S}$:

[TABLE] 3. 3.

Queries to unitaries with gate complexity $O\left(\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon^{\prime}_{\rm enc})\right)$, where $\epsilon_{\rm enc}^{\prime}\in\Theta\left(\frac{\operatorname{\Delta}_{k}^{2}}{a^{2}N}\right)$:

[TABLE] 4. 4.

Number of qubits: $O(\log(N)+b+\log(1/\gamma)),$ 5. 5.

Other elementary gates: $O\left(b\cdot(Q_{1}+Q_{2})\right)$.

Proof.

The algorithm executes the following steps. Initially, it sets $j=2$, and $\operatorname{\Delta}=1/2$. At each iteration $j=2,\ldots$, the goal is to search for all gaps that are larger than $\operatorname{\Delta}^{j}=\frac{1}{2^{j}}$.

In the first step of each iteration, a uniform grid with points $h_{i}=-\frac{1}{2}+i\frac{\operatorname{\Delta}^{j}}{2}$, $i=0,\ldots,\frac{2}{\operatorname{\Delta}^{j}}$, such that $|h_{i}-h_{i+1}|=\operatorname{\Delta}^{j}/2$. If there exists any gap with width at least $\operatorname{\Delta}^{j}$, then there exists at least one point $h_{i}$ such that $\sigma_{\min}(\boldsymbol{\mathrm{H}}-h_{i})\geq\operatorname{\Delta}^{j}/4$. For each $h_{i}$ in the grid, the algorithm proceeds as follows:

Use Lemma A.1 to create an $(a+|h_{i}|,b+1,\epsilon_{\rm enc})$-QUBE $\boldsymbol{\mathrm{U}}_{i}$ of $\boldsymbol{\mathrm{H}}-h_{i}$, where we set $\epsilon_{\rm enc}:=(\frac{\operatorname{\Delta}^{j}}{8})^{2}\frac{1}{4a}$. The cost of $\boldsymbol{\mathrm{U}}_{i}$ is $\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon_{\rm enc})+O(1)$. 2. 2.

Execute $\mathsf{QSMIN}$ with parameters $\epsilon_{\sigma}=\frac{\operatorname{\Delta}^{2j}}{64}$ and $\delta_{\sigma}=\frac{\delta_{j}}{2}$ to test if $\sigma_{\min}(\boldsymbol{\mathrm{H}}-h_{i})$ satisfies the desired bound.

The requirement $\epsilon_{\rm enc}\leq\frac{\epsilon_{\sigma}}{4a}$ of Lemma 2.2 is satisfied, and therefore $\mathsf{QSMIN}$ returns a value $\widetilde{\sigma}_{i}^{2}$ which satisfies

[TABLE]

This means that, for each $h_{i}$, the algorithm can test if the computed $\widetilde{\sigma}_{i}^{2}$ satisfies $\widetilde{\sigma}_{i}^{2}\geq\operatorname{\Delta}^{2j}/32$, and, if this holds for any $i$, then it is ensured that $\sigma^{2}_{\min}(\boldsymbol{\mathrm{H}}-h_{i})\geq\widetilde{\sigma}^{2}_{i}-\epsilon_{\sigma}>\frac{\operatorname{\Delta}^{2j}}{64}$, or, in other words, $\sigma_{\min}(\boldsymbol{\mathrm{H}}-h_{i})>\frac{\operatorname{\Delta}^{j}}{8}$.

Now, this lower bound allows us to count the eigenvalues smaller than $h_{i}$ efficiently using $\mathsf{QCOUNT}$:

Use Lemma A.1 to construct an $(a+|h_{i}|,b+1,\epsilon^{\prime}_{\rm enc})$-QUBE $\boldsymbol{\mathrm{U}}^{\prime}_{i}$ for $\boldsymbol{\mathrm{H}}-h_{i}$ with $\epsilon^{\prime}_{\rm enc}:=\tfrac{\Delta^{2j}}{64\cdot 32C_{\rm sgn}a^{2}N}$. Here $C_{\rm sgn}$ is the constant from Lemma 2.3 (note that this value of $\epsilon^{\prime}_{\rm enc}$ satisfies the requirement of the lemma). The cost of $\boldsymbol{\mathrm{U}}^{\prime}_{i}$ is $\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon^{\prime}_{\rm enc})+O(1)$. 2. 4.

By executing $\mathsf{QCOUNT}(\boldsymbol{\mathrm{U}}^{\prime}_{i},\frac{\operatorname{\Delta}^{j}}{8},\frac{\delta_{j}}{2})$, we obtain $z_{i}$, the true number of eigenvalues that are smaller that $h_{i}$, with probability at least $1-\delta_{j}/2$.

We can now calculate the total costs for each iteration. By replacing $\epsilon_{\sigma}=\frac{\operatorname{\Delta}^{2j}}{64}$ and $\delta_{\sigma}=\delta_{j}/2=\frac{\delta\operatorname{\Delta}^{j}}{2^{j+2}}$ in Lemma 2.2, each call to $\mathsf{QSMIN}$ costs:

Queries to $\boldsymbol{\mathrm{U}}_{i}$:

[TABLE] 2. 2.

Queries to random state-preparation circuits $\boldsymbol{\mathrm{U}}_{S}$:

[TABLE] 3. 3.

Number of qubits:

[TABLE] 4. 4.

Other elementary gates:

[TABLE]

From Lemma 2.3, each call to $\mathsf{QCOUNT}\left(\boldsymbol{\mathrm{U}}^{\prime}_{i},\frac{\Delta^{j}}{8},\frac{\delta\Delta^{j}}{2^{j+1}}\right)$ costs:

Queries to $\boldsymbol{\mathrm{U}}^{\prime}_{i}$: $O\left(\tfrac{a^{2}N}{\operatorname{\Delta}^{j}}\log\left(N\right)\log(\tfrac{2^{j}}{\operatorname{\Delta}^{j}\delta})\right)$, 2. 2.

Number of qubits: $O(\log(N)+b)$, 3. 3.

Other elementary gates: $O\left(b\tfrac{a^{2}N}{\operatorname{\Delta}^{j}}\log\left(N\right)\log(\tfrac{2^{j}}{\operatorname{\Delta}^{j}\delta})\right)$, 4. 4.

Total cost: $O\left(\left(\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon^{\prime}_{\rm enc})+b\right)\cdot\tfrac{a^{2}N}{\operatorname{\Delta}^{j}}\log\left(N\right)\log(\tfrac{2^{j}}{\operatorname{\Delta}^{j}\delta})\right)$.

The number of iterations can be bounded as follows. First, note that a sufficient termination criterion is $\operatorname{\Delta}^{j}\leq\operatorname{\Delta}_{k}$. If this holds, there exists at least one point $h_{i}$ in the grid that satisfies $\sigma_{\min}(\boldsymbol{\mathrm{H}}-h_{i})\geq\frac{\Delta_{k}}{4}>\frac{\Delta^{j}}{4}$. Subsequently, the approximate singular value returned by $\mathsf{QSMIN}$ will satisfy $\widetilde{\sigma}_{i}^{2}\geq\sigma_{i}^{2}-\operatorname{\Delta}^{2j}/64>\frac{1.5\operatorname{\Delta}^{2j}}{32}$, which means that the algorithm will certainly execute $\mathsf{QCOUNT}$ for that point $h_{i}$. Ultimately, $\mathsf{QCOUNT}$ will verify that $h_{i}$ is inside the correct gap, and the algorithm will terminate.

We can obtain a necessary termination criterion $\frac{\operatorname{\Delta}^{j}}{4}\leq\operatorname{\Delta}_{k}$ in a similar way. To see this, assume that the algorithm terminated after finding a point that satisfies $\widetilde{\sigma}_{i}^{2}>\frac{\operatorname{\Delta}^{2j}}{32}$. This implies that $\sigma_{i}^{2}>\frac{\operatorname{\Delta}^{2j}}{64}\implies\sigma_{i}>\frac{\operatorname{\Delta}^{j}}{8}$. But since $h_{i}$ is inside the correct $k$-th gap, it means that $\operatorname{\Delta}_{k}\geq 2\sigma_{i}\geq\frac{\operatorname{\Delta}^{j}}{4}$, as initially claimed.

Since we are halving $\operatorname{\Delta}^{j}$ in every iteration, we conclude that the algorithm terminates after $m:=\lceil\log(1/\operatorname{\Delta}_{k})\rceil+C_{\rm gap}$ iterations, where $C_{\rm gap}\in[0,2]$. This means that $\operatorname{\Delta}^{m}=(1/2)^{m}\in\Theta(\operatorname{\Delta}_{k})$. The most expensive iteration is the last one, where the corresponding block-encodings have costs:

• •

$Q:=\mathcal{C}(\boldsymbol{\mathrm{U}}_{i};\epsilon_{\rm enc})\in\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}},\epsilon_{\rm enc})+O(1)$, where

[TABLE]

• •

$Q^{\prime}:=\mathcal{C}(\boldsymbol{\mathrm{U}}^{\prime}_{i};\epsilon^{\prime}_{\rm enc})\in\mathcal{C}(\boldsymbol{\mathrm{U}}^{\prime}_{\boldsymbol{\mathrm{H}}},\epsilon^{\prime}_{\rm enc})+O(1)$, where

[TABLE]

The total cost of the last iteration is therefore:

[TABLE]

Let us bound each additive term separately. By replacing $\operatorname{\Delta}=1/2$ we have:

[TABLE]

For the second term we have

[TABLE]

After simplifying some constants, the final total cost is bounded by the sum of the following terms:

• •

$O\left(\left(\mathcal{C}(\boldsymbol{\mathrm{U}}_{i};\epsilon_{\rm enc})+b\right)\frac{a^{2}}{\gamma\operatorname{\Delta}_{k}^{2}}\log(\frac{1}{\gamma})\log(\frac{a}{\operatorname{\Delta}_{k}})\log(\frac{\log(a/\operatorname{\Delta}_{k})}{1-2\delta_{\rm osp}})\frac{\log(1/\delta)}{1-2\delta_{\rm osp}}\log(\frac{1}{\operatorname{\Delta}_{k}})\right)$,

• •

$O\left(\mathcal{C}(\boldsymbol{\mathrm{U}}_{S})\frac{1}{\gamma}\log(\frac{a}{\operatorname{\Delta}_{k}})\log(\frac{\log(a/\operatorname{\Delta}_{k})}{1-2\delta_{\rm osp}})\frac{\log(1/\delta)}{1-2\delta_{\rm osp}}\log(\frac{1}{\operatorname{\Delta}_{k}})\right)$,

• •

$O\left(\left(\mathcal{C}(\boldsymbol{\mathrm{U}}^{\prime}_{i};\epsilon^{\prime}_{\rm enc})+b\right)\frac{a^{2}N}{\operatorname{\Delta}_{k}}\log\left(N\right)\log(\frac{1}{\operatorname{\Delta}_{k}\delta})\log(\frac{1}{\operatorname{\Delta}_{k}})\right),$

where $\mathcal{C}(\boldsymbol{\mathrm{U}}_{i};\epsilon_{\rm enc})$ and $\mathcal{C}(\boldsymbol{\mathrm{U}}^{\prime}_{i};\epsilon^{\prime}_{\rm enc})$ depend on the values of $\epsilon_{\rm enc}$ and $\epsilon_{\rm enc}^{\prime}$, respectively, in the final iteration, and the trailing $\log(\tfrac{1}{\operatorname{\Delta}_{k}})$ factor is the number of iterations.

Next, we argue about the correctness of the algorithm and the approximation guarantees of the final results. Once we find some $h_{i}$, for which the result of $\mathsf{QCOUNT}$ satisfies $z_{i}=k$, we know that the point $h_{i}$ is inside the correct gap $\operatorname{\Delta}_{k}$. The algorithm then returns $\widetilde{\mu}_{k}\leftarrow h_{i}$ and $\widetilde{\operatorname{\Delta}}_{k}\leftarrow\Delta^{j}$, which satisfy the following bounds. In the above, we already argued that $\operatorname{\Delta}_{k}\geq\operatorname{\Delta}^{j}/4$ is necessary for the algorithm to terminate. But we can also see that $\operatorname{\Delta}^{j}\geq\operatorname{\Delta}_{k}/2$ holds as well. Indeed, the algorithm cannot terminate when $\operatorname{\Delta}^{j}<\operatorname{\Delta}_{k}/2$, otherwise it would have terminated in the previous step. The combination of these two facts implies that the returned $\operatorname{\Delta}^{m}$ satisfies:

[TABLE]

Combining this with the termination criterion of the algorithm, we also conclude that the returned midpoint $\widetilde{\mu}_{k}$ satisfies:

[TABLE]

where in the last we replaced $\lambda_{k+1}=\mu_{k}-\frac{\operatorname{\Delta}_{k}}{2}$ and $\lambda_{k}=\mu_{k}+\frac{\operatorname{\Delta}_{k}}{2}$.

Finally, we bound the success probability. There are a total of at most $l_{j}:=1+2/\operatorname{\Delta}^{j}$ calls to $\mathsf{QCOUNT}$ and $\mathsf{QSMIN}$ in each iteration $j$. Each call fails with probability at most $p_{j}:=\delta_{j}/2$. Thus, the total success probability for iteration $j$ is at least

[TABLE]

Taking a union bound over all iterations, we have

[TABLE]

∎

With some additional steps the solution accuracy can be increased to any $\epsilon>0$.

Theorem 3.1**.**

Let $\boldsymbol{\mathrm{H}}$ be an $N\times N$ Hermitian matrix, $\|\boldsymbol{\mathrm{H}}\|\leq 1/2$, $k\in[N-1]$ be an integer, and $\epsilon,\delta\in(0,1)$ an accuracy and failure probability parameter, respectively. Let

• •

$\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon_{\rm enc})$* denote the cost of the unitary $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}$ which forms an $(a,b,\epsilon_{\rm enc})$-QUBE of $\boldsymbol{\mathrm{H}}$.*

• •

$\boldsymbol{\mathrm{\Pi}}$* be a $(\gamma,\delta_{\rm osp})$-OSP distribution.*

Algorithm 2 returns $\widetilde{\lambda}_{k},\widetilde{\lambda}_{k+1},\widetilde{\mu}_{k},\widetilde{\operatorname{\Delta}}_{k}$ that satisfy:

[TABLE]

In addition to the costs induced by calling Algorithm 1 in the first step, it requires

Queries to unitaries with complexity $O(\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon_{\rm enc}))$, where $\epsilon_{\rm enc}\in\Theta(\operatorname{\Delta}_{k}^{2}/a)$:

[TABLE] 2. 2.

Queries to OSP circuits $\boldsymbol{\mathrm{U}}_{S}\sim\boldsymbol{\mathrm{\Pi}}$:

[TABLE] 3. 3.

Queries to unitaries with complexity $Q_{2}:=O\left(\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon^{\prime}_{\rm enc})\right)$, where $\epsilon^{\prime}_{\rm enc}\in\Theta(\epsilon^{2}\operatorname{\Delta}_{k}^{2}/(a^{2}N))$:

[TABLE] 4. 4.

Number of qubits: $O\left(\log(N)+b+\log(\frac{1}{\gamma})\right)$, 5. 5.

Other elementary gates: $O\left(b\cdot(Q_{1}+Q_{2})\right)$.

If we use Lemma 2.1 for $\boldsymbol{\mathrm{\Pi}}$, after some simplifications, the total circuit complexity (of both Alg. 1 and 2) is upper bounded by

[TABLE]

where $\epsilon_{1}\in\Theta\left(\frac{\epsilon^{2}\operatorname{\Delta}_{k}^{2}}{a^{2}N}\right)$.

Proof.

Let us call $\widehat{\mu}_{k}$ and $\widehat{\operatorname{\Delta}}_{k}$ the initial solutions returned by Algorithm 1, i.e., the solutions that satisfy the constant factor accuracy of Lemma 3.1. We know that $\widehat{\operatorname{\Delta}}_{k}\in[0.5\operatorname{\Delta}_{k},4\operatorname{\Delta}_{k}]$. By setting $\epsilon_{k}:=\epsilon\widehat{\operatorname{\Delta}}_{k}/4$, it means that $\epsilon_{k}\leq\epsilon\operatorname{\Delta}_{k}$.

We now approximate both eigenvalues $\lambda_{k}$ and $\lambda_{k+1}$ to accuracy $\epsilon_{k}/2$ each as follows. Consider first $\lambda_{k+1}$. We make a grid on the left of $\widehat{\mu}_{k}$ with step $\epsilon_{k}/4$, consisting of the points $l_{i}=\widehat{\mu}_{k}+i\epsilon_{k}/4$, for $i=-\lceil 2\widehat{\operatorname{\Delta}}_{k}/\epsilon_{k}\rceil,\ldots,-2,-1$. Note that, there exists at least one point $l^{*}>\lambda_{k+1}$ such that $\epsilon_{k}/4\leq|\lambda_{k+1}-l^{*}|<\epsilon_{k}/2,$ or, equivalently, $\sigma_{\min}(\boldsymbol{\mathrm{H}}-l^{*})\in[\epsilon_{k}/4,\epsilon_{k}/2]$.

To find such a point, we first set the parameters:

[TABLE]

and prepare two block-encodings:

• •

$\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}\leftarrow\mathsf{QENC}_{a,b}(\boldsymbol{\mathrm{H}},\epsilon_{\rm enc})$,

• •

$\boldsymbol{\mathrm{U}}^{\prime}_{\boldsymbol{\mathrm{H}}}\leftarrow\mathsf{QENC}_{a,b}(\boldsymbol{\mathrm{H}},\epsilon^{\prime}_{\rm enc})$.

Then we execute the steps 10-18 for each $i$, similar to Algorithm 1, but with the aforementioned parameters. We can verify that the returned values $\widetilde{\lambda}_{k},\widetilde{\lambda}_{k+1},\widetilde{\mu}_{k},\widetilde{\operatorname{\Delta}}_{k},$ satisfy the advertised bounds.

It remains to bound the complexity. There are a total of $O(\epsilon^{-1})$ points in the two grids, which gives the total number of calls to the subroutines $\mathsf{QSMIN}$ and $\mathsf{QCOUNT}$. At each grid point, we call $\mathsf{QENC}$ with parameter $\epsilon_{\rm enc}=3\epsilon_{k}^{2}/256a\in\Theta(\operatorname{\Delta}_{k}^{2}/a)$, followed by $\mathsf{QSMIN}$ with parameters:

[TABLE]

The cost for each $\mathsf{QSMIN}$ call is

$O\left(\frac{a^{2}}{\gamma\epsilon^{2}\operatorname{\Delta}_{k}^{2}}\log(\frac{1}{\gamma})\log(\frac{a}{\epsilon\operatorname{\Delta}_{k}})\log(\frac{\log(a/(\epsilon\operatorname{\Delta}_{k}))}{1-2\delta_{\rm osp}})\frac{\log(1/(\epsilon\delta))}{1-2\delta_{\rm osp}}\right)$ queries to $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}$, 2. 2.

$O\left(\frac{1}{\gamma}\log(\frac{a}{\epsilon\operatorname{\Delta}_{k}})\log(\frac{\log(a/(\epsilon\operatorname{\Delta}_{k}))}{1-2\delta_{\rm osp}})\frac{\log(1/(\epsilon\delta))}{1-2\delta_{\rm osp}}\right)$ queries to OSP circuits $\boldsymbol{\mathrm{U}}_{S}\sim\boldsymbol{\mathrm{\Pi}}$, 3. 3.

$O\left(\log(N)+b+\log(\frac{1}{\gamma})\right)$ qubits, 4. 4.

$O\left(\frac{a^{2}}{\gamma\epsilon^{2}\operatorname{\Delta}_{k}^{2}}\log(\frac{1}{\gamma})\log(\frac{a}{\epsilon\operatorname{\Delta}_{k}})\log(\frac{\log(a/(\epsilon\operatorname{\Delta}_{k}))}{1-2\delta_{\rm osp}})\frac{\log(1/(\epsilon\delta))}{1-2\delta_{\rm osp}}\right)$ other elementary gates.

In the worst case, $\mathsf{QCOUNT}$ will also be called in every iteration. Each call takes as input an $(a,b,\epsilon^{\prime}_{\rm enc})$–QUBE, accuracy $\epsilon^{\prime}_{\rm enc}=\epsilon^{2}_{k}/(64a^{2}C_{\rm sgn}N)\in\Theta(\epsilon^{2}\operatorname{\Delta}_{k}^{2}/(a^{2}N)),$ and failure probability $\delta_{\sigma}=\Theta(\epsilon\delta)$. Substituting in Lemma 2.3, each $\mathsf{QCOUNT}$ call costs:

[TABLE]

calls to the block-encoding circuit $\boldsymbol{\mathrm{U}}^{\prime}_{\boldsymbol{\mathrm{H}}}$. ∎

4 Complexity in QRAM

In this section we specialize the complexity analysis in the QRAM model [58]. We describe a specific technique to construct block-encodings via the so-called quantum-accessible data structures, studied in [41, 49, 50, 65]. We assume a dense $N\times N$ matrix $\boldsymbol{\mathrm{A}}\in\mathbb{C}$, whose elements are known (i.e. pre-computed), and stored as an array in a classical memory. Given this classical array, a data structure of the form of [41] can be constructed in an offline fashion and stored in a classical memory in $O(N^{2})$ arithmetic operations, and it requires $O(N^{2}\log(N))$ space, i.e. words, assuming each matrix element requires one word to store.

For dense matrices, the data structure consists of $N$ binary trees with $O(\log(N))$ depth and $N$ leaves each. If we were to prepare a quantum circuit for the access oracle of this data structure, we would need $O(N^{2})$ gates to implement it. To elaborate more, each binary tree $B_{i}$ can be obtained by constructing a circuit that rotates the state $\ket{0}^{\otimes\log(N)}$ to the state $\ket{B_{i}}$, where the latter contains the $N$ leaves of the tree $B_{i}$. Such a circuit can be constructed using $O(N)$ one and two qubit gates (rotations and CNOT gates) [66, 67], and we need such a circuit for every $i\in[N]$. However, if the quantum computer has the ability to query the elements of the classical data structure in superposition, which is possible in the QRAM model, the gate complexity can be reduced to $O(\operatorname{polylog}(N))$ [41]. In particular, Theorem 10 of [41] showed that this construction enables efficient circuits implementing the maps

[TABLE]

where $\ket{\boldsymbol{\mathrm{A}}_{i}}$ is the normalized quantum state corresponding to the $i$-th row of $\boldsymbol{\mathrm{A}}$ and $\braket{i|\widetilde{\boldsymbol{\mathrm{A}}}}=\norm{\boldsymbol{\mathrm{A}}_{i}}$ encodes the row norms. These circuits can be implemented to precision $\epsilon$ with gate complexity $O\left(\operatorname{polylog}(N/\epsilon)\right)$ [41, 50]. The following lemma can be used to construct a block-encoding from a quantum-accessible data structure in QRAM.

Lemma 4.1** (Imported result from [50, 49]).**

Let $\boldsymbol{\mathrm{A}}\in\mathbb{C}^{N\times N}$ be a matrix stored in a quantum-accessible data structure. There exists a $\mathsf{QENC}_{a,b}(\boldsymbol{\mathrm{A}},\epsilon)$ subroutine with $a=\|\boldsymbol{\mathrm{A}}\|_{F},b=\lceil\log(N)\rceil+1,$ and cost $\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}};\epsilon)=O(\operatorname{polylog}(N/\epsilon))$.

By using the specific block-encodings from Appendix A and state preparation circuits from Section 2, we can bound the complexity of the algorithm as follows.

Corollary 4.1**.**

Using Lemma˜4.1 to construct block-encodings, and Lemma˜2.1 for state preparation, the total complexity of Theorem 3.1 becomes:

[TABLE]

and the algorithm succeeds with probability at least $1-\delta$.

Proof.

It suffices to replace the costs of the block-encodings in Theorem 3.1. For all block-encodings we can use Lemma˜4.1, which means that $a=\|\boldsymbol{\mathrm{H}}\|_{F}\in O(\sqrt{N})$. The cost has a polylogarithmic dependence on the accuracy. Replacing the costs in Eq. (9) we have

[TABLE]

Combining everything together, Eq. (9) becomes

[TABLE]

where we absorbed all polylogarithmic factors in the trailing term. ∎

5 A lower bound in the Black-Box model

In this section we derive a lower bound for the problem of deciding the existence of a spectral gap of a binary matrix, in the so-called Black-Box access model. Here we assume a unitary $\boldsymbol{\mathrm{O}}_{\boldsymbol{\mathrm{A}}}$ through which we can query elements of the matrix $\boldsymbol{\mathrm{A}}$ via:

[TABLE]

Several lower bounds are known for related problems. A straightforward query lower bound for determining the midpoint $\mu_{k}$ can be obtained from the median problem [68]. Let $\mathcal{A}$ be a quantum algorithm which, given $\boldsymbol{\mathrm{O}}_{\boldsymbol{\mathrm{A}}}$, and an integer $k$, it returns any (classical) value $\widetilde{\mu}_{k}$ such that $\widetilde{\mu}_{k}\in(\lambda_{k},\lambda_{k+1})$. Then $\mathcal{A}$ can be directly used to solve the median problem, where, given oracle access to a sequence of $N$ numbers $\{x_{1},\ldots,x_{N}\}$, where $N$ is even, we seek to find a number that lies strictly inside $x\in(x_{N/2},x_{N/2}+1)$. Such a sequence can be seen as the elements of a diagonal matrix, and the eigenvalues of such a matrix are the diagonal elements themselves. [68] proved a $\Omega(N)$ query lower bound for the median problem, assuming that we have oracle access to a quantum circuit implementing the sequence. Indeed, the same bound holds even if we only want to verify that $\widetilde{\mu}_{k}$ is inside the correct interval, since this generalizes the counting problem [68, 69, 70].

Here we obtain a stronger lower bound for deciding whether there exists a gap between two consecutive eigenvalues, inspired by the reductions of [56, 57]. We remark that the construction below is inherently non-Hermitian. Whether analogous black-box lower bounds hold for Hermitian matrices remains an open question.

Theorem 5.1**.**

Given a matrix $\boldsymbol{\mathrm{A}}\in\{0,1\}^{N\times N}$, a unitary $\boldsymbol{\mathrm{O}}_{\boldsymbol{\mathrm{A}}}:\ket{i}\ket{j}\ket{0}\mapsto\ket{i}\ket{j}\ket{\boldsymbol{\mathrm{A}}_{i,j}}$, and an integer $k$, any quantum algorithm which answers TRUE if $\operatorname{\Delta}_{k}:=|\lambda_{k}(\boldsymbol{\mathrm{A}})-\lambda_{k+1}(\boldsymbol{\mathrm{A}})|>0$, and FALSE otherwise, must execute $\Omega(N^{2})$ queries to $\boldsymbol{\mathrm{O}}_{\boldsymbol{\mathrm{A}}}$. Here the eigenvalues of $\boldsymbol{\mathrm{A}}$ are order such that $Re(\lambda_{1})\geq Re(\lambda_{2})\geq\ldots\geq Re(\lambda_{N})$, and if $Re(\lambda_{k})=Re(\lambda_{k+1})$, then $Im(\lambda_{k})\geq Im(\lambda_{k+1}).$

Proof.

The proof closely follows the techniques of [56, 57] for other similar problems (rank, determinant, matrix power, and inverse), essentially by reducing the problem to counting the number of triangles in a directed graph. Assume that we are given a set $X\in\{0,1\}^{N^{2}}$ and the value $a=\lceil N^{2}/2\rceil$, and we ask whether $\sum_{i=1}^{N^{2}}X_{i}=a$. The query complexity of this problem is $\Omega(N^{2})$ [56].

For the reduction the first step is to create a couple of graphs. We create a bipartite graph $G$, with $N$ sources $s_{i}$ and $N$ targets $t_{j}$, and we insert an edge between $s_{i}$ and $t_{j}$ iff $X_{iN+j}=1$. We also create a utility node $v$, which has $N$ outgoing edges to every one of the sources, and $N$ incoming edges from every one of the targets.

Thereafter, we construct a second bipartite graph $G^{\prime}$ which also has $N$ sources $s^{\prime}_{i}$, $N$ targets $t^{\prime}_{j}$, and a utility node $v^{\prime}$. We create exactly $a$ edges between the sources and the targets, in a round-robin manner, so that the edges are “evenly distributed” among the two vertex sets (for each $s_{i}$, we add roughly $\lceil a/N\rceil$ outgoing edges). Finally, we add $N$ edges from $v^{\prime}$ to every $s^{\prime}_{i}$, and $N$ edges from every $t^{\prime}_{j}$ to $v^{\prime}$.

Note that the two graphs are disjoint. Consider the larger graph $F=(V\cup V^{\prime},E\cup E^{\prime})$, and its adjacency matrix $\boldsymbol{\mathrm{A}}$. $F$ has $4N+2$ nodes. If we choose the node labels of $F$ appropriately, $\boldsymbol{\mathrm{A}}$ consists of two diagonal blocks: one block which corresponds to the adjacency of $G$, and the other to the one of $G^{\prime}$. More specifically,

[TABLE]

We will show that $\boldsymbol{\mathrm{A}}_{G}$ and $\boldsymbol{\mathrm{A}}_{G^{\prime}}$ have only three nonzero eigenvalues.

The elements of the matrix $\boldsymbol{\mathrm{A}}^{3}$ contain the number of paths of length three starting from every node $i$ to every other node $j$. Since the only length-three paths from $v$ and $v^{\prime}$ are towards themselves, the matrix $\boldsymbol{\mathrm{A}}^{3}$ has the following structure:

[TABLE]

where $a$ is located in position $\boldsymbol{\mathrm{A}}^{3}_{v^{\prime},v^{\prime}}$ and $\sigma:=\sum_{i=1}^{N^{2}}X_{i}$ in position $\boldsymbol{\mathrm{A}}^{3}_{v,v}$. It is not hard to see that $\sigma$ and $a$ are eigenvalues, since $\det(\sigma\boldsymbol{\mathrm{I}}-\boldsymbol{\mathrm{A}}^{3})=\det(a\boldsymbol{\mathrm{I}}-\boldsymbol{\mathrm{A}}^{3})=0$.

The matrices $\boldsymbol{\mathrm{B}}$ and $\boldsymbol{\mathrm{B}}^{\prime}$ have a special structure. Let us analyze $\boldsymbol{\mathrm{B}}$ first. From each source vertex $s_{i},i\in[N]$, there are exactly $d_{i}=outdeg(s_{i})$ paths of length three to every source $s_{j},$ including itself. This means that the row $i$ of $\boldsymbol{\mathrm{B}}$ contains exactly $N$ non-zeros, with the same value $d_{i}$, at the positions corresponding to the sources. This gives rise to a submatrix $\boldsymbol{\mathrm{B}}_{s}$. Looking at the targets $t_{l}$, note that there are exactly $f_{l}=indeg(t_{l})$ paths from every node $t_{m},m\in[N]$, to $t_{l}$. This means that the $l$-th column of the matrix $\boldsymbol{\mathrm{B}}$ has exactly $N$ nonzeros, equal to $f_{l}$, at the positions that correspond to target vertices. This gives rise to a submatrix $\boldsymbol{\mathrm{B}}_{t}$, where $\boldsymbol{\mathrm{B}}_{s}$ and $\boldsymbol{\mathrm{B}}_{t}$ are disjoint. With this description, we conclude that the matrix $\boldsymbol{\mathrm{B}}$ has the following structure:

[TABLE]

where both $\boldsymbol{\mathrm{B}}_{s}$ and $\boldsymbol{\mathrm{B}}_{t}$ have size $N\times N$. It turns out that they are rank-one matrices. Indeed, we can write $\boldsymbol{\mathrm{B}}_{s}=\boldsymbol{\mathrm{d}}_{s}\boldsymbol{\mathrm{e}}^{\top}$, where $\boldsymbol{\mathrm{e}}$ is the all-ones vector, $\boldsymbol{\mathrm{d}}_{s}^{\top}=\begin{pmatrix}d_{1}&\ldots&d_{N}\end{pmatrix}$, and $\boldsymbol{\mathrm{B}}_{t}=\boldsymbol{\mathrm{e}}\boldsymbol{\mathrm{f}}_{t}^{\top}$, where $\boldsymbol{\mathrm{f}}_{t}^{\top}=\begin{pmatrix}f_{1}&\ldots&f_{N}\end{pmatrix}$. The only nonzero eigenvalues of $\boldsymbol{\mathrm{B}}_{s}$ and $\boldsymbol{\mathrm{B}}_{t}$ are equal to $\lambda_{s}=\boldsymbol{\mathrm{d}}_{s}^{\top}\boldsymbol{\mathrm{e}}=\sum_{i=1}^{N}d_{i}$ and $\lambda_{f}=\boldsymbol{\mathrm{f}}_{t}^{\top}\boldsymbol{\mathrm{e}}=\sum_{i=1}^{N}{f_{i}}$. We can see that

[TABLE]

Next, we shift our attention to $\boldsymbol{\mathrm{B}}^{\prime}$. Since it has very similar structure to $\boldsymbol{\mathrm{B}}$, we can argue accordingly that the only two non-zero eigenvalues of $\boldsymbol{\mathrm{B}}^{\prime}$ are equal to $a$. Now, it follows that

[TABLE]

The above imply that $\boldsymbol{\mathrm{A}}^{3}$ is diagonalizable. The eigenvalues of $\boldsymbol{\mathrm{A}}$ are the third-roots of the eigenvalues of $\boldsymbol{\mathrm{A}}^{3}$. Since $\boldsymbol{\mathrm{A}}$ is real, its complex eigenvalues come in conjugate pairs. The spectrum of $\boldsymbol{\mathrm{A}}$ is equal to

[TABLE]

where $\{a_{1},a_{2},a_{3}\}$ are either all equal to $|a^{1/3}|$ or equal to $|a^{1/3}|\times\{z_{1},z_{2},z_{3}\}$ where $z_{1},z_{2},z_{3}$ are the three third-roots or unity. Similarly, $\{\sigma_{1},\sigma_{2},\sigma_{3}\}$ are either all equal to $|\sigma^{1/3}|$ or $|\sigma^{1/3}|\times\{z_{1},z_{2},z_{3}\}$.

Let $\{\lambda_{1},\lambda_{2},\ldots,\lambda_{6}\}$ be the eigenvalues of $\boldsymbol{\mathrm{A}}$ sorted in descending order of their real parts, where conflicts are resolved with respect to the imaginary part, i.e.

[TABLE]

We now derive the final argument.

$(\implies):$ If $a=\sigma$, we have the following four cases for the eigenvalues of $\boldsymbol{\mathrm{A}}$:

[TABLE]

$(\impliedby):$ On the opposite direction, we prove each case separately:

• •

If $\lambda_{i}=\lambda_{i+1}$, for $i=1,\ldots,5$, all eigenvalues are equal and hence $a=\sigma$.

• •

$\lambda_{1}=\lambda_{2}=\lambda_{3}=\lambda_{4}\neq\lambda_{5}\neq\lambda_{6}$. We already argued that $a_{i}$ are either all equal to each other or equal to the scaled roots of unity. The same for $\sigma_{i}$. The only way that the top-three $\lambda_{i}$ are equal is if either all $a_{i}$ are equal and real, or all $\sigma_{i}$. Since they are real, and equal to $\lambda_{4}$, it means that $\lambda_{4}$ is also real. We conclude that at least one of $\lambda_{1},\ldots,\lambda_{4}$ is equal to $a$, and at least one of them is equal to $\sigma$, therefore $a=\sigma$.

• •

The only way that top-6 eigenvalues are equal in pairs, is if both $a_{i}$ and $\sigma_{i}$ are the scaled roots of unity, and, in this case, with the same scaling factor, which also implies that $a=\sigma$.

We conclude that $a=\sigma$ if and only if one of the following relation holds:

• •

$\lambda_{i}=\lambda_{i+1},i=1,\ldots,5.$

• •

$\lambda_{1}=\lambda_{2}=\lambda_{3}=\lambda_{4}\neq\lambda_{5}\neq\lambda_{6}$.

• •

$\lambda_{1}=\lambda_{2}\neq\lambda_{3}=\lambda_{4}\neq\lambda_{5}=\lambda_{6}$.

Therefore, one can call the gap algorithm five times, for $k=1,\ldots,5$, and revisit the three cases above in order to verify whether $a=\sigma$. ∎

6 Conclusion

In this work we studied algorithms for computing spectral gaps between eigenvalues of Hermitian matrices and the corresponding midpoints. One of the main subroutines is a randomized state preparation method, which achieves at least $1/N$ overlap with any fixed target state with probability at least $3/5$, and requires $O(N)$ classical random bits to prepare. This can be used to efficiently approximate the smallest singular value of a block-encoded matrix with high probability, and to count the eigenvalues of a Hermitian matrix that are smaller than a given threshold. In the quantum-accessible data structure model, the total query complexity of our main algorithm scales at most quadratically with respect to the problem size $N$. This indicates the potential for a polynomial quantum speed-up, since the complexity of the best-known classical algorithms require $\widetilde{O}(N^{\omega})$ arithmetic operations. In the black-box access model, we also proved an $\Omega(N^{2})$ query lower bound for deciding whether there is a gap between two consecutive eigenvalues of a binary matrix.

\c@NAT@ctr

Appendix A Preliminaries on block-encodings and eigenvalue transformation

In this section we recall some basic definitions and known results for block-encodings. Block-encoding refers to the derivation of a unitary operator $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}}$ which can be described by a quantum circuit and “encodes” a given matrix $\boldsymbol{\mathrm{A}}$, which was stated in Definition 1.1.

Constructing block-encodings with low cost is not a trivial task. A common approach in the literature is to assume a “black-box” unitary oracle that gives access to a (sparse) matrix $\boldsymbol{\mathrm{A}}$. See eg. [49, 71]. Recent works take a step further, describing explicit quantum circuits for the block-encoding of special classes of matrices [72, 73, 50, 74]. In the main algorithm we assume the existence of a general subroutine that can prepare block-encodings, defined as follows.

Definition A.1**.**

A subroutine $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}}\leftarrow\mathsf{QENC}_{a,b}(\boldsymbol{\mathrm{A}},\epsilon)$ takes as input a matrix $\boldsymbol{\mathrm{A}}$ of size $N$ and an accuracy $\epsilon$, and prepares an $(a,b,\epsilon)$-QUBE of $\boldsymbol{\mathrm{A}}$. $a$ and $b$ might be functions of $N$ or $\|\boldsymbol{\mathrm{A}}\|_{\xi}$, for some norm $\xi$.

It is straightforward to obtain a block-encoding after a diagonal shift.

Lemma A.1** ($\mathsf{QSHIFT}$, follows from [49]).**

Let $h\in[-1/2,1/2]$ be a scalar and $\boldsymbol{\mathrm{H}}$ be a Hermitian matrix with $\|\boldsymbol{\mathrm{H}}\|\leq 1/2$, and $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}$ be an $(a,b,\epsilon_{\rm enc})$-QUBE for $\boldsymbol{\mathrm{H}}$ with cost $\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon_{\rm enc})$. We can construct a $(a+|h|,b+1,\epsilon_{\rm enc})$-QUBE for $\boldsymbol{\mathrm{H}}-h$ with a subroutine

[TABLE]

where $\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}-h};\epsilon_{\rm enc})=\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon_{\rm enc})+O(1)$.

Proof.

We use [49, Lemma 29] in a similar way as in [38]. First we construct a couple of $2$-by-$2$ unitaries $\boldsymbol{\mathrm{P}}_{L}\ket{0}=c_{1}\ket{0}+c_{2}\ket{1}$ and $\boldsymbol{\mathrm{P}}_{R}\ket{0}=d_{1}\ket{0}+d_{2}\ket{1}$, as follows. Set $\beta=1+|h|/a$. The values $c_{1},d_{1},c_{2},d_{2}$ are chosen such that $\beta(c_{1}^{*}d_{1})=1$ and $\beta(c_{2}^{*}d_{2})=-h/a$. Let $\boldsymbol{\mathrm{W}}=\ket{0}\bra{0}\otimes\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}+\ket{1}\bra{1}\otimes\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}$. One can verify that $(\boldsymbol{\mathrm{P}}_{L}^{\dagger}\otimes\boldsymbol{\mathrm{I}})\boldsymbol{\mathrm{W}}(\boldsymbol{\mathrm{P}}_{R}\otimes\boldsymbol{\mathrm{I}})$ is an $(a+|h|,b+1,\epsilon_{\rm enc})$-QUBE of $\boldsymbol{\mathrm{H}}-h$. ∎

A.1 Quantum eigenvalue transformation and matrix sign function

We can now describe how to prepare a block-encoding of the matrix sign function $\operatorname{sgn}(\boldsymbol{\mathrm{H}}-h)$ using Quantum Eigenvalue Transformation (QET). We follow the methodology of [38] which is convenient for our analysis, but we also note that many works have studied similar approximations [71, 49, 75, 76]. E.g., an elegant an practical approach is mentioned [75], using a recursive QET based on Padé approximations for the sign function, with sharp complexity bounds.

Corollary A.1** (Adapted from [38, 49]).**

Let $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}$ be an $(a,b,\epsilon_{\rm enc})$-QUBE of a Hermitian matrix $\boldsymbol{\mathrm{H}}$, $\|\boldsymbol{\mathrm{H}}\|\leq 1/2$, with gate complexity $\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon_{\rm enc})$. Let $h$ be a scalar in $[-1/2,1/2]$, and assume that $\Lambda\left(\boldsymbol{\mathrm{H}}-h\right)\subseteq[-1,-\Delta]\cup[\Delta,1]$, for some (known) $\Delta\in(0,1)$. Let $\epsilon_{\rm sgn}\in(0,1)$ be an input accuracy parameter. If $\epsilon_{\rm enc}\leq\epsilon_{\rm sgn}\frac{\Delta^{2}}{C_{\rm sgn}(a+|h|)^{2}}$, for some absolute constant $C_{\rm sgn}>1$, then we can construct a unitary $\boldsymbol{\mathrm{U}}_{\rm sgn}$ which is an $(1,b+2,\epsilon_{\rm sgn})$-QUBE of $\operatorname{sgn}(\boldsymbol{\mathrm{H}}-h)$. The cost of $\boldsymbol{\mathrm{U}}_{\rm sgn}$ is:

Queries to $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}-h}$: $O\left(\frac{a+|h|}{\Delta}\log(\frac{1}{\epsilon_{\rm sgn}})\right)$, with $\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}-h};\epsilon_{\rm enc})=\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon_{\rm enc})$, 2. 2.

Number of qubits: $\log(N)+b+2$, 3. 3.

Other elementary gates: $O\left(\frac{a+|h|}{\Delta}b\log(\frac{1}{\epsilon_{\rm sgn}})\right)$ 4. 4.

Total complexity: $C({\boldsymbol{\mathrm{U}}_{\rm sgn};\Delta,\epsilon)}=O\left((\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}};\epsilon_{\rm enc})+b)\cdot\frac{a+|h|}{\Delta}\log(\frac{1}{\epsilon_{\rm sgn}})\right)$.

Proof.

Since $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}}$ is an $(a,b,\epsilon_{\rm enc})$-QUBE for $\boldsymbol{\mathrm{H}}$, then, by definition, it is also an $(1,b,0)$-QUBE for a matrix $\widetilde{\boldsymbol{\mathrm{H}}}$, where $\|\boldsymbol{\mathrm{H}}-a\widetilde{\boldsymbol{\mathrm{H}}}\|\leq\epsilon_{\rm enc}$. We use Lemma A.1 to prepare an $(a+|h|,b+1,\epsilon_{\rm enc})$-QUBE $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}-h}$ for $\boldsymbol{\mathrm{H}}-h$ with the same cost (up to constants). Again, by definition, $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}-h}$ is an $(1,b+1,0)$-QUBE for a matrix $\widetilde{\boldsymbol{\mathrm{H}}}^{\prime}$ such that $\|(\boldsymbol{\mathrm{H}}-h)-(a+|h|)\widetilde{\boldsymbol{\mathrm{H}}}^{\prime}\|\leq\epsilon_{\rm enc}$.

Using standard techniques from holomorphic functional calculus (see e.g. [19, Lemma B.1]) we know that as long as $\epsilon_{\rm enc}\leq\frac{\epsilon_{\rm sgn}}{2}\cdot\frac{\sigma_{\min}^{2}(\boldsymbol{\mathrm{H}}-h)}{C_{\rm sgn}(a+|h|)}$, for some constant $C$, then it holds that $\left\|\operatorname{sgn}\left(\frac{\boldsymbol{\mathrm{H}}-h}{a+|h|}\right)-\operatorname{sgn}(\widetilde{\boldsymbol{\mathrm{H}}}^{\prime})\right\|\leq\epsilon_{\rm sgn}/2$. A sufficient upper bound for $\epsilon_{\rm enc}$ is to set $\epsilon_{\rm enc}\leq\epsilon_{\rm sgn}\frac{\Delta^{2}}{2C_{\rm sgn}(a+|h|)}$. It also holds that $\operatorname{sgn}\left(\frac{\boldsymbol{\mathrm{H}}-h}{a+|h|}\right)=\operatorname{sgn}\left(\boldsymbol{\mathrm{H}}-h\right)$ since $\operatorname{sgn}$ is scale invariant.

It remains to approximate $\operatorname{sgn}(\widetilde{\boldsymbol{\mathrm{H}}}^{\prime})$ with a polynomial (QET). From Weyl’s inequality, $|\lambda_{i}(\boldsymbol{\mathrm{H}}-h)-(a+|h|)\lambda_{i}(\widetilde{\boldsymbol{\mathrm{H}}}^{\prime})|\leq\epsilon_{\rm enc}$. Therefore, it follows that $\sigma_{\min}(\widetilde{\boldsymbol{\mathrm{H}}}^{\prime})\geq\frac{\sigma_{\min}(\boldsymbol{\mathrm{H}}-h)-\epsilon_{\rm enc}}{a+|h|}\geq\frac{\Delta-\epsilon_{\rm enc}}{a+|h|}\gg\frac{\Delta}{2(a+|h|)}$, where the last bound holds due to the aforementioned value of $\epsilon_{\rm enc}\leq\epsilon_{\rm sgn}\frac{\Delta^{2}}{2C(a+|h|)}\ll\Delta/2$. From [38, Lemma 3], we can find a polynomial $p(\cdot)$ with degree $l=O\left(\frac{2(a+|h|)}{\Delta}\log(\frac{2}{\epsilon_{\rm sgn}})\right)$, such that $|p(x)|\leq 1$ for all $|x|\leq 1$, and $|p(x)-\operatorname{sgn}(x)|\leq\epsilon_{\rm sgn}/2$ for all $x\in[-1,-\Delta/2(a+|h|)]\cup[\Delta/2(a+|h|),1]$. Consequently, $\left\|p(\widetilde{\boldsymbol{\mathrm{H}}}^{\prime})-\operatorname{sgn}(\widetilde{\boldsymbol{\mathrm{H}}}^{\prime})\right\|\leq\epsilon_{\rm sgn}/2$.

From [49, Thm. 2] (or [38, Thm. 1]), we can obtain an $(1,b+2,0)$-QUBE of $p(\widetilde{\boldsymbol{\mathrm{H}}}^{\prime})$ using $l$ queries of $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}-h},\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{H}}-h}^{\dagger}$, and $O((b+2)l)$ elementary gates. This is $\boldsymbol{\mathrm{U}}_{\rm sgn}$, and it is easy to see that it forms an $(1,b+2,\epsilon_{\rm sgn})$-QUBE for $\operatorname{sgn}(\boldsymbol{\mathrm{H}}-h)$, since:

[TABLE]

∎

A.2 Computing the trace of a block-encoded matrix

The main algorithm requires a subroutine to compute the trace of a block-encoded matrix.

Definition A.2** ([61]).**

Let $\rho$ be a density operator on $\mathcal{H}_{S}$ and let $\ket{\mathbf{0}}$ be some easy-to-prepare state in $\mathcal{H}_{S}$. A unitary $\boldsymbol{\mathrm{U}}_{\rho}$ on $\mathcal{H}_{S}\otimes\mathcal{H}_{A}$ is an $R$-preparation unitary of $\rho$ if we have

[TABLE]

where $\ket{\rho}=\boldsymbol{\mathrm{U}}_{\rho}\left(\ket{\mathbf{0}}\ket{\mathbf{0}}\right)$ and $\boldsymbol{\mathrm{U}}_{\rho}$ is implementable using $R$ elementary gates.

Proposition A.1**.**

Let $\boldsymbol{\mathrm{I}}_{n}$ denote the $2^{n}\times 2^{n}$ identity matrix. Then there exists a $2n$-preparation of $\frac{1}{2^{n}}I_{n}$ using $n$ CNOT gates and $n$ single-qubit gates on $\mathcal{H}\otimes\mathcal{H}$.

Proof.

Let $\ket{i}_{n}$ denote the $i^{\text{th}}$ basis state on the $2^{n}$-dimensional Hilbert space. Then we can prepare the state

[TABLE]

by first applying $\boldsymbol{\mathrm{H}}^{\otimes n}\otimes I_{n}$, where $\boldsymbol{\mathrm{H}}$ denotes the Hadamard gate, and then CNOT gates between the $j^{\text{th}}$ and the $(j+n)^{\text{th}}$ qubits for $j=0,\cdots,n-1$. Then, we have

[TABLE]

∎

Corollary A.2**.**

Let $\boldsymbol{\mathrm{A}}$ be an $N\times N$ matrix, and $\delta_{\rm tr},\epsilon_{\rm tr}\in(0,1)$. Let $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}}$ be an $(a,b,\tfrac{\epsilon_{\rm tr}}{2N})$-QUBE of $\boldsymbol{\mathrm{A}}$ with circuit complexity $\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}};\frac{\epsilon_{\rm tr}}{2N})$ such that

[TABLE]

where $\widetilde{\boldsymbol{\mathrm{A}}}:=(\bra{0}^{\otimes b}\otimes\boldsymbol{\mathrm{I}})\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}}(\ket{0}^{\otimes b}\otimes\boldsymbol{\mathrm{I}})$. There exists a quantum algorithm that produces an estimate $\xi$ of $\tr(\boldsymbol{\mathrm{A}})$ such that

[TABLE]

with probability at least $(1-\delta_{\rm tr})$. The algorithm has complexity:

Queries to $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}}$: $O\left(\tfrac{aN}{\epsilon_{\rm tr}}\log(\tfrac{1}{\delta_{\rm tr}})\right)$, 2. 2.

Number of qubits: $2\log(N)+b$, 3. 3.

Other elementary gates: $O\left(\log(N)\tfrac{aN}{\epsilon_{\rm tr}}\log(\tfrac{1}{\delta_{\rm tr}})\right)$ 4. 4.

Total complexity: $O\left((\log(N)+\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}};\tfrac{\epsilon_{\rm tr}}{2N}))\tfrac{aN}{\epsilon_{\rm tr}}\log(\tfrac{1}{\delta_{\rm tr}})\right)$.

Proof.

From Proposition A.1, there exists an $R$-preparation of $\rho=\frac{1}{N}I_{n}$ with cost $R=O(\log(N))$ elementary gates. From Lemma 5 of [61], we know that if $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}}$ is an $(a,b,0)$-QUBE of $\boldsymbol{\mathrm{A}}$ with circuit complexity $\mathcal{C}(\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}})$, and $\rho$ has an $R$-preparation unitary, then for every $\epsilon,\delta>0$ there exists an algorithm that produces an estimate $\xi$ of $\tr(\rho\boldsymbol{\mathrm{A}})$ such that

[TABLE]

with probability at least $(1-\delta)$. The algorithm uses $O\left(\frac{a}{\epsilon}\log(\frac{1}{\delta})\right)$ applications of a Grover operator that requires four uses of the R-preparation unitary and two queries to $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}}$. Note that $\boldsymbol{\mathrm{U}}_{\boldsymbol{\mathrm{A}}}$ is an $(1,b,0)$-QUBE of $\widetilde{\boldsymbol{\mathrm{A}}}$. Using [61, Lemma 5] on $\widetilde{\boldsymbol{\mathrm{A}}}$ with error $\epsilon:=\tfrac{\epsilon_{\rm tr}}{2aN}$, we can compute an estimate $\xi^{\prime}$ such that

[TABLE]

with probability at least $1-\delta_{\rm tr}$ and circuit complexity

[TABLE]

Then we can set $\xi=aN\xi^{\prime}$, which satisfies

[TABLE]

∎